移动粒子半隐式法核函数特征对压力求解稳定性的影响

张凯,孙中国,席光

(西安交通大学能源与动力工程学院,710049,西安)

移动粒子半隐式法(MPS)是一种无网格数值方法,相比于传统网格剖分方法,移动粒子半隐式法采用拉格朗日描述,在处理大变形问题时,只需追踪粒子携带的位置、压力和速度等信息,避免了传统网格方法容易产生网格畸变以及需要频繁网格重构的情况,因而在处理具有大变形复杂边界或移动边界的问题时具有明显优势[1]。近年来,利用MPS法模拟楔形体入水、多孔介质流动等自由表面流体流动现象的方法得到了广泛的研究与应用[2-4]。流体压力采用隐式求解,求解精度和稳定性将直接影响MPS算法的整体精度,压力求解的数值振荡容易导致模拟精度下降甚至不收敛,因此有效抑制MPS法在流动模拟中的压力振荡对提升算法精度与稳定性具有重要意义。

国内外学者已从多方面开展了抑制MPS算法压力振荡的研究,包括压力梯度修正、表面粒子判定、添加泊松方程源项等。Khayyer等通过在压力梯度以及泊松方程中添加源项并允许轻微的可压缩性,提出了改进MPS法用来提升预测沿海结构波浪冲击压力的精度[5]。Shibata等提出了在自由表面上布置虚拟粒子的方法来抑制压力振荡,该方法同时可模拟负压现象[6]。Debadatta等通过修改压力泊松方程源项,并采用混合表面粒子判定技术以及改进的固液边界条件,获得较稳定的压力求解,模拟了剧烈晃动现象[7]。

上述学者对压力振荡的改进主要集中在压力梯度以及泊松方程方面,对核函数的研究相对较少,实际上,核函数数值特征对模拟的稳定性具有重要影响。Ataie-Ashtiani等对不同核函数提高溃坝模拟稳定性做了相关研究,发现紧支域半径以及粒子大小对模拟稳定性具有较大影响[8]。潘徐杰等发现在粒子相撞时具有有限值的核函数可减缓压力振荡,并研究了5种常见核函数对模拟的适用性[9-10]。Shobeyri发现,与具有陡峭形状的核函数相比,扁平形状的核函数可提高模拟的稳定性[11]。朱跃等采用两种核函数模拟了溃坝问题,与传统核函数相比,采用高斯型核函数模拟得到的粒子数密度偏差波动较小,压力振荡幅度较小[12]。

本文分析了MPS法核函数形状特征对压力求解稳定性的影响,通过归纳总结主要特性,提出了一种可有效抑制压力振荡的新核函数。针对静水压力问题的静压求解以及液体晃动问题的动力学压力计算,采用典型算例进行验证。

1 数值方法

1.1 移动粒子半隐式法

移动粒子半隐式法中,流体在流动过程中遵循能量守恒方程和动量守恒方程,矢量形式为

(1)

(2)

式中:ρ为密度;t为时间;u为速度矢量;p为压力;μ为动力黏度;f为重力;为梯度算子;2为拉普拉斯算子。算子式分别为

(3)

(4)

(5)

式中:d为维度;〈〉i为粒子i物理量的近似值;xij为粒子i与j之间的x向距离;wij为核函数;r=|ri-rj|,ri、rj为粒子i、j的坐标;re为核函数的影响半径,一般情况取2.1l0,l0为单个流体粒子的直径。

粒子i的粒子数密度可表示为

(6)

对于不可压缩流体,粒子数密度保持不变,取初始布置时的粒子数密度n0。

粒子压力可通过求解压力泊松方程得到

(7)

为了使计算结果更稳定,调整泊松方程的源项[13],该方法可获得更稳定的压力场,即

(8)

式中γ为(0,1)之间的常值。

1.2 核函数基本要素与构造思路

在MPS算法中,粒子间的相互作用用核函数加权表示,传统核函数和加权形式参见式(5)(6)。当粒子间距r>re时,核函数值为0,忽略粒子间相互作用,核函数满足紧支性条件并可提高计算效率。计算中物理参数值可通过核函数加权计算获得,核函数满足非负性条件可防止出现非物理现象。为了描述粒子间相互作用的实际情况,即粒子间距离越小,粒子间的相互作用越大,核函数须满足单调递减条件。

传统核函数在粒子间距非常小时值为无穷大,导致在压力方程建立和求解过程中容易出现数值振荡现象。众多学者在模拟过程中会采用不同形状核函数以减轻模拟过程中的压力振荡现象,一种常见的核函数处理方式是在粒子相互碰撞时,将核函数的值改为有限值[9-11];另外,核函数的形状特征对压力求解稳定性的影响是本文关注的重点。基于上述分析,为便于把握曲线的形状特征,本文构造有限值核函数及以指数函数为基础的多项式形式核函数,重点研究核函数的极值和形状特征,并分别通过流体静力学和流体动力学中两个典型算例,即静水压力求解和液体晃动压力求解,来分析并验证指数多项式型核函数对抑制压力振荡现象的效果。

2 静水压力求解

构造两种有限值核函数,基于静水压力问题模拟,对比分析不同核函数在模拟过程中对压力求解稳定性的影响。分别从核函数是否为有限值、核函数表征形状特点以及核函数真实形状特征3个方面,分析归纳抑制压力振荡的核函数特性。本文算例为二维模型,重力加速度g=-9.8 m/s2。

2.1 算例形状及初始布置

算例模拟区域和粒子初始布置如图1所示,其中水池底面中点p1为压力监测点,l0=0.008 m,液体部分离散为2 700个粒子,壁面粒子布置为3层。p1点压力的解析解为

p=ρgh

(9)

式中h为水池深度。

(a)几何模型 (b)粒子初始离散模型图1 静水问题模型

2.2 讨论及结果分析

2.2.1 有限值核函数 分别采用3种不同核函数模拟静水压力问题,其中w1(r)为传统核函数,在粒子间距为0时核函数值为无穷大,w2(r)为有限值核函数[14],w3(r)为本文构造的指数型多项式核函数,为有限值且在整个支持域内无奇点,即

(10)

(11)

3种核函数的曲线如图2所示。3种不同核函数模拟静水问题时,监测点p1点的压力随时间变化的曲线如图3所示。由图3可知:传统核函数w1(r)压力振荡比较明显,最终计算不收敛;核函数w2(r)虽在粒子相互碰撞时采用了有限值,但未能明显抑制压力振荡;新构造核函数w3(r)压力振荡得到明显抑制,且模拟结果与解析解吻合较好。

图2 核函数w1(r)、w2(r)、w3(r)曲线

图3 3种核函数模拟静水问题时p1点压力随时间的变化

采用w1(r)时,当相邻粒子间距较小时,其核函数值将迅速增大为较大数值(趋近无穷大),造成相邻粒子压力巨大,发生粒子飞溅弹开的数值现象,这与实际物理过程不符。当采用有限值核函数w2(r)时,压力振荡仍然较严重,这是由于当粒子间距很小时,核函数值随粒子间距变化依然比较敏感(大斜率),这使得粒子间的作用关系在间距较小时呈现出不稳定的趋势,成为诱发压力振荡的重要因素。

本文提出变化率相对较小的w3(r),在粒子间距较小时,其核函数值增大至合理范围,与周边粒子及前一时刻自身的核函数值相差较小,不会产生过大压力,粒子运动相对稳定,避免了剧烈飞溅。结果分析表明,压力求解稳定的核函数须在满足粒子间距与核函数值成反比例关系的基础上,同时应具备在粒子间距较近,即r/re为0～0.1时为有限值且变化率较小的形状特征。

2.2.2 核函数表征形状特点 为便于探究核函数曲线的形状特点,尤其是[0,1]区域两端曲线斜率对压力振荡的影响,本文构造核函数采用w(r)=exp(-a(r/re)2)-exp(-a)的形式,其中a为常数,当r/re=1时核函数值为0,保证了核函数的连续性。a取不同的值时,核函数曲线具有不同的形状特点,a越小核函数形状越呈现出扁平状,核函数值在0.9

采用核函数w4(r)(a=3)、w5(r)(a=6)和w6(r)分别进行数值模拟,如图4所示,监测点p1历时10 s,压力随时间变化的曲线如图5所示。核函数w4(r)中a=3时,压力振荡较明显;核函数w5(r)中a=6时,压力与解析解一致,随时间分布均匀且较稳定。

图4 核函数w4(r)、w5(r)、w6(r)曲线

图5 p1点压力随时间的变化

由图4、图5对比可知,核函数值在r/re=1附近变化率较大时,容易出现压力求解不稳定。若粒子间距较大时核函数值较大,粒子相互接近的趋势将会受到抑制,因此在r/re=1附近,核函数值应始终维持在0附近,核函数变化率不应过大。压力求解稳定的核函数具有在r/re=0处为有限值,且在r/re=0,1两端处应尽量保持平缓变化的形状特点。

2.2.3 核函数真实形状特征 构造核函数采用wm(r)=bexp(-6(r/re)2)-bexp(-6)的形式,其中b取不同值时,核函数wm(r)在粒子间距为0时有限值的大小不同。当b分别取1、5、9时,核函数w7(r)、w8(r)、w9(r)曲线如图6所示,核函数与其对应粒子数密度常数的比值wn(r)=wm(r)/n0。可知p1点压力与b无关,b取不同数值时,p1点压力完全一致。

图6 核函数w7(r)～w12(r)曲线

SPH算法中光滑函数的构造需要满足归一性条件[15],而MPS算法并不需要强制满足,物理量参数的近似值可表示为

(12)

MPS算法中粒子数密度作为规范化因子,对粒子各作用模型进行规范化处理。

当b取3种不同数值时,wn(r)所表达的3条曲线完全重合,模拟得到包含p1点在内的压力值分布完全一致。核函数与粒子数密度常数的比值在压力求解时起重要作用,尽管核函数曲线形状和极值不同,但wn(r)的曲线能真实反映粒子间的相互作用。

3 液体晃动求解

3.1 初始布置

液体晃动动压问题模拟区域和初始布置如图7所示,其中液箱在水平方向上做正弦运动,液箱运动规律为vx=0.02×5.39sin(5.39t),p2点为液箱右侧壁一点,p2点与液箱底部距离为55 mm,粒子直径l0=0.005 m,布置8 000个液体粒子,壁面粒子布置为3层。

3.2 讨论及结果分析

通过静水压力模拟发现,核函数在[0,1]区间曲线两端附近数值变化较平缓时,可有效改善压力求解振荡现象;当核函数在r/re=0.8附近数值过小时,系统动力学计算会出现偏差。为了兼顾压力求解稳定性与保证系统的动力学性能,对核函数做进一步改进,可提高r/re=0.8附近核函数值,新的分段指数多项式型核函数为

(13)

核函数改进前后典型时刻液箱内自由液面模拟结果与实验的对比如图8所示。改进前流体压力值计算存在偏差且有瞬态波动,如图中A、B、C区域所示,同时流体内部出现较多压力零散点;改进后压力振荡明显减弱,压力分层分布更加合理,压力零散点基本消失,如图中A′、B′、C′区域所示。与w1(r)相比,w6(r)在压力求解过程中具有更好的稳定性,液体在晃动过程中自由液面的变化规律与实验结果更加吻合。

(a)几何模型

(b)粒子初始离散模型图7 液体晃动模型

改进前后液体晃动算例监测点p2的压力计算结果与实验结果[16]对比如图9所示。采用改进前传统核函数,模拟结果存在较大振荡,造成较大误差,双峰现象较模糊;改进后压力振荡被抑制且成功捕获双峰现象,与实验结果吻合较好。压力峰值计算较为准确,谷值较实验值略大,需进一步分析液体晃动中数值耗散产生的影响。指数多项式型核函数在整个区间内连续且最大值为有限值,函数构造简单易于实现,可有效抑制压力振荡现象,与传统核函数相比,压力求解具有更好的稳定性,可广泛适用于普通流体的模拟。

t=nT+0.2T t=nT+0.5T t=nT+0.7T(a)实验结果[16]

t=6.06 s t=6.41 s t=6.64 s(b)w1(r)模拟得到的压力分布云图

t=6.06 s t=6.41 s t=6.64 s(c)w6(r)模拟得到的压力分布云图图8 不同时刻液箱内自由液面实验结果与模拟结果对比

图9 改进前后p2点压力变化的计算结果与实验结果[16]比较

4 结 论

本文分析并研究了核函数曲线形状特点对MPS法流体压力求解稳定性的影响,构造了指数多项式型核函数,通过模拟静水压力及液体晃动典型问题,验证了所提核函数可准确模拟液体的压力分布及监测点的压力变化,并可有效抑制求解过程中的压力振荡现象。研究发现,核函数的最大值为有限值且r/re在[0,1]区间曲线两端附近数值变化平缓时,可有效改善压力求解振荡现象,在r/re=0.8附近核函数值过小会影响系统的动力学性能。此外,核函数的极值及曲线形状仅是表征特点,核函数与粒子数密度的比值可真实表征粒子间的相互作用,在压力分析稳定性中起重要作用。