ICVI工艺致密化C/C复合材料的数值模拟①

王建武，李克智，张守阳，李 伟

(西北工业大学碳/碳复合材料工程技术研究中心，凝固技术国家重点实验室，西安 710072)

0 引言

C/C复合材料不仅具有密度小、比模量高、比强度大、摩擦性能好等优异的力学性能，而且还保持了炭材料所固有的耐高温、耐热冲击、抗烧蚀、热膨胀系数低等良好的热物理性能，被广泛应用于航空航天领域［1－2］。等温化学气相渗透(ICVI)工艺是制备 C/C复合材料最主要的工艺，但在沉积过程中，为了避免制件表面过早结壳，通常在较低的温度和压力下进行沉积，导致整个制备周期过长，生产成本很高。因此，寻求最优的沉积工艺，实现C/C复合材料的快速沉积是实际生产中亟待解决的问题之一。

影响ICVI工艺的参数众多，主要包括沉积温度、压强、气体流量，若对各个工艺参数的影响规律设计实验进行探索，需要消耗大量的时间和物力。通过计算机模拟技术对ICVI工艺过程进行模拟，可研究沉积环境中的流场和温度场分布，显示预制体致密化的过程，对C/C复合材料ICVI工艺的研究具有十分重要的意义［2－3］。目前，ICVI工艺的计算机模拟受到越来越多的关注，包括预制体结构建模和致密化过程仿真、计算域内物理场的分布计算及前驱体反应的化学反应动力学研究。Huttinger等［4］认为，烯烃、芳香烃和大分子烃是主要成炭的物质，以毛细管作为孔隙简化模型，建立了以甲烷为前驱体的ICVI工艺模型，并进行了实验验证。李爱军［5－6］建立了双孔隙模型和“平行-串式”反应模型，分析了预制体致密化过程中孔隙拓扑结构的变化规律，以及均相反应和非均相反应的竞争关系与热解炭织构的关系。孙国岭、赵建国等［7］直接运用NS方程描述了前驱体在自由空间的流动，但并未给出预制体内部流场的分布规律。焦妍琼、边国军［8－9］等分别建立了2D轴对称瞬态模型，利用热传导、对流和扩散、传质、沉积反应的多场耦合计算，模拟了热梯度CVI中流场、温度场的分布及密度的演变规律。

本文针对ICVI工艺，根据实验用的沉积装置，建立简化的几何模型，依据质量、动量和能量守恒定律建立各个区域的偏微分控制方程，计算沉积炉内流场和温度场的分布，分析前驱体进入沉积炉后的流动速率以及升温过程;分析整个沉积过程预制体密度演变规律，并对比整体平均密度随时间变化的计算值和实验值，验证模型的可靠性。

1 几何模型建立和模型假设

ICVI的工艺实验装置如图1(a)所示。前驱体从底部通入，流经窄缝区时通过扩散方式进入预制体，在预制体内发生热解反应生成热解炭，并沉积在纤维表面，反应副产物随通道中未反应的气体从上方的出口排出。对沉积装置进行简化，建立了如图1(b)所示的2D轴对称几何模型。其中，Ω1为预制体区，Ω2为自由流动区，Ω1，1、Ω1，2分别表示预制体内、外侧边界，Ωin、Ωout分别为气体的入口和出口，Ωw为炉壁。模拟采用的沉积工艺参数为沉积温度为1 393 K，压强为10 kPa，天然气在入口处的流速为4 m/s;预制体初始孔隙率为 0．675，初始密度为 0．575 kg/m3。

另外，对模型作进一步假设:

(1)不考虑碳氢气体热解反应对系统热量的影响;

(2)反应气体为理想气体，且不可压缩，预制体为各向同性体;

(3)热解反应均为一级反应，热解的中间产物仅考虑 C2H4、C2H2、C6H6、H24 种主要的组分，反应采用简化的多步反应模型，如图2所示。其中，k1～k4、K1～K4分别为气相反应和气固反应的反应速率常数，其值见表2，C∞表示热解炭;

(4)预制体区的比表面积远大于自由流动区的比表面积，故仅考虑发生在预制体区内的沉积反应。

图1 ICVI工艺原理图及简化后的2D轴对称几何模型Fig．1 diagram of ICVI furnace and simplified 2D axis-symmetric geometry model

图2 简化的甲烷热解的多步化学反应动力学模型Fig．2 Simplified multi-step chemical reaction kinetic model of methane pyrolysis

2 数学模型的建立

2．1 传热方程

ICVI过程中炉壁温度设为沉积温度Tw，热量以前驱体和预制体为媒介通过热传导和对流的方式从炉壁向内部传递，根据能量守恒定律，分别建立Ω1区和Ω2区的热传导方程:

其中

式中 ρ、ρp分别为甲烷和预制体的密度，kg/m3;cp和cp，p分别为甲烷和预制体纤维的常压比定压热容，J/(kg·K);k、kp分别为甲烷和预制体纤维的热导率，W/(m·K);θp为预制体中纤维的体积分数;u为速度矢量。

炉壁的温度为沉积温度，出口处仅考虑热对流，得到以下边界条件和初始条件:

在流场和温度场的耦合计算中，并未考虑甲烷的裂解。所以，仅需要计算甲烷的比定压热容和热导率。按照多项式拟合，得到［11］:

预制体纤维的比定压热容和热导率:

2．2 流动方程

气体在Ω1区和Ω2区内的流动方式不同，根据动量守恒原理，分别建立如下的弱可压缩气体流场方程和连续方程［1，9］:

式中 ρ为混合气体的密度，kg/m3;ε为预制体的孔隙率;η为混合气体的粘度系数，Pa·s;κ为预制体的渗透率，设为 0．25［9］;p 为炉内压强，Pa;I为单位矢量。

入口边界和出口边界分别设为速度和压强，边界条件和初始条件如下:

按照混合法则，混合气体密度计算方法如下:

式中 ρi为第 i组分气体的密度，kg/m3;Ci为第 i组分气体的浓度，mol/m3;n为气体组分数。

其中，纯组分气体的密度根据理想气体状态方程计算:

对于混合气体的粘度，采用Sutherland模型进行估算［11］:

式中 ηi和yi是纯组分i的粘度和分子分率［12］。

式(8)中的关联系数φij可由简单的wilke近似式计算:

2．3 传质方程

前驱体的温度迅速升高到裂解温度后，会发生分解，生成其他组分的气体，见图2。混合气体通过扩散方式向预制体内传输的同时，发生热解反应生成热解炭，沉积在纤维的表面。根据质量守恒原理，对Ω1区和Ω2区分别建立浓度变化的偏微分方程:

式中 Ci为第i种反应物的浓度，mol/m3;Dieff、Diab分别为第i种组分在Ω1区和Ω2区的扩散系数，m2/s;R1，i、R2，i分别为第 i种反应物在 Ω1区和 Ω2区的反应速率，mol/(m3·s)，见表 1。

其边界条件和初始条件为

综合考虑Fick扩散和Knudsen扩散，得到有效扩散系数表达式［11］:

式中 ε为预制体孔隙率;τ表示预制体内孔隙结构的曲折因子，取值一般为3～7。

Dik和Diab的值见表2。

表1 致密化过程中各组分在Ω1和Ω2区的反应速率Table 1 Reaction rate of species in Ω1and Ω2during the densification process

表2 各组分的Fick扩散系数和Knudsen扩散系数的值及各反应的反应速率常数ki和KiTable 2 Value of and and reaction rate constant kiand Ki

表2 各组分的Fick扩散系数和Knudsen扩散系数的值及各反应的反应速率常数ki和KiTable 2 Value of and and reaction rate constant kiand Ki

Di ab/k/组分Di(m2/s)(m2/s)ki Ki CH4 5．642 ×10－3 9．33 ×10－6 1．0454 2．44 ×10－6 C2H4 2．357 ×10－3 7．05 ×10－6 5．6011 5．43 ×10－6 C2H2 2．466 ×10－3 7．32 ×10－6 2．4585 2．50 ×10－5 C6H6 1．497 ×10－3 4．23 ×10－6 — 2．00 ×10－4 H2 8．679 ×10－3 2．64 ×10－5— —

2．4 孔隙率的变化方程

前驱体及其裂解产物在预制体区发生热解反应，生成热解炭并填充在孔隙中，使孔隙率不断降低，预制体孔隙结构及比表面计算均采用“双孔隙模型”，其变化过程可用如下方程表示:

式中 Mc、ρc分别为热解炭的摩尔质量和密度，kg/mol，kg/m3;R为生成热解炭的反应速率，也即热解炭的沉积速率，mol/(m3·s);它与预制体的比表面积Sv有关。

设定所有的边界条件为绝缘边界，初始条件为

3 计算求解

本文基于以上建立的气体流动、热量传递和孔隙率变化的偏微分方程，考虑到各个物理场之间的耦合关系，通过多物理场耦合软件(COMSOL Muitiphysics)，采用有限元方法进行迭代耦合计算，具体的计算过程如图3所示。

图3 ICVI致密化过程模拟流程图Fig．3 Flow diagram of densification simulation for ICVI process

4 结果与讨论

4．1 流场分析

图4(a)显示了t=5 min时预制体区和自由流动区内速度场的分布。可看出，前驱体在Ω2中间位置的速率远大于Ω1区内的速率，高出约3个数量级，且整个Ω1区域的速率接近于零。这是由于在预制体内有大量的微孔，整个区域的比表面积达到了2×105m－1，气态前驱体受到相当大得粘滞阻力。因此，前驱体的流动受到很大的限制，流动速率接近于零。

为了进一步了解在不同高度z处速度沿半径r方向的变化情况，在 z= －0．0 1m，z=0和 z=0．0 1m 处各取一个横截面，其速度-位移(v-r)曲线如图4(b)所示。该图显示在不同高度z上速度沿径向r的分布基本一致，即Ω2区为未完全发展的层流，窄缝中心部位流动速度较快，往炉壁一侧速度逐渐降为零，在靠近预制体一侧速度降低到0．5 m/s左右;Ω1区内速度接近于零，但并不为零。不同z值的3条曲线也说明，炉内流场的分布随高度的增加其变化并不显著，流场分布较为稳定。从曲线中还可看出，在Ω1区与Ω2区的交界面Ω1，2处流速的变化是连续的，这与实际流动是相符的。

图4 t=5 min时沉积炉内速度场的分布及不同高度z处速率沿r向的变化曲线Fig．4 Distribution of velocity in the furnace and velocity curves as a function of r at different z after 5 min

4．2 温度场分析

图5给出了t=5 min时整个区域内温度场的分布，以及不同高度z处温度沿径向r的变化曲线。从图5(a)可看出，入口附近前驱体的温度还较低，但前驱体与壁面经过短时间的热量交换后，其温度迅速升高，已接近设定的沉积温度，而且整个沉积炉内的温度分布基本一致，符合等温沉积的条件。图5(b)显示从下到上气体温度逐渐升高，这是由升温时间不同造成的，但上下温度差在2 K以内。z=－0．01 m处的曲线在靠近界面位置(0．03 m＜r＜0．04 m)温度有下降趋势。这是由于气体受热时间短，未被充分加热便传输到预制体内，且在同高度的Ω2区相对较冷的气体在z方向的快速流动还带走了部分热量;z=0和z=0．01 m位置的温度都随着r增加而升高，因为越靠近炉壁热量传输的距离越短，温度也越高。此外，通过炉壁温度和预制体中心温度的比较可看出，不同高度z位置内外的温度差都在3 K左右，温度相差都很小，进一步验证了前驱体的沉积过程是在等温条件下进行的。

图5 t=5 min时沉积炉内温度场的分布及不同高度z处温度沿r向的变化曲线Fig．5 Distribution of temperature in the furnace and temperature curves as a function of r at different z value after 5 min

4．3 致密化过程分析

图6给出了沉积温度为1 393 K，压强为10 kPa，天然气在入口处的流速为4 m/s时，预制体内密度分布的演变规律。致密化前期0～50 h，预制体中心位置(0～0．02 m之间)，热解炭的沉积很快，密度最大的区域出现在0．02～0．03 m之间，而在靠近窄缝的区域内(0．03～0．04 m 之间)，热解炭的沉积相对较少，如图6(a)。致密化中期50～100 h，预制体内密度的分布规律与前期相同，但密度最大的区域向外侧移动，靠近窄缝的密度较小的区域也在逐渐缩小，见图6(b)、(c)。致密化后期100～150 h，预制体内整体的密度分布较为均匀，只有在靠近窄缝附近的一个狭窄区域内密度偏低，见图6(d)。

图6 不同沉积时间预制体内的密度分布Fig．6 Density distribution of the preform at different deposition time

图6中还显示，在不同沉积阶段，靠近窄缝位置的预制体密度总是很低。这是因为靠近窄缝的区域，气体在z向的流动比较强，将未来得及反应的气体带走，反应的气体在该区域内的滞留时间很短。所以，致密化的效果较其他位置较差。

从密度最大值的变化也可反映出整个预制体致密化的进程。致密化前期0～50 h，最大密度增加为0．945 9 g/cm3，增加近一倍;致密化中后期50～100 h和100～150 h，最大密度增加分别为0．408 1 g/cm3和0．082 4 g/cm3，增加幅度明显降低。由此可很明显地看出，预制体密度的增加主要发生在致密化前期，占整个密度增加的66%，越到沉积后期，密度增加越困难，随着沉积的不断进行，致密化效率也在不断降低。这是由于沉积前期，预制体内孔隙率和比表面积都很大，前驱体及其裂解气体不仅可在预制体内很好地进行传递，及时补充因沉积而消耗的气体，而且热解炭的沉积反应也比较剧烈;沉积进行到后期，孔隙已经被大量的热解炭填充，气体传输的通道被堵塞，沉积的气体得不到及时补充，沉积的速率明显降低。

图7给出了致密化150 h过程中预制体整体平均密度随时间的变化曲线及不同致密化时间预制体密度的实验值。

图7 致密化150 h预制体整体平均密度随沉积时间变化的模拟曲线和实验结果比较Fig．7 Comparison of average density of preform between simulated and experimental results after 150h densification

沉积140 h，预制体平均密度的计算值为1．871 g/cm3，实验值为 1．845 g/cm3，误差为 0．026 3 g/cm3，误差百分比为1．42%，说明模拟值与实验值相差很小。7个实验值与对应的模拟值之间的平均误差为0．055 g/cm3，误差百分数为2．98%，尽管实验值与模拟值略有偏差，但偏差很小，从一定程度上验证了模型的可靠性。从密度变化曲线可得出，在致密化前期，曲线的斜率较大，说明前期预制体致密化速率很快，密度从0．35 g/cm3增加到1 g/cm3以上;致密化中期，预制体致密化速率有所降低，50 h内密度仅增加0．4 g/cm3;致密化后期，曲线趋于平稳，50 h内密度增加仅为0．13 g/cm3，尤其在130 h以后，预制体密度基本不变，其致密化效率非常低。

图8为在预制体区内z=0高度沿r向取的3点(0，0)、(0，0．01)和(0，0．015)处密度的计算值随时间的变化曲线。

图8 预制体区内不同点处的密度值随时间的变化曲线Fig．8 Curve of density at different points in the preform as a function of time

图8显示，3点处密度的变化规律与预制体整体平均密度的变化规律一致，每个点的密度都是逐渐增加，最终保持不变，与实际的致密化过程吻合，说明整个耦合计算在这3点处的收敛性很好。密度分布云图则显示，计算在整个区域内的收敛性也很好，说明对3个控制方程的耦合求解是可行的，计算结果是可靠的。

5 结论

(1)不考虑化学反应，计算了整个区域的流场和温度场，模拟结果显示，自由流动区内气体流速比预制体区内的高出3个数量级，预制体区内流速基本为零;反应炉内温度场分布很均匀，预制体上下温度相差很小，炉壁和预制体中心温度相差3 K左右，认为沉积反应在等温条件下进行。

(2)确定沉积温度为1 393 K，对质量、动量守恒方程和孔隙率变化方程进行耦合计算，不同时间预制体密度分布云图说明，随沉积进行，最高密度区域从预制体中心向两侧移动，靠近窄缝的预制体区域，材料密度偏低。

(3)将整体平均密度的计算结果与实验值进行比较，结果相差很小，且计算值和实验值的变化趋势相同，且预制体区内沿r向所取的3点密度的变化曲线也是收敛的，验证了模型的可靠性。从密度变化曲线也表明，预制体在致密化开始50 h内，致密化速率很快，50 h以后，致密化速率逐渐减小，直到密度保持不变，这与实验过程也是相吻合的，说明建立的耦合模型可描述ICVI工艺过程。

［1］ 杨锦，李贺军，李克智，等．ICVI制备炭/炭复合材料流场模拟［J］．科学技术与工程，2008，8(11):2786-2791．

［2］ 向巧，罗瑞盈，章劲草．炭/炭复合材料等温CVI工艺计算机模拟的应用［J］．炭素技术，2009，28(1):40-43．

［3］ 姜开宇，李贺军，侯向辉，等．轴对称C/C复合材料件等温CVI过程的数值模拟研究［J］．西北工业大学学报，2000，18(9):665-668．

［4］ Zhang W G，Hu Z J，Huttinger K J．Chemical vapor infiltration of carbon fiber felt:optimization of densification and carbon microstructure［J］．Carbon，2002，40:2529-2545．

［5］ Li He-jun，Li Ai-jun，Bai Rui-cheng，et al．Numerical simulation of chemical vapor infiltration of propylene into C/C composites with reduced multi-step kinetic models［J］．Carbon，2005，43(14):2937-2950．

［6］ Li Ai-jun，Koyo Norinaga，Zhang Wei-gang，et al．Modeling and simulation of materials synthesis:Chemical vapor deposition and infiltration of pyrolytic carbon［J］．Composites Science and Technology，2008，68(5):1097-1104．

［7］ 孙国岭，李贺军，齐乐华，等．CVI工艺建模研究进展［J］．材料工程，2006(增刊1):441-444．

［8］ 焦妍琼，李贺军，李克智．TCVI工艺制备C/C复合材料的多物理场耦合模拟［J］．中国科学 E辑:技术科学．2009，52(11):3173-3179．

［9］ 边国军，齐乐华，宋永善，等．基于数值模拟的热梯度CVI制备C/C复合材料致密化行为［J］．复合材料学报，2011，28(4):29-33．

［10］ 李爱军．碳/碳复合材料想能预测与ICVI工艺系统虚拟设计［D］．西安:西北工业大学，2004．

［11］ 孙国岭，李贺军，齐乐华，等．C/C复合材料热梯度CVI工艺致密化过程的非稳态温度场分析［J］．金属学报，2006，42(10):1046-1050．

［12］ 时钧，汪家鼎，余国琮，等．化学工程手册［M］．北京:化学工业出版社，1996．

［13］ 侯向辉，李贺军，李克智，等．CVI制备碳/碳复合材料致密化行为模拟研究［J］．兵器材料科学与工程，1999，22(2):28-33．