基于多場(chǎng)耦合碳/碳復(fù)合材料傳熱及燒蝕響應(yīng)

孫學(xué)文，楊海波，米 濤

1) 北京科技大學(xué)機(jī)械工程學(xué)院，北京 100083 2) 北京科技大學(xué)流體與材料相互作用教育部重點(diǎn)實(shí)驗(yàn)室，北京 100083

3) 東莞材料基因高等理工研究院，東莞 523808 4) 北華航天工業(yè)學(xué)院材料工程學(xué)院，廊坊 065000

隨著高超聲速技術(shù)的發(fā)展，飛行器服役環(huán)境越來(lái)越惡劣，對(duì)熱防護(hù)材料提出了更高的要求[1?3].碳/碳復(fù)合材料具有較高的化學(xué)潛熱，在高溫環(huán)境下仍保持較高的強(qiáng)度，被廣泛地應(yīng)用于飛行器的熱防護(hù)系統(tǒng)[4?6]. 飛行器服役過(guò)程中周圍的熱化學(xué)非平衡流場(chǎng)會(huì)影響碳/碳復(fù)合材料的傳熱及燒蝕，導(dǎo)致防熱層的厚度和形狀發(fā)生變化，材料燒蝕產(chǎn)物又與周圍的高溫氣體發(fā)生化學(xué)反應(yīng)，從而改變流場(chǎng)的溫度及組元濃度等特性，流場(chǎng)特性的改變反過(guò)來(lái)又會(huì)影響飛行器防熱層的燒蝕[7?10]. 因此，高超聲速熱化學(xué)非平衡流場(chǎng)與碳/碳復(fù)合材料之間存在著強(qiáng)烈的耦合作用，準(zhǔn)確的預(yù)測(cè)碳/碳復(fù)合材料的內(nèi)部溫度分布以及燒蝕響應(yīng)，對(duì)熱防護(hù)系統(tǒng)的設(shè)計(jì)及優(yōu)化具有重要的意義.

對(duì)熱防護(hù)材料燒蝕性能的預(yù)測(cè)一直是發(fā)展高超聲速飛行器所面臨的關(guān)鍵問(wèn)題之一，近年來(lái)國(guó)內(nèi)外相關(guān)學(xué)者對(duì)其開展了一定的研究[11?13]. Martin與Boyd[14]建立了來(lái)流氣體與熱防護(hù)材料的耦合模型，主要分析了來(lái)流環(huán)境下熱防護(hù)材料的溫度響應(yīng)以及防熱材料燒蝕產(chǎn)物對(duì)外部流場(chǎng)的影響，并以IRV-2為計(jì)算模型對(duì)數(shù)值方法進(jìn)行了驗(yàn)證.Cross與Boyd[15]針對(duì)火箭噴管環(huán)境下防熱材料的響應(yīng)進(jìn)行研究，建立了噴管內(nèi)流場(chǎng)與碳/酚醛材料的耦合反應(yīng)模型，對(duì)噴管流場(chǎng)環(huán)境以及材料響應(yīng)進(jìn)行了研究. Mortensen與Zhong[16]針對(duì)防熱材料表面燒蝕對(duì)高超聲速邊界層的影響，考慮了材料表面的燒蝕和真實(shí)氣體效應(yīng)，建立了流場(chǎng)的熱化學(xué)非平衡模型，并以鈍椎體為例進(jìn)行驗(yàn)證，分析了表面燒蝕對(duì)高超聲速邊界層的影響. Chen等[17?18]針對(duì)碳化材料燒蝕響應(yīng)預(yù)測(cè)，開發(fā)了全隱式的燒蝕熱響應(yīng)程序，用于模擬防熱材料熱解氣體流動(dòng)、熱化學(xué)燒蝕以及外形變化等，最后采用三組算例對(duì)計(jì)算程序進(jìn)行了驗(yàn)證. Kumar[19]建立了防熱材料燒蝕與外部流場(chǎng)耦合的數(shù)值模型，重點(diǎn)分析了碳化材料熱解過(guò)程，以及熱解氣體與外部流場(chǎng)之間的相互影響關(guān)系. Li等[20]針對(duì)碳化復(fù)合材料，建立了燒蝕面后退的非線性熱解層模型，對(duì)其燒蝕響應(yīng)進(jìn)行預(yù)測(cè). Candler等[21]針對(duì)碳基復(fù)合材料的氧化燒蝕過(guò)程，比較了有限速率反應(yīng)模型以及化學(xué)平衡兩種氣?固反應(yīng)模型對(duì)材料響應(yīng)的預(yù)測(cè)，并采用球錐模型進(jìn)行驗(yàn)證. Qin等[22]針對(duì)火箭發(fā)動(dòng)機(jī)內(nèi)部的熱環(huán)境，建立碳/碳復(fù)合材料的多尺度熱化學(xué)燒蝕模型，考慮基體及纖維的反應(yīng)速率，對(duì)復(fù)合材料的燒蝕響應(yīng)進(jìn)行預(yù)測(cè). Yin 等[23]，Meng 等[24]，Chen[25]考慮來(lái)流與固體材料的單向耦合作用，對(duì)碳/碳復(fù)合材料的燒蝕響應(yīng)進(jìn)行計(jì)算，預(yù)測(cè)燒蝕體的外形變化和溫度分布等情況. 現(xiàn)有防熱材料響應(yīng)模型通常將氣動(dòng)熱載荷直接作用于材料模型表面進(jìn)行材料響應(yīng)的單向耦合預(yù)測(cè)，未同時(shí)考慮由于材料溫度及燒蝕外形變化等對(duì)外部流場(chǎng)的影響.

基于以上分析，本文在熱化學(xué)燒蝕理論的基礎(chǔ)上，采用流?熱?燒蝕多場(chǎng)耦合策略，考慮外部流場(chǎng)熱化學(xué)非平衡效應(yīng)、材料傳熱以及表面燒蝕等因素，建立了碳/碳復(fù)合材料在高超聲速環(huán)境下的雙向耦合模型，對(duì)其傳熱及燒蝕響應(yīng)進(jìn)行預(yù)測(cè)，分析不同時(shí)刻材料模型的溫度、燒蝕速率以及燒蝕外形的變化.

1 控制方程及多場(chǎng)耦合策略

碳/碳復(fù)合材料燒蝕過(guò)程的數(shù)值模擬需要考慮復(fù)雜的物理及化學(xué)過(guò)程，包括氣動(dòng)加熱、表面燒蝕、材料熱響應(yīng)以及燒蝕邊界移動(dòng)等. 本文通過(guò)Fluent對(duì)外部流場(chǎng)進(jìn)行建模計(jì)算，得到材料外部的氣動(dòng)熱載荷，采用Abaqus計(jì)算材料的瞬態(tài)熱響應(yīng)，考慮有限速率燒蝕模型用來(lái)預(yù)測(cè)碳/碳復(fù)合材料的燒蝕后退，并用網(wǎng)格移動(dòng)策略對(duì)燒蝕邊界進(jìn)行追蹤，最后利用Mpcci實(shí)現(xiàn)流場(chǎng)與固體材料之間的數(shù)據(jù)傳遞.

1.1 外部流場(chǎng)控制方程

針對(duì)多組元的化學(xué)反應(yīng)氣體混合物，可壓縮黏性熱化學(xué)非平衡流動(dòng)的Navier-Stokes控制方程組在直角坐標(biāo)系中表達(dá)形式如下：

式中，Q為守恒變量，E、F分別為x、y方向的無(wú)黏通量，Ev、Fv分別為x、y方向的黏性通量，S為反應(yīng)源項(xiàng)體現(xiàn)化學(xué)非平衡的影響.

式中，u、v分別為x、y方向的運(yùn)動(dòng)速度，ρ、p分別為氣體壓力和密度，E、H為單位質(zhì)量氣體的總能量和總焓，ρi、i、Yi分別為氣體組元i的密度、化學(xué)反應(yīng)源項(xiàng)、質(zhì)量分?jǐn)?shù)，Dim、hi為 氣體組元i的擴(kuò)散系數(shù)和單位質(zhì)量絕對(duì)焓，eve、ve表示單位質(zhì)量氣體的振動(dòng)能及振動(dòng)能源項(xiàng)，qx、qvex為氣體在x方向的平動(dòng)?轉(zhuǎn)動(dòng)熱流和振動(dòng)熱流，qy、qvey為 氣體在y方向的平動(dòng)?轉(zhuǎn)動(dòng)熱流和振動(dòng)熱流，τxx、τxy、τyy、τyx為剪切應(yīng)力張量分量，ns為氣體組元數(shù)量.

采用Gupta化學(xué)動(dòng)力學(xué)模型計(jì)算化學(xué)反應(yīng)引起的組元變化，具體化學(xué)反應(yīng)模型見參考文獻(xiàn)[26].

1.2 氣?固界面及燒蝕速率模型

來(lái)流氣體與固體材料之間需要滿足能量守恒，圖1為氣?固界面能量傳遞示意，其滿足關(guān)系式：

式中，q為 傳入材料模型內(nèi)部的熱流，qconv為外部流場(chǎng)對(duì)材料壁面的氣動(dòng)加熱熱流，qrad?in為外部流場(chǎng)對(duì)材料的輻射加熱量，chcs為壁面化學(xué)反應(yīng)產(chǎn)生的熱量，hcs為 燒蝕質(zhì)量損失帶走的熱量，qrad?out為材料對(duì)外部流場(chǎng)的輻射散熱量.

圖1 氣?固界面能量傳遞示意圖Fig.1 Energy transfer at the gas?solid interface

采用有限速率燒蝕模型來(lái)模擬碳/碳復(fù)合材料的表面燒蝕，不考慮由于機(jī)械剝蝕及材料表面熔化等物理因素引起的質(zhì)量損失，只考慮由于表面氣固化學(xué)反應(yīng)引起的質(zhì)量變化.

碳/碳復(fù)合材料表面燒蝕主要包括碳的氧化和氮化，表面化學(xué)反應(yīng)機(jī)制及質(zhì)量損失率如下[24]：

（1）碳的氧化反應(yīng)：O+Cs→CO

質(zhì)量損失速率為：

（2）碳的氧化反應(yīng)：O2+2Cs→2CO

質(zhì)量損失速率為：

（3）碳的氮化反應(yīng)：N+Cs→CN

質(zhì)量損失速率為：

材料表面總的質(zhì)量損失率為：

1.3 材料熱響應(yīng)模型及燒蝕邊界追蹤

固體材料的熱傳導(dǎo)會(huì)影響材料表面的溫度分布、氧化屬性以及燒蝕速率等，因此要準(zhǔn)確預(yù)測(cè)材料的燒蝕響應(yīng)，必須考慮材料熱傳導(dǎo). 材料內(nèi)部的熱傳導(dǎo)遵循傅里葉導(dǎo)熱定律和能量守恒定律，材料內(nèi)部熱傳導(dǎo)的控制方程在直角坐標(biāo)系下可寫為：

式中：t為時(shí)間，T為溫度，ρs為 材料密度，cp為材料的定壓比熱容，λ為材料的熱傳導(dǎo)系數(shù)，q為施加在固體材料邊界上的熱載荷.

材料燒蝕為動(dòng)態(tài)過(guò)程且燒蝕面不斷發(fā)生變化，為準(zhǔn)確預(yù)測(cè)燒蝕響應(yīng)，需要捕捉燒蝕面的位置. 材料的燒蝕速率可按下式計(jì)算得到.

材料表面節(jié)點(diǎn)的移動(dòng)通過(guò)燒蝕速率及時(shí)間步長(zhǎng) ?t可計(jì)算得到，節(jié)點(diǎn)移動(dòng)方向垂直于邊界，節(jié)點(diǎn)移動(dòng)位移為：

式中，δx、δy分別為燒蝕表面節(jié)點(diǎn)沿x、y方向的位移，nx、ny為燒蝕面的內(nèi)法線在x、y方向的分量.

建模中通過(guò)Abaqus網(wǎng)格運(yùn)動(dòng)算法實(shí)現(xiàn)燒蝕表面后退的模擬，利用用戶自定義接口函數(shù)Umeshmotion更新燒蝕表面節(jié)點(diǎn)位置. 網(wǎng)格節(jié)點(diǎn)的移動(dòng)可能導(dǎo)致單元發(fā)生巨大變形，采用ALE（Arbitrary lagrangian-eulerian）網(wǎng)格自適應(yīng)技術(shù)對(duì)模型內(nèi)部網(wǎng)格進(jìn)行重劃分，進(jìn)而避免網(wǎng)格的畸形，最終實(shí)現(xiàn)燒蝕表面后退過(guò)程的模擬.

1.4 流?熱?燒蝕耦合分析策略

高超聲速流動(dòng)與材料的耦合建模采用分區(qū)法實(shí)現(xiàn)，根據(jù)物理空間分為流體部分和固體部分. 流體域與固體域的耦合實(shí)質(zhì)是流體氣動(dòng)加熱問(wèn)題、固體內(nèi)部熱傳導(dǎo)以及壁面燒蝕通過(guò)耦合界面發(fā)生相互作用的物理化學(xué)過(guò)程. 圖2（a）中 ?f和 ?s分別為流體域和固體域，Γ為耦合界面，在耦合界面上，固體向流體提供壁面溫度及位移邊界，而流體向固體提供氣動(dòng)熱載荷. 流場(chǎng)與材料模型的計(jì)算數(shù)據(jù)在耦合界面上反復(fù)交換，通過(guò)Mpcci實(shí)現(xiàn)兩個(gè)區(qū)域非匹配網(wǎng)格間的數(shù)據(jù)傳遞，圖2（b）為非匹配網(wǎng)格間的數(shù)據(jù)傳遞示意.

圖2 流固耦合示意及界面數(shù)據(jù)傳遞. （a）流固耦合示意圖；（b）非匹配網(wǎng)格間的數(shù)據(jù)傳遞Fig.2 Fluid structure coupling and data transfer at the interface:(a) fluid structure coupling; (b) data transfer between unmatched grids

圖3 耦合計(jì)算流程Fig.3 Flow of coupled computing

流?熱?燒蝕耦合分析的具體計(jì)算流程如圖3所示. 首先計(jì)算熱化學(xué)非平衡流場(chǎng)，獲得壁面處的熱流以及各化學(xué)組元的質(zhì)量分?jǐn)?shù)，然后以流場(chǎng)計(jì)算結(jié)果作為邊界條件，進(jìn)行材料內(nèi)部傳熱和壁面燒蝕的計(jì)算，計(jì)算后將壁面溫度、表面位置作為邊界條件反饋回流場(chǎng)，并計(jì)算該時(shí)刻的流場(chǎng)分布情況. 通過(guò)耦合面數(shù)據(jù)的反復(fù)迭代，獲得不同時(shí)刻材料溫度及壁面燒蝕的分布情況.

2 碳/碳復(fù)合材料流?熱?燒蝕耦合模型

對(duì)零攻角碳/碳復(fù)合材料前緣模型燒蝕性能進(jìn)行數(shù)值計(jì)算. 該算例模型選自文獻(xiàn)[24]，模型尺寸可滿足后續(xù)材料燒蝕試驗(yàn)的研究，模型半錐角為10°，前緣半徑為0.025 m，其幾何模型如圖4 所示. 來(lái)流馬赫數(shù)為8，來(lái)流環(huán)境取20 km高空環(huán)境，靜溫217 K，靜壓 5475 Pa，來(lái)流組元 N2、O2質(zhì)量分?jǐn)?shù)分別為77%及23%，碳/碳復(fù)合材料的性能如表1所示.

采用二維模型進(jìn)行計(jì)算，利用Ansys ICEM對(duì)流場(chǎng)模型進(jìn)行網(wǎng)格劃分，網(wǎng)格為四邊形結(jié)構(gòu)化網(wǎng)格，流場(chǎng)壁面邊界層中的網(wǎng)格質(zhì)量直接影響到氣動(dòng)熱環(huán)境的數(shù)值模擬和熱流計(jì)算，對(duì)模型壁面附近網(wǎng)格進(jìn)行加密，以滿足壁面邊界層對(duì)網(wǎng)格的要求，外部流場(chǎng)模型的網(wǎng)格數(shù)量為5251. 材料模型采用Abaqus中的四邊形CPE4T熱力耦合單元，網(wǎng)格數(shù)量為2355. 流場(chǎng)與碳/碳前緣模型網(wǎng)格分布如圖5所示. 初始時(shí)刻材料前緣模型溫度為300 K，耦合計(jì)算時(shí)間步長(zhǎng)為 ?t= 0.01，耦合計(jì)算總時(shí)長(zhǎng)為 30 s，以此來(lái)分析材料模型溫度分布、壁面燒蝕速率以及燒蝕外形的變化.

3 計(jì)算結(jié)果及分析

圖4 碳/碳復(fù)合材料前緣模型Fig.4 Leading edge of carbon/carbon composite

表1 碳/碳材料性能參數(shù)Table 1 Performance parameters of carbon/carbon materials

通過(guò)計(jì)算，可以得到不同時(shí)刻外部流場(chǎng)以及材料模型響應(yīng)的結(jié)果. 在t= 20 s時(shí)，外部流場(chǎng)馬赫數(shù)及溫度的云圖分布如圖6所示，在模型頭部前端形成激波，在駐點(diǎn)區(qū)附近，馬赫數(shù)較小，接近于0，駐點(diǎn)處流場(chǎng)溫度最高可達(dá)到2860 K.

圖7為材料前緣模型燒蝕20 s后與初始狀態(tài)對(duì)比的位置云圖，材料表面產(chǎn)生了明顯的燒蝕后退，由圖可以看出在駐點(diǎn)區(qū)發(fā)生了較為嚴(yán)重的燒蝕，在側(cè)面部分燒蝕量較小，該現(xiàn)象從圖8燒蝕深度沿壁面的分布情況也可看出. 產(chǎn)生這種現(xiàn)象的主要原因是由于是駐點(diǎn)區(qū)材料模型表面溫度較高，化學(xué)反應(yīng)更為活躍，進(jìn)而導(dǎo)致燒蝕氧化速度更快.

圖5 外部流場(chǎng)及材料模型網(wǎng)格劃分. （a）流場(chǎng)網(wǎng)格示意圖；（b）材料前緣模型網(wǎng)格示意Fig.5 External flow field grids and material model:(a) external flow field grids; (b) leading edge model grids

圖6 流場(chǎng)馬赫數(shù)（a）和溫度（b）云圖分布Fig.6 Mach (a) and temperature (b) distribution of the flow field

圖7 材料模型的燒蝕外形Fig.7 Material model ablation profile

圖8 材料模型壁面燒蝕深度分布Fig.8 Distribution of the ablation depth along the wall

通過(guò)計(jì)算可得到不同時(shí)刻熱流密度的分布情況，如圖9所示，駐點(diǎn)區(qū)熱流密度最大，遠(yuǎn)離駐點(diǎn)區(qū)熱流密度隨之減小，隨著燒蝕時(shí)間的推進(jìn)，駐點(diǎn)區(qū)熱流密度減小. 這是由于熱流密度的大小與壁面附近溫度梯度的大小有關(guān)，溫度梯度越大，熱流值也越大. 在初始時(shí)刻，外部流場(chǎng)與材料模型發(fā)生熱交換，流場(chǎng)溫度較高，由于熱傳導(dǎo)存在一定的延遲性，材料壁面溫度較低，此時(shí)壁面附近的溫度梯度較大，導(dǎo)致熱流密度也較大，隨著熱量不斷向材料模型內(nèi)部進(jìn)行傳遞，材料模型的溫度場(chǎng)也隨之逐漸升高，使得溫度梯度變小，熱流密度也就變小. 同時(shí)也說(shuō)明了如果采用非耦合模型，將無(wú)法預(yù)測(cè)到熱流密度的變化，從而高估熱載荷的大小.圖10為材料模型壁面溫度分布情況，由圖可以看出隨著燒蝕時(shí)間的推進(jìn)，駐點(diǎn)區(qū)壁面溫度升高較快，而遠(yuǎn)離駐點(diǎn)區(qū)的壁面溫度升高較慢，同時(shí)也驗(yàn)證了壁面熱流變化的原因.

圖9 不同時(shí)刻壁面熱流分布Fig.9 Heat flow distribution

圖10 不同時(shí)刻外壁面的溫度分布Fig.10 Distribution of wall temperature

圖11為預(yù)測(cè)不同時(shí)刻材料模型的燒蝕外形，隨著燒蝕時(shí)間的推移，模型燒蝕后退深度逐漸增大，10、20、30 s時(shí)刻材料駐點(diǎn)所對(duì)應(yīng)的燒蝕深度分別為 3.42、9.86、17.47 mm，材料模型頭部半徑也逐漸變大，耦合模型考慮了外形變化對(duì)氣動(dòng)熱環(huán)境的影響. 圖12為駐點(diǎn)處燒蝕速率隨著時(shí)間的變化情況，駐點(diǎn)處在開始階段燒蝕速率變化最大，而在之后變化逐漸緩慢，10、20、30 s時(shí)刻材料駐點(diǎn)所對(duì)應(yīng)的燒蝕速率分別為 0.548、0.725、0.795 mm·s?1，這是由于開始階段的高熱流導(dǎo)致材料表面溫度迅速升高，進(jìn)而導(dǎo)致燒蝕速率也迅速變大，隨著燒蝕時(shí)間的推移，材料表面溫度上升，溫度梯度減小，導(dǎo)致熱流密度減小，進(jìn)而表面溫度趨于平穩(wěn)，材料的燒蝕速率也將逐漸趨于平穩(wěn).

圖11 燒蝕外形預(yù)測(cè)Fig.11 Ablation profile prediction

圖12 駐點(diǎn)處燒蝕速率隨時(shí)間的變化Fig.12 Rate of recession at the point of stagnation

通過(guò)以上分析可以看出，在高超聲速氣動(dòng)熱環(huán)境下碳/碳復(fù)合材料前緣模型的頭部區(qū)域?qū)a(chǎn)生一定的燒蝕后退，從而導(dǎo)致外部流場(chǎng)發(fā)生變化，使得氣動(dòng)熱載荷發(fā)生變化. 已有的相關(guān)模型未考慮流固雙向耦合，忽略燒蝕后退對(duì)氣動(dòng)熱載荷的影響，將導(dǎo)致氣動(dòng)熱載荷預(yù)測(cè)產(chǎn)生較大的誤差，從而影響碳/碳復(fù)合材料模型燒蝕響應(yīng)的有效預(yù)測(cè).

4 結(jié)論

（1）考慮熱化學(xué)非平衡效應(yīng)、固體材料傳熱以及材料表面熱化學(xué)燒蝕，建立了高超聲速流場(chǎng)與碳/碳復(fù)合材料燒蝕響應(yīng)的雙向流?熱?燒蝕多場(chǎng)耦合模型，并對(duì)碳/碳復(fù)合材料前緣模型的傳熱及燒蝕響應(yīng)進(jìn)行了預(yù)測(cè).

（2）碳/碳復(fù)合材料前緣模型初始階段駐點(diǎn)區(qū)熱流值最大，隨著燒蝕時(shí)間的推移，材料壁面溫度逐漸升高，駐點(diǎn)區(qū)溫度梯度變小，熱流值也減小.壁面溫度和熱流隨時(shí)間都發(fā)生了顯著的變化，因此采用多場(chǎng)耦合模型可更好的預(yù)測(cè)流場(chǎng)及材料響應(yīng)的變化過(guò)程.

（3）在高超聲速氣動(dòng)熱服役環(huán)境下，碳/碳復(fù)合材料前緣模型駐點(diǎn)區(qū)的溫度較高，材料表面反應(yīng)活躍，燒蝕最為嚴(yán)重，而模型側(cè)面只發(fā)生少量燒蝕，燒蝕前后材料模型外形產(chǎn)生一定的變化，前緣半徑增大.