代曉松,王一白，劉 宇,常 桁,趙宇輝

(1.北京航空航天大學 宇航學院， 北京 100191; 2.北京機電工程研究所， 北京 100074)

固體軌控發(fā)動機是推力可調(diào)針栓噴管(Pintle Nozzle)技術(shù)的具體應(yīng)用，是一種喉栓式固體火箭發(fā)動機[1]。喉栓式固體火箭發(fā)動機能夠?qū)崿F(xiàn)固體發(fā)動機的可控，即可以根據(jù)實際需要調(diào)節(jié)推力的大小，從多方面提高了飛行器的性能，但隨著高能推進劑和高壓強發(fā)動機技術(shù)的發(fā)展，與常規(guī)固體火箭發(fā)動機一樣，喉栓式固體火箭發(fā)動機仍然面臨著嚴峻的噴管燒蝕問題。(

碳基材料(石墨和C/C)噴管密度低，升華溫度較高，具有優(yōu)異的抗熱振性能，其燒蝕的物理化學過程極為復雜，包含熱化學燒蝕、燃氣機械沖刷及凝相粒子侵蝕等[2]。其中最主要的是熱化學燒蝕，為噴管內(nèi)壁面碳與燃氣氧化性組分異相化學反應(yīng)所致的噴管內(nèi)壁面退移，其速率由異相反應(yīng)動力學速率與燃氣氧化性組分向噴管內(nèi)壁的擴散速率共同決定[3]。

20世紀80年代，Keswani與Kuo[4,5]提出了較為完備的熱化學燒蝕理論模型，包括氣相守恒方程、固相導熱方程及氣-固界面邊界條件。近年來，Acharya與Kuo[6,7]發(fā)展了該模型，應(yīng)用異相反應(yīng)動力學的研究成果開發(fā)出GNEM程序(Graphite Nozzle Erosion Minimization)。Thakre與Yang[2]則采用標準k-ε兩方程湍流模型結(jié)合雙層模型的壁面函數(shù)法模擬氣相湍流流動，并取實際燒蝕率為兩種控制機制下燒蝕率的最小值，克服了對反應(yīng)級數(shù)為1的限制性要求。

為了準確預(yù)示固體軌控C/C復合材料針栓噴管在固體高溫高壓燃氣環(huán)境的傳熱及燒蝕規(guī)律，本文綜合考慮了燃氣流動、燃氣向固壁傳熱、固相導熱、氧化性組分向壁面擴散、壁面異相化學反應(yīng)等過程，通過Fluent流體計算軟件進行了流固耦合仿真，得到了固體軌控發(fā)動機碳基材料針栓噴管的瞬態(tài)溫度分布，研究了噴管與閥桿壁面燒蝕率分布規(guī)律以及燃燒室壓強和氧化性組分對燒蝕率的影響規(guī)律。

1 數(shù)學物理模型

1.1 控制方程

假設(shè)燃氣為凍結(jié)流，組分為理想氣體，組分擴散遵守Fick定律，不考慮凝相粒子的作用，不考慮輻射換熱的影響，不存在體積力。

1) 氣相控制方程與湍流模型[8]

組分擴散方程：

(1)

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

(2)

湍流模型采用Realizablek-ε兩方程模型，近壁處理采用增強型壁面函數(shù)。

2) 固相控制方程

固相材料二維軸對稱瞬時導熱方程：

(3)

3) 氣-固界面邊界條件

質(zhì)量平衡：

(4)

組分平衡：

(5)

能量平衡：

(6)

4) 熱化學反應(yīng)速率

一般固體推進劑燃氣中的主要組分為H2O、CO2、HCl、N2和H2，此外還有少量的OH-、H+，而O+2原子和O2的成分由于燃氣基本屬于富燃氣體可以忽略。由于在任何合理的表面溫度下碳的蒸汽壓都很低，因此燒蝕反應(yīng)發(fā)生在噴管內(nèi)壁面上，為表面異相化學反應(yīng)。綜上所述，主要是C與燃氣氧化性組分H2O、CO2、OH反應(yīng)形成CO[9]，燒蝕反應(yīng)方程式及其動力參數(shù)見表1。

表1 燒蝕反應(yīng)方程式及其動力參數(shù)

反應(yīng)速率根據(jù)Arrhenius 定律計算，由反應(yīng)所致的壁面退移速率就是燒蝕率：

(7)

耦合求解方程(1)～(7)，即可獲得湍流邊界層內(nèi)各物理量的分布，并求得燒蝕率。

1.2 算例驗證

為驗證以上計算方法的有效性與準確性，針對Geisler的70-lb BATES 發(fā)動機石墨噴管進行了流固耦合傳熱燒蝕計算，并與試驗數(shù)據(jù)進行了比較。

圖1為70-lb BATES 發(fā)動機。其噴管為石墨材料，推進劑AP/HTPB/Al，端-側(cè)面燃燒藥柱，內(nèi)徑117 mm，外徑305 mm，長508 mm，燃燒時間不大于5 s。噴管喉徑50.8 mm，收斂角45°，擴張角15°，擴張比9.5。石墨材料物性參數(shù)ρ=1830 kg/m3，cp=1 050 J/(kg·K)，λ=70 W/(m·K)。

圖1 70-lb BATES 發(fā)動機結(jié)構(gòu)

圖2為計算得到的噴管燒蝕率沿軸向的分布和Geisler 的喉部燒蝕率測試結(jié)果[10]。

圖2 70-lb BATES發(fā)動機噴管燒蝕率沿軸向分布

計算值與試驗值見表2，二者吻合較好，驗證了燒蝕計算方法的有效性。

表2 喉部燒蝕率計算值與試驗值

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

為減小計算量，本文僅考慮噴管開度為100%時的非穩(wěn)態(tài)燒蝕，將計算模型簡化為二維軸對稱模型。簡化后的計算模型如圖3所示。推進劑選擇HTPB/Al，噴管喉徑6 mm，收斂角40°，擴張角15°，擴張比13，燃燒時間7 s，時間步長取1×10-4s。C/C復合材料物性參數(shù)ρ=1 900 kg/m3，cs=1 182 J/(kg·K)，λs=98.67 W/(m·K)。噴管入口燃氣溫度 3 266 K，壓強8 MPa。

圖3 噴管二維軸對稱傳熱-燒蝕模型

對以上二維軸對稱模型進行網(wǎng)格劃分，如圖4。

使流場近壁面第一層網(wǎng)格y+<1，以精確模擬湍流邊界層內(nèi)的傳熱傳質(zhì)。噴管入口設(shè)為壓力入口、噴管出口設(shè)為壓力出口，噴管內(nèi)壁及閥桿壁面設(shè)為傳熱耦合壁面及異相反應(yīng)發(fā)生壁面，噴管側(cè)壁、外壁設(shè)為絕熱。

圖4 噴管傳熱-燒蝕計算網(wǎng)格

2 計算結(jié)果及分析

2.1 噴管計算結(jié)果分析

圖5為不同時刻噴管內(nèi)壁溫度沿軸向的分布?？梢姡跍卦谏仙倪^程中，上升速率逐漸減小，這是由于隨噴管內(nèi)壁溫度升高，燃氣向固相材料的對流換熱不斷減弱。

圖5 不同時刻噴管內(nèi)壁溫度

圖6為不同時刻噴管幾何喉部區(qū)域(X=2.5 mm)溫度隨工作時間的變化曲線。由于噴管固體壁面較厚，因此在 7 s的工作時間內(nèi)尚未完全達到穩(wěn)態(tài)傳熱。

圖6 噴管喉部溫度隨時間的變化曲線

圖7、圖8為噴管內(nèi)壁熱流密度曲線和燒蝕率沿軸向的分布曲線。可以看出，燒蝕率曲線與熱流密度曲線形狀基本一致，表明燒蝕率的分布遵循熱流密度的分布規(guī)律。分析可知，燒蝕率和熱流密度在噴管節(jié)流面區(qū)域較高，在遠離節(jié)流面區(qū)域之后則迅速下降。這也再次說明位于噴管收斂段與閥桿頭部之間的節(jié)流面，才是推力器的實際喉部位置。燃氣在沿著噴管軸線方向向下游流動的過程中，繼續(xù)擴壓增速，速度持續(xù)加大，壓強和溫度快速降低，噴管內(nèi)壁面燒蝕率和熱流密度隨之迅速降低。

圖7 噴管內(nèi)壁面熱流密度曲線

圖8 噴管內(nèi)壁面燒蝕率隨時間的變化曲線

另由圖8可看出，燒蝕率的峰值出現(xiàn)在環(huán)狀節(jié)流面的上游位置，而不是等效喉部的錐形截面上。燒蝕率的這一分布特點是由噴管湍流邊界層流動及傳熱傳質(zhì)特性所決定的。從收斂段開始，邊界層將逐漸減薄，傳熱速率不斷增大；而在節(jié)流面上游某處，邊界層外緣將達到聲速，該處邊界層最薄，燃氣質(zhì)量流率、熱流密度、組分質(zhì)量擴散流率均達到峰值，材料的熱化學燒蝕率也最大；之后在擴張段，邊界層又將迅速增厚，傳熱速率也急劇下降。

圖9為燒蝕最大值位置(X=4 mm)的燒蝕速率隨工作時間的變化曲線。隨喉襯內(nèi)壁溫度升高，燃氣向喉襯的對流換熱不斷減弱，使得內(nèi)壁溫及燒蝕率的上升速度逐漸減小。與溫度變化相似，由于噴管徑向厚度較大，導致當工作時間達到7 s時燒蝕率仍未達到穩(wěn)態(tài)變化。分析整個傳熱燒蝕過程，燒蝕率隨傳熱過程的進行逐漸增大，當時間為7 s時，達到工作期間的最大燒蝕率1.375 mm/s。

2.2 閥桿計算結(jié)果分析

圖10為不同時刻閥桿內(nèi)壁面溫度沿軸向的分布曲線。與噴管溫度變化相似，閥桿壁面溫度在上升的過程中呈現(xiàn)出上升速率逐漸減小的特點。閥桿內(nèi)壁面溫度的最大值位于桿頭和節(jié)流面處，一方面因為桿頭處流場存在的回流區(qū)(見圖11)，溫度較高，另一方面是由于節(jié)流面流動速度快，加劇了熱對流，使熱流密度增大。

圖9 喉部入口處燒蝕率隨時間的變化曲線

圖10 不同時刻閥桿壁面溫度分布曲線

圖11 閥桿頭部回流區(qū)示意圖

圖12和圖13分別為閥桿內(nèi)壁熱流密度和燒蝕率沿軸線的分布曲線。分析可知，由于閥桿頭部頂端的回流區(qū)存在著高溫高壓，因此該區(qū)域的燒蝕率較為嚴重(1.79 mm/s)；而與噴管燒蝕相似，在節(jié)流面上游位置，燃氣速度增加，傳熱速率不斷增大，材料的熱化學燒蝕率也隨之不斷增大。從圖中也可以看出，閥桿的壁面燒蝕率在t=7 s時達到1.79 mm/s，可見7 s的工作時間內(nèi)，C/C材料閥桿將由于燒蝕而使得型面失真，喪失原有功能，這是設(shè)計中必須考慮的問題。

2.3 熱化學燒蝕影響因素分析

本文主要研究了燃燒室壓強和燃氣的氧化組分含量對熱化學燒蝕率的影響。由于噴管和閥桿的計算結(jié)果所反映出的基本規(guī)律相同，下文僅就噴管固相材料的熱化學燒蝕展開分析。

1) 壓強對熱化學燒蝕的影響

本文計算了pc分別為2、4、6、8 MPa時的噴管燒蝕率及內(nèi)壁熱流密度，結(jié)果見圖14、圖15。由于發(fā)動機燃燒室壓強變化對燃燒室溫度及燃氣組分質(zhì)量分數(shù)的影響有限，因此對上述工況的計算中皆采用相同的入口溫度、組分質(zhì)量分數(shù)。

由圖中看出，入口總壓導致流場內(nèi)壓強分布變化是影響燒蝕率的重要因素。噴管壓強升高時，流量密度增加，對流換熱及質(zhì)量傳輸增強；流動雷諾數(shù)增加，邊界層變薄，燃氣向壁面?zhèn)鳠?、氧化性組分向壁面擴散的阻力降低；壁面異相化學反應(yīng)速率增大。這些因素都將導致較高的固相材料燒蝕率。由圖可以看出，內(nèi)壁面熱流密度與燒蝕率分布趨勢始終保持一致。隨著噴管內(nèi)流場燃氣壓強的增大，燒蝕率和壁面熱流密度保持線性增長。

圖12 閥桿壁面熱流密度分布曲線

圖13 閥桿壁面燒蝕率分布曲線

圖14 不同壓強下噴管壁面熱流密度

圖15 不同壓強下噴管壁面熱化學燒蝕率

圖16展示了噴管節(jié)流面上游位置(最大燒蝕點)的熱化學燒蝕率隨壓強的變化情況，可以看出，燒蝕率與壓強近似成線性關(guān)系。

圖16 壓強對噴管喉部燒蝕率的影響

2) 氧化性組分對熱化學燒蝕的影響

圖17為噴管幾何喉部截面燒蝕反應(yīng)相關(guān)組分質(zhì)量分數(shù)沿徑向的分布。可以看出，由于壁面化學反應(yīng)的存在，C氧化使CO和H2的質(zhì)量分數(shù)增加，而氧化性組分H2O、CO2的質(zhì)量分數(shù)則因反應(yīng)消耗而降低。此外，由于C與H2O燒蝕反應(yīng)的指前因子遠大于CO2，因此H2O的反應(yīng)消耗也更大，是決定燒蝕的主要氧化性組分。

圖17 組分質(zhì)量分數(shù)沿喉部的徑向分布

3 結(jié)論

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

1) 碳基材料針栓噴管的燒蝕主要是熱化學燒蝕，噴管內(nèi)壁及閥桿壁面的溫度及燒蝕率分布均遵循壁面熱流密度的分布規(guī)律，并且在等效節(jié)流面上游處達到最大值。

2) 閥桿頭部位置存在高溫高壓回流區(qū)，導致該區(qū)域的燒蝕率較高。

3) H2O是決定噴管內(nèi)壁及閥桿燒蝕的主要氧化性組分，因為其在燃氣中質(zhì)量分數(shù)較高，并且與C的異相反應(yīng)速率更大。

4) 噴管內(nèi)壁及閥桿燒蝕率隨燃燒室壓強升高而增加，并且與壓強近似成線性關(guān)系。

[1] NAPIOR J,GARMY V.Controllable Solid Propulsion for Launch Vehicle and Spacecraft Application [C].57thInternational Astronautical Congress.Spain,2006.

[2] THAKRE P.Chemical Erosion of Graphite and Refractory Metal Nozzles and Its Mitigation in Solid-Propellant Rocket Motors[D].The Pennsylvania State University,2008,5-6,10.

[3] THAKRE P,YANG V.Chemical Erosion of Carbon-Carbon/Graphite Nozzles in Solid-Propellant Rocket Motors[J].Journal of Propulsion and Power,2008,24(4):822-833.

[4] KUO K K,KESWANI S T.A Comprehensive Theoretical Model for Carbon-Carbon Composite Nozzle Recession[J].Combustion Science and Technology,1985,42:145-164.

[5] KESWANI S T,KUO K K.Validation of an Aerothermochemical Model for Graphite Nozzle Recession and Heat-Transfer Process[J].Combustion Science and Technology,1986,47:177-192.

[6] ACHARYA R,KUO K K.Effect of Pressure & Propellant Composition on Graphite Rocket Nozzle Erosion Rate[J].Journal of Propulsion and Power,2007,23(6):1242-1254.

[7] ACHARYA R,KUO K K.Numerical Simulation of Graphite Nozzle Erosion with Parametric Analysis[R].AIAA 2010-6846.

[8] 周力行.湍流氣粒兩相流動和燃燒的理論與數(shù)值模擬[M].北京：科學出版社,1994:154-161.

[9] GEISLER R L.The Relationship Between Solid Propellant Formulation Variables and Nozzle Recession Rates[C].JANNAF(15th)Rocket Nozzle Technology Subcommittee Workshop,1978.

[10] GEISLER R L,BECKMAN C W.The History of the BATES Motors at the Air Force Rocket Propulsion Laboratory [R].AIAA 98-3981.