精确可靠的模拟表面张力对液滴变形与运动的影响在环境工程、微纳米工程以及医药工程等领域有着十分重要的应用. 采用基于网格的数值方法模拟此类问题往往需要复杂的界面追踪算法, 典型的如流体体积函数法^[1,2]和水平集方法^[3]等. 因此, 越来越多的研究者采用拉格朗日粒子类方法模拟表面张力相关现象, 利用粒子本身实现对边界的自动追踪, 极大地简化了算法, 其中最常用的数值方法是光滑粒子动力学(smoothed particle hydrodynamics, SPH)^[4,5]方法. 在SPH方法中, 主要有两种方法来考虑表面张力的影响: 粒子间相互作用力(inter-particle interaction force, IIF)方法^[6-10]和连续力方法(continuum surface force, CSF) ^[11]. IIF方法来源于如下思想, 即表面张力来自于液体表面分子之间和内部分子之间作用力的差异. 受此启发, 在SPH方法中, 可以通过在宏观粒子之间加入额外远程吸引力和近程排斥力来近似模拟表面张力. 这种方法简单灵活, 粒子之间的额外作用力形式非常多. 这种方法的缺点是表面张力系数不是输入参数, 而是需要先模拟某个标准算例, 然后测量与这种宏观粒子之间吸引力等价的表面张力系数. 此外, 这种方法往往存在不收敛的问题. CSF方法通过给不同的流体赋以不同的色值(color)来区分不同的流体, 利用不连续的色函数(color function)求解表面法向以及曲率来计表面张力, 然后将只作用在表面的力转换为连续的体积力. SPH方法中, 界面法向向量及其曲率的计算误差一般较大. 无论哪种表面张力模拟方法, SPH方法都存在如下固有缺点: 1) SPH形函数不能严格满足1阶精度(特别是粒子分布不均匀或者自由界面处), 导致收敛性不好, 常常需要对形函数及其导数进行修正^[12,13]; 2)形函数不满足Kronecker delta性质, 因此位移边界条件处理复杂, 一般需要通过各种虚粒子处理边界条件. 因此, 本文尝试采用一种新型无网格方法模拟表面张力问题, 即最优输运无网格(optimized transportation meshfree, OTM)方法^[14].

OTM方法采用粒子系统离散连续介质, 利用局部最大熵(local maximum entropy, LME)^[15-17]形函数构造连续场, 并结合数学中的最优输运理论计算系统的动能积分. 相比于其他无网格方法, OTM方法在精度、收敛性、边界处理以及数学理论基础上均具有优势. OTM方法具有无网格方法能够处理大变形的优点, 同时又继承了有限元方法的部分优点, 如形函数严格满足一阶精度、形函数在边界处满足弱Kronecker delta性质等. 当前, OTM方法已经成功应用于金属侵彻^[18]、超高速碰撞^[19-21]、裂纹扩展^[22]以及热传导等问题^[23], 但是将其应用于表面张力模拟还未见文献报道.

本文从拉格朗日方程出发, 从能量的角度考虑表面张力的影响, 即将表面能加入到拉格朗日函数中, 得到的表面张力广义力精确地作用在流体表面, 而且表面张力系数是唯一的输入参数. 为了减少计算量, 推导了轴对称形式的OTM离散格式, 用于三维算例数值模拟, 并给出了详细的计算流程. 通过模拟泊肃叶流并与解析解对比, 验证了本文基本离散格式和轴对称处理的正确性. 通过模拟静止液滴, 并与Young-Laplace公式对比, 验证了表面张力处理的正确性. 模拟了液滴的周期振荡问题并与解析解对比, 严格验证了液滴中速度分布的空间收敛性. 最后, 模拟了液滴在重力、表面张力以及界面共同作用下的形变问题, 考察了液滴中的压力分布并与理论解对比, 进一步验证了本文OTM离散格式.

OTM方法中的连续介质离散介于拉格朗日有限元与SPH方法之间. OTM方法将有限元网格中每个小网格转换为1个粒子, 位于小网格形心. 对于三角形, 形心为三条中线的交点. 粒子存储所有的物理量, 与SPH方法类似. 同时, 保留所有的网格节点, 其中存储了速度和位置, 如图1所示. 初始时刻, 边界上只有节点, 且节点随着流场运动. 边界节点很好地表征了连续介质的边界位置, 与拉格朗日有限元方法类似. 此外, 由于局部最大熵形函数满足弱Kronecker delta性质, 边界条件处理简单. 抛弃有限元方法中单元和节点之间的固定连接关系, 而是采用近邻粒子搜索方式动态确定, 从而可以处理大变形, 这与SPH方法类似.

最优输运无网格方法与更新拉格朗日有限元方法最大的1个区别是采用LME形函数, 因此不需要网格, 属于无网格方法. 给定1个节点集$ {\boldsymbol{x}}_{a}(a= $$ 1, 2, \cdots , n) $, 每个节点的形函数记为$ {N_a}({\boldsymbol{x}}) $. 利用形函数可以从离散数值插值出1个连续场:

其中, $ {u^h}({\boldsymbol{x}}) $是插值出来的连续场, $ {u_a} $是物理量u在节点a处的值. 在有限元和SPH等方法中, 形函数及其导数都具有简单解析表达式, 但是LME形函数没有显式的解析表达式, 而是需要数值求解.

如果要求(1)式可以精确重构常数场, 那么有$\displaystyle \sum\nolimits_a {{N_a}({\boldsymbol{x}}) = 1}$, 即形函数满足归一化条件. 此外, 对于固定的x, 如果形函数大于零, 那么形函数$ {N}_{a}(a=1, 2, \cdots , n) $可以看成是某个事件的概率分布, 其熵为

其中, 由于x是固定值所以将其省略. 为了减少计算量, 形函数一般都是局域的, 即当$ \left| {{\boldsymbol{x}} - {{\boldsymbol{x}}_a}} \right| $大于某一临界值时, 形函数$ {N_a}({\boldsymbol{x}}) $为零. 因此, 可以定义如下形函数影响域平均大小:

形函数影响域越大, 计算量越大. LME形函数的基本思路是让形函数熵最大, 同时影响域又尽量小, 因此, 构造了如下优化问题:

其中, β是1个可调参数, 用于调节两个目标的权重. 最后1个约束条件表示形函数严格满足一阶精度. 上述优化问题解的存在性、唯一性以及求解方法文献[15-17]中有详细论述, 本文只给出最后结果. LME形函数可以写成如下形式:

其中Z函数的定义为

$ {{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}) $是以下无约束优化问题的解:

定义:

那么, 可以采用如下迭代方法求解$ {{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}) $:

初始时刻, $ {{\boldsymbol{\lambda }}^{(0)}} $一般取零. 流场发生变化后, $ {{\boldsymbol{\lambda }}^{(0)}} $一般取上1个时刻的$ {{\boldsymbol{\lambda }}^*} $. 当β等于零时, 得到的最大熵形函数是全局的, 不适合数值计算. 当β趋于无限大时, 即不考虑熵最大, 只要求形函数影响域最小, 那么LME形函数可以连续过渡到有限元线性形函数. LME形函数具有以下特点: 1)满足一阶精度; 2)不依赖网格; 3)边界处满足弱Kronecker delta性质; 4)形函数非负.

得到$ {{\boldsymbol{\lambda }}^*}({\boldsymbol{x}}) $与形函数后, 形函数导数为

首先考虑无黏、可压缩以及等温流体. 当初始状态已知, 节点的位置描述了流场的形变的流动, 实际上完全描述了流场, 因此, 节点的位置可以当成系统的广义坐标. 本文下标p表示粒子, 用下标a, b, c表示节点. 由于无黏, 这是1个保守系统, 拉格朗日方程为

其中, E是系统动能, V是系统势能, $ {{\boldsymbol{f}}_a} $是广义节点力, $ {{\boldsymbol{v}}_a} $是节点速度, $ {{\boldsymbol{x}}_a} $是节点位置, N是节点总数. 系统动能和势能为所有粒子动能和内能之和:

其中$ {m_p} $表示粒子质量, $ {{\boldsymbol{v}}_p} $表示粒子速度, $ {e_p} $表示粒子单位质量内能, $ {\rho _p} $表示粒子密度.

粒子内能和压力的关系通过物态方程(equation of state, EOS)描述:

与有限元一样, 粒子位置和速度可以通过节点位置和速度插值得到:

其中$ {\text{δ}} {{\boldsymbol{x}}_{p}} $和$ {\text{δ}} {{\boldsymbol{x}}_a} $分别为粒子和节点的虚位移, $ {N_a}({{\boldsymbol{x}}_{p}}) $表示节点a处形函数在粒子p处的值, $ \nabla {N_a}({{\boldsymbol{x}}_p}) $表示形函数梯度, 求和是对粒子的所有近邻节点. 粒子密度变化和位置变化的关系由连续性条件给出:

将(14e)式代入(15)式, 得到

利用求导链式规则, 并且利用(13)式和(16)式, 得到

利用(12a)式、(14c)式和(14d)式得到

在(18)式中交换求和顺序, 拉格朗日方程变成:

其中, M为相容质量矩阵

(19)式与更新拉格朗日有限元格式有两点区别: 1)有限元格式中的单元积分对应了OTM中的粒子积分; 2)有限元中使用基于网格的形函数, 而OTM方法中使用无网格的局部最大熵形函数. 注意到局部最大熵形函数可以连续过渡到有限元线性形函数, 因此, OTM方法可以看成是一种特殊的单点积分拉格朗日有限元方法.

在有限元中方法, 为避免矩阵求逆, 常常将相容质量矩阵按行求和, 得到集中质量矩阵, 大大减少了计算量和存储量, OTM方法中同样如此. 利用形函数的单位分解性质, 可以将粒子的质量分配给其相邻的节点, 得到如下节点质量定义:

注意到粒子质量和局部最大熵形函数都是严格大于零的, 因此, 节点质量也是严格大于零的. 这样连续介质既可以看成是粒子系统, 也可以看成是节点系统. 很明显, 这两种系统的总质量和总动量都是严格相等的.

系统的总势能通过粒子系统计算. 系统的总动能既可以通过粒子系统计算, 也可以通过节点系统近似计算:

将(23)式代入(11)式, 得到如下离散方程:

从(20)式和(21)式可知, 节点质量刚好就是相容质量矩阵按行求和. 利用节点广义力和节点质量可以更新节点的速度和位置:

用(25b)式更新节点位置后, 粒子位置根据(14a)式更新. 注意到节点位置变化代表了流场的形变, 因此两个相邻时刻的连续映射$\varphi :{{\boldsymbol{x}}^k} \to {{\boldsymbol{x}}^{k + 1}}$可以通过更新的节点位置插值得到

上述映射的梯度即为变形梯度张量:

变形梯度张量的行列式在粒子位置处的值表征了粒子的密度变化

(28)式可以用于更新粒子密度. 粒子密度更新后, 根据物态方程更新粒子压力.

以上推导没有考虑黏性, 因而是保守系统. 处理非保守力的基本想法将其看成外力, 然后利用虚功原理推导对应的广义力. 假设$ {{\boldsymbol{f}}_p} $是作用在粒子p上的外力, 那么此力的虚功为

由虚位移的独立性, 得到相应的广义力为

例如, 作用在1个粒子上的重力$ {{\boldsymbol{f}}_p} = {m_p}{\boldsymbol{g}} $. 根据(30)式得到广义力${{\boldsymbol{f}}_a} = \displaystyle \sum\nolimits_p {{m_p}{\boldsymbol{g}}{N_{a, p}}}$. 考虑到重力实际上是保守力, 因此也可以将重力势能加入拉格朗日函数中, 得到重力对应的广义力. 简单计算可以证明两者是一致的.

材料的很多行为(如弹塑性, 黏弹性等)都可以统一用1个应力张量描述, 而应力张量的影响也可以看成外力, 如下式中的动量守恒方程:

应力张量散度的物理意义为单位体积受到的力, 因此, 应力张量虚功为

其中用到了高斯积分公式, 并且假设了边界上虚位移或者应力张量为零, 这在固壁边界以及自由表面边界上都成立. 同理, 得到广义力为

如果只考虑牛顿流体的黏性, 那么黏性应力张量为

其中μ为黏性系数.

同理, (31)式也能通过势能和拉格朗日方程得到. 从(33)式可以看到, 1个粒子在无限短时间内受到的力只和应力张量本身有关, 而和应力应变关系无关. 因此, 可以临时假设1个线性应力应变关系, 从而得到简单的弹性势能. 对弹性势能求变分, 就得到了相应的广义力, 详细证明见附录A. 两种方法的结果是一致的.

考虑二维情况, 液体边界被边界节点表征, 如图2所示.

表面势能正比于表面积:

其中$ {r_{ab}} = \left| {{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_b}} \right| $, $ {r_{ac}} = \left| {{{\boldsymbol{x}}_a} - {{\boldsymbol{x}}_c}} \right| $为节点之间的间距. γ是表面张力系数. 对(35)式求变分, 得到广义力

可以看到, (36)式非常简洁, 可以解释为节点a同时受到两个相邻节点b和c的吸引力, 吸引力大小刚好等于表面张力系数. (36)式中没有可调参数, 表面张力系数是唯一的输入参数. 此外, 表面张力对应的广义力只作用在表面节点上, 即精确的作用在流体表面.

记(r, θ, h)为柱坐标系, 并且假设环向没有位移, 而且所有物理量都和环向坐标无关, 即$ {\text{δ}} {u_\theta } = 0 $, $ \partial /\partial \theta = 0 $. 定义:

那么, (11)—(14)式依然成立. 速度矢量散度在柱坐标下的表达式为

因此, 粒子密度变化和粒子位置变化关系式(15)需要修改为

于是

其中$ {{\boldsymbol{n}}_r} $是径向单位向量. 于是, 轴对称条件下的广义力为

(18)—(26)式仍然成立. 注意(27)式中的变形梯度张量的行列式表征的是粒子在(r, h)平面上的面积的变化:

因此, 可以先更新粒子在(r, h)平面上的面积$ {S^{k + 1}} $, 然后再更新密度:

很明显, 外力对应的广义力(30)式仍然成立. 黏性力对应的广义力也仍然成立, 详细证明见附录B.

三维轴对称条件下, 表面势能为

其中, $ {\tilde r_{ab}} = {r_a} + {r_b} $, $ {\tilde r_{ac}} = {r_a} + {r_c} $. $ {r_k}(k = a, b, c) $为节点的径向坐标. 对(44)式求变分, 得到表面张力对应的广义力:

式中, 前两项为节点b和节点c对节点a的吸引力, 后两项为对称轴对边界节点的吸引力. 可见, 笛卡尔坐标系下的OTM计算格式只需要经过微小改动, 就能计算三维轴对称问题.

笛卡尔坐标系下, 采用OTM方法模拟表面张力的详细步骤总结如下.

1)初始化: 给定节点的位置和速度, 给定粒子的位置和其他物理量.

2)近邻粒子搜索: 找到所有粒子-节点对.

3)对于所有粒子-节点对, 计算形函数$ {N_{a, p}} $和形函数导数$ \nabla {N_{a, p}} $.

4)计算节点质量、粒子速度梯度、黏性应力张量和广义节点力:

5)更新节点速度和位置:

6)根据边界条件调整节点位置, 然后更新粒子位置:

7)计算变形梯度张量并更新粒子密度和压力:

8)完成1个时间步长, 回到步骤2).

所有的初始化借助于有限元网格. 网格节点变成节点, 网格形心放置粒子, 粒子体积即为网格大小. 物态方程如下:

其中, $ {\rho _0} $表示材料初始密度, ρ表示当前密度, $ {c_0} $表示声速, p表示压力.

两块无限长的二维平板中间充满了静止液体, 液体在平行于平板方向的体积力作用下开始运动, 称之为泊肃叶(Poiseuille)流. 如果是无限长的三维圆管, 那么称之为Hagen-Poiseuille流. 对于二维泊肃叶, 参数和文献[24]一样, 解析解见文献[24]. 将平板距离变成圆管半径, 即为Hagen-Poiseuille流, 解析解见文献[25]. OTM模拟结果见图3.

OTM模拟结果与理论解吻合很好. 为了模拟无限长管道, 用到了周期边界条件.

二维液滴的半径R = 0.2 m, 密度ρ = 1 kg/m³, 表面张力系数γ = 1 Pa·m, 黏性系数η = 0.5 Pa·s, 声速c₀ = 50 m/s. 静止状态下液滴的理论压力值可以根据Young-Laplace 公式得到ρ = γ/R = 5 Pa. 如果是三维液滴, 那么其理论压力值为ρ = 2γ/R = 10 Pa. OTM模拟的平均压力分别为 5.004 Pa和 9.967 Pa. 压力分布如图4所示.

由于弱可压假设, 压力场中有误差. 但是液滴中的平均压力非常接近理论值(误差0.5%以内), 证明了本方法的精度.

在上1个算例的基础上, 将黏性减小到$ \eta = $$ 5 \times {10^{ - 3}} $ Pa·s, 并且附加1个散度为0的初始速度场$ {v_x} = x $, $ {v_y} = - y $给所有节点. 二维液滴的理论振荡周期为$ T = 2{\text{π }}\sqrt {\rho {R^3}/6\gamma } $ ^[26]. 为了验证OTM方法的收敛性, 采用不同的粒子总数模拟这个问题, 并且检查两条相互垂直线$ x = 0 $和$ y = 0 $上的物理量分布, 时刻$ t = T/2 $, 如图5所示.

可以看到, 压力分布的收敛不明显, 但是速度分布的收敛很明显. 改变表面张力系数的大小, 得到的周期与理论值的对比如图6所示. 周期是根据液滴右上部分质心轨迹测量得到, 与文献[26]一样.

对于轴对称的三维液滴, 散度为零的初始速度场取为$ {v_{\text{r}}} = r $, $ {v_h} = - 2 h $, 其他物理量不变. 振荡周期的理论解为$ T = 2\pi \sqrt {\rho {R^3}/8\gamma } $^[27]. 模拟结果与理论的对比如图7所示. 模拟的周期根据对称轴端点的位置轨迹测量. 可以看到, 两种情况下OTM计算结果和理论非常吻合.

将液滴半径扩大到 1 m, 表面张力系数扩大到10 Pa·m, 加上大小为10 m/s²的重力. y = –1 m处放置1个光滑的固定平板. 液滴在重力、表面张力以及光滑平板的共同作用下产生变形, 如图8所示.

在图8(d)中, 压力等值线与重力垂直, 和理论相符. 液滴内部的压力场需要满足两个条件, 一是压力沿Y轴的线性分布规律$ {\text{δ}} p = \rho g{\text{δ}} Y $, 另一个是自由表面附近的压力满足Young-Laplace方程$ p = \gamma /R $. 这两个条件实际上决定了液滴的形状, 即液滴自由表面部分的微分方程为

其中H和p₀分别为液滴最高点的Y坐标和压力. 上述理论曲率半径在实际计算时误差太大, 因此, 本文采用如下方法计算曲率半径: 对于每个边界节点, 找到与其相邻的4个表面节点, 然后用1个圆来拟合这5个节点, 即$\mathop {\min }\limits_{{x_0}, r} \displaystyle \sum\nolimits_{i = 1}^5 {{{(\left| {{x_i} - {x_0}} \right| - r)}^2}}$, 自由变量为圆的中心坐标和半径. 圆的半径即为此节点的曲率半径.

液体中自由边界的右边部分如图9(a)所示. 采用曲率半径计算的边界压力和根据EOS计算的压力比较如图9(b)所示.

两种方法得到压力值大致相等, 说明模拟的自由表面形状和理论相符. 靠近固壁处误差较大, 经检查是因为固壁附近节点的分布更加无序, 导致搜索算法对初始参数敏感. 压力在Y轴方向的线性分布和理论相符, 拟合直线的斜率也和理论$ {\text{d}}p/{\text{d}}Y = - \rho g = - 10 $ Pa/m相符.

本文基于拉格朗日方程, 通过将表面张力势能加入拉格朗日函数, 得到了OTM方法框架下的表面张力效应表达式, 并且给出了轴对称的处理方式. 模拟了泊肃叶流, 静止和振荡的液滴, 液滴在重力、表面张力以及光滑平板共同作用下的形变等典型算例, 通过与解析解对比, 验证了OTM方法在模拟表面张力现象中具有如下优势: 1)表面边界节点精确的表征了边界的位置, 保证了表面张力势能的准确计算, 而且表面张力对应的广义力精确的作用在流体表面; 2)表面张力对应的广义力形式简洁, 具有直观的物理意义, 而且表面张力系数是唯一的输入参数, 不需要任何可调参数; 3)可以保证速度分布的空间收敛性.

本文在连续介质离散时, 节点和节点、节点和粒子之间没有任何固定联系, 而是通过近邻粒子搜索动态确定. 但是, 在计算表面张力势能时, 为了方便计算表面积, 假设了表面节点之间的连接关系不变. 因此, 所有表面节点可以看成是1个表面有限元网格. 表面有限元网格与内部节点在拉格朗日方程框架下相互作用和协调运动. 下一步, 将考虑表面拓扑结构发生变化的情况, 如多个液滴的融合过程.

下面给出(33)式的势能推导方法. 将应力和应变张量写成向量形式${:} $

应变与位移的关系为

其中L是微分矩阵

u是位移. 广义胡克定律为

其中c是1个对称矩阵:

以下4个等式成立:

将(16)式写成矩阵形式:

其中约定了形函数导数$ \nabla {N_{a, p}} $为列向量. 对于固定的应力张量$ {{\boldsymbol{\sigma}}_p} $, 由于c可以取得足够大, 使得应变为小变形, 于是弹性势能为

对弹性势能求变分, 得到

将向量形式改写为常用的张量形式, 得到

(A12)式第二项中有应变, 因此是小量, 可以忽略

(A13)式和 (33)式一致.

两个张量的双点乘是1个标量, 因此, 可以在柱坐标下求解

考虑到轴对称, 得到

利用定义(37a)式和(37c)式, 虚功为

(B3)式与(32)式相同.