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

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

Том 8, № 3, 2003

ОБ ОДНОРОДНЫХ ЧИСЛЕННЫХ АЛГОРИТМАХ ПРИ СОПРЯЖЕНИИ ПОЛНОЙ И ГИДРОСТАТИЧЕСКОЙ МОДЕЛЕЙ ГИДРОТЕРМИКИ ВОДОЕМОВ*

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

The problems on the conjunction of full and hydrostatic models for describing viscous incompressible stratified fluid flows in deep water reservoirs are investigated. The approach to the construction of two-dimensional conjugate model within the projective numerical algorithm is suggested.

При численном моделировании течений вязкой несжимаемой стратифицированной жидкости, как правило, хорошо работает гидростатическое приближение уравнений Навье — Стокса (т. е. распределение давления в области течения подчиняется гидростатическому закону). Однако считается, что такая модель не вполне подходит для описания процессов, протекающих в глубоких водоемах в области подъема и опускания вод, когда вертикальная неоднородность процессов становится существенной. В этом случае следует использовать полную модель. Но поскольку гидростатическая модель более экономична в вычислительном отношении, для эффективности расчетов в таких водоемах желательно применять сопряженные модели, а именно сопряжение гидростатической и негидростатической моделей с учетом особенностей течений в разных областях водоема. Сложность заключается в том, что численные реализации этих моделей существенно различны, в особенности при использовании условия твердой крышки на свободной поверхности водоема. Условие несжимаемости жидкости сильно усложняет численную реализацию полной модели. Для решения этой проблемы стали применять проекционные методы.

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

Запишем двумерную систему уравнений Навье — Стокса, рассматриваемую в приближении Буссинеска [1], в безразмерных переменных, учитывая возможную анизотропию турбулентной вязкости:

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

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

Q L

x = Lx, z = Hz, U =—, T = —, t = Tt,

H U

e = h, u = UZ, w = Qw = Uewt, u = —ü , Ф = Qz,

L

L

H

Р = ро + Арр, Р = родН(р - г), Щх = ^ = .

Здесь х — горизонтальная координата (скорость и); г — вертикальная координата (скорость Q — расход; р — плотность, р0 — некоторая характерная плотность, Ар — предполагаемое изменение плотности; Р — давление; д — ускорение свободного падения; ш — завихренность; Ф — функция тока; — коэффициенты горизонтальной и вертикаль-

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

2 1

Р = Р - р(11-

^ р

х

Тильдой отмечены безразмерные переменные, но в дальнейшем ее будем опускать. В итоге получаем следующую систему уравнений:

ut + uux + wuz = div AteVu — - jv Pxdz,

p J

wt + uwx + wwz = div AteVw —

ux + wz = 0.

Pz

Fr2e2'

(1)

Здесь Fr = U / \[gH; Fp = U/^/gHAp/p o; Ate — диагональная матрица с элементами vtx/Rex, vtz/Reze2 на главной диагонали, Rexi = LU/vt0Xi• Пусть

V = (u,w), f

Fp2 p

. _ du 2 dw

Pxdz, 0 , rot2eV = ---£ —, A£

dz dx

d2

+ £

d 2

dz2 dx2'

Тогда Ф-проекционный алгоритм Чорина — Тимухина (проекция на подпространство со-леноидальных векторных полей осуществляется с использованием функции тока) для (1) имеет вид [2]

^ — * — n —

v —v + (vnV)v* = div Ate Vv* + fn,

u = rot2ev*, АеФ = u,

un+1 = Фz, wn+1 = —Фx.

1

z

1

1

2

z

Давление при необходимости можно восстановить из уравнений

др

п+1

дх

дрп+1 дг

^ (и* - ип+1),

¿Г2 (т* - ^п+1).

Здесь т — шаг по времени; п — номер временного слоя. Переход к гидростатическому приближению в системе (1) соответствует исключению малых членов при условиях е ^ 1, Ие^е2 ~ 1, Кех ~ 1, что приводит к следующей системе:

щ + иих + шиг = ¿¡у Аг У и — -Рх2 + /1,

гГ

дР=о,

дг ' их + = 0,

где А4 — диагональная матрица с элементами и ^ на главной диагонали. В этой системе из первого уравнения определяется компонента скорости и, из второго уравнения — давление р, из третьего — компонента скорости т. Проводя аналогичную процедуру с системой (2), придем к следующему алгоритму:

и* — ип

+ (ьпУ)и* = ¿¡у Аг У и* + /

и

п+1

и

фп+1

и

п+1

(3)

и

п+1 = фп+1

т

п+1 = _фга+1

Таким образом, для численного расчета гидростатического приближения в рамках проекционного метода (2) достаточно исключить расчет компоненты т* ив операторах го12е и Д£ положить е = 0. Это легко осуществить в процессе расчетов по данному алгоритму, что, в частности, позволяет на одних участках расчетной области реализовать полную модель, а на других — гидростатическую, т. е. осуществить сопряжение двух моделей в процессе выполнения расчетов. Тем самым в рамках одного численного алгоритма можно реализовать три модели: полную, гидростатическую и сопряженную. Исключение из расчета уравнения для вертикальной компонеты скорости т* дает заметную экономию.

Можно также рассмотреть р-проекционный алгоритм [3], который имеет вид

V* - V™

+ (мпУ)м* = ¿¡у А4£ Ум* + /

¡д2 рп+1 + 5д 2рп+1 дх2

дг2

е2гг2

Нр-^У V*,

(4)

и

п+1

и ^ 2 Рх

Гг2

т

п+1 = _ * 5т р = т — е2Гг2 р •

Здесь 5 — вспомогательный параметр, пока равный единице. После соответствующего предельного перехода V* будем искать из первого уравнения системы (4). Для определения

рп+1 и мга+1 решаем уравнения

е

рп+ = ^т ¿¡у V*,

ип+1 = и* — рх, (5)

гг2

тп+1

Легко видеть, что ¿¡у = 0 в силу первого уравнения системы (5). Уравнение для

давления содержит здесь обе компоненты скорости V*, которые приходится рассчитывать по полным конвективно-диффузионным уравнениям. Таким образом, в рамках р-проек-ционного численного алгоритма переход к гидростатическому приближению можно осуществить, если положить 5 = 0 в системе (4). Однако этот подход не дает существенных упрощений и экономии ресурсов, поэтому его реализация не представляется целесообразной.

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

Наибольшие различия между полной и гидростатической моделями, по-видимому, должны наблюдаться при наличии достаточно мощных вертикальных потоков. Для учета особенностей геометрической и сеточной моделей будем использовать два параметра: отношение глубины области к ее протяженности и отношение вертикального и горизонтального шагов сетки. При расчетах течений в особо крупных водоемах (как в случае озера Байкал) мы, будучи ограничены в вычислительных ресурсах, вынуждены прибегать к довольно грубой сетке, особенно в горизонтальном направлении, так как протяженность водоема измеряется многими десятками километров. В связи с этим мы провели исследования на нескольких модельных областях, чтобы определить зависимость результатов расчетов от соотношения ес = Нх/кх. Во всех расчетах область была прямоугольной, сетка по горизонтали и по вертикали строилась равномерная, но с разным соотношением шагов кх и кх. На вертикальных границах области выделены зоны втекания (в верхней части левой границы) и вытекания (в нижней части правой границы) с заданными равными расходами воды. Мы не будем останавливаться на постановке начально-краевой задачи и деталях численной реализации, которые достаточно стандартны и подробно описаны в [2, 4]. Отметим только, что модель дополнялась уравнениями переноса тепла и уравнением состояния Чена — Миллеро р(р, Т) [5]. Численные эксперименты были выполнены:

1) для области с соотношением горизонтального и вертикального размеров Н/Ь = 1/4, Н = 50 ми равными шагами сетки кх = кх = Н/50;

2) на той же области, но при соотношении шагов кх = Н/50, кх = Юк^;

3) на области с соотношением размеров Н/Ь = 2, Н = 100 м равными шагами по вертикали и горизонтали кх = кх = Н/100, выходное отверстие расположено на глубине Н = 60 м;

4) на той же области, но при соотношении шагов кх = Н/50, кх = 5кг.

Температура втекающей воды была равна 6 °С, скорость — 0.5 м/с, ух = = 10-4 м2/с,

поток тепла через свободную поверхность —100 Вт/м2 (охлаждение поверхности). Продолжительность расчетного периода составляла 5 ч. Во всех четырех экспериментах строилось начальное поле данных по полной модели (расчетный период 1 ч, начальное состояние —

т*

покой и полная гомотермия с температурой 7 °С). Далее расчеты выполнялись по разным моделям.

Результаты расчетов (рис. 1-5) подтвердили, что различия между полной и гидростатической моделями проявляются тем сильнее, чем больше отношение глубины водоема к его протяженности: на рис. 1 картины течений очень близки, а на рис. 3 существенно различаются в зоне вертикальной струи у правой стенки. Видно также, что полная модель дает более широкую струю (см. рис. 3). Здесь как раз и может проявиться влияние параметра ес через горизонтальную схемную вязкость. Действительно, в эксперименте 1 (при одинаковых шагах сетки) мы наблюдаем небольшие отличия у правой стенки, а в эксперименте 2 (см. рис. 2) схемная диссипация практически стирает отличия между моделями и, кроме того, уничтожает вихри ниже основной струи. В экспериментах 3 и 4 (е = 4) различия в расчетах по разным моделям при ес = 1 видны отчетливо (см. рис. 3), а при ес = 0.2 они слабые (см. рис. 4).

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

Н, м

20 40 60 80 100 120 140 160 180 200 L, м

Н, м

sn --т--==-

20 40 60 80 100 120 140 160 180 200 Ь, м

Рис. 1. Картины течений, полученные в эксперименте 1: а — гидростатическая модель, Фтш = -0.37, Фтах = 2.50; б —полная модель, ФтЫ = -0.39, Фтах = 2.51; -0.39, -0.18, -0.02, 0, 0.13, 1.25, 2.5 — линии 1-7 соответственно.

Я,м

20 40 60 80 100 120 140 160 180 200 L, м

Н, м

20 40 60 80 100 120 140 160 180 200 L, м

Рис. 2. Картины течений, полученные в эксперименте 2: a — гидростатическая модель, б — полная модель; Фтш = 0; Фтах = 2.5; 0, 0.63, 1.25, 1.88, 2.50, 3.06 — линии 1-6 соответственно.

Н, м я, м

10 20 30 40 50 Ь,м 10 20 30 40 50L, м

Рис. 3. Картины течений, полученные в эксперименте 3: a — гидростатическая модель, Фт1п = -1.18, Фтах = 3.0; б — полная модель, Фт;п = -0.78, Фтах = 3.0; -0.06, 0, 0.15, 1.5, 3.0 — линии 1-5 соответственно.

Я, м а Я, м

10 20 30 40 50L, м ю 20 30 40 50L, м

Рис. 4. Картины течений, полученные в эксперименте 4: a — гидростатическая модель, Фт1п = -0.43, Фтах = 3.0; б — полная модель, Фт;п = -0.56, Фтах = 3.0; -0.05, 0, 0.15, 1.5, 3.0 — линии 1 - 5 соответственно.

СО РАН [6-8]. В районе мыса Лиственичный глубина поперечного сечения озера достигает 900 м. Уклон дна здесь достаточно крутой, а расчеты, выполненные нами ранее [2, 4], показали, что возникающие присклоновые течения могут достигать глубинных слоев. В

н,л

Рис. 5. Картины течений, полученные по трем моделям: а — гидростатическая модель, Фтш = -0.25, Фтах = 4.0; б — полная модель, Фтт = -0.26, Фтах = 6.7; в — сопряженная модель, Фтт = -0.38, Фтах = 6.8; -0.1, 0, 0.05, 2.0, 4.0, 6.0 — линии 1-6 соответственно.

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

Выполнены два численных эксперимента: при отсутствии и при наличии притока. Относительное отклонение давления от гидростатического, рассчитанное с использованием уравнения для компоненты скорости w [4], лежит в диапазоне от 1.1 • 10-9 до 2 • 10-7, H/L = 0.1. В первом случае картины течений, полученные по трем моделям, качественно не различались, но имели некоторые количественные отличия. Отметим, что этот эксперимент выполнен с учетом сжимаемости воды, т.е. использовалось уравнение состояния Чена — Миллеро р = p(p, T), и в условиях, приближенных к натурным. Во втором случае расчет проведен в модельных условиях: был подключен приток со скоростью втекающей воды uin = 0.02 м/с и температурой Tin = 6 °C. Кроме того, в уравнении состояния принято p = 0, чтобы максимально приблизиться к условиям расчетов в [2]. На рис. 5 приведены картины течений на пятые сутки от начала расчетного периода, полученные по трем моделям. Здесь уже видны качественные отличия между тремя расчетами, которые требуют дополнительного анализа. На первый взгляд гидростатическая модель в некоторой степени препятствует развитию вихрей и их продвижению в горизонтальном направлении. Это можно объяснить тем, что в этой модели отсутствует часть горизонтальной турбулентной диффузии, которая отвечает за горизонтальное растяжение вихрей. Видно, что сопряженная модель (полагаем е = 0 при x > 4.7 км) дает в некотором смысле промежуточные результаты между полной и гидростатической моделями, как и ожидалось при ее построении.

Еще один эксперимент с использованием сопряженной модели в условиях, близких к натурным, проведен при расчетах течений в Байкале в районе дельты реки Селенги. Этот район интересен по нескольким причинам: во-первых, глубины здесь достигают 300 м, но при этом имеется мелководная шельфовая зона со средним уклоном дна порядка 0.01; во-вторых, приток создает дополнительное воздействие на водообмен. В качестве расчетной обасти использовалось поперечное сечение озера Харауз — Красный Яр. В качестве граничных условий заданы скорость течения в притоке uin = 0.05 м/с и температура втекающей воды Tin = 10 °C. Расчеты по полной модели показали, что относительное отклонение давления от гидростатического лежит в диапазоне от 1.2• 10-11 до 1.8• 10-9, H/L = 0.01. Результаты расчетов, выполненных в этих условиях по полной и гидростатической моделям, оказались практически идентичными. Хотя глубины позволяют развить достаточные вертикальные ускорения, пологий рельеф их гасит, а в результате существенных отклонений от гидростатики в задаче не возникает.

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

[1] Плсконов В.М., ПОЛЕЖАЕВ В.И., Чудов Л.А. Численное моделирование процессов тепло- и массообмена. М.: Наука, Гл. ред. физико-мат. лит., 1984. 288 с.

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

[3] БЕЛОЦЕРКОВСКИй О.М. Численное моделирование в механике сплошных сред. М.: Наука, 1984. 520 с.

[4] Бочаров О.Б., Овчинникова Т.Э. О термогравитационной конвекции в прибрежной зоне глубокого озера в период весенне-летнего прогревания // Вычисл. технологии. 1998. Т. 3, № 4. С. 3-11.

[5] CHEN C.N., Millero F.J. Precise thermodynamic properties for natural waters covering only the lymnologics range // Limnol. Oceanogr. 1986. Vol. 31(3). P. 1207-1230.

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

[7] Течения в Байкале. Новосибирск: Наука, 1977. 160 с.

[8] Шимараев М.Н. Элементы теплового режима озера Байкал. Новосибирск: Наука, 1977. 149 с.

Поступила в редакцию 4 ноября 2002 г., в переработанном виде —13 февраля 2003 г.