王建武,李克智,張守陽,李 偉
(西北工業(yè)大學(xué)碳/碳復(fù)合材料工程技術(shù)研究中心,凝固技術(shù)國家重點(diǎn)實(shí)驗(yàn)室,西安 710072)
C/C復(fù)合材料不僅具有密度小、比模量高、比強(qiáng)度大、摩擦性能好等優(yōu)異的力學(xué)性能,而且還保持了炭材料所固有的耐高溫、耐熱沖擊、抗燒蝕、熱膨脹系數(shù)低等良好的熱物理性能,被廣泛應(yīng)用于航空航天領(lǐng)域[1-2]。等溫化學(xué)氣相滲透(ICVI)工藝是制備 C/C復(fù)合材料最主要的工藝,但在沉積過程中,為了避免制件表面過早結(jié)殼,通常在較低的溫度和壓力下進(jìn)行沉積,導(dǎo)致整個制備周期過長,生產(chǎn)成本很高。因此,尋求最優(yōu)的沉積工藝,實(shí)現(xiàn)C/C復(fù)合材料的快速沉積是實(shí)際生產(chǎn)中亟待解決的問題之一。
影響ICVI工藝的參數(shù)眾多,主要包括沉積溫度、壓強(qiáng)、氣體流量,若對各個工藝參數(shù)的影響規(guī)律設(shè)計(jì)實(shí)驗(yàn)進(jìn)行探索,需要消耗大量的時(shí)間和物力。通過計(jì)算機(jī)模擬技術(shù)對ICVI工藝過程進(jìn)行模擬,可研究沉積環(huán)境中的流場和溫度場分布,顯示預(yù)制體致密化的過程,對C/C復(fù)合材料ICVI工藝的研究具有十分重要的意義[2-3]。目前,ICVI工藝的計(jì)算機(jī)模擬受到越來越多的關(guān)注,包括預(yù)制體結(jié)構(gòu)建模和致密化過程仿真、計(jì)算域內(nèi)物理場的分布計(jì)算及前驅(qū)體反應(yīng)的化學(xué)反應(yīng)動力學(xué)研究。Huttinger等[4]認(rèn)為,烯烴、芳香烴和大分子烴是主要成炭的物質(zhì),以毛細(xì)管作為孔隙簡化模型,建立了以甲烷為前驅(qū)體的ICVI工藝模型,并進(jìn)行了實(shí)驗(yàn)驗(yàn)證。李愛軍[5-6]建立了雙孔隙模型和“平行-串式”反應(yīng)模型,分析了預(yù)制體致密化過程中孔隙拓?fù)浣Y(jié)構(gòu)的變化規(guī)律,以及均相反應(yīng)和非均相反應(yīng)的競爭關(guān)系與熱解炭織構(gòu)的關(guān)系。孫國嶺、趙建國等[7]直接運(yùn)用NS方程描述了前驅(qū)體在自由空間的流動,但并未給出預(yù)制體內(nèi)部流場的分布規(guī)律。焦妍瓊、邊國軍[8-9]等分別建立了2D軸對稱瞬態(tài)模型,利用熱傳導(dǎo)、對流和擴(kuò)散、傳質(zhì)、沉積反應(yīng)的多場耦合計(jì)算,模擬了熱梯度CVI中流場、溫度場的分布及密度的演變規(guī)律。
本文針對ICVI工藝,根據(jù)實(shí)驗(yàn)用的沉積裝置,建立簡化的幾何模型,依據(jù)質(zhì)量、動量和能量守恒定律建立各個區(qū)域的偏微分控制方程,計(jì)算沉積爐內(nèi)流場和溫度場的分布,分析前驅(qū)體進(jìn)入沉積爐后的流動速率以及升溫過程;分析整個沉積過程預(yù)制體密度演變規(guī)律,并對比整體平均密度隨時(shí)間變化的計(jì)算值和實(shí)驗(yàn)值,驗(yàn)證模型的可靠性。
ICVI的工藝實(shí)驗(yàn)裝置如圖1(a)所示。前驅(qū)體從底部通入,流經(jīng)窄縫區(qū)時(shí)通過擴(kuò)散方式進(jìn)入預(yù)制體,在預(yù)制體內(nèi)發(fā)生熱解反應(yīng)生成熱解炭,并沉積在纖維表面,反應(yīng)副產(chǎn)物隨通道中未反應(yīng)的氣體從上方的出口排出。對沉積裝置進(jìn)行簡化,建立了如圖1(b)所示的2D軸對稱幾何模型。其中,Ω1為預(yù)制體區(qū),Ω2為自由流動區(qū),Ω1,1、Ω1,2分別表示預(yù)制體內(nèi)、外側(cè)邊界,Ωin、Ωout分別為氣體的入口和出口,Ωw為爐壁。模擬采用的沉積工藝參數(shù)為沉積溫度為1 393 K,壓強(qiáng)為10 kPa,天然氣在入口處的流速為4 m/s;預(yù)制體初始孔隙率為 0.675,初始密度為 0.575 kg/m3。
另外,對模型作進(jìn)一步假設(shè):
(1)不考慮碳?xì)錃怏w熱解反應(yīng)對系統(tǒng)熱量的影響;
(2)反應(yīng)氣體為理想氣體,且不可壓縮,預(yù)制體為各向同性體;
(3)熱解反應(yīng)均為一級反應(yīng),熱解的中間產(chǎn)物僅考慮 C2H4、C2H2、C6H6、H24 種主要的組分,反應(yīng)采用簡化的多步反應(yīng)模型,如圖2所示。其中,k1~k4、K1~K4分別為氣相反應(yīng)和氣固反應(yīng)的反應(yīng)速率常數(shù),其值見表2,C∞表示熱解炭;
(4)預(yù)制體區(qū)的比表面積遠(yuǎn)大于自由流動區(qū)的比表面積,故僅考慮發(fā)生在預(yù)制體區(qū)內(nèi)的沉積反應(yīng)。
圖1 ICVI工藝原理圖及簡化后的2D軸對稱幾何模型Fig.1 diagram of ICVI furnace and simplified 2D axis-symmetric geometry model
圖2 簡化的甲烷熱解的多步化學(xué)反應(yīng)動力學(xué)模型Fig.2 Simplified multi-step chemical reaction kinetic model of methane pyrolysis
ICVI過程中爐壁溫度設(shè)為沉積溫度Tw,熱量以前驅(qū)體和預(yù)制體為媒介通過熱傳導(dǎo)和對流的方式從爐壁向內(nèi)部傳遞,根據(jù)能量守恒定律,分別建立Ω1區(qū)和Ω2區(qū)的熱傳導(dǎo)方程:
其中
式中 ρ、ρp分別為甲烷和預(yù)制體的密度,kg/m3;cp和cp,p分別為甲烷和預(yù)制體纖維的常壓比定壓熱容,J/(kg·K);k、kp分別為甲烷和預(yù)制體纖維的熱導(dǎo)率,W/(m·K);θp為預(yù)制體中纖維的體積分?jǐn)?shù);u為速度矢量。
爐壁的溫度為沉積溫度,出口處僅考慮熱對流,得到以下邊界條件和初始條件:
在流場和溫度場的耦合計(jì)算中,并未考慮甲烷的裂解。所以,僅需要計(jì)算甲烷的比定壓熱容和熱導(dǎo)率。按照多項(xiàng)式擬合,得到[11]:
預(yù)制體纖維的比定壓熱容和熱導(dǎo)率:
氣體在Ω1區(qū)和Ω2區(qū)內(nèi)的流動方式不同,根據(jù)動量守恒原理,分別建立如下的弱可壓縮氣體流場方程和連續(xù)方程[1,9]:
式中 ρ為混合氣體的密度,kg/m3;ε為預(yù)制體的孔隙率;η為混合氣體的粘度系數(shù),Pa·s;κ為預(yù)制體的滲透率,設(shè)為 0.25[9];p 為爐內(nèi)壓強(qiáng),Pa;I為單位矢量。
入口邊界和出口邊界分別設(shè)為速度和壓強(qiáng),邊界條件和初始條件如下:
按照混合法則,混合氣體密度計(jì)算方法如下:
式中 ρi為第 i組分氣體的密度,kg/m3;Ci為第 i組分氣體的濃度,mol/m3;n為氣體組分?jǐn)?shù)。
其中,純組分氣體的密度根據(jù)理想氣體狀態(tài)方程計(jì)算:
對于混合氣體的粘度,采用Sutherland模型進(jìn)行估算[11]:
式中 ηi和yi是純組分i的粘度和分子分率[12]。
式(8)中的關(guān)聯(lián)系數(shù)φij可由簡單的wilke近似式計(jì)算:
前驅(qū)體的溫度迅速升高到裂解溫度后,會發(fā)生分解,生成其他組分的氣體,見圖2?;旌蠚怏w通過擴(kuò)散方式向預(yù)制體內(nèi)傳輸?shù)耐瑫r(shí),發(fā)生熱解反應(yīng)生成熱解炭,沉積在纖維的表面。根據(jù)質(zhì)量守恒原理,對Ω1區(qū)和Ω2區(qū)分別建立濃度變化的偏微分方程:
式中 Ci為第i種反應(yīng)物的濃度,mol/m3;Dieff、Diab分別為第i種組分在Ω1區(qū)和Ω2區(qū)的擴(kuò)散系數(shù),m2/s;R1,i、R2,i分別為第 i種反應(yīng)物在 Ω1區(qū)和 Ω2區(qū)的反應(yīng)速率,mol/(m3·s),見表 1。
其邊界條件和初始條件為
綜合考慮Fick擴(kuò)散和Knudsen擴(kuò)散,得到有效擴(kuò)散系數(shù)表達(dá)式[11]:
式中 ε為預(yù)制體孔隙率;τ表示預(yù)制體內(nèi)孔隙結(jié)構(gòu)的曲折因子,取值一般為3~7。
Dik和Diab的值見表2。
表1 致密化過程中各組分在Ω1和Ω2區(qū)的反應(yīng)速率Table 1 Reaction rate of species in Ω1and Ω2during the densification process
表2 各組分的Fick擴(kuò)散系數(shù)和Knudsen擴(kuò)散系數(shù)的值及各反應(yīng)的反應(yīng)速率常數(shù)ki和KiTable 2 Value of and and reaction rate constant kiand Ki
表2 各組分的Fick擴(kuò)散系數(shù)和Knudsen擴(kuò)散系數(shù)的值及各反應(yīng)的反應(yīng)速率常數(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— —
前驅(qū)體及其裂解產(chǎn)物在預(yù)制體區(qū)發(fā)生熱解反應(yīng),生成熱解炭并填充在孔隙中,使孔隙率不斷降低,預(yù)制體孔隙結(jié)構(gòu)及比表面計(jì)算均采用“雙孔隙模型”,其變化過程可用如下方程表示:
式中 Mc、ρc分別為熱解炭的摩爾質(zhì)量和密度,kg/mol,kg/m3;R為生成熱解炭的反應(yīng)速率,也即熱解炭的沉積速率,mol/(m3·s);它與預(yù)制體的比表面積Sv有關(guān)。
設(shè)定所有的邊界條件為絕緣邊界,初始條件為
本文基于以上建立的氣體流動、熱量傳遞和孔隙率變化的偏微分方程,考慮到各個物理場之間的耦合關(guān)系,通過多物理場耦合軟件(COMSOL Muitiphysics),采用有限元方法進(jìn)行迭代耦合計(jì)算,具體的計(jì)算過程如圖3所示。
圖3 ICVI致密化過程模擬流程圖Fig.3 Flow diagram of densification simulation for ICVI process
圖4(a)顯示了t=5 min時(shí)預(yù)制體區(qū)和自由流動區(qū)內(nèi)速度場的分布??煽闯?,前驅(qū)體在Ω2中間位置的速率遠(yuǎn)大于Ω1區(qū)內(nèi)的速率,高出約3個數(shù)量級,且整個Ω1區(qū)域的速率接近于零。這是由于在預(yù)制體內(nèi)有大量的微孔,整個區(qū)域的比表面積達(dá)到了2×105m-1,氣態(tài)前驅(qū)體受到相當(dāng)大得粘滯阻力。因此,前驅(qū)體的流動受到很大的限制,流動速率接近于零。
為了進(jìn)一步了解在不同高度z處速度沿半徑r方向的變化情況,在 z= -0.0 1m,z=0和 z=0.0 1m 處各取一個橫截面,其速度-位移(v-r)曲線如圖4(b)所示。該圖顯示在不同高度z上速度沿徑向r的分布基本一致,即Ω2區(qū)為未完全發(fā)展的層流,窄縫中心部位流動速度較快,往爐壁一側(cè)速度逐漸降為零,在靠近預(yù)制體一側(cè)速度降低到0.5 m/s左右;Ω1區(qū)內(nèi)速度接近于零,但并不為零。不同z值的3條曲線也說明,爐內(nèi)流場的分布隨高度的增加其變化并不顯著,流場分布較為穩(wěn)定。從曲線中還可看出,在Ω1區(qū)與Ω2區(qū)的交界面Ω1,2處流速的變化是連續(xù)的,這與實(shí)際流動是相符的。
圖4 t=5 min時(shí)沉積爐內(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
圖5給出了t=5 min時(shí)整個區(qū)域內(nèi)溫度場的分布,以及不同高度z處溫度沿徑向r的變化曲線。從圖5(a)可看出,入口附近前驅(qū)體的溫度還較低,但前驅(qū)體與壁面經(jīng)過短時(shí)間的熱量交換后,其溫度迅速升高,已接近設(shè)定的沉積溫度,而且整個沉積爐內(nèi)的溫度分布基本一致,符合等溫沉積的條件。圖5(b)顯示從下到上氣體溫度逐漸升高,這是由升溫時(shí)間不同造成的,但上下溫度差在2 K以內(nèi)。z=-0.01 m處的曲線在靠近界面位置(0.03 m<r<0.04 m)溫度有下降趨勢。這是由于氣體受熱時(shí)間短,未被充分加熱便傳輸?shù)筋A(yù)制體內(nèi),且在同高度的Ω2區(qū)相對較冷的氣體在z方向的快速流動還帶走了部分熱量;z=0和z=0.01 m位置的溫度都隨著r增加而升高,因?yàn)樵娇拷鼱t壁熱量傳輸?shù)木嚯x越短,溫度也越高。此外,通過爐壁溫度和預(yù)制體中心溫度的比較可看出,不同高度z位置內(nèi)外的溫度差都在3 K左右,溫度相差都很小,進(jìn)一步驗(yàn)證了前驅(qū)體的沉積過程是在等溫條件下進(jìn)行的。
圖5 t=5 min時(shí)沉積爐內(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
圖6給出了沉積溫度為1 393 K,壓強(qiáng)為10 kPa,天然氣在入口處的流速為4 m/s時(shí),預(yù)制體內(nèi)密度分布的演變規(guī)律。致密化前期0~50 h,預(yù)制體中心位置(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,預(yù)制體內(nèi)密度的分布規(guī)律與前期相同,但密度最大的區(qū)域向外側(cè)移動,靠近窄縫的密度較小的區(qū)域也在逐漸縮小,見圖6(b)、(c)。致密化后期100~150 h,預(yù)制體內(nèi)整體的密度分布較為均勻,只有在靠近窄縫附近的一個狹窄區(qū)域內(nèi)密度偏低,見圖6(d)。
圖6 不同沉積時(shí)間預(yù)制體內(nèi)的密度分布Fig.6 Density distribution of the preform at different deposition time
圖6中還顯示,在不同沉積階段,靠近窄縫位置的預(yù)制體密度總是很低。這是因?yàn)榭拷p的區(qū)域,氣體在z向的流動比較強(qiáng),將未來得及反應(yīng)的氣體帶走,反應(yīng)的氣體在該區(qū)域內(nèi)的滯留時(shí)間很短。所以,致密化的效果較其他位置較差。
從密度最大值的變化也可反映出整個預(yù)制體致密化的進(jìn)程。致密化前期0~50 h,最大密度增加為0.945 9 g/cm3,增加近一倍;致密化中后期50~100 h和100~150 h,最大密度增加分別為0.408 1 g/cm3和0.082 4 g/cm3,增加幅度明顯降低。由此可很明顯地看出,預(yù)制體密度的增加主要發(fā)生在致密化前期,占整個密度增加的66%,越到沉積后期,密度增加越困難,隨著沉積的不斷進(jìn)行,致密化效率也在不斷降低。這是由于沉積前期,預(yù)制體內(nèi)孔隙率和比表面積都很大,前驅(qū)體及其裂解氣體不僅可在預(yù)制體內(nèi)很好地進(jìn)行傳遞,及時(shí)補(bǔ)充因沉積而消耗的氣體,而且熱解炭的沉積反應(yīng)也比較劇烈;沉積進(jìn)行到后期,孔隙已經(jīng)被大量的熱解炭填充,氣體傳輸?shù)耐ǖ辣欢氯?,沉積的氣體得不到及時(shí)補(bǔ)充,沉積的速率明顯降低。
圖7給出了致密化150 h過程中預(yù)制體整體平均密度隨時(shí)間的變化曲線及不同致密化時(shí)間預(yù)制體密度的實(shí)驗(yàn)值。
圖7 致密化150 h預(yù)制體整體平均密度隨沉積時(shí)間變化的模擬曲線和實(shí)驗(yàn)結(jié)果比較Fig.7 Comparison of average density of preform between simulated and experimental results after 150h densification
沉積140 h,預(yù)制體平均密度的計(jì)算值為1.871 g/cm3,實(shí)驗(yàn)值為 1.845 g/cm3,誤差為 0.026 3 g/cm3,誤差百分比為1.42%,說明模擬值與實(shí)驗(yàn)值相差很小。7個實(shí)驗(yàn)值與對應(yīng)的模擬值之間的平均誤差為0.055 g/cm3,誤差百分?jǐn)?shù)為2.98%,盡管實(shí)驗(yàn)值與模擬值略有偏差,但偏差很小,從一定程度上驗(yàn)證了模型的可靠性。從密度變化曲線可得出,在致密化前期,曲線的斜率較大,說明前期預(yù)制體致密化速率很快,密度從0.35 g/cm3增加到1 g/cm3以上;致密化中期,預(yù)制體致密化速率有所降低,50 h內(nèi)密度僅增加0.4 g/cm3;致密化后期,曲線趨于平穩(wěn),50 h內(nèi)密度增加僅為0.13 g/cm3,尤其在130 h以后,預(yù)制體密度基本不變,其致密化效率非常低。
圖8為在預(yù)制體區(qū)內(nèi)z=0高度沿r向取的3點(diǎn)(0,0)、(0,0.01)和(0,0.015)處密度的計(jì)算值隨時(shí)間的變化曲線。
圖8 預(yù)制體區(qū)內(nèi)不同點(diǎn)處的密度值隨時(shí)間的變化曲線Fig.8 Curve of density at different points in the preform as a function of time
圖8顯示,3點(diǎn)處密度的變化規(guī)律與預(yù)制體整體平均密度的變化規(guī)律一致,每個點(diǎn)的密度都是逐漸增加,最終保持不變,與實(shí)際的致密化過程吻合,說明整個耦合計(jì)算在這3點(diǎn)處的收斂性很好。密度分布云圖則顯示,計(jì)算在整個區(qū)域內(nèi)的收斂性也很好,說明對3個控制方程的耦合求解是可行的,計(jì)算結(jié)果是可靠的。
(1)不考慮化學(xué)反應(yīng),計(jì)算了整個區(qū)域的流場和溫度場,模擬結(jié)果顯示,自由流動區(qū)內(nèi)氣體流速比預(yù)制體區(qū)內(nèi)的高出3個數(shù)量級,預(yù)制體區(qū)內(nèi)流速基本為零;反應(yīng)爐內(nèi)溫度場分布很均勻,預(yù)制體上下溫度相差很小,爐壁和預(yù)制體中心溫度相差3 K左右,認(rèn)為沉積反應(yīng)在等溫條件下進(jìn)行。
(2)確定沉積溫度為1 393 K,對質(zhì)量、動量守恒方程和孔隙率變化方程進(jìn)行耦合計(jì)算,不同時(shí)間預(yù)制體密度分布云圖說明,隨沉積進(jìn)行,最高密度區(qū)域從預(yù)制體中心向兩側(cè)移動,靠近窄縫的預(yù)制體區(qū)域,材料密度偏低。
(3)將整體平均密度的計(jì)算結(jié)果與實(shí)驗(yàn)值進(jìn)行比較,結(jié)果相差很小,且計(jì)算值和實(shí)驗(yàn)值的變化趨勢相同,且預(yù)制體區(qū)內(nèi)沿r向所取的3點(diǎn)密度的變化曲線也是收斂的,驗(yàn)證了模型的可靠性。從密度變化曲線也表明,預(yù)制體在致密化開始50 h內(nèi),致密化速率很快,50 h以后,致密化速率逐漸減小,直到密度保持不變,這與實(shí)驗(yàn)過程也是相吻合的,說明建立的耦合模型可描述ICVI工藝過程。
[1] 楊錦,李賀軍,李克智,等.ICVI制備炭/炭復(fù)合材料流場模擬[J].科學(xué)技術(shù)與工程,2008,8(11):2786-2791.
[2] 向巧,羅瑞盈,章勁草.炭/炭復(fù)合材料等溫CVI工藝計(jì)算機(jī)模擬的應(yīng)用[J].炭素技術(shù),2009,28(1):40-43.
[3] 姜開宇,李賀軍,侯向輝,等.軸對稱C/C復(fù)合材料件等溫CVI過程的數(shù)值模擬研究[J].西北工業(yè)大學(xué)學(xué)報(bào),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ìn)展[J].材料工程,2006(增刊1):441-444.
[8] 焦妍瓊,李賀軍,李克智.TCVI工藝制備C/C復(fù)合材料的多物理場耦合模擬[J].中國科學(xué) E輯:技術(shù)科學(xué).2009,52(11):3173-3179.
[9] 邊國軍,齊樂華,宋永善,等.基于數(shù)值模擬的熱梯度CVI制備C/C復(fù)合材料致密化行為[J].復(fù)合材料學(xué)報(bào),2011,28(4):29-33.
[10] 李愛軍.碳/碳復(fù)合材料想能預(yù)測與ICVI工藝系統(tǒng)虛擬設(shè)計(jì)[D].西安:西北工業(yè)大學(xué),2004.
[11] 孫國嶺,李賀軍,齊樂華,等.C/C復(fù)合材料熱梯度CVI工藝致密化過程的非穩(wěn)態(tài)溫度場分析[J].金屬學(xué)報(bào),2006,42(10):1046-1050.
[12] 時(shí)鈞,汪家鼎,余國琮,等.化學(xué)工程手冊[M].北京:化學(xué)工業(yè)出版社,1996.
[13] 侯向輝,李賀軍,李克智,等.CVI制備碳/碳復(fù)合材料致密化行為模擬研究[J].兵器材料科學(xué)與工程,1999,22(2):28-33.