，

(1.中交第二航務(wù)工程局有限公司，武漢 430040；2.中交公路長大橋建設(shè)國家工程研究中心有限公司，北京 250073;3.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240)

1 引 言

人工堆積的尾礦壩和粉砂質(zhì)邊坡等土工構(gòu)筑物極容易發(fā)生滑移失穩(wěn)等災(zāi)害，往往會造成極大的人身和財產(chǎn)損失。準確模擬這類土體大變形問題，對工程設(shè)計和災(zāi)害防御具有重大意義。土力學(xué)中常見的數(shù)值模擬方法如FEM(有限元法)和DEM(離散元法)在求解這類問題時都存在缺陷。FEM常因網(wǎng)格束縛導(dǎo)致計算結(jié)果不準確；DEM則由于選取計算參數(shù)困難，較難得出精確的數(shù)值模擬結(jié)果[1]。光滑粒子流體動力學(xué)SPH(Smoothed Particle Hydrodynamics)方法是一種拉格朗日型無網(wǎng)格粒子方法，起初主要用于求解天體物理學(xué)問題[2,3]。SPH方法將問題域離散成有限個攜帶物質(zhì)信息的粒子，每個粒子都在周圍粒子影響下遵守守恒定律運動，十分適用于求解大變形問題。目前SPH已廣泛應(yīng)用于流體力學(xué)及固體力學(xué)等各個領(lǐng)域，成功模擬了爆炸、穿甲及潰壩等大變形問題。近年來，越來越多的學(xué)者開始應(yīng)用SPH方法模擬土體大變形問題。Bui等[4]將SPH方法應(yīng)用在土體大變形問題中，采用Drucker-Prager屈服模型，模擬了不同內(nèi)摩擦角的土坡滑移現(xiàn)象。黃雨等[5]基于等效牛頓黏度系數(shù)的概念，采用Bingham模型描述土體大變形，模擬了土體流動。Chen等[6]將各向異性的Drucker-Prager模型與修正Kondner和Zelasko模型應(yīng)用到SPH中，用以分析土體的非線性特性及邊坡的動力反應(yīng)和滑移特征。Wang等[7,8]基于混合物原理，建立了兩相流模型，模擬了射流沖刷和水下滑坡問題。

在SPH方法中，粒子的場函數(shù)由支持域內(nèi)其他粒子場函數(shù)的加權(quán)平均而得，其中權(quán)函數(shù)是具有偶函數(shù)性質(zhì)的光滑函數(shù)。因此，若想獲得較高的計算精度，需支持域內(nèi)具有足夠多的粒子，且粒子整體呈對稱分布。但這兩條性質(zhì)在邊界附近都很難滿足[9]。該邊界附近粒子缺失的問題，需要通過邊界處理方法來解決。常用的邊界處理方法包括排斥力法和虛粒子法，排斥力法易受初始擾動的影響[10]，而虛粒子法則難以處理復(fù)雜邊界[11]，且存在邊界零粒子層問題[12]。近年來，一種新興的邊界處理方法正引起學(xué)者們的注意。Ferrand等[13]提出了能夠應(yīng)用于任意形狀邊界的邊界處理方法，即統(tǒng)一半解析壁面邊界條件USAW(unified semi-analytical wall boundary conditions)處理方法，成功模擬了潰壩流問題。Leroy等[14]改進了USAW邊界處理方法，并將其與ISPH方法結(jié)合。Mayrhofer等[15]應(yīng)用USAW方法成功模擬了三維潰壩流問題。Cercos-Pita等[16]提供了一套基于GPU并行計算的USAW-SPH開源代碼。文獻[17]基于USAW方法，提出了兩種不規(guī)則粒子分布，更好地模擬了復(fù)雜邊界潰壩問題。USAW邊界處理方法在流體動力學(xué)上的成功，給土體大變形問題的邊界處理提供了新思路。本文首次應(yīng)用USAW邊界處理方法模擬土體大變形問題。

以往USAW方法中，研究者一般生成規(guī)則的初始粒子分布，其中內(nèi)部粒子質(zhì)量相同，而邊界粒子質(zhì)量取決于該點在邊界處所形成的角度[13]。但在這種方法中，邊界粒子的質(zhì)量會影響問題域的模擬，并且會出現(xiàn)邊界零粒子層，影響模擬精度。結(jié)合USAW基本理論，本文提出了邊界粒子無體積和無質(zhì)量的做法。

2 土體大變形理論模型

土體的控制方程由質(zhì)量守恒和動量守恒方程構(gòu)成：

(1，2)

定義土體為理想彈塑性材料，且服從Drucker-Prager屈服準則。在該模型中，屈服函數(shù)可表示為

(3)

(4，5)

本文算例是大變形問題，選取與剛體旋轉(zhuǎn)無關(guān)的Jaumann應(yīng)力率,

(6)

(7)

彈塑性材料可將材料的應(yīng)力應(yīng)變關(guān)系分解為彈性部分和塑性部分。彈性部分采用胡克定律求解，塑性部分則遵守塑性流動法則，最后獲得本構(gòu)方程[7]：

(8)

(9)

(10)

式中Ψ為塑性勢函數(shù)中的剪脹角，參考文獻[4]取值為0，本文選用的塑性勢函數(shù)g為

(11)

3 SPH理論和邊界處理方法

在SPH方法中，將問題域離散成多個粒子，每個粒子具有體積和質(zhì)量等屬性，并攜帶應(yīng)力、速度和位置等信息，能夠在力的作用下運動。粒子a的場函數(shù)fa(如壓強或速度)及其導(dǎo)函數(shù)fa，均可由其支持域內(nèi)各粒子的場函數(shù)插值近似得到，

(12)

Wa b

(13)

式中Vb為粒子b的體積，P為支持域內(nèi)所有粒子點的集合，本文符號較多，為避免混淆，將其列入表1。在SPH方法中，目標粒子的場函數(shù)通過周圍粒子對應(yīng)函數(shù)的加權(quán)平均求得，見式(12)。式(13)是SPH方法最基礎(chǔ)的導(dǎo)函數(shù)近似式，該式將導(dǎo)函數(shù)的加權(quán)平均轉(zhuǎn)換成了核函數(shù)導(dǎo)數(shù)的加權(quán)平均，方便了導(dǎo)函數(shù)的求解。Wa b=W(xa b,h)為光滑核函數(shù)，是SPH方法加權(quán)近似的權(quán)函數(shù)，具有緊支性和歸一性，其大小隨粒子a與粒子b之間相對距離的增大而減??；xa b=xa-xb，xa和xb分別指粒子a和b的位置矢量，h為光滑長度。本文選取Wendland五階核函數(shù)：

(14)

式中Ra b=xa b/h；αd為保證光滑核函數(shù)歸一性的量，在二維算例中，其值為7/(4πh2)。更詳細的SPH理論推導(dǎo)可參考文獻[7]。

表1 符號含義

Tab.1 Symbol definition

符號 含義a,b粒子(含邊界粒子)e僅邊界粒子s邊界上兩個粒子之間的線段P支持域所有粒子的集合Ω整個計算域F內(nèi)部粒子的集合(除邊界)S邊界上所有線段的集合

嚴格來講，式(12,13)只適用于遠離邊界處的粒子，當(dāng)粒子靠近邊界時，由于邊界的截斷，粒子的支持域內(nèi)將缺少部分粒子。為了彌補邊界粒子缺失，USAW邊界處理方法在插值近似時引入了修正因子(renormalization factor)：

(15)

式中γa為粒子a的修正因子，定義為

(16)

式中Ω為粒子a的支持域，W(xa-x,h)為粒子a的光滑函數(shù)。當(dāng)粒子a遠離壁面時，支持域完整，γa=1；靠近壁面時，由于邊界的截斷，支持域破缺，導(dǎo)致γa<1。γa可采用時間步進法[13]或解析方法[17]求解，本文采用后者。在計算函數(shù)導(dǎo)數(shù)時，與式(13)不同，USAW邊界處理方法保留了邊界積分項，

(17)

(18)

式中l(wèi)為線段s的長度。

邊界粒子的密度可通過SPH插值近似得到，

(19)

(20)

而粒子a處速度的偏導(dǎo)，可由式(21)離散得到。

(21)

式(21)采用常用的SPH變化技巧[7]，具有一定的反對稱性，符合作用力與反作用力原理。粒子a上其他場函數(shù)的偏導(dǎo)數(shù)可做類似處理。

傳統(tǒng)的SPH方法大多假定邊界粒子與內(nèi)部粒子屬性一致，使用相同的方法來求解密度和壓強；在USAW方法中，則采取不同的方法。

內(nèi)部粒子的密度，通過修正后的密度求和法[13]來求解：

(22)

式中 上標n指的是時間步。式(22)是運用傳統(tǒng)的密度求和法，將兩個相鄰時間步的密度求解方程相減得到[13]。同時，為了減少數(shù)值振蕩，本文在方程中添加了δ-SPH 技術(shù)[18]，即式中最后一項，其中，δ為控制修正大小的一個值，本文取0.1。

為了減少數(shù)值振蕩，除上文介紹過的δ-SPH 技術(shù)外，本文還采用了Marrone等[19]提出的人工粘性技術(shù)。其不僅能夠?qū)幽苻D(zhuǎn)變?yōu)闊崮?，同時也提供了沖擊的耗散，避免了粒子互相靠近時的穿透現(xiàn)象。在模擬過程中發(fā)現(xiàn)，固體粒子在滑坡過程中容易出現(xiàn)聚集成塊的現(xiàn)象，即所謂的張力不穩(wěn)定現(xiàn)象。為了減小這種張力不穩(wěn)定，本文采取了由Monaghan[20]提出的人工應(yīng)力法。在數(shù)值模擬中，有時應(yīng)力會躍出屈服面外，這是不允許的。本文采用張裂處理(Tension crack treatment)和比例拉回(scaling back procedure)的辦法將屈服面外的應(yīng)力映射回屈服面上[4]，并采用二階精度的蛙跳法(Leap -Frog)進行時間步進。具體可參見文獻 [7,8]。

圖1 USAW中的邊界定義

Fig.1 Boundary definition in USAW method

關(guān)于SPH方法的粒子分布，一般先將問題域劃分為若干個相鄰的矩形體積單元，再取矩形單元的中心位置為粒子位置，粒子的體積等于對應(yīng)單元的體積。除問題域外，還需布置若干粒子來模擬邊界線，這些粒子是位于問題域外的虛粒子。虛粒子雖相對于邊界靜止，但同樣參與計算，如圖2(a)所示。圖中黑色實心點為內(nèi)部粒子，其體積之和為模擬的問題域，灰色粒子點為邊界粒子，虛線為問題域的邊界。該方法由于邊界粒子無法與壁面貼合，較難模擬復(fù)雜壁面形狀。

USAW邊界處理方法可以處理任意形狀邊界，部分原因是該方法直接假定邊界粒子的連線作為壁面，如圖2(b)所示。圖中邊界粒子僅在問題域內(nèi)部分上色，是為了反應(yīng)傳統(tǒng)USAW方法中邊界粒子的質(zhì)量分布[13]，即邊界處粒子質(zhì)量一般是內(nèi)部粒子的一半，直角拐點處粒子質(zhì)量是內(nèi)部粒子質(zhì)量的1/4。在這種分布方法中，邊界粒子體積占據(jù)了一定內(nèi)部問題域的體積，即相當(dāng)于在邊界處固定了一排靜止不動的內(nèi)部粒子，從而導(dǎo)致問題域小于需模擬的尺寸，是不合理的。

本文發(fā)現(xiàn)USAW邊界處理方法中，邊界粒子可以不具有體積和質(zhì)量。USAW方法在求解導(dǎo)函數(shù)時，將積分項分成了體積積分項和邊界上的面積分項，在邊界的面積分項中，僅邊界粒子的密度參與了計算，見式(17)。同時在邊界處密度和壓強的求解式(25,26)中，邊界粒子的質(zhì)量或體積都不參與計算。故在邊界處采用這種無體積和無質(zhì)量的邊界粒子十分合理。

圖2 有質(zhì)量邊界粒子的分布

Fig.2 Distribution of boundary particles with mass in different methods

本文提出邊界粒子無體積和無質(zhì)量的處理方法，如圖3所示，其中空心粒子為邊界粒子。USAW邊界處理方法將場函數(shù)導(dǎo)數(shù)分為體積分和面積分的累加，而這種邊界粒子無質(zhì)量的分布方法，也對應(yīng)地將粒子區(qū)分為有體積和無體積兩部分，具有更高的精度，也更符合USAW邊界處理方法的本質(zhì)。使用無體積和無質(zhì)量的邊界粒子，既能保證邊界粒子與壁面完美貼合，正確模擬復(fù)雜形狀，又能準確地求解導(dǎo)函數(shù)中的面積分。

4 數(shù)值計算結(jié)果分析

模擬沙土滑移，基本模型如圖4所示，在一個矩形槽左下角，存在一個矩形沙柱，初始時刻沙粒靜止，當(dāng)快速拔除右邊的擋沙板之后，沙柱在重力作用下垮塌。其中灰色矩形是初始時刻的沙土，hi和di分別為初始高和寬，高寬比為k；虛線是滑坡終止時的沙面形狀，其中h∞和d∞分別為滑坡結(jié)束時沙土的最大高度和最大寬度，δd為滑坡終止時滑移的最遠距離。

本文模擬不同初始高寬比k下的滑坡運動，分析初始高寬比對滑移最遠距離和最終沙粒高度的影響，并與Lube等[21]測量的實驗值進行比較，驗證算例的可靠性。在SPH模擬中，采用與實驗所用的粗石英砂相同的參數(shù)，即內(nèi)摩擦角為31°，粘聚力為0。長寬比k分別為0.5，1和3，初始長度di恒為0.2 m，各算例初始粒子間距均為0.002 m，光滑長度均取初始粒子間距的1.2倍，3個模型的總粒子數(shù)分別為5350，10650和31950，計算時間步長為5×10-5s，計算結(jié)果如圖5所示。

圖3 無質(zhì)量邊界粒子分布

Fig.3 Distribution of boundary particles without mass

圖4 潰沙模型和基本符號

Fig.4 Landslide model and basic symbols

3個算例的滑移最遠距離d∞分別為0.35 m，0.51 m和1.13 m，最終沙粒高度h∞分別為0.1 m，0.19 m和0.29 m。Lube等[21]分析了滑坡的實驗數(shù)據(jù)，總結(jié)出滑移距離δd與初始寬度di的無量綱比值ζ(k)可表示為

(23)

式中c1=1.6,c2=2.2；當(dāng)k位于1.8～2.8之間時，沒有準確函數(shù)可表達。滑移后左端高度h∞與初始寬度di比值ξ(k)可表示為

(24)

數(shù)值計算值和實驗值的比較結(jié)果列入表2，相對誤差計算方法為實驗值與模擬值之差的絕對值除以實驗值。由表2可知，模擬結(jié)果與實驗值的ζ(k)和ξ(k)相對誤差較小，在5%左右，屬工程誤差范圍之內(nèi)，故可認為本文方法能夠較好地模擬潰沙過程。

圖5 不同高寬比下的潰沙終止形態(tài)

Fig.5 Landslide simulations with different aspect ratios

從圖6(a)可以看出，內(nèi)摩擦角增大時，滑移距離減小。從圖6(b)可以看出，黏聚力增大時，滑移距離減?。火ぞ哿^續(xù)增大時，沙體的表面將變得粗糙，沙柱右上角在滑動過程中在沙面上形成一個尖點，滑動終止時也無法消失，且黏聚力越大，最后保留的尖點也越大。

SPH方法模擬運動學(xué)問題時，常出現(xiàn)邊界零粒子層現(xiàn)象[12]，在運動介質(zhì)的舌尖與邊界之間，經(jīng)常會出現(xiàn)一個大約為1個粒子大小的空隙。該零粒子層會影響邊界層內(nèi)的數(shù)值模擬精度，而且粒子數(shù)越少，單個粒子所占體積越大，零粒子層就越厚。本文采用三種邊界處理方法，模擬了圖5(a)的潰沙問題，并觀察其舌尖與邊界之間是否存在零粒子層現(xiàn)象，如圖7所示。

表2 模擬值與實驗值之間的相對誤差

Tab.2 Relative error between simulated and experimental values

Kζ(k)模擬值ζ(k)實驗值相對誤差/%ξ(k)模擬值ξ(k)實驗值相對誤差%0.50.750.86.250.50.5011.551.63.130.951534.654.61.11.351.556.4

圖6 不同參數(shù)下的模擬結(jié)果

Fig.6 Simulation results under different soil parameters

可以看出，使用虛粒子法模擬時，粒子舌尖與邊界之間明顯存在零粒子層現(xiàn)象；在傳統(tǒng)的USAW邊界處理方法中，依然可以觀察到邊界零粒子層的存在，但較虛粒子法稍薄，主要原因是其邊界粒子質(zhì)量僅為內(nèi)部粒子質(zhì)量的1/2；圖7(c)顯示，本文提出的邊界粒子無質(zhì)量的USAW方法有效消除了邊界零粒子層現(xiàn)象，而且其舌尖位置更加接近實驗位置。

探索滑坡沖擊楔形體的運動過程，并分析滑坡對楔形體的沖擊力。如圖8(a)所示，在問題域的左下角，存在矩形沙土，沙土高為0.6 m，寬為0.2 m，底部距原點0.3 m處存在楔形體凸起，楔形體呈等邊直角三角形形狀，斜邊長為0.18 m，方向朝下。SPH模擬的算例共由8400個粒子組成，其余參數(shù)均與上文保持一致。

滑動前0.8 s內(nèi)的沙體輪廓如圖8(a)所示?；麻_始0.14 s后，滑坡才抵達楔形體底部，隨后越來越多沙土覆蓋楔形體左側(cè)，楔形體左側(cè)壓力相應(yīng)急劇上升；t=0.22 s時，沙土剛好全部覆蓋楔形體；此后，部分沙土沖過楔形體，形成一個勾形舌尖；t=0.34 s時，舌尖落入地面，與楔形體形成一個空穴，并隨后填滿；然后沙土繼續(xù)向前滑行，速度不斷減慢，約進行到0.8 s時，停止滑行，最遠距離達到1.1 m，整個過程中有少數(shù)粒子脫離舌尖，以較大速度拋物前行。

圖7 邊界零粒子層

Fig.7 Zero particle layer near boundary

圖8 楔形體潰沙輪廓圖

Fig.8 Profiles of landslide on a wedge

沙土沖擊建筑物所造成的壓力是工程設(shè)計關(guān)注的重點。該算例中，楔形體左側(cè)壓力隨時間變化如圖9所示。可以看出，滑坡開始0.14 s后，沙土接觸楔形體并對楔形體施壓，楔形體左側(cè)壓力對應(yīng)急劇上升；0.22 s時，沙土剛好全部覆蓋楔形體，楔形體左側(cè)所受壓力達到最大值，約為15.5 kN；此后，沙土沖過楔形體，對應(yīng)的楔形體左側(cè)壓力急劇減小，在0.4 s時，沙粒與楔形體右側(cè)壁面相互作用，從而導(dǎo)致左側(cè)壁面壓強形成一個小波峰；然后，楔形體左側(cè)沙土減少，且沙土流過楔形體的速度逐漸降低，左側(cè)壓強緩慢減小，最終靜止時壓力值約為3 kN。邊界壓力的強烈振蕩一直是限制SPH方法發(fā)展的重要阻礙，本文方法模擬的邊界壓力振蕩較小。

圖9 楔形體左側(cè)壓力隨時間變化情況

Fig.9 Time evolution of the pressure on the left ide of the wedge

5 結(jié) 論

SPH方法是拉格朗日型無網(wǎng)格粒子法，十分適用于處理大變形問題。本文應(yīng)用SPH方法，結(jié)合改進后的USAW邊界處理方法，較精確地模擬了滑坡問題，具體結(jié)論如下。

(1) 探討并改進了USAW邊界處理方法，首次使用梯形法求解γa s，并將傳統(tǒng)的SPH數(shù)值技術(shù)如 δ-SPH技術(shù)、人工粘性以及人工應(yīng)力與USAW邊界處理方法相結(jié)合，改善了USAW方法的計算效果。

(2) 提出了邊界粒子無質(zhì)量的處理方法，該方法適用于模擬復(fù)雜邊界問題，并提高了USAW邊界處理方法的精度，還能避免邊界零粒子層現(xiàn)象。

(3) 首次應(yīng)用USAW邊界處理方法模擬滑坡問題，證明了該方法模擬土體大變形問題的可行性，分析了內(nèi)摩擦角和黏聚力等土體物理特性參數(shù)對滑坡過程的影響，并計算模擬了滑坡沖擊楔形體過程，研究了楔形體所受壓力的變化規(guī)律。

本文將改進的SPH方法應(yīng)用于滑坡問題計算，通過比較計算結(jié)果與實驗結(jié)果，驗證了方法的有效性。但對于土力學(xué)中大量的其他問題并沒有涉及，這將是今后的研究內(nèi)容之一。

:

[1] 潘建平，曾慶筠.SPH方法在土工大變形分析中的應(yīng)用研究進展[J].安全與環(huán)境學(xué)報，2015，15(5):144-150.(PAN Jian-ping，ZENG Qing-yun.Research review over the advances of the application of SPH method to the analysis of the serious geotechnical deformation[J].JournalofSafetyandEnvironment，2015，15(5):144-150.(in Chinese))

[2] Lucy L B.A numerical approach to the testing of the fission hypothesis[J].TheAstronomicalJournal，1977，82:1013-1024.

[3] Gingold R A，Monaghan J J.Smoothed particle hydrodynamics:theory and application to non-spherical stars[J].MonthlyNoticesoftheRoyalAstronomicalSociety，1977，181(3):375-389.

[4] Bui H H，F(xiàn)ukagawa R，Sako K，et al.Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics，2008，32(12):1537-1570.

[5] 黃 雨，郝 亮，謝 攀，等.土體流動大變形的 SPH 數(shù)值模擬[J].巖土工程學(xué)報，2009，31(10):1520-1524.(HUANG Yu，HAO Liang，XIE Pan，et al.Numerical simulation of large deformation of soil flow based on SPH method[J].ChineseJournalofGeotechnicalEngineering，2009，31(10):1520-1524.(in Chinese))

[6] Chen W，Qiu T.Simulation of earthquake -induced slope deformation using SPH method[J].InternationalJournalforNumericalandAnalyticalMe-thodsinGeomechanics，2014，38(3):297-330.

[7] Wang C，Wang Y，Peng C，et al.Smoothed particle hydrodynamics simulation of water-soil mixture flows[J].JournalofHydraulicEngineering，2016，142(10):04016032.

[8] Wang C，Wang Y，Peng C，et al.Two -fluid smoothed particle hydrodynamics simulation of submerged granular column collapse[J].MechanicsResearchCommunications，2017，79:15-23.

[9] 周 杰，徐勝利.SPH方法的運動邊界條件研究[J].計算力學(xué)學(xué)報，2016，33(3):412-417.(ZHOU Jie，XU Sheng-li.Research on the moving boundary condition in SPH methods[J].ChineseJournalofComputationalMechanics，2016，33(3):412-417. (in Chinese))

[10] 上官子檸，馬慶位，韓端鋒.SPH 邊界處理中兩種排斥力模型的對比分析[J].船舶力學(xué)，2014，18(1):37-44.(SHANGGUAN Zi-ning，MA Qing-wei，HAN Duan-feng.Comparisons of two repulsive models for boundary treatment in SPH [J].JournalofShipMechanics，2014，18(1):37-44.(in Chinese))

[11] 艾孜海爾·哈力克，熱合買提江·依明，開依沙爾·熱合曼，等.不同質(zhì)量粒子分布 SPH 方法及其應(yīng)用[J].振動與沖擊，2015，34(22):62-67.(AZHAR Halik，RAHMATJAN Imin，KAYSAR Rahman，et al.SPH method with different mass particle distributions and its application[J].JournalofVibrationandShock，2015，34(22):62-67.(in Chinese))

[12] Ni X Y，F(xiàn)eng W B.Numerical simulation of wave overtopping based on DualSPHysics [J].AppliedMechanicsandMaterials，2013，405:1463-1471.

[13] Ferrand M，Laurence D R，Rogers B D，et al.Unified semi-analytical wall boundary conditions for inviscid，laminar or turbulent flows in the meshless SPH method[J].InternationalJournalforNumericalMethodsinFluids，2013，71(4):446-472.

[14] Leroy A，Violeau D，F(xiàn)errand M，et al.Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[J].JournalofComputationalPhysics，2014，261:106-129.

[15] Mayrhofer A，F(xiàn)errand M，Kassiotis C，et al.Unified semi-analytical wall boundary conditions in SPH:analytical extension to 3-D[J].NumericalAlgorithms，2015，68(1):15-34.

[16] Cercos-Pita J L.AQUAgpusph，a new free 3D SPH solver accelerated with OpenCL[J].ComputerPhy-sicsCommunications，2015，192:295-312.

[17] 駱 釗，汪 淳.統(tǒng)一半解析邊界處理SPH方法與不規(guī)則粒子分布[J].水動力學(xué)研究與進展A輯，2017，32(2):189-197.(LUO Zhao，WANG Chun.Unified semi-analytical wall boundary treatment in SPH and irregular particle distribution[J].ChineseJournalofHydrodynamics，2017，32(2):189-197.(in Chinese))

[18] Leroy A，Violeau D，F(xiàn)errand M，et al.Unified semi-analytical wall boundary conditions applied to 2-D incompressible SPH[J].JournalofComputationalPhysics，2014，261:106-129.

[19] Marrone S，Antuono M，Colagrossi A，et al.δ-SPH model for simulating violent impact flows[J].ComputerMethodsinAppliedMechanicsandEnginee-ring，2011，200(13):1526-1542.

[20] Monaghan J J.SPH without a tensile instability[J].JournalofComputationalPhysics，2000，159(2):290-311.

[21] Lube G，Huppert H E，Sparks R S J，et al.Collapses of two -dimensional granular columns[J].PhysicalReviewE，2005，72(4):041301.