close

Вход

Забыли?

вход по аккаунту

?

Расчет течений вязкой несжимаемой жидкости с помощью схем повышенной точности.

код для вставкиСкачать
НАУЧНЫЙ ВЕСТНИК МГТУ ГА
2012
№ 184
УДК 519.63
РАСЧЕТ ТЕЧЕНИЙ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ
С ПОМОЩЬЮ СХЕМ ПОВЫШЕННОЙ ТОЧНОСТИ
В.И. КОСТИН
Статья представлена доктором физико-математических наук, профессором Лобановым А.И.
В статье рассматриваются вопросы численного моделирования течений несжимаемой стратифицированной
жидкости. Применяется неявный алгоритм решения уравнений Навье-Стокса, записанных в приближении Буссинеска, основанный на применении компактных аппроксимаций третьего порядка точности для конвективных членов уравнений. Диффузионные члены уравнений аппроксимируются с четвёртым порядком точности. Приводятся
результаты расчётов.
Ключевые слова: компактные аппроксимации, стратифицированная среда, несжимаемая жидкость.
Введение
Одним из универсальных методов расчёта течений несжимаемой жидкости на протяжении
длительного времени является метод МАС [1] и целый ряд его модификаций [2-4].
В основе данной схемы расчётов для уравнений, записанных в переменных скорость - давление, лежит идея определения давления [2] (или поправки к давлению [5]) при условии, что на
каждом временном слое обращается в ноль величина дивергенции скорости.
Первоначальные попытки получения решения были связаны с применением явных алгоритмов вычисления независимых переменных на новых временных слоях [1; 2]. Однако при
проведении расчётов течений с малой вязкостью жесткие ограничения на критерий устойчивости разностных схем приводили к дополнительным сложностям при реализации вычислительных алгоритмов. Обладающие большим запасом устойчивости неявные разностные схемы [3-7]
в значительной мере снимают эти проблемы. Другой возможностью получить решение удовлетворительного качества (с хорошими свойствами аппроксимации и монотонности) является выбор надлежащей аппроксимации исходных уравнений.
Для этого часто используют ориентированные разности на многоточечном шаблоне [9], что
при использовании неявных схем приводит к существенному усложнению алгоритма расчёта.
Альтернативой этому выступает схема расчётов, основанная на использовании трехточечного шаблона для аппроксимации конвективных и диффузионных членов, входящих в исходные уравнения [4-8; 10]. Её преимущество заключается в том, что она обладает хорошими дисперсионными и диссипативными свойствами и имеет высокий порядок аппроксимации: от
третьего и выше [10-12], а получающаяся система линейных уравнений имеет трехдиагональную форму, которая легко решается методом прогонки.
В работе [8] реализован алгоритм расчёта, основанный на аппроксимации конвективных и
диффузионных членов уравнений с третьим порядком точности, позволяющий проводить исследования в широком диапазоне определяющих параметров течений неоднородной несжимаемой жидкости.
В настоящей работе рассматривается алгоритм, который содержит аппроксимацию диффузионных членов уравнений четвёртого порядка точности.
Расчет течений вязкой несжимаемой жидкости …
63
1. Постановка задачи и основные уравнения
Рассматривается течение полностью перемешанной жидкости в стратифицированной среде
в рамках вязкой несжимаемой жидкости с линейной стратификацией плотности (или температуры) по глубине.
Представленная задача описывается системой уравнений Навье-Стокса, записанных для переменных скорость - давление в приближении Буссинеска, когда отличие плотности среды от ее
невозмущенного значения учитывается только в силах плавучести.
Таким образом, изучаемое явление рассматривается в декартовой системе координат
( 0, x, y, z ) или ( 0, x1 , x2 , x3 ) , причем начало координат совпадает с центром возмущённой области. Ось Ox направлена вдоль потока, Oz - вертикально вверх, а Oy дополняет систему координат до правой тройки.
Стратификация среды в невозмущенном состоянии определяется распределением плотности ρ = ρ 0 ( z ) (или температуры T = T0 ( z ) ). Как известно, степень стратификации характеризуется частотой Брента-Вяйсяйля, определяемой следующим соотношением
g ∂ρ
N 2 = − ⋅ 0 = α g,
(1)
ρ* ∂z
где ρ* = ρ 0 ( 0 ) , a g ускорение свободного падения.
r
Если через u1 , u2 , u3 обозначить проекции вектора скорости V на соответствующие оси координат, то согласно [2-5, 13, 14] течение вязкой несжимаемой жидкости описывается следующей системой уравнений
D ρ ui
∂p
=−
+ ρ v∆ui − ρ g i ,
(2)
Dt
∂xi
∂ui
(3)
= 0.
∂xi
Здесь g1 = g 2 = 0, g 3 = g - ускорение свободного падения; ∆ - оператор Лапласа, а в формуле
Df
(3) - суммирование по повторяющемуся индексу.
- полная производная по времени.
Dt
В приближении Буссинеска система (2) примет вид (4)
Dui
1 ∂ρ
ρ
=−
+ v∆ui − g i .
(4)
Dt
ρ* ∂xi
ρ*
Для определения величины плотности ρ к уравнениям (3-4) необходимо добавить следующее соотношение [3-5, 15]
D( ρ − ρ0 )
∂ρ
= −ui 0 .
(5)
Dt
∂xi
Исследуемая перемешанная жидкость представляет собой бесконечный цилиндр радиусом L
с образующей – осью Ox, таким образом рассматривается нестационарное течение, генерируемое участком перемешанной жидкости ( y 2 + z 2 ≤ L2 ) при наличии устойчивой стратификации
ρ = ρ0 ( z ) невозмущенного потока с частотой Брента-Вяйсяйля (1).
Считается, что в начальный момент времени жидкость покоится, а плотность внутри перемешанной области ("пятна") постоянна
v(0, y, z ) = w(0, y, z ) = 0 ρ (0, y, z ) = ρ* = const.
(6)
Для линейной стратификации
ρ0 ( z ) = ρ* (1 − α z ).
(7)
64
В.И. Костин
1
выбрать в качестве характерных значений масштабов длины, скорости и
N
времени рассматриваемого явления, то в уравнениях (3)-(5) можно перейти к безразмерным величинам
u
x
vi = i , yi = i , τ = Nt.
(8)
U
L
Безразмерное возмущение плотности при этом вводится по формуле
ρ − ρ0
ρ=
.
(9)
∂ρ0
L
∂z
Величину градиента давления обычно рассматривают в виде

∂π
1  ∂p
= 
+ ρ0 gi  .
(10)
∂yi ρ*  ∂yi

Следовательно, проведя процедуру обезразмеривания исходных уравнений, легко получить
следующие соотношения
∂vi
= 0,
(11)
∂yi
Dvi
∂π
1
ρ
=−
+
∆vi + 2 δ i 3 ,
(12)
∂yi Re
Dτ
Fr
Dρ
= −viδ i 3 .
(13)
Dτ
Здесь безразмерные числа Рейнольдса и Фруда определены следующими соотношениями (N период Брента-Вяйсяйля (1))
UL
U
Re =
, Fr =
.
(14)
v
NL
Поскольку, в соответствии с (6), в начальный момент времени (τ = 0) вне "пятна" жидкость
покоится, то в этой области p = 0 , а внутри "пятна" (в соответствии с формулами (1), (6)-(7), (9)
p = − z . Так что справедливо
Если L,U и
 0 при y 2 + z 2 ≥ L2
(15)
2
2
2
− z при y + z ≤ L .
Давление, определяемое (10), при этом задаётся интегрированием величины ρ ( y, z ) таким
образом, чтобы оно соответствовало гидростатическому соотношению во всей расчётной области [3-5]
1
π ( 0, y, z ) = 2 ∫ ρ ( y, z )dz.
(16)
Fr
Таким образом, сформулирована задача об эволюции "пятна" в устойчиво стратифицированной среде (7) с начальными условиями (6), (15), (16), которая описывается системой уравнений (11)-(13).
ρ =
2. Алгоритм решения
Для численного решения полученной системы (11)-(13) в расчётной области вводится разностная сетка ωh с узлами, имеющими координаты ( tm , yi , z j )
ωh : tm = (m − 1)τ ; yi = (i − 1)hy ; z j = ( j − 1)hz .
(17)
65
Расчет течений вязкой несжимаемой жидкости …
Значения всех функций, входящих в рассматриваемые уравнения, задаются в узлах сетки
ωh (17)
f ijm = f ( tm , yi , z j ) .
(18)
Решение системы (11)-(13) основано на определении поправки к давлению из уравнения
Пуассона с некоторой правой частью, которая получается из условия несжимаемости divV = 0
на каждом временном слое t = tm [3-5].
Согласно [3-5, 13-15] уравнения (10) можно записать в виде
∂vi ∂vi vk
∂π
1
ρ
+
=−
+
∆vi − 2 δ i 3 .
(19)
∂τ
∂yk
∂yi Re
Fr
Если ввести оператор
∂f v
1
ρ
L ( fm ) = m k −
∆f m − 2 δ m 3 ,
(20)
∂yk
Re
Fr
то соотношения (10) примут вид
∂vi
∂π
+ L ( vi ) = −
,
(21)
∂τ
∂yi
или для каждой из компонент вектора скорости в поперечном сечении
∂v
∂π
+ L (v) = −
,
(22)
∂τ
∂y
∂w
∂π
+ L ( w) = −
.
(23)
∂τ
∂z
Каждое из уравнений (22) и (23) решается в два шага, например, для (22) имеем
v% − v m
∂π m
I.
+ L ( v% ) = −
;
(24)
∂y
τ
II .
v m +1 − v m
τ
+ L ( v% ) = −
∂π m +1
.
∂y
(25)
∂π m
— некоторый разностный аналог производной от давления, а v% — промежуточное
∂y
значение скорости на новом временном слое.
Из (24) и (25) следует
∂ (π m +1 − π m )
v m +1 = v% − τ
.
(26)
∂y
Здесь
Аналогично, применяя процедуру (24)–(25) к уравнению (23), легко получить
∂ (π m +1 − π m )
m +1
w = w% − τ
.
∂z
(27)
Теперь, дифференцируя уравнения (26) и (27) по y и z соответственно и складывая полученный результат, можно получить уравнение для искомой поправки к давлению δπ = π m +1 − π m
∂ 2 (δπ )
∂ 2 ( δπ )
∂ν m +1
∂wm +1
∂ν%
∂w%
.
(28)
+
=
+
−τ
−τ
∂y
∂z
∂y
∂z
∂y 2
∂z 2
Учитывая то, что жидкость несжимаема, на каждом временном шаге должно выполняться уравнение (11), то есть справедливо
66
В.И. Костин
∂ν m +1
∂wm +1
+
= 0,
∂y
∂z
(29)
следовательно, из (28) следует
∂ 2 (δπ )
∂y
2
+
∂ 2 (δπ )
∂z
2
=
1  ∂ν%
∂w% 
+

.
τ  ∂y
∂z 
(30)
Это и есть искомое уравнение Пуассона для поправки к давлению.
Таким образом, алгоритм решения исходной системы уравнений следующий (вычисление
неизвестных величин на временном слое t = tm +1 по известным значениям функций на m-м слое
t = tm ).
1. Определяются из соотношения типа (24) компоненты скорости ν% и w% .
2. Из уравнения Пуассона (30) находится поправка к давлению δπ .
3. Определяются точные значения скорости ν m +1 и wm +1 из формул (26) и (27).
4. Из уравнения (13) находится значение плотности p m +1.
3. Аппроксимация исходных уравнений
Достаточно подробное обсуждение вопроса аппроксимации разностного оператора L, введенного формулой (20), было проведено в работе [5]. Там было отмечено, что для задач рассматриваемого типа, когда требуются описания течений, которые содержат сильные неоднородности, используются разностные схемы, обладающие большим запасом устойчивости и
имеющие хорошие свойства монотонности и малую схемную диссипацию. Как показано в целом ряде работ [4-8, 10-12], такими свойствами обладают компактные аппроксимации повышенной точности (выше второго порядка). В [8] реализован алгоритм, основанный на аппроксимации как конвективных, так и диссипативных членов уравнений Навье-Стокса с третьим
порядком точности. Вместе с тем, достаточно несложно обобщить этот алгоритм на случай, когда диссипативные члены уравнений Навье-Стокса учитываются ещё более точно - с четвёртым
порядком аппроксимации.
Для этого оператор L, заданный согласно (20)
df ν dfw 1 d 2 f
p g
L( f ) =
+
−
+ 2 i,
(31)
2
dy
dz Re dy
Fr g
необходимо представить в виде
p g
Lf = Ly f + Lz f + 2 i .
(32)
Fr g
Разностная аппроксимация данного оператора проводится следующим образом. Например,
для оператора Ly f применяется следующий подход
 3 h4 
1 −1
B
f
o
∆
+
h +
,
y
2y
h
Re
Re 

где аппроксимация второй производной выполняется с четвёртым порядком точности
d2 f
d2 f
−1
r=
=
B
∆
=
+ O ( h4 ) .
y
2 yf
2
2
dy h
dy
Ly f
= Ay−1∆ y uf −
Операторы By и ∆ 2 y имеют следующий вид
By f = ( E +
1
1
1
10
1
∆ 2 ) f = f i + ( f i +1 − 2 f i + fi −1 ) =
fi +1 −
fi +
fi −1 ,
12
12
12
12
12
(33)
(34)
67
Расчет течений вязкой несжимаемой жидкости …
∆2 y f =
1
( f i +1 − 2 fi + f i −1 ) .
h2
Причем, в соответствии с (34)
By r = ∆ 2 y f .
(35)
Аппроксимация конвективных членов осуществляется парой следующих трехточечных
операторов
df
df
= Ay−1 ∆ y f =
+ O ( h3 ) .
(36)
dy h
dy
Введённые разностные операторы вычисляются по следующим формулам
s
1
Ay f = ( A0 − ∆ 0 ) f ; ∆ y f =
(∆ 0 − s∆ 2 ) f .
4
2h
Разностные операторы A0 , ∆ 0 и ∆ 2 имеют следующий вид
1
4
1
A0 f i = f i −1 + f i + f i +1 ;
6
6
6
∆ 0 f i = f i +1 − fi −1 ;
∆ 2 f i = fi +1 − 2 f i + f i −1 .
В качестве s берётся s = signV .
Следовательно, применяя компактные разности (33), (34), (36) к уравнениям (11)-(13) в
форме (21), можно получить разностную схему их решения. Так, для скорости v из (22) получается следующее соотношение
ν −ν m
+ Ly (ν ) + Lz (ν ) = −∇ y p m .
(37)
τ
Если ввести обозначение (38), то (37) примет вид (39)-(40), поскольку будет справедливо
ν =ν m +τ ∈
ν −ν m
∈=
,
(38)
τ
E ∈ + Ly (υ + τ ∈) + Lz (υ m + τ ∈) = −∇ y p m ,
(39)
E ∈ +τ Ly (∈) + τ Lz (∈) = ϕ ,
(40)
m
m
ϕ = − Ly (ν
m
m
) − L (ν ) − ∇
m
z
m
y
p .
(41)
Соотношения (40)-(41) удобно представить в виде
[ E + τ Ly + τ Lz ] ∈= ϕ m .
(42)
Последнее соотношение допускает факторизацию (с погрешностью О(τ2))
[ E + τ Ly ][ E + τ Lz ] ∈= ϕ m .
(43)
Полученное уравнение решается в два приема: вводится новая переменная
ω = [ E + τ Ly ] ∈
(44)
и относительно неё решается уравнение (43), принявшее вид (45)
[ E + τ Lz ]ω = ϕ m .
(45)
После определения из (45) величины ω находится искомое значение ∈ (38) из соотношения
(44) и далее - переменная v.
Более подробно процедура вычисления неизвестной величины ∈ выглядит следующим
образом.
68
В.И. Костин
С учётом того, что
1 −1
Bz ∆ 2 z f ,
Re
1 −1
Ly f = Ay−1∆ yν f −
By ∆ 2 y f ,
Re
уравнения (44) и (45) будут иметь вид
Lz f = Az−1∆ z wf −
[ E + τ Az−1∆ z w −
τ
Re
[ E + τ Ay−1∆ yν −
τ
Bz−1∆ 2 z ]ω = ϕ m ,
(46)
By−1∆ 2 y ] ∈= ω.
(47)
Re
Теперь, c учётом соотношения (35), из (46) следует (48)-(49)
[ A z +τ∆ z w]ω −
τ
Azϑ = Az [ϕ m ],
Re
∆ 2 z ω − Bzϑ = 0
Полученные соотношения можно записать в матричном виде
τ


m
 [ A z +τ∆ z w] − Re Az   ω  =  Az [ϕ ]  ,




  ϑ   0 
∆
−
B
2z
z 

(48)
(49)
(50)
и, если ввести обозначения
τ


[ A +τ∆ z w] −
Az  r  ω 
G = z
, x =  ,
Re


ϑ 
∆2z
− Bz 

r  A [ϕ m ] 
f = z
,
 0 
(51)
то получится следующая система линейных уравнений
Gx=f
(52)
с известными матрицей G и вектором правой части f. Поскольку все операторы A и ∆ , как показано выше, вычисляются на трехточечном шаблоне, то матрица G будет иметь трехдиагональный вид и, следовательно, соотношение (52) можно записать в каждой точке разностной
сетки ωh в следующем виде
r
r
r
r
Aj x j −1 + B j x j +1 − C j x j = f j , j = 2,..., N − 1 .
(53)
Очевидно, к полученным уравнениям необходимо добавить условия на границах рассматриваемой области решения
r
r r
x1 = Ψ1 x2 + µ1 ;
r
r
r
xN = Ψ N xN −1 + µ N .
Полученная система уравнений решается методом векторной прогонки [16], который подробно рассмотрен в [8].
После определения величины ω из системы уравнений (53) описанным выше методом определяется неизвестное значение ∈ из уравнения (47).
4. Результаты расчётов
С помощью предложенного алгоритма было проведено численное исследование рассмотренной выше задачи об эволюции перемешанной жидкости ("пятна"). Приведённые результаты
расчётов относятся к случаю, когда число Фруда Fr = 1, а величина числа Рейнольдса варьировалась от 402 до 405.
69
Расчет течений вязкой несжимаемой жидкости …
Рис. 1
Рис. 3
Рис. 2
Рис. 4
На рис. 1 - 4 представлено распределение скорости ν ( y ) на оси "пятна" ( z = 0) в моменты
времени t = 3 и t = 9 для числа Re=402. На следующих двух рисунках изображены распределения скорости ν ( y ) на оси "пятна" ( z = 0) в моменты времени t = 6, t = 12 для числа Re=405.
Приведённые данные показывают, что конфигурация поперечной скорости течения в рассматриваемом диапазоне числа Рейнольдса на протяжении всего времени исследования носит
монотонный характер. Наблюдается хорошее совпадение результатов расчетов по схемам
третьего и четвёртого порядка точности.
На следующих двух рисунках представлено распределение вертикальной скорости по глубине w( z ) в момент времени t = 12 для числа Рейнольдса Re=402. Здесь рассмотренная характеристика течения даётся как в центре возмущенной области (рис. 5), так и в области максимальных амплитуд (рис. 6). Приведённые данные свидетельствуют о том, что схемы третьего и
четвёртого порядка точности дают очень близкие результаты.
Рис. 5
Рис. 6
70
В.И. Костин
На рис. 7, 8 представлены конфигурации поперечной скорости w( z ) в момент времени t = 12
в области максимальных амплитуд и в центре возмущения среды для числа Re=405. Здесь также
наблюдается очень хорошее совпадение результатов расчётов для обеих рассматриваемых схем.
Рис. 7
Рис. 8
Распределение функции плотности в области максимальных амплитуд для чисел Рейнольдса
Re=402 и Re=405 показано на рис. 9, 10.
Рис. 9
Рис. 10
Как и для уже рассмотренных характеристик исследуемого течения - величин скорости,
приведённые данные говорят о монотонном поведении функции плотности на протяжении всего представленного временного интервала T ≤ 12 во всём диапазоне числа Рейнольдса
3.102 ≤ Re ≤ 3.105 ,
Приведённые результаты расчётов показывают, что применение компактных разностных
схем повышенного порядка точности для аппроксимации уравнений Навье-Стокса даёт вполне
удовлетворительное описание процессов, содержащих сильные неоднородности. Использование схемы четвёртого порядка точности для аппроксимации диффузионных членов уравнений
приводит к практически тем же результатам, что были получены для алгоритма с третьим порядком точности в работе [8].
ЛИТЕРАТУРА
1. Harlow F., Welch J. Numerical calculation of time-dpendent viscous incompressible flow of fluid with free surface. - Physics of Fluids, 1965, v.8, №12, pp.2182-2189.
2. Белоцерковский О.М., Гущин В.А., Щенников В.В. Метод расщепления в применении к решению задач
динамики вязкой несжимаемой жидкости // Журнал вычислительной математики и математической физики. - 1975.
- №1. - Т.15. - С. 197-207.
71
Расчет течений вязкой несжимаемой жидкости …
3. Даниленко А.Ю., Костин В.И., Толстых А.И. О неявном алгоритме расчёта течений однородной и неоднородной жидкости. - М.: ВЦ АН СССР, 1985.
4. Костин В.И. Численное моделирование некоторых задач океанологии на основе компактных схем третьего
порядка. - М.: ВЦ АН СССР, 1991.
5. Костин В.И. Применение компактных аппроксимаций третьего порядка к решению некоторых задач неоднородной несжимаемой жидкости: Деп. в ВИНИТИ № 630-В90 // Труды XV конф. молодых учёных МФТИ. - М., 1990.
- Ч. 1. - С. 135 - 148.
6. Копысов А.Н., Костин В.И. Излучение звука турбулентной струей в спутном потоке при малых числах
Маха: Деп. в ВИНИТИ № 630-В90 // Труды XV конф. молодых учёных МФТИ. - М., 1990. - Ч. 1. - С.120 - 134.
7. Копысов А.Н., Костин В.И. Численное моделирование излучения звука турбулентным течением в стратифицированной среде. - М.: ВЦ АН СССР, 1991.
8. Костин В.И. Об одном методе расчета течений вязкой несжимаемой жидкости // Речной транспорт (XXI век).
- 2011. - № 4 (52). - С. 64-70.
9. Рыков В.В. Численное моделирование пространственных течений несжимаемой вязкой жидкости. - М.: ВЦ
АН СССР, 1983.
10. Толстых А.И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. - М.: Наука, 1990.
11. A.I. Tolstykh. Development of arbitrary-order multioperators-based schemes for parallel calculations. 1: Higherthan-fifth-order approximations to convection terms // J. Comput. Phys. 2007. V.225. P.2333-2353.
12. Липавский М.В., Толстых А.И., Чигерёв Е.Н. О схеме с мультипликаторными аппроксимациями 9-го
порядка для параллельных вычислений и её применении в задачах прямого численного моделирования // Журнал
вычислительной математики и математической физики. - 2006. - Т. 46. - № 8. - С.1433-1452.
13. Монин А.С., Яглом А.М. Статистическая гидромеханика. - М.: Наука, 1964. - Ч.1.
14. Кочин Н.Е., Кобель И.А, Розе Н.В. Теоретическая гидромеханика. - М.: Физматгиз, 1963. - Ч. 1, 2.
15. Океанология. Физика океана. - Т. 1. Гидрофизика океана; под ред. А.С. Монина. - М.: Наука, 1978.
16. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978.
ON THE NUMERICAL METHOD FOR CALCULATING VISCOUS INCOMRESSIBLE FLUIDS
Kostin V.I.
This paper is devoted on the numerical simulations of incompressible stratified fluid. The implicit algorithm for solving the Navier-Stokes in Boussinesq approach equations based on the use of compact approximations of the third order for
convective terms for equations is used. The diffusion terms of equations are approximated by the fourth order of accuracy.
The results of the calculations.
Key words: compact approximations, stratified medium, incompressible fluid.
Сведения об авторе
Костин Владимир Иванович, 1959 г.р., окончил МФТИ (1982), Финансовую академию при Правительстве РФ (2000), кандидат физико-математических наук, доцент, ректор МГАВТ, автор более 30 научных работ, область научных интересов – численное моделирование задач аэрогидродинамики, вопросы
анализа рисков и экономической безопасности, организация и управление образовательными процессами.
Документ
Категория
Без категории
Просмотров
8
Размер файла
2 427 Кб
Теги
точности, помощь, расчет, вязкой, повышенных, схема, жидкости, течение, несжимаемой
1/--страниц
Пожаловаться на содержимое документа