林鵬智，劉 鑫

(1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川成都 610065;2.河海大學(xué)港口海岸與近海工程學(xué)院，江蘇南京 210098)

1 光滑粒子水動(dòng)力學(xué)簡(jiǎn)介

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

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

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

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

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

然后再通過(guò)求解壓力Poisson方程(pressure Poisson equation，PPE)獲得壓強(qiáng)，進(jìn)而獲得速度場(chǎng)。

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

1.2 計(jì)算域的SPH離散

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

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

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

圖1 SPH方法中的粒子離散

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

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

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

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

c．基于拉格朗日法特性，SPH避免了對(duì)非線性對(duì)流項(xiàng)的離散操作，簡(jiǎn)化了求解過(guò)程并且避免了相應(yīng)的數(shù)值誤差。

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

本文將圍繞SPH的上述特點(diǎn)和優(yōu)勢(shì)，著眼于水利與海洋工程中各類實(shí)際問(wèn)題，在對(duì)其適用的問(wèn)題進(jìn)行分類的基礎(chǔ)上，從工程應(yīng)用的角度論述SPH在水動(dòng)力學(xué)領(lǐng)域的研究進(jìn)展，并對(duì)SPH未來(lái)的發(fā)展方向進(jìn)行展望。

2 SPH在純水流問(wèn)題中的應(yīng)用

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

2.1 波浪模擬

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

2.1.1 造波方法

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

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

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

2.1.2 波浪破碎模擬

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

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

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

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

2.2 明渠流模擬

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

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

2.3 淺水流動(dòng)模擬

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

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

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

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

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

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

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

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

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

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

3.2.1 動(dòng)量傳遞法

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

3.2.2 積分力法

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

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

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

4 SPH在泥沙輸運(yùn)問(wèn)題中的應(yīng)用

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

4.1 SPH泥沙顆粒兩相流模型

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

4.2 SPH底床高度沖刷模型

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

5 SPH氣液兩相流模擬應(yīng)用研究

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

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

兩相流問(wèn)題中，雖然可以忽略水的壓縮性，但氣體的壓縮性是不能忽略的，如水中大氣泡的運(yùn)動(dòng)、破碎波中卷積的氣囊等。這就需要采用不同的模型模擬水體和氣體，目前這方面的工作還很少，只在近兩年，新加坡國(guó)立大學(xué)Luo等通過(guò)粒子法實(shí)現(xiàn)了對(duì)可壓縮氣體與不可壓縮液體的模擬。另一個(gè)水動(dòng)力學(xué)中常見(jiàn)的問(wèn)題是當(dāng)水面紊動(dòng)很強(qiáng)時(shí)，會(huì)在水面發(fā)生自摻氣，可摻入的氣泡又會(huì)在紊流變?nèi)鯐r(shí)逸出水面。這個(gè)問(wèn)題的模擬尚未見(jiàn)報(bào)道。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改進(jìn))方法改進(jìn)了模型，使之能夠處理晃蕩問(wèn)題中的氣腔形變，模擬結(jié)果和實(shí)驗(yàn)吻合良好。

6 研究展望

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

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

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

［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］高睿，任冰，王國(guó)玉，等.孤立波淺化過(guò)程的SPH數(shù)值模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展: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］劉謀斌，宗智，常建忠.光滑粒子動(dòng)力學(xué)方法的發(fā)展與應(yīng)用［J］.力學(xué)進(jìn)展，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方法及其對(duì)二維非線性水波的模擬［J］.計(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］王亞萍，光滑粒子流體動(dòng)力學(xué)(SPH)方法在潰壩問(wèn)題中的應(yīng)用研究［D］.天津:天津大學(xué)，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算法改進(jìn)及在晃蕩與入水中的應(yīng)用［D］.大連:大連理工大學(xué)，2014.

［68］李大鳴，陳海舟，張建偉，等.基于SPH法的二維矩形液艙晃蕩研究［J］.計(jì)算力學(xué)學(xué)報(bào)，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)，龔凱，等.二維矩形水槽晃蕩過(guò)程的SPH方法模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展: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］//吳有生，劉樺，許維臨，等.第九屆全國(guó)水動(dòng)力學(xué)學(xué)術(shù)會(huì)議暨第二十二屆全國(guó)水動(dòng)力學(xué)研討會(huì)文集.北京:海洋出版社，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.