电场作用下气泡上升行为特性的数值计算研究1)

王悦柔 王军锋 刘海龙

（江苏大学能源与动力工程学院，江苏镇江 212013）

引言

气液多相流不仅存在于自然现象中，在工业生产中也有广泛的应用[1-7]，如核池沸腾、液体分离、废水处理和化学反应等.对流体中上升气泡的运动特性及规律的认识对相关工业装置的设计和运行十分重要.为此研究人员已经对黏性流体中的气泡动力学[8-20]进行了的大量研究.Clift 等[8]最初对气泡上升过程的研究进行了回顾，表明气泡的运动可以用无量纲数描述，即雷诺数Re、莫顿数Mo和欧特沃斯数Eo，后者也可以称为邦德数Bo.并且他认为在低Re和Bo下气泡能维持球状，随着Re和Bo的增加会产生不同的气泡形态，如椭圆帽状、球帽状、椭球状以及摆动形状.之后的研究表明[9-11]，在Re和Bo足够大的情况下，气泡可能会变为环形.Tripathi 等[12]在Mo较低的情况下发现了一种新的气泡破裂模式，气泡变为球状并产生卫星气泡.王焕然等[13]发现在Mo较低，Re较高的情况下，气泡形状表现出周期振荡变化.此外，其他参数也可能影响气泡上升机制.Chen 等[14]研究了雷诺数、邦德数、密度和黏度比对上升气泡的影响.结果表明，Re和Bo以及密度比的增大会增强气泡下方液体射流的影响，从而导致气泡中心穿孔.Amaya-Bower[15]分析表明，黏度比对气泡的终端速度和形状的影响更大，随着黏度比的降低，剪切应力增大，气泡变形程度和终端速度减小.Alizadeh 等[16]在二维和三维数值研究的基础上验证了表面张力对气泡形变行为的影响符合拉普拉斯定律.张洋等[17]通过数值研究表明当域约束比Cr8 时可忽略壁面对裙带气泡特性的影响.

目前，多数研究集中在对非电场作用下的气泡动力学特性上.将电场应用于不同介质的多相体系中，能够有效地控制气泡的运动及变形[21-27]，使气泡获得较好的流动状态和上升特性.早期，Taylor[21]和Garton 等[22]对这一问题进行了经典研究，受电场作用的液滴和气泡会沿电场方向发生变形，在电场方向上呈现出一个扩展的椭球体形状.陈凤等[23-24]在蠕动流近似的基础上分析计算了具有不同介电常数工质的气泡在外电场力作用下的受力状况以及气泡内外的速度场分布.计算结果表明电场力的水平分力使气泡沿场强方向伸长，使得气泡内部的流体运动加剧.并且液体和气体相对介电常数之比越大的工质所受的电场力越大.Yang 等[25]的研究表明电场形式对气泡运动也有影响，垂直电场加速了气泡的上升运动，而水平电场阻碍了上升运动.Rahmat 等[26]研究了电毛细管数与域限制比对气泡（特别是环形气泡）上升的影响，结果表明，当电力足够强时，气泡就会被穿透.对于较小的域限制比率，需要更大的电毛细管数来穿透气泡.Andalib 等[27]通过实验研究了不同电场形式下单个粗气泡形态及运动路径的变化过程.然而实验条件的限制使得多数情况下的气泡行为特性难以通过实验捕捉，目前的多数研究也是依靠数值计算手段实现.且由于静电学与流体动力相互作用的复杂性，相关物性参数对受电场力作用的气泡形态及上升速度的影响的研究尚不充分，对于电场作用下气泡的运动特性及电场的作用机制仍需开展深入系统的研究.

目前已有多种数值方法被开发应用于气液两相流体动力学研究中，包括边界积分法、前跟踪法、流体体积（volume of flui，VOF）法、水平集（level set，LS）法和光滑粒子流体力学（smoothed particle hydrodynamics，SPH）法等.边界积分法[28]和前跟踪法[11,29]是精确计算体积力的显式方法，但实现起来非常复杂，且需要消耗大量的计算时间.VOF 法、LS 法和相场法是一种隐式方法，可以对界面的破裂、聚并等拓扑变化进行常规处理.与相场法[25]相比，VOF 法和LS 法在电流体动力学中得到了广泛的应用.VOF 方法[14,26,30]很好地解释了界面的拓扑变化，满足质量守恒，但不能准确计算界面附近的参数.LS 方法[32-34]在保持良好质量守恒的前提下，能准确捕捉两相流体界面.另外，一种利用具有密度和黏度等物理性质的质点来描述流体的SPH方法也可被用于多相流数值研究中.这种方法非常适用于模拟大变形、移动边界和多相混合的动态现象.然而从已有文献可知，目前SPH 法[35]还没有在电流体动力学中得到广泛的验证与应用.

本文借助数值手段对外加电场作用下单个气泡在流体中上升的动力行为进行了模拟研究.利用水平集方法，通过求解电场方程与Navier-Stokes 方程，捕获气液两相界面.该方法能够精确计算界面法向和曲率，平滑界面周围的物理性质，准确地捕捉气液两相流的拓扑变化.之后从网格无关性及与前人研究对比两方面验证了该模型的有效性，通过雷诺数、邦德数和电邦德数等不同参数研究了液体黏度、表面张力和电场力对气泡运动及变形的影响，并给出了不同情况下的气泡形态图.

1 数学描述

本文中的气液两相流体被认为是具有不同介电常数的介质，在没有外加电场的情况下，在流体本体区或界面处没有自由电荷.用于本文的二维模型如图1 所示，初始半径为R的气泡位于高H=16R，宽W=8R的矩形计算域中，从气泡中心距离底边h1=3R处释放，气泡上升至距离顶边h1处停止计算.这里选定R=1 进行模拟，模拟中所用的单位为mm.计算域的四周边界上应用无滑移壁面边界条件.在竖直方向上施加电场强度E∞=Ut/H的均匀电场.电势在顶部边界上为Ut，在底部边界处为0.

图1 均匀电场中气泡上升的二维模型示意图Fig.1 Schematic diagram of a two-dimensional model of a rising bubble under a uniform electric fiel

1.1 水平集方法

水平集方法通过跟踪水平集函数φ 的等值线来确定流体界面.水平集φ=0.5 决定界面的位置.控制φ 的传递和重新初始化的方程为

式中，u为流体速度，γ 和ε 为重新初始化参数.参数ε 确定界面周围的层厚度，其中φ 从0 变为1.γ 参数决定重新初始化的数量.

研究分为相位初始化和时间相关性研究两个阶段.相位初始化部分求解初始接口到计算域节点Dwi的距离.时间相关步骤根据以下表达式使用水平集函数的初始条件

式中，正号用于界面内部初始域，负号用于界面外部初始域.

由于水平集函数是一个平滑阶跃函数，因此通过以下公式确定全局密度和动力黏度

式中，ρb，ρl，µb和µl分别表示气体和周围流体的恒定密度和黏度.

1.2 流体流动方程

对于电场作用下的流体流动，流体流动受惯性力、黏性力、电场力等控制.在不可混相不可压缩牛顿流体的假设下，将N-S 方程表示为

式中，ρ 表示密度，u是速度，g为重力加速度，t为Fst表示表面张力，Fes表示外加体积力，这里即代表电场力.

表面张力由下式定义

式中，σ 为表面张力系数，I为单位矩阵，n为界面单位法向，δ 为狄拉克δ 函数，该函数仅在流体界面为非零值.界面法向的计算公式如下

水平集参数φ 还用于通过由下式定义的平滑函数来近似δ 函数

1.3 电场方程

在电流体力学中，由于磁感应效应是可以忽略的，所以电场强度E被假定为无旋的 (∇×E=0).对于理想的介质材料，电场的控制方程可以写为

式中，ε0为真空的介电常数，εr为相对介电常数.在没有任何时变磁场的情况下，电场的卷曲为零.电场可以用电势梯度表示

因此，式(10)可以根据电势V重写

作用在两相流体中的外加体积力Fes即为电场力，由麦克斯韦应力张量的形式给出

式中，D为电位移场，由下式给出

1.4 无量纲参数定义

式中，R是气泡半径，E∞是未扰动电场强度.下标b 和l 分别表示气体和周围液体.

2 模拟结果

2.1 模型验证

为了研究结果的网格依赖性，采用不同的网格数模拟ρ∗=0.001，µ∗=0.01，ε∗=0.04，Re=20，Bo=0.3，Boe=1 条件下的上升气泡.图2 比较了不同网格数下T=5.74 时的气泡的形状及位置.这里用∆表示最大网格尺寸，可以发现，随着R/∆的增大，气泡形状及位置的差异逐渐减小，R/∆=70 之后，气泡形状的差异较小，误差在5%以内，考虑到计算精度和计算机资源成本，这里选用R/∆=70 的网格尺寸条件进行接下来的模拟.如果区域的宽度远大于气泡的半径，则可以忽略区域的边界效应.根据研究表明[36]，当域高和域宽分别小于气泡直径的8 倍和3 倍时，边界效应对气泡的影响较大，这里设定域高H=16R，宽W=8R.此外，当前2D模拟结果与文献[19]的实验结果及对应的3D 模拟结果进行了比较.气泡最终形态与终端雷诺数Re∗的对比结果分别如表1 和图3 所示.

图2 不同网格大小下气泡形态及位置Fig.2 Bubble geometry and position with different grid sizes

表1 当前模拟结果与前人实验结果的对比Table 1 Comparison between current numerical results and previous experimental results

图3 数值模拟预测的最终气泡形态与前人实验结果[19]的比较Fig.3 Comparison of the bubble terminal deformation predicted by numerical simulation with previous experimental results[19]

可以看出，2D 仿真预测的气泡最终形态与前人实验观测结果及3D 预测的气泡形态基本一致，且与气泡终端速度相关的终端雷诺数Re∗误差在5%以内，因此2D 模型足以用于当前对气泡形态变化及上升速度的数值研究，且比3D 模型节省了更多的计算量，验证了该模型的有效性.

2.2 非电场作用下的气泡上升行为

初始气泡形状假设为球形，在本研究中，这一假设将继续成立.接下来，讨论了流体黏度和表面张力对气泡动力学的孤立影响.这里采用ρ∗=0.001，µ∗=0.01 条件进行接下来的模拟，通过调整黏度和表面张力系数获取不同Re和Bo，预测的气泡在h/R=13 处的最终形状如图4 所示.

一般情况下，当Re<40，Bo<30 时气泡沿直线上升且最终形成沿直线对称的完整的球状或扁椭球状[12].随着Re和Bo的增大，气泡可能会发生破碎与穿透为环状，或者出现路径不稳定的上升的现象，在这种情况下，需要一个三维模型来更好地预测气泡的速度和形状，而不是本文所使用的二维模型.气泡形状主要受惯性力(即重力)、表面张力和黏性力的控制.表面张力和黏性力阻止了气泡的变形，并努力使气泡保持球形.然而，惯性力导致气泡变形.当惯性力大于黏性力和表面张力时，气泡发生变形.在初始阶段，气泡的变形和加速上升是由于表面张力，惯性力和黏性力的不平衡.在达到力平衡后，气泡形状波动微小并稳定上升.这里通过改变气泡周围流体的黏度及表面张力获得Re<40，Bo<30 范围内的不同Re和Bo.从图4 可以发现，在Re和Bo均较小的情况下，气泡更偏近球状，随着Re和Bo增大，气泡最终呈扁长状.即液体黏度和气液两相间的表面张力越大，气泡在上升过程中克服惯性力，使得气泡维持球形.随着黏性力及表面张力的减小，气泡不再能维持球状，受惯性力影响发生变形.

图4 根据不同Re 和Bo 预测的上升气泡的最终形状Fig.4 The terminal shape of the rising bubble predicted by different Re numbers and Bo numbers

2.3 电场强度对气泡形状及速度的影响

在电场作用下，介电液滴在介质流体中的两相界面上受到电场力的作用，会沿电场方向发生拉伸变形[37].定义长径比Ar=Wb/Hb来反映气泡的变形，Wb和Hb分别表示气泡的宽度和高度.通过改变Boe来研究电场强度对气泡运动和变形的影响.可以知道，当Ar>1 时气泡发生横向变形，Ar越大，形变程度越大；反之，则发生纵向变形，即沿电场方向拉伸变形，Ar越小，变形程度越大.图5 给出了Re=20，Bo=0.3 条件时不同Boe下气泡长径比Ar随T值的变化，气泡初始形状均设为球状.当Boe=0 时，随着T值的增大，Ar值逐渐增大，并趋于稳定，Ar>1，气泡沿水平方向变形，最终变形为扁椭球状.随着Boe增大，Ar减小，气泡被克服表面张力的电力沿电场力方向拉伸，电场越强，Ar越小，气泡形状拉伸程度越大，且Ar波动变化，较难达到稳定值.分别选取Ar达到幅值时的不同特征时刻T下的气泡形态图，如图6 所示.电场力和表面张力之间的相互作用主导了上升的早期阶段，气泡沿电场呈明显拉伸状态，Ar达到一个最小值，之后气泡形状出现竖直和水平方向交替的摇摆不定的振荡形式，且Boe越大，振荡幅度越大.随着气泡上升，由于黏滞效应，摇摆运动的振幅随时间而减小，当气泡上升时，惯性力效应进入并与其他机制竞争.此外可以发现，在施加了电场的情况下，随着气泡垂直上升，不同Boe下，当前计算工况下的气泡均发生了六次振荡，气泡Ar达到幅值处的特征时间也相差不大，气泡形状发生振荡的时间节点基本一致，不因电场强度的改变而发生过大改变.

图5 气泡长径比Ar 随T 的变化图：Re=20，Bo=0.3Fig.5 Temporal evolution of the aspect ratio Ar of the gas bubble：Re=20，Bo=0.3

图6 不同Boe下气泡形态演化过程：(a)Boe=1，(b)Boe=2，(c)Boe=3，(d)Boe=4Fig.6 Bubble shape evolution under different Boe：(a)Boe=1，(b)Boe=2，(c)Boe=3，(d)Boe=4

与此同时，不同Boe下气泡上升特征速度U值随T值变化的曲线如图7 所示.在没有电场的情况下，上升速度主要受两种力控制：密度差引起的浮力和黏性效应引起的阻力.然而，流动阻力与气泡的形状密切相关[11,38].其中，当Ar较小时，流动阻力较小，上升速度较快，反之亦然.根据本研究结果可以发现，这种情况也适用于加电情况.气泡的上升速度是通过气泡质心位置的时间差值得到的.由图7 可见，施加垂直方向的电场导致气泡速度增大.电场越强时气泡上升得更快，原因主要是由于由电场引起的垂直拉伸导致Ar减小，引起的气泡所经受的阻力减小，气泡速度增大.由于气泡的不稳定振荡变形带动气泡周围流体的流动，阻力发生变化，上升速度值也随着气泡的不稳定抖动变形出现振荡形式，从图中可见速度U发生振荡时间节点与纵横比Ar振荡点几乎一致.

图7 气泡上升速度U随T的变化图：Re=20，Bo=0.3Fig.7 Temporal evolution of the vertical velocity U of the gas bubble：Re=20，Bo=0.3

2.4 表面张力对气泡形状及速度的影响

保持其他参数值不变，改变表面张力，探究不同Bo下气泡受电场力作用的上升行为.Bo象征重力和表面张力的比值，较小的Bo对应着较大的表面张力，相应的浮力效应较小，不加电场的情况下，Bo越小，气泡形状越圆整.外加电场时，气泡长径比变化如图8 所示.对于所有的Bo值，气泡上升初始阶段响应电场力，沿电场方向呈明显拉伸状态，Ar均小于1.Bo为0.1∼0.5 时，表面张力效应较大，Ar一直小于1，气泡呈长椭球状，Bo越小即表面张力越大的情况下，随着气泡的上升，气泡的长径比Ar相对越小，即气泡沿电场方向的变形程度越大；此时，Bo的减小使得气泡Ar值达到最小值的时间节点提前，Bo越低时，气泡所受的各种作用力间的相互竞争效应更强，气泡形状振荡越频繁，不稳定性更高.可以看出，当前计算域内，Bo=0.1 时气泡Ar达到振幅的振荡次数几乎是Bo=0.3 情况下的两倍.此外，如前文所述，表面张力效应增大导致相应的Ar值变小，阻力减小，气泡速度U越大（如图9）,且均呈不稳定振荡状态，Ar的低峰值点和高峰值点分别对应相同时间节点的U的高峰值点和低峰值点.

图8 气泡长径比Ar 随T 的变化图：Re=20，Boe=2Fig.8 Temporal evolution of the aspect ratio Ar of the gas bubble：Re=20，Boe=2

图9 气泡速度U 随T 的变化图：Re=20，Boe=2Fig.9 Temporal evolution of the vertical velocity Uof the gas bubble：Re=20，Boe=2

Bo=1 和Bo=10 时，气泡最终长径比Ar>1，此时气泡惯性力效应更强，表面张力效应较低，气泡在上升初期响应电场力发生竖直方向拉伸后，受惯性力作用迅速拓展发生水平方向拉伸变形，并且Ar逐步趋于稳定值，没有再发生明显的振荡情况，可以看出，Bo=1 时气泡后期呈现碗状，Bo=10 时，气泡在T=4.504 处呈现为盖状.气泡速度U值逐渐增大后也达到相对稳定的值，不再发生明显振荡.并且，Bo越小即表面张力越大的情况下，气泡Ar值越小，气泡沿电场方向发生拉伸变形程度越大，速度U越小.

2.5 液体黏度对气泡形状及速度的影响

黏性力也是影响气泡运动和变形的重要因素.Re表征惯性力与黏性力的比值，这里通过改变流体黏度获取不同大小的Re，Re对气泡形状和速度的影响如图10 和图11 所示.较大Re状态下，Re在10∼30 范围内时，由于黏性力较小，气泡对电力响应得更快，并因此沿电场方向伸长更多，Ar达到一个最小值，气泡上升初期的变形程度最大，随后由于各作用力相互竞争效应，气泡发生不稳定振荡变形，Re越大，即黏性力越小的情况下，振荡程度越大，之后变形波动程度逐步减小，并达到相对稳定状态，气泡初期与后期变形程度也相差更大.在Re=30 时，电场使得气泡沿电场方向发生一定程度的拉伸，气泡在上升后期由于惯性力的主导作用，出现Ar>1 的状态，即气泡发生横向变形.

图10 气泡长径比Ar随T的变化图：Bo=0.3，Boe=2Fig.10 Temporal evolution of the aspect ratio Ar of the gas bubble：Bo=0.3，Boe=2

图11 气泡速度U随T的变化图：Bo=0.3，Boe=2Fig.11 Temporal evolution of the vertical velocity U of the gas bubble：Bo=0.3，Boe=2

而当黏性力较大时，如Re=1 和5 的情况下，气泡首先响应电场力发生拉伸变形后，黏滞力阻碍了惯性力作用，使得气泡维持稳定的长椭球状，不再发生振荡.随着Re的增大，气泡形态发生不稳定振荡，然而Ar达到幅值的时间点接近，相比于表面张力对气泡行为的影响，黏性力的改变似乎不会大幅度影响气泡发生振荡的频率.气泡上升特征速度U值受气泡形状变化的影响也发生改变，相同Re下，速度U达到幅值的时间点与Ar幅值的时间点相对应.然而，不同Re下的U值大小差异并没有出现前文所述的类似情况，此时Re越大，气泡Ar值越大，U值也越大，这里考虑是黏性力作用影响，Re较大的情况下黏性力较小，气泡更容易受惯性力影响快速上升，Re较小时，周围液体较大的黏性力阻碍了气泡的上升.

3 结论

本研究对外加电场作用下单个气泡在流体中上升的动力行为进行了数值模拟研究.利用水平集方法捕捉气液两相界面.通过Re，Bo和Boe与气泡长径比及上升速度之间的关系研究了Re<40，Bo<30 范围内的液体黏度、表面张力和电场力对气泡运动及变形的影响.研究得到以下结论.

未加电场的情况下，当液体黏度和表面张力较大时，气泡几乎维持球状；反之，由于表面张力，惯性力和黏性力的不平衡，气泡在初始阶段的变形和加速上升.在达到力平衡后，气泡形状波动微小并稳定上升.

外加电场对气泡形状和上升速度有显著影响.气泡沿着电场方向发生拉伸变形，Boe越大，电场强度越大，沿电场方向拉伸更明显，变形程度更大，且由于施加电场力后各作用力之间相互竞争，气泡形状出现竖直和水平方向交替的摇摆不定的振荡形式，且Boe越大，振荡幅度越大.随着气泡上升，气泡振幅随时间而减小.此外，发现竖直方向的匀强电场加速了气泡的上升.电场作用下，保持其他参数不变，随着液体黏度的减小，Re增大，气泡上升初期对电场的响应速度更快，气泡初期与后期变形程度相差更大，且气泡上升速度更大，气泡形态出现更为明显的振荡形式，更难达到稳定状态.液体黏度较大时，气泡响应电场力发生变形后很快维持稳定形态，基本不发生波动，液体较强的黏性力也减缓了气泡的上升速度.同一Re下，随着表面张力系数增大，Bo减小，气泡不稳定振荡更频繁.气泡上升速度受气泡形状变化的影响，发生相同频率的不稳定振荡，且长径比Ar越小，气泡上升速度U越大.