CC BY
32
УДК 519.62
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ УРАВНЕНИИ НАВЬЕ-СТОКСА
В СЛОЖНОЙ ОБЛАСТИ
© 2007 г. А.Б. Усов
It is offered the semi-implicit scheme of the method of the final differences for the numerical research of the nonlinear Navie-Stoks equations in curvilinear physical field, depending from time. The case of the compressed liquid motion is explored. The algorithm of the count on this final differences scheme is specified. Its conditions of stability is deduced. The comparative analysis between this scheme and others used scheme is organized.
Вопросам численного решения нелинейных уравнений движения вязкой жидкости - уравнений На-вье-Стокса - посвящено большое количество работ [1-5]. С развитием вычислительной техники все чаще численно решаются полные уравнений Навье-Стокса. Наиболее распространенными разностными схемами решения являются схемы переменных направлений [3], дробных шагов [5], вариант схемы Лакса-Вендрофа, предложенный Томменом [1], схемы Мак-Комака [4], Бима-Уорминга [2] и др. Однако остается множество проблем при численной реализации алгоритмов решения конкретных задач, в которых выбор расчетной схемы проводится с учетом целей исследования, входных данных, рассматриваемого промежутка времени, математической модели, в состав которой входят решаемые уравнения, а также рельефа дна водной системы, разбиения на камеры, длины каждой камеры и многими другими факторами.
Постановка задачи
Рассмотрим задачу о плоском движении вязкой сжимаемой жидкости в области криволинейной формы с твердыми стенками (рисунок). Введем декартову прямоугольную систему координат 0x2; О - точка пересечения свободной поверхности жидкости и левой границы исследуемого участка реки; ось Ох направлена вдоль реки, ось 02 - вертикально вниз. Изотермическое движение жидкости в 0x2 описывается системой уравнений Навье-Стокса [6]:
8V
Г Г г г
-+vx — + v, — = —Vp + /Ar + -/Wdv V)+F, (1)
dt дх dz р 3
dt dx ' dz р dp d(vxp) d(vzp)
dt
dx
dz
= 0 , p = D p.
Здесь р - плотность жидкости; V=(vx, vj - вектор скорости частиц жидкости; р - гидродинамическое давление; и - коэффициент динамической вязкости жидкости; F=(Fx, Fz) - вектор динамической составляющей внешних массовых сил, действующих на жидкость и отнесенных к единице массы; D = const -квадрат скорости звука в жидкости; А - оператор Лапласа на плоскости; V - оператор Набла; div V - дивергенция вектора скорости жидкости.
Физическая область
Граничные и начальные условия: на свободной поверхности (при 2 = £ (х, 0)
1 dv„
1 8v„
3 dx
3 dz
dv dv
г vx
(nxrz + nzrx)
8vz ^ dvz V 8z dx
. г c'v. dv
dz
dx
= 8 2-
d<Z dt
А =.
1 +
dx д<; дх 1
nz =--; пх=- ,
А дх дх А
тг = -п,
условия «прилипания» на дне (при г = Н(х, I)) на левой и правой границах (при х=0 и х=Ь)
V = V = 0;
(3)
(4)
(5)
zu = In
x0 =ln ■
Обратно
zj = 0,5
Д + 2zt -1 ß\ ~2zj -1
ß2 + 2xx -1
/In
fl+1 А
(8)
/ln
vth-Л
-Ä
в начальный момент времени
ух(х,г,0) = /0(х,г);
vг(x,z,0) = р{х,г,0) = р0(х,г) . (6)
Здесь L - длина рассматриваемой области; функции f0, /¡, р0 определяют поле скоростей и плотность жидкости в начальный момент времени; р* - плотность жидкости в невозмущенном состоянии; п = (пх ,п2); т= (тх ,т2) -вектора внешней к свободной поверхности жидкости нормали и касательной в момент времени t.
Задача (1) - (6) решается методом конечных разностей.
Численная схема метода конечных разностей
Процесс численного решения уравнений Навье-Стокса значительно упрощается с помощью специально построенной расчетной сетки. Используя приведенные в [7] преобразования координат, получим равномерную, не зависящую от времени сетку в вычислительной плоскости, хотя узлы сетки в физической плоскости подвижны и расположены неравномерно.
Изучаемая область в физической плоскости имеет вид, изображенный на рисунке. При расчетах учитываются возвышение свободной поверхности жидкости и рельеф дна.
Расчетная область является криволинейной и зависит от времени, поэтому при переходе с одного временного слоя на другой необходимо пересчитывать значения искомых функций на пространственных слоях. Избежать такого пересчета позволяет переход к переменным (хь 2\), переводящим переменную (по времени) криволинейную физическую область В в постоянную вычислительную область <2 .
*!=(*-«-)/(#-?); х1 =х/ь. (7)
Таким образом, физическая область С, < г < Н; 0 < х < Ь переходит в вычислительную 0 < х, < 1; 0 <гг <1.
Вблизи свободной поверхности жидкости, дна, боковых границ происходит резкое изменение скорости жидкости и ее плотности. Для того чтобы уловить эти изменения, необходима мелкая пространственная сетка. Считать во всей вычислительной области на такой сетке нельзя из-за ограниченности памяти ЭВМ, поэтому счет проводится на неравномерной разностной сетке.
Введем в вычислительной области 2 неравномерную сетку (хь 2Х), которая сгущается вблизи границ. Для измельчения разностной сетки вблизи границ сделаем замену переменных, переводящую 2 в область Я с равномерной сеткой
1 + <д+1)/(А-1)~
Xj = 0,5
(ß2+i)<ß2+i)/(ß2-i)y -l-ß:
/
\+iß2+i)/(ß2-iy
Параметры Д, /?2 задаются вычислителем и влияют на густоту разностной сетки вблизи боковых границ, дна и свободной поверхности. Чем они ближе к единице, тем сильнее сгущается сетка вблизи границ.
Переходя в уравнениях Навье-Стокса к координатам О х02°, получим: уравнение неразрывности
.0 Я „ я „О я/
dp | dz" dp | gx» 8(/»x) |
dt dt dzö dx dx0 , dz0 8(pvx) , dz0 d(pvz)
dx dz0
dz dz
= 0;
(9)
уравнение движения в проекции на ось Ох
+ v„
dz0 8vx
о
dt dt dz dz0 dvv
- + v„
dx d vv
dx dx0
dx dz
dz0 dv
~дГа70
= -D
dx^_dp_ dz0 dp dx дх° +
dz 0
dz
4 + —
3
(10)
4 + —
3
dx0
dx
d2vr „ dz0 dx0 d2vv dLx0 8vr
J„0 я,, ^
dv
0
n 0 ^ 0 ^2
dz ex с v7
d( x0)2 fd2 z0
5z 2
+ 2-
z x z 0 x 0 x 2 x
4S2 z
3 8 x2
z x z 0 x0
0 0 2 dz dz ovz
dz dx d(z0)2
20 д z dvz
z x z _0
Уравнение движения жидкости по оси 02 выписывается аналогично, уравнение состояния сохраняет свой вид. Граничные и начальные условия запишутся в виде: на неподвижных границах (при х° = 1; х° = 0 и при 2°=1)
vx =0;v,=0 ;
(11)
на свободной поверхности жидкости (при г =0)
2
= v
z
1
Т. = пх
/
0
/
/
0
0
x
v
x
2
2
0
+
(
D{p-p*) = 2/j.
(nl-1)
(dx° дVv 5z°dV„1
dx dx0 dx dz0
/2 OZ
L dz0 dv,
(12)
+ nxnz
dz dv, dx dv, dz 8vr
■ +
dx dz0 dx dx° dz dz0
(«xTz + nzTx)
dxö 8V7 dz0 5V
■> J
0
fdz0 дv,
x z 0 x x0 z z 0
+ 2
dz dv. dz dz
(dz0 dvr dx° dvr )
dx dz0
dx dx°
d£ dt
= v- T = .
T
1 + d£ dx°
8x° dx
d£ dx 1
dx° dx T
(13)
tk=k At; x,°= (/ - 0,5)/M; Z}°= (y—0,5)/ N2; i,j=QX---N!(N2); ^=0,1,2,...,
(14)
где А/ - шаг сетки по времени; Л'/, Ы2 - число узлов разностной сетки по координатам х0, z0 соответственно.
Используется неразнесенная сетка - все неизвестные функции ищутся в узлах (х0, z0, ), определяемых (14). Для удовлетворения граничных условий узлы разностной сетки вблизи границ располагаются на расстояниях 1/(2 N1) и 1/(2 N2) от них, вводятся фиктивные узлы: (хд,г°); (х°,г^2+1);
(х°,гд); г,у = 0Д,...,ЛГ1+1(ЛГ2+1) и используются
соотношения /,
'j
-f (x0, z 0, tk ) = 05(ftrn, ,-ftn2,j ) =
-0,5 (fij+1 / 2 ftj-1 / 2 ) :
аппроксимации по пространственным переменным. Такой подход позволяет удовлетворить условия на границах, не вводя дополнительных граничных условий для метода конечных разностей. Наличие адвективных слагаемых в уравнениях (9), (10)
дvr 8у7 8 р др ч (V. —-: уг —- V.-; -), заставляет использо-
dz dx dz dx
вать разностные схемы, аппроксимирующие эти слагаемые «против потока». Используется неразнесенная сетка, производные аппроксимируются разностными соотношениями с первым порядком аппроксимации по пространственным и временной координатам.
Заменяя производные в (9)-(13) разностными соотношениями, получим конечно-разностный аналог задачи. В разностных уравнениях движения в проекции на ось Ох значения ух берутся на (к+1)-м временном слое, все остальные функции - на к-м; в проекции на ось Oz значения функции р - на к-м слое по време-
ни, все остальные функции - на (к+1)-м; уравнение неразрывности аппроксимируется по неявной схеме метода конечных разностей и рассматривается не только во внутренних узлах разностной сетки, но и на границах расчетной области. Разностные уравнения движения в проекции на ось Ох примут вид (Ах, Аг -шаги сетки по пространственным координатам):
Производные от функций z0 (х, х, 1), х0 (х) по х ,х, / определяются из (7), (8). Методом конечных разностей на равномерной сетке решается задача (9) - (13). Введем разностную сетку
/ , 1к ) _ 0,5(Л+1/2,/ Ji-1/2,j ) имеющие второй порядок
Л Рг,3
f f xk+l / \k f xk+l / xk+l
Оx)i,j "OxXj , .к (vx),,j ~(vx)i-ij
At
" + (Vx tj
Ax
dx
Я 0
dz
dt
dz0
dz
z0
dx \
(Vx )u +
(15)
(Vz )k
k
z->lJ
(v )k+1-(V )k+\ v x'l, j V x'l, j-1
z
= -D
f J- J-Рг,} -Рг-l.j
Ax
i,j
J- J-Pu~Pu-\
Az
i,j
+Fk +
i,j
+ iu
8
+3
4 (v x )l+1, j
-2(V
x 'l, j + (Vx )г- 1j
(Ax)2
dz
dz
Л+1
к 1 0 dx
8x
/ \k-i-i / \k+1 / \k+1 . / \k+1 (v*)i+l,/+l ~(Vx)z+l,/-l ~(Vx)z-l,/+l + (Vx)
x Л- 1,j-1
Га 0\2
dz
v /
5x v У
4AxAz
k
/ \k+1 о/ \k+1 , / \k+1 (Vx)i,j+1~ 2(VxXj +(Vx j
',j J
82 x0
5x2
k (V )k+}-(V )k*.
x
( л -,2„
4 dz dz
^^k (V )^^(V )k+
3 5 x2
dz
2
M
z
1
Ö2 z 0
^ dxdz
i,j
к (Vz tj- (Vz )k,j-1
z
Ji,j
(Vz )k+1,;+^(Vz )k+ 1,^^(Vz )k-1,^^(Vz )k
Я 0 " dz k я 0 8x
_dz dx
j
z )l- 1,j-1
OZ
az
az
dx
^A^Az
Разностные уравнения движения в проекции на ось Ох и уравнения неразрывности выписываются аналогично. При построении разностной схемы использовались идеи метода расщепления: система разностных уравнений расщепилась на три подсистемы (для определения функций V,.. . р). которые решаются независимо друг от друга. Уравнения движения (10) были записаны в неконсервативной форме, чтобы обеспе-
0
k
x
k
k
nzTz
k
2
k
k
0
0
x
z
1
x
x
2
k
0
= "x
= ~nz ■
x
3
x
к
+
2
3
4
3
k
k
k
чить возможность расщепления исходной системы разностных уравнений на независимые подсистемы.
Условия на свободной поверхности жидкости в разностной форме примут вид
\п2Л-ИЪ) *
f Г Ч&+1 Г >.к+\ , S xk+1 , Ч&+1
(vx )г,1 "(vx )г-1,0 + (vx )г,0 -(vx );-1,1
2Ax
(16)
ox
dz0
8z
z \к+1 / \к+1
(vx )k 1 "(vx )i, 0
¿,1/2 к
tSz
0 dz k
x i,1 / 2
¿,1/2
(v^^Jö-i)+
ff,. \клЛ.
+(nxnz)f
k+1
ВДГ "(Vj-O
Az
я и öz
~s7
¿,1/2
/ \k+1 / \k+1 , / 4k+ 1 / \k+1 (vz\i -(vJi-io+CvJyo -(vJi-u
8xv dx
i, 1/2
+ 2 ( т *
(vx )k1 "(vx )k-l 10"(vx )k-1 \+(vx )
2Ax
8x 8x
¿,1/2
(vx ^-(vx )k0 0 z0 k
Az x i,1/2
(17)
+2(«zrz)f
5z 0 5z
(vz ^"^z )k,0
,1/2
z
(t,. \k+1
k+1
(vx )y
Az
я 0 öz
az
,1/2
8x
dx
, (vz )g-(vz )k-1^-(vz )k-14+(vz )k0 +
¿,1/2
(vJu "(V2)*o
Az
я о dz
dx
i, 1/2
Tk 2 ¿1/2'
*-k+l _ rk
2j_~ ¿
&t
dx
(vx )k,,1+(vx )k0
1/2 к . r \k
2
Ax
(18)
rpk _
' 1
i+
3x°
Sx
,1/2
rk rk bi ~hi-1
x
(«x)f =
gx» 3x
¿,1/2
1 .
2Ax Tk '
(гЛ=(пЛ- (т,)к = (п,)к .
Условия «прилипания» в разностной форме принимают вид
(Ух)кгт-~(Ух)кгт+1, / = 0,1,...,^; (ух)к,=-(ух)к ■
(vz)o,; =<vztr j = 0X-,N2, {vz)4m
\k
x )1,j; к
'ZJ l,j> \k
(Уг)т.1 =~(yz)k
~ (vz)i,W2+l = z)N\+\,j'>
(Vх ) ОТ
/ = 0,1,...,^2. (19)
Предложенная разностная схема - полунеявная; часть слагаемых в разностных уравнениях берется на к-м временном слое, часть - на (к+1)-м.
Алгоритм вычислений по схеме (15)-(19).
1) Во всех входных функциях (Е, g1, g2 и других) и начальных условиях проводится переход от переменных (х,г) к (х0 по формулам (7), (8).
2) Задается число узлов сетки по пространственным переменным х°, г" (формулы (14)), шаг по времени А/.
3) На текущем временном слое определяются возвышение свободной поверхности жидкости по формулам (18), преобразование координат по формулам (7), (8), сеточные функции
öz
(л о\к dz
dz1
( я о\ ох
Ji,j
dx
Л,
dxL
Ji,i
v St ,
( л o\k dz
dx
f Я2 0\k О z
dxz
fЯ2 0 \k О z
JUj
dxdz
/ = 1,2,...,М1;к = 1,2,....
В (17) значения функции ух берутся на (к+1)-м слое по времени, уг - на к-м слое, в (16) значения функций V,.. - на (Л'+1 )-м слое, плотности р- на к-м.
Кинематические условия расписываются по явной схеме метода конечных разностей, т.е.
Производные функции £ (х,0 по переменным х и / находятся по формулам (16), (18).
4) Разностные уравнения движения (15) ( I, ] = 1, 2, ..., N (N2)) условия (17) (I = 1, 2, ..., Щ), (19) образуют систему (М+2)(Щ2+2) линейных разностных уравнений с (М+2)(Щ2+2) неизвестными для определения компоненты вектора скорости ух на (к+1)-м
слое по времени.
5) Из разностных уравнений движения по оси Oz, условий (16) ( 1= 1, 2, ..., Щ1), (19) определяется сеточная функция на (к+1)-м слое по времени.
6) Разностные уравнения, полученные из (9) и рассматриваемые внутри области и на ее границах, служат для определения плотности жидкости р на (&+1)-м слое по времени.
При определении функций . . р на очередном слое по времени решаются системы (Л/Г1+2)(Л'2+2) линейных разностных уравнений с (М+2)(Щ2+2) неизвестными.
7) Переходим к следующему шагу по времени (к = к + 1) и, если не достигнут окончательный момент времени, то возвращаемся к пункту 3 алгоритма. При расчете на первом временном слое используются начальные условия задачи.
к
+
k
0
x
к
к
+
k
+
k
k
k
к
k
2
2
k
При построении разностной схемы использовались 9-, 6- и 4-точечные шаблоны. Получающиеся в пунктах 4-6 алгоритма системы линейных уравнений решаются модифицированным методом Гаусса с выбором главного элемента по строкам или одним из итерационных методов в зависимости от размера матриц. Благодаря расщеплению исходной системы разностных уравнений на три независимые подсистемы ресурсов современных ЭВМ достаточно для устойчивого счета по предложенной разностной схеме. Программа, реализующая метод конечных разностей, написана в среде обработки Ое1рЫ-6.
Условия устойчивости предложенной разностной схемы выводятся аналогично [7].
Предложена новая разностная схема для решения нелинейных уравнений Навье-Стокса сжимаемой жидкости. К ее достоинствам относятся не слишком жесткие ограничения на шаг по времени в разностной схеме, проведение расчетов на равномерной сетке, расщепление системы уравнений Навье-Стокса на три системы разностных уравнений, решаемых независимо друг от друга. Разностная схема была оттестирована на модельных примерах, для которых известны аналитические решения: движение поршня
бесконечной длины в вязкой жидкости и движение безграничной жидкости со свободной поверхностью при различных нагрузках на ней. Результаты, полученные аналитически и численно, разнятся менее чем на 8 % на промежутке времени до 5 сут. Разработанная схема метода конечных разностей позволяет проводить численное исследование нелинейных уравнений Навье-Стокса с учетом рельефа дна и возвышения свободной поверхности жидкости.
Литература
1. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. М., 1990. Т. 1. Т. 2.
2. Бим Р.М., Уорминг Р.Ф. // Ракетная техника и космонавтика. 1978. Т. 16. № 4. С. 145-156.
3. Годунов С.К., Рябинький В.С. Разностные схемы. Введение в теорию. М., 1977.
4. МаккормакР.В. // Аэрокосмическая техника. 1983. Т. 1. № 4. С. 114-123.
5. Рихтмайер Р.Д., Мортон К. Разностные методы решения краевых задач. М., 1972.
6. Лойцянский Л.Г. Механика жидкости и газа. М., 1970.
7. Усов А.Б. // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2007. № 5. С. 15-18.
Ростовский государственный университет_30 ноября 2006 г.
