二维喷射液滴在重力场中的运动模型研究

刘展位,薄涵亮

(清华大学 核能与新能源技术研究院,北京 100084)

欧拉与拉格朗日是流体力学描述的两种方法,拉格朗日方法着眼于描述流体质点的位置随时间变化的规律,迹线即为流体质点在空间运动所描绘的曲线。欧拉方法着眼于空间点,在空间每一点上面描述流体运动随时间变化的情况,采用流线来描述流场。欧拉与拉格朗日方法的解耦问题一直是流体力学研究的热点和难点问题。根据Yang[1]对多相流模拟的总结,多相流模型包括欧拉-欧拉(EE)模型[2]、MUSIG模型、非齐次MUSIG模型[3]。其中EE模型计算代价较小但模拟精度相对比较有限。非齐次MUSIG模型扩展了齐次的MUSIG模型,能够使不同尺寸的离散相具有不同的速度场,主要通过求解多组气相的Navier-Stokes方程,以及改进的连续性方程,从而获得不同的离散相速度场,其缺点是增加了计算成本,且一般速度组为2个,因此对离散相速度组的描述受到一定限制。

基于粒子的方法包括欧拉-拉格朗日方法以及其他粒子法。欧拉-拉格朗日方法包括离散相模型(DPM)[4]、稠密离散相模型(DDPM)[5]、离散元法(DEM)[6]、光滑粒子法(SPH)[7]等。其他粒子法包括格子玻尔兹曼方法(LBM)[8]、耗散粒子动力学(DPD)[9]。本课题组采用欧拉-拉格朗日方法的整体框架,对汽水分离器内液滴行为进行建模,从探究汽水分离器中液滴的微观机理到获取流场宏观参数,从而优化设计汽水分离器结构的研究路线,并通过自编程序与商用CFD软件结合开发出一整套计算方法。这些行为包括:液滴的产生、运动、碰撞、聚合、破碎、蒸发、冷凝和消亡等[10]。欧拉-拉格朗日方法虽然较欧拉-欧拉法更能精确地模拟粒子运动,但是计算量随着粒子数目增加而明显增加。因此,开发了欧拉网格逼近方法[11],欧拉网格逼近算法已被应用于描述小液滴在重力空间中的运动特征,同时该方法也被应用于多液滴的情况并验证了计算的可行性[11]。由于实际蒸汽发生器的汽水分离以及其他的汽水分离装置空间内具有部分大液滴,因此有必要将该方法拓展到大液滴的工况。

本文提出一种根据欧拉网格逼近方法来计算汽水分离器中大液滴运动的方法。介绍大液滴的计算公式和密度定义,验证欧拉网格逼近方法应用于大液滴粒子计算的有效性,最后将欧拉网格逼近方法应用于实际AP1000重力分离空间中大液滴的计算,分析3种尺寸的喷射液滴在重力分离空间中的行为,为汽水分离等装置中的大液滴分析提供有效的计算分析工具。

1 理论分析

1.1 喷射液滴形成机理

气泡在自由表面处破裂主要过程为:首先,当气泡从水室上浮到蒸汽表面时,气泡上半部(液帽)突出汽空间,水膜不断减薄最终破裂形成一些细小的水滴;之后,由于气液两相压力差及流体表面张力的作用,气泡周围的水向中心聚集以填补空虚位置形成1个向上射流。射流以一定速度向上升起且在峰顶断裂生成几个大水滴。因此,可以归纳气泡破裂产生液滴的两种主要类型(图1):第一类气泡液膜厚度变薄后破裂产生的细小液滴被称为膜液滴;第二类由于气泡破裂的一瞬间在自由表面留下空洞,同时由于瑞利不稳定性的作用,周围液体涌入填补空缺位置便产生向上的不稳定瑞利射流,进而在表面张力作用下产生射流断裂形成的液滴叫做喷射液滴[12]。文献[13]的实验已观察到这两类液滴,关于喷射液滴的产生文献[14]已进行了大量的实验研究与数值模拟。

图1 气泡破裂产生液滴的两种类型[12]Fig.1 Two types of bubble rupture to produce droplet[12]

文献[15]对汽水分离器入口膜液滴的尺寸分布进行了实验研究,其尺寸分布如图2所示,而文献[16]等对瑞利不稳定性产生的喷射液滴开展了研究,得到了平均喷射液滴直径为1.5 mm。

图2 蒸汽发生器内膜液滴尺寸的分布[15]Fig.2 Distribution of steam generator inner film droplet size[15]

1.2 基本假设

为了简化计算,假设如下:1) 外流场处于稳态;2) 单个液滴为刚性圆球,其质量和大小都是固定的;3) 不考虑粒子的碰撞、聚并、破碎等作用;4) 液滴和蒸汽之间的传热和传质可以忽略不计;5) 单向耦合假设,认为粒子对于流场的作用力可以忽略,只考虑流场对粒子的作用力。

2 数学描述

2.1 数学模型

图3为网格中大液滴示意图,大液滴的直径大于网格尺寸,一部分网格能被大液滴完全覆盖,而另一部分处于液滴边缘的网格只能被大液滴部分覆盖。认为液滴每步的运动都是质点运动,大液滴覆盖到的所有网格节点的速度均一致,而对于这两种不同的网格需要采用不同的密度定义。其中D为大液滴直径,d为网格长度。

图3 网格中大液滴示意图Fig.3 Schematic diagram of large droplet in grid

假设液滴能完全覆盖的网格节点为Nc,不完全覆盖的网格节点为Ne,可以得到完全覆盖的网格节点所占的面积Ac为:

Ac=Ncd2

(1)

大液滴所占面积Ao为:

(2)

每个网格面积Am为:

Am=d2

(3)

不完全覆盖的网格节点所占面积Ae为:

(4)

不完全覆盖的网格节点的系数k为:

(5)

其中,k的范围为(0,1)。

2.2 密度公式

欧拉网格逼近方法的基本方程参考文献[11],认为在二维情况下完全覆盖的网格节点的密度是一致的,而不完全覆盖的网格节点的密度需要乘以1个修正系数。因此,对于大于网格直径的液滴喷射密度定义如下。

1) 在完全覆盖的网格中

当运动速度|v(i)|≤Df时,则:

ρl(i)=ρf

(6)

当运动速度|v(i)|>Df时,则:

(7)

2) 在非完全覆盖的网格中

当运动速度|v(i)|

ρl(i)=kρf

(8)

当运动速度|v(i)|>Df时,则:

(9)

其中:ρl为网格节点密度;ρf为水密度;md为液滴质量;N为液滴数目;ΔV为网格体积;v为液滴运动速度;f为液滴入射频率。

2.3 方程求解

大液滴的方程求解与小液滴的方程求解类似,主要求解步骤如下。

1) 计算x和y方向路程

对于x方向速度vx和y方向速度vy均为正方向,且xacc≥0、yacc≥0时,则:

Δx(i)=dx-xacc(i)

Δy(i)=dy-yacc(i)

(10)

对于vx为正方向、vy为负方向,且xacc≥0、yacc≤0时,则:

Δx(i)=dx-xacc(i)

Δy(i)=-dy+yacc(i)

(11)

其中,xacc、yacc分别为x方向和y方向的路程。

2) 计算运动时间

(1) 求解x方向方程和y方向方程得到tx和ty

对于vx和vy均为正方向,则:

(12)

其中,ax和ay分别为x和y方向的加速度。

(2) 判断第i步计算的tx(i)与ty(i)是否为虚根,修正tx(i)和ty(i)

当粒子运动到y方向最高点时,可能会出现ty等于虚根的情况,原因在于粒子不能到达y正方向下一个网格节点,而是经历了节点间的上升和下降过程。因此需要采用新的y方向求解方程,以保证ty为实根。对于vx为正方向、vy为正方向,且xacc(i)≥0、yacc(i)≥0时,y方向产生了虚根,此时修正方程为式(13),得到的tx(i)和ty(i)均为实根。

Δx(i)=dx-xacc(i)

Δy(i)=-yacc(i)

(13)

(3) 比较tx和ty,并将其中较小值作为时间参数t

t(i)=min{tx(i),ty(i)}

(14)

3) 更新速度、累积量、密度和位移

时间参数t(i)后的vx和vy为:

vx(i+1)=vx(i)+ax(i)t(i)

vy(i+1)=vy(i)+ay(i)t(i)

(15)

对于时间参数t=ty(i),则:

xacc(i+1)=xacc(i)+vx(i)t(i)+

0.5ax(i)t(i)2

yacc(i+1)=0

(16)

对于时间参数t=tx(i),则:

yacc(i+1)=yacc(i)+vy(i)t(i)+

0.5ay(i)t(i)2

xacc(i+1)=0

(17)

4) 更新液滴覆盖的节点密度和液滴覆盖的节点速度

大液滴运动时可能会出现网格节点重叠的问题,需要比较此步节点密度和原有节点密度,并取其中的较大值作为节点密度。其中,ρn为此步大液滴的节点密度,ρo为节点原密度。同时,如果节点密度发生改变,节点速度也改变为对应的值。

ρl(i+1)=max(ρn(i+1),ρo)

(18)

3 模型应用

3.1 模型验证与比较

为了展现模型在工程应用上的意义,本文以AP1000蒸汽发生器的汽水分离系统为研究对象。AP1000汽水分离器的压力p=5.76 MPa,外流场为纯重力场且蒸汽流速v=1.45 m/s。喷射液滴的速度、半径、高度信息由文献[16]得到。

首先根据压力确定临界特征液滴半径,液滴受到的力包括重力、浮力、曳力、附加质量力、Magnus力、Saffman力[17]。图4为不同尺寸喷射液滴的网格节点示意图,其中黑色节点代表完全被液滴覆盖的节点,红色节点代表不完全被液滴覆盖的节点。由图4可知,液滴尺寸越大,其覆盖的黑色节点越多,且总的覆盖节点也更多。液滴的受力达到平衡时,即得到悬浮液滴的半径r=1 356 μm,此时液滴为临界尺寸液滴,临界尺寸液滴减速为0且悬浮在某一位置。

液滴尺寸:a——临界尺寸(r=1 356 μm);b——小于临界尺寸(r=1 250 μm);c——大于临界尺寸 (r=1 500 μm)图4 网格节点示意图Fig.4 Schematic diagram of grid node

分别对3种类型尺寸的液滴进行了对比验证,两种方法的速度相对误差达到0.1%(图5),密度相对误差达到0.6%(图6)。

液滴尺寸:a——临界尺寸(r=1 356 μm);b——小于临界尺寸(r=1 250 μm);c——大于临界尺寸(r=1 500 μm)图5 欧拉网格逼近方法和欧拉-拉格朗日方法的速度对比Fig.5 Comparison of velocity between Euler grid approximation and Euler-Lagrange methods

液滴尺寸:a——临界尺寸(r=1 356 μm);b——小于临界尺寸(r=1 250 μm);c——大于临界尺寸(r=1 500 μm)图6 欧拉网格逼近方法和欧拉-拉格朗日方法的密度对比Fig.6 Comparison of density between Euler grid approximation and Euler-Lagrange methods

在计算效率方面,采用传统欧拉-拉格朗日方法在时间步长为0.000 1 s时对于单液滴粒子的计算需要2 000 s,而采用欧拉网格逼近算法计算时间为10 s,能将计算效率提升1～2个数量级。

3.2 模型应用

根据文献[12]计算得到单位面积、单位时间产生的气泡数目。下部重力空间的直径为5 m,本文采用的网格尺寸为0.5 mm。根据二维截面得到气泡产生总频率f=250 s-1,即每隔20 mm的节点产生气泡频率f=1 s-1。假设每个气泡能产生8个喷射液滴,即可得到每个喷射液滴产生的频率。

图7～9为不同尺寸喷射液滴局部图,可看出,小于临界尺寸的液滴会随蒸汽流场向上流动,而大于临界尺寸的液滴会先向上做减速运动到最高点,然后下落回自由表面。

图7 半径1 250 μm的喷射液滴局部图Fig.7 Local view of jet droplet with radius of 1 250 μm

图8 半径1 356 μm的喷射液滴局部图Fig.8 Local view of jet droplet with radius of 1 356 μm

图9 半径1 500 μm的喷射液滴局部图Fig.9 Local view of jet droplet with radius of 1 500 μm

计算得到的重力空间中不同尺寸喷射液滴速度场和密度场如图10、11所示。

液滴半径:a——1 250 μm;b——1 356 μm;c——1 500 μm 图10 重力空间中的不同尺寸喷射液滴速度场Fig.10 Velocity field of jet droplet with different sizes in gravity space

计算得到的重力分离空间中的不同尺寸喷射液滴合并密度场如图12、13所示,密度较大的位置是可能的碰撞点位置,实际重力分离空间的液滴总密度图需要再叠加膜液滴的密度分布。

图12 不同尺寸喷射液滴合并密度局部图Fig.12 Local plot of combined density for jet droplet with different sizes

图13 不同尺寸喷射液滴合并密度图Fig.13 Combined density diagram for jet droplet with different sizes

4 结论

本文基于欧拉网格逼近理论基础提出了一套用于大液滴计算的模型,采用欧拉-拉格朗日方法进行了模型验证,两种方法的速度相对误差达到0.1%,密度相对误差达到0.6%。同时欧拉网格逼近方法的计算效率相较于欧拉-拉格朗日方法能提升1～2个数量级。本文将该模型应用于AP1000汽水分离器的下部重力分离空间中喷射液滴的模拟,分别得到达到悬浮的临界尺寸喷射液滴、小于临界尺寸和大于临界尺寸喷射液滴的运动图,且给出了合成密度图,合成密度图中密度较大的位置处是可能的碰撞点区域。本工作可为汽水分离等装置中的大液滴分析提供有效的计算分析工具。