周义仁，石先德，张 茜

(太原理工大学水利科学与工程学院，太原 030024)

灌区精准量水是提高农业用水利用率的基础所在[1]，本文针对灌区量水过程中存在的测流精度低、水头损失大等问题提出了一种新型量水装置----摆杆式明渠测流装置，其工作原理需建立摆杆摆动角度和渠道流量之间关系进行测流，探究该装置在不同类型渠道和测流摆杆下的测流特性具有重要意义。随着数值模拟技术的发展，计算流体力学所能计算的领域越来越多[2]，通过数值模拟方法建立数学模型解决实际工程问题具有投入少、研究周期短、可调参数多等优点。FLUENT是目前应用最为广泛的计算流体力学软件，其计算方法及后处理功能十分先进[3]。摆杆式明渠测流装置进行测流任务时，测流摆杆由于水流冲击而摆动，是属于固体边界在流域中产生的被动运动，在FLUENT中，动网格模型可以用来模拟由于流域边界运动引起流域形状随时间变化的流动情况，也可以用来模拟由于流场改变而导致固体边界发生被动运动情况[4,5]。利用FLUENT软件对摆杆式明渠测流装置测流特性进行研究，可以提高该装置的测流精度及适用性。

1 摆杆式明渠测流装置简介

摆杆式明渠测流装置机械结构如图1所示。

1-渠道边坡；2-调平旋钮；3-可伸缩不锈钢支架；4-超声波探头；5-液晶显示器；6-供电及数据处理模块；7-水平仪；8-角度传感器；9-测流摆杆图1 摆杆式明渠测流装置

圆柱形摆杆放置于流动的水中时受到渠道水流阻力和升力作用，会沿水流方向摆动，通过建立渠道流量和摆杆摆动角度之间的关系，测得摆杆旋转角度便可以得知渠道瞬时流量。

2 数值模拟研究

2.1 计算模型及控制方程

计算采用VOF自由液面模型，对于水气两相流动，VOF模型对其进行描述时，先定义函数qw(x,y,z,t)和qa(x,y,z,t)，它们分别表示水与空气在某一时刻的某一单元格内的体积分数。而在每一个单元格内[6]：

qw+qa=1

(1)

通过公式(1)能够得到，当qw=1时，表示在该单元格内水的体积分数为1，也就是这个单元内充满了水；反之，qa=1时，则表示这个单元内全都是空气；若0

(2)

式中：t为时间；ui为速度分量；xi为坐标分量。

湍流模型选用Realizablek-ε模型，连续性方程、动量方程、湍流动能k方程、湍流耗散率ε方程如下所示[8]：

(3)

(4)

Gk+Gb-ρε-Ym+Sk

(5)

(6)

2.2 几何建模及网格划分

本文模拟及试验所用渠道为D50型U形渠道，其主要参数为：渠深0.5 m，渠道底弧半径0.25 m，底弧圆心角163°，渠道倾角8.5°，渠顶宽度0.58 m，渠道底坡为1/1 000。为了充分模拟U形渠道水流流态并尽量减少计算机模拟时间，渠道长度设为4 m，入口处1 m设置长方体稳流池，确保水流从入口处平稳过渡到U形渠道中。设备放置于U形渠道中央，距出口2 m位置，取摆杆顶端圆心为坐标原点，水流方向为Z轴正方向，竖直向上为Y轴正方向，顺水流方向左岸为X轴正方向，如图2所示。

图2 坐标轴示意图

结构化网格与非结构化网格联合进行求解，可以很好地解决计算缓慢以及网格质量差等问题，使得模拟更加精准和高效[9,10]。计算模型在最大尺度以及最小尺度之间存在较大差异(渠道长达4 m，测流摆杆直径为5 mm)，所以在摆杆上下游0.5 m以内采用非结构化网格，在摆杆运动范围之外采取结构化网格。设置全局网格最大尺寸为4 cm，在摆杆周围进行非结构化网格划分，并对摆杆进行局部加密处理，摆杆表面最大网格尺寸为1 mm，网格总数约为36万个，网格划分如图3所示。

图3 网格划分示意图

2.3 边界条件及求解方法

为了更加精确的模拟渠道水流流动对摆杆式明渠测流精度的影响，上游设置进口为速度入口(VELOCITY-INLET)，利用平均流速来控制渠道瞬时流量大小；下游出口设置为压力出口(PRESSURE-OUTLET)；渠道顶部设为压力进口(PRESSURE-INLET)，受到大气压强作用，顶部参考压强为101 kPa；摆杆以及渠道侧面全部设为无滑移固体边壁，渠道壁面为不锈钢材质，设置壁面粗糙度高度为0 cm，粗糙度常数为0.5。整个计算域中最开始时充满了空气，通过VOF方法进行迭代，自动生成水气交界面。

本文主要采用有限体积法进行计算，主要步骤是进行计算区域划分后，进行区域的离散化与控制方程的离散，计算区域的离散化实质上便是上述所提到的网格划分。控制方程的离散化可以利用如下所示的通用形式表示[11]：

(7)

式中：ρ为流体密度；u为液体流速；t为时间；φ为广义变量，可以为速度、温度或者浓度等一些待求物理量；Γ是相应于φ的广义扩散项；Sk是广义源项。

建立了可以进行计算的代数方程组之后，需要对各个未知量的求解顺序以及方式进行特殊处理，对于瞬态计算，采用PISO算法比SIMPLE或者SIMPLEC算法更加合适。

2.4 动网格设定

在FLUENT中动网格更新方法主要有3种，包括弹簧近似光滑模型(Spring-based Smoothing)、动态分层模型(Dynamic Layering)以及局部网格重构模型(Local Remeshing)[12]。本文主要采用弹簧光顺与局部网格重构相结合进行计算。

动网格流场控制模型中，计算区域是变化的，其守恒方程为[13]：

(8)

FLUENT中6DOF(六自由度)模型主要是用于模拟计算域中的刚体构件受到流体的作用之后运动的轨迹及状态，通过模拟流体的运动对刚体产生的作用力，来推算刚体的运动轨迹。其中刚体的质量、位移、旋转以及转动惯量可以通过UDF(用户自定义函数)来进行编写录入[14]。

本文主要研究3种不同材质的测流摆杆在D50型U形渠道下的测流特性，不同材质的摆杆参数如表1所示。

3 结果分析

3.1 水面线分析

对明渠进行三维数值模拟时，水面线是一项基本研究对象，通过试验与模拟结果的比较，可以验证数值模拟的可行性。

表1 不同材质测流摆杆参数

对不同流量下D50型U形渠道水流进行模拟可以得到测流装置下方(Z=0 m)水位随流量变化图，如图4所示。

图4 水位与流量关系图

当渠道瞬时流量为300 m3/h时，采用铝材测流摆杆，t=0 s时开始从入口注水，渠道水流从入口到出口水面线随时间变化瞬态图如图5所示，图中蓝色区域为空气，红色区域为水流，红蓝相间区域为水气交界面，白色区域为摆杆。观察图5可以发现，随着入口处水流不断涌入渠道，渠道中水位逐渐上升，当t=6.87 s以后，水位趋于平缓。稳流池与渠道相接处水面略有壅高，这是由于稳流池出口截面较小，水流撞击壁面后造成水位上升。测流设备上下游1 m水流平缓，在渠道出口处水位有所下降。由图4与图5可以发现，利用数值模拟对U形渠道进行仿真模拟得到的渠道水位与试验结果相近，所以利用VOF方法以及Realizablek-ε湍流模型对渠道水流进行数值模拟可以得到与实际相符的结果。

图5 300 m3/h流量下渠道沿程水面线变化图

3.2 流速分布

为了解U形渠道过水断面流速变化特征，对测流设备下方(Z=0 m)测流断面进行流速分析，如图6所示为不同流量下，渠道过水断面沿水流方向(Z轴正方向)水流流速分布云图。

由图6可以发现，在底坡为1/1 000的D50型U形渠道中，随着渠道内流量的增加，过水断面面积增大，水位逐渐升高，水流流速逐渐增大，当流量为150 m3/h时，最大流速为0.96 m/s，当流量为600 m3/h时，最大流速达到1.37 m/s。对于渠道断面沿水流方向的流速分布而言，在靠近渠道边壁的流速急剧下降，在渠道边壁区域接近于0 m/s，其最大流速位于水面下方，在气液交界面处，流速较缓，这与渠道实际流速分布规律相吻合。

3.3 动网格变化分析

对于动网格的计算方法，本文主要采用弹簧光顺与局部网格重构相结合进行计算。弹簧光顺将网格节点视为可变形的弹簧，固体边界发生位置变化时，将会挤压或者拉伸网格，通过求解弹簧的应力平衡方程得到新网格的节点位置，但是其只能处理网格微小变形的情况，当网格变形较大时，则需要利用局部网格重构技术。针对网格尺寸以及畸变率，当网格出现较大程度的变形时，将会对该区域的网格进行网格重构，通过守恒定律以及插值映射从旧网格中得到物理量的大小。如图7所示，以不锈钢材质摆杆在不同流量的水流作用下产生网格重构的变化为例，对动网格变化进行分析。图7蓝色区域为空气，红色区域为水流，红蓝相间区域为水气交界面。

图6 不同流量下断面流速分布云图

图7 不锈钢材质摆杆局部网格变化

由图7可以发现，当摆杆受到水流作用向下游摆动时，摆杆周围的网格发生变形，下游网格受到摆杆挤压，体积变小，而摆杆上游位置的网格体积增大。随着流量的增大，摆杆摆动的角度随之增大，动网格更新的区域也随之增大，且之前摆杆经过的区域网格变得更加密集。在整个迭代计算过程中，网格体积始终没有出现负体积。通过模拟研究发现，铝材以及有机玻璃两种材质的测流摆杆与图7中不锈钢材质摆杆周围网格变化规律相似，说明采用FLUENT中6DOF动网格模型模拟摆杆式测流装置测流特性是可行的。

3.4 测流公式及精度分析

摆杆式明渠测流装置测得流量Q与以下因素有关：水的密度 、重力加速度g、测流摆杆密度 、测流摆杆长度L、测流摆杆直径D、测流摆杆被水推动的水平位移d(d=Lsinθ)、水中泥沙浓度S、水动力黏滞系数v、渠道糙率n以及渠道底坡i。摆杆式明渠测流装置的测流计算公式可以通过量纲分析法推导得到。选取d，g，ρ杆3个相互独立的量作为研究基本物理量，最终可以得到以下测流公式：

(9)

(10)

Q=Ksin2.5θ

(11)

对D50型U形渠道下不同材质的测流摆杆测流特性进行数值模拟，得到3种材质摆杆在不同流量下相应摆动角度，利用线性拟合得到式(11)中系数K的大小，便得到了摆杆式明渠测流装置在该渠道下的测流公式，如表2所示。

表2 不同材质测流摆杆测流公式

图8为模拟值与试验值对比图，主坐标轴反映了不同材质的测流摆杆sin2.5θ与渠道流量Q之间的关系，次坐标轴反映了利用表2算得流量与实际流量之间的相对误差。由图8可以发现，不同材质的测流摆杆在U形渠道下，渠道流量Q与摆杆摆动角度sin2.5θ之间成正比例函数关系，模拟值与试验值变化一致。利用测流公式得到渠道流量与真实流量之间相对误差α较小，其中不锈钢材质的测流摆杆测得流量最大相对误差为5.33%，平均相对误差为2.5%；铝材测流摆杆测得流量最大相对误差为5.55%，平均相对误差为1.87%；有机玻璃材质测流摆杆测得流量最大相对误差为5.26%，平均相对误差为2.22%。

4 结 语

利用FLUENT对摆杆式明渠测流装置在底坡为1/1 000的D50型U形渠道下测流特性进行数值模拟研究，通过数值模拟以及试验数据对比分析得到以下结论。

图8 不同材质摆杆流量与角度关系图

(1)利用FLUENT中VOF模型与Realizablek-ε湍流模型对D50型U形渠道下水流特性进行模拟研究是可行的，模拟渠道沿程水面线与试验相符。渠道过水断面上流速分布规律为：在靠近渠道边壁的区域水流流速急剧下降，在边壁旁趋近于零，而水流最大流速并不是位于水面位置，而是在水面下方的位置，这也与国内外其他学者研究结论一致。

(2)利用6DOF动网格模型对测流摆杆摆动特性进行模拟是可行的，通过水流与摆杆之间的耦合计算可以得到摆杆的摆动角度与渠道流量大小之间的变化规律，随着渠道内瞬时流量的增大，摆杆摆动的角度也逐渐增大，摆杆周围非结构化网格产生重构。

(3)利用数值模拟对摆杆式测流装置测流精度进行研究，发现量纲分析得到的测流公式Q=Ksin2.5θ中系数K大小随着摆杆密度的减小而减小，测得流量与实际流量之间平均相对误差小于3%，可以满足明渠测流精度要求。

本文对摆杆式明渠测流装置测流特性进行数值模拟研究，模拟结果与实际情况一致，充分验证了该测流装置的可行性，对提高该装置测流精度有重要意义，为灌区精准量水的发展提供了新思路。