气泡夹带对壁面油膜流动特性的影响

2022-12-19童宝宏张国涛胡晓磊

计算力学学报 2022年6期

童宝宏，陈武，张国涛，胡晓磊

(1.安徽工业大学机械工程学院,马鞍山 243032；2.特种重载机器人安徽省重点实验室,马鞍山 243032)

1 引言

含气泡液体流动普遍存在于自然界及工程领域中，如水轮机工作时造成的空化现象[1]、水面破碎产生的气泡波浪[2]和进行废水处理的气泡[3]等。已有研究表明气泡的存在具有一定的危害性，尤其在油-气润滑系统工作过程中。油-气润滑系统工作时，干燥压缩空气将微米级的油滴传输至润滑点形成润滑油膜。相关研究表明[4]，在高速气流输送润滑油过程中，润滑油与管壁之间容易形成微气室，气室内的气体在压力作用下破裂为许多气泡，游离在润滑油膜内的气泡会降低油膜承载能力，破坏油膜的完整性，造成润滑不良并最终损坏摩擦表面[5]。因此，考察含气泡油膜沿壁面流动的动态演化过程，对于深入研究油膜的润滑性能和查明含气泡油膜流动力学机制具有重要的意义。

油膜的形成与流动过程本身是一个复杂的非定常行为[6]。广大学者对沿倾斜壁面上的液膜流动特性及形态演变过程开展了一系列研究工作，通过实验观测和数值模拟考察了壁面倾斜角度、壁面微结构和雷诺数等对液膜波动特性的影响。研究发现，增大壁面倾斜角度和改变壁面形状会降低液膜的稳定性[7-9]；减少入口流量会导致液膜出现溪流[10]和滴流[11]等现象，破坏液膜的完整性。气流是油-气润滑过程中不可缺少的组成，已有研究采用数值模拟方法深入考察了气流对液膜厚度及流动速度的影响。结果表明，逆向气流趋向于降低和平坦化区域中的液体速度，同向气流减薄液膜厚度[12-14]。

在油膜实际流动过程中，气泡的形成、破裂及射流总是伴随着与周围流体强烈的耦合作用[15-18]。对气泡靠近自由液面运动的研究发现，初始气泡压力和位置对自由液面的形态及破裂形式具有不同程度的影响[19,20]。气泡破裂过程会改变周围流体的流动形式，气泡周围流体的压力分布也更为复杂。相关学者对液体中气泡破裂过程产生的压力变化开展研究，研究结果表明，气泡初始位置影响了壁面压力大小和压力作用范围[18]。

国内外学者已对液膜流动行为进行了系统且深入的研究，并对液膜的流动速度、形态、厚度及不同物理量对其影响进行相应考察，取得了丰富的研究成果，但涉及含气泡油膜流动的动态行为研究鲜有报道。因此，为了深入探究气泡夹带影响下润滑油膜的动态演化过程及内在流动机理，本文采用VOF方法对建立的含气泡油膜沿倾斜壁面流动的模型进行三维数值仿真模拟，考察气泡及其破裂行为对油膜的形态和流动特性的影响，进一步查明含气泡润滑油膜流动的内在机理。

2 数值方法及计算域

2.1 控制方程

倾斜壁面上含气泡油膜流动是典型的气液两相流现象，自由液面的形状一般不是规整的几何面，且含气泡油膜流动过程中，其相界面的位置随着空间和时间不断发生变化。而求解此类问题的关键在于如何将相界面进行离散化，以及如何精确追踪相界面随时间的变化。目前流体体积法VOF(volume of fluid)广泛应用到追踪自由相界面的数值模拟中。关于VOF模型计算的控制方程为

连续性方程

(1)

(2)

VOF方法框架中，在网格上定义了一个额外的场变量，即体积分数函数γ(x,y,z)，为了指示每一个计算单元内一种流体的体积分数,取值范围在0～1。当函数值为1时，表示计算单元内仅有一种流体填充，γ为0时，则由另一种流体填充，介于0～1时，表示相界面必须穿过该单元。体积分数方程为

∂γ/∂t+·(υγ)=0

(3)

在目标域中任意点的流体密度和粘度为两种流体体积分数的加权平均值，即

(4,5)

本文同时考虑到气相剪切力和表面张力的影响，在数值计算求解过程中将气相剪切力和表面张力合并成同一个动量源相添加到动量方程中，即

(6)

式中σ为表面张力系数，ki为界面曲率。

固液之间因表面张力而产生的力为壁面粘附力，CSF模型中[21]，对这种力的处理方法不是通过设置壁面处的边界条件，而是直接将其并入到表面张力动量源相中，此时表面法向量为

(7)

2.2 计算方法可行性验证

为检验VOF方法模拟含气泡油膜沿壁面流动的可行性，根据文献[22]的实验方法和物性参数，建立无气泡液膜沿壁面流动的三维模型。图1为物理模型和边界条件。考虑到液膜沿倾斜壁面流动呈明显的对称性，本文计算时仅建立了1/2模型，对称轴为y=0的xoz平面，并将其设置为对称边界条件，模型与水平面夹角固定在45°，采用无滑移壁面边界条件，油液从进口流入，出口流出，液相入口设定为速度入口边界条件。

图1 几何模型及边界条件

本文选取文献[22]的部分工况进行CFD计算，表1为根据文献[22]的液膜参数而采用的数值参数。

表1 数值模拟中液膜物理参数

表1中液膜Re计算公式为

Re=4Q/(bυ)

(8)

式中Q为体积流量(m3/s)，b为板宽(m)，υ为运动粘度(m2/s)。

数值计算时，采用有限体积法对控制方程进行离散，压力速度耦合的求解选用PISO方法，压力求解采用PRESTO!方法，气液界面插值选取几何重构方法，动量求解选择二阶迎风格式。为了获得更好的收敛性，使用10-5s～10-4s的可变时间步长和恒定全局库朗数0.25进行瞬态模拟。收敛标准是绝对残差，监控的变量是出入口液相质量流量和，残差波动在不同的时间步之间没有明显变化且出入口液相质量流量之和小于10-4kg/s时，表示计算达到准稳态，此时计算模拟停止。

取数值模型中部分区域作为研究区域，将数值计算出的液膜平均厚度δs与文献[22]实验测量的平均膜厚δe对比。表2为对比结果。模拟结果和试验之间的相对误差为

δ=(δe-δs)/δe×100%

(9)

表2 平均膜厚

由表2可知，平均液膜厚度随着Re数增加而增加，数值模拟结果与实验结果的吻合性较好，这一规律与文献[22-23]描述的液膜厚度与Re之间的变化规律相同。由表2可知，雷诺数在107附近时，数值结果与实验数据误差最小，Re为15，24和202时，数值模拟结果与实验值的相对误差大于6%，其原因可能是模拟条件和试验条件之间存在细微差别，此外，文献[22]也指出，通过测量数据计算出的入口流量Qave与实际入口流量Q之间有着0.02 LPM～0.3 LPM的误差，进而导致数值模拟时设置的入口流速存在一定的偏差。总体来看，VOF方法能较为准确地模拟无气泡油膜沿壁面流动的动态过程，在此基础上，本文将采用VOF方法进一步数值模拟含气泡油膜沿倾斜壁面流动的动态行为。

2.3 参数选择及模型确定

图2 流动油膜添加气泡

2.4 网格无关性验证

为精确捕捉流动油膜的气-液界面变化和提高计算效率，对气泡附近区域网格进行无关性讨论。网格划分过程中，采用六面体结构网格划分整个计算域，然后对气泡及其周围区域进行局部加密，网格最小尺寸分别为0.07 mm,0.06 mm,0.05 mm,0.04 mm和0.03 mm。定义气泡刚添加到油膜内为初始时刻，即t=0 ms。气泡贯穿整个油膜流动过程，其直径变化在模拟过程中能够直观体现，此外气泡直径也能反映气泡的气-液界面变化过程，因此可以用不同网格尺寸下气泡直径来表明网格质量。

图3给出了气泡破裂前，不同网格尺寸下沿板宽(y)方向xoz平面上气泡最大直径d随时间的变化情况。

图3 不同网格尺寸下气泡最大直径变化对比(y =0 mm)

可以看出，网格尺寸不同，气泡最大直径的变化规律基本相似。网格尺寸对气泡破裂时间有显著影响，其中网格尺寸为0.07 mm和0.06 mm的最大直径和破裂时间与网格尺寸0.04 mm和 0.03 mm 相比存在较大差距。而网格尺寸为0.04 mm和0.03 mm气泡的最大直径随时间变化基本相同，且破裂时间也非常接近。总体来看，为保证数值计算结果具有参考意义的同时考虑计算效率，本文选择最小网格尺寸为0.04 mm的网格对气泡流动区域的油膜进行划分。

3 结果与讨论

数值研究过程中，将气泡中心设置在沿油膜流动方向x=8 mm、板宽方向y=0 mm及垂直于壁面方向z=0.9 mm，此时气泡中心与液膜表面距离为0.93 mm。

3.1 含气泡油膜流动形态分析

图4展示了含气泡油膜流动过程中油膜形态的演化过程。初始时刻，气泡与流体之间的相对速度为零。气泡受粘性剪切力及自身浮力的作用，沿着壁面向下游运动，同时产生远离壁面的上浮运动。油膜受气泡上浮作用影响，在气泡上游的油膜表面出现液面隆起现象。随后气泡与油膜表面接触发生破裂，气泡破裂后，在油膜中形成空腔，空腔在表面张力和浮力作用下逐渐减小，直至形成射流。这与文献[16,25]描述的现象一致，也进一步证明了本文数值模拟的准确性。

图4 含气泡油膜形态演化

图5展示了t=10 ms时有无气泡夹带油膜流动过程中气液界面高度对比。气泡的运动和形变诱发油膜的气液界面发生不同幅度的波动，在气泡上游区域油膜的气液界面波动幅度明显大于下游，另外，气泡上游油膜气液界面有外凸趋势，界面高度增加了0.042 mm，出现波峰；气泡下游油膜气液界面有内凹趋势，界面高度减小了0.012 mm，出现波谷。

图5 油膜流动过程的气液界面(t =10 ms)

为进一步分析气泡存在对油膜形态带来的影响规律，本文统计了气泡破裂前气泡上游油膜气液界面峰值A随时间的变化。如图6所示，纵坐标(A-H1)代表不同时刻含气泡油膜上游气液界面峰值与初始油膜气液界面的差值。随着气泡在油膜内部发生形变与迁移，气泡受到周围流体的粘性剪切力作用；与此同时，气泡对周围流体也具有相同的反作用力，使得气泡上游油膜的气液界面产生波动，油膜气液界面峰值开始逐渐增加；在t=12.5 ms时，有无气泡油膜气液界面差值达到最大值0.048 mm；t>12.5 ms后，气液界面差值没有发生显著变化，而是处于一定幅值范围内的细微波动状态，直至气泡发生破裂，这主要是因为隆起的油膜受表面张力和重力作用需要保持原有的形态。

图6 油膜气液界面峰值随时间变化

3.2 气泡对油膜流动速度的影响

油膜流动表速度和油膜内速度是反映油膜形成速率和传热特性的两个重要参数。图7展示了有无气泡夹带时油膜表速度分布，其中区域I-III分别表示气泡的上游、顶部和下游的油膜。可以看出，气泡邻域内的油膜流动表速度受气泡影响产生扰动，在不同区域内存在速度梯度分布，这主要是因为随着气泡向下游运动，油膜与气泡之间存在速度差，致使流动油膜受到气泡表面粘性阻力作用，阻碍了气泡周围流体的运动。在I-III区域中，气泡对区域III内的油膜流动速度阻碍程度最大，速度峰值减小了0.021 m/s；区域I内的油膜阻碍程度次之,速度峰值减小了0.018 m/s；区域Ⅱ内的油膜流动速度阻碍作用最小，尤其靠近气泡中心位置，其影响程度最小。以上分析说明气泡的存在能够延迟油膜表面流动时间，不利于摩擦表面润滑油膜快速形成。

图7 含气泡油膜表速度分布(y =10 mm，t =18.7 ms)

图8展示了油膜内某点速度随时间变化曲线，可以看出，当6.5 ms

图8 含气泡油膜内某点速度时变曲线(x =17.5 mm，y =0 mm，z =1.65 mm)

图9展示了t=18.7 ms时，沿流动方向上、区域I-III有无气泡夹带时油膜流动各方向的速度分布情况。可以看出，区域I-III油膜减小的表速度主要体现在方向x和方向z的速度Vx和Vz，根据3.1节所述，油膜中夹带的气泡因其上浮运动诱发油膜表面出现隆起现象，受粘性剪切力作用致使油膜流动受阻，原本沿方向x和方向z的速度减小。气泡阻碍的油膜因气泡上浮运动而产生方向z的运动趋势，导致气泡周围的油膜表面高于远处油膜，流经气泡附近的流体发生分流现象，引起沿板宽(y)方向的分速度，沿板宽(y)方向的速度Vy略有增加，这也进一步解释了3.1节所提到的，气泡上游油膜的气液界面峰值在增加到最大值后没有继续增加的原因。

图9 含气泡油膜表面不同方向速度(y =10 mm，t =18.7 ms)

3.3 气泡破裂阶段对油膜压力的影响

气泡在壁面附近破裂时，会引起周围流体产生复杂的物理现象，这将影响着油膜的润滑性能。图10 展示了气泡破裂过程油膜压力分布云图，其中黑色曲线为气液相界面曲线。气泡发生破裂后，油膜中的空腔在表面张力和浮力作用下逐渐减小，直至形成射流。从图10(a,b)可以看出，在气泡破裂前期，气泡内部压力由原来的69.8 Pa减少到4.3 Pa，气体迅速释放，气泡空腔体积快速减小，空腔外围流体随着空腔壁向中心运动，但外围流体由于惯性不易迅速跟上，致使在其邻域形成负压区域[26]，负压最小值为-48.5 Pa。随后空腔受惯性效应减弱，空腔体积收缩速度放缓，如图10(c)所示，空腔下方的负压区域也进一步缩小，而负压区域主要集中在空腔下游，上下游油膜压力差值为38.5 Pa，这主要是由于空腔上下游非对称性流动引起的，空腔上游流体运动方向与空腔收缩速度方向相同，空腔下游流体则相反。在气泡破裂后期，如图10(d)所示，t=32.3 ms时，空腔底部宽度基本消失，空腔类似倒置的圆锥状，此时空腔上下游流体向空腔底部聚集，空腔底部流体形成局部高压区域，该高压区域也是空腔形成射流的主要原因。

图10 气泡破裂过程油膜压力分布

为定量描述气泡破裂阶段对油膜压力的影响情况，选取气泡下方壁面压力来间接分析油膜压力变化规律。图11展示了含气泡油膜破裂过程中气泡下方壁面压力分布，其中压力正负号代表壁面所受压力方向(压力值为负，代表壁面受z轴正方向的力)。t=30.4 ms时，气泡为未发生破裂壁面压力分布；t=30.7 ms时，空腔底部壁面形成一个负压区；t=30.9 ms时，区域负压值继续增大，出现第一个压力峰值，数值为-30.8 Pa；t=31.8 ms时，壁面压力减小；t=31.9 ms时，壁面呈现正负压同时存在的现象；t=32 ms时，与图11(e)中壁面压力相比，壁面处的压力有所增加，壁面所受负压力全部消失；t=32.3 ms时，壁面压力出现第二个压力峰值，数值为 60.3 Pa，这表明气泡破裂过程可能对摩擦表面造成损害；t=32.7 ms时，随着空腔射流演化，壁面压力值逐渐降低。