胡文才

（沂沭泗水利管理局，江蘇 徐州 221009）

0 引言

近年來(lái)由于水利工程的大量修建，改變了河流與湖泊特性，在做湖泊洪水預(yù)報(bào)調(diào)度中發(fā)現(xiàn)，單純采用水文學(xué)的方法做洪水預(yù)報(bào)調(diào)度，其精度完全不能滿足實(shí)際防洪調(diào)度需求，因此不斷有專家學(xué)者將水力學(xué)和水文學(xué)的方法結(jié)合起來(lái)研制洪水預(yù)報(bào)調(diào)度系統(tǒng)。20 世紀(jì) 80 年代初意大利的 E.Todini 教授與我國(guó)合作進(jìn)行了淮河王家壩至正陽(yáng)關(guān)河道洪水演進(jìn)與滯洪區(qū)洪水調(diào)度，采用水力學(xué)與水文學(xué)集合的方法；2003 年河海大學(xué)李致家教授等在南四湖流域分別采用一維和二維水力學(xué)方法進(jìn)行應(yīng)用研究，取得了一定的成果[1-4]。

本次在駱馬湖洪水預(yù)報(bào)調(diào)度中區(qū)間產(chǎn)流采用經(jīng)驗(yàn)法計(jì)算至入湖口，南四湖和沂河來(lái)水采用馬斯京根法計(jì)算至駱馬湖入湖口。沂河港上、中運(yùn)河運(yùn)河站和房亭河土山站的入流分別作為網(wǎng)格入流處理。采用水量平衡的方法進(jìn)行駱馬湖洪水調(diào)度，以確定各個(gè)閘的放水流量過(guò)程，作為水力學(xué)的下邊界條件，然后根據(jù)二維水力學(xué)方法演算駱馬湖內(nèi)的洪水過(guò)程以確定各個(gè)網(wǎng)格點(diǎn)的流量和水位。整個(gè)駱馬湖的初始水位場(chǎng)由駱馬湖內(nèi) 9 個(gè)水位站的實(shí)時(shí)水位進(jìn)行資料插值獲得。

1 基本原理

駱馬湖二維水動(dòng)力學(xué)模型，采用平面二維水流數(shù)學(xué)模型進(jìn)行模擬計(jì)算，為反映駱馬湖內(nèi)外的河道及岸線邊界，采用邊界擬合網(wǎng)格，建立正交曲線坐標(biāo)系下的二維水流運(yùn)動(dòng)數(shù)學(xué)模型，運(yùn)用有限差分法離散水流運(yùn)動(dòng)基本方程組，通過(guò)交替隱式解法（ADI）法求解離散方程[5]，得到整個(gè)計(jì)算區(qū)域上的水位和流速。

1.1 二維水動(dòng)力學(xué)原理

1.1.1 水流運(yùn)動(dòng)基本方程

假定壓力沿水深按靜水壓力分布，在正交曲線坐標(biāo)系 ξ-η下的二維不恒定流連續(xù)性方程為：

式中：u，v分別為ξ和η方向的流速分量；Z 為水位；H 為水深，H = Z + h，h 為基面以下水深；q 為包括雨量等旁側(cè)流量；g 為重力加速度；C 為謝才系數(shù)，C =H1/6，n 為糙率；εξ、εη為紊動(dòng)粘性系數(shù)；ρ、ρα分別為水和空氣的密度；vf為水面上 10 m 處風(fēng)速；CD為風(fēng)剪切應(yīng)力系數(shù)，CD= (1.1 + 0.0536vf)×10-3。

1.1.2 ADI 法

采用有限差分法離散方程，通過(guò) ADI 法求解方程。

圖1 為正交曲線坐標(biāo)下的流速和水深點(diǎn)的交錯(cuò)布置，對(duì)以水位、水深點(diǎn)為中心的陰影計(jì)算單元，流速 u 布置在 Z 方向的邊上（圖中橫箭頭），流速v布置在 h 方向的邊上（圖中豎箭頭），圓點(diǎn)處為水位、水深和糙率等的位置，可以避免將變量布置在同一位置而引起的計(jì)算波動(dòng)。

圖1 交錯(cuò)網(wǎng)格布置

根據(jù) ADI 法，在 n ～（n + 1/2）時(shí)間步長(zhǎng)內(nèi)對(duì)式 (1)、(2) 聯(lián)立求解，變量 v以 n 時(shí)刻值代入，其差分方程為：

式中：i，j 為空間變量；A，B，C，D 是對(duì)原來(lái)連續(xù)性和動(dòng)量方程經(jīng)過(guò)這個(gè)交替隱式解法差分后的系數(shù)，無(wú)具體物理意義。

數(shù)學(xué)計(jì)算流程如圖 2 所示。

1.2 馬斯京根法

馬斯京根法是一個(gè)流行的集總式河道演算方法，它假設(shè)了 1 個(gè)可變化的流量～蓄量方程：

取水量平衡方程和槽蓄方程差分解，可得流量演算方程：

圖2 數(shù)學(xué)模型計(jì)算流程

式中：W 為槽蓄量；k 為蓄量常數(shù)；x 為求槽蓄量 W時(shí)，河段楔蓄量所占的權(quán)重；C0、C1、C2為河道演算系數(shù)。

1.3 降雨徑流相關(guān)法

區(qū)間產(chǎn)流計(jì)算采用經(jīng)驗(yàn)方法——降雨徑流相關(guān)法，即 P + Pa～R。

2 具體應(yīng)用

2.1 流域概況及水文特征

駱馬湖位于沂河末端，中運(yùn)河?xùn)|側(cè)，是以防洪、蓄水為主，結(jié)合航運(yùn)、發(fā)電、水產(chǎn)養(yǎng)殖等綜合利用的多功能湖泊。駱馬湖承接沂河干流、南四湖、邳蒼地區(qū) 5.1 萬(wàn) km2面積的來(lái)水，調(diào)蓄后主要由新沂河排入黃海。入湖主要河道有沂河、中運(yùn)河，出湖主要河道有新沂河、中運(yùn)河和六塘河。

駱馬湖南北長(zhǎng) 20 km，東西寬 16 km，周長(zhǎng) 70 km。一般湖底高程 20.0 m，最低 19.0 m。正常蓄水位 23.0 m 時(shí)，湖面面積 375 km，容積 9.0 億 m3；退守宿遷大控制后，當(dāng)達(dá)到設(shè)計(jì)洪水位 25.0 m 時(shí)，包括駱馬湖與中運(yùn)河之間三角地帶面積 18 km2，湖面面積為 450 km2，容積 15 億 m3；校核洪水位 26.0 m 時(shí)，總庫(kù)容 19.0 億 m3[6]。

2.2 入湖水量計(jì)算

駱馬湖入湖水量主要由以下 3 部分組成：1）沂河來(lái)水，沂河臨沂站來(lái)水除去彭道口閘和江風(fēng)口閘分泄流量，剩余水量全部進(jìn)入駱馬湖；2）運(yùn)河來(lái)水，中運(yùn)河來(lái)水包括南四湖下泄水量和運(yùn)駱邳蒼區(qū)間產(chǎn)水；3）房亭河來(lái)水，房亭河劉集閘站以上流域面積 933 km2的區(qū)間產(chǎn)水。各控制斷面位置如圖 3 所示。

圖3 駱馬湖入湖水量控制斷面示意圖

圖中 1、2、3 分別為 3 個(gè)入湖控制斷面，其中斷面 1 流量過(guò)程由臨沂站過(guò)程采用馬斯京根法演算得出，斷面 2、3 流量過(guò)程分別由 P + Pa～R 計(jì)算和馬斯京根法演算結(jié)果疊加得出，馬斯京根法河段演算系數(shù)如表1 所示。

表1 馬斯京根法河段演算系數(shù)表

2.3 模型應(yīng)用范圍

根據(jù)地形和工程的情況，以駱馬湖湖區(qū)范圍建立模型，建立如圖 4 所示的駱馬湖二維水動(dòng)力模型網(wǎng)格的概化地形圖，共有 30×38 個(gè)，網(wǎng)格尺寸190～1200 m。駱馬湖二維水動(dòng)力模型湖區(qū)地形采用 2001年1:10000 地形圖，高程采用廢黃河零點(diǎn)基面。模型的邊界條件為新沂河嶂山閘、六塘河洋河灘閘、中運(yùn)河皂河閘、中運(yùn)河運(yùn)河鎮(zhèn)、沂河入湖口的流量過(guò)程。湖區(qū)庫(kù)容的計(jì)算范圍在駱馬湖湖區(qū)內(nèi)。

2.4 初始條件

初始條件如下：

圖4 駱馬湖二維水動(dòng)力模型概化地形

式中：u0、v0分別為初始流速在 ξ 和 η 上的分量，計(jì)算時(shí)取流速 u0= 0 和 v0= 0；Z0為初始水位，為一給定值，由駱馬湖現(xiàn)有的 9 個(gè)水位站（皂河閘、洋河灘閘、嶂山閘、苗圩、窯灣、袁場(chǎng)、張宅、新店、曉店）預(yù)報(bào)開始時(shí)刻的實(shí)測(cè)資料插值獲得整個(gè)駱馬湖的初始水位場(chǎng)。模型經(jīng)過(guò)一定時(shí)間的運(yùn)行，初始條件的影響將會(huì)消除。

2.5 動(dòng)邊界處理

在河湖地區(qū)，灘地和沙洲隨著水位的漲落淹沒和露出，在數(shù)模計(jì)算中要進(jìn)行動(dòng)邊界處理，根據(jù)水位的變化連續(xù)不斷地修正邊界，模擬灘地水邊線的移動(dòng)。本文采用富裕水深法，即在計(jì)算中判斷每個(gè)網(wǎng)格的水深，當(dāng)網(wǎng)格水深大于某一給定的小水深時(shí)，將單元開放，作為計(jì)算水域；反之，將單元關(guān)閉，置流速于零。

2.6 邊界條件

模型開邊界一般用水位或流量過(guò)程控制。在駱馬湖二維水動(dòng)力數(shù)學(xué)模型中，出、入湖邊界條件均用流量控制，出、入湖方向?yàn)?3 進(jìn) 4 出，數(shù)學(xué)模型邊界條件如圖 5 所示。模型的邊界條件為新沂河嶂山閘、六塘河洋河灘閘、中運(yùn)河皂河閘、中運(yùn)河運(yùn)河鎮(zhèn)、沂河入湖口、老沂河分洪道入湖口的流量過(guò)程。閉邊界采用不可入條件，即取邊界的外法線方向流速 vn= 0。

2.7 計(jì)算調(diào)度

駱馬湖湖區(qū)洪水演進(jìn)預(yù)報(bào)采用 96 h的時(shí)段長(zhǎng)度，與水文學(xué)模型相結(jié)合，根據(jù)實(shí)測(cè)降雨過(guò)程作洪水預(yù)報(bào)，通過(guò)水文學(xué)模型計(jì)算駱馬湖流域的產(chǎn)匯流狀況，得出駱馬湖湖區(qū)的入流邊界條件，分別是沂河、中運(yùn)河、駱馬湖區(qū)間等入流過(guò)程及駱馬湖湖面產(chǎn)流過(guò)程。然后再根據(jù)不同的洪水調(diào)度方案，得出駱馬湖湖區(qū)的出流邊界條件，最后模擬預(yù)測(cè) 96 h 內(nèi)駱馬湖湖區(qū)水位、流量的變化狀況。

計(jì)算時(shí)首先根據(jù)洋河灘閘的 t 時(shí)刻的水位 Zt，利用駱馬湖水位-庫(kù)容曲線，求出 t 時(shí)刻的駱馬湖的庫(kù)容 Vt；再利用駱馬湖預(yù)報(bào)模型中產(chǎn)匯流及湖區(qū)降雨量得出的 △t 時(shí)間內(nèi)的入湖量 △V，求出 (t + △t) 時(shí)刻駱馬湖的庫(kù)容 Vt+△t，Vt+△t=Vt+△V；

最后通過(guò)駱馬湖水位-庫(kù)容曲線，預(yù)估出 (t + △t)時(shí)刻駱馬湖洋河灘閘的水位 Zt+△t。

圖5 駱馬湖數(shù)學(xué)模型邊界條件

根據(jù)上述駱馬湖洪水調(diào)度原則和方案，結(jié)合預(yù)估出的 (t + △t) 時(shí)刻駱馬湖洋河灘閘的水位 Zt+△t，確定具體的駱馬湖湖區(qū)洪水的調(diào)度規(guī)則，具體調(diào)度規(guī)則如表2 所示，再通過(guò)駱馬湖二維水動(dòng)力數(shù)學(xué)模型計(jì)算，預(yù)測(cè)當(dāng)前駱馬湖各種調(diào)度規(guī)則下，湖區(qū)的洪水水位及流量的變化過(guò)程。

表2 駱馬湖湖區(qū)水動(dòng)力洪水調(diào)度規(guī)則

2.8 實(shí)例應(yīng)用

選擇 2006、2008 和 2009 年，共 5 場(chǎng)洪水進(jìn)行模擬預(yù)報(bào)調(diào)度，預(yù)報(bào)調(diào)度結(jié)果如表3 所示。由于篇幅問題示例圖僅列 2009年2 場(chǎng)洪水預(yù)報(bào)調(diào)度過(guò)程圖，具體如圖 6、7 所示。

表3 實(shí)際應(yīng)用成果統(tǒng)計(jì)表

圖6 2009-07-17 洪水預(yù)報(bào)調(diào)度結(jié)果

圖7 2009-08-17 洪水預(yù)報(bào)調(diào)度結(jié)果

模擬預(yù)報(bào)結(jié)果顯示，水位預(yù)報(bào)誤差最大為 0.13 m。預(yù)報(bào)最高水位出現(xiàn)時(shí)間，與實(shí)際最高水位出現(xiàn)時(shí)間相比，有的提前有的滯后，最大相差 1 d。出現(xiàn)這一結(jié)果主要有 3 個(gè)方面的原因：1）入湖過(guò)程滯后，模型計(jì)算港上站至駱馬湖洪水傳播時(shí)間為 6 h，但是實(shí)際傳播時(shí)間為 12 h 左右；2）實(shí)際與計(jì)算調(diào)度過(guò)程不可能完全一致；3）預(yù)報(bào)入湖過(guò)程存在誤差，入湖過(guò)程預(yù)報(bào)采用經(jīng)驗(yàn)方法，其預(yù)報(bào)結(jié)果與實(shí)際入湖過(guò)程存在一定的誤差。以 2009年7月17日洪水為例，預(yù)報(bào)入湖時(shí)間為 22日 5:00，但是實(shí)際入湖時(shí)間為 22日 11:00 左右，滯后 6 h；預(yù)報(bào)入湖洪峰流量為 5835 m3/s，但實(shí)際最大入湖洪峰流量為 4600 m3/s左右，預(yù)報(bào)入湖洪峰流量比實(shí)際入湖洪峰大1235 m3/s。由于實(shí)際入湖過(guò)程比預(yù)報(bào)過(guò)程滯后，導(dǎo)致計(jì)算結(jié)果最高水位出現(xiàn)時(shí)間比實(shí)際時(shí)間提前，預(yù)報(bào)最大入湖洪峰流量大于實(shí)際入湖洪峰流量，導(dǎo)致預(yù)報(bào)最高水位高于實(shí)際最高水位。

2.9 優(yōu)缺點(diǎn)

該系統(tǒng)采用水力學(xué)與水文學(xué)結(jié)合的方法進(jìn)行預(yù)報(bào)調(diào)度計(jì)算，湖內(nèi)洪水演進(jìn)采用水力學(xué)的方法，考慮了動(dòng)庫(kù)容的問題，解決了傳統(tǒng)水文學(xué)方法計(jì)算時(shí)由于動(dòng)庫(kù)容帶來(lái)的誤差問題，使預(yù)報(bào)精度提高了。該系統(tǒng)水力學(xué)計(jì)算采用 Fortran 語(yǔ)言編程，計(jì)算速度提高了，目前做 1 次預(yù)報(bào)時(shí)間大約僅需要 2 min，與傳統(tǒng)預(yù)報(bào)系統(tǒng)相比時(shí)間縮短了。

其他相似湖泊要應(yīng)用該系統(tǒng)，需要考慮入湖流量和時(shí)間誤差，以及實(shí)際出湖過(guò)程與預(yù)報(bào)過(guò)程不一致的情況。如果考慮上述幾個(gè)因素帶來(lái)的影響，進(jìn)行預(yù)報(bào)結(jié)果修正，將會(huì)帶 來(lái)更好的預(yù)報(bào)成果。

3 結(jié)語(yǔ)

本文將二維水流運(yùn)動(dòng)數(shù)學(xué)模型，有限差分法離散水流運(yùn)動(dòng)基本方程組，ADI 法求解離散方程的方法和水文學(xué)的方法結(jié)合應(yīng)用于駱馬湖洪水預(yù)報(bào)調(diào)度中；地形資料采用了 2001年1:10000 地形圖，與實(shí)際地形基本吻合；采用 9 個(gè)水位控制站的實(shí)時(shí)水位，通過(guò)插值獲取計(jì)算初始時(shí)刻駱馬湖的水位場(chǎng)。在實(shí)際應(yīng)用中發(fā)現(xiàn)二維水力學(xué)進(jìn)行湖泊水力學(xué)演算，如果地形資料準(zhǔn)確，湖泊內(nèi)水位觀測(cè)站點(diǎn)代表性較好，采用該方法能夠較好地模擬各個(gè)測(cè)點(diǎn)的水位。在實(shí)時(shí)洪水預(yù)報(bào)和調(diào)度中，采用水文學(xué)和水力學(xué)結(jié)合的方法進(jìn)行洪水預(yù)報(bào)調(diào)度，能取得令人滿意的預(yù)報(bào)成果。

[1]李致家，董增川，梁忠民，等. 大流域洪水預(yù)報(bào)與洪水調(diào)度管理方法研究[J]. 河海大學(xué)學(xué)報(bào)，2004, 30 (1): 12-15.

[2]李致家. 通用一維河網(wǎng)不恒定流軟件的研究[J]. 水利學(xué)報(bào)，1998 (8): 14-18.

[3]汪德灌. 計(jì)算水力學(xué)-理論與應(yīng)用[M]. 南京：河海大學(xué)出版社，1989: 64-101.

[4]李致家，包紅軍，孔祥光，等. 水文學(xué)與水力學(xué)相結(jié)合的南四湖洪水預(yù)報(bào)模型[J]. 湖泊科學(xué)，2005, 17: 299-304.

[5]徐峰俊，劉俊勇. 伶仃洋海區(qū)二維不平衡非均勻輸沙數(shù)學(xué)模型[J]. 水力學(xué)報(bào)，2003 (7): 16-23.

[6]鄭大鵬. 沂沭泗防汛手冊(cè)[M]. 徐州：中國(guó)礦業(yè)大學(xué)出版社，2003: 95-97.