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

УДК 536-12+544.4

Б. Г. Трусов

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

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

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

— экспериментальное исследование кинетики химических превра-ений;

— построение структурной кинетической модели (механизма реакций);

— решение обратной задачи для определения констант скоростей элементарных стадий процесса;

— решение прямой задачи применительно к новым условиям;

— проверка адекватности модели.

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

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

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

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

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

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

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

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

Для идеальной смеси, состоящей из к компонентов — индивидуальных веществ при температуре Т и давлении р = ^ рг, энтропия может

г

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

5 = V Б^Пг = V ($ - Д 1п ^ Пг, (1)

г г г г г 0 р0 г

где пг — содержание ¿-го газообразного компонента в системе, выра-

женное в молях на единицу массы; Б(Рг) — энтропия компонента, отнесенная к одному молю вещества, при том парциальном давлении рг, которое он будет иметь в системе; Бг0 — энтропия г-го компонента при стандартном давлении р0, являющаяся функцией только температуры Т [4]; Я0 — универсальная газовая постоянная.

Если р0 = 1, то с учетом уравнения состояния идеального газа для г-го компонента имеем

5 = ЕЙ0 - ^ 1п «г, (2)

г ^ '

где V — удельный объем системы.

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

Тогда для производной энтропии по времени будет справедливо выражение

dS V^ I оО Г> 1 R0Tnl \ . ^ Cvini

- = £ [So- Ro ln —— Ro)n i + £ —

l

где От — удельная теплоемкость г-го вещества при постоянном объеме.

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

В то же время, второй закон термодинамики утверждает, что производство энтропии, связанное с изменениями внутри систем , всегда неотрицательно: ¿Б > 0. О но вызывается естественными необратимыми процессами, к числу которых могут быть отнесены и химические реакции. ри этом в ней, да е если система изолирована и ее полная внутренняя энергия неизменна, за счет работ химических термодинамических сил имеет место выделение "нескомпенсированной теплоты" [5]. Эту величину, по аналогии с механикой, иногда называют таюке диссипативной функцией [1]:

¿Б 0 ¿Б >> 0 (4)

1ъ = -ТЫ > 0 или И = ы > 0' (4)

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

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

w =

Y LXi J, (5)

V т ¿п где Хг — термодинамическая сила, Jг = —--соответствующий по-

аЬ

ток ¿-го вещества, Ьг — феноменологические коэффициенты.

ценим вид зависимостей для термодинамических сил, воспользовавшись основным уравнением термодинамики (уравнением Гиббса):

таБ = аи + рвю — ^^ ^гапг, = (аи) , ] = ¿, (6)

где и — внутренняя энергия систем .

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

Де Донде Т. [8] на основе этого же уравнения ввел два новых понятия: сродство химической реакции А и степень ее завершенности £. Сродство определяется разностью химических потенциалов реагентов и продуктов реакции, а степень завершенности имеет смысл координаты процесса. С этой точки зрения в работах [5,8,10] сродство трактуется как термодинамическая сила химического процесса, а производная степени завершенности по времени — как скорость рассматриваемой химической реакции:

dS sr-^ Ai d£i

~dt ~ ^ TTt

Y LiXi Ji. (7)

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

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

концентраций всех компонентов системы J = dn/dt. В этом случае мерой химических термодинамических сил в изолированной системе в соответствии с уравнением (6) может стать градиент химического потенциала этого компонента:

x — i d^i = 1 d^i dn (8)

г T dt T dn dt

В справочной литературе для индивидуальных веществ обычно вместо значений химических потенциалов приводятся температурные зависимости приведенной энергии Гиббса G* = — ^¡/T. Если все вещества, входящие в систему, находятся в газообразном состоянии и подчиня тся уравнени состояния идеального газа, то

G = (s°(T) — R ln ^T^) — T (H(T) + А/H(To)), (9)

где H0(T), S0(T) — стандартная энтальпия и стандартная энтропия i-го вещества, А/H0(To) — энтальпия образования вещества, т. е. мера его химической энергии [4].

dG* 1 d^

Вычислив для всех компонентов системы —— =----—, получим

dni T dni

выражения для термодинамических сил и диссипативной функции:

X = — dG* dni = Ro dn w = y^ a Ro i dn^|2

i — dni dt ni dt i ni dt

г 4 '

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

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

l2

т = D, (11)

где l — длина свободного пробега частиц, D,t — коэффициент самодиффузии частиц сорта i.

римем, что длина свободного пробега и коэффициент диффузии равны соответственно [9]:

l = , D, - ^^

4па2и0' 1 ра2П(1'1)*'

где а — сечение соударений, п0 — число частиц в единице объема, Т — температура, шг — молярная масса, р — давление, П(1,1)* — интеграл столкновений.

Тогда для диссипативной функции может быть предложен вид

- = £ *Д (^)2.

г 4 7

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

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

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

У^йдпг = о, з = 1,ш, (13)

г

где ад — стехиометрические коэффициенты, ш — число химических элементов, образующих систему;

б) уравнение состояния идеального газа в дифференциальной форме

^(рь) = ^ ^ДоТ £ п^ , или ьр+Ьр = ^Яо £ п^ Т+ДоТ £ пг;

1 1 1 (14)

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

и = £ итпг = £(иг + Д/И°)пг, (15)

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

С учетом соотношений (4) и (12) закон сохранения энергии представим в виде

£ итиг + £ СтигГ + £ Ьги2 =0, (17)

г г г

где Ьг = ктг —0--феноменологические коэффициенты.

иг

Процедуру нахождения максимальной скорости возрастания энтропии с учетом ограничений (13), (14), (17) выполним, используя метод неопределенных множителей [11]. Для этого составим функцию Ла-гранжа:

Л = £ - —0 1П - —0Пг + £ Т+

гг

+ Хр ^ир + Ьр - —0 £ игТ - —оТ £ и^ + + У^ У^ а^Щ + Хи I £ и„гиг + £ С'огПгТ + £ ЬгП2 I . (18)

] г \ г г г /

осле дифференцирования функции Л по всем переменным (иг, Т, р, и, Xj, Хр, Хи), приравнивания полученные уравнения ну-л и используя частичные подстановки, получаем расчетну систему уравнений

2Ьг ¿и —0Тиг ^ \ • ГТ

= - —о1п и +1^ , г = l,k,

j

у~] а^ги = 0, з = 1, т,

г ( V (19)

ир + Ьр - I —0 £ иг I Т - —0 Ту^иг = 0,

£ итпг + £ Cmn{T + £ Ьгп2 = 0.

г

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

Алгоритм решения этой системы уравнений был реализован в виде компьютерной программы, позволяющей исследовать системы произвольного химического состава благодаря подключению базы данных, содержащей термодинамические и термохимические свойства индивидуальных веществ, а таюке параметры потенциальной функции Леннард-Джонса [12]. Программа позволяет исследовать как адиабатически изолированные, так и изотермические системы, находящиеся при постоянном объеме или давлении.

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

По этой причине в настоящей работе сделана попытка определить степень практической пригодности разработанного метода путем сравнения с результатами расчетов, полученными по какой-либо признанной и проверенной кинетической модели. В качестве такой была выбрана модель горения углеводородных топлив GRI-Mech, опубликованная в сети Internet на странице университета Беркли [16]. Там же приведена информация о ее экспериментальной проверке. Необходимые расчеты по этой кинетической модели проводились с использованием компь -терной программы Kintecus [17].

а рис. 1 представлены результаты моделирования кинетического процесса химических превращений в системе 0,01%CH4+O,02%02+ +99,97%Ar для условий T = 2000 K и p = 0,0984 МПа по описанному выше методу, а на рис. 2 приведены концентрационные кривые, вычисленные с использованием формальной кинетической модели GRI-Mech. Сопоставление результатов по абсолютным значениям концентраций (моль/см3) и по характеру их последовательного образования позволяет сделать вывод о принципиальной работоспособности метода. При этом, конечно, необходимо принять во внимание тот факт, что в описываемой термодинамической модели время можно рассматривать только как условный параметр, однонаправленный с реальным временем.

В процессе моделирования использовались термодинамические и термохимические свойства индивидуальных веществ, опубликованные в работах [4,13]. Параметры потенциальной функции Леннард-

lgcL -8

-9

-10

-11

-11

Ог H« Н20 \

--

СН4 -н ^s. со

С02 он

/ / /// У

0

сн2о

1E-06 2E-06 ЗЕ-06 Условное время

Рис. 1. Результаты моделирования кинетики превращений в системе

CH4 - Ö2 - Ar

0,00Е+00 ЩЕ-ОЬ в,46Е-04- 1J1E-0J 1Щ-03 %29Е-02 trc

Рис. 2. Расчет кинетики превращений в системе СН4 — О2 — Аг с использованием модели ОМ-МесИ

Джонса, использовавшиеся для вычисления термодинамических сил, первоначально принимались по данным работы [14]. Значения этих параметров, в частности диаметры столкновений а, обычно рассчитываются по данным о теплопроводности, вязкости или по р-у—Т-измерениям для реальных газов. Разные источники могут приводить для них отличающиеся значения [9], поэтому при сопоставлении данных, приведенных на рис. 1 и рис. 2, значения а варьировались в пределах 20% для того, чтобы оценить чувствительность результатов к изменению величин термодинамических сил процесса.

Наилучшее качественное совпадение результатов моделирования и кинетических расчетов по формальной модели ОЫ-МесЬ, как видно, достигается на начальном этапе процесса. В дальнейшем поведения концентрационных кривых различаются. Можно предположить две причины для объяснения этого обстоятельства:

а) "условное время" моделирования течет намного быстрее реального;

б) кинетическая модель ОЫ-МесЬ не предназначена для адекватного описания процесса вблизи равновесия.

снованием для второго предполо ения могут слу ить результаты расчета концентраций в условиях равновесия для рассмотренного примера [7]. В табл. 1 приведены равновесные концентрации основных компонентов, полученные с помощью представленной модели (они строго соответству т условиям равновесия) и в результате установления решения по кинетической модели ОЫ-МесЬ.

Таблица 1

Компонент Равновесные концентрации, 1013, см 3

Моделирование Расчет GRI-Mech

O 1,59 3,25

Ö2 26,37 36,5

H 3,05 6,34

OH 12,53 17,8

H2 17,4 25,02

CO 37,94 53,16

CO2 31,02 29,47

H2O 669,5 658,5

На упомянутой вше странице Internet [16], посвященной кинетической модели GRI-Mech, приведены также некоторые экспериментальные данные, позволяющие количественно оценить точность предлагаемого термодинамического метода моделирования. Для уже рассмотренной системы CH4—O2—Ar была выполнена серия расчетов в температурном диапазоне 2000...2500K. На рис.3 представлены рассчитанные при температуре 2000 K концентрационные кривые, моль/см3, от начала процесса до установления равновесия. Видно, что концентрации CH3 и CO имеют максимум. На рис. 4 сопоставлены результата этой серии расчетов для CH3 и соответствующие опытные значения, полученные экспериментально в ударной трубе [15]. Видно, что описанный термодинамический подход к моделированию

Igq

-7 -8

-9

-10

110 1,0Е-Об Щ-Од Щ-06 Щ-06 Ь Условное время

Рис. 3. Установление равновесия в системе СН4-02-Лг

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

Вторая группа расчетов для сопоставления с экспериментальными данными была выполнена применительно к системе 0,4%СН4 + 5%02 + 94,6%Аг для условий Т = 1941 К,

р = 0,258 МПа и Т = 1821 К, Рис. 4. Максимальные значения концен- _ тт _

трации СНз: Р = 0,1696 МПа. На рис. 5 пока-

прямая — результат моделирования, зан характер изменения концен-

точки — эксперимент [15] трации СО для данных условий,

а на рис. 6 — соответствующие экспериментальные значения, полученные в ударной трубе [15]. Если снова принять во внимание, что условное расчетное время не пропорционально реальному времени, то из сравнения результатов видно, что абсолютная погрешность моделирования не превышает 30%.

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

CHjmax,PPm

0,40 0,М ОМ 0,US 0,48 1000¡Т

Ссо-108, моль/см3

\ 7

г

//

О Щ-06 10Е-05 1Д-05 2, ОЕ-05 2,5Е~0,5 t Условное время

Рис. 5. Графики изменения концентрации CO:

1 — T = 1941K, p = 0,258 МПа, 2 — T = 1821K, p = 0,1696МПа

С 7

6

5 4 J

г 1

О 50 100 150 ZOO 2JO 300t,nc

Рис. 6. Максимальные значения концентрации CO по данным работы [15]:

• — T = 1941K, p = 0,258 МПа, о — T = 1821K, p = 0,1696МПа

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

1. Циглер Г. Экстремальные принципы термодинамики необратимых процессов и механика сплошной среды. - М.: Мир, 1966. - 134 с.

2. Пригожин И. Введение в термодинамику необратимых процессов. - Ижевск: НИЦ "Регулярная и хаотическая динамика", 2001. - 160 с.

со-108,МОЛЬ/СМ3

• •

• •«

• • • • О

• • *| О Г)

N о с

• о

- • • •• п О о

3. Гленсдорф П., П р и г о ж и н И. Термодинамическая теория структуры, устойчивости и флуктуаций. - М.: Мир, 1973. - 280 с.

4. Термодинамические свойства индивидуальных веществ: Справоч. в 4-х т. / Л.В. Гурвич, И.В. Вейц, В.А. Медведев и др. - М.: Наука, 1978-1982.

5. Пригожин И., Д е ф э й Р. Химическая термодинамика. - Новосибирск: Наука, 1966. - 509 с.

6. Де Гроот С., Мазур П. Неравновесная термодинамика. - М.: Мир, 1964.-456 с.

7. Применение ЭВМ для термодинамических расчетов металлургических процессов / Г.Б. Синярев, Н.А. Ватолин, Б.Г. Трусов, Г.К. Моисеев. - М.: Наука, 1982.-263 с.

8. Де Донде Т., Ван Риссельберг П. Термодинамическая теория сродства (книга принципов): Пер. с англ. / Под ред. В.М. Глазова. - М.: Металлургия, 1984. -136 с.

9. Гиршфельдер Дж., Кертисс Ч., Берд Р. Молекулярная теория газов и жидкостей. - М.: ИЛ, 1961. - 929 с.

10. Пригожин И., Кондепуди Д. Современная термодинамика. От тепловых двигателей до диссипативных структур: Пер. с англ. - М.: Мир, 2002. - 461 с.

11. Бертсекас А. Условная оптимизация и методы множителей Лагранжа. - М.: Радио и связь, 1987. - 400 с.

12. Трусов Б. Г. Программная система TERRA для моделирования фазовых и химических равновесий в плазмохимических системах // 3 Межд. симпозиум по теор. и прикл. плазмохимии. - Иваново, 2002. - С. 217-220.

13. Иориш В. С., Белов Г. В., Юнгман В. С. Программный комплекс ИВТАНТЕРМО для WINDOWS и его использование в прикладном термодинамическом анализе: Препринт № 8-415. - М.: ОИВТ РАН, 1998. - 56 с.

14. S v e h l a R. A. Viscosities and Thermal Conductivities of Gases at High Temperatures // NASA TR R-132. Washington,Gov.Print.Off.,1962. - 25 p.

15. Yu C. L., Wang C., Frenklach M. Chemical Kinetics of Methyl Oxidation by Molecular Oxygen // J. Phys. Chem. - 1996. - V. 99. - P. 14377.

16. http://www.me.berkeley.edu/gri_mech/.

17. I a n n i C. J. - http://www.kintecus.com/.

Статья поступила в редакцию 25.04.2005

Борис Георгиевич Трусов родился в 1940 г., окончил в 1963 г. МВТУ им. Н.Э. Баумана. Д-р техн. наук, профессор, заведующий кафедрой "Программное обеспечение ЭВМ и информационные технологии" МГТУ им. Н.Э. Баумана. Специализируется в области компьютерного моделирования высокотемпературных процессов с фазовыми и химическими превращениями, компьютерного моделирования в плазмохимии, энергетике, двигателестроении.

B.G. Trusov (b. 1940) graduated from the Bauman Moscow Higher Technical School in 1963. D. Sc. (Eng.), professor, head of "Software and Information Technologies" department of the Bauman Moscow State Technical University. Specializes in the field of computer-aided simulation of high-temperature processes with phase and chemical conversions, computer-aided simulation in plasma chemistry, power engineering, engines building.