高超声速飞行器热流密度/分层温度/碳化层研究

2018-07-23，，，

中国空间科学技术 2018年3期

，，，

北京航天长征飞行器研究所高超声速飞行器防隔热中心，北京 100076

高超声速飞行器在临近空间长时间、高机动飞行时，高温真实气体效应、化学非平衡效应、激波/边界层干扰，以及转捩等特征明显，给气动热环境预示带来一系列复杂问题,特别是更难以准确预测边界层转捩。地面风洞试验目前仅能做到小于Ma12的模拟，因此，直接测量高超声速飞行器表面热流密度这一关键参数，对摸清真实飞行条件下高超声速飞行器热环境和飞行器热防护设计都有着极为重要的意义。国外一些研究者在该方面开展了许多研究[1-5]。典型的是德国SHEFEX-Ⅱ测得Ma9.3飞行热流数据，受传感器使用温度不能超过700℃的限制，仅测得飞行器尾部热流。国内在地面风洞试验技术和热流测量技术方面开展了许多研究[6-9]，获得了Ma10以下大量风洞热流测量数据。树脂基材料由于密度小、导热系数低，被广泛应用于高超声速飞行器热防护结构。在高温气动加热条件下，通过材料的烧蚀/热解/碳化及对边界层的“热阻塞”作用，实现对高超声速飞行器热防护目的。国内外学者针对树脂基材料、蜂窝材料温度场计算方法进行了大量的理论和试验研究[10-23]，但未研究树脂基复合材料导热微分方程中加入导热系数随温度变化的烧蚀温度场计算方法，也未定量给出碳化层厚度及计算结果与实测结果的比较。

本文采用在飞行器表面开孔安装长时耐高温热流传感器直接测量热流密度方法，国内首次获得Ma12以上飞行器表面临近空间热流密度真实变化历程和边界层转捩特征，验证了热环境预示方法的正确性。在考虑热解吸热的树脂基常规固体导热微分方程中加入材料导热系数随温度变化项后的导热微分方程，计算了飞行器热结构分层温度和烧蚀碳化层厚度，并与实测飞行器热结构烧蚀碳化量和分层温度进行了比较，验证了修改计算方法的合理性和正确性。

1 锥面热流密度

再入飞行器以高超声速在大气层飞行时，受气体粘性作用影响飞行器表面相对速度为零，在紧贴近飞行器表面薄层里，速度由零迅速增至来流速度，此薄层即为边界层。在边界层内气体动能部分转化为热能，存在很大的温度梯度，由此温度梯度形成了对高超声速飞行器表面的气动加热。采用参考焓方法计算有攻角再入飞行器锥面热流密度。在转捩区内热流密度用层流值和湍流值加权平均方法。

轴向X位置的锥面热流密度qw：

(1)

式中：qx为不考虑焓值影响的锥面热流，按照层流、湍流、转捩3种状态分别计算；hre为恢复焓；hw为壁焓。层流状态下：

(2)

(3)

式(2)中带“*”参考值由当地物面压力和参考焓确定，h*参考焓和hre恢复焓分别为：

h*=0.19hre+0.23he+0.58hw

(4)

(5)

湍流状态下：

当Re<107，

(6)

当Re≥107，

(7)

h*=0.22hre+0.28he+0.50hw

(8)

转捩状态下：

qxtr=(1-Γ)qxL+ΓqxT

(9)

式中：qxtr为转捩热流密度；qxL为层流热流密度；qxT为湍流热流密度；Γ为加权函数。

2 热结构烧蚀碳化热解模型

2.1 烧蚀表面控制方程

高超声速飞行器热防护通常采用烧蚀热防护方式，其烧蚀防热机理是在气动热环境作用下，依靠材料自身物理、化学变化吸收气动加热，同时烧蚀产物蒸发气体和树脂热解气体向边界层引射，起到对边界层高温气体的阻塞，从而减少了部分气动加热。常用热结构防热材料为硅基类复合材料。该类材料特点是在一定高温条件下飞行器表面形成一层熔融的液态层，其烧蚀表面控制方程包括：

(1)表面能量平衡方程

fcΔhc+αSiO2fSiO2Δhv]

(10)

式中：Ψ为引射因子；ε为材料辐射系数；σ为斯蒂芬-波尔兹曼常数；Tw为表面温度；T0为初始温度；Cp为材料比热；mt为质量烧蚀速度；fp为树脂含量；Δhp为树脂分解热；fc为碳化分数；Δhc为碳燃烧热；αSiO2为二氧化硅蒸发分数；fSiO2为二氧化硅质量分数；Δhv为二氧化硅蒸发热。左边第一项为表面气动加热，包括质量引射对气动加热的引射因子Ψ，第二项为表面向环境辐射热；右边第一项为材料初始时刻升至烧蚀温度时本身吸收热量，第二项为树脂热解吸热，第三项为碳燃烧热，第四项为SiO2蒸发吸热。

(2)质量守恒方程

(11)

式中：B为积分因子；ρ为材料密度；δL为液态层厚度；P为边界层外缘压力；μw为液态层粘度；τ为剪切力。

(3)引射因子

(12)

式中：层流下β=0.62，湍流下β=0.2；层流下η=0.26，湍流下η=0.33；Ma为空气分子量；Mj为引射气体分子量。Ψ反映了质量引射时气动热流密度与无质量引射时气动热流密度之比,质量烧蚀速度mt可通过迭代方法求得。

2.2 热结构温度响应/热解碳化多元耦合控制方程

高超声速飞行器热结构温度场通常采用一维固体内部导热微分方程计算。求解区域为除烧蚀层的热结构内部温度和承力结构内壁温度。常规固体导热微分方程求解树脂基热结构内部温度响应时不考虑材料物性参数随温度变化和树脂热解碳化过程的吸热，计算的承力结构内壁温度往往偏高。偏差主要原因是气动热通过表面传向内部时，防热材料中的树脂在高温作用下会发生热解碳化，一般300℃时开始分解，800℃左右分解碳化完毕，树脂基材料密度由100%降至60%左右。这期间材料热物理性能发生了明显变化，热物性参数已不是常温状态下的热物性参数，同时材料热解碳化过程吸收一定热量，特别是长时间飞行条件下，该作用更为明显。为合理计算高温树脂基热结构内部温度响应，应在考虑了材料热解吸热的常规导热微分方程中加入材料热物性参数随温度变化项。

考虑材料动态热解吸热、材料热物性参数随温度变化条件下，飞行器热结构内部温度响应的控制方程为：

(13)

根据修正的阿仑尼乌斯(S.Arrhenius)热动力学方程，热解气体质量生成率为：

(14)

式中：ρ0为原始材料密度；ρc为材料完全碳化的极限密度；K0为频率因子；E为活化能；R为气体常数；n为反应级数。

材料导热系数随温度变化关系为：

λ=a+bT+cT2

(15)

式中：a，b，c为通过试验回归计算所得系数。对式(13)所示的控制方程采用时间和空间有限差分离散，以式(10)烧蚀移动方程求得的Tw为边界条件，通过迭代方法即可求得考虑物性参数变化和热解吸热的飞行器热结构内部温度响应和碳化层厚度。

3 计算结果与讨论

3.1 实测表面热流与理论计算比较

图1为国内首次实测Ma12以上飞行器0°子午线(迎风母线)，X/L=0.53处(X为测点距飞行器原点距离，L为飞行器总长)热流密度真实变化历程和理论预示结果的对比。由图1可知，实测飞行器表面热流密度变化趋势与理论预示热流密度变化趋势一致，实测热流密度略高于理论计算值，结合飞行器热流传感器地面试验结果分析认为属于热流传感器系统测量误差。总的来看理论预测热流密度结果与实测值偏差较小，两者偏差小于20%。热环境工程设计允许的误差在30%。

图1 热流密度实测值和计算值比较Fig.1 Comparison between measured and calculated values of the heat flux density

3.2 理论计算转捩热流与实测比较

飞行器在临近空间飞行时，随着高度的降低，Re逐渐增大，当飞行器降低到一定高度，Re达到某个值时，飞行器表面的边界层出现转捩，流态由层流逐渐变为湍流。一直以来，边界层转捩问题在飞行器热防护设计中备受关注，因为边界层一旦发生转捩，将对气动力、气动热产生重要影响。图2给出了实测飞行器0°子午线，X/L=0.53处边界层转捩真实历程和理论计算结果的比较。由图2可知热流密度在某一时刻出现了突增，表明在该时刻边界层出现转捩。理论计算和实测热流密度的转捩起始时间点基本一致，曲线变化规律一致且实测值略高于理论值，主要原因是转捩是非常复杂的流态变化过程，准确计算转捩热流值目前还很困难，计算方法有待于完善。

图2 转捩热流实测值和计算值比较Fig.2 Comparison between the measured and calculated values of the transition heat flux density

3.3 热结构内部分层温度实测结果与理论计算比较

图3、图4是飞行器热结构X/L=0.53，0°子午线处距离防热层表面4 mm、8 mm处理论计算值与实测分层温度的比较。由图3、图4可知，理论计算防热结构分层温度响应规律与遥测结果一致，距离防热结构表面越近则温度变化速率越大，主要是防热结构表面热传递是非稳态过程，距表面越近则温度梯度越大，沿厚度方向温度梯度逐渐减小，表现温度变化率趋于减缓。从计算值和实测值对比还可知，不考虑树脂热解吸热和材料物性参数变化，理论计算值明显高于实测值，图3最大偏差275℃，图4是320℃，而考虑热解吸热和物性参数随温度的变化，则理论计算值与实测值最大偏差小于70℃，可见在树脂基常规导热微分方程中考虑树脂热解吸热和导热系数随温度变化因素更符合高温气动加热条件下树脂基防热材料物理化学变化过程，计算结果更为合理。

3.4 热结构烧蚀碳化实测结果与理论计算比较

图5为飞行器0°和45°子午线，X/L=0.53两个测点处防热层碳化层厚度实测值和理论计算结果比较。可知，理论计算和实测结果变化趋势一致，理论值略高于实测值，0°迎风子午线防热层烧蚀碳化厚度均高于45°子午线，是由于0°迎风子午线的热环境高于45°子午线，热环境高的防热层表面温度也高，表面温度越高，烧蚀碳化越严重。