CC BY
41
Дармаев Тумэн Гомбоцыренович, кандидат физико-математических наук, доцент, заведующий лабораторией вычислительных и геоинформационных технологий Научно-образовательного инновационного центра системных исследований и автоматизации Института математики и информатики Бурятского госуниверситета, тел. (+73012) 221215, [email protected]
Khandarov Fedor Vladimirovich, research associate, scientific and educational innovation center of system research and automation, Institute of Mathematics and Computer Science, Buryat State University.
Darmaev Tumen Gombotsyrenovich, candidate of physical and mathematical sciences, associate professor, head of laboratory of calculation and geoinformational technologies, scientific and educational innovation center of system research and automation, Institute of Mathematics and Computer Science, Buryat State University.
УДК 004.942
© Т.В. Якубайлик
АДАПТАЦИЯ И ВЕРИФИКАЦИЯ ТРЕХМЕРНОГО ЧИСЛЕННОГО АЛГОРИТМА ДЛЯ РАСЧЕТА ТЕЧЕНИЙ В НЕГЛУБОКИХ ЗАМКНУТЫХ СТРАТИФИЦИРОВАННЫХ ВОДОЕМАХ
Работа выполнена при финансовой поддержке междисциплинарного интеграционного проекта №56,
грант РФФИ №13-05-00853
Проведена адаптация известной модели GETM [1] для расчета течений в неглубоких замкнутых стратифицированных водоемах. Данная программа была написана для моделирования течений в морях и заливах, но в силу универсальности позволяет рассчитывать течения практически в любых водоемах. Для этого требуется соответствующая настройка модели, так как процесс моделирования конкретного объекта подразумевает выбор из большого набора параметров. При помощи GETM проводились расчеты для Северного моря [2], Балтийского моря [3], озера Альпах (Швейцария) [4], водоемов-охладителей Шатурской ГРЭС [5] и многих других водоемов [6].
Ключевые слова: численный алгоритм, неглубокий водоем, стратифицированная жидкость.
T.V. Yakubaylik
ADAPTATION AND VERIFICATION OF THREE-DIMENSIONAL NUMERICAL ALGORITHM FOR CALCULATION OF FLOWS IN SHALLOW CLOSED STRATIFIED RESERVOIRS
The model GETM [1] has been adapted for calculation of flows in shallow closed stratified reservoirs. This program has been written for simulation of flows in seas and bays, but because of its versatility it allows to calculate the flow in various basins. It requires a corresponding adjustment of the model, as the simulation of specific object implies a choice of large set ofparameters. GETM was used for the calculations in the North Sea [2], in the Baltic Sea [3], in Lake Alpnach (Switzerland) [4], for cooling ponds of Shaturskaya power station [5] and in many other water bodies [6].
Keywords: numerical algorithm, shallow water body, stratified fluid.
Введение
Для моделирования течений в неглубоких соленых замкнутых стратифицированных озерах использовалась модель GETM. Сейчас это наиболее динамично развивающаяся модель, специализированная для расчета течений в различных водоемах (в отличие, например, от универсальной MITgsm, которая позволяет наряду с водоемами рассчитывать и движение в атмосфере и за счет этого имеет некоторые ограничения в вариативности). В GETM возможен выбор из четырех вариантов сетки для дискретизации по вертикали, достаточно большой выбор разностных аппроксимаций для расчета адвекции (по сравнению, например, с POM). Также она позволяет подключать много разных моделей вертикального турбулентного перемешивания (для расчета коэффициента вертикального турбулентного обмена).
1. Описание алгоритма
Система уравнений, используемая в этой программе, является модификацией хорошо известной и часто используемой для расчетов российскими и зарубежными учеными, системой уравнений гидро-
физики, основанной на приближениях Буссинеска и гидростатики. По сравнению с системой уравнений, выписанной, например, в [7], из уравнений не исключается вертикальная молекулярная вязкость.
Для построения численного алгоритма используется метод расщепления на моды [8], который позволяет рассчитать свободную поверхность с наименьшими затратами машинного времени, вычисляя скорости переноса отдельно от трехмерного расчета скорости и термодинамики. Уравнения для внешней моды скорости получаются интегрированием по глубине уравнений для внутренней моды. Вычисление по явным разностным схемам ограничено шагами дискретизации по пространству Ах и времени АА. Уравнения импульса, усредненного по глубине, вычисляются с маленьким шагом по времени во внешней моде и полные уравнения рассчитываются с более длинными шагами по времени при расчете внутренней моды. Для внешней моды решаются уравнения мелкой воды, которые описаны для модели в а - координатах в [8]. Модель, использующая расщепление на моды, способна охватить более широкий спектр пространственных и временных масштабов при одинаковой стоимости вычислительных ресурсов.
Для дискретизации по вертикали используются общие вертикальные координаты. Для пространственной дискретизации используется ступенчатая С-сетка Аракавы. Дискретизация проводится на основе метода конечных объемов. Члены горизонтальной и вертикальной диффузии аппроксимируются схемами второго порядка с центральными разностями. Пакет GETM дает большую возможность выбора схем для членов переноса. Для расчета адвекции во внешней моде запрограммировано девять конечно-разностных схем: от простой схемы первого порядка с разностями против потока до нелинейных монотонных схем третьего порядка аппроксимации. Те же самые схемы доступны для расчета горизонтального переноса во внутренней моде. Для расчета вертикальной адвекции во внутренней моде есть возможность выбора только из пяти схем (первого, второго и третьего порядка аппроксимации), большинство из которых являются монотонными.
Коэффициент вертикального турбулентного обмена рассчитывается с использованием моделей турбулентного замыкания разного уровня сложности. Существует возможность расчета с постоянным коэффициентом вертикального турбулентного обмена. Запрограммировано несколько вариантов модели. Также присутствует возможность расчета по моделям второго порядка. Коэффициенты горизонтальной турбулентной диффузии берутся постоянными. Также присутствует возможность выбора из нескольких различных вариантов уравнения состояния, включающих линейное уравнение состояния и уравнение состояния по формуле ЮНЕСКО.
2. Выбор сеток и разностных аппроксимаций членов переноса на основе тестовых расчетов
Для тестирования программы были проведены расчеты при постоянном коэффициенте вертикального турбулентного обмена в случае однородной жидкости, в процессе которых определились параметры сеток по горизонтали и вертикали. Расчеты проводились для бассейна прямоугольной формы с ровным дном глубиной 50 м и 25 м при западном ветре 5 м/с. На рис. 1 представлен годограф скорости в центре бассейна, (а) полученный для аналитического решения [9] и (б) полученный в результате расчета в программе GETM для бассейна глубиной 25 м. Расчеты проводились до выхода на стационар. Получено не только качественное, но и неплохое количественное совпадение при шаге по пространству 200 м как по оси Ох, так и по оси Оу. По вертикали оптимальное количество слоев оказалось равным шестидесяти.
0,004 - 1 \
0^02 - 0.004 ^ \ \ \ \ \
роо 0,005 0,010 0,015 -\ - 1
* и, м/с -0.005 \ 0 0.005 0.01 0.015 0.02
-0,00?- \ и, м/с
-0,004 - -0.004 — \ \
-0,006 - 15 -0.008
-0,008 - V, м/с ---
Рис. 1. Сравнение аналитического и численного решений для однородной жидкости: а) модель с учетом бокового обмена [9], б) численное решение
В целях экономии машинного времени для расчета переноса в уравнениях движения была взята наиболее простая схема первого порядка с разностями против потока. В озерах с устойчивой стратификацией в летний период профили температуры и солености несильно изменяются на небольших масштабах времени (менее полусуток). Так как стратификация в летний период достаточно сильная, профили имеют вид «ступеньки». Для адекватного описания данного физического процесса была выбрана схема, наилучшим образом сохраняющая профили со временем. Для расчета полей температуры и солености после ряда тестов была выбрана схема TVD-P2-PDM (третьего порядка, монотонная), наиболее адекватно отражающая картину распределения рассчитываемых величин по глубине. Более простые схемы «размазывали» профиль со временем. Более сложные схемы сильно увеличивали время счета.
3. Адаптация алгоритма на тестовых расчетах в бассейнах различной геометрии
Для адаптации численного алгоритма был проведен ряд тестовых расчетов для модельных водоемов простой формы: прямоугольного бассейна с ровным дном, цилиндрического бассейна с ровным дном, параболического бассейна и бассейна с береговой линией реального озера (с линейными размерами примерно 4 км на 9 км) и ровным дном. Известно из натурных наблюдений, что при постоянном ветре скоростью 5 м/с скорости течения воды на поверхности водоема имеют величину на два порядка меньше (около 5 см/с). Не должно быть сильного перемешивания, так как предполагается устойчивая стратификация. При постоянном ветре с наветренной стороны должен наблюдаться апвел-линг (поднятие на поверхность холодной воды). При переменном ветре ожидается появление внутренних волн. Все эти эффекты были определяющими при адаптации пакета GETM к расчету течений в неглубоком стратифицированном замкнутом водоеме.
На рис. 2 представлены профили восточной и северной компонент вектора скорости при расчете в прямоугольном бассейне длиной 10 км, шириной 7 км, глубиной 25 м, при постоянном ветре 6 м/с. В этом случае на поверхности величина вектора горизонтальной скорости составляет около 5,6 см/с. Сохранение стратификации подтверждается приведенным ниже рис. 3. Эффект апвеллинга проиллюстрирован на рис. 4 (при северо-западном и юго-западном ветре). Как и ожидалось, подъем холодной воды происходит с подветренной стороны.
zonal vel.[ttiK=3arlatc=13,lonc=25]
Рис. 2. Профили скорости в прямоугольном бассейне после шести часов расчета с северо-западным ветром величиной 6 м/с: а) восточная скорость, б) северная скорость
teniperature[time=30,latc=20,lonc=27] salinity[time=30,latc=20,lonc=27]
1
5 10 20 13 14 IS 17 18 1
Рис. 3. Профиль температуры (слева) и солености (справа) в середине прямоугольного бассейна после 6 часов расчета (квазистационарное течение)
На рис. 5 приведено распределение горизонтального вектора скорости на поверхности бассейна при северо-западном и юго-восточном ветре.
90.11 90 16 90.18 93.20 90.12 9124 30.26 90.28 17 85 Ц| до.16 90.18 90.20 90 22 90.24 90.26 90.2В
Рис. 4. Эффект апвеллинга: слева при северо-западном ветре, справа при юго-западном ветре
Рис. 5. Распределение горизонтального вектора скорости на поверхности бассейна: слева при северозападном, справа при юго-западном ветре 6 м/с
Рис. 6. Распределение горизонтального вектора скорости на поверхности бассейна при постоянном западном ветре 7 м/с, действующем на всей поверхности квадратного бассейна 100 на 100
км глубиной 100 м, по данным расчетов из [10]
Эти результаты сравнивались с расчетами из работы [10], выполненными для прямоугольного бассейна 100 на 100 км с ровным дном, глубиной 100 м, при постоянном западном ветре 7 м/с (рис. 6). Картина поля скоростей на поверхности бассейна в наших расчетах аналогична картине, полученной в расчетах других авторов.
Далее приведены результаты расчетов в прямоугольном бассейне с ровным дном при переменном ветре. На присутствие внутренних волн указывает и рис. 7, на котором представлено изменение изотерм в фиксированной точке бассейна с течением времени.
Рис. 7. Изотермы в фиксированной точке бассейна с течением времени (первый момент времени -начало штиля после переменного ветра): слева - в центральной и справа - в прибрежной
На рис. 8, 9 представлены изотермы на поверхности и в зоне термоклина (на глубине 6 м) при расчете в прямоугольном бассейне при переменном ветре. Ветер был запрограммирован по следующему сценарию: северо-западный, сначала 6 м/с в течение шести часов, затем 11 м/ с в течение часа или двух, затем полное отсутствие ветра (штиль) в течение двенадцати часов. Изотермы приведены, начиная с момента отключения ветра (начало штиля) с интервалом в два часа. Видно, что при отсутствии ветра температура на поверхности быстро выравнивается, в то время как в области термоклина продолжается движение теплых и холодных водных масс, что говорит о возможности образовании внутренних волн.
Для адаптации алгоритма к непрямоугольной береговой линии проведены расчеты в цилиндрическом бассейне с ровным дном и в бассейне с береговой линией реального озера и ровным дном.
■5-и7е1 йтрега1:иге[йте=0,г=б0]
■5-и7е1 1етрега1иге[1:1те=1О,г=б0]
•5.447е1 1втрега1иге[0те=20^=60]
□
гетрега1иге[ите-30,г-б0]
Рис. 8. Изотермы в прямоугольном бассейне глубиной 25 м на поверхности при переменном северо-западном ветре (изменение со временем с интервалом в 2 часа)
93.14 ».16 90.18 ».20 ».22 ».24 Я.26 90.28 ||Ьк1е СБбН
ЁтрегаШге[Ите=30,г=41]
□
□
Рис. 9. Изотермы в прямоугольном бассейне глубиной 25 м на глубине 6 м при переменном северо-западном
ветре (изменение со временем с интервалом в 2 часа)
На рис. 10 приведено распределение горизонтального вектора скорости на поверхности цилиндрического бассейна радиусом 2 км, глубиной 25 м при постоянном северном и юго-восточном ветре. Эффект апвеллинга проиллюстрирован на рис. 11 (при северном (слева) и юго-восточном (справа) ветре). Как и ожидалось, подъем холодной воды происходит с подветренной стороны.
0.01 0.02 0.03 0.04 0.05 0.06 0.07
0.02 0.03 0.04 0.05 0.06
Рис. 10. Распределение горизонтального вектора скорости на поверхности цилиндрического бассейна: слева при северном, справа при юго-восточном постоянном ветре 6 м/с
Рис. 11. Эффект апвеллинга в цилиндрическом бассейне: слева при северном ветре, справа при юго-
восточном ветре
На рис. 12 представлено распределение горизонтального вектора скорости в водоеме с реальной береговой линией и ровным дном, глубиной 25 м при постоянном ветре. На рис. 13, 14 представлены изотермы на поверхности и в зоне термоклина (на глубине 6 м) при расчете в бассейне с реальной береговой линией при переменном ветре.
90.14 90.16 90.18 90.20 90.22 90.24 90.26 90.2
90.14 90.16 90.18 90.20 90.22 90.24 90.26 90.2
Рис. 12. Распределение горизонтального вектора скорости на поверхности бассейна: слева при северозападном, справа при юго-восточном постоянном ветре 6 м/с
Рис. 13. Изотермы на поверхности при переменном северо-западном ветре (изменение со временем с интервалом в 2 часа)
Ветер был запрограммирован по сценарию, описанному выше. Изотермы приведены, начиная с момента отключения ветра (начало штиля) с интервалом в два часа. Так же, как и в случаях прямоугольного и цилиндрического бассейна, наблюдается эффект постепенного выравнивания температуры на поверхности при отсутствии ветра, в то время как в области термоклина продолжается циркуляционное движение, которое в случае произвольно изрезанной береговой линии имеет более сложный характер.
90.14 90.16 90.18 90.20 90.22 90.24 90.26 93.28 1опдШе (°ЕаЯ)
Рис. 14. Изотермы на глубине 6 м при переменном северо-западном ветре (изменение со временем с интервалом в 2 часа)
В случае достаточно сложных очертаний береговой линии алгоритм работоспособен. В бассейне с параболическим дном проводились тестовые расчеты для адаптации алгоритма к неровному дну. Бассейн на поверхности имеет диаметр четыре километра, максимальная глубина - 25 метров. Перемешивания слоев не наблюдается и в данном тесте (рис. 15). Видно влияние неровного
дна на изменение профилей температуры и солености со временем. Происходит более заметное их сглаживание со временем.
temperature[time=30,latc=20, lonc=27]
sal inity [time=30,latc=20, lonc=27]
temperature{time=30,latc=20,lorK=27] (°C)
ilimty[time=30,ratc=20,lonc=27] (PSU)
Рис. 15. Профили температуры (слева) и солености (справа) в параболическом бассейне после шести часов расчета с северо-западным ветром 6 м/с
На рис. 16 приведено распределение горизонтального вектора скорости на поверхности бассейна при северо-западном и юго-восточном ветре. Эффект влияния дна виден в изменении поведения вектора скорости на поверхности по сравнению с расчетом в цилиндрическом бассейне (рис. 10).
+5 içeal vel.[tlme=29,z=60], meridional vel.[time=29,z=60]
t5 îîs^I vei.[time=29,z=60], meridional vel.[time=29,z=60]
ч \ чч » Ч Ч \ 4 \ ч \ Ч 4
0.040 -0.035 -
ч N44 ^ > * * 4 \ \ v \ V *
0.097559832036»
0.01 0.02 0.03 0.04 0.05 0.06 0.07 longitude ("East) ^
Рис. 16. Распределение горизонтального вектора скорости на поверхности параболического бассейна: слева - при северо-западном, справа - при юго-восточном постоянном ветре величиной 6 м/с
Эффект апвеллинга для данного бассейна лучше виден на поперечном разрезе (рис. 17, 18).
Рис. 17. Изотермы в параболическом бассейне при северо-западном ветре на среднем разрезе по долготе:
слева - в начале расчета , справа - в конце расчета
Рис. 18. Изотермы в параболическом бассейне при юго-восточном ветре на среднем разрезе по долготе: слева - в начале расчета, справа - в конце расчета
На рис. 19 изображены изотермы при переменном ветре в фиксированных точках бассейна (после прекращения ветра, в течение 24 часов): в центральной - самой глубокой - и в более мелких в южном, восточном, северном и западном направлениях.
temperature[latc=20,lonc=25]
temperature[latc=12Jor>c=Z5]
0:00 11:00 11:00 17:00 2D:00 23:00 02:00
а) центральная глубокая точка
temperature[latc=20,1опс=33]
05:00 11:00 14:00 17:00 20:00 23:00 02:00
в) северная точка глубиной 11м
В:00 11:00 11:00 17:00 20:00 23 0D 02:00
б) южная точка глубиной 11м
В:00 11:00 14:00 17:00 20:00 23:0D 02:00
г) западная точка глубиной 11м
temperatu^tlatc=20Jonc=20]
OS:00 11:00 14:00 17:00 20:0D 23:00 02:00
д) восточная точка глубиной 14 м
Рис. 19. Изотермы в нескольких точках параболического бассейна после прекращения переменного северо-западного ветра
Заключение
В результате пакет GETM (GOTM) был адаптирован для расчета течений в неглубоких замкнутых водоемах с устойчивой стратификацией. Причем данный пакет позволяет проводить расчеты в водоемах с произвольной береговой линией и не только для ровного, но и для неровного дна.
Литература
1. GETM [Электронный ресурс]. _ Режим доступа: http://getm.eu/
2. Van Leeuwen S.M., van der Molen J., Ruardij P., Fernand L., Jickells T. Modelling the contribution of deep chlorophyll maxima to annual primary production in the North Sea // Biogeochemistry. 2012. Vol. 113. P. 137-152.
3. Gr'awe U., Burchard H. Storm surges in the Western Baltic Sea: the present and a possible future // Climate Dynamics. 2012. Vol. 39. P. 165-183.
4. Becherer J., Umlauf L. Boundary mixing in lakes. Part I. modeling the effect of shearinduced convection // Journal of Geophyscal Research: Oceans. 2011. Vol. 116.
5. Дебольская Е.И., Масликова О.А., Исаенков А.Ю. Использование программы GETM для математического моделирования термического режима Шатурских озер-охладителей // Динамика и термика рек, водохранилищ и прибрежной зоны морей: тр. VII конф. / РФФИ. ИВП РАН. 2009. С. 168-176.
6. Список публикаций, связанных с GETM [Электронный ресурс]. Режим доступа: http://getm.eu/index.php?option=com_include&Itemid=51
7. Марчук Г.И., Саркисян А.С. Математическое моделирование циркуляции океана. - М.: Наука, 1988.
8. Blumberg A.F., Mellor G.L. A description of a coastal ocean circulation model // Three dimensional ocean models, edited by N. S. Heaps. Washington: American Geophysical Union, D.C., 1987 P. 1-16.
9. Компаниец Л.А., Якубайлик Т.В., Гаврилова Л.В., Володько О.С. Аналитические решения для задач стационарного ветрового движения жидкости. - Красноярск: Изд-во СФУ, 2012. 112 c.
10. Добровольская З.Н., Епихов Г.П., Корявов П.П., Моисеев Н.Н. Математические модели для расчета динамики и качества сложных водных систем // Водные ресурсы. 1981. № 3. С. 33-51.
Якубайлик Татьяна Валерьевна, младший научный сотрудник ИВМ СО РАН, тел. (391-2) 498811, e-mail: [email protected]
Yakubaylik Tatyana Valerevna, junior researcher, ICM SB RAS, тел. (391-2) 498811, e-mail: [email protected]
