基于长程力的SPH方法固壁边界处理

周学君，陈 丁，黄文雄

(1.河海大学力学与材料学院，江苏 南京 210098； 2.黄冈师范学院数理学院，湖北 黄冈 438000)

基于长程力的SPH方法固壁边界处理

周学君1，2，陈 丁1，黄文雄1

(1.河海大学力学与材料学院，江苏 南京 210098； 2.黄冈师范学院数理学院，湖北 黄冈 438000)

将近场动力学(PD)中的长程力引入光滑粒子流体动力学(SPH)的固壁边界处理模型中，给出一种适合SPH方法边界处理的新方法。通过液柱坍塌、Poiseuille流、水铸造充型仿真和土体崩塌4个算例对该方法的有效性进行验证。结果表明：该方法能很好地解决粒子非物理穿透边界；能客观反映真实的边界作用；对于复杂几何边界和具有材料强度的动力学问题也相当有效。

光滑粒子流体动力学(SPH)；长程力；近场动力学(PD)；固壁边界处理

光滑粒子流体动力学(smoothed particle hydrodynamics，简称SPH)是一种纯Lagrangian形式的无网格粒子计算方法，由Lucy[1]和Gingold等[2]分别独立提出。经过不断地研究和改进，SPH方法的应用领域已由早期的天体物理学和流体力学相关的问题，延伸到具有材料强度的固体力学和生物力学等范畴。SPH方法在模拟大变形瞬时动力学问题时，不存在网格类算法(如FEM)中因出现网格畸变或缠绕而导致的算法精度降低和失败等难题[3]。

尽管SPH方法在许多问题的应用中已经取得巨大成功，但对于边界条件的处理仍然存在着一定的困难。这是因为SPH方法不具有Kronecker性质[4]，即在边界上或者邻近边界处进行积分时会被边界截断，导致计算不能完全覆盖整个区域。为了解决SPH方法的边界处理问题，人们基于不同的思路提出了一些解决方案[5-14]。目前，SPH方法施加边界条件，或是需要在边界上布置镜像虚粒子，或是在边界上布置边界虚粒子。前者虽然守恒性比较好，但镜像粒子的生成比较麻烦，特别是固壁边界形状不太规则时更是不易；后者虽然不受固壁边界形状的影响，容易实现，但虚粒子对实粒子的排斥力是人为给定的，没有统一的排斥力模型，并且参数也不好确定。笔者针对后一种情况，建议在布置边界虚粒子的基础上，利用近场动力学(peridynamics，简称PD)中描述物质点之间相互作用的本构力概念，将非局部的长程作用力模型引入到边界条件中。

PD方法是由美国Sandia国家实验室的Silling[15]提出的一种新兴的非局部、无网格物质点类方法[16-17]，已经在固体材料和结构的静动力变形以及破坏分析中成功应用[18-20]。PD方法中，用来描述“近场”范围内两物质点间相互作用的长程力已经被众多学者进行了深入细致的研究，已取得一些有效的模型[15，21-24]。另一方面，SPH方法的固壁边界处理中边界虚粒子与实粒子之间的作用力，类似于PD方法中物质点间的本构力。因此，笔者将长程作用力的概念和模型引入SPH方法的边界处理中，并通过编程和具体算例对该方法处理SPH固壁边界的有效性进行验证。

1 SPH方法的基本思想和控制方程

1.1 基本思想

SPH方法将问题域离散成有限个具有独立质量和体积的粒子，每个粒子都携带有速度、密度、应力等材料特性，且通过其影响域(支持域)内邻近粒子进行插值近似来计算它们的材料特性。SPH方法近似一般分2步进行，第1步是将函数按照积分加权的形式进行积分近似，第2步是将问题域和积分进行离散化(粒子近似)。

对于定义在问题域Ω内任意点x处的场函数f(x)，采用如下积分表达：

(1)

其中

式中：x、x′∈Ω——位置坐标向量；W——光滑函数或者核函数，采用三次样条函数[3]为光滑函数；h——光滑长度，表征光滑函数的影响区域；ad——正则化因子；R——点x和x′处两粒子之间的相对距离。

在将问题域采用有限个粒子离散化后，在粒子i处的场函数f(xi)的SPH积分表达采用以下离散化后的粒子近似式：

(2)

其中

Wij=W(xi-xj，h)

式中：mj、ρj——粒子j的质量和密度；N——粒子i支持域内的粒子总数。

1.2 控制方程

Lagrangian描述下的无黏性流体动力学控制方程包括连续性方程、动量方程和运动方程，应用SPH核近似和粒子近似后，离散化的SPH公式为

(3)

其中

式中：v、p[6]——速度和压力；g——体力；Πij——人工黏度项，它的作用是消除SPH方法在模拟流体动力学问题时产生的数值不稳定性，本文算例采用了文献[25]中建议的人工黏度表达式；ρ0——流体初始密度；γ——参数，取7；vmax——流体的最大速度。

SPH方法模拟固体大变形问题时，动量方程要采用如下形式：

(4)

式中：σ——总应力张量；l——单位张量

对于土体，采用基于Drucker-Prager失效条件的理想弹塑性本构模型描述，SPH离散后的公式为[26]

(5)

其中

2 边 界 处 理

对于固壁边界的处理，先在边界上布置一定数量的虚粒子，然后将PD方法中的本构力引入模型中，用作边界虚粒子对实粒子的作用力，达到施加边界条件的目的。

在PD方法中，为了描述在“近场”范围内物质点之间相互作用力，给出本构力的概念。本构力是一种非局部长程作用力，Huang等[23]提出的本构力表达式为

(6)

其中

式中各项意义见文献[23-24]。

图1 边界虚粒子对实粒子的作用力

(7)

其中

式中：rij——i与j之间的距离；a——常数，表征边界力的作用范围与光滑长度之间的数量关系，常取0.5～1。

式(7)作用力的方向与(ξ+η)相同，固壁边界虚粒子i的位置可以认为保持不变，其位移u=0，所以f的方向为[(x′+u′)-x]的方向，如图1。当边界虚粒子的影响域内有多个实粒子，它们之间的排斥力可以直接叠加得到。

边界排斥力的表达式(式(7))形式简单、没有过多的参数，但传统SPH方法中边界条件施加过程中罚函数的参数较多且很难有实际意义，参数值的确定受人为因素影响较大。

3 数 值 算 例

通过4个算例来验证本文所讨论的边界处理方法的有效性，其中前3个算例为流体动力学模型，最后一个为固体动力学模型。在前3个流体模型中，假设流体具有弱可压缩性，流体密度为103kg/m3，边界排斥力表达式(7)中，体积弹性模量K可以通过式(3)中p的表达式和Navier-Stokes连续性方程得到。在固体模型中，以土体崩塌为算例，并与镜像虚粒子法的模拟结果比较，此时式(7)中参数的取值直接采用土体相关的材料参数。4个算例中，模型均视为平面应变问题，设定式(7)中的δ与光滑长度h、粒子的初始间距相同，时间积分采用跳蛙法。

众所周知，SPH方法积分近似在边界上会因支持域的截断导致计算精度降低。事实上，当固壁边界仅采用单层粒子分布时， 会造成边界附近粒子缺失，从而导致积分近似计算的精度降低；本文采用设置多层(3层)边界粒子的办法来最大限度地减小其数值影响，即在边界粒子的外侧再增加设置两排虚粒子，且这些边界虚粒子的密度、速度等信息，分别通过其邻近的实粒子的信息插值计算。

3.1 液柱坍塌

为了验证本文提出的边界处理方法对阻止粒子非物理穿透固壁时的表现，选取液柱坍塌算例进行模拟研究。如图2所示，尺寸为0.4 m×0.6 m的液柱放置在无限长的边界上，在自身重力作用下，液柱向下坍塌。流体粒子和边界粒子的初始间距均为0.02 m，共需要600个流体粒子，时间步长为1×10-5s。因为液柱受重力作用的最大速度vmax=3.43m/s，所以参数B=1.68×105，式(7)中β=1.0、a=1。

图3展示了t=0.3 s时的模拟结果，可以看出坍塌后流体粒子完全附着在边界之上，未见流体粒子有穿透边界的情况发生。在其他的排斥力边界处理方法中，常用减小边界虚粒子间距的办法来阻止粒子非物理穿透固壁[13]，但在本文所研究的模型中只需要求边界虚粒子分布与实粒子的初始粒子分布保持一致，即可保证模拟效果的有效性和正确性，且控制了边界虚粒子的数目，提高了运算效率。

图2 液柱坍塌的初始SPH 粒子分布

图3 t=0.3 s时液柱坍塌的SPH粒子分布

图4 Poiseuille流的初始SPH 粒子分布

3.2 Poiseuille流

Poiseuille流是流体动力学的经典模型，可以用来比较边界方法中排斥力与真实边界力的一致性。如图4所示，在该模型中，初始静止的流体受到外力F的作用在2块间距为l的平行且无限大的固定平板间流动，最终达到稳定状态。流体的尺寸是0.000 5 m×0.001 m，流体粒子的数量为20×40。

Poiseuille流任一点的流动速度可以通过求解Navier-Stokes动量方程得到，Morris等[27]已经推导出Poiseuille流的流体速度的理论级数解。由Morris等[27]的数解可得流体前沿的最大速度vmax=1.25×10-5m/s，于是相应的Re=1.25×10-2，B=2.23×10-6，式(7)中β=1.0，a=1。

在SPH模拟中，时间步长为1×10-5s，以流体最前端一列间隔选取的20个粒子的速度为参照，经过5 000步计算后流体的速度达到稳定状态。图5给出分别由SPH方法和Morris等[27]的理论级数解得到的在t=0.01 s、0.1 s、0.2 s和最终稳定状态t=∞时坐标y与速度的函数图像，可以看出两者吻合的程度很高，且本文SPH模拟值与理论值的最大相对误差为1.8%，比文献[27]中给出的两者最大相对误差2%略低，表明本文研究的边界方法中的排斥力模型能客观地反映真实边界力效果，且能有效处理平行于壁面的流动问题。

图5 Poiseuille流的速度分布(Re=1.25×10-2)

图6 “弓”字型型腔示意图(单位：mm)

3.3 水模拟铸造充型

水模拟铸造充型来自Schmid等[28]的试验，用来测试本文边界施加方法对于复杂几何边界的处理能力。几何尺寸如图6所示，整体上来看该型腔为“弓”字型，含有5个直角转角，2个弧形转角，底部入口处有一垂直延伸段以便水能够被腔壁约束顺利入腔。试验过程中水从型腔底部入口以8.7 m/s的速度充入腔内，不计重力影响。

SPH模拟中共需要设置4 290个虚粒子和27 000个水粒子，粒子的初始间距为0.001 m，时间步长为2×10-7s，参数B=1.08×106、β=1.0、a=1。

图7 水铸造充型模拟的试验结果与模拟结果对比

试验、文献[12]模拟和本文SPH模拟结果如图7所示，共选取4个时间节点的型腔内的流场形态进行比较，分别是7.15 ms、25.03 ms、39.34 ms和53.64 ms。由图7可知，虽然腔壁边界的几何形状复杂，但文献[12]和本文所研究的边界处理方法有效地施加了边界条件，使得流体完全控制在腔内沿腔壁运动，未出现流体非物理穿透腔壁的现象。本文SPH模拟结果比较精确地预测流体在不同时刻的流动状态及各个转角处空腔的形成和消失过程。另一方面，由于试验与模拟过程中某些客观条件无法完全一致，如空腔内气体压力和流体自身重力影响等，所以流体形态在细节方面与试验存在着一定差异，但是这种差异在其他的边界处理方法的SPH仿真水模拟铸造充型结果中也存在[12，14]。文献[12]中利用罚函数施加边界条件，罚函数中参数较多且难以设定，而本文提出的边界处理方法中模型简单且参数少、易于设定。

为进一步验证本文边界处理方法对流场的影响，图8描述了在SPH模拟水铸造充型过程中粒子的最大速度随时间变化的图像，可以看出全程粒子的最大速度出现在5.7 ms时刻，大致上在水流通过第1个弧形转角后获得，图中出现的最大速度峰值几乎都出现在水流通过转角时刻，这些预测与文献[14]的结果吻合。

图8 SPH模拟水铸造充型过程中的粒子最大速度随时间变化曲线

3.4 土体崩塌

SPH方法模拟流体动力学问题时，本文提出的排斥力模型能有效地模拟流体与固壁间的边界力作用。本算例将该模型应用到固体大变形问题上，用来验证边界排斥力方法在具有材料强度的动力学问题的有效性。

土体崩塌是SPH方法应用到固体大变形问题的常用模型，Bui等[26]采用文献[7]中的镜像虚粒子的方法来处理边界，其模拟效果已经得到广泛认可，这里只需将本文所研究的边界处理方法的模拟结果与之作比较即可。

土体的原始尺寸为4 cm×6 cm，需要2 400个土体粒子，在边界上布置共计795个虚粒子，粒子的初始间距为0.001 m，粒子的密度为1 850 kg/m3。土体泊松比为0.3，内摩擦角φ=27.5°，剪胀角ψ=0，式(5)中其他参数K=8.3×105Pa、G=3.9×105Pa、αφ=0.15；式(7)中β=0.1、a=0.5；时间步长设为2×10-6s。

为了更好地观察模拟结果的细节，将土粒子平均分成3层，分别用不同颜色标注。图9给出了2种边界处理方法的土体崩塌最终效果，土块在历时0.4 s后完全停止运动，2种模拟结果土体崩塌后的形态相似，上端自由表面基本一致，静止角的大小几乎相等，均略小于内摩擦角φ，这也与文献[26]中的结论吻合。另外，从图9(b)来看，未见土体粒子非物理穿透边界，且土体实粒子分布规则有序，充分说明本文研究的边界处理方法对于土体的变形问题是有效的。

图9 SPH模拟土体崩塌

在镜像虚粒子法中，虚粒子的生成比较麻烦，特别是对于具有复杂几何形状的边界，计算量大，难以推广到2种变形体之间的边界力处理，而本文提出的边界处理方法易于实施，且预测效果真实，具有很好的推广性。

4 结 语

利用SPH方法模拟流体和固体大变形，重点研究固壁边界处理问题。在固壁上布置边界虚粒子的基础之上，考虑利用PD中描述物质点之间相互作用的本构力概念，将非局部的长程作用力模型引入到边界条件中。以液柱坍塌、Poiseuille流、水铸造充型仿真和土体崩塌为算例，验证本文所研究的固壁边界处理方法的有效性。

PD方法中的长程作用力与SPH方法中虚、实粒子的边界排斥力所需要满足的条件相似，且其表达式简单，含有的参数数目较少且易于给定。从数值模拟的结果来看，本文所采用的边界处理模型能够很好地解决粒子非物理穿透边界，能较客观地反映真实的边界作用，对于复杂几何边界和具有材料强度的动力学问题也相当有效。理论的契合和仿真的吻合，同时说明长程力作为SPH固壁边界处理的新思路是切实可行的，值得推广。

[ 1 ]LUCYLB.Anumericalapproachtothetestingofthefissionhypothesis[J].AstronomicalJournal， 1977， 82(12)： 1013-1024.

[ 2 ]GINGOLDRA，MONAGHANJJ.Smoothedparticlehydrodynamics:theoryandapplicationtononsphericalstars[J].MonthlyNoticesoftheRoyalAstronomicalSociety， 1977， 181(2)： 375-389.

[ 3 ]LIUMB，LIUGR.Smoothedparticlehydrodynamics(SPH)：anoverviewandrecentdevelopments[J].ArchComputMethodsEng,2010，17： 25-76.

[ 4 ]LIUMB，LIUGR.Smoothedparticlehydrodynamics：ameshfreeparticlemethod[M].Singapore：WorldScienicPublishingCoPteLtd， 2003.

[ 5 ]CAMPBELLPM.Somenewalgorithmsforboundaryvaluesproblemsinsmoothedparticlehydrodynamics[R].Albuquerque：MissionResearchCorp， 1989.

[ 6 ]MONAGHANJJ.SimulationfreesurfaceflowswithSPH[J].JournalofComputationalPhysics， 1994， 110： 399-406.

[ 7 ]LIBERSKYLD，PETSCHCKAG，CARNEYTC，etal.HighstrainLagrangianhydrodynamics：athree-dimensionalSPHcodefordynamicmaterialresponse[J].JournalofComputationalPhysics， 1993， 109： 67-75.

[ 8 ]RANDLESPW，LIBERSKYLD.Smoothedparticlehydrodynamics：somerecentimprovementsandapplications[J].ComputerMethodsinAppliedMechanicsandEngineering， 1996， 138： 375-408.

[ 9 ]LIUMB，LIUGR，ZONGZ.Anoverviewonsmoothedparticlehydrodynamics[J].InternationalJournalofComputerMethods， 2008， 5(1)：135-188.

[10]LIUMB，LIUGR，LAMKY.Investigationsintowatermitigationsusingameshlessparticlemethod[J].ShockWaves， 2002， 12(3)：181-195.

[11] 龚凯，刘桦，王本龙.SPH固壁边界处理方法的改进[J]. 力学季刊， 2008， 29(4)， 507-514. (GONGKai，LIUHua，WANGBenlong.AnimprovedboundarytreatmentapproachforSPHmethod[J].ChineseQuarterlyofMechanics， 2008， 29(4):507-514. (inChinese))

[12] 强洪夫，韩亚伟，王坤鹏，等. 基于罚函数SPH新方法的水模拟充型过程的数值分析[J]. 工程力学， 2011， 28 (1)： 245-250.(QIANGHongfu，HANYawei，WANGKunpeng，etal.NumericalsimulationofwaterfillingprocessbasedonnewmethodofpenaltyfunctionSPH[J].EngineeringMechanics， 2011， 28 (1)： 245-250. (inChinese))

[13] 韩亚伟，强洪夫，赵玖玲，等. 光滑粒子流体动力学方法固壁处理的一种新型排斥力模型[J]. 物理学报，2013， 62(4)：326-336. (HANYawei，QIANGHongfu，ZHAOJiuling，etal.Anewrepulsivemodelforsolidboundaryconditioninsmoothedparticlehydrodynamics[J].ActaPhysicaSinica，2013， 62(4)：326-336.(inChinese))

[14] 刘虎，强洪夫，陈福振，等. 一种新型光滑粒子动力学固壁边界施加模型[J]. 物理学报，2015， 64(9)：375-388. (LIUHu，QIANGHongfu，CHENFuzhen，etal.Anewboundarytreatmentmethodinsmoothedparticlehydrodynamics[J].ActaPhysicaSinica，2015， 64(9)：375-388.(inChinese))

[15]SILLINGSA.Reformulationofelasticitytheoryfordiscontinuitiesandlong-rangeforces[J].JournaloftheMechanicsandPhysicsofSolids， 2000， 48(1)： 175-209.

[16]SILLINGSA，EPTONM，WECKNERO，etal.Peridynamicstatesandconstitutivemodeling[J].JournalofElasticity， 2007， 88(2)： 151-184.

[17] 黄丹， 章青， 乔丕忠， 等. 近场动力学方法及其应用[J]. 力学进展， 2010， 40(4)： 448-459. (HUANGDan，ZHANGQing，QIAOPizhong，etal.Areviewonperidynamicsmethodanditsapplication[J].AdvanceinMechanics， 2010， 40(4)： 448-459. (inChinese))

[18]KILICB.Peridynamictheoryforprogressivefailurepredictioninhomogeneousandheterogeneousmaterials[D].Tucson：TheUniversityofArizona， 2008.

[19] 胡祎乐， 余音， 汪海. 基于近场动力学理论的层压板损伤分析方法[J]. 力学学报， 2013， 45(4)： 624-628. (HUYile，YUYin，WANGHai.Damageanalysismethodforlaminatesbasedonperidynamictheory[J].ChineseJournalofTheoreticalandAppliedMechanics， 2013， 45(4)： 624-628. (inChinese))

[20] 章青， 顾鑫， 郁杨天. 冲击荷载作用下颗粒材料动态力学响应的近场动力学模拟[J]. 力学学报， 2016，48(1)：56-63. (ZHANGQing，GUXin，YUYangtian.Peridynamicssimulationfordynamicsresponseofgranularmaterialsunderimpactloading[J].ChineseJournalofTheoreticalandAppliedMechanics， 2016，48(1)：56-63.(inChinese))

[21]WECKNERO，ABEYARATNER.Theeffectoflong-rangeforcesonthedynamicsofabar[J].JournalofMechPhysSolids， 2005， 53： 705-728.

[22]SILLINGSA，ASKARIE.Ameshfreemethodbasedontheperidynamicmodelofsolidmechanics[J].ComputStruct， 2005， 83：1526-1535.

[23]HUANGDan,LUGuangda,QIAOPizhong.Animprovedperidynamicapproachforquasi-staticelasticdeformationandbrittlefractureanalysis[J].InternationalJournalofMechanicalSciences， 2015，94/95：111-122.

[24]HUANGDan,LUGuangda,WANGChongwen.etal.Anextendedperidynamicapproachfordeformationandfractureanalysis[J].EngineeringFractureMechanics，2015， 141：196-211.

[25]LATTANZIOJC，MONAGHANJJ，PONGRACICH，et.al.Controllingpenetration[J].SIAMJournalonScienicandStatisticalComputing， 1986， 7(2)： 591-598.

[26]BUIHH，FUKAGAWAR，SAKOK，etal.Lagrangianmeshfreeparticlesmethod(SPH)forlargedeformationandfailureflowsofgeomaterialusingelastic-plasticsoilconstitutivemodel[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics，2008， 32(12)：1537-1570.

[27]MORRISJP，PATRICKJF，ZHUY.ModelinglowreynoldsnumberincompressibleflowsusingSPH[J].JournalofComputationalPhysics， 1997， 136： 214-226.

[28]SCHMIDM，KLEINF.Fluidflowindiecavities-experimentalandnumericalsimulation[R].Indianapolis：NADCA， 1995.

Solid boundary treatment for SPH method based on long-range force

ZHOU Xuejun1，2， CHEN Ding1， HUANG Wenxiong1

(1.CollegeofMechanicsandMaterials,HohaiUniversity,Nanjing210098,China;2.CollegeofMathematicsandPhysics,HuanggangNormalUniversity,Huanggang438000,China)

A new method for boundary treatment of the smoothed particle hydrodynamics (SPH) method is presented through introduction of the long-range force in peridynamics (PD) to the solid boundary treatment of SPH. The reliability of the method is verified by four numerical simulation examples, including a liquid column break, Poiseuille flow, water casting filling, and soil collapse. The results show that the method can solve the problem of non-physical particle penetration boundary and objectively reflect actual boundary action, and is effective for dynamic problems with complex geometric boundies and material strength.

smoothed particle hydrodynamics (SPH); long-range force; peridynamics (PD); solid boundary treatment

10.3876/j.issn.1000-1980.2017.02.009

2016-04-13

国家自然科学基金(11172088)；江苏省普通高校研究生科研创新计划项目(KYZZ16_0268)；黄冈师范学院高级别培育项目(201617603)

周学君(1981—)，男，湖北黄冈人，讲师，博士，主要从事计算力学与工程仿真研究。E-mail：zhouxj@hhu.edu.cn

O242

A

1000-1980(2017)02-0153-08