蔡國飆 賀碧蛟

（北京航空航天大學(xué)宇航學(xué)院）

1 引言

飛船、空間站等大型航天器具有質(zhì)量大、艙段多的特點，其上需要布置多塊太陽能帆板才能滿足航天器對電力的需求。由于帆板數(shù)量較多，對與部分發(fā)動機(jī)距離較近的帆板，發(fā)動機(jī)真空羽流會對其產(chǎn)生較大的氣動力和力矩作用。為了對發(fā)動機(jī)真空羽流效應(yīng)進(jìn)行數(shù)值模擬，必須針對不同的流動狀態(tài)建立相應(yīng)的數(shù)學(xué)物理模型，并采用與之相應(yīng)的數(shù)值模擬方法。目前能夠比較成功地模擬真空羽流效應(yīng)的研究方法是直接模擬蒙特卡羅方法[1]（DSMC）。

國外針對發(fā)動機(jī)羽流問題的研究開展得比較早，開發(fā)了一系列基于DSMC的模擬計算軟件，其中比較典型的通用計算軟件有：由美國約翰遜航天中心（Johnson Space Center）的 LeBeau[2-4]等人開發(fā)的DAC（DSMC Analysis Code）軟件，由俄羅斯理論與應(yīng)用機(jī)械研究所（the Institute of Theoretical and Applied Mechanics）的 Ivanov[5]等人開發(fā)的 SMILE（Statistical Modeling In Low-density Environment）軟件，由美國Cornell大學(xué)的Dietrich和Boyd等人開發(fā)的MONACO 計算軟件[6，7]等。

文章所介紹的基于 DSMC的 PWS（Plume WorkStation）羽流計算軟件是依據(jù)模塊化的原則，適當(dāng)引入面向?qū)ο蟮木幊趟枷?，特別設(shè)計具有較強(qiáng)通用性的粒子邊界處理模塊，實現(xiàn)PWS軟件在羽流模擬計算的通用化，并在國外“CUBRC”試驗研究的基礎(chǔ)上開展試驗驗證工作。使用PWS軟件就神舟飛船120N平移發(fā)動機(jī)羽流對太陽能帆板氣動力效應(yīng)進(jìn)行了數(shù)值模擬研究，并且對比了由不同帆板角度對作用在帆板中心點的力和力矩影響。

2 羽流數(shù)值模擬模型和方法

流場依據(jù)Knudsen（Kn）數(shù)劃分為三大領(lǐng)域，即連續(xù)流區(qū)領(lǐng)域（Kn＜0.1）、過渡領(lǐng)域（0.1＜Kn＜10）和自由分子流領(lǐng)域（Kn＞10）[8]，Kn= λ/L，其中 λ 是平均分子自由程，L為流動特征長度。在Kn數(shù)較大，即氣體變稀薄時，基于連續(xù)流介質(zhì)的N-S方程失效。

真空羽流流場包括了全部流態(tài)，即連續(xù)介質(zhì)流、過渡領(lǐng)域流和自由分子流等[1]。羽流核心區(qū)為連續(xù)流，外圍為過渡領(lǐng)域流和自由分子流，羽流流動極為復(fù)雜，很難用一種方法來描述。對于連續(xù)流動，對應(yīng)的數(shù)學(xué)方程是N-S方程；對于自由分子流，對應(yīng)的數(shù)學(xué)方程是無碰撞項的Boltzmann方程；而對于過渡領(lǐng)域流的研究，目前尚無完善的理論，與之相對應(yīng)的數(shù)學(xué)方程是完全的Boltzmann方程，它是一個非常復(fù)雜的非線性積分-微分方程，除了極其簡單的情形，它的理論求解極為困難。

為了對羽流流場進(jìn)行數(shù)值模擬，必須針對不同的流動狀態(tài)建立相應(yīng)的數(shù)學(xué)物理模型，并采用與之相應(yīng)的數(shù)值模擬方法。目前能夠比較成功地模擬自由分子流和過渡領(lǐng)域流的數(shù)值方法是DSMC。DSMC方法是20世紀(jì)60年代初由G.A.Bird提出并發(fā)展的直接模擬方法，在數(shù)學(xué)上已被證明與Boltzmann方程的相容。

G.A.Bird[1]提出的DSMC直接從物理實際出發(fā)，利用少量的模擬分子代替真實流場內(nèi)數(shù)目眾多的氣體分子，用計算模擬代替物理過程，經(jīng)統(tǒng)計平均獲得宏觀流動參數(shù)，達(dá)到求解稀薄氣體問題的目的。典型的DSMC計算流程圖如圖1所示。

圖1 DSMC計算流程圖

由于羽流場DSMC計算域的密度值較低，分析時常作如下假設(shè)：（1）流場中分子的碰撞為二體碰撞；（2）由于氣體溫度不高，僅考慮分子的轉(zhuǎn)動內(nèi)能，忽略分子的振動內(nèi)能和分子化學(xué)非平衡效應(yīng)；（3）氣體流動為定常流動；（4）將分子視為變徑硬球分子（VHS）。

3 PWS軟件及驗證

DSMC使用大量的“模擬分子”以模擬真實的氣體，為此需要較大的存儲空間和較高計算速度，這極大地限制了DSMC的應(yīng)用。DSMC的計算主程序主要有6個主模塊：粒子模塊（粒子參數(shù)存儲）、網(wǎng)格模塊、粒子進(jìn)入模塊、粒子運(yùn)動模塊、邊界模塊（主要指固體邊界）和粒子碰撞模塊?？紤]到計算速度的需要，在6個主模塊的基礎(chǔ)上再增加一個并行模塊。在羽流效應(yīng)計算時，DSMC的各模塊信息交換如圖2所示。

圖2 DSMC計算模塊信息交換關(guān)系圖

根據(jù)模塊化和面向?qū)ο蟮能浖O(shè)計方法，PWS軟件對DSMC計算程序進(jìn)行如下改進(jìn)：

（1）使用松散的軟件架構(gòu)保證各個計算模塊的相對獨立性，即改進(jìn)其中一個模塊時，基本上不需要改動其他模塊；

（2）對各個模塊的功能實現(xiàn)采用通用性設(shè)計，即在使用單/多組分分子、不同分子模型（碰撞模型和內(nèi)能模型）、不同固體邊界條件時不需要改動源代碼。

如圖2所示，基于DSMC的PWS羽流計算通用軟件主要包括網(wǎng)格模塊、粒子參數(shù)模塊、粒子進(jìn)入模塊、粒子運(yùn)動及邊界模塊、粒子碰撞模塊、并行模塊模塊6大模塊。在這6大基本模塊的基礎(chǔ)上，針對各種實際問題（單組分、多組分、多種碰撞模型、多種邊界條件），PWS羽流計算軟件可以方便地搭建數(shù)值模擬應(yīng)用平臺。

軟件計算結(jié)果的驗證是軟件編制不可缺少的一環(huán)。本文使用了美國Holden[9]等人做的“CUBRC”第7次試車數(shù)據(jù)作為軟件計算驗證的對比數(shù)據(jù)?！癈UBRC”實驗是一個超音速稀薄氣流對雙圓錐體的氣動力和氣動加熱效應(yīng)的測量試驗，試驗用的雙圓錐體構(gòu)型如圖3所示。圖4為利用PWS軟件對雙圓錐體表面壓強(qiáng)分布的計算結(jié)果與試驗測量數(shù)據(jù)對比圖，圖5為利用PWS軟件對雙圓錐體表面熱流密度分布的計算結(jié)果與試驗測量數(shù)據(jù)對比圖，軸向距離X指從雙圓錐體的頂點為原點，計算點在雙圓錐體的中心線上投影的坐標(biāo)值（m）。

圖3 “CUBRC”試驗雙圓錐體結(jié)構(gòu)圖（單位：mm）

圖4 雙圓錐體表面壓強(qiáng)分布計算值與試驗值對比

圖5 雙圓錐體熱流密度分布計算值與試驗值對比

從圖4、圖5可以看出雙圓錐體表面壓強(qiáng)和熱流密度計算值與試驗值相符，最大值和最小值的位置一致；除了2個點之外，其他各點的壓強(qiáng)計算值相對試驗值的偏差都在±20%以內(nèi)；熱流密度計算值相對試驗值的偏差都在-30%～+10%范圍以內(nèi)，特別是在雙圓錐連接處附近符合得相當(dāng)好。經(jīng)過與實驗值比對，可以認(rèn)為PWS計算軟件具有較高的置信度。

4 太陽能帆板羽流效應(yīng)計算

神舟飛船平移發(fā)動機(jī)推力相對較小，故可直接將發(fā)動機(jī)出口截面作為DSMC計算的模擬分子入口條件，使用PWS軟件直接對發(fā)動機(jī)三維羽流場進(jìn)行數(shù)值模擬，得到發(fā)動機(jī)單機(jī)/雙機(jī)工作時羽流對太陽能帆板在0°、45°和90°三種狀態(tài)下產(chǎn)生的氣動力效應(yīng)值，具體計算算例如表1所示。發(fā)動機(jī)與太陽能帆板的相對位置如圖6所示。

表1 對比算例

圖6 發(fā)動機(jī)與太陽能帆板相對位置

4.1 計算條件

（1）粒子入口條件

粒子入口條件為：以發(fā)動機(jī)出口截面作為PWS軟件計算外噴流和后流區(qū)的粒子入口截面條件；“模擬粒子”選取燃燒產(chǎn)物主要組分 N2、H2O、CO、CO2和H2等5種氣體分子，質(zhì)量分?jǐn)?shù)如表2所示。

表2 發(fā)動機(jī)燃燒產(chǎn)物主要組分質(zhì)量分?jǐn)?shù)

粒子入口截面物理相對值分布如圖7所示，圖中各參考量的的定義為：參考密度ρ0=0.005kg/m3，參考壓強(qiáng)P0=2000Pa參考溫度T0=3000K，參考軸向速度u0=3500m/s，參考徑向速度v0=500m/s，發(fā)動機(jī)出口半徑R0=0.039m。

圖7 粒子入口物理參數(shù)分布圖

（2）邊界條件

真空邊界：粒子經(jīng)過真空邊界后逃逸（即刪除）；

固體邊界：粒子撞擊到固體邊界按照Maxwell反射[1]處理，Maxwell反射需確定固體邊界壁面溫度和熱適應(yīng)系數(shù)，發(fā)動機(jī)壁面溫度設(shè)定為800 K，飛船和太陽能帆板表面溫度設(shè)定為300 K，熱適應(yīng)系數(shù)統(tǒng)一設(shè)定為1。

4.2 計算結(jié)果與分析

圖8至圖 13為算例1至算例6太陽能帆板所受的發(fā)動機(jī)羽流氣動壓強(qiáng)分布圖，圖中P指的是太陽能帆板表面正壓強(qiáng)（Pa）。

圖8 算例1太陽能帆板受到的羽流氣動壓強(qiáng)分布圖

圖9 算例2太陽能帆板受到的羽流氣動壓強(qiáng)分布圖

圖10 算例3太陽能帆板受到的羽流氣動壓強(qiáng)分布圖

圖11 算例4太陽能帆板受到的羽流氣動壓強(qiáng)分布圖

圖12 算例5太陽能帆板受到的羽流氣動壓強(qiáng)分布圖

從圖8至圖 13可以看出T1和T2平移發(fā)動機(jī)羽流對太陽能帆板的氣動壓強(qiáng)分布：T2發(fā)動機(jī)單機(jī)工作時羽流對0°狀態(tài)太陽能帆板氣動壓強(qiáng)最大值達(dá)到12Pa，且壓強(qiáng)值在1Pa以上的作用面積較小；T2發(fā)動機(jī)單機(jī)工作時羽流對45°狀態(tài)太陽能帆板氣動壓強(qiáng)最大值達(dá)到10Pa，且壓強(qiáng)值在1Pa以上的作用面積較大；T2發(fā)動機(jī)單機(jī)工作時羽流對90°狀態(tài)太陽能帆板氣動壓強(qiáng)最大值僅達(dá)到2.4Pa，但壓強(qiáng)值在1Pa以上的作用面積較大；T1和T2發(fā)動機(jī)雙機(jī)工作時羽流對0°狀態(tài)太陽能帆板氣動壓強(qiáng)最大值依然為12Pa，且壓強(qiáng)值在1Pa以上的作用面積與單機(jī)工作相比有所增加，但還是較??；T1和T2發(fā)動機(jī)雙機(jī)工作時羽流對45°狀態(tài)太陽能帆板氣動壓強(qiáng)最大值達(dá)到10Pa，且壓強(qiáng)值在1Pa以上的作用面積很大；T1和T2發(fā)動機(jī)雙機(jī)工作時羽流對90°狀態(tài)太陽能帆板氣動壓強(qiáng)最大值達(dá)到7.5Pa，且壓強(qiáng)值在1Pa以上的作用面積很大。這些說明，最大氣動力效應(yīng)工況應(yīng)該出現(xiàn)在雙機(jī)工作的45°或90°位置。表 3為六種工況下作用點為太陽能帆板的中心點的力和力矩的值。

圖13 算例6太陽能帆板受到的羽流氣動壓強(qiáng)分布圖

從表3可以看出，發(fā)動機(jī)羽流對太陽能帆板影響最大的工況出現(xiàn)在算例5，即雙機(jī)工作時帆板處于45°角位置時，其作用在帆板中心點上的合力為約為24N，合力矩約為7.7Nm。

5 結(jié)束語

PWS軟件是一款基于DSMC的羽流計算通用軟件，并且具備了可擴(kuò)展特性，松散的架構(gòu)為軟件的進(jìn)一步完善提供了方便的接口，軟件的計算值與國外“CUBRC”試驗測量值符合得很好，認(rèn)為PWS軟件具有較高的置信度。

使用PWS軟件對神舟飛船平移發(fā)動機(jī)羽流效應(yīng)進(jìn)行計算表明：發(fā)動機(jī)羽流對著太陽能帆板羽流氣動壓強(qiáng)最大值出現(xiàn)在的算例1和算例4，即單機(jī)和雙機(jī)工作時帆板處于0°角位置時，但作用力不大；發(fā)動機(jī)羽流對太陽能帆板影響最大的工況出現(xiàn)在算例5，即雙機(jī)工作時帆板處于45°角位置時，其作用在帆板中心點上的合力約為24N，合力矩約為7.7Nm。

平移發(fā)動機(jī)對處于羽流直接撞擊區(qū)的太陽能帆板有較大的力效應(yīng)作用，需要結(jié)合飛船質(zhì)心進(jìn)行進(jìn)一步分析，數(shù)值模擬結(jié)果可以為神舟飛船平移發(fā)動機(jī)使用模式提供一定參考作用。 ◇

［1］Bird G A.Molecular gas dynamics［M］.Oxford:Clarendon Press，1976

［2］Wilmoth R G，Carlson A B，LeBeau G J.DSMC grid methodologies for computing low-density.hypersonic flows about reusable launch vehicles，AIAA 1996-1812［R］

［3］Boyles K A，LeBeau G J，Lumpkin III F E.The use of virtual subcells in DSMC analysis of orbiter aerodynamics at high altitudes upon reentry.AIAA 2003-1030［R］

［4］LeBeau G J，Boyles K A，Lumpkin III F E.Virtual sub-cells for the direct simulation Monte Carlo method.AIAA 2003-1031［R］

［5］Ivanov M S，Markelovf G N，Gimelshein S F.Statistical simulation of reactive rarefied flows:numerical approach and applications.AIAA 1998-2669［R］

［6］Dietrich S，Boyd I D.A scalar optimized parallel implementation of the DSMC method.AIAA 1994-355［R］

［7］Dietrich S，Boyd I D.Parallel implement at ion on the IBM SP-2 of the Direct Simulation Monte Carlo method.AIAA 1995-2029[R]

［8］沈青.稀薄氣體動力學(xué)［M］.北京：國防工業(yè)出版社，2003

［9］Holden M，Wadhams T.Code validation study of laminar shock/boundary layer and shock/shock interactions in hypersonic flow.Part A:experimental measurements，AIAA 2001-1031［R］