Modeling of Stefan-type problems under conditions of thermal decomposition of binders in thermal protection composite materials

Full Text

Введение

Композиционные материалы (КМ) находят широкое применение в системах тепловой защиты элементов конструкций скоростных летательных аппаратов (ЛА) с 1970-х годов. Это обусловлено особенностями их состава, которые позволяют при высоких температурах отводить значительное количество теплоты благодаря эндотермическому разложению органических связующих с образованием газовой фазы и пористого коксового остатка. Через этот остаток происходит фильтрация газов, которая сопровождается поглощением теплоты за счет конвекции, вдува газов в газодинамический пограничный слой, оттеснения его от стенки, существенного снижения тепловых потоков, а также за счет объемной теплоемкости материала [1–6].

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

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

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

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

Следует отметить, что до настоящего времени не решена не только изложенная комплексная проблема, но и задача моделирования генерации массы газа в зоне разложения, пригодная для произвольного теплозащитного КМ. Существующее моделирование осуществляется с привлечением химической кинетики связующего конкретного КМ, что требует знания не только точного химического состава связующих и их остатков после разложения, но и всех термодинамических характеристик: температур разложения, энергий активации для каждой реакции, коксовых чисел, коэффициентов газификации и т. п. Такой подход к моделированию газообразования является тупиковым, поскольку не позволяет в процессе проектирования теплозащиты из КМ проводить сравнительный анализ эффективности теплоотвода для различных материалов [7–9].

Ранее задачи типа Стефана с одной подвижной границей без учета газообразования и фильтрации рассматривались в ряде работ, например [1, 10]. В последнее время опубликовано значительное число работ, посвященных решению задач типа Стефана в процессах сушки и задачам теплоизоляции, связанным с сохранением теплоты, как с учетом, так и без учета фильтрации паров [11, 12].

Для преодоления трудностей, связанных с формализацией химической кинетики разложения связующих теплозащитных КМ, в работах [13, 14] предложена идея использования плотностей и температур начала ($\rho _{b}$, $T_{b}$) и окончания ($\rho _{e}$, $T_{e}$) разложения связующих, которые имеются в паспорте каждого теплозащитного КМ. Это позволяет сформировать характеристики закона разложения связующих КМ в виде зависимости плотности КМ от температуры в зоне разложения, ограниченной двумя нестационарно подвижными границами с температурами $T_{b}$ и $T_{e}$. В настоящей статье осуществляется развитие данного подхода, которое заключается в формулировке характеристик закона разложения и их интеграции в комплексную физико-математическую модель тепломассопереноса.

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

1. Математическая формулировка задачи

Рассматривается задача получения аналитического решения нелинейной проблемы тепломассопереноса, заключающейся в определении нестационарного температурного распределения $T(x, x_{b}(t), x_{e}(t), t)$ в теплозащитном КМ. Учитывается термическое разложение связующих в подвижной зоне, ограниченной границами $x_{b}(t)$ и $x_{e}(t)$ с температурами $T_{b}(x,t)$ и $T_{e}(x,t)$ соответственно (при $x_{e} < x_{b}$, $T_{b} < T_{e}$), а также определяется скорость движения $\dot{x}(t)$ этой зоны на основе баланса подводимых и отводимых тепловых потоков. Модель также включает газообразование и неизотермическую фильтрацию в пористом коксовом остатке. При формулировке и решении задачи приняты следующие допущения:

• границы начала $x_{b}(t)$ и окончания $x_{e}(t)$ разложения связующих первоначально определяются интерполяцией температурного распределения, полученного из решения задачи теплопроводности без фазовых превращений;
• фазовые превращения происходят только внутри зоны, ограниченной границами $x_{e}(t)$ и $x_{b}(t)$;
• гетерогенные химические реакции в пористом коксовом остатке отсутствуют;
• фильтрация газов, обусловленная перепадом давления между зоной разложения и внешней границей тела (давлением в пограничном слое), является стационарной;
• ввиду малой скорости газов в зоне разложения их давление и плотность принимаются равными давлению и плотности торможения;
• хотя область разложения нестационарна и подвижна, в каждый момент времени границы $x_{b}$ и $x_{e}$ начала и окончания разложения движутся с одинаковой скоростью $\dot{x}(t)$, что позволяет считать распределение температур и плотности газа в этой области квазистационарным.

Формулируются следующие задачи.

1. Задача определения нестационарного распределения температур $T_{2}(x, t)$ в пористом коксовом остатке (фаза «2») с учетом неизотермической фильтрации газов. Для этого зафиксируем в текущий момент времени $t_{f}$ границы $x_{fb}$ и $x_{fe}$ и времена $t_{e}$, $t_{b}$ достижимости этих границ.
   Область определения $0 < x < x_{fe}$:
   \[\begin{equation}
   \tag{1}
   \frac{\partial T_{2}}{\partial t} = a_{2} \frac{\partial^{2} T_{2}}{\partial x^{2}} - \frac{(\rho v c_{p})_{g}}{c_{2} \rho_{2}} \frac{\partial T_{2}}{\partial x}, \quad t > 0;
   \end{equation}\]
   \[\begin{equation}
   \tag{2}
   \!\!\!\!\!\!\!\!\!\!
   \lambda_{2} \frac{\partial T_{2}(x,t)}{\partial x}\Bigr|_{x=0} = q_{w} \exp \Bigl(-(\rho v c_{p})_{g}^{2} \frac{t}{4a_{2}} \Bigr) - (\rho v c_{p})_{g} T_{2}(0,t), \quad t > 0;
   \end{equation}\]
   \[\begin{equation}
   \tag{3}
   T_{2}(x_{fe}, t) = T_{fe} , \quad t > 0.
   \end{equation}\]
2. Задача теплопроводности в области (фаза «1») для определения распределения температур $T_{1}(x, t)$.
   Область определения $x_{fb} < x < \infty$:
   \[\begin{equation}
   \tag{4}
   \frac{\partial T_{1}}{\partial t} = a_{1} \frac{\partial^{2} T_{1}}{\partial x^{2}}, \quad t > 0;
   \end{equation}\]
   \[\begin{equation}
   \tag{5}
   T_{1}(x_{fb}, t) = T_{fb}, \quad t > 0;
   \end{equation}\]
   \[\begin{equation}
   \tag{6}
   \lim\limits_{x \to \infty} T_{1}(x, t) = T_{0}, \quad t > 0;
   \end{equation}\]
   \[\begin{equation}
   \tag{7}
   T_{1}(x, 0) = T_{0}. 
   \end{equation}\]
3. Задача определения квазистационарного распределения температур $T(x)$ в зоне разложения связующего КМ (фаза «3») с учетом фильтрации и эндотермических реакций с тепловым эффектом $Q^{*}$.
   Область определения $x_{fe} < x < x_{fb}$:
   \[\begin{equation}
   \tag{8}
   \frac{\partial^{2} T}{\partial x^{2}} + \frac{\rho c_{p} \dot{x}}{\lambda} \frac{\partial T}{\partial x} = \frac{\dot{\rho} Q^{*}}{\lambda};
   \end{equation}\]
   \[\begin{equation}
   \tag{9}
   T(x_{fe}) = T_{fe};
   \quad 
   T(x_{fb}) = T_{fb}.
   \end{equation}\]
   Здесь $\dot{\rho}(t)$ можно принять пропорциональной линейной скорости $\dot{x}(t)$ движения зоны разложения:
   \[\begin{equation*}
   \dot{\rho}(t) = \frac{\rho_{b} - \rho_{e}}{x_{fb} - x_{fe}} \cdot \dot{x}(t).
   \end{equation*}\]
   Изменение плотности $\rho(x(t))$ КМ в зоне разложения описывается обыкновенным дифференциальным уравнением первого порядка [15]:
   \[\begin{equation}
   \tag{10}
   \frac{d\rho(x(t))}{dt} = A \rho(x(t)) \frac{E}{\bar{R}T(x)}, \quad x_{fe}(t) < x < x_{fb}(t).
   \end{equation}\]
   Условие Стефана на подвижных границах $x = x_{fe}$ и $x = x_{fb}$:
   \[\begin{equation}
   \tag{11}
   \lambda_{2} \frac{\partial T_{2}(x_{fe}, t)}{\partial x} - \lambda_{1} \frac{\partial T_{1}(x_{fb}, t)}{\partial x} = \dot{m} Q^{*}, \quad \dot{m} = \rho_{1} \cdot \dot{x}.
   \end{equation}\]
4. Задача неизотермической фильтрации газов в пористом коксовом остатке включает следующие соотношения:

• стационарное уравнение неразрывности
  \[\begin{equation}
  \tag{12}
  \frac{d(\rho v)_{g}}{dx} = 0, \quad 0 < x < x_{fe};
  \end{equation}\]
• закон нелинейной фильтрации [1, 2]
  \[\begin{equation}
  \tag{13}
  v_{g}(x) = \frac{k}{\mu_{g} (1 + \Pi \cdot \mathsf{Re}_{g})} \frac{dp_{g}(x)}{dx};
  \end{equation}\]
• уравнение состояния
  \[\begin{equation}
  \tag{14}
  p_{g}(x) = \rho_{g} \bar{R} \cdot T_{2}(x).
  \end{equation}\]

Объединяя выражения (12)–(14), приходим к нелинейному квазистационарному уравнению фильтрации
\[\begin{equation}
\tag{15}
\frac{d}{dx} \left[ \frac{p_{g}}{\bar{R}T_{2}(x)} \cdot \frac{k}{\mu_{g}(1 + \Pi \cdot \mathsf{Re}_{g})} \frac{dp_{g}}{dx} \right] = 0, \quad 0 < x < x_{fe}
\end{equation}\]
с граничными условиями
\[\begin{equation}
\tag{16}
p_{g}(0) = p_{w} , \quad p_{g}(x_{fe}) = p_{0},
\end{equation}\]
где $k$ — коэффициент проницаемости, $\mu_{g}$ — динамическая вязкость газа, $\Pi$ — пористость коксового остатка, $\mathsf{Re}_{g} =  {\rho_{g} v_{g} d}/{\mu_{g}}$ — фильтрационное число Рейнольдса, $d$ — усредненный диаметр капилляра. Решение задачи (15), (16) дает распределение давления фильтрации $p_{g}$ в пористой области; из (13) определяется распределение скорости $v_{g}(x)$ фильтрации, а из (14) — плотности $\rho_{g}(x)$.

В соотношениях (1)–(10) приняты следующие обозначения: $c_{p}$ — теплоемкость при постоянном давлении; $x$, $t$ — пространственная переменная и время; $v_{g}$, $(\rho v)_{g}$ — линейная и массовая скорости фильтрации; $a = \lambda / (c\rho) $ — температуропроводность; $\lambda$, $c$, $\rho $ — теплопроводность, теплоемкость и плотность соответственно; $Q^{*}$ — эндотермический эффект разложения связующего КМ; $q_{w}$ — плотность теплового потока на границе $x = 0$; $E$, $\bar{R}$ — энергия активации и усредненная газовая постоянная; $A$ — размерный множитель. Индексы: $w$ — наружная граница; $b$, $e$ — начало и конец разложения связующего соответственно.

2. Метод решения

Решение комплексной проблемы (1)–(16) осуществляется в следующей последовательности.

1. Фиксируются текущее время $t=t_{f}$ и границы зоны разложения $x=x_{fe}$ и $x=x_{fb}$ с температурами $T_{e}$ и $T_{b}$ соответственно, а также времена $t_{e}$, $t_{b}$ достижимости этих границ.

2. Определяется распределение плотности КМ в зоне разложения путем интегрирования уравнения (10). Получено выражение
\[\begin{equation}
\tag{17}
\rho(x(t))=C \exp[(B/T(x))t], \quad B=A E/\bar{R},
\end{equation}\]
в котором коэффициент $B$ и постоянная интегрирования $C$ определяются из условий
\[\begin{equation}
\tag{18}
\begin{array}{ll}
x=x_{fe} =t_{e} \dot{x}: & C \exp[(B/T_{e})t_{e}]=\rho _{e}; \\
x=x_{fb} =t_{b} \dot{x}: & C \exp[(B/T_{b})t_{b}]=\rho _{b},
\end{array}
\end{equation}\]
где $\dot{x}$ — линейная скорость движения зоны разложения.

Объединяя (17) и (18) и учитывая равенства $t_{e} =x_{fe} /\dot{x}$ и $t_{b} =x_{fb} /\dot{x}$, получим искомый закон распределения плотности КМ в зоне разложения:
\[\begin{equation}
\tag{19}
\rho(x(t))=\rho _{b} \bigl( {\rho _{e} }/{\rho _{b} } \bigr)^{T_{e} \big(\frac{T_{b} }{T(x) } x(t)-x_{fb} \big)\big/(T_{b} x_{fe} -T_{e} x_{fb})}.
\end{equation}\]

Одновременно определяются распределение плотности $\rho _{g}(x)$ газов в зоне разложения, а также плотность $\rho _{g0}$ и давление $p_{g0}$ торможения:
\[\begin{equation}
\tag{20}
\rho _{g}(x)=\rho _{b} -\rho(x(t)), 
\;\; 
\rho _{g0} =\frac{1}{x_{b} -x_{t} } \int _{x_{b} }^{x_{e} }    \rho _{g} (x(t))dx, 
\;\; 
p_{g0} =\rho _{g0} \bar{R}\bar{T},
\end{equation}\]
где $\bar{T}$ — усредненная температура, например $\bar{T}=(T_{b} +T_{e})/2$.

3. Распределение температуры $T(x)$ в (19) определяется из решения задачи (8), (9):
\[\begin{multline}
\tag{21}
T(x)=T_{b} -\frac{\dot{\rho }Q^{*} }{\rho c_{p} \dot{x}}(x_{fb} -x)
\Bigl[(T_{e} -T_{b})+ {}
\\
{}+\frac{\dot{\rho }Q^{*} (x_{fb} -x_{fe})}{\rho c_{p} \dot{x}} \Bigr]
\frac{1-\exp  ( {\dot{x}}/{\bar{a}} )(x_{fb} -x)}{1-\exp ({\dot{x}}/{\bar{a}} )(x_{fb} -x_{fe})},
\end{multline}\]
где $\bar{a}$, $\rho c_{p}$ — усредненные значения в зоне разложения.

Для определения гидравлических характеристик фильтрации в пористом остатке «2» решается задача (15), (16) относительно давления фильтрации $p_{g}(x)$. Если в первом приближении распределение температур $T_{2}(x)$ в пористом остатке принять линейным
\[
T_{2}(x)=T_{f}[T_{2}(0)-T_{f}]({x_{f} -x})/{x_{f} }, \quad
x_{f} =(x_{fb} +x_{fe})/2, \quad
T_{f} =(T_{b} +T_{e})/2,
\]
что соответствует близкому распределению (см. результаты расчетов), то распределение давления фильтрации имеет квадратичный вид:
\[\begin{equation}
\tag{22}
\frac{p_{g}^{2}(x)-p_{gw}^{2} }{p_{g0}^{2}(t)-p_{gw}^{2} } =\frac{T_{2}(0,t)x-[T_{2}(0)-T_{f}]x^{2} /2}{T_{2}(0)x_{f} -[T_{2}(0)-T_{f}]x_{f}^{2} /2}.
\end{equation}\]

Это распределение давления в области «2» можно теперь использовать для нахождения линейной скорости фильтрации $v_{q}(x)$ и массовой скорости $(\rho v)_{g} = \text{const}$ из (12).

4. Решается задача Стефана (1)–(11), в которой подвижная зона заменена подвижной границей $x_{f} =(x_{fb} +x_{fe})/2$ с температурой $T_{f} =(T_{b} +T_{e})/2$, так как каждая точка зоны движется с одинаковой линейной скоростью $\dot{x}$.

При решении уравнения (1) в области «2» для исключения конвективного члена применяется подстановка [16]
\[
T_{2}(x,t)=u(x,t)\exp \Bigl(\frac{bx}{2a_{2} } -\frac{b^{2} t}{4a_{2} } \Bigr), \quad
b= \frac{(\rho vc_{p})_{g} }{c_{2} \rho _{2}},
\]
с помощью которой граничное условие третьего рода (2) трансформируется в граничное условие второго рода, а в уравнении теплопереноса относительно функции $u(x,t)$ конвективный член исчезает. В результате решением задачи (1)–(3) является функция [17]
\[\begin{multline}
\tag{23}
T_{2}(x,t)=\Bigl\{\Bigl[ T_{f} +\frac{q_{w} }{\vartheta c_{2} \rho _{2} } \operatorname{erf}\Big(\frac{x_{f} }{2\sqrt{a_{2} t_{f} } } \Bigr)\Bigr] 
\frac{1+\frac{b}{2\vartheta } \operatorname{erf}\bigl(\frac{x}{2\sqrt{a_{2} t} } \bigr)}{1+\frac{b}{2\vartheta } \operatorname{erf}\bigl(\frac{x_{f} }{2\sqrt{a_{2} t_{f} } } \bigr)} - {}
\\
{}- \frac{q_{w} }{\vartheta c_{2} \rho _{2} } \operatorname{erf}\Bigl(\frac{x}{2\sqrt{a_{2} t} } \Bigr)
\Bigr\}\exp \Bigl(\frac{bx}{2a_{2} } -\frac{b^{2} }{4a_{2} } \Bigr), \quad x\in(0,x_{f}), \; t>0,
\end{multline}\]
где $\vartheta =\sqrt{a_{2} /\pi t}$.  

Решением задачи (4)–(7) является функция
\[\begin{equation}
\tag{24}
T_{1}(x,t)=T_{0} +(T_{f} -T_{0})\frac{\operatorname{erfc}\bigl(\frac{x}{2\sqrt{a_{1} t} } \bigr)}{\operatorname{erfc}\bigl(\frac{x_{f} }{2\sqrt{a_{1} t_{f} } } \bigr)}, \quad x\in[x_{f} ,\infty], \; t>0.
\end{equation}\]

Так как аргументы функций ошибок должны быть безразмерными, естественно находить новое значение координаты подвижной границы $x_{f}$ в виде равенства
\[\begin{equation}
\tag{25}
x_{f} =2\chi \sqrt{a_{1} t_{f}},
\end{equation}\]
где коэффициент $\chi$ определяется из условия Стефана (11), в котором вместо $x_{fb}$ берется $x_{fe}$.

Подставляя (23)–(25) в (11), получим трансцендентное уравнение относительно $\chi$:
\[\begin{multline}
\tag{26}
\chi =\frac{1}{\rho _{1} Q^{*}} \biggl[
\frac{\big(q_{w} -T_{f} (c_{p} \rho v)_{g} /2\big) \sqrt{\frac{a_{2} }{\pi a_{1}}} \exp \big(- \chi ^{2} {a_{1} } / {a_{2}}  \big)}{\frac{v_{g} \bar{V}}{\sqrt{\pi}} +\frac{v_{g}}{2\bar c} \operatorname{erf}\big(\chi \sqrt{ {a_{1} } / {a_{2}}} \big)} - \\
-\frac{\lambda _{1}(T_{f} -T_{0})}{\sqrt{\pi}} \frac{\exp(-\chi ^{2})}{1-\operatorname{erf}(\chi)} \biggr],
\end{multline}\]
где
\[
\bar{V}=\frac{\sqrt{a_{2} /t}}{v_{g}} , \qquad \bar{c}= \frac{ c_{2} \rho _{2} } {(c_{p} \rho)_{g}}.
\]
Уравнение (26) решается численно. Полученное значение $\chi$ подставляется в (25), а новое значение $x_{f}$ — в соотношения (23), (24), в результате чего окончательно получаем решение комплексной проблемы Стефана:
\[\begin{multline}
\tag{27}
T_{2}(x,t)=\biggl\{ 
\frac{T_{f} +\frac{q_{w} }{\vartheta c_{2} \rho _{2}} \operatorname{erf}\big(\chi \sqrt{ {a_{1} }/{a_{2}}} \big)}{1+\frac{b}{2\vartheta} \operatorname{erf}\big(\chi \sqrt{ {a_{1} } / {a_{2}}} \big)} 
\Bigl[1+\frac{b}{2\vartheta} \operatorname{erf}\Bigl(\frac{x}{2\sqrt{a_{2} t}} \Bigr)\Bigr]- {} \hspace{2.5cm} 
\\
\hspace{1.5cm} -\frac{q_{w} }{\vartheta c_{2} \rho _{2}} \operatorname{erf}\Bigl(\frac{x}{2\sqrt{a_{2} t}} \Bigr)\biggr\}
\exp \Bigl(\frac{bx}{2a_{2}} -\frac{b^{2} t}{4a_{2}} \Bigr), \quad x\in[0; x_{f}], \; t>0, 
\end{multline}\]
\[\begin{equation}
\tag{28}
T_{1}(x,t)=T_{0} +(T_{f} -T_{0})\frac{\operatorname{erfc}\big(\frac{x}{2\sqrt{a_{2} t}} \big)}{\operatorname{erfc}(\chi)}, \quad x\in[x_{f} ;\infty], \; t>0.
\end{equation}\] 

5. Определяются новые значения границ зоны разложения $x_{e}$ и $x_{b}$ путем интерполяции распределений температур (27) по температуре $T_{e}$ и (28) — по температуре $T_{b}$ соответственно: $T_{2}(x_{e}, t)=T_{e}$ и $T_{1}(x_{b}, t)=T_{b}$. С этими значениями границ зоны разложения весь алгоритм метода решения повторяется с п. 1 при новом значении времени $t$.

3. Результаты и их анализ

На рис. 1–4 представлены результаты расчетов по формулам (19), (21), (27), (28) и (22). Входные данные принимали следующие значения: $q_{w} =10^{4}$ Вт/м$^2$, $Q^{*} =800$ Дж/кг, $T_{b} =500$ K, $T_{e} =1000$ K, $\rho _{b} =1400$ кг/м$^3$, $\rho _{e} =1300$ кг/м$^3$, $a=10^{-7}$ м$^2$/c, $\lambda _{1} =0.1$ Вт/(м${}\cdot{}$K), $\lambda _{2} =0.3$ Вт/(м${}\cdot{}$K). Для некоторого фиксированного момента времени $t_{f}$: $x_{b} =0.004$ м, $x_{e} =0.002$ м, $\dot{x}=0.00012$ м/с.

На рис. 1 показано стационарное распределение температуры в зоне разложения связующего КМ, рассчитанное по формуле (21) и используемое для определения распределения плотности КМ в (19). Температурный профиль имеет выраженную вогнутость (значительную величину второй производной $\partial ^{2} T/\partial x^{2}$), что обусловлено поглощением теплоты эндотермическими реакциями разложения $Q^{*}$ и конвективным переносом (см. уравнение (8)).

Рис. 1. Стационарное распределение температуры в зоне разложения связующего композиционного материала
[Figure 1. Steady-state temperature distribution in the binder decomposition zone of the composite material]

На рис. 2 приведены результаты расчетов плотности КМ $\rho(x)$ в зоне разложения по формуле (19) и распределение плотности генерированных газов $\rho _{g}(x)$ в соответствии с выражением (20). Эти зависимости имеют точки перегиба, в которых вторые производные меняют знак: с положительного на отрицательный для $\rho(x)$ и с отрицательного на положительный для $\rho _{g}(x)$. При этом изменение плотностей на границах $x_{b}$ и $x_{e}$ происходит скачкообразно, а не плавно: 
\[
\frac{d\rho }{dx}\Big|_{\begin{subarray}{l} {x=x_{b} } \\ {x=x_{e} } \end{subarray}} \ne 0 
\quad \text{и} 
\quad 
\frac{d\rho _{g}} {dx}\Big|_{\begin{subarray}{l} {x=x_{b} } \\ {x=x_{e} } \end{subarray}} \ne 0.
\] 
Из рисунка также видно, что ввиду малой абсолютной величины кривизны используемые в расчетах приближенно линейные распределения $\rho(x)$ и $\rho _{g}(x)$ являются обоснованными.

Рис. 2. Распределение плотности композиционного материала $\rho (x)$ и газовой фазы $\rho _{g} (x)$ в зоне разложения связующего
[Figure 2. Density distribution of the composite material $\rho (x)$ and gas phase $\rho _{g} (x)$ in the binder decomposition zone]

На рис. 3 показаны распределения температур в различные моменты времени $t=10$ c, 25 c, 50 c, 100 c и 150 c, полученные по формулам (25)–(28). На рисунке видны координаты подвижной срединной границы $x_{f} =(x_{b} +x_{e})/2$ с температурой $T_{f} =(T_{b} +T_{e})/2$, соответствующие точкам излома температурных кривых.

Принятое ранее предположение о линейном распределении температур в пористом коксогазовом остатке (кривые левее точек излома) подтверждается данными расчетами.

Рис. 3. Эволюция температурных полей $T(x,t)$ в пористом остатке и области неповрежденного материала в различные моменты времени: 1 — $t =10$ с, 2 — $t=20$ с, 3 — $t=50$ с, 4 — $t=100$ с, 5 — $t=150$ с
[Figure 3. Evolution of temperature fields $T(x,t)$ in the porous residue and virgin material region at different time instants: 1—$t=10$ s, 2—$t=20$ s, 3—$t=50$ s, 4—$t=100$ s, 5—$t=150$ s]

На рис. 4 приведены результаты расчетов давления фильтрации $p_{g}(x,t)$ в пористом коксовом остатке в различные моменты времени $t_{1} =100$ c, $t_{2} =200$ c, $t_{3} =300$ c, рассчитанные по формуле (22). Кривые имеют параболический вид с выпуклостью вверх ($\partial ^{2} p_{g}(x,t)/\partial x^{2} <0$) и ограничены подвижной границей $(x_{b} +x_{e})/2$ фазовых превращений. В этих точках давление газов принималось равным давлению торможения в зоне разложения, а на левой границе — атмосферному давлению $p_{g}(0, t)=10^{5}$ Па.

Рис. 4. Эволюция распределения давления фильтрации $p_g(x,t)$ в пористом коксовом остатке во времени: 1 — $t=100$ с, 2 — $t =200$ с, 3 — $t=300$ с
[Figure 4. Evolution of filtration pressure $p_g(x,t)$ distribution in the porous coke residue over time: 1—$t=100$ s, 2—$t=200$ s, 3—$t=300$ s]

Выводы

На основе представленных результатов можно сделать следующие выводы.

1. Сформулирована комплексная физико-математическая модель тепломассопереноса в теплозащитных композиционных материалах (КМ) в условиях фазовых превращений связующих. В рамках модели получены аналитические решения, основанные на приближении задач типа Стефана.
2. Идентифицирован закон разложения связующих КМ, использующий экспериментальные данные по температурам и плотностям начала и окончания разложения. Предложенный подход позволяет избежать трудностей, связанных с формализацией химической кинетики, что делает его применимым для решения задач тепломассопереноса в широком классе теплозащитных КМ.
3. Проведенный анализ демонстрирует существенно нелинейный характер распределения температур в зоне разложения, что указывает на некорректность использовавшихся ранее предположений о срединной температуре. Установлено, что распределения плотности КМ и газов, несмотря на наличие точек перегиба, характеризуются незначительной кривизной, обуславливающей их близость к линейным зависимостям. Также показано, что изменение плотностей на границах зоны разложения происходит скачкообразно.
4. Распределение температур в пористом коксовом остатке близко к линейному, а распределение давления фильтрации в этой области имеет параболический вид и ограничено нестационарно подвижной границей фазовых превращений. 

Конкурирующие интересы. Автор заявляет об отсутствии конфликта интересов в отношении данной публикации.
Авторский вклад и ответственность. Автор несет полную ответственность за содержание и подготовку окончательной версии рукописи. Окончательная версия рукописи одобрена автором.
Финансирование. Работа выполнена при финансовой поддержке Российского научного фонда (грант   № 22-19-00420-П, https://rscf.ru/project/22-19-00420/), выданного Московскому авиационному институту.

About the authors

Vladimir F. Formalev

Author for correspondence.
Email: formalev38@yandex.ru
ORCID iD: 0000-0003-2135-0926
SPIN-code: 4502-3688
Scopus Author ID: 6602473764
ResearcherId: T-1483-2018
https://www.mathnet.ru/rus/person30958

Dr. Phys. & Math. Sci., Professor; Professor; Dept. of Computational Mathematics and Programming

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Figure 1. Steady-state temperature distribution in the binder decomposition zone of the composite material

Download (44KB)

Indexing metadata

3. Figure 2. Density distribution of the composite material $\rho (x)$ and gas phase $\rho _{g} (x)$ in the binder decomposition zone

Download (74KB)

Indexing metadata

4. Figure 3. Evolution of temperature fields $T(x,t)$ in the porous residue and virgin material region at different time instants: 1—$t=10$ s, 2—$t=20$ s, 3—$t=50$ s, 4—$t=100$ s, 5—$t=150$ s

Download (74KB)

Indexing metadata

5. Figure 4. Evolution of filtration pressure $p_g(x,t)$ distribution in the porous coke residue over time: 1—$t=100$ s, 2—$t=200$ s, 3—$t=300$ s

Download (72KB)

Indexing metadata