輸水明渠壅水狀態(tài)斷面流速特性三維仿真研究

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

（1.武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢430072；2.湖北省國(guó)際灌排研究培訓(xùn)中心，武漢430072）

0 引言

灌區(qū)量測(cè)水設(shè)施是實(shí)施計(jì)劃用水、合理配水的重要手段；是農(nóng)業(yè)用水計(jì)量收費(fèi)的重要依據(jù)和保證；也是水資源管理三條紅線的重要技術(shù)支撐，是建設(shè)現(xiàn)代化灌區(qū)重要的一環(huán)。傳統(tǒng)的量水方法有堰槽測(cè)流、流速儀測(cè)流、超聲波測(cè)流、水工建筑物測(cè)流等。由于我國(guó)灌區(qū)量水的需求與日俱增，近些年來(lái)涌現(xiàn)了大量的新型量水技術(shù)和方法，如雷達(dá)波流速儀測(cè)流、多點(diǎn)超聲波測(cè)流等。

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

自然界中大多數(shù)的明渠水流都是紊流，對(duì)于明渠紊流的運(yùn)動(dòng)規(guī)律和水力特性許多學(xué)者進(jìn)行了一定的探索和研究。Keu?lengan［3］將平板邊界層的研究與明渠水流相結(jié)合，提出明渠斷面垂線上沿渠道方向的流速滿足對(duì)數(shù)分布。后來(lái)有學(xué)者［4］認(rèn)為明渠水流可沿水深分成內(nèi)區(qū)和外區(qū)。內(nèi)區(qū)由黏性底層、過(guò)渡層和對(duì)數(shù)區(qū)組成，外區(qū)用Coles［5］提出的尾流函數(shù)對(duì)對(duì)數(shù)流速分布進(jìn)行修正。Prandtl 等［6］和Afzal［7］給出了指數(shù)形式的流速公式，認(rèn)為主流區(qū)的流速梯度和水分子的黏性呈漸進(jìn)指數(shù)關(guān)系。由一些學(xué)者的研究成果［8-10］知，矩形斷面和梯形斷面順?biāo)鞣较蛄魉?u分布如圖1所示，最大流速發(fā)生在水面以下；矩形明渠過(guò)水?dāng)嗝嫔纤搅魉?v和垂直流速-w分布圖如圖2所示，可以看出有表面渦流和底部渦流，即有二次流的存在。

圖1 典型斷面順?biāo)鞣较蛄魉俜植紙DFig.1 Velocity distribution diagram of typical cross-section along the flow direction

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

現(xiàn)有的明渠流速分布規(guī)律是針對(duì)恒定均勻流的，但實(shí)際渠道中糙率的不均勻分布，泥沙淤積，下游壅水等都會(huì)影響渠道中的流速分布。由于輸配水渠道系統(tǒng)中常用節(jié)制閘進(jìn)行水位、流量調(diào)節(jié)，渠道在臨近節(jié)制閘處會(huì)產(chǎn)生不同程度的壅水，尤其是在渠道高水位低流量運(yùn)行時(shí)。目前少有文獻(xiàn)研究下游壅水高度對(duì)過(guò)水?dāng)嗝娴牧魉俜植嫉挠绊?，因此本文選取漳河三干渠三分干，基于FLOW-3D 三維仿真，對(duì)輸水明渠在下游壅水時(shí)的斷面流速分布進(jìn)行研究。

1 數(shù)學(xué)模型

1.1 控制方程

FLOW-3D中采用的控制方程為N-S方程，并且在求解N-S方程時(shí)采用雷諾平均法，在控制方程中分別加入面積分?jǐn)?shù)和體積分?jǐn)?shù)［11］。

其連續(xù)性方程為：

動(dòng)量方程為：

式中：i、j=1、2、3；ui為x、y、z方向的速度；Ax、Ay、Az為計(jì)算單元x、y、z方向面積；VF為各計(jì)算單元內(nèi)液體的體積分?jǐn)?shù)；ρ為液體密度；P為壓強(qiáng)；gi為重力；fi為雷諾應(yīng)力［12，13］。

1.2 湍流模型

RNGk-ε模型是基于重整化群（Renormalization Group）的理論提出來(lái)的標(biāo)準(zhǔn)模型的修正方程。RNGk-ε模型相對(duì)于k-ε模型，可以更好地處理高應(yīng)變率及流線彎曲程度較大的流動(dòng)［14-16］。因此本文采用RNGk-ε模型。RNGk-ε模型的湍動(dòng)能k和湍動(dòng)能耗散率ε方程分別為：

湍動(dòng)能方程：

湍動(dòng)能耗散率方程：

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

1.3 自由表面

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

2 三維流場(chǎng)建模

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

2.1 模型范圍及網(wǎng)格劃分

2.1.1 模型范圍

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

2.1.2 網(wǎng)格無(wú)關(guān)性檢驗(yàn)

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

圖3 不同網(wǎng)格寬度穩(wěn)定后沿程水深Fig.3 The steady water depth under different mesh sizes

由表1可知：當(dāng)網(wǎng)格尺寸大于0.05 m時(shí)，網(wǎng)格寬度對(duì)模擬結(jié)果影響較大，此時(shí)網(wǎng)格尺寸越小，模擬結(jié)果越接近于恒定均勻流，但相應(yīng)運(yùn)行時(shí)間會(huì)大幅增加；當(dāng)網(wǎng)格尺寸小于0.05 m 時(shí)，網(wǎng)格寬度對(duì)模擬結(jié)果影響較小。因此出于時(shí)間和精確度的考慮，本文網(wǎng)格寬度取0.05 m。

表1 10 m長(zhǎng)渠道不同網(wǎng)格寬度模擬結(jié)果Tab.1 Simulation results of a 10 m long channel under different mesh sizes

2.1.3 網(wǎng)格劃分

當(dāng)模型渠道長(zhǎng)度為100 m，網(wǎng)格寬度為0.05 m 時(shí)，總網(wǎng)格數(shù)可達(dá)到千萬(wàn)以上，因此考慮用相連網(wǎng)格塊，在加密網(wǎng)格的同時(shí)盡可能減少成指數(shù)倍增長(zhǎng)的網(wǎng)格數(shù)。且模擬結(jié)果表明，從水面到渠底依次用0.05、0.1、0.2 m 的網(wǎng)格寬度與計(jì)算區(qū)域均用0.05 m的網(wǎng)格寬度沿程水深及流速分布大致相同。

因此計(jì)算區(qū)域用3 個(gè)連接的網(wǎng)格塊來(lái)劃分，靠近渠底的網(wǎng)格塊單元尺寸為0.2 m×0.2 m×0.2 m，渠道中部網(wǎng)格塊單元尺寸為0.1 m×0.1 m×0.1 m，靠近水面處網(wǎng)格塊尺寸為0.05 m×0.05 m×0.05 m，網(wǎng)格總數(shù)為500～600萬(wàn)，具體劃分見(jiàn)圖4。

圖4 網(wǎng)格劃分示意圖Fig.4 Sketch of mesh discretization

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

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

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

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

2.3 單位、流體特性及時(shí)間設(shè)置

選擇SI 單位制，流體為20 ℃的水。初始時(shí)間步長(zhǎng)設(shè)置為0.001 s，最小時(shí)間步長(zhǎng)設(shè)置為10-7s。時(shí)間步長(zhǎng)的控制方法使用Stability and convergence 和stability，步長(zhǎng)會(huì)自適應(yīng)調(diào)節(jié)。

2.4 模型驗(yàn)證

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

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

3 仿真及結(jié)果分析

3.1 渠道流速分布

國(guó)內(nèi)外學(xué)者根據(jù)梯形斷面渠道、矩形斷面渠道、復(fù)式斷面河槽等的室內(nèi)試驗(yàn)和現(xiàn)場(chǎng)資料證實(shí)，渠道中心自由水面處的流速不是最大值，最大流速發(fā)生在自由水面以下。一些學(xué)者提出是由于過(guò)水?dāng)嗝娲嬖诙瘟鳎?2］。

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

圖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方向）的流速分布圖。結(jié)果表明：①邊壁附近的水流有向中心聚集的趨勢(shì)；②渠道中垂線兩側(cè)流速分布大致對(duì)稱相等，且渠道中垂線附近流速大致為零；③在流量相同，下游水深不同的情況下，隨著下游水深增大，即壅水程度的增大，此斷面上水平方向流速最大值點(diǎn)的位置相對(duì)水面向下降，且最大流速值減小。

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

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

3.2 點(diǎn)流速測(cè)流的精度分析

使用雷達(dá)技術(shù)測(cè)流時(shí)，需要將雷達(dá)測(cè)得的表面點(diǎn)流速乘以水面流速系數(shù)k得到斷面平均流速，然后根據(jù)過(guò)水?dāng)嗝婷娣e得到流量：

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

通常情況下，在使用雷達(dá)測(cè)流之前會(huì)通過(guò)率定的方式確定一個(gè)固定的水面流速系數(shù)，設(shè)為k0。但渠道在實(shí)際運(yùn)行過(guò)程中，由于受到壅水作用的影響，在流量不變的情況下，隨著下游水深變化，水面流速系數(shù)和斷面流速分布均會(huì)發(fā)生變化。

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

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

表2 Δh不同時(shí)使用固定的k0產(chǎn)生的測(cè)流誤差Tab.2 Discharge measurement error using constant k0 under different values of Δh

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

3.3 水面流速系數(shù)

3.3.1 水面流速系數(shù)變化規(guī)律

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

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

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

圖9 各個(gè)斷面上水面流速系數(shù)與水深之間的關(guān)系Fig.9 The relation between water depth and surface velocity coefficient at different cross-sections

因此在渠道實(shí)際運(yùn)行過(guò)程中，Δh在20 cm 范圍內(nèi)時(shí)，用正常水深情況下率定的水面流速系數(shù)測(cè)流比較可靠，Δh大于20 cm 時(shí)，實(shí)際水面流速系數(shù)大于率定的水面流速系數(shù)，測(cè)量的流量會(huì)小于實(shí)測(cè)流量，造成一定的流量偏差。

3.3.2 考慮水深變化對(duì)水面流速系數(shù)進(jìn)行修正

本文選取三次多項(xiàng)式對(duì)距上游邊界50 m 斷面處的水面流速系數(shù)和下游邊界處的水深關(guān)系進(jìn)行擬合，得到修正后水面流速系數(shù)k'與下游水深hd的表達(dá)式：

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

由表3 可見(jiàn)，在本文所采用的工況下，k'值所推算出來(lái)的斷面流量與實(shí)際流量誤差小于1.4%。

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

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

3.3.3 多點(diǎn)流速測(cè)量

雷達(dá)測(cè)流采用的是無(wú)線電多普勒效應(yīng)。得到的水面數(shù)據(jù)的個(gè)數(shù)與水面波浪和漂浮物有關(guān)。因此灌區(qū)實(shí)際測(cè)流并不能像數(shù)值模擬一樣得到水面很多點(diǎn)的數(shù)據(jù)，只能測(cè)出水面有限點(diǎn)的流速。

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

表4 多點(diǎn)流速測(cè)流產(chǎn)生的誤差Tab.4 Discharge measurement error using multiple point measurement

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

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

4 結(jié)論

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

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

（2）當(dāng)所取渠段下游水深邊界變化Δh在20 cm 范圍內(nèi)時(shí)，用點(diǎn)流速和斷面平均流速的關(guān)系率定好的水面流速系數(shù)測(cè)流比較可靠，Δh大于20 cm 時(shí)，實(shí)際水面流速系數(shù)會(huì)大于率定的水面流速系數(shù)，測(cè)量的流量會(huì)小于實(shí)測(cè)流量，造成一定的流量偏差。

（3）考慮壅水導(dǎo)致的水深變化而對(duì)水面流速系數(shù)k修正時(shí)，測(cè)流誤差會(huì)比用固定水面流速系數(shù)明顯減小，最大測(cè)流誤差可從16.7%減小到1.4%。

（4）實(shí)際用雷達(dá)技術(shù)測(cè)表面流速時(shí)，設(shè)置的水面流速測(cè)點(diǎn)越多，流量測(cè)量結(jié)果越準(zhǔn)確。若可測(cè)水面7個(gè)點(diǎn)的流速時(shí)，測(cè)流誤差在9%以內(nèi)，在渠道量水規(guī)范［23］推薦的精度可接受范圍內(nèi)。

本文只選取了一個(gè)典型渠段且渠段長(zhǎng)度僅有100 m，該結(jié)論是否具有尺度效應(yīng)還需選取不同渠道特性的更長(zhǎng)的渠段進(jìn)行研究，在此基礎(chǔ)上還可對(duì)水面流速系數(shù)與糙率、斷面尺寸比、底坡等的關(guān)系進(jìn)行進(jìn)一步的研究?！?/p>