高超聲速稀薄流的氣粒多相流動(dòng)DSMC算法建模研究

李 潔，石于中，徐振富，王小虎

（國(guó)防科技大學(xué) 航天與材料工程學(xué)院，湖南 長(zhǎng)沙410073）

0 引 言

在高超聲速稀薄流的氣體－顆粒多相流動(dòng)中，由于流場(chǎng)的稀薄氣體效應(yīng)顯著，連續(xù)介質(zhì)假定受到限制，傳統(tǒng)的基于Navier－Stokes方程的CFD方法將不再適用［1］。隨著航天飛行器在高空的機(jī)動(dòng)飛行和制動(dòng)的需求發(fā)展，基于高空氣體－顆粒耦合流場(chǎng)的多相稀薄流動(dòng)研究日益受到重視。國(guó)外早期Epstein［2］和Baines等［3］學(xué)者基于 Maxwellian或 Chapman－Enskog分子速度分布函數(shù)，研究了稀薄氣體中的單個(gè)球形顆粒所受的作用力和熱。這些理論研究工作為開展多相稀薄流的數(shù)值模擬奠定了良好的基礎(chǔ)。通常稀薄流動(dòng)中的顆粒直徑較小，一般是微米量級(jí)甚至更?。?］，類似于氣體分子，顆粒具有質(zhì)量、體積、位置、速度和內(nèi)能，若將顆?？醋鞒赡撤N類型的模擬分子，則可以將成功運(yùn)用于解決稀薄流動(dòng)問題的DSMC方法［5］延伸到模擬顆粒的運(yùn)動(dòng)。鑒于此想法，在前述學(xué)者的工作基礎(chǔ)上，Gallis等［6］利用Green函數(shù)，建立了適用于DSMC方法求解的在不同分子速度分布的稀薄氣相流場(chǎng)中顆粒所受到的作用力和傳熱模型。同樣，Zhang J等［7］利用DSMC方法研究Io星體上的火山羽流流動(dòng)。上述研究工作均是假定顆粒是稀疏的，忽略顆粒的布朗運(yùn)動(dòng)，并且只考慮氣相對(duì)顆粒相的作用，忽略顆粒相對(duì)氣相的影響。此后，為克服上述不足，Burt［8］和Gimelshein［8］利用Gallis［6］推導(dǎo)的單個(gè)球形顆粒熱力學(xué)模型，基于分子動(dòng)力學(xué)理論，構(gòu)造了氣粒雙向耦合作用的熱力學(xué)模型，為多相稀薄流建立了較好的數(shù)學(xué)描述方法，但這些方法存在某些缺陷，主要是由于氣體分子和顆粒的相間作用處理模式各不相同，導(dǎo)致動(dòng)量和能量的守恒僅是從時(shí)間平均層面上實(shí)現(xiàn)，不能保證單個(gè)時(shí)間步長(zhǎng)內(nèi)的氣粒相間碰撞的動(dòng)量守恒和能量守恒，誤差將會(huì)隨著氣體分子和顆粒的數(shù)密度差距增大而增大［1］。國(guó)內(nèi)對(duì)于氣粒多相稀薄流的研究工作起步較晚，樊菁等［9］圍繞空間氣液、氣固多相羽流場(chǎng)展開研究，由于局限于空間真空環(huán)境下，未考慮氣粒兩相的雙向耦合作用。以上國(guó)內(nèi)外的研究工作都未考慮液態(tài)顆粒間的碰撞，因而忽略了氣粒多相稀薄流中的顆粒碰撞、合并現(xiàn)象，而這些現(xiàn)象廣泛存在于多相流中，對(duì)流場(chǎng)的流動(dòng)特性產(chǎn)生一定影響。為此，本文在文獻(xiàn)［10－11］工作基礎(chǔ)上，構(gòu)造了適用于DSMC算法的固態(tài)和液態(tài)顆粒碰撞、聚合和分離模型，發(fā)展稀薄條件下雙向耦合作用的氣粒多相流的直接模擬 Monte Carlo（DSMC）算法，在此基礎(chǔ)上實(shí)現(xiàn)高超聲速稀薄流環(huán)境中的氣粒多相噴流流場(chǎng)數(shù)值模擬，嘗試為解決稀薄流過渡區(qū)氣粒多相復(fù)雜流動(dòng)問題提供一種新的研究思路。

1 氣粒雙向耦合作用的熱力學(xué)模型

利用Gallis［6］的單個(gè)氣體模擬分子對(duì)單個(gè)顆粒的力和熱作用公式［8］：

式中Fnum表示氣體分子權(quán)因子，即一個(gè)模擬分子所代表的真實(shí)分子個(gè)數(shù)。V是網(wǎng)格體積。τ是顆粒表面熱適應(yīng)系數(shù)。Rp是顆粒半徑。Tp是顆粒溫度。m是氣體分子質(zhì)量。ur是氣體分子相對(duì)于顆粒的速度，Cr是ur的模值。ζ是氣體分子的轉(zhuǎn)動(dòng)能自由度。erot是單個(gè)氣體分子的轉(zhuǎn)動(dòng)能。

在一個(gè)時(shí)間步長(zhǎng)Δt內(nèi)，每個(gè)網(wǎng)格中依據(jù)公式（1）和（2）累加計(jì)算出所有氣體模擬分子對(duì)單個(gè)顆粒的力和熱作用，逐個(gè)求出網(wǎng)格內(nèi)每一顆粒所受到的相間作用力和熱傳遞，累計(jì)網(wǎng)格內(nèi)所有顆粒的作用力和熱傳遞，得到網(wǎng)格內(nèi)氣相所受到的總相間力和熱傳遞。

接下來，將網(wǎng)格總的氣粒相互作用體現(xiàn)在單個(gè)的氣相模擬分子上。若不考慮化學(xué)反應(yīng)等傳質(zhì)現(xiàn)象，氣粒相互作用后，兩相的宏觀速度和溫度將發(fā)生變化，氣相的宏觀速度是氣體模擬分子的平均速度，而氣相的溫度是氣體模擬分子熱運(yùn)動(dòng)速度和內(nèi)能的度量，從而氣粒相互作用改變了氣體模擬分子的平均速度、熱運(yùn)動(dòng)速度和內(nèi)能，并且其改變量要滿足相間力和熱的傳遞公式［11］。

首先，氣粒相間相互作用傳遞動(dòng)量。由于動(dòng)量傳遞使得單個(gè)網(wǎng)格中氣相分子的平均速度發(fā)生變化，而動(dòng)量傳遞并未改變氣相分子能量分布，所以氣相分子的平動(dòng)能不變。再者，氣粒相間相互作用傳遞能量。從宏觀上講，能量傳遞改變了氣相的溫度分布。因?yàn)闅庀嗟臏囟扔蓺庀嗄M分子的平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度構(gòu)成，其中平動(dòng)溫度是分子熱運(yùn)動(dòng)速度的度量，所以相間能量傳遞改變了單個(gè)網(wǎng)格中氣相分子的熱運(yùn)動(dòng)速度和內(nèi)能，而能量傳遞并未改變氣相分子平均速度，故氣相分子的平均速度在能量傳遞過程中保持不變［11］。

因此，根據(jù)相間動(dòng)量和能量傳遞守恒原則和氣體分子能量按自由度均布原則，可以構(gòu)造出相間作用后的氣體和顆粒速度、溫度等物理量分布模型［11］。

2 顆粒碰撞、聚合和分離模型

由于稀薄流中顆粒的直徑通常是微米的量級(jí)，這種細(xì)微顆粒的碰撞若采用傳統(tǒng)的決定論方法，不僅導(dǎo)致記錄參數(shù)數(shù)目的增大而且對(duì)網(wǎng)格的分辨率要求高，這需要巨大的計(jì)算機(jī)內(nèi)存和機(jī)時(shí)，以致難以實(shí)現(xiàn)數(shù)值模擬，因此有必要采用隨機(jī)論觀點(diǎn)，通過顆粒的碰撞幾率函數(shù)來確定發(fā)生碰撞的顆粒，根據(jù)顆粒碰撞的力學(xué)機(jī)理和散射模型確定碰撞后各變量值，以此構(gòu)造細(xì)微顆粒的碰撞模型［10］。

由于采用隨機(jī)論觀點(diǎn)，顆粒可以看作成某種類型的模擬粒子，滿足混沌假設(shè)，且顆粒間碰撞是在瞬時(shí)完成的。類似于氣體分子，顆粒具有質(zhì)量、直徑、位置、速度和內(nèi)能，滿足粒子的輸運(yùn)特性，從而可以構(gòu)造出顆粒的碰撞幾率函數(shù)，以此判斷顆粒發(fā)生碰撞的可能性。然后采用Bird的無時(shí)間記數(shù)器方法，選取合適的碰撞對(duì)。根據(jù)碰撞動(dòng)力學(xué)的動(dòng)量守恒和能量守恒定律得出碰撞后各參數(shù)值。由于固態(tài)顆粒和液態(tài)顆粒表面物性不相同，導(dǎo)致碰撞參數(shù)的處理方式不相同，需區(qū)分對(duì)待。

固態(tài)顆?？煽醋饔睬颍舨挥?jì)表面摩擦和轉(zhuǎn)動(dòng)的影響，光滑圓球顆粒的碰撞遵循質(zhì)量和動(dòng)量守恒定律，碰撞前后的質(zhì)心速度不變。同時(shí)，考慮顆粒間碰撞的非彈性，在碰撞過程中有能量損失，記e為顆粒間碰撞的彈性恢復(fù)系數(shù)，碰撞前后兩顆粒的相對(duì)運(yùn)動(dòng)速度為

然后通過硬球顆粒散射模型確定相對(duì)速度方向，假定顆粒的尺寸較小，顆粒的稠密程度不高，可以忽略顆粒的體積對(duì)碰撞幾率的影響，利用球?qū)ΨQ散射模型，類似于氣體分子可確定出g′的方向。記n＝（n1，n2，n3）為碰撞后顆粒相對(duì)運(yùn)動(dòng)速度g′的方向單位矢量，由于碰撞的隨機(jī)性，n為

式中：χ1和χ2為Euler角，由隨機(jī)抽樣確定。其抽樣方法為

式中：R1、R2為隨機(jī)數(shù)。

以上是固體顆粒的碰撞模型。對(duì)于液態(tài)顆粒而言，由于液態(tài)顆粒存在表面張力，碰撞后顆粒運(yùn)動(dòng)形態(tài)將不同于固態(tài)顆粒。根據(jù)文獻(xiàn)［12］～文獻(xiàn)［15］，液態(tài)顆粒碰撞后將發(fā)生合并、摩擦分離（streching separation）或反射分離（reflexive separation）。在摩擦分離過程中，兩個(gè)碰撞粒子重疊區(qū)域小，顆粒起初合并，合并后的液滴由于所具有的內(nèi)部動(dòng)能過量，液滴不能穩(wěn)定存在，在經(jīng)歷短暫的合并后又發(fā)生分離。液態(tài)顆粒碰撞后究竟發(fā)生合并還是分離取決于以下三個(gè)因素：顆粒直徑比值Δ、韋伯?dāng)?shù)We和兩顆粒相撞時(shí)的無量綱偏心距χ。其中，顆粒直徑比值Δ定義為碰撞對(duì)中的小顆粒直徑與大顆粒直徑比值。碰撞韋伯?dāng)?shù)We定義如下：

式中di是小顆粒直徑，v是碰撞顆粒的相對(duì)速度，σd是液滴表面張力。

Ashgriz和Poo［15］通過理論分析建立了合并－反射分離的邊界模型，經(jīng)計(jì)算與試驗(yàn)結(jié)果的對(duì)比，定義了發(fā)生反射分離的臨界韋伯?dāng)?shù)。Benson等［16］定義了摩擦分離臨界韋伯?dāng)?shù)。

本文采用Schmidt D P和Rutland C J提出的NTC（no－time－counter）模型［17］，類似于氣體分子，僅當(dāng)位于同一網(wǎng)格內(nèi)的液滴樣本才可能發(fā)生碰撞，發(fā)生碰撞的次數(shù)Mcol定義為：

其中Np為落在網(wǎng)格內(nèi)的液滴樣本的數(shù)目。Fnump為液滴的權(quán)因子，即一個(gè)液滴樣本所代表的真實(shí)液滴的個(gè)數(shù)。Vc為網(wǎng)格的體積。cr是兩液滴相對(duì)速度。σT是兩液滴碰撞截面：

這里，下標(biāo)p1和p2代表兩個(gè)液滴樣本。

從液滴樣本中采用置換方法隨機(jī)選取Mcol個(gè)碰撞對(duì)，每個(gè)碰撞對(duì)發(fā)生碰撞的頻率為：

首先從（0－1）間產(chǎn)生隨機(jī)數(shù)XX，如果XX≥Pcollision，則液滴發(fā)生碰撞，其碰撞特性可以通過計(jì)算臨界碰撞參數(shù)bcr得到。

其中

其次從（0－1）間產(chǎn)生隨機(jī)數(shù)YY，計(jì)算碰撞參數(shù)b＝（rsmall＋rbig），如果b大于臨界碰撞參數(shù)bcr，那么兩液滴樣本中的液滴發(fā)生摩擦分離。對(duì)于這種碰撞，兩液滴樣本內(nèi)的液滴個(gè)數(shù)不發(fā)生變化，僅僅是液滴的速度發(fā)生變化。碰撞后兩液滴速度分別為：

式中Zz＝

如果b小于臨界碰撞參數(shù)bcr，那么液滴發(fā)生聚合。發(fā)生聚合后的液滴參數(shù)如下：

3 算例和結(jié)果分析

計(jì)算區(qū)域?yàn)椋?，0）、（0，0.025m）、（0.2m，0）、（0.2m，0.025m）四點(diǎn)構(gòu)成的矩形域。噴流中心點(diǎn)坐標(biāo)為（0，0），噴口半徑為4mm。噴流周圍是真空環(huán)境。由于整個(gè)流場(chǎng)關(guān)于x軸呈對(duì)稱分布，故只計(jì)算半個(gè)流場(chǎng)。

噴流為氣粒兩相流。噴流中氣體的組分為H2、N2和CO，摩爾比例分別為0.38、0.31和0.31。氣相速度為3100m／s，密度為0.00011kg／m3，溫度分兩種情況考慮，分別為1000K和1700K。噴流中顆粒由液態(tài)顆粒和固態(tài)顆粒組成，顆粒的材料密度為3970kg／m3，平均脈動(dòng)速度分為25m／s和50m／s，顆粒直徑分布分兩組：0.3μm～1.0μm和0.03μm～0.1μm，即兩組之間差一個(gè)量級(jí)，不同直徑的顆粒物性參數(shù)見表1。

表1 噴流出口顆粒相參數(shù)分布Table 1 Particle properties at the nozzle exit

圖1描述的是粒徑為0.03μm和0.1μm的顆粒在不同脈動(dòng)速度下的位置分布圖。由圖可見，顆粒聚集在中心線附近逐漸向外稀疏，呈扇形分布，脈動(dòng)速度越大，顆粒擴(kuò)散越顯見。

圖1 不同顆粒脈動(dòng)速度下的顆粒位置分布圖（Tg＝1700K）Fig.1 Location of particles in different peculiar particle velocities

圖2描述的是顆粒發(fā)生碰撞、聚合和分離后，粒徑dp在不同范圍內(nèi)顆粒位置分布圖。不難發(fā)現(xiàn)，液態(tài)顆粒聚合更多發(fā)生在離噴口遠(yuǎn)些的地方，顆粒脈動(dòng)速度越大，碰撞聚合越顯見。由式（13）可知，液體顆粒發(fā)生碰撞、聚合和分離主要受三個(gè)因素的影響：碰撞韋伯?dāng)?shù)、無量綱偏心距和液滴半徑比。隨著噴口噴出的顆粒遠(yuǎn)離噴口，不同粒徑顆粒間的相對(duì)速度進(jìn)一步加大，引起韋伯?dāng)?shù)的增大，加劇碰撞聚合的發(fā)生。由式（21）可知，圖2（a）和2（b）粒徑介于0.03μm 與0.04μm之間的顆粒是由兩個(gè)粒徑均為0.03μm的顆粒碰撞聚合而成。同理，圖2（c）和2（d）所示顆粒是由兩種碰撞對(duì)產(chǎn)生：粒徑為0.03μm 的顆粒與0.04μm的顆粒碰撞、兩個(gè)粒徑皆為0.04μm的顆粒碰撞，而圖2（e）和2（f）所示顆粒由三種碰撞對(duì)產(chǎn)生：粒徑為0.03μm顆粒與0.06μm的顆粒碰撞、粒徑為0.04μm的顆粒與0.06μm的顆粒碰撞、兩個(gè)粒徑均為0.06μm的顆粒碰撞。故在圖2（e）和2（f）中，由更多種碰撞對(duì)產(chǎn)生的顆粒數(shù)相應(yīng)更多。

圖2 不同粒徑液態(tài)顆粒的位置分布圖（Tg＝1700K）Fig.2 Location of liquid particles with variable diameters

通過上述算例可見，本文構(gòu)造的適用于DSMC算法的固態(tài)和液態(tài)顆粒碰撞、聚合和分離模型能夠較好地反映稀薄流過渡區(qū)氣粒多相復(fù)雜流動(dòng)問題。

4 結(jié) 論

本文在氣粒雙向耦合作用的熱力學(xué)模型基礎(chǔ)上，建立適用于DSMC算法的固態(tài)和液態(tài)顆粒碰撞、聚合和分離模型，發(fā)展稀薄條件下雙向耦合作用的氣粒多相流的直接模擬Monte Carlo（DSMC）算法，在此基礎(chǔ)上實(shí)現(xiàn)高超聲速稀薄流環(huán)境中的氣粒多相噴流流場(chǎng)數(shù)值模擬。

通過二維真空噴流算例計(jì)算結(jié)果表明，顆粒聚集在噴流中心線附近逐漸向外稀疏，呈扇形分布，脈動(dòng)速度越大，顆粒擴(kuò)散越顯見。液態(tài)顆粒聚合更多發(fā)生在離噴口稍遠(yuǎn)的地方，顆粒脈動(dòng)速度越大，碰撞聚合越顯見。

由于本文涉及DSMC算法和顆粒動(dòng)力學(xué)模型的結(jié)合和數(shù)值處理，算法建模的可靠性需要比較驗(yàn)證，而受目前實(shí)驗(yàn)條件的限制，國(guó)內(nèi)外多相稀薄流的實(shí)驗(yàn)數(shù)據(jù)極少，有待今后同相關(guān)實(shí)驗(yàn)結(jié)果進(jìn)行比較。在現(xiàn)有顆粒動(dòng)力學(xué)研究成果基礎(chǔ)上，本文初步實(shí)現(xiàn)了顆粒碰撞、聚合和分離行為的Monte Carlo模擬，為解決稀薄流過渡區(qū)氣粒多相復(fù)雜流動(dòng)問題提供了一種新的研究思路。

［1］BURT J M，BOYD I D.Development of a two－way coupled model for two phase rarefied plume flows［R］.AIAA－2004－1351.

［2］EPSTEIN P S.On the resistance experienced by spheres in their motion through gases［J］.PhysRev，1924，24：710－733.

［3］BAINES M J，WILLIAMS I P，ASEBIOMO A S.Resistance to the motion of a small sphere moving through agas［J］.MonNotRAstronSoc，1965，130：63－74.

［4］BURT J M，BOYD I D.A Monte Carlo radiation model for simulating rarefied multiphase plume flows［R］.AIAA－2005－4691.

［5］BiRd G A.Molecular gas dynamics and the direct simulation of gas flows［M］.Clarendon Press，Oxford，1994.

［6］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］.PhysicsofFluids，2001，13（11）：3482－3492.

［7］ZHANG J，GOLDSTEIN D B，GIMELSHEIN N E，GIMELSHEIN S F，LEVIN D A，VARGESE P L.Modeling low density sulfur dioxide jets：application to volcanoes on Jupiter＇s moon Io［R］.AIAA－2001－2767.

［8］BURT J M，BOYD I D.Development of a two－way coupled model for two phase rarefied plume flows［R］.AIAA－2004－1351.

［9］樊菁，劉宏立，蔣建政，等.火箭剩余推進(jìn)劑排放過程的分析與模擬［J］.力學(xué)學(xué)報(bào)，2004，36（2）：129－139.

［10］李潔，尹樂，顏力.稀薄流氣粒兩相耦合作用的熱力學(xué)模型［J］.國(guó)防科技大學(xué)學(xué)報(bào).2009，31（3）：6－10.

［11］李潔，任兵，陳偉芳.稀薄流過渡區(qū)氣固兩相噴流的建模與數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào).2005，23（4）：484－489.

［12］WILLIS K D，ORME M E.Experiments on the dynamics of droplet collision in a vacuum［J］.Exp.Fluids，2000，29 （4）：347－358.

［13］QIAN J，LAW C K.Regimes of coalescene and separation in droplet collision［J］.J.FluidMech.，1997，331：59－80.

［14］ORME M.Experiments on droplet collisions，bounce，coalescene and disruption［J］.Prog.EnergyCombust.Sci，1997，23（1）：65－79.

［15］ASHGRIZ N，POO J Y.Coalescene and separation in binary collision of liquid drops［J］.J.FluidMech.，1990，221：183－204.

［16］BENSON C M，ZHONG J，GIMELSHEIN S F.A general model for the simulation of aerosol droplets in a high－temperature environment［R］.AIAA－2002－3181.

［17］SCHMIDT D P，RUTLAND C J.A new droplet collision algorithm［J］.JournalofComputationalPhysics，2000，64（1）：62－80.