杜 濤， 沈 丹，李凰立，蘇 虹，徐珊姝

(北京宇航系統(tǒng)工程研究所，北京 100076)

高超聲速多體分離的Monte Carlo模擬仿真分析研究

杜 濤， 沈 丹，李凰立，蘇 虹，徐珊姝

(北京宇航系統(tǒng)工程研究所，北京 100076)

在臨近空間飛行器飛行試驗任務(wù)中，運載器與飛行器之間的分離多數(shù)在臨近空間發(fā)生。飛行器與箭體之間的氣動特性，對分離過程的安全性有重要影響。然而，分離過程的動態(tài)特性和飛行器外形的復(fù)雜性給氣動特性的正確預(yù)測造成了巨大困難，地面試驗也很難驗證。為此，發(fā)展了Monte Carlo打靶分析方法，在誤差帶寬內(nèi)隨機生成氣動特性，通過增大氣動特性的覆蓋性，大大提高了氣動特性對分離過程的安全性影響評估的正確性和可靠性。

高超聲速；分離;分離氣動特性;Monte Carlo打靶仿真

0 引言

傳統(tǒng)宇航飛行任務(wù)中，在嚴(yán)格的大氣動壓限制條件下，頭體分離、整流罩分離等過程選擇在大氣稀薄區(qū)域，氣動特性對分離的影響忽略不計。隨著高超聲速技術(shù)的發(fā)展，諸多高超聲速飛行試驗任務(wù)要求在大氣層內(nèi)高速飛行條件下完成分離過程，例如X-43的飛行器與運載器分離高度為28km，而分離速度達(dá)到7馬赫[1]。大氣層內(nèi)的分離問題，最大的困難在于分離氣動特性的復(fù)雜性，難以預(yù)測，因為分離對象的氣動特性與分離過程的狀態(tài)高度耦合，互相作用，分離狀態(tài)受分離氣動特性影響。而分離氣動特性除與來流和外形相關(guān)外，還與分離體之間的幾何位置和運動狀態(tài)相關(guān)。要正確評估飛行器與運載器之間分離的安全性在飛行器設(shè)計上是一件非常困難的事情[2]。

預(yù)測分離過程氣動特性，目前主要依靠風(fēng)洞試驗和數(shù)值仿真。模擬分離過程的風(fēng)洞試驗，包括軌跡捕獲法(CTS)或者風(fēng)洞自由飛技術(shù)[3-4]。前者根據(jù)分離彈道模擬結(jié)果確定飛行體之間的幾何位置關(guān)系，確定模型在風(fēng)洞中的位置關(guān)系，從而測量獲得分離體的氣動特性。該方法的缺點是采用了準(zhǔn)定常假設(shè)，只能模擬分離體之間的靜態(tài)幾何關(guān)系，無法模擬分離過程的動態(tài)效應(yīng)，同時模型支持方式對分離過程的氣動特性也有影響。風(fēng)洞自由飛能夠克服軌跡捕獲方法的缺點[5]，模型在風(fēng)洞試驗過程無約束，但是自由飛行的模型尺度一般比較小,需要同時準(zhǔn)確模擬幾何外形和質(zhì)量特性，加工困難，而且是一次性，成本比較高，同時分離的初始條件的準(zhǔn)確實現(xiàn)比較困難[6]。

隨著數(shù)值模擬技術(shù)和計算機能力的提升，數(shù)值仿真技術(shù)在飛行器設(shè)計上發(fā)揮著越來越大的作用[7]。采用算流體動力學(xué)(CFD)和彈道方程耦合的方式可以開展頭體分離過程的分析，模擬動態(tài)分離過程，甚至考慮控制律[8]。但是數(shù)值仿真一次只能對單一工況進行模擬，動態(tài)模擬耗時和成本仍然比較高。最大困難是目前缺乏對數(shù)值計算本身可信度的驗證(verification)和確認(rèn)(validation)，給飛行器設(shè)計工作帶來了風(fēng)險[9-10]。

臨近空間飛行器的分離問題，還有自身的特殊性，高馬赫數(shù)帶來了化學(xué)反應(yīng)、非平衡效應(yīng)對氣動特性的影響[11]，跨流域問題涉及流動是連續(xù)流、過渡流還是稀薄流動的復(fù)雜性[12-13]，即使采用CFD方法或者地面試驗手段開展臨近空間飛行的飛行器的靜態(tài)氣動特性的模擬，無論理論研究，還是開展工程研制都是困難和有挑戰(zhàn)性的問題[14]。同時分離過程的動態(tài)特性，采用準(zhǔn)定常假設(shè)，可能會導(dǎo)致計算不準(zhǔn)確，甚至計算錯誤。因此，如何正確評估出氣動特性對臨近空間分離的影響顯得非常重要。

綜上所述，目前并沒有一種方法可以為大氣層內(nèi)高速分離問題的評估提供準(zhǔn)確的模擬。設(shè)計上更關(guān)心的不是對分離過程的正確模擬，更關(guān)注分離的安全裕度和魯棒性，在允許的偏差范圍內(nèi)，分離是否仍然安全。Monte Carlo打靶分析方法是一種有效的方法，可以克服上述方法的不足，給分離的安全性評估提供有價值的結(jié)論。Monte Carlo打靶方法的基本原理是在已有分離氣動數(shù)據(jù)庫的基礎(chǔ)上，通過引進一個滿足某種分布的隨機變量，產(chǎn)生某些在允許偏差范圍內(nèi)變動的隨機參數(shù)，作為單次分析的輸入量，通過大量的分析，從統(tǒng)計的角度評估這些參數(shù)變動對分離過程的影響。Monte Carlo分析方法的優(yōu)點是計算量小，結(jié)果覆蓋性強，可以分離出多種因素對分離結(jié)果的影響。Tartabini等[14]用Monte Carlo打靶分析方法分析了飛行姿態(tài)和飛行參數(shù)相關(guān)的分離條件在3σ范圍內(nèi)變化對X-43A在馬赫數(shù)為10的條件下分離的安全性評估，陳波等[15]用Monte Carlo方法分析了正態(tài)分布的隨機量來模擬彈射力、導(dǎo)彈質(zhì)量、慣矩、分離高度等偏差引起的分離邊界。本文將利用Monte Carlo打靶分析方法開展氣動特性偏差對頭體分離的安全性開展評估，能夠以比較少的計算量，通過覆蓋偏差條件，實現(xiàn)對頭體分離安全性的全面和準(zhǔn)確的評估。

1 Monte Carlo打靶分析方法

分析由初始分離條件開始，根據(jù)飛行器在每一時刻的位置和飛行器與箭體之間的相對幾何關(guān)系，從氣動特性數(shù)據(jù)庫中插值確定當(dāng)前時刻的氣動力。然后，隨機地在誤差范圍內(nèi)確定誤差值，從而得到包含誤差的氣動力的特性。最后由六自由度運動方程確定下一時刻的位置，以此推進，直到頭體距離足夠安全，仿真分析結(jié)束。分離過程中的每一個步驟都要對分離的安全性進行判斷，一旦發(fā)現(xiàn)飛行器與箭體發(fā)生相撞，宣告分離失敗，中止仿真。上述分析仿真步驟通常反復(fù)進行上千次，每次仿真的氣動特性都由隨機變量獲取，以保證仿真結(jié)果具有充分的可靠性和覆蓋性。

1.1 標(biāo)準(zhǔn)氣動特性數(shù)據(jù)庫

參考CFD數(shù)值仿真與彈道方程耦合仿真的結(jié)果，確定分離體之間的幾何位置狀態(tài)的網(wǎng)格化狀態(tài)，完整覆蓋耦合仿真結(jié)果。通過CFD數(shù)值仿真計算產(chǎn)生上述狀態(tài)下的準(zhǔn)定常氣動力特性，構(gòu)建出一個完整的氣動特性數(shù)據(jù)庫。氣動力數(shù)據(jù)庫要求能夠覆蓋所有分離過程中可能的分離姿態(tài)。

分離過程持續(xù)時間很短，耦合分析表明，整個分離過程中飛行器速度變化不超過3m/s，因此盡管這個速度變化對分離過程很重要，但是對氣動特性變化并不重要，在氣動特性預(yù)測上可以不考慮飛行器速度變化的影響。同時，假定飛行器和箭體分離發(fā)生在一個平面上，暫不考慮橫側(cè)向氣動特性，以簡化問題。耦合仿真表明，分離過程的法向位置變化不到0.2m，上述假設(shè)是合理的。

氣動數(shù)據(jù)庫的自變量為：兩個物體之間的距離dL，飛行器姿態(tài)角α1，運載器姿態(tài)角α2。

氣動特性數(shù)據(jù)庫是一個三維數(shù)據(jù)庫，每個維度包含6個狀態(tài)，狀態(tài)根據(jù)耦合仿真結(jié)果確定，數(shù)據(jù)庫共計由216個狀態(tài)組成。采用N-S方程求解準(zhǔn)定常狀態(tài)獲得，精度較高，狀態(tài)如表1所示。已知分離體的位置及姿態(tài)后，采用三線性插值，就可以得到該時刻的氣動力特性。

表1 計算狀態(tài)參數(shù)

1.2 偽隨機數(shù)的生成

素數(shù)模乘同余發(fā)生器獲得滿足均勻分布的偽隨機數(shù)：

Xi+1=aXimodm

δ=ri=Xi/m

ri為服從均勻分布的偽隨機數(shù)。

由偽隨機數(shù)生成滿足各種分布的隨機數(shù)的過程稱為隨機變量的抽樣，本文采用正態(tài)分布的抽樣，公式為：

1)產(chǎn)生r～U(0,1);

c0=2.515517,c1=0.802853,c2=0.010328,

d1=1.432788,d2=0.189269,d3=0.001308,

產(chǎn)生的X～N(0，1)。

4)令Y=σX+μ,則Y服從均值為μ,方差為σ的正態(tài)分布。

1.3 氣動特性的修正

通過氣動數(shù)據(jù)庫插值得到的是標(biāo)準(zhǔn)的氣動特性，在每次仿真中，借助偽隨機數(shù)，以氣動偏差的形式對標(biāo)準(zhǔn)氣動特性進行修正。修正公式為：

Fmodify=F+δ·Δ·F

其中，F(xiàn)是氣動數(shù)據(jù)庫插值得到的氣動力(或力矩)，δ是服從正態(tài)分布的隨機數(shù)，Δ是氣動力誤差寬帶或力矩誤差寬帶。在本文中，氣動力特性誤差帶寬設(shè)定為50%，俯仰力矩誤差帶寬設(shè)定為100%。單次的分離過程仿真過程中偽隨機數(shù)δ保持不變。通過監(jiān)測子樣的分布情況，可以確定分離仿真計算子樣的數(shù)量足夠多，保證隨機數(shù)的取值能夠覆蓋整個誤差帶寬范圍。經(jīng)過分析，3000次仿真能夠滿足覆蓋性要求。修正后的氣動特性參與飛行器的動力學(xué)計算。

1.4 碰撞判據(jù)

每一步迭代計算后，需要判斷飛行器與箭體之間是否發(fā)生碰撞，一旦飛行器與運載器之間在空間中發(fā)生重疊或相交，即宣布本次仿真分離過程失敗，計算中止。采用一個更加偏嚴(yán)格的判據(jù)，對箭體和飛行器是否發(fā)生碰撞進行識別。當(dāng)滿足如下關(guān)系時候，即宣告兩者發(fā)生碰撞，仿真中止：

δ<(δ1+δ2)

其中，δ表示飛行器質(zhì)心與箭體質(zhì)心之間的距離，δ1表示飛行器質(zhì)心到后端面的距離，δ2表示箭體質(zhì)心到前端面的距離。當(dāng)飛行器質(zhì)心與箭體質(zhì)心之間的距離大于δ1和δ2之和時，飛行器是安全的，不會與箭體發(fā)生碰撞。

2 仿真結(jié)果

本文對兩個例子進行了飛行器和箭體分離的Monte Carlo打靶分析仿真分析，每一例子共計進行了3000次打靶仿真。兩個例子均是馬赫數(shù)為23，不同的是第一個例子的飛行高度為70km，攻角為0°，第二個例子的飛行高度為80km，初始攻角為10°。下面分別給出兩個例子的分析結(jié)果。

2.1 分離高度為70km和初始攻角為0°

圖2給出了飛行器質(zhì)心位置變化，圖中每一條線代表一次仿真結(jié)果。圖中標(biāo)識的兩條線分別表示采用CFD與動力學(xué)耦合仿真在相同高度和攻角下馬赫數(shù)分別為20和25的仿真結(jié)果，后面圖中表示相同。圖3～圖5依次給出了運載器質(zhì)心位置變化，飛行器攻角變化和運載器攻角變化歷程。

圖6給出了飛行器與運載器之間距離的變化歷程。經(jīng)過統(tǒng)計，發(fā)現(xiàn)0.3%的分離過程失敗。因此，70km高度條件下的分離過程，存在分離過程發(fā)生碰撞的可能性。但是發(fā)生碰撞的概率非常低，采取一些措施，能夠消除發(fā)生碰撞的可能性。

2.2 分離高度為80km和初始攻角為10°

飛行器整體尺寸為米級，在80km高空高馬赫數(shù)飛行的Kn數(shù)比較小，整體處于連續(xù)流狀態(tài)，但是局部區(qū)域已經(jīng)有過渡流動特征，相關(guān)研究可以參見文獻(xiàn)[16]。

圖7給出了飛行器質(zhì)心位置變化，圖中每一條線代表一次仿真結(jié)果。圖中標(biāo)識的兩條線分別表示采用CFD與動力學(xué)耦合仿真在相同高度和攻角下馬赫數(shù)分別為20和25的仿真結(jié)果，后面圖中表示相同。圖8～圖10依次給出了運載器質(zhì)心位置變化，飛行器攻角變化和運載器攻角變化歷程。

圖6給出了飛行器與運載器之間距離的變化歷程。經(jīng)過統(tǒng)計，發(fā)現(xiàn)所有的分離過程都不會發(fā)生碰撞。因此，80km高度條件下的分離過程，是比較安全的，即使在氣動特性的偏差條件下，也不會發(fā)生由于氣動特性引起的碰撞的可能性。

3 結(jié)論

采用Monte Carlo打靶方法評估了分離氣動特性對臨近空間飛行器與運載器在大氣層內(nèi)分離的安全性。研究表明，Monte Carlo打靶方法中引入了誤差的隨機分布，可以有效識別不確定性因素對飛行安全性的影響，從而提高分析結(jié)果的可靠性和工程可用性，并指導(dǎo)工程設(shè)計的改進。在無法完全評估偏差因素對分離過程的影響下，Monte Carlo打靶方法在飛行器設(shè)計上可以提供覆蓋性的評估結(jié)果。該方法還可用在多種因素的影響評估上發(fā)揮作用。在后續(xù)研究中，將進一步利用Monte Carlo打靶方法分離單一因素對分離過程的影響評估，開展多重獨立不確定性干擾源作用下，對分離過程的安全性影響研究。

[1] Joyce P J，Pomroy J B，Grindle L. The Hyper-X launch vehicle: challenges and design considerations for hypersonic flight testing[R].AIAA 2005-3333，2005.

[2] Bermdez L M, Gladden R D, Jeffries M S, et al . Aerodynamic characterization of the Hyper-X launch vehicle[R].AIAA 2003-7074, 2003.

[3] 王元靖, 吳繼飛, 陶洋,等. 高超聲速多體干擾與分離試驗[J].航空動力學(xué)報,2010,25(4):902-906.

[4] 王飛，祝汝松，張俊生.卡爾曼濾波在某跨聲速風(fēng)洞CTS機構(gòu)連續(xù)動態(tài)軌跡捕獲試驗技術(shù)中的應(yīng)用[J].實驗流體力學(xué)，2015,29(1):103-108.

[5] 楊益農(nóng), 賈區(qū)耀. 多體分離拋撒初條件與分離特征參數(shù)[J].宇航學(xué)報,2005, 26(2):140-142.

[6] 賈區(qū)耀, 楊益農(nóng), 蔣增輝.風(fēng)洞自由飛實驗結(jié)果的精度、準(zhǔn)度[J].宇航學(xué)報,2005, 30(6):2082-2085.

[7] 王巍，劉君，白曉征，等. 非結(jié)構(gòu)動網(wǎng)格技術(shù)及其在超聲速飛行器頭罩分離模擬中的應(yīng)用[J].空氣動力學(xué)報,2008, 26(1):131-135.

[8] Clifton J D, Ratcliff C J, Bodkin D J, et al. Determining the stability and control characteristics of high-performance maneuvering aircraft using high-resolution CFD simulation with and without moving control surfaces[R]. AIAA 2013-0972, 2013.

[9] 鄧小剛，宗文剛，張來平，等.計算流體力學(xué)中的驗證與確認(rèn)[J]. 力學(xué)進展,2007, 35(2):279-288.

[10] American Institute of Aeronautics and Astronautics Staf.AIAA guide for the verification and validation of computational fluid dynamics simulations[C].AIAA G-007-1998,1998.

[11] James N M， Vincent C J R , Simmonds A L. Nonequilibrium effects for hypersonic transitional flows[R]. AIAA 1987-0404, 1987.

[12] 王智慧, 鮑麟, 童秉綱. 高超聲速尖頭體駐點熱流從連續(xù)態(tài)過渡到稀薄態(tài)的變化特征和橋函數(shù)研究[J]. 中國科學(xué) G, 2009，39(8):1134-1140.

[13] Wang W L, Boyd I D. Continuum breakdown in hypersonic viscous flows[R].AIAA 2002-0651, 2002.

[14] Tartabini P V, Bose D M, Thornblom M N,et al. Mach 10 stage separation analysis for the X43-A[C].AIAA 2006-1038, 2006.

[15] 陳波，劉剛，肖涵山，等.基于蒙特卡羅模擬的導(dǎo)彈分離邊界計算方法研究[J].空氣動力學(xué)學(xué)報，2012，30(4):508-513.

[16] 李志輝，蔣新宇，吳俊林，等. 求解Boltzmann模型方程高性能并行算法在航天跨流域空氣動力學(xué)應(yīng)用研究[J]. 計算機學(xué)報, 2016, 39(9):1801-1811.

Research on Monte Carlo Simulation for HypersonicMulti-body Separation

DU Tao, SHEN Dan, LI Huang-li, SU Hong, XU Shan-shu

(Beijing Institute of Astronautical Systems Engineering,Beijing 100076,China)

In the hypersonic vehicle flight experiment, separation between hypersonic vehicle and launch vehicle will take place in the near space. The separation aerodynamic characteristic plays an important role on security in separation process. However, because of complex configuration of vehicles and dynamic characteristic for separation process, the aerodynamic characteristic for separation process is difficult to be simulated with numerical method or in wind tunnel experiment. The Monte Carlo simulation method is developed to improve the validity and reliability for evaluating the separation process through stochastic aerodynamic characteristic being incumbent on error bound.

Hypersonic; Separation; Separation aerodynamic characterisitic; Monte Carlo simulation

2017-01-13；

2017-04-24

杜濤(1973-)，男，博士，研究員，主要從事高超聲速空氣動力學(xué)和氣動布局方面研究。E-mail:dutao_calt@yahoo.com

V417

A

2096-4080(2017)01-0054-06