Прикладная математика и механика, 2022, T. 86, № 5, стр. 765-778

1. Введение. Во второй половине XX века было обнаружено, что большая часть турбулентного атмосферного пограничного слоя (АПС) занята роллами – упорядоченными спиралевидными вихрями (валами) с горизонтальной осью, по направлению близкой к среднему направлению геострофического ветра (см., например, [1]). На спутниковых снимках пограничного слоя роллы выглядят, как “облачные улицы” – вытянувшиеся на сотни километров параллельные ряды с периодом в несколько километров, сохраняющиеся в течение нескольких суток. Облака в таком случае формируются в области восходящих движений между роллами при соответствующих термодинамических условиях (см. рис. 1). Возникновение роллов, как правило, является следствием развития гидродинамических неустойчивостей в АПС [1, 2], в том числе и конвективной неустойчивости при умеренных ветрах (2–3.5 м/с), развивающейся как в умеренных и южных широтах, так и в высоких при холодных вторжениях [3–7].

В безоблачных условиях спиралевидные вихри также хорошо наблюдаются средствами дистанционного зондирования в других диапазонах волн [3, 4, 8, 9]. На рис. 2 представлено пространственно-периодическое распределение аэрозоля (мелкая фракция), зафиксированное с помощью лидарных измерений коэффициента обратного рассеяния на пыли в пустынях Калмыкии [10], где из-за низкой влажности образование облачности в области локализации вихревых структур не происходит. Пространственный масштаб квазипериодичности зафиксированных аэрозольных слоев составляет 2–5 км, что характерно для валиковой циркуляции, высота роллов – не более 2 км [11]. Периодичность в распределении аэрозоля отмечается на всех высотах лидарного зондирования – от 300 м до 1000 м.

Согласно оценкам, вклад мезомасштабных роллов в тепломассоперенос через АПС составляет от 20 до 60%, причем вклад в вертикальный перенос в высоких широтах превалирует над турбулентным [6, 7, 12]. Роллы играют важную роль в индуцировании формирования новых частиц и, соответственно, образования ядер облачной конденсации [13]. В океаническом пограничном слое наблюдаются схожие движения (ленгмюровская циркуляция) [14]. Совокупность этих факторов и того, что оценки мощности антропогенных выбросов аэрозолей в атмосферу составляют порядка 10⁷–10⁸ т. в год [15, 16], делают задачу моделирования распределения концентраций мелкодисперсных аэрозолей в роллах существенной для моделирования процессов переноса и химической трансформации антропогенных примесей в АПС в целом.

Особый интерес представляет задача переноса аэрозолей в случае, когда тяжелый аэрозоль вступает в реакцию с внешней средой или распадается. Распространяясь в атмосфере, он диффундирует и под действием силы тяжести опускается на землю с постоянной скоростью, которая предвычисляется из задачи Стокса [17, 18].

В настоящей работе предлагаются два обоснованных метода расчета полей концентраций мелкодисперсных аэрозолей в вихревых структурах, использующих результаты моделирования поля скоростей в вихрях [19]. На основе методов теории возмущений [20] получено приближенное решение стационарной пространственно-периодической сингулярно возмущенной задачи реакция–диффузия–адвекция, моделирующей распределений аэрозоля, оценен остаточный член и предложен метод численного решения задачи нулевого приближения. В качестве альтернативного подхода к задаче численного моделирования поля концентраций рассмотрена реализация метода эволюционной факторизации, на основе которого исследован вопрос о влиянии величин скоростей распада и оседания аэрозоля на распределение концентрации этой примеси в роллах. Использование этих двух подходов строго обосновано с привлечением методов и результатов [20–23]. Таким образом, описано единственное устойчивое стационарное состояние системы путем описания распределения концентрации аэрозоля, отвечающего этому состоянию. С использованием модельных данных получена оценка количества переносимого аэрозоля.

В заключение этого параграфа отметим другие примеры эффективного использования методов асимптотического анализа при решении прямых и обратных задачах атмосферной диффузии малых примесей [24–26], связанных с получением оценок эмиссий локальных и распределенных источников загрязняющих примесей и параметров турбулентной диффузии для устойчивого или нейтрального АПС, а также при решении обратных задач по восстановлению параметров моделей реакция–диффузия–адвекция других приложений [27–32].

2. Постановка задачи. Изменение концентрации $C({{x}_{1}},{{x}_{2}},{{x}_{3}},t)$ аэрозоля при миграции через АПС описывается уравнением типа реакция–диффузия–адвекция [17]:

где $\left\langle k \right\rangle $ – среднее значение горизонтального коэффициента турбулентной диффузии, ${{k}_{3}}(x)$ – вертикальный коэффициент турбулентной диффузии [17], ${{u}_{i}}(x)$, $i = \overline {1,3} $ – составляющие средней скорости переноса примеси вдоль соответствующих осей, ${{{v}}_{g}}$ – абсолютная величина скорости оседания частиц примеси, $F(C,x)$ – мощность стока вещества.

В качестве масштабов длины, скорости и времени выберем толщину экмановского слоя $L = {{\left( {K{\text{/}}f} \right)}^{{1/2}}}\sim 1000$ м, скорость геострофического ветра $G$ ~ 10 м/с [2] и $1{\text{/}}f$, где $f = \Omega \sin \theta $, Ω – угловая скорость вращения Земли, $\theta - $ широта, K – характерная вертикальная турбулентная вязкость (подробнее см. [1, 2, 33, 34]). Поскольку характерные горизонтальные масштабы развивающихся структур кратны высоте экмановского слоя, то в качестве характерного пространственного масштаба по трем измерениям используется высота слоя Экмана (см., например, [2, 11, 19]).

Заметим, что детальный анализ данных позволяет говорить о сосуществовании в АПС наряду со структурами, горизонтальный масштаб $\lambda $ которых в несколько раз превосходит вертикальный $h$, где ${\lambda \mathord{\left/ {\vphantom {\lambda h}} \right. \kern-0em} h} \geqslant 2$, структур, вертикальный масштаб которых равен и меньше горизонтального ${\lambda \mathord{\left/ {\vphantom {\lambda h}} \right. \kern-0em} h} \leqslant 1$ [6, 35].

Определим безразмерные переменные посредством соотношений: $x{\kern 1pt} ' = \frac{x}{L}$, $t{\kern 1pt} ' = \frac{t}{T}$, $C{\kern 1pt} ' = \frac{C}{{\tilde {C}}}$, где T – характерное время существования валиковых структур (несколько суток), $\tilde {C}$ – характерная концентрация, например, фоновое значение концентрации (или ее среднее на некоторой высоте). Тогда в безразмерных переменных уравнение (2.1) принимает вид:

Согласно [18] на высоте более 100 м над поверхностью Земли изменчивостью коэффициента обмена можно пренебречь. Следовательно

где Pr_D – турбулентное диффузионное число Прандтля, Re – число Рейнольдса. В левой части уравнения (2.2) при производной от концентрации по времени в качестве множителя содержится отношение характерного временного масштаба $L{{G}^{{ - 1}}}$ к характерному времени существования валиковых структур. Этот масштаб соответствует времени установления экмановского профиля при изменении внешних условий (напр., при изменении направления и скорости ветра), а также характерному времени развития субмезомасштабных структур вследствие сдвиговых неустойчивостей [36], обеспечивающих основной первичный вынос аэрозоля с поверхности. Поскольку $T \gg L{{G}^{{ - 1}}}$, то влиянием нестационарности можно пренебречь и рассматривать процесс как стационарный. Сохранив прежние обозначения для безразмерных концентраций, скорости переноса и оседания примеси приходим к уравнению в безразмерных переменных:

Заметим, что валиковыми структурами переносятся в основном мелкодисперсные аэрозоли, для которых скорость оседания ∼1 см/с [34]. В таком случае безразмерная скорость оседания ${{{v}}_{g}}$ ~ ${{10}^{{ - 3}}}$.

Пусть в уравнении (2.1) $F(C,x) = - \gamma (C - \tilde {C})$, где $\gamma > 0$ – скорость распада аэрозоля (величина, обратно пропорциональная времени жизни примеси), а слагаемое $F(C,x)$ описывает распад вещества за счет столкновения частиц при миграции аэрозоля. Тогда уравнение (2.3) примет следующий вид:

где под обозначением ${\mathbf{\bar {u}}}(x)$ понимается вектор с компонентами ${{u}_{1}}(x)$, ${{u}_{2}}(x)$, ${{u}_{3}}(x) - {{{v}}_{g}}$, ${{{v}}_{g}}$ и $\gamma : = \gamma L{{G}^{{ - 1}}} - $ безразмерные параметры модели.

Поскольку поле скоростей ${\mathbf{\bar {u}}}(x)$ в уравнении (2.4) предполагается гладким, то решение краевой задачи для уравнения в ограниченной области с гладкой границей, включающей в себя гряду роллов, и граничным условием смешанного типа [17], которое соответствует заданию концентрации или нулевого потока вещества на различных частях границы, существует и единственно [21, 22], а также глобально устойчиво по Ляпунову [23] как стационарное решение соответствующей начально-краевой задачи для уравнения:

где $L{{(TG)}^{{ - 1}}}\sim {{({{\Pr }_{D}}\operatorname{Re} )}^{{ - 1}}}$. Эти свойства стационарного решения ниже будут использованы для обоснования численного алгоритма расчета распределения концентрации аэрозоля в вихревых структурах.

С целью выбора рациональных средств численного моделирования распределения аэрозоля внутри вихревых структур, перейдем к эквивалентной (в области локализации вихревых структур) пространственно-периодической задаче для уравнения (2.4) в полосе:

где $\varepsilon : = {{({\text{P}}{{{\text{r}}}_{D}}{\text{Re}})}^{{ - 1}}}\sim {{10}^{{ - 3}}} > 0$ – малый параметр, ${{L}_{1}}$, ${{L}_{2}}$, $a$, $b$ – некоторые положительные числа (параметры расчетной области), функции ${{u}_{i}}(x) = {{u}_{i}}({{x}_{2}},{{x}_{3}})$, $i = \overline {1,3} $ – ${{L}_{2}}$-периодические по переменной ${{x}_{2}}$, непрерывно дифференцируемые в области ${{\mathbb{R}}^{1}} \times [a,b]$, ${{C}^{a}}({{x}_{1}},{{x}_{2}})$ – непрерывно дифференцируемая ${{L}_{1}}$-периодическая по переменной ${{x}_{1}}$ и ${{L}_{2}}$-периодическая по переменной ${{x}_{2}}$ функция, n – внутренняя нормаль к границе ${{x}_{3}} = b$. Ось ${{x}_{1}}$ ориентирована по оси симметрии вихревых структур.

3. Численное моделирование распределения концентрации аэрозоля. Существуют разные возможности (в смысле выбора средств) в моделировании распределения аэрозоля внутри вихревых структур. Рассмотрим два возможных подхода к решению этой задачи.

Алгоритм построения численного решения задачи (2.6) основан на использовании свойства единственности и глобальной устойчивости решения этой задачи как стационарного решения соответствующей начально-краевой задачи для уравнения (2.5) и реализуется путем стационирования решения задачи для уравнения (2.5) к решению задачи (2.6). В свою очередь, численное решение задачи для уравнения (2.5) выполняется с использованием метода эволюционной факторизации [37] в расчетной области, представляющей собой прямоугольный параллелепипед $\Pi = \{ ({{x}_{1}},{{x}_{2}},{{x}_{3}})$: $0 \leqslant {{x}_{1}},{{x}_{2}} \leqslant 49$, $0.3 \leqslant {{x}_{3}} \leqslant 12\} $, в котором введена равномерная сетка, состоящая из ${{N}_{1}} \times {{N}_{2}} \times {{N}_{3}}$ узлов. Безразмерный шаг сетки соответствует пространственному расстоянию в 0.1 км. Пространственные шаги сетки одинаковы по каждому из трех направлений: ${{h}_{x}} = {{h}_{y}} = {{h}_{z}}$ = = 0.1 (в безразмерных переменных). Безразмерный шаг по времени выбирается равным $\tau = 0.01$. Точность метода эволюционной факторизации при таком выборе пространственного и временного шага составляет $O\left( {{{\tau }^{2}} + h_{x}^{2} + h_{y}^{2} + h_{z}^{2}} \right)$ [37]. Поэтому численные расчеты проводились с точностью 0.01. По вертикали расчетная область начинается с уровня ${{x}_{3}} = 0.3$, поскольку вынос аэрозоля возникает на субмезомасштабных движениях, которые “не разрешимы” в используемом приближении. К структурам, обеспечивающим такое движение относятся, например, спонтанные вихри с вертикальной осью – пыльные дьяволы и термики. Они тянут частицы аэрозоля вверх.

В ходе реализации метода эволюционной факторизации трехмерный оператор Лапласа заменяется на произведение трех одномерных, что позволяет решить уравнение в три этапа, на каждом из которых решение уравнения сводится к решению ${{N}_{i}} \times {{N}_{j}}$, $i,j = 1,2,3$ обыкновенных дифференциальных уравнений методом прогонки. Метод прогонки реализуется в соответствии с алгоритмом из работы [38].

На рис. 3–4 а)–в) представлены компоненты u₁, u₂, u₃ вектора скорости переноса [33], рассчитанные при условиях нейтральной стратификации для достаточно характерных для АПС значениях ${\text{Re}} = 200$, угле поворота по отношению к геострофическому ветру 15°. На рис. 3–4 г)–е) представлены результаты по расчету безразмерной концентрации $C(x)$, полученные при численном решении задачи (2.6) при различных значениях безразмерных параметров ${{{v}}_{g}}$, $\gamma $ и $\varepsilon = 0.004$ (в сечении вертикальной плоскостью). Рисунки 3–4 ж) или з) – разности безразмерных концентраций, которые представлены распределениям г) (уменьшаемое) и д) (вычитаемое) или д) (уменьшаемое) и е) (вычитаемое). Расчеты выполнены при различных граничных значениях концентрации ${{C}^{a}}$, равных 10 и 100.

Отметим, что согласно использованным данным [33] на верхней границе ${{x}_{3}} = 12$ (в безразмерных единицах) величина компоненты скорости u₃ равна нулю, а рассчитанная концентрация постоянна и совпадает с фоновой. В связи с этим на рисунках представлены только те области изменения соответствующих физических величин, на которых наблюдаются распределения, отличные от постоянных.

В работах [13, 34] показано, что значительное количество аэрозоля попадает в приосевую область вихрей. Согласно численным расчетам, представленным на рис. 3, 4, распределение концентрации аэрозоля имеет пространственно-периодическую структуру (что соответствует данным наблюдений [10]), согласованную с периодической структурой поля скоростей [19, 33]. Основная масса увлекаемого вихрями аэрозоля удерживается ими и переносится в горизонтальной плоскости. Смещение области локализации вещества к нижней границе (по сравнению с результатами моделирования [13, 34]) связано с учетом распада и оседания вещества в модели (2.6). Учет процесса оседания аэрозоля приводит к незначительным изменениям в распределении концентрации этой примеси, не превышающим $1\% $ от граничного значения ${{C}^{a}}$ (по данным моделирования), причем наибольшее отличие величин концентраций при ${{{v}}_{g}} = 0$ и ${{{v}}_{g}} = 0.001$ приходится на области с наибольшими значениями компонент скорости переноса ${{u}_{2}}$ и ${{u}_{3}}$, которые определяют структуру распределения примеси по вертикали (см. рис. 3 б), в), ж) и рис. 4 б), в), ж)). На эти же области приходится наибольшее отличие величин концентраций при различных скоростях распада: $\gamma = 0.1$ и $\gamma = 1$ (см. рис. 3 д), е), з) и рис. 4 д), е), з)). Покажем, что при таком выборе параметров модели (2.6) процессы распада и переноса аэрозоля превалируют над диффузией и определяют структуру поля концентраций аэрозоля в вихрях.

Рассмотрим следующую задачу Коши

которая получается из задачи (2.6) при $\varepsilon = 0$ без учета условия на верхней границе при ${{x}_{3}} = b$. Задача (3.1) разрешима единственным образом в классе функций ${{C}^{1}}\left( {{{\mathbb{R}}^{2}} \times (a,b)} \right)$ ∩ $C\left( {{{\mathbb{R}}^{2}} \times [a,b]} \right)$, так как одна из горизонтальных компонент скорости переноса всегда отлична от нуля [39, с. 46]. Заметим, что две другие компоненты скорости переноса могут обращаться в нуль. Получим оценку точности приближенного решения $\hat {C}(x,\varepsilon ) = {{\bar {C}}_{0}}(x)$ задачи (2.6) при $\gamma = 1$, используя метод верхних и нижних решений (см., например, [20]).

Определение. Упорядоченная пара ${{L}_{1}}$-периодических по переменной ${{x}_{1}}$ и ${{L}_{2}}$-периодических по переменной ${{x}_{2}}$ функций $\beta (x,\varepsilon )$, $\alpha (x,\varepsilon )$ ∈ ${{C}^{1}}\left( {{{\mathbb{R}}^{2}} \times [a,b]} \right)$ ∩ ${{C}^{2}}\left( {{{\mathbb{R}}^{2}} \times (a,b)} \right)$ называется соответственно верхним и нижним решениями задачи (2.6), если при $0 < \varepsilon \ll 1$ выполняются следующие дифференциальные неравенства:

Известно, что если существуют верхнее и нижнее решения задачи (2.6), то решение $C(x,\varepsilon )$ этой задачи заключено между верхним и нижним решениями (см. напр., [20]):

Верхнее и нижнее решения задачи (2.6) получим путем модификации приближенного решения $\hat {C}(x,\varepsilon )$ определенным образом:

где ${{\bar {C}}_{0}}(x)$ – решение задачи (3.1), $\varepsilon = 0.004$, $\delta = 2 \times {{10}^{{ - 2}}}$, $B = {{10}^{{ - 2}}}$, $\xi = \left( {{{x}_{3}} - b} \right){\text{/}}\varepsilon $.

Лемма. Функции $\beta (x,\varepsilon )$, $\alpha (x,\varepsilon )$, определяемые равенствами (3.3), являются верхним и нижним решениями задачи (2.6).

Для доказательства леммы достаточно проверить выполнение неравенств 1)–3). Подставляя функции (3.3) в неравенства 1)–2), имеем:

При ${{x}_{3}} = a$ и ${{x}_{3}} = b$ имеем неравенства:

выполнение которых обеспечивается за счет соответствующего выбора параметра модели b. Лемма доказана. Итак, имеет место двойное неравенство (3.2).

Получим оценку остаточного члена. Поскольку $\beta (x,\varepsilon ) - \alpha (x,\varepsilon ) = O(\mu )$ и β(x, ε) – ‒ ${{\bar {C}}_{0}}(x)$ = $O(\mu )$, где $\mu \sim {{10}^{{ - 2}}}$, то в силу неравенства треугольника:

в равномерной норме.

Теорема. Существует единственное классическое решение $C(x,\varepsilon )$ задачи (2.6), удовлетворяющее при $\gamma = 1$ оценке (3.4), где ${{\bar {C}}_{0}}(x)$ – решение задачи (3.1).

Заметим, что в случае $\gamma = 0.1$ справедлива аналогичная теорема с оценкой (3.4), в которой $\mu \sim {{10}^{{ - 1}}}$. Доказательство теоремы легко получить, если положить $\delta $ = $2 \times {{10}^{{ - 1}}}$.

Таким образом, в качестве модельной задачи может быть использована задача (3.1), причем допущенная при этом ошибка в решении определяется неравенством (3.4). Для численного решения задачи (3.1) можно использовать метод конечных элементов [40], реализуемый в FEniCS Project (вычислительная платформа с открытым исходным кодом (LGPLv3) для решения уравнений в частных производных) через API для Python, дающий результат вычислений близкий (в соответствии с теоремой) к результату, представленному на рис. 3, 4.

4. Обсуждение результатов. Оценка количества переносимого аэрозоля. Из общих теорем теории линейных уравнений 2-го порядка в частных производных [21–23] следует, что состояние системы, которому соответствует стационарное распределение концентрации аэрозоля, изменившись по причине внешнего возмущения, возвращается в исходное стационарное состояние в силу устойчивости и единственности стационарного решения, описывающего это состояние. Таким образом, стационарное распределение концентрации аэрозоля отвечает наиболее вероятному состоянию этой системы.

В результате численного моделирования установлено, что поле концентрации аэрозоля в вихревых структурах имеет пространственно-периодическую структуру, согласованную с периодическим изменением поля скоростей [19, 33], что также согласуется с данными наблюдений [10]. Учет распада и оседания аэрозоля приводит к смещению области локализации вещества по вертикали в направлении нижней границы расчетной области (по сравнению с данными моделирования [13, 34]).

Количество вещества Q, переносимое через сечение $S \approx (50 \times 4)$ км² за 1 мин определяется формулой (в размерных переменных):

где ${{m}_{0}}$ – масса молекулы аэрозоля, $S{\text{'}}$ – прямоугольная область: $({{x}_{2}},{{x}_{3}})$ ∈ ∈ $[0,\;45] \times [0.3,\;4]$.

В разд. 3 на рис. 3, 4 приведены вычисленные значения интеграла

для поля скоростей [19, 33] и различных значений параметров модели (2.6).

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 18-29-10080 и Фонда развития теоретической физики и математики “Базис” в рамках гранта № 19-2-6-178-1.

Авторы выражают глубокую благодарность проф. В.Ф. Бутузову за внимание к работе и полезные рекомендации.