液滴在固体表面上的铺展现象广泛存在工业应用中，如打印机将墨水喷到纸面^[1]、农药喷洒在蔬菜表面^[2]、表面喷涂^[3]、喷雾冷却^[4]等。液滴撞击壁面大致会出现3种现象：沉积、溅射和反弹^[5]。很多工业应用涉及到液滴在壁面上的沉积问题^{[2-3, 6]}，为抑制液滴溅射和反弹来提高液滴在壁面上的沉积率，通常有两种途径^[6]：1)加入添加剂改变流体的流变特性；2)改变固壁表面的性质。第一种方法用途广泛，例如在液体中添加聚合物使液滴呈现非牛顿性，抑制液滴在壁面上的溅射和反弹，提高药物的使用效率^[2]；铁轨上的摩擦改进剂，通过水溶剂中添加聚合物复合材料改变流变特性，提高沉积率，以降低机车燃料消耗和轨轮的磨损^[3]等等。在流体中添加不同的成分得到的流变特性不同^[7-11]，例如添加弹性颗粒会促进液滴在壁面上的反弹^[11]，而在牛顿流体中添加高分子聚合物会使液体产生幂律流变特性^{[9, 12-14]}。幂律流体在控制液体沉积问题中有着重要的应用，目前对于幂律流体的动态湿润理论研究虽然已经取得一些进展^[14-16]，但数值模拟的研究仍然较少^[15]。

与传统数值方法相比，格子Boltzmann方法作为一种介观的数值方法，凭借其独特的优势在研究非牛顿流体方面得到了广泛的关注^[17-21]。在格子Boltzmann方法中，最常用的是单松弛模型，由于该模型的数值稳定性较差，笔者给出了非牛顿流体的多松弛格子Boltzmann模型。目前，研究两相流界面运动常有的数值方法有水平集法、流体体积函数法(VOF)^[22]、相场法^[23]等，相对于其他方法，基于自由能理论的相场法在研究多组分系统中介观结构的演化等复杂现象时具有特有的优势，如多相流、接触线以及涉及复杂界面演化和相变问题；与此同时，相场法既有水平集、VOF等方法的捕捉大变形界面的能力，又无需涉及界面的重构过程，相界面的位置捕捉以及厚度大小的控制均可通过一个对流扩散的控制方程来实现，基于该方法追踪得到的两相界面较光滑且连续性好。因此，相场法是一种可以有效模拟液滴撞击固壁问题的方法。

笔者以格子Boltzmann方法为基础，基于相场理论并借助Cahn-Hilliard方程，建立了非牛顿幂律流体与牛顿流体的两相流模型，模拟了假塑性流体在不同幂律指数n(0.5≤n≤1.0)和韦伯数We(5≤We≤45)下液滴撞击壁面的铺展，该模拟结果对于通过改变流变特性控制液滴的沉积提供了参考。

1 理论 1.1 幂律(Power-law)模型

对于不可压缩流体，切应力张量τ与变形速度张量$ \tilde {\mathit{\boldsymbol{S}}} $之间的关系为：$ {\mathit{\pmb{\tau}}} = 2\eta \tilde {\mathit{\boldsymbol{S}}} $，其中$ \tilde {\mathit{\boldsymbol{S}}} = [\nabla {\mathit{\boldsymbol{u}}} + {(\nabla {\mathit{\boldsymbol{u}}})^{\rm{T}}}]/2 $，η为动力学黏性系数。对于幂律模型^[24]，黏性系数为剪切速率$ {\dot \gamma } $的函数，即

$ \eta (\dot \gamma ) = {\eta _0}|\dot \gamma |n - 1, n > 0, $

(1)

式中：η₀为渐进黏度；$ \dot \gamma = 2\sqrt {{D_\prod }} $，二维情况下$ {D_\prod } = \sum\nolimits_{\alpha , \beta = 1}^2 {{{\tilde S}_{\alpha \beta }}} {{\tilde S}_{\alpha \beta }} $，是应变速率张量$ \tilde {\mathit{\boldsymbol{S}}} $的第二不变量；n为幂律指数，当n>1时，式(1)描述的是膨胀性流体，也称剪切增稠流，例如玉米面糊就属于膨胀性流体；当n=1时，则是牛顿流体；当n＜1时，则是假塑性流体，也称剪切稀化流，包括大部分高分子聚合物的溶液。

1.2 Call-Hilliard相场模型

Call-Hilliard对流扩散方程^[25-27]：

$ \left\{ \begin{array}{l} \frac{{\partial C}}{{\partial t}} + {\mathit{\boldsymbol{u}}} \cdot \nabla C = M{\nabla ^2}{\mu _{\rm{c}}}, \\ {\mu _{\rm{c}}} = {C^3} - 1.5{C^2} + 0.5C - {\varepsilon ^2}{\nabla ^2}C, \end{array} \right. $

(2)

式中：C表示体积分数，u是速度。移动系数M=M₀|C(1－C)|，μ_c是化学势。界面的厚度取为$ D = 4\sqrt 2 \varepsilon $。引入静态接触角θ^eq后，在固壁处的边界条件为^[28]：$ {{\mathit{\boldsymbol{n}}}_{\rm{w}}} \cdot \nabla {\mu _{\rm{c}}}{|_{\rm{s}}} = 0 $和$ {{\mathit{\boldsymbol{n}}}_{\rm{w}}} \cdot \nabla C = - |\nabla C|\cos {\theta ^{{\rm{eq}}}} $，n_w为固壁上的单位法向量。在本文中，参数均取为$ \varepsilon = 1/\sqrt 2 $，M₀=0.2。

对于等温两相流流动的控制方程^[25-26]：

$ \nabla \cdot {\mathit{\boldsymbol{u}}} = 0, $

(3)

$ \rho (\frac{{\partial {\mathit{\boldsymbol{u}}}}}{{\partial t}} + {\mathit{\boldsymbol{u}}} \cdot \nabla {\mathit{\boldsymbol{u}}}) = - \nabla p + \nabla \cdot [\eta (\nabla {\mathit{\boldsymbol{u}}} + \nabla {{\mathit{\boldsymbol{u}}}^{\rm{T}}})] + {{\mathit{\boldsymbol{F}}}_{\rm{s}}}, $

(4)

式中：密度ρ=Cρ₁+(1－C)ρ₂，黏度η=Cη₁+(1－C)η₂，下标1和2代表不同的组分。p表示压力，$ {{\mathit{\boldsymbol{F}}}_{\rm{s}}} = \frac{{6\sqrt 2 \sigma }}{\varepsilon }\mu \nabla C $为体积表面张力，这里σ是表面张力。

1.3 两相流格子Boltzmann方法

下面采用格子Boltzmann方法^[25-27]来求解上述两相流流动的控制方程(3)和(4)。对于D2Q9模型，为了得到压力和速度，引入了一个变量$ {g_\alpha } = {f_\alpha }c_{\rm{s}}^{\rm{2}} + (p - \rho c_{\rm{s}}^{\rm{2}}){\Gamma _\alpha }(0) $，其中c_s是声速，f_α是粒子分布函数，Γ_α(u)是一个跟宏观速度以及f_α的平衡态f_α^eq有关的函数：

$ f_\alpha ^{{\rm{eq}}} = \rho {\Gamma _\alpha }({\mathit{\boldsymbol{u}}}) = \rho {\omega _\alpha }[1 + \frac{{{{\mathit{\boldsymbol{e}}}_\alpha } \cdot {\mathit{\boldsymbol{u}}}}}{{c_{\rm{s}}^{\rm{2}}}} + \frac{{{{({{\mathit{\boldsymbol{e}}}_\alpha } \cdot {\mathit{\boldsymbol{u}}})}^2}}}{{2c_{\rm{s}}^4}} - \frac{{{\mathit{\boldsymbol{u}}} \cdot {\mathit{\boldsymbol{u}}}}}{{2c_{\rm{s}}^{\rm{2}}}}], $

(5)

式中e_α为α方向的微观粒子的速度; ω_α是相对应的权系数。马赫数Ma=|u|/c_s，在低马赫数的情况下，宏观的水动力方程可由式(6)g_α的离散Boltzmann方程恢复得到：

$ \frac{{\partial {g_\alpha }}}{{\partial t}} + {e_\alpha } \cdot \nabla {g_\alpha } = - {{\mathit{\pmb{\Lambda}}} _{\alpha \beta }}({g_\beta } - g_\beta ^{{\rm{eq}}}) + ({{\mathit{\boldsymbol{e}}}_\alpha } - {\mathit{\boldsymbol{u}}}) \cdot [\nabla \rho c_{\rm{s}}^{\rm{2}}({\Gamma _\alpha } - {\Gamma _\alpha }(0)) + {\mathit{\boldsymbol{F}}}{\Gamma _\alpha }], $

(6)

式中采用的是多松弛时间模型。Λ_αβ为粒子碰撞矩阵，分子间作用力为$ {\mathit{\boldsymbol{F}}} = \nabla \rho c_{\rm{s}}^{\rm{2}} - (\nabla p - {{\mathit{\boldsymbol{F}}}_{\rm{s}}}) $，其中g_α的平衡态函数：$ g_\alpha ^{{\rm{eq}}} = c_{\rm{s}}^{\rm{2}}f_\alpha ^{{\rm{eq}}} + (p - \rho c_{\rm{s}}^{\rm{2}}){\Gamma _\alpha }(0) $。对式(6)进行Crank-Nicolson积分，再引入一个新的粒子分布函数$ {{\bar g}_\alpha } $，其具体表达式为：$ {{\bar g}_\alpha } = g_\alpha ^{{\rm{eq}}} + 1/2{S_{\alpha \beta }}({g_\alpha } - g_\alpha ^{{\rm{eq}}}) - 1/2 \cdot \delta t({{\mathit{\boldsymbol{e}}}_\alpha } - {\mathit{\boldsymbol{u}}}) \cdot [\nabla \rho c_{\rm{s}}^{\rm{2}}({\Gamma _\alpha } - {\Gamma _\alpha }(0)) + {{\mathit{\boldsymbol{F}}}_{\rm{s}}}{\Gamma _\alpha }] $，式(6)转化为：

$ {{\bar g}_\alpha }({\mathit{\boldsymbol{x}}} + {{\mathit{\boldsymbol{e}}}_\alpha }\delta t, t + \delta t) - {{\bar g}_\alpha }({\mathit{\boldsymbol{x}}}, t) = - ({S_{\alpha \beta }} + 2{I_{\alpha \beta }})({{\bar g}_\beta } - \bar g_\beta ^{{\rm{eq}}}) + \\\;\delta t({{\mathit{\boldsymbol{e}}}_\alpha } - {\mathit{\boldsymbol{u}}}) \cdot [\nabla \rho c_{\rm{s}}^{\rm{2}}({\Gamma _\alpha } - {\Gamma _\alpha }(0)) + {{\mathit{\boldsymbol{F}}}_{\rm{s}}}{\Gamma _\alpha }], $

(7)

式中粒子碰撞矩阵S_αβ+2I_αβ= $ M_{\alpha \gamma }^{ - 1}{{\hat S}_{\gamma \varepsilon }}{M_{\varepsilon \beta }} $，其中M_αβ为正交转换矩阵。$ {{\hat {\mathit{\boldsymbol{S}}}}_{\alpha \beta }} = {M_{\alpha \gamma }}({S_{\gamma \lambda }} + 2{I_{\gamma \lambda }}){\mathit{\boldsymbol{M}}}_{\lambda \beta }^{ - 1} $是一个对角松驰矩阵，可以表达为$ {\hat {\mathit{\boldsymbol{S}}}} $=(s₀, s₁, s₂, s₃, s₄, s₅, s₆, s₇, s₈)。其中s₀，s₃，s₅表示守恒量(质量和动量)，通常取值为s₀=s₃=s₅；由于对称性的要求s₄=s₆，此外系数s₇=s₈，均和黏度有关。所以独立变量有3个s₁，s₂和s₄，这3个取值变化可以用来改善模型的稳定性。所以，经过比较参考文献[25]和[27]的数值，选择取值为$ {{\hat {\mathit{\boldsymbol{S}}}}_{\alpha \beta }} $=diag(1, 1, 1, 1, 1.7, 1, 1.7, s₇, s₈)。其中松弛参数s₇跟黏性系数η的关系为：$ \eta = \rho c_{\rm{s}}^{\rm{2}}(\frac{1}{{{s_7}}} - \frac{1}{2})\delta t $。

这样，宏观速度和压力可以通过$ {{\bar g}_\alpha } $的零阶矩和一阶矩得到：

$ \begin{array}{l} p = \sum\limits_{\alpha = 0}^8 {{{\bar g}_\alpha }} + \frac{{\delta t}}{2}{\mathit{\boldsymbol{u}}} \cdot \nabla \rho c_{\rm{s}}^{\rm{2}}, \\ {\mathit{\boldsymbol{u}}} = \frac{1}{{\rho c_{\rm{s}}^{\rm{2}}}}(\sum\limits_{\alpha = 0}^8 {{{\bar g}_\alpha }} {{\mathit{\boldsymbol{e}}}_\alpha } + \frac{{\delta t}}{2}{{\mathit{\boldsymbol{F}}}_{\rm{s}}})。\end{array} $

(8)

由于非牛顿流体的黏度与变形速度张量$ {\tilde {\mathit{\boldsymbol{S}}}} $有关，在每个当前时间步得到速度u后，利用关系式$ \tilde {\mathit{\boldsymbol{S}}} = [\nabla {\mathit{\boldsymbol{u}}} + {{(\nabla {\mathit{\boldsymbol{u}}}})^{\rm{T}}}]/2 $并采用差分格式离散求得变形速度张量$ {\tilde {\mathit{\boldsymbol{S}}}} $，通过$ {\tilde {\mathit{\boldsymbol{S}}}} $求得$ {\dot \gamma } $后，根据幂律流体的黏度与剪切速率$ {\dot \gamma } $的关系，可得出下一个时间步的黏性系数^{[17, 29]}。为了比较各种参数，如表面张力、黏性力、惯性力对液滴撞击固体表面的影响，对于该两相流流动涉及的量纲一的量有奥内佐格数(Ohnesorge number) $ Oh = \eta /\sqrt {\rho R\sigma } $和韦伯数(Weber number)We=ρU²R/σ(U为特征速度)，此外雷诺数$ {\mathit{Re}} = \sqrt {We} /Oh $。本文中所涉及的算例均在格子单位下进行的，时间步和空间步均设为1。

2 模型和程序的验证 2.1 方腔流

为了验证格子Boltzmann模型模拟非牛顿幂律流体的能力和程序的可靠性，将两相流程序中的体积分数C全部设置为1，模拟了单相幂律流体在幂律指数n=0.5时的方腔流流动。其中，Re=ρU^2－nLⁿ/η₀= 100, U=0.05, 计算域为[0, 200]×[0, 200]。将本文的数值结果与Neofytou的文献结果^[30]相比较，从图 1中看出，两者吻合。

2.2 Laplace定律

验证了基于相场法的格子Boltzmann模型模拟幂律流体、牛顿流体两相流的能力和程序的可靠性。根据Laplace定律，液滴放置在母液中，在未受其他外力的作用下，液滴达到稳定状态时，液滴内外压差与其半径以及表面张力系数存在一定的关系。对二维情况有：Δp=p_in－p_out=σ/R，其中R为液滴半径，p_in是液滴的压力，p_out是液滴外母液的压力。初始时刻液滴位于计算区域中心，计算为[0, 6R]×[0, 6R]，四边均为周期性边界条件。参数n=0.5，ρ₁=1.0，ρ₂=1.0，η₁=η₂=0.1。下标1表示液滴，下标2表示液滴周围的液体，以下所涉及算例的下标含义均是如此。图 2中，在2种表面张力参数σ=1.0 E－3及σ=2.0 E－3下，数值解和理论解的相对误差均在2.53 %以内，且随着半径R的增大，相对误差随之减小。最大伪势速度在液滴界面附近，大小为10^－7量级。

2.3 液滴在壁面上的动态接触

验证了在引入考虑接触角影响的边界条件后的两相流模型模拟液滴动态接触的能力和程序的正确性。首先，考察了固壁上液滴达到平衡态时，静态接触角是否跟设置值一致。初始时刻，在固壁上放置一个半圆形液滴，由于表面张力的作用，液滴会变形，直至达到平衡状态。计算域[0, 200]×[0, 80]，上下边界为完全反弹边界条件，左右边界施加周期性边界条件。液滴半径R=30，其他参数为n=0.5，σ=10^－3，η₁=0.01，η₂=0.001和ρ₁=ρ₂=0.1。平衡态接触角是固壁附近处界面两层网格线性拟合的直线以顺时针旋转到壁面的夹角。图 3为液滴在不同静态接触角下达到平衡状态后的形状。从表 1可以看出，数值模拟得到的平衡态接触角与理论值的相对误差的最大绝对值为3.3 %(即静态接触角为45°)，稳态时的液滴面积与初始液滴面积的相对误差的最大绝对值为2.98 %。

考察了移动接触线的变化规律。初始时刻，液滴刚接触固壁，在表面张力作用下液滴会变形，接触线也随之变化。液滴半径R=20，设液滴与壁面上的接触线半径为r，则量纲一的铺展长度为r/R，时间采用$ \sqrt {\rho {R^3}/\sigma } $处理得到量纲一的时间t^*。计算域[0, 200]×[0, 80]，幂律指数n=1.0，静态接触角θ^eq=45°。从图 4可以看出，在液滴铺展的初始阶段r/R与(t^*)^1/2成正比，这一结论与文献[31]相一致，验证了所建模型的可靠性。

3 液滴撞击固壁上的铺展

图 5为单液滴以一定速度u₀撞击固壁，D是液滴的铺展长度，R为液滴半径。以下算例中均不考虑重力效应, 且黏度比均为10。上下两边施加完全反弹边界条件，左右两边界为周期性边界。计算域[0, 200]×[0, 200]，其他参数R=20，θ^eq=90°，σ=10^－3，η₁=0.01，η₂=0.001，ρ₁=ρ₂=0.1。韦伯数We=ρ₁u₀²R/σ，奥内佐格数$ Oh = {\eta _1}/\sqrt {{\rho _1}R\sigma } $。

首先，模拟了在We=45及Oh=0.158时幂律指数n=0.5~1.0，液滴撞击壁面上的铺展。从图 6可知，对于幂律指数n≥0.7的液滴，撞击壁面后大致可分为3个阶段：快速铺展、回缩、稳定。在快速铺展阶段，液滴迅速地铺展到最大长度，幂律指数越大，最大铺展长度越大；在回缩阶段，液滴开始由最大铺展长度向中心收缩，幂律指数越大，回缩的速度越大；最后，液滴回缩至其铺展长度达到一个恒定值，与回缩阶段相对应，液滴回缩的速度越快，越能更快达到稳定状态。对于幂律指数n=0.5的液滴，铺展过程中无回缩现象，液滴的铺展速度由初始最大值逐渐减小，当液滴铺展到最大铺展长度时，液滴即到达稳定状态。

从整体上看，对于幂律指数n≤1.0的液滴，最大铺展长度随着幂律指数n的减小而减小，且幂律指数越大，液滴越趋近于牛顿流体；幂律指数越小，液滴的非牛顿特性就越明显，表现为回缩现象不明显，达到稳态所需要的时间也变长。这些结果表明，在碰撞速度相同下，与牛顿流体相比，假塑性流体更易黏附在壁面，且n越小，流体非牛顿特性越强。

图 7为选取上述算例中幂律指数n=0.5以及n=1.0(牛顿流体)时液滴撞击壁面后形态的变化。对于牛顿液滴，初期在惯性的作用下迅速铺展，液滴两端凸起直至铺展到最大铺展长度，然后在以表面张力为主导的作用下，液滴向中心收缩，当液滴的动能全部转化为内能时，液滴最终以半圆的形态静止在壁面上。相对于牛顿流体，幂律指数n=0.5的假塑性液滴在撞击壁面后，黏性迅速增大导致液滴的铺展、变形受到抑制，液滴的拓扑结构没有发生较大变化，仍然大致呈一个圆形；铺展后期，液滴的能量耗散速度变慢，经过较长一段时间后，液滴铺展到最大铺展长度时，并以半圆形的形态静止在壁面上，液滴达到稳态。上述分析表明，相对于牛顿液滴，假塑性流体的液滴撞击壁面后的形态更稳定。

然后，模拟了在Oh=0.158下幂律指数分别为n=0.5以及n=1.0时，韦伯数We=5~45液滴撞击壁面上的铺展。从图 8中可以看出，对于n=0.5的假塑性幂律流体，液滴铺展到最大铺展长度时液滴即到达稳态，稳态跟韦伯数无关。当韦伯数等于5和20时，液滴均相对较缓慢地到达稳态，到达稳态的所需时间均较长；当韦伯数继续增大至45时，液滴铺展到最大铺展长度所需要的时间急剧减小。对于n=1.0的牛顿流体，液滴的稳态也与韦伯数无关，且随着速度的增大，液滴到达最大铺展长度的所需时间变长，最大铺展长度也增大；液滴到达稳态的时间随速度增加而增长但不显著。此外，假塑性幂律流体和牛顿流体在稳态时的铺展长度均相同，碰撞速度相等时前者到达稳态的时间要比后者长，随着速度的增大，这种差异会急剧减小，例如，在We=45时假塑性液滴和牛顿液滴到达稳态的时间已基本上相同，且稳定后的铺展长度均为54.5格子长度。以上结果均表明，在高碰撞速度下，与牛顿流体相比，假塑性流体会抑制液滴的铺展，更有利于液滴在壁面上的沉积。

4 结论

在两相流格子Boltzmann模型中引入流体的非牛顿特性，建立了基于相场法的牛顿流体、非牛顿幂律流体两相流格子Boltzmann模型，并在不考虑重力作用的条件下，对不同幂律指数和不同韦伯数下液滴撞击壁面的铺展进行了数值模拟。当n＜1时幂律流体的非牛顿特性会抑制液滴的铺展，且n越小，幂律流体表现的非牛顿特性就越强；韦伯数We越大，液滴到达稳定的速度越快。这些结果表明，假塑性流体的非牛顿特性有利于液滴在壁面上的沉积。