Research on the Dynamic Contact Angle Model for the Droplet Impact Process
-
摘要: 基于计算流体力学(CFD)模拟液滴冲击壁面,对于理解液滴在固体壁面铺展的动力学行为有重要的意义,可以为超疏水结构设计及防除冰涂层开发提供技术支撑,其中的难点在于如何在模型中准确刻画接触线及动态接触角的演化过程. 总结了四种典型的动态接触角模型,从理论上分析了其应用范围,借助FLUENT中的UDF功能,将动态接触角模型应用于壁面边界条件. 首先对液滴冲击光滑壁面的动力学过程进行了数值模拟研究,通过定量分析液滴形态的各项参数变化并与实验结果对比表明,Seebergh动态接触角模型更适用于模拟低毛细数下液滴的运动,Kistler模型与Jiang模型应用范围更广并且可以较准确地描述高毛细数下液滴的运动. 随后基于Kistler动态接触角模型,对液滴在微结构表面的冲击与铺展过程进行了仿真研究,发现应用动态接触角模型会导致液滴内部流场在表面张力起主导作用的阶段内发生变化,并且在平衡状态下液滴接触角的模拟值与理论值相近.Abstract: The simulation of droplet-wall impact process based on computational fluid dynamics (CFD) is of great significance for understanding the dynamic behavior of droplets spreading on the solid wall, and can provide technical support for the design of superhydrophobic structures and the development of anti-icing coating. The difficulty lies in how to accurately describe the evolution process of the contact line and the dynamic contact angle in the model. Herein, 4 typical dynamic contact angle models were summarized, and their application ranges were analyzed theoretically. With the UDF function in FLUENT the dynamic contact angle model was applied to the wall boundary conditions, and the dynamic process of droplet impact on smooth wall was numerically simulated. The quantitative analysis of the changes of droplet shape parameters and the comparison with the experimental results show that, the Seebergh dynamic contact angle model is more suitable for simulating the motion of droplets with lower capillary numbers. The Kistler model and the Jiang model are more widely used and can accurately describe the motions of droplets with higher capillary numbers. Then, based on the Kistler dynamic contact angle model, the impact and spreading processes of droplets on the microstructure surface were simulated. It is found that, the application of the dynamic contact angle model will lead to the change of the internal flow fields of droplets with the surface tension playing a dominant role, and the simulated droplet contact angle value in equilibrium is close to the theoretical value.
-
0. 引言
低能耗、高效率一直是航空防除冰设备研发的目标,采用超疏水微结构表面与其他加热方式结合的复合防除冰设备具有非常高的研究价值[1]. 对液滴冲击固体壁面的物理过程进行研究可以为航空复合式防除冰设备的开发提供参考,因此保证数值模拟的准确性是很有必要的[2-3]. 液滴在固体表面的动力学行为复杂多样,包括撞击、反弹、振动等已有较多研究[4-6]. 当液滴在空气环境中与壁面接触时,会形成一条固体、液体和气体三相相交的接触线,接触线在衬底上运动,直到总能量达到最小值时停止,即热力学平衡状态,在该动态条件下,接触线被定义为动态接触线[7-9]. 液-气界面和液-固界面在动态接触线处的夹角称为动态接触角θD[10]. 当动态接触线达到平衡状态时,动态接触角转变为静态接触角θ0[11].
为了准确描述液滴在壁面上的运动,国内外已有许多学者对液滴冲击问题中接触线的运动展开了实验和理论研究. 在早期研究中,人们通过大量的实验去寻找接触线运动速度和宏观接触角之间的关系,并拟合出了经验公式,一直沿用至今,通常用毛细数Ca=vμ/σ来表征黏性力与表面张力对液体流动情况的影响,其中v为接触线速度,μ为黏性系数, σ为表面张力系数. Hoffman[12]通过实验测量了在黏性力和界面力做主导的条件下,5种液体以10-5 < Ca < 10在毛细管中流动时所对应的接触角,其中大部分实验的毛细数处于10-5~10-1之间,而对应的接触角涵盖了0°~180°. 根据实验结果,他提出了连续的Hoffman函数,以此来描述动态接触角和毛细数之间的关系. Kistler以Hoffman方程为基础,考虑接触角滞后效应,通过理论分析修正并完善了Hoffman函数,提出了应用较为广泛的Kistler模型[13]. Jiang和Slattery等[14]在Hoffman实验的基础上,对其实验数据进行最小二乘法拟合得到了经验模型. Bracke等[15]通过实验测量了将固体条带拉入液体池的过程中动态接触角的变化,并根据实验数据得出了经验公式. Seebergh等[16]总结前人工作,通过实验测量了低毛细数与平衡接触角下动态接触角的演变规律,提出了经验模型. 此外,也有不少学者通过力学分析建立了动态接触角的理论模型. Blake等[17]从分子动力学的角度出发,通过描述分子在接触线附近的运动,考虑了表面张力与摩擦力的作用,用分子动力学参数表示接触线速度与接触角之间的关系. Cox[18]则应用流体力学的理论进行分析,建立了滑移长度与动态接触角的关系式. 然而,理论模型部分参数的确定存在难度,如微观与宏观的滑移长度以及分子动力学参数等,相较于宏观经验模型来说应用更加复杂.
随着计算流体力学(CFD)的不断发展,越来越多的学者选择通过数值模拟的方式研究液滴冲击过程. 模拟自由表面流动的两个主要挑战是规定固体表面上的边界条件,特别是移动接触线上的边界条件,以及考虑润湿效应,即将接触角包含在模型方程中. 在大多数研究工作中考察的无量纲参数包括Weber数We=ρD0V02/σ,Reynolds数Re=ρD0V0/μ,其中ρ为液体密度,D0为液滴初始半径,V0为液滴初始撞击速度. Šikalo等[19]使用VOF方法,引入Kistler动态接触角模型,模拟了液滴冲击壁面的过程,模拟得到的水滴形态变化与实验观察结果吻合较好. Xie等[20]使用局部网格动态加密技术,模拟了液滴以不同速度冲击亲水与疏水表面的过程,研究了液滴铺展系数、高度系数、接触角以及表面积的变化. 然而,少有学者提及模拟计算中动态接触角模型的选取方法以及各个模型的特点.
对于液滴冲击超疏水微结构表面的研究同样受到了国内外学者的广泛关注. 最经典的粗糙表面接触理论模型主要有Wenzel模型[21]与Cassie-Baxter模型[22]两种,前者认为液滴最终完全浸润到微结构内部,后者认为微结构与液滴之间存在空气,两种模型均为多功能微结构表面的设计提供了参考. Tuteja等[23]在粗糙微结构的基础上考虑表面曲率,设计出了对低表面张力液体呈现疏水性的微结构表面. Liu等[24-25]设计了圆锥形微结构表面,并且在高We数下较为稳定地实现了液滴的煎饼弹跳,大幅减少了液滴与表面的接触时间. Du等[26]通过VOF方法对液滴冲击环形槽状微结构表面的过程进行了研究,考察液滴最大铺展状态并建立了最大铺展系数的预测模型. Yang等[27]对液滴冲击方柱表面的过程进行了数值模拟,通过研究液滴内部流场与压力的分布情况定性分析了液滴出现铺展与回缩反弹行为的原因. 采用已有的数值模拟研究多数应用静态接触角作为壁面边界条件,然而对于形式各样的微结构表面,何种工况下该应用动态接触角模型并未被讨论.
本文在前人研究的基础上,采用数值模拟的方法,使用Kistler、Jiang、Bracke和Seebergh四种动态接触角模型,对液滴冲击动力学过程进行了模拟研究. 对比分析了不同参数条件下动态接触角模型对于液滴形态的影响,总结了不同模型的特点及使用条件,并给出了模型改进的建议. 此外,使用动态接触角代替本征接触角对液滴冲击微结构壁面进行数值模拟,分析了液滴运动过程中内部流场的变化,对比分析了应用动态接触角模型与静态接触角模型计算结果的差异.
1. 动态接触角模型
1.1 宏观模型
本小节首先介绍Kistler、Jiang、Bracke和Seebergh四种动态接触角模型,然后从理论上对比了这四种模型描述动态接触角的差异性.
模型1 Kistler模型[13]
Kistler模型考虑接触角滞后效应,即接触线速度大于0时采用前进角θadv,小于0时采用后退角θrec,将动态接触角定义为接触线速度的函数,其中自变量为毛细数Ca和Hoffman函数对应前进或后退角的反函数值,具体形式如下:
$$ \theta_{\mathrm{D}} =f_{\mathrm{H}}\left[C a+f_{\mathrm{H}}^{-1}\left(\theta_{0}\right)\right], $$ (1) $$f_{\mathrm{H}} =\arccos \left\{1-2 \tanh \left[5.16\left(\frac{x}{1+1.31 x^{0.99}}\right)\right]^{0.706}\right\}, \quad x=C a+f_{\mathrm{H}}^{-1}\left(\theta_{0}\right), $$ (2) $$\begin{align*} \theta_{0} & = \begin{cases}\theta_{\text {adv }}, & v>0, \\ \theta_{\text {rec }}, & v <0\end{cases} .\end{align*} $$ (3) 模型2 Jiang模型[14]
由于Hoffman的实验主要考虑黏性力和界面力的作用,Jiang模型基于Hoffman的实验数据,通过最小二乘法拟合得到,因此二者适用范围相同. 所以理论上在不考虑重力和惯性力的条件下,Jiang模型适用于任何几何形状,动态接触角大小仅取决于所涉及的材料的性质和接触线的位移速度. 具体形式如下:
$$ \begin{equation*} \frac{\cos \theta_{0}-\cos \theta_{\mathrm{D}}}{\cos \theta_{0}+1}=\tanh \left(4.96 C a^{0.702}\right) \end{equation*}. $$ (4) 在模型计算结果与实验数据的对比中,对于低速(We < 10-3)和小接触角(θ0 < 70°)的实验数据误差在10%以内. 然而对于惯性力更高的情况该模型的表现较差.
模型3 Bracke模型[15]
Bracke等在实验中将固体条带拉入液体池,观测并记录了接触线速度与接触角的关系,得到了接触角余弦值与条带移动速度的平方根之间的线性关系,通过采用不同材料进行实验,确定了斜率和静态接触角余弦值的线性关系. 针对较低的毛细数范围(Ca < 10-1)总结出与Jiang模型形式类似的经验公式,具体如下:
$$ \begin{equation*} \frac{\cos \theta_{0}-\cos \theta_{\mathrm{D}}}{\cos \theta_{0}+1}=2 C a^{0.5} \end{equation*}. $$ (5) 模型4 Seebergh模型[16]
Seebergh等总结了Jiang等和Bracke等的工作,整理出了动态接触角经验模型的结构形式,并针对更低的毛细数(Ca < 10-3)和平衡接触角(θ0 < 70°)改进参数,得到了以下模型:
$$ \begin{equation*} \frac{\cos \theta_{0}-\cos \theta_{\mathrm{D}}}{\cos \theta_{0}+1}=2.24 C a^{0.54} . \end{equation*} $$ (6) 在上述四种模型中,除Kistler模型以外,其余模型均在Hoffman实验的基础上进行细化,在不同的毛细数范围内修正拟合参数,但是并未包含接触角滞后信息.
1.2 模型对比
为了对比四种动态接触角模型的差异,本文计算了平衡接触角为10°,90°以及160°情况下动态接触角的计算值随毛细数的变化规律,如图 1所示. 考虑实验工况以及经验公式的定义域,Bracke模型只考虑Ca < 0.1的情况,Seebergh模型考虑Ca < 0.01的情况,Jiang与Kistler模型则考虑了更广的毛细数范围. 从计算结果可以看出,动态接触角随着毛细数的增加而增大,并且对于高毛细数的情况,即使平衡接触角不同,动态接触角都将趋于180°. 在平衡接触角为10°且毛细数较小的情况下,Seebergh和Bracke模型计算结果相近并且大于Kistler与Jiang模型. 但随着平衡接触角的增大,四种模型的计算结果在小毛细数时越来越接近,且Kistler模型的动态接触角略小于其他三种模型. 当Ca>0.1时,只有Kistler和Jiang模型可以计算出结果,可以看出,Kistler模型计算的动态接触角要小于Jiang模型,并且Jiang模型在Ca>0.01后,动态接触角会更迅速地增大到最大值.
2. VOF方法及计算设置
2.1 VOF模型
本文采用Hirt和Nichols[28]提出的VOF模型来追踪界面. 假设气相和液相都是不可压缩的,则质量和动量守恒方程以及体积分数的输运方程如下:
$$ \nabla \cdot \boldsymbol{u}=0, $$ (7) $$ \frac{\partial}{\partial t}(\rho \boldsymbol{u})+\nabla \cdot(\rho \boldsymbol{u} \boldsymbol{u})=-\nabla p+\nabla \cdot\left[\mu\left(\nabla \boldsymbol{u}+\nabla \boldsymbol{u}^{\mathrm{T}}\right)\right]+\rho \boldsymbol{g}+f_{\mathrm{vol}}, $$ (8) $$ \frac{\partial \alpha_{q}}{\partial t}+\boldsymbol{u} \cdot \nabla \alpha_{q}=0, $$ (9) 其中p为压力;g为重力加速度;fvol为考虑表面张力影响的源项;αq为相体积分数,q=1为气相,q=2为液相. 式中流体参数是指控制体内的加权有效物性参数,如下所示:
$$ \rho =\sum\limits_{q} \alpha_{q} \rho_{q}, $$ (10) $$ \mu =\sum\limits_{q} \alpha_{q} \mu_{q} . $$ (11) 采用Brackbill等[29]提出的连续表面张力(continuum surface force,CSF)模型将表面张力转换为体积力,由此在动量方程中引入体积力的源项为
$$ \begin{equation*} f_{\text {vol }}=\sigma \kappa \nabla \alpha_{2} \frac{\rho}{\left(\rho_{1}+\rho_{2}\right) / 2} , \end{equation*} $$ (12) 其中$\nabla \alpha_{2}$是液体体积分数的梯度,κ为曲面曲率,由单位法线的差值定义:
$$ \kappa=-\nabla \cdot \boldsymbol{n}, $$ (13) $$\boldsymbol{n}=\frac{\nabla \alpha_{2}}{\left|\nabla \alpha_{2}\right|} . $$ (14) 与壁面相邻单元的表面法线表示为
$$ \begin{equation*} \boldsymbol{n}=n_{\mathrm{w}} \cos \theta_{\mathrm{W}}+t_{\mathrm{w}} \sin \theta_{\mathrm{W}}, \end{equation*} $$ (15) 其中nw是与壁面垂直的单位向量,tw是与壁面相切的单位向量,θW为壁面处的接触角.
2.2 计算设置
本文采用FLUENT中的VOF模型分别对液滴冲击平面和微结构表面的动力学过程进行模拟,计算域的物理模型及网格如图 2所示. 为了在保证计算精度的同时尽可能节约计算资源,采用二维旋转轴对称模型进行模拟. 设置Y轴为对称轴,X轴为壁面,其余边界设置为压力出口. 整个计算域的尺寸为8 mm×8 mm的正方形区域,采用结构化网格划分,并且对贴近对称轴与壁面的空气区域进行加密,取液滴运动区域网格尺寸为20 μm×20 μm.
数值模拟设置选择瞬态计算过程,采用Pressure-Based求解器,压力方程采用PISO算法求解,体积力加权离散. 动量和能量控制方程采用二阶迎风格式,气液界面的追踪采用Geo-Reconstruct格式离散,瞬态格式采用一阶迎风格式求解. 时间步长设置为1 μs,计算液滴运动总时长为15 ms.
动态接触角模型通过FLUENT中UDF功能实现,其中接触线速度根据单元格的实际速度计算:
$$ \begin{equation*} \boldsymbol{u}=\left(u_{\text {cell }} \cdot n_{t}\right) \frac{n_{t}}{\left|n_{t}\right|} . \end{equation*} $$ (16) 接触线的运动方向由式(17)计算:
$$ \begin{equation*} \mathcal{\varepsilon}_{\mathrm{sign}}=\boldsymbol{n}_{\mathrm{cl}} \cdot \boldsymbol{n}, \end{equation*} $$ (17) 当εsign为正时,接触线回缩;当εsign为负时,接触线铺展.
本文所涉及到的计算工况如表 1所示,为了包含工程问题中常见的毛细数范围,使用水和甘油两种液滴作为研究对象,所研究的毛细数范围为10-3 < Ca < 101,并考察了101 < Re < 103,100 < We < 102范围内的液滴冲击现象. 水滴的参数为ρ=998 kg · m-3,μ=10-3 N · s · m-2,σ=0.073 N · m-1,甘油液滴的参数为ρ=1 220 kg · m-3,μ=0.116 N · s · m-2,σ=0.063 N · m-1.
condition wall liquid diameter D0/mm impact velocity v/(m/s) contact angle θ/(°) We Re Ca 1 plane water 2.45 1.64 105,95 90.0 4 010 0.023 2 glycerol 2.45 1.41 97,90 94.0 36 2.590 3 water 2.47 0.21 63.4 1.5 523 0.003 4 water 2.47 0.21 100 1.5 523 0.003 5 water 2.47 0.21 160 1.5 523 0.003 6 glycerol 2.45 1.41 17,13 94.0 36 2.590 7 glycerol 2.45 1.41 160 94.0 36 2.590 8 micro structure water 2.00 1.00 165,155 27.0 2 000 0.014 8 water 2.00 0.50 165,155 13.5 1 000 0.007 10 water 2.00 1.00 95,85 27.0 2 000 0.014 11 water 2.90 0.34 165,155 4.7 984 0.005 3. 计算结果及分析
3.1 低毛细数
选取Šikalo等[30]和Chen等[31]的实验结果作为数值计算的参考. 图 3展示了在工况1中应用Kistler模型计算所得结果与实验的对比情况. 左图为实验结果,从图中可以看出模拟所得液滴形态随时间的变化与实验结果吻合度较好. 液滴冲击过程可以分为铺展、回缩和反弹三个阶段. 在液滴接触壁面后,随着时间的推移,液滴在惯性力的作用下开始铺展,同时表面张力和黏性力也在抑制三相线的移动,直到达到最大铺展后,液滴边缘形成突起. 随后在表面张力的主导作用下,三相线后退,液滴开始回缩,表面能转化为向上的动能,使液滴产生反弹的趋势,当动能足够大时,液滴发生反弹.
为进一步对比不同动态接触角模型的计算结果,对液滴形态进行量化分析,其中铺展系数定义为当前液滴直径与液滴初始直径之比,即D/D0,高度系数定义为当前液滴轴线上液膜的高度与液滴初始直径之比,即H/D0,对时间进行无量纲化处理为τ=tV0/D0. 图 4展示了在工况1中应用不同接触角模型计算得到的液滴铺展系数和高度系数随时间变化的情况,在此工况下,液滴冲击壁面的最大毛细数Ca=0.023,这时Seebergh模型不适用,故只对Kistler、Jiang与Bracke模型进行验证. 从图 4中可以看出,三种模型的模拟结果在液滴铺展阶段基本相同,这是因为三种模型均考虑了前进角的影响,并且在该毛细数下,三种模型对应动态接触角的变化范围相近. 计算所得铺展系数略大于实验测得数据,但绝对误差值处于0.1至0.3之间,相对误差值在τ>0.1时保持在10%以内. 差异主要体现在液滴回缩反弹阶段,从回缩阶段开始,Kistler模型计算所得铺展系数略大于Jiang与Bracke模型. 由于Kistler模型针对回缩阶段的后退角进行修正,使得液滴的接触角更小,因此回缩速度更慢,铺展系数也就相对较大. 并且由于会出现如图 3中t=7.80 ms所示的液滴分离后重聚的情况,导致在回缩结束前的短时间内Kistler模型计算所得的高度系数为0. 此外,我们可以发现虽然三种模型在回缩结束时的铺展系数是不同的,但其最大高度系数基本相同,这是由撞击过程中的惯性力决定的,在该工况下黏性损耗占比较小. 其中Jiang模型模拟的液滴达到最大高度的时间要早于Bracke和Kistler模型,而后两者在反弹阶段的模拟结果几乎是相同的,由于Bracke和Kistler模型在毛细数处于10-2~10-1计算所得接触角更接近,所对应的液滴形态也就更相似.
图 5为在工况3中应用四种动态接触角模型模拟所得的铺展系数与高度系数随时间变化的情况. 由于此工况下毛细数很小(Ca=0.003),所以Seebergh模型也可以计算出结果. 从计算结果可以看出,在更低毛细数的情况下,液滴达到最大高度系数时,Jiang模型与Bracke模型的计算结果相近,都要早于并且大于Kistler模型. 同时我们可以发现,随着毛细数的降低,前三种模型计算结果在液滴回缩与反弹阶段更加大于实验值,即误差增加. 但是Seebergh模型的计算结果与实验值吻合得比较好,在达到最大铺展系数时误差约为12%左右. 同时高度系数也与实验值表现出较高的符合性,其最小高度系数和最大高度系数出现时间更早. 相比于其余三种模型,Seebergh模型针对毛细数处于10-5~10-3的范围内接触角与三相线速度的关系式进行拟合参数的优化,可以更准确地反应低毛细数下表面粗糙度引起的黏滑效应,因此Seebergh模型更加适用于低毛细数下的液滴动力学行为模拟.
图 6展示了四种模型在工况4下的计算结果,相较于工况3增大了平衡接触角,计算了液滴冲击疏水表面的过程. 从图中我们可以看出,四种模型计算所得铺展系数和高度系数的变化趋势与撞击亲水表面的情况一致. 由于壁面的平衡接触角增大,液滴最大铺展系数和高度系数出现的时间点前移,铺展系数降低,高度系数增加. 虽然此工况下没有实验数据作为参照,但我们根据工况1与工况3计算的结果可以推断在低毛细数的情况下,四种模型计算所得的动态接触角是略小于实际值的,而图中所示Seebergh模型的计算结果相较于Kistler、Jiang、Bracke模型来说铺展系数更小,高度系数更大. 因此我们可以判断,即使在较大平衡接触角的情况下,Seebergh模型对于低毛细数的液滴运动模拟结果也更加贴近实际情况.
为了进一步了解平衡接触角对于低毛细数液滴冲击运动的影响,图 7展示了工况3、4、5下四种模型的计算结果. 从图中可以看出,四种模型模拟得到的铺展系数和高度系数随时间变化趋势基本一致. 随着平衡接触角的增加,液滴的最大铺展系数逐渐降低,最大高度系数与最小高度系数均逐渐增加,并且最小高度系数的发生时间也有所提前. Kistler、Jiang、Bracke模型模拟结果基本相同,平衡接触角变化导致的液滴形态改变程度相差不大. Seebergh模型对于平衡接触角改变的响应程度相较于前三种模型更高,增大接触角导致液滴铺展系数的差值达到0.3至0.4.
3.2 高毛细数
选用了黏度更高的甘油液滴作为研究对象,模拟了其冲击不同平衡接触角壁面的过程. 在工况2、6、7中,毛细数Ca=2.59相对较大,Bracke和Seebergh模型并不适用,因此只对Kistler和Jiang模型进行验证. 图 8展示了在工况2、6、7下计算所得的接触角随时间的变化过程. 从图中可以看出,在冲击亲水表面的情况下,Kistler和Jiang模型的模拟结果基本相同,与实验值的变化趋势相同,但是数值上存在误差,并且接触角变化率更低. 当接触角增大后,二者的差异逐渐体现出来,在相同的平衡接触角条件下,动态接触角变化趋势基本相同,但是Jiang模型计算的数值要大于Kistler模型,在工况2中尤为明显. 当平衡接触角进一步增大,壁面表现为超疏水,差值减小. 当平衡接触角较小时,从图 1中我们可以了解到二者通过理论计算所得的平衡接触角差距很小,这也就解释了为何模拟结果相近. 当平衡接触角增大后理论值差距增加,在90°时,两者对应毛细数Ca=2的动态接触角差值达到10°左右,在160°时,这个差值仍然存在但降低到了5°左右. 从结果来看,平衡接触角越小,计算所得的动态接触角数值随时间降低的速度越快,对于超疏水表面,动态接触角在整个冲击过程中的变化趋势为在平衡接触角附近波动.
图 9展示了Kistler和Jiang模型在工况6下计算的甘油液滴铺展系数和高度系数随时间变化的情况. 从图中可以看出在低平衡接触角、高毛细数的情况下,两种模型计算的结果基本相同,由于甘油液滴黏度较高,且壁面为超亲水壁面,液滴并未出现明显回缩与反弹的情况. 在τ=1之前,模拟所得铺展系数略大于实验值,这与图 8中动态接触角的变化规律相符,在液滴铺展阶段实验测得的动态接触角大于计算值,表面更加疏水,这也就导致了液滴的铺展系数更小. 当液滴铺展到最大半径的时候,实验值减小并且低于计算值,因此所表现出的铺展系数也趋于稳定并最终高于模拟值. 综上所述,对于高毛细数液滴冲击问题的模拟,Kistler和Jiang模型的表现近乎相同,计算所得动态接触角都与实验值有一定偏差,并且接触角的偏差导致的液滴形状的改变存在一定的滞后性,绝对误差值小于0.2,可以较准确地模拟高毛细数液滴冲击问题.
3.3 微结构表面
3.3.1 二维轴对称环形槽状微结构
本文通过工况8、9、10分析了应用动态接触角模型后,液滴撞击环形槽状微结构过程中液滴内部流场的变化. 边界条件与微结构的尺寸如图 10所示,采用二维轴对称模型,其中微结构的宽度a、间隔b与高度h均为200 μm. 根据Cassie-Baxter模型,表观接触角可表示为$\cos \theta_{0}=\phi\left(1+\cos \theta_{\mathrm{Y}}\right)-1$,其中$\phi$是固相分数,定义为界面投影区域中的液-固接触面积的分数. 本文所选微结构参数对应$\phi=0.511, \theta_{\mathrm{Y}}=160^{\circ}$,以此来模拟经过化学处理后具有微观二级结构的超疏水表面的表观接触角,并将其作为宏观微结构表面的本征接触角. 根据理论公式计算得到处于Cassie-Baxter状态的液滴表观接触角为165.7°,应用动态接触角模型计算液滴稳定时的表观接触角为168°,与理论值相近.
经过模拟发现,是否应用动态接触角模型主要影响其内部流场的变化情况. 图 11为工况8中应用动态接触角模型后在5 ms内从铺展阶段到回缩阶段液滴内部流场的变化. 可以看出在液滴接触壁面后,接触点的速度迅速减小,边缘三相线速度较快,并向外沿壁面铺展. 在2 ms时可观察到液滴边缘壁面处的流速接近0,并且产生向内的涡流,说明此时液滴外沿已开始回缩,而靠近对称轴的速度矢量则表示液滴中心处仍在向外流动,同时中心低速区在扩大,这就导致了在3 ms时出现了速度差较大的局部凹陷区域,而在此区域的下方无壁面支撑,表面张力不足以维持液滴的形态,因此在4 ms时液滴分离为三部分. 由于分离前液滴外沿部分已存在回缩并向上的速度,因此在5 ms时液滴重新聚集并呈现出反弹的趋势.
应用固定接触角作为边界条件与应用动态接触角模型模拟得到的液滴形状与运动情况基本相同,然而内部流场却存在偏差,这种偏差会影响液滴内部的传热与传质过程[32]. 图 12(a)为应用静态接触角模拟的3 ms时液滴内部流场情况,可以发现位于微结构间隙处的液滴局部存在涡流. 在铺展阶段结束时,液滴外缘靠近壁面处流场趋于静止,远离壁面处流速较快并向液滴中心流动. 在回缩阶段开始后,由表面张力主导液滴外缘向中心移动,壁面处流速逐渐增大,而此时液滴中心还存在向外铺展的区域,对向流动相互作用,因此在间隙处形成涡流. 然而当应用动态接触角模型进行模拟后,靠近壁面处的流动变化相对平缓,高速区与低速区过渡平滑,间隙处的涡流明显减弱,如图 12(b)所示.
同样的规律可以在降低We数后更明显地观察到. 图 13为工况2模拟得到的3 ms时的液滴内部流场,同样可以观察到液滴靠近壁面间隙处存在涡流,并且在应用动态接触角模型时有所减小. 在维持We相同的情况下降低本征接触角至90°,从图 14中同样可以观察到由于接触角滞后效应导致近壁面处的流速减缓,阻碍了涡流的形成.
综上,针对二维轴对称环形槽构型,应用动态接触角模型会对液滴冲击微结构壁面的模拟结果产生影响,主要体现在内部流场的变化上. 在固相分数较大的情况下,应用动态接触角模型会导致壁面处的流速变化更加平缓,微结构间隙处不易形成涡流.
3.3.2 三维方柱微结构
然而,针对方柱微结构和圆锥微结构等复杂三维构型,固相分数降低后,动态接触角模型对液滴形态与液滴内部流动的影响尚不明确,本文参考文献[24-25]中的实验工况,通过数值模拟进行了验证,建立了如图 15所示的三维仿真模型. 采用了1/4对称模型进行计算,其中宽度a=100 μm,间隔b=200 μm,高度h=800 μm,固相分数为0.4. 计算了文献[24]中的触发液滴煎饼弹跳的工况,对模型与方法的准确性进行验证. 图 16为不同时刻下液滴形态与实验结果的对比,可以看到,应用Kistler动态接触角模型可以较为准确地模拟三维复杂构型中的液滴形态演化过程.
随后针对文献[25]中的实验工况,考虑到Jiang模型、Bracke模型与Seebergh模型都未包含接触角滞后信息,而Jiang模型的适用范围更广、鲁棒性更好,因此分别采用静态接触角、Kistler动态接触角模型与Jiang动态接触角模型进行了数值模拟研究. 液滴撞击初始We=4.7,壁面本征接触角为160°,前进与后退角分别为165°与155°,详细信息见表 1中工况11.
图 17展示了液滴形态的量化结果,从铺展系数与高度系数随无量纲时间的变化过程中我们可以看出,在针对固液接触面积更小、流固耦合作用更突出的复杂构型时,动态接触角模型对于液滴宏观形貌的演化影响较小. 此外,由于计算工况对应的毛细数小于0.01,且本征接触角较大,此时动态接触角的变化区间更小,并且不同模型在相同毛细数下对应的接触角数值相近,因此体现到液滴形态上的差异不明显.
图 18展示了不同接触角模型在4 ms铺展阶段以及7.5 ms回缩阶段的液滴轮廓以及内部流场信息. 可以观察到无论是在惯性效应占主导的铺展阶段,还是表面张力占主导的回缩阶段,液滴整体轮廓以及内部流动并未因接触角模型的不同产生明显差异. 这可能归因于高本征接触角对应的动态接触角数值变化幅度较小,同时较低的固相分数进一步弱化了壁面接触角的影响,二者同时作用导致液滴的动力学行为对小范围的接触角变化不敏感.
综上,针对具有复杂构型的微结构壁面,是否选用动态接触角模型需要根据壁面的固相分数、液滴撞击速度以及本征接触角的大小进行判断. 当固相分数较大时,液滴与微结构壁面的相互作用更类似于水平表面,应用动态接触角模型对液滴的运动过程将产生较为明显的影响. 此外,滞后角作为壁面重要参数对于液滴动力学行为的影响较为复杂,不同壁面条件下接触角模型的适用性还需进一步验证.
4. 结论
对于低毛细数的情况(Ca < 0.1),四种模型得到的液滴形态变化的过程具有相同的趋势,均表现出与实验结果相符的特征. 其中Kistler、Jiang与Bracke模型计算得到的结果基本相同,均略大于实验结果. Seebergh模型在低毛细数的工况下与实验结果符合较好. 随着平衡接触角的增大,Seebergh模型表现出更加明显的变化,铺展系数的降幅显著高于其他三个模型,说明动态接触角随平衡接触角的变化呈非线性的关系. Seebergh模型相较于其他三个模型来说更适用于低毛细数的情况,Kistler、Jiang与Bracke模型仍需要根据更多实验数据进行毛细数相关的系数修正,以此满足低毛细数情况下的应用.
对于高毛细数的情况(Ca>0.1),Kistler和Jiang模型可以进行模拟,二者的模拟结果与实验相比均存在一定偏差,这与其高毛细数下对应的动态接触角的非线性变化有关. 在毛细数Ca>0.1后二者都迅速增大到接近180°,虽然动态接触角的理论值之间存在小于10°的偏差,但体现在液滴形态上的差距并不明显. 接触角变化所带来的液滴形态的变化并非实时的,需要进一步研究与这种滞后性相关的影响因素,才能为现有模型的改进提供理论支撑.
将动态接触角模型应用到液滴冲击微结构壁面的模拟中. 针对固相分数较大的二维轴对称简单构型,应用Kistler动态接触角模型计算所得液滴形态变化与应用静态接触角模型的模拟结果相近,然而内部流场却有所区别,应用动态接触角模型会导致壁面处的涡流减弱,可以预见这种改变将影响液滴运动过程中的传热与传质过程. 而针对三维复杂构型采用动态接触角模型,由于滞后角较小、本征接触角较大、固相分数较小等原因,液滴形态与内部流动均未受到明显影响.
-
表 1 计算工况
Table 1. Calculation conditions
condition wall liquid diameter D0/mm impact velocity v/(m/s) contact angle θ/(°) We Re Ca 1 plane water 2.45 1.64 105,95 90.0 4 010 0.023 2 glycerol 2.45 1.41 97,90 94.0 36 2.590 3 water 2.47 0.21 63.4 1.5 523 0.003 4 water 2.47 0.21 100 1.5 523 0.003 5 water 2.47 0.21 160 1.5 523 0.003 6 glycerol 2.45 1.41 17,13 94.0 36 2.590 7 glycerol 2.45 1.41 160 94.0 36 2.590 8 micro structure water 2.00 1.00 165,155 27.0 2 000 0.014 8 water 2.00 0.50 165,155 13.5 1 000 0.007 10 water 2.00 1.00 95,85 27.0 2 000 0.014 11 water 2.90 0.34 165,155 4.7 984 0.005 -
[1] ZHU Y T, WANG Z L L, LIU X L, et al. Anti-icing/de-icing mechanism and application progress of bio-inspired surface for aircraft[J]. Transactions of Nanjing University of Aeronautics and Astronautics, 39(5): 542-554. [2] 林贵平, 卜雪琴, 申晓斌, 等. 飞机结冰与防冰技术[M]. 北京: 北京航空航天大学出版社, 2016.LIN Guiping, BU Xueqin, SHEN Xiaobin, et al. Aircraft Icing and Anti-icing Technology[M]. Beijing: Beihang University Press, 2016. (in Chinese) [3] 张旋. 过冷水滴的结冰与碰撞及其耦合特性研究[D]. 北京: 清华大学, 2019.ZHANG Xuan. Research on freezing and impact processes of supercooled water droplet and their coupling characteristics[D]. Beijing: Tsinghua University, 2019. (in Chinese) [4] 王凯宇, 庞祥龙, 李晓光. 超疏水表面液滴的振动特性及其与液滴体积的关系[J]. 物理学报, 2021, 70(7): 076801.WANG Kaiyu, PANG Xianglong, LI Xiaoguang. Oscillation properties of water droplets on a superhydrophobic surface and their correlations with droplet volume[J]. Acta Physica Sinica, 2021, 70(7): 076801. (in Chinese) [5] PANG X L, DUAN M, LIU H, et al. Oscillation-induced mixing advances the functionality of liquid marble microreactors[J]. ACS Applied Materials & Interfaces, 2022, 14: 11999. [6] 严裕, 娄钦, 陈家豪. 双液滴在具有接触角滞后性微通道内的运动行为研究[J]. 应用数学和力学, 2023, 44(3): 304-318. doi: 10.21656/1000-0887.430165YAN Yu, LOU Qin, CHEN Jiahao. Lattice Boltzmann study on the motion of dual droplets in microchannels with contact angle hysteresis[J]. Applied Mathematics and Mechanics, 2023, 44(3): 304-318. (in Chinese) doi: 10.21656/1000-0887.430165 [7] 焦云龙, 刘小君, 逄明华, 等. 液滴平壁铺展过程中的滞后效应及力学机制研究[J]. 应用数学和力学, 2016, 37(1): 14-26. doi: 10.3879/j.issn.1000-0887.2016.01.002JIAO Yunlong, LIU Xiaojun, PANG Minghua, et al. Study of contact angle hysteresis at moving contact lines based on CFD simulation and mechanical analysis[J]. Applied Mathematics and Mechanics, 2016, 37(1): 14-26. (in Chinese) doi: 10.3879/j.issn.1000-0887.2016.01.002 [8] LI X G, WANG Y Q, YANG Y, et al. Dynamic behavior of droplets under interfacial jamming of nanoparticles[J]. Applied Physics Letters, 2018, 113: 133702. doi: 10.1063/1.5045775 [9] MOHAMMAD K A, SUSZYNSKI W J. Physics of dynamic contact line: hydrodynamics theory versus molecular kinetic theory[J]. Fluids, 2022, 7(10): 1-19. [10] GANESAN S. On the dynamic contact angle in simulation of impinging droplets with sharp interface methods[J]. Microfluidics and Nanofluidics, 2013, 14(3/4): 615-625. [11] YOUNG T. An essay on the cohesion of fluids[J]. Philosophical Transactions of the Royal Society of London, 1805, 95: 65-87. doi: 10.1098/rstl.1805.0005 [12] HOFFMAN R L. A study of the advancing interface[J]. Journal of Colloid and Interface Science, 1975, 50(2): 228-241. doi: 10.1016/0021-9797(75)90225-8 [13] VOINOV O V. Hydrodynamics of wetting[J]. Fluid Dynamics, 1977, 11(5): 714-721. doi: 10.1007/BF01012963 [14] JIANG T S, SOO-GUN O H, SLATTERY J C. Correlation for dynamic contact angle[J]. Journal of Colloid and Interface Science, 1979, 69(1): 74-77. doi: 10.1016/0021-9797(79)90081-X [15] BRACKE M, DE VOEGHT F, JOOS P. The kinetics of wetting: the dynamic contact angle[J]. Progress in Colloid & Polymer Science, 1989, 79: 142-149. [16] SEEBERGH J E, BERG J C. Dynamic wetting in the low capillary number regime[J]. Chemical Engineering Science, 1992, 47(17/18): 4455-4464. [17] BLAKE T D, BRACKE M, SHIKHMURZAEV Y D. Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[J]. Physics of Fluids, 1999, 11(8): 1995-2007. [18] COX R G. The dynamics of the spreading of liquids on a solid surface, part 1: viscous flow[J]. Journal of Fluid Mechanics, 1986, 168: 169-194. [19] ŠIKALO Š, WILHELM H D, ROISMAN I V, et al. Dynamic contact angle of spreading droplets: experiments and simulations[J]. Physics of Fluids, 2005, 17(6): 062103. [20] XIE P, DING H B, INGHAM D B, et al. Analysis and prediction of the gas-liquid interfacial area for droplets impact on solid surfaces[J]. Applied Thermal Engineering, 2020, 178: 115583. [21] WENZEL R N. Resistance of solid surface to wetting by water[J]. Industrial & Engineering Chemistry, 1936, 28(8): 988-994. [22] CASSIE A B D, BAXTER S. Wettability of porous surfaces[J]. Transactions of the Faraday Society, 1944, 40(10): 546-551. [23] TUTEJA A, CHOI W, MA M L, et al. Designing superoleophobic surfaces[J]. Science, 2007, 318(5856): 1618-1622. [24] LIU Y H, MOEVIUS L, XU X P, et al. Pancake bouncing on superhydrophobic surfaces[J]. Nature Physics, 2014, 10: 515-519. [25] MOEVIUS L, LIU Y H, WANG Z K, et al. Pancake bouncing: simulations and theory and experimental verification[J]. Langmuir, 2014, 30: 13021-13032. [26] DU J Y, WANG X, LI Y Z, et al. Maximum spreading of liquid droplets impact on concentric ring-textured surfaces: theoretical analysis and numerical simulation[J]. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 2021, 630: 127647. [27] YANG C J, CAO W R, YANG Z. Study on dynamic behavior of water droplet impacting on super-hydrophobic surface with micro-pillar structures by VOF method[J]. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 2021, 630: 127634. [28] HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39(1): 201-225. [29] BRACKBILL J U, KOTHE D B, ZEMACH C. A continuum method for modeling surface tension[J]. Journal of Computational Physics, 1992, 100(2): 335-354. [30] ŠIKALO Š, MARENGO M, TROPEA C, et al. Analysis of impact of droplets on horizontal surfaces[J]. Experimental Thermal and Fluid Science, 2002, 25(7): 503-510. [31] CHEN B, ZHANG Y H, DAI Z F, et al. Experimental research on the dynamics of a train of droplets impacting, from droplets to liquid film, continuity and inheritance[J]. Energy, 2022, 256: 124670. [32] 章振宇, 张宸玮, 张鹏. 小韦伯数下液滴撞击光滑壁面的数值模拟[J]. 工程热物理学报, 2021, 42(12): 3296-3303.ZHANG Zhenyu, ZHANG Chenwei, ZHANG Peng. Numerical simulation of droplet impacting on free slip wall under small Weber number[J]. Journal of Engineering Thermophysics, 2021, 42(12): 3296-3303. (in Chinese)