二維地表水動(dòng)力模型與SWMM耦合研究及應(yīng)用

王自明，查良瑜，王 斌

（浙江省水利河口研究院，浙江 杭州 310020）

1 問題的提出

城市雨洪的計(jì)算主要涉及產(chǎn)流計(jì)算、匯流計(jì)算以及地下管網(wǎng)水力計(jì)算[1-2]。目前對(duì)于產(chǎn)匯流的計(jì)算主要有2種途徑:一種是采用水文學(xué)方法；另一種是水動(dòng)力方法。水文學(xué)方法多將匯水區(qū)域當(dāng)做黑箱或者灰箱系統(tǒng)，建立輸入和輸出之間的關(guān)系以模擬坡面流[3]，其中較為代表的有SWMM[4]、Info Works CS[5]等。傳統(tǒng)水文學(xué)方法所建立的城市雨洪模擬多數(shù)應(yīng)用于城市規(guī)劃和設(shè)計(jì)階段，具有模型輕量、計(jì)算速度快的特點(diǎn)，但在實(shí)際城市內(nèi)澇模擬中難以表征地表澇水產(chǎn)生的過程和機(jī)理，模擬精度有限。水動(dòng)力學(xué)方法直接通過求解二維淺水方程來表征地表水流產(chǎn)匯流的物理過程，可為對(duì)城市內(nèi)澇產(chǎn)生過程、機(jī)理等方面的研究提供更直觀有力的支撐，具備更高的科研與使用價(jià)值。

但由于城市的復(fù)雜性，高精度的城市雨洪模型容易造成數(shù)值算法的不穩(wěn)定且計(jì)算耗時(shí)太長，實(shí)用性較差。解以揚(yáng)[6]等以基于無結(jié)構(gòu)網(wǎng)格上二維水力學(xué)模型為主體結(jié)構(gòu)，利用不同方法預(yù)報(bào)的降雨量作為系統(tǒng)的降雨邊界條件以及模擬城市排水系統(tǒng)產(chǎn)生的結(jié)果作為控制方程中連續(xù)方程的匯源項(xiàng)驅(qū)動(dòng)澇災(zāi)的形成和發(fā)展，但模型只用了2 281、678、646個(gè)網(wǎng)格對(duì)天津、南京、南昌3個(gè)大城市進(jìn)行離散，模型簡(jiǎn)化程度較大，難以表征樓房、道路等構(gòu)筑物的小尺度地形特征。與此同時(shí)，多數(shù)使用地表二維水動(dòng)力模型模擬城市內(nèi)澇的研究中，往往忽略地下管網(wǎng)對(duì)內(nèi)澇產(chǎn)生的影響或簡(jiǎn)單地采用水量扣除法概化管網(wǎng)對(duì)地表澇水的匯集與排除功能。近年來出現(xiàn)的一、二維耦合計(jì)算模型如MIKE -SWMM等，具備了一定耦合計(jì)算二維坡面產(chǎn)匯流和管網(wǎng)非恒定流的能力，但其精度和穩(wěn)定性難以匹配城市洪澇模型復(fù)雜的下墊面條件，且冗長的計(jì)算時(shí)間難以滿足城市內(nèi)澇預(yù)警、實(shí)時(shí)展現(xiàn)等更復(fù)雜應(yīng)用的需求。

本文構(gòu)建了高精細(xì)化、穩(wěn)定性能好的城市二維計(jì)算模型，通過引入GPU并行計(jì)算技術(shù)，解決了地表水動(dòng)力模塊的計(jì)算時(shí)間問題；同時(shí)耦合SWMM的地下管網(wǎng)計(jì)算結(jié)果作為地表連續(xù)方程的源匯項(xiàng)，建立更為準(zhǔn)確的城市內(nèi)澇模型，并將該模型應(yīng)用于東南某城鎮(zhèn)區(qū)域的內(nèi)澇計(jì)算和預(yù)報(bào)工作。

2 地表二維水動(dòng)力模型

文中使用的高性能高分辨率的二維水動(dòng)力模型，可以直接在城市的高分辨DEM上進(jìn)行模擬計(jì)算，為解決傳統(tǒng)模型在復(fù)雜地形條件下計(jì)算時(shí)間步長小，容易發(fā)散等問題，程序中將底坡源項(xiàng)作為單元界面上的通量，以便于達(dá)到全穩(wěn)條件，對(duì)摩阻源項(xiàng)采用半隱式算法，大大提高了模型的穩(wěn)定性。同時(shí)為解決在高分辨DEM上計(jì)算耗時(shí)較長的問題，模型中采用GPU加速技術(shù)，根據(jù)GPU和CPU構(gòu)造的不同，模型基于CUDA語言給GPU和CPU分配了不同的任務(wù)。GPU的計(jì)算任務(wù)主要為通過MUSCL線性插值以及HLLC近似黎曼解分別獲取單元體在預(yù)測(cè)步和校正步的單元界面通量，其并行計(jì)算能力可以同時(shí)完成多個(gè)單元體的通量求解過程；CPU負(fù)責(zé)模型的主程序，它的主要任務(wù)是利用GPU求解得到的界面通量對(duì)單元體的水深、流速等水力特征量進(jìn)行迭代計(jì)算。GPU和CPU的協(xié)作運(yùn)算，保證模型可以高效完成整個(gè)求解過程。模型的基本方程如下:

二維淺水方程，在地形復(fù)雜條件下無法保持內(nèi)力項(xiàng)與源項(xiàng)之間的平衡。為嚴(yán)格保持平衡并保證涉及復(fù)雜地形模擬計(jì)算時(shí)的精度，采用如下改進(jìn)形式的控制方程[7]:

式中:η為水面高（m）；η=h+zb，h為水深（m）；zb為河床高程（m）；t代表時(shí)間（s）；x，y分別代表平面2個(gè)方向的坐標(biāo)（m）；u，f，g，S分別代表包含著流體變量，x方向通量，y方向通量以及源項(xiàng)的矢量；q是考慮降雨等因素的點(diǎn)源項(xiàng)；u，v代表x，y方向的速度（m/s）；分別代表x、y方向上底床斜率；Cf代表著底床糙率系數(shù)。淺水下河床地形示意見圖1。

圖1 淺水下河床地形示意圖

模型使用Godunov型有限體積法求解上述控制方程。借助2步MUSCL - Hancock和龍格庫塔方法來保證其空間和時(shí)間的二階精度。自適應(yīng)時(shí)間步長，即借用CFL條件，由CFL條件確定的時(shí)間步長為:

式中:C為Courant數(shù)，取值在（0，1）之間，通常建議C 值設(shè)為 0.5 ～ 0.8。

3 地表及SWMM耦合

地表二維水動(dòng)力模型與SWMM模型耦合計(jì)算的示意見圖2，其中二維模型使用的是正方形網(wǎng)格，網(wǎng)格編號(hào)順序?yàn)閺淖笙陆情_始逐列從下到上所有的有效網(wǎng)格見圖3。網(wǎng)格左下角為坐標(biāo)原點(diǎn)，因此可以根據(jù)點(diǎn)坐標(biāo)以及正方形網(wǎng)格的大小，將點(diǎn)索引到二維模型所在的網(wǎng)格編號(hào)中。此外地表模型中通過2個(gè)文件來與地下管網(wǎng)模型來進(jìn)行數(shù)據(jù)交換，一個(gè)是網(wǎng)格編號(hào)及源匯文件名的mask.dat文件，文件中包含2列信息，分別為網(wǎng)格編號(hào)以及對(duì)應(yīng)源匯文件名稱，如網(wǎng)格編號(hào)10對(duì)應(yīng)的源匯文件名稱為1.dat；另一個(gè)是每個(gè)管網(wǎng)節(jié)點(diǎn)的源匯項(xiàng)文件，給出源匯強(qiáng)度的時(shí)間序列如網(wǎng)格編號(hào)10讀取的即為文件名為1.dat所對(duì)應(yīng)的源匯強(qiáng)度的時(shí)間序列。由此只需提取出SWMM模型中各個(gè)管網(wǎng)節(jié)點(diǎn)的溢流強(qiáng)度隨時(shí)間變化即可將地表二維與SWMM模型耦合起來。

SWMM生成的結(jié)果文件是一個(gè)二進(jìn)制文件，其中保存每個(gè)對(duì)象在輸出步長中各水力要素的結(jié)果，因此在讀取結(jié)果時(shí)需要調(diào)用動(dòng)態(tài)鏈接庫中的GetSwmmResult（iType，iIndex，vIndex，period，Value）函數(shù)，其中iType為結(jié)果中的對(duì)象類型（0 - 子流域，1 - 節(jié)點(diǎn)，2 - 管線，3 - 系統(tǒng)）；iIndex為被查詢對(duì)象的序號(hào)；vIndex為對(duì)應(yīng)iType的不同水力要素的結(jié)果，分為子流域、節(jié)點(diǎn)、管網(wǎng)，其中與地表二維模型交換過程中較為關(guān)注的為節(jié)點(diǎn)的變量，包括:0 - 節(jié)點(diǎn)底部高程以上的水深，1 - 水頭，2 - 蓄滯量，3 - 橫向入流，4 - 總?cè)肓鳎? - 淹沒溢出流量，需要與地表維進(jìn)行交換的數(shù)據(jù)主要為5 - 淹沒溢出流量；period為查詢時(shí)刻；Value為查詢結(jié)果的返回值，可得到相應(yīng)的溢流強(qiáng)度隨時(shí)間變化的文件，完成模型耦合。

一般情況下，地下管網(wǎng)的出流作為地表水的內(nèi)邊界，需要通過水力學(xué)計(jì)算公式如（堰流公式）得到。為了在不修改地表水模型的前提下實(shí)現(xiàn)二者的耦合，對(duì)此過程進(jìn)行一定的概化處理。認(rèn)為管網(wǎng)計(jì)算得到的每個(gè)節(jié)點(diǎn)流量，只在二維模型所在的一個(gè)單元格上釋放，再將流量通過換算成該單元格的降雨量，因此最終以降雨的形式，將二者耦合。如管網(wǎng)節(jié)點(diǎn)的出水量為0.01 m3/s，單元格的面積為1 m2，則認(rèn)為該單元格上的降雨強(qiáng)度為0.01 m/s。

圖2 地表模型和管網(wǎng)模型耦合示意圖

圖3 網(wǎng)格編號(hào)示意圖

4 耦合模型的應(yīng)用

某城鎮(zhèn)地區(qū)，三面環(huán)山一面朝海，地勢(shì)總體東北高，西南低，呈簸箕狀，城鎮(zhèn)下墊面類型復(fù)雜，同時(shí)擁有居民區(qū)、水庫、農(nóng)田、河道以及外海等不同的區(qū)域，區(qū)域總面積約25 km2。模型研究范圍見圖4。研究區(qū)域具有明顯的海洋性氣候特征，8月中旬至9月為多雨期，期間雨量的多少受臺(tái)風(fēng)影響較大。利用耦合模型對(duì)該區(qū)域在2009年“莫拉克”臺(tái)風(fēng)期間的城市內(nèi)澇情況進(jìn)行模擬。其中，地表水動(dòng)力模型所采用的網(wǎng)格直接由ArcGIS的柵格文件導(dǎo)出，其網(wǎng)格大小為3 m×3 m。管網(wǎng)主要分布于主城區(qū)的2條主干道上，其余區(qū)域管網(wǎng)較少，共分布有管網(wǎng)節(jié)點(diǎn)660個(gè)，其網(wǎng)格DEM見圖5。降雨強(qiáng)度見圖6。

圖4 研究區(qū)域圖

圖5 地表模型 DEM 及管網(wǎng)示意圖

圖6 “莫拉克”期間降雨強(qiáng)度圖

二維模型中，完全采用水動(dòng)力方法模擬地表產(chǎn)匯流過程，圖7為計(jì)算流域內(nèi)水流向水庫匯集的過程圖。實(shí)際調(diào)查中該城鎮(zhèn)老城區(qū)暴雨時(shí)內(nèi)澇較為嚴(yán)重，采用耦合模型對(duì)圖5中老城區(qū)進(jìn)行模擬計(jì)算，結(jié)果見圖8，模擬產(chǎn)生的老城區(qū)受淹范圍與調(diào)查情況基本一致，區(qū)域內(nèi)積水深度50 ～ 100 cm，與實(shí)際調(diào)查結(jié)果也較為接近。

圖7 產(chǎn)匯流過程圖

圖8 老城中的淹沒范圍及水深圖

分析老城區(qū)造成積水的原因主要有以下幾個(gè)方面:①地勢(shì)低洼；②城市建設(shè)導(dǎo)致原先的河道被大量占用甚至完全填埋（見圖9），上游山區(qū)河道保留較好，到達(dá)城區(qū)后，河道變?yōu)榘登?，且河道寬度和深度急劇減小，造成城區(qū)道路及房屋受淹較為嚴(yán)重； ③城區(qū)道路中排水管網(wǎng)系統(tǒng)薄弱，造成暴雨后雨水無處可排。

圖9 河道侵占造成城區(qū)淹沒圖

5 結(jié) 論

本文采用高性能二維模型對(duì)城鎮(zhèn)區(qū)域房屋、道路、河流等進(jìn)行精細(xì)化的模型，同時(shí)建立SWMM管網(wǎng)模型，將二者耦合起來，并將計(jì)算與實(shí)際現(xiàn)場(chǎng)調(diào)查的淹沒范圍及水深情況進(jìn)行對(duì)比。結(jié)果表明，耦合模型基本可以對(duì)城鎮(zhèn)的內(nèi)澇情況進(jìn)行準(zhǔn)確的捕捉。同時(shí)采用數(shù)據(jù)交換的方式將二維模型與SWMM模型耦合起來的方法無時(shí)間尺度問題，數(shù)據(jù)交換方便，可以為相似的工程提供借鑒。