Развитие полного лагранжева подхода для моделирования течений разреженных дисперсных сред (обзор)

Полный текст

ВВЕДЕНИЕ

В различных разделах физики и механики используются континуальные бесстолкновительные модели сред, лишенных собственного давления [1]. Пренебрежение “напряжениями” при введении континуального описания среды, состоящей из дискретных элементов (частиц), справедливо в случае, когда флуктуационными скоростями соседних частиц можно пренебречь по сравнению с их среднемассовой скоростью. Такие модели успешно применялись для описания динамики самогравитирующих дискретных масс (крупномасштабное распределение вещества во Вселенной, “блины Зельдовича” [2], формирование планетных систем, волны плотности и спиральная структура галактик [3]), “холодной” и пылевой компонент в плазме [4], коллективного движения птиц, насекомых и микроорганизмов [5], моделирования транспортных потоков [6] и ряда других систем.

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

Можно упомянуть также работы [8, 9], в которых рассматривались среды без давления (либо с давлением, зависящим только от времени [10]) с условием “слипания” частиц среды при пересечении их траекторий. Иногда в литературе такое условие, применимое при значительных объемных концентрациях частиц, называется условием запрещенного обгона [1]. При таком ограничении возникает необходимость вводить в рассмотрение сильные разрывы [9], несущие конечную массу, импульс и энергию поверхностной фазы, движущейся вдоль разрыва [8].

Течения разреженных двухфазных сред типа “газ/жидкость — инерционные частицы/капли/пузырьки”, как правило, описываются в рамках двухконтинуального подхода, при этом несущая фаза, в общем случае, описывается уравнениями Навье–Стокса с источниковыми членами, учитывающими обратное влияние частиц, а дисперсная фаза моделируется континуумом, лишенным собственных напряжений [11–16]. Такой континуум является чрезвычайно “cжимаемым”. Более того, как отмечалось ранее, в разреженной “холодной” среде возможно возникновение зон пересекающихся траекторий частиц, на границах которых возникают каустики — огибающие траекторий. Вблизи каустик, а также точек или линий “сборки” континуума частиц происходит резкое возрастание числовой концентрации (осредненной плотности) среды частиц в силу локального схлопывания трубок тока дисперсной фазы [1, 7, 17].

В литературе для расчета параметров дисперсной фазы часто используются конечно-разностные методы, основанные на эйлеровом описании среды частиц. Такой подход предполагает однозначность полей континуальных параметров дисперсной фазы, поэтому он неприменим для исследования течений с пересекающимися траекториями частиц. Другой известный подход — использование лагранжева метода типа “частицы-в-ячейках” [18], основанного на отслеживании многих траекторий отдельных частиц без учета уравнения неразрывности дисперсной фазы. Концентрация дисперсной фазы при этом вычисляется путем суммирования всех частиц, приходящихся на одну эйлерову ячейку несущей фазы. Иногда малый лагранжев объем дисперсной среды заменяют “крупной частицей”, представляющей группу соседних частиц. Масса “крупной частицы” считается неизменной в процессе движения.

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

Альтернативой является полный лагранжев подход (ПЛП), предложенный автором ранее — см., например, работы [17, 19, 20]. Этот подход использует лагранжево описание полей скорости и числовой плотности дисперсной фазы. При этом уравнение неразрывности среды частиц, записанное в лагранжевой форме, решается путем привлечения дополнительных уравнений для нахождения компонент якобиана перехода от эйлеровых к лагранжевым переменным. Эти уравнения решаются совместно с уравнениями динамики частиц на выбранных траекториях дисперсной фазы.

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

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

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

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

Разделы 9–11 посвящены развитию бессеточных, полностью лагранжевых подходов для расчета параметров несущей и дисперсной фазы при исследовании нестационарных вихревых газодисперсных течений с несжимаемой вязкой несущей фазой. Описан бессеточный лагранжев метод, основанный на комбинации ПЛП для дисперсной фазы и метода вязких вихревых доменов для несущей фазы. Обсуждаются примеры расчета плоских вихревых двухфазных течений, течений с вихревыми кольцами и импульсных газодисперсных струй, в том числе — с учетом полидисперсности и испарения дисперсной фазы.

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

1. О МОДЕЛИ РАЗРЕЖЕННОЙ ДВУХФАЗНОЙ СРЕДЫ

Рассматривается смесь, состоящая из несущей фазы (жидкости или газа) и сферических включений, имеющих радиус σ и массу m. В общем случае при учете фазовых переходов радиус и масса частиц могут изменяться в процессе движения среды. Температура вещества частиц T_s предполагается однородной по объему частицы. Объемная концентрация дисперсной фазы считается пренебрежимо малой, столкновениями частиц и броуновскими эффектами пренебрегается. Режим обтекания отдельных частиц может варьировать в широких пределах: от континуального до свободномолекулярного. В общем случае при наличии распределения частиц по размерам, скоростям и температурам осредненное поведение дисперсной фазы может быть описано с помощью бесстолкновительного кинетического уравнения для одночастичной функции распределения f ⁽¹⁾(t, r, c, T_s, σ):

$\frac{\partial f^{\left(1\right)}}{\partial t}+\mathbf{c}\frac{\partial f^{\left(1\right)}}{\partial \mathbf{r}}+\frac{\partial f^{\left(1\right)}(\mathbf{f}/m)}{\partial \mathbf{c}}+\frac{\partial f^{\left(1\right)}(q/c_{s}m)}{\partial T_s}+\frac{\partial f^{\left(1\right)}j}{\partial \sigma }=0.$ (1.1)

Здесь t — время, c — скорость, r — радиус-вектор, f — сила, действующая на пробную частицу, q — приток тепла к частице со стороны несущей фазы, c_s — теплоемкость вещества частицы, j — скорость изменения радиуса частицы (за счет испарения, горения или конденсации). Частные производные по векторным переменным для краткости обозначают дивергенцию по компонентам соответствующих векторов. В общем случае в число фазовых переменных могут быть добавлены и другие параметры, например, угловая скорость вращения частиц, параметр формы и пр.

Числовая концентрация n_s и макроскопическая скорость V_s дисперсной фазы определяются стандартным образом:

$n_s(t,\mathbf{r})=\int \begin{array}{cc}{f}^{\left(1\right)}d\mathbf{c}d\sigma d{T}_s,& n_s{\mathbf{V}}_s=\int f^{\left(1\right)}\mathbf{c}d\mathbf{c}d\sigma d{T}_s\end{array}.$

Индекс s здесь и далее относится к макропараметрам среды частиц.

На макромасштабе несущая фаза в общем случае описывается уравнениями Навье–Стокса с дополнительными источниковыми членами Q₁, Q₂ и Q₃, которые учитывают межфазный обмен массой, импульсом и энергией. Если массообмен отсутствует ( j = Q₁ = 0), источниковые члены в правых частях уравнения импульса и уравнения энергии, записанного в форме уравнения притока тепла, имеют вид

$Q_2=-\int f^{\left(1\right)}\mathbf{f}d\mathbf{c}d\sigma d{T}_s,Q_3=-\int f^{\left(1\right)}qd\mathbf{c}d{T}_{s}d\sigma +\int f^{\left(1\right)}\mathbf{f}(\mathbf{c}-\mathbf{V})d\mathbf{c}d{T}_{s}d\sigma .$

Здесь V — скорость несущей фазы. В последнем соотношении первый член описывает приток тепла от частиц к несущей фазе, а второй — работу сил трения на относительном перемещении фаз. В частном случае отсутствия фазовых переходов и распределения частиц по температурам кинетико-континуальное описание двухфазной среды использовалось, например, в работах [21–23]. В абсолютном большинстве публикаций по моделированию дисперсных сред используются дальнейшие упрощения: считается, что фазовые переходы отсутствуют, а среду частиц можно разбить на конечное число сортов частиц с концентрациями n_sk. Частицы каждого сорта имеют одинаковый размер, одинаковые локальные скорости V_sk и температуру T_sk (k — номер сорта частиц). Здесь, как и ранее, T_sk — температура материала частиц, а не термодинамическая температура континуума частиц, которая для каждого сорта частиц в принятых предположениях равна нулю. Предполагается, что функция распределения по скоростям и температурам для частиц k-го сорта имеет вид дельта-функции Дирака

$f^{\left(1\right)}=n_{sk}(t,\mathbf{r})\delta (\mathbf{c}-{\mathbf{V}}_{sk})\delta (T_s-T_{sk}).$

Для такой функции распределения из кинетического уравнения для k-го сорта частиц стандартным образом, после интегрирования по скоростям и температурам, получаются уравнения “холодного” континуума (индекс k в дальнейшем опущен)

$\frac{\partial n_s}{\partial t}+\mathrm{div}\left(n_s{\mathbf{V}}_s\right)=\begin{array}{cc}\mathrm{0,}& m\end{array}\frac{d{\mathbf{V}}_s}{dt}=\begin{array}{cc}{\mathbf{f}}_s,& c_s\end{array}m\frac{d{T}_s}{dt}=q_s.$ (1.2)

При формулировке уравнений “холодной” среды вместо числовой плотности частиц n_s часто используют массовую плотность ρ_s = mn_s и переписывают уравнения (1.2) в более привычном для механики сплошной среды виде:

$\frac{\partial {\rho }_s}{\partial t}+\mathrm{div}\left({\rho }_s{\mathbf{V}}_s\right)=\mathrm{0,}{\rho }_s\frac{d{\mathbf{V}}_s}{dt}={\mathbf{F}}_s,c_s{\rho }_s\frac{d{T}_s}{dt}=Q_s,{\mathbf{F}}_s=n_s{\mathbf{f}}_s,Q_s=n_s{q}_s.$

В данном случае скорость континуума частиц V_s совпадает с реальной скоростью пробной частицы, f_s — осредненная сила, действующая на пробную частицу (в зависимости от условий обтекания частиц f_s может задаваться различными выражениями и включать составляющие негидродинамической природы), q_s — осредненный тепловой поток к пробной частице. В качестве f_s и q_s обычно используют выражения для силы и притока тепла для одиночной сферы в вязком потоке несущей фазы, параметры которого (V и T ) находятся из осредненных уравнений несущей фазы. Выражения для f_s и q_s зависят от локальных условий обтекания частицы, т. е., в основном, от локальных значений чисел Рейнольдса и Маха обтекания частиц:

${\mathrm{Re}}_s=\frac{2\sigma \rho \left|\mathbf{V}-{\mathbf{V}}_s\right|}{\mu }\begin{array}{cc},& {\text{M}}_s\end{array}=\frac{\left|\mathbf{V}-{\mathbf{V}}_s\right|}{a}.$

Здесь ρ, μ и a — плотность, динамическая вязкость и скорость звука несущей фазы. Например, при обтекании малой частицы (капли, пузырька) вязкой средой при малых числах Рейнольдса суммарная сила, действующая на частицу со стороны несущей фазы, имеет вид [24]

$$\begin{gathered} {\mathbf{f}}_s={\mathbf{f}}_{St}+{\mathbf{f}}_{Ar}+{\mathbf{f}}_{vm}+{\mathbf{f}}_{BB};{\mathbf{f}}_{St}=6\pi \sigma \mu (\mathbf{V}-{\mathbf{V}}_s)\begin{array}{cc},& {\mathbf{f}}_{Ar}\end{array}=\frac{4}{3}\pi {\sigma }^3\rho \left(\frac{d\mathbf{V}}{dt}-\mathbf{g}\right), \\ {\mathbf{f}}_{vm}=\frac{2}{3}\pi {\sigma }^3\rho \left(\frac{d\mathbf{V}}{dt}-\frac{d_s{\mathbf{V}}_s}{dt}\right);{\mathbf{f}}_{BB}=6{\sigma }^2\sqrt{\pi \rho \mu }\underset{0}{\overset{t}{\int }}\frac{d_s(\mathbf{V}-{\mathbf{V}}_s)/d\tau }{\sqrt{t-\tau }}/d\tau \text{.} \end{gathered}$$ (1.3)

Здесь f_St, f_Ar, f_vm и f_BB — силы Стокса, Архимеда, присоединенных масс и Бассэ–Буссинеска соответственно. Индекс s при производной здесь обозначает дифференцирование вдоль траектории частицы. Для газовзвесей и аэрозолей первая сила (сила аэродинамического сопротивления при стационарных условиях обтекания) значительно превосходит все остальные. При конечных числах Рейнольдса обтекания частицы ее принято представлять в виде

${\mathbf{f}}_d=\frac{1}{2}{C}_d\rho \left|\mathbf{V}-{\mathbf{V}}_s\right|(\mathbf{V}-{\mathbf{V}}_s)\pi {\sigma }^2.$

Здесь C_d — коэффициент сопротивления, в общем случае зависящий от чисел Re_s и M_s обтекания частицы. При малых числах Маха для C_d часто используется формула Л.С. Клячко [25], удовлетворительно аппроксимирующая стандартную кривую сопротивления сферы до Re_s ~ 10³:

$C_d=\frac{24}{{\mathrm{Re}}_s}\left(1+\frac{1}{6}{\mathrm{Re}}_s^{2/3}\right).$

В условиях обтекания частицы в режиме сплошной среды при конечных M_s для C_d широко используются аппроксимационная формула из работы [26]. Выражение для C_d при свободномолекулярном режиме обтекания частиц можно найти, например, в работе [27]. При учете вращения частиц и движении частиц в узких областях с большими градиентами продольной скорости (пограничных слоях, слоях смешения, узких каналах) важную роль играют подъемные силы, действующие на частицы (см. примеры в разд. 8). Среди сил негидродродинамической природы, влияющих на динамику частиц, часто учитываются сила тяжести mg (здесь g — в общем случае, градиент гравитационного потенциала) и сила Лоренца (в системе СИ), действующая на заряженную частицу с зарядом e во внешнем электромагнитном поле с напряженностями E и B [23, 28]

${\mathbf{f}}_L=e\left[\mathbf{E}+(\mathbf{V}-{\mathbf{V}}_s)\mathbf{B}\right].$

Для теплового потока к частице чаще всего используется аппроксимационная формула Ранца–Маршалла [29]:

$q_s=4\pi \sigma \lambda (T-T_s)G({\mathrm{Re}}_s,\mathrm{Pr},{\text{M}}_s),G=1+0.3{\mathrm{Re}}_s^{1/2}{\mathrm{Pr}}^{1/3}.$

В случае конечных чисел Маха обтекания частиц для поправочной функции G можно использовать аппроксимационную формулу из [26].

В уравнениях несущей фазы источниковые члены, соответствующие одному сорту частиц, принимают вид

$Q_2=-n_s{\mathbf{f}}_s,Q_3=-n_s{q}_s+n_s{\mathbf{f}}_s({\mathbf{V}}_s-\mathbf{V}).$ (1.4)

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

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

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

При использовании уравнений типа (1.2) основные математические трудности связаны с очень высокой “сжимаемостью” континуума невзаимодействующих частиц. Более того, в некоторых областях течения траектории частиц могут пересекаться в одних и тех же точках пространства, формируя “сборки” и “складки”, на границах которых (каустиках) концентрация дисперсной фазы неограниченно возрастает [1, 7, 17].

Вопрос о границах применимости модели невзаимодействующих частиц при наличии зон накопления или пересечения траекторий частиц обычно решается на основе дополнительного микроструктурного анализа. Примеры такого анализа приведены в работе [17], где с использованием элементов теории вероятности, в предположении пуассоновского распределения числа частиц в конечном объеме, получены выражения для среднего значения (математического ожидания) расстояния между частицами M(l) и его дисперсии D(l) в точках неограниченного роста числовой плотности частиц на каустиках:

$\frac{\sigma }{M\left(l\right)}=\frac{\text{s}}{\Gamma (1/s)}{\left[\frac{3k}{4\pi }{\phi }_{\text{s}\infty }{\left(\frac{L}{\sigma }\right)}^{3-s}\right]}^{1/s}\text{},\frac{D\left(l\right)}{M^2\left(l\right)}=\frac{2s\Gamma (2/s)}{{\Gamma }^2(1/s)}-1.$ (1.5)

Здесь Γ — гамма-функция Эйлера, L — макромасштаб длины для рассматриваемой задачи, на котором среда частиц описывается континуальными уравнениями, ${\phi }_{\text{s}\infty }=(4/3)\pi {\sigma }^3{n}_{s\infty }$ — характерная объемная доля среды частиц вдали от сингулярности концентрации, а положительные константы (интегральные характеристики сингулярности) k и s определяются из следующих соображений.

Рассмотрим точку x₀, лежащую на каустике, в которой числовая плотность среды частиц неограниченно (но интегрируемым образом) возрастает, и окружим эту точку сферой W_r радиуса r. Предположим, что при малых r/L среднее число частиц в сфере W_r имеет следующее асимптотическое представление:

$N_s(r/L)=\underset{{\Omega }_r}{\int }{n}_{s}dv=k{\left(\frac{r}{L}\right)}^s\text{}+o{\left(\frac{r}{L}\right)}^s\text{}.$

Из выражения (1.5) следует, что, несмотря на неограниченный рост числовой плотности на каустике, среднее расстояние между частицами остается конечным, поскольку конечный лагранжев объем среды частиц “схлопывается” не в точку, а в конечный элемент поверхности с конечной поверхностной плотностью среды частиц. Оценки показывают [17], что для большинства типичных сингулярностей плотности среды частиц, возникающих при моделировании течений газовзвесей в условиях, представляющих интерес для приложений, среднее расстояние между частицами в точках сингулярности числовой плотности остается много больше размера частиц, т. е. бесстолкновительная модель дисперсной фазы остается применимой.

На каустиках, являющихся огибающими траекторий частиц в точках разворота траекторий и лежащих в области конечных значений скорости несущей фазы, числовая плотность частиц имеет интегрируемую сингулярность типа $1/\sqrt{|\mathbf{x}-{\mathbf{x}}_0|}$ [30], где $|\mathbf{x}-{\mathbf{x}}_0|$ — расстояние до точки сингулярности x₀. При учете небольшого распределения частиц по скоростям (малого отличия исходной функции распределения частиц по скоростям от дельта-функции) числовая плотность среды частиц вблизи каустики становится конечной, оставаясь достаточно большой [31].

Имеющиеся в литературе оценки границ применимости бесстолкновительной среды для моделирования течений с локальными зонами пересекающихся траекторий частиц показывают, что в условиях, представляющих интерес для большинства приложений, для инерционных частиц (размером более 1 мкм) бесстолкновительная модель остается применимой при исходной объемной концентрации, не превосходящей величину ~ 10^–4–10^–5 (см., например, работы [31–32]).

С увеличением концентрации дисперсной фазы возрастает вероятность столкновений между частицами, движущимися по пересекающимся траекториям, например, вблизи тела, обтекаемого запыленным газом с инерционными частицами, отражающимися от обтекаемой поверхности. Существует довольно много публикаций, посвященных построению полуэмпирических столкновительных моделей дисперсных систем (см., например, [33–35]), в которых используются статистические и кинетические подходы и предположения, во многом аналогичные теории газовых смесей. В частности, как правило, предполагается быстрая “максвеллизация” функции распределения сталкивающихся частиц по скоростям.

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

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

Ниже рассматривается предельная модель полностью бесстолкновительной среды частиц, справедливая для достаточно разреженных дисперсных систем.

2. ОСНОВНЫЕ ИДЕИ ПОЛНОГО ЛАГРАНЖЕВА ПОДХОДА

Для вычисления концентрации частиц в зонах пересекающихся траекторий (“складок”) частиц и вблизи каустик естественно перейти от эйлерова к лагранжеву описанию среды частиц. В качестве лагранжевых координат для простоты используем значения декартовых координат выбранной лагранжевой частицы x₀, y₀, z₀ в некоторый момент, принятый за начало отсчета t = 0. Рассматривая радиус-вектор выбранной частицы r_s, скорость V_s и температуру T_s как функции лагранжевых координат r₀ = (x₀, y₀, z₀) и времени t, уравнения (1.2) можно записать в лагранжевой форме:

$\frac{\partial {\mathbf{r}}_s}{\partial t}=\begin{array}{cc}{\mathbf{V}}_s,& m\end{array}\frac{\partial {\mathbf{V}}_s}{\partial t}=\begin{array}{cc}{\mathbf{f}}_s,& c_s\end{array}m\frac{\partial T_s}{\partial t}=q_{s,}{n}_s({\mathbf{r}}_0,t)\left|\det \left(J\right)\right|=n_{s0}\left({\mathbf{r}}_0\mathrm{,0}\right).$ (2.1)

Здесь компоненты матрицы Якоби J имеют вид J_ij = ∂x_i / ∂x_0j (i, j = 1, 2, 3; x₁ = x, x₂ = y, x₃ = z). При фиксированном значении r₀ первые три уравнения превращаются в систему обыкновенных дифференциальных уравнений, позволяющих рассчитать данную траекторию частиц и распределение скорости и температуры дисперсной фазы вдоль этой траектории. Если бы компоненты якобиана J_ij на траектории были известны, то из последнего соотношения (2.1) можно было бы найти распределение концентрации частиц вдоль данной траектории. В большинстве моделей разреженных дисперсных сред межфазная сила f_s считается известной функцией координат, скоростей и температур фаз (а также производных параметров несущей фазы по эйлеровым координатам и времени), но не зависит явно от концентрации частиц. Последнее условие позволяет вывести замкнутую систему обыкновенных дифференциальных уравнений для компонент J_ij на траектории частиц. Продемонстрируем вывод на примере стоксовой межфазной силы. Пусть

${\mathbf{f}}_s=\frac{m}{\tau }(\mathbf{V}-{\mathbf{V}}_s\begin{array}{cc}),& \tau \end{array}=\frac{m}{6\pi \sigma \mu }=\frac{2{\sigma }^2{\rho }_{\text{s}}^{\text{0}}}{\text{9}\mu }.$

Здесь τ — время скоростной релаксации частиц, ρ_s⁰ — плотность вещества частиц. Дифференцируя первые два векторных уравнения (2.1) по лагранжевым координатам x_0j и меняя порядок дифференцирования по времени и x_0j, получаем систему

$\frac{\partial J_{ij}}{\partial t}=\begin{array}{cc}{\Omega }_{ij},& \frac{\partial {\Omega }_{ij}}{\partial t}\end{array}=\frac{1}{\tau }\left(\sum _k\frac{\partial v_i}{\partial x_k}{J}_{kj}-{\Omega }_{ij}\right),{\Omega }_{ij}=\frac{\partial v_{si}}{\partial x_{0j}}.$ (2.2)

Таким образом, если скорость несущей фазы и ее производные по эйлеровым координатам известны, уравнения (2.1) совместно с (2.2) составляют замкнутую систему обыкновенных дифференциальных уравнений (ОДУ) и конечного соотношения, позволяющих определить все параметры дисперсной фазы, включая концентрацию, на выбранной траектории частиц.

Для трехмерных неустановившихся течений порядок системы ОДУ для расчета одной траектории, компонент скорости, температуры и вспомогательных функций J_ij, Ω_ij, необходимых для расчета концентрации частиц вдоль этой траектории, равен 25, для двумерных — 13. Начальные условия для указанной системы ОДУ получают естественным образом из граничных условий для рассматриваемого течения. Например, если расчет траектории начинается из области однородного потока частиц, направленного вдоль оси x, то при t = 0 имеем r = r₀, T_s = T_s0, n_s = n_s0, u_s = u_s0, v_s = 0, w_s = 0, J_ii = 1, J_ij (i ≠ j) = 0, Ω_ij = 0. При использовании других моделей межфазного взаимодействия уравнения для вспомогательных функций Ω_ij имеют более сложный вид (см., например, [23, 27, 38–41], но описанная выше процедура вывода этих уравнений остается неизменной.

Предложенный метод расчета концентрации частиц устраняет проблему пересекающихся траекторий, так как различным траекториям, пересекающимся в одной и той же точке физического пространства, соответствуют различные значения лагранжевых координат r₀ и t. Границы областей пересекающихся траекторий (каустики) легко определяются в процессе расчета, поскольку на них якобиан обращается в ноль.

Расчет концентрации можно продолжить за точку пересечения траекторий частиц, при этом для автоматического учета возникновения “складки” определитель в уравнении неразрывности среды частиц, в общем случае, следует взять по модулю, как это сделано в выражении (2.1), поскольку за точкой пересечения траекторий меняется направление обхода лагранжева элемента. Если траектория частиц пересекает поверхность сильного разрыва параметров несущей фазы (например, ударную волну [41]), правые части ОДУ (2.1)–(2.2) и переменные J_ij, Ω_ij, также претерпевают разрывы. При этом для продолжения расчетов параметров дисперсной фазы после разрыва на выбранной траектории следует учитывать непрерывность координат и скоростей среды частиц, а также их производных вдоль поверхности разрыва.

3. РЕАЛИЗАЦИЯ ПОЛНОГО ЛАГРАНЖЕВА ПОДХОДА В КРИВОЛИНЕЙНЫХ КООРДИНАТАХ

Если для описания движения дисперсной смеси используются криволинейные координаты ξ₁, ξ₂, ξ₃, то уравнение неразрывности среды частиц в лагранжевой форме имеет вид [42]

$n_s({\xi }_{10},{\xi }_{20},{\xi }_{30},t)\sqrt{g}=n_s({\xi }_{10},{\xi }_{20},{\xi }_{30}\mathrm{,0})\sqrt{g_0},g=\det ||g_{ij}||.$ (3.1)

Здесь g_ij — компоненты метрического тензора лагранжевой системы координат, индекс ноль отмечает значение координат и определителя при t = 0. Используя векторные уравнения из (2.1), можно вывести замкнутую систему ОДУ на фиксированной траектории для определения текущего значения n_s. В общем случае эта процедура слишком громоздка, поэтому приведем ее для криволинейных ортогональных координат, наиболее часто встречающихся в приложениях. В этом случае уравнение неразрывности среды частиц (2.1) принимает вид

$n_s({\xi }_{10},{\xi }_{20},{\xi }_{30},t)H_1{H}_2{H}_3\det ||\partial {\xi }_i/\partial {\xi }_{i0}||=n_{s0}({\xi }_{10},{\xi }_{20},{\xi }_{30}\mathrm{,0})H_{10}{H}_{20}{H}_{30}.$ (3.2)

Здесь H_i — коэффициенты Ламе. Введем в рассмотрение физические компоненты скорости v_si и ускорения a_si частицы. Воспользуемся известным представлением компонент ускорения материальной точки единичной массы через ее кинетическую энергию K:

$a_{si}=\frac{1}{H_i}\left[\frac{d}{dt}\left(\frac{\partial K}{\partial {\dot{\xi }}_i}\right)-\frac{\partial K}{\partial {\xi }_i}\right]\begin{array}{cc},& K\end{array}=\frac{|V{|}_s^2}{2}.$ (3.3)

Здесь точка над ξ_i обозначает производную по времени. С использованием выражения (3.3) можно записать уравнения траектории и импульса частиц в проекциях на оси криволинейной ортогональной системы координат следующим образом:

$\frac{\partial {\xi }_i}{\partial t}=\frac{v_{si}}{H_i}\begin{array}{cc},& \frac{\partial v_{si}}{\partial t}\end{array}=\sum _k\frac{v_{sk}}{H_i{H}_k}\left(v_{sk}\frac{\partial H_k}{\partial {\xi }_i}-v_{si}\frac{\partial H_i}{\partial {\xi }_i}\right)+\frac{f_{si}}{m}.$ (3.4)

В уравнениях (3.4) суммирование по i отсутствует, искомые функции ξ_i и v_si считаются зависящими от времени t и лагранжевых координат ξ_j 0, f_si — соответствующая компонента силы, действующей на частицу. Дифференцируя уравнения (3.4) по лагранжевым координатам ξ _j 0 и изменяя порядок дифференцирования, получаем систему

$$\begin{gathered} \frac{\partial J_{ij}}{\partial t}=\frac{{\Omega }_{ij}}{H_i}-\frac{v_i}{{H_i}^2}\sum _k\frac{\partial H_j}{\partial {\xi }_k}{J}_{kj},J_{ij}=\frac{\partial {\xi }_i}{\partial {\xi }_{j0}},{\Omega }_{ij}=\frac{\partial v_{si}}{\partial {\xi }_{j0}}, \\ \frac{\partial {\Omega }_{ij}}{\partial t}=\frac{\partial }{\partial {\xi }_{j0}}\left\{\sum _k\frac{v_{sk}}{H_i{H}_k}\left(v_{sk}\frac{\partial H_k}{\partial {\xi }_i}-v_{si}\frac{\partial H_i}{\partial {\xi }_i}\right)\right\}+\frac{1}{m}\frac{\partial f_{si}}{\partial {\xi }_{j0}}. \end{gathered}$$ (3.5)

Для заданной модели силового взаимодействия фаз ( f_si — известные функции параметров несущей фазы и скоростей частиц) и фиксированной системы координат (H_i — известные функции координат) дифференцирование выражения в фигурных скобках и f_si по лагранжевым координатам следует проводить по правилам дифференцирования сложных функций. Полученные в результате уравнения (3.2), (3.4)–(3.5) при ξ_i0 = const превращаются в замкнутую систему ОДУ для ξ_i, v_si, J_ij, Ω_ij и конечного соотношения для n_s, позволяющих рассчитывать концентрацию одновременно с расчетом траектории частиц.

В случае очень сложных выражений для силы, действующей на частицы (например, при учете подъемных сил Магнуса, Сэфмана, термофоретической силы и др.), для расчета компонент якобиана вместо дифференцирования и вывода дополнительных уравнений (3.5) может быть применена численная процедура. Например, можно одновременно рассчитывать нескольких соседних траекторий частиц (у которых начальные координаты отличаются на малую величину) и вычислять частные производных $\partial x_i/\partial x_{i0},\partial V_{si}/\partial x_{j0}$ с помощью конечно-разностной аппроксимации. Такая процедура была предложена в работе [40].

4. ВАЖНЫЕ ЧАСТНЫЕ СЛУЧАИ

Установившееся течение несущей фазы. В этом случае траектории среды частиц стационарны, что позволяет выразить часть компонентов якобиана перехода от эйлеровых к лагранжевым переменным через компоненты скорости частиц и существенно сократить число искомых функций при вычислении концентрации дисперсной фазы. Выберем траекторию частиц r = r(r₀, t). В силу ее стационарности, имеем

$\mathbf{r}({\mathbf{r}}_0,t+dt)=\mathbf{r}({\mathbf{r}}_0+d{\mathbf{l}}_0,t),d{\mathbf{l}}_0={\mathbf{V}}_{s0}dt.$

Здесь V_s0 — скорость частиц в точке r₀. Следовательно, производная радиуса-вектора вдоль лагранжева направления l₀ может быть вычислена как

$\frac{\partial \mathbf{r}}{\partial l_0}({\mathbf{r}}_0,t)=\frac{1}{V_{s0}}\frac{\partial \mathbf{r}}{\partial t}=\frac{{\mathbf{V}}_s}{V_{s0}}.$ (4.1)

Отсюда можно выразить три компоненты якобиана преобразования от эйлеровых к лагранжевым переменным через остальные компоненты якобиана и компоненты вектора скорости. Например, при использовании декартовых координат из (4.1) получаем соотношения

$\frac{\partial x_i}{\partial x_0}=\frac{v_{si}}{V_{s0}{n}_1}-\frac{n_2}{n_1}\frac{\partial x_i}{\partial y_0}-\frac{n_3}{n_1}\frac{\partial x_i}{\partial z_0},n_i=\cos ({\mathbf{V}}_{s0},{\mathbf{e}}_i),(x_1=x,x_2=y,x_3=z),(i=\mathrm{1,2,3}).$ (4.2)

Здесь e_i — орты осей координат. После подстановки (4.2) в выражение для якобиана уравнение неразрывности для среды частиц упрощается

$n_s({\mathbf{r}}_0,t)\det \left|\begin{array}{ccc}{u}_s& \partial x/\partial y_0& \partial x/\partial z_0\\ v_s& \partial y/\partial y_0& \partial y/\partial z_0\\ w_s& \partial z/\partial y_0& \partial z/\partial z_0\end{array}\right|=n_{s0}\left({\mathbf{r}}_0\mathrm{,0}\right)u_{s0}.$ (4.3)

Соответственно, порядок системы ОДУ, которую необходимо решать на выбранной траектории частиц, сокращается на шесть уравнений.

Стационарность траекторий частиц позволяет также значительно упростить процедуру нахождения нестационарного поля концентрации в случае, когда граничные условия для концентрации известным образом зависят от времени (n_s0(r₀, t_g) = G(t_g), t_g — “глобальное” время) (например, при исследовании эволюции неоднородностей концентрации в стационарном поле скоростей несущей фазы [43]). В этом случае, как следует из выражения (4.3), мгновенное значение концентрации на выбранной траектории в момент t (здесь t — время движения частицы по траектории от точки r₀) выражается через функцию G и значение концентрации n_st(r₀,t), найденное для стационарного граничного условия n_s0 = 1: 

$n_s(r_0,t)=G(t_g-t)n_{st}(r_0,t). n_s(r_0,t)=G(t_g-t)n_{st}(r_0,t).$

Течения с рождением или исчезновением частиц. В случае заданных источников/стоков количества частиц уравнение неразрывности в лагранжевой форме принимает вид

$\frac{\partial }{\partial t}\left[n_s({\mathbf{r}}_0,t)\det \left(J\right)\right]=j_n\det \left(J\right).$ (4.4)

Здесь j_n — мощность объемного источника/стока числа частиц. Таким образом, в данном случае вместо конечного соотношения (2.1) на выбранной траектории частиц следует решать дополнительное обыкновенное дифференциальное уравнение (4.4).

Эволюция материальных поверхностей и линий. При исследовании эволюции материальных поверхностей или линий удобно вводить поверхностную n_sur или линейную n_l числовую плотность частиц. Закон сохранения числа частиц для элемента поверхности (линии) имеет вид

$n_{sur}d\Sigma =n_{sur0}d{\Sigma }_0,n_{l}dl=n_{l0}d{l}_0.$

Если использовать декартовы координаты и ввести параметризацию поверхности (I) (линии (II)) в начальный момент

$$\begin{gathered} \begin{array}{cc}\left(\text{I}\right)& z_0=s(x_0,y_0),\end{array} \\ \begin{array}{cc}\left(\text{II}\right)& y_0=g\left(x_0\right),\end{array}{z}_0=h\left(x_0\right), \end{gathered}$$

то уравнения неразрывности в лагранжевой форме можно записать в виде

$$\begin{gathered} \begin{array}{cc}\left(\text{I}\right)& n_{sur}\end{array}(x_0,y_0,t)\sqrt{|J{|}^2+|J{|}_1^2+|J{|}_2^2}=n_{sur0}(x_0,y_0\mathrm{,0})\sqrt{1+{s^{\text{'}}}_{x_0}^2+{s^{\text{'}}}_{y_0}^2}, \\ \left|J\right|=\det \left|\begin{array}{cc}\partial x/\partial x_0& \partial y/\partial x_0\\ \partial x/\partial y_0& \partial y/\partial y_0\end{array}\right|,{\left|J\right|}_1=\det \left|\begin{array}{cc}\partial z/\partial x_0& \partial y/\partial x_0\\ \partial z/\partial y_0& \partial y/\partial y_0\end{array}\right|,{\left|J\right|}_2=\det \left|\begin{array}{cc}\partial x/\partial x_0& \partial z/\partial x_0\\ \partial x/\partial y_0& \partial z/\partial y_0\end{array}\right|; \\ \begin{array}{cc}\left(\text{II}\right)& n_l\end{array}(x_0,t)\sqrt{{\left(\frac{\partial x}{\partial x_0}\right)}^2+{\left(\frac{\partial y}{\partial x_0}\right)}^2+{\left(\frac{\partial z}{\partial x_0}\right)}^2}=n_{l0}\left(x_0\mathrm{,0}\right)\sqrt{1+g^{\text{'}}{}_{x_0}^2+h^{\text{'}}{}_{x_0}^2}. \end{gathered}$$

В данные соотношения входит лишь часть компонентов якобиана. Таким образом, при исследовании эволюции материальных поверхностей и линий число ОДУ, которые необходимо решать на выбранной траектории для одновременного вычисления концентрации, существенно сокращается: в случае (I) порядок системы — 18, а в случае (II) — 12.

5. ВЫЧИСЛЕНИЕ ДИФФЕРЕНЦИАЛЬНЫХ ХАРАКТЕРИСТИК ПОЛЕЙ КОНТИНУАЛЬНЫХ ПАРАМЕТРОВ ДИСПЕРСНОЙ ФАЗЫ

Вдоль выбранных траекторий частиц можно вычислять и дифференциальные характеристики континуальных параметров дисперсной фазы, содержащие их производные по эйлеровым координатам. При известном законе деформирования среды частиц в декартовых координатах $\mathbf{r}({\mathbf{r}}_0,t)$ выражения для градиентов концентрации, температуры, дивергенции и ротора скорости среды частиц в эйлеровых и лагранжевых координатах связаны следующим образом:

$$\begin{gathered} \nabla n_s(\mathbf{r},t)=J_T^{-1}{\nabla }_0{n}_s({\mathbf{r}}_0,t),\nabla T_s(\mathbf{r},t)=J_T^{-1}{\nabla }_0{T}_s({\mathbf{r}}_0,t),\text{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}div(}{\mathbf{V}}_s(\mathbf{r},t))=\nabla {\mathbf{V}}_s(\mathbf{r},t)=J_T^{-1}{\nabla }_0{\mathbf{V}}_s({\mathbf{r}}_0,t), \\ {⃐}_s(\mathbf{r},t)=\left[\nabla \times {\mathbf{V}}_s(\mathbf{r},t)\right]=\left[J_T^{-1}{\nabla }_0\times {\mathbf{V}}_s({\mathbf{r}}_0,t)\right],\nabla =\left\{\mathbf{i}\frac{\partial }{\partial x},\mathbf{j}\frac{\partial }{\partial y},\mathbf{k}\frac{\partial }{\partial z}\right\},{\nabla }_0=\left\{\mathbf{i}\frac{\partial }{\partial x_0},\mathbf{j}\frac{\partial }{\partial y_0},\mathbf{k}\frac{\partial }{\partial z_0}\right\}. \end{gathered}$$ (5.1)

Здесь J_T^-1 — матрица, обратная транспонированной матрице Якоби; выражения в левых частях содержат частные производные по эйлеровым переменным $(x,y,z)$, а в правых — по лагранжевым $(x_0,y_0,z_0).$ Правые части в уравнениях для завихренности и дивергенции скорости в системе (5.1) содержат компоненты матриц J_ij и Ω_ij, которые при стоксовом законе сопротивления частиц вычисляют на выбранной траектории частиц из решения обыкновенных дифференциальных уравнений (2.2).

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

$q_s=4\pi \sigma \lambda (T-T_s).$

Дифференцируя третье уравнение в выражении (1.2) по лагранжевым координатам, можно получить векторное уравнение переноса для градиента температуры частиц по лагранжевым переменным:

$\frac{\partial }{\partial t}\left({\nabla }_0{T}_s\right({\mathbf{r}}_0,t\left)\right)=\theta \left\{J_T\nabla T\left[\mathbf{r}\right({\mathbf{r}}_0,t\left)\right]-{\nabla }_0{T}_s\right\},\theta =\frac{6\pi \sigma \lambda }{m{c}_s}.$ (5.2)

Градиент температуры по эйлеровым координатам пересчитывается с использованием второй формулы из выражения (5.1).

Для вычисления градиента числовой плотности частиц на выбранной траектории применим оператор ∇₀ к уравнению неразрывности среды частиц

$n_s({\mathbf{r}}_0,t)\det \left(J\right)=n_{s0}\left({\mathbf{r}}_0\mathrm{,0}\right)$.

В результате получим

${\nabla }_0{n}_s({\mathbf{r}}_0,t)=\frac{{\nabla }_0{n}_{s0}\left({\mathbf{r}}_0\mathrm{,0}\right)-n_s({\mathbf{r}}_0,t){\nabla }_0\det \left(J\right)}{\det \left(J\right)}.$ (5.3)

Видно, что уравнение (5.3) содержит градиент якобиана перехода от эйлеровых переменных к лагранжевым. Его вычисление на выбранной траектории частиц требует знания вторых производных ${\partial }^2{x}_i/\partial x_{0j}\partial x_{0k},$ т. е. компонент гессиана $H_{ijk}$ перехода от эйлеровых к лагранжевым переменным.

При стоксовском законе сопротивления частиц, дифференцируя уравнения (2.2) по лагранжевым координатам и меняя порядок дифференцирования, получим дополнительные уравнения для вычисления компонент гессиана H_ijl на выбранной траектории частиц:

$\frac{\partial H_{ijl}}{\partial t}={\Psi }_{ijl},\frac{\partial {\Psi }_{ijl}}{\partial t}=\frac{1}{\tau }\left(\sum _k\frac{\partial v_i}{\partial x_k}{H}_{kjl}+\sum _{k,m}\frac{{\partial }^2{v}_i}{\partial x_k\partial x_m}{J}_{kj}{J}_{ml}-{\Psi }_{ijl}\right),H_{ijl}=\frac{{\partial }^2{x}_i}{\partial x_{j0}\partial x_{l0}}\begin{array}{cc},& {\Psi }_{ijl}\end{array}=\frac{{\partial }^2{v}_{si}}{\partial x_{0j}\partial x_{0l}}.$ (5.4)

Как и ранее, при использовании уравнений (5.1)–(5.4) совместно с системой (2.1)–(2.2) на выбранной траектории частиц предполагается, что поля температуры и скорости известны (вплоть до вторых частных производных по координатам). Начальные условия $H_{ijl}\left({\mathbf{r}}_0\mathrm{,0}\right),{\Psi }_{ijl}\left({\mathbf{r}}_0\mathrm{,0}\right)$ для уравнений (5.4) в различных задачах должны задаваться в соответствии с граничными условиями на входном участке течения. Следует отметить, что компоненты гессиана могут быть использованы для вычисления относительной деформации малого лагранжева элемента на каустике, где якобиан матрицы J равен нулю. Примеры такого подхода можно найти в работах [44, 45].

Для несжимаемых стратифицированных течений с малоинерционной дисперсной примесью, которая движется фактически со скоростью несущей фазы (τ = 0, V_s = V), якобиан det(J) становится равным единице и процедура вычисления градиента концентрации значительно упрощается. Например, первое соотношение из выражения (5.1) принимает простой вид:

$\nabla n_s(\mathbf{r},t)={J_T}^{-1}(\mathbf{r},t){\nabla }_0{n}_{s0}\left({\mathbf{r}}_0\mathrm{,0}\right).$

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

Рис. 1. Пример расчета ∂ρ/∂x при обтекании сферы вертикально стратифицированным потоком в вертикальной срединной плоскости лагранжевым методом (а) и теневая фотография из работы [47] (б).

Для качественного сравнения на рис. 1б приведена теневая фотография аналогичного поля продольной компоненты градиента плотности при обтекании вихревого кольца стратифицированным потоком из работы [47]. Внешнее течение в таком потоке (за исключением следа) качественно соответствует невязкому обтеканию сферы.

На рис. 2 приведен пример расчета с помощью лагранжева подхода мгновенной картины модуля градиента плотности $|\nabla \rho |$ в горизонтально стратифицированном потоке при потенциальном обтекании круглого цилиндра единичного радиуса. Предполагается, что безразмерный градиент плотности в набегающем потоке изменяется по гармоническому закону $\nabla \rho =\left\{\sin \left(kx\right)\mathrm{,0}\right\}.$

Рис. 2. Мгновенная картина модуля градиента плотности при обтекании цилиндра горизонтально стратифицированным потоком с гармоническими колебаниями плотности.

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

6. ПОЛНЫЙ ЛАГРАНЖЕВ ПОДХОД ДЛЯ РЕШЕНИЯ КИНЕТИЧЕСКОГО УРАВНЕНИЯ БЕССТОЛКНОВИТЕЛЬНОЙ СРЕДЫ С ФАЗОВЫМИ ПЕРЕХОДАМИ

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

В качестве примера рассмотрим случай, когда частицы полидисперсны, размер частиц изменяется в силу фазовых переходов, но частицы одинакового размера имеют одинаковые скорости и температуру. Тогда одночастичная функция распределения частиц, удовлетворяющая кинетическому уравнению (1.1), имеет вид

$f^{\left(1\right)}=N(t,\mathbf{r},\sigma )\delta (\mathbf{c}-{\mathbf{V}}_s)\delta (T-T_s).$

Для простоты выкладок рассмотрим одномерный случай и введем лагранжевы координаты x₀, σ₀ (индекс 0 отмечает значения соответствующих величин в начальный момент). Уравнение (1.1) в лагранжевых переменных примет вид

$$\begin{gathered} N(t,x_0,{\sigma }_0)\det \left(J\right)=N(\mathrm{0,}{x}_0,{\sigma }_0), \\ {J}_{11}=\partial x/\partial x_0\begin{array}{cc},& J_{12}\end{array}=\partial x/\partial {\sigma }_0,J_{21}=\partial \sigma /\partial x_0,J_{22}=\partial \sigma /\partial {\sigma }_0. \end{gathered}$$

Это уравнение позволяет найти функцию N, если значения компонентов якобиана известны. Как и в разд. 2, уравнения для этих компонентов могут быть выведены путем дифференцирования уравнений движения и испарения частицы по лагранжевым координатам. В результате на фиксированной “траектории” в фазовом пространстве получаем систему ОДУ:

$$\begin{gathered} \frac{dx}{dt}=\begin{array}{cc}{u}_s,& \frac{d{u}_s}{dt}\end{array}=\begin{array}{cc}F(x,\sigma ,u_s),& \frac{d\sigma }{dt}\end{array}=\begin{array}{cc}j(x,\sigma ),& F\end{array}=\frac{f_{sx}}{m}; \\ \begin{array}{cc}\frac{d{J}_{11}}{dt}=A_1,& \frac{d{J}_{12}}{dt}\end{array}=\begin{array}{cc}{A}_2,& \frac{d{J}_{21}}{dt}\end{array}=\frac{\partial j}{\partial x}{J}_{11}+\frac{\partial j}{\partial \sigma }\begin{array}{cc}{J}_{21},& \frac{d{J}_{22}}{dt}\end{array}=\frac{\partial j}{\partial x}{J}_{12}+\frac{\partial j}{\partial \sigma }{J}_{22}; \\ \frac{d{A}_1}{dt}=\frac{\partial F}{\partial x}{J}_{11}+\frac{\partial F}{\partial \sigma }\begin{array}{cc}{J}_{21},& \frac{d{A}_2}{dt}\end{array}=\frac{\partial F}{\partial x}{J}_{12}+\frac{\partial F}{\partial \sigma }{J}_{22}. \end{gathered}$$

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

7. БАЗОВАЯ ЭЙЛЕРОВО-ЛАГРАНЖЕВА МОДЕЛЬ СРЕДЫ “ГАЗ / ЖИДКОСТЬ — ИНЕРЦИОННЫЕ ЧАСТИЦЫ” И ПАРАМЕТРЫ ПОДОБИЯ

Сформулируем замкнутую систему уравнений для смешанного эйлерово-лагранжева описания газодисперсных сред с вязкой, сжимаемой несущей фазой. В общем случае в эйлеровых координатах несущая фаза описывается уравнениями Навье–Стокса с источниковыми членами, учитывающими обратное влияние дисперсной фазы, а среда частиц описывается моделью “холодного” континуума с помощью полного лагранжева подхода. Тогда эйлерово-лагранжевы уравнения двухфазной среды в декартовых координатах x, y, z и x₀, y₀, z₀ принимают вид

$$\begin{gathered} \frac{\partial \rho }{\partial t}+\mathrm{div}\left(\rho \mathbf{V}\right)=\mathrm{0,}\rho \left(\frac{\partial \mathbf{V}}{\partial t}+(\mathbf{V}\nabla )\mathbf{V}\right)=-\nabla p+{\nabla }_i{\tau }^{ij}{\mathbf{e}}_j-\sum _k{n}_{sk}{\mathbf{f}}_{sk},c_v\rho \left(\frac{\partial T}{\partial t}+(\mathbf{V}\nabla )T\right)= \\ =-p\text{div}\left(\mathbf{V}\right)+\text{div}(\lambda \nabla T)-\sum _k\left[n_{sk}{q}_{sk}-n_{sk}{\mathbf{f}}_{sk}(\mathbf{V}-{\mathbf{V}}_{sk})\right],n_s(t,{\mathbf{r}}_0)\left|\text{\hspace{0.17em}det}\right(J\left)\right|=n_{s0}\left(\mathrm{0,}{\mathbf{r}}_0\right), \\ \frac{\partial \mathbf{r}(t,{\mathbf{r}}_0)}{\partial t}={\mathbf{V}}_s,m\frac{\partial {\mathbf{V}}_s(t,{\mathbf{r}}_0)}{\partial t}={\mathbf{f}}_s,c_{s}m\frac{\partial T_s(t,{\mathbf{r}}_0)}{\partial t}=q_s,\frac{\partial J_{ij}(t,{\mathbf{r}}_0)}{\partial t}={\Omega }_{ij}, \\ m\frac{\partial {\Omega }_{ij}(t,{\mathbf{r}}_0)}{\partial t}=\frac{\partial f_{si}}{\partial x_{0j}},{\tau }^{ij}=2\mu e^{ij}-\frac{2}{3}\mu {\delta }^{ij}\text{div}\left(\mathbf{V}\right),e^{ij}=\frac{1}{2}\left(\frac{\partial V_i}{\partial r_j}+\frac{\partial V_j}{\partial r_i}\right),J_{ij}=\frac{\partial r_i(t,{\mathbf{r}}_0)}{\partial r_{0j}}, \\ p=\begin{array}{cc}\rho RT,& \frac{\mu }{{\mu }_0}\end{array}=\frac{\lambda }{{\lambda }_0}={\left(\frac{T}{T_0}\right)}^{\kappa },\mathbf{r}=\{x,y,z\}\begin{array}{cc},& {\mathbf{r}}_0\end{array}=\{x_0,y_0,z_0\}. \end{gathered}$$ (7.1)

Здесь δ^ij — символы Кронекера, R — газовая постоянная, e_j — базисные векторы системы координат, κ — показатель степенной зависимости вязкости и теплопроводности несущей фазы от температуры, индекс 0 относится к параметрам, заданным граничными условиями, суммирование по k означает суммирование по числу “складок” среды частиц в областях пересекающихся траекторий. Заранее неизвестные границы этих областей находятся в процессе решения уравнений дисперсной фазы в лагранжевых координатах. В случае несжимаемой несущей фазы вместо уравнения состояния совершенного газа имеем ρ = const, при этом система (7.1) значительно упрощается: уравнение неразрывности несущей фазы примет вид div(V) = 0, а выражение для касательных напряжений τ ^ij = 2μe^ij.

Система (7.1) представляет собой более общую формулировку классической модели “запыленного газа” [14], учитывающую возможность многозначности полей континуальных параметров среды частиц при возникновении “складок” и “сборок” в дисперсной фазе. Для замыкания модели (7.1) требуется конкретизировать выражения для сил f_s и притока тепла к частице q_s, вычислить частные производные $\partial f_{si}/\partial x_{0j}$(дифференцируя $f_{si}\left(\mathbf{V}\right[\mathbf{r}({\mathbf{r}}_0,t)],{\mathbf{V}}_s({\mathbf{r}}_0,t\left)\right]\mathrm{...})$ как сложные функции), а также задать начальные и граничные условия.

Для параметров несущей фазы V, p, T граничные и начальные условия задаются так же, как в случае течения чистой фазы, а для среды частиц, описываемой системами обыкновенных дифференциальных уравнений на выбранных траекториях частиц, следует задавать все необходимые начальные значения V_s, n_s, T_s, J_ij, Ω_ij в области втекания дисперсной фазы. Формально не требуется задавать граничные условия для среды частиц на обтекаемых твердых поверхностях, однако при различных условиях взаимодействия частиц/капель со стенкой (формирование жидкой пленки, регулярный либо диффузный отскок и пр.) следует вводить в рассмотрение пристеночные структуры (пленку жидкости, плотный поверхностный слой дисперсной фазы, слой отраженных частиц).

При формулировке уравнений модели в безразмерном виде, в дополнение к стандартным параметрам подобия для уравнений Навье–Стокса (числам Маха M, Рейнольдса Re, Прандтля Pr, показателю адиабаты γ и показателю степени κ), добавляется ряд параметров, характеризующих дисперсную фазу.

При стоксовском законе сопротивления частиц такими параметрами являются (см., например, работу [43]) относительная массовая концентрация частиц α = mn_s0 / ρ₀, параметр инерционности частиц β, равный отношению макромасштаба рассматриваемой задачи L к характерной длине скоростной релаксации частиц mV₀/6πσμ₀ (в некоторых работах авторы вводят параметр, обратный β и называемый числом Стокса частиц Stk), а также отношение теплоемкости несущей фазы при постоянном давлении к теплоемкости вещества частиц χ = c_p/c_s.

Параметр инерционности частиц характеризует важность двухскоростных эффектов в рассматриваемой задаче. Так, при b → 0 (или Stk → ∞) очень инерционные частицы практически не реагируют на несущую фазу, сохраняя исходные значения своих параметров на рассматриваемом масштабе, а при b → ∞ (или Stk → 0) малоинерционные частицы “вморожены” в среду и движутся без скоростного проскальзывания с несущей фазой. В последнем случае двухконтинуальная модель (7.1) вырождается в односкоростную модель “эффективного газа” [14], характеризуемого измененными (эффективными) значениями параметров подобия M_ef, Re_ef, Pr_ef и γ_ef [14]:

${\mathrm{Re}}_{ef}=(1+\alpha )\mathrm{Re},M_{ef}^2=M^2\frac{(1+\alpha )(1+\alpha \gamma \chi )}{(1+\alpha \chi )},{\mathrm{Pr}}_{ef}=\mathrm{Pr}\frac{(1+\alpha \chi )}{(1+\alpha )},{\gamma }_{ef}=\gamma \frac{1+\alpha \chi }{1+\alpha \gamma \chi }.$

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

В случае более сложных выражений для межфазной силы (при учете нестоксовского режима обтекания частиц, силы тяжести, подъемных сил, заряда частицы и др.) в число параметров подобия добавляются характерные числа обтекания частиц Рейнольдса, Маха и Кнудсена [27, 41], число Фруда, безразмерные коэффициенты при подъемной силе Сэфмана [39], силе Лоренца [23] и др. Дополнительные параметры подобия для несущей фазы могут возникать при задании тепловых граничных условий на неадиабатической стенке, а также при учете других внешних силовых полей.

При конечных массовых концентрациях дисперсной фазы α, численные расчеты по модели (7.1), как правило, осуществляются с помощью итераций по источниковым членам, учитывающим обратное влияние частиц. На каждом шаге итерации при расчете параметров дисперсной фазы параметры несущей фазы считаются известными на эйлеровой сетке. Такие численные алгоритмы использовались, например, в работах [49–51]. При этом в исследовании [51] учитывались также эффекты полидисперсности и испарения дисперсной фазы.

8. ПРИМЕРЫ РАСЧЕТОВ, ИЛЛЮСТРИРУЮЩИЕ ВОЗМОЖНОСТИ И ПРЕИМУЩЕСТВА ПОЛНОГО ЛАГРАНЖЕВА МЕТОДА

За последние два десятилетия полный лагранжев подход приобрел значительную популярность при расчетах распределения концентрации инерционной примеси в газодисперсных потоках применительно к различным приложениям. В ряде работ вычислительного характера (см., например, [51–57]) на основе сравнения различных численных методов подчеркивались значительные преимущества ПЛП по сравнению с другими методами нахождения поля концентрации инерционной дисперсной фазы. В литературе уже появились примеры успешного внедрения вариантов процедуры полного лагранжева подхода в большие вычислительные пакеты OpenFoam и Ansys Fluent [51, 58].

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

Далее приведен ряд характерных примеров течений с зонами пересекающихся траекторий частиц и сингулярностями их числовой плотности вблизи каустик. Распределение концентрации среды частиц в этих течениях рассчитывалось с помощью ПЛП. Нахождение полей концентрации дисперсной фазы в таких течениях с помощью стандартных эйлеровых подходов невозможно в силу многозначности полей континуальных параметров, а методы типа “частицы-в-ячейках” или Монте-Карло требуют расчета намного большего числа траекторий частиц.

Перехлест траекторий инерционных частиц при опрокидывании волны сжатия и в зонах сходящихся линий тока несущей фазы. Самый простой пример течения, в котором возникает каустика в среде частиц — одномерная задача о распространении волны сжатия, т. е. возмущения скорости частиц, созданного затухающим импульсом несущей фазы. В такой задаче “задние быстрые” частицы, получившие большее ускорение, в силу их инерционности со временем догоняют “передние медленные”, и в среде частиц возникает “складка”, на фронте которой появляется каустика с бесконечным значением числовой плотности частиц [1, 7, 17, 31].

В обычном газе, с ненулевым давлением, при опрокидывании волны сжатия вместо перехлеста траекторий частиц среды появляется сильный разрыв — ударная волна. В разреженной дисперсной среде ничто не препятствует возникновению пересечений траекторий частиц и формированию “складки” среды. Характерным примером такого поведения среды частиц могут служить задачи о распространении взрывных волн в запыленном газе, в которых за фронтом расширяющейся и затухающей ударной волны формируется выраженный слой повышенной концентрации частиц [59, 60].

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

Рис. 3. Траектории частиц (пунктир) и изолинии безразмерной концентрации частиц (сплошные линии с цифрами) в потенциальном течении вблизи пробоотборника аэрозоля; стенки пробоотборника — жирные линии; K = U_∞ /U_a = 0.2, Stk = 1.

В работе [61] поля концентрации частиц (в том числе при появлении каустик) находились с помощью полного лагранжева метода в окрестности и внутри плоского тонкостенного пробоотборника при различных значениях определяющих параметров (числа Стокса частиц и отношения скорости набегающего потока к скорости аспирации U_∞/U_a).

Слоистые структуры в сталкивающихся потоках [62, 63]. Вблизи критических точек, возникающих при столкновении противоположно направленных потоков (рис. 4), структура поля концентрации примеси может быть очень сложной. Это обусловлено возникновением локальных зон многократного пересечения траекторий инерционных частиц, связанного с последовательными разворотами траекторий частиц, попадающих во встречный поток. Распределение концентрации частиц в окрестности критической точки достаточно общего вида (сталкиваются два потока под разными углами, потоки имеют различные плотности и вязкости, частицы примеси имеют различные параметры инерционности) исследовано в [62].

Рис. 4. Течение в окрестности “косой” критической точки [62]. Синие линии — линии тока несущей среды, красные — траектории частиц, пунктир — огибающие траектории частиц (каустики).

Рисунки 4 и 5 соответствуют случаю, когда частицы присутствуют только в верхнем потоке (индекс 1), причем потоки имеют одинаковые плотности, но разные коэффициенты вязкости: вязкость нижнего потока в десять раз превосходит вязкость верхнего потока.

Рис. 5. Распределения вертикальных скоростей несущей фазы (синие линии), частиц (красные линии) (а) и безразмерной концентрации частиц (б) в окрестности критической точки [62].

В работе [62] поле скорости несущей фазы найдено численно из решения автомодельной задачи для окрестности неортогональной критической точки. Масштабы скорости и длины: V = (μ₁C₁/ρ₁)^1/2, L = (μ₁/C₁ρ₁)^1/2, где C₁ — “скорость растекания” (продольный градиент скорости вдали от критической точки) верхнего потока. Инерционные частицы неоднократно пересекают контактную поверхность, при этом возникают многослойные “складки” в среде частиц, а на огибающих траекторий (каустиках) и на контактной поверхности возникают узкие зоны накопления частиц. Таким образом, в поле суммарной концентрации формируются слоистые структуры, расчет которых другими методами вряд ли возможен.

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

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

Формирование каустик и зон накопления частиц при обтекании затупленных тел запыленным потоком. Другой наглядный пример формирования каустик и локальных зон накопления частиц — обтекание затупленного тела запыленным газом. Публикаций, в которых использовался ПЛП для расчета распределения концентрации инерционных частиц в задачах двухфазного обтекания тел до- и сверхзвуковыми газодисперсными потоками, довольно много, здесь ограничимся упоминанием лишь трех трудов [32, 40, 64], в которых продемонстрированы основные особенности поведения дисперсной фазы вблизи поверхности обтекаемого тела.

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

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

Во втором случае, при зеркальном отражении частиц от обтекаемой поверхности, на огибающей траектории отраженных частиц возникает каустика с интегрируемой сингулярностью концентрации частиц [30]. Внутри слоя отраженных частиц в различных сечениях могут возникать разрывы концентрации частиц, соответствующие границам зон отраженных частиц, содержащим различное число “складок”, т. е. различное количество отраженных траекторий, приходящих в одну точку пространства. Количество “складок” увеличивается с приближением к оси симметрии.

Указанные особенности распределения концентрации частиц продемонстрированы на рис. 6 из работы [40] (Stk = 1/β — число Стокса частиц), в которой авторы с помощью полного лагранжева подхода рассчитывали поля параметров дисперсной фазы при сверхзвуковом обтекании плоского цилиндра невязким запыленным газом. Параметры несущей фазы при этом рассчитывались на эйлеровой сетке.

Рис. 6. Картины траекторий частиц (левые рисунки) и изолиний безразмерной концентрации частиц (правые рисунки) при обтекании цилиндра невязким сверхзвуковым потоком [64]. Верхние рисунки соответствуют режиму отсутствия инерционного осаждения частиц, а нижние — режиму с отражением частиц.

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

В работе [32] на основе вычисления вероятности столкновения падающих и отраженных частиц на оси симметрии было показано, что для частиц размером более 1 мкм и сферы радиуса 1 м столкновениями частиц можно пренебречь, когда объемная доля частиц в набегающем потоке не превосходит величин ~ 10^–5.

Фокусировка частиц за точкой пересечения ударных волн. Как отмечалось ранее, возникновения зон пересекающихся траекторий инерционных частиц следует ожидать в локальных областях сходящихся линий тока несущей фазы. Расчеты траекторий и поля концентрации частиц в окрестности точки взаимодействия косых стационарных ударных волн показывают [41], что как при регулярном, так и при маховском взаимодействии волн за точкой взаимодействия происходит аэродинамическая фокусировка инерционных частиц, формируется узкая область, в которой траектории частиц пересекаются, а концентрация частиц резко возрастает. Примеры расчета траекторий и профилей суммарной концентрации частиц, формирующихся в области, удаленной вниз по потоку, приведены на рис. 7. Число Маха волн равно пяти, для симметричного регулярного взаимодействия угол между волнами равен 120°. Концентрация частиц отнесена к ее значению в невозмущенном потоке, единица длины — длина релаксации частицы при стоксовском законе сопротивления, вычисленная для параметров невозмущенного потока.

Рис. 7. Картины траекторий частиц (верхние рисунки) и соответствующие профили суммарной концентрации частиц (нижние рисунки), формирующиеся вдали за точкой взаимодействия волн; y_tr на последнем рисунке отсчитывается по нормали к контактной поверхности.

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

Существует оптимальный режим фокусировки частиц, при котором конечный объем дисперсной фазы “схлопывается” в поверхность, создавая коллимированный пучок частиц [41]. Формирование таких пучков может быть использовано в технологиях газодисперсного напыления, но в некоторых случаях этот эффект может представлять серьезную опасность. Например, в работе [66] проведено численное исследование сверхзвукового обтекания плоского круглого цилиндра слабозапыленным потоком в условиях, когда на головную ударную волну падает косой скачок уплотнения (рис. 8, 9). Условия течения несущей фазы, соответствовали экспериментальной работе [67], число Маха набегающего потока равно шести.

Рис. 8. Теневая фотография из работы [67] и расчеты распределения числа Маха в ударном слое из работы [66] для случая, когда на боковую поверхность цилиндра приходит контактный разрыв.

Рис. 9. Типичные картины траекторий частиц различной инерционности и поля безразмерной концентрации дисперсной фазы, соответствующие расчетам из работы [66] (отраженные частицы не учитывались).

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

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

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

Фокусировка частиц в сдвиговых потоках за счет действия подъемных сил. Еще один важный механизм фокусировки частиц связан с действием на частицы боковых (подъемных) сил. Их возникновение обусловлено двумя основными механизмами — неоднородностью потока, обтекающего частицу, и возможным вращением частиц. “Шнурование” частиц (т. е. формирование сгустков частиц) на оси дозвуковой струи, истекающей из узкого канала, отмечалось еще в экспериментах [69]. Этот эффект обычно связывают с действием подъемной силы Магнуса, обусловленной большой скоростью вращения, которую достаточно крупные инерционные частицы приобретают за счет столкновений со стенками канала, из которого истекает струя.

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

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

${\mathbf{F}}_{Saf}=6.46{\sigma }^2\sqrt{\mu \rho \left|\frac{\partial u}{\partial y}\right|}(u-u_s)\text{\hspace{0.17em}sign}\left(\frac{\partial u}{\partial y}\right)\mathbf{j}.$ (8.1)

Здесь u и u_s — продольные скорости газа и частиц, j — единичный вектор, направленный перпендикулярно потоку. Как видно из выражения (8.1), направление действия силы Сэфмана зависит от направления градиента поперечной скорости потока, обтекающего частицу в системе отсчета, связанной с частицей.

В работах [39] и [50] исследована структура пограничного слоя запыленного газа, формирующегося за ударной волной, движущейся с постоянной скоростью вдоль плоской поверхности. В межфазном обмене импульсом, кроме силы аэродинамического сопротивления, учитывалась подъемная сила Сэфмана (8.1). Частицы могут поступать в пограничный слой как из внешнего потока (случай движения волны в запыленном газе вдоль твердой стенки на рис. 10а), так и с подстилающей поверхности (случай движения волны над эродирующим слоем дисперсного осадка на рис. 10б). В обоих случаях в пристеночной зоне за ударной волной возникают области многократного пересечения траекторий частиц и протяженные зоны повышенной концентрации дисперсной фазы.

Рис. 10. Траектории частиц (линии 2) и каустики (жирные линии 3) в системе отсчета, связанной с движущейся ударной волной вдоль стенки [39] (граница пограничного слоя — линия 1); (а) — частицы поступают в пограничный слой из внешнего потока; (б) — частицы поступают с поверхности разрушающегося слоя дисперсной фазы, лежащего на стенке.

На рис. 10а, 10б приведены типичные примеры расчетов “фонтанирующих” и “сальтирующих” траекторий частиц в движущейся системе отсчета, связанной с фронтом волны. Продольная координата отнесена к длине скоростной релаксации частиц, а поперечная — к характерной толщине пограничного слоя на этой длине. На каустиках (жирные линии 3) концентрация частиц неограниченно, но интегрируемым образом возрастает, а под нижней ветвью каустики на рис. 10а концентрация частиц равна нулю.

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

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

1) частицы монотонно мигрируют по направлению к стенкам канала;

2) частицы движутся по направлению к оси симметрии канала (рис. 11а);

3) траектории частиц многократно пересекают ось симметрии канала с постепенно уменьшающейся амплитудой (рис. 11б).

Рис. 11. Типичные траектории частиц различной инерционности в системе отсчета, связанной с ударной волной, движущейся в узком канале (поперечная координата отнесена к радиусу канала, продольная — к длине скоростной релаксации частиц).

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

Для каналов диаметром порядка 10^–3 м и слабых ударных волн эффект фокусировки имеет место для частиц с диаметром порядка 10^–5 м.

Фокусировка частиц за счет силы Сэфмана происходит также при стационарном течении аэрозоля с несжимаемой вязкой несущей фазой в профилированном сужающемся микроканале. На рис. 12 приведен пример типичного поведения траекторий инерционных частиц в осесимметричном сужающемся канале, при этом канал можно спрофилировать таким образом, чтобы на выходе получить идеально коллимированный пучок частиц. Подобный эффект фокусировки и создания коллимированных пучков частиц используется в технологиях прямого нанесения вещества (direct-write) [74].

Рис. 12. Типичная картина траекторий частиц в цилиндрических координатах в осесимметричном сужающемся микроканале. Каустики и области накопления частиц формируются на огибающих траектории частиц.

Еще один пример исследования фокусировки частиц под действием боковых сил сдвиговой природы содержится в работе [75], где рассмотрено горизонтальное ламинарное течение разреженной суспензии в вертикальной ячейке Хеле–Шоу. Течения такого типа представляют интерес для описания поведения частиц проппанта в трещинах гидроразрыва нефтяных пластов.

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

В данном случае, в отличие от решения Сэфмана [71], частица движется перпендикулярно сдвигу скорости несущей фазы. Для такого течения коэффициент подъемной силы был вычислен в [76]. В работе [75] был найден характерный продольный масштаб длины, на котором частицы мигрируют поперек течения на расстояние порядка полуширины канала. Эволюция полей числовой концентрации и скорости среды частиц вдоль ячейки исследована с помощью полного лагранжева метода. В зависимости от параметра инерционности частиц обнаружены различные режимы поперечной миграции: с пересечением и без пересечения частицами средней линии канала.

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

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

Каустики в стационарных вихревых течениях. Движение одиночных инерционных включений в стационарных модельных вихревых течениях, поля скорости которых заданы аналитическими формулами, исследовали во многих публикациях (см., например, работы [77–79] и обзор [80]). В [77] показано, что легкие включения (например, пузырьки газа в жидкости) под действием сил Архимеда и присоединенных масс собираются в центре вихрей. Более тяжелые частицы под действием центробежной силы, наоборот, стремятся покинуть центр вихря и собраться на его периферии [78].

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

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

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

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

Расчеты показали, что в случае частиц, превосходящих по плотности несущую фазу вдали от оси вихря, существует разделяющая (“критическая”) высота начала трубки тока дисперсной фазы. Характерное поведение траекторий частиц, начинающихся на внешней границе расчетной области выше и ниже этой высоты в отсутствие силы тяжести, показано на рис. 13а и 13б соответственно.

Рис. 13. Типичные восходящие (а) и нисходящие (б) траектории “тяжелых” частиц в течении типа торнадо [84].

Траектории из некоторого конечного по толщине слоя, стартующие выше критического уровня, приближаются к поверхности накопления частиц, имеющей форму осесимметричной “чаши”, внутри которой частиц нет. Траектории, стартующие ниже критического уровня, стремятся к основанию вихря, при этом сечения соответствующих трубок тока частиц плоскостью (r, z) многократно пересекаются между собой (рис. 14).

Рис. 14. Картина трубок тока несущей фазы (синие пунктирные линии), трубок тока среды “тяжелых” частиц (зеленые сплошные линии), каустик (зеленые пунктирные линии) и изолиний величины, обратной концентрации дисперсной фазы (разноцветные сплошные линии с числами), в плоскости r, z [84].

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

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

Рис. 15. Формирование спиралевидной верхней части “чашеобразной” зоны накопления частиц в течении типа торнадо [84]. Зеленые сплошные линии — трубки тока частиц, цветные пунктирные линии с числами — изолинии величины, обратной концентрации дисперсной фазы.

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

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

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

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

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

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

В рамках использованной модели был обнаружен эффект скопления пыли в тонких слоях, расположенных на расстоянии 0.03–1 а.е. (для частиц 0.41 мкм) и 1–10 а.е. (для частиц 0.17 мкм) по обе стороны от токового слоя. Эти плотные слои в распределении пыли, совпадающие с каустиками на границах зон пересекающихся траекторий частиц, формируются, в основном, из-за фокусирующего воздействия солнечного магнитного поля.

Пример расчета траекторий движения частиц межзвездной пыли размером 0.41 мкм в срединной плоскости, перпендикулярной токовому слою, показан ни рис. 16а. Трубки тока частиц имеют множественные пересечения, из-за чего формируются “складки”, на границах которых (каустиках) числовая плотность частиц принимает бесконечно большие значения. Структура взаимного расположения каустик очень сложна. Пример картины зон накопления частиц в срединной плоскости, перпендикулярной токовому слою, показан на рис. 16б. Обнаруженный эффект формирования плотных слоев космической пыли в гелиосфере представляет интерес при планировании будущих запусков космических аппаратов для исследования структуры гелиосферы.

Рис. 16. Траектории (а) и каустики (б) космической пыли в срединной плоскости гелиосферы, перпендикулярной токовому слою магнитного поля. Солнце находится в начале координат, координаты измеряются в астрономических единицах.

В работе [86] в гидродинамическом приближении рассматривается центробежная сепарация инородных включений (частиц) во вращающемся сферическом объеме самогравитирующей жидкой или газообразной среды. Отдельно рассмотрен случай наличия твердого ядра у вращающегося объема. В данном случае существенное влияние на движение частиц оказывают переменная по радиусу сферы сила гравитации, а также центробежная сила и сила Кориолиса.

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

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

В общем случае картины распределения дисперсной фазы достаточно сложны, они зависят от трех безразмерных параметров, характеризующих инерционность частиц, вклад центробежной силы и вклад силы гравитации. Рисунок 17 из работы [86] иллюстрирует формирование складок в экваториальной плоскости при различном вкладе центробежной силы для частиц с одинаковыми инерционными свойствами. На рис. 18 показана типичная пространственная траектория инерционной частицы и сечения трубок тока дисперсной фазы вертикальной плоскостью. Огибающие трубок тока частиц (каустики) представляют собой вложенные друг в друга эллиптические поверхности, число которых возрастает с приближением к оси симметрии. Внутри области накопления частиц в экваториальной плоскости формируются концентрические окружности с неограниченно растущей числовой плотностью частиц. Для “легких” включений, перемещающихся к периферии объема, концентрация максимальна вблизи оси симметрии.

Рис. 17. Траектории частиц одинаковой инерционности в экваториальной плоскости при различных вкладах центробежной силы; (а) — меньший вклад центробежной силы, (б) — больший вклад центробежной силы.

Рис. 18. Типичная пространственная траектория инерционной частицы (сплошная линия) (а) и сечения трубок тока частиц вертикальной плоскостью (б) из [86].

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

9. КОМБИНИРОВАННЫЙ БЕССЕТОЧНЫЙ МЕТОД РАСЧЕТА ВИХРЕВЫХ НЕСТАЦИОНАРНЫХ ГАЗОДИСПЕРСНЫХ ТЕЧЕНИЙ

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

Для преодоления проблемы пересчета параметров фаз в ряде работ [87–89] был предложен полностью лагранжев (бессеточный) метод расчета параметров обеих фаз. Указанный метод основан на сочетании и модификации двух лагранжевых подходов: вихревого метода для решения уравнений Навье–Стокса/Эйлера для несущей фазы и полного лагранжева метода для расчета параметров дисперсной среды.

В предыдущие годы активно развивались бессеточные лагранжевы подходы к численному моделированию течений чистой жидкости. Это, в частности, методы гидродинамики сглаженных частиц (smoothed-particle hydrodynamics, SPH) [90–91], в которых используют стандартные переменные скорость–давление, а также различные варианты вихревых методов, использующие переменные скорость–завихренность и обобщающие методы дискретных вихрей на случай вязкой среды (см., например, книги [92–93], статьи [94–97] и цитированную в них литературу).

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

В развиваемом в работах [87–89] бессеточном методе расчета параметров вязкой несжимаемой несущей фазы используют вариант вихревого метода, основанный на введении “диффузионной скорости завихренности” и дивергентой форме записи уравнения переноса завихренности в вязкой среде [97–98]. Для нахождения параметров инерционной дисперсной фазы используется полный лагранжев метод.

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

Следует отметить, что в литературе появляются публикации (см., например, работы [98–99] и цитируемую в них литературу), посвященные развитию методов, сочетающих различные варианты вихревых методов для несущей фазы и лагранжев траекторный подход для моделирования примеси. Однако в этих работах рассчитывают лишь траектории частиц и отсутствует аккуратный расчет поля концентрации примеси, учитывающий пересечения траекторий частиц и существенную деформацию элементарного фазового объема дисперсной среды.

Лагранжев вихревой метод для несущей фазы [87–89]. В качестве основы математического описания двухфазной среды используется модель (7.1) с несжимаемой вязкой несущей фазой. В межфазном взаимодействии учитываются силы Стокса, присоединенных масс, Архимедова сила и сила тяжести. Если не оговорено противное, массовая концентрация примеси считается малой, и влиянием частиц на течение несущей фазы пренебрегается. Сделанные предположения позволяют решать задачу о расчете вихревого течения несущей фазы независимо от течения дисперсной среды.

Для расчета параметров несущей фазы используется вариант вихревого метода на основе введения диффузионной скорости завихренности [94, 95]. Рассматриваются нестационарные плоские или осесимметричные течения без закрутки в неограниченной области. В этом случае имеется единственная ненулевая компонента ротора скорости ω (называемая далее завихренностью), направленная перпендикулярно плоскоcти течения в плоском случае и по направлению изменения азимутального угла j в осесимметрчном случае. В переменных скорость–ω уравнения Навье–Стокса переписываются в дивергентной форме в декартовых r = (x, y, z) (I) либо цилиндрических координатах $\mathbf{r}=\left(r,\phi ,z\right)$ (II), где ось симметрии совпадает с осью Oz:

$$\begin{gathered} \frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}=\mathrm{0,}\frac{\partial \omega }{\partial t}+\frac{\partial }{\partial x}\left[\left(u_x+u_{dx}\right)\omega \right]+\frac{\partial }{\partial y}\left[\left(u_y+u_{dy}\right)\omega \right]=0\left(I\right), \\ \frac{\partial }{\partial r}\left(u_{r}r\right)+\frac{\partial }{\partial z}\left(u_{z}r\right)=\mathrm{0,}\frac{\partial \omega }{\partial t}+\frac{\partial }{\partial r}\left[\left(u_r+u_{dr}\right)\omega \right]+\frac{\partial }{\partial z}\left[\left(u_z+u_{dz}\right)\omega \right]=0\left(II\right). \end{gathered}$$   (9.1)

Здесь $\mathbf{v}=(u_x,u_y\mathrm{,0})\begin{array}{c}\end{array}\left(I\right)$ и $\mathbf{v}=\left(u_r\mathrm{,0,}{u}_z\right)\begin{array}{c}\end{array}\left(II\right)$ — конвективная скорость несущей фазы в декартовых и цилиндрических координатах; ${\mathbf{v}}_d=(u_{dx},u_{dy}\mathrm{,0})\begin{array}{c}\end{array}\left(I\right)$ и ${\mathbf{v}}_d=\left(u_{dr}\mathrm{,0,}{u}_{dz}\right)\begin{array}{c}\end{array}\left(II\right)$ — соответствующая диффузионная скорость завихренности, которая в рассматриваемых случаях имеет вид

$$\begin{gathered} {\mathbf{v}}_d=-\frac{1}{\mathrm{Re}}\frac{\nabla \omega }{\omega },\begin{array}{c}\\ \end{array}\omega =\left(\frac{\partial u_x}{\partial y}-\frac{\partial u_y}{\partial x}\right),\begin{array}{c}\\ \end{array}\nabla \omega =\left(\frac{\partial \omega }{\partial x},\frac{\partial \omega }{\partial y}\mathrm{,0}\right)\begin{array}{c}\\ \end{array}\left(I\right), \\ {\mathbf{v}}_d=-\frac{1}{\mathrm{Re}}\frac{\nabla \left(r\omega \right)}{r\omega },\begin{array}{c}\\ \end{array}\omega =\left(\frac{\partial u_r}{\partial z}-\frac{\partial u_z}{\partial r}\right),\begin{array}{c}\\ \end{array}\nabla \left(r\omega \right)=\left(\frac{\partial \left(r\omega \right)}{\partial r},\begin{array}{c}\\ \end{array}\mathrm{0,}\begin{array}{c}\\ \end{array}\frac{\partial \left(r\omega \right)}{\partial z}\right)\begin{array}{c}\\ \end{array}\left(II\right). \end{gathered}$$ (9.2)

В отличие от идеальной жидкости, в вязкой завихренность не является “вмороженной” в среду, а переносится со скоростью, равной сумме скорости среды и диффузионной скорости. Здесь и далее все уравнения записаны в безразмерном виде, для масштабирования использованы некоторые характерные значения скорости U, длины L, времени L/U, завихренности U/L и циркуляции UL. Число Рейнольдса имеет вид Re = LUr/m. Для конкретных примеров характерные масштабы длины и скорости будут пояснены далее.

Для конкретных начальных распределений завихренности поле завихренности в текущий момент рассчитываем из уравнения переноса завихренности (9.1). При найденной завихренности конвективная скорость v в плоском течении может быть восстановлена из интеграла Био–Савара в форме

$\mathbf{v}\left(R\right)=\frac{1}{2\pi }\underset{\Omega }{\int }\omega \left(r\right){\mathbf{e}}_z\frac{\left(\mathbf{R}-\mathbf{r}\right)}{{\left|\mathbf{R}-\mathbf{r}\right|}^2}dS.$ (9.3)

Здесь R — радиус-вектор рассматриваемой точки течения, Ω — область ненулевой завихренности, dS — дифференциал площади в этой области.

В случае осесимметричного течения без закрутки из формулы Био–Савара имеем следующие выражения для компонент скорости [100]:

$u_r\left(r\right)=-\frac{1}{2\pi r}\iint \omega (R,Z)f^{\text{'}}\left(k\right)\frac{\partial k}{\partial z}{\left(rR\right)}^{1/2}dRdZ,$        (9.4)

$u_z\left(r\right)=\frac{1}{2\pi r}\iint \omega (R,Z)\left(\frac{f\left(k\right)}{2r}+f^{\text{'}}\left(k\right)\frac{\partial k}{\partial r}\right){\left(rR\right)}^{1/2}dRdZ,$

$f\left(k\right)=\left(\frac{2}{k}-k\right)K\left(k\right)-\frac{2}{k}E\left(k\right), k(\mathbf{r},\mathbf{R})={\left(\frac{4rR}{{(r+R)}^2+{(z-Z)}^2}\right)}^{1/2}.$

Здесь радиус-векторы r и R рассматриваются на плоскости и имеют координаты (r, z) и (R, Z). Интегрирование по Z ведется от минус до плюс бесконечности, а по R — от нуля до бесконечности, K(k) и E(k) — полные эллиптические интегралы первого и второго рода:

$K\left(k\right)=\text{}\underset{0}{\overset{\pi /2}{\int }}{(1-k^2{\cos }^2\theta )}^{-1/2}d\theta ,E\left(k\right)=\text{}\underset{0}{\overset{\pi /2}{\int }}{(1-k^2{\cos }^2\theta )}^{1/2}d\theta .$

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

Опишем дискретный вариант процедуры нахождения параметров несущей фазы в рамках бессеточного вихревого метода. Область начальной ненулевой завихренности W в плоскости (x, y) (в плоском случае) или (r, z) (в осесимметричном случае) разбивается на совокупность N малых лагранжевых элементов с площадями ∆S_i и координатами центров (x_i, y_i) и (r_i, z_i) соответственно (i = 1,…,N).

В осесимметричном случае каждому такому элементу в пространстве можно поставить в соответствие тороидальное вихревое кольцо, характеризуемое радиусом r_i, координатой вдоль оси симметрии z_i, малым радиусом сечения $a_i\ll r_i$ и циркуляцией скорости  равной потоку завихренности через ∆S_i. Радиус сечения кольца удовлетворяет соотношению πa_i² = ∆S_i.

В силу дивергентной формы уравнения переноса завихренности в (9.1) вихревые лагранжевы элементы, движущиеся с суммарной скоростью v + v_d, сохраняют поток завихренности, а значит, и свое значение циркуляции скорости. Поэтому циркуляции выбранных элементов  постоянны, при этом $\underset{i}{\Sigma }{\gamma }_i=\Gamma ,$ где Г — суммарная циркуляция области завихренности W. В лагранжевых переменных t и ${\mathbf{r}}_{i0}={\mathbf{r}}_i(t=0)$ уравнения движения центров выбранных малых вихревых элементов можно записать в виде

$\frac{\partial {\mathbf{r}}_i\left(t,{\mathbf{r}}_{i0}\right)}{\partial t}=\mathbf{v}\left(t,{\mathbf{r}}_i\right)+{\mathbf{v}}_d\left(t,{\mathbf{r}}_i\right).$ (9.5)

Здесь r_i — радиус вектор i-го лагранжева элемента в плоскости (x, y) (в плоском случае) или (r, z) (в осесимметричном случае). Значения компонент скорости среды в точках r_i находятся из дискретного аналога интеграла Био–Савара. В плоском случае (9.3) имеем

$\mathbf{v}\left({\mathbf{r}}_i\right)=\frac{1}{2\pi }\sum _{k=1}^N{\gamma }_k{\mathbf{e}}_z\frac{{\mathbf{r}}_i-{\mathbf{r}}_k}{|{\mathbf{r}}_i-{\mathbf{r}}_k{|}^2}.$ (9.6)

В осесимметричном случае из выражения (9.4) получаем

$u_r\left({\mathbf{r}}_i\right)=-\frac{1}{2\pi r_i}\sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^N{\gamma }_j{f}^{\text{'}}\left(k_{ij}\right)\frac{\partial k_{ij}}{\partial z_i}\sqrt{r_i{r}_j}, k_{ij}={\left(\frac{4r_i{r}_j}{{(r_i+r_j)}^2+{(z_i-z_j)}^2}\right)}^{1/2}\text{},$  (9.7)

$u_z\left({\mathbf{r}}_i\right)=\frac{{\gamma }_i}{4\pi r_i}\left(\ln \frac{8r_i}{a_i}-\frac{1}{4}\right)+\frac{1}{2\pi r_i}\sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^N{\gamma }_j\left(\frac{f\left(k_{ij}\right)}{2r_i}+f^{\text{'}}\left(k_{ij}\right)\frac{\partial k_{ij}}{\partial r_i}\right)\sqrt{r_i{r}_j}.$

Здесь первое слагаемое в выражении $u_z\left({\mathbf{r}}_i\right)$ отвечает за самоиндуцированную скорость i-го вихревого кольца [101].

Для нахождения завихренности и ее градиента по эйлеровым координатам, необходимых для вычисления диффузионной скорости ${\mathbf{v}}_d$ (9.5) в произвольной точке r, используется интерполяция, основанная на идеях метода SPH [90–92].

$\omega \left(\mathbf{r}\right)=\int \omega \left(\mathbf{R}\right)\delta \left(\mathbf{r},\mathbf{R}\right)d\mathbf{R}\approx \int \omega \left(\mathbf{R}\right){\delta }_{\epsilon }\left(\mathbf{r},\mathbf{R}\right)d\mathbf{R}\approx \sum _{i=1}^N{\gamma }_i{\delta }_{\epsilon }\left(\mathbf{r},{\mathbf{r}}_i\right),$ (9.8)

$\nabla \omega \left(r\right)\approx \sum _{i=1}^N{\gamma }_i\nabla {\delta }_{\epsilon }\left(\mathbf{r},{\mathbf{r}}_i\right).$ (9.9)

Здесь d — дельта-функция, δ_e — приближение дельта-функции “шапочкой”, обеспечивающее слабую сходимость к дельта-функции при e → 0. Для улучшения точности аппроксимации значения малого параметра ε_i выбирают различными для окрестностей различных вихревых элементов (с индексом i), поскольку они зависят от расстояния до центра ближайшего элемента. Вид “шапочки” может быть выбран достаточно произвольно при соблюдении необходимого баланса между громоздкостью вычислений и требуемой точностью. В различных расчетах использовались “шапочки” как первого порядка аппроксимации

${\delta }_{\epsilon }\left(\mathbf{r},{\mathbf{r}}_i\right)=\frac{1}{{\epsilon }_i}\exp \left(-\frac{\left|\mathbf{r}-{\mathbf{r}}_i\right|}{{\epsilon }_i}\right),$

так и более высоких порядков, например

${\delta }_{\epsilon }\left(\mathbf{r},{\mathbf{r}}_i\right)=\frac{1}{3\pi {\epsilon }_i^2}\left(4-\frac{{\left|\mathbf{r}-{\mathbf{r}}_i\right|}^2}{{\epsilon }_i^2}\right)\exp \left(-\frac{{\left|\mathbf{r}-{\mathbf{r}}_i\right|}^2}{{\epsilon }_i^2}\right).$

В осесимметричном случае дополнительно также используются соотношения для малого переменного радиуса сечения вихревого элемента $a_i=a_i\left(t\right):a_i^2{\omega }_i=a_{i0}^2{\omega }_{i0}={\text{const}}_i.$  Эти соотношения получены из условия постоянства потока завихренности через выбранное сечение лагранжева вихревого элемента, центр которого движется с суммарной скоростью v + v_d. В соответствии с рекомендациями [92] величина “носителя” шапочки ε_i берется пропорциональной расстоянию до ближайшего вихревого элемента l_i(t) : ε_i(t) = Cl_i(t), где C — константа больше единицы.

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

Описанный ранее вихревой метод сводит задачу решения уравнений Навье–Стокса для нахождения поля скорости несущей фазы к решению системы обыкновенных дифференциальных уравнений (9.5) для координат вихревых элементов. С учетом постоянства циркуляций скорости для указанных элементов γ_i из выражения (9.9) (при заданном выражении для “шапочки”) получаются явные выражения для диффузионной скорости завихренности, а скорость несущей фазы на каждом шаге по времени восстанавливается с помощью вычисления конечных сумм типа (9.6)–(9.7). Технические детали реализации алгоритма можно найти в работах [87–89].

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

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

Отметим важное преимущество предлагаемого комбинированного полностью лагранжева метода по сравнению со стандартными лагранжево-эйлеровыми методами, в которых поле скорости несущей фазы рассчитывается в фиксированных узлах эйлеровой сетки. В общем случае узлы сетки не совпадают с текущими положениями дисперсных частиц ${\mathbf{r}}_s$, поэтому в точках $\mathbf{r}={\mathbf{r}}_s$ необходимо интерполировать значения требуемых функций с эйлеровой сетки. Это приводит к дополнительной погрешности, порядок которой зависит от порядка интерполяционного многочлена.

Часто используется квадратичная интерполяция (см., например, работу [53]). В предложенном комбинированном бессеточном подходе значения скоростей в любой точке находятся без внесения дополнительной погрешности.

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

10. ПРИМЕРЫ ПРИМЕНЕНИЯ КОМБИНИРОВАННОГО ЛАГРАНЖЕВА ПОДХОДА ДЛЯ РАСЧЕТА НЕСТАЦИОНАРНЫХ ВИХРЕВЫХ ДВУХФАЗНЫХ ТЕЧЕНИЙ

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

$\omega =\frac{\mathrm{Re}}{4\pi t}\exp \left(-\mathrm{Re}\frac{r^2}{4t}\right),v_{\phi }=\frac{1}{2\pi r}\left[1-\exp \left(-\mathrm{Re}\frac{r^2}{4t}\right)\right].$

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

$l_v=\sqrt{\frac{m\Gamma }{6\pi \sigma \mu }},\begin{array}{c}\end{array}U=\frac{\Gamma }{l_v},\begin{array}{c}\end{array}T=\frac{l_v^2}{\Gamma },\begin{array}{c}\end{array}\Omega =\frac{\Gamma }{l_v^2}.$

Здесь Γ — циркуляция скорости; l_v — длина скоростной релаксации частицы. В межфазном обмене импульсом, кроме силы Стокса, учитываются силы присоединенных масс, Архимеда и сила тяжести (1.3). Задача зависит от трех безразмерных параметров: числа Рейнольдса Re = Γρ/μ, числа Фруда Fr² = U ²/gl_v и отношения плотностей фаз η = ρ/ρ_s0.

В типичном расчете участвует порядка 10³ вихревых доменов и столько же пробных траекторий частиц дисперсной фазы. Расчет начинается от безразмерного момента t = 1, начальная безразмерная концентрация частиц однородна в конечном круге и равна единице. В расчетах было получено хорошее соответствие теоретических и численных значений скорости и производных скорости несущей фазы по времени и пространству, а также завихренности. Сравнивались значения концентрации частиц, вычисленные на основе полного лагранжева подхода с использованием точных формул для параметров несущей фазы и расчетов методом вихревых доменов. Примеры расчетов концентрации и траекторий частиц приведены на рис. 19–21.

Рис. 19. Вихревые домены (серые точки) и пробные частицы дисперсной фазы (цветные точки, цвет соответствует локальной концентрации частиц) в вихре Ламба–Озеена для “тяжелых” стоксовских частиц в моменты t = 2, 3, 4 при Re = 100, Fr = 1, η = 0.

Видно, что за счет силы тяжести облако оседает и становится несимметричным, при этом за счет центробежной силы на границе облака формируется область накопления частиц. Поведение “легких частиц”, для которых важны силы Архимеда и присоединенных масс, в вихре оказывается гораздо более сложным. Уже на начальном этапе в облаке формируются зоны пересекающихся траекторий частиц и каустики (рис. 20).

Рис. 20. Распределения концентрации “легких” частиц в начальные моменты времени при Re = 100, Fr = ∞, η = 1.2. Сплошные линии соответствуют использованию аналитических формул для параметров несущей фазы, крестики — расчету с помощью метода вихревых доменов.

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

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

Рис. 21. Типичные траектории “легких” частиц в вихре Ламба–Озеена при Re = 100, Fr = ∞, η = 1.2. Сплошные линии — траектории до начала расширения облака, пунктир — после начала расширения облака. Огибающие показаны точечной линией, начала траекторий показаны стрелочками.

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

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

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

Примеры численных расчетов распределения концентрации в струе в различные моменты времени показаны на рис. 22–23.

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

Рис. 22. Вихревые домены (серые точки) и пробные частицы дисперсной фазы (цветные точки, цвет соответствует местной концентрации частиц) в плоской импульсной струе; Re = 100, β = 10.

Параметр инерционности стоксовых частиц также посчитан для этих параметров. Безразмерная длительность впрыска Δt = 0.2. Для более инерционных частиц картина распределения концентрации примеси изменяется (рис. 23).

Рис. 23. То же, что на рис. 22. Re = 100, β = 1.

Инерционные частицы сначала собираются в головной части вихревой зоны, а с течением времени — в крестообразной области позади вихрей. Интересно поведение во времени лагранжевой поверхности, состоящей из одних и тех же инерционных частиц (рис. 24).

Рис. 24. Типичное поведение лагранжевой поверхности, состоящей из одних и тех же частиц, в плоской импульсной струе.

Как видно, лагранжева поверхность претерпевает многочисленные самопересечения, что свидетельствует о формировании множественных складок в среде частиц. Этот пример, как и картина траекторий частиц в вихре Ламба–Озеена на рис. 21, иллюстрируют невозможность использования стандартных эйлеровых моделей для описания подобных задач. В то же время предложенный комбинированнный лагранжев подход успешно справляется с расчетом распределения концентрации в таких течениях.

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

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

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

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

В работе [102] численно моделировалось поведение инерционной примеси в задаче о столкновении двух вязких вихревых колец при умеренных числах Рейнольдса, посчитанных по начальной циркуляции скорости в вихревом кольце. Общая постановка задачи показана на рис. 25.

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

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

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

Рис. 26. Вихревые элементы (серые точки) и распределение безразмерной концентрации инерционных частиц (цветные точки) в задаче о столкновении двух вихревых колец. Re = 100, β = 1; безразмерное время t = 1 (a) и t = 2 (б).

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

11. РАЗВИТИЕ КОМБИНИРОВАННОГО ЛАГРАНЖЕВА ПОДХОДА НА СЛУЧАЙ НЕИЗОТЕРМИЧЕСКИХ ТЕЧЕНИЙ С ФАЗОВЫМИ ПЕРЕХОДАМИ

Метод, аналогичный методу вязких вихревых доменов, может быть использован для расчета поля температуры несжимаемой несущей фазы в неизотермических газодисперсных течениях. В работах [103–104] такой метод назван методом тепловых доменов. Уравнение конвективной теплопроводности несущей фазы переписывается в безразмерном дивергентном виде:

$\frac{\partial T}{\partial t}+\text{div}\left[T\left(\mathbf{V}+{\mathbf{V}}_{TD}\right)\right]=\begin{array}{cc}\mathrm{0,}& {\mathbf{V}}_{TD}\end{array}=-\frac{\gamma }{\mathrm{PrRe}}\frac{\nabla T}{T}.$

Здесь V_TD — диффузионная скорость тепловых доменов, γ — показатель адиабаты. Далее исходная область течения разбивается на N малых лагранжевых элементов (тепловых доменов), центры которых движутся со скоростью V_i + V_TDi.

Для вычисления температуры и градиента температуры по эйлеровым координатам, входящего в определения диффузионной скорости, применяется процедура метода сглаженных частиц, аналогичная процедуре, использованной в методе вязких вихревых доменов (9.8)–(9.9). Вместо циркуляций скорости Γ_i в данном случае в формулах будут фигурировать сохраняющиеся в процессе движения интегралы от температуры по площади тепловых доменов Θ_i.

Таким образом, задача определения поля температур несущей фазы сводится к системе обыкновенных дифференциальных уравнений движения центров тепловых доменов. Эти уравнения должны решаться совместно с системами уравнений комбинированного лагранжева подхода, описанного в предыдущем разделе. С помощью такого подхода в работах [103–104] исследовались поля концентрации и изменение размеров испаряющихся капель в неизотермической импульсной двухфазной струе с холодными каплями, впрыскиваемой в область горячего пара.

12. ПРИМЕНЕНИЕ ПОЛНОГО ЛАГРАНЖЕВА ПОДХОДА В СЛУЧАЕ ПУЛЬСИРУЮЩИХ И ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ

В силу большого круга разнообразных прикладных задач (от экологии и физики облаков до ядерных технологий и медицины) проблеме неоднородного распределения дисперсной примеси в турбулентных потоках посвящено много десятков статей (см., например, обзоры [105–108]). Особый интерес в минувшие годы вызывало выяснение физических механизмов формирования локальных зон предпочтительной аккумуляции и высоких градиентов концентрации частиц в турбулентных потоках. Далее упомянем лишь несколько наиболее важных публикаций, посвященных указанной тематике, а основное внимание уделим попыткам применения полного лагранжева подхода для моделирования поведения концентрации частиц в пульсирующих и турбулентных газодисперсных течениях.

При описании распределения дисперсной примеси на больших пространственных масштабах (перенос загрязнений в атмосфере и океане, течения промышленных аэрозолей в системах вентиляции и теплоэнергетических установках) широко применяют инженерные модели на базе уравнения конвективной диффузии с полуэмпирическими коэффициентами турбулентной диффузии частиц [109]. Неоднородность осредненных характеристик турбулентности при этом учитывается за счет переменности коэффициента турбулентной диффузии, а также так называемого эффекта турбофореза — осредненной миграции частиц из более турбулентных зон в менее турбулентные [110].

При детальном анализе структуры поля концентрации дисперсной примеси на мезомасштабах (масштабах когерентных структур в неоднородной турбулентости, колмогоровском масштабе в изотропной турбулентности) оказывается, что мгновенные картины распределения примеси крайне неоднородны даже при очень малой инерционности частиц. Этот факт подтверждается многочисленными расчетами и экспериментами [111].

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

К аналогичному выводу приходят и авторы исследования [113] на основе прямого численного моделирования динамики миллиона стоксовых частиц различной инерционности в трехмерном поле изотропной турбулентности.

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

В работе [115] на основании численных расчетов траекторий многих частиц в случайном поле скорости, моделирующем изотропную турбулентность, показано, что с течением времени траектории малоинерционных частиц, первоначально равномерно распределенных в объеме, стремятся в фазовом пространстве к аттрактору, проекция которого на физическое пространство имеет фрактальную размерность. Это означает, что малоинерционные частицы в однородной турбулентности с течением времени собираются в областях фрактальной структуры. Этот эффект имеет место для частиц, инерционность которых ниже некоторого порогового значения, за счет диссипативного характера их движения в вязкой среде. Более инерционные частицы при перемешивании продолжают оккупировать все физическое пространство.

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

В некоторых работах (например, [116]) отмечается возможность образования каустик на огибающих мгновенных линий тока инерционных частиц в однородной турбулентности. Типичная качественная картина распределения дисперсных частиц в поле однородной турбулентности показана на рис. 27.

Рис. 27. Схематичная картина мгновенного распределения концентрации малоинерционных частиц и каустик [116] в поле однородной турбулентности. Каустики показаны синим цветом.

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

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

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

В связи с работами [117–120] в диссертации [121] с помощью полного лагранжева подхода была исследована модельная задача о нестационарной (пульсирующей) критической точке скорости в газопылевой среде. Точки нулевого ускорения должны двигаться с постоянными скоростями, поэтому можно перейти в инерциальную систему отсчета, связанную с такой точкой. В этой системе отсчета скорость несущей фазы в начале координат равняется нулю, а локальное поле скоростей несущей фазы в простейшем случае будет линейно зависеть от координат с коэффициентом (скоростью растекания), пульсирующим во времени. Взяв одну гармонику во временном законе изменения скорости растекания, локальное двумерное поле безразмерной скорости несущей фазы можно записать в виде

$x=x\text{\hspace{0.17em}sin}(\omega t\begin{array}{cc}),& y\end{array}=-y\text{\hspace{0.17em}sin}(\omega t).$

В работе [121] рассматривалось поведение конечного облака дисперсной фазы радиуса L в окрестности такой точки, масштабы скорости и времени при этом были AL и A^–1, где А — максимальная размерная скорость растекания в критической точке. На основании параметрических расчетов поведения облака частиц с помощью полного лагранжева метода было показано, что в зависимости от значений параметра инерционности частиц (числа Стокса) и безразмерной частоты колебания ω с течением времени облако частиц либо стягиваются в точку в начале координат, либо неограниченно расширяется. При этом существует критическое значение инерционности частиц Stk_с, такое что менее инерционные частицы притягиваются к особой точке при любых частотах колебаний. В обоих случаях происходят гармонические колебания облака частиц с периодическим появлением неограниченных значений концентрации частиц (рис. 28), при этом внутри облака концентрация частиц остается однородной и зависит только от времени. При переходе через ноль якобиана det(J) облако частиц как бы выворачивается наизнанку.

Рис. 28. Примеры поведения траекторий (а) и концентрации (б) частиц в окрестности пульсирующей особой точки в случае, когда точка является притягивающей (сплошные линии 1) или отталкивающей (пунктир 2).

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

Полный лагранжев метод был впервые применен для исследования аккумуляции частиц в турбулентном потоке в работе [122], где авторы рассматривали развитое турбулентное течение газовзвеси в плоском канале при малой концентрации частиц.

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

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

Далее в ряде работ [55, 123, 124] с помощью полного лагранжева подхода исследовалась кластеризация частиц в модельных двумерных случайных полях скорости, заданных отрезками рядов Фурье (всего порядка 200 гармоник) со случайными коэффициентами. Распределения коэффициентов выбирали из условия удовлетворения законам подобия однородных турбулентных течений. Такие поля являются удобным инструментом для моделирования однородной изотропной турбулентности.

В работах [125, 126] также использовался полный лагранжев подход для дисперсной фазы, но поле однородной турбулентности несущей фазы находили из прямого решения трехмерных уравнений Навье–Стокса псевдоспектральным методом. Распределением частиц по скоростям в малом локальном объеме в этих работах также пренебрегалось.

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

Авторы установили также, что существует критическое значение числа Стокса частиц (Stk_cr ~ 0.7), выше которого происходит смена знака производной осредненного по времени значения якобиана перехода от эйлеровых переменных к лагранжевым (как характеристики сжимаемости среды частиц).

При Stk < Stk_cr осредненный по времени локальный объем среды частиц непрерывно уменьшается, т.е. со временем происходит кластеризация частиц. При Stk > Stk_cr осредненный по времени локальный объем дисперсной фазы непрерывно расширяется. При этом по-прежнему могут формироваться как мгновенные зоны кластеризации, так и зоны разрежения частиц.

Рис. 29. Области притяжения (внутри кривой) и отталкивания (вне кривой) рассматриваемой особой точки в пространстве безразмерных определяющих параметров, точка на кривой соответствует критическому числу Stk_c.

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

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

1. в осредненном (или фильтрованном) поле скорости несущей фазы, полученном из решения осредненных уравнений (уравнений Рейнольдса или LES) осредненная динамика среды, состоящей из достаточно инерционных частиц, может быть описана уравнениями “холодного” континуума;
2. при осредненном описании переноса массы среды частиц в лагранжевой форме учитывается дополнительный поток массы за счет мелкомасштабных (подсеточных) флуктуаций скорости частиц, вызванных мелкомасштабными флуктуациями скорости несущей фазы;
3. этот дополнительный поток массы описывается законом Фика с использованием одной из известных моделей для коэффициента турбулентной диффузии частиц;
4. при расчете числовой плотности частиц вдоль выбранной траектории пространственные производные плотности числа частиц (в диффузионном слагаемом) находятся путем проецирования решения, получаемого с помощью полного лагранжева подхода, на эйлерову сетку, что приводит к гибридному лагранжево-эйлерову алгоритму.

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

Итак, в работе [44] предполагается, что осредненная по пространству реализаций (или каким-то другим способом) скорость среды частиц в турбулентном потоке есть <V>_s(r, t), соответствующие этой скорости осредненные траектории частиц носят детерминированный характер и являются достаточно гладкими (рис. 30).

Рис. 30. Осредненная траектория среды частиц (жирная линия) и реальные траектории частиц (ломаные линии) в турбулентном потоке.

Выбрав замкнутый лагранжев объем дисперсной фазы Ω_s(t) c границей Σ_s(t), перемещающийся с осредненной скоростью <V>_s(r, t), и учитывая предположения 1)–3), можно записать закон сохранения массы среды частиц в интегральной и в дифференциальной лагранжевой форме следующим образом:

$$\begin{gathered} \frac{d}{dt}\underset{{\Omega }_s}{\int }{n}_{s}d{\Omega }_s=-\underset{{\Sigma }_s}{\int }{\mathbf{j}}_{sD}\mathbf{n}d{\Sigma }_s=\underset{{\Sigma }_s}{\int }(D_{sT}\nabla n_s)\mathbf{n}d{\Sigma }_s=\underset{{\Omega }_s}{\int }\text{div}(D_{sT}\nabla n_s)d{\Omega }_s, \\ \frac{d}{dt}\left[n_s\right({\mathbf{r}}_0,t\left)\left|\det J\right|\right]=\text{div}(D_{sT}\nabla n_s)\left|\det J\right|. \end{gathered}$$ (12.1)

Здесь j_sD — вектор диффузионного потока массы на границе лагранжева объема Ω_s(t) за счет подсеточных флуктуаций скорости среды частиц, D_s — коэффициент турбулентной диффузии частиц, |det(J)| — модуль якобиана перехода от эйлеровых переменных к лагранжевым. Якобиан, как и ранее в выражении (2.1), взят по модулю для автоматического учета смены ориентации объема при пересечении осредненных траекторий среды частиц.

При стоксовском законе сопротивлении частиц для осредненных траекторий и скорости среды частиц имеем уравнения, аналогичные (2.1):

$\frac{\partial <\mathbf{r}>({\mathbf{r}}_0,t)}{\partial t}=\begin{array}{cc}<\mathbf{V}{>}_s,& \frac{\partial <\mathbf{V}{>}_s({\mathbf{r}}_0,t)}{\partial t}\end{array}=\frac{1}{\tau }(<\mathbf{V}>-<\mathbf{V}{>}_s).$ (12.2)

Чтобы осредненное уравнение импульса частиц (12.2) было справедливо, значения корреляций флуктуационных скоростей частиц $<{v^{\text{'}}}_{si}{v^{\text{'}}}_{sj}>$ должны быть пренебрежимо малы по сравнению с величинами <V_si> <V_sj>. Если поле осредненной (фильтрованной) скорости несущей фазы <V> (r, t) известно, то можно применить полный лагранжев подход и вычислить вдоль выбранной осредненной траектории частиц скорость и компоненты якобиана перехода от эйлеровых к лагранжевым переменным из уравнений, аналогичных (2.2).

Отличие от ранее рассмотренных случаев состоит в том, что для нахождения текущего значения концентрации частиц на выбранной траектории теперь следует решать дифференциальное уравнение (второе уравнение в выражении (12.1)). Чтобы замкнуть такую постановку задачи, необходимо выбрать модель для коэффициента турбулентной диффузии частиц D_sT и алгоритм вычисления на выбранной траектории частиц производных по пространственным координатам от n_s, входящих в правую часть второго уравнения (12.1). Например, если фильтрованное поле скоростей несущей фазы < V > (r, t) находится с помощью метода LES, то истинное поле скорости можно представить в виде

$\mathbf{V}(\mathbf{r},t)=<\mathbf{V}>(\mathbf{r},t)+{\mathbf{v}}^{\mathbf{\text{'}}}(\mathbf{r},t),$

где v′ — скорость флуктуаций, определяемая масштабом пространственного фильтра Δ. В работе [44] для вычисления коэффициента турбулентной диффузии частиц была принята часто используемая модель

$D_{sT}=\Lambda {\left(\frac{2}{3}K\right)}^{1/2}\text{}\begin{array}{cc}\Delta ,& K\end{array}=\frac{<{\mathbf{v}}^{\mathbf{\text{'}}}{\mathbf{v}}^{\mathbf{\text{'}}}>}{2}\begin{array}{cc},& \Lambda \end{array}=\frac{U_s}{U}.$

Здесь Λ — функция числа Стокса частицы, равная отношению масштабов флуктуационных скоростей фаз U_s и U; K — кинетическая энергия пульсаций скорости несущей фазы.

В случае очень инерционных частиц коэффициент турбулентной диффузии обнуляется, и модель среды частиц превращается в модель “холодного” континуума, движущегося на фоне осредненного поля скорости несущей фазы. В противоположном случае малоинерционных частиц модель переходит в уравнение конвективной диффузии пассивной примеси. В работе [44] при вычисления диффузионного члена в уравнении для концентрации частиц (12.1) значения концентрации на лагранжевой сетке пересчитывались на эйлерову сетку, используемую для вычисления параметров несущей фазы. Таким образом, алгоритм полного лагранжева подхода дополнялся процедурой лагранжево-эйлерова пересчета параметров дисперсной фазы.

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

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

$U(x,t)=\text{sin}\left(\frac{\pi t}{10}\right)\text{sin}\left(\pi x\right)+\frac{1}{16}\text{sin}\left(\frac{\pi t}{160}\right)\text{sin}\left(\frac{\pi x}{16}\right).$

В качестве характерных масштабов при обезразмеривании взяты максимальная скорость основного колебания и линейный масштаб L, равный четверти начальной высоты слоя частиц.

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

Вычислялись также отношение продольной компоненты тензора пульсационных напряжений в среде частиц к потоку импульса, переносимого со среднемассовой скоростью, и отношение диффузионного потока частиц к среднему потоку массы. Показано, что для достаточно инерционных частиц (число Stk = 5) вкладом пульсационных скоростей частиц в поток импульса вполне можно пренебречь, но при этом следует учитывать диффузионный вклад в перенос массы частиц. Типичный расчет траекторий и концентрации частиц на фоне “фильтрованного” поля скоростей несущей фазы показан на рис. 31.

Рис. 31. Расчет траекторий и концентрации частиц в “фильтрованном” пульсирующем поле скоростей в вертикальной трубе; ось абсцисс — безразмерное время, ось ординат — безразмерная координата, цветом показана величина безразмерной концентрации частиц.

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

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

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

Рис. 32. Рассчитанные мгновенные поля безразмерной концентрации частиц (цвет соответствует величине log(n_s)) различной инерционности в различные моменты безразмерного времени: (а) Stk = 0.01, t = 1; (б) Stk = 0.1, t = 1; (в) Stk = 1, t = 1; (г) Stk = 1, t = 23.

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

В работе [44] был также предложен так называемый полный лагранжев метод второго порядка, основанный на вычислении, наряду с компонентами якобиана преобразования от эйлеровых переменных к лагранжевым, компонент гессиана соответствующего преобразования. Такой подход позволяет вычислять среднее значение числа частиц в малом, но конечном лагранжевом объеме вблизи каустик, где якобиан стремится к нулю. В недавно опубликованной работе [127] этот метод был применен к расчету концентрации микрокапель в трехмерных вихревых образованиях, формирующихся при кашле больного COVID-19.

В завершение обзора отметим недавнюю публикацию [128], в которой авторы предложили эффективный алгоритм пересчета концентрации дисперсной примеси, рассчитанной на выбранных пересекающихся лагранжевых траекториях частиц, на фиксированную эйлерову сетку. Этот алгоритм объединяет идеи полного лагранжева подхода и метода SPH: предлагается выбирать величину носителя “шапочки”, аппроксимирующей дельта-функцию в методе SPH, пропорциональной значению якобиана перехода от эйлеровых переменных к лагранжевым. Такой подход позволил авторам резко сократить число рассчитываемых траекторий частиц при сохранении требуемой точности вычисления концентрации частиц.

ЗАКЛЮЧЕНИЕ

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

Главные достоинства метода — его простота и экономичность, обусловленные возможностью построения полей концентрации частиц в сложных неодномерных и нестационарных потоках на основе расчета небольшого числа траекторий. Во многих публикациях подчеркивается значительная (на порядки) экономия времени вычисления концентрации частиц при использовании полного лагранжева подхода. Например, в работе [53] отмечено, что для повышения относительной точности расчета концентрации дисперсной фазы в зонах пересекающихся траекторий частиц на 0.1% с помощью стандартных лагранжевых подходов требуется расчет порядка 10⁴ траекторий частиц, приходящихся на одну ячейку эйлеровой сетки, на которой рассчитаны параметры несущей фазы. В то же время при использовании полного лагранжева подхода достаточно расчета одной траектории частиц из каждого слоя “складки”.

Дальнейшее развитие полного лагранжева подхода, по-видимому, будет связано с учетом фазовых переходов, полидисперсности и влияния частиц на параметры несущей фазы, моделированием коллективного поведения активных (“самодвижущихся”) частиц, анализом механизмов кластеризации частиц в трехмерных вихревых потоках и когерентных структурах турбулентных двухфазных потоков, а также объединением идей полного лагранжева подхода с идеями метода SPH.

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

БЛАГОДАРНОСТИ

Автор выражает глубокую благодарность кандидатам физико-математических наук И.В. Голубкиной, Н.А. Лебедевой и О.Д. Рыбдыловой, а также доктору А. Папоуцакису (A. Papoutsakis) (Великобритания) и профессору Ван Бо-И (Wang Bo-Yi) (КНР), принимавшим активное участие в совместных работах по развитию полного лагранжева подхода. Автор благодарен также профессорам Ю.М. Циркунову и В.В. Измоденову, любезно предоставившим иллюстрации из своих работ.

ФИНАНСИРОВАНИЕ

Работа выполнена по открытому плану МГУ им. М.В. Ломоносова.

КОНФЛИКТ ИНТЕРЕСОВ

Автор заявляет, что у него нет конфликта интересов.

Об авторах

А. Н. Осипцов

Автор, ответственный за переписку.
Email: osiptsov@imec.msu.ru

Научно-исследовательский институт механики

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

Дополнительные файлы

Доп. файлы

Действие

1. JATS XML

Скачать

2. Рис. 1. Пример расчета ∂ρ/∂x при обтекании сферы вертикально стратифицированным потоком в вертикальной срединной плоскости лагранжевым методом (а) и теневая фотография из работы [47] (б).

Скачать (315KB)

Метаданные

3. Рис. 2. Мгновенная картина модуля градиента плотности при обтекании цилиндра горизонтально стратифицированным потоком с гармоническими колебаниями плотности.

Скачать (394KB)

Метаданные

4. Рис. 3. Траектории частиц (пунктир) и изолинии безразмерной концентрации частиц (сплошные линии с цифрами) в потенциальном течении вблизи пробоотборника аэрозоля; стенки пробоотборника — жирные линии; K = U∞ /Ua = 0.2, Stk = 1.

Скачать (133KB)

Метаданные

5. Рис. 4. Течение в окрестности “косой” критической точки [62]. Синие линии — линии тока несущей среды, красные — траектории частиц, пунктир — огибающие траектории частиц (каустики).

Скачать (396KB)

Метаданные

6. Рис. 5. Распределения вертикальных скоростей несущей фазы (синие линии), частиц (красные линии) (а) и безразмерной концентрации частиц (б) в окрестности критической точки [62].

Скачать (167KB)

Метаданные

7. Рис. 6. Картины траекторий частиц (левые рисунки) и изолиний безразмерной концентрации частиц (правые рисунки) при обтекании цилиндра невязким сверхзвуковым потоком [64]. Верхние рисунки соответствуют режиму отсутствия инерционного осаждения частиц, а нижние — режиму с отражением частиц.

Скачать (511KB)

Метаданные

8. Рис. 7. Картины траекторий частиц (верхние рисунки) и соответствующие профили суммарной концентрации частиц (нижние рисунки), формирующиеся вдали за точкой взаимодействия волн; ytr на последнем рисунке отсчитывается по нормали к контактной поверхности.

Скачать (564KB)

Метаданные

9. Рис. 8. Теневая фотография из работы [67] и расчеты распределения числа Маха в ударном слое из работы [66] для случая, когда на боковую поверхность цилиндра приходит контактный разрыв.

Скачать (225KB)

Метаданные

10. Рис. 9. Типичные картины траекторий частиц различной инерционности и поля безразмерной концентрации дисперсной фазы, соответствующие расчетам из работы [66] (отраженные частицы не учитывались).

Скачать (395KB)

Метаданные

11. Рис. 10. Траектории частиц (линии 2) и каустики (жирные линии 3) в системе отсчета, связанной с движущейся ударной волной вдоль стенки [39] (граница пограничного слоя — линия 1); (а) — частицы поступают в пограничный слой из внешнего потока; (б) — частицы поступают с поверхности разрушающегося слоя дисперсной фазы, лежащего на стенке.

Скачать (219KB)

Метаданные

12. Рис. 11. Типичные траектории частиц различной инерционности в системе отсчета, связанной с ударной волной, движущейся в узком канале (поперечная координата отнесена к радиусу канала, продольная — к длине скоростной релаксации частиц).

Скачать (212KB)

Метаданные

13. Рис. 12. Типичная картина траекторий частиц в цилиндрических координатах в осесимметричном сужающемся микроканале. Каустики и области накопления частиц формируются на огибающих траектории частиц.

Скачать (157KB)

Метаданные

14. Рис. 13. Типичные восходящие (а) и нисходящие (б) траектории “тяжелых” частиц в течении типа торнадо [84].

Скачать (168KB)

Метаданные

15. Рис. 14. Картина трубок тока несущей фазы (синие пунктирные линии), трубок тока среды “тяжелых” частиц (зеленые сплошные линии), каустик (зеленые пунктирные линии) и изолиний величины, обратной концентрации дисперсной фазы (разноцветные сплошные линии с числами), в плоскости r, z [84].

Скачать (167KB)

Метаданные

16. Рис. 15. Формирование спиралевидной верхней части “чашеобразной” зоны накопления частиц в течении типа торнадо [84]. Зеленые сплошные линии — трубки тока частиц, цветные пунктирные линии с числами — изолинии величины, обратной концентрации дисперсной фазы.

Скачать (175KB)

Метаданные

17. Рис. 16. Траектории (а) и каустики (б) космической пыли в срединной плоскости гелиосферы, перпендикулярной токовому слою магнитного поля. Солнце находится в начале координат, координаты измеряются в астрономических единицах.

Скачать (302KB)

Метаданные

18. Рис. 17. Траектории частиц одинаковой инерционности в экваториальной плоскости при различных вкладах центробежной силы; (а) — меньший вклад центробежной силы, (б) — больший вклад центробежной силы.

Скачать (95KB)

Метаданные

19. Рис. 18. Типичная пространственная траектория инерционной частицы (сплошная линия) (а) и сечения трубок тока частиц вертикальной плоскостью (б) из [86].

Скачать (232KB)

Метаданные

20. Рис. 19. Вихревые домены (серые точки) и пробные частицы дисперсной фазы (цветные точки, цвет соответствует локальной концентрации частиц) в вихре Ламба–Озеена для “тяжелых” стоксовских частиц в моменты t = 2, 3, 4 при Re = 100, Fr = 1, η = 0.

Скачать (356KB)

Метаданные

21. Рис. 20. Распределения концентрации “легких” частиц в начальные моменты времени при Re = 100, Fr = ∞, η = 1.2. Сплошные линии соответствуют использованию аналитических формул для параметров несущей фазы, крестики — расчету с помощью метода вихревых доменов.

Скачать (133KB)

Метаданные

22. Рис. 21. Типичные траектории “легких” частиц в вихре Ламба–Озеена при Re = 100, Fr = ∞, η = 1.2. Сплошные линии — траектории до начала расширения облака, пунктир — после начала расширения облака. Огибающие показаны точечной линией, начала траекторий показаны стрелочками.

Скачать (184KB)

Метаданные

23. Рис. 22. Вихревые домены (серые точки) и пробные частицы дисперсной фазы (цветные точки, цвет соответствует местной концентрации частиц) в плоской импульсной струе; Re = 100, β = 10.

Скачать (284KB)

Метаданные

24. Рис. 23. То же, что на рис. 22. Re = 100, β = 1.

Скачать (300KB)

Метаданные

25. Рис. 24. Типичное поведение лагранжевой поверхности, состоящей из одних и тех же частиц, в плоской импульсной струе.

Скачать (226KB)

Метаданные

26. Рис. 25. Постановка задачи о перемешивании дисперсной примеси двумя сталкивающимися вихревыми кольцами в вязкой среде.

Скачать (112KB)

Метаданные

27. Рис. 26. Вихревые элементы (серые точки) и распределение безразмерной концентрации инерционных частиц (цветные точки) в задаче о столкновении двух вихревых колец. Re = 100, β = 1; безразмерное время t = 1 (a) и t = 2 (б).

Скачать (426KB)

Метаданные

28. Рис. 27. Схематичная картина мгновенного распределения концентрации малоинерционных частиц и каустик [116] в поле однородной турбулентности. Каустики показаны синим цветом.

Скачать (387KB)

Метаданные

29. Рис. 28. Примеры поведения траекторий (а) и концентрации (б) частиц в окрестности пульсирующей особой точки в случае, когда точка является притягивающей (сплошные линии 1) или отталкивающей (пунктир 2).

Скачать (130KB)

Метаданные

30. Рис. 29. Области притяжения (внутри кривой) и отталкивания (вне кривой) рассматриваемой особой точки в пространстве безразмерных определяющих параметров, точка на кривой соответствует критическому числу Stkc.

Скачать (63KB)

Метаданные

31. Рис. 30. Осредненная траектория среды частиц (жирная линия) и реальные траектории частиц (ломаные линии) в турбулентном потоке.

Скачать (49KB)

Метаданные

32. Рис. 31. Расчет траекторий и концентрации частиц в “фильтрованном” пульсирующем поле скоростей в вертикальной трубе; ось абсцисс — безразмерное время, ось ординат — безразмерная координата, цветом показана величина безразмерной концентрации частиц.

Скачать (161KB)

Метаданные

33. Рис. 32. Рассчитанные мгновенные поля безразмерной концентрации частиц (цвет соответствует величине log(ns)) различной инерционности в различные моменты безразмерного времени: (а) Stk = 0.01, t = 1; (б) Stk = 0.1, t = 1; (в) Stk = 1, t = 1; (г) Stk = 1, t = 23.

Скачать (1MB)

Метаданные