1.   引　言

双重乳液是一种分散相液滴中包裹着更小液滴的结构化流体, 其独特的壳-核结构使之具有更高的灵活性和可控性, 在化学反应^[1,2]、功能材料合成^[3-5]、药物输运^[6,7]等领域有着广泛的应用前景. 自20年前问世以来, 微流控技术获得了长足发展, 可以实现对多相流体在微尺度上的精确操控, 并且具有装置高度集成、实验试剂消耗少等优势. 因此, 微流控技术已经成为制备双重乳液的一种有效手段, 并引起了众多学者的广泛关注^[8-10]. 然而, 为了保证对多相流体的精确和稳定操控, 微流控装置内各相工质通常处于流速较低的稳定层流状态, 这在一定程度上限制了微流控乳液制备的产率. 为了实现乳液的高通量制备, 研究者们多利用微流控分裂结构对制备生成的乳液进行分裂倍增, 从而实现乳液产率的提升^[11,12]. 因此, 深入认识微流控分裂结构中乳液分裂过程的流体动力学机理十分必要, 这不仅对完善乳液多相流体动力学理论具有重要的科学意义, 并且可为工程中微流控乳液制备产率的提升提供关键的技术支撑.

为实现乳液尺寸的精确操控和乳液产量的倍增, 研究者们设计出多种微流控分裂流道结构, 如十字通道^[13]、T型通道^[14-23]、Y型通道^[24-29]、带阻碍结构通道^[14,30]等. Link等^[14]首次将T型微通道用于液滴被动破裂, 并将一系列T型通道相连形成树状网络, 可在不改变高分散相体积分数的情况下实现小液滴的大量制备. 自Link以后, 众多学者继续致力于通过实验研究T型微通道内液滴破裂特性和机理. 研究结果表明: 液滴在T型微通道内的被动破裂存在两种流型, 分别为隧道破裂流型和阻塞破裂流型^[16]; T型微通道中液滴破裂与否, 主要取决于液滴初始长度和连续相流体的毛细数, 两者之间存在指数关系, 以此表征破裂与不破裂之间的流型临界线^[14,16,21]. 考虑到已有的微流控装置制备双重乳液存在产率低的问题, 并且T型微通道不易实现双重乳液分裂, Abate和Weitz^[26]将Link^[14]等的工作进行了拓展, 将Y型通道作为树状网络的单元结构, 实现了液滴和双重乳液的高通量制备. 其研究结果表明, 主通道中连续相流速增大, 双重乳液形变加快, 有助于双乳液分裂. 最近, 不同于传统的2D光刻技术, 新型的3D对称毛细管被用于液滴和双重乳液的分裂, 可实现高质量的乳液可控倍增^[31]. 然而, 针对双重乳液在特征几何结构微通道内的流动破裂特性实验研究还未见到. 由于双重乳液相比单液滴具有更复杂的核壳结构^[32], 因此有必要对双重乳液变形破裂过程中内外液滴界面的相互影响开展深入讨论与分析.

与实验研究相比, 计算流体动力学方法为深入探究双重乳液多相流动特性提供了另一个研究途径, 此方法可以详细给出实验研究中较难获得的速度分布、压力梯度以及界面形变的实时情况^[33]. 而充分了解这些流场和压力场的演化规律将有助于优化微流控装置的设计. 目前, 相场方法已被用于模拟研究T型分叉处的液滴流动破裂特性^[15,23]. 研究表明, 液滴是否破裂与连续相和离散相黏度比、毛细数和液滴长度有着重要关联^[15,23], 在大黏度比和大毛细数下液滴更易于破裂, 而液滴尺寸减小使得液滴不易破裂^[23]. 此外, Carlson等^[25]运用该方法研究Y型分叉微通道中液滴破裂特性, 其工况分为破裂工况和不破裂工况, 两种工况的临界线与液滴长度以及毛细数相关, 该研究认为Rayleigh-Plateau不稳定性是液滴破裂的驱动机制. 梁宏等^[27]也采用此方法模拟了液滴在Y型分叉微通道中的迁移过程. 结果表明: 当毛细数足够小时, 液滴则滞留在分叉口处, 不发生破裂. 除了相场方法, 已有基于体积分数法(VOF)方法研究对称T型微通道中液滴破裂和不破裂特性的工作, 与实验研究结果类似, 发现了三种不同破裂机制, 分别为隧道破裂、不连续阻塞破裂和阻塞破裂^[18,19,22], 液滴破裂特性与毛细数以及液滴长度密切相关^[17-19,22], 毛细数增大会加快液滴破裂^[20,22]. 随后, 该方法还被用于研究不同分叉角度微通道中的液滴破裂特性, 结果表明, 两个子通道之间的夹角越小, 液滴分裂得更快^[29]. 此外, 已有文献基于数值方法研究双重乳液在不同制备结构中的生成过程^[34-39]、在剪切拉伸流中的变形破裂过程^[32,40,41]、在剪切拉伸流中双重乳液的碰撞过程^[42,43]等, 但尚缺乏采用数值方法全面深入了解双重乳液在Y型分叉微通道内流动破裂特性的研究.

综上所述, 具有Y型分叉结构的微通道为精确操控乳液尺寸和高通量制备微小乳液提供了有效途径. 与T型微通道相比, Y型微通道更利于单乳液分裂, 并已被成功应用于双重乳液分裂实验. 因此, 充分了解Y型微通道内双重乳液破裂这一过程有着重要的科学意义和工程应用价值. 然而, 现有研究主要集中于单乳液在Y型分叉通道中流动破裂特性的研究. 此外, 与单乳液相比, 双重乳液中内液滴的存在增加了内外液滴间的相互作用^[32,40,41]. 目前, 关于内液滴对双重乳液在Y型微通道中流动破裂特性影响的研究还有所不足, 其内在机理有待进一步阐明. 为此, 本文基于VOF方法模拟研究了双重乳液在Y型微通道中破裂的机理, 详细分析了双重乳液流经Y型微通道时的流场信息以及双重乳液形变参数的变化, 并给出了Y型微通道内双重乳液运动相图.

2.   数学模型

本文建立了三相流体流动的二维数学模型以研究双重乳液在Y型分叉微通道内的动力学行为特性. 计算区域如图1所示, 通道尺寸为w₀ = 200 ${\text{μ}} {\rm{m}}$, w₁ = 0.75w₀, l₀ = 12w₀, l₁ = 10w₀, l₂ = 8w₀, $\theta$ = 30°. 内相流体(密度${\rho _{\rm{i}}}$, 黏度${\mu _{\rm{i}}}$)、中间相流体(密度${\rho _{\rm{m}}}$, 黏度${\mu _{\rm{m}}}$)和外相流体(密度${\rho _{\rm{o}}}$, 黏度${\mu _{\rm{o}}}$)由主通道入口处注入并且在表面张力的作用下形成双重乳液, 在外相流体的带动下向下游运动. 内液滴和外液滴的尺寸由内相和中间相流体的注入时间控制.

2.1   控制方程

不相混且不可压缩的三相流体流动可用VOF表征. 在流场中定义一个体积率函数${\alpha _j}$. 计算网格中${\alpha _j}$ = 0, 表示网格内不存在j相流体; 网格中${\alpha _j}$ = 1, 表示网格内充满j相流体; 0 < ${\alpha _j}$ < 1表示网格内为相界面. 此方法通过研究网格中每一相的体积率函数${\alpha _j}$来确定相界面, 能捕捉各种不规则形状的自由表面, 可较好地模拟双重乳液在流场中的生成^[34-39]、变形^[32,40,41]、碰撞^[42,43]等复杂的变化. 在每个网格中, 各相流体的体积率函数总和为1, 即${\alpha _{\rm{i}}} + {\alpha _{\rm{m}}} + {\alpha _{\rm{o}}} = {\rm{ }}1$. 在整个计算区域内满足以下控制方程:

体积率函数方程

$$\frac{{\partial {\alpha _j}}}{{\partial t}} + \nabla \cdot \left( {{{U}} \cdot {\alpha _{\rm{j}}}} \right) = 0,$$

(1)

连续性方程

$$\nabla \cdot {{U}} = 0,$$

(2)

动量方程

$$\begin{split} \frac{{\partial \left( {\rho {{U}}} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho {{UU}}} \right) = &- \nabla p + \nabla \cdot \mu ( \nabla {{U}} \\ &+ \nabla {{{U}}^{\rm{T}}} )+ \rho {{g}} + {{{F}}_{\text{σ}}} ,\end{split}$$

(3)

式中t为时间(s), U为速度(m/s), p为压力(Pa), $\rho$为密度(kg/m³), $\mu$为黏度(Pa·s), g为重力加速度(m/s²), ${{{F}}_{\text{σ}}}$为表面张力(N/m³). 因为通道尺度为微米级, 所以重力加速度的影响可忽略不计. $\rho$和$\mu$可由下式计算:

$$\rho = {\alpha _{\rm{i}}}{\rho _{\rm{i}}} + {\alpha _{\rm{m}}}{\rho _{\rm{m}}} + {\alpha _{\rm{o}}}{\rho _{\rm{o}}},$$

(4)

$$\mu = {\alpha _{\rm{i}}}{\mu _{\rm{i}}} + {\alpha _{\rm{m}}}{\mu _{\rm{m}}} + {\alpha _{\rm{o}}}{\mu _{\rm{o}}}.$$

(5)

源项${{{F}}_{\text{σ}}}$采用连续表面力(CSF)方法^[44]计算,

$${{{F}}_{\text{σ}}} = \sigma \kappa \hat n{\delta _{\rm{s}}},$$

(6)

其中$\sigma$为表面张力系数, $\kappa$为界面平均曲率, ${\delta _{\rm{s}}}$为狄拉克方程. 另外包含壁面的网格中$\hat {{n}}$与接触角${\theta _{\rm{w}}}$有关,

$${\hat{ n}} = {{\hat{ n}}_{\rm{w}}}\cos {\theta _{\rm{w}}} + {{\hat{ t}}_{\rm{w}}}\sin{\theta _{\rm{w}}},$$

(7)

其中${{\hat{ n}}_{\rm{w}}}$和${{\hat{ t}}_{\rm{w}}}$分别为垂直于和平行于壁面的单位向量. 本文中通道壁面对中间相流体是不可润湿的, 故${\theta _{\rm{w}}}$ = 180°.

本文中各相流体的物性参数如表1所列. 内液滴界面(${\sigma _{\rm{i}}}$)和外液滴界面(${\sigma _{\rm{o}}}$)的表面张力系数分别为0.0073和0.0226 N/m. 外相流体毛细数($Ca =$${\mu _{\rm{o}}}{u_0}/{\sigma _{\rm{o}}}$, u₀为主通道入口处流体速度)的范围为0.005 ≤ Ca ≤ 0.03. 外液滴与内液滴之间的体积比为V_oi = V_o/V_i, 主通道中双重乳液的无量纲初始长度定义为l^* = l_e/w₁. 采用无量纲参数定量描述双重乳液在分叉处的形貌演化过程, 分别为图1中虚线框内所示的内、外液滴颈部厚度$\delta _{{\rm{in}}}^* = {\delta _{{\rm{in}}}}/{w_1}$, $\delta _{{\rm{out}}}^* = {\delta _{{\rm{out}}}}/{w_1}$. 内、外液滴前端在分支通道内的运动距离$\Delta l_{{\rm{in}}}^* = {\rm{ }}({l_{{\rm{in}}}} - {l_{{\rm{in}},0}})/{w_1}$, $\Delta l_{{\rm{out}}}^* = {\rm{ }}({l_{{\rm{out}}}} - {l_{{\rm{out}},0}})/{w_1}$. 此外定义了无量纲时间t^* = u₀(t – t₀)/w₁, t₀为初始时刻.

2.2   边界条件

因为数值模拟中雷诺数小于60, 所以本文采用层流模型. 主通道入口处采用周期性边界条件以控制各相流体注入通道的体积. 如图2所示, 主通道入口处流体速度为匀速u₀. 双重乳液形成的具体步骤如下: 步骤(a), 当0 < t < t_m, 中间相流体流入; 步骤(b), 当t_m < t < t_m+i, 内相和中间相流体同时流入(w_in = 0.8w₀); 步骤(c), 当t_m+i < t < t_m+i+m, 中间相流体流入使外液滴包裹内液滴, 从而形成一个双重乳液; 步骤(d), 当t_m+i+m < t, 外相流体流入.

通道出口处为一个大气压(1 atm = 1.013 × 10⁵ Pa). 通道的内壁面采用无滑移边界条件

$${u_\varGamma } = {v_\varGamma } = {w_\varGamma } = {\rm{0}},$$

(8)

其中$\varGamma$表示内壁面.

2.3   数值求解方法

采用有限体积差分法同时离散求解VOF模型中的控制方程. 求解时, 采用二阶迎风差分格式离散动量方程, 压力场与速度场之间的耦合采用SIMPLE算法, 压力场的离散采用体积力分数法, 相界面采用分段线性界面重构法(PLIC)^[45]进行几何重构. 控制方程离散后, 采用Gauss-Seidel方法进行迭代求解. 为了保证计算的收敛性和收敛速度, 各项松弛因子设为: 0.2(压力), 0.3(密度), 0.3(体积力)和0.2(动量). 采用相对准则判断收敛, 即当每个计算步长内各计算单元中所有变量的相对残差之和小于初始值0.1%时判定计算收敛. 计算中的时间步长依据全局库朗数不大于0.2的规则进行自动调整.

本文采用四边形网格对图1所示的计算区域进行网格划分. 对近壁面附近的网格进行局部加密以确保能捕捉到双重乳液与壁面之间的液膜. 计算前采用四套不同的网格数进行网格独立性检验. 如图3所示, 当计算单元总数大于315610时, 随着网格数的增加, 双重乳液形貌基本重合. 考虑计算成本, 本文采用的网格数为315610.

2.4   模型验证

为验证本文所建立的数学模型的正确性, 依据文献[41]中的实验, 基于上文的二维数学模型数值模拟了如图4所示的剪切流场中双重乳液形变, 并且与实验结果进行了比对. 模拟中内、中、外三相流体分别为去离子水(20 ℃下${\mu _{\rm{i}}}$ = 0.001 Pa·s, ${\rho _{\rm{i}}}=$ 998 kg/m³)、Ucon润滑油(20 ℃下${\mu _{\rm{m}}}$ = 0.125 Pa·s, ${\rho _{\rm{m}}}$ = 991 kg/m³)和不同分子量的硅油混合物(20 ℃下${\mu _{\rm{o}}}$ = 5.2 Pa·s, ${\rho _{\rm{o}}}$ = 989 kg/m³). 表面张力系数${\sigma _{\rm{o}}}$ = 0.024 N/m, ${\sigma _{\rm{i}}}$ = 0.003 N/m. 计算区域的几何尺寸为W × H = 20R_o × 8R_o, R_o为外液滴半径. 如图5所示, 对于双重乳液的形貌和内、外液滴的Taylor^[46]稳态形变参数D(D_i = (L_i – B_i)/(L_i + B_i), D_o = (L_o – B_o)/(L_o + B_o)), 数值模拟结果与实验结果吻合较好, 这表明本文的数学模型能有效地预测Y型微通道内双重乳液的形变.

3.   结果分析

基于上述多相流模型, 本文开展了双重乳液在Y型分叉微通道内流动特性的研究. 在文献[16, 22, 23]中, 液滴在流经T型分叉微通道时会呈现三种典型流型, 分别为阻塞破裂(obstructed breakup)、隧道破裂(tunnel breakup)和不破裂(non-breakup)流型. 与液滴情况类似, 双重乳液流经Y型分叉微通道时, 同样存在这三种典型的流型. 本文针对这三种典型流型的动力学信息(界面形貌、压力分布、速度分布等)瞬时演化规律进行详细分析, 以探究双重乳液破裂过程中的驱动和阻碍作用变化情况, 揭示双重乳液破裂的机理.

3.1   阻塞破裂

当双重乳液流经Y型分叉处时, 上游压力驱动双重乳液形变, 此时双重乳液前端以及尾部界面形貌发生变化(图6(a)), 而双重乳液界面张力将阻碍该过程的发生, 该过程中速度场演化如图6(b)所示. 为阐明上游驱动压力以及双重乳液界面张力之间相互关系, 图6—8分别给出了阻塞破裂工况中流场整体压力分布情况、双重乳液界面张力以及上游压力瞬时演化曲线. 如图6(a)所示, 可将双重乳液流经Y型分叉处分为三个阶段: entering阶段、squeezing阶段以及post-breakup阶段.

阻塞破裂工况中的entering阶段(图6(a)中$t_{{\rm{OB}}1}^*$—$t_{{\rm{OB}}3}^*$时间段)定义为从双重乳液前端离开主通道进入分叉子通道开始, 到乳液前端接触Y型分叉尖角时刻为止. 在阻塞破裂工况entering阶段中, 双重乳液前端进入主通道与分叉子通道连接处, 该过渡空间中的通道宽度逐渐增加, 双重乳液受限程度降低. 此时, 双重乳液前端界面在前进过程中其曲面半径逐渐增大, 即前端界面张力$\Delta {p_{\sigma ,{\rm{front}}}}$逐渐减小(见图7(a)中$t_{{\rm{OB}}1}^*$—$t_{{\rm{OB}}3}^*$时间段). 双重乳液前端界面张力$\Delta {p_{\sigma ,{\rm{front}}}}$方向与流向相反, 并且阻碍双重乳液形变, 因此$\Delta {p_{\sigma ,{\rm{front}}}}$减小意味着流阻减小. 在阻塞破裂工况entering阶段中, 双重乳液尾部一直处于主通道中, 界面形貌没有变化, 故其界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$维持不变(见图7(b)中$t_{{\rm{OB}}1}^* - t_{{\rm{OB}}3}^*$时间段). 双重乳液尾部界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$方向沿主流方向, 即$\Delta {p_{\sigma ,{\rm{tail}}}}$增加意味着流阻减小. 图7(c)给出了$\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$的演化曲线以定量描述双重乳液流动所受的阻碍作用. 在阻塞破裂工况entering阶段中, $\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$呈下降趋势, 在$t_{{\rm{OB}}3}^*$时刻到达最低点. 这与图8上游压力p_inlet演化规律一致, 且在该阶段中, $\Delta {p_{\sigma ,{\rm{front}}}} -$$\Delta {p_{\sigma ,{\rm{tail}}}}$与p_inlet下降压差均为60 Pa.

进入squeezing阶段(图6(a)中$t_{{\rm{OB}}3}^*$—$t_{{\rm{OB}}5}^*$时间段), 双重乳液前端分别进入两个子通道. 在开始阶段(图6(a)中$t_{{\rm{OB}}3}^*$—$t_{{\rm{OB}}5}^*$时间段), 前端界面曲面半径逐渐减小, 前端界面张力$\Delta {p_{\sigma ,{\rm{front}}}}$增大(图7(a)中$t_{{\rm{OB}}3}^*$—$t_{{\rm{OB}}5}^*$时间段). 此后, 由于双重乳液阻塞子通道, 前端界面曲面半径约等于子通道半径, 故$\Delta {p_{\sigma ,{\rm{front}}}}$维持不变(图7(a)中$t_{{\rm{OB}}5}^*$—$t_{{\rm{OB}}8}^*$时间段). 双重乳液尾部在该阶段先进入主通道与分叉子通道的连接处, 然后随着阻塞破裂过程发展, 尾部界面一分为二, 完全进入子通道. 在进入过渡区域的初始阶段(图7(b)中$t_{{\rm{OB}}4}^*$—$t_{{\rm{OB}}5}^*$时间段), 尾部界面在前进过程中其曲面半径迅速增加, 界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$陡降. 此后, 在$t_{{\rm{OB}}5}^*$—$t_{{\rm{OB}}7}^*$时间段, 尾部界面沿着渐扩通道贴近壁面前进, 该过程中的尾部界面曲面半径缓慢增加, 界面张力$\Delta {p_{\sigma {\rm{tail}}}}$逐渐减小. 当尾部界面接近分叉尖角, 如$t_{{\rm{OB}}7}^*$时刻, 乳液尾部近乎完全进入子通道, 但此时双重乳液还未完全破裂, 特别是乳液尾部在子通道中拉伸, 界面趋于平直, 因此界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$在$t_{{\rm{OB}}7}^*$—$t_{{\rm{OB}}8}^*$时间段迅速下降趋向0 Pa. 结合$\Delta {p_{\sigma ,{\rm{front}}}}$和$\Delta {p_{\sigma ,{\rm{tail}}}}$演化规律, 纵观squeezing阶段, 流动受阻参数的$\Delta {p_{\sigma ,{\rm{front}}}} -$$\Delta {p_{\sigma ,{\rm{tail}}}}$先增大至$t_{{\rm{OB}}5}^*$时刻, 随后缓慢增加至$t_{{\rm{OB}}7}^*$时刻, 最后阶段出现阶跃变化, 该阶段中, $\Delta {p_{\sigma ,{\rm{front}}}}-$ $\Delta {p_{\sigma ,{\rm{tail}}}}$增加压力324 Pa. 其演化曲线与图8上游压力p_inlet演化规律一致, 且该阶段中, $\Delta {p_{\sigma ,{\rm{front}}}} -$$\Delta {p_{\sigma ,{\rm{tail}}}}$与p_inlet上升压差为369 Pa, 两者变化数值接近. 综上所述, 阻塞破裂工况中上游压力p_inlet与乳液两端界面张力差$\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$正相关.

当双重乳液完全破裂, 进入post-breakup阶段, 双重乳液尾部回缩, 界面曲面半径减小, 尾部界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$增加, 从而流阻减小, 上游压力p_inlet相应降低. 该预测与图8中p_inlet变化规律相一致.

为深入理解双重乳液破裂的机理, 采用无量纲特征参数定量描述squeezing阶段中双重乳液形貌(包括内液滴以及外液滴)的演化过程(见图9), 分别为外液滴颈部厚度$\delta _{{\rm{out}}}^*$、外液滴前端运动距离$\Delta l_{{\rm{out}}}^*$、内液滴颈部厚度$\delta _{{\rm{in}}}^*$和内液滴前端运动距离$\Delta l_{{\rm{in}}}^*$. $t_{{\rm{out}},0}^{\rm{*}}$以及$t_{{\rm{in}},0}^{\rm{*}}$分别代表外液滴和内液滴前端接触分叉尖角的时刻, 而$\delta _{{\rm{out}},0}^{\rm{*}}$和$\delta _{{\rm{in}},0}^{\rm{*}}$分别为所对应时刻的外液滴和内液滴颈部厚度. 如图9(a)所示, 外液滴颈部厚度随时间减小直至双重乳液破裂, 插图中给出了$\delta _{{\rm{out}},0}^{\rm{*}} - \delta _{{\rm{out}}}^*$和${t^*} - t_{{\rm{out}},0}^{\rm{*}}$的对数坐标图. 根据指数标度律可判定squeezing阶段包含两个子过程. 子过程I中$\delta _{{\rm{out}},0}^{\rm{*}} - \delta _{{\rm{out}}}^*$和${t^*} - t_{{\rm{out}},0}^{\rm{*}}$呈线性关系, 而在子过程II中两者呈指数关系$\delta _{{\rm{out}},0}^{\rm{*}} - \delta _{{\rm{out}}}^*$ ~ (${t^*} - t_{{\rm{out}},0}^{\rm{*}}$)^4/5, 这与T型微通道中液滴阻塞破裂的特性类似^[23]. 两个子过程中颈部厚度随时间的变化规律不同, 这主要是由于外液滴尾部界面张力发生了改变. 与颈部厚度变化不同的是, 外液滴前端运动距离随时间线性增加, 有$\Delta l_{{\rm{out}}}^*$ ≈ $2/3({t^*} -$$t_{{\rm{out}},0}^{\rm{*}})$ (见图9(b)). 此时, 外液滴前端以2u₀/3恒定速度运动, 这主要由于该工况中子通道完全被外液滴阻塞所导致. 图9(c)和图9(d)分别给出了内液滴颈部厚度$\delta _{{\rm{in}}}^*$和内液滴前端运动距离$\Delta l_{{\rm{in}}}^*$随时间的演化曲线, 其演化特性与外液滴一致. 内液滴颈部厚度随时间缩小, 也分为两个子过程, $\delta _{{\rm{in}},0}^{\rm{*}} - \delta _{{\rm{in}}}^*$和${t^*} - t_{{\rm{in}},0}^{\rm{*}}$分别呈线性关系和指数关系$\delta _{{\rm{in}},0}^{\rm{*}} - \delta _{{\rm{in}}}^*$ ~ (${t^*} - t_{{\rm{in}},0}^{\rm{*}}$)^4/5; 内液滴前端运动距离$\Delta l_{{\rm{in}}}^*$也随时间线性增加, 有$\Delta l_{{\rm{in}}}^*$ ≈ 0.7(${t^*} - t_{{\rm{in}},0}^{\rm{*}}$). 此时, 内液滴前端以0.7u₀匀速前进.

3.2   隧道破裂

从流动行为角度分析, 隧道破裂流型在entering和post-breakup阶段与阻塞破裂流型非常类似, 但在squeezing阶段两者有明显区别(见图10). 在squeezing阶段$t_{\rm {TB4}}^*$时刻, 隧道破裂流型中乳液与子通道壁面间开始出现两条隧道. 从压力演化角度分析, 隧道破裂流型在entering和squeezing阶段均展现出一定的不同(见图11和图12).

在entering阶段(图10(a)中$t_{{\rm{TB}}1}^*$—$t_{{\rm{TB}}3}^*$时间段), 当双重乳液前端开始进入渐扩过渡通道时, 由于其受限程度降低, 双重乳液前端界面在前进过程中其曲面半径逐渐增大, 即前端界面张力$\Delta {p_{\sigma ,{\rm{front}}}}$逐渐减小(见图11(a)中$t_{{\rm{TB}}1}^*$—$t_{{\rm{TB}}2}^*$时间段). 值得注意的是, 由于该工况中双重乳液长度较短, 在entering阶段, 双重乳液尾部已离开主通道进入渐扩过渡通道(图11(b)中$t_{{\rm{TB}}2}^*$—$t_{{\rm{TB}}3}^*$时间段), 在前进过程中尾部界面曲面半径迅速增加, 界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$陡降, 在$t_{{\rm{TB}}3}^*$时刻变化趋于平缓. 因此, 不同于阻塞破裂工况, 在隧道破裂entering阶段, $\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$先下降后上升, 在$t_{{\rm{TB}}2}^*$时刻为最低点(见图11(c) $t_{{\rm{TB}}1}^*$—$t_{{\rm{TB}}3}^*$时间段), 这与图12上游压力p_inlet演化规律一致. 在该阶段中, $\Delta {p_{\sigma ,{\rm{front}}}} -$ $\Delta {p_{\sigma ,{\rm{tail}}}}$先下降37 Pa后上升70 Pa, 而p_inlet先下降24 Pa后上升37 Pa. 两者之间略有差异, 这可能是由于在隧道破裂工况中, 双重乳液并不完全阻塞渐扩过渡通道(见图10(a) $t_{{\rm{TB}}3}^*$时刻). 在entering阶段, 隧道破裂流型速度场与阻塞破裂流型非常类似, 如图10(b)所示.

进入squeezing阶段(图10(a)中$t_{{\rm{TB}}3}^*$—$t_{{\rm{TB}}5}^*$时间段), 双重乳液前端分别进入两个子通道. 在开始阶段(图10(a)中$t_{{\rm{TB}}3}^*$—$t_{{\rm{TB}}4}^*$时间段), 前端界面曲面半径逐渐减小, 前端界面张力$\Delta {p_{\sigma ,{\rm{front}}}}$迅速增大(图11(a)中$t_{{\rm{TB}}3}^*$—$t_{{\rm{TB}}5}^*$时间段). 此后不同于阻塞破裂工况, 由于隧道的存在, 双重乳液在拉伸破裂过程中隧道宽度增加, 双重乳液前端半径会继续减小. 因此, 从$t_{{\rm{TB}}4}^*$时刻开始, 前端界面张力$\Delta {p_{\sigma ,{\rm{front}}}}$仍会增大, 其增速有所降低, 略为平缓(图11(a)中$t_{{\rm{TB}}4}^*$—$t_{{\rm{TB}}7}^*$时间段). 在squeezing阶段, 双重乳液尾部继续在渐扩过渡通道中前进, 直至尾部界面在分叉尖角处破裂, 尾部界面曲线在该过程中基本保持不变, 界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$略有下降(图11(b)中$t_{{\rm{TB}}3}^*$—$t_{{\rm{TB}}7}^*$时间段). 值得注意的是, 由于当尾部界面接近分叉尖角时, 隧道的存在使得尾部界面仍然保持一定弧度, 没有被拉伸趋于平直. 因此, 在隧道破裂工况squeezing阶段的最后时刻, $\Delta {p_{\sigma ,{\rm{tail}}}}$没有出现陡降现象, 只是略有下降. 在squeezing阶段, 隧道破裂流型由于隧道的出现, 连续相流体从隧道流过, 此时隧道处速度明显增大, 如图10(b)所示.

结合$\Delta {p_{\sigma ,{\rm{front}}}}$和$\Delta {p_{\sigma ,{\rm{tail}}}}$的演化规律, 在隧道破裂工况squeezing阶段, $\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$先迅速上升至$t_{{\rm{TB}}4}^*$时刻, 随后上升速度略有放缓, 在最后阶段增加速率略有加快. 该阶段$\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$增加了227 Pa. 由图12可知, $\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$演化曲线与上游压力p_inlet演化规律基本一致, 在该阶段, p_inlet上升压差为305 Pa, 大于$\Delta {p_{\sigma ,{\rm{front}}}} -$$\Delta {p_{\sigma ,{\rm{tail}}}}$变化数值. 两者出现差异的原因在$t_{{\rm{TB}}3}^*$—$t_{{\rm{TB}}5}^*$阶段, 可以发现p_inlet演化曲线从$t_{{\rm{TB}}3}^*$时刻开始以较高速率上升, 且没有在$t_{{\rm{TB}}4}^*$时刻减缓, 一直维持到$t_{{\rm{TB}}5}^*$时刻. 这是由于在$t_{{\rm{TB}}4}^* - t_{{\rm{TB}}5}^*$阶段, 双重乳液内液滴被挤压变形, 外液滴界面逐渐接近内液滴, 这使得外液滴尾部内侧压力上升, 这对其变形起到了阻碍作用. 在阻塞破裂工况中, squeezing阶段p_inlet上升压差也高过$\Delta {p_{\sigma ,{\rm{front}}}} - \Delta {p_{\sigma ,{\rm{tail}}}}$的变化45 Pa, 由于在阻塞破裂工况中内外液滴距离较大, 因此相互作用影响程度小于隧道破裂工况.

与阻塞破裂类似, 进入post-breakup阶段, 双重乳液将完全破裂, 其尾部回缩, 界面曲面半径减小, 尾部界面张力$\Delta {p_{\sigma ,{\rm{tail}}}}$增加, 流阻减小, 上游压力p_inlet相应降低(图11(b)中$t_{{\rm{TB}}3}^*$—$t_{{\rm{TB}}7}^*$时间段).

为深入探究“隧道效应”对双重乳液破裂行为的影响, 图13给出了隧道破裂工况中典型几何参数(包括外液滴颈部厚度$\delta _{{\rm{out}}}^*$、外液滴前端运动距离$\Delta l_{{\rm{out}}}^*$、内液滴颈部厚度$\delta _{{\rm{in}}}^*$和内液滴前端运动距离$\Delta l_{{\rm{in}}}^*$)在squeezing阶段随时间演化的情况. 与阻塞破裂类似, $\delta _{{\rm{out}},0}^{\rm{*}} - \delta _{{\rm{out}}}^*$随时间演化分为两个子过程, 子过程I中$\delta _{{\rm{out}},0}^{\rm{*}} - \delta _{{\rm{out}}}^*$和${t^*} - t_{{\rm{out}},0}^{\rm{*}}$呈线性关系, 而在子过程II中两者呈指数关系(见图13(a)插图). 不同之处在于, 子过程II的指数约为3.75/5, 小于阻塞破裂的4/5. 这表明隧道的出现使得双重乳液外液滴颈部收缩速率降低, 这也间接说明隧道内润滑流动引起的剪切作用并没促进外液滴的挤压形变. 图13(b)给出了外液滴前端运动距离随时间的演化曲线, 如图所示, 在隧道破裂流型中, $\Delta l_{{\rm{out}}}^*$演化趋势前期为线性变化, 随后变为指数变化, 指数为4/5. 这表明隧道的出现也降低了外液滴沿流向的拉伸速率.

图13(c)和图13(d)分别给出了内液滴颈部厚度$\delta _{{\rm{in}}}^*$和内液滴前端运动距离$\Delta l_{{\rm{in}}}^*$随时间的演化曲线, 其演化特性与该工况中外液滴的类似. 不同之处在于特征参数的指数关系分别为$\delta _{{\rm{in}},0}^{\rm{*}} - \delta _{{\rm{in}}}^*$ ~ (${t^*} - t_{{\rm{in}},0}^{\rm{*}}$)^4/5, $\Delta l_{{\rm{in}}}^*$ ~ (${t^*} - t_{{\rm{in}},0}^{\rm{*}}$)^3.75/5, 其指数分别为4/5和3.75/5, 与外液滴的指数3.75/5和4/5均有所不同. 经讨论分析, 隧道破裂工况中内液滴颈部厚度演化曲线的指数与阻塞破裂中的指数一致, 略大于外液滴的指数.

3.3   不破裂

当双重乳液尺寸较小时, Y型分叉微通道中出现不破裂流型, 此时乳液将经历三个阶段, 依次为: entering阶段、sliding阶段、和non-breakup阶段. 其中, entering阶段与隧道破裂类似, 后两个阶段是其特有的(见图14和图15). 在entering阶段, 双重乳液进入渐扩过渡通道, 由于此时双重乳液长度小于通道宽度, 外液滴界面没有显著形变, 因此上游压力p_inlet维持不变. 进入sliding阶段, 当双重乳液前端前进撞击分叉尖角, 双重乳液发生形变, 由圆形被挤压为椭圆形. 在$t_{{\rm{NB}}2}^*$—$t_{{\rm{NB}}4}^*$阶段, 双重乳液内部压力增加, 同时上游压力增加, 在$t_{{\rm{NB}}4}^*$时刻到达最大值. 在$t_{{\rm{NB}}4}^*$—$t_{{\rm{NB}}5}^*$阶段, 由于不对称, 双重乳液滑向左侧子通道, 此时双重乳液内存在不对称涡流, 双重乳液内部压力不再增加, 上游压力也开始回落. 值得注意的是, 由于双重乳液滑向单侧子通道, 两个子通道出口压力不再完全相同, 出现略微差异. 当双重乳液完全越过主通道中心线, 双重乳液内部压力以及上游压力也回落到初始水平.

3.4   相　图

图16给出了Y型微通道内双重乳液和单乳液运动相图. 图中x轴坐标为主通道中乳液无量纲长度, y轴坐标为主通道毛细数, 所得相图结构与前人研究结果类似. 图中采用幂函数${l^*} = \beta C{a^b}$描述隧道破裂和不破裂工况临界线, 其中双重乳液和单乳液运动相图临界线分别为l^* = 0.09Ca^–0.49和l^* = 0.07Ca^–0.49; 采用线性关系${l^*} = \alpha$描述隧道破裂和阻塞破裂工况临界线, 其中双重乳液和单乳液运动相图临界线分别为l^* = 1.66和l^* = 1.53. 由图16可知, 当保持Ca不变, 随着l^*的减小, 流型由隧道破裂逐渐转变为不破裂. 这是由于随着乳液尺寸减小, 乳液与壁面间易出现隧道, 从而会降低乳液形变速率, 最终使得乳液不破裂. 而当保持l^*不变, 随着Ca减小, 流型由隧道破裂逐渐转变为不破裂. 这是由于在相同表面张力情况下毛细数越大即外流体流速越大, 从而对应更大的上游压力, 而后者正是Y型微通道中乳液产生形变的主要驱动力, 因此毛细数的增加利于乳液破裂发生. 此外, 对比两个相图可知, 当单乳液变为双重乳液, 各工况的分界线均向右移动, 不破裂工况与隧道破裂工况分界线位置右移表明双重乳液比单乳液更不易破裂; 隧道破裂与阻塞破裂分界线右移则是由于相同工况下双重乳液上游压力水平高于单乳液工况, 从而更容易出现隧道.

4.   结　论

本文基于VOF液/液相界面追踪方法, 建立Y型微通道中双重乳液流动破裂非稳态理论模型, 数值模拟研究了Y型微通道内双重乳液破裂情况. 为揭示双重乳液破裂流型的内在机理, 详细分析了双重乳液流经Y型微通道时的流场信息以及双重乳液形变参数演化, 并给出了Y型微通道内双重乳液运动相图. 研究结果表明:

1) 在Y型微通道中, 双重乳液运动有三种流型, 即阻塞破裂流型、隧道破裂流型和不破裂流型; 双重乳液阻塞破裂或隧道破裂过程均可以分为entering, squeezing和post-breakup 三个阶段; 而双重乳液不破裂过程经历三个阶段为entering阶段、sliding阶段和non-breakup阶段;

2) 在阻塞破裂和隧道破裂工况的entering阶段以及squeezing阶段, 双重乳液均受上游压力驱动产生形变, 形变过程中乳液前端与尾部界面张力之差阻碍双重乳液形变破裂, 且两者正相关; 与阻塞破裂工况相比, 隧道破裂工况中隧道的出现将减缓双重乳液外液滴颈部收缩速率以及沿流向拉伸速率, 并减缓了内液滴沿流向拉伸速率, 对于内液滴颈部收缩速率影响不大; 当减小毛细数或乳液初始长度, 隧道破裂工况将过渡到不破裂工况, 此时乳液不再破裂;

3) 隧道破裂和不破裂工况临界线可以采用幂律关系式${l^*} = \beta C{a^b}$进行预测, 隧道破裂和阻塞破裂工况临界线可以采用线性关系${l^*} = \alpha$描述; 与单乳液运动相图相比, 双重乳液运动相图各工况的分界线关系式系数$\alpha$和$\beta$均相应增大.