ICVI工藝致密化C/C復合材料的數(shù)值模擬①

王建武，李克智，張守陽，李 偉

(西北工業(yè)大學碳/碳復合材料工程技術研究中心，凝固技術國家重點實驗室，西安 710072)

0 引言

C/C復合材料不僅具有密度小、比模量高、比強度大、摩擦性能好等優(yōu)異的力學性能，而且還保持了炭材料所固有的耐高溫、耐熱沖擊、抗燒蝕、熱膨脹系數(shù)低等良好的熱物理性能，被廣泛應用于航空航天領域［1－2］。等溫化學氣相滲透(ICVI)工藝是制備 C/C復合材料最主要的工藝，但在沉積過程中，為了避免制件表面過早結殼，通常在較低的溫度和壓力下進行沉積，導致整個制備周期過長，生產(chǎn)成本很高。因此，尋求最優(yōu)的沉積工藝，實現(xiàn)C/C復合材料的快速沉積是實際生產(chǎn)中亟待解決的問題之一。

影響ICVI工藝的參數(shù)眾多，主要包括沉積溫度、壓強、氣體流量，若對各個工藝參數(shù)的影響規(guī)律設計實驗進行探索，需要消耗大量的時間和物力。通過計算機模擬技術對ICVI工藝過程進行模擬，可研究沉積環(huán)境中的流場和溫度場分布，顯示預制體致密化的過程，對C/C復合材料ICVI工藝的研究具有十分重要的意義［2－3］。目前，ICVI工藝的計算機模擬受到越來越多的關注，包括預制體結構建模和致密化過程仿真、計算域內(nèi)物理場的分布計算及前驅體反應的化學反應動力學研究。Huttinger等［4］認為，烯烴、芳香烴和大分子烴是主要成炭的物質，以毛細管作為孔隙簡化模型，建立了以甲烷為前驅體的ICVI工藝模型，并進行了實驗驗證。李愛軍［5－6］建立了雙孔隙模型和“平行-串式”反應模型，分析了預制體致密化過程中孔隙拓撲結構的變化規(guī)律，以及均相反應和非均相反應的競爭關系與熱解炭織構的關系。孫國嶺、趙建國等［7］直接運用NS方程描述了前驅體在自由空間的流動，但并未給出預制體內(nèi)部流場的分布規(guī)律。焦妍瓊、邊國軍［8－9］等分別建立了2D軸對稱瞬態(tài)模型，利用熱傳導、對流和擴散、傳質、沉積反應的多場耦合計算，模擬了熱梯度CVI中流場、溫度場的分布及密度的演變規(guī)律。

本文針對ICVI工藝，根據(jù)實驗用的沉積裝置，建立簡化的幾何模型，依據(jù)質量、動量和能量守恒定律建立各個區(qū)域的偏微分控制方程，計算沉積爐內(nèi)流場和溫度場的分布，分析前驅體進入沉積爐后的流動速率以及升溫過程;分析整個沉積過程預制體密度演變規(guī)律，并對比整體平均密度隨時間變化的計算值和實驗值，驗證模型的可靠性。

1 幾何模型建立和模型假設

ICVI的工藝實驗裝置如圖1(a)所示。前驅體從底部通入，流經(jīng)窄縫區(qū)時通過擴散方式進入預制體，在預制體內(nèi)發(fā)生熱解反應生成熱解炭，并沉積在纖維表面，反應副產(chǎn)物隨通道中未反應的氣體從上方的出口排出。對沉積裝置進行簡化，建立了如圖1(b)所示的2D軸對稱幾何模型。其中，Ω1為預制體區(qū)，Ω2為自由流動區(qū)，Ω1，1、Ω1，2分別表示預制體內(nèi)、外側邊界，Ωin、Ωout分別為氣體的入口和出口，Ωw為爐壁。模擬采用的沉積工藝參數(shù)為沉積溫度為1 393 K，壓強為10 kPa，天然氣在入口處的流速為4 m/s;預制體初始孔隙率為 0．675，初始密度為 0．575 kg/m3。

另外，對模型作進一步假設:

(1)不考慮碳氫氣體熱解反應對系統(tǒng)熱量的影響;

(2)反應氣體為理想氣體，且不可壓縮，預制體為各向同性體;

(3)熱解反應均為一級反應，熱解的中間產(chǎn)物僅考慮 C2H4、C2H2、C6H6、H24 種主要的組分，反應采用簡化的多步反應模型，如圖2所示。其中，k1～k4、K1～K4分別為氣相反應和氣固反應的反應速率常數(shù)，其值見表2，C∞表示熱解炭;

(4)預制體區(qū)的比表面積遠大于自由流動區(qū)的比表面積，故僅考慮發(fā)生在預制體區(qū)內(nèi)的沉積反應。

圖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 數(shù)學模型的建立

2．1 傳熱方程

ICVI過程中爐壁溫度設為沉積溫度Tw，熱量以前驅體和預制體為媒介通過熱傳導和對流的方式從爐壁向內(nèi)部傳遞，根據(jù)能量守恒定律，分別建立Ω1區(qū)和Ω2區(qū)的熱傳導方程:

其中

式中 ρ、ρp分別為甲烷和預制體的密度，kg/m3;cp和cp，p分別為甲烷和預制體纖維的常壓比定壓熱容，J/(kg·K);k、kp分別為甲烷和預制體纖維的熱導率，W/(m·K);θp為預制體中纖維的體積分數(shù);u為速度矢量。

爐壁的溫度為沉積溫度，出口處僅考慮熱對流，得到以下邊界條件和初始條件:

在流場和溫度場的耦合計算中，并未考慮甲烷的裂解。所以，僅需要計算甲烷的比定壓熱容和熱導率。按照多項式擬合，得到［11］:

預制體纖維的比定壓熱容和熱導率:

2．2 流動方程

氣體在Ω1區(qū)和Ω2區(qū)內(nèi)的流動方式不同，根據(jù)動量守恒原理，分別建立如下的弱可壓縮氣體流場方程和連續(xù)方程［1，9］:

式中 ρ為混合氣體的密度，kg/m3;ε為預制體的孔隙率;η為混合氣體的粘度系數(shù)，Pa·s;κ為預制體的滲透率，設為 0．25［9］;p 為爐內(nèi)壓強，Pa;I為單位矢量。

入口邊界和出口邊界分別設為速度和壓強，邊界條件和初始條件如下:

按照混合法則，混合氣體密度計算方法如下:

式中 ρi為第 i組分氣體的密度，kg/m3;Ci為第 i組分氣體的濃度，mol/m3;n為氣體組分數(shù)。

其中，純組分氣體的密度根據(jù)理想氣體狀態(tài)方程計算:

對于混合氣體的粘度，采用Sutherland模型進行估算［11］:

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

式(8)中的關聯(lián)系數(shù)φij可由簡單的wilke近似式計算:

2．3 傳質方程

前驅體的溫度迅速升高到裂解溫度后，會發(fā)生分解，生成其他組分的氣體，見圖2。混合氣體通過擴散方式向預制體內(nèi)傳輸?shù)耐瑫r，發(fā)生熱解反應生成熱解炭，沉積在纖維的表面。根據(jù)質量守恒原理，對Ω1區(qū)和Ω2區(qū)分別建立濃度變化的偏微分方程:

式中 Ci為第i種反應物的濃度，mol/m3;Dieff、Diab分別為第i種組分在Ω1區(qū)和Ω2區(qū)的擴散系數(shù)，m2/s;R1，i、R2，i分別為第 i種反應物在 Ω1區(qū)和 Ω2區(qū)的反應速率，mol/(m3·s)，見表 1。

其邊界條件和初始條件為

綜合考慮Fick擴散和Knudsen擴散，得到有效擴散系數(shù)表達式［11］:

式中 ε為預制體孔隙率;τ表示預制體內(nèi)孔隙結構的曲折因子，取值一般為3～7。

Dik和Diab的值見表2。

表1 致密化過程中各組分在Ω1和Ω2區(qū)的反應速率Table 1 Reaction rate of species in Ω1and Ω2during the densification process

表2 各組分的Fick擴散系數(shù)和Knudsen擴散系數(shù)的值及各反應的反應速率常數(shù)ki和KiTable 2 Value of and and reaction rate constant kiand Ki

表2 各組分的Fick擴散系數(shù)和Knudsen擴散系數(shù)的值及各反應的反應速率常數(shù)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 孔隙率的變化方程

前驅體及其裂解產(chǎn)物在預制體區(qū)發(fā)生熱解反應，生成熱解炭并填充在孔隙中，使孔隙率不斷降低，預制體孔隙結構及比表面計算均采用“雙孔隙模型”，其變化過程可用如下方程表示:

式中 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時預制體區(qū)和自由流動區(qū)內(nèi)速度場的分布?？煽闯觯膀岓w在Ω2中間位置的速率遠大于Ω1區(qū)內(nèi)的速率，高出約3個數(shù)量級，且整個Ω1區(qū)域的速率接近于零。這是由于在預制體內(nèi)有大量的微孔，整個區(qū)域的比表面積達到了2×105m－1，氣態(tài)前驅體受到相當大得粘滯阻力。因此，前驅體的流動受到很大的限制，流動速率接近于零。

為了進一步了解在不同高度z處速度沿半徑r方向的變化情況，在 z= －0．0 1m，z=0和 z=0．0 1m 處各取一個橫截面，其速度-位移(v-r)曲線如圖4(b)所示。該圖顯示在不同高度z上速度沿徑向r的分布基本一致，即Ω2區(qū)為未完全發(fā)展的層流，窄縫中心部位流動速度較快，往爐壁一側速度逐漸降為零，在靠近預制體一側速度降低到0．5 m/s左右;Ω1區(qū)內(nèi)速度接近于零，但并不為零。不同z值的3條曲線也說明，爐內(nèi)流場的分布隨高度的增加其變化并不顯著，流場分布較為穩(wěn)定。從曲線中還可看出，在Ω1區(qū)與Ω2區(qū)的交界面Ω1，2處流速的變化是連續(xù)的，這與實際流動是相符的。

圖4 t=5 min時沉積爐內(nèi)速度場的分布及不同高度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時整個區(qū)域內(nèi)溫度場的分布，以及不同高度z處溫度沿徑向r的變化曲線。從圖5(a)可看出，入口附近前驅體的溫度還較低，但前驅體與壁面經(jīng)過短時間的熱量交換后，其溫度迅速升高，已接近設定的沉積溫度，而且整個沉積爐內(nèi)的溫度分布基本一致，符合等溫沉積的條件。圖5(b)顯示從下到上氣體溫度逐漸升高，這是由升溫時間不同造成的，但上下溫度差在2 K以內(nèi)。z=－0．01 m處的曲線在靠近界面位置(0．03 m＜r＜0．04 m)溫度有下降趨勢。這是由于氣體受熱時間短，未被充分加熱便傳輸?shù)筋A制體內(nèi)，且在同高度的Ω2區(qū)相對較冷的氣體在z方向的快速流動還帶走了部分熱量;z=0和z=0．01 m位置的溫度都隨著r增加而升高，因為越靠近爐壁熱量傳輸?shù)木嚯x越短，溫度也越高。此外，通過爐壁溫度和預制體中心溫度的比較可看出，不同高度z位置內(nèi)外的溫度差都在3 K左右，溫度相差都很小，進一步驗證了前驅體的沉積過程是在等溫條件下進行的。

圖5 t=5 min時沉積爐內(nèi)溫度場的分布及不同高度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時，預制體內(nèi)密度分布的演變規(guī)律。致密化前期0～50 h，預制體中心位置(0～0．02 m之間)，熱解炭的沉積很快，密度最大的區(qū)域出現(xiàn)在0．02～0．03 m之間，而在靠近窄縫的區(qū)域內(nèi)(0．03～0．04 m 之間)，熱解炭的沉積相對較少，如圖6(a)。致密化中期50～100 h，預制體內(nèi)密度的分布規(guī)律與前期相同，但密度最大的區(qū)域向外側移動，靠近窄縫的密度較小的區(qū)域也在逐漸縮小，見圖6(b)、(c)。致密化后期100～150 h，預制體內(nèi)整體的密度分布較為均勻，只有在靠近窄縫附近的一個狹窄區(qū)域內(nèi)密度偏低，見圖6(d)。

圖6 不同沉積時間預制體內(nèi)的密度分布Fig．6 Density distribution of the preform at different deposition time

圖6中還顯示，在不同沉積階段，靠近窄縫位置的預制體密度總是很低。這是因為靠近窄縫的區(qū)域，氣體在z向的流動比較強，將未來得及反應的氣體帶走，反應的氣體在該區(qū)域內(nèi)的滯留時間很短。所以，致密化的效果較其他位置較差。

從密度最大值的變化也可反映出整個預制體致密化的進程。致密化前期0～50 h，最大密度增加為0．945 9 g/cm3，增加近一倍;致密化中后期50～100 h和100～150 h，最大密度增加分別為0．408 1 g/cm3和0．082 4 g/cm3，增加幅度明顯降低。由此可很明顯地看出，預制體密度的增加主要發(fā)生在致密化前期，占整個密度增加的66%，越到沉積后期，密度增加越困難，隨著沉積的不斷進行，致密化效率也在不斷降低。這是由于沉積前期，預制體內(nèi)孔隙率和比表面積都很大，前驅體及其裂解氣體不僅可在預制體內(nèi)很好地進行傳遞，及時補充因沉積而消耗的氣體，而且熱解炭的沉積反應也比較劇烈;沉積進行到后期，孔隙已經(jīng)被大量的熱解炭填充，氣體傳輸?shù)耐ǖ辣欢氯?，沉積的氣體得不到及時補充，沉積的速率明顯降低。

圖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，誤差百分數(shù)為2．98%，盡管實驗值與模擬值略有偏差，但偏差很小，從一定程度上驗證了模型的可靠性。從密度變化曲線可得出，在致密化前期，曲線的斜率較大，說明前期預制體致密化速率很快，密度從0．35 g/cm3增加到1 g/cm3以上;致密化中期，預制體致密化速率有所降低，50 h內(nèi)密度僅增加0．4 g/cm3;致密化后期，曲線趨于平穩(wěn)，50 h內(nèi)密度增加僅為0．13 g/cm3，尤其在130 h以后，預制體密度基本不變，其致密化效率非常低。

圖8為在預制體區(qū)內(nèi)z=0高度沿r向取的3點(0，0)、(0，0．01)和(0，0．015)處密度的計算值隨時間的變化曲線。

圖8 預制體區(qū)內(nèi)不同點處的密度值隨時間的變化曲線Fig．8 Curve of density at different points in the preform as a function of time

圖8顯示，3點處密度的變化規(guī)律與預制體整體平均密度的變化規(guī)律一致，每個點的密度都是逐漸增加，最終保持不變，與實際的致密化過程吻合，說明整個耦合計算在這3點處的收斂性很好。密度分布云圖則顯示，計算在整個區(qū)域內(nèi)的收斂性也很好，說明對3個控制方程的耦合求解是可行的，計算結果是可靠的。

5 結論

(1)不考慮化學反應，計算了整個區(qū)域的流場和溫度場，模擬結果顯示，自由流動區(qū)內(nèi)氣體流速比預制體區(qū)內(nèi)的高出3個數(shù)量級，預制體區(qū)內(nèi)流速基本為零;反應爐內(nèi)溫度場分布很均勻，預制體上下溫度相差很小，爐壁和預制體中心溫度相差3 K左右，認為沉積反應在等溫條件下進行。

(2)確定沉積溫度為1 393 K，對質量、動量守恒方程和孔隙率變化方程進行耦合計算，不同時間預制體密度分布云圖說明，隨沉積進行，最高密度區(qū)域從預制體中心向兩側移動，靠近窄縫的預制體區(qū)域，材料密度偏低。

(3)將整體平均密度的計算結果與實驗值進行比較，結果相差很小，且計算值和實驗值的變化趨勢相同，且預制體區(qū)內(nèi)沿r向所取的3點密度的變化曲線也是收斂的，驗證了模型的可靠性。從密度變化曲線也表明，預制體在致密化開始50 h內(nèi)，致密化速率很快，50 h以后，致密化速率逐漸減小，直到密度保持不變，這與實驗過程也是相吻合的，說明建立的耦合模型可描述ICVI工藝過程。

［1］ 楊錦，李賀軍，李克智，等．ICVI制備炭/炭復合材料流場模擬［J］．科學技術與工程，2008，8(11):2786-2791．

［2］ 向巧，羅瑞盈，章勁草．炭/炭復合材料等溫CVI工藝計算機模擬的應用［J］．炭素技術，2009，28(1):40-43．

［3］ 姜開宇，李賀軍，侯向輝，等．軸對稱C/C復合材料件等溫CVI過程的數(shù)值模擬研究［J］．西北工業(yè)大學學報，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］ 邊國軍，齊樂華，宋永善，等．基于數(shù)值模擬的熱梯度CVI制備C/C復合材料致密化行為［J］．復合材料學報，2011，28(4):29-33．

［10］ 李愛軍．碳/碳復合材料想能預測與ICVI工藝系統(tǒng)虛擬設計［D］．西安:西北工業(yè)大學，2004．

［11］ 孫國嶺，李賀軍，齊樂華，等．C/C復合材料熱梯度CVI工藝致密化過程的非穩(wěn)態(tài)溫度場分析［J］．金屬學報，2006，42(10):1046-1050．

［12］ 時鈞，汪家鼎，余國琮，等．化學工程手冊［M］．北京:化學工業(yè)出版社，1996．

［13］ 侯向輝，李賀軍，李克智，等．CVI制備碳/碳復合材料致密化行為模擬研究［J］．兵器材料科學與工程，1999，22(2):28-33．