GPU加速的SPH方法在溢洪道水流模擬中的應(yīng)用

摘要：基于拉格朗日描述的光滑粒子動力學(xué)方法（SPH）擅長于處理自由面劇烈變化的水流現(xiàn)象，十分適合水利工程中泄洪等問題的數(shù)值模擬。然而，SPH方法通常采用均勻分布的粒子對流體計算域進(jìn)行空間離散，對于工程問題而言需要的粒子數(shù)量較多、計算量大。為了突破SPH方法在實際大規(guī)模計算中的適用范圍，采用C++和CUDA混合編程的技術(shù)，借助GPU實現(xiàn)了對SPH方法的并行加速。通過WES三圓弧段組成的光滑溢洪道過流問題，驗證了GPU加速的SPH方法的計算精度和可靠性，計算效率相對原始的SPH仿真過程提高了61.8倍。最后，將GPU加速的SPH方法應(yīng)用于水利工程的溢洪道泄流問題，分別模擬了光滑溢洪道和臺階式溢洪道流動特性，通過自由面的演化過程及泄流沿程截面上的速度分布狀態(tài)，對比分析了臺階對泄流現(xiàn)象的影響。

關(guān)鍵詞：光滑粒子動力學(xué)方法;GPU加速;臺階式溢洪道;消能率

中圖法分類號：TV512

文獻(xiàn)標(biāo)志碼：A

DOI：10.16232/j.cnki.1001-4179.2019.03.038

1研究背景

臺階式溢洪道是將傳統(tǒng)光滑溢洪道的泄流槽做成臺階式，水流在流經(jīng)臺階時與每級臺階均產(chǎn)生劇烈的碰撞，形成水流的旋滾及內(nèi)部的紊動剪切作用，促使水流表面破碎，進(jìn)而能顯著增加溢洪道的泄流消能率，有助于減小下游消力池的規(guī)模[1-2]。所以，目前臺階式溢洪道在國內(nèi)外許多工程上得以應(yīng)用[3-5]。然而，水流與臺階之間的劇烈作用將對溢洪道的安全提出挑戰(zhàn)，有必要對臺階式溢洪道的流動現(xiàn)象及機(jī)理進(jìn)行系統(tǒng)研究，以確保溢洪設(shè)施的安全。盡管世界各國水利工程技術(shù)人員對此流動問題開展了大量的試驗研究[6]，但受限于尺度效應(yīng)、測量手段、經(jīng)費(fèi)等原因，難以獲取流動機(jī)理分析所必備的詳細(xì)數(shù)據(jù)。相較之下，采用數(shù)值模擬的手段對該泄流問題進(jìn)行研究能夠得到豐富的可視化結(jié)果，故而深受研究人員的青睞[7-8]。

在既往的研究中，水利工作者多采用以有限體積法為代表的網(wǎng)格類方法開展溢洪過程中流場的數(shù)值模擬。此類方法雖然通過與VOF、Level-set等自由面捕捉技術(shù)的結(jié)合，實現(xiàn)了對一系列水力學(xué)問題的模擬，但仍難以真實再現(xiàn)自由面的翻卷、破碎、飛濺等劇烈流動現(xiàn)象。近年來具有拉格朗日特性的SPH方法在自由表面流動模擬中取得了較大進(jìn)展[9-10]，能夠較為真實地捕捉自由面的翻卷、破碎等非線性現(xiàn)象，并成功應(yīng)用于入水沖擊[11]、液艙晃蕩[12]、波浪與結(jié)構(gòu)物的相互作用等復(fù)雜流動問題[13]。由于臺階式溢洪道的泄流問題與此類問題具有相似的特性，例如水流的翻滾、碰撞等，故而本文嘗試采用SPH方法對臺階式溢洪道的泄洪過程進(jìn)行數(shù)值模擬。

傳統(tǒng)的SPH方法是通過粒子的運(yùn)動來表示流動現(xiàn)象，為了精確模擬強(qiáng)非線性變化的自由面，需要采用大量的粒子對流體空間進(jìn)行離散，同時需要采用較小的計算時間步長，故而對計算資源的消耗大，計算效率通常不高。近年來，GPU（GraphicsProcessingUnit，圖形處理器）硬件的計算核心數(shù)和存儲能力迅速提升，國外一些學(xué)者逐步采用基于GPU加速的技術(shù)來提高SPH方法的計算效率。例如，Crespo等[14]，Xia和Liang[15]，Mokos等[16]分別采用CUDA語言編寫程序調(diào)用GPU硬件環(huán)境實現(xiàn)對SPH方法的加速計算，這些研究成果表明GPU對SPH方法計算效率的提高有明顯的幫助。然而，國內(nèi)在SPH方法中采用GPU加速技術(shù)的相關(guān)研究較少，本文將在此方面作初步嘗試，并應(yīng)用于水利工程中的溢洪道流動問題研究。

2SPH數(shù)值方法

2.1控制方程

本文基于拉格朗日描述的SPH方法進(jìn)行流場數(shù)值模擬.相應(yīng)流體控制方程為

公式

式中，ρ為流體密度，P為壓力，V為速度向量，g為重力加速度向量，y是運(yùn)動黏性系數(shù)。式（1）和式（2）的時間導(dǎo)數(shù)項是以物質(zhì)導(dǎo)數(shù)的形式給出的。在粒子法中，粒子的位置和其他物理量都是基于拉格朗日描述法表達(dá)的，因此不需要計算對流項。

2.2粒子作用模型

2.2.1核函數(shù)

在SPH方法中，需要借助核函數(shù)對空間任意場函數(shù)f（r）進(jìn)行近似積分表示：

公式

式中，r為任意粒子的空間矢量，r'為目標(biāo)粒子的空間矢量，Ω為r的積分域，h為積分域的光滑半徑，W（r-r'，h）為核函數(shù)。本文采用的核函數(shù)表達(dá)式如下：

公式

式中，q=r/h。

2.2.2密度模型

在SPH方法中，流體是弱可壓的，流體粒子的密度可通過其作用域內(nèi)所有粒子的密度作加權(quán)平均得到，這里的加權(quán)函數(shù)即為，上述核函數(shù)：

公式

式中，pi為粒子i的密度，mj為粒子j的質(zhì)量，rij為粒子i與j之間的距離。

2.2.3梯度模型

由于控制方程中存在壓力梯度項，本文SPH方法采用粒子模型形式的梯度表達(dá)式如下

公式

式中，Pi為粒子i的壓力。

2.3狀態(tài)方程

為實現(xiàn)流場壓力的求解，本文SPH方法采用了描述弱可壓流體的狀態(tài)方程建立粒子壓力與流體密度間的關(guān)系：

公式

式中，ρo為參考密度，常數(shù)γ=7，系數(shù)B=ρoc2/γ，聲速c通常取流場最大速度的10倍。

3GPU加速在SPH方法中的實現(xiàn)

GPU在設(shè)計之初是為了滿足具有高度并行化特征的圖像處理工作，相較于CPU具備更多的計算核心。例如，普通的臺式電腦配備的CPU擁有8個計算核心，而配備的圖形顯卡則擁有，上千個計算核心，這些計算核心能并發(fā)執(zhí)行較多的任務(wù)，從而能使一臺普通電腦完成一臺高性能工作站的計算任務(wù)。

為了充分利用GPU設(shè)備的并行計算能力，需要對SPH計算程序進(jìn)行改進(jìn)，使之能夠同時調(diào)用CPU和GPU的異構(gòu)硬件資源。本文在SPH程序的實現(xiàn)過程中，采用C++語言編寫主流程，采用CUDA語言（ComputeUnifiedDeviceArchitecture，通用并行計算架構(gòu)）調(diào)用GPU設(shè)備并發(fā)執(zhí)行計算任務(wù)，改造后的SPH計算流程如圖1所示。對于SPH方法，鄰居粒子的搜尋和粒子間相互作用力的計算是整個流程中計算量最大的模塊部分，由于該部分是顯示計算過程，非常適合采用并發(fā)任務(wù)執(zhí)行策略，故而本文將該部分計算任務(wù)分配于GPU設(shè)備端。同時，為了減少GPU設(shè)備和CPU之間數(shù)據(jù)通信帶來的計算延遲現(xiàn)象，本文將每個時間步內(nèi)粒子位置的更新任務(wù)同樣分配于GPU設(shè)備端。

4GPU加速的SPH方法數(shù)值驗證

在采用SPH方法數(shù)值模擬臺階式溢洪道的流動問題之前，首先采用GPU加速對該方法的精度和效率進(jìn)行驗證。本節(jié)對WES曲線（Xxl.85=2.0Hd0.85Y）表達(dá)的三圓弧光滑式溢流道的過流問題進(jìn)行模擬，并基于Michels的實驗數(shù)據(jù)對本文GPU加速的SPH方法進(jìn)行驗證[17]。

本文采用的WES光滑式溢流道模型如圖2所示。仿真時水的密度ρ=1000kg/m3，流體的運(yùn)動黏性系數(shù)v=1.01x10-6m2/s，重力加速度g=9.81m/s2。采用粒子對圖中模型進(jìn)行空間離散，初始粒子間距l(xiāng)為0.04m，模型粒子總數(shù)為42.15萬，計算時間步長為0.0005s。本文的仿真工作均采用表1所示的計算環(huán)境。

圖3為利用SPH方法模擬得到的溢洪道上各水平位置處流體粒子的最大速度與實驗測定流速的對比。其中，L表示相對壩頂?shù)乃骄嚯x，h為水面相對溢洪道的垂直距離。由圖可知，在整個溢洪道上，SPH方法得到的流速與實驗結(jié)果基本一致。基于此，SPH方法亦能夠獲得與實驗一致的溢流水面形態(tài)，二者結(jié)果的對比如圖4所示。通過本文結(jié)果與Michels實驗結(jié)果的對比可知，SPH方法在處理此類自由面問題時具有較好的可靠性，可進(jìn)一步應(yīng)用于臺階式溢洪道流動等水利問題的數(shù)值模擬研究

本文分別在不調(diào)用和調(diào)用GPU加速的工況下開展仿真模擬，兩種工況下計算5000個時間步的機(jī)器運(yùn)行時間如圖5所示，在SPH方法中采用GPU加速技術(shù)的計算效率相對無加速時提高了61.8倍。由此可見，在SPH方法不但能夠?qū)σ绾榈佬沽鲉栴}進(jìn)行精確地模擬，通過結(jié)合GPU加速技術(shù)，還能夠顯著提高該方法的數(shù)值仿真效率，有助于SPH方法在工程問題中的推廣應(yīng)用。

5臺階式溢洪道泄流問題數(shù)值模擬

截至目前，水利工作者對溢洪道的泄流問題開展了較多的數(shù)值研究，但受到計算規(guī)模和效率的限制。采用SPH方法對該問題研究時，通常在蓄水池處分配了較小的計算域，難以保持泄洪流動狀態(tài)的穩(wěn)定性。本文采用的GPU加速技術(shù)有助于提高SPH方法的計算規(guī)模和效率，故而可采用較寬的計算域?qū)υ撔沽鲉栴}進(jìn)行數(shù)值模擬。

5.1計算模型及工況

為對比分析溢洪道泄流水動力特性，本文分別開展了光滑溢洪道和臺階式溢洪道的泄流數(shù)值模擬研究，兩種溢洪道的斜升角相同，二者幾何模型尺寸如圖6所示。模型左側(cè)為蓄水池，中部為過渡段及溢洪道，右側(cè)為出流區(qū)域，其中臺階式溢洪道共設(shè)置了10個臺階，每級臺階的長和高分別為0.1m和0.05m。采用初始間距為0.005m的粒子對圖6幾何模型進(jìn)行空間離散，SPH方法的計算工況參數(shù)為：流體密度為1000kg/m3，運(yùn)動黏性系數(shù)為1.01x10-6m2/s;粒子間距為0.0025m，粒子總數(shù)為275529，光滑長度系數(shù)為1.0，時間步長為0.00005s，模擬時長為6s。

5.2數(shù)值結(jié)果及分析

本文首先采用GPU加速的SPH方法對光滑式溢洪道的溢流過程進(jìn)行數(shù)值模擬，并與無加速時的計算速度進(jìn)行了對比，兩者運(yùn)行的硬件環(huán)境如表1所示，運(yùn)行時間的對比如圖7所示。由該圖可見，在GPU加速狀態(tài)下對27.5萬粒子數(shù)規(guī)模的溢流問題進(jìn)行仿真，計算運(yùn)行時間相較無加速工況下單個CPU核心計算時間提高了73倍。

圖8為臺階式溢洪道和光滑溢洪道的自由水面演化過程，包括開始泄流的瞬間、水流與溢洪道相互作用的初始階段、水流在溢洪道上的穩(wěn)定泄流階段等。由圖8可見，泄流開始時刻（t=0.5s），兩種溢洪道的流態(tài)高度相似，上游水平堰頂?shù)乃罹鶠?.11m，水頭越過水平堰頂?shù)挠叶讼蛞绾榈佬沽?。在t=1.0s時刻，光滑溢洪道上的水流表面較為光順;在臺階式溢洪道上，水流因與臺階的碰撞而產(chǎn)生明顯的翻卷、融合現(xiàn)象，自由面的形狀粗糙。在t=2.0s和t=3.6s時，光滑溢洪道上水流的速度和自由面形態(tài)基本相似，泄流過程中溢洪道的底部水流速度相對，上部明顯增高，達(dá)到3.5m/s。該階段的泄流過程得到充分發(fā)展，流動狀態(tài)穩(wěn)定。在t=2.0s時刻，更多的水流粒子沖擊到臺階上，并因局部的水流翻卷而形成大量的粒子聚集。在t=3.6s時刻，臺階式溢洪道上水流狀態(tài)達(dá)到穩(wěn)定狀態(tài)，并形成水流截面上的流速分層現(xiàn)象，即臺階近壁面上流速較低，水流表面處速度較高，但溢洪道上水流最大速度為3.29m/s，低于光滑溢洪道上的流動速度。

此外，在流動的穩(wěn)定階段，臺階式溢洪道截面，上水面厚度相對光滑溢洪道上水面厚度較厚。通過兩種溢洪道上流動狀態(tài)的定性對比可知，盡管本文模擬的溢洪道泄流高度較小，仍可觀測到階梯式溢洪道具有較好的泄洪消能作用。

圖9通過矢量圖的形式展示了臺階式溢洪道上流體粒子的運(yùn)動速度。可見，在各級臺階近壁面處均出現(xiàn)了渦旋。由圖10的局部放大矢量圖可見，臺階的整個夾角區(qū)域均由流體渦旋占據(jù)，夾角區(qū)域外側(cè)與光滑溢洪道上流動狀態(tài)相似，流體質(zhì)點以滑移流的方式向下游運(yùn)動，但在臺階，上滑移流和渦旋之間存在動量的轉(zhuǎn)換，進(jìn)而導(dǎo)致水流動能的沿程耗散。

圖11為兩種溢洪道泄流的沿程水深對比。由圖可見，兩種溢洪道上游泄流段水深高度基本一致，均約為0.081m，表明其來流狀態(tài)相同。隨著向下游的泄流過程，溢洪道上的水深逐步降低。其中，光滑溢洪道上的水深沿程降低更為迅速，在x=0.9m處水深高度約為0.0397m;臺階式溢洪道上的水深沿程降低幅度相對較小，在x=0.9m處水深高度約為0.05426m。

圖12～13分別展示了光滑溢洪道和臺階式溢洪道不同沿程截面處沿水深方向的流體粒子速度分布。其中，光滑溢洪道上沿程流速逐步增加，各截面處在水深方向上近似保持相對恒定的流速，而臺階式溢洪道上沿程流速以及各截面沿水深方向上的流速均逐步增加。溢洪道的上半段（例如x=0.1～0.3m段），近臺階處流速較低，在水深d=0.02m處流速增大至較大值，水深d>0.02m處流速近似恒定。溢洪道的下半段（例如x=0.7～0.9m段），雖然近臺階處流速較低，但在較深的水深處（d>0.04m）流速才接近穩(wěn)定，這表明在溢洪道的下半段臺階內(nèi)的渦旋對水深的影響范圍更大。

為定量比較臺階式溢洪道的泄洪消能效果，本文通過消能率對兩種溢洪道的消能效果進(jìn)行比較。消能率η的計算公式如下：

公式

式中，E為溢洪道壩頂I-I截面處水流總能量，由流經(jīng)該截面處的流體粒子動能與勢能之和組成;Ez為溢洪道壩下II-II截面處水流總能量。

表2對兩種溢洪道的消能率進(jìn)行了對比，本文采用的臺階式溢洪道雖然具有較小的壩體高度并設(shè)置了較少的臺階數(shù)量，但消能率仍達(dá)到了38.79%，是光滑式溢洪道的14.01%的消能率的2.77倍。美國Rice等[18]研究結(jié)論表明，臺階式溢洪道的能量損失是光滑溢洪道的能量損失的2～3倍，本文結(jié)果與該結(jié)論基本一致。

6結(jié)語

采用結(jié)合了GPU加速技術(shù)的SPH方法對臺階式溢洪道的泄流過程進(jìn)行了數(shù)值研究。首先，通過對WES三圓弧段組成的光滑溢洪道標(biāo)準(zhǔn)問題進(jìn)行數(shù)值模擬，驗證了SPH方法的精度和效率。數(shù)值結(jié)果表明本文GPU加速的SPH方法獲得了與已發(fā)表實驗數(shù)據(jù)相吻合的結(jié)果，同時計算效率相對原始的SPH仿真過程提高了61.8倍。隨后，對光滑溢洪道的泄流過程進(jìn)行了數(shù)值模擬。在GPU加速狀態(tài)下對27.5萬粒子數(shù)規(guī)模的泄流問題進(jìn)行仿真，計算運(yùn)行時間相較無加速工況下單個CPU核心計算時間提高了73倍，粒子數(shù)較多時GPU加速技術(shù)能夠更有效地提高SPH方法的計算效率。最后，對臺階式溢洪道的泄流過程進(jìn)行了數(shù)值研究，并與光滑溢洪道的流動現(xiàn)象進(jìn)行了對比分析。通過自由面的演化過程可見，SPH方法能夠成功捕獲泄流過程水頭與臺階的沖擊、臺階上水體的翻卷等現(xiàn)象。通過水粒子的速度矢量圖可見，泄流過程中臺階上存在明顯的渦旋，渦旋和臺階外側(cè)滑移流之間的動量交換使得泄流過程水流動能的耗散。對比了兩種溢洪道的沿程水深變化以及各沿程截面處的流速分布情況。兩種溢洪道沿程水深均有所降低，但臺階式溢洪道的降低幅度較小?？傊?，本文作為初步嘗試，實現(xiàn)了SPH方法在臺階式溢洪道流動問題的應(yīng)用，并展示了GPU加速技術(shù)對于SPH方法計算速度提高的明顯效果。

參考文獻(xiàn)：

[1]張峰，劉韓生，張為法.臺階式溢洪道滑行水流消能特性研究[J].長江科學(xué)院院報，2014，31（6）：37-40.

[2]任雨，王承恩，劉斌.階梯式溢洪道水平面上時均壓強(qiáng)試驗研究[J].人民黃河，2011，33（3）：123-124.

[3]鐘曉鳳，張法星，孫寧，等.某水庫階梯溢洪道水力特性研究[J].人民黃河，2016，38（6）：115-118.

[4]Chakib B.Numerical Computation of Inception Point Location for Flat-sloped Stepped Spillway[J].International Journal of Hydraulic Engineering，2013，2（3）：47-52.

[5]Meireles I，Matos J.Skimming flow in the nonaerated region of steppedspillways over embankment dams[J].Jourmal of Hydraulic Engineering，2009，135（8）：685-689.

[6]張志昌，曾東洋，劉亞菲.臺階式溢洪道滑行水流水面線和消能效果的試驗研究[J].應(yīng)用力學(xué)學(xué)報，2005，22（1）：30-35.

[7]高夢露，劉亞坤，孫洪亮，等.階梯式溢洪道的數(shù)值模擬研究[J].水利與建筑工程學(xué)報，2015（2）：28-32.

[8]趙相航，解宏偉，郭馨.臺階式溢洪道水力特性數(shù)值研究[J].青海大學(xué)學(xué)報，2016，34（4）：98-103.

[9]Monaghan J，Kocharyan A.SPH simulation of multi-phase flow[J].Computer Physics Communications，1995，87（1-2）：225-235.

[10]Amicarelli A Albano R，Mirauda D，et al.A smoothed particle hydrodynamics model for 3D solid body transport in free surface flows[J].Computers & Fluids，2015（116）：205-228.

[11]Oger G，Doring M，Alessandrini B，et al.Two dimensional SPH simulations of wedge water entries[J].Journal of Computational Physics，2006（213）：803-822.

[12]Cao X，Ming F，Zhang A，et al.Sloshing in a rectangular tank basedon SPH simulation[J].Applied Ocean Research，2014（29）：241-254.

[13]Altomare C，Crespo A，Dominguez J，et al.Applicability of smoothedparticle hydrodynamics for estimation of sea wave impact on coastalstructures[J].Coastal Engineering，2015（31）：1-12.

[14]Crespo A，Dominguez J，Rogers B，et al.Dualsphysics：open-sourceparallel CFD solver based on smoothed particle hydrodynamics[J].Computer Physics Communications，2015（187）：204-216.

[15]Xia X，Liang Q.A GPU-accelerated smoothed particle hydrodynamics model for the shallow water equations[J].Environmental Modeling& Software，2016（75）：28-43.

[16]Mokos A，Rogers B，Stansby P，et al.Multi-phase SPH modelling ofviolent hydrodynamics on GPUs[J].Computer Physics Communications，2015（196）：304-316.

[17]Michels V，Lovely M.Some prototype observations of air entrainedflow[C]// Proceedings of the Minuesota hydraulies convention，PartIV，1953.

[18]Rice C，Kadavy K.Model study of a roller compacted concretestepped spillway[J].Jourmal of Hydraulic Engineering，1996，122（6）：292-297.

引用本文：王巍.GPU加速的SPH方法在溢洪道水流模擬中的應(yīng)用[J].人民長江，2019，50（3）：216-221.

Application of GPU-based accelerated SPH method in flow simulation of spillway

WANG Wei

（Jilin Province Water Resource and Hydropouver Consultative Company，Changchun 130021，China）

Abstract：Smoothed particle hydrodynamics（SPH），a smooth particle dynamics method based on Lagrangian description isgood at numerical simulation of flood discharge with drastic changed free surface.However，the SPH method usually uses uniformly distributed particles to discretize the computational domain of the fluid in space.For engineering problems，the particlesnumbers and calculation complexity are large.In order to break through the application scope of SPH method in practical large-scale computing issues，the parallel acceleration of SPH method was realized by the mixed programming technology ofC++ andCUDA and the help of GPU hardware equipment.Through calculating the flow passing issue of smooth spillway composed of threearc segments of WES，the accuracy and reliability of GPU-based SPH method in this paper are verified.The computational efficiency is 61.8 times higher than that of the original SPH simulation process.Finally，the GPU-based accelerated SPH methodwas applied to the discharge simulation of hydraulic engineering，the flow characteristics of smooth spillway and stepped spillwaywere simulated respectively.By observing the variation process of the free surface and the velocity distribution along the dischargesection，the effect of the step on the discharge phenomenon is analyzed.

Key words：smoothed particle hydrodynamics;GPU acceleration;stepped spillway;energy dissipation ratio