CC BY
46
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2011 Математика и механика № 1(13)
УДК 519.6:532.516
Б.О. Цыденов, А.В. Старченко ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЭФФЕКТА ТЕРМОБАРА В ОЗЕРЕ БАЙКАЛ В ПЕРИОД ВЕСЕННЕ-ЛЕТНЕГО ПРОГРЕВАНИЯ
Разработана математическая модель и построен вычислительный алгоритм, позволяющий воспроизвести процесс формирования и дальнейшего развития термобара в озере Байкал с учетом основных физических факторов. Полученные численные результаты демонстрируют возникновение циркуляционного течения вблизи берега и его смещение к центру озера с течением времени.
Ключевые слова: термический бар, температура максимальной плотности, приближение Буссинеска, конвекция, численный эксперимент.
Проблема «чистой воды» - одна из важнейших проблем, которые ставит перед человечеством научно-технический прогресс. По оценкам ученых через несколько десятилетий чистая пресная вода станет важнейшим ресурсом, поскольку она незаменима в отличие от других природных богатств Земли. К истощению водных ресурсов ведет не рост расходуемой воды, а её загрязнение. Озеро Байкал является самым крупным хранилищем пресной воды на планете (около 20 % мировых запасов).
Под термобаром понимается узкая зона в глубоком озере умеренных широт, в которой происходит погружение имеющей наибольшую плотность воды от поверхности до дна. С физической точки зрения, причиной формирования термобара является так называемый эффект уплотнения при смешении вод, т. е. аномальное изменение плотности воды. Известно, что плотность воды определенного солевого состава, находящейся при фиксированном давлении, достигает максимума при некоторой температуре - температуре максимальной плотности (ТМП). Поэтому если смешиваются две водные массы, имеющие общую боковую границу, и температура одной выше ТМП, а другой - ниже ТМП, то в результате получается смесь, которая будет тяжелее как первого, так и второго объемов. Естественно, более тяжелая вода должна опускаться, вследствие чего в месте смешения образуется как бы барьер для горизонтального перемещения воды. То есть термобар препятствует обмену водных масс между прибрежными и центральным районами озера, являясь в то же время зоной конвергенции этих масс (гидрологическим фронтом).
Для сохранения уникальности Байкала и его экосистемы необходимо понимание всех физических механизмов, участвующих в процессах водообмена и формирования качества его вод. С одной стороны, важность изучения термобара как явления, которое может оказать существенное влияние на процессы распространения загрязнения, состоит в том, что интенсивные нисходящие течения, возникающие между двумя конвективными ячейками, могут привести к быстрому распространению загрязнения из поверхностных слоев до очень больших глубин. С другой стороны, установлено, что придонные воды Байкала «моложе» и богаче кислородом, чем воды основного глубинного ядра. Все это вызывает интерес к исследованию термобара.
Явление термического бара впервые было обнаружено на Женевском озере швейцарским физиком и географом Ф. Форелем, который положил начало науке об изучении озер, образовавшей отдел океанографии — лимнологии. Честь вторичного открытия термического бара и, что самое главное, широкая разработка этой проблемы и освещение ее огромного значения при объяснении многих процессов в жизни озер принадлежит А.И. Тихомирову. Тщательные исследования, проведенные на Ладоге в 1957 - 1962 гг., позволили всесторонне оценить значение термобара для различных сторон жизни водоема. Заинтересовавшись статьей Тихомирова, переведенной на английский язык, американский ученый Роджерс вскоре исследовал и описал термический бар на озере Онтарио. Первая попытка исследования общих закономерностей водообмена в Байкале с помощью моделей была осуществлена в Государственном гидрологическом институте в 1966 - 1967 гг. Однако ограничение модели только южной частью озера не позволило получить достаточно общих результатов. Ввиду сложности создания таких моделей, невозможности моделирования на них эффектов, связанных с влиянием вращения Земли и вертикальной плотностной стратификацией, в дальнейшем этот вид исследования полей течений не получил развития [1, с. 95-99]. В 1995 г. С.Дж. Уолкер и др. [2] предложили трехмерную численную модель для проверки гипотезы о том, что глубинный водообмен инициируется плотностной неоднородностью, вызванной штормами. Однако в своих численных экспериментах они использовали параметры, не характерные для Байкала, поэтому осталось не совсем ясным, что же происходит в реальных условиях в озере. П. Д. Киллворт и др. [3], используя двумерную численную модель гидротермических процессов в водоеме в гидростатическом приближении, пытались проверить гипотезу о переменных ветровых воздействиях, способных вызвать глубинный водообмен. Но, как показали данные наблюдений, ни самые сильные шторма (ноябрь 1994 г.), ни осенние ветры не приводят к перемешиванию слоев глубже 80 - 100 м. Большинство современных исследователей связывает вентиляцию глубинных вод с эпизодически возникающими опусканиями дискретных объемов воды верхнего слоя, происходящими из-за так называемой термобарической неустойчивости столба жидкости. Механизм такой неустойчивости предполагает наличие внешней силы, которая может преодолеть потенциальный барьер, создаваемый архимедовыми силами. В работах Е.А. Цветовой [4 - 7] используются как двумерные, так и трехмерные модели с учетом сжимаемости воды, сил Кориолиса, адиабатического градиента температуры и других факторов. Расчеты выполнялись на прямоугольных и трапециевидных областях. Однако в этих работах нет анализа влияния отдельных физических факторов на процессы перемешивания вод озера. В работе [8] предложена математическая модель термобара в глубоком озере, согласно которой глубинная конвекция в окрестности фронта термобара может быть обусловлена более высокой минерализацией прибрежных вод озера по сравнению с минерализацией вод его основной части. С использованием двумерной гидродинамической модели в работе [9] численно исследуются плотностные течения, возникающие под влиянием минерализованного притока в области больших уклонов дна. Для оценки влияния отдельных факторов, таких, как температурная стратификация, различие в минерализации вод озера и притока, расход притока, выполнена серия численных экспериментов с целью выявления возможности проникновения вод притока в придонные области озера.
Целью данной работы является разработка математической модели и вычислительного алгоритма, основанного на численном методе высокого порядка точно-
сти, а также воспроизведение динамической картины формирования и дальнейшего развития явления термического бара с учетом физических факторов, свойственных климатическим условиям южного бассейна озера Байкал.
Натурные исследования показывают, что основные изменения происходят в направлении, перпендикулярном фронту, то есть от берега к центру. При этом характеристики в направлении, параллельном берегу, достаточно однородны. На этом основании полагают, что модель, в которой исключены все градиенты в направлении, параллельном берегу, должна качественно и правильно описывать физический процесс, а явление термического бара считают квазидвухмерным [7]. Модель, представленная в данной работе, также использует гипотезу о двухмер-ности, при этом не учитывается влияние силы Кориолиса, связанной с вращением Земли. Рассматривается взятый из работы [10] прибрежный профиль озера, соответствующий реальным условиям южного бассейна озера Байкал (рис. 1). Протяженность расчетной области (Ьх = 10 км) намного больше глубины, а глубина (Н = 900 м) примерно соответствует средним глубинам южного бассейна Байкала [11]. В качестве уравнения состояния используется уравнение, связывающее плотность воды с температурой. На свободной поверхности ставятся граничные условия типа «твердой крышки» (отсутствие ветровых напряжений) и задается поток тепла. На дне, помимо условия непроницаемости, задается связь касательных напряжений с придонной скоростью (квадратичный закон трения), а также условие отсутствия теплообмена с дном. Начальные условия соответствуют состоянию покоя и заданным полям температуры (температура однородна по горизонтали, но переменна по вертикали), взятым на основе натурных наблюдений в мае месяце [10] (рис. 2).
АТМОСФЕРА
У, м 800
600
400
200
Рис. 1. Профиль южного бассейна озера Байкал
Рис. 2. Начальное распределение температуры
Математическая постановка задачи основана на двумерной негидростатической модели в приближении Буссинеска для конвективного течения:
.2
ди ди дыу
ді дх ду
1 др д ( ди \ д ( ди
------— +—I К — 1+—I К —
Р0 дх дх V х дх) ду V у ду
0
dv+duv+dv2^.di+—{K dv)+—(K dvV -(р-ро).
дґ дх ду р0 ду дх V дх) ду V у ду) р0
ди ду л — + — = 0;
дх ду
дТ + диТ + дуТ = д ( дТ А 5 ( дТ
дґ дх ду дх Vх дх) ду Vу ду где р = p + gyp0'; р = р0 (-у(Т-Тт )2), ) = 8.572628-10-6 град-2; и, V - составляющие скорости по осям х и у соответственно, Т - температура, Тт и 4 °С - температура максимальной плотности (ТМП), р0 = 998,2 кг/м3 - характерная плотность воды, р - давление, g - ускорение свободного падения, коэффициенты Кх, Ку и Бх, Бу характеризуют интенсивность диффузионного переноса импульса и тепла в соответствующем направлении и рассчитываются по формуле Обухова [12]:
Ку = (0,05£0 )2^(|у] ~р~ду (формула Обухова);
Кх = Ку; Бх = 0,7Кх; Бу = 0,7К ; Ь0 - характерный размер области перемешивания.
Начальные условия задаются в виде
ґ = 0: и = V = 0; Т = Т0(у).
Профиль Т0 (у) представлен на рис. 2.
Граничные условия имеют вид
дТ
и = 0; V = 0; — = 0 - на дне;
дп
дТ
х = Ьх : и = 0; V = 0; — = 0 - на правой границе;
х дх
г du n n п дТ
У = Ly : — = 0; v = 0; cwaterРоDy — = Q, - га шверхшсш, ду ду
где Q = 210 Вт/м2 - поток тепла через свободную поверхность, cwater - теплоемкость воды.
Решение конвективно-диффузионных уравнений основано на конечно-разностном методе конечного объёма [13]. Перед применением метода конечного объема для получения дискретизации исходной системы строится шахматная сетка: давление и температура определяются в узловых точках построенной сетки, а компоненты скорости рассчитываются в точках, расположенных на гранях конечного объема. Дифференциальное уравнение интегрируют по каждому конечному объему. Следующий шаг состоит в аппроксимации полученных интегральных соотношений конечно-разностными, и в результате получается совокупность разностных уравнений.
Численный алгоритм нахождения поля течения и температуры опирается на разностную схему Кранка - Николсона. Конвективные слагаемые в уравнениях аппроксимируются по противопотоковой схеме Леонарда Quick 2-го порядка. Для согласования поля скорости и давления использована процедура SIMPLE Патан-кара [13]. Метод SIMPLE (Semi-Implicit Method for Pressure-Linked Equations - по-
лунеявный метод для связывающих давление уравнений) основан на циклической последовательности операций «предположение - коррекция» при решении уравнений. Используя некоторое начальное поле давления, сначала вычисляют компоненты скорости по уравнению движения. Затем давление и компоненты скорости корректируются так, чтобы удовлетворять уравнению неразрывности. Этот процесс продолжается до тех пор, пока итерационный процесс не сойдется.
Системы разностных уравнений на каждом шаге по времени решаются методом нижней релаксации или явным методом Булеева. Для вычисления компонент скорости и, V применяется метод релаксации. Для нахождения температуры и давления используется явный метод Булеева.
Тестовые эксперименты проводились на классической задаче о тепловой гравитационной конвекции в квадратной каверне с использованием результатов исследований В. И. Полежаева и др. [14, с. 26-31]. Расчетная область покрывалась прямоугольной сеткой 65x65. Начальное условие расчета - покой. В двумерной негидростатической модели в приближении Буссинеска коэффициенты интенсивности диффузионного переноса импульса и тепла положим равными
где V - коэффициент кинематической вязкости, а - коэффициент температуропроводности.
Вводя масштабы для искомых величин и независимых переменных, исходную систему уравнений в приближении Буссинеска можно привести к безразмерному виду. Тогда в таких уравнениях возникают безразмерные параметры: вг = gpЬ3Д7Уv2 - число Грасгофа (в - коэффициент теплового изменения плотности, Ь = Ьх = Ьу), определяющее интенсивность тепловой конвекции; Рг = v/a -число Прандтля, представляющее отношение толщин динамического и теплового пограничных слоев; Яа = вгхРг - число Рэлея. В данных расчетах Рг=0,7, Яа=104. Эксперименты, проведенные в замкнутой квадратной области при изотермических боковых границах и с условием теплоизоляции горизонтальных границ, показывают, что полученные результаты совпадают с картиной течения, приведенной в работе [14] (рис. 3 и 4).
Кх Ку V, Бх Бу а,
1
0.8-
0.6-
0.4-
0.2-
0 0.2 0.4 0.6 0.8 1
Рис. 3. Изотермы, полученные в ходе Рис. 4. Изотермы, приведенные в [14, с. 27]
вычислительного эксперимента
Результаты тестовых расчетов показывают, что физическая картина течения в каверне состоит в том, что нагреваемые у стенки слои жидкости поднимаются вверх, а охлаждаемые опускаются вниз, в связи с чем в прямоугольной области образуется циркуляционное течение, перемешивающее жидкость и переносящее тепло от нагретой стенки к холодной [15].
В численных расчетах термобара сделана попытка воспроизвести реальные условия озера. Путем блокировки (делая равными нулю скорости в выключенной зоне за счет использования очень больших значений коэффициента вязкости в этой зоне и задавая нулевое значение скорости на фиктивной границе) некоторых конечных объемов прямоугольной неравномерной сетки расчетная область была приближена к прибрежному профилю озера. Используется неравномерная ортогональная сетка 126x90 с измельчением шагов у берега (х = 0). В результате Их меняется от 25 до 200 м, Ну = 10 м. Шаг по времени Дґ = 60 с. Начальные условия соответствуют состоянию покоя и заданным полям температуры (температура однородна по горизонтали, но переменна по вертикали), соответствующим измеренному его вертикальному распределению в мае, которое изменялось от 1,7 градусов на водной поверхности до максимального значения 3,8 на глубине 900 м [10]. На поверхности задан поток тепла Q = 210 Вт/м2 [16]. Ветровое трение отсутствует.
На рис. 5 - 8 представлены результаты расчетов линий тока, проведенных в условиях, близких к естественным. На рис. 5 показано образование циркуляционного течения вблизи берега, которое с течением времени продвигается к открытой границе (рис. 6 и 7), и после 30 суток от начала вычислительного эксперимента система циркуляционных течений за счет явления термобара продвигается на расстояние более 9 км от берега (рис. 8).
Рис. 6. Линии тока через 15 суток
Рис. 8. Линии тока через 30 суток
Отметим, что рисунки здесь представлены в искаженном масштабе: отношение масштаба по горизонтали к масштабу по вертикали на рисунках уменьшено в 10 раз по сравнению с реальным их отношением для рассматриваемого участка озера, и поэтому соответственно в 10 раз уменьшена в них горизонтальная компонента скорости.
В процессе постепенного перехода воды в бассейне через ТМП возможны три фазы развития водообмена [17]:
1) формирование вдольсклонового потока;
2) образование подповерхностной струи;
3) трансформация одного типа циркуляции в другой при переходе через ТМП.
Как видно из рисунков, за промежуток времени 30 суток за счет явления термобара поверхностные воды проникают примерно до глубин свободной конвекции 150 - 250 м и затем имеют тенденцию продвигаться в горизонтальном направлении в сторону центральной части озера.
Наблюдается хорошее качественное соответствие структуры течений, полученной в ходе численного моделирования, с результатами натурных наблюдений [10]. В области, где температура воды ниже ТМП, наблюдаются вдольсклоновые потоки; в области, где температура воды выше ТМП, существует активная подповерхностная струя теплых вод, в среднем слое наблюдается компенсационное течение, направленное к берегу. Опускание вод происходит между ячейками, а компенсирующий подъем - на периферии, причем подъем в прибрежной ячейке даже нарушает общее направление движения вниз по склону. Модель наглядно демонстрирует, что глубоководные слои вовлечены в вертикальное перемешивание в течение долгого времени после прогрева верхнего слоя - до тех пор, пока хоть какая-то часть вод с температурой, меньшей ТМП, находится над наклонным дном
и прогревается сверху, получая тепло не от воздуха, а от верхнего тёплого слоя воды. В частности, это может быть причиной порождения холодноводных придонных интрузий, недавно обнаруженных при проведении натурных исследований на озере Байкал: авторы связывают их происхождение с наличием термобара, однако «реальный механизм водообмена остаётся неясным» [17].
Типичные абсолютные значения скоростей, полученные в численных экспериментах, были следующими: компоненты и - 0,08-0,46 см/с, v - 0,09-0,45 см/с. Интересно отметить, что вертикальная и горизонтальная компоненты скорости изменяются в почти одинаковом диапазоне.
Приведенные результаты показывают, что в условиях, близких к естественным, происходит проникновение поверхностных вод в глубокие слои при весеннем прогреве. Это согласуется с описаниями натурных наблюдений [10].
Полученные результаты математического моделирования и построенный вычислительный алгоритм являются фундаментальными исследованиями и могут найти широкое применение в изучении механизмов, управляющих процессами перемешивания и водообмена в природных водоемах, являющегося важным направлением гидрофизических исследований. Результаты исследований, полученные для озера Байкал, могут быть также использованы рыбопромысловыми организациями для предсказания районов высокой концентрации рыбы, приуроченных к существованию подобного фронта в весенний период в озере. Выводы представляют интерес при решении таких важных практических задач, как распространение загрязнений в водоемах и рациональное использование их природных ресурсов.
ЛИТЕРАТУРА
1. Афанасьева Э.Л., Бекман М.Ю., Безрукова Е.В. Путь познания Байкала. Новосибирск: Наука, 1987. 301 с.
2. Walker S.J., Watts R.G. A three-dimensional numerical model of deep ventilation in temperate lakes // J. Geophys. Res. 1995. V. 100. P. 22711-22731.
3. Killworth P.D., Carmack E.C., Wiess R.F., Mateas R. Modelling deep-water renewal in Lake Baikal // Limnol. Oceanogr. 1996. V. 41. No. 7. P. 1521-1538.
4. Цветова Е.А. Математическое моделирование Байкальского термобара // Математические проблемы экологии: тр. второй Всерос. конф. Новосибирск, 1994. C. 44-49.
5. Tsvetova E.A. Convective currents assotiated with the thermal bar of Lake Baikal // Advanced Mathematics: Coputations and Applications. NCC Publisher, 1995. P. 386-393.
6. Цветова Е.А. Специфичесие проявления конвекции в глубоких озерах // Математические проблемы экологии: тр. Третьей Междунар. конф. Новосибирск, 1996. C. 181-189.
7. Цветова Е.А. Численная модель термобара в озере Байкал // Метеорология и гидрология. 1997. № 9. C. 58-68.
8. Квон В.И., Квон Д.В. Численный анализ механизма глубокого проникновения поверхностных вод в прибрежной зоне озера в период весенне-летнего термобара // Вычисл. технологии. 1997. Т. 2. № 5. C. 46-56.
9. Овчинникова Т.Э., Бочаров О.Б. Сезонное влияние вод притока на водообмен в глубоком озере в условиях больших уклонов дна // Вычисл. технологии. 2007. Т. 12. № 6. C. 59-72.
10. Shimaraev M.N., Verbolov V.I., Granin N.G., Sherstyankin P.P. Physical Limnology of Lake Baikal: a Review. Irkutsk - Okayama, 1994. 81 с.
11. Бочаров О.Б., Овчинникова Т.Э. Численное моделирование явления термобара в озере Байкал // Вычисл. технологии. 1996. Т. 1. № 3. C. 21-28.
12. Белолипецкий В.М., Шокин Ю.И. Математическое моделирование в задачах охраны окружающей среды. Новосибирск: ИНФОЛИО-пресс, 1997. 240 с.
13. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. -М.: Энергоатомиздат, 1984. С. 96-98.
14. Полежаев В.И., Бунэ А.В., Верезуб Н.А., и др. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье - Стокса. М.: Наука, 1987. 271 с.
15. Цыденов Б.О. Численное моделирование конвективных течений в каверне // Перспективы развития фундаментальных наук: тр. VI Междунар. конф. Томск: Изд-во ТПУ, 2009. Т. 2. С. 673-676.
16. Блохина Н.С., Соловьев Д.А. Влияние ветра на динамику развития термобара в период весеннего прогрева водоема // Вестник Моск. ун-та. 2006. Сер. 3. № 3. C. 59-63.
17. Демченко Н.Ю. Исследование структуры и динамики термобара в пресных и солоноватых водоемах: автореф. дис. ... канд. физ.-мат. наук. Калининград, 2008. 24 с.
Статья принята в печать 28.01.2011 г.
Tsydenov B.O., Starchenko A.V. NUMERICAL MODELLING OF THE THERMAL BAR EFFECT IN LAKE BAIKALIN A SPRING-SUMMER WARMING PERIOD. The developed mathematical model and numerical algorithm make it possible to reproduce the process of the thermal bar dynamics in Lake Baikal with allowance for main physical factors. The results presented here show the appearance of a circulation flow off shore and its movement towards the central part of the lake with time.
Keywords: thermal bar, temperature of maximum density, Boussinesq approximation, convection, numerical experiment.
TSYDENOV Bair Olegovich (Tomsk State University)
E-mail: ba1r@sibmail.com.
STARCHENKO Alexander Vasilievich (Tomsk State University)
E-mail: starch@math.tsu.ru
