董添文，柳和生，黃益賓，余　忠

(上饒師范學院，江西 上饒 334001)

前言

1977年，Lucy[1]， Gingold與Monaghan[2]進行天體物理領域的計算時，提出了一種無網(wǎng)格粒子方法——光滑粒子流體動力學(Smoothed Particle Hydrodynamics，SPH)。當時有限元法非常流行，因此，SPH沒有受到重視。近年來，計算力學界和工程領域界學者研究發(fā)現(xiàn)，SPH在計算一些復雜問題，例如：自由表面流、大變形、運動交界面等問題時，較傳統(tǒng)有限元或有限差分法有很大優(yōu)勢。伴隨著SPH方法的廣泛應用，經(jīng)典SPH的一些缺點也逐漸被發(fā)現(xiàn)，隨后，學者紛紛提出了各種對傳統(tǒng)SPH方法的改進措施，其理論體系也不斷得到完善和發(fā)展。本文將圍繞SPH方法的理論體系及其應用，主要包含計算流體力學領域、材料加工中的計算問題、彈塑性材料的變形與失效三個主要領域進行綜述。

1　SPH方法的原理

1.1　核函數(shù)近似

如圖1所示，SPH方法是純基于Lagrange的粒子法，用大量粒子來表達計算域，它們攜帶密度，壓力等信息。例如，位于r處的粒子i的場信息函數(shù)f(r)是通過它的支持域內(nèi)的其它粒子的核函數(shù)插值來計算的，相當于有限元中的形函數(shù)插值。

(1)

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

〈〉——核近似算子

W(r-r′,h)——核函數(shù)

h——光滑域半徑，是W緊支性的度量

利用高斯散度原理可以得到函數(shù)導數(shù)的核近似表達式

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

(2)

圖1　核函數(shù)近似示意圖

1.2　離散原理

經(jīng)典SPH方法采用所謂粒子近似法(屬于一種配點法)來離散連續(xù)積分式(1)和式(2)。

(3)

f(ri)?iWij

(4)

同樣，對于函數(shù)空間二次導數(shù)的粒子近似式可以表示為

(5)

式中，△iWij——W對空間向量rij的二階導數(shù)。

1.3　核函數(shù)

核函數(shù)一般是圓形或球形函數(shù)，應具有歸一性、緊支性、非負性等屬性。常用的核函數(shù)有：高斯核函數(shù)[2]、三次樣條核函數(shù)[3]、分段四次樣條核函數(shù)、分段五次樣條核函數(shù)[4]等。通常，核函數(shù)越光滑(高階可導)，越不會對粒子分布不規(guī)則敏感，故更穩(wěn)定。在文獻[5]中，Liu給出了一種核函數(shù)的構造方法。

1.4　Navier-Stokes方程的SPH表達式

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

連續(xù)方程

(6)

式中，ρ——密度；υ——速度矢量；t——時間

動量方程

(7)

式中，p——壓力；τ——粘性剪切力張量；F——體積力

(8)

式中，μ——粘度系數(shù)；I——單位矩陣；ε——剪應變速率張量

(9)

能量方程

(10)

式中，q為熱流密度。將以上公式右端進行SPH核近似和離散就可以得到N-S方程的SPH形式，SPH的離散格式很多[6,7]。

常用的連續(xù)方程離散格式有：

(11)

(12)

在邊界附近，若采用式(11)，計算結果會出錯，因此在自由表面流的計算中常用式(12)。

壓力項的對稱離散格式為

(13)

粘性力項的對稱離散格式為

對于不可壓縮粘性流，粘性力計算公式為

(14)

最早，Takeda[8]采用直接對核函數(shù)求二次導數(shù)的方法，即式(5)的方法離散式(14)，然而，Brookshaw[9]研究發(fā)現(xiàn)二階導數(shù)直接近似公式對粒子不規(guī)則性敏感，會引起計算不穩(wěn)定，因此又出現(xiàn)了兩種方法：一種是Flebbe等[10]提出來的嵌套核函數(shù)插值法，將二次導數(shù)分別用兩次一階導數(shù)核函數(shù)插值的嵌套乘積方式完成；另一種是，將其中一次求導用差分計算，再進行核函數(shù)導數(shù)插值的表達式[11,12]。研究表明該方法在計算粘度項時更加穩(wěn)定，且運算量小，因此在實際應用中使用最多[12-16]，其離散表達式為

(15)

利用以上離散法則，可以得到動量控制方程和能量控制方程的SPH形式：

(16)

(17)

1.5　狀態(tài)方程

因“壓力-速度”耦合，不可壓縮流的計算較可壓縮流困難更大。經(jīng)典SPH方法中，通過引入“人工壓縮率”的概念(理論上不可壓縮流體實際上是可壓縮的)，建立狀態(tài)方程求解壓力[17,18]，因此，經(jīng)典SPH也被稱作弱可壓縮SPH，(Weakly compressible SPH ，WCSPH)。常用的狀態(tài)方程有：

p=ρc2

(18)

(19)

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

表面流時使用式(19)。一般聲速不能取真實物理值，這是因為，真實聲速非常大，時間步長則將受聲速限制而太小，會影響計算效率，但當聲速過小，則會出現(xiàn)密度波動大，進一步引起壓力振蕩，造成不穩(wěn)定，Monaghan建議，c至少要為最大粒子速度的10倍，在Lee[14]的研究建議，c=100υmax，具體問題中，還需要調(diào)試。

2　經(jīng)典SPH存在的問題

2.1　粒子插值的不連續(xù)

連續(xù)性意味著插值形函數(shù)再生多項式的能力[19]，且與收斂性，穩(wěn)定性密切相關。Bonet[20]，F(xiàn)ang[21]等指出：角動量或線動量的不嚴格守恒會導致長時間模擬的不穩(wěn)定。線動量守恒只要求粒子插值具有0階連續(xù)性，而角動量守恒，則要求粒子插值近似至少具有一階連續(xù)性。Liu[22]給出函數(shù)及其前兩階導數(shù)具有n階核函數(shù)近似精度必須滿足的條件。0階連續(xù)性條件即是核函數(shù)歸一化條件，1階連續(xù)性則要求

(20)

然而，它們是由積分推導出的，不能保證離散化后方程的連續(xù)性。如圖2所示：當粒子位于邊界附近時，它的支持域就被邊界截斷而不完整，這時是不能保證離散后的方程滿足插值連續(xù)的。另外，當粒子的分布極不規(guī)則時，也不能保證插值連續(xù)性[23]。

圖2　核函數(shù)被邊界截斷

2.2　張力不穩(wěn)定

早期學者發(fā)現(xiàn)，即使?jié)M足CFL條件，當粒子之間的應力為負(即為張力)的時候，粒子會相互吸引，形成大量粒子“團聚”，最終造成計算失敗。最近研究表明，這種現(xiàn)象也可能在應力都是正的情況中出現(xiàn)，在彈性體的模擬中則表現(xiàn)為非物理的斷裂[24]。Swegle等[25,26]采用Von Neumann判據(jù)對張力不穩(wěn)定進行了深入研究。發(fā)現(xiàn)造成張力不穩(wěn)定的根本原因是由于核插值改變了粒子之間的本構關系，從而，改變了原始偏微分方程的性質。并提出了張力不穩(wěn)定性判據(jù)：

△WT>0

(21)

很多學者提出了解決張力不穩(wěn)定性問題的方法。Morris[4]研究發(fā)現(xiàn)高階樣條核函數(shù)與高斯型核函數(shù)能顯著降低不穩(wěn)定性的增長速度。然而，高階核函數(shù)在緊支域的某些地方為負值，將產(chǎn)生非物理結果[7]。Johnson[27]提出二次光滑函數(shù)來模擬高速沖擊問題，并稱可解決壓力不穩(wěn)定，但是因△W恒大于零，故不能解決拉伸失穩(wěn)。

Swegle等[28]等提出守恒光滑法，通過平滑2h波長的振動替代常用的人工粘性力來消除張力不穩(wěn)定性。Vignjevic[29]稱CSA法與施加人工排斥力法會影響材料強度，故使用時要引起注意。Dyka[30,31]等提出了“應力點法”的思想是在SPH粒子外設置應力點粒子來求應力和其它變量信息，而SPH粒子自身的運動學變量信息例如，位移、速度、質量等不變。該應力點方法和EFG方法中的“高斯積分點”類似[32]。Vignjevic[29]利用正則化核函數(shù)將“應力點法”擴展到二維問題。Monaghan[24]提出一種人工排斥力的方法來消除張力不穩(wěn)定性。

Randles與Libersky[33]提出“雙粒子動力學”方法來解決張力不穩(wěn)定性。另外，Liu[34]提出的RKPM，Chen等[35-38]提出了CSPM等方法，能改善甚至完全消除張力不穩(wěn)定性。Hu[39,40]在他的計算中采用了“參考壓力法”來消除張力不穩(wěn)定性。

2.3　邊界條件實施困難

因核函數(shù)被邊界截斷， SPH在邊界條件的處理上較其它數(shù)值方法更困難。近期主要使用的邊界施加方法有：

1) 排斥力方法

排斥力方法理論來源于計算分子力的Lennard-Jones方程[17]。該方法是在固壁邊界設置一組虛粒子，對靠近邊界的內(nèi)部粒子施加一定的排斥力，防止它穿透邊界。后來研究發(fā)現(xiàn)，早期的模型會導致粒子沿邊界“游走”，與實際情況不符。為此，Monaghan[41-43]對排斥力模型進行了改進。該方法的優(yōu)點是不受邊界形狀的影響，實施方便，缺點是模型的參數(shù)依然憑經(jīng)驗設置。

2) 鏡像粒子法

鏡像粒子方法是在每個時間步，距離邊界附近一定范圍內(nèi)，將內(nèi)部粒子以邊界為對稱面生成相應的鏡像虛粒子[44]。Takeda[8]，Morris[18]分別對鏡像粒子方法進行改進用來處理曲邊界。虛粒子邊界的守恒性好，但對于復雜形狀邊界(比如鋸齒狀邊界)或兩相交界面，這種方法實施起來有困難[43]。

3) 靜態(tài)粒子邊界

該方法與鏡像粒子法不同，其中虛粒子是固定的，虛粒子要參與內(nèi)部粒子的運算，虛粒子層可看作固壁的加厚[14,15,45]。Dalrymple與Knio[46]提出了被稱作“Dynamic Boundary"的一種雙層交錯粒子來模擬邊界。與Koshizuka的“DUMMY”粒子相似，邊界粒子與內(nèi)部粒子一樣要參與N-S方程計算，但是，位置固定不變或另外施加(例如：造浪器)。該方法被證明守恒性好，邊界實施容易[47]。

3　對經(jīng)典SPH方法的改進

3.1　NSPH( Normalized SPH)

Randles等[48]在Johnson等[49]的工作的基礎上，基于歸一化的思想提出了對核函數(shù)及導數(shù)的修正方法

(22)

(23)

式(22)與(23)分別保證了核函數(shù)插值的0階與1階連續(xù)。Bonet等[20]利用變分原理推導得到相同的核函數(shù)導數(shù)修正公式，保證角動量守恒。Belytschko[19]認為，該方法雖然能保證修正的核函數(shù)插值線性連續(xù)性，但修正的核函數(shù)導數(shù)不具有可積性，會導致采用Galerkin法建立的離散格式不能通過分片試驗。

3.2　RKPM(Reproducing kernel method，RKPM)

RKPM是在SPH方法的基礎上發(fā)展起來的一種基于核近似的無網(wǎng)格方法。由于SPH存在兩個不足：粒子插值難滿足線性一致性條件；當總粒子數(shù)目相對較少時，計算精度降低。Liu等針對這些問題對SPH方法進行改進，提出RKPM[34]。RKPM利用修正的核函數(shù)來保證粒子插值連續(xù)性，構造修正的核函數(shù)：

(24)

由n階再生條件，建立離散線性方程組，解方程組可以得到修正的核函數(shù)。因為高階再生核函數(shù)計算量大，在實際應用中線性再生核用得較廣泛。RKPM采用的是全局Galerkin離散格式。RKPM方法不僅解決了SPH法在邊界上的不連續(xù)性，而且完全消除了張力不穩(wěn)定性。Liu還將它與小波分析技術結合，使得RKPM具有多尺度自適應分析能力。

RKPM法已經(jīng)有大量的應用研究：如結構動力學[50]、應力集中[51]、大變形[52]、剪切帶的形成[53]等。

3.3　CSPM(Corrective smoothed particle method)

CSPM方法是Chen等[35-38]提出的。以二維直角笛卡兒坐標系為例，將函數(shù)f在點(xi,yi)處展開，并在等式兩邊同時乘以核函數(shù)W，并在求解域Ω內(nèi)積分得

忽略所有求導項，則得到核近似和粒子求和近似公式

(26)

明顯，式(26)與式(22)是等價的。通常，若要得到函數(shù)的一階導數(shù)近似，只需要式(25)中的核函數(shù)具有反對稱性(通常采用核函數(shù)一階導數(shù)W,x)即可消除0次導數(shù)項與二次導數(shù)項，高階導數(shù)求解依此類推。該方法中使用的核函數(shù)不要求與經(jīng)典SPH方法的核函數(shù)具有相同的性質，例如：經(jīng)典光滑核函數(shù)要求非負，在CSPM方法中則不要求，此方法的核函數(shù)選擇更靈活，但要求保證聯(lián)立方程組的系數(shù)矩陣不會奇異。CSPM在求解域內(nèi)具有二階精度，在邊界附近也具有一階精度。該方法求解函數(shù)及其導數(shù)的公式與NSPH方法相同，但這兩種方法的思想是不同的，前者更接近數(shù)學，后者更接近物理。CSPM方法具有三個優(yōu)點：理論上它可以方便的求解任意階導數(shù)；它在邊界條件的實施上較經(jīng)典SPH方法簡單，邊界粒子的參數(shù)可以直接設置；它在將非線性偏微分方程離散化時，較經(jīng)典SPH方法更直接。Zhang等[54]在推導函數(shù)及其導數(shù)的計算公式時，均保留高階導數(shù)項，建立聯(lián)立方程組，統(tǒng)一求解。這種方法精度較高，但是計算量更大。Fang[21]為了平衡精度與效率，將CSPM方法與經(jīng)典SPH方法耦合，在邊界附近采用CSPM方法求解以提高精度與穩(wěn)定性，在計算域內(nèi)部采用SPH方法來提高計算效率。

3.4　MLSPH(Moving least squares particle hydrodynamics)

MLS(Moving least squares)近似最早是由數(shù)學家設計用于數(shù)據(jù)擬合以及表面構造[55]，利用MLS構造的形函數(shù)可積，適當選擇權函數(shù)可以保證MLS形函數(shù)滿足全局連續(xù)性[19]。MLS已經(jīng)被廣泛用來構造無網(wǎng)格法的形函數(shù)，如EFG，F(xiàn)PM等。Dilts[56,57]提出移動最小二乘SPH方法MLSPH。

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

(27)

(28)

式中，p(rj)是基函數(shù)向量在rj處的值，以二次多項式基為例

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

(29)

pT是計算點ri的鄰域內(nèi)各點的基函數(shù)向量組成的矩陣。

MLS形函數(shù)可以看作對經(jīng)典SPH形函數(shù)，在粒子分布不規(guī)則時的一種修正，當粒子規(guī)則分布時有

(30)

MLSPH方法能提高精度和消除張力不穩(wěn)定性，然而，Monaghan[24]指出它的計算量大約是經(jīng)典SPH方法的8倍。RKPM，NSPH與MLSPH均是對經(jīng)典對稱非負核函數(shù)及其導數(shù)進行修正或重構來提高精度，而重構

的核函數(shù)在緊支域中可能出現(xiàn)負值，這樣就產(chǎn)生了一些場變量出現(xiàn)非物理值，比如：負密度，負能量，這樣會造成計算崩潰。而CSPM不改變經(jīng)典SPH的核函數(shù)，不會出現(xiàn)核函數(shù)出現(xiàn)負值情況。

3.5　ISPH(Incompressible SPH)

Cummins[58]綜合移動粒子半隱式方法(MPS)與投影方法(Projection)的優(yōu)點，提出了一種投影SPH方法。為了區(qū)別于WCSPH，該類方法后來被稱作“ISPH"。該方法屬于分數(shù)步法，在中間臨時步，可以先不考慮不可壓縮條件，即先假設粒子不受壓力而只受粘性力或外力，當粒子運動后，它的密度或速度散度會發(fā)生變化，再根據(jù)不可壓縮條件可以推導出壓力Possion方程

(31)

求解方程得到壓力值，最后，再考慮壓力的影響，相當于把速度和位置都改變的粒子再“壓”回去。與Projection方法采用有限差分離散不同，ISPH函數(shù)的離散是通過光滑核函數(shù)插值離散完成的。近年來的研究表明ISPH在收斂性，穩(wěn)定性，計算精度，計算效率，壓力值求解諸方面均較經(jīng)典SPH方法好[14-16]。近期有不同學者對ISPH和經(jīng)典SPH的應用情況進行了比較[59-62]。發(fā)現(xiàn)，如果對經(jīng)典WCSPH算法計算的密度每隔20步進行MLS修正，計算結果與ISPH相當。

4　SPH的應用

4.1　可壓縮流體計算

SPH方法最早是由天文學家在計算天體物理中提出來的，因為天體運動可以被看作是可壓縮流[1,2]。Springel[63]基于SPH方法開發(fā)了一套用于計算天體的軟件，開展了大量計算，例如星系碰撞、宇宙結構的形成、千年計算等。SPH在可壓縮流計算領域的另一個重要應用是對爆炸的模擬。因為在爆炸過程中，存在著極大變形和自由表面，傳統(tǒng)有網(wǎng)格方法難以勝任。特別是對一些細節(jié)，例如，破碎、飛濺、融合等的模擬，無網(wǎng)格粒子法表現(xiàn)更佳。Liu[64]等將SPH方法應用到水下爆炸沖擊問題中。

4.2　不可壓縮流體計算

1994年，Monaghan首次將SPH方法應用到自由表面流的模擬，也為SPH方法在不可壓縮流中的計算開了先河[17]；隨后,世界各地學者紛紛加入到這個領域中來，成立一個SPH歐洲研究興趣組織(SPH European Research Interest Community, SPHERIC)。SPHERIC的成員開展了大量研究，其中應用最廣泛的領域是在海洋工程，例如波浪對海岸(或石油鉆井平臺)的沖擊、船艙中的液體晃蕩。在工業(yè)領域應用較為成功的是鑄造充型過程的模擬。Cleary等[65]采用SPH方法對高壓壓鑄充型過程進行了研究，與VOF方法計算結果與實驗進行比較，表明SPH可以有效地計算這種高速流動現(xiàn)象。國內(nèi)，曹文炅等[66]采用類似方法對高壓壓鑄充型過程進行了分析。SPH方法在工業(yè)領域復雜加工中的應用有很多，例如，采用ISPH算法模擬聚合物螺桿擠出展開流道[67]、模擬磨料水射流中單磨粒加速過程[68]、激光水下加工中，水與熔融金屬之間的相互作用[69]。長期阻礙其應用發(fā)展的兩大主要因素有：復雜加工領域中的計算域幾何形態(tài)復雜，且都是動態(tài)變化的；SPH作為一種無網(wǎng)格粒子法，其計算量遠遠大于傳統(tǒng)網(wǎng)格方法，特別是在粒子搜索階段，每個時間步，都需要重新進行鄰域粒子搜索配對，而有網(wǎng)格法，一旦網(wǎng)格固定，節(jié)點與節(jié)點的連接是固定的。最近，SPHERIC成員提出了一種前處理器，首先利用已經(jīng)非常成熟的商業(yè)網(wǎng)格法生成器進行復雜計算域的網(wǎng)格劃分，然后再將網(wǎng)格區(qū)域用粒子代替；他們采用了基于GPU的多核并行計算方案，計算效率比單CPU版的程序提高了60倍[70]。

4.3　材料的變形與沖擊斷裂

材料加工過程中，存在材料的大變形。Libersky[44,48,71]等將SPH應用于鍛壓過程的計算。近年來研究者開始將SPH法應用到切削過程的模擬中．Heinstein等[72]證明了SPH法在直角切削模擬中的有效應用．Limido等[73]采用SPH法對高速切削過程進行模擬，分析了切屑的形成并進行了切削力預報。Spreng[74]采用采用自適應SPH算法模擬切屑過程，避免了數(shù)值斷裂。我國學者也開展了很多研究。例如，郭曉光[75]等將SPH方法應用于超精密切削過程仿真。

Johnson[27, 48,71]等最早將它應用到高速碰撞問題的研究中。徐金中等[76]對初始光滑長度和粒子間距對計算結果的影響進行了深入研究。張志春等[77]采用基于SPH-FEM耦合的算法對子彈沖擊高強鋼板進行了數(shù)值計算。

5　展望

(1) 工程領域的應用基本都是三維問題，有時甚至是多相情況下的三維問題，邊界實施較二維更為復雜。采用排斥力設置邊界的優(yōu)點是對復雜形態(tài)邊界的適應性強，缺點是參數(shù)的設置一直沒有一套科學的，簡單的，普適性的方案。值得進一步深入探索。

(2) 對于無網(wǎng)格粒子法，計算時間會隨著粒子數(shù)量增加急劇增加，這是無網(wǎng)格粒子法在大型復雜加工過程中的應用的一個瓶頸。除了依靠并行算法以外，是否能有從算法本身進行新的改進措施，特別是對于每個時間步都要重新進行的粒子鄰域搜索，是計算耗時一個關鍵因素，有待深入探索。

(3) SPH方程在不同工業(yè)領域的應用中不能完全通用，需要結合實際情況，進行修正，例如，在聚合物領域中的應用，要考慮流體粘彈性問題，在金屬切屑加工中的應用要考慮金屬剝離后的粒子鄰域搜索。

參考文獻：

[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. 光滑粒子流體動力學——一種無網(wǎng)格粒子法[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方法建模及數(shù)值模擬[J]. 華南理工大學學報(自然科學版), 2011, 39(1): 100-105.

[67] 董添文, 江順亮, 黃興元, 柳和生, 基于不可壓縮光滑粒子流體動力學的單螺桿擠出三維流動數(shù)值模擬[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].強度與環(huán)境, 2009, 36(1):1-7.

[77] 張志春, 強洪夫, 傅學金等. 基于SPH-FEM 轉換算法的7. 62 mm步槍彈沖擊30CrMnSiA 鋼板的數(shù)值計算[J]. 計算物理，2012, 29(1):73-81.