李有為

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

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

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

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

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

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

1.1 宜昌—大通河段

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

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

表1 日均區(qū)間流量占下游最近水文站流量比對應(yīng)的時間

圖1 2018 年枝城—沙市站日均流量過程統(tǒng)計

圖2 2018 年各水文站點區(qū)間流量差值統(tǒng)計

1.2 大通—瀏河口河段

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

圖3 2005—2018 年大通站與徐六涇站區(qū)間年際、年內(nèi)水量變化量及變化率

1.2.1 年際變化

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

1.2.2 年內(nèi)變化

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

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

2 數(shù)學(xué)模型建立

2.1 基本控制方程

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

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

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

2.2 數(shù)值解法

2.2.1 數(shù)值離散

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

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

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

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

2.2.2 通量計算

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

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

當(dāng)系數(shù)α＝1 時模擬最小限制通量，當(dāng)系數(shù)α＝2時模擬最大限制通量，本項目中取α＝1。則對空間上進行數(shù)據(jù)重構(gòu)，上標R，L 分別代表i節(jié)點左右兩側(cè):

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

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

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

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

2.3 河網(wǎng)計算

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

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

圖4 河網(wǎng)計算單元概化

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

3 模型應(yīng)用

3.1 河道概化

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

圖5 長江宜昌至瀏河口段一維河網(wǎng)全局布置

3.2 參數(shù)率定

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

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

3.3 驗證成果

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

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

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

4 結(jié)論

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

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