Моделирование трехмерной гидродинамики мелкого моря в приближении возмущенной плотности Текст научной статьи по специальности «Физика»

УДК 519.6

В.С. Васильев

МОДЕЛИРОВАНИЕ ТРЕХМЕРНОЙ ГИДРОДИНАМИКИ МЕЛКОГО МОРЯ В ПРИБЛИЖЕНИИ ВОЗМУЩЕННОЙ ПЛОТНОСТИ

Безразмерные переменные. Введем характерные масштабы lp , Lt, lp , L(p , lq , l^ , l^ , Lu, Lv, Lw, Lx, Ly, Lz плотности p, времени t, давления p, гравитационного потенциала ф , угловой скорости вращения Земли q , первого и второго коэффициентов вязкости ц и ц', компонент u, v и w вектора скорости v вдоль координатных направлений х, y и z соответственно. Физические постоянные (g и, возможно, ц и ц') примут в выбранных масштабах соответствующие числовые значения. Соотношение масштабов выберем таким, чтобы числа Струхала Sh х = lxl-1l-1, Sh y = LyL-1L-1 и Sh z = llL 1 были равны между

собой Shx = Shy = Shz = Sh и равны единице Sh = 1. Потребуем также выполнения

для чисел Эйлера Euх = LpL-V, Euy = LpL-^L-2, Euz = LpLp1L-2 и Фруда Frx = Lф1L2u ,

Fry = l^lV , Frz = l^lI соотношений Euх = Fr-1, Euy = Fr-1, Euz = Fr;1. Тогда система

уравнений вязкой сжимаемой жидкости [1] в безразмерных переменных во вращающейся с постоянной угловой скоростью системе координат примет вид

p' + div(pv )= 0, (1)

(pv)t' + i a div(pva v )= -R2 (grad(p-X'div v)+p grad ф-i a div(dv/ йха))+

+ ia div(xR2 gradva)+ 2R П R-1pv . (2) где v = iavа (парные строчные греческие буквы традиционно для тензорного исчисления означают суммирование от 1 до размерности пространства n = з); ij, vj, 1 < j < n - декартовы единичные орты и компоненты вектора v вдоль соответствующих координатных направлений X, 1 < j < n соответственно; r = díagL-1, L~y,l-1}; x и X' - коэффициенты турбулентного обмена импульсом; п -

матрица Кориолисова ускорения.

Физические компоненты вектора скорости содержит вектор R-1v .

Однородное поле. В однородном поле тяжести gradф = ^ (g - ускорение свободного падения) появляется выделенное направление (вертикальное). При выборе направлений х, y, z на восток, на север и вертикально вверх соответственно матрица п имеет вид

0 sin S - cos Sl

П = Q

v

- sin S 0 cos S 0

(3)

где ^ - северная широта.

Если потребовать равенства единице чисел Эйлера и Фруда в уравнениях для горизонтальных компонент вектора скорости Ейх = Ейу = Ей = 1, Ргх = Ргу = Рг = 1,

то Ейг = я2Еи = я2, Ргг = я~2Рг = я~2 и я = а1аё{1,1,я], где я = ьк/ь2 - отношение масштабов по горизонтали ьк = ьх = ьу и вертикали ьг.

0

0

Для получения уравнения для возвышения свободной поверхности е = е((, х,у) проинтегрируем уравнения системы (1), (2) по вертикали от поверхности дна ё = ё((, х, у) до е

( e A (

+ div

Jpdz

V d Jt

(2)

Jpvdz + Pe (v e, n e )-e^ + Pd (v d , n d )+ d't)= 0,

(4)

Jpvdz

A (

+ ia div,

V d

(2)

H

Vd A

Jpva vdz +Pe v e ( ,n e )-e;)+Pd V d (v d >П d )+ dt') =

Vd

(e A

grad(2) Jndz + R 2

V d J

-п-ene -nd n d + gJpdz + i a div(2) Jxdv/ dxa dz

\A

yj

+ ia div(2)

|XR2 gradvadz + т(R)ne + т(R)nd + 2R П R-1 Jpvdz , (5)

V d

i e Л

где J ft'dz = J fdz - fee't+ fdd't \ J diV a dz = diV(2) J a dz +(ae, n e ) + (ad , n d ),

( e A

a

V d J

( e A

diV (2)( + ayj + azk )= (ax )x +(ay ^ ; J grad(f)dz = grad(2 ) J fdz

+ fen e + fdn d ,

V d J

grad(2) f = fXi + fy j+0k ; ne =-eXi - ey j+k, n d = d'x i+d'y j - k - внешние (по отношению к столбу жидкости) нормали; п = p -X'div v ;

2К x(My + vX ) x(r 2u'z + w'x )h

t(r ) =

x(uy + vx )

2Xv'„

Ш V + w'„

2и'г + w'x) фЧ + w'y) 2ХЯ2<

Мелкое море. В мелком море я составляет 104... 105, то есть я >> 1. В этом случае наличие w в уравнениях (2) для и и V, вносимое слагаемыми R21 аdiv(^5^/5xа), 2Rй R-1ру, не существенно из-за коэффициентов я~2 и яЛ В то же время наличие и и V с коэффициентами я2 и я1 в уравнении для w, вносимое этими же слагаемыми, ведет к численной неустойчивости. Причем сгущение сетки в горизонтальных направлениях (уменьшение я ) проблему не устраняет. Причина в масштабах процессов, которые таковы, что различие в поле скоростей на расстоянии 1 м по вертикали может быть существеннее различия на расстоянии 10 км по горизонтали (длинноволновая гидродинамика).

Наличие компоненты V в уравнении для и и компоненты и в уравнении для V, вносимого слагаемым R 21 а &у(хду/дха), хотя и не имеет коэффициентов я2 и я1, но оставляет возможность «перекачивания» энергии между компонентами и и V (строго диссипативным является не каждое уравнение в отдельности, а всего лишь система в целом). Вследствие этих подсеточных процессов развивается еще один тип неустойчивости (или, по крайней мере, замедление сходимости).

Гидростатическое приближение. Если в давлении п выделяется слагаемое а, не зависящее от координаты г, то содержащееся в (5) выражение

i e \

E = !

grad(2) Jndz + R2 ^ne +^dnd - gJpdz

V d J

преобразуется к виду

+

t

E = grad(2) Ha + J (n - a)dz V d

где h = e - d - полная глубина.

В гидростатическом приближении

J(n- a)dz + R2 пеПе nd - gJpdz

d

p = a + J pgdz .

P'z + pg = 0,

где a((, x, y) - атмосферное давление, и

f ее Л f е Л е

E = H grad(2) a + g grad^) JdzJpds + g Jpdz grad(2) d - R2 Jgrad(X'div v)dz

V d z J V d J d

е

В случае постоянной плотности р = const Jpgds = pg(е - z) и

E =

е

H (grad(2) a + pg grad(2) е) - R2 J grad(X' div v )dz .

В соленоидальном поле div v = 0

E = H(grad(2) a + Pg grad(2) е) является двумерным и не содержит r 2.

С другой стороны, если r2ia div(Xdv/dxa)= r2 grad(X div v) (например, в случае X = const), то (2) преобразуется к виду

(pv)' + iа div(pvаv - XR2 grad vа) = -R2 (grad(p - (X + X')div v) - pg) + 2R fi R-1pv ,

а в соленоидальном поле - к виду

(pv)' + iа div(pvаv - XR2 gradvа)= -R2(gradp - pg) + 2R fi R-1pv . (6)

В этом случае слагаемые, вызывающие численную неустойчивость, содержит только 2R fi R-1pv .

Приближение возмущенной плотности. В приближении возмущенной плотности

p = po +8 ,

где |s/p0| << 1, p0 - некоторое среднее значение плотности,

P'z +pg - 0,

но наличие коэффициента r 2 может делать вклад выражения p'z + pg сравнимым с вкладом остальных слагаемых.

С другой стороны, для устойчивого разрешения системы (1), (6) матрица (3) вынужденно преобразуется к виду (длинноволновая гидродинамика)

f 0 sin Q 0Л fi' =Q - sin Q 0 0

v 0 0 0J

а (6) - к виду

(pv)' + ia div(pvav - XR2 grad va)- 2fi' pv = -R2 (grad p - pg).

(7)

Для разрешения системы уравнений гидродинамики мелкого моря (1), (7) (вязкой сжимаемой в приближении возмущенной плотности жидкости) предлагается следующая модификация МАС-техники [1]

рТ = ру + т(а div(xR2 gradVа - руау)+ 2П'ру), (8)

~ = с +тdiv((grad~-су), Т = Т +тdiv(кgradТ-Ту), р = р(Т,~, (9)

Н = Н + Т2ШУ

(2)

(н^Ц2)(ро 1а + я«г))+ Гн ,

Т2а1у(я2 (( р -р

ру = ( ( гС Р~ + ((г - 2)Ру )С ) - тК2 (( р - р§)"

(10)

гр, (11)

Л (12)

где р~ - промежуточное, не удовлетворяющее уравнению неразрывности поле динамической скорости; т - шаг по времени; с - концентрация примеси; в -коэффициент диффузии рассматриваемой примеси; т -температуры среды; к -коэффициент теплопроводности (в правой части опущено слагаемое с тепловыделением из-за трения, так как для рассматриваемой задачи оно не существенно);

-(е ) + ШУ(2) т(юе + Vй)-р-11ру ёг

ее ^^ е , . ) е

|руёг = |руёг + т 1а а1у(2)1 J(xR2grad уа-руау)г + т(я П + т^ ^ + 2П' |руёг

)

Л

юе и юё - интенсивности испарения-осадков (юе > о - испарение) и фильтрации (юё > о - в грунт) в кинематических условиях

е!+®е = (уе , Пе ), + (у ё , )=юё ;

гр = р - ( I-1 (г^ Р + ((г - г)р)?' ) + т div((гZ ^ (г?рр + ( - г^)) )] •

Введение криволинейных координат (§, п, с), подвижных вслед за движениями свободной поверхности (не среды), изложено в [2].

Уравнение (10) получается подстановкой уравнения (5) в форме

е е / ~ |руёг = | р~ёг - т(( grad(2)(a + +

ё ё в уравнение (4) в форме

л )

+ р( е + ®ё )

ёшё

))

| рёг = | рёг - т div(2) | руёг

ё ё ( (ё )

Аналогично, уравнение (11) получается подстановкой (12) в сеточное уравнение (1)

р = ( I-1 (г?р + ((г - г)р)с )- т^(ру) •

Тем самым, на сеточном уровне обеспечивается выполнение уравнения неразрывности.

Возможны модификации (8)-(12), когда с учетом нового положения узлов рассчитываются с и г, а р в правой части гр не отрабатывается по р, а

рассчитывается из алгебраического уравнения состояния р = р(с, Т). Отличие (8)-(12) от [3] в том, что не требуется закон стратификации плотности. В отличие от [4] система (8)-(12) содержит одновременно с и т, и расширяется на большее количество параметров. В отличие от [5] модель (8)-(12) рассчитана на существенное непостоянство Н и на колебания уровня е , сравнимые с Н .

Ги = X

н

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах: Пер. с англ. - М.: Мир, 1991.

2. Васильев В.С. Энергетически нейтральная строго диссипативная аппроксимация двумерной гидродинамической системы на подвижной криволинейной сетке в вертикальном поле // Известия ТРТУ, 2003, № 5, с. 184-189.

3. Белоцерковский С.А., Гущин В.А., Коньшин А.Г. Метод расщепления для исследования течений стратифицированной жидкости со свободной поверхностью // Ж. вычисл. матем. и матем. физики, 1987, т. 27, № 4, с. 594609.

4. Иевлев В.М. Численное моделирование турбулентных течений. М.: Наука, 1990. - 216 с.

5. Марчук Г.И., Кочергин В.П. Саркисян А.С. и др. Математические модели циркуляции в океане. - Новосибирск: 1980, - 285 с.