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

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

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

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

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

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

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

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

0 引 言

固體火箭發(fā)動機(jī)為了達(dá)到提高比沖的目的，廣泛采用含金屬復(fù)合推進(jìn)劑。燃燒后生成大量的固體顆粒，燃?xì)馍淞魍ǔＪ怯啥喾N氣體成分、固體顆粒形成典型的兩相流動，這些顆粒既對燃?xì)馀艑?dǎo)裝置產(chǎn)生巨大的沖刷侵蝕作用，又對燃?xì)馍淞髁鲌霎a(chǎn)生直接影響，形成復(fù)雜的相互作用干擾流場［1-3］。

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

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

1 數(shù)值模擬算法

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

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

考慮固體顆粒對周圍氣體的影響，首先要確定在每個(gè)時(shí)間步長內(nèi)哪個(gè)仿真分子將與顆粒進(jìn)行碰撞。對Bird的非時(shí)間計(jì)數(shù)方法進(jìn)行修正，來確定與所選的顆粒可能發(fā)生碰撞的計(jì)算分子數(shù)ns。

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

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

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

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

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

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

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

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

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

2 算例與討論

2.1 算例驗(yàn)證

為了驗(yàn)證本文建立兩相流DSMC雙向耦合算法可靠性，對文獻(xiàn)［10］的算例進(jìn)行了計(jì)算。計(jì)算區(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，每種粒子的質(zhì)量流量相等；在入口邊界上，粒子的速度為1 200m/s，溫度為2 200K，總的質(zhì)量流量為13.33kg/s·m2，材料密度為3 970kg/m3，熱容為765J/kg·K。表面熱適應(yīng)系數(shù)設(shè)為0.89。計(jì)算區(qū)域劃分為5 000個(gè)網(wǎng)格，時(shí)間步長為1.5×10-9s。流場穩(wěn)定后，約有270 000個(gè)氣體仿真分子和10 000個(gè)固體仿真顆粒。

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

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

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

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

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

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

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

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

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

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

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

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

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

式中，為顆粒尺寸的平均半徑；σ為標(biāo)準(zhǔn)偏差，表示顆粒尺寸分布的分散程度。

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

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

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

溫度的分布比較復(fù)雜（圖7），沿軸向，隨著氣相的迅速膨脹，溫度迅速下降，能量從固體顆粒向氣相轉(zhuǎn)化，顆粒的溫度也逐漸下降。沿徑向，軸線附近的溫度比外圍的溫度稍低，更外圍的顆粒溫度更低，這是因?yàn)?種顆粒與氣相相互作用時(shí)能量交換不同（參見圖8），導(dǎo)致顆粒的溫度不同；同時(shí)動量交換也不相同，顆粒擴(kuò)散的程度因此會有所差異，這些因素導(dǎo)致溫度的分布比較復(fù)雜。不過整個(gè)流場溫度分布差別不大，最大差異在300K左右。軸向速度表明，外圍顆粒被氣相加速較大，這是因?yàn)橥鈬饕切〕叽绲念w粒，尺寸顆粒越小，越容易被氣相加速。

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

圖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 結(jié)束語

本文在DSMC方法的基礎(chǔ)上，采用雙向耦合技術(shù)，建立了適用于DSMC方法的兩相流的數(shù)值模擬方法，通過相間作用的解耦處理，實(shí)現(xiàn)了氣固兩相動量和能量相互作用的模擬。通過對有關(guān)算例的計(jì)算表明，本文建立的方法雖然在每一個(gè)時(shí)間步長內(nèi)不能保證相間動量和能量交換時(shí)守恒，但是在大樣本統(tǒng)計(jì)平均條件下，能夠保證相間動量和能量交換時(shí)守恒。同時(shí)模擬也表明，固體顆粒尺寸、氣體稀薄程度等條件對相間的相互作用有很大的影響。本文的方法為稀薄兩相流的計(jì)算提供了一種數(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ā)動機(jī)高空羽流流場干擾效應(yīng)的DSMC方法研究［J］.空氣動力學(xué)學(xué)報(bào)，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］.空氣動力學(xué)學(xué)報(bào)，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

國家自然科學(xué)基金項(xiàng)目（91016027、11325212）；國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（2014CB744100）

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

李中華，李志輝，彭傲平，等.一種稀薄兩相流動的數(shù)值模擬方法［J］.空氣動力學(xué)學(xué)報(bào)，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.