基于數(shù)字流域系統(tǒng)的漳衛(wèi)河系洪水演進(jìn)模型建立

邢 斌，朱曉春，李永慶

（海河水利委員會科技咨詢中心，天津 300170）

洪水的運(yùn)動規(guī)律復(fù)雜多變，威脅人們正常的生產(chǎn)、生活，帶來嚴(yán)重的生命財產(chǎn)損失和安全隱患。利用現(xiàn)代計算機(jī)技術(shù)，洪水演進(jìn)模型能夠根據(jù)上游干支流的洪水過程及區(qū)間來水，模擬計算下游任意斷面的洪峰流量、出現(xiàn)時間、最高水位及洪水過程，同時能夠模擬水流在蓄滯洪區(qū)內(nèi)的運(yùn)動規(guī)律，為流域防洪規(guī)劃、洪水預(yù)報、控制工程調(diào)度等提供科學(xué)的決策依據(jù)。

目前，洪水演進(jìn)模擬方法主要有水文學(xué)法、歷史洪災(zāi)調(diào)查法、一維或二維水力學(xué)法等。其中，二維水力學(xué)模擬方法能夠分別從時間與空間角度模擬洪水演進(jìn)過程，提供更為直觀的洪水淹沒范圍變化，特別是一維、二維、零維耦合的水力學(xué)模型能在大流域尺度內(nèi)更好地為防洪調(diào)度的決策提供服務(wù)。

隨著水力學(xué)模型技術(shù)日趨成熟，建模平臺可以使用丹麥的MIKE 系列、荷蘭的De1ft3D 等商業(yè)化軟件，也可以采用我國自主開發(fā)的如天津大學(xué)、河海大學(xué)、清華大學(xué)等高校研制的洪水演進(jìn)模型。其中，河海大學(xué)王船海開發(fā)的數(shù)字流域系統(tǒng)應(yīng)用開發(fā)軟件，因其可視化的建模、需求方案的可定制、計算成果的在線查詢分析、動態(tài)顯示等功能，在太湖、長江、淮河等流域有很多成功應(yīng)用案例，如太湖流域洪水預(yù)報系統(tǒng)、三峽水庫動庫容預(yù)報調(diào)度系統(tǒng)、江蘇省里下河河網(wǎng)水動力模型系統(tǒng)等。筆者采用數(shù)字流域系統(tǒng)建立海河流域漳衛(wèi)河系洪水演進(jìn)模型，依據(jù)漳衛(wèi)河系的現(xiàn)行調(diào)度規(guī)劃模擬了漳衛(wèi)河系發(fā)生50年一遇洪水的演進(jìn)過程。

1 流域概況

漳衛(wèi)河系是海河流域最南部的防洪骨干水系，流經(jīng)山西、河南、河北、山東、天津四省一市，流域面積3.76 萬km2，西以太岳山為界，南接黃河、徒駭馬頰河，北界滏陽河，東達(dá)渤海。流域地形總的趨勢西高東低，地面坡度山區(qū)丘陵區(qū)為10‰～0.5‰、平原為0.1‰～0.3‰。

漳衛(wèi)河系由漳河、衛(wèi)河、衛(wèi)運(yùn)河、漳衛(wèi)新河、南運(yùn)河組成。流域上游有漳河和衛(wèi)河兩大支流。漳、衛(wèi)兩河于徐萬倉匯合后稱衛(wèi)運(yùn)河，至四女寺樞紐，以下為漳衛(wèi)新河和南運(yùn)河。漳衛(wèi)新河是漳衛(wèi)河水系的入海尾閭河道。主要的蓄滯洪區(qū)有良相坡、長虹渠、白寺坡等。

2 模型建立及參數(shù)率定

漳衛(wèi)河流域面積較大，河流眾多，蓄滯洪區(qū)分布其中。根據(jù)流域?qū)嶋H情況，洪水演進(jìn)模型需要建立一維、二維及零維模型。

2.1 模型原理

2.1.1 一維模擬

河道洪水波傳播屬不恒定流，其演變規(guī)律符合圣維南方程組，并采用偏微分方程的數(shù)值解法進(jìn)行水流的數(shù)值模擬，公式如下：

式中：q 為旁側(cè)入流（m3/s）；Q，A，B，Z 分別為河道斷面流量（m3/s），過水面積（m2），河寬（m）和水位（m）；Vx為旁側(cè)入流流速在水流方向上的分量，一般可以近似為0；K為流量模數(shù)，反映河道的實(shí)際過流能力；α為動量校正系數(shù)，是反映河道斷面流速分布均勻性的系數(shù)；g 為重力加速度（m/s2）；t 為時間（s）；x 為距上一斷面的距離（m）。

2.1.2 二維模擬

蓄滯洪區(qū)的水流采用二維淺水波方程來描述：

式中：Z為水位（m）；u，v分別為x與y方向上的流速（m/s）；U，V分別為x與y方向上的單寬流量（m3/s）為單寬流量的矢量，為它的模，；q為考慮降雨等因素的源項（m/s）；g 為重力加速度（m/s2）；c為謝才系數(shù)；f為柯氏力系數(shù)；wxt，wyt分別為風(fēng)應(yīng)力沿x和y方向的分量；t為時間（s）；h為水深（m）；ρ為空氣密度（kg/m3）；x，y為自變量坐標(biāo)。

對上述二維淺水波方程直接求解有一定的困難。因此，采用破開算子法將該方程分成如下兩分步方程組，然后分別對其采用合適的方法進(jìn)行求解。

第一分步為：

第二分步為：

對上面兩分步方程組的數(shù)值求解，采用直角坐標(biāo)系下非均勻矩形規(guī)則網(wǎng)格的控制體積法，具體方法如圖1所示。

圖1 二維差分網(wǎng)格示意

蓄滯洪區(qū)邊界條件一般有堤，四周封閉，在大水情況下通過口門、漫堤、閘堰等與河道干流進(jìn)行水量交換，其交換水量對行洪區(qū)講是邊界條件，對整個模型計算講是內(nèi)邊界條件。

由單元I流進(jìn)單元J的流量為QX=δX(ZI-ZJ)+βX，單元K流到單元J的流量為QY=δY(ZK-ZJ)+βY，對式（4）離散可得：

式中：A 為單元J 的面積（m2）；∑Qi表示包括降雨在內(nèi)單位時間內(nèi)流進(jìn)單元J水量的代數(shù)和（m3/s）；Z為單元J的水位（m）；Δt為計算步長（s）。

2.1.3 零維模擬

計算軟件對無地形資料的蓄洪區(qū)、湖泊等對象處理為零維區(qū)域，主要考慮水流在該單元區(qū)域內(nèi)的調(diào)蓄作用，至于其中的運(yùn)動規(guī)律則不考慮，僅僅考慮區(qū)域內(nèi)的水位變化。水位的變化規(guī)律遵循水量平衡原理，流入?yún)^(qū)域的凈水量等于區(qū)域內(nèi)的蓄量增量，即進(jìn)出蓄洪區(qū)內(nèi)的水量等于其內(nèi)的蓄量變化。其計算公式為：

式中： )(zA為蓄洪區(qū)內(nèi)水面面積（m2），一般為水位的某種函數(shù)關(guān)系；∑Q為蓄變量（m3/s）；Z為t時刻的水位（m）；t為時間（s）。

對式（6）離散后為：

式中：∑Q為蓄變量（m3/s）； )(zA為蓄洪區(qū)內(nèi)水面面積（m2）；Z0為初始水位（m）；Z 為t 時刻的水位（m）；t為時間（s）。

2.1.4 模型耦合

流域洪水運(yùn)動模擬由零維、一維、二維模擬所組成，各部分模擬必須耦合聯(lián)立才能求解，而各部分模擬的耦合是通過“聯(lián)系”來實(shí)現(xiàn)的?！奥?lián)系”要素是指流域中控制水流運(yùn)動的堰、閘及口門等，其過流流量一般以寬頂堰流量公式計算。

寬頂堰上的水流可分為自由出流、淹沒出流兩種流態(tài)，不同流態(tài)采用不同的計算公式。當(dāng)出流為自由出流時，其計算公式為：

當(dāng)出流為淹沒出流時其計算公式為：

式中：B 為堰寬（m）；Z1為堰上水位（m）；Z2為堰下水位（m）；H0為堰上總水頭（m），H0=Z1-Zd；hs為堰頂水深（m），hs=Z2-Zd；g為重力加速度（m/s2）；m為自由出流系數(shù)，一般取m=0.325～0.385；?m為淹沒出流系數(shù)，一般取?m＝1.0～1.18；Zd為堰頂高程（m）。

對式（8）—（9）離散后，自由出流流量Q=δZ1Z1+βZ1，淹沒出流流量Q=δZ2(Z1-Z2)。其中，δZ1，δZ2，1Zβ為與Z1，Z2有關(guān)的系數(shù)，一般常采用時段初水位來計算，有時為了提高計算精度，可采用迭代法計算。

2.2 模型建立

模型的建立主要包括地形及網(wǎng)格剖分、人工控制工程概化、控制條件設(shè)置、邊界條件確立等。

2.2.1 模擬范圍

本次涉及的模擬范圍主要包括漳河、衛(wèi)河出山口以下的河道及蓄滯洪區(qū)，即漳河（京廣鐵路橋下—徐萬倉）、衛(wèi)河（淇門—徐萬倉）、衛(wèi)運(yùn)河（徐萬倉—四女寺）、漳衛(wèi)新河（四女寺—入海口）4 條主要河道，還包括淇河、湯河、安陽河3 條衛(wèi)河的主要支流；大名泛區(qū)、良相坡、長虹渠（含柳圍坡）、白寺坡、共渠西行洪區(qū)、小灘坡、任固坡7個蓄滯洪區(qū)。計算范圍概化示意，如圖2所示。

圖2 洪水演進(jìn)模型建模范圍示意

2.2.2 數(shù)字地形與網(wǎng)格剖分

一維模型的地形數(shù)據(jù)是以1∶2 000 河道實(shí)測大斷面資料為基礎(chǔ)，將兩堤坐標(biāo)、起點(diǎn)距、高程導(dǎo)入斷面編輯器，參與模型計算。斷面編輯器中可以對河道進(jìn)行主槽、灘地分界，設(shè)置糙率。河道斷面編輯器界面，如圖3所示。

圖3 河道斷面編輯器界面

良相坡、長虹渠等蓄滯洪區(qū)地形數(shù)據(jù)是以1∶10 000地形圖提取的散點(diǎn)高程為基礎(chǔ)，以蓄滯洪區(qū)圍堤為固定邊界，將x，y，z文件導(dǎo)入軟件，生成不規(guī)則三角網(wǎng)，得到模擬區(qū)域的地形文件，結(jié)果如圖4所示。

本次二維計算采用單元中心的有限體積法進(jìn)行數(shù)值模擬，計算區(qū)域采用直角坐標(biāo)系下非均勻矩形網(wǎng)格進(jìn)行地形的剖分。網(wǎng)格間距根據(jù)模擬區(qū)域面積大小、地形變化情況、模擬精度等因素綜合分析后設(shè)置，當(dāng)關(guān)注某一區(qū)域的水文要素時可以設(shè)置網(wǎng)格加密，分為單獨(dú)橫向加密和單獨(dú)縱向加密以及局部單元格加密。模型模擬區(qū)域網(wǎng)格剖分結(jié)果，如圖5所示。

小灘坡和任固坡無近期地形資料，利用水位－蓄量關(guān)系建立零維模型，以水位節(jié)點(diǎn)的形式參與計算，忽略其內(nèi)部流場。零維要素編輯器界面，如圖6所示。

圖4 蓄滯洪區(qū)（以大名泛區(qū)為例）地形高程渲染

圖5 蓄滯洪區(qū)（以大名泛區(qū)為例）地形剖分

圖6 零維要素編輯器界面

2.2.3 工程模擬

工程設(shè)施主要是指流域中控制水流運(yùn)動的堰、水閘及蓄滯洪區(qū)分洪口門等。這些設(shè)施的過流滿足水力學(xué)上的堰流公式，其原理與模型耦合一致。模型中將工程概化為“聯(lián)系”。河道上相鄰的斷面之間、不同河道的斷面之間、河道斷面與二維網(wǎng)格之間、不同的二維網(wǎng)格之間均可根據(jù)水流方向建立“聯(lián)系”，分別模擬河道攔河閘、穿堤涵洞、分洪（閘）口門等水利工程。模擬這些工程時，首先根據(jù)本身實(shí)際坐標(biāo)找到在模型中的位置，然后將閘底板高程、閘孔寬度、閘孔高度輸入模型，如果是寬頂堰型，則只需輸入堰頂高程和寬度。為了提高計算穩(wěn)定性，根據(jù)剖分網(wǎng)格尺寸大小考慮將較寬的口門分為若干“聯(lián)系”進(jìn)行模擬，如網(wǎng)格邊長200 m、口門寬度700 m可用3個寬度為200 m、1個寬度為100 m的“聯(lián)系”組合模擬。

2.2.4 控制條件設(shè)置

水閘、口門等工程的具體控制調(diào)度方法在模型中是以控制條件的形式實(shí)現(xiàn)的。本次計算的控制條件依據(jù)《漳衛(wèi)河洪水調(diào)度方案》編寫。臨時破堤型口門，當(dāng)水位或者流量達(dá)到調(diào)度條件時被扒開。由于破堤口門被炸開分洪時并不是瞬間達(dá)到最大寬度，模型中考慮了口門寬度發(fā)展速度控制，即模擬口門自扒開至最大口門寬度這個過程隨時間變化的實(shí)際情況。

2.2.5 邊界條件確立

模型中共有5個洪水入口，分別是漳河入流、良相坡入流、淇河入流、湯河入流和安陽河入流，其設(shè)計洪水流量過程作為模型的上邊界條件；1 個洪水出口，即漳衛(wèi)新河河口，其設(shè)計潮位過程作為模型的下邊界條件。模型邊界條件示意，如圖7所示。

圖7 模型邊界條件示意

2.2.6 露灘問題處理

對于水流還沒有流到的露灘區(qū)域，需進(jìn)行特殊處理。目前，比較合理的方法是采用動邊界方法確定出計算域中有水和無水區(qū)域的界限進(jìn)行模擬，但是在程序編譯處理過程中比較復(fù)雜困難。本模型作近似簡化處理，通過“凍結(jié)法”和最小水深假設(shè)的方法把露灘問題處理為固定邊界。

2.3 參數(shù)率定

計算模型的主要參數(shù)包括用于計算河道流量模數(shù)（K）的曼寧糙率系數(shù)（n）、計算分洪口門流量的堰流系數(shù)（m）。其中，分洪口門一般按照自由出流考慮，堰流系數(shù)一般取0.325～0.385，考慮到分洪口門沒有控制閘門，計算中一般按低限考慮，可不設(shè)定。因此，模型計算時需要設(shè)置的主要參數(shù)是糙率。

天然河道的糙率一般與斷面形狀、河床特性等因素有關(guān)，變化范圍較大。由于沒有完整的實(shí)測資料來率定模型糙率參數(shù)，目前采用漳衛(wèi)河系防洪規(guī)劃成果中的設(shè)計糙率。模型糙率取值，見表1。

表1 模型糙率（n）取值

3 模擬計算

運(yùn)行模型之前需要給定初始條件，對于流域模型而言，上游至下游主要洪水流路必須通暢，有條件控制的分洪口門及蓄滯洪區(qū)不參加初始條件計算。給定一個高于最高地面高程的水位值進(jìn)行初始條件計算，當(dāng)整個模型的水位變幅為0 即初始條件計算達(dá)到恒定時開始正式計算。

選取50年一遇設(shè)計洪水過程，利用之前制作的地形及模型參數(shù)，進(jìn)行漳衛(wèi)河系洪水演進(jìn)模擬，不僅可以獲得蓄滯洪區(qū)洪水的淹沒要素，還可以輸出河道水面線、控制站水位（流量）過程、蓄滯洪區(qū)控制點(diǎn)水位、蓄滯洪區(qū)最大入流、出流等。主要計算結(jié)果，見表2。

表2 50年一遇洪水模擬主要成果

洪水演進(jìn)過程支持實(shí)時查詢，即在地圖上點(diǎn)擊任意位置，就可獲得該時刻該位置的水位、水深、流速、水量、地形高程等信息。以洪水入流后160 h為例，蓄滯洪區(qū)對應(yīng)的淹沒范圍、河道水面線、流量過程、水位過程，如圖8—10所示。

圖8 長虹渠淹沒范圍

圖9 衛(wèi)運(yùn)河水面線

圖10 四女寺樞紐流量、水位過程（當(dāng)發(fā)生50年一遇設(shè)計洪水時）

經(jīng)過對比，本次研究洪水演進(jìn)的成果與防洪規(guī)劃成果有所不同，主要不同點(diǎn)見表3—4。

表3 河道計算成果對比水位：m；流量：m3/s

表4 蓄滯洪區(qū)計算成果對比水位：m；流量：m3/s

本次研究洪水演進(jìn)成果與防洪規(guī)劃成果之所以有所不同，是由于兩次洪水演進(jìn)采用的口門規(guī)模、運(yùn)用規(guī)則不同。本次研究按照現(xiàn)行調(diào)度方案編寫控制條件。如，規(guī)劃中王灣口門寬度采用800 m、淇門160 m、圈里口門200 m，本次研究分別采用370、380、140 m，底高也不盡相同（見表5）；本次研究長虹渠僅通過淇門入流，而規(guī)劃先運(yùn)用宋村進(jìn)行分洪；本次研究白寺坡王灣入流大于規(guī)劃，是由于共渠西上部李橋水位高于63.471 m 時即通過邢固扒口分洪入白寺坡，白寺坡最大入流為邢固、王灣合計最大，而規(guī)劃中邢固口門未啟用?？傮w上講，本次研究洪水演進(jìn)的結(jié)果較為合理，與規(guī)劃報告的趨勢保持一致。

表5 規(guī)劃與本次研究采用口門規(guī)模及控制條件對比

4 結(jié)論

筆者以漳衛(wèi)河系洪水演進(jìn)為例研究了數(shù)字流域系統(tǒng)模擬一維、二維、零維水力學(xué)的建模過程，得到以下主要結(jié)論。

（1）數(shù)字流域系統(tǒng)中的水動力學(xué)模塊能夠較好地應(yīng)用于大型流域河網(wǎng)一維、蓄滯洪區(qū)二維及零維洪水演進(jìn)數(shù)值模擬。本次研究以漳衛(wèi)河山區(qū)以下發(fā)生50年一遇洪水時按照《漳衛(wèi)河洪水調(diào)度方案》進(jìn)行調(diào)度為模擬對象，與以往成果進(jìn)行對比驗證，結(jié)果顯示模型的模擬精度滿足工作需要，可用于平原河網(wǎng)地區(qū)、大型流域的洪水演進(jìn)分析。

（2）二維模型采用隱格式有限體積法計算，在滿足模型穩(wěn)定性和精度要求的條件下計算步長可達(dá)600 s以上，計算效率較高。較快的計算速度使模型能夠應(yīng)用于洪水預(yù)報、調(diào)度方案調(diào)整、大型流域防洪規(guī)劃等領(lǐng)域。

（3）由于海河流域缺乏完整的實(shí)測洪水資料，本次計算僅使用以往規(guī)劃成果所采用的模型參數(shù)進(jìn)行模型調(diào)試，再與規(guī)劃成果進(jìn)行對比，分析各項條件的異同，判斷計算結(jié)果及其趨勢是否合理，因此對模型驗證結(jié)果的可靠性支撐略顯欠缺。未來，隨著流域洪水記錄的積累，可對模型進(jìn)行更為全面的驗證，以進(jìn)一步提高模擬精度。

[1]田景環(huán)，李芳芳.河道洪水演進(jìn)淺析及一維數(shù)學(xué)模型的建立[J].中國水運(yùn)，2007，（8）：94-95.

[2]魏凱，梁忠民，王軍.基于MIKE21 的濛洼蓄滯洪區(qū)洪水演算模擬[J].南水北調(diào)與水利科技，2013，（6）：16-19.

[3]王船海，李光熾.行蓄洪區(qū)型流域洪水模擬[J].成都科技大學(xué)學(xué)報，1995，（2）：6-14.

[4]周潮洪，杜津媛，常守權(quán).蓄滯洪區(qū)二維洪水演進(jìn)系統(tǒng)的建立[J].海河水利，2002，（6）：33-34.

[5]海河水利委員會.漳衛(wèi)河系防洪規(guī)劃[R].天津：海河水利委員會，2008.