1 前言

毛细管液相色谱是一种高效的色谱分离技术，具有固定相用量少，流动相消耗少，样品消耗少，环境污染少等优点^[1]。毛细管液相色谱的核心部件是毛细管柱，其管径通常为毫米级与微米级，由于表面张力与壁面粘附力占主导地位，毛细管中流体特性与大尺寸管道中流体具有非常大的区别^[2]。

弹状流是毛细管中常见两相流流型之一，弹状流的特点是：某一相均匀子弹形状的液滴或气泡被另一相分隔，液滴或气泡与管壁之间有薄薄的液膜^[3]。气泡的形状、气泡速度、液弹长度及液膜厚度是非常重要的流体力学参数^[4]。毛细管中弹状流的压降、传质与液弹长度有关^[5]，且液膜厚度也会影响传热速率与传质速率，因为液膜的存在增加了界面面积，因而在传质过程中起决定性作用^[6]。然而不能像其他参数一样(如表观速度和流体的物理性质)预先确定液膜厚度、液弹的长度及液弹的速度等重要的流体力学参数。对于微通道中两相流的流体力学已有很多实验研究^[7~11]。目前虽有一些气液弹状流的实验与模拟研究^{[5, 12~17]}和液液弹状流的实验研究^[18~21]，大部分研究液液弹状流的压降^[18]、内部循环^[19]与传质^[21]问题，但对于微通道中液液弹状流的流体力学的分析还不是很清楚。且大多数文献研究的是气液弹状流，由于分散相的黏度，液液弹状流的界面剪切力显著，因此，液液弹状流的液膜厚度、液弹速度与液弹长度可能不同于气液弹状流^[8]。所以有必要系统地研究液液弹状流的流体力学特性。本文以T型毛细管为模型，使用计算流体力学的方法研究毛细管中液液弹状流的流体力学特性。研究不同操作条件下液弹长度、液膜厚度与液弹速度。模拟结果用Dore等^[19]的实验结果进行验证。

2 数值模型与算法 2.1 控制方程

使用计算流体学(CFD)Fluent(16.2)软件进行数值模拟，选择由Hirt和Nichols^[22]提出的VOF多相流模型追踪不相溶两液体的界面。VOF多相流模型的控制方程如下：

连续性方程：

$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \vec v} \right) = 0 $

(1)

动量方程：

$ \frac{{\partial \left( {\rho \overrightarrow v } \right)}}{{\partial t}} + \nabla \cdot \left( {\rho \overrightarrow v \overrightarrow v } \right) =- \nabla P + \nabla \cdot \left[{\mu \left( {\nabla \overrightarrow v + {{\left( {\nabla \overrightarrow v } \right)}^T}} \right)} \right] + \rho \overrightarrow g + \overrightarrow F $

(2)

体积分数方程：

$ \frac{{\partial {\alpha _\text{d}}}}{{\partial t}} + \overrightarrow v \cdot \nabla {\alpha _\rm{d}} = 0 $

(3)

流体物理性质方程：

$ \rho = {\alpha _{\rm{c}}}{\rho _{\rm{c}}} + {\alpha _{\rm{d}}}{\rho _{\rm{d}}}\;\;\;\;\;\;\;\;\mu = {\alpha _{\rm{c}}}{\mu _{\rm{c}}} + {\alpha _{\rm{d}}}{\mu _{\rm{d}}} $

(4)

毛细管中流体流动为层流，界面张力与壁面吸附力占主导作用，Fluent中界面张力模型是由Brackbill等^[23]提出的连续表面力(CSF)模型，VOF模型中附加的界面张力导致了动量方程中的源项。因此，动量方程中的源项为：

$ \bar F = \sigma \left[ {\frac{{\rho \kappa \hat n}}{{\left( {1/2} \right)\left( {{\rho _{\text{c}}} + {\rho _{\text{d}}}} \right)}}} \right] $

(5)

$ \hat n$为界面法线，κ为局部界面曲率，计算公式为：

$ \hat n = \frac{n}{{\left| n \right|}}\;\;\;\;\;\;\;n = \nabla {\alpha _{\rm{d}}} $

(6)

$ \kappa = \nabla \cdot \hat n $

(7)

2.2 物理模型

以往研究中均有使用二维与三维模型，当管道横截面为矩形时通常使用三维模型研究管道中两相流的流动情况^{[24, 25]}。且Qian等^[5]和Cherlo等^[25]的研究均发现分散相弹长在二维和三维模型中没有很大差异，但三维模型中的压降大约是二维模型的两倍。本文主要讨论圆形管道中液弹长度，液膜厚度与液弹速度，所以为了节省运算时间，采用二维模型。使用GAMBIT软件构造几何模型并进行网格划分。图 1为使用T型管的几何模型简图。

2.3 网格及边界条件

两相入口为速度入口边界条件，出口处设为压力出口边界^[26]且为标准大气压。壁面为无滑移边界条件，壁面接触角为180°。分散相入口α_d=1，即该处全为分散相，连续相入口α_c=1即该处全为连续相。连续相先充满管道，当t=0时，两相液体同时流入管道。所用流体的物理性质如表 1。

采用四种均匀正方形网格划分方式，边长分别为20、12.5、10、8 μm，对应网格数分别为4100、10496、16400、25650.

2.4 数值方法

使用基于压力的非稳态求解器，不考虑重力作用的影响。连续相为[C4mim][NTf2]离子液体，分散相为水。空间离散化梯度选择更准确的基于节点的格林-高斯方法^[16]。压力项关联采用交错压力格式(PRESTO!)，流场采用压力与速度耦合的PISO方法，动量方程选择二阶迎风格式。时间步长取1×10^-5 s, 为了保证计算收敛，连续性方程的残差绝对收敛标准在10^-5以内，速度残差绝对收敛标准为10^-4。

3 结果与讨论

T型管道直径d=0.2 mm，为了形成液弹且使流动充分发展，管长与直径比值的最小值应为40^[27]。因此，管道长度设为40 d。分散相流量从0.13~1.2 mL·h^-1变化，连续相流量从0.2~0.8 mL·h^-1变化，根据Tsaoulidis等^[28]的两相流型图可知，流动在弹状流范围内。

3.1 网格独立性检查

图 2表示不同网格数时分散相体积分数图。由图 2可知，仅在界面处发生显著的差异，分辨率越高，界面越清晰。发现网格数为16400时交界面分布良好，过渡自然，没有明显差异。网格数为4100、10496、16400、25650(Q_d=Q_c=0.8 mL·h^-1)对应的液弹速度分别为0.0153、0.0159、0.0159、0.0159 m·s^-1，对应的液弹长度分别为0.463、0.467、0.490、0.484 mm。液弹长度与液弹速度几乎与网格分辨率无关。本文使用的网格数为16400，网格划分如图 3所示。

3.2 液弹长度分析 3.2.1 毛细管数的影响

Garstecki等^[29]发现液弹的形成有“挤压”和“剪切”两种机理，两种机理的临界毛细管数(Ca)约为10^-2，本文Ca包含临界Ca，即低毛细管数时为挤压机理，刚产生的液弹进入主通道且不断增长直到几乎占满通道横截面，连续相的压力在液弹上游积聚，使得液弹在下游延长。当毛细管数增大时为剪切机理，分散相的剪切力也促使了液弹的断裂。

图 4为三种不同连续相流量分数下无量纲液弹长度与毛细管数(Ca)的关系，由图可知低速时，无量纲液弹长度随Ca的增大略有不同，但随Ca进一步增大趋于减小。且无量纲液弹长度随连续相流量分数的增大而减小。模拟结果与Dore等^[19]的实验结果总体一致。

3.2.2 分散相与连续相流率的影响

连续相流率不变，改变分散相流率，观察液弹长度变化。图 5为分散相流率对液弹长度的影响，当连续相流率为常数时，液弹长度随分散相流率的增加而变长。且由3.2.1节知，液弹长度随连续相流率增加而减小。

3.2.3 关系式

关于弹长的关系式大部分来源于气液弹状流^{[5, 17, 30]}。有些关于液液弹状流液弹长度^{[20, 29, 31]}的研究。Tsaoulidis等^[20]提出的液弹长度关系式只适用于较大毛细管数(0.03~0.18)，该毛细管数远远大于本文所用毛细管数。实际上，液弹长度受两方面影响：一是受两相流率比Q_d/Q_c影响的界面断裂，二是受连续相毛细管数影响的界面张力与连续相剪切力之间的平衡^[31]。因此，根据模拟结果提出液弹长度的估测式：

$ \frac{{{L_\rm{d}}}}{R} = \left( {1.263 + 1.981 \times \frac{{{Q_\rm{d}}}}{{{Q_\rm{c}}}}} \right)C{a_\rm{c}}^{-0.0891} $

(8)

为了验证模拟结果的正确性，将式(8)的计算结果与Dore等^[19]的实验值进行比较如图 6, 计算值与实验值的平均误差在20%左右。偏差较大的原因可能是实验中接触角不明确，且模拟中假设壁面光滑无滑移，与实验条件有一定的差异所致。

3.3 液膜厚度

如图 7所示，在液液弹状流中，分散相液弹几乎占据了整个毛细管的横截面积，这使连续相形成薄液膜将分散相液弹与管壁分隔。分散相前端和后端为半球形，液膜沿分散相液弹长度分布不均匀，分为均匀液膜区和过渡区，本文研究均匀液膜区液膜厚度。采用FLUENT后处理软件，选取分散相体积分数为0.1的界面，即该界面处几乎全部为水相，计算液膜厚度。

Bretherton^[13]发现液膜厚度与毛细管数(Ca)有关，且随毛细管数增大而变厚。图 8为模拟的无量纲液膜厚度与Ca的关系，且模拟数值与文献中的液膜计算公式进行比较。可以发现模拟液膜厚度与Dore等^[19]的模型吻合很好，平均误差约为11%。但Fairbrother等^[14]、Bretherton^[13]、Aussillous等^[15]和Irandoust等^[12]模型计算的液膜厚度远远小于模拟数据。这是可以理解的，因为这些模型是在气液弹状流的基础上总结的，气液弹状流与液液弹状流主要差异是两相黏度差，在气液弹状流中，分散相气体的黏度可以忽略，因此对界面的剪切力也可以忽略。然而，在液液弹状流中两相黏度差减小意味着两相界面剪切力的增大，这导致液膜区不可忽略的流动，因此液弹速度、液膜厚度与气液弹状流不同^[32]。

3.4 液弹速度

由于液弹速度与停留时间有关，所以它是传质过程中一个重要参数。毛细管中的弹状流，由于连续相液膜的存在，液弹或气泡并不能完全占据管道横截面积，这使得液弹速度总比混合速度大^[20]。图 9为模拟液弹速度与混合速度关系图。由图可知模拟液弹速度与Dore等^[19]的实验值整体上吻合很好。

本文模拟了不同混合速度下的液弹速度，因此，液弹的速度可由下式进行估算：

$ {U_{\rm{d}}} = 1.157{U_{\operatorname{mix} }} $

(9)

为了验证该式的准确性，将式(9)的计算值与Dore等^[19]的实验值对比如图 10，实验值与式(9)计算值整体吻合很好，平均误差在17.7%左右。且液弹速度与混合速度差随混合速度的变大而增大。

4 结论

使用计算流体力学(CFD)的方法，采用VOF多相流模型，使用水和[C4mim][NTf2]离子液体两相流体研究毛细管中液液弹状流的流体力学特性，得出结论如下：

(1) 液弹长度与两相流率和毛细管数有关，液弹长度随分散相流率的增大而增大，随连续相流率的增大而减小。不同分散相流量分数下，低速时，无量纲液弹长度随Ca的增大略有不同，但随Ca进一步增大趋于减小。

(2) 液膜厚度受毛细管数的影响且随毛细管数的增大而变厚，液膜厚度与现有的关联式符合较好，平均误差大约为11%。

(3) 由于液膜的存在，液弹速度高于混合速度，且液弹速度与混合速度的差异随混合速度的增大而增大。根据模拟数值提出关于计算液弹长度与液弹速度的关联式与实验值吻合性较好，平均误差分别约为20%，17.7%。

符号说明：

Ca	—毛细管数(${\mu _{\text{c}}}{U_{\operatorname{mix} }}/\sigma $)	Ud	—液弹速度，m·s-1
Cac	—连续相毛细管数(${\mu _{\rm{c}}}{U_{\rm{c}}}/\sigma $)	Umix	—混合速度，m·s-1
d	—毛细管直径，m	α	—体积分数
F	—体积力，N	εc	—连续相流量分数(${Q_{\rm{c}}}/\left( {{Q_{\rm{c}}} + {Q_{\rm{d}}}} \right) $)
g	—重力加速度，m·s-2	κ	—曲率，m-1
h	—液膜厚度，m	μ	—黏度，Pa·s
Ld	—液弹长度，m	ρ	—密度，Kg·m-3
$\hat n $	—界面法线	σ	—界面张力，N·m-1
P	—压力，Pa	下标
Q	—体积流率，mL·h-1	c	—连续相
R	—毛细管半径，m	d	—分散相
t	—时间，s	mix	—混合物