摘    要：本文利用開源計算流體力學軟件OpenFOAM搭建二維數值模型研究管道-河道耦合水動力的作用機制，通過設置不同的管道上下游邊界水位和管道坡度探究不同降雨時期河道水位和管道坡度對雨水管道末端排水能力的影響。研究發現下游河道水位上升不一定降低管道排水能力，在一定強度降雨條件和管道坡度下，河道水位淹沒雨水管道末端一定程度時能促進管道提前形成滿管流，進而提高管道排水能力，而且坡度小的管道更容易在更低的河水位下形成滿管流。研究結果有助于更全面地揭示雨水管道末端與河道的水動力相互作用機制，為三維模擬與實驗研究的開展提供依據，也可為城市雨水管道改造及河道水位調控提供指導。

關鍵詞：城市內澇;雨水管道;耦合;數值模擬;OpenFOAM;

城市洪水是最常見和最具破壞性的自然災害之一[1]。極端暴雨的頻發、城市化進程的地面硬化、雨水管道的排水能力不足等都是導致城市內澇的重要因素。

在過去的幾十年里，國內外許多研究者通過理論分析、模型試驗和數值模擬的方法，對城市內澇的致災機理及解決方案進行了大量研究。在極端氣候方面，鄭飛飛等[2]分析了澳大利亞69個雨量計在1966—2012年間的數據，結果表明極端降水在不同時間尺度上發生變化會對洪水風險產生巨大的影響；呂洪等[3]在對鄭州進行城市洪澇災害預測時，利用極端降雨事件的影響構建新的風險指標體系，得到了更全面準確的預測結果；M. Bermúdez等[4]在西班牙西北部沿海河段進行的測試表明，氣候變化對河流流量的影響在具有復合洪澇潛力的地區具有非常重要的影響。在地表徑流方面，陳潔云等[5]總結了城市雨水徑流估算分析模型的發展，證明封閉式分析模型能高效準確地估算城市雨水徑流量；Martin Bruwier等[6]系統分析了9個城市特征在洪泛情況下對地表徑流的影響，表明洪水嚴重程度主要受建筑面積的影響；Francisco Peňa等[7]采用 FLO-2D 和MODFLOW-2005提出了一種基于物理、松耦合的建模框架，對降雨過程中地表水和地下水的相互作用導致的城市洪澇進行細化研究。

對于地表徑流和管道流耦合的過程，Noh Seong Jin等[8]在城市淹沒分析中研究了不同形狀雨水口的流量系數；Ricardo Martins等[9]采用計算流體力學(computational fluid dynamics, CFD)和試驗相結合的方法對不同排水條件下溝渠的水力特性進行了研究；陳倩等[10]對不同工況雨水口的泄流能力進行了試驗研究；還有一些研究人員在城市洪水耦合模型中對進水口進行了不同的模擬分析，并通過試驗進行了驗證[11,12]。在管道流動方面，Heekyung Park等[13]早在1998年就使用水動力學模型對環狀和樹狀管道系統進行了詳細的水動力評估；Pachaly Robson Leo等[14]評估了在雨水管理模型5.1(Storm Water Management Model 5.1,SWMM)中加入Preissmann槽算法后，在復雜、高度動態流入場景中的性能；葉家強等[15]研究在管道建模過程中通過入滲效率來估算管道排水量，其中入滲效率和降雨高潮位都有關系。

在雨水管末端區域，河道水位的頂托效應對其排水性能具有重大影響。Félix L等[16]綜述了前人多復合事件驅動因素耦合方法，并提出未來的研究應集中于開發一種緊密耦合的程序，以準確地反映風暴潮和降雨徑流之間復雜的物理相互作用；張兆祥等[17]采用水動力模擬，研究淹沒出流時水位、坡度、管徑對內澇的影響，結果表明，河道水位上升嚴重制約雨水管道排澇能力；羅鳴等[18]也研究了河道邊界水位對管道排水能力的影響，探究不同降雨重現期下出水口管道分別為自由出流、半淹沒出流和完全淹沒出流3種情形的管道相對排水能力變化。

針對城市內澇的成因，研究人員對降雨、地表徑流、地表和管道耦合、管道流動、管道和受納水體耦合的過程進行了細致的分析。盡管已有不少文獻對下游高水位和極端降雨的復合事件造成的城市內澇進行研究，但是大部分文獻僅是從宏觀上描述高潮位會阻礙管道排水從而導致城市內澇，沒有量化高潮位對管道排水能力的影響程度。且大多研究采用SWMM進行模擬，而SWMM在描述管道內水動力過程方面存在嚴重的不足，無法精準刻畫下游頂托情況下管道內的水流復雜情況。

本文主要通過二維模型對末端雨水管道進行模擬，以初步探究不同的降雨強度、水位高度、管道坡度對管道排水能力的影響。采用二維模型的主要原因是其可以高效地捕捉到管道中的明滿流交替現象，進而闡釋雨水管道末端與河道之間的水動力相互作用趨勢和規律，為復雜三維模擬及實驗研究的開展提供依據。而且，本文的二維模型模擬結果通過經驗公式驗證精度可靠。

1 數學模型

1.1 控制方程

OpenFOAM中的interFoam求解器可以求解2個不可壓縮相的雷諾平均Navier-Stokes方程，包括連續性方程(公式1)和動量守恒方程(公式2),并通過VOF(volume of fluid)方法捕捉自由表面的運動。在笛卡爾坐標系中，Navier-Stokes控制方程由以下幾個部分組成：

式中U表示速度矢量，ρ是流體密度，p*是偽動壓，g是重力加速度，X是位移矢量，μeff是有效動力黏度，μ是水或空氣的動力黏度，σ是表面張力系數，κ是界面的曲率。此外，為封閉求解雷諾時均N-S方程，本文采用k-omega SST湍流模型。

對于VOF方法，必須求解一個額外的方程來描述相的運動。相函數α被定義為每個單元中水體積的比例，需滿足一個經典的平流方程，如下式所示：

本文采用的VOF方法與壓縮技術相結合，有助于抑制數值擴散，保留自由表面的銳度。為了實現上述期望，Weller[19]在公式(3)中增加了一個人工壓縮項，最終表達式如下：

其中Ur是只在界面上工作的水相和空氣相之間的壓縮速度。此外，為了實現公式(4)的約束性，顯式求解的多維通用限制器算法(MULES)用來確保α的計算結果總是約束在0和1之間。上式中，ρ和μ都是通過使用相位函數對水和空氣進行加權計算：

ρ=αρwater+(1-α)ρair (5)

μ=αμwater+(1-α)μair (6)

1.2 數值方法

在求解控制方程(公式1～3)的過程中，采用了有限體積法(finite volume method, FVM)進行數值離散化。具體而言，計算域被離散成一系列的小網格，并使用Rhie和Chow[20]開發的拼合網格方法，將所有的流場信息都存儲在各個網格的中心。

采用OpenFOAM開展研究的優勢在于控制方程中具有多種不同項的離散化方法和插值方法。基于前人研究基礎[21]:本文的瞬態項采用隱式歐拉方法，平流項采用Gauss limitedLinearV1方法，黏性擴散項采用線性校正方法；除上述之外，所有的剩余項都采用了線性內插法；在求解壓力和速度耦合問題時，采用PIMPLE算法進行處理。PIMPLE是一種混合了壓力關聯方程(semi-implicit method for pressure linked equations, SIMPLE)算法的壓力隱式算子分裂法(pressure implicit with splitting of operators, PISO)和半隱式方法，主回路繼承自PISO。雖然本文中沒有使用欠松弛，但其允許方程欠松弛以保證所有方程在每個時間步收斂[22]。

2 模型驗證

2.1 數值模型

本文采用的二維數值模型如圖1所示。由于數值模擬要求上下游水位邊界水位盡量穩定，所以將水塔和水槽的長度設置得盡量長以滿足模擬要求，如圖1所示，水塔長度為350 m, 高度為1 m, 水槽長度為100 m, 高度為1 m。除此之外，排水管的直徑D為0.1 m, 水平方向的長度L為8 m, 坡度i取1%、0.5%和0%,管道下緣距離水槽底部0.25 m, G1,G2和G3是流量監測面。這些參數的選擇都是根據工程經驗以及實驗平臺構建方案決定。

圖1 二維數值模型示意圖 

2.2 網格收斂性

對于基于有限體積法的數值模擬，其計算精度不僅取決于離散化方案和計算算法，還與網格質量和網格大小有關，因此需要對網格的靈敏度進行驗證。

為了檢查網格的靈敏度，首先比較了管道中4種不同網格大小的數值模擬結果。為了捕捉自由表面，水塔中從排水管上沿至上方0.7 m處為加密區域，網格的縱向尺寸為2×10-3 m; 同理，水槽中從排水管上沿至上方0.6 m處為加密區域，網格縱向尺寸為1.85×10-3 m。由于管道進出口處水流狀態復雜，在水塔和水槽中靠近管道的水平方向設置10 m的網格加密段，該段網格水平方向的最小尺寸與管道水平方向最小尺寸相同。由于本數值模擬中主要關注的是管道部分，所以管道中需要有一個足夠詳細的網格尺寸來捕捉相關的流場特征。如圖2所示，在研究中，考慮了4種不同離散度的網格，分別是(a)粗網格、(b)中網格、(c)細網格、(d)最細網格。管道縱向分別劃分30、40、50、60份，最小網格尺寸分別為3.3×10-3、2.5×10-3、2×10-3 、1.7×10-3 m。網格的長徑比最大(Δx/Δz)為5,最小為1,在靠近管道兩端接口處最小。從兩端到管道中間各有水平方向1 m的網格長度漸變區。數值模擬采用時間步長為1×10-3 s, 計算至80 s, 取穩定后的流量平均值作為管道排水量。圖3比較了不同網格尺寸在測試工況下的管道排水量結果。粗網格和中網格的結果最大偏差接近1.01%,而對于較細和更細網格的結果偏差可以忽略不計。此外，細網格可以節省12%的計算量。因此考慮到計算精度和成本，在接下來的章節中，都采用了管道縱向網格尺寸最小為2×10-3 m的細網格。

圖2 4種不同離散度的網格用于數值收斂性研究 

圖3 不同尺寸網格的數值收斂研究結果 

2.3 經驗結果對比

本文模型在管道自由出流情況下可類比閘孔出流，閘孔出流流量計算經驗公式可參照Brater E F等[23]編制的水力學手冊，如公式(7)所示：

式中，h是自由出流下閘孔中心點水頭，a是閘孔長度，g是重力加速度，C=0.6是排放系數[24]。

為了描述方便，參考前人方法[18],針對排水管出口處的水位情況定義無量綱的淹沒度S,S指的是水槽的水位和排水管出口下緣的高度之差與排水管直徑的比值，具體如下：

其中，hc是水槽的水位；h0是排水管出口下緣的高度，h0=0.25 m; D是排水管的直徑。

選取管道坡度1%,Δh=0.30 m(Δh是水塔和水槽的水位差),管道兩端水位同時上漲的一組工況，對比數值結果和經驗結果之間的差異。其中由于只能在下游頂托較低時將模型看作自由出流，進而采用經驗公式計算流量值，故只計算S=0和S=0.25 2個工況下的經驗流量進行對比。如圖4所示，經驗結果與數值結果很接近，前者略大是因為經驗公式對應的是閘孔出流，閘孔后沒有管壁存在，而二維模型包含管道管壁對水流的阻尼效應。根據以上分析推斷本文中的二維模型可開展雨水管道末端與河道水動力相互作用規律的探析研究，模擬結果滿足工程精度要求。

圖4 數值和經驗結果的比較 

3 排水管道-河道耦合水動力過程

3.1 降雨初期

降雨初期，上下游水位同時上升。在數值模擬中假設水塔和水槽的水位差保持不變以簡化問題。降雨初期管道排水量隨頂托變化的模擬結果如圖5所示。

圖5(a～c)展示了當坡度為1%(工程中常見的雨水管道坡度),不同Δh條件下管道排水量隨S的增加先變大后趨于平穩。以Δh=0.20 m工況中管道流態的變化分析原因。如圖6所示，S=0.25時管道未形成滿管流，S=0.5時管道接近滿管流，S=0.75時管道已呈滿管流。管道中平均流速主要受Δh影響，Δh不變，S從0.25變化到0.75的過程中管道過流面積不斷增加，因此排水量增加；當 S>0.75 時，管道中均為滿管流，此時過水面積不變，在Δh相同的條件下，管道排水量基本保持不變。

圖5(d～f)展示了不同坡度管道在Δh=0.30 m 情況下排水量隨S的變化特征。i=0%和i=0.5%的管道在S=0.25后隨頂托的增加排水量保持不變，i=1%的管道在S=0.75后排水量保持不變，說明管道坡度越小越容易在較低的河道水位下形成滿管流。

圖5 降雨初期管道排水量與S的關系 

圖6 穩定后排水管入口處的水流狀態 

3.2 降雨中期

降雨中期，排水管上游來水基本穩定，下游河道水位持續上漲。在數值模擬中假定水塔水位固定，水槽水位上漲。降雨中期管道排水量隨頂托變化的模擬結果如圖7所示，增加管道排水量隨著下游頂托的呈現先增加后下降的趨勢，這與工程中普遍認為的管道排水量隨頂托增加而下降不同。

圖7 降雨中期管道穩定排水量與S的關系 

圖8 降雨末期管道穩定排水量與h水塔的關系 

圖7(a～c)展示了1%坡度的管道排水量隨S的增大先基本穩定，再突然增大，最后逐漸減小。在S=0到S=0.35之間排水量不變，這是因為隨著頂托的增加，水位差減少，導致流速降低；但同時下游水位上升，造成過流面積增加，在這兩者耦合影響下，最終流量變化不顯著。在S=0.5左右管道形成滿管流，此時排水量達到最大值，對于水塔水位為0.55、0.65 m和0.7 m的情況，管道排水量比自由出流時分別增加23%、35%和25%。在此之后，管道都是滿管流，過流面積不變，頂托增加導致水位差減小，流速減小，流量減小，水塔水位為0.55、0.65 m和0.7 m時，在S=1.5的情況下排水量比排水量最大時分別減少20%、28%和17%。

水塔水位不變，不同坡度管道排水量隨S的變化趨勢如圖7(d～f)所示。坡度為0的管道在S=0.05左右就已經達到排水量最大值，排水量比自由出流時增加8%,坡度為1%和0.5%的管道分別在S=0.25和S=0.4達到排水量最大值，排水量比自由出流分別增加17%和28%,這也證明了小坡度管道更容易形成滿管流。

3.3 降雨末期

降雨末期，排水管上游水位下降，下游河道水位基本穩定。在數值模擬中假定水槽水位固定，水塔水位下降。圖8展示了降雨末期管道排水量隨水塔水位變化的數值模擬結果。從結果可知，降雨末期管道頂托不變，排水量隨著水塔水位的降低而減小，這與普遍的認知一致。

4 結語

本文主要發現在降雨初期，頂托增加的同時上游來水量也在增大，導致管道排水量先增加，直到管道形成滿管流后排水量趨于穩定；在降雨中期隨著頂托的增加排水量會先達到一個最大值，達到排水量最大值的頂托程度和管道坡度有關；在降雨末期，頂托不變，管道排水量隨水塔水位的降低而減小。總的來說，管道坡度越小越容易在較小的頂托程度下形成滿管。研究結果可以為城市雨水管道改造及河道水位調控提供指導，對后續的三維模擬與試驗研究的開展具有重要指導意義。

文章來源：科技通報