近場動力學(xué)方法在C/C復(fù)合材料燒蝕問題中的應(yīng)用

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

(1.武漢理工大學(xué)力學(xué)系，湖北 武漢 430070;2.武漢理工大學(xué)理學(xué)院新材料力學(xué)理論與應(yīng)用湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430070)

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

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

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

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

1 C/C復(fù)合材料的熱化學(xué)燒蝕計算

1. 1 表面化學(xué)反應(yīng)

C/C復(fù)合材料的熱化學(xué)燒蝕主要包括碳的氧化反應(yīng)、升華、碳氮反應(yīng)以及空氣中氮?dú)夂脱鯕獾姆纸夥磻?yīng)。

熱力學(xué)分析表明，C/C復(fù)合材料在發(fā)生熱化學(xué)燒蝕時的主要化學(xué)反應(yīng)包括：

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

(1)

O2→2[O]

(2)

N2→2[N]

(3)

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

(4)

(5)

(6)

上式中：C2N、CN、CO為相應(yīng)反應(yīng)生成物。

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

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

(7)

(8)

(9)

(10)

(11)

(12)

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

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

(13)

1. 2 相容性條件

在熱化學(xué)燒蝕條件下，對于具有i種組元的反應(yīng)系統(tǒng)，表面元素組元的質(zhì)量守恒方程的形式為

(14)

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

(1+B)ciw=cie

(15)

式中：e代表反應(yīng)邊界；ci為對應(yīng)組元的質(zhì)量濃度；B為無因次質(zhì)量損失率，其表達(dá)式為

(16)

1. 3 表面能量平衡

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

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

由圖1可以寫出C/C復(fù)合材料表面熱化學(xué)燒蝕的能量守恒關(guān)系式，即凈熱流表達(dá)式:

(17)

1. 4 表面燒蝕率

C/C復(fù)合材料的燒蝕特性是通過燒蝕率來表示的，當(dāng)C/C復(fù)合材料表面在進(jìn)行熱化學(xué)燒蝕時，會伴隨著機(jī)械剝蝕的發(fā)生，根據(jù)試驗(yàn)測定的經(jīng)驗(yàn)關(guān)系式為

(18)

(19)

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

(20)

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

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

2.1 基于第一近鄰非常規(guī)態(tài)的PD溫度場的計算

2014年，Oterkus等[14]推導(dǎo)出一種拉格朗日形式的基于常規(guī)態(tài)的近場動力學(xué)瞬態(tài)熱傳導(dǎo)方程：

(21)

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

若假設(shè)當(dāng)物質(zhì)點(diǎn)鄰域半徑較小時，物質(zhì)點(diǎn)x鄰域內(nèi)各處的溫度梯度G近似相等，同時認(rèn)為只有最近鄰的物質(zhì)點(diǎn)間存在熱量的傳遞[18]，如圖2所示，Δ為粒子間距，藍(lán)色區(qū)域內(nèi)粒子僅與紅色區(qū)域內(nèi)粒子發(fā)生熱傳遞。

圖2 近場動力學(xué)瞬態(tài)熱傳導(dǎo)理論的質(zhì)點(diǎn)鄰域Fig.2 Particle neighborhood of peridynamics transient heat conduction theory

依據(jù)經(jīng)典的熱傳導(dǎo)平衡方程,可以得到第一近鄰非常規(guī)態(tài)的近場動力學(xué)瞬態(tài)熱傳導(dǎo)平衡方程如下：

(22)

式中：hs為物質(zhì)點(diǎn)x的單位時間熱生成量。

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

(23)

如果邊界Rt被施加了熱流qb，邊界上質(zhì)點(diǎn)的熱流密度為

qi=qb(i∈Ri)

(24)

由此,即可計算模型的瞬態(tài)溫度場。

2. 2 移動邊界的實(shí)現(xiàn)

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

(25)

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

若物質(zhì)點(diǎn)在燒蝕區(qū)內(nèi)(即粒子yi坐標(biāo)值大于Hi)，則時間相關(guān)函數(shù)為

(26)

則當(dāng)μ=0時，粒子脫落，燒蝕發(fā)生，實(shí)現(xiàn)邊界移動；當(dāng)μ=1時，粒子不脫落，不發(fā)生燒蝕。

3 C/C復(fù)合材料熱化學(xué)燒蝕與溫度場的耦合分析

C/C復(fù)合材料的燒蝕涉及到熱化學(xué)燒蝕計算和溫度場的計算，而傳熱計算為熱化學(xué)燒蝕計算提供了反應(yīng)溫度Tw，即

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

(27)

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

qy=qN

(28)

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

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

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

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

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

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

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

4 算例計算與分析

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

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

導(dǎo)熱率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方向的導(dǎo)熱率[W/(m2·K)];T為溫度(K)。

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

C=997+2.930×10-2T

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

q=qN

(29)

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

其他邊界設(shè)定為自由邊界條件,即

q=0

(30)

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

表1 算例條件及數(shù)值模擬計算結(jié)果與試驗(yàn)結(jié)果的對比

圖4 計算模型Fig.4 Calculation model

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

圖5 離散后的數(shù)值模型Fig.5 Discreted numerical model

4. 2 溫度場分布結(jié)果與分析

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

由圖6可見，PD方法可以實(shí)現(xiàn)對燒蝕邊界移動的自然描述。

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

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

由圖7可見，隨著燒蝕時間的增加，溫度逐漸趨于線性分布，C/C復(fù)合材料的燒蝕率隨著反應(yīng)的進(jìn)行也趨于一個常數(shù)。

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

4. 3 熱化學(xué)燒蝕計算結(jié)果與分析

圖8為不同時刻進(jìn)入C/C復(fù)合材料模型內(nèi)部的凈熱流計算結(jié)果。

圖8 不同時刻進(jìn)入C/C復(fù)合材料模型內(nèi)部的凈熱流Fig.8 Net heat flow into the C/C composites model at different times

由圖8可見，隨著燒蝕時間的增加，C/C復(fù)合材料的表面溫度不斷增加，材料與空氣的反應(yīng)更加充分，帶走更多的熱量，導(dǎo)致進(jìn)入C/C復(fù)合材料內(nèi)部的凈熱流逐漸降低，其與文獻(xiàn)[19]中的試驗(yàn)結(jié)果進(jìn)行對比，兩者結(jié)果相符。

圖9為不同時刻C/C復(fù)合材料的壁面溫度計算結(jié)果。

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

由圖9可見，燒蝕開始時，C/C復(fù)合材料表面溫度迅速上升，隨著化學(xué)反應(yīng)的進(jìn)行，帶走部分熱量，導(dǎo)致溫度上升速率逐漸變??；后期隨著燒蝕時間的增加，化學(xué)反應(yīng)、壁面溫度以及進(jìn)入C/C復(fù)合材料內(nèi)部的凈熱流達(dá)到動態(tài)平衡，材料表面溫度會逐漸趨于穩(wěn)定在2 431 K左右，與表1中試驗(yàn)結(jié)果進(jìn)行對比，兩者結(jié)果是相符的。

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

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

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

5 結(jié) 論

通過以上算例計算結(jié)果與試驗(yàn)結(jié)果的對比，可得到如下結(jié)論：

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

(2) C/C復(fù)合材料的燒蝕與壁面溫度、各組元熱焓、分壓、加載熱流及材料性質(zhì)有關(guān)，是多種因素共同作用、相互影響的結(jié)果。隨著燒蝕時間的增加，各種作用達(dá)到一個動態(tài)平衡的狀態(tài)，壁面溫度與壁面各組元熱化學(xué)反應(yīng)速率趨于穩(wěn)定， 從而導(dǎo)致C/C復(fù)合材料的燒蝕率趨于一個恒定值。