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

Вычислительные технологии

Том 11, № 5, 2006

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВЕТРОВЫХ ТЕЧЕНИЙ В СТРАТИФИЦИРОВАННЫХ ВОДОЕМАХ МЕТОДОМ РАСЩЕПЛЕНИЯ. ГИДРОСТАТИЧЕСКОЕ ПРИБЛИЖЕНИЕ*

В.М. Белолипецкий, П. В. Белолипецкий Институт вычислительного моделирования СО РАН, Красноярск, Россия

e-mail: belolip@icm.krasn.ru

Considered here is a numeric algorithm for solving the equations of an incompressible stratified fluid in the hydrostatic approximation, which is based on the splitting method. Results of calculations based on hydrostatic and none hydrostatic models are compared.

Введение

Рассматривается численный алгоритм для исследования ветровых течений в замкнутых водоемах, основанный на уравнениях движения стратифицированной жидкости в приближениях Буссинеска, гидростатики, "твердой крышки", на уравнениях конвекции-диффузии для температуры и солености, использующий метод расщепления по физическим процессам, метод конечных элементов и схему с разностями против потока. Выводятся граничные условия для этапов расщепления. Построено аналитическое решение стационарной задачи для прямоугольной области, которое использовалось в качестве теста. Приводятся результаты сравнения расчетов двумерных в вертикальной плоскости течений по гидростатической и негидростатической моделям.

1. Постановка задачи.

Исходные уравнения и граничные условия

Математическая модель ветровых течений в стратифицированных водоемах (негидростатическая модель) основывается на уравнениях турбулентных течений неоднородной жидкости в приближениях Буссинеска [1]:

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований — ККФН (грант № 05-01-97700-р-енисей), РФФИ-НВО (грант № 05-05-8902 НВО), Министерства образования и науки РФ и Американского фонда гражданских исследований и развития (грант № RUX0-002-KR-06), междисциплинарного интеграционного проекта СО РАН № 24.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

dU 1 dp д dU д dU д dU

+ IV;

(1)

dt р0 дх дх x дх ду у ду dz z dz

(V 1 др д дУ д тдУ д „ дУ 1ТТ

— =--7Г + 7Г+ ТТКу1Г + 7Г- Ш5

ро ду дх дх ду ду дх дх

(Ъ ¿Ш

1 др ро дх

ддШ д

дШ

--^ + 7Г Кх1Г + 7Г Ку^~ + — К— + - д;

дх дх ду ду

д_

д х^

дШ

дх

ро

ди дУ дШ

дх ду дх

0;

(Т — АК дТ АК дТ дТк дТ;

(Ъ дх хТ дх ду уТ ду дх хТ дх ' (С д дС д дС д дС

—¡Т — ТТ~ КхСт;--+ --+ 5

аъ дх дх ду ду дх дх р — Ро(1 +Рг(Т,С)).

(2)

(3)

(4)

(5)

(6)

Здесь и, V, Ш — составляющие скорости течения воды в направлениях Ох, Оу, Ох соответственно, ось Ох направлена на восток, ось Оу — на север, ось Ох направлена вниз; Ъ — время; р — плотность воды; р0 — характерное значение плотности воды; р — давление; д — ускорение свободного падения; Т — температура воды; С — соленость воды; Кх, Ку, Кх, Кхт, Кут, Кхт, Кхс, Кус, Кхс — коэффициенты турбулентного обмена; I —

( д тт д т, д д г

параметр Кориолиса, —- — ——+ и ——+ V ——+ Ш —--субстанциональная производная.

аЪ дЪ дх ду дх

Уравнения (1)-(6) дополняются начальными и граничными условиями. Начальные условия

и — 0, V — 0, Ш — 0, Т — То(х), С — Со(х). Граничные условия на водной поверхности х — я(Ъ,х,у)

К. 8~и

дх

Тх г, дУ _ Ту --, — —--,

ро

дУ_

х Т

КгТ^Т — -

х

ро Г

1 п

сРро

Ш — ^ + и* + У^я,

Ъ х у

р — Ра

К

хС'

дС

х

(7)

~Гг

— -ГС

Граничные условия на твердой поверхности:

и — 0, V — 0, Ш — 0,

, N ^ дТ . дТ . . ^ дТ Рдн 008(п, х) КхТдх + С0в(п, у) КуТду + 008(п, х) КхТдх — ср ■■

С С С

008(и, х) КхС-ТТ- + 0О8(и, у) КхС+ 0ОВ(п, х) — РдС.

х у х

(8)

Здесь тх, ту — составляющие напряжения трения ветра; гра — давление на водной поверхности; Гп — полный поток тепла через свободную поверхность; Гдн — теплообмен с ложем водоема; гс — поток соли через свободную поверхность; ГдС — массообмен с ложем водоема, п — внешняя нормаль к поверхности; То(х), Со(х) — заданные функции.

Для движений, горизонтальный масштаб которых много больше вертикального, применимо гидростатическое приближение. Вместо уравнения (3) рассматривается уравнение статики [1, 2]:

р

— — дро (1 + р1). (9)

Упрощение уравнений

Проинтегрируем уравнение (9) от q до z с учетом граничного условия на водной поверхности:

z z q

p = gj Po(1 + Pi) dz+pa = gpoj (1 + Pi) dz - gpoj (1 + pi) dz + pa ^

q 0 0

z

~ gpoj (1 + Pi) dz - gpoS + pa-0

Для pa = const из предыдущего соотношения следует

z

—V2p = g / V2pidz - gV2q, (10)

Po

o

г = / _d_ д 2 \дХ дуy

Проинтегрируем уравнение неразрывности по z от q до H:

н

Г idU dV\

J \дХ + -8y)dz + W|z=H - Wlz=q = 0- (11)

q

Так как выполняются условия на свободной поверхности (7) и прилипания на дне, из (11) получаем

н н

/ Udz+-^ [Vdz - % = 0. дх J д у J дt

qq

В приближении "твердой крышки" это соотношение принимает вид

н н

UUiz+hiViz=0- (12)

oo

Итак, в приближениях гидростатики (10) и "твердой крышки" (12) для исследования стратифицированных течений в водоемах получена следующая система уравнений:

z

dU = &U 9U &U lv dl ?dpidz (13)

dt дх x дх ду y ду dz z dz gдх j дх Z;

o

z

dv = ov + dlk dV + ^kov_iU + di_ [ дpi

dt дх x дх ду y ду dz z дz ^ду j ду

o

н н

^ f Udz+д I Vdz = 0; (15)

д х ду

oo

-Ш

дг

ди — У_

дх ду

(16)

Особенностью гидростатического приближения является то, что для определения вертикальной составляющей скорости течения достаточно одного граничного условия (либо на дне, либо на водной поверхности). Система уравнений (13)-(16) дополняется граничными и начальными условиями: на водной поверхности (г = 0)

ди

Тх

г дг

на дне и боковой границе начальными условиями

К — = - —, К

д V

Ту

Ро дг

и = 0, V = 0; и = 0, V = 0, Ш = 0.

, Ш = 0;

Ро

(17)

(18) (19)

2. Численный алгоритм

Сформулированную задачу решаем методом расщепления [3], реализованным в работах

[4, 5].

Этап 1

и - ип д ип д ип ди + ип -— + Уп -— + Шп

А г V - Vп А г

дх д Vп

+ ип + V'

дх

ду дV -у

+ Шп

дг дV п дг

= IVп,

-нип.

Этап 2

и * - и д ^ ди * д ^ ди * д ^ ди * ^ —К* — + —Ку— + - Р?

А г

дх дх ду ду дг дг

V* - V д дV* д „ дV* д „ дV* -= — Кх^— + — Ку+ — К*-

А г

дх дх ду ду дг дг

г г

РГ = 9 [^г, Р2п = д [ —Рп¿г.

- Рп,

-х

о

о

-у

Этап 3

- (Н~) + -—- (н—^-

-х \ -х) -у \ -у

1

н н

л+\7Т [ и*йг + V д А г \ — у и

оо

Этап 4

ип+1 = и* + дА г-^, Vn+1 = V * + дА г-^;

-х -у

г

Шп+1 = (—ип+1 + —¿г.

х у

о

(20)

(21)

(22)

(23)

(24)

Здесь А г — шаг по времени. Уравнение (22) получается в результате подстановки соотношений (23) в уравнение (12). Начальными условиями для этапа 1 являются условия (19).

Вывод граничных условий для этапов схемы расщепления

Водоем аппроксимируется сеточной областью в виде прямоугольных параллелепипедов. В результате граница области состоит из прямоугольных площадок, расположенных в горизонтальных или вертикальных плоскостях. Вертикальное сечение водоема представляется в виде горизонтальной линии (невозмущенная поверхность воды) и ломаной линии, отрезки которой параллельны соответствующим осям координат. Контур границы водоема аппроксимируется ломаной линией, отрезки которой параллельны оси координат Ох или

Oy.

Из соотношений (17) и (23) находятся условия на водной поверхности:

dU * dz

т- ^ d V *

—-, Kz

Po

dz

- Tv, W = 0. P0

(25)

На горизонтальных участках твердой границы z = Hi = const по условию прилипания Un+1 = 0, Vn+1 = 0, тогда из (23) следует

U * = -gA tp-, V* = -gA t

dx dy

(26)

из (13) и (14) находим

д _ dUn

g dx dz z dz

дя

dr dVn

+ Fl , g dy dzKz dz

+ Fn,

следовательно, на горизонтальных площадках границы

dr dUn

U* = Atl—Kz _ d z dz

- F1n

V*

At [ —K

ddVn

dz"dz

Fn

(27)

Аналогично выводятся условия на вертикальных площадках границы:

U *|

A t — K.

dr dUn

d x - dx

Fn

d dUn

U U = At[^Kv^ - Fn

dy dy

V *\

V *\

A t —K,

dr dVn

dx - dx

F

= A 'Ф'' ж- F

(28)

Граничные условия для уравнения (22) определяются из условий на границе контура водоема (18) и соотношений (23):

dx

1 U *l --— U x=-i , ——

g A t z=o dy

v=vj

JKt V *\ V=0f •

(29)

Обычно коэффициенты турбулентного обмена Кх, Ку, КхТ, КуТ, КхС, КуС считают постоянными. Коэффициент вертикального турбулентного обмена Кх вычисляется по формуле Прандтля — Обухова [6]:

Kz

(0.05h)2vB + Kn K

1 min

при B > 0, при B < 0,

где

dU 2 dV

Bчiz) чт,)-%1;

g dP

x=x

x=x

2

к — глубина квазиоднородного слоя, определяемая по первой от поверхности расчетной точке, в которой выполняется условие

(0.05г^)УВ\г=гк < Кт1Е

Кт-т — фоновое значение коэффициента вертикального турбулентного обмена.

После определения составляющих вектора скорости течения вычисляются температура, соленость и плотность воды.

Замечание 1. Уравнение (22) для определения я можно привести к другому виду. Выполним дифференцирование по горизонтальным координатам:

дА г

д\ —2_Л д_Н-к -н -я

х2 у2 х х у у

н

о

_-Н *. _—Нт/*|

дх и 1г=н -у У 1г=н,

отсюда с учетом (23) следует

-и* -V*'. ,

-7,--+ -7— I аг-

х у

н

-2я + —2я2 = К-и* + —Й ¿г. (30)

дх2 -у2 дН А г } \ -х -у

о

Замечание 2. Задача Неймана (30), (29) имеет единственное (с точностью до постоянной) решение, если выполняется условие [7]

+ / = °- (31)

о

Здесь Г — граница области С; п — нормаль к границе Г;

н

= 1 [{-и * -V *\

1 = дНКг] \ех + ~—у)

о

Так как по формуле Гаусса — Остроградского

!{аа = дж^ДЙ ¿г = дА;Дп* ¿Г

о о \г / г

1 -V*

и по (29)

—К

-п г дАг

условие (31) выполнено.

Замечание 3. Из соотношений (23), (24), (30) получаем формулу для вычисления вертикальной составляющей скорости Шп+1:

^ н

ттт+1 Г (-и* дV*\ , г Г (-и* дV*\ ,

Ш = - А + 1у)аг+н]Ы + -ду)аг- (32)

оо

В гидростатическом приближении вертикальная составляющая скорости Жп+1 обращается в нуль только на водной поверхности (либо на дне). На дне (или на водной поверхности) Жп+1, вообще говоря, не равна нулю. Можно из предварительно продифференцированного по вертикали уравнения неразрывности построить решение для Жп+1, удовлетворяющее двум граничным условиям, но тогда уравнение неразрывности внутри области будет выполняться приближенно [2]. В предлагаемом алгоритме согласно (32) Жп+1и = 0 и Жп+%=н = 0.

Итак, численный алгоритм для исследования трехмерных ветровых течений в стратифицированных водоемах представляется соотношениями (20), (21), (23), (25), (27)-(30), (32).

Замечание 4. В качестве граничных условий на твердых участках границы для уравнений (21) можно использовать соотношения (26) с известными (на п-м временном шаге) правыми частями. Для уравнения (22) граничные условия можно получить из проинтегрированных по г от 0 до Н уравнений (13), (14) с учетом условий на границе берегового контура водоема (18) и соотношений (23).

3. Численный алгоритм для двумерных в вертикальной плоскости течений

Численный алгоритм для двумерных в вертикальной плоскости течений является частным случаем алгоритма для пространственных задач. Этап 1

Этап 2

и - ип ТТ„ д ип ттгп дип

+ ип—- + =0.

д г

дх

дг

(33)

и * - Ц = дЦ* + _ТК дии^_ рп

Д г дх х дх дг г дг 1 '

т=>Г£

К

дЦ *

д г

--, Ц* 1*=н = д М«тК

г=0

Ро

ВВЦ п

дг дг

-К

(34)

Этап 3

* д дЦп п

Ц |х=х,- ДМ дхКх дх ^1

дх \ дх

д Д г дх ^

н

Ц *йг,

(35)

Этап 4

дх

х=хо

1 Ц *1 ^ --— Ц х=хо , ——

д Дг *=0 дх

х=х£

-1- Ц * I х=хь

д Д г 4=0

г н

Цп+1 = Ц* + дДЖп+1 = - [ ^¿г + г [ ^¿г.

дх д х Н дх

о о

Для двумерного варианта этапы 3 и 4 упрощаются. Проинтегрируем уравнение относительно я по х и с учетом граничных условий для непроточного водоема из (36) получим этап 3а

н г н

ип+ = и* - }_! иШп+ = ^ + (37)

0 0 0 Таким образом, в данном варианте исключается этап численного решения задачи (35), компоненты вектора скорости определяются в результате последовательного выполнения этапов 1, 2 и 3а.

3.1. Аналитическое решение двумерной стационарной задачи

Рассмотрим двумерную стационарную задачу:

ТТ-и тттдШ д2и -я ТЛ, ,

и^Г + = 7ГТ + 9д- - Р1(х,г),

х г г2 х

—и+—ш=0 _

х г х

иаг = 0,

-и

г

г=0

-32(2 + а)(£) (1 - £) , и\г=1 = 0

и\х=0 = 0, и\х=ь = 0, 0 < х < Ь, 0 < г < 1, а > 1. Эта задача имеет аналитическое решение:

и = В0(х) А^г), Ш = -В1(х) А0(г),

(38)

(39)

где

В0(х) = 1б (хГ (1 - Ь)\ в,(х) = | (Ь) (1 - Ь) (1 - 2

А0(г) = г(1 - г)2(1 - аг), А^г) = (1 - г) [1 - (3 + 2а)г + 4аг2]. Правая часть первого уравнения (38) вычисляется по формуле

д—х - Р1 = А1В0В1 - А0В0В1 ¿А- - В0 ¿¿А. -х йг йг2

(40)

(41)

Решение задачи (38) численно будем искать методом установления с помощью алгоритма (33), (34), (37). Предварительно необходимо определить функцию Р1(х,г). Проинтегрируем первое уравнение (38) по г от 0 до 1:

-и

г

г=1

-и

г

г=0

1 1 1

я и и

+ д-^ - Р ¿г = IU—йг+ W—йг.

х х г

Тогда из первого уравнения (38) получаем интегродифференциальное уравнение для определения Р1(х,г):

2и и Р = Р ¿г + — - —

г2 г

г=1

—и

г

г=0

1 1

и и и и

- и---Ш— - и—¿г - Ш—^.

х г х г

1

Решение этого уравнения имеет вид

Fi(x,z) = ^ bn

(42)

где

bi

Ьз

b5

- 24aB0 + 2(2 + a)B0B1, (8 + 16a + 8a2)B0B1, 6a(1 + 2a)B0B1,

b2 b4

be

-(8 +8a + 2a2)B0B1, -(3 + 16a + 14a2)B0B1, -4a2 B0B1.

В частности, F1\z=0 = 0, F1\

'1\z=1

B0(B1 - 24a).

3.2. Тестовая задача для нестационарных течений

Рассмотрим упрощенную задачу для ветровых течений однородной жидкости в прямоугольной области [0 < x < L, 0 < z < H] при Kx = const, Kz = const, условиях непротекания на вертикальных стенках и условии проскальзывания на дне. Для нахождения U и я получаем задачу

dU ^ d2U ^ d2U 1Я

— = + Kz— + g—, (43)

dt

dx2 н

dz2 dx

Udz = 0;

(44)

z = 0: Kzd-U = -т; dz

dU

z = H : Kz— = 0;

dz

x = 0, x = L : U = 0.

(45)

(46)

(47)

Проинтегрируем уравнение (43) по z от нуля до H:

н н

d Г d2 Г dU

— Udz = K^—— Udz + Kz — dt J dx2 J dz

00

Отсюда в силу условий (44)-(46) находим

z=H

- к. dU

z dz

м1я gH^~ = -т•

dx

Для определения U получаем уравнение

dU = K d!U K d2U т

dt x dx2 z dz2 H'

z=0

+gH£■

Таким образом, задача (43)-(47) свелась к задаче (45)-(48).

n

z

4. Результаты расчетов

Стационарная задача. При а = 1 аналитическое решение (39) описывает течение неоднородной жидкости с одной циркуляционной зоной, а для а > 1 — с двумя циркуляционными зонами. Решение задачи (38) предлагаемым методом позволяет проверить правильность постановки граничных условий для этапов расщепления. На рис. 1, а приведен пример численного решения стационарной задачи методом установления на прямоугольной сетке 24 х 10 для а = 1.8, L = 25. На грубой расчетной сетке численное решение существенно отличается от аналитического. Однако уже на сетке 32 х 20 результат расчета графически совпадает с аналитическим решением (рис. 1, б).

Нестационарная задача. Для нестационарной задачи (43)-(47) проведены расчеты как для т = const, так и для т = 0.0125 | Wa| Wa, где Wa = 10sin (At), Л = 20п. Численные решения задачи (43)-(47) методом расщепления графически совпадали с численными ре-

а б

0.311-0.180 -0.049 0.082 0.213 0.344 0.476 0.607 0.738 0.869 1.000 -0.295 -0.165 -0.035 0.095 0.225 0.355 0.485 0.616 0.746 0.876 1.006

Рис. 1. Решения стационарной задачи (38) (Ь = 25, х = Ь/2): сплошная линия — аналитическое решение, штриховая — численное решение; а — сетка 24 х 10; б — сетка 32 х 20.

-4.070 -3.072 -2.074 -1.076 -0.077 0.921 1.919 2.917 3.916 4.914 5.912

Рис. 2. Профили горизонтальной составляющей скорости течения в глубоководной области водоема: сплошная линия — негидростатическая модель, штриховая — гидростатическая модель.

шениями задачи (45)-(48) (на сетках 30 х 30 при одинаковых шагах по времени). Расчеты показали, что использование в методе расщепления краевых условий с первым порядком приводит к снижению устойчивости численной схемы.

Сравнение гидростатической и негидростатической моделей. Сравнение гидростатической и полной моделей гидротермики глубоких водоемов в рамках проекционных методов для двумерных течений выполнено в работе [8]. В настоящей работе проведена серия расчетов двумерных в вертикальной плоскости ветровых течений в непроточном водоеме по гидростатической и негидростатической моделям. На рис. 2 представлены рассчитанные профили горизонтальной составляющей скорости в глубоководной области водоема (глубина водоема 20 м, длина — 4800 м, вода пресная) для ветра 5 м/с и вертикального распределения температуры Т0(г):

Для ветра 18 м/с формируется одна циркуляционная зона. Переход от двухциркуляци-онной картины течения к одноциркуляционной происходит при близких значениях ветра по обеим моделям (для рассмотренного примера при ветре 17.5 м/с).

Сравнение результатов численного моделирования двумерных течений в вертикальной плоскости по гидростатической и негидростатической моделям показало, что обе модели дают достаточно близкие значения скоростей ветровых течений (отличия максимальных значений скоростей не превышают 10%). Поэтому для исследования гидротермических режимов слабопроточных водоемов можно использовать гидростатическое приближение.

Список литературы

[1] Филатов Н.Н. Динамика озер. Л.: Гидрометеоиздат, 1983. 166 с.

[2] Саркисян А.С. Моделирование динамики океана. СПб.: Гидрометеоиздат, 1991. 292 с.

[3] БЕлоцЕрковский О.М. Численное моделирование в механике сплошных сред. М.: Наука,

[4] БЕлолипЕцкий П.В. Численное моделирование ветровых течений в стратифицированных водоемах методом конечных элементов // Тр. VI конф. "Динамика и термика рек, водохранилищ и прибрежной зоны морей". М., 2004. С. 26-29.

[5] Белолипецкий П.В. Численное моделирование двумерных в вертикальной плоскости ветровых течений в стратифицированных водоемах методом расщепления // Вычисл. техноло-

гии. 2005. Т. 10, № 5. С. 19-28. [6] Математические модели циркуляции в океане / Г.И. Марчук и др. Новосибирск: Наука,

[7] Михайлов В.П. Дифференциальные уравнения в частных производных. М.: Наука. 1976.

[8] Бочаров О.Б., Овчинникова Т.Э. Численное моделирование явления термобара в озере Байкал // Вычисл. технологии. 1996. Т. 1, № 3. С. 21-28.

1984. 520 с.

1980. 288 с.

391 с.

Поступила в редакцию 15 ноября 2005 г., в переработанном виде — 20 февраля 2006 г.