搜索

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

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

张晓林 黄军杰

引用本文:
Citation:

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

张晓林, 黄军杰

Study on wetting and spreading behaviors of compound droplets on wedge by lattice Boltzmann method

Zhang Xiao-Lin, Huang Jun-Jie
PDF
HTML
导出引用
  • 固壁上液滴的润湿铺展行为是自然界中普遍存在的现象, 针对楔形体上的复合液滴, 采用基于相场理论的格子Boltzmann方法对其润湿铺展行为进行探究. 通过理论分析和数值模拟, 发现液滴润湿面积随接触角、楔形体顶角的减小而增大, 液滴也越容易分裂. 处于理论分裂临界状态附近的液滴, 在一定密度比、黏度比条件下将沿楔形体壁面分裂. 基于模拟结果生成以密度比、黏度比为坐标的液滴分裂状态相图, 比较发现相同条件下初始状态为平衡态的复合液滴更不易发生分裂. 另外模拟还表明非对称界面张力及非对称运动黏度比也是影响液滴分裂结果的重要因素.
    The wetting and spreading of droplets on solid walls are commonly seen in nature. The study of such a phenomenon can deepen our understanding of solid-liquid interaction and promote the development of relevant cutting-edge technological applications. In this work, the lattice Boltzmann method based on phase field theory is used to investigate the wetting and spreading of a compound droplet on a wedge. This method combines the finite-difference solution of the Cahn-Hilliard equations for ternary fluids to capture the interface dynamics and the lattice Boltzmann method for the hydrodynamics of the flow. Symmetric compound droplets with equal interfacial tensions on a wedge are considered first. Through theoretical analysis and numerical simulation, it is found that the wetted area on the wedge increases with the decrease of the contact angle of the wedge surface and the wedge apex angle. Depending on these two factors, the droplet may or may not split on the wedge. We also find that the droplet near the critical state predicted not to split by static equilibrium analysis could split during the spreading along the wall of the wedge under certain density and viscosity ratios. Based on the simulation results, a phase diagram of the droplet splitting state is generated with the density ratio and viscosity ratio as the coordinates. As the density ratio and kinematic viscosity ratio increase, the inertia effect becomes more prominent in the wetting and spreading process and the droplet is more likely to split. By comparing the phase diagrams in different initial conditions, it is found that under the same conditions, the compound droplet with an equilibrium initial state is less likely to split than that with an unequilibrium initial state, which is possibly because the initial total energy of the former is relatively small. Our study also shows that the kinematic viscosity ratio between the left half and the right half droplet may affect the results of droplet splitting. The increase of such a viscosity difference is conducive to the splitting of the compound droplet. Besides, asymmetric compound droplets with unequal interfacial tensions are also simulated, and it is found that the greater the wrapping degree between the left half and right half, the more difficult it is to separate the compound droplet.
      通信作者: 黄军杰, jjhuang@cqu.edu.cn
    • 基金项目: 国家自然科学基金(批准号: 11972098)资助的课题.
      Corresponding author: Huang Jun-Jie, jjhuang@cqu.edu.cn
    • Funds: Project supported by the National Natural Science Foundation of China (Grant No. 11972098).
    [1]

    Latthe S S, Sutar R S, Kodag V S, et al. 2019 Prog. Org. Coat. 128 52Google Scholar

    [2]

    Woerthmann B M, Totzauer L, Briesen H 2022 Powder Technol. 404 117443Google Scholar

    [3]

    Eres M H, Schwartz L W, Roy R V 2000 Phys. Fluids 12 1278Google Scholar

    [4]

    Dai Q W, Huang W, Wang X L, Khonsari M M 2021 Tribol. Int. 154 106749Google Scholar

    [5]

    Yang Y, Li X J, Zheng X, Chen Z Y, Zhou Q F, Chen Y 2018 Adv. Mater. 30 1704912Google Scholar

    [6]

    Young T 1805 Philos. Trans. R. Soc. London 95 65

    [7]

    Sui T, Wang J D, Chen D R 2011 J. Colloid Interface Sci. 358 284Google Scholar

    [8]

    Li Y Q, Wu H A, Wang F C 2016 J. Phys. D Appl. Phys. 49 085304Google Scholar

    [9]

    Han Z Y, Duan L, Kang Q 2019 AIP Adv. 9 085203Google Scholar

    [10]

    Wang F, Schiller U D 2021 Soft Matter 17 5486Google Scholar

    [11]

    Herminghaus S, Brinkmann M, Seemann R 2008 Ann. Rev. Mater. Res. 38 101Google Scholar

    [12]

    Chang F M, Hong S J, Sheng Y J, Tsao H K 2010 J. Phys. Chem. C 114 1615Google Scholar

    [13]

    Zhou L M, Yang S M, Quan N N, et al. 2021 ACS Appl. Mater. Interfaces 13 55726Google Scholar

    [14]

    Ma B J, Shan L, Dogruoz B, Agonafer D 2019 Langmuir 35 12264Google Scholar

    [15]

    Courbin L, Bird J C, Reyssat M, Stone H A 2009 J. Phys. Condes. Matter 21 464127Google Scholar

    [16]

    Frank X, Perre P 2012 Phys. Fluids 24 042101Google Scholar

    [17]

    Lee Y, Matsushima N, Yada S, Nita S, Kodama T, Amberg G, Shiomi J 2019 Sci. Rep. 9 7787Google Scholar

    [18]

    Ben Said M, Selzer M, Nestler B, Braun D, Greiner C, Garcke H 2014 Langmuir 30 4033Google Scholar

    [19]

    Weyer F, Ben Said M, Hotzer J, Berghoff M, Dreesen L, Nestler B, Vandewalle N 2015 Langmuir 31 7799Google Scholar

    [20]

    Zhang C Y, Ding H, Gao P, Wu Y L 2016 J. Comput. Phys. 309 37Google Scholar

    [21]

    He Q, Li Y J, Huang W F, Hu Y, Wang Y M 2020 Phys. Rev. E 101 033307Google Scholar

    [22]

    Li S, Lu Y, Jiang F, Liu H H 2021 Phys. Rev. E 104 015310Google Scholar

    [23]

    Huang J J 2021 Phys. Fluids 33 072105Google Scholar

    [24]

    Chen S Y, Doolen G D 1998 Annu. Rev. Fluid Mech. 30 329Google Scholar

    [25]

    Jacqmin D 1999 J. Comput. Phys. 155 96Google Scholar

    [26]

    Huang J J, Wu J, Huang H B 2018 Eur. Phys. J. E 41 1Google Scholar

    [27]

    Liang H, Chai Z H, Shi B C, Guo Z L, Zhang T 2014 Phys. Rev. E 90 063311Google Scholar

    [28]

    Lee T 2009 Comput. Math. Appl. 58 987Google Scholar

    [29]

    Bouzidi M, Firdaouss M, Lallemand P 2001 Phys. Fluids 13 3452Google Scholar

    [30]

    Lallemand P, Luo L S 2000 Phys. Rev. E 61 6546Google Scholar

    [31]

    Guo Z L, Shi B C, Zheng C G 2011 Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 369 2283Google Scholar

    [32]

    Carlson A, Do-Quang M, Amberg G 2011 J. Fluid Mech. 682 213Google Scholar

  • 图 1  物理模型图示

    Fig. 1.  Physical model illustration.

    图 2  不同网格密度下液滴的平衡形态

    Fig. 2.  Equilibrium morphology of the droplet at different mesh densities.

    图 3  不同网格密度下液滴接触线位置的演化 (a) $ {r_{{\rho _{13}}}} = $$ {r_{{\nu _{13}}}} = 1 $; (b) $ {r_{{\rho _{13}}}} = {r_{{\nu _{13}}}} = 10 $

    Fig. 3.  Evolution of the contact line position of droplet under the different mesh densities: (a) $ {r_{{\rho _{13}}}} = {r_{{\nu _{13}}}} = 1 $;(b) $ {r_{{\rho _{13}}}} = $$ {r_{{\nu _{13}}}} = 10 $.

    图 4  不同计算域$ {L_x} \times {L_y} $下液滴的平衡形态

    Fig. 4.  Equilibrium morphology of droplets in different computational domains of $ {L_x} \times {L_y} $.

    图 5  不同计算域$ {L_x} \times {L_y} $下液滴接触线位置的演化

    Fig. 5.  Evolution of the contact line position of the droplet in different computational domains of $ {L_x} \times {L_y} $.

    图 6  不同接触角下的液滴平衡形态

    Fig. 6.  Equilibrium morphology of droplets at different contact angles.

    图 7  h随接触角的变化曲线

    Fig. 7.  Variation of h with the contact angle.

    图 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.

    图 10  初始状态为非平衡态Janus状液滴的润湿铺展过程及速度场分布 (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 1 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $

    Fig. 10.  Wetting and spreading process and velocity field distribution of Janus-like droplet with non-equilibrium initial state: (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 1 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $.

    图 11  初始状态为非平衡态Janus状液滴的润湿分裂状态相图 (a) rv13 - rρ13; (b) rη13 - rρ13

    Fig. 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  初始状态为非平衡态时系统动能最大值随运动黏度比的变化曲线

    Fig. 12.  Variation of the maximum kinetic energy of the system with the kinematic viscosity ratio under the non-equilibrium initial state.

    图 13  初始状态为平衡态复合液滴的润湿铺展过程及速度场分布 (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 10 $

    Fig. 13.  Wetting and spreading process and velocity field distribution of compound droplet with equilibrium initial state: (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 10 $.

    图 14  初始状态为平衡态复合液滴的润湿分裂状态相图 (a) rv13 - rρ13; (b) rη13 - rρ13

    Fig. 14.  Split/non-split phase diagram of compound droplets with equilibrium initial state: (a) rv13 vs. rρ13; (b) rη13 vs. rρ13

    图 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.

  • [1]

    Latthe S S, Sutar R S, Kodag V S, et al. 2019 Prog. Org. Coat. 128 52Google Scholar

    [2]

    Woerthmann B M, Totzauer L, Briesen H 2022 Powder Technol. 404 117443Google Scholar

    [3]

    Eres M H, Schwartz L W, Roy R V 2000 Phys. Fluids 12 1278Google Scholar

    [4]

    Dai Q W, Huang W, Wang X L, Khonsari M M 2021 Tribol. Int. 154 106749Google Scholar

    [5]

    Yang Y, Li X J, Zheng X, Chen Z Y, Zhou Q F, Chen Y 2018 Adv. Mater. 30 1704912Google Scholar

    [6]

    Young T 1805 Philos. Trans. R. Soc. London 95 65

    [7]

    Sui T, Wang J D, Chen D R 2011 J. Colloid Interface Sci. 358 284Google Scholar

    [8]

    Li Y Q, Wu H A, Wang F C 2016 J. Phys. D Appl. Phys. 49 085304Google Scholar

    [9]

    Han Z Y, Duan L, Kang Q 2019 AIP Adv. 9 085203Google Scholar

    [10]

    Wang F, Schiller U D 2021 Soft Matter 17 5486Google Scholar

    [11]

    Herminghaus S, Brinkmann M, Seemann R 2008 Ann. Rev. Mater. Res. 38 101Google Scholar

    [12]

    Chang F M, Hong S J, Sheng Y J, Tsao H K 2010 J. Phys. Chem. C 114 1615Google Scholar

    [13]

    Zhou L M, Yang S M, Quan N N, et al. 2021 ACS Appl. Mater. Interfaces 13 55726Google Scholar

    [14]

    Ma B J, Shan L, Dogruoz B, Agonafer D 2019 Langmuir 35 12264Google Scholar

    [15]

    Courbin L, Bird J C, Reyssat M, Stone H A 2009 J. Phys. Condes. Matter 21 464127Google Scholar

    [16]

    Frank X, Perre P 2012 Phys. Fluids 24 042101Google Scholar

    [17]

    Lee Y, Matsushima N, Yada S, Nita S, Kodama T, Amberg G, Shiomi J 2019 Sci. Rep. 9 7787Google Scholar

    [18]

    Ben Said M, Selzer M, Nestler B, Braun D, Greiner C, Garcke H 2014 Langmuir 30 4033Google Scholar

    [19]

    Weyer F, Ben Said M, Hotzer J, Berghoff M, Dreesen L, Nestler B, Vandewalle N 2015 Langmuir 31 7799Google Scholar

    [20]

    Zhang C Y, Ding H, Gao P, Wu Y L 2016 J. Comput. Phys. 309 37Google Scholar

    [21]

    He Q, Li Y J, Huang W F, Hu Y, Wang Y M 2020 Phys. Rev. E 101 033307Google Scholar

    [22]

    Li S, Lu Y, Jiang F, Liu H H 2021 Phys. Rev. E 104 015310Google Scholar

    [23]

    Huang J J 2021 Phys. Fluids 33 072105Google Scholar

    [24]

    Chen S Y, Doolen G D 1998 Annu. Rev. Fluid Mech. 30 329Google Scholar

    [25]

    Jacqmin D 1999 J. Comput. Phys. 155 96Google Scholar

    [26]

    Huang J J, Wu J, Huang H B 2018 Eur. Phys. J. E 41 1Google Scholar

    [27]

    Liang H, Chai Z H, Shi B C, Guo Z L, Zhang T 2014 Phys. Rev. E 90 063311Google Scholar

    [28]

    Lee T 2009 Comput. Math. Appl. 58 987Google Scholar

    [29]

    Bouzidi M, Firdaouss M, Lallemand P 2001 Phys. Fluids 13 3452Google Scholar

    [30]

    Lallemand P, Luo L S 2000 Phys. Rev. E 61 6546Google Scholar

    [31]

    Guo Z L, Shi B C, Zheng C G 2011 Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 369 2283Google Scholar

    [32]

    Carlson A, Do-Quang M, Amberg G 2011 J. Fluid Mech. 682 213Google Scholar

计量
  • 文章访问数:  251
  • PDF下载量:  8
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-07-22
  • 修回日期:  2022-10-22
  • 上网日期:  2022-10-27
  • 刊出日期:  2023-01-20

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

  • 1. 重庆大学航空航天学院, 重庆 400044
  • 2. 重庆大学, 非均质材料力学重庆市重点实验室, 重庆 400044
  • 通信作者: 黄军杰, jjhuang@cqu.edu.cn
    基金项目: 国家自然科学基金(批准号: 11972098)资助的课题.

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

English Abstract

    • 由于固液、液液和气液分子之间的不同相互作用, 固壁上的流动界面存在丰富独特的自然现象, 在固体表面的众多性质中, 浸润性是一个非常重要的特性. 探究固壁上流体的润湿和铺展, 能促进人们对固液相互作用的理解, 这对于相关基础科学研究和前沿技术应用的发展都具有特殊意义. 尤其是随着仿生材料制造技术的进步, 已经实现了独特的液滴湿润和铺展行为, 推动着自清洁表面[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]研究了圆形和椭圆形界面上液滴的润湿平衡形态, 并讨论了圆柱纤维上两液滴的吞噬包裹过程. 由于液滴的润湿铺展行为与固壁形态存在很强关联, 而目前还没有楔形体上复合液滴润湿铺展行为的相关研究, 因此我们将就其开展探究, 以提高对固壁上液滴微操作(如复合微液滴的分裂)的认知.

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

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

      $\begin{split} &{\mathcal{F}}\left({\phi}_{1},{\phi}_{2},{\phi}_{3},\nabla {\phi}_{1},\nabla {\phi}_{2},\nabla {\phi}_{3}\right)\\=\; &{{\displaystyle \int }}_{V}\left[\varPsi \left({\phi}_{1},{\phi}_{2},{\phi}_{3}\right)+\frac{1}{2}\sum _{i=1}^{3}{\kappa }_{i}{\left|\nabla {\phi}_{i}\right|}^{2}\right]\text{d}V,\end{split}$

      其中${\phi _i}\left( {i = 1, 2, 3} \right)$为相序参数, 变化范围为[–1, 1], 满足$\displaystyle \sum \nolimits_{i=1}^{3}{\phi }_{i}=-1 $, 故只有两个独立的相序参数. (1)式右侧积分中的第一项$ \varPsi \left( {{\phi _1}, {\phi _2}, {\phi _3}} \right) $为体自由能密度, 第二项为界面能密度, 本文中取

      $ \varPsi \left({\phi }_{1},{\phi }_{2},{\phi }_{3}\right)=\sum _{i=1}^{3}{a}_{i}{\left({\phi }_{i}^{2}-1\right)}^{2} \text{, } $

      其中系数$ {a_i} $为常量. $ {\sigma _{ij}} $为流体$ i $和流体$ j $之间的界面张力, γ1, γ2, γ3为与界面张力相关的参数:

      $\begin{split} &{\gamma _1} = \frac{1}{2}\left( {{\sigma _{12}} + {\sigma _{13}} - {\sigma _{23}}} \right), \\& {\gamma _2} = \frac{1}{2}\left( {{\sigma _{12}} + {\sigma _{23}} - {\sigma _{13}}} \right), \\& {\gamma _3} = \frac{1}{2}\left( {{\sigma _{13}} + {\sigma _{23}} - {\sigma _{12}}} \right).\end{split} $

      系数$ {a_i} $, $ {\kappa _i} $与参数$ {\gamma _i} $、界面厚度$ W $之间的关系为: ${a_i} = {{3{\gamma _i}}}/{{(4 W)}}$, $ {\kappa _i} = {{3{\gamma _i}W}}/{8}$. 化学势可表示为

      $ {\mu _i} = \frac{{\delta}\varPsi }{{\delta}{\phi _i}} + {\beta _\phi } \text{, } $

      其中,

      $\begin{split} & {\beta }_{\phi }=-\sum _{j=1}^{3}\frac{4{\gamma }_{\text{T}}}{W{\gamma }_{j}}\frac{\partial \varPsi }{\partial {\phi }_{j}},\\ & {\gamma _{\text{T}}} = \frac{{3{\gamma _1}{\gamma _2}{\gamma _3}}}{{{\gamma _1}{\gamma _2} + {\gamma _1}{\gamma _3} + {\gamma _2}{\gamma _3}}},\end{split} $

      可进一步得到

      $\begin{split} {\mu _i} =\; & 4{a_i}\left[ {\phi _i}\left( {\phi _i^2 - 1} \right) - {\delta _i}\left( {1 + {\phi _1}} \right)\right.\\ & \times\left.\left( {1 + {\phi _2}} \right)\left( {1 + {\phi _3}} \right) \right] - {\kappa _i}{\nabla ^2}{\phi _i} , \end{split}$

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

      $ \frac{{\partial {\phi _i}}}{{\partial t}} + {\boldsymbol{u}} \cdot \nabla {\phi _i} = \nabla \cdot \left( {{M_i}\nabla {\mu _i}} \right)\quad (i = 1,2) \text{, } $

      其中常数$ {M_i} = \dfrac{{{M_0}}}{{{\gamma _i}}} $是扩散率. $ {\mu _3} $可以通过关系式$ \displaystyle\sum\limits_{i = 1}^3 {\dfrac{{{\mu _i}}}{{{\gamma _i}}}} = 0 $得到. 取参考尺寸$ {L_{\text{r}}} $, 参考速度$ {U_{\text{r}}} $, 参考时间$ {T_{\text{r}}} = {{{L_{\text{r}}}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}} {{U_{\text{r}}}}}} \right. } {{U_{\text{r}}}}} $, 将(5)式无量纲化得到:

      $ \frac{{\partial {\phi _i}}}{{\partial {t^ * }}} + {{\boldsymbol{u}}^ * } \cdot {\nabla ^ * }{\phi _i} = \frac{1}{{Pe}}{\left( {{\nabla ^ * }} \right)^2}\mu _i^ * \quad (i = 1,2) \text{, } $

      其中$ Pe = {U_{\text{r}}}L_{\text{r}}^{\text{2}}/{M_0} $, $ \mu _i^ * = {{{L_{\text{r}}}{\mu _i}} \mathord{\left/ {\vphantom {{{L_{\text{r}}}{\mu _i}} {{\gamma _i}}}} \right. } {{\gamma _i}}} $.

      三相流体交界处的界面角$ {\varphi _i} $ (见补充材料A图S1)满足

      $ \sum\limits_{i = 1}^3 {{\varphi _i}} = 2{\text{π }} \text{, } \frac{{\sin {\varphi _3}}}{{{\sigma _{12}}}} = \frac{{\sin {\varphi _2}}}{{{\sigma _{13}}}} = \frac{{\sin {\varphi _1}}}{{{\sigma _{23}}}} . $

      在固壁上对于CHE, 化学势$ {\mu _i} $应满足无通量边界条件

      $ {{\boldsymbol{n}}_{\text{S}}} \cdot \nabla {\mu _i} = \frac{{\partial {\mu _i}}}{{\partial {{\boldsymbol{n}}_{\text{S}}}}} = 0,$

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

      $\begin{split} & {\theta _1} = \frac{{\left( {1 + {\phi _3}} \right){\theta _{13}} + \left( {1 + {\phi _2}} \right){\theta _{12}}}}{{\left( {1 + {\phi _3}} \right) + \left( {1 + {\phi _2}} \right)}} ,\\ & {\theta _2} = \frac{{\left( {1 + {\phi _3}} \right){\theta _{23}} + \left( {1 + {\phi _1}} \right){\theta _{21}}}}{{\left( {1 + {\phi _3}} \right) + \left( {1 + {\phi _1}} \right)}},\end{split} $

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

      $ \sin {\phi _2}\cos {\theta _{13}} - \sin {\phi _3}\cos {\theta _{12}} - \sin {\phi _1}\cos {\theta _{23}} = 0 . $

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

      $ {\left[ {\nabla \phi ({{\boldsymbol{x}}},t)} \right]^{(2)}} = \frac{3}{{c{\delta _x}}}\sum\limits_{l = 1}^8 {{w_l}{{{\boldsymbol{e}}}_l}\phi ({{\boldsymbol{x}}} + {{{\boldsymbol{e}}}_l}{\delta _t},t)} . $

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

      $\begin{split} &{\left[ {\nabla \phi ({{\boldsymbol{x}}},t)} \right]^{(4)}} \\=\;& \frac{1}{{2c{\delta _x}}}\sum\limits_{l = 1}^8 {{w_l}{{{\boldsymbol{e}}}_l}\left[ {8\phi ({{\boldsymbol{x}}} + {{{\boldsymbol{e}}}_l}{\delta _t},t) - \phi ({{\boldsymbol{x}}} + 2{{{\boldsymbol{e}}}_l}{\delta _t},t)} \right]} .\end{split}$

      ${\nabla ^2}\phi $的计算在全域均采用二阶精度:

      $\begin{split} & {\left[ {{\nabla ^2}\phi ({{\boldsymbol{x}}},t)} \right]^{(2)}} \\=\;& \frac{6}{{\delta _x^2}}\sum\limits_{l = 1}^8 {\left[ {{w_l}\phi ({{\boldsymbol{x}}}+{{{\boldsymbol{e}}}_l}{\delta _t},t) - (1 - {w_0})\phi ({{\boldsymbol{x}}},t)} \right]} ,\end{split} $

      其中${\delta _x}$为网格尺寸, ${\delta _t}$为时间步长, 格子速度为$ c = {\delta _x}/{\delta _t} $, ${{{\boldsymbol{e}}}_l}$为格子速度矢量, ${w_l}$为权系数.

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

      $\begin{split} & {f_l}\left( {{\boldsymbol{x}} + {{\boldsymbol{e}}_l}{\delta _t},t + {\delta _t}} \right) - {f_l}({{\boldsymbol{x}}},t) \\ =\; & - \frac{1}{\tau }\left( {{f_l} - f_{_l}^{{\text{eq}}}} \right) + \left( {1 - \frac{1}{{2\tau }}} \right)\left( {{{\boldsymbol{e}}_l} - {\boldsymbol{u}}} \right) \\ &\times \left[ {\nabla \rho c_{\text{s}}^{\text{2}}\left[ {{\varGamma _l} - {\varGamma _l}(0)} \right] + {{\boldsymbol{F}}_{{\text{ST}}}}{\varGamma _l}} \right]{\delta _t},\end{split} $

      其中$ {c_{\text{s}}} $表示格子声速$\left( {c_{\text{s}}} = {c}/{{\sqrt 3 }} \right)$, $ \tau $为松弛时间, ${{\boldsymbol{F}}}_{\text{ST}}=\displaystyle\sum\nolimits_{i=1}^{3}{\mu }_{i}\nabla {\phi }_{i}$表征界面张力作用, $ {f_l} $$ f_{_l}^{{\text{eq}}} $分别为沿着格子速度矢量${{{\boldsymbol{e}}}_l}$的分布函数和平衡分布函数,

      $\begin{split} f_l^{{\text{eq}}} =\; & {w_l}\left[ p + \rho c_{\text{s}}^{\text{2}}\left( \frac{1}{{c_{\text{s}}^{\text{2}}}}{{e}_{l\alpha }}{u_\alpha }\right.\right.\\ &\left.\left. + \frac{1}{{2c_{\text{s}}^{\text{4}}}}\left( {{{e}_{l\alpha }}{{e}_{l\beta }} - c_{\text{s}}^{\text{2}}{\delta _{\alpha \beta }}} \right){u_\alpha }{u_\beta } \right) \right].\end{split} $

      ${\varGamma _l}$为动水压力,

      $ {\varGamma _l} = {w_l}\left[ {1 + \frac{1}{{c_{\text{s}}^{\text{2}}}}{{e}_{l\alpha }}{u_\alpha } + \frac{1}{{2c_{\text{s}}^{\text{4}}}}\left( {{{e}_{l\alpha }}{{e}_{l\beta }} - c_{\text{s}}^{\text{2}}{\delta _{\alpha \beta }}} \right){u_\alpha }{u_\beta }} \right] . $

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

      $ \rho ({\phi }_{i})=\sum _{i=1}^{3}\frac{{\phi }_{i}+1}{2}{\rho }_{i} ,~ \nu ({\phi }_{i})={\left[\sum _{i=1}^{3}\frac{{\phi }_{i}+1}{2}\frac{1}{{\nu }_{i}}\right]}^{-1}, $

      其中$ {\rho _i} $$ {\nu _i} $分别为流体$ i $的密度和运动黏度(相应的动力黏度$ {\eta _i} = {\rho _i}{\nu _i} $). 宏观物理量的计算为

      $\begin{split} & p=\sum _{l=0}^{8}{f}_{l}+\frac{1}{2}{\delta }_{t}\left({\boldsymbol{u}}\cdot \nabla \rho {c}_{\text{s}}^{\text{2}}\right) , \\ & \rho {\boldsymbol{u}}=\frac{1}{{c}_{\text{s}}^{\text{2}}}\sum _{l=1}^{8}{f}_{l}{{\boldsymbol{e}}}_{l}+\frac{1}{2}{\delta }_{t}{{\boldsymbol{F}}}_{\text{ST}} . \end{split} $

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

      $ \frac{{\partial p}}{{\partial t}} + \rho c_{\text{s}}^{\text{2}}\nabla \cdot {\boldsymbol{u}} = 0 \text{, } $

      $ \rho \left( {\frac{{\partial {\boldsymbol{u}}}}{{\partial t}} + {\boldsymbol{u}} \cdot \nabla {\boldsymbol{u}}} \right) = - \nabla p + {{\boldsymbol{F}}_{{\text{ST}}}} + \nabla \cdot {\boldsymbol{\varPi}} , $

      其中${\boldsymbol{\varPi}} = \rho \nu [ {\nabla {\boldsymbol{u}} + {{(\nabla {\boldsymbol{u}})}^{\text{T}}}} ]$为黏性应力张量. 松弛参数和运动黏度的关系为$\nu = c_{\text{s}}^{\text{2}}( {\tau - 0.5} ){\delta _t}$. 取参考尺寸$ {L_{\text{r}}} $, 参考速度$ {U_{\text{r}}} $, 参考时间$ {T_{\text{r}}} = {{{L_{\text{r}}}}}/{{{U_{\text{r}}}}} $, 参考密度$ {\rho _{\text{r}}} $以及参考运动黏度$ {\nu _{\text{r}}} $, 将(18)式无量纲化得:

      $ {\rho ^ * }\left( {\frac{{\partial {{\boldsymbol{u}}^ * }}}{{\partial {t^ * }}} + {{\boldsymbol{u}}^ * } \cdot {\nabla ^ * }{{\boldsymbol{u}}^ * }} \right) = - {\nabla ^ * }{p^ * } + {\boldsymbol{F}}_{{\text{ST}}}^ * + \frac{1}{{Re}}{\nabla ^ * } \cdot {{\boldsymbol{\varPi}} ^ * } , $

      其中,

      $ {\rho ^ * } = \frac{\rho }{{{\rho _{\text{r}}}}} \text{, }~~ {p^ * } = \frac{p}{{{\rho _{\text{r}}}U_{\text{r}}^{\text{2}}}} \text{, } $

      $ {{\boldsymbol{F}}}_{\text{ST}}^{\text{*}}=\frac{{{\boldsymbol{F}}}_{\text{ST}}}{{\rho }_{\text{r}}{U}_{\text{r}}^{\text{2}}/{L}_{\text{r}}}=\frac{1}{We}\sum _{i=1}^{3}{\gamma }_{i}^{\ast }{\mu }_{i}^{\ast }{\nabla }^{\ast }{\phi }_{i} , $

      $ We = \frac{{{\rho _{\text{r}}}U_{\text{r}}^{\text{2}}{L_{\text{r}}}}}{{{\sigma _{12}}}} \text{, } ~~\gamma _i^ * = \frac{{{\gamma _i}}}{{{\sigma _{12}}}} \text{, } $

      $ Re = \frac{{{\rho _{\text{r}}}{U_{\text{r}}}{L_{\text{r}}}}}{{{\eta _{\text{r}}}}} = \frac{{{U_{\text{r}}}{L_{\text{r}}}}}{{{\nu _{\text{r}}}}} \text{, } $

      $ {{\boldsymbol{\varPi}} ^ * } = {\rho ^ * }{\nu ^ * }\left[ {{\nabla ^ * }{{\boldsymbol{u}}^ * } + {{\left( {{\nabla ^ * }{{\boldsymbol{u}}^ * }} \right)}^{\text{T}}}} \right] ,$

      $ {\nu ^*} = {\nu }/{{{\nu _{\text{r}}}}}. $

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

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

      图  1  物理模型图示

      Figure 1.  Physical model illustration.

    • 模拟中将$ {L_{\text{r}}} $$ {N_{\text{L}}} $个网格离散, $ {T_{\text{r}}} $$ {N_{\text{t}}} $个时间步离散, 可得网格尺寸${\delta _x} = {{{L_{\text{r}}}}}/{{{N_{\text{L}}}}}$、时间步长${\delta _t} = {{{T_{\text{r}}}}}/{{{N_{\text{t}}}}}$. 相场模拟中有无量纲参数$Cn = {W}/{{{L_{\text{r}}}}} = ({{W/{\delta _x}}})/{{{N_{\text{L}}}}}$, 在给定$ {N_{\text{L}}} $的条件下, 应保证界面厚度$ W $足够小, 但$ W $过小时又无法精确捕捉界面, 同时考虑到计算资源的消耗, 取$W/{\delta _x} = {\text{4}}$. 为减小棋盘效应的影响[31], 计算域采用奇数网格离散. 为计算稳定, 取$Pe = 2 \times {10^4}$.

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

      图  2  不同网格密度下液滴的平衡形态

      Figure 2.  Equilibrium morphology of the droplet at different mesh densities.

      图  3  不同网格密度下液滴接触线位置的演化 (a) $ {r_{{\rho _{13}}}} = $$ {r_{{\nu _{13}}}} = 1 $; (b) $ {r_{{\rho _{13}}}} = {r_{{\nu _{13}}}} = 10 $

      Figure 3.  Evolution of the contact line position of droplet under the different mesh densities: (a) $ {r_{{\rho _{13}}}} = {r_{{\nu _{13}}}} = 1 $;(b) $ {r_{{\rho _{13}}}} = $$ {r_{{\nu _{13}}}} = 10 $.

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

      图  4  不同计算域$ {L_x} \times {L_y} $下液滴的平衡形态

      Figure 4.  Equilibrium morphology of droplets in different computational domains of $ {L_x} \times {L_y} $.

      图  5  不同计算域$ {L_x} \times {L_y} $下液滴接触线位置的演化

      Figure 5.  Evolution of the contact line position of the droplet in different computational domains of $ {L_x} \times {L_y} $.

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

      图  6  不同接触角下的液滴平衡形态

      Figure 6.  Equilibrium morphology of droplets at different contact angles.

      图  7  h随接触角的变化曲线

      Figure 7.  Variation of h with the contact angle.

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

      图  8  不同顶角楔形体上液滴的平衡形态

      Figure 8.  Equilibrium morphology of droplets on wedges with different vertex angles.

      图  9  液滴分裂/不分裂理论临界界线

      Figure 9.  Theoretical critical boundary of droplet splitting/non-splitting state.

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

      图  10  初始状态为非平衡态Janus状液滴的润湿铺展过程及速度场分布 (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 1 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $

      Figure 10.  Wetting and spreading process and velocity field distribution of Janus-like droplet with non-equilibrium initial state: (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 1 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $.

      图  11  初始状态为非平衡态Janus状液滴的润湿分裂状态相图 (a) rv13 - rρ13; (b) rη13 - rρ13

      Figure 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为不同密度比液滴在润湿铺展过程中所达到的系统动能最大值${E_{{\rm{k}}, \max }}$随运动黏度比的变化曲线(小圆圈代表分裂, “+”符号代表未分裂), 其中插图为演化过程中系统动能的变化曲线, 其表达式为

      图  12  初始状态为非平衡态时系统动能最大值随运动黏度比的变化曲线

      Figure 12.  Variation of the maximum kinetic energy of the system with the kinematic viscosity ratio under the non-equilibrium initial state.

      $ {E_{\rm{k}}} = \iint {\frac{1}{2}\rho ({u^2} + {v^2})}{\text{d}}x{\text{d}}y \text{, } $

      单位为$ {\rho _1}U_{\text{r}}^{\text{2}}L_{\text{r}}^{\text{2}} $. 从插图可知在润湿铺展前期, 动能增加, 惯性逐渐起主导作用, 之后黏性耗散发挥主要作用, 动能逐渐减小[32]. 对于密度比较大的情况($ {r_{{\rho _{13}}}} \geqslant 10 $), 环境流体密度较小, 其动能可以忽略, 故系统动能可近似看作复合液滴的动能. 通过比较, 发现液滴分裂均发生在动能峰值较大的情况, 说明密度比、运动黏度比越大时, 润湿铺展过程中的惯性效应越大, 液滴就更容易分裂.

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

      图  13  初始状态为平衡态复合液滴的润湿铺展过程及速度场分布 (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 10 $

      Figure 13.  Wetting and spreading process and velocity field distribution of compound droplet with equilibrium initial state: (a) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 5 $; (b) $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 10 $.

      图  14  初始状态为平衡态复合液滴的润湿分裂状态相图 (a) rv13 - rρ13; (b) rη13 - rρ13

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

      图  15  左右侧液滴接触线位置的演化

      Figure 15.  Evolution of the position of the left and right droplet contact lines.

      另外我们还关注了运动黏度比不对称的情形, 即$ {r_{{\nu _{12}}}} \ne 1 $. 初始时刻液滴处于平衡形态, $ {\theta _{13}} = {\theta _{23}} = 80^\circ $, 其余参数与4.2节保持一致. 图16给出了$ {r_{{\rho _{12}}}} = 1 $, $ {r_{{\rho _{13}}}} = 50 $, $ {r_{{\nu _{13}}}} = 7 $, 而$ {r_{{\nu _{12}}}} $= 1.0, 1.4, 2.0时左右侧液滴接触线位置随时间的演化, 此时分裂发生的临界值为$ {d_{{\text{CL}}}} = 2.23 $. 同样图16也给出了$ t $ = 0, 3, 30时刻液滴的形态, 可见在$ {r_{{\nu _{12}}}} $ = 1的情况下液滴没有分裂, 平衡时又恢复了对称状态,而在$ {r_{{\nu _{12}}}} $= 1.4, 2.0的情况下, 液滴发生了分裂, 平衡时左侧壁面上的液滴位置高于右侧壁面上的液滴. 同时随$ {r_{{\nu _{12}}}} $增大, 液滴分裂的程度也越大(即分裂后的液滴离楔形体顶点的距离越远), 表明随左右侧液滴黏度差的增加, 液滴更容易分裂.

      图  16  左右侧液滴接触线位置的演化

      Figure 16.  Evolution of the position of the left and right droplet contact lines.

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

参考文献 (32)
补充材料:

目录

    /

    返回文章
    返回