CC BY
24
Серия «Математика»
Том 1 (2007), № 1, С. 42-51
Онлайн-доступ к журналу: http://isu.ru/izvestia
ИЗВЕСТИЯ
Иркутского государственного университета
УДК 556.55
Численное моделирование течений и переноса примесей в стратифицированных озерах
В. К. Аргучинцев, А. В. Аргучинцева (arg@math.isu.ru)
Иркутский государственный университет
Аннотация. В данной работе предлагается негидростатическая нестационарная трехмерная нелинейная модель для описания течений, полей температуры и плотности в озерах и водохранилищах с учетом неоднородностей рельефа дна и использованием уравнения состояния, связывающего давление, температуру и плотность. Система уравнений гидротермодинамики решается численно с применением метода фиктивных областей. Используется неявная конечно-разностная схема покомпонентного расщепления. Приводятся результаты расчетов преобладающих течений в верхних слоях озера Байкал в безледный период. Найденные на основе гидротермодинамической модели скорости движения и турбулентные характеристики используются для расчета переноса примесей от Байкальского целлюлозно-бумажного комбината.
Ключевые слова: моделирование, уравнения гидротермодинамики, озеро, течения, примесь
1. Введение
Существующая сеть гидрологических наблюдений слишком редка для экспериментального изучения мезомасштабных процессов в водоемах. Одним из перспективных методов их изучения является математическое моделирование, которое позволяет объяснять с теоретической точки зрения количественные и качественные закономерности, решать диагностические и прогностические задачи. Заметим, что большинство работ посвящено двумерным моделям. Однако мезомасштабные процессы в гидросфере почти всегда существенно трехмерны. Поэтому для практического применения нужны пространственные модели, разработка которых наталкивается на трудности, связанные с необходимостью построения эффективных вычислительных алгоритмов для решения сложной системы уравнений гидротермодинамики (Цветова, 1974, 1977;
Марчук, Кочергин, Цветова, 1978; Квон, 1979; Гурина, Демин, Филатов, 1984; Астраханцев, Руховец, 1986; Brugge, Jones, Marshall, 1991; Walker, 1994; Walker, Watts, 1995; Пененко, Цветова, 1998). Одним из способов упрощения уравнений гидротермодинамики является фильтрация волн, скорость распространения которых велика по сравнению со скоростями основных процессов. Наличие таких волн затрудняет разработку эффективных численных методов интегрирования уравнений. Применение гидростатического приближения позволяет отфильтровать звуковые волны. В настоящее время большинство созданных трехмерных моделей основано на гидростатическом приближении. Однако гидростатический закон выполняется с достаточной степенью точности, если горизонтальный масштаб изучаемых движений существенно преобладает над вертикальным.
Разработка негидростатических моделей наталкивается на трудности, связанные с необходимостью построения эффективных вычислительных алгоритмов. При использовании негидростатических уравнений для отфильтровывания звуковых волн обычно пренебрегают локальной производной от плотности в уравнении неразрывности. При этом искажаются основные типы движений и усложняется метод решения из-за отсутствия эволюционного уравнения для давления.
В этой статье предлагается негидростатическая модель с учетом сжимаемости воды, которая является более универсальной для описания мезомасштабных процессов в стратифицированных водоемах.
2. Основные уравнения
В качестве основы математического моделирования мезомасштабных процессов в стратифицированных водоемах примем уравнения геофизической гидродинамики, выражающие основные законы сохранения энергии, импульса и массы стратифицированных сплошных сред.
В систему дифференциальных уравнений нестационарной трехмерной нелинейной модели включаются: уравнение движения —
йу 1
— = — gradp —2ш х у + д + Бу, (2.1)
уравнение неразрывности
dt + р div v = 0 (2.2)
уравнение притока тепла
§ - ^ = DT + MT (2.3)
dt Cpp dt
уравнение переноса солености
йд
^ = Бд + И,, (2.4)
уравнение состояния, записанное в общем виде —-
р = р(р,Т,д), (2.5)
полуэмпирическое уравнение турбулентной диффузии примеси —-
дв дв — + div(вV) - Шд — = Бв — а в + К + Г, (2.6)
дЬ дг
соотношение А.М. Обухова (1946, 1988)
, I (0.05п)2\/С + ко, если С > 0 , .
к* = п / п (2.7)
где С = (&) + (&) + ,
соотношение Смагоринского (Smagorinsky, 1963)
2 / (дт ди\2 /дм д^2 . .
где
й д д д д ^ , д . дф д . дф д . дф "гг = ог+ и— + + ш—; Бф = — к^— + — к^— + — кг — йЬ дЬ дх ду дг дх дх ду ду дг дг
Здесь ф — любая из функций рассматриваемой задачи; Ь — время; и,у — горизонтальные, ш — вертикальная компоненты вектора V скорости движения среды вдоль осей декартовой прямоугольной системы координат (х,у,г), х,у — горизонтальные координаты, а ось г направлена вертикально вверх; р — плотность воды; р — давление; Т — температура воды; и — вектор угловой скорости вращения Земли, направленный параллельно оси Земли к Северному полюсу; д — сила тяжести; ср — удельная теплоемкость при постоянном давлении; а = — 1 дт — коэффициент термического расширения; д — соленость; Ит — скорость притока тепла от внешних источников; И, — интенсивность источников или стоков субстанции; кх,ку,кг — коэффициенты турбулентного обмена соответственно по горизонтали и вертикали; т — масштаб, пропорциональный шагу горизонтальной сетки Дв; п — толщина верхнего квазиоднородного слоя, ко — эмпирическая константа, в -— концентрация рассматриваемой субстанции, а и К учитывают трансформацию и взаимодействие различных субстанций между собой; Г(х, у, г, Ь) — функция, описывающая распределение и мощности
источников рассматриваемой субстанции; wg — скорость гравитационного осаждения, определяемая по формуле Стокса (Фукс, 1955): wg = 2pngr2/(9pv), где pn и r — соответственно плотность и радиус частицы взвеси; v — вязкость воздуха.
Для воды используется эмпирическое уравнение состояния (Chen, Millero, 1986), связывающее плотность, температуру, давление и соленость:
p(q, T, p) = po(q, T)/(1 - p/K(q, T,p)),
где po(q, T) — плотность воды при стандартном атмосферном давлении, K — объемный модуль упругости.
Систему уравнений (1)-(8) будем рассматривать в параллелепипеде Ü{x,y,z : 0 < x < X, 0 < y < Y, 0 < z < h}, где x = 0, x = X,y = 0,y = Y — границы области счета. Преобразуя уравнения (2)-(3) с помощью уравнения состояния, получим эволюционные уравнения для T и p. На свободной поверхности водоема задаются касательное трение ветра, поток солнечной радиации и теплообмен водной поверхности с атмосферой. На дне водоема и боковой поверхности ставятся условия прилипания или непротекания с квадратичным законом трения и задается теплообмен с дном. На границах втекания ставятся условия первого рода, а на границах вытекания — второго рода.
Из-за отсутствия необходимой информации о гидрометеорологических полях начальные условия заменяются значениями, полученными на основе крупномасштабной модели или находятся из решения соответствующих стационарных задач.
3. Метод решения
Уравнения модели интегрируются численно в декартовой прямоугольной системе координат с применением метода фиктивных областей. Введение таких областей позволяет проводить расчеты с произвольной функцией, описывающей рельеф суши и дна водоемов.
Сложность решения рассматриваемой системы уравнений обусловлена наличием физических процессов с различными характерными временными масштабами. Поэтому численный алгоритм решения уравнений строится на основе метода расщепления по физическим процессам (Марчук, 1974).
Решение задачи на каждом временном шаге осуществляется в два основных этапа: 1) перенос субстанций вдоль траекторий и турбулентный обмен; 2) процесс согласования гидрологических полей. Такой подход позволяет в принципе использовать разные шаги по времени на каждом этапе.
На первом этапе для каждой искомой функции рассматривается эволюционное уравнение вида:
% + ^ = <•■
где Ь = ^2^=1 Ьт-
Редукция сложных задач к более простым возможна в тех случаях, когда исходный положительно полуопределенный оператор задачи представим в виде суммы положительно полуопределенных простейших операторов.
Аппроксимация по времени строится на основе метода покомпонентного расщепления по геометрическим переменным (Марчук, 1974), который состоит в разложении сеточного оператора Ь1 > 0 на более простые операторы Ь^п > 0. Операторы Ь^п > 0 аппроксимируем со вторым порядком точности по координатам.
Введем неравномерную сетку с основными узловыми точками
xi = гАх (г = 0,1,...,1 + 1); у = ЦАу (з = 0,1,...,,] + 1);
гк = кАгк (к = 0,1,...,К + 1); Ьп = иАЬ Ц = 0,1,...)
и шагами сетки Ах, Ау, Ахк, АЬ.
Будем также использовать вспомогательные точки х+\/2,У3+1/2, Хк+\/2, расположенные в серединах основных интервалов. Обозначим:
Ф?,3, к = Ф^'Уз ,гк ,Ьи),
Щ±1/2,3, к = (Щ±1/2,3, к + Щ ,3 , к)/2, (г = 1, 2,...,1),
щ ,з±1/2, к = Щ ,з±1/2, к + Vi ,3, к )/2, (з = 1, 2,...,]),
Щ,3,к±1/2 = (™^з,к±1/2 + wi,j,k)/2,
Ак = (Агк+1 +Агк )/2 (к = 1, 2,..., К).
Приведем разностные аналоги операторов:
Ь'ФХзк =
= уП+1 /23к и1-1 1,3,к 3,к = 2Ах
АХ [к**+1/2 *. к (^+1,3,к - ) - ,, к к (Аз,к - ф-1,з,к )],
ГтН 1\ 3+1/2,к'Фь3+1,к - 3-1/2,к^,3-1,к (Ь2 У)%3,к = -
2Ау
А? К,к (^,3+1,к - А3,к) - кк - Ф^-1,к)],
^3к+1/2Фь3,к+1 - К,3,к-1/2Фь3,к-1
3 '
= ^ 2Ак
_кп фг,],к+1 — фг,],к + кп фг,3,к — фг,],к-1
Х - • 1 , -. / о 4 4 ХА
Дгк+Д ДгкДк
(¿4 ф)г,3,к =
2Дк
Используя на каждом дробном шаге [Ьп, Ьп+1 ] схему Кранка-Никол-сона, алгоритм расщепления примет вид:
фп+т/3 — фп+(т-1)/3 ^ фп+т/3 + фп+(т-1)/3
ДЬ + ¿т 2 =0 т = 1 2'
Для повышения точности расчетов используется двуциклическая перестановка этапов расщепления.
На втором этапе для фильтрации звуковых волн использовалась неявная разностная аппроксимация первого порядка точности по времени - схема "естественного фильтра" (Марчук, Дымников, Залесный и др., 1984). Эволюционное уравнение для давления решалось методом покомпонентного расщепления по координатам. При реализации алгоритма использовалась немонотонная прогонка (Самарский, Николаев, 1978). Построенные конечно-разностные схемы абсолютно устойчивы, имеют первый порядок аппроксимации по времени и второй — по координатам.
4. Численные эксперименты
Для иллюстрации возможностей модели приведем результаты расчетов течений и переноса примесей в оз. Байкал.
Рисунок 1 характеризует общую циркуляцию вод в верхних слоях озера Байкал в летне-осенний период.
Вследствие действия силы Кориолиса циркуляция имеет циклонический характер, как и во всех крупных и глубоких озерах Северного полушария. При этом образуются вихри в каждой котловине в пределах единого циклонического круговорота. Течения в центральных областях характеризуются пониженными скоростями. Приведенная схема течений образуется постепенно в течение всего летне-осеннего периода. В конце безледного периода сформировавшиеся в поверхностном слое преобладающие течения проникают на более значительную глубину. Полученные схемы поверхностных течений согласуются с экспериментальными данными (Айнбунд, 1988; Shimaraev е! а1., 1994).
На рисунке 2 представлены результаты расчетов преобладающих течений в котловине Южного Байкала. Найденные на основе гидротермодинамической модели скорости движения и турбулентные характеристики используются для расчета переноса примесей от Байкальского целлюлозно-бумажного комбината (БЦБК) (рис. 3).
Рис. 1. Общая циркуляция вод озера Байкал в летне-осенний период
Номер изолинии на рисунке 3, умноженный на 10, означает концентрацию примеси в процентах от сбрасываемой БЦБК. Из этих рисунков видно, что особенности циркуляции вод Южного Байкала препятствуют распространению примесей в другие котловины, уменьшая тем самым способности Южной котловины к самоочищению и приводя к необратимым опасным последствиям. Учитывая уникальность Байкала, необходимо или ликвидировать комбинат, или перевести его на замкнутый цикл водоснабжения.
Рис. 3. Распределение сбросов БЦБК при преобладающем течении
в Южном Байкале
В. К. АРГУЧИНЦЕВ, А. В. АРГУЧИНЦЕВА 5. Заключение
Предложенная модель для стратифицированных водоемов, с одной стороны, позволяет описать крупномасштабные процессы (например, действие силы Кориолиса) в озерах, имеющих большие размеры, а с другой, - мезомасштабные явления, например, такие как термический бар, формирующий вертикальный водообмен весной и осенью вследствие различной стратификации воды в прибрежных и центральных частях водоема. Использование полного уравнения состояния с учетом минерализации позволяет моделировать плотностные течения в окрестности геотермальных источников.
Учет сжимаемости позволяет с физической точки зрения описать адекватно разнообразные локальные процессы в стратифицированных водоемах, а с математической — вследствие эволюционного типа уравнений построить эффективную конечно-разностную схему для их решения.
Модель может быть применена и для решения различных экологических задач, требующих знания детальных гидрологических характеристик.
Список литературы
1. Айнбунд М.М. Течения и внутренний водообмен в озере Байкал / М.М.Айнбунд. - Л.: Гидрометеоиздат, 1988. - 247 с.
2. Астраханцев Г.П. Дискретная гидротермодинамическая модель климатической циркуляции глубокого озера / Г.П.Астраханцев, Л.А.Руховец // Вычислительные процессы и системы. - М.: Наука, 1986. - Вып. 4. - С. 135-178.
3. Гурина А.М. Нелинейная диагностическая модель течений глубокого озера (на примере Ладожского озера) / А.М.Гурина, Ю.Л.Демин, Н.Н.Филатов // Моделирование переноса вещества и энергии в природных системах. - Новосибирск: Наука, 1984. - С. 77-89.
4. Квон В.И. Температурно-стратифицированное течение в проточном водоеме / В.И.Квон // Метеорология и гидрология.- 1979.- №6. - С. 74-79.
5. Марчук Г.И. Численное решение задач динамики атмосферы и океана / Г.И.Марчук. - Л.: Гидрометеоиздат, 1974. - 303 с.
6. Математическое моделирование общей циркуляции атмосферы и океана / Г.И.Марчук, В.П.Дымников, В.Б.Залесный и др. - Л.: Гидрометеоиздат, 1984. -320 с.
7. Марчук Г.И. Численное моделирование динамики вод озера Байкал / Г.И.Марчук, В.П.Кочергин, Е.А.Цветова // Математическое моделирование качества воды водоемов. - М.: Наука, 1978. - С. 43-51.
8. Обухов А.М. Турбулентность в температурно-неоднородной атмосфере / А.М.Обухов //Тр. Ин-та теоретической геофизики АН СССР. - 1946. - Т.1. - С. 95-115.
9. Обухов А.М. Турбулентность и динамика атмосферы / А.М.Обухов. - Л.: Гидрометеоиздат, 1988. - 413 с.
10. Пененко В.В. Структура комплекса моделей для исследования взаимодействий в системе "озеро Байкал-атмосфера региона" / В.В.Пененко, Е.А.Цветова // Оптика атмосферы и океана.-1998.-Т.11, №6.-С. 586-593.
11. Самарский А.А. Методы решения сеточных уравнений / А.А.Самарский, Е.С.Николаев. - М.: Наука, 1978. - 592 с.
12. Фукс Н.А. Механика аэрозолей / Н.А.Фукс. - М.: АН СССР, 1955. - 351 с.
13. Цветова Е.А. Нестационарные ветровые течения в озере Байкал / Е.А.Цветова // Численные методы расчета океанических течений. - Новосибирск: ВЦ СО АН СССР, 1974. - С. 115-128.
14. Цветова Е.А. Математическое моделирование циркуляций вод озера / Е.А.Цветова // Течения в Байкале. - Новосибирск: Наука, 1977. - С. 63-81.
15. Brugge R. Non Hydrostatic modelling for studies of open-ocean deep convec-tion / R.Brugge, H.L.Jones, J.C.Marshall // Deep Convection and Deep Water For-mation in the Oceans: edited by P.C. Chu and J.C. Gascard, Elsevier publishing. - 1991. -P.325-340.
16. Chen C.T. Precise thermodynamic properties for natural waters covering only the limnological range / C.T.Chen, F.J.Millero // Limnol. Oceanogr. - 1986. - V. 31, No.3. - P. 657-662.
17. Physical limnology of lake Baikal: a review / M.N.Shimaraev, V.I.Verbolov, N.G.Granin et. al. - Irkutsk-Okayama: Baikal International Center for Ecological Research, 1994. - 81 p.
18. Smagorinsky J. General circulation experiments with the primitive equations: 1. The basic experiment / J.Smagorinsky // Mon. Weather Rev. - 1963. - V. 91, No.2.
- P. 99-164.
19. Walker S.J. A three-dimensional numerical model of deep water renewal in temperate lake / S.J.Walker // Abstr. for the degree of doctor of philosophy.- Tulan University. - 1994. - 119 p.
20. Walker S.J. A three-dimensional numerical model of deep ventilation in tem-perate lakes /S.J.Walker, R.G.Watts // J. of Geophysical Research. - 1995. - V.100, No.C11.
- P. 22711-22731.
V. K. Arguchintsev, A. V. Arguchintseva
Numeric modeling of flows and admixture transfer in stratified lakes
Abstract. The given scientific work considers a nonhydrostatic nonstationary three-dimensional nonlinear model for the description of flows, temperature fields and density fields in lakes and reservoirs. The model takes into account inhomogeneities of a bottom relief and the state equation for water connecting pressure, temperature and density. The system of hydrothermodynamic equations is solved numerically with application of the method of fictitious areas. The implicit finite difference scheme of component splitting is used. The calculations results of prevailing flows in the top layers of Lake Baikal in the period without ice are presented. The motion velocities found on the basic of the hydrothermodynamic model and turbulent characteristics are used for the calculation of admixture transfer from the Baikal pulp-and-paper combine at prevailing flow in hollow of Southern Baikal.
