Large-eddy simulation of a droplet-laden air-flow over a waved water surface

全文:

1. Введение

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

Воздушный пограничный слой над взволнованной водной поверхностью в натурных и лабораторных условиях как правило характеризуется наличием капель брызг. Результаты экспериментов показывают, что инжекция капель в воздух обусловлена многими факторами, такими как всплывающие к поверхности и лопающиеся пузырьки, обрушение (или микрообрушение) поверхностных волн, а также гидродинамическая неустойчивость индуцируемого ветром дрейфового течения, приводящая к росту возмущений на водной поверхности с их последующей атомизацией [Bortkovskii, 1987; Andreas et al., 2010; Troitskaya et al., 2017, 2023]. Результаты экспериментов [Andreas et al., 2010; Troitskaya et al., 2017] говорят о том, что основной вклад в объемную концентрацию водной фазы привносят капли с диаметром в диапазоне от 100 до 300 микрон.

Определение вклада капель в процессы обмена в лабораторных и натурных экспериментах сопряжено со значительными трудностями. Поэтому в моделях прогноза, как правило, используются оценки вклада капель в потоки импульса и тепла, основанные на балковых параметрах задачи, таких как скорость ветра на высоте 10 м от поверхности воды, значительная высота волн, разность температур воды и воздуха, а также средняя скорость гравитационного оседания капель. Подобные оценки существенно зависят от задаваемых параметров и предположений о динамике капель и могут приводить к значительному количественному разбросу модельных предсказаний (см., например, обзор в [Andreas et al., 2015]).

Более детальную картину динамики капель дают лагранжевы стохастические модели (ЛСМ) [Edson et al., 1994; Mueller, Veron, 2014; Troitskaya et al., 2016]. В ЛСМ решаются уравнения, где поле скорости несущего потока для каждой капли задается в виде суммы средней скорости, определяемой профилем согласно предсказанию теории подобия атмосферного погранслоя, и гауссовых флуктуаций. К недостаткам ЛСМ можно отнести неспособность искусственно задаваемого стохастического форсинга учитывать влияние на динамику капель когерентных вихревых структур, наличие которых в погранслойных течениях хорошо известно из наблюдений в многочисленных экспериментах (см., например, обзор [Robinson, 1991]).

В последнее время для исследования процессов обмена в атмосферном погранслое все чаще используется метод прямого численного моделирования (“Direct numerical simulation” или сокращенно DNS) [Peng, Richter, 2019; Druzhinin et al., 2017]. При этом интегрируются полные уравнения несущей и дискретной фаз без привлечения каких-либо процедур осреднения и гипотез замыкания. В таких численных экспериментах процессы обмена полностью контролируются на всех физически значимых масштабах, вплоть до масштаба вязкой диссипации, и вклады капель определяются явно. Основным недостатком DNS, однако, является ограничение на число Рейнольдса несущего воздушного потока (основанное на характерных значениях длины волны, скорости ветра и кинематической вязкости воздуха) с максимальными рассматриваемыми значениями порядка Re ~ 10⁴, обусловленное конечными объемами вычислительных ресурсов.

Ограничение на число Рейнольдса в численных расчетах может не быть столь существенным при использовании метода вихреразрешающего моделирования (или LES, “Large Eddy Simulation”). В этом случае интегрируются уравнения, пространственно отфильтрованные по подсеточным (т.е., не разрешаемым заданной сеткой) масштабам течения. При этом в фильтрованных уравнениях движения присутствуют слагаемые, соответствующие подсеточным турбулентным напряжениям, расчет которых требует привлечения процедуры замыкания относительно разрешаемых, фильтрованных полей. Существует множество моделей замыкания, обладающих как достоинствами, так и недостатками (см., например, обзоры в [Piomelli, Balaras, 2002; Глазунов, 2009]).

В настоящей работе за основу взята модель LES, первоначально разработанная для исследования динамики однофазного атмосферного пограничного слоя над взволнованной водной поверхностью в [Sullivan et al., 2008]. Позднее эта модель была адаптирована в [Richter et al., 2019] для моделирования динамики капель в воздушном потоке над взволнованной водной поверхностью в канале со свободной границей. Однако, в отличие от [Richter et al., 2019], где стационарное течение в канале поддерживается постоянным градиентом давления, в настоящей работе используется постановка задачи, где стационарный, турбулентный воздушный поток создается постоянным средним напряжением сдвига, поддерживаемого движением верхней плоской границы с заданной постоянной скоростью (см. рис. 1 и обсуждение ниже). В этой постановке, средний градиент давления равен нулю, что обусловливает постоянное, однородное по вертикали, значение полного потока импульса воздушного потока [Druzhinin et al., 2017].

Рис. 1. Схема численного эксперимента (см. текст).

Модель, рассматриваемая в настоящей работе, основана на решении трехмерных, отфильтрованных по подсеточным флуктуациям, уравнений движения воздушной (несущей) фазы в эйлеровой формулировке и лагранжевых уравнений движения отдельных капель. Рассматриваются достаточно мелкие капли (с диаметром до 300 микрон), что позволяет пренебречь их деформацией. В качестве нижней границы рассматривается стационарная поверхностная волна, задаваемая неявно с помощью отображения координат (обсуждаемого ниже). Для расчета подсеточных турбулентных напряжений в уравнении для несущей фазы используется концепция турбулентной вязкости, определяемой кинетической энергией подсеточных пульсаций, для которой решается отдельное прогностическое уравнение, и пространственным масштабом, задаваемым локальным размером ячеек расчетной сетки. Влияние капель на несущий воздушный поток не учитывается, поверхностная волна считается заданной; не рассматриваются эффекты, связанные с испарением капель. Верификация модели осуществляются сравнением фазово-осредненных профилей скорости и потока импульса несущей фазы и концентрации капель с результатами прямого численного моделирования полных уравнений динамики несущей и дисперсной фаз для различных сценариев инжекции в поток и размеров капель. Ниже в разделе 2 представлены уравнения движения несущей и дисперсной фаз и способы их численного решения; в разделе 3 верификация и калибровка модели на основе сравнения с результатами прямого численного моделирования [Druzhinin, et al., 2017] при числе Рейнольдса воздушного потока Re = 1.5 × 10⁴. Далее в разделе 4 обсуждаются результаты вихреразрешающего моделирования двухфазного воздушного потока над взволнованной водной поверхностью с числом Re = 10⁵.

2. Математическая модель

Постановка задачи схематично представлена на рис. 1 и в целом эквивалентна ранее использовавшейся постановке расчетов прямого численного моделирования (DNS) в [Druzhinin et al., 2017]. Рассматривается воздушный поток над взволнованной водной поверхностью, несущий капли. Верхняя граница области счета плоская и движется в направлении распространения поверхностной волны со скоростью U₀, связанной со средней, характерной скоростью ветра в лабораторных или натурных условиях. Нижняя граница совпадает с водной поверхностью, на которой распространяется двумерная, стационарная волна с амплитудой a, длиной λ и фазовой скоростью c; крутизна волны ka = 2πa/λ задается равной ka = 0.1. Значения U₀ и λ определяют масштабы скорости и длины задачи, на которые далее нормируются все переменные.

Фильтрованные уравнения Навье-Стокса для воздуха решаются численно в безразмерных переменных в виде:

$$\begin{gathered} \frac{\partial U_i}{\partial t}+\frac{\partial \left(U_i{U}_j\right)}{\partial x_j}=-\frac{\partial P}{\partial x_j}+\frac{1}{Re}\frac{{\partial }^2{U}_i}{\partial x_j\partial x_j}-\frac{\partial {\tau }_{ij}}{\partial x_j}, \\ \frac{\partial U_j}{\partial x_j}=0. \end{gathered}$$ (1)

В уравнениях (2) и далее по повторяющимся индексам проводится суммирование, $x_i=x,y,z(i=\mathrm{1,}\mathrm{2,}3);U_i^{}$, и P — фильтрованные поля скорости и давления; число Рейнольдса определяется масштабами длины и скорости,

$Re=\frac{U_0\lambda }{v}$ , (2)

где ν — кинематическая вязкость воздуха (≈ 0.15 см²/сек). τ_ij — тензор подсеточных напряжений поля скорости воздуха, для которого используется параметризация в виде [Sullivan, et al., 2008]:

$$\begin{gathered} {\tau }_{ij}=-v_t\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right),v_t=c_k{e}^{1/2}l \\ l={\left(\Delta x_1\Delta x_2\Delta x_3\right)}^{1/3} \end{gathered}$$ (3)

В выражении (3) $v_t$ — турбулентная вязкость, определяемая полем кинетической энергии подсеточных пульсаций e и пространственным масштабом l, задаваемым локальным размером ячеек расчетной сетки. Поле e вычисляется интегрированием отдельного прогностического уравнения [Sullivan et al., 2008]:

$$\begin{gathered} \frac{\partial e}{\partial t}+\frac{\partial \left(U_{j}e\right)}{\partial x_j}=-\frac{{\tau }_{ij}}{2}\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)+ \\ +2\frac{\partial }{\partial x_j}\left(v_t\frac{\partial e}{\partial x_j}\right)-c_D\frac{e^{3/2}}{l} \end{gathered}$$ (4)

В уравнения (3) и (4) входят модельные константы $c_k=0.1$ и $c_D=0.93$. Значения этих констант подобраны в предположении о наличии протяженной части инерционного интервала на подсеточном масштабе. Применимость этого предположения для рассматриваемой конфигурации потока верифицируется сравнением результатов расчетов LES с результатами прямого численного моделирования при Re = 1.5 × 10⁴ (см. раздел 3 ниже).

Одновременно с уравнениями (1) и (4) для эйлеровых полей несущей фазы решаются уравнения для лагранжевых скоростей и координат капель [Druzhinin et al., 2017]:

$$\begin{gathered} \frac{d{r}_i^n}{dt}=V_i^n \\ \frac{d{V}_i^n}{dt}=\frac{f\left({Re}_n\right)}{{\tau }_d^n}\left(U_i^n-V_i^n\right)-g^{\text{'}}{\delta }_{iz} \end{gathered}$$ (5)

В уравнениях (5) $r_i^n,V_i^n(i=x,y,z)$ — координата и скорость n–й капли, $U_i^n$ — локальная, мгновенная скорость воздуха; производная по времени вычисляется вдоль траектории капли; $g^{\text{'}}$ — безразмерное, приведенное ускорение силы тяжести (обсуждаемое ниже). Поскольку отношение плотности воздуха к плотности воды мало (ρ_a /ρ_w порядка 10^–3), в правой части уравнения (5) фигурируют только силы вязкого трения и тяжести, существенно превалирующие над остальными силами (градиента давления, Бассе, Сэфмана и присоединенной массы). Сила вязкого трения включает фактор

$f\left({\mathrm{Re}}_n\right)=1+0.15{\mathrm{Re}}_n^{0.687}$, (6)

учитывающий конечность числа Рейнольдса капли, равного

 ${Re}_n=\frac{d_n\left|U^n-V^n\right|}{v}$. (7)

В силу трения входит также безразмерное время релаксации капли ${\tau }_d^n$,

${\tau }_d^n=\frac{d_n^2}{18v}\frac{{\rho }_w}{{\rho }_a}\frac{U_0}{\lambda }$ . (8)

В локальную скорость воздуха $U_i^n$ в правой части уравнения (5) входят как фильтрованная скорость воздуха $U_i\left(r_{}^n,t\right)$, так и вклад подсеточных пульсаций скорости воздуха, окружающего каплю, $u_i^n\left(t\right)$,

$U_i^n=U_i\left(r_{}^n,t\right)+u_i^n\left(t\right)$. (9)

В общем случае вклад $u_i^n\left(t\right)$ определяется решением стохастического уравнения в виде [Pozorski, Apte, 2009; Michałek, et al., 2013]:

$\frac{d{u}_i^n}{dt}=-\frac{u_i^n}{T_L}+C{\left(\frac{e}{T_L}\right)}^{1/2}{w}_i+$

$+{\delta }_{iz}\exp \left(-\frac{{\tau }_d{e}^{1/2}}{l}\right)\frac{\partial e}{\partial \eta }\left[1+\frac{{\left(u_i^n\right)}^2}{e}\right],$ (10)

где δ_ij — символ Кронекера,  $w_i$— дельта-коррелированный, гауссовский случайный процесс с единичной дисперсией, а время корреляции подсеточных пульсаций вдоль лагранжевой траектории капли определяется как

$T_L=\frac{l}{e^{1/2}}{\left[1+\frac{3W_s^2}{2e}\right]}^{-1/2}$. (11)

В уравнение (10) входят поправки, обусловленные инерцией капель и зависимостью распределения кинетической энергии подсеточных пульсаций от расстояния (η) от местоположения капли до водной поверхности. Физический смысл аппроксимации (11) для времени корреляции состоит в том, что время корреляции пульсаций скорости воздуха вдоль траектории капли уменьшается с увеличением скорости ее оседания $W_s$ благодаря эффекту «пересечения траекторий» [Wells, Stock, 1983]; $W_s$ находится численным решением методом Ньютона уравнения для установившегося движения капли (5) в покоящемся воздухе:

$W_{s}f\left(\frac{W_{s}d}{v}\right)={\tau }_d{g}^{\text{'}}$. (12)

Последнее слагаемое в правой части выражения (10) учитывает эффекты, связанные с неоднородностью распределения кинетической энергии подсеточных пульсаций воздуха с расстоянием до водной поверхности [Michałek, et al., 2013].

В случае однородной турбулентности в пределе больших значений Re значение константы С в правой части уравнения (10) связано с значением константы Колмогорова структурной функции скорости в инерционном интервале спектра. В рассматриваемом случае турбулентного воздушного потока над взволнованной поверхностью с $\text{Re = O}\left({\text{10}}^{\text{4}}{\text{-10}}^{\text{5}}\right)$ значение С ≈ 2.3 определяется согласием результатов предсказаний модели для фазово-осредненного профиля концентрации капель с результатами прямого численного интегрирования полных уравнений динамики воздушной и дисперсной фаз для различных сценариев инжекции капель в поток (см. ниже обсуждение в разделе 3 и рис. 4, 5).

Следует отметить, что, как показывают результаты вихреразрешающего моделирования двухфазных турбулентных потоков в каналах, вклад подсеточных пульсаций в динамику частиц может существенным образом влиять на распределение поля их концентрации (см. напр. [Michałek, et al., 2013]). Однако в [Richter, et al., 2019] вклад подсеточных пульсаций не учитывался, поскольку в этих расчетах их кинетическая энергия быстро спадала с расстоянием до водной поверхности и оказывалась пренебрежимо малой по сравнению с кинетической энергией пульсаций фильтрованного поля скорости. В настоящей работе вклад $u_i^n\left(t\right)$ не отбрасывается с целью проверки обоснованности предположения о его малости путем сравнения результатов LES с результатами прямого численного моделирования.

Расчеты проводятся в системе отсчета, движущейся с фазовой скоростью волны. При этом используется конформное преобразование координат декартовой системы в вертикальной плоскости, x(ξ,η) и z(ξ,η), в виде [Druzhinin, et al., 2017]:

$$\begin{gathered} x=\xi -a\exp (-k\eta )\sin k\xi \\ z=\eta +a\exp (-k\eta )\cos k\xi \end{gathered}$$ (13)

где безразмерное волновое число k = 2π. Преобразование (13) переводит область с волнообразной нижней границей в область, ограниченную нижней плоской границей при η = 0. Форма водной поверхности $z_s(x,t)$ задается неявно преобразованием (13) и с точностью до членов второго порядка малости по $ka$ совпадает с решением, описывающим волну Стокса [Gent, Taylor, 1976]. Используется также дополнительное преобразование вертикальной координаты, обеспечивающее сгущение узлов расчетной сетки по вертикали и разрешение больших градиентов в пограничных слоях, формирующихся вблизи верхней и нижней границ области счета.

Следует отметить, что сгущение узлов сетки у верхней и нижней границ позволяет разрешать с достаточной точностью поля обтекания в окрестности поверхностных волн, а также течение в вязком подслое. Такое разрешение в особенности необходимо для применимости в расчетах условий инжекции капель, которая, согласно данным недавних лабораторных экспериментов [Troitskaya et al., 2017], происходит вблизи водной поверхности благодаря развитию гидродинамической неустойчивости возмущений поля скорости в приповерхностном водном слое.

Уравнения (1) решаются численно с использованием метода Адамса–Бэшфорта второго порядка точности в области 0 ≤ ξ ≤ 6, –2 ≤ y ≤ 2, 0 ≤ η ≤1 на разнесенной сетке. Используется метод расщепления, и уравнение Пуассона для давления решается методом итераций с помощью преобразования Фурье в горизонтальной плоскости и методом Гаусса по вертикальной координате. Локальная скорость воздуха в местоположении каждой капли определяется с помощью Эрмитовой интерполяции [Fletcher, 1991; Белоцерковский, 1984]. Все поля являются периодическими по горизонтальным координатам. На верхней и нижней (движущихся) границах ставятся условия Дирихле для полей скорости воздуха и нулевое условие для кинетической энергии подсеточных пульсаций. На нижней границе задается двумерное поле скорости, определяемое потенциальным решением для волны Стокса на глубокой воде второго порядка точности по крутизне волны [Филлипс, 1980].

При инициализации расчетов задается поле скорости слабо возмущенного ламинарного течения Куэтта (с амплитудой флуктуаций порядка 0.1 %). По истечении переходного процесса турбулентное течение становится статистически стационарным, и в воздушный поток случайным образом инжектируются капли с изначально равномерным пространственным распределением, диаметром в диапазоне от 100 до 300 µm, и со скоростью, равной мгновенной скорости окружающего воздуха. Капли, достигающие верхней границы или падающие на водную поверхность, ре-инжектируются, так что полное число капель поддерживается постоянным (равным 10⁶) на протяжении всего времени счета. Рассматриваемое число капель достаточно велико, и погрешности, связанные с усреднением и определением статистических характеристик дисперсной фазы, пренебрежимо малы.

Поскольку известные данные натурных и лабораторных наблюдений не позволяют определить параметры инжекции капель в воздушный поток, в математических моделях рассматриваются различные сценарии [Peng, Richter, 2019; Druzhinin et al., 2017; Richter et al., 2019; Druzhinin 2021]. В настоящей работе рассматриваются сценарии, когда капли инжектируются либо вблизи гребней волн, либо с однородным пространственным распределением по фазе волны. Частицы инжектируются на границе вязкого подслоя и в буферной области пограничного слоя со скоростями, равными либо скоростям соответствующих Лагранжевых частиц в поверхностной волне, либо локальной скорости воздуха.

Известные данные экспериментов и наблюдений говорят о том, что для капель с диаметром около 200 µm отношение произведения скорости оседания капли в покоящейся среде (равной произведению τ_dg ׳, где τ_d определяется в (8)) к произведению постоянной Кармана и скорости трения, Ƙu_*, порядка единицы. Соответственно, приведенное значение ускорения гравитации g ׳ определяется таким образом, что это условие выполняется для рассматриваемых параметров воздушного потока.

По достижении стационарного режима течения рассчитываются статистические характеристики поля течения воздушного потока и динамики капель. Фазово-осреднённые поля находятся с помощью осреднения по поперечной координате, времени и «оконным» (определяемым длиной волны) осреднением по продольной координате. Вертикальные профили получаются дополнительным осреднением фазово-осреднённых полей по продольной координате.

3. Верификация модели

С целью верификации вихреразрешающей (LES) модели, представленной выше, проводились расчеты для числа Рейнольдса воздушного потока (3) равного Re = 1.5 × 10⁴ на сетке, состоящей из 180 × 120 × 90 узлов. С одной стороны, рассматриваемое Re достаточно велико (соответствующее число Рейнольдса, основанное на скорости трения u_* ≈ 0.03U₀, равно Re_* = u_*λ/ν ≈ 500) и позволяет достичь воздушному потоку состояния полностью развитой турбулентности. С другой стороны, в этом случае результаты можно сравнить с соответствующими результатами прямого численного (DNS) моделирования, проводившегося для тех же конфигурации и параметров течения, но с использованием сетки из 360 × 240 × 180 узлов [Druzhinin, et al., 2017]. Во всех рассматриваемых случаях задавалась крутизна поверхностной волны $ka$ = 0.1. Расчеты проводились для капель с диаметром в диапазоне от 100 до 300 µm, как для сценариев однородной, так и неоднородной инжекции капель со скоростью, равной скорости лагранжевых частиц в поверхностной волне. Кроме того, чтобы оценить вклад в динамику капель подсеточных флуктуаций, проводились также расчеты и без использования стохастического форсинга. Сравнение расчетов DNS и LES представлено на рис. 2–5.

Рис. 2. Профили средней скорости U_x (a) и ее среднеквадратичных флуктуаций  $U_i^{\text{'}}(i=x,y,z, \mathrm{b}-\mathrm{d})$. Результаты расчетов LES и DNS черным и синим цветом; η – расстояние от водной поверхности. Асимптотика (14) — коротким штрихом (а); среднеквадратичная скорость подсеточных флуктуаций — длинным штрихом (b,c,d).

Рис. 3. Распределение фазово-осредненных полей горизонтальной (a,b) и вертикальной (c,d) компонент скорости, давления (e,f) и потока импульса (g,h) в расчетах DNS (a,c,e,g) и LES (b,d,f,h) в приповерхностном воздушном слое (z < 0.1). Инкремент изолиний 0.02 (a,b), 0.001 (c,d), 0.2 (e,f), и 5 × 10^–5 (g,h). Отрицательные значения показаны штрихом.

Рис. 4. Распределение фазово-осредненных полей концентрации капель в случаях однородной (a,b) и неоднородной (c,d) инжекции в расчетах DNS (a,c) и LES (b,d). Инкремент изолиний 4 × 10^–5 (a,b) и 10^–4 (c,d).

Рис. 5. Профили средней концентрации капель C в расчетах DNS (длинным штрихом) и LES для случаев неоднородной (a) и однородной (b) инжекции. Коротким штрихом обозначены профили в LES без применения стохастического форсинга.

На рис. 2 показаны профили скорости воздушного потока $U_x$(a) и компонент её среднеквадратичных флуктуаций ${U^{\text{'}}}_i^{}$ ($i=x,y,z,\text{ b, c, d}$), полученные в расчетах LES и DNS осреднением по поперечной и продольной координатам y и ξ и времени. Профили средней скорости также сравниваются с асимптотикой

$U_{\log }=\frac{1}{\kappa }\ln \frac{\eta }{{\eta }_0}$ (14)

где Ƙ = 0.4 – постоянная Кармана; η₀ — параметр шероховатости водной поверхности, определяемый как вязкостью воздуха, так и амплитудой волн [Zeng et al., 1998]. В рассматриваемом случае этот параметр задается в виде [Druzhinin, 2021]

${\eta }_0=\frac{0.11v}{u_{*}}+0.025a$. (15)

Сравнение результатов расчетов LES и DNS на рис. 2 показывает хорошее согласие для вертикальных профилей средней скорости (рис. 1а, различие порядка 5 %), и неплохое согласие для профилей флуктуаций скорости (рис. 1b,c,d, различие порядка 10–15 %). Из рисунка видно также, что амплитуда подсеточных пульсаций сравнима с амплитудой флуктуаций фильтрованной скорости вблизи водной поверхности (при η ≈ 0.01), и становится пренебрежимо малой на достаточном удалении от поверхности (при η > 0.05). Наблюдается также неплохое согласие профилей фазово-осредненных горизонтальной и вертикальной компонент скорости, давления и потока импульса

null, (16)

представленных на рис. 3.

Как видно из рис. 2 (b,c,d), кинетическая энергия подсеточных флуктуаций скорости мала по сравнению с амплитудой фильтрованной турбулентности на достаточном удалении от водной поверхности (при η > 0.1). Профиль средней скорости здесь следует логарифмической асимптотике (14). Этот же комментарий можно отнести и к случаю с $\text{Re }={10}^5$ (см. рис. 6 ниже). Таким образом, подсеточный поток импульса в этой области мал по сравнению с потом импульса фильтрованной компоненты турбулентности. С другой стороны, рисунок 3 (g,h) показывает, что в приповерхностном воздушном слое (η < 0.1) распределение фазово-осредненного потока импульса фильтрованной турбулентности в LES довольно хорошо согласуется с турбулентным потоком импульса в расчетах DNS. Небольшие отличия могут быть связаны с влиянием подсеточной компоненты турбулентности в области η < 0.01, где ее амплитуда сравнима с амплитудой фильтрованной компоненты. По-видимому, согласие расчетов LES и DNS на рис. 2 и 3 может быть улучшено благодаря коррекции замыкания (3) и модельных констант для рассматриваемой конфигурации потока при значениях Re = O(10⁴ – 10⁵). Этот вопрос требует дополнительного исследования, выходящего за рамки настоящей работы.

Рис. 6. Профили средней скорости U_x (a), кинетической энергии пульсаций фильтрованного поля скорости и подсеточных пульсаций (b), и концентрации капель (c) в расчетах LES с числом Рейнольдса Re = 10⁵. Расстояние от водной поверхности, равное высоте вязкого подслоя (η = 5η_*) и амплитуде волны (η = a) отмечено коротким штрихом (а). Штрихом отмечена асимптотика (14) (а), кинетическая энергия подсеточных пульсаций (b), и концентрация капель в LES без использования стохастического форсинга (с).

На рис. 4 и 5 представлены фазово-осредненный и средний вертикальный профили поля концентрации капель для различных сценариев их инжекции в поток. На рис. 5 также приводятся профили концентрации в расчетах LES без использования стохастического форсинга. Мгновенное поле концентрации определяется согласно

$C\left(r\right)=\frac{\pi }{6}\sum _{n=1}^{N_d}\frac{d_n^{3}w\left(r^n,r\right)}{{\Omega }_g^n},$ (17)

где $w(r^n,r)$ — геометрический весовой коэффициент, обратно пропорциональный расстоянию между n-й каплей с координатами $r^n=(x^n,y^n,z^n)$ и узлом сетки с координатами $r=(x,y,z);{\Omega }_g^n$  – объем рассматриваемой ячейки сетки. Для каждой капли коэффициенты $w(r^n,r)$ (их всего восемь – по числу ближайших, окружающих каплю, узлов ячейки) нормируются таким образом, что их сумма равна единице. Таким образом, значения концентрации капель экстраполируются на узлы эйлеровой расчетной сетки. Расчеты фазово-осредненного и среднего профилей осуществляются аналогично этой процедуре для полей скорости.

Как видно из рис. 4, результаты расчетов LES и DNS для концентрации капель довольно хорошо согласуются во всех рассмотренных случаях. При этом, как показывает рис. 5, применение форсинга (для учета подсеточных флуктуаций в динамике капель) незначительно улучшает согласие профилей концентрации в случае неоднородной инжекции (~ на 10%) и практически не влияет на профили C в случае однородной инжекции.

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

Результаты вихреразрешающего моделирования при Re = 105

Модель LES, обсуждаемая выше, была задействована для расчетов двухфазного потока над водной поверхностью с числом Рейнольдса Re = 10⁵, ${\mathrm{Re}}_{*}\approx 3000.$ Расчеты проводились с использованием сетки из 360 × 240 × 180 узлов. Остальные параметры модели задавались такими же, как и в расчетах LES при Re = 1.5 × 10⁴ (Re ≈ 500). Результаты представлены на рис. 6 и 7.

Рис. 7. Распределение фазово-осредненных полей горизонтальной и вертикальной компонент скорости (a,b), давления (c) и потока импульса (d) воздушного потока, и концентрации (e,f) и силы воздействия капель на воздух (g,h) в расчетах LES с Re = 10⁵ в приповерхностном воздушном слое (z < 0.1). Отрицательные значения показаны штрихом.

На рис. 6 приводятся профили средних полей скорости воздушного потока (6a), кинетической энергии фильтрованного поля скорости воздуха и подсеточных пульсаций (6b), и концентрации капель (как с учетом, так и без учета стохастического форсинга) (6c). Как видно из рис. 6а, сгущение узлов сетки у водной поверхности обеспечивает необходимое разрешение поля скорости воздуха как области вблизи границы вязкого подслоя (η ≈ 5η_* ≈ 1.8 × 10^–3), так и в буферном слое и на расстоянии порядка амплитуды волны (в слое η ≈ a ≈ 0.016), позволяя таким образом отслеживать форму водной поверхности. На достаточном удалении от поверхности профиль скорости следует логарифмическому закону и хорошо описывается асимптотикой (14). Из рис. 6b также видно, что в области вне вязкого подслоя (η > 0.02 ) кинетическая энергия турбулентности фильтрованного поля скорости E растет и значительно превышает значение энергии подсеточных пульсаций e, что обусловливает слабое влияние стохастического форсинга на распределение концентрации капель (рис. 6c).

Фазово-осреднённые поля скорости, давления и потока импульса воздушного потока (на рис. 7a–d) качественно мало отличаются от соответствующих распределений полей в расчетах при Re = 1.5 × 10⁴, что дает основание полагать, что поле средней скорости при Re = 10⁵ по-прежнему остается безотрывным [Druzhinin, et al., 2017]. На рис. 7e,f и 7g,h показаны профили концентрации капель (17) и силы $F\left(r\right),$ с которой капли воздействуют на воздушный поток при учете их взаимодействия с воздухом (обозначения те же, что и в (17)), определяемой согласно

$$\begin{gathered} F\left(r\right)=3\pi \tau \sum _{n=1}^{N_d}{d}_n\left(V_x^n-U_x^n\right)f\left({Re}_n\right)\frac{w\left(r^n,r\right)}{{\Omega }_g^n}= \\ =\sum _{n=1}^{N_d}\frac{\pi d_n^3}{6}\frac{{\rho }_w}{{\rho }_a}\frac{1}{{\tau }_d^n}\left(V_x^n-U_x^n\right)f\left({Re}_n\right)\frac{w\left(r^n,r\right)}{{\Omega }_g^n} \end{gathered}$$ (18)

Расчеты профилей проводились для двух различных сценариев ре-инжекции падающих в воду капель — со скоростями частиц в поверхностной волне (рис. 7e,g) и со скоростями, равными локальной скорости воздуха (рис. 7f,h). В обоих случаях инжекция производилась в окрестности гребней волн. Рисунки показывают, что сценарий инжекции мало влияет на профиль концентрации, но коренным образом, качественно влияет на характер воздействия частиц на поток: капли, инжектируемые со скоростью поверхностной волны, в среднем тормозят поток на наветренном склоне волны, т.к. в этой области ${⟨F⟩}_x<0$, в то время как капли, инжектируемые со скоростью воздуха, ускоряют поток на подветренном сколе, т.к. здесь ${⟨F⟩}_x>0$ . Этот результат также согласуется с выводами по расчетам DNS, проводившимся в [Druzhinin et al., 2017] при Re = 1.5 × 10⁴.

5. Выводы

Представлена математическая модель и проведено вихреразрешающее моделирование динамики турбулентного воздушного потока над взволнованной водной поверхностью, несущего мелкодисперсные капли. Математическая модель основана на решении трехмерных, отфильтрованных по подсеточным флуктуациям, уравнений движения воздушной (несущей) фазы в эйлеровой формулировке и уравнений движения отдельных капель в лагранжевой формулировке. Для замыкания подсеточных напряжений поля скорости воздуха использована концепция турбулентной вязкости, где кинетическая энергия неразрешаемых расчетной сеткой пульсаций определяется решением прогностического уравнения. Рассмотрены достаточно мелкие капли (с диаметром до 300 микрон), что позволяет пренебречь их деформацией, не учитываются также столкновения капель между собой, их теплообмен с воздухом и испарение. Массовая концентрация капель задается достаточно малой, что позволяет пренебречь их воздействием на поток; поверхностная волна считается заданной, неизменяющейся под действием капель и(или) ветра. Проведена верификация модели сравнением с результатами прямого численного интегрирования полных уравнений динамики воздушной и дисперсной фаз для различных сценариев инжекции в поток и размеров капель. Получены фазово-осредненные профили скорости, потока импульса воздуха и концентрации капель для значений числа Рейнольдса несущего потока до 10⁵, которое на порядок превышает типичные значения Re, задаваемые при проведении прямого численного моделирования.

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

Источник финансирования

Работа выполнена при поддержке проекта РНФ № 24-27-00216.

作者简介

O. Druzhinin

编辑信件的主要联系方式.
Email: druzhinin@ipfran.ru
俄罗斯联邦, Nizhny Novgorod ул. Ульянова 46, 603950

参考

补充文件

附件文件

动作

1. JATS XML

下载

2. Fig. 1. Schematic diagram of the numerical experiment (see text).

下载 (17KB)

索引源数据

3. Fig. 2. Profiles of the mean velocity (a) and its root-mean-square fluctuations ( ). Results of LES and DNS calculations in black and blue; η is the distance from the water surface. Asymptotics (14) is shown by a short dash (a); the root-mean-square velocity of subgrid fluctuations is shown by a long dash (b,c,d).

下载 (36KB)

索引源数据

4. Fig. 3. Distribution of phase-averaged fields of the horizontal (a,b) and vertical (c,d) components of velocity, pressure (e,f), and momentum flux (g,h) in DNS (a,c,e,g) and LES (b,d,f,h) calculations in the near-surface air layer (z < 0.1). The increment of the isolines is 0.02 (a,b), 0.001 (c,d), 0.2 (e,f), and 5 × 10–5 (g,h). Negative values are shown with a dash.

下载 (107KB)

索引源数据

5. Fig. 4. Distribution of phase-averaged droplet concentration fields for uniform (a,b) and non-uniform (c,d) injection cases in DNS (a,c) and LES (b,d) calculations. The increment of the isolines is 4 × 10–5 (a,b) and 10–4 (c,d).

下载 (37KB)

索引源数据

6. Fig. 5. Profiles of the average droplet concentration C in DNS (long dash) and LES calculations for non-uniform (a) and uniform (b) injection cases. Short dashes indicate profiles in LES without stochastic forcing.

下载 (15KB)

索引源数据

7. Fig. 6. Profiles of the mean velocity (a), kinetic energy of the filtered velocity field pulsations and subgrid pulsations (b), and droplet concentration (c) in LES calculations with the Reynolds number Re = 105. The distance from the water surface equal to the height of the viscous sublayer (η = 5η*) and the wave amplitude (η = a) is marked with a short dash (a). The dash indicates the asymptotics (14) (a), kinetic energy of the subgrid pulsations (b), and droplet concentration in LES without using stochastic forcing (c).

下载 (22KB)

索引源数据

8. Fig. 7. Distribution of phase-averaged fields of the horizontal and vertical components of velocity (a,b), pressure (c), and momentum flux (d) of the air flow, and concentration (e,f) and the force of droplet impact on air (g,h) in LES calculations with Re = 105 in the near-surface air layer (z < 0.1). Negative values are shown as dashes.

下载 (79KB)

索引源数据