同轴微通道中双乳液滴生成的数值模拟研究

Numerical Simulation Study on the Generation of Dual Emulsion Droplets in Coaxial Microchannels

摘要: 本文基于COMSOL Multiphysics 6.2与水平集(Level Set)法追踪油‑水‑油界面演化过程，构建了二维轴对称同轴微通道模型，探究聚焦孔直径与出口流道宽度之比(w_n)、颈长(L₃)和开角(θ₃)对双乳液滴生成的影响。结果表明：该模型中双乳液滴存在滴流与射流两种生成模式，滴流模式下界面张力主导、单分散性好；射流模式下粘性剪切强、生成频率高。几何参数分析表明：增大w_n可减小外液滴体积、增加生成频率，过大易导致液滴合并；增加L₃主要延长射流长度，对液滴体积影响较小；减小θ₃可提升剪切力，减小外液滴体积、增加生成频率，内液滴体积变化不明显。综上，几何参数可有效调控双乳液滴尺寸、频率等，为微流控双乳液滴在药物载体、功能胶囊等领域的应用提供了数值依据和设计参考。

Abstract: This study used COMSOL Multiphysics 6.2 and the Level Set method to track the evolution process of oil-water-oil interfaces. A two-dimensional axisymmetric coaxial microchannel model was constructed to investigate the effects of the ratio of the diameter of the focusing hole to the width of the outlet channel (w_n), the neck length (L₃), and the opening angle (θ₃) on the generation of double emulsion droplets. The results show that in this model, there are two generation modes for double emulsion droplets: droplet flow and jet flow. In the droplet flow mode, the interfacial tension dominates and the monodispersity is good; in the jet flow mode, the viscous shear is strong and the generation frequency is high. Geometric analysis shows that raising w_n reduces outer droplet volume and increases generation frequency; an excessively large w_n is prone to cause droplet merging; increasing L₃ mainly extends the jet length, with a relatively small impact on the droplet volume; reducing θ₃ can enhance the shear force, reduce the volume of the outer droplet, increase the generation frequency, and the volume change of the inner droplet is not significant. In summary, the geometric parameters can effectively control the size and frequency of double emulsion droplets, providing numerical basis and design references for the application of microfluidic double emulsion droplets in drug carriers, functional capsules, etc.

1. 引言

微流控技术因其具有在微尺度下精确操控流体的优势，被广泛应用于生物医学、药物递送和材料合成等领域[1]-[3]。微液滴技术作为微流控的重要分支，可高效制备单分散性高的液滴，并进一步用于生成双乳液滴(如W/O/W或O/W/O)，克服了传统方法中液滴尺寸不均一的难题[4] [5]。双乳液滴微球具有独特的核壳结构[6]，在封装活性物质、可控释放和多相反应方面具有巨大潜力[7]。然而，双乳液滴的制备过程涉及多相流体复杂作用，其生成机制和稳定性等仍面临挑战。目前，生成双乳液滴的微流控装置结构主要包括双T型通道[8] [9]、双十字聚焦通道[10] [11]，同轴毛细管结构[12] [13]等。在双乳液滴生成中，双T型结构因非对称剪切作用易产生多个单分散且不稳定的液滴；双十字聚焦通道则在交汇处易受扰动，导致液滴尺寸不均一且操作复杂。而同轴毛细管微通道因嵌套式结构和同心圆设计，液流稳定、可控性高，成为当前研究的热点。Umbanhowar等[14]最早利用该装置研究了高度单分散的单乳液滴生成规律，发现液滴尺寸与毛细管直径、两相流速和界面张力密切相关。Abate等[15]通过实验证明同轴流动聚焦装置可一步生成多重乳液滴，但其尺寸均一性和稳定性仍受流体参数(如毛细管数Ca、韦伯数We)的显著影响。Zhang等[16]利用微尺度3D打印了一种双同轴结构装置，能够稳定生成单、双乳液滴，双乳液滴的制备提供了一种新的途径。

液滴微流控系统在搭建及实验环节工艺要求高、成本高，多相界面动态行为、局部剪切力分布等信息难以通过高速摄像完整捕捉，仅通过实验难以深入探究微观液滴形成机理[5] [17]。基于有限元的流体力学模型可精准捕捉流体界面瞬时运动，适于多参数优化，并能耦合温度场、电场等揭示其复杂作用机制[18] [19]。Chen等[20]运用VOF法模拟了流动聚焦装置中双乳液滴形成过程，揭示了粘度比、流速比和界面张力比对双乳液滴形态的影响。Nabavi等[21]采用VOF-CSF模型建立了结合流动聚焦与共流的二维轴对称装置，研究了流速、流体性质及几何形状对双乳液滴特性的影响。结果表明：增大毛细管孔口尺寸可有效抑制喷射流；而低流量下增加注入口与收集口间距，则会出现复合液滴与单液滴交替生成的现象。

尽管已有研究对流动聚焦微通道中双乳液滴的形成进行了较深入的数值模拟分析，针对同轴微通道中多相流动特性及其对双乳液滴生成行为影响的系统研究仍相对有限。现有的模拟多采用VOF法结合CSF处理界面，其在处理界面拓扑变化及多相相互作用方面虽有一定优势，但在界面清晰度与质量守恒方面仍面临一定挑战。为此，本文基于COMSOL Multiphysics 6.2仿真平台，采用水平集法建立了二维轴对称数值模型，对同轴微通道中油包水包油(O/W/O)型双乳液滴的生成过程进行系统模拟。本研究重点通过系统调节聚焦孔直径、聚焦孔颈长度与聚焦孔开角三类几何结构参数，并基于三维体积分数图、压力场以及速度场分析液滴的生成特性和机理，深入探究几何结构对双乳液滴体积、生成频率及流动形态的影响规律，为未来实现双乳液滴的精确调控提供参考依据。

2. 数值分析方法及可行性验证

2.1. 控制方程

微观尺度低流速通道中，假设流体为不可压缩的层流，忽略重力影响，采用二维轴对称模型进行描述，液滴的自由界面由“层流–水平集”法定义与捕捉。根据质量和动量守恒原理，流动过程由以下连续性方程(1)与不可压缩的Navier-Stokes方程(2) [22]描述：

$\rho \nabla \cdot u=0$ (1)

$\rho \left(\frac{\partial u}{\partial t}+u\cdot \nabla u\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla u+\nabla u^T\right)\right]+F_{st}$ (2)

式中：u为流体速度矢量(m·s⁻¹)；ρ为密度(kg·m⁻³)；μ为动力黏度(Pa·s)；$\rho \nabla$ 为单位体积流体所受的压强梯度；$\mu \left(\nabla u+\nabla u^T\right)$ 为流体因剪切变形引起的黏性应力；$F_{st}$ 为表面张力引起的体积力(N·m⁻³)，通过水平集函数耦合至动量方程。

水平集函数$\varnothing$ 用于追踪油–水–油三相界面，其在计算域内满足输运方程[23]：

$\frac{\partial \varnothing }{\partial t}+u\cdot \nabla \varnothing =\gamma \nabla \cdot \left({ϵ}_{ls}\nabla \varnothing -\varnothing \left(1-\varnothing \right)\frac{\nabla \varnothing }{\left|\nabla \varnothing \right|}\right)$ (3)

式中：$\varnothing$ 为水平集变量，取值在0 (壳层油相)到1 (内水相/外连续油相)之间平滑过渡；${ϵ}_{ls}$ 为界面厚度参数，控制相界面扩散尺度；$\gamma$ 为界面重新初始化参数，用于保持水平集函数的符号距离特性。密度与黏度在界面附近由水平集函数平滑过渡，由以下公式计算：

$\rho ={\rho }_1+\left({\rho }_2-{\rho }_1\right)\varnothing$ (4)

$\mu ={\mu }_1+\left({\mu }_2-{\mu }_1\right)\varnothing$ (5)

式中，${\rho }_1$ 表示分散相密度，${\rho }_2$ 表示连续相密度，${\mu }_1$ 表示分散相粘度，${\mu }_2$ 表示连续相粘度。

表面张力$F_{st}$ 通过以下连续表面力(CSF)模型方程引入：

$F_{st}=\sigma \kappa \delta (\varnothing )n$ (6)

$\kappa =\nabla \cdot \left(\frac{\nabla \varnothing }{\left|\nabla \varnothing \right|}\right)$ ，$n=\frac{\nabla \varnothing }{\left|\nabla \varnothing \right|}$ (7)

其中：σ为油–水界面张力系数(N·m⁻¹)；$\kappa$ 为界面曲率；n为界面单位法向量；$\delta (\varnothing )$ 为Dirac分布函数，用于将表面张力集中于界面附近。

2.2. 几何模型及边界条件

本文构建的用于模拟双乳液滴生成的二维轴对称同轴聚焦微通道几何模型如图1所示，利用轴对称特性可显著降低计算资源，同时维持界面演变高保真度。几何结构主要由三个同轴区域组成：内相入口1通道、中间相入口2通道及最外层的连续相入口3通道，均汇集于聚焦孔处，该结构由聚焦孔直径(D₃)、颈长(L₃)及开角(θ₃)三个关键几何参数定义。模型的尺寸如下：D₁ = 30 μm、D₂ = 90 μm、D₃ = 300 μm，L₁ = 900 μm，L₂ = 960 μm，总长L₃ = 1000 μm。本文采用三角形网格离散计算域，以确保水平集方法精确捕捉界面变化，并完成了网格无关性验证。模拟流体体系如下：1,6-己二醇二丙烯酸酯作为内相(油相)与外相(连续油相)，以2 wt%的聚乙烯醇(PVA)水溶液作为中间相(水相)，初始控制内相流速为0.05 m/s，中间相流速为0.05 m/s，外相流速为0.04 m/s，在所有模拟中，时间步长固定为0.001 s。内相、中间相与外相依次对应为流体1、2、3，其在25℃下的物性参数详见表1所示[24]。

Figure 1. Geometric dimensions of the coaxial microchannel model

图1. 同轴微通道模型几何尺寸

Table 1. Physical properties of three-phase fluids

表1. 三相流体的物理特性

流体	密度(kg/m^3)	粘度(Pa·s)	界面张力(N/m)
内相流体	1000	0.00671	0.005
中间相流体	1000	0.00195	0.005
外相流体	1000	0.00671	0.005

2.3. 网格可行性分析

为了研究网格细化程度对于生成双乳液滴的准确性影响，同时节省计算资源，选择了如图2所示四种不同细化程度的网格(截取自聚焦孔区域)进行可行性分析。结果显示：细化程度越高，三角形网格数量越多，差异主要体现于流体聚焦前沿区域，其余流型稳定部分则无明显区别。体积分数图可看出网格过粗会导致计算精度不足、界面失真，甚至引发液滴破裂现象；随着网格逐渐精细化，计算可捕捉到卫星液滴[25]，且相界面变得更精细。由表2可看出，内外液滴的体积比随着网格数量的增加而变化，当网格数从12,823增加到90,487时，液滴的体积比从0.04556增加到0.11872，增加了161%，当网格数从90,487增加到340,828时，液滴的体积比增加了7.2%，而计算时长从533 min增加到2725 min。因此考虑到计算资源成本，选取平均单元网格质量为0.9663、单元数量为90,487的网格3用于后续研究。

Figure 2. Local diagrams and phase diagrams of different unit size networks

图2. 不同单元尺寸网格局部图及其相图

Table 2. Data of different grids

表2. 不同网格的数据

网格类型	三角形网格数量	平均单元质量	内外液滴体积比	计算时长(min)
网格1	12,823	0.9535	0.04773	79
网格2	19,877	0.9552	0.06257	91
网格3	90,487	0.9663	0.13471	533
网格4	340,828	0.9639	0.14585	2725

3. 结果与讨论

3.1. O/W/O型双乳液滴典型流态模式

为验证本模型在几何结构与物理参数设置上的合理性，本文参考了Sattari等[26]在双共轴微通道中进行双乳液滴生成的实验研究，该实验采用一步法生产双乳液滴，本文模拟结果与Sattari等实验中所观测到的滴流与射流模式及其转变条件高度吻合。因此，本文成功预测了同轴微通道中生成双乳液滴的2种典型流态模式——滴流式和射流式，并揭示了其生成机制：内相流体在聚焦口处先被中间相包裹形成初级液滴，在外相剪切作用下发生颈缩并断裂，最终形成双乳液滴。整个过程可分为初级液滴的生长、颈缩和断裂三个阶段。

3.1.1. 滴流模式

Figure 3. Drop-flow mode of the double emulsion droplet experiment [26]

图3. 双乳液滴实验滴流模式[26]

Figure 4. Phase diagram under droplet flow mode (left) and pressure and velocity contour maps (right)

图4. 滴流模式下相图(左)及压力和速度云图(右)

如图4所示，左部分的体积分数图呈现出双乳液滴的滴流式生成过程，与图3实验滴流模式结果一致，该模式通常发生在低流速情况下，即内相毛细管数和外相毛细管数较低时，此过程周期性强，液滴单分散性强。在此模式下，界面张力主导双乳液滴的断裂过程，中间相在聚焦孔处呈现出周期性的“颈缩–夹断”动力学行为，内界面流体形状从椭圆形逐渐转变为圆形，液滴在颈部处收缩后自然脱落，形成尺寸均匀的双液滴。从图4的压力和速度图可看出，在液滴形成初期，内相区域的压力和速度最高，且内相和中间相的压力均显著高于外相区域，该压力驱动其突破界面张力与粘性剪切力，驱使中间相流体均匀地包裹内相。液滴生长时，颈部形成逆向压力梯度，协同Marangoni效应使液膜变薄拉长，引发界面失稳与收缩。随液滴增大，当颈部达临界尺寸，压力与速度达到峰值，内外界面依次断裂，完成液滴脱离。整个过程表明，滴流模式下的液滴生成是压力、界面张力与粘性力相互竞争耦合的结果。该同轴微通道可实现三相流体径向对称流动，有效避免液滴偏心与形态失稳，从而易于生成高单分散性、球形度佳的双乳液滴。相较于双T型或双十字型结构，该方法适用性更广且一步即可完成制备。

3.1.2. 射流模式

Figure 5. Jet mode of the double emulsion droplet experiment [26]

图5. 双乳液滴实验射流模式[26]

Figure 6. Phase diagram under jet flow mode (left) and pressure and velocity contour maps (right)

图6. 射流模式下相图(左)及压力和速度云图(右)

随着中间相或外相毛细数的增加，粘性剪切力开始发挥更显著的作用，双乳液滴的生成模式由滴流向射流模式转变，即如图5所示的实验结果，本文模拟结果图6所示与图5实验吻合良好。降低中间相表面张力会导致粘性力作用增强、毛细管数增大，使中间相包裹内相的界面更易被连续相拉伸，断裂位置会向下游移动，形成稳定的射流。根据图6可知，射流区压力分布均匀，无滴流模式下的局部高压集中，仅在射流末端出现周期性高低压变化。内外相流速均高于中间相，中间相包裹内相进入通道后被高速外相流拉伸，于下游形成长射流。下游Rayleigh-Plateau不稳定性主导的周期性压力扰动，使得射流出现颈缩现象，主液滴断裂后其连接颈部分可能无法被主液滴完全吸收，从而导致卫星液滴的生成。此模式下，粘性剪切力代替界面张力成为主导，压力场转化为界面不稳定性的媒介。这种机制导致射流模式下双乳液滴生成频率显著高于滴流模式，且液滴尺寸小。且外相流速越高，射流稳定性越强，断裂位置处越远。此机制对于需要高通量但不要求极高单分散性的乳液制备应用具有重要意义。

3.2. 聚焦孔直径对双乳液滴生成的影响

Figure 7. Volume fraction diagrams of the double emulsion droplets for different focusing hole diameters (left) and pressure and velocity contour maps (right)

图7. 不同聚焦孔直径的双乳液滴体积分数图(左)及压力和速度云图(右)

微流道装置的聚焦结构对双乳液滴的生成至关重要，本研究在三维模型中考察了聚焦孔直径(D₃)对双乳液滴生成的影响。设置D₃为120、150、180、210、240 μm，对应其与出口流道宽度之比w_n = 0.4、0.5、0.6、0.7、0.8，其余参数保持不变。相关结果见图7、8。如图7所示，随着聚焦孔直径增大，射流长度缩减，断裂长度减小，逐渐转变为滴流，当w_n ≤ 0.6时，聚焦孔较窄，流场的压力、速度梯度较弱，仅能形成单一大体积双乳液滴，当w_n ≥ 0.7时，聚焦孔变宽，流场压力及速度峰值区显著增多，易形成多个小体积双乳液滴。

由图8可知，随着w_n由0.4扩大至0.7，外液滴体积持续减小，从2.08 μL降至1.10 μL；双乳液滴的生成频率持续升高，从114.03 Hz升至212.31 Hz，因为w_n增大至0.7使流场扰动即压力和速度梯度增强，液滴生成频率加快，导致毛细管数增大。当w_n = 0.8时，外液滴体积明显回升，到1.75 μL，这是因为聚焦孔过宽，流场对液滴的约束减弱，小液滴易发生合并，导致外液滴体积变大；生成频率显著下降至151.06 Hz，体积回升与频率下降均源于液滴合并对生成节奏的打破。而在整个w_n变化范围内，内液滴体积略有减小但变化不大，说明w_n对双乳液滴“内层核心液滴”的体积调控作用很弱。因此，改变聚焦孔的直径w_n ≤ 0.6时，易获得稳定性强的单一大体积双乳液滴，其外层包裹液体积大、载药量高，实现药物的缓慢释放，选择适中的w_n，如0.7时，易获得高频率小体积双乳液滴，适于批量生产。

Figure 8. (a) Volume of inner and outer droplets; (b) Frequency of formation of double emulsion droplets

图8. (a) 内外液滴体积；(b) 双乳液滴生成频率

3.3. 聚焦孔颈长对双乳液滴生成的影响

本节研究了聚焦孔颈长L₃ (设为30、60、90、120、150 μm，聚焦孔直径固定为150 μm，其余参数同3.2节)对双乳液滴生成的影响。如图9体积分数图所示，所有工况均为射流模式，且射流长度随L₃增加而增长。当L₃在30 μm∼120 μm范围内时，射流在出口处保持稳定，形成连续均匀的射流束并断裂为完整外液滴。相关模拟结果见图9、图10。当L₃ = 150 μm时，此时聚焦孔颈过长，射流在流动过程中失去有效约束，伴随出现卫星液滴，由于卫星液滴的尺寸较小，运动速度高于双乳液滴，一段时间后追赶上双乳液滴并与之融合。从右侧压力和速度云图可知，随L₃增大，聚焦孔区域的流速梯度变平缓，同时压力分布的集中程度降低。

从图10分析可知，随着聚焦孔颈长L₃在30 μm~150 μm范围内的增大，内液滴体积基本保持稳定，说明颈长变化对包裹的内相流体体积影响有限，外液滴体积则呈现出缓慢上升的趋势，主要是因为L₃增大，削弱了聚焦孔结构对射流的约束使得射流直径略有增加，同时液流在孔颈内的停留时间延长，导致射流在断裂前能够积累更多连续相流体，从而缓慢增大了最终液滴体积；双乳液滴的生成频率先略升高后略降低，在L₃ = 90 μm时，射流稳定、液滴体积可控、生成频率最高。总体而言，聚焦孔长度的改变对双乳液滴的生成过程影响不大，效应主要体现在射流段的拉长，控制颈长小于150 μm，可避免卫星液滴生成。

Figure 9. Volume fraction diagrams of the two emulsion droplets with different neck lengths of the focusing holes (left) and pressure and velocity cloud diagrams (right)

图9. 不同聚焦孔颈长的双乳液滴体积分数图(左)及压力和速度云图(右)

Figure 10. (a) Volume of inner and outer droplets; (b) Frequency of formation of double emulsion droplets

图10. (a) 内外液滴体积；(b) 双乳液滴生成频率

3.4. 聚焦孔开角对双乳液滴生成的影响

聚焦孔开角(θ₃)是影响孔口附近连续相流动收敛与压力分布的关键几何参数。本节固定其他工艺参数，仅改变θ₃ (取值90˚、75˚、60˚、40˚、30˚)，研究其对双乳液滴生成的影响，相关体积分数分布、压力速度云图及统计结果见图11与图12。由图11可知，当θ₃减小至30˚时，流道收敛加剧，显著提高了局部剪切率与毛细管数，从而使Marangoni效应增强，压力梯度更显著，Marangoni效应与Rayleigh-Plateau不稳定性耦合，聚焦孔流道收缩更剧烈，流体通过聚焦孔时的局部流速更高；当θ₃增大至90˚时，Marangoni效应减弱，界面张力梯度弱，压力与流速分布更均匀，流道收缩更平缓，界面与流场的耦合作用被削弱。

Figure 11. Volume fraction diagrams of the double emulsion droplets for different focusing hole diameters (left) and pressure and velocity cloud diagrams (right)

图11. 不同聚焦孔直径的双乳液滴体积分数图(左)及压力和速度云图(右)

由图12统计结果可知，随聚焦孔开角θ₃从90˚减小至30˚，外液滴体积显著减小，液滴生成频率逐渐增加，原因：聚焦孔开角越小，毛细管数越大，Marangoni效应越明显，加速外相流体的剪切力，加快液滴分离，因此外液滴体积减小，液滴生成频率增加，内液滴体积基本稳定。若实际应用中需小尺寸、高产量、薄壳的双乳液滴，可选择较小的聚焦孔开角(如30˚)的模型以强化流场约束与剪切力，需大尺寸、厚壳度的液滴，则可增大开角。

Figure 12. (a) Volume of inner and outer droplets; (b) Frequency of formation of double emulsion droplets

图12. (a) 内外液滴体积；(b) 双乳液滴生成频率

4. 结论

本文基于水平集法对同轴微通道中双乳液滴的生成与演化进行了数值模拟，成功复现了滴流和射流两种典型模式，并从压力场与速度场揭示了其生成机制。重点考察了聚焦孔直径、颈长及开角等几何参数对液滴尺寸与生成频率的影响：较小直径、较短颈长和较小开角有利于形成外液滴体积小、壳层薄、生成频率高的模式，适用于快速释放型药物载体；而较大直径、较长颈长和较大开角则对应外液滴体积大、壳层厚、生成频率低的模式，更适合高负载功能微胶囊等应用。二维轴对称模型基于完美轴对称假设，无法考察实际加工中通道偏心、截面非圆度等几何偏差对流动及液滴一致性的影响。本文模拟处于低雷诺数层流区，惯性力可忽略；若流速提高进入惯性主导区，则可能出现三维失稳、涡脱落等现象，超出该模型的预测能力。因此，为更精确预测实际复杂工况，未来需开展全三维数值模拟。此研究为微流控双乳液滴的定制化制备提供了结构设计依据，后续可通过多参数协同调控，进一步提升液滴性能的精准可控性，拓展其在生物医学、功能材料领域等方面的应用。

NOTES

^*通讯作者。