周志壇，丁逸夫，樂貴高，梁曉揚(yáng)，孫培杰，盛英華

(1.南京理工大學(xué)機(jī)械工程學(xué)院，南京 210094；2.上海宇航系統(tǒng)工程研究所，上海 201109)

0 引 言

近年來，隨著載人登月、空間站組建、深空探測(cè)等一系列重大航天活動(dòng)相繼開展，我國(guó)加大了對(duì)重型運(yùn)載火箭及其動(dòng)力系統(tǒng)的研究工作，確定了液氧煤油發(fā)動(dòng)機(jī)和液氫液氧發(fā)動(dòng)機(jī)是運(yùn)載火箭動(dòng)力系統(tǒng)的最佳選擇[1-3]。液體運(yùn)載火箭發(fā)動(dòng)機(jī)在工作時(shí)，由于外界環(huán)境壓力過低，燃?xì)饬鬟M(jìn)入外界環(huán)境后急劇膨脹，在火箭底部形成回流，并對(duì)箭體底部形成對(duì)流加熱效應(yīng)，同時(shí)，高溫燃?xì)庵袊娚錈霟岬腃O2和H2O等混合氣流對(duì)火箭底部形成輻射加熱效應(yīng)[4]。火箭底部相當(dāng)于背風(fēng)面，容易受發(fā)動(dòng)機(jī)噴流回流造成的對(duì)流加熱與輻射加熱耦合作用，導(dǎo)致溫度迅速升高。對(duì)火箭底部熱環(huán)境估計(jì)過低會(huì)給火箭整體的安全性帶來極大威脅，甚至誘發(fā)爆炸等重大事故造成飛行失敗，而估計(jì)過高又會(huì)導(dǎo)致熱防護(hù)設(shè)備保守設(shè)計(jì)，增加發(fā)射成本[5-7]，因此，對(duì)液體運(yùn)載火箭底部進(jìn)行熱環(huán)境分析研究成為當(dāng)前的重中之重。

國(guó)外已有大量學(xué)者對(duì)火箭底部熱環(huán)境開展研究工作。1970年，Kramer[8]對(duì)泰坦III底部熱環(huán)境進(jìn)行試驗(yàn)測(cè)試，結(jié)果顯示在剛起飛1 s內(nèi)輻射熱流是總熱流的主要來源，隨著飛行高度上升，輻射熱流減少，對(duì)流熱流升高；1994年，Reijasse等[9]對(duì)阿里安5進(jìn)行數(shù)值模擬和縮比試驗(yàn)對(duì)比研究，得出影響對(duì)流熱流的主要因素是飛行高度和噴管數(shù)量以及噴管排列方式；2007年，Negishi等[10]對(duì)H-IIA底部熱環(huán)境進(jìn)行數(shù)值模擬，結(jié)果顯示底部熱流峰值在30 km 高空附近。

近年來國(guó)內(nèi)也有部分學(xué)者開展了噴流試驗(yàn)以及數(shù)值仿真方面的工作。史亞男[11]采用有限體積法求解RANS方程，湍流模型選取RNGk-ε兩方程模型對(duì)火箭底部熱環(huán)境進(jìn)行仿真；周一鵬等[12]利用反向蒙特卡羅法計(jì)算火箭發(fā)動(dòng)機(jī)噴管內(nèi)的輻射熱流；戴欣怡[13]采用標(biāo)準(zhǔn)k-ε兩方程模型分析了液體發(fā)動(dòng)機(jī)噴流對(duì)底部熱環(huán)境的影響。

相比于國(guó)外對(duì)火箭底部熱環(huán)境研究程度，國(guó)內(nèi)目前還處于早期階段，尤其是我國(guó)液體運(yùn)載火箭飛行試驗(yàn)數(shù)據(jù)極少，國(guó)內(nèi)學(xué)者對(duì)火箭底部熱環(huán)境多為數(shù)值模擬方面的研究，缺乏對(duì)比分析導(dǎo)致仿真方法的有效性難以得到檢驗(yàn)。并且由于計(jì)算資源有限，所以大多數(shù)學(xué)者對(duì)數(shù)值模型進(jìn)行了不同程度簡(jiǎn)化，如幾何模型只考慮四分之一、網(wǎng)格無(wú)邊界層、熱輻射計(jì)算中采用了精度較低的P-1模型等，這些都在很大程度上影響了計(jì)算精度。

隨著計(jì)算機(jī)運(yùn)算速度和內(nèi)存規(guī)模飛速發(fā)展，為開展含邊界層網(wǎng)格的結(jié)構(gòu)化全箭三維流動(dòng)數(shù)值模擬創(chuàng)造了有利條件。本文對(duì)液體運(yùn)載火箭底部熱環(huán)境進(jìn)行數(shù)值模擬研究，在與試驗(yàn)測(cè)試完成對(duì)比驗(yàn)證基礎(chǔ)上，通過大量研究分析不同飛行高度環(huán)境下，火箭底部熱流的變化規(guī)律。

1 飛行試驗(yàn)

試驗(yàn)測(cè)試對(duì)象為某三級(jí)液體運(yùn)載火箭，火箭共分為三個(gè)子級(jí)，由于液氧煤油發(fā)動(dòng)機(jī)密度比沖高，因此一子級(jí)推進(jìn)劑選擇為液氧和煤油。在火箭底部放置兩個(gè)熱流傳感器，傳感器位置如圖1所示。

測(cè)試結(jié)果為該火箭按照實(shí)際軌道飛行過程中，箭體底部?jī)蓚€(gè)熱流傳感器反饋的數(shù)據(jù)。

圖1 火箭底部熱流傳感器位置Fig.1 The heat flux sensor of rocket tail compartment

2 計(jì)算模型

2.1 計(jì)算網(wǎng)格

采用多重網(wǎng)格方法對(duì)三維模型進(jìn)行網(wǎng)格劃分，該方法通過在不同疏密網(wǎng)格下求解原方程組，平順高、低頻誤差，加速收斂過程。

為確保激波能夠完整描述，對(duì)噴管及噴流區(qū)域網(wǎng)格進(jìn)行加密處理，整個(gè)計(jì)算域均采用結(jié)構(gòu)網(wǎng)格以保證網(wǎng)格的正交性及光順性。在數(shù)值計(jì)算同時(shí)，為防止邊界處的溫度、壓力等物理量變化劇烈，本文對(duì)邊界層計(jì)算域的網(wǎng)格進(jìn)一步加密，確保壁面附近網(wǎng)格Y+值小于3，具體如圖2所示。

圖2 計(jì)算網(wǎng)格Fig.2 Computational mesh

2.2 數(shù)值方法

2.2.1控制方程

基于液體運(yùn)載火箭多組分燃?xì)鉂M足連續(xù)介質(zhì)假設(shè)和理想氣體狀態(tài)方程的特點(diǎn)，建立三維多組分可壓縮Navier-Stokes方程，包括質(zhì)量、動(dòng)量與能量輸運(yùn)的控制方程，形式如下。

多組分輸運(yùn)方程為：

(1)

式中：Yl為組分l的質(zhì)量分?jǐn)?shù)，Rl為組分l在化學(xué)反應(yīng)后的凈生成率，Sl為自定義源項(xiàng)的離散相引起的生成率。

組分?jǐn)U散通量Jl可表示為：

(2)

式中：Dl，m為組份l的質(zhì)量耗散系數(shù)，DT，l為組份l的熱擴(kuò)散系數(shù)。

在直角坐標(biāo)系中，單一組分l的可壓縮N-S方程的守恒形式為：

(3)

(4)

(5)

(6)

式(3)～(6)中:U為流動(dòng)變量；F,G,H為氣流通量矢量;Fv,Gv,Hv為黏性通量矢量;K為熱傳導(dǎo)系數(shù)；T為環(huán)境溫度；p,ρ,e,τ,μ分別為壓力、密度、比動(dòng)能、應(yīng)力、黏性系數(shù);u,v,w分別為速度在x,y,z方向上的分量。

2.2.2空間離散格式

本運(yùn)載火箭采用了大推力發(fā)動(dòng)機(jī)設(shè)計(jì)方案，燃燒室總壓、總溫高，噴流速度為典型超聲速，高空階段箭體以超聲速飛行，來流與噴流強(qiáng)烈干擾而形成復(fù)雜的發(fā)動(dòng)機(jī)噴流效應(yīng)，二階迎風(fēng)TVD格式收斂速度快，數(shù)值穩(wěn)定性高，由于設(shè)計(jì)了魯棒性強(qiáng)的通量限制器，不僅能求得高分辨率的復(fù)雜波系結(jié)構(gòu)，而且抑制強(qiáng)間斷解的非物理震蕩，擬選用二階迎風(fēng)TVD格式計(jì)算全箭這類大型噴流模型。

將含箭體的整個(gè)計(jì)算域劃分為六面體的結(jié)構(gòu)化網(wǎng)格，采用有限體積法對(duì)每個(gè)體網(wǎng)格單元進(jìn)行數(shù)值積分，流動(dòng)變量U在體單元中心的體積平均為：

(7)

式中:Vi,j,k為計(jì)算體網(wǎng)格單元，采用二階迎風(fēng)TVD格式對(duì)對(duì)流通量進(jìn)行離散。

(8)

其中，數(shù)值流通量分別為：

(9)

式中:λk，μk，νk為線性化替換矩陣的特征值;αk,βk,γk為線性化替換矩陣的展開項(xiàng)系數(shù);ek(A),ek(B),ek(C)為線性化替換矩陣的特征向量，黏性通量采用二階中心差分離散。

2.2.3湍流模型

湍流模型采用Realizablek-ε兩方程模型，與標(biāo)準(zhǔn)k-ε兩方程模型相比較，湍動(dòng)耗散率ε方程中的產(chǎn)生項(xiàng)不再包含有湍動(dòng)能k方程中的產(chǎn)生項(xiàng)Gk，并且湍動(dòng)黏度μt中的系數(shù)Cμ不是常數(shù)，而是與應(yīng)變速率相關(guān)，這樣的形式能更好地表示能量轉(zhuǎn)換。

湍動(dòng)能k方程為：

(10)

式中:μ為混合物黏性，Gk為平均速度梯度引起的湍動(dòng)能k的產(chǎn)生項(xiàng)，常數(shù)系數(shù)σk=1.0。

湍動(dòng)耗散率ε方程為：

(11)

式中:常數(shù)系數(shù)σε=1.2；C1=1.44，C2=1.9，E為平均特征應(yīng)變率。

2.2.4輻射模型

對(duì)輻射熱流采用離散坐標(biāo)法(DOM)進(jìn)行求解，該方法具有易于處理散射問題，易于和流動(dòng)方程聯(lián)立求解以及計(jì)算精度較高的特點(diǎn)。該方法最早由Chandrasekhar[14]提出，由輻射傳遞方程沿著s方向進(jìn)行離散得到：

在r位置沿著s方向的輻射傳遞方程為：

(12)

由式(12)可得光譜強(qiáng)度Iλ(r,s)的輻射傳遞方程為：

(13)

式中:r為位置矢量，s為方向矢量，s′為散射方向矢量，δ為氣體層厚度，a為吸收系數(shù)，n為折射率，σ為黑體輻射常數(shù)，σs為散射系數(shù)，I為光譜輻射強(qiáng)度，T為黑體的熱力學(xué)溫度，Φ為散射相函數(shù)，Ω′為立體角，λ為波長(zhǎng)，Ibλ為普朗克方程給定的黑體強(qiáng)度。

氣體的折射率n會(huì)受環(huán)境溫度Tδ以及波長(zhǎng)λ影響，具體如下：

(14)

式中:T0為標(biāo)準(zhǔn)溫度，n0為氣體在標(biāo)準(zhǔn)溫度下的折射率，Tδ表示環(huán)境溫度，η為線脹系數(shù)，可由文獻(xiàn)[15]查得。

吸收系數(shù)a及散射系數(shù)σδ可表示為：

(15)

式中:kλ為光譜減弱系數(shù)，βλ為散射衰減系數(shù)，ρ為氣體密度。

對(duì)散射角為θ的散射相函數(shù)定義為Φ(θ)，離散坐標(biāo)法求解散射光強(qiáng)正是利用勒讓德多項(xiàng)式展開，對(duì)Φ(θ)做近似處理，展開如下：

(16)

式中:Φg是第g階勒讓德函數(shù)，cg是第g階展開項(xiàng)的系數(shù)，可表示如下：

(17)

由式(14)可知，要獲得精確的散射相函數(shù)需要無(wú)窮多階勒讓德多項(xiàng)式展開，在實(shí)際計(jì)算時(shí)，由于只在前向周圍很窄的區(qū)域內(nèi)存在顯著的前向峰，因此，將散射相函數(shù)近似表示為前向部分和截?cái)嘞嗪瘮?shù)之和，即：

(18)

式中:f為截?cái)嘁蜃?，代表前向散射所占比例，M為近似的階數(shù)，cg*為第g階近似展開項(xiàng)系數(shù)。

3 數(shù)值結(jié)果

液體運(yùn)載火箭從地面發(fā)射后，高度不斷提升，隨海拔變化的主要參數(shù)包括環(huán)境壓強(qiáng)、環(huán)境溫度以及飛行馬赫數(shù)。

為研究液體運(yùn)載火箭底部熱環(huán)境變化趨勢(shì)，本文選取8個(gè)典型高度的工況進(jìn)行數(shù)值模擬，表1為這8個(gè)工況的自由流條件參數(shù)，噴管喉部總溫為3800 K，總壓為18.66 MPa，氣體組分包括H2O,CO2,CO,H2以及空氣等。

3.1 不同高度下的馬赫數(shù)場(chǎng)分布

由圖3可知，在5 km～25 km高度，噴管尾焰中馬赫波系結(jié)構(gòu)較為明顯，特別是在5 km和10 km 時(shí)，燃?xì)饬鲌?chǎng)中可見多個(gè)噴流波節(jié)。由于發(fā)動(dòng)機(jī)總壓遠(yuǎn)大于外部環(huán)境壓強(qiáng)，因此燃?xì)庠趪姽艹隹诟浇纬赏靶螇嚎s波，在馬赫數(shù)云圖中可觀察到清晰的尾焰噴流核心區(qū)、發(fā)展段以及射流邊界。

表1 自由流條件參數(shù)Table 1 Parametric free stream conditions

圖3 不同高度下的馬赫數(shù)云圖Fig.3 Contour of Mach number at different altitudes

隨著高度逐漸爬升，大氣環(huán)境壓強(qiáng)降低，噴流對(duì)外界環(huán)境的擾動(dòng)區(qū)域變大，噴流馬赫數(shù)峰值增大，噴流波節(jié)數(shù)減少，桶形激波與噴管軸線的交點(diǎn)坐標(biāo)向下游移動(dòng)。在35 km～40 km，噴流膨脹半徑及膨脹距離劇增，桶形激波與噴管軸線不相交。

燃?xì)庠诶郀枃姽苤刑幱诔R界工作狀態(tài)，在擴(kuò)張段形成超聲速膨脹波并進(jìn)一步加速，該膨脹波在噴管尾部由于受到大氣壓縮作用呈桶形，在噴流激波和噴流邊界下游與噴管軸線相交。隨著高度上升，環(huán)境壓強(qiáng)和溫度呈下降趨勢(shì)，外部伴隨氣流對(duì)噴流的壓縮作用逐漸減弱，噴流的膨脹程度繼續(xù)加強(qiáng)，這與圖3情況一致。

3.2 不同高度下的底部溫度場(chǎng)分布

選取4個(gè)典型工況的底部溫度場(chǎng)進(jìn)行對(duì)比展示，圖4為10 km,20 km,30 km和40 km這4個(gè)工況的箭體底部及噴管附近溫度場(chǎng)等值線分布圖。

圖4 不同高度下的溫度等值線圖Fig.4 Contour of temperature at different altitudes

從圖4中火箭底部溫度等值線的顏色及疏密程度可以看出，隨著飛行高度提升，火箭底部溫度整體為逐漸上升的趨勢(shì)。隨著海拔增加，噴管尾部噴流膨脹加強(qiáng)，膨脹后的噴流邊界與超聲速來流相互作用形成剪切層，噴管噴出的高溫燃?xì)庠诩羟袑恿鲃?dòng)作用下回流至火箭底部，形成對(duì)流熱流，自由來流、發(fā)動(dòng)機(jī)噴流激波與箭體底部及噴管壁面構(gòu)成一個(gè)封閉區(qū)域，對(duì)流熱流在封閉的剪切層區(qū)域內(nèi)形成再循環(huán)系統(tǒng)，對(duì)箭體底部持續(xù)進(jìn)行對(duì)流傳熱，造成底部對(duì)流熱流值在前期隨著海拔上升逐漸升高。但隨著海拔繼續(xù)上升，尾部噴流激波膨脹劇烈，致使所構(gòu)成的封閉區(qū)域減小，回流至箭體底部的高溫燃?xì)饬繙p小，這也是后期箭體底部對(duì)流熱流下降的原因。

3.3 不同高度下的底部對(duì)流熱流分布

選取4個(gè)典型工況的底部熱流云圖進(jìn)行對(duì)比展示，從圖5可以看出，20 km處對(duì)流熱流大于10 km處及30 km處，而40 km處對(duì)流熱流值最小。

圖5 不同高度下的對(duì)流熱流云圖Fig.5 Contour of convection heatflux at different altitudes

隨著海拔增加，環(huán)境壓力下降，噴流膨脹加強(qiáng)，引起更多回流，火箭底部的對(duì)流換熱增加，另一方面，海拔升高時(shí)，前方來流馬赫數(shù)(飛行速度)增加，對(duì)箭體底部回流的阻滯效果加強(qiáng)，因而減弱了火箭底部的對(duì)流加熱。在低海拔時(shí)，環(huán)境壓強(qiáng)對(duì)對(duì)流熱流的影響大于來流速度，所以隨著海拔提升對(duì)流熱流表現(xiàn)為增加的趨勢(shì)；在高海拔地區(qū)，環(huán)境壓強(qiáng)對(duì)對(duì)流熱流的影響小于來流速度，所以隨著海拔提升對(duì)流熱流表現(xiàn)為降低的趨勢(shì)。

4 數(shù)值模擬與飛行試驗(yàn)對(duì)比

在對(duì)某三級(jí)液體運(yùn)載火箭數(shù)值模擬過程中，選擇了與熱流傳感器相同的位置布置兩個(gè)監(jiān)測(cè)點(diǎn)，通過數(shù)值計(jì)算得到了監(jiān)測(cè)點(diǎn)的熱流參數(shù)，圖6顯示了兩個(gè)監(jiān)測(cè)點(diǎn)的數(shù)值計(jì)算總熱流值與飛行實(shí)際測(cè)量總熱流值變化曲線，可以看出兩者趨勢(shì)以及數(shù)值大小吻合度都較高?？偀崃鞒霈F(xiàn)了兩次峰值，分別是72 s 和99 s，對(duì)應(yīng)的正是10 km附近和20 km附近。

圖6中顯示的計(jì)算總熱流值可分解為圖7中的對(duì)流熱流和輻射熱流。圖7表明,72 s附近的總熱流峰值主要來源于輻射加熱，99 s附近的總熱流峰值主要來源于對(duì)流加熱。對(duì)流熱流在50～99 s逐漸增加，并在99 s達(dá)到最大值，然后在99～133 s逐漸減小，輻射熱流在剛開始達(dá)到最大值，然后開始降低。

圖6 總熱流對(duì)比Fig.6 Comparison of total heating rate

圖7 數(shù)值計(jì)算的對(duì)流熱流以及輻射熱流Fig.7 Computed convection and radiation heating rate

從圖7還能看出，在較低飛行高度時(shí)，輻射熱流占總熱流的比重較大，而在更高的海拔時(shí)，輻射熱流逐漸減小，對(duì)流熱流在火箭底部熱環(huán)境中扮演更加重要的角色。因?yàn)橐后w運(yùn)載火箭輻射傳熱介質(zhì)主要為H2O和CO2，隨著飛行高度上升，環(huán)境含氣量降低，只能由燃?xì)馀欧诺腍2O和CO2作為介質(zhì)進(jìn)行輻射傳熱，所以輻射熱流逐漸減小，同時(shí)，隨著海拔升高，火箭燃?xì)鈬娏鬓D(zhuǎn)變?yōu)楦叨扰蛎?，來自噴管附近的高溫氣體回流至火箭底部，對(duì)底部產(chǎn)生對(duì)流加熱，因而對(duì)流熱流增加，當(dāng)海拔繼續(xù)爬升時(shí)，外界來流馬赫數(shù)增加，來流對(duì)回流的阻滯效果加強(qiáng)，削弱了火箭底部的對(duì)流換熱，致使對(duì)流熱流下降。

5 結(jié) 論

本文對(duì)液體運(yùn)載火箭的高空底部熱流問題進(jìn)行了數(shù)值模擬，并與實(shí)際飛行數(shù)據(jù)進(jìn)行了對(duì)比分析，主要得出以下結(jié)論：

1)基于合理假設(shè)，建立含燃?xì)?空氣多組分Navier-Stokes方程、Realizablek-ε兩方程模型和熱輻射模型的火箭燃?xì)鈬娏髂Ｐ停捎酶咝У亩ATVD格式離散對(duì)流通量，離散坐標(biāo)法求解熱輻射方程，中心差分求解黏性引起的耗散項(xiàng)。對(duì)液體運(yùn)載火箭高空底部熱流進(jìn)行數(shù)值研究，仿真結(jié)果與測(cè)試數(shù)據(jù)吻合良好，校驗(yàn)了本文所采用的數(shù)值方法的精度和有效性。

2)本文的三級(jí)液體運(yùn)載火箭底部總熱流出現(xiàn)了兩次峰值，第一次峰值在10 km附近，主要受到輻射熱流影響所致；第二次峰值在20 km附近，主要由對(duì)流熱流所致。

3)分別對(duì)8個(gè)典型飛行高度的工況的熱流變化規(guī)律進(jìn)行研究，火箭底部的輻射熱流在剛起飛階段達(dá)到最大值，然后隨著海拔上升逐漸降低，對(duì)流熱流表現(xiàn)為先上升后下降的趨勢(shì)，并在20 km附近達(dá)到最大值。

4)在分析液體運(yùn)載火箭底部的熱環(huán)境時(shí)，需要考慮到對(duì)流/輻射耦合換熱的影響，本文所形成的火箭高空噴流模型和熱流數(shù)值模擬方法對(duì)液體運(yùn)載火箭底部熱防護(hù)安全性設(shè)計(jì)具有一定的參考價(jià)值。