一种新型SPH固壁边界处理的排斥力模型

周学君，陈 丁，唐 轶

(1.黄冈师范学院 数理学院，湖北 黄冈 438000; 2.河海大学 力学与材料学院，南京 210098；3.云南民族大学 数学与计算机科学学院，昆明 650500)

一种新型SPH固壁边界处理的排斥力模型

周学君1，2，陈 丁2，唐 轶3

(1.黄冈师范学院 数理学院，湖北 黄冈 438000; 2.河海大学 力学与材料学院，南京 210098；3.云南民族大学 数学与计算机科学学院，昆明 650500)

边界排斥力法是光滑粒子流体动力学(SPH)固壁边界处理的方法之一，但由于缺乏统一的排斥力模型而制约其广泛应用。考虑将近场动力学(Peridynamics，PD)中描述颗粒间接触作用的短程排斥力引入到固壁边界处理模型中，提出一种新型SPH方法边界排斥力模型。通过Couette流和溃坝2个算例验证了排斥力模型的有效性。排斥力表达式简单，参数易于给定，为SPH方法中固壁边界处理提供新思路。

光滑粒子流体动力学(SPH)；排斥力模型；近场动力学(PD)；固壁边界; Couette流

1 研究背景

光滑粒子流体动力学(Smoothed Particle Hydrodynamics，SPH)是一种产生最早、发展最成熟的纯Lagrangian形式的无网格计算方法。SPH方法最初是由Lucy(1977)[1]，Gingold等[2]提出用来解决三维开放空间的天体物理学问题，之后经过不断的改进和完善，已经广泛应用到流体力学、固体力学和生物力学等领域。SPH方法利用有限数量的粒子来离散问题域，粒子之间通过核函数来建立联系，在模拟大变形问题时，不存在网格类算法(如FEM)中因出现网格畸变或缠绕而导致的算法精度降低和失败等难题[3]。

SPH方法在模拟自由表面流问题时，由于粒子的Lagrangian特性能够自动地获得自由表面边界，而不需要额外地添加边界条件，非常方便。但对于固壁边界条件的处理仍然存在着一定的困难，这是因为在边界上或者邻近边界处，计算会被边界截断导致不能完全覆盖整个区域，产生较大误差。目前使用较为广泛的SPH固壁边界处理方法有边界排斥力法、镜像虚粒子法和耦合边界粒子法，但各有优缺点。

边界排斥力法最早于1994年由Monaghan[4]提出，通过在固壁边界上布置一定数量的虚粒子，使之对邻近边界的真实粒子产生适当的排斥力，从而防止邻近边界的真实粒子非物理穿透边界。边界排斥力法虽然不受固壁边界形状的影响，容易实现，但排斥力是人为给定的，没有统一的排斥力模型，并且参数也不好确定。Libersky等(1993)[5]首创在边界上镜像布置虚粒子的方法来施加边界条件；在此基础之上，Randles等(1996)[6]对镜像虚粒子法进行改进，要求所有镜像虚粒子具有边界处等值的变量值，采用粒子近似将边界条件施加到边界粒子上。镜像虚粒子法虽然守恒性比较好，但镜像粒子的生成比较麻烦，特别是固壁边界形状不太规则时更是不易。Liu等(2002[7],2003[8])提出耦合边界粒子法来处理固壁边界，在固壁处布置多层的边界虚粒子，并对这些边界虚粒子赋予与实粒子相同的质量、密度、压力等参数。在计算过程中，耦合边界虚粒子参与到实粒子的守恒方程求解过程中，从而实现固壁边界条件。耦合边界粒子法的优点是可以降低SPH 方法的边界效应，但有时不足以防止实粒子非物理穿透边界。

本文针对边界排斥力法的不足，将近场动力学(Peridynamics，PD)中描述颗粒间接触作用的短程排斥力引入到固壁边界处理模型中。

PD方法是由美国Sandia国家实验室的Silling[9]提出的一种非局部、无网格物质点类方法[10-11]，已经在固体材料和结构的静动力变形以及破坏分析中成功应用[12-18]。在PD方法中，单颗粒视为由许多物质点组成，单颗粒内部物质点间与分属不同颗粒的物质点间均存在非局部作用力，该作用力形式简单，参数易于给定。边界排斥力法中边界虚粒子与实粒子之间的作用力，与颗粒间短程排斥力相似。本文主要工作是将短程排斥力模型引入到SPH方法的固壁边界处理中，并通过编程和具体算例的计算分析，证明边界排斥力法处理SPH固壁边界的有效性。

2 SPH方法的近似和控制方程的离散

2.1 SPH方法的近似

在SPH方法中，将连续体离散成有限个粒子，通过对某个粒子的支持域内其他粒子的加权求和获得该点的数值解。SPH方法近似分2步进行：第1步是积分近似，第2步是粒子近似。

函数f(x)在问题域Ω内的积分近似表示式为

(1)

式中：x，x′为包含在问题域Ω内点的坐标向量；dx′为x处无穷小体元；h为光滑长度；W为光滑函数(核函数)，本文采用三次样条函数[3]为光滑函数。

式(1)可写成离散化的粒子近似式，即

(2)

式中：mj，ρj分别表示粒子j的质量和密度，j=1,2,…,N；N为在x处粒子的支持域内的粒子总数。

粒子i处场函数的粒子近似式可写成

(3)

(4)

2.2 控制方程的SPH离散

Lagrangian描述下的流体动力学控制方程应用SPH近似后，离散化的SPH公式为：

(5)

(6)

(7)

式中：α表示Cartesian坐标分量x，y和z；在同一项中重复出现的下标(i或j)表示Einstein求和约定；ρ，v，p分别为密度、速度和压力；Πij为人工黏度项；f为体力；D/Dt代表物质导数。

式(5)和式(6)中采用对称结构，可以减少由于粒子不一致问题产生的误差[19]；另外，式(6)中Πij的作用是消除SPH方法在模拟流体动力学问题时产生的数值不稳定性，本文采用的人工黏度为文献[20]中的形式。式(6)中压力p是由弱可压缩流体的状态方程[6]得出，即

(8)

3 固壁边界施加模型

利用边界排斥力法来处理固壁边界，需要先在边界上布置一定数量的虚粒子，然后构造合理的排斥力公式来模拟真实边界力。考虑将PD方法中的颗粒间接触作用的短程排斥力引入到模型之中，作为边界虚粒子对实粒子的作用力，达到施加边界条件的目的。

PD理论认为物质点间的相互作用力具有非局部特征，这种作用力强化了近距离物质点间的相互作用，弱化了远距离物质点的相互作用，恰好吻合边界排斥力的分布特征。对于无黏性颗粒材料，PD方法采用短程排斥力描述分属不同颗粒的物质点间相互作用，本文将这种短程排斥力引入边界排斥力模型中，并采用文献[14]中弹性排斥力形式,即

(9)

式中：x，x′为物质点的位置矢量；d为物质点x和x′之间的短程排斥力作用的临界距离，表征当不同的物质点间距离

(10)

式中β为常量参数，取值在1.0左右。

若将布置在边界上的虚粒子和实粒子看作不同的物质点，d的大小等于光滑长度h的a倍，则边界虚粒子会对在其影响域内的实粒子沿着两粒子中心线方向上产生排斥力，该力的大小会随着两粒子的距离增大而减小，直到两粒子的距离超过d后作用力消失。

如图1所示，当实粒子B成为边界虚粒子A影响域内的粒子时，则会在沿着两粒子的中心线处虚粒子A对实粒子B产生一个作用力，即

(11)

图1 边界虚粒子对实粒子的排斥力Fig.1 The repulsive force of virtual particles on the boundary to real particles

传统SPH方法处理固壁的边界排斥力模型复杂且参数值的确定受人为因素影响较大，而上述改进的边界排斥力公式形式简单，参数值易于给定。

4 算例分析

本节中通过2个流体动力学的经典算例——Couette流和溃坝，来验证本文所研究的边界处理模型的有效性。

这里值得一提的是，对于初始粒子的数量和间距的设置，需要根据问题域的几何尺寸合理设定。过多的粒子虽然能保证精度但需要更多的计算量，过少的粒子又会造成计算精度降低。对于本文算例中初始粒子的设定，课题组均进行了不同数量粒子的数值试验；本文所采用SPH代码是由FORTRAN语言编写，在CPU为Intel Core i7、内存大小为32 G的台式电脑上运行。通过比较数值试验的CPU花费时间和计算精度，找到最佳的初始粒子配置。

4.1 Couette流

在Couette流模型中，初始静止的流体夹在2块固定且水平放置于间距为l的无限大平板中，流体由于上平板突然以恒定速度V0水平方向运动而产生流动，最终达到稳定状态。流体的尺寸是0.5 mm×1 mm，在SPH模拟中设置流体粒子的数量为20×40，初始间距为2.5×10-5m，流体密度ρ=1.0×103kg/m3，平板间距l=1.0×10-3m，V0=1.25×10-5m/s，见图2。

Couette流中流体水平速度Vx与时间t相关的理论级数表达式[21]为

(12)

式中：y为流体粒子的竖向位置；μ为运动黏度，这里取1.0×10-6m2/s，Reynolds数Re=1.25×10-2。在SPH模拟中，时间步长为1×10-5s，以流体最前端一列间隔选取的20个粒子的速度为参照，经过50 000步计算后流体的速度达到稳定状态。边界排斥力公式(式(11))中，h=2.5×10-5m，a=1.0，β=1.0。平板采用非滑移边界条件，并设定当流体沿平板切线方向运动时，边界虚粒子具有与流体实粒子大小相等方向相反的速度。

图3 SPH模拟与理论解的Couette流前端速度分布对比Fig.3 Comparison between the simulation solution in SPH and theory solution of front velocity distribution for Couette flow

图3给出SPH方法和理论级数解(式(12))得到的在t=0.01，0.1，0.2 s和最终稳定状态t=∞时流体前端速度分布对比，可以看出两者相当吻合，经计算得出SPH模拟值的最大相对误差为0.78%，表明本文研究的边界方法中的排斥力模型能较客观地反映真实边界力效果。

图4 溃坝的初始SPH粒子分布Fig.4 The initial particles distribution of dam breakin SPH simulation

图5 H=300 mm时试验[22]与SPH模拟在不同时刻的流场形态对比Fig.5 Comparison between experimental data[22] and SPH simulation result of flow field at different moments when H=300 mm

图6 H=600 mm时试验[22]与SPH模拟在不同时刻的流场形态对比Fig.6 Comparison between experimental data[22] and SPH simulation result of flow field at different moments when H=600 mm

4.2 溃 坝

图5和图6分别给出了H=300mm和600mm时，在选取的时间节点的试验和SPH模拟结果的流场形态比较。两图中SPH模拟结果都采用压强云图，从整体上看，试验和模拟结果对于自由表面运动的描述比较吻合，粒子压强分布也符合实际情况；即使在最容易发生实粒子非物理穿透边界的时刻[23-24]，如图5(f)和图6(f)，SPH模拟结果也没有出现粒子穿透现象，且粒子在整个过程中分布有序，说明本文所研究的边界力法能在不对流场产生明显干扰的情况下有效处理固壁边界。

对某些流场形态如水流垂直爬升高度、水流前端翻转的空腔的位置和大小，模拟结果与试验结果有些细节差异。其主要原因是SPH模拟中边界粒子对水粒子的影响较大，以及近似处理的流体粒子湍流问题与真实流体存在一定的差异；另外，SPH模拟中器壁简化为自由滑移边界，没有考虑气压等影响，这与试验环境难以保持完全一致，也是差异产生的可能原因。这些差异在文献[23-24]的SPH模拟溃坝算例中也同样存在。

图7 SPH模拟与试验[22]得到的溃坝水流前端位置比较Fig.7 Comparison between experimental data[22] and SPH simulation of propagation of the surge front position of dam-break water flow

表>1时水流前端归一化的平均速度的试验值与模拟值比较

注：相对误差=(试验值-模拟值)/试验值×100%

综合试验值和SPH模拟值的流场形态、水流前端位置和速度的比较，可以说明本文所研究的固壁边界处理方法在自由表面流问题中是有效的。

5 结 论

针对SPH固壁边界处理的边界排斥力法，本文提出一种新的排斥力模型。该模型的思路来源于PD方法中描述颗粒间接触作用的短程排斥力，提出的排斥力模型简单，参数易于设定，利于SPH固壁边界条件的施加。从Couette流和溃坝的数值算例结果来看，本文提出的排斥力模型能够很好地解决粒子非物理穿透边界，能较客观地反映真实的边界作用，粒子在排斥力的作用下运动规律，分布有序，表明该排斥力模型能够在对流场不产生明显干扰的情况下有效地处理固壁边界难题。

[1] LUCY L B. A Numerical Approach to the Testing of the Fission Hypothesis[J]. The Astronomical Journal，1977，82(12): 1013-1024.

[2] GINGOLD R A，MONAGHAN J J. Smoothed Particle Hydrodynamics：Theory and Application to Nonspherical Stars[J]. Monthly Notices of the Royal Astronomical Society，1977，181(3): 375-389.

[3] LIU M B，LIU G R. Smoothed Particle Hydrodynamics (SPH): An Overview and Recent Developments[J]. Archives of Computational Methods in Engineering，2010，17(1):25-76.

[4]MONAGHAN J J. Simulation Free Surface Flows with SPH[J]. Journal of Computational Physics， 1994，110(2): 399-406.

[5] LIBERSKY L D，PETSCHCK A G，CARNEY T C，etal. High strain Lagrangian Hydrodynamics: A Three-dimensional SPH Code for Dynamic Material Response[J]. Journal of Computational Physics，1993，109(1): 67-75.

[6] RANDLES P W，LIBERSKY L D. Smoothed Particle Hydrodynamics: Some Recent Improvements and Applications[J]. Computer Methods in Applied Mechanics and Engineering， 1996，139(1): 375-408.

[7] LIU M B, LIU G R. Smoothed Particle Hydrodynamics: A Meshfree Particle Method[M]. Singapore: World Scientific Publishing Co. Pte. Ltd., 2003: 138-141.

[8] LIU M B，LIU G R，LAM K Y. Investigations into Water Mitigations Using a Meshless Particle Method[J]. Shock Waves，2002，12(3):181-195.

[9] SILLING S A. Reformulation of Elasticity Theory for Discontinuities and Long-range Forces[J]. Journal of the Mechanics and Physics of Solids， 2000， 48(1): 175-209.

[10]SILLING S A， EPTON M， WECKNER O，etal. Peridynamic States and Constitutive Modeling[J]. Journal of Elasticity， 2007， 88(2): 151-184.

[11]黄 丹，章 青，乔丕忠，等. 近场动力学方法及其应用[J]. 力学进展，2010，40(4)：448-459. [12]KILIC B. Peridynamic Theory for Progressive Failure Prediction in Homogeneous and Heterogeneous Materials[D]. Tucson,USA: The University of Arizona， 2008.

[13]胡祎乐，余 音，汪 海. 基于近场动力学理论的层压板损伤分析方法[J]. 力学学报，2013，45(4): 624-628. [14]章 青， 顾 鑫， 郁杨天. 冲击荷载作用下颗粒材料动态力学响应的近场动力学模拟[J]. 力学学报，2016，48(1):56-63.

[15]WECKNER O，ABEYARATNE R. The Effect of Long-range Forces on the Dynamics of a Bar[J]. Journal of the Mechanics & Physics of Solids，2005，53(3): 705-728.

[16]SILLING S A，ASKARI E. A Meshfree Method Based on the Peridynamic Model of Solid Mechanics[J]. Computers & Structures，2005，83(17/18):1526-1535.[17]KILIC B. Peridynamic Theory for Progressive Failure Prediction in Homogeneous and Heterogeneous Materials[D]. Tucson,USA: The University of Arizona，2008.[18]HUANG Dan，LU Guang-da，QIAO Pi-zhong. An Improved Peridynamic Approach for Quasi-static Elastic Deformation and Brittle Fracture Analysis[J]. International Journal of Mechanical Sciences，2015，(94/95):111-122.

[19]MONAGHAN J J. An Introduction to SPH[J]. Computer Physics Communications，1988，48(1): 89-96.

[20]LATTANZIO J C，MONAGHAN J J，PONGRACIC H，etal. Controlling Penetration[J]. SIAM Journal on Scientific and Statistical Computing，1986，7(2): 591-598.

[21]MORRIS J P，PATRICK J F，ZHU Y. Modeling Low Reynolds Number Incompressible Flows Using SPH[J]. Journal of Computational Physics，1997，136(1): 214-226.

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

[24]LIU Hu, QIANG Hong-fu, CHEN Fu-zhen，etal. A New Boundary Treatment Method in Smoothed Particle Hydrodynamics[J]. Acta Physica Sinica, 2015, 64(9):094701.

(编辑：黄 玲)

A Repulsive Model for Solid Boundary Treatment inSmoothed Particle Hydrodynamics

ZHOU Xue-jun1，2， CHEN Ding2，TANG Yi3

(1.College of Mathematics and Physics，Huanggang Normal University，Huanggang 438000， China；2.College of Mechanics and Materials，Hohai University，Nanjing 210098，China；3.College of Mathematics and Computer Science，Yunnan University of Nationalities，Kunming 650500，China)

Boundary repulsive method is one of the methods for solid boundary treatment in smoothed particle hydrodynamics (SPH)， but the method is difficult to be widely applied due to the lack of unified repulsive model. The short-range repulsive force which describes the acting force between granules in Peridynamic (PD) is introduced to solid boundary treatment model to build a new boundary repulsive model in the framework of SPH. The reliability of the method is verified by two numerical simulation examples including Couette flow and dam-break. Moreover， the repulsive formulation is simple and the parameters are easy to be given. Therefore, the present method provides a new alternative for solid boundary treatment in SPH.

smoothed particle hydrodynamics (SPH); repulsive model; peridynamic (PD); solid boundary; Couette flow

2016-04-20；

2016-06-22

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

周学君(1981-)，男，湖北蕲春人，讲师，博士研究生，主要从事计算力学与工程仿真研究，(电话)13477625972(电子信箱)zhouxj@hhu.edu.cn。

10.11988/ckyyb.20160375

2017,34(7):54-59

O242

A

1001-5485(2017)07-0054-06