不同質(zhì)量粒子分布SPH方法及其應(yīng)用

不同質(zhì)量粒子分布SPH方法及其應(yīng)用

艾孜海爾·哈力克1, 2，熱合買(mǎi)提江·依明2，開(kāi)依沙爾·熱合曼2，買(mǎi)買(mǎi)提明·艾尼1

(1. 新疆大學(xué)機(jī)械工程學(xué)院，烏魯木齊830047;2. 新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，烏魯木齊830046)

摘要：光滑粒子流體動(dòng)力學(xué)(SPH)方法是一種新發(fā)展起來(lái)的純拉格朗日數(shù)值方法，其最大的優(yōu)點(diǎn)是不需要對(duì)自由表面進(jìn)行特殊處理，有利于模擬飛濺、破碎等高度非線性流。SPH方法的一個(gè)缺點(diǎn)是計(jì)算的時(shí)間復(fù)雜性高，因?yàn)樾枰罅苛Ｗ釉敿?xì)定義流體，使之計(jì)算時(shí)間相應(yīng)增加。在核半徑保持不變的情況下為了提高計(jì)算效率和精度提出不同質(zhì)量粒子分布的方法，即在關(guān)鍵計(jì)算區(qū)域分布較小質(zhì)量粒子，而在其它區(qū)域分布較大質(zhì)量粒子。首先通過(guò)靜水?dāng)?shù)值模擬，初步驗(yàn)證了此方法的有效性。之后，對(duì)半浸泡振蕩圓筒所產(chǎn)生的小振幅水波進(jìn)行了數(shù)值模擬。將數(shù)值模擬結(jié)果與已有的實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。分析結(jié)果表明此方法能夠明顯提高計(jì)算精度的同時(shí)可以節(jié)省大量計(jì)算時(shí)間。

關(guān)鍵詞：光滑流體動(dòng)力學(xué)方法；波浪；振蕩圓筒；不同質(zhì)量粒子；流體與結(jié)構(gòu)物相互作用

中圖分類(lèi)號(hào):TP391文獻(xiàn)標(biāo)志碼：A

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助(50975036，51275064)；中央高校基本科研業(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金(3132013063)；遼寧省自然科學(xué)基金(2015020136)

收稿日期：2014-11-21修改稿收到日期：2015-04-27

SPH method with different mass particle distributions and its application

AZHARHalik1,2,RAHMATJANImin2,KAYSARRahman2,MAMTIMINGeni1(1. School of Mechanical Engineering, Xinjiang University, Urumqi 830047, China; 2. Schools of Mathematics and System Science, Xinjiang University, Urumqi 830046, China)

Abstract:The smoothed particle hydrodynamic (SPH) method is a newly developed pure Lagrangian numerical method, its major advantage is that no special treatment of a free surface is required, it is beneficial for simulating highly non-linear flows to run up and splash around bodies. Its one drawback is that it has a high computational time complexity associated with a large number of particles desirable for good flow definition in order to improve the computation accuracy without substantially increasing computational time. Here, a method with different mass particle distributions was used for a pre-selected area where the high resolution was desirable. The effectiveness of the proposed method was tested and verified with a static water problem. Then, the surface waves generated by an oscillating cylinder with different wave periods and strokes were numerically simulated, and the simulation results were compared with the experimental data. The results showed that using this method brings a significant improvement in the computation accuracy without substantially increasing computational time.

Key words:smoothed particle hydrodynamic (SPH) method; wave; oscillating cylinder; different mass particle

通過(guò)歐拉計(jì)算流體動(dòng)力學(xué)方法對(duì)自由表面流體問(wèn)題進(jìn)行處理是比較困難的，主要是因?yàn)榱黧w不僅復(fù)雜多變，而且運(yùn)動(dòng)時(shí)所受力的數(shù)學(xué)模型多為空間和時(shí)間函數(shù)的偏微分方程，編程實(shí)現(xiàn)難度較大。SPH方法是一種純拉格朗日數(shù)值方法。SPH方法首先由Lucy等[1-2]分別提出并應(yīng)用于天體物理學(xué)中。此方法中每個(gè)粒子卷帶了質(zhì)量、位置、速度、內(nèi)能等各種物理量。SPH方法最大的優(yōu)點(diǎn)是，不需要對(duì)自由表面進(jìn)行特殊處理，有利于模擬飛濺、破碎等高度非線性流。Monaghan[3]通過(guò)SPH方法模擬了剛體與水的相互作用，Hopkins等[4]模擬了流體混合問(wèn)題。國(guó)內(nèi)則有強(qiáng)洪夫等[5]通過(guò)大密度差多相流 SPH 方法模擬了二維液滴碰撞，劉棟等[6]模擬了液滴變形及表面張力。

SPH方法的一個(gè)缺點(diǎn)是計(jì)算效率低。通過(guò)CPU和GPU并行計(jì)算可以一定程度提高計(jì)算效率[7]，但還是需要大量粒子詳細(xì)定義流體，使計(jì)算時(shí)間隨之增加。為了減少粒子數(shù)量，先前的SPH方法采用的是可變光滑長(zhǎng)度法[8]或粒子動(dòng)態(tài)分布法[9]。可變光滑長(zhǎng)度法需要更多的計(jì)算時(shí)間來(lái)處理額外的時(shí)間和空間函數(shù)，而粒子動(dòng)態(tài)分布法則存在守恒問(wèn)題[10]。為了提高計(jì)算效率和精度,本文提出不同質(zhì)量粒子分布的思想，即在關(guān)鍵計(jì)算區(qū)域分布較小質(zhì)量粒子，而在其它區(qū)域分布較大質(zhì)量粒子。本文首先通過(guò)靜水問(wèn)題的數(shù)值模擬，初步驗(yàn)證了此方法的有效性。之后，對(duì)半浸泡振蕩圓筒產(chǎn)生的小振幅水波進(jìn)行了數(shù)值模擬。將SPH數(shù)值模擬結(jié)果與文獻(xiàn)[11]中的實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。最后通過(guò)不同質(zhì)量粒子分布法對(duì)各種波幅進(jìn)行了詳細(xì)比較，驗(yàn)證了此方法的有效性。

1SPH方法基本原理

1.1SPH方法

SPH方法中，在插值理論的基礎(chǔ)上, 首先采用核函數(shù)近似方法將函數(shù)及其梯度用核函數(shù)的積分表示，把偏微分形式的連續(xù)體方程轉(zhuǎn)化為積分形式的方程，隨后采用粒子近似方法將連續(xù)形式的積分方程轉(zhuǎn)化成離散形式的方程[12]。位于ri處的粒子i 的場(chǎng)函數(shù)φ的SPH離散形式可寫(xiě)為

(1)

式中:WPs表示流體粒子，φ可為場(chǎng)域中的某個(gè)矢量或標(biāo)量，如速度或密度等，ωj為j粒子在rj處的體積，Wi j=W(ri j,h)為核函數(shù)或被稱為光滑函數(shù)，h為核半徑，ri j=|ri-rj|。用分部積分、高斯定理和核函數(shù)性質(zhì)，梯度▽?duì)盏腟PH離散形式為

(2)

1.2狀態(tài)方程

在SPH方法中不可壓縮流體是當(dāng)做弱可壓縮流體來(lái)處理的[13]。以下?tīng)顟B(tài)方程用來(lái)計(jì)算流體壓力。

(3)

1.3核函數(shù)

由于三次樣條函數(shù)是模擬水槽中水流演進(jìn)過(guò)程的最佳選擇[14]，本文選擇它作為核函數(shù)。

1.4邊界條件

SPH方法中常用的邊界條件有虛粒子法[15]、固定粒子法[16]、排斥力法[17]以及動(dòng)態(tài)邊界條件法[18]等。這些方法都有相應(yīng)的優(yōu)缺點(diǎn)，如虛粒子法在邊界有尖角或高曲率曲面時(shí)比較難以實(shí)現(xiàn)，固定邊界法則可能產(chǎn)生較大的非物理邊界層。本文選擇排斥力邊界條件法，此方法比較容易實(shí)現(xiàn)復(fù)雜邊界以及粒子間的相互作用。排斥力邊界條件法由Monaghan提出[17]。根據(jù)此方法Rogers和Dalrymple等[19]對(duì)垂直于固壁邊界移動(dòng)的流體粒子的受力計(jì)算進(jìn)行了修正[19]。

1.5時(shí)間積分

本文開(kāi)發(fā)模擬程序采用的是預(yù)測(cè)校正時(shí)間積分法，該積分法預(yù)測(cè)的時(shí)間演化公式為

(5)

式中:n為時(shí)間步，對(duì)這些值在半個(gè)時(shí)間步時(shí)應(yīng)力作用下做如下修正

(6)

時(shí)間步結(jié)束時(shí)再做如下修正

(7)

2數(shù)值模擬

2.1不同質(zhì)量粒子分布法

為了減少粒子數(shù)量，先前的SPH方法采用的是可變光滑長(zhǎng)度法[8]或粒子動(dòng)態(tài)分布法[9]?？勺児饣L(zhǎng)度法需要更多計(jì)算時(shí)間來(lái)處理額外的時(shí)間和空間函數(shù)，而粒子動(dòng)態(tài)分布法則存在守恒問(wèn)題[10]。本文采用的不同質(zhì)量粒子分布法中核半徑保持不變，只在關(guān)鍵區(qū)域內(nèi)質(zhì)量為m的粒子由4個(gè)質(zhì)量為m/4的粒子所取代細(xì)化，見(jiàn)圖1。

圖1　質(zhì)量為m的粒子由4個(gè)質(zhì)量為m/4的粒子所取代 Fig.1 Replacing a particle with a mass m with four lighter ones with a mass m/4

本文模擬案例中選取核半徑h=1.3Δx,Δx為較重粒子的初始間距。此長(zhǎng)度被證明具有很好的魯棒性[20]。由于較輕粒子保持與較較重粒子相同的核半徑，因此兩類(lèi)粒子核半徑內(nèi)都有可能包含不同質(zhì)量粒子。較輕粒子相互作用的粒子數(shù)會(huì)更多。只在所選區(qū)域分布較輕粒子可以在不降低精確度的條件下節(jié)省大量計(jì)算時(shí)間。

2.2不同質(zhì)量分布法在靜水模擬中的應(yīng)用

在將此方法運(yùn)用于動(dòng)態(tài)問(wèn)題之前，首先在靜水條件下對(duì)其準(zhǔn)確性進(jìn)行了初步評(píng)估。模擬的初始設(shè)置如圖2所示，水槽寬度為2 m，水深為1 m，將計(jì)算域分解成兩個(gè)區(qū)域，在區(qū)域1中分布較輕粒子(Δx1=0.01 m)，在區(qū)域2中分布較重粒子(Δx2=0.02 m)較輕粒子分布對(duì)稱于中心垂軸。

圖2　不同質(zhì)量粒子分布區(qū)域 Fig.2 The regional distribution of particles with different mass

圖3給出了t=5 s時(shí)靜水壓力的精確值與SPH計(jì)算結(jié)果的比較，可以看出兩個(gè)結(jié)果較吻合。在兩種粒子交界處產(chǎn)生了很小的誤差，這是由于運(yùn)用三次樣條函數(shù)時(shí)不同質(zhì)量粒子核半徑保持不變，因此在不同質(zhì)量粒子交界面可能產(chǎn)生拉伸失穩(wěn)現(xiàn)象。Monaghan[21]提出，可以通過(guò)壓力校正避免此類(lèi)問(wèn)題。在z→0處由于核函數(shù)支持域缺損而產(chǎn)生一定的誤差。圖4給出了t=5 s時(shí)兩種粒子分布圖，可以看出此時(shí)沒(méi)有出現(xiàn)明顯的粒子聚集現(xiàn)象，粒子分布關(guān)于中心軸對(duì)稱，不同粒子交界面沒(méi)有出現(xiàn)非物理混合現(xiàn)象，保持了良好的靜水狀態(tài)。

圖5(b)給出了t=10 s時(shí)的粒子分布圖。可以看出此時(shí)也沒(méi)有出現(xiàn)明顯的粒子聚集現(xiàn)象，粒子分布關(guān)于中心軸對(duì)稱，不同粒子交界面沒(méi)有出現(xiàn)非物理混合現(xiàn)象，保持了良好的靜水狀態(tài)。在自由表面附近由于核函數(shù)支持域缺損，出現(xiàn)了很小的移動(dòng)或變形。

圖3　t=5 s時(shí)靜水壓力的精確值與SPH計(jì)算結(jié)果的比較 Fig.3 Comparisons between the exact and SPH hydrostatic pressure at t=5 s

圖4　t=5 s時(shí)兩種粒子分布圖 Fig.4 Distribution of two kinds of particles at t=5 s

圖6給出的模擬中較輕粒子分布于斜角區(qū)域，粒子質(zhì)量比同樣為1∶4?？梢钥闯隽Ｗ颖３衷?，沒(méi)有出現(xiàn)不規(guī)則交界面引起的移動(dòng)或混合現(xiàn)象。以上兩個(gè)數(shù)值模擬中的時(shí)間步(Courant number,Cr=0.2)均為原時(shí)間步(Cr=0.4)的一半。

為檢查質(zhì)量比的極限，圖7中給出的模擬中分布的粒子質(zhì)量比為1∶16?？梢钥闯鰐=10 s時(shí)，由于較輕粒子初始間距很小，使之核梯度遠(yuǎn)遠(yuǎn)小于之前案例，從而在不同粒子交界區(qū)域出現(xiàn)了明顯的拉伸失穩(wěn)和粒子聚集現(xiàn)象。通過(guò)引入較大的人工壓力修正此失穩(wěn)現(xiàn)象時(shí)，可能會(huì)引起靜水壓力誤差明顯增大等其他數(shù)值問(wèn)題。

圖5 粒子質(zhì)量比為1∶4時(shí)SPH計(jì)算結(jié)果Fig.5SPHresultsforthemassratioof1∶4圖6 粒子分布于斜角區(qū)域,粒子質(zhì)量比為1∶4Fig.6Diagonallyslanteddistributionofparticlesforthemassratioof1∶4圖7 粒子質(zhì)量比為1∶16時(shí)SPH計(jì)算結(jié)果Fig.7Particledistributionforthemassratioof1∶16

為了評(píng)估方法在靜水案例中的收斂性，可以定義如下全局相對(duì)誤差[20]

(8)

式中：n表示時(shí)間步，xi為粒子i在x方向上的坐標(biāo)。本模擬中xi始終大于零，因此避免了xi→0時(shí)Rx的雜散值。圖8給出了圖6中的斜角分布粒子的橫坐標(biāo)誤差。當(dāng)Rx<4×10-5時(shí)，可以認(rèn)為模擬達(dá)到穩(wěn)定狀態(tài)[22]。通過(guò)上述案例可以看出，粒子質(zhì)量比為1∶4時(shí)，不同質(zhì)量粒子分布法可以成功運(yùn)用到靜水案例中。

圖8　斜角分布粒子的橫坐標(biāo)誤差 Fig.8 Error of the horizontal particle positions for diagonally slanted mass distribution

2.3SPH方法模擬振蕩圓筒產(chǎn)生的小幅波

Yu等[11]測(cè)量了由振蕩圓筒產(chǎn)生的表面波,提出了無(wú)限寬水槽中半浸泡振蕩圓筒產(chǎn)生的水波問(wèn)題的線性理論，得到了表面波運(yùn)動(dòng)的理論公式,并在寬度為30m的水槽中進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果與理論結(jié)果相當(dāng)吻合。

圖9給出了振蕩圓筒SPH數(shù)值模擬的初始設(shè)置，水槽寬度為12 m，水深為0.58 m，振蕩圓筒置于中點(diǎn)位置。水波被兩側(cè)寬度為1.5 m的海綿邊界層所吸收衰減。一開(kāi)始圓筒為半浸泡狀態(tài)，之后開(kāi)始縱向振動(dòng)，作小振幅簡(jiǎn)諧運(yùn)動(dòng)，圓筒的震蕩振幅為1.23 cm，震蕩周期為0.76 s。參加運(yùn)算粒子數(shù)為18 500個(gè)，粒子初始間距和圓筒直徑分別為Δx=0.02 m，D=0.15 m。圖10給出了兩個(gè)不同時(shí)刻的粒子分布圖。

圖9　振蕩圓筒SPH數(shù)值模擬的初始設(shè)置 Fig.9 Schematic of SPH simulation of a heaving cylinder

圖10　振蕩圓筒SPH模擬壓力分布圖 Fig.10 Pressure distributions of SPH simulation of a heaving cylinder

圖11給出了t=10 s時(shí)的表面波形理論結(jié)果和SPH模擬結(jié)果的比較。兩種結(jié)果較吻合。與實(shí)驗(yàn)結(jié)果進(jìn)行比較時(shí)采用了[11]中提出的通過(guò)波幅比進(jìn)行比較的方法。將距圓筒三個(gè)波長(zhǎng)遠(yuǎn)處的波幅比與公式(9)中的無(wú)限遠(yuǎn)處理論波幅比RA進(jìn)行比較。SPH模擬波幅比為RA,SPH=0.58，這與實(shí)驗(yàn)值RA,實(shí)驗(yàn)值=0.543較接近。

(9)

圖11　t=10 s時(shí)的表面波形理論 結(jié)果和SPH模擬結(jié)果的比較 Fig.11 Comparison of surface profile between analytical and SPH results at t=10 s

表1給出了圓筒在不同震蕩周期和振幅條件下的SPH波幅比和實(shí)驗(yàn)波幅比，可以看出結(jié)果相當(dāng)吻合。粒子細(xì)化分布(Δx=0.01 m)時(shí)，可以得到更好的結(jié)果(RA,細(xì)化分布=0.54)，但是粒子數(shù)大量增多(71 000個(gè))，使計(jì)算時(shí)間大幅增加。

表1　圓筒在不同震蕩周期和

2.4不同質(zhì)量粒子分步法模擬震蕩圓筒

圖12給出了粒子的初始分布，由于圓筒周?chē)鷧^(qū)域?yàn)槟M中的關(guān)鍵區(qū)域，因此在此區(qū)域分布較輕粒子(Δx1=0.01 m)，其它區(qū)域分布較重粒子(Δx2=0.02 m) ，圓筒震蕩周期為0.76 s，振幅為1.23 cm。圖13給出了t=10 s時(shí)表面波形的理論結(jié)果和不同質(zhì)量粒子分布SPH法模擬結(jié)果的比較，可以看出結(jié)果較吻合。

圖12　不同質(zhì)量粒子初始分布圖 Fig.12 Initial distribution of different mass particles

數(shù)值模擬的波幅比為RA,SPH=0.53，這與實(shí)驗(yàn)結(jié)果相當(dāng)吻合。此模擬中用到的粒子數(shù)為22 300個(gè)，這比圖10中的18 500個(gè)粒子只多出了3 800個(gè)粒子，但是比粒子細(xì)化分布(Δx=0.01 m)中的71 000個(gè)粒子少了48 700個(gè)粒子。細(xì)化粒子分布單機(jī)運(yùn)算時(shí)間為40 h。不同質(zhì)量粒子分布法的單機(jī)運(yùn)算時(shí)間為22 h。粗略粒子分布(Δx=0.02 m)的運(yùn)算時(shí)間為16 h。可以看出此方法在計(jì)算效率和精確度方面均有所優(yōu)勢(shì)。

圖13　t=10 s時(shí)的表面波形理論結(jié)果和 不同質(zhì)量粒子分布SPH法模擬結(jié)果的比較 Fig.13 Comparison of surface profile between analytical and different mass distribution SPH results at t=10 s

圖14給出了幾種方法通過(guò)公式(10)計(jì)算的誤差率，可以看出不同質(zhì)量粒子分布法在保證粒子數(shù)沒(méi)有大幅增多的條件下達(dá)到了接近于細(xì)化粒子分布時(shí)的誤差率。

(10)

圖14　4種不同的粒子分布測(cè)試中誤差率的比較 Fig.14 Comparison of 4 convergence tests with different kinds of distribution

表2給出了三種不同的粒子分布條件下模擬振幅為1.23 cm，周期為0.76 s的振蕩圓筒案例時(shí)的粒子數(shù)、波幅比及計(jì)算時(shí)間。可以看出不同質(zhì)量粒子分步法在提高了精度的同時(shí)節(jié)省了大量計(jì)算時(shí)間。

表2　采用3種不同的粒子分布計(jì)算結(jié)果的比較

表3給出了不同振幅和周期的振蕩圓筒案例通過(guò)三種粒子分布所得到的波幅比的數(shù)值結(jié)果?？梢钥床煌|(zhì)量粒子分步法在控制粒子數(shù)的前提下達(dá)到了接近于細(xì)化粒子分布時(shí)的計(jì)算精度。

表3　不同振幅和周期的振蕩圓筒案例通過(guò)三種粒子分布所得到的波幅比的數(shù)值結(jié)果

3結(jié)論

為了提高SPH數(shù)值模擬方法的計(jì)算效率和精確度，本文提出了不同質(zhì)量粒子分布法，即在關(guān)鍵計(jì)算區(qū)域分布較小質(zhì)量粒子，其他區(qū)域分布較大質(zhì)量粒子，從而控制參加計(jì)算的粒子數(shù)量，提高計(jì)算效率。此方法既避免了細(xì)化分布法中粒子數(shù)過(guò)多的問(wèn)題，也可使核半徑保持不變。首先通過(guò)靜水模擬對(duì)此方法進(jìn)行了初步評(píng)估。之后，對(duì)半浸泡振蕩圓筒產(chǎn)生的小振幅水波進(jìn)行了數(shù)值模擬。將SPH數(shù)值模擬結(jié)果與Ursell和Yu的實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。分析結(jié)果表明此方法能夠明顯提高計(jì)算精度的同時(shí)可以節(jié)省大量計(jì)算時(shí)間。今后將進(jìn)一步研究通過(guò)Monaghan提出的人工壓力引入法修正不同質(zhì)量粒子交界面可能產(chǎn)生拉伸失穩(wěn)現(xiàn)象。

參考文獻(xiàn)

[1]Lucy L B. A numerical approach to the testing of fusionprocess [J]. Astronomical Journal,1977, 88:1013-1024.

[2]Gingold R A, Monaghan J J. Smoothed particle hydrodynamics [J]. Monthly Notices of the Royal Astronomical Society,1977, 235:911-934.

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

[4]Hopkins P F. A general class of Lagrangian smoothed particle hydrodynamics methods and implications for fluid mixing problems[J]. Monthly Notices of the Royal Astronomical Society, 2013, 428(4): 2840-2856.

[5]強(qiáng)洪夫, 石超, 陳福振, 等. 基于大密度差多相流 SPH 方法的二維液滴碰撞數(shù)值模擬[J].物理學(xué)報(bào),2013,62(21): 214701-214701.

QIANG Hong-fu,SHI Chao,CHEN Fu-zhen,et al.Simulation of two-dimensional droplet collisions based on SPH method of multi-phase flows with large density differences [J]. Acta Phys. Sin,2013, 62(21): 214701-214701.

[6]劉棟, 郭印誠(chéng), 林文漪. 液滴變形及表面張力模擬的光滑粒子動(dòng)力學(xué)方法[J]. 清華大學(xué)學(xué)報(bào):自然科學(xué)版,2013 (3):384-388.

LIU Dong, GUO Yin-cheng, LIN Wen-yi. Droplet deformation and surface tension modeling using the smoothed particle hydrodynamics method[J]. J Tsinghua Univ.:Sci.& Tech.,2013(3):384-388.

[7]Domínguez J M, Crespo A J C, Gómez-Gesteira M. Optimization strategies for CPU and GPU implementations of a smoothed particle hydrodynamics method[J]. Computer Physics Communications, 2013, 184(3): 617-627.

[8]Imin R, Geni M. Stress analysis of gear meshing impact based on SPH method[J]. Mathematical Problems in Engineering, 2014,2014:328216.1-328216.7.

[9]Barcarolo D A, Le Touzé D, Oger G, et al. Adaptive particle refinement and derefinement applied to the smoothed particle hydrodynamics method[J]. Journal of Computational Physics, 2014:273(5):640-657.

[10]Feldman J, Bonet J. Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems[J]. International Journal for Numerical Methods in Engineering,2007, 72:295-324.

[11]Yu Y S, Ursell F. Surface wave generated by an oscillating circular cylinder on water of finite depth: theory and experiment[J]. Journal of Fluid Mechanics,1961,11:529-551.

[12]Monaghan J J. Smoothed particle hydrodynamics [J]. Reports on Progress in Physics, 2005, 68:1703-1759.

[13]Dalrymple R A, Rogers B D. Numerical modeling of water waves with the SPH method. Coastal Engineering, 2006, 53:141-147.

[14]Guilcher P M, Ducorzet G, Alessandrini B, et al. Water wave propagation using SPH models[J]. Proceedings of Second International SPHERIC Workshop, Spain, 2007:119-124.

[15]Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics[J]. Journal of Computational Physics,2003,191:448-475.

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

[17]Monaghan J J, Kos A. Solitary waves on a Cretan Beach[J]. Journal of Waterway, Port, Coastal and Ocean Engineering,1999, 125:145-154.

[18]Crespo A J C,Gomez-Gesteira M, Dalrymple R A, Boundary conditions generated by dynamic particles in SPH methods. Computers, Materials & Continua,2007,5(3):173-184.

[19]Rogers B D, Dalrymple R A. SPH modelling of tsunami waves[J]. In Advances in Coastal and Ocean Engineering, 2008,10: 75-100.

[20]Quinlan N J, Basa M, Lastiwka M. Truncation error in mesh-free particle methods[J]. International Journal for Numerical Methods in Engineering,2006, 66:2064-2085.

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

[22]Zhou J, Causon D, Mingham C, et al. The surface gradient method for the treatment of source terms in the shallow water equations[J]. Journal of Computational Physics,2001,168:1-25.

第一作者夏冬生男，博士，副教授，1973年12月生

通信作者張會(huì)臣男，博士，教授，1965年10月生