多因素影响的梯形明渠三维流动规律

白 伦, 吴永婷, 孙洪广*

(1.河海大学力学与材料学院，水文水资源与水利工程科学国家重点实验室，南京 210098；2.山东鲁轻安全评价技术有限公司，济南 250100)

梯形明渠是水工建筑中常见的引水渠道，广泛应用于农田水利灌溉和城市给排水工程。一方面，明渠流量的测量是水资源配置优化和科学管理的基础。要测量明渠断面流量，无论是采用经验公式、半经验半理论公式，还是理论方法，均需先掌握明渠流动的特性。另一方面，明渠的施工、监测和维护在实际工程中也是重要的课题，包括最少工程量、明渠使用寿命、边坡稳定及渗流等问题，这些都需要进一步地了解多因素情形下明渠的流动规律。

目前关于明渠流动的研究，主要集中在矩形明渠的流动。矩形明渠实验发现，预期在水面的最大流速点发生在水面以下，它是一组较强的自由面涡与较弱的底层涡相互作用的结果[1]。进一步的实验结果表明：二次元流引起的主流区流速波动分别为10%[2]、5%[3]和2%[4]，同时也导致了横截面上切应力的波动。Rodriguez等[5]通过实验证明，在宽浅明渠中，尽管二次元流会影响流动的横向动量传输，从而进一步影响流速分布，但二维流动的流向流速分布依然遵守对数律公式。Bonakdari等[6]将矩形渠道划分出一个梯形中心区和两个三角形壁面区来考虑壁面影响。Mahananda等[7]对窄明渠流动的研究发现，对于发展中的湍流，宽深比的变化会影响整个水深的流动特性，而对于充分发展的湍流，流动特性的影响只表现在外区。应用流速公式预测流量时，Rickenmann等[8]通过实验数据评估了6种广泛应用的流阻方程和结合了两种幂律流阻方程的Ferguson法的预测效果，后者表现更佳。Johnson等[9]通过使用由纵向流速和横向流速相关关系得到的积分尺度这一局部尺度，再结合对数律流速公式，来远程监测河道流量。

而对于梯形明渠，其坡度的存在使得湍流结构较矩形明渠更加复杂，在自由液面和侧壁之间产生了另一组涡[10]。张长高[11]用镜像法将梯形断面拓展，应用Fourier双重级数得到了流速分布的解析解，但公式较为复杂。胡云进等[12]对12种不同底坡的梯形明渠进行了数值仿真，认为梯形明渠内区流速符合对数律分布，而外区流速符合幂律分布。在梯形明渠均匀流水力计算中，王正中等[13]、赵延风等[14]、黄朝煊[15]的梯形明渠正常水深计算公式都是通过迭代提高精度。对于广义梯形中的自由溢流，Vatankhah[16]认为当渠道的坡度为负、零或平缓时，上游断面的将是临界流量，端截面深度仅由引水渠的形状及其临界深度决定，并提出了一个基于亚临界流动出口水深的直接流量公式。Abrari等[17-19]结合能量方程和连续性方程，分别利用垂向上3个不同点的速度[17]，耦合流线曲率和边缘压力分布[18]，利用端界面几何中心的速度[19]，来计算梯形明渠临界流动条件下的流量。从香农熵的角度，Sheikh等[20]提出了一种梯形明渠边界剪应力的新计算方程。梯形明渠流动内部规律的研究仍相对缺乏。

从工程设计方面考虑，王开民[21]推导出了最佳断面边坡系数计算公式，可以使挖砌方量最少。González等[22]结合曼宁公式和牛顿迭代法来推求最优梯形渠道。但实际施工中，还要结合边坡的材料、受力和渗流[23-24]，柴春岭等[25]利用变水位静水法，构建梯形渠道渗漏强度函数。Roushangar等[26]综合考虑了衬砌成本、土方成本、渗漏和蒸发损失等方面的影响，来进行渠道结构的优化设计。

由于现有湍流理论的不完善，未能充分解释天然明渠流动的性质，实验和仿真都是流体力学中不可缺少的研究手段，被应用于水力计算公式的建立。目前流行计算流体动力学(computational fluid dynamics，CFD)软件Fluent有着丰富的物理模型，能够模拟和分析复杂几何体内流体的流动，被广泛应用于灾害预测、科研和工程施工[27-31]。在上述研究成果的基础上，结合数值模拟，进一步研究多因素影响下梯形明渠的流速分布规律。先对一特定明渠流动进行分析，再综合考虑不同坡度、不同流速的情形，总结出流动规律，为进一步的梯形渠道设计、施工提供参考资料。

1 理论框架

采用计算流体动力学软件Fluent对梯形渠道进行数值模拟。控制方程采用的是Realizablek-ε模型，引入了湍动能(k)和耗散率(ε)方程，具有可靠、收敛性好、内存需求低的优点，适用于模拟流动分离和复杂二次流[32]。相比于其k-ε模型，Realizablek-ε模型湍动黏度的计算引入了旋转和曲率，ε方程倒数第二项不再具有奇异性。

连续性方程为

(1)

动量方程为

(2)

湍动能方程为

(3)

耗散率方程为

(4)

2 几何建模及边界条件

参照某灌区实地测量的参数[33]，设计模型长12 m，底部宽2 m，水深0.8 m。本文使用六面体结构化网格，垂向划分60层，横向划分80层，网格数约为393 000个。断面网格划分如图1所示。

图1 梯形断面网格划分Fig.1 Grid division of trapezoidal section

由于研究的是稳态流动，对自由水面采用刚盖假定处理，即假设流动稳定时的自由水面是恒定的。上游进口为速度进口，下游出口为自由出流，自由水面为对称边界，壁面及床面边界为无滑移壁面，用壁面函数近似处理。速度进口条件选择湍流强度I和水力直径D。记梯形明渠底宽为b/m，边坡坡度为1:n，全水深为h/m，则水力直径表达式为

(5)

湍流强度表达式为

(6)

式(6)中：U为平均流速；ν为运动黏滞系数。模型的具体计算参数见表1，使用SIMPLE算法进行求解，当计算各变量残差小于10-5或随迭代次数不变时，计算收敛。

表1 速度进口参数设定值Table 1 Set values of velocity inlet

3 计算结果及分析

对本次模拟中边坡坡度为1的梯形明渠进行分析，取对称面x=0，流速的沿程分布，i点流速以w=im/s为原点，如图2所示(横向流速为u，垂向流速为v，纵向流速为w，分别对应x、y、z轴)。

图2 垂向速度的沿程分布Fig.2 Vertical velocity distributions along the flow path

流动稳定后，梯形明渠对称面x=0上的速度分布图形沿程变化不明显。笔者使用多种函数拟合并对比误差，发现底部的流速符合对数分布律，与胡云进等[12]结果类似，但范围为y<0.2 m，而外区(0.2 m

图3 y=0.7 m处横向各点速度的沿程变化Fig.3 The variations of velocity at several transverse points of the height y=0.7 m along the flow path

图4 y=0.7 m处横向各点静压的沿程变化Fig.4 The variations of static pressure at several transverse points of the height y=0.7 m along the flow path

考虑梯形明渠边壁上不同高度的壁面切应力(τ)的沿程变化，结果如图6所示。

从图6可知，壁面切应力的值在入口处急剧增大，在z=1 m附近有极大值，随着流程的增大，总体上各点的值先减小后增大，但变化幅度在0.5 Pa以内。在出口处，除y=0.8 m外的点的切应力大致相等。壁面切应力的最大值出现在边壁y=0.8 m、z=1 m附近，其次是边壁y=0 m、z=1 m附近。

接着考虑梯形明渠边壁上不同高度总压(P)的分布及出口与入口的差值，结果如图7和图8所示。

图5 y=0.7 m处横向各点总压的沿程变化Fig.5 The variations of total pressure at several transverse point of the height y=0.7 m along the flow path

图6 边壁上不同高度的壁面切应力的沿程变化Fig.6 The variations of wall shear stress of different heights on the walls along the flow path

图7 边壁上不同高度的总压沿程分布Fig.7 The distributions of total pressure of different heights on the walls along the flow path

图8 边壁上不同高度的入口和出口的总压差Fig.8 Total pressure differences at inlet and outlet at different heights on the walls

从图7可以看到，梯形明渠边壁上高度为0 m的总压最大，因此最容易发生渗漏。根据图8可知，总压差在高度为y=0.1 m和y=0.4 m处分别有一个极大值和极小值，而且在y=0.4 m处取最小值，y=0.8 m处取最大值，说明在边壁上0.1 m和0.8 m附近，机械能损失更大，也说明在0.1 m和0.8 m附近流动对边壁的侵蚀更加严重。

接着，对比分析不同边坡坡度(1:n)的梯形明渠中的流动，取流速为1 m/s，壁面切应力(τ)和总压(P)的结果如图9和图10所示，可以看到，当1/n≤1.5时，边壁和床面的壁面切应力的最大值τwm和τbm随坡度变化平缓；当1/n>1.75时，边壁和床面的壁面切应力的最大值τwm和τbm随坡度增大，且趋于相等。而边壁和床面的壁面切应力平均值τwa和τba随坡度的变化均不明显。

图9 边壁和床面的壁面切应力随坡度的变化Fig.9 The variation trends of wall shear stress on the walls and bed with slopes

对于边壁和床面总压P，当1/n>1.5时，边壁和床面的总压的最大值Pwm和Pbm随坡度增大，且趋于相等。边壁和床面的总压的平均值Pwa和Pba随坡度的变化均不明显。

接着，对比分析梯形明渠中不同平均流速(U)的情形，取坡度为1，壁面切应力和总压(P)的结果如图11和图12所示。

图10 边壁和床面的壁面总压随坡度的变化Fig.10 The variation trends of total pressure on the walls and bed with slopes

图12 边壁和床面的总压随平均流速的变化Fig.12 The variation trends of total pressure on the walls and bed with mean velocity

可以看到，边壁和床面的壁面切应力随速度单调递增，当U>3.5 m/s时，边壁的最大值τwm趋近于床面的最大值τbm。二次项拟合得到公式为

τwm=1.56U2-0.13U+2.34

(7)

τbm=0.54U2+4.28U+0.68

(8)

τwa=0.99U2+0.80U+0.24

(9)

τba=0.57U2+3.09U+0.49

(10)

同时，边壁和床面的总压P也随速度单调递增，当U>1.5 m/s时，边壁的最大值趋近于壁面的最大值。二次项拟合得到公式为

Pwm=533.6U2-196.4U+8 073.7

(11)

τbm=581.3U2-534.5U+8 657.0

(12)

Pwa=316.6U2-8.7U+3 891.6

(13)

Pba=315.9U2-2.7U+7 988.1

(14)

4 结论

基于刚盖假定，使用Fluent对梯形明渠中的流动进行了三维数值模拟。先对边坡坡度比为1、平均流速为1 m/s的流动进行分析，接着对比研究了不同坡度1:n、不同平均流速U情形下边壁和床面切应力τ和总压P，得出如下结论。

(1)对称面上速度分布的沿程变化不明显，底区(y<0.2 m)流速满足对数律分布，外区流速(0.2 m

(2)同一水平面上的各点，速度值和总压随流程的减小趋势，离壁面越近表现得越明显，而其静压的减小趋势则大致保持相同。

(3)边壁上不同高度的各点，壁面切应力在1 m处出现极值，随后减小后增大，变化幅度在0.5 Pa以内，壁面切应力最大值出现在边壁y=0.8 m，z=1 m附近，同时出口和入口的总压差也在y=0.8 m处取最大值，说明此处流动对边壁的侵蚀更加严重。

(4)平均流速为1 m/s的情形，当坡度大于1.75时，边壁和床面的壁面切应力的最大值随坡度增大，且趋于相等；当坡度大于1.5时，边壁和床面的总压的最大值随坡度增大，且趋于相等。边壁和床面的壁面切应力、总压的平均值随坡度变化不明显。

(5)坡度为1的情形，边壁和床面的壁面切应力、总压随速度单调递增，当平均速度大于3.5 m/s时，边壁的壁面切应力最大值趋近于床面的最大值；当平均速度大于1.5 m/s时，边壁的总压最大值趋近于壁面的最大值。用二次项拟合的公式，较符合模拟流动的壁面切应力、总压随速度的变化。