光滑粒子水動力學在水利與海洋工程中的應用研究進展

林鵬智，劉 鑫

(1.四川大學水力學與山區(qū)河流開發(fā)保護國家重點實驗室，四川成都 610065;2.河海大學港口海岸與近海工程學院，江蘇南京 210098)

1 光滑粒子水動力學簡介

近年來，隨著計算機硬件技術(shù)的不斷進步，各種新的數(shù)值模擬方法逐漸被科研和工程領域重視并得到了廣泛應用，其中光滑粒子水動力學(smoothed particle hydrodynamics，SPH)是最具代表性的一種。不同于傳統(tǒng)的基于歐拉網(wǎng)格的數(shù)值計算方法，SPH基于拉格朗日系統(tǒng)建立控制方程，在整個模擬過程可以完全不依賴計算網(wǎng)格［1］，這使得SPH具有一些網(wǎng)格法所不具有的優(yōu)勢，如適于模擬具有復雜或者大變形自由面的水流問題和處理具有復雜流固交界面甚至是運動結(jié)構(gòu)物的流固耦合問題。自SPH在1977年被兩位學者［2-3］獨立提出開始，其應用范圍已經(jīng)從最初的天體物理學發(fā)展到目前的流體動力學問題。

1.1 拉格朗日描述的控制方程

在計算流體運動時，SPH通過求解 Navier-Stokes方程(N-S方程)或它的其他演化形式(如水深平均的淺水方程等)解決問題。拉格朗日描述下的一般形式的N-S方程如下:

式中:u為速度;p為壓強;g為重力加速度;ν0為分子運動黏性系數(shù);ρ為密度。相比傳統(tǒng)網(wǎng)格法中采用歐拉描述的控制方程，方程(2)中不存在對流項，網(wǎng)格之間的對流效應信息由拉格朗日粒子的自身速度攜帶。

在對方程(1)(2)進行數(shù)值求解時，SPH數(shù)值模型可細分為弱可壓縮模型(weakly compressible SPH，WCSPH)和不可壓縮模型(incompressible SPH，ISPH)兩大類。WCSPH模型中，方程(1)包含流體密度變化項，輔以虛擬的流體狀態(tài)方程直接將密度變化和壓強建立聯(lián)系，可以在不求解矩陣的情況下顯式獲得壓強。但為了使時間步長不至于太小，往往會采用較小的人工彈性模量，在計算高速流動時可能產(chǎn)生較大的誤差。ISPH模型則采用不可壓假設，將方程(1)簡化為速度散度為零的形式，即

然后再通過求解壓力Poisson方程(pressure Poisson equation，PPE)獲得壓強，進而獲得速度場。

一般講，ISPH模型能夠采用更大的時間步長進行計算，得到比WCSPH模型更加光滑穩(wěn)定的壓力場，在模擬劇烈的自由面變化等問題時壓力場更加穩(wěn)定，但對邊界條件的處理要求也更加嚴格［4］。WCSPH模型則更適合并行計算的要求［5］，為了獲取更加穩(wěn)定的壓力場，一般可以通過引入人工黏性等數(shù)值處理方法。當前學術(shù)界對二者的研究成果都很豐富。

1.2 計算域的SPH離散

根據(jù)狄拉克函數(shù)的數(shù)學性質(zhì)，計算域中任意點x0處的變量(如密度、速度、壓強)可以精確表示成其周圍點的積分形式:

式(4)被稱為場變量的積分表示，而δ即δ函數(shù)，或稱狄拉克函數(shù)，Ω為積分域。在SPH模擬中，使用離散的拉格朗日點離散上述連續(xù)積分，并用多項式函數(shù)W代替超越函數(shù)δ，將上述積分寫成離散形式，得到積分表示的粒子離散形式(圖1):

式中:i為當前參考點的變量值;j為積分域內(nèi)離散的拉格朗日點的變量;W為核函數(shù)，也有學者稱其為權(quán)函數(shù)或者光滑函數(shù)，是用以代替狄拉克函數(shù)的多項式函數(shù);h為光滑長度，κh為搜索域半徑，系數(shù)κ根據(jù)模型不同而不同，通常取2～3;V為積分自變量，在二維情況下是圓的面積，而在三維情況下則是球的體積。關于粒子離散的詳細介紹可見文獻［6］，此處不再贅述。具體到水動力學問題中，拉格朗日點成為具有一定質(zhì)量、體積、速度、壓強等屬性，并攜帶位置信息的粒子，而j則表示搜索域Ω內(nèi)i的所有相鄰粒子。這時，根據(jù)式(5)，粒子i的SPH表達式為

圖1 SPH方法中的粒子離散

式中m為粒子質(zhì)量。

總體來說，與一般的網(wǎng)格法相比，SPH方法有以下優(yōu)勢和特性:

a．由于SPH方法是一種純拉格朗日粒子法，使得其可以方便地描述水流的自由表面，而不需要傳統(tǒng)歐拉方法中對自由面的捕捉和重建操作，這也是SPH非常適合處理具有復雜自由面問題的原因。

b．由于SPH的無網(wǎng)格特性，使其可以靈活地處理流固耦合問題而不需要采用傳統(tǒng)歐拉方法中復雜的網(wǎng)格重建或者網(wǎng)格切割操作，大大簡化了流固耦合求解流程，使其非常適合模擬流固耦合問題。

c．基于拉格朗日法特性，SPH避免了對非線性對流項的離散操作，簡化了求解過程并且避免了相應的數(shù)值誤差。

d．不同于歐拉網(wǎng)格法中的處理，單相SPH模型不需要對計算域內(nèi)非流體區(qū)域劃分網(wǎng)格，減少了對計算資源的額外消耗，對于非連續(xù)流體(如飛濺水流)問題的模擬更加高效。

本文將圍繞SPH的上述特點和優(yōu)勢，著眼于水利與海洋工程中各類實際問題，在對其適用的問題進行分類的基礎上，從工程應用的角度論述SPH在水動力學領域的研究進展，并對SPH未來的發(fā)展方向進行展望。

2 SPH在純水流問題中的應用

SPH被大規(guī)模應用在水利和海洋工程研究領域是從Monaghan開始的，他在文獻［7］中詳細討論了SPH的自由面邊界處理方法、固體邊界處理方法以及N-S方程的求解流程等問題，為以后SPH在水利工程中的應用研究奠定了基礎。SPH在水動力學問題中最廣泛的應用領域是純水流問題，即河流、海洋中的各種水體流動問題。

2.1 波浪模擬

單純的波浪問題是海洋工程中最基本的水動力學問題，也是迄今為止SPH研究的熱點問題之一。波浪屬于典型的自由面大變形問題，傳統(tǒng)歐拉網(wǎng)格法中對其模擬需要依賴高質(zhì)量的界面追蹤算法如VOF(volume of fluid)模型［8］、Level set方法等，這些方法均要求對自由面附近網(wǎng)格進行特殊處理，重建自由面位置并保證自由面位置清晰，而SPH基于無網(wǎng)格離散，流體自由面被拉格朗日粒子自動描述，省去了自由面重建操作，簡化了模擬［9-11］。SPH波浪問題的應用研究主要集中在粒子法中造浪的方法及波浪與海岸結(jié)構(gòu)物相互作用的模擬方面。

2.1.1 造波方法

SPH作為一種粒子法，其造波技術(shù)也有其相應的特點和要求。Monaghan［7］最早進行了相關工作，他根據(jù)孤立波的理論波形設置了初始粒子分布，并給定粒子群孤立波的理論速度場從而簡單地產(chǎn)生波浪。很顯然，該方法很難輸出需要反復輸入動量的周期波。Monaghan在文獻［7］中也提出了推波板造波法，即使用一個類似物理固體邊界的粒子墻，按照實驗造波理論給定其運動規(guī)律，進而產(chǎn)生波浪。該方法應用非常廣泛，且涉及的波浪形式也多種多樣，如Gotoh等［12］在其MPS模型中(一種和 SPH非常類似的模型)使用了造波板產(chǎn)生橢圓余弦波(cnoidal wave)并研究了波浪淺化破碎的特性。在后續(xù)研究中，周期波［13］、孤立波［14］甚至是隨機波［15］都被成功模擬。另外，還曾有學者對運動結(jié)構(gòu)物運動產(chǎn)生的波浪［16］以及滑坡體涌浪［17］進行了一些研究，這些內(nèi)容已經(jīng)涉及了一些流固耦合的問題。

如前面所述，推板造波取得了巨大的成功，但也存在著無法克服的缺陷，即存在二次反射問題，當計算域中存在結(jié)構(gòu)物時，結(jié)構(gòu)物的反射波重新與固體造波板作用形成二次反射而擾亂入射波參數(shù)。一個比較簡單的回避方式是增加計算域延后二次反射到達的時間［18］，但僅能用于短時間波浪模擬情況。Frigaard等［19］提出了一種主動吸波造波板理論，造波板運動的一個額外的分量恰好吸收掉反射波而避免了二次反射，但這種方法非常依賴對反射波的監(jiān)測。Hayashi等［20］在 MPS模型中使用主動吸波造波板實現(xiàn)了長時間波浪破碎和越頂模擬;而Didier等［21］則在SPH中整合了無反射和黏性吸波層，研究了波浪的相互作用問題。

另一種方法被稱為“內(nèi)造波法”(internal wave maker)，它通過源項代替造波板產(chǎn)生波浪。源項只作用于源項區(qū)域(source region)內(nèi)部的流體粒子，通過周期變化源項驅(qū)動源項區(qū)流體發(fā)生運動，產(chǎn)生水面波動，形成向外傳播的行進波。結(jié)構(gòu)物的反射波會穿過源項區(qū)域然后被消波層耗散，避免了二次反射。Larsen等［22］是最早使用該方法進行造波的人，他們基于Boussinesq方程實現(xiàn)了質(zhì)量源項造波。之后的學者把質(zhì)量源項拓展到了緩坡方程［23］和N-S方程［24］，又在其基礎上發(fā)展了動量源項造波［25-27］。而由于SPH方法的拉格朗日特性，使得其難于進行質(zhì)量源項造波法中的質(zhì)量吞吐操作，所以動量源項造波法相比之下更加適合。基于動量源，Liu等［28］改進了內(nèi)造波算法并將其引入到ISPH模型中，系統(tǒng)地測試了SPH內(nèi)造波算法對各種源項區(qū)域尺寸、波浪參數(shù)的敏感性，并將其應用到典型的波浪結(jié)構(gòu)物相互作用問題中，展示了ISPH內(nèi)造波方法的應用潛力。

2.1.2 波浪破碎模擬

波浪破碎、水體抨擊等屬于典型的自由面劇烈變形、水體嚴重不連續(xù)的物理現(xiàn)象，是傳統(tǒng)網(wǎng)格法比較難于處理的方向，也是SPH的優(yōu)勢所在。由于拉格朗日特性，SPH不需要進行FDM、FVM中的追蹤自由面操作，算法簡潔有效。對于破碎波的模擬，技術(shù)問題主要包含兩方面內(nèi)容，一是自由面處的邊界處理，二是有效的紊流模型。

SPH不進行自由面重建操作，但對自由面處的積分域截斷問題仍需要進行額外的處理，最傳統(tǒng)的簡單方法是根據(jù)自由面處粒子數(shù)密度(number density)的變化判斷自由面粒子身份，然后對自由面粒子通過鏡像附近粒子的方式補全積分域，施加邊界條件。該方法作為SPH邊界處理的基本方法，被現(xiàn)今幾乎所有的SPH模型采用或者改進［29-30］。一種比較典型的改進方式是在上述理論基礎上增加粒子對稱性條件，判斷粒子周圍水粒子分布對稱性，輔助判斷自由面身份。該方法由Khayyer等［31］提出，采用相鄰粒子坐標積分來加強判定準確率，類似的改進方法則被Liu等［30］用于流固耦合問題模擬，采用臨近粒子坐標的SPH積分代替直接求和，并發(fā)現(xiàn)該方法對結(jié)構(gòu)物附近自由面粒子誤判問題改善顯著。

基于這些技術(shù)，有大量的SPH破碎波模擬研究成果發(fā)表。使用 WCSPH 模型，Colagrossi等［32-33］分別研究了復雜形狀結(jié)構(gòu)物附近的波浪破碎問題，而Shao［34］則進行了一系列的工作，探索ISPH在波浪破碎模擬上的表現(xiàn)。Khayyer等［35］通過改進上述自由面判別條件和處理，增加了粒子角動量守恒，較大地提高了自由面判定精度。

對破碎波模擬的另一重要方面是紊流模型的建立。SPH被提出至今已經(jīng)近40年，即便是其在流體動力學中的大范圍應用也已經(jīng)有20多年了，但仍然沒有很好地解決紊流問題。學術(shù)界普遍認為的原因是由于傳統(tǒng)SPH方法離散精度比網(wǎng)格法稍低，存在一些無法消除的數(shù)值誤差而掩蓋了紊流的脈動速度和壓強。這使得在更高階精度的SPH格式成熟應用之前的紊流模擬研究都不會達到很高的精度，這在一定程度上限制了SPH紊流模型的發(fā)展。2002年，SPH的鼻祖Monaghan提出了一種求解N-S方程的垂向二維SPH紊流模型，雖然模擬結(jié)果相對較好，但隨后就被指出其計算效率太低。Violeau等［36］提出了兩種改進的代數(shù)紊流模型和一種基于雷諾應力模型的SPH紊流處理方法，在模擬類似潰壩的問題中確實能夠表現(xiàn)出更高的精度，但仍然弱于同等情況下網(wǎng)格法中的精度。同時，Lo等［37］也提出了基于大渦模擬(LES)和k-ε模型的ISPH紊流模型，研究了破碎波復雜自由面附近紊動能分布。Shao等使用這些模型模擬了一系列的波浪、水流問題［38］，其中LES模型在后續(xù)的SPH中應用相對比較廣泛［39］。

2.2 明渠流模擬

上述波浪模型有一個共同的特點，即在模擬區(qū)域并無流體進出，計算域整體的質(zhì)量也是守恒的或者是方便進行拉格朗日描述的。由于SPH的拉格朗日特性，能夠保證進口和出口區(qū)域粒子無進出且自動守恒，而對于需要在進出口設置入流和出流等歐拉描述的邊界條件的情況，SPH則會遇到一些困難。在傳統(tǒng)歐拉網(wǎng)格法中，進出口邊界被固定地設置在一系列(通常是直線或者平面)網(wǎng)格的邊界處，邊界另一側(cè)布置一層虛擬網(wǎng)格，用來給定入流或者出流速度、VOF等信息從而施加入流、出流邊界條件，相對比較方便和直接。而在拉格朗日模型下施加這類歐拉邊界條件則需要一些額外的處理。在水利工程中，常常需要進行一段河道或者部分實驗水槽等實際問題的模擬，對這類入流出流邊界處理的依賴較大，故有一些學者做了這方面的嘗試。

為了解決邊界處粒子進出問題，KAJTAR等［40］采用對進出口流體粒子的直接引入或刪除來簡易地實現(xiàn)進出流，但他們并未給出進出流處理的細節(jié)，而是主要著眼于流固耦合問題。一般認為，直接刪除這些粒子導致的邊界處的截斷誤差(truncation kernel error［41］)會非常大，且對于壓力場的誤算往往影響很大一塊區(qū)域。Lastiwka等［42］稍后提出了一個更加細致的處理方法，但這種方法僅僅適用于氣體動力學一類無自由表面的問題，顯然對于明渠流這類自由面大變形問題并不適用。Federico等［43］于2012年給出了詳細的進出流邊界條件處理細節(jié)，在他們的模型中，進出流邊界外布置了數(shù)層(據(jù)核函數(shù)而定)inflow/outflow粒子，進口邊界粒子被施以指定的速度、壓強并流入計算域，而出口邊界粒子則保持原有速度和壓力場不變。這些邊界粒子只用于避免截斷誤差和施加邊界條件，本身并不通過N-S方程求解，所以流出內(nèi)部粒子支持域的邊界粒子隨即被刪除而不會影響計算域內(nèi)任何粒子。Federico使用改進的模型進行了明渠流的測試，隨后應用于水躍問題，該方法能夠較準確地實現(xiàn)恒定速度進口和自由出流邊界條件。類似地，Bouscasse等［44］模擬了明渠流中圓柱后方流速場，而de Padova等［45］則模擬了三維的水躍問題?？偟膩碇v，對于明渠流問題的模擬研究仍處于起步階段，有很多的問題尚未完全解決。

2.3 淺水流動模擬

在模擬河流甚至是流域尺度問題時，基于N-S方程的數(shù)值模型將面臨運算量過大的問題。對網(wǎng)格法而言，此時可以選擇基于水深平均的方程求解，如淺水方程等。這些方程中流體的垂向參數(shù)被在水深方向積分，也就不存在自由面的捕捉問題，但一般認為SPH方法在淺水方程中的應用優(yōu)勢在于無需處理非線性對流項，簡化了求解流程，且拉格朗日方法良好的守恒性在大尺度模擬時也更有優(yōu)勢。淺水SPH模型具有一些不同于N-S方程模型的獨特性質(zhì)，雖然拉格朗日水粒子體積在模擬過程中不會變化，但由于粒子顯含水深維度的參數(shù)并隨時間變化，故單個水粒子的水平方向參數(shù)(x或x、y)必然發(fā)生變化，導致會牽涉到一系列諸如粒子光滑長度的實時改變甚至是粒子分裂合并的問題，因而淺水SPH模型面臨一些不同于N-S方程模型的困難。但由于大尺度問題對數(shù)值模擬技術(shù)的需求不斷增加，所以相關淺水SPH模型的研究工作也得到了很大的發(fā)展。對于淺水SPH，較早的工作始于Bonet等［46］，他們首次提出了SPH-SWE(shallow water equation)模型，隨后de Leffe等［47］也進行了類似研究。近年來，很多淺水SPH研究集中在對諸如開放邊界［48］、固體結(jié)構(gòu)物邊界、源項［49］等問題上，逐步地實現(xiàn)了對真實河流地形和洪水的模擬［48］。Vacondio等［50］指出淺水SPH模型不同于傳統(tǒng)N-S方程SPH模型的特性之一是在模擬中粒子的光滑長度往往要根據(jù)每個粒子的情況隨位置發(fā)生變化。處理河流問題的淺水SPH與基于N-S方程的傳統(tǒng)SPH具有很大區(qū)別，同時也有更多的問題有待解決。

3 SPH在水流-結(jié)構(gòu)物相互作用問題中的應用研究

3.1 水流-固定結(jié)構(gòu)物相互作用

3.1.1 固定剛性不透水結(jié)構(gòu)物

除了對自由面的模擬外，另外一個SPH擅長的領域是對結(jié)構(gòu)物的模擬。傳統(tǒng)的歐拉網(wǎng)格法中對于拉格朗日描述的結(jié)構(gòu)物的處理常常通過網(wǎng)格重建［51］或切割網(wǎng)格［52］等方法實現(xiàn)結(jié)構(gòu)物、網(wǎng)格之間的信息傳遞;而SPH中對流體和結(jié)構(gòu)物統(tǒng)一采用拉格朗日描述的特性避免了上述問題，拉格朗日流體粒子直接從結(jié)構(gòu)物邊界粒子處獲得信息，施加邊界條件。處理固體結(jié)構(gòu)物的重點在于解決流體粒子在邊界處的截斷誤差和施加邊界條件。在早期模型中采用對邊界附近粒子施加斥力［7］的辦法阻止粒子溢出，但是邊界斥力的值完全是根據(jù)粒子間距計算的，并無物理依據(jù)，這導致了邊界處雖然不發(fā)生粒子溢出，但粒子分布與內(nèi)部差異很大?，F(xiàn)在比較流行的邊界處理方法被稱為虛擬邊界粒子法，包括固定虛粒子和鏡像粒子。前者采用幾層固定于邊界附近的虛粒子阻止流體溢出并施加邊界條件。Randles等［53］是較早改進固定邊界粒子用以施加邊界條件的研究者之一，該技術(shù)隨后被多位學者采用［29，54］。從理論上說，由于固定邊界粒子攜帶了質(zhì)量、壓強、密度等比邊界斥力更多的物理信息，所以比斥力法具有更大的靈活性，也能夠獲得更好的邊界附近流體粒子分布，該方法至今還被廣泛采用。而鏡像粒子法則是從固定邊界粒子法改進而來，鏡像粒子位置根據(jù)邊界另一側(cè)的流體粒子實時更新獲得，同時具有一定的速度，攜帶比固定邊界粒子更多的信息。Morris等［55］是最早使用鏡像粒子概念的研究人員，Cummins等［56］則在其成果中描述了詳細的鏡像規(guī)則，Yildiz等［57］將鏡像粒子邊界處理用于相對復雜的結(jié)構(gòu)物表面的模擬，但在結(jié)構(gòu)物角點處鏡像粒子相互重疊，發(fā)生過度鏡像的問題。Liu等［18］改進了這個方法，防止了角點過度鏡像的發(fā)生，并發(fā)現(xiàn)其在模擬潰壩問題上能夠獲得更加穩(wěn)定光滑的壓力場。

被 SPH 模型研究最多的問題是潰壩問題［13，18，29，58］。實驗中的潰壩水流撞擊一側(cè)擋墻后發(fā)生翻卷水舌并破碎，期間存在大量的不連續(xù)和液滴飛濺現(xiàn)象，這種情況在傳統(tǒng)網(wǎng)格法中模擬起來相對困難，而SPH方法能夠輕易地捕捉到水舌的翻卷和不連續(xù)的液滴飛濺。SPH壓力場計算結(jié)果相對穩(wěn)定光滑，對抨擊點高壓區(qū)的捕捉也比較準確。另一類SPH應用熱點則集中在周期性的波浪和海工建筑物(防波堤、海墻等)相互作用問題，如Shao［59］模擬了波浪與幕墻式防波堤的相互作用，得到了波浪的形變和各種水動力學參數(shù);而Gotoh等［60］則通過模擬波浪和非淹沒式防波堤的相互作用，指出防波堤對波浪的主要耗散作用來自于附近漩渦的生成和脫落;相比之下，三維SPH模型也得到了一定程度的發(fā)展，且通常能夠應用在更加符合物理實際情況的模擬，如Crespo等［61］開發(fā)了3D-SPH模型模擬更大尺度的長波與海岸結(jié)構(gòu)物的作用，而Altomare等［62］則使用并行化的DualSPHysics模型模擬了armour block防波堤的波浪作用特性。在Altomare等的模擬中，由于在塊體間會摻入大量的流體粒子，所以該模擬需要巨大的粒子數(shù)和極大的內(nèi)存消耗，計算中，總計2146095個流體粒子被用于波浪模擬，有效的并行計算框架是保證這樣規(guī)模的模擬能夠順利進行的必要條件。

3.1.2 多孔介質(zhì)結(jié)構(gòu)物

上述波浪與結(jié)構(gòu)物相互作用模擬僅限于不透水剛體結(jié)構(gòu)物，還有部分研究人員著眼于多孔介質(zhì)與水流的相互作用，這些模擬技術(shù)在沙床、多孔防波堤以及植物區(qū)水流模擬等領域應用需求較大。廣義上的多孔介質(zhì)可以是錯綜排列的微小結(jié)構(gòu)物(如植物區(qū))，同網(wǎng)格法中的處理類似，SPH并不直接離散多孔結(jié)構(gòu)物細節(jié)，而是通過求解空間平均后的N-S方程來反應多孔介質(zhì)對水流平均流速的影響。這方面較早的工作始于日本京都大學的 Gotoh等［12］的MPS模型，他們通過在N-S方程中引入一個簡單的摩擦項來描述多孔介質(zhì)影響，模擬了孤立波在透水沙坡上的變形問題。通過改進上述方法，很多相關研究［63-65］對SPH多孔介質(zhì)模型的建立進行了系統(tǒng)的驗證。Shao［66］通過在多孔介質(zhì)區(qū)域內(nèi)外使用不同方程的方法提出了ISPH多孔介質(zhì)模型，并模擬了周期波在淺灘破碎等問題。

3.2 水流-剛體結(jié)構(gòu)物耦合運動

對于結(jié)構(gòu)物運動和水流耦合計算的處理，分為運動學耦合和動力學耦合兩大類，前者指運動結(jié)構(gòu)物速度等參數(shù)事先給定，典型的如液箱晃蕩問題等［67-69］;相比之下，動力學耦合模型除了需要精確的邊界條件處理方法和運動邊界的處理方法外，還需要計算運動結(jié)構(gòu)物受力的方法。在傳統(tǒng)的網(wǎng)格法中，對于拉格朗日結(jié)構(gòu)物的描述需要一整套額外的自適應網(wǎng)格或者切割網(wǎng)格理論而顯得復雜，尤其是當運動結(jié)構(gòu)物處于自由面附近甚至穿透自由面時，處理起來會遇到更多的困難，實際工程中港口船舶的錨固、各種浮式結(jié)構(gòu)物的波浪響應、水庫滑坡體入水涌浪等問題卻都涉及到運動結(jié)構(gòu)物穿破自由面的情況，對數(shù)值模型提出了更高的要求。在現(xiàn)有的大部分SPH模型中，邊界條件都是由虛擬邊界粒子(固定粒子或鏡像粒子)施加的，加上流體直接用拉格朗日粒子描述，所以它更加適合對上述流固耦合問題的模擬，通過引入簡單的運動邊界處理和加速度求解模塊，即可實現(xiàn)流固耦合。流固耦合計算的主要難點是準確地求解物體受力及運動軌跡，對此，已有的SPH算法一般可分為動量傳遞法和積分表面壓力法兩大類。

3.2.1 動量傳遞法

早期的SPH模型(尤其是MPS模型)由于各種原因，對流體壓強計算時會產(chǎn)生較大的數(shù)值壓力震蕩。在進行流固耦合模擬時需要積分結(jié)構(gòu)物表面的壓強從而獲得壓力的合力，這個合力也會發(fā)生震蕩而影響后續(xù)加速度計算。MPS方法的創(chuàng)始人Koshizuka提出過一個巧妙的辦法回避壓強積分這個操作［70］，這個方法隨后被 Shao 等［71-72］引入 SPH方法用于運動結(jié)構(gòu)物的受力及運動模擬。在他們的模型中，結(jié)構(gòu)物被離散成和固體密度一樣的固體粒子群，其初始間距和流體粒子相同，在計算中，固體粒子與流體粒子一樣參與求解N-S方程并獲得速度，在N-S方程求解步驟之后，按實際固體密度積分所有固體粒子的動量，計入結(jié)構(gòu)物運動的總動量，再根據(jù)結(jié)構(gòu)物剛體假設重新分配這部分動量。在這樣處理之后，計算中的總動量并沒有丟失而保持守恒，結(jié)構(gòu)物的剛體假設也能夠保持。在下一時刻初始，經(jīng)過剛體假設更新的固體粒子的位置會對周圍流體施加作用。這個方法的應用并不廣泛，能夠查詢到的文獻只有B?ckmann等發(fā)表的，他們進行了一系列流固耦合問題的模擬研究［73］。這個方法整個流程中并未涉及對壓強的積分，也就避免了因為壓強震蕩而導致的誤差，它通過動量守恒特性，先使用固體粒子參與計算“收集”由流體傳來的動量，隨后根據(jù)剛體假設和獲得的總動量等物理特性重新分配而更新結(jié)構(gòu)物參數(shù)。但這個方法的一個缺點是結(jié)構(gòu)物附近流體“不可壓”特性會被破壞，且由于粒子質(zhì)量的差異，結(jié)構(gòu)物附近粒子分布也會更不規(guī)則［30］。

3.2.2 積分力法

相比之下，另一種方法應用更加廣泛，通過積分結(jié)構(gòu)物附近(表面)粒子的壓強(或結(jié)構(gòu)物邊界粒子壓強)求得壓力，結(jié)合結(jié)構(gòu)物受力情況依據(jù)牛頓第二定律計算結(jié)構(gòu)物加速度和角加速度，時間積分平動速度和轉(zhuǎn)動速度及位置。理論上講，牛頓第二定律本質(zhì)也是動量守恒，故和第一種方法相比具有相同的理論基礎，不過對流體力直接積分的方法更加符合流固耦合的物理實際。考慮到計算流程的簡潔，大部分SPH流固耦合模型建立在WC-SPH下，如Vandamme等［74］采用固定邊界粒子(見 3.1.1節(jié))方法的WCSPH模型模擬了楔形體入水和圓柱入水問題;而Oger等［75］則基于鏡像粒子繼續(xù)了流固耦合的研究，他們引入了鏡像邊界附近壓強積分的算法，詳細介紹了流固耦合的處理并給出了一系列運動結(jié)構(gòu)物的模擬結(jié)果，展示了模型的應用。Liu等［30］綜合前人成果改進了一套用于處理流固耦合的ISPH模型，模擬了液箱晃蕩和方柱入水等問題。液箱晃蕩問題中，他們使用了真實的箱體運動外激勵速度晃蕩箱體，代替在網(wǎng)格法中普遍使用的非慣性坐標系［76］，模擬結(jié)果和實驗對比良好，隨后他們將其應用于更加復雜的方柱入水問題中，該問題結(jié)構(gòu)物運動是包含平動和轉(zhuǎn)動的全自由運動，ISPH結(jié)果較好地還原了方柱各個自由度的運動和自由面的變化，模擬結(jié)果和FLUENT軟件計算結(jié)果相差不大。在水利工程實際問題方面，筆者也將該模型拓展到對滑坡體涌浪的產(chǎn)生、傳播和翻壩問題的模擬上，并通過實驗驗證了SPH結(jié)果的精度，展示了ISPH在水利工程問題中的應用潛力。

3.3 水流-柔性結(jié)構(gòu)物模擬

上述運動結(jié)構(gòu)物的處理算法都僅限于剛體結(jié)構(gòu)物，對于可變形的彈性/柔性結(jié)構(gòu)物來說，情況則更加復雜一些。如前面所述，基于網(wǎng)格的數(shù)值方法在處理柔性物體時會面臨比剛體更嚴重的困難。實際計算中對彈性結(jié)構(gòu)物的模擬通常使用FEM方法或者類似的拉格朗日法，這使得原本就基于拉格朗日描述的SPH具有更大的理論優(yōu)勢。最早使用SPH模擬上述問題的是Müller等［77］，在他們的方法中，基于Lennard-Jones勢引入邊界力來模擬彈性結(jié)構(gòu)物，但他們的模擬結(jié)果在彈性邊界處會發(fā)生嚴重的粒子分布不均勻現(xiàn)象，精度較差。Lenaerts等［78］則模擬了柔性外殼與水的相互作用，他們采用了邊界力和位置修正聯(lián)合作用的辦法防止固體與流體粒子的互相穿透，但也指出粒子穿透現(xiàn)象在材料彈性模量較小時很難避免。Amini等［79］則系統(tǒng)研究了多種柔性/彈性結(jié)構(gòu)物的流固耦合問題，這些類似問題在傳統(tǒng)網(wǎng)格法中是很難模擬的。

4 SPH在泥沙輸運問題中的應用

水流中的泥沙輸運問題是水利和海洋工程中極其重要的一大類問題。早期的關于泥沙的數(shù)值模擬研究都是基于網(wǎng)格法的，如基于歐拉網(wǎng)格的有限差分法(FDM)［80］和有限體積法(FVM)［81］，以及基于拉格朗日網(wǎng)格的有限元法(FEM)［82］。另外還有一種特殊的網(wǎng)格法，即格子玻爾茲曼方法(LBM)［83］。由于SPH的粒子特性，能夠基于泥沙顆?；蛘吡Ｗ拥状驳雀拍顏砻枋黾澳M泥沙問題，具有一定的優(yōu)勢，目前SPH方法對泥沙的處理模型大致分為泥沙顆粒兩相流模型和底床高度沖刷模型。

4.1 SPH泥沙顆粒兩相流模型

利用拉格朗日粒子概念描述非連續(xù)的泥沙顆粒具有一定的優(yōu)勢，它可以實現(xiàn)泥沙-水流的“雙向耦合”，即在泥沙被水流沖蝕的同時，水流也會受到泥沙的反作用，且這一切是由SPH水粒子和泥沙粒子混合物自動完成的，不需太多的額外處理。Monaghan等［84］最早將這樣的處理加入到SPH方法中，而Shakibaeinia等［85］使用了類似SPH的方法改進MPS方法研究了潰壩水流導致的泥沙輸運和底床變形。在這些研究中，底床被可沖蝕的泥沙粒子代替，當?shù)状布袅_到設定的臨界值時，泥沙粒子開始啟動，啟動后的泥沙顆粒作為大密度的流體粒子參與流體計算，泥沙顆粒的沉降也完全按照其真實受力特性計算。

4.2 SPH底床高度沖刷模型

基于泥沙和水流的兩相流模型具有處理簡單、邊界條件易于施加等特點，且實踐證明其在具有底床大變形和水沙混摻劇烈的問題中(如3.1.1節(jié)中的潰壩)模擬精度良好、可靠性高。但單個泥沙顆粒的啟動模型存在一定的缺陷:首先，在實際模擬中SPH泥沙粒子體積遠遠大于真實的沙粒，故必須等待底床沖蝕積累到一定量之后才能啟動，在模擬中造成了顯著的非連續(xù)底床變動，對于精細且緩慢的動床變化計算誤差較大;其次，游離在水中的單個泥沙粒子的水動力學特性和它所表征的一團細顆粒泥沙差別較大，而將其作為一個整體通過自重沉降回底床的做法可能帶來很大誤差。從現(xiàn)有研究結(jié)果來看，大多數(shù)SPH泥沙輸運模擬均只關注了比較劇烈的諸如潰壩、抨擊等問題的泥沙輸運，可能也與上述缺陷有一定關系。因此，Kri?tof等［86］嘗試在 SPH 中引入類似網(wǎng)格法的連續(xù)底床沖刷模型來實現(xiàn)連續(xù)的底床變形和推移質(zhì)輸運。在他們的方法中，底床仍然用泥沙粒子來表征形變和施加邊界條件，但該泥沙顆粒并不啟動，而是根據(jù)水流條件實時計算沖蝕并由推移質(zhì)輸沙率和懸起沉降率共同保證質(zhì)量守恒。同時，該方法還依據(jù)類似網(wǎng)格法中的處理加入了基于拉格朗日描述的懸移質(zhì)對流擴散方程，將懸移質(zhì)離散到水粒子系統(tǒng)中并通過拉格朗日粒子懸移質(zhì)濃度概念求解對流擴散方程，實現(xiàn)真實的泥沙輸運和底床變形模擬。

5 SPH氣液兩相流模擬應用研究

SPH方法在流體模擬中的另一個應用領域是對于多相流、特別是氣液兩相流的模擬［87］。由于傳統(tǒng)單相SPH認為空氣壓強為常數(shù)，造成了當水流一側(cè)壓強過大時會將自由面向空氣一側(cè)推動，這在描述水中氣泡時顯然是錯誤的。氣液兩相SPH模型的最初并未局限于氣液兩相而是指廣義上的兩種不同密度流體的耦合模擬，Monaghan等［88］最早將兩相流概念引入SPH，他們使用了兩種不同粒子描述不同流相，相間的相互作用根據(jù)N-S方程計算而兩相界面處則增加一個非物理的作用力實現(xiàn)耦合，用來模擬dust gas流。后續(xù)研究將其應用在很多領域［32，36，84］。這些早期的 SPH 方法均基于微可壓縮模型，在處理密度差很大的兩相流交界面時會遇到一些困難。如Colagrossi等［32］曾指出WCSPH在模擬大密度比的兩相流時，需要在交界面處做非常大的限制(一個密度額外修正處理、一個施加在低密度流體上的增大表面張力的處理、一個速度光滑處理等)才能保證精度，并且需要在兩種流體間設置差異極大的聲速來均衡WCSPH中聲速帶來的影響。這些做法除了嚴重增加算法的非物理限制外，還大大增加了運算量，WCSPH方法中過多的人工耗散和修正的弊端在兩相流問題面前暴露無遺。

不可壓縮模型因為不需要采用人工聲速，可以較好地處理水氣兩相流，其中比較有代表性的是Hu等［89-90］的一系列研究成果，其ISPH模型被成功地用于模擬大密度差的兩相流問題。Shao［91］則提出了類似的ISPH模型，將其同時用于大密度差和小密度差兩相流的模擬，展示了ISPH模型良好的通用性和穩(wěn)定性。在Shao的模型中，依然采用兩種粒子來分別模擬兩相(如水和氣)，各相分別獨自求解一套N-S方程描述流體粒子的運動，而交界面上則使用另一相粒子位置充當邊界粒子、施加邊界條件，這樣的處理類似多孔介質(zhì)［66］甚至是運動結(jié)構(gòu)物中的處理。這種ISPH模型在計算過程中不發(fā)生密度變化，保證了交界面兩側(cè)各自流體密度的恒定?，F(xiàn)有的大部分ISPH模型對于氣體相都處理成不可壓流體，這在某些情況下并不符合物理實際，需要改進模型。

兩相流問題中，雖然可以忽略水的壓縮性，但氣體的壓縮性是不能忽略的，如水中大氣泡的運動、破碎波中卷積的氣囊等。這就需要采用不同的模型模擬水體和氣體，目前這方面的工作還很少，只在近兩年，新加坡國立大學Luo等通過粒子法實現(xiàn)了對可壓縮氣體與不可壓縮液體的模擬。另一個水動力學中常見的問題是當水面紊動很強時，會在水面發(fā)生自摻氣，可摻入的氣泡又會在紊流變?nèi)鯐r逸出水面。這個問題的模擬尚未見報道。Luo等①LUO M，KOH C G，GAO M，et al.A particle method for two-phase flows with large density difference［J］.International Journal for Numerical Methods in Engineering，2015(已錄用，待刊).基于兩相CPM(consistent particle method，基于 MPS改進)方法改進了模型，使之能夠處理晃蕩問題中的氣腔形變，模擬結(jié)果和實驗吻合良好。

6 研究展望

通過分析發(fā)現(xiàn)，SPH具有穩(wěn)定簡潔的特點，在交界面大變形問題的模擬上具有優(yōu)勢，而在未來則會有更多的關于泥沙輸運和水流摻氣的SPH研究成果出現(xiàn)。而復雜的多物理過程的研究也將成為熱點，如需要考慮材料特性的柔性體變形模擬，需要考慮破壞土體特性的滑坡體入水問題，以及需要考慮固液氣三相的泥石流運動問題等。

在對純水流問題與流固耦合問題的模擬中，SPH已經(jīng)取得了很大的進展，計算方法相對成熟穩(wěn)定，已經(jīng)可以與傳統(tǒng)的網(wǎng)格計算法分庭抗禮，而在交界面大變形(自由面和結(jié)構(gòu)物表面)問題中，SPH還略占優(yōu)勢。但在泥沙模擬和氣液兩相流模擬方面，尤其是動態(tài)摻氣過程的模擬中，SPH方法才剛剛起步，尚有大量工作需要逐步推進?？紤]到SPH方法在處理多相流方面具有自身的優(yōu)勢，可以預見在未來會有越來越多的關于泥沙輸運和水流摻氣的SPH研究成果出現(xiàn)。同時采用SPH方法模擬復雜的多物理過程問題的研究也將出現(xiàn)，這些問題包含需要考慮材料特性的柔性體變形模擬，需要考慮破壞土體特性的滑坡體入水問題，以及需要考慮固液氣三相的泥石流運動問題等。

未來SPH在水利與海洋工程研究領域的發(fā)展將主要集中在以下幾個方面:①SPH計算效率與精度的提高。目前大部分成熟的水動力學SPH算法均基于使用各向同性均勻分布粒子的核函數(shù)，在計算多尺度物理問題時，面臨計算效率低的問題。在結(jié)合N-S方程求解方法的基礎上，開發(fā)出高階精度的具有橢圓支持域的非均勻各向異性核函數(shù)的SPH模型，可以有效提高模型性能，用于模擬狹長計算域問題并提高局部模擬精度(如長波與結(jié)構(gòu)物作用問題)，是 SPH未來的重要的理論創(chuàng)新方向之一。②在高精度SPH模型開發(fā)基礎上，發(fā)展適用于SPH特點的紊流模型，如大渦模擬，更好地模擬各類實際問題，可以大大拓寬SPH的工程應用地位，擺脫目前SPH多用于算法研究的困境。③為了克服SPH相比網(wǎng)格法更加耗時的弱點，對三維SPH模型實現(xiàn)并行計算是必要的發(fā)展趨勢。在這些研究成功的基礎上，SPH方法將獲得更為廣闊的發(fā)展空間，并展現(xiàn)出其獨特的優(yōu)勢與迷人的前景。

［1］LIU M B，LIU G R.Smoothed particle hydrodynamics(SPH):an overview and recent developments［J］.Archives of Computational Methods in Engineering，2010，17:25-76.

［2］GINGOLD R A，MONAGHAN J J.Smoothed particle hydrodynamics:theory and application to non-spherical stars［J］.Montly Notices of the Royal Astronomical Society，1977，181(3):375-389.

［3］LUCY L B.A numerical approach to the testing of fission hypothesis［J］.The Astronomical Journal，1977，82(12):1013-1024.

［4］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:8417-8436.

［5］MORRIS J P，ZHU Y，F(xiàn)OX P J.Parallel simulation of pore-scale flow through porous media［J］.Computers and Geotechnics，1999，25:227-246.

［6］LIU G R，LIU M B.Smoothed particle hydrodynamics:a meshfree particle method［M］.Singapore:World Scientific，2003.

［7］MONAGHAN J J.Simulating free surface flows with SPH［J］.Journal of Computational Physics，1994，110(2):399-406.

［8］LIN P Z，XU W L.NEWFLUME:a numerical water flume for two-dimensional turbulent free surface flows［J］.Journal of Hydraulic Research，2006，44(1):79-93.

［9］高睿，任冰，王國玉，等.孤立波淺化過程的SPH數(shù)值模擬［J］.水動力學研究與進展:A 輯，2010，25(5):620-629.(GAO Rui，REN Bing，WANG Guoyu，et al.SPH model of solitary waves shoaling on a mild sloping beach［J］.Chinese Journal of Hydrodynamics，2010，25(5):620-629.(in Chinese))

［10］劉謀斌，宗智，常建忠.光滑粒子動力學方法的發(fā)展與應用［J］.力學進展，2011，41(2):217-234.(LIU Moubin，ZONG Zhi，CHANG Jianzhong.Development and application of smoothed particle hydrodynamics［J］.Advances in Mechanics，2011，41(2):217-234.(in Chinese))

［11］鄭興，段文洋.K2_SPH方法及其對二維非線性水波的模擬［J］.計算物理，2011，28(5):659-666.(ZHENG Xing，DUAN Wenyang.K2_SPH method and application for 2D nonlinear water wave simulation［J］.Chinese Journal of Computational Physics，2011，28(5):659-666.(in Chinese))

［12］GOTOH H，SAKAI T，Lagrangian simulation of breaking waves using particle method［J］.Coastal Engineering Journal，1999，41(3/4):303-326.

［13］KHAYYER A，GOTOH H.Modified moving particle semiimplicit methods for the prediction of 2D wave impact pressure［J］.Coastal Engineering，2009，56(4):419-440.

［14］DAO MH，XU H，CHAN ES，TKALICH P.Modelling of tsunami-like wave run-up，breaking and impact on a vertical wall by SPH method［J］.Natural Hazards and Earth System Sciences，2013，13:3457-3467.

［15］ZHENG J H，SOE M M，ZHANG C，et al.Numerical wave flume with improved smoothed particle hydrodynamics［J］.Journal of Hydrodynamics:Ser B，2010，22(6):773-781.

［16］GALLATI M，BRASCHI G，F(xiàn)ALAPPI S.SPH simulations of the waves produced by a falling mass into a reservoir［J］.Nuovo Cimento Societa Intalian di Fisica C:Geophysics and Space Physics，2005，28:129.

［17］DU X，W U W，GONG K，et al.Two dimensional SPH simulation of waterwaves generated by underwater landslide［J］.Journal of Hydrodynamics:Ser A，2006，21(5):579-586.

［18］LIU X，XU H H，SHAO S D，et al.An improved incompressible SPH model for simulation of wave-structure interaction［J］.Computers and Fluids，2013，71:113-123.

［19］FRIGAARD P，BRORSEN M.A time domain method for separating incident and reflected irregular waves［R］.Aalborg，Danish:Hydraulics ＆ CoastalEngineering Laboratory，Department of CivilEngineering，Aalborg University，1993.

［20］HAYASHI M，GOTOH H，MEMITA T，et al.Gridless numerical analysis of wave breaking and overtopping at upright seawall［C］//Coastal Engineering Conference.San Francisco:Asce American Society of Civil Engineers，2001:2100-2113.

［21］DIDIER E，NEVES M G.A semi-infinite numerical wave flume using smoothed particle hydrodynamics［J］.International Journal of Offshore and Polar Engineering，2012，22(3):193-199.

［22］LARSEN J，DANCY H.Open boundaries in short wave simulations-a new approach［J］.Coastal Engineering，1983，7(3):285-297.

［23］LEE C，SUH K D.Internal generation of waves for timedependent mild slope equations［J］.Coastal Engineering，1998，34(1):35-57.

［24］LIN P，LIU P L F.Internal wave-maker for Navier-Stokes equations models［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，1999，125(4):207-215.

［25］WEI G，KIRBY J T，SINHA A.Generation of waves in Boussinesq models using a source function method［J］.Coastal Engineering，1999，36(4):271-299.

［26］CHOI J，YOON S B. Numerical simulations using momentum source wave-maker applied to RANS equation model［J］. Coastal Engineering，2009，56 ( 10 ) : 1043- 1060.

［27］HA T，LIN PZ，CHO Y S.Generation of 3D regular and irregular waves using Navier-Stokes equations model with an internal wave maker［J］.Coastal Engineering，2013，76:55-67.

［28］LIU X，LIN PZ，SHAOSD.ISPH wave simulation by using an internal wave maker［J］.Coastal Engineering，2015，95:160-170.

［29］SHAO S D，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.

［30］LIU X，LIN P，SHAO S D.An ISPH simulation of coupled structure interaction with free surface flows［J］.Journal of Fluids and Structures，2014，48:46-61.

［31］KHAYYERA，GOTOH H，SHAOSD.Enhanced predictions of wave impact pressure by improved incompressible SPH methods［J］.Applied Ocean Research，2009，31:111-131.

［32］COLAGROSSI A，LANDRINI M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics［J］.Journal of Computational Physics，2003，191(2)，448-475.

［33］MONAGHAN J J，KOS A，ISSA N.Fluid motion generated by impact［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2004，129:250-259.

［34］SHAO S D.Incompressible SPH simulation of wave breaking and overtopping with turbulence modelling［J］.International Journal for Numerical Methods in Fluids，2006，50(5):597-621.

［35］KHAYYERA，GOTOH H，SHAO SD.Corrected Incompressible SPH method for accurate water-surface tracking in breaking waves［J］.Coastal Engineering，2008，55(3):236-250.

［36］VIOLEAU D，BUVAT C，ABED-MERA?M K，et al.Numerical modelling of boom and oil spill with SPH［J］.Coastal Engineering，2007，54(12):895-913.

［37］LO E Y M，SHAO S D.Simulation of near-shore solitary wave mechanics by an incompressible SPH method［J］.Applied Ocean Research，2002，24(5):275-286.

［38］SHAO S D，JI C，GRAHAM D I，et al.Simulation of wave overtopping by an incompressible SPH model［J］.Coastal Engineering，2006，53(9):723-735.

［39］DALRYMPLE R A，ROGERS B D.Numerical modeling of waterwaves with the SPH method［J］.Coastal Engineering，2006，53(2):141-147.

［40］KAJTAR J，MONAGHAN JJ.SPH simulationsof swimming linked bodies［J］.Journal of Computational Physics，2007，227:8568-8587.

［41］ANDREA A，JEAN-CHRISTOPHE M，F(xiàn)RANCIS L，et al.SPH truncation error in estimating a 3D function［J］.Computers and Fluids，2011，44:279-296.

［42］LASTIWKA M，BASA M，QUINLAN N J.Permeable and non-reflecting boundary conditions in SPH ［J］.International Journal for Numerical Methods in Fluids，2009，61:709-724.

［43］FEDERICO I，MARRONE S，COLAGROSSI A，et al.Simulating 2D open-channel flows through an SPH model［J］.European Journal of Mechanics-B/Fluids，2012，34:35-46.

［44］BOUSCASSE B，COLAGROSSI A，MARRONE S，et al. Viscous flow past a circular cylinder below a free surface ［C］/ /ASME 33rd International Conference on Ocean，Offshore and Arctic Engineering. San Francisco: merican Society of Mechanical Engineers ( AMSE ) ，2014: V002T08A080-V002T08A080.

［45］de PADOVA D，MOSSA M，SIBILLA S，et al.3D SPH modelling of hydraulic jump in a very large channel［J］.Journal of Hydraulic Research，2013，51(2):158-173.

［46］BONET J，KULASEGARAM S，RODRIGUEZ-PAZ M X，Variational formulation for the smooth particle hydrodynamics(SPH)simulation of fluid and solid problems［J］.Computer Methods in Applied Mechanics and Engineering，2004，193(12):1245-1256.

［47］de LEFFE M，LE TOUZé D，ALESSANDRINI B.SPH modeling of shallow-water coastal flows［J］.Journal of Hydraulic Research，2010，48(Sup1):118-125.

［48］VACONDIO R，ROGERS B D，STANSBY P K，et al.SPH modeling of shallow flow with open boundaries for practical flood simulation［J］.Journal of Hydraulic Engineering，2011，138(6):530-541.

［49］VACONDIO R，ROGERS B D，STANSBY P K，et al.Shallow water SPH for flooding with dynamic particle coalescing and splitting［J］.Advances in Water Resources，2013，58:10-23.

［50］VACONDIO R，ROGERS B D，STANSBY P K.Accurate particle splitting for smoothed particle hydrodynamics in shallow water with shock capturing［J］.International Journal for Numerical Methods in Fluids，2012，69(8):1377-1410.

［51］HASSAN O，PROBERT E J，MORGAN K，et al.Unsteady flow simulation using unstructured meshes［J］.Computer Methods in Applied Mechanics and Engineering，2000，189:1247-1275.

［52］LAI M C，PESKIN C S.An immersed boundary method with formal second-order accuracy and reduced numerical viscosity［J］.Journal of Computational Physics，2000，160:705-719.

［53］RANDLES P W，LIBERSKY L D.Smoothed particle hydrodynamics:some recent improvements and applications［J］.Computer Methods in Applied Mechanics and Engineering，1996，139:375-408.

［54］LIU G R，GU Y T.A local radial point interpolation method(LRPIM)for free vibration analyses of 2D solids［J］.Journal of Sound and Vibration，2001，246(1):29-46.

［55］MORRIS J P，F(xiàn)OX P J，ZHU Y.Modeling low Reynolds number incompressible flows using SPH［J］.Journal of Computational Physics，1997，136:214-226.

［56］CUMMINS SJ，RUDMAN M，An SPH projection method.Journal of Computational Physics［J］.1999，152:584-607.

［57］YILDIZ M，ROOK R A，SULEMAN A.SPH with the multiple boundary tangentmethod［J］.International Journal for Numerical Methods in Engineering，2009，77:1416-1438.

［58］王亞萍，光滑粒子流體動力學(SPH)方法在潰壩問題中的應用研究［D］.天津:天津大學，2012.

［59］SHAO S D.SPH simulation of solitary wave interaction with a curtain-type breakwater［J］.Journal of Hydraulic Research，2005，43(4):366-375.

［60］GOTOH H，SHAO S D，MEMITA T.SPH-LES model for numerical investigation of wave interaction with partially immersed breakwater［J］.Coastal Engineering Journal，2004，46(1):39-63.

［61］CRESPO A J C，GóMEZ-GESTEIRA M，DALRYMPLE R A.3D SPH simulation of large waves mitigation with a dike［J］.Journal of Hydraulic Research，2007，45(5):631-642.

［62］ALTOMARE C，CRESPO A J C，ROGERS B D，et al.Numerical modelling of armour block sea breakwater with smoothed particle hydrodynamics［J］.Computers ＆Structures，2014，130:34-45.

［63］TARTAKOVSKY A M，MEAKIN P.Pore scale modeling of immiscible and miscible fluid flowsusing smoothed particle hydrodynamics［J］.Advances in Water Resources，2006，29(10):1464-1478.

［64］TARTAKOVSKY A M，MEAKIN P，SCHEIBE T D，et al.Simulations of reactive transport and precipitation with smoothed particle hydrodynamics［J］.Journal of Computational Physics，2007，222(2):654-672.

［65］TARTAKOVSKY A M. Langevin model for reactive transport in porous media［J］. Physical Review E，2010， 82( 2) : 026302.

［66］SHAO S D.Incompressible SPH flow model for wave interactions with porous media［J］.Coastal Engineering，2010，57(3):304-316.

［67］陳臻.SPH算法改進及在晃蕩與入水中的應用［D］.大連:大連理工大學，2014.

［68］李大鳴，陳海舟，張建偉，等.基于SPH法的二維矩形液艙晃蕩研究［J］.計算力學學報，2010，27(2):369-374.(LI Daming，CHEN Haizhou，ZHANG Jianwei，et al.A study of 2D liquid sloshing in rectangle tanks based on SPH method［J］.ChineseJournal of Computational Mechanics，2010，27(2):369-374.(in Chinese))

［69］崔巖，吳衛(wèi)，龔凱，等.二維矩形水槽晃蕩過程的SPH方法模擬［J］.水動力學研究與進展:A 輯，2008，23(6):618-624.(CUI Yan，WU Wei，GONG Kai，et al.Numerical simulation of sloshing in two dimensional rectangular tanks with SPH ［J］.Chinese Journal of Hydrodynamics，2008，23(6):618-624.(in Chinese))

［70］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.

［71］SHAO S D，GOTOH H，Simulating coupled motion of progressive wave and floating curtain wall by SPH-LES model［J］.Coastal Engineering Journal，2004，46(2):171-202.

［72］SHAO S D.Incompressible SPH simulation of water entry of a free-falling object［J］.International Journal for numerical methods in fluids，2009，59(1):91-115.

［73］B?CKMANN A，SHIPILOVA O，SKEIE G.Incompressible SPH for free surface flows［J］.Computers ＆ Fluids，2012，67:138-151.

［74］VANDAMME J，ZOU Q P，REEVE D E.Modeling floating object entry and exit using smoothed particle hydrodynamics［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2011，137(5):213-224.

［75］OGER G，DORING M，ALESSANDRINI B，et al.Twodimensional SPH simulations of wedge water entries［J］.Journal of Computational Physics，2006，213:803-822.

［76］XUE M A，LIN P Z.Numerical study of ring baffle effects on reducing violent liquid sloshing［J］.Computers ＆Fluids，2010，52:116-129.

［77］MüLLER M，SCHIRM S，TESCHNER M，et al.Interaction of fluids with deformable solids［J］.Computer Animation and Virtual Worlds，2004，15(3/4):159-171.

［78］LENAERTS T，DUTRé P.Unified SPH model for fluidshell simulations［R］.Leuven，Belgium:Department Computer Science，Katholieke Universiteit Leuven，2008.

［79］AMINI Y，EMDAD H，F(xiàn)ARID M，A new model to solve fluid-hypo-elastic solid interaction using the smoothed particle hydrodynamics(SPH)method［J］.European Journal of Mechanics-B/Fluids，2011，30(2):184-194.

［80］ROACHE P J. Computational fluid dynamics ［M］.Albuquerque，New Mexico: Hermosa Publishers， 1998.

［81］EYMARD R，GALLOUET T，HERBIN R.Finite volume methods［M］.Wroclaw，Poland:University of Wroclaw，2008.

［82］GLOWINSKI R.Finite element methods for incompressible viscous flow［M］.Amsterdam:Elsevier Science B.V.，2003.

［83］THORNE D T，MICHAEL C.Lattice Boltzmann modeling:an introduction for geoscientists and engineers［M］.Berlin:Springer，2006.

［84］MONAGHAN J J，CAS RAF，KOS A M，et al.Gravity currents descending a ramp in a stratified tank［J］.Journal of Fluid Mechanics，1999，379:39-69.

［85］SHAKIBAEINIA A，JIN Y C.A mesh-free particle model for simulation of mobile-bed dam break［J］.Advances in Water Resources，2011，34(6):794-807.

［86］KRI?TOF P，BENE? B，KRˇIVáNEK J，et al.Hydraulic erosion using smoothed particle hydrodynamics［J］.Computer Graphics Forum，2009，28(2):219-228.

［87］龔凱，劉樺，王本龍.兩相 SPH 方法［C］//吳有生，劉樺，許維臨，等.第九屆全國水動力學學術(shù)會議暨第二十二屆全國水動力學研討會文集.北京:海洋出版社，2009:340-345.

［88］MONAGHAN J J，KOCHARYAN A.SPH simulation of multi-phase flow［J］.Computer Physics Communications，1995，87(1):225-235.

［89］HU X Y，ADAMS N A.An incompressible multi-phase SPH method［J］.Journal of Computational Physics，2007，227(1):264-278.

［90］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.

［91］SHAO S D. Incompressible smoothed particle hydrodynamics simulation of multifluid flows［J］.International Journal for Numerical Methods in Fluids，2012，69(11):1715-1735.