光滑粒子流體動力學方法固壁處理的一種新型排斥力模型*

韓亞偉 強洪夫 趙玖玲 高巍然

(第二炮兵工程大學601室,西安 710025)

(2012年8月10日收到;2012年9月16日收到修改稿)

1 引言

近年來,光滑粒子流體動力學(SPH)方法作為一種純Lagrangian型的無網(wǎng)格算法得到了很快發(fā)展,該算法最早是由Lucy[1]與Gingold[2]分別提出,最初用于解決天體物理中無邊界三維流動問題,現(xiàn)已推廣應用到計算力學諸多領域的數(shù)值仿真之中.該算法的特點是在計算中不需要生成任何輔助網(wǎng)格,因此在處理大變形問題時不存在Lagrangian網(wǎng)格算法中網(wǎng)格纏繞的技術瓶頸,理論上可以處理任意的變形問題;并且由于算法的Lagrangian特性,SPH離散粒子自然地追蹤流體物質(zhì)的運動及邊界的變化,不需要像Eluerian網(wǎng)格算法中采取復雜的界面追蹤技術,不存在Eluerian網(wǎng)格算法中因數(shù)值耗散造成界面追蹤精度的降低.由于SPH算法的以上特點,使得它特別適合模擬自由表面、運動界面及大變形問題[3?6].

SPH方法發(fā)展還尚未成熟,尤其是固壁邊界條件的施加方法,一直是困擾SPH方法發(fā)展的難點之一.早期Monaghan[7]模擬高速自由液面流動時引入了Lennard-Jones對勢的強排斥力,其中固壁邊界用虛粒子來描述,該模型中虛粒子對鄰近邊界的粒子強制施加排斥力,從而阻止鄰近邊界的粒子非物理穿透邊界.但是該方法動量不守恒,對鄰近邊界的粒子的描述與實際不符.Libersky等[8]和Randles等[9]首次提出鏡像粒子的方法來施加固壁邊界,在邊界之外布置關于邊界對稱的鏡像粒子,該粒子與真實粒子密度相同,速度相反,并參與計算,該方法守恒性好但鏡像粒子生成復雜,無法處理復雜壁面邊界的問題.Gotoh等[10]建議固壁邊界采用邊界粒子平衡內(nèi)部壓強來防止粒子穿透邊界,界面粒子采用多層鏡像粒子來實現(xiàn)滑移邊界條件,但仍難以處理復雜邊界問題.文獻[11,12]采用兩種類型的虛粒子來處理固定邊界條件,即一種虛粒子設置在固定邊界上,與Monaghan[7]所用的相似,另一種虛粒子分布在邊界的周圍,這與Libersky等[8]的鏡像粒子相似,即同時使用兩種虛粒子來阻止實粒子穿透邊界,但是這一方法也繼承了兩類虛粒子方法的不足.為了能夠處理復雜壁面邊界問題,近期國內(nèi)外學者把研究的重點放在如何更好地施加排斥力上,先后提出了幾種排斥力公式.Rogers等[13]提出了一種排斥力方法,該方法能很好地防止流體穿過壁面,但不能保證動量守恒;強洪夫等[14]基于Galerkin形式的SPH方程推導的罰函數(shù)邊界力能較好地保證動量守恒,但該方法含有未知參數(shù),需要經(jīng)過試算才能確定其取值,且邊界附近的流體粒子存在運動不穩(wěn)定問題;Monaghan等[15]和Liu等[16]也分別提出了各自的排斥力施加方法,這兩種方法類似,都能處理復雜邊界,但都需要將邊界粒子尺寸設為流體粒子尺寸的一半或更小,從而導致邊界粒子數(shù)目的增加,降低了計算效率,且這兩種方法不能保證動量守恒.

為了既保持排斥力方法能處理復雜邊界的優(yōu)勢,又能克服上述幾種方法的不足,本文提出一種新的排斥力模型來施加固壁邊界條件.分別通過靜止液柱算例和液柱坍塌算例,驗證了各種排斥力方法的動量守恒特性及邊界粒子尺寸對排斥力施加效果的影響;當僅采用邊界上單層粒子計算邊界力時存在邊界粒子缺失問題,本文采用了Liu等[16]提出的耦合邊界法中的虛粒子處理方法,以提高邊界附近流體粒子的計算精度,將該法應用到容器中液體的靜止算例及潰壩算例.對容器中液體的靜止算例,計算得到了容器底部壓強變化過程,通過與Monaghan等[15]與及強洪夫[14]所提出的方法進行對比,本文方法可有效克服邊界附近的流體粒子存在的運動不穩(wěn)定問題;通過潰壩算例驗證了本文方法處理復雜邊界問題的有效性.

2 SPH的基本方程及離散

2.1 控制方程

拉格朗日描述的流體動力學方程為:

參數(shù)P0可通過P0=100ρ0V/γ獲得,其中ρ0為流體初始密度,γ=7,Vmax為流體的最大速度.為保證流體的弱可壓縮性,一般取聲速c=10Vmax.

2.2 流體動力學方程的離散

在SPH方法中,連續(xù)的流場離散成為一系列相互作用的粒子,通過核函數(shù)估計技術在這些粒子上離散控制方程組,得到一組描述各粒子物理量隨時間變化的常微分方程組,即SPH基本方程組,再對這組方程采用相應的常微分方程組求解方法來推進時間進程的求解.

與耗散粒子動力學方法類似[17],求解域?被離散為一系列質(zhì)點,其中每一個質(zhì)點i記錄了流場在該點處的物理量,如密度ρi,速度vi,位置xi等,并且給予一個統(tǒng)計意義的質(zhì)量mi,每一個粒子同時擁有一個影響域S,影響域內(nèi)所包含的其他SPH粒子 j為該粒子的鄰近粒子,對于系統(tǒng)中任意點x處的物理量及其空間導數(shù)通過核函數(shù)估計得到,詳見文獻[3,4].

這樣,流體運動的控制方程(1)—(3)可離散為

2.3 時間積分

形成虛粒子,從而導致其缺乏通用性,很難推廣到復雜邊界區(qū)域;排斥力法一般僅需要在邊界處設置邊界粒子,避免了在計算區(qū)域外設置虛粒子的問題,因此具有很強的適應性,但需要構造合理的排斥力公式來保證邊界力的計算精度及穩(wěn)定性;耦合邊界法是排斥力法與虛粒子法的耦合,因此,合理的排斥力模型是耦合邊界法的基礎之一.本文著重研究以下5種典型的排斥力法,流體粒子與邊界粒子的相互作用如圖1所示.

本文采用leap-frog格式的時間積分方法,即

?表示密度ρ或速度v,xi是粒子i的位置坐標.為了使計算過程穩(wěn)定,本文采用Monaghan[7]給出的分別考慮具有黏性耗散和外力作用的時間步長表達式:

圖1 流體粒子與邊界粒子相互作用

方法1

Monaghan[7]最先提出用一組固定在邊界上的粒子對鄰近的流體粒子施加排斥力,從而防止流體粒子非物理穿透邊界,該排斥力類似于計算分子力時的Lennard-Jones方程,表達式如下:

其中f是作用在單位質(zhì)量上的力,ci代表當?shù)芈曀?

2.4 傳統(tǒng)排斥力模型

自由滑移邊界條件的定義如下:

為了提高計算的精度及穩(wěn)定性,國內(nèi)外學者提出了多種方法來施加壁面邊界條件.這些方法大致可分為三類:1)排斥力法;2)虛粒子法;3)耦合邊界法.其中虛粒子法需要將流體粒子映射到邊界外部

其中參數(shù)n1和n2一般取值為12和4,D一般取為最大速度的平方.r0為截止半徑,若其取值過大,會使鄰近邊界的流體粒子受到較大的邊界力,導致計算失敗;若取值較小,就很難阻止流體粒子穿透邊界.

方法2

Rogers等[13]提出的排斥力表達式為

其中nj為邊界粒子 j處的邊界法向,Ψ為流體粒子i到邊界的距離,ξ為xij在粒子 j處切線方向的投影,u⊥為粒子i的速度在nj上的投影.排斥力函數(shù)R(Ψ)的表達式為

其中?b為兩鄰近邊界粒子的距離,P(ξ)可保證粒子平行于邊界運動時,排斥力保持不變.ε(z,u⊥)是修正函數(shù),可根據(jù)流體的深度及速度調(diào)整邊界力的大小:

其中

式中h0為初始流體深度,z為流體局部深度,c0為初始聲速.

方法3

為克服Lennard-Jones方程計算邊界力時的不足,Monaghan等[15]又提出了一種邊界力公式:

式中K取值為gh0,β一般取4,mi為粒子i的質(zhì)量.

方法4

強洪夫等[14]基于Galerkin形式的SPH方程推導出了新的排斥力公式:

其中Aj是與邊界粒子 j有關的邊界權函數(shù),ε為罰參數(shù),其選取是個較困難的問題,一般需要若干次試算.由于該方法施加的邊界力大小與流體相對于邊界法向的速度大小有關,當流體粒子靠近邊界且相對速度很小時,會導致邊界力振蕩,從而導致靠近邊界的流體粒子速度、壓強和密度的不穩(wěn)定.

方法5

與Monaghan提出的邊界力公式相似,Liu等[16]也提出了一種排斥力方法:

上述5種方法各有優(yōu)缺點,下文將通過數(shù)值算例進行說明.

3 新型排斥力模型

如引言所述,國內(nèi)外學者已提出多種排斥力方法,但都存在各自的不足之處.為了克服上述問題,本文基于Galerkin加權余量法并結合上述排斥力方法,提出了一種新型排斥力公式(方法6)來施加壁面邊界條件.

Lagrangian描述的動量方程的張量形式如下:

vα是速度分量,σαβ是總應力張量分量,gα為單位質(zhì)量體積力分量,以上公式中的重復上標符合Einstein求和約定.應用Galerkin加權余量法,動量方程(27)式可寫為如下的弱形式:

式中右側第一項為內(nèi)力項,第二項為表面力項即邊界力項,第三項為體力項,其中?代表整個計算域空間,Γ為計算域邊界區(qū)域,nα為邊界的外法線方向分量,δα代表權函數(shù)(test function)分量,代表試函數(shù)(trial function)分量.

Galerkin方法的權函數(shù)及其梯度的表達式如下:

其中δ(x?xj)為Diracδ函數(shù).

方程(28)右側的第二項表示邊界力對動量方程的貢獻,應用罰函數(shù)法對該項進行變換.首先根據(jù)滑移接觸條件(14)來定義邊界的約束罰勢ΠB,其張量形式的表達式為

式中ε為罰參量,罰勢的一階變分形式為

將該罰勢的變分形式引入到動量方程的弱形式(28)式中,并代替邊界作用項,得到如下的形式:

將權函數(shù)和試函數(shù)的近似表達式(29),(30)和(31)代入(34)式,得到如下形式:

由于δvα的任意性,因此對于每個粒子必須滿足下式:

應用核函數(shù)的性質(zhì)并對上式進行點積分變化得到如下表達式:

通過將(39)式與(24)式對比可確定罰參數(shù)ε及Hij的取值.為使(39)式計算所得的邊界力與(24)式具有相同的數(shù)量級,本文取

當流體靠近壁面時,(39)式所得的邊界力與(vi?j)·nj成正比,為了避免方法4產(chǎn)生的壓力振蕩問題,本文對該項也進行了修正,最終得到的邊界排斥力的定義如下:

4 數(shù)值算例

本節(jié)通過4個數(shù)值算例來驗證本文方法的有效性,并與其他方法進行了對比.為了保證計算結果的可對比性,每種方法都采用相同的粒子配置,核函數(shù)及其影響域也保持一致.

分別通過靜止液柱算例及液柱坍塌算例,對比了傳統(tǒng)排斥力公式及本文新型排斥力公式的優(yōu)缺點,證明了本文方法的有效性.

對容器中液體靜止算例及潰壩算例,當固壁粒子僅采用邊界上的單層粒子分布時,會造成邊界附近粒子缺失,從而導致邊界附近的流體粒子計算精度降低,本文采用Liu等[16]提出的虛粒子施加方法來解決該問題.即在邊界粒子外側設置兩排虛粒子,邊界粒子的信息通過其鄰近的流體粒子插值得到,虛粒子的信息通過邊界粒子和流體粒子的信息插值獲得,如(41)—(43)式:

通過液體靜止算例,驗證了本文方法所得的邊界力與真實值的一致性;通過潰壩算例,驗證了本文方法處理復雜邊界問題的有效性.

4.1 靜止液柱算例

為了驗證邊界力對流體粒子初始分布的擾動特性,本節(jié)選取靜止液柱算例進行數(shù)值分析.

該算例中,液柱不受重力作用.其尺寸為0.4 m×0.6 m,邊界長度為1 m.核函數(shù)統(tǒng)一采用三次樣條函數(shù),粒子i的影響域hi=1.0?di,?di為粒子i的初始尺寸.初始密度ρ0=1000.0 kg/m3,邊界粒子與流體粒子尺寸大小相同,取為0.02 m,時間步長為1×10?5s,初始粒子分布如圖2,最底部的一層粒子為邊界粒子,其余粒子為流體粒子.由于流體靜止且不受外力作用,要保證系統(tǒng)動量守恒,在計算過程中流體粒子應始終保持靜止.

t=0.3 s時,采用不同方法模擬的結果如圖3,圖中箭頭方向為粒子的運動方向.從圖3可知,方法4及本文方法能保證流體粒子始終靜止.

圖2 靜止液柱算例初始SPH粒子分布(單位為m)

圖3 t=0.3 s時用不同方法模擬靜止液柱獲得的SPH粒子分布及速度分布(單位為m) (a)方法1;(b)方法2;(c)方法3;(d)方法4;(e)方法5;(f)方法6

方法1的模擬結果和截止半徑r0及影響域h的取值有關.若r0＜h/2,則流體粒子保持靜止,否則流體粒子將受到邊界力的擾動.本文取r0=1.1?d＞h/2,鄰近邊界的流體粒子由于受到邊界粒子排斥力的作用而產(chǎn)生運動;若增大r0的取值,流體粒子所受的排斥力增加,流體粒子的運動速度也會增大.選取四次樣條核函數(shù)及高斯核函數(shù)也存在同樣的問題.因此,截止半徑r0的取值問題不利于方法1在實際問題中的應用.

方法2、方法3及方法5中,無論選取何種核函數(shù),與邊界粒子鄰近的流體粒子受到的邊界排斥力均不為0,邊界力的作用使流體粒子產(chǎn)生運動,因此不能保證動量守恒;且方法2中邊界排斥力最大,流體粒子運動速度明顯大于其他方法所得的流體粒子運動速度.

方法4的罰參數(shù)取值為1.0.在方法4及方法6中,排斥力大小與流體和壁面之間的相對速度有關,僅當流體靠近壁面時,才受到邊界力的作用.這就使初始靜止流體能始終保持靜止,保證了鄰近邊界的流體粒子分布與實際相符.

若核函數(shù)統(tǒng)一采用高斯核函數(shù)或四次樣條函數(shù),影響域取初始粒子間距的1.0—1.5倍,計算結果與上述結果相似,所得的結論保持不變.

4.2 液柱坍塌算例

液柱在重力作用下會向下運動,排斥力的存在將阻止其穿透壁面,本算例對各種排斥力阻止流體粒子非物理穿透壁面的特性進行研究.

核函數(shù)及其影響域、液柱尺寸、初始設置與上節(jié)相同,但受重力作用,重力加速度g=9.81 m/s2,沿y軸負方向,

t=0.3 s時的計算結果如圖4所示.從圖4中可知:當邊界粒子和流體粒子尺寸大小相同時,方法1、方法3、方法4及方法5并不能阻止流體粒子穿透邊界,究其原因,是由于這4種方法的排斥力大小與邊界粒子和流體粒子之間的距離成反比.要防止流體粒子非物理穿透邊界,有兩種方法:一種是將邊界粒子尺寸設為流體粒子的1/2或更小,如文獻[15],但是這樣會增加邊界粒子數(shù)目,降低計算效率,尤其是三維問題;另一種是使排斥力大小與流體粒子到邊界粒子的法向距離成反比,如方法2和方法6,這種處理方法可減小邊界力精度對邊界粒子尺寸的依賴,提高計算效率.

方法4中,罰參數(shù)ε選取需要經(jīng)過若干次試算才能獲得.若罰參數(shù)過大,流體粒子在向下運動的過程中會受到過大的邊界力,從而導致該粒子加速遠離邊界運動,擾亂整個流場的粒子秩序,致使模擬失效;若罰參數(shù)過小,邊界排斥力很難阻止流體粒子穿透邊界.本算例經(jīng)過3次試算,分別取罰參數(shù)為0.1,1.0,10.0,最終確定罰參數(shù)取1.0.

4.3 容器中液體靜止算例

為了驗證本文方法所得的排斥力與邊界力的真實值是否一致,選取文獻[13]的容器中液體靜止算例進行檢驗.本節(jié)僅選取方法3、方法4與本文方法進行對比,來驗證本文方法的有效性.

圖4 t=0.3 s時,不同方法計算液柱坍塌得到的SPH粒子分布 (a)方法1;(b)方法2;(c)方法3;(d)方法4;(e)方法5;(f)方法6

容器長2 m,高1 m,水深H=0.9 m.方法3中,計算參數(shù)與文獻[15]一致,為防止流體粒子穿透邊界,邊界粒子尺寸設為流體粒子的1/2;方法4中,罰參數(shù)設為2.0,流體粒子尺寸設為0.02,邊界粒子尺寸設為0.015 m;本文方法中所有粒子的尺寸都設為0.02 m.3種方法的時間步長同取2.0×10?5s,所有流體粒子的初始壓強都設為0,初始密度ρ0=1000.0 kg/m3,重力加速度g=9.81 m/s2,沿y軸負方向

在重力作用下,流體粒子向下運動,但邊界力的存在阻礙了其進一步運動,這樣流體粒子就會在重力和邊界力的共同作用下振蕩運動,隨著時間的增加,流體粒子所受的力會達到平衡,容器底部的流體粒子的壓強理論值應為ρgH,H為液體深度.

t=1 s時,粒子分布狀態(tài)如圖5.方法3的計算結果如圖5(a)所示,可知高度越高,靠近邊界的流體粒子離邊界越遠,且該距離大于初始粒子間距,這就導致左上角及右上角處流體粒子的聚集及粒子秩序的紊亂.方法4及本文方法中,邊界力與粒子間的相對速度有關,從而克服了方法3存在的問題,鄰近邊界的流體粒子具有更好的分布秩序.

圖5 t=1.0 s時三種不同方法模擬的SPH粒子分布狀態(tài) (a)方法3;(b)方法4;(c)方法6

圖6 容器底部中央的SPH粒子相對壓強隨時間變化過程

以容器底部靠近中央的流體粒子為研究對象,記錄其壓強隨時間的變化過程如圖6.圖6中壓強為其相對值,即p/ρgH,液體趨于靜止時,相對壓強的理論值為1.0.從圖6可知,方法3及本文方法在t=0.6 s時壓強都達到了穩(wěn)定狀態(tài),相對壓強都趨于1.0.本文方法計算所得的相對壓強略小于1.0,主要是由人工黏性產(chǎn)生的耗散所導致的.采用方法4時,計算所得的相對壓強持續(xù)振蕩,t=0.6 s時壓強基本趨于穩(wěn)定,但t=1.02 s時,壓強又開始劇烈振蕩.其原因是:當鄰近邊界的流體粒子壓強趨于穩(wěn)定時,其運動速度趨于0,從(16)式可知,邊界力亦趨于0,流體粒子在壓力作用下會產(chǎn)生較大的加速度及速度向邊界靠近,此時在(16)式作用下又會產(chǎn)生較大的邊界力使流體粒子遠離邊界,這樣鄰近邊界的流體粒子就會產(chǎn)生振蕩的運動速度,由(5)式及(4)式可知,粒子運動速度的劇烈變化會導致粒子密度及壓強的大幅振蕩.本文方法則克服了方法4存在的缺陷.

圖7 采用本文方法得到的不同時刻的相對壓強分布 (a)t=0.8s;(b)t=1.0s;(c)t=1.2s

為了考察壓力場的空間分布特性,圖7給出了0.8,1.0,1.2s三個不同時刻的相對壓強分布狀態(tài).由圖7可知,本文方法計算所得的壓力場比較光滑;容器底部壓力場穩(wěn)定后的不同時刻,壓力場的空間分布一致且不隨時間變化.

4.4 潰壩算例

潰壩是經(jīng)典的自由表面流動問題,常用來驗證壁面邊界施加方法的有效性.計算模型如圖8,水柱位于容器左側,右側有一擋板,初始時刻水柱靜止,然后將擋板迅速抽出.水柱高度H=2L,寬度為L=0.146m,容器尺寸為2L×4L.

圖8 潰壩幾何模型

本文對潰壩過程進行了模擬,并與Koshizuka等[18]的實驗以及方法4、方法5所得的結果進行了對比.計算參數(shù)為:初始密度ρ0=1000.0 kg/m3,邊界粒子數(shù)目為805,邊界粒子初始尺寸為1.5mm,流體粒子數(shù)目為5000,流體粒子初始間距為 2.92mm,時間步長 1×10?5s,Vmax=方法4的罰參數(shù)經(jīng)試算后取0.5.核函數(shù)統(tǒng)一采用三次樣條函數(shù),粒子i的影響域hi=1.0?di.為克服邊界粒子缺失問題,本算例的3種方法都采用了Liu等[16]提出的虛粒子法來提高邊界附近的計算精度.

3種邊界施加方法所得的模擬結果如圖9,第1行為0.2,0.6s時的實驗結果;第2行為方法4所得的結果;第3行為方法4所得的結果;第4行為采用本文方法所得的模擬結果.圖9中箭頭方向表示粒子的運動方向,其長短表示粒子運動速度的大小.

從圖9中可知,3種方法所得的自由表面運動過程與實驗結果都比較符合.但是,采用方法4施加自由滑移邊界條件時,靠近邊界處部分流體粒子的速度明顯大于其周圍區(qū)域流體粒子的速度;隨著時間的增加,這種非物理的速度振蕩并不會消失,且邊界流體粒子速度失真會導致粒子密度和壓強的不穩(wěn)定,最終導致邊界附近流場的數(shù)值模擬失真.采用方法5及本文方法時,靠近邊界的流體粒子始終保持平行壁面邊界運動,粒子速度場分布光滑且秩序良好,很好地克服了方法4所導致的速度振蕩問題.

5 結論

本文提出一種新型排斥力模型,該模型不需要確定未知參數(shù),可較好地防止流體粒子非物理穿透壁面,能克服流體相對壁面速度較小時速度不穩(wěn)定問題.通過液柱靜止算例、液柱坍塌算例及容器中液體靜止算例對比并驗證了本文方法的動量守恒性、防止流體粒子非物理穿透邊界的特性及邊界力的精確度;最后,將本文方法應用到潰壩算例,進一步驗證了該方法適應復雜邊界問題的能力.

本文方法可處理復雜形狀邊界,不需要確定未知參數(shù),可較好地施加固壁邊界條件.本文方法還得出:邊界排斥力的大小要與流體粒子到邊界的法向距離成反比,而不是與流體粒子到邊界粒子的距離成反比,此時才能保證邊界力不受邊界粒子尺寸的影響,有效防止流體粒子非物理穿透壁面邊界;排斥力模型中需要考慮流體與壁面的相對速度對排斥力的影響,否則會導致邊界附近流體粒子速度失真.

圖9 三種不同邊界施加方法所得的潰壩速度場及與實驗結果的對比 (a)實驗結果,t=0.2 s;(b)實驗結果,t=0.6 s;(c)方法4,t=0.2 s;(d)方法4,t=0.6 s;(e)方法5,t=0.6 s;(g)本文方法,t=0.2 s;(h)本文方法,t=0.6 s

[1]Lucy L B 1977 Astron.J.82 1013

[2]Gingold R A,Monaghan J J 1977 Mon.Not.R.Astron.Soc.181 375

[3]Zhang A M 2008 Chin.Phys.B 17 927

[4]Sun Z H,Han R J 2008 Chin.Phys.B 17 3185

[5]Zhang A M,Yao X L 2008 Acta Phys.Sin.57 339(in Chinese)[張阿漫,姚熊亮2008物理學報57 339]

[6]Monaghan J J 2005 Rep.Prog.Phys.68 1703

[7]Monaghan JJ 1994 J.Comput.Phys.110 399

[8]Libersky L D,Petscheck A G,Carney T C 1993 J.Comput.Phys.109 67

[9]Randles P W,Libersky L D 1996 Comput.Methods Appl.Mech.Eng.138 375

[10]Gotoh H,Sakai T 1999 Coast.Eng.41 303

[11]Liu G R,Gu Y T 2001 Struct.Eng.Mech.246 29

[12]Liu M B,Chang JZ 2010 Acta Phys.Sin.59 3654(in Chinese)[劉謀斌,常建忠2010物理學報59 3654]

[13]Rogers B,Dalrymple R 2007 Adv.Num.Model Simul.Tsun.Wave Runup.10 75

[14]Qiang H F,Han Y W,Wang K P,Gao W R 2011 Eng.Mech.28 245(in Chinese)[強洪夫,韓亞偉,王坤鵬,高巍然2011工程力學28 245]

[15]Monaghan J J,Kajtar J B 2009 Comput.Phys.Commun.180 1811

[16]Liu M B,Shao J R 2011 Sci.China:Technol.Sci.10 1

[17]Chang J Z,Liu M B,Liu H T 2008 Acta Phys.Sin.57 3954(in Chinese)[常建忠,劉謀斌,劉漢濤2008物理學報57 3954]

[18]Koshizuka S,Oka Y 1996 Nucl.Sci.Eng.123 421