改性C/C-ZrC-SiC复合材料热化学烧蚀数值分析 ①

梅宗书，范学领，石成英，晋小超

(1.火箭军工程大学 核工程学院，西安 710025；2.西安交通大学 航天航空学院 机械结构强度与振动国家重点实验室，西安 710049)

0 引言

碳/碳(C/C)复合材料因其突出的耐烧蚀、高比强度、高比模量、低热膨胀系数、耐热冲击等一系列优异性能，广泛应用于航天、航空领域，如再入飞行器鼻锥、固体火箭喷管和航天飞机机翼前缘[1]。但C/C材料在高温下易氧化，从而降低了其抗烧蚀性能。基体改性通过在碳基质中加入改性抑制剂或密封剂封闭活性点来降低氧化速率从而提高抗烧蚀性能。改性材料根据引入相的不同，可分为SiC陶瓷、UHTCs、难熔金属等，其中UHTCs颗粒掺杂物主要是难熔金属Zr、Hf和Ta的碳化物及硼化物，研究较多的有ZrC、ZrB2、ZrB2-SiC、HfB2-SiC、ZrB2-ZrC-SiC等，具有较强的应用前景，可在2200 ℃以上高温对烧蚀结构进行有效防护[2]。

在数值计算方面，碳基复合材料主要从表面热化学烧蚀和内部材料热响应两个方面，对其烧蚀表面气体产物浓度、质量烧蚀率、线烧蚀率等性能参数进行计算求解。YIN Tiantian等[3]根据C/C热化学烧蚀在低、高温下分别由反应动力学速率和扩散速率控制的特点，研究了C/C复合材料前缘结构在不均匀压力、温度分布下的表面退移率，同时探讨了基体改性对材料烧蚀性能的影响。杨德军等[4]、WANG Chen[5]根据表面能量平衡、质量守恒及化学平衡方程构建表面烧蚀模型，并通过虚拟失效、重新构建网格的方法，实现了热化学烧蚀、烧蚀表面退缩及温度场耦合作用下C/C复合材料烧蚀性能变化规律有限元分析。有限速率表面烧蚀模型[6-8]是C/C、碳/酚醛复合材料表面热化学烧蚀计算的又一重要而有效的途径，建立了一种基于表面元素质量守恒的有限速率烧蚀表面边界条件公式群，可预测再入空间飞行器的热罩表面的气动加热、化学物质浓度和碳质材料烧蚀速率。目前，改性C/C复合材料的研究主要是通过试验的方法[9-11]，数值计算大多针对未掺杂改性碳/碳复合材料而言[12-14]，改性C/C复合材料的数值分析则较少[11,15]，本文通过热化学平衡常数和表面元素质量守恒得到ZrC-SiC改性C/C复合材料烧蚀表面产物浓度，并构建CCZS表面氧化模型。通过烧蚀试验结果数据和热-力-化多场耦合数值分析计算结果比较，对构建的改性C/C复合材料氧化模型进行验证。

1 改性复合材料烧蚀行为数值分析

利用流体动力学和结构有限元数值分析方法，结合流-固耦合多场耦合方案，对CCZS复合材料的烧蚀行为和质量烧蚀率进行数值分析，将多场耦合数值分析方法的范围从C/C复合材料扩展到基体改性C/C复合材料。CCZS复合材料的烧蚀是一个多组分、多相态(固态、液态、气态)的复杂热化学反应过程，其最终产物的确定目前仍是基于化学平衡的研究方法。本章将基于热化学平衡原则和最小自由能原理，对烧蚀体系中的各组分分压和浓度进行理论推导和数值分析。

1.1 热化学表面氧化机制

CCZS复合材料在高温气流作用下表现出复杂的热化学烧蚀行为，涉及多种化学反应和气体扩散。根据最小自由能判断，CCZS复合材料在氧化烧蚀过程中可能发生的化学反应如式(1)～式(10)所示，热化学烧蚀包括碳纤维、碳基体、碳化锆和碳化硅的氧化反应、碳氮反应及碳的升华等。

3C(s)=C3(g)

(1)

SiO2(s)=SiO2(1)

(2)

C(s)+O2(g)=CO2(g)

(3)

2C(s)+O2(g)=2CO(g)

(4)

2C(s)+N2(g)=2CN(g)

(5)

4C(s)+N2(g)=2C2N(g)

(6)

ZrC(s)+2O2(g)=ZrO2(s)+CO2(g)

(7)

SiC(s)+2O2(g)=SiO2(s)+CO2(g)

(8)

本文试验设定了三种烧蚀温度，分别为1400、1600、1750 ℃，在这三种不同烧蚀温度下，温度越高，发生的热化学反应越多，本节以烧蚀温度为1750 ℃的工况进行组分浓度表达式的推导，在此温度下，式(8)中产生的固态SiO2会发生相变，呈液态，会对碳材料的氧化反应起到阻隔氧化性气体的作用。但因其含量较少，在进行氧化模型构建时，忽略SiO2熔化对氧化反应的影响，即SiO2熔化后对氧化性气体的阻隔作用忽略不计。CCZS复合材料在烧蚀过程中的表面热化学反应，除了式(1)～式(8)中所列出的材料固体相与气体的反应，还应包括气体自身的离解。

O2(g)=2O(g)

(9)

N2(g)=2N(g)

(10)

1.2 热化学反应模型

当反应达到平衡时，根据化学平衡计算原理，边界反应层中产物分压和反应物分压之比为一个以温度为变量的平衡常数，其表达式为

(11)

其中，pp,i、pr,i分别为式i(i=1～10)中产物分压和反应物分压；α,β分别为式i中产物和反应物的系数。式(1)～式(10)所对应的平衡常数系数a、b经过查表[16]后拟合可得表1。

表1 平衡常数系数

组元浓度Ci和组元分压pi之间的关系式为

(12)

(13)

根据元素质量守恒和相容性条件，可得元素浓度和组元浓度之间的关系为

(15)

(16)

(17)

(18)

式中f1为复合材料中ZrC质量分数；f2为复合材料中SiC质量分数；CO,e、CN,e分别为边界层外缘氧元素和氮元素质量分数；B为无因次质量损失率。

联立式(14)～式(16)，可得

(19)

从而得到无因次质量损失率为

(20)

1.3 产物浓度计算

根据式(1)～式(10)得到的反应平衡常数关系式和式(17)、式(18)、式(20)可得以下各组元浓度表达式：

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

上列方程组需迭代求解，先给定B和T，由式(21)～式(30)可求出全部壁面组分质量分数，将算出的浓度代入式(20)，可求出一个新的B值。若给定的B与求出的B之差满足精度要求，则给定的B为所求的解；若不满足，需要重新给出B，反复迭代，表面温度T由表面能量平衡方程确定。经过数值计算，得到CCZS复合材料表面烧蚀层的各组元浓度。

2 CCZS试样烧蚀试验

2.1 试样及烧蚀平台

采用聚合物浸渍裂解(Precursor Infiltration and Pyrolysis，PIP)结合化学气相沉积(Chemical Vapor Deposition，CVD)制备CCZS材料。为研究材料在三种试样表面温度下的抗氧化烧蚀性能，制备了3颗形状及材料相同的锥体烧蚀试样。编号1～3的试样分别对应表面温度为1400、1600、1750 ℃的三种工况。由于制备工艺的限制，三颗烧蚀试样从尺寸上有微小的差别，但这种外形尺寸上的微小差别对地面烧蚀流场和材料烧蚀行为的影响程度是有限的。因此，在做烧蚀试验及试验后的分析时，可认为三颗试样是完全相同的。

烧蚀体几何尺寸和基本物理属性见表2，试样的密度约为(2.795±0.014) g/cm3，掺杂浓度(15±0.29)%，高度(11.2±0.04) mm，底部直径(10±0.03) mm，锥角为42°±0.5°。由于三颗试样的物理属性差别是很小的，在后面的分析中，不考虑三颗试样本身的差异对烧蚀性能的影响，而主要讨论不同工况下试样的烧蚀行为及性能的差异。在下文的数值模拟建模中，为了简化建模过程，烧蚀体的物理参数设为密度2.795 g/cm3，高度11.2 mm，底部直径10 mm，锥角42°。

采用氧-丙烷热冲击试验平台进行地面烧蚀试验，通过调节O2和C3H8的体积流量进行烧蚀温度调节。图1为烧蚀试验设备及不同工况下的试样表面温度和气体流量随时间的变化曲线。图1中，温度在设定的目标温度上下有微小的波动，温度的稳定区域可满足烧蚀时长的需要。因此，工况是满足试验要求的。

(a)Oxygen-propane thermal shock test platform (b)1400 ℃

采用氧-丙烷火焰对1～3号试样分别在1400、1600、1750 ℃的高温下烧蚀90 s，火焰喷口离试样驻点200 mm，火焰方向与试样轴线平行。试样烧蚀过程中，可观察到在1750 ℃时试样表面有亮斑，推断为SiO2，推断基于以下两点：一是排除为ZrO2，因其熔点为2677 ℃，且具有低的蒸汽压，SiO2的熔点为1670 ℃，且具有较高的蒸汽压，在1750 ℃的气流下，只可生成熔融的SiO2；二是亮斑的位置大多在锥体中间和尾端部位，应是由于气流的冲刷作用导致的，符合液态SiO2流动性好的特点。

2.2 烧蚀形貌分析

三种工况下，三颗试样烧蚀后的宏观形貌如图2所示。根据图2(a)、(b)可知，在较低温度(≤1600 ℃)下，试样的宏观形貌变化并不明显，但试样表面形成蓬松易脱落的物质。

(a)1400 ℃ (b)1600 ℃ (c)1750 ℃

通过扫描电子显微镜(SEM)和能谱仪(EDS)，对试样表面的微观形貌进行观察。研究发现，表面蓬松物质主要为ZrO2和SiO2呈颗粒形式的堆积物。根据图2(c)可知，在高温下试样表面形成熔融SiO2，并在局部位置聚集成微小颗粒附着在试样表面，烧蚀结束冷却后在试样表面呈现出白色颗粒。观察图2(c)试样的微观形貌，可观察到明显的白色线状物质，即试样表面形成的熔融的SiO2，由于SiC含量较少，所以没有形成成片的SiO2保护膜，但仍对保护内部碳质材料和氧化形成的ZrO2起到了重要的作用。

在1400、1600、1750 ℃烧蚀90 s工况下，CCZS试样质量损失速率较低，分别为0.838、0.875、0.881 mg/s，表明CCZS复合材料具有良好的抗烧蚀性能。根据化学反应速率可知，热化学烧蚀一般随温度的增加呈指数型增加，CCZS试样的质量损失速率虽然随着温度的增加提高，但相比于由1400 ℃到1600 ℃，由1600 ℃到1750 ℃时烧蚀速率增加并不明显。根据上述微观形貌的分析结果可推测，试样表面的熔融SiO2能够起到隔绝氧气的作用，对内部材料具有较好的保护作用，并使疏松的ZrO2不易被气流冲走，且随着温度的升高，生成的SiO2量增多，从而降低了随着温度升高而引起的质量损失。

3 多物理场耦合分析

3.1 烧蚀流场数值分析

利用Fluent中的非预混燃烧模型表征氧气与丙烷的燃烧反应，其中丙烷作为燃料流，氧气作为氧化剂流。非预混燃烧计算中不使用有限速率化学反应模型，而是用统一的混合物浓度作为未知变量进行求解。因此，无需计算代表组元生成或消失的源项。

对燃气发生器喷嘴和氧-丙烷焰燃烧状态进行仿真，建立二维氧-丙烷焰热流场物理模型，根据物理模型对称，建模时只考虑其1/2模型，如图3所示。模型尺寸为50 mm×520 mm，针对发生器的3层小孔，构建3个1 mm混合气体出口，每个出个间隔2 mm。沿着发生器喷口方向提供速度为1 mm/s的流动空气进行替代。构建烧蚀试样壁面，在壁面、燃气出口等区域进行网格加密处理。

图3 氧-丙烷热流场计算域网格划分

对模型边界条件进行设置，三个燃气入口均采用质量入口(Mass flux inlet)，入口参数需要给定气体在出口处的质量流率和初始静压。对于燃烧过程在开放环境下的自由射流场，出口边界选择压力出口边界条件(Pressure outlet)。对于喷嘴部分，除了出口通气孔外，其他部位均为壁面边界条件。上边界同样选为出口压力边界，下边界选为轴对称(Axisymmetric)。

对三种工况下的热流场进行稳态数值计算，流场温度分布云图如图4(a)所示。由图4(a)可知，随着燃气流量的增加，流场最高温度并没有明显的升高，推测热流场温度应与燃气中氧气和丙烷的混合比例有关。三种工况1400、1600、1750 ℃对应的试样驻点处的温度分别为1471.78、1677.57、1931.13 ℃，流场对称轴上的温度变化曲线如图4(b)所示。烧蚀试样驻点处的温度基本上与测量得到的主控温度一致，说明流场模拟的准确性，可作为多场耦合数值分析的流场部分使用。

(a)Cloud map of flow field temperature distribution at T=1750 ℃ (b)Variation curves of temperature

3.2 烧蚀试样建模与设置

利用Abaqus对烧蚀试样进行有限元建模分析，结构的尺寸和坐标与Fluent中的一致。试样各项尺寸参数根据2.1节中的设置，分别为密度2.795 g/cm3，高度11.2 mm，底部直径10 mm，锥角42°，结构建模如图5所示。网格采用结构化的四边形网格和轴对称传热-变形耦合单元CAX4T，单元数为4875，节点数为5061，在烧蚀最为剧烈的端头部分，进行了网格加密处理。

图5 结构有限元模型网格划分

(31)

式中 下标w为固体表面；下标e为来流边界。

表3 材料热力学性能

(32)

式中fr为机械剥蚀因子，一般取1或1.5，本文取为1.5。

总的材料质量损失率和线烧蚀率为

(33)

(34)

将烧蚀速率嵌入到表面烧蚀退移子程序Umeshmotion中，对CCZS复合材料烧蚀试样热-力-化多场耦合模型进行计算。

3.3 数值计算结果分析

图6为三种烧蚀温度下烧蚀试样内部温度分布，根据表面温度可知，模型可较准确地模拟烧蚀试验的温度分布。其中，最大温度处驻点的温度分别为1384、1566、1734 ℃。得到的表面温度均比试验测得的主控温度1400，1600，1750 ℃要低，推理可能是如下两种因素导致：

(1)主控温度的测量采用的是红外测温仪，火焰可能会对红外测温仪造成影响，从而产生误差；

(2)流场能量项没有考虑表面反应热，从而对流场热量计算时偏低。

(a)T=1400 ℃ (b)T=1600 ℃

由于热膨胀的作用，且热膨胀导致的计算节点正位移要比氧化烧蚀导致的负位移要大，即膨胀速率大于烧蚀退移速率，计算得到的试样尺寸随着烧蚀时间逐渐变大，尤其是轴向的长度变化较大。随着烧蚀温度的升高，热膨胀占主导作用的节点正位移逐渐提高，但单位温度增长下提高的温度幅度明显降低。原因是随着温度的升高，热膨胀率有减低趋势；另外，温度越高，氧化导致的负位移越高，抵消了部分热膨胀位移。

图7为三种工况下，试样驻点温度随时间的变化，线型数据表示通过数值计算得到的数据，符号型数据代表试验测得的试验数据。

图7 驻点温度变化曲线

由图7可知，在烧蚀的前10 s内，驻点温度上升最快，之后随着烧蚀时间增加，温度上升速率逐渐降低，在烧蚀约40 s后，温度已经基本达到了主控温度，之后略微有所上升。试验数据和计算数据在烧蚀时长为10 s后吻合度较高，10 s之前由于温度处于急速上升阶段，存在一定的误差。整体而言，试验数据与计算数据变化趋势一致，表明所构建模型的有效性。

图8为烧蚀90 s时长中，试样的轴向表面位置出的质量损失率。由图8可知，在前30 s内，质量损失率增长较快，之后增长趋于平缓，烧蚀60 s后，质量损失率逐渐趋于稳定。稳定后的三种工况下，表面温度由低到高，试样驻点的质量损失率分别为2.41×10-3、2.83×10-3、3.23×10-3kg/(m2·s)。比较驻点处和尾部的质量损失率可知，两者之间的质量损失率大小相差较小，氧化烧蚀在试样空间位置的差异不是很明显，原因应是试样尺寸较小，流场参数在试样表面的几乎相同。

(a)T=1400 ℃ (b)T=1600 ℃

将得到的单位面积的质量损失率通过乘以试样受氧化面积，得到单位时间试样整体的质量损失率，其单位变为mg/s。图9曲线数据为三种工况下，试样整体质量烧蚀率随时间的变化。

图9中的矩形数据为烧蚀试验得到的不同工况下烧蚀时长90 s下的平均质量损失率。可知，随着温度的升高，烧蚀质量损失率增加趋势变缓，这与上文中通过分析掺杂超高温陶瓷颗粒烧蚀机理时得到的结论是一致的。同时，与试验测量得到的平均烧蚀率相比，也说明计算数据与试验数据得到的质量损失速率较为一致，进一步证明了模型的有效性和可靠性。

图9 烧蚀质量损失率

4 结论

利用热化学平衡原则进行改性复合材料氧化烧蚀行为的表征，推导并验证了基体改性C/C复合材料的氧化模型。仿真与试验验证了模型的有效性，为发展改性材料的氧化模型提供了理论参考。文中的理论方法是一种简化近似计算方法，没有对液相产物在阻隔氧气扩散和降低覆盖点位处材料剥蚀率等方面的影响，因此接下来需要对液相产物的物理尺寸及附着在材料表面的分布位置和稀疏度等分布属性进行研究，从而获得液相产物更为精确的防热机制和表征方法。