李林博，王于剛，龔秀秀

（1.陜西水環(huán)境設(shè)計(jì)集團(tuán)有限公司,陜西 西安 710018；2.中國市政工程西北設(shè)計(jì)研究院有限公司,甘肅 蘭州730000）

運(yùn)用數(shù)學(xué)模型對河道水動(dòng)力場及涉河建筑物的模擬具有廣泛的應(yīng)用。宋昊明等人[1]采用數(shù)學(xué)模型方法分析了無為縣白茆取水工程對河道行洪的影響；任杰等人[2]采用模型定量計(jì)算了防浪林對河道行洪的影響；梁艷慧等人[3]構(gòu)建平面二維數(shù)學(xué)模型,模擬了河道在變交汊角條件下分流比特性,開展了多因素影響分析；寇爾丹[4]利用數(shù)值計(jì)算方法對曲線疊水低堰進(jìn)行了模擬；郭鳳清等人[5]構(gòu)建洪水演進(jìn)模型對洪水進(jìn)行動(dòng)態(tài)量化的估算,為防洪調(diào)度和防洪規(guī)劃具有指導(dǎo)意義。楊填堰工程壩軸線與河道主流斜交,平面布置具有其特殊性,附近水動(dòng)力特性與以往工程有所不同。本文利用MIKE軟件中二維水動(dòng)力模塊構(gòu)建楊填堰所在河段平面二維數(shù)學(xué)模型,分析不同工況下該河段水動(dòng)力特性。

1 研究背景

楊填堰灌區(qū)始建于漢代,其工程建設(shè)歷史悠久,灌區(qū)工程體系經(jīng)歷次修繕,目前灌區(qū)灌溉面積共計(jì)2.13 萬畝。2017年,該堰與山河堰、五門堰三堰被授予世界灌溉工程遺產(chǎn),對當(dāng)?shù)剞r(nóng)業(yè)經(jīng)濟(jì)的發(fā)展及水利文化的弘揚(yáng)作用更加重要。

2021 年9月,湑水河流域普降大到暴雨,26 日楊填堰上游升仙村水文站形成最大洪峰2660 m3/s。本次洪水造成楊填堰堰體水毀110 m溢流壩、10 m擋水壩。經(jīng)調(diào)查,該堰歷史上有記錄的修復(fù)重建共計(jì)十余次,僅建國以來,楊填堰水毀發(fā)生共計(jì)六次。

表1 建國以來?xiàng)钐钛甙l(fā)生水毀情況統(tǒng)計(jì)

2 計(jì)算模型

2.1 研究區(qū)域概況

湑水河系漢江左岸一級支流,發(fā)源于陜西省周至縣境內(nèi)秦嶺光頭山,流域面積2307 km2,干流全長165.5 km,河道平均比降5.59%。本次研究區(qū)域位于楊填堰所在的湑水河下游平原河段。

楊填堰堰體走向：右岸向上游傾斜,與湑水河河道主流大致成45°,樞紐總長321 m,主要建筑物包括：溢流壩252 m,右岸擋水壩69 m。堰壩左岸設(shè)沖沙閘3孔,灌溉引水閘1孔。樞紐設(shè)計(jì)引水流量3.0 m3/s,最大引水流量4.0 m3/s。

圖2 楊填堰平面布置

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

MIKE21 FM水動(dòng)力模塊是基于數(shù)值解的二維淺水方程,可用于模擬各種作用力下產(chǎn)生的水位、水流變化及模擬二維自由表面流,其采用非結(jié)構(gòu)化網(wǎng)格,數(shù)值計(jì)算方法是單元中心的有限體積法。模型的控制方程是由質(zhì)量守恒和垂向積分的動(dòng)量方程組成。

質(zhì)量守恒方程：

X方向動(dòng)量方程：

Y方向動(dòng)量方程：

式中：ζ為水位,m；t為時(shí)間,s；x、y、表示平面坐標(biāo),m；d為隨時(shí)間變化的水深,m；h為水深（h=ζ-d,m）；p、q分別為單寬流量在x與y方向上的分量,m3/（s·m）,其中p=uh,q=νh；u、ν分別為x與y方向上的垂向平均流速,m/s；C為謝才系數(shù),m1/2/s；g為重力加速度,m/s2；f為風(fēng)摩擦系數(shù)；V、Vx和Vy分別為風(fēng)速及其在x與y方向上的分量,m/s；為科氏（Coriolis）系數(shù),s-1；Pa為大氣壓強(qiáng),Pa；ρw為水的密度,kg/m3；yy、xy、τyy為有效切應(yīng)力分量,N/m2。

2.3 基礎(chǔ)數(shù)據(jù)及模型建立

本次建模原始數(shù)據(jù)采用2018年實(shí)測的楊填堰上、下游各1 km的1∶1000河道帶狀地形圖及研究區(qū)域內(nèi)23條橫斷面數(shù)據(jù)。模型劃分采用非結(jié)構(gòu)化網(wǎng)格,為保證生成網(wǎng)格質(zhì)量,設(shè)定網(wǎng)格最大角度30°。本次共創(chuàng)建網(wǎng)格節(jié)點(diǎn)64688個(gè),網(wǎng)格128100 個(gè),其中楊填堰處加密網(wǎng)格最大面積0.5 m2,其余河段網(wǎng)格最大面積20 m2。對網(wǎng)格進(jìn)行平滑處理后內(nèi)插高程點(diǎn),形成地形文件,見圖3。

圖3 模型地形

2.4 參數(shù)選取

研究區(qū)域位于UTM-49（墨卡托投影）分區(qū),計(jì)算考慮科氏力影響。模型設(shè)置計(jì)算時(shí)間步長為0.01 s,共計(jì)1.08×106步,模擬時(shí)間3 h。模型上游入口設(shè)置為流量邊界,下游出口為水位-流量邊界,水位-流量關(guān)系曲線見圖4。為避免過強(qiáng)淺水效應(yīng)的影響,對過淺水區(qū)域進(jìn)行忽略處理,其中定義干水深0.05 m,淹沒水深0.1 m,濕水深0.15 m。模型計(jì)算初始水深設(shè)置參照計(jì)算的實(shí)際情況設(shè)置為常水位692.00 m。底床阻力設(shè)置為曼寧系數(shù),參照《水力計(jì)算手冊》同時(shí)結(jié)合《城固縣湑水河右岸呂家村段防洪工程初步設(shè)計(jì)》中相關(guān)成果,綜合確定曼寧數(shù)取值為28.5 m1/3/s。

圖4 出口邊界水位—流量關(guān)系曲線

2.5 計(jì)算工況

本文對研究區(qū)域在工況一湑水河多年平均流量（34.69 m3/s）、工況二“9·26”洪水洪峰流量（2660 m3/s）、工況三20年一遇洪水洪峰流量（3360 m3/s）共3 個(gè)不同工況進(jìn)行模擬計(jì)算。

3 參數(shù)率定與模型驗(yàn)證

數(shù)值計(jì)算模型參數(shù)率定參照2021年“9·26洪水”洪痕調(diào)查資料進(jìn)行。經(jīng)率定,模型綜合曼寧系數(shù)取值為28.5 m1/3/s較為合理。選取計(jì)算區(qū)域內(nèi)5 個(gè)洪痕調(diào)查點(diǎn)（見圖5）,數(shù)值計(jì)算結(jié)果與洪痕對比見表2。模擬結(jié)果與洪痕結(jié)果最大差值17 cm,模擬精度滿足本次計(jì)算要求,數(shù)值模型對楊填堰段河道水動(dòng)力特性模擬是可行的。

表1 測點(diǎn)數(shù)值計(jì)算水位與洪痕比較表 單位：m

圖5 調(diào)查洪痕位置及結(jié)果分析斷面

4 計(jì)算結(jié)果與分析

采用前述建立的模型對各工況進(jìn)行數(shù)值計(jì)算,本次分析楊填堰附近河段水動(dòng)力特征,選取局部計(jì)算成果進(jìn)行分析,模擬得到的水位分布及流速場見圖6～圖8。為分析楊填堰上下游水位、流速變化規(guī)律,在楊填堰溢流壩及擋水壩段分別布設(shè)特征斷面S1、S2,斷面分布見圖5。

圖6 工況一（34.69 m3/s）水位分布及流速場

圖7 工況二（2660 m3/s）水位分布及流速場

圖8 工況三（3360 m3/s）水位分布及流速場

4.1 水位變化情況

由圖6～圖8分析可知,湑水河多年平均流量時(shí),因流量較小,楊填堰僅溢流壩段過流,堰后水流繞行下游河心灘；工況二、三流量相對較大,溢流、擋水壩段均過流,下游全斷面過流。

圖9、圖10為以距離為橫坐標(biāo),水位為縱坐標(biāo)的楊填堰附近水位變化圖,圖中橫坐標(biāo)距離以堰體上游面為0向上為負(fù),向下為正。從圖中可以看出：（1）溢流堰在流量較小時(shí),堰面水流比較平順,隨著流量增大,堰面水流發(fā)生波動(dòng),在堰面斜坡上水位先減小后增大,考慮在堰面斜坡上發(fā)生水躍；（2）隨著流量增大,溢流堰堰體上游及其堰面水面比降逐漸加大,而在堰后平臺及河道水面比降變化不大,堰體上下游水位差逐漸減小,楊填堰壅水效果相對減弱,分析是由于楊填堰兩側(cè)擋水壩及灘面同時(shí)過流,過流斷面加大所致；（3）楊填堰擋水壩上下游水位沒有明顯波動(dòng),相對平穩(wěn)；上下游水位差隨流量增大而減小,與溢流堰段水位變化趨勢相同。

圖1 楊填堰所在位置

圖9 溢流壩S1斷面水位變化

圖10 擋水壩S2斷面水位變化

4.2 流速場變化情況

由圖6～圖8流速云圖分析可以看出：各種工況下,（1）堰體上下游流速因堰體本身走向的影響,其右岸流速較左岸流速偏大；（2）隨著流量增大,堰面最大流速在各工況下分別為1.35 m/s、8.4 m/s、9.6 m/s,呈現(xiàn)逐漸增大趨勢；（3）在擋水壩與溢流壩段交界處形成偏向左岸的楔形大流速區(qū),楔形區(qū)域內(nèi)最大流速區(qū)的長度、面積也隨流量增大而增大,楔形核心區(qū)流速為4.8 m/s的堰后長度從10 m增加到29 m。

5 結(jié)論

本文通過MIKE軟件構(gòu)建湑水河楊填堰段平面二維數(shù)學(xué)模型,主要分析了不同工況下楊填堰附近河段的水位、流速等水動(dòng)力特征分布規(guī)律,主要研究結(jié)論如下：

1）從水位分布變化情況來看,楊填堰溢流壩及擋水壩上下游水位差隨著流量增大而減小,楊填堰壅水效果相對減弱；隨著流量增大,將在堰面發(fā)生水躍。

2）從流速場變化情況可以看出,堰體上下游流速因堰體走向影響,其右岸流速較左岸流速大；堰面最大流速隨著流量增大而逐漸增大,20年一遇洪水時(shí)堰面最大流速為9.6 m/s。

3）因右岸擋水壩導(dǎo)流作用,在溢流壩與擋水壩段交界處形成楔形大流速區(qū),且該區(qū)域長度、面積隨流量增大而增大。

綜合以上分析,當(dāng)來流量較大時(shí),在楊填堰溢流壩與擋水壩段交界處形成楔形大流速區(qū),考慮將首先在此處堰體下游造成沖刷,形成沖坑從而造成堰體不穩(wěn)發(fā)生水毀,堰體多次從此處發(fā)生險(xiǎn)情也驗(yàn)證了本次數(shù)值分析成果。本文的計(jì)算成果可為楊填堰的水毀修復(fù)設(shè)計(jì)提供參考,對類似工程的設(shè)計(jì)具有借鑒意義。