宋永善，齊樂(lè)華，張守陽(yáng)，張佳平，李逸仙，李賀軍

(1.西北工業(yè)大學(xué) 機(jī)電學(xué)院，西安 710072；2.西北工業(yè)大學(xué) C/C復(fù)合材料研究中心，西安 710072)

0 引言

C/C復(fù)合材料具有一系列優(yōu)異的性能[1-3]，如高強(qiáng)度、高模量，以及其強(qiáng)度隨溫度升高不降反升，此外還具有良好的抗熱沖擊、抗燒蝕、熱膨脹系數(shù)低及良好的導(dǎo)熱性能，是迄今為止惰性氣氛下最為理想的高溫結(jié)構(gòu)材料之一，被廣泛用于固體火箭發(fā)動(dòng)機(jī)喉襯、飛行器頭錐等航空航天結(jié)構(gòu)部件。然而，在高溫、高速?zèng)_擊等極端苛刻的服役環(huán)境下，C/C復(fù)合材料不可避免地會(huì)發(fā)生燒蝕現(xiàn)象，進(jìn)而導(dǎo)致材料性能急劇下降。

目前，關(guān)于C/C復(fù)合材料的燒蝕問(wèn)題已經(jīng)開展了很多實(shí)驗(yàn)研究，但受燒蝕過(guò)程中的高溫限制，實(shí)驗(yàn)研究所能測(cè)量的參數(shù)只有線燒蝕率和質(zhì)量燒蝕率，燒蝕機(jī)理的研究也大多基于對(duì)燒蝕后材料表面的微觀形貌進(jìn)行分析，得到的信息十分有限。數(shù)值模擬因此成為進(jìn)一步深入研究C/C復(fù)合材料燒蝕機(jī)理不可或缺的手段。Piyush Thakre等[4]和Daniele Bianchi等[5]基于不同假設(shè)分別計(jì)算了C/C復(fù)合材料發(fā)動(dòng)機(jī)噴管的熱化學(xué)燒蝕過(guò)程。王臣等[6]通過(guò)對(duì)碳基復(fù)合材料熱化學(xué)燒蝕機(jī)理分析，根據(jù)質(zhì)量守恒、能量守恒及系統(tǒng)化學(xué)反應(yīng)的熱化學(xué)平衡原理，建立了相應(yīng)的熱化學(xué)和熱力學(xué)燒蝕模型，模擬了材料內(nèi)部熱應(yīng)力場(chǎng)分布。楊德軍等[7]采用虛擬失效、重新構(gòu)建網(wǎng)格部件的方法實(shí)現(xiàn)燒蝕表面的退縮，建立了燒蝕表面退縮下瞬態(tài)溫度場(chǎng)的有限元模型，分析了熱化學(xué)燒蝕、燒蝕表面退縮及溫度場(chǎng)耦合作用下 C/C復(fù)合材料燒蝕性能的變化規(guī)律。解惠貞等[8]研究表明，纖維分布方向?qū)/C復(fù)合材料的燒蝕性能有重要影響，垂直燃?xì)饬鞣较虻睦w維與平行燃?xì)饬鞣较虻睦w維的線燒蝕率和質(zhì)量燒蝕率區(qū)別非常大。汪海斌等[9]針對(duì)軸編C/C復(fù)合材料的結(jié)構(gòu)形式和燒蝕機(jī)理，建立了噴管喉襯燒蝕的多尺度分析方法。通過(guò)宏觀-微觀的漸進(jìn)分析，獲得了噴管喉襯的燒蝕率和燒蝕形貌。王臣[10]建立了C/C復(fù)合材料的表面氧化模型來(lái)研究材料的燒蝕過(guò)程和熱環(huán)境參數(shù)，并與實(shí)驗(yàn)結(jié)果吻合良好。上述研究主要針對(duì)C/C復(fù)合材料的宏觀燒蝕機(jī)理分析，未考慮C/C復(fù)合材料中炭纖維和碳基體在密度和性能上的差異，因而很難揭示微尺度下C/C復(fù)合材料燒蝕過(guò)程中的形貌變化。實(shí)際上，C/C復(fù)合材料微尺度下的形貌變化會(huì)導(dǎo)致其宏觀尺度的表面粗糙化，進(jìn)而加劇材料表面的傳熱、傳質(zhì)，影響材料表面的流場(chǎng)。Lachaud J等[11]通過(guò)理論計(jì)算結(jié)合實(shí)驗(yàn)測(cè)定了C/C復(fù)合材料及其組分(纖維、基體)的固有氧化速率，建立了穩(wěn)態(tài)下C/C復(fù)合材料氧化模型，推導(dǎo)出的函數(shù)能直接計(jì)算氧化達(dá)到穩(wěn)態(tài)時(shí)纖維的形貌，但未能給出燒蝕過(guò)程中的微觀形貌演變規(guī)律。

本文建立微觀尺度燒蝕過(guò)程的形貌演變瞬態(tài)模型，研究C/C復(fù)合材料微觀燒蝕過(guò)程中纖維和基體的相對(duì)質(zhì)量損失規(guī)律，以期為研究C/C復(fù)合材料微觀燒蝕機(jī)理提供借鑒。

1 燒蝕過(guò)程建模

1.1 基本假設(shè)

在承受高溫高速?zèng)_刷的苛刻服役環(huán)境中，C/C復(fù)合材料會(huì)因氧化、升華、等離子腐蝕和機(jī)械剝離而被逐漸破壞，導(dǎo)致不能忽略的表面退移，該現(xiàn)象被稱為燒蝕。同時(shí)，燒蝕還會(huì)產(chǎn)生其他作用[5]：(1)增大反應(yīng)表面，進(jìn)一步加劇發(fā)生在材料表面的傳熱和傳質(zhì)；(2)表面粗糙化導(dǎo)致C/C復(fù)合材料的有效燒蝕速率不能簡(jiǎn)單地用其組分的燒蝕速率算術(shù)平均代替；(3)加劇機(jī)械剝離；(4)導(dǎo)致層流向湍流轉(zhuǎn)變。

但材料燒蝕率較低，對(duì)氣體流動(dòng)的連續(xù)性影響較小，故本文對(duì)燒蝕模擬過(guò)程做了簡(jiǎn)化，做出如下假設(shè)：(1)氣源位于離初始材料表面足夠遠(yuǎn)的一個(gè)距離；氣體流動(dòng)狀態(tài)為層流；(2)氧化反應(yīng)為一階反應(yīng)。則材料表面的氧化速率可定義為

v=-kCVmol

(1)

式中v為材料的氧化速率，m/s；k為氣體反應(yīng)速率常數(shù)，m/s；C為局部反應(yīng)氣體濃度，mol/m3；Vmol為材料的摩爾體積，m3/mol。

1.2 水平集法

本文采用水平集法來(lái)追蹤燒蝕過(guò)程中材料的表面移動(dòng)。材料表面由水平集函數(shù)定義：

(2)

其中，d為到材料表面的最短距離，φ=0代表材料表面。在定義了φ的初值及合適的速度場(chǎng)v后，可以通過(guò)計(jì)算水平集方程(3)來(lái)的追蹤材料的表面移動(dòng)。

(3)

但燒蝕速度v只在材料表面有物理意義，針對(duì)這種情況，必須通過(guò)外推構(gòu)造合適的速度場(chǎng)。本文選用方程(4)所示方法[12]進(jìn)行外推：

(4)

式中φ為速度v的任意一個(gè)分量。

方程(4)中符號(hào)函數(shù)S(φ)和單位法向量n分別定義為

(5)

(6)

在計(jì)算方程(3)一個(gè)步長(zhǎng)后，φ不再為符號(hào)距離函數(shù)，必須重新初始化。通過(guò)求解方程(7)直至穩(wěn)態(tài)可以完成重新初始化。

(7)

方程(7)中sign(φ0)被定義為

(8)

式中τ為虛時(shí)間；φ0為求解水平集方程(3)一個(gè)步長(zhǎng)以后得到的水平集函數(shù)。

1.3 擴(kuò)散方程

氣體區(qū)域和氣-固界面分別采用擴(kuò)散方程(9)和(10)描述氧化過(guò)程中反應(yīng)氣體的質(zhì)量傳輸：

(9)

-DCn=-kf/mC

(10)

式中D為反應(yīng)氣體的擴(kuò)散系數(shù)；kf/m為由材料組分決定的非均相反應(yīng)速率常數(shù)。

根據(jù)邊界條件(7),可得

v?DCnvmol

(11)

將式(11)等號(hào)右側(cè)的r向z向分量分別代入方程(4)外推，可得到外推后的速度場(chǎng)vext。

氧化過(guò)程中，材料表面不斷移動(dòng)，這意味著每計(jì)算一個(gè)步長(zhǎng)，都要重新設(shè)定邊界條件。為了避免這種情況，本文選擇文獻(xiàn)[13]中的近似方法。對(duì)邊界條件(10)作如下處理。

氧化過(guò)程中，反應(yīng)氣體在材料表面處的消耗SC可表示為

(12)

函數(shù)δ(φ)表示為

(13)

其中，ε為正比于網(wǎng)格大小的參數(shù)。函數(shù)δ(φ)的特點(diǎn)為在區(qū)間[0,2ε]內(nèi)的積分1,且在φ=ε處取最大值。

經(jīng)過(guò)處理后的邊界條件(12)可整合為方程中的源項(xiàng)：

(14)

擴(kuò)散系數(shù)被定義為φ的函數(shù)：

(15)

氣體中擴(kuò)散系數(shù)為D，為了避免材料表面附件的數(shù)值震蕩，材料一側(cè)氣體擴(kuò)散系數(shù)被定義為一個(gè)非常小的數(shù)εD。

1.4 計(jì)算方法

在COMSOL中建立如圖1(a)所示的2D軸對(duì)稱模型。如圖1(b)所示，模型中纖維半徑為r，基體厚度為t，材料初始長(zhǎng)度為h1，氣源距材料表面的初始距離為h2。模型通過(guò)映射方式劃分正方形網(wǎng)格，網(wǎng)格尺寸為10-5mm，氣源處邊界條件為給定的初始濃度。

求解C/C復(fù)合材料燒蝕形貌的具體計(jì)算流程圖如圖2所示。模擬中通過(guò)Comsol中的系數(shù)型PDE模塊寫入擴(kuò)散方程(14)，定義相應(yīng)的初始值和邊界條件。其中，C和φ均通過(guò)COMSOL with MATLAB以插值函數(shù)形式分別導(dǎo)出和導(dǎo)入COMSOL。

(a)模型基本框架 (b)局部網(wǎng)格

圖2 C/C復(fù)合材料燒蝕過(guò)程計(jì)算流程Fig.2 Flow digram of computation for C/C composites during ablation

不同時(shí)刻的相對(duì)質(zhì)量Δm的計(jì)算通過(guò)COMSOL后處理中的體積積分功能實(shí)現(xiàn)，具體過(guò)程如下：

(1)在COMSOL后處理模式中計(jì)算各時(shí)間步長(zhǎng)下的纖維體積Vf。計(jì)算公式為

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

(16)

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

(17)

式中r為2D軸對(duì)稱坐標(biāo)下的橫坐標(biāo)變量；rf為纖維半徑值；H(φ)為定義的Heaviside函數(shù)。

通過(guò)引入Heaviside函數(shù)可以很好地平滑界面處的不連續(xù)，得到更準(zhǔn)確的結(jié)果。

(2)記初始時(shí)刻的纖維體積為Vf0，其余時(shí)刻的纖維體積為Vfi。則不同時(shí)刻纖維的相對(duì)質(zhì)量Δmi=Vfi/Vf0，相應(yīng)的基體的相對(duì)質(zhì)量Δmi=Vmi/Vm0。

2 燒蝕過(guò)程相對(duì)質(zhì)量損失速率分析

2.1 穩(wěn)態(tài)結(jié)果驗(yàn)證

本文建立了半徑為3.5 μm的單根纖維及包裹在纖維外圍的厚度為0.5 μm的界面層組成的C/C復(fù)合材料模型。燒蝕的溫度大約為3000 K，設(shè)定此時(shí)的氣體擴(kuò)散系數(shù)為7×10-4m2/s[14]；此時(shí)界面層的氧化速率常數(shù)km=182 m/s[15]。圖3(a)為D=7×10-4m2/s，kf=20 m/s，km=180 m/s時(shí)計(jì)算得到的穩(wěn)態(tài)形貌，纖維形狀呈不規(guī)則的圓錐形，纖維母線為曲線；圖3(b)是單向C/C復(fù)合材料氧乙炔焰后的表面形貌，燒蝕溫度2500 ℃，O2和C2H2氣體流量分別為0.24 L/s和0.10～0.20 L/s，燒蝕時(shí)間為60 s，火焰距試樣表面距離為10 mm,焰流密度為2.38 MW/m2，與圖3(a)中的仿真形貌吻合良好。

2.2 燒蝕過(guò)程相對(duì)質(zhì)量損失速率

在氧化過(guò)程中，材料形狀發(fā)生了劇烈變化，增大了纖維的露出表面，增加了基體距離氣源的相對(duì)距離。為了準(zhǔn)確分析這種變化對(duì)纖維和基體質(zhì)量損失速率的影響，本文對(duì)纖維和基體微觀形貌達(dá)到穩(wěn)態(tài)之前的質(zhì)量損失過(guò)程進(jìn)行了分析。

Jean Lchaud提出了C/C復(fù)合材料微尺度的穩(wěn)態(tài)燒蝕模型[16]，模型定義了Sherwood數(shù)Sh和反應(yīng)速率比A，本文以這兩個(gè)參數(shù)作為變量，計(jì)算了不同反應(yīng)/擴(kuò)散機(jī)制下C/C復(fù)合材料的微觀形貌演變過(guò)程。由于材料的形貌演變過(guò)程十分緩慢，只保存形貌演變過(guò)程中的特定時(shí)間的數(shù)據(jù)。特定時(shí)間的選擇基于Vignoles G L[17]提出的參考時(shí)間公式:

τ0=rf/(C0υskf)

(18)

圖4為A=5，Sh分別取0.01、0.1、1時(shí)的纖維和基體相對(duì)質(zhì)量損失曲線；圖5為Sh=0.1，A分別取2、4、8時(shí)的纖維和基體相對(duì)質(zhì)量損失曲線。其中，Δm為纖維當(dāng)前的相對(duì)質(zhì)量。圖6為A=5時(shí)反應(yīng)氣體的濃度變化云圖。

(a)計(jì)算結(jié)果 (b)實(shí)驗(yàn)結(jié)果

(a)A=5，Sh=0.01時(shí),纖維 (b)A=5，Sh=0.1時(shí),纖維 (c)A=5，Sh=1時(shí),纖維

(d)A=5，Sh=0.01時(shí),基體 (e)A=5，Sh=0.1時(shí),基體 (f)A=5，Sh=1時(shí),基體

由圖4(a)可見，當(dāng)Sh=0.01，A=5時(shí)，在C/C復(fù)合材料微觀形貌演變過(guò)程中，纖維質(zhì)量損失速率迅速增大，在t=τ0時(shí)達(dá)到最大值。

由圖6(a)可以發(fā)現(xiàn)，在整個(gè)形貌演變過(guò)程中，材料表面反應(yīng)氣體濃度值基本一致，大致保持0.85C0的高濃度，這是因?yàn)榉磻?yīng)氣體的擴(kuò)散系數(shù)較高。但觀察形貌演變過(guò)程中5條濃度等值線(0.85C0～0.9C0)的變化，可以發(fā)現(xiàn)，最初濃度等值線完全緊貼材料表面，與纖維的輪廓相仿，露出纖維表面的濃度一致，但隨著材料形貌演變的進(jìn)行，濃度等值線不再保持之前的形狀，露出纖維底部的濃度低于纖維頂部的濃度，基體表面的反應(yīng)氣體濃度降低。因此，纖維質(zhì)量損失速率并沒(méi)有隨著纖維露出表面的增大而持續(xù)增大，在之后受擴(kuò)散的制約，纖維露出表面的增大不能加劇纖維的質(zhì)量損失。

圖4(b)、(e)為當(dāng)A=5，Sh=0.1時(shí)計(jì)算得到的燒蝕過(guò)程，τ02為參考時(shí)間?？梢园l(fā)現(xiàn)與圖4(a)類似的現(xiàn)象，在形貌演變初期(t=0～τ02),基體被大量消耗，使纖維露出表面增大，纖維質(zhì)量損失速率在燒蝕初期就達(dá)到了最大值，但質(zhì)量損失速率的增幅已不太明顯。隨著燒蝕的進(jìn)行，由于反應(yīng)氣體擴(kuò)散系數(shù)較低，纖維底部與纖維頂部出現(xiàn)了明顯的濃度差(圖6(b))，但露出纖維表面的增大還能夠略微增大纖維的質(zhì)量損失速率。但對(duì)基體來(lái)說(shuō)，其質(zhì)量損失速率迅速減小，保持穩(wěn)定(圖4(f))。

(a)Sh=0.1,A=2時(shí),纖維 (b)Sh=0.1,A=4時(shí),纖維 (c)Sh=0.1,A=8時(shí),纖維

(d)Sh=0.1,A=2時(shí),基體 (e)Sh=0.1,A=4時(shí),基體 (f)Sh=0.1,A=8時(shí),基體

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

由圖4(c)可見，當(dāng)Sh=1，A=5時(shí)，纖維質(zhì)量損失速率幾乎保持恒定，說(shuō)明燒蝕過(guò)程為純擴(kuò)散機(jī)制控制。由圖4(f)可見，當(dāng)Sh=1，A=5時(shí)，基體的質(zhì)量損失速率迅速減小，并很快達(dá)到穩(wěn)態(tài)。

由圖6可知，當(dāng)反應(yīng)速率比一定時(shí)，隨著Sh的增大，擴(kuò)散反應(yīng)堆氧化反應(yīng)的抑制逐漸增大，當(dāng)Sh=0.01時(shí)，燒蝕主要受氧化反應(yīng)控制；當(dāng)Sh=1時(shí)，燒蝕過(guò)程則主要受擴(kuò)散反應(yīng)控制。

由公式(18)可知，τ06=2τ05=4τ04，則隨著反應(yīng)速率比A的增大，基體質(zhì)量損失基體相對(duì)于纖維的相對(duì)質(zhì)量損失的變化進(jìn)一步增大，且相對(duì)質(zhì)量損失速率達(dá)到穩(wěn)態(tài)所需的時(shí)間增加，如圖5所示。

3 結(jié)論

(1)基于水平集算法，耦合擴(kuò)散方程建立了C/C復(fù)合材料在纖維尺度下的2D軸對(duì)稱燒蝕模型，計(jì)算了C/C復(fù)合材料的微觀形貌演變過(guò)程。

(2)燒蝕至穩(wěn)態(tài)的模擬結(jié)果與C/C復(fù)合材料燒蝕后的微觀形貌吻合較好，驗(yàn)證了本文模型的可靠性。

(3)分析了Sh數(shù)和反應(yīng)速率比A對(duì)纖維和基體的相對(duì)質(zhì)量損失速率的影響規(guī)律。隨著Sh的降低時(shí)纖維的相對(duì)質(zhì)量損失速率會(huì)迅速增大并在燒蝕形貌達(dá)到穩(wěn)態(tài)前穩(wěn)定，基體的相對(duì)質(zhì)量損失速率會(huì)迅速降低并在燒蝕形貌達(dá)到穩(wěn)態(tài)前穩(wěn)定；纖維和基體的相對(duì)質(zhì)量損失速率達(dá)到穩(wěn)態(tài)所需時(shí)間隨反應(yīng)速率比A的增大而增大。