动态润湿现象广泛地存在于自然界和工业应用中,比如自然界中荷叶上的水滴滴落[1],工业中的增材打印[2]、喷墨印刷[3]等,对大自然和人类的生活生产有着深刻的影响.润湿性是材料表面的重要特性之一,体现了固体与液体表面的相互作用.度量润湿性的一个主要参数是接触角,它是指液体相界面切向与固体壁面切向所形成的夹角.长期以来,为了实现对液滴的操控,研究液滴在固体表面的润湿行为特性,是力学、生物、材料等学科致力研究的重要方向.由于该问题的跨尺度特性和实验技术的限制,润湿现象一直是相关研究工作的难点,发展动态润湿现象的数值模拟技术具有重要的学术和应用价值.
接触角的研究有着很长的历史,早在1805年,Young[4]就提出了描述静态接触角的经典Young方程.但是在随后相关的研究中,学者们遇到应力奇异性问题,这是移动接触线与固壁无滑移边界条件的矛盾所致,又称为Huh-Scriven佯谬[5].为了解决应力奇异性的问题,滑移模型[6]、前驱膜模型[7]和扩散界面模型[8]被相继提出.基于这些模型,学者进一步对微观接触角和宏观接触角的关系开展了理论研究,得出了与实验相契合的接触角模型[9-12],精确的接触角描述成为可能.
润湿现象具有明确的相界面,可以采用多相流模型进行描述.目前,针对多相流动数值模拟,学者们提出了多种界面捕捉方法,包括front tracking方法[13]、volume of fluid方法[14]、level set方法[15]和相场方法[16].基于自由能理论,相场方法在描述相界面时,不仅考虑了相界面的拓扑变化,而且考虑了界面应力的效应,使其在表面张力的计算上更具优势.在移动接触线的模拟中,相场方法也具有优势.在相场方法中,接触线的移动可视为相界面的扩散所导致,其相当于隐式施加滑移边界条件[17], 从而解决接触线附近的奇异性.基于上述特点, 相场方法已经得到了快速的发展和广泛的应用, 如热毛细流模拟[18]、 非Newton液滴的撞击铺展问题模拟[19]等.
综上所述,动态润湿现象模拟的难点主要是动态接触角的合理描述和接触线附近的奇异性.对于前者,本文引入Yokoi动态接触角模型[20]对接触角进行描述;对于后者,本文采用相场方法来解决接触线的奇异性问题.本文工作总结为:基于相场方法,结合动态接触角边界条件,发展考虑动态润湿的两相流模型,并通过OpenFOAM开源平台对多相流模型进行编程实现;然后对动态润湿现象进行模拟研究,并对比研究不同接触角模型的润湿行为特征.
相场方法是一种典型的扩散界面方法,它具有一定宽度的界面厚度,这个厚度用人工给定的参数ε来表征,本文取为网格尺度.对于两相不相溶流体A和B,流场的相分布用任意一相流体的体积分数C来表示.若A的相分布为C,则流体B的相分布则为1-C.相分布随时间的演化由对流Cahn-Hilliard方程[21]来控制:
·(uC)=·(Mφ),
(1)
其中,u为流动速度,M=M0|C(1-C)|为迁移率,M0为迁移率系数.φ为化学势,其定义式如下:
(2)
对于无质量流入的边界,化学势φ的边界条件如下:
nw·φ=0,
(3)
其中,nw为边界的法向方向矢量,指向单元内.对于一维问题,当相分布处于平衡态时,φ=0,进而可计算得相的分布为
(4)
在自由能理论中,表面张力的作用可以视为单位表面面积的额外自由能,对于平衡态下的界面可以得到等式:
(5)
其中,根据以上两个公式,不考虑切向力的表面张力可表述为势形式[22],表达式为
C.
(6)
基于相场方法,固壁上的接触角润湿边界条件[23]可表示为
nw·
(7)
对于静态润湿条件,接触角θ为常数;对于动态润湿现象,其接触角θ随接触线速度ucl的变化而变化.一种简单的接触角与速度的关系是准动态接触角模型(quasi-dynamic contact angle model,QDCA model):
(8)
其中,ucl为接触线在壁面上的移动速度.这个模型考虑了动态润湿现象中液滴具有不同的前进角θadv和后退角θrec.该模型在早期的研究中被广泛地使用[24-25],但是其难以准确地反映真实的动态湿润现象.为了精确描述动态接触角,文献[9-10,12]基于水动力学理论,得出了动态接触角的关系θ3~Ca,Ca=μu/σ为毛细数.进一步,Yokoi等[20]结合实验,提出了一种动态接触角模型 (Yokoi’s model):
(9)
其中,μ为黏度系数,θe,θmda和θmdr分别代表静态接触角、最大前进角和最小后退角,ka和kr分别是前进接触角和后退接触角的变化系数.
对于等温不可压两相流动,控制方程分布为不可压缩的N-S方程和连续性方程:
P+·(μ(u+uT))+Fs+Fg,
(10)
·u=0,
(11)
其中,Fg=ρg为重力项,Fs为表面张力项,表达式为方程(6).对于两相流体,其相界面的运动控制方程采用Cahn-Hilliard方程,即方程(1).方程中的平均密度ρ和平均黏度μ定义如下:
ρ=ρ1C+ρ2(1-C), μ=μ1C+μ2(1-C).
(12)
本文的数值模型基于OpenFOAM开源平台进行程序实现.在OpenFOAM中,输运偏微分方程采用有限体积法进行离散,求解变量储存在控制体中心.由于界面拓扑变化与速度场之间的强耦合效应,迭代方法的选择对该模型的数值求解至关重要.本文采用以下迭代过程进行求解,从时间步tn到tn+1的求解步骤如下.
1) 采用二阶Runge-Kutta方法求解Cahn-Hilliard方程(方程(1))来更新相分布场C,求解过程中的速度场u为时间步tn的值:
(13)
其中 Z(t,C,u)≡-·(uC)+·(Mφ).
2) 得到Cn+1后,利用方程(12)对物性参数进行更新,得到ρn+1和μn+1,以及对表面张力项Fs进行计算,得到
3) 基于新的物性参数和表面张力项,采用PISO迭代算法[26]对N-S方程组(方程(10)和方程(11))进行求解,得到新的速度场un+1和压力场pn+1.
图1 接触角边界条件示意图
Fig. 1 The schematic diagram of the contact angle boundary condition
上述所涉及的离散和边界条件的实施细节可参见文献[26].
下面对动态接触角边界条件施加的相关细节进行说明,主要是关于接触线速度的计算(图1).
在接触线问题的数值模拟中,接触线速度的计算是重要因素.对于固壁边界,一般采用无滑移边界条件,即ub=0,所以移动接触相的运动描述不能直接采用边界上的速度值ub,需要采用其他方法对接触线的速度进行计算.
本文采用接触线附近控制单元的速度值ucell来估计接触线的速度值:
(14)
由于相场方法是一种扩散界面方法,所以具有变量连续过度的一定宽度的相界面.基于这点,本文认为在接触线附近的变量场也具有连续变化的特点,也就是接触线的速度与附近网格的速度值应该是相近的,所以可以利用邻近网格的速度值估计接触线的速度值.
首先,基于静态接触角问题验证程序接触角边界条件(方程(7))引入的正确性,随后程序的动态接触角引入方法相同,仅仅是接触角大小的区别.对放置在平直壁面上的具有静态接触角的液滴进行模拟,并将模拟结果与相应理论解进行对比.
图2 液滴静态接触角示意图
Fig. 2 The schematic diagram of the static angle of a droplet
首先,考虑忽略重力(g=0)的情况下,液滴在表面张力的作用下,逐渐变形至平衡状态,这个平衡状态将由液滴与壁面的接触角θ决定.如图2所示,半球形液滴初始时刻半径为R0,其平衡状态下高度和铺展半径分别为H和L.根据质量守恒和几何关系,可得铺展前后几何形状关系的解析表达式:
(15)
(16)
数值模拟中,取轴对称坐标系下边长为80的正方形计算域.液滴的初始半径R0为30,其相关的液滴和环境空气的无量纲物性参数为ρd=ρg=1,μd=μg=0.1,σ=0.001.迁移率系数取为M0=0.001.模拟中,接触角θ分别取为30°,60°,90°,120°和150°,计算结果如图3所示.从图中的液滴形状可以看到,随着接触角逐渐增大,液滴从亲水状态转变为疏水状态,其对应的铺展半径变小.图4展示了液滴铺展的理论解与计算结果的对比,可以看出两者基本一致,其数值误差在1%左右.
图3 不同接触角下,模拟所得的平衡态液滴形状 图4 数值解与理论值对比
Fig. 3 The equilibrium shapes of droplets Fig. 4 Comparison between the numerical with different contact angles and theoretical solutions
图5 θ=60°时,不同Eo数下,液滴高度与 图6 θ=120°时,不同Eo数下,液滴高度与
理论解对比图 理论解对比图
Fig. 5 Comparison of numerical and theoretical Fig. 6 Comparison of numerical and theoretical solutions of the droplet height as a function solutions of the droplet height as a function of the Eo number for θ=60° of the Eo number for θ=120°
接下来,考虑具有重力(g≠0)的情况下液滴的铺展状态.在重力的作用下,除了接触角外,液滴的平衡状态还由重力和表面张力之间的竞争关系决定.为了表征重力与表面张力的相对大小,这里引入无量纲Eötvös数Eo=ρgR2/σ.当Eo数远远小于1时,表面张力的作用对液滴的铺展起主导作用,液滴呈现圆弧的形状.当Eo数远远大于1时,则重力的作用对液滴的铺展起主导作用,液滴此时不能维持圆弧的形状.对于这两种情况,该问题具有理论解[27],分别为
(17)
(18)
在此算例中, 取直角坐标系下的计算域为[0,4R0]×[0,2R0].液滴的初始半径为R0=0.01,其相关的无量纲物性参数为ρd=1 000,ρg=1,μd=μg=0.1,σ=0.05.迁移率取为M0=0.000 1.模拟中,接触角θ分别取为60°和120°,以表征液滴的亲水状态和疏水状态.改变重力加速度g的大小,Eo数也随之变化.计算结果如图5、6所示.
图中结果展示了液滴平衡状态下的高度和Eo数之间的关系,图5和图6分别是接触角为60°和120°的情况,求解的Eo数范围为0.01~10 .当Eo≪0.1时,液滴高度与忽略重力下的理论解方程(17)相吻合;当Eo≫8时,此时重力起主导作用,液滴高度与方程(18)的理论解相吻合.以上结果表明,本文基于相场法引入静态接触角模型及相应程序可正确有效地处理静态润湿问题.
本小节针对文献[20]中的实验进行数值模拟,从而验证本文基于相场法引入动态接触角模型及相应程序的有效性和准确性,同时对比分析不同接触角模型结果的差异.
参照文献[20]中的实验,对具有一定初速度V0的液滴下落撞击壁面的过程进行模拟.在这个过程中,液滴与壁面的接触角θ与接触线的移动速度ucl相关,呈现动态润湿的特性.如图7所示,一个液滴以V0=1 m/s的初速度撞击壁面,初始半径为R0=1.14 mm,相关的物性参数为ρd=1 000 kg/m3,ρg=1.25 kg/m3,μd=1.0×10-3 Pa·s,μg=1.82×10-5 Pa·s,σ=0.072 N/m,重力加速度为g=9.8 m/s2.采用Yokoi动态接触角模型,即方程(9),相关参数为静止接触角θe=90°,最大前进角θmda=114°,最大后退角θmdr=77°,前进角系数ka=9.0×10-9,后退角系数kr=9.0×10-8.
图7 液滴撞击壁面示意图图8 采用不同网格计算所得的铺展半径变化图
Fig. 7 The schematic diagram of a droplet Fig. 8 The evolution of the spreading radius impacting on a wall for different meshes
首先,进行网格有效性验证.采用轴对称坐标系,计算域是边长为0.004 m的正方形区域,分别采用100×100,160×160,200×200和240×240四种网格,计算结果如图8所示.
如图8所示,当网格密度达到200×200后,计算结果已趋于收敛,其最大相对误差小于1%.综合考虑精度和计算资源因素,以下计算基于200×200网格进行数值模拟.
对比文献[20]实验和本文数值模拟结果如图9、10所示.
图9展示了实验与数值模拟结果,在不同时刻具有相同的液滴形态特征,其定性地验证了模型和程序的有效性.液滴铺展半径的对比(图10)定量地表征了模型和程序的有效性,数值模拟结果与实验结果在图9、10都表现出很高的吻合度.在铺展半径的时间演化图中,液滴的铺展情况与实验高度一致,但最大的误差出现在0.025 s时刻附近.结合两图可以看到,液滴在撞击壁面初期(t<0.004 s),铺展半径迅速增大,动能转化为内能和势能,接触角呈现前进角(钝角)的状态.在0.004 s附近,液滴的动能被转化为势能以及被黏性消耗,在表面张力的作用下液滴呈现回缩的趋势,此时接触角从前进角(钝角)转化为后退角(锐角)状态.随后,液滴快速回缩,并经历一段时间的振荡后恢复至平衡状态,在这个过程中接触角也随接触线速度的振荡而振荡,最后稳定在静态接触角附近.
图9 不同时刻下,实验与计算的液滴形态对比图,左侧截图源于文献[20]
Fig. 9 Comparison of droplet shapes between experiment and simulation at different moments
with the left-side photos caputred from ref. [20]
图10 200×200网格的模拟结果与实验[20]的铺展半径随时间的变化对比图
Fig. 10 Comparison of time evolutions of the droplet radius between experiment[20] and simulation for a 200×200 mesh
在上一个算例的基础上,本文进一步采用不同的接触角模型(静态接触角模型和准动态接触角模型)进行模拟和对比,分析不同接触角模型对液滴动态浸润行为的影响.
首先,对比采用静态接触角模型的模拟结果,接触角为θ=90°,结果如图11所示.
图11为静态接触角模型与Yokoi接触角模型(方程(9))的铺展半径时间演化对比图,从图中可看出,静态接触角模型下的接触线在初始铺展阶段(t<0.004 s)和回缩阶段(t>0.004 s)都具有更快的速度,并且在t=0.019 s时,发生了液滴的分离(图12).这些结果说明,在静态接触角模型(θ=90°)下,液滴在整个铺展过程中具有更小的能量耗散,并且更大的动能导致了液滴在回缩过程中出现了液滴的分离现象.此结果与实验现象存在很大的差异,说明静态接触角模型无法正确描述液滴在壁面铺展的实际情况.
图11 不同接触角模型的铺展半径随时间的变化对比图
Fig. 11 Comparison of time evolutions of the droplet radius between different contact angle models
图12 t=0.019 s时,静态接触角 图13 不同动态接触角模型的铺展半径
模型的模拟结果 随时间的变化对比图
Fig. 12 At t=0.019 s, the simulation result of Fig. 13 Comparison of time evolutions of the droplet radius the static contact angle model between different dynamic contact angle models
接下来,对比准动态接触角模型的模拟结果,其中前进角θadv和后退角θrec分别等于上一算例中的最大前进角θmda和最大后退角θmdr.模拟的对比结果如图13所示.
由铺展半径时间演化对比图(图13)所示, 两种接触角模型的模拟结果在初始的铺展阶段一致,但在液滴回缩的过程中具有较大的差异, 并且准动态接触角模型下的液滴铺展经历了较短时间的振荡过程便达到了相对平衡的状态.而通过液滴的形态对比图(图14)可知,液滴在0.004 s时刻具有一致的形态,这与图13中的初始铺展段曲线(t<0.004 s)相对应.在液滴回缩阶段(图14中的t=0.010 s时刻图),准动态接触角模型下的液滴具有更小的接触角,所以在接触线附近表面张力对液滴的牵拉作用更大,具有更显著抑制液滴运动的效果,最终体现为液滴的回缩相对更慢.由于回缩速度较小,动能耗尽所需时间更短,所以液滴所经历的振荡过程时间也较少.该算例结果显示,准动态接触角模型与实验现象依然具有一定的差异,无法准确描述动态润湿现象.这些结果表明,采用合理的接触角模型对动态湿润现象的模拟至关重要.
(a) t=0.004 s (b) t=0.010 s
(c) t=0.015 s (d) t=0.030 s
图14 不同时刻下的液滴形态对比图
Fig. 14 Comparison of droplet shapes between different contact angle models at different moments
本文基于相场方法和动态接触角模型,建立了能够模拟动态润湿行为的两相流数值模型,解决了动态湿润现象中的接触线附近的应力奇异性问题,并实现了接触角的准确描述.基于该模型,本文利用OpenFOAM开源平台进行了程序实现,通过对平直壁面上的液滴撞击模拟,并对比了理论值,验证了模型与程序的有效性和正确性.随后,基于该模型,对实验中液滴的动态润湿行为进行了模拟,模拟结果与实验数据取得了较高的吻合度.最后,通过不同的接触角模型的模拟结果对比,发现不同模型间的结果具有较大差异.在静态接触角边界条件下,模拟出现了液滴分离现象,与实验存在显著的差异;与准静态接触角模型的对比中,差异主要在液滴回缩过程中的速度,且回缩过程中较小的接触角能减少液滴达到稳定状态的时间,这是由于较小的接触角对接触线的移动具有更大的牵制作用.以上结果说明了准确的动态接触角模型对动态润湿现象的准确模拟具有至关重要的意义.
[1] FENG L, LI S, LI Y, et al. Super-hydrophobic surfaces: from natural to artificial[J]. Advanced Materials, 2002, 14(24): 1857-1860.
[2] YARIN A L. Drop impact dynamics: splashing, spreading, receding, bouncing...[J]. Annual Review of Fluid Mechanics, 2006, 38(1): 159-192.
[3] HUE P L. Progress and trends in ink-jet printing technology[J]. Journal of Imaging Science and Technology, 1998, 42: 49-62.
[4] YOUNG T. An essay on the cohesion of fluids[J]. Philosophical Transactions of the Royal Society of London, 1805, 95: 65-87.
[5] HUH C, SCRIVEN L E. Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[J]. Journal of Colloid and Interface Science, 1971, 35(1): 85-101.
[6] DUSSAN E B. On the spreading of liquids on solid surfaces: static and dynamic contact lines[J]. Annual Review of Fluid Mechanics, 1979, 11(1): 371-400.
[7] DE GENNES P G. Wetting: statics and dynamics[J]. Reviews of Modern Physics, 1985, 57(3): 827-863.
[8] JACQMIN D. Contact-line dynamics of a diffuse fluid interface[J]. Journal of Fluid Mechanics, 2000, 402: 57-88.
[9] 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(1): 169-194.
[10] VOINOV O V. Hydrodynamics of wetting[J]. Fluid Dynamics, 1977, 11(5): 714-721.
[11] HOCKING L M. The spreading of a thin drop by gravity and capillarity[J]. The Quarterly Journal of Mechanics and Applied Mathematics, 1983, 36(1): 55-69.
[12] TANNER L H. The spreading of silicone oil drops on horizontal surfaces[J]. Journal of Physics D: Applied Physics, 1979, 12(9): 1473-1484.
[13] UNVERDI S O, TRYGGVASON G. A front-tracking method for viscous, incompressible, multi-fluid flows[J]. Journal of Computational Physics, 1992, 100(1): 25-37.
[14] 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.
[15] OSHER S, SETHIAN J A. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[J]. Journal of Computational Physics, 1988, 79(1): 12-49.
[16] JACQMIN D. Calculation of two-phase Navier-Stokes flows using phase-field modeling[J]. Journal of Computational Physics, 1999, 155(1): 96-127.
[17] DING H, SPELT P D M. Inertial effects in droplet spreading: a comparison between diffuse-interface and level-set simulations[J]. Journal of Fluid Mechanics, 2007, 576: 287. DOI: 10.1017/S0022112007004910.
[18] QIAO L, ZENG Z, XIE H, et al. Modeling thermocapillary migration of interfacial droplets by a hybrid lattice Boltzmann finite difference scheme[J]. Applied Thermal Engineering, 2018, 131: 910-919.
[19] 周平, 曾忠, 乔龙. 假塑性流体液滴撞击壁面上的铺展的格子Boltzmann模拟[J]. 重庆大学学报, 2018, 41(12): 1-9.(ZHOU Ping, ZENG Zhong, QIAO Long. Simulation of shear-thinning droplets impact on solid surfaces by using lattice Boltzmann method[J]. Journal of Chongqing University, 2018, 41(12): 1-9.(in Chinese))
[20] YOKOI K, VADILLO D, HINCH J, et al. Numerical studies of the influence of the dynamic contact angle on a droplet impacting on a dry surface[J]. Physics of Fluids, 2009, 21(7): 072102. DOI: 10.1063/1.3158468.
[21] CAHN J W, HILLIARD J E. Free energy of a nonuniform system Ⅲ: nucleation in a two-component incompressible fluid[J]. The Journal of Chemical Physics, 1959, 31(3): 688-699.
[22] CHELLA R, VIALS J. Mixing of a two-phase fluid by cavity flow[J]. Physical Review E, 1996, 53(4): 3832-3840.
[23] DING H, SPELT P D. Wetting condition in diffuse interface simulations of contact line motion[J]. Physical Review E, 2007, 75(4): 046708. DOI: 10.1103/PhysRevE.75.046708.
[24] PASANDIDEH-FARD M, AZIZ S D, CHANDRA S, et al. Cooling effectiveness of a water drop impinging on a hot surface[J]. International Journal of Heat and Fluid Flow, 2001, 22(2): 201-210.
[25] OLSSON E, KREISS G, ZAHEDI S. A conservative level set method for two phase flow[J]. Journal of Computational Physics, 2007, 225: 785-807.
[26] JASAK H. Error analysis and estimation for finite volume method with applications to fluid flow[D]. PhD Thesis. London: Imperial College London, 1996.
[27] DUPONT J B, LEGENDRE D. Numerical simulation of static and sliding drop with contact angle hysteresis[J]. Journal of Computational Physics, 2010, 229(7): 2453-2478.
曾忠(1968—),教授,博士,博士生导师(通讯作者. E-mail: zzeng@cqu.edu.cn).