楔形体上复合液滴润湿铺展行为的格子Boltzmann 方法研究*

张晓林 黄军杰2)†

1) (重庆大学航空航天学院,重庆 400044)

2) (重庆大学,非均质材料力学重庆市重点实验室,重庆 400044)

固壁上液滴的润湿铺展行为是自然界中普遍存在的现象,针对楔形体上的复合液滴,采用基于相场理论的格子Boltzmann 方法对其润湿铺展行为进行探究.通过理论分析和数值模拟,发现液滴润湿面积随接触角、楔形体顶角的减小而增大,液滴也越容易分裂.处于理论分裂临界状态附近的液滴,在一定密度比、黏度比条件下将沿楔形体壁面分裂.基于模拟结果生成以密度比、黏度比为坐标的液滴分裂状态相图,比较发现相同条件下初始状态为平衡态的复合液滴更不易发生分裂.另外模拟还表明非对称界面张力及非对称运动黏度比也是影响液滴分裂结果的重要因素.

1 引言

由于固液、液液和气液分子之间的不同相互作用,固壁上的流动界面存在丰富独特的自然现象,在固体表面的众多性质中,浸润性是一个非常重要的特性.探究固壁上流体的润湿和铺展,能促进人们对固液相互作用的理解,这对于相关基础科学研究和前沿技术应用的发展都具有特殊意义.尤其是随着仿生材料制造技术的进步,已经实现了独特的液滴湿润和铺展行为,推动着自清洁表面[1]、涂装[2,3]、微液滴定向操作[4]和油水分离[5]等相关领域的快速发展.

液滴在固体壁面上的润湿铺展行为研究,主要涉及两类表面: 一类是光滑表面,如一般的平面和曲面;另一类是含奇异点的非光滑表面,如简单的楔形/矩形槽面、台面和复杂的微结构表面.对于一般光滑平面来说,对应的就是杨氏方程[6]所描述的理想情况,在忽略重力时液滴平衡呈球帽形状.Sui 等[7]从理论和实验方面分析流体与球形固体之间的润湿结果,发现接触角与杨氏接触角一致,液滴或气泡体积无穷大时,存在极限润湿位置.Li 等[8]研究了圆锥状基底的曲率比对其上液滴自发定向运动的影响,当曲率比超过临界值时,液滴就会停滞.Han 等[9]研究了二维液滴在两侧壁接触角不同的圆角V 形槽上的润湿行为,发现液滴的形态除与槽体开角、接触角相关外,还依赖于液滴体积和圆角半径.Wang 等[10]采用格子Boltzmann 方法(LBM)对两平行圆柱纤维上液滴的形态进行了模拟,确定了不同纤维间距和液滴体积下的完整形态图,存在三种可能的平衡构型: 桶形液滴、液滴桥和液柱,并提出了液滴形态转变的能量和力的解释.对于非光滑的表面,Herminghaus 等[11]讨论了液滴在楔形、矩形槽上的润湿,发现同一结构表面可以存在多种液滴形态.Chang 等[12]从理论和实验两方面验证了锥台表面角点对液滴润湿铺展所起的抑制作用.Zhou 等[13]探究了液滴在狭窄矩形台面上的润湿铺展.Ma 等[14]则考虑了液滴在圆形、三角形以及方形柱上的钉扎现象.对于其他微结构表面上液滴的润湿铺展研究[15−17],不再详述.

上述均是固壁上单液滴的润湿铺展研究,而复合液滴的润湿铺展属于三相流问题,相较于单液滴的二元流情形存在更多流动界面,同时流体和固体之间相互作用也更加复杂,因此相关数值研究更具挑战性.Said 等[18]对平壁面上复合液滴的润湿铺展行为进行了实验和数值探究.Weyer 等[19]研究了圆柱纤维上油水复合液滴的形态及液滴分离的相关影响因素.Zhang 等[20]分析了平壁面和毛细管中复合液滴的平衡构型.He 等[21]和Li 等[22]模拟了圆形界面上复合液滴的润湿铺展.Huang[23]研究了圆形和椭圆形界面上液滴的润湿平衡形态,并讨论了圆柱纤维上两液滴的吞噬包裹过程.由于液滴的润湿铺展行为与固壁形态存在很强关联,而目前还没有楔形体上复合液滴润湿铺展行为的相关研究,因此我们将就其开展探究,以提高对固壁上液滴微操作(如复合微液滴的分裂)的认知.

2 数值方法及理论模型

LBM 在诞生的30 多年里已经取得了长足的发展,具有许多优点,诸如算法的简单性、并行性,同时能够处理复杂的边界条件[24].由于LBM 的介观特性,在处理多相流问题时展现出了相当的优势.本文采用混合LB-有限差分方法[23]进行数值模拟研究,流场方程采用LBM 求解,基于Cahn-Hilliard 相场理论的界面演化方程,则采用二阶有限差分法进行空间离散和二阶Runge-Kutta 法进行时间推进.

2.1 数值方法

2.1.1 三元流体相场理论

对于三元流系统,自由能泛函定义为

其中系数ai为常量.σij为流体i和流体j之间的界面张力,γ1,γ2,γ3为与界面张力相关的参数:

系数ai,κi与参数γi、界面厚度W之间的关系为:ai=3γi/(4W) ,κi=3γiW/8 .化学势可表示为

其中,

可进一步得到

其中δi=γT/γi.Cahn 和Hilliard 将扩散通量近似为与化学势梯度成正比的量[25],得到Cahn-Hilliard方程(CHE).对于三元流体,有两个独立方程:

三相流体交界处的界面角φi(见补充材料A图S1)满足

在固壁上对于CHE,化学势µi应满足无通量边界条件

nS为固壁上指向流体的单位法向量;而润湿边界条件用于求固壁内与流体点相邻的虚拟点处的相序参数ϕ,以计算∇ϕ(x,t) 和∇2ϕ(x,t) .采用文献[26]发展的几何润湿边界条件,假设扩散界面附近ϕ的等值线相互平行,同时接触线处界面切线与接触角一致,随后建立与该切线平行的过虚拟点的特征线,并利用相序参数在界面法向满足双曲正切函数分布的特征来计算此点的ϕ值.该方法在曲面附近无需采用复杂的插值计算,也不用判断界面和网格的相对构型,实施起来较为方便.另外上述润湿边界条件使用的是文献[20]中提出的加权接触角:

其中θij(ij) 表示流体i-j界面与固壁之间的接触角,处于流体i一侧,具体见补充材料图S1.接触角θij和θji互为补角,即θij=π−θji.而且三个接触角不相互独立,由静力学平衡满足如下关系:

对于∇ϕ的计算,基于LBM 中D2Q9 模型使用的各向同性离散格式,在楔形体边界处采用二阶精度:

在离开楔形体边界的地方采用四阶精度[27]:

而∇2ϕ的计算在全域均采用二阶精度:

其中δx为网格尺寸,δt为时间步长,格子速度为c=δx/δt,el为格子速度矢量,wl为权系数.

2.1.2 LBM 求解Navier-Stokes 方程

使用单松弛碰撞模型时的格子Boltzmann 方程为[28]

通过相序参数可以计算流体密度和运动黏度:

其中ρi和νi分别为流体i的密度和运动黏度(相应的动力黏度ηi=ρiνi).宏观物理量的计算为

通过Chapman-Enskog 分析,(13)式在宏观尺度近似下列方程:

对于任意形状的固壁边界,采用插值回弹方法[29]求解未知的分布函数fl.为提高稳定性,模拟采用多松弛碰撞模型[30].

2.2 物理模型

本文主要研究楔形体上二维复合液滴(实际上对应第三维度无限长的复合液柱)的润湿铺展问题.在液滴尺寸较小、表面张力相对较大的情况下可以忽略重力的影响.物理模型如图1所示,计算域为Lx×Ly的矩形,楔形体顶角为χ,两侧壁长b.初始状态的复合液滴为Janus 状液滴,由关于过楔形体顶点的竖直线对称的两个半圆形液滴组成,半径为r0,其中左侧液滴为流体1 (Fluid1),右侧液滴为流体2 (Fluid2),周围环境为流体3 (Fluid3).计算域四周均采用周期边界条件.对称参数条件下σ12:σ13:σ23=1:1:1,接触角θ13=θ23(根据(9)式,θ12=90◦),流体1 与流体2 的密度比=1,运动黏度比=1,流体1 与流体3 的密度比为,运动黏度比为,动力黏度比为.文中取ρr=ρ1,νr=ν1,Lr=r0,Ur=σ12/η1,则有Tr=r0η1/σ12,Re=,模拟中Re=100.部分非对称情况的参数见4.4 节.

图1 物理模型图示Fig.1.Physical model illustration.

3 网格密度和计算域大小检验

模拟中将Lr用NL个网格离散,Tr用Nt个时间步离散,可得网格尺寸δx=Lr/NL、时间步长δt=Tr/Nt.相场模拟中有无量纲参数Cn=W/Lr=(W/δx)/NL,在给定NL的条件下,应保证界面厚度W足够小,但W过小时又无法精确捕捉界面,同时考虑到计算资源的消耗,取W/δx=4 .为减小棋盘效应的影响[31],计算域采用奇数网格离散.为计算稳定,取P e=2×104.

3.1 网格密度检验

网格数量的选择需要考虑计算资源消耗和计算精度两方面因素,我们在接触角θ13=105◦,楔形体顶角χ=90◦的情况下,考察网格密度对液滴润湿铺展模拟结果的影响.液滴润湿平衡形态的理论解推导见补充材料A.对于静态问题,计算域大小不影响液滴的最终平衡形态,选用Lx×Ly=5×5的计算域,网格密度为NL=30.2,40.2,60.2,80.2,100.2,相应的Nt=302,402,602,802,2004.不同网格密度下液滴平衡形态的模拟结果见图2(对称性图中只给出一侧液滴的情形,后同),可见除三相点和接触线位置附近的差别稍大外,界面其余位置的结果基本一致.我们还监测了接触线位置随时间的变化过程,由于初始时刻接触线处于楔形体顶点处,因此以接触线到楔形体顶点的距离dCL来表征接触线的位置,相应的计算结果如图3(a)所示.在液滴未分裂时,dCL还反映了其在楔形体表面的润湿面积.相较NL=30.2,40.2,在NL=60.2,80.2,100.2 时,接触线位置随时间的变化比较接近,而在趋近平衡时,5 种网格密度下液滴的接触线位置基本一致,综合考虑,在静力学计算中采用NL=40.2.对于动力学问题,针对Lx×Ly=5×5的计算域考察了=10,=10的液滴的动态润湿铺展过程,此时还采用了更大的网格密度NL=120.2,相应的各种网格分别对应Nt=302,402,1204,1604,2004,2404.液滴接触线位置的演化如图3(b)所示,可以看出在NL=100.2,120.2时的计算结果很接近,总体趋于收敛.但采用较大的网格密度对计算资源消耗巨大,同时接触线位置演化的误差主要存在于部分时刻(dCL为极值附近),且NL=60.2,120.2 时的相对误差约为2.1%,故后面动力学问题中的计算采用NL=60.2.

图2 不同网格密度下液滴的平衡形态Fig.2.Equilibrium morphology of the droplet at different mesh densities.

图3 不同网格密度下液滴接触线位置的演化 (a) ==1;(b) ==10Fig.3.Evolution of the contact line position of droplet under the different mesh densities: (a) ==1 ;(b) ==10 .

3.2 计算域大小检验

对于液滴分裂动力学问题,计算域大小将会影响液滴分裂的模拟结果,因此我们关注=10和=10的液滴在χ=90◦和θ13=80◦的楔形体上的润湿铺展,考察了5 种计算域大小(Lx×Ly=5×5,7.5×5.0,10×5,12.5×5.0,10×10,对应的网格数量分别为Nx×Ny=301×301,451×301,601×301,751×301,601×601)对模拟结果的影响.图4 为上述5 种计算域中模拟得到的液滴最终平衡形态,比较发现计算域为Lx×Ly=7.5×5.0,10×5,12.5×5.0,10×10 时的液滴形态一致,且较Lx×Ly=5×5 的结果更准确.图5 为上述不同计算域尺寸下液滴润湿铺展过程中接触线位置随时间演化的曲线,可见对于液滴动态演化过程也有类似的结果.因此取计算域x轴方向尺寸Lx≥7.5,同时y轴方向尺寸Ly≥5 对动力学模拟是合适的,本文采用Lx×Ly=10×5 的计算域.

图4 不同计算域 Lx×Ly 下液滴的平衡形态Fig.4.Equilibrium morphology of droplets in different computational domains of Lx×Ly .

图5 不同计算域Lx×Ly 下液滴接触线位置的演化Fig.5.Evolution of the contact line position of the droplet in different computational domains of Lx×Ly .

4 结果与分析

4.1 楔形体顶角、壁面润湿性对液滴平衡形态的影响

本节关注液滴润湿的平衡态问题,主要讨论楔形体顶角及壁面润湿性对未分裂复合液滴平衡形态的影响.因流体物性对本节考虑的液滴最终平衡形态无影响,为简单起见,模拟中选用的流体密度比、运动黏度比均为1∶1∶1 (即=1 和=1),计算域大小为Lx×Ly=5×5,其他相关参数为NL=40.2,Nt=402,C n=0.1 .首先探究楔形体顶角一定时,不同接触角下液滴平衡形态的差异.图6 为楔形体顶角χ=90◦时随接触角变化的液滴平衡形态图,能看到此时液滴的润湿特性与平壁面类似,即接触角越大,壁面越疏水,液滴润湿面积越小,三相点到楔形体顶点处的距离就越大.在此给出各种接触角下液滴三相点到楔形体顶点处的距离hnum,并与对应的理论值htheo(见补充材料A)进行比较,发现接触角θ13接近临界值2π/3−χ/2和 2π/3+χ/2 时,计算误差稍大,而其余情况则比较吻合(图7).

图6 不同接触角下的液滴平衡形态Fig.6.Equilibrium morphology of droplets at different contact angles.

图7 h 随接触角的变化曲线Fig.7.Variation of h with the contact angle.

需要注意的是,接触角为理论临界分裂值时(2π/3−χ/2=75◦),得到hnum>0,此时液滴并没有分裂.我们还讨论了楔形体顶角对液滴平衡形态的影响,图8 为θ13=90◦时三种不同顶角 (χ=75◦,90◦,105◦) 楔形体上液滴的平衡形态,对应的液滴润湿面积s分别为1.82,1.72,1.66,可知在壁面润湿特性一定时,楔形体顶角越小,被润湿的固壁面积就越大.根据理论分析我们在各种楔形体顶角、接触角下给出了液滴的临界分裂界线(图9),当液滴处于临界分裂界线左下方时对应分裂状态,处于临界分裂界线右上方时对应不分裂状态.(90◦,80◦)表示楔形体顶角χ=90◦,接触角θ13=80◦,此时理论预测液滴未分裂.分析表明接触角、楔形体顶角越小,液滴越容易处于分裂状态.

图8 不同顶角楔形体上液滴的平衡形态Fig.8.Equilibrium morphology of droplets on wedges with different vertex angles.

图9 液滴分裂/不分裂理论临界界线Fig.9.Theoretical critical boundary of droplet splitting/non-splitting state.

4.2 密度比、黏度比对液滴润湿分裂行为的影响

从4.1 节的分析可知,接触角和楔形体顶角是影响液滴润湿行为的两个重要因素,而对于动态润湿铺展过程,当接触角处于理论临界分裂值附近时,其动力学结果可能不同,故本节探究楔形体顶角χ=90◦,接触角θ13=80◦下密度比、黏度比对复合液滴润湿分裂行为的影响.计算域大小为Lx×Ly=10×5,取NL=60.2,Nt=1204,则C n=0.07,Nx×Ny=601×301.图10 给出两种典型的液滴不分裂和分裂的润湿铺展过程,分别对应=50,=1和=50 ,=5 .从图10 可以看到在演化的初期,液滴界面发生了扭曲变形,当三相点接近楔形体顶点时,在=1 的情况下液滴发生了回缩,最终不能分裂,而=5 时的液滴则发生了分裂.从速度场分布图10 中可知,在演化过程中液滴界面附近产生了涡,随着演化趋于平衡,涡逐渐远离液滴,并随黏性耗散而衰弱,直至消失.考虑液滴分裂与不分裂两种结果,在图11(a)建立了不同密度比、运动黏度比下楔形体上复合液滴润湿铺展的分裂状态相图,发现在运动黏度比≥3时,考虑的大部分密度比下的液滴发生了分裂,只是在小密度比=1 时,液滴没有分裂,分裂发生在运动黏度比稍大的情况下.可见当较小时,液滴惯性较小,不利于实现液滴分裂,同时液滴密度一定时,在较大黏度比的情况下,环境流体黏度较小,对液滴润湿铺展的阻碍作用就越小,则液滴更容易分裂.图11(a)稍做调整便得到了密度比、动力黏度比条件下的液滴分裂状态相图,如图11(b)所示.

图10 初始状态为非平衡态Janus 状液滴的润湿铺展过程及速度场分布 (a) =50 ,=1 ;(b) =50 ,=5Fig.10.Wetting and spreading process and velocity field distribution of Janus-like droplet with non-equilibrium initial state:(a) =50,=1 ;(b) =50 ,=5 .

图11 初始状态为非平衡态Janus 状液滴的润湿分裂状态相图 (a) rv13 -rρ13;(b) rη13 -rρ13Fig.11.Split/Non-split phase diagram of Janus-like droplets with non-equilibrium initial state: (a) rv13 vs. rρ13;(b) rη13 vs. rρ13.

在上述探究参数范围内,我们还考虑液滴润湿演化过程中的能量变化.图12 为不同密度比液滴在润湿铺展过程中所达到的系统动能最大值Ek,max随运动黏度比的变化曲线(小圆圈代表分裂,“+”符号代表未分裂),其中插图为演化过程中系统动能的变化曲线,其表达式为

图12 初始状态为非平衡态时系统动能最大值随运动黏度比的变化曲线Fig.12.Variation of the maximum kinetic energy of the system with the kinematic viscosity ratio under the nonequilibrium initial state.

4.3 液滴初始形态对其润湿分裂行为的影响

前面讨论的是初始条件为非平衡态Janus 状液滴的情形,本节将探究初始形态为平衡态的复合液滴在楔形体上的润湿铺展行为.平衡态的复合液滴由非平衡态Janus 状液滴演化而来,其半径为r=0.7884r0,具体推导见补充材料B.初始时刻平衡态复合液滴的下三相点位于楔形体顶点处,其余参数与4.2 节一致,我们仍关注该情况下密度比、黏度比对液滴润湿分裂行为的影响.图13 为此初始条件下液滴分裂、不分裂的两种典型过程,分别对应.与图10 比较,发现=5时初始形态为非平衡的液滴发生了分裂,而初始形态为平衡态的液滴却没有分裂,分裂发生在=10 .以密度比、黏度比为坐标轴生成的液滴分裂状态相图如图14 所示,比较图11(a)和图14(a)可知,此时液滴分裂发生在运动黏度比更大的情况;比较图11(b)和图14(b),发现复合液滴的分裂界线相较于初始状态为非平衡态的情况更加缓和,说明在相同密度比的条件下,初始状态为平衡态的复合液滴的分裂发生在动力黏度比更大的情况.总之,对于初始形态为平衡态的液滴,其界面能较低,系统总势能较小,在相同的条件下更不容易突破能量壁垒而发生分裂.要使液滴分裂,应该提高润湿铺展前期的惯性作用,可以通过减小环境流体(或增大液滴)的密度或黏度来实现.

图13 初始状态为平衡态复合液滴的润湿铺展过程及速度场分布 (a) rρ13=50 ,rν13=5 ;(b) rρ13=50 ,rν13=10Fig.13.Wetting and spreading process and velocity field distribution of compound droplet with equilibrium initial state:(a) r ρ13=50,r ν13=5 ;(b) r ρ13=50 ,r ν13=10 .

图14 初始状态为平衡态复合液滴的润湿分裂状态相图 (a) rv13 -rρ13;(b) rη13 -rρ13Fig.14.Split/non-split phase diagram of compound droplets with equilibrium initial state: (a) rv13 vs. rρ13;(b) rη13 vs. rρ13.

4.4 非对称条件下复合液滴的润湿铺展行为

前面讨论的均是对称参数下复合液滴的润湿铺展行为,本节研究非对称参数下复合液滴的润湿铺展.在研究分裂问题时,根据4.3节的结果,将讨论的参数选取在分裂状态相图分界线附近.考虑界面张力不对称时的两种情形,分别对应σ12:σ13:σ23=1.000:1.115:0.816(ϕ1=135◦,ϕ2=105°,ϕ3=120◦)以及σ12:σ13:σ23=1.000:1.155:0.577 (ϕ1=150◦,ϕ2=90◦,ϕ3=120◦),平衡时两侧液滴具有形状上的非对称性,相关推导见补充材料B.在此我们采用稍小的接触角θ13=θ23=75◦,根据(9)式可得θ12分别为 85.564◦和 81.406◦.初始时刻平衡态复合液滴的下三相点位于楔形体顶点处,其余参数与4.2 节一致.在对称和非对称界面张力的条件下,我们观察了rρ13=50 和rν13=7 液滴的润湿铺展过程,其左右侧接触线位置随时间的演化如图15 所示(图例中LS 表示左侧液滴接触线,RS 表示右侧液滴接触线).在该条件下发生分裂的临界值为dCL=2.35,当dCL>2.35 时说明液滴发生了分裂.图15 还给出t=0,2,4,30 时刻液滴的形态,可以看到在σ12:σ13:σ23=1:1:1 和σ12:σ13:σ23=1.000:1.115:0.816两种情况下,液滴均发生了分裂,但相较于对称情形,非对称界面张力下液滴的分裂程度较低;在σ12:σ13:σ23=1.000:1.155:0.577时,液滴没有分裂.根据上述三种界面张力下的模拟结果,我们发现随某一界面角的减小,一侧液滴被另一侧液滴包裹的程度更大,相应地就越不容易发生分裂.

图15 左右侧液滴接触线位置的演化Fig.15.Evolution of the position of the left and right droplet contact lines.

图16 左右侧液滴接触线位置的演化Fig.16.Evolution of the position of the left and right droplet contact lines.

5 结论

本文对楔形体上复合液滴的润湿铺展行为进行了数值研究.对于固壁上液滴润湿的静态问题,理论分析表明接触角、楔形体顶角越小,液滴越容易分裂.在不发生分裂的情况下,楔形体顶角一定时,液滴的润湿特性与平壁面类似,即接触角越大,壁面越疏水,液滴润湿面积越小;壁面润湿特性一定时,楔形体顶角越大,被润湿的固壁面积就越小.对液滴润湿铺展的动态问题,当液滴处于理论临界分裂接触角附近时,在一定的密度比、黏度比条件下会沿楔形体壁面发生分裂.基于模拟结果生成以密度比、运动黏度比为坐标的液滴分裂状态相图,通过能量分析表明密度比、运动黏度比越大时,润湿铺展过程中的惯性效应越大,液滴就更容易分裂.对比液滴的分裂状态相图,发现在相同的条件下初始形态为平衡态的复合液滴更不易发生分裂.对于非对称情况,我们发现左右侧液滴运动黏度差的增加有利于液滴分裂,同时界面张力不对称导致的液滴被包裹程度越大,液滴越不易分裂.本研究只讨论了对称参数和部分非对称参数下楔形体上二维复合液滴的润湿铺展行为,对于三维问题,将在以后的工作中进一步探讨.