金 溪，王 芳，張翔凌

（1.武漢理工大學(xué)土木工程與建筑學(xué)院，湖北武漢430070；2.武漢市規(guī)劃研究院，湖北武漢430014）

排水管網(wǎng)底泥的淤積和沖刷會造成許多問題。隨著城市化的發(fā)展這些問題表現(xiàn)得越來越突出，尤其是在使用合流制管網(wǎng)的地區(qū)這些問題的影響尤其嚴(yán)重。近幾年來排水系統(tǒng)底泥輸運造成的環(huán)境問題成為了相關(guān)領(lǐng)域的研究熱點［1-2］。

目前，相當(dāng)多的研究聚焦于將一維或二維水力模型與底泥輸運模型進行耦合，實現(xiàn)對于底泥沉積及沖刷過程的模擬［3-4］。但是研究的內(nèi)容集中于個別管段在單場降雨中導(dǎo)致的污染物釋放問題，無法在管網(wǎng)的規(guī)模上進行考慮底泥輸運問題的水力和水質(zhì)模擬，并且研究內(nèi)容缺少對于底泥沉積過程的考慮，因此無法在長周期上對底泥淤積量的變化進行預(yù)測。Seco等［5］提出了一種基于SWMM（storm water management model）與底泥污染物釋放模型的模擬方法，可以在管網(wǎng)規(guī)模上進行考慮底泥污染物沖刷的水質(zhì)模擬，但是模擬過程中SWMM與底泥污染物釋放模型是完全分開計算的，因此并不是一種耦合模擬過程。

為了實現(xiàn)考慮底泥輸運問題的排水管網(wǎng)水力、水質(zhì)模擬，提出一種基于SWMM及底泥輸運模型（sediment transport model，STM）耦合模擬的計算方法，在考慮底泥淤積和沖刷的情況下進行水力、水質(zhì)模擬。

1 控制方程

在管渠底泥輸運模型中，其控制方程常用Saint-Venant-Exner（SWExner）方 程 組 進 行 描 述。SWExner方程組的形式如下所示。連續(xù)性方程為

動量方程為

其中：x為流行距離，m；t為時間，s；A為管渠過流斷面面積，m2；Q為流量，m3·s-1；H為水頭，m；Sf為水力坡度；g為重力加速度，m·s-2；ρs為底泥干密度，kg·m-3；e為底泥空隙率；C為水中懸移質(zhì)濃度，kg·m-3；Ad為管渠底部沉積物橫斷面積，m2。

在求解SWExner方程組時可以利用經(jīng)典的計算流體力學(xué)方法以及Exner方程求解方法對其進行耦合求解。許多研究者利用該思路對SWExner方程組的求解進行了研究，并驗證了其可行性［6-9］。

SWExner方程組的求解方法可以分為兩大類：完全耦合求解方法（1-step approach）以及弱耦合求解方法（weak coupling approach）。對于完全耦合求解方法，在每一個時間步長的求解過程中，需要考慮底泥輸運對于管道流動的影響，并且管道的流動按照非恒定流進行考慮，因此SWExner方程組中的方程需要同時求解。而在弱耦合求解方法中，每一個時間步長中，管段的流動假定為恒定流動，并且假定底泥輸運對于水在管道中的流動沒有影響。通過上述假設(shè)，在弱耦合求解方法中可以先進行Saint-Venant方程組的求解，再將Saint-Venant方程組的水力計算結(jié)果作為已知條件代入Exner方程進行求解。多名研究者對完全耦合求解及弱耦合求解方法進行了綜合比較［10-11］，總體上來說完全耦合求解方法適用的場合更廣泛，尤其在流量及懸移質(zhì)濃度劇烈變化導(dǎo)致底泥形變或運動速度較快的場景下完全耦合求解方法可以給出更好的模擬效果，但是該方法也具有計算量大、耦合機制復(fù)雜的缺點。相對應(yīng)的，弱耦合求解方法只能用于底泥形變或運動速度相對于水流速度較低的情境下，否則其前提條件將無法滿足，影響模擬效果。但是弱耦合求解方法具有求解簡單、對Exner方程的形式適應(yīng)性更強等特點。Cunge等［12］認為，對于大多數(shù)底泥輸運問題，底泥的運動速度一般遠低于水流的速度，因此弱耦合求解方法可以應(yīng)對絕大多數(shù)底泥輸運問題的模擬場景。多名研究學(xué)者將弱耦合求解方法用于排水管段底泥輸運問題的求解，取得了較好的效果。驗證了弱耦合求解方法在排水管道場景應(yīng)用的可行性［13-14］。通過上述比較以及對排水系統(tǒng)中底泥輸運特點的考慮，采用弱耦合求解方法對SWExner方程組進行求解。

在耦合求解的過程中，SWMM5［15］被用來求解SWExner方程組中的Saint-Venant方程組部分。在SWMM5的水力、水質(zhì)模型求解過程中，管段被看作最基本的計算單元，忽略了各個變量數(shù)值沿管長的變化。因此當(dāng)Exner方程與SWMM5進行耦合計算時，其方程中的（?（QC））/?x部分被忽略掉，從而簡化為方程（4）所示的形式：

求解Exner方程的目的是獲得Ad與C的值。由于該方程中具有2個未知量，因此需要補充一個方程使其封閉。這里補充如公式（5）所示的河床形變方程：

其中：α為恢復(fù)飽和系數(shù)，表示在沖刷（或者沉積）過程中實際沖刷（或沉積）的底泥質(zhì)量占潛在的可能被沖刷（或沉積）的底泥質(zhì)量的比例；ω為底泥顆粒的沉降速率，m·s-1，沉降速率與底泥顆粒粒徑ds（mm）具有很強的相關(guān)性，本研究中ω與ds之間的數(shù)值關(guān)系采用Chen［16］總結(jié)的方程進行描述；B為底泥橫斷面表層寬度，m；C*為水流的攜沙能力，kg·m-3，這里采用張麗春等［17］中總結(jié)的攜沙能力計算公式進行描述：

其中：U為管渠過流斷面平均流速，m·s-1；R為管渠過流斷面水力半徑，m；k、m為經(jīng)驗參數(shù)，其取值與U3/（gRω）的數(shù)值有關(guān)，這里采用張小峰等［18］總結(jié)的k、m與U3/（gRω）的數(shù)值關(guān)系進行k、m值的計算。

2 求解方法

在耦合求解過程中SWMM5負責(zé)計算水力、水質(zhì)變量數(shù)值，包括Q、H、U、C，然后利用上述變量值代入Exner方程獲得Ad與C的值。在上述過程中C被更新了2次，SWMM5計算獲得的C為不考慮底泥顆粒沉積和沖刷因素情況下獲得的底泥顆粒在水中的濃度。而Exner方程的求解則是在考慮底泥輸運的情況下對底泥顆粒的濃度以及管底底泥的體積進行更新。通過上述耦合求解方法可以實現(xiàn)在考慮底泥輸運情況下的管網(wǎng)水力、水質(zhì)以及管底底泥淤積情況的模擬。該求解過程可以被描述為下述步驟：

（1）SWMM5演進一個時間步長，獲得Q、H、U、C數(shù)值。

（2）Exner方程求解的準(zhǔn)備。獲得ω、R、B的數(shù)值，根據(jù)U3/（gRω）的數(shù)值獲得k、m數(shù)值，最后計算C*的數(shù)值。

（3）Exner方程求解。聯(lián)合求解方程（4）、（5）獲得Ad與C的值。

（4）Exner方程求解結(jié)果的校驗與修正。Ad與C的計算結(jié)果沒有考慮耦合求解的數(shù)值穩(wěn)定性，因此需要對結(jié)果進行校驗與修正，其方法詳見2.2節(jié)。

（5）管段水力橫斷面的更新。根據(jù)Ad的數(shù)值對管段的水力橫斷面進行更新。返回步驟（1）或者結(jié)束模擬。

在上述求解過程中水中懸移質(zhì)濃度C可以利用SWMM5模型中定義的某一種污染物來代表，例如定義TSS（total suspended solid）代表水中懸移質(zhì)濃度。

2.1 管段水力橫斷面的更新

在耦合求解過程的第（5）步需要對管渠的橫斷面依據(jù)Ad進行更新，但是在SWMM5中沒有斷面類型可以在模擬過程中進行調(diào)整，因此需要對SWMM5的代碼進行適當(dāng)修改，增加一種可以在模擬過程中更新幾何尺寸的斷面類型。以SWMM5中的可淤積圓管（filled circular）斷面類型為基礎(chǔ)創(chuàng)建一種可以在模擬過程中改變斷面幾何參數(shù)的斷面類型?？捎俜e圓管中的淤積深度（filled depth）屬性被用來描述底泥淤積深度，并且該屬性需要在每一個時間步長中根據(jù)Ad的計算結(jié)果進行更新，更新方法如下所示：

（1）根據(jù)Ad計算結(jié)果更新淤積深度。

（2）根據(jù)淤積深度更新淤積表面的寬度以及滿流條件下的水力半徑。

（3）根據(jù)淤積深度更新管段進出口的偏移高度。

2.2 耦合求解過程中的數(shù)值穩(wěn)定性

耦合求解過程中需要將時間離散為較短的時間步長，而時間步長長度的確定對于求解過程的穩(wěn)定性具有非常重要的影響。通過上述耦合過程的描述可以看出，SWMM5模型為耦合求解過程的整體框架，Exner方程的求解被嵌入在SWMM5的步長中，因此耦合的時間步長長度是由SWMM5的時間步長確定機制控制的。SWMM5的時間步長是利用Courant-Friedrichs-Lewy（CFL）原則確定的，CFL原則確定的時間步長僅能滿足Saint-Venant方程組求解對數(shù)值穩(wěn)定性的要求，而無法保證Exner方程求解的數(shù)值穩(wěn)定性。

Exner方程求解出現(xiàn)數(shù)值不穩(wěn)定的情況主要是因為在求解過程中沒有考慮水流攜沙能力對于沉積物濃度變化范圍的限制。實際上在一個時間步長內(nèi)如果沉積物的濃度（C）達到了水流的攜沙能力確定的濃度（C*）后，沉積物的沉積與沖刷會達到一個平衡狀態(tài)，從而使得沉積物濃度停留在當(dāng)前時間步長內(nèi)水流攜沙能力確定的沉積物濃度數(shù)值。因此在Exner方程求解過程中，在獲得了C與Ad的計算結(jié)果以后需要利用C*進行校驗與修正。校驗與修正的過程為：

（1）首先通過比較沉積物濃度初始值（C0）與攜沙能力（C*）判斷當(dāng)前時間步長內(nèi)是沖刷過程還是沉積過程。

（2）如果為沉積過程，比較沉積物濃度計算結(jié)果C與C*的數(shù)值，如果C＜C*則賦值C=C*，并利用公式（4）重新計算Ad；如果為沖刷過程，比較C與C*的數(shù)值，如果C＞C*則賦值C=C*，并利用公式（4）重新計算Ad。

上述修正過程的流程可由圖1表示。

通過上述修正方法可以在不改變SWMM5時間步長確定機制的條件下，實現(xiàn)耦合過程的數(shù)值穩(wěn)定性。但是這種方法在處理高濃度、高沉降性（如粘滯性性小、粒徑較大）的沉積物顆粒時，仍然可能出現(xiàn)相鄰時間步長C、Ad計算結(jié)果差異較大的情況，雖然不會出現(xiàn)數(shù)值震蕩的情況，但是可能增加模擬過程的連續(xù)性誤差。這種情況下可以通過縮小SWMM5的時間步長（routing step）來減小連續(xù)性誤差。

圖1 沉積物濃度計算結(jié)果驗證與修正流程Fig.1 Flowchart of validation and correction sedi?ment mass concentration

3 案例研究

將上述耦合計算方法應(yīng)用2個不同的案例模型中。模型1為一個簡單的虛擬排水管網(wǎng)，其管網(wǎng)結(jié)構(gòu)如圖2所示，參數(shù)見表1。

圖2 案例1管網(wǎng)結(jié)構(gòu)Fig.2 Sewer network of case1

管網(wǎng)中各圖元的屬性如表1所示。該案例的作用是對于耦合算法的計算過程的數(shù)值穩(wěn)定性及連續(xù)性誤差問題進行分析。

表1 案例1管網(wǎng)模型基本參數(shù)值Tab.1 Parametersof features in sewer network of case1

案例2為巢湖市某區(qū)域排水管網(wǎng)系統(tǒng)，利用4場降雨的監(jiān)測數(shù)據(jù)對耦合方法的計算結(jié)果與SWMM5模型的計算結(jié)果進行對比，以評估耦合模型的模擬效果。案例2研究區(qū)域的管網(wǎng)圖如圖3所示。監(jiān)測降雨分別為當(dāng)?shù)貢r間2011-06-24，2011-07-12，2011-07-17以及2011-07-26的降雨。案例2管網(wǎng)模型的基本參數(shù)值如表2所示。

表2 案例2管網(wǎng)模型基本參數(shù)值Tab.2 Parameters of features in combined sewer network of case 2

研究區(qū)域內(nèi)管網(wǎng)在4場降雨過程中的水力模擬的準(zhǔn)確性是水質(zhì)模擬準(zhǔn)確性的基礎(chǔ)，因此必須首先對模型進行水力模擬校核，通過調(diào)整匯水子面積相關(guān)屬性的參數(shù)使模型的模擬結(jié)果與監(jiān)測流量數(shù)據(jù)達到較好的模擬效果。水力校核的結(jié)果如圖4所示。

4 結(jié)果與討論

4.1 案例1模擬結(jié)果與分析

為了模擬管網(wǎng)底泥在長周期運行中淤積深度的變化，對案例1模型利用弱耦合求解方法進行了4年的長周期模擬，其中降雨數(shù)據(jù)的時間步長為1h，管段底泥初始淤積深度為0.1m。管段底泥淤積深度以及水中TSS濃度的SWMM5模擬結(jié)果以及弱耦合求解方法模擬結(jié)果的對比如圖5所示。

圖4 研究區(qū)域管網(wǎng)水力模擬校核結(jié)果Fig.4 Calibrated result of SWMMhydraulic simulation of case 2

通過圖5 a、5b可以看出，通過耦合模擬方法可以模擬出管底淤積以及管網(wǎng)水質(zhì)情況隨降雨而產(chǎn)生的變化，在降雨密集的時間段，由于管內(nèi)流量的增加造成底泥沖刷，從而降低了管底的淤積深度，同時管道水流的TSS濃度顯著增高。而在降雨較少的時期，晴天流量中的TSS會在管底造成淤積，出現(xiàn)管底淤積深度增加而管道內(nèi)水流TSS濃度降低的情況。相對應(yīng)的，由于SWMM5模擬中沒有考慮管底淤積物的沖刷和淤積對于淤積深度和污染物濃度的影響，其管道淤積情況不會隨著時間發(fā)生變化，并且水中TSS的濃度沒有考慮雨天底泥沖刷的影響，因此濃度值普遍偏低。

圖5 案例1SWMM5模擬結(jié)果與耦合模擬結(jié)果對比Fig.5 Comparison of simulation results with and without coupling process of case 1

在耦合求解過程中數(shù)值穩(wěn)定性以及連續(xù)性誤差的大小是衡量求解方法質(zhì)量的重要指標(biāo)。在求解過程中，輸入的TSS質(zhì)量（Minlet）、排出的TSS質(zhì)量（Moutlet）以及由于沖刷或淤積導(dǎo)致的水中TSS質(zhì)量變化（Mdf）之間需要滿足質(zhì)量守恒原則，即滿足關(guān)系式：Minlet=Moutlet+Mdf，其中Mdf在沖刷情況下取負值，在淤積情況下取正值。在理想情況下這個關(guān)系式應(yīng)該被嚴(yán)格遵守，但是由于采用時間離散的方式進行求解不可避免會產(chǎn)生誤差，如果誤差較大則說明求解方法存在問題或者某些參數(shù)值設(shè)置不合理。圖6顯示了案例1的4年模擬過程中某個降雨沖刷時段（6h）的連續(xù)性誤差時間過程線。

圖6 案例1沉積物連續(xù)性誤差時間過程線Fig.6 Time evolution of sediment mass error of case1

從圖6可以看出，在整個沖刷過程中連續(xù)性誤差的值在很低的范圍內(nèi)波動。誤差的量化指標(biāo)可以通過相對誤差（E）來表示，相對誤差的計算公式如公式（7）所示：

在所截取的6h時間段內(nèi)，TSS連續(xù)性誤差的值為-3.18%，在整個4年的模擬時段中TSS的連續(xù)性誤差值1.83%。上述結(jié)果說明弱耦合求解方法在進行排水管網(wǎng)底泥輸運模擬中可以保持數(shù)值的穩(wěn)定性及較低的連續(xù)性誤差。

為了驗證耦合求解過程在污染物變化劇烈的條件下是否仍然能保持計算的穩(wěn)定性及低連續(xù)性誤差，將案例1中晴天流量TSS濃度提高至5 000mg·L-1，同時為了比較不同恢復(fù)飽和系數(shù)（α）及不同時間步長（Δt）對于數(shù)值穩(wěn)定性及連續(xù)性誤差的影響，對于沉降顆粒的恢復(fù)飽和系數(shù)及時間步長取不同的數(shù)值進行多次模擬。恢復(fù)飽和系數(shù)取值為0.010、0.025、0.050、0.10、0.250、0.500、1.000。時間步長取值為0.5s、1s、2s、5s、10s、30s、60s。α與Δt共有49種組合，對每種參數(shù)組合逐一進行模擬，圖7展示了所有參數(shù)組合對應(yīng)的模擬結(jié)果中TSS連續(xù)性誤差以及流量連續(xù)性誤差。

圖7 案例1不同α、Δt參數(shù)組合情況下的TSS及流量連續(xù)性誤差Fig.7 Flow and mass error of case1 at different combinations ofαandΔt

通過圖7可以看出所有α、Δt參數(shù)組合情況下耦合求解方法的流量連續(xù)性誤差均小于0.5%，時間Δt的不同取值幾乎對流量連續(xù)性誤差沒有影響。α取值對于流量連續(xù)性誤差有比較明顯的影響，這主要是因為連續(xù)性誤差的大小主要受到耦合過程中管段過流斷面變化的影響，而較高的α值會造成管段過流斷面變化速率較快（較高的沉積或沖刷速率會導(dǎo)致沉積物的快增加或降低），從而造成相對較大的流量連續(xù)性誤差。

對于TSS連續(xù)性誤差，α、Δt取值對其均有比較明顯的影響，并且在α、Δt均取較高值的情況下出現(xiàn)了很高的TSS連續(xù)性誤差。說明在底泥顆粒濃度高且易于沉降或沖刷的情況下較長的時間步長會導(dǎo)致連續(xù)性誤差的增加。但是可以看出隨著時間步長的縮短，TSS連續(xù)性誤差快速降低至較低的數(shù)值，說明縮小時間步長可以有效應(yīng)對高濃度且易于沉降或沖刷的底泥顆粒的輸運模擬?？梢钥闯鲈谂潘到y(tǒng)中，雖然在高濃度、高恢復(fù)飽和系數(shù)的底泥輸運情況下底泥的形變及運動速率有了較大的提升，但是其數(shù)值相對于水流速率仍然很低，因此當(dāng)時間步長長度較短的情況下將每個時間步長內(nèi)的流動近似為恒定流的前提條件仍然是成立的，因此弱耦合求解方法用于管道底泥輸運問題的求解問題的模擬仍然可以滿足數(shù)值穩(wěn)定性及低連續(xù)性誤差的要求。

4.2 案例2模擬結(jié)果與分析

在案例2的SWMM5模擬結(jié)果與弱耦合求解方法模擬結(jié)果對比的過程中，為了使2種不同的模擬方法均達到最佳模擬效果，2種不同模擬方法中的水質(zhì)模擬相關(guān)參數(shù)采用不同的取值。不同模擬方法的參數(shù)取值如表3所示。

表3 SWMM5及STM模型在耦合及非耦合模擬中的參數(shù)取值Tab.3 Parameters in SWMM5 and STM with and without coupling process

2種模擬方法獲得的TSS濃度模擬結(jié)果對比如圖8所示。通過圖8可以看出考慮底泥輸運的耦合求解方法的模擬結(jié)果與實測數(shù)據(jù)吻合得更好。SWMM5模擬結(jié)果的問題主要體現(xiàn)在模擬的后半段TSS濃度模擬值普遍較實測值偏低。這個現(xiàn)象在2011-07-12、2011-07-17以及2011-07-26三天的降雨體現(xiàn)得比較明顯，這主要是由于在沒有考慮底泥輸運情況下，當(dāng)?shù)孛嫖廴疚镌诒挥晁疀_刷完后，管道水流中的污染物量大大減少，造成TSS濃度非常低。而實際情況下水流對于管道底泥的沖刷會造成水中底泥顆粒仍然維持在一定的濃度。耦合模擬結(jié)果由于考慮水流對于管底沉積物的沖刷，因此模擬結(jié)果更加接近實測值。

2011-06-24的降雨事件中監(jiān)測的TSS濃度曲線中有2個峰值，分別出現(xiàn)在0：50以及2：00這2個時刻。SWMM5模擬結(jié)果中只模擬出2：00的TSS峰值，而耦合求解方法的模擬結(jié)果中2個濃度峰值都有較好的體現(xiàn)。這主要是因為在SWMM5模擬過程中，污染物只來自于地面污染物的沖刷，當(dāng)降雨強度較低的時候，各個匯水子面積上的地表徑流量很小，對地面污染物的沖刷能力很弱，因此地表徑流沖刷引入的污染物量幾乎對于管道中的TSS濃度沒有影響。耦合模擬過程中管道中的流量會在流動的過程中從上游至下游逐漸累加，因此管道內(nèi)的水流在從起端流至監(jiān)測點的過程中沖刷能力逐漸變強，并不斷對管底沉積物進行沖刷，從而形成比較明顯的濃度提升，形成第一個TSS濃度峰值。

圖8 案例2耦合求解方法與SWMM5模擬結(jié)果對比Fig.8 Comparison of simulation results with and without coupling process

為了定量檢驗2種方法模擬結(jié)果的質(zhì)量，對2種方法的模擬結(jié)果的納什效率系數(shù)（Nash-Sutcliffe efficiency coefficient，NSE）進行了計算。表4顯示了2種方法水質(zhì)模擬結(jié)果的納什效率系數(shù)及流量、污染物連續(xù)性誤差。

表4 2種求解方法的納什系數(shù)及連續(xù)性誤差對比Tab.4 Numerical accuracy evaluation of simulation results for each storm event

通過表4可以看出耦合求解方法對應(yīng)的NSE數(shù)值明顯較SWMM5的NSE值更高，說明耦合求解方法具有更好的模擬效果，耦合求解方法的連續(xù)性誤差雖然較SWMM5的連續(xù)性誤差高，但是仍然在較低且可接受的范圍內(nèi)。

5 結(jié)論

以SWMM5軟件為基礎(chǔ)框架，將Exner方程的求解方法嵌入SWMM5的水力、水質(zhì)計算過程中，實現(xiàn)SWExner方程組的弱耦合求解。利用該耦合求解方法可以實現(xiàn)在考慮管道底泥輸運情況下的管網(wǎng)水力、水質(zhì)模擬，從而實現(xiàn)在考慮管道底泥輸運的情況下，對排水管網(wǎng)的水力、水質(zhì)工況進行更加準(zhǔn)確的模擬和評估。

在該耦合求解過程的每一個時間步長中，Exner方程求解的結(jié)果需要利用水流攜沙能力進行校驗與修正，從而使得SWMM5原有的時間步長確定機制可以用來確定耦合過程的時間步長。通過對SWMM5中填充圓（filled circular）斷面類型的改造，實現(xiàn)了在模擬過程中更新斷面參數(shù)。目前的耦合求解方法中只能對于圓形斷面實現(xiàn)斷面參數(shù)的更新。本文2個案例的管段斷面均為圓形，因此可以利用耦合求解方法進行模擬，下一步工作會將SWMM5中的多個斷面類型進行改造實現(xiàn)對于沉積物斷面在模擬過程中的更新。

通過2個模型案例的模擬，驗證了該管道底泥輸運模型弱耦合求解方法的數(shù)值穩(wěn)定性及模擬結(jié)果的準(zhǔn)確性。通過模擬結(jié)果可以看出該方法的水流、污染物連續(xù)性誤差均較小，而且其模擬結(jié)果與單純的SWMM模擬結(jié)果相比，與實測的水力、水質(zhì)數(shù)據(jù)具有更高的吻合度。