APP下载

广义守恒相场简化多相流格子Boltzmann 方法

2022-07-13章诗婷肖鸿威周锦翔牛小东

空气动力学学报 2022年3期
关键词:液滴三相透镜

章诗婷,肖鸿威,周锦翔,牛小东,*

(1. 汕头大学 工学院,汕头 515063;2. 汕头大学 智能制造技术教育部重点实验室,汕头 515063)

0 引 言

多相流动现象不仅与自然界、人们的日常生活息息相关,也普遍存在于工业生产领域[1-2]。对于计算流体力学(computational fluid dynamic,CFD)研究人员来说,由于多相流固有的非线性、物理性质(如密度和黏度)的巨大差异,以及因质量、动量和能量同时进行交换而产生移动界面的复杂拓扑变化[3],因此开发一个简单有效的数值模型来研究这些复杂现象在实践中面临着巨大的挑战。到目前为止,已经有许多研究学者对多相流进行了数值仿真,并针对多相流问题发展了一些数值模拟方法,包括相场方法(phase-field method, PFM)[4-6]、 流 体 体 积 法( volume of fluid,VOF)[7-8]、水平集方法(level set method,LSM)[9-10]、前沿跟踪法(front-tracking method,FTM)[11-12]、光滑粒子流 体 动 力 学 方 法( smooth particle hydrodynamics method,SPHM)[13]。其中,相场方法是研究这类流动最常用的方法之一。在相场方法中,尖锐界面被视为有限厚度的过渡区,物理参数在其上平滑变化。

相场方法作为扩散界面方法中研究复杂界面动力学的一种,已被许多学者应用在模拟多相流问题上。相场方法主要分为两种:一种是用包含空间四阶导数的Cahn-Hilliard(C-H)方程[14]来捕捉界面变化,另一种是具有二阶精度的Allen-Cahn(A-C)界面控制方程[15]。当前,大多数描述多相流的相场方法都是基于C-H 理论,其方程源自梯度流法[16]。Lee 等[17-19]开发了一种基于总自由能假设的非整体相场方法;Kim[20]进一步改进了该方法,通过引入新的离散delta 函数提出了广义连续表面张力公式。与非整体相场方法不同,Elliort 和Luckhaus[21]提出了一种具有修正自由能公式的整体相场方法,对其引入了对称正定矩阵,其中的系数与不同流体之间的成对表面张力有关。Dong[2,22-24]使用整体相场方法进行了一系列应用。然而,我们观察到,由于总自由能随时间减少,C-H 方程在保持各相体积方面的性能较差[25]。此外,当构造非线性势时,在两两表面张力不均匀的情况下很难满足还原一致性条件。

格 子Boltzmann 方 法(lattice Boltzmann method,LBM)作为一种介观方法,自20 世纪80 年代末以来,已发展成为一种流行的CFD 替代工具。因为它的简单性、在并行计算机上的可扩展性以及易于处理复杂几何图形,被广泛应用于模拟流体动力学[26-27]和非线性方程组[28-29]。LBM 使用基于粒子的密度分布函数来演化宏观系统。最近,基于相场理论,许多研究者结合LBM 来仿真研究多相流问题。Liang 等[30]提出一个LB 模型,通过巧妙地构建流场LB 方程中的力项来进一步发展C-H 动力学。Reza 等[31]开发了另一个PF-LB 模型,使用C-H 方程跟踪三种不同流体之间的界面。然而,上述两种方法仍受到体积守恒和质量守恒的影响。在A-C 方程的动力学建模中,Reza 等[32]提出了一种描述三相流的守恒相场LB 模型,该模型可保持良好的质量守恒性,并且满足还原一致性条件。Zheng 等[33]针对不混溶多组分流体开发了另一种满足还原一致性条件的守恒PF-LB 方法。Hu 等[34]提出一种广义守恒PF-LB 模型,其中每一相相变量满足广义eikonal 方程和体积分数约束。

在本文中,我们在单松弛LB 框架下提出一种简化的多相流LB 方法,用于仿真具有大密度比、大黏度比的多相流问题。该模型基于Hu 等[34]提出的广义守恒PF 理论,其优点在于可以处理不可压、不混溶多相流,并具有良好的还原一致性以及质量和体积守恒性。此外,该方法保持了我们早期开发的用于两相流大密度比的SMLBM 的良好数值稳定性[35]。与用于模拟多相流的常用单松弛和多松弛LB 方法不同,简化多相流格子Boltzmann 方法(simplified multiphase lattice Boltzmann method,SMLBM)在单松弛LB 方法的框架内通过预测-校正策略来模拟流体系统和跟踪界面演化,其计算过程中仅需要考虑平衡分布函数的演化。因此,SMLBM 大大减少了计算中对虚拟内存的需求,简化了物理边界条件的实施,而这两点往往是限制传统单松弛和多松弛LB 方法应用的主要因素。最近,Li 等[36]构造了一种统一的SMLBM 方法来模拟大密度比、大黏度比的磁流体两相流问题,验证了此方法模拟多相磁流体流动的效率性和准确性。Khan 等[37]运用SMLBM模拟研究了均匀磁场下磁流体液滴变形和气泡在磁流体中合并的动力学过程。随后,Chen 等[38]发展了一种三相流SMLBM 方法,并将其应用于中等雷诺数剪切流中的复合液滴数值研究,但研究中的各相流体之间的密度比和黏度比依旧较小。因此,在模拟多相流系统中,解决大密度比、大黏度比多相流问题还较为困难,尤其是三相流及以上的多相流问题。但在现实生活中,气-液多相流中的各相流体之间的密度比普遍都达到了1 000 以上。所以为了提高模型对大密度比多相流问题数值仿真的鲁棒性,并克服单松弛模型对黏度的敏感性,需要提高相场方程界面捕捉的稳定性和精确性。因此,我们将广义守恒PF 方程中扩散项构造的源分布函数吸收到相场平衡分布函数中,这使得该方法的实现更加简单。

通过拉普拉斯定律和具有不同表面张力比的液滴透镜扩展两个基准问题的数值仿真,我们首先验证了本方法的稳定性和准确性。进一步,通过模拟不同黏度比的三相泊肃叶流问题以及不同密度比复合液滴在平板上的扩散算例,将数值结果与解析解进行比较,展示了本模型的精确性和处理大密度比、大黏度比复杂多相流问题的能力。

1 数值方法

1.1 广义守恒相场模型

相场模型实质上就是一类特殊的扩散界面模型。为了准确地描述整个流体系统中多种非混溶流体共存的状态以及他们在相互作用下的运动,在相场模型中将任意位置上多种非混溶流体所占的体积分数定义为该点的相变量φi,显然存在如下约束:

1.2 广义守恒相场简化格子Boltzmann 模型

显然,上述SMLBM 的整个实施过程只涉及平衡分布函数,而平衡分布函数又可以直接从宏观变量计算得到。与传统的单松弛和多松弛LB 方法相比,不再需要计算流场的碰撞演化过程,提高了模型的计算效率。此外,由广义守恒相场方程(4)中扩散项构造的源分布函数被吸收到式(20)中的相场平衡分布函数中,这使得该方法更加简单和稳定。

2 模拟算例

2.1 拉普拉斯定律

本节首先进行拉普拉斯定律的基准试验,以此来验证模型的准确性。我们将两个不同密度的液滴(流体1、流体2)同时浸入流体3 中,如图1 所示。在表面张力的影响下,液滴内部的压力大于液滴外部的压力,在达到平衡时,液滴会以圆形保持静止。同时液滴内外界面压力差Δpij与表面张力σij之间满足以下公式:

图1 两个液滴浸没到另一种流体中示意图Fig. 1 Diagram of two droplets immersed in another fluid

首先为了验证网格的收敛率,分别在计算域为x×y= 100×50、200×100、400×200 和600×300 时进行仿真。我们设定参考长度Lc= 2Rij,并将其分别设置为30、60、120 和180。根据下述公式计算了三种相变量在不同网格大小下的数值误差:

三种相变量的数值误差如图2 所示。可以发现三种相变量的误差呈现相同的趋势,即网格越小,相变量的误差越小,说明模型的收敛精度与网格的大小有关。当计算域为100×50 时,网格最大,其收敛精度为1×10-3;当计算域为600×300 时,网格最小,收敛精度达到1×10-4,均在二阶精度以内。从图3 中可以看出,400×200 和600×300 两种网格的误差非常接近,可以判断在这个网格数量区间内,模型的精度不会随着网格数量的增加而有较大变化,进一步证明了网格独立性。为了保证计算精度,同时尽可能地减少计算时间,我们选择计算域为400×200 的网格大小。

图2 当前模型的收敛误差Fig. 2 Convergence errors of the present model

图3 在不同网格大小下 E (φ3)随时间变化Fig. 3 The time histories of E (φ3) under different grid sizes

图4 中列出了界面压力差在两个液滴不同的表面张力(σ13、σ23)下和不同半径(Rij= 20、30、40、50、60)下模拟值与理论值的对比。图中可以清楚地看出模拟值和理论值吻合得非常好。我们对比了两个液滴在不同半径下相变量φ在水平中央线的模拟值和初始值,如图5 所示。从图中可以看出,在液滴半径为Ri j=20、60时,相变量的模拟值和理论值依旧重合得很好。显然,模型能够准确地追踪界面的变化且能够保持较高的稳定性。

图4 液滴1 和液滴2 界面压差的模拟值和理论值对比Fig. 4 Comparison of the pressure difference across the interfaces of droplet 1 and droplet 2 between the numerical results and theoretical solutions

图5 相变量在水平中央界面分布对比图(R ij=20、60)Fig. 5 Comparison of the simulated profile of the order parameter with the theoretical solutions (R ij=20, 60)

2.2 液滴透镜

本节研究了一个透镜在两层不相溶流体间的扩散运动这一经典算例,液滴透镜扩散问题在文献[10,32,41,42]中被广泛用于验证模拟多相流的数值方法。液滴透镜如果一开始存在于两相不混溶的流体之间,在表面张力的作用下透镜会发生形变,其最终的平衡状态是上下两层流体间相互作用的结果。透镜的平衡状态由三个表面张力系数决定,根据诺伊曼定律可以得出以下关系[43]:

图6 液滴透镜参数示意图Fig. 6 Diagram of the spreading of a liquid lens

表1 液滴透镜算例中表面张力系数的设定Table 1 Surface tension coefficients used in the simulation of spreading of a liquid lens

4 种情况下的透镜平衡状态如图7 所示,显然,不同的表面张力系数比会导致透镜最后的平衡形状不同。当σ23的值增加时流体2 的作用力会增大;同理,当σ13的值减小时,流体1 的作用力也会相应减小。正因为上下作用力的改变形成压力差导致了透镜形变的不同。我们将平衡时透镜两个三相点之间的距离d,以及h1和h2的模拟值和理论值进行了对比,如表2 所示。从表中可知,模拟值和理论值吻合得较好,最大的相对误差值为5.6%。这种微小的误差是由于在界面厚度为零的极限情况下计算了理论值。

图7 不同表面张力系数下液滴透镜的平衡状态Fig. 7 The equilibrium shape of the liquid lens for different surface tension coefficients

图8 展示了四组不同表面张力系数下液滴透镜质量随时间的变化。在图中,以液滴透镜的初始质量(m0)为评判标准,若质量守恒,液滴透镜的初始质量和实时计算质量(m1)的比值为1。从图中可以看出,四种方案的液滴透镜在演变过程中的质量比曲线与准确质量比1 的直线基本重合,质量耗散误差都在1%以内,表明本方法能够保持良好的质量守恒性。

图8 四种方案液滴透镜的初始质量和计算质量比值随时间的变化Fig. 8 Variations of liquid lens mass ratio with time for four cases

2.3 三相泊肃叶流

泊肃叶流动问题主要用于验证黏度比对模型计算精度的影响。这里我们模拟了三相泊肃叶流动问题。假设各相流体之间的界面是光滑的,对于相场的相变量初始化使用了如下函数:

其中,h为矩形槽道的高度,槽道长度为2h,水平槽道内充满了三层互不混溶的流体,如图9 所示。

图9 三相泊肃叶流示意图Fig. 9 Diagram of Poiseuille flow with three phases

在三相泊肃叶流动中,流体由一个体力G=ρg驱动,g是x方向的恒定加速度。因为各相的流体密度为常量,整个系统的N-S 方程表述如下:

表3 三相泊肃叶流算例中密度和运动黏度的设定Table 3 Densities and kinematic viscosities used in three-phase Poiseuille flow example

我们绘制了速度分量ux在y轴上的分布曲线,并将当前结果与Hu 等的结果进行对比,如图10 所示,发现三种情况下两者的曲线吻合较好。图10(a)中,由于三种流体之间的黏度比为1,沿y轴的速度剖面较为光滑。而在方案2 和方案3 中,黏度比增大,图10(b、c)非常好地捕捉到了因黏度变化而引起的相界面处明显的速度分量ux的跳跃。

图10 当前结果与Hu 等结果[34]之间速度剖面的比较Fig. 10 Comparison of velocity profiles between present results and Hu et al.[34]

2.4 复合液滴铺展

为了验证模型是否能够克服大密度比多相流交互界面处的较大压力梯度问题,这里采用一个二维复合液滴在光滑平板上的铺展问题来进行数值模拟。在这个算例中充分考虑了湿润边界条件,图11 中流体1、流体2 和流体3 为非混溶的三种流体,与固体壁面从左至右形成的接触角分别为θ13、θ12和θ23,由几何关系可知,θ21= π -θ12。为了满足三相表面张力系数平衡条件,流体i和流体j与光滑平面形成的接触角θij可根据Young’s 法则定义如下:

图11 三种流体与固体壁面接触角示意图Fig. 11 Diagram of a compound droplet on a solid substrate

为了解决三相流系统与固体平面接触面上接触角的不连续性问题,在湿润边界处理上设置加权接触角θ1和θ2来判定湿润度,保证界面均匀光滑地过渡。加权接触角的定义如下:

湿润边界的实施还需要人为在固体壁面下侧添加一层虚拟网格,如图12 所示。使用中心有限差分格式离散单元格(i,1)后,通过相邻网格的数值近似获得虚拟网格(i,0)的值[45]:

图12 湿润边界条件实施示意图Fig. 12 Diagram of the implementation of wetting boundary conditions

复合液滴的初始形态如图13 所示,流体1 和流体2 分别为一个四分之一圆,它们共同组成一个半圆复合液滴,周围区域为流体3。此算例中,三相流体的运动黏度 υ相同。为了验证模型在不同密度比下的准确性和稳定性,设置了四组不同密度比的复合液滴铺展算例进行对比,如表4 所示。复合液滴初始直径2R选作特征长度,计算域设定为4R×2R,网格划分成400×200 的 均 匀 网 格,界 面 厚 度 ξ = 4,参 考 速 度Uc= 0.01,流体区域四周边界均采用无滑移边界条件。

表4 复合液滴铺展中三种流体的密度Table 4 Densities used in the simulation of spreading of the compound droplet

图13 复合液滴在平面上铺展初始形态Fig. 13 Initial shape of the compound droplet

在该算例中将湿润边界条件接触角设置为:θ13=90°、θ12= 120°、θ23= 60°。假设任意两相界面间的表面张力系数相同,其值为0.01,则有三相角φ1=φ2=φ3= 120°。整个复合液滴不考虑重力因素,主要在表面张力和黏性力作用下推动液滴形态变化。这里用无量纲参数Re和We来表征液滴的动力学特性,设置Re=30,Weij=2。

图14 表示了在方案1 情况下,复合液滴在光滑平板上随时间发展的状态图,不同颜色的线分别表示不同时间节点复合液滴所处状态,这里时间节点分别取t= 0、0.10、0.25、0.50、0.75、2.00。从图中可看到复合液滴随着时间的演变,慢慢地向外铺展,当复合液滴到达平衡状态时其形态会保持不变,平衡时复合液滴的形状主要由初始设置的接触角和三相角给定。这里,我们令流体1 和流体2 与固体平面接触的润湿长度分别为L1和L2,其总润湿长度为L,如图15 所示。其中润湿长度的理论值是根据Zhang 等[46]计算的结果,并与模拟值进行了对比,如表5 所示。从表中可以发现误差值均小于1%,理论值与数值模拟结果非常吻合,验证了该模型的精确性。

图14 方案1 的复合液滴随时间演化过程示意图Fig. 14 Time evolution of the compound droplet of Case 1

图15 方案1 的复合液滴平衡状态Fig. 15 Equilibrium state of a compound droplet of Case 1

表5 方案1 的润湿长度理论值与模拟值对比Table 5 Comparison of the numerical results and analytical solutions of wetting length in Case 1

方案2~方案4(如图16 所示)的复合液滴的平衡形状与方案1 相似。但是,随着密度比的增加,复合液滴达到平衡状态所需的时间逐渐增加。这是因为密度比的增加导致压力梯度和界面上黏性力的增加,从而进一步延迟复合液滴扩散到平衡状态。然而,即使对于密度比高达1 200 的方案4,本模型仍能给出稳定解。四种情况下接触角的数值结果与表6 中的理论值进行了比较,其中 Δθij表示各接触角的相对误差。结果表明,误差均小于2%,进一步证明了模型的准确性。

表6 不同密度比复合液滴平衡状态下接触角的模拟值与理论值对比Table 6 Comparison of the numerical results and analytical solutions of the compound droplet for different cases under equilibrium state

图16 不同情况下复合液滴的平衡形态Fig. 16 The equilibrium shape of the compound droplet for different cases

为了更清楚地观察不同密度比下复合液滴的平衡形态,将四种情况下复合液滴平衡态图绘制在一张图中,如图17 所示。我们发现了一个有趣的现象,随着密度比的增加,流体1 和流体2 之间的交界线逐渐向左移动。此外,在较大的密度比下,润湿边界总长度L也增加,三相边界点不仅向左移动,而且随着密度比的增加向上移动。显然,在方案4 中,复合液滴中流体1 和流体2 的形状不像方案1 那样规则,从而增加了总润湿长度。有趣的是,四种情况下设置的接触角是相同的,并且由于密度比的增加,复合液滴的形状并不规则,无法计算润湿长度的理论值,但是最终平衡状态下接触角仍与预设的平衡接触角一致。

图17 不同密度比下复合液滴平衡状态图Fig. 17 Equilibrium shape of compound droplets at different density ratios

为了更加直观地观察三相交界点的变化规律,绘制了在总润湿长度L的变化过程中,不同密度比下三相交界点高度H的变化,如图18 所示。首先在图中可以明显地看出,密度比增大的同时,复合液滴三相交界点的高度随之增高,并在复合液滴演化的过程中一直保持这种趋势。其次,在复合液滴到达平衡状态时,密度比越大,总湿润长度也会越长,其中方案4 的L是增长最显著的。这一现象可归因于两种流体之间的作用力和压力差,它们随着密度比的增加而增加,因此复合液滴扩散得更多。

图18 不同密度比下复合液滴三相交界点高度随总润湿长度的变化曲线Fig. 18 Curves of H varying with L for different cases

3 结 论

本文提出了一种广义守恒相场简化多相流格子Boltzmann 方法。该方法继承了Hu 等推导出的广义守恒相场方程的还原一致性和质量体积守恒的优点,并结合早前提出的SMLBM 用于处理大密度比、大黏度比问题,其中广义守恒相场方程的源分布函数被吸收入相场平衡分布函数中,以保持方法的良好稳定性。此外,本方法大大降低了内存的需求,这通常是限制传统单松弛或多松弛LB 方法应用的一个主要因素。为了验证本方法在处理多相流时的网格独立性和精确性,我们首先模拟了拉普拉斯定律和液体透镜两个基准算例,并获得了令人满意的结果。其次,对不同黏度比下的三相泊肃叶流进行了仿真,验证了此方法能够模拟大黏度比(最高可达500)问题。最后,对不同密度比下的复合液滴在润湿基底上的扩散进行了模拟,数值结果表明,该方法能够模拟大密度比多相流问题,其中最大密度比可达1 200。对此还发现在较大的密度比下,复合液滴轮廓变得不规则,导致总润湿长度增加,并且还找到了复合液滴在不同密度比下三相交界点的变化规律。以此我们可以得出结论,本文所提出的方法在模拟多相流的准静态问题上有较好的表现。

猜你喜欢

液滴三相透镜
涡扇发动机内外涵道水量分布比例数值仿真
透镜(部级优课)
射流过程中主液滴和伴随液滴的形成与消除研究
钠火分析程序-NACOM介绍
三相异步电动机在维修电工中的难点研究
未来的眼镜或许比纸还薄
油田三相分离器效能分析
油田三相分离器效能分析
采用并联通用内模的三相APF重复控制策略
透镜及其应用专题复习