杨伟斌， 朱庆勇

(中山大学航空航天学院,广东深圳 518107)

引 言

在航空航天领域,烧蚀是高超声速飞行器最重要的热防护形式.当烧蚀材料受热时,材料发生热解,释放出热解气体.烧蚀过程中,热解气体进入边界层,降低碳化层热流,起到了热阻塞的作用;同时热解气体流经碳化层引起温度的降低,起到冷却剂的作用.在这类热防护系统设计过程中,温度场的预测格外重要.因温度场预测不准而不得不采用过大的防热冗余,将严重影响飞行器性能;而不足的防热设计将有可能导致飞行器损毁.

烧蚀材料通常用作热防护材料,可以保护航天器免受大气摩擦产生的大量热量[1].针对材料的烧蚀机理以及高超声速飞行器的烧蚀流场特性,国内外展开了广泛的研究[2-3].Gokcen等用一维模型计算了碳化物在烧蚀过程中的内部热响应[4].李鸿权等对热解气在碳化层内流动进行了研究，发现热解气对热防护结构的温度分布有很大的影响[5].刘骁等研究了碳化烧蚀材料的三维温度场[6].多孔材料碳化过程中热解气体产生的机理已被许多研究者广泛研究[7-8].用质量守恒、动量守恒和能量守恒方程描述物质守恒,可以分析热解气体在模型中的压力分布和温度分布,气体在多孔介质中的流动遵循Darcy定律,朱庆勇等通过从Darcy定律导出的数学模型来探讨流体在多孔介质中的扩散,建立了毛细管内流体以质量和能量平衡形式扩散的控制方程[9-10].

考虑到热解气体的运动,气体在碳化物中的渗透率对烧蚀过程和热防护有着重要的影响,对于具有自相似性的多孔介质,可以用分形维数来评价其渗透率[11-14].分形渗透率模型与3个参数密切相关,即多孔介质的迂曲度分形维数、孔径分形维数和最大孔径.Yu等建立了渗透率模型,指出分形分析对评价多孔介质的输运性质具有重要意义,通过对输运性质的研究,得到了双重多孔介质中气体输运模型[13].Li等建立了考虑气体压缩性和滑移效应的分形致密多孔介质中的气体输运模型[15].王登科等针对气体渗流中存在的Klinkenberg效应研究了煤体瓦斯中的渗流规律[16].Cao等采用二元混合格子Boltzmann模型(lattice Boltmann method,LBM)模拟了分数Brown运动模型产生的随机分形多孔介质中的扩散过程[17].

渗透率随烧蚀过程的变化而变化,为了准确评估烧蚀过程,文章结合全局坐标和局部坐标的关系,推导出了材料的渗透率张量,并将新模型引入烧蚀模型中,得到了模型温度随烧蚀过程的变化,并将温度分布与实验数据进行了比较,验证了方法的有效性.模拟计算了不同孔径分形维数和不同迂曲度分形维数材料三维球头结构的烧蚀过程.

1 烧蚀模型

1.1 物理模型

烧蚀结构受到外部热流的影响,发生化学热解反应,并同时受扩散因素的影响.当温度较低时,烧蚀速率由化学动力学(Arrhenius方程)和扩散速度共同决定.根据壁面组元质量守恒和相容条件,表面烧蚀率可计算为[18]

(1)

热解气包含O2,O,N2,N,CO,CN,C2N,C1,C2,C3,C4,C5等13种组分.各组分之间反应如下[19-20]

(2)

(3)

其中,Mi是组分i的分子质量,Mav是气体平均分子质量.烧蚀率可由式 (1)～ (3)计算得出.

1.2 渗透率模型

图1为扫描电子显微镜下的纤维材料表面形貌.图像显示碳化纤维材料的横截面,经过图像处理运用盒子计数法可以得到如图2所示的分形维数、最大孔径和孔隙率等材料属性,从而得出热解气体在碳化层的渗透率[21].

图1 烧蚀材料的扫描电镜图像Fig.1 SEM image of ablative material

在计算多孔介质微孔中的气体扩散时,须考虑气体的滑移效应.对于气体流动,可以分为黏滞流、分子流和黏滞-分子流[22].哪种流动起主要作用,主要由Knudsen数决定.

其中,λ为微孔尺寸,l为气体平均分子自由程.平均分子自由程可以表示为

其中,μ为动力黏性系数,P为压强,M为分子质量,T为材料温度.热解气的平均分子自由程可以计算为1.643×10-8m.烧蚀材料的孔径大约为1 μm,进而计算得出Knudsen数为0.033,在计算中要考虑气体的滑移效应.

朱庆勇等综合考虑模型倾斜角度、材料分形特征、气体滑移边界条件等,推导出气体在多孔介质中的渗透率张量形式[22]为

其中

在实际工程当中,热流往往沿着一个特定的方向,因此可仅考虑气体在一个方向上的扩散,方程表示为

其中，DpT和DtL分别表示材料模截面上的孔径分形维数和长度方向的迂曲度分形维数.

(a)Tortuosity fractal dimensions

1.3 烧蚀模型[18]

当发生热解之前,热传递方式主要是材料的热传导.温度达到一定程度,材料热解为热解气和碳化结构.热解气经由碳化层的孔隙扩散到边界层.

通过上述分析,文章做出以下简化和假设[2]:

(1)在烧蚀过程中,烧蚀材料由外到内可分为碳化层-热解面-原始材料3层结构;

(2)热解反应仅发生在热解面上;

(3)热解气体与碳化层之间换热充分,二者无温差;

材料的热解过程主要发生在一定的温度范围内,一般在热解面附近.首先,材料受到外部热流的影响逐渐被加热,其中表面温度最先升高,同时由于材料自身热传导的影响,原始材料层温度也会升高.当温度达到热解温度后,材料发生热解吸热,形成热解气体和多孔碳化层.同时由于航天器使用环境使得边界层热流密度一般很大,热解层很薄,厚度可以忽略不计,近似热解面.热解气通过多孔碳化层注入热防护机构边界层带走部分热量.

图3 烧蚀结构示意图Fig.3 Schematic diagram of ablation structure

原始层的传热方程

碳化层的传热控制方程

烧蚀边界边界条件

热解面边界条件

其中，qN为热流密度;ρ表示密度,cp表示比热,κ表示导热系数,h表示汽化潜热.下标1,2分别表示原始层和碳化层

热解面上的热解气体质量流量由热解温度决定,为

其中，下标Ep表示热解面所在单元.

碳化层为多孔介质,当热解气通过碳化层流动时服从气体连续方程和Darcy定律[6]

其中，ε是孔隙率,[K]是渗透率张量,μ是热解气动力黏性系数.碳化层外边界条件为

P=Pe

其中，Pe为外壁面压强.热解面边界条件为

2 计算结果和分析

根据前面阐述的方法进行系列数值计算,将实验结果与参考文献[24]的实验结果进行对比.计算条件为:热流密度为279.7 kW/m2,材料厚度为30 mm,烧蚀时间为600 s.材料物理参数和初始条件如表1所示，材料迂曲度分形维数为1.1，孔径分形维数为1.8.

表1 材料的物理参数Table 1 Detailed parameters of calculation model and test material

图4给出不同位置的温度变化情况.对比可以看出,数值计算的温度变化与实验结果基本一致.实验结束时,近壁点温度为1 206 K，与实验温度 1 236 K 拟合良好.从中可以看出本文模型能够很好地计算地热防护机构的防护效果.

图4 烧蚀结构内部温度计算结果与实验对比Fig.4 Comparison of temperature increase at different thickness locations with test results

建立三维球头模型计算材料孔径分形维数和迂曲度分形维数不同对烧蚀的影响.球头模型的半径为1.9 cm.具体模型参见参考文献[23].

2.1 不同孔径分形维数材料的烧蚀

图5给出了不同孔径分形维数材料制成的热防护结构烧蚀70 s时的材料密度分布情况.图6显示的是不同孔径分形维数的热防护机构烧蚀70 s时的热解气密度分布情况.结合温度分布可以看出,在球头位置,温度最高,相应地在密度上反映的是材料密度减少,形成热解气,造成球头位置热解气密度较高.70 s 时,球头位置密度分别是1.719 g/cm3(Dp=1.7)，1.717 g/cm3(Dp=1.8)和1.714 g/cm3(Dp=1.9).可以看出随着孔径分形维数的增加,固体密度变小.相应地,孔径分形维数为1.7,1.8和1.9的球头结构在 70 s 时的气体密度分别是0.187，0.189，0.192 g/cm3.

(a)Dp=1.7

(a)Dp=1.7

分别计算3种孔径分形维数材料热防护层的防护效果.具体温度分布如图7所示:计算70 s时,孔径分形维数1.7,1.8和1.9的材料组成的球头前点温度分别是1 901,1 909和1 919 K,而相应地热防护层内层的温度分别是388,393和398 K.

(a)Dp=1.7

从图8可以明显看出,模拟计算的前50 s,材料背面的温度基本没有变化,相当于初始温度.50 s 时,3种不同孔径分形维数材料热防护结构的背面温度均为300.14 K,50 s仅上升0.14 K.50 s后球头背面温度快速上升,到70 s时,3种材料(孔径分形维数分别为1.7,1.8和1.9)的温度已到达388,393和398 K,分别上升了88,93和98 K.可以明显看出孔径分形维数越小,热防护效果越好.

图8 烧蚀70 s时不同孔径分形维数材料烧蚀结构背面的温度变化Fig.8 Temperature change in 70 s of ablative material with different pore area fractal dimensions

2.2 不同迂曲度分形维数材料的烧蚀

图9给出了不同迂曲度分形维数的材料制成的球头烧蚀70 s时的材料密度分布情况,图10显示的是70 s时的热解气体密度分布情况.在球头位置,材料受到热流的直接影响,温度最高,相应的材料热解程度最高,形成更多的热解气,造成球头位置热解气密度较高.

(a)Dt=1.05

(a)Dt=1.05

70 s时,3种不同材料球头位置密度分别是1.707 g/cm3(Dp=1.05),1.717 g/cm3(Dp=1.1)和1.724 g/cm3(Dp=1.15).可以看出随着迂曲度分形维数的增加,固体密度变大.相应地,3种材料的球头结构在70 s时的气体密度分别是0.200,0.189,0.181 g/cm3.

图11，12分别为不同迂曲度材料烧蚀结构温度场分布和背面的温度变化.从图11，12可以看出:球头温度分布受到烧蚀材料孔隙迂曲度分形维数的影响.当烧蚀材料的迂曲度分形维数分别是1.05,1.1和1.15时,计算70 s时球头结构的前端温度分别是1 952,1 909 和1 878 K.相应地,球头背面温度分别是418,393和378 K.

(a)Dt=1.05

图12 不同迂曲度分形维数材料烧蚀结构背面的温度变化Fig.12 Temperature change in 70 s with different tortuosity fractal dimensions

几种不同迂曲度的材料在烧蚀前期,热防护效果差别不大.在计算50 s时,球头背面温度均为300.1 K,没有差别.随着烧蚀的进行,差距逐渐出现,计算70 s结束时,三者(迂曲度分形维数分别是1.05,1.1和1.15的球头背面)的温度分别是418,393和378 K,温差达到40 K,出现明显的差别.

3 结论

本文在三维球头烧蚀模型中引入多孔介质的分形渗透率张量模型,通过与实验对比验证了模型的可靠性.并分别研究了孔径分形维数和迂曲度分形维数对热解气在三维球头烧蚀过程中扩散的影响,得出如下结论:(1)孔径分形维数越小的材料热防护效果越好.烧蚀70 s时,分形维数分别为 1.7,1.8和1.9材料制成球头的内部温度分别为388，393 和398 K.(2)迂曲度分形维数越大的材料热防护效果越好.烧蚀70 s时,分形维数分别为 1.05,1.1和1.15材料制成球头的内部温度分别为418,393和378.通过研究不同烧蚀材料对结构热防护效果的影响,有助于更好地研究发展航空航天器的热防护材料和热防护结构.