Журнал вычислительной математики и математической физики, 2019, T. 59, № 11, стр. 1899-1914

ВВЕДЕНИЕ

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

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

На основе модельного кинетического уравнения в [6], [7] была построена регуляризованная (квазигазодинамическая, КГД) система уравнений для описания течений бинарной смеси газов (подробнее о КГД моделях однокомпонентных газов см. [7], [8]). Предполагалось, что каждая компонента газовой смеси имеет свои плотность, скорость и температуру. В отсутствие химических реакций (приближение холодного газа) обменные слагаемые в уравнениях баланса масс компонент отсутствовали, а в уравнениях баланса импульса и полной энергии каждой компоненты присутствовали и описывали взаимодействие между компонентами смеси. При этом за основу брался вид обменных слагаемых для кинетических уравнений и предполагалось, что газы являются одноатомными. Тем самым такие КГД уравнения для смесей были ориентированы на описание простейших одноатомных газов. Эта система уравнений была протестирована в расчетах одномерных течений разреженного газа на примере задач о неподвижной ударной волне и течении газов между двумя резервуарами. Результаты сравнивались с полученными методом Монте-Карло (DSMC) [9] и были достаточно точными. Дополнительные приложения были даны в [10].

На основе этих уравнений в [7] было также построено гомогенное (одножидкостное) приближение в предположении равенства скоростей и температур отдельных компонент смеси. При этом предполагалось равенство показателей адиабат $\gamma $ и чисел Прандтля ${{\alpha }_{{{\text{Pr}}}}}$ компонент и понятие внутренней энергии смеси явным образом не вводилось. Работоспособность такой модели была также проверена на примере тех же двух задач.

В дальнейшем в [11] обменные слагаемые в КГД системе были обобщены на случай течений многоатомных газов, с произвольными $\gamma $ и ${{\alpha }_{{{\text{Pr}}}}}$. Кроме того, сами КГД уравнения для смеси газов были записаны в другой, более стандартной для газовой динамики форме, особенно удобной для их последующей дискретизации. Там же было показано, что КГД система для смеси газов удовлетворяет законам сохранения суммарной полной энергии и неубывания суммарной энтропии. Физический смысл предложенного обобщения обменных членов был выяснен в [12].

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

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

Однако попытки численного моделирования на основе такой модели гомогенных бинарных смесей неразреженных газов в задаче типа Рэлея–Тейлора показали наличие слишком большой численной диффузии в уравнениях баланса массы компонент из-за больших скачков давления компонент. Поэтому для подобных смесей газов было решено существенно модифицировать регуляризующие слагаемые всех уравнений движения и тем самым, во-вторых, построить еще одну регуляризованную систему уравнений гомогенной бинарной смеси, для неразреженных газов. Отметим, что эта модификация имеет некоторые схожие элементы с существенно более простой, чем в данной статье, квазигидродинамической регуляризацией в задачах о течении бинарной смеси (иными словами, двухкомпонентной жидкости) при малых изменениях суммарной плотности, но с дополнительным учетом межфазных эффектов (описываемых уравнениями Навье–Стокса–Кана–Хилларда для сжимаемой жидкости) [13], [14]. Вторая построенная в данной статье система уравнений в основном сохраняет отмеченные выше свойства первой, но содержит только суммарные регуляризующие диффузионные потоки. Подчеркнем, что снова обеспечивается справедливость уравнения баланса полной энтропии смеси с неотрицательным производством энтропии (теорема 2), при этом само уравнение существенно отличается от соответствующего ему для первой системы уравнений. Выполнение указанного свойства имеет принципиальное значение, поскольку обеспечивает физическую непротиворечивость модели, и добиться этого для подобной регуляризации для смесей оказалось непросто.

Для анализа свойств математических и численных моделей смесей в качестве тестовой часто используется задача моделирования гравитационной неустойчивости Рэлея–Тейлора, возникающей в поле силы тяжести на границе раздела между расположенным сверху тяжелым газом и расположенным снизу легким газом. Другим типом тестовой задачи является задача о неустойчивости Рихтмайера–Мешкова, возникающей при прохождении ударной волны через границу раздела двух газов. Имеется большая недавняя обзорная статья [15] по обоим типам неустойчивости. Численное исследование и детальный анализ результатов по неустойчивости Рэлея–Тейлора для сжимаемых газов приведены в том числе в [16]. В данной работе мы сосредоточимся на задаче о гравитационной неустойчивости, включая примеры типа имеющихся в [16]. При этом существенным моментом является численное моделирование различных бинарных смесей, во всем диапазоне $(0,1)$ числа Атвуда, характеризующего разброс молекулярных весов компонент смеси.

1. РЕГУЛЯРИЗОВАННЫЕ УРАВНЕНИЯ ДВИЖЕНИЯ ОБЩЕЙ И ГОМОГЕННОЙ БИНАРНЫХ СМЕСЕЙ УМЕРЕННО-РАЗРЕЖЕННЫХ ГАЗОВ

1.1. Регуляризованная КГД система уравнений бинарной смеси умеренно-разреженных газов в форме [11] состоит из следующих уравнений баланса массы, импульса и полной энергии для газов $\alpha = a,b$

Здесь и ниже ${{\partial }_{t}}$ и ${{\partial }_{i}}$ – частные производные по $t$ и ${{x}_{i}}$, операторы $div$ и $\nabla $ берутся по пространственным координатам $x = ({{x}_{1}},\; \ldots ,\;{{x}_{n}})$, где $n = 1,2,3$, причем дивергенция тензора берется по его первому индексу, а $ \otimes $ и $ \cdot $ – знаки тензорного и скалярного произведения векторов.

Основными искомыми функциями являются ${{\rho }_{\alpha }} > 0$, ${{{\mathbf{u}}}_{\alpha }} = ({{u}_{{1\alpha }}},\; \ldots ,\;{{u}_{{n\alpha }}})$, ${{\theta }_{\alpha }} > 0$ – плотность, скорость, абсолютная температура газа $\alpha $, зависящие от $(x,t)$. Участвуют также полная энергия, давление и удельная внутренняя энергия

газа $\alpha $, который считается совершенным политропным. Здесь присутствуют постоянные ${{R}_{\alpha }} = {{R}_{0}}{\text{/}}{{M}_{\alpha }} > 0$, причем ${{R}_{0}}$ – универсальная газовая постоянная, а ${{M}_{\alpha }}$ – молекулярный вес, и ${{c}_{{V\alpha }}}$ — удельная теплоемкость при постоянном объеме. Используется также другая форма записи давления: где ${{\gamma }_{\alpha }}$ – показатель адиабаты, а ${{c}_{{p\alpha }}}$ – удельная теплоемкость при постоянном давлении.

В этих уравнениях ${{\Pi }_{\alpha }} = \Pi _{\alpha }^{{NS}} + \Pi _{\alpha }^{\tau }$ – тензор вязких напряжений. При этом $\Pi _{\alpha }^{{NS}}$ – это классический тензор вязких напряжений Навье–Стокса

с коэффициентами динамической ${{\mu }_{\alpha }} = {{\mu }_{\alpha }}({{\rho }_{\alpha }},{{\theta }_{\alpha }}) > 0$ и объемной вязкости ${{\lambda }_{\alpha }} = {{\lambda }_{\alpha }}({{\rho }_{\alpha }},{{\theta }_{\alpha }}) \geqslant 0$ и единичным тензором $\mathbb{I}$ (порядка $n$). Кроме того, есть регуляризующий тензор. Отметим, что в последней формуле ${{\gamma }_{\alpha }}{{p}_{\alpha }} = c_{{s\alpha }}^{2}{{\rho }_{\alpha }}$, где $c_{{s\alpha }}^{2} = {{\gamma }_{\alpha }}({{\gamma }_{\alpha }} - 1){{\varepsilon }_{\alpha }}$ – квадрат скорости звука компоненты $\alpha $. Величины ${{{\mathbf{F}}}_{\alpha }}$ и ${{Q}_{\alpha }} \geqslant 0$ – заданные плотность массовой силы и мощность теплового источника, а $\tau > 0$ – параметр релаксации, который может зависеть от всех искомых функций ${{\rho }_{a}}$, ${{{\mathbf{u}}}_{a}}$, ${{\theta }_{a}}$ и ${{\rho }_{b}}$, ${{{\mathbf{u}}}_{b}}$, ${{\theta }_{b}}$.

Величины ${{{\mathbf{S}}}_{{u,\alpha }}}$, ${{S}_{{E,\alpha }}}$ – обменные слагаемые, зависящие от всех искомых функций и связывающие между собой уравнения для газов $a$ и $b$. Они таковы, что

Конкретный вид ${{{\mathbf{S}}}_{{u,\alpha }}}$, ${{S}_{{E,\alpha }}}$ дан в [11] (см. также физическую мотивацию в [12]); ниже он не потребуется и поэтому здесь опущен.

Регуляризующие скорости ${{{\mathbf{w}}}_{\alpha }}$ и ${{{\mathbf{\hat {w}}}}_{\alpha }}$ таковы

при этом ${{\rho }_{\alpha }}({{{\mathbf{u}}}_{\alpha }} - {{{\mathbf{w}}}_{\alpha }})$ – плотность потока массы газа $\alpha $.

Тепловой поток ${{{\mathbf{q}}}_{\alpha }}$ можно выразить формулами

с коэффициентом теплопроводности ${{\varkappa }_{\alpha }} = {{\varkappa }_{\alpha }}({{\rho }_{\alpha }},{{\theta }_{\alpha }}) > 0$ и регуляризующим тепловым потоком ${\mathbf{q}}_{\alpha }^{\tau }$.

В регуляризованных уравнениях могут использоваться как физические коэффициенты вязкости и теплопроводности, так и искусственные; для последних чаще всего используются формулы (см. [7]):

где ${{\alpha }_{{S\alpha }}} > 0$ и ${{\alpha }_{{{\text{Pr}}\alpha }}} > 0$ – числа Шмидта и Прандтля компоненты $\alpha $.

Приведенная система регуляризованных уравнений бинарной смеси является достаточно сложной и содержит $2(n + 2)$ искомых скалярных функций ${{\rho }_{\alpha }}$, ${{{\mathbf{u}}}_{\alpha }} = ({{u}_{{1\alpha }}},\; \ldots ,\;{{u}_{{n\alpha }}})$, ${{\theta }_{\alpha }}$, $\alpha = a,b$. Поэтому практический интерес представляют и упрощенные модели движения бинарной смеси.

1.2. Для гомогенной бинарной смеси полагается ${{{\mathbf{u}}}_{a}} = {{{\mathbf{u}}}_{b}} = {\mathbf{u}}$, ${{\theta }_{a}} = {{\theta }_{b}} = \theta $ (см. [1]). Пусть также ${{{\mathbf{F}}}_{a}} = {{{\mathbf{F}}}_{b}} = {\mathbf{F}}$. Чтобы вывести уравнения движения такой смеси, выполним агрегирование выписанных выше уравнений по $\alpha = a,b$. Уравнение баланса массы (1.1) просто перепишем в виде

положив величина ${{{\mathbf{\hat {w}}}}^{{(\alpha )}}}$ потребуется чуть ниже.

Введем суммарные плотность $\rho = {{\rho }_{a}} + {{\rho }_{b}}$ и давление $p = {{p}_{a}} + {{p}_{b}}$. Сложив по $\alpha = a,b$ уравнения баланса импульса вида (1.2) с ${{{\mathbf{u}}}_{a}} = {{{\mathbf{u}}}_{b}} = {\mathbf{u}}$, ${{\theta }_{a}} = {{\theta }_{b}} = \theta $, получим агрегированное уравнение баланса импульса гомогенной смеси

Здесь суммарные регуляризующие скорости и слагаемые тензора вязких напряжений $\Pi = {{\Pi }^{{NS}}} + {{\Pi }^{\tau }}$ таковы с $Q: = {{Q}_{a}} + {{Q}_{b}}$, причем учтено важное (хотя и очевидное) свойство – аддитивность регуляризующих импульсов

Введем суммарные удельную внутреннюю энергию и теплоемкость при постоянном объеме

Сложив по $\alpha = a,b$ уравнения баланса полной энергии вида (1.3) с ${{{\mathbf{u}}}_{a}} = {{{\mathbf{u}}}_{b}} = {\mathbf{u}}$, ${{\theta }_{a}} = {{\theta }_{b}} = \theta $, получим агрегированное уравнение баланса полной энергии гомогенной смеси Здесь суммарные полная энергия, тепловой поток и его регуляризующая составляющая выражаются формулами причем учтены аддитивность $\rho $, $\rho \varepsilon $, ${{c}_{V}}\rho $, $\rho {\mathbf{w}}$. Обменные члены в уравнениях (1.13) и (1.19) сократились в силу свойств (1.7).

Выведенная система регуляризованных уравнений гомогенной бинарной смеси гораздо проще исходной и содержит только $n + 3$ искомые скалярные функции ${{\rho }_{a}}$, ${{\rho }_{b}}$, ${\mathbf{u}} = ({{u}_{1}},\; \ldots ,\;{{u}_{n}})$, $\theta $.

Можно дополнительно ввести суммарные теплоемкости при постоянном давлении, газовую постоянную, показатель адиабаты по формулам

тогда $\gamma - 1 = R{\text{/}}{{c}_{V}}$ и поэтому верны другие естественные формулы для давления смеси Во избежание недоразумений подчеркнем, что при внешней похожести в отличие от уравнений состояния однородного совершенного политропного газа в формулах для гомогенной бинарной смеси величины $R$, ${{c}_{p}}$, ${{c}_{V}}$, $\gamma $ не постоянные, а функции концентраций компонент ${{\rho }_{a}}{\text{/}}\rho $ и ${{\rho }_{b}}{\text{/}}\rho = 1 - {{\rho }_{a}}{\text{/}}\rho $.

В (1.16) и (1.19) можно использовать формулы

Здесь $c_{s}^{2}$ играет роль квадрата скорости звука смеси.

Обратим внимание на то, что при выполненном агрегировании одни величины берутся аддитивно, т.е. просто суммируются (как ${{\rho }_{\alpha }}$, ${{p}_{\alpha }}$ и т.д.), а другие берутся как линейные комбинации с весами-концентрациями ${{\rho }_{\alpha }}{\text{/}}\rho $ (как ${{{\mathbf{w}}}^{{(\alpha )}}}$, ${{\varepsilon }_{\alpha }}$ и т.д.).

Для контроля результата отметим, что в простейшем случае газов с одинаковыми ${{c}_{{V\alpha }}}$ и ${{c}_{{p\alpha }}}$ имеем

Как следствие, $\varepsilon = {{\varepsilon }_{\alpha }} = {{c}_{V}}\theta $ и с учетом формул (1.17) уравнение баланса полной энергии и выражение для ${{{\mathbf{q}}}^{\tau }}$ принимают стандартный для однокомпонентного газа вид ср. с (1.3), (1.10) (здесь в $\Pi $ и ${\mathbf{q}}$ стоят суммарные коэффициенты вязкости и теплопроводности).

Далее, при $\tau = 0$ (но без использования формул (1.11)) полученная система уравнений (1.12), (1.13), (1.19) переходит в систему уравнений типа Навье–Стокса для бинарной смеси

с введенными выше $p$, ${{\Pi }^{{NS}}}$, $Q$ и $ - {\mathbf{q}} = ({{\varkappa }_{a}} + {{\varkappa }_{b}})\nabla \theta $. Если же и ${{\mu }_{\alpha }} = {{\lambda }_{\alpha }} = {{\varkappa }_{\alpha }} = 0$, то получается система уравнений типа Эйлера для бинарной смеси

Выведем из полученных регуляризованных уравнений набор следствий стандартного типа. Сложение уравнений (1.12) по $\alpha $ с учетом аддитивности $\rho $ и $\rho {\mathbf{w}}$ дает естественное уравнение баланса суммарной плотности

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

Скалярное умножение уравнения баланса импульса (1.13) на ${\mathbf{u}}$ с использованием формул дифференцирования произведения функций и уравнения (1.29) приводит к уравнению баланса кинетической энергии

Вычитание из уравнения баланса полной энергии (1.19) последнего уравнения с учетом формулы $\rho {\mathbf{w}} = \tau div(\rho {\mathbf{u}}){\mathbf{u}} + \rho {\mathbf{\hat {w}}}$ приводит к уравнению баланса внутренней энергии

где $\nabla {\mathbf{u}} = \left\{ {{{\partial }_{i}}{{u}_{j}}} \right\}_{{i,j = 1}}^{n}$ и $:$ означает скалярное произведение тензоров.

Введем энтропии компоненты $\alpha $ и смеси

где ${{S}_{{0\alpha }}} = {\text{const}}$. Пусть ${{\left| {\mathbb{D}({\mathbf{u}})} \right|}^{2}} = \mathbb{D}({\mathbf{u}}):\mathbb{D}({\mathbf{u}})$.

Теорема 1. Верно следующее регуляризованное уравнение баланса энтропии гомогенной бинарной смеси умеренно-разреженных газов:

с производством энтропии ${{\mathcal{P}}^{{NS}}} + \mathcal{P}_{a}^{\tau } + \mathcal{P}_{b}^{\tau }$, где причем $\mathcal{P}_{\alpha }^{\tau } \geqslant 0$ при условии $\tau ({{\gamma }_{\alpha }} - 1){{Q}_{\alpha }} \leqslant 4{{p}_{\alpha }}$, $\alpha = a,b$.

Доказательство. Ясно, что с учетом аддитивности $\rho s$ и уравнения баланса массы компоненты (1.12) имеем

Используя определение ${{s}_{\alpha }}$ и еще дважды применяя уравнение (1.12), получаем Используя уравнение баланса внутренней энергии смеси (1.31), формулу и разагрегирование некоторых введенных выше суммарных величин, выводим где

Как нетрудно видеть, величину ${{\mathcal{P}}^{{NS}}}$ можно переписать в форме (1.34). Кроме того, величину $\mathcal{P}_{\alpha }^{\tau }$ можно упростить так

Известно, что последнее выражение можно преобразовать к виду (1.35); при ${{Q}_{\alpha }} = 0$ (см., например, [7]) и краткий вывод в [17, утверждение 2]. Случай ${{Q}_{\alpha }} \ne 0$ охватывается так же, как в [18] (см. также ниже доказательство теоремы 2).

Выведенное уравнение баланса энтропии делает построенную модель физически корректной. При этом оно выглядит так, как если бы было получено суммированием уравнений баланса энтропии компонент с ${{{\mathbf{S}}}_{{u,\alpha }}} = 0$ при ${{{\mathbf{u}}}_{a}} = {{{\mathbf{u}}}_{b}} = {\mathbf{u}}$, ${{\theta }_{a}} = {{\theta }_{b}} = \theta $ (ср. с [11, теорема 1]). Разумеется, на самом деле это не так, поскольку для гомогенной бинарной смеси уравнения баланса импульса, полной энергии и энтропии для отдельных компонент уже не выполняются.

Обсудим характер регуляризации в построенных уравнениях. В уравнении баланса массы компоненты (1.12) участвуют неагрегированные регуляризующие импульсы ${{\rho }_{\alpha }}{{{\mathbf{w}}}^{{(\alpha )}}}$, зависящие от плотности и давления только той же компоненты. В уравнение же баланса импульса (1.13), как и в уравнение баланса суммарной массы (1.29), входят, напротив, только суммарные регуляризующие импульсы $\rho {\mathbf{w}}$, $\rho {\mathbf{\hat {w}}}$. При этом в формуле для ${{\Pi }^{\tau }}$ (см. (1.16) и (1.24)) можно считать, что агрегируются также квадраты скоростей звука компонент в сочетании с использованием суммарной плотности. В конвективные слагаемые уравнения баланса полной энергии (1.19) (а также уравнения баланса полной энтропии (1.33)) опять входят неагрегированные регуляризующие импульсы ${{\rho }_{a}}{{{\mathbf{w}}}^{{(a)}}}$, ${{\rho }_{b}}{{{\mathbf{w}}}^{{(b)}}}$ (см. также (1.25)). При этом в формуле (1.21) для ${{{\mathbf{q}}}^{\tau }}$ агрегируются величины ${{c}_{{V\alpha }}}$ (и дают ${{c}_{V}}$) и суммируются величины ${{R}_{\alpha }}\nabla {{\rho }_{\alpha }}$.

Отметим, что описанная выше процедура построения регуляризованных уравнений гомогенной бинарной смеси отнюдь не единственно возможная. Вместо агрегирования регуляризованных уравнений (общей) бинарной смеси можно пытаться строить такие уравнения на основе уравнений типа Навье–Стокса бинарной смеси (1.26)–(1.28) в том числе по аналогии с процедурами из [19] или [20] в случае однокомпонентного газа. Они даже приводят к более простому, чем (1.25), регуляризующему выражению в левой части уравнения баланса полной энергии смеси (1.19), однако не позволяют в итоге обеспечить в уравнении баланса энтропии смеси неотрицательность производства энтропии.

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

где $(\rho ,{\mathbf{u}},\theta )$ и $(\hat {\rho },{\mathbf{\hat {u}}},\hat {\theta })$ – значения искомых функций на текущем и следующем слоях по времени, а $\hat {E} = \tfrac{1}{2}\hat {\rho }{{\left| {{\mathbf{\hat {u}}}} \right|}^{2}} + {{c}_{V}}\hat {\rho }\hat {\theta }$, $\Delta t$ – шаг по времени, ${\mathbf{d}}$ – вектор из параметров газа $R$, ${{c}_{V}}$, ${{\alpha }_{{{\text{Sc}}}}}$, ${{\alpha }_{{{\text{Pr}}}}}$ и функции $Q$. Выражения ${{A}_{1}}[{\mathbf{d}}](\rho ,{\mathbf{u}},\theta )$, ${{A}_{2}}[{\mathbf{d}}](\rho ,{\mathbf{u}},\theta )$, ${{A}_{3}}[{\mathbf{d}}](\rho ,{\mathbf{u}},\theta )$ – это аппроксимации на текущем слое по времени всех слагаемых уравнений баланса массы, импульса и полной энергии соответственно, кроме слагаемых с производными по времени.

Тогда явный метод Эйлера для решения построенных уравнений гомогенной бинарной смеси можно записать в виде

где $\rho = {{\rho }_{a}} + {{\rho }_{b}}$, а ${{{\mathbf{d}}}_{\alpha }}$ – значение вектора параметров для газа $\alpha $. Неизвестные на новом слое по времени $(\mathop {\hat {\rho }}\nolimits_a ,\mathop {\hat {\rho }}\nolimits_b ,{\mathbf{\hat {u}}},\hat {\theta })$ легко последовательно находятся из выписанных формул. Здесь $E = \tfrac{1}{2}\rho {{\left| u \right|}^{2}} + \left( {{{c}_{{Va}}}{{\rho }_{a}} + {{c}_{{Vb}}}{{\rho }_{b}}} \right)\theta $, $\hat {E} = \tfrac{1}{2}\left( {{{{\hat {\rho }}}_{a}} + {{{\hat {\rho }}}_{b}}} \right){{\left| {{\mathbf{\hat {u}}}} \right|}^{2}} + \left( {{{c}_{{Va}}}{{{\hat {\rho }}}_{a}} + {{c}_{{Vb}}}{{{\hat {\rho }}}_{b}}} \right)\hat {\theta }$.

Отметим, что сложение первых двух формул дает

что соответствует аппроксимации уравнения для суммарной плотности (1.29).

2. РЕГУЛЯРИЗОВАННЫЕ УРАВНЕНИЯ ГОМОГЕННОЙ БИНАРНОЙ СМЕСИ НЕРАЗРЕЖЕННЫХ ГАЗОВ

Приведенные в начале разд. 1 уравнения для описания течений негомогенной бинарной смеси (1.1)–(1.10) позволяют рассчитывать течения умеренно-разреженного газа, см. [6], [7], [10]. Но, как показывает практика численного моделирования, с ростом плотности газов, сопровождающимся увеличением числа взаимных столкновений молекул отдельных компонент, обменные слагаемые в этой системе стремительно возрастают, что приводит к неустойчивости численного алгоритма. Для таких течений удается использовать построенные в разд. 1 уравнения гомогенной бинарной смеси (1.12)–(1.16), (1.18)–(1.21). Однако эти уравнения также теряют свою адекватность с уменьшением числа Кнудсена. Последнее связано с увеличением градиентов давлений отдельных компонент в областях границ между разными компонентами смеси.

Поэтому для моделирования течений неразреженных (достаточно плотных) гомогенных газовых бинарных смесей в данном разделе строится альтернативная регуляризованная модель, в которой в отличие от двух предыдущих систем уравнений из разд. 1 регуляризаторы включают в себя давления и плотности всей смеси, а не отдельных ее компонент. Отметим, что использование единственной регуляризующей скорости смеси ${\mathbf{\hat {w}}}$ было предложено в [13] и успешно реализовано в 3D изотермических расчетах [14] для существенно более простой, чем в данной статье, квазигидродинамической (КГидроД) регуляризации в задачах о течении бинарной (двухкомпонентной) смеси при малых изменениях суммарной плотности, но с дополнительным учетом межфазных эффектов (описываемых уравнениями Навье–Стокса–Кана–Хилларда для сжимаемой жидкости).

Построим регуляризованные (КГД) уравнения гомогенной бинарной смеси неразреженных газов. Уравнение баланса массы компоненты изменим посредством замены ${{{\mathbf{w}}}^{{(\alpha )}}}$ на ${\mathbf{w}}$ в (1.12). Такую же замену сделаем в уравнении баланса полной энергии смеси (1.19). Кроме того, в регуляризующих слагаемых уравнений баланса импульса и полной энергии смеси внесем такие дополнительные изменения, чтобы в дальнейшем было выполнено уравнение баланса энтропии смеси с неотрицательным производством энтропии (это совсем нетривиально в подобной модели). А именно, предлагаемые регуляризованные уравнения гомогенной бинарной смеси неразреженных газов таковы:

В них по-прежнему $\rho = {{\rho }_{a}} + {{\rho }_{b}}$, $p = {{p}_{a}} + {{p}_{b}}$, $E = \tfrac{1}{2}\rho {{\left| {\mathbf{u}} \right|}^{2}} + \rho \varepsilon $ и используются формулы (1.18) для $\varepsilon $. Обратим внимание на введение в уравнение (2.3) нового регуляризующего слагаемого с разностью газовых постоянных $\delta R: = {{R}_{a}} - {{R}_{b}}$ и $C = {{\rho }_{a}}{\text{/}}\rho $ – концентрацией компоненты $a$. Кроме того, используются прежние регуляризующие скорости смеси а также тензор вязких напряжений $\Pi = {{\Pi }^{{NS}}} + {{\tilde {\Pi }}^{\tau }}$ и поток тепла $ - {\mathbf{q}} = ({{\varkappa }_{a}} + {{\varkappa }_{b}})\nabla \theta - {{{\mathbf{\tilde {q}}}}^{\tau }}$ со следующими модифицированными слагаемыми-регуляризаторами: где ${{c}_{V}}$, $R$, $\gamma $ по-прежнему определяются формулами (1.18), (1.22). Модификации подвергнуты все три слагаемых с множителем–тензором $\mathbb{I}$ в (2.5) и слагаемое $R\theta \nabla \rho $ в (2.6). В соответствии с модифицированным слагаемым $\gamma pdiv{\mathbf{u}} = \tilde {с}_{s}^{2}\rho div{\mathbf{u}}$ скорость звука в смеси теперь вычисляется по формуле привычного вида ${{\tilde {с}}_{s}} = \sqrt {\gamma (\gamma - 1)\varepsilon } $. В частном случае ${{R}_{a}} = {{R}_{b}}$ слагаемое c $\delta R$ в уравнении (2.3) исчезает и, как обычно, оказывается $R{\mathbf{u}} \cdot \nabla (p{\text{/}}R) = {\mathbf{u}} \cdot \nabla p$ в (2.5).

Сложение по $\alpha = a,b$ новых уравнений баланса плотностей компонент (2.1) приводит к прежнему уравнению баланса суммарной плотности (1.29). Поэтому и уравнение баланса кинетической энергии смеси (1.30) сохраняет свой вид.

Его вычитание из уравнения баланса полной энергии смеси (2.3) приводит к такому уравнению баланса внутренней энергии смеси

Энтропия смеси по-прежнему задается формулами (1.32).

Теорема 2. Верно следующее регуляризованное уравнение баланса энтропии гомогенной бинарной смеси неразреженных газов

с производством энтропии ${{\mathcal{P}}^{{NS}}} + {{\mathcal{P}}^{\tau }}$, где ${{\mathcal{P}}^{{NS}}}$ дается формулой (1.34), а ${{\mathcal{P}}^{\tau }}$ — формулойпричем ${{\mathcal{P}}^{\tau }} \geqslant 0$ при условии $\tau (\gamma - 1)Q \leqslant 4p$.

Доказательство. В силу аддитивности $\rho s$ и нового уравнения баланса массы компоненты (2.1) в полной аналогии с (1.36), (1.37) (если угодно, опуская верхние индексы у ${{{\mathbf{w}}}^{{(a)}}}$, ${{{\mathbf{w}}}^{{(b)}}}$ и ${{{\mathbf{w}}}^{{(\alpha )}}}$) имеем

Отсюда с учетом аддитивности $p$ и $\rho \varepsilon $ и с помощью уравнения баланса внутренней энергии смеси (2.7) выводим где в силу формулы $pdiv({\mathbf{u}} - {\mathbf{w}}) - pdiv{\mathbf{u}} + div(p{\mathbf{w}}) = \nabla p \cdot {\mathbf{w}}$ имеем

С помощью формул

где учтено, что $p = R\rho \theta $ (см. (1.23)), выполним цепочку следующих преобразований по частичному выделению квадратичных слагаемых:

Далее, используя последовательно формулы

где учтено, что $R = {{c}_{V}}(\gamma - 1)$ (см. (1.22)), а также формулы и $p = {{c}_{V}}(\gamma - 1)\rho \theta $ (см. (1.23)), завершаем выкладки по выделению квадратичных слагаемых следующим образом: Формула (2.8) получена, и теорема доказана.

Остановимся также на существенно более простой КГидроД регуляризации (в том числе для сопоставления с [13]), которая также может использоваться в приложениях, хотя и для существенно более узкого класса течений, при умеренных числах Маха. Соответствующая система уравнений использует только одну регуляризующую скорость ${\mathbf{\hat {w}}}$, см. (2.4), и принимает вид

В ней $\Pi = {{\Pi }^{{NS}}} + {{\tilde {\Pi }}^{\tau }}$, но регуляризующий тензор вязкости намного проще: ${{\tilde {\Pi }}^{\tau }} = \rho {\mathbf{u}} \otimes {\mathbf{\hat {w}}}$, а поток тепла регуляризующего слагаемого не содержит вовсе: $ - {\mathbf{q}} = ({{\varkappa }_{a}} + {{\varkappa }_{b}})\nabla \theta $.

Для этой системы уравнение баланса энтропии смеси выводится и выглядит намного проще:

снова с производством энтропии ${{\mathcal{P}}^{{NS}}} + {{\mathcal{P}}^{\tau }} \geqslant 0$.

3. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ГОМОГЕННЫХ БИНАРНЫХ СМЕСЕЙ НЕРАЗРЕЖЕННЫХ ГАЗОВ В ЗАДАЧЕ ТИПА РЭЛЕЯ-ТЕЙЛОРА

Рассматривается задача о моделировании неустойчивости гомогенной бинарной смеси неразреженных газов в рамках уравнений из разд. 2 в двумерной постановке в координатах $(x,z)$ в поле силы тяжести ${\mathbf{F}} = (0, - g)$ (где $g = 9.82$ м/сек² – ускорение свободного падения) в области шириной $L = 0.125$ м и высотой $H = 0.5$ м. Граница раздела двух газов изначально находится вдоль прямой $z = 0$, при этом в верхней части области $0 < z < {{H}_{2}} = 0.2$ располагается тяжелый газ, в нижней ${{H}_{1}} = - 0.3 < z < 0$ – легкий газ.

Берутся газы с параметрами, указанными в табл. 1. Коэффициенты вязкости и теплопроводности таковы $\mu = {\text{const}}$, $\lambda = 0$, $\varkappa = (\gamma {\text{/}}(\gamma - 1)){{R}_{0}}\mu {\text{/}}(M{{\alpha }_{{Pr}}})$. Напомним, что ${{R}_{0}}$ – универсальная газовая постоянная, а $M$ – молекулярный вес газа.

Разностная аппроксимация уравнений строилась согласно [7]. Параметр регуляризации и шаг по времени вычислялись по формулам

где ${{\tilde {с}}_{s}} = \sqrt {\gamma (\gamma - 1)\varepsilon } $ – скорость звука смеси, ${{h}_{x}}$ и ${{h}_{z}}$ – шаги сетки по $x$ и $z$, а $min$ из двух величин берется по всем узлам сетки по пространству. Здесь $\alpha $ и $\beta $ (число Куранта) – параметры из $(0,1)$.

Начальные условия. Начальное состояние газа в момент времени $t = 0$ определяется гидростатическим равновесием в предположении адиабатичности: ${{\partial }_{z}}p = - g\rho $, ${\mathbf{u}} = 0$. Для него зависимости плотности и давления компонент газа от высоты $z$ задаются с помощью формул

где $\gamma _{\alpha }^{'}: = {{\gamma }_{\alpha }}{\text{/}}({{\gamma }_{\alpha }} - 1)$; температура $\theta $ может быть вычислена по ним в соответствии с уравнениями состояния газов. При этом для $z = 0$ задаются параметры нормальной атмосферы, а именно, ${{p}_{a}} = {{p}_{b}} = {{p}_{0}} = 101\,325$ Па, $\theta $ = 273 K, а плотности ${{\rho }_{{a0}}}$ и ${{\rho }_{{b0}}}$ вычисляются по ним в соответствии с уравнениями состояния газов. Отметим, что значения ${{p}_{a}}(z)$, ${{p}_{b}}(z)$, а также $\theta (z)$ слабо отличаются от указанных их значений при $z = 0$.

В той части области, где в начальный момент один из газов практически отсутствует, его плотность берется малым значением ${{10}^{{ - 8}}}$.

Для инициализации процесса в окрестности границы раздела газов задается осциллирующее симметричное по $x$ и большое по амплитуде возмущение вертикальной скорости ${{\left. {{{u}_{z}}} \right|}_{{0 \leqslant z \leqslant {{z}_{0}}}}} = 10cos(4\pi x{\text{/}}L)$, где ${{z}_{0}} = 0.0025$.

Отметим, что значения молекулярных весов $M$ необходимо брать с достаточно высокой точностью, чтобы при $t = 0$ с надлежащей точностью выполнялись условия гидростатического равновесия для обоих газов. Иначе слои газа сразу приходят в движение.

Граничные условия. На вертикальных (боковых) и нижней границах расчетной области ставятся условия непротекания и скольжения

где ${\mathbf{u}} = ({{u}_{x}},{{u}_{z}})$. Для обеспечения условия непротекания потока массы на нижней стенке дополнительно берется ${{\left. {\mathbf{F}} \right|}_{{z = - {{H}_{1}}}}} = 0$.

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

где ${{p}_{{{{H}_{2}}}}}$ найдено при $z = {{H}_{2}}$ согласно (3.1). Такие условия позволяют возмущениям решения покидать область расчета через верхнюю границу.

Результаты расчетов. Расчеты выполнены для трех смесей (пар) газов, часто рассматриваемых в литературе (см., например, [15], [16]) и характеризуемых значениями чисел Атвуда $At = ({{M}_{a}} - {{M}_{b}}){\text{/}}({{M}_{a}} + {{M}_{b}})$, равными 0.091, 0.57, 0.947. То есть охвачены все характерные значения $At$, включая наиболее сложный режим $At \approx 1$.

На фиг. 1 представлены распределения суммарной плотности $\rho $ для смеси Ar–O₃ с $At = 0.091$ на три характерных момента времени $t = 0.1,0.2,0.3$ с. Распределения температуры смеси, а также плотностей отдельных компонент, в значительной мере повторяют распределения плотности смеси и поэтому для краткости не приводятся.

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

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

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

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

В силу неустойчивости задачи Рэлея–Тейлора результаты ее моделирования очень сильно зависят от многих параметров, включая число Атвуда, показатели адиабат ${{\gamma }_{a}}$ и ${{\gamma }_{b}}$, форму и вид начального возмущения данных, конкретный момент времени и т.д. Поэтому имеется возможность только качественно сопоставить полученные результаты с результатами других авторов. Тем не менее можно отметить, что полученные распределения полей плотности характерны для такого типа задач и, в частности, аналогичны полям, представленным в [16] для числа $At = 0.1$.

Базовые расчеты проведены на сетке 125$ \times $400 узлов при параметрах $\alpha = 0.2$, $\beta = 0.3$. Отметим, что оптимизация этих величин в данных расчетах не выполнялась и возможно, что коэффициент $\alpha $, определяющий величину искусственной диссипации, может быть уменьшен, а число Куранта увеличено.

На фиг. 2 представлены результаты расчетов задачи о смеси Ar–SF₆ с $At = 0.57$ на моменты времени $t = 0.1,0.15,0.2$ с. Сопоставление с предыдущим расчетом показывает, что в целом общая картина развития процесса неустойчивости сохраняется. Однако с увеличением $At$ процесс формирования пузырей-грибов происходит более интенсивно, величина акустических возмущений возрастает, проваливающиеся вниз пузыри тяжелого газа SF₆ становятся более узкими и со временем формируют тонкие струи, верхние части которых стремятся оторваться от основного потока и сформировать отдельные капли. При этом всплывающие пузыри легкого газа становятся более широкими. Такая форма пузырей соответствует известным данным [5], [15], [16].

Дополнительные результаты, полученные в расчетах смеси Воздух–SF₆ с $At = 0.66$, оказались в целом близки к последним и поэтому здесь опущены.

Результаты расчетов задачи для смеси He–SF₆ с $At = 0.947$ представлены на фиг. 3. Наблюдается бурное развитие неустойчивости, формирование грибообразных пузырей всплывающего газа и зон опускания тяжелого газа происходит с существенно большей скоростью, чем для смесей с меньшими числами $At$, поэтому распределения плотности смеси представлены на более ранние моменты времени $t = 0.06,0.10,0.12$ с. Отметим, что расчеты подобных задач представляют серьезные трудности из-за огромной разницы в плотностях газов и большой скорости протекающих процессов.

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

На изолиниях температуры наблюдается стремительный нагрев гелия за счет его сжатия под действием тяжелого газа, играющего роль поршня. Для бó льших $t$ происходит взаимодействие образовавшихся пузырей тяжелого газа с нижней границей области, и картина течения перестает быть характерной для задачи гравитационной неустойчивости. Эти расчеты проведены снова на сетке 125 $ \times $ 400 узлов, при $\alpha = 0.4$ и меньшем $\beta = 0.03$.

Для анализа зависимости численного решения от шага сетки по пространству на фиг. 4, 5 показаны результаты для первого варианта смеси, т.е. Ar–O₃ с $At = 0.091$, на те же моменты времени $t = 0.1,0.2,0.3$ с, для сеток 60 $ \times $ 200 и 250 $ \times $ 800 узлов соответственно. Результаты наглядно показывают, что даже на грубой сетке характерные особенности течения разрешаются, за исключением специфической загибающейся формы границы грибов, поскольку на эту границу приходится всего одна или несколько ячеек сетки. При сгущении сетки выявляются дополнительные более тонкие эффекты формирования грибообразных структур. Кроме того, полученные в расчетах на подробной сетке структуры закрученных течений оказываются более округлыми и обладают большим числом витков спирали, чем удается наблюдать на более грубой сетке.

4. ЗАКЛЮЧЕНИЕ

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

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

На основе последней системы уравнений реализован численный метод. Он позволил выполнить расчеты задачи о неустойчивости типа Рэлея–Тейлора во всем реальном диапазоне числа Атвуда, начиная с небольших и вплоть до близких к его максимальному значению 1. Полученные результаты качественно соответствуют известным в литературе. Численное решение обладает надлежащей симметрией. Показано улучшение его качества с измельчением сетки по пространству.

Численный алгоритм уже реализован и для расчета 3D течений; более того, описанные выше 2D расчеты выполнены по 3D коду с минимальным числом узлов по третьему из направлений.

Первая версия алгоритма на основе регуляризованных уравнений однокомпонентного газа недавно была включена в состав открытого программного комплекса OpenFOAM [21]. Разработанный в данной статье алгоритм для моделирования течений газовых смесей также может быть включен в этот комплекс.

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