Математическое моделирование химически реагирующих течений Текст научной статьи по специальности «Математика»

УДК 536.4:621.454

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ХИМИЧЕСКИ РЕАГИРУЮЩИХ ТЕЧЕНИЙ

КОРЕПАНОВ М.А.

Институт прикладной механики УрО РАН, Ижевск, Россия, ipm@udman.ru

АННОТАЦИЯ. Описана система уравнений и метод ее решения, позволяющие определять состав химически реагирующей газовой смеси с учетом скорости химических реакций и максимально возможного числа веществ. Приведены результаты численных исследований течений в соплах ракетных двигателей.

КЛЮЧЕВЫЕ СЛОВА: математическое моделирование, метод решения, химический состав, химическая реакция

Список обозначений:

I - энтальпия, кДж/кг;

Кр - константа равновесия химической реакции;

К - константа скорости обратной химической реакции;

К - число "кинетически" учитываемых веществ;

I - число молекулярных компонентов газовой смеси; т - число атомарных компонентов газовой смеси; п - общее рассматриваемых химических реакций; р - давление, МПа;

Я - газовая постоянная, Дж/(кг К), Я0 =8,31441 Дж/(моль К);

Т - температура, К;

V - объем, м3;

и - внутренняя энергия, кДж/кг;

- количество атомов ьго химического элемента в j-ом веществе;

ЬгТ - количество атомов ьго химического элемента в топливе; q - общее количество компонентов смеси;

а - коэффициент избытка окислителя; д - мольная масса, кг/моль; р - плотность, кг/м3;

Уу, Уу ’ - стехиометрические коэффициенты химических реакций; г, ], к, s - индексы.

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

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

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

Основным недостатком большинства моделей является отсутствие данных

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

и равновесный, основанный на концепции химического равновесия [4, 5]. Первый сложен с точки зрения описания системы химических реакций и дальнейшего численного решения, а второй не учитывает время, необходимое для протекания химических реакций.

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

Таким образом, процессы в НРИС можно разделить на три группы:

1) изменение состава и свойств рабочего тела за счет протекания химических реакций;

2) изменение состава и количества вещества в реакторе за счет массообмена с соседними реакторами или с окружающей средой:

- в форме расходов (из-за перепада давлений);

- в виде диффузии (из-за градиента концентраций);

3) теплообмен:

- изменение внутренней энергии за счет массообмена;

- изменение внутренней энергии за счет теплообмена со стенками реактора;

- нестационарные тепловые процессы в стенках реактора.

V = ад

с

Рис.1. Нестационарный реактор идеального смешения

ИЗМЕНЕНИЕ СОСТАВА И СВОЙСТВ РАБОЧЕГО ТЕЛА ЗА СЧЕТ ХИМИЧЕСКИХ РЕАКЦИЙ

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

Другим широко распространенным методом является описание процессов в реагирующих смесях с точки зрения формальной химической кинетики [1, 2, 3, 6]. В этом случае для нахождения состава реагирующей смеси в любой момент времени необходимо решать систему дифференциальных уравнений, обладающую "жесткостью" [6, 7]. При этом во внимание принимаются только те вещества, для которых есть данные по их термодинамическим свойствам и по протекающим с их участием химическим реакциям. Количество веществ, принимаемых в рассмотрение, в этом случае обычно существенно меньше, чем при равновесной модели.

Для устранения указанных сложностей можно использовать концепцию "стесненного равновесия", рассмотренную в работе [8]. На основании этой концепции было предложено ввести в метод расчета равновесного состава реагирующей смеси вместо уравнений диссоциации уравнения химической кинетики для определяющих веществ [9]. Это позволяет наложить на систему уравнений ограничения, определяемые конечным временем протекания химических реакций, при сохранении максимально возможного набора учитываемых веществ.

Для удобства численного решения все уравнения должны быть записаны в одинаковых переменных. В [6] отмечается, что состав химически реагирующей смеси целесообразно характеризовать относительными величинами: мольными г. или весовыми долями gi. При изменении параметров состояния р, V, Т эти величины не изменяются, в отличие от используемых в уравнениях равновесной термодинамики парциальных давлений р. и в уравнениях химической кинетики мольных концентраций С..

В [4, 6] предлагается при расчетах использовать в качестве неизвестных логарифмы искомых величин (мольных долей). Это позволяет, во-первых, линеаризовать уравнения закона действующих масс для равновесной термодинамики, а во-вторых - увеличивает устойчивость расчетов, т.к. исключает появление в процессе расчета отрицательных значений искомых величин, не имеющих физического смысла.

Учитывая то, что г. = pi / р , уравнения для определения равновесного состава газовой

смеси [4, 5] подстановкой 1пр. = 1п(р • г. ) = 1пр + 1п г. могут быть преобразованы

к следующему виду:

- уравнение диссоциации:

- уравнение сохранения вещества:

( ч Л

1п ^а. • г. - 1пМТ - 1пЬТ + 1пр = 0, I = 1,2,...,т ;

V .= )

- уравнение закона Дальтона заменяется уравнением нормировки:

(2)

(3)

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

drl

k

dt

ч

| £сР Л -1

V *оТ У

_ ( ч \ ( ч ^

Кр} • рм+-м • ехр| Е у& \щ - ехр| Е у&'\пг

У

У.

£ = 1,2............К,

(4)

где М + =ЕУу - число молекул, участвующих в прямой реакции;

,=1

ч

М- = Еуу' - число молекул, участвующих в обратной реакции;

,=1

п - общее число рассматриваемых химических реакций.

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

ч ч

у • М =ЕУу'Мг ,

2=1 2=1

где Mi - молекула /-го вещества;

у - номер химической реакции.

Решение системы уравнений (1 - 4) целесообразно вести относительно логарифмов мольных долей 1пг/. Для удобства записи уравнений необходимо ввести переменную

= 1пг/, дифференциал которой равен dyi = (1/ г,) • dri. Уравнения (1 - 4) после замены переменной будут выглядеть следующим образом:

- уравнение диссоциации

т ( т \

У]- Е ау • у, +1пКр + 1 - Е ау •1п^ =0, у = 1,2,.,/ - К;

,=1 V ,=1 У

(5)

уравнение сохранения вещества

г

1п

Е ау • ехр(^у) - 1пМт - 1пЬ,т + 1п^ = 0, / = 1,2,. ..,т ;

(6)

- уравнение нормировки

1п

Л

Е ехр(г,)

V2=1

=0;

(7)

уравнение химической кинетики

d У £ dt

( ч

= Е (У У-У ^ • К 1(Т) • еХ^ ЕУ £, '•у,-^,

У=1

V ,=1

^ ( Кр_^

V *оТ У

м--1

ехр

чч

/пК] +(м + -м-)• /пр + ЕVу • у, - Е'VуЧ

,=1

У у

,=1 У

-1

£ = 1,2 К

(8)

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

Изменение состава и количества вещества, а также изменение внутренней энергии в НРИС описывается системой обыкновенных дифференциальных уравнений [5, 6].

ч

п

X

X

Для замыкания системы уравнений, описывающих процессы в НРИС, ее необходимо дополнить уравнением состояния

p = р • R • T (9)

и зависимостью внутренней энергии газовой смеси от температуры

U = f(T, /ь ^.^ rq). (10)

Полученная система уравнений, описывающая процессы в НРИС, состоит из

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

РЕШЕНИЕ СИСТЕМЫ УРАВНЕНИЙ ДЛЯ ОПРЕДЕЛЕНИЯ СОСТАВА

И ТЕМПЕРАТУРЫ РЕАГИРУЮЩЕЙ ГАЗОВОЙ СМЕСИ

Система уравнений для определения состава высокотемпературной реагирующей газовой смеси при p, T = const состоит из уравнений (5 - 8). При этом уравнение кинетики (8) является дифференциальным уравнением и описывает неравновесный релаксационный процесс. В [1] отмечается, что подобные уравнения вблизи равновесия являются, как правило, уравнениями с малым параметром при старшей производной, что существенно

усложняет их численное интегрирование. В [1, 7] отмечается также, что для кинетических

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

В [1, 6, 7] отмечается, что для численного интегрирования систем уравнений химической кинетики лучше всего подходят методы, основанные на использовании неявных разностных схем. В [1, 6] для решения уравнения вида:

d^7 = f (x,yi,y2,...,yq ) (11)

at

предлагается использовать следующую конечно-разностную схему :

Уп+1 -Уп -h• [аf(xn,yn) + (1 -a>f(xn+1 ,yn+1 )]= ^ (12)

где n - номер шага по времени; h - шаг по времени;

а - параметр неявной схемы, в [1, 6] рекомендуется а = 0,4 .

При записи кинетических уравнений (8) в виде (12) для решения системы уравнений определения состава реагирующей газовой смеси можно использовать метод Ньютона, при использовании которого исходная система нелинейных уравнений относительно yi сводится

к системе линейных уравнений относительно поправок A/i.

Для увеличения эффективности метода Ньютона в [6, 7] рекомендуется аналитическое вычисление частных производных. Линеаризованные относительно Ayi уравнения диссоциации (5), сохранения вещества (6) и нормировки (7) примут следующий вид: уравнение диссоциации:

m

AYj -Z -ДУг =-5j , j = 1,2,..У - K , (13)

i =1

mm

где = yj aj Yi + XnKj + 1 aij • Ыр;

i =1

уравнение сохранения вещества:

9

Ё аь- •Ау]. • ехр() - В - А 1пМТ = -5; • Д, г = 1,2,. ..,да , (14)

1=1

9

гДе Д = Ё аа -ехР(^;),

1=1

5І = 1п Д. - 1пМТ - 1пЬгТ + 1п^; уравнение нормировки:

9 9

2 еХР(^ ) • АУг=-5г- Ё еХР(^ ), (15)

г=1 г =1

( 9

где 5Г = 1п

Ё ехР(^) I.

V г=1 )

Уравнения химической кинетики в линеаризованной форме записываются в виде

9 гдО

Ё -Аг, = -вкп , = 1,2,...Д, (16)

к=1 )

где 5Кг = - у - к • [а • /(х„, у) +(1 - а) • /(х„+1, ^)].

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

записанных в форме (12), имеют вид: д/<

п /9 ^

X

дУ,

= 5к -1 -а)• Ё(V„ -Vj(T)•exp Ёv/• У* -У

1 =1

V*=1 )

/ \М--1

г кр_ ^

V *оТ )

- (V1’-5 ,Н, г = 1,2,...,К; к = 1,2................9, (17)

V г

+Ёv * •г * -Ёv */ • у *

5=1 1=1

здесь 5„ - символ Кронекера,

, Г1, г = к,

5 = Г [о, г ф к.

Для удобства решения системы линеаризованных уравнений (13 - 16) матрица коэффициентов при поправках (рис. 2) структурирована по аналогии с [4, 5]:

столбцы:

1) поправки, относящиеся к молекулярным продуктам, учитываемым по уравнениям диссоциации;

2) поправки, относящиеся к атомарным продуктам;

3) поправки, относящиеся к молекулярным продуктам, учитываемым по уравнениям химической кинетики;

4) поправка к числу МТ ;

строки:

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

2) уравнения сохранения вещества в той же последовательности, что и атомарные продукты в столбцах;

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

4) уравнение нормировки.

X

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

1

2

3

4

2

3

4 {

апч

«21*1

«12*2

а22г2

а\1-КГ1-К

а21-КГ1-К

VI ат2г2 • * ат1-КГ1-К 3/1 0 5/1

ду1 ду2 ’ * &/1-К дУж

ъ. ъ. ' ' д/к Як

дух ду2 ’ ’ ЗУ ,_к дУж

Т\ Г2 • • П-К ГА1

-а

11

- а

т\

-а

12

-а

т2

-а

11-К

-а

пй-К

ГА1

Ат

Г4

Ж.

^ Ат

к

ду

Ат

1*

'Ат

аиг1

«21*1

VI

дА

дУ\

ъ,

5у1

атКг1

дУк

к

ду

к

'к

О

О

О

а\кгк В\

а2КГК ~ ^2

в„

»1

52

5/-аг

^А2^2

^Ат^т

>К 1

’КК

-8#У

1

Рис. 2. Матрица коэффициентов системы линейных уравнений

Уточнение приближенных значений неизвестных у I производится по формуле

п) = у\п-1) + п),

до получения заданной точности.

Температура газовой смеси определяется из уравнения (10) по методике, описанной в [1, 4, 5].

Дифференциальные уравнения баланса массы и энергии в НРИС и уравнение нестационарной теплопроводности решаются известными методами решения систем обыкновенных дифференциальных уравнений.

МОДЕЛИРОВАНИЕ ХИМИЧЕСКИ НЕРАВНОВЕСНЫХ ТЕЧЕНИЙ В СОПЛАХ

При рассмотрении неравновесного течения по соплу большое внимание уделяется определению потерь удельного импульса на химическую неравновесность %н. В [4] приведены результаты расчетов по определению %н для широкого класса топлив и различных условий. В литературе [1, 5, 6] обращается также внимание на существенное отклонение состава продуктов сгорания при неравновесном течении от равновесного. Для проведения расчетов было выбрано топливо на основе ^04+ЩК-Ы(СН3)2, т.к. в этом случае имеется возможность сравнения результатов не только по потерям удельного импульса, но и по составу продуктов сгорания [1, 6].

Стационарное течение реагирующей смеси рассмотрено на основе концепции НРИС (рис.3).

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

1) по результатам расчетов равновесного течения определяется средний показатель изэнтропы:

Пер = ln

f Pk_^

P J/

г

ln

Pk_ к Р

RT

л

RkTk J

(18)

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

Рис.3. Моделирование течения в сопле

При решении обратной задачи теории сопла движение НРИС описывается уравнением:

1 dp

dL

к 2 j

+---------- = 0.

p dL

(19)

Алгоритм решения основан на методике, описанной в [1] и ведется итерационным методом в следующей последовательности:

1) при известных термогазодинамических параметрах в точке Ьт-1 определяется положение новой точки Ьт = Ьт1 + (Жт-1 + Жт )■ Лt / 2, давление рт в которой определяется по

заданному распределению по длине;

2) затем определяется скорость движения НРИС:

W2 = W2 + R

Wm nm-\^ R0

T_ .PV)

.pv)

m-1

■(pm - Pm-1 ) ;

(20)

3) определяется состав и свойства смеси в точке Lm при p,T = const;

4) на основании уравнения сохранения энергии IT = IПС + W2 / 2 определяется

поправка по температуре [5]

ATm =-5t / Cp, £j = IT - Inc - W /2.

(21)

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

AT < є-j,.

Значение £н определялось по формуле

jn - jn EH = -ид-100%,

Hn

J л і А

(22)

/п т п

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

m

Значение пустотного импульса определяется из выражения:

1п = Жа + R0Ta /(МаЖа), (23)

полученного на основании формулы тяги R = GWa + Fa (ра - рн ), уравнений неразрывности

G = pWF и состояния идеального газа р = pR0T / и.

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

- Н2, СО, Н20, N2 при а < 0,8;

- СО2, Н2О, N при 0,8 < а < 1,2;

- О2, СО2, Н2О, N при а > 1,2.

В [1] отмечается, что а оказывает значительное влияние на величину %н. На рис.4

показаны зависимости %н = f (а) для сопла с диаметром критического сечения 50 мм при различной степени расширения потока в сравнении с аналогичными данными из [4]. На графиках точками показаны данные из [4], а сплошными линиями - результаты численных исследований. Характерной особенностью является наличие максимума потерь вблизи стехиометрического значения а = 1, что объясняется [1, 6] существенной

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

температур, при а< 0,5, течение практически заморожено, что приводит к резкому

возрастанию потерь удельного импульса на неравновесность. Некоторое уменьшение потерь импульса при а < 0,3 можно объяснить образованием конденсата и выделением при этом значительного количества тепловой энергии.

Рис.4. Зависимость £,н = f (а)

При определении состава продуктов сгорания обращается внимание [1, 4, 5] на существенное отличие его от равновесных значений. С точки зрения образования токсичных веществ наиболее интересным является отклонение состава продуктов сгорания на срезе сопла от равновесных значений, особенно для таких веществ, как оксиды азота N0 и углерода СО. На рис.5 для сравнения результатов расчетов равновесного и неравновесного течения приведены графики мольных долей N0 и СО.

толстая линия - неравновесный расчет; тонкая линия - равновесный расчет.

Рис.5. Концентрации токсичных веществ в зависимости от а

В [1] отмечается, что при фиксированном а мольные доли веществ Н2, СО, Н20, N2 СО2 изменяются лишь на 5 - 10% при варьировании давления рк в камере сгорания, диаметра критического сечения dkr и относительного радиуса сопла г / гкг в широких пределах: 2,5 МПа<рк <25 МПа, 25 мм<dkг <250 мм, 1<г/гкг <15. На рис.6 приведена

"универсальная" зависимость мольных долей Н2, СО, Н2О, N2 СО2 для продуктов сгорания топлива М2О4+Н2М-МСН3)2. На рис.7 показаны аналогичные зависимости, полученные при рк =15 МПа, dkr =50 мм.

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

0,25 0,50 1,0 а

Рис.6. "Универсальные" зависимости мольных долей от а

0,2 0,4 0,6 0,8 1,0 1,2 а

г / г'к(г =4 - сплошная линия; г / гкг =8 - точки.

Рис.7. Мольные доли основных компонентов смеси от а

СПИСОК ЛИТЕРАТУРЫ:

1. Пирумов У.Г., Росляков Г.С. Газовая динамика сопел. М.: Наука, 1990. 368 с.

2. Турбулентные течения реагирующих газов: Пер. с англ. // Под ред. П.Либби, Ф.Вильямса. М.: Мир, 1983. 328 с.

3. Кондратьев В.Н., Никитин Е.Е. Кинетика и механизм газофазных реакций. М.: Наука, 1974. 512 с.

4. Алемасов В.Е., Дрегалин А.Ф., Тишин А.П., Худяков В.А. Термодинамические и теплофизические свойства продуктов сгорания: Справочник в 10 томах / Под ред. В.П.Глушко. М.: ВИНИТИ, 1971 - 1980. Т 1. С.10.

5. Алемасов В.Е., Дрегалин А.Ф., Тишин А.П. Теория ракетных двигателей / Под ред. В.Е.Алемасова. - Изд. 2-е, перераб. и доп. М.: Машиностроение, 1969. 548 с.

6. Алемасов В.Е., Дрегалин А.Ф., Крюков В.Г., Наумов В.И. Математическое моделирование высокотемпературных процессов в энергосиловых установках. М.: Наука, 1989. 256 с.

7. Полак Л.С., Гольденберг М.Я., Левицкий А.А. Вычислительные методы в химической кинетике. М.: Наука, 1984. 280 с.

8. Каретто Л.С. Математическое моделирование образования загрязняющих веществ // Образование и разложение загрязняющих веществ в пламени: Пер. с англ. / Ред. Н.А.Чигир. М.: Машиностроение, 1981. С. 84 - 137.

9. Корепанов М.А. Математическое моделирование процессов термической переработки элементов РДТТ: дис. ... канд. техн. наук. Ижевск, 1997. 162 с.

SUMMARY. The set of equations and method of its solutions allowing to determine composition a chemically reacting mixed gas in view of velocity of chemical reactions and the greatest possible number of substances is described. The results of numerical examinations of currents in nozzles of rocket engines are given.