基于GPU并行的改進SPH方法對粘性流場的模擬

金善勤，鄭興，段文洋

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

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

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

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

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

1 SPH的基本理論

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

式中:h為光滑長度，N為i粒子支持域內的相鄰粒子數，fj表示j粒子點的場函數值，Vj表示j粒子的體積，W ri-rj，h

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

2 δ-SPH的控制方程

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

式中:p、ρ、v分別表示粒子點的壓力、密度、速度;f為外力項;σ為總應力張量，由壓力p和粘性應力τ組成;α、β表示坐標方向。

對于牛頓流體，粘性剪切應力和剪切應變成正比，于是有

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

其中，

采用SPH方法模擬粘性流場的流動既是一個邊值問題，又是一個初始條件問題。采用不同的邊界處理方法，會對計算的精度和穩定性造成不同的影響。所以，有必要對邊界處理方法進行相應的研究，尋找最優的邊界處理方法。目前常用的邊界處理方法有:基于排斥力的固壁粒子方法、鏡像粒子邊界處理方法和周期性邊界處理方法，關于其詳細介紹請參考文獻［13］。

3 基于粒子對的GPU并行算法

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

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

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

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

采用粒子對并行計算方法時，多個線程在不同流處理器中對相應的粒子對進行計算，此時會存在多個線程同時對某一粒子的變量內存進行讀取和寫入操作的情況，這樣會導致數據疊加錯誤，必須對線程操作進行排序，才能得到正確的計算結果。為了保證每個線程操作結果的正確，本文采用了原子操作。原子操作能夠控制對同一變量內存進行讀取寫入操作的線程數量，使得基于粒子對并行模式的GPU計算方法得以實現。

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

4 數值模擬和分析

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

4.1.1 腔內剪切流動問題數值模型

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

4.1.2 腔內剪切流數值結果與分析

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

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

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

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

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

圖5 6 400個粒子計算模型最終流場及不同粒子點計算模型水平中心線處無量綱垂向速度分布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 流動模擬

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

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

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

圖6 不同時刻Poiseuille流動速度和加速度分布Fig.6 Velocity and acceleration of Poiseuille flow at different moments

4.3 Couette 流動模擬

Couette流動是t=0時刻，上層平板突然以一個平行于x軸的恒定速度v0運動。Morris等［20］給出其流場速度隨時間變化的解析表達式為

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

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

圖7 不同時刻Couette流動速度和加速度分布Fig.7 Velocity and acceleration distribution of Couette flow at different moments

4.4 孤立波砰擊模擬

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

模擬過程中造波板運動方程采用Ma［21］所給出的計算公式:

圖8 孤立波砰擊計算模型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點壓力模擬結果與實驗值Fig.10 Simulation pressure results and experimental pressure value at point P0and P1

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

4.5 GPU加速性能分析

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

Couette流動和Poiseuille流動采用的是類似的計算模型，整個流場的流域實粒子、邊界鏡像粒子、周期性粒子共計1 222個粒子。采用鏡像法進行邊界處理的腔內剪切流動粒子數分別為:2 116、7 396、11 236個粒子。孤立波砰擊模擬中整個流場的粒子總計28 168個粒子點。表3中的數據為各計算模型采用不同計算模式的單步時間消耗。

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

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

圖11 不同計算模式的時間消耗柱狀圖Fig.11 Histograms of time consumption of different calculation

5 結論

本文對δ-SPH方法的GPU并行計算方法進行了研究，將其應用于對粘性流場問題的模擬。通過對相關數值模擬可以得到以下結論:

1)δ-SPH方法在模擬粘性流場問題時，采用傳統的鏡像粒子法對邊界進行處理，由腔內剪切流的計算結果表明該方法的結果精度上較傳統SPH有所提高。

2)在大規模計算時，采用基于粒子對模式的GPU并行計算方法比基于單個粒子的GPU并行計算方法的加速性能要好，其計算速度更是CPU串行計算不可比擬的。計算結果并不會因GPU流處理器在表示浮點數方面能力較差而帶來精度損失。

3)隨著粒子數量的增加，基于粒子對的GPU并行計算方法的加速效果增加迅速，體現了其在并行方面的強大優勢。

［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強非線性水動力學數值模型的改進與應用［D］.大連:大連理工大學，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方法及其對二維非線性水波的模擬［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］徐鋒.基于眾核架構的并行 SPH算法的研究與實現［D］.上海:上海交通大學，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］朱小松.粒子法的并行加速及在液體晃蕩研究中的應用［D］.大連:大連理工大學，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方法改進研究及其在自由面流動問題中的應用［D］.哈爾濱:哈爾濱工程大學，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通用計算平臺上的SPH流體模擬［J］.計算機輔助設計與圖形學學報，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，FLEBBE 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高性能運算之CUDA［M］.北京:中國水利水電出版社，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，FOX 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.