МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
2019 год, том 31, номер 5, стр. 56-68
ПЕРЕОБУСЛАВЛИВАНИЕ ХИМИЧЕСКИХ ИСТОЧНИКОВ
В УРАВНЕНИЯХ ТИПА
ДИФФУЗИЯ - КОНВЕКЦИЯ - ХИМИЧЕСКАЯ КИНЕТИКА
2019 г.
А.С. Петрусёв
НИЦ «Курчатовский институт»
petrusev_as@nrcki.ru
ГУ «Московский физико-технический институт»
DOI: 10.1134/S0234087919050046
Рассмотрены причины вырождения конечно-разностных уравнений неразрывности
компонентов системы уравнений многокомпонентной гидродинамики химически
активного газа. Предложен метод переобуславливания, позволяющий преодолеть
подобное вырождение уравнений. Метод отработан при одномерном моделирова-
нии горения газовых смесей и оказался работоспособным в случае наличия ано-
мально быстрых химических реакций.
Ключевые слова: химическая кинетика, вырождение уравнений, переобуславлива-
ние, квазилинейное исключение.
THE PRECONDITIONING OF CHEMICAL SOURCES TERMS IN EQUATIONS
WITH DIFFUSION, CONVECTION AND CHEMICAL KINETICS
A.S. Petrusev
NRC «Kurchatov Institute»
SU «Moscow Institute of Physics and Technology»
The reasons of degeneration have been considered for finite-differences equations of
multi-components hydrodynamics of chemically active gas. To prevent the equations de-
generation, the preconditioning method is developed. The method is adjusted for model-
ing of one-dimension gas mixture combustion in the case of extremely fast chemical re-
actions present.
Key words: chemical kinetics, equations degradation, preconditioning, quasi-linear
elimination.
1. Введение
Численное интегрирование жёстких уравнений химической кинетики,
а также включающих их более сложных уравнений многокомпонентной
Переобуславливание химических источников в уравнениях типа …
57
гидродинамики, требует использования неявных численных методов. Один
из распространённых способов решения систем уравнений в частных про-
изводных состоит в том, что после дискретизации по пространству они сво-
дятся к системам обыкновенных дифференциальных уравнений (ОДУ) вы-
сокого порядка с помощью метода прямых. Для решения полученной сис-
темы ОДУ используются широко распространённые неявные одношаговые
или многошаговые методы. Так устроены, в частности, многие современ-
ные программные пакеты [1,2] для моделирования одномерного горения.
Подавляющее большинство неявных численных методов высокого по-
рядка точности, используемых для интегрирования ОДУ, основаны в той
или иной форме на решении систем линейных алгебраических уравнений
(СЛАУ). Так, многошаговые методы Гира, Адамса или им подобные, а так-
же неявные одношаговые методы Рунге-Кутты основаны на решении сис-
темы нелинейных уравнений методом Ньютона, использующим решение
СЛАУ. Полуявные одношаговые методы Розенброка не используют метод
Ньютона, но непосредственно включают решение СЛАУ. Таким образом, и
те, и другие требуют решения СЛАУ, для чего обычно используется метод
Гаусса.
Таким образом, работоспособность указанных неявных методов суще-
ственно зависит от разрешимости системы уравнений вида (упрощённо):
(1 J)u f ,
(1)
где f - вектор правой части ОДУ, J Df / Du - матрица Якоби правой час-
ти системы ОДУ, - временной шаг. Причём для нормальной работы неяв-
ного алгоритма требуется разрешимость (1) при больших , вплоть до
. А для этого желательна невырожденность матрицы Якоби J.
Поясним сказанное на примере. Для системы ОДУ типа
x
(
x y
),
(2)
(
x y
),
где <<-1 - параметр, матрица Якоби всегда вырождена, что, казалось бы,
не препятствует численному решению системы. Как известно, после прохо-
да погранслоя, на квазистационарном участке, интегрирование жёстких
ОДУ неявными методами ведётся с большим временным шагом
J
1.
При этом условие
1J
0 остаётся существенным. Но для системы (2)
вычисление определителя в условиях машинной арифметики сопровожда-
58
А.С. Петрусёв
ется потерей точности, т.к. содержит разность больших близких чисел. Ес-
ли обозначить через точность представления машинного числа (примерно
108 для чисел одинарной точности, 1016 для двойной и т.д.), то
2 2
1
12
12
,
2 2
2 2
1
J
12
2 2
1
0
12
2 2
В случае когда
12
, указанный определитель будет равен
нулю, т.е. матрица (1 J) будет восприниматься компьютером как вырож-
денная. Это обстоятельство накладывает ограничение на временной шаг:
1
(
)
(3)
В случае больших
такое ограничение может быть весьма обремени-
тельным. Причём даже в случае выполнения (3) решение системы ОДУ мо-
1
жет иметь значительную погрешность, а численная матрица
(1
J
)
, по
сравнению с точной, совершенно другой спектр и утратить свои стабилизи-
рующие свойства для слабо жёстких мод. На практике оба эти обстоятель-
ства ведут к нарушению нормальной работы неявных алгоритмов, вызывая
уменьшение временного шага и резкое замедление счёта.
В то же время, введя замену переменных
(u,v) ( y x, y x) ,
(4)
получим эквивалентную систему уравнений
u
2
u,
(5)
0.
Для (5)
1J
12 0, причём ограничение типа (3) отсутствует.
Ограничения, аналогичные (3), возможны и в случае более сложных систем
ОДУ, а угадать нужную подстановку (4) становится затруднительно.
Численное вырождение матрицы Якоби сопровождается ростом её
числа обусловленности. Однако число обусловленности (обычно опреде-
1
ляемое как
M(J)
J
J
и много обсуждаемое в литературе) ничего не
говорит о возможности численного решения (1). Например, при 0 1,
для эквивалентных СЛАУ (6), (7) и (8)
x y
0,
(6)
x
(1)
y
,
Переобуславливание химических источников в уравнениях типа …
59
x y
0,
(7)
y
,
x
1,
(8)
y
числа обусловленности велики ( 1/ ) и близки друг другу. Тем не менее
20
при
10
системы (7) и (8) решаются численно методом Гаусса без ви-
димых трудностей, в то время как попытка решения системы (6) приводит к
делению на ноль.
Разрешимость системы уравнений можно характеризовать «числом
Адамара» матрицы, характеризующим ортогональность её векторов-строк:
1
n
n
2
H(J)
det(J)
J
1
ij
i
1
j
1
Обычно СЛАУ с H ~ 1 (H ) разрешима численно без сколько-
нибудь заметного сокращения знаков, а численное решение системы с
H 1 (H ) приводит к большой погрешности результата или делению на
ноль. Для систем (6)-(8) числа Адамара примерно равны , 0.7 и 1 соответ-
ственно. Использованная выше замена переменных (4) приводит к ортого-
нализации векторов-строк матрицы Якоби. Для (5) соответствующее число
Адамара равно 1.
Причина такого поведения численных алгоритмов вполне понятна: при
использовании машинных чисел с плавающей точкой сокращение знаков
происходит при вычитании близких чисел, а при делении - нет. При деле-
нии на малую величину в СЛАУ (7) и (8) сокращения знаков не происхо-
дит, но происходит при вычислении определителя СЛАУ (6) из-за того, что
векторы-строки матриц этих систем уравнений почти коллинеарны, а это
ведёт к вычислению разности близких чисел.
Ограничение H можно преодолеть, перейдя к использованию длин-
ных машинных чисел высокой точности. Существуют математические па-
кеты, предоставляющие подпрограммы альтернативной арифметики, опе-
рирующие с вещественными числами длиной 16 байт и более. Однако автор
полагает, что это не тот способ, которым следует решать проблему.
Предлагаемый ниже метод предназначен для уравнений, включающих
химическую кинетику с быстрыми реакциями. Отличительной особенно-
стью этого случая является резко выраженное вырождение уравнений (час-
60
А.С. Петрусёв
то ещё на этапе вычисления их коэффициентов), что вызывает затруднения
при их решении. В то же время, используя вид химического источникового
члена, удаётся построить эффективный алгоритм переобуславливания.
2. Переобуславливание для уравнений химической кинетики
Основные особенности задачи проявляются уже при решении 0D урав-
нений (химическая кинетика без конвекции и молекулярного транспорта):
y
W,
(9)
t
где y и W - векторы концентраций химических компонентов и химических
источников соответственно размерности m, m - число веществ в системе.
Трудности могут возникать при приближении к равновесному составу, ко-
гда быстрые прямые и обратные реакции почти компенсируют друг друга, а
их малая разность примерно равна вкладу медленных реакций. При этом (9)
переходит в систему уравнений химического равновесия, которые подвер-
жены численному вырождению, известному как проблема «малых концен-
траций» [3]. Вычислительные трудности подобного рода также могут при-
сутствовать в системе (9) и вдали от равновесия при наличии в используе-
мом кинетическом механизме химических реакций с аномально высокой
скоростью. Аномально быстрые реакции могут встречаться, например, при
построении кинетического механизма методом обратной подгонки, когда
константы скоростей реакций подбираются по поведению решения прямой
кинетической задачи. В этом случае константы скоростей не имеют непо-
средственного физического смысла и могут превосходить газокинетиче-
ские. Эти обстоятельства ведут к вырождению получаемой СЛАУ типа (1),
что накладывает довольно жёсткие требования на численный алгоритм ре-
шения прямой кинетической задачи.
Предположим, что имеется реакция с такой большой скоростью
w
,
j
1
что для всех остальных реакций
w
w
, при
j j
. Тогда численно
w
j
j
1
j
1
1
w
w
и Wi пропорциональны
w
для всех веществ i, входящих в j1-ю
j
j
1
j1
реакцию. Такое преобладание
w
приводит как к вырождению матрицы
1
j
Якоби (из-за пропорциональности нескольких строк), так и к неверному
вычислению правых частей (медленные реакции маскируются j1-й, их ско-
рости не попадают в разрядную сетку машинного числа). В этом случае
система (9) будет подобна (2), которая, как отмечалось, неудобна для чис-
ленного решения. Данная ситуация вполне аналогична упоминавшейся про-
Переобуславливание химических источников в уравнениях типа …
61
блеме «малых концентраций», возникающей при расчёте термодинамиче-
ского равновесия [3]. Для преодоления таких трудностей в [4] был предло-
жен алгоритм квазилинейного исключения, позволяющий эффективно пре-
одолеть численное вырождение системы уравнений элементного баланса.
Аналогичный метод может быть применён и для системы (9). Уравнения (9)
нелинейны относительно y, однако они линейны поwj ; химические источ-
никовые члены W разлагаются в сумму по скоростям отдельных реакций:
W
Sw. Здесь
=
W
W
- вектор химических источников размерности m,
i
=[
]
w
w
- вектор скоростей реакций размерности n, (n - число реакций),
j
[
]
S
s
- стехиометрическая матрица размерности m n . Все реакции счи-
ij
таются прямыми (необратимыми); обратимые реакции представляются в
виде двух необратимых. При этом всегда wj>0.
Аналогично случаю уравнений элементного баланса [4], благодаря
квазилинейной структуре W, можно перейти от (9) к преобразованной сис-
теме по записи, аналогичной (8), применяя квазилинейное исключение ли-
дирующих элементов. Целью такого исключения является присутствие
скорости
w
быстрой j1-й реакции только в одном из уравнений (9). Про-
1
j
цесс квазилинейного исключения является некоторым аналогом исключе-
ния Гаусса с выбором главного элемента. В данном случае лидирующие
элементы следует выбирать по величине модулей произведений sijwj. Как и
в случае уравнений элементного баланса, квазилинейное исключение не
позволяет привести систему (9) к диагональному виду, аналогичному (8),
поскольку число реакций обычно больше числа веществ. Однако это имеет
место по отношению к первым m n , наиболее быстрым реакциям. Как
альтернатива, (9) можно приводить к форме, аналогичной треугольной (7).
Использование такой псевдотреугольной формы не менее эффективно для
численного решения, зато более экономично по числу арифметических
операций. Везде ниже используется именно этот вариант исключения.
Квазилинейное исключение предлагается производить следующим об-
разом:
Начинаем с первой строки, i=1.
В i-й строке матрицы S выбирается j-й «лидирующий элемент» sij, для
которого произведение sijwj максимально по модулю.
i-я строка последовательно прибавляется ко всем нижеследующим
строкам, k i 1,m , с таким коэффициентом, чтобы j-й элемент k-й
строки обратился в нуль.
Переходим к (i+1)-й строке.
62
А.С. Петрусёв
Таким образом, при квазилинейном исключении преобразуется только
стехиометрическая матрица S : S S . После процедуры исключения все
элементы skj равны нулю при k>i, где j - номер лидирующего элемента i-й
строки матрицы S. Благодаря тому, что все стехиометрические коэффици-
енты по порядку величины равны единице, квазилинейное исключение не
сопровождается сколько-нибудь существенной потерей точности (сокраще-
нием знаков). Принципиально важно, чтобы на месте коэффициентов skj,
обращаемых в нуль при исключении, стояли точные нули, а не малые вели-
чины, связанные с приближённостью машинной арифметики.
Алгоритм квазилинейного исключения требует довольно подробной
информации о структуре источникового члена W, что затрудняет числен-
ное нахождение преобразованной матрицы Якоби правой части уравнений
(9), её необходимо вычислять аналитически.
Алгоритм остаётся применимым при записи уравнений (9) относи-
тельно массовых либо смешанных долей химических компонентов. Ис-
пользование смешанных долей по сравнению с массовыми удобно тем, что
стехиометрические коэффициенты sij оказываются целыми и безразмерны-
ми, а в коэффициенты квазилинейного исключения не входит отношение
молекулярных весов.
Описанное квазилинейное исключение можно представить в векторно-
матричном виде S QS , где Q - матрица размером m m . Она имеет тре-
угольный вид: в ней заполнен нижний треугольник, а по диагонали стоят
единицы. В верхнем треугольнике все элементы нулевые. Аналогично про-
изводится и преобразование матрицы Якоби правой части.
Скорости реакций в стандартном виде вычисляются как:
m a
jl
w
K
(T)
y
,
j
j
l
1
l
где Kj(T) - константа скорости j-й реакции, ajl - эмпирические коэффициен-
ты, для простых реакций совпадающие со стехиометрическими slj. Для
сложных реакций коэффициенты ajl могут быть нецелыми и не совпадать с
slj. Отсюда
n
w
n
a
n
W
j
jl
i
J
s
s
w
s
r
;
i,l
1,m,
il
ij
ij
j
ij jl
y
y
y
l
j1
l
j1
l
j1
или J=SR, где R=[rjl], rjl =w
a
/
y
- прямоугольная матрица Якоби скорос-
j jl
l
Переобуславливание химических источников в уравнениях типа …
63
тей реакций размером nm. После квазилинейного исключения J=SR=QSR.
При численном решении (9) от системы уравнений типа
(1 J)ξ b(Sw) ,
(10)
переходим к системе уравнений
(Q J)ξ b(S )
w .
(11)
Можно показать, что описанное выше вырождение системы уравнений
(10), вызванное присутствием быстрой реакции, в преобразованной системе
уравнений (11) отсутствует. В самом деле, большие величины
w
и
r
бу-
j
j
k
1
1
дут преобладающими только в одном уравнении системы (11). В силу по-
строения алгоритма квазилинейного исключения, если скорость
w
доми-
j
1
нирует в i1-м уравнении, то в k-м уравнении она будет отсутствовать при
k>i1. При k<i1 скорость
w
либо отсутствует, либо присутствует, но не до-
j1
минирует, т.к. не была выбрана в k-м уравнении лидирующим элементом
(иначе
w
отсутствовала бы в i1-м уравнении).
j
1
3. Переобуславливание для одномерных конечно-разностных уравне-
ний химическая кинетика - конвекция - диффузия
При моделировании широкого ряда физико-химических процессов
возникает необходимость численного решения систем уравнений типа кон-
векция - диффузия - химия:
y
i
v
y
divj
W
,
i
1,m,
(12)
i
i
i
t
где y - вектор смешанных долей химических компонентов, - плотность
смеси, ji - мольный диффузионный поток i-го компонента, v - среднемас-
совая скорость смеси. Использование смешанных долей вместо массовых
связано с причинами, отмеченными выше. В частности, такие уравнения
встречаются при численном моделировании процессов горения газовых
смесей. На рассмотрении этого примера в одномерном случае мы и остано-
вимся.
Трудность численного решения (12) существенно зависит от числа
Дамкёлера, характеризующего отношение величин химических источников
и диффузионных потоков. Наиболее прост для численного решения случай
малых и умеренных чисел Дамкёлера, когда величины диффузионных по-
токов относительно велики. Это соответствует большой скорости горения.
В случаях медленного горения, например, когда состав исходной смеси да-
64
А.С. Петрусёв
лёк от стехиометрического, в уравнении (12) преобладает член Wi. Как пра-
вило, такое преобладание связано с тем, что локальный состав смеси близок
к равновесному, прямые и обратные химические реакции почти полностью
компенсируют друг друга, а их малая разность примерно равна вкладу кон-
вективных и диффузионных процессов. Этот случай известен как «химиче-
ски равновесные течения» [5]. В случае таких течений, в [5] отмечается по-
теря точности при решении уравнений, содержащих химические источни-
ковые члены. Там же для моделирования течений подобного типа рекомен-
дуется либо полуявный учёт этих членов с записью уравнений относитель-
но приращений величин концентраций на шаге интегрирования, либо ис-
пользование специального алгоритма с заменой уравнений неразрывности
компонентов на уравнения химического равновесия Гульдберга-Вааге.
Первый способ действительно улучшает обусловленность решаемых урав-
нений, однако не устраняет сокращение знаков при вычислении их правых
частей. Кроме того, для около равновесных течений он обычно приводит к
чрезвычайно медленной сходимости алгоритма. Очевидным недостатком
второго подхода является грубость описания диффузионных процессов, что
ведёт к заметному дисбалансу в уравнениях неразрывности компонентов.
Помимо этого, как отмечалось выше, сами уравнения химического равно-
весия также подвержены вырождению.
Большие числа Дамкёлера могут встречаться и в случае сравнительно
быстрых пламён, если присутствуют аномально быстрые реакции. В обоих
случаях может происходить вырождение уравнений, полностью аналогич-
ное описанному выше для (9).
В одномерном случае при использовании трёхточечной аппроксима-
ции дискретный аналог (12) можно записать как
A
C
J
ξ
A
ξ
C
ξ
b
Sw
,
(13)
k
k
k k
k k
1
k k
1
k
k
где =1/, k - номер узла сетки,
ξ
y
y
- приращения неизвестных на
k
k
k
временном шаге. Здесь
ξ
, b
- векторы размера m;
A
, C
,
J
- матрицы
k k
k
k k
mm. Матрицы
A
, C
и векторbk соответствуют конвективным и диффу-
k
k
зионным членам;Jk , как и ранее, матрица Якоби химических источников.
Как и ранее, система уравнений (13) линейна поξk , за исключением
правой части, которая линейна по химическим скоростямwk , что позволя-
ет применить алгоритм квазилинейного исключения в полной аналогии со
случаем системы (10). Преобразование уравнений сводится к их почленно-
му домножению в каждом узле сетки на переобуславливающую матрицу
Переобуславливание химических источников в уравнениях типа …
65
Qk, после чего уравнения принимают вид
A
C
Q
J
ξ
A
ξ
C
ξ
Q
b
S
w ,
(14)
k
k
k
k k
k k
1
k k
1
k k
k k
A
Q
A
;
C
Q
C
;
J
Q
J
Q
SR
;
S
Q
S.
k
k k
k
k k
k
k k
k
k
k
k
МатрицаQk строится в каждом узле сетки аналогично описанному ра-
нее, но наличие конвекции и диффузии препятствует вырождению уравне-
ний. Вырождение возможно только если
s
w
b
. В обратном случае
ij
jk
ik
квазилинейное исключение в i-м уравнении не требуется. Данное обстоя-
тельство несколько сокращает потребное число арифметических операций.
Вычисление преобразованных химических источниковых членов в
правой части (14) следует производить как (Q
S)w
, а не
Q
(Sw
). Хотя
k
k
k
k
второй вариант требует меньше арифметических операций, но при этом те-
ряется эффект переобуславливания. Аналогично, k
J следует вычислять как
(Q
S)
R
, а не
Q
(SR
). Поэтому алгоритм требует вычисления и хране-
k
k
k
k
ния в памяти громоздкой матрицыS
Q
S в каждом узле сетки. Эти об-
k
стоятельства приводят к дополнительным затратам как оперативной памя-
ти, так и арифметических операций.
Конкретный вид
A
,
C
,b зависит как от используемой модели мно-
k
k k
гокомпонентной диффузии, так и от её численной реализации. Эти матрицы
могут быть диагональными или плотно заполненными, вычисляться ана-
литически или численно (как это делается в [1, 2]). Алгоритм квазилинейно-
го исключения остаётся применимым во всех этих случаях, но всегда тре-
бует аналитического вычисления матрицы Якоби химических источников.
4. Уравнение энергобаланса
Если в задаче требуется рассчитывать поле температуры, то возникает
необходимость решения уравнения энергобаланса. Одномерное уравнение
в недивергентной форме имеет вид
m
m
T
T
T
c
c
v
c
j
H
W
0
(15)
p
p
pi i
i i
t
z
z
z
i1
i1
Здесь cp - теплоёмкость смеси, T - её температура, v - среднемассовая ско-
рость смеси, - коэффициент теплопроводности смеси,
c
c
(T),
pi
pi
H
H
(T)
- молярные теплоёмкость и энтальпия i-го компонента смеси
i
i
соответственно.
66
А.С. Петрусёв
Наличие членов Wi приводит к потере точности, аналогично уравне-
нию (12), и дискретный аналог (15) требует переобуславливания. Однако
присутствие множителей Hi в сумме нарушает подобие слагаемых - хими-
ческих источников, имеющееся в (12), что затрудняет процедуру квазили-
нейного исключения. Более удобен другой подход. При записи уравнения
энергобаланса в дивергентной или полудивергентной форме слагаемые ви-
W отсутствуют:
m
H
H
T
v
H
j
0
i i
t
z
z
z
z
i1
или
m
m
T
T
T
y
y
j
i
i
i
c
c
v
c
j
H
v
0,
(16)
p
p
pi i
i
t
z
z
z
t
z
z
i1
i
1
где H H ( y,T ) - энтальпия смеси.
Уравнение (16) не содержит явно химических источников, поэтому его
конечно-разностный аналог не требует переобуславливания, независимо
как от используемых моделей диффузии и теплопроводности, так и от их
численной реализации.
5. Численный эксперимент
Для проверки и отработки алгоритма было проведено несколько проб-
ных расчётов скорости горения смеси ацетилен - воздух с использованием
кинетической схемы высокой жёсткости из 42 веществ и 256 реакций, сре-
ди которых имелись аномально быстрые.
При расчётах для сравнения использовались три программы, отли-
чающиеся способом решения стационарной системы конечно-разностных
уравнений. Все программы использовали процедуру итерационного реше-
ния на основе одноэтапного метода Розенброка. Первая и третья програм-
мы использовали метод Розенброка с вещественными коэффициентами. Во
второй использовался метод Розенброка с комплексными коэффициентами
и улучшенной устойчивостью. Третья программа включала также описан-
ную выше процедуру квазилинейного исключения. Во всех случаях расчёт
вёлся на установление, что позволило использовать упрощённую процеду-
ру автоматического выбора временного шага (основанную на ограничении
максимального изменения решения на шаге), а также вычислять, переобу-
славливать (только в третьей программе) и обращать матрицу Якоби один
Переобуславливание химических источников в уравнениях типа …
67
раз на 10 итераций. При вычислении матрицы Якоби химические источни-
ки учитывались точно, а диффузионные - приближённо, что существенно
упростило их вычисление. В используемой модели термо- и бародиффузия
не учитывались, а для вычисления диффузионных потоков химических ком-
понентов применялся полуявный алгоритм решения уравнений Стефана-
Максвелла [5].
Расчёты проводились на персональном компьютере с одноядерным
процессором Intel тактовой частотой 3.1 ГГц. Использовались числа двой-
ной точности (64-битовые) и неравномерная адаптивная сетка, содержащая
100 узлов. На начальном этапе расчёта во всех трёх программах производи-
лось ограничение констант скоростей реакций сверху, что позволило «от-
ключить» аномально быстрые реакции и построить хорошее начальное
приближение. При этом переобуславливание в третьей программе не про-
водилось. На окончательном этапе константы скорости реакций не ограни-
чивались, а в третьей программе выполнялось переобуславливание. И на
первом, и на втором этапах счёт начинался с мелким шагом 1010с.
Как и следовало ожидать, на первом этапе работа всех трёх программ
протекала практически одинаково, вычисленные начальные приближения
совпали с хорошей точностью. На втором этапе первая и вторая программы
после 500 итераций выбирали чрезвычайно мелкий шаг по времени (1015
1014с), при котором решение практически не изменялось (и погрешность
не убывала) на протяжении 2500 итераций, после чего программы были ос-
тановлены. Третья программа без видимых затруднений наращивала шаг до
значения около 0.01с. Решение сошлось с заданной относительной точно-
стью (103) за 170 итераций, что заняло около 20с.
Отметим, что распространённые программные пакеты для моделиро-
вания одномерного горения [1, 2] в случае данной кинетической схемы
также оказались неработоспособными: вычисления сопровождались зацик-
ливанием или авостом.
6. Выводы
Описанная процедура квазилинейного исключения позволила постро-
ить алгоритм решения уравнений (12), (16), устойчивый к численному вы-
рождению указанного выше типа. В отличие от упоминавшегося метода
расчёта химически равновесных течений, преобразованные уравнения (14)
полностью натуральны (физически эквивалентны исходным (12)) и не при-
водят к дополнительным погрешностям в описании процессов как диффу-
зии или конвекции, так и химической кинетики.
68
А.С. Петрусёв
С использованием этого алгоритма был построен блок моделирования
одномерного горения, использованный в программном комплексе «Хими-
ческий Верстак» [6], многократно тестировавшийся на модельных расчётах
различных ламинарных газовых пламён и сопоставлением результатов рас-
чётов с пакетами [1, 2].
Вместе с тем, алгоритму присущи определённые недостатки. Преодоле-
ние вырождения достигается довольно высокой ценой. Прежде всего, это
дополнительные трудности при программировании вычисления химичес-
ких источников, их матрицы Якоби и коэффициентов уравнений (14). Кро-
ме того, для решения (14) и (16) необходимо использовать векторную про-
гонку, обладающую сравнительно большой трудоёмкостью. Экономичный
полуявный алгоритм решения уравнений неразрывности компонентов [5] и
энергобаланса [7], реализуемый раздельными скалярными прогонками, в
данном случае неприменим, т.к. ни (14), ни дискретные аналоги (16) не
имеют необходимого диагонального преобладания.
СПИСОК ЛИТЕРАТУРЫ
1.
2.
corporate/resourcelibrary/brochure/chemkin-pro-brochure.pdf
3.
Г.В. Белов. Термодинамическое моделирование. - М.: Научный Мир, 2000.
G.V. Belov. Termodinamicheskoe modelirovanie. - M.: Nauchniy Mir, 2000.
4.
А.С. Петрусёв. О проблеме «малых концентраций» при расчёте термодинамиче-
ского равновесия // Теплофизика высоких температур, 2005, т.43, №4, с.526-532.
англ. пер.: A.S. Petrusev. The Problem of "Low Concentrations" in the Calculation of
Thermodynamic Equilibrium // High Temperature, 2005, v.43, № 4, p.522-529.
5.
Ю.В. Лапин, М.Х. Стрелец. Внутренние течения газовых смесей. - М.: Наука, 1989.
Iu.V. Lapin, M.H. Streles. Vnutrennie techeniya gazovyh smesey. - M.: Nauka, 1989.
6.
7.
А.С. Петрусёв Численное моделирование ламинарного факела / Материалы 51-й
конференции МФТИ, 2008, с.73-76.
A.S. Petrusev. Chislennoe modelirovanie laminarnogo fakela / Materialy 51 Konferentsii
MFTI, 2008, s.73-76.
Поступила в редакцию 04.06.18
После доработки 22.11.18
Принята к публикации 10.12.18
