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

ВВЕДЕНИЕ

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

Реалистичные химические модели, как правило, включают большое число компонент и реакций, причем последние нередко являются жесткими. Это значит, что их временные масштабы много меньше характерных времен процессов конвекции. Кроме того, например, в детонационных волнах (ДВ) химические реакции протекают в узких зонах за этими волнами. Аккуратное разрешение всех этих масштабов при численном моделировании требует использования очень подробных сеток, подразумевающих чрезвычайно большие (иногда неприемлемые) объемы вычислений в многомерном случае. Примечательно, что эта проблема обусловлена в большей степени жесткостью химических реакций и в меньшей степени их числом и количеством компонент газа [1], [2]. Из сказанного следует, что численная схема для динамики идеальных и вязких газовых смесей с химическими реакциями должна корректно воспроизводить главные структуры течения, к примеру, те же ДВ, на относительно грубых сетках.

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

В подавляющем большинстве численных схем для уравнений Эйлера и Навье–Стокса с химическими источниками применяется расщепление по физическим процессам на конвекцию, диффузию (при наличии вязкости) и реакции. Для расчета конвективной части чаще всего используются явные схемы сквозного счета. Со времен работы [1] известно, что численная диссипация таких схем вкупе с жесткостью реакций и грубостью сетки приводит к физически неверной расчетной эволюции ДВ. Вместо одной истинной сильной ДВ схема порождает две ложные волны: ударную волну (УВ) без реакций и ослабленную ДВ (прекурсор), распространяющуюся со скоростью порядка $h{\text{/}}\tau $, где $h$ и $\tau $ – шаги по пространству и времени соответственно. Для борьбы с этим нежелательным явлением было разработано немало разных приемов: переключение [3] и рандомизация [4], [5] температуры реакции; кусочно-постоянная реконструкция MinMax [6] с разрывом внутри ячейки; метод подсеточного разрешения [7], [8], основанный на оценке положения фронта ДВ и экстраполяции температуры и массовых долей реагентов из соседних ячеек; метод выделения равновесных состояний в ячейке [9]; метод пороговых значений [10], явно ограничивающий локальную скорость численной ДВ через приближенную оценку скорости реальной ДВ; монотонная и свободная от ограничителей реконструкция на гиперболических тангенсах [11]. Детальный обзор разнообразных методов обеспечения правильной скорости ДВ и их хронология изложены в работе [7].

Отметим также работы [12], [13], развивающие оригинальный подход к моделированию движений единичной плоской ДВ. Его сутью является явное выделение фронта ДВ в качестве границы вычислительной области. Скорость и координата фронта находятся в процессе счета. Хотя подход [12], [13] уступает схемам сквозного счета в общности, он превосходит их в точности отслеживания положения фронта и его скорости: в [12], [13] эти величины могут быть посчитаны с любой точностью, в то время как схемы сквозного счета принципиально не могут локализовать фронт ДВ с погрешностью меньше $h$.

Как уже было указано выше, наибольшей популярностью при расчетах конвекции пользуются именно явные схемы сквозного счета. Им свойственна численная диссипация, провоцирующая преждевременное срабатывание реакций, поэтому ее стараются уменьшить за счет малости шага по времени. Однако, чем меньше отношение $\tau {\text{/}}h$, тем активнее и быстрее ложная ДВ (прекурсор) стремится оторваться от истинной ДВ (скорость прекурсора имеет порядок $h{\text{/}}\tau $ [1]). Это явление пытаются подавить при помощи достаточно искусственных приемов. Разумно задаться следующим вопросом: почему бы не считать конвекцию по неявным схемам сквозного счета при больших шагах по времени? По-видимому, в этом случае необходимость в процедурах подавления отпадет, если, конечно, численная диссипация не будет слишком значительной.

В последние годы существенное развитие получил новый класс высокоточных схем – бикомпактные схемы [14]–[17]. Применительно к обсуждаемой проблеме они обладают двумя актуальными достоинствами: во-первых, они неявные ($A$-устойчивые и жестко-точные по пространственным переменным [18] и могут быть $L$-устойчивыми и жестко-точными по времени) и при этом экономичные [15], во-вторых, имеют относительно невысокую численную диссипацию даже при числах Куранта, близких к единице (примеры см. в [16]). Кроме того, пространственный шаблон этих схем является компактным и занимает одну ячейку сетки [14], [18], бикомпактные схемы обладают уникальными дисперсионными свойствами [18]. Бикомпактные схемы также допускают счет на декартовых сетках с адаптацией к решению (см. [17]). Адаптивное измельчение позволяет, с одной стороны, экономично разрешать мелкие масштабы, связанные с химией, и, с другой, увеличивать локальное отношение $\tau {\text{/}}h$.

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

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

1. МОДЕЛЬ

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

где ${{X}_{1}}$ и ${{X}_{2}}$ – условные обозначения химических формул первой и второй компонент газа соответственно. Реакция (1) представляет собой простейшую модель горения.

Одномерные течения такого газа описываются системой уравнений Эйлера:

где искомый вектор консервативных переменных ${\mathbf{Q}} = ({{Q}_{1}}, \ldots ,{{Q}_{m}}) = {\mathbf{Q}}(x,t)$, вектор потоков ${\mathbf{F}}({\mathbf{Q}})$ и вектор источников ${\mathbf{S}}({\mathbf{Q}})$ определяются в виде Символ ${{\partial }_{x}} \equiv \partial {\text{/}}\partial x$. Через $\rho $, $v$, $p$, $E$, ${{Z}_{1}}$, $T$ обозначены соответственно плотность, скорость, давление, энергия единицы объема, массовая доля первой компоненты и температура газа. Последний считается совершенным, параметр $\gamma = {\text{const}}$ – показатель адиабаты. Параметр ${{q}_{1}} = {\text{const}}$ равен суммарному тепловому эффекту от участия компоненты ${{X}_{1}}$ во всех реакциях, протекающих в газе; в нашем случае ${{q}_{1}}$ совпадает с тепловым эффектом реакции (1).

Реакция (1) характеризуется своей скоростью ${{K}_{1}}(T)$, для задания которой мы применяем так называемую “кинетику Хэвисайда”:

где ${{B}_{1}} = {\text{const}}$ – параметр с размерностью ${{K}_{1}}$, $\theta ( \cdot )$ – функция Хэвисайда, ${{T}_{1}}$ – температура реакции (1) (температура воспламенения). Важно уточнить, что реакция (1) полагается жесткой: ее число Дамкелера где ${{V}_{0}}$ – характерная скорость процессов конвекции.

Отметим, что хотя “кинетика Хэвисайда” (3) не является полностью корректной с физической точки зрения, она почти не отличается от кинетики Аррениуса на ДВ большой амплитуды. Это объясняется двумя факторами. Во-первых, перепад температур на такой ДВ значителен, он сосредоточен в очень узкой зоне вблизи самой волны и носит резкий, скачкообразный характер. Следовательно, от экспоненты Аррениуса остаются практически одни лишь ее предельные значения 0 и 1. Во-вторых, даже промежуточные значения экспоненты Аррениуса, будучи умноженными на достаточно большой параметр ${{B}_{1}}$, все равно дают высокие скорости реакции, удовлетворяющие неравенствам вида (4). Поскольку предметом настоящей работы будет расчет ДВ большой амплитуды, мы считаем возможным ограничиться более простой кинетикой (3).

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

где искомый вектор консервативных переменных ${\mathbf{Q}} = ({{Q}_{1}}, \ldots ,{{Q}_{m}}) = {\mathbf{Q}}(x,y,t)$, векторы потоков ${\mathbf{F}}({\mathbf{Q}})$ и ${\mathbf{G}}({\mathbf{Q}})$, вектор источников ${\mathbf{S}}({\mathbf{Q}})$ выражаются по формулам Здесь ${\mathbf{v}} = ({{v}_{x}},{{v}_{y}})$ – вектор скорости.

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

2. ЧИСЛЕННАЯ СХЕМА

Применим для численного решения системы уравнений (5) схему расщепления по физическим процессам Марчука–Стрэнга второго порядка [19], [20]. Нетрудно видеть, что система (5) состоит из конвективной части

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

Для расчета конвективной части (6) мы воспользуемся одной бикомпактной схемой четвертого порядка аппроксимации по $x,y$ и третьего порядка аппроксимации по $t$ [14]. Интегрирование по $t$ в ней осуществляется при помощи $L$-устойчивого жестко-точного трехстадийного SDIRK-метода [21, теорема 5]. Выбранная высокоточная бикомпактная схема задействует три вычислительных приема (метода).

1. Симметризованное локально-одномерное расщепление по $x,y$ [19], [22], [23] (см. также [24]).

2. Метод консервативной монотонизации бикомпактных схем [16]. В качестве монотонной схемы-партнера для высокоточной бикомпактной схемы мы возьмем “неявный левый уголок”.

3. Декартовы сетки с адаптацией к решению. Реализация бикомпактных схем на таких сетках исчерпывающе изложена в работе [17].

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

Зададим сразу несколько параметров, значения или способ нахождения которых не будем варьировать в зависимости от задачи: ${{\omega }_{0}} = 2$, ${{\omega }_{1}} = 1$, ${{\omega }_{2}} = 0.1$ (как в [17]); $\sigma = 0.1$ (как в [16]); $C_{2}^{x}$, $C_{2}^{y}$, $\tau $ вычисляются автоматически исходя из данных на предыдущем слое, (см. [17]), $\delta = 0.2$; ${\text{rtol}} = {{10}^{{ - 7}}}$. Значения ${{N}_{x}}$, ${{N}_{y}}$, ${{R}_{{max}}}$, ${{C}_{1}}$, числа Куранта $\kappa $ мы будем указывать непосредственно при описании постановок расчетов. Обратим внимание на то, что, в отличие от работы [17], параметр ${{C}_{1}}$ модифицируется в зависимости от ранга ячейки; во всех формулах ${{C}_{1}}$ заменяется на ${{C}_{1}} \cdot {{2}^{R}}$.

Обозначим оператор временного послойного перехода бикомпактной схемы, по которой рассчитывается конвективная часть (6) системы (5), как ${\text{A}}(\tau )$.

Расчет химической части (7) можно выполнять в каждом узле пространственной сетки независимо от других узлов. Для интегрирования уравнения (7) по $t$ мы будем применять $L$-устойчивый жестко-точный двухстадийный DIRK-метод второго порядка с таблицей Бутчера [21, теорема 5]:

Заметим, что источники отличны от нуля только в уравнениях относительно $\rho {{Z}_{i}}$. Следовательно, при счете химической части можно оставлять постоянными значения консервативных переменных $\rho $, $\rho {\mathbf{v}}$, $E$ и итерировать только переменные $\rho {{Z}_{i}}$.

Обозначим оператор послойного перехода метода (8), по которой рассчитывается химическая часть (7) системы (5), как ${\text{R}}(\tau )$.

Запишем теперь результирующий оператор перехода схемы расщепления Марчука–Стрэнга:

Как можно видеть из формулы (9), при расчете химической части шаг $\tau $ дробится на ${{N}_{r}}$ частей. Это делается для того, чтобы не привязывать шаг при счете конвективной части к мелкому шагу, который необходим для достижения разумной точности при счете химической части. Хотя последняя рассчитывается по формально $L$-устойчивому методу, это отнюдь не значит, что шаг для химии можно задавать сколь угодно большим, так как это негативно отразится на фактической точности решения.

Величину ${{N}_{r}}$ будем вычислять автоматически по формуле

где ${{z}_{{{\text{ODE}}}}}$ – задаваемый безразмерный параметр схемы интегрирования уравнения (7).

3. ВЫЧИСЛИТЕЛЬНЫЕ ПРИМЕРЫ

Ниже мы рассмотрим несколько задач Римана о распаде произвольного разрыва для систем уравнений (2) и (5). Всюду далее показатель адиабаты $\gamma = 1.4$.

Задача 1. Постановку этой задачи мы берем из работы [10]. Решается система одномерных уравнений (2) при ${{x}_{{max}}} = 30$, ${{t}_{{max}}} = 1.5$. Реакция (1) жесткая и имеет такие параметры: ${{B}_{1}} = {{10}^{4}}$, ${{T}_{1}} = 2$, ${{q}_{1}} = 20$. Начальное условие задается так:

Граничные условия постоянные: ${\mathbf{Q}}(0,t) = {\mathbf{Q}}(0,0),{\mathbf{Q}}({{x}_{{max}}},t) = {\mathbf{Q}}({{x}_{{max}}},0).$

Распределение (10) имеет простой смысл: левее координаты $x = 10$ находится полностью сгоревший газ при температуре $T = 10 > {{T}_{1}}$, а правее – полностью не сгоревший газ при температуре $T = 1 < {{T}_{1}}$. Известно (см. [10]), что в результате распада начального профиля (10) образуются три структуры: волна разрежения (ВР), контактный разрыв (КР) и ДВ. ВР движется налево, а КР и ДВ – направо.

Посчитаем эту задачу по схеме, описанной в разд. 2. Расчеты проведем как на равномерной сетке, так и сетке с адаптацией к решению. Для обеих сеток ${{N}_{x}} = 300$ (шаг ${{h}_{x}}(0) = 0.1$, как у схем в [10]), ${{C}_{1}} = 2$, ${{z}_{{{\text{ODE}}}}} = 10$. Вычисления на равномерной сетке (${{R}_{{max}}} = 0$) выполним при $\kappa = 2$, а на адаптивной – при ${{R}_{{max}}} = 2$ и $\kappa = 4$ (число Куранта в ячейках нулевого ранга равно 1). Отметим, что мы в полной мере используем временную неявность и $L$-устойчивость бикомпактной схемы для конвективной части и численного метода (8) для химической части, так как число Куранта в разы больше 1, а параметр ${{z}_{{{\text{ODE}}}}}$ велик.

На фиг. 1 приведены посчитанные профили температуры в конечный момент времени. Сплошной кривой нарисовано решение, вычисленное на подробной равномерной сетке с ${{N}_{x}} = 6000$, с ним сопоставляются решения на сетке с ${{N}_{x}} = 300$, изображенные пунктиром и маркерами (маркеры ставятся на каждом пятом целом узле ячеек нулевого ранга). Видно, что схема расщепления (9) на основе бикомпактной схемы обеспечивает хорошее качество решения. Адаптация сетки до ранга 2 несколько улучшает разрешение узкого пика ЗНД (Зельдовича–Неймана–Деринга) за ДВ, КР и “полочки” между КР и ДВ, но в целом она не приводит к визуально заметному повышению точности, поскольку решение составлено преимущественно из констант и скачков. В течение расчетов на равномерной сетке $\tau \approx 0.02$, ${{N}_{r}} \approx 25$, а на адаптивной – $\tau \approx 0.01$, ${{N}_{r}} \approx 12$.

Результат, полученный в задаче 1, позволяет сравнить неявную схему расщепления (9) на основе бикомпактной схемы и явные схемы расщепления на основе таких популярных схем, как MUSCL и WENO5 (см., например, [8], [10]). Из этого сравнения мы можем заключить, что наш подход обладает двумя существенными преимуществами:

1. Явные схемы расщепления на основе MUSCL и WENO5 неверно передают динамику ДВ на не слишком подробных сетках: истинная ДВ распадается на ложную УВ и сильно опережающую ее ложную слабую ДВ (прекурсор). Из-за этого данные схемы требуют либо значительного сгущения сетки, нежелательного в многомерном случае, либо специальных процедур, связанных с оценками скорости фронта ДВ, его явным выделением, экстраполяцией температуры и массовых долей компонент газа из соседних ячеек и тому подобных приемов (см. Введение). Неявная схема расщепления на основе бикомпактной схемы позволяет сразу получить корректное решение, не требуя каких бы то не было коррекций и поправок, приспособленных лишь к узкой конкретной постановке.

2. Шаг по времени в неявной схеме (9) больше в 100–200 раз ($0.01$ или $0.02$ против ${{10}^{{ - 4}}}$ в [10]) при том же ${{N}_{x}}$, что говорит о более высокой скорости счета даже с поправкой на неявность.

Задача 2. Ее постановка отличается от задачи 1 только одним числом: начальная скорость $v(x,0)$ при $x < 10$ равна не 2, а 4 (см. [10]). Начальный разрыв распадается на УВ, движущуюся налево, а также на КР и ДВ, движущиеся направо.

Посчитав эту задачу по схеме (9) на тех же сетках при тех же параметрах, мы получим профили температуры в конечный момент времени, изображенные на фиг. 2. Очевидно, бикомпактная схема в составе схемы расщепления (9) выдает корректное решение, весьма близкое к точному. Значения $\tau $ и ${{N}_{r}}$ получаются такими же, как в задаче 1.

Интересно испытать нашу схему, варьируя параметры реакции (1) в рамках задачи 2.

Начнем со скорости реакции: сделаем дополнительные расчеты при ${{B}_{1}} = {{10}^{2}},\;{{10}^{3}},\;{{10}^{5}}$ при тех же параметрах схемы. Результаты этих расчетов в виде профилей температуры в конечный момент времени показаны на фиг. 3. Расчеты при меньших ${{B}_{1}}$ позволяют лишний раз убедиться в правильности решения, посчитанного ранее при ${{B}_{1}} = {{10}^{4}}$, а расчет при большем ${{B}_{1}}$ (на порядок) – в том, что при увеличении жесткости системы уравнений наша схема не перестает вести себя правильно.

Попробуем теперь поменять тепловой эффект реакции (1): посчитаем задачу 2 при ${{q}_{1}} = 10$ и ${{q}_{1}} = 30$, оставив ${{B}_{1}} = {{10}^{4}}$; параметры схемы возьмем те же, что и раньше. Искомые профили температуры в конечный момент времени приведены на фиг. 4. Ясно, что неявная схема расщепления (9) на основе бикомпактной схемы проходит и это испытание.

Наконец, подвергнем нашу схему последней, самой сложной проверке. Положив параметры ${{B}_{1}} = {{10}^{4}}$ и ${{q}_{1}} = 20$ (исходные значения), понизим температуру реакции: ${{T}_{1}} = 1.1$. При такой ${{T}_{1}}$ даже небольшое повышение температуры по сравнению с состоянием перед ДВ спровоцирует срабатывание реакции. Не меняя параметров схемы, мы получим неправильный результат: ДВ все же распадется на УВ и опережающую слабую ДВ. Однако этот распад будет незначительным: УВ и слабая ДВ к моменту времени $t = 1.5$ будут находиться друг от друга на расстоянии примерно 2.5. Для сравнения, в явной схеме расщепления на основе MUSCL это расстояние составит больше 30 [10]. Увеличив число Куранта κ до 3 на равномерной сетке и до 8 на адаптивной сетке (в ячейках нулевого ранга оно составит 2), мы получим правильное решение. В этом можно убедиться по профилям температуры в конечный момент времени, изображенным на фиг. 5.

Задача 3. Рассмотрим двумерную версию задачи 1 с цилиндрической симметрией. Решается система уравнений (5) при ${{x}_{{max}}} = {{y}_{{max}}} = 30$, ${{t}_{{max}}} = 2$. Параметры реакции (1) такие же, как и в задачах 1, 2. Начальное условие имеет вид

На границах $x = 0,\;y = 0$ ставятся условия симметрии, на границах $x = {{x}_{{max}}},\;y = {{y}_{{max}}}$ – постоянные условия, или фон.

Вычисления выполним на адаптивной сетке с ${{N}_{x}} = {{N}_{y}} = 75$ (${{h}_{x}}(0) = {{h}_{y}}(0) = 0.4$) при ${{R}_{{max}}} = 1$ и $2$; параметры ${{C}_{1}} = 2$, ${{z}_{{{\text{ODE}}}}} = 10$. Число Куранта $\kappa = 2$ для ${{R}_{{max}}} = 1$ и $\kappa = 4$ для ${{R}_{{max}}} = 2$ (число Куранта в ячейках нулевого ранга равно 1 в обоих случаях).

На фиг. 6, 7 приведены двумерные графики температуры и изображения сетки в конечный момент времени при ${{R}_{{max}}} = 1,2$ соответственно. Значения температуры показаны цветовой заливкой, а изображения сетки составлены только из расчетных ячеек. Хорошо видно, что схема расщепления (9) на основе бикомпактной схемы рассчитывает динамику ДВ без проблем. Последняя не распадается на ложные волны и, кроме того, имеет правильную круговую форму. Заметим, что более глубокая адаптация сетки к решению отчетливее проявляет узкий ЗНД-пик за ДВ (темная полоса, примыкающая к ДВ), за которым следует некоторый круговой слой немонотонностей. Впрочем, они достаточно быстро подавляются схемой. Сетка адаптируется в тех слоях, где локализованы ВР и ДВ, что представляется разумным. Укажем дополнительно, что в течение расчетов параметры $\tau \approx 0.05$, ${{N}_{r}} \approx 50$.

ЗАКЛЮЧЕНИЕ

Для системы двумерных газодинамических уравнений Эйлера с жесткими химическими источниками построена неявная схема расщепления. Эта схема использует классическое расщепление Марчука–Стрэнга по физическим процессам. Конвективная часть системы рассчитывается по бикомпактной схеме четвертого порядка аппроксимации по пространству и третьего порядка аппроксимации по времени. Эта высокоточная бикомпактная схема является $L$-устойчивой по времени; в ней применяется метод консервативной монотонизации и адаптивное измельчение (декартовой) сетки в зависимости от локального поведения решения. Химическая часть системы рассчитывается по $L$-устойчивой SDIRK-схеме второго порядка аппроксимации.

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

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