三峽庫區(qū)枝狀河網(wǎng)水動力過程實時模擬

，,b，，

(華中科技大學(xué) a.水電與數(shù)字化工程學(xué)院； b.數(shù)字流域科學(xué)與技術(shù)湖北省重點實驗室, 武漢 430074)

1 研究背景

為快速高效地模擬河網(wǎng)水流過程，反映河流真實形態(tài)，大量學(xué)者已經(jīng)研究開發(fā)了多種一維河網(wǎng)水流數(shù)學(xué)模型，建模過程中常見的離散格式主要為六點中心Abbott差分格式及Preissmann四點隱式差分格式等[1-5]。這些模型在計算中具有較好的穩(wěn)定性；計算結(jié)果精度較高、可靠性好、適用范圍廣、可擴展性良好。但是在計算速度與效率方面略有不足，在模擬較長時段水流運動過程時需花費較長時間，實時性較差。而在現(xiàn)有諸多實際工程問題中，尤其在汛期防洪問題及水庫實時調(diào)度問題中可能需要利用水力學(xué)方法進行實時模擬計算，在此背景條件下，亟需開發(fā)一種高效的河網(wǎng)水流數(shù)學(xué)模型滿足工程需求。為此，以提高計算效率為目標，引入粒子追蹤思想，構(gòu)建通用河網(wǎng)拓撲關(guān)系，研究開發(fā)一種顯隱式結(jié)合的河網(wǎng)水流數(shù)學(xué)模型。并在此基礎(chǔ)上，建立了長江干流朱沱—三峽壩址河段的河網(wǎng)模型，計算檢驗了模型的穩(wěn)定性、準確性與高效性，可進行水動力過程的實時模擬。

2 高效求解的河網(wǎng)水流數(shù)學(xué)模型

2.1 模型原理

2.1.1 控制方程

針對一維河網(wǎng)水力學(xué)計算問題，基于定床假設(shè)，假定波動水面呈漸變趨勢，斷面水面線無橫比降，過水斷面壓力在垂線上滿足靜水壓力分布規(guī)律，河(渠)道底坡足夠小，僅考慮水頭沿程損失，且非恒定流阻力可借用恒定流阻力公式進行計算。在上述假定條件的基礎(chǔ)上，可使用圣維南方程組作為水流數(shù)學(xué)模型的控制方程組，其基本形式如式(1)、式(2)所示，其中式(1)為連續(xù)方程，式(2)為動量守恒方程。

B(?η/?t)+?Q/?x=0 ；

(1)

?Q/?t+gA(?η/?x)+?(Q2/A)/?x+

gn2Q|Q|/(AR4/3)=0 。

(2)

式中：B為斷面水面寬度；η為斷面中心水位；Q為斷面流量；A為過水斷面面積；g為重力加速度；n為河道綜合糙率值；R為斷面水力半徑；t為時間；x為距離。

2.1.2 離散格式

研究采取以斷面為中心的方式構(gòu)建控制體，其示意圖如圖1(a)所示，圖中豎實線代表斷面，虛線所圍成的矩形構(gòu)成了水位控制體，豎虛線表示水位控制體的前后流量介面。

圖1 控制體及離散格式示意圖Fig.1 Diagram of control body and discretization format

在保證模型計算穩(wěn)定、準確的同時，亦考慮模型的高效求解，采用θ半隱的離散格式[6]對動量守恒方程的水位梯度項進行離散，如圖1(b)所示。得到簡化的動量守恒方程離散形式如式(3)所示，相應(yīng)得到連續(xù)方程的θ半隱離散形式如式(4)所示，對阻力項則采用顯式格式進行離散，其形式如式(5)所示。

(3)

(4)

(5)

式中：i為水位控制體編號，取值為0，1，2，…；N為當前時刻，取值為0，1，2，…；N+1為下一時刻；Qi±1/2為第i個水位控制體前后介面中心流量；Ai±1/2為第i個水位控制體前后介面處過水斷面面積；θ為隱式因子，一般取值在0.5～1.0之間；ηi為第i個水位控制體水位；Δxi+1/2表示第i個斷面到i+1個斷面的距離；Li為第i個水位控制體長度；Ri±1/2表示第i個水位控制體前后介面處水力半徑。

2.1.3 歐拉-拉格朗日粒子追蹤法

為準確、快速得到各介面流速大小，引入歐拉-拉格朗日粒子追蹤法[7](ELM)計算對流項。對于一維水流，通常以斷面中心的狀態(tài)來描述整個斷面的情況，故可將某一斷面概化為其中心處的單一粒子，推求該粒子在時間、空間內(nèi)運動的變化。假設(shè)從tN時刻到tN+1時刻，某粒子從位置x運動到位置x′，若不考慮質(zhì)量損失，則粒子運動的速率一定，若想獲取該粒子tN+1時刻的運動狀態(tài)，則只需獲得其在tN時刻所在的空間位置與速率；根據(jù)河道水流方向及當前位置的前后斷面中心點坐標計算可得到粒子的運動方向。任一粒子隨時間、空間變化的運動狀態(tài)可用式(6)來描述。

(6)

式中：uj為第j個粒子的速率；degj為第j個粒子的運動方向，即與水平線所呈的夾角大?。籜i，Yi分別為第i個斷面中心點的橫坐標和縱坐標。

2.1.4 枝狀河網(wǎng)預(yù)測校正算法

計算將河網(wǎng)劃分成多個河段，對每個河段運用預(yù)測校正法[8]，將單步求解分成2步，能夠在保證計算精度的同時提高計算速度。

(1)預(yù)測步：每一計算時段內(nèi)，對每一計算河段，根據(jù)邊界條件、上一時段末狀態(tài)與水力參數(shù)求解各控制體水位。

(2)校正步：將預(yù)測步計算得到的各河段末斷面的水位作為新的下邊界條件，根據(jù)預(yù)測步的上邊界條件與參數(shù)再次求解各斷面的水力要素，由此計算得到各斷面各時刻較為準確的水位值。

2.2 模型求解步驟

基于上述離散原理，結(jié)合實測數(shù)據(jù)資料，依照地形數(shù)據(jù)，構(gòu)建一維河網(wǎng)水流數(shù)學(xué)模型，其求解的主要步驟如下所述。

(1)構(gòu)建拓撲關(guān)系。根據(jù)斷面基本信息數(shù)據(jù)，建立河道、斷面、控制體及計算介面間的拓撲關(guān)系，依據(jù)斷面平面坐標，構(gòu)建斷面空間位置結(jié)構(gòu)，在此基礎(chǔ)上，規(guī)定計算正方向。

(2)初始化參數(shù)。讀取控制體糙率值、初始水位，初始化各斷面及各計算介面上的水位、流量等水力要素。

3 實例計算與驗證

以長江干流朱沱—三峽大壩壩址河段為研究對象，計算驗證模型的合理性、準確性，測試模擬時間，論證實時模擬的可能性。概化計算河段如圖2所示，朱沱—三峽壩址河段長約760 km，河段自上而下沿程設(shè)有朱沱、寸灘、清溪場、萬縣、廟河5個水文站，壩址附近鳳凰山水位站為壩前水位代表站，2大支流嘉陵江和烏江的入庫控制站分別為北碚站和武隆站，干流河道庫面寬度一般為700～1 700 m，丘陵峽谷交替變化，地貌奇特，地形復(fù)雜，具有一定的代表性[9]。模型上邊界條件為朱沱站流量過程，下邊界條件為鳳凰山站水位過程，主要中間邊界條件為北碚站和武隆站流量過程。

圖2 計算河段水系概化圖Fig.2 Generalizationofcalculatedriversystem

選取2005年6—10月實測數(shù)據(jù)進行參數(shù)率定，選取2006年6—9月實測數(shù)據(jù)進行驗證，率定得到計算河段的河道綜合糙率值為0.034～0.045。具有代表性的寸灘站、清溪場站的率定與驗證結(jié)果如圖3所示。各站點水位、流量過程的驗證結(jié)果良好，水位計算誤差一般 ≤ 10 cm，最大誤差約為20 cm；流量計算誤差一般 ≤ 5%，最大誤差約為 10%，數(shù)學(xué)模型計算結(jié)果與實測資料符合較好，實現(xiàn)了高精度的水流模擬，模型構(gòu)建合理。

圖3 水位和流量率定驗證結(jié)果Fig.3 Results of calibration and verification of water level and discharge

在此基礎(chǔ)上，以1 a為計算總時段，測試模型計算效率。計算環(huán)境搭建在Windows 7操作系統(tǒng)下，運用Intel Xeon E5-2697a v4(clock speed: 2.6 GHz)處理器進行單核計算，統(tǒng)計50次計算平均耗時約為23.7 s。由此可見，模型計算效率較高，多次運算時間差別不大，穩(wěn)定性較好，可以滿足三峽庫區(qū)水動力過程實時模擬的要求，從運算耗時的角度來看，模型能用于三峽水庫實時調(diào)度研究。

4 結(jié) 語

以實時模擬三峽庫區(qū)水動力過程為研究目標，引入歐拉-拉格朗日法處理動量守恒方程的對流項，采取θ半隱離散格式求解水位梯度項，運用有限體積法離散連續(xù)方程，在此基礎(chǔ)上建立一維河網(wǎng)水流數(shù)學(xué)模型，并基于預(yù)測校正法利用追趕法求解以河段為單元離散構(gòu)建的三對角矩陣系統(tǒng)方程組?；谒Ｐ?，模擬實測數(shù)據(jù)下三峽庫區(qū)朱沱—壩址河段水流演進過程，率定驗證了模型具有較好的穩(wěn)定性與精度，并通過大量年時段計算，檢驗了模型具備實時模擬的能力，實現(xiàn)了研究目標。

參考文獻：

[1] 謝鑒衡. 河流模擬[M]. 北京: 水利電力出版社, 1990:68-72.

[2] 汪德顴. 計算水力學(xué)理論與應(yīng)用[M]. 北京: 科學(xué)出版社, 2011:139-154.

[3] 衣秀勇, 關(guān)春曼, 果有娜,等. DHI MIKE FLOOD洪水模擬技術(shù)應(yīng)用與研究[M]. 北京: 中國水利水電出版社, 2014:4-5.

[4] 夏軍強, 張曉雷, 鄧珊珊,等. 黃河下游高含沙洪水過程一維水沙耦合數(shù)學(xué)模型[J]. 水科學(xué)進展, 2015, 26(5):686-697.

[5] 陳煉鋼, 施 勇, 錢 新,等. 閘控河網(wǎng)水文-水動力-水質(zhì)耦合數(shù)學(xué)模型——Ⅰ.理論[J]. 水科學(xué)進展, 2014, 25(4):534-541.

[6] HU De-chao, ZHONG De-yu, WANG Guang-qian,etal. A Semi-implicit Three-dimensional Numerical Model for Non-hydrostatic Pressure Free-surface Flows on an Unstructured, Sigma Grid[J]. International Journal of Sediment Research, 2013, 28(1):77-89.

[7] HU De-chao, ZHANG Hong-wu, ZHONG De-yu. Properties of the Eulerian-Lagrangian Method Using Linear Interpolators in a Three-dimensional Shallow Water Model Usingz-level Coordinates[J]. International Journal of Computational Fluid Dynamics, 2009, 23(3):271-284.

[8] HU De-chao, ZHONG De-yu, ZHANG Hong-wu,etal. Prediction-Correction Method for Parallelizing Implicit 2D Hydrodynamic Models. I: Scheme[J]. Journal of Hydraulic Engineering, 2015, 141(8): 04015014.

[9] 黃仁勇. 長江上游梯級水庫泥沙輸移與泥沙調(diào)度研究[M]. 北京: 科學(xué)出版社，2017:14-16.