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

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

Том 10, № 5, 2005

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДВУМЕРНЫХ В ВЕРТИКАЛЬНОЙ ПЛОСКОСТИ ВЕТРОВЫХ ТЕЧЕНИЙ В СТРАТИФИЦИРОВАННЫХ ВОДОЕМАХ МЕТОДОМ РАСЩЕПЛЕНИЯ*

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

Красноярск, Россия e-mail: belolip@icm.krasn.ru

A numerical algorithm for solving equations of an uncompressed stratified flow is considered. The algorithm relies on the method of splitting by physical processes, as well as the finite element method and finite difference counter flow scheme. The examples of test calculations are presented.

Введение

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

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

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

ди ТТди Т1 ди 1 др дди дди

"ТГТ + и— + =--+ — Кг— + — кх —,

дъ д х д г р0 дх д г д г дх дх

* Работа выполнена при поддержке Красноярского краевого фонда науки (грант № 12F003M), Российского фонда фундаментальных исследований — ККФН (грант № 05-01-97700-р-енисей). Research described in this publication was made possible in part by Award N. KY-002-X1 of the U.S. Civilian Research & Development Foundation for the Independent States of the Former Soviet Union (CRDF).

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

дШ ттдШ ттгдШ + и—--+ Ш-

дЬ

дх

д г

1 др д тдШ д дШ р

р0 д г д г д г дх дх р0

ди дШ _ 0

дх дг

(1)

(2)

дТ ттдТ ттгдТ д тдТ д дТ + и— + Ш— _ — КхТ— + — КхТ—, дЬ дх дг дг дг дх дх

дС ттдС ВС д тдС д дС "яГ + и— + Ш— _ — К^ + КжС —, дЬ дх дг дг дг дх дх

р _ Ро(1+ Рг(Т,С)).

(3)

Здесь и,Ш — составляющие скорости течения воды в направлениях х и г; ось г направлена вниз; Ь — время; р — плотность воды; р0 — характерная плотность воды; р — давление; 9 — ускорение свободного падения; Т — температура воды; С — соленость воды; Кх, Кхт, Кхс, Кг, Кгт, Кгс — коэффициенты турбулентного обмена. Вводится функция р1:

Р1 _ Р - 9Ро р\^г,

(4)

тогда уравнения (1) приводятся к виду

ди ди ди + и---+ Ш-

дЬ

дх

дг

1 д р 1 д д и д д и д

ро д х д г д г д х д х д х

о

дШ ттдШ ТТГдШ + и—--+ Ш-

дЬ

дх

дг

1 др1 д тдШ д „ дШ

ро д г д г д г д х д х

(1а)

Уравнения (1а), (2), (3) дополняются начальными и граничными условиями. Граничные условия на водной поверхности (г _ 0):

К

Т

К

дТ

дг

ди

дг

--, Ш

ро

0,

1 п

СрРо

К

с

дС

дг

-Р

с,

(5)

РД

СрРо

на дне и на боковой поверхности (г _ Н или х _ х^):

и _0, Ш _0,

дТ дТ

СО$(и, х)КхТ^--+ 008(и, г)КгТ-тг~ _

дх дг

дС дС

дх дг

Здесь и — внешняя нормаль к поверхности; т — напряжение трения ветра; Еп — полный поток тепла через водную поверхность; Рдн — плотность потока тепла на твердой границе области; Рс — поток соли через водную поверхность; Рдс — плотность потока соли на твердой границе области; ср — удельная теплоемкость воды.

до-

(6)

г

Для определения коэффициента вертикального турбулентного обмена применяется формула Прандтля — Обухова [2, 3]:

Кх

II ди\2 д дР I V (ди\2 9дР-^(л

(ди\2 9дР^п

Кт1п при —--— < 0,

\дг ; ро дг

(7)

где Ктп — фоновое значение коэффициента вертикального турбулентного обмена; Н — глубина верхнего квазиоднородного слоя, которая определяется по первой от поверхности расчетной точке, где выполняется условие

(?Е)2 - < 0.

\ дг ) р0 дг

Если значение Н больше половины максимальной глубины водоема, то Н приравнивается этой половине.

Коэффициенты горизонтального турбулентного обмена Кх, Кхт, Кхс считаются постоянными; Кхт = атКг, Кхс = асКг, где ат, ас — турбулентные числа Прандтля. Уравнение состояния для пресной воды [1]

Р1 = -0.68 ■ 10-5(Т - 4)2, ро = 1.0 г/см3. (8)

Уравнение состояния для соленой воды берется в приближении Буссинеска [4]:

р1 = -1.8573 ■ 10-4Т + 0.8 ■ 10-3С, р0 = 1.0000726 г/см3. (9)

Соотношения (8), (9) позволяют оценить различие плотностной стратификации для пресной и соленой воды.

Напряжение трения на водной поверхности определяется по формуле Н. А. Давтян [5]:

т =1.25 ■ 10-6(0.69 + 1.07 ■ 10-3 Ша) Wa2 [г/(см3 ■ с2)], (10)

где Wa — скорость ветра, см/с.

Подобная задача в работе [6] численно решалась для уравнений, записанных относительно переменных функция тока — вихрь, в работе [7] численный алгоритм основывался на так называемом ^-проекционном методе, в работе [2] рассматривалась задача в приближении медленных течений и без учета горизонтального турбулентного обмена.

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

Сечение водоема приближенно представляется в виде горизонтальной линии и кусочно-линейного дна, на котором отрезки соединяют известные значения глубины. На эту схематизацию накладывается прямоугольная сетка, из которой удалены элементы, полностью не лежащие или большая часть которых не лежит внутри области, в результате получается сетка для расчетов. Написана подпрограмма, реализующая построение нерегулярной сетки. Граница водоема аппроксимируется ломаной, каждое звено которой параллельно одной из осей координат. Для лучшего описания течений вблизи берегов водоема большой протяженности используется сетка со сгущением в прибрежной области (рис. 1).

к

А-

\

ч

ч

\

S 1

/

1

ч

ч.

ч

■v /

Ч. /

1

-- /

/

Рис. 1. Схематизация водоема.

Численный алгоритм основан на методе расщепления по физическим процессам [8, 9]. Поля скоростей и давления находятся из решения задачи (1а), (2), (5), (6) по известным значениям плотности. Затем из решения задачи (3), (5), (6) с учетом (8) или (9) определяются температура, соленость и плотность воды.

Вначале на временном интервале Ьп < Ь < Ьп+1 (Ьп+1 — Ьп _ ДЬ) рассчитывается конвективный перенос по явной схеме:

U - Un + Un d-U! + WndU

At

W - Wn

At

dx dz

dWn dWn + un -— + Wn

ox dz

0,

0,

(11)

где ип _ и(Ьп, х, г), Шп _ Ш(Ьп,х,г). Так как (Уи) _ 0 на всей границе (V _ (и,Ш)), для уравнений (11) не требуется граничных условий.

Далее на этом же временном интервале решаются уравнения, полученные с помощью метода расщепления по физическим процессам [8]:

U * - U At

д / dU* dx\ x дх

д / dU* dz\ z dz

dp

dx

1 dz,

W * - W At

JL (K dW*

dx V x dx

d (K dW *

dz \ z dz

+ g;

Api = PL (dU* + dW *

At V dx dz

Un+1 = U * -

At dp1

Wn+1 = w * -

At dp1

(12)

(13)

(14)

po dx} po dz

В рассматриваемых задачах возможны следующие типы границ: горизонтальная и вертикальная твердые стенки, граница раздела вода — атмосфера (водная поверхность). Сформулируем граничные условия для уравнения (13). Пусть x = xi = const — вертикальный участок твердой стенки. Из уравнения движения в проекции на ось x (1а) и условия

z

g

непротекания находится соотношение

1 dP = «к, дЦ. - д[М. d, (15)

р0 дх дх дх J дх

о

Так как на горизонтальных участках границы области (z = const) W = 0, из проекции уравнения движения на ось z (1а) следует условие

1 dpi 0K0Wn

— + 9- (16)

р0 дz дz дz

Рассмотрим вывод граничных условий для уравнений (12). На горизонтальных участках границы области (z = const) U = 0, W = 0, тогда из (14) и проекции уравнения движения на ось х получается соотношение

z

U' = Д 4Kz ^ - Д ^ "'Z, (17)

о

а из (14) и проекции уравнения движения на ось z следует условие

д дWn

W * = Д t—Kz д— + Д tg. (18)

dz dz

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

z

д dUn Г д nn д dWn

U* = Д t-d-Kx^ - Дt9 (р1 "z, W* = Дt-d-K,^ + Дtg- (19)

дх дх дх дх дх о

На водной поверхности при z = 0 для W* ставится условие (18). Из (5), (14) и (16) находится условие для U* при z = 0:

г, (U* т Л ( д „ dWn\

= - ро + 5кД ^gzKz-^). (2°)

Здесь 5 — вспомогательный параметр, равный единице.

Таким образом, схема определения компонент вектора скорости течения в момент времени tn+1 представляется в виде четырехэтапной схемы расщепления:

1) по известным в начальный (или предыдущий) момент времени полям скорости и плотности из уравнений (11) находятся U, W;

2) из решения задач (12), (17) - (20) определяются U*,W*;

3) из решения задачи (13), (15), (16) находится p1;

4) по формулам (14) вычисляются значения составляющих скорости Un+1,Wn+1. Расчеты температуры и солености воды выполняются в два этапа. На первом этапе

рассчитывается конвективный перенос из решения уравнений

T - Tn тт,п дТn ттт dTn ГТ" + u + W —— = 0,

Дt дх dz

£_£! + и-^ + w d-f = 0. (21)

Дt дх dz

z

На втором этапе решаются уравнения Тп+1_т

А г

(к (к

" хТ дх ) дг\ гТ дг )

дх

сп+1 _ с = дс п+Л + (к дсп+1 ^

А г дх \ хС дх ) дг \ хС дг )

с учетом граничных условий (5), (6).

По найденным значениям Тп+1, Сп+1 с помощью соотношения (8) или (9) определяется плотность воды рп+1.

Таким образом, исходная задача свелась к последовательному решению задач для уравнений гиперболического, параболического и эллиптического типа. В работе [10] отмечается, что в задачах с преимущественным направлением течения схемная вязкость в явной схеме с разностями против потока незначительна и характеристики течения рассчитываются достаточно точно. Поэтому для решения уравнений (11) и (18) применяется схема с разностями против потока [10]:

-п+1 _ -п я- яр

-3_3 + ип — + Шп—

Аг + °гз дх + ' гз дг

0.

(22)

Здесь

д- -п - п

дх АХг

д- - п - п - г+13 - г3

дх Ах г +1

д- - п -п - г3 - г3-1

дг Агз

д- -п - п - г3+1 - г3

дг Аг3+1

при Щ > 0, при Щ < 0, при > 0,

при < 0

(23)

-п = - (гп, хг, ), хг = хг-1 + А хг, = Х^-1 + А .

Ограничения на шаг по времени при решении уравнений (11), (21) по схеме (22), (23) получаются из условия [10]:

Аг <

\и |

Ах

шах + \Ш \

Аг

-1

Эллиптические и параболические уравнения решаются с помощью метода конечных элементов [9, 11]. Неизвестные функции представляются в виде разложения по первым элементам некоторого базиса. Коэффициенты при базисных функциях подлежат определению. После подстановки в соответствующие уравнения разложений вместо неизвестных функций определяются невязки, зависящие от неизвестных коэффициентов. Из равенства нулю интегралов по области, в которой ищется решение задачи, от произведения невязки на базисные функции для определения коэффициентов получается система линейных алгебраических уравнений.

Базисные функции выбираются следующим образом: с каждым узлом сетки связывается одна базисная билинейная функция, равная единице в этом узле и нулю на всех элементах, кроме тех, которые содержат данный узел. На каждом четырехугольном элементе

Рис. 2. Четырехугольный элемент сетки: Ах, Аг — размеры элемента, х € [0, Ах], г € [0, Аг] — локальные координаты.

не равными нулю являются четыре базисные функции, которые в системе координат, связанной с данным элементом, имеют следующий вид (рис. 2):

(Дж — х)(Дг — г) х(г — Дг) хг ¡ (Ах — х)г

ДхДг ' ДхДг ' ДхДг' ДхДг

Затем формула Грина

= —I дС!(П + I ¡С сов(и, г) (Г

дг I д г I

применяется к членам уравнений, содержащим вторые производные. Цель этого шага — уменьшить порядок производной в уравнении, что является стандартным приемом в методе конечных элементов.

Граничные условия учитываются следующим образом. Если на части границы поставлено условие Неймана, то в соответствующей части интеграла по границе в уравнениях производная функции по нормали заменяется на заданное значение. Если на части границы поставлено условие Дирихле, то соответствующая часть интеграла по границе в уравнениях не учитывается. При формировании системы уравнений в строках, соответствующих неизвестным, лежащим на этой границе, все элементы, кроме диагонального, зануляются, диагональный берется равным единице, а в правую часть подставляется значение известное из условий Дирихле.

В результате получаются системы линейных алгебраических уравнений относительно неизвестных коэффициентов. Матрицы этих систем линейных уравнений являются сильно разреженными и поэтому хранятся в компактной форме [12]. Хранится только следующая информация: ненулевые элементы матрицы, номера их столбцов, указатели на первые элементы каждой строки, для каждого элемента указатель на следующий элемент в строке. При таком виде хранения системы алгебраических уравнений можно решать итерационными методами, такими как простой итерации, Якоби, Зейделя, верхней релаксации, наискорейшего спуска, минимальных невязок [13, 14]. В работе использовался метод Зейделя.

Введение переменной р\ позволяет легко использовать построенный численный алгоритм для гидростатического приближения.

3. Примеры расчетов

Для проверки схемы расщепления и выведенных граничных условий используется следующий тест. Рассматривается течение неоднородной жидкости в замкнутой квадратной области [0 < х < 1 ,0 < г < 1], вызываемое заданным напряжением на водной поверхности; остальные границы области — твердые стенки.

Исходные уравнения (1), (2) записываются в переменных вихрь ш — функция тока —, которые связаны с исходными переменными соотношениями

дШ ди тт д— д—

ш = ^--ТГ, и = Ш = —

дх дг дг дх

Уравнения (1), (2) принимают вид

дш + ^дш + ^дш

д£ дх дг

д дш д дш 9 др!

= —К--1--К--1- ———

дг дг дх дх р0 дх

Д— = —ш.

(24)

Точное стационарное решение уравнений (24) в области [0,1] х [0,1] для Кг = Кх = 1

при

9 дР!

ро д х

1024х(1 — х) (1 — 2х)г(1 — г)3(1 — 3х + 3х2 — 3г + 9хг + 27х2 г) + +2048х3(1 — х)3(1 — 2х)(1 — г)(1 — 3г + 3г2 + 384г(1 — г )2+

+ 128(1 — 6х + 6х2)(3г — 2) с граничными условиями прилипания на твердых стенках

и = 0, Ш = 0

(25)

и условиями на верхней границе

ди

дг

—64х2(1 — х)2

Ш = 0

имеет вид

— = 16х2(1 — х)2г(1 — г)2.

(26)

Численное решение сформулированной задачи находится по описанному алгоритму методом установления с начальными условиями и = 0, Ш = 0. Расчеты показали, что наибольшее отличие численных данных от аналитического решения в окрестности верхней границы г = 0. Были проведены расчеты на сетках: Дх = Дг = 0.1; Дх = Дг = 0.02; на неравномерной сетке с шагами Дх = Дг = 0.1 внутри области со сгущением вблизи границ до Дх = Дг = 0.02. При выбранных параметрах разница между точным и численным решениями почти неразличима.

Численные эксперименты показали, что результаты расчетов для граничных условий (20) при 5 = 0 и 5 =1 практически не отличаются. Поэтому на водной поверхности можно использовать условие (20) при 5 = 0.

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

Картина ветровых течений зависит от силы ветра, стратификации и геометрии водоема. Для однородной жидкости в водоеме формируется одна циркуляционная зона, причем

с увеличением длины водоема возрастает максимальное значение скорости течения воды. Влияние размеров водоема, силы ветра и стратификации было изучено на следующих примерах. Рассматривалось пресноводное озеро длиной 500 м и максимальной глубиной 20 м с линейным перепадом температуры воды от 7 до 4 °С в верхнем восьмиметровом

-4.529 -2.901 -1.4ЭЭ 0.114 1.6С2 3.209 4.757 6,304 7Л52 9.Э99 10.947 ' ' ' 10.000

201.000

400.000

Рис. 3. Профиль горизонтальной компоненты скорости и (г) в вертикальном сечении х = 200 м для течения с одной циркуляционной зоной.

-2.453 1.774 1.096 -0.417 0.2С2 0.941 1.620 2.299 2Л77 3.656 4.335

1 1 1 1 Голой-

200.000

\

1400.900

. I

кшо.ооо

/

1800.900

яша.аоо

Рис. 4. Профиль горизонтальной компоненты скорости и (г) в вертикальном сечении х = 200 м для течения с двумя циркуляционными зонами.

слое и ниже до дна постоянным значением температуры 7 °C. При скорости ветра 9 м/с формируется одна циркуляционная зона (рис. 3). При малых скоростях ветра в водоеме образуются две циркуляционные зоны. На рис. 4 приведен профиль горизонтальной составляющей скорости течения воды для скорости ветра 5 м/с.

Для соленой воды при том же перепаде температуры и постоянной солености и для скорости ветра 20 м/с формируются две циркуляционные зоны. Для пресноводного водоема длиной 5000 м, максимальной глубиной 20 м и тем же профилем температуры при скорости ветра 5 м/с формируется одна циркуляционная зона с небольшими по размеру вихрями у дна. Для соленой воды такого же водоема формируются две циркуляционные зоны.

Разработанный программный комплекс будет применяться для исследования гидрофизических и химико-биологических параметров замкнутых водоемов.

Автор выражает благодарность В.В. Шайдурову, В.М. Белолипецкому и С.Ф. Пятаеву за постановку задачи, ценные советы и замечания.

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

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

[2] Белолипецкий В.М. Численное моделирование ветровых течений в стратифицированных водоемах // Водные ресурсы. 2001. Т. 28, № 2. С. 133-137.

[3] Математические модели циркуляции в океане / Г.И. Марчук и др. Новосибирск: Наука, 1980. 288 с.

[4] Кочергин В.П. Теория и методы расчета океанических течений. М.: Наука, 1978. 128 с.

[5] Судольский А.С. Динамические явления в водоемах. Л.: Гидрометеоиздат, 1991. 262 с.

[6] Генова С.Н., Лукавенко П.Н. Двумерная в вертикальной плоскости модель гидротермического режима непроточного водоема // Вычисл. технологии. 2002. Т. 7, № 4. С. 9-17.

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

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

[9] Флетчер К. Вычислительные методы в динамике жидкости. Т. 1, 2. М.: Мир, 1991.

[10] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.

[11] Флетчер К. Численные методы на основе метода Галёркина. М.: Мир, 1988. 352 с.

[12] Писсанецки С. Технология разреженных матриц. М.: Мир, 1988. 410 с.

[13] КАлиткин Н.Н. Численные методы. М.: Наука, 1987. 512 с.

[14] Самарский А.А. Введение в численные методы. М.: Наука, 1987. 288 с.

Поступила в редакцию 5 мая 2004 г., в переработанном виде — 2 ноября 2004 г.