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

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

（國防科技大學 航天與材料工程學院，湖南 長沙410073）

0 引 言

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

1 氣粒雙向耦合作用的熱力學模型

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

式中Fnum表示氣體分子權因子，即一個模擬分子所代表的真實分子個數。V是網格體積。τ是顆粒表面熱適應系數。Rp是顆粒半徑。Tp是顆粒溫度。m是氣體分子質量。ur是氣體分子相對于顆粒的速度，Cr是ur的模值。ζ是氣體分子的轉動能自由度。erot是單個氣體分子的轉動能。

在一個時間步長Δt內，每個網格中依據公式（1）和（2）累加計算出所有氣體模擬分子對單個顆粒的力和熱作用，逐個求出網格內每一顆粒所受到的相間作用力和熱傳遞，累計網格內所有顆粒的作用力和熱傳遞，得到網格內氣相所受到的總相間力和熱傳遞。

接下來，將網格總的氣粒相互作用體現在單個的氣相模擬分子上。若不考慮化學反應等傳質現象，氣粒相互作用后，兩相的宏觀速度和溫度將發生變化，氣相的宏觀速度是氣體模擬分子的平均速度，而氣相的溫度是氣體模擬分子熱運動速度和內能的度量，從而氣粒相互作用改變了氣體模擬分子的平均速度、熱運動速度和內能，并且其改變量要滿足相間力和熱的傳遞公式［11］。

首先，氣粒相間相互作用傳遞動量。由于動量傳遞使得單個網格中氣相分子的平均速度發生變化，而動量傳遞并未改變氣相分子能量分布，所以氣相分子的平動能不變。再者，氣粒相間相互作用傳遞能量。從宏觀上講，能量傳遞改變了氣相的溫度分布。因為氣相的溫度由氣相模擬分子的平動溫度、轉動溫度和振動溫度構成，其中平動溫度是分子熱運動速度的度量，所以相間能量傳遞改變了單個網格中氣相分子的熱運動速度和內能，而能量傳遞并未改變氣相分子平均速度，故氣相分子的平均速度在能量傳遞過程中保持不變［11］。

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

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

由于稀薄流中顆粒的直徑通常是微米的量級，這種細微顆粒的碰撞若采用傳統的決定論方法，不僅導致記錄參數數目的增大而且對網格的分辨率要求高，這需要巨大的計算機內存和機時，以致難以實現數值模擬，因此有必要采用隨機論觀點，通過顆粒的碰撞幾率函數來確定發生碰撞的顆粒，根據顆粒碰撞的力學機理和散射模型確定碰撞后各變量值，以此構造細微顆粒的碰撞模型［10］。

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

固態顆粒可看作硬球，若不計表面摩擦和轉動的影響，光滑圓球顆粒的碰撞遵循質量和動量守恒定律，碰撞前后的質心速度不變。同時，考慮顆粒間碰撞的非彈性，在碰撞過程中有能量損失，記e為顆粒間碰撞的彈性恢復系數，碰撞前后兩顆粒的相對運動速度為

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

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

式中：R1、R2為隨機數。

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

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

Ashgriz和Poo［15］通過理論分析建立了合并－反射分離的邊界模型，經計算與試驗結果的對比，定義了發生反射分離的臨界韋伯數。Benson等［16］定義了摩擦分離臨界韋伯數。

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

其中Np為落在網格內的液滴樣本的數目。Fnump為液滴的權因子，即一個液滴樣本所代表的真實液滴的個數。Vc為網格的體積。cr是兩液滴相對速度。σT是兩液滴碰撞截面：

這里，下標p1和p2代表兩個液滴樣本。

從液滴樣本中采用置換方法隨機選取Mcol個碰撞對，每個碰撞對發生碰撞的頻率為：

首先從（0－1）間產生隨機數XX，如果XX≥Pcollision，則液滴發生碰撞，其碰撞特性可以通過計算臨界碰撞參數bcr得到。

其中

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

式中Zz＝

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

3 算例和結果分析

計算區域為（0，0）、（0，0.025m）、（0.2m，0）、（0.2m，0.025m）四點構成的矩形域。噴流中心點坐標為（0，0），噴口半徑為4mm。噴流周圍是真空環境。由于整個流場關于x軸呈對稱分布，故只計算半個流場。

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

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

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

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

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

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

通過上述算例可見，本文構造的適用于DSMC算法的固態和液態顆粒碰撞、聚合和分離模型能夠較好地反映稀薄流過渡區氣粒多相復雜流動問題。

4 結 論

本文在氣粒雙向耦合作用的熱力學模型基礎上，建立適用于DSMC算法的固態和液態顆粒碰撞、聚合和分離模型，發展稀薄條件下雙向耦合作用的氣粒多相流的直接模擬Monte Carlo（DSMC）算法，在此基礎上實現高超聲速稀薄流環境中的氣粒多相噴流流場數值模擬。

通過二維真空噴流算例計算結果表明，顆粒聚集在噴流中心線附近逐漸向外稀疏，呈扇形分布，脈動速度越大，顆粒擴散越顯見。液態顆粒聚合更多發生在離噴口稍遠的地方，顆粒脈動速度越大，碰撞聚合越顯見。

由于本文涉及DSMC算法和顆粒動力學模型的結合和數值處理，算法建模的可靠性需要比較驗證，而受目前實驗條件的限制，國內外多相稀薄流的實驗數據極少，有待今后同相關實驗結果進行比較。在現有顆粒動力學研究成果基礎上，本文初步實現了顆粒碰撞、聚合和分離行為的Monte Carlo模擬，為解決稀薄流過渡區氣粒多相復雜流動問題提供了一種新的研究思路。

［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］.力學學報，2004，36（2）：129－139.

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

［11］李潔，任兵，陳偉芳.稀薄流過渡區氣固兩相噴流的建模與數值模擬［J］.空氣動力學學報.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.