基于淺水波方程的充排水過程模擬研究

王志良，徐立洲，丁志宏

（中水北方勘測設計研究有限責任公司，天津 300222）

修建輸水系統是實現區域間水資源均衡配置、解決供需水矛盾的重要手段。近年來，我國建設了多項長距離引水工程，有效促進了我國社會經濟發展。與此同時，大量關于引水系統充排水過程的研究被逐步開展。長距離管道輸水過程是涉及明滿流交替和氣液兩相流[1]的變化過程，楊開林等[2]提出了模擬無壓隧洞充水過渡過程的“虛擬流動法”，求解了萬家寨引黃入晉工程洞內初始無水條件下的充水問題。楊敏等[3]基于圣維南方程、假想窄縫法和虛擬流量法，建立了長距離串聯輸水管線的連續充水數學模型，并應用于南水北調天津干線有壓系統充水模擬。王克忠等[4]采用Fluent軟件對長距離無壓引水洞岔洞的水力特性進行了模擬分析。

基于有限體積法和二維淺水波方程，搭建了二維水動力數學模型，對無壓輸水系統的充排水過程進行了數值仿真分析，研究成果可為復雜城市輸水系統及長距離引水工程充排水過程的研究提供參考。

1 數值模型

1.1 數學模型

模型控制方程是服從靜水壓力分布假設的二維淺水波方程，這是描述非恒定漸變流的基本微分方程，其向量形式如下：

式中：h為水位（m）；u為x方向流速（m/s）；v為y方向流速（m/s）；g為重力加速度（m/s2）；Sox和Soy分別為x方向和y方向的地形坡度；z為地形高程（m）；Sfx和Sfy分別為x方向和y方向的摩擦力，本文暫不考慮。

1.2 數值方法

有限體積法是將計算域劃分為若干控制體，通過計算進出控制體邊界的通量，基于質量和動量守恒定律，得到時段末各控制體上的物理量分布[5]。該方法物理意義明確，具有較好的守恒性，擬采用該方法結合Rusanov格式求解淺水波方程。計算域被劃分為形狀規則的矩形網格，首先將方程（1）離散并在控制體內積分，x和y方向的空間步長分別為Δx和Δy，時間步長為Δt，離散方程為：

式中：i和j分別代表某一控制體在x和y方向的節點中心坐標；n代表某一時刻。

有限體積法的核心是構造邊界通量，Rusanov格式的具體形式為：

式中：λ1、λ2、λ3為方程（1）的Jacobi矩陣的特征值。

2 模型驗證

2.1 理想地形靜止算例

該算例中計算域的長和寬均為25 m，空間步長Δx=Δy=0.25 m，時間步長Δt=0.125 s，模擬時長為100 s，初始條件和地形方程為：

在無出入流條件下，計算域內應始終保持靜止。計算時段末的水位分布和流速分布分別如圖1和圖2所示，可以看出，計算時段末整個計算域水面仍保持在0.1 m，且保持靜止，即數值結果可很好地保持穩態，與實際情況相符。

圖1 水位分布

圖2 流速分布

2.2 干床潰壩算例

為了準確模擬輸水系統的充排水過程，數值模型需具備處理干濕界面的能力，該模型中設計水位小于10-6m時，即為干河床。干床潰壩算例的計算域長度為10 m，寬度為1 m，空間步長Δx=Δy=0.05 m，時間步長Δt=0.02 s，模擬時長為6 s。壩體位于x=5 m處，初始條件和地形方程為：

y方向上的物理量變化率為0，計算時段末x方向的水位分布和流速分布分別如圖3和圖4所示，可以看出，潰壩發生后，落水波傳播至上游，上游水位下降，漲水波傳播至下游，下游水位上升，且潰壩波未傳播到的區域仍保持開始狀態。特別地，下游落水波未傳播到的區域可保持干床狀態，這表明該模型可有效處理干濕界面。為了說明該模型的數值計算精度，將數值解與精確解進行對比，可以看出，計算水位與精確解誤差較小，而計算流速除下游間斷處外誤差較小。

圖3 水位分布

圖4 流速分布

使用誤差計算公式（10）-（11）進行計算，得到h和u的計算誤差分別為0.004和0.1045，滿足計算要求。

式中：eu和eh分別代表流速和水位的計算誤差；unumi,j和分別代表（i，j）處的控制體的流速的數值解和精確解；和分別代表（i，j）處的控制體的水位的數值解和精確解。

3 數值實驗

3.1 單向排水過程模擬

該算例上游入口為Wall邊界，下游出流為20 m3/s，計算域長度為2000 m，寬度為30 m，系統坡度為0.001%，空間步長Δx=Δy=2 m，時間步長Δt=0.05 s，初始條件和地形為：

排水100、1800、3600 s時的水位分布和流速分布如圖5和圖6所示，可以看出，排水前期落水波逐漸向上游傳播，出口處流速最大，上游水位為6 m，且保持靜止；排水后期水面線為直線，越往下游流速越大。由圖5和圖6可知，出流100 s后，落水波傳播至距入口1183 m處，出口處流速為5.6 m/s；排水1800 s后，系統水位自入口0.57 m降為出口0.41 m，出口處流速為2.12 m/s；排水3600 s后，系統水位自入口0.19 m降為出口0.14 m，且出口流速為1.28 m/s。

圖5 水位分布

圖6 流速分布

3.2 雙向充水過程模擬

該算例為雙向同時充水問題，上、下游充水流量均為20 m3/s，計算域長度為2000 m，寬度30 m，坡度為0.001%，空間步長Δx=Δy=2 m，時間步長Δt=0.05 s，初始條件和地形為：

充水1800 s和3600 s時的水位分布和流速分布如圖7和圖8所示，可以看出，上下游同時充水過程中，漲水波由系統兩端向內部傳播，匯合處水位最高，流速最小。由圖7和圖8可知，充水1800 s后，匯合段長162 m，水位0.03 m，匯合段上、下游側流速分別為3.15、3.12 m/s；充水3600 s后，匯合段長356 m，水位為0.04 m，匯合段上、下游側流速仍分別為3.15、3.12 m/s。與匯合處相比，系統兩端的水位波動較大。

圖7 水位分布

圖8 流速分布

3.3 有障礙物的單向充水過程研究

輸水系統的底面形狀會影響水流的流態，以下算例模擬了有凸起的輸水系統中的充水過程。入口邊界施加水深0.7 m，出口為自由出流邊界。計算域的長度為25 m，寬度為5 m，計算域內部有一圓錐形凸起，其中心位于（12 m，2.5 m）處，高度為1 m，其余區域平坦。空間步長Δx=Δy=2 m，時間步長Δt=0.05 s，模擬時長為6 s，初始條件為：

計算時段末的水位分布和流速分布如圖9、圖10和圖11所示，可以看出，凸起前水位逐漸降落，凸起周圍出現繞流現象，流態較為復雜。水流流至凸起處受到阻力而涌起，凸起上游側及兩側水位驟升，水流的重力勢能迅速轉換為動能，流速增加。由圖10和圖11可知，4 s時繞過凸起的兩股水流在下游側有合并的趨勢，6 s時凸起前側的流速增大，水流涌起區域擴展，且兩股水流在下游側合并。

圖9 水位分布（t=4 s）

圖10 水位分布

圖11 流速分布

4 結論

（1）基于淺水波方程搭建了二維水動力數值模型，采用有限體積法和Rosanov格式進行求解。靜止算例驗證了該模型可有效保持靜水狀態和適應復雜地形，干床潰壩算例驗證了該模型可有效處理干濕界面，且計算精度較高。

（2）采用該模型模擬了具有坡度的輸水系統中的單向排水過程和雙向充水過程，不同時刻的水位、流速分布可反映出輸水系統中的落水波和漲水波傳播規律；采用該模型模擬了具有凸起的輸水系統中的充水過程，模擬結果較好地反映出了障礙物周圍的繞流現象，均符合實際物理過程。

實際輸水系統的布置較為復雜，系統的形狀、糙率及出入流狀態都會影響水流流態。因此，后期需進一步在模型中考慮摩擦，并結合實際工程開展數值研究。