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

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

Том 7, № 4, 2002

ДВУМЕРНАЯ В ВЕРТИКАЛЬНОЙ ПЛОСКОСТИ МОДЕЛЬ ГИДРОТЕРМИЧЕСКОГО РЕЖИМА НЕПРОТОЧНОГО ВОДОЕМА*

С. Н. ГЕНОВА Институт вычислительного моделирования СО РАН

Красноярск, Россия e-mail: sv@icm.krasn.ru П. Н. ЛУКАВЕНКО Красноярский государственный университет, Россия

Two-dimensional in vertical plane wind flow model in non-flowing reservoirs is considered. Numerical algorithm is suggested in Boussinesq approximation and approximation of boundary layer for variable coefficient of vertical turbulence exchange. New variables: stream-function and vorticity are used. The examples of calculations of water temperature and wind flows in stratified reservoir are given.

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

В настоящей работе рассматривается численный алгоритм для двумерной негидростатической модели с переменным коэффициентом вертикального турбулентного обмена. Приводятся примеры расчетов ветровых течений в водоемах произвольной геометрии.

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

1.1. Двумерная в вертикальной плоскости модель

Основные допущения: приближение Буссинеска, приближение твердой крышки. Опускаются члены горизонтального турбулентного обмена. С учетом перечисленных допущений основные уравнения записываются в виде

dU TTdU TxrdU 1 dp д BU

ж + Udx + =- po dx +

dW TTdW BW 1 dp д dW p

"Т5Г + U-JT + W— =--/ + — K— + ^ g, (1)

dt dx dz p0 dz dz dz p0

dU dW = 0

dx dz

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проект №99-05-64695, Министерства образования РФ и СИБЕ, грант №ИЕС-002. © С.Н. Генова, П.Н. Лукавенко, 2002.

Уравнение переноса и диффузии для температуры Т:

дТ идТ щдТ = д.К —

дЬ дх дг дг т дг'

(2)

Здесь Ь — время; х, г — прямоугольные координаты, ось г направлена вертикально вниз; и,Ш — компоненты вектора скорости; р — давление; р = р(Т) — плотность воды; р0 — характерная плотность воды; д — ускорение свободного падения; К, Кт — коэффициенты вертикального турбулентного обмена.

Уравнения (1), (2) дополняются начальными и граничными условиями. В начальный момент задаются распределения всех искомых переменных. На боковых границах ставятся условия непротекания и отсутствия потоков тепла и примесей.

Граничные условия на свободной поверхности (г = 0)

W = 0, К

дг

Кт

Ро

дТ

дг

Би

еро'

(3)

Граничные условия на дне (г = Н)

и = 0, W = 0, Кт ^Т = — •

дг еро

(4)

Здесь т — напряжение ветра; Би — тепловой поток через свободную поверхность; с — удельная теплоемкость воды; Бы — поток тепла через донную поверхность. Теплообмен с атмосферой зависит от метеоданных и определяется по формуле

Би — Бд + Бси + Бо

где Бд — радиационный теплообмен; — затраты тепла на испарение; Бс — конвективный теплообмен. Для определения составляющих теплообмена используются известные формулы [6].

В переменных функция тока ф и вихрь и система уравнений (1) записывается в виде

дш ттдш ТТГдш д [ ди\ д др д (дК дW дК ди —+ и--+ W— = — К— + -—- +-------

дЬ дх дг дг \ дг ) р0 дх дг \ дх дг дг дг

(5)

и

д2ф д2ф

дх2 дг2 '

дф дг

W = -

дф дх

(6) (7)

Начальные условия: ф(х,г, 0) = ф0(х,г), ш(х,г, 0) = ш0(х,г). Граничные условия для системы (5), (6) следуют из их аналогов для системы (1):

ф = 0,

и

КоРо

при г = 0,

ф = 0,

дф дг

д 2ф

Ш = - дг2

На боковых границах для непроточного водоема ф = 0.

(8)

при г = Н.

9

Для оценки коэффициента вертикального турбулентного обмена К(ж,г,£), учитывающего сдвиговые эффекты и стратификацию, применяется формула Прандтля — Обухова [6]:

Г (0.05П)2^В при В > 0, К = \ 1 } (10)

[ Кт;п при В < 0 или К < Кт;п,

^ [ди\2 д др , .

где В = ----——; п — глубина квазиоднородного слоя, находящаяся по первой от

\дг ) ро дг

поверхности расчетной точке, в которой выполняется условие

(0.05zfc)УB\z=zk < Kmin, (11)

Kmin — фоновое значение коэффициента вертикального турбулентного обмена. Уравнение состояния воды берется в виде

р =1 - 1.5 ■ 10-4T.

Напряжение ветра находится по формуле

т/р0 = 1.25 ■ 10-2V2 [см2/с2], (12)

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

Для медленных течений в уравнениях движения (1), (5) можно пренебречь нелинейными инерционными членами.

1.2. Одномерная модель температурного режима водоема

Рассмотрим простейший случай формирования вертикальной структуры температурного режима непроточного водоема в одномерном приближении. Распределение температуры воды находится из решения задачи (2) - (4) при U = 0 и W = 0.

На формирование турбулентной структуры водоемов в основном влияют обрушение поверхностных и внутренних волн, диффузия энергии турбулентности, конвективные процессы, сдвиговые эффекты. Для оценки коэффициента вертикального турбулентного обмена применяется формула (10), в которой

B =( дМ2 + С м2 _ Jgdp

\0z) \dz) р0 dz

(u, v — горизонтальные составляющие скорости ветровых течений).

Для оценки коэффициента K определим вертикальные градиенты скорости из приближенного решения Экмана для ветровых течений:

g az

uek = --— [(тх + ту) cos(az) + (ту _ тх) sin(az)]

2pgazKo (13)

2poaKo

vek = 777"[(ту _ тх) cos(az) + (тх + ту) sin(az)] ,

где тх, ту — составляющие напряжения трения ветра; К0 — характерное значение коэффициента вертикального турбулентного обмена; а = у/ТКо, / — параметр Кориолиса

(f = 2w sin w — угловая скорость вращения Земли, ^ — широта местности). Из (13) следует

Ш+Ш ==ш*- т=(14)

Глубину верхнего квазиоднородного слоя можно оценить из следующего условия: толщина слоя, в котором К0 ^ ) + (~д~) ^ уменьшается в еп раз, равна толщине верхнего квазиоднородного слоя h1. Тогда из (14) получается hi = К0/(2f). Предполагая, что Ко = K(0), находим

К = (0-05п)2т . (15)

jW 4f2 + (0.05n)4gp* |z=o

(T Л 2 g др

- e-2az--—.

poKoJ po dz

В описанной методике интенсивность вертикального турбулентного обмена определяется действием ветра.

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

Рассматривается разностная сетка, равномерная по глубине и неравномерная по длине: Xi+i = Xi + Axi, zj+i = Zj + Az (i = 1, ...,ii, j = 1, ..., jji) (рис. 1).

Первые производные функций по x аппроксимируются следующим образом:

/ dF(x)\ = 1 / Axi-i () _ Axi \ Axi-i - Axi

V dx ).. = (Axi + Axi-iA Axi F (x)i+1'j Axi-1F (x)i-j AxiAxi-i F (x)i,j'

Численный алгоритм основан на методе расщепления по физическим процессам. Поле скорости находится из решений уравнений (5) - (7) по известным значениям плотности. Затем определяются температура и плотность. Предполагается, что Кт = К.

На первом этапе находится коэффициент вертикального турбулентного обмена К по формуле (10). Глубина квазиоднородного слоя h находится из условия (11). Далее рассчитывается вихрь скорости w. Уравнение (5) описывает два различных физических процесса — конвективный и диффузионный. Для построения численного алгоритма используется

Рис. 1. Разностная сетка для водоема произвольной формы.

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

Конвективный процесс описывается следующей задачей:

дЬ

Начальное приближение

дШ1 + и^ + = 0,

дх

дг

Ш1

t—

(х, г).

:1б)

и0 = 0, Ш0 = 0, ш0 = 0.

о _

0

Вертикальный турбулентный обмен с учетом напряжения трения ветра на свободной поверхности в каждом створе х = х^ описывается задачей

д^2 ~Ж

(к—

дг \ дг

д др д (дК дШ дК ди

т

г—0 Коро

Ш2

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

д2ф = — ТТ^) ш2

г—щ дг2 t—1„

:17)

где Щ = Н (х^) — глубина водоема в створе х = х^.

Численный алгоритм решения уравнения для ш, основанный на методе расщепления, представляется в следующем виде. Рассматривается элементарный временной интервал Ьп < Ь < Ьп+1, на котором вначале решается задача (16) с помощью схемы с разностями

против потока

где

. ,2+1 = п т гп+1/2 ,* ШП+1/2, >** ш1ад = ш2 1,3 — и 1,3 ш1 г,з - Wi,j ш1 г,з ,

:18)

ш

1 i,j

АЬ АЬ

Ш1 ¿,З - ш

п

1 i 1,3

) при ип+1 > 0

Ах К+и - шп3 при ип+1 < 0

:19)

ш

1 ¿,З

АЬ л п п А— 1ш1 i,j - ш1 ¿,3-1

АЬ Л ,п ,п

А" 1ш1 ¿,З+1 - ш1 ¿,З

при шп+1/2 > 0, при ш+1/2 < 0.

(20)

Здесь Щ1- = F(xi, г3- ,Ьп), ^п+1/2 = + п+1). Схема с разностями против потока яв-

и АЬ

ная и устойчивая при условии с =

Ах

< -.

Затем на этом же временном интервале в каждом створе х = x¿ решается задача (17) с помощью неявной схемы:

шп+1 - шп = д (Кдшп+1\ д дрп ^/дК^ дК^ дип

АЬ дг \ дг ) р0 дх дг\ дх дг дг дг

Для уравнения (21) применяется разностная схема:

/ ,п+1 , , Ш0,- — Ш.

2 ¿,3

= К,3+1( Ш2г+3+1 — Шп++1) — КЦ-1( Шп++1 — Ш2+3-1)

АЬ = (Аг)2

+ /*п 1 ^ ¿,З

*

**

(Аг)2

+ £

п

где

г = д ^КТ^- дШ, - дК, дЩ + ддр

п г,3

дг \ дх дг дг дг ) р0 дх Полученная система разностных уравнений решается методом прогонки.

На третьем этапе по найденному вихрю находятся функция тока ф и составляющие скорости и, Ш:

д 2фп+1 д2фп+!

дх2

+

дг2

-и

п+1

ип+1 = Шп+1 = - дфп+!

дг

Разностный аналог уравнения (22) имеет вид

дх

(22)

(23)

Ф

2-

г+1,3

,■ - Фгз -

Ахг

АХ~Г

(Фг,3 - Фг-1,3) Ф. - 2Ф + Ф 1

АХг(АХг + АХг-1)

+

Аг2

= -и

г,3 ,

(24)

Фг,1 = Фг,33 = 0.

Ф1,3 = Фгг,3

Эта система решается методом последовательной верхней релаксации:

Фгк+1 = (1 - 0)Ф& + Ш 1г ^Ф^ + Ш2гФк+1,з + Фк,+-1 + Аг2 и,

Ш 1г

2 1 +

Аг2

АхгАхг-1

Ш 2 г

2Аг2

Ахг(Ахг + Ахг-1)

Ш 3г

2Аг2

Ахг-1 (Ахг + Ахг-1)

в — параметр релаксации (1 < в < 2). Составляющие вектора скорости (и,Ш) определяются по (23).

На четвертом этапе по найденным составляющим скорости находится распределение температуры. Уравнение (2) описывает два различных физических процесса: перенос температуры течением и процесс, связанный с диффузией. Для построения численного алгоритма используется метод расщепления по физическим процессам с применением схемы с разностями против потока и неявной схемы для диффузионных процессов.

Перенос тепла течением описывается следующим образом:

^ + и1^ + Ш1^ = 0, г,

t—

ТП(х,г).

25)

дг дх дг

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

К Т2

дг

¿—о

СрРо'

дт1 = дг =

к Г2

дг

д_

дг

¿—н

К

дТ-2

дг СрРо'

(26)

Т2

t—^п

ТП+1(хг,г)

где Нг = Н(хг) — глубина водоема в створе х = хг.

Задачи (25), (26) решаются с помощью описанных выше схемы с разностями против потока и неявной схемы для диффузионных процессов.

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

Были проведены модельные расчеты по определению термоклина в оз. Шира по одномерной модели (в глубоководной части озера). Использовались реальные метеоданные для этой местности. Результаты расчетов согласуются с натурными данными (рис. 2). На рис. 2 также приведены результаты расчетов распределения температуры по двумерной в вертикальной плоскости модели, которые приближаются к натурным данным. Результаты расчетов температуры воды в глубоководной области с учетом и без учета конвективных членов в уравнении (5) практически совпадают. Из приведенных данных следует, что для оценки температурного режима в глубоководной области водоема можно применять одномерную модель.

Выполнена серия расчетов для водоема прямоугольной формы и водоема произвольной геометрии, схематизирующих оз. Шира в поперечном сечении. Для исследования влияния конвективных членов в уравнениях (1), (5) были выполнены расчеты с учетом и без учета конвективных членов.

Линии тока при скорости ветра 10 м/с приведены на рис. 3, а — 3, в, при скорости ветра 5 м/с — на рис. 3, г; ДТ = 16 °С — рис. 3, а, 3, б и 3, г, ДТ = 26 °С — рис. 3, в (ДТ = Ттах — Тт[п). Воздействие стратификации на картину течения демонстрируется на рис. 3, а и 3, в. Существенное влияние на картину течения в водоемах оказывает скорость ветра (рис. 3, а, и 3, г). На формирование течений в водоемах влияет также и профиль дна водоема, что видно из сравнения результатов, показанных на рис. 3, а и 3, б. Для сложного профиля дна (рис. 3, б) в водоеме образуются придонные циркуляционные зоны, что не наблюдается в водоеме прямоугольной формы (рис. 3, а).

На рис. 4 приведены профили горизонтальной составляющей скорости в среднем створе. Показано сравнение результатов расчетов в непрямоугольной области с учетом (сплошная линия) и без учета (пунктирная линия) конвективных членов. Как видно, отличие результатов расчетов по двум моделям незначительное.

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

Рис. 2. Распределение температуры по глубине: натурные и расчетные профили; Т\нат — начальное распределение температуры воды, Т2нат — распределение температуры воды через 2 мес., Т\ расч и Т2расч — расчетные профили, полученные по одномерной и двумерной моделям.

Рис. 3. Линии тока в водоеме прямоугольной (а, в, г) и произвольной (б) формы.

Рис. 4. Скорость течения воды в среднем створе водоема непрямоугольной формы при различных значениях скорости ветра: а — V = 5 м/с, б — V = 9 м/с, в — V = 12 м/с с учетом (сплошная линия) и без учета (пунктирная линия) конвективных членов.

Авторы выражают благодарность профессору В. М. Белолипецкому за постановку задачи и внимание к работе.

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

[1] Архипов Б. В., Солбаков В. В., Шапочкин Д. А. Двумерная вертикальная модель температурного режима водоема-охладителя // Вод. ресурсы. 1995. Т. 22, № 6. С. 653-666.

[2] Васильев О. Ф., Квон В. И. О теоретическом описании гидротермических явлений в водоемах-охладителях // Проблемы теплофизики и физической гидродинамики. Новосибирск: Наука, 1974. С. 100-111.

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

[4] Бочаров О. Б., Овчинникова Т. Э. Влияние стратификации на гидротермический режим водоема и требования к вычислительному алгоритму // Вычисл. технологии: Сб. науч. тр. / Сиб. отд-ние РАН; ИВТ. 1993. Т. 2, № 6. С. 155-162

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

[6] Белолипецкий В. М., Генова С. Н., Туговиков В. Б., Шокин Ю. И. Численное моделирование задач гидроледотермики водотоков. Новосибирск: Ин-т вычисл. технологий СО РАН, Вычисл. центр в г. Красноярске, 1994. 138 с.

[7] Belolipetskii V. M., Genüva S. N. Investigation of hydrothermal and ice refimes in hydropoewr Station bays // IJFCD. 1998. Vol. 10. P. 151-158.

Поступила в редакцию 4 декабря 2001 г., в переработанном виде — 11 марта 2002 г.