CC BY
47
Термогидродинамика океана
УДК 504.3+551.51
М.В. Шокуров*, С.Ю. Артамонов*, И.Н. Эзау**
¿ДО-модель турбулентного атмосферного планетарного пограничного слоя: описание и тестовые расчеты
Приведено описание вихреразрешающей модели ЬЕБЖС с динамической смешанной схемой турбулентного подсеточного замыкания. Детально рассмотрены использованные в модели конечно-разностные схемы и их свойства. Выполнены тестовые расчеты нейтрально стратифицированного атмосферного планетарного пограничного слоя над однородной шероховатой поверхностью. Результаты показывают, что модель ЬЕБМС правильно воспроизводит основные характеристики этого слоя - угол поворота ветра, коэффициент сопротивления, профили скорости. На основе тестовых расчетов показано, что данная модель может быть использована для целенаправленных численных исследований планетарного пограничного слоя.
Ключевые слова: атмосферный пограничный слой, численное моделирование, вихрераз-решающие модели.
1. Введение
Исследования атмосферного планетарного пограничного слоя (III 1С) важны для гидрофизики океана. Для моделирования циркуляции в океанах и морях необходимо знать граничные условия на поверхности (потоки импульса, тепла и влаги), для исследования ветровых волн - не только напряжения, но и всю структуру атмосферного ППС. Для климатических исследований важны данные о потоках парниковых газов и других веществ через поверхность океана, они необходимы для расчетов будущих сценариев изменения климата с помощью климатических моделей. Информация такого рода может быть получена при изучении атмосферного ППС. Таким образом, исследования этого слоя - как теоретические, экспериментальные, так и численные, -актуальны для геофизической гидродинамики.
Атмосферный ППС изменяется с течением времени и в пространстве. Его структура определяется вращением Земли, горизонтальными градиентами давления и температуры, потоками тепла и влаги на поверхности, стратификацией плотности (зависящей от стратификации виртуальной потенциальной температуры), интенсивностью турбулентности, наличием облачности, а также шероховатостью, наклоном и неоднородностью подстилающей поверхности.
Наиболее простые теории и модели построены для пограничного слоя, в котором многие из вышеперечисленных факторов не учитываются. А именно, наиболее простым является режим нейтрально стратифицированного (эк-мановского) ППС, возникающего при течении стационарного однородного потока воздуха с постоянной виртуальной потенциальной температурой над горизонтальной поверхностью с однородной шероховатостью при отсутствии
© М.В. Шокуров, С.Ю. Артамонов, И.Н. Эзау, 2013 0233-7584. Мор. гидрофиз. журн., 2013, № 1
3
потоков тепла и влаги на поверхности, горизонтальных градиентов температуры и облачности. В нейтральном ППС отсутствуют эффекты плавучести вследствие нагревания или охлаждения подстилающей поверхности.
При таких условиях ветер в верхней части атмосферы над ППС является геострофическим, его скорость направлена вдоль изобар и не зависит от высоты. Трение на поверхности нарушает геострофический баланс. Поскольку число Рейнольдса Re в атмосферном ППС велико, порядка 109, движение в нем всегда турбулентно. Турбулентность осуществляет перенос импульса по вертикали, возникает поток импульса с верхних уровней к поверхности, скорость ветра у поверхности уменьшается и появляется агеострофический поток поперек изобар, в котором выполняется баланс между силой Кориолиса, градиентом давления и турбулентными напряжениями.
В нижней части нейтрального ППС, в поверхностном нейтральном подслое, влиянием силы Кориолиса и градиента давления можно пренебречь, турбулентные напряжения не меняются с высотой, профиль скорости имеет логарифмический вид.
Над поверхностным подслоем сила Кориолиса становится важна, с увеличением высоты она приводит к повороту направления средней скорости по часовой стрелке (в Северном полушарии). Турбулентные напряжения с высотой уменьшаются.
С точки зрения теории для определения структуры нейтрального ППС необходимо знать турбулентные напряжения. Уравнения для напряжений включают, в свою очередь, третьи моменты, уравнения для третьих моментов включают четвертые моменты и т. д. Для обрыва этой бесконечной цепочки требуется вводить гипотезы замыкания, выражающие старшие моменты через младшие. Простейшим замыканием является связь турбулентных напряжений с вертикальным градиентом скорости, имеющим коэффициент пропорциональности, называемый коэффициентом турбулентной вязкости.
Предположение о постоянстве турбулентной вязкости с увеличением высоты дает наиболее простую теорию нейтрального ППС - экмановскую теорию. Однако такая теория далека от реальности - она приводит к слишком большому повороту направления ветра с высотой, равному 45°, а также в поверхностном подслое - к линейному профилю скорости вместо наблюдаемого логарифмического.
Более реалистичной является двухслойная модель, в нижнем слое которой турбулентная вязкость линейно растет с высотой, а в верхнем слое она постоянна. Существует класс моделей, в которых коэффициент турбулентной вязкости выражается через турбулентную кинетическую энергию (ТКЭ) и масштаб турбулентности. В свою очередь для определения ТКЭ необходимо решить уравнение ее переноса, в которое входят третьи моменты, для них нужны гипотезы замыкания.
Первой работой по изучению нейтрально стратифицированного ППС была работа [1], авторы которой предложили выразить интегральные свойства пограничного слоя (скорость трения на поверхности и* и угол а между направлением ветра на поверхности и геострофическим ветром) через «внешние» параметры - геострофическую скорость ветра G, параметр Кориолиса f и высоту шероховатости поверхности z0. Кроме скорости трения традиционно 4 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
используется также геострофический коэффициент сопротивления
Cg = u* / G.
Тот же результат был получен позже в [2] с использованием теории размерности. Авторы взяли за основу идею о сшивании двух решений - для внутреннего логарифмического подслоя и внешнего подслоя нейтрального
ППС в области высот z0 << z << , где оба эти решения должны быть справедливы. Реализация этой идеи дала так называемые законы сопротивления -универсальную зависимость основных интегральных параметров нейтрального ППС (полный угол поворота ветра и коэффициент сопротивления) от од-
G
ного безразмерного числа - поверхностного числа Россби Ro =-. Сама
А
идея и многие термины в этой работе были заимствованы по аналогии из общей теории турбулентного пограничного слоя над гладкой поверхностью без вращения, в которой сшивались внутренний вязкий подслой с внешним подслоем. Законы сопротивления позволяют выразить коэффициент сопротивления и угол поворота ветра аналитически через число Россби. Неопределенными при этом остаются две константы А и В, которые можно определить, например, из эксперимента, описанного в [3].
Начиная с 1970-х годов появилась еще одна возможность исследования ППС - численное моделирование. Вследствие большой величины числа Рей-нольдса (Re = 109) численное решение самих уравнений Навье - Стокса (так называемое прямое численное моделирование, DNS - Direct Numerical Simulation) невозможно даже на современных компьютерах. В настоящее время прямое численное моделирование турбулентных пограничных слоев выполняется при числах Рейнольдса порядка 2500 [4]. Однако этого все еще недостаточно для обеспечения ключевого критерия моделирования природных ППС - подобия характеристик ППС и их независимости от числа Рейнольдса. Минимальное значение числа Рейнольдса, при котором этот критерий выполняется, равно 40 000.
Вместо DNS используются модели больших вихрей, или вихреразре-шающие модели (LES - Large Eddy Simulation), в которых вихри размером, большим шага сетки, описываются явно, а действие более мелких вихрей параметризуется. Основная роль мелких вихрей заключается в «получении» кинетической энергии от крупномасштабных вихрей и превращении ее за счет вязкой диссипации в тепло. Однако в некоторых условиях, например вблизи твердой поверхности, подсеточная турбулентность играет роль также в переносе импульса. В настоящее время разработаны различные методы параметризации подсеточной турбулентности в LES-моделях атмосферного ППС. В следующем разделе будет дано их описание.
В отличие от рассмотренных выше моделей ППС с постоянной турбулентной вязкостью К, заданной турбулентной вязкостью K(z), турбулентной вязкостью, определяемой через ТКЭ, а также в отличие от универсальных законов сопротивления численное моделирование с использованием LES-мо-делей дает полную информацию о структуре ППС.
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
5
Кроме вертикальных профилей средней скорости, ТКЭ и турбулентных напряжений, можно рассчитать третьи, четвертые моменты и моменты более высокого порядка. Численное моделирование позволяет определить значимость различных слагаемых - ТКЭ и потоков импульса - в эволюционных уравнениях для вторых моментов и на основе этой информации построить замыкания, пренебрегая малыми членами и выражая старшие моменты через младшие. LES-модели позволяют определить спектральный состав турбулентности, роль разных масштабов в разных процессах - переносе импульса, ТКЭ и др. Из экспериментов и численного моделирования известно, что турбулентность вблизи твердой стенки анизотропна, перенос импульса осуществляется вихрями с характерной пространственной структурой в виде «шпилек». Численное моделирование позволяет «увидеть», идентифицировать и количественно описать такие или подобные «когерентные» структуры.
Первые численные расчеты с применением LES-моделей выполнялись для неустойчивого конвективного ППС [5]. Однако эти же модели вскоре были использованы и для моделирования нейтрального ППС. В [6] приводится сравнение четырех разных LES-моделей с разными численными схемами и подсеточными параметризациями. Было показано, что основные характеристики нейтрального ППС воспроизводятся даже при грубом разрешении и что они слабо зависят от выбора численной схемы и метода параметризации подсеточной турбулентности. В последнее время выполнены расчеты нейтрального слоя с более высоким разрешением и использованием более современных схем замыкания и численных схем более высокого порядка точности [7 - 13].
Естественным шагом было бы применить LES-модель для моделирования атмосферного ППС над морем. С этой целью в Морском гидрофизическом институте НАН Украины была инсталлирована LES-модель LESNIC (Large Eddy Simulation NERSC Improved Code) с динамическим смешанным подсе-точным замыканием, разработанная в NERSC (Nansen Environmental and Remote Sensing Center - Международный центр по дистанционному исследованию окружающей среды им. Ф. Нансена) одним из авторов.
Цель настоящей работы - описание модели и численных схем, представление результатов тестовых расчетов, а также обсуждение проблем, возникающих при моделировании. В следующем разделе приведено описание модели, в третьем разделе дается описание использованных конечно-разностных схем и их свойств, в четвертом разделе кратко описаны тестовые численные эксперименты и их результаты.
2. Модель больших вихрей LESNIC
LES-модель LESNIC [7] была создана для исследования турбулентного атмосферного ППС. Модель основана на уравнениях Навье - Стокса для движения несжимаемой жидкости в приближении Буссинеска при больших числах Рейнольдса. Подсеточные мелкомасштабные турбулентные движения параметризуются с помощью современной динамической смешанной схемы замыкания.
Техника LES постепенно становится популярным инструментом для исследования турбулентных течений в атмосфере и океане. LES-модели имеют
6 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, N 1
несколько важных преимуществ. Во-первых, они намного дешевле и проще для использования, чем натурные измерения. Во-вторых, они обеспечивают полную структуру трехмерных гидродинамических полей при полностью контролируемых внешних условиях. В-третьих, они более доступны в том смысле, что их можно использовать в условиях, в которых натурные измерения выполнить трудно. И наконец, ¿ДО-модели очень полезны для концептуального понимания атмосферных процессов. Действительно, тот или иной физический процесс может быть легко исключен из рассмотрения в модели, следовательно, роль различных процессов может быть идентифицирована.
Для определения характеристик турбулентного движения несжимаемой жидкости при очень больших числах Рейнольдса решается стандартный набор уравнений гидродинамики в приближении Буссинеска:
ды1 д I \ др ды
—'- =--(ыы )—— - а1, —L = 0, (1)
д1 дxJ ' 3 дх1 ' дх1
где х3 = (х,у,2) - оси декартовой системы координат, направленные на восток, север и вертикально вверх соответственно; ы1 = ы1 (х1,х2,х3,1 ) = (ы,V,ч>) -компоненты скорости; р - мгновенное значение давления. Подразумевается суммирование по повторяющимся индексам. Члены, содержащие молекулярную вязкость, опущены в силу их малости при очень больших числах Рей-нольдса. При описании нейтрального ППС плотность однородна, поэтому сила плавучести отсутствует.
Вертикальная и горизонтальные компоненты скорости вращения Земли учитываются в параметрах Кориолиса f = 20 ътр, /' = 20 со8(, где
0 = 7,29 10-5с-1 - угловая скорость вращения Земли, р - широта. В традиционном приближении в геофизической гидродинамике пренебрегают компонентами силы Кориолиса /м> и /и , пропорциональными локальной горизонтальной компоненте угловой скорости вращения Земли 0,у = Осо8ф, учитывают только локальную вертикальную компоненту = ^тф. Однако в результате недавних исследований было выяснено, что горизонтальная компонента может быть важна в атмосферном ППС [8, 14]. Поэтому в данной работе использовалось полное точное выражение для ускорения Кориолиса:
ах =- /V + /' м>,
ау = /и, (2)
а2 =- / Ы.
Наиболее естественным подходом к численному решению рассматриваемой системы дифференциальных уравнений является прямое численное моделирование ОЫБ, в котором все масштабы движения считаются явно. Однако, исходя из теории Колмогорова [15], минимальный размер вихрей, присутствующих в турбулентном потоке, есть величина порядка Яе-34. Таким образом, чтобы явно описать такие вихри в численной модели, шаг сетки должен быть равен к = о(Яе34) и, следовательно, число узлов сетки N = о(Яе9 4).
0233-7584. Мор. гидрофиз. журн., 2013, № 1
7
Для планетарного пограничного слоя Яе = 109, т. е. сеточная модель должна содержать примерно 1021 узлов. Современные компьютеры не достигли производительности, при которой возможно использование прямого численного моделирования.
Согласно гипотезе Колмогорова [15], при достаточно больших числах Рейнольдса в спектре трехмерной изотропной турбулентности существует инерционный интервал, в котором энергия не продуцируется и не диссипиру-ет, а только передается от больших волновых чисел к малым. Спектр энергии пульсаций в инерционном интервале задается соотношением
5 (к ) = ске21ъ к "5/3,
где ск - постоянная Колмогорова; е - скорость диссипации турбулентной энергии во всем диапазоне частот; к - волновое число.
Если вышеприведенные предположения выполняются, тогда, чтобы адекватно описывать статистические характеристики крупномасштабных вихрей, достаточно построить численную модель, верно воспроизводящую генерацию турбулентных движений в низкочастотном диапазоне и перераспределение энергии во всем или хотя бы в части инерционного интервала (рис. 1).
Считаются явно Параметризуются
-►
Р и с. 1. Качественное представление о разделении масштабов в вихреразрешающих моделях
Таким образом, в LES-моделях движения с пространственными масштабами, которые больше шага сетки, описываются явно, подсеточные же масштабы параметризуются. За время использования LES-моделей методы параметризации подсеточных процессов претерпели развитие от простой параметризации Смагоринского, использованной в [5], до наиболее современных [7, 13].
Ниже приводится краткое описание модели LESNIC. Математическая техника LES основана на свертке уравнений Навье - Стокса с некоторым
8
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
пространственным фильтром, что приводит к разделению параметров течения на крупномасштабную и мелкомасштабную составляющие. В качестве фильтра обычно выбирается гауссов фильтр
gSxj - xj ) =
6(xj ~ xj)2
J сТПл7
6
rexp
(3)
где А J = (Ах,Ау,Аг) - шаг сетки; А - ширина пропускания фильтра, определяющая наименьший размер вихря, остающегося после применения фильтра. Тогда скорость может быть представлена в виде суммы крупномасштабной (разрешенной) и подсеточной составляющих: и = и1 + и", где
u' = \u
b
(x'j )ga (Xj - x'j ]dx'j, a = Xj - A. /2, b = Xj+ A . / 2. (4)
Осредненное давление определяется аналогичным образом.
После применения процедуры пространственной фильтрации (3) к системе уравнений (1) получим:
ди' д ! , , \ др' ди' п ...
—^ =--(и'и'.+ т..)-—, —L = 0, (5)
дГ дх/'' дх/ дх, У'
т. = (и¡и.)' - (и,)' (и..)', (6)
где т. - тензор турбулентных напряжений.
Для определения давления к уравнению Навье - Стокса применяется оператор дивергенции. Учитывая уравнение неразрывности, получим уравнение Пуассона для осредненного давления:
Э2p' Э Э ( 1 1 \
*ТЭТ lu'uj + TJ). (7)
Система (5) является незамкнутой. Поэтому необходимо использовать некоторые физические предположения, чтобы получить замкнутую систему для средних составляющих турбулентных характеристик. Центральной проблемой вихреразрешающего моделирования является выбор способа параметризации турбулентных потоков t., обусловленных взаимодействием мелкомасштабных (подсеточных) и крупномасштабных (описываемых явно) вихрей. При этом необходимо выразить интересующие нас потоки в терминах известных в модели фильтрованных величин. На данный момент существует большое число подходов к решению этой проблемы. Такие параметризации называются подсеточными турбулентными замыканиями. Универсальный подход к построению замыканий для LES-моделей пока не найден. Опишем лишь динамическую смешанную схему замыкания DMM (Dynamic Mixed Subgrid Closure Mode'), которая используется в модели LESNIC.
a
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
9
Подставляя ы = ы1 + и* в (6), получим разложение тензора турбулентных напряжений тт на три части:
т.. (ы.,ы.) = Ь. + Я + С. , (8)
3У '' 3' '3 '3 '3 , у '
где Ь3 = (ы1м3)' - (ы[)' (ы3)' - тензор Леонарда, эта часть напряжений создается разрешенными крупномасштабными вихрями; = (ы*ы*) - (ы*) (ы*) - тензор Рейнольдса, турбулентные напряжения этого тензора создаются неразрешенными подсеточными вихрями; С 3 = (ы\ы*)' - (ы\)'(ы*)' + (ы*ы3)' - (ы*)' (ы3)' -
эта часть напряжений создается как разрешенными, так и отфильтрованными вихрями.
Исторически первым замыканием для турбулентных напряжений была модель вихревой вязкости Смагоринского [16]. В этом замыкании предполагается, что анизотропную часть тензора турбулентных напряжений можно представить в виде
т -=-2С^|, (9)
здесь CS - постоянная Смагоринского; S. = —
1 ( du du.
■ + -dxj dx.
тензор скоростей
деформации; |S| =^2SjS~j - норма тензора скоростей деформации. Заметим,
что для несжимаемой жидкости след тензора турбулентных напряжений равен нулю и уравнение (9) может быть записано в более простой форме:
t =-2Cs2ä2|Sl\S' . (10)
Модели вихревой вязкости достаточно верно описывают суммарную диссипацию крупных вихрей и правильно воспроизводят спектр диссипации в диапазоне волновых чисел, разрешаемых моделью. Однако корреляция компонент тензора турбулентных напряжений, рассчитанных на основе модели вихревой вязкости, с компонентами этого тензора, которые рассчитаны по результатам прямого численного моделирования, мала (0 - 0,25).
Второй класс замыканий основан на соображениях подобия, в этих моделях тензор турбулентных напряжений аппроксимируется путем формальной замены неотфильтрованной скорости на известную фильтрованную:
t=с (uu )-u ) u)), (11)
где С - эмпирический коэффициент пропорциональности, часто принимаемый равным единице. Вторая фильтрация с верхним индексом L проводится с привлечением либо того же самого фильтра, либо более широкого фильтра.
Модели, основанные на соображениях подобия, достаточно хорошо воспроизводят корреляцию между вторыми моментами, однако значительно занижают диссипацию, что может приводить к численной неустойчивости в
10 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
моделях, использующих бездиссипативные аппроксимации адвективных слагаемых. В практических расчетах модели, построенные из соображения подобия, используются в комбинации с моделями вихревой вязкости. Такие модели принято называть смешанными [17].
Третий класс замыканий основан на так называемом динамическом подходе [18]. Предполагается, что одна и та же подсеточная модель без какого-либо изменения может использоваться как при более грубом, так и при более мелком пространственном разрешении. Применяя более широкий фильтр к (8), определим
Т =((ии ) )-(и'и') . (12)
С другой стороны, тензор турбулентных напряжений Т., полученный после последовательного применения двух фильтров к системе уравнений Навье -Стокса, имеет вид
Т =((ии ) ^-(и'Щ^ . (13)
Вычитая (12) из (13), получим
Т.-Т =Ш)-(и')(и.) . (14)
Последнее тождество уже не содержит нефильтрованных переменных. Тогда, согласно предположению, используемому в динамическом подходе, можно применить одну и ту же подсеточную модель для вычисления т. и Т.. Причем в данном случае постоянная Смагоринского будет играть роль свободного параметра, вычисляемого на каждом шаге по времени и используемого для минимизации погрешности.
В динамической смешанной схеме замыкания тензор турбулентных напряжений представляется в виде
тт. =(и'и.)-(и]) и ) - 201^8% , (15)
где тензор Леонарда параметризуется с использованием замыкания, построенного из соображений подобия, а вторая часть параметризуется с использованием модели вихревой вязкости с постоянной Смагоринского, вычисляемой на каждом шаге как в динамическом методе. В данной работе применялась динамическая смешанная схема замыкания для турбулентных напряжений.
3. Численная схема
Уравнение Навье - Стокса и уравнение неразрывности были представлены в конечно-разностной форме на сдвинутой С-сетке Аракавы. Давление и температура задаются в центрах ячеек, компонента скорости и - в центре левой грани, компонента скорости V - в центре передней грани, вертикальная компонента скорости w - в центре верхней грани ячейки.
Схема для вычисления адвективных членов основана на центральных разностях [19], реализованы два варианта адвекции - в дивергентной и в ко-
0233-7584. Мор. гидрофиз. журн., 2013, № 1
11
сосимметричной форме, в настоящей работе для расчетов использовалась ко-сосимметричная форма. Схема адвекции полностью консервативна, сохраняет массу, импульс и кинетическую энергию, эти свойства консервативности необходимы для предотвращения нелинейной неустойчивости. Далее приведены конечно-разностные выражения для кососимметричной формы.
Кососимметричная форма адвективного ускорения имеет вид
Эи . 1 Э(и,и.) 1 Эи.
—'- =--— + — и.—-
2 J Эх,.
Эх 2 Эх
(16)
В конечно-разностной форме на сдвинутой С-сетке Аракавы кососимметричная форма адвективного ускорения записывается следующим образом:
1 SUU;) 1 ---—+ —
2 Sxt 2
-х, S и; —'J Sx
J J
где оператор осреднения по , -и координате
-х 1 фх' = —
2
f f Ах j^..., х, +
2
+
Ах, ' 2
а оператор первой центральной разности по , -и координате
Sj
Sх.
Ах
С ( Ах j|..., х, + 2
• V
-j
, х — -
Ах
_I
2
(17)
(18)
(19)
Такая схема имеет второй порядок точности.
К недостаткам данной схемы можно отнести ее немонотонность и, как следствие, - вычислительную дисперсию: в области резких градиентов скорости или температуры генерируется вычислительная «рябь» на самых малых масштабах. Однако поскольку модель предназначена для исследования турбулентности, имеющаяся в этой модели подсеточная турбулентная вязкость эффективно подавляет эту «рябь».
Динамическое давление вычислялось методом дробного шага [20, 21]. Идея этого метода состоит во введении «промежуточной» скорости й между моментами времени п и п + 1 и «расщеплении» уравнения для скорости на два уравнения:
и — и
At
= Adv + V's .
n+1 Л
и — и
At
= —grad( pn+l),
div^1) = 0,
(20)
(21)
(22)
здесь Adv и V's - адвективные и вязкие члены; grad и div - конечно-разностные выражения для градиента и дивергенции. Из первого уравнения определяется промежуточная скорость и, она подставляется во второе уравнение.
12 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
и
\
х
1
Взяв конечно-разностную дивергенцию от второго уравнения и воспользовавшись (22), получим уравнение Пуассона для давления, в правой части которого стоит дивергенция промежуточной скорости:
Ар- = ^ = р. (23)
А?1
Для решения этого уравнения при условии периодичности по горизонтальным координатам использовалось двойное преобразование Фурье по х, у:
к2 -к2
рп+ = р , (24)
где ~п+' и р - преобразование Фурье от рп+1 и р. Получившееся обыкновенное дифференциальное уравнение по г второго порядка после дискретизации сводилось к линейной системе с трехдиагональной матрицей, которая решалась методом прогонки. После решения трехдиагональной системы давление в физическом пространстве рп+1 получалось из ~рп+1 обратным преобразованием Фурье.
После решения уравнения для давления скорость в момент времени п + 1 определялась из (21):
ип+1 = и -А? • егаа( рп+'). (25)
Турбулентные вязкие напряжения т. рассчитываются в центрах ячеек как результат применения схемы подсеточного замыкания. Производные от турбулентных напряжений по пространственным координатам аппроксимируются первыми центральными разностями.
Для интегрирования по времени использовалась явная схема Рунге - Кут-та четвертого порядка в формулировке [22]. Четвертый порядок точности позволяет увеличить шаг по времени в 4 - 5 раз по сравнению со стандартными схемами второго порядка - «чехарда» и схемой Адамса - Бэшфорта.
Поскольку для вязких членов использовалась явная схема, шаг по времени ограничивался не только адвективным критерием Куранта - Фридрихса -Леви, но и коэффициентом турбулентной вязкости. Величина шага по времени выбиралась автоматически и рассчитывалась по формуле
(,.' . / л 1 1 У\
А? С€РЬ
111
" + —2 + ТТ
+12 8
Ах ^ " V Ах> V™ ^ — /у
Ах2 Ау2 Аг
(26)
Для реализации турбулентного замыкания использовались трехточечные фильтры Симпсона - для базового и более широкого тестового фильтра:
и' (х) = а(и. (х + Ах) + и. (х - Ах) + Ьи. (х)), (27)
где коэффициенты а и Ь равны 1/24 и 22 для основного гауссова фильтра и 1/6 и 4 - для более широкого фильтра.
Для задания нижнего граничного условия на твердой границе (подстилающей поверхности) предполагается, что с нижнего уровня модели (на вы-
0233-7584. Мор. гидрофиз. журн., 2013, № 1 13
соте Az/2), на котором заданы горизонтальные компоненты скорости, к поверхности импульс переносится неразрешенными турбулентными вихрями, а сама величина потока импульса (напряжения) определяется из предположения логарифмического профиля скорости в слое от поверхности до Az/2:
и* . (28)
ln| — I
I 2 Zo J
Тогда подсеточное турбулентное напряжение на поверхности определяется следующим образом:
= и2 и(х,*Az/2) , т = и*2 ,v(х,yAz/2) . (29)
хг | v( х, y, Az/2)| | v( х, y, Az/2)|
На верхней границе ставилось условие скольжения, то есть нулевые напряжения t = 0 , tyz = 0 . На боковых границах использовались традиционные
для ZES-моделей периодические граничные условия по х и у для всех переменных:
j(х + Ьх,y) = j(х,y), j(х,y + Ly) = j(х,y). (30)
Такие граничные условия позволяют использовать преобразование Фурье при решении уравнения Пуассона для определения давления в прямоугольной области.
4. Результаты тестовых численных экспериментов
4.1. Описание тестовых экспериментов
Согласно существующей теории нейтрального III 1С, его структура определяется только одним безразмерным числом - поверхностным числом Росс-
G
би Ro =-. Геострофическая скорость ветра G имеет типичное значение 5 -
fz0
10 м/c, в штормовых условиях она может достигать 20 м/с. В настоящей работе для тестовых экспериментов была выбрана геострофическая скорость 5 м/с из соображений экономии вычислительных ресурсов, при меньшей скорости шаги по времени для численной схемы можно сделать больше. Параметр Кориолиса во внетропических широтах также меняется незначительно, для умеренных широт (j = 50°) было выбрано типичное значение f = 1,17 10 4 c-1. Следует отметить, что для экваториальной области определение структуры нейтрального ППС представляет значительные трудности в силу его большой толщины.
Наиболее изменчивым параметром, входящим в определение числа Росс-би, является высота шероховатости поверхности. Ее значение изменяется на Земле на четыре порядка в пределах от 1 м для поверхности суши, покрытой лесом, до 10-4 м для поверхности океана при слабых ветрах. Для тестовых экспериментов выбрано значение z0 = 0,1 м. Таким образом, число Россби было равно 4,27-105.
14 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
Для тестовых численных экспериментов выбрана область в виде бокса длиной Ьх = 4 км, шириной Ьу = 2 км и высотой Ьг = 1,5 км. Число точек сетки Ых = 128, Ыу = 64, М = 128. Таким образом, шаг сетки по горизонтали Ах = Ау = 31,25 м, по вертикали Аг = 11,72 м.
Источником турбулентности в нейтрально стратифицированном ППС является сдвиг скорости. В начальный момент времени профиль скорости однородный, сдвиг равен нулю, поэтому для инициализации развития турбулентности в начальные условия к вертикальной компоненте скорости добавлялась случайная компонента небольшой амплитуды 0,1 м/с на нижних 10 уровнях, то есть от поверхности до 117 м. В противном случае на развитие турбулентности потребуется большее время.
В проведенных тестовых численных экспериментах граничные условия на боковых границах были периодическими, на верхней границе задавалось условие скольжения, на нижней границе напряжение трения определялось «пришиванием» логарифмического профиля к нижнему уровню (см. раздел 3).
В тестовых экспериментах интегрирование выполнялось на двое суток, шаг по времени рассчитывался автоматически по формуле (26), типичная величина шага составляла 0,5 с.
При выполнении расчетов начальный профиль выбран равномерным по г, равным геострофической скорости ветра: и = О = 5 м/с, V = 0. Это приводит к возникновению инерционных колебаний с периодом 2п// = 18 ч. Введем количественную меру нестационарности ППС следующим образом: проинтегрируем по вертикали от поверхности до верхней границы уравнения движения
Тогда в предположении отсутствия напряжения на верхней границе получим
В левой части этих уравнений стоит ускорение полного (интегрального по г) потока, в правой части - силы, действующие на всю толщу пограничного слоя. Первое слагаемое - сумма силы Кориолиса и градиента давления, второе слагаемое - напряжение трения, действующее на нижней поверхности. В стационарном состоянии эти два члена полностью компенсируют друг друга, поэтому определим для каждой компоненты скорости в качестве меры неста-
4.2. Инерционные колебания
(31)
(32)
0233-7584. Мор. гидрофиз. журн., 2013, № 1
15
ционарности отношение напряжения трения на нижнеи поверхности к сумме силы Кориолиса и градиента давления:
г Lz
C = - — Г(V-V (и w )0 0
0 0 (33)
г
С = —— \(и - и
(vw)o 0
Зависимость величин Си(0, С„(0 от времени приведена на рис. 2 для Яо = = 4,27-105, видны постепенно затухающие инерционные колебания. Для установившегося режима величины Си(1), С„(0 должны быть равны единице. Время затухания инерционных колебании и установления равновесного состояния оказалось достаточно большим. Хотя к концу вторых суток инерционные колебания и продолжались, далее при обсуждении результатов для осреднения по времени выбирался интервал с 24 до 48 ч.
Р и с. 2. Зависимость коэффициентов нестационарности Cu(t), Cv(t) от времени
4.3. Профили скорости, спираль Экмана, логарифмический профиль
Профили компонент скорости и напряжение на поверхности определялись осреднением в горизонтальной плоскости и осреднением по времени с 24 до 48 ч.
На рис. 3 приведены профили компонент скорости u(z), v(z), на рис. 4 та же информация представлена в параметрической форме в виде экмановской спирали. Видно, что, в отличие от классической экмановской спирали с постоянным коэффициентом вязкости и углом поворота ветра 45°, в турбулентном экмановском пограничном слое с переменным коэффициентом турбулентной вязкости экмановская спираль сильнее «прижимается» к направлению геострофического ветра. При Ro = 4,27-105 угол поворота равен 15,6°, это значение согласуется с полученными ранее экспериментальными и модельными значениями [3, 6, 7, 9 - 11].
16 ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
гШ
1,4
0,6
0,2
-0,2
0,2
0,4
0,6
0,8
1 и/в у/в
Р и с. 3. Профили компонент скорости ветра, обезразмеренных на геострофическую скорость G = 5 м/с. По вертикальной оси - безразмерная высота 2/ / и*
у/в
-0.2
0 0.2 0.4 0.6 0.8
Р и с. 4. Экмановская спираль средней скорости ветра
и/в
Скорость трения на поверхности и* рассчитывается в модели, для данного тестового численного эксперимента она равна 0,246 м/с. Геострофический коэффициент сопротивления О = и* / G = 0,049, это значение также согласуется с полученными ранее значениями [3, 6, 7, 9 - 11].
Скорость трения определяет логарифмический профиль скорости в нижнем поверхностном подслое. На рис. 5 представлены профили модуля скорости Щг) = (и2(г)+ у2(г))1/2. Видно, что логарифмический профиль очень хорошо воспроизводится моделью. Тангенс угла наклона профиля скорости вбли-
0233-7584. Мор. гидрофиз. журн., 2013, № 1
17
зи поверхности определяется скоростью трения. Прямой линией изображен
теоретический логарифмический профиль и (г) =—1п
к
( г ^
V го У
Р и с. 5. Логарифмический профиль скорости ветра в поверхностном слое. Квадратиками показаны реальные профили, прямыми линиями - логарифмические. По вертикальной оси - безразмерная высота в логарифмическом масштабе, по горизонтальной - модуль скорости ветра
5. Заключение
В данной работе с использованием модели больших вихрей ЬЕБШС с конечно-разностными аппроксимациями высокого порядка точности и современной схемой подсеточного турбулентного замыкания выполнены тестовые численные расчеты для нейтрально стратифицированного атмосферного планетарного пограничного слоя. Результаты показывают, что модель ЬЕБШС правильно воспроизводит основные характеристики нейтрального ППС -угол поворота ветра, коэффициент сопротивления, профили скорости. В процессе установления стационарного режима нейтрального ППС возникают инерционные колебания большой амплитуды, которые медленно затухают. Это требует дополнительных вычислительных ресурсов. На основе тестовых расчетов можно сделать вывод, что данная модель может быть использована для целенаправленных численных исследований ППС.
Более детальное описание структуры нейтрального ППС и ее зависимости от поверхностного числа Россби на основе серии численных экспериментов с разными значениями Яо будет дано в следующей статье.
Работа выполнена при частичной поддержке европейской программы ЕР-7 (проект РБЬ-РМЕ^).
18
!ББМ 0233-7584. Мор. гидрофиз. журн., 2013, № 1
СПИСОК ЛИТЕРАТУРЫ
1. Rossby C. G., Montgomery R.B. The layers of frictional influence in wind and ocean currents // Phys. Oceanogr. Meteorol. - 1935. - 3, № 3. - P. 1 - 101.
2. Казанский А.Б., Монин А.С. О турбулентном режиме выше приземного слоя воздуха // Изв. АН СССР. Сер. геофиз. - 1960. - 1. - С. 165 - 168.
3. Зилитинкевич С.С. Динамика пограничного слоя атмосферы. - Л.: Гидрометеоиздат, 1970. - 292 с.
4. Iwamoto K., Kasagi N., Suzuki Y. Direct numerical simulation of turbulent channel flow at Re = 2320 // Proc. 6th Symp. Smart Control of Turbulence, Tokyo, March 6 - 9. - 2005. -P. 327 - 333.
5. Deardorff W.J. Numerical investigation of neutral and unstable planetary boundary layers // J. Atmos. Sci. - 1972. - 29. - C. 91 - 115.
6. Andren A., Brown A.R., Graf J. et al. Large-eddy simulation of a neutrally stratified boundary layer: A comparison of four computer codes // Quart. J. Roy. Meteorol. Soc. - 1994. - 120. -P. 1457 - 1484.
7. Esau I. Simulation of Ekman boundary layers by large eddy model with dynamic mixed subfilter closure // Environ. Fluid Mech. - 2004. - 4. - P. 273 - 303.
8. Esau I. The Coriolis effect on coherent structures in planetary boundary layers // J. Turbul. -2003. - 4, № 1. - P. 1 - 19.
9. Zilitinkevich S., Esau I. On integral measures of the neutral barotropic planetary boundary layer // Bound.-Layer Meteorol. - 2002. - 104, № 3. - P. 371 - 379.
10. Zilitinkevich S., Esau I. The effect of baroclinicity on the depth of neutral and stable planetary boundary layers // Quart. J. Roy. Meteorol. Soc. - 2003. - 129. - P. 3339 - 3356.
11. Zilitinkevich S.S., Esau I.N. Resistance and heat transfer laws for stable and neutral planetary boundary layers: Old theory, advanced and re-evaluated // Ibid. - 2005. - 131. - P. 1863 -1892.
12. Глазунов А.В. Моделирование нейтрально стратифицированного турбулентного потока воздуха над шероховатой поверхностью // Известия РАН. ФАО. - 2006. - 42, № 3. -C. 307 - 325.
13. Глазунов А.В. Вихреразрешающее моделирование турбулентности с использованием смешанного динамического локализованного замыкания. Часть 1. Формулировка задачи, описание модели и диагностические численные тесты // Там же. - 2009. - 45, № 1. -C. 7 - 28.
14. Глазунов А.В. О влиянии направления геострофического ветра на турбулентность и квазиупорядоченные крупномасштабные структуры в пограничном слое атмосферы // Там же. - 2010. - 46, № 6. - C. 786 - 807.
15. Колмогоров А.Н. Локальная структура турбулентности в несжимаемой жидкости при очень больших числах Рейнольдса // Докл. АН СССР. - 1941. - 30, № 44. - С. 99 - 102.
16. Lilly D.K. The representation of small scale turbulence in numerical simulation experiments // Proc. IBM Sci. Computing Symp. Environmental Sci. - 1967. - 320 - 1951. - P. 195 - 210.
17. Zang Y., Street R.L., Koseff J. A dynamic mixed subgrid scale model and its applications to turbulent recirculation flows // Phys. Fluids A. - 1993. - 5, № 12. - P. 3186 - 3196.
18. Germano M., Piomelli U., Moin P. et al. A dynamic subgridscale eddy viscosity model // Phys. Fluids. - 1991. - 3, № 7. - P. 1760 - 1765.
19. Morinishi Y., Lund T.S., Vasiliev O.V. et al. Fully conservative higher order finite difference schemes for incompressible flows // J. Comput. Phys. - 1998. - 43. - P. 90 - 124.
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
19
20. Kim J., Moin P. Application of fractional step method to incompressible Navier-Stokes equations // Ibid. - 1985. - 59. - P. 308 - 323.
21. Armfield S., Street R. The fractional step method for the Navier-Stokes equations on staggered grids: The accuracy of three variations // Ibid. - 1999. - 153. - P. 660 - 665.
22. Jameson A., Schmidt W., Turkel E. Numerical simulation of the Euler equations by finite volume method using Runge-Kutta time stepping schemes // AIAA. - Paper 1981 - 1259. -P. 1 - 17.
*Морской гидрофизический институт НАН Украины, Материал поступил
Севастополь в редакцию 20.09.11
E-mail: shokurov.m@gmail.com, После доработки 08.11.11 sergei.artamonov@gmail.com
"Международный центр по дистанционному исследованию окружающей среды им. Ф. Нансена, Норвегия, г. Берген E-mail: igor.ezau@nersc.no
АНОТАЦ1Я Наведено опис вихоророздшьноТ моделi LESNIC з динашчною змшаною схемою турбулентного шдаткового замикання. Детально розглянуп використаш в моделi шнцево-рiзницевi схеми та Тх властивосп. Виконаш тестовi розрахунки нейтрально стратифжованого атмосферного планетарного прикордонного шару над однорвдною шорсткою поверхнею. Результата показують, що модель LESNIC правильно вщтворюе основш характеристики цього шару - кут повороту в^у, коефщент опору, профш швидкостГ На основi тестових розрахун-шв показано, що ця модель може бути використана для цшеспрямованих чисельних досль джень планетарного прикордонного шару.
Ключовi слова: атмосферний прикордонний шар, чисельне моделювання, вихоророздiльнi модель
ABSTRACT Description of the eddy-resolving model LESNIC with dynamic mixed scheme of turbulent sub-grid closure is given. The finite-difference schemes used in the model and their features are examined in details. Test calculations of the neutrally stratified atmospheric planetary boundary layer over a homogeneous rough surface are performed. The results show that the model LESNIC correctly reproduces main characteristics of this layer, i.e. the angle of wind rotation, the resistance coefficient and velocity profiles. Based on the test calculations it is shown that the model can be used for targeted numerical investigations of the planetary boundary layer.
Keywords: atmospheric boundary layer, numerical modeling, large eddy simulation models.
20
ISSN 0233-7584. Мор. гидрофиз. журн., 2013, № 1
