王 玨,張海榮,孫振宇,張明亮,3
(1.大連海洋大學海洋科技與環境學院,遼寧 大連 116023;2.中國長江電力股份有限公司水資源研究中心,湖北 宜昌 443133;3.大連理工大學海岸和近海工程國家重點實驗室,遼寧 大連 116025)
洪水是由暴雨、急劇融冰化雪、風暴潮等引起的江河湖泊水量激增、水位上漲的一種自然災害。通過對歷史資料分析發現:擴大耕地、圍湖造田等人類活動不斷破壞地表下墊面形態,致使洪水頻率升高,致災程度也不斷加深。在強降雨期間,短時間內流域水量突增、匯聚,形成洪水在河網水系內傳播。針對河網內突發的洪水問題,國內外諸多學者開展了大量的研究工作。早期學者多采用瞬時流態法[1]、馬斯京根法[2-3],由于上述方法過度簡化圣維南方程,在精度上往往不能滿足實際需求。考慮到河網水動力模型的核心是對圣維南方程組的求解,Preissmann和Cunge[4]提出四點隱式差分格式,經過不斷完善,因其收斂快、效率高并且穩定性好,被廣泛應用于實際工程中[5-8]。一般來說,河網水系常表現為樹狀或環狀,其相關研究主要著重于一維問題,但是需要考慮分汊點處的水流銜接情況,從而使研究問題復雜化,由于三級聯解法穩定性好、計算速度快,是當前一種成熟的主流算法[9-10]。
三峽水庫屬于大型河道水庫,在防洪、航運以及發電等方面產生了積極作用[11-12]。目前,國內學者多采用MIKE11、HEC-RAS等模型對三峽庫區水動力過程模擬[13-16],在水庫調度[17-20]及洪水傳播規律[21-23]等方面取得較多成果。通過文獻分析發現,先前研究區域主要針對寸灘至壩址河段,沒有全面考慮各支流入流和區間入流對庫區洪水演進的影響。此外,MIKE、HEC等軟件功能齊全、較為穩定,但是無法對模型進行修改并耦合水文預報模型,進而更好的應用于三峽庫區水動力模擬與水情實時預報,因此開發適用于三峽庫區的一維河網水動力模型具有重要意義。
基于此,本文采用Preissmann四點隱式差分格式離散圣維南方程,應用三級聯解法構建了一種高效、精準和可靠的水動力數值模型。針對三峽庫區復雜河網系統,構建了長江干流朱沱至壩址河段以及包括嘉陵江、烏江在內共21條支流的三峽庫區水動力模型,并結合水文預報模型數據,充分考慮支流入流和區間入流的影響,模擬三峽庫區洪水演進過程并探究其傳播規律,為三峽庫區水文預報及洪水調度提供指導。此外本模型可以用于其他河網系統,對河網水動力相關問題的研究具有積極推動作用。
圣維南方程是描述水流在河道中運動的一維非恒定流基本方程組,其包含連續性方程(1)和動量守恒方程(2),形式如下:
(1)
(2)

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

圖1 Preissmann四點隱式差分格式
若以f(x,t)表示式(1)、(2)中待計算的物理量Q(x,t)和z(x,t),采用fn+1=fn+Δf的形式表示相鄰時間步長的因變量函數值,則式(1)、(2)中函數和導數的表達式見式(3)—(5):
(3)
(4)
(5)
式中n——時間層;j——空間層;Δt——時間步長;Δx——空間步長;θ——加權因子,由穩定性決定,0.5≤θ≤1。
采用Preissmann四點隱式差分格式在時間和空間上離散圣維南方程組,以時段初的系數值代替時段平均值線性化方程組可得:
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時間內水位增量;ΔQj、ΔQj+1——第j、j+1斷面在Δt時間內流量增量;a1j、b1j、c1j、d1j、e1j、a2j、b2j、c2j、d2j、e2j——時間步長Δt內河段斷面j的差分方程系數,具體表達式如下。
(8)
根據外河道首末斷面關系,設有如下的線性方程:
ΔQj=FjΔzj+Gj
(9)
Δzj=HjΔQj+1+IjΔzj+1+Jj
(10)
式(8)、(9)中的追趕系數表達式如下:
(11)

針對樹狀、環狀等復雜河網,在式(9)、(10)的基礎上,采用三級算法進行求解,具體步驟如下:①將河網劃分成河道與節點,并以河道-節點-河道的形式構建計算模型,每個河道皆有若干計算斷面組成,針對每個計算斷面采用有限差分法離散圣維南方程組,進而得到水位和流量2個自變量的單一河道差分方程組;②考慮河道上下游邊界條件和節點連接情況,形成封閉的節點水位方程組,通過超松弛迭代法求解方程組進而得到各節點水位;③將各節點水位回帶至相對應河道,最后求解河道各斷面的水位和流量。
根據河道上下游關系,可分為上游邊界和下游邊界,涉及的邊界條件有3種,分別是水位邊界、流量邊界和水位流量關系邊界,求解1.2節循環計算方程(9)、(10)時需要給出相應的邊界條件,邊界條件不同時,將有不同的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)
三峽水庫是典型的大型河道型水庫,從朱沱至壩址約750 km,庫區沿程斷面及測站位置分布見圖2。本文對始于朱沱站止于壩址的長江干流,以及嘉陵江、烏江等21條主要支流構建三峽庫區的一維水動力模型,模型共計654處斷面。上游邊界在有水文站的斷面給定實測流量過程,如朱沱、北碚、武隆等,其他沒有水文站的斷面給定水文預報模型的流量數據;下游邊界為壩址區域,給定實測水位過程線。此外,區間入流流量數據由水文預報模型提供,根據河道劃分分段接入水動力模型。

圖2 三峽庫區河網示意
三峽庫區河道較長,沿程地形變化明顯,河道糙率變化大。為保證模型計算的高效性和準確性,在模型中分段設置河道糙率參數,選擇三峽庫區沿程朱沱、寸灘、長壽、清溪場、忠縣、萬縣、奉節、巫山、巴東等重要站點的實測水位數據對庫區干流河道糙率進行率定。首先在朱沱至壩址所有干流河道給定0.04的固定糙率值,在模型率定過程中,不斷改變各個河道的糙率值直至達到最優解,本文所率定河道糙率值與相關文獻[21]中2套糙率對比見表1。

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

圖3 2020年8—9月模擬結果與實測值對比

圖4 2021年8—9月模擬結果與實測值對比

圖5 三峽水庫沿程水面線變化
恒定流計算是河網數值模擬中對輸入邊界的理想化處理,幾乎無法出現在實際工程中,但是一定程度上可以反映出三峽庫區洪水的演進過程,可為水庫調度人員了解庫區的洪水演進規律提供參考。本文針對不同組合流量和水位條件對庫區洪水傳播規律進行了探討。
以朱沱、北碚、武隆和區間來水的恒定流量組合代表不同來水情景,模擬不同壩址水位條件下各來水組合的庫區水面線,具體見圖6。圖6a給出了上游入流流量恒定時,控制下游壩址水位條件下庫區水面線變化情況,結果顯示:隨著壩址水位的升高,庫區水面線不斷抬升,近壩區間水位抬升幅度大于庫尾區域,庫區水面線比降從上游至下游呈減小趨勢,壩址水位對寸灘站水位的影響已經較小,對朱沱站的影響甚至可以忽略。圖6b—6e給出了下游壩址水位恒定時,控制上游入流流量條件下庫區水面線變化情況,結果顯示:當朱沱流量增大時,清溪場以上河段水位的抬升最為明顯,且抬升幅度相對一致,清溪場以下河段水位抬升幅度逐步減少,水面線比降不斷增大;當北碚流量增大時,庫區水面線的中部有所抬升,其中寸灘站到清溪場站抬升最為明顯,清溪場站下游水面線比降增大;當武隆入流流量增大時,庫區水面線的中部有所抬升,其中清溪場站抬升最為明顯,清溪場站上游水面線比降不斷減小,其下游水面線比降不斷增大;當控制區間入流流量時,庫區水面線的變化規律與控制武隆、北碚來水情景類似,但總體變化幅度較小。

a)控制壩址水位

b)控制朱沱流量

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

d)控制武隆流量

e)控制區間流量續圖6 各情景庫區沿程水面線變化
總體而言,從水面線的沿程分布來看,壩址水位對庫區沿程各站存在頂托作用,最遠影響至寸灘站附近;入流流量增大的影響主要體現在中上游河段水位的抬升,朱沱至清溪場河段水位抬升最為明顯,水面線形態較陡且比降減小,清溪場至壩址河段水面線逐步平緩,且水面線比降隨入流流量增加而增大,其中忠縣至萬縣河段水面線較為平緩,巴東至壩址河段水位基本呈水平漲落。
為研究庫區洪水傳播時間規律,根據實況洪水過程模擬朱沱站來水情況:設定洪水周期為14 d,起漲至洪峰歷時7 d,起漲流量為5 000 m3/s,洪峰流量分別為10 000、20 000、30 000、40 000、50 000、60 000 m3/s,同時考慮壩址水位對傳播時間的影響,從低水位145 m至高水位170 m每5 m劃分水位級,共模擬36種洪水演進場景。
根據模擬結果,統計不同洪峰流量與壩址水位條件下,朱沱站來水至壩址洪水傳播時間,結果見圖7。洪峰流量為10 000 m3/s時,壩址水位為145 m時洪水傳播時間為42 h,壩址水位為170 m時傳播時間為23 h;朱沱洪峰流量為60 000 m3/s時,壩址水位為145 m時洪水傳播時間為44 h,壩址水位為170 m時傳播時間為31 h。結果表明:在相同洪峰流量條件下,壩址水位越高傳播時間越短,其原因在于水深隨水位增加而不斷增加,動力波波速與水深呈正比,因此洪水傳播時間縮短;當壩址水位相同時,洪峰流量越大傳播時間越長,其原因在于流量增大導致運動波和動力波的傳播速度都增大,但以運動波特性傳播距離更長,從而增加了洪水傳播時間。通過查閱相關文獻[23],本文結論與文獻結論相一致,可為庫區洪水預報和防汛工作提供參考依據。

圖7 三峽庫區洪水傳播時間
實際工程中流量測驗通常采用流量計法、容積法、浮標法等方法,水位測驗通常采用水尺或水位計。水位測驗工作簡單,連續的水位數據獲取較為容易,相較而言流量測驗工作技術復雜、耗資較大且難以獲取連續的流量數據。在水文預報、水動力計算過程中,常常通過水位流量關系將連續的水位數據轉換、推算成連續的流量數據。三峽庫區沿程測站除寸灘站有連續的流量數據外,僅清溪場站、萬縣站有少量不連續的流量數據,其余各站皆無流量數據。寸灘站位于重慶市江北區三家灘,距朱沱站約150 km,距三峽壩址約605 km,是長江上游最重要的控制站之一,該站數據時間序列長、測驗精度高,可代表重慶市區水情變化。因此,以寸灘站實測數據為支撐,分析該站水位流量關系,有助于探求其他測站的水位流量關系,并進一步推算流量過程。
采用2020年8—9月、2021年8—9月兩場洪水過程模擬的寸灘站水位、流量數據繪制水位流量關系,并與實測值進行比較,結果見圖8。寸灘站汛期水位一般在165~190 m變化,本文選取不同的水位區間,對其水位流量關系進行分析:2020年9月1日0時至4日12時,水位變化為170、175、170 m,流量變化為23 000、36 000、19 000 m3/s,漲水階段用時24 h,退水階段用時60 h,為逆時針繩套曲線;2021年9月2日0時至6日0時,水位變化為175、183、175 m,流量變化為29 000、48 000、22 000 m3/s,漲水階段用時37 h,退水階段用時83 h,為逆時針繩套曲線;2020年8月17日0時至22日6時。水位變化為180、190、180 m,流量變化為49 000、81 000、3 800 m3/s,漲水階段用時70 h,退水階段用時56 h,為逆時針繩套曲線。通過分析結果,寸灘站水位在170~175 m時水位流量關系十分復雜,其他水位時則表現出較為簡單的水位流量對應關系,隨著水位升高,流量增大趨勢更為明顯。此外,寸灘站受洪水漲落影響較大,多表現為逆時針繩套曲線。總體而言,寸灘站的模擬結果與實際水位流量關系較為吻合,因此通過構建三峽庫區一維水動力模型能精確地計算三峽庫區沿程測站的水位流量關系,可直觀、準確地為庫區汛期的水文計算及預報提供服務。

a)2020年8—9月

b)2021年8—9月圖8 寸灘站水位流量關系對比
本文構建三峽庫區一維水動力模型,對庫區2020年8—9月、2021年8—9月洪水的流量和水位進行模擬驗證,各站點水位平均絕對誤差小于0.3 m,流量平均相對誤差不超過6%,計算的庫區沿程水面線與實測值基本一致。采用模型對洪水傳播規律進行探究,主要得出以下結論。
a)以朱沱站、北碚站、武隆站和區間流量組成上游來水,結合不同壩址水位模擬庫區沿程水面線,結果顯示:壩址水位頂托作用影響至寸灘站附近;上游高洪來水使庫尾水位抬升明顯,巴東至壩址河段水位呈水平漲落。
b)結合實際洪水過程,分別考慮朱沱站來水情況與壩址水位條件,共模擬36種洪水演進場景,相同條件下,水位越高庫區洪水傳播時間越短,流量越大則庫區洪水傳播時間越長,整體傳播時間為23~44 h。
c)寸灘站是三峽水庫的干流入庫控制站,汛期水位在165~190 m,受洪水漲落影響較大,多表現為逆時針繩套曲線。寸灘站模擬所得水位流量關系與實際水位流量關系較為吻合,因此本模型可模擬庫區沿程各站點水位流量關系,并為各站點水文預報及洪水調度提供依據。