楊亮亮
(廣晟昊興勘測設計有限公司甘肅分公司,甘肅 蘭州 730000)
黃河干流在甘肅境內全長913km,分別流經甘南、臨夏、蘭州、白銀4個市(州)的9個縣6個區。近年來,為完善白銀市城市防洪體系,避免已建防洪堤帶病運行,提升重要城市防洪減災能力,對黃河干流白銀段進行防洪治理。本次黃河干流白銀段防洪治理工程中一個關鍵問題是確定各洪水頻率工況下的河道水面線位置,為堤防和護岸工程的加高加固提供理論依據。采用能量方程逐段試算法(伯努利方程法)對白銀市段局部區域內的水面線進行計算,但由于能量方程逐段試算法是一種一維方法,且在適用條件和計算方法上存在諸多的局限,使計算結果在可靠性方面還存在一定的問題。為對計算結果進行進一步的校核,本次水面線復核計算采用二維水動力數值模型,對多種洪水頻率工況下黃河干流白銀市水川段的水面線高程進行模擬計算,并對能量方程逐段試算法計算結果進行復核和評價。
黃河干流上游白銀段洪水相對漲落較緩,洪水歷時和間隔較長,洪次較少,水面線推算時認為河段設計洪水下的水流流態是明渠恒定漸變流。在明渠恒定漸變流條件下,工程設計時選擇一維能量方程來推求河道不同斷面的水位[1]。天然河道水面線計算基本方程為
(1)

為克服能量逐段試算法在河道水面線計算上的局限和不足,本文對計算結果進行進一步的校核優化,采用二維水動力數值模型,對多種洪水頻率工況下黃河干流白銀市水川段水面線高程進行模擬計算,并與一維能量方程推求的河道水面線計算結果進行比較,對能量方程逐段試算法計算結果進行復核和評價。
采用二維淺水方程來描述典型河段的水流流動問題,并選擇Poisson方程作為坐標轉換方程,把物理域上的不規則區域轉化為計算域上矩形區域,因為不規則的邊界在離散時可能會出現顯著的誤差[2-3]。曲線坐標系下的二維淺水方程包括:
a.沿水深積分的連續方程:
(1)
式中:u、v分別為ξ、η方向上的水流速度,m/s;Gξξ、Gηη分別為曲線坐標系轉換成直角坐標的轉換系數;d為基準水位下的水深,m;ζ為基準水位(z=0)以上的水位,m。
b.動量方程:
(2)
(3)
式中:u、v、Gξξ、Gηη含義同上;fu、fv分別為柯氏力系數(地理緯度北緯33°),1/s;Fξ、Fη分別為ξ、η方向上的紊動動量通量,kg/s2。
c.定解條件:

邊界條件:出、入流邊界上,給定水位、流速或流量過程,固壁邊界采用無滑移邊界條件。

本次模擬計算的數值離散方法采用有限差分法,計算方法是基于交錯網格的交替方向顯、隱式混合格式,即ADI法。ADI方法將一個時間步長分為兩個時間層,各由半個時間步長組成。在每一時間層中,模型方程的所有項通過在空間上至少有二階精度的相同方法求解。ADI方法的優點是計算穩定性好、計算精度高,缺點是靈活性差,不適用于具有混合偏導數的情況。在二維模擬計算中,ADI方法是無條件穩定的[4]。

3.1.1 斷面布置
黃河干流在白銀市境內流經白銀區、靖遠縣、平川區、景泰縣,自西南的大峽水電站庫區流入,經白銀區,至景泰縣紅柳后拐向西北,經過景泰縣五佛后折向東北流出甘肅境內,轄區內河道彎曲,平面上大致呈“S”形,境內河長254km,主要為天然河道段。根據河道分段及斷面布置原則,共劃分244條橫斷面,平均斷面間距500m。本次計算范圍從橫1斷面至橫37斷面,河道長度約21km,計算范圍內包含38個實測橫斷面,上游邊界為大峽水電站;烏金峽水電站在計算區域下邊界下游約9km處,其尾水影響范圍至橫23斷面;安寧渡水文站位于計算范圍下邊界下游約70km處。圖2為黃河干流白銀市水川段計算區域示意圖。
圖2中藍色線表示本次防洪治理工程設計狀況計算范圍,與天然狀況下計算范圍不同的地方用黃色線標示。由圖2可知,天然狀況和設計狀況計算區域的主要差異位于橫11斷面至橫30斷面之間。
3.1.2 網格劃分
由于天然狀況和設計狀況計算區域不同,在計算中須生成2套網格分別對天然條件下和設計條件下的工況進行計算。在本次模擬計算中所采用的天然狀況和設計狀況下的網格點數分別為29×761和22×756。白銀市水川段天然狀況和設計狀況計算網格見圖3。

圖3 白銀市水川段設計狀況計算網格
3.1.3 地形插值
為在模擬計算中真實反映計算段河道的地理信息,采用地形插值方法對白銀市水川段計算區域內的高程散點進行插值計算生成地形文件[5]。插值計算所涉及的高程散點總數為1403370個,其中橫斷面散點數為2514個,等高線散點數為1393638個,地形高程散點數為7218個。白銀市水川段插值地形圖見圖4。

圖4 白銀市水川段天然狀況插值地形圖
白銀市水川段計算分10年一遇洪水、2012年實測最大洪水、5年一遇枯水期施工洪水和冬季實測流量4種工況,具體計算工況見表1。

表1 白銀市水川段計算工況
邊界條件:兩種方法在河道水面線計算時,工況4的下游水位邊界條件采用實測值。對于工況1~工況3,由于缺乏水位實測資料,采用能量方程計算值作為下游邊界條件。
由白銀市安寧渡水文站實測糙率-流量關系(見圖5)可知,糙率實測值在流量為1000m3/s上下時有明顯區別:當流量大于1000m3/s時,糙率值變化不大,介于0.025~0.032之間;當流量小于1000m3/s時,糙率值隨流量減小而顯著增大,介于0.030~0.060之間。因此,在本次模擬計算中,對工況1~工況3(Q>1000m3/s),采用同一套槽-灘-植被的率定糙率組合值,對工況4(Q=600m3/s),由于流量小,水流不上灘,因此在計算區域內不設邊灘、植被糙率,取安寧渡水文站實測糙率值0.045(可由圖5得到)。在對工況1~工況3的模擬計算中進行比選的糙率組合有:?河道主槽糙率0.03,邊灘糙率0.05,植被糙率0.06~0.15;?河道主槽糙率0.03,邊灘和植被糙率0.04;?河道主槽糙率0.03,邊灘和植被糙率0.06;?河道主槽糙率0.03,邊灘糙率0.05,植被糙率0.10。
采用上述糙率組合對白銀市水川段工況進行模擬計算,將其計算結果與能量方程計算值進行比較,選擇擬合結果最優的第一組糙率組合為工況1~工況3模擬計算的糙率組合值,糙率值的具體分布見圖5(天然狀況下)。

圖5 安寧渡水文站實測糙率-流量關系散點圖
2012年12月1—3日對白銀市水川段進行了水位流量現場觀測(工況4)。利用該實測資料對模型進行率定,由于冬季流量較小,水流不上灘,因此在整個計算范圍內采用單一糙率值0.045。白銀市水川段沿程各斷面上的實測值、二維計算結果和兩者之間的差值見表2。工況4的二維計算值與實測值的對比見圖6。由于橫23斷面后的區域內水位受下游烏金峽水電站尾水影響,其尾水影響無法在二維數值模型中進行模擬,因此本文僅比較分析橫1斷面至橫23斷面間區域的結果。由表2可知,二維計算結果與實測值之間的最大差值為1.50m,出現在橫19斷面。
由圖6可知,二維計算值與實測值之間差值主要出現在橫3斷面至橫8斷面、橫15斷面至橫17斷面以及橫19斷面至橫22斷面。具體原因分析如下:
a.由圖6可知,二維計算值在橫3斷面至橫8斷面間低于實測值,其主要原因是二維計算值在橫2斷面和橫3斷面間存在突降現象。由橫2斷面至橫4斷面的地形圖和實地踏勘可知:橫2斷面至橫3斷面之間距離較大(約1.4km),遠大于水川段的平均橫斷面距離(500m),水流在橫2斷面以后先經過一個近90°的彎道,且伴隨有明顯的河床縮窄現象。天然工況1條件下橫1斷面至橫3斷面流速分布見圖7,由圖7可知,水流在出彎以后呈現向對岸折沖的趨勢,主流整體偏向右岸,左岸局部出現小范圍回流,經過一段距離的調整,主流在橫3斷面處回到河道中心。理論上,以上地形與河勢條件會引起河道內水流的能量損耗,使橫2斷面上游局部出現壅水,橫2斷面和橫3斷面間出現較明顯的水位跌落現象。本次模擬計算中進行了多組針對糙率與地形的數值試驗,結果表明橫2斷面與橫3斷面間較大的水位跌落不可避免。由于橫2斷面至橫3斷面間較大的水位跌落現象,造成了平面二維計算水位值在橫3斷面至橫8斷面間整體低于實測值。
b.二維計算值在橫15斷面至橫17斷面間低于實測值,其原因也與橫14斷面至橫17斷面之間的地形和河勢有關。該河段為45°彎道,且在橫15斷面存在河心洲,受河心洲的擠壓,橫15斷面附近水面寬度僅為橫14斷面水面寬度的一半。受地形和河勢的影響,二維計算值在橫14斷面和橫15斷面間出現跌落現象。同樣,在橫17斷面下游約50m處,河道中再次出現阻水洲灘,使過水斷面再次縮窄,水面寬度約為橫17斷面處水面寬度的3/5左右,且河床底高程出現負坡,造成橫16斷面至橫17斷面水面坡降放緩。
c.二維模擬值在橫19斷面至橫22斷面間高于實測值。結合實測地形可知:一方面,橫18斷面至橫20斷面間,以及橫25斷面至橫27斷面間河道中央各有一個洲,在橫20斷面和橫24斷面,河道有兩處明顯束窄;另一方面,橫18斷面至橫19斷面河床為負坡,橫19斷面至橫24斷面,河床坡度較緩。綜合兩方面因素,造成橫19斷面至橫22斷面的壅水現象,水面坡降放緩。
由于橫23斷面后的區域內水位受下游烏金峽水電站尾水影響,其尾水影響無法在二維數值模型中進行模擬,本文僅分析各工況下橫1斷面至橫23斷面間區域的計算結果。篇幅所限,10年一遇洪水(天然狀況計算結果)見表3,列出兩個工況(工況2和工況3)下能量方程逐段試算法與二維數值計算值進行比較分析,比較結果見表4,二維計算值與能量方程計算值的對比見圖8。

表4 二維計算水位高程與能量方程計算值比較成果(工況2和工況3) 單位:m

圖8 白銀市水川段水面線沿程分布
由兩種計算方法比較可知:在下游邊界取能量方程逐段試算法計算值的情況下,對于工況1的天然狀況和設計狀況,二維計算值與能量方程計算值的最大差值分別為1.77m和1.44m,均出現在橫16斷面;對于工況2、工況3和工況4,二維計算值與能量方程計算值的最大差值分別為1.68m、1.28m和1.20m,分別出現在橫20斷面、橫19斷面和橫19斷面。
由圖8可知:
a.二維計算值與能量方程計算值的最大差值為1.77m,發生在天然工況1的橫16斷面。
b.在大流量工況(工況1和工況2)下,橫6斷面至橫10斷面間的二維計算值大于能量方程計算值,且在該河段以較大比降下降;隨著流量降低,橫6斷面至橫10斷面間二維計算值開始減小,部分或全部低于能量方程計算值,但仍維持較大比降。這是因為橫6斷面至橫10斷面間河道有兩個大角度的轉角,且伴隨有河床束窄的現象。隨著河道中水位降低,地形的影響開始占據主導地位。
c.工況1的天然狀況和設計狀況下,二維計算值在橫13斷面至橫23斷面間均低于能量方程計算值,這是由于橫13斷面至橫16斷面間的水面線較大幅度的降落引起的。橫12斷面至橫16斷面間為一個近120°的彎道,彎道內側為顧家善,是一塊灘地。在10年一遇洪水工況下,水流上灘,使橫12斷面至橫16斷面間的糙率增加,再加上彎道的作用,使橫13斷面至橫16斷面間的水面線產生突降。一維模型中采用綜合糙率,難以考慮水流上灘所引起的糙率增大現象和彎道的影響,可能致使二維計算值在橫13斷面至橫23斷面間低于能量方程計算值。隨著流量的降低,工況2~工況4中顧家善區域內均無水流上灘現象,因此在二維計算中沒有造成糙率的突然增大,因此對比能量方程計算值,沒有出現水面線的突降。隨著河道中水位的降低,地形的影響開始占據主導地位,橫19斷面至橫22斷面間的壅水現象開始控制橫18斷面至橫23斷面的水位值,使這段河道中的二維計算值大于能量方程計算值。
利用能量方程和二維水動力數值模型,在10年一遇洪水、2012年實測最大洪水、5年一遇枯水期施工洪水和冬季實測流量4種工況下對黃河干流白銀市水川段水面線高程進行計算,并對兩種計算結果進行比較分析,得到以下結論:
a.一維能量方程是計算河道水面線的基本理論與方法,是一種適用于恒定流的方法,適用于河道順直、河道平面形態變化不大且兩推算斷面間過水斷面相對均勻的河段,能量方程逐段試算法基本滿足防洪工程一般河道水面線的需要,是一種適用廣泛、簡便易行的工程計算分析方法。
b.二維淺水方程適用于研究海岸、河道入海口、湖泊、大型水庫等具有廣闊水域的地區,對河道彎道、斷面突變、河道水位突降、壅水以及糙率系數變化的河段計算結果準確,當用于河流水面線計算時,適合于研究具有大片灘地的河段。在洪水到來時,如整個灘地上的水流不能假定是沿同一方向的,則應建立二維模型,充分地重演洪水的側向演進。
總之,工程設計中可將一維、二維模型嵌套計算,建議對于相對順直,有堤防約束的河段,可用能量方程計算河道水面線;對于局部地形復雜,水流在灘、槽之間游移變動的河段,以及存在大片灘地的天然狀況或者無堤防約束的洪水漫流等,可利用二維淺水方程進行建模計算。同時兩種計算方法在河道沿線堤防進行提前預報、預警方面有重要的意義,針對預警結果采取相應的應急措施,能夠快速準確地進行防洪。