JFM：载有颗粒的界面在胶体液滴剪切诱导变形中的作用

2.1 问题描述

如图1所示，我们考虑一个孤立的载胶体液滴悬浮在一个恒定剪切速率γ的简单剪切流中。初始半径为R的球形液滴位于长方体计算域的中心。液滴周围流体和液滴本身都是不可压缩的牛顿流体，具有相同的密度ρ，它们的粘度分别为μout和μin。在长方体区域内的流动剖面被上下运动板块保持为稳定的线性剪切流动，而液滴在时间框架t = 0时突然释放。值得注意的是，液滴承载着一组总数目为N、半径为a的胶体粒子，它们在t = 0时随机但均匀地悬浮在液滴内部或分布在液滴表面。由于布朗扩散，自由悬浮的颗粒可能靠近液滴表面，然后被液滴表面吸收。被吸附的颗粒也可能因颗粒间的排斥相互作用而被迫分离(Gu & Botto 2018)。在这项研究中，特别感兴趣的是揭示颗粒负载界面在剪切诱导的胶体液滴变形中的作用。

2.2 两相流的水平集方法

采用水平集法(Sethian & Smereka 2003)捕获液滴表面并跟踪液滴的变形。通常，使用一组连续性和Navier-Stokes方程来解决滴液内外的流场。

这里， u 是速度向量；t是时间；p是压力。除非另有规定，所有方程都通过按R缩放长度、按γR缩放速度、按1/γ缩放时间、按ργ²R²缩放压力和按ργ²R缩放体力项来无量纲化。在(2.2)中，μ(x)在滴液外部和内部分别是1和λμ = μin/μout; Re是雷诺数，定义为ργR²/μout; 体力f由滴液表面张力fσ和悬浮粒子施加的力fp引起。滴液表面表示为水平集函数φ的零等值线，这是到滴液表面的有符号距离函数。因此，滴液表面根据水平集方程（Gibou, Fedkiw & Osher 2018）移动：

这里，n是滴液表面的向外单位法向量。由于φ在求解(2.3)后不保留有符号距离函数的性质，因此使用重新初始化方程将水平集函数φ0转换为有符号距离函数φ：

这里，τ是虚拟时间，sign是平滑的符号函数。

为了避免数值不稳定性，通过以下Heaviside函数对液滴表面流体粘度的急剧变化（即(2.2)中的μ(x)）进行平滑处理：

在这里，通常将半平滑宽度ε设置为求解流场网格尺寸的1.5倍，通常认为这是流体表面的半厚度。计算液滴表面向外单位法向n和平均曲率κ

然后，由于液滴表面张力而产生的体积力fσ可以按照以下方式计算：

式中，Ca为净毛细数，定义为μoutγ R/σc， σc为净表面张力;δε(φ) =∂Hε(φ)/∂φ是光滑的狄拉克函数。

2.3 粒子运动的控制方程

拉格朗日点粒子法用于跟踪粒子运动，该方法已广泛应用于分散多相流的高效模拟(Balachandar& Eaton 2010). & Eaton 2010)。根据牛顿第二定律，粒子运动的控制方程为

为了便于理解，该方程中的所有变量都是有量纲的。在这里，ρp是粒子密度；up是粒子速度矢量；g是重力加速度；Fdrag是由于粒子-流体相互作用而产生的流体动力阻力；Fpp代表粒子间相互作用力；Fip是由于粒子-界面相互作用而产生的附着力。

为简单起见，本研究中省略了粒子的自旋，因此通过Stokes阻力模型考虑粒子-流体相互作用，如下:

这里，uf是质点中心的流体速度;Fd为阻力系数，对于自由悬浮在散装流体中的颗粒，通常为1。对于浸入流体界面的颗粒，当液滴黏度与周围流体不同时，其阻力系数fd变化显著。根据先前的分析研究(Pozrikidis 2007;Dörr 等人 2016)，我们将阻力系数fd修改如下:

这里，d是粒子中心到

流体界面之间的距离，对于悬浮在液滴内部的粒子来说，d是负值。

其次，简要讨论了惯性力项和重力、浮力项。将式(2.8)变换为无因次形式，将等式两边除以μoutγ R2

所有变量都是无量纲的，但为了简单起见，符号保持不变。其中，λρp = ρp/ρ为粒子与流体的密度比;ω = a/R为颗粒-液滴尺寸比;Fr为弗鲁德数，定义为γ 2R/|g|。在本研究中，λρp， ω， Re和Fr分别被设置为近似<10,10−2,10−1和100。因此，与水动力阻力项相比，(2.11)中的惯性力项和重力浮力项可以忽略不计，即它们的比值小于10−4。

那么，粒子运动的控制方程可以进一步简化为:

其中rp是粒子位置向量。

粒子布朗运动通过随机游走模型考虑，该模型可以建模为Wiener过程(Ma 等人 2012)。因此，得出了以下随机微分方程，作为粒子运动控制方程的最终形式：

根据Gu & Botto(2016, 2018, 2020)提出的技术，通过对浸没在流体界面中的颗粒施加附着力Fip来有效地考虑颗粒-界面相互作用，其计算公式为

这里，d是从粒子中心到液滴表面的最小距离的矢量。在该模型中，粒子表面的物理和化学性质对粒子-界面相互作用的影响被包含在无量纲参数fip中。

这里，θp是粒子与流体界面的接触角。

对于粒子-粒子相互作用，本研究考虑了排除的体积相互作用、粒子碰撞和流体动力润滑，因此(2.8)中的力Fpp实际上由以下三部分组成:

遵循先前研究中广泛使用的技术(Ma 等人 2012;Liu 等人 2016;Zhao & Yong 2017, 2018)，摩尔斯势的排斥部分用于两个粒子之间的排除体积相互作用，吸引部分用于弱吸引相互作用(如果有的话)。量纲形式的摩尔斯电势由下式给出。

在这里，r是两粒子之间的距离；εrep和εatt分别是排斥力和吸引力的Morse相互作用强度；λm衡量的是相互作用势发生显著变化时的长度尺度；rm是Morse相互作用力变为零时的特定距离。为避免粒子重叠，采用赫兹势来描述两粒子之间的接触力（Frijters 等人2012; Xie & Harting 2018; Vowinckel 等人2019），其维数形式为：

在这里，εh是赫兹相互作用强度，rh表示接触力变为零时的特定距离。因此，由于Morse势引起的粒子间相互作用力Fpp,m可以用无量纲形式表示如下：

在此方程中，无量纲参数定义如下:ε∗re rep = εrep/(σca2)， ε∗att = εatt/(σca2)， λ∗m λm = λm/a, r∗m m = rm/a。考虑粒子碰撞时，由赫兹势引起的接触力Fpp,h可由下式给出：

这里，ε∗h∗h = εha1/2/σc, r∗h r∗h = rh/a。在下面的讨论中，为方便起见，去掉了表示无量纲参数的上标“*”。值得注意的是，当两个接近的颗粒之间的距离小于流体网格尺寸时，需要一个润滑模型来适当考虑它们之间的水动力相互作用。与两个粒子的相对速度ui - uj相反的法向润滑力Fpp,l由((Frijters 等人2012; Xie & Harting 2018; Vowinckel 等人2019)下式给出：

这里ε是一个临界长度，它远小于流体网格尺寸he。当粒子间间隙r−2a小于ε时，润滑力保持恒定以避免数值不稳定(Vowinckel 等人 2019)。流体动力润滑力在很大程度上取决于颗粒间隙内流体的粘度，因此系数flu的形式与式(2.10)中的fd相同。值得注意的是，在我们目前的方法中省略了切向润滑力，因为它通常比正常的要小得多，正如以前对颗粒悬浮液的数值研究所指出的那样（Lambert 等人2013; Suzuki & Hayakawa 2019)。

2.4 模拟粒子对液滴表面张力的影响

值得注意的是，吸附颗粒完全不影响流体界面的固有张力，但颗粒间的排斥性相互作用产生正表面压力，从而降低了有效表面张力(manikanantan & Squires 2020)。在我们目前的方法中，粒子-流体相互作用是以双向耦合的方式考虑的，因此考虑了粒子对流场的反向影响。为此，在动量守恒方程中加入了流体动力阻力Fdrag的反作用力和作用在粒子上的附着力Fip。

因此，(2.2)中粒子施加的物体力项f p可以用无因次形式给出：

其中，δ(rf−rp)为光滑的狄拉克函数(Luo, Shang & Bai 2019a)，其中rf和rp分别代表流体网格点和粒子的位置。Gu & Botto(2020)已经证实，通过在(2.2)中加入体力项fp即(2.22)，可以正确地再现粒子间相互作用对粒子负载界面有效表面张力的影响。

2.5 数值实现

控制方程是通过我们自己的三维有限差分代码在一个规则并置的均匀笛卡尔网格上求解的。采用四步投影法求解流体速度和压力解耦。采用二阶中心差分格式对(2.1)和(2.2)中的所有空间导数进行离散化。对于时间离散，扩散项采用了Crank-Nicolson半隐式技术，对流和体力项采用了三级龙格-库塔技术。为了降低计算成本，压力泊松方程只在最后一个龙格-库塔阶段求解，并采用四层多重网格技术和交替方向隐式格式进一步提高计算效率。流求解器在我们之前的工作中有详细的讨论(Luo & Bai 2018)。对于水平集平流和重初始化方程(2.3)和(2.4)，分别采用五阶加权基本非振荡格式和三阶总变差递减龙格-库塔格式进行空间和时间离散化。为了推进胶体颗粒的位移，采用两阶段龙格-库塔法求解颗粒运动的控制方程，即(2.13)。有关数值实现的更多细节见于附录A。

2.6 无量纲参数总结

通过对所有控制方程的总结，我们发现了15个可能影响剪切流动中胶体液滴动力学的无量纲参数。表1总结了本研究中使用的控制无量纲参数及其值。值得注意的是，总粒子数N通常用粒子在液滴表面的表面覆盖率φ = Nω2/4或液滴内部的体积粒子浓度cnp = Nω3来代替。

参考以前的实验研究(Mulligan & Rothstein 2011;Mei et al . 2016;Kaganyuk & Mohraz 2020)，雷诺数通常远小于1，因此惯性效应可以忽略不计。根据我们之前的数值研究(Luo, He & Bai 2015;Luo et al . (2019b)， Re = 0.1足够低，可以消除惯性效应，而进一步降低Re会显著增加计算时间。如§2.3所述，由于Re和ω都很小，λρ和Fr的影响可以忽略不计。在我们的模拟中使用的几个重要参数Ca = 0.05-0.3， ω = 0.01-0.04， θ p = 45°-90°，Pep = 102-104， ϕ = 0-0.4和cnp = 0-0.05的值涵盖了以往实验中提出的典型范围。先前的实验无法获得与粒子间相互作用相关的参数，但在粒子悬浮液的几项数值研究中采用了莫尔斯和赫兹势(Liu 等人 2016;Gu & Botto 2018, 2020;Xie & Harting 2018;Zhao & Yong 2018)。在本研究中，我们发现εh = 100和rh = 0.01足以防止粒子重叠。另外，λm和rm参照以往数值研究取4，εrep = 0.05 ~ 0.5表示粒子间排斥性相互作用对液滴动力学的影响。注意，在接下来的讨论中，为了便于理解，εrep被εm代替。

表1。本研究中使用的无量纲参数及其值。在第3.2节中，对于颗粒覆盖的液滴，εm， ϕ和Ca是变化的，而其他参数则是固定的，除非另有说明。在第3.3节中，对于含有胶体颗粒的液滴，只有Pep、θ p和cnp是变化的。

3.1 验证：净液滴的变形

在我们之前对剪切流动中胶囊或液滴动力学的研究中，流动求解器已经得到了充分的验证，其中使用前缘跟踪方法进行界面跟踪(Luo 等人 2015;Luo & Bai 2018;Luo 等人 2019b)。在这里，我们通过模拟简单剪切流中清洁液滴的变形，进一步验证了水平集函数的求解器。根据我们之前工作的验证，2.5 ×10−4的时间步长，10R×6R×10R的计算域大小和120 × 72 × 120的欧拉分辨率足以消除它们对液滴动力学的影响，因此我们在本研究中使用了它们。实际上，当欧拉分辨率从120 × 72 × 120提高到200 × 120 × 200时，颗粒覆盖液滴的变形变化被限制在1%以下。与以往研究类似，液滴动力学主要由Taylor参数表征:(i)变形指数D = (L−B)/(L + B)，其中L和B分别为变形椭球状液滴的半长轴和半短轴;(ii)变形液滴长轴与流动方向的倾角θ。

通过预测简单剪切流中洁净液滴的变形，将我们的水平集方法与Kennedy, Pozrikidis & Skalak(1994)的边界积分方法、Komrakova等人(2014)的点阵玻尔兹曼方法和我们自己的锋面跟踪方法(Luo 等人2019b)进行了定量比较。如图2所示，稳态变形指数D和倾角θ是毛细管数Ca的函数。从定性的角度来看，变形液滴的三维形状轮廓类似于椭球体，这与之前的研究观察结果一致。从定量的角度来看，我们的方法和Komrakova等人的方法预测的液滴变形略高于Kennedy等人的方法，这可能是由于雷诺数（Reynolds number）的差异造成的，即我们方法中Re=0.1，而Kennedy等人的方法中Re=0。值得注意的是，我们的方法与Komrakova等人方法之间的最大相对误差限制在4.5%以下。而我们自己的水平集（level-set）方法和前沿追踪（front tracking）方法之间的相对误差更小，小于2.5%，这两种方法都采用了相同的流动求解器。此外，导致液滴破碎的临界毛细管数Cacr大致等于我们水平集方法预测的0.375，这与Kennedy等人的结果非常接近，即Cacr=0.37。

图2。通过简单剪切流中清洁液滴变形的建模验证我们的水平集方法。将稳态变形指数D (a)和倾角θ (b)与之前报告的结果进行比较，包括我们自己使用前沿跟踪方法的工作(Luo等人2019b)， Komrakova等人(2014)使用晶格玻尔兹曼方法和Kennedy等人(1994)使用边界积分方法。图(c)中绘制了四个示例的稳态形状轮廓。

3.2 颗粒覆盖的液滴变形

粒子引起的液滴表面有效张力的变化被认为是胶体粒子对整个液滴变形产生潜在影响的主要原因。因此，在本小节中，我们将重点关注简单剪切流中颗粒覆盖液滴的变形。所有颗粒初始均以随机而均匀的模式排列在液滴表面，没有颗粒自由悬浮在液滴体相中。此外，通过设置相当大的强度，即fip = 10 in(2.14)，为粒子-界面相互作用产生的附着力，省去了粒子分离和再附着的动力学。

我们首先通过计算在静止流体中颗粒覆盖液滴表面上的压力跃变，来研究吸附颗粒如何影响有效表面张力。如图3(a)所示，当达到稳态时，液滴表面上的吸附颗粒排列成有序的六边形微结构，这主要是由于颗粒间的排斥相互作用。一般来说，如图3(b)所示，吸附颗粒会导致整个液滴的有效表面张力σeff降低。这种有效的降低主要归因于颗粒间的排斥作用，它会产生一个二维表面压力来抵消流体界面张力。我们的结果与Gu和Botto（2016）之前研究中的报告结果在性质上相似。由于颗粒间排斥力的有限范围，因此需要最小颗粒覆盖率ϕmin才能对有效表面张力产生非零影响。然而，当ϕ>ϕmin时，σeff以近似指数形式减小，因为颗粒间的排斥作用变得显著。此外，随着εm的增大，σeff的减小更加显著，因为颗粒间的排斥作用得到了增强。

图3. 通过计算悬浮在静止流体中被颗粒覆盖的液滴表面的压差，研究了吸附颗粒对有效表面张力的影响。被吸附的颗粒最初以随机、均匀的模式分布在整个液滴表面。(a)绘制了稳定状态下颗粒单层的微观结构，用于不同的颗粒覆盖φ。绿色填充点表示吸附颗粒。(b)在不同粒子间相互作用强度εm时，负载粒子界面的有效表面张力σeff是φ的函数。

图4. 剪切流中被颗粒覆盖的液滴的典型例子。最初，所有的颗粒均匀地覆盖在液滴表面，不允许它们分离。对于不同的颗粒表面覆盖度，即图(a)和(b)中的ϕ = 0.2和0.4，给出了不同时间框架下变形液滴的三维形状轮廓，绿色填充点表示吸附颗粒。蓝色标记的点用于表示吸附颗粒沿液滴表面的旋转运动。绘制了(c)变形指数D和(D)倾角θ在不同φ下的时间演化图。其中，毛细管数Ca固定为0.15，粒子间相互作用强度εm固定为1。

在图4中，我们给出了几个典型的剪切流动中颗粒覆盖的液滴的例子，其中绘制了不同颗粒表面覆盖下的变形指数D和倾角θ的时间演变图，即φ = 0, 0.2, 0.3和0.4。与干净的液滴类似，被颗粒覆盖的液滴也会变形为椭球形状，并最终保持这种稳定的形状，其表面不断旋转。然而，在整个瞬态时间内，颗粒覆盖的液滴通常比清洁的液滴表现出更大的变形，并且在颗粒表面覆盖率较高时，液滴变形能力的增强更为明显。因此，颗粒覆盖的液滴沿流动方向呈现较小的倾斜角。此外，被颗粒覆盖的液滴达到稳态的总时间也比清洁的液滴要长得多，并且随着颗粒表面覆盖率的增加，达到稳态的时间会进一步延长，这与变形能力越强的液滴更偏向于流动的预期是一致的。

上述几个例子表明，吸附颗粒通常会促进液滴的变形。接下来，我们进行了系统的模拟，以揭示含颗粒界面在剪切作用下对颗粒覆盖液滴变形的影响。为了直接说明这一点，图5展示了在不同颗粒间相互作用强度εm和颗粒表面覆盖率ϕ组合下，变形液滴在稳态下的三维形状轮廓。一般来说，颗粒覆盖的液滴呈近椭圆形。值得注意的是，图中还展示了液滴表面上吸附颗粒的分布情况。原则上，由于颗粒间的排斥作用，吸附颗粒不能聚集在一起，因此可以在整个液滴表面上自由移动。因此，吸附颗粒的分布可能看起来不均匀。在所有展示的情况下，颗粒在液滴尖端附近的堆积密度都高于其他表面积。实际上，尽管吸附颗粒会沿着液滴表面连续旋转，但这种特定的分布不会随时间变化（数据未显示）。这种特定的分布主要是由表面对流引起的，表面对流也被认为是表面活性剂分子类似分布的原因（Luo et al，2019b）。当然，随着ϕ的增加，颗粒分布的这种不均匀性会明显减少，因为更强的颗粒间排斥作用使吸附颗粒的自由移动性降低。

图5。在Ca = 0.1的条件下，对颗粒间相互作用强度εm和颗粒表面覆盖φ的不同组合，绘制了剪切流中颗粒覆盖液滴的稳态形状。彩色轮廓表示颗粒间距离dr，以表示吸附颗粒在液滴表面上的不均匀分布，每个颗粒(例如标记为i)计算为最接近颗粒i的六个颗粒的平均dr。

图6.剪切流中颗粒覆盖液滴的稳态变形。(a)对于不同颗粒间相互作用强度εm，稳定变形指数D作为颗粒表面覆盖φ的函数，Ca固定为0.1。(b) εm = 0.5时，不同毛细数Ca下的D vs φ。

接下来，我们从定量的角度研究吸附颗粒对液滴变形的影响。为此，我们绘制了稳态变形指数D随颗粒表面覆盖率ϕ变化的图，如图6所示，图中考虑了不同的颗粒间排斥强度εm或不同的毛细管数Ca。吸附颗粒确实会导致D和θ发生显著变化，这在固定Ca时可以明显观察到。一般来说，液滴的变形能力显著提高，这最初可以归因于颗粒引起的有效表面张力的降低（如图3所示）。更具体地说，D随ϕ显著增加，且在更大的εm下增加得更快。例如，当ϕ从0增加到0.4时，在εm分别为0.1和1的情况下，D分别从0.12增加到0.13和0.43。实际上，我们还通过增加颗粒与液滴的尺寸比ω（从0.01增加到0.04）来研究了其影响，但发现只要ϕ和εm固定，ω的影响就可以忽略不计（数据未显示）。值得注意的是，这一结果并不意味着颗粒尺寸a对液滴变形没有影响，实际上，它已经被包含在颗粒表面覆盖率ϕ = Na2/4R2和无量纲颗粒间相互作用强度εm = εrep/(σca2)的影响中。总之，吸附颗粒通常会增加液滴的变形能力，并且这种效果会随着颗粒间排斥相互作用的增强而显著增强。

图7. 剪切流动中颗粒覆盖液滴的稳态三维形状。(a)在Ca = 0.1条件下，不同粒子间相互作用强度εm时，涡度方向上的液滴宽度W/R作为粒子表面覆盖φ的函数;(b) W/R随变形指数d绘制，此处使用与图6相同的参数值。实线为我们之前在表面活性剂负载液滴研究中发现的W/R = 1−D2曲线(Luo et al, 2019b)。(c)液滴表面积A与D的关系，实心曲线代表了L/R = 1/(1 - D)、W/R = 1 - D2、B/R = 1/(1 + D)三个主轴的椭球体。(D) W/R = 1 - D2曲线与前人实验数据吻合良好(Mei et al .2016)。

简单剪切流涡度方向的液滴变形对于掌握变形液滴的完整三维形状也很重要，因此我们分析了图7中的稳定宽度W/R。显然，被吸附的颗粒也显著增加了涡度方向的液滴变形，W/R明显减小。这种W/R的降低也通过增加φ或εm得到加强。在我们之前对剪切流中清洁液滴和表面活性剂负载液滴的研究(Luo et al, 2019b)中，我们发现了意想不到的W/R与D之间的关系，即W/R = 1−D2。如图7(b)所示，被颗粒覆盖的液滴的所有数据点也坍缩到这条曲线上。此外，如图7(c)所示，液滴表面积A与L/R = 1/(1−D)、W/R = 1−D2和B/R = 1/(1 + D)这三个主轴椭球的表面积A也非常吻合。此外，我们将W/R = 1−D2曲线与图7(D)中报道的颗粒覆盖液滴的实验数据(Mei et al 2016)进行了比较，两者也非常吻合。结果表明，在简单剪切流中，被颗粒覆盖的液滴仍保持着L/R = 1/(1−D)、W/R = 1−D2和B/R = 1/(1 + D)三个主轴的完美椭球形状。

图8. 稳定变形D是有效毛细数Caeff = μoutγ R/σeff的函数，其中σeff为加载颗粒界面的有效表面张力，如图3所示。图(a)中，σeff的取值为:ϕeff = Nπa2/ A，其中A为变形液滴的表面积。图(b)中，Caeff−max = μoutγ R/σeff−min，其中σeff−min为颗粒堆积最密集的水滴尖端附近出现的最小σeff。实线表示清洁液滴的变形，用于比较。值得注意的是，对于同一符号，随着D的增加，颗粒表面覆盖率从0.05增加到0.4，即ϕ= 0.05,0.1,0.15,0.2,0.25,0.3,0.35和0.4。

颗粒引起的有效表面张力的降低是解释剪切流动中颗粒覆盖液滴变形扩大的一个明显而重要的机制，但它可能不是唯一的机制。为了验证这一假设，我们在图8(a)中绘制了颗粒覆盖的液滴的稳定变形D作为有效毛细数Caeff = μoutγ R/σeff的函数，其中σeff通过使用图3中的数据与给定的ϕ和εm值来评估。显然，大多数数据点显示出与清洁液滴曲线的正偏差，特别是在高颗粒覆盖率时。注意，图8(a)中的Caeff表示整个液滴表面σeff的平均值，它是用平均粒子覆盖率计算得出的，即:ϕeff = Nπa2/ A(A为变形液滴的表面积)。考虑到颗粒分布的明显不均匀性，我们在图8(b)中绘制了D与基于液滴表面σeff最小值计算的Caeff−max的关系图，由于颗粒堆积最密集，通常出现在液滴尖端。对于干净的液滴，大多数数据点都低于曲线，但少数数据点(即Ca = 0.2和ϕ = 0.4)仍然高于曲线。由此可以得出结论，除了整个液滴表面有效表面张力的整体降低外，吸附颗粒的不均匀分布是导致颗粒变形能力增加的另一个重要原因。更具体地说，与有效表面张力均匀恒定的液滴相比，不均匀的颗粒分布使得液滴末端的表面张力较低，而在液滴赤道处的表面张力较高。因此，为了满足法向粘性应力平衡，液滴末端和赤道的形状需要更尖锐和更平坦，最终导致整个液滴的变形更大。实际上，在之前的研究中也观察到类似的现象(Stone & Leal 1990;Li & Pozrikidis 1997)对于含有表面活性剂的液滴具有极高的psamet数，即O(1000)，而在低psamet数，即O(0.1) ~ O(10)时，这种现象不明显，甚至消失。

值得注意的是，在图8中，我们观察到了一些令人惊讶的现象，这些现象与颗粒覆盖液滴在相对较高颗粒覆盖率下的变形有关。首先，导致液滴破裂的有效毛细管数的临界阈值通常低于干净液滴的临界阈值。其次，颗粒覆盖的液滴可能会表现出极大的变形，即D > 0.65，这在干净液滴中并未观察到，即D < 0.6。第三，在颗粒表面覆盖率较高（例如ϕ > 0.35）的情况下，有效毛细管数Caeff不再随颗粒表面覆盖率的增加而增加，但D在高覆盖率下仍然显著增加。控制参数ϕ = Nω2/4是根据最初球形液滴的表面积计算的，它实际上仅取决于颗粒数N（在ω固定的情况下）。然而，Caeff是根据有效颗粒覆盖率ϕeff计算的，它除了取决于颗粒数N外，还取决于变形液滴的表面积A。因此，对于在高ϕ下高度变形的液滴，ϕeff可能不再随N的增加而增加，因为表面积A也在增加。

我们还分析了吸附颗粒对液滴表面速度Us大小的影响，如图9所示。显然，Us在相位角方向上表现出很大的变化，尤其是在中间剪切平面上。一般来说，Us在液滴赤道附近达到最大值，在液滴尖端附近达到最小值。液滴表面上的这种非恒定速度会引起表面对流效应，这是液滴尖端颗粒堆积更密集的最根本原因。基本上，Us取决于液滴的形状，如图9(b)和9(c)所示，通过降低最大表面速度Us-max和最大与最小表面速度之差Us-max - Us-min来表示。更重要的是，吸附颗粒显著降低了液滴的表面流动性；即，即使通过将它们与D作图来排除液滴形状的影响，Us-max和Us-max - Us-min都会减小。这种降低的表面流动性与含有表面活性剂的液滴所观察到的相似，它是由切向马兰戈尼应力引起的（Li & Pozrikidis 1997）。颗粒引起的Us减小主要归因于颗粒阻塞效应（Luo & Bai 2020）。因此，随着粒子间斥力的增加，Us-max−Us-min的变化幅度减小，因为粒子干扰随着εm或φ的增加而增强。反过来，由于表面对流效应减弱，表面迁移率降低导致颗粒分布的不均匀性降低(图6)。值得注意的是，液滴表面附近的局部剪切速率随着表面迁移率的降低而增加，如图9(a)所示。因此，周围流体作用于液滴表面的剪切力可能会被放大。液滴尖端上增加的剪切力可能是将液滴拖到更长的形状的重要因素。因此，颗粒表面迁移率的降低可能是另一个潜在的机制，即使在相同的毛细数下，颗粒覆盖的液滴也比干净的液滴表现出更高的变形能力(图8)。

图9.吸附颗粒对液滴表面附近流场和液滴表面流场的影响。(a)给出了液滴表面速度大小Us的三维轮廓和中间剪切平面局部剪切速率∂u/∂z +∂w/∂x的二维轮廓。(b)最大表面速度Us-max和(c)最大和最小表面速度之差Us-max - Us-min作为液滴变形指数d的函数绘制，参数值与图8相同。

3.3 胶体液滴变形

图10. 剪切流中胶体液滴的示例。在初始时刻，所有颗粒都在液滴内部随机、均匀地悬浮，且没有颗粒设置在液滴表面上。在图(a–c)中，展示了不同时间框架下的液滴形状轮廓，其中绿色实心点代表吸附在液滴表面的颗粒，分别对应于不同的颗粒Péclet数，即Pep = 104、103和102。蓝色标记点表示吸附在液滴表面上的几个颗粒的旋转运动。在图(a)中的三个示例中，绘制了(d)液滴变形D和(e)吸附颗粒总数Nnp-s随时间的变化。此处，Ca = 0.1, εm = 0.5, ω = 0.02, θp = 90° 和 cnp = 0.05。

在实际应用中，胶体颗粒通常悬浮在体积相中，体积相与液滴表面之间的颗粒交换可能以脱落和再附着的形式发生。在本节中，我们进一步考虑了剪切流中胶体液滴的变形，这种液滴由一个包含颗粒悬浮液的颗粒负载表面组成。在图10中，我们展示了不同颗粒Péclet数（即Pep = 104、103和102）下液滴变形D和浸入液滴表面的颗粒总数Nnp-s随时间演变的典型例子。一般来说，这些胶体液滴也会变形成椭球形。对于Pep = 104，颗粒的布朗扩散相当弱，因此一旦颗粒吸附到液滴表面，就很少会脱落。此外，如图10(a)中标出的颗粒所示，由于表面对流，吸附的颗粒显然会连续旋转。在Pep = 103时，吸附的颗粒仍然很少脱落，但它们不再随着液滴表面连续旋转。相比之下，由于扩散效应的增加，吸附的颗粒可能会在液滴表面以随机模式移动。随着Pep进一步降低到102，布朗运动导致吸附的颗粒在界面上波动，甚至频繁脱落和再附着。因此，对于Pep = 102，Nnp-s表现出明显的振荡，但对于Pep = 103和104，Nnp-s则缓慢且连续地增加，没有振荡。因此，在较大的Péclet数（即Pep = 104）下，胶体液滴表现出与颗粒覆盖液滴相似的变形行为。然而，在Pep = 102时，胶体液滴的变形比具有相同数量吸附颗粒的颗粒覆盖液滴更大。

图11. 液滴内部不同颗粒浓度下的胶体液滴变形。(a) 对于不同的颗粒浓度cnp，绘制了吸附颗粒总数Nnp-s随时间的变化图，其中Pep = 103，θp = 90°。(b) 吸附颗粒沿界面厚度的数密度分布图，其中Nnp是浸没在厚度为0.1a的薄流体层中的吸附颗粒数，这些颗粒与液滴表面的平均距离为d。注意，负d表示液滴内部的位置。(c) 稳态下的液滴形状轮廓图，其中绿色和红色实心点分别表示吸附在液滴表面上的颗粒和在体积流体中自由悬浮的颗粒。

显然，液滴表面吸附的颗粒表面浓度取决于液滴内部的体积颗粒浓度。在图11中，展示了不同体积颗粒浓度（即cnp = 0.01–0.05）下Nnp-s随时间的变化，其中颗粒Péclet数Pep固定为103。对于所有给定的cnp值，自由悬浮的颗粒会扩散到液滴表面附近并吸附在液滴表面上，因此Nnp-s随时间连续增加并最终达到平稳状态。因此，随着吸附颗粒数量的增加，液滴变形D也持续增加（数据未显示）。值得注意的是，在所有时间框架内，Nnp-s均随cnp的增加而显著增加。此外，图11(b)展示了不同cnp或不同Peb下吸附颗粒沿界面厚度的数密度分布。由于液滴表面的附着力，吸附颗粒的数密度沿垂直于界面向两侧延伸的法线方向逐渐减小。随着cnp的减小，吸附颗粒的数密度分布变化不大，但总体上有所减小。更有趣的是，随着Peb的减小，数密度分布变得更加平坦，这表明由于明显的扩散效应，吸附颗粒的脱落变得更加频繁。值得注意的是，界面附近体积流体中自由悬浮的颗粒需要克服来自吸附颗粒的排斥相互作用力，才能最终吸附到液滴表面上。当Peb从103增加到104时，尽管颗粒脱落减弱，但吸附颗粒的总数反而减少（图10e）。

图12. 颗粒表面润湿性对剪切流中胶体液滴变形的影响。对于不同的颗粒接触角，即θp = 90°、60°和45°，在Pep = 103和cnp = 0.04的条件下，绘制了(a)液滴变形D和(b)吸附颗粒总数Nnp-s随时间的变化图。在(c)和(d)图中，绘制了不同颗粒接触角（即θp = 90°、60°和45°）下，cnp对稳态时D和Nnp-s的影响。图(c)中的虚线表示Ca = 0.1时纯净液滴的变形。

除了载有颗粒的界面外，自由悬浮的颗粒也可能影响液滴变形，因为它们可能会增加体相的有效粘度。在图12中，我们研究了由颗粒接触角θp表征的颗粒表面润湿性对液滴变形的影响。随着θp的减小，沉浸在液滴表面的颗粒总数Nnp-s显著减少，因为防止颗粒脱落的附着力减弱，如(2.14)式所示。因此，无论是瞬态还是稳态，液滴变形D都随着θp的减小而显著减小。值得注意的是，在颗粒接触角较小的情况下，即θp = 45°时，液滴表面的颗粒覆盖率低于0.05，因此载有颗粒的界面对液滴变形的影响可忽略不计。在这种情况下，发现体相中的颗粒浓度对液滴变形的影响很小。这是因为本研究中使用的自由悬浮颗粒的体相浓度，即cnp = 0–0.05，太低而无法对体相的有效粘度产生明显影响。更具体地说，正如Guazzelli和Pouliquen（2018）所指出的，当cnp低于约0.05时，可以使用爱因斯坦的线性依赖定律很好地预测稀颗粒悬浮液的有效粘度，即μeff = 1 + 2.5cnp。因此，在我们的模拟中，cnp = 0–0.05时，颗粒引起的内部到外部粘度比的有效增加被限制在0.125以下，根据Rallison（1984）的预测D = (19λμ + 16)/(16λμ + 16)Ca，这只能引起液滴变形的很小变化（即<1％）。因此，对于本研究中使用的参数值，载有颗粒的界面对剪切流中颗粒引起的液滴变形变化起主要作用。然而，在未来的工作中，探索载有颗粒的界面与体相中浓密颗粒之间的耦合相互作用对胶体液滴变形的影响将是非常有趣的，因为浓密悬浮液可能会导致有效流体粘度的显著增加，甚至产生其他复杂的非牛顿流变学特性（Guazzelli & Pouliquen 2018; Morris 2020）。

在这项工作中，我们开发了一种数值方法来模拟胶体液滴（即载有胶体颗粒的液滴）的流动动力学。我们特别关注许多颗粒对液滴动力学的集体效应，采用点粒子法来追踪颗粒运动，这已被证明是在可实现的数值精度和可承受的计算资源之间取得平衡的一个绝佳选择（Balachandar & Eaton 2010）。已经开发了基于格子玻尔兹曼方法（Usta 等人2009；Zhao & Yong，2017，2018；Xie & Harting，2018）的几种方法，用于模拟含有内部胶体颗粒的液滴，但尚未详细讨论吸附在液滴表面上的颗粒及其对表面张力的影响。Gu和Botto（2016，2018，2020）开发了一种用于载有颗粒的流体界面的快速数值方法，其中有效处理了颗粒-界面相互作用及其对有效界面张力的影响。然而，现有的点粒子方法没有考虑悬浮颗粒的物理体积和布朗运动。我们目前的方法进一步考虑了这些特殊效应，因此，它可能成为研究流动条件下胶体液滴动力学的有用工具。值得注意的是，与以往研究中使用的格子玻尔兹曼方法和相场方法不同，我们采用了水平集方法来捕获流体界面，因为它在计算局部界面特性（如平均曲率和点到界面的距离向量）方面具有优势，这些特性对于计算（2.28）中的体力项fpe至关重要。

通过系统模拟，我们证明了在低体积分数(即本研究中使用的cnp < 0.05)下自由悬浮在体相中的胶体颗粒对液滴变形的影响相当有限，这与以往的实验和数值数据一致(Usta 等人 2009;Ammel et al 2020)。这一现象可以解释为，在低颗粒浓度下，颗粒引起的附加粘性应力相当小。然而，相同数量的颗粒一旦吸附在液滴表面并保持在液滴表面，则可能对液滴变形产生显著影响。对于表面润湿性为中性的颗粒，即颗粒接触角θp固定为90°时，我们的结果表明了这种效应。此外，在典型的高粒子psamclet数下，即Pep = 104，粒子一旦吸附到液滴表面就会表现出罕见的分离。特别是由于弱布朗扩散，吸附颗粒的数量增加缓慢，这解释了之前实验中观察到的液滴变形缓慢增加的原因(Mei 等人2016)。这些结果表明，只要颗粒倾向于停留在液滴表面，颗粒负载界面对整个液滴的变形起着至关重要的作用。实际上，研究粒子布朗运动对液滴变形的影响将是非常有趣的，特别是在小的psamclet数(例如Pep ~ O(1))下。沿液滴表面的强扩散和频繁的分离和再附着可能使颗粒分布更加均匀，从而可能抑制D-Caeff曲线与干净液滴的偏离。然而，离散随机漫步模型的时间步长随着Pep成比例地减小，因此求解需要大量的计算时间(2.13)。显然，为了进一步研究Pep低至0(1)时胶体液滴的变形，迫切需要一种高速并行方案来提高数值实现的效率。

为了将我们的注意力集中在载有颗粒的界面所起的作用上，我们对简单剪切流中颗粒覆盖液滴的变形进行了系统研究。特别令人感兴趣的是确定颗粒间相互作用和颗粒表面覆盖率对液滴变形的详细作用。一般来说，吸附的颗粒会促进液滴变形，这与之前关于低颗粒覆盖率液滴的实验数据一致（Mei 等人 2016；Kaganyuk & Mohraz，2020）。在这些实验中，还发现当颗粒接近完全覆盖时，吸附的颗粒会抑制液滴变形，而这一点在我们的当前研究中并未研究。Frijters 等人（2012）的数值模拟也预测了吸附颗粒会增加液滴变形，颗粒覆盖率范围为0–0.55。这种增加效果归因于在液滴尖端附近积聚的吸附颗粒产生的额外惯性力。通过设置较低的雷诺数，我们的研究展示了颗粒增强液滴变形性的完全不同的机制，即颗粒引起的有效表面张力降低。事实上，在Mei et al.（2016）和Kaganyuk & Mohraz（2020）的实验中，液滴的雷诺数均限制在10−1以下。因此，我们认为颗粒降低的宏观表面张力可能是先前实验中观察到的具有相对较低颗粒覆盖率的液滴变形增加的真正机制。

此外，通过绘制液滴变形与有效毛细管数的关系图，我们发现，即使在相同的毛细管数下，颗粒覆盖的液滴仍然表现出比干净液滴更大的变形。因此，颗粒引起的整体表面张力降低不能完全解释载有颗粒的界面在颗粒覆盖液滴变形中的促进作用。我们讨论了其他两种潜在机制，即液滴表面颗粒分布的不均匀性和颗粒引起的液滴表面迁移率降低。首先，有效毛细管数实际上是根据整个液滴表面颗粒覆盖率的平均值来评估的。然而，由于表面对流，吸附的颗粒通常表现出明显的不均匀分布。特别是，颗粒通常在液滴尖端处表现出更密集的堆积，这在Kaganyuk & Mohraz（2020）对部分涂层液滴的实验中也观察到了。显然，使用有效毛细管数无法正确描述这种不均匀颗粒分布对液滴变形的影响。其次，液滴表面速度被吸附的颗粒降低，特别是在液滴尖端附近颗粒堆积更密集的区域。我们之前对微通道中颗粒覆盖液滴的研究（Luo & Bai 2020）也证明了表面迁移率的降低。由于表面迁移率的降低，作用在液滴表面上的周围剪切力变得更大，这可能会将两个液滴尖端向相反方向拖动，从而对增加液滴变形产生潜在贡献。

我们的数值方法中仍有几个悬而未决的问题。首先，由于三重接触线的存在，胶体颗粒在流体界面上的布朗扩散与在散装流体中的布朗扩散有很大不同，因此吸附颗粒沿界面或垂直界面的扩散系数可能不同(Boniello 等人2015;Dörr 等人 2016)。因此，在一些实际应用中，需要考虑吸附颗粒的界面扩散。其次，粒子间相互作用模型是另一个值得关注的问题。我们目前的方法包括流体动力润滑和接触力，这通常需要粒子悬浮液的点粒子方法(Frijters 等人 2012;Xie & Harting 2018;Vowinckel et al 2019)。我们还考虑了由于颗粒表面的聚合物涂层而产生的额外空间或重叠排斥，这在实际系统中经常用于稳定胶体颗粒(Garbin 等人 2015;Israelachvili 2015)。然而，许多其他种类的相互作用力(特别是短程相互作用力)需要进一步仔细考虑，例如静电、范德华和毛细相互作用(Deshmukh et al. 2015;Sharifi-Mood, Liu & Stebe 2015)，这对于密集堆积的颗粒单层尤其重要。还应该注意的是，通过修改胶体相互作用的经典理论来定量描述这些粒子间力可能在纳米尺度上失败(Batista, Larson & Kotov 2015)。此外，对于不同类型的粒子，用于粒子间相互作用的参数可能有很大不同的值，但在实验中很难测量(Kim 等人 2019b)。第三，前人的实验研究(Mei et al . 2016;Kaganyuk & Mohraz(2020)提出假设，粒子之间的弱吸引相互作用可以通过表面弹性来解释抑制液滴变形。此外，研究发现，改变吸引-排斥相互作用的相对大小可以在具有密集堆积颗粒的流体界面上诱导复杂的粘弹性(Rahman等，2019)。更重要的是，吸引粒子间的相互作用赋予了形成界面裂纹的可能性，即颗粒单层断裂(Kim et al . 2019a;Chai 等人 2020)，在此之后，新的界面区域暴露出来，因此流体界面呈现出更不均匀的表面张力。显然，这些复杂的界面特性可能会导致颗粒负载的液滴产生更丰富的动力学，特别是在极端变形的情况下(Garbin 2019)，这可以通过使用我们的方法通过修改(2.19)中的吸引粒子间相互作用来揭示。

附录A 数值实现

一旦流体速度通过求解(2.1)和式(2.2)得到更新，水平集函数就可以通过求解(2.3)和式(2.4)得到平流。水平集平流方程和重初始化方程在空间上通过五阶加权基本非振荡格式离散化，在时间上通过三阶总变分递减龙格-库塔格式离散化。更多细节可以在Gibou等人(2018)最近的综述中找到。众所周知，在多相流模拟的水平集方法中存在质量损失。一个重要的原因是重新初始化过程中的数值误差。我们使用Sussman & Fatemi(1999)提出的技术在重新初始化过程后对φ进行校正，如下所示:

其中λ_ijk是流体网格单元ijk的常数，计算公式为:

其中Ωijk表示[(x, y,z)|xi−1/2 < x < xi+1/2, yj−1/2 < y < yj+1/2, zk−1/2 < z < zk+1/2]的定义域。

虽然在重初始化过程中可以保持质量守恒，但由于求解水平集平流方程时的数值误差，仍然存在质量损失。在我们目前的方法中，每个时间步长的质量损失被限制在10−6左右。在本研究中大多数模拟中，总时间步长约为5 × 104，少数情况下可能达到5 × 105。如果在每个时间步不进行校正，液滴体积在最终稳定状态的总损失可能高达5% - 50%。为了进一步消除液滴体积中的这种误差，我们使用Ge等人(2018)提出的接口校正方法。通过沿液滴表面向外法线方向移动零能级集，将液滴体积修正为初始值，通过求解下式实现:

其中δV为液滴在时间间隔δt内的体积变化，A为液滴的表面积。采用五阶加权基本非振荡格式对对流项进行离散化。通过求解(A3)对每个时间步的液滴体积进行校正，最终将所有模拟的液滴体积变化限制在10−3左右。值得注意的是，在每个时间步长的液滴体积修正过程中，液滴表面的位移实际上很小，在10−5a左右，因此这种对液滴体积的人工修正对吸附颗粒的运动影响很小。

对于粒子运动控制方程即(2.13)的求解，采用两阶段龙格-库塔法:

其中，△t为时间步长;N和N +1是时间步长数;N +1/2表示中间变量。粒子位置处的流体速度uf通过狄拉克函数从相邻的流体网格u插值，如下所示:

这里，狄拉克函数δ的离散形式由下式给出：

为了计算由于粒子-界面相互作用而产生的附着力Fip(即(2.14))，还通过狄拉克函数从水平集函数中插值出粒子与液滴表面之间的距离|d|和力方向d/|d|

在求解动量方程（即（2.2））时，需要强调的一件事是，颗粒引起的体力项fp和fpe（即（2.22）和（2.24））需要分布在多个网格上，以避免点源引起的数值不稳定性，这也可以通过使用离散狄拉克δ函数来实现，如（A7）所示。

（来源：多相流在线）