退水渠道水力特性模拟的SPH方法及工程应用

高仁祖，马晓泉，李向东，顾声龙

(1.青海省水利水电勘测设计研究院有限公司，青海 西宁 810001；2.青海大学水利电力学院，青海 西宁 810016；3.黄河上游生态保护与高质量发展实验室，青海 西宁 810016)

0 引言

实现生态文明、绿色发展是我国西北社会经济发展的重要目标。随着兰-西城市群发展规划理念的提出，加快推进“一优两高”战略以及“大西宁”的规划构建，“引黄济宁”工程作为解决区域发展不平衡不充分问题、实施扶贫开发和改善民生具有重要意义[1]。“引黄济宁”工程分为引水工程和供水工程两部分，工程规划从黄河干流上游龙羊峡水库引水，经隧洞穿越拉脊山自流输入到湟水南岸，向湟水干流及南岸、城乡生活、工业园区、农业生产、生态环境建设保护供水，以提高水安全保障能力。

在供水工程中，为保证渠道及重要建筑物的运行安全和渠道分段检修，满足渠道安全运行等需要，结合实际地形沟道设置相应的退水建筑物。由于退水渠的特殊性，退水流量较大，坡度较陡，为防止水流对退水出口处地面造成冲刷，退水渠末端设置消力池进行消能。因此，退水渠的水力特性研究对工程设计和运行管理具有十分重要的意义。通常规模以上重要退水渠的水力特性研究一般以物理模型试验[2]和数值模拟[3]为主，由于资金有限以及场地的限制，所以本次研究使用数值模拟方法。本文基于传统的拉格朗日流体力学光滑质点流体动力学(SPH)方法，对退水渠的水力特性进行数值模拟。

SPH(Smoothed Particles Hydrodynamics)方法首次被Gingold和Monaghan提出[4]，近30多年的发展，在流体动力学领域应用广泛，比如，溃坝流[5]、自由表面流动[6]、多相流[7]等，相应成熟的SPH程序软件有DualSPHysics[8]、SPHinXsys[9]、GPUSPH[10]等。DualSPHysics开源程序包已得到国内外许多学者的验证[11-12]以及应用[13-14]，2019年张云云[15]等数值模拟陡坡水跃，结果与文献中的试验结果基本吻合，验证了SPH出入流方法对陡坡水跃的数值模拟可行性。然而，国内外对退水渠开展的数值模拟以及物理模型试验甚少，相关成果不多，尤其通过SPH方法数值模拟退水渠并未给出相应合理的研究成果。

因此，本文使用DualSPHysics 5.0开源程序代码对退水渠进行三维数值模拟，在4种流量的工况下，分析退水渠水流的水面线变化、流速变化以及渠道末端消力池的消能效果，并得出一些趋势性的结论，更好地指导工程设计和运行管理。

1 模拟方法

1.1 SPH方法的基本方程

SPH方法认为，如果场变量在计算区域中是连续且光滑的，则粒子上的场变量可以由周围粒子的场变量值核估计而来。SPH方程的建立有2个重要的步骤。第一步是核计算，第二步是粒子估计。核估计完成场变量或场变量梯度的插值，而粒子估计则完成对核估计积分表达式中的粒子离散。

1.2 核估计与粒子估计

对于一个连续光滑函数f(x)，在定义域Ω上一点的函数值可以表达为：

(1)

其中，

(2)

式中，x—空间位置矢量；δ(x-x′)—狄拉克δ函数。

若δ函数用函数W(x-x′，h)替代[16]，则f(x)可近似为：

(3)

式中，W(x-x′，h)—SPH中的核函数，它的值取决于两点之间的距离|x-x′|和光滑长度h，它与光滑因子k一同决定了光滑函数影响域的大小。

对于(3)式，进行粒子离散化为：

(4)

式中，mj—配对粒子j的质量；ρj—配对粒子j的密度。

对粒子i处场函数的粒子估计最终为：

(5)

其中，

Wij=W(xi-xj，h)

(6)

粒子i处场函数空间导数的粒子估计最终[17]为：

(7)

其中，

(8)

rij=|xi-xj|

(9)

式中，rij—粒子间距；xi—i粒子的矢量位置；xj—j粒子的矢量位置。

1.3 控制方程

在流体动力学中，流体运动的控制方程(Navier-Stokes方程)表示为：

(10)

(11)

式中，ρ、v—流体的密度、速度；P、fs—流体的静水压力、表面张力；μ—流体的动力黏性系数；t—时间。

在不可压缩流动问题中，利用弱可压缩[18]状态方程表示如下：

(12)

(13)

式中，B—参考压强；γ—常数，在水中γ=7；ρ0—水的密度1000kg/m3；c0—声速，应在流体最大速度的10倍以上，来保证模拟流场的不可压缩性。

SPH方法通过核函数和粒子近似，将拉格朗日形式的N-S方程式(10)和式(11)离散为：

(14)

(15)

式中，Wij—核函数；Πij—人工黏性项[19]。

2 模型布置与方法验证

2.1 几何模型布置

根据“引黄济宁”工程中的众多退水渠原型，如图1所示建立了1∶1的退水渠数值几何模型，其中主干渠是“引黄济宁”工程中的供水渠道，分水池的建立是为了调取部分主干渠的水量对周边县级区域起到灌溉作用，本文是研究退水渠的水力特性，因此不再考虑分水池调水灌溉这一作用。图1中节制闸与渠道入口的距离为10.6m，退水渠的渠道中心与渠道入口的垂直距离为5.95m，圆形分水池的半径为3.8m，坡度为1∶6，消力池底坎高度为1m，其它渠段尺寸见表1。在数值模拟中，粒子间距dx为0.1m，时间步长为1.5×10-4s，模拟总时长80s，使用出入流开放边界条件[18]来稳定流量，入流面积为2m×2.1m(宽×高)，4种退水流量分别为5.7m3/s(最小流量)、5.99m3/s(设计流量)、7.29m3/s(加大流量)和8.4m3/s(限制流量)。节制闸的开闭状态在数值模拟时为完全闭合。

表1 几何模型各部分的尺寸大小

图1 退水渠三维几何模型布置图及坐标原点分布

2.2 本文的SPH方法在陡坡水跃中的验证

为了验证DualSPHysics方法在数值模拟陡坡水力特性方面的准确性及可行性，选取文献中陡坡水跃研究的方案，文献中葛旭峰[20]研究不同坡度不同流量的陡坡水跃特性，给出了试验值与流体软件Fluent模拟值一致的结论。如图2所示，参考了文献[20]中的几何模型尺寸建立了数值几何模型，选取文献[20]中坡度为30°、陡坡长度为0.6m、流量Q分别为42.17、25.18L/s的工况，且临界水深hk分别为0.1043、0.074m。在模拟过程中选择出入流开放边界条件[15]，粒子间距dx=0.01m，时间步长为1.5×10-4s，模拟时长为50s。图2中x表示从坡底的折坡点开始x方向上的距离，即离消力池起点的绝对距离。

图2 陡坡三维几何模型布置图

如图3所示，给出了不同流量的试验值[20]与DualSPHysics方法的模拟值的水面线变化曲线图，x/hk表示离消力池起点的相对位置。从图3中分析得出，模拟值变化与试验值变化趋势一致，然而此处为水流从急变流向缓流变化的转折处，水面线起伏较大，流态处于复杂且随机的湍流状态，因此，此处水面线的模拟值与试验值之间的误差较大，其误差大于其余段的水面线误差。当流量Q=25.18L/s时，模拟值与试验值[20]的平均误差为2.47%；流量Q=42.17L/s时，模拟值与试验值[20]的平均误差为-8.7%。在水跃旋滚剧烈处，水面波动明显，数值模拟情况较差，但水面线总体趋势一致，通过水面线的平均误差分析，得出DualSPHysics方法在数值模拟陡坡水力特性方面，具有一定的准确性以及可行性。

图3 水面线对比曲线图

3 计算结果及分析

图4 四种流量下，沿程水面线变化曲线图

如图5所示，给出了设计流量5.99m3/s和加大流量7.29m3/s的流速分布云图，对退水渠的流场进一步分析，从图5中分析得出：整个渠道中的水流流态稳定，无严重的水流冲击、飞溅等对工程不利的现象。整个退水渠流速最大处分布在陡坡段末，然而，流速越大的位置，受到水流的冲刷越大，因此，在具体工程中应该对冲刷严重的部位加强巩固。通过整个渠道流速色值来分析，陡坡色值表较大，而水流进入消力池后，色值变低，从而更加反映了消力池削减水流流速的效果。如图6所示，给出了4种流量下沿退水渠的流速变化曲线，从图6中分析得出：流速整体沿陡坡增大，图6中色值明显增加；在陡坡段末流速受到消力池旋滚区域的影响，水流流态处于湍流，水流之间存在着许多个随机的大大小小的漩涡结构，漩涡结构之间的剪切作用使得流速趋势变得急剧下降，最后沿消力池段变得平缓。最大流速6.53m/s发生在陡坡段末x=35m的位置。

图5 两种流量下，流速分布三维云图

图6 四种流量下，流速分布曲线图

最后，如图7所示，对沿退水渠变化的消能率进行分析，从图7中分析得出：消能率沿退水渠增大，在消力池段，整体消能率的增长变得缓慢。在消力池段，消能率随着流量的增大，消能效果反而下降。在x=35m至x=45m之间，由于水流流态复杂，处于湍流旋滚区域，蓝点线和红点线消能率曲线呈先减小后增大趋势。然而，总体消能率的趋势沿程呈增长趋势，总体消能率在55%左右，消能效果显著。其中消能率的计算公式[22]如下：

图7 四种流量下，消能率分布曲线图

(16)

式中，Ei—断面单位重量液体所具有的总能量；i=1，2，分别表示跃前、后断面；zib—渠道底面面标高；vi—断面平均流速；hi—断面平均水深。

4 讨论与结论

退水渠是渠系建筑物不可缺少的部分，在“引黄济宁”的供水工程各段分布众多，退水渠承担着排泄灌溉渠道内剩余水量或入渠洪水的责任，不同流量的退水渠会产生不同的水力特性，这些变量最终将会影响整个退水渠排水的功能，影响着下游原始河道防冲消能，因此，对退水渠的数值模拟其意义重大。DualSPHysics 5.0开源程序代码是基于SPH粒子法比较成熟的流体软件，在国内外应用广泛[8，11-12]，在自由表面流方面的数值模拟技术趋于成熟，其次三维的数值模拟相比于以往的二维流体仿真[11-12]，在模拟整个流场方面会更加真实以及数值结果精度会提高很多。本文基于DualSPHysics方法对退水渠进行三维数值模拟，通过四种不同流量的工况，分析研究退水渠水流沿程水深、水流流速以及消能率的变化规律，得出以下结论：

(1)当退水渠道体型一定，入流流量越大，渠道沿程水深在消力池中也随之增大。根据GB 50288—2018《灌溉与排水工程设计标准》[21]对分水池段的水面线超高验算，得出加大流量为7.29m3/s时退水渠的边墙高度设计合理。

(2)对退水渠中水流流速分析，在陡坡段末x=35m处的渠道水流流速最大，此位置受水流冲刷影响较大，需在具体工程中对冲刷严重的部位加强巩固。通过陡坡段以及消力池段的流速云图发现，色值变化明显，极大的反应了消力池的效能效果，同时消力池中的流态为湍流，表现出流态稳定且无明显不利工程运行的问题。

(3)退水渠道消能率随着入流流量的减小而增大，即消力池在小流量的工况下消能效果更好，总体消能率达到了55%左右，消能效果明显。

(4)本文只是对流量变化下的退水渠的水力特性进行了数值分析，得出了一些趋势性的结论，但是对于不同尺寸退水渠的数值分析没有系统的论证，希望在今后的时间里做出合理的成果。