基于LevelSet方法的气_液_固三相流动模型与模拟_李彦鹏

第25卷　第4期应用力学学报Vol.25　No.4 2008年12月CHINESE JOURNAL OF APPL IE D MECHANICS Dec.2008文章编号:100024939(2008)0420578205基于Level Set方法的气2液2固三相流动模型与模拟3李彦鹏1　王焕然2(长安大学　710064　西安)1　(西安交通大学　710049　西安)2摘要:采用基于Level Set方法与离散颗粒模型相结合的方法,建立了一个用于描述气液固三相流动的新模型。

模型耦合了颗粒与气泡、颗粒与液相以及气泡与液相之间的相互作用。

应用该模型对液固悬浮液中的典型现象———气泡的单孔及多孔形成过程以及颗粒夹带现象进行了三维模拟,检验了其可行性。

并进一步研究了颗粒的存在对气泡的形成与上升过程的影响以及气泡诱导的液相流动对颗粒行为的影响。

研究结果表明,所提出的模型能够真实地预测三相流中气泡与颗粒分散相的特征,为研究多尺度的三相流动提供了一种新途径。

关键词:气液固三相流;Level Set方法;离散颗粒模型;数值模拟中图分类号:O35911 文献标识码:　A1　引　言气2液2固三相流动过程广泛存在于石油、化工、冶金、能源及环境等工程实际中,比如石油开采中的油2气2砂流动,化工反应装置中的气液固流态化和能源工程中的高钙煤灰颗粒的增湿脱硫过程。

三相流的流体动力学特性是由每一相的运动以及各相之间的相互作用共同决定的。

因此,三相流动问题远比单相或两相流动复杂,对其的模化及模拟已成为流体力学领域极具挑战性的问题之一。

目前国内外研究者对气2液2固三相流的模拟主要采用欧拉2欧拉(Eulerian2Eulerian)方法[122]和欧拉2拉格朗日(Eulerian2Lagrangian)方法[324]两种。

然而,这两种模型只能描述三相流动的宏观特性(诸如时间平均的流型及气含率的空间分布等),不能反映诸如气泡的变形与运动这样的微观流动状态。

而在气泡及颗粒为分散相的三相流动系统中,气泡和固体粒子的行为对深入理解气泡、颗粒与流体的相互作用,进而对理解整个三相流动现象都具有重要意义。

直至最近,随着计算机技术的发展,一些研究者开始尝试利用界面追踪技术对三相流中气泡的行为进行直接模拟,比如Li等人[5]利用CFD2VOF2 DPM模型对二维气2液2固三相流化床内气泡尾涡的运动进行了数值模拟。

Annaland等人[6]采用Front t racking方法结合离散颗粒模型模拟了气2液2固三相流中气泡的上升行为。

由于VO F(Vol2 ume of Fluid)方法和Front t racking方法在气液交界面重构上的复杂性,限制了其在三维体系下的应用。

相比较而言,Level Set方法在模拟三维相界面发生拓扑结构变化的过程具有简单方便的优势。

因此,本文以Level Set方法为基础,结合离散颗粒模型,建立一个新的三相流模型,并对液固悬浮液中气泡形成和上升过程进行数值模拟,分析其气泡与颗粒行为,揭示分散相运动的细节。

2　三相流模型211　气液相界面运动的Level Set方法Level Set方法是由Osher和Fedkiw[7]提出的3基金项目:国家自然科学基金(20706007);教育部科学技术研究重点项目(108116)　来稿日期:2007201216　修回日期:2007207216第一作者简介:李彦鹏,男,1973年生,长安大学环境科学与工程学院,副教授;研究方向———多相流。

E2m ail:liyanp01@一种隐式界面追踪方法,已经在许多复杂的界面追踪问题中应用成功。

在此方法中,假设在计算区域中气、液两相占据的区域分别为Ω1(t )、Ω2(t ),两相分界面为Γ(t )。

构造Level Set 函数<(x ,t ),使得在任意时刻t ,Γ(t )恰是<(x ,t )的零等值面,即要求Γ(t )={x |<(x ,t )=0}(1)函数<(x ,t )的初值应满足在Γ(t )附近为法向单调,在Γ(t )上为零。

一般可取<(x ,0)为x 点到界面Γ(0)的符号距离<(x ,0)=d (x ,Γ(0)), x ∈Ω10, x ∈Γ(0)-d (x ,Γ(0)),　x ∈Ω2(2) 随着时间的推进,Level Set 函数<(x ,t )的流体动力学行为由下述方程控制5<5t+V ・ <=0(3)式中V 为流场速度矢量。

在任意时刻,只要求出<(x ,t ),就可以通过取该函数在三维空间的零等值面来隐式确定气液相界面的位置。

关于Level Set 方法的具体细节可参见文献[8]。

212　气液两相流动的控制方程在定义了Level Set 函数后,由气液相界面分隔的不相混的气液两相就可以处理为单一流体体系。

这样,考虑颗粒相存在的影响,气液两相的控制方程可由单相流的质量与动量守恒方程采用体积平均技术得到5εf ρ5t+ ・(εf ρV )=0(4)5εf ρV 5t+ ・(εf ρVV )=-εf p +εf ・τ+εf ρg +f sf +f bf(5)式中:ρ、V 、p 分别表示流体的密度、速度矢量、压力;τ为粘性应力张量;εf 表示流体的体积分数,它满足εf +εs =1,其中εs 为固相的体积分数;f sf 和f bf 分别表示单位体积的固体和气泡(气液相界面)对流体的作用力,其表达式将在后面给出。

需要注意的是,在上述两相流单一场控制方程中的流体密度与粘性在空间上是变化的,其定义为ρ(<)=ρg +(ρ1-ρg )H (<)(6)μ(<)=μg +(μ1-μg )H (<)(7)式中:下标l 和g 分别表示液相和气相,H (<)是一个Heaviside 函数。

213　颗粒运动的离散颗粒模型对于三相流中固体粒子的运动,通过在Lagrangian 坐标系中分析每一个颗粒的运动而进行描述。

本文将颗粒的运动分解为受颗粒作用的碰撞过程和受流体作用的悬浮过程两部分。

1)　在颗粒的悬浮过程,每个颗粒的运动遵循Newton 第二定律。

由于本文所考虑的粒子直径较小,远小于计算网格的尺寸,忽略粒子在这个过程中的旋转运动,则其运动方程为d x p d t =V p (8)m p d V p d t =m p g -π6d 3p p +F f s +F f b (9)式中:x p 、V p 分别表示颗粒位置、速度矢量;d p 为颗粒直径。

方程(9)右边各项分别表示重力、颗粒位置处的流体压力梯度力、流体2颗粒作用力和气泡2颗粒作用力。

由于考虑的颗粒尺寸较小,本模型中流体-颗粒作用力F f s 只包括了曳力和附加质量力,其中颗粒所受的曳力F d 采用Wen 等人[9]的计算公式F d =18πd 2p C ′d ε2f ρf |V -V p |(V -V 0)(10)这里c ′d 表示有效曳力系数。

对于球形颗粒,附加质量力F am 的表达式为F am =112πd 3p ρ(d V d t -d V p d t )(11) 当颗粒接近气泡时,表面张力会通过液膜作用于粒子。

由于计算网格远大于液膜厚度,采用气泡诱导力模型(BIF )[5]计算气泡-颗粒作用力F f b F f b =π6d 3p σK (<)δ(<) <(12)式中:σ、κ(<)、δ(<)、<分别表示表面张力系数、界面曲率、Dirac Delta 函数、Level Set 函数。

2)　对于颗粒间的碰撞运动,采用Hoomans 等人[10]的硬球模型处理。

该模型假设碰撞是二元双向瞬时碰撞,接触点为一点,颗粒为球形的准刚体。

因此可根据碰撞动力学中的动量守恒定律确定碰撞后颗粒的速度变化规律。

214　相间耦合在控制方程(5)中气泡对流体的作用力f bf 实际上是由气液相界面的表面张力所诱导。

本文采用连续表面力模型(CSF )[11]将此表面力转化为界面附近的体积力,即有f bf =σκ(<)δ(<) <(13) 根据Newton 第三定律,作用在颗粒上的流体力将会施加给流体反作用力。

因此,方程(5)中的流体2颗粒作用力项f sf 可由下式给出f sf =-ΣF n f s ΔV n ij k,　(x n p ∈Ωij k )(14)975第4期 李彦鹏,等:基于Level Set 方法的气2液2固三相流动模型与模拟式中:下标ij k 定义了计算网格的位置;Ω和ΔV 表示网格的区域与体积;x n p 是颗粒n 的位置向量。

215　数值方法基于上述模型,我们自编了三相流数值模拟代码。

其中,气液主流场控制方程采用三维中心网格的有限容积方法进行离散,采用AL E (Arbit rary Lagrangian Eulerian )算法计算流场。

离散中,采用二阶通量受限的Barton 格式计算对流相,采用标准的二阶中心格式离散扩散相。

压力泊松方程采用鲁棒性好的ICC G 算法进行求解。

对于Level Set 方程的求解,采用三阶ENO 迎风格式计算对流项V ・ <,采用三阶TVD Runge 2Kutta 格式推进时间步。

颗粒运动方程采用简单的一阶显式积分格式求解。

计算中所用到的欧拉体系下的流体相参数,比如流体相体积分数εf ,压力与速度均采用体积平均的方法进行插值计算。

对于颗粒与颗粒以及颗粒与固壁之间的碰撞,采用Hoomans 等人提出的事件驱动算法来判断颗粒间相互作用的先后顺序。

3　结果与讨论为了检验所建立的三相流模型的正确性,对研究比较成熟的单个气泡在液固悬浮液中的上升速度进行数值计算,并与Luo 等人的实验数据[12]进行比较,其结果如图1所示。

数值模拟采用与实验相同的参数:气泡直径10mm ,从静止状态开始上升,空气的密度和黏度分别为111kg ・m -3和118×10-5Pa ・s ;甘油水溶液的密度和黏度分别为1206kg ・m -3和010529Pa ・s ,表面张力系数是010629N ・m -1;1000个粒径为1mm 、密度为215kg ・m -3的玻璃球随机布置在同上的计算区域。

从图1中可以看到,模拟结果与实验结果吻合较好,二者的相对误差最大为6%,这表明本文的数值模型是合理的。

图1　液固悬浮液中单气泡上升速度的比较311　单孔气泡在液固悬浮液中的生成下面,对气泡通过单孔在液固悬浮液中的生长现象进行模拟。

计算模型为一尺寸为4cm ×4cm ×8cm的方截面柱体,容器中的液体介质为Paratherm NF 传热流体,底部正中有一直径为014cm 的曝气孔,以较低的气体流量供气。