山區(qū)多沙水庫干支流泥沙淤積數(shù)學(xué)模型及應(yīng)用

陳,舒彩文,武 見,梅志宏,談廣鳴

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,江蘇南京 210098;2.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430072;3.黃河勘測規(guī)劃設(shè)計(jì)有限公司,河南鄭州 450003;4.中國水電顧問集團(tuán)昆明勘測設(shè)計(jì)研究院,云南昆明 650051)

我國西部地區(qū)水能資源豐富,大力開發(fā)水電資源實(shí)現(xiàn)西電東送,是我國西部大開發(fā)戰(zhàn)略的重要組成部分.西部地區(qū)部分河流水土流失較為嚴(yán)重,含沙量較高,以元江為例,干流多年平均含沙量為3.5～4.3kg/m3,部分河段最大含沙量超過15kg/m3,在元江等多沙河道修建水庫后泥沙問題將較為突出,因而預(yù)測建庫后庫區(qū)泥沙淤積對于工程的優(yōu)化設(shè)計(jì)與運(yùn)行具有重要的意義.目前預(yù)測水庫泥沙淤積常采用一維水沙數(shù)學(xué)模型[1-5],國外比較有影響的一維模型有HEC-6模型[6]、FLUVIAL-12模型[7]、MIKE11模型[8]等,這些模型常用于低含沙沖積河流的沖淤計(jì)算,對于存在急流的山區(qū)多沙河流的沖淤計(jì)算不太適用.本文針對山區(qū)多沙河道水沙運(yùn)動的特點(diǎn),建立了適用于該類河流的一維非均勻非飽和輸沙數(shù)學(xué)模型,將模型應(yīng)用于元江上游戛灑江水電站庫區(qū)泥沙淤積計(jì)算,計(jì)算結(jié)果可為水庫設(shè)計(jì)與運(yùn)行提供依據(jù).

1 模型建立

1.1 基本方程

一維非均勻非飽和輸沙數(shù)學(xué)模型的基本方程如下:

水流連續(xù)方程

水流運(yùn)動方程

泥沙連續(xù)方程

河床變形方程

式中:x——沿程距離;Q——流量;Z——水位;g——重力加速度;B——河寬;t——時間;ql,ul,Sl——側(cè)向入流的單寬流量、流速和含沙量;A——過水?dāng)嗝婷娣e;R——斷面水力半徑;β——動能修正系數(shù);n——糙率系數(shù);S——斷面平均含沙量;S*——水流挾沙力;u——斷面平均流速;ω——泥沙顆粒沉速;ρ——泥沙干密度;α——恢復(fù)飽和系數(shù);Gb——推移質(zhì)輸沙率.

1.2 方程離散及求解

求解一維非恒定水流基本方程組采用有限體積法,此方法的特點(diǎn)是水量保持守恒.模型針對山區(qū)河道存在急流的特點(diǎn),對運(yùn)動方程的對流項(xiàng)采用一階迎風(fēng)格式進(jìn)行離散,由于該格式求解非線性對流擴(kuò)散方程時能獲得相當(dāng)高的精度[9],能夠較好地模擬山區(qū)河道為急流時的水流運(yùn)動,離散后的方程為

水流連續(xù)方程

式中上角標(biāo)為“0”的量表示上時刻的值.

1.3 補(bǔ)充方程

1.3.1 水流挾沙力公式

采用張瑞瑾水流挾沙力公式[10]:

分組水流挾沙力S*k采用竇國仁模式[11]:

其中

式中:S*——水流挾沙力;S*k——第k組沙水流挾沙力;K,M——水流挾沙力系數(shù)和指數(shù);U——斷面平均流速;——泥沙平均沉速;P*k,Pk——第k組沙水流挾沙力級配和懸沙級配;ωk——第k組沙對應(yīng)的沉速;N——非均勻沙粒徑組數(shù);β——指數(shù).

1.3.2 推移質(zhì)輸沙率

單寬推移質(zhì)輸沙率采用輸沙經(jīng)驗(yàn)關(guān)系式[12]:

其中

式中:vd——近底流速;d,d50——單組床沙粒徑和非均勻沙中值粒徑;h——平均水深.

1.4 支流入?yún)R處理

考慮到支流入?yún)R對干流的影響,方程(1),(2)和(3)均在等式右側(cè)增加了側(cè)向入流項(xiàng)的處理,而對于支流入?yún)R后干流含沙量的計(jì)算,采用以下公式:

式中:Qm,Sm,Qc,Sc——干、支流入?yún)R前上游最近斷面的流量和含沙量;Q0,S0——入?yún)R后首個干流斷面的流量和含沙量.

2 模型應(yīng)用

2.1 戛灑江水電站簡介

戛灑江一級水電站位于元江干流上游戛灑江河段,壩址流域面積為1.81萬km2,庫區(qū)由上游石羊江段和支流綠汁江組成,干流石羊江庫區(qū)平均比降約0.298%,支流綠汁江庫區(qū)平均比降約0.3%,庫區(qū)河道平面圖見圖1.大壩為面板堆石壩,最大壩高135 m,裝機(jī)容量200MW.水庫總庫容14.79億m3,有效庫容為6.72億m3,正常蓄水位為675m,死水位為660m.壩址處多年平均流量135m3/s,多年平均入庫沙量為2840萬t,多年平均含沙量為6.06kg/m3,歷史最大含沙量超過18.0kg/m3.戛灑江具有徑流小、含沙量高的特點(diǎn).本文應(yīng)用所建多沙水庫干支流泥沙淤積數(shù)學(xué)模型,對戛灑江水電站的泥沙淤積問題進(jìn)行研究.

2.2 計(jì)算范圍及邊界條件

2.2.1 計(jì)算范圍與斷面布置

計(jì)算范圍包括干流石羊江長約60km河段(進(jìn)口斷面S1到壩址處),以及支流綠汁江長約50km河段(進(jìn)口斷面L1到支流入?yún)R點(diǎn)),圖1為水電站庫區(qū)平面示意圖.其中石羊江共布置47個斷面,平均間距1.3km,支流綠汁江共布置41個斷面,平均間距1.2km.

2.2.2 邊界條件

模型進(jìn)口邊界條件為水庫干支流進(jìn)口斷面的流量和含沙量值,出口邊界為壩址處水位.由于模型進(jìn)口斷面處無水文站,其水沙過程通過壩址處元江水文站提供的水沙資料分配推算.其中,進(jìn)口采用1988—1997年共10a的水文泥沙資料作為計(jì)算系列,年平均流量為121.94m3/s,年平均含沙量為7.09kg/m3,屬于小水大沙的年份系列,入庫懸沙中值粒徑0.028mm.出口壩址正常蓄水位為675m,6—8月汛限水位為655m,9月汛限水位為665m.

模擬水庫100a泥沙淤積過程時,選取1988—1997年的水沙系列為邊界條件,前后連續(xù)計(jì)算10次共1200個月,模型計(jì)算時間步長為60s.

圖1 戛灑江一級水電站庫區(qū)示意圖Fig.1 Reservoir area of Jiashajiang Hydropower Station

2.3 有關(guān)計(jì)算參數(shù)的選取

由于缺乏實(shí)測資料對模型挾沙力系數(shù)率定,采用平衡輸沙的方法,假定河道建庫前多年平均含沙量和水流挾沙力基本相等來確定系數(shù).率定結(jié)果為:石羊江K平均值為0.93;綠汁江K平均值為1.16,M平均值為0.92.干支流糙率值通過實(shí)測水面線資料率定,糙率范圍為0.035～0.060.恢復(fù)飽和系數(shù) α在淤積時取為0.25,沖刷時取為1.0[13].

2.4 模型驗(yàn)證

圖2為2006年6月25日至9月30日石羊江干流25號斷面(壩上游33km)實(shí)測水位和模型計(jì)算值對比圖,從圖中可見兩者最大誤差不超過10cm,滿足工程要求.由于缺乏水庫泥沙沖淤驗(yàn)證資料,采用與戛灑江水電站相似的陡灘口水庫泥沙淤積資料進(jìn)行驗(yàn)證.圖3為陡灘口水庫淤積平衡后的沿程計(jì)算沖淤量與模型試驗(yàn)?zāi)M沖淤量對比圖,兩者誤差較小.所建一維非均勻非飽和輸沙數(shù)學(xué)模型計(jì)算水庫泥沙沖淤具有較好的精度.

圖2 水位實(shí)測與計(jì)算結(jié)果對比Fig.2 Comparison between observed and simulated water levels

圖3 陡灘口水庫沿程沖淤量驗(yàn)證Fig.3 Verification of deposition amount in Doutankou reservoir

2.5 水庫淤積計(jì)算成果分析

2.5.1 河床縱橫剖面變化

水電站運(yùn)行100a期間,庫區(qū)河道干流和支流河床縱剖面變化過程見圖4.由圖4可知,庫區(qū)干支流泥沙均呈現(xiàn)三角洲淤積形態(tài),三角洲淤積頂點(diǎn)隨著水庫運(yùn)行年限的增加逐漸向壩前推移,但淤積速度逐漸變緩.由于干流庫區(qū)泥沙淤積速度大于支流庫區(qū),當(dāng)水庫運(yùn)行40 a后,干流庫區(qū)泥沙淤積到支流入?yún)R口附近,干流部分泥沙進(jìn)入支流綠汁江庫區(qū)發(fā)生淤積,使得綠汁江淤積加快,支流庫區(qū)在運(yùn)行80a后達(dá)到?jīng)_淤平衡.

2.5.2 水庫庫區(qū)泥沙淤積及排沙比分析

表1為戛灑江庫區(qū)干流、支流及整個庫區(qū)泥沙累積淤積總量及排沙比計(jì)算結(jié)果.由表1可見,在計(jì)算年限100a內(nèi),水庫一直處于淤積狀態(tài),干流庫區(qū)運(yùn)行100a后未達(dá)到?jīng)_淤平衡,總淤積量為13.75億t,排沙比達(dá)到85.78%,支流庫區(qū)運(yùn)行80a后達(dá)到?jīng)_淤平衡,總淤積量為8.8億t,排沙比達(dá)到100%.

2.5.3 庫容變化

水庫初始運(yùn)行時正常蓄水位675m,水庫總庫容、死庫容和有效庫容分別為14.79億m3、8.07億m3和6.72億m3.水庫運(yùn)行后泥沙淤積導(dǎo)致庫容損失,表2為戛灑江一級電站運(yùn)行后庫容變化情況及庫容損失率情況.從表2可以看出,水庫泥沙不斷發(fā)生淤積,庫容損失不斷增加.水庫運(yùn)行40a后,死庫容損失率接近50%,運(yùn)行至100a,死庫容損失率達(dá)到100%,有效庫容損失率為6.53%,總庫容損失率為57.53%.

圖4 庫區(qū)河道深泓線變化Fig.4 Evolution of talweg in reservoir area

表1 戛灑江庫區(qū)泥沙累積淤積總量及排沙比計(jì)算結(jié)果Table 1 Simulated results of total sediment deposition amount and delivery rate in Jiashajiang reservoir area

表2 建庫后戛灑江水庫庫容損失特征值Table 2 Values of loss rate of storage in Jiasajiang reservoir

3 結(jié) 論

針對我國西南部山區(qū)多沙河流水沙運(yùn)動的特點(diǎn),建立了一維非恒定非飽和輸沙數(shù)學(xué)模型用于水庫泥沙淤積預(yù)測,模型對運(yùn)動方程對流項(xiàng)采用一階迎風(fēng)格式進(jìn)行離散,使得模型能夠?qū)ι絽^(qū)河道急流流態(tài)進(jìn)行模擬.最后將模型用于元江上游戛灑江水電站庫區(qū)泥沙淤積計(jì)算,結(jié)果表明:干支流庫區(qū)均呈三角洲淤積形態(tài),三角洲淤積頂點(diǎn)隨運(yùn)行年限的增加逐漸向壩前推移,但淤積速度逐漸變緩.干流庫區(qū)運(yùn)行100a后未達(dá)到?jīng)_淤平衡,支流庫區(qū)運(yùn)行80a后達(dá)到?jīng)_淤平衡.干流庫區(qū)運(yùn)行100a后總淤積量為13.75億t,排沙比達(dá)到85.78%,支流庫區(qū)運(yùn)行80a后總淤積量為8.8億t,排沙比達(dá)到100%.水庫運(yùn)行100a后,干支流累計(jì)有效庫容損失率為6.53%,死庫容損失率為100%,總庫容損失率為57.53%.

[1]湯立群,陳界仁,陳國祥,等.白沙水庫一維泥沙數(shù)學(xué)模型及電廠取水方案計(jì)算[J].河海大學(xué)學(xué)報(bào),1998,26(5):34-38.(TANG Li-qun,CHEN Jie-ren,CHEN Guo-xiang,et al.One dimensional model of reservoir sedimentation and its application incalculation for water intake of a power plant[J].Journal of Hohai University,1998,26(5):34-38.(in Chinese))

[2]楊國錄,吳衛(wèi)民,陳振虹,等.SUSBED-1動床恒定非均勻全沙模型[J].水利學(xué)報(bào),1994(4):1-10.(YANG Guo-lu,WU Weiming,CHEN Zhen-hong,et al.SUSBED-1 movable bed modeling of graded sediments[J].Journal of Hydraulic Engineering,1994(4):1-10.(in Chinese))

[3]張小峰,謝葆玲,許全喜.三峽水庫蓄水初期泥沙淤積預(yù)測模型研究[C]//中國水力發(fā)電工程學(xué)會水文泥沙專業(yè)委員會第六屆學(xué)術(shù)討論會論文集.成都:時代出版社,2005.

[4]黃仁勇,黃悅.三峽水庫干支流河道一維非恒定水沙數(shù)學(xué)模型初步研究[J].長江科學(xué)院院報(bào),2009,26(2):9-13.(HUANG Ren-yong,HUANG Yue.Preliminary study on 1-D numerical simulation of unsteady flow and sediment transport in mainstream and tributaries of three gorges reservoir area[J].Journal of Yangtze River Scientific Research Institute,2009,26(2):9-13.(in Chinese))

[5]黃永健,談廣鳴,呂平,等.黏性泥沙數(shù)學(xué)模擬中幾個關(guān)鍵問題的研究進(jìn)展[J].水利水電科技進(jìn)展,2009,29(3):90-93.(HUANG Yong-jian,TAN Guang-ming,LV Ping,et al.Research advancesof severalkey problems in the numerical simulation of cohesive sediment[J].Advances in Science and Technology of Water Resources,2009,29(3):90-93.(in Chinese))

[6]FELDMAN A D.HEC models for water resources system simulation:theory and experience[J].Advances in Hydroscience,1981(12):297-423.

[7]CHANG H H,HANISON L L,LEE W,et al.Numerical modeling for sediment-pass-through reservoirs[J].Journal of Hydraulic Engineer,ASCE,1996,122(7):381-388.

[8]王領(lǐng)元.應(yīng)用MIKE對河流一、二維的數(shù)值模擬[D].大連:大連理工大學(xué),2007.

[9]張小峰,中川一,許全喜.一階迎風(fēng)差分格式的精度問題[J].武漢大學(xué)學(xué)報(bào):工學(xué)版,2001,34(1):6-10.(ZHANG Xiao-feng,NAKAGAWA H,XU Quan-xi.On accuracy of first order upwind scheme[J].Journal of Wuhan University of Hydraulic and Electric Engineering,2001,34(1):6-10.(in Chinese))

[10]張瑞瑾.河流泥沙動力學(xué)[M].北京:中國水利水電出版社,1989.

[11]竇國仁,趙士清,黃亦芬.河道二維全沙數(shù)學(xué)模型研究[J].水利水運(yùn)工程學(xué)報(bào),1987(2):1-12.(DOU Guo-ren,ZHAO Shi-qing,HUANG Yi-fen.Study on two-dimensional total sediment transport mathematic model[J].Hydro-Science and Engineering,1987(2):1-12.(in Chinese))

[12]國務(wù)院三峽工程建設(shè)委員會辦公室泥沙專家組.長江三峽泥沙問題研究:第五卷:1996—2000[M].北京:知識產(chǎn)權(quán)出版社,2002.

[13]韓其為.水庫淤積[M].北京:科學(xué)出版社,2003.