光滑粒子流体动力学方法固壁处理的一种新型排斥力模型*

2013-12-12韩亚伟强洪夫赵玖玲高巍然

物理学报 2013年4期

韩亚伟强洪夫赵玖玲高巍然

(第二炮兵工程大学601室,西安 710025)

(2012年8月10日收到;2012年9月16日收到修改稿)

1 引言

近年来,光滑粒子流体动力学(SPH)方法作为一种纯Lagrangian型的无网格算法得到了很快发展,该算法最早是由Lucy[1]与Gingold[2]分别提出,最初用于解决天体物理中无边界三维流动问题,现已推广应用到计算力学诸多领域的数值仿真之中.该算法的特点是在计算中不需要生成任何辅助网格,因此在处理大变形问题时不存在Lagrangian网格算法中网格缠绕的技术瓶颈,理论上可以处理任意的变形问题;并且由于算法的Lagrangian特性,SPH离散粒子自然地追踪流体物质的运动及边界的变化,不需要像Eluerian网格算法中采取复杂的界面追踪技术,不存在Eluerian网格算法中因数值耗散造成界面追踪精度的降低.由于SPH算法的以上特点,使得它特别适合模拟自由表面、运动界面及大变形问题[3−6].

SPH方法发展还尚未成熟,尤其是固壁边界条件的施加方法,一直是困扰SPH方法发展的难点之一.早期Monaghan[7]模拟高速自由液面流动时引入了Lennard-Jones对势的强排斥力,其中固壁边界用虚粒子来描述,该模型中虚粒子对邻近边界的粒子强制施加排斥力,从而阻止邻近边界的粒子非物理穿透边界.但是该方法动量不守恒,对邻近边界的粒子的描述与实际不符.Libersky等[8]和Randles等[9]首次提出镜像粒子的方法来施加固壁边界,在边界之外布置关于边界对称的镜像粒子,该粒子与真实粒子密度相同,速度相反,并参与计算,该方法守恒性好但镜像粒子生成复杂,无法处理复杂壁面边界的问题.Gotoh等[10]建议固壁边界采用边界粒子平衡内部压强来防止粒子穿透边界,界面粒子采用多层镜像粒子来实现滑移边界条件,但仍难以处理复杂边界问题.文献[11,12]采用两种类型的虚粒子来处理固定边界条件,即一种虚粒子设置在固定边界上,与Monaghan[7]所用的相似,另一种虚粒子分布在边界的周围,这与Libersky等[8]的镜像粒子相似,即同时使用两种虚粒子来阻止实粒子穿透边界,但是这一方法也继承了两类虚粒子方法的不足.为了能够处理复杂壁面边界问题,近期国内外学者把研究的重点放在如何更好地施加排斥力上,先后提出了几种排斥力公式.Rogers等[13]提出了一种排斥力方法,该方法能很好地防止流体穿过壁面,但不能保证动量守恒;强洪夫等[14]基于Galerkin形式的SPH方程推导的罚函数边界力能较好地保证动量守恒,但该方法含有未知参数,需要经过试算才能确定其取值,且边界附近的流体粒子存在运动不稳定问题;Monaghan等[15]和Liu等[16]也分别提出了各自的排斥力施加方法,这两种方法类似,都能处理复杂边界,但都需要将边界粒子尺寸设为流体粒子尺寸的一半或更小,从而导致边界粒子数目的增加,降低了计算效率,且这两种方法不能保证动量守恒.

为了既保持排斥力方法能处理复杂边界的优势,又能克服上述几种方法的不足,本文提出一种新的排斥力模型来施加固壁边界条件.分别通过静止液柱算例和液柱坍塌算例,验证了各种排斥力方法的动量守恒特性及边界粒子尺寸对排斥力施加效果的影响;当仅采用边界上单层粒子计算边界力时存在边界粒子缺失问题,本文采用了Liu等[16]提出的耦合边界法中的虚粒子处理方法,以提高边界附近流体粒子的计算精度,将该法应用到容器中液体的静止算例及溃坝算例.对容器中液体的静止算例,计算得到了容器底部压强变化过程,通过与Monaghan等[15]与及强洪夫[14]所提出的方法进行对比,本文方法可有效克服边界附近的流体粒子存在的运动不稳定问题;通过溃坝算例验证了本文方法处理复杂边界问题的有效性.

2 SPH的基本方程及离散

2.1 控制方程

拉格朗日描述的流体动力学方程为:

参数P0可通过P0=100ρ0V/γ获得,其中ρ0为流体初始密度,γ=7,Vmax为流体的最大速度.为保证流体的弱可压缩性,一般取声速c=10Vmax.

2.2 流体动力学方程的离散

在SPH方法中,连续的流场离散成为一系列相互作用的粒子,通过核函数估计技术在这些粒子上离散控制方程组,得到一组描述各粒子物理量随时间变化的常微分方程组,即SPH基本方程组,再对这组方程采用相应的常微分方程组求解方法来推进时间进程的求解.

与耗散粒子动力学方法类似[17],求解域Ω被离散为一系列质点,其中每一个质点i记录了流场在该点处的物理量,如密度ρi,速度vi,位置xi等,并且给予一个统计意义的质量mi,每一个粒子同时拥有一个影响域S,影响域内所包含的其他SPH粒子 j为该粒子的邻近粒子,对于系统中任意点x处的物理量及其空间导数通过核函数估计得到,详见文献[3,4].

这样,流体运动的控制方程(1)—(3)可离散为

2.3 时间积分

形成虚粒子,从而导致其缺乏通用性,很难推广到复杂边界区域;排斥力法一般仅需要在边界处设置边界粒子,避免了在计算区域外设置虚粒子的问题,因此具有很强的适应性,但需要构造合理的排斥力公式来保证边界力的计算精度及稳定性;耦合边界法是排斥力法与虚粒子法的耦合,因此,合理的排斥力模型是耦合边界法的基础之一.本文着重研究以下5种典型的排斥力法,流体粒子与边界粒子的相互作用如图1所示.

本文采用leap-frog格式的时间积分方法,即

ϕ表示密度ρ或速度v,xi是粒子i的位置坐标.为了使计算过程稳定,本文采用Monaghan[7]给出的分别考虑具有黏性耗散和外力作用的时间步长表达式:

图1 流体粒子与边界粒子相互作用

方法1

Monaghan[7]最先提出用一组固定在边界上的粒子对邻近的流体粒子施加排斥力,从而防止流体粒子非物理穿透边界,该排斥力类似于计算分子力时的Lennard-Jones方程,表达式如下:

其中f是作用在单位质量上的力,ci代表当地声速.

2.4 传统排斥力模型

自由滑移边界条件的定义如下:

为了提高计算的精度及稳定性,国内外学者提出了多种方法来施加壁面边界条件.这些方法大致可分为三类:1)排斥力法;2)虚粒子法;3)耦合边界法.其中虚粒子法需要将流体粒子映射到边界外部

其中参数n1和n2一般取值为12和4,D一般取为最大速度的平方.r0为截止半径,若其取值过大,会使邻近边界的流体粒子受到较大的边界力,导致计算失败;若取值较小,就很难阻止流体粒子穿透边界.

方法2

Rogers等[13]提出的排斥力表达式为

其中nj为边界粒子 j处的边界法向,Ψ为流体粒子i到边界的距离,ξ为xij在粒子 j处切线方向的投影,u⊥为粒子i的速度在nj上的投影.排斥力函数R(Ψ)的表达式为

其中∆b为两邻近边界粒子的距离,P(ξ)可保证粒子平行于边界运动时,排斥力保持不变.ε(z,u⊥)是修正函数,可根据流体的深度及速度调整边界力的大小:

其中

式中h0为初始流体深度,z为流体局部深度,c0为初始声速.

方法3

为克服Lennard-Jones方程计算边界力时的不足,Monaghan等[15]又提出了一种边界力公式:

式中K取值为gh0,β一般取4,mi为粒子i的质量.

方法4

强洪夫等[14]基于Galerkin形式的SPH方程推导出了新的排斥力公式:

其中Aj是与边界粒子 j有关的边界权函数,ε为罚参数,其选取是个较困难的问题,一般需要若干次试算.由于该方法施加的边界力大小与流体相对于边界法向的速度大小有关,当流体粒子靠近边界且相对速度很小时,会导致边界力振荡,从而导致靠近边界的流体粒子速度、压强和密度的不稳定.

方法5

与Monaghan提出的边界力公式相似,Liu等[16]也提出了一种排斥力方法:

上述5种方法各有优缺点,下文将通过数值算例进行说明.

3 新型排斥力模型

如引言所述,国内外学者已提出多种排斥力方法,但都存在各自的不足之处.为了克服上述问题,本文基于Galerkin加权余量法并结合上述排斥力方法,提出了一种新型排斥力公式(方法6)来施加壁面边界条件.

Lagrangian描述的动量方程的张量形式如下:

vα是速度分量,σαβ是总应力张量分量,gα为单位质量体积力分量,以上公式中的重复上标符合Einstein求和约定.应用Galerkin加权余量法,动量方程(27)式可写为如下的弱形式:

式中右侧第一项为内力项,第二项为表面力项即边界力项,第三项为体力项,其中Ω代表整个计算域空间,Γ为计算域边界区域,nα为边界的外法线方向分量,δα代表权函数(test function)分量,代表试函数(trial function)分量.

Galerkin方法的权函数及其梯度的表达式如下:

其中δ(x−xj)为Diracδ函数.

方程(28)右侧的第二项表示边界力对动量方程的贡献,应用罚函数法对该项进行变换.首先根据滑移接触条件(14)来定义边界的约束罚势ΠB,其张量形式的表达式为

式中ε为罚参量,罚势的一阶变分形式为

将该罚势的变分形式引入到动量方程的弱形式(28)式中,并代替边界作用项,得到如下的形式:

将权函数和试函数的近似表达式(29),(30)和(31)代入(34)式,得到如下形式:

由于δvα的任意性,因此对于每个粒子必须满足下式:

应用核函数的性质并对上式进行点积分变化得到如下表达式:

通过将(39)式与(24)式对比可确定罚参数ε及Hij的取值.为使(39)式计算所得的边界力与(24)式具有相同的数量级,本文取

当流体靠近壁面时,(39)式所得的边界力与(vi−j)·nj成正比,为了避免方法4产生的压力振荡问题,本文对该项也进行了修正,最终得到的边界排斥力的定义如下:

4 数值算例

本节通过4个数值算例来验证本文方法的有效性,并与其他方法进行了对比.为了保证计算结果的可对比性,每种方法都采用相同的粒子配置,核函数及其影响域也保持一致.

分别通过静止液柱算例及液柱坍塌算例,对比了传统排斥力公式及本文新型排斥力公式的优缺点,证明了本文方法的有效性.

对容器中液体静止算例及溃坝算例,当固壁粒子仅采用边界上的单层粒子分布时,会造成边界附近粒子缺失,从而导致边界附近的流体粒子计算精度降低,本文采用Liu等[16]提出的虚粒子施加方法来解决该问题.即在边界粒子外侧设置两排虚粒子,边界粒子的信息通过其邻近的流体粒子插值得到,虚粒子的信息通过边界粒子和流体粒子的信息插值获得,如(41)—(43)式:

通过液体静止算例,验证了本文方法所得的边界力与真实值的一致性;通过溃坝算例,验证了本文方法处理复杂边界问题的有效性.

4.1 静止液柱算例

为了验证边界力对流体粒子初始分布的扰动特性,本节选取静止液柱算例进行数值分析.

该算例中,液柱不受重力作用.其尺寸为0.4 m×0.6 m,边界长度为1 m.核函数统一采用三次样条函数,粒子i的影响域hi=1.0∆di,∆di为粒子i的初始尺寸.初始密度ρ0=1000.0 kg/m3,边界粒子与流体粒子尺寸大小相同,取为0.02 m,时间步长为1×10−5s,初始粒子分布如图2,最底部的一层粒子为边界粒子,其余粒子为流体粒子.由于流体静止且不受外力作用,要保证系统动量守恒,在计算过程中流体粒子应始终保持静止.

t=0.3 s时,采用不同方法模拟的结果如图3,图中箭头方向为粒子的运动方向.从图3可知,方法4及本文方法能保证流体粒子始终静止.

图2 静止液柱算例初始SPH粒子分布(单位为m)

图3 t=0.3 s时用不同方法模拟静止液柱获得的SPH粒子分布及速度分布(单位为m) (a)方法1;(b)方法2;(c)方法3;(d)方法4;(e)方法5;(f)方法6

方法1的模拟结果和截止半径r0及影响域h的取值有关.若r0＜h/2,则流体粒子保持静止,否则流体粒子将受到边界力的扰动.本文取r0=1.1∆d＞h/2,邻近边界的流体粒子由于受到边界粒子排斥力的作用而产生运动;若增大r0的取值,流体粒子所受的排斥力增加,流体粒子的运动速度也会增大.选取四次样条核函数及高斯核函数也存在同样的问题.因此,截止半径r0的取值问题不利于方法1在实际问题中的应用.

方法2、方法3及方法5中,无论选取何种核函数,与边界粒子邻近的流体粒子受到的边界排斥力均不为0,边界力的作用使流体粒子产生运动,因此不能保证动量守恒;且方法2中边界排斥力最大,流体粒子运动速度明显大于其他方法所得的流体粒子运动速度.

方法4的罚参数取值为1.0.在方法4及方法6中,排斥力大小与流体和壁面之间的相对速度有关,仅当流体靠近壁面时,才受到边界力的作用.这就使初始静止流体能始终保持静止,保证了邻近边界的流体粒子分布与实际相符.

若核函数统一采用高斯核函数或四次样条函数,影响域取初始粒子间距的1.0—1.5倍,计算结果与上述结果相似,所得的结论保持不变.

4.2 液柱坍塌算例

液柱在重力作用下会向下运动,排斥力的存在将阻止其穿透壁面,本算例对各种排斥力阻止流体粒子非物理穿透壁面的特性进行研究.

核函数及其影响域、液柱尺寸、初始设置与上节相同,但受重力作用,重力加速度g=9.81 m/s2,沿y轴负方向,

t=0.3 s时的计算结果如图4所示.从图4中可知:当边界粒子和流体粒子尺寸大小相同时,方法1、方法3、方法4及方法5并不能阻止流体粒子穿透边界,究其原因,是由于这4种方法的排斥力大小与边界粒子和流体粒子之间的距离成反比.要防止流体粒子非物理穿透边界,有两种方法:一种是将边界粒子尺寸设为流体粒子的1/2或更小,如文献[15],但是这样会增加边界粒子数目,降低计算效率,尤其是三维问题;另一种是使排斥力大小与流体粒子到边界粒子的法向距离成反比,如方法2和方法6,这种处理方法可减小边界力精度对边界粒子尺寸的依赖,提高计算效率.

方法4中,罚参数ε选取需要经过若干次试算才能获得.若罚参数过大,流体粒子在向下运动的过程中会受到过大的边界力,从而导致该粒子加速远离边界运动,扰乱整个流场的粒子秩序,致使模拟失效;若罚参数过小,边界排斥力很难阻止流体粒子穿透边界.本算例经过3次试算,分别取罚参数为0.1,1.0,10.0,最终确定罚参数取1.0.

4.3 容器中液体静止算例

为了验证本文方法所得的排斥力与边界力的真实值是否一致,选取文献[13]的容器中液体静止算例进行检验.本节仅选取方法3、方法4与本文方法进行对比,来验证本文方法的有效性.

图4 t=0.3 s时,不同方法计算液柱坍塌得到的SPH粒子分布 (a)方法1;(b)方法2;(c)方法3;(d)方法4;(e)方法5;(f)方法6

容器长2 m,高1 m,水深H=0.9 m.方法3中,计算参数与文献[15]一致,为防止流体粒子穿透边界,边界粒子尺寸设为流体粒子的1/2;方法4中,罚参数设为2.0,流体粒子尺寸设为0.02,边界粒子尺寸设为0.015 m;本文方法中所有粒子的尺寸都设为0.02 m.3种方法的时间步长同取2.0×10−5s,所有流体粒子的初始压强都设为0,初始密度ρ0=1000.0 kg/m3,重力加速度g=9.81 m/s2,沿y轴负方向

在重力作用下,流体粒子向下运动,但边界力的存在阻碍了其进一步运动,这样流体粒子就会在重力和边界力的共同作用下振荡运动,随着时间的增加,流体粒子所受的力会达到平衡,容器底部的流体粒子的压强理论值应为ρgH,H为液体深度.

t=1 s时,粒子分布状态如图5.方法3的计算结果如图5(a)所示,可知高度越高,靠近边界的流体粒子离边界越远,且该距离大于初始粒子间距,这就导致左上角及右上角处流体粒子的聚集及粒子秩序的紊乱.方法4及本文方法中,边界力与粒子间的相对速度有关,从而克服了方法3存在的问题,邻近边界的流体粒子具有更好的分布秩序.

图5 t=1.0 s时三种不同方法模拟的SPH粒子分布状态 (a)方法3;(b)方法4;(c)方法6

图6 容器底部中央的SPH粒子相对压强随时间变化过程

以容器底部靠近中央的流体粒子为研究对象,记录其压强随时间的变化过程如图6.图6中压强为其相对值,即p/ρgH,液体趋于静止时,相对压强的理论值为1.0.从图6可知,方法3及本文方法在t=0.6 s时压强都达到了稳定状态,相对压强都趋于1.0.本文方法计算所得的相对压强略小于1.0,主要是由人工黏性产生的耗散所导致的.采用方法4时,计算所得的相对压强持续振荡,t=0.6 s时压强基本趋于稳定,但t=1.02 s时,压强又开始剧烈振荡.其原因是:当邻近边界的流体粒子压强趋于稳定时,其运动速度趋于0,从(16)式可知,边界力亦趋于0,流体粒子在压力作用下会产生较大的加速度及速度向边界靠近,此时在(16)式作用下又会产生较大的边界力使流体粒子远离边界,这样邻近边界的流体粒子就会产生振荡的运动速度,由(5)式及(4)式可知,粒子运动速度的剧烈变化会导致粒子密度及压强的大幅振荡.本文方法则克服了方法4存在的缺陷.

图7 采用本文方法得到的不同时刻的相对压强分布 (a)t=0.8s;(b)t=1.0s;(c)t=1.2s

为了考察压力场的空间分布特性,图7给出了0.8,1.0,1.2s三个不同时刻的相对压强分布状态.由图7可知,本文方法计算所得的压力场比较光滑;容器底部压力场稳定后的不同时刻,压力场的空间分布一致且不随时间变化.

4.4 溃坝算例

溃坝是经典的自由表面流动问题,常用来验证壁面边界施加方法的有效性.计算模型如图8,水柱位于容器左侧,右侧有一挡板,初始时刻水柱静止,然后将挡板迅速抽出.水柱高度H=2L,宽度为L=0.146m,容器尺寸为2L×4L.

图8 溃坝几何模型

本文对溃坝过程进行了模拟,并与Koshizuka等[18]的实验以及方法4、方法5所得的结果进行了对比.计算参数为:初始密度ρ0=1000.0 kg/m3,边界粒子数目为805,边界粒子初始尺寸为1.5mm,流体粒子数目为5000,流体粒子初始间距为 2.92mm,时间步长 1×10−5s,Vmax=方法4的罚参数经试算后取0.5.核函数统一采用三次样条函数,粒子i的影响域hi=1.0∆di.为克服边界粒子缺失问题,本算例的3种方法都采用了Liu等[16]提出的虚粒子法来提高边界附近的计算精度.

3种边界施加方法所得的模拟结果如图9,第1行为0.2,0.6s时的实验结果;第2行为方法4所得的结果;第3行为方法4所得的结果;第4行为采用本文方法所得的模拟结果.图9中箭头方向表示粒子的运动方向,其长短表示粒子运动速度的大小.

从图9中可知,3种方法所得的自由表面运动过程与实验结果都比较符合.但是,采用方法4施加自由滑移边界条件时,靠近边界处部分流体粒子的速度明显大于其周围区域流体粒子的速度;随着时间的增加,这种非物理的速度振荡并不会消失,且边界流体粒子速度失真会导致粒子密度和压强的不稳定,最终导致边界附近流场的数值模拟失真.采用方法5及本文方法时,靠近边界的流体粒子始终保持平行壁面边界运动,粒子速度场分布光滑且秩序良好,很好地克服了方法4所导致的速度振荡问题.

5 结论

本文提出一种新型排斥力模型,该模型不需要确定未知参数,可较好地防止流体粒子非物理穿透壁面,能克服流体相对壁面速度较小时速度不稳定问题.通过液柱静止算例、液柱坍塌算例及容器中液体静止算例对比并验证了本文方法的动量守恒性、防止流体粒子非物理穿透边界的特性及边界力的精确度;最后,将本文方法应用到溃坝算例,进一步验证了该方法适应复杂边界问题的能力.

本文方法可处理复杂形状边界,不需要确定未知参数,可较好地施加固壁边界条件.本文方法还得出:边界排斥力的大小要与流体粒子到边界的法向距离成反比,而不是与流体粒子到边界粒子的距离成反比,此时才能保证边界力不受边界粒子尺寸的影响,有效防止流体粒子非物理穿透壁面边界;排斥力模型中需要考虑流体与壁面的相对速度对排斥力的影响,否则会导致边界附近流体粒子速度失真.