CN106508040B - 一种多孔热解材料传热过程计算方法 - Google Patents

Abstract

一种多孔热解材料传热过程计算方法，包括推导复合材料内部各层质量、动量及能量守恒方程，实现多孔介质结构内部固相与气相的质量、动量及能量方程的求解，得到固相温度T、密度ρ、导热系数k_s、气相密度ρ_f、导热系数k_f、孔隙率φ、热解气体压力p、分解气体溢出速度u_n。本发明能对不同种类的碳化热解材料传热过程进行分析，能三维显示材料内部热解气体热解过程，得到热解气体压力、速度分布，以及热解后形成的多孔介质结构温度分布云图。通过计算多孔热解材料传热过程，解决了高超声速飞行器复合材料温度及舱内环境温度难于准确计算的问题。

Description

一种多孔热解材料传热过程计算方法

技术领域

本发明涉及一种多孔热解材料传热过程计算方法，特别是一种基于fluent的复合材料多孔介质结构传热过程计算方法。

背景技术

在火箭发动机内部热防护或高超声速飞行器外部气动热防护中，硅基复合材料或碳基复合材料是常用的热防护材料，准确预示这类绝热材料烧蚀过程中的温度分布，对于发动机及飞行器结构设计意义重大。其制作工艺有多种，可以将SIO2玻璃纤维浸渍酚醛树脂缠绕固化，或用玻璃纤维织成的布浸渍酚醛树脂叠层固化，或用乱玻璃纤维调入酚醛树脂后模压成形。其烧蚀过程分为四个阶段，如果材料为含有玻璃纤维等的硅基复合材料，其烧蚀过程中除了有碳化层、热解层和基体层外，表面还存在一定的熔融物。这类材料有一个共同点，即在热解碳化过程中，材料内部释放出热解气体，结构逐渐呈现出多孔介质结构。其烧蚀物理过程具体如下：

(1)发动机点火开始后，壁面温度不高，尚未达到材料的热解温度，材料机械强度较大，燃气与粒子流对绝热材料表面的烧蚀作用较弱，除了微弱的化学烧蚀外，在绝热层碳化之前，没有机械侵蚀发生。在绝热层内只有热传导发生，且绝热层内只存在基体层。

(2)随着烧蚀表面温度的不断升高，材料中的高聚物逐渐开始分解，释放出热解气体，材料密度随之降低。由于材料热解过程中仍具有一定的结构强度，绝热层表面的退移量很小。绝热层内存在热解层和基体层，热解气体的产生使得传导和对流同时存在；

(3)随着热解进行完毕，表面温度持续升高，材料内逐渐出现碳化物质。本阶段，由内向外分别存在基体层、热解层和碳化层。碳化层物质结构强度明显降低，在碳化层内，孔隙率较大，外界氧化性气体扩散到碳化层内与碳反应，消耗碳化层质量，使得孔隙率逐渐变大，碳化层密度逐渐降低，当密度低于临界值时，由于此部分碳化物质强度太低，将被燃气冲刷烧蚀掉，表面开始产生退移。热解层向内部退移，烧蚀表面也在不断向内推进。由于颗粒撞击，机械侵蚀也发生作用，加速了表面退移率，但是颗粒在表面的沉积对碳化层也起到了一定的隔热保护作用。绝热层内传导和对流同时存在。

由于材料内部受热后存在热解、碳化等过程，内部传热方式比较复杂，准确分析材料内部温度分布及热解气体压力分布具有一定难度。传统的计算这类材料导热过程的方法，未考虑材料内部热解过程吸热，常常导致计算得到的温度过高，误差较大；或者是考虑了热解吸热，并将材料热解碳化后的多孔特性通过孔隙率来体现，但是却未考虑热解气体及由此产生的气体分压对传热过程的影响，且计算多是针对一维或二维结构。通过试验测量材料内部温度来拟合出温度与导热系数的经验关系式，也是一种常用的方法，但是这类方法所得到的经验公式仅对特定的材料有效，对于新材料而言，需要通过试验重新建立关系式，通用性不强。

为掌握碳化热解材料内部传热过程计算方法，获得材料内部温度分布，本发明基于多孔介质理论，推导了材料碳化热解过程中质量、动量及能量守恒关系式，基于fluent软件二次开发了UDF程序，对碳化热解过程中的传热传质过程进行求解，得到了材料内部三维温度分布及热解气体速度、压力分布。

发明内容

为了解决现有热解材料传热过程计算无法得到材料内部三维温度分布及热解气体速度、压力分布的问题，本发明提出了一种多孔热解材料传热过程计算方法。

本发明通过以下技术方案来完成发明任务：

第一步，建立碳化层质量、动量、能量守恒方程：

公式(1)中，ρ_f为热解气体密度，为待求量；ρ_s为固相密度，待求量；φ为孔隙率，为待求量；V为热解气体速度，待求量。

其中，ρ为材料总密度，待求量；c_p为材料总比热容，待求量；k为材料总导热系数，待求量；φ为孔隙率，待求量；ρ_s为固相密度，待求量；ρ_f为热解气体密度，待求量；c_ps为固相比热容，待求量；c_pf为热解气体比热容，待求量；k_f为热解气体导热系数，待求量；k_s为固相导热系数，待求量；为材料烧蚀率，已知量；Q_reaction为碳化材料反应热，已知量。

在材料表面热解气体的溢出速度u＝u₀＝const。

第二步，建立热解层质量、动量、能量守恒方程：

其中u_i、u_j为热解气体流速，待求量；T为材料温度，待求量；T₀为材料初始温度，已知量；p为热解气体压力，待求量；μ为粘性系数，已知量；k为材料渗透率，已知量；x_j为坐标变量，已知量；Q_s为热解吸热量，已知量。

第三步，求解碳化层及热解层材料密度、比热容、导热系数及孔隙率，计算公式如下：

ρ＝(1-φ)ρ_s+φρ_f (7)

k＝φk_f+(1-φ)k_s (9)

ρ_s由确定；

第四步，计算碳化层及热解层内热解气体与固相之间热交换能量源项，表达式如下：

碳化层：

热解层：

第五步，计算热解层与碳化层内热解气体溢出速度与压力，计算公式如下：

p＝φρ_fRT (14)

n为法向，n_i为热解指数，已知量；R为通用气体常数，已知量；u_n为热解气体速度，待求量；A为计算单元表面积，已知量；E为能量系数，已知量。

第六步，设置整个流场内固相温度T、密度ρ、导热系数k_s、气相密度ρ_f、导热系数k_f、孔隙率φ、热解气体压力p初值，设置材料表面加载热流边界条件及其他各表面温度边界条件，根据公式(1)～(14)进行计算，求解固相温度T、密度ρ、导热系数k_s、气相密度ρ_f、导热系数k_f、孔隙率φ、热解气体压力p、分解气体溢出速度u_n。

第七步，将计算结果进行三维显示，得到材料内部温度、热解气体压力、速度分布云图。

通过以上计算，可以得到热解材料多孔介质结构内部三维温度、热解气体压力、速度三维分布云图。为了验证本发明计算方法的正确性，对高硅氧酚醛复合材料加热过程进行了计算，计算得到的温度分布如图1所示，热解气体溢出速度分布如图2所示，材料内气体压力分布如图3所示，材料孔隙率分布如图4所示，计算得到的冷壁面温升曲线与试验测量结果对比如图5所示，由此可见，本发明所建立的计算方法可正确模拟复合材料多孔介质结构传热过程。解决了高超声速飞行器复合材料温度及舱内环境温度难于准确计算的问题。

附图说明

图1热解材料温度分布。

图2多孔介质材料内部热解气体溢出速率分布。

图3多孔介质材料内部热解气体压力分布。

图4多孔介质材料孔隙率分布。

图5计算结果与试验结果对比。

具体实施方式：

推导碳化层、热解层及基体层内气相及固相的质量、动量及能量守恒方程：

根据雷诺输运定理，物理量的守恒定理可写成：

公式(15)中，为待求解物理量，V为热解气体速度，t为时间，为源项。

(1)碳化层：

碳化层内除了存在有碳化物残渣，还有分解产生的热解气体流过，属于多孔介质。因此，在分析其传热传质过程时，要分固体和分解气体两部分进行。

在碳化层内任取一微元控制体，控体制内既包含固体，也包含分解气体。

质量守恒方程：

对于碳化层中的固相部分，虽然热解过程已进行完毕，不再有可分解成份存在，但是，碳化物质与扩散过来的氧化气体反应，可使碳化层质量减少，密度降低，使得孔隙率变大，因此有：

公式(16)中，S_m为质量源项，ρ_s为固相密度，为待求量，为碳化物质烧蚀速率，为已知量。

因此，整个控制体内的质量守恒方程可写为：

公式(1)中，ρ_f为热解气体密度，为待求量，φ为孔隙率，为待求量。

动量守恒方程：

标准流体的动量守恒方程可写为：

公式(17)中p为热解气体压力，为待求量，τ为热解气体剪切力，g为重力加速度，为已知量。

在公式(17)中，粘性项变为由于流速很小，很小，且认为流速分布均匀，则多孔介质中流体相的动量守恒方程可写为：

能量守恒方程：

对于固体部分，控制体总能量变化率等于单位时间内外界加入的热量与内部热源所增加的热量之和。固相总能量变化率可写为：

公式(18)中，e为材料内能，为待求量，Ω为微元控制体。

公式(18)中所以

公式(19)中c_ps为固相材料比热容，为已知量，ΔT为温差。

单位时间内外界加入的热量为：

公式(20)中k_s为固相导热系数，为待求量。此项为由于热传导进入固体内的热量。

碳化层内不存在热分解过程，没有热解吸热；但存在化学反应热。因此，内部热源产生的热增量其中，q_s为化学反应产生的热流率，Q_reaction为化学反应热。因此，固体碳化层内的能量守恒方程可写为：

微分形式为：

忽略定压比热容随时间的变化，并认为k_s与位置无关，则公式(22)可写为：

此即为固相能量守恒方程的最终表达式。

对于分解气体所属的流体相，其控制体内总能量的变化率为：

单位时间内外界加入的热量有体积力以及表面力所做的功W、控制体外流体通过体系表面对控制体内流体导热Q_c以及辐射传热Q_r。

其中，

公式(25)中，A为微元表面积，σ为表面力，n为微元表面法向。

由于碳化层内温度相对不高，因此忽略辐射影响。

传导传热量

公式(26)中k_f为热解气体导热系数，为待求量。

因此，流体相能量守恒方程可写为：

其微分形式为：

式(28)左端第一项可化为：

公式(29)中，c_pf为热解气体比热容，为待求量。

将公式(1)代入公式(29)，上式可化为：

式(28)中左端第二项可化为：

将公式(30)、(31)代入式(28)中，可得

对于压力做功可如下处理：

由公式(2)，可得

则：

所以，公式(32)可化为：

由于烧蚀开始时的非稳态阶段时间较短，流速较低，且气体压缩性不大，因此忽略以及项的影响，所以有：

即：

如果忽略粘性力做功，则公式(37)可写成：

考虑到多孔介质孔隙率影响以及连续性方程，公式(38)的最终表达式为：

综上，碳化层中整个控制体内的能量守恒方程为：

即：

(2)热解层：

热解层各相关控制方程推导过程与碳化层类似，这里不作详细说明，仅给出具体表达式如下：

质量守恒方程：

由上面的分析，由于热解层内存在热解过程，即源项并且假设分解气体质量流均匀溢出，因此，质量守恒方程为：

热解层动量守恒方程参见公式(5)。

能量守恒方程：

分析过程同前，能量守恒方程基本形式同式，所不同的是流体与固体的密度随时间的变化率不再为零，并且热解过程要吸收大量的热量，因此源项q_h不为零，

因此固相能量守恒方程变为：

流体相能量守恒方程变为：

即：

考虑到孔隙率影响以及连续性方程，热解层控制体内能量守恒方程的最终表达式为：

(3)基体层：

能量守恒方程：

碳化层与热解层内材料密度、比热容及导热系数公式如下：

为确定分解气体溢出速率，在基体层与分解层交界面处取一微元控制体，认为此微元内只有溢出气体，没有流进气体，且流速沿法向流向内表面。则有如下方程组：

p＝φρ_fRT (14)

联立式(12)～(14)进行迭代求解，即可得到分解气体溢出速度u_n

设置整个流场内固相温度T、密度ρ、导热系数k_s、气相密度ρ_f、导热系数k_f、孔隙率φ、热解气体压力p初值，设置材料表面加载热流边界条件及其他各表面温度边界条件，根据公式(1)～(14)进行计算，求解固相温度T、密度ρ、导热系数k_s、气相密度ρ_f、导热系数k_f、孔隙率φ、热解气体压力p、分解气体溢出速度u_n。

Claims (1)

1.一种多孔热解材料传热过程计算方法，其特征是包括以下步骤：

第一步，建立碳化层质量、动量、能量守恒方程：

$\frac{\∂(1-\mathrm{\φ}){\mathrm{\ρ}}_s}{\∂t}+\frac{\∂{\mathrm{\φ\ρ}}_f}{\∂t}+\▿\·\left({\mathrm{\ρ}}_{f}V\right)=0---\left(1\right)$

公式(1)中，ρ_f为热解气体密度，为待求量；ρ_s为固相密度，待求量；φ为孔隙率，为待求量；V为热解气体速度，待求量；

$\frac{\∂\left({\mathrm{\ρ}}_{f}V\right)}{\∂t}=-\▿p+\mathrm{\ρ}g-\frac{\mathrm{\μ}}{k}V---\left(2\right)$

$\begin{array}{c}\frac{\∂(1-\mathrm{\φ}){\mathrm{\ρ}}_s}{\∂t}\left(c_{ps}\right(T-T_0)-c_{pf}(T-T_0)-Q_{reaction})+\frac{\∂T}{\∂t}\left(\right(1-\mathrm{\φ}){\mathrm{\ρ}}_s{c}_{ps}+{\mathrm{\φ\ρ}}_f{c}_{pf})\\ =({\mathrm{\φk}}_f+(1-\mathrm{\φ}\left)k_s\right)\frac{{\∂}^{2}T}{\∂x_j^2}-{\mathrm{\ρ}}_f{c}_{pf}{u}_j\frac{\∂T}{\∂x_j}+\frac{\mathrm{\μ}}{k}{u}_j{u}_j-p\frac{\∂u_i}{\∂x_j}\end{array}---\left(3\right)$

其中，ρ为材料总密度，待求量；c_p为材料总比热容，待求量；k为材料总导热系数，待求量；φ为孔隙率，待求量；ρ_s为固相密度，待求量；ρ_f为热解气体密度，待求量；c_ps为固相比热容，待求量；c_pf为热解气体比热容，待求量；k_f为热解气体导热系数，待求量；k_s为固相导热系数，待求量；Q_reaction为碳化材料反应热，已知量；

在材料表面热解气体的溢出速度u＝u₀＝const；

第二步，建立热解层质量、动量、能量守恒方程：

$\frac{\∂(1-\mathrm{\φ}){\mathrm{\ρ}}_s}{\∂t}+\frac{\∂{\mathrm{\φ\ρ}}_f}{\∂t}+\▿\·\left({\mathrm{\ρ}}_{f}V\right)=0---\left(4\right)$

$\frac{\∂\left({\mathrm{\ρ}}_{f}V\right)}{\∂t}=-\▿p+\mathrm{\ρ}g-\frac{\mathrm{\μ}}{k}V---\left(5\right)$

$\begin{array}{c}\frac{\∂(1-\mathrm{\φ}){\mathrm{\ρ}}_s}{\∂t}\left(c_{ps}\right(T-T_0)-c_{pf}(T-T_0)-Q_s)+\frac{\∂T}{\∂t}\left(\right(1-\mathrm{\φ}){\mathrm{\ρ}}_s{c}_{ps}+{\mathrm{\φ\ρ}}_f{c}_{pf})\\ =({\mathrm{\φk}}_f+(1-\mathrm{\φ}\left)k_s\right)\frac{{\∂}^{2}T}{\∂x_j^2}-{\mathrm{\ρ}}_f{c}_{pf}{u}_j\frac{\∂T}{\∂x_j}+\frac{\mathrm{\μ}}{k}{u}_j{u}_j-p\frac{\∂u_i}{\∂x_j}\end{array}---\left(6\right)$

其中u_i、u_j为热解气体流速，待求量；T为材料温度，待求量；T₀为材料初始温度，已知量；p为热解气体压力，待求量；μ为粘性系数，已知量；k为材料渗透率，已知量；x_j为坐标变量，已知量；Q_s为热解吸热量，已知量；

第三步，求解碳化层及热解层材料密度、比热容、导热系数及孔隙率，计算公式如下：

ρ＝(1-φ)ρ_s+φρ_f (7)

$c_p=\frac{(1-\mathrm{\φ}){\mathrm{\ρ}}_s{c}_{ps}+{\mathrm{\φ\ρ}}_f{c}_{pf}}{(1-\mathrm{\φ}){\mathrm{\ρ}}_s+{\mathrm{\φ\ρ}}_f}---\left(8\right)$

k＝φk_f+(1-φ)k_s (9)

ρ_s由确定，为材料烧蚀率，为已知量；

第四步，计算碳化层及热解层内热解气体与固相之间热交换能量源项，表达式如下：

碳化层：

热解层：

第五步，计算热解层与碳化层内热解气体溢出速度与压力，计算公式如下：

$\frac{\∂p}{\∂n}=-\frac{\mathrm{\μ}}{k}{u}_n---\left(12\right)$

$\frac{\∂(1-\mathrm{\φ}){\mathrm{\ρ}}_s}{\∂t}=-{\mathrm{\ρ}}_{0}A{\left(\frac{{\mathrm{\ρ}}_s-{\mathrm{\ρ}}_f}{{\mathrm{\ρ}}_0}\right)}^{n_1}\exp (-\frac{E}{RT})=-\frac{\∂{\mathrm{\φ\ρ}}_f}{\∂t}-\▿\·{\mathrm{\φ\ρ}}_f{u}_n---\left(13\right)$

p＝φρ_fRT (14)

n为法向，n₁为热解指数，已知量；R为通用气体常数，已知量；u_n为热解气体速度，待求量；A为计算单元表面积，已知量；E为能量系数，已知量；

第六步，设置整个流场内材料温度T、密度ρ、导热系数k_s、热解气体密度ρ_f、导热系数k_f、孔隙率φ、热解气体压力p初值，设置材料表面加载热流边界条件及其他各表面温度边界条件，根据公式(1)～(14)进行计算，求解材料温度T、密度ρ、导热系数k_s、热解气体密度ρ_f、导热系数k_f、孔隙率φ、热解气体压力p、热解气体速度u_n；

第七步，将计算结果进行三维显示，得到材料内部温度、热解气体压力、速度分布云图。