基于改進垂向流量交換的城市內(nèi)澇模擬方法

金 溪，周鵬飛，張翔凌，劉承宇

(武漢理工大學(xué)土木工程與建筑學(xué)院，湖北 武漢 430070)

近年來隨著極端暴雨災(zāi)害的頻繁發(fā)生，越來越多的城市開展了排水防澇規(guī)劃與建設(shè)工作。在這個過程中城市內(nèi)澇水文水動力模型得到了越來越廣泛的應(yīng)用[1-4]。早期水文水動力模型側(cè)重于地面產(chǎn)流過程及地下管網(wǎng)水動力過程的模擬，無法模擬地面內(nèi)澇[5]。為彌補這一缺陷，Hsu等[6]采用SWMM模型與地表二維模型耦合實現(xiàn)了城市內(nèi)澇耦合模擬。但該方法中模型間的水流交換是單向的，當(dāng)管網(wǎng)超載時，水流從檢查井溢出到地表；當(dāng)?shù)叵鹿芫W(wǎng)有足夠過流能力時，地表水流不能回流至管網(wǎng)中。隨著研究不斷深入，逐漸出現(xiàn)一批能夠?qū)崿F(xiàn)地下地表水流雙向交換的耦合模型[7-9]。耦合模型雖然具有結(jié)果詳細(xì)、擅長處理方向不定水流模擬的優(yōu)勢，但計算過程復(fù)雜、計算效率低。為了達到模擬精度與計算效率之間的平衡，部分研究者開始探索既能避開求解繁瑣物理方程又具有良好精度的內(nèi)澇模擬方法，例如基于GIS 的淹沒分析算法、改進填洼算法的洪澇計算以及部分簡單元胞自動機模型[10-12]。但這些方法只能確定淹沒范圍和水深，不能模擬洪澇演進的水動力過程，因此應(yīng)用范圍有限。

目前耦合模型是城市內(nèi)澇模擬的主要方法。垂向流量交換是耦合水文水動力模型的關(guān)鍵環(huán)節(jié)，對模擬結(jié)果的準(zhǔn)確性及水量平衡會產(chǎn)生重要影響[13-15]。目前的垂向流量交換主要以一維模型節(jié)點水位(H1D)、二維模型地面水位(H2D)之間的差值為作用水頭，利用孔口或堰流公式計算一維模型與二維模型之間的交換流量。該方法計算過程復(fù)雜并且獨立于一維模型，沒有考慮垂向交換流量與管道流量間的相互影響，容易造成計算不穩(wěn)定，需要通過設(shè)置最大允許回流量等措施實現(xiàn)耦合計算的穩(wěn)定性。

本文將節(jié)點進出水量平衡原則引入垂向流量交換計算過程，對耦合模型進行改進，簡化垂向流量交換計算過程，提高耦合模擬的穩(wěn)定性。采用理論案例及實際案例對改進方法進行驗證，分析改進模型的可靠性和模擬精度，以期為耦合模型垂向流量交換提供一種新方法。

1 水文水動力耦合模型

對于耦合模型的耦合方式，國內(nèi)外已有學(xué)者提出了一些較為符合實際的耦合理論并構(gòu)建了耦合模型。耦合模型中一維、二維模型的時間同步以及水量交換是耦合方法的關(guān)鍵。對于一維管網(wǎng)模型與地表二維模型耦合，主要存在垂向水流交換，即水流在檢查井、雨水篦處的交換；對于一維河網(wǎng)模型與地表二維模型耦合，水流交換主要存在水平方向上的正向交換和側(cè)向交換。

本研究著重討論一維管網(wǎng)模型與地表二維模型間的耦合過程。耦合模型中一維模型采用SWMM作為水文水動力計算引擎進行地表水文及管道系統(tǒng)水動力計算；二維模型采用守恒型淺水方程作為控制方程，利用自主開發(fā)的二階Godunov格式有限體積法求解模塊進行求解。

2 改進垂向交換流量計算方法

節(jié)點垂向流量交換的3種典型情況如圖1所示。圖1(a)、圖1(b)表示節(jié)點超載情況(節(jié)點井內(nèi)承壓且節(jié)點水位(HJ)大于節(jié)點相連管段的管頂高程點最大值(HP)，由于一維模型一般忽略節(jié)點井筒內(nèi)空間，此時也表示一維節(jié)點水位超過地面高程(H0))，此種情況下一維模型管段流量的求解應(yīng)該考慮節(jié)點位置處二維模型地表水位，并且通過管段流量的計算結(jié)果獲得節(jié)點的入流量(Qbf)及出流量(Qof)。出流量為節(jié)點出流管段的流量之和，入流量包括節(jié)點入流管段的流量、地表徑流流量、其他進流量等，將入流量及出流量的差值作為垂向交換流量，垂向交換流量大于0則代表溢流，小于0代表回流，等于0代表沒有垂向流量交換。圖1(c)表示節(jié)點未超載情況(HJ

圖1 節(jié)點垂向流量交換的典型情況Fig.1 Three typical scenarios of node vertical exchange flow

2.1 垂向交換流量計算

將一維模型節(jié)點與其所在的二維模型地表單元格進行關(guān)聯(lián)，方法是將地表單元格面積作為節(jié)點的積水面積，展示節(jié)點積水面積與地表單元格之間的關(guān)系，如圖2所示。圖2中，深藍色單元格為節(jié)點所在位置對應(yīng)的地表單元格，耦合模擬開始前利用該單元格面積更新對應(yīng)節(jié)點的積水面積，該單元格除參與二維模型計算外還參與垂向流量交換的計算。其余單元格為普通單元格，只參與二維模型計算。

圖2 節(jié)點積水面積與地表單元格關(guān)系Fig.2 Relationship between node ponded area and cells of 2-D mesh grid

2.1.1 節(jié)點超載情況下垂向交換流量的計算

當(dāng)節(jié)點處于超載狀態(tài)時，相連管道流量的計算需要考慮地表水位的影響，在一維水力模擬每個時間步長動力波計算開始之前利用節(jié)點所在地表單元格水位更新節(jié)點水位，然后執(zhí)行動力波水力計算，從而實現(xiàn)考慮二維模型中地表水位變化情況下的一維模型水力模擬。利用節(jié)點進出水量守恒原則計算節(jié)點的垂向交換流量。根據(jù)水量守恒原則，節(jié)點積水面積內(nèi)水體積變化與節(jié)點相連管段流量之間存在如下關(guān)系：

(1)

式中：Vold為時間步長開始時刻節(jié)點積水面積內(nèi)積水體積，m3；Vnew為時間步長結(jié)束時刻節(jié)點積水面積內(nèi)積水體積，m3；Δt1D為一維模型時間步長，s；Qi為與節(jié)點相連管道的流量，m3/s，入流為正，出流為負(fù)；Qlat包含地表徑流量、旱季流量等其他進流量之和，m3/s；m為節(jié)點相連管道的數(shù)量。

一維模型中節(jié)點積水面積內(nèi)的水體積變化為一維、二維模型之間流量交換的體積，即

(2)

式中：Qexf為垂向交換流量，m3/s，正值表示節(jié)點溢流，負(fù)值表示節(jié)點回流。

2.1.2 節(jié)點未超載情況下垂向交換流量的計算

當(dāng)節(jié)點不超載但地面有積水時，相連管道流量的計算不需要考慮節(jié)點位置處地表水位的影響，因此不進行節(jié)點水位的更新，僅計算節(jié)點回流流量。本文采用底孔口流量公式進行回流流量的計算，回流過程中地表水位會發(fā)生變化，采用式(3)所示的變水頭流量公式作為回流流量計算的基礎(chǔ)公式：

(3)

式中：t為回流時間，s；S為積水面積，m2；H1為回流時間t起始時刻節(jié)點水位，m；H2為回流時間t結(jié)束時刻節(jié)點水位，m；μ為孔口流量系數(shù)，取值為0.62；A為節(jié)點井筒或雨水口面積，m2。

假設(shè)H2=0，利用式(3)求出節(jié)點所在單元格地面積水全部回流所需時間t。如果t>Δt1D，說明當(dāng)前SWMM時間步長內(nèi)單元格內(nèi)積水沒有全部回流，利用式(4)、式(5)計算當(dāng)前時間步長結(jié)束時刻地面積水水位及回流流量：

(4)

(5)

如果t≤Δt1D，說明當(dāng)前SWMM時間步長內(nèi)單元格內(nèi)積水全部回流，利用式(6)計算當(dāng)前時間步長結(jié)束時刻回流流量：

(6)

如果在時間步長結(jié)束時節(jié)點為超載狀態(tài)，說明節(jié)點在該時間步長內(nèi)由回流狀態(tài)轉(zhuǎn)化為溢流狀態(tài)，此時節(jié)點的垂向交換流量采用進出水量守恒原則利用式(2)進行計算。

2.1.3 修正節(jié)點水位更新造成的連續(xù)性誤差

圖3 節(jié)點水位更新引起水量不平衡原因示意Fig.3 Diagram of continuity error caused by node head updating

(7)

2.2 模型耦合實現(xiàn)方式

采用改進垂向流量交換計算方法耦合模型中的SWMM模型需要開啟節(jié)點積水功能(啟用Allow Ponding)。在Δt1D對應(yīng)的耦合過程中具體計算步驟如下：

(1) 根據(jù)展示的3種情況進行節(jié)點是否超載判斷；

(2) 超載節(jié)點進行水位更新操作，未超載節(jié)點進行回流流量的計算；

(3) 一維模型演進一步；

(4) 根據(jù)水量守恒原則計算超載節(jié)點垂向交換流量；

(5) 連續(xù)性誤差消除；

(6) 垂向交換流量傳遞給二維模型；

(7) 二維模型演進至Δt1D結(jié)束時刻。

由于SWMM模塊具有計算節(jié)點積水體積變化的功能，節(jié)點超載情況下的垂向交換流量可以通過SWMM節(jié)點積水體積的變化直接獲得。上述計算步驟中的(1)—(5)步可通過向SWMM模塊添加相應(yīng)的函數(shù)在SWMM模塊內(nèi)完成，降低一維、二維模型間耦合度，使得2個模塊更加獨立。

3 案例分析

3.1 理論案例

采用如圖4所示的理論案例對改進方法進行驗證。管網(wǎng)系統(tǒng)由6個節(jié)點和6根管道組成，地表區(qū)域為一個邊長為200 m的方形閉合平原區(qū)域。平原區(qū)域覆蓋的節(jié)點2—節(jié)點5均允許垂向流量交換，管道和節(jié)點的屬性信息如表1所示。平原高程為0 m，糙率為0.025。水流通過節(jié)點1進入系統(tǒng)，入流量在模擬開始的10 min內(nèi)由0逐漸增加至1.0 m3/s，然后保持不變至第20 h，然后降為0；水流通過節(jié)點6離開系統(tǒng)，出口設(shè)置為自由出流。系統(tǒng)初始為干，即管道和平原的初始水深均為0 m。

圖4 理論案例示意[16]Fig.4 Diagram of theoretical case

表1 管網(wǎng)要素屬性表

利用Infoworks ICM軟件建立理論案例模型作為對照。ICM模型采用堰流公式進行垂向流量交換計算，堰流系數(shù)取值0.36。在改進垂向流量交換方法中節(jié)點所在單元格的面積會作為節(jié)點的積水面積，因此，網(wǎng)格劃分對于耦合計算產(chǎn)生的影響需要進行討論。采用2種單元格面積(4、16 m2)進行理論案例網(wǎng)格劃分。

圖5(a)展示了單元格面積16 m2網(wǎng)格劃分下改進模型的水量平衡情況。在垂向流量交換過程中節(jié)點溢流量應(yīng)該等于地面積水量與節(jié)點回流量之和，通過圖5(a)可以看出改進模型在垂向流量交換過程中保持了水量守恒，節(jié)點累積溢流曲線與“節(jié)點累積回流+地表積水”曲線幾乎完全重合。在一維、二維耦合模擬過程中，進入耦合模型的總水量應(yīng)該等于流出系統(tǒng)的水量及在模型中積水水量之和。通過圖5(b)可以看出模型整體水量平衡也得到了很好的滿足。

圖5 改進模型水量平衡分析Fig.5 Analytical diagram of water balance of improved model

圖6以節(jié)點2為典型溢流節(jié)點，節(jié)點3為典型回流節(jié)點，A點及節(jié)點3位置作為地面淹沒深度代表點，展示了2種網(wǎng)格劃分(S4(網(wǎng)格面積4 m2)、S16(網(wǎng)格面積16 m2))情況下改進模型與ICM模型計算結(jié)果的對比。利用納什效率系數(shù)(Nash-Sutcliffe efficiency coefficient，ENS)對2個模型模擬結(jié)果的吻合程度進行量化評價。通過對比可以看出，采用不同網(wǎng)格劃分情況下改進垂向流量交換方法的計算結(jié)果與ICM模型模擬結(jié)果吻合度較高。在米級、十米級地表單元格劃分條件下，采用改進垂向流量交換方法獲得的耦合計算結(jié)果與ICM模型獲得結(jié)果基本一致。通過圖6(c)看出，網(wǎng)格劃分的不同對于地面積水深度的模擬結(jié)果有一定影響，但是在不同網(wǎng)格劃分的情況下改進方法均可以獲得與ICM模型幾乎吻合的地面水深變化曲線。

圖6 改進模型與ICM模型計算結(jié)果的對比Fig.6 Comparison of improved model and ICM model results

3.2 城市內(nèi)澇模擬案例

采用武漢市武昌區(qū)東沙匯水子系統(tǒng)作為城市內(nèi)澇模擬案例，研究區(qū)域管網(wǎng)匯水范圍為82.73 km2。匯水區(qū)域內(nèi)地勢比較平坦，南側(cè)為沙湖與東湖，湖泊為雨水的主要受納水體，由于湖泊水位常年保持在較高水位，對于雨水入湖造成一定的影響，在局部低洼地區(qū)容易發(fā)生內(nèi)澇。

利用Infoworks ICM及改進模型分別建立研究區(qū)域耦合模型。ICM模型地面網(wǎng)格采用非結(jié)構(gòu)網(wǎng)格，單元格平均大小為450 m2，單元格數(shù)量為21萬個。改進模型地面網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格，單元格尺寸為20 m×20 m，單元格數(shù)量為24萬個。降水?dāng)?shù)據(jù)采用當(dāng)?shù)?016年6月18日及7月6日的2場暴雨(以下簡稱20160618場次、20160706場次)，2場暴雨的累計降水量分別為143.5、125.6 mm，降水強度過程線如圖7所示。

圖7 暴雨降水強度過程線Fig.7 Research region storm system network and rainfall intensity curves of selected storms

采用改進垂向流量交換方法的耦合模型對于2場降水進行模擬后獲得的地面淹水情況如圖8所示。20160618場次降水量大、降水歷時較短，導(dǎo)致了比較嚴(yán)重的城市內(nèi)澇，模擬結(jié)果可以看出地面積水點較多、積水深度較深，積水范圍與當(dāng)日實際記錄的內(nèi)澇點基本吻合。20160706場次雖然降水量也比較大，但是降水歷時較長，并且在降水過程中雨量分布比較均勻，峰值降水強度約為20160618場次峰值雨強的一半，因此形成的內(nèi)澇較輕，當(dāng)日記錄的內(nèi)澇點較少。通過模擬結(jié)果可以看出，采用改進垂向流量交換方法的耦合模型可以比較準(zhǔn)確地模擬地面內(nèi)澇淹沒范圍及深度。

圖8 2場暴雨地面淹水情況模擬結(jié)果Fig.8 Flood areas of simulation results of two selected storms

圖9 不同模型地面積水體積時間過程線對比Fig.9 Comparison of flood volume-time curves of different models

圖9展示了改進方法模型與ICM模型在2場降水情況下地面積水體積隨時間變化的情況?？梢钥闯?，改進方法模型與ICM模型的垂向流量交換計算結(jié)果幾乎一致，地面積水體積隨時間變化的曲線重合度較高，說明改進方法節(jié)點垂向流量交換計算結(jié)果與ICM模型計算結(jié)果吻合度較高。

在城市內(nèi)澇模擬場景中超載節(jié)點的垂向流量交換數(shù)量占了總體垂向流量交換數(shù)量的絕大多數(shù)。通過對模擬結(jié)果的統(tǒng)計，在2場暴雨模擬過程中超載節(jié)點垂向流量交換數(shù)量占總垂向流量交換節(jié)點數(shù)量的比例分別為81.43%和88.66%，采用改進方法可以在絕大多數(shù)情況下利用SWMM模塊通過水量守恒原則直接計算出節(jié)點的垂向流量交換結(jié)果，僅在少量的非超載節(jié)點回流情況下需要對回流流量采用變水頭孔口流量公式進行計算。從而大大減少了采用孔口或堰流公式進行垂向流量交換的數(shù)量，提高了耦合模擬的計算效率及計算穩(wěn)定性。

4 結(jié) 論

以SWMM 模型與自主開發(fā)的二維水動力模型為基礎(chǔ)，構(gòu)建城市內(nèi)澇水文水動力耦合模型，在此基礎(chǔ)上采用節(jié)點水量平衡原則對垂向流量交換方法進行改進。采用理論管網(wǎng)及城市排水管網(wǎng)作為案例，對改進方法從模擬精度、可靠性及水量平衡等方面進行了驗證。主要結(jié)論如下：

(1) 通過水量守恒原則計算節(jié)點垂向流量交換可以極大減少孔口或堰流相關(guān)的計算，計算過程可以由SWMM計算模塊完成，簡化了耦合模型的垂向流量交換的計算復(fù)雜程度。

(2) 非超載情況下的節(jié)點回流量計算通過增加相關(guān)函數(shù)的方式融合在SWMM計算模塊中。SWMM模塊在計算過程中采用迭代方法計算節(jié)點水頭及管段流量，并以節(jié)點水頭收斂為迭代停止條件，實現(xiàn)了節(jié)點垂向流量交換與管段流量的計算穩(wěn)定性，避免出現(xiàn)節(jié)點溢流、回流交替出現(xiàn)的數(shù)值震蕩現(xiàn)象。

(3) 通過與主流耦合模型軟件infoworks ICM模擬結(jié)果進行對比，表明采用改進垂向流量交換方法的耦合模型具有良好的精度和可靠性，在米級、十米級、百米級地表網(wǎng)格劃分情況下與infoworks ICM模擬結(jié)果高度吻合，可應(yīng)用于城市內(nèi)澇數(shù)值模擬計算。