羽流中固體顆粒在真空環(huán)境下分布的數(shù)值仿真①

李中華，李志輝，吳俊林，彭傲平

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

羽流中固體顆粒在真空環(huán)境下分布的數(shù)值仿真①

李中華，李志輝，吳俊林，彭傲平

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

通過對氣-固兩相間動量和能量相互作用解耦處理，建立了一種適于模擬真空環(huán)境氣固兩相混合物羽流的DSMC雙向耦合算法與固體顆粒空間輸運(yùn)變化特性的TPMC計算技術(shù)。仿真了固體火箭發(fā)動機(jī)兩相羽流流場和固體顆粒在遠(yuǎn)離發(fā)動機(jī)噴口數(shù)十公里空間擴(kuò)散運(yùn)動分布特性，通過將計算結(jié)果與典型文獻(xiàn)結(jié)果及理論分析比較確認(rèn)，證實(shí)本文方法的準(zhǔn)確可靠性。結(jié)果表明，固體顆粒對氣相的擴(kuò)散有一定的阻滯作用;僅在離發(fā)動機(jī)噴口一定距離，氣相對固體顆粒有較大影響，致顆粒溫度下降、速度增加;顆粒溫度隨發(fā)動機(jī)噴口距離增大而減小，一直要在遠(yuǎn)離噴口上百公里顆粒溫度才隨軸向位置趨于平衡，且不同尺寸的顆粒溫度差別較大，對指導(dǎo)外層空間高真空環(huán)境氣固兩相羽流傳輸影響工程研制具有重要意義。

稀薄氣體;DSMC方法;兩相流;雙向耦合;TPMC方法;高空羽流

0 引言

固體火箭發(fā)動機(jī)廣泛采用含金屬復(fù)合推進(jìn)劑，燃燒后會生成大量的固體顆粒。燃?xì)馍淞魍ǔＪ怯啥喾N氣體成分、固體顆粒形成典型的兩相多組分混合物流動，產(chǎn)生復(fù)雜的相互作用干擾流場［1－6］。如何準(zhǔn)確可靠地模擬真空環(huán)境氣固兩相混合物羽流傳輸及固體顆粒遠(yuǎn)離噴口擴(kuò)散運(yùn)動軌跡、分布變化規(guī)律，一直是國內(nèi)外學(xué)術(shù)界與工程應(yīng)用部門關(guān)心的問題。

在高空，一般采用DSMC(Direct Simulation Monte Carlo)方法［7］模擬羽流流場。該方法采用有限數(shù)目的仿真分子模擬實(shí)際流場中數(shù)目巨大的真實(shí)分子。通過跟蹤流場中仿真分子的運(yùn)動和分子間的碰撞達(dá)到流場模擬的目的。該方法廣泛用于模擬高空火箭或太空船推進(jìn)器的羽流，在描述這些流動中的氣體特性時有很高的精確度［8－9］。為了模擬氣固兩相流問題，Gallis提出了一種改進(jìn)的DSMC模擬方案［10］，利用Green函數(shù)發(fā)展的DSMC方法適用于求解在任意分子速度分布的氣相流場中顆粒所受的力和熱，可模擬包括稀薄和化學(xué)惰性固體顆粒相在內(nèi)的稀薄流動。但Gallis方法只考慮氣相對固體顆粒的作用，忽略顆粒相作用于氣體的影響。這種假定在固體推進(jìn)劑火箭羽流中是有一定缺陷的，因?yàn)樵诠腆w火箭發(fā)動機(jī)羽流中，很多模擬范圍內(nèi)相間動量和能量的傳遞對氣體流動特性有重大的影響。后來經(jīng)Burt、Boyd等人的發(fā)展，建立了一種適用于DSMC方法的雙向耦合(Two-way coupled)算法［5，11－13］，既考慮氣相對固相的力和熱的作用，又考慮固相顆粒對氣相的作用，能較準(zhǔn)確地描述固相顆粒在稀薄過渡流中的輸運(yùn)過程。

對于在大氣層外真空環(huán)境工作的發(fā)動機(jī)，其羽流流場中含有大量Al2O3粒子，且大小不一。這些粒子會在空間形成一團(tuán)粒子流云團(tuán)，對相關(guān)的紅外觀測設(shè)備造成污染影響，而且粒子云團(tuán)的存在還會嚴(yán)重干擾紅外觀測效果。對于長時間高真空環(huán)境工作的發(fā)動機(jī)，會出現(xiàn)羽流中固體顆粒如Al2O3粒子形成約數(shù)公里長的粒子流，隨后又運(yùn)動較長的時間。Al2O3粒子在數(shù)公里長的粒子流中的數(shù)密度、溫度分布難以按照常規(guī)宏觀流體力學(xué)數(shù)值方法和關(guān)于空間位置、速度的飛行力學(xué)數(shù)值推進(jìn)計算確定。

為了研究大氣層外高真空環(huán)境羽流場中固體顆粒空間分布、變化特性，本文在前期研究［9，14－16］基礎(chǔ)上，提出發(fā)展適于氣固兩相羽流DSMC雙向耦合算法與TPMC(Test Particle Monte Carlo)模擬技術(shù)，采用DSMC和TPMC相結(jié)合的方法開展研究，擬分兩步進(jìn)行計算。第一步，首先考慮發(fā)動機(jī)固定不動，把羽流視為一個定常問題，利用DSMC方法計算得到穩(wěn)定噴流中流場溫度、粒子數(shù)密度及粒子速度。在DSMC方法中，采用雙向耦合算法計算固體顆粒與噴流氣體分子之間的相互作用。第二步，在第一步計算的基礎(chǔ)上，分析得到發(fā)動機(jī)噴流徑向剖面上的流場溫度、粒子數(shù)密度及粒子速度基本不再受噴流氣體影響的剖面計算參數(shù)(包括粒子數(shù)密度、溫度、流動速度)，以此數(shù)據(jù)為基礎(chǔ)，采用TPMC方法計算粒子在遠(yuǎn)離發(fā)動機(jī)噴口的空間位置分布，得到不同時刻粒子的溫度和數(shù)密度分布。

1 氣-固兩相羽流DSMC雙向耦合算法

DSMC方法是羽流流場計算行之有效的模擬方法。依托DSMC方法，采用變剛球分子模型，能量交換使用Larsen-Borgnakke模型［17］;網(wǎng)格采用二級直角網(wǎng)格，碰撞網(wǎng)格根據(jù)密度自適應(yīng)［9，14－15］，碰撞對的選取限制在碰撞網(wǎng)格內(nèi)，發(fā)展考慮氣-固兩相相互作用影響的DSMC雙向耦合算法。為了提高計算效率，采用了基于MPI的DSMC并行計算技術(shù)，將計算區(qū)域分解為N個子區(qū)域，分別分配給N個節(jié)點(diǎn)。各節(jié)點(diǎn)進(jìn)行獨(dú)立且完整的DSMC仿真。仿真分子穿過邊界區(qū)域時，計算區(qū)域之間進(jìn)行仿真數(shù)據(jù)信息交換。兩相流DSMC模擬運(yùn)算準(zhǔn)則是基于相間動量和能量從微粒特性的瞬時變化的傳遞進(jìn)行解耦，分為兩個方面:

(1)考慮氣相對固體顆粒的作用［10］。假設(shè)固體顆粒處于當(dāng)?shù)刈杂煞肿恿鞯臓顟B(tài)，不考慮多原子氣體的振動激發(fā)，在同一個網(wǎng)格里氣相的每個DSMC仿真分子作用到一個固體顆粒上的能量和動量可表示為

式中 Rp為等效顆粒半徑;Ng為每個仿真分子所模擬的真實(shí)氣體分子數(shù)目;τ為顆粒表面熱適應(yīng)系數(shù);Vc為網(wǎng)格體積;m為單個氣體分子質(zhì)量;ur為氣體分子與相關(guān)顆粒的相對速度，cr是ur的絕對值;k為Boltzmann常數(shù);Tp為顆粒溫度;Λ為氣體轉(zhuǎn)動自由度數(shù);erot為單個分子轉(zhuǎn)動能。

(2)考慮固體顆粒對周圍氣體的影響，首先要確定在每個時間步長內(nèi)哪個仿真分子將與顆粒進(jìn)行碰撞。通過對Bird的非時間計數(shù)方法［7］進(jìn)行修正，來確定與所選的顆?？赡馨l(fā)生碰撞的仿真分子數(shù)ns。

式中 Np為一個仿真顆粒所表示的實(shí)際固體顆粒數(shù)量;ng為與固體顆粒在同一網(wǎng)格里的氣體仿真分子數(shù);Δt為時間步長;(cr)max為網(wǎng)格內(nèi)采樣到的分子-顆粒對碰撞前最大相對速度。

一個與這個顆粒發(fā)生碰撞的給定的氣體仿真分子，要么為以概率等于顆粒熱適應(yīng)系數(shù)τ的等溫漫反射碰撞，要么為概率為1-τ的鏡面反射。如果發(fā)生鏡面反射，則相對速度cr在碰撞中不發(fā)生改變，碰撞后的相對速度可通過cr與單位矢量?n相乘得到。如果發(fā)生漫反射，碰撞后相對速度圍繞初始相對速度的方位角ε在［0，2π］上等概率分布。在漫反射碰撞中，碰撞后相對速度能假定為等于初始相對速度cr，而是需要通過使用“取舍”法，從如下分布函數(shù)來確定的值:

式中 β為氣體基于顆粒溫度的最可幾熱運(yùn)動速度的倒數(shù)，β=［m/(2kTp)］1/2。

2 真空環(huán)境固體顆粒分布TPMC計算

大氣層外真空環(huán)境發(fā)動機(jī)長時間工作，固體顆粒會隨羽流擴(kuò)散運(yùn)動到離發(fā)動機(jī)噴口很遠(yuǎn)的距離，這時再采用上文介紹的DSMC和雙向耦合算法來模擬顆粒運(yùn)動分布，顯然不再合適。為了跟蹤固體顆粒長時間的運(yùn)動軌跡，本文建立了一種粒子跟蹤的TPMC方法。發(fā)動機(jī)羽流在真空環(huán)境膨脹到一定程度后，流場中氣流密度會變得很低，氣相和固相不再相互作用。這時，固體顆粒保持勻速直線運(yùn)動。在這樣的羽流場，模擬固體顆粒的TPMC方法基本步驟可設(shè)計為

(1)在羽流場入口處產(chǎn)生一個樣本粒子(試驗(yàn)粒子)。根據(jù)DSMC方法得到羽流場的結(jié)果，在TPMC界面上隨機(jī)產(chǎn)生一個固體顆粒。

(2)跟蹤顆粒的運(yùn)動。顆粒在起始點(diǎn)以本身所具有的速度做勻速直線運(yùn)動。在給定時間步長內(nèi)，運(yùn)動到新的位置。運(yùn)動一定時間后，記錄顆粒的位置。

(3)在大氣層外真空環(huán)境中運(yùn)動，固體顆粒要吸收太陽輻射能量，本身也要向外輻射能量，溫度會發(fā)生變化，要進(jìn)行溫度變化的處理。

(4)重復(fù)(1)～(3)的過程，直到記錄到足夠多顆粒為止。

還需要建立在高真空環(huán)境考慮太陽輻射影響的可計算模型。假設(shè)太陽輻射是均勻的，單位時間固體顆粒受到太陽輻射的能量為

式中 Es為固體顆粒受到的太陽輻射能量功率;fa為太陽常數(shù);在外層空間高真空環(huán)境可設(shè)置fa=1 368 W/m2;Sp為顆粒等效受熱面積;α為發(fā)射系數(shù)和吸收系數(shù)。

在Δt時間步長內(nèi)，顆粒受太陽輻射引起的溫度變化可表示為

式中 cp為顆粒材料的比重壓熱容;Mp為顆粒質(zhì)量。

顆粒由自身輻射引起的溫度變化應(yīng)為

式中 σ為Stefan-Boltzmann常數(shù);cp為顆粒材料的熱容;ρp為顆粒材料密度。

在每個時間步長內(nèi)，對式(8)和式(9)進(jìn)行迭代求解，可給出顆粒的溫度變化。

顆粒數(shù)分布按顆粒半徑的連續(xù)分布密度，可近似地由對數(shù)正態(tài)分布曲線來擬合［18］。即

3 算例與結(jié)果討論

3.1 氣-固相互作用模型驗(yàn)證

為了比較確認(rèn)本文發(fā)展的兩相流雙向耦合算法的準(zhǔn)確可靠性，基于文獻(xiàn)［11］提供兩相稀薄流中一個氣-固相互作用影響模擬的雙向耦合模型算例，擬定0.1 mm寬和20 mm長的矩形區(qū)域(圖1)進(jìn)行計算。

模擬中，氣相混合物為H2、CO和N2組成的混合氣體，分子數(shù)密度分別設(shè)置為 2×1023、1×1023、1×1023m－3;在入口處，氣體速度為2 000 m/s，溫度為1 000 K。固相有2 種鋁粒子:直徑分別為 3×10－6、6×10－6m，每種粒子的質(zhì)量流量相等;在入口邊界上，粒子的速度為1 200 m/s，溫度為 2 200 K。

圖1 氣固兩相流計算區(qū)域與邊界條件示意圖Fig.1 Sketch of domain and boundary types

圖2和圖3繪出了流場中不同截面上動量與能量傳遞速率隨位置的變化關(guān)系。圖中顯示，氣-固兩相之間有較強(qiáng)的動量和能量傳遞。

圖2 兩相流中動量傳遞速率隨位置變化比較Fig.2 Comparison of momentum transfer rates

圖3 兩相流中能量傳遞速率隨位置變化比較Fig.3 Comparison of energy transfer rates

通過比較氣相和固相的變化梯度，可看出兩者的梯度大小大致相等，符號相反。這表明在模擬計算中，氣相作用在固相的動量和能量，等于固相反作用于氣相的動量和能量，即在相間相互作用時，動量和能量是守恒的。在前面介紹的雙向耦合算法中，氣相對固相的作用和固相對氣相的作用過程是解耦的，兩個過程獨(dú)立計算。在計算固相對氣相的作用時，其過程是一個隨機(jī)統(tǒng)計計算過程，不能保證在每一步的計算中嚴(yán)格地遵守動量和能量守恒，但在大樣本數(shù)情況下統(tǒng)計平均時，可保證動量和能量守恒。通過圖2、圖3將本文仿真模擬方法結(jié)果與典型文獻(xiàn)結(jié)果［11］進(jìn)行定量化對比分析，表明本文所建立的氣-固兩相流動DSMC雙向耦合算法用于計算仿真固體火箭發(fā)動機(jī)羽流中氣相混合物與固體顆粒輸運(yùn)過程是可靠的。

3.2 固體發(fā)動機(jī)兩相羽流流場計算

本文對某固體發(fā)動機(jī)真空環(huán)境產(chǎn)生的兩相羽流流場進(jìn)行模擬，其中發(fā)動機(jī)出口馬赫數(shù)1.5，出口溫度2 960 K，出口速度設(shè)為均勻流動，固體顆粒的質(zhì)量流量占總流量的30%。羽流中氣相組元的摩爾分?jǐn)?shù)見表1。

表1 羽流氣相組元摩爾分?jǐn)?shù)Table 1 Mole fraction of plume gas species

計算中，背景環(huán)境設(shè)為真空邊界條件。關(guān)于固體顆粒Al2O3的模擬，采用有限分組方法，根據(jù)式(10)，把Al2O3分成不同尺寸和質(zhì)量的5組，視為5種顆粒來處理。每一組顆粒具有相同的尺寸和質(zhì)量，顆粒直徑在20～200 μm之間，具體尺寸見表2。

表2 固體顆粒尺寸Table 2 Particles size

圖4給出了計算得到的兩相羽流流場數(shù)密度分布等值線云圖。從圖4(a)的氣相流場可看出，由于噴流出口馬赫數(shù)較低，只比音速稍高，噴流向真空膨脹十分迅速。沿流動方向，數(shù)密度迅速下降，溫度下降也很快。沿徑向，數(shù)密度下降更快。固相顆粒的擴(kuò)散是氣相擴(kuò)散引起的，圖4(b)中，顯示出固相顆粒的擴(kuò)散是有界的。當(dāng)氣相膨脹到密度很低時，對固相顆粒的影響會變得很小，顆粒會保持自己的運(yùn)動狀態(tài)做勻速直線運(yùn)動。從數(shù)密度分布來看，固體顆粒比較集中地分布在核心區(qū)域。為了驗(yàn)證固體顆粒對氣相羽流的影響，還進(jìn)行了不包含固體顆粒的單相羽流對比計算。圖5給出了考慮固體顆粒對氣相影響的兩相羽流與單相羽流2種情況下模擬得到的羽流軸線上氣相溫度與流動速度分布比較?？煽闯?，單相羽流模擬得到的氣體流動速度較雙相羽流模擬值普遍偏大，說明在沒有固體顆粒的情況下，氣相流場膨脹得更快，速度增加更多，在x－=6處(以噴流半徑為參考無量綱化)2種情況模擬得到的速度相差約160 m/s，證實(shí)固體顆粒對氣相的擴(kuò)散有一定的阻滯作用。兩相情況下，氣體不斷從顆粒相吸收能量，溫度較高，圖中顯示出在x－=6處，兩相羽流模擬得到的氣相溫度高于單相羽流模擬值約40 K。

圖6給出了該發(fā)動機(jī)兩相羽流場中固體顆粒溫度分布等值線云圖，呈現(xiàn)出復(fù)雜流場結(jié)構(gòu)，沿軸向，隨著氣相的迅速膨脹，溫度迅速下降，能量從固體顆粒向氣相轉(zhuǎn)化，顆粒的溫度也逐漸下降。沿徑向，軸線附近顆粒的溫度比外圍的溫度稍低，更外圍的顆粒溫度更低，這是因?yàn)楸?設(shè)置的5種顆粒與氣相相互作用時能量交換不同，導(dǎo)致顆粒的溫度不同;同時，動量交換也不相同，顆粒擴(kuò)散的程度因此會有所差異，這些因素導(dǎo)致溫度的分布比較復(fù)雜。不過整個流場中顆粒溫度分布差別不大，最大差異在300 K左右。由圖7所示的顆粒軸向速度分布表明，外圍顆粒被氣相加速較大，這是因?yàn)橥鈬饕切〕叽绲念w粒，顆粒尺寸越小，越容易被氣相加速。

圖4 兩相羽流流場數(shù)密度分布Fig.4 Number density distribution in two-phase plume flow field

圖5 單相和兩相羽流軸線上溫度與流動速度分布比較Fig.5 Comparison on temperature and flow velocity distribution of one-and two-phase plume flows

圖6 羽流中固體顆粒溫度分布Fig.6 Particles temperature distribution in two-phase plume flow field

圖7給出5種顆粒沿羽流軸線不同位置的溫度、速度分布，顯示出在距離噴口x－=5以內(nèi)的羽流場中，氣相對固體顆粒有較大影響，顆粒溫度下降和速度增加較快，由此計算證實(shí)了式(1)、式(2)表征的氣相分子數(shù)對固相顆粒動量和能量的作用影響很大;但在x－=5處，氣相數(shù)密度比出口處下降1個量級，此時氣相對固體顆粒的影響已經(jīng)變得很小。x－＞5之后，不同位置的顆粒溫度下降和速度增加較小，基本上各種顆粒的溫度、速度分布均保持在相應(yīng)量值。

圖7 不同固體顆粒沿軸線溫度、速度分布計算比較Fig.7 Comparison on temperature and velocity distribution of particles along the axis of the plume flow

為了揭示該發(fā)動機(jī)長時間真空環(huán)境工作所產(chǎn)生兩相羽流場固體顆粒運(yùn)動軌跡與分布特性，利用所建立TPMC方法，對發(fā)動機(jī)工作70 s后的顆粒分布進(jìn)行了模擬計算。由于考慮太陽輻射和自身輻射時，發(fā)射系數(shù)對顆粒溫度變化過程有一定影響。根據(jù)文獻(xiàn)［19］中的結(jié)果，Al2O3在 800～1 300 K 之間、波長為 1.89 μm的發(fā)射系數(shù)比較接近0.5，而且隨溫度變化不大。因此，在計算中假設(shè)，Al2O3在所有溫度范圍、各波段內(nèi)的發(fā)射系數(shù)均為0.5。圖8給出了發(fā)動機(jī)工作70 s后，第1種顆粒在中心對稱面的數(shù)密度、溫度等值線云圖分布。計算表明，此時在羽流軸線方向，粒子擴(kuò)散到離發(fā)動機(jī)噴口63～115 km范圍的扇形區(qū)域內(nèi)。在周向，基本上保持軸對稱分布，由此顆粒分布是在一個近似圓臺的范圍內(nèi)，擴(kuò)散范圍很大。圖9展示了本文計算設(shè)置的5種顆粒分別在離發(fā)動機(jī)噴口74、115 km的擴(kuò)散半徑(R)與顆粒直徑(dp)變化關(guān)系。由圖9分析，顆粒直徑越大，擴(kuò)散半徑越小，兩者基本上呈反比關(guān)系。

圖10繪出了該發(fā)動機(jī)羽流軸對稱中心線上不同顆粒的溫度隨遠(yuǎn)離噴口軸向位置60 km＜x＜110 km變化關(guān)系。圖10中顯示出，顆粒溫度隨遠(yuǎn)離噴口位置增大而減小;顆粒尺寸越大，起始溫度越高，在整個擴(kuò)散過程中，也保持較高的溫度。在發(fā)動機(jī)工作達(dá)70 s時，尺寸最小的第1種顆粒溫度隨位置x基本保持不變而趨近平衡溫度，其他顆粒仍隨軸向位置x增大而減小，還未到達(dá)平衡溫度。這種現(xiàn)象說明，受太陽輻射和自身輻射的影響，顆粒溫度隨著遠(yuǎn)離發(fā)動機(jī)噴口的擴(kuò)散運(yùn)動而逐漸下降，一直要在遠(yuǎn)離發(fā)動機(jī)噴口上百公里顆粒溫度，才隨軸向位置x逐漸趨于平衡，且不同尺寸的顆粒溫度也各不相同。這對指導(dǎo)固體發(fā)動機(jī)高真空環(huán)境氣固兩相羽流空間環(huán)境影響工程研制具有重要意義。

圖8 發(fā)動機(jī)工作70 s后顆粒1參數(shù)分布Fig.8 Number density and temperature distribution for Particle 1 after 70 s with engine working

圖9 顆粒擴(kuò)散半徑與顆粒直徑的關(guān)系Fig.9 Relation between diffused radius and particle diameters

圖10 不同顆粒溫度分布比較Fig.10 Comparison on temperature distribution for different types of particles

4 結(jié)論

(1)采用雙向耦合技術(shù)，建立了適于氣固兩相多組元混合物羽流問題DSMC仿真方法，通過對相間相互作用進(jìn)行解耦處理，實(shí)現(xiàn)了氣-固兩相間動量和能量相互作用的模擬。研究證實(shí)固體顆粒對氣相的擴(kuò)散有一定的阻滯作用，會使氣相溫度增加。顆粒溫度下降，且沿徑向軸線附近顆粒溫度比外圍溫度稍低。

(2)在距離發(fā)動機(jī)噴口一定位置，如x－=5以內(nèi)的羽流場中，氣相對固體顆粒有較大影響，顆粒溫度下降和速度增加，而在x－＞5之后不同位置，氣相對固體顆粒的影響甚微，顆粒溫度和速度逐漸趨于相應(yīng)平衡值。

(3)研制了一套適于真空環(huán)境固體顆粒空間分布運(yùn)動特性模擬的TPMC方法。通過對發(fā)動機(jī)長時間工作(如70 s)所產(chǎn)生顆粒分布進(jìn)行模擬計算，結(jié)果表明，顆粒擴(kuò)散集聚到離發(fā)動機(jī)噴口63～115 km范圍的近似圓臺區(qū)域;在一定氣體數(shù)密度范圍內(nèi)，羽流中氣-固兩相有較強(qiáng)的相互作用，顆粒尺寸對顆粒的擴(kuò)散范圍和溫度分布有很大影響。顆粒直徑越大，擴(kuò)散半徑越小，兩者基本上呈反比關(guān)系。

(4)受太陽輻射和自身輻射的影響，顆粒溫度隨著遠(yuǎn)離發(fā)動機(jī)噴口的擴(kuò)散運(yùn)動而逐漸下降，直至遠(yuǎn)離發(fā)動機(jī)噴口上百公里顆粒溫度，才隨軸向位置逐漸趨于平衡，且不同尺寸的顆粒溫度也各不相同。對指導(dǎo)固體發(fā)動機(jī)高真空環(huán)境氣固兩相羽流空間環(huán)境影響工程研制具有重要作用。

(5)通過DSMC和TPMC方法相結(jié)合，本文方法為稀薄兩相流的計算、大尺度范圍內(nèi)計算高空羽流中顆粒輸運(yùn)過程提供了一種數(shù)值模擬手段。研究體會到，固體顆粒尺寸對顆粒的擴(kuò)散范圍和溫度分布有很大影響。在一定氣體數(shù)密度范圍內(nèi)，羽流中氣-固兩相有較強(qiáng)的相互作用。在氣體數(shù)密度很小時，氣-固兩相作用很弱，可不予考慮。

［1］Dettleff G.Plume flow and plume impingement in space technology［J］.Progress in Aerospace Sciences，1991，28(1).

［2］Siebenhaar A，Bulman M J.The strutjet engine:the overlooked option for space launch［R］.AIAA 95-3124.

［3］于勝春，湯龍生.固體火箭發(fā)動機(jī)噴管及羽流流場的數(shù)值分析［J］.固體火箭技術(shù)，2004，27(2).

［4］李東霞，徐旭，蔡國飆，等.火箭發(fā)動機(jī)氣體-顆粒兩相流雙流體模型研究［J］.固體火箭技術(shù)，2005，28(4).

［5］Marris H，Deboudt K，Augustin P，et al.Fast changes in chemical composition and size distribution of fine particles during the near-field transport of industrial plumes［J］.Science of The Total Environment，2012，427-428:126-138.

［6］Johnson G R，Jayaratne E R，Lau J，et al.Remote measurement of diesel locomotive emission factors and particle size distributions［J］.Atmospheric Environment，2013，81:148-157.

［7］Bird G A.Molecular gas dynamics and the direct simulation of gas flows［M］.London:Oxford Univ.Press，1994.

［8］Ivanov M S，Khotyanovskyy D V，Kudryavtsev A N，et al.Numerical study of backflow for nozzle plumes expanding into Vacuum［R］.AIAA 2004-2687.

［9］李志輝，李中華，楊東升，等.衛(wèi)星姿控發(fā)動機(jī)混合物羽流場分區(qū)耦合計算研究［J］.空氣動力學(xué)學(xué)報，2012，30(4):483-491.

［10］Gallis M A，Rader D J，Torczynski J R.DSMC simulations of the thermophoretic force on a spherical macroscopic particle［R］.AIAA 2001-2890，American Institute of Aeronautics and Astronautics，Washington，D C，2001.

［11］Burt J M，Boyd I D.Development of a two-way coupled model for two phase rarefied flows［C］//42nd AIAA Aerospace Sciences Meeting and Exhibit，Reno，Nevada，5-8 Jan.2004，AIAA 2004-1351.

［12］Burt J M，Boyd I D.Monte carlo simulation of a rarefied multiphase plume flow［C］//43rd AIAA Aerospace Sciences Meeting and Exhibit 10-13 January 2005，Reno，Nevada.2005，AIAA 2005-964.

［13］Harris A J L，Ripepe M，Hughes E A.Detailed analysis of particle launch velocities，size distributions and gas densities during normal explosions at Stromboli［J］.Journal of Volcanology and Geothermal Research，2012，231-232:109-131.

［14］李中華.飛行器材料高空放氣污染蒙特卡羅數(shù)值模擬［C］//中國宇航學(xué)會首屆學(xué)術(shù)年會，2005.

［15］李中華，李志輝，李海燕，等.過渡流區(qū)NS/DSMC耦合算法研究［J］.空氣動力學(xué)學(xué)報，2013，31(3).

［16］Li Z H，Wu J L，Peng A P.Coupled navier-stokes/direct simulation monte carlo simulation of multicomponent mixture plume flows［J］.Journal of Propulsion and Power，2014，30(3):672-689.

［17］Borgnakke C，Larsen P S.Statistical collision model for Monte Carlo simulation of polyatomic gas mixtures［J］.J.of Comput.Phys.，1975，18:405-420.

［18］李潔，任兵，陳偉芳.稀薄流過渡區(qū)氣固兩相噴流的建模與數(shù)值模擬［J］.空氣動力學(xué)學(xué)報，2005，23(4).

［19］Freeman G N，Ludwing C B，Malkmus W，et al.Development and validation of standardized infrared radiation model(SIRRM)［R］.AFRPL 79-55.

(編輯:崔賢彬)

Numerical simulation on particle distribution of the plume flow in vacuum

LI Zhong-hua，LI Zhi-hui，WU Jun-lin，PENG Ao-ping
(China Aerodynamics Research and Development Center，Mianyang 621000，China)

A numerical simulation approach was presented to simulate rarefied two-phase mixture plume flow by applying twoway coupling technique of DSMC method.The interaction between rarefied gas and solid particles was dealt with decoupling computation of momentum and energy exchange between phases.The TPMC computing technology was also developed to solve the space transport characteristics of solid particles in distant field space.The two-phase plume flow field and particle distribution characteristics far away from the engine nozzle were simulated by the DSMC and TPMC method.It is indicated that the present methods and models are correct by the comparison of the computed results with the typical confirmed results from the reference and theoretical analysis.There exists some retarded elasticity in gas phase diffusion of solid particles to some extent.From the engine nozzle distance，gas and solid particles have a greater impact caused with temperature drop and increases of the speed of particles.The granular temperature away from the nozzle exit tends to balance，and the particle temperature is very different with different size of particles.It is shown that there are more intensely interaction between phases in case of higher number density，and particle size affects diffusing range and temperature distribution.The results are of great significance to guide the vacuum gas-solid two-phase plume transmission.

rarefied gas;DSMC method;two-phase flow;two-way coupled method;TPMC method;plume flow

V438

A

1006-2793(2014)06-0797-07

10.7673/j.issn.1006-2793.2014.06.011

2014-03-03;

2014-08-12。

973計劃(2014CB744100);國家自然科學(xué)基金(91016027、11325212);國防基礎(chǔ)科研基金(51313030104)。

李中華(1970—)，男，高工，研究方向?yàn)橄”怏w動力學(xué)。E-mail:lzhzcb@aliyun.com

李志輝。E-mail:zhli0097@x263.net