GPU加速的SPH方法在溢洪道水流模擬中的應用

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

關鍵詞：光滑粒子動力學方法;GPU加速;臺階式溢洪道;消能率

中圖法分類號：TV512

文獻標志碼：A

DOI：10.16232/j.cnki.1001-4179.2019.03.038

1研究背景

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

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

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

2SPH數值方法

2.1控制方程

本文基于拉格朗日描述的SPH方法進行流場數值模擬.相應流體控制方程為

公式

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

2.2粒子作用模型

2.2.1核函數

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

公式

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

公式

式中，q=r/h。

2.2.2密度模型

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

公式

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

2.2.3梯度模型

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

公式

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

2.3狀態方程

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

公式

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

3GPU加速在SPH方法中的實現

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

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

4GPU加速的SPH方法數值驗證

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

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

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

本文分別在不調用和調用GPU加速的工況下開展仿真模擬，兩種工況下計算5000個時間步的機器運行時間如圖5所示，在SPH方法中采用GPU加速技術的計算效率相對無加速時提高了61.8倍。由此可見，在SPH方法不但能夠對溢洪道泄流問題進行精確地模擬，通過結合GPU加速技術，還能夠顯著提高該方法的數值仿真效率，有助于SPH方法在工程問題中的推廣應用。

5臺階式溢洪道泄流問題數值模擬

截至目前，水利工作者對溢洪道的泄流問題開展了較多的數值研究，但受到計算規模和效率的限制。采用SPH方法對該問題研究時，通常在蓄水池處分配了較小的計算域，難以保持泄洪流動狀態的穩定性。本文采用的GPU加速技術有助于提高SPH方法的計算規模和效率，故而可采用較寬的計算域對該泄流問題進行數值模擬。

5.1計算模型及工況

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

5.2數值結果及分析

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

圖8為臺階式溢洪道和光滑溢洪道的自由水面演化過程，包括開始泄流的瞬間、水流與溢洪道相互作用的初始階段、水流在溢洪道上的穩定泄流階段等。由圖8可見，泄流開始時刻（t=0.5s），兩種溢洪道的流態高度相似，上游水平堰頂的水深均為0.11m，水頭越過水平堰頂的右端向溢洪道泄流。在t=1.0s時刻，光滑溢洪道上的水流表面較為光順;在臺階式溢洪道上，水流因與臺階的碰撞而產生明顯的翻卷、融合現象，自由面的形狀粗糙。在t=2.0s和t=3.6s時，光滑溢洪道上水流的速度和自由面形態基本相似，泄流過程中溢洪道的底部水流速度相對，上部明顯增高，達到3.5m/s。該階段的泄流過程得到充分發展，流動狀態穩定。在t=2.0s時刻，更多的水流粒子沖擊到臺階上，并因局部的水流翻卷而形成大量的粒子聚集。在t=3.6s時刻，臺階式溢洪道上水流狀態達到穩定狀態，并形成水流截面上的流速分層現象，即臺階近壁面上流速較低，水流表面處速度較高，但溢洪道上水流最大速度為3.29m/s，低于光滑溢洪道上的流動速度。

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

圖9通過矢量圖的形式展示了臺階式溢洪道上流體粒子的運動速度?？梢姡诟骷壟_階近壁面處均出現了渦旋。由圖10的局部放大矢量圖可見，臺階的整個夾角區域均由流體渦旋占據，夾角區域外側與光滑溢洪道上流動狀態相似，流體質點以滑移流的方式向下游運動，但在臺階，上滑移流和渦旋之間存在動量的轉換，進而導致水流動能的沿程耗散。

圖11為兩種溢洪道泄流的沿程水深對比。由圖可見，兩種溢洪道上游泄流段水深高度基本一致，均約為0.081m，表明其來流狀態相同。隨著向下游的泄流過程，溢洪道上的水深逐步降低。其中，光滑溢洪道上的水深沿程降低更為迅速，在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）流速才接近穩定，這表明在溢洪道的下半段臺階內的渦旋對水深的影響范圍更大。

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

公式

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

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

6結語

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

參考文獻：

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

[2]任雨，王承恩，劉斌.階梯式溢洪道水平面上時均壓強試驗研究[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].應用力學學報，2005，22（1）：30-35.

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

[8]趙相航，解宏偉，郭馨.臺階式溢洪道水力特性數值研究[J].青海大學學報，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方法在溢洪道水流模擬中的應用[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