陈铁军，房 媛，田 辉*

(1.河北省承德市兴隆县政府办公室,河北 承德 067300；2.河北石油职业技术大学 机械工程系,河北 承德 067000)

液滴下落过程普遍存在于多相流系统、喷墨打印、涂装工艺及降雨过程等领域，在此过程中的运动、变形甚至破碎现象一直深受学者的关注。Eggers[1]基于系统的理论分析获得了液滴下落变形及破碎的非线性动力学规律。Stone[2]着重分析了表面张力影响下低雷诺数时下落液滴的变形规律。其研究显示：低雷诺数时液滴下落过程满足连续性方程及斯托克斯方程，液滴变形的程度主要受两相间表面张力作用影响。Leal[3]发现液滴的尾迹受惯性力影响显著表现出非线性特点。

由于下落过程中液滴与周围空气间存在悬殊的密度、粘性差异，使得相界面变形等迁移特性的数值研究困难重重，国内外专家学者相继提出了一系列求解算法，如MAC[4]方法、Front Tracking[5]法、VOF[6]法、Level Set[7]法等。然而以上单一算法中仍存在系统质量守恒性差、捕捉界面不够准确(税利的界面细节被算法抹平)等问题。CLSVOF[8]、Particle Level Set[9]等混合方法的提出有效的改善了以上问题。本文在经典粒子水平集法(Particle Level Set)的基础上通过改进距离函数修正方式实现高精度两项界面迁移过程的捕捉并保证系统质量守恒性，为解决液滴下落等复杂工程问题提供有效途径。

1 无量纲控制方程组

本文通过一套无量纲不可压缩流体控制方程组[10]求解气液两相流动(如方程组1-4所示)，不考虑气液两相间质量传递；气液两相相应区域以不同的物性参数表征，无质量虚拟粒子随流场迁移模拟相界面的迁移运动。计算中均使用无量纲参数，则描述变密度不可压缩流动问题的无量纲控制方程组可描述为张量形式：

(1)

(2)

(3)

(4)

2 改进粒子水平集法

粒子水平集法将拉格朗日思想及欧拉思想良好的结合在一起，为多相流相界面迁移运动问题的准确求解提供了有力的工具。然而算法中粒子对距离函数的修正策略局限于网格内部，当待修正网格点、逃逸粒子及界面不在同一网格内的时候，现有方法所获得的修正值存在一定偏差[11]。

本文不再采用根据网格点距离函数的正/负来对局部距离函数的求解进行分类。这里借助如图1所示的虚线l来判定待修正网格点局部距离函数的求解方式。虚线l过投影点P′且其斜率由此处界面切向量决定。l将计算域分为两部分，若待修正点与原始水平集法确定的界面位于同一侧，则局部距离函数可通过方程(5)获得；若待修正点与原始水平集法确定的界面分别在虚线l的两侧，则局部距离函数通过方程(6)获得。

φP(x)=Sp(rp+|x-xp′|)

(5)

φP(x)=Sp(rp-|x-xp′|)

(6)

(7)

3 控制方程离散求解方式

4 液滴下落过程的数值模拟

液滴在重力作用下的下落及浮力作用气泡的上升过程受以下控制参数影响：两相的密度ρ液和ρ气，两相的动力粘性系数μ液和μ气，重力加速度g，两相间表面张力系数σ，空间位置xi以及时间T。基于π定理，通常选取以上8个变量所构成的5个独立变量来描述此问题[14]。本文采用的控制参数为：(1)密度比ρ液/ρ气；(2)动力粘性系数比μ液/μ气；(3)无量纲时间T；(4)Re数(Reynold number)；(5)Eo数(Eotvos number)。

(8)

(9)

(10)

液滴下落过程的数值模拟在1×2的矩形计算域中进行，采用结构化均分网格，计算域四边均为壁面无滑移边界条件，CFL数设为0.1。液体密度设为ρ液=1 000 kg·m-3，动力粘性系数设为μ液=1.137×10-3Pa·s保持与水一致，初始全场速度均为0。

首先对后续数值模拟采用的网格进行无关性验证。分别采用三种结构化网格：64×128、128×256、256×512，液滴直径为D，初始于(0.5,1.7)，Re=200，We=100。图2给出了三种情况下液滴质量(面积)守恒特性。采用64×128网格时液滴的质量亏损严重，液滴发生破裂后质量亏损超过10%，计算过程CPU耗时325 s；采用128×256网格时液滴的质量亏损较小，液滴发生破裂前后质量亏损始终保持在1%以内，图3给出了三种网格下液滴下落过程典型形态。图中可见三种情况液滴变形情况极其相似，且网格数的增加使液滴表面更光顺。图3(a)中可见液滴边缘粗糙的棱角。计算过程CPU耗时4 167 s；采用256×512网格时液滴的质量亏损几乎可以忽略，但消耗大量计算资源，CPU耗时129 527 s。兼顾求解精度及计算效率，本文均采用128×256的网格对液滴的相关流动过程进行数值研究。

在壁面无滑移边界条件的影响下，近壁面的流动区域存在明显的速度梯度，流动表现为剪切流。如图2所示，液滴破碎后在近壁面区域受壁面作用，下落过程出现明显的左右摆动现象。本文通过数值模拟到壁面不同距离液滴的下落过程，分析总结壁面作用对下落液滴运动、变形的影响。

图4为到壁面不同距离的三个液滴下落过程，到壁面距离分别为D、1.5D、2D。图中可见壁面对液滴下落过程影响显著，越靠近壁面此作用越突出。三种情况液滴下落初期均首先向靠近壁面一侧摆动，但很快便开始向另一侧摆动，到一定极限时摆动方向调转，如此反复直至下落到底。以图4(a)的情况为例，液滴初始时刻圆心到壁面的距离为D，图5给出了此过程中三个典型时期液滴变形、相对流线及周围压力的分布(红色表示高压区域，蓝色表示低压区域)。计算域内绝对速度场减去液滴的平均速度可得到液滴的相对速度场，相对流线便是在此相对速度场中绘出的。液滴下落初期水平及垂直速度均较小，液滴左右两侧变形差异很小整体呈现为椭圆形，相对流线紧贴液滴下侧，而上侧出现两个强度较小且旋向相反的两个涡，液滴左右两侧压力基本一致。随着液滴的下落流场趋于复杂，液滴下落中期相对速度场中出现了三个影响范围较大的涡。液滴呈现出类似鸡蛋的形状且大头在左侧，其周围流场完全被旋向相反的两个涡控制。此时液滴顶部左右两侧压力不平衡，左侧略低。液滴下落后期情况与下落中期恰好相反，鸡蛋形的液滴大头在右侧，液滴顶部右侧压力略低，另外相对速度场中出现了四个涡。

5 结论

本文在推导求解气液两相流无量纲控制方程组的基础上，针对传统粒子水平集法距离函数修改方式缺陷，提出了一种基于逃逸粒子投影法判断待修正点与界面位置进而完善距离函数修正的方法。在此基础上针对壁面作用下液滴在空气环境中的下落过程进行了数值模拟，数值计算结果显示：下落液滴在近壁面剪切流影响下发生摆动，而液滴的左右摆动诱发其尾迹区内产生涡；另外强度及影响范围不断增加的涡又会影响液滴的后续运动；涡产生后其强度及影响范围不断发生变化然而涡的中心稳定于液滴晃动方向发生调转的区域。本文所发展的数值求解方法具有较高的界面捕捉精度，适用于求解存在高密度差、高粘性差、互不相容的多相流界面迁移运动。