基于二維淺水方程的城市地面洪水演進數值模擬研究

冀永鵬，張洪興，王運濤，李 晉,張明亮

(1.大連海洋大學 海洋科技與環境學院，遼寧 大連 116023；2.大連理工大學 港口、海岸與近海工程國家重點實驗室，遼寧 大連 116025；3.大連理工大學 水資源與防洪研究所，遼寧大連 116025；4.盤錦鴛鴦溝國家級海洋公園管理辦公室，遼寧 盤錦 124010)

1 研究背景

近年來隨著城市化速度的加快，城市及周邊區域自然狀態下的森林植被、湖泊沼澤等區域逐步被城市道路、建筑等基礎設施所取代，其下墊面的不透水面積不斷增加。城區下墊面性質和功能的變化致使城市的水循環途徑發生了改變，某種程度上增加了城市雨洪的地表徑流量[1]。再加上城市建設中排水設施存在諸多隱患和問題，導致城市雨洪和內澇災害頻繁發生，造成城區居民人員傷亡以及巨大經濟損失?；诖耍P于城市內澇成因、地表洪水演進預報及內澇防治等關鍵問題的研究受到越來越多學者的關注。作為研究城市內澇問題的重要技術手段，城市雨洪數值模擬技術是城市防洪減災的關鍵技術之一，也是當前研究的熱點問題。目前，國內外已經發展了很多城市雨洪模型，城市雨洪模型包括地面雨洪產流、地面雨洪匯流以及地面雨洪管流的計算，其中產、匯流計算是模擬城市雨洪過程的關鍵。美國環保署研發的SWMM模型[2]、美國工程師協會水文工程中心的Storm模型以及英國Wallingford水力學研究機構開發的 InfoWorks CS軟件模型等[3]是根據水文學原理建立的城市雨洪模型。還有很多學者依據水動力學原理建立了城市雨洪模型，一維和二維水動力模型是應用較為廣泛的模型。例如丹麥水資源及水環境(DHI)研究的MIKE FLOOD模型[4]；Cea等[5]在二維淺水水動力方程的基礎上提出了解耦離散模型(DHD)，用于求解陸地徑流和降雨徑流的轉換，并通過一些實驗室案例對模型進行了驗證；余江順等[6]應用以水動力學為理論基礎的水利軟件(HydroInfo)對成都繞城公路區域雨洪演進進行了模擬；Xia等[7]基于淺水動力方程提出了一種新的一階Godunov型有限體積法模型 (SRM)，使用HLLC格式計算界面通量，對城市雨洪的理想案例以及實驗室案例進行了驗證。由于城市地形變化大，水動力學模型的計算比水文學模型繁雜，且初始值以及邊界條件的設置較為復雜，所以目前對水動力學模型的應用相對較少，但由于水動力學模型的模擬準確度更高，越來越多的學者開始關注水動力學模型的應用。也有學者建立水文學與水動力學耦合模型對城市洪水進行模擬[8-10]。另外，相比國外模型，國內模型的二次開發能力弱，模型功能也較為單一，應用范圍較為狹窄，可視化水平也有待提高[11]，并且大多模型在城市雨洪模擬中并未考慮到對干濕邊界水流交替進行特殊的處理。

本文基于平面二維淺水方程，在非結構網格下采用Godunov型有限體積法建立二維水動力模型，以Roe格式的近似黎曼解計算三角網格的界面通量，并且針對復雜地形進一步對干濕邊界以及負水深問題進行了處理，利用第一個算例對模型的動邊界處理進行了驗證，利用第二個算例對模型的靜水和諧性及穩定性進行了驗證，并將模型應用于城市街道及附近區域的地表雨洪演進模擬，驗證模型在城市復雜地形上進行洪水預測的準確性，最后使用模型對渾太胡同內的洪水演進進行了模擬預測。

2 水動力模型

2.1 平面二維淺水模型

平面二維淺水模型是對Navier-Stokes方程進行垂向平均近似，并遵循淺水Boussinesq假設所得到的控制方程。本文忽略了二維淺水方程中的風應力、科氏力和二階擴散項，二維淺水控制方程最終可以表達為：

(1)

式中：U為守恒性變量；F和G分別為x方向和y方向對流通量變量；S為源項。

(2)

式中：h為水深，m；u和v分別為x和y方向的流速，m/s；Ri為降雨強度，mm/h；g為重力加速度，m/s2；η為水位，m；τbx和τby分別為x和y方向的摩阻項，τbx=n2u(u2+v2)1/2h-4/3，τby=n2v(u2+v2)1/2h-4/3；n為曼寧系數。

2.2 數值離散

本文采用有限體積法，以非結構三角網格為控制單元對控制方程進行離散，運用 Green 公式對控制體進行積分整理可得：

(3)

本文模型為顯格式計算模式，下一時刻的離散值由上一時刻求得，其表達式為：

(4)

2.3 邊界通量計算

由于對方程進行深度平均積分，物理變量在每個三角網格單元中均為常數，在整個求解域內構成一系列分片函數，非結構網格在每個控制單元邊界有間斷現象，形成了局部的Riemann問題。本文采用Roe格式求解Riemann問題，數值通量表達式為：

(5)

2.4 變時間步長

本文控制方程求解采用顯格式計算，因此用CFL條件來控制其穩定性。在非結構三角網格下，應用CFL條件來確定時間步長[12]，其計算方法為：

(6)

式中：Cr為庫朗數，取值范圍是0

3 動邊界及負水深處理

3.1 動邊界處理

由于城市下墊面并不平坦，降水后城市下墊面會出現局部積水現象，從而導致非結構網格單元出現干濕交替或者部分淹沒的情況，形成干濕界面即動邊界。

本文針對復雜地形引入了干濕界面處理技術[13]，設置了最小限制水深h0，水深大于h0即為濕邊界，小于h0即為干邊界。4種干濕邊界分類情況如圖1所示。圖1 (a)中為濕邊界，邊界兩側水深均大于h0；圖1(b)、1(c)均為半濕邊界，其中圖1(b)中邊界有通量，圖1 (c)中邊界無通量；圖1(d)為干邊界，邊界兩側水深均小于h0。

圖1 4種干濕邊界分類示意圖

根據干濕邊界的判斷，可進一步確定網格單元的干濕情況：網格邊界均為濕邊界或有通量的半濕邊界則為濕網格；網格邊界均為干邊界或者無通量的半濕邊界則為干網格；其余均為半濕網格。在無通量半濕邊界，由于ηL與ηR并不相等，形成的非物理流量會產生質量不守恒問題，所以將網格單元的ηR降到與左側網格單元的ηL相等的位置(圖1(c)虛線所示)，但這種處理仍會導致偽流速問題的出現。本文針對無通量半濕邊界處的偽流速問題做了進一步處理，由于網格單元的水位梯度最終由網格節點水位決定，因此將無通量半濕邊界處節點的水位降至與左側網格單元的ηL水位相等，保證了質量守恒的同時也避免了偽流速的形成，提高了模型的靜水和諧性和穩定性。

3.2 負水深處理

在基于平面二維淺水方程進行模擬的過程中，邊界通量的計算由時間步長、單寬通量以及邊界長度決定。若時間步長設置過小，則會極大地降低計算效率，而在小水深情況下，時間步長設置較大就會使得網格邊界通量大于上游網格的水量，導致邊界上游網格水深成為負值，形成負水深，最終造成質量不守恒。

針對負水深問題，本文將周圍網格單元最小限制水深h0以上的水量填補入負水深網格單元中，直至負水深成為最小限制水深，從而保持質量守恒，也能保證周圍網格不會因水量填補而再次出現負水深。

4 數值模擬驗證

4.1 拋物型無阻力自由水面模擬

本算例模擬了在無地面阻力情況下，拋物型地形中自由水面的震蕩過程。Thacker[14]給出了該情況下自由水面震蕩情況的解析值。該算例是驗證非平底地形條件下模型處理干濕界面交替時水流運動精度的經典案例。拋物型地形公式為：z(x，y)=h*/a2(x2+y2)，其中h*和a為大于0的常數。當σ

σ)+h*-z(x,y)]

(7)

u=-σwsin(wt)v=-σwcos(wt)

(8)

式中：σ為周期性運動的幅度；w為周期性運動的頻率。在本算例中a=8 025.5 m，h*=10 m，σ=a/10，水流運動周期為3 600 s。該算例中模擬了3h，即3個周期的水流運動；采用了41 276個網格,共有20 889個節點；模擬區域為矩形，x與y方向范圍均為-10 000～10 000 m，四周邊界為固定邊界。模擬的最小限制水深為0.01 m，曼寧系數設置為0。

模擬了y=0沿線在t=5 400 s、t=7 200 s時刻的水位，另外模擬了研究區內坐標為(-7 900，0)、(-5 500，0)兩個測點的水深和流速，模擬的數值解與解析值的比較見圖2～3。由圖2可以看出，每個時刻y=0沿線的水位模擬值與解析值基本吻合，水位震蕩幅度在8～14 m之間,且在無阻力情況下，隨著時間的增加,水位的震蕩幅度并未減??；由圖3可以看出，水深和流速模擬結果良好，坐標為(-5 500，0)的測點處于完全淹沒狀態，坐標為(-7 900，0)的測點位于干濕交替處，水深和流速模擬值在干濕交接界面處沒有明顯波動，準確地模擬了復雜干濕界面的水流運動情況，說明模型對動邊界的處理較為合理。

圖2 不同時刻y=0 沿線水位模擬值與解析值對比

圖3 研究區內兩個測點的水深及流速模擬值與解析值對比

4.2 洪泛漫灘洼地填充模擬

本算例是模擬漫灘洼地區域洪水運動的理想算例，涉及到復雜地形上低動量流的運動情況，主要驗證模型的動邊界處理能力、靜水和諧能力及穩定性[15]。圖4所示為模擬區域等高線、測點以及入流位置分布，模擬區域為2000 m×2000 m的矩形區域，由4×4的洼地矩陣組成，各個洼地形狀一致，地形變化平滑，沿西北向東南對角線地形高程下降2 m。測點均分布于洼地中心，洼地序號(亦為測點編號)如圖4所示。入流分布于西北角向南延伸100 m，入流流量為20 m3/s，從第10 min開始持續到86 min，其他邊界均為固壁邊界，沒有流量流出。該算例模擬區域采用了26 508個網格，初始狀態為干河床，時間步長為變化時間步長，最小限制水深為0.003 m，曼寧系數設定為0.015，共模擬了48 h。

圖5為不同模型計算的部分測點水位隨時間變化的比較。所參照對比的模型包括MIKE FLOOD[4]模型、InfoWorks ICM[10]模型以及ISIS Fast Dynamic[15]模型。MIKE FLOOD模型是基于一維和二維水動力模型進行動態耦合的建模系統。InfoWorks ICM模型是基于有限體積法的二維水動力淺水方程建立的模型，模型具有模擬自然界水流運動以及工程排水系統的能力，并且可以捕捉快速流動的水流變化。ISIS Fast Dynamic模型是利用曼寧均勻流動定律，基于體積擴散法建立的模型，模型可以預測洪水的最后淹沒程度以及類似于二維淺水在低動量流情況下的水流動態，但不適于預測洪水的淹沒速度。從圖5中可以看出，本模型模擬結果與其他模型結果均大致相符，雖然結果有時存在有一定的數值差異，但水位的變化趨勢是基本一致的。當入流流量開始進入時，在入流點附近測點(測點4、7)的水位可以觀測到瞬時水位峰值。入流消失后，附近測點的水位會逐漸下降，直至洼地邊緣最高處洪水無法流出洼地則水位達到平衡。離入流點越遠的測點瞬時水位峰值會逐漸變小甚至消失。測點5、10的水位變化波動較大，比其他模型的水位到達時間要早大約5～9 h，兩個測點的水位相比其他模型也較高。這主要是由于洼地之間區域最小限制水位以上的水深較淺，而限制水位以上的水深有著微小的差異就會導致水位的大小以及水位隨時間的變化表現出顯著的不同。

圖4 漫灘洼地模擬區域等高線、測點及入流位置示意圖

圖5 漫灘洼地區域洪水演進算例不同模型計算的部分測點水位隨時間變化的比較

圖6為第6、12 h的水深淹沒范圍與流場分布。分析圖6可知，水流沿著地勢方向向下流動，在計算區域東西方向地形平坦,水流速度緩慢，南北方向地形有較大起伏，水流速度較快，水流平穩后的洼地四周幾乎沒有流速的分布也說明了模型對偽流速處理的合理性。

綜上所述，該模型對本算例的模擬較為良好，模型的動邊界處理較為合理，靜水和諧性較好，可以應用于模擬和預報城市雨洪形成的地表徑流演進過程。

圖6 漫灘洼地區域第6、12 h的水深淹沒范圍及流場分布

4.3 城市降水及點源徑流洪水演進模擬

本算例的模擬區域位于英國格拉斯哥的科肯齊街及周圍街道，忽略了街道以外的下墊面以及建筑[15]。模擬區范圍面積約為0.4 km×0.96 km，地形及道路分布、點源徑流的入流點、測點分布位置見圖7，地形數據分辨率為0.5 m。城市洪水的來源有兩種：一是降雨徑流，二是地面的點源徑流。模擬降水時間為1～4 min，降雨強度為400 mm/h，模擬入流流量從第23 min開始逐漸增加，在37～39 min達到5 m3/s，隨后流量逐漸減小，在53 min時停止入流。模擬采用了21 261個網格，區域四周邊界均為固壁邊界，模擬初始狀態為干河床，時間步長為變化時間步長，共模擬了5 h的洪水演進過程。道路路面的曼寧系數為0.02，其余下墊面曼寧系數為0.05。圖8(a)顯示了1、2、3、6測點本模型模擬的水位與其他模型模擬水位的比較，對比模型同樣包括有4.2中所述的各類模型。從圖8(a)中可以看出本模型模擬水位與其他模型水位變化趨勢基本一致，都呈現“雙峰”的形態，這是由于在t=1 min到t=4 min時發生了強降雨，而點源入流量峰值在t=37 min時出現。由于測點3位于下游聚洼地，其由降水形成的水位峰值持續了近30 min，而點源入流形成的峰值也出現了滯后且水位逐漸保持在一定高度。本模型模擬的流速與其他模型模擬流速的比較如圖8(b),由于ISIS Fast Dynamic模型無法預測速度，因此對比模型僅包括MIKE FLOOD以及InfoWorks ICM模型。圖8(b)顯示，各模型預測的速度變化趨勢大致相同，同樣都呈現“雙峰”狀態，但數值之間存在較大差異，這可能是由于各模型采用的計算方法和格式不同，或者地形影響、水深極淺區域處理方式的不同所引起的。總體上來說，本模型能夠處理干濕交替的復雜水流運動，能夠較為準確地模擬城市雨洪的演進過程。

圖7 研究區地形及道路、入流點、測點分布

4.4 渾、太河潰堤洪水在渾太胡同的演進

渾太胡同位于鞍山市溫香鎮，左鄰渾河右鄰太子河，兩條河流共同匯入形成大遼河，由大遼河河口流入遼東灣海域。每年夏季強降雨會對渾河、太子河汛期的水位產生極大的影響，據歷史資料統計，渾河及太子河曾頻繁發生潰堤洪澇災害，在從1798年到1997年近200年時間里，一共發生了114 次不同程度的水災，發生頻率為每1.75 a一次。其中1898年到1997 年間，發生了20次大面積洪水和內澇災害，遼、渾、太3河均包括在受災區域內。另外，連續多年發生水災的次數也很頻繁，近200年里，連續 2 年以上發生水災有22次，而且有些年發生的是較大洪水，對渾太胡同內的城區造成了極大的損害[16]。

本文以渾太胡同為研究區域，利用本文模型對渾、太河內100年一遇的特大降雨造成的潰堤洪水進行了數值模擬。模擬以Roe格式求解二維淺水控制方程，基于改進的動邊界處理以及負水深處理技術，提高了模型在復雜地形下模擬的穩定性和精度。馬家堡子在渾太胡同歷年發生的洪澇災害中多次形成潰口，如1949和1950 年發生的潰堤洪水災害。因此，本文首先選擇馬家堡子作為潰口位置，并且選擇了同在上游的遲坨子作為另一個潰口，首先模擬單個潰口情況下特大洪水在分布有城區和村莊的渾太胡同中地表水流的運動狀況；然后將兩個潰口進行組合，模擬兩個潰口同時潰堤形成的特大洪水在渾太胡同中地表水流的運動狀況。模擬區域非結構三角網格共12 659個，潰口位置處設置為入流邊界，潰口入流邊界網格長度較小，設置為80～100 m，為了減少計算量，其他區域的空間網格劃分尺度較大。特大降雨形成的徑流由入流邊界進入，其余邊界均為閉邊界。模擬中曼寧系數設置為0.05，采用變化的時間步長，模擬方案總時長為118 h。圖9為研究區內不同工況下潰堤洪水的淹沒范圍以及流場分布。

圖8 城市降水及點源徑流洪水演進算例不同模型計算的部分測點水位及流速比較

圖9 渾河和太子河潰堤洪水模擬不同工況下潰堤洪水在渾太胡同的水流演進

從圖9中可以看出，洪水由東北方向向西南方向演進，與渾太胡同地勢東北高西南低相符合。隨著潰口入流流量的進入，淹沒面積不斷擴大；在入流流量逐漸減弱后，上游水深明顯下降，洪水不斷向下游聚集，最終在下游處出現最大水深，馬家堡子、遲坨子以及組合潰口的最大淹沒水深分別為6.25、3.94以及7.44 m。同時也能看出，在入流徑流減弱后的洪水運動末期水流流速有了明顯下降，水流逐漸穩定。

圖10給出了幾種潰口工況下不同水深(H>0.5 m、H>1.5 m、H>2.5 m)淹沒面積的比較。由圖10可以看出，入流總量與各深度淹沒面積成正比，組合潰口情況下各深度淹沒面積均為最大。本文計算采用CPU單核串行模式，處理器型號為Intel Core CPU i7-6700 CPU@3.4GHz，模擬方案時段長度為118 h，計算耗時統計結果如下：馬家堡子、遲坨子以及組合潰口分別在19、12以及21 min內完成了模擬。結果表明，本文模型效率較高，能夠在較短的時間內完成大尺度洪泛區及城區的洪水預報。在模擬潰堤時，本文模型可同時模擬任意潰口、多個潰點洪水的實時演進預報?？傮w而言，本文模型能夠較為準確地模擬城區復雜地面上洪水的演進過程，并且可以實時地為城市洪水的演進進行預測。

圖10 渾河和太子河潰堤洪水不同工況下淹沒面積比較

5 結 論

基于有限體積法以Roe格式求解Riemann解并對方程進行離散，改進對動邊界以及負水深的處理技術，建立了非結構網格下的平面二維淺水模型。本文驗證了模型的干濕邊界處理能力以及靜水和諧性，并將模型應用于渾太胡同潰堤洪水模擬，最終得出以下結論：

(1)模擬拋物型地形下無阻力自由水面的運動，模擬的水位、水深以及流速值與解析值吻合，特別是干濕邊界處測點的模擬值基本沒有明顯的數據波動，因此模型具有較好的干濕邊界處理能力；模擬了洪泛區多個洼地的填充過程，洼地附近的水流基本平穩后，流速明顯降低，流場分布范圍減小，因此模型具有良好的靜水和諧性，模型能夠穩定地模擬復雜地形區域的洪水運動；

(2)模擬英國格拉斯哥城市雨洪的演進過程，結果很好地反映了在降水和點源徑流兩種城市洪水來源下的水位與流速變化，模擬結果具有較高的精度和準確性，模型可以針對復雜地形區域進行降雨和洪水演進的模擬；

(3)在不同潰口工況下，針對鞍山市渾太胡同100年一遇的特大洪水的模擬結果表明：洪水演進方向與地勢走向相符，洪水演進過程中流場分布合理，洪水流量與淹沒面積呈正相關。模型在21 min以內完成了不同工況下118 h的大面積洪水演進，具有較高的模擬效率。本文模型能夠在較短時間內對長時間、大范圍區域的洪水進行實時地預報，為預判洪水淹沒情況以及及時做出洪水治理決策提供科學的依據，對降低洪水災害損失和人員傷亡具有重要的意義。

在實際情況中，城市的降雨洪水在城市下墊面有一定的下滲現象，所以為使模擬更加精確，在未來研究中應將降水的下滲作用加入模型中。此外，在實際的城市雨洪中，暴雨洪水造成的污染物擴散是影響城市環境的重要因素，因此城市雨洪水動力學模型與水質模型的耦合也是本模型研究發展的一個重要方向。