在模拟具有复杂界面形状和剧烈界面变形的流动时，传统的基于网格的方法存在网格畸变和数值扩散等问题，计算的准确性和稳定性难以保证。粒子法因各相采用不同类型的粒子表示，可避免重构界面几何形状的复杂过程，在模拟这类问题时具有明显优势。移动粒子半隐式(moving particle semi-implicit, MPS)方法^[1]和光滑流体动力学(smoothed particle hydrodynamics, SPH)方法^[2-3]是目前计算流体力学中常用的粒子法，采用完全Lagrange描述定义流体运动，将连续的流体介质离散为运动的粒子，并将流体控制方程中的微分算子转换为粒子间相互作用模型。

表面张力作为流体的重要性质，在小尺度的流动如气泡和液滴的破裂或聚合等过程中起着至关重要的作用。粒子法中常用的表面张力模型可分为两类：一类是粒子间相互作用力(inter-particle interaction force, IIF)模型，通过构建粒子之间的相互作用力来模拟粒子受到的表面张力，但模型中的经验参数通常难以确定，Nugent等^[4]最早在SPH中用这种方法模拟Van der Waals液体；另一类是连续表面力(continuum surface force, CSF)模型，它最早由Brackbill等^[5]提出并用于网格法，其原理是将作用在界面上的表面张力转换为界面周围流体的体积力加在Navier-Stokes方程中，可避免表面力计算中重构界面的复杂过程。Morris^[6]首次将CSF模型应用于SPH方法，其中曲率的准确计算是CSF模型的难点，Morris根据法线的长度对计算所得法线进行筛选，再进行曲率的计算。Adami等^[7]在多相流动SPH方法的CSF模型中采用了一种新的再生散度近似格式来计算曲率。Zhang^[8]和Duan等^[9]在粒子法中通过构造出界面或等值线的几何形状，计算流体界面的法线与曲率，对于复杂形状的界面其过程较为复杂。在早期MPS方法中，Nomura等^[10]采用一个简单的余弦模型来近似CSF模型中的曲率。

在传统的粒子法CSF模型中，界面附近某些内部粒子的法线计算错误，使界面曲率的计算产生较大误差。而表面张力是一种表面力，它仅取决于界面的形状。粒子法中，在准确识别界面粒子的基础上，表面张力可根据界面粒子计算得到而不依赖于内部粒子。据此，本文对CSF模型进行改进，首先采用几何法准确识别界面粒子，再通过计算界面粒子单位法线的面散度得到曲率，且表面张力仅作用于界面粒子。改进后的模型可避免原模型中内部粒子带来的误差，且格式简单易于计算。随后本文通过计算圆和椭圆的曲率验证改进后的曲率计算方法，并结合MPS法模拟方形液滴振荡和液滴碰撞过程的二维单相流动过程，以验证改进后的CSF模型在动态算例中的效果。

1 MPS方法

不可压缩流体的Lagrange型连续性方程和动量方程为

$ \frac{{D\rho }}{{Dt}} = 0, $

(1)

$ \rho \frac{{D\mathit{\boldsymbol{u}}}}{{Dt}} = - \nabla {p} + \mu {\nabla ^2}\mathit{\boldsymbol{u}} + \mathit{\boldsymbol{F}}. $

(2)

其中：ρ、p、μ、u分别表示流体的密度、压力、动力黏度和速度；F为质量力，包括重力和表面张力等。

粒子间的相互作用关系在MPS方法中通过权函数实现，

$ w\left( r \right) = \left\{ \begin{array}{l} \frac{{{r_{\rm{e}}}}}{r} - 1,\;\;\;\;0 < r < {r_{\rm{e}}};\\ 0,\;\;\;\;\;\;\;\;\;\;{r_{\rm{e}}} \le r. \end{array} \right. $

(3)

其中：r=| r_i-r_j|为2个粒子之间的距离，r_e为粒子作用域的半径。粒子数密度的计算式为

$ {n_i} = \sum\limits_{j \ne i} {w\left( {\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|} \right)} . $

(4)

其中粒子j是在粒子i作用域内的粒子。对于不可压缩流体，粒子数密度为常数并等于初始粒子数密度n₀，即n_i=n₀。

控制方程中的微分算子通过权函数转化为粒子之间的相互作用。任意标量ϕ的梯度和Laplace算子的离散格式分别为

$ {\left\langle {\nabla \phi } \right\rangle _i} = \frac{d}{{{n_0}}}\sum\limits_{j \ne i} {\left[ {\frac{{{\phi _j} - {\phi _i}}}{{{{\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|}^2}}}\left( {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right)w\left( {\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|} \right)} \right]} , $

(5)

$ {\left\langle {{\nabla ^2}\phi } \right\rangle _i} = \frac{{2d}}{{\lambda {n_0}}}\sum\limits_{j \ne i} {\left( {{\phi _j} - {\phi _i}} \right)w\left( {\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|} \right)} . $

(6)

其中d为空间维度，

$ \lambda = \frac{{\sum\limits_{j \ne i} {{{\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|}^2}w\left( {\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|} \right)} }}{{\sum\limits_{j \ne i} {w\left( {\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|} \right)} }}. $

(7)

MPS中每个时间步内的计算分2步进行：第1步显式计算动量方程中除压力以外的力对粒子的作用引起的Δu_i^*，得到粒子的临时速度u_i^*和位置r_i^*,

$ \mathit{\boldsymbol{u}}_i^ * = \mathit{\boldsymbol{u}}_i^n + \Delta \mathit{\boldsymbol{u}}_i^ * , $

(8)

$ \mathit{\boldsymbol{r}}_i^ * = \mathit{\boldsymbol{r}}_i^n + \mathit{\boldsymbol{u}}_i^ * \Delta t. $

(9)

第2步中为满足不可压缩条件，求解离散后的压力泊松方程^[11]

$ {\nabla ^2}{p_i} = \left( {1 - \gamma } \right)\frac{\rho }{{\Delta t}}\nabla \cdot \mathit{\boldsymbol{u}}_i^ * + \gamma \frac{\rho }{{\Delta {t^2}}}\frac{{{n^0} - n_i^n}}{{{n^0}}}. $

(10)

其中γ为常数，本文取γ=1×10^-3。求解全体粒子的压力泊松方程组采用不完全Cholesky分解共轭梯度法，之后利用更新后的压力场求解速度修正值：

$ \Delta {{\mathit{\boldsymbol{u'}}}_i} = - \frac{{\Delta t}}{\rho }\nabla {p_i}. $

(11)

其中计算压力梯度时使用的格式为^[12]

$ \nabla {p_i} = \frac{D}{{{n^0}}}\sum\limits_{j \ne i} {\left[ {\frac{{{p_j} + {p_i}}}{{{{\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|}^2}}}\left( {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right)w\left( {\left| {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{r}}_i}} \right|} \right)} \right]} . $

(12)

相比于式(5)，该格式既能保证粒子间压力为推力，又满足动量守恒。

最后再次修正粒子的速度和位置,

$ \mathit{\boldsymbol{u}}_i^{n + 1} = \mathit{\boldsymbol{u}}_i^ * + \Delta {{\mathit{\boldsymbol{u'}}}_i}, $

(13)

$ \mathit{\boldsymbol{r}}_i^{n + 1} = \mathit{\boldsymbol{r}}_i^ * + \Delta {{\mathit{\boldsymbol{u'}}}_i}\Delta t. $

(14)

2 表面张力

在经典CSF模型中，不同计算区域的粒子被赋于不同的色函数值，色函数的梯度方向即界面法线方向。微分几何中曲率可表示为

$ \kappa = - {\nabla _{\rm{S}}} \cdot \mathit{\boldsymbol{\hat n}}. $

(15)

其中$ \mathit{\boldsymbol{\hat n}}$为单位法矢量。在正交曲面坐标系(t₁, t₂, n)下，由于$ \mathit{\boldsymbol{\hat n}}$总是单位矢量，则有

$ {\nabla _{\rm{N}}} \cdot \mathit{\boldsymbol{\hat n}} = \mathit{\boldsymbol{\hat n}} \cdot \frac{{\partial \mathit{\boldsymbol{\hat n}}}}{{\partial n}} = \mathit{\boldsymbol{\hat n}} \cdot {\bf{0}} = 0. $

(16)

结合式(15)和(16)有

$ \kappa = - \left( {\nabla - {\nabla _{\rm{N}}}} \right) \cdot \mathit{\boldsymbol{\hat n}} = - \nabla \cdot \mathit{\boldsymbol{\hat n}}. $

(17)

理论上式(15)和(17)均可用于计算曲率。网格法的经典CSF模型中，在界面附近一定厚度的过渡区域内的流体元处采用式(17)计算曲率，可将表面张力转化为体积力而无需构造出曲面的形状。然而针对离散后的流体，式(16)仅在结构化、规则或正交的网格或粒子分布时成立。在粒子法中，虽然不可压缩流体在宏观上分布均匀，但局部粒子在一定程度上呈现不规则的随机分布；此时式(16)并不成立，故采用式(17)计算曲率则会引入误差或造成计算不稳定。而粒子法中界面粒子的识别较易实现，故计算表面张力时可不必转化为体积力的形式，而是在界面粒子处采用式(15)计算曲率，以避免∇_N· $\mathit{\boldsymbol{\hat n}}$≠0带来的误差。本文中关于表面张力的计算过程为：1)精准识别出界面粒子；2)对于界面粒子，计算单位法线及其面散度得到曲率和表面张力；3)在动态算例中求解压力泊松方程时将表面张力值作为压力边界条件。

2.1 界面粒子的识别

粒子法中识别界面粒子的基本原理为：根据目标粒子作用域被其他粒子充塞情况，界定内部粒子和界面粒子；如其作用域被其他粒子完全充塞，则界定为内部粒子，反之，界定为界面粒子。所采用的方法可划分为体积法和几何法两类。体积法^[1]是通过一个与体积有关的量来判断目标粒子作用域被其他粒子充塞情况，计算简单但准确度不高。几何法^[13-14]是通过粒子间的几何关系判断目标粒子作用域被其他粒子充塞情况，计算复杂但准确度较高。经典MPS方法中常采用反映目标粒子作用域充塞情况的体积量——粒子数密度界定界面粒子和内部粒子^[1]，方法简单但误判率很高。Dilts^[13]提出了用于SPH的3种几何法，其原理为将检查目标粒子作用域充填，转化为检查其作用域边界被其他粒子作用域的覆盖情况来识别界面粒子，其中弧度法准确度最高，但在某些粒子分布稀疏之处仍可能出现少量误判粒子。

本文采用一种直观的几何法来识别界面粒子，如图 1所示，以目标粒子为圆心作扇形，使扇形在粒子作用域内扫描一周，若扫描至某一位置处扇形内部没有充填任何一个粒子，则认为目标粒子为界面粒子，如图 1中A粒子。在动态算例中，由于粒子分布不均匀，可能出现某处内部粒子在某个方向上分布稀疏的情况，如图 1中的B粒子。故应设置合适的扫描角度和半径，在尽可能识别出小角度凹面上粒子的同时避免内部粒子被误判为界面粒子。本文选取扇形的圆心角为90°，半径为2.5l₀(l₀为初始粒子间距)。

2.2 表面张力的计算

假设表面张力系数为常量，故表面梯度为零，表面张力的方向为表面法线方向。考虑表面张力仅存在于表面粒子处，CSF模型可表示为

$ {\mathit{\boldsymbol{F}}_{{\rm{st}}}} = \sigma \kappa \mathit{\boldsymbol{\hat n}}{\delta _{\rm{s}}}. $

(18)

其中：F_st为表面张力矢量，σ为表面张力系数。表面张力仅作用于界面粒子，即对界面粒子有δ_s=1，内部粒子δ_s=0。引入色函数：

$ {c_i} = \sum\limits_j {\frac{m}{\rho }{W_{5{\rm{p}},ij}}} . $

(19)

单位内法线为

$ {{\mathit{\boldsymbol{\hat n}}}_i} = \frac{{\nabla {c_i}}}{{\left| {\nabla {c_i}} \right|}}. $

(20)

色函数的梯度计算采用SPH方法中一阶离散的对称格式^[6]

$ \nabla {c_i} = \sum\limits_j {m\left( {\frac{{{c_i} + {c_j}}}{\rho }} \right)\nabla {W_{5{\rm{p}},ij}}} . $

(21)

曲率的计算格式为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\kappa }}_i} = - \left( {{\nabla _{\rm{s}}} \cdot {{\mathit{\boldsymbol{\hat n}}}_i}} \right) = }\\ { - \sum\limits_{j \in B} {\frac{m}{{\rho _B}}\left( {{{\mathit{\boldsymbol{\hat n}}}_j} - {{\mathit{\boldsymbol{\hat n}}}_i}} \right)} \cdot \nabla {W_{5{\rm{p}},ij}}.} \end{array} $

(22)

计算单位法线面散度时仅选择界面粒子进行加权求和，故

$ {\rho _B} = \sum\limits_{j \in B} {m{W_{5{\rm{p}},ij}}} . $

(23)

其中B为界面粒子集合。由于SPH常用的核函数具有良好的归一性，本节计算中的式(19)、(21)—(23)均采用了SPH的离散格式及5次样条核函数^[15]

$ {W_{5{\rm{p}},ij}} = {\alpha _{\rm{d}}}\left\{ \begin{array}{l} {\left( {3 - s} \right)^5} - 6{\left( {2 - s} \right)^5} + 15{\left( {1 - s} \right)^5},\;\;\;\;\;\;0 \le s < 1;\\ {\left( {3 - s} \right)^5} - 6{\left( {2 - s} \right)^5},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;1 \le s < 2;\\ {\left( {3 - s} \right)^5},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;1 \le s < 3;\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;s \ge 3. \end{array} \right. $

(24)

其中：s=r_ij/h，h为光滑长度，α_d在一维、二维和三维空间内分别为120/h、7478πh²和3359πh³。

2.3 自由表面压力边界条件

传统的MPS方法在求解压力泊松方程时将自由表面上的粒子的压力设置为零，即假定的大气压力，这在表面张力可以忽略的情况下是合理的。当表面张力不可忽略时，气液界面间压力满足Laplace方程。本文在求解压力泊松方程时将界面粒子压力赋值为表面张力的大小且其符号与曲率一致

$ {p_{i \in B}} = {p_0} + \sigma \mathit{\boldsymbol{\kappa }}. $

(25)

其中大气压力p₀=0。

3 算例 3.1 曲率计算的验证

本节以圆和椭圆为例验证曲率计算的准确性。曲率的理论值为

$ {\kappa _0} = \frac{{\left| {y''} \right|}}{{{{\left( {1 + {{y'}^2}} \right)}^{3/2}}}}. $

(26)

取圆方程x²+y²=1，即圆形界面的曲率为1。考察粒子分辨率和光滑长度对计算准确性的影响：当粒子间距取0.1和0.05时，圆形界面粒子数分别为N_b=63和N_b=126，图 2为2种分辨率下计算所得表面张力分布；在2种分辨率下且光滑长度分别取h=1.0l₀、h=1.25l₀和h=1.5l₀时，圆周上所有界面粒子的曲率计算结果与理论值的偏差如表 1所示。可见采用改进的CSF模型后，不同分辨率和光滑长度下曲率计算误差均很小，表面张力计算结果准确。

为考察高曲率和曲率随空间变化的情况，计算椭圆x²/4+y²=1在第一象限内界面粒子处的曲率。根据式(26)可知当0≤x≤2，y≥0时曲率理论值为0.25≤κ₀≤2。当粒子间距为0.1和0.05时，1/4椭圆上对应的界面粒子数为N_b=25和N_b=49；光滑长度取h=1.0l₀和h=1.5l₀。表面张力计算结果如图 3所示，椭圆曲率计算的相对误差如图 4所示，其中横坐标为椭圆的长半轴方向。在x较小处椭圆的曲率和曲率变化均较小，计算误差很小。在x较大处，曲率迅速增加处作用域内的高曲率界面粒子参与计算，使所得曲率值偏大；曲率最大值附近，作用域内的较低曲率界面粒子使得计算曲率值偏小。光滑长度较大时通常误差较大，这是由于计算曲率时更大的作用域内的粒子被考虑在内。当分辨率增加时误差均显著减小，这是由于作用域减小，曲率计算更集中于局部粒子。高分辨率下光滑长度对曲率的影响也有所减弱。动态算例中发现，当光滑长度较小时，尽管计算的精确度更高，但其光滑性较差。考虑精确度和光滑性两方面的影响，本文随后的方形液滴振荡和液滴碰撞算例中选取h=1.5l₀。

3.2 方形液滴振荡

初始速度为零的方形液滴受表面张力和压力作用，其振荡过程存在分析解。本节通过模拟该经典算例考察动态过程中表面张力计算的准确性。取乙醇液滴边长为0.01 m，密度为797.88 kg/m³，表面张力系数为0.023 61 N/m，动力黏度为1.49×10^-3 Pa·s。液滴振荡频率的理论值为^[16]

$ \omega = \frac{{s\left( {{s^2} - 1} \right)\sigma }}{{\rho R_0^3}}. $

(27)

其中：s为振荡模数，液滴为方形时取s=4；R₀为液滴的等效半径。根据以上参数可得液滴的理论振荡周期为T₀=6.3×10^-2 s。图 5为粒子数为1 600时不同时刻下的液滴形状，其中黑色表示内部流体粒子，灰色表示界面粒子；计算所得液滴振荡周期为T=6.9×10^-2 s。初始时刻方形的四角受表面张力作用向内运动，液滴在前半个周期内演变为旋转45 °后的方形；后半个周期开始时四角受表面张力作用向内运动并演变回原始形状。图 6为粒子数分别为400、1 600和6 400时液滴的动能随时间的变化，其中动能变化的周期对应液滴振荡的半周期。可见随着分辨率的增加，液滴运动表现出更明显的周期性。当粒子数较少时，液滴仅呈现一个较明显的振荡周期，且振荡周期偏长。粒子数增加时，模拟过程中表面张力和压力的计算更精确，能实现更多的振荡周期，且周期更接近理论值。粒子数为1 600和6 400时振荡周期相差很小，此时表面张力的计算受分辨率的影响很小。

3.3 液滴正面碰撞

为进一步验证改进后的CSF模型，本节对受惯性力、黏性力、表面张力和压力作用的二维液滴碰撞过程进行模拟。2个相距0.25 m、相对速度为0.2 m/s的液滴正面碰撞。取液滴直径为0.1 m，密度为1 000 kg/m³，表面张力系数为0.5 N/m，动力黏度为0.2 Pa·s。对应无量纲参数Re=100，We=8。初始粒子间距为2.5×10^-2 m，2个液滴的总粒子数为2 626。图 7为不同时刻下的液滴形状。2个液滴碰撞后结合，结合处的凹液面在表面张力作用下向外凸出，而液滴的外侧在初速度的作用下向内运动并形成凹面。随后表面张力使凹面向外运动演变为凸面，之后的过程中液滴在表面张力作用下其形状在横椭圆与竖椭圆之间变化，且振幅越来越小。针对真空下二维液滴碰撞过程的SPH模拟，文[17]在计算表面张力时使用的是IIF模型；文[8]则采用了表面张力计算的几何方法。本文模拟液滴碰撞结合后形状演变过程与文[8]和文[17]相似；此外尽管本文中液滴Re更大，但并未出现文[17]中粒子分散现象，液滴表面更光滑，表明本文的表面张力计算结果更为合理。

4 结论

本文对适用于粒子法的CSF模型进行了改进。采用了一种几何法通过扇形扫描邻近区域的粒子充塞情况识别界面粒子，曲率由作用域内界面粒子单位法线的面散度计算得到，并且仅在界面粒子处存在表面张力作用。动态算例中求解压力泊松方程时将表面张力值作为压力边界条件。

通过计算圆和椭圆曲率对改进后的CSF模型进行验证，发现曲率是常数的圆形在不同分辨率和光滑长度下计算结果均非常准确，而椭圆在分辨率很低时存在一定误差，其中曲率较大处计算值偏小，曲率迅速增加处计算值偏大。多数情况下光滑长度的增加会引起误差的增加；当分辨率增加时误差显著减小，且光滑长度对计算误差的影响也明显减弱。随后在MPS方法中采用改进后的CSF模型进行了2组动态算例的计算：1)表面张力为主要驱动力的二维方形液滴振荡过程的数值模拟，得到的形状和振荡周期与理论值接近；当分辨率增加时振荡的周期性更明显，振荡周期也更准确。2)二维单相液滴碰撞过程的模拟，不同时刻下液滴形状模拟结果与已有研究结果相似，而在液滴表面光滑程度方面表现更好，表明表面张力计算结果更为合理。在2组动态算例的计算过程中同样可见界面粒子的识别效果也很好，未出现内部粒子被误识别的现象。