基于平面二維淺水方程的海河干流行洪模擬研究

李 婉 亭

(中國市政工程華北設計研究總院有限公司，天津 300074)

天津市地處海河流域最下游，承擔著流域洪水宣泄入海的任務，海河流域洪水的預見期短，一旦暴發(fā)洪水，天津市極易形成洪災[1]。由于天津市地處特殊的地理位置，受到渤海的頂托作用，流經天津市入海的河道極易受到洪、澇、潮的共同影響。應對平原河網地區(qū)洪澇潮的方法，包括工程措施與非工程措施，其中非工程措施方面，重在對工程措施的合理調度。董炳江[2]等建立了在非結構非結構混合網格上求解水深平均的二維淺水方程的模型。潘存鴻[3]等采用特殊的底坡源項處理技術，建立了三角形網格下二階精度的KFVS和諧格式。劉剛[4]等對經典的Roe格式計算數值通量中的靜水壓力項進行了修正，并證明了修正后的Roe格式具備和諧性。許棟[5]等提出了一種河道灘地洪水淹沒分析的多分辨率處理方法，有效將二維淺水模擬的低分辨率網格和局部灘地高分辨率地形結合。本研究利用平面二維淺水方程，及經驗證的城市暴雨洪水管理模型SWMM[6]，以海河干流段為研究對象，進行不同重現期設計暴雨下2種不同情景下的河道行洪模擬，旨在為城市河道防洪提供技術支持。

1 平面二維淺水方程基本理論

1.1 平面二維淺水方程的適用性

二維淺水方程適用于水平尺度遠大于垂直尺度、無明顯垂直環(huán)流、平面大范圍的自由表面流動，本文所研究海河干流河長約73 km，河道寬度為100～400 m，水深5 m左右，且河道不存在明顯環(huán)流，屬表面自由流動，水平尺度上遠大于垂直尺度，可以引入靜水壓力假設，并沿水深方向進行積分來簡化方程[7,8]。

1.2 控制方程

連續(xù)性方程：

(1)

X方向動量方程：

(2)

Y方向動量方程：

(3)

式中：h為總水深，m；η為河底高程，m；d為靜水深，m；t為時間，s；u、v分別為x、y方向的平均流速；S為源匯項，m3/s；us、vs為源匯項水流流速；f為科式力系數；ρ0為水的相對密度；pa為大氣壓強；Txx、Txy、Tyy由渦黏性公式計算：

1.3 數值解法

(1)空間離散。采用有限體積法離散求解[9,10]，黎曼近似解通量差分裂格式計算通過各單元邊的水流法向數值通量[11]，對每個離散體沿邊界的輸入、輸出通量，利用控制方程，求解出水深、流速等參數。

在笛卡爾坐標系下，二維淺水方程組寫為：

(4)

式中：I、V分別為無黏性和黏性的通量，各項分別為：

運用高斯定理，對式(4)第i個單元進行積分，得出：

(5)

式中：A為控制體的單元面積；L為控制體的單元邊界；ds為沿邊界的積分變量。則方程(5)可寫為：

(6)

式中：Ui為單元i的U的平均值；Si為單元i的S的平均值；N為單元的邊界數；ΔLi為單元i的邊界長度。

(2)時間積分。

二維淺水方程的求解方法分為低階方法與高階方法兩種。

低階方法為顯式的Euler方法：

Un+1=Un+tG(Un)

高階方法為以如下形式的使用了二階的Runge Kutta方法：

(3)邊界條件。開邊界：上下游邊界處具體的流量過程、水位過程等設置為開邊界；閉合邊界：設置陸地邊界為閉合邊界；干濕邊界：判斷網格單元處的水深，將計算單元判定為干單元、半干單元、濕單元分別計算相應質量與動量。

2 水動力模型的建立與求解

2.1 模型計算區(qū)域及離散

(1)建立工作區(qū)域與生成邊界線。模型的計算范圍由北運河下游、子牙河下游、海河干流組成。利用衛(wèi)星遙感圖像在GIS中提取模型的各個邊界線見圖1。

圖1 研究區(qū)域邊界線

(2)網格化分與地形插值。天然水體的邊界及水下地形很不規(guī)則，本研究采用無結構的網格劃分，因其最為靈活，普適性好[12]，河道順直段劃分網格采用四邊形網格，在河道交匯處及河道變化較大處則采用三角形網格，最終得到研究區(qū)劃分單元數5 306個，節(jié)點6 722個。地形信息需通過掌握的河道斷面資料進行散點插值，在已劃分好的網格上布設散點高程信息，得到最終進行模型計算的基礎地形文件。

2.2 水動力模型參數率定與驗證

模型計算率定的主要參數為河床糙率，通過調節(jié)糙率大小使模擬結果與實測值相吻合，選擇結果吻合度最好的河床糙率對模型進一步驗證[13]。

2.2.1 參數率定

(1)邊界條件：模型邊界將上游西河閘、屈家店閘設為流量邊界，下游海河閘邊界設為水位邊界。采用2012年天津市主汛期(7.21～8.10)的實測水位及流量資料；模型的旁側入流利用SWMM模型模擬計算的各泵站排水結果。

(2)閘門調度：模型的內部控制邊界(二道閘、海河閘)采用2012年的水閘啟閉調度資料控制。

(3)初始條件：初始水位采用2012年7月21日8時的實測水位，初始流量取為0 m3/s，河床糙率設置為0.025，并根據計算結果進行調整。

水動力模型整體上符合要求，然而由于海河干流實際情況較模型復雜許多，加之實測值在測量中也存在一定誤差，因此模擬值并不能夠跟實測值完全相同，但模擬值與實測值水位、流量變化趨勢基本一致[13]，最終得率定結果如圖2所示。

圖 2 二道閘閘上水位、二道閘過閘流量、海河閘過閘流量模擬值與實測值對比圖

2.2.2 參數驗證

根據現有資料并結合天津市2012年之后的各年水資源公報，有2013年、2014年、2017年屬偏枯水年，2016年屬偏豐水年，2015年屬平水年，為了參數驗證的合理性，選擇平水年即2015年進行驗證，同模型參數率定的步驟類似，其驗證結果如圖3所示。

圖3 二道閘閘上水位、二道閘過閘流量、海河閘過閘流量模擬值與實測值對比圖

3 不同重現期設計暴雨下河道匯流計算

在海河干流現有控制性工程的基礎上，設計在下游海河閘出口處增設泵站，旨在研究現狀工程系統(tǒng)、增設泵站抽排水的工程系統(tǒng)對海河干流行洪的不同影響，為海河流域減輕行洪壓力提供技術支持。具體情景設計為：情景1：不同重現期設計暴雨下的現狀河道行洪計算；情景2：情景1的基礎上，在下游海河閘處增設抽水排澇泵站的河道行洪計算。

3.1 計算條件

(1)初始條件與邊界條件：考慮到2012年屬近30年來的罕見洪澇狀況，根據2012年相應資料，設置模型初始條件。邊界條件上游為流量過程(0 m3/s)，下游為水位過程(2012年汛期出現的最高潮位記錄)[14,15]。旁側入流采用SWMM模型計算的流量過程，另外在情景2計算時，在海河閘處需增設一泵站(流量控制)。

(2)閘門調度：流域中控制水流運動的堰、閘等的過流流量滿足水力學上的堰流公式等[16]，本研究中涉及的閘門的啟閉通過閘上、閘下的水位差控制：大于0.01 m時，閘門開啟，小于等于0.01 m時，閘門關閉。

3.2 計算結果及分析

根據對現狀河道資料的分析，選取行洪最不利的北運河段上的北洋橋斷面、內部邊界二道閘作為研究對象，模擬計算結果表明，不同重現期設計暴雨下(降雨最大值出現在9 h)，北洋橋處斷面的水位變化在9 h前后顯著加快，當遭遇3/a設計暴雨時，該處堤防尚能防止洪水漫堤，一旦遭遇3/a以上的暴雨時，該處堤防則會發(fā)生洪水漫堤。

圖4 不同重現期設計暴雨下北洋橋處水位、二道閘閘上水位、二道閘過閘流量過程

在海河閘處增設抽水排澇泵站，研究其對所選取的研究斷面對象(北洋橋、二道閘閘上)的影響，由于篇幅限制僅列出3 a一遇下的2處研究對象的水位變化圖。

圖5 3 a一遇設計暴雨不同情境下北洋橋處、二道閘閘上水位變化過程

表1 二道閘閘上最高水位

由上面計算結果可知，抽水排澇泵站可在一定程度上減輕河道的行洪壓力，所選取的研究斷面據下游的距離從大到小分別為：北洋橋、二道閘、海河閘，其中：不同重現期設計暴雨下,北洋橋斷面處的最高水位的減小幅度較??；二道閘閘上最高水位的最大變化量(抽水情景與無泵站進行對比)分別為0.147、0.148、0.150、0.151 m，平均為0.149 m；海河閘閘上最高水位的最大變化量(抽水情景與無泵站進行對比)分別為0.210、0.210、0.196、0.207 m，平均最大變化量為0.201 m。

表2 海河閘閘上最高水位

4 結 語

(1)平面二維淺水方程模擬河道洪水演進的結果較好，但是模擬計算結果可信度是建立在資料充分的情況下，河道地形數據、水文氣象數據等是模型建立的重要前提，且對模型的率定驗證也是保證模型后續(xù)應用準確性的重要步驟。

(2)不同重現期設計暴雨下，海河干流的堤防薄弱處存在洪水漫堤的風險，設計通過在下游海河閘處增設抽水泵站，研究其對各控制斷面最高水位的影響，結果表明受洪澇潮影響的城市河道行洪壓力在下游增設抽水排澇泵站的情景下可得到有效緩解。

(3)雖然增設泵站可以在一定程度上減輕河道的行洪壓力，起到降低最高水位的作用，但是對于距離下游泵站較遠的上游堤防薄弱段，仍存在洪水漫堤的危險，這說明由于天津市城市地區(qū)河道行洪的特殊性，單一的工程措施已經很難保障其河道的行洪安全，需要通過更多樣化的措施(如海綿措施、及時的分洪措施、智慧水系的系統(tǒng)調控措施等)來進行城市地區(qū)的防洪保護。

□