然而，对于更加一般、不能归结为准静态问题的情形，如亲水区域的长度与宽度相近且均小于lc的情形，限域吸附过程中气-液界面和三相接触线经历复杂的动力学演化过程，式(1)、(4)理论上不再适用。Brasjen等[17]对方形亲水区域限域吸附液滴最大高度进行了初步的理论和实验研究，实验结果表明，边长小于毛细长度的正方形亲水区域限域吸附产生液滴最大高度与毛细数之间的关系与式(4)存在明显差别，而且理论计算结果与实验结果呈现的趋势并不吻合。因此，目前对于亲水区域的长度与宽度相近且均小于lc的情形，仍然未能建立足够准确的hmax理论模型，导致目前仍然缺乏对利用此类亲水区域生成的液滴微阵列中液滴最大高度的系统、精准调控方法。
本文对长度、宽度均小于毛细长度的亲水区域限域吸附液滴最大高度进行了理论建模和数值研究。首先通过相场动力学和润滑近似方法分别对亲/疏水图案化基底限域吸附液体过程进行了数值建模，然后基于测量限域吸附溶液干燥后固体残留形貌的方法对限域吸附液滴最大高度进行了实验研究，最后对比了理论和实验结果。结果表明，本文提出的基于相场动力学的模拟方法比基于润滑近似的方法能更准确地描述亲水区域限域吸附生成液滴最大高度随毛细数的变化规律，揭示了限域液体吸附过程中黏性力对形成液滴的作用减弱的现象。
1 相场动力学方法
1.1 域内方程
求解域为接触线附近的一个半圆柱形区域，如图 1所示。求解域内采用Cahn-Hilliard形式的相场模型模拟气液两相界面的动态演化过程。这一模型的基本思想是：互不相溶的气液两相的体相内仅有气体或者液体单相，而在两相的界面处存在一个两相混合的过渡层，称为扩散界面。引入相场变量ϕ来描述这一过渡现象，该变量在气相取值为-1，在液相取值为1，而在两相的过渡界面区域取-1和1之间的值，对应两相界面过渡区内气相的体积分数和液相的体积分数分别为(1-ϕ)/2和(1+ϕ)/2，ϕ为相场变量。由此可以定义该过渡区域内的物性：
其中：ρ和μ分别为过渡区域内一点处的物质密度和动力黏度；ρg和μg分别为气相的密度和动力黏度，而ρl和μl分别为液相的密度和动力黏度。一般将ϕ=0的等值面定义为气液两相的实际界面所在位置。
两相之间的扩散界面会因两相内的流动和化学势梯度导致的界面扩散通量而发生演化，这一演化过程由Cahn等[22]提出的方程描述：
其中：u为局部流动速度；G为化学势；Γ为迁移率调整参数，表征界面扩散通量与化学势梯度之间的比例关系；化学势和迁移率调整参数中的ε为扩散界面厚度，表征两相之间混合过渡区的厚度；χ为迁移率调整参数和扩散界面厚度平方之间的比例常数；λ为混合能量密度。λ与γ之间可以建立如下联系：
为利用式(7)求解ϕ的分布，进而给出气液相界面的演化过程，还需要给出求解域内的u，因此除求解式(7)外，还需要求解不可压缩Newton流体的Navier-Stokes方程和连续性方程[23]：
其中：p为压力，g为重力加速度矢量。
1.2 边界与初始条件
图 1中AEBH和CFDG面采用对称边界条件：
其中：n为边界外法向，I为二阶单位张量。
ABCD面采用无滑移润湿壁边界条件：
其中：ni为界面法向，由液相指向气相；t为平行于壁移动速度方向的界面切向；θ为壁的接触角；γ▽G表示界面扩散通量。
其中：β=α-φ，α为该处界面切线与壁之间的夹角，φ为界面切线与楔形区域内某一点和原点连线的夹角；p0为气相压力；c和d为常数。该边界上指定流体正应力，正应力数值由式(21)推导得出：
该边界相场的边界条件设置为固定界面交线条件，即认为气液界面(ϕ=0等值面)与该边界的交线始终位于EF直线上：
其中：tanh函数为双曲正切函数，Dwi为初始界面到边界某点的距离。
对于圆弧面ADEF，该边界为环境大气边界，流场边界设置为压力边界条件：
相场边界同样设置为固定界面交线条件：
由于求解域尺度较小，仅为接触线附近区域，因此重力的影响可以忽略不计，流场的初始条件采用零流速等压力条件，即：
相场初始条件为：
1) 指定面BHGC、面EFGH和圆弧面BCEF围成区域初始为液相：
2) 指定面AHGD、面EFGH和圆弧面ADEF围成区域为气相，相场变量取值为
1.3 几何与物性参数设置、网格划分与求解
相场动力学方法的数值求解域中，半圆柱体的直径AB长度设置为200 μm，半圆柱体长度(HG长度)设置为200 μm。由于AEBH、CFDG面均采用对称边界条件，因此模拟中亲水图案沿垂直提拉方向的间距等于半圆柱体长度的2倍，为400 μm，与实验中亲水图案沿垂直提拉速度方向的间距保持一致。相场动力学方法涉及气体、液体两相内流场和相界面的计算，因此需要指定两相的物性参数，计算中使用的物性参数如表 1所示。其中液相物性参数与实验中使用的0.5 mol/L氧化锌前驱体溶液的物性参数保持一致，气相的物性参数为25 ℃、一个标准大气压(1.013×105 Pa)下空气的物性参数。计算中采用的U取值范围为5 000~20 000 μm/s，对应毛细数取值范围为6.5×10-4~2.6×10-3。
表 1 相场动力学方法中使用的物性参数 |
| 物性 | 液相 | 气相 |
| 密度/(kg·m-3) | 965 | 1.29 |
| 动力黏度/(mPa·s) | 4.38 | 0.01 |
| 表面张力/(mN·m-1) | 33.44 | 33.44 |
使用COMSOL Multiphysics 6.0进行数值求解。在气液界面可能运动到达的区域采用非结构化四面体网格，并且对网格进行了加细处理，而在其他区域采用结构化六面体网格以减小整体网格数，提高计算效率。细化的非结构化四面体网格区域四面体的最长边长为4.6 μm。
2 润滑近似方法
2.1 域内方程
其中：t为时间；pa表示环境压力；ρl表示液体的密度；Δ=∂2/(∂x2)+∂2/(∂y2)为Laplace算符；П(h)为分离压，用于模拟三相接触线相对有限接触角表面的运动。参考文[17]，将分离压取为
其中：n、m为给定整数，本文中取n=3、m=2；h0为前驱膜厚度，本文中取h0=50 nm。
2.2 边界条件与初始条件
参考图 3，AB、CD边界采用对称边界条件：
对于BC边界，该处为前驱膜，边界条件设置为：
根据渐进匹配的思想，外部静态弯月面解在接触线附近应该能够平滑过渡到数值计算求解域内的解，即对于通过式(10)求解得到的弯月面轮廓与在数值求解域内求解得到的液膜轮廓，在过渡连接区AD边界应该有相同的曲率。当大气压不变时，由Young-Laplace方程可知气液界面液相一侧的压力与曲率之间为线性关系，因此AD边界静态弯月面解给出的压力与在数值求解域内求解得到的压力应该相同。
对于外部静态弯月面解，AD边界对应h→0，∂h/∂y→-sinθ，此时有
静态弯月面区域内压力为流体静压，又y→0时静态弯月面趋向于水平液面，此处曲率为0，因此p=pa，则AD边界压力可写为
故AD边界边界条件为：
其中hres表示过渡区液膜厚度，本文中取hres=100 μm。
本节中数值求解采用的初始条件为仅在y方向发生变化而满足AD边界条件的二次函数曲面。该曲面在h=h0处截断，小于h0的取值全部取为h0，表达式如下：
其中y0为人为指定的截断点，要求0 < y0 < ymax，ymax为求解域内y坐标最大值。
2.3 几何与物性参数设置、网格划分与求解
润滑近似方法的求解域中，沿y轴方向的矩形边长设置为800 μm，沿x轴方向的矩形边长设置为200 μm。由于边AB、CD均设置为对称边界条件，因此模拟中亲水图案间距等于矩形计算域沿x轴边长的2倍，为400 μm，与实验条件保持一致。润滑近似方法设置的密度和表面张力与实验中使用的溶液物性参数保持一致，数值如表 1所示。润滑近似方法的计算中，固定U=1 000 μm/s，通过改变μ来改变计算结果对应的Ca。
使用COMSOL Multiphysics 6.0进行数值求解，求解域内采用全局自由三角形网格划分，在可能存在接触线的区域进行网格细分。
3 实验与材料
3.1 亲/疏水图案化基底制备
以氧化层厚度为100 nm的氧化硅片作为基底，对基底进行氧气等离子体处理，处理功率为70 W，通入氧气流量为100 mL/min，处理时间为3 min。在基底表面以转速3 000 r/min旋涂NR9-3000PY负性光刻胶，旋涂时间为40 s。旋涂结束后将带光刻胶基底静置于120 ℃恒温热板烘烤90 s，取下冷却后用汞灯光刻机(URE-2000，中科院光电所)曝光40 s，然后再次将基底静置于120 ℃恒温热板烘烤90 s。取下基底冷却后用RD6显影液显影35 s，获得部分区域被光刻胶覆盖的基底。
对带有光刻胶图案的基底进行氧气等离子体处理3 min(处理功率100 W，通入100 mL/min O2)，然后将基底置于培养皿中并同时滴加20 μL三甲氧基(1H, 1H, 2H, 2H-十七氟癸基)硅烷，抽真空后将培养皿置于200 ℃恒温热板加热10 min。冷却后取出基底，浸入丙酮中并在25 ℃下超声5 min，完全去除光刻胶后获得带有亲/疏水图案的基底。本文使用的基底带有边长为40 μm×80 μm的亲水矩形图案阵列，亲水图案长40 μm的边平行于提拉速度方向。亲水图案沿平行提拉方向的间距为320 μm，沿垂直提拉方向的间距为400 μm。
3.2 溶液配制与物性表征
用于生成氧化锌(ZnO)固体薄膜的0.5 mol/L前驱体溶液是通过将2.1951 g二水合乙酸锌(Zn(CH3COO)2·2H2O，纯度99.99%)和0.6108 g乙醇胺(C2H7NO，纯度99%)溶解在20 mL乙二醇单甲醚(2-MOE，纯度99%)中，然后电磁搅拌12 h制成。使用旋转流变仪Anton Paar Physica MCR302测试配制的前驱体溶液的动力黏度，使用Dataphysics DCAT21表面界面张力仪测试配制的溶液的表面张力，结果分别为4.38 mPa·s和33.44 mN/m。
3.3 浸涂与液滴最大高度表征
使用浸涂机(SYDC-100ZT，上海三研科技有限公司)将带有亲/疏水图案的基底以2 000 μm/s速度浸入0.5 mol/L氧化锌前驱体溶液中，浸泡1 min，然后在1 000~10 000 μm/s取定一个速度，以该速度匀速将基底提拉出溶液，在25 ℃、10%相对湿度环境下干燥5 min，然后将基底先后在150 ℃和400 ℃下加热5 min和1 h，使得限域吸附产生的前驱体溶液液滴完全固化为氧化锌固体薄膜。使用Bruker Multimode 8原子力显微镜测试图案化氧化锌固体薄膜表面形貌，通过MATLAB数值积分获得固体薄膜体积。
3.4 图案化固体薄膜体积与亲水区域生成液滴最大高度的关系
亲水图案区域生成液滴的体积为
其中：ρZnO表示生成固体薄膜的密度，MZnO表示生成固体薄膜的摩尔质量，Vsolid表示生成固体薄膜的体积，csolution表示前驱体溶液中锌离子的浓度。取ρZnO=5.6 g/cm3，MZnO=81.4 g/mol，csolution=0.5 mol/L，可以由固体薄膜的体积换算得到液滴的体积。
对于亲水图案区域生成的限域液滴，由于亲水区域的长度和宽度均远小于毛细长度，表征重力与表面张力作用之比的无量纲Bond数Bo=ρgw2/γ≈0.002≪1，因此重力对液滴形状的影响可以忽略不计，在平衡态下液滴表面为给定体积和铺展区域条件下的表面积最小曲面。利用COMSOL 6.0求解该带有积分约束的变分问题，可由换算得到的Vfilm和亲水图案形状求得生成液滴的最大高度。
3.5 高速摄像机拍摄限域吸附过程
使用高速摄像机(X190，合肥中科君达视界技术股份有限公司)拍摄亲/疏水图案化基底限域吸附液滴的动态过程，拍摄帧率为10 000帧/s。
4 结果与讨论
4.1 限域吸附动力学过程
基于润滑近似方法和相场动力学方法分别在Ca=10-4条件下求解并绘制40 μm×80 μm亲水区域限域吸附动力学过程，并与高速摄像机拍摄的40 μm×80 μm亲水图案区域限域吸附过程(U=10 000 μm/s，对应毛细数Ca≈1.31×10-3)对比，结果如图 4所示。由于不同毛细数条件下高速摄像机拍摄得到限域吸附动力学过程的特征不存在显著的定性差异，而较低毛细数条件下生成卫星滴较小，因此为获得更好的比较效果，采用更高毛细数条件下的高速摄像机拍摄得到限域吸附过程与仿真计算结果对比。可以看出，亲水区域的限域吸附动力学过程包括2个阶段：第1个阶段，亲水图案区域边缘钉扎接触线，导致该区域接触线随亲水图案一起上升，然而由于亲水图案宽度有限，图案两侧接触线并未随亲水图案一起上升，仅仅在靠近亲水图案的区域略微抬升，在远离图案的区域仍然保持接近平直的状态；第2个阶段，当亲水图案被完全提拉出弯月面后，仍有部分液体存在于图案和弯月面之间形成液桥，随着图案继续上升，由于图案下部疏水区域不能钉扎接触线，液桥两侧接触线将向图案中部收缩，致使液桥宽度逐渐减小，最终发生断裂并在断裂后产生若干卫星滴[26]。由图 3c的高速摄像拍摄结果可知，液桥首先在靠近亲水图案一侧发生断裂，然后在极短的时间内失稳并分裂为若干卫星滴。
将2种数值计算方法得到的结果与实验结果对比，可以发现在限域吸附动力学过程的第1个阶段，2种数值方法得到的结果与实验结果保持一致；在图案浸涂的第2个阶段，润滑近似方法给出的液桥断裂和卫星滴形成过程与高速摄像拍摄结果保持高度一致，然而相场动力学方法并未给出卫星滴形成过程，这是由于液桥断裂后形成卫星滴尺寸过小，与相场方法中的扩散界面厚度在同量级甚至更小，超出了相场方法的相界面计算分辨率。加细网格可以降低计算中使用的扩散界面厚度，提高方法的相界面分辨率，但同时也会增加计算量。由于算力受限，本文并未加细网格尝试分辨卫星滴，但后续液滴最大高度随毛细数变化规律的结果表明，即使不能分辨出卫星滴，相场动力学方法也能给出足够准确的亲水区域限域吸附产生液滴的最大高度随毛细数变化规律。
4.2 液滴最大高度变化规律
5 结论
本文通过相场动力学和润滑近似方法分别对亲/疏水图案化基底限域吸附液体过程进行了数值建模，然后基于测量限域吸附溶液干燥后固体残留形貌的方法对限域吸附液滴最大高度进行了实验研究，最后对比了理论和实验结果，主要结论如下：
1) 基于润滑近似的模拟方法和基于相场动力学的方法均能模拟出亲水区域限域吸附的2个阶段的动力学过程。相比相场方法，润滑近似方法能够更准确地刻画液桥断裂和卫星滴形成过程。
2) 相场方法相比润滑近似方法给出了与实验结果更加吻合的吸附液滴最大高度与毛细数之间的关系，揭示了亲水区域限域吸附液体过程中黏性力对形成液膜的作用减弱的规律，有望用于指导限域吸附液滴微阵列中液滴最大高度的精准调控。
