液滴碰撞壁面是常见的两相流现象，广泛存在于应急喷淋冷却、汽水分离再热、电路印刷等各种工业科技领域。液滴在壁面上的运动行为对许多物理过程有重要影响，例如在喷淋冷却中，需要增大液滴碰壁后的铺展直径D，并减小液滴中心厚度h_c，从而扩大换热面积，促进流体蒸发，利用汽化潜热加强散热^[1-2]。液滴碰撞壁面后的运动状态，主要受流体性质、壁面润湿性、碰撞速度、入射角度等因素的影响^[3-8]。通过研究上述因素与液滴行为的内在联系，有效控制液滴运动，揭示其动力学机理，是研究和应用的共同需要。在近20年里，随着高速摄像技术的应用和普及，液滴碰撞动力学得到了快速发展^[9-12]。

为描述初始直径D₀的液滴以速度U₀垂直碰撞壁面后的形态变化，研究者们用铺展因子β=D/D₀表征液滴直径变化，用无量纲中心厚度h^*=h_c/D₀表征液滴厚度变化，观测β和h^*随无量纲时间τ=t·U₀/D₀变化的情况，重点关注最大铺展因子β_max和最小中心厚度h_min^*。针对上述参数，学者们开展了大量研究。

根据β随τ变化情况，液滴碰撞分为运动、铺展、弛豫和平衡4个阶段^[13]，如图 1所示。在运动阶段(τ < 0.1)，惯性力占主导地位，液滴上部仍保持球状，下部铺展成楔形边缘，文[13, 19-22]指出该阶段可忽略流体性质影响，β随τ呈指数上升，即β ∝τ^b，其中指数b介于0.45到0.57之间，一般取b=0.5；在铺展阶段，受黏性力或毛细力的主导作用，β逐渐增至β_max，然后液滴开始回撤收缩，进入驰豫阶段^[14]；在平衡阶段，液滴可能稳定铺展，也可能发生反弹^[15-16]，或者碎裂飞溅^[17-18]。本文主要针对运动阶段和铺展阶段进行研究。

对于直径小于其毛细长度的液滴碰撞(即ρgD₀² < σ，其中ρ和σ分别是液滴的密度和表面张力，g为重力加速度；本文研究体系均满足该条件)，一般认为β_max主要受We和Re这2个特征数控制^[22]。We=ρD₀U₀2/σ，Re=ρD₀U₀/μ，其中μ是液滴的黏度。

1996年，Pasandideh-Fard等^[23]对水滴碰撞进行实验和数值研究，认为β_max由We、Re和前进接触角θ_a共同决定：β_max=[(We+12)/(3-3cosθ_a+4We/Re^0.5)]^0.5。2009年，Roisman等^[24]通过理论推导，认为在We、Re较大的情况下，有β_max=0.87Re^0.2-0.40Re^0.4We^-0.5。上述预测模型认为We和Re都对β_max产生影响，而Bartolo等^[3]、Clanet等^[14]、Eggers等^[26]利用碰撞因子P=We/Re^b(0 < b < 1)判定液滴铺展由哪种力主导，认为当毛细力主导铺展时，β_max由We决定；当黏性力主导铺展时，β_max由Re决定。

h^*随τ变化分3个阶段^[27]：降低阶段、平台阶段、回升阶段。在降低阶段的初期，受惯性力主导，h^*随τ线性降低^[22]。当U₀较大时，h_min^*接近黏性边界层厚度，h^*保持在h_min^*一段时间即进入平台阶段。随后，受表面张力作用，h^*进入回升阶段。

对于h_min^*，2010年Schroll等^[28]数值模拟了硅油碰撞平面，忽略空气影响，发现h_min^*由边界层厚度决定。2010年，Eggers等^[26]通过数值模拟，认为当Re较大时存在规律h_min^*∝Re^-0.4，但对于Re较小的情况，可以忽略双曲流影响，有h_min^*∝Re^-0.5。2012年，Lagubeau等^[29]通过实验测量水和甘油溶液液滴碰撞平面，认为在100 < Re < 10 000的范围内，h_min^*与Re^-0.4成正比。2017年，Zhu等^[30]通过数值模拟液滴撞击球面，认为当We较小时，表面张力也显著影响液滴厚度。

β_max和h_min^*的多种预测模型如表 1所示。从表 1可以看出，在黏性力区，已有的理论、实验和模拟研究普遍认为β_max主要由Re控制，与Re^0.2成正比，但目前实验数据和模拟验证均较少，已有实验主要针对硅油、甘油溶液，黏度在300 mPa·s以内，缺乏对于更高黏度液体的研究；另外，Eggers等^[26]的模拟研究不针对真实流体，且忽略空气影响和假设壁面完全疏水，与现实情况具有差异。在毛细力区，尽管Eggers等^[26]的模拟结果与Fedorchenko等^[25]的理论相符合，都认为β_max∝We^0.5，但Clanet等^[14]的实验结果却表明β_max∝We^0.25，实验所得We指数大小仅为理论的一半。因此，仍需对β_max进行更多的实验、模拟研究。

针对h_min^*的实验和模拟研究结论主要基于水和甘油溶液，研究流体种类较少；同时，Lagubeau等^[29]的实验结论h_min^*∝Re^-0.4和Du等^[22]的数值结果h_min^*∝Re^-0.5也存在一定差异，尽管2种结论所基于的Re范围有所不同，但Re是否对指数有直接影响以及Re是否为h_min^*的唯一影响因素，有待进一步探讨。

综上所述，针对β_max和h_min^*的研究，不同学者所得结论具有差异，研究的流体种类也较为单一，We和Re的研究范围也各不相同。因此本文采用相场法，对9组真实流体碰撞不同润湿性平面进行数值模拟，将μ和σ的研究范围分别扩展为1~970 mPa·s和20~500 mN/m，并通过改变U₀，获得1≤Re≤18 600、1≤We≤4 660的测试体系，探索流体性质、壁面润湿性等因素在液滴碰撞过程中的作用机理，补充了低Re液滴碰撞数据，探索已有理论的适用性，尝试揭示更具普适性的液滴碰撞机制。

1 相场法模型介绍

多相流数值模拟的技术关键点之一是对相界面进行描述和追踪，常用的界面捕捉技术有VOF法^[31]、水平集法^[32]、界面追踪法^[33]、相场法^[34]等。与VOF法和水平集法相比，相场法对高密度比两相流模拟具有更好的数值稳定性，且能更加真实地体现表面张力的影响^[35]。同时，相场法用有限厚度的扩散层来描述相界面，用界面扩散来体现三相接触线移动，解决了接触线“奇点”问题^[36]。

基于上述优点，相场法在液滴碰撞模拟研究中得到了广泛应用^[36-38]。本文采用有限元计算软件COMSOL Multiphysics^®中的相场法模块，建立2-D轴对称的液滴碰撞模型。

1.1 相场法原理

本文采用的相界面控制方程为Cahn-Hilliard方程^[39]：

$ \partial \phi / \partial t+\boldsymbol{u} \cdot \nabla \phi=\nabla \cdot(\gamma \nabla G). $

(1)

其中：ϕ为相函数，在单相流体内部为常数，而在相间交界层内线性变化；u为流场速度矢量；γ是迁移率，表征相间扩散的时间尺度；G是化学势，代表多相流系统自由能对相函数的变化率^[40]。

$ G=\lambda\left[-\nabla^{2} \phi+\phi\left(\phi^{2}-1\right) / \varepsilon^{2}\right] . $

(2)

其中：λ为混合能密度；ε是相界面厚度，取最小网格尺寸h的一半，即ε=h/2。

界面张力F_st由相函数ϕ的梯度和化学势G共同决定：

$ \boldsymbol{F}_{\mathrm{st}}=G \nabla \phi. $

(3)

将F_st作为体积力加入到控制流场的Navier-Stokes方程中去，并结合连续性方程，计算不可压缩流体的运动情况：

$ \rho(\partial \boldsymbol{u} / \partial t+\boldsymbol{u} \cdot \nabla \boldsymbol{u})=-\nabla p+\mu \nabla^{2} \boldsymbol{u}+\boldsymbol{F}_{\mathrm{st}}+\rho \boldsymbol{g}, $

(4)

$ \nabla \cdot \boldsymbol{u}=0 $

(5)

其中：ρ和μ分别是流体的密度和动力黏度，p是压强，g是重力加速度。

通过求解计算式(1)—(5)，利用相函数获得新的界面位置，得到各相流体的体积分数，对流体密度、黏度等参数进行更新，实现两相流的瞬态模拟计算。

1.2 模型验证

为保证模型的数值稳定性，h和D₀需满足关系：h/D₀ < 0.01^[41]。为研究网格敏感性，设置D₀=2 mm，分别取h为0.014和0.010 mm，建立非结构化网格。用2种网格设置的模型模拟水滴碰撞疏水壁面，所得β~τ曲线如图 2所示，2种网格设置的计算结果相差约0.03%。为节省计算时间，本文取h=0.014 mm。

为验证模型可靠性，建立与Šikalo等^[42]的水滴碰撞实验相同的数值模型：D₀=2.45 mm，Re=4 010，We=90，壁面接触角θ_e=100 °。相场法模拟结果与Šikalo等的实验数据、VOF法模拟结果的对比如图 3所示，其中灰色为实验照片，橘色为VOF模拟结果，蓝色为本文相场法模拟结果。可以看出，相场法模型能较为准确地呈现液滴形态的变化过程。

图 4给出了上述碰撞过程的β~τ曲线。在运动阶段，相场法模型所得β值比实验数据略小；从τ>0.1到液滴铺展至β_max，相场法模型与实验数据整体符合良好，相场法对β_max的模拟精度略高于VOF法；在驰豫阶段，相场法模拟结果与实验数据基本吻合。

2 液体性质和壁面润湿性对液滴碰撞的影响

本文的研究体系如表 2所示，包括9组流体：水、3种质量分数的甘油水溶液、4种型号的PDMS硅油，以及镓基合金。镓基合金是一种常温液态金属，密度约为水的6倍，表面张力高达500 mN/m，黏度与水接近，其特殊的物性为在25℃常温下研究低黏度、超高表面张力和超高密度液滴碰撞提供了依据。PDMS硅油的表面张力约为水的1/3，是一种低表面张力流体，其黏度范围非常广，可以低至与水同一量级，也可以是水的上千倍。

所研究液滴D₀均为2 mm，通过改变U₀(0.1~ 7 m/s)，获得了We变化范围覆盖3个数量级、Re变化范围覆盖4个数量级的模拟结果。同时，为分析壁面润湿性的影响，设置θ_e为90°、150°、180°。

本文采用Clanet等^[14]提出的P=We/Re^0.8来判断液滴铺展的主导力，认为当P>1时，黏性力起主导作用；当P < 1时，毛细力起主导作用。将表 2中的离散数据点绘制到图 5中。可以看到，因为镓基合金和水的黏度很低，这2种流体的碰撞体系完全处于毛细力区；甘油溶液、PDMS硅油作为水基和油基的黏性流体，大部分处于黏性力区。

2.1 铺展因子β随时间变化

数值模拟表 2中部分流体碰撞θ_e为150 °的壁面，运动阶段的β~τ曲线如图 6所示，Re>100的数据用空心点表示，Re < 100的数据用实心点表示，黑色实线为Re>100数据的拟合结果β=2.17τ^0.48(拟合方差R²=0.945)，符合Rioboo等^[13]的结论。

从数据点分布情况可以看出，对于Re>100的液滴碰撞，尽管包含物性比较特殊的镓基合金，其分布也较为集中，基本遵循一致的指数变化规律。Re < 100的液滴碰撞结果则与指数规律β=2.17τ^0.48相差较多，且数据点较为分散。这是因为低Re液滴的黏性力接近于惯性力，而β ∝τ^0.5规律的成因是惯性力主导铺展，所以该规律对Re < 100的情况并不完全适用。

2.2 最大铺展因子β_max 2.2.1 毛细力主导区

如表 1所示，不同学者给出的毛细力区(P < 1)β_max模型存在差异，但都认为与We密切相关，而与Re无关。Fedorchenko等^[25]通过理论推导、Eggers等^[26]利用数值模拟，认为β_max∝We^0.5，但Eggers是在假设壁面完全疏水(θ_e=180°)的前提下得出该结论。Bartolo等^[3]、Antonini等^[5]、Clanet等^[14]通过实验研究，得出β_max∝We^0.25；Clanet等和Bartolo等认为无论壁面润湿性如何，该规律均适用，但Antonini等则认为当We < 200时, 壁面疏水性会导致β_max偏离该规律。

本文在毛细力区研究的流体包括水、60%质量分数甘油溶液、PDMS 10硅油和镓基合金，β_max随We变化的情况如图 7所示。对不同润湿性壁面的模拟结果进行拟合，如图 7中虚线所示：当θ_e=90°时，β_max∝We^0.26(R²=0.921)；当θ_e=180°时，β_max∝We^0.40(R²=0.990)；当θ_e=150°时，β_max∝We^0.29(R²=0.918)，介于前两者之间。因为真实表面无法绝对疏水，所以Clanet等^[14]实验研究结果β_max∝We^0.25更接近本文θ_e=90°的模拟结果，而Eggers等^[26]针对完全疏水表面提出的β_max∝We^0.5更适用于θ_e=180°的情况。

虽然图 7中的4种流体物性差异较大，但模拟结果整体分布在以We为变量的幂函数曲线附近，说明在毛细力区，We是β_max的主要控制因素。从图 7中3条拟合曲线的指数变化可以看出，随着壁面疏水性增强，黏性力对液滴铺展的阻碍作用减弱，液滴的β_max随之增大，导致θ_e=180°时β_max∝We^b规律的指数b最大。

2.2.2 黏性力主导区

如表 1所示，文[14, 22, 25-26]指出在黏性力主导区(P>1)存在规律β_max ∝Re^0.2。本文在黏性力区的研究采用3种质量分数的甘油溶液、3种型号的PDMS甘油，黏度变化范围为10~970 mPa·s，跨越近3个数量级。其中，黏度极高的PDMS 1 000硅油在低速碰撞壁面时有Re=1。

β_max随Re变化情况如图 8所示，图中虚线为碰撞不同润湿性壁面的拟合结果：当θ_e=90°时，β_max∝0.93Re^0.19(R²=0.980)；当θ_e=150°时，β_max∝Re^0.27(R²=0.953)。可以看出，θ_e=90°的数据整体落在了θ_e=150°数据的上方，但随着Re增大，二者的差异逐渐减小。由此可见在黏性力区，液滴碰撞亲水壁面的β_max更大；但随着Re增加，黏性力影响减弱，与黏性力相关联的壁面润湿性影响随之减弱^[43]，此时疏水壁面上β_max迅速增加，导致β_max∝Re^b规律的指数b更大。

当Re趋近于1时，数据点明显偏离拟合曲线，这是因为高黏度液滴低速度碰撞壁面时，入射动能很快通过黏性力耗散殆尽，液滴铺展基本不受惯性力影响，近似于自然铺展，导致β_max与Re关联性减弱。

2.3 液滴中心厚度h^*随时间变化

图 9给出了表 2中4种流体在θ_e为150°的壁面上的h^*~τ曲线。当τ < 0.1时，无论We、Re如何变化，均有规律h^*=1-τ，说明该阶段h^*受惯性力主导。当We较大时，液滴的高入射速度和低表面张力使h^*出现平台期。对于水和镓基合金，当We < 10时，h^*随时间变化具有明显波动性，这说明在低速入射的情况下，高表面张力会导致液滴厚度产生振荡，这与Wang等^[44]在水滴碰撞实验中观测到的现象一致。

如表 1所示，文[22, 24-26, 29]认为h_min^*仅与Re相关。从图 9中甘油溶液和PDMS硅油的曲线可以看出，二者Re相近，h_min^*几乎相同，与已有结论一致；但对于We相近、Re相差较大的水和镓基合金，二者也有相似的h_min^*。因此，h_min^*是否仅与Re相关，需要进一步探讨。

2.4 最小中心厚度h_min^*

如表 1所示，Roisman等^[24]和Fedorchenko等^[25]通过理论推导，Eggers等^[26]通过数值模拟，发现当Re较大时有h_min^*∝Re^-0.4，Lagubeau等^[29]通过实验测量支持了该结论。当Re < 1 000时，Du等^[22]和Eggers等^[26]通过模拟研究，认为此时可以忽略双曲流影响，有h_min^*∝Re^-0.5。也有研究^[30]认为当We较小时，表面张力显著影响液滴厚度。

由于h^*达到h_min^*时，液滴的中心比边缘要薄，因此很难通过拍摄手段获得准确数据^[29]，并且已有结论主要基于水和甘油溶液，体系覆盖范围有限。因此，本文对表 2中9组流体进行碰撞模拟，得到h_min^*随Re变化的情况如图 10所示(水和镓基合金低速碰撞θ_e为90°的壁面时，液滴中心剧烈振荡无稳定模拟结果，故剔除这部分数据)。

本文所研究的碰撞体系大部分满足Re < 1 000，因此与文[22]中规律h_min^*∝Re^-0.5比较。对We≥10、Re < 1 000的数据拟合得到h_min^*=1.20Re^-0.50(R²=0.955)。可以看出空心数据点大多落在该拟合曲线上，本文模拟结果与文[22]的结论一致。从2种形状的空心数据点在拟合曲线上基本重合可以看出，壁面润湿性对We≥10液滴的h_min^*~Re关系无显著影响。

We < 10的实心数据点整体分布在拟合曲线上方，且不同θ_e对应的h_min^*相差较大。当We=1时，液滴碰撞θ_e为150°的壁面的h_min^*比拟合曲线h_min^*=1.20Re^-0.50高出197.53%。这一现象证实了Zhu等^[30]的研究结论：We较小时，表面张力和壁面润湿性对h_min^*产生明显影响。因此本文认为，h_min^*主要受Re控制的理论仅适用于We较大的情况。

此外，当Re接近于1时，h_min^*明显偏离h_min^*∝Re^-0.5规律，其原因与上文所述一致，即高黏度液滴低速碰撞壁面时，液滴更多的表现出自然铺展特性。

3 结论

本文利用相场法，数值模拟了9组流体碰撞不同润湿性壁面的过程，拓展了液滴碰撞模拟研究的流体物性范围，主要结论如下：

1) 在碰撞初期(τ < 0.1)，发现文[13]提出的规律β∝τ ^0.5[13]主要适用于Re>100的液滴碰撞，这是因为低Re液滴的黏性力较强，惯性力不再处于主导地位；h^*的模拟结果符合已有结论h^*=1-τ，说明在碰撞初期，h^*主要受惯性力影响。

2) 在毛细力区，本文模拟结果显示β_max∝We^b，符合前人研究结论；随着壁面疏水性增强，指数b逐渐增大，说明当毛细力主导铺展时，壁面越疏水，黏性力对铺展的阻碍作用越小，β_max越大。在黏性力区，数值结果满足规律β_max∝Re^b，与已有理论相符；当液滴Re较小时，亲水壁面上的β_max更大，但随着Re增大，惯性力增强，壁面润湿性的影响随之减弱，疏水壁面上的β_max迅速增加，导致β_max∝Re^b规律的指数b较大。

3) 当We≥10、Re < 1 000时，h_min^*模拟结果与文[22]提出的规律h_min^*∝Re^-0.5相符，但We < 10时会偏离该规律，壁面润湿性和表面张力对h_min^*产生明显影响。其机理在于低We液滴的毛细力较强，高表面张力和壁面疏水性减弱了液滴的铺展能力，所以不能再仅依据Re去预测h_min^*。

4) 对于以往鲜有研究的高黏度PDMS 1 000硅油和超高表面张力的镓基合金，尽管这2种流体物性较为特殊，但其碰撞过程仍主要受Re和We影响。但当PDMS 1 000硅油以低速碰撞壁面时，Re接近于1，此时液滴的入射动能很快通过黏性力耗散殆尽，惯性力作用可忽略不计，液滴在壁面上近似于自发铺展，导致β_max和h_min^*主要由液滴初始动能和壁面润湿性决定，偏离结论2和3中与We、Re之间的幂函数关系。