CFD数值模拟技术在液滴微流控多相流特性研究的应用进展

王炳捷，李辉，杨晓勇，白志山

（华东理工大学机械与动力工程学院，上海，200237）

过程强化是实现绿色工艺的关键技术，通过小尺寸和小型化设备的发展，降低单位能耗及副产品生成，并最终达到提高生产效率、降低生产成本、提升安全性和减少环境污染的目的。然而，小尺寸紧凑型设备的发展也给其中流体动力学特性和运动行为的研究带来了新的挑战，特别是多相流体系[1-3]。与传统的宏观大尺寸系统相比，借助微通道控制的两相流体具有比表面积大、传递距离短、混合速度快等优点，可以减少传质限制，以获得更好的性能[4]。因此，液液微流控体系在化学反应[5]、液液萃取[6]、生物分析[7]、结晶过程调控[8]以及结构性材料制备[9]等领域取得了重大进展。为了分析操作条件、流体性质以及通道结构对多相流的影响，当前已经从试验角度开展了许多有益的研究[2-3,10-13]。

根据以往的试验研究，微通道内两相流间的相互作用及动力学行为受到许多参数的影响，如两相流速、两相黏度、界面张力、微通道结构及壁面润湿特性等。同时，预测微通道装置中液滴尺寸及成型机制的理论分析仍然比较复杂，且准确性不足[14-15]。因此，仅靠试验研究或理论分析很难对微通道中两相流动现象和特性有全面的了解。与此同时，微流控系统中参数的测量和流场特性的获取也是极为困难的工作。由于这两个原因，目前微流控系统理论多为经验研究所得。

此外，液滴和气泡在复杂几何结构中的运动对于许多科学和工程应用而言非常重要。在过去的几十年里，无约束介质中气泡和液滴的运动引起了科学界的极大关注[16]。相反，关于气泡/液滴在复杂几何结构通道内演变的研究内容非常有限。一方面微通道壁面可以压缩气泡和液滴，在某种程度上，液滴/气泡的形状由微通道结构所决定。另一方面，气泡/液滴可以在无界介质中自由演化，但在微通道结构的限制下运动。为了解决这两类问题，试验过程通常将基本现象分离出来，从宏观上研究第二相的运动。而液滴或气泡在复杂几何中演化的问题则通过简化控制方程的方法来获得一些近似的解析解[17]。因此，开发能够实现微通道内液滴成型、液滴聚并/融合、液滴溶解、颗粒聚焦等复杂传质传热过程的数值模拟方法极为必要。

与液滴微流控技术相关的现有文章的综述重心主要集中在微通道结构和微流控设备的发展[18-19]、微流控技术在获得结构性/功能性材料的应用[20-21]以及微流控技术及其产物在不同领域的拓展[22-23]，却很少把目光集中于数值模拟方法对于液滴微流控多相流特性研究的重要辅助作用。本文将从数值模拟所涉及的液滴微流控装置结构及演变（第2 节）、液滴微流控模拟方法及优化（第3 节）、微通道内两相流作用过程及原理（第4节）三个方面，系统综述当前数值模拟方法在液滴微流控技术及应用方面的研究进展（图1所示），并在第5节针对当前数值模拟方法所存在的瓶颈进行重点阐述，展望未来数值模拟技术在液滴微流控探索应用的发展方向。

1 微通道结构

微通道限制了微流体边界，因此通道尺寸结构也会对流体特性产生影响。利用黏性剪切力来产生液滴的基础微通道结构可分为：同轴流（co-flow）、 交 错 流（cross-flow） 和 聚 焦 流（focusing-flow）。微流控技术的发展也使得微通道结构以上述三种结构为基础，不断向复杂组合或优化结构拓展和演变。而针对不同的研究对象（液滴成型、液滴聚并、气泡溶解等）所采用的微通道结构尺寸和模拟方法以及所考虑的影响因素也不尽相同（见表1和表2）。

1.1 同轴流结构

1.1.1 基础同轴流结构

用于液-液介质研究的同轴流微通道最早由Umbanhowar 等[53]提出，其结构特征在于分散相流体与连续相流体在平行流中相遇。同轴流结构可以是准二维（2D）平面的或三维（3D）同轴的。其中，前者通常采用多板层组合获取封闭的同轴流通道结构，而后者则采用不同尺寸的圆形管（或方形管/槽）进行同轴内配合组装获得。Lan 等[27]借助CFD 数值模拟方法研究了准二维（2D）平面同轴流结构中的液滴成型过程。通过分析液滴分离位置附近的局部压力和流速，揭示了滴流成型和射流成型的本质区别［图2(a)］。Xiong等[26]结合micro-PIV测试和数值模拟发现了液滴形成过程中的内部循环流，并探究了两相流速、两相黏度和界面张力对液滴内部循环流特性的影响［图2(b)］。而Deng等[45]、Guillaument 等[47]和Nacimirad 等[24]则探究了3D 同轴流微流控结构中两相流速、两相黏度、界面张力以及微通道壁面润湿特性等参数对液滴成型方式和尺寸的影响规律。

1.1.2 优化同轴流结构

在传统同轴流微通道结构的基础上，通过对分散相或者连续相通道结构优化，可改变微通道内部流体的流动特性。Heuberger 等[31]通过在分散相通道尖端开设方形槽道或布置特定结构的开孔挡板，构建了独特的喷射式同轴流液滴成型结构，其本身具有双相流体振荡器特征，可在高速率和高长径比条件下连续产生单分散性液滴。随后的数值模拟帮助确定了上游质量和动量传递在激发液-液界面振荡中的重要作用，而试验结果则进一步证明了该种优化结构在提升液滴成型的稳定性和产量方面的潜力。而Wang[46]和Chekifi[32]等则改变两相流体接触段的几何结构，构建了锥形同轴流结构微流控系统，并利用数值模拟方法研究了流动参数和锥形结构对液滴成型过程中液滴尺寸和生成频率的影响。如图3所示。

图1 数值模拟方法在液滴微流控技术中的应用

图2 基础同轴流微通道结构中液滴成型过程

1.2 交错流结构

1.2.1 基础交错流结构

交错流结构是指分散相流体和连续相流体以一定角度θ（0°＜θ≤180°）在微通道中相遇[54]。在传统交错流结构中，分散相和连续相流体通常是在T形结构通道内以正交方式相遇的。由于T形通道兼具结构简单性及产物的高单分散性（粒径偏差CV＜2%）特征，从而被广泛应用于液滴/气泡的生成[55]。Yeh等[51]模拟研究了交错流微通道结构中的多相流流动特性，揭示了单分散海藻酸钙颗粒形成过程及尺寸调控方法。而Padoin等[43]借助数值模拟方法评估亲疏水壁面润湿特性变化对微通道中气液流动模式的影响。模拟结果表明，表面张力主导了微通道内部气液流动特性，具体表现为在相似的边界条件下，微通道壁面亲疏水特性的改变（由亲水壁面θ=25°变为疏水壁面θ=105°），造成其中气液流动模式的改变（由泰勒流变为分层流）。如图4所示。

表1 模拟方法及对象

1.2.2 优化交错流结构

微通道的几何形状对微通道的动态特性和总压降有显著影响，通过改变主通道结构会对其中液滴/气泡的成型产生显著的影响。Chandra 等[35]通过在两相流体交汇处下游主通道内侧等间距布置不同高度的三角形阻塞单元，首先构建了优化交错流结构微流控模型。随后，采用数值模拟方法分别研究了优化结构对液滴成型、尺寸、形貌和压降的影响，证明了优化结构对流体传质传热的有益作用。研究结果表明，在所形成的液滴内部及界面处的混合效果显著提升。此外，三角形阻塞单元间距的变化也会影响液滴成型、液滴长度和压降。而Zhu等[36]更是通过改变主通道横截面形状（梯形、三角形、圆底长方形、半圆形等），构建了优化交错流结构微流控模型，探究了不同结构对液滴分离时间和尺寸形貌的作用关系。结果表明，对于不同界面的微通道，液滴分离时间、分离直径和去除时间服从以下规律：三角形截面＜梯形截面＜底部弯曲矩形截面＜矩形截面＜倒立梯形截面。此外，已有报道在微流控装置周边布置特殊场（离心场、磁场、电场等），借助外加场作用克服流体原本的物理化学性质（黏性力和表面张力），强化和调控微通道内部流体的流动形式，进而实现流体的输送、混合等过程。Ren 等[56]模拟研究了离心场作用下并联交错流微通道中液滴成型和融合过程，揭示了离心场转速与液滴尺寸、液滴产率、融合速率之间的构效关系，获得了液滴成型方式转变的调控机制。如图5所示。

表2 微通道结构与特征尺寸

1.2.3 交错流结构阵列

实际应用中，单通道交错流微流控结构不能满足液滴高通量可控合成的要求。考虑到微通道并联不存在放大效应，因此可以通过并行交错流结构中的分散相通道，共用连续相通道，构成交错流结构阵列。在这种结构中，分散相流体被迫同时进入多个交错流结构微通道，然后在两相交汇处被连续相流体截断，同时形成多个拟球型液滴。Kobayashi等[40]采用数值模拟方法探究了交错流结构阵列中通道尺寸和分散相流速对液滴成型过程的影响，并获取了该过程中动态界面张力的变化规律。研究表明，最终的液滴尺寸与通道尺寸成正比，而液滴生成速率与通道尺寸成反比。van Dijke等[42]在单分散油滴高通量合成过程的模拟研究中发现，液滴形成点之间的距离由分散相通道末端平台的高度和施加的压力决定，而与单分散液滴生成过程中的液滴或气泡大小无关。如图6所示。

图3 优化同轴流微通道结构中液滴成型过程数值模拟［31-32，46］

图4 基础交错流微通道结构中两相流体相互作用过程数值模拟［38，43，51］

图5 基础交错流微通道结构中两相流体相互作用过程数值模拟［35-36］

图6 交错流结构阵列中两相流体相互作用过程数值模拟［40，42］

1.3 聚焦流结构

在聚焦流结构微通道中，两相流体通过流体动力学进行聚焦，然后以延伸流状态通过收缩的微通道，继而获得尺寸更小的液滴。与同轴流结构微流控装置类似，聚焦流微流控装置也分为准二维（2D）平面和三维（3D）轴对称类型。Lan 等[28]在模拟聚焦流结构通道中液滴形成过程时发现，液滴形成过程的总传质比可达50%。在液滴形成过程的早期，由于连续相的剪切作用，出现了较强的内环流。这增强了两相之间的传质，并在液滴顶部形成了一个浓度缺口。在液滴形成过程的后期，液滴运动速度明显增加，液滴内部不再有循环流动，传质速度变慢。Dziubiński[52]通过建立准二维（2D）平面聚焦流结构微流控模型，探究了两侧连续相不对等流速条件下，聚焦流流型和位置的变化规律，提出了描述聚焦流形态和位置的简单理论模型。该模型可在误差高达25%的情况下确定聚焦流的位置和形状，有助于推动单颗粒流变学发展，使得借助水动力聚焦来直接控制颗粒受力和运动成为可能。Dang 等[50]在两相微通道夹角θ为30°的聚焦流微通道结构中，系统研究了通道壁面的接触角、表面张力和流体黏度对泰勒气泡长度、体积和形状的影响。研究发现，在最高表面张力和最高流体黏度的条件下，随着接触角的增大，气泡长度大幅度减小。这是气泡端部形状由凸形向凹形的变化和气泡周围液膜体积减小的共同作用结果。针对高黏度流体，液滴粒径控制和成型将更加困难。为解决这一问题，Li等[30]通过对聚焦流结构微流控装置外加电场作用，模拟研究了电场和电荷在流体界面上的相互作用特性，实现了电场力对高黏度液滴成型动力学的有效调控。如图7所示。

2 模拟方法

由于微尺度装置中存在着较大的比表面积，因此，流体界面在许多物理过程中扮演着重要的角色，例如液滴碰撞、合并和破裂[57]。在具有移动边界和界面流动问题的模拟过程中，存在基于固定网格的离散Navier-Stokes（NS）方程，并以显式或隐式来表达界面演化的两种方法，即界面追踪方法和界面捕捉方法。其中，移动边界的界面追踪基于计算网格对界面的显式描述，具有较高的数值计算精度，但其适用范围仅限于Stokes流。常用的前沿追踪法（front-tracking method）[58]属于界面追踪方法。与界面追踪方法不同，界面运动可以简单地通过界面捕捉方法中对应的相位函数来得到。界面捕获方法使用单独的相位函数，以隐式表示界面。用于解决不可压缩两相问题的常用界面捕捉方法包括水平集法（level setmethod，LS）[28,59]、流体体积 法 （volume of fluid method， VOF）[60-61]和格子-玻尔兹曼法（lattice Boltzmann method，LB）[62-64]。以上方法对比见表3。

图7 聚焦流微通道结构中两相流体相互作用过程数值模拟［28，50，52］

表3 模拟方法对比

2.1 LS方法

2.1.1 保守LS方法

保守LS 方法的关键在于附加维度的引入[65]，是基于空间曲面的隐函数表达，最初由Osher 和Sethian[66]于1988 年提出，随后在1994 年Sussman、Smereka 和Osher 等将其进一步发展。该方法可以简单地追踪移动的界面和形状，精确求解具有复杂拓扑变化的不可压缩两相流问题，并准确地表示界面法向和曲率等界面变量[25]。目前已成为模拟不可压缩多相流问题的一种有效方法，原理如图8 所示。此外，由于LS 方法采用的是光滑的距离函数来捕捉相界面，各个物理量可以在界面上光滑连续地过渡，相界面的捕捉效果好。Bashir 等[67]采用保守两相LS 方法模拟交错流微通道内的液滴成型过程，重点研究了微通道壁面润湿性对液滴成型的影响。研究结果表明两相流中液体之间的竞争润湿特性会导致不稳定的流型，从而破坏均匀液滴的正常成型。特别是对于疏水润湿区域而言，连续相流体毛细数Cacont小于0.02 时，接触角的影响作用越来越显著，后续试验也验证了该模拟结果。Naeimirad 等[24]同样采用保守LS 方法模拟受驱动力影响的不同界面形貌，探究液滴成型方式转变与系统量纲为1 参数（流速比、黏度比、毛细管数等）变化的构效关系。结果表明，增大流速和降低界面张力会减小液滴尺寸，从而导致液滴成型方式由滴流向喷射转变。

2.1.2 优化LS方法

采用水平集函数的经典控制方程，在经过多步计算后，流场的畸变将导致水平集函数的梯度变得过大或者过小，造成水平集函数本身的失真。该问题被称为“质量损失”，是保守LS 方法所固有的，会导致界面的法线、平均曲率和高斯曲率的不准确近似，并且不能通过提高计算精度而解决。为了克服“质量损失”问题，保持水平集函数为一个真正的符号距离函数，需要对它进行重新初始化，此过程被称作“重定距”[28]。然而，函数初始化过程总伴随着界面位置的移动，造成质量的损失，导致质量不守恒。此外，改善初始化步骤来矫正质量守恒又会增加计算时间，提升计算成本。近年来，通过优化控制方程、添加新控制方程或耦合其他模型来获取优化LS方法，已成为解决保守LS方法“质量损失”问题的新途径。Lan 等[27]通过在经典水平集控制方程中引入两个额外的附加项，来补偿水平集函数的失真，数值模拟所得液滴尺寸的预测值与试验测量值保持较高的一致性。在随后的研究中，Xiong 等[26]同样借鉴该方法来研究同轴流微通道内的液滴成型过程中流场变化和旋流强度。Lan 等[28]通过在水平集函数中设置一个新的控制方程来克服“质量损失”，并将同轴流微通道中液滴形成过程的模拟结果与已有报道的试验结果相比较，验证了优化LS 方法的可行性和准确性。而Marchandise等[68-69]则采用不连续的Galerkin（DG）方法来执行水平集输运方程，确保了质量守恒。Olsson 等[70-71]利用双曲正切函数代替传统的符号距离函数，进而减少质量守恒误差。Wong 等[29]将保守LS 方法与Carreau-Yasuda 应力模型结合，用以监测交错流通道内液滴成型过程中的尺寸变化。通过经验获得的流体流变和物理特性数据服从该应力模型，为模拟计算提供了理论依据。模拟结果表明，剪切成型液滴的尺寸变化由被测流体的流速、黏度等物理特性控制。随着流速比的增大，数值模拟计算对液滴平均直径的预测越准确。在测试的最高流速比（0.125）下，水滴平均直径预测值与测量值的差异约为11%，而在较低的流速比下，则显示出相当好的一致性。Li 等[30]通过在保守LS 方法中耦合静电模型，模拟研究了外加静电场对液滴成型过程影响。数值模拟表明，电场和电荷相互作用产生的电场力，对液滴形成动力学起着重要的控制作用。Gutiérrez等[17]在保守LS方法的基础上，通过引入拉格朗日-欧拉函数和浸入边界方法，提出了一种研究复杂几何形状中液滴和气泡演变问题的新方法。具体而言，首先利用保守LS 方法来处理多相域，同时控制质量守恒。随后引入拉格朗日-欧拉函数来优化仿真区域，在此过程中控制体积的大小被均匀化，网格在重要区域得到细化。因此，移动网格将跟随气泡的运动，从而减少了计算域大小、提高了网格质量，显著降低了计算资源的消耗。最后，利用浸入边界方法处理复杂的几何结构，并在拉格朗日-欧拉框架内重新生成内部边界。该研究方法能够处理完整的非结构化网格，增加了该模型的适用性。

图8 LS方法原理

2.2 VOF方法

2.2.1 VOF基础方法

使用最为广泛的界面捕捉模型是由Hirt 和Nichols[72]于1981 年首次提出的VOF 方法。该方法采用几何重建策略来构造流体界面，利用VOF 函数追踪不同相的体积分数，并采用分段线性界面计算（piecewise-linear interface calculation，PLIC）方法来逼近多相界面[73-74]，原理如图9 所示。VOF 方法已被用于探究微通道内的液滴成型、聚并、融合、溶解等各种传热传质并存的过程[48]。Padoin等[43]和Kashid等[39]利用VOF方法分别探究了交错流微通道壁面润湿特性对等温气-液流动模式（泰勒流和分割流）和柱塞状液滴成型机制的影响规律。Chen 等[37]采用VOF 方法研究了交错流微通道壁面润湿特性对液滴成型过程的影响，获取了毛细管数Ca和通道壁面润湿角与液滴成型方式和液滴尺寸控制的相互作用关系。Sontti 等[34]采用VOF 方法研究了交错流微通道中非牛顿流体液滴的成型过程和方式，并获得了两相流速、界面张力等因素对液滴成型机制的影响规律。Glatzel等[75]则针对打印机喷头墨滴生成的实际问题，基于三维模型和VOF 方法模拟研究了墨滴的生成过程及影响因素。

此外，还有众多学者采用VOF 方法模拟研究特殊结构微通道中的液体动力学行为。针对优化的同轴流微通道结构，Heuberger 等[31]采用VOF 方法研究了分散相通道尖端开设方形槽道的特殊同轴流微通道装置内部的流体动力学行为以及液滴的成型方式。Chekifi[32]基于VOF方法，研究了在低雷诺数Re和毛细管数Ca条件下，锥形同轴流结构微流控装置中液滴的生成过程，确定了液滴成型方式、液滴形状尺寸、产生频率与流动条件（速度比、黏性效应）和通道几何尺寸的关系。针对优化的交错流微通道结构，Chandra 等[35]借助VOF 方法研究了主通道齿形壁面结构（齿形间距、大小等）对通道内部柱栓状液滴成型、大小、形状以及齿间压降的影响，获取了通道壁面膜和压降等关键特征参数。Zhu等[36]则结合三维模型和VOF 方法模拟研究了主通道横截面几何形状变化对微液滴演化、移动等动力学行为的影响，验证了液态水对低温燃料电池中阴极气体微通道几何形状的敏感性。Ren 等[44]利用VOF 方法模拟研究了离心场作用下，并联交错流微通道中液滴/气泡的成型、聚并以及分裂行为。针对交错流微通道阵列结构，Kobayashi 等[40-41]和van Dijke等[42]均采用VOF方法模拟研究了交错流结构阵列中多液滴成型过程以及速度场分布和压力梯度的变化规律。

图9 VOF方法原理

而在微通道内两相传热传质方面，Dai 等[33]基于VOF 方法，开发出一个用来分析和预测微通道内部液-液泰勒流的流动特性和传热行为的通用模型。同时，借助该模型证实了传热系数与流动条件之间的强烈依赖关系，解释了试验过程中液-液流动传热测量的巨大不确定性和困难。此外，模拟的准确性也被试验所验证（大部分Nusselt 试验测量值都落在数值模拟预测值±10%的范围以内）。而Li等[38]结合micro-PIV 和VOF 方法，探究了子弹状液滴的生成过程，预测了液滴的尺寸和流速，获取了液滴内外流场特性。相关研究发现将有助于优化微通道中柱塞状/子弹状液滴结构特性，从而提高多相流传质传热和混合效果。

2.2.2 优化VOF方法

尽管VOF 方法能够模拟需要追踪界面的多相问题，但VOF 函数中的体积分数在空间上是一个不连续的阶梯函数，造成了两个相邻网格的界面也是不连续的，且相关物理量在通过界面时也是不连续的，这个现象称为寄生流动（parasitic current）。这些非物理速度是由于界面曲率离散化的数值误差而产生的，常被用来计算表面张力。然而表面张力计算中所带入的数值误差被传递到动量方程中，此时，则需要速度项来平衡该误差，进而导致了伪速度。伪速度和寄生流动在微观尺度上更加严重，因为毛细管力占主导地位，这可能会影响微观尺度上模拟的准确性[76]。目前优化VOF方法的主要途径就是缓解数值方法造成的伪速度和寄生流动现象[48]。最近，多种用以最小化VOF 模型中的伪速度，进而缓解伪速度问题的创新方法已经被开发，并成功用于研究涉及孔尺度和微流体的微尺度问题[76-79]。

Wang[46]采用VOF 方法和连续表面力（continuum surface force）模型分别捕捉气-液界面和模拟流体表面张力，对锥形同轴流结构微流控装置中流体状态参数和通道外尺寸参数与泰勒气泡成型方式、频率和尺寸的构效关系进行了数值研究。类似地，Deng 等[45]同样采用VOF 模型和连续表面力模型探究了可控尺寸的单分散微液滴的成型过程。模拟结果表明，连续相流体的黏滞阻力、分散相流体的惯性力和液滴界面张力之间的平衡在很大程度上决定了液滴的成型方式及尺寸大小。并在模拟研究的基础上提出了预测液滴直径的经验关系式，同时得到了相关试验的验证（液滴直径的预测值与试验数据高度吻合，最大误差为9.2%，平均误差仅为3.8%）。Guillaument 等[47]使用VOF 模型和单流体（one-fluid）模型组合来追踪二维微通道中两相流界面，探究了同轴流微通道结构中界面张力和微通道壁面润湿特性对超临界CO2流体成型方式及相位反转机制的影响规律。Soh 等[48]利用VOF 方法捕捉多相物理量，并在表面张力计算中应用平滑操作以最小化伪速度。同时，耦合多相/多组分模型对不同组分的迁移进行追踪，模拟研究了交错流微通道结构中CO2气泡成型和CO2气泡在连续相流体硅油中的溶解过程。该优化模型有助于推进微通道中液滴溶解、血管中涉及传质的药物运输以及CO2强化采油率等过程的研究。

2.3 LB方法

2.3.1 LB基础方法

格子-玻尔兹曼（LB）方法是统计物理学和计算科学的最新发展[49]。该方法是局部的，几乎不需要进行集群、多线程或GPU 并行处理，与顺序处理相比，加速比最高可达20。因此，该方法特别适用于不可压缩的流动和多物理场（热流体、电流体动力学、磁流体动力学、传质和反应流）等流动模拟问题[65]，原理如图10 所示。与LS 方法不同，LB 方法并不会直接追踪界面。相反，该模型使用网格来定义整个求解域的节点矩阵，并将流体视为有质量无体积的粒子。粒子在网格的节点间流动，并在节点处相互碰撞。通过模拟确定这些粒子在各个方向上的分布函数，可从这些分布函数中计算得到诸如密度、速度等宏观参数。该分布函数有4种常用模型：着色模型、Shan-Chen模型、自由能模型和He 模型。着色模型以不同的颜色区分相位，并根据相邻粒子之间的相互作用对界面进行追踪。Shan-Chen模型则将多相流体视为非理想流体，可在局部动量不守恒的情况下，很好地解决相变问题。自由能模型使用总密度和密度差来确定两种流体的密度。He 模型主要用于不可压缩流体，通过主函数来追踪界面，其中界面张力被分子间的相互作用所取代。Frisch 等[80]在其研究中证明，可以使用有限差分近似将LB 模型视为Navier-Stokes 方程的稳定形式，对于复杂几何结构的模拟具有较高的计算效率和精度。Wang等[81]在研究液-固传质特性的数值模拟过程中发现，使用LB 方法可在更小网格单元尺寸和更短计算时间的条件下，获得相同精度的计算结果。

图10 LB方法原理

2.3.2 优化LB方法

然而在众多模拟研究中发现，宏观尺度模型在简单几何结构中的应用要明显优于LB模型[82-84]。这是因为，对于LB 方法而言，当使用反弹边界条件时，模拟的准确性可能会降低。为了获得较高精度而采用半反弹策略，则会增加计算时间。因此，针对简单几何结构，提高效率并保持高精度是LB 方法优化的重要方向[65]。Riaud 等[49]开发了一种模拟两相微流体相互作用的优化LB 方法。该方法使用Latva-Kokko 再分配策略实现相分离，而界面张力则通过Shan-Chen模型的表达式求解。通过引入粒子碰撞后的再分配步骤，提高算法的稳定性，同时使界面厚度保持恒定，不受界面张力的影响。因此，该优化方法提供了独立的界面张力控制，受黏度影响很小。与连续表面力模型相比，界面曲率更容易获取，而伪流动可以保持在与着色模型相当的低水平。随后，他们采用该模型研究了交错流微通道结构中气泡成型过程，并结合以空气/乙醇为工作介质的T形结构微通道中气泡成型试验，通过对界面性质、界面张力以及剪切力与界面张力之间竞争关系的研究，验证了该模型的有效性。

此外，多相/多组分流体系统在微流体领域中广泛存在，涉及液滴/气泡在微通道内的成型、聚并、分裂等过程。由于体系具有较大的比表面积，微通道壁面润湿特性会对微流体流动行为产生较大的影响。在LB 方法中，相界面问题的处理可以简单地转化为对流体粒子间以及流体与微通道壁面间相互作用力的控制[85]。通过耦合多相/多组分模型，则可实现微通道中多组分追踪，以及流体性质、微通道特征参数等对液滴成型过程作用机制的模拟研究。Dupin 等[62]结合LB 方法（自由能模型）和多相/多组分模型，在低雷诺数Re和低毛细管数Ca条件下，模拟研究了二维交错流微通道结构中多相流动力学行为，并对流体与微通道壁面相互作用、相界面伪速度等进行调控。Yang等[86]以交错流微通道结构中液滴成型过程为基础，将浸入式边界条件引入LB 方法（Shan-Chen 模型），模拟研究了细胞在封装过程中的三维旋转和形变特性。

2.4 LS+VOF耦合方法

LS方法和VOF方法是用以模拟研究复杂界面两相流相互作用规律及动力学行为最为广泛的两种方法。在LS 方法中，通过水平集函数追踪和捕捉界面。由于水平集函数具有连续性和光滑性特征，可以精确地计算空间梯度。然而，LS方法在保持体积守恒方面存在固有缺陷。相比而言，VOF方法本质上体积守恒，因为该模型中计算和追踪的是每个单元中特定相的体积分数，而不是界面本身。VOF方法的不足在于空间导数的计算，因为VOF函数本身在界面上不连续。由于不当方法致使表面张力离散化和表面曲率近似化，都会产生伪速度。为了克服LS 方法和VOF 方法的各自缺点，ANSYS FLUENT（ANSYS Inc.，美国）在第14.0 版更新中提供了LS和VOF的耦合方法（CLSVOF）。该方法采用分段线性界面构造（PLIC）的几何重建来执行重新初始化过程。曲率和界面法向由水平集函数计算获得，而界面的准确位置则通过平衡每个单元中的体积来调节，从而实现由VOF模型计算获得体积分数。该方法在重定距水平集函数的同时加强了质量守恒，而流体的表面张力和物理特性则用LS模型相类似的方法计算求解[87]。与VOF 方法相比，CLSVOF 耦合方法可以产生更精确的气液界面，特别是在气泡与分散相流体分离阶段，可获得的与试验更加吻合的结果。Dang等[50]在异构化聚焦流微通道中泰勒气泡成型过程的数值模拟研究中发现，相较于VOF方法，采用CLSVOF耦合方法模拟所得的前后两个气泡的模拟长度差异更小，更加吻合试验结果）。此外，研究结果还表明，在泰勒气泡成型过程中，通道壁面接触角和流体黏度对气泡的形状产生了较大影响，表面张力对其的影响则相对较小。

3 研究对象

微通道中微流体间相互作用的结果多样，借助于CFD数值模拟的方法，通过研究流体间的传热传质过程，可获取微通道中液滴成型、液滴聚并/融合、液滴溶解以及颗粒聚焦等过程原理和调控机制。

3.1 液滴成型

3.1.1 量纲为1参数

液滴成型是微流控研究的基础，对于微流体应用至关重要。在液滴成型过程中，从注射泵或者压力控制器引入的能量部分转化为界面能，促使液-液界面失稳，进而引发离散的液滴从分散相流体中分离[88]。微流体中可以发生各种性质的流体运动，这通常由相互竞争的物理效应所决定，例如力的平衡。两相流体相互作用形成液滴过程中液滴本身所受到平衡力可以归纳为两大类：分离力和附着力。尽管微通道结构复杂、尺寸各异，但液滴成型过程中所涉及的物理学特性和潜在成型机制相似（图11）。多种力（重力、黏性力、惯性力、界面张力等）相互作用的结果是液滴以不同的流动方式成型，具体而言包括：挤压成型（squeezing）、滴流成型（dripping）和喷射成型（jetting）[89]。液滴生成过程与界面流动有关，通常涉及3 个主要步骤：首先，分散相和连续相流体在两相交界处相遇，并形成非混溶界面；随后，非混溶界面发生大变形，进入不稳定状态；最后，不稳定界面自发破碎，衰变为相互独立的液滴。因此，界面张力在液滴成型过程中起着关键作用，而与界面张力有关的两个量纲为1 参数（毛细管数Ca和韦伯数We）在液滴成型过程中尤为重要。

（1）毛细管数Ca毛细管数Ca是黏性力与毛细管压力的比值，见式(1)。由于惯性力本身不取决于通道尺寸，而通道尺寸的减小将导致黏性力和毛细管压力的增大以及重力效应的削减。因此，如果通道尺寸足够小，黏性力和毛细管压力就会占主导地位。这一原因使得毛细管数Ca成为描述微流体中液滴产生过程最常用的量纲为1参数。在微流体中，毛细管力的取值范围为10-3～10。

式中，uj为流体流速；μj为流体黏度；γ为两相间的界面张力。其中下角标j既可代表分散相流体（j=disp），也可代表连续相流体（j=cont）。

（2）韦伯数We尽管对于大多数微流体流动而言，流体的惯性可以忽略。但是对于高流速条件下的射流和非线性气泡形成以及液滴和气泡从分散相中的脱离，流体惯性很重要。惯性力和毛细管压力的相互竞争的结果是产生了韦伯数We。对于大多数微流体流动，We＜1，见式(2)。

式中，ρj为流体密度；d为通道尺寸。

3.1.2 挤压（squeezing）成型

在毛细管数Ca较小的条件下，此时黏性力被通道壁的限制所取代，液滴以两个阶段的挤压方式成型，即液滴的出现和生长、与分散相流体分离。在第1阶段，分散相流体在连续相流体施加的侧向剪切力作用下，在流体流动方向的主通道内不断积聚。随着分散相流体前端液滴的不断增大，液相界面和通道壁面的间隙不断减小，并在液滴上形成压力梯度。当此压力梯度足以克服分散相液滴内部压力时，分散相界面挤压变形，形成颈缩。在第2阶段，在动态界面张力所驱动的Rayleigh-Plateau 不稳定作用下，分散相流体前端液滴在颈缩处分离，形成长方形或者子弹型（plug-shaped）的柱栓状液滴/气泡。挤压方式产生液滴过程中，液滴在脱落前几乎完全阻塞两相交界区域，液滴的大小与毛细管数无关，而由通道的尺寸和两相流速比决定[34,37,89]。Wang[46]和Sontti 等[34]以挤压方式在锥形同轴流结构微流控装置中分别获得了子弹型的柱栓状气泡和液滴。而Chen 等[37]则以挤压方式在交错流结构微流控模型中获得了长方形柱栓状液滴。

图11 液滴成型过程模拟

3.1.3 滴流（dripping）成型

当毛细管数Ca继续增大（更高黏度和/或流速），液滴成型方式由挤压成型转变为滴流成型，其中拖曳界面破裂的黏性力占主导地位，而界面张力作用则确保初生的液滴稳定不破裂。在同轴流和聚焦流结构微通道中，液滴滴流成型过程同样分为两个阶段：液滴的出现和生长；与分散相流体分离。与挤压方式成型略有不同，由于较大的黏性剪切力，分散相流体前段不断增大的液滴在完全阻塞微通道之前，即被连续相流体在颈缩处剪断，因此，所生成的液滴尺寸小于微通道尺寸，保持拟球型形貌，具有高单分散性特征[37]。在这种情况下，液滴尺寸既取决于流速，也取决于两相黏度比。高黏度的连续相流体在液滴与通道壁面之间将会遇到更大的流体阻力，颈缩位置与分散相流体前沿液滴距离将进一步缩短，其结果是液滴分离过程速度更快，所获得的液滴尺寸更小。而对于交错流结构微通道而言，液滴/气泡在从分散相流体分离后会首先附着在微通道壁面上。因此，存在液滴与微通道壁面分离，并在自身界面张力作用下形成拟球型的自由液滴的第3阶段。在此阶段，液滴成型的持续时间与通道壁面润湿特性高度相关。

3.1.4 喷射（jetting）成型

通过增加分散相或者连续相流速，毛细管数Ca将进一步增大，液滴由滴流方式成型向喷射方式成型转变[90]。在同轴流和聚焦流结构微通道中，液滴喷射成型过程依然分为两个阶段，即液滴的出现和生长、与分散相流体分离。但与挤压和滴流方式成型均不相同，在液滴喷射成型过程中，由于界面（Rayleigh-Plateau）不稳定性，当连续相流体和分散相流体的惯性力所施加的黏性力大于界面张力时，两相流体在交汇处形成连续的射流，直观可见分散相流体具有明显的“拖尾”现象，此射流在下游末端断裂成拟球型液滴[37]。由于射流状态的不稳定性，且液滴从分散相流体脱离过程受两相黏度和流速等多方面原因的影响，致使喷射方式所获得液滴具有多分散性特征。此外，液滴形成位置向下游移动，液滴成型位置与液滴之间的距离随着连续相流体Ca值的增大而增大。当Ca值超过临界值时，通道内将不能产生液滴，两相流体流动变为稳定的平行流动。为了避免液滴喷射成型，在高流速下应采用低黏度的连续相。而在交错流结构微通道中，液滴成型过程仍然包含液滴与微通道壁面分离的第3阶段，且分离时间与通道壁面状态和润湿特性高度相关。已有众多学者[27,32,45-46]结合试验观测和CFD 模拟方法，系统地研究了液滴/气泡在不同微通道结构中的滴流和射流成型过程，并从中获取两种成型方式的转换机制和调控方法。而Chen 等[37]和Guillaument等[47]更是探究了微通道壁面润湿特性等因素对液滴/气泡成型的影响。

3.1.5 成型方式转换

随着毛细管力的增加，分散相流体前端的液滴部分地、暂时地阻塞了连续相流体的流动，从而形成了液滴从挤压成型到滴流成型的过渡状态。由于此状态下，分散相流体局部或间歇地阻塞连续相流体流动，剪切应力和挤压压力的共存，因此，过渡状态下的液滴尺寸既取决于毛细管数Ca，也取决于两相流速比[37,89]。此外，大量试验和数值模拟结果表明，3种液滴被动成型方式可以通过改变与界面张力有关的量纲为1 参数（毛细管数Ca和韦伯数We），从而实现相互转换。然而，相较于分散相毛细管数Cadisp和连续相韦伯数Wecont，连续相的毛细管数Cacont对液滴成型方式有较大的影响，但不同剪切形式所对应的Ca值也存在较大差异（见表4）。这是由于液滴成型方式是由两相组分、两相黏度、两相流速、微通道尺寸结构、边壁润湿特性等多方面因素共同作用的结果[91-92]。

3.2 液滴的聚并/融合/溶解

3.2.1 液滴聚并

CFD数值模拟为微通道内液滴之间的碰撞过程提供了更为直观的研究方法，有助于研究在复杂微通道结构或者多场作用下相同或不同属性液滴间的相关作用关系。同属性液滴在微通道中的相互作用过程会发生破碎、聚并和再分裂（图12）。Anandan 等[93]通过在气泡逻辑门模型中引入交错流微通道结构作为气泡生成通道，模拟研究了气泡的生成过程，以及气泡以不同操作条件在两相邻微通道切线交汇处的聚并和分裂过程。而Ren等[44]模拟了离心力作用下并联交错流微通道中液滴的成型过程，以及液滴在通道交汇处的碰撞、聚并和再分裂过程。结果表明，随着离心力的增大，液滴成型频率增大，液滴尺寸减小。对于并联交错流微通道而言，微通道与转动中心间距的增大，致使液滴在成型过程中所受到的离心力作用增大，其结果是距离转动中心最远的微通道中所获得液滴尺寸最小，但液滴成型频率最快。此外，离心力的作用同样会加剧并联微通道所生成的液滴在汇集通道中的碰撞和聚并过程。

3.2.2 液滴融合

不同属性液滴在微通道中的融合属于两相传质问题，是一个相对复杂的过程（图13）。Trivedi等[25]借助于试验观测和CFD数值模拟研究了交错流微通道结构中不同属性液滴间的相互融合过程。研究发现，不同属性液滴间的融合可以分为接触、混合、颈缩和脱离4 个连续过程。具体而言，在第1阶段，当液滴A被连续相流体输送到两相相互作用处时，表面力驱使液滴A界面与分散相流体界面接触连接，在该过程中毛细作用过程自然发生。在第2阶段，液滴A注入分散相中，与分散相流体相互混合，促进了混合液滴B在分散相流体前沿的出现和长大。在第3阶段，受到连续相流体施加的横向剪切力，分散相流体不断向前端积聚，致使分散相前端液滴不断长大，分散相后端出现“颈缩”。在第4阶段，当液滴增长到临界尺寸，混合液滴B从分散相颈缩处断裂，并逐渐脱离微通道顶端壁面，形成自由液滴B。CFD模拟研究表明，在适当的操作条件下，该过程可以自动同步。

表4 液滴/气泡成型方式及转换条件

3.2.3 气泡溶解

微通道中气泡在连续相流体中的溶解过程同样属于两相传质问题，试验方法对于该过程的研究往往存在极大的难度，仅仅停留在针对气泡尺寸变化的直观测量。而借助CFD 模拟方法，通过研究气泡在连续相流体中溶解阶段的内外流场变化特性，从而获取气泡微通道溶解过程机理和调控方法（图14）。Soh 等[48]通过耦合传质模型和VOF 模型，在二维微通道中模拟研究了微通道内CO2气体剪切成型和硅油中的溶解过程，获取了CO2气泡溶解过程中长度和体积的变化规律，导出了溶解气泡长度和体积随时间变化的解析解，并与恒定传质系数模型等3 种模型的模拟结果以及试验结果进行比较，进而证明了模型的准确性和可靠性。在此过程中，流体体积模型用于捕获多相物理量，并在表面张力计算中应用平滑操作以最小化伪速度。对于相迁移的多组分追踪，采用了α因子和驱逐操作方法，以确保在相应的阶段追踪正确的对象。这项研究所开发的耦合模型可用于研究微通道中液滴的溶解、血管内药物的运输等传质过程。

3.3 颗粒聚焦

水动力聚焦是微流体中应用最广泛的技术之一，已被广泛应用于化学/生物分析，包括流式细胞术、单分子检测以及用于快速化学和酶动力学研究的层流混合器[94]。通过在弯曲微流体通道中引入Dean流，Mao等[95]基于三维流体动力聚焦技术，开发了“微流控漂移”技术，实现了单层平板微流体装置中细胞群或颗粒群的三维流体动力聚焦（图15）。借助于CFD模拟，颗粒群在微通道内的三维流动聚焦过程得以研究，具体分为3 个阶段：第1阶段，包含颗粒的流体（以下简称为“颗粒流”）和鞘层流分别从平行进口引入微流控芯片中，并在进入90°曲线通道前进行混合；在第2 阶段，在曲线通道处，所引入的Dean 流在垂直方向上以双涡的形式对颗粒进行横向加速，将颗粒从通道顶部和底部扫向通道中心平面，并进一步将颗粒整体拖向通道中心区域，这一步被称为“微流控漂移”，有效地实现了颗粒在垂直方向上的流动聚焦；在第3阶段，颗粒流受到垂直于流动方向的两侧对称旁路鞘层流的水平剪切作用，促使原本聚焦在微通道中心区域的颗粒进一步在中心轴线聚焦，实现单颗粒在中心轴线方向的依次排列。该研究成果有助于实现细胞在微通道中的流动聚焦检测。

图12 液滴聚并过程模拟［44，93］

图13 液滴融合过程模拟［25］

图14 气泡溶解过程模拟［48］

4 结语与展望

在受限微流控通道中操作两相流体获得液滴，进而实现微通道液滴聚并/融合、液滴溶解、颗粒聚焦等复杂传质传热过程，已被广泛应用于许多科学领域。当下数值模拟方法的快速发展为预测和分析微流控通道内的多相流问题提供了更为直观有效的帮助，但是随着微通道结构的复杂化、多相流体系的复杂化以及多场耦合微通道内多相流相互作用的复杂化，也给数值模拟方法带来了前所未有的挑战，同时促使着更为完备的数值模拟技术的革新。在未来，数值模拟技术在液滴微流控探索应用的发展方向将会集中于以下几个方面。

图15 颗粒聚焦过程模拟［94-95］

（1）数值模拟方法的准确性和适用性的提升。随着微通道结构复杂化和小型化集成，现有基础数值模拟方法，包括LS 方法、VOF 方法和LB 方法等，由于本身函数/模型的局限性或缺陷，难以避免试验结果与数值模拟结果的吻合性较差。此前，一些研究人员试图通过耦合其他模型或优化本构方程的方式，克服LS方法的“质量损失”、VOF方法的伪速度和寄生流动以及LB 方法的（半）反弹边界问题，进而提升数值模拟结果的准确性，但涉及组合方法或者优化方法的适用范围仍需进一步扩展。

（2）模拟对象从牛顿流体体系向非牛顿流体体系的拓展。微通道中液滴成型的研究大多采用传统的牛顿流体体系，而对非牛顿流体的研究很少。非牛顿流体具有剪切力和变形速率非线性相关的特性，不同的力，包括表面张力和黏性力，控制着液滴动力学。当涉及非牛顿流体时，流体的流变学可能很重要。然而，当前没有建立单独的本构方程来描述所有流体的流变图。了解微流控通道中的非牛顿流体液滴成型的动力学机制对于确保液滴尺寸、形态和产量等特性要求至关重要。此前，一些研究人员专注于分析具有黏弹性特性的液滴的动力学和形变，但对于其他非牛顿流体体系产生可控液滴背后的相关物理学性质的深入了解仍然有限。

（3）液滴输运过程动力学行为及内部流场特性研究的深入。微通道受限空间内的多相流相互作用过程以及液滴/气泡在复杂几何结构的动力学行为对于许多科学和工程应用而言非常重要。此前，这两方面的研究常常通过简化控制方程的方式来寻求一些近似的解析解。随着光学观测技术日新月异的进步，为揭示复杂化/小型化微通道内多相流相互作用过程中，液滴内部流场特性及变化规律提供了有利的支撑。而数值模拟软件的发展为定量化预测液滴在受限空间内的动力学行为（自转、公转、跳跃）提供了可能，加深了对微通道内传质传热过程强化的认知。