CC BY
62
Вычислительные технологии
Том 4, № 6, 1999
РАСЧЕТ ДИНАМИКИ ТЕЧЕНИЯ МАЛЫХ ОЗЕР НА ПРИМЕРЕ ОЗЕРА ШИРА*
Л.В. ГАВРИЛОВА Красноярский государственный технический университет, Россия
Л. А. КОМПАНИЕЦ Институт вычислительного моделирования СО РАН
Красноярск, Россия e-mail: kla@cc.krascience.rssi.ru
For the numerical modeling of the wind flows in the lake Shira a two-stage procedure of obtaining velocities versus the depth is employed. At the first stage from equations of linear shallow water averaged velocities are found and the free surface elevation that will be further used for solving the diffusion equation. It is shown that for small lakes such as Shira the Coriolis forces can be neglected which enables us to simplify the numerical algorithm.
Введение
Озеро Шира в силу уникальности своих лечебных свойств имеет огромное значение для всего сибирского региона. Использование лечебных свойств оз. Шира началось в давние времена, но в последнее время сильно затруднено в связи с ухудшением экологической обстановки на курортном озере и в его окрестности. Возник ряд медико-биологических и экологических проблем курортного комплекса “Озеро Шира”, и одной из важных из них является проблема изучения качества воды в озере. При исследовании характеристик качества воды необходимо совместное изучение всей совокупности процессов — физических, химических, биологических и гидрофизических.
Основным фактором, связывающим абиотические и биотические процессы в водных экосистемах, является тепломассоперенос. Гидрофизические модели, описывающие основные физические процессы, обеспечивают химико-биологические модели информацией о температуре воды в характерных зонах и водо- и массообмене между зонами. При совместном исследовании гидрофизических и химико-биологических процессов существенно увеличивается сложность задачи и практический интерес представляют упрощенные гидрофизические модели.
* Работа выполнена при финансовой поддержке Федеральной целевой программы “Интеграция”, проект №73, 1998, Красноярского краевого фонда науки, грант 7F0036, гранта “Исследование поверхностных и внутренних гравитационных волн в жидкости” (в рамках Интеграционной программы фундаментальных исследований СО РАН) и Российского фонда фундаментальных исследований, грант №99-05-64695.
© Л. В. Гаврилова, Л. А. Компаниец, 1999.
1. Обоснование выбора математических моделей
При выборе математических моделей и методов расчета необходимо учитывать все особенности объекта, для которого проводятся расчеты. Так, для расчета конвекции глубоких озер типа Байкал используется квазидвумерная модель, в которой учитывается сжимаемость воды и два параметра Кориолиса [14].
Для морей и глубоких озер разработаны модели на основе уравнений гидродинамики турбулентного движения и распространения солености и температуры с учетом влияния сил Кориолиса, ветрового напряжения на поверхности и трения на дне [2, 3].
Для вытянутых в продольном направлении бассейнов (Кондопожская губа) используется двумерная в вертикальной плоскости x,z модель с двуслойной температурной стратификацией [7, 9].
Двумерные (в горизонтальной плоскости) так называемые уравнения мелкой воды применяются в случае, если горизонтальные размеры водоемов гораздо больше его глубины (Онежское озеро, Каспийское море, Чудское озеро). Получаемые при этом усредненные по глубине значения скоростей позволяют определить горизонтальную циркуляцию вод
[1, 4, 8, 12].
В то же время в ряде работ [3, 6, 7, 9, 11] рассматриваются методы, позволяющие находить скорости в зависимости от глубины, используя предварительно вычисленные величины усредненных по глубине скоростей и возвышения свободной поверхности.
Оз. Шира представляет собой бессточное озеро без островов размером 9x5 км, максимальной глубины 21 м (рис. 1). В озеро впадает одна речка Сон, но в силу малости притока все ее влияние сосредоточено в приустьевой зоне; основной тип движения — ветровые течения, и при математическом моделировании ветровых течений влиянием реки можно пренебречь. Для расчета пространственной картины течения использовалась методика [6, 7, 9]. Коэффициент турбулентной диффузии считался по формуле Прандтля — Обухова или принимался постоянным [5], хотя возможно его определение по двупараметрической модели [6, 10].
Рассматривается система уравнений
ди д'ц д/ ди
dt - lv = -ддХ + dZ [KzЖ
dv дп д( dv
д +lu = -gThj + a~z (KzШ
ди дv д'ш , .
дх + ду + дz ’
где х, у, z — пространственные переменные, при этом ось х направлена на восток, у — на север, z — вверх, U = (u,v,w) — вектор скорости течения, l = 2и sin p — параметр Кориолиса, и — угловая скорость вращения Земли, p — широта места, g — ускорение силы тяжести, п — отклонение поверхности жидкости от равновесного положения, Kz — коэффициент диффузии, в общем случае Kz = K(z). На боковой поверхности бассейна (вертикальные стенки) ставятся условия непротекания.
В зависимости от постановки граничных условий на свободной поверхности и дне получаем различные модификации для определения средних по глубине скоростей и возвышения свободной поверхности.
Рис. 1. Схема батиметрии оз. Шира.
Следуя [3, 9], для системы уравнений (1) рассмотрим граничные условия при г = п:
К 94 К Ра \Ш\
Кг— = КХ — Щ \ШХ, дг рь
ЪТ ду ЪТ Ра \Ш\
К ді = Кі Рь ^ к'
дп
дг =
а при г = — Н(х,у) условие проскальзывания в виде
Кгди/дг = К2IV|и, Кгду/дг = К2\У|У, ш = 0,
где ра — плотность атмосферы, рь — плотность воды, Щ(г) = (шх(х,у,1), шу(х,у,1)) — вектор скорости ветра, К, К2 — безразмерные коэффициенты касательного напряжения и придонного трения, V = (и, у) — вектор осредненных по глубине скоростей.
Интегрируя соотношения (1) от дна до поверхности с учетом граничных условий и предполагая малость п/Н, получим обычные уравнения мелкой воды с учетом ветрового воздействия и трения на дне [3, 7]:
ди , дп 1 ^ и^К ^ РаШх|Щ 1
И + 9вх: =1у — К2^ + К1^^~,
ду , дП , К I , К РаШу \Ш |
дг + % = —1и — К2~Н + К1^Н~,
ш + Ч(ИУ) = о, (2)
где V = (д/дх, д/ду). Начальные условия для систем уравнений (1), (2) брались нулевыми.
Ввиду отсутствия притоков на всех берегах озера (вертикальные стенки) ставятся условия непротекания.
2. Колебание жидкости в бассейне под действием постоянного ветра
Система уравнений (2) может быть решена независимо от системы (1) и позволяет определить поведение усредненных по глубине скоростей и возвышения свободной поверхности. Для решения этой системы уравнений методом конечных разностей на плоскости стандартным образом строится разностная сетка с шагом Ах = А у, при этом береговая линия заменяется ломаной линией со звеньями, параллельными осям координат. На береговой линии используется разностная аппроксимация условия непротекания Уп = 0, где Уп — проекция вектора У на нормаль к боковой поверхности. Численные алгоритмы для системы уравнений (2) известны, например, в [4] выписаны формулы неявной разностной схемы.
Нами на первом этапе исследования использовалась явная условно устойчивая двушаговая разностная схема второго порядка аппроксимации Мак-Кормака.
Для определения основных параметров движения в оз. Шира под воздействием ветра был проведен ряд расчетов в одномерном случае.
В тестовых расчетах были взяты 53 точки и шаг расчетной сетки равнялся 138 м, что соответствует длине оз. Шира. Глубина бассейна равнялась 10 м, что соответствует средней глубине озера. Из литературы известно, что наиболее часто применяемое значение ветрового коэффициента меняется от 1.2 • 10-3 до 2.6 • 10-3, и из одномерных тестовых расчетов получалось, что при коэффициенте К = 1.2 • 10-3 и скорости ветра 10 м/с имеем колебание поверхности с полупериодом около 700 с. Рис. 2 показывает возвышение свободной поверхности через 70, 150, 360 и 430 шагов по времени (расстояние в метрах). По оси х отложены расчетные точки. Расчеты проведены при А1 = 10 с, д = 10 м/с,
Рис. 2. Возвышение свободной поверхности в за- Рис. 3. Мареограмма в последней точке зада-даче о колебании жидкости в бассейне под дей- чи о колебании жидкости в бассейне под действием постоянного ветра (одномерный случай). ствием постоянного ветра (одномерный слу-
т
0.010 -і
0.03 -|
-0.005 -
0.005 -
0 -
-0.010
о
10
20
30
40
50
60
1000 2000 3000 4000 5000 6000
X
і
чай, ровное дно).
Рис. 3 показывает мареограмму в последней расчетной точке. По оси I отложены шаги по времени.
Уменьшение амплитуды колебаний может быть связано с выходом на стационарное решение; подобный эффект наблюдался при расчетах двумерных в вертикальной плоскости течений, индуцированных ветром, в Кондопожской губе [9].
Увеличение длины области расчетов в 2 и 4 раза приводит к почти линейному увеличению амплитуды и периода решения.
При одномерных расчетах в бассейне с вертикальными стенками, полого переходящими в наклонное дно в соответствии с батиметрией оз. Шира (рис. 4), получаем мареограмму, показанную на рис. 5; колебательный тип движения сохраняется, характер течения существенно зависит от аппроксимации граничных условий для скорости на промежуточном шаге в разностной схеме Мак-Кормака.
В двумерных расчетах для модельного прямоугольного бассейна использовалась расчетная сетка 55 х 19 с шагом 138 м по обоим направлениям, глубина бассейна полагалась равной 10 м. Для расчета применялась разностная схема Мак-Кормака. В случае, когда ветер дует под углом в 450 к оси х и скорость ветра составляет 10 м/с, картина течения изображена на рис. 6. Здесь также наблюдается колебательный характер течения.
Оз. Шира по своим размерам можно отнести к малым озерам, и, естественно, возникает вопрос о влиянии параметра Кориолиса на течение в численных расчетах. С этой целью были проведены расчеты для системы уравнений (2) с I = 0 и I = 0. При этом считалось, что параметр Кориолиса — постоянная величина, и = 7.292 • 10-5 рад/с, ^ = 540.
Расчеты проводились на равномерной разностной сетке 34 х 27, расстояние между расчетными точками 276 м, шаг по времени 10 с, д = 10 м/с, К = 1.2 • 10-3,
К2 = -2^-, где С = — Н1/6, коэффициент Маннинга кт = 0.02 м 1/3 с, W = (10, 0), С 2Н кт
число шагов по времени 200. Рис. 7 дает качественную картину скоростей при расчете по
0.03
О 1000 2000 3000 4000 5000 6000
£
Рис. 4. Батиметрия модельного бассейна в задаче Рис. 5. Мареограмма в последней точке зада-о колебании жидкости в бассейне с переменным чи о колебании жидкости в бассейне под дей-дном. ствием постоянного ветра (одномерный слу-
чай, переменное дно). 1 — аппроксимация граничного условия для скорости на промежуточном шаге в разностной схеме МакКормака иN = 0, 2 - то же, иN = — UN-1. N -последняя расчетная точка.
Рис. 6. Возвышение свободной поверхности в задаче о колебании жидкости в бассейне под действием постоянного ветра. а — 70 шагов по времени (двумерный случай), б—150 шагов по времени (двумерный случай).
Рис. 7. Плановая картина течения при расчете по разностной схеме Мак-Кормака без учета (а) и с учетом (б) сил Кориолиса.
схеме Мак-Кормака без учета и с учетом сил Кориолиса. Видно, что второй случай почти не меняет картину течения качественно, что наблюдается даже для более крупных озер (типа Чудского [6]), и количественно, что связано с небольшими размерами озера. Так, разница между величинами скоростей, полученными с учетом и без учета сил Кориолиса по формуле [6]
(Е"=, ((/!1) - /Г); + (/21) - л(2>)2)<)1/2 (£"=,((/Г )2 + (/21))2М1/2 '
где /,, /2 — компоненты скорости в направлении х и у, индекс 1 соответствует расчету с учетом параметра Кориолиса, составляет 5%.
3. Пространственная картина течения
Для определения изменения скоростей в зависимости от глубины использовалась система
(1). При решении первых двух уравнений системы (1) величины дцдх, дц/ду, V считаются известными из решения системы уравнений (2) на каждом временном слое, вторые производные по z аппроксимируются неявным образом, иначе появляется сильное ограничение на шаг по времени. Это приводит к необходимости применять для решения первых двух уравнений из (1) метод матричной прогонки. Но так как влияние параметра l несущественно, то полагая его равным 0, можно вместо матричной прогонки применять две скалярные. Для решения третьего уравнения применяется схема Эйлера.
Относительно выбора Kz в многочисленной литературе нет определенных рекомендаций [13]. Применение двупараметрической модели турбулентности сильно усложняет алгоритм, и при отсутствии точных данных о величине параметра Kz в уравнениях (1) возникает необходимость выбора его приближенного значения. Коэффициент Kz был принят следующим:
1. Kz = K(z) в соответствии с формулой Прандтля — Обухова
Kz = (0.05г)2\т |/(рь Ko )e"“z = Ke-az,
Ko = (0.05п)2|т\/(2Pbl), r = wy/Ko/(2l), a = у/ЩЩ, (3)
т — ветровое напряжение, т = K\paW\W\;
2. Kz = Kcp = const как усредненное значение по глубине для Kz, посчитанного по формуле Прандтля — Обухова
Kcp = K/(Ha)(1 - e"“H);
3. Kz = const = 10-2 м2/с.
■000995 0.0« -0.0216 0.08S8
Рис. 8. Распределение скоростей по глубине в прямоугольном бассейне.
На рис. 8 показана зависимость скорости и от глубины для точки (10,10) модельного двумерного прямоугольного бассейна постоянной глубины 10 м при различных Kz и коэффициентах ветрового напряжения (а — K = 1.2 • 10-3, б — K = 2.6 • 10-3).
Анализ численных расчетов показывает, что при Kz = 10-2 = const мы имеем на поверхности большее значение скорости, чем при Kcp и в соответствии с (3). При этом на дне ситуация обратная (на рис. 8 по вертикали отложены расчетные точки).
4. Расчеты для озера Шира
Рис. 9. Распределение скоростей и и V по глубине для двух точек оз. Шира. а — 70, б — 220 шагов по времени.
После методических расчетов, позволивших убедиться в работоспособности численных алгоритмов, были проведены расчеты для оз. Шира. Алгоритм оставался прежним: сначала рассчитывалось возвышение свободной поверхности по разностной схеме для уравнений линейной мелкой воды, затем находились скорости в зависимости от глубины.
Графики рис. 9 соответствуют двум точкам центрального сечения оз. Шира вдоль оси х, ближе к левому берегу и в центре, К = Кср, Кх = 2.6 • 10-3. На рис. 9, аизображены скорости через 12 мин, на рис. 9, б — через 37 мин после начала ветрового воздействия. В этом случае на поверхности озера в центральной части его скорость больше, чем у берегов, и скорость в направлении х при западном ветре 10 м/с достигает 6 см/с.
5. Выводы
Для оценки ветровых течений в оз. Шира в приближении однородной жидкости численно реализована упрощенная математическая модель, достоверность расчетов по которой проверена качественно в сравнении с расчетами для других водных бассейнов. Модель позволяет оценить некоторые особенности течения в озере и, следовательно, правильно выбрать модель и численный алгоритм, учитывающие неоднородность жидкости, что существенно при расчете течений в стратифицированных озерах, каким является оз. Шира.
Список литературы
[1] Андреев О. А., Воробьева Л. В. Ветровая и стоковая циркуляция Ладожского озера (численные эксперименты). В “Моделирование и экспериментальные исследования гидрологических процессов в озерах”. Наука, Ленингр. отд-ние, Л., 1986, 17-21.
[2] Архипов Б. В., Корявов П. П. Трехмерная нестационарная модель динамики вод и изменения солености в водоеме. Там же, 10-13.
[3] Астраханцев Г. П., Оганесян Л. А., Руховец Л. А. Трехмерная математическая модель гидротермики замкнутого водоема. Там же, 13-16.
[4] Баклановская В. Ф., Пальцев Б. В., Чечель И. И. О краевых задачах для системы уравнений Сен-Венана на плоскости. ЖВМиМФ, 19, №3, 1979, 708-725.
[5] Белолипецкий В. М., Генова С. Н., Туговиков В. Б., Шокин Ю. И. Численное моделирование задач гидроледотермики водотоков. Новосибирск, 1994.
[6] Волкова Г. Б., Квон В. И., Филатова Т.Н. Численное моделирование ветровых течений в Чудском озере. Водные ресурсы, №3, 1981, 91-99.
[7] Добровольская З. Н., Епихов Г. П., Корявов П.П., Моисеев Н. Н. Математические модели для расчета динамики и качества сложных водных систем. Там же, 32-51.
[8] Добровольская З.Н., Корявов П.П., Симонов А. И. Расчет течений в Онежском озере с учетом антропогенного воздействия. Там же, 100-104.
[9] Добровольская З. Н., Симонов А. И. Математическое моделирование течений в стратифицированном водоеме. В “Моделирование и экспериментальные исследования гидрологических процессов в озерах”. Наука, Ленингр. отд-ние, Л., 1986, 6-10.
[10] Игнатова Г. Ш., Квон В. И. Одномерная модель сезонного термоклина в озерах. Водные ресурсы, №6, 1979, 119-126.
[11] Кочергин В. П. Теория и методы расчета океанических течений. Наука, М., 1978.
[12] Перминов С. М., Чечель И. И. Исследование течений, горизонтальной циркуляции вод и водообмена между восточной и западной частями Северного Каспия. В “Гидрофизика Северного Каспия”. Наука, М., 1986, 6-10.
[13] Филатов Н. Н. Динамика озер. Гидрометеоиздат, Л., 1983.
[14] Цветова Е. А. Влияние сил Кориолиса на конвекцию в глубоком озере: вычислительный эксперимент. ПМТФ, 39, №4, 1981, 127-134.
Поступила в редакцию 3 декабря 1998 г., в переработанном виде 31 марта 1999 г.
