互不相溶的液−液两相体系广泛存在于自然界以及人类的生产、生活过程中，在能源热转换、石油化工、材料制备和生物医学等领域均扮演着重要角色。对由多个不混溶相组成的高黏度流体组分进行加工，将一种液体分散到另一种不混溶液体中，是药物、化妆品、药品和食品乳液制造过程中的关键步骤^［1-2］。因此，有必要研究剪切流中液滴的流体力学，以细致了解多相分散体的流变特性，准确掌握剪切流动中液滴的流变特性及其内在机理，这对研究液滴动力学行为的理论研究具有重要的学术意义。

从20世纪30年代Taylor^［3-4］的开创性研究开始，许多学者对液滴变形和破裂的研究做出了贡献。在实验研究方面，部分学者研究了一些无量纲数对液滴动力学行为的影响。张红光等^［5］通过实验验证了剪切流场中液滴的变形收缩系数和拉伸率之间的非线性关系，通过收缩系数的变化判定剪切对液滴形变的影响；林长志等^［6］在实验中研究了牛顿液滴在不相溶连续体相中的变形断裂，结果表明，液滴的断裂取决于毛细数和黏度比，并通过三种机理得以实现，分别为颈缩、端部夹断和毛细不稳定性；Salkin等^［7］结合实验和理论，研究了可变形物质（如液滴和气泡）对线性微障碍物的破碎动力学，通过黏度对比讨论了其对临界毛细数的影响。相比于实验，利用数值模拟方法能对液滴运动的一些内在机理与影响规律进行细致研究。Rallison等^［8］在二维剪切流的研究中，研究了黏度比和毛细数对液滴变形和破碎的影响；Stone等^［9］研究了低雷诺数黏性流场中液滴变形和破碎的动力学过程，分析了流动条件和物性参数对液滴尺寸分布的影响；Li等^［10］对浸没在低雷诺数和中等雷诺数的简单剪切流中的液体球体的内部、外部和表面流动进行了数值研究，结果表明，雷诺数Re=0的流线是液球表面上的闭合Jeffery轨道，也是液球内外的闭合曲线；王程遥等^［11］基于流体体积（VOF）液/液相界面追踪方法，建立了不可压缩水/油单乳液液滴动力学模型并进行数值求解，模拟了剪切流场条件下两个相同体积的液滴在碰撞过程中的相互作用及变形行为；骆政园等^［12］研究了流体惯性影响下线性剪切流场中弹性界面液滴的变形动力学行为，从而揭示了雷诺数对弹性界面液滴变形特性的影响规律及其内在机理。

本文基于格子Boltzmann（lattice Boltzmann，LB）颜色梯度模型，对剪切流场中液滴的动力学行为进行了数值模拟与详细分析。为了深入研究液滴界面对其内部不同位置流体的作用，研究中引入示踪粒子并对其运动位置和运动形态进行追踪，探索液滴在不同毛细数Ca、雷诺数和流体黏度比下的变形现象，在Re−Ca相图中区分了不同的变形及破裂模式，并对液滴内部不同初始位置的示踪粒子运动进行研究。

1 数学方法

在Liu等^［13］的两相颜色梯度LB模型中，两种不混溶流体分别表示为红色流体和蓝色流体。颜色梯度模型由碰撞、重着色和流动三个步骤组成，流场中混合流体分布函数定义为${f_i}\left( {{\boldsymbol{x}},t} \right) = \displaystyle \sum\limits_i {f_i^k}$（k=r或b），有

$\qquad {f}_{i}{}^{k}({\boldsymbol{x}} + {\boldsymbol{e}}_{i}\Delta t,t + \Delta t)={f}_{i}{}^{k}({\boldsymbol{x}},t)-{\varOmega }_{i}({f}_{i}{}^{k}({\boldsymbol{x}},t)) $

(1)

式中：x为粒子当前的坐标位置；$ {f_i}^k $为k相流体在位置x处第i个速度方向上的总分布函数；t为时间；Δt为时间步长；$ {{\boldsymbol{e}}_i} $为第i个方向上的格子离散速度；$ {\varOmega _i}\left( {f_i^k\left( {{\boldsymbol{x}},t} \right)} \right) $为碰撞项，由碰撞算子$\varOmega _i^1\left( {f_i^k\left( {{\boldsymbol{x}},t} \right)} \right)$、扰动算子$\varOmega _i^2\left( {f_i^k\left( {{\boldsymbol{x}},t} \right)} \right)$与重着色算子$\varOmega _i^3$组成。

$ {\varOmega _i}({f_i}^k({\boldsymbol{x}},t)) = \varOmega _i^3({\varOmega _i^1}({f_i}^k({\boldsymbol{x}},t)) + {\varOmega _i^2}({f_i}^k({\boldsymbol{x}},t))) $

(2)

平衡态分布函数$f_i^{{\rm{eq}}}\left( {{\boldsymbol{x}},t} \right)$表示为

$\qquad {f_i}^{{\rm{eq}}}(x,t) = {\omega _i}\rho \left[ {1 + \frac{{{{\boldsymbol{e}}_i} \cdot {\boldsymbol{u}}}}{{{c_{\rm{s}}}^2}} + \frac{{{{({{\boldsymbol{e}}_i} \cdot {\boldsymbol{u}})}^2}}}{{2{c_{\rm{s}}}^4}} - \frac{{{{\boldsymbol{u}}^2}}}{{2{c_{\rm{s}}}^2}}} \right]$

(3)

式中：$ {\omega _i} $为权重系数；$\rho = \displaystyle \sum\limits_k {{\rho _k}}$为流体总密度，$ {\rho _k} $为k相流体的密度；${c_{\rm{s}}} = c/\sqrt 3$为声速；c=dx/dt 为格子速度；局部速度由$\rho {\boldsymbol{u}} = \displaystyle \sum\limits_k {\displaystyle \sum\limits_i {f_i^k{{\boldsymbol{e}}_i}} }$来确定。

$\qquad {\omega _i} = \left\{ \begin{gathered} \frac{4}{9},\;\;i = 0 \\ \frac{1}{9},\;\;i = 1,2,3,4 \\ \frac{1}{36},\;\;i = 5,6,7,8 \\ \end{gathered} \right. $

(4)

$\qquad {{\boldsymbol{e}}_i} = \left\{ {\begin{array}{*{20}{l}} {\begin{array}{*{20}{l}} {(0,0),}&{i = {\rm{0}}} \end{array}}\\ {\begin{array}{*{20}{l}} {c\left( {{\rm{cos}}\left[ {\dfrac{\left( i-1 \right){\text{π}} }{2}} \right],{\rm{sin}}\left[ {\dfrac{\left( i-1 \right){\text{π}} }{2}} \right]} \right),}&\\{i = 1,2,3,4} \end{array}}\\ {\begin{array}{*{20}{l}} {\sqrt {\rm{2}} c\left( {{\rm{cos}}\left[ {\dfrac{\left( 2i-9 \right){\text{π}} }{4}} \right],{\rm{sin}}\left[ {\dfrac{\left( 2i-9 \right){\text{π}} }{4}} \right]} \right),}&\\{i = {\rm{5}},{\rm{6}},{\rm{7}},{\rm{8}}} \end{array}} \end{array}} \right. $

(5)

在Bhatnagar−Grass−Krook(BGK)碰撞算子中，一个计算时间步长内，粒子由一个位置结点运动到相邻的结点，并在该结点上与其他流体粒子发生碰撞，即

$\qquad {\varOmega _i^1}({f_i}^k({\boldsymbol{x}},t)) = \left( { - \frac{1}{\tau }} \right)\left[ {{f_i}^k({\boldsymbol{x}},t) - {f_i}^{k({\rm{eq}})}({\boldsymbol{x}},t)} \right] $

(6)

式中，$\tau $为弛豫时间。

考虑到在不互溶流体混合区域不同流体的黏度不相等，通过调和平均值来得到混合流体的黏度，即

$\qquad \frac{1}{\nu } = \frac{1}{\rho }\displaystyle \sum\limits_k {\frac{{{\rho _k}}}{{{\nu _k}}}} $

(7)

式中：$\nu $为混合流体运动黏度；${{\nu }_{k}} $为k相流体的运动黏度。

k相流体的弛豫时间$\tau $与其运动黏度有关。

$\qquad \tau = \frac{{{\nu _k}}}{{{c_{\rm{s}}}^2\Delta t}} + 0.5 $

(8)

Liu等^［13］推导出两相模拟中微扰算子的表达式，并采用连续表面力的概念来表示界面张力以及质量守恒和动量守恒的约束。

$\qquad \varOmega _i^2({f_i}^k({\boldsymbol{x}},t)) = \frac{{{A_k}}}{2}\left| {\boldsymbol{G}} \right|\left[ {{\omega _i}\frac{{{{({{\boldsymbol{e}}_i} \cdot {\boldsymbol{G}})}^2}}}{{{{\left| {\boldsymbol{G}} \right|}^2}}} - {B_i}} \right] $

(9)

式中：$ {B_i} $为流体界面特性参数，B₀=−4/27，B₁~B₄均为2/27，B₅~B₈均为5/108；G为颜色梯度；A_k为控制界面张力σ参数。

梯度算子计算式为

$\qquad {\boldsymbol{G}} = \frac{{{\rho _{\rm{r}}}}}{\rho }\nabla \frac{{{\rho _{\rm{b}}}}}{\rho } - \frac{{{\rho _{\rm{b}}}}}{\rho }\nabla \frac{{{\rho _{\rm{r}}}}}{\rho } $

(10)

式中：$\nabla $为梯度算子；ρ_r为红色流体密度；ρ_b为蓝色流体密度。

$\qquad \nabla \varphi \left( {{\boldsymbol{x}},t} \right) = \frac{1}{{c_{\rm{s}}^2{\Delta_{\rm{t}}}}}\sum\limits_i {{\omega _i}\varphi \left( {{\boldsymbol{x}} + {{\boldsymbol{e}}_i}{\Delta_{\rm{t}}}} \right)} {{\boldsymbol{e}}_i} $

(11)

σ的表达形式为

$\qquad \sigma = \frac{1}{9}\left( {{A_{\rm{r}}} + {A_{\rm{b}}}} \right){\tau} $

(12)

为了促进异相的分离并维持相界面的稳定，采用重着色方法减少不同流体的混合。在D'Ortona等^［14］的开创性工作基础上，Latva−Kokko等^［15］提出了一些改进，各分布函数定义为

$ \Omega _i^3({f_i}^k({\boldsymbol{x}},t)) = \frac{{{\rho _k}}}{\rho }{f^ * }\left( {{\boldsymbol{x}},t} \right) + \sum\limits_{l,l \ne k} {{\beta _0}{\omega _i}\frac{{{\rho _k}{\rho _l}}}{\rho }} {{\boldsymbol{n}}_{kl}} \cdot {{\boldsymbol{e}}_i} $

(13)

式中：$ {f^*}\left( {{\boldsymbol{x}},t} \right) $为通过扰动步计算后得到的新的混合流体的分布函数；ρ_l为l相流体的密度；n_kl为界面处的法向向量；β₀为异相界面厚度有关的分离系数。

2 物理模型及验证 2.1 物理模型

本文对剪切流动中纯液滴的变形及破裂过程进行了研究，同时对液滴中示踪粒子的运动进行了讨论。图1为液滴初始分布和示踪粒子示意图，其中：H为上、下边界距离；L为左、右边界距离。红色液滴浸没在蓝色流体中，初始状态时，液滴被置于整个流场的中心，流场的上、下边界以大小相同方向相反的速度U移动。流场的L=400，H=200，产生的剪切速率为$\dot \gamma $=2U/H，在该剪切流场下，液滴会产生形变甚至破碎。整个流场的水平方向上应用周期性边界条件，同时在上、下壁面应用Zou−He速度边界条件^［16］。同时，在水平方向左侧红、蓝相界面处放置一个示踪粒子跟踪界面处流体的运动情况。

液滴在流场中受到的剪切作用为黏性力、惯性力和表面张力的共同作用。定义无量纲时间${t^ * } = t\sigma /({\mu _k}\sqrt A )$，其中：k为b、r时分别表示蓝色、红色流体的动力黏度；A为液滴的面积。影响液滴变形的因素可以用雷诺数、毛细数、流体黏度比等无量纲数来描述，各参数表达式分别为：毛细数$Ca = \dot \gamma {\mu _{k}}R/\sigma$，R为液滴半径，雷诺数$Re = \dot \gamma \rho {R^2}/{\mu _{k}}$，流体黏度比$\lambda = {\mu _{{\rm{r}}}}/{\mu _{\rm{b}}}$。采用变形参数D和偏转角度θ定量表征液滴运动情况，$D = \left( {{L_1} - {B_1}} \right)/ \left( {{L_1} + {B_1}} \right)$，其中：L₁、B₁分别为液滴长轴和短轴的长度；θ为液滴变形后长轴与y方向的夹角。需要指出的是，若无特别说明，本文中所有物理量均取格子单位。图2为液滴变形参数示意图。

2.2 Young−Laplace定律验证

首先，为验证模型方法的正确性，模拟了不同半径红色液滴悬浮在蓝色流体中的情况。计算区域选取200 × 200，两种流体初始分布情况为

$ \left\{ {\begin{array}{l}{\rho }_{{\rm{r}}}=1,{\rho }_{{\rm{b}}}=0\qquad{\left(x-100\right)}^{2} + {\left(y-100\right)}^{2}\leqslant {R}^{2}\\ {\rho }_{{\rm{b}}}=1,{\rho }_{{\rm{r}}}=0\qquad\qquad\qquad\qquad{\text{其他}}\end{array}} \right. $

(14)

界面张力和流体黏度分别取σ=0.01和ν=0.1，并且在x和y方向上使用周期性边界条件。根据杨拉普拉斯定律，当系统达到平衡状态时穿过界面的压差∆p与界面张力σ的关系为

$\qquad \Delta p = \frac{\sigma }{R} $

(15)

图3为Young−Laplace定律验证结果。由图中可知，相界面的内外压差与液滴半径的倒数1/R成正比，其斜率为0.010 1，与本文中给定的界面张力（0.01）一致，符合Young−Laplace定律。

2.3 网格无关性验证

本文采用网格密度为200 × 400进行模拟，并选取100 × 200和400 × 800两种网格密度与其进行对比验证，计算中除网格密度外其他参数均一致。结果如表1所示，其中ε为误差。在三种不同的网格密度下，液滴在剪切流场中变形参数和偏转角度的差别均很小。

3 结果与讨论 3.1 毛细数Ca对液滴变形的影响

首先，在Re=0.8、R/H=1/2、λ=1下，研究毛细数对剪切流中液滴变形的影响，结果如图4～5所示。将毛细数限制在0.05～0.55内时可以控制液滴的变形程度且使其不会发生破裂。从图中可以看到，随着毛细数的变化，液滴的形态发生改变，液滴的变形参数和偏转角度相应发生变化。在剪切力作用下，液滴从静止开始沿剪切力作用方向不断拉伸和倾斜，逐渐偏离椭圆形，并在达到最大形变时保持稳定。当Ca较小时，液滴受到的黏性剪切力不足以抵抗表面张力，液滴的变形程度较小；当Ca不断增大时，黏性剪切力对液滴的影响开始占主导地位，液滴的变形程度加大。同时，从图5(b)中也可以发现，液滴的变形参数随Ca近似线性变化。

相界面处示踪粒子在不同Ca下的运动路程和剪切作用如图6所示，其中：S为示踪粒子运动路程；Δu_x为示踪粒子速度在x方向上的分量与其所处位置的剪切速度之差，可反映剪切流场对示踪粒子产生的作用。图6(a)中，低Ca下示踪粒子运动速度较小，其运动路程也较小。随着Ca增加，示踪粒子运动速度增大，运动路程增加；当Ca增大到一定程度时，示踪粒子的变化趋于稳定。图6(b)中，当Ca较小时示踪粒子受剪切作用的影响较大，速度变化较大，当Ca增大时示踪粒子受剪切作用的影响以及相应速度的变化均不大。

3.2 雷诺数Re对液滴变形的影响

在Ca=0.35、R/H=1/2、λ=1下，雷诺数Re对剪切流动中液滴变形的影响如图7～8所示，其中将Re控制在0.2～1.2之间。低Re下黏性力占主导地位，液滴不易发生变形，随着Re增大，液滴的变形参数随之增大；Re增大意味着液滴受到的惯性力越大，液滴在水平方向受到的剪切作用增强，液滴的变形效果更加强烈。从图8中可以看出，随着Re增加，液滴的偏转角度变化幅度较小，说明液滴的偏转角度已处于一个较稳定的范围。

图9为Re对示踪粒子运动路程和剪切作用的影响。图9(a)中界面处示踪粒子的运动路程呈逐渐增大的趋势。图9(b)中，Re越大，示踪粒子水平方向速度与剪切速度的差值越大；随着Re的增大，示踪粒子受到的惯性力增大，示踪粒子的运动速度和运动路程均增大。在Re增大时，液滴的变形也随之加剧。

3.3 红蓝相黏度比λ对液滴变形的影响

在Ca=0.35、Re=0.8、R/H=1/2下，红蓝相黏度比λ对液滴动力学行为的影响以及对液滴变形参数和偏转角度的影响分别如图10、11所示。随着λ的增加，液滴的变形参数先增大后减小，其偏转角度先增大后保持稳定。这说明当λ较小时，λ的增大会使得施加在液滴上的黏性剪切力增加，使液滴在剪切作用下的变形加剧，从而往水平方向逐渐偏转；当λ较大时，液滴相对于外部连续相表现得更像固体，流场对其施加剪切力的作用效果较小，液滴变形效果降低，其朝水平方向的偏转程度变化不大。

图12为λ对示踪粒子运动路程和剪切力作用的影响。从图中可以看出，当λ大于0.25时，不同λ对界面处示踪粒子的运动影响较小，改变液滴的黏性而保持流场区域流体黏性不变对示踪粒子的驱动力作用不明显，示踪粒子受到的剪切作用效果变化不大，其运动速度和运动路程的变化均较小。

3.4 液滴变形和破裂的Re−Ca相图

图13为剪切流场中液滴不同变形和破裂的Re−Ca相图。可以看到，当将Ca控制在0.05～0.35时，液滴会发生类型1变形，且在变形稳定后其内部会形成一个涡；此时增大Re或Ca会使液滴的变形加剧，并逐渐被拉伸成哑铃状，内部形成两个涡；如果再进一步增大Re或Ca会使得液滴达到变形的极限，哑铃状液滴的中间连接部分被拉断，从而破裂形成两个独立的液滴。

3.5 液滴不同位置示踪粒子的运动情况

在Ca=0.35、Re=0.8、R/H=1/2、λ=1下，示踪粒子与液滴圆心初始距离对其运动路程以及剪切作用的影响。如图14所示，其中p为示踪粒子初始位置距离液滴圆心的距离。图14(a)、(b)分别反映了液滴发生变形生成单个涡和两个涡时的情况。离液滴圆心越远的示踪粒子运动路程越大，直到示踪粒子接近界面时，由于表面张力的存在它受到来自界面的约束作用增强，运动路程才有所减少。

从图14（c）可以看到，当液滴内部生成一个涡时，初始位置离液滴圆心越近的示踪粒子受到的剪切作用越小，其速度越小，运动路程也越短。图14（d）中液滴在变形过程中产生两个涡，离圆心较近的示踪粒子会在变形后向哑铃状液滴的一侧运动，距离较远的示踪粒子则会在整个液滴区域内作绕流运动，同样离圆心距离越远的示踪粒子受到的剪切作用较小。

图15(a)为单个涡中示踪粒子与界面距离Q随时间的变化。在变形过程中示踪粒子与界面的距离越来越小，而后保持稳定。图15（b）为在两个涡中示踪粒子与界面距离随时间的变化，相比于图15（a），示踪粒子最终与界面的距离更小，此时液滴被拉伸成哑铃状，总体呈细长状态。

4 结　论

应用格子玻尔兹曼颜色梯度模型研究了液滴在二维剪切流动中的动力学行为。讨论了雷诺数Re、毛细数Ca和流体黏度比对复合液滴拉伸变形的影响和示踪粒子的运动规律，对液滴变形和破裂的Re−Ca相图进行分析，得到以下结论：

（1）在剪切流场中，随着Ca的增大，液滴逐渐拉伸为椭圆形，且变形程度逐渐增大，示踪粒子的运动路程随之增大，其受到剪切作用的影响先增大后趋于稳定。

（2）随着Re的增加，液滴的变形参数增大，偏转角度变化不大，示踪粒子的运动路程增大，其受到的剪切作用也越大。

（3）随着λ的增加，液滴的变形参数先增大后减小，偏转角度先增大后保持稳定，示踪粒子随黏度比的变化，其运动路程和受到剪切作用的变化变化均不大。

（4）在Re−Ca相图中可以看到两种不同的变形类型和液滴破裂形态的分布，随着Re与Ca的增大，液滴形态由变形转为开始破裂。

（5）离液滴圆心初始距离越近（离界面越远）的示踪粒子运动路程越大，受剪切作用的影响越大。