1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 61, № 1, 2021

Журнал вычислительной математики и математической физики, 2021, T. 61, № 1, стр. 124-135

Численное моделирование газовых смесей в рамках квазигазодинамического подхода на примере взаимодействия ударной волны с пузырьком газа

Т. Г. Елизарова 1, *, Е. В. Шильников 1, **

1 ИПМ им. М.В. Келдыша РАН
125047 Москва, Миусская пл., 4, Россия

* E-mail: telizar@mail.ru
** E-mail: shilnikov@imamod.ru

Поступила в редакцию 21.04.2020
После доработки 02.07.2020
Принята к публикации 18.09.2020

Полный текст (PDF)

Аннотация

Представлен новый численный алгоритм для моделирования течений нереагирующих газовых смесей в трансзвуковых режимах. Алгоритм основан на методе конечного объема, записанного для регуляризованных, или квазигазодинамических, уравнений. Уравнения для описания течения смеси выведены феноменологически на базе существующей регуляризованной системы для однокомпонентного газа и классических уравнений для газовой смеси. Примеры численного моделирования включают в себя расчет задачи о нестационарном взаимодействии газового потока с тяжелой и легкой каплями газа. Библ. 23. Фиг. 9. Табл. 1.

Ключевые слова: смесь газов, квазигазодинамические уравнения, метод конечного объема, трансзвуковые течения, газовые капли.

ВВЕДЕНИЕ

Полученные более тридцати лет назад квазигазодинамические (КГД), или регуляризованные, уравнения для описания течений газов различной природы успешно применялись для построения численных алгоритмов как для традиционных вычислительных систем, так и для систем с массивным параллелизмом (см. [1]–[3]). Оценки показывают, что КГД-подход особенно эффективен для моделирования нестационарных или неустановившихся течений газа. В настоящее время распространение этого подхода дополнительно поддерживается тем фактом, что КГД-алгоритм включен в открытый программный комплекс OpenFOAM и доступен широкому кругу пользователей как в России, так и за рубежом (см. [4]–[6]).

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

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

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

1. СИСТЕМА УРАВНЕНИЙ ДЛЯ ОПИСАНИЯ ТЕЧЕНИЯ ГАЗОВОЙ СМЕСИ

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

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

Такая модель описания течений смеси газов является весьма распространенной. В частности, она применяется в ряде работ Р. Абгралла и С. Карни, например, в [7]–[10]. При использовании этого вида моделей предпринималось много усилий для подавления нефизичных осцилляций на границе раздела флюидов при расчетах по схемам повышенного порядка аппроксимации. Для этого применялась так называемая Double flux модификация широко распространенных численных методов, таких как TVD-MUSCL схемы (см. [9]), разрывный метод Галеркина (см. [11]), схемы Годунова (см. [12], [13]). Отметим, что в отличие от описанной выше модели, в [12] рассматривается многожидкостная модель, когда каждый флюид имеет не только свою собственную плотность, но и свои скорости и давления. Это приводит к появлению обменных членов в правых частях уравнений. Такая модель хорошо справляется с моделированием течений смеси жидкостей и газов. Однако наличие обменных членов приводит к необходимости интегрирования на каждом шаге, вообще говоря, жесткой системы ОДУ, что сильно усложняет алгоритм. В случае моделирования течений реагирующих газов жесткую систему ОДУ приходится решать по причине использования в расчетах сильно различающихся скоростей химических реакций (см. [13], [14]), но при отсутствии химических реакций представляется более желательным обойтись без этого усложнения.

Для краткости рассмотрим смесь двух газов a и b. В этом случае описанная выше модель представляется в виде следующей системы уравнений:

(1)
$\frac{{\partial {{\rho }_{a}}}}{{\partial t}} + \operatorname{div} ({{\rho }_{a}}{\mathbf{u}}) = {\text{0,}}$
(2)
$\frac{{\partial {{\rho }_{b}}}}{{\partial t}} + \operatorname{div} ({{\rho }_{b}}{\mathbf{u}}) = {\text{0,}}$
(3)
$\frac{{\partial \rho {\mathbf{u}}}}{{\partial t}} + \operatorname{div} (\rho ({\mathbf{u}} \otimes {\mathbf{u}})) + \nabla p = \operatorname{div} {{\Pi }_{{{\text{NS}}}}} + \rho {\mathbf{F}}{\text{,}}$
(4)
$\frac{{\partial E}}{{\partial t}} + \operatorname{div} ((E + p){\mathbf{u}}) = - \operatorname{div} {{{\mathbf{q}}}_{{{\text{NS}}}}} + \operatorname{div} ({{\Pi }_{{{\text{NS}}}}} \cdot {\mathbf{u}}) + \rho ({\mathbf{u}} \cdot {\mathbf{F}}) + Q.$

В приведенной системе использованы общепринятые обозначения газодинамических величин. Дополнительно введены вектор F и скаляр Q, которые обозначают удельную внешнюю силу и внешний источник или сток энергии соответственно. Последнее будет существенно для случая реагирующих газов, который пока не включен в рассмотрение. Тензор вязких напряжений Навье–Стокса ПNS и вектор теплового потока qNS в форме закона Фурье имеют традиционный вид

${{\Pi }_{{{\text{NS}}}}} = \mu \left( {(\nabla \otimes {\mathbf{u}}) + {{{(\nabla \otimes {\mathbf{u}})}}^{{\text{T}}}} - \frac{2}{3}{\text{I}}\operatorname{div} {\mathbf{u}}} \right),\quad {\mathbf{q}}_{{{\text{NS}}}}^{{}} = - \kappa \nabla T,$
где I – единичная матрица. В выписанной модели предполагается, что газовая смесь имеет единую скорость u и температуру T, а плотность смеси, ее давление и удельная полная энергия определяются через параметры ее компонент как

$\rho = {{\rho }_{a}} + {{\rho }_{b}},\quad p = {{p}_{a}} + {{p}_{b}},\quad E = \rho \varepsilon + {{\rho {{{\mathbf{u}}}^{2}}} \mathord{\left/ {\vphantom {{\rho {{{\mathbf{u}}}^{2}}} 2}} \right. \kern-0em} 2}.$

Кроме того, газодинамическая смесь удовлетворяет обычным соотношениям для идеального политропного газа

$p = \rho RT = \rho \varepsilon (\gamma - 1),$
где γ – показатель адиабаты смеси, R – газовая постоянная и ε – удельная внутренняя энергия смеси:

$R = \frac{{{{R}_{a}}{{\rho }_{a}} + {{R}_{b}}{{\rho }_{b}}}}{\rho } = {{c}_{{\text{p}}}} - {{c}_{{\text{V}}}},\quad \gamma = \frac{{{{c}_{{\text{p}}}}}}{{{{c}_{{\text{V}}}}}},\quad \gamma - 1 = \frac{R}{{{{c}_{{\text{V}}}}}},$
$\varepsilon = \frac{{{{\varepsilon }_{a}}{{\rho }_{a}} + {{\varepsilon }_{b}}{{\rho }_{b}}}}{\rho } = {{c}_{{\text{V}}}}T,\quad {{c}_{{\text{V}}}} = \frac{{{{c}_{{\text{V}}}}_{a}{{\rho }_{a}} + {{c}_{{\text{V}}}}_{b}{{\rho }_{b}}}}{\rho }.$

Подчеркнем, что выписанные выше термодинамические параметры для смеси сp, ${{с}_{{\text{V}}}}$, а следовательно, и R, γ, в отличие от аналогичных величин для идеального политропного газа уже не являются постоянными и определяются через взвешенные значения параметров смеси. Они зависят от плотности каждой компоненты газа в каждой пространственно-временной точке (x, t).

Скорость звука для смеси может быть вычислена с использованием соотношения

$\rho c_{{\text{s}}}^{2} = {{\gamma }_{a}}{{p}_{a}} + {{\gamma }_{b}}{{p}_{b}}.$

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

2. РЕГУЛЯРИЗОВАННАЯ СИСТЕМА УРАВНЕНИЙ ДЛЯ ОПИСАНИЯ ТЕЧЕНИЯ ГАЗОВОЙ СМЕСИ

Регуляризованные (или КГД) уравнения для описания течения вязкого политропного идеального газа имеют следующий вид (см., например, [2], [3]):

(5)
$\frac{{\partial \rho }}{{\partial t}} + \operatorname{div} (\rho ({\mathbf{u}} - {\mathbf{w}})) = {\text{0,}}$
(6)
$\frac{{\partial \rho {\mathbf{u}}}}{{\partial t}} + \operatorname{div} (\rho ({\mathbf{u}} - {\mathbf{w}}) \otimes {\mathbf{u}}) + \nabla p = \operatorname{div} \Pi + {\text{(}}\rho - \tau \operatorname{div} (\rho {\mathbf{u}})){\mathbf{F}}{\text{,}}$
(7)
$\frac{{\partial E}}{{\partial t}} + \operatorname{div} ((E + p)({\mathbf{u}} - {\mathbf{w}})) = - \operatorname{div} {\mathbf{q}} + \operatorname{div} (\Pi \cdot {\mathbf{u}}) + \rho ({\mathbf{u}} - {\mathbf{w}}) \cdot {\mathbf{F}} + Q.$

Эта система уравнений включает в себя регуляризирующие ее добавки вида

(8)
${\mathbf{w}} = \frac{\tau }{\rho }\left( {\operatorname{div} (\rho {\mathbf{u}} \otimes {\mathbf{u}}) + \nabla p - \rho {\mathbf{F}}} \right),$
(9)
${\mathbf{\hat {w}}} = \frac{\tau }{\rho }\left( {\rho ({\mathbf{u}}\nabla ){\mathbf{u}} + \nabla p - \rho {\mathbf{F}}} \right),$
(10)
${\mathbf{q}} = {{{\mathbf{q}}}_{{{\text{NS}}}}} + {{{\mathbf{q}}}^{\tau }},$
(11)
${\text{П}} = {{{\text{П}}}_{{{\text{NS}}}}} + \rho {\mathbf{u}} \otimes {\mathbf{\hat {w}}} + \tau ({\mathbf{u}}\nabla p + \gamma p\operatorname{div} {\mathbf{u}} - (\gamma - 1)Q),$
(12)
$ - {{{\mathbf{q}}}^{\tau }} = \tau \rho {\mathbf{u}}\left( {{\mathbf{u}}\nabla \varepsilon + p({\mathbf{u}}\nabla )\left( {\frac{1}{\rho }} \right) - \frac{Q}{\rho }} \right).$

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

(13)
$\tau = {l \mathord{\left/ {\vphantom {l {{{c}_{{\text{s}}}}}}} \right. \kern-0em} {{{c}_{{\text{s}}}}}},$
где cs – скорость звука и l – некоторый характерный размер, определяемый рассматриваемой задачей. В большинстве вычислительных задач в качестве характерного размера удобно выбрать величину шага пространственной сетки h:
(14)
$\tau = {{\alpha h} \mathord{\left/ {\vphantom {{\alpha h} {{{c}_{{\text{s}}}}}}} \right. \kern-0em} {{{c}_{{\text{s}}}}}},$
с численным коэффициентом $\alpha \leqslant 1$, который выбирается в процессе решения задачи из соображений точности и устойчивости численного алгоритма. Коэффициент τ определяет величину подсеточной диссипации алгоритма и тесно связан с условием устойчивости алгоритма.

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

По аналогии с КГД-уравнениями для однокомпонентного газа (5)–(7) выпишем регуляризованный аналог системы для смеси. Таким же образом, как в системе (1)–(4), положим ρ = ρa + ρb и расщепим уравнение неразрывности (5) на уравнения для каждой из компонент ρa и ρb. Соответствующим образом расщепим и регуляризатор, входящий в исходное уравнение для общей плотности смеси. Принимая во внимание, что компоненты смеси имеют одинаковую скорость, положим, что и регуляризирующие добавки к скоростям для компонент смеси одинаковы: wa = = wb = w. Тогда регуляризованная (или КГД) система уравнений для смеси газов будет иметь следующий вид:

(15)
$\frac{{\partial {{\rho }_{a}}}}{{\partial t}} + {\text{div}}({{\rho }_{a}}({\mathbf{u}} - {\mathbf{w}})) = {\text{0,}}$
(16)
$\frac{{\partial {{\rho }_{b}}}}{{\partial t}} + {\text{div}}({{\rho }_{b}}({\mathbf{u}} - {\mathbf{w}})) = {\text{0,}}$
(17)
$\frac{{\partial \rho {\mathbf{u}}}}{{\partial t}} + {\text{div}}(\rho ({\mathbf{u}} - {\mathbf{w}}) \otimes {\mathbf{u}}) + \nabla p = \operatorname{div} \Pi + {\text{(}}\rho - \tau \operatorname{div} (\rho {\mathbf{u}})){\mathbf{F}}{\text{,}}$
(18)
$\frac{{\partial E}}{{\partial t}} + \operatorname{div} ((E + p)({\mathbf{u}} - {\mathbf{w}})) = - \operatorname{div} {\mathbf{q}} + \operatorname{div} {\text{(}}\Pi \cdot {\mathbf{u}}{\text{) + }}\rho ({\mathbf{u}} - {\mathbf{w}}) \cdot {\mathbf{F}} + Q.$

Эта система дополняется регуляризирующими добавками вида (8)–(12). Для вычисления коэффициента τ через скорость звука cs для смеси используются соотношения (13) или (14).

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

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

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

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

Близкая по своей структуре более простая квазигидродинамическая модель течения смеси была построена в [20] для моделирования медленных двухфазных течений в пористых средах с учетом межфазных взаимодействий. Вид этой системы уравнений имеется в [19].

3. ПОСТАНОВКА ЗАДАЧИ О ВЗАИМОДЕЙСТВИИ ПУЗЫРЬКА ГАЗА С УДАРНОЙ ВОЛНОЙ

Тестовая задача состоит в моделировании взаимодействия плоской ударной волны, движущейся в воздухе, с цилиндрическим пузырем другого газа. Такой эксперимент описан в [21] и численно исследован в [7], [8], [22], [23]. Рассматривается прямоугольная двумерная область, заполненная воздухом (фиг. 1). Плоская ударная волна, проходя через воздух, падает на цилиндрический пузырь из гелия или хладоагента R22 (CHClF2). Пузырь с радиусом $R = 0.025$ помещается в воздух с центром пузыря в точке $({{x}_{c}},{{y}_{c}}) = (0.32,\;0)$. Пренебрегая граничными эффектами, т.е. считая цилиндр бесконечно длинным, можно моделировать такое течение в двумерной постановке. В центральном сечении течение является практически двумерным и краевыми эффектами можно пренебречь, если длина цилиндра не мала по сравнению с высотой установки.

Фиг. 1.

Схема расчетной области.

Все газовые компоненты считаются идеальными газами. Газовая постоянная $R = {{{{R}_{{{\text{univ}}}}}} \mathord{\left/ {\vphantom {{{{R}_{{{\text{univ}}}}}} m}} \right. \kern-0em} m},$ где ${{R}_{{{\text{univ}}}}}$ – универсальная газовая постоянная, m – молекулярная масса газа. Мы пренебрегаем физической вязкостью µ и теплопроводностью κ газов и используем уравнения смеси в формулировке Эйлера. Внешняя сила F и источник тепла Q считаются равными нулю. Начальные динамически равновесные параметры газов в расчетной области приведены в табл. 1.

Таблица 1.  

Начальные параметры газов. Размерности всех величин соответствуют системе СИ

Газ ρ u ${v}$ p γ m R cs
Воздух 1.0 0.0 0.0 105 1.4 28.96 287.1 374.16
Гелий 0.182 0.0 0.0 105 5/3 4.003 2077 915.75
R22 3.1538 0.0 0.0 105 1.249 86.47 96.15 199.0

На правой границе задается условие притока воздуха с параметрами, лежащими за ударной волной, движущейся справа налево через воздух со скоростью, соответствующей числу Маха ${{M}_{{\text{s}}}} = 1.22$:

$\left( {\rho ,u,{v},p,\gamma } \right)\left| {_{{{\text{right}}}}} \right. = \left( {1.3764,\; - {\kern 1pt} {\text{124}}{\text{.82414,}}\;0.0,\;{\text{156983}}{\text{.9256,}}\;1.4} \right).$

Все остальные границы области рассматриваются как твердые непроницаемые стенки с граничными условиями скольжения. Расчеты проводятся на прямоугольных сетках: грубая сетка, состоящая из $1300 \times 178$ ячеек с пространственным шагом $h = 5 \times {{10}^{{ - 4}}}$ в обоих направлениях, и в два раза более подробная сетка, состоящая из $2600 \times 356$ ячеек.

Для численной реализации системы QGD уравнений (15)–(18) мы используем явную по времени схему, построенную на основе метода конечных объемов, с аппроксимацией всех пространственных производных центральными разностями второго порядка. Все газодинамические переменные отнесены к центрам ячеек. Их значения в центрах граней ячеек рассчитываются с использованием линейной интерполяции. Условие устойчивости для этой схемы имеет вид условия Куранта, и шаг по времени определяется по формуле

(19)
$\Delta t = \beta \mathop {\min }\limits_i \frac{{{{h}_{i}}}}{{{{c}_{{{\text{s}}i}}} + \left| {{{{\mathbf{u}}}_{i}}} \right|}},$
где минимум берется по всем ячейкам сетки, i – номер ячейки, β – числовой коэффициент (число Куранта), который не зависит от шага пространственной сетки. Конечное время расчета принимается равным ${{t}_{{{\text{fin}}}}} = 1.4 \times {{10}^{{ - 3}}}$ для пузырька гелия и ${{t}_{{{\text{fin}}}}} = 1.1 \times {{10}^{{ - 3}}}$ для пузырька R22. Время измеряется в секундах и отсчитывается от момента начала движения ударной волны от правой границы расчетной области.

4. АНАЛИЗ РЕЗУЛЬТАТОВ РАСЧЕТОВ

Расчеты для случая пузыря хладоагента R22 проводились с параметром схемы $\alpha = 0.4$ и числом Куранта $\beta = 0.5$. На фиг. 2 численные Шлирен-образы приведены в последовательные моменты времени (сверху вниз): $t = 7.61 \times {{10}^{{ - 4}}}$, $t = 8.47 \times {{10}^{{ - 4}}}$, $t = 1.01 \times {{10}^{{ - 3}}}$, $t = 1.1 \times {{10}^{{ - 3}}}$.

Фиг. 2.

Пузырь R22. Численные Шлирен-образы в последовательные моменты времени (сверху вниз): $t = 7.61 \times {{10}^{{ - 4}}}$, $t = 8.47 \times {{10}^{{ - 4}}}$, $t = 1.01 \times {{10}^{{ - 3}}}$, $t = 1.1 \times {{10}^{{ - 3}}}$.

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

После того как падающая ударная волна наталкивается на пузырь, последний начинает двигаться налево и деформируется (фиг. 2). На верхней картинке, соответствующей времени $t = 7.61 \times {{10}^{{ - 4}}}$, падающая ударная волна состоит из двух вертикальных фрагментов сверху и снизу от пузыря. Внутри пузыря видна преломленная ударная волна, которая вследствие меньшей скорости распространения искривляется. Это связано с тем, что скорость звука в R22 почти в два раза меньше, чем в воздухе (табл. 1). При этом внешние концы преломленной волны движутся вместе с внутренними концами фрагментов падающей волны. Изогнутая отраженная волна, более слабая, чем падающая и преломленная, перемещается вправо от пузыря к правой границе области. Отраженная от пузыря волна соединена с фрагментами падающей волны двумя слабыми волнами, отраженными от горизонтальных стенок. Дифракция волн вокруг пузыря приводит к его деформации. На следующих рисунках на фиг. 2 видны формирование и развитие двух завитков в верхней и нижней частях пузыря, в которых возникает завихренность. Через некоторое время фрагменты падающей ударной волны, миновав пузырь, соединяются и в дальнейшем двигаются влево как целое. Отстающая вначале точка соединения фрагментов падающей волны постепенно догоняет их внешние концы, и фронт волны выпрямляется (нижнее изображение на фиг. 2). Все эти эффекты находятся в хорошем согласии с экспериментом и численными результатами других авторов.

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

Фиг. 3.

Прохождение ударной волны через пузырь R22.

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

Сравнение формы пузырька в последний момент времени, полученное в наших расчетах, с численными из [7] и экспериментальными из [21] результатами показано на фиг. 4. Здесь приведена форма пузыря R22 в момент окончания расчета $t = 1.1 \times {{10}^{{ - 3}}}$: (а) – наши результаты, (б) – результаты из [7], (в) – экспериментальные данные из [21]. Следует отметить, что как эволюция формы пузырька, так и общая структура сложного потока с большим количеством волн разных типов совпадают с результатами этих и других авторов.

Фиг. 4.

Форма пузыря R22 в момент окончания расчета $t = 1.1 \times {{10}^{{ - 3}}}$: (а) – наши результаты, (б) – результаты из [7], (в) – экспериментальные данные из [21].

Расчеты для случая гелиевого пузыря проводились с параметром $\alpha = 0.4$ и числом Куранта β = 0.2. Они представлены на фиг. 5. Здесь приведены численные Шлирен-образы в последовательные моменты времени (сверху вниз): $t = 7.04 \times {{10}^{{ - 4}}}$, $t = 7.52 \times {{10}^{{ - 4}}}$, $t = 9.13 \times {{10}^{{ - 4}}}$, $t = 1.342 \times {{10}^{{ - 3}}}$.

Фиг. 5.

Пузырь гелия. Численные Шлирен-образы в последовательные моменты времени (сверху вниз): $t = 7.04 \times {{10}^{{ - 4}}}$, $t = 7.52 \times {{10}^{{ - 4}}}$, $t = 9.13 \times {{10}^{{ - 4}}}$, $t = 1.342 \times {{10}^{{ - 3}}}$.

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

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

Скорость преломленной волны настолько велика, что на фиг. 5 даже не видно, как она проходит через пузырь. Уже на верхнем рисунке фиг. 5 эта волна достигла левой границы пузыря. Чтобы увидеть процесс движения преломленной волны внутри пузыря, на фиг. 6 приведены фрагменты картины течения в более ранние близкие моменты времени. Моменты времени составляют (слева направо): $t = 6.86 \times {{10}^{{ - 4}}}$, $t = 6.95 \times {{10}^{{ - 4}}}$, $t = 7.03 \times {{10}^{{ - 4}}}$. Видно, что за то время, когда преломленная волна прошла весь пузырь, падающая волна почти не продвинулась.

Фиг. 6.

Прохождение преломленной волны через пузырь гелия. Моменты времени (слева направо): $t = 6.86 \times {{10}^{{ - 4}}}$, $t = 6.95 \times {{10}^{{ - 4}}}$, $t = 7.03 \times {{10}^{{ - 4}}}$.

Еще одно отличие от предыдущего случая состоит в том, что отраженная волна здесь является не ударной волной, а волной разрежения. Это хорошо видно на фиг. 7. На ней изображены профили давления на момент окончания расчета вдоль центрального сечения z = 0 (линии 1) и вдоль сечения z = 0.4 (линии 2), проходящей между стенкой и пузырем для обоих вариантов расчета.

Фиг. 7.

Профили давления вдоль линии z = 0 (линии 1) и z = 0.4 (линии 2): (а) – для пузыря R22, (б) – для пузыря гелия.

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

Сравнение формы пузырька на момент окончания расчета, полученное в наших расчетах с теми же результатами, что и для случая R22, показано на фиг. 8. Здесь показана форма пузыря гелия в момент времени $t = 1.342 \times {{10}^{{ - 3}}}$: (а) – наши результаты, (б) – результаты из [7], (в) – экспериментальные данные из [21].

Фиг. 8.

Форма пузыря гелия в момент времени $t = 1.342 \times {{10}^{{ - 3}}}$: (а) – наши результаты, (б) – результаты из [7], (в) – экспериментальные данные из [21].

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

Сравнение результатов, полученных на разных пространственных сетках, представлены на фиг. 9. На нем приведены результаты расчетов с каплей R22 на момент времени окончания расчета на разных пространственных сетках, а именно, показаны численные Шлирен-образы на момент времени $t = 1.1 \times {{10}^{{ - 3}}}$, полученные на сетке $1300 \times 178$ ячеек (слева) и $2600 \times 356$ ячеек (справа). Видно, что при дроблении сетки общая структура течения сохраняется, а его качество улучшается. Волны становятся тоньше и более четкими, проявляются более тонкие особенности течения, особенно внутри капли. Уменьшение параметра регуляризации α влияет на решение примерно так же, как и дробление сетки. Однако при этом ухудшается устойчивость схемы и требуется уменьшение шага по времени. Поэтому расчеты проводились при минимальном значении α, обеспечивающем устойчивость при разумном значении числа Куранта β.

Фиг. 9.

Численные Шлирен-образы на момент времени $t = 1.1 \times {{10}^{{ - 3}}}$, (а) – полученные на сетке $1300 \times 178$ ячеек и (б) – для $2600 \times 356$ ячеек.

Отметим также, что в работе [7], с которой сравнивались наши результаты, использовалась гораздо более сложная разностная схема на динамически адаптирующейся встроенной сетке, эквивалентная равномерной сетке $16\,000 \times 800$ ячеек, т.е. на порядок более подробная, чем в настоящей работе. Для оценки эффективности работы алгоритма укажем, что, например, расчет задачи о тяжелой капле на подробной сетке $2600 \times 356$ с числом точек порядка миллиона ячеек занимал около 10 ч процессорного времени на персональном компьютере.

ЗАКЛЮЧИТЕЛЬНЫЕ ЗАМЕЧАНИЯ

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

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

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

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

  1. Четверушкин Б.Н. Кинетические схемы и квазигазодинамическая система уравнений. М.: МАКС Пресс, 2004.

  2. Елизарова Т.Г. Квазигазодинамические уравнения и методы расчета вязких течений. М.: Научный мир, 2007. Перевод: Elizarova TG. Quasi-Gas Dynamic Equations. Springer, Berlin, 2009.

  3. Шеретов Ю.В. Динамика сплошных сред при пространственно-временном осреднении. Москва–Ижевск: РХД, 2009. Перевод названия: Sheretov Yu. Continuum Dynamics under Spatiotemporal Averaging. Regular and Chaotic Dynamics, Moscow-Izhevsk, 2009 (in Russian).

  4. Kraposhin M.V., Ryazanov D.A., Smirnova E.V., Elizarova T.G., Istomina M.A. Development of OpenFOAM solver for compressible viscous flows simulation using quasi-gas dynamic equations // IEEE eXplore. 2017. V. 29. № 6. P. 117–124.

  5. Kraposhin M.V., Smirnova E.V., Elizarova T.G., Istomina M.A. Development of a new OpenFOAM solver using regularized gas dynamic equations // Comput. Fluids. 2018. V. 166. P. 163–175.

  6. OpenFOAM User Guide, version 7. https://cfd.direct/openfoam/user-guide/

  7. Quirk J., Karni S. On the dynamics of a shock-bubble interaction // J. Fluid Mech. 1996. V. 318. P. 129–163.

  8. Abgrall R. How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach // J. Comput. Phys. 1996. V. 125. P. 150–160.

  9. Abgrall R., Karni S. Computations of compressible multifluids // J. Comput. Phys. 2001. V. 169. P. 594–623.

  10. Billet G., Abgrall R. An adaptive shock-capturing algorithm for solving unsteady reactive flows // Comput. Fluids. 2003. V. 32. P. 1473–1495.

  11. Billet G., Ryan J. A Runge–Kutta discontinuous Galerkin approach to solve reactive flows: The hyperbolic operator // J. Comput. Phys. 2011. V. 230. P. 1064–1083.

  12. Saurel R., Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows // J. Comput. Phys. 1999. V. 150. P. 425–467.

  13. Borisov V.E., Rykov Yu.G. Modified Godunov method for multicomponent flow simulation // J. Phys.: Conf. Ser. 2019. N1250. P. 012006.

  14. Lian Y.S., Xu K. A Gas-kinetic scheme for multimaterial flows and its application in chemical reactions // J. Comput. Phys. 2000. V. 163. P. 349–375.

  15. Shilnikov E.V., Elizarova T.G. Numerical simulation of gravity instabilities in gas flows by use of the quasi gas dynamic equation system // IOP Conf. Ser.: Materials Science and Engng. 2019. V. 657. P. 012035.

  16. Elizarova T.G., Graur I.A., Lengrand J.C. Two-fluid computational model for a binary gas mixture // Europ. J. Mech. B. 2001. V. 3. P. 351–369.

  17. Елизарова Т.Г., Злотник А.А., Четверушкин Б.Н. О квазигазо- и гидродинамических уравнениях бинарных смесей газов // Докл. АН. 2014. Т. 459. № 4. С. 395–399.

  18. Elizarova T.G., Lengrand J.-C. Free parameters for the modelization of nonmonatomic binary gas mixtures // J. Thermophys. Heat Transfer. 2016. V. 30. № 3. P. 695–697.

  19. Elizarova T.G., Zlotnik A.A., Shil’nikov E.V. Regularized equations for numerical simulation of flows of homogeneous binary mixtures of viscous compressible gases // Comput. Math. and Math. Phys. 2019. V. 59. № 11. P. 1832–1847.

  20. Балашов В.А., Савенков Е.Б. Многокомпонентная квазигидродинамическая модель для описания течений многофазной жидкости с учетом межфазного взаимодействия // Прикл. механ. техн. физ. 2018. Т. 59. № 3. С. 57–68.

  21. Haas J.F., Sturtevant B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities // J. Fluid Mech. 1987. V. 181. P. 41–76.

  22. Иванов И.Э., Крюков И.А. Численное моделирование течений многокомпонентного газа с сильными разрывами свойств среды // Матем. моделирование. 2007. Т. 19. № 12. С. 89–100.

  23. Борисов В.Е., Рыков Ю.Г. Численное моделирование течений многокомпонентных газовых смесей с использованием метода двойного потока // Матем. моделирование. 2020. Т. 32. № 10. С. 3–20.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал