李國(guó)良，楊云軍，龔安龍，劉 周

(中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京100074)

0 引 言

導(dǎo)彈飛行過(guò)程中發(fā)動(dòng)機(jī)噴流與來(lái)流繞流的作用會(huì)改變彈體的氣動(dòng)特性，尤其是對(duì)底阻特性的影響較大。底阻受到彈體長(zhǎng)度、邊界層狀態(tài)、尾部形狀、發(fā)動(dòng)機(jī)噴流參數(shù)、飛行高度和馬赫數(shù)等因素的影響,在高空高馬赫數(shù)狀態(tài)下底部有時(shí)會(huì)出現(xiàn)正推力。探索飛行器底部復(fù)雜流動(dòng)結(jié)構(gòu)及規(guī)律，是近年來(lái)研究的熱點(diǎn)也是難點(diǎn)之一。在不考慮發(fā)動(dòng)機(jī)噴流狀態(tài)下，超聲速底部流動(dòng)包含膨脹波、自由剪切層、分離渦、激波等多種流動(dòng)結(jié)構(gòu)，分離區(qū)中的大渦結(jié)構(gòu)與激波的相互干擾導(dǎo)致流動(dòng)呈現(xiàn)高頻非定常性。在考慮底部發(fā)動(dòng)機(jī)噴流狀態(tài)下，底部高溫高速噴流與來(lái)流相互干擾，二次回流區(qū)在噴口附近形成了環(huán)帶，因此與無(wú)噴流狀態(tài)相比，底部發(fā)動(dòng)機(jī)噴流改變了底部流場(chǎng)結(jié)構(gòu)及飛行器的氣動(dòng)特性。

底部流動(dòng)研究涉及兩大類問(wèn)題：無(wú)噴流時(shí)的底部大分離流動(dòng)和發(fā)動(dòng)機(jī)噴流干擾下的摻混流動(dòng)。目前國(guó)內(nèi)外研究學(xué)者對(duì)無(wú)噴流狀態(tài)底部流動(dòng)開展了大量研究。肖志祥等[1]采用RANS/LES混合方法研究超聲速底部流動(dòng)，對(duì)包含豐富流動(dòng)結(jié)構(gòu)和復(fù)雜流動(dòng)機(jī)理的導(dǎo)彈類超聲速底部流場(chǎng)進(jìn)行數(shù)值分析。高瑞澤等[2]采用LES/RANS混合方法及Smagorinsky-BL及Vreman-BL模型對(duì)超聲速底部流動(dòng)開展了計(jì)算，獲得了較好的速度型及底面壓力分布結(jié)果。朱德華等[3]采用二階總變差不增格式(TVD)及二階Runge-Kutta法模擬了高超聲速鈍錐底部尾跡流場(chǎng)，觀察到了底部出現(xiàn)了非定常周期性流動(dòng)，呈現(xiàn)流動(dòng)的不穩(wěn)定性。國(guó)外的研究學(xué)者Soshi等[4]、Pramod等[5]、Franck等[6]、Krishnendu[7]、Soo等[8]對(duì)底部大分離流動(dòng)開展了計(jì)算方法的適用性研究，并與試驗(yàn)進(jìn)行了比較。

[4-8]的計(jì)算結(jié)果表明，RANS方法不能完全準(zhǔn)確模擬底部流動(dòng)大分離特性。但本文的計(jì)算模型底部均帶有噴管，噴管的存在實(shí)際削弱了底部流動(dòng)的大分離特性，呈現(xiàn)了定常的特性。根據(jù)參考文獻(xiàn)[9-11]，RANS方法適用于此類問(wèn)題的數(shù)值模擬。Bannink等[9]、Bakker等[10]采用了冷噴方法開展鈍錐旋成體底部流場(chǎng)計(jì)算，計(jì)算方法為RANS結(jié)合不同的湍流模型，包括S-A、k-ω、剪切壓力傳輸湍流模型(SST)等，通過(guò)與試驗(yàn)結(jié)果的比較得出采用SST湍流模型獲得的結(jié)果與試驗(yàn)一致性最好。林敬周等[11]采用二階無(wú)波動(dòng)無(wú)自由參數(shù)的耗散差分(NND)格式，B-L代數(shù)模型對(duì)文獻(xiàn)[9-10]中的外形開展了底部噴流對(duì)來(lái)流干擾的研究，結(jié)果表明有無(wú)噴流對(duì)底部流場(chǎng)結(jié)構(gòu)影響很大，同時(shí)噴流對(duì)底部壓力分布有顯著影響。田耀四等[12]采用單組分變比熱比二維軸對(duì)稱方程對(duì)發(fā)動(dòng)機(jī)噴流開展了數(shù)值研究，對(duì)不同飛行高度及飛行速度的尾噴流激波流場(chǎng)特性進(jìn)行了詳細(xì)分析。

依據(jù)文獻(xiàn)調(diào)研發(fā)現(xiàn)，國(guó)內(nèi)外在考慮底部發(fā)動(dòng)機(jī)噴流干擾下的摻混流動(dòng)數(shù)值模擬研究較少，而且多是采用冷噴方法，即認(rèn)為來(lái)流與噴流為同一介質(zhì)，通常為空氣。本文建立兩組分熱噴方法，即來(lái)流為空氣，發(fā)動(dòng)機(jī)燃?xì)鉃榱硪环N介質(zhì)。這種計(jì)算方案能夠更真實(shí)地反映燃?xì)獾呐蛎浶?yīng)，獲得可信度更高的底部壓力分布，更好地應(yīng)用到工程問(wèn)題中。

本文采用RANS結(jié)合SST湍流模型，并與試驗(yàn)數(shù)據(jù)進(jìn)行比較，驗(yàn)證方法的準(zhǔn)確性。同時(shí)針對(duì)本文選用的飛行器外形比較冷噴與熱噴兩種處理方法對(duì)計(jì)算結(jié)果的差異。最后采用熱噴方法對(duì)其超聲速底部流動(dòng)開展計(jì)算研究，并與監(jiān)測(cè)點(diǎn)試驗(yàn)數(shù)據(jù)進(jìn)行比較。選定來(lái)流馬赫數(shù)2.5，開展不同飛行高度及發(fā)動(dòng)機(jī)總壓對(duì)燃?xì)鈹U(kuò)散及底部流場(chǎng)的影響研究。

1 計(jì)算方法與算例驗(yàn)證

1.1 方法介紹

湍流采用SST兩方程模型，Roe的Riemann近似解算器計(jì)算無(wú)黏通量。時(shí)間推進(jìn)采用上下對(duì)稱的高斯-賽德爾(LU-SGS)方法。熱噴計(jì)算時(shí)采用雙組分的三維Navier-Stokes方程，如下所示：

(1)

(2)

(3)

其中，

(4)

式中：V為逆變速度，u，v，w為三個(gè)速度分量，ni為單位矢量在坐標(biāo)系下的分量;p為壓力，E為總能，ρ為密度，H為總焓;yi為第i種氣體的質(zhì)量分?jǐn)?shù),Hi為第i種氣體的總焓，Di為第i種氣體的擴(kuò)散系數(shù);ρD=μl/Scl+μt/Sct,μl為層流黏性系數(shù)，μt為湍流黏性系數(shù)，Scl為層流施密特?cái)?shù)，Sct為湍流施密特?cái)?shù);τij為應(yīng)力張量項(xiàng)，qi為熱傳導(dǎo)項(xiàng)。

在考慮冷噴計(jì)算時(shí)，噴流與來(lái)流為同一種介質(zhì)，即空氣，熱噴控制方程褪化為常規(guī)Navier-Stokes方程。

1.2 噴流算例驗(yàn)證

選取經(jīng)典試驗(yàn)開展方法驗(yàn)證，試驗(yàn)?zāi)Ｐ蜑檩S對(duì)稱鈍錐外形[9-10]，底部噴管出口與喉部面積比為10.7584，進(jìn)口與喉部面積比為9,如圖1所示。

來(lái)流馬赫數(shù)Ma∞=1.96，靜壓p∞=28300 Pa，攻角α=0。圖2所示為無(wú)噴流時(shí)底部流線，在底部存在兩對(duì)分離渦。圖3所示為無(wú)噴流時(shí)模型柱體部分壓力系數(shù)計(jì)算與試驗(yàn)[9]的比較?？v坐標(biāo)壓力系數(shù)為測(cè)點(diǎn)壓力與來(lái)流壓力的比值(下同)。從圖3可以看出，計(jì)算值與試驗(yàn)值一致性較好，總體誤差控制在0.5%以內(nèi)。

在考慮噴流條件下，計(jì)算了壓力比(噴管進(jìn)口總壓pin,t與來(lái)流靜壓p∞之比)為12.9和44.0兩種情況，噴流介質(zhì)為高壓空氣。圖4為無(wú)噴流及兩種壓力比三種情況下的計(jì)算值與試驗(yàn)值的比較。橫坐標(biāo)表示為底部位置徑向坐標(biāo)與模型半徑的比值。在無(wú)噴流時(shí)，底部計(jì)算壓力系數(shù)與試驗(yàn)值的總體誤差控制在8%以內(nèi)，壓力比12.9時(shí)誤差控制在2.5%以內(nèi)，壓力比44.0時(shí)誤差控制在5.5%以內(nèi)。從數(shù)值比較來(lái)看，本文所建立的來(lái)流干擾條件下發(fā)動(dòng)機(jī)底部噴流數(shù)值方法是準(zhǔn)確的。同時(shí)計(jì)算與試驗(yàn)都表現(xiàn)出隨著壓力比的不斷增加，底部的壓力系數(shù)逐漸下降的趨勢(shì)。這表明隨著壓力比增加，噴流氣體對(duì)底部流場(chǎng)的引射作用增強(qiáng)。還可以看出，三種情況下，底部壓力系數(shù)不隨徑向位置的改變發(fā)生變化，基本保持恒定。圖5為壓力比44.0時(shí)底部流場(chǎng)流線，從圖5可以看出，高速噴流氣體削弱了底部分離區(qū)域的大小，同時(shí)改變了底部流場(chǎng)的渦結(jié)構(gòu)。

圖6為流場(chǎng)當(dāng)?shù)豴與p∞比值的等值線云圖，此時(shí)的壓力比增加到349，底部流動(dòng)出現(xiàn)了明顯的流動(dòng)特征，即出現(xiàn)了桶狀激波(圖中標(biāo)注為1)與羽流激波(圖中標(biāo)注為2)。這與參考文獻(xiàn)[9-10]的試驗(yàn)結(jié)果與計(jì)算結(jié)果一致。

2 飛行器外形與冷/熱噴計(jì)算方法對(duì)比

2.1 飛行器外形參數(shù)

計(jì)算模型為圖7所示的三維軸對(duì)稱尖錐外形飛行器，底部為發(fā)動(dòng)機(jī)噴管。模型總長(zhǎng)為10 m，柱體直徑為1.45 m,計(jì)算網(wǎng)格采用分區(qū)多塊結(jié)構(gòu)方式。監(jiān)測(cè)傳感器位置位于模型底部xy對(duì)稱面上，距離噴管中心線0.65 m。

2.2 冷噴與熱噴計(jì)算結(jié)果對(duì)比

對(duì)于發(fā)動(dòng)機(jī)噴流問(wèn)題處理方式分為等效冷噴與熱噴。熱噴方法中來(lái)流與發(fā)動(dòng)機(jī)燃?xì)饩鶠榱繜嵬耆珰怏w。冷噴等效通常有兩種方法：1)出口動(dòng)量、壓力比相等，馬赫數(shù)不同；2)出口動(dòng)量、馬赫數(shù)相等，壓力不同。

以計(jì)算來(lái)流條件：來(lái)流溫度T∞=227 K，Ma∞=1.8，P∞=29400 Pa，α=0，Pin,t=6.0×106Pa，燃燒室總溫Tin,t=3500 K，燃?xì)鈿怏w常數(shù)Rc=284 J/(kg·K)為例來(lái)說(shuō)明兩種冷噴等效方法的差異。

等熵流Ma和面積A之間的關(guān)系式：

(5)

喉道Mat=1，發(fā)動(dòng)機(jī)燃?xì)猞?1.1，出口面積Aout與喉部面積At比11.6，半徑比3.40。得到出口Maout=3.09。

出口靜壓Pout與Pin,t的關(guān)系式：

(6)

式中:Pout=81970 Pa。

出口靜溫Tout與Tin,t的關(guān)系式：

(7)

式中:Tout=2369 K。

1)動(dòng)量相等，壓力比相等，馬赫數(shù)變化

2)動(dòng)量相等，馬赫數(shù)相等，壓力變化

出口馬赫數(shù)相等，等效出口MaA=Maout=3.09。出口動(dòng)量相等，γAPA=γPout，空氣γA=1.4，得到出口PA=64405 Pa，與實(shí)際Pout相差21.4%。得到空氣下的出口與喉部面積比4.61，半徑比2.15。在此條件下根據(jù)式(3)得Pin,t=2706116 Pa，總溫不變。

綜上所述，動(dòng)量、壓力比相等時(shí)，出口馬赫數(shù)與實(shí)際值差11.3%；動(dòng)量、馬赫數(shù)相等時(shí)，出口壓力與實(shí)際值差21.4%。因此選取第1種相似條件，即動(dòng)量、壓力比相等更為合適。在此相似條件下噴管進(jìn)口Main=0.49，靜壓Pin=1722166 Pa， 靜溫Tin=3340 K。

圖8為噴管原始形線、壓力比相等形線、馬赫數(shù)相等形線示意圖。圖9為采用冷噴與熱噴方法計(jì)算的噴管出口的馬赫數(shù)分布曲線，與理論計(jì)算的結(jié)果一致性較好，熱噴方法計(jì)算的噴管出口馬赫數(shù)平均值為3.10，與理論值的誤差為0.3%。冷噴方法計(jì)算的噴管出口馬赫數(shù)平均值為2.80，與理論值的誤差為2.2%。

采用冷噴方法與熱噴方法開展計(jì)算結(jié)果的對(duì)比分析。圖10為冷噴方法與熱噴方法的底部流線圖，從圖中可以看出，兩種方法流場(chǎng)結(jié)構(gòu)基本相同，在底部凹腔內(nèi)形成一對(duì)主分離渦及小分離渦。圖11為兩種方法計(jì)算值與試驗(yàn)值的定量比較。與第1.2節(jié)算例驗(yàn)證類似，底部壓力系數(shù)隨著y坐標(biāo)的變化基本保持恒定，表明底部的壓力相對(duì)均勻。監(jiān)測(cè)點(diǎn)的試驗(yàn)壓力系數(shù)-0.147，冷噴方法監(jiān)測(cè)點(diǎn)壓力系數(shù)-0.236，與試驗(yàn)誤差為60.5%。熱噴方法監(jiān)測(cè)點(diǎn)的壓力系數(shù)-0.162，與試驗(yàn)誤差10.2%。造成誤差差異的原因?yàn)橥ㄟ^(guò)熱噴方法計(jì)算保持了原來(lái)噴管的外形，噴管出口位置的傾斜角度為14.75°。冷噴方法改變了噴管外形，噴管出口位置的傾斜角度減小為9.27°。阻礙了噴管出口的高壓氣體向空腔內(nèi)擴(kuò)散，因此冷噴方法計(jì)算的壓力系數(shù)值較熱噴方法要低。

采用熱噴方法計(jì)算得到的壓力數(shù)據(jù)與試驗(yàn)值一致性較好，表明本文采用的兩組分熱噴RANS方法是可靠的。燃?xì)獾慕M分分布無(wú)法直接測(cè)量，但是燃?xì)庀虻撞繑U(kuò)散，直接影響底部壓力的分布。從壓力比較的結(jié)果可知本文采用的方法適用于燃?xì)饨M分分布的數(shù)值模擬。

圖12為同一算例采用熱噴方法計(jì)算得到的發(fā)動(dòng)機(jī)燃?xì)馕擦髻|(zhì)量分?jǐn)?shù)空間分布，燃?xì)馕擦鞒视馉罘植贾饾u向周圍擴(kuò)散。

3 飛行器底部流場(chǎng)分析

采用熱噴方法對(duì)飛行器開展Ma∞=2.5典型狀態(tài)不同飛行高度H及Pin,t對(duì)噴流及底部流場(chǎng)的影響研究。

3.1 進(jìn)口總壓Pin,t變化

H=15 km不變，改變Pin,t，5種條件分別為6.0 MPa，7.2 MPa，8.4 MPa，9.6 MPa，10.8MPa，研究發(fā)動(dòng)機(jī)燃?xì)鈪?shù)對(duì)底部擴(kuò)散及壓力分布的影響。

圖13為底部壓力系數(shù)隨Pin,t變化曲線，隨著Pin,t提高，底部壓力系數(shù)逐漸增加，同時(shí)5種條件下底部壓力分布均各自保持基本恒定。

圖14為底部燃?xì)赓|(zhì)量分?jǐn)?shù)隨Pin,t變化曲線，隨著Pin,t提高，壁面處的燃?xì)赓|(zhì)量分?jǐn)?shù)變化較小，在0.155～0.165之間。在底部凹腔x=9.9 m壁面處燃?xì)赓|(zhì)量分?jǐn)?shù)同樣基本保持不變，在0.163左右。其余位置燃?xì)赓|(zhì)量分?jǐn)?shù)逐漸提高。從圖14可以看出，發(fā)動(dòng)機(jī)燃?xì)庀虻撞堪记粌?nèi)擴(kuò)散，壁面處的燃?xì)鉂舛茸罡撸疫_(dá)到一個(gè)額定值后，隨Pin,t提高基本不變。

3.2 飛行高度H變化

Pin,t=7.2 MPa不變，改變H，6種條件分別為11 km，13 km，15 km，17 km，19 km，21 km。從圖15可以看出，Pin,t不變，底部壓力系數(shù)隨著H的增加單調(diào)增加。在H=21 km時(shí)底部壓力系數(shù)為0.0147，即底部出現(xiàn)正推力。正推力的出現(xiàn)有可能改變飛行器的飛行彈道，在設(shè)計(jì)階段需要提前開展評(píng)估。

圖16為Ma∞=2.5，H=11 km，13 km，15 km，17 km，19 km，21 km時(shí)噴管出口下游中心線馬赫數(shù)分布。圖中坐標(biāo)軸x代表中心線上一點(diǎn)與噴管出口中心點(diǎn)的距離，x=0代表噴管出口中心點(diǎn)處。從圖16可以看出,來(lái)流馬赫數(shù)一定，噴管出口中心線上發(fā)動(dòng)機(jī)燃?xì)庾畲篑R赫數(shù)隨著飛行高度增加不斷后移且數(shù)值不斷增大。

圖17、圖18為Ma∞=2.5，H=11 km，21 km時(shí)底部馬赫數(shù)云圖，從圖中可以看出,H=11 km時(shí)發(fā)動(dòng)機(jī)燃?xì)庾畲篑R赫數(shù)為4.7，H=21 km時(shí)發(fā)動(dòng)機(jī)燃?xì)庾畲篑R赫數(shù)為5.9，同時(shí)發(fā)動(dòng)機(jī)燃?xì)庠贖=21 km時(shí)膨脹效應(yīng)更強(qiáng)，高馬赫數(shù)的范圍更廣。

圖19為H=11 km在距噴管出口0.2 m,0.5 m,1.0 m,2.0 m處沿y方向的燃?xì)赓|(zhì)量分?jǐn)?shù)曲線。圖20為H=21 km在距噴管出口0.2 m,0.5 m,1.0 m,2.0 m處沿y方向的燃?xì)赓|(zhì)量分?jǐn)?shù)曲線。H=11 km，在x=1.0 m處燃?xì)獾呐蛎浄秶?.2 m，在x=2.0 m處燃?xì)獾呐蛎浄秶?.35 m。H=21 km，在x=1.0 m處燃?xì)獾呐蛎浄秶?.7 m，在x=2.0 m處燃?xì)獾呐蛎浄秶?.1 m。由于本文選用的飛行器長(zhǎng)度較長(zhǎng)，底部的環(huán)境壓力恢復(fù)為來(lái)流壓力，飛行高度增加，來(lái)流壓力減小，底部環(huán)境壓力對(duì)發(fā)動(dòng)機(jī)燃?xì)獾募s束減弱，發(fā)動(dòng)機(jī)燃?xì)獾呐蛎浶?yīng)增強(qiáng)，高速區(qū)整體向后移動(dòng)。

4 結(jié) 論

本文建立了考慮發(fā)動(dòng)機(jī)噴流的三維冷噴與兩組分熱噴數(shù)值模擬方法，在相同的來(lái)流與發(fā)動(dòng)機(jī)參數(shù)條件下，進(jìn)行了冷噴與熱噴方法的對(duì)比分析，兩組分熱噴方法獲得的數(shù)據(jù)與試驗(yàn)值一致性較好，較冷噴方法精度有明顯提升。采用熱噴方法，對(duì)飛行器開展了飛行高度及噴管進(jìn)口總壓對(duì)噴流及底部流場(chǎng)的影響研究。隨著噴管進(jìn)口總壓增加，底部壓力系數(shù)逐步提高。最高燃?xì)赓|(zhì)量濃度位于底部壁面處，而且達(dá)到一個(gè)額定值后，不隨總壓提高而變化。隨著飛行高度的增加，底部壓力系數(shù)單調(diào)增加，噴管出口中心線上發(fā)動(dòng)機(jī)燃?xì)庾畲篑R赫數(shù)不斷后移且數(shù)值不斷增加，在一定的飛行高度，底部壓力系數(shù)由負(fù)轉(zhuǎn)正，即底部會(huì)出現(xiàn)正推力。本文計(jì)算方案把燃?xì)庾鳛橐环N均質(zhì)流體，因此不能完全表征燃?xì)獾膶?shí)際流動(dòng)特性。未來(lái)需要開展發(fā)動(dòng)機(jī)燃?xì)舛嘟M分?jǐn)?shù)值方法的研究工作。

參考文獻(xiàn)

[1] 肖志祥, 符松. 用RANS/LES混合方法研究超聲速底部流動(dòng)[J]. 計(jì)算物理, 2009, 26(2):221-230. [Xiao Zhi-xiang, Fu Song. Study on supersonic base flow using RANS/LES method [J]. Chinese Journal of Computational Physics, 2009, 26(2): 221-230.]

[2] 高瑞澤, 閆超. LES/RANS混合方法對(duì)超聲速底部流動(dòng)的應(yīng)用[J]. 北京航空航天大學(xué)學(xué)報(bào), 2011, 37(9):1095-1099. [Gao Rui-ze, Yan Chao. LES/RANS hybird method for supersonic axisymmetric base flow[J]. Journal of Beijing University of Aeronautics and Astronautics, 2011, 37(9):1095-1099. ]

[3] 朱德華, 沈清. 鈍頭體高超聲速繞流底部失穩(wěn)特征數(shù)值模擬[J]. 力學(xué)學(xué)報(bào), 2012, 44(3):465-472. [Zhu De-hua, Shen Qing. Numerical study of the stability of hypersonic base flow over a blunt body and Apollo command module[J]. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(3): 465-472. ]

[4] Soshi K, Kozo F. Computational analysis of the characteristics of subsonic, transonic and supersonic base flows[C]. The 35th AIAA Fluid Dynamics Conference and Exhibit, Toronto, Canada，June 6-9, 2005.

[5] Pramod S, Graham V C. Numerical investigation of supersonic base flows using DES[C]. The 43rd Aerospace Science Meeting and Exhibit Conference, Reno, USA, January 10-13, 2005.

[6] Francs S, Sebastein D, Philippe G, et al. RANS-LES simulations of supersonic base flow[C]. The 44th Aerospace Science Meeting and Exhibit Conference, Reno, USA, January 9-12, 2006.

[7] Krishnendu S. Eeffect of Reynolds number on detached eddy simulation of hypersonic base flow[C]. The 45th AIAA Aerospace Science Meeting and Exhibit Conference, Reno, USA, January 8-11, 2007.

[8] Soo H P, Sumanta A. An improved formulation ofk-εturbulence models for supersonic base flow[C]. The 17th AIAA Computational Fluid Dynamics Conference, Toronto, Canada, June 6-9, 2005.

[9] Bannink W J, Bakker E M. Base flow/underexpanded exhaust plume interaction in a supersonic external flow[C]. The 8th AIAA International Space Planes and Hypersonic Systems and Technology Conference, Norfolk, USA, April 27-30, 1998.

[10] Bakker P G, Bannink W J. CFD validation for base flows with and without plume interaction[C]. The 40th Aerospace Science Meeting and Exhibit Conference , Reno, USA, January 14-17, 2002.

[11] 林敬周, 李樺. 超聲速底部噴流干擾流場(chǎng)數(shù)值模擬[J]. 空氣動(dòng)力學(xué)報(bào), 2005, 23(4):516-520. [Lin Jing-zhou, Li Hua. Numerical investigation of base jet interaction in a supersonic external flow [J]. ACTA Aerodynamics Sinica, 2005, 23(4): 516-520.]

[12] 田耀四, 蔡國(guó)飆, 朱定強(qiáng), 等. 固體火箭發(fā)動(dòng)機(jī)噴流流場(chǎng)數(shù)值仿真[J]. 宇航學(xué)報(bào), 2006, 27(5):876-880. [Tian Yao-si, Cai Guo-biao，Zhu Ding-qiang，et al. Exhaust plume simulation of solid rocket[J]. Journal of Astronautics, 2006, 27(5):876-880.]