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