CC BY
32
Вычислительные технологии
Том 7, № 5, 2002
РАСЧЕТ ПАРАМЕТРОВ СЛОЯ ПЕРЕМЕШИВАНИЯ В ВОДОЕМЕ С ИСПОЛЬЗОВАНИЕМ РАЗЛИЧНЫХ СПОСОБОВ АППРОКСИМАЦИИ АДВЕКЦИИ В ВИХРЕРАЗРЕШАЮЩЕЙ МОДЕЛИ *
В. А. Шлычков Институт водных и экологических проблем СО РАН Новосибирск, Россия e-mail: slav@ad-sbras.nsc.ru
The results of numerical experiments with eddy-resolving models are presented. The monotonia and centered-difference advection schemes are used. The equivalence of the schemes is shown with the condition of optimal spectral characteristics of model. Numerical calculations and full-scale data are compared.
Введение
При численном исследовании геофизических пограничных слоев широко используется LES-подход (Large Eddy Simulation) [1], согласно которому в спектре турбулентных движений выделяется крупновихревая часть, подлежащая явному описанию на основе уравнений гидротермодинамики. Мелкомасштабная часть спектра считается подсеточной и параметризуется с привлечением различных гипотез турбулентного замыкания. Крупные вихри в природных водоемах формируются, например, в условиях плотностной неустойчивости, когда происходит выхолаживание свободной поверхности и в верхних слоях развиваются процессы термической конвекции. Наблюдения показывают [2], что конвективный обмен реализуется в виде неупорядоченного ансамбля вертикальных струй и термиков, переносящих охлажденные массы воды в глубину.
Формирование собственных пространственно-временных масштабов конвективных крупных вихрей в значительной мере обусловлено механизмами нелинейной адвекции. Поэтому качество численного решения будет зависеть от способа дискретной аппроксимации инерционных членов. Так, при анализе проникающей конвекции важна энергетика процесса [3], и с этой точки зрения выбор пространственной аппроксимации следует сделать в пользу схем со свойствами энергетической сбалансированности [4]. К использованию другого класса схем приводит требование транспортивности — совпадения направления переноса с ориентацией вектора скорости (монотонные схемы с разностями против потока)
* Работа поддержана грантом №00-05-98542 "Ведущие научные школы".
© В. А. Шлычков, 2002.
[5]. В настоящее время успешно применяются ТУБ-схемы с ограниченной полной вариацией, которые сочетают качества консервативных и монотонных схем [6]. Сопоставительный анализ адвективных схем для атмосферной ЬЕБ-модели проводился в [7].
В настоящей работе изучается влияние способов сеточной аппроксимации слагаемых адвекции на численное описание стохастического ансамбля конвективных крупных вихрей в водоеме.
В качестве исходной примем систему уравнений гидротермодинамики турбулентной жидкости, записанную в приближении Буссинеска [8]. Математическая постановка задачи базируется на выделении среднего крупномасштабного потока и крупновихревого компонента поля турбулентности с индивидуальным воспроизведением когерентных структур [9]. Размеры водоема предполагаются достаточно большими, а область решения представляет выделенную в воде колонну квадратного сечения со стороной Ь, причем среднее течение внутри колонны считается горизонтально однородным. Уравнения гидротермодинамики для компонентов скорости и,У и температуры Т среднего течения имеют вид
ди „г д ди д дУ 1ТТ д дУ д "ттг = -1У + 7Т К— - — Ш, — = -и +— К— - — Ш,
дЬ дх дх дх дЬ дх дх дх
дТ = ±Кт дТ - (1)
дЬ дх т дх дх ' где х — вертикальная координата; Ь — время; I — параметр Кориолиса; Кт — коэффициент вертикального турбулентного обмена подсеточного масштаба; черта сверху означает горизонтальное осреднение; и,у,'Ш,Т' — компоненты скорости и отклонение температуры крупновихревого течения. Уравнения для и,у,'Ш,Т' получим вычитанием (1) из исходной системы:
'и ди 1 др , ^ д тлди д_
— + т— = —— + IV + ПХу и + — К— + — ит, оЬ дх р дх дх дх дх
¿V дУ 1 др , ^ д т^дv д_
-г- + т— = —---1и + DxyV + — К— + — vw,
оЬ дх р ду дх дх дх
'т 1 др „ . ^ д тдт . ,
л = - + 'лТ' + ^т + д:ках- (2)
+ «Тт = Dxy Т' + ±К, £ + Рг,
оЬ дх дх дх дх
ди дv дт ^ дх ду дх '
' д д д д д д д д где — = — + (и + и)— + (У + V) — + т—; Dxy = тт"Кх— + — Ку —; вт — коэффициент
'Ь дЬ дх ду дх дх дх ду ду
термического расширения; д — ускорение силы тяжести; р — давление; р — плотность
воды.
Направим ось х вверх и совместим уровень х = 0 с поверхностью водоема. Для уравнений (1) поставим следующие граничные условия:
ди дУ „ дТ
— = — = 0, рерКт дх = -Во при х = 0, (3)
дТ
и = У = 0, — = 7н при х = Н, (4)
дх
где B0 — заданный поверхностный поток тепла; H — нижняя граница расчетной области; Yh — устойчивая температурная стратификация водоема на глубине. Краевые условия для системы (2) запишем в виде
du dv
— = — = 0, w = 0, T' = T'o(t, x, y) при z = 0, (5)
dz dz
f - = 0 при z = H, (6)
где ф = (u,v,w,T'), T'o — случайные малые возмущения температуры; C — скорость внутренних гравитационных волн. В качестве краевых условий по горизонтали примем предположение о периодичности процессов. В начальный момент зададим состояние покоя и линейное убывание температуры с глубиной.
Задача (1) - (6) решалась методом конечных разностей. Для интегрирования по времени разработан полунеявный метод расщепления в варианте [9], близкий к схеме "предиктор — корректор". Метод обеспечивает второй порядок точности для операторов адвективного переноса и строго неявную реализацию слагаемых, описывающих распространение быстрых гравитационных волн.
Пространственная дискретизация полей проводится в целые и полуцелые узлы с проекцией компонентов скорости на смежные грани элементарного объема в соответствии с уравнением неразрывности в (2). Рассмотрим два варианта конечно-разностных аналогов исходной задачи. Вариант A второго порядка точности по пространству получен из консервативной формы уравнений с использованием центральных разностей при аппроксимации адвективных членов. В варианте B для этой цели применяются односторонние разности с общим понижением порядка точности схемы.
Проблема корректного разделения спектра турбулентности на крупные вихри и мелкомасштабные пульсации связана со способом описания подсеточных процессов, т. е. применяемой моделью турбулентной вязкости. Аппроксимация B приводит к возникновению дополнительной "схемной" вязкости [5] и увеличивает уровень диссипации, обусловленный влиянием диффузионных членов в уравнениях. В отличие от турбулентной вязкости схемная вязкость не поддается прямой физической интерпретации и не связана явно с параметрами физической системы. В связи с этим целесообразно провести оценку влияния эффектов схемной вязкости на численное решение, задавая "контролируемый" механизм турбулентной диссипации.
Эксперимент 1. Рассмотрим простейшую модель турбулентного замыкания с постоянными значениями Kx = Ky = 10 см2/с, K =1 см2/с, характерными для рассматриваемых масштабов явления [10]. Для значений L = 100 ми H = —8 м принято сеточное разрешение 128 х 128 х 80 узлов, шаг по времени задан равным 15 с. Моделировалось развитие турбулентной проникающей конвекции в стратифицированном водоеме при возрастании B0 от нуля до максимального значения 300 Вт/м2, типичного для теплоотдачи в ночное время суток осенью [11]. По мере развития конвективного обмена в верхних слоях водоема формируется слой перемешивания с безразличной стратификацией, причем его толщина Hi растет со временем. Возрастает также средняя кинетическая энергия системы KcOnv = p(u2 + v2 + w2), где угловые скобки означают осреднение по области.
Кривые 1, 2 на рис. 1 иллюстрируют ход кинетической энергии решений, полученных по схемам A, B соответственно. Разница в решениях становится заметной при достижении "конечных" амплитуд крупных вихрей и усилении влиянии нелинейных эффектов.
Ясапу > эрг/м3
Рис. 1. Изменение со временем кинетической энергии конвективных движений при аппроксимации А (кривые 1, 3) и В (кривые 2, 4).
м3/с2 Бы, м3/с2
Рис. 2. Модельные энергетические спектры в экспериментах 1 (слева) и 2 (справа) при аппроксимации А (кривые 1, 3) и В (кривые 2, 4).
Отчетливо проявляется диссипативный характер монотонной схемы (кривая 2), выражающийся в более слабых темпах нарастания Ксопу по сравнению с центрально-разностной схемой (кривая 1). Вместе с тем, статистические характеристики полей и структура энергетических переходов практически не меняются, а относительная разница в конвективном потоке тепла юТ' не превышает 5%.
Анализ внутренней структуры решения проведем в терминах энергетических спектров. Согласно представлениям [12], оптимальное задание спектрального окна и разрешения сеточной модели должно обеспечивать детальное описание интервала энергии и начальной зоны инерционного интервала спектра турбулентных пульсаций. Мгновенное распределение спектральной функции по волновым числам к в логарифмическом масштабе иллюстрирует рис. 2, а. Отметим однотипность конфигурации спектров в моделях А, В (кривые 1, 2 разнесены по вертикали для удобства анализа). Энергетический пик расположен вблизи значения ^ к = 0, что соответствует оценкам масштаба плавучести в океане (Ьъ ~ 6 м) [10]. Согласно рис. 2, а, в диапазоне ^ к > 0.5 спектральные кривые убывают примерно по линейному закону; эту область спектра можно интерпретировать как инерционный интервал. Отметим, что наклон кривой 2 (составляющий значение около -9) слегка превышает наклон кривой 1 (~ -8), что соответствует большей диссипативности модели В. Транспортивность модели А реализуется за счет малости амплитуд коротких волн, вносящих основную фазовую ошибку при описании адвекции.
Проведенный анализ показывает близость схем А и В по воспроизведению ансамбля крупных вихрей при условии оптимальности спектральных характеристик модели. Последнее условие является существенным, так как отказ от оптимальности обусловливает появление значительных отличий в свойствах схем. Это можно показать путем расчета задачи с увеличенным размером области, задавая, например, Ь = 1000 м при неизменных значениях прочих параметров. Воспроизводимая область волновых чисел расположится левее указанной на рис. 2, а на значение Д ^ к = 1, так что коротковолновая часть спектра совместится с интервалом энергоснабжения, а разрешение инерционного интервала будет ограничено малым числом последних гармоник, недостаточным для оптимальности. Интегрирование задачи по схеме А привело к нереалистичному результату в виде пространственной "пилы", обусловленной вкладом энергетически активных старших компонентов Фурье. Решение, полученное по модели В, оказалось более гладким, спектральное распределение имеет максимум, правее которого подавление гармоник происходит в основном за счет фиктивной схемной диссипации. Это качество монотонных схем используется в [12] для фильтрации гармоник малых масштабов в процедуре "рационального осреднения". Как следует из приведенных результатов, в модели пограничного слоя с "явным" описанием турбулентности эффекты схемной диссипации невелики, однако становятся существенными, когда величина турбулентной вязкости не обеспечивает необходимое коротковолновое усечение спектра.
Эксперимент 2 служит для проверки полученных выводов с использованием более сложной модели турбулентного замыкания. Расчет Кх, Ку проводился по модели Смаго-ринского [4], а для описания вертикального подсеточного обмена использовалась двухпа-раметрическая модель
дЬ д т^дЬ де д т^де б т е2 Ь2
■77Т = Ктг + КЗ - е, — = К— + е^-КЗ - еъ —, К = ел —, (7)
дЬ дх дх дЬ дх дх Ь Ь е
где Ь — кинетическая энергия турбулентности, е - скорость диссипации, З = и2 + У,2 -двтТ; ai, еi — эмпирические постоянные. Краевые условия для системы (7) заданы из предположения о затухании турбулентных движений на границах области.
Изменение кинетической энергии в данной постановке показывают кривые 3, 4 на рис. 1, построенные на основе схем А, В соответственно. Диссипативные свойства монотонной схемы (кривая 4) здесь практически не проявляются — уровень энергетического насыщения примерно одинаков в обоих решениях, близки временные масштабы флуктуаций Ксопу, обусловленные внутренней динамикой ансамбля крупных вихрей. Незначительно отличаются размеры структурных элементов ансамбля, определенные по среднеэнергети-ческому волновому числу: 4.8 и 5.1 м при Ь = 8 ч. Анализ коэффициентов турбулентного обмена показывает, что их средние значения в консервативной схеме на 3-5% больше, чем в монотонной. То есть "эффективная" вязкость в методе В, обусловленная совместным действием турбулентной и схемной вязкостей, близка в среднем к расчетной турбулентной вязкости в методе А, что, по сути, делает эти схемы эквивалентными в смысле близости воспроизводимых статистических и интегральных характеристик крупных вихрей. Авторы [7] пришли к аналогичному выводу, рассматривая ЬЕБ-модели конвекции в атмосфере на основе монотонной и центрально-разностной аппроксимаций.
Спектральное распределение энергии в данном эксперименте иллюстрируют кривые 3, 4 на рис. 2, б. Прямые линии показывают теоретический наклон спектра в инерционном интервале [10]. Видно, что спектры решений А, В, незначительно различаясь в деталях,
Рис. 3. Распределение по глубине суммарного (кривая 1) и конвективного (кривая 2) потоков тепла: точки — данные измерений [11].
идентичны в целом, правильно воспроизводят интервал энергии и отражают асимптотику затухания турбулентных пульсаций в диапазоне больших волновых чисел.
Эксперимент 3. Тестирование модели (1)-(7) проводилось с помощью данных полевого эксперимента на водохранилище Веллингтон [11]. В период ночного охлаждения в верхнем слое водоема зарегистрирована конвекция. Поверхностный поток тепла и начальный профиль средней температуры задавались по данным наблюдений. Начальные распределения скорости дрейфа и полей турбулентности определялись с помощью пред-расчетного решения уравнений (1), (7) с использованием информации об измеренных характеристиках течения [11]. Период интегрирования задачи соответствовал 00-06 ч локального времени.
Модель правильно воспроизводит эволюцию слоя перемешивания, расчетное положение суточного термоклина близко к фактическому. Разница между модельной и измеренной температурами поверхности составила в среднем 0.07 0С. Варианты расчета по схемам А, В дали практически одинаковые профили средней температуры. На рис. 3 представлен расчетный профиль полного потока тепла (кривая 1) в сравнении с фактическими значениями (отдельные точки). Отметим неплохое количественное согласие данных. Кривая 2 иллюстрирует вклад конвективной составляющей рсрТ'ю.
Заключение
Численное исследование вихреразрешающих моделей, построенных на основе симметричной и односторонней аппроксимации адвективных членов, показало приемлемость обеих схем для прямого воспроизведения ансамбля турбулентных конвективных структур в водоеме. Необходимым условием инвариантности решений относительно способа аппроксимации является согласованность границ спектрального окна, сеточного разрешения и модели турбулентного замыкания, обеспечивающая корректное описание интервала энергии и каскадного переноса в инерционном интервале спектра. Схема с направленными разностями приводит к дополнительной диссипации энергии, не предусмотренной в исходных дифференциальных уравнениях. В моделях с расчетом турбулентной вязкости этот эффект компенсируется формированием более низких значений коэффициентов обмена по сравнению с консервативной схемой, так что характер убывания спектров в обеих схемах практически одинаков и соответствует универсальным колмогоровским законам. В целом представленная модель крупных вихрей в водоеме позволяет численно воспроизводить
конвективный ансамбль с устойчивыми статистическими и средними характеристиками,
удовлетворительно совпадающими с данными наблюдений.
Список литературы
[1] Stevens B., Lenschow D.H. Observations, experiments, and large eddy simulation // Bull. Am. Met. Soc. 2001. Vol. 82, No. 2. P. 283-294.
[2] Raash S., Etling D. Modelling deep ocean convection: large eddy simulation in comparison with laboratory experiments //J. Phys. Oceanogr. 1998. Vol. 28. No. 9. P. 1786-1802.
[3] Ван Мигем Ж. Энергетика атмосферы. Л.: Гидрометеоиздат, 1977. 327 с.
[4] Пененко В. В., Алоян А. Е. Модели и методы для задач охраны окружающей среды. Новосибирск: Наука, 1985. 256 с.
[5] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.
[6] Пинчуков В.И., Шу Ч.-В. Численные методы высоких порядков для задач газовой динамики. Новосибирск: Изд-во Сиб. отд-ния РАН, 2000. 231 c.
[7] Brown A. R., MacVean M. K., Mason P. J. TI The effects of numerical dissipation in large eddy simulations //J. Atmos. Sci. 2000. Vol. 57, No. 19. P. 3337-3345.
[8] Белолипецкий В. М., Костюк В. Ю., Шокин Ю. И. Математическое моделирование течений стратифицированной жидкости. Новосибирск: Наука, 1991. 175 с.
[9] Шлычков В. А., Пушистов П. Ю. Подобие структур конвективных пограничных слоев атмосферы и водоема: результаты численных экспериментов с вихреразрешаю-щими моделями // Вычисл. технологии. 2002. Т. 7, №2. С. 113-122.
[10] Озмидов Р. В. Диффузия примесей в океане. Л.: Гидрометеоиздат, 1985. 280 с.
[11] Imberger J. The diurnal mixed layer // Limnol. Oceanogr. 1985. Vol. 30, No. 4. P. 737770.
[12] Белоцерковский О. М., Опарин А. М. Численный эксперимент в турбулентности: от порядка к хаосу. М.: Наука, 2000. 223 с.
Поступила в редакцию 5 апреля 2002 г.
