一種基于Boltzmann方程的碳基復(fù)合材料氧化燒蝕微觀界面的追蹤方法

袁經(jīng)超， 朱子鏑， 張方舟， 張 丹， 張偉剛， 李愛軍

(1. 上海大學(xué) 材料科學(xué)與工程學(xué)院，上海 200072；2. 上海大學(xué) 機(jī)電工程與自動(dòng)化學(xué)院，上海 200072； 3. 中國科學(xué)院過程工程研究所 多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100190)

1 前言

炭纖維增強(qiáng)碳基體(C/C)復(fù)合材料的超高溫?zé)g機(jī)理受到化學(xué)動(dòng)力學(xué)和對流-擴(kuò)散機(jī)制的雙重控制，在燒蝕過程中會(huì)依次出現(xiàn)氧化和升華[1-4]。但碳基復(fù)合材料的纖維相和基體相具有不同的抗氧化與升華能力，因此在高溫環(huán)境下材料的表面會(huì)產(chǎn)生不均勻的退讓，引起表面粗糙度的變化。進(jìn)一步，在熱氣流的沖擊下，露在外面的纖維受剪切力和旋渦阻力的作用會(huì)產(chǎn)生剝落，造成不可預(yù)期的表層型面變化。降低了結(jié)構(gòu)的穩(wěn)定性和可靠性，增加了結(jié)構(gòu)失效的風(fēng)險(xiǎn)[5,6]。為了對機(jī)械剝蝕更好地進(jìn)行預(yù)測，細(xì)觀尺度上化學(xué)燒蝕形貌演化過程的建模和計(jì)算尤為重要。

碳基復(fù)合材料燒蝕過程的實(shí)質(zhì)是一個(gè)多相界面的演變過程，影響因素眾多。當(dāng)前已進(jìn)行了大量的實(shí)驗(yàn)研究，如電弧噴射試驗(yàn)[7]、熱重法[8,9]和氧乙炔燒蝕試驗(yàn)[10-12]等。不少研究者在試驗(yàn)觀察碳基復(fù)合材料的燒蝕形貌的基礎(chǔ)上，采用了不同尺度的計(jì)算方法，如有限元[13]、VOF[14-16]和Level-Set[17,18]方法對燒蝕過程中界面的演變進(jìn)行了數(shù)值模擬，為理解碳基復(fù)合材料的燒蝕行為提供了有效的途徑。其中，Vignoles等[19]提出的基于布朗運(yùn)動(dòng)的隨機(jī)行走模型，成功預(yù)測了碳基復(fù)合材料單胞燒蝕的形貌變化。在這個(gè)模型中，假設(shè)氣相組分在氣固界面處被迅速消耗，濃度呈線性變化；高溫氧化反應(yīng)被簡化為粘附事件概率，并與反應(yīng)速率常數(shù)直接相關(guān)聯(lián)。模型的不足之處在于，粘附概率的選取對于該模型的有效性起到?jīng)Q定性作用，然而其物理意義卻并不明確，并且主要還依賴于經(jīng)驗(yàn)系數(shù)，必須經(jīng)過反復(fù)的優(yōu)化方能確定。同時(shí)布朗運(yùn)動(dòng)模型為概率模型，需對單個(gè)粒子運(yùn)動(dòng)進(jìn)行模擬，所需計(jì)算資源較大，且每次計(jì)算均引入隨機(jī)誤差，計(jì)算效率低，當(dāng)應(yīng)用兩相以上情況時(shí)模型復(fù)雜程度會(huì)急劇提高。

筆者采用Boltzmann方法模擬含纖維/基體界面相的碳基復(fù)合材料燒蝕界面的演化過程。作為一種介觀數(shù)值模擬方法，Boltzmann同時(shí)具有微觀與宏觀方法的計(jì)算優(yōu)勢，易應(yīng)用于復(fù)雜的幾何域[20]，易追蹤多相界面[21]，并且具有良好的并行性[22]。同時(shí)參考相場模型處理多相問題的方法，引入空間占有率概念，在Euler網(wǎng)格下，針對三相單胞模型實(shí)現(xiàn)了氣相組分濃度的預(yù)測和燒蝕界面的自捕捉，彌補(bǔ)了前人模型的不足。

2 單胞燒蝕模型

當(dāng)前，單胞燒蝕模型是被眾多研究者所認(rèn)可的[17,19,23,24]，可以反映碳基復(fù)合材料介觀尺度上燒蝕行為的一個(gè)理想簡化模型。如圖1(a)所示，一根圓柱狀炭纖維被埋在了不同材料構(gòu)成的基體里面，而炭纖維相與基體相之間有一層非常薄的中間層。

圖 1 (a) 纖維/界面相/基體的單胞模型和(b) 燒蝕過程中單胞模型變化示意圖

在初始時(shí)刻，C/C復(fù)合材料的表面是平滑的，炭纖維相與基體相完美貼合，高度一致。燒蝕過程中，高溫氧氣氣流在材料表面高速流過，并形成一定厚度的邊界層，如圖1(b)所示。在邊界層內(nèi)部氧氣幾乎是靜止的，邊界層內(nèi)部的主要傳質(zhì)是由濃度差引起的純擴(kuò)散過程[8]，可由Fick擴(kuò)散定律進(jìn)行描述：

(1)

式中,C表示氣體中擴(kuò)散物質(zhì)的濃度；D表示擴(kuò)散物質(zhì)的擴(kuò)散系數(shù)。在邊界層外部，對流傳質(zhì)將占主導(dǎo)地位，相較于邊界層內(nèi)部，邊界層外部的氧氣濃度可以看作是不變的，可采用Dirichlet邊界條件[19]：

C=C0

(2)

整體材料可以認(rèn)為是單胞模型的重復(fù)，所以橫向邊界采取周期性邊界條件。

溫度在898 K左右時(shí)，氧化反應(yīng)為體系中的主要化學(xué)反應(yīng)[8]，這時(shí)碳基復(fù)合材料的燒蝕是以碳的氧化反應(yīng)為主要反應(yīng)?？珊喕癁槿缦禄瘜W(xué)反應(yīng)方程式：

Cs+O2(g)CO2(g)

(3)

假設(shè)這是一個(gè)一級動(dòng)力學(xué)反應(yīng)[5]，氧氣在材料表面與材料發(fā)生反應(yīng)，可以認(rèn)為氧氣和材料一起從這個(gè)體系中移除，材料表面處根據(jù)反應(yīng)速率和局部反應(yīng)物濃度得到流出通量：

DC·n=k·C

(4)

其中,n為表面單位法矢量；k為表面反應(yīng)速率常數(shù)。隨著反應(yīng)的進(jìn)行，材料表面會(huì)發(fā)生退讓。盡管控制方程的形式很簡單，但對于非均勻材料，表面的這種退讓會(huì)運(yùn)動(dòng)更加復(fù)雜，由于不同相具有不同的反應(yīng)速率常數(shù)，使得燒蝕表面會(huì)出現(xiàn)類似波動(dòng)的變化，這極大地提升了求解難度。尤其對于兩相以上的復(fù)雜情況，難以得到解析。

3 數(shù)值模擬方法

格子Boltzmann方法在求解復(fù)雜流場問題上有很好的優(yōu)勢，且可方便地拓展多路徑的化學(xué)反應(yīng)過程?；镜腂GK格子Boltzmann方程如下[22]：

+ΔtwiS

(5)

通過Chapman-Enskog展開可以建立宏觀的擴(kuò)散系數(shù)與微觀的碰撞頻率之間的關(guān)系[22]：

(6)

式中,Δx表示極小的空間步長，各個(gè)方向的空間步長相同,d表示維度。

對于二維軸對稱問題，需額外添加一項(xiàng)以應(yīng)對徑向空間的變化[25]。

(7)

(8)

式中r表示圓柱的徑向,j表示徑向的離數(shù)，f1和f3方向由圖1(b)給出。式(8)較(7)更為精確[25]，式(7) 或者(8)可以直接添加進(jìn)源項(xiàng)。

盡管帶有源項(xiàng)的BGK格子Boltzmann方程可以很好地描述含有均相體反應(yīng)的擴(kuò)散-反應(yīng)體系，但碳基復(fù)合材料的不均勻性，使得數(shù)值算法必須面對非均相表面反應(yīng)體系。為此，參考相場模型的處理方法，構(gòu)建一個(gè)新的變量Φi(0≤Φi<1)用于區(qū)分不同的固體與氣體，0表示不存在固相i，1表示完全為固相i。Φi代表的是某一種固相i在空間的占有率。引入Φi之后，氧氣與固相i反應(yīng)的消耗量Ri，即式(4)中的流出通量，可以通過含Φi的式子表示：

Ri=ΦikiC

(9)

固相i對應(yīng)的燒蝕量可以通過氧氣的消耗量Ri得到，因此Φi可以表示為：

(10)

式中,Mi表示固相i的分子質(zhì)量；ρi表示固相i的密度。Φi值的求解可確定材料內(nèi)不同固體相所占比例。

4 結(jié)果討論

4.1 結(jié)果驗(yàn)證

Vignoles等[5]在一定假設(shè)下對只包含基體和纖維相的單胞模型求得解析解。基體相較纖維擁有更高的反應(yīng)活性，即更大的反應(yīng)速率常數(shù)[8]。解析解表明：由于基體的反應(yīng)速率大于纖維，基體表面率先發(fā)生衰退，使得纖維逐漸露出，進(jìn)而纖維側(cè)向開始被氧化，使得纖維逐漸變細(xì)，最終變?yōu)獒樞?；兩個(gè)重要的參數(shù)可反映出穩(wěn)態(tài)下單胞的表面粗糙度，即基體與纖維的反應(yīng)活性對比A以及舍伍德數(shù)Sh。

(11)

(12)

圖 2 兩相單胞模型的數(shù)值計(jì)算結(jié)果與解析解的對比(A=5)

圖 3 兩相單胞模型濃度場和質(zhì)量傳輸流線 (A=5)

Vignoles等[5]在計(jì)算解析解時(shí)假設(shè)氧氣濃度延垂直方向呈線性變化，而水平方向無變化。本文的模型無需此假設(shè)，可同時(shí)計(jì)算出氣相的濃度場。圖3為燒蝕達(dá)到穩(wěn)態(tài)時(shí)氧氣濃度的等值線，顯然，在遠(yuǎn)離燒蝕界面處水平方向上不存在濃度梯度，假設(shè)成立，橫向的濃度傳遞可以忽略。而在燒蝕界面附近，等值線沿著表面發(fā)生曲折，表明水平方向上也存在有濃度差。圖3進(jìn)一步表明，擴(kuò)散限制模式下(Sh遠(yuǎn)大于1)，固體表面極大的反應(yīng)速率，使得固體表面濃度幾乎為0，氧氣從邊界層頂部源源不斷向表面擴(kuò)散；反應(yīng)限制模式(Sh遠(yuǎn)小于1)，由于反應(yīng)速率遠(yuǎn)小于擴(kuò)散速率，所以整個(gè)邊界層內(nèi)部濃度均達(dá)到邊界層外部濃度C0。僅在反應(yīng)限制模式下，不需要計(jì)算濃度場，可將材料非均相的性質(zhì)解耦出來計(jì)算。

4.2 纖維和基體界面相

碳基復(fù)合材料在炭纖維與基體之間往往存在一個(gè)界面相，是增強(qiáng)相與基體相連接的橋梁，也是應(yīng)力的傳遞者。通常界面相較于基體有更為活潑的化學(xué)性質(zhì)，即擁有更大的反應(yīng)速率常數(shù)[5]。因此包含中間相的單胞模型更接近C/C復(fù)合材料的真實(shí)狀態(tài)，同時(shí)也可以進(jìn)一步檢驗(yàn)算法在更為復(fù)雜的多相材料燒蝕問題求解的穩(wěn)定性。本文在基體和纖維不變的情況下，加入了中間相，假設(shè)反應(yīng)速率常數(shù)ki=8km。

圖 4 不同舍伍德數(shù)下三相單胞燒蝕模型的界面燒蝕形貌

圖4為不同舍伍德數(shù)下，包含纖維-界面相-基體的單胞燒蝕模型的計(jì)算結(jié)果。如圖4所示，在擴(kuò)散限制模式下，材料的非均相性無法體現(xiàn)，燒蝕形貌與圖2基本一致；而在反應(yīng)限制模式下，材料的不均勻性明顯體現(xiàn)，纖維露出高度，不論相對于基體還是界面相，都明顯提升。由于界面相的反應(yīng)速率常數(shù)大于基體的反應(yīng)速率常數(shù)，所以體系中界面衰退最快的是界面相。界面相將基體與纖維分成了兩部分，界面相與基體構(gòu)成一個(gè)體系，與纖維構(gòu)成另一個(gè)體系，分別與最初單胞模型一致。因此，對于每個(gè)體系主要的控制參數(shù)依然為舍伍德數(shù)Sh與反應(yīng)活性比A。若將中間相的參數(shù)作為Sh和A的分子項(xiàng)，纖維—界面相系統(tǒng)參數(shù)可寫為：

(13)

(14)

基體—界面相系統(tǒng)參數(shù)可寫為：

(15)

(16)

圖5中左圖和右圖分別是將包含基體-界面相-纖維的三相體系拆分成界面相-纖維和界面相-基體進(jìn)行計(jì)算所得的結(jié)果。左圖和右圖相互疊加可以獲得與中間圖及其相似的結(jié)果。這表明對單胞燒蝕模型而言，多相材料的粗糙度取決于相鄰兩相之間的狀態(tài)。相鄰兩相之間的作用不會(huì)受到其它相作用的影響。并且與相鄰兩相的厚度無關(guān)。但同時(shí)也表明：中間相作為化學(xué)反應(yīng)最活潑的相，會(huì)增加纖維的露出高度，增大材料整體的表面粗糙度。

圖 5 三相模型與兩相模型對比

4.3 燒蝕深度

圖 6 反應(yīng)活性比(A2)與燒蝕深度(Lp)的關(guān)系(A=5)

燒蝕極限Lpm深度在Rf確定時(shí)僅與Sh相關(guān)，對此對Lpm進(jìn)行無量綱轉(zhuǎn)換，并且做出與Sh相關(guān)的圖7。圖7表明了Sh與Rf/Lpm的線性關(guān)系，黑色點(diǎn)為圖6中各個(gè)Sh下的極限燒蝕深度Lpm，實(shí)線為線性的擬合函數(shù)，虛線是根據(jù)現(xiàn)有點(diǎn)數(shù)據(jù)所構(gòu)成的誤差邊界，菱形點(diǎn)為隨機(jī)選取Sh模擬結(jié)果，用來驗(yàn)證預(yù)測模型的可靠性。

圖 7 舍伍德數(shù)(Sh)與極限燒蝕深度(Lpm)之間的關(guān)系(A=5)

顯然驗(yàn)證點(diǎn)很好的落在了預(yù)測邊界以內(nèi)，這表明預(yù)測的直線具有很好的準(zhǔn)確性。此外，這條直線近乎于通過(0,0)點(diǎn)，Sh趨向于0的時(shí)候，即基體和纖維的反應(yīng)系數(shù)趨向于0，理論上燒蝕深度應(yīng)該趨向于，然而這里存在穩(wěn)態(tài)的概念，一般燒蝕輪廓不發(fā)生變化時(shí)即為穩(wěn)態(tài)，而(0,0)處不存在穩(wěn)態(tài)的情況，所以Rf/Lpm受到擴(kuò)散與達(dá)到穩(wěn)態(tài)的時(shí)間影響，具有一定深度。另一方面，直線的斜率代表了氣相擴(kuò)散進(jìn)入燒蝕孔內(nèi)的傳質(zhì)阻力，氣相一邊與基體與纖維反應(yīng)，一邊擴(kuò)散進(jìn)入燒蝕孔內(nèi)部，達(dá)到穩(wěn)態(tài)時(shí)，擴(kuò)散最深的狀態(tài)便是燒蝕的最大深度。若假設(shè)截距為0，斜率可由任意Sh下的一組數(shù)據(jù)估計(jì)，由此可以估算出其他Sh下的Lpm，即可預(yù)測任意體系下的燒蝕深度或者露出高度?？蓴M合出如下關(guān)系式：

(17)

5 結(jié)論

基于Boltzmann方法，通過構(gòu)建空間占有率參數(shù)Φi，提出了擴(kuò)散—表面非均相反應(yīng)體系的新算法。該算法應(yīng)用于纖維—基體兩相單胞模型時(shí)，計(jì)算結(jié)果與解析解具有很好的一致性，驗(yàn)證了算法的可靠性。并且可以獲得以往數(shù)值算法無法得到的濃度場的信息。氣相濃度的線性化假設(shè)在遠(yuǎn)離界面處具有合理性，但在界面附近明顯呈非線性分布。進(jìn)一步，將算法用于求解纖維-基體-界面相的三相單胞模型。氧化燒蝕的形貌僅取決于相鄰兩相之間的狀態(tài)。含界面相的氧化燒蝕，存在最大燒蝕深度。并且通過數(shù)學(xué)變換可以得到燒蝕深度與各參數(shù)之間的關(guān)系，提出了預(yù)測燒蝕深度的數(shù)學(xué)表達(dá)式。在求解過程中，提出的算法在多相結(jié)構(gòu)計(jì)算中表現(xiàn)出了良好的穩(wěn)定性。