于勇，張航維

（北京理工大學 宇航學院，北京 100081）

固體火箭發(fā)動機通過噴管將高壓高溫燃氣中的熱能轉換為動能而產(chǎn)生推力[1].因為具有低密度、高比強度、耐熱沖擊和耐燒蝕等一系列優(yōu)異性能[2]，碳基材料(石墨和C/C)常用作固體火箭發(fā)動機的噴管材料.在發(fā)動機工作過程中，噴管長時間處在高溫高壓的流場環(huán)境下，同時受到化學反應燒蝕和凝相粒子的侵蝕作用，引起碳基材料的損失，導致噴管型面衰退[3].其中化學燒蝕是導致噴管燒蝕的重要原因，主要表現(xiàn)為燃氣中的H2O、CO2、OH 等氧化性組分和噴管中碳基材料的化學反應[4].喉部區(qū)域型面的變化直接影響著發(fā)動機性能，因此喉部的燒蝕性能是噴管考核的重要指標[5-6].

在針對化學燒蝕的仿真研究中，BIANCHI 和NASUTI[7]建立了固體火箭發(fā)動機碳基噴管的化學燒蝕模型，證明噴管燒蝕率與燃燒室壓力成線性關系.目前研究認為噴管化學燒蝕主要受化學動力學和組分擴散控制，并在此基礎上發(fā)展出了多種控制機制模型.1985 年KOU 和KESWANI[8]率先提出了計算噴管化學燒蝕速率的雙控制模型.考慮到雙控制模型的局限性，THAKRE 和YANG[9]于2008 年提出了最小機制控制模型，該模型綜合考慮了化學反應速率和燃氣擴散速率的影響，分析了不同Al 質(zhì)量分數(shù)工況下的燒蝕率，結果與試驗吻合較好.張斌等[10]利用多種化學燒蝕機制分析噴管燒蝕過程，發(fā)現(xiàn)計算H2O 引起的燒蝕率時, 應綜合考慮壁面溫度和燃氣組分中H2O 濃度；而計算CO2引起的燒蝕率時,可直接選用化學反應速率控制機制，其研究工作同時表明溫度升高導致燒蝕率增加主要是通過影響化學反應速率實現(xiàn).針對燃氣中參與反應的氧化性組分，張曉光等[11]的研究結果表明H2O 在噴管化學燒蝕過程中占主要地位.魯文濤[12]在模擬噴管化學燒蝕的過程中考慮了H2O 和CO 之間的可逆氣相反應，發(fā)現(xiàn)在預測C/C 喉襯退移規(guī)律時，可以忽略氣相反應帶來的影響.鞏倫坤等[13]以半潛入噴管的收斂段為研究對象，分析碳/酚醛材料的燒蝕規(guī)律與氣流方向的關系，認為當氣流與材料表面平行時，氣相對流引起的化學燒蝕起主要作用.

此外，國內(nèi)外的研究利用試驗分析了影響噴管燒蝕過程的因素.EVANS 等[14]使用固體火箭發(fā)動機模擬裝置研究了金屬推進劑和非金屬推進劑兩種情況下石墨噴管的喉部瞬態(tài)燒蝕特性，發(fā)現(xiàn)采用非金屬推進劑測試的噴管喉部表面更粗糙，被認為是由更顯著的非均相反應造成的.陳博等[15]進行了燃氣發(fā)生器條件下C/C 復合材料噴管的燒蝕性能試驗，結果表明壓強和流速影響了氧化性組分通過邊界層的質(zhì)量傳輸，氧化性組分濃度和溫度的升高均會導致噴管燒蝕率上升.蘇君明等[16]對工程應用的5 種C/C 復合材料和2 種石墨喉襯進行了地面點火試驗，證明噴管熱環(huán)境下的燃氣氧化性組分濃度和燃燒室壓強是影響碳基材料喉襯燒蝕率的主要原因.

在材料燒蝕導致邊界移動方面的研究，國內(nèi)外也取得了一系列成果.CHEN 和MILOS[17]利用TITAN程序和動網(wǎng)格技術成功模擬了流固熱耦合作用下熱防護材料的型面衰退.MENG 等[18]、ZHANG 等[19]利用Abaqus 中的任意拉格朗日?歐拉(ALE)自適應網(wǎng)格分析了多物理場作用下的表面燒蝕現(xiàn)象.薛海峰等[20]以碳/酚醛燃氣舵為研究對象，建立了化學燒蝕模型，并利用彈性光順法和網(wǎng)格重構法模擬了由熱化學燒蝕引起的材料表面退移過程.

綜上所述，目前國內(nèi)外對噴管化學燒蝕的研究多集中在噴管內(nèi)壁面上的化學反應機理，以流場影響噴管內(nèi)部傳熱的單向流固耦合為主，忽略了二者的雙向耦合過程，即噴管型面變化和流場變化在固體火箭發(fā)動機工作期間的相互作用.本文建立了動態(tài)燒蝕耦合模型，綜合考慮了傳熱，化學燒蝕和型面衰退對噴管燒蝕率的影響，實現(xiàn)了噴管型面衰退和噴管流場的雙向流固耦合，著重分析了動態(tài)燒蝕過程中噴管內(nèi)壁面和流場的變化情況.

1 數(shù)學物理模型

1.1 基本假設

本文研究基于以下假設：

①燃氣組分為理想氣體；

②不考慮燃氣中各組分氣體之間的化學反應；

③不考慮凝相顆粒導致的機械侵蝕的影響；

④忽略輻射換熱，忽略重力等體積力的影響.

1.2 控制方程

①氣相控制方程與湍流模型.

噴管中多組分氣體湍流流動的Reynolds 時均方程組通用形式[21]如下：

式中：?為通用變量；T?為輸運系數(shù)；S?為各方程源項.對于連續(xù)、組分、動量和能量方程，上述變量具有特定的形式.

理想氣體狀態(tài)方程為

式中：Yi為某氣體組分i的質(zhì)量分數(shù)；Wi為某氣體組分i的摩爾質(zhì)量.

湍流模型采用Realizable k-epsilon 兩方程模型，采用增強型壁面函數(shù)處理近壁面的流動.

② 固相控制方程.

固相材料瞬時導熱方程為

式 中：ρs為 固 相 材 料的 密 度；Cs為 固 相 材 料 的 比 熱；Ts為固相材料的溫度；λs為固相材料的導熱系數(shù).

③ 氣?固界面邊界條件.

質(zhì)量平衡：

組分平衡：

能量平衡：

式中：ρg為氣體混合物的密度；r˙為固相材料消耗率；Di,m為 某 氣 體組 分i的 擴 散 系 數(shù)；ω˙i為某 氣 體 組 分i的生成或消耗速率；λg為氣體混合物的導熱系數(shù)；qi為某氣體組分i的熱流密度；下標 w表示氣?固交界面.

④ 表面異相反應.

噴管固相材料和燃氣在噴管內(nèi)壁面上發(fā)生反應，該反應的反應速率遵循Arrhenius 定律，見式(7).反應中固相材料的消耗率見式(8).

式中：kj為反應速率；Aj為反應指前因子；Tw為噴管內(nèi)壁面溫度；βj為溫度指數(shù)；Ej為反應活化能 ；R為通用氣體常數(shù)；pj,w為某氧化性組分j在噴管內(nèi)壁面處的分壓；nj為反應壓強指數(shù)；ρs為噴管固體材料密度.

根據(jù)式(7)利用用戶自定義函數(shù)UDF 和DEFINE_SR_RATE 宏定義化學反應速率，建立化學反應模型，完成噴管內(nèi)壁面處的化學反應計算.

⑤化學燒蝕模型.

化學燒蝕模型表征噴管燒蝕率的表現(xiàn)形式，本文采用化學動力學控制模型，見式(9)，該公式表示噴管由于化學反應導致的固相材料消耗率即為噴管燒蝕率.

式中r˙erosion為噴管燒蝕率.

1.3 物理模型和參數(shù)

本文以70-lb BATES 發(fā)動機[22]為原型，模擬計算其噴管內(nèi)流場以及噴管內(nèi)壁面的化學燒蝕情況.發(fā)動機幾何結構如圖1 所示，噴管喉徑為50.8 mm，收斂半角為45°，擴張半角為15°，擴張比為9.5，燃燒時間為5 s.噴管材料為石墨，其物性參數(shù) ρs=1 830 kg/m3，Cs=1 050 J/(kg·K)，λs=70 W/(m·K).

圖1 70-lb 發(fā)動機幾何結構Fig.1 70-lb BATES motor configuration

為精確模擬噴管內(nèi)壁面處的傳熱傳質(zhì)現(xiàn)象，在網(wǎng)格劃分時需保證近壁面處網(wǎng)格的y+<1，近壁面第一層網(wǎng)格高度設定為5×10?7m.

結構網(wǎng)格在保證近壁面處y+值大小的同時，其動網(wǎng)格更新往往會導致遠場處的網(wǎng)格變形較大，難以滿足后續(xù)的計算需求.為適應網(wǎng)格變形并減少網(wǎng)格數(shù)量，本文采用混合網(wǎng)格進行計算，即在流體區(qū)域近噴管內(nèi)壁面處采用結構網(wǎng)格，其余流體區(qū)域和固體區(qū)域采用非結構網(wǎng)格，如圖2 所示.

圖2 噴管計算網(wǎng)格Fig.2 Computational grid for nozzle

動網(wǎng)格更新方式為彈性光順法和網(wǎng)格重構法，利用用戶自定義函數(shù)UDF 和DEFINE_GRID_MOTION宏定義動網(wǎng)格，控制噴管內(nèi)壁面網(wǎng)格節(jié)點的移動，節(jié)點的移動量由計算得到的燒蝕量給出.

圖3 為采用動網(wǎng)格和不采用動網(wǎng)格兩種計算方式下15% Al 質(zhì)量分數(shù)噴管內(nèi)壁面處的y+值分布對比，可以看出混合網(wǎng)格劃分方案在保證網(wǎng)格運動的同時基本滿足了y+值的需求.

圖3 噴管壁面y+Fig.3 y+ value of nozzle wall

噴管入口設為壓力入口、噴管出口設為壓力出口，噴管內(nèi)壁面設為傳熱耦合壁面及表面化學反應發(fā)生面，噴管側壁、外壁設為絕熱，噴管固體域初始溫度設定為300 K.

本文化學燒蝕模型主要考慮C 與氧化性組分H2O、CO2的反應，其反應方程式及動力學參數(shù)見表1[23].

表1 表面反應動力學參數(shù)Tab.1 Kinetic parameters of chemical reaction

1.4 耦合機理和方法

在噴管的化學燒蝕過程中，噴管固體域中的石墨和燃氣發(fā)生化學反應，噴管型面發(fā)生衰退，該過程中流場和噴管型面產(chǎn)生雙向流固耦合作用.耦合計算流程如圖4 所示:

圖4 耦合計算流程Fig.4 Couping calculation scheme

①在初始時刻，給定流場計算的邊界條件；

②求解流場區(qū)域，獲得噴管內(nèi)壁面附近的溫度、組分濃度和組分分壓等參數(shù)；

③根據(jù)步驟②中的參數(shù)計算化學反應速率和燒蝕率，得到噴管內(nèi)壁面上的燒蝕率分布，當前計算時間為t，計算時間步長為 ?t；

④利用動網(wǎng)格技術，將噴管內(nèi)壁面上的網(wǎng)格節(jié)點以對應的燒蝕率進行移動；

⑤重復步驟②～④，直至完成所有計算.

2 計算結果及分析

2.1 Al 質(zhì)量分數(shù)對燒蝕的影響

表2 為不同Al 質(zhì)量分數(shù)時的噴管入口條件，包括燃氣壓強、溫度及組分質(zhì)量分數(shù)[9].圖5 為計算得到的t=5 s 時不同Al 質(zhì)量分數(shù)工況下噴管內(nèi)壁面燒蝕率沿軸向的分布結果，其中喉部的燒蝕率與試驗值[22]吻合較好，具體結果見表3.

表2 不同Al 質(zhì)量分數(shù)工況下噴管入口條件Tab.2 Nozzle inlet conditions with different aluminum mass fraction

表3 不同Al 質(zhì)量分數(shù)工況下喉部燒蝕率與試驗值對比Tab.3 Comparison between calculated and measured throat erosion rates

圖5 不同Al 質(zhì)量分數(shù)工況下噴管壁面燒蝕率分布Fig.5 Erosion rate of under different aluminum mass fraction

結果表明，噴管燒蝕率隨推進劑中Al 質(zhì)量分數(shù)的增加而降低，其主要原因是燃氣中參與化學反應的氧化性組分H2O 和CO2濃度的降低.噴管燒蝕主要分為化學燒蝕和顆粒侵蝕兩部分，推進劑中Al 質(zhì)量分數(shù)的增加導致燃氣中Al2O3顆粒濃度增加，進而增強了Al2O3顆粒對噴管內(nèi)壁面的機械侵蝕強度，而實際的噴管燒蝕率卻不升反降，由此證明了噴管燒蝕中顆粒侵蝕所占比重較小，由此證明了化學燒蝕才是導致噴管燒蝕的主要原因[9].

圖6 為15%的Al 質(zhì)量分數(shù)工況下燒蝕率隨時間變化的分布情況.隨著燒蝕過程的進行，燒蝕率的上升變慢，這主要是由于噴管內(nèi)壁面溫度上升，使得噴管內(nèi)壁面和燃氣之間的熱交換趨于平衡，導致燃氣向固相材料的對流換熱減弱.

圖6 不同時刻下15%Al 質(zhì)量分數(shù)的噴管燒蝕率分布Fig.6 Erosion rate with 15% aluminum mass fraction at different times

圖7 和圖8 分別顯示了15%的Al 質(zhì)量分數(shù)工況下t=5 s 時的噴管壓力分布和速度分布.圖9 對比了15%的Al 質(zhì)量分數(shù)工況下不同時刻的噴管溫度分布，顯示了噴管內(nèi)部的傳熱過程.

圖7 15%Al 質(zhì)量分數(shù)的噴管壓力分布 (單位：MPa)Fig.7 Pressure distribution of nozzle with 15% aluminum mass fraction(unit: MPa)

圖8 15%Al 質(zhì)量分數(shù)的噴管馬赫數(shù)分布Fig.8 Mach number distribution of nozzle with 15% aluminum mass fraction

圖9 15%Al 質(zhì)量分數(shù)的噴管溫度分布(單位：K)Fig.9 Temperature distribution of nozzle with 15% aluminum mass fraction (unit: K)

2.2 動態(tài)燒蝕過程

在噴管工作過程中，由于內(nèi)壁面發(fā)生化學燒蝕，造成了固體域材料損失和噴管型面變化，研究表明噴管喉部面積增加5%就會對發(fā)動機的性能造成巨大影響[24]，因此想要深入了解噴管的化學燒蝕過程，需結合噴管型面衰退進一步分析.非定常工況下噴管受到化學燒蝕導致流場和內(nèi)壁面相互耦合作用，在該作用下流場和內(nèi)壁面發(fā)生變化的過程即為噴管的動態(tài)燒蝕過程.本文通過耦合化學燒蝕模型和動網(wǎng)格模型，實現(xiàn)了噴管的動態(tài)燒蝕.

2.2.1 噴管燒蝕率對比

選用表2 中的工況分別模擬噴管的動態(tài)燒蝕過程.圖10～12 分別是15% Al、18% Al、21% Al 質(zhì)量分數(shù)工況下t=5 s 時動態(tài)燒蝕對噴管內(nèi)壁面燒蝕率的影響，并與未采用動網(wǎng)格的計算結果進行對比，具體結果見表4.圖13 是計算得到的t=5 s 時不同Al 質(zhì)量分數(shù)工況下動態(tài)燒蝕噴管內(nèi)壁面燒蝕率沿軸向的分布結果.

表4 動態(tài)燒蝕工況下喉部燒蝕率與試驗值對比Tab.4 Comparison between calculated and measured throat erosion rates under dynamic erosion

圖10 15% Al 質(zhì)量分數(shù)的噴管燒蝕率分布對比Fig.10 Comparison of erosion rate with 15% aluminum mass fraction

圖11 18% Al 質(zhì)量分數(shù)的噴管燒蝕率分布對比Fig.11 Comparison of erosion rate with 18% aluminum mass fraction

圖12 21% Al 質(zhì)量分數(shù)的噴管燒蝕率分布對比Fig.12 Comparison of erosion rate with 21% aluminum mass fraction

圖13 不同Al 質(zhì)量分數(shù)工況下動態(tài)燒蝕噴管壁面燒蝕率分布Fig.13 Dynamic erosion rate of under different aluminum mass fraction

可以看出，3 種工況下采用動網(wǎng)格的喉部燒蝕率都更貼近試驗值，驗證了本文動態(tài)燒蝕計算方法的有效性.同時3 種工況下采用動網(wǎng)格的噴管燒蝕率峰值都較高.圖14 為15%Al 質(zhì)量分數(shù)工況下2 種計算方式的噴管內(nèi)壁面溫度分布.對比不采用動網(wǎng)格的計算方式，可以看出采用動網(wǎng)格的噴管由于噴管型面和流場的變化，導致噴管喉部及喉部上游區(qū)域溫度較高，兩種計算方式的溫度差最高可達40 K，溫度的升高增強了化學反應速率，導致了較高的峰值燒蝕率.

圖14 15% Al 質(zhì)量分數(shù)的噴管內(nèi)壁面溫度分布對比Fig.14 Comparison of temperature distribution of nozzle with 15% aluminum mass fraction

2.2.2 噴管型面變化和溫度分布對比

圖15(a)和圖15(b)分別為截取喉部上游區(qū)域和下游區(qū)域的噴管內(nèi)壁面型面變化趨勢.由于燒蝕率在噴管上的分布，使得噴管上游型面變化最大.可以看出t=1 s 的初始時刻噴管型面變化較小，隨后逐漸加大.這是由于初始時刻噴管內(nèi)壁面溫度較低，化學反應速率和燒蝕率較小，隨著噴管內(nèi)壁面溫度的上升，化學反應速率加快，進而增大了單位時間內(nèi)的燒蝕量.噴管內(nèi)壁面溫度變化見圖16.

圖15 15% Al 質(zhì)量分數(shù)的噴管喉部上游型面及下游型面變化Fig.15 Variation of upstream and downstream profile of throat with 15%aluminum mass fraction

圖16 不同時刻下15% Al 質(zhì)量分數(shù)的噴管內(nèi)壁面溫度分布Fig.16 Temperature distribution of nozzle with 15% aluminum mass fraction at different times

2.2.3 噴管壓力變化

圖17 和圖18 分別是15% Al 質(zhì)量分數(shù)工況下沿噴管內(nèi)壁面和沿喉部徑向的壓力分布.兩種計算方式下噴管收斂段的壓力差異較小，而相比于未采用動網(wǎng)格的計算方式，動態(tài)燒蝕的噴管在喉部及下游區(qū)域的壓力更高，沿喉部徑向的最大壓力差約為0.3 MPa.這種壓力變化的差異可以根據(jù)圖10 的燒蝕率分布分析得出，圖10 中15% Al 質(zhì)量分數(shù)工況下兩種計算方式在收斂段的燒蝕率曲線較為吻合，采用動網(wǎng)格的噴管燒蝕率在喉部及下游區(qū)域較高.噴管型面衰退和流場變化的耦合作用導致了兩種計算方式下燒蝕率的差異，進而產(chǎn)生了不同的壓力分布.

圖17 15% Al 質(zhì)量分數(shù)的噴管內(nèi)壁面壓力分布對比Fig.17 Comparison of pressure distribution of nozzle with 15% aluminum mass fraction

圖18 15% Al 質(zhì)量分數(shù)的噴管喉部徑向壓力分布對比Fig.18 Radial distribution of pressure at the nozzle throat with 15% aluminum mass fraction

2.2.4 噴管喉部組分變化

圖14 的溫度分布表明，考慮型面衰退影響的噴管在喉部的溫度更高，導致喉部化學反應速率加快，同時在喉部壁面附近形成了組分濃度梯度，組分的濃度分布同時受到化學反應速率和擴散速率的影響.圖19 是15% Al 質(zhì)量分數(shù)工況下沿噴管喉部徑向的組分分布，包括參與反應的H2O 和CO2.可以看出，動態(tài)燒蝕的噴管喉部壁面附近H2O 和CO2的濃度偏低，質(zhì)量分數(shù)分別為0.001 7 和0.003 3，而不采用動網(wǎng)格計算得到的H2O 和CO2的質(zhì)量分數(shù)分別為0.001 6 和0.002 5.這是由于高溫狀態(tài)下H2O 和CO2的擴散速率大于化學反應的消耗速率，2 種反應物的濃度分布主要受擴散速率的影響.

圖19 15% Al 質(zhì)量分數(shù)的噴管喉部徑向組分分布對比Fig.19 Radial distribution of species at the nozzle throat with 15% aluminum mass fraction

3 結 論

通過本文的研究，可以得出以下結論：

①隨著推進劑中Al 質(zhì)量分數(shù)的增加，噴管燃氣中參與化學反應的氧化性組分H2O 和CO2的濃度減小，噴管燒蝕率隨之降低；

②建立的動態(tài)燒蝕模型綜合考慮了傳熱、化學反應和噴管衰退等現(xiàn)象.相比于未采用動網(wǎng)格的計算方式，采用動態(tài)燒蝕模型得到的計算結果更貼近試驗值，能夠有效模擬噴管的動態(tài)化學燒蝕過程；

③噴管內(nèi)壁面的衰退使得噴管流場發(fā)生變化，導致采用動網(wǎng)格技術計算的噴管峰值燒蝕率更高，喉部附近的溫度和壓力更大.喉部附近的溫度差最高可達40 K，沿喉部徑向的最大壓力差約為0.3 MPa；

④受高溫狀態(tài)下化學反應速率和擴散速率的影響，采用動態(tài)燒蝕模型計算的噴管在喉部壁面附近的H2O 和CO2的濃度更高.