梁 杰，李志輝，杜波強，方 明

(中國空氣動力研究與發(fā)展中心超高速所，綿陽621000)

0 引言

大型航天器如空間站或衛(wèi)星在壽命即將結(jié)束時，地面控制中心利用剩余的推進(jìn)劑進(jìn)行多次剎車離軌控制，將其軌道高度逐漸減低，然后在無控狀態(tài)下再入大氣層。當(dāng)再入到120～100 km左右高度時，在大氣的作用下開始?xì)鈩蛹訜帷?00 km高度左右，太陽翼電池帆板、大型天線等設(shè)備在氣動熱和氣動力的雙重作用下開始撕裂，并與航天器本體分離。隨著航天器表面溫度持續(xù)升高，航天器局部材料出現(xiàn)軟化、結(jié)構(gòu)失效，艙體內(nèi)外部設(shè)備陸續(xù)與航天器發(fā)生分離。解體后的部件在墜落過程中熔化或燃燒，質(zhì)量逐漸減小，一些較大的部件則會殘留一部分，直到撞擊地面。所有的碎片沿航天器的軌跡方向散布在一個狹長的區(qū)域［1］。如果解體后的碎片散落在人口密集區(qū)，則會對人類的安全造成破壞［2］。這類危險目標(biāo)的隕落跟蹤與落點預(yù)報均會引起世界各國的重視，如過去的美國天空實驗室、原蘇聯(lián)宇宙1402號核動力衛(wèi)星、我國失控衛(wèi)星以及俄羅斯“和平”號空間站的隕落，均受到世界各國的重視，提前進(jìn)行了再入風(fēng)險評估［3］。國外開展空間碎片再入預(yù)測及地面風(fēng)險評估的研究已有十多年的歷史，并利用工程預(yù)測方法形成了相對成熟的軟件系統(tǒng)，如美國的碎片評估軟件(Debris Assessment Software，DAS)［4］和目標(biāo)再入存活分析工具(Object Reentry Survival Analysis Tool，ORSAT)［5］，歐洲的航天器大氣再入與氣動加熱解體(Spacecraft Atmosphere Reentry and Aerothermal Breakup，SCARAB)［6］等。國內(nèi)的清華大學(xué)利用工程預(yù)測手段初步建立了空間碎片再入分析方法和評估軟件［7］。我國的天宮一號目標(biāo)飛行器由實驗艙和資源艙組成，長約7．8 m的太陽翼電池帆板位于資源艙的兩側(cè)，兩個艙體的總長約10．5 m，艙體最大直徑約3．4 m，設(shè)計在軌壽命2年。自2011年發(fā)射升空后，已經(jīng)超期服役即將墜入大氣層損毀，需要預(yù)先對其隕落的風(fēng)險進(jìn)行詳細(xì)的評估，同時還要不斷完善這種大型航天器的再入解體和地面風(fēng)險評估手段，為后續(xù)的貨運飛船、空間站、大型衛(wèi)星的隕落風(fēng)險預(yù)報做準(zhǔn)備。

關(guān)于太陽翼在軌運行的熱分析，已有一些文獻(xiàn)開展過相關(guān)研究［8－9］，但對于大型航天器離軌隕落解體過程中太陽翼帆板的力熱分析國內(nèi)尚屬空白。而在航天器再入隕落的初期，流動處于稀薄過渡流區(qū)域，求解N－S方程的連續(xù)流計算方法已經(jīng)失效，需要采用基于分子碰撞理論的直接模擬蒙特卡洛(Direct Simulation Monte Carlo，DSMC)方法［10］。文獻(xiàn)［11］采用DSMC方法對鈍體返回艙外形在稀薄流域的高溫真實氣體效應(yīng)進(jìn)行了計算分析，文獻(xiàn)［12］則采用并行DSMC技術(shù)開展了過渡流區(qū)復(fù)雜噴流干擾流動的模擬，尚未針對帶太陽翼帆板的大型航天器存在復(fù)雜激波干擾環(huán)境下的氣動特性開展相關(guān)研究，因此本文在構(gòu)建模擬大型復(fù)雜航天器高溫氣體熱化學(xué)非平衡特性的DSMC數(shù)值方法的基礎(chǔ)上，對壽命末期離軌隕落的低軌航天器再入稀薄流域120～90 km高度的氣動力、熱特性進(jìn)行了模擬分析，重點分析了作用在水平和垂直放置的單側(cè)太陽翼帆板上的力、力矩和峰值壓力、熱流，為下一步的飛行器結(jié)構(gòu)解體分析提供數(shù)據(jù)支撐。

1 DSMC數(shù)值模擬方法

DSMC方法是通過模擬仿真分子的運動、與邊界的碰撞以及分子之間的碰撞隨時間的變化，并存貯仿真分子的位置坐標(biāo)、速度分量以及內(nèi)能，最后通過統(tǒng)計網(wǎng)格內(nèi)仿真分子的運動狀態(tài)實現(xiàn)對真實氣體流動問題的模擬。每個仿真分子代表1012～1014個真實氣體的分子，該方法的關(guān)鍵是在時間步長Δt內(nèi)將分子的運動與碰撞解耦。在Δt時間內(nèi)讓所有分子運動一定的距離并考慮在邊界的反射，然后計算此Δt內(nèi)具有代表性的分子間的碰撞。

1．1計算網(wǎng)格及網(wǎng)格自適應(yīng)

在DSMC方法中，流場中的網(wǎng)格是用來選取可能的碰撞分子對以及對宏觀流動參數(shù)取樣?！疤鞂m一號”目標(biāo)飛行器結(jié)構(gòu)復(fù)雜，涉及到局部凸起物、薄壁和細(xì)支架等需要精細(xì)描述的復(fù)雜壁面，同時計算區(qū)域需要覆蓋飛行器本體、展開的太陽翼帆板以及激波干擾引起的復(fù)雜流動區(qū)域，因此網(wǎng)格的計算效率和高可靠性至關(guān)重要。在早期的研究中作者采用兩級直角網(wǎng)格結(jié)構(gòu)對美國“獵戶座”飛船在稀薄流域的氣動力系數(shù)進(jìn)行了計算，并與國際知名軟件DS3V和DAC的計算結(jié)果一致［13］。針對“天宮一號”復(fù)雜外形，本文采用計算效率高的流場直角網(wǎng)格和對飛行器幾何外形描述精確的表面三角形非結(jié)構(gòu)網(wǎng)格(如圖1)，建立了混合網(wǎng)格結(jié)構(gòu)的DSMC數(shù)值模擬方法。在描述物面幾何形狀的非結(jié)構(gòu)網(wǎng)格建立以后，直接將其嵌入到直角網(wǎng)格的流場中，使DSMC計算對流場網(wǎng)格的依賴程度大大降低。同時通過判斷分子運動軌跡方程和物面三角形面元上任一點的位置方程，唯一確定出分子與物面的碰撞點坐標(biāo)，解決了這種混合網(wǎng)格流場分子運動與物面碰撞的難題。對分子在物體三角形面元上碰撞、反射前后的流場參數(shù)進(jìn)行統(tǒng)計取樣就可以獲得飛行器的整體氣動力特性以及表面力、熱載荷分布。

圖1 計算網(wǎng)格及網(wǎng)格自適應(yīng)Fig．1 Computational grid and self-adaption

為了解決流場中因激波壓縮、激波干擾以及氣體膨脹后出現(xiàn)的密度梯度急劇變化的流動特征，計算中采用了網(wǎng)格自適應(yīng)的策略［14］。即在背景網(wǎng)格的基礎(chǔ)上，根據(jù)流場中密度梯度的變化分別對碰撞網(wǎng)格和取樣網(wǎng)格進(jìn)行細(xì)化(如圖1)。由于流場的梯度沿各個方向的變化是不同的，梯度變化大的方向網(wǎng)格細(xì)分的更密一些，因此沿三個坐標(biāo)方向是各自獨立地進(jìn)行自適應(yīng)，碰撞分子則是在自適應(yīng)后最小的亞網(wǎng)格內(nèi)選取，保證了計算的空間精度。

1．2分子轉(zhuǎn)動和振動松弛模型

天宮飛行器再入隕落的速度達(dá)到7．6 km/s，需要考慮氣體分子的轉(zhuǎn)動和振動激發(fā)以及化學(xué)反應(yīng)。連續(xù)流氣體分子的轉(zhuǎn)動松弛碰撞數(shù)和振動松弛碰撞數(shù)需要經(jīng)過內(nèi)自由度的修正才可用于DSMC的模擬［15］，修正如下:

式中:ζt、ζr和ζv分別是分子平動、轉(zhuǎn)動和振動自由度。Parker［16］給出了連續(xù)流轉(zhuǎn)動松弛碰撞數(shù)ZCr的表達(dá)式

式中:T*是氣體分子作用勢的特征溫度，(Zr)∞是實驗測定的極限值。

連續(xù)流的振動松弛碰撞數(shù)定義為

式中:τc為分子平均碰撞時間，τMW為振動松馳時間，τP為高溫修正項。

由Landau-Teller理論結(jié)果給出

Millikan和White［17］在高溫激波管內(nèi)通過干涉儀觀測了氣體分子振動松馳過程，在104K的溫度范圍內(nèi)，給出了如下擬合參數(shù)

式中:mr為分子折合質(zhì)量，θv是振動特征溫度。Park［18］給出了高溫修正的 τP為

式中:n是數(shù)密度，σv是分子振動碰撞截面(=10－20m2)，ˉc是分子平均熱運動速度。

1．3化學(xué)反應(yīng)模型及非彈性碰撞模擬策略

通常反應(yīng)速率常數(shù)可以寫成下面的方程形式:

式中:Ea是反應(yīng)中需要的活化能，k是Boltzmann常數(shù)，Λ和η是常數(shù)，由實驗定出。

化學(xué)反應(yīng)幾率可以推導(dǎo)為下面的表達(dá)式［10］

式中:ε是對稱因子，對于不同類分子ε=1，對同類分子ε=2?！ウ螢榕鲎卜肿拥钠骄鶅?nèi)自由度，ωAB為碰撞對的粘性溫度指數(shù)。當(dāng)碰撞中的總碰撞能量Ec＞Ea時，反應(yīng)截面σR與總碰撞截面σT的比值就是反應(yīng)發(fā)生的概率，這種化學(xué)反應(yīng)模型也稱為TCE(Total Collision Energy)模型。

在DSMC模擬過程中為簡化碰撞算法，假設(shè)非彈性碰撞的化學(xué)反應(yīng)過程和內(nèi)能交換是相互獨立的事件。由于碰撞分子發(fā)生化學(xué)反應(yīng)的幾率在所有非彈性碰撞事件中是最低的，所以首先判斷其是否發(fā)生化學(xué)反應(yīng)。離解反應(yīng)、交換反應(yīng)等不同的化學(xué)反應(yīng)過程在滿足反應(yīng)幾率的條件下都有可能發(fā)生。在無法獲知反應(yīng)分子遵循何種反應(yīng)軌跡生成新的反應(yīng)生成物時，可以假定每種反應(yīng)過程都是獨立的。對所有可能發(fā)生化學(xué)反應(yīng)的幾率求和得到總的反應(yīng)幾率，如果Pr＞Rf，碰撞分子發(fā)生化學(xué)反應(yīng)［19］。具體發(fā)生的是何種反應(yīng)類型，則同樣根據(jù)幾率/Pr與隨機(jī)數(shù)Rf的比較來選取，然后按照相應(yīng)的化學(xué)反應(yīng)類型產(chǎn)生新的生成物(分子或原子)。

對于不發(fā)生化學(xué)反應(yīng)的碰撞分子，即Rf＞Pr時，則根據(jù)分子振動抽樣概率，判斷是否發(fā)生振動自由度的激發(fā);對于未發(fā)生振動激發(fā)的分子，再根據(jù)分子轉(zhuǎn)動抽樣概率，判斷是否發(fā)生轉(zhuǎn)動自由度的激發(fā);對于發(fā)生內(nèi)自由度激發(fā)的分子，按相應(yīng)的能量交換模型進(jìn)行能量的再分配，否則分子按彈性碰撞處理。

由于本文模擬的高度較高，僅考慮五組元空氣(O2，N2，O，N，NO)的離解反應(yīng)和交換反應(yīng)，燒蝕的影響暫不作考慮。

1．4 DSMC并行算法

DSMC并行算法采用區(qū)域分解的策略，根據(jù)計算的處理器數(shù)(或CPU核心數(shù))的多少將計算區(qū)域劃分為等量的子區(qū)域，每個處理器在其分配的子區(qū)域內(nèi)部獨立地計算模擬分子的碰撞和遷移，離開子區(qū)域的模擬分子把攜帶的信息傳遞給對應(yīng)子區(qū)域的處理器。并行計算總的時間包括每個處理器計算碰撞、遷移的時間、處理器之間的通訊時間以及各處理器之間為同步而等待的時間。提高并行計算的效率主要是通過減少通訊和同步等待的時間來實現(xiàn)。

為了減少不同處理器計算時間的差別，通常采用負(fù)載平衡技術(shù)［20－22］，本文采用靜態(tài)隨機(jī)負(fù)載平衡方法用于解決不同處理器之間的計算時間同步問題［12］，該方法基于概率近似原理，將計算區(qū)域的全部網(wǎng)格平均分配給指定的所有處理器。由于采用相同的概率隨機(jī)選取流場網(wǎng)格，按照均分后的數(shù)量分配給每個處理器，因此當(dāng)計算區(qū)域的網(wǎng)格數(shù)量較多時，每個處理器包含近似相等的物體邊界、高密度流動區(qū)域和稀薄氣體區(qū)域的網(wǎng)格數(shù)，這樣每個處理器的計算負(fù)載也非常接近。

2 計算結(jié)果分析

本文選取天宮一號目標(biāo)飛行器的簡化模型作為研究對象，飛行器上的太陽翼帆板朝向考慮兩種情況，分別是板面平行于和垂直于飛行器中心軸線。計算的來流速度為7．6 km/s，高度范圍為120 km～90 km，每隔10 km一個狀態(tài)，攻角(Angle of Attack，AOA)為0°和20°，壁面溫度設(shè)置為500 K、壁面反射模型為Maxwell模型(80%的漫反射、20%鏡面反射)，僅計算對稱的半個流場，計算中考慮了分子的轉(zhuǎn)動、振動能量激發(fā)，90 km高度時還考慮了五組元空氣的離解反應(yīng)和交換反應(yīng)。

2．1太陽翼帆板平行于目標(biāo)飛行器

圖2給出了太陽翼帆板平行于天宮一號目標(biāo)飛行器情況下，0°攻角下不同再入高度時的壓力等值線分布云圖。從圖中可以看出120 km、110 km高度的大氣密度相對較低，盡管飛行速度達(dá)到了7．6 km/s，但飛行器的頭部并未形成明顯的脫體激波，而是逐漸增強的壓縮過程，氣體在整個太陽翼帆板的側(cè)面也產(chǎn)生明顯的壓縮。圖中也顯示出，飛行器頭部簡化的對接機(jī)構(gòu)形狀對飛行器前方壓縮氣流的擾動影響，而隨著飛行高度的降低，這種影響不斷減弱，直至形成較強的頭部脫體激波。100 km和90 km的頭部激波作用在太陽翼帆板上產(chǎn)生了比駐點值高幾倍的峰值壓力和峰值熱流。

圖2 太陽翼水平放置時0°攻角不同高度壓力分布Fig．2 Pressure distribution of 0°AOA in different altitude with horizontal solar array

圖3給出了120 km和90 km在20°攻角下的壓力分布云圖。圖中顯示出，由于攻角的存在，除了在飛行器的頭部產(chǎn)生壓縮以外，整個飛行器艙體的迎風(fēng)面也會產(chǎn)生較強的壓縮，使艙體和太陽翼帆板的迎風(fēng)受力面增大，作用在太陽翼帆板上的力和力矩都有較大增加。

圖3 太陽翼水平放置時20°攻角不同高度壓力分布Fig．3 Pressure distribution of 20°AOA in different altitude with horizontal solar array

表1、表2分別給出了0°攻角和20°攻角不同高度情況下，作用在單側(cè)太陽翼帆板上的軸向力、偏航力矩和峰值壓力、峰值熱流。由于太陽翼帆板平行于軸向方向，0°攻角時只有帆板的側(cè)面為主要受力面，因此作用在帆板上面的軸向力和偏航力矩(力矩作用點取帆板與艙體的連接處)相對較小，但在100 km和90 km由于頭部激波的干擾作用，使干擾區(qū)域的峰值熱流非常高，100 km時的峰值熱流為150 kW/m2左右，而到90 km時則達(dá)到了1．2 MW/m2。20°攻角時，盡管作用在帆板上的峰值壓力和峰值熱流與0°攻角時相差不大，但作用的軸向力和偏航力矩明顯增大。因此，當(dāng)太陽翼帆板平行放置時，激波干擾引起較強的氣動加熱可能最先破壞太陽翼帆板的結(jié)構(gòu)，然后才是在氣動力的作用下將其撕毀并脫離飛行器艙體。

表1 作用在水平太陽翼上的氣動力和熱流峰值(0°攻角)Table 1 Aerodynamics and peak heat flux acting on horizontal solar array at 0°AOA

表2 作用在水平太陽翼上的氣動力和熱流峰值(20°攻角)Table 2 Aerodynamics and peak heat flux acting on horizontal solar array at 20°AOA

2．2太陽翼帆板垂直于目標(biāo)飛行器

圖4給出了太陽翼帆板垂直于目標(biāo)飛行器0°攻角、不同再入高度時的壓力等值線分布云圖。圖中顯示出的流場結(jié)構(gòu)與圖2相似，只是此時太陽翼帆板有較大的迎風(fēng)面，其前方的壓縮層更厚(120 km、110 km)，激波/激波、激波/邊界層的干擾更強、更復(fù)雜(100 km、90 km)，同時也與艙體產(chǎn)生了較強的相互作用，該區(qū)域的激波干擾使資源艙表面的壓力和熱流上升，尤其是在太陽翼帆板與艙體的連接處附近，更容易使太陽翼帆板遭到破壞，另外作用在帆板上的力和力矩明顯比水平放置時大得多，因此垂直放置的太陽翼帆板會比水平放置的太陽翼帆板更早地在氣動力、氣動熱的作用下被撕裂并脫離飛行器艙體。

圖5給出了90 km高度、0°攻角下的壓力等值線分布，可以更加詳細(xì)地分析復(fù)雜激波干擾下的流動特征。頭部脫體激波與太陽翼帆板的脫體激波相互作用并產(chǎn)生一道正激波對帆板表面形成非常強烈的壓縮，使壓力急劇上升。而在艙體表面附近區(qū)域，由于氣體的膨脹，壓力迅速下降，因此在頭部脫體激波、艙體和太陽翼帆板表面之間的空域形成較大的逆壓梯度，引起流動分離，同時產(chǎn)生的分離激波也與頭部脫體激波相互干擾。圖中還可以看出，從形成的正激波作用在帆板表面的位置開始，沿帆板左右各形成一股“噴流”，朝向艙體的“噴流”在較大的壓力梯度下加速至超聲速并直接撞擊在資源艙的表面，形成局部的脫體激波。由于在較大的區(qū)域存在著逆壓梯度，形成一個大的分離渦，而在帆板脫體激波影響的區(qū)域僅有微弱的流動分離，如圖6所示。

圖4 太陽翼垂直放置時0°攻角不同高度壓力分布Fig．4 Pressure distribution of 0°AOA in different altitude with perpendicular solar array

圖5 90 km、0°攻角的壓力等值線分布Fig．5 Pressure contours distribution of 90 km at 0°AOA

圖7、圖8分別給出了90 km高度、0°和20°攻角下的太陽翼帆板表面的壓力和熱流分布的對比，0°攻角作用在帆板表面的高壓力和高熱流區(qū)域近似為條狀，而20°攻角時則近似為弧形，并且0°攻角時作用在帆板上的壓力和熱流強度高于20°攻角情況。從表3和表4給出的作用在帆板上的力和力矩表明，兩種攻角下的力和力矩相差并不明顯。在相同高度，作用在帆板的力和力矩要比水平放置時大一個量級。盡管100 km高度帆板表面的峰值熱流比帆板水平放置時有較大降低，但由于100 km高度時太陽翼帆板所受的力和力矩是帆板水平放置時在90 km高度時的近3倍，超過了600 N和2700 N·m，到90 km時更是達(dá)到了3950 N和17000 N·m。因此太陽翼帆板垂直放置時，在100 km高度就有可能在氣動力、氣動熱的雙重作用下被撕裂并脫離飛行器本體。

圖6 90 km、0°攻角的流線分布Fig．6 Stream lines distribution of 90 km at 0°AOA

圖7 太陽翼帆板表面壓力分布Fig．7 Pressure distribution of solar panels surface

圖8 太陽翼帆板表面熱流分布Fig．8 Heat flux distribution of solar panels surface

表3 作用在垂直太陽翼上的氣動力和熱流峰值(0°攻角)Table 3 Aerodynamics and peak heat flux acting on perpendicular solar array at 0°AOA

表4 作用在垂直太陽翼上的氣動力和熱流峰值(20°攻角)Table 4 Aerodynamics and peak heat flux acting on perpendicular solar array at 20°AOA

3 結(jié)論

通過上面的計算分析，得出以下幾點初步結(jié)論:

1)100 km以上高度，來流氣體對目標(biāo)飛行器及太陽翼帆板是漸進(jìn)的壓縮過程，作用在太陽翼帆板上的力和力矩、峰值壓力和熱流相對較低。

2)100 km及以下高度，由于飛行器頭部脫體激波與太陽翼帆板產(chǎn)生比較強的干擾，帆板所受的氣動力、力矩以及表面上的峰值壓力、熱流隨高度降低明顯增大，且增大的速率加快。

3)太陽翼帆板水平放置時，激波干擾引起的側(cè)緣峰值熱流達(dá)到了兆瓦量級，較高的氣動加熱可能最先引起帆板結(jié)構(gòu)的破壞，然后在氣動力的共同作用下將其撕毀并脫離飛行器本體。

4)太陽翼帆板垂直放置時，飛行器頭部脫體激波與帆板脫體激波產(chǎn)生更強烈、更復(fù)雜的激波/激波、激波/邊界層的干擾。在相同高度，作用在帆板上的氣動力和力矩要比水平放置時大一個量級，因此太陽翼帆板垂直放置時要比水平放置時更快地在氣動力、氣動熱的雙重作用下被撕裂并脫離目標(biāo)飛行器。

［1］ 朱魯青．美國“高層大氣研究衛(wèi)星”隕落事件分析［J］．國際太空，2011，11:27－33．［Zhu Lu-qing．Analysis of UARS reentry event［J］．Space International，2011，11:27－33．］

［2］ 都亨，張文祥，龐寶軍，等．空間碎片［M］．北京:中國宇航出版社，2007:1－128．

［3］ 吳連大，王昌彬．危險目標(biāo)再入跟蹤預(yù)報研究［J］．紫金山天文臺臺刊，2000，19(2):100－105．［Wu Lian-da，Wang Chang-bin．Research on tracking and prediction of risk object reentry［J］．Publications of Purple Mountain Observatory，2000，19(2):100－105．］

［4］Opiela JN，Hillary E，Whitlock D O，et al．DASuser's guide［R］．Houston:Orbital Debris Program Office in NASA Johnson Space Center，TX，2007

［5］Bouslog SA，Ross B P，Madden C B．Space debris reentry risk analysis［R］．AIAA 1994－591，1994．

［6］Fritsche B，Klinkrad H，Kashkovsky A，et al．Spacecraft disintegration during uncontrolled atmospheric re-entry［J］．Acta Astronaut，2000，47(2－9):513－522．

［7］Wu Z N，Hu R F，Qu X，et al．Space debris reentry analysis method and tools［J］．Chin J．Aeronaut，2011，24:387－395．

［8］Li JL，Yan SZ．Thermally induced vibration of composite solar array with honeycomb panels in low earth orbit［J］．Applied Thermal Engineering，2014，71:419－432．

［9］Li JL，Yan SZ，Cai R Y．Thermal analysis of composite solar array subjected to space heat flux［J］．Aerospace Science and Technology，2013，27(1):84－94．

［10］Bird G A．Molecular gas dynamics and the direct simulation of gas flows［M］．Oxford:Clarendon Press，1994．

［11］ 黃飛，趙波，程曉麗，等．返回器高空稀薄氣動特性的真實氣體效應(yīng)研究［J］．宇航學(xué)報，2014，35(3):283－290．［Huang Fei，Zhao Bo，Cheng Xiao-li，et al．Real gas effects of reentry vehicle aerodynamics under hypersonic［J］．Journal of Astronautics，2014，35(3):283－290．］

［12］ 梁杰，閻超，楊彥廣，等．過渡區(qū)側(cè)向噴流干擾的并行DSMC數(shù)值模擬研究［J］．宇航學(xué)報，2011，32(5):1012－1018．［Liang Jie，Yan Chao，Yang Yan-guang，et al．Parallel DSMC simulation of lateral jet interaction in rarefied transitional region［J］．Journal of Astronautics，2011，32(5):1012－1018．］

［13］ 梁杰，閻超，杜波強．基于兩級直角網(wǎng)格結(jié)構(gòu)的三維DSMC算法研究［J］．空氣動力學(xué)學(xué)報，2010，28(4):466－471．［Liang Jie，Yan Chao，Du Bo-qiang．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．］

［14］ 梁杰，閻超，李志輝，等．稀薄過渡流區(qū)橫向噴流干擾效應(yīng)數(shù)值模擬研究［J］．空氣動力學(xué)學(xué)報，2013，31(1):27－33．［Liang Jie，Yan Chao，Li Zhi-hui，et al．Numerical investigation of lateral jet interaction effects in rarefied transition flow regime［J］．Acta Aerodynamica Sinica，2013，31(1):27－33．］

［15］Lumpkin F E，Hass B L，Boyd I D．Resolution of differences between collision number definition in particle and continuum simulations［J］．Physics of Fluids，1991，A3:2282－2284．

［16］Parker J G．Rotational and vibrational relaxation in diatomic gases［J］．Physics of Fluids，1959，2(4):449－462．

［17］Millikan R C，White D R．Systematics of vibrational relaxation［J］．The Journal of Chemical Physics，1963，39(12):3209－3213．

［18］Park C．Problems of rate chemistry in the flight regimes of aeroassisted orbital transfer vehicles［C］．Progress in Astronautics and Aeronautics，AIAA，New York，USA，Jule9－14，1985:511－537．

［19］Gimelshein N E，Gimelshein S F，Levin D A，et al．Reconsideration of DSMC models for internal energy transfer and chemical reactions［C］．23rd International Symposium of Rarefied Gas Dynamics，Whistler，Canada，July 20－25，2002:349－357．

［20］Ivanov M，Markelov G，Taylor S，et al．Parallel DSMCstrategies for 3D computation［C］．Proc．Parallel CFD’96，Capri，Italy，July 13－16，1996:485－492．

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

［22］Wu J S，Tseng K C．Parallel DSMC method using dynamic domain decomposition［J］．Int．J．Numer．Meth．Engng．2005，63:37–76．