Стационарная модель твердофазного горения с учетом термических напряжений. Асимптотический анализ
А.Г. Князева
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия
Получено асимптотическое решение стационарной задачи о распространении фронта реакции в твердой деформируемой среде в приближении малых деформаций. Выявлена область параметров модели, в которой существуют быстрые режимы превращения. Построено аналогичное решение задачи в более полной постановке, учитывающей уравнение неразрывности и уравнение движения в общем виде и не использующей условие малости деформации. Показано, что в этом случае задача сводится к поиску массовой скорости горения, а область параметров, где существуют быстрые режимы превращения, существенно расширяется. Обнаружены, по крайней мере, три типа быстрых режима превращения в твердой фазе.
1. Введение
Исследование физико-химических превращений в твердой фазе показало, что их характеристики в существенной степени зависят от внутренних напряжений и деформаций, т.е. от напряжений и деформаций, сопровождающих превращения. С одной стороны, напряжения и деформации могут непосредственно влиять на скорость превращения. С другой стороны, их влияние может быть обусловлено связанным характером различных физических процессов, в том числе теплопереноса и деформирования. Если реакция в твердой фазе является экзотермической, процесс послойного распространения превращения по веществу называют твердофазным горением. К числу таких процессов принадлежат само-распространяющийся высокотемпературный синтез, низкотемпературные радикальные реакции, многие реакции твердофазной полимеризации, реакции экзотермического разложения в твердой фазе и др. Если зона реакции, распространяющаяся по веществу, характеризуется большими градиентами температуры, то напряжения и деформации, сопровождающие реакцию, имеют, в основном, термическую природу. Для описания таких превращений следует использовать связанные модели. Под структурой волны твердофазного горения понимают характер распределения в зоне реакции температуры, степени превращения, плотности и скорости среды.
В ряде работ показано, что в случае экзотермических превращений связанные модели допускают существование различных режимов [1-5]. Скорость медленных стационарных режимов превращения определяется переносом энергии за счет теплопроводности. Скорость быстрых режимов превращения определяется переносом энергии волной механических возмущений. Учет особенностей среды, в которой возможны превращения, структурных неоднородностей, деформаций и напряжений различных типов, изменения плотности и других свойств в ходе реакции приводит к появлению новых закономерностей.
В настоящей работе проводится сравнение стационарной модели горения, построенной для случая малых деформаций, с моделью твердофазного горения, где условие малости деформаций не используется.
2. Простейшая модель твердофазного горения
Математическая постановка простейшей задачи твердофазного горения с учетом только термических напряжений и в приближении малых деформаций имеет вид [1]:
дТ .
СсР------- = ЛТ ------
Е дг Т дх2
д Т - 3КаТТ — + Т дг
+ак ГФ1( у)Ф2(Т),
(1)
© Князева А.Г, 2004
дt
= кг Фі( у)Ф2(т X
р
де Э а,
дt2
Эх2
(2)
(3)
а її — (Л + 2ц)е — ЗКат (Т — Ту), а 22 = азз = Ле — ЗКат (Т — Т0),
Ф2(т)=ехр[^--¿Тт
х —+<~: т = т>, у = 0,
(4)
(5)
х : Т = Ть, ^ = 1, (6)
где Т — температура; у — степень превращения (или массовая доля продукта реакции); е — компонента тензора деформации в направлении распространения фронта; ст^ — компоненты тензора напряжений (отличные от нуля); х — пространственная координата; t — время; АТ — коэффициент теплопроводности; р — плотность; се — теплоемкость при постоянной деформации; А, ц — коэффициенты Ламэ; К — изотермический модуль всестороннего сжатия (К =А + 2ц/ 3); аТ — линейный коэффициент теплового расширения; ф1 (у) — кинетическая функция; кг, Еа, аг — формально-кинетические параметры; Я — универсальная газовая постоянная; Т0 — температура недеформированного состояния; Ть — температура продуктов реакции, которую требуется определить. Условие для температуры в области продуктов эквивалентно условию дТ/ дх = 0.
В такой модели предполагается, что экзотермическая реакция может быть описана простой суммарной схемой
А ^ В,
фронт реакции — плоский, термические напряжения в волне горения преобладают над остальными (концентрационными, структурными). В работе [1] решение этой задачи в форме бегущей волны найдено с использованием качестве масштаба температур температуры продуктов реакции для данной модели (т.е. с учетом работы напряжений). Вид решения в этом случае оказывается не очень удобным для анализа, так как требует многочисленных преобразований переменных. Оказывается, что явное введение температуры продуктов реакции Ть, отличной от температуры продуктов тепловой теории горения, более удобно для построения аналитического решения.
В системе координат, связанной с фронтом, движущимся вправо со скоростью уп, вместо (1)-(3) имеем
дт дт
сер|----------vn-------
е 1 ді п дх
= Л
д 2т
дх 2
де де Л
— 3Катт 1 — - Vn "дх I + бкгф1(у)ф2 (т)’
Iу — V п Iу = к г Фї( у)ф2(т),
ді ох
ґ д 2е
2п Л
~ о2е + 2 д2е
—т- — 2Vn---------+ Vn —“
ді2 дхді дх
д а1
~дх2
Перейдем в задаче к безразмерным переменным
е = . т — т0
тЬ0 — т0
•и х і ап
£= —, т = -, 5 = ^ (7)
х і* а*
и введем новое определение для деформации (которая и так является безразмерной величиной) е = —, где
е»
а* , —1
е* = -——, і* = кг ехр
( Е Л
^ть0у
Л + 2ц
= у]Кті* , а* = 3Кат(ть0 — т0),
тЬ0 = т0 +
-б-.
сер
Тогда стационарная задача, которая получается из нестационарной при отбрасывании производных по времени, принимает вид
.¿е
d2е
- ¿е
¿у
’
— + ®(Є + а)^— — ^,
— V ^? = Ф1( У)Ф2(е),
(8)
Ф2 (е) = ехр
= ехр
1 — е
в(Є+ а)
(1 — е)(1 + а) у (е + а)
¿2е ,
2’
где
^ —— +^: е = 0, у = 0,
(ЗКат )2 ть0 — т> ю = ^————----------у =
(9)
(10)
(11)
(12)
Л + 2ц сер
се ¿ть0
Еа бг
а=
ть0 + т0
а=
рх*
в= ¿ть0
V = -
і*2 (Л + 2ц) х*/ і
В области продуктов реакции (на бесконечном удалении от зоны превращения) все возмущения затухают:
е ¿е
£ ——^: у = 1^ -г = °
(13)
так что температура продуктов 0 = 0Ь является искомой величиной,
р
еь= ть — т0
тЬ0 — т0
Уравнение движения (10) сразу интегрируется
е
1 —а2.2
и приводит к уравнению теплопроводности
— V
1+
ю(е + а)
1 — а2.2
¿е = 4е — ¿у ¿£ = ¿^2 '
(14)
Перейдем в стационарной задаче (9), (14) к новой переменной ^1 = Ьу, что эквивалентно введению нового пространственного масштаба, т.е.
^1 ^
кт/ ^
В результате придем к уравнениям
1+
ю(е + а)
1—а
ае = ^6— (15)
— тУ~г = Ф1(у) Ф2(е),
2 —1 2 2
где т = .у ; а1 = а V с теми же граничными условиями:
^1 — +^: е = 0, у = 0,
¿е
^1 — —^. у =1, —- = 0. ^1
¿е
(16)
Введем р =-----, что позволяет понизить порядок
^1
системы дифференциальных уравнений (15)
1+
ю(е + а)
1 — а1
— ¿у ¿е — ¿е ,
¿у
— ту17Г Р = Ф1( у)Ф2 (е).
¿е
Аналитическое решение задачи проведем методом сращиваемых асимптотических разложений.
Задача во внешней области (инертная задача)
1+
ю (е + а) 1—а
¿Р ¿е ,
е — 0: р — 0, у = 0. (17)
Решение задачи (17) имеет вид
Р=
¿е
1+
юа
1 — а1
е—-
юе2 1 — а1 2
у=0
(18)
1+
юа
1—а
е—
юе2 1 — а1 2
Не представляет труда записать общий интеграл этого уравнения. Причем постоянная интегрирования будет находиться из условия сращивания с внутренним реше-
нием. Выражение для температуры не приведено, так как оно далее не используется.
В частном случае ю = 0 решение внешней задачи имеет вид
р = -0 и 0 = С1 ехр(-^1).
Для решения внутренней задачи переходим к внутренней переменной
1 -0
і = -
У
, ¿е = —уйі.
(19)
Тогда задача во внутренней области (в зоне химических реакций) принимает вид
-у
1+
ю(1 — уі + а)
1 — а1
= — ¿Р+¿у , (20)
¿і ¿і
¿у
т р = ф1( у)ф2(і), ¿і
1—е
Ь .
где
Ф2(і) = ехр
: Р — 0, у = 1,
(1 + а) і
1 + а — уі
(21)
(22)
(23)
Представление функции (23) в виде ряда Тейлора по степеням малого параметра у с точностью до слагаемых второго порядка малости будет иметь вид
Ф2(і, У) = Ф2(і, 0)
1—
+.
(24)
1+а
Решение внутренней задачи (20)-(22) ищем в виде
2
Р = Р0 +УР1 + У Р2 + .->
у = У0+уу1+у2 у2+...,
2
(25)
т = т0 + ут1 + у т2 +....
Подставляя (25) в (20), найдем для нулевого прибли-
жения
Фр = ¿у0 ¿і ¿і ’
т0Р0^У° =(1 — Уo)фo,
для первого
1+
¿і
ю(а +1) 1 — а,
¿Р1 + ¿у1 ¿і ¿і
¿у1 ¿у0
т0 Р0~ГГ +(Р1т0 + Р0т1)~1-¿і ¿і
= (1 — У0)ф1 — У1ф o,
для второго юі
¿Р2 + 4у! ¿і 4і
У
Рис. 1. Зависимость температуры и скорости для низкоскоростного режима от параметров в модели с малыми деформациями: а < 1; ст = 0.2 (1),
0.4 (2), 0.6 (3), 0.8 (4)
т0 Ро
—+(Р1^0 + Р0т)^уу-+
dt
+ (Р2 т0 + Р1т1 + Р0т2>
&
4у0
= (1 -У0)Ф2 -У1ф1 -у2Ф0 и для последующих
- + 4^ = 0,
йг йг
т0 Р0<йг +(Р1т0 + Р0 т)
+ (Р2 т0 + Р1да1 + Р0 «2)
йУк-1
йг
йУк - 2
йг
+ (Р3 т0 + Р2 т1 + Р0т3 + т2 Р1)-йУк=3
+... +
+ (т0 Рк + т1 Рк-1 + т2 Рк - 2 + ... + ткР 0)
йУ0
йг
: (1 - У0)Ф к - У1Ф к-1 - У 2Ф к - 2 -...- У к Ф 0>
где
Ф 0 =ф 2(г, у = 0); Ф к = д к ф2(г ’7)
ду к
у=0
Решение внутренней задачи в нулевом приближении с учетом условия
1 -0
Ь
У
имеет вид
Р0 = у0 -1 т0(1 - у0) = е -ь - е-.
(26)
1
Переходя в (26) к внешнему пределу г ^ — и при-
У
равнивая затем внутреннее решение к внешнему, записанному для г ^ -ь, придем к системе равенств
юст
1 - а-«
0 Ь + ---0 Ь -1 = 0
2 1 - а^
т0
' ехр
1-0
1 - 0Ь
ехр
У
(28)
Уравнение (27) имеет два типа положительных решений 0Ь < 1 для а1 << 1 и 0Ь2 > 0Ь1 > 1 для а > 1, зависящие от параметров ю, ст, а! [1]. Уравнение (28) дает для этих значений температуры соответственно т0 < 1 и т0 > 1. Зависимости температуры продуктов и скорости от коэффициента связанности для низкоскоростного режима показаны на рис. 1.
В случае высокоскоростного режима из (27) имеем
0 Ь1,2 =
1 -Ю1СТ±д/(1 - Ю^СТ)2 - 2Ю!
ю1
где ю
ю
а
1
. Оба значения температуры действи-
тельные, если выполняется неравенство
Рис. 2. Кривые, разделяющие область параметров модели на подобласти, в которых существуют вещественные (выше ю+ и ниже ю_) корни уравнения (27). Ниже кривой ю* корни положительные. а > 1
У
Рис. 3. Зависимость температуры реакционной зоны (а), температуры продуктов (б) и массовой скорости горения (в) от ст. ю1 = 0.1 (1), 0.15 (2), 0.2 (3); а1 > 1
Ю>2ст2 - 2Ю!(ст + 1) +1 > 0.
Этому условию удовлетворяют значения параметров, лежащие выше кривой ю+ или ниже кривой юна рис. 2. Хотя бы один корень уравнения (27) в этом
случае положительный, если ю1 < ю* = —. Следова-
ст
тельно, высокоскоростные режимы превращения существуют только в области параметров ниже кривой
ю-, т.е. при малых ю (что соответствует реальной области изменения этого параметра) и больших а1.
Из двух значений температуры продуктов при а > 1 ниже кривой ю- для расчета скорости выбираем меньшее значение. Как показывает численное исследование задачи [1], обе температуры 0Ь2 >0Ь1 > 1 (для а> 1) имеют физический смысл, что фактически указывает на двухтемпературную структуру фронта горения. Аналогичный результат получается и в более сложных случаях [1-3]. Обе температуры и массовая скорость горения существенно зависят от параметров задачи, что показано на рис. 3.
Очевидно, при ю = 0 имеем 0Ь = 1 и т0 = 1, что совпадает с известными решениями.
3. Обобщение модели на случай конечных деформаций
Рассмотренная выше связанная модель твердофазного горения имеет место, если деформации, скорости, ускорения считаются малыми, что позволяет не рассматривать уравнение неразрывности в явном виде. Примем, что деформации, ускорения, перемещения не малы. В этом случае одномерная модель твердофазного горения включает [4, 5] уравнение неразрывности
йр йи .
— + Р— = 0,
й йх
уравнение движения йи дстт 1
р^“= 3 ,
йг дх
уравнение баланса внутренней энергии в форме уравнения теплопроводности
Р се
йТ = А д 2Т йе +
— АТ 3атТ +
йг Т дх2 йг
+ & к г Ф1( У )Ф 2(Т X
уравнение химической кинетики для суммарной реакции
= кг ф1( У )ф 2(Т)
и линейные соотношения между компонентами тензоров напряжений и деформаций, записанные в приращениях:
йсту = 2цйеу +5у (Айекк - 3КйТ),
так что
йст. 1 ^ ч йе ^ йТ
----= (А + 2ц)—- 3Ка т —.
йг йг йг
В этой системе уравнений и — скорость точек среды.
Следовательно, в системе координат, связанной с фронтом реакции, движущимся вправо со скоростью ип, имеем стационарную модель
dp 2 dv
-да — = -p2--------------,
dx dx
dv da,,
- да— = -
- да ce
dx dx
dr
dx
_A,
d2r
dx2
„ _ К йе Qг йу
+ 3атТ—т------------т —,
Р йх Р йх
- т йУ = ркгф1(У)ф2(Т),
йх
= (А + 2ц)—-3КаТ —, йх йх йх
где т = р(ип - и) — массовая скорость «горения».
В рассматриваемом случае можно принять
(29)
^kk
;е = Ро/Р-1,
что позволяет замкнуть и частично проинтегрировать систему (29) с заданными условиями в реагентах
(X — +го)
т = т,, р = р0, V = 0, е = 0, ап = 0, у = 0 и с условием затухания возмущений в продуктах реакции (х — —^). Действительно, так как
de
dx
Ро dp
p
2 dx
то из второго и пятого уравнений системы (29) находим
йи р0 йТ
т— = —2 (А + 2ц)---------+ 3КаТ —.
йх р2 йх йх
Используя для нахождения производной плотности уравнение неразрывности (первое уравнение системы (29)), придем к уравнению
dv
dx
3даКат
dr
да -p о (А + 2ц) dx
Аналогично имеем
3Kar
dT
_1Ф _________________________
p2 dx да2-p0(А + 2ц) dx
(30)
(31)
Следовательно, уравнение теплопроводности (третье уравнение системы (29)) заменяется уравнением
dr А
- да ce--------------------_ А
dx
d2T
dx2
(3KaT) Т даpо dr Qr dy
p(да2-po(A + 2ц)) dx p dx
(32)
В безразмерных переменных 0, p_p/p0, v _ T/t* , £ из (30)-(32) имеем
V e*
d0
dV ________________
d£ _ 1 - (aV)2 d£’
1 dp _ _ e* d0 p2 d£ _ 1 - (aV)2 d^,
(33)
(34)
- V
^1+^£+4 p-1A 1 - (aV)2
d20
d!j
d0
d^
2 + Ф1( У)ф 2(0)
где все параметры аналогичны предыдущей задаче, но включают начальную плотность вещества
ю _
V _
(3arKГ Тьо - To A + 2ц CeРo
да
a2 _-
e* _
pDyl Kr ft*
3KaT (Тьо - То)
pp Kr
A + 2ц t*
A + 2ц
«Масштаб» деформаций е» становится теперь параметром и представляет собой произведение термической деформации aT(Tb0 -T0) и отношения скоростей распространения объемной и продольной механических волн 3K/ (A + 2ц).
Уравнение кинетики (9) не изменяется. Соотношение между напряжением s и деформацией e также остается прежним (11).
Полагая, что е» ~ const, ю ~ const, из (33), (34) с учетом условий в реагентах находим
_ 1 _ e»V
р =--------------, v = -
1+
e»0
1 - (aV)2
0,
(36)
1 - ^ )2
т.е. плотность и скорость являются функциями температуры.
Следовательно, вместо уравнения (35) имеем
-V
1+
ю (0 + a)
1 - (aV)2
(
1 + -
e*0
Л'
1 - (aV)2
d0
dt
d20
2 +Ф1( У)Ф2(0)-
(37)
Некоторые свойства этого уравнения (которое в частном случае ф1 ( у) = 1 представляет собой уравнение Льенара), позволяющие судить о качественном характере решений, описаны в [5].
Для нахождения скорости фронта, как и выше, в задаче (37), (9) с условиями (12), (13) используем метод сращиваемых асимптотических разложений.
После перехода к переменной ^1, а затем к переменной р внешняя задача будет включать уравнение
1 +
ю (0 + a) 1 - a1
(
1+
e*0 1 - a1
Y
dp P _— P d0
и условия, аналогичные задаче (17). Решение внешней задачи с учетом условия в реагентах имеет вид
Рис. 4. Зависимость температуры продуктов 8* для дозвукового режима, а < 1, а = 0.3. 1 — ю = 0.5, а = 0.1; 2 — ю = 1.5, а = 0.1; 3 — ю = = 1.5, а = 0.7
Рис. 5. Зависимость температуры продуктов для сверхзвукового режима, а > 1
Р = —
+ -
1 + •
юа
1 — а1
0 +
1 +
8*а
1 — аі
1 — а і 2
ю8* 03
(1 — а0 3
(38)
Внутренняя задача включает уравнения и условия в продуктах реакции
-у-1 +
ю
1 — а1
(а +1 — у*)
1+
8*
1 — а
-(1 —у*)
= _ & + ^
& & ’ т ^Ур = ФД у )ф 2(і X
(39)
1 — 0 ь У
Р — 0, У = 1.
Еще одним условием является условие сращивания внешнего и внутреннего решений, которое позволяет найти скорость фронта (в данном случае массовую скорость горения) и температуру продуктов реакции. Из системы уравнений (39) видно, что в нулевом приближении решение внутренней задачи не зависит от параметра е*. Более того, так как, по определению, е* << 1, то второе приближение также не будет зависеть от этого параметра, так как слагаемыми порядка е*у в решении можно пренебречь как величинами второго порядка малости.
Аналогично предыдущему для скорости горения и температуры продуктов в нулевом приближении найдем уравнения (28) и
1+
юа
1 — а1
0ь +
1+
8*а
1 — а1
ю 0 ь
1 — а1 2
+
ю8* 0
(40)
В зависимости от величины и знака дискриминанта этого уравнения, оно имеет одно или три действительных решения. Анализ показывает, что в реальной области изменения параметров это уравнение имеет либо один положительный корень, либо два положительных корня, либо три. Если а1 < 1, то значения 0Ь < 1 дают массовую скорость горения для дозвукового режима, причем как значения температуры, так и значения скорости слабо зависят от параметра е* (рис. 4). Иная ситуация наблюдается в сверхзвуковой области. Если а1 мало превышает единицу, то существует единственный положительный корень уравнения (40), который может быть как больше, так и меньше единицы. Этот режим горения является однотемпературным, но температура 0 Ь может быть как больше, так и меньше единицы (рис. 5). Условно, если 0Ь < 1, имеем низкотемпературную низкоскоростную твердофазную детонацию. В противном случае т >> 1. При а1, превышающем некоторое значение, и е*, меньшем некоторой величины (конкретные значения зависят от других параметров задачи), уравнение (40) имеет три положительных корня, т.е. мы получаем «трехтемпературный» режим горения, как и в представленной выше задаче. Меньшее значение температуры слабо зависит от е*, большее—резко возрастает с ростом е* и с уменьшением ю. Так, при е* = = 0.05, ст = 0.3, ю = 1, а1 = 10 имеем 0Ь1 ~ 1.164 и 0Ь2 ~ 17.63. При увеличении е* до 0.1 находим 0Ь1 ~ 1.16 и 0Ь2 ~ 19.3. Третий корень 0Ь3 >> 1 уменьшается с ростом е*. Если е* = 0.2, то 0Ь1 ~ 1.01, но второй и третий корни исчезают. Существует область параметров, где все три температуры имеют конечные
е*
значения. Так, при ст = 0.3, е*1 =-------= 0.1, ю1 =
8
а1 —1
а1 — 1
0.45 имеем 0ь1 = 2.12, 0ь2 = 3.44, 0
ь3
1.97. Область изменения ю1, где существуют трехтем-
+
* — * ь =
+
пературные режимы, расширяется с уменьшением е*. Для ст = 0.3 и е* > 0.11 существует только однотемпературный режим твердофазного горения.
Вопрос об устойчивости тех или иных режимов пока остается открытым и требует специального исследования. Физический смысл двух или трех температур 0Ь, по-видимому, можно будет выяснить при численном решении задачи.
Как было показано выше, характер изменения плотности и скорости в волне горения определяется характером изменения температуры (см. формулы (36)).
В задаче с малыми деформациями изменение плотности (или изменение удельного объема) определяется изменением первого инварианта тензора деформаций
Переходя к безразмерным переменным и учитывая решение стационарной задачи, получим тот же результат, что следует из (36).
Таким образом, в работе найдены приближенные аналитические решения задач твердофазного горения
для случая малых деформаций и при явном изменении плотности в зоне реакции. Показано, что и в том, и в другом случаях существуют как дозвуковые, так и сверхзвуковые стационарные режимы превращения.
Работа выполнена при финансовой поддержке фонда РФФИ, грант № 03-01-00074.
Литература
1. ТимохинА.М., Князева А.Г. Режимы распространения фронта реакции в связной термомеханической модели твердофазного горения // Хим. физика. - 1995. - Т. 15. - № 10. - С. 85-100.
2. Князева А.Г, Дюкарев Е.А. О режимах твердофазного разложения
одиночных кристаллов инициирующих взрывчатых веществ // Физ. мезомех. - 2000. - Т. 3. - № 3. - С. 97-106.
3. КнязеваА.Г., ДюкаревЕ.А. Термомеханическая модель автовол-нового распространения низкотемпературных реакций с учетом разрушения // Физика горения и взрыва. - 1998. - Т. 34. - № 4. -С. 84-93.
4. Князева А.Г Об одной причине существования быстрых режимов твердофазных превращений // Матер. XII Симпозиума по горению и взрыву «Химическая физика процессов горения и взрыва». -Черноголовка, 2000. - Ч. III. - С. 89-91.
5. Князева А.Г. Решение задачи термоупругости в форме бегущей вол-
ны и его приложение к анализу возможных режимов твердофазных превращений // ПМТФ. - 2003. - Т. 44. - № 2. - С. 14-26.
_—-1.
akk _о
p
Stationary model of solid-phase combustion with regard to thermal stresses. Asymptotic analysis
A.G. Knyazeva
Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia
The asymptotic solution of a stationary problem about the reaction front propagation in a solid medium under deformation is obtained in the small strain approximation. The region of model parameters, where fast modes of transformation exist, is revealed. A similar problem solution is constructed in a more complete statement, which accounts for the equation of continuity and equation of motion in the general form and uses no condition of small strain. It is shown that in this case the problem reduces to finding the mass rate of combustion, while the parameter region, where fast modes of transformation exist, expands considerably. At least three types of fast modes of transformations in the solid phase are revealed.