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

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

Том 10, № 5, 2005

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕМПЕРАТУРНО-СТРАТИФИЦИРОВАННЫХ ТЕЧЕНИЙ В СИСТЕМАХ ГЛУБОКИХ

ВОДОЕМОВ*

О. Ф. Васильев Институт водных и экологических проблем СО РАН Новосибирск, Россия e-mail: vasiliev@ad-sbras.nsc.ru

А. Ф. Воеводин, В. С. Никифоровская Институт гидродинамики СО РАН, Новосибирск, Россия e-mail: voevodin@hydro.nsc.ru

The numerical model is proposed for calculation of thermally-stratified unsteady flows in the systems of water bodies of oblong form. The basic equations governing this kind of flows are obtained by averaging the 3D equations of hydrodynamics. In addition to the initial and boundary conditions the junctions conditions are formulated in the places of stream confluence. Some numerical results are presented.

Введение

Математическая модель разработана для расчета гидротермических процессов в слабопроточных стратифицированных по плотности узких глубоких водоемах. Основные уравнения двумерной продольно-вертикальной модели получены путем поперечного осреднения трехмерных уравнений гидродинамики и предположения о гидростатичности давления. Полученная таким образом начально-краевая задача решается для системы двумерных областей со свободной границей. Кроме начальных и граничных условий во входных и выходных створах системы на свободной и донной поверхности для узловых вертикалей формулируются условия сопряжения, обеспечивающие баланс водных потоков и потоков тепла [1, 2]. Численный метод разработан на основе неявных разностных схем и метода дробных шагов. При этом для решения системы разностных уравнений предложен эффективный алгоритм, который является обобщением метода прогонки на случай областей, топологическая структура которых описывается графом "дерево".

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 03-05-65399) и Совета Программы поддержки ведущих научных школ (грант № 22.2003.5). © Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

1. Основные уравнения

Если поперечные размеры русла (или водоема) достаточно малы по сравнению с его протяженностью, то трехмерная гидродинамическая модель может быть аппроксимирована двумерной продольно-вертикальной моделью. В частности, такой тип модели может применяться к стратифицированным слабопроточным водоемам вытянутой формы. Основные уравнения двумерной продольно-вертикальной модели вертикального и продольного распределения плотности и скорости могут быть получены путем поперечного осреднения трехмерных динамических уравнений и уравнения неразрывности [1-4]. После осреднения основные уравнения в гидростатическом приближении выглядят следующим образом:

ди ди ди д . 1

~К7 + + = С + —

дЪ дх дг дх \ р0 J

д (Ьи) + д (Ьт)

,1 1 д (1 ди\ к + шг*аг) — Ьт;

где

к =

2

Е

г=1

дх

, дЬг

1+1 ТТ дх

+

дг

дк

дг

(1)

(2)

1/2

Л1 1 т = ^1и1и-

Здесь х,г — продольная (горизонтальная) и вертикальная декартовы координаты; Ъ — время; и(х, г, Ъ), /ш(х, г, Ъ) — горизонтальная и вертикальная составляющие скорости соответственно; Т(х, г, Ъ) — температура; Б(х, г, Ъ) — соленость; р(Т, Б) — плотность жидкости; р0 — характерное значение плотности; д — гравитационное ускорение; Ь = Ь1 + Ь2 — ширина русла; Ь1(х, г), Ь2(х, г) характеризуют положение и контур средней поверхности канала; ^(х,Ъ) — уровень свободной поверхности; V* — коэффициент турбулентной вязкости; т — сопротивление трения (касательное напряжение) на боковых поверхностях; Л — коэффициент трения.

Уравнения переноса тепла и солености получаются путем осреднения соответствующих трехмерных уравнений:

дТ идТ тдТ = 1 д / дТ4

дЪ дх дг Ьдг ( Т дг

дБ = 1 д (ЬК дБ

дЪ дх дг Ьдг \ 3 дг

(ькт«П + -±-(0, — »

V дг; роСо V Ь /

(3)

(4)

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

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

4 +1° = 0- О

дЪ дх

С

Ьи<г, ( — г0 = Н,

(5)

¿0

где В(х) = Ь(х, () — ширина русла на свободной поверхности; ° — расход; г0(х) — отметка дна; Н — глубина потока.

0

2

2

Используя уравнение (5), мы удовлетворяем кинематическому условию на свободной поверхности. Поэтому это уравнение часто называют кинематическим условием.

Для решения системы (1)-(5) необходимо определить коэффициент турбулентной вязкости vt и коэффициенты диффузии Kt , Ks . Самый простой способ — применение модифицированной теории турбулентных течений Прандтля для описания процессов перемешивания в стратифицированных потоках. Пример такого использования можно найти в статье А.А. Блумберга [1]. Однако в настоящее время для замыкания системы часто используются модели турбулентности на базе двух уравнений. Возможно, самой популярной из них является так называемая e — е (или k — е)-модель, включающая уравнения переноса энергии турбулентности e (или k) и скорости ее диссипации е. Также начинают широко применяться модели турбулентности со многими уравнениями, использующие турбулентные напряжения, а не турбулентную вязкость [3].

Для вычисления плотности используется уравнение состояния [5]: р = p(S,T).

1.1. Граничные условия

Кроме кинематического условия (5) на свободной поверхности задаются другие условия. 1. При z = Z

ди tw д( д(

vt^r = —, w = — + u—, dz po dt dx

К дТ Ф K dS 0 kt aZ = POCO, Ks dZ = °,

где tw — ветровое напряжение; Ф — поток тепла через единицу площади свободной поверхности.

При отсутствии ледяного покрова теплообмен через единицу площади свободной поверхности вычисляется как сумма [4-7]:

Ф = Фsn + Ф^ - Фbr - Фе - Фк - Фev + Фрг,

где Ф5П — поглощенная солнечная радиация; Ф^ — поглощенная длинноволновая атмосферная радиация; Фьг — длинноволновая радиация водной поверхности; Фе — затраты тепла на испарение или выделение тепла при конденсации; Фк — поток тепла, обусловленный конвекцией; — конвективный отток тепла с испаряющейся или приток с конденсирующейся водой; Фрг — конвективный приток тепла с атмосферными осадками.

Конвективный отток тепла с испаряющейся водой и приток с осадками и конденсирующимся паром определяются исходя из того, что температура осадков считается равной температуре воздуха, а температура испаряющейся воды — температуре моря. Если температура воздуха выше температуры замерзания воды, то атмосферные осадки считаются жидкими. В противном случае они считаются твердыми, т. е. снегом, и в соответствующей формуле учитывается скрытая теплота плавления снега, падающего в воду. При наличии ледяного покрова теплообмен вычисляется как сумма

Ф = -Ф„ + Ф^ - Ф

f,

где Фи — теплообмен между водой и льдом; — прошедшая через лед и поглощенная водой солнечная радиация; Ф^ — отток тепла с замерзающей или приток с тающей водой.

2. При z = z0

du dzo

V — = К |u|u, w =

^ dT ^ dS

= 0, = 0.

dz dz

На входе в водохранилище задаются расход воды Qin(t), температура воды T(t), соленость S (t) или плотность p(S, T) ив предположении, что глубина h здесь известна, может быть задано также вертикальное распределение продольной скорости uin(z,t). В случае эстуария такие условия могут быть заданы в сечении реки, расположенном достаточно далеко вверх по течению от устья. Аналогичные проблемы возникают и в выходном створе. С точки зрения кинематики более точным здесь было бы предположение о том, что на

т „ dzo

границе x = L вертикальная компонента скорости w линейно изменяется от w = u—— на

dx

дне до w = — на свободной поверхности. Это предположение в сочетании с уравнением (2)

позволяет проинтегрировать уравнение (1) в этом створе и определить uin. Таким образом, условия на границе при x = L будут следующими [3, 4]:

Z (L,t) = Zin(t),

w = dZin (z - zo) + u dzo (Z - z) dt h dx h '

O r-T I Sin(t),Tin(t), если u < 0,

S,T =i (6) не задается, если u > 0,

где Zin(t) — функция, определяющая колебания уровня в зависимости от времени; щ = u(L, z0, t).

При постановке условия на горизонтальную компоненту скорости вверх по потоку (x = 0) в некоторых типичных для стратифицированных течений случаях предпочтительно не задавать заранее профиль скорости. Поэтому, используя в этом случае тот же подход, что и при постановке условий на границе вниз по потоку, рассматриваются следующие условия в сечении вверх по потоку:

Q = Qin(t),

ST = J Sin(t),Tin(t) если u> 0, не задается, если u < 0,

и "мягкое" кинематическое условие

w = dZin (z - zo) + u d^o (Z - z) (7)

dt h dx h '

1.2. Условия сопряжения

При моделировании стратифицированных течений в естественных водоемах со сложной конфигурацией в плане целесообразно использовать их как систему рукавов. Схематически такие водоемы могут быть представлены (или составлены) в виде системы двумерных

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

Применяя основные законы сохранения в интегральной форме и делая некоторые предположения о динамике и кинематике потока, получим следующие условия, которые должны выполняться для узловой вертикали (или вертикали ветвления) [4, 8]:

1. Равенство уровней свободной поверхности на узловой вертикали для всех рукавов:

(i(t) = idem, i = 1, 2,...

2. Равенство давлений в рукавах в каждой точке узловой вертикали:

pi(z, t) = idem, i = 1, 2,...

3. Баланс водных потоков:

= 0, i = 1, 2,... ,

1 , если секция примыкает к узлу правой границей, 1 , если секция примыкает к узлу левой границей.

4. Баланс потоков тепла и солености:

у^ nibiTiUi = 0, i = 1, 2,..., "Y^nibiSiUi = 0, i = 1, 2,...

5. Будем предполагать, что вертикальное распределение вертикальной же компоненты скорости w линейно во всех ветвях на узловой вертикали (так называемое "мягкое" кинематическое условие):

' z - z0 \ dZ

к ) от

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

1. Горизонтальные скорости отсутствуют, а вертикальные определяются аналогично (7) при отсутствии уклона дна:

■т* ' '

к у ^

где х* — средняя отметка дна емкости; £* — уровень свободной поверхности емкости.

2. Уровни свободной поверхности в основном рукаве слева и справа от емкости совпадают с уровнем поверхности в емкости:

Сл С* Спр-

3. Интегральный баланс водных потоков:

^л - ^пр + П(С*) ^ = 0.

4. Баланс водных потоков в каждой точке узловой вертикали:

М. - (И„Р -I од ¡+(= 0.

5. Баланс потоков тепла и солености:

(ЪиТ). - (ЪиТ)Пр -(ЪиБ). - (ЪиБ)пр -

Ж аь

д_

дг

дТ*

^^Т" + 7Т" ( Пы*Т* - К*

дБ* д П—+ — ( Пы*Б* - К

„дБ*

дь ' ^ дг

где Т* = Т*(г,Ь), Б* = Б*(г,Ь) — температура воды и концентрация примеси в емкости

0

0

2. Численный метод

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

2.1. Построение сетки

Дискретная математическая модель одним из необходимых элементов включает построение разностной сетки. Исходная область для каждой секции при Ь > 0

0 < х < Ь^, 0 < г < ((,1)(х,Ь)

преобразуется в прямоугольную область, где г — номер секции,

0 < х < Ь({), 0 < п < 1,

с помощью замены независимых переменных:

«

г -

х = х, п = -тт •

с - 4г)

Такое преобразование позволяет построить сетку с одинаковым шагом по оси п и шагом Ах^ = Ь^/М различным для каждой области. В результате для каждой секции получим сеточную область

П(г) : I хп = пАх«; п = 0,...,М н \ Пт = тАп; т = 0,... М

Шаг по оси Ь является общим для всех секций и равняется АЬ. Считаем, что при Ь = 0 решение задано, т.е. известны для каждой секции

и(х, п, 0) = ит(х, п), ы(х, п, 0) = ыт(х, п), С(х, 0) = (т(х). Решение ищется на момент времени Ьк+х = + АЬ, к = 0,...

2.2. Расчет поля скоростей и поля плотности

Алгоритм решения задачи на каждом шаге по времени складывается из трех этапов.

1. Расчет свободной поверхности путем аппроксимации уравнений (5) и уравнения (1), проинтегрированного по глубине:

2. Расчет поля скоростей путем численного решения уравнений (1), (2).

3. Расчет поля температур и концентрации примесей путем численного решения уравнений (3), (4). Расчет поля плотности с помощью уравнения (6).

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

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

3.1. Пример 1

Рассматривалась система открытых русел, состоящая из трех секций (рукавов). При этом первая и третья секции (основной водоток) имели длину Ь\ = 500 м, Ь3 = 500 м и одинаковую ширину Б\ = Б3 = 100 м, дно считалось горизонтальным. Вторая секция (подводной канал) имела длину Ь2 = 200 м и ширину В2 = 10 м, дно так же было горизонтальным.

В качестве начального состояния при £ = 0 рассматривался установившийся режим, который был получен при граничных условиях: Q[)1) = 2 м3/с, Q02^> = 2.5 м3/с и QL = f (к). Кривая связи f (к) в конечном створе третьей секции задавалась так, что при к^3) = 2 м QL3) = 4.5 м3/с. Температура воды при £ = 0 для всех секций задавалась Т(х, 0) = 10°С. Теплообмен с окружающей средой не учитывался.

Нестационарный режим в системе моделировался следующим образом. В начальном створе второй секции за 5 мин температура воды поднималась до 25°С и поддерживалась постоянной в течение 20 мин и далее за 5 мин снова опускалась до 10°С. Таким образом имитировался аварийный сброс подогретой воды в основной водоток.

На рис. 1 показано распределение температуры по длине второй и третьей секций в различные моменты времени. В первой секции температура воды была постоянной до момента £ = 60, и при £ > 60 подогретая вода начинала проникать и в первую секцию.

3.2. Пример 2. Расчет термически-стратифицированного течения в непроточном водоеме

В непроточный водоем вытянутой формы (длина 500 км, уровень воды г = 33.0 м, общая площадь зеркала П = 22.75 тыс. км2, объем наполнения Ж = 141 км3), представленный

У ктйг

С

в виде водотока, состоящего из двух сообщающихся друг с другом участков в 250 км с параметрами, представленными в таблице.

В качестве исходных данных была использована информация, соответствующая предполагаемой морфометрии Аральского моря по источникам [13-15]. В течение длительного времени (£ = 1.5 года) в его мелководную часть поступает сток речной воды.

Гипотетический гидрограф (ежемесячный расход воды), поступающий в восточную часть водоема (восточный участок), изменение температуры воды, температуры атмосферного воздуха и солености воды во времени приведены на рис. 2. На рис. 3 представлено распределение температуры по длине участков в различные моменты времени. Изменение концентрации примеси было незначительным.

Рис. 1. Распределение температуры по длине второй (а) и третьей (б) секций (для £ = 20, 40, 60 мин).

Участок 1 Участок 2

Исходные данные восточный западный

мелководный глубоководный

Глубина, м со 4 ... 45

Ширина, км 15 ...110 17 ...45

Площадь зеркала, тыс.км2 16.374 6.376

Объем наполнения, км3 58 83

Температура, 0С 0 0

Соленость,

Рис. 2. Изменение расхода воды, температуры воды, температуры атмосферного воздуха и солености воды.

Г, °с

0,0 3,0 6,0 9,0 12,0 15,0 18,0 21 ,0 24,0 27,0 30,0

Восточный участок

Западный участок

Рис. 3. Распределение температуры для восточного и западного участков водоема для £ =120 и 400 суток.

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

[1] Blumberg A.F. Numerical model of estuarial circulation //J. Hydraulic Division, ASCE. 1977. Vol. 103, N 3. P. 295-310.

[2] Васильев О.Ф., Квон В.И., Чернышева Р.Т. Температурно-стратифицированное тече-

ние в водоеме вытянутой формы // Гидротехническое строительство. 1974. № 4. C. 35-43.

[3] Vasiliev O.F., Dumnov S.V. A two-dimensional mathematical model for salt water intrusion in an estuary // Proc. of the 20th IAMR Congress, Moscow, USSR. 1983. Vol. 2. P. 10-19.

[4] Vasiliev O.F. Vertical two-dimensional hydrodynamic models of water bodies: the state of the art and current issues // Proc. of the 5th Intern. Conf. on Hydro-Science and Engineering, Warsaw, Poland. 2002. Vol. 5. P. 3.

[5] Астраханцев Г.П., Руховец Л.А. Дискретная гидродинамическая модель климатической циркуляции глубокого озера // Вычисл. процессы и системы. 1986. Т. 4. С. 135-178.

[6] Wake A., Rumer R.R. Modelling ice regime of lake Eria // J. Hydraulic Division. ASCE. 1979. Vol. 105, N 7. P. 827-844.

[7] Квон Д.В., Квон В.Н., Семчуков А.Н. Численный расчет продольно-вертикальной термической структуры Телецкого озера в годовом цикле // Вычисл. технологии. 2000. Т. 5, № 3. C. 29-45.

[8] Васильев О.Ф., Воеводин А.Ф., Никифоровская В.С. Численное моделирование стратифицированных течений в системах открытых русел и водоемах разветвленной формы // Вычисл. технологии. 2004. Т. 9, № 2. C. 26-41.

[9] Марчук Г.И., Дымников В.П., Залесный В.Б. и др. Математическое моделирование общей циркуляции атмосферы и океана. Л.: Гидрометеоиздат, 1984.

[10] Самарский А.А., Вавишевич П.Н., Матус П.П. Разностные схемы с операторными множителями. Минск: Ин-т мат. моделирования РАН, Ин-т математики НАНБ, 1988.

[11] Васильев О.Ф., Атавин А.А., Воеводин А.Ф. Методы расчета неустановившихся течений в системах открытых русел и каналов // Численные методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1975. Т. 6, № 4. C. 21-30.

[12] Воеводин А.Ф., Шугрин С.М. Численные методы решения одномерных систем. Новосибирск: Наука, 1981.

[13] Иванова Л.В. Гидрологические проблемы Аральского моря // Водные ресурсы. 1992. № 2. C. 39-49.

[14] Гидрометеорология и гидрохимия морей СССР. Т. VII: Аральское море. Л.: Гидрометеоиз-дат, 1990.

[15] Гидрологический ежегодник 1978. Ч. 1 и 2, Т. II, вып. 3 (Аральское море, устье р. Аму-Дарьи). Алма-Ата, 1979.

Поступила в редакцию 24 мая 2005 г.