Известия РАН. Механика жидкости и газа, 2021, № 3, стр. 100-120

1.1. Турбулентность на входе

Для построения надежного метода LES с целью моделирования турбулентного атмосферного пограничного слоя требуется задание реалистичных турбулентных пульсаций на входе в область, которые будут далее эволюционировать во всей расчетной области. С теоретической точки зрения эти пульсации должны удовлетворять нескольким критериям, а именно (а) они должны быть стохастически переменными на масштабах вплоть до пространственных и временных масштабов фильтра; (б) они должны быть совместимы с уравнениями Навье–Стокса; (в) они должны состоять из когерентных вихрей во всем диапазоне пространственных масштабов вплоть до длины фильтра; (г) они должны позволять легкое определение турбулентных свойств; (д) они должны быть легко реализуемы [10].

Двумя из наиболее распространенных подходов к генерации входных турбулентных пульсаций в методе LES являются синтетический и предваряющий методы. В синтетическом методе на входе строятся случайные поля, а в предваряющем методе необходимые пульсации определяются путем дополнительного моделирования. Последние методы более точны, но требуют бóльших вычислительных затрат и труднее в исполнении [10]. Их обзор можно найти в статье [10]; далее в этой работе они не рассматриваются. Синтетический метод, выбранный в настоящей работе для генерации входной турбулентности, кратко обсуждается ниже. Более детальное описание синтетических методов можно найти в работах [9, 10].

В работе [11] развита синтетическая модель, первоначально предложенная в [12]. В ней входные турбулентные пульсации создаются изменением масштаба поля скоростей в точке, находящейся ниже по течению, и заданием этого поля в качестве граничного условия на входе. Таким образом, экономным образом пограничный слой развивается во времени и пространстве. По сравнению с простейшими методами задания случайных возмущений на входе, данный синтетический метод сокращает расстояние адаптации вниз по течению до десятикратной высоты пограничного слоя [11]. Другой распространенный синтетический метод – это вихревой метод, предложенный в [13] и позднее усовершенствованный в [14–16]; в этом методе на входной границе задаются случайные двумерные вихри, которые затем эволюционируют в расчетной области. Эти вихри определяются реалистическими значениями временнóго и пространственного масштабов, а также величиной завихренности, вытекающими из информации о среднем течении и расстояниях между узлами сетки. В работе [17] разработан метод, основанный на синтезировании случайных бездивергентных турбулентных скоростей с учетом спектральных и когерентных функций, соответствующих статистическим данным о течениях в атмосферных пограничных слоях. Эта схема обеспечивает выполнение спектра турбулентности и когерентной функции.

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

1.2. Сеточное разрешение

Несмотря на наличие мелкомасштабных вихрей в устойчивых пограничных слоях, в некоторых работах, выполненных методом LES, расчеты проводились на сетках с малым разрешением [2, 6]. В работе [18] выполнен ряд расчетов таких слоев методом крупных вихрей в широком диапазоне шагов сетки. При этом было обнаружено, что как средние, так и турбулентные параметры течения существенно зависят от шага сетки, даже до такой степени, что от его величины зависит сама сходимость решения. В другом расчете устойчивого пограничного слоя методом LES [19] использованы изотропные сетки с переменным шагом и показано, что величина подсеточной ТКЭ существенно зависит от шага сетки. На грубых сетках в численном решении доминировали подсеточные потоки с весьма малыми значениями ТКЭ, вплоть до не поддающихся разрешению. В работе [19] сообщается, что увеличение разрешения сетки ограничивает доминирование подсеточных потоков очень узким пристеночным слоем, тогда как в остальной части области турбулентные вихри могут быть успешно разрешены. Увеличение уровня устойчивости пограничного слоя требует еще большего сеточного разрешения и бóльших вычислительных ресурсов, так как доминирующие вихри становятся мельче и в большей степени обладают свойством перемежаемости [18].

1.3. Подсеточная модель

Наиболее широко используемая подсеточная модель – это модель замыкания при помощи вихревой вязкости, предложенная Смагоринским [20]. В этой модели перенос импульса в неразрешенном поле скорости осуществляется эффективной вязкостью [5]. Эта модель была успешно использована во многих работах по LES-моделированию нейтральных и конвективных атмосферных пограничных слоев [11, 14]. К другим моделям подсеточной вязкости принадлежат адаптированная к стенке модель локальной вихревой вязкости (Wall-Adapted Local Eddy-viscosity, WALE) [21] и получившая недавно популярность подсеточная модель турбулентной кинетической энергии [22, 23].

Для устойчивых пограничных слоев имели место попытки исследовать потенциальное влияние устойчивости пограничного слоя на показатели работы модели Смагоринского. Так, в работе [24] выполнены расчеты устойчивого пограничного слоя методом крупных вихрей с использованием модели Смагоринского, дополненной функциями, описывающими затухание на стенке и вводящими поправки на устойчивость. Недостатки моделирования были объяснены свойствами подсеточной модели и чрезмерными требованиями к сеточному разрешению при нехватке вычислительных ресурсов. В [25] при решении аналогичной задачи использованная подсеточная модель состояла из двух частей; при этом возникали численные неустойчивости, приписанные как численной процедуре, так и самой подсеточной модели. В работе [26] выполнено сравнение двух подсеточных моделей, зависимой и не зависимой от масштаба, при расчете устойчивого пограничного слоя методом крупных вихрей и исследована зависимость этих моделей от устойчивости и высоты над землей. Сделан вывод, что модель, зависимая от масштаба, дает правильное значение коэффициента Смагоринского с поправкой на устойчивость, а модель, не зависимая от масштаба, занижает это значение. В более широком LES-исследовании [27] изучено влияние подсеточных потоков и диссипаций на мелкомасштабную турбулентность в устойчивом пограничном слое с привлечением данных специально разработанного полевого эксперимента. Обнаружено, что устойчивость не изменяет долю подсеточных потоков в их общем числе, но коэффициенты подсеточной модели значительно меняется в зависимости от уровня устойчивости. Несмотря на эти результаты, во многих случаях имело место успешное LES-моделирование устойчивых пограничных слоев с использованием простой подсеточной модели Смагоринского [18, 19, 28].

1.4. Моделирование приповерхностного слоя

Для большей экономичности численных расчетов при моделировании приповерхностных течений широко используются функции стенки. Их использование основано на гипотезе о подобии течений в пограничных слоях на гладкой и шероховатой стенках [9]. При моделировании атмосферных пограничных слоев обычно используется функция стенки, основанная на линейном масштабе аэродинамической шероховатости z₀ [29]. Функции стенки обычно подчиняются логарифмическому закону, т.е. имеет место линейное соотношение между ${{U}^{ + }} = U{\text{/}}{{u}_{\tau }}$ и логарифмом величины ${{z}^{ + }} = z{{u}_{\tau }}{\text{/}}\nu $, где ${{u}_{\tau }}$ – скорость трения, а $\nu $ – кинематическая вязкость.

Аналогично функциям стенки для скорости, для сокращения вычислительных затрат широко используются функции стенки для температуры. Наиболее известная функция стенки для температуры основана на подобии между полями импульса и теплопереноса [30]. Имеется ряд уточненных функций стенки для температуры, соответствующих различным режимам термической устойчивости. На основе асимптотической комбинации естественной и вынужденной конвекции в работе [31] были предложены функции стенки для течений в вертикальном канале между параллельными пластинами, подходящие для расчета смешанной конвекции при всех значениях чисел Ричардсона. В работе [32] также была предложена температурная функция стенки, параметрически зависящая от числа Ричардсона; она использовалась в расчетах естественной и вынужденной конвекции в атмосферном пограничном слое при обтекании тупых тел, установленных на поверхностях. Отмечено, что использование этой функции приводило к более точным результатам, чем при использовании функции стенки, соответствующей вынужденной конвекции. Наряду с этими достижениями имеется и много примеров успешного LES-моделирования атмосферного пограничного слоя при различных уровнях температурной стратификации на основе классической температурной функции стенки, описанной в [30] (см., например, [23, 33]).

1.5. Гибридные методы крупных вихрей и очень крупных вихрей

Ввиду больших вычислительных затрат, связанных с применением методов крупных вихрей, особенно в областях больших размеров, что свойственно метеорологическим приложениям, разрабатывались также гибридные методы крупных вихрей с целью сократить потребные расходы при сохранении точности расчетов. В литературе описаны многие гибридные методы, хотя терминология в этой области достаточно нечеткая, так же как и отнесение тех или иных методов к гибридным. Наиболее известная группа гибридных методов известна как гибридные RANS-LES методы; в этих методах модель RANS используется в части расчетной области [34]. Эти модели разделяют на два основных класса, а именно, зональные и незональные методы, в соответствии со стратегией определения границы между областями, описываемыми моделями RANS и LES [35]. В зональных методах расчетная область разбивается на две зоны, разделяемые границей, определяемой как ряд фиксированных геометрических поверхностей. RANS- и LES-моделирования проводятся по отдельности в соответствующих зонах, а затем проводится стыковка полученных решений. Основной недостаток этого метода состоит в необходимости задания сложных условий стыковки на границе раздела [36]. В незональном подходе граница между областями, управляемыми RANS и LES, устанавливается заданием специальных граничных условий. Например, в незональном методе отсоединенных вихрей (detached-eddy simulation, DES) подход RANS применяется во всем присоединенном пограничном слое или в его большей части, тогда как в отрывных зонах выполняется LES-моделирование [37]. Незональный подход более прост в реализации, причем имеется множество его версий, обладающих повышенной точностью. В группе незональных методов, названных методами моделирования стенки (wall-modeling in LES, WMLES), вклад RANS ограничен весьма тонким пристеночным слоем. Именно этот подход используется в настоящей работе.

Наряду с численными методами, используемыми в гибридных LES, метод очень крупных вихрей (VLES) также представляет собой мощный инструмент экономии ресурсов. Концепция VLES, первоначально предложенная в работе [38], представляет собой один из самых ранних гибридных методов. Основное различие между VLES и стандартным LES состоит в определении ширины фильтра относительно шага сетки. В чистом LES ширина фильтра связана с шагом сетки, тогда как во VLES она может быть задана произвольно, имея любое значение от шага сетки до характерного линейного размера течения [35, 37]. Увеличение ширины фильтра снижает точность моделирования [35, 39]. Это объясняется тем, что увеличение ширины фильтра во VLES уменьшает отношение количества разрешенных и моделируемых вихрей. Влияние роста ширины фильтра на вычислительные затраты неясно, так как во VLES ширина фильтра может меняться без изменения шага сетки. Согласно данному определению, метод VLES становится просто LES, когда ширина фильтра берется как нижний предел шага сетки. В работе [40] предложено численное определение различия между LES и VLES. Согласно этой работе, LES на достаточно подробной сетке и с малой шириной фильтра должен разрешать более, чем 80% кинетической энергии турбулентности во всей области, за исключением приповерхностных зон, где используются функции стенки. Напротив, VLES определяется как модель с достаточно грубой сеткой и длиной фильтра, которая разрешает менее 80% кинетической энергии турбулентности в расчетной области. В настоящей работе приводятся данные о долях разрешаемой и моделируемой кинетической энергии турбулентности с целью дать оценку численным характеристикам модели.

1.6. Проблемы и цели

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

Работа построена следующим образом. В разд. 2 дано описание метода VLES. В этом разделе обсуждаются детали трех основных составных частей метода, а именно, подсеточная модель, метод синтетических вихрей и функции стенки. Основные особенности численных схем и методологии метода VLES обсуждаются в разд. 3, а полученные результаты приведены в разд. 4. Наконец, в Заключении подведены итоги проведенных исследований и даны рекомендации для будущих работ.

2.1. Формулирование и реализация подсеточной модели

Рассматривается несжимаемое, турбулентное, термически неоднородное течение с однопараметрической подсеточной моделью турбулентности. В этих предположениях течение в рассматриваемой модели VLES описывается системой четырех уравнений для безразмерных переменных ${{U}_{i}}$, $\Theta $ и ${{k}_{{sgs}}}$. Эти уравнения включают уравнение неразрывности, уравнение переноса импульса, уравнение переноса тепла или энергии и уравнение для кинетической энергии подсеточной турбулентности. Ниже эти уравнения (уравнения (2.1)–(2.4)) приводятся в безразмерной форме, причем за характерный линейный размер выбрана высота пограничного слоя $\delta $, а за масштабы скорости и температуры взяты их значения в набегающем потоке U₀ и T₀

Хотя все члены в этих уравнениях подробно обсуждаются в других работах [19, 22, 23], ниже приводится их краткое описание. Черта над переменными указывает на их пространственную разрешенность; $\bar {p} = \bar {p}{\text{*}} + \tfrac{1}{3}{{\tau }_{{ii}}}$ есть пространственно разрешенное модифицированное кинематическое давление, нормализованное по постоянной плотности, где $\bar {p}{\text{*}}$ – пространственно разрешенное статическое давление; ${{k}_{{sgs}}}$ – кинетическая энергия подсеточной турбулентности; ${{\tau }_{{ij}}} = \overline {{{U}_{i}}{{U}_{j}}} - {{\bar {U}}_{i}}{{\bar {U}}_{j}} = - 2{{\nu }_{T}}{{\bar {S}}_{{ij}}}$ – поток подсеточного импульса, где ${{\bar {S}}_{{ij}}} = \left( {\tfrac{{\partial {{{\bar {U}}}_{i}}}}{{\partial {{x}_{j}}}} + \tfrac{{\partial {{{\bar {U}}}_{j}}}}{{\partial {{x}_{i}}}}} \right)$ – скорость деформации; π_i = $\overline {{{U}_{i}}\Theta } - {{\bar {U}}_{i}}\bar {\Theta } = - {{\nu }_{\theta }}\tfrac{{\partial{ \bar {\Theta }}}}{{\partial {{x}_{i}}}}$ – подсеточный кинематический тепловой поток; Re = ${{U}_{0}}\delta {\text{/}}\nu $ – число Рейнольдса; ${\text{R}}{{{\text{e}}}_{T}} = {{U}_{0}}\delta {\text{/}}{{\nu }_{T}}$ – число Рейнольдса в модели подсеточной турбулентности; ${\text{Pr}} = \nu {\text{/}}\alpha $ – ламинарное число Прандтля; ${\text{Ri}} = g\delta \Delta \Theta {\text{/}}(\bar {\Theta }U_{0}^{2})$ – число Ричардсона; $P = - {{\tau }_{{ij}}}{{\bar {S}}_{{ij}}}$ – производство сдвигового напряжения; $B = - g{{\nu }_{\theta }}\tfrac{{\partial{ \bar {\Theta }}}}{{\partial z}}$ – производство плавучести; ε = ${{C}_{\varepsilon }}k_{{sgs}}^{{3/2}}{\text{/}}l$ – скорость диссипации и δ_ij – дельта-функция Кронекера.

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

где ${{C}_{\Delta }}$ – определяющий параметр масштаба l и, следовательно, всей подсеточной модели. Модель турбулентности замыкается заданием остальных параметров, включая турбулентную вязкость ${{\nu }_{T}} = {{C}_{k}}k_{{sgs}}^{{1/2}}l$ и турбулентную тепловую вязкость ${{\nu }_{\theta }} = {{\nu }_{T}}{\text{/P}}{{{\text{r}}}_{T}}$. Значение ${{C}_{k}}$ берется равным 0.094, ${{C}_{\varepsilon }}$ = 1.048 и турбулентное число Прандтля ${\text{P}}{{{\text{r}}}_{T}}$ = 0.85.

2.2. Метод синтетических вихрей

В настоящей модели VLES для генерации турбулентности во входном сечении расчетной области используется метод синтетических вихрей. Первоначально метод был разработан в [13] и с тех пор непрерывно совершенствовался в работах [14–16]. Основная концепция данного метода заключается в том, чтобы генерировать пульсации скорости в виде синтетических вихрей, исходя из осредненной статистической информации о поле течения как функции пространства (высоты над землей) и времени. Определяющими параметрами являются количество вихрей, размер каждого вихря, завихренность (или эквивалентное поле скорости, характеризующее каждый вихрь) и время жизни вихрей [15]. Вихри двумерны, а вектор завихренности параллелен направлению течения. Схематически генерация вихрей изображена на рис. 1a. Согласно этой теории, полностью разработанной в литературе [13—16], поле пульсаций скорости на данном шаге по времени определяется функцией

где ${\mathbf{u}}$ – возмущение скорости на входе, которое затем накладывается на среднюю скорость на входе, ${\mathbf{x}}$ – положение вектора точки на входной границе, $N$ – количество вихрей, вводимых на входной границе (в настоящей работе принято $N$= 200), $i$ – индекс текущего вихря, ${{\Gamma }_{i}}$ – циркуляция текущего вихря, ${{{\mathbf{x}}}_{i}}$ – вектор, определяющий положение центра текущего вихря, s – единичный вектор в направлении потока и σ_i(x_i) – характерный радиус текущего вихря. Предполагая, что поток направлен в направлении +x, а направление, нормальное к стенке, есть $ + z$, зададим степенной профиль средней скорости на входе в виде [43] где ${{z}_{{ref}}}$ – характерная высота, ${{U}_{{ref}}}$ – характерная скорость и $\alpha = 1{\text{/}}ln({{z}_{{ref}}}{\text{/}}{{z}_{0}})$ – показатель степени, зависящий от ${{z}_{{ref}}}$ и длины аэродинамической шероховатости z₀ [5].

Размер вихрей, несущих энергию, задается на основе теории длины пути смешения [9] в виде

где $\kappa = 0.41$ – постоянная Кармана, а ${{\sigma }_{{max}}}$ – размер наибольшего из несущих энергию вихрей, отнесенный к характерному размеру границы $L$ [9]. Таким образом, ${{\sigma }_{{max}}}$ и L связаны соотношением ${{\sigma }_{{max}}} = {{a}_{\sigma }}L$, где постоянная a_σ будет определена ниже.

Используя масштабный анализ, можно оценить характерное время несущих энергию вихрей. Для наибольших несущих энергию вихрей с характерной скоростью U₀ и линейным размером ${{\ell }_{0}} = {{\sigma }_{{max}}}$ характерное время можно получить в виде

где $\varepsilon $ – скорость диссипации, определенная и подробно описанная в [9]. Этот временной масштаб характерен не для всех несущих энергию вихрей, но лишь для наибольших из них. Временной масштаб, характерный для всех вихрей, можно определить соотношением $\tau = {{a}_{\tau }}{{\tau }_{0}}({{\ell }_{0}})$, где постоянная ${{a}_{\tau }}$ будет определена ниже.

Температура на входе задается следующим степенным профилем

где ${{\Theta }_{\infty }}$ – температура на верхней границе пограничного слоя, ${{\Theta }_{s}}$ – температура поверхности, ${{z}_{{max}}}$ – высота области и α – показатель степени, зависящий от длины аэродинамической шероховатости, как определено выше уравнением (2.7). Ввиду недостатка места, здесь не приводится описание профиля интенсивности турбулентности и циркуляции вихрей на входной границе; соответствующие данные можно найти в [9].

2.3. Задание функции стенки

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

где $\kappa = 0.41$ – постоянная Кармана. В библиотеке OpenFOAM эта функция стенки имеет название nutkAtmRoughWallFunction.

Температурная функция стенки в настоящей модели VLES первоначально взята из [5]; она связывает величину ${{\Theta }^{ + }} = ({{\Theta }_{s}} - \Theta )\rho {{c}_{p}}{{u}_{\tau }}{\text{/}}{{q}_{s}}$, где q_s – тепловой поток к поверхности, и логарифм z⁺ линейным соотношением вида

где ${{\kappa }_{\theta }} = 0.48$ – тепловая постоянная Кармана, а ${{B}_{\theta }} = 3.9$ – константа модели турбулентности. Однако более принято задавать ${{\Theta }^{ + }}$ как функцию ${{U}^{ + }}$ в виде где ${\text{P}}{{{\text{r}}}_{T}} = 0.85$ – турбулентное число Прандтля, а величина Pr_f  описана в [30]. В библиотеке OpenFOAM эта функция стенки известна как alphatJayatillekeWallFunction.

Для кинетической энергии турбулентности использована функция стенки, известная в OpenFOAM как kqRWallFunction [41]

где ${{C}_{\mu }} = 0.09$ – постоянная.

3.1. Расчетная область и генерация сетки

Схематически расчетная область изображена на рис. 1б. Длина, ширина и высота области составляют $X = 5$ м, $Y = 1.5$ м и Z = 1.2 м соответственно. Поперечное сечение области соответствует сечению аэродинамической трубы, использованной в экспериментах [44]. Поток воздуха протекает в направлении $ + x$ с характерными скоростью и температурой, заданными во входном сечении, как показано на рис. 1б. На этом же рисунке представлены типичные профили скорости и температуры в пограничном слое. Результаты моделирования анализируются по данным, полученным на шести вертикальных линиях, на каждой из которых расположено большое количество датчиков. Типичный профиль показан на рис. 1б. Все шесть линий расположены в середине области $y = 0$, занимая всю область по высоте от z = 0 до z = 1.2 м; в продольном направлении они занимают положения от $x = 0$ до $x = 5$ м с шагом в 1 м.

Для исследования влияния шага сетки на результаты, полученные с помощью модели VLES, рассмотрены четыре сетки с различным разрешением. Их параметры представлены в табл. 1; их разрешение меняется от очень точного с 1 000 000 контрольных объемов до весьма грубого с 62 500 объемов. Результаты исследования влияния сеточного разрешения представлены в разделе 4.1.1. Сетка в направлениях $x$ и $y$ равномерная, а в направлении z шаг сетки меняется, со сгущением в окрестности поверхности. Сетка генерируется при помощи программы blockMesh, содержащейся в библиотеке OpenFOAM. Высота ячеек сетки у поверхности жестко контролируется и меняется отдельно от разбиения во внутренней части области. Поэтому влияние подсеточной модели и функций стенки обсуждается отдельно. Значения величины z⁺ в первом слое приводятся в разделе 4.2, где используются функции стенки; в остальных случаях, когда течение у стенки определяется численным решением (раздел 4.1), ${{z}^{ + }} < 1.5$.

3.2. Граничные условия

Во входном сечении расчетной области два отдельных степенных профиля для скорости и температуры задаются уравнениями (2.7) и (2.10). Что касается турбулентных пульсаций, метод синтетических вихрей используется лишь для генерации пульсаций скорости в соответствии с уравнением (2.6). Как показано ниже, ниже по потоку течение становится турбулентным относительно пульсаций как скорости, так и температуры. В выходном сечении и на верхней границе области ставятся условия нулевых производных для всех переменных. Также для всех переменных ставятся периодические граничные условия на передней и задней границах расчетной области. На дне области ставятся условия прилипания для скорости и условие постоянной температуры.

Для подсеточной турбулентной кинетической энергии на входе ставится условие, задаваемое процедурой atmBoundaryLayerInletK. Это условие предполагает, что вся входная граница представляет собой инерционный поверхностный слой атмосферного пограничного слоя, так что скорость трения и турбулентная кинетическая энергия не зависят от высоты [3]. В этом граничном условии сначала в предположении логарифмического закона рассчитывается скорость трения

а затем равномерная подсеточная турбулентная кинетическая энергия в виде ${{k}_{{sgs}}} = u_{\tau }^{2}{\text{/}}C_{\mu }^{{1/2}}$, где ${{C}_{\mu }} = 0.09$ – постоянная. На выходе для подсеточной турбулентной кинетической энергии ставится условие нулевой производной. На поверхности возможна постановка двух граничных условий, а именно либо условие нулевой кинетической энергии при численном определении течения у стенки, либо условие, задаваемое процедурой kqRWallFunction, при использовании стандартных функций стенки для шероховатых поверхностей. В работе используются оба эти граничных условия с целью исследования влияния функций стенки на результаты, получаемые при помощи модели VLES (раздел 4.2).

Для турбулентной вязкости условие нулевой нормальной производной ставится во входном и выходном сечениях. На стенке возможна постановка двух условий. Если течение у поверхности рассчитывается численно, то ставится условие нулевой нормальной производной. В противном случае может быть поставлено граничное условие для шероховатой поверхности, реализуемое процедурой nutkAtmRoughWallFunction. Это условие модифицирует турбулентную вязкость у поверхности таким образом, что

где $\bar {E} = (z + {{z}_{0}}){\text{/}}{{z}_{0}}$. Температурная функция стенки, описываемая уравнением (2.13), модифицирует турбулентную тепловую диффузию у стенки следующим образом где α – коэффициент молекулярной тепловой диффузии. В библиотеке OpenFOAM это граничное условие известно как alphatJayatillekeWallFunction.

3.3. Численные схемы

В расчетах используется неявная обратная схема второго порядка по времени; все градиентные схемы основаны на гауссовом интегрировании второго порядка с линейной интерполяцией. Все лапласовы схемы основаны на уточненном гауссовом интегрировании с линейной интерполяцией, которое обеспечивает неограниченное консервативное численное поведение второго порядка. Дивергентные схемы основаны на гауссовом интегрировании с линейной либо противопоточной интерполяцией, в зависимости от конкретной переменной [41].

Во всех расчетах шаг по времени выбирается таким образом, что максимальное число Куранта удовлетворяет условию $Co = \Delta t\left| {\bar {U}} \right|{\text{/}}\Delta x < 1$. Разрешение матрицы давления, предобусловленной согласно неполному диагональному методу Холески, проводится при помощи предобусловленного сопряженного градиентного солвера. Предобусловливание других переменных проводится диагональным неполным нижне-верхним методом, а решение проводится при помощи пред-обусловленного двусопряженного градиентного солвера. Уравнения, связанные с давлением, решаются гибридным методом, включающим два алгоритма, а именно неявный метод расщепления оператора и полунеявный метод [41].

Имеется большое число работ по исследованию влияния численных схем на моделирование с разрешением масштабов [45–50]. В работе [45] изучено влияние нескольких разностных схем, в том числе центральных схем второго и четвертого порядка, а также спектрального метода, на результаты расчетов нестационарного слоя смешения методом крупных вихрей. Сделан вывод, что спектральная схема дает несколько лучшие результаты, чем разностные, но при этом она требует бóльших вычислительных ресурсов. В [46] обсуждаются детали двух широко используемых численных схем, применяемых при LES-моделировании сложных конфигураций, а именно, метод погруженных границ и неструктурированные сетки. Взаимодействие между численной дискретизацией и подсеточным моделированием изучено в [48]. В работе [48] использованы неявные подсеточные модели, основанные на конечно-объемной дискретизации, что создает возможность полного объединения дискретизации и подсеточного моделирования. Предложены схема нелинейной дискретизации для подсеточной модели в областях турбулентного течения и схема второго порядка в областях ламинарного течения. Влияние критериев сходимости и шага по времени на расчеты методом LES некоторого течения, ограниченного стенкой, изучено в [49]. Обнаружено, что моделирование с более мягким критерием сходимости более эффективно с вычислительной точки зрения. В работе [50] изучены ошибки численного счета и моделирования при использовании явных и неявных динамических разностных схем при расчетах методом LES трехмерного вихревого течения Тэйлора–Грина. Отмечается, что динамические разностные схемы обеспечивали оптимальную точность на всех разрешенных масштабах в потоке в любой момент времени и создавали меньшие вычислительные ошибки, чем стандартные асимптотические конечно-разностные схемы. Численные схемы, используемые в настоящей работе, обеспечивают приемлемые результаты для течений в различных условиях. Поэтому влияние различных разностных схем на результаты моделирования в настоящей работе не исследуется и интересующийся читатель отсылается по этому вопросу к исследованиям, содержащимся в работах [45–50].

4.1. Моделирование устойчивого пограничного слоя с разрешением у стенки

В этом разделе содержатся результаты моделирования с непосредственным решением у стенки для четырех случаев устойчивости, описанных в табл. 2. Сначала в разделах 4.1.4–4.1.3 анализируется чувствительность численной модели для варианта 1 в табл. 2. В соответствии с минимизационным подходом, принятым в данной работе, результаты этого анализа позволят свести модель к ограниченному числу входных параметров, что будет использовано при анализе устойчивости в разделах 4.1.4 и 4.1.5.

4.1.1. Исследование сеточного эффекта. Предполагается, что модель VLES должна с хорошей точностью предсказывать средние характеристики течения, а также турбулентные флуктуации и потоки, на грубых сетках. Таким образом, одним из важнейших тестов модели VLES является ряд расчетов на различных сетках, от мелких до очень грубых. Показатели работы модели VLES на четырех сетках, приведенных в табл. 1, представлены на рис. 2а–2в. Ввиду недостатка места, на рисунках представлены лишь некоторые переменные. Поскольку проведение теста на характеристики сетки для всех вариантов устойчивости потребовало бы чрезмерных вычислительных ресурсов, такой тест был проведен лишь для одного варианта. На рисунках также приведены экспериментальные данные [44]. Можно видеть, что результаты численного моделирования находятся в приемлемом согласовании с экспериментом, начиная с сетки уровня III. Для очень грубой сетки уровня IV имеет место большое отличие от экспериментальных данных, особенно, в том что касается турбулентных переменных. Эти результаты показывают, что уровень сетки III представляет собой предельный по грубости сетки уровень, при котором еще возможно правильное определение средних и турбулентных параметров пограничного слоя. Поэтому данная сетка использована в дальнейших численных расчетах настоящей работы.

Помимо описанного выше исследования, нужно было убедиться, что использованный шаг сетки способен адекватно разрешить турбулентные структуры. С этой целью проведено сравнение шага сетки с размером наибольшего из несущих энергию вихрей, который является интегральным линейным масштабом турбулентности. Оценка интегрального линейного масштаба турбулентных вихрей внутри расчетной области ($z > 0.05\delta $) для четырех рассмотренных вариантов устойчивости дает величину от 0.15 до 0.7 м, причем эта величина возрастает с удалением от стенки. Процедура оценки интегральных пространственных и временных масштабов заимствована из работ [5, 51]. Сетка, использованная в данных расчетах, имела постоянные шаги $\Delta x$ = = 0.05 м и $\Delta z = 0.03$ м, а размер шага в направлении z менялся со сгущением сетки у поверхности. Так, шаг сетки $\Delta z$ был равен 0.001 м у стенки и возрастал до 0.095 м в области свободного течения. Из этих данных можно заключить, что использованный шаг сетки адекватен моделированию турбулентного течения. Что касается временных масштабов, использовались шаг по времени 0.005 с и время осреднения различных переменных 100 с. Для четырех рассмотренных вариантов устойчивости интегральный масштаб времени внутри расчетной области ($z > 0.05\delta $) изменялся от 0.14 до 0.76 с. Имея в виду, что шаг по времени был в от 28 до 152 раз меньше интегрального масштаба времени, а выбранное время осреднения было по крайней мере в 100 раз дольше, чем интегральный масштаб времени, можно заключить, что использованные в данных расчетах шаг по времени и время осреднения адекватны для разрешения турбулентных структур.

4.1.2. Зависимость от параметров синтетического метода вихрей. Как сказано выше, входная турбулентность генерировалась синтетическим методом вихрей, который позволяет минимизировать количество определяющих параметров лишь двумя параметрами ${{a}_{\sigma }}$ и ${{a}_{\tau }}$. В этом разделе исследуется чувствительность модели к этим двум параметрам. Сначала исследуется влияние на модель числового параметра ${{a}_{\sigma }}$, определяющего размеры вихрей на входе. Далее исследуются показатели работы модели в зависимости от численного параметра ${{a}_{\tau }}$, определяющего, как часто появляются на входе новые вихри.

На рис. 2г–2е показано, что изменение ${{a}_{\sigma }}$ не влияет на среднюю скорость. С другой стороны, возрастание ${{a}_{\sigma }}$ от 1 до 5 существенно увеличивает турбулентную статистику. Это объясняется тем фактом, что более высокие значения ${{a}_{\sigma }}$ соответствуют более крупным вихрям, которые несут в себе бóльшую энергию и вносят больший вклад в уровень турбулентности моделируемого течения. Хотя и трудно подобрать идеальное значение ${{a}_{\sigma }}$, которое давало бы хорошее соответствие с экспериментальными данными по всей высоте пограничного слоя, представляется, что значение ${{a}_{\sigma }} = 3$ дает наилучшее согласование с экспериментом.

Анализ чувствительности синтетического метода по отношению к параметру ${{a}_{\tau }}$, представленный на рис. 2ж–2и, показывает те же тенденции, что и в случае ${{a}_{\sigma }}$, т.е. пренебрежимо малое влияние на средние переменные и существенное влияние на турбулентную статистику. Ввиду того факта, что ${{a}_{\tau }}$ определяет время жизни наибольших энергонесущих вихрей на входе, можно ожидать, что бóльшие значения ${{a}_{\tau }}$ будут соответствовать бóльшим значениям турбулентной статистики. В целом по данным, приведенным на рис. 2ж–2и, можно заключить, что наилучшее совпадение с трубным экспериментом имеет место при ${{a}_{\tau }} = 0.01$.

4.1.3. Зависимость от параметров подсеточной модели. Тестирование чувствительности подсеточной модели осуществляется изменением постоянной ${{C}_{\Delta }}$, определяющей линейный масштаб l модели (см. уравнение (2.5)), как показано на рис. 3. Можно видеть, что изменение ${{C}_{\Delta }}$ практически не влияет на средние переменные (рис. 3а) и турбулентные флуктуации (рис. 3б). Однако влияние ${{C}_{\Delta }}$ на турбулентные потоки импульса (рис. 3в) и тепла (рис. 3г и 3д) больше. Это объясняется тем фактом, что данный параметр управляет явлениями переноса на подсеточных масштабах, определяя подсеточную длину пути смешения. Из данных на рис. 3г и 3д можно заключить, что все три значения ${{C}_{\Delta }}$ завышают значение вертикального турбулентного теплового потока в придонной части пограничного слоя, тогда как значения горизонтального турбулентного теплового потока рассчитываются с приемлемой точностью. Это затрудняет правильный выбор данного параметра. В этом выборе может быть полезным анализ отношений значения турбулентной кинетической энергии, полученного при моделировании, к полному значению этого параметра. График этой величины, построенный на рис. 3е, показывает, что с увеличением ${{C}_{\Delta }}$ это отношение возрастает внутри расчетной области. Оно велико вблизи стенки, чего и следует ожидать, так как здесь течение по характеру близко к вязкому подслою. Обнаружено, что для успешной работы метода VLES не более 20% турбулентной кинетической энергии внутри расчетной области должно быть получено моделированием, а более 80% должно быть получено разностным решением [40]. Имея в виду этот критерий (рис. 3е) вместе с информацией, приведенной на рис. 3а–3д, можно заключить, что наиболее подходящим значением в численной модели, описывающей приповерхностные течения, является значение ${{C}_{\Delta }} = 1$.

4.1.4. Варианты тепловой устойчивости. После определения надлежащих числовых значений параметров модели, что было сделано в предыдущих разделах, метод VLES был применен для изучения эффекта тепловой стратификации атмосферного пограничного слоя. На рис. 4 представлены результаты моделирования с численным определением течения у стенки для четырех вариантов устойчивости, приведенных в табл. 2, а также сравнения с экспериментальными результатами [44]. Эти расчеты выполнены на сетке уровня III с параметрами метода синтетических вихрей ${{a}_{\sigma }} = 3$ и ${{a}_{\tau }} = 0.01$ и параметром подсеточной модели ${{C}_{\Delta }} = 1$.

По профилям средней скорости, представленным на рис. 4а, можно видеть, что с увеличением уровня устойчивости от слабо устойчивого варианта ($R{{i}_{\delta }} = 0.12$) до сильно устойчивого ($R{{i}_{\delta }}$ = = 0.74) дефицит скорости возрастает под влиянием тепловой стратификации. Профили средней температуры на рис. 4б сдвигаются в сторону бóльших значений, причем этот рост незначителен для трех первых вариантов устойчивости, но заметен для сильно устойчивого случая с $R{{i}_{\delta }} = 0.74$.

Что касается турбулентных пульсаций горизонтальной скорости, вертикальной скорости и температуры, представленных на рис. 4в–4д, очевидно, что рост тепловой устойчивости приводит к снижению турбулентных пульсаций. Это связано со стабилизирующим эффектом сил плавучести, обусловленным тепловой стратификацией: они подавляют турбулентность в пограничном слое. В особенности эти силы снижают пульсации вертикальной компоненты скорости, как указывалось еще в [1]. Это можно видеть и по результатам настоящей работы, на рис. 4в и 4г, где турбулентные пульсации компоненты $w$ примерно в два раза меньше турбулентных пульсаций компоненты $u$ в рассмотренных течениях в пограничном слое.

Что касается турбулентного импульса и тепловых потоков, стабилизирующий эффект тепловой стратификации при увеличении числа Ричардсона хорошо виден на рис. 4е–4з. Еще одна характерная особенность обеих групп турбулентных пульсаций на рис. 4в–4д и турбулентных потоков на рис. 4е–4з – это различное поведение параметров течения в слабо устойчивых случаях ($R{{i}_{\delta }} = 0.12$ и 0.24) и в сильно устойчивых случаях ($R{{i}_{\delta }} = 0.40$ и 0.74). В двух слабо устойчивых случаях параметры течения имеют существенно ненулевые значения в придонной половине пограничного слоя и начинают стремиться к нулю лишь на очень малой высоте над поверхностью, примерно при $z = 0.05\delta $. В двух сильно устойчивых случаях это стремление к нулю начинается на бóльших высотах, где значения параметров турбулентности очень близки к нулю. Это является признаком подавления турбулентности при сильной тепловой стратификации, что можно наблюдать также на рис. 5, где мгновенные контуры $y$-компоненты завихренности построены в центральной плоскости $x - z$ для всех четырех вариантов устойчивости. Можно видеть, что в двух слабо устойчивых случаях ($R{{i}_{\delta }} = 0.12$ и 0.24) в большей части области течения доминируют крупные вихри, о чем свидетельствуют высокие уровни завихренности. Напротив, в двух сильно устойчивых случаях ($R{{i}_{\delta }} = 0.40$ и 0.74) области высокой завихренности ограничены лишь небольшими приповерхностными зонами, тогда как внутри расчетной области турбулентные движения подавлены.

В целом по данным, приведенным на рис. 4, можно сказать, что согласование между данными численными результатами и экспериментальными данными [44] достаточно хорошее как для средних, так и для турбулентных профилей. Хотя некоторые расхождения и имеют место, модель VLES в целом хорошо предсказывает поведение пограничного слоя при различных уровнях устойчивости.

4.1.5. Спектральный анализ при численном моделировании течения у стенки. Хотя оценка гидродинамических вычислительных моделей в основном проводится в соответствии с полученными на их основе предсказаниями средних и турбулентных параметров течения, успешный спектральный анализ турбулентных пульсаций может придать большее доверие полученным результатам. Спектральный анализ четырех вариантов устойчивости, рассчитанных при помощи VLES с численным моделированием течения у стенки, представлен на рис. 6а–6в для турбулентных пульсаций горизонтальной скорости, вертикальной скорости и температуры соответственно. Коспектральный анализ турбулентного напряжения $\overline {uw} $, вертикального теплового потока $\overline {w\theta } $ и горизонтального теплового потока $\overline {u\theta } $ приведен на рис. 6г–6е для тех же вариантов устойчивости. Детали спектрального анализа, выполненного для результатов настоящего моделирования, можно найти в [3]. Для анизотропного течения в атмосферном пограничном слое предполагается, что в инерционном интервале для пульсаций скорости и температуры наклон кривой спектральной плотности энергии $E(\kappa )$ по отношению к волновому числу $\kappa $ в логарифмической шкале равен –5/3 (спектр Колмогорова) [3, 52]. Наклон кривой коспектральной плотности энергии $C(\kappa )$ для пульсаций скорости в продольном (x) и вертикальном (z) направлениях по отношению к волновому числу $\kappa $ в логарифмической шкале равен приблизительно –7/3 [52]. Для коспектральной плотности энергии турбулентных вертикального и горизонтального тепловых потоков этот наклон в инерционном интервале равен –7/3 и –5/2 соответственно [52]. Для сравнения эти наклоны изображены на рис. 6 для соответствующих величин.

Можно заметить, что, как показывают рис. 6а–6в, спектральный состав турбулентных пульсаций скорости и температуры разрешается моделью VLES для волновых чисел, изменяющихся почти на два порядка величины. Инерционный диапазон частично matched, в то время как меньшие масштабы инерционного диапазона моделируются (а не разрешаются), что проявляется в резком падении и усечении спектров. Кроме того, спектральные данные в инерционном диапазоне довольно хорошо следуют наклону –5/3. Из данных на рис. 6а–6в также следует, что рост тепловой устойчивости приводит к уменьшению амплитуды плотности спектральной энергии, что объясняется подавлением турбулентности силами плавучести.

Коспектры турбулентного импульса и тепловых потоков также успешно разрешаются для волновых чисел, изменяющихся почти на два порядка величины, и частично соответствуют инерционному интервалу (рис. 6г–6е). Из этих графиков также следует хорошее совпадение коспектральных данных с наклоном в логарифмической шкале в инерционном интервале. По этим же графикам видно, что коспектральная плотность между двумя соответствующими параметрами убывает с увеличением числа Ричардсона для всех четырех рассмотренных вариантов устойчивости.

4.2. Решение методом VLES с использованием функций стенки

В настоящем разделе обсуждаются результаты расчетов методом VLES течения над шероховатой поверхностью с использованием функций стенки для скорости и температуры с целью более экономного проведения расчетов. Для этого исследования выбрана сетка уровня III, а также используются числовые параметры синтетического метода вихрей и подсеточной модели, выбранные в предыдущем разделе для случая численного моделирования течения у стенки (${{a}_{\sigma }} = 3$, ${{C}_{\Delta }}$ = 1). Однако параметр, определяющий временной масштаб вихрей, имеет большее значение, ${{a}_{\tau }}$ = 0.05, в связи с тем, что при использовании функций стенки требуется большее значение характерного времени образования крупных вихрей. Это объясняется тем фактом, что при использовании функций стенки производство турбулентности у поверхности моделируется, а не разрешается численно, и поэтому образование вихрей на некотором расстоянии от стенки происходит с бóльшим значением временнóй константы. Соответственно, и перенос турбулентной кинетической энергии от стенки во внешние слои требует больше итераций по времени, прежде чем на входе возникнут новые вихри.

При численном моделировании с использованием функций стенки очень большое значение имеет выбор высоты первого слоя сетки. Поэтому расчеты методом VLES проводились при четырех различных значениях этой высоты. На рис. 7а–7г приведены данные по чувствительности модели к изменению ${{z}_{p}}$, половины высоты первой вычислительной ячейки для варианта устойчивости 1. Значения z⁺, соответствующие четырем высотам первого слоя ячеек, равны 87, 278, 325 и 411. Процедура расчета z⁺ описана в [9]. По данным на рис. 7а–7г можно видеть, что результаты настоящего моделирования находятся в удовлетворительном соответствии с экспериментальными данными [44] для средних и турбулентных параметров в областях, далеких от поверхности, хотя некоторые расхождения имеют место. Вблизи от стенки, где функции стенки моделируют турбулентность, два наименьших значения z⁺ (87 и 278), по-видимому, обеспечивают лучшее совпадение с экспериментом. Таким образом, можно заключить, что для успешного моделирования атмосферного пограничного слоя с температурной стратификацией величина z⁺ должна быть меньше 278. Ухудшение турбулентной статистики вблизи стенок при использовании функций стенки является ожидаемым и даже желательным фактом, поскольку использование этих функций приводит к моделированию, а не численному разрешению пульсаций. Поэтому данное обстоятельство не следует интерпретировать как слабое место модели. На самом деле, чтобы преодолеть это обстоятельство, следует лучшим образом подбирать входные параметры синтетического метода, но это находится за рамками настоящей работы. Тем не менее ошибки, получаемые при использовании функций стенки с двумя наименьшими значениями z⁺, имеют тот же порядок, что и при численном разрешении приповерхностного течения.

Спектральный анализ результатов расчетов методом VLES с использованием функций стенки для расчетной области, соответствующей z⁺ = 278, и варианта устойчивости 1 приведен на рис. 7д–7ж. На этих же рисунках изображены наклоны зависимостей различных параметров в логарифмической шкале, обсуждаемые в разделе 4.1.5. По спектрам и коспектрам скорости и температуры на рис. 7д–7ж видно, что модель VLES способна разрешать каскады энергии в диапазонах волновых чисел протяженностью в два порядка. Инерционный интервал частично соответствует представленным спектральным и коспектральным кривым. Также можно отметить удовлетворительное согласование этих спектров и коспектров с соответствующими наклонами в инерционном интервале.