長江干線宜昌—瀏河口段一維河網水流數學模型

李有為

(長江航道規劃設計研究院，湖北 武漢 430040)

20 世紀50 年代以來，國內外相繼開展了水流數值模擬的研究與應用，其中，一維水流數學模型經過數十年的研究已較為成熟，可普遍用于研究大型復雜河道長距離、長時段的水動力變化，其數值計算方法主要有特征線法、有限差分法[1]等。長江干線一維數學模型根據河道水流特點、計算精度及效率，通常以大通為界將宜昌以下河道分為宜昌—大通[2]、大通—瀏河口[3]兩段分別建立徑流和潮流模型，而在當前新水沙條件和新型復雜江湖河網關系下，傳統的兩分段模型難以整體性解決長河段復雜河網的水流模擬問題，也無法實現數字航道建設要求的一體化模擬。同時，針對圣維南方程組的隱式差分格式解法中壓縮存儲消元法、雙追趕法均存在著絕對值較小的數作除數而引起計算中斷或數值不穩定的隱患，特別是在超大型河網水力計算中將顯著加劇數值震蕩進而降低計算精度。多年以來，人們一直在探求如何壓縮系數矩陣的儲存信息、用較少的單位來儲存河網矩陣，用什么方法能夠較快地求解超大型河網方程組并提高計算穩定性，這一直是河網水力計算中的核心問題。

為了解決上述問題，針對長江宜昌—瀏河口段長達1 660 km 的超長河網，本文研發了穩定性強、精度高的基于有限體積法離散控制方程的長江航道一維河網水流數學模型。該模型采用具有TVD 特性的SLIC(slope limited centred scheme)格式計算數值通量，在處理河網交匯方面與基于三級聯解算法[4]的有限差分法有所不同，其以控制體積為計算單元，實質上不需要對汊道進行額外處理，對于大型復雜河網也不用進行復雜系數矩陣的求解程序，即可以實現自動對存在支流的河段進行水力要素的計算，該模型物理過程更為清晰明了、計算更加方便快捷。

本模型可為長江航道干支聯通規劃、工程方案設計、維護疏浚設計、水位預測、數字航道等做技術支撐，能促進超長距離、超大范圍航道綜合治理工程模型試驗研究，加速實現“暢通、高效、平安、綠色” 的黃金水道以及交通強國目標。

1 宜昌—瀏河口段各水文站區間來流特點

長江中下游流域面積達80 萬km2，覆蓋區域廣闊，沿線氣候條件復雜多變，依靠數量有限的控制水文站難以準確計算出站與站之間的未知區間來流量，而超長一維河網模型建立的重點、難點是能否平衡同期沿線水量、率定出合理糙率等因子。因此，在建模之前需要通過主要水文站點的實測數據分析宜昌—瀏河口段各水文站區間來流特點及水量平衡情況。

1.1 宜昌—大通河段

宜昌—大通河段主要支流為清江、松滋口、太平口、江漢運河、藕池口、漢江等6 條，主要入匯湖泊有洞庭湖、鄱陽湖，對應支流湖泊下游最近的水文站點分別有枝城站(清江)、沙市站(松滋口、太平口、江漢運河)、監利站(藕池口)、螺山站(洞庭湖)、漢口站(漢江)、大通站(鄱陽湖)。

從2018 年全年各主要水文站點的流量過程統計(圖1)、水文站點區間流量統計(圖2)來看:區間出入流對干線相鄰站流量影響比較明顯，對峰現時間則影響不大；在考慮主要支流湖泊出匯情況后，上下游最近水文站點的流量基本平衡，但仍然存在未經統計的次級水源出入(如毛細支流、局部降雨產匯流、人工取排水工程、蒸發下滲等)，具體而言，該部分水流占下游最近水文站流量比絕對值在10%以內的河道在全年90%以上時間的為宜昌—螺山段，螺山—漢口段、九江—大通段未知水流占比較多，特別是螺山—漢口段水流占下游最近水文站流量比絕對值均大于2%，九江—大通段水流占下游最近水文站流量比絕對值在2%以內的河道占全年時間僅25.2%，見表1。可見，螺山—漢口段、九江—大通段的區間未知水量需要重點考慮。

表1 日均區間流量占下游最近水文站流量比對應的時間

圖1 2018 年枝城—沙市站日均流量過程統計

圖2 2018 年各水文站點區間流量差值統計

1.2 大通—瀏河口河段

大通站以下較大的入江支流，北岸主要有裕溪河、滁河、淮河等，南岸主要有青弋江、水陽江、秦淮河以及太湖水系等。近年來長江下游兩岸各地在長江干流沿江修建了大量的引排水涵閘，阻斷了各支流的天然通江狀態。本節以大通站和徐六涇站之間的年內、年際來水量關系分析區間來流特點。大通站與徐六涇站年際年內徑流量比較見圖3。

圖3 2005—2018 年大通站與徐六涇站區間年際、年內水量變化量及變化率

1.2.1 年際變化

徐六涇站較大通站徑流量偏大，歷年平均偏大量為313.3 億m3，平均增幅為3.64%。歷年偏大變化率在0.22%～7.70%，增幅最大的為2017 年，增幅最小的為2013 年。兩站之間的水量變化主要受沿江通江河流的引排水影響，兩站各年的水量變化不一，其中，2005、2015 和2017 年區間水量變化比較大，分別為556 億、676 億、722 億m3，水量變化在500 億m3以上；2007、2008、2009 和2013 年變化比較小，在100 億m3以下。

1.2.2 年內變化

歷年年內各月只有5 月和6 月這兩個月徐六涇站的水量小于大通站，其它各月均是徐六涇站水量大于大通站；其中，兩站區間水量變化最大的是主汛期8 月，水量變化量達65.6 億m3；各月水量變化率最大的為10 月，變化率為8.16%。

綜上，宜昌—大通河段應重點考慮螺山—漢口段、九江—大通段的區間未知水源，大通—瀏河口河段年內、年際間區間來水量占總量比均小于10%，不需考慮區間未知水源。

2 數學模型建立

2.1 基本控制方程

根據水流的質量和動量守恒定律，建立以面積流量為變量的圣-維南方程組:

式中:t為時間；x為水流方向；g為重力加速度；A為過水斷面面積；Q為流量；zb為河底高程；H為水深；β為動量修正系數；Sf為阻力坡度，為綜合糙率，根據Cao 等[5]求解，R為水力半徑；ql為旁側入流的單寬流量；ulx表示單位流程上的側向出流流速在主流方向的分量。

該方程的水動力學部分較全面地考慮了河道不規則形態、河床與岸坡組成對水流動量通量的影響，并引入動量修正系數β處理。該模型適用于非恒定流，當時間偏導項為零時，退化為恒定流，因此本模型可普遍應用于恒定、非恒定流的計算。

2.2 數值解法

2.2.1 數值離散

為求解方程(1)(2)，采用Cao 等提出的數值解法，即首先對方程進行算子分裂，再采用有限體積SLIC 數值計算格式求解雙曲型方程組。方程(1)(2)可以寫成矩陣形式:

式中:U為變量矩陣；F為通量矩陣；Sb為除阻力項的源項矩陣；Sd為阻力的源項矩陣。

采用算子分裂分別對方程(3)進行離散，其中除阻力項的源項矩陣Sb采用二階龍格庫塔法計算源項，阻力的源項矩陣Sd采用全隱格式求解。

式中:Δt為時間步長；Δx為空間步長；i為空間節點；n為時間節點；Fi＋1∕2和Fi-1∕2為有限體積的SLIC數值格式求解的交界面通量；U**、U*為中間變量矩陣。對于公式(8)中的阻力源項矩陣Sd采用全隱格式[6]，該方法能夠處理在小水深條件下由于阻力過大引起的數值震蕩，提高計算的穩定性。

2.2.2 通量計算

在計算Fi-1∕2、Fi＋1∕2相鄰單元交界面上的通量時，采用SLIC 的數值格式求解，該數值格式方法能夠捕捉激波，自動處理部分河段出現急流的情況，同時亦可適用于恒定流的求解計算。該方法分為3 個步驟:

1)第一步數據重構。為避免寄生振蕩，在數據重構時采用limited slopes作為約束條件。在相鄰單元Ii＝[xi-1∕2，xi＋1∕2]、Ii＋1＝[xi＋1∕2，xi＋3∕2]中，已知位于時間節點n和空間節點i-1、i、i＋1 和i＋2 的變量可以得到

當系數α＝1 時模擬最小限制通量，當系數α＝2時模擬最大限制通量，本項目中取α＝1。則對空間上進行數據重構，上標R，L 分別代表i節點左右兩側:

2)第二步數據重構。在每個單元Ii上，取時間t＝1∕2·Δt，對時間上進行數據重構得

因此可求出相鄰單元交界面的通量Fi＋1∕2，代入式(6)進行求解，即可求出下一步各變量的值。

為了使計算穩定，對于顯格式時間步長需要滿足克朗條件:

即克朗數Cr＜1，λmax為由雅格比矩陣求得的最大特征值。

2.3 河網計算

與基于三級聯解算法的有限差分法不同，本文以控制體積為計算單元，實質上不需要對汊道進行額外處理，對于復雜河網也不用進行復雜系數矩陣的求解，程序即可實現自動對存在支流的河段進行水力要素的計算。

如圖4，對于由abcd組成的第i點的計算控制體單元Ii＝[xi-1∕2，xi＋1∕2]，從物理上而言，控制體中i節點上n＋1 時刻過流面積等于n時刻該節點上的過流面積加上該時間段內控制體水量的變化量。此時若有支流流入或流出，則支流流量也屬于控制體中流入或流出的水量，因此可得到其過流面積及水位值，這可由控制方程的連續性方程中體現出來。同理，控制體中i節點上流量的變化，也可以根據其動量守恒過程，考慮支流水流流入在水流方向上的投影對控制體動量的貢獻，從而可得到n＋1 時刻i節點的流量。

圖4 河網計算單元概化

在模型計算時，只需要輸入河段中支流的位置及其流入或流出的流量大小和流速方向，即可對有支流的情況進行求解，該求解過程與沒有支流時一致，因此并不需要額外增加節點導致增加計算量。

3 模型應用

3.1 河道概化

本文研究的模擬河道概化范圍見圖5，上起宜昌水文站、下至瀏河口(航道里程25.3 km)，包含漢江、北支兩條主要支流。由于河段較長，以經度帶為界，按石首、九江、安慶節點將長江干線劃分為4 個主干河段，其余分汊、支流河道以相互連接方式并入主干河段，進而形成完整的一維河網。模型對河網的26 條水道進行間距800～3 000 m的斷面剖分，共剖分出1 128 條斷面。模型外部邊界上游為宜昌站、仙桃站流量，下游為瀏河口附近的楊林站、北支的連興港站潮位，支流清江入匯采用高壩洲同期實測逐日流量過程，松滋口、太平口、藕池口三口分流分別采用新江口、沙道觀、彌陀寺、藕池(康、管)同期實測逐日流量過程，洞庭湖、鄱陽湖入匯采用城陵磯、湖口水文站實測逐日流量過程，江漢運河分流采用引江濟漢工程渠道設計流量固定值，螺山—漢口段分流、九江—大通段匯流采用前述區間水量平衡計算成果在區間水域沿線以線源方式設置相應補充水量過程，合計設置了15 處邊界條件。模型驗證起算地形以2018 年實測地形為主。

圖5 長江宜昌至瀏河口段一維河網全局布置

3.2 參數率定

本模型主要根據水位對模型中的綜合糙率系數進行率定。根據文獻資料搜集沿程河床泥沙組成特性，同時考慮到對于同一斷面的非恒定流過程其糙率與水位的關系并非單一，會受到漲水和退水、人工建筑物等影響。通過與枝城、監利、螺山、漢口、黃石港、九江、彭澤和大通等站實測水位數據的對比，率定成果見表2。

表2 宜昌—瀏河口段沿程綜合糙率驗證取值

3.3 驗證成果

長江干線為非恒定流，傳統工程上應用的一維模型通常將非恒定水流過程梯級化，每個梯級用恒定流模型計算，這種方法與直接采用非恒定流模型的計算結果差異仍然不是十分清楚。本模型直接采用非恒定模型計算，無需將非恒定來流過程概化為分級恒定流，從物理上減少了模型計算誤差，提高了模型精度。值得注意的是，在本模型中，當非恒定流模型的時間偏導項為零時，模型將退化為恒定流模型，同樣也適用于恒定水流過程的計算，因此模型具有良好的通用性；由于模型的數值計算方法使其具備了能夠自動捕捉激波和處理小水深的能力，對于長江干線局部的急流或急緩流交界的情況，能夠自行進行計算并保證模型計算的穩定性。

從模型過程驗證(圖6)可以看出:長江干線枯、洪期分界明顯，年內最高水位一般在7 月中旬，水位與流量成正比關系；大通及以上徑流河段實測水位與計算水位之間誤差均在10 cm 以內，誤差分布較均勻，大通站計算值在枯水期偏低、洪水期偏高；大通以下潮流河段徐六涇站在枯水期2018-02-22T00∶00—2018-02-28T23∶00 時的高低潮位計算誤差均在10 cm 以內、相位差為0 h，誤差分布較均勻；流量實測值與計算值走勢基本一致，兩者相對誤差在10%以內。因此，本次研究建立的新型一維河網數學模型計算穩定性強、精度較高，能較好地模擬超長網狀河段的水動力情況。

圖6 2018 年主要站點水位、潮位、流量過程驗證

4 結論

1)自主研發了基于有限體積法離散控制方程的一維河網水流數學模型，采用具有TVD 特性的SLIC 格式計算數值通量，在處理河網交匯方面以控制體為計算單元，可對存在支流的河段進行水力要素的自動快速計算，避免了傳統河網計算中對復雜系數矩陣的求解，提高了計算精度。模型可適用于恒定和非恒定流、自動處理急緩流交界處和小水深條件下的水流計算問題，能保證模型計算穩定性，具有較為廣泛的應用前景。

2)2018 年宜昌—瀏河口段區間來流計算結果表明應重點考慮螺山—漢口段、九江—大通段的區間未知水源，其余區段可暫不考慮。將建立的一維河網模型應用于該河段，計算結果表明:水位流量計算值與實測值基本吻合，能較準確反映超長河網復雜非恒定河道的水動力情況，模型具有較高的模擬精度。