山地灌溉管道水力特性的数值模拟
刘家宏1, 周晋军1, 2, 王浩1, 3, 吕宏兴2    
1. 中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室, 北京 100038;
2. 西北农林科技大学 水利与建筑工程学院, 杨凌 712100;
3. 水利部 水资源与水生态工程技术研究中心, 北京 100044
摘要:有压管道的水力特性是管道设计工作的依据和前提。为研究山地灌溉输水管道水力特性,该文建立山地灌溉输水管道水动力学数学模型,运用计算流体动力学数值模拟方法,对有压灌溉管道中凸起管段的水力特性进行模拟。该文以输入流量为控制变量,计算管道水流能量损失,分析了Re与管道凸起段135°弯头局部阻力系数的关系;研究输水管道凸起段水压和流速分布特性。结果表明:当Re<5.0×104时弯头局部阻力系数随着Re增大而迅速减小,当Re>7.0×104时弯头局部阻力系数趋于稳定,当2.3×104 <Re <7.0×104时,弯头Ⅰ、弯头Ⅱ、弯头Ⅲ、弯头Ⅳ局部阻力系数变幅范围为4.12~0.37;弯头处静压分布均表现出外侧压力大,内侧压力小;速度值的大小均表现为弯头外侧速度小,内侧速度大。
关键词山地灌溉管道    局部阻力系数    Reynolds数    计算流体动力学    
Numerical simulation of the hydraulic characteristics of hilly irrigation systems
LIU Jiahong1, ZHOU Jinjun1, 2, WANG Hao1, 3, LV Hongxing2    
1. State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China;
2. College of Water Resources and Architectural Engineering, Northwest Agriculture and Forestry University, Yangling 712100, China;
3. Research Center for Water Resources and Hydro-ecological Engineering, Ministry of Water Resources, Beijing 100044, China
Abstract: The hydraulic characteristics of pressurized pipelines are important for irrigation system designs. The hydrodynamics irrigation pipes in a hilly terrain were analyzed here using a computational fluid dynamics model. The energy losses in the pipeline were calculated with the input discharge as the control variable. The relationship between the Reynolds number and flow resistance coefficient of a convex section (a 135° elbow) was analyzed. The pressure and velocity distributions are presented for the pipeline. The results show that:1) when the Reynolds number is less than 5.0×104, the flow coefficient of the elbow decreases rapidly with increasing Reynolds number; 2) when the Reynolds number is more than 7.0×104, the flow coefficient is nearly constant; 3) when the Reynolds number is in the range of 2.3×104 and 7.0×104, the flow coefficient ranges in 4.12~0.37. The pressure on the backside of the elbow is high, while the inside pressure is low. The velocity distribution is just the opposite with a low velocity near the backside of the elbow and higher velocities near the inside.
Key words: hilly irrigation systems    flow coefficient    Reynolds number    computational fluid dynamics    

随着灌溉技术在山地、 丘陵地区的推广和发展,山地灌溉面积不断扩大。但由于山地区域地势起伏不断,输水管道布置出现凸起或下凹,导致管道水力特性不同于一般平铺管道。国内外有大量学者研究了有压管道的水力特性问题。刘竹溪和刘光临[1]对泵站水锤的防护及计算进行研究。索丽生[2]对锥管水击计算的特征线法进行研究。杨玉思等[3]对有压供水管道水锤防护进行了研究。万五一[4]对长距离输水系统的非恒定流特性进行了研究。文俊等[5]对90°圆形弯管三维紊流进行数值模拟研究。陈江林等[6]对T型三通管水力特性进行数值模拟与试验研究。李文全等[7]对井筒式潜水轴流泵出水管道水力特性进行数值模拟研究。王梦婷等[8]对正虹吸管道水力特进行试验研究。严继松等[9]对有压管道充水过程水力特性进行数值模拟研究。郑文玲等[10]对异形岔管道水力特性进行数值模拟研究。石喜[11]对灌溉管网非恒定流计算及应用进行了研究。周晋军等[12]对山地灌溉管道含气囊运动水力特性进行了试验研究。Strowger等[13] 和Wood[14]对水电站瞬变过程中调压室漏空压力钢管的气液两相流问题进行研究,同时结合工程实际对水电站系统非恒定流的规律和计算进行了研究。Jayaraj等[15]和 Afshar 等[16]对管道瞬变流、 水锤等问题进行模拟研究。Estrada[17]研究加压管网水力解算器在灌溉系统中的应用。

这些研究主要是集中在有压管道水锤问题和城市供水有压管道的水力特性研究,关于山地灌溉输水管道水力特性的研究较少,特别是在不同流速下管道弯头局部阻力系数的变化问题方面少有报道。为保障山地灌溉输水管道合理布设和优化设计工作,需要对管道输水水头损失进行计算,因此有必要对山地灌溉输水管道的水力特性进行研究。本文通过计算流体动力学数值模拟方法,以流量为控制变量,运用FLUENT软件模拟,研究135°弯头局部阻力系数随Re的变化规律,分析山地灌溉输水管道凸起段水压分布规律和水流速度矢量分布规律。

1 模型建立与网格划分

采用FLUENT前处理软件Gambit建立三维管道模型。坐标原点为管道进口断面圆心处,x正向为管道中水流流动方向,并进行网格划分。图 1为山地灌溉管道凸起部位模型计算域网格划分,模型中管道内径d=65 mm,模拟管道材质是有机玻璃管,其糙率系数n=0.007,模型中有4个135°弯头,弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ位置如图 1所示。弯头Ⅰ和弯头Ⅱ距离为1.35 m,弯头Ⅱ和弯头Ⅲ距离为1.5 m,模型管道爬坡角为45°,模型管道上下坡对称。假设模型是封闭的,管道中不含气囊。模型弯头是由两边管道和管道同心球体嵌套而成,弯头外侧圆滑,内侧稍显棱角,假设弯头不会影响管道光滑度。弯头部位网格划分见图 2

图 1 管段凸起部位管段网格划分图  
图 2 弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ网格划分图  
2 模型求解方法

方程采用有限体积法(FVM)离散,应用QUICK格式,速度场与压力的耦合计算采用SIMPLEC算法[18]。模拟过程边界条件设置为速度进口,即模型输入为速度。

1) 连续性方程:

$ \frac{{\partial \left( {\rho {u_i}} \right)}}{{\partial {x_i}}} = 0. $ (1)

2) 动量方程:

$ \begin{array}{*{20}{c}} {\frac{\partial }{{\partial {x_i}}}\left( {\rho {u_i}{u_j}} \right) = - \frac{{\partial P}}{{\partial {x_i}}} + }\\ {\frac{\partial }{{\partial {x_i}}}\left[{\mu \left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) - \frac{2}{3}\mu {\delta _{ij}}\frac{{\partial {u_k}}}{{\partial {x_k}}}} \right] + \frac{\partial }{{\partial {x_j}}}\left( { - \rho \overline {{u_i}'{u_j}'} } \right).} \end{array} $ (2)
其中:ρ是液体密度,kg/m3uiuj分别是流速分量,m/s; P是流体的压力,N/m2; μ流体动力黏度,N·s/m2δij是“Kronecker delta”符号。当ij时,δij=1; 当ij时,δij=0 。

在山地灌溉输水管道弯头及凸起部位水流存在漩涡流动,因此选用RNG(renormalization group) k-ε湍流模型如式(3)和式(4)所示,即

$ \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho k{u_i}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_i}}}\left[{{\alpha _k}{\mu _{{\rm{eff}}}}\frac{{\partial k}}{{\partial {x_j}}}} \right] + {G_k} + \rho \varepsilon , $ (3)
$ \begin{array}{*{20}{c}} {\frac{{\partial \left( {\rho \varepsilon } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho \varepsilon {u_i}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left[{{\alpha _\varepsilon }{\mu _{{\rm{eff}}}}\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right] + }\\ {\frac{{C_{1\varepsilon }^ * \varepsilon }}{k}{G_k} + {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k}.} \end{array} $ (4)
其中:μeff是修正的湍动黏度;Gk是平均速度梯度引起的湍动能k的产生项; C1ε*是修正系数。 各项物理量的计算如式(5)所示,即

$ \left. {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\mu _{{\rm{eff}}}} = \mu + {\mu _t},{\mu _t} = \rho {C_\mu }\frac{{{k^2}}}{\varepsilon },{C_\mu } = 0.0845;}\\ {{G_k} = {\mu _t}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)\frac{{\partial {u_i}}}{{\partial {x_j}}},{a_k} = {a_{\rm{\varepsilon }}} = 1.39;}\\ {C_{1\varepsilon }^ * \varepsilon = {C_{1{\rm{\varepsilon }}}} - \frac{{\eta \left( {1 - \eta /{\eta _0}} \right)}}{{1 + \beta {\eta ^3}}};}\\ {{C_{1{\rm{\varepsilon }}}} = 1.42,{\eta _0} = 4.377;} \end{array}}\\ {\eta = {{\left( {2{E_{ij}} \cdot {E_{ij}}} \right)}^{1/2}}\frac{k}{\varepsilon };}\\ {{E_{ij}} = \frac{1}{2}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right);}\\ {{C_{2{\rm{\varepsilon }}}} = 1.68;\beta = 0.012.} \end{array}} \right\} $ (5)

凸起段的入口断面满足边界条件为:

$ u = {\rm{const}},v = w = 0,\frac{{\partial p}}{{\partial x}} = 0. $
其中:uvw的物理意义依次是xyz方向的动量;p表示大气压强;const表示常数。

管道出口断面满足边界条件如式(6)所示,即

$ \frac{{\partial u}}{{\partial x}} = \frac{{\partial v}}{{\partial x}} = \frac{{\partial w}}{{\partial x}} = 0,p = {\rm{const}}{\rm{.}} $ (6)
3 数值模拟结果分析

将网格模型导入FLUENT软件,在前面所述假设条件下进行数值模拟计算,数值模拟选择RNG k-ε模型对管道流场进行计算。考虑到弯头间距离远大于管道直径,在局部阻力系数计算过程中考虑了沿程水头损失。

根据Bernoulli方程[19],弯头前后断面间能量守恒方程为式(7),即

$ {z_i} + \frac{{{p_i}}}{\gamma } + \frac{{{a_i}v_i^2}}{{2g}} = {z_j} + \frac{{{p_j}}}{\gamma } + \frac{{{a_j}v_j^2}}{{2g}} + {h_{\rm{w}}}. $ (7)
其中: zizj的物理意义是位能; $ \frac{{{p_i}}}{\gamma } $和$ \frac{{{p_j}}}{\gamma } $物理意义是压能;$ \frac{{{a_i}v_i^2}}{{2g}} $和$ \frac{{{a_j}v_j^2}}{{2g}} $的物理意义是动能;hw的物理意义是水流能量损失。 hw由式(8)可得,即

$ {h_{\rm{w}}} = \sum {{h_{\rm{f}}} + \sum {{h_{\rm{j}}}} } . $ (8)
其中: hf是沿程水头损失; hj是局部水头损失。hfhj由式(9)和式(10)可得,即

$ {\left. {\begin{array}{*{20}{c}} {{h_{\rm{f}}} = \lambda \frac{l}{d}\frac{{{v^2}}}{{2g}},}\\ {{h_{\rm{j}}} = \xi \frac{{{v^2}}}{{2g}}.} \end{array}} \right\}} $ (9)

Re计算方法如式(10)所示,即

$ {\mathop{\rm Re}\nolimits} = \frac{{vd}}{v}. $ (10)
其中: Re为Reynolds数; v代表流速; d代表管道内直径; ν代表液体的运动粘滞系数,ν=1.306×10-3 cm2/s

图 3给出了Re与弯头局部阻力系数的对数相关关系图,纵坐标为弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ的对数,横坐标为Re对数。可以看出,弯头Ⅰ 局部阻力系数与Re的对数相关性稍差,其余的较好。

图 3 Re与弯头局部阻力系数的对数相关关系  

图 4给出了管段凸起部位弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ局部阻力系数随Re的变化而变化的数值模拟结果。Re的区间是2.3×104~1.4×105,对应的速度区间为0.46~2.83 m/s。

图 4 弯头Ⅰ、Ⅱ、Ⅲ、Ⅳ局部阻力系数与Re关系  

图 4可以看出,Re<5.0×104时,弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ局部阻力系数都随着Re的增大而迅速减小; 当Re>7.0×104后,局部阻力系数逐渐趋于稳定; 当Re在2.3×104~7.0×104之间变化时,弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ局部阻力系数变幅范围为4.12~0.37。

图 5给出了山地灌溉管道在入口流速为 1.5 m/s (Re=74 655)工况下的凸起部位各管段静压分布图。选取这个工况的原因是局部阻力系数在此趋于稳定。从图 5中可以看出,管段入口处压力最大,出口处压力最小,从入口到出口管道输水静压逐渐减小; 坡前段、 上坡段的静水压力明显大于下坡段和坡后段压力。

图 5 管段凸起部位管段静压分布图  

图 6给出了入口流速为1.5 m/s工况下,管段凸起部位弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ静压分布图。从图 6可以看出,弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ处静压分布均表现出外侧压力最大,内侧压力最小,中间压力居中的特点; 且在每个弯头的内外侧都分别至少有3个压力分区。就弯头外侧(内侧)的压力来看,压力从大到小的弯头依次是为弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ,即与图 5中的压力变化趋势相同。弯头Ⅰ外侧的静压是凸起部位管段中除进水口处压力最大的区域,弯头Ⅳ的内侧静压是凸起部位管段中输水常态下唯一的负压区,弯头Ⅱ和弯头Ⅲ的压力值介于弯头Ⅰ和弯头Ⅳ的压力值之间,处于较为安全的位置。

图 6 弯头Ⅰ、弯头Ⅱ、弯头Ⅲ、弯头Ⅳ静压分布图  

图 7给出了进水口流速为1.5 m/s工况下弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ数值模拟显示的速度矢量图。从图 7可以看出弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ处速度值的大小均表现为弯头外侧速度小,内侧速度大,与图 6中静压大小分布特点正好相反,符合能量守恒定律。同时,在弯头处速度矢量图中表现为速度矢量的方向是从弯头内侧指向外侧,呈现扩散的流态,在弯头外侧均会出现漩涡区,在弯头内侧偏下游处也会有狭长的漩涡区出现。

图 7 弯头Ⅰ、弯头Ⅱ、弯头Ⅲ、弯头Ⅳ速度矢量图  
4 结 论

本文通过计算流体动力学数值模拟方法,运用FLUENT 软件进行模拟,对山地灌溉管道水力特性进行了研究,得到如下主要结论:

1) 通过数值模拟,计算分析了弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ局部阻力系数与Re的关系,得出弯头局部阻力系数随Re增大而减小,最后基本趋于稳定,模拟结果符合管道水动力学理论。

2) 通过数值模拟研究发现: 输水常态下,管段凸起部位管段弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ处静压分布均表现出外侧压力最大,内侧压力最小; 速度值的大小均表现为弯头外侧速度小,内侧速度大。模拟结果符合能量守恒定律。

3) 本次数值模拟结果表明,当Re在 2.3×104~7.0×104 变化时(即管道水流流速在 0.46~1.40 m/s 变化时),弯头Ⅰ、 弯头Ⅱ、 弯头Ⅲ、 弯头Ⅳ的局部阻力系数的变幅范围为4.12~0.37,变幅较大。研究结果表明: 不同设计流速下,管道的局部水头损失系数取值变化大。建议在山地灌溉管道设计中予以关注。

参考文献
[1] 刘竹溪, 刘光临. 泵站水锤的防护措施及其简易计算(上)[J]. 农田水利与小水电, 1984, 59(5):32-37.LIU Zhuxi, LIU Guanglin. Pumping water hammer protection measures and simple calculation[J]. Irrigation and Water Conservancy and Hydropower, 1984, 59(5):32-37. (in Chinese)
[2] 索丽生. 锥管水击计算的特征线法[J]. 水力发电学报, 1997, 58(3):61-68.SUO Lisheng. Method of characteristics for computation of water hammer in conical tubes[J]. Journal of Hydroelectric Engineering, 1997, 58(3):61-68. (in Chinese)
[3] 杨玉思, 徐艳艳, 羡巨智. 长距离高扬程多起伏输水管道水锤防护的研究[J]. 给水排水, 2009, 45(4):108-111.YANG Yusi, XU Yanyan, XIAN Juzhi. Research on water hammer prevention in high-lift, hilly and long distance water transmission pipeline[J]. Water & Wastewater Engineering, 2009, 45(4):108-111. (in Chinese)
[4] 万五一. 长距离输水系统的非恒定流特性研究[D]. 天津:天津大学, 2004.WAN Wuyi. Study on Unsteady Flow in Long-Distance Water Diversion Projects[D]. Tianjin:Tianjin University, 2004. (in Chinese)
[5] 文俊, 刁明军, 李斌华, 等. 90°圆形弯管三维紊流是指模拟[J]. 四川水力发电, 2008, 27(2):111-113.WEN Jun, DIAO Mingjun, LI Binhua, et al. Numerical simulation on three-dimensional turbulent of 90° circular bends[J]. Sichuan Water Power, 2008, 27(2):111-113. (in Chinese)
[6] 陈江林, 吕宏兴, 石喜, 等. T型三通管水力特性的数值模拟与实验研究[J]. 农业工程学报, 2012, 28(5):73-77.CHEN Jianglin, LV Hongxing, SHI Xi, et al. Numerical simulation and experimental study on hydrodynamic characteristics of T-type[J]. Transactions of the Chinese Society of Agricultural Engineering, 2012, 28(5):73-77. (in Chinese)
[7] 李文全, 杨祖欣, 游强强, 等. 井筒式潜水轴流泵出水管道水力特性数值模拟研究[J]. 水电能源科学, 2013, 31(1):150-153.LI Wenquan, YANG Zuxin, YOU Qiangqiang, et al. Numerical simulation of hydraulic characteristics of water outflow pipe of well-type submerged axial flow pump[J]. Water Resources and Power, 2013, 31(1):150-153. (in Chinese)
[8] 王梦婷, 李琳, 谭义海, 等. 正虹吸管道水力特性试验研究, [J]. 水电能源科学, 2014, 32(12):87-91.WANG Mengting, LI Lin, TAN Yihai, et al. Hydraulic model test research of new pressure regulating device for small and medium diversion type hydropower station[J]. Water Resources and Power, 2014, 32(12):87-91. (in Chinese)
[9] 严继松, 廖国玲. 有压管道充水过程水力特性三维数值模拟[J]. 水利水电技术, 2015, 46(3):110-114. YAN Jisong, LIAO Guoling. 3-D numerical simulation on hydraulic characteristics of water filling process of pressure pipeline[J]. Water Resources and Hydropower Engineering, 2015, 46(3):110-114. (in Chinese)
[10] 郑文玲, 张耀哲, 杨石磊, 等. 异形岔管水力特性的数值模拟[J]. 西北农林科技大学学报:自然科学版, 2014, 42(11):183-190.ZHENG Wenling, ZHANG Yaozhe, YANG Shilei, et al. Numerical simulation of hydraulic characteristics in heterotypic bifurcated pipe[J]. Journal of Northwest Agriculture and Forestry University:Natural Science Edition, 2014, 42(11):183-190. (in Chinese)
[11] 石喜. 灌溉管网非恒定流计算及应用研究[D]. 杨凌:西北农林科技大学, 2013.SHI Xi. Research on Caculation and Application of Unsteady Flow in Irrigation Network[D]. Yangling:Northwest Agriculture and Forestry University, 2013. (in Chinese)
[12] 周晋军, 吕宏兴, 朱德兰. 山地灌溉管道含气囊运动的水力特性研究[J]. 人民黄河, 2013, 35(12):101-103.ZHOU Jinjun, LV Hongxing, ZHU Delan. Research on hydraulic characteristics of mountainous irrigation pipe with air movement[J]. Yellow River, 2013, 35(12):101-103. (in Chinese)
[13] Strowger E B, Derr S L. Speed changes of hydraulic turbines for sudden changes of load[J]. Journal of Turbomachinery, 1926, 48:209-262.
[14] Wood F M. Discussion of speed changes of hydraulic turbines for sudden changes of load[J]. Journal of Turbomachinery, 1926, 48:56-68.
[15] Jayaraj K, Ganesan N, Padmanabhan C. A new finite element formulation based on the velocity of flow for water hammer problems[J]. International Journal of Pressure Vessels and Piping, 2005, 82:1-14.
[16] Afshar M H, Rohani M. Water hammer simulation by implicit method of characteristic[J]. International Journal of Pressure Vessels and Piping, 2008, 85:851-859
[17] Estrada C, Gonzalez C, Aliod R, et al. Improved pressurized pipe network hydraulic solver for applications in irrigation systems[J]. Journal of Irrigation and Drainage Engineering, 2009, 135:421-430.
[18] 王福军. 计算流体动力学分析——CFD软件原理与应用[M]. 北京:清华大学出版社, 2004. WANG Fujun. Computational Fluid Dynamics Analysis——Software of CFD Principles and Applications[M]. Beijing:Tsinghua University Press, 2004. (in Chinese)
[19] 吕宏兴, 裴国霞, 杨玲霞. 水力学[M]. 北京:中国农业出版社, 2002. LV Hongxing, PEI Guoxia, YANG Lingxia. Hydraulics[M]. Beijing:China Agriculture Press, 2002. (in Chinese)