Scientific journal
Fundamental research
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,685

EFFECTIVE NUMERICAL METHOD FOR SIMULATION OF CHEMICALLY NON-EQUILIBRIUM FLOWS

Molchanov A.M. 1
1 Moscow Aviation Institute (National Research University)
This article is devoted to the development of an effective numerical method for solving a system of equations describing the gas flow with nonequilibrium chemical reactions. The main attention is paid to the stiffness of the species mass conservation equations. Several numerical schemes are considered: a purely implicit scheme in which convection, diffusion and a source are presented in implicit form; an explicit-implicit scheme in which convection and diffusion are presented in explicit form, and the source in the implicit form; a scheme with splitting into physical processes; an explicit-implicit scheme with alternation on half-steps. The proposed explicit-implicit numerical scheme is fully coupled. In this case, there is no need for repeated circulation at each time step of the 7-diagonal matrix. Thus, there is no the most laborious part of the solution of the basic system of equations that appears with the implicit convection-diffusion representation and the presence of a large number of additional equations. Test calculations were carried out, which showed the best efficiency of an explicit-implicit scheme with alternation on half-steps. It is also shown that when using a scheme with splitting into physical processes, there are problems with convergence. Analysis of the results obtained shows the greatest effectiveness of the scheme with a change in half-steps. The results of the calculation are exactly the same as the calculations obtained on a purely implicit scheme and on an explicitly implicit scheme, and the calculation time is much shorter.
numerical method
chemically reacting flows
stiff differential equations
Navier-Stokes equations
implicit method
explicitly implicit method
splitting method

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

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

mol01.wmf, (1)

где τG – характерное газодинамическое время; τchem – характерное время протекания химических реакций; L – характерный размер потока (например, диаметр трубы); V – скорость потока. Как правило, за τchem принимают то время, за которое происходит полное превращение исходных реагентов в продукты реакции.

Наибольшую сложность представляет расчет в случае, когда характерное газодинамическое время сопоставимо со временем протекания химических реакций. Такое течение называется химически неравновесным. Его диапазон, как правило, принимается в пределах примерно: 0,01 < Da < 100. В этом случае необходимо решать уравнения неразрывности всех компонентов газовой смеси и учитывать в них скорости образования этих компонентов в результате всех протекающих химических реакций.

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

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

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

Основная система уравнений

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

mol02.wmf (2)

где ρ – плотность газовой смеси; ui – компоненты скорости; p – давление; τij – тензор вязких напряжений; Cs – массовая доля компонента s; mol03.wmf – скорость образования компонента s в результате химических реакций; NC – количество компонентов газовой смеси; gs, j – диффузионная скорость; Et – термодинамическая часть полной энергия; qj – плотность теплового потока, обусловленного теплопроводностью и диффузией; mol04.wmf – химическая энергия компонента s.

Термодинамическая часть полной энергии E отличается от последней тем, что в нее не входит химическая составляющая:

mol05.wmf. (3)

Численный метод

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

1. Векторная форма записи основной системы уравнений. Система (2) в декартовой системе координат может быть представлена в векторной форме:

mol06.wmf (4)

где

mol07.wmf (5)

mol08.wmf (6)

В новой криволинейной системе координат (ξ, η, ζ) система (4) примет вид

mol09.wmf (7)

где потоки объединены mol10.wmf и:

mol11.wmf (8)

Здесь J – якобиан преобразования координат. Конечно-объемное представление уравнения имеет вид

mol12.wmf (9)

где половинные индексы i + 1/2, i – 1/2, j + 1/2 и т.д. означают потоки на соответствующих гранях контрольных объемов.

2. Неявный метод решения системы в случае отсутствия источников. При отсутствии химических реакций mol13.wmf. Неявные численные методы решения системы (9) описаны во многих работах, в частности в работах [1, 2].

В работе [1] для представления конвективного потока на гранях использовался метод расщепления Стегера –Уорминга [3], а в монографии [2] – метод Рое [4, 5].

При неявном представлении конвективных и вязких потоков на каждом шаге по времени в каждом узле сетки возникает необходимость многократного обращения матриц Якоби mol14.wmf и т.п., а также проведения операций перемножения таких матриц.

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

mol15.wmf, (10)

где матрицы-блоки имеют следующие размеры: A11 – 5х5; A21 – Naх5; D – диагональная матрица размером Na. Здесь Na – число дополнительных уравнений; Na + 5 = Neq – общее число уравнений. Матрица A11 относится к основным уравнениям Навье – Стокса (уравнения неразрывности, количества движения и энергии).

Основные достоинства матрицы вида (10): 1) ее обращение сводится к однократному обращению матрицы A11 и тривиальным операциям перемножения матриц; 2) любые необходимые операции с такой матрицей не меняют ее форму, т.е. блок A12 остается нулевым, а блок A22 остается диагональной матрицей.

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

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

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

Уравнение (9) можно представить в виде

mol16.wmf, (11)

где mol17.wmf – конечно-объемное представление конвекции и диффузии. Здесь оно представлено в явной форме. Источник представлен в неявной форме. Применяем разложение в ряд Тейлора:

mol18.wmf. (12)

Тогда из (11) получается

mol19.wmf, (13)

где mol20.wmf – явное приращение за счет конвекции и диффузии.

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

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

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

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

mol21.wmf (14)

mol22.wmf (15)

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

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

Более перспективным является подход, в котором учитываются оба физических процесса, но происходит чередование их явного и неявного представления. Для этого разбиваем текущий шаг по времени Δt на два полушага Δt/2 и решаем систему в два этапа:

Первый полушаг:

mol23.wmf. (16)

Второй полушаг:

mol24.wmf, (17)

mol25.wmf. (18)

На первом полушаге mol26.wmf представляется в явной форме, а источник mol27.wmf – в неявной. Определяется промежуточное значение вектора mol28.wmf. Затем явное и неявное представления меняются местами, и находится окончательное значение вектора на (n + 1)-ом шаге mol29.wmf.

Можно ввести и более общую форму уравнений на полушагах:

mol30.wmf, (19)

mol31.wmf, (20)

где на втором полушаге: mol32.wmf.

Здесь α – дополнительный параметр. При α = 1 получается чисто неявная схема; значение α > 1 дает завышенную релаксацию, но это может улучшать сходимость метода.

При ϕ = 1, ψ = 1 из (19), (20) получается схема с расщеплением на физические процессы (уравнения (14), (15)), а при ϕ = 0,5, ψ = 0,5 – схема с полушагами (уравнения (16), (17)).

Уравнение (19) решается так же, как в пп. 3:

mol34.wmf. (21)

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

Из (19) следует:

mol36.wmf. (22)

Подставляя это выражение в (20) и используя формулы (18), получаем

mol37.wmf. (23)

Здесь использовалась линейность представления конвекции и диффузии:

mol38.wmf. (24)

Это уравнение решается так же, как в пп. 2.

Результаты расчетов

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

Таблица 1

Параметры на срезе сопла

pa, Па

Ta, K

Ua, м/с

Ra, м

XH2

XH2O

XCO

XCO2

XN2

1х105

1330

2500

1,0

0,4

0,05

0,15

0,05

0,35

Здесь pa, Ta, Ua – соответственно давление, температура и скорость потока на срезе сопла; Xs – мольные доли компонентов; Ra – радиус сопла. Параметры спутного воздушного потока: pe = 0,25х105 Па, Te = 217 K, Ue = 300 м/с.

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

mol1.tif

Распределение температуры вдоль оси струи при использовании различных численных схем: 1 – явно-неявная схема (пп. 2.3); 2 – схема с переменой на полушагах (пп. 2.4 ψ = 0,5, ϕ = 0,5); 3 – схема с расщеплением по физическим процессам (ϕ = 1, ψ = 1)

Таблица 2

Сравнение времени расчетов на сетке 120х60

Схема

t, сек

Схема

t, сек

Чисто неявная

350

С переменой на полушагах (пп. 4)

130

Явно-неявная (пп. 3)

750

С расщеплением по физическим процессам

270

Анализ полученных результатов показывает наибольшую эффективность схемы с переменой на полушагах (пп. 4 ψ = 0,5, ϕ = 0,5); результаты расчета в точности совпадают расчетами, полученными на чисто неявной схеме (на рисунке не приведена) и на явно-неявной схеме, а время расчета существенно меньше. Для схемы с расщеплением вообще не удалось получить результат с удовлетворительной сходимостью.