稀薄過渡流區(qū)橫向噴流干擾效應數(shù)值模擬研究

梁 杰，閻 超，李志輝，李緒國

（1．北京航空航天大學 航空科學與工程學院，北京 100191；2．中國空氣動力研究與發(fā)展中心，四川 綿陽 621000）

0 引 言

飛行器在飛行過程中因機動需要快速改變飛行姿態(tài)或軌道，通常采用反作用控制系統(tǒng)（RCS）來實現(xiàn)，分布在飛行器表面不同位置的RCS發(fā)動機提供的推力完成俯仰、偏航和滾轉等姿態(tài)變化以及軌道的轉移［1］。RCS發(fā)動機工作時產生的噴流氣體與周圍氣流發(fā)生相互作用，對推力和飛行器的壓心位置會產生附加影響，并改變飛行器表面的熱環(huán)境。在連續(xù)流區(qū)域，經過幾十年持續(xù)的風洞試驗和數(shù)值模擬研究［2－4］，對不同類型的噴流干擾及流動機理有了非常細致的了解。但在稀薄大氣環(huán)境下，尤其是飛行器處于稀薄過渡流域時，由于連續(xù)介質假設失效，飛行器表面出現(xiàn)了滑移效應，噴流與前方邊界層（粘性層）相互作用引起的流動分離和流場結構均出現(xiàn)了不同程度的變化。無論是地面風洞試驗所用的相似準則，還是數(shù)值求解連續(xù)介質假設下的NS方程，都已無法滿足精確模擬稀薄過渡流區(qū)橫向噴流干擾特征的要求。另外，在高空稀薄氣流中噴流會形成羽流，它所涉及的問題還包括對周圍環(huán)境的污染［5］，影響目標的識別。在文獻［6］中，我們詳細研究了三維平板模型在不同的噴流壓力比環(huán)境下音速噴流與稀薄大氣干擾的流動規(guī)律，計算得到的干擾流場結構、流動分離特征和分離長度與低密度風洞輝光放電和油流測量結果有較好的一致性，確認了建立的直接模擬Monte Carlo（DSMC）算法的可靠性和計算精度。由于旋成體外形的飛行姿態(tài)不同，可能存在攻角，也可能還存在側滑角，飛行器外流場的復雜性導致旋成體上橫向噴流干擾流場較平板的噴流干擾要復雜得多。更為重要的差別在繞旋成體流動的分離邊界，很可能從上表面繞到下表面，產生所謂的包裹效應［4］。這些復雜的橫向噴流流動現(xiàn)象難以進行風洞試驗模擬，都需要依靠數(shù)值模擬開展細致的研究分析。本文在對平板模型噴流干擾研究的基礎上，重點開展了細長鈍雙錐外形高超聲速稀薄來流與超聲速／高超聲速橫向噴流干擾效應的研究，計算分析了不同高度、不同速度、不同攻角、不同噴流推力下復雜流場結構和對氣動力特性的影響規(guī)律?？疾炝薘CS噴管出口參數(shù)不同（均勻／非均勻）對噴口附近分離渦和分離長度的影響。計算中采用了研究建立的兩級直角網格及碰撞和取樣網格自適應算法［7］，以及MPI并行環(huán)境下的基于靜態(tài)隨機負載平衡技術的DSMC并行算法。

1 DSMC數(shù)值方法

1．1 DSMC方法

DSMC方法［8］是用若干個模擬分子代替真實氣體的分子，并在計算機上存貯模擬分子的位置坐標、速度分量以及內能，這些量隨模擬分子的運動、與邊界的碰撞以及分子之間的碰撞而改變，最后通過統(tǒng)計網格內模擬分子的運動狀態(tài)實現(xiàn)對真實氣體流動問題的模擬。該方法的關鍵是在時間步長Δt內將分子的運動與碰撞解耦。在Δt時間內讓所有分子運動一定的距離并考慮在邊界的反射，然后計算此Δt內具有代表性的分子間的碰撞。

1．2 計算網格及網格自適應方法

在DSMC方法中，流場中的網格是用來選取可能的碰撞分子對以及對宏觀流動參數(shù)取樣。流場采用均勻的直角坐標網格，追蹤分子的效率非常高，計算區(qū)域內的模擬分子可以直接根據(jù)分子的位置坐標來確定分子所屬的網格，而不必跟蹤分子從一個網格運動到另一個網格［7］。其缺點是無法精確地描述物面邊界。本文結合直角網格計算的高效率以及表面非結構網格對飛行器幾何外形的精細描述，建立了混合網格結構的DSMC數(shù)值模擬方法。在描述物面幾何形狀的非結構網格建立以后，直接將其嵌入到直角網格的流場中，使DSMC計算對流場網格的依賴程度大大降低。同時通過判斷分子運動軌跡方程和物面三角形面元上任一點的位置方程，唯一確定出分子與物面的碰撞點坐標，解決了這種混合網格流場分子運動與物面碰撞的難題。對分子在物體三角形面元上碰撞、反射前后的流場參數(shù)進行統(tǒng)計取樣就可以獲得飛行器的整體氣動力特性以及表面力、熱載荷分布。

稀薄過渡流區(qū)的噴流干擾流場結構非常復雜，流場中的氣體密度、壓力有幾個量級的變化。噴流核心區(qū)和激波干擾區(qū)的密度非常高，其它大部分區(qū)域的密度相對較低。為了解決流場中密度急劇變化的流動特征，采用了網格自適應的方法。即在背景網格的基礎上，根據(jù)流場中密度梯度的變化分別對碰撞網格和取樣網格進行細化。由于流場的梯度沿各個方向的變化是不同的，梯度變化大的方向網格細分的更密一些，因此沿三個坐標方向是各自獨立地進行自適應，碰撞分子則是在自適應后最小的亞網格內選取，保證了計算的空間精度。

碰撞網格自適應的目的是使亞網格的尺度與當?shù)仄骄杂沙桃恢拢泳W格自適應是為了捕捉當?shù)氐牧鲃犹荻?，正常情況下取樣網格的數(shù)量遠遠小于碰撞網格的數(shù)量。為了快速定位自適應后分子所屬的網格，采用如下的搜索方法：

對于流場中均勻的矩形網格，模擬分子在背景網格中的位置可以用下面的編號表示：

式中ij是在（x，y，z）三個坐標方向的序號，rj是分子坐標，Δj＝Lj／Nj是網格間距，Lj是沿相應坐標方向的流場大小，Nj是沿相應坐標方向劃分的背景網格數(shù)，j0是背景網格的起始坐標。則模擬分子在背景網格中的序號為：

每個背景網格根據(jù)當?shù)氐牧鲃犹荻妊厝齻€坐標方向分別劃分成lx，ly，lz個二級碰撞或取樣網格，二級網格在一級（背景）網格中的起始編號為：

因為二級網格也是均勻的矩形網格，因此二級網格的序號也可以按照公式（2）同樣的方式獲得：

1．3 隨機負載平衡的DSMC并行算法

采用DSMC方法模擬三維RCS噴流與稀薄大氣的干擾流場，由于噴流核心區(qū)的氣流密度相當高，導致模擬分子數(shù)和網格劃分引起的計算量仍然非常大，需要開展并行計算。DSMC并行算法采用區(qū)域分解，根據(jù)計算的處理器數(shù)（或CPU核心數(shù)）的多少將計算區(qū)域劃分為等量的子區(qū)域，每個處理器在其分配的子區(qū)域內部獨立地計算模擬分子的碰撞和遷移，離開子區(qū)域的模擬分子把攜帶的信息傳遞給對應子區(qū)域的處理器。并行計算總的時間包括每個處理器計算碰撞、遷移的時間、處理器之間的通訊時間以及各處理器之間為同步而等待的時間。提高并行計算的效率主要是通過減少通訊和同步等待的時間來實現(xiàn)。

如果模擬分子在一個時間步長之后，離開原來所在的網格遷移進入相鄰的網格，并且改變后的網格隸屬于其它的處理器，這樣的分子及其攜帶的所有信息就被記錄下來并放入緩沖數(shù)組，當所有的模擬分子都被確定以后，放入緩沖數(shù)組的分子就在相應的處理器之間進行信息交換。

DSMC方法的一個顯著特征是模擬分子在流場中的分布是不斷變化的，計算時間依賴于子區(qū)域內的分子個數(shù)、分子的碰撞次數(shù)以及分子與物面的碰撞次數(shù)。因為模擬前不可能預測到不同子區(qū)域內的分子個數(shù)，因此在區(qū)域分解時將各子區(qū)域的計算負載劃分得比較均衡就特別困難。如果簡單地將流場均勻劃分成若干同等大小的子區(qū)域，就有可能造成處理器之間的負載分布不均勻，負載最少的處理器就要化費較長的時間等待負載最大的處理器，造成同步等待的時間顯著增加。這種方式并不適合DSMC并行模擬。

為了減少不同處理器計算時間的差別，通常采用負載平衡技術來實現(xiàn)［9－10］，本文根據(jù)稀薄大氣條件下飛行器橫向噴流干擾特點，發(fā)展了一套靜態(tài)隨機負載平衡方法用于解決不同處理器之間的計算時間同步問題，該方法基于概率近似原理，將計算區(qū)域的全部網格平均分配給指定的所有處理器。由于采用相同的概率隨機選取流場網格，按照均分后的數(shù)量分配給每個處理器，因此當計算區(qū)域的網格數(shù)量較多時，每個處理器包含近似相等的物體邊界、高密度流動區(qū)域和稀薄氣體區(qū)域的網格數(shù)，這樣每個處理器的計算負載也非常接近。當用于并行計算的CPU數(shù)量不多時（Np＜32），這種方式得到的并行算法加速比和并行效率都非常高［6］。

2 計算結果分析

2．1 平板模型噴流干擾計算與試驗對比

為了驗證本文算法的精確度，計算模擬了三維平板模型橫向噴流干擾的分離距離并與低密度風洞試驗結果進行了對比。風洞試驗氣體為氮氣，前室總溫593K、總壓400kPa，流場馬赫數(shù)12．99。音速噴流的氣體為氦氣，噴嘴前室總溫為300K，模擬的三個試驗狀態(tài)的總壓變化見表1。在DSMC數(shù)值模擬中，采用VHS分子模型，L－B內能交換模型，壁面采用完全漫反射模型（下同）。平板模型長120mm，寬68mm，高15mm。音速噴嘴出口直徑2mm，在平板中心距前端70mm的位置。

噴流干擾的分離長度是通過測量噴管出口中心到主分離線的徑向距離得到的。用DSMC方法計算的三個試驗狀態(tài)的分離長度，試驗測量的分離長度都列在了表1中。計算與試驗總體上吻合較好，只是計算得到的噴流干擾流場分離長度比試驗測量的分離長度偏小些，最大偏差3mm（Case3），最大相對偏差10%（Case1）。這是由于噴嘴尺寸較小，為簡化起見，計算中音速噴嘴出口參數(shù)設置為均勻流，即出口只有軸向速度，徑向速度為0。這樣會減弱噴流沿徑向的膨脹速度，可能會造成由于逆向流動減弱而引起計算的分離距離偏小，噴管出口參數(shù)不同對分離長度的影響在下節(jié)的算例中還會進一步分析說明。

表1 DSMC計算模擬的試驗狀態(tài)及分離距離Table 1 Test conditions and separation distance

模擬的激波／邊界層干擾、流動分離特征以及表面極限流線與低密度風洞輝光放電和表面油流圖譜的對比分析詳見文獻［6］。這些結果驗證了本文算法的可靠性以及對稀薄過渡流區(qū)復雜橫向噴流干擾流場的模擬能力。

2．2 鈍雙錐橫向噴流與稀薄大氣干擾特性數(shù)值模擬

下面選取總長度2．1m的細長鈍雙錐外形作為研究對象，噴流中心距頭部1m，方向與Y軸正向相同。這樣就可以考察頭部斜激波與噴流弓形主激波干擾引起的復雜流動結構。RCS發(fā)動機分別采用推力為150N和10N的型面噴管，150N噴管的出口馬赫數(shù)近似為4．4，燃氣組分為 H2、N2、NH3和 H2O，10N噴管的出口馬赫數(shù)6．3，燃氣組分同上。模擬的大氣參數(shù)組分為N2和O2。

2．2．1 飛行高度對噴流干擾流場的影響

圖1分別給出了噴流推力150N時85km、90km和100km高度、飛行速度6000m／s的干擾流場密度分布等值線云圖。飛行高度較低時，飛行器頭部產生較強的斜激波，頭激波與噴流弓形主激波相交產生激波／激波空間干擾，噴流與前方邊界層（粘性層）內氣體相互作用形成三維分離流動區(qū)。噴流弓形激波后為高壓力和高密度的強剪切層區(qū)，由于噴流阻礙了來流，在噴流的下游存在較寬的低密度和低壓區(qū)。

隨著飛行高度的增加，來流氣體的稀薄度增大，頭部斜激波變弱，噴流與稀薄來流的干擾也越來越弱，噴流自身的影響區(qū)域卻越來越大，在100km以上高空，噴流的影響已經覆蓋整個飛行器。如果噴流的燃氣組分中含有較強污染效應的成分（如NH3和H2O），會粘附在飛行器表面，從而影響傳感器的正常工作。另外，隨著飛行高度的增加，噴流與粘性層內氣體的干擾減弱，引起的分離也逐漸減小，100km的計算狀態(tài)已經沒有分離的形成。

圖1 不同高度150N推力橫噴干擾流場密度分布Fig．1 Density distribution of lateral jet interaction flowfield of 150Nthrust in different altitudes

2．2．2 噴管出口參數(shù)分布對噴流干擾流場影響

為方便起見，前面的計算中噴管出口參數(shù)采用的是均勻分布，沿徑向的氣體速度為零，但實際上燃氣在噴管內流動到達出口時分布并不是均勻的。圖2、圖3分別給出了采用真實噴管非均勻出口參數(shù)（求解NS方程獲得）和均勻分布計算的90km高度噴流干擾分離區(qū)附近的流線局部放大圖，不論是主分離渦還是馬蹄渦均有明顯差別，采用真實噴管型面出口參數(shù)計算的分離區(qū)明顯大于采用均勻出口參數(shù)計算的分離區(qū)。

另外，從圖4、圖5分別繪出的采用噴管出口非均勻與均勻參數(shù)計算得到的表面極限流線分布，也可以看出兩者的不同。在圖4中數(shù)字1代表主分離線，2代表再附線，3是馬蹄渦引起的二次分離線，與圖2中的空間流線相對應。圖5中僅可看到主分離線和再附線，同時主分離線和再附線之間的分離區(qū)范圍也減小了。因此，噴管出口參數(shù)分布對分離區(qū)的大小和形狀有明顯的影響。

2．2．3 來流速度及噴管推力對干擾流場的影響

為了比較飛行器來流速度不同對飛行器橫向噴流干擾流場的影響，圖6給出了來流速度4000m／s、高度85km、噴流推力150N的干擾流場密度分布，與圖1（a）相比，來流速度的降低使頭部斜激波減弱，噴流主激波向上游擴展增強，噴流產生的逆向流動引起的分離距離和空間分離渦更大一些。圖7、圖8繪出了噴流推力為10N時來流速度分別為4000m／s和6000m／s的噴流干擾流場密度分布，飛行高度為85km，噴流馬赫數(shù)6．3。從圖中可以看出，盡管噴流馬赫數(shù)較高，但由于推力小，噴流產生的逆流與粘性層內氣體的干擾較弱，并沒有引起明顯的流動分離。而且，對于小推力噴流來說，來流速度的增加，使噴流主激波向下游傾斜得更大，對噴口下游區(qū)域的再壓縮增強，使下游飛行器表面壓力增大，對力矩產生的附加影響有所加大。計算表明，對于推力比較小的噴流，在來流參數(shù)不變的條件下，來流對噴流的壓制作用要明顯高于大推力噴流。

2．2．4 攻角變化對噴流干擾流場的影響

旋成體細長雙錐外形的飛行姿態(tài)不同，如正攻角和負攻角所對應的干擾流場特性是不同的。10°攻角時（高度90km、來流速度6000m／s、噴流推力150N），由于頭激波遠離噴口，對噴口影響較小，與噴流弓形激波的作用位置上移，粘性層內氣體的流動減弱，噴流的逆向流動引起的空間分離渦相比0°攻角時有所增強（圖9）。同時下表面的氣流向上方卷起，使分離產生的包裹效應減弱，表面分離形成的高壓區(qū)減小。

圖6 150N推力橫噴干擾流場密度分布（V＝4000m／s）Fig．6 Density distribution of lateral jet interaction flowfield of 150Nthrust（V＝4000m／s）

圖7 10N推力橫噴干擾流場密度分布（V＝4000m／s）Fig．7 Density distribution of lateral jet interaction flowfield of 10Nthrust（V＝4000m／s）

圖8 10N推力橫噴干擾流場密度分布（V＝6000m／s）Fig．8 Density distribution of lateral jet interaction flowfield of 10Nthrust（V＝6000m／s）

圖9 10°攻角時噴口上游流線及分離渦Fig．9 Stream line and seperation vortex of 10°AOA

圖10 －10°攻角時噴口上游流線分布Fig．10 Stream line and seperation vortex of－10°AOA

圖10 是－10°攻角的情況。這時的頭激波靠近噴口，并產生較強的壓制作用。噴口附近的當?shù)貕毫Ρ葴p小，噴流干擾引起的噴口上游空間的分離基本消失，但由于氣流是向下表面方向流動，空間流動引起的包裹效應在表面引起的高壓力區(qū)域有所增大。

2．2．5 噴流干擾對氣動力特性的影響

對于旋成體上的側向噴流干擾引起的流動分離，很可能從上表面繞到下表面跨過φ＝90°的分界線（定義φ＝0°為上表面迎風母線），假如噴流干擾的壓力正增量區(qū)（在主分離線之后的相當大區(qū)域）包裹在旋成體上，那么在φ＜90°范圍內的正增量投影將放大噴管推力，對應φ＞90°的干擾區(qū)，雖然表面壓力仍具有正增量，但其在推力方向的投影為負增量，因此產生所謂的包裹效應。通常用放大系數(shù)來定義干擾引起的力特性變化量［4］：

推力放大系數(shù) KF＝（T＋ΔT）／T

力矩放大系數(shù) KM＝（M0＋ΔM）／M0

ΔT與ΔM表示由于噴流干擾引起的推力與力矩增量，本文中力矩參考點為頭部頂點。

表2給出了噴管推力為150N時噴流干擾引起氣動力系數(shù)的增量以及計算的推力和力矩放大系數(shù)，來流速度為6000m／s。表中數(shù)據(jù)表明，在相同的飛行狀態(tài)下，隨著飛行高度的增加，噴流干擾引起的氣動力系數(shù)的增量在數(shù)值上是不斷增大的，但推力和力矩放大系數(shù)卻是減小的趨勢，主要是由于來流密度降低使動壓減小，造成推力和力矩的增量減小。另外，正攻角引起推力和力矩放大系數(shù)的減小，負攻角引起推力放大系數(shù)的增大，力矩放大系數(shù)的減小。噴管推力10N時噴流干擾對氣動力特性的影響，與150N推力噴流干擾的變化規(guī)律類似，只是當推力比較小時，推力和力矩放大系數(shù)有所增大，主要是處于分母的推力和力矩比較小造成的，與氣動中心的Φ1m高超聲速風洞的側向噴流干擾試驗有相同的結論。

表2 噴流干擾對氣動力特性的影響（150N噴管）Table 2 The effects of jet interaction on aerodynamic properties（150Nthrust）

3 結 論

綜合分析上述數(shù)值模擬結果，可以得出以下幾點初步研究結論：

（1）在噴流參數(shù)保持不變的情況下，隨著來流稀薄度的增加，因噴流干擾引起的空間分離渦不斷減弱直至消失，稀薄來流對噴流的影響越來越弱，而噴流自身的影響區(qū)域卻越來越大。

（2）采用真實RCS噴管出口參數(shù)計算的空間分離渦和表面極限流線與采用均勻出口參數(shù)有明顯的不同，噴管出口參數(shù)分布對分離區(qū)的大小、形狀以及表面壓力分布都有一定的影響。

（3）高超聲速來流對小推力的RCS噴流有較強的壓制作用，噴口附近的干擾流場結構及分離特征與當?shù)亓鲃訔l件密切相關。

（4）來流速度的增加，使噴流主激波更向下游傾斜，并對噴口下游的流動再壓縮增強。

（5）飛行姿態(tài)不同，如正攻角和負攻角所對應的干擾流場特性明顯不同，正攻角引起推力和力矩放大系數(shù)的減小，負攻角引起推力放大系數(shù)的增大，力矩放大系數(shù)的減小。

（6）在相同飛行狀態(tài)下，隨著飛行高度的增加，噴流干擾引起的氣動力系數(shù)的增量在數(shù)值上是不斷增大的，但推力和力矩放大系數(shù)卻是減小的趨勢。

（7）在稀薄流域，RCS噴流干擾引起的推力和力矩放大系數(shù)總體上是比較小的。

［1］ GIMELSHEIN S F，ALEXEENKO A A，LEVIN D A．Modeling of chemically reacting flows from a side jet at high altitudes［R］．AIAA 2002－0212，2002．

［2］ BRANDEIS J，GILL J．Experimental investigation of side jet steering for missiles at supersonic and hypersonic speeds［R］．AIAA 95－0316，1995．

［3］ SRIVASTAVA B．CFD analysis and validation of lateral jet control of a missile［R］．AIAA 96－0288，1996．

［4］ 李素循．激波與邊界層主導的復雜流動［M］．北京：科學出版社，2006：162－190．（LI suxun．Shock wave and boundary layer dominant complicated flow［M］．Beijing：Science Press，2006：162－190（in Chinese））

［5］ GIMELSHEIN S F，ALEXEENKO A A，LEVIN D A．Modeling of the interaction of a side jet with a rarefied atmosphere［R］．AIAA 2001－0503，2001．

［6］ 梁杰，閻超，楊彥廣，等．過渡區(qū)側向噴流干擾的并行DSMC數(shù)值模擬研究［J］．宇航學報，2011，32（5）：1012－1018．（LIANG Jie，YAN Chao，YANG Yanguang，et al．Parallel DSMC simulation of lateral jet interaction in rarefied transitional region［J］．Journal of Astronautics，2011，32（5）：1012－1018（in Chinese））

［7］ 梁杰，閻超，杜波強．基于兩級直角網格結構的三維DSMC算法研究［J］．空氣動力學學報，2010，28（4）：466－471：466－471．（LIANG Jie，YAN Chao，DU Boqiang．An algorithm study of three－dimensional DSMC simulation based on two－level Cartesian coordinates grid structure［J］．Acta Aerodynamica Sinica，2010，28（4）：466－471（in Chinese））

［8］ BIRD G A．Molecular gas dynamics and the direct simulation of gas flows［M］．Oxford：Clarendon Press，1994．

［9］ IVANOV M，MARKELOV G，TAYLOR S，et al．Parallel DSMC strategies for 3Dcomputation［A］．Proc．Par－allel CFD＇96［C］，1996：485－492．

［10］LEBEAU G J．A parallel implementation of the direct simulation Monte Carlo method［J］．Comput．Methods Appl．Mech．Engrg．，1999，174：319－337．