改進(jìn)的SPH邊界處理方法與土體大變形模擬

，

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

1 引 言

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

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

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

2 土體大變形理論模型

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

(1，2)

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

(3)

(4，5)

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

(6)

(7)

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

(8)

(9)

(10)

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

(11)

3 SPH理論和邊界處理方法

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

(12)

Wa b

(13)

式中Vb為粒子b的體積，P為支持域內(nèi)所有粒子點(diǎn)的集合，本文符號(hào)較多，為避免混淆，將其列入表1。在SPH方法中，目標(biāo)粒子的場函數(shù)通過周圍粒子對(duì)應(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之間相對(duì)距離的增大而減??；xa b=xa-xb，xa和xb分別指粒子a和b的位置矢量，h為光滑長度。本文選取Wendland五階核函數(shù)：

(14)

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

表1 符號(hào)含義

Tab.1 Symbol definition

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

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

(15)

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

(16)

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

(17)

(18)

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

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

(19)

(20)

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

(21)

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

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

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

(22)

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

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

圖1 USAW中的邊界定義

Fig.1 Boundary definition in USAW method

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

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

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

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

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

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

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

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

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

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

Fig.3 Distribution of boundary particles without mass

圖4 潰沙模型和基本符號(hào)

Fig.4 Landslide model and basic symbols

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

(23)

式中c1=1.6,c2=2.2；當(dāng)k位于1.8～2.8之間時(shí)，沒有準(zhǔn)確函數(shù)可表達(dá)?；坪笞蠖烁叨萮∞與初始寬度di比值ξ(k)可表示為

(24)

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

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

Fig.5 Landslide simulations with different aspect ratios

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

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

表2 模擬值與實(shí)驗(yàn)值之間的相對(duì)誤差

Tab.2 Relative error between simulated and experimental values

Kζ(k)模擬值ζ(k)實(shí)驗(yàn)值相對(duì)誤差/%ξ(k)模擬值ξ(k)實(shí)驗(yàn)值相對(duì)誤差%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

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

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

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

圖7 邊界零粒子層

Fig.7 Zero particle layer near boundary

圖8 楔形體潰沙輪廓圖

Fig.8 Profiles of landslide on a wedge

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

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

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

5 結(jié) 論

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

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

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

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

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

:

[1] 潘建平，曾慶筠.SPH方法在土工大變形分析中的應(yīng)用研究進(jìn)展[J].安全與環(huán)境學(xué)報(bào)，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] 黃 雨，郝 亮，謝 攀，等.土體流動(dòng)大變形的 SPH 數(shù)值模擬[J].巖土工程學(xué)報(bào)，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方法的運(yùn)動(dòng)邊界條件研究[J].計(jì)算力學(xué)學(xué)報(bào)，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 邊界處理中兩種排斥力模型的對(duì)比分析[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].振動(dòng)與沖擊，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].水動(dòng)力學(xué)研究與進(jìn)展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.