近年来，微型制造技术的快速发展促进了微流控技术的广泛应用^［1−3］。对微流控技术来说，整个工艺的小型化使得能够用更自动化、廉价和便携的方式来快速分析微量的样品。随着在微流体装置中产生液滴（气泡）技术的迅速发展，目前已经形成了基于微液滴（气泡）的微流体化学处理的集中应用，包括微粒分散^［4］、蛋白质结晶^［5］、纳米颗粒合成和微量化学分析等。使用微液滴（气泡）技术时，由于样品或试剂被限制在微型液滴中，因此可以避免由于泰勒分散效应^［6］引起的样品稀释，并且可以改善混合性能。除此之外，液膜的存在避免了样品与容器表面的相互作用，从而消除了表面吸附导致的样品交叉污染。

T型微通道是一种结构相对简单的液滴生成装置，分散相流体受到连续相流体的横向剪切作用使分散相最终发生断裂而生成液滴。大量实验研究^［7−8］发现，改变两相流体的表观速度时在微通道中会产生多种两相流流型。Thorsen等^［9］在T型微通道中通过不断改变两相进口速度比，得到了不同尺寸的液滴。任桂香^［10］对T型微米级通道内微液滴形成过程进行了研究，发现液滴生成经历液滴生长、液滴从离散相分离、液滴从壁面分离三个过程。刘志鹏等^［11］用实验的方法研究了T型微通道内液滴流型分布，并对液滴界面的不稳定性进行了理论分析，结果表明，在液滴形成过程中，较大的黏性剪切力可以促使界面失稳以致破碎，并得出界面张力的改变引起的Mrangoin效应^［11］加剧了界面的不稳定性。郭书恒^［12］利用“流体体积”（volume of fluid，VOF）模型模拟微通道内液滴大小，并对液滴尺寸与连续相毛细数Ca和两相流量比Q之间的关系作出定量分析。由前人的研究可知，微流通道中液滴尺寸和流态的变化与Ca、Q、界面张力、两相流体的黏性以及通道壁面的浸润性有很大关系。其中液滴生成尺寸与毛细数和两相流体的流量比之间存在一定的函数关系。

基于格子玻尔兹曼方法(lattice Boltzmann method，LBM)颜色梯度模型，模拟T型微通道内液滴的生成过程，研究Ca、Q对液滴生成尺寸以及频率的影响，着重分析液滴尺寸与Ca、Q的函数关系。并以T型微通道生成的弹状液滴为基础，研究液滴在微流控梯状网络中的相位同步情况，定量分析Ca、液滴初始相位长度和旁通管数目对液滴同步程度的影响。

1 数值方法与模型验证 1.1 数值方法

本文采用LBM颜色梯度模型和D2Q9模型^［13］来进行二维T型微通道内液滴生成过程和微流控梯状网络中液滴相位同步行为的数值模拟。在颜色梯度模型中，红、蓝两种分布函数 $f_i^{\rm R}$ 和 $f_i^{\rm B}$ 分别表示两种不互溶流体， $i$ 为格子速度方向，计算域中混合流体分布函数定义为： $ {f_i} =$ $ f_i^{\rm R} + f_i^{\rm B}$ ，每种颜色的流体都分别执行碰撞步和流动步，即

$ \begin{split} \quad\quad & f_i^k\left( {{{x}} + {{{e}}_i}{\delta _t},t + {\delta _t}} \right) =\\ &f_i^k\left( {{{x}},t} \right) - \frac{{f_i^k\left( {{{x}},t} \right) - f_i^{k,{\rm eq}}\left( {{{x}},t} \right)}}{{{\tau _k}}} + {{\varOmega }}_i^k\end{split}$

(1)

式中：k = R或B表示流体颜色； ${\delta _t}$ 为时间步长，且 ${\delta _t} = 1$ ； ${{{e}}_i}$ 为 $i$ 方向上的格子速度； ${{\varOmega }}_i^k$ 为 $i$ 方向上 $k$ 相流体的碰撞算子（即扰动步）； $f_i^k\left( {{{x}},t} \right)$ 为k相流体在t时刻 ${{x}}$ 位置处的分布函数； ${\tau _k}$ 为k相流体的弛豫时间，每相流体动力黏度 ${\nu _k}$ 与 ${\tau _k}$ 的关系为 ${\nu _k} = c_{\rm s}^2\left( {{\tau _k} - 1/2} \right)$ ， ${{{c}}_{\rm{s}}}$ 为格子声速， ${{c}}_{\rm{s}}^{\rm{2}} = 1/3$ ； $f_i^{k,{\rm eq}}$ 为平衡态分布函数，即

$\quad\quad f_i^{k,{\rm eq}} = \rho {\omega _i}\left[ {1 + \frac{{{{{e}}_i} {\text{·}} {{u}}}}{{c_{\rm s}^2}} + \frac{{{{\left( {{{{e}}_i} {\text{·}} {{u}}} \right)}^2}}}{{2c_{\rm s}^4}} - \frac{{{{{u}}^2}}}{{2c_{\rm s}^4}}} \right]$

(2)

式中： ${\omega _i}$ 为权系数， ${\omega _0} = 4/9$ ， ${\omega _{1 \sim 4}} = 1/9$ ， ${\omega _{5 \sim 8}} = $ 1/36； ${{u}}$ 为流体宏观速度； $\rho $ 为计算区域内混合流体密度。

每相流体均满足质量守恒和动量守恒，即

$\quad\quad{\rho _k} = \sum\limits_i {f_i^k} = \sum\limits_i {f_i^{k,{\rm eq}}} $

(3)

$\quad\quad {\rho _k}{{u}} = \sum\limits_i {f_i^k{{{e}}_i}} = \sum\limits_i {f_i^{k,{\rm eq}}{{{e}}_i}} $

(4)

式中， ${\rho _k}$ 为k相流体的宏观密度，则计算区域内混合流体密度 $\rho = {\rho _{\rm R}} + {\rho _{\rm B}}$ 。

${{\varOmega }}_i^k$ 有助于相界面区域混合并产生界面张力，在二维LBM中其表达式为

$\quad\quad{{\varOmega }}_i^k = \frac{{{A_k}}}{2}\left| {{G}} \right|\left[ {{\omega _i}{{\frac{{\left( {{{{e}}_i} {\text{·}} {{G}}} \right)}}{{{{\left| {{G}} \right|}^2}}}}^2} - {B_i}} \right]$

(5)

式中： ${{G}}$ 表示颜色梯度， ${{G}}\left( {{{x}},t} \right) = \nabla {\rho ^N}$ ， $\rho^N $ 为两相相场； ${B_i}$ 为界面特性参数， ${B_0} = - 4/27$ ， ${B_{1 \sim 4}} = $ $ - 2/27$ ， $ {B_{5 \sim 8}} = - 5/108$ ； ${A_k}$ 为控制界面张力 $\sigma $ 的自由参数。

$\quad\quad{\rho ^N}\left( {{{x}},t} \right) = \frac{{{\rho _{\rm R}}\left( {{{x}},t} \right) - {\rho _{\rm B}}\left( {{{x}},t} \right)}}{{{\rho _{\rm R}}\left( {{{x}},t} \right) + {\rho _{\rm B}}\left( {{{x}},t} \right)}},\;\;{\rho ^N} \in \left[ { - 1,1} \right]$

(6)

界面张力 $\sigma $ 表示为

$\quad\quad\sigma = \frac{2}{9}\left( {{A_{\rm R}} + A{}_{\rm B}} \right)\tau $

(7)

式中：取 ${A_{\rm R}} = A{}_{\rm B}$ ，由给定的表面张力确定 ${A_k}$ ； $\tau $ 为混合流体的弛豫时间。

为了准确计算两相流体的弛豫时间^［14］，采用的计算式为

$\quad\quad\frac{1}{{\tau - 0.5}} = \frac{{1 + {\rho ^N}}}{{2\left( {{\tau _{\rm R}} - 0.5} \right)}} + \frac{{1 - {\rho ^N}}}{{2\left( {{\tau _{\rm B}} - 0.5} \right)}}$

(8)

扰动算子虽然产生了界面张力，但不能保证两种流体的不混溶性。为了促进相分离并维持相界面，应用Latva−Kokko等^［15］提出的改进后的重着色算法。该算法允许红、蓝相流体在界面切线处适度混合，同时保持颜色梯度分布对称。因此，该算法可以减少界面处假速度并消除前人原始着色步算法^［16］产生的格子聚集问题。根据该算法，红、蓝相流体经过重着色算子后，分布函数定义分别为

$\quad\quad f_i^{\rm R} = \frac{{{\rho _{\rm R}}}}{\rho }{f^*} + \beta \frac{{{\rho _{\rm R}}{\rho _{\rm B}}}}{\rho }{\omega _i}\cos {{\varphi _i}} $

(9)

$\quad\quad f_i^{\rm B} = \frac{{{\rho _{\rm B}}}}{\rho }{f^*} + \beta \frac{{{\rho _{\rm R}}{\rho _{\rm B}}}}{\rho }{\omega _i}\cos {{\varphi _i}} $

(10)

式中： ${f}^*$ 为经过扰动步后，在重着色预分离之前的混合流体分布函数；β为分离系数，且一般取0.7^［17］； ${\varphi _i}$ 为相场梯度（颜色梯度）和 ${{{e}}_i}$ 之间的夹角，有

$\quad\quad\cos {{\varphi _i}} = \frac{{{{{e}}_i} {\text{·}} \nabla {\rho ^N}}}{{\left| {{{{e}}_i}} \right|\left| {\nabla {\rho ^N}} \right|}}$

(11)

1.2 模型验证 1.2.1 Young−Lapace定律验证

由Young−Lapace方程可知，液滴相界面内外的压差与其曲率半径 ${R^*}$ 的关系为： ${p_{\rm{i}}} - {p_{\rm{o}}} = $ $ {\sigma /{{R^*}}}$ ，其中 ${p_{\rm{i}}}$ 、 $ {p_{\rm{o}}} $ 分别表示液滴界面内、外压力。在颜色梯度模型中，界面张力 $\sigma $ 如式（7）所示，由 ${A_k}$ 所决定，本文取 $\sigma = 0.03$ 。为了验证该模型方法的准确性，将半径不同的红相流体液滴放入计算域为 $200 \times 200$ 个格子单位的蓝相流体中，计算域边界为周期性边界条件。给定不同液滴半径模拟计算 ${p_{\rm{i}}} - {p_{\rm{o}}}$ 与 $1/{R^*}$ 的关系，如图1所示。由图中可知，两者呈正比例关系，斜率约为0.030 4，与本文给定的界面张力 $\sigma = 0.03$ 一致。这证明本文的模型符合Young−Lapace定律。

1.2.2 浸润性验证

在颜色梯度模型中，固体壁面的浸润性是通过给一个固定的相场 $\rho _{\rm s}^N$ 来实现， $\rho _{\rm s}^N \in $ $ \left[ { - {\rm{1,1}}} \right]$ ，不同 $\rho _{\rm s}^N$ 对应不同的接触角。为了验证该方法的可靠性，选取长、宽分别为100、50个格子单位长度的计算域。计算域左、右两侧是周期性边界条件，上、下为固体壁面。采用半步长反弹格式，在下方固体壁面上初始设置一个方形的红色相流体区域，其他区域为蓝色相流体，分别取 $\;\rho _{\rm s}^N = -0.75{\text{、}}$ $-0.50{\text{、}}0{\text{、}}0.50{\text{、}}0.75 $ 来模拟最终液滴接触角的变化。液滴与固体壁面接触角结果如图2所示。由图中可以看出，固体壁面的相场 $\rho _{\rm s}^N$ 与接触角 $\theta $ 呈线性关系，这说明本文采用的浸润性方法可行。当 $\rho _{\rm s}^N \leqslant - {\rm{0}}{\rm{.75}}$ 时，液滴与固体壁面的接触角已经大于140°，即液滴基本处于疏水状态，因此本文取 $\rho _{\rm s}^N = - 0.75$ 。

2 几何模型

本文采用的二维T型微通道几何结构如图3所示，图中：微通道总长度 ${L_{\rm{1}}}{\rm{ = 800}}$ 个格子单位长度；主、侧通道宽度 $d{\rm{ = 30}}$ 个格子单位长度；连续相入口距离分散相通道长度 ${L_{\rm{2}}}{\rm{ = 50}}$ 个格子单位长度；分散相入口距离主通道长度 ${L_{\rm{3}}}{\rm{ = 50}}$ 个格子单位长度。

连续相流体从水平主通道注入，流动方向为从左到右，分散相流体从旁通道垂直注入，与连续相交汇；进、出口均采用非平衡反弹格式^［18］，两相流体入口均采用速度入口，主通道出口采用压力出口。通道固体壁面采用半步长反弹格式^［19］获得非滑移边界条件。通道内所生成液滴长度用 $L$ 表示，液滴无量纲长度可表示为 $L/d$ 。

3 结果与分析 3.1 网格无关性验证

本文分别对 $6 \times 10^4$ 和 $1 \times 10^5$ 个网格进行了模拟计算。结果发现，两种不同网格密度下的液滴长度和液滴间距几乎没有差别，两种网格密度下的液滴无量纲长度最大误差仅为 $4.7\% $ 。为了节省计算时间，本文采用 $6 \times 10^4$ 网格模型来进行模拟计算。

3.2 T型微通道内液滴生成尺寸与影响因素

本文中T型微通道内流动数值模拟时采用液液两相流，并通过连续相对分散相的剪切来产生液滴，主要研究连续相毛细数 $Ca$ 和两相流量比 $Q$ 两个无量纲数对分散相液滴长度及液滴生成频率的影响。 $Ca$ 与黏性和表面张力的比值成正比，它反映了黏性力和表面张力对流体流动的影响。 $Ca$ 、 $Q$ 定义分别为

$\quad\quad Ca =\frac{{U_{\rm c}}{\nu _{\rm c}}}{\sigma }$

(12)

$\quad\quad Q = \frac{{U_{\rm f}}d}{{U_{\rm c}}d} =\frac {U_{\rm f}}{{U_{\rm c}}}$

(13)

式中： $U$ 为宏观速度； $\nu $ 为动力黏度；下标“ ${\rm{c}}$ ”和“ $\rm f$ ”分别表示“连续相”和“分散相”。

3.2.1 毛细数对液滴生成的影响

保持其他物性参数不变，通过改变连续相初始速度改变连续相毛细数。在液滴形成初期，随着分散相体积的增大，液滴头部增大并逐渐完全阻塞主通道。连续相流体对分散相流体的挤压作用使得液滴颈部区域内分散相压力升高，最终在T型通道拐角处挤断分散相使其生成液滴。断裂后的分散相受两相表面张力的影响逐渐收缩，宽度最终几乎与通道宽度相同。连续相对分散相的挤压力和表面张力是液滴形成的主要因素，黏性力的作用可忽略不计。图4为连续相毛细数对液滴尺寸的影响。随着毛细数不断增大，连续相对分散相的挤压作用逐渐增大，使得连续相从注入主通道到被夹断生成液滴的时间逐渐变短，因此液滴长度也会随之不断变小，液滴的生成频率不断升高。

图5（a）为魏丽娟等^［20］采用实验方法得出的液滴无量纲长度与毛细数关系。根据前人研究^{［21−22］}，液滴无量纲长度与毛细数之间存在幂率关系，其一般通式为 $L/d = AC{a^B}$ ，式中A、B均为拟合参数。

本文所得液滴无量纲长度与毛细数的定量关系如图5（b）所示，由图中可见，在不同的分散相流量 $Q_{\rm f} $ 下，随着毛细数的增大，液滴无量纲长度不断减小，其函数关系为

$\quad\quad\left\{ \!\!{\begin{array}{*{20}{l}} {L/d = {\rm{0}}{\rm{.121}}C{a^{ - {\rm{0}}{\rm{.681}}}},\begin{array}{*{20}{l}} {}&{} \end{array}{Q_{\rm f}} = 0.050} \\ {L/d = {\rm{0}}{\rm{.084}}C{a^{ - {\rm{0}}{\rm{.861}}}},\begin{array}{*{20}{l}} {}&{} \end{array}{Q_{\rm f}} = 0.075} \end{array}} \!\!\!\!\right.\!\!\!\!\!$

(14)

在不同的分散相流量下，毛细数对液滴的无量纲生成频率的影响如图6所示。由图中可见，液滴的无量纲生成频率 ${f_0}$ 随着连续相毛细数 $Ca$ 的增大而逐渐增大。

3.2.2 流量比对液滴生成的影响

保持其他参数不变，改变分散相入口速度，以探究液滴无量纲长度 $L/d$ 与流量比之间的关系。图7为 $Ca = 0.03{\rm{2}}$ 下流量比与液滴长度的关系。从图中可以看出，随着流量比的增大，单位时间内流入主通道内的分散相不断增加，并且此时连续相对分散相的挤压作用力还不足以“剪断”分散相，因此最终液滴长度也会随之不断增大。

本文模拟了不同连续相毛细数工况下流量比对液滴长度和液滴生成频率的影响，结果如图8所示。从图中可以看出，液滴无量纲长度与流量比之间有明显的线性增长关系，且在不同毛细数工况下，变化趋势基本一致。两者拟合式为

$\quad\quad\left\{ {\begin{array}{*{20}{l}} {L/d = 0.815 + 2.899Q,\begin{array}{*{20}{l}} {}&{} \end{array}Ca = 0.048} \\ {L/d = 0.742 + 2.174Q,\begin{array}{*{20}{l}} {}&{} \end{array}Ca = 0.064} \end{array}}\!\!\!\! \right.\!\!\!\!\!\!\!$

(15)

根据文献［23］的研究，流量比Q与液滴无量纲长度之间的一般通式为 $L/d = M + NQ$ ，其中M和N为拟合参数，与本文模拟结果一致。在不同Ca下，液滴无量纲生成频率随着流量比的增加逐渐增加，但增加的趋势逐渐变缓，两个流量比变化趋势大致相同。

3.3 微流控梯状网络内弹状液滴行为 3.3.1 几何模型

图9为二维微流控梯状网络的几何模型，图中：L₄为水平长度；L₅为两平行通道间距；W₁为通道宽度；W₂为旁通宽度；Q₀为进口流量；p为出口压强。该模型为平行双通道，L₄为1 100个格子单位，L₅为60个格子单位，W₁为20个格子单位。两通道间由若干个相互平行的旁通管相连。为防止液滴流入旁通管内，W₂为15个格子单位，两相流体的表面张力和黏度为定值，流动过程中两通道Q₀均保持恒定不变，即Ca恒定。两个尺寸相同的弹状液滴在主通道内的初始相位长度为X_i，经过梯状旁通管后，最终相位长度为X_o。本文采用X_o/X_i表示液滴同步程度，T为无量纲时间，定义为 $T = {t_{\rm s}}/{t_0}$ ， ${t_0}$ 为固定时间步长， ${t_{\rm s}}$ 为实际时间步，本文取 ${t_0} = 10\;000$ 。

微通道内液滴的生成尺寸受连续相毛细数和两相流量等因素的影响。弹状液滴因尺寸大小合适，且纵向尺寸受限于微通道宽度，从而广泛用于各种微通道内两相流的研究。但弹状液滴只能在一定范围的连续相毛细数和两相流量比等条件的控制下产生，所以此处通过控制一定的变量参数来生成尺寸合适的弹状液滴，研究初始相位不同的弹状液滴在微流控梯状网络中的同步情况，并定量分析旁通管数目n、主通道内Ca和初始相位X_i对液滴同步程度X_o/X_i^［24］的影响。

3.3.2 旁通管数目对液滴相位同步的影响

控制Ca = 0.021、X_i 为 90个格子单位，分别模拟了n = 1～17时液滴同步情况。图10为n = 10的结构中不同无量纲时刻的液滴同步情况，可以看出，液滴经过梯状网络后，液滴相位有所缩短，但缩短长度不明显。

图11为旁通管数目n对液滴同步程度的影响。从图中可看出，随着n的增加，液滴同步程度X_o/X_i逐渐变小，最终趋近于0.8。该工况下两液滴没有实现完全同步，即使不断增加旁通管数目，也不会出现液滴完全同步的情况，X_o/X_i最终会无限趋近于0.8。综上所述，增加旁通管数目可以促使液滴同步，但对其影响有限。

3.3.3 液滴初始相位长度对液滴相位同步的影响

控制Ca和n恒定，研究X_i对液滴同步程度的影响，结果如图12所示。在Ca = 0.005 3、n = 16时，X_o/X_i随着X_i增加而变大，同步效果变差。当X_i大于140个格子单位时，X_o/X_i开始在0.9附近波动。保持n不变，Ca = 0.002 7，此时不同初始相位长度的液滴均实现相位同步。

图13为Ca = 0.002 7和n = 16时不同初始相位长度液滴同步的无量纲相位时间。从图中可以看出，随着液滴初始相位长度的增大，液滴同步的无量纲相位时间T不断增长，两者呈现明显的线性关系。

3.3.4 初始进口速度对液滴相位同步的影响

控制X_i和n恒定，研究连续相毛细数对液滴相位同步的影响。图14为Ca = 0.005 3和n = 16时不同无量纲时刻液滴的相位同步状况。从图中可以看出，液滴经过梯状网络后，液滴的相位逐渐实现同步。

图15为连续相毛细数对液滴相位同步程度的影响。从图中可以看出，Ca = 0.002 6～0.005 3内X_o/X_i越来越大，即相位同步效果越来越差，在 $Ca{\text{≈}}0.005\;3 $ 时X_o/X_i达到峰值。随后，随着Ca的增大同步效果逐渐好转，在Ca = 0.53时又基本实现液滴相位的同步。由此可见，Ca对液滴相位同步过程有着显著影响。

4 结　论

基于LBM颜色梯度模型对液滴在T型微通道的生成过程和初始相位长度不同的弹状液滴在微流控梯状网络内的相位同步过程进行了数值模拟，主要结论有：

（1）液滴无量纲长度随着连续相毛细数增大而变短，两者存在明显的幂率关系，且液滴无量纲生成频率与连续相毛细数呈正相关。这与文献[20]中的实验结果吻合较好。

（2）两相流量比与液滴无量纲长度呈线性关系，液滴的无量纲生成频率随着两相流量比增大而加快，但上升趋势逐渐变缓。

（3）在梯状网络结构中，增加旁通管数目和缩短液滴初始相位长度在一定程度上可以改善液滴相位同步程度，但不能实现液滴相位完全同步；连续相毛细数对液滴相位同步影响明显，且在连续相毛细数小于0.002 7和大于0.005 3时，基本可以实现液滴相位的同步。