李中華，黨雷寧，李志輝，2，李緒國

(1. 中國空氣動力研究與發(fā)展中心超高速空氣動力研究所， 綿陽621000；2. 北京航空航天大學(xué) 北京前沿創(chuàng)新中心國家計算流力學(xué)實驗室，北京100191)

1 引言

天宮飛行器隕落再入大氣層過程會經(jīng)歷稀薄過渡流區(qū)。 在強(qiáng)烈的氣動力/熱作用下會發(fā)生解體并產(chǎn)生大量碎片，其中一部分碎片在空中被完全燒毀，其余存活下來的殘骸碎片會落到地面。 過渡流區(qū)的氣動力、熱特性的準(zhǔn)確預(yù)測對航天器解體及其殘骸碎片隕落的預(yù)報至關(guān)重要。

過渡流區(qū)的流動無論試驗還是數(shù)值計算均難于處理[1-2]。 在數(shù)值仿真方面，基于連續(xù)流的CFD 技術(shù)[3-4]和廣泛應(yīng)用于稀薄流的DSMC 方法在過渡流區(qū)均會遇到各自的問題[5-7]。 在求解稀薄效應(yīng)很重的流動非平衡問題時，連續(xù)流方程會失效，CFD 的結(jié)果可能會在流場和物面力/熱參數(shù)等計算方面產(chǎn)生比較大的誤差；DSMC 方法受網(wǎng)格尺寸和時間步長的限制，計算效率很低，難以實現(xiàn)有效模擬[8-9]。 過去在計算機(jī)不夠先進(jìn)的情況下，過渡流區(qū)氣動特性的計算主要以當(dāng)?shù)鼗椒ǖ墓こ逃嬎銥橹鳌?/p>

近年來，隨著高性能計算機(jī)與數(shù)值計算技術(shù)的發(fā)展，在過渡流區(qū)，基于Boltzmann 模型方程的氣體動理論統(tǒng)一算法已研究建立[10-11]，并取得了較好的效果。 但其計算量巨大，應(yīng)用受到一定的限制。 另外一種途徑是連續(xù)流/稀薄流的耦合算法[12]。 基于CFD 和DSMC 方法的N-S/DSMC 耦合算法，能夠發(fā)揮各自的優(yōu)點(diǎn)，擴(kuò)展CFD 和DSMC方法的應(yīng)用范圍，提高DSMC 方法在模擬過渡區(qū)流動的計算效率[13-14]。 但是，對于耦合方法的研究，多以二維/軸對稱簡單外形為主，研究耦合算法的機(jī)理問題，對于三維復(fù)雜流動問題，研究文獻(xiàn)很少[15]。

本文基于MPC(Modular Particle-Continuum)耦合處理方法[16-18]，建立一種適于三維復(fù)雜外形飛行器過渡流區(qū)高超聲速繞流問題的N-S/DSMC耦合算法；并利用該耦合計算方法，對天宮外形體在低密度不同克努森數(shù)(Knudsen，Kn)條件下的繞流進(jìn)行數(shù)值計算研究，通過與相關(guān)風(fēng)洞實驗和計算結(jié)果的比較，驗證所建N-S/DSMC 耦合模型算法的有效性。

2 數(shù)值模擬方法

本文采用MPC 耦合處理技術(shù)，不需要對CFD和DSMC 程序做改動，直接在2 個獨(dú)立的程序模塊外部加入網(wǎng)格和信息交換的計算模塊，以實現(xiàn)N-S/DSMC 的耦合計算[17-18]。

把整個流場分為連續(xù)流和稀薄流區(qū)(圖1)。在連續(xù)流區(qū)域，采用時間一階、空間二階的隱式NND 格式求解式(1)所示三維非定?？蓧嚎sN-S方程[3]:

式中， U 為守恒變量矢量；E、F、G 為無粘通量矢量；Ev、Fv、Gv為粘性通量矢量。

圖1 N-S/DSMC 耦合算法分區(qū)示意圖Fig.1 Sketch of N-S/DSMC hybrid method

在稀薄流區(qū)域，采用DSMC 方法求解非平衡流場[8，19-20]。 采用式(2)所示當(dāng)?shù)豄n 數(shù)判斷連續(xù)流方程的失效問題[11，13]:

式中，λ 為當(dāng)?shù)亓鲌鰵怏w分子平均自由程；Q為當(dāng)?shù)亓鲌龊暧^流動參數(shù)，可以是壓力、密度或溫度。 取梯度最大的流動參數(shù)來計算Knl。 當(dāng)Knl≥0.02 時，在耦合邊界上，流動已經(jīng)開始出現(xiàn)非平衡現(xiàn)象，此時速度分布函數(shù)已經(jīng)不再符合Maxwell 平衡分布，認(rèn)為連續(xù)流方程失效，可以采用Chapman-Enskog 分布來描述。

對于偏離Maxwell 分布的非平衡分布，其一階展開如式(3)[21]:

式中，C 為無量綱化的分子熱運(yùn)動速度， C ＝c/(2kT/m)1/2。 Г(C)為偏離度，滿足式(4):

本文基于Maxwell 分布的速度采樣方法，在邊界面上給出符合Chapman-Enskog 分布的分子速度[21]，根據(jù)當(dāng)?shù)氐姆瞧胶舛龋捎谩熬芙^-接受”的隨機(jī)篩選準(zhǔn)則，產(chǎn)生符合Chapman-Enskog 分布的分子速度。 步驟如下:

2)設(shè)置幅度參數(shù)A ＝1 ＋30B；

3)根據(jù)Maxwell 分布產(chǎn)生一個速度矢量Ctry；

4)產(chǎn)生一個隨機(jī)數(shù)Rf。 如果ARf≤Γ(C)，接受Ctry；否則，返回到第3 步，產(chǎn)生新的Ctry；

5)給出分子速度c ＝(2kT/m)1/2Ctry＋u。

由于DSMC 方法存在較大的統(tǒng)計波動，這些統(tǒng)計波動如果傳遞到CFD 求解區(qū)域，可能會導(dǎo)致CFD 求解不穩(wěn)定。 為了抑制DSMC 方法的統(tǒng)計波動對CFD 計算的影響，本文發(fā)展了基于亞松弛的收斂處理技術(shù)。 在每個時間點(diǎn)j 上，某個網(wǎng)格中由DSMC 方法模擬得到的宏觀流動參數(shù)為Qj(包括ρ、u、v、w、T)，傳遞給CFD 的參數(shù)與上一步的參數(shù)進(jìn)行式(4)所示加權(quán)平均:

式中，θ 為網(wǎng)格中新參數(shù)的比較小的權(quán)重(本文取θ＝0.02)；相應(yīng)地，(1-θ)是比較大的權(quán)重。

3 算例結(jié)果與討論

3.1 天宮兩艙體低密度風(fēng)洞試驗耦合計算驗證

采用前述N-S/DSMC 耦合模型算法計算了天宮兩艙結(jié)構(gòu)體外形在稀薄過渡流區(qū)不同來流條件下的試驗狀態(tài)。 低密度風(fēng)洞試驗狀態(tài)見表1。 試驗氣體為氮?dú)?N2)。 每個試驗狀態(tài)分別計算-2°～25°迎角繞流流場。

表1 兩艙體外形過渡流區(qū)試驗狀態(tài)Table 1 Test cases of two-module configuration

兩艙體0°迎角典型流場壓力分布云圖如圖2。 頭部端面形成主激波，主激波打在艙體肩部，形成次激波，在肩部產(chǎn)生僅次于駐點(diǎn)區(qū)的高壓和高熱流集聚效應(yīng)，因兩艙體外形的巨大尺度，在背風(fēng)區(qū)尾部出現(xiàn)較大范圍的低壓真空區(qū)。

圖2 N-S/DSMC 計算天宮兩艙體流場壓力分布(狀態(tài)1)Fig.2 Pressure distribution of two-module body calculated by N-S/DSMC (case1)

氣動力特性耦合計算與試驗結(jié)果的比較情況如圖3。 對于軸向力系數(shù)CA，兩個試驗狀態(tài)差別比較明顯。 在試驗參數(shù)范圍內(nèi)，CA隨著稀薄效應(yīng)的增大而增加。 法向力系數(shù)CN和俯仰力矩系數(shù)(絕對值)Cm0隨迎角增加而增加，而隨稀薄度變化不大。 壓心系數(shù)隨稀薄度的變化也較小。 這些均符合過渡流區(qū)氣動力特性變化規(guī)律。 圖3 所示本文算法計算與低密度風(fēng)洞實驗2 種結(jié)果的對比表明，耦合算法結(jié)果與試驗結(jié)果吻合很好，迎角越小，耦合算法結(jié)果與低密度風(fēng)洞實驗吻合越好。

圖3 不同狀態(tài)下兩艙體氣動特性計算與風(fēng)洞試驗比較Fig.3 Comparison of aerodynamic characteristics between the present computation and the windtunnel experiment for various cases of twomodule configuration

3.2 天宮多次解體殘骸兩體間干擾流場耦合計算

天宮多次解體后，會產(chǎn)生各種類型的殘骸碎片，相互之間的流場會有干擾。 簡單起見，本文計算了球(D ＝100 mm)和球柱(D ＝100 mm、L ＝200 mm)2 種外形在表2 所示來流條件決定的流場干擾特性。 兩體的位置為頭部并齊，平行放置，軸線間距離為d。

表2 兩體干擾流場狀態(tài)計算參數(shù)Table 2 Cases of interfere flow field

不同干擾距離條件下兩個球體干擾壓力場等值線云圖分布如圖4。 在d/D ＝1.25 條件下，激波在對稱面上形成干擾激波，干擾激波很強(qiáng)，打在球體表面后形成反射激波，并在對稱面上再次反射，形成復(fù)雜波系。 在d/D ＝1.5 條件下，反射激波與尾流區(qū)的膨脹波相互作用，在尾流區(qū)形成復(fù)雜波系。

兩體干擾對氣動特性的影響如圖5。 從圖中可以看出，隨著稀薄度的增大，軸向力系數(shù)增加。在Ma ＝10 的條件下，當(dāng)d/D ＝1.25 時，兩體之間的干擾對氣動特性有一定影響，兩種稀薄度條件下，CA均有一定下降，這是因為干擾激波打在球體的下、后表面，在下、后表面產(chǎn)生高壓區(qū)域(參見圖4 所示)，向前產(chǎn)生一定推力，使得CA下降。這個高壓區(qū)域同時產(chǎn)生法向力和俯仰力矩，并使壓心系數(shù)后移。 從CN系數(shù)隨不同來流Kn 數(shù)變化來看，稀薄度越大，干擾的影響越大。 當(dāng)d/D≥1.5 時，干擾對氣動特性基本不產(chǎn)生影響。 這為空氣動力融合軌道數(shù)值預(yù)報提供實時計算工程設(shè)計準(zhǔn)則。

圖4 天宮解體殘骸球體繞流不同間距干擾流場(Case2)Fig.4 Interference flow field between spheres of varies distances (case2)

圖5 球體干擾對氣動特性的影響Fig.5 Effects of interference on aerodynamic characteristics around two side-by-side spheres

圖6 是天宮多次解體殘骸簡化物兩個并排球柱體繞流的流場結(jié)構(gòu)。 d/D ＝1.25 時，干擾激波在對稱面和柱體間進(jìn)行多次反射，殘骸間繞流干擾嚴(yán)重；d/D＝1.5 時，干擾激波在對稱面和柱體間只有一次反射，殘骸間繞流干擾迅速減小。

圖6 球柱體繞流不同間距干擾流場(Case2)Fig.6 Interference flow field between sphere-poles of varies distances (case2)

圖7 繪出干擾對氣動特性的影響。 從圖中可以看出，與球體類似，隨著稀薄度的增大，軸向力系數(shù)增加。 在Ma ＝10 的來流條件下，當(dāng)d/D ＝1.25 時，兩體之間的干擾對氣動特性有一定影響，兩種稀薄度條件下，干擾距離越小，兩體間繞流干擾越嚴(yán)重。 殘骸間干擾使CA、CN、Cm增加，壓心位置Xcp后移；當(dāng)d/D＝2 時，干擾基本不對氣動特性產(chǎn)生影響。

圖7 球柱體干擾對氣動特性的影響Fig.7 Effects of interference on the aerodynamic characteristics of sphere-poles

對以上結(jié)果分析表明:物體的形狀、殘骸間距離、來流稀薄條件等因素會對殘骸碎片氣動特性產(chǎn)生不同的干擾影響，基于表2 計算條件，一般在兩體相距1.5 倍之外，殘骸碎片間繞流干擾在再入解體飛行航跡數(shù)值預(yù)報設(shè)計中，可以忽略不計其對氣動特性的影響。

4 結(jié)論

1) 通過對CFD 和DSMC 方法計算程序進(jìn)行基于網(wǎng)格與信息交換的雙向耦合設(shè)計，發(fā)展了求解過渡流區(qū)高超聲速三維繞流問題的N-S/DSMC耦合算法。

2) 采用耦合算法，對天宮飛行器兩艙體在低密度風(fēng)洞試驗條件下的氣動力特性開展了仿真，與試驗結(jié)果進(jìn)行對比，驗證了本文所建立的耦合算法對數(shù)值求解類天宮飛行器服役期滿再入解體過程所產(chǎn)生殘骸碎片多體繞流干擾的可靠可行性。

3) 通過對天宮解體殘骸碎片簡化外形的兩體(球體、球柱體)干擾流場進(jìn)行仿真計算，結(jié)果表明，對于球、球柱這類簡化外形，在一定的距離范圍內(nèi)，兩體干擾會對氣動特性產(chǎn)生很強(qiáng)的影響，隨著兩體間隔距離的增大，如簡單球體間隔1.5倍、球柱體間隔2 倍直徑，殘骸碎片間繞流干擾引起的氣動特性影響可以忽略不計，為天宮飛行器實際再入解體氣動融合結(jié)構(gòu)/軌道數(shù)值預(yù)報提供重要設(shè)計規(guī)范標(biāo)準(zhǔn)。