光滑粒子流体动力学的研究进展及应用

董添文，柳和生，黄益宾，余　忠

(上饶师范学院，江西 上饶 334001)

前言

1977年，Lucy[1]， Gingold与Monaghan[2]进行天体物理领域的计算时，提出了一种无网格粒子方法——光滑粒子流体动力学(Smoothed Particle Hydrodynamics，SPH)。当时有限元法非常流行，因此，SPH没有受到重视。近年来，计算力学界和工程领域界学者研究发现，SPH在计算一些复杂问题，例如：自由表面流、大变形、运动交界面等问题时，较传统有限元或有限差分法有很大优势。伴随着SPH方法的广泛应用，经典SPH的一些缺点也逐渐被发现，随后，学者纷纷提出了各种对传统SPH方法的改进措施，其理论体系也不断得到完善和发展。本文将围绕SPH方法的理论体系及其应用，主要包含计算流体力学领域、材料加工中的计算问题、弹塑性材料的变形与失效三个主要领域进行综述。

1　SPH方法的原理

1.1　核函数近似

如图1所示，SPH方法是纯基于Lagrange的粒子法，用大量粒子来表达计算域，它们携带密度，压力等信息。例如，位于r处的粒子i的场信息函数f(r)是通过它的支持域内的其它粒子的核函数插值来计算的，相当于有限元中的形函数插值。

(1)

式中，W(r-r′,h)=αdK(R),R=|r-r′|/h

〈〉——核近似算子

W(r-r′,h)——核函数

h——光滑域半径，是W紧支性的度量

利用高斯散度原理可以得到函数导数的核近似表达式

f(r)≅W(r-r′,h)dr′

(2)

图1　核函数近似示意图

1.2　离散原理

经典SPH方法采用所谓粒子近似法(属于一种配点法)来离散连续积分式(1)和式(2)。

(3)

f(ri)≅iWij

(4)

同样，对于函数空间二次导数的粒子近似式可以表示为

(5)

式中，△iWij——W对空间向量rij的二阶导数。

1.3　核函数

核函数一般是圆形或球形函数，应具有归一性、紧支性、非负性等属性。常用的核函数有：高斯核函数[2]、三次样条核函数[3]、分段四次样条核函数、分段五次样条核函数[4]等。通常，核函数越光滑(高阶可导)，越不会对粒子分布不规则敏感，故更稳定。在文献[5]中，Liu给出了一种核函数的构造方法。

1.4　Navier-Stokes方程的SPH表达式

Lagrange形式的Navier-Stokes方程如下：

连续方程

(6)

式中，ρ——密度；υ——速度矢量；t——时间

动量方程

(7)

式中，p——压力；τ——粘性剪切力张量；F——体积力

(8)

式中，μ——粘度系数；I——单位矩阵；ε——剪应变速率张量

(9)

能量方程

(10)

式中，q为热流密度。将以上公式右端进行SPH核近似和离散就可以得到N-S方程的SPH形式，SPH的离散格式很多[6,7]。

常用的连续方程离散格式有：

(11)

(12)

在边界附近，若采用式(11)，计算结果会出错，因此在自由表面流的计算中常用式(12)。

压力项的对称离散格式为

(13)

粘性力项的对称离散格式为

对于不可压缩粘性流，粘性力计算公式为

(14)

最早，Takeda[8]采用直接对核函数求二次导数的方法，即式(5)的方法离散式(14)，然而，Brookshaw[9]研究发现二阶导数直接近似公式对粒子不规则性敏感，会引起计算不稳定，因此又出现了两种方法：一种是Flebbe等[10]提出来的嵌套核函数插值法，将二次导数分别用两次一阶导数核函数插值的嵌套乘积方式完成；另一种是，将其中一次求导用差分计算，再进行核函数导数插值的表达式[11,12]。研究表明该方法在计算粘度项时更加稳定，且运算量小，因此在实际应用中使用最多[12-16]，其离散表达式为

(15)

利用以上离散法则，可以得到动量控制方程和能量控制方程的SPH形式：

(16)

(17)

1.5　状态方程

因“压力-速度”耦合，不可压缩流的计算较可压缩流困难更大。经典SPH方法中，通过引入“人工压缩率”的概念(理论上不可压缩流体实际上是可压缩的)，建立状态方程求解压力[17,18]，因此，经典SPH也被称作弱可压缩SPH，(Weakly compressible SPH ，WCSPH)。常用的状态方程有：

p=ρc2

(18)

(19)

式中：γ=7,B=c2ρ0/γ，ρ0是参考密度，c是声速。通常，式(18)在计算非自由表面流时使用，在计算自由

表面流时使用式(19)。一般声速不能取真实物理值，这是因为，真实声速非常大，时间步长则将受声速限制而太小，会影响计算效率，但当声速过小，则会出现密度波动大，进一步引起压力振荡，造成不稳定，Monaghan建议，c至少要为最大粒子速度的10倍，在Lee[14]的研究建议，c=100υmax，具体问题中，还需要调试。

2　经典SPH存在的问题

2.1　粒子插值的不连续

连续性意味着插值形函数再生多项式的能力[19]，且与收敛性，稳定性密切相关。Bonet[20]，Fang[21]等指出：角动量或线动量的不严格守恒会导致长时间模拟的不稳定。线动量守恒只要求粒子插值具有0阶连续性，而角动量守恒，则要求粒子插值近似至少具有一阶连续性。Liu[22]给出函数及其前两阶导数具有n阶核函数近似精度必须满足的条件。0阶连续性条件即是核函数归一化条件，1阶连续性则要求

(20)

然而，它们是由积分推导出的，不能保证离散化后方程的连续性。如图2所示：当粒子位于边界附近时，它的支持域就被边界截断而不完整，这时是不能保证离散后的方程满足插值连续的。另外，当粒子的分布极不规则时，也不能保证插值连续性[23]。

图2　核函数被边界截断

2.2　张力不稳定

早期学者发现，即使满足CFL条件，当粒子之间的应力为负(即为张力)的时候，粒子会相互吸引，形成大量粒子“团聚”，最终造成计算失败。最近研究表明，这种现象也可能在应力都是正的情况中出现，在弹性体的模拟中则表现为非物理的断裂[24]。Swegle等[25,26]采用Von Neumann判据对张力不稳定进行了深入研究。发现造成张力不稳定的根本原因是由于核插值改变了粒子之间的本构关系，从而，改变了原始偏微分方程的性质。并提出了张力不稳定性判据：

△WT>0

(21)

很多学者提出了解决张力不稳定性问题的方法。Morris[4]研究发现高阶样条核函数与高斯型核函数能显著降低不稳定性的增长速度。然而，高阶核函数在紧支域的某些地方为负值，将产生非物理结果[7]。Johnson[27]提出二次光滑函数来模拟高速冲击问题，并称可解决压力不稳定，但是因△W恒大于零，故不能解决拉伸失稳。

Swegle等[28]等提出守恒光滑法，通过平滑2h波长的振动替代常用的人工粘性力来消除张力不稳定性。Vignjevic[29]称CSA法与施加人工排斥力法会影响材料强度，故使用时要引起注意。Dyka[30,31]等提出了“应力点法”的思想是在SPH粒子外设置应力点粒子来求应力和其它变量信息，而SPH粒子自身的运动学变量信息例如，位移、速度、质量等不变。该应力点方法和EFG方法中的“高斯积分点”类似[32]。Vignjevic[29]利用正则化核函数将“应力点法”扩展到二维问题。Monaghan[24]提出一种人工排斥力的方法来消除张力不稳定性。

Randles与Libersky[33]提出“双粒子动力学”方法来解决张力不稳定性。另外，Liu[34]提出的RKPM，Chen等[35-38]提出了CSPM等方法，能改善甚至完全消除张力不稳定性。Hu[39,40]在他的计算中采用了“参考压力法”来消除张力不稳定性。

2.3　边界条件实施困难

因核函数被边界截断， SPH在边界条件的处理上较其它数值方法更困难。近期主要使用的边界施加方法有：

1) 排斥力方法

排斥力方法理论来源于计算分子力的Lennard-Jones方程[17]。该方法是在固壁边界设置一组虚粒子，对靠近边界的内部粒子施加一定的排斥力，防止它穿透边界。后来研究发现，早期的模型会导致粒子沿边界“游走”，与实际情况不符。为此，Monaghan[41-43]对排斥力模型进行了改进。该方法的优点是不受边界形状的影响，实施方便，缺点是模型的参数依然凭经验设置。

2) 镜像粒子法

镜像粒子方法是在每个时间步，距离边界附近一定范围内，将内部粒子以边界为对称面生成相应的镜像虚粒子[44]。Takeda[8]，Morris[18]分别对镜像粒子方法进行改进用来处理曲边界。虚粒子边界的守恒性好，但对于复杂形状边界(比如锯齿状边界)或两相交界面，这种方法实施起来有困难[43]。

3) 静态粒子边界

该方法与镜像粒子法不同，其中虚粒子是固定的，虚粒子要参与内部粒子的运算，虚粒子层可看作固壁的加厚[14,15,45]。Dalrymple与Knio[46]提出了被称作“Dynamic Boundary"的一种双层交错粒子来模拟边界。与Koshizuka的“DUMMY”粒子相似，边界粒子与内部粒子一样要参与N-S方程计算，但是，位置固定不变或另外施加(例如：造浪器)。该方法被证明守恒性好，边界实施容易[47]。

3　对经典SPH方法的改进

3.1　NSPH( Normalized SPH)

Randles等[48]在Johnson等[49]的工作的基础上，基于归一化的思想提出了对核函数及导数的修正方法

(22)

(23)

式(22)与(23)分别保证了核函数插值的0阶与1阶连续。Bonet等[20]利用变分原理推导得到相同的核函数导数修正公式，保证角动量守恒。Belytschko[19]认为，该方法虽然能保证修正的核函数插值线性连续性，但修正的核函数导数不具有可积性，会导致采用Galerkin法建立的离散格式不能通过分片试验。

3.2　RKPM(Reproducing kernel method，RKPM)

RKPM是在SPH方法的基础上发展起来的一种基于核近似的无网格方法。由于SPH存在两个不足：粒子插值难满足线性一致性条件；当总粒子数目相对较少时，计算精度降低。Liu等针对这些问题对SPH方法进行改进，提出RKPM[34]。RKPM利用修正的核函数来保证粒子插值连续性，构造修正的核函数：

(24)

由n阶再生条件，建立离散线性方程组，解方程组可以得到修正的核函数。因为高阶再生核函数计算量大，在实际应用中线性再生核用得较广泛。RKPM采用的是全局Galerkin离散格式。RKPM方法不仅解决了SPH法在边界上的不连续性，而且完全消除了张力不稳定性。Liu还将它与小波分析技术结合，使得RKPM具有多尺度自适应分析能力。

RKPM法已经有大量的应用研究：如结构动力学[50]、应力集中[51]、大变形[52]、剪切带的形成[53]等。

3.3　CSPM(Corrective smoothed particle method)

CSPM方法是Chen等[35-38]提出的。以二维直角笛卡儿坐标系为例，将函数f在点(xi,yi)处展开，并在等式两边同时乘以核函数W，并在求解域Ω内积分得

忽略所有求导项，则得到核近似和粒子求和近似公式

(26)

明显，式(26)与式(22)是等价的。通常，若要得到函数的一阶导数近似，只需要式(25)中的核函数具有反对称性(通常采用核函数一阶导数W,x)即可消除0次导数项与二次导数项，高阶导数求解依此类推。该方法中使用的核函数不要求与经典SPH方法的核函数具有相同的性质，例如：经典光滑核函数要求非负，在CSPM方法中则不要求，此方法的核函数选择更灵活，但要求保证联立方程组的系数矩阵不会奇异。CSPM在求解域内具有二阶精度，在边界附近也具有一阶精度。该方法求解函数及其导数的公式与NSPH方法相同，但这两种方法的思想是不同的，前者更接近数学，后者更接近物理。CSPM方法具有三个优点：理论上它可以方便的求解任意阶导数；它在边界条件的实施上较经典SPH方法简单，边界粒子的参数可以直接设置；它在将非线性偏微分方程离散化时，较经典SPH方法更直接。Zhang等[54]在推导函数及其导数的计算公式时，均保留高阶导数项，建立联立方程组，统一求解。这种方法精度较高，但是计算量更大。Fang[21]为了平衡精度与效率，将CSPM方法与经典SPH方法耦合，在边界附近采用CSPM方法求解以提高精度与稳定性，在计算域内部采用SPH方法来提高计算效率。

3.4　MLSPH(Moving least squares particle hydrodynamics)

MLS(Moving least squares)近似最早是由数学家设计用于数据拟合以及表面构造[55]，利用MLS构造的形函数可积，适当选择权函数可以保证MLS形函数满足全局连续性[19]。MLS已经被广泛用来构造无网格法的形函数，如EFG，FPM等。Dilts[56,57]提出移动最小二乘SPH方法MLSPH。

(rj)=PT·A-1·p(rj)Wij

(27)

(28)

式中，p(rj)是基函数向量在rj处的值，以二次多项式基为例

P(r)=[1,x,y,x2,xy,y2]T,m=6

(29)

pT是计算点ri的邻域内各点的基函数向量组成的矩阵。

MLS形函数可以看作对经典SPH形函数，在粒子分布不规则时的一种修正，当粒子规则分布时有

(30)

MLSPH方法能提高精度和消除张力不稳定性，然而，Monaghan[24]指出它的计算量大约是经典SPH方法的8倍。RKPM，NSPH与MLSPH均是对经典对称非负核函数及其导数进行修正或重构来提高精度，而重构

的核函数在紧支域中可能出现负值，这样就产生了一些场变量出现非物理值，比如：负密度，负能量，这样会造成计算崩溃。而CSPM不改变经典SPH的核函数，不会出现核函数出现负值情况。

3.5　ISPH(Incompressible SPH)

Cummins[58]综合移动粒子半隐式方法(MPS)与投影方法(Projection)的优点，提出了一种投影SPH方法。为了区别于WCSPH，该类方法后来被称作“ISPH"。该方法属于分数步法，在中间临时步，可以先不考虑不可压缩条件，即先假设粒子不受压力而只受粘性力或外力，当粒子运动后，它的密度或速度散度会发生变化，再根据不可压缩条件可以推导出压力Possion方程

(31)

求解方程得到压力值，最后，再考虑压力的影响，相当于把速度和位置都改变的粒子再“压”回去。与Projection方法采用有限差分离散不同，ISPH函数的离散是通过光滑核函数插值离散完成的。近年来的研究表明ISPH在收敛性，稳定性，计算精度，计算效率，压力值求解诸方面均较经典SPH方法好[14-16]。近期有不同学者对ISPH和经典SPH的应用情况进行了比较[59-62]。发现，如果对经典WCSPH算法计算的密度每隔20步进行MLS修正，计算结果与ISPH相当。

4　SPH的应用

4.1　可压缩流体计算

SPH方法最早是由天文学家在计算天体物理中提出来的，因为天体运动可以被看作是可压缩流[1,2]。Springel[63]基于SPH方法开发了一套用于计算天体的软件，开展了大量计算，例如星系碰撞、宇宙结构的形成、千年计算等。SPH在可压缩流计算领域的另一个重要应用是对爆炸的模拟。因为在爆炸过程中，存在着极大变形和自由表面，传统有网格方法难以胜任。特别是对一些细节，例如，破碎、飞溅、融合等的模拟，无网格粒子法表现更佳。Liu[64]等将SPH方法应用到水下爆炸冲击问题中。

4.2　不可压缩流体计算

1994年，Monaghan首次将SPH方法应用到自由表面流的模拟，也为SPH方法在不可压缩流中的计算开了先河[17]；随后,世界各地学者纷纷加入到这个领域中来，成立一个SPH欧洲研究兴趣组织(SPH European Research Interest Community, SPHERIC)。SPHERIC的成员开展了大量研究，其中应用最广泛的领域是在海洋工程，例如波浪对海岸(或石油钻井平台)的冲击、船舱中的液体晃荡。在工业领域应用较为成功的是铸造充型过程的模拟。Cleary等[65]采用SPH方法对高压压铸充型过程进行了研究，与VOF方法计算结果与实验进行比较，表明SPH可以有效地计算这种高速流动现象。国内，曹文炅等[66]采用类似方法对高压压铸充型过程进行了分析。SPH方法在工业领域复杂加工中的应用有很多，例如，采用ISPH算法模拟聚合物螺杆挤出展开流道[67]、模拟磨料水射流中单磨粒加速过程[68]、激光水下加工中，水与熔融金属之间的相互作用[69]。长期阻碍其应用发展的两大主要因素有：复杂加工领域中的计算域几何形态复杂，且都是动态变化的；SPH作为一种无网格粒子法，其计算量远远大于传统网格方法，特别是在粒子搜索阶段，每个时间步，都需要重新进行邻域粒子搜索配对，而有网格法，一旦网格固定，节点与节点的连接是固定的。最近，SPHERIC成员提出了一种前处理器，首先利用已经非常成熟的商业网格法生成器进行复杂计算域的网格划分，然后再将网格区域用粒子代替；他们采用了基于GPU的多核并行计算方案，计算效率比单CPU版的程序提高了60倍[70]。

4.3　材料的变形与冲击断裂

材料加工过程中，存在材料的大变形。Libersky[44,48,71]等将SPH应用于锻压过程的计算。近年来研究者开始将SPH法应用到切削过程的模拟中．Heinstein等[72]证明了SPH法在直角切削模拟中的有效应用．Limido等[73]采用SPH法对高速切削过程进行模拟，分析了切屑的形成并进行了切削力预报。Spreng[74]采用采用自适应SPH算法模拟切屑过程，避免了数值断裂。我国学者也开展了很多研究。例如，郭晓光[75]等将SPH方法应用于超精密切削过程仿真。

Johnson[27, 48,71]等最早将它应用到高速碰撞问题的研究中。徐金中等[76]对初始光滑长度和粒子间距对计算结果的影响进行了深入研究。张志春等[77]采用基于SPH-FEM耦合的算法对子弹冲击高强钢板进行了数值计算。

5　展望

(1) 工程领域的应用基本都是三维问题，有时甚至是多相情况下的三维问题，边界实施较二维更为复杂。采用排斥力设置边界的优点是对复杂形态边界的适应性强，缺点是参数的设置一直没有一套科学的，简单的，普适性的方案。值得进一步深入探索。

(2) 对于无网格粒子法，计算时间会随着粒子数量增加急剧增加，这是无网格粒子法在大型复杂加工过程中的应用的一个瓶颈。除了依靠并行算法以外，是否能有从算法本身进行新的改进措施，特别是对于每个时间步都要重新进行的粒子邻域搜索，是计算耗时一个关键因素，有待深入探索。

(3) SPH方程在不同工业领域的应用中不能完全通用，需要结合实际情况，进行修正，例如，在聚合物领域中的应用，要考虑流体粘弹性问题，在金属切屑加工中的应用要考虑金属剥离后的粒子邻域搜索。

参考文献：

[1] Lucy L. B. A numerical approach to the testing of the fission hypothesis [J] . Astrophys. J. 1977,82, 1013-1024.

[2] Gingold R.A. and Monaghan JJ, Smoothed particle hydrodynamics: theory and application to non-spherical stars [J]. Mon. Not. R.Astron. Soc. 1977, 181, 375-389.

[3] Monaghan J. J. Particle methods for hydrodynamics[J]. Computer Physics Report, 1985, 3(2): 71-124.

[4] Morris J. P.A study of the stability properties of smooth particle hydrodynamics[J]. Publ.Astron.Soc.Aust, 1996, 13(1): 97-102.

[5] Liu G. R. and Liu M. B. 光滑粒子流体动力学——一种无网格粒子法[M]. 韩旭, 杨刚, 强洪夫, 译. 长沙: 湖南大学出版社, 2005.

[6] Monaghan J. J. Smoothed particle hydrodynamics[J]. Annual Review Astronomics and A strophysics, 1992, 30: 543-574.

[7] Fulk D. A. A numerical analysis of smoothed particle hydrodynamics[D]. Air Force Institute of Technology, 1994.

[8] Takeda H., Miyama S. M., Sekiya M. Numerical simulation of viscous flow by smoothed particle hydrodynamics[J]. Progress of Theoretical Physics, 1994, 92(5): 939-960.

[9] Brookshaw L. A method of calculating radiative heat diffusion in particle simulations[C]. Proceedings of the Astronomical Society of Australia, 1985, 6(2): 207-210.

[10] Flebbe O., Munzel S., Herold H. et al. Smoothed particle hydrodynamics- physical viscosity and the simulation of accretion disks[J]. The Astrophysical Journal, 1994, 431(2): 754-760.

[11] Dilts. G. A. Some recent developments for moving least squares particle methods[C]. Proceedings on Computational Fluid and Solid Mechanics, Cambridge MA, 2001: 1-5.

[12] Shao S. D., Edmond Y. M. Lo. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[J]. Advances in Water Resources, 2003, 26(7): 787-800.

[13] Cummins S. J. Rudman. M. An SPH projection method[J]. Journal of Computational Physics. 1999, 152(2): 584-607.

[14] Lee E. S., Moulinec C., Xu,R. et al. Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method[J]. Journal of Computational Physics, 2008, 227(18): 8417-8436.

[15] Khayyer A., Gotoh H., Shao. S. D. Enhanced predictions of wave impact pressure by improved incompressible SPH methods[J]. Applied Ocean Research, 2009, 31(2): 111-131.

[16] Xu R., Stansby P., Laurence. D. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[J]. Journal of Computational Physics, 2009, 228(18): 6703-6725.

[17] Monaghan J. J. Simulating Free Surface Flows with SPH[J]. Journal of Computational Physics, 1994, 110(2): 399-406.

[18] Morris J. P., Fox P. J. Y. Zhu. Modeling low reynolds number incompressible flows using SPH[J]. Journal of Computational Physics, 1997, 136(1): 214-226.

[19] Belytschko T., Krongauz Y., Dolbow J. etal. On the completeness of the meshfree particle[J]. International Journal for Numerical

Methods in Engineering, 1998, 43(5): 785-819.

[20] Bonet J., Lok. T-S L. Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations[J]. Comput. Methods Appl. Mech. Engrg, 1999, 180(1-2): 97-115.

[21] Fang J. N. Parriaux, A. Rentschler,M. et al. Improved SPH methods for simulating free surface flows of viscous fluids[J]. Applied Numerical Mathematics, 2009, 59(2): 251-271.

[22] Liu M. B., Liu, G. R. Lam. K. Y. Constructing smoothing functions in smoothed particle hydrodynamics with applications[J]. Journal of Computational and Applied Mathematics, 2003, 155(2): 263-284.

[23] Belytschko T., Krongauz Y., Organ, D. et al. Meshless methods: An overview and recent developments[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4): 3-47.

[24] Monaghan. J. J. SPH without a tensile instability[J]. Journal of Computational Physics, 2000, 159(2): 290-311.

[25] Balsara D. S. Von Neumann stability analysis of smoothed particle hydrodynamics suggestions for optimal algorithms[J]. Journal of Computational Physics, 1995, 121(2): 357-372.

[26] Swegle J. W., Hicks, D. L. and Attaway. S. W. Smoothed Particle Hydrodynamics stability analysis[J]. Journal of Computational Physics, 1995, 116(1): 123-134.

[27] Johnson G. R., Petersen, E. H. and Stryk. R. A. Incorporation of an SPH option into EPIC code for a wide range of high velocity impact computations[J]. Int. J. Impact Eng, 1993, 14 (1-4): 385-394.

[28] Guenther C., Hicks D. L., Swegle. J .W. Conservative smoothing versus artificial viscosity[R]. SANDIA REPORT, SAND94-1853 UC-705, 1994.

[29] Vignjevic. R. Review of development of the SPH method predictive modeling of dynamic processes[M]. US: Springer, 2009, 367-396.

[30] Dyka C. T., Inge. R. P. Addressing tension instability in SPH Methods[R]. Technical Report NRL/MR/6384, NRL, 1994.

[31] Dyka C. T., Inge. R. P. An approach for tension instability in smoothed particle hydrodynamics (SPH) [J]. Computers & Structures, 1995, 57(4): 573-580.

[32] Li S. F., Liu. W. K. Meshfree and particle methods and their applications[J]. Appl Mech Rev, 2002, 55(1): 1-34.

[33] Randles P. W., Libersky.L. D. Normalized SPH with stress points[J]. International Journal for Numerical Methods in Engineering, 2000, 48(10): 1445-1462.

[34] Liu W. K., Jun S., Zhang. Y. F. Reproducing kernel particle methods[J]. Int. J. Numer. Meth. Fluids. 1995, 20(8-9): 1081-1106.

[35] Chen J. K., Beraun.J. E. A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems[J]. Comput. Methods Appl. Mech. Engrg, 2000, 190(1-2): 225-239.

[36] Chen J. K., Beraun,J. E. and Jih. C. J. An improvement for tensile instability in smoothed particle hydrodynamics[J]. Computational Mech, 1999, 23(4): 279-287.

[37] Chen J. K., Beraun,J. E. and Jih. C. J. A corrective smoothed particle method for boundary value problems in heat conduction[J]. International Journal for Numerical Methods in Engineering, 1999, 46(2): 231-252.

[38] Chen J. K., Beraun J. E., and Jih. C. J. Completeness of corrective smoothed particle method for linear elastodynamics[J]. Computational Mechanics, 1999, 24(4): 273-285.

[39] Hu X. Y., Adams. N. A. A constant-density approach for incompressible multi-phase SPH[J]. Journal of Computational Physics, 2009, 228(6): 2082-2091.

[40] Hu X. Y., Adams. N. A. An incompressible multi-phase SPH method[J]. Journal of Computational Physics, 2007, 227(1): 264-278.

[41] Monaghan J. J., Kos. A. Solitary waves on a cretan beach[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1999, 125(3): 145-154.

[42] Monaghan J. J., Kos A. and Issa. N. Fluid motion generated by impact[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2003, 129(6): 250-260.

[43] Monaghan J. J., Kajtar. J. B. SPH particle boundary forces for arbitrary boundaries[J]. Computer Physics Communications, 2009, 180(10): 1811-1820.

[44] Libersky L. D., Petschek A. G., Carney T. C., et al. High strain Lagrangian hydrodynamics a three-dimensional SPH code for dynamic material response[J]. J. Comput. Phys, 1993, 109(1): 67-75.

[45] Koshizuka S, Nobe A, Oka Y. Numerical analysis of breaking waves using the moving particle semi-implicit method[J]. International Journal for Numerical Methods in Fluids, 1998, 26(7): 751-769.

[46] Dalrymple, R. and Knio, O. SPH modelling of waves[C]. In ASCE, editor, Proc. Coastal Dynamics,Lund, Sweden, 2001: 779-787.

[47] Crespo A. J. C. Application of the smoothed particle hydrodynamics model SPHysics to free-surface hydrodynamics[D]. UNIVERSIDADE DE VIGO DEPARTAMENTO DE FISICA APLICADA. 2008.

[48] Randles P. W. and Libersky. L. D. Smoothed particle hydrodynamics: Some recent improvements and applications[J]. Comput. Methods Appl. Mech. Eng, 1996, 139(1-4): 375-408.

[49] Johnson G. R. and Beissel. S. R. Normalized smoothing functions for SPH impact computations[J]. Int. J. Numer. Methods Eng, 1996, 39(16): 2725-2741.

[50] Liu W. K., Jun S., Li S. F.et al. Reproducing kernel particle methods for structural dynamics[J]. Int. J. Numer. Methods Engrg. 1995, 38(10): 1655-1679.

[51] Lee S. H., Kim H. J.,Jun. S. Two scale meshfree method for the adaptivity of 3D stress concentration problems[J]. Comput. Mech. 2000, 26(4): 376-387.

[52] Chen J. S., Pan C., Wu C. T. et al. Reproducing kernel particle methods for large deformation analysis of non-linear structures[J]. Comput. Methods Appl. Mech. Engrg. 1996, 139(1-4): 195-227.

[53] Li S. F., Hao W., Liu. W. K. Meshfree simulations of shear banding in large deformation[J]. International Journal of Solid and Structures, 2000, 37(48-50): 7185-7206.

[54] Zhang G. M., Batra. R. C.Modified smoothed particle hydrodynamics method and its application to transient problems[J]. Computational Mechanics, 2004, 34(2): 137-146.

[55] Lancaster P. and Salkauskas. K. Surface generated by moving least squares methods[J]. Mathe Comp, 1981, 37(155): 141-158.

[56] Dilts G. A. Moving least-square particle hydrodynamics I: Consistency and stability[J]. Int. J. Numer. Methods Eng, 1999, 44(8): 1115-1155.

[57] Di B G. A. Dilts. Moving least-square particle hydrodynamics II: Conservation and boundaries[J]. Int. J. Numer. Methods Eng, 2000, 48(10): 1503-1524.

[58] Cummins S. J., Rudman. M. An SPH projection method[J]. Journal of Computational Physics. 1999, 152(2): 584-607.

[59] Hughes J. P., Graham. D. I. Comparison of incompressible and weakly-compressible SPH models for free-surface water flows[J]. Journal of Hydraulic Research. 2010, 48, s(1):105-117.

[60] Szewc K., Pozorski J. and Minier. J. P. Analysis of the incompressibility constraint in the smoothed particle hydrodynamics method[J]. Int. J. Numer. Meth. Engng. 2012, 92(4): 343-369.

[61] Rafiee A., Cummins S., Rudman M. etal. Comparative study on the accuracy and stability of SPH schemes in simulating energetic free-surface flows[J]. European Journal of Mechanics B/Fluids. 2012, 36: 1-16.

[62] Shadloo M. S., Zainali A., Yildiz M. et al. A robust weakly compressible SPH method and its comparison with an incompressible SPH[J]. Int. J. Numer. Meth. Engng. 2012, 89(8): 939-956.

[63] Springel V. The cosmological simulation code GADGET-2[J]. Mon. Not. R. Astron. Soc. 2005,364, 1105-1134.

[64] Liu M. B., Liu G. R., Lam K. Y. Investigations into water mitigation using a meshless particle method [J]. Shock Waves, 2002, 12(3): 181-195.

[65] Cleary P.,Prakash M., Ha. J. Novel applications of smoothed particle hydrodynamics (SPH) in metal forming[J]. Journal of Materials Processing Technology, 2006, 177(3): 41-48.

[66] 曹文炅, 周照耀, 何毅等. 压铸充型过程的SPH方法建模及数值模拟[J]. 华南理工大学学报(自然科学版), 2011, 39(1): 100-105.

[67] 董添文, 江顺亮, 黄兴元, 柳和生, 基于不可压缩光滑粒子流体动力学的单螺杆挤出三维流动数值模拟[J], 机械工程学报, 2012, 48(22): 80-86.

[68] 王建明, 余丰, 刘飞宏. SPH和FEM 耦合法模拟磨料水射流中单磨粒加速过程[J]. 山东大学学报(工学版). 2011, 41 (5): 114-120.

[69] Yan Y. Z. Li L. Sezer K. etal. CO2 laser underwater machining of deep cavities in alumina[J]. Journal of the European Ceramic Society. 2011, 31(15): 2793-2807.

[70] Valdez-Balderas D., Dominguez J., Rogers B.D., Crespo A.J.C. Towards accelerating Smoothed Particle Hydrodynamics simulations for free-surface flows on multi-GPU clusters[J]. Computer Physics Communications, 2013,184(3): 617-627.

[71] Libersky LD and Petschek AG, Smooth particle hydrodynamics with strength of materials, In: Advances in the Free-Lagrange Method[J]. Springer, New York, 1991:248-257.

[72] Heinstein M, Segalman D. Simulation of or thogonal cutting with smooth particle hydrodynamics[R]. Sandia National Laboratories Report. New Mexico: Sandia National Laboratories, 1997.

[73] Limidoa J. Espinosaa C. Salaüna M., Lacome J.L. SPH method applied to high speed cutting modeling[J]. 2007, 49(7): 898-908.

[74] Fabian S. Dirk S. Alexandra M. etal. Smoothed Particle Hydrodynamics with Adaptive Discretization[C]. 9th international SPHERIC workshop Paris, France, June, 2014.

[75] 郭晓光, 魏延军, 张小冀 等. 基于FEM-SPH耦合方法的AISI4340钢超精密切削过程仿真[J]. 大连理工大学学报, 2013, 53(4): 526-531.

[76] 徐金中, 汤文辉. SPH方法在层裂损伤模拟中的应用[J].强度与环境, 2009, 36(1):1-7.

[77] 张志春, 强洪夫, 傅学金等. 基于SPH-FEM 转换算法的7. 62 mm步枪弹冲击30CrMnSiA 钢板的数值计算[J]. 计算物理，2012, 29(1):73-81.