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

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

Том 3, № 4, 1998

О ТЕРМОГРАВИТАЦИОННОЙ КОНВЕКЦИИ В ПРИБРЕЖНОЙ ЗОНЕ ГЛУБОКОГО ОЗЕРА В ПЕРИОД ВЕСЕННЕ-ЛЕТНЕГО ПРОГРЕВАНИЯ *

О. Б. БОЧАРОВ, Т. Э. ОВЧИННИКОВА Институт водных и экологических проблем СО РАН Новосибирск, Россия e-mail: [email protected]

The formation and evolving of convective currents in a deep lake in a spring-summer warming period is modelled numerically. The calculations are performed under the conditions approximating the natural situation of the South Baikal. The effect of the water state equation on the development of global vertical circulations in the lake is investigated.

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

В предыдущих работах авторов [1, 2] с помощью двумерной гидротермической модели несжимаемой жидкости в приближении Буссинеска было показано, что действительно в мае — июне формируются глобальные циркуляции, достигающие дна. Однако использовалось уравнение состояния Кнудсена р = р(Т), не учитывающее влияния давления на плотность, в то время как натурные данные показывают, что температура максимальной плотности воды Tm зависит от глубины, что, вообще говоря, может существенно повлиять на картину течений.

1. Анализ уравнения состояния

На рис. 1 приведены зависимости Tm от давления (1 бар ^ 10 м водного столба) по натурным данным (А. Б. Каплун, M. K. Штрем [3]) для уравнения состояния Объединенной комиссии ЮНЕСКО по океанографическим таблицам и стандартам [4] и для уравнения

*По материалам доклада, сделанного на конференции "Математические модели и численные методы механики сплошной среды" (Новосибирск, Академгородок, 27 мая - 2 июня 1996 г.), посвященной памяти академика Н.Н. Яненко.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-05-66101.

© О. Б. Бочаров, Т. Э. Овчинникова, 1998.

Рис. 1. Температура максимальной плотности: 1 — по данным Гилла, Чена — Миллеро, 2 — Каплуна, 3 — Штрема.

состояния, часто используемого лимнологами [5]. Столь различные результаты свидетельствуют о том, что вопрос о зависимости температуры максимальной плотности от давления требует дополнительного исследования.

Нами проведены численные эксперименты с квадратичным уравнением состояния Мар-кофски — Харлемана, в которое была введена экспериментальная функция Тт (г) по данным А. Б. Каплуна

р = 999.975[1 - 6.8 • 10-6 (Т - Тт (г ))2 ].

Использовались также и полные уравнения р(Т, р, 5) Гилла [4] и Чена —Миллеро [5] с постоянной соленостью 5.

Р(Т)-Р0

ал -

0.2 .................................................. гр 0р

0 2 4 б 8

Рис. 2. Зависимость плотности от температуры: 1 — по данным Гилла, 2 — Чена —Миллеро, 3 —Харлемана; ро = 999.6, S = 0, р = 5.

Сопоставляя профили плотности в диапазоне температур 0 — 10 °С и давлений 0 — 150 бар, что почти соответствует максимальным глубинам о. Байкал и диапазону температур в мае — июле, можно увидеть, что в наиболее интересной для нас области параметров уравнения [4, 5] дают мало различимые результаты. Если в уравнение Маркофски — Харлемана не вводить зависимость температуры максимальной плотности от давления, то на поверхности все уравнения дают близкие значения, а с глубиной расхождение с [4, 5] растет. Для примера на рис. 2 приведены графики плотности при давлении 5 бар. Следует заметить, что в окрестности температуры максимальной плотности уравнение состояния можно представить следующим образом [6]:

Р = Рт(р){1 — S(р)[Т — Тт(р)]2} . (1)

В таком виде оно удобно для качественного анализа устойчивости равновесия стратифицированной жидкости.

2. Математическая модель и численный алгоритм

В основу исследований положена известная двумерная вертикальная модель расчета осред-ненных характеристик турбулентного течения, получаемая из уравнений Навье — Стокса — Рейнольдса с использованием гипотезы Буссинеска для напряжений и приближения Обе-бека — Буссинеска при учете плотностной неоднородности. Достаточно распространенная гипотеза о гидростатическом распределении давления не использовалась. Система уравнений при этом имеет следующий вид [7]:

дv д д 1 _ д ^ д д ^ д р

— + — п\ + — -V + — Ур = — Кх— V + — Кг—V + — g,

дъ дх дг р0 дх дх дг дг р0

дд

—и + — и = 0, (2)

дх дх

дТ + д + д д^д + д п дг

— + — иТ + — иТ = — Вх— Т + — и2 — . дъ дх дг дх дх дг дг

Р = Р(Р,Т )•

Здесь и, и — компоненты вектора скорости V вдоль осей координат х и г (направленных перпендикулярно берегу и по вертикали вверх соответственно), Т — температура, р0 — характерная плотность воды, коэффициенты х, Кх и Ох, характеризуют интенсивность турбулентного переноса импульса и тепла в соответствующем направлении (Д^ = 1.4КЯ,

q = х, г).

Граничные и начальные условия задаются следующим образом. На свободной поверхности г = Н ставится кинематическое граничное условие типа абсолютно гладкой неподвижной "твердой крышки" и задается атмосферное давление

ди

и = 0 , Кх — = 0 , р = Ра(х,Ъ) ,

а также задается поток тепла Ф, зависящий от потока солнечной радиации, температуры воздуха и воды, облачности и влажности, рассчитываемый по методике, описанной в [8]:

дТ

= ф.

дг

На дне г = гь(х) задаются касательные напряжения, пропорциональные квадрату скорости, и ставится условие отсутствия теплообмена с грунтом

0ут _ | дТ 0

—— = -КЬТ | _ I , —— = 0 .

дИк дИл

Здесь vТ — касательная составляющая скорости, к — коэффициент трения о дно, д/дЫа, (а = к,й) — производная по конормали, определяемая соответствующим тензором турбулентного переноса:

д д д Кх ео8(_,,х) ——+ К ео8(п,

ONk 'Ox 'Ox

n — вектор внешней нормали к границе области.

Что касается граничных условий при x = 0 (у берега) и при x = L (открытая граница), то в работе [1] анализировались различные их виды. Здесь мы остановились на варианте, который представляется нам наиболее физичным и дающим наибольшую устойчивость счета, — твердая стенка при x = 0 и условие симметрии вихрей при x = L:

x = 0: u = 0, w = 0, Tx = 0;

x = L : u = 0, wx = 0, Tx = 0 .

Начальные условия при t = 0 соответствуют состоянию покоя u = w = 0 и заданному полю температуры T0(x,z), по которому определяется гидростатическое распределение давления с учетом уравнения состояния.

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

Численный алгоритм определения поля скорости основан на так называемом Ф-про-екционном методе Чорина — Тимухина [9], который в общих чертах выглядит следующим образом. На первом шаге рассчитывается вспомогательное поле скоростей v* из уравнений с исключенным градиентом давления:

dv* д д д д д д о

^ + дик* + -^wv* = дКхд v* + v* + Р g. (3)

дt ox Oz дx дx Oz дz p0

По этому полю определяется завихренность

rot v* = rot v = и = u*z — w* . Затем функция тока определяется из уравнения

ДФ = и , (4)

а окончательное соленоидальное поле скоростей — из уравнений связи по функции тока:

дФ дФ

u = -7Г- , w = —— .

OZ OX

В работе [1] в ходе численных экспериментов с задачами рассматриваемого типа было установлено, что относительное отклонение давления от гидростатического, характеризуемое величиной

p0 (Ow д д д Ow д Ow\

— + Т" uw + ТГ- uw — — K^---— Kz— , (5)

pg \ Ot Ox Oz Ox Ox Oz Oz J

имеет порядок 10-9 и только в узкой области нисходящего течения достигает порядка 10-4. Следовательно, максимальное по области отклонение давления от гидростатического составляет сотую долю процента. В связи с этим для упрощения алгоритма при расчете плотности по уравнению состояния р = р(р, Т) давление определялось из гидростатического уравнения

1 = —р(р,Т ^

Аппроксимация уравнений выполнена на неравномерной ортогональной сетке с измельчением шагов у свободной поверхности (г = Н) и берега (х = 0). В результате Нг меняется от 1 до 33 м, а кх — от 33 до 1350 м. Используются консервативные неявные разностные схемы. Конвективные члены аппроксимируются с помощью разностей против потока, операторы второго порядка — с помощью центральных разностей. В целом алгоритм имеет первый порядок аппроксимации как по времени, так и по пространству. Линеаризация конвективных членов осуществляется путем использования значений скоростей с предыдущего временного слоя. Алгебраические системы решаются методом верхней релаксации. Подробный анализ схем такого класса для задач конвекции проведен в [10]. Здесь хотелось бы подчеркнуть, что последовательный расчет скорости и температуры привносит в алгоритм свойства явных схем. В результате при больших шагах по времени (более 300 с в нашем случае) возникает неустойчивость, проявляющаяся в развитии хаотических вихревых образований. Подходящий шаг по времени при выполнении расчетов подбирался экспериментально из условий приемлемой скорости сходимости итерационных процессов (20-30 итераций на шаг) и стабильности характеристик течения при изменении шага.

3. Численные эксперименты

В численных расчетах мы попытались по возможности воспроизвести реальные условия озера Байкал. Прибрежный профиль озера взят из работы [11], а глубина Н = 900 м примерно соответствует средним глубинам южного бассейна Байкала. Протяженность расчетной области Ь = 10 км. Начальное вертикальное распределение температуры, график которого приведен на рис. 3, основано на данных натурных наблюдений в южной части Байкала в июне [11]. Следует заметить, что от поверхности до глубины ^ 25 м температура падает (Тг > 0), далее до глубины ^ 300 м растет (Тг < 0), а ниже отметки 300 м до самого дна опять идет медленное падение температуры (Тх > 0). Метеоусловия, соответствующие маю — июню, взяты из работы [12]. В расчетах использовались натурные данные для коэффициента вертикального турбулентного обмена [13] (рис. 4) а горизонтальная турбулентная вязкость определялась по закону 4/3 Ричардсона: Кх = Кг (Нх/Нг)4/г.

В случае зависимости р(р, Т) условие гидростатической устойчивости имеет вид [4]

—^ I+£ * 0 , с2=* (6)

Медленность вертикальных движений в океане и крупных водоемах часто позволяет применять гидростатическое приближение для давления. Применимость этого приближения в изучаемой проблеме была обоснована численными расчетами в [1]. В гидростатическом приближении условие устойчивости (6) приобретает вид

Г = > 0 (7)

Л = — дТдг > 0 • (7)

Рис. 3. Начальный профиль температуры: 1 — Т(К), 2 — Тгт по данным Миллеро, 3 — Каплуна.

900

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

При несколько других предположениях аналогичное условие устойчивости было получено и в [14]. Применяя для уравнения состояния представление (1), из (7) получим условие устойчивости:

дТ

и • 104 = 2рт (р)Б (р)[Т - Тт (р)] ^ > 0.

В силу положительности функций рт (р) и 5 (р) оценка устойчивости теперь сводится к анализу отклонения температуры от Тт (р) и знака дТ/дг.

На рис. 5. приведены графики поведения функции устойчивости /8(г), рассчитанные по начальному распределению температуры и трем уравнениям состояния. Видно, что верхняя зона до глубины порядка 30 м обладает неустойчивой стратификацией, далее до глубины 200 м имеет место устойчивое распределение температуры. Ниже 200 м, с учетом

■0,3 -0.2 0,0 0.2 0.3 0,5

_1М 1Ц.ИМ -Ш-М /V / ......

! || 1 и 3 ! |

н ¡1 л

! -1 ..!________1

! 1 ! 1 : 1 1 1 ; | /

! 1 ; 1 ! 1

Рис. 5. Функция устойчивости: 1 — по данным Гилла, 2 — Чена — Миллеро, 3 —Харлемана.

Рис. 6. Динамика функции устойчивости: 1 — 10-е, 2 — 20-е, 3 — 30-е сутки.

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

Расчеты, выполненные с использованием трех уравнений состояния, привели к качественно одинаковым результатам. На рис. 6. показана динамика устойчивости стратификации для уравнения состояния Чена — Миллеро [5]. Сравнивая с начальным распределением (см. рис. 5.), можно заключить, что на первой стадии прогревание и конвективное перемешивание приводят к понижению устойчивости верхних слоев (¡8 < 0 до глубины

Рис. 7. Линии тока на 10-е, 20-е, 30-е сутки при Ф = —0.001 (1), 0.001 (2).

примерно 300 м и по модулю существенно меньше, чем раньше). В дальнейшем устойчивость верхних слоев повышается, и стратификация начинает приближаться к характерному летнему распределению с формированием сильно устойчивой зоны термоклина. Устойчивость глубинных слоев (ниже 400 м) практически не меняется. На рис. 7. приведены картины линий тока, рассчитанные с использованием уравнения состояния Чена — Миллеро [5] через 10, 20 и 30 суток от начала расчетного периода, причем первая из них показывает в укрупненном виде прибрежную зону.

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

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

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

[2] Бочаров О. Б., Васильев О.Ф., Квон В. И., Овчинникова Т. Э. Математическое моделирование термобара в глубоком озере. Докл. РАН, 349, №4, 1996, 530-532.

[3] Проблемы Байкала. Труды ЛИН, 16, №36, Наука, Новосибирск, 1978.

[4] Гилл А. Динамика атмосферы и океана. Т. 2, Мир, М., 1986.

[5] CHEN C. T., Millero F. J. Precise thermodynamic properties for natural waters covering only the limnologies range. Limnol. Oceanogr. 31, No. 3, 1986, 657-662.

[6] Gebhart B., Mollendorf J. C. A new density relation for pure and saline water. Deep Sea Res., 24, 1977, 831-848.

[7] Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная механика и теплообмен. Т. 2, Мир, М., 1990.

[8] Архипов Б. В., СолБАков В. В. Расчет термогидродинамического режима водоема по двумерной модели. Изв. АН ФАО, 30, №5, 1994, 671-685.

[9] Тимухин Г. И. О построении некоторых схем приближенного решения уравнений динамики вязкой несжимаемой жидкости. Числ. методы механики сплошн. среды. 3, №2, 1972, 77-95.

[10] Полежаев В. И., Бунэ А. В., Верезуб Н. А. и др. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье — Стокса. Наука, М., 1987.

[11] Shimaraev M. N., Verbolov V. I., Granin N.G., Sherstyankin P.P Physical Limnology of Lake Baikal: a Review. Irkutsk — Okayama, 1994.

[12] Верболов В. И., Сокольников В.М., Шимараев М. Н. Гидрометеорологический режим и тепловой баланс озера Байкал. Наука, М.—Л., 1965.

[13] ШиМАРАЕВ М. Н. Элементы теплового режима озера Байкал. Наука, Новосибирск, 1977.

[14] КАМЕНКОВИЧ В. М. Основы динамики океана. Гидрометеоиздат, Л., 1973.

Поступила в редакцию 8 декабря 1997 г.