蔡國(guó)飆 賀碧蛟

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

1 引言

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

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

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

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

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

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

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

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

圖1 DSMC計(jì)算流程圖

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

3 PWS軟件及驗(yàn)證

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

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

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

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

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

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

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

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

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

圖5 雙圓錐體熱流密度分布計(jì)算值與試驗(yàn)值對(duì)比

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

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

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

表1 對(duì)比算例

圖6 發(fā)動(dòng)機(jī)與太陽(yáng)能帆板相對(duì)位置

4.1 計(jì)算條件

（1）粒子入口條件

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

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

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

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

（2）邊界條件

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)束語(yǔ)

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

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

平移發(fā)動(dòng)機(jī)對(duì)處于羽流直接撞擊區(qū)的太陽(yáng)能帆板有較大的力效應(yīng)作用，需要結(jié)合飛船質(zhì)心進(jìn)行進(jìn)一步分析，數(shù)值模擬結(jié)果可以為神舟飛船平移發(fā)動(dòng)機(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］沈青.稀薄氣體動(dòng)力學(xué)［M］.北京：國(guó)防工業(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］