金善勤，鄭興，段文洋

(哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001)

光滑粒子水動(dòng)力學(xué)(smoothed particle hydrodynamics，SPH)方法是一種無網(wǎng)格的拉格朗日粒子算法，其目前已成為水動(dòng)力學(xué)中數(shù)值模擬的主要方法之一。1994年Monaghan［1］首先將光滑粒子水動(dòng)力學(xué)方法(SPH)應(yīng)用于自由表面流動(dòng)的模擬，基于弱可壓縮流體假設(shè)引入狀態(tài)方程，以二維潰壩模型為研究對(duì)象，討論了SPH方法在計(jì)算水動(dòng)力學(xué)問題中的可行性。隨后，Colagrossi等［2-3］采用 SPH方法，對(duì)任意形狀物體入水后的自由表面形狀及強(qiáng)非線性效應(yīng)進(jìn)行了研究。

傳統(tǒng)SPH方法在壓力計(jì)算和物理粘性模擬等方面存在不足，國(guó)內(nèi)外許多學(xué)者對(duì)其進(jìn)行了改進(jìn)。Monaghan等［4］為了減少粒子在運(yùn)動(dòng)中出現(xiàn)的非物理震蕩，防止粒子在運(yùn)動(dòng)過程中由于距離太小而發(fā)生非物理穿透，其在動(dòng)量方程中加入了人工粘性修正項(xiàng)。針對(duì)弱可壓縮SPH方法在壓力求解過程中出現(xiàn)的不穩(wěn)定現(xiàn)象，Marrone等［5-6］基于黎曼解的壓力修正思想提出了δ-SPH方法，并以二維潰壩問題為研究對(duì)象，證實(shí)了該方法的有效性。國(guó)內(nèi)的許多學(xué)者也在SPH算法改進(jìn)方面進(jìn)行研究，高睿等［7］采用黎曼求解N-S方程，改進(jìn)了SPH方法，并研究了規(guī)則波與結(jié)構(gòu)物的相互作用，鄭興等［8-9］針對(duì)傳統(tǒng)SPH方法二階核近似后存在精度不高的問題，提出了改進(jìn)的K2_SPH方法，并將其應(yīng)用于強(qiáng)非線性自由表面問題的研究。

由于傳統(tǒng)SPH方法計(jì)算效率不高，提高 SPH計(jì)算效率也是研究人員關(guān)注的一個(gè)重點(diǎn)。采用GPU加速SPH算法是近幾年來發(fā)展的一種新的SPH方法計(jì)算模式。借助于GPU的眾核優(yōu)勢(shì)，對(duì)SPH算法并行優(yōu)化，能夠顯著的提高計(jì)算效率。目前國(guó)內(nèi)外許多研究人員已經(jīng)在該方面取得了一些成果，例如Xiong等［10］采用自適應(yīng)粒子分裂和合并技術(shù)在單GPU和多GPU上分別實(shí)現(xiàn)了SPH方法的并行計(jì)算，徐鋒等［11］對(duì)SPH方法在GPU上的物理存儲(chǔ)模式和線程使用方式進(jìn)行了較為詳細(xì)的論述，朱小松等［12］對(duì)GPU和CPU上的MPS方法和SPH方法的并行計(jì)算進(jìn)行研究，并借助此方法對(duì)液艙晃蕩問題進(jìn)行了分析。雖然目前基于GPU的SPH算法發(fā)展迅速，但是仍然存在一些問題有待解決，如并行計(jì)算模式單一、壓力計(jì)算結(jié)果不穩(wěn)定等。

針對(duì)以上不足，本文提出一種基于粒子對(duì)的并行計(jì)算方法，并將其與δ-SPH方法結(jié)合，通過對(duì)粘性流場(chǎng)的模擬，研究不同邊界處理方法對(duì)計(jì)算結(jié)果的影響，不同粒子加速方法對(duì)計(jì)算精度和加速效果的影響。

1 SPH的基本理論

SPH方法是一個(gè)在粒子系統(tǒng)中進(jìn)行插值的計(jì)算方法，流域被離散成許多的有限體積單元，這些有限體積單元即為粒子。每個(gè)粒子點(diǎn)都具有自身的物理屬性如:質(zhì)量、密度、速度等。在i粒子處的場(chǎng)函數(shù)和場(chǎng)函數(shù)的梯度可用核近似的形式來表示:

式中:h為光滑長(zhǎng)度，N為i粒子支持域內(nèi)的相鄰粒子數(shù)，fj表示j粒子點(diǎn)的場(chǎng)函數(shù)值，Vj表示j粒子的體積，W ri-rj，h

( )表示核函數(shù)值。當(dāng)支持域半徑趨向于0時(shí)，核函數(shù)W ri-rj，h( )便變成了狄拉克函數(shù)。

2 δ-SPH的控制方程

δ-SPH方法是在傳統(tǒng)SPH方法的基礎(chǔ)上，基于密度修正和黎曼解的壓力修正思想提出的一種改進(jìn)SPH方法，其控制方程可描述為

式中:p、ρ、v分別表示粒子點(diǎn)的壓力、密度、速度;f為外力項(xiàng);σ為總應(yīng)力張量，由壓力p和粘性應(yīng)力τ組成;α、β表示坐標(biāo)方向。

對(duì)于牛頓流體，粘性剪切應(yīng)力和剪切應(yīng)變成正比，于是有

式(3)中，Ψij和πij分別為密度修正項(xiàng)和人工粘性項(xiàng)，h為平均光滑長(zhǎng)度，c為平均聲速，ω和ε分別為密度修正項(xiàng)和人工粘性項(xiàng)系數(shù)。Ψij表達(dá)式為

其中，

采用SPH方法模擬粘性流場(chǎng)的流動(dòng)既是一個(gè)邊值問題，又是一個(gè)初始條件問題。采用不同的邊界處理方法，會(huì)對(duì)計(jì)算的精度和穩(wěn)定性造成不同的影響。所以，有必要對(duì)邊界處理方法進(jìn)行相應(yīng)的研究，尋找最優(yōu)的邊界處理方法。目前常用的邊界處理方法有:基于排斥力的固壁粒子方法、鏡像粒子邊界處理方法和周期性邊界處理方法，關(guān)于其詳細(xì)介紹請(qǐng)參考文獻(xiàn)［13］。

3 基于粒子對(duì)的GPU并行算法

借助GPU的眾核架構(gòu)優(yōu)勢(shì)，加速SPH算法是近年來發(fā)展的一種全新的計(jì)算方法。傳統(tǒng)SPH方法良好的并行特性，為該方法能夠順利在GPU上的實(shí)現(xiàn)并行計(jì)算提供了基礎(chǔ)。此外相比于完全不可壓縮SPH(ISPH)，本文采用的δ-SPH方法與傳統(tǒng)SPH方法一樣具有良好的并行特性。目前基于GPU的SPH算法大多采用基于單個(gè)粒子的并行計(jì)算模式，單個(gè)線程計(jì)算單個(gè)粒子的流體控制方程，得到單個(gè)粒子點(diǎn)的密度變化率和加速度等信息。在GPU的單指令多線程(SIMT)并行模式下，多個(gè)線程可同時(shí)發(fā)射、同時(shí)運(yùn)行，由于SPH方法各個(gè)粒子之間不存在較強(qiáng)的相關(guān)性，這樣多個(gè)粒子可同時(shí)在GPU上開展并行計(jì)算，關(guān)于該方面的詳細(xì)介紹可參考文獻(xiàn)［11，14］。研究表明隨著粒子數(shù)量的不斷增加，采用基于單個(gè)粒子并行方法的加速能力提高緩慢。為此，本文提出了一種基于粒子對(duì)的并行計(jì)算方法，該方法的全部計(jì)算任務(wù)都在GPU上完成，中間計(jì)算過程不需要進(jìn)行GPU和CPU的數(shù)據(jù)傳遞，因此該方法的計(jì)算效率得到了較大提高。圖1給出了采用基于單個(gè)粒子和粒子對(duì)2種不同并行計(jì)算方法的加速比隨著粒子數(shù)量變化的曲線。由圖1可得，采用基于粒子對(duì)的并行計(jì)算方法優(yōu)勢(shì)明顯。傳統(tǒng)SPH計(jì)算方法大多采用逐個(gè)求解單個(gè)粒子控制方程的計(jì)算模式，Riffert等［15］提出了成對(duì)相互作用法，通過分析粒子的成對(duì)相互作用，可以較大的提高粒子的計(jì)算效率和存儲(chǔ)空間。

圖1 加速比隨粒子數(shù)增加的變化曲線Fig.1 Speed ratio curve along with increase of the particle number

在CUDA架構(gòu)下，每個(gè)線程塊中的線程都有系統(tǒng)賦予其唯一的標(biāo)識(shí)符［16］。雖然不同計(jì)算能力的GPU所能使用的最大線程數(shù)量不同，當(dāng)在所有流處理器都被占滿時(shí)，由于計(jì)算指令處于等待模式，這樣使得粒子數(shù)量再多，同樣可以使得每個(gè)線程和粒子對(duì)編號(hào)一一對(duì)應(yīng)。計(jì)算開始時(shí)，由CPU主機(jī)端將流場(chǎng)粒子的初始物理信息全部拷貝至GPU設(shè)備端，根據(jù)δ-SPH的求解順序，每個(gè)線程控制對(duì)應(yīng)粒子依次進(jìn)行粒子對(duì)查找，密度計(jì)算，修正項(xiàng)密度計(jì)算，加速度計(jì)算等。然后采用跳蛙法求出下一步的流場(chǎng)初始信息，不斷循環(huán)至流場(chǎng)達(dá)到穩(wěn)定狀態(tài)，基于粒子對(duì)的并行計(jì)算流程如圖2所示。

在GPU上大量的晶體管被用于執(zhí)行單元，而不像CPU那樣作為控制單元和緩存，因此GPU對(duì)選擇、循環(huán)等分支結(jié)構(gòu)的處理能力并不如CPU那樣出色。在SPH算法中，邊界處理和鄰近粒子的查找需要許多分支語句，因此，這類操作在整個(gè)并行過程中耗時(shí)最多。為了提高計(jì)算效率，本文2種并行算法均采用基于網(wǎng)格哈希值的鏈表搜索方法，關(guān)于其詳細(xì)介紹請(qǐng)參考文獻(xiàn)［14］。

采用粒子對(duì)并行計(jì)算方法時(shí)，多個(gè)線程在不同流處理器中對(duì)相應(yīng)的粒子對(duì)進(jìn)行計(jì)算，此時(shí)會(huì)存在多個(gè)線程同時(shí)對(duì)某一粒子的變量?jī)?nèi)存進(jìn)行讀取和寫入操作的情況，這樣會(huì)導(dǎo)致數(shù)據(jù)疊加錯(cuò)誤，必須對(duì)線程操作進(jìn)行排序，才能得到正確的計(jì)算結(jié)果。為了保證每個(gè)線程操作結(jié)果的正確，本文采用了原子操作。原子操作能夠控制對(duì)同一變量?jī)?nèi)存進(jìn)行讀取寫入操作的線程數(shù)量，使得基于粒子對(duì)并行模式的GPU計(jì)算方法得以實(shí)現(xiàn)。

圖2 基于粒子對(duì)的GPU并行流程Fig.2 GPU parallel process based on particle pair

4 數(shù)值模擬和分析

4.1 多種邊界處理方法的比較

4.1.1 腔內(nèi)剪切流動(dòng)問題數(shù)值模型

腔內(nèi)剪切流動(dòng)是一種經(jīng)典的流體運(yùn)動(dòng)，流體放置在正方形容器內(nèi)，容器的頂部平板以恒定速度vtop做水平運(yùn)動(dòng)，容器的其余三邊為固定邊界，通過頂部邊界的運(yùn)動(dòng)來驅(qū)動(dòng)容器內(nèi)流體的運(yùn)動(dòng)，流場(chǎng)穩(wěn)定時(shí)能在靠近頂部處形成回流［17-18］。這一計(jì)算模型能夠較好地檢驗(yàn)SPH算法對(duì)運(yùn)動(dòng)邊界和固定邊界的處理。計(jì)算數(shù)值模型相關(guān)參數(shù)為:正方形容器邊長(zhǎng)l=0.001 m ，頂部邊界速度vtop=0.001 m/s，動(dòng)粘系數(shù)ν=10-6m2/s，密度ρ=103kg/m3。在此情況下雷諾數(shù)為1。在正方形容器內(nèi)分布有1 600(40×40)個(gè)粒子，分別采用基于排斥力的固壁粒子邊界處理法、鏡像粒子處理法、固壁粒子加鏡像混合處理方法3種不同方法對(duì)正方形邊界進(jìn)行處理。

4.1.2 腔內(nèi)剪切流數(shù)值結(jié)果與分析

在模擬腔內(nèi)剪切流動(dòng)問題時(shí)，壓力狀態(tài)方程的c0取0.01，固壁粒子半徑取初始粒子間距的一半。圖3為采用不同邊界處理方法，流場(chǎng)右上角粒子的位置及速度分布情況。

為了驗(yàn)證模擬結(jié)果，圖 4 給出與 G.R.Liu［19］中有限差分結(jié)果對(duì)比圖。由圖4可知，垂向粒子點(diǎn)的水平速度和FDM的計(jì)算結(jié)果吻合較好，但是在水平中心線處的垂向速度分布上，δ-SPH方法與FDM計(jì)算結(jié)果存在偏差，均存在回流中心左移的現(xiàn)象。從計(jì)算效率和結(jié)果精度兩方面考慮，鏡像粒子點(diǎn)法完全能夠滿足粘性流場(chǎng)的需求。

為了研究GPU并行程序隨粒子數(shù)量增加的加速效果及數(shù)值結(jié)果的收斂性，本文采用鏡像粒子邊界處理方法分別對(duì)流域內(nèi)布置80×80個(gè)粒子點(diǎn)及100×100個(gè)粒子點(diǎn)2種計(jì)算模型進(jìn)行了計(jì)算。圖5(a)中給出6 400(80×80)個(gè)粒子點(diǎn)的最終流場(chǎng)分布。在靠近正方形頂部邊界處可以看見明顯的回流。圖5(b)為不同粒子點(diǎn)使水平中心線處無量綱垂向速度分布。由圖可知隨著流場(chǎng)內(nèi)實(shí)粒子數(shù)量的增加，水平中心線處的垂向速度基本相同。

圖3 不同邊界處理方法的最終流場(chǎng)分布Fig.3 Final flow pattern of different boundary processing

圖4 腔內(nèi)剪切流流場(chǎng)水平中心線處無量綱垂向速度分布及垂向中心線處無量綱水平速度Fig.4 Dimensionless vertical velocity of horizontal center line and dimensionless horizontal velocity of vertical center line in cavity shear flow

圖5 6 400個(gè)粒子計(jì)算模型最終流場(chǎng)及不同粒子點(diǎn)計(jì)算模型水平中心線處無量綱垂向速度分布Fig.5 Eventually flow pattern of model with 6 400 particles and dimensionless horizontal velocity on vertical center line of different particle number module

4.2 Poiseuille 流動(dòng)模擬

流體在兩塊無限長(zhǎng)的固定平行板(y=0和y=l)之間的流動(dòng)稱之為Poiseuille流動(dòng)。流域內(nèi)流體在平行于x軸的外力F作用下在平板間逐漸流動(dòng)并最終達(dá)到穩(wěn)態(tài)。Morris等［20］給出了Poiseuille流動(dòng)隨時(shí)間(t＞0)變化的解析表達(dá)式:

動(dòng)粘系數(shù)ν=10-6m2/s，傳動(dòng)力F=2×10-4m/s2，計(jì)算域l=10-3m，密度ρ=103kg/m3。通過解析式可知流場(chǎng)中的最大速度為2.5×10-5m/s，相應(yīng)的雷諾數(shù)Re=2.5×10-2。流場(chǎng)計(jì)算域?yàn)?0.000 5 m×0.001 m，粒子點(diǎn)數(shù)量20×39，在上下邊界處各布置20個(gè)粒子點(diǎn)，具體細(xì)節(jié)可參考文獻(xiàn)［9］。計(jì)算選用五次樣條核函數(shù)，時(shí)間步長(zhǎng) dt=0.000 1 s。表1 給出t=0.6 s時(shí)整個(gè)流場(chǎng)左側(cè)第1列39個(gè)粒子的誤差比較。表1的計(jì)算結(jié)果表明，δ-SPH方法在CPU上運(yùn)算結(jié)果的誤差范圍大約是0.1%～0.5%。在 GPU的最大誤差不超過1.5%。圖6為在GPU和CPU平臺(tái)上不同時(shí)刻Poiseuille流動(dòng)水平方向速度和加速度分布情況。

表1 t=0.6 s GPU和CPU與Poiseuille流動(dòng)的解析解速度結(jié)果比較Table 1 Velocity comparison of GPU and CPU with analytical solutions when t=0.6 s for Poiseuille flow

圖6 不同時(shí)刻Poiseuille流動(dòng)速度和加速度分布Fig.6 Velocity and acceleration of Poiseuille flow at different moments

4.3 Couette 流動(dòng)模擬

Couette流動(dòng)是t=0時(shí)刻，上層平板突然以一個(gè)平行于x軸的恒定速度v0運(yùn)動(dòng)。Morris等［20］給出其流場(chǎng)速度隨時(shí)間變化的解析表達(dá)式為

計(jì)算域和 Poiseuille 流動(dòng)，v0=2.5×10-5m/s，雷諾數(shù)為Re=2.5×10-2，采用五次樣條核函數(shù)，時(shí)間步長(zhǎng)dt=0.000 1 s。表 2 給出t=0.6 s時(shí)，采用 δ-SPH 方法在CPU和GPU上計(jì)算所得x方向計(jì)算速度與解析解的誤差分布比較。圖7為在GPU和CPU平臺(tái)上不同時(shí)刻Couette流動(dòng)的速度和加速度。通過采用δ-SPH方法在CPU和GPU上對(duì)低雷諾數(shù)的Poiseuille流動(dòng)和Couette流動(dòng)的模擬結(jié)果得到，在GPU上和CPU上計(jì)算結(jié)果沒有太大差異，沒有因GPU流處理器在表示浮點(diǎn)數(shù)方面能力較差而帶來精度損失。

表2 t=0.6 s GPU和CPU與Coeutte流動(dòng)的解析解速度結(jié)果比較Table 2 Velocity comparison of GPU and CPU with analytical solutions when t=0.6 s for Coeutte flow

圖7 不同時(shí)刻Couette流動(dòng)速度和加速度分布Fig.7 Velocity and acceleration distribution of Couette flow at different moments

4.4 孤立波砰擊模擬

為了進(jìn)一步驗(yàn)證本文方法的通用性，本節(jié)還給出孤立波砰擊計(jì)算結(jié)果。其計(jì)算模型如圖8所示。其中，水槽長(zhǎng)10 m，水深0.25 m。在水槽末端壁面上距水池底部0.05 m及0.15 m處布置P0、P12個(gè)壓力測(cè)點(diǎn)，在距離水槽左端2 m和8 m處分別布置浪高儀。粒子間距dx=0.01 m，粒子數(shù)為25 000(25×1 000)。其中，水槽末端粒子分布如圖9所示。

模擬過程中造波板運(yùn)動(dòng)方程采用Ma［21］所給出的計(jì)算公式:

圖8 孤立波砰擊計(jì)算模型Fig.8 Calculating model of the solitary wave slamming

圖9 2 m處孤立波波形Fig.9 Wave elevation of solitary wave at 2 m

圖10 P0和P1點(diǎn)壓力模擬結(jié)果與實(shí)驗(yàn)值Fig.10 Simulation pressure results and experimental pressure value at point P0and P1

從圖9和圖10可以看出，采用δ-SPH方法可以準(zhǔn)確的模擬孤立波，數(shù)值模擬結(jié)果與實(shí)驗(yàn)值吻合十分完好。從P0與P1點(diǎn)的壓力變化情況可以看出δ-SPH方法有效的減少了壓力的非物理振蕩，壓力模擬結(jié)果與實(shí)驗(yàn)值誤差較少，但是在孤立波與水槽末端直壁碰撞過后的壓力模擬結(jié)果仍然與實(shí)際存在偏差。

4.5 GPU加速性能分析

與傳統(tǒng)SPH方法相比，δ-SPH方法增加了密度修正項(xiàng)，需要求出修正后的密度梯度，其中涉及大量的矩陣運(yùn)算，因此在計(jì)算時(shí)間上付出了一定的代價(jià)。本文中采用基于單個(gè)粒子和基于粒子對(duì)2種不同的并行設(shè)計(jì)方案實(shí)現(xiàn)了δ-SPH的GPU并行計(jì)算，所有的并行計(jì)算程序都是在NVIDIA的CUDA環(huán)境中開發(fā)，版本為5.5，模擬在一臺(tái)PC機(jī)上運(yùn)行，其配置如下:Intel Core(TM)i5-3470 CPU 3.20GHZ，4.0 GB內(nèi)存和GTX660顯卡，內(nèi)存2G。

Couette流動(dòng)和Poiseuille流動(dòng)采用的是類似的計(jì)算模型，整個(gè)流場(chǎng)的流域?qū)嵙Ｗ?、邊界鏡像粒子、周期性粒子共計(jì)1 222個(gè)粒子。采用鏡像法進(jìn)行邊界處理的腔內(nèi)剪切流動(dòng)粒子數(shù)分別為:2 116、7 396、11 236個(gè)粒子。孤立波砰擊模擬中整個(gè)流場(chǎng)的粒子總計(jì)28 168個(gè)粒子點(diǎn)。表3中的數(shù)據(jù)為各計(jì)算模型采用不同計(jì)算模式的單步時(shí)間消耗。

表3 不同計(jì)算模式的單步時(shí)間消耗Table 3 Single step time consumption with different calculation models

圖11為各計(jì)算模型不同計(jì)算模式的時(shí)間消耗柱狀圖。從圖中比較可得，不論是在CPU平臺(tái)還是在GPU平臺(tái)，采用粒子對(duì)進(jìn)行計(jì)算總能得到很好的計(jì)算效率。隨著粒子數(shù)的增加，GPU加速效果不斷提高。采用粒子對(duì)并行模式的GPU計(jì)算方法明顯優(yōu)于傳統(tǒng)的基于單個(gè)粒子的并行計(jì)算模式，雖然在計(jì)算的過程中大量的使用了原子操作，但隨著粒子數(shù)量的增加，計(jì)算量的迅速增加有效的掩蓋了原子操作所造成的時(shí)間等待。

圖11 不同計(jì)算模式的時(shí)間消耗柱狀圖Fig.11 Histograms of time consumption of different calculation

5 結(jié)論

本文對(duì)δ-SPH方法的GPU并行計(jì)算方法進(jìn)行了研究，將其應(yīng)用于對(duì)粘性流場(chǎng)問題的模擬。通過對(duì)相關(guān)數(shù)值模擬可以得到以下結(jié)論:

1)δ-SPH方法在模擬粘性流場(chǎng)問題時(shí)，采用傳統(tǒng)的鏡像粒子法對(duì)邊界進(jìn)行處理，由腔內(nèi)剪切流的計(jì)算結(jié)果表明該方法的結(jié)果精度上較傳統(tǒng)SPH有所提高。

2)在大規(guī)模計(jì)算時(shí)，采用基于粒子對(duì)模式的GPU并行計(jì)算方法比基于單個(gè)粒子的GPU并行計(jì)算方法的加速性能要好，其計(jì)算速度更是CPU串行計(jì)算不可比擬的。計(jì)算結(jié)果并不會(huì)因GPU流處理器在表示浮點(diǎn)數(shù)方面能力較差而帶來精度損失。

3)隨著粒子數(shù)量的增加，基于粒子對(duì)的GPU并行計(jì)算方法的加速效果增加迅速，體現(xiàn)了其在并行方面的強(qiáng)大優(yōu)勢(shì)。

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

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

［3］COLAGROSSI A，LUGNI C，BROCCHINI M.A study of violent sloshing wave impacts using an improved SPH method［J］.Journal of Hydraulic Research，2010，48(S1):94-104.

［4］MONAGHAN J J.Particle methods for hydrodynamics［J］.Computer Physics Reports，1985，3(2):71-124.

［5］MARRONE S，ANTUONO M，COLAGROSSI A，et al.δ-SPH model for simulating violent impact flows［J］.Computer Methods in Applied Mechanics and Engineering，2011，200(13):1526-1542.

［6］MOLTENI D，COLAGROSSI A.A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH［J］.Computer Physics Communications，2009，180(6):861-872.

［7］高睿.SPH強(qiáng)非線性水動(dòng)力學(xué)數(shù)值模型的改進(jìn)與應(yīng)用［D］.大連:大連理工大學(xué)，2011:20-78.GAO Rui.Correction and application of high nonlinear SPH hydrodynamic model［D］.Dalian:Dalian University of Technology，2011:20-78.

［8］鄭興，段文洋.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.

［9］ZHENG Xing，DUAN Wenyang.Numerical simulation of dam breaking using smoothed particle hydrodynamics and viscosity behavior［J］.Journal of Marine Science and Application，2010，(9):34-41.

［10］XIONG Qingang，LI Bo，XU Ji.GPU-accelerated adaptive particle splitting and merging in SPH［J］.Computer Physics Communications，2013，184(7):1701-1707.

［11］徐鋒.基于眾核架構(gòu)的并行 SPH算法的研究與實(shí)現(xiàn)［D］.上海:上海交通大學(xué)，2013:35-65.XU Feng.Research and implementation of the smoothed particle hydrodynamics algorithm based on multi-core architecture［D］.Shanghai:Shanghai Jiao Tong University，2013:35-65.

［12］朱小松.粒子法的并行加速及在液體晃蕩研究中的應(yīng)用［D］.大連:大連理工大學(xué)，2011:74-109.ZHU Xiaosong.Parallel acceleration of particle method and its application on the study of liquid sloshing［D］.Dalian:Dalian University of Technology，2011:74-109.

［13］鄭興.SPH方法改進(jìn)研究及其在自由面流動(dòng)問題中的應(yīng)用［D］.哈爾濱:哈爾濱工程大學(xué)，2010:55-89.ZHENG Xing.An investigation of improved SPH and its application for free surface flow［D］.Harbin:Harbin Engineering University，2010:55-89.

［14］溫嬋娟，歐嘉蔚，賈金原.GPU通用計(jì)算平臺(tái)上的SPH流體模擬［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2010，22(3):406-411.WEN Chanjuan，OU Jiawei，JIA Jinyuan.GPU-based smoothed particle hydrodynamic fluid simulation［J］.Journal of Computer-Aided Design ＆ Computer Graphics，2010，22(3):406-411.

［15］RIFFERT H，HEROLD H，F(xiàn)LEBBE O，et al.Numerical aspects of the smoothed particle hydrodynamics method for simulating accretion disks［J］.Computer Physics Communications，1995，89(1-3):1-16.

［16］張舒，褚艷利.GPU高性能運(yùn)算之CUDA［M］.北京:中國(guó)水利水電出版社，2009:34-56.

［17］LIU G R，LIU M B，LI Shaofan.Smoothed particle hydrodynamics—a meshfree method［J］.Computational Mechanics，2004，33(6):491.

［18］LIU M B，LIU G R，LAM K Y，et al.Smoothed particle hydrodynamics for numerical simulation of underwater explosion［J］.Computational Mechanics，2003，30(2):106-118.

［19］LIU M B，LIU G R，ZONG Z，et al.Numerical simulation of incompressible flows by SPH［C］//International Conference on Scientific and Engineering Computational.Beijing，China，2001:3-6.

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

［21］MA Q W，ZHOU Juntao.MLPG-R method for numerical simulation of 2D breaking waves［J］.Computer Modeling in Engineering and Sciences，2009，43(3):277-304.