光滑粒子流体动力学方法在拱坝中孔泄流冲击水垫塘模拟中的应用

2022-04-12张建伟主攀陈海舟杜宇

华北水利水电大学学报(自然科学版) 2022年1期

张建伟，主攀，陈海舟，杜宇

(1.华北水利水电大学水利学院，河南郑州 450046； 2.上海道盾科技股份有限公司，上海 200131)

对拱坝中孔泄流过程进行分析，对优化泄流方式、减少坝体空蚀和磨蚀的产生、削弱河床和岸坡冲刷具有重要意义。拱坝中孔泄流时具有水头高、流速大、下泄功率大的特点[1]，物理模型试验获取准确的观测资料难度较大、投资大、周期长，每次物理模型试验只能针对较少的组合工况进行模拟，具有局限性[2]。数值模拟具有建模方便、仿真时间短、成本低的特点，适用于大坝泄流模拟，可作为物理模型试验研究的补充手段。2009年,孙建等[3]在刚盖假定的前提下，采用紊流数学模型模拟了溪洛渡高拱坝水垫塘内三维流场变化过程；2019年,周振兴等[4]采用UG-FLUENT软件对二河新泄洪闸泄流时的下游流场进行了模拟。由于泄流过程自由液面会发生大变形，传统的基于网格的数值方法如有限差分法(FDM)和有限体积法(FVM)只能对下游水垫塘进行模拟，很难对泄流过程进行仿真，寻找新的数值模拟方法已成为目前泄流仿真研究的热点。

光滑粒子流体动力学(Smoothed Particle Hydrodynamics，SPH)方法是一种基于拉格朗日形式的纯粒子法，最初用于解决天体物理学中的三维问题，后由MONAGHAN J[5]对其进行修改，用于解决流体计算中的自由表面流问题。SPH方法可以用来处理自由表面大变形问题，能够模拟出流场的整体形态和局部的破碎、飞溅、卷曲等细节。目前,SPH方法已被广泛应用于溃坝、晃荡、液滴撞击和滑坡涌浪等问题研究中。李大鸣等[6]应用SPH方法并使用边壁加载和相对加速度加载两种加载方式模拟二维矩形容器内液体的纵荡，发现在不同水深下边壁加载和相对加速度加载两种方式下的计算结果较一致；陈海舟[7]应用SPH方法模拟了溃坝问题、容器内液体晃荡问题和液滴冲击液面问题，在SPH基本方程的基础上推导出了适用于广义流体力学的SPH控制方程，为SPH方法在流体力学中的广泛应用奠定了基础；张建伟等[8]应用SPH方法模拟了溃坝冲击水垫塘的过程，对冲击过程中的流速和压强进行了分析，提出SPH方法未来可应用于消能研究。

文中应用SPH方法模拟了拱坝中孔泄流冲击水垫塘过程，并对泄流过程中流速、压强等进行分析，同时研究了下游水垫塘对泄流消能的影响。

1 SPH方法原理

1.1 基本原理

SPH方法是一种纯拉格朗日粒子法，用该方法做相关分析时先进行核近似，然后再对粒子近似[9]。通过积分表示法对场函数进行近似，核近似方程为：

(1)

式中：f为场函数；Ω为计算域；x为坐标向量；h为光滑长度；W为核函数。在SPH方法中习惯用角括弧对核近似算子进行标记。

应用粒子近似法对核近似方程式(1)进一步近似，近似后的函数表达式为：

(2)

式中：ρ为密度；m为质量；i、j分别为粒子序号。

1.2 控制方程

采用无黏假定将二维Navier-Stokes方程离散为SPH形式的流体控制方程。

(3)

动量守恒方程为：

(4)

能量守恒方程为：

(5)

式中:P为压强；v为速度；vij=vi-vj；e为能量；N为粒子支持域中的粒子数。

1.3 人工黏度

SPH算法在模拟冲击波问题时会产生非物理振荡，采用MONAGHAN型人工黏度Πij能有效防止粒子相互接近时的非物理穿透。MONAGHAN人工黏度表达式为：

(6)

其中：

全国已建基层水利站29040个，12个省份实行了村级水管员制度。全国80%的县成立了县级农村饮水安全工程管理机构。在全国27个省55个县深入开展农业水价综合改革示范。取消4项行政审批事项，到2015年年底拟将水利部行政审批事项由48项减少到29项。55个全国示范县市、150个省级试点县市开展小型工程管理体制改革试点工作。

设置αΠ分别取1.0、1.2、1.4、1.6、1.8和2.0，通过试算发现，当αΠ取1.8时，既不会造成水流的间断，同时也能很好地模拟水体大变形状态，故αΠ设置为1.8。

1.4 密度近似法

在SPH方法中,光滑长度的变化和粒子的分配与粒子密度变化具有直接的关系。目前,关于粒子密度近似方面应用最广泛的方法是密度求和法。任意粒子i的密度ρi的计算式为：

(7)

式中Wij为粒子j对粒子i影响的光滑函数。

1.5 边界处理

SPH方法属于无网格方法，只需要考虑固壁边界条件的处理。为防止实粒子穿透边壁，采用第一种虚粒子与第二种虚粒子相结合的方法来布置边壁粒子。

第一种虚粒子称为固壁虚粒子，可以对边壁附近的实粒子施加一个强排斥力，该排斥力的表达式为：

(8)

式中：Fij为固壁粒子j对边壁附近实粒子的排斥力；n1、n2分别取12和4；r0为截止半径；D是由具体问题而定的参数，一般取与速度最大值的平方相等的量级值；rij为粒子i与粒子j之间的距离。

第二种虚粒子称为镜像虚粒子，其密度、光滑长度和速度大小与对应的实粒子相同，速度方向相反。边壁粒子布置方案如图1所示。

图1 边壁粒子布置方案示意图

2 数值计算

2.1 模型验证

拱坝中孔泄流仿真难点在于对自由液面大变形进行模拟，为验证SPH方法是否能准确模拟流体大变形问题，选取文献[10]中溃坝的算例进行验证。采用FORTRAN语言编写了SPH方法的计算程序，应用所开发的程序模拟了溃坝过程。模型为长0.584 m、高0.350 m的单宽矩形槽，上游水位设置为0.293 m，通过快速提闸模拟溃坝过程[11]，模型尺寸如图2所示，计算参数见表1。

图2 溃坝水位模型

表1 模拟溃坝过程的参数条件

图3为模拟所得的不同时刻水体流态与文献[10]的对比结果。

图3 不同时刻的水体流态与文献[10]的结果对比

选取了0.2、0.4、0.8 s三个时刻，绘制出物理模型试验结果中的自由液面曲线和等压线图与本文的模拟结果进行对比，如图4所示。

图4 自由液面比较曲线

在图4中物理模型试验结果中的自由表面曲线和数值模拟结果中的自由表面曲线上分别取若干个点，计算点的位置距离差后与物理模型试验值进行对比，得到各点的位置相似度，取平均值后得到自由液面曲线的相似度，结果见表2。

表2 自由液面的相似度

文中的模拟结果和文献[10]的物理模型试验结果十分接近，自由液面曲线综合相似度为91.7%。相对于文献[10]中MPS半隐式运动粒子法模拟结果，SPH方法模拟得到的液面流态更接近物理模型试验的结果。同时，在图4中的等压线图中可以看到，相对压强为零的等压线与自由液面曲线基本重合，说明压强分布规律合理，验证了文中所编写的计算程序的准确性。在0.8 s时，溃坝洪水形成了漩涡，自由液体表面产生了大变形，说明SPH方法可以准确模拟自由液面大变形问题。

2.2 拱坝泄流冲击水垫塘过程数值模拟

为模拟拱坝中孔泄流冲击水垫塘过程，布置了一个长12.2 m、高3.5 m的单宽矩形水槽模型。距水槽左端4 m处设置一挡水建筑物，该挡水建筑物高2.5 m、宽0.2 m。距其底部1.1 m位置处设置一直径为0.1 m的泄水孔用来模拟中孔泄流。下游布置有水垫塘。上游水位设置为2.5 m，下游水位设置为0.4 m。模型尺寸如图5所示，计算参数见表3，初始粒子布置如图6所示。初始粒子密度分布如图7所示。

图5 中孔泄流水位模型

表3 模拟拱坝泄流冲击水垫塘过程的参数条件

图6 粒子的初始布置位置示意图

图7 粒子的密度分布情况

中孔泄流冲击水垫塘过程中的几个关键时刻的压强图(左)和流速图(右)如图8所示。由图8可以看出：上游水体从泄水孔流出后，由于中孔位置处受到的压强较大，使得下泄水流的流速较快；在0.16 s时，上游孔口处粒子受到的最大压强达到14 kPa，孔口的最大流速为3.2 m/s，出射水股具有极大的动能；从孔口流出后，下泄水体在空中形成射流，这一过程中粒子的势能逐步转化为动能；在0.46 s时，可以观察到粒子的压强进一步减小，流速增大，粒子的最大流速达到5.1 m/s；在0.49 s时，上游水体开始冲击下游水垫塘，发生碰撞的一瞬间，由于受到下游水垫塘粒子的作用力，上游水体粒子的压强突然增大到24 kPa，由于动量守恒，下泄水体粒子的速度减小为3.8 m/s；随后在0.76～1.10 s时，由于上游水体的持续冲击，下游水垫塘形成水舌，并逐渐扩大为涌浪，上下游水体粒子混掺在一起向下游边界处运动；在1.32 s时，由于水体的持续冲击，下游水垫塘再次出现水舌并持续翻滚形成漩涡；在1.49 s时可以清楚地观察到漩涡中心处出现空穴，随后上下游水体粒子再次混掺在一起形成涌浪，推高下游水位，说明出射水股动能逐渐转化为水垫塘内水体的势能。

图8 拱坝中孔泄流冲击水垫塘过程中水流的压强图(左)和流速图(右)

为了更清晰地展示漩涡部分的流场分布，绘制出漩涡形成过程中几个关键时刻的流速矢量图(左)和局部放大图(右)，如图9所示。图9中：在1.32 s时，由于上游水体的流速较大，受到冲击的下游水垫塘形成水舌；随后下游水垫塘的水舌进一步发展，上下游水体混掺后翻滚形成漩涡，在1.49 s时，在漩涡的中心处出现明显的空穴，由于空穴处水体粒子相互挤压，空穴处流速较大；随后空穴部分逐渐被下泄水体填充，空穴处流速减小，实现消能，在1.73 s时，漩涡演变为涌浪，在这一过程中水股动能逐渐转化为内能，具有一定的消能效果。