王 玨,張海榮,孫振宇,張明亮,3

(1.大連海洋大學(xué)海洋科技與環(huán)境學(xué)院,遼寧 大連 116023;2.中國(guó)長(zhǎng)江電力股份有限公司水資源研究中心,湖北 宜昌 443133;3.大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧 大連 116025)

洪水是由暴雨、急劇融冰化雪、風(fēng)暴潮等引起的江河湖泊水量激增、水位上漲的一種自然災(zāi)害。通過對(duì)歷史資料分析發(fā)現(xiàn):擴(kuò)大耕地、圍湖造田等人類活動(dòng)不斷破壞地表下墊面形態(tài),致使洪水頻率升高,致災(zāi)程度也不斷加深。在強(qiáng)降雨期間,短時(shí)間內(nèi)流域水量突增、匯聚,形成洪水在河網(wǎng)水系內(nèi)傳播。針對(duì)河網(wǎng)內(nèi)突發(fā)的洪水問題,國(guó)內(nèi)外諸多學(xué)者開展了大量的研究工作。早期學(xué)者多采用瞬時(shí)流態(tài)法[1]、馬斯京根法[2-3],由于上述方法過度簡(jiǎn)化圣維南方程,在精度上往往不能滿足實(shí)際需求?？紤]到河網(wǎng)水動(dòng)力模型的核心是對(duì)圣維南方程組的求解,Preissmann和Cunge[4]提出四點(diǎn)隱式差分格式,經(jīng)過不斷完善,因其收斂快、效率高并且穩(wěn)定性好,被廣泛應(yīng)用于實(shí)際工程中[5-8]。一般來說,河網(wǎng)水系常表現(xiàn)為樹狀或環(huán)狀,其相關(guān)研究主要著重于一維問題,但是需要考慮分汊點(diǎn)處的水流銜接情況,從而使研究問題復(fù)雜化,由于三級(jí)聯(lián)解法穩(wěn)定性好、計(jì)算速度快,是當(dāng)前一種成熟的主流算法[9-10]。

三峽水庫(kù)屬于大型河道水庫(kù),在防洪、航運(yùn)以及發(fā)電等方面產(chǎn)生了積極作用[11-12]。目前,國(guó)內(nèi)學(xué)者多采用MIKE11、HEC-RAS等模型對(duì)三峽庫(kù)區(qū)水動(dòng)力過程模擬[13-16],在水庫(kù)調(diào)度[17-20]及洪水傳播規(guī)律[21-23]等方面取得較多成果。通過文獻(xiàn)分析發(fā)現(xiàn),先前研究區(qū)域主要針對(duì)寸灘至壩址河段,沒有全面考慮各支流入流和區(qū)間入流對(duì)庫(kù)區(qū)洪水演進(jìn)的影響。此外,MIKE、HEC等軟件功能齊全、較為穩(wěn)定,但是無法對(duì)模型進(jìn)行修改并耦合水文預(yù)報(bào)模型,進(jìn)而更好的應(yīng)用于三峽庫(kù)區(qū)水動(dòng)力模擬與水情實(shí)時(shí)預(yù)報(bào),因此開發(fā)適用于三峽庫(kù)區(qū)的一維河網(wǎng)水動(dòng)力模型具有重要意義。

基于此,本文采用Preissmann四點(diǎn)隱式差分格式離散圣維南方程,應(yīng)用三級(jí)聯(lián)解法構(gòu)建了一種高效、精準(zhǔn)和可靠的水動(dòng)力數(shù)值模型。針對(duì)三峽庫(kù)區(qū)復(fù)雜河網(wǎng)系統(tǒng),構(gòu)建了長(zhǎng)江干流朱沱至壩址河段以及包括嘉陵江、烏江在內(nèi)共21條支流的三峽庫(kù)區(qū)水動(dòng)力模型,并結(jié)合水文預(yù)報(bào)模型數(shù)據(jù),充分考慮支流入流和區(qū)間入流的影響,模擬三峽庫(kù)區(qū)洪水演進(jìn)過程并探究其傳播規(guī)律,為三峽庫(kù)區(qū)水文預(yù)報(bào)及洪水調(diào)度提供指導(dǎo)。此外本模型可以用于其他河網(wǎng)系統(tǒng),對(duì)河網(wǎng)水動(dòng)力相關(guān)問題的研究具有積極推動(dòng)作用。

1 水動(dòng)力模型構(gòu)建

1.1 控制方程

圣維南方程是描述水流在河道中運(yùn)動(dòng)的一維非恒定流基本方程組,其包含連續(xù)性方程(1)和動(dòng)量守恒方程(2),形式如下:

(1)

(2)

1.2 控制方程的離散及求解

圣維南方程屬于非線性雙曲型偏微分方程,該方程只有在理想情況下才可求得解析解,因此在實(shí)際應(yīng)用中常采用數(shù)值近似解[24-25]。本文采用Preissmann四點(diǎn)隱式差分格式離散方程,其網(wǎng)格點(diǎn)結(jié)構(gòu)見圖1。

圖1 Preissmann四點(diǎn)隱式差分格式

若以f(x,t)表示式(1)、(2)中待計(jì)算的物理量Q(x,t)和z(x,t),采用fn+1=fn+Δf的形式表示相鄰時(shí)間步長(zhǎng)的因變量函數(shù)值,則式(1)、(2)中函數(shù)和導(dǎo)數(shù)的表達(dá)式見式(3)—(5):

(3)

(4)

(5)

式中n——時(shí)間層;j——空間層;Δt——時(shí)間步長(zhǎng);Δx——空間步長(zhǎng);θ——加權(quán)因子,由穩(wěn)定性決定,0.5≤θ≤1。

采用Preissmann四點(diǎn)隱式差分格式在時(shí)間和空間上離散圣維南方程組,以時(shí)段初的系數(shù)值代替時(shí)段平均值線性化方程組可得:

a1jΔzj + 1+b1jΔQj +1+c1jΔzj+d1jΔQj=e1j

(6)

a2jΔzj+1+b2jΔQj+1+c2jΔzj+d2jΔQj=e2j

(7)

式中 Δzj、Δzj+1——第j、j+1斷面在Δt時(shí)間內(nèi)水位增量;ΔQj、ΔQj+1——第j、j+1斷面在Δt時(shí)間內(nèi)流量增量;a1j、b1j、c1j、d1j、e1j、a2j、b2j、c2j、d2j、e2j——時(shí)間步長(zhǎng)Δt內(nèi)河段斷面j的差分方程系數(shù),具體表達(dá)式如下。

(8)

根據(jù)外河道首末斷面關(guān)系,設(shè)有如下的線性方程:

ΔQj=FjΔzj+Gj

(9)

Δzj=HjΔQj+1+IjΔzj+1+Jj

(10)

式(8)、(9)中的追趕系數(shù)表達(dá)式如下:

(11)

針對(duì)樹狀、環(huán)狀等復(fù)雜河網(wǎng),在式(9)、(10)的基礎(chǔ)上,采用三級(jí)算法進(jìn)行求解,具體步驟如下:①將河網(wǎng)劃分成河道與節(jié)點(diǎn),并以河道-節(jié)點(diǎn)-河道的形式構(gòu)建計(jì)算模型,每個(gè)河道皆有若干計(jì)算斷面組成,針對(duì)每個(gè)計(jì)算斷面采用有限差分法離散圣維南方程組,進(jìn)而得到水位和流量2個(gè)自變量的單一河道差分方程組;②考慮河道上下游邊界條件和節(jié)點(diǎn)連接情況,形成封閉的節(jié)點(diǎn)水位方程組,通過超松弛迭代法求解方程組進(jìn)而得到各節(jié)點(diǎn)水位;③將各節(jié)點(diǎn)水位回帶至相對(duì)應(yīng)河道,最后求解河道各斷面的水位和流量。

1.3 邊界條件

根據(jù)河道上下游關(guān)系,可分為上游邊界和下游邊界,涉及的邊界條件有3種,分別是水位邊界、流量邊界和水位流量關(guān)系邊界,求解1.2節(jié)循環(huán)計(jì)算方程(9)、(10)時(shí)需要給出相應(yīng)的邊界條件,邊界條件不同時(shí),將有不同的F1和G1的初始值,具體見式(12)—(15):

定解條件:z(x,t)|x=0=z(0,t)

Q1(x,t)|x=k=Q(k,t) 或z(x,t)|x=k=z(k,t)

或Q(x,t)|x=k,0=f(z(x,t))

(12)

(13)

(14)

(15)

2 三峽庫(kù)區(qū)水動(dòng)力模型構(gòu)建與驗(yàn)證

2.1 模型構(gòu)建及率定

三峽水庫(kù)是典型的大型河道型水庫(kù),從朱沱至壩址約750 km,庫(kù)區(qū)沿程斷面及測(cè)站位置分布見圖2。本文對(duì)始于朱沱站止于壩址的長(zhǎng)江干流,以及嘉陵江、烏江等21條主要支流構(gòu)建三峽庫(kù)區(qū)的一維水動(dòng)力模型,模型共計(jì)654處斷面。上游邊界在有水文站的斷面給定實(shí)測(cè)流量過程,如朱沱、北碚、武隆等,其他沒有水文站的斷面給定水文預(yù)報(bào)模型的流量數(shù)據(jù);下游邊界為壩址區(qū)域,給定實(shí)測(cè)水位過程線。此外,區(qū)間入流流量數(shù)據(jù)由水文預(yù)報(bào)模型提供,根據(jù)河道劃分分段接入水動(dòng)力模型。

圖2 三峽庫(kù)區(qū)河網(wǎng)示意

三峽庫(kù)區(qū)河道較長(zhǎng),沿程地形變化明顯,河道糙率變化大。為保證模型計(jì)算的高效性和準(zhǔn)確性,在模型中分段設(shè)置河道糙率參數(shù),選擇三峽庫(kù)區(qū)沿程朱沱、寸灘、長(zhǎng)壽、清溪場(chǎng)、忠縣、萬縣、奉節(jié)、巫山、巴東等重要站點(diǎn)的實(shí)測(cè)水位數(shù)據(jù)對(duì)庫(kù)區(qū)干流河道糙率進(jìn)行率定。首先在朱沱至壩址所有干流河道給定0.04的固定糙率值,在模型率定過程中,不斷改變各個(gè)河道的糙率值直至達(dá)到最優(yōu)解,本文所率定河道糙率值與相關(guān)文獻(xiàn)[21]中2套糙率對(duì)比見表1。

表1 三峽庫(kù)區(qū)各河段糙率參數(shù)

2.2 模擬驗(yàn)證

本文對(duì)三峽庫(kù)區(qū)2020年8—9月、2021年8—9月兩場(chǎng)洪水過程進(jìn)行模擬,各河道糙率系數(shù)采用率定所得結(jié)果,庫(kù)區(qū)干流上游邊界采用朱沱水文站實(shí)測(cè)流量數(shù)據(jù),支流上游邊界有水文站斷面采用各自上游水文站實(shí)測(cè)流量數(shù)據(jù),水文預(yù)報(bào)模型給定未有水文站的支流斷面流量數(shù)據(jù)以及干支流區(qū)間入流流量數(shù)據(jù),下游邊界采用壩址實(shí)測(cè)水位數(shù)據(jù),并對(duì)庫(kù)區(qū)干流沿程各個(gè)站點(diǎn)模擬的水位和流量進(jìn)行驗(yàn)證分析。具有代表性的站點(diǎn)(寸灘站、清溪場(chǎng)站、萬縣站)2020年8—9月、2021年8—9月模擬結(jié)果和實(shí)測(cè)值的對(duì)比見圖3、4,庫(kù)區(qū)沿程水面線模擬與實(shí)測(cè)結(jié)果對(duì)比見圖5。分析庫(kù)區(qū)沿程各站點(diǎn)的模擬結(jié)果發(fā)現(xiàn):水位平均絕對(duì)誤差在0.1～0.3 m,根據(jù)模擬水位值所繪制庫(kù)區(qū)水面線與實(shí)際水面線吻合良好,流量平均相對(duì)誤差在2%～6%。結(jié)果顯示本文模型能準(zhǔn)確模擬庫(kù)區(qū)水面線、水情要素的變化以及洪水傳播過程,說明三峽庫(kù)區(qū)河道糙率率定的結(jié)果準(zhǔn)確,模型構(gòu)建合理,具有較高的模擬精度。

圖3 2020年8—9月模擬結(jié)果與實(shí)測(cè)值對(duì)比

圖4 2021年8—9月模擬結(jié)果與實(shí)測(cè)值對(duì)比

圖5 三峽水庫(kù)沿程水面線變化

3 三峽庫(kù)區(qū)洪水傳播規(guī)律分析

恒定流計(jì)算是河網(wǎng)數(shù)值模擬中對(duì)輸入邊界的理想化處理,幾乎無法出現(xiàn)在實(shí)際工程中,但是一定程度上可以反映出三峽庫(kù)區(qū)洪水的演進(jìn)過程,可為水庫(kù)調(diào)度人員了解庫(kù)區(qū)的洪水演進(jìn)規(guī)律提供參考。本文針對(duì)不同組合流量和水位條件對(duì)庫(kù)區(qū)洪水傳播規(guī)律進(jìn)行了探討。

3.1 庫(kù)區(qū)水面線變化規(guī)律分析

以朱沱、北碚、武隆和區(qū)間來水的恒定流量組合代表不同來水情景,模擬不同壩址水位條件下各來水組合的庫(kù)區(qū)水面線,具體見圖6。圖6a給出了上游入流流量恒定時(shí),控制下游壩址水位條件下庫(kù)區(qū)水面線變化情況,結(jié)果顯示:隨著壩址水位的升高,庫(kù)區(qū)水面線不斷抬升,近壩區(qū)間水位抬升幅度大于庫(kù)尾區(qū)域,庫(kù)區(qū)水面線比降從上游至下游呈減小趨勢(shì),壩址水位對(duì)寸灘站水位的影響已經(jīng)較小,對(duì)朱沱站的影響甚至可以忽略。圖6b—6e給出了下游壩址水位恒定時(shí),控制上游入流流量條件下庫(kù)區(qū)水面線變化情況,結(jié)果顯示:當(dāng)朱沱流量增大時(shí),清溪場(chǎng)以上河段水位的抬升最為明顯,且抬升幅度相對(duì)一致,清溪場(chǎng)以下河段水位抬升幅度逐步減少,水面線比降不斷增大;當(dāng)北碚流量增大時(shí),庫(kù)區(qū)水面線的中部有所抬升,其中寸灘站到清溪場(chǎng)站抬升最為明顯,清溪場(chǎng)站下游水面線比降增大;當(dāng)武隆入流流量增大時(shí),庫(kù)區(qū)水面線的中部有所抬升,其中清溪場(chǎng)站抬升最為明顯,清溪場(chǎng)站上游水面線比降不斷減小,其下游水面線比降不斷增大;當(dāng)控制區(qū)間入流流量時(shí),庫(kù)區(qū)水面線的變化規(guī)律與控制武隆、北碚來水情景類似,但總體變化幅度較小。

a)控制壩址水位

b)控制朱沱流量

c)控制北碚流量圖6 各情景庫(kù)區(qū)沿程水面線變化

d)控制武隆流量

e)控制區(qū)間流量續(xù)圖6 各情景庫(kù)區(qū)沿程水面線變化

總體而言,從水面線的沿程分布來看,壩址水位對(duì)庫(kù)區(qū)沿程各站存在頂托作用,最遠(yuǎn)影響至寸灘站附近;入流流量增大的影響主要體現(xiàn)在中上游河段水位的抬升,朱沱至清溪場(chǎng)河段水位抬升最為明顯,水面線形態(tài)較陡且比降減小,清溪場(chǎng)至壩址河段水面線逐步平緩,且水面線比降隨入流流量增加而增大,其中忠縣至萬縣河段水面線較為平緩,巴東至壩址河段水位基本呈水平漲落。

3.2 庫(kù)區(qū)洪水傳播時(shí)間分析

為研究庫(kù)區(qū)洪水傳播時(shí)間規(guī)律,根據(jù)實(shí)況洪水過程模擬朱沱站來水情況:設(shè)定洪水周期為14 d,起漲至洪峰歷時(shí)7 d,起漲流量為5 000 m3/s,洪峰流量分別為10 000、20 000、30 000、40 000、50 000、60 000 m3/s,同時(shí)考慮壩址水位對(duì)傳播時(shí)間的影響,從低水位145 m至高水位170 m每5 m劃分水位級(jí),共模擬36種洪水演進(jìn)場(chǎng)景。

根據(jù)模擬結(jié)果,統(tǒng)計(jì)不同洪峰流量與壩址水位條件下,朱沱站來水至壩址洪水傳播時(shí)間,結(jié)果見圖7。洪峰流量為10 000 m3/s時(shí),壩址水位為145 m時(shí)洪水傳播時(shí)間為42 h,壩址水位為170 m時(shí)傳播時(shí)間為23 h;朱沱洪峰流量為60 000 m3/s時(shí),壩址水位為145 m時(shí)洪水傳播時(shí)間為44 h,壩址水位為170 m時(shí)傳播時(shí)間為31 h。結(jié)果表明:在相同洪峰流量條件下,壩址水位越高傳播時(shí)間越短,其原因在于水深隨水位增加而不斷增加,動(dòng)力波波速與水深呈正比,因此洪水傳播時(shí)間縮短;當(dāng)壩址水位相同時(shí),洪峰流量越大傳播時(shí)間越長(zhǎng),其原因在于流量增大導(dǎo)致運(yùn)動(dòng)波和動(dòng)力波的傳播速度都增大,但以運(yùn)動(dòng)波特性傳播距離更長(zhǎng),從而增加了洪水傳播時(shí)間。通過查閱相關(guān)文獻(xiàn)[23],本文結(jié)論與文獻(xiàn)結(jié)論相一致,可為庫(kù)區(qū)洪水預(yù)報(bào)和防汛工作提供參考依據(jù)。

圖7 三峽庫(kù)區(qū)洪水傳播時(shí)間

3.3 庫(kù)區(qū)測(cè)站水位流量關(guān)系分析

實(shí)際工程中流量測(cè)驗(yàn)通常采用流量計(jì)法、容積法、浮標(biāo)法等方法,水位測(cè)驗(yàn)通常采用水尺或水位計(jì)。水位測(cè)驗(yàn)工作簡(jiǎn)單,連續(xù)的水位數(shù)據(jù)獲取較為容易,相較而言流量測(cè)驗(yàn)工作技術(shù)復(fù)雜、耗資較大且難以獲取連續(xù)的流量數(shù)據(jù)。在水文預(yù)報(bào)、水動(dòng)力計(jì)算過程中,常常通過水位流量關(guān)系將連續(xù)的水位數(shù)據(jù)轉(zhuǎn)換、推算成連續(xù)的流量數(shù)據(jù)。三峽庫(kù)區(qū)沿程測(cè)站除寸灘站有連續(xù)的流量數(shù)據(jù)外,僅清溪場(chǎng)站、萬縣站有少量不連續(xù)的流量數(shù)據(jù),其余各站皆無流量數(shù)據(jù)。寸灘站位于重慶市江北區(qū)三家灘,距朱沱站約150 km,距三峽壩址約605 km,是長(zhǎng)江上游最重要的控制站之一,該站數(shù)據(jù)時(shí)間序列長(zhǎng)、測(cè)驗(yàn)精度高,可代表重慶市區(qū)水情變化。因此,以寸灘站實(shí)測(cè)數(shù)據(jù)為支撐,分析該站水位流量關(guān)系,有助于探求其他測(cè)站的水位流量關(guān)系,并進(jìn)一步推算流量過程。

采用2020年8—9月、2021年8—9月兩場(chǎng)洪水過程模擬的寸灘站水位、流量數(shù)據(jù)繪制水位流量關(guān)系,并與實(shí)測(cè)值進(jìn)行比較,結(jié)果見圖8。寸灘站汛期水位一般在165～190 m變化,本文選取不同的水位區(qū)間,對(duì)其水位流量關(guān)系進(jìn)行分析:2020年9月1日0時(shí)至4日12時(shí),水位變化為170、175、170 m,流量變化為23 000、36 000、19 000 m3/s,漲水階段用時(shí)24 h,退水階段用時(shí)60 h,為逆時(shí)針繩套曲線;2021年9月2日0時(shí)至6日0時(shí),水位變化為175、183、175 m,流量變化為29 000、48 000、22 000 m3/s,漲水階段用時(shí)37 h,退水階段用時(shí)83 h,為逆時(shí)針繩套曲線;2020年8月17日0時(shí)至22日6時(shí)。水位變化為180、190、180 m,流量變化為49 000、81 000、3 800 m3/s,漲水階段用時(shí)70 h,退水階段用時(shí)56 h,為逆時(shí)針繩套曲線。通過分析結(jié)果,寸灘站水位在170～175 m時(shí)水位流量關(guān)系十分復(fù)雜,其他水位時(shí)則表現(xiàn)出較為簡(jiǎn)單的水位流量對(duì)應(yīng)關(guān)系,隨著水位升高,流量增大趨勢(shì)更為明顯。此外,寸灘站受洪水漲落影響較大,多表現(xiàn)為逆時(shí)針繩套曲線?？傮w而言,寸灘站的模擬結(jié)果與實(shí)際水位流量關(guān)系較為吻合,因此通過構(gòu)建三峽庫(kù)區(qū)一維水動(dòng)力模型能精確地計(jì)算三峽庫(kù)區(qū)沿程測(cè)站的水位流量關(guān)系,可直觀、準(zhǔn)確地為庫(kù)區(qū)汛期的水文計(jì)算及預(yù)報(bào)提供服務(wù)。

a)2020年8—9月

b)2021年8—9月圖8 寸灘站水位流量關(guān)系對(duì)比

4 結(jié)論

本文構(gòu)建三峽庫(kù)區(qū)一維水動(dòng)力模型,對(duì)庫(kù)區(qū)2020年8—9月、2021年8—9月洪水的流量和水位進(jìn)行模擬驗(yàn)證,各站點(diǎn)水位平均絕對(duì)誤差小于0.3 m,流量平均相對(duì)誤差不超過6%,計(jì)算的庫(kù)區(qū)沿程水面線與實(shí)測(cè)值基本一致。采用模型對(duì)洪水傳播規(guī)律進(jìn)行探究,主要得出以下結(jié)論。

a)以朱沱站、北碚站、武隆站和區(qū)間流量組成上游來水,結(jié)合不同壩址水位模擬庫(kù)區(qū)沿程水面線,結(jié)果顯示:壩址水位頂托作用影響至寸灘站附近;上游高洪來水使庫(kù)尾水位抬升明顯,巴東至壩址河段水位呈水平漲落。

b)結(jié)合實(shí)際洪水過程,分別考慮朱沱站來水情況與壩址水位條件,共模擬36種洪水演進(jìn)場(chǎng)景,相同條件下,水位越高庫(kù)區(qū)洪水傳播時(shí)間越短,流量越大則庫(kù)區(qū)洪水傳播時(shí)間越長(zhǎng),整體傳播時(shí)間為23～44 h。

c)寸灘站是三峽水庫(kù)的干流入庫(kù)控制站,汛期水位在165～190 m,受洪水漲落影響較大,多表現(xiàn)為逆時(shí)針繩套曲線。寸灘站模擬所得水位流量關(guān)系與實(shí)際水位流量關(guān)系較為吻合,因此本模型可模擬庫(kù)區(qū)沿程各站點(diǎn)水位流量關(guān)系,并為各站點(diǎn)水文預(yù)報(bào)及洪水調(diào)度提供依據(jù)。