王志良，徐立洲，丁志宏

（中水北方勘測設(shè)計(jì)研究有限責(zé)任公司，天津 300222）

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

基于有限體積法和二維淺水波方程，搭建了二維水動(dòng)力數(shù)學(xué)模型，對(duì)無壓輸水系統(tǒng)的充排水過程進(jìn)行了數(shù)值仿真分析，研究成果可為復(fù)雜城市輸水系統(tǒng)及長距離引水工程充排水過程的研究提供參考。

1 數(shù)值模型

1.1 數(shù)學(xué)模型

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

式中：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 數(shù)值方法

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

式中：i和j分別代表某一控制體在x和y方向的節(jié)點(diǎn)中心坐標(biāo)；n代表某一時(shí)刻。

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

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

2 模型驗(yàn)證

2.1 理想地形靜止算例

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

在無出入流條件下，計(jì)算域內(nèi)應(yīng)始終保持靜止。計(jì)算時(shí)段末的水位分布和流速分布分別如圖1和圖2所示，可以看出，計(jì)算時(shí)段末整個(gè)計(jì)算域水面仍保持在0.1 m，且保持靜止，即數(shù)值結(jié)果可很好地保持穩(wěn)態(tài)，與實(shí)際情況相符。

圖1 水位分布

圖2 流速分布

2.2 干床潰壩算例

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

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

圖3 水位分布

圖4 流速分布

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

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

3 數(shù)值實(shí)驗(yàn)

3.1 單向排水過程模擬

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

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

圖5 水位分布

圖6 流速分布

3.2 雙向充水過程模擬

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

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

圖7 水位分布

圖8 流速分布

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

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

計(jì)算時(shí)段末的水位分布和流速分布如圖9、圖10和圖11所示，可以看出，凸起前水位逐漸降落，凸起周圍出現(xiàn)繞流現(xiàn)象，流態(tài)較為復(fù)雜。水流流至凸起處受到阻力而涌起，凸起上游側(cè)及兩側(cè)水位驟升，水流的重力勢(shì)能迅速轉(zhuǎn)換為動(dòng)能，流速增加。由圖10和圖11可知，4 s時(shí)繞過凸起的兩股水流在下游側(cè)有合并的趨勢(shì)，6 s時(shí)凸起前側(cè)的流速增大，水流涌起區(qū)域擴(kuò)展，且兩股水流在下游側(cè)合并。

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

圖10 水位分布

圖11 流速分布

4 結(jié)論

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

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

實(shí)際輸水系統(tǒng)的布置較為復(fù)雜，系統(tǒng)的形狀、糙率及出入流狀態(tài)都會(huì)影響水流流態(tài)。因此，后期需進(jìn)一步在模型中考慮摩擦，并結(jié)合實(shí)際工程開展數(shù)值研究。