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

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

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

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

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

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

高超聲速；分離;分離氣動(dòng)特性;Monte Carlo打靶仿真

0 引言

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

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

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

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

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

1 Monte Carlo打靶分析方法

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

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

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

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

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

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

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

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

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

Xi+1=aXimodm

δ=ri=Xi/m

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

由偽隨機(jī)數(shù)生成滿足各種分布的隨機(jī)數(shù)的過(guò)程稱為隨機(jī)變量的抽樣，本文采用正態(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 氣動(dòng)特性的修正

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

Fmodify=F+δ·Δ·F

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

1.4 碰撞判據(jù)

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

δ<(δ1+δ2)

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

2 仿真結(jié)果

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

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

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

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

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

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

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

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

3 結(jié)論

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

[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] 王元靖, 吳繼飛, 陶洋,等. 高超聲速多體干擾與分離試驗(yàn)[J].航空動(dòng)力學(xué)報(bào),2010,25(4):902-906.

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

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

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

[7] 王巍，劉君，白曉征，等. 非結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)及其在超聲速飛行器頭罩分離模擬中的應(yīng)用[J].空氣動(dòng)力學(xué)報(bào),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] 鄧小剛，宗文剛，張來(lái)平，等.計(jì)算流體力學(xué)中的驗(yàn)證與確認(rèn)[J]. 力學(xué)進(jìn)展,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] 王智慧, 鮑麟, 童秉綱. 高超聲速尖頭體駐點(diǎn)熱流從連續(xù)態(tài)過(guò)渡到稀薄態(tài)的變化特征和橋函數(shù)研究[J]. 中國(guó)科學(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ì)算方法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2012，30(4):508-513.

[16] 李志輝，蔣新宇，吳俊林，等. 求解Boltzmann模型方程高性能并行算法在航天跨流域空氣動(dòng)力學(xué)應(yīng)用研究[J]. 計(jì)算機(jī)學(xué)報(bào), 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-)，男，博士，研究員，主要從事高超聲速空氣動(dòng)力學(xué)和氣動(dòng)布局方面研究。E-mail:dutao_calt@yahoo.com

V417

A

2096-4080(2017)01-0054-06