高超声速稀薄流的气粒多相流动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.