微流控技术依托MEMS(micro-electro-mechanical system)工艺可以制备出10~500 μm宽度的微沟道, 用于生物细胞、蛋白质、DNA或病毒的输运、分离、富集和提纯^[1-2].这种微型化的生物处理器件可极大提高精准识别和快速诊断生物样本的能力.因此, 针对生物芯片实验室(lab-on-a-chip)或微全分析系统(micrometer-scale total analysis systems, μTASs)的现代生化检测技术由于微流控技术的发展得到了越来越多学者的广泛关注和推广.在微流控技术多元化的研究过程中, 液滴微流控(droplet-based microfluidics)是其发展和创新的一个研究分支.液滴微流控芯片可连续性地产生皮升(pL)到纳升(nL)量级的液珠, 用于单细胞培育^[3]、基因识别^[4]和药性毒理研究等领域^[5].以单细胞基因测序为例, 要进行单个活体细胞的基因测序, 必须进行目标细胞和免疫磁珠的空间孤立, 防止其他细胞的特异性干扰, 阻碍基因标记物的特征提取.因此, 液珠包裹和封装的生物细胞、蛋白质、病毒和反应物可以有效降低外界对生物样本的生化污染和信息干扰.每个液珠作为一个微型反应容器, 为个体生物样本提供很好的研究平台.由此可以极大提高子代细胞的基因序列识别, 减少种群中光学、电学、声学或应力等造成的宏观信息干扰, 使分子生物学更能目标化地测量单个生物的个体分化、生殖和基因突变^[6].

当前, 液滴微流控芯片结构主要分为正交型^[7]和同向型^[8].不相溶两相流体的黏滞系数、流速、接触角度, 以及毛细数(capillary number, C_a)等多维变量共同决定了液滴大小和几何均匀性.比如：增加连续相的液体速度, 易导致离散相液体没有足够时间在沟道伸展形成均匀球形液滴, 实验观测为喷射(jetting)现象.反之, 离散相流体速度高于连续相流体, 液滴容易在两相交界面发生冻结(squeezing)现象.因此在实验前期发展一种理论模型显得尤为必要.合理的理论模型可以引导实验开展, 并且丰富实验内涵.已经有大量学者开展了T形液滴微流控相关建模和仿真^[9], 得到与实验相符的理论解析集.国内方面, 已有大量的学者参与了液滴形状的关键参数研究, 并提出了相关参数的优化结果^[10-12].然而在液滴形成与细胞封装的理论研究上, 鲜有学者进行相关报道.尤其在细胞随液体流动的轨迹和运动特征方面, 缺少相关的动态仿真数据.因此发展一种水平集与动力学结合的数值模型, 用于T形微流控细胞封装具有一定的研究价值.

1 理论模型

当前的液滴微流控模拟方法主要分为VOF法(volume-of-fluid method)、格子玻尔兹曼法(lattice Boltzmann method)和水平集方法(level set method, LS).就前两种方法而言, VOF法的计算精度不高, 不能平滑得到离散相液滴平均曲率边界, 容易造成表面张力误差增大^[9].格子玻尔兹曼方法对离散相和连续相的界面张力很难找到合理的控制参数, 用以适应液滴滑移带来的收敛问题.由于LS方法满足质量守恒约束, 并且具有较好的计算精度, 因此本文选择该方法分析T形微沟道的液滴成形和细胞封装.T形微流控芯片的几何结构如图 1所示.其中：L和H分别表示连续相液体流动长度和T形微槽长度; W_d和W_c分别表示离散相和连续相的微槽宽度; U_d和U_c分别是离散相和连续相的流体速度.两相流体的分布和流向服从Navier-Stokes方程, 如式(1)所示.

(1)

式中：u为液体速度(m/s); t为时间(s); ρ为液体密度(kg/m³); η为液体黏滞系数(Pa·s); p为压强(Pa); F_sf是两相流之间表面张力(N/m³), 其表达式如式(2)所示：

(2)

式中, σ, δ, κ分别为液体之间的表面张力系数(N/m)、两相流界面函数、两相流界面的曲率函数.它们的函数表达式如式(3)和式(4)所示：

(3)

(4)

其中, n为单位向量, 与液滴界面的法线方向一致, 其表达式如式(5)所示：

(5)

式中, Φ为相函数, 取值范围在[0, 1]之间, 代表了两种不相溶液体的浓度变化关系, 满足式(6):

(6)

式中, γ和ε分别为网格重置参数(m/s)和界面厚度控制参数(m).两相流密度与时间、流速的函数关系, 如式(7)所示：

(7)

根据式(8), 流体的密度和黏滞系数可与势函数Φ联立, 求解方程(1), (6)和(7)即可获得不同时间和位置对应的流体浓度变化.

(8)

如图 1所示, 在Φ等于0.5时, 其闭合曲线为离散相液滴与连续相流体的交界面.因此通过闭合曲线积分(如式(9)所示)可得到相应液滴的有效尺寸d(m),

(9)

在两相流中, 毛细数C_a反映了表面张力对液体流动的影响.它与黏性力和表面张力的比值有关, 如式(10)所示：

(10)

不失一般性, 本文设定的两相流流体中, 连续相势函数Φ=0, 离散相势函数Φ=1.此外液滴的形状还与表面所受应力有关, 具体为流体与湿壁剪切摩擦F_fr有关(F_fr=-ηu/β), β为滑移长度(m), 与离散网格的具体尺寸有关.当液滴在微槽移动过程中, 连续相、离散相与湿壁三相交点处所作的液-液界面的切线与固体交界面之间的夹角称为接触角θ(rad), 如图 2所示.求解上述方程可得相关液滴的形成.

本文主要以红细胞为封装的目标粒子, 其受力满足Stokes应力, 如式(11)所示：

(11)

式中：m_p和v分别为细胞质量和运动速度; τ_p等于ρ_pd_p²/18η(ρ_p和d_p分别为细胞密度和直径).细胞与墙壁保持弹性碰撞, 满足能量守恒.同时细胞密度接近培养溶液密度, 并且尺寸为微米量级.因此忽略相关重力、浮力和布朗运动力等因素干扰, 能节约计算时间, 提高仿真效率.

2 水平集仿真与结果

根据上述理论模型的简要说明, 本文将油/水(O/W)两相流作为仿真目标, 其中连续相为正十二烷(n-dodecane oil), 离散相为水溶液(water).它们的基本物理量如表 1所示.为更好研究毛细数的影响, 离散相水的流速U_d=0.012 m/s, 且为常数.因此主要改变连续相油的流速, 其初始数值为U_c=0.037 m/s.因此毛细数C_a等于0.01;W_d和W_c分别为50 μm和80 μm; L和H分别为700 μm和100 μm.油水相的表面张力σ为5.0 mN/m, 且接触角θ为165°.仿真平台采用COMSOL有限元软件, 其提供的Laminar Two Phase Flow, Level Set (tpf)和Particle Tracing for Fluid Flow (fpt)可以很好地进行物理耦合, 实现上述偏微分方程的动态计算.求解结果如图 3所示.

对T形微流控通道而言, 液滴成形主要与离散相和连续相的流体运动速度有关.改变连续相的流体速度U_c, 研究不同C_a数值对液滴大小的影响, 如图 4所示.

随着毛细数C_a逐渐增大, 固定时间内, T形微流控通道产生液滴的频率逐渐变高.在0~30 ms时间内, C_a等于0.008或0.01时仅产生一个液滴, 但是C_a为0.04, 0.08和0.1时分别产生4, 7和9个液滴, 该仿真结果也与文献[9, 11-12]一致.通过上述的液滴形态和流体分布进行红细胞封装的理论研究, 本文仿真的红细胞ρ_p和d_p分别为1 050 kg/m³和5 μm.忽略细胞运动对流体速度场的影响, 根据式(11), 结合牛顿经典应力公式, 即可得到红细胞运动轨迹, 如图 5所示.

红细胞在微流通道内首先受鞘流作用, 使其保持在中线位置随流体运动.图 5的红细胞自上向下运动到两相流的交界面, 受流体黏滞力作用, 被封装在了液滴内部, 跟随液滴向输出通道方向流动.本文根据上述的C_a数值进一步将细胞动力模型与流体运动模块相互耦合, 在30 ms时间内得到了不同连续相流速场内部的细胞的运动轨迹如图 6所示.该结果暗示液滴直径越大, 相应封装细胞的能力越强, 细胞在液滴内部也存在一定的抖动.反之, 较小的液滴尺寸, 高的输出频率使得细胞的包裹能力相对较弱, 细胞极易从离散相进入连续相中单独运动.细胞在进入两相流交界面位置时, 由于黏滞力大于界面张力, 容易驱使细胞贴近微槽边缘移动, 由此降低了细胞封装效率.低的C_a数值虽然能够保证细胞的自动封装, 但其频率较低.单个细胞与细胞之间需要保持足够距离, 否则单个液滴将会封装多个细胞.倘若开展细胞内的信使RNA测序实验, 多个细胞与单个编码磁珠配对容易造成基因误码率增大.因此, C_a数值与细胞间距对于封装的研究是一个多参量优化问题.除此之外, 液滴的接触角、细胞与墙壁非弹性碰撞都需要进一步的研究学习.通过对细胞或粒子进行荧光染色, 实验观测可采用高速CCD摄像机连续抓拍液滴内部的微粒, 掌握其相应的运动轨迹.

3 结论

1) 利用水平集方法仿真T形微流控通道的液滴成形和尺寸参数, 结果表明：离散相速度固定在0.012 m/s时, 毛细数C_a越大, 连续相速度越大会使液滴尺寸越小, 输出频率越高.

2) 通过Stokes动力学模型耦合水平集方法, 得到红细胞在两相流的运动轨迹.毛细数C_a越大, 黏滞力越大, 细胞越容易紧贴墙壁呈直线运动.但是封装效率迅速降低.若黏滞系数过小, 细胞容易在液滴内部发生抖动, 影响细胞与液滴的同步输运, 该问题需要多目标优化共同实现.