輸水明渠壅水狀態斷面流速特性三維仿真研究

王靜茹，管光華，靳偉榮，熊 驥

（1.武漢大學水資源與水電工程科學國家重點實驗室，武漢430072；2.湖北省國際灌排研究培訓中心，武漢430072）

0 引言

灌區量測水設施是實施計劃用水、合理配水的重要手段；是農業用水計量收費的重要依據和保證；也是水資源管理三條紅線的重要技術支撐，是建設現代化灌區重要的一環。傳統的量水方法有堰槽測流、流速儀測流、超聲波測流、水工建筑物測流等。由于我國灌區量水的需求與日俱增，近些年來涌現了大量的新型量水技術和方法，如雷達波流速儀測流、多點超聲波測流等。

目前在實際工程中使用標準斷面進行量水，常忽略了渠道壅水高度對于水位流量關系的影響，即認為標準斷面的水位-流量是固定的［1，2］。而使用雷達技術測流只能得到渠道水面流速，要得到斷面平均流速還需乘以一定的水面流速系數，因此水面流速系數的取值對于這一類流速儀的精度具有決定性的作用。因此，本文主要研究明渠在不同下游壅水高度條件下水位-流量關系和水面流速系數的變化規律，結合我國灌區內渠道主要斷面形式，本文選取的渠段均為具有梯形斷面的緩坡明渠。

自然界中大多數的明渠水流都是紊流，對于明渠紊流的運動規律和水力特性許多學者進行了一定的探索和研究。Keu?lengan［3］將平板邊界層的研究與明渠水流相結合，提出明渠斷面垂線上沿渠道方向的流速滿足對數分布。后來有學者［4］認為明渠水流可沿水深分成內區和外區。內區由黏性底層、過渡層和對數區組成，外區用Coles［5］提出的尾流函數對對數流速分布進行修正。Prandtl 等［6］和Afzal［7］給出了指數形式的流速公式，認為主流區的流速梯度和水分子的黏性呈漸進指數關系。由一些學者的研究成果［8-10］知，矩形斷面和梯形斷面順水流方向流速-u分布如圖1所示，最大流速發生在水面以下；矩形明渠過水斷面上水平流速-v和垂直流速-w分布圖如圖2所示，可以看出有表面渦流和底部渦流，即有二次流的存在。

圖1 典型斷面順水流方向流速分布圖Fig.1 Velocity distribution diagram of typical cross-section along the flow direction

圖2 矩形明渠和分布圖（二次流）Fig.2 Velocity distribution diagram at horizontal direction and vertical direction

現有的明渠流速分布規律是針對恒定均勻流的，但實際渠道中糙率的不均勻分布，泥沙淤積，下游壅水等都會影響渠道中的流速分布。由于輸配水渠道系統中常用節制閘進行水位、流量調節，渠道在臨近節制閘處會產生不同程度的壅水，尤其是在渠道高水位低流量運行時。目前少有文獻研究下游壅水高度對過水斷面的流速分布的影響，因此本文選取漳河三干渠三分干，基于FLOW-3D 三維仿真，對輸水明渠在下游壅水時的斷面流速分布進行研究。

1 數學模型

1.1 控制方程

FLOW-3D中采用的控制方程為N-S方程，并且在求解N-S方程時采用雷諾平均法，在控制方程中分別加入面積分數和體積分數［11］。

其連續性方程為：

動量方程為：

式中：i、j=1、2、3；ui為x、y、z方向的速度；Ax、Ay、Az為計算單元x、y、z方向面積；VF為各計算單元內液體的體積分數；ρ為液體密度；P為壓強；gi為重力；fi為雷諾應力［12，13］。

1.2 湍流模型

RNGk-ε模型是基于重整化群（Renormalization Group）的理論提出來的標準模型的修正方程。RNGk-ε模型相對于k-ε模型，可以更好地處理高應變率及流線彎曲程度較大的流動［14-16］。因此本文采用RNGk-ε模型。RNGk-ε模型的湍動能k和湍動能耗散率ε方程分別為：

湍動能方程：

湍動能耗散率方程：

式中：k為湍動能，m2/s2；ε為湍動能耗散率，kg·m2/s2；μ為流體動力黏滯系數，N·s/m2；μt為流體的湍動黏度，μt=ρCμk2/ε，N·s/m2；αε，αk，C1ε和C2ε為經驗常數，βη3），其中=0.012，C1ε=1.42；常數C2ε=1.68；Gk為由平均速度梯度引起的紊動能產生項，可由式定義［14，17］。

1.3 自由表面

在數值模擬中，自由表面就是氣體和液體的交界面，氣體只對液體施加壓力。FLOW-3D采用truVOF方法追蹤自由表面流動。在水氣兩相流中，定義函數αw和αa分別為計算單元中水和氣的體積分數。每個單元中，水和氣的體積分數之和為1，即αw+αa=1。當αw=0 時，計算單元內全為氣相；當αw=1 時，計算單元內全為液相；當0<αw<1 時，計算單元內部分是水，部分是氣，包含自由表面［18］。

2 三維流場建模

漳河灌區三干渠三分干全長36.8 km，本次建模主要選擇從三干渠三分干分水閘到下游三分干卞廟節制閘的一段，全長9 300 m。該渠段的幾何參數及水力特性為：底坡i=1/8 000，底寬b=4 m，邊坡系數m=1.75，糙率n=0.015，渠道設計流量10 m3/s，渠道加大流量為13 m3/s。

2.1 模型范圍及網格劃分

2.1.1 模型范圍

模型范圍取其中的100 m，按原型1∶1建立幾何實體模型。

2.1.2 網格無關性檢驗

網格質量的好壞直接影響到數值計算的穩定性和精度。選擇合適的網格數目對于數值計算的準確性以及縮短計算時間有非常重要的影響［19-21］。本文取10 m 長渠道，對0.2，0.1，0.05，0.03 m 四種網格寬度的均勻網格進行無關性檢驗。網格形狀為正方體。數值模擬結果如表1所示。不同網格寬度穩定后沿程水深如圖3所示。

圖3 不同網格寬度穩定后沿程水深Fig.3 The steady water depth under different mesh sizes

由表1可知：當網格尺寸大于0.05 m時，網格寬度對模擬結果影響較大，此時網格尺寸越小，模擬結果越接近于恒定均勻流，但相應運行時間會大幅增加；當網格尺寸小于0.05 m 時，網格寬度對模擬結果影響較小。因此出于時間和精確度的考慮，本文網格寬度取0.05 m。

表1 10 m長渠道不同網格寬度模擬結果Tab.1 Simulation results of a 10 m long channel under different mesh sizes

2.1.3 網格劃分

當模型渠道長度為100 m，網格寬度為0.05 m 時，總網格數可達到千萬以上，因此考慮用相連網格塊，在加密網格的同時盡可能減少成指數倍增長的網格數。且模擬結果表明，從水面到渠底依次用0.05、0.1、0.2 m 的網格寬度與計算區域均用0.05 m的網格寬度沿程水深及流速分布大致相同。

因此計算區域用3 個連接的網格塊來劃分，靠近渠底的網格塊單元尺寸為0.2 m×0.2 m×0.2 m，渠道中部網格塊單元尺寸為0.1 m×0.1 m×0.1 m，靠近水面處網格塊尺寸為0.05 m×0.05 m×0.05 m，網格總數為500～600萬，具體劃分見圖4。

圖4 網格劃分示意圖Fig.4 Sketch of mesh discretization

2.2 物理模型、邊界條件、初始條件

除了紊流模型選擇RNGk-ε模型外，還需要在軟件中選擇重力模型和黏性流動。重力模型中設定重力加速度為9.8 m/s2。

上游邊界條件為固定流量邊界條件（不設定流體高度，默認流體從整個邊界開放區域流入，流動方向與邊界垂直）；下游邊界為固定水位的壓力出流邊界條件，水深hd=h0+Δh，其中h0是正常水深1.76 m，Δh為壅水高度，分別取0、10、20、30、40、50 cm；渠道兩側及底部為無滑移的wall邊界條件；上部為壓力為0的壓力邊界條件。

本文選取沿程1.76 m 的正常水深作為初始條件，用CAD 建立水深模型，輸出為stl格式，再將其作為流體導入。

2.3 單位、流體特性及時間設置

選擇SI 單位制，流體為20 ℃的水。初始時間步長設置為0.001 s，最小時間步長設置為10-7s。時間步長的控制方法使用Stability and convergence 和stability，步長會自適應調節。

2.4 模型驗證

下游壅水高度Δh=0 時，模擬的是恒定均勻流，數值模擬的結果與用HEC-RAS 計算和曼寧公式計算的沿程水面線對比如圖5所示。由于上游邊界條件的影響，用FLOW-3D 模擬出來的沿程水深在進口處會有0.2 m 的降落，很快穩定在正常水深1.76 m 附近，HEC-RAS 和曼寧公式計算出沿程水深均為1.76 m。FLOW-3D 得到的水面線誤差小于1 mm，因此FLOW-3D 數值計算的結果具有一定可信度。

圖5 數值模擬水深和HEC-RAS恒定均勻流水深Fig.5 Simulation results of water depth by flow3d and HEC-RAS

3 仿真及結果分析

3.1 渠道流速分布

國內外學者根據梯形斷面渠道、矩形斷面渠道、復式斷面河槽等的室內試驗和現場資料證實，渠道中心自由水面處的流速不是最大值，最大流速發生在自由水面以下。一些學者提出是由于過水斷面存在二次流［22］。

圖6（a）～（f）為距離上下游邊界均50 m 的斷面沿水流方向（y方向）的流速分布圖。本文同樣得到最大流速發生在自由水面以下。同時對過水斷面流速等值線圖分析可得：①渠道中垂線兩側流速分布大致對稱相等；②在流量相同，下游水深不同的情況下，隨著下游水深增大，即壅水程度的增大，最大流速點的位置相對水面向下降，且最大流速值減小；③渠道邊壁處流速等值線圖除靠近水面處向渠道中心彎曲外，其他地方近似與邊壁輪廓平行。

圖6 距離上下游邊界均50 m的斷面沿水流方向（y方向）的流速分布圖Fig.6 Distribution of velocity along y-axis at the cross-section 50 m away from upstream and downstream end

圖7（a）～（f）為距離上下游邊界均50 m 的斷面沿水平方向（x方向）的流速分布圖。結果表明：①邊壁附近的水流有向中心聚集的趨勢；②渠道中垂線兩側流速分布大致對稱相等，且渠道中垂線附近流速大致為零；③在流量相同，下游水深不同的情況下，隨著下游水深增大，即壅水程度的增大，此斷面上水平方向流速最大值點的位置相對水面向下降，且最大流速值減小。

圖7 距離上下游邊界均50 m的斷面沿水平方向（x方向）的流速分布圖Fig.7 Distribution of velocity along x-axis at the cross-section 50 m away from upstream and downstream end

這可以用二次流來解釋，由于邊壁和自由水面的影響，在水面附近存在一個指向中心的二次流，將邊壁附近的低速水流帶到中心，這也是造成表面流速降低的原因之一［22］。

3.2 點流速測流的精度分析

使用雷達技術測流時，需要將雷達測得的表面點流速乘以水面流速系數k得到斷面平均流速，然后根據過水斷面面積得到流量：

式中：vc表示斷面平均流速；vs表示水面平均流速。

通常情況下，在使用雷達測流之前會通過率定的方式確定一個固定的水面流速系數，設為k0。但渠道在實際運行過程中，由于受到壅水作用的影響，在流量不變的情況下，隨著下游水深變化，水面流速系數和斷面流速分布均會發生變化。

以本文所取漳河典型渠段為例，由于其底坡較小，因此下游邊界水深改變導致的壅水影響范圍較大。為研究相同流量下，下游水深不同時，相同斷面處的流速分布和水面流速系數變化規律，取渠段內距上下游邊界均50 m 的斷面，研究不同下游壅水條件下采用不壅水時的水面流速系數k0測流導致的流量誤差。

表2 中h表示斷面處水深，A表示過水面段面積，Q1表示實際流量，由于仿真中設置入口處為固定流量10m3/s 并且水流為恒定流，因此各斷面處流量均為10 m3/s，Q2表示根據k0計算得到的流量：

表2 Δh不同時使用固定的k0產生的測流誤差Tab.2 Discharge measurement error using constant k0 under different values of Δh

由表2 可知，由于渠道底坡較小，下游水深邊界改變，此斷面水深幾乎與下游邊界相等，只小于下游水深1～5 mm，且相較于下游水深的減小值隨下游水深的增大而增大。在該渠段內，如果忽略壅水作用對水面流速系數的影響，而采用固定的水面流速系數，所測量得到的流量誤差隨壅水高度的增加而增大，在壅水高度達到40 cm 時誤差大于10%，不滿足量水精度的要求。

3.3 水面流速系數

3.3.1 水面流速系數變化規律

本文對不同下游邊界條件恒定非均勻流時離上游邊界40、50、60、70 m 斷面的水面流速系數進行了研究，如圖8所示。為減小誤差，每個水面平均流速的值均是選取穩定后3 個時刻數據的平均值。由圖8 可知，下游水深條件改變時渠道沿水流方向各斷面水面流速系數并不為一定值。下游水深變化量Δh=0、10、20 cm 時沿水流方向各斷面的水面流速系數沿程略有增加，大致為1；Δh=30 cm 時沿程水面流速系數基本不變，為1.05；Δh=40 cm 或50 cm 時沿程水面流速系數從1.2～1.21 變小到1.1～1.12。

圖8 距上游邊界不同距離的斷面的水面流速系數Fig.8 The values of surface velocity coefficient at different cross-sections

考慮到實際渠道運行過程中，下游水深邊界改變，離上游邊界距離的不同，最直接影響的物理量是水深。如圖9所示，水深不變或增加10、20 cm 時，離上游邊界40、50、60、70 m 的斷面處水面流速系數基本維持在1附近，之后水深再增加時，水面流速系數也增大，水深從1.96 m 增大到2.16 m 時，水面流速系數增大地較快，水深從2.16 m再增大時，水面流速系數基本不變。

圖9 各個斷面上水面流速系數與水深之間的關系Fig.9 The relation between water depth and surface velocity coefficient at different cross-sections

因此在渠道實際運行過程中，Δh在20 cm 范圍內時，用正常水深情況下率定的水面流速系數測流比較可靠，Δh大于20 cm 時，實際水面流速系數大于率定的水面流速系數，測量的流量會小于實測流量，造成一定的流量偏差。

3.3.2 考慮水深變化對水面流速系數進行修正

本文選取三次多項式對距上游邊界50 m 斷面處的水面流速系數和下游邊界處的水深關系進行擬合，得到修正后水面流速系數k'與下游水深hd的表達式：

擬合曲線與三維仿真結果的決定系數R2=0.99，故可采用該式計算不同下游水深下該斷面的水面流速系數，進而求斷面流量。

由表3 可見，在本文所采用的工況下，k'值所推算出來的斷面流量與實際流量誤差小于1.4%。

表3 Δh不同時同一斷面的水面平均流速和流量Tab.3 The average surface velocity and flow rate at the same section under different Δh

表3 中h表示斷面處水深，A表示過水面段面積，Q1表示實際流量，由于仿真中設置入口處為固定流量10 m3/s并且水流為恒定流，因此各斷面處流量均為10 m3/s，表示根據k'計算得到的流量：

3.3.3 多點流速測量

雷達測流采用的是無線電多普勒效應。得到的水面數據的個數與水面波浪和漂浮物有關。因此灌區實際測流并不能像數值模擬一樣得到水面很多點的數據，只能測出水面有限點的流速。

本文對測點的個數與測流誤差之間的關系進行了一定的討論。測量水面上三點（離水面中心點距離0 m、±2m）、五點（離水面中心點距離0 m、±2 m、±4 m）、七點（離水面中心點距離0 m、±1.5 m、±3 m、±4.5 m）的流速得到的水面平均流速分別為v3、v5、v7，則測流產生的誤差如表4所示。

表4 多點流速測流產生的誤差Tab.4 Discharge measurement error using multiple point measurement

表中Q3，Q5，Q7為根據v3、v5、v7和k'計算得到的流量，以Q3為例：

由表4 可知，水面測量的點流速越多，測流產生的誤差越小。只測量水面3 個點的流速時，在6 種工況下產生的誤差約在10%～19%范圍內；五點流速的測流誤差在5%～12%范圍內；七點流速的測流誤差在0%～9%范圍內。因此測量表面流速的個數為7個時，測流誤差可以接受，且壅水高度Δh<20 cm 時，測流誤差在5%以內，誤差相對較小。

4 結論

本文采用三維仿真軟件對輸水明渠不同下游邊界條件下渠道流速分布進行了三維數值模擬，得出結論如下。

（1）不同壅水高度下渠道同一橫斷面水面流速系數會顯著改變。若忽略渠道壅水特性，認為水面流速系數固定，由測量出的水深計算出過流面積，再得到流量可能會造成顯著的測量誤差。根據本文計算的算例，測流方法引入的誤差最大可為16.7%。

（2）當所取渠段下游水深邊界變化Δh在20 cm 范圍內時，用點流速和斷面平均流速的關系率定好的水面流速系數測流比較可靠，Δh大于20 cm 時，實際水面流速系數會大于率定的水面流速系數，測量的流量會小于實測流量，造成一定的流量偏差。

（3）考慮壅水導致的水深變化而對水面流速系數k修正時，測流誤差會比用固定水面流速系數明顯減小，最大測流誤差可從16.7%減小到1.4%。

（4）實際用雷達技術測表面流速時，設置的水面流速測點越多，流量測量結果越準確。若可測水面7個點的流速時，測流誤差在9%以內，在渠道量水規范［23］推薦的精度可接受范圍內。

本文只選取了一個典型渠段且渠段長度僅有100 m，該結論是否具有尺度效應還需選取不同渠道特性的更長的渠段進行研究，在此基礎上還可對水面流速系數與糙率、斷面尺寸比、底坡等的關系進行進一步的研究。□