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

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

一種稀薄兩相流動的數(shù)值模擬方法

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

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

從描述粒子運動的微觀層次出發(fā)，采用雙向耦合技術，建立了一種適用于稀薄條件下兩相流動的DSMC數(shù)值模擬方法。對相間相互作用進行解耦處理，實現(xiàn)了氣固兩相間動量和能量相互作用的模擬。相關算例驗證了方法的可靠性。結果表明：在大樣本統(tǒng)計平均條件下，該模擬方法能夠保證相間動量和能量交換時守恒，固體顆粒尺寸、氣體稀薄程度等條件對相間的相互作用有很大的影響。

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

0 引 言

固體火箭發(fā)動機為了達到提高比沖的目的，廣泛采用含金屬復合推進劑。燃燒后生成大量的固體顆粒，燃氣射流通常是由多種氣體成分、固體顆粒形成典型的兩相流動，這些顆粒既對燃氣排導裝置產生巨大的沖刷侵蝕作用，又對燃氣射流流場產生直接影響，形成復雜的相互作用干擾流場［1-3］。

對于高空羽流，一般采用有效仿真稀薄氣體流動的DSMC（Direct Simulation Monte Carlo）方法模擬。該方法采用有限數(shù)目的仿真分子模擬實際流場中數(shù)目巨大的真實分子。通過跟蹤流場中仿真分子的運動和分子間的碰撞達到流場模擬的目的。該方法廣泛地用于模擬高空火箭或太空船推進器的羽流，在描述這些流動中的氣體特性時有很高的精確度［4-7］。Gallis提出了一種改進的DSMC方法［8］，利用Green函數(shù)發(fā)展的DSMC方法適用于求解在任意分子速度分布的氣相流場中顆粒所受的力和熱，可以模擬包括稀薄和化學惰性固體顆粒相在內的稀薄流動。Gallis方法只考慮氣相對固體顆粒的作用，忽略顆粒相作用于氣體的影響。這種假定在固體推進劑火箭羽流中是有一定缺陷的，因為在固體火箭發(fā)動機羽流中很多模擬范圍內相間動量和能量的傳遞對氣體流動特性有重大的影響。后來經Burt、Boyd等人的發(fā)展，建立了一種適用于DSMC方法的雙向耦合算法，既考慮氣相對固相的力和熱的作用，又考慮固相顆粒對氣相的作用，能夠準確描述固相顆粒在稀薄過渡流中的輸運過程［9-15］。

本文采用雙向耦合技術，建立了模擬稀薄兩相流的DSMC方法。對相關文獻上的算例計算，模擬了氣固兩相羽流，驗證了計算方法，為稀薄條件下的氣固兩相流動問題提供了一種新的研究手段。

1 數(shù)值模擬算法

兩相流DSMC模擬運算法則是基于相間動量和能量從微粒特性的暫時變化的傳遞的解耦。

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

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

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

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

式中，β為氣體在顆粒溫度處最可幾速度的倒數(shù)，β＝［m/（2kTp）］1/2。

對于漫反射多原子分子氣體，碰撞后轉動能erot也必須改變。漫反射雙原子氣體分子的轉動能可計算如下：

式中，Rf為（0，1）之間的一個隨機數(shù)。

整體坐標系下ur的分量為ur、vr、wr。采用Bird二元彈性碰撞，相對速度的各分量可以由ur、vr、wr、cr和、角δ和ε計算得到。有：

式中，δ為碰撞偏轉角，定義為-ur與碰撞后的相對速度矢量之間的夾角，這里分別為下面碰撞時氣體分子和顆粒的絕對速度。

應用SPSS 22.0軟件，計數(shù)資料比較采用卡方檢驗，采用Spearman相關分析法分析癌組織不同NF-κB與HMGB1表達與其T分期和N分期的關系，P<0.05為差異有統(tǒng)計學意義。

由于增加了大量的固體顆粒，并需要計算氣-固兩相的相互作用，計算量會大大增加。為了提高計算效率，本文采用了基于MPI的并行計算技術。將計算區(qū)域分解為N個子區(qū)域，分別分配給N個節(jié)點。各節(jié)點進行獨立且完整的兩相DSMC仿真。仿真粒子（包含氣相、固相）穿過區(qū)域邊界時，計算區(qū)域之間進行數(shù)據(jù)交換。

2 算例與討論

2.1 算例驗證

為了驗證本文建立兩相流DSMC雙向耦合算法可靠性，對文獻［10］的算例進行了計算。計算區(qū)域確定為0.1mm寬和20mm長的矩形區(qū)域（圖1），上下為鏡面反射邊界，在該邊界上不會發(fā)生能量和動量交換，這樣保證流動過程中只在氣體分子與固體顆粒之間發(fā)生能量和動量交換。入口邊界為均勻流動，出口邊界為超聲速邊界。模擬中氣相混合物為H2、CO和N2組成的混合氣體，分子數(shù)密度分別為2×1023/m3、1×1023/m3、1×1023/m3；在入口處，氣體的速度為2 000m/s，溫度為1 000K。固相有兩種鋁粒子：直徑分別為3×10-6m、6×10-6m，每種粒子的質量流量相等；在入口邊界上，粒子的速度為1 200m/s，溫度為2 200K，總的質量流量為13.33kg/s·m2，材料密度為3 970kg/m3，熱容為765J/kg·K。表面熱適應系數(shù)設為0.89。計算區(qū)域劃分為5 000個網格，時間步長為1.5×10-9s。流場穩(wěn)定后，約有270 000個氣體仿真分子和10 000個固體仿真顆粒。

圖2、圖3是本文對兩相流中固體顆粒與氣相宏觀流動參數(shù)模擬結果與文獻［10］計算結果的比較情況。對固體顆粒，在入口邊界上速度比氣相速度低，會受到氣體的推動，速度逐漸增加（圖2）；而初始溫度比氣相溫度高，熱能會逐漸向氣相傳遞，溫度逐漸下降。由于固體顆粒的尺寸和質量不一樣，由前面介紹的方法可知，在相同條件下，固體顆粒所受到的力與顆粒半徑的平方成正比，而顆粒質量與顆粒半徑的立方成正比，綜合起來，顆粒的加速度與顆粒的半徑成反比。所以小尺寸的顆粒加速度比較大。在流場中，不同尺寸的顆粒會逐漸分離，其數(shù)密度會逐漸下降。圖2表明顆粒相的速度、數(shù)密度、溫度流向變化基本上屬于線性變化。對氣相流動，受較低速度的固相顆粒的阻擋，氣體的速度沿流動方向是下降的（圖3），數(shù)密度逐漸上升。由于不斷從溫度較高的固相顆粒獲得能量，氣相的溫度沿流動方向逐漸升高。相應地，氣相流動的參數(shù)變化基本上也是線性的。從圖2、圖3顯示出，本文對氣固兩相宏觀流動參數(shù)模擬結果與文獻研究結果吻合很好。

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

圖2 氣固兩相流中固體顆粒宏觀流動參數(shù)計算比較Fig.2 Comparison of parameters of particle in two phase flow

圖3 氣固混合物兩相流中氣體宏觀流動參數(shù)計算比較Fig.3 Comparison of parameters of gas in two phase flow

圖4和圖5繪出流場中不同截面上動量與能量傳遞速率隨位置的變化關系。能量流率和動量流率通過在計算中分別采樣統(tǒng)計通過某個截面的不同相的粒子動量和能量得到。

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

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

圖中顯示出，不論是氣相還是固相，流動中的動量流率和能量流率基本上呈線性變化過程，通過比較氣相和固相的變化梯度，可以看出，兩者的梯度大小大致相等，符號相反。這表明在模擬計算中，氣相作用在固相的動量和能量，等于固相反作用于氣相的動量和能量，即在相間相互作用時，動量和能量是守恒的。在前面介紹的雙向耦合算法中，氣相對固相的作用和固相對氣相的作用過程是解耦的，兩個過程獨立計算。在計算固相對氣相的作用時，其過程是一個隨機統(tǒng)計計算過程，不能保證在每一步的計算中嚴格地遵守動量和能量守恒，但是在大樣本數(shù)情況下統(tǒng)計平均時，可以保證動量和能量守恒。

通過以上的對比模擬，表明本項目所建立的氣-固兩相流動DSMC雙向耦合算法用于計算仿真固體火箭發(fā)動機羽流中氣相混合物與固體顆粒輸運過程是可靠的。

2.2 發(fā)動機兩相羽流流場計算

采用本文的方法，對某固體發(fā)動機兩相羽流流場進行了模擬。羽流中氣相組元的摩爾分數(shù)見表1。

表1 羽流組元摩爾分數(shù)Table 1 Gas Mole fraction of plume

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

式中，為顆粒尺寸的平均半徑；σ為標準偏差，表示顆粒尺寸分布的分散程度。

本文采用有限分組方法，把Al2O3固體顆粒分成5組不同尺寸和質量的顆粒，作為5種顆粒來處理，5組顆粒數(shù)密度相等。每一組顆粒具有相同的尺寸和質量，顆粒直徑在20μm～200μm之間。具體尺寸見表2。

表2 固體顆粒尺寸Table 2 Size of particles

發(fā)動機出口馬赫數(shù)為1.5，出口溫度為2 960K。出口速度設為均勻流動，其中固體顆粒的質量流量占總流量的30%。計算中，環(huán)境設為真空邊界條件。圖6是兩相羽流的流場數(shù)密度分布。從氣相流場看，由于噴流出口馬赫數(shù)較低，只比聲速稍高，噴流向真空膨脹十分迅速，數(shù)密度和溫度均迅速下降。固相顆粒的擴散是氣相擴散引起的，圖中顯示，固相顆粒的擴散是有界的。當氣相膨脹到密度很低時，對固相顆粒的影響很小，顆粒會保持自己的運動狀態(tài)做勻速直線運動。從數(shù)密度分布來看，固體顆粒比較集中地分布在核心區(qū)域，這表明氣流對固體顆粒的擴散作用是有限的。

溫度的分布比較復雜（圖7），沿軸向，隨著氣相的迅速膨脹，溫度迅速下降，能量從固體顆粒向氣相轉化，顆粒的溫度也逐漸下降。沿徑向，軸線附近的溫度比外圍的溫度稍低，更外圍的顆粒溫度更低，這是因為5種顆粒與氣相相互作用時能量交換不同（參見圖8），導致顆粒的溫度不同；同時動量交換也不相同，顆粒擴散的程度因此會有所差異，這些因素導致溫度的分布比較復雜。不過整個流場溫度分布差別不大，最大差異在300K左右。軸向速度表明，外圍顆粒被氣相加速較大，這是因為外圍主要是小尺寸的顆粒，尺寸顆粒越小，越容易被氣相加速。

圖8是噴流軸線上幾種粒子的參數(shù)分布比較。由于氣相膨脹很快，在一定距離之外，氣相密度變得很低，對固體顆粒的影響已經很小。本文的計算中，在＝5（以噴流半徑為參考無量綱化）以內，氣相對固體顆粒有較大的影響，顆粒溫度下降較快，＝5以外，溫度下降較小，基本上保持在一個溫度上。顆粒尺寸越小，溫度變化越大。在＝5以內，顆粒的速度快速增加，在＝5以外，顆粒速度基本上保持在同一個速度上。顆粒尺寸越小，速度變化越大。

圖6 兩相羽流流場數(shù)密度分布Fig.6 Number density distribution of two phase plume

圖7 羽流中固體顆粒溫度分布Fig.7 Temperature distribution of particle in plume

圖8 不同固體顆粒沿軸線參數(shù)分布Fig.8 Particle′s parameter distribution along axis

3 結束語

本文在DSMC方法的基礎上，采用雙向耦合技術，建立了適用于DSMC方法的兩相流的數(shù)值模擬方法，通過相間作用的解耦處理，實現(xiàn)了氣固兩相動量和能量相互作用的模擬。通過對有關算例的計算表明，本文建立的方法雖然在每一個時間步長內不能保證相間動量和能量交換時守恒，但是在大樣本統(tǒng)計平均條件下，能夠保證相間動量和能量交換時守恒。同時模擬也表明，固體顆粒尺寸、氣體稀薄程度等條件對相間的相互作用有很大的影響。本文的方法為稀薄兩相流的計算提供了一種數(shù)值模擬手段。

［1］ Geisler R L.A global view of the use of aluminum fuel in solid rocket motors［R］.AIAA 2002-3748.38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference.Indianapolis，2002.

［2］ Dettleff G.Plume flow and plume impingement in space technology［J］.Progress in Aerospace Sciences，1991，28（1）：1-71.

［3］ Simmons S.Rocket exhaust plume phenomenology［M］.EI Segundo，California：The Aerospace Press，2000.

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

［5］ Gimelshein S F，Boyd I D，Ivanov M S.Modeling of internal energy transfer in plume flows of polyatomic［R］.AIAA 99-0738，1999.

［6］ 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.

［7］ Huang Lin，Nie Wansheng，Chen Weifang.Studying of multiplume interference effects for attitude-control thruster with DSMC method［J］.Acta Aerodynimca Sinica，2003，21（1）：104-108.（in Chinese）黃琳，聶萬勝，陳偉芳.姿控發(fā)動機高空羽流流場干擾效應的DSMC方法研究［J］.空氣動力學學報，2003，21（1）：104-108.

［8］ Gallis M A，Rader D J，Torczynski J R.DSMC simulations of the thermophoretic force on a spherical macroscopic particle［R］.AIAA 2001-2890.

［9］ Gallis M A，Torczynski J R，Rader D J.An approach for simulating the transport of spherical particles in a rarefied gas flow via the direct simulation Monte Carlo method［J］.Physics of Fluids，2001，13（11）：3482-3492.

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

［11］Burt J M，Boyd I D.Monte Carlo simulation of a rarefied multiphase plume flow［C］.43rdAIAA Aerospace Sciences Meeting and Exhibit.Reno，Nevada，2005，AIAA 2005-964.

［12］Sergey F Gimelshein，Alina A Alexeenko，Dean C Wadsworth，et al.The influence of particulates on thruster plume/shock layer interaction at high altitudes［C］.43rdAIAA Aerospace Sciences Meeting and Exhibit，Reno，Nevada，2005.AIAA 2005-766.

［13］Burt J M，Boyd I D.Particle rotation effects in rarefied twophase plume flows［C］.24thInternational Symposium on Rarefied Gas Dynamics.Monopoli，Italy，2004.

［14］Gosse S，Sarou Kanian V，Veron E，et al.Characterization and morphology of alumina particles in solid propellant subscale rocket motor plumes［C］.36thAIAA Thermophysics Conference.Orlando，F(xiàn)L，2003.AIAA 2003-3649.

［15］Li Jie，Ren Bin，Chen Weifang.Modeling and numerical simulation of gas-solid jet in the transitional-rarefied regime［J］.Acta Aerodynamica Sinica，2005，23（4）：484-489.（in Chinese）李潔，任兵，陳偉芳.稀薄流過渡區(qū)氣固兩相噴流的建模與數(shù)值模擬［J］.空氣動力學學報，2005，23（4）：484-489.

A numerical method for simulating rarefied two phase flow

Li Zhonghua，Li Zhihui，Peng Aoping，Wu Junlin
（HypervelocityAerodynamicsInstituteofChinaAerodynamicsResearchandDevelopmentCenter，Mianyang621000，China）

A numerical simulation approach is presented to simulate rarefied two phase flow by applying two-way coupling technical in DSMC method.The interaction between rarefied gas and solid particles is dealt with decoupled approximation to compute momentum and energy exchange between phases.Results of simulation examples validate the methods and the models by the comparison of the present results with that of typical confirmed results from references，and indicate that it can ensure momentum and energy conservation between phases even in case of large sample number.A two-phase plume flow field from solid rocket engine is simulated using this present method.Away from the engine nozzle，gas and solid particles have a greater impact caused by temperature drop and increases of the speed of particles.The particles temperature away from the nozzle exit tends to be balanced.Simulation shows that the particle size and gas rarefaction would affect the interaction evidently.

rarefied gas；DSMC method；two phase flow；two-way coupled method；plume

V211.3

：Adoi：10.7638/kqdlxxb-2013.0028

0258-1825（2015）02-0266-06

2013-3-6；

：2013-5-8

國家自然科學基金項目（91016027、11325212）；國家重點基礎研究發(fā)展計劃（2014CB744100）

李中華＊（1970-），男，河南商丘人，高級工程師，主要從事稀薄氣體動力學研究.E-mail：lzhzcb＠aliyun.com

李中華，李志輝，彭傲平，等.一種稀薄兩相流動的數(shù)值模擬方法［J］.空氣動力學學報，2015，33（2）：266-271.

10.7638/kqdlxxb-2013.0028 Li Z H，Li Z H，Peng A P，et al.A numerical method for simulating rarefied two phase flow［J］.Acta Aerodynamica Sinica，2015，33（2）：266-271.