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

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

Том 3, № 6, 1998

ОЦЕНКА ВЛИЯНИЯ “ХОЛОДНОЙ ПОВЕРХНОСТНОЙ ПЛЕНКИ” ПРИ МОДЕЛИРОВАНИИ КОНВЕКЦИИ В ВОДОЕМЕ *

В. А. Шлычков Институт водных и экологических проблем СО РАН Новосибирск, Россия e-mail: slav@iwep.nsk.su

The problem of direct numerical modeling of thermal boundary layer at liquid surface with intensive evaporation was formulated. The results of the problem solution are used for the calibration of parametrization formulae which relate the heat flux and the temperature gradient in a cold film. The formulae are tested on the problem of water mesoscale perturbation development under convective instability.

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

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

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

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

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

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

© В. А. Шлычков, 1998.

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

1. Моделирование холодной пленки

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

При постановке задачи будем исходить из уравнений Навье — Стокса и переноса тепла, записанных для водной среды [4]. Выделим среднее (фоновое) поле температуры Т, которое будем считать функцией времени Ь и вертикальной координаты г. При выводе уравнения переноса тепла воспользуемся подходом [5], основанным на методе расщепления с введением конвергентных слагаемых, отражающих влияние конвективных движений на средний поток

дТ д2Т д—*

аГ =к' -Щ2 - а;—т- (1)

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

Для описания возмущений используем линейный аналог уравнений гидротермодинамики в приближении Буссинеска

дй др „ д2й

дЬ дх + Д дх2 ’

% — -!-.^ф «

дТ ~ _ д2Т

+ ь— — д

дЬ дх2 ’

дй д— дх + дг 0

где х — горизонтальная координата, й, р — возмущения полей горизонтальной скорости и аналога давления, Ь — дТ/дг — параметр фоновой стратификации, Д — коэффициент кинематической вязкости, вт — коэффициент термического расширения воды, д — ускорение силы тяжести. Запись уравнений в форме двумерных в плоскости (х,г) не внесет существенных погрешностей в решение, т. к., согласно [4], начальные стадии развития микроконвекции характеризуются образованием вытянутых валов с двумерной структурой.

Перейдем к постановке краевых условий. По горизонтали будем предполагать периодичность искомой вектор-функции ф — (й,—,Т,р). Совместим уровень г — 0 со свободной поверхностью воды и обозначим через Н глубину, на которой затухают конвективные пульсации. Краевые условия по вертикали для уравнения (1) запишем в виде

Р'шсркш дг расрка дг + Ра^шка при г °, (3)

Т = Тн при г = Н, (4)

где рад, ср — плотность и теплоемкость воды, 0, Q — потенциальная температура и удельная влажность воздуха, ка — коэффициент турбулентной диффузии температуры в воздухе, ра, ср — плотность и теплоемкость воздуха, Ьш — удельная теплота испарения, Тц — заданная температура воды на нижней границе. Уравнение (3) выражает баланс тепловых потоков на поверхности водоема. Слагаемые в правой части (3), описывающие атмосферные потоки тепла и влаги, рассчитывались с помощью модели слоя постоянных потоков

[3], согласно которой

д 0 80 ка 8^ = С@иь(0Н — Т), ка = С& иь(0Н — 0*) ,

(5)

где ин, 0н, Qh — заданные значения скорости, температуры и влажности на верхней границе слоя г = Н, с© — коэффициент теплообмена воздушной массы с водой, Q* — насыщающее значение влажности вблизи водной поверхности.

Для системы (2), учитывая, что она имеет второй порядок дифференцирования по координате г, на границах поставим условия “твердой крышки”:

гш

0 при г = 0, г = Н.

(6)

Начальный профиль температуры зададим в виде линейной функции от г со значением на поверхности Т = 20 ° Си убыванием с глубиной 10 ° С/м (кривая 1 на рис. 1).

Рис. 1. Вертикальные профили средней температуры при Ь = 0, 50, 60, 70, 80, 90 с физического времени (кривые 1, 2, 3, 4, 5, 6 соответственно). Точечные линии — нормированные профили

пульсаций вертикальной скорости.

В силу линейности задачи (2), (6), ее решение будем искать в виде одной гармоники Фурье ф ~ ехр(ікх) с волновым числом к = 2п/Ь, где Ь — масштаб возмущения по горизонтали. В этом предположении из системы (2) получим уравнения для синус- и косинус-компонент искомой вектор-функции, зависящих от г и Ь. Эти уравнения совместно с (1) интегрировались численно, для аппроксимации по времени применялась схема Кранка — Николсона. Аппроксимация по вертикали проведена на сдвинутых сетках так,

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

Размер шага по времени ДЬ = 0.1 с задан исходя из критерия сходимости итераций по нелинейности; в линейной задаче физические временные масштабы допускают увеличение ДЬ на порядок. Дискретизация по г проводилась из соображений детализации процессов с характерным вертикальным размером несколько миллиметров. Принято Дг = 0.25 мм, количество узлов сетки равнялось 200. Расчеты по оптимизации параметров дискретизации показали, что во всех случаях имеется запас точности, позволяющий увеличить Дг в 2-4 раза. Для принятых значений сеточных параметров измельчение шага сетки практически не влияет на результаты. Оценка точности схемы проводилась дополнительно путем сопоставления численного решения с точным, полученным в рамках линейного приближения. Относительная погрешность в определении сд (см. ниже) не превышала 2 %.

В расчетах принималось Н = —5 см, Н =10 м, ин = 5 м/с, 0н = 20 °С, значение Qh рассчитывалось по упругости водяного пара при температуре 0н и относительной влажности, последняя взята равной 85 %. Остальные параметры заданы табличными значениями.

Амплитуда возмущений в (2) при Ь = 0 задавалась достаточно малой, чтобы исключить влияние нелинейностей в (1) на начальном этапе интегрирования. В первые 50 с физического времени температура водной поверхности падает за счет потерь тепла на испарение со скоростью примерно 0.9 °С/мин. В результате формируется профиль температуры с неустойчивой стратификацией 5' < 0 в приповерхностном слое (кривая 2 на рис. 1).

За этот период поля возмущений, еще не достигшие конечных амплитуд, перестраиваются со взаимной адаптацией и выходят на асимптотический режим линейной нормальной моды, характеризующейся экспоненциальным ростом амплитуды со временем, ф ~ ехр(сдЬ). Показатель роста сд определялся численно по формуле

= 1 {к (ь + дь)) д = 2ДЬ {к (Ь)) ’

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

Распределение сд по спектру волновых возмущений показывает, что свободная конвекция обладает свойством избирательности по отношению к горизонтальным масштабам пульсаций: возмущения длиной менее 1 см являются устойчивыми, а доминирующими оказываются возмущения с длиной волны Ь ^ 2.5 см, (кривая 1 на рис. 2, а), что соответствует результатам лабораторных экспериментов [1]. Время роста амплитуды в е раз Ье = 1 /сд составляет около 6 с.

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

Дальнейшее интегрирование проводилось с Ьь = 2.5 см. При Ь > 50 с амплитуда пуль-сационной скорости возрастает до значений порядка 2 • 10-3 м/с, которые являются предельными в нелинейной задаче. Влияние конвергентного слагаемого в (1) обусловливает

а б

0 1 2 3 4 5 6 7 8 X, см 0 5 10 15 20 £,м

Рис. 2. Распределение показателя скорости роста по спектру а) в задаче микроконвекции при £ = 50, 60, 70, 80, 90 с (кривые 1, 2, 3, 4, 5 соответственно); б) в задаче мезомасштабной конвекции при £ = 50 мин с учетом и без учета параметризации (7) (кривые 1, 2 соответственно) и при £ = 75, 100, 150 мин (кривые 3, 4, 5). Пунктирная кривая 6 получена при Ьъ = 6 м на момент

£ = 150 мин.

деформацию профиля Т(г) с образованием температурной волны, проникающей в глубинные слои жидкости (рис. 1, кривые 3-6). Слой приповерхностного перепада температуры сокращается до 2-3 мм; ниже формируется слой перемешивания. Линейный анализ профилей 3-6 показал, что при сохранении избирательного характера спектра неустойчивости длина доминирующей моды с течением времени смещается в сторону больших Ь (кривые 2-5 на рис. 2, а иллюстрируют результаты линейного анализа соответственно профилей 3-6 на рис. 1). Этим обстоятельством можно объяснить явление увеличения пространственных масштабов конвективных ячеек после стадии образования холодной пленки, получившее интерпретацию в [1] как результат взаимодействия микроконвекции с крупномаштабными внутренними волнами.

Конвективный теплообмен с глубинными массами теплой воды приводит к увеличению поверхностной температуры. Перестройка профиля Т(г) разрушает первоначально неустойчивую стратификацию и подавляет рост амплитуды возмущений. Решение становится квазистационарным, а величина д стабилизируется вблизи значения 180 Вт/м2, которое является характерным для поверхностного потока тепла при испарении в природных условиях. Перепад температуры АТ в холодной пленке составил 0.13 °С, а температурный градиент на поверхности Б = 150 °С/м.

Полученные результаты в целом хорошо согласуются с данными работ [1, 4]. Это делает правомерным применение данного подхода для оценок точности интегральных моделей, параметрически описывающих влияние холодной пленки. В качестве рабочих аппроксимаций рассмотрим два соотношения, взятых из работы [2]:

АТ = аг(дрт Сррт кт/д) 1/4д3/4,

(7)

АТ = а,2м/кт (т/рш )1/2, (8)

где кт — коэффициент молекулярной теплопроводности воды, т — касательное напряжение ветра, а1, а2 — безразмерные коэффициенты, подлежащие определению. Подбор а1, а2 проводился путем сопоставления левой и правой частей в формулах (7), (8) при подстановке туда численных значений АТ и д, полученных в результате решения задачи (1) - (6). Величина т определялась из модели слоя постоянных потоков (5), где Ц, варьировалась в различных вариантах от 2 до 15 м/с.

Из расчетов получены оптимальные значения а1 = 0.31, а2 = 7. При этом конвективные движения развиваются только при Ц, > 3 м/с. В этом диапазоне скоростей относительная ошибка формулы (7) не превышала 5%; погрешность соотношения (8) оказалась значительно больше и составила при Ц, = 15 м/с около 30%. Далее в расчетах будем применять параметризацию (7) с указанным значением а1.

2. Моделирование мезомасштабной конвекции

Для оценки влияния эффектов холодной пленки при ее параметрическом описании рассмотрим задачу о развитии мезомасштабных возмущений в водоеме в условиях конвективной неустойчивости. Нас будут интересовать возмущения с характерными горизонтальными размерами 10 м и временем развития порядка нескольких минут, наблюдаемые в верхних слоях океана [6, 4]. Генерация конвективных структур с такими масштабами возможна лишь при сравнительно большом перепаде температур между водой и атмосферой, который не может быть достигнут только лишь за счет испарения. Поэтому введем дополнительный механизм выхолаживания приводного слоя атмосферы, например, за счет адвекции, допускающий интерпретацию в виде натекания холодного воздушного потока на более теплый водоем. Зададим в (5) 0и = 15 °С, Т = 20 °С.

Уравнения для описания средней температуры и возмущений в водоеме имеют прежний вид (1), (2), где теперь под кш , д будем понимать соответственно коэффициенты турбулентной температуропроводности и турбулентной вязкости, считая жидкость турбули-зованной за счет процессов подсеточного масштаба.

Краевое условие (3) для уравнения (1) необходимо преобразовать с целью учета параметризации (7). Для этого, обозначив левую часть (3) д = радсркшдТ/дг, заменим ее выражением, полученным обращением формулы (7). В итоге краевое условие примет вид

д0 дЦ

(АТ/а 1 )4/3(двтсрршкт/д)1/3 = расРка — + раЬшка — при г = 0. (9)

При численной реализации условия (9) величина АТ представляется как разность температур на поверхности и втором по глубине уровне дискретной сетки по г.

Укажем значения основных параметров в данной задаче: А£ =15 с, Н = -12 м, кш = 10-4 м2/с, закон изменения горизонтальной диссипации задан в виде Д = 10-3Ь4/3 м2/с. Равномерная по г сетка содержала 200 уровней. Наклон профиля Т(г) в начальный момент времени задан величиной 0.1 ° С/м.

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

появления диапазона неустойчивости. В данной задаче этому моменту соответствует физическое время £ = 50 мин.

Анализ устойчивости возмущений проводился при фиксированном поле средней температуры путем численного интегрирования уравнений (2) по “собственному” времени начиная с £ = 0 до выхода на асимптотический режим установления сд для каждого Ь из рассматриваемого участка спектра. Таким образом, положение спектра неустойчивости как функции волнового числа менялось со временем в соответствии с изменением фоновой стратификации.

Кривая 1 на рис. 2, б иллюстрирует распределение показателя скорости роста сд, полученное в задаче (1)-(2), (4)-(6), (9), т. е. с параметризацией термической пленки. Видно, что размер максимально неустойчивого возмущения соответствует Ь 3 м, а волны с Ь < 1.5 м являются нейтральными или затухающими со временем. Кривая 2 на рис. 2, б построена для задачи (1) - (6) без учета параметризации на тот же момент времени, при этом величина дТ/дг в левой части (3) аппроксимировалась обычными конечными разностями. Сопоставляя кривые 1 и 2, приходим к выводу, что при схожести свойств спектральной избирательности обеих задач фильтрация поверхностных эффектов занижает скорость развития возмущений на порядок. Так, если максимум кривой 1 дает значение сд = 5.6 • 10-3 с-1 что соответствует £е = 2.9 мин, то для кривой 2 наиболее быстрорастущая мода имеет значение £е = 22 мин. При этом, если температуры поверхности в разных подходах отличаются незначительно, то значение д в модели с учетом холодной пленки превышает тепловой поток, вычисленный прямым способом (3), в 3-4 раза.

Пренебрежение эффектами высоких температурных градиентов в холодной пленке приводит к ослаблению мощности слоя неустойчивой стратификации и, следовательно, к понижению запаса потенциальной энергии неустойчивости. В результате формируются неадекватные временные и пространственные масштабы конвективных пульсаций (как следует из формы кривой 2 на рис. 2, б, доминирующее возмущение имеет длину Ь ^ 7 м). Величина теплового потока д = 220 Вт/м2, полученная в модели с параметризацией (7), представляется более реалистичной для принятых параметров метеорологического режима в приводном слое. Поэтому дальнейшие результаты и выводы будут касаться только модели с учетом холодной пленки.

При £ > 50 мин возмущения из диапазона неустойчивости начинают интенсивно развиваться, причем скорость возрастания амплитуды зависит от дисперсионного соотношения, представленного кривой 1 на рис. 2, б. Согласно приведенным выше соображениям, зададим базовый масштаб Ьъ = 3 ми рассмотрим результаты дальнейшего интегрирования по времени. Основные закономерности формирования вертикальной структуры средней температуры носят характер, подобный полученному при изучении микроконвекции (см. рис. 1). Отметим лишь отсутствие выраженного температурного фронта на нижней границе слоя перемешивания. Толшина слоя, охваченного конвекцией, растет со скоростью примерно 2.5 см/мин. Одновременно происходит качественное изменение конфигурации кривых, характеризующих спектр неустойчивости. Нелинейная перестройка термической структуры приводит к подавлению всех возмущений, горизонтальные размеры которых меньше базового значения Ьъ “несущей” волны. Это иллюстрирует кривая 3 на рис. 2, б, построенная по результатам анализа спектра собственных значений линейных нормальных мод для профиля Т на момент £ = 75 мин. Видим, что ранее неустойчивый интервал 1.5 < Ь < 3 м утерял это свойство, а доминантная частота переместилась к значению Ь ^ 8 м. С развитием конвективных движений тенденция увеличения масштаба максимально неустойчивой моды сохраняется; так, экстремум кривой 5, построенной при £ = 150 мин,

находится уже при Ь ^ 10 м.

Таким образом, задавая исходный размер возмущения Ь = Ьь, получим исходящий из этой точки спектра пучок дисперсионных кривых, описывающих характеристики спектра неустойчивости в разные моменты времени. Коротковолновое усечение масштабов с Ь < Ьь означает, в сущности, отрицательную обратную связь со свойствами избирательности в процессе взаимодействия среднего течения и пульсаций, когда развитие пульсаци-онной компоненты обусловливает реакцию среднего потока, направленную на погашение возмущений базового и более мелкого масштабов и стимулирования роста более крупных возмущений. Такая противокаскадная передача энергии через поле средней температуры сохраняется и для больших Ьь: пунктирная кривая 6 на рис. 2, б, построеннная для Ьь = 6 м, имеет коротковолновое усечение соответственно при Ь 6 м и максимум сд при Ь ^ 14 м (не зависящий от Ьь первоначальный спектр неустойчивости как и ранее представлен кривой 1, рис. 2, б).

Эти выводы справедливы для участка спектра, включающего значение Ьь = 30 м. Дальнейшее увеличение Ьь сопровождается принципиальным изменением характера дисперсионных кривых. Нелинейное взаимодействие среднего потока и возмущений с горизонтальными размерами более 30 м не приводит к подавлению коротковолнового экстремума, а формирует спектр неустойчивости, подобный изображенному кривой 1 на рис. 2, б с максимумом при Ь ^ 3 ми монотонным убыванием сд с возрастанием Ь. Значение частоты сд в точке максимума увеличивается с увеличением базовой длины волны и при Ьь > 50 м превышает величину показателя скорости роста линейной задачи (кривая 1, рис. 2, б). Отсюда можно сделать вывод, что эффекты взаимодействия длинноволновых возмущений со средним потоком приводят в действие механизм дополнительного стимулирования развития коротких волн, которые, в свою очередь, являются неустойчивыми и трансформируются со временем в возмущения больших масштабов. Эти эффекты, очевидно, могут быть описаны лишь в рамках более полной модели с учетом разномасштабности конвективных движений и детализированным описанием нелинейных взаимодействий.

Резюмируя сказанное, отметим, что начальный этап формирования средней температуры в водоеме обусловлен в основном механизмом турбулентной диффузии. По мере образования достаточно мощного слоя с неустойчивой стратификацией происходит дестабилизация части волнового спектра с характерным максимумом дисперсионного соотношения при Ь ^ 3 ми коротковолновым усечением при Ь ^ 1.5 м. Эволюция вертикальной структуры при конечных значениях амплитуды возмущений сопровождается подавлением наиболее неустойчивых коротких волн и возрастанием горизонтальных размеров конвективных ячеек. Одновременно заглубляется граница перемешанного слоя и, следовательно, происходит увеличение вертикального масштаба конвекции. Этот процесс ограничен масштабом по горизонтали Ь = 30 м, превышение которого ведет к деблокированию коротковолновых возмущений и возможной интенсификации их роста.

Отсюда следует, что при рассмотренном наборе физических параметров масштабы свободной конвекции в водоеме лежат в диапазоне 3 < Ь < 30 м. Можно предположить, что процесс формирования конвективных структур не является строго детерминированным а носит динамико-стохастический характер, обусловленный взаимодействием разномасштабных возмущений и среднего течения.

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

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

[1] Бунэ А. В., Гинзбург А. И., Полежаев В. И., Федоров К. Н. Численное и лабораторное моделирование развития конвекции в охлаждающемся с поверхности слое воды. Изв. АН СССР, Физика атм. и океана, 21, №9, 1985, 956-963.

[2] Гинзбург А. И., Федоров К.Н. Термическое состояние пограничного слоя охла-ждающейси воды при переходе от свободной конвекции к вынужденной. Там же, 14, №7, 1978, 778-785.

[3] Лыкосов В. Н. Параметризация пограничного слоя атмосферы в моделях крупномасштабной циркуляции. В “Вычисл. процессы и системы”. Наука, М., вып.10, 1993, 64-95.

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

[5] Пушистов П.Ю., Васкевич Л. А. Численное моделирование конвективного пограничного слоя атмосферы. Изв. АН СССР, Физика атм. и океана, 24, №11, 1988, 1142-1154.

[6] Федоров К.Н. О многообразии физических режимов верхнего слоя океана. Препринт №4, ОВМ АН СССР. М., 1981.

[7] Шлычков В. А., Пушистов П. Ю., Бочаров О. Б. Математическая модель для описания взаимодействия конвективных возмущений в пограничном слое атмосферы и водоеме. Тез. докл. междунар. конф. “Математические модели и методы их исследования”, Красноярск, 1997.

Поступила в редакцию 19 ноября 1997 г., в переработанном виде 19 мая 1998 г.