近場動力學方法在C/C復合材料燒蝕問題中的應用

陳明偉,劉齊文，2,劉立勝，2

(1.武漢理工大學力學系，湖北 武漢 430070;2.武漢理工大學理學院新材料力學理論與應用湖北省重點實驗室，湖北 武漢 430070)

C/C復合材料具有抗燒蝕、抗沖擊、高比強度、高比模量等特點，在航天器和武器等飛行器的熱防護安全領域得到了廣泛的應用，目前正逐步取代傳統防熱抗燒蝕材料，成為飛行器結構設計中防熱結構的主要材料。飛行器在高溫、高焓及高熱通量等惡劣環境下，如航天器進入大氣層時，經受7 000 K左右超高溫、幾十兆瓦熱通量、100個標準重力加速度的過載，此時通過材料自身熱化學燒蝕以及機械剝蝕引起質量損失，帶走大部分的熱量，阻止外部熱量向飛行器內部傳遞，從而保護內部工作人員的安全與內部結構的正常工作[1]。鑒于C/C復合材料在飛行器安全領域中應用的重要性，目前對C/C復合材料的研究也愈加廣泛。通過對C/C復合材料熱化學燒蝕性能的研究，可增加對其內部結構及其服役性能的了解，以為飛行器熱防護材料的設計和研究提供理論支持。

C/C復合材料在高熱通量條件下，由于材料本身與高熱氣流的物理、化學和力學作用等，會發生質量損失。目前，針對C/C復合材料的燒蝕研究，已經提出了一些有效的燒蝕計算模型，并可以計算在一定燒蝕條件下材料的線燒蝕率[2]。現有主要的燒蝕模型有熱化學燒蝕模型、機械剝蝕模型等[3-4]。根據熱化學燒蝕理論可以對C/C復合材料的熱化學燒蝕行為進行建模，并對材料的質量損失與線燒蝕率的變化規律進行理論上的預測[5]。

C/C復合材料的燒蝕過程是一個包括質量損失和結構失效的復雜過程，是一個非常典型的非線性、不連續問題[6-7]。目前對此類問題的研究主要有實驗和數值模擬兩種方法，但實驗方法條件苛刻、經濟成本過高，而數值模擬方法則恰好彌補了實驗方法的缺點，又可滿足一定的精度。數值模擬方法近年來發展迅速，已被廣泛應用于生物、建筑、機械等多個領域。如徐遲等[8]和陳超等[9]利用數值模擬方法分別模擬計算了三峽庫區某滑坡庫水位與降雨聯合作用滲流應力耦合和汽車油罐車爆炸燃燒特性，獲得了較好的計算精度。此外，在處理不連續問題時，傳統的有限單元法、有限差分法等一些宏觀方法的計算結果主要依賴于網格尺寸及質量[10]。燒蝕問題包含邊界移動，在燒蝕發生后，因材料的質量損失也必須重新劃分網格，計算結果也具有強烈的網格依賴性。而近場動力學(peridynamics,PD)理論的提出與發展突破了經典理論的限制，它基于非局域思想建模，采用空間積分方程代替偏微分方程來描述物質點的運動，自然而然地得到對不連續現象的描述，從而可以避免傳統計算方法的奇異性和計算效率的問題[11-12]。PD方法經過多年的發展，已成功運用于許多研究領域。如Ha等[13]采用基于鍵型PD理論研究了動態脆性裂紋問題；Oterkus等[14]采用基于常規態的PD方法研究了熱擴散的問題；Bobaru等[15]采用基于鍵型PD方法研究了瞬態熱傳導問題;廖洋等[16]采用PD強度折減法分析了邊坡的穩定性問題；王超聰等[17]采用PD方法對復合材料的燒蝕溫度場進行了模擬。

本文在C/C復合材料熱化學燒蝕理論的基礎上，建立了C/C復合材料在熱化學燒蝕與溫度場相互耦合條件下的PD求解模型。首先，根據熱化學燒蝕理論計算進入C/C復合材料模型內部的凈熱流；然后，根據凈熱流進一步由PD傳熱理論計算移動邊界條件下C/C復合材料的燒蝕溫度場。在對燒蝕表面退縮問題求解時，滿足燒蝕判據的粒子自動脫落，進而實現了燒蝕面的移動，可以直觀地表征燒蝕移動邊界條件下瞬態溫度場的分布。

1 C/C復合材料的熱化學燒蝕計算

1. 1 表面化學反應

C/C復合材料的熱化學燒蝕主要包括碳的氧化反應、升華、碳氮反應以及空氣中氮氣和氧氣的分解反應。

熱力學分析表明，C/C復合材料在發生熱化學燒蝕時的主要化學反應包括：

iC(s)→Ci(g)(i=1,2,3)

(1)

O2→2[O]

(2)

N2→2[N]

(3)

上式中：C(s)為固態碳；Ci為升華碳對應分子：C1、C2、C3，且式(1)中包含三個反應式；O2為氧氣分子；[O]為氧離子；N2為氮氣分子；[N]為氮離子；s表示組元為固態；g表示組元為氣態。

(4)

(5)

(6)

上式中：C2N、CN、CO為相應反應生成物。

與反應方程式(1)～(6)相對應的平衡常數k分別為

kpci=pci(i=1,2,3)

(7)

(8)

(9)

(10)

(11)

(12)

上式中：pC1、pC2、pC3、pCO、pO2、pC2N、pN2、pCN、pO、pN為對應的各組元相對分壓。

平衡常數k僅僅與溫度相關，可通過查閱《物理化學手冊》得到各反應的平衡常數與溫度的關系式；公式(7)～(12)有8個方程和10個未知數，因此需要補充質量守恒方程來封閉方程組。引入質量濃度ci(單位為1)，ci與組元分壓的關系(可直接由物態方程導出)為

(13)

1. 2 相容性條件

在熱化學燒蝕條件下，對于具有i種組元的反應系統，表面元素組元的質量守恒方程的形式為

(14)

由熱化學邊界層理論，在Le=Pr=1(Le為路易斯數，Pr為普朗特數)的條件下，有以下壁面質量守恒條件，即相容性條件：

(1+B)ciw=cie

(15)

式中：e代表反應邊界；ci為對應組元的質量濃度；B為無因次質量損失率，其表達式為

(16)

1. 3 表面能量平衡

在進行熱化學燒蝕計算時，首先根據上個時刻的熱流密度條件利用PD傳熱理論得到在下一個時刻的壁面溫度， 然后計算該溫度條件下相應的各個物理化學參數，最后通過C/C復合材料表面能量平衡方程求得進入C/C復合材料壁面內部的凈熱流，至此可完成一個完整的循環，反應壁面的能量平衡見圖1。

圖1 反應壁面的能量平衡示意圖Fig.1 Schematic diagram of wall energy balance

由圖1可以寫出C/C復合材料表面熱化學燒蝕的能量守恒關系式，即凈熱流表達式:

(17)

1. 4 表面燒蝕率

C/C復合材料的燒蝕特性是通過燒蝕率來表示的，當C/C復合材料表面在進行熱化學燒蝕時，會伴隨著機械剝蝕的發生，根據試驗測定的經驗關系式為

(18)

(19)

求解出材料總的固體質量損失率后，即可由下式計算材料的表面燒蝕率v-∞(m/s)：

(20)

式中：ρ為材料的密度(kg/m3)。

2 移動邊界下模型溫度場的計算

2.1 基于第一近鄰非常規態的PD溫度場的計算

2014年，Oterkus等[14]推導出一種拉格朗日形式的基于常規態的近場動力學瞬態熱傳導方程：

(21)

式中：ρ為材料的密度(kg/m3);cH為材料的比熱容[J/(kg·K)];h為熱流態;hs(x,t)為質點x的單位時間熱生成量。

若假設當物質點鄰域半徑較小時，物質點x鄰域內各處的溫度梯度G近似相等，同時認為只有最近鄰的物質點間存在熱量的傳遞[18]，如圖2所示，Δ為粒子間距，藍色區域內粒子僅與紅色區域內粒子發生熱傳遞。

圖2 近場動力學瞬態熱傳導理論的質點鄰域Fig.2 Particle neighborhood of peridynamics transient heat conduction theory

依據經典的熱傳導平衡方程,可以得到第一近鄰非常規態的近場動力學瞬態熱傳導平衡方程如下：

(22)

式中：hs為物質點x的單位時間熱生成量。

熱流密度的散度表示如下：

(23)

如果邊界Rt被施加了熱流qb，邊界上質點的熱流密度為

qi=qb(i∈Ri)

(24)

由此,即可計算模型的瞬態溫度場。

2. 2 移動邊界的實現

采用PD理論實現燒蝕邊界的移動并不需要刪除網格、重建網格，這是由于模型是由物質點組成的，在由熱化學燒蝕計算得到燒蝕率v-∞后，可由時間積分得到模型的燒蝕深度。由每個時間步燒蝕率疊加得到模型的燒蝕深度，即可得到燒蝕安全區為

(25)

式中:H為材料的原始厚度(m)；Hi為燒蝕后材料的厚度(m)。

若物質點在燒蝕區內(即粒子yi坐標值大于Hi)，則時間相關函數為

(26)

則當μ=0時，粒子脫落，燒蝕發生，實現邊界移動；當μ=1時，粒子不脫落，不發生燒蝕。

3 C/C復合材料熱化學燒蝕與溫度場的耦合分析

C/C復合材料的燒蝕涉及到熱化學燒蝕計算和溫度場的計算，而傳熱計算為熱化學燒蝕計算提供了反應溫度Tw，即

h(x′,t)]dAx′+hs(x,t)

(27)

公式(17)關于C/C復合材料表面熱化學燒蝕計算為傳熱計算提供了熱流邊界條件qN，對于燒蝕層粒子而言，其熱流密度為

qy=qN

(28)

因此,可以對C/C復合材料的燒蝕和溫度場進行耦合求解。具體耦合求解流程如下(見圖3)：

(1) 設置初始條件，生成鄰域，利用PD理論對溫度場進行求解，得到當前時刻的壁面溫度Tw。

(2) 由壁面溫度Tw代入熱化學燒蝕計算公式，求得該時間步的表面燒蝕率v-∞和進入材料內部的凈熱流qN，并判斷燒蝕安全區Hi。

(3) 增加一個時間步Δt，判斷下一時刻t是否大于總的計算時間t總，然后決定結束計算或者進行下一個時間步的計算。

(4) 搜索粒子形成鄰域，將凈熱流qN作為邊界條件施加到更新模型的燒蝕層粒子上,并耦合PD熱傳導方程,求解當前時間步的壁面溫度Tw。

(5) 重復上述步驟(2)～(4)。

圖3 C/C復合材料熱化學燒蝕與溫度場的耦合計算流程圖Fig.3 Flow chart for coupling calculation of ablation and temperature field of C/C composites

4 算例計算與分析

4. 1 計算模型與算例條件

參考文獻[19]的計算模型，并將其簡化為二維計算模型，采用本文中的計算方法對模型溫度場、燒蝕率以及進入C/C復合材料內部的凈熱流等進行計算，并將數值模擬計算結果與文獻[19]中的試驗結果進行了對比。計算模型采用文獻[19]中試驗測得的變物性參數，具體如下：

導熱率K為

Kx=172.528-0.236T+1.994×10-5T2-7.792×10-9T3+1.14×10-11T4

Ky=82.785-0.103T+9.846×10-6T2-4.372×10-9T3+7.315×10-13T4

式中：Kx、Ky分別為x、y方向的導熱率[W/(m2·K)];T為溫度(K)。

比熱容C[J/(kg·K)]為

C=997+2.930×10-2T

根據能量守恒，給定熱流的邊界條件為

q=qN

(29)

其中,qN由熱化學燒蝕計算以及壁面能量平衡方程計算得到，當y=0時，溫度邊界為Ty=0=300 K。

其他邊界設定為自由邊界條件,即

q=0

(30)

表1為算例條件以及數值模擬計算結果與試驗結果的對比。本文在駐點壓力Ps為0.598 MPa、熱焓Hs為8 051 kJ/kg的條件下，對C/C復合材料的熱化學燒蝕性能進行分析。計算模型的尺寸為0.050 5 m(長)×0.016 m(寬)，計算y方向與熱流方向平行，見圖4。

表1 算例條件及數值模擬計算結果與試驗結果的對比

圖4 計算模型Fig.4 Calculation model

圖5為離散后的數值模型，該模型中粒子數為m(506)×n(161)=81 466，粒子間距為0.000 1 m；背景單元大小為0.000 1 m×0.000 1 m，數目為80 800，這里所說的背景單元只是粒子分布的方式，實際上并不存在也不參與計算。

圖5 離散后的數值模型Fig.5 Discreted numerical model

4. 2 溫度場分布結果與分析

圖6為不同時刻(7 s、14 s、21 s和28 s)模型燒蝕溫度場分布。

由圖6可見，PD方法可以實現對燒蝕邊界移動的自然描述。

圖6 不同時刻模型燒蝕溫度場分布Fig.6 Temperature field distribution of the ablation model at different times

與圖6相對應的溫度場分布曲線見圖7，即為不同時刻模型x=0.008時沿y方向的溫度分布。

由圖7可見，隨著燒蝕時間的增加，溫度逐漸趨于線性分布，C/C復合材料的燒蝕率隨著反應的進行也趨于一個常數。

圖7 不同時刻模型x=0.008時沿y方向的溫度分布Fig.7 Temperature distribution of the model along y-direction at different times

4. 3 熱化學燒蝕計算結果與分析

圖8為不同時刻進入C/C復合材料模型內部的凈熱流計算結果。

圖8 不同時刻進入C/C復合材料模型內部的凈熱流Fig.8 Net heat flow into the C/C composites model at different times

由圖8可見，隨著燒蝕時間的增加，C/C復合材料的表面溫度不斷增加，材料與空氣的反應更加充分，帶走更多的熱量，導致進入C/C復合材料內部的凈熱流逐漸降低，其與文獻[19]中的試驗結果進行對比，兩者結果相符。

圖9為不同時刻C/C復合材料的壁面溫度計算結果。

圖9 不同時刻C/C復合材料的壁面溫度Fig.9 Wall temperature of the C/C composites at different times

由圖9可見，燒蝕開始時，C/C復合材料表面溫度迅速上升，隨著化學反應的進行，帶走部分熱量，導致溫度上升速率逐漸變小；后期隨著燒蝕時間的增加，化學反應、壁面溫度以及進入C/C復合材料內部的凈熱流達到動態平衡，材料表面溫度會逐漸趨于穩定在2 431 K左右，與表1中試驗結果進行對比，兩者結果是相符的。

圖10為不同時刻C/C復合材料的線燒蝕率(v-∞)計算結果。

圖10 不同時刻C/C復合材料的線燒蝕率Fig.10 Line ablation rate of the C/C composites at different times

由圖10可見，熱化學燒蝕剛開始時，C/C復合材料的線壁面溫度較低，反應緩慢，燒蝕率低；隨著燒蝕時間的增加，壁面溫度上升，化學反應速度增大，反應迅速，C/C復合材料的線燒蝕率升高；最后整個過程達到動態平衡，C/C復合材料的線燒蝕率趨于穩定，且與表1中試驗結果相符。

5 結 論

通過以上算例計算結果與試驗結果的對比，可得到如下結論：

(1) 本文提出的基于熱化學燒蝕理論的PD傳熱方法，在引入損傷函數后，實現了對燒蝕邊界的自然描述，不需要引入其他復雜判定條件。利用該方法對C/C復合材料的壁面溫度、凈熱流、燒蝕率進行求解，其計算結果與試驗結果基本吻合，說明PD方法可用于C/C復合材料熱化學燒蝕與溫度場的耦合分析。

(2) C/C復合材料的燒蝕與壁面溫度、各組元熱焓、分壓、加載熱流及材料性質有關，是多種因素共同作用、相互影響的結果。隨著燒蝕時間的增加，各種作用達到一個動態平衡的狀態，壁面溫度與壁面各組元熱化學反應速率趨于穩定， 從而導致C/C復合材料的燒蝕率趨于一個恒定值。