C/C复合材料微观尺度烧蚀过程质量损失速率模拟*

宋永善，齐乐华，张守阳，张佳平，李逸仙，李贺军

(1.西北工业大学 机电学院，西安 710072；2.西北工业大学 C/C复合材料研究中心，西安 710072)

0 引言

C/C复合材料具有一系列优异的性能[1-3]，如高强度、高模量，以及其强度随温度升高不降反升，此外还具有良好的抗热冲击、抗烧蚀、热膨胀系数低及良好的导热性能，是迄今为止惰性气氛下最为理想的高温结构材料之一，被广泛用于固体火箭发动机喉衬、飞行器头锥等航空航天结构部件。然而，在高温、高速冲击等极端苛刻的服役环境下，C/C复合材料不可避免地会发生烧蚀现象，进而导致材料性能急剧下降。

目前，关于C/C复合材料的烧蚀问题已经开展了很多实验研究，但受烧蚀过程中的高温限制，实验研究所能测量的参数只有线烧蚀率和质量烧蚀率，烧蚀机理的研究也大多基于对烧蚀后材料表面的微观形貌进行分析，得到的信息十分有限。数值模拟因此成为进一步深入研究C/C复合材料烧蚀机理不可或缺的手段。Piyush Thakre等[4]和Daniele Bianchi等[5]基于不同假设分别计算了C/C复合材料发动机喷管的热化学烧蚀过程。王臣等[6]通过对碳基复合材料热化学烧蚀机理分析，根据质量守恒、能量守恒及系统化学反应的热化学平衡原理，建立了相应的热化学和热力学烧蚀模型，模拟了材料内部热应力场分布。杨德军等[7]采用虚拟失效、重新构建网格部件的方法实现烧蚀表面的退缩，建立了烧蚀表面退缩下瞬态温度场的有限元模型，分析了热化学烧蚀、烧蚀表面退缩及温度场耦合作用下 C/C复合材料烧蚀性能的变化规律。解惠贞等[8]研究表明，纤维分布方向对C/C复合材料的烧蚀性能有重要影响，垂直燃气流方向的纤维与平行燃气流方向的纤维的线烧蚀率和质量烧蚀率区别非常大。汪海斌等[9]针对轴编C/C复合材料的结构形式和烧蚀机理，建立了喷管喉衬烧蚀的多尺度分析方法。通过宏观-微观的渐进分析，获得了喷管喉衬的烧蚀率和烧蚀形貌。王臣[10]建立了C/C复合材料的表面氧化模型来研究材料的烧蚀过程和热环境参数，并与实验结果吻合良好。上述研究主要针对C/C复合材料的宏观烧蚀机理分析，未考虑C/C复合材料中炭纤维和碳基体在密度和性能上的差异，因而很难揭示微尺度下C/C复合材料烧蚀过程中的形貌变化。实际上，C/C复合材料微尺度下的形貌变化会导致其宏观尺度的表面粗糙化，进而加剧材料表面的传热、传质，影响材料表面的流场。Lachaud J等[11]通过理论计算结合实验测定了C/C复合材料及其组分(纤维、基体)的固有氧化速率，建立了稳态下C/C复合材料氧化模型，推导出的函数能直接计算氧化达到稳态时纤维的形貌，但未能给出烧蚀过程中的微观形貌演变规律。

本文建立微观尺度烧蚀过程的形貌演变瞬态模型，研究C/C复合材料微观烧蚀过程中纤维和基体的相对质量损失规律，以期为研究C/C复合材料微观烧蚀机理提供借鉴。

1 烧蚀过程建模

1.1 基本假设

在承受高温高速冲刷的苛刻服役环境中，C/C复合材料会因氧化、升华、等离子腐蚀和机械剥离而被逐渐破坏，导致不能忽略的表面退移，该现象被称为烧蚀。同时，烧蚀还会产生其他作用[5]：(1)增大反应表面，进一步加剧发生在材料表面的传热和传质；(2)表面粗糙化导致C/C复合材料的有效烧蚀速率不能简单地用其组分的烧蚀速率算术平均代替；(3)加剧机械剥离；(4)导致层流向湍流转变。

但材料烧蚀率较低，对气体流动的连续性影响较小，故本文对烧蚀模拟过程做了简化，做出如下假设：(1)气源位于离初始材料表面足够远的一个距离；气体流动状态为层流；(2)氧化反应为一阶反应。则材料表面的氧化速率可定义为

v=-kCVmol

(1)

式中v为材料的氧化速率，m/s；k为气体反应速率常数，m/s；C为局部反应气体浓度，mol/m3；Vmol为材料的摩尔体积，m3/mol。

1.2 水平集法

本文采用水平集法来追踪烧蚀过程中材料的表面移动。材料表面由水平集函数定义：

(2)

其中，d为到材料表面的最短距离，φ=0代表材料表面。在定义了φ的初值及合适的速度场v后，可以通过计算水平集方程(3)来的追踪材料的表面移动。

(3)

但烧蚀速度v只在材料表面有物理意义，针对这种情况，必须通过外推构造合适的速度场。本文选用方程(4)所示方法[12]进行外推：

(4)

式中φ为速度v的任意一个分量。

方程(4)中符号函数S(φ)和单位法向量n分别定义为

(5)

(6)

在计算方程(3)一个步长后，φ不再为符号距离函数，必须重新初始化。通过求解方程(7)直至稳态可以完成重新初始化。

(7)

方程(7)中sign(φ0)被定义为

(8)

式中τ为虚时间；φ0为求解水平集方程(3)一个步长以后得到的水平集函数。

1.3 扩散方程

气体区域和气-固界面分别采用扩散方程(9)和(10)描述氧化过程中反应气体的质量传输：

(9)

-DCn=-kf/mC

(10)

式中D为反应气体的扩散系数；kf/m为由材料组分决定的非均相反应速率常数。

根据边界条件(7),可得

v≅DCnvmol

(11)

将式(11)等号右侧的r向z向分量分别代入方程(4)外推，可得到外推后的速度场vext。

氧化过程中，材料表面不断移动，这意味着每计算一个步长，都要重新设定边界条件。为了避免这种情况，本文选择文献[13]中的近似方法。对边界条件(10)作如下处理。

氧化过程中，反应气体在材料表面处的消耗SC可表示为

(12)

函数δ(φ)表示为

(13)

其中，ε为正比于网格大小的参数。函数δ(φ)的特点为在区间[0,2ε]内的积分1,且在φ=ε处取最大值。

经过处理后的边界条件(12)可整合为方程中的源项：

(14)

扩散系数被定义为φ的函数：

(15)

气体中扩散系数为D，为了避免材料表面附件的数值震荡，材料一侧气体扩散系数被定义为一个非常小的数εD。

1.4 计算方法

在COMSOL中建立如图1(a)所示的2D轴对称模型。如图1(b)所示，模型中纤维半径为r，基体厚度为t，材料初始长度为h1，气源距材料表面的初始距离为h2。模型通过映射方式划分正方形网格，网格尺寸为10-5mm，气源处边界条件为给定的初始浓度。

求解C/C复合材料烧蚀形貌的具体计算流程图如图2所示。模拟中通过Comsol中的系数型PDE模块写入扩散方程(14)，定义相应的初始值和边界条件。其中，C和φ均通过COMSOL with MATLAB以插值函数形式分别导出和导入COMSOL。

(a)模型基本框架 (b)局部网格

图2 C/C复合材料烧蚀过程计算流程Fig.2 Flow digram of computation for C/C composites during ablation

不同时刻的相对质量Δm的计算通过COMSOL后处理中的体积积分功能实现，具体过程如下：

(1)在COMSOL后处理模式中计算各时间步长下的纤维体积Vf。计算公式为

Vf=(1-H(φ))dVindomainr≤rf

(16)

Vm=(1-H(φ))dVindomainr>rf

(17)

式中r为2D轴对称坐标下的横坐标变量；rf为纤维半径值；H(φ)为定义的Heaviside函数。

通过引入Heaviside函数可以很好地平滑界面处的不连续，得到更准确的结果。

(2)记初始时刻的纤维体积为Vf0，其余时刻的纤维体积为Vfi。则不同时刻纤维的相对质量Δmi=Vfi/Vf0，相应的基体的相对质量Δmi=Vmi/Vm0。

2 烧蚀过程相对质量损失速率分析

2.1 稳态结果验证

本文建立了半径为3.5 μm的单根纤维及包裹在纤维外围的厚度为0.5 μm的界面层组成的C/C复合材料模型。烧蚀的温度大约为3000 K，设定此时的气体扩散系数为7×10-4m2/s[14]；此时界面层的氧化速率常数km=182 m/s[15]。图3(a)为D=7×10-4m2/s，kf=20 m/s，km=180 m/s时计算得到的稳态形貌，纤维形状呈不规则的圆锥形，纤维母线为曲线；图3(b)是单向C/C复合材料氧乙炔焰后的表面形貌，烧蚀温度2500 ℃，O2和C2H2气体流量分别为0.24 L/s和0.10～0.20 L/s，烧蚀时间为60 s，火焰距试样表面距离为10 mm,焰流密度为2.38 MW/m2，与图3(a)中的仿真形貌吻合良好。

2.2 烧蚀过程相对质量损失速率

在氧化过程中，材料形状发生了剧烈变化，增大了纤维的露出表面，增加了基体距离气源的相对距离。为了准确分析这种变化对纤维和基体质量损失速率的影响，本文对纤维和基体微观形貌达到稳态之前的质量损失过程进行了分析。

Jean Lchaud提出了C/C复合材料微尺度的稳态烧蚀模型[16]，模型定义了Sherwood数Sh和反应速率比A，本文以这两个参数作为变量，计算了不同反应/扩散机制下C/C复合材料的微观形貌演变过程。由于材料的形貌演变过程十分缓慢，只保存形貌演变过程中的特定时间的数据。特定时间的选择基于Vignoles G L[17]提出的参考时间公式:

τ0=rf/(C0υskf)

(18)

图4为A=5，Sh分别取0.01、0.1、1时的纤维和基体相对质量损失曲线；图5为Sh=0.1，A分别取2、4、8时的纤维和基体相对质量损失曲线。其中，Δm为纤维当前的相对质量。图6为A=5时反应气体的浓度变化云图。

(a)计算结果 (b)实验结果

(a)A=5，Sh=0.01时,纤维 (b)A=5，Sh=0.1时,纤维 (c)A=5，Sh=1时,纤维

(d)A=5，Sh=0.01时,基体 (e)A=5，Sh=0.1时,基体 (f)A=5，Sh=1时,基体

由图4(a)可见，当Sh=0.01，A=5时，在C/C复合材料微观形貌演变过程中，纤维质量损失速率迅速增大，在t=τ0时达到最大值。

由图6(a)可以发现，在整个形貌演变过程中，材料表面反应气体浓度值基本一致，大致保持0.85C0的高浓度，这是因为反应气体的扩散系数较高。但观察形貌演变过程中5条浓度等值线(0.85C0～0.9C0)的变化，可以发现，最初浓度等值线完全紧贴材料表面，与纤维的轮廓相仿，露出纤维表面的浓度一致，但随着材料形貌演变的进行，浓度等值线不再保持之前的形状，露出纤维底部的浓度低于纤维顶部的浓度，基体表面的反应气体浓度降低。因此，纤维质量损失速率并没有随着纤维露出表面的增大而持续增大，在之后受扩散的制约，纤维露出表面的增大不能加剧纤维的质量损失。

图4(b)、(e)为当A=5，Sh=0.1时计算得到的烧蚀过程，τ02为参考时间。可以发现与图4(a)类似的现象，在形貌演变初期(t=0～τ02),基体被大量消耗，使纤维露出表面增大，纤维质量损失速率在烧蚀初期就达到了最大值，但质量损失速率的增幅已不太明显。随着烧蚀的进行，由于反应气体扩散系数较低，纤维底部与纤维顶部出现了明显的浓度差(图6(b))，但露出纤维表面的增大还能够略微增大纤维的质量损失速率。但对基体来说，其质量损失速率迅速减小，保持稳定(图4(f))。

(a)Sh=0.1,A=2时,纤维 (b)Sh=0.1,A=4时,纤维 (c)Sh=0.1,A=8时,纤维

(d)Sh=0.1,A=2时,基体 (e)Sh=0.1,A=4时,基体 (f)Sh=0.1,A=8时,基体

(a)Sh=0.01 (b)Sh=0.1 (c)Sh=1

由图4(c)可见，当Sh=1，A=5时，纤维质量损失速率几乎保持恒定，说明烧蚀过程为纯扩散机制控制。由图4(f)可见，当Sh=1，A=5时，基体的质量损失速率迅速减小，并很快达到稳态。

由图6可知，当反应速率比一定时，随着Sh的增大，扩散反应堆氧化反应的抑制逐渐增大，当Sh=0.01时，烧蚀主要受氧化反应控制；当Sh=1时，烧蚀过程则主要受扩散反应控制。

由公式(18)可知，τ06=2τ05=4τ04，则随着反应速率比A的增大，基体质量损失基体相对于纤维的相对质量损失的变化进一步增大，且相对质量损失速率达到稳态所需的时间增加，如图5所示。

3 结论

(1)基于水平集算法，耦合扩散方程建立了C/C复合材料在纤维尺度下的2D轴对称烧蚀模型，计算了C/C复合材料的微观形貌演变过程。

(2)烧蚀至稳态的模拟结果与C/C复合材料烧蚀后的微观形貌吻合较好，验证了本文模型的可靠性。

(3)分析了Sh数和反应速率比A对纤维和基体的相对质量损失速率的影响规律。随着Sh的降低时纤维的相对质量损失速率会迅速增大并在烧蚀形貌达到稳态前稳定，基体的相对质量损失速率会迅速降低并在烧蚀形貌达到稳态前稳定；纤维和基体的相对质量损失速率达到稳态所需时间随反应速率比A的增大而增大。