非恒定量法在三門峽庫區(qū)水沙輸移中的應(yīng)用

黃河三門峽水庫位于河南省石門市，控制面積68.84萬。km2，占總湖區(qū)面積的91.5%。從1919年到1985年，交叉口站多年的平均降水量和沙量分別為464億株和沙量。t3和15.6g，占總河流總數(shù)的89%和98%。由此可知,三門峽水庫控制黃河洪水和泥沙的作用很大,研究三門峽庫區(qū)的泥沙問題具有重要的意義。三門峽水庫自1960年9月建成并投入運(yùn)用40年以來,改變了水庫上、下游天然河道的水沙輸移狀況,從而引起了水庫上、下游河道沖淤所導(dǎo)致的河床調(diào)整。對(duì)于黃河河道沖淤計(jì)算的一維泥沙數(shù)學(xué)模型經(jīng)過國家“七五”、“八五”重點(diǎn)攻關(guān)項(xiàng)目的資助已經(jīng)取得了豐碩的成果。僅就水動(dòng)力學(xué)泥沙數(shù)學(xué)模型而言,在參考文獻(xiàn)中就包括了清華大學(xué)王士強(qiáng)的不平衡輸沙模型、黃河水利科學(xué)研究院曲少軍和張啟衛(wèi)的水庫泥沙數(shù)學(xué)模型、下游河道泥沙數(shù)學(xué)模型以及武漢水利電力大學(xué)韋直林的分組挾沙計(jì)算模型等等。這些數(shù)學(xué)模型在黃河相關(guān)河段的計(jì)算預(yù)測工作中已得到了較為滿意的結(jié)果。但在以往研究中,由于受計(jì)算機(jī)發(fā)展水平和黃河問題本身復(fù)雜性的限制,大部分的泥沙輸移數(shù)學(xué)模型還是以恒定水流和恒定輸沙為主,往往將非恒定過程概化為恒定流梯級(jí),即忽略了水流運(yùn)動(dòng)方程中?A/?t、?Q/?t和泥沙運(yùn)動(dòng)方程中?S/?t項(xiàng)。本文建立的一維泥沙輸移模型希望考慮水沙輸移過程中非恒定項(xiàng)的影響,一方面該模型能夠回答長時(shí)段河床變形的規(guī)律,另一方面也能回答在整個(gè)洪水過程中泥沙輸移和河床變形規(guī)律的變化,為日益受關(guān)注的數(shù)字化河流工程中實(shí)時(shí)調(diào)度系統(tǒng)的建立打下基礎(chǔ)。在計(jì)算結(jié)果中還比較了使用恒定流輸沙和非恒定流輸沙兩種不同的模型對(duì)計(jì)算結(jié)果的影響。1基本公式(1)t+xq+gahx+sf-s0b連續(xù)方程?Q?x+?A?t=q(1)動(dòng)量方程?Q?t+??x(Q2A)+gA(?h?x+Sf-S0)=0(2)式中Q為流量;h為水深;A為過水面積;q為匯入或分出流量;Sf為摩阻坡降;S0為河床的底坡降。(2)非均勻沙組成的含沙量變化泥沙非恒定輸移方程?(QS)?x+?(AS)?t+αωB(S*-S)=0(3)河床變形方程?Ζ0?t+αωB(S*-S)=0(4)式中S為斷面平均含沙量;S*為斷面平均挾沙力;Z0為河床高程;B為河寬;ω為泥沙沉降速度;α為泥沙恢復(fù)飽和系數(shù)。對(duì)于水流運(yùn)動(dòng)方程,本文采用四點(diǎn)偏心格式和追趕法進(jìn)行離散和求解。而對(duì)于泥沙運(yùn)動(dòng)方程,在以往研究中,求解方程式(3)時(shí),一般不考慮時(shí)變項(xiàng)?(AS)/?t的影響,而將它簡化為恒定流輸沙,應(yīng)用十分普遍的含沙量計(jì)算公式為S=S*+(S0-S0*)e-αωLQb+(S0-S0*)QbαωL(1-e-αωLQb)(5)式中Qb=Q/B為單寬流量;L為河段長度;下標(biāo)0為上個(gè)河段變量。如果不忽略時(shí)變項(xiàng),則推導(dǎo)非恒定輸沙式(3)可得到:?S?t+QA?S?x+(q-αωB)SA+αωBS*A=0(6)用特征線法解方程式(6)可得:S=Aq-αωBef(x-QAt)q-αωBAt-αωBq-αωBS*(7)由t=0時(shí)刻的初始條件解出:[S]tx=[S*]tx+[S-S*]0x-QAteαωB-qAt+[qq-αωBS*]0x-QAteαωB-qAt-[qq-αωBS*]tx(8)從式(8)可看出,出口斷面的含沙量[S]tx決定于出口斷面的挾沙力[S*]tx、初始斷面的含沙量和挾沙力的差值[S-S*]0x-QAt以及匯入流量q。假設(shè)式中出口斷面的挾沙力[S*]tx與進(jìn)口斷面的挾沙力[S*]0x-QAt相等,則有[S]tx=[S*]tx+[S-S*]0x-QAteαωB-qAt+[1αωB-q(?Q?x+?A?t)S*]tx(1-eαωB-qAt)(9)對(duì)?Q?x、?A?t使用離散格式后得到:[S]tx=[S*]tx+[S-S*]0x-QAteαωB-qAt+[1αωB-q(QΔx+AΔt)S*]tx(1-eαωB-qAt)-[1αωB-q(QΔx)S*]tx-Δx(1-eαωB-qAt)-[1αωB-q(AΔt)S*]t-Δtx(1-eαωB-qAt)(10)式(10)僅適用于均勻沙,即流速ω為常數(shù),將該式推廣應(yīng)用于非均勻沙時(shí),應(yīng)求分組含沙量的沿程變化,即設(shè)將非均勻沙分組,每組泥沙的含沙量計(jì)算均可利用式(10)。根據(jù)各組泥沙總量百分比pi的含義,應(yīng)有Si=piSS*i=piS*(i=1?2?3???n)(i=1?2?3???n)將式(10)對(duì)i求和,即可得到非均勻沙總含量的沿程變化:[S]tx=[S*]tx+[S-S*]0x-QAt[n∑i=1pieαωiB-qAt]0x-QAt-[AS*Δtn∑i=1pi1-eαωiB-qAtαωiB-q]t-Δtx+[(QΔx+AΔt)S*n∑i=1pi1-eαωiB-qAtαωiB-q]tx-[QS*Δxn∑i=1pi1-eαωiB-qAtαωiB-q]tx-Δx(11)式中上標(biāo)t和0分別為該時(shí)段與初始時(shí)段的計(jì)算值,下標(biāo)x為不同位置的計(jì)算值。對(duì)挾沙力的計(jì)算可采用張瑞瑾挾沙力公式和韓其為方法:S*=Κ[U3gR(n∑i=1piωi)10.92]m(12)式中U為斷面平均速度;K和m為經(jīng)驗(yàn)系數(shù)。2泥沙模型及網(wǎng)格黃河是多泥沙河流,由于淤積或沖刷的影響,黃河的河床斷面形狀每年都在發(fā)生變化,本文選用黃河1965年初的河床斷面數(shù)據(jù),計(jì)算時(shí)段從1965年初到1966年底兩年共730d。原始資料中河床的斷面形狀按近似階梯形處理,從潼關(guān)到三門峽全程113.7km的河道共有23個(gè)已知斷面,其中包括潼關(guān)、史家灘等水文站,每個(gè)水文站有水位及含沙量等驗(yàn)證資料。研究模型中使用網(wǎng)格的距離步長是不等間距的,泥沙計(jì)算時(shí)間步長選取Δt=3800s。進(jìn)口邊界選擇潼關(guān)以上,水流計(jì)算的邊界為流量隨時(shí)間的變化過程,泥沙計(jì)算的邊界為含沙量隨時(shí)間的變化過程,出口邊界選為三門峽壩址,根據(jù)水位調(diào)度作為下游邊界。3比較模型和計(jì)算結(jié)果(1)控制水力學(xué)的恒流模型模型驗(yàn)證分水力學(xué)和泥沙兩種模型進(jìn)行驗(yàn)證。水力學(xué)模型驗(yàn)證選取了代表性潼關(guān)水位和史家灘水位作為驗(yàn)證的數(shù)據(jù),驗(yàn)證結(jié)果見表1。結(jié)果表明水力學(xué)的恒定流計(jì)算模型和非恒定流計(jì)算模型都能取得了較為滿意的結(jié)果。泥沙模型的驗(yàn)證選取了整個(gè)河段的沖淤總量作為驗(yàn)證的數(shù)據(jù),在模型中分別采用了包含?S/?t項(xiàng)的非恒定流計(jì)算方法以及去掉?S/?t項(xiàng)的恒定流計(jì)算方法,從表2所示的對(duì)比結(jié)果可以看出,非恒定流模型計(jì)算的驗(yàn)證結(jié)果較恒定流模型計(jì)算的驗(yàn)證結(jié)果在精度上有了較大的提高。(2)降低了對(duì)于黃河沙量的側(cè)積和沙量與河底高程的沖刷經(jīng)本模型計(jì)算,選取其中2個(gè)典型斷面,分別繪出它們的流量、含沙量過程線如圖1、圖2所示。圖1為1965-1966年史家灘水文站的計(jì)算流量過程。在每年的6月到10月是洪水季節(jié),相應(yīng)的史家灘水文站流量增加,其中計(jì)算得到1965年最大洪峰流量為3984m3/s,發(fā)生在1965年的7月下旬,與實(shí)測流量3540m3/s能較好地符合。計(jì)算得到的1966年兩次洪峰流量分別為5942m3/s和5347m3/s,與實(shí)測洪峰流量4810m3/s和4870m3/s也能較好地符合。圖2為1965-1966年北村水文站的沙量隨時(shí)間變化的過程,由圖2計(jì)算過程可見,黃河高含沙河流的水土流失主要集中在汛期。1965年6月中旬最大含沙量約225kg/m3,1966年7月初和7月底兩次最大含沙量分別高達(dá)764kg/m3和523kg/m3,與發(fā)生的洪水過程相對(duì)應(yīng)。根據(jù)沙量平衡的關(guān)系,本模型還計(jì)算了在1965年初到1966年底長時(shí)段內(nèi),三門峽庫區(qū)河底高程的變化,如圖3所示。1965年1月到10月各斷面的河底高程降低較為明顯,表明在該時(shí)段內(nèi)黃河沖刷較為嚴(yán)重;從1965年10月到1966年6月,黃河沖刷較小,但還是有一些沖刷;從1966年6月到10月該河段基本沒有沖刷,相對(duì)還有少量淤積。這與表2實(shí)測的沖淤量相比大致相符。4非測定泥沙輸移曲線本文針對(duì)黃河流域的水沙輸移規(guī)律,應(yīng)用圣維南方程組以及非恒定泥沙連續(xù)方程建立了黃河三門峽庫區(qū)一維非恒定非均勻泥沙輸移數(shù)學(xué)模型,主要結(jié)論如下:(1)對(duì)水力學(xué)計(jì)算運(yùn)用圣維南方程和Preissmann四點(diǎn)偏心差分格式建立了一維非恒定流水力學(xué)模型;對(duì)非恒定泥沙輸移計(jì)算,應(yīng)用特征線方法建立了含沙量濃度的計(jì)算公式。并用黃河流域三門峽庫區(qū)1965年初到1966年底的實(shí)測資料對(duì)模型進(jìn)行了驗(yàn)證,計(jì)算得到的結(jié)果與實(shí)測資料符合良好,驗(yàn)證結(jié)果表明,此模型在非恒