孫培杰，嚴(yán)立，李雙菲，楊眉，陸辰昱，樂貴高，張亮，張衛(wèi)東

（1.上海宇航系統(tǒng)工程研究所，上海 201109；2.南京理工大學(xué)機(jī)械工程學(xué)院，江蘇 南京 210094；3.上海航天技術(shù)研究院，上海 201109）

0 引言

火箭發(fā)動機(jī)燃燒室通常以高壓、高溫的能量密度模式工作。由于燃?xì)饪倝号c外部大氣環(huán)境壓強(qiáng)差梯度很大，因此在火箭起飛階段，發(fā)射臺附近區(qū)域形成很強(qiáng)的斜激波、馬赫盤正激波、膨脹波，以及接觸面間斷和激波附面層干擾，對于固液混合運載火箭，還形成了燃?xì)夂虯l2O3顆粒兩相射流，引起流線或流體運動界面不穩(wěn)定以及引起化學(xué)反應(yīng)等復(fù)雜現(xiàn)象，這給起飛過程噴流研究工作帶來很大困難。20 世紀(jì)60 年代開始，國內(nèi)外通過大量的工程計算［1-3］、數(shù)值分析［4-5］、地面試驗［6］等對底部噴流和熱環(huán)境的規(guī)律、機(jī)理進(jìn)行了深入研究。國內(nèi)也采用數(shù)值計算方法對新一代運載火箭底部熱環(huán)境進(jìn)行預(yù)示分析，指導(dǎo)總體設(shè)計［7-10］。目前，國內(nèi)外對發(fā)射臺的熱環(huán)境研究報導(dǎo)甚少，公開文獻(xiàn)主要趨向于外形相對簡單的平板、導(dǎo)流器或?qū)Я鞑坌泵娴娜細(xì)饬鳑_擊效應(yīng)原理性研究［11］，具有大型運載火箭發(fā)射臺外形的研究僅局限于射流之間干擾分布狀況，而對固液捆綁火箭發(fā)射臺的力熱環(huán)境計算無文獻(xiàn)參考。

DETTLEFF等［13］設(shè)計了一套平板熱流計測量發(fā)動機(jī)尾部羽流場的熱流密度，通過熱電偶測量被沖擊平板一點的溫度變化來計算得到熱流密度。TSUTSUMI等［14］對火箭發(fā)動機(jī)布局對發(fā)射臺內(nèi)部流場的影響進(jìn)行了數(shù)值分析，選擇了具有典型發(fā)動機(jī)和發(fā)射臺的日本H3 運載火箭進(jìn)行分析。當(dāng)芯級有3 個液體發(fā)動機(jī)時，相鄰噴管的射流沖擊導(dǎo)流槽斜面后相互作用，會產(chǎn)生沿斜面向上的側(cè)向流動，而在側(cè)向流動區(qū)域布置一個發(fā)動機(jī)會減小噴管之間的相互作用。ALLGOOD等［15］對ARES V 概念級運載火箭進(jìn)行了數(shù)值計算，主要研究了發(fā)動機(jī)羽流對導(dǎo)流槽的沖擊特性，定量分析了導(dǎo)流槽斜面的壓力以及溫度分布情況。KIRIS等［16］使用美國宇航局Ames研究中心開發(fā)的升空與飛行器空氣動力學(xué)（lava）求解器，計算了導(dǎo)流槽斜面的瞬時熱流密度，并建立了塔架模型，研究了發(fā)動機(jī)羽流對塔架的沖擊狀況。趙業(yè)輝等［17］利用一套雙熱電偶測溫裝置和熱流測量計的組合測量裝置，測量了發(fā)動機(jī)噴管外高溫燃燒產(chǎn)物的溫度和熱流密度，并通過數(shù)值模擬的方法對固體火箭發(fā)動機(jī)噴管的外流場進(jìn)行了對比驗證，仿真結(jié)果與測量結(jié)果相似。陳勁松等［18］利用發(fā)射燃?xì)饬鲃恿W(xué)瞬態(tài)數(shù)值模擬方法，研究了雙面排導(dǎo)型發(fā)射臺承受的燃?xì)饬鳑_擊作用力隨時間變化特性。吳地勇等［19］對火箭發(fā)射時的燃?xì)饬鲌鲞M(jìn)行了仿真分析，監(jiān)測了發(fā)射臺上表面所受的最大壓力、最大總壓和最大靜溫，提取了相應(yīng)的發(fā)射臺表面壓力和靜溫分布矢量。周國儀等［20］采用計算流體力學(xué)方法，數(shù)值模擬了火箭發(fā)射時與帶孔平板形成的沖擊流場，分析了流場結(jié)構(gòu)及特性。盛文成等［21］使用有限體積法求解三維N-S 方程，通過數(shù)值模擬揭示導(dǎo)彈發(fā)射過程中導(dǎo)流器表面的壓強(qiáng)、溫度、馬赫數(shù)分布，以及燃?xì)饬鲗?dǎo)流器沖擊載荷變化規(guī)律。

以上研究均為火箭熱環(huán)境的評估提供了良好的參考依據(jù)。然而，對于固液捆綁火箭，由于固液組合推進(jìn)系統(tǒng)的燃?xì)鈱儆诙嘟M分、氣-固兩相超聲速流，其對發(fā)射臺和導(dǎo)流裝置的沖擊流動建模復(fù)雜，求解難度大，目前未見相關(guān)的研究報導(dǎo)。為此，本文擬建立燃?xì)?顆粒介質(zhì)的歐拉離散相可壓縮N-S 方程及熱輻射方程，并結(jié)合捕捉激波魯棒性強(qiáng)的Roe 迎風(fēng)差分格式，求解熱輻射方程效率高的離散坐標(biāo)法（Discrete-Ordinates Method，DOM），開展多簇噴流強(qiáng)沖擊發(fā)射臺和導(dǎo)流裝置的力熱效應(yīng)數(shù)值模擬。

1 控制方程與計算方法

1.1 控制方程

假設(shè)液體燃料火箭燃?xì)舛嘟M分均滿足連續(xù)介質(zhì)假設(shè)和理想氣體狀態(tài)方程，建立三維多組分混合燃?xì)獾目蓧嚎sNavier-Stokes 方程組，具體形式如下：

式中：U為守恒變量；Q為源項；F、G、H分別為x、y、z方向上的無黏通量；Fv、Gv、Hv分別為x、y、z方向上的黏性通量，可以表示為

式中：ρ為流體密度；E為能量密度；p為流體所受的正應(yīng)力；τ為流體黏性力；K為熱傳導(dǎo)系數(shù)；T為流體溫度；u、v、w為速度在x、y、z3 個方向的分量。

為求解黏性力τij，可以建立應(yīng)力與形變速率之間的關(guān)系，根據(jù)本構(gòu)方程可以得到剪應(yīng)力和法向應(yīng)力的表達(dá)式如下：

式中：μ為流體黏性系數(shù)。

1.2 湍流模型

為了計算提高黏性引起的湍流耗散的計算精度，Launder和Spalding最早提出并經(jīng)Yakhot和Orzag修正的k-ε兩方程湍流模型即RNGk-ε模型［22-24］。

1.3 輻射模型

采用DOM，該方法易于處理散射問題，易于和流動方程聯(lián)立求解，而且計算精度較高。由輻射傳遞方程沿著s方向進(jìn)行離散得到：

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

可得光譜強(qiáng)度Iλ（r，s）的輻射傳遞方程為

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

1.4 離散相模型

歐拉離散項模型將顆粒和氣體處理為2 種相互貫穿的連續(xù)相，空間各點都有2 種連續(xù)相各自不同的速度、溫度和密度。兩相之間相互滲透，各自具有不同的體積分?jǐn)?shù)，相互作用并且相互之間有滑移。因此顆粒物在空間中有連續(xù)的速度、溫度及體積分?jǐn)?shù)的分布。歐拉離散相模型的控制方程如下

式中：ρp為離散相密度；up為離散相速度矢量為質(zhì)點內(nèi)能。相間拖拽力FD和傳熱項Qp定義如下

式中：v為流體的速度矢量；fD為斯托克斯修正系數(shù)；τu為斯托克斯松弛系數(shù)；Cp為離散相的定壓比熱；fN為努塞爾修正系數(shù)；τT為溫度松弛時間；TP為離散相溫度。

1.5 差分格式

由于火箭燃?xì)馍淞鲗儆诟咚倏蓧嚎s流，流場中有強(qiáng)射流激波和馬赫盤正激波的流動變量間斷現(xiàn)象，需要激波捕捉能力強(qiáng)的數(shù)值格式，這里選擇二階對流迎風(fēng)差分格式［25］，該方法經(jīng)過大量計算檢驗，具有魯棒性強(qiáng)，收斂于物理解數(shù)值求解。流動方程組中耗散項采用二階中心差分離散。

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

2.1 數(shù)值方法驗證

實驗裝置為某標(biāo)準(zhǔn)試驗發(fā)動機(jī)，采用常規(guī)含鋁丁羥（HTPB）推進(jìn)劑。實驗裝置見文獻(xiàn)［17］，在距離發(fā)動機(jī)出口3 m 的軸線上放置傳感器安裝平臺。傳感器主要包括熱電偶和熱流測量計，用以測量發(fā)動機(jī)尾焰的溫度和平臺表面熱流密度。射流穩(wěn)定狀態(tài)下，溫度約為500 ℃（773.15 K），熱流密度為2.320 MW/m2。對該標(biāo)準(zhǔn)試驗發(fā)動機(jī)尾焰進(jìn)行數(shù)值模擬。穩(wěn)定工作下兩相燃?xì)饬鞲鹘M分的質(zhì)量分?jǐn)?shù)見表1。

表1 兩相燃?xì)饬鞲鹘M分質(zhì)量分?jǐn)?shù)Tab.1 Mass fractions of each component

傳感器平臺附近燃?xì)饬骱蜔崃鳒囟惹闆r如圖1和圖2 所示。從圖1 和圖2 可以看出，數(shù)值模擬得到的傳感器平臺附近燃?xì)饬鳒囟燃s為790 K，傳感器平臺表面的熱流密度為2.234 MW/m2。而測試的燃?xì)饬鳒囟燃s為773.15 K，熱流密度為2.32 MW/m2，兩者溫度相對誤差為2.2%，熱流密度相對誤差為3.7%。計算與實驗吻合良好，驗證了計算模型的精確性。

圖1 傳感器平臺附近燃?xì)饬鳒囟惹闆rFig.1 Temperature of the gas flow near the sensor platform

圖2 傳感器平臺表面熱流密度Fig.2 Surface heat flux of the sensor platform

2.2 火箭起飛階段流場與熱環(huán)境數(shù)值分析

2.2.1 計算模型

火箭地面發(fā)射系統(tǒng)包括發(fā)射臺、導(dǎo)流孔、導(dǎo)流槽、支撐臂、服務(wù)塔、箭體等，其中運載火箭由芯級和4 個助推器捆綁而成，芯級為液體燃料，裝有2 臺發(fā)動機(jī)，助推器為固體燃料，每個助推器安裝1 臺發(fā)動機(jī)?；鸺l(fā)射臺幾何模型如圖3 所示。

圖3 火箭發(fā)射臺幾何模型Fig.3 Geometric model of the rocket launch pad

2.2.2 計算網(wǎng)格

由于火箭發(fā)射裝置幾何外形十分復(fù)雜，因此采用分塊化網(wǎng)格生成方法將整個計算域分成多個子塊，將各子塊分別生成六面體結(jié)構(gòu)網(wǎng)格，以保證網(wǎng)格的正交性及光順性。發(fā)射臺、支撐臂、服務(wù)塔、導(dǎo)流孔和噴管局部的網(wǎng)格如圖4 所示。

圖4 火箭發(fā)射臺各子塊表面計算網(wǎng)格Fig.4 Iso-surface meshes of the rocket launch pad

2.2.3 邊界條件

①噴管入口條件：采用總壓總溫邊界條件，芯級發(fā)動機(jī)入口壓力17.7 MPa，溫度3 805.34 K；助推發(fā)動機(jī)入口壓力6.3 MPa，溫度3 235.88 K；在助推器噴管出口燃?xì)鉃?74 kg 的情況下，出口凝相質(zhì)量Al2O3為152.248 8 kg。②外部邊界定義：大氣環(huán)境值。③壁面條件：噴管內(nèi)壁采用絕熱壁面條件，而發(fā)射臺、導(dǎo)流孔、發(fā)射臂、箭體底部和噴管外壁面定義為恒溫壁面邊界條件。

2.2.4 計算狀態(tài)

分別以0.2、3.0、20.0、50.0 m 4 種典型高度，X方向（靠近塔架方向，即平行于服務(wù)塔方向）、Y方向（平行于塔架方向指向?qū)Я鞑?，即垂直于服?wù)塔方向）和45°方向（支撐臂方向）3 種漂移方向的飛行狀態(tài)為對象，研究火箭點火和起飛時燃?xì)鈬娏魉矐B(tài)過程及燃?xì)饬鳡顟B(tài)。漂移參數(shù)見表2。

表2 各高度漂移參數(shù)Tab.2 Drift parameters at each height

2.3 火箭起飛階段分析

過芯級噴管中心截面的馬赫數(shù)分布如圖5 所示，由圖可見隨著箭體上升和Y方向漂移，芯級噴管右側(cè)燃?xì)庥鹆髦饾u沖擊在發(fā)射臺面，羽流反濺程度逐漸加強(qiáng)。過助推器其中2 個噴管中心截面的馬赫數(shù)分布如圖6 所示，由圖可見隨著火箭上升和Y方向漂移，靠近右側(cè)的羽流越來越多地沖擊在發(fā)射臺面上并發(fā)生反濺?；鸺仙⒀刂鳼方向或者X方向漂移時，助推器排出的Al2O3顆粒質(zhì)量濃度分布如圖7 所示，由圖可見當(dāng)羽流遇到發(fā)射臺后，同樣表現(xiàn)一定的沖擊和反濺?；鸺痫w并沿著Y方向漂移時燃?xì)饬鳑_擊在發(fā)射臺面的壓強(qiáng)分布如圖8 所示，由圖可見隨著漂移量增加，臺面的沖擊范圍增大，當(dāng)上升20.0 m 時，發(fā)射臺面壓強(qiáng)最大，達(dá)到0.885 MPa，隨后最大值減小。

圖5 芯級中心截面的馬赫數(shù)分布Fig.5 Mach number distribution through the core state central section

圖6 過助推器截面的馬赫數(shù)分布Fig.6 Mach number distribution through the boosting rocket central section

圖7 過助推器截面的Al2O3顆粒質(zhì)量濃度分布Fig.7 Concentration distribution of Al2O3 particles through the core state central section

圖8 發(fā)射臺表面壓力分布Fig.8 Pressure distribution on the rocket launch pad surface

支撐壁表面的總熱流密度分布如圖9 所示，由圖可見當(dāng)火箭升至20.0 m 時，靠近Y方向側(cè)壁面的總熱流密度更大，峰值達(dá)到8.63 MW/m2，隨著火箭高度增加，Y軸方向一側(cè)的發(fā)射臂面熱流密度較高的區(qū)域增大，但最大總熱流密度逐漸下降。火箭底部熱流隨著高度的增加，熱流逐漸減低，在起飛時刻熱流最大，為1.50 MW/m2，如圖10 所示。火箭起飛階段沿Y軸漂移時噴管外壁面的總熱流密度分布如圖11 所示，由圖可見在3.0 m 高度噴管總熱流密度最大，峰值為1.27 MW/m2，因為初始階段燃?xì)饬鳑_擊導(dǎo)流槽和發(fā)射臺發(fā)生一定的反濺，隨后逐漸減小。

圖9 支撐臂表面總熱流密度Q 分布Fig.9 Total heat flux distribution of the support arm surface

圖10 箭體底板總熱流密度Q 分布Fig.10 Total heat flux distribution of the rocket base plate

圖11 噴管外壁面總熱流密度Q 分布Fig.11 Total heat flux distribution on the nozzle outer wall

芯級和助推發(fā)動機(jī)燃?xì)饬鲗Πl(fā)射臺和噴管外壁面產(chǎn)生的最大總熱流密度、最大總壓力和熱流密度集中區(qū)域的燃?xì)饬鳒囟冉y(tǒng)計見表3～表5。

表3 助推發(fā)動機(jī)最大總熱流密度Q 匯總表Tab.3 Maximum total heat flux of the booster engine單位：MW/m2

表4 助推發(fā)動機(jī)最大總壓力匯總表Tab.4 Maximum total pressure of the booster engine單位：Pa

表5 助推發(fā)動機(jī)熱流密度集中區(qū)域附近燃?xì)饬鳒囟葏R總表Tab.5 Temperature of the combustion gas flow near the heat flux concentration area of the booster engine單位：K

助推噴管在發(fā)射臺表面形成的熱流密度集中區(qū)域最大總熱流密度值雖然僅為芯級的1/5，但其熱流密度較高的區(qū)域主要分布在支撐臂、服務(wù)塔等設(shè)備附近。因此，進(jìn)行熱防護(hù)需要著重考慮固體發(fā)動機(jī)燃?xì)獾挠绊憽?/p>

3 結(jié)束語

在火箭點火和起飛階段力熱環(huán)境大型數(shù)值模擬過程中，主要取得了以下結(jié)論：

1）基于連續(xù)介質(zhì)假設(shè)，建立適用于含燃?xì)?顆粒多組分的兩相射流模型及其數(shù)值方法，基于黑體輻射理論建立了燃?xì)饬鲌鰺彷椛浞匠蹋x取標(biāo)準(zhǔn)火箭噴管為驗證對象，開展了固體火箭燃?xì)馍淞鳑_擊流場的數(shù)值仿真，得到溫度和熱流傳感器平臺數(shù)值結(jié)果與實驗數(shù)據(jù)進(jìn)行對比，溫度、熱流密度兩者的相對誤差分別為2.2%和3.7%，驗證了本文數(shù)值方法的精度和計算手段的有效性。

2）在火箭起飛初始階段，由于燃?xì)饬鞒跏紱_擊波傳播范圍逐漸變大，以及箭體漂移，有部分燃?xì)馍淞髋c發(fā)射臺面碰撞、反濺，對箭體底部有一些影響，隨著箭體進(jìn)一步上升，反濺逐漸減弱，演變?yōu)槎啻貒娏飨嗷ジ蓴_的形態(tài)影響箭體底部，這種影響比燃?xì)饬鞣礊R的影響要小。

3）從馬赫數(shù)分布云圖看，芯級射流半徑要大于助推器射流，這是因為：其一，助推器為固體燃料，射流中含有一定濃度的Al2O3，不利于射流的快速膨脹；其二，助推器噴管總壓、總溫均比芯級噴管的低。因此，燃?xì)馍淞髑放蛎洺潭炔煌炯壍母鼜?qiáng)。

4）通過大量數(shù)值模擬，得到火箭起飛和漂移狀態(tài)下固液捆綁火箭多簇噴管氣固兩相燃?xì)馍淞鳑_擊力熱效應(yīng)，包括發(fā)射臺面、支撐臂、箭體底部和噴管外壁面產(chǎn)生的最大總熱流密度、最大總壓力和熱流密度集中區(qū)域的燃?xì)饬鳒囟确逯到y(tǒng)計，其結(jié)果可為發(fā)射臺和箭體底部等重點部位的防熱設(shè)計提供指導(dǎo)。