長(zhǎng)江干線(xiàn)宜昌—瀏河口段一維河網(wǎng)水流數(shù)學(xué)模型

李有為

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

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

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

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

1 宜昌—瀏河口段各水文站區(qū)間來(lái)流特點(diǎn)

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

1.1 宜昌—大通河段

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

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

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

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

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

1.2 大通—瀏河口河段

大通站以下較大的入江支流，北岸主要有裕溪河、滁河、淮河等，南岸主要有青弋江、水陽(yáng)江、秦淮河以及太湖水系等。近年來(lái)長(zhǎng)江下游兩岸各地在長(zhǎng)江干流沿江修建了大量的引排水涵閘，阻斷了各支流的天然通江狀態(tài)。本節(jié)以大通站和徐六涇站之間的年內(nèi)、年際來(lái)水量關(guān)系分析區(qū)間來(lái)流特點(diǎn)。大通站與徐六涇站年際年內(nèi)徑流量比較見(jiàn)圖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 月這兩個(gè)月徐六涇站的水量小于大通站，其它各月均是徐六涇站水量大于大通站；其中，兩站區(qū)間水量變化最大的是主汛期8 月，水量變化量達(dá)65.6 億m3；各月水量變化率最大的為10 月，變化率為8.16%。

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

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

2.1 基本控制方程

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

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

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

2.2 數(shù)值解法

2.2.1 數(shù)值離散

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

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

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

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

2.2.2 通量計(jì)算

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

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

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

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

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

為了使計(jì)算穩(wěn)定，對(duì)于顯格式時(shí)間步長(zhǎng)需要滿(mǎn)足克朗條件:

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

2.3 河網(wǎng)計(jì)算

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

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

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

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

3 模型應(yīng)用

3.1 河道概化

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

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

3.2 參數(shù)率定

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

表2 宜昌—瀏河口段沿程綜合糙率驗(yàn)證取值

3.3 驗(yàn)證成果

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

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

圖6 2018 年主要站點(diǎn)水位、潮位、流量過(guò)程驗(yàn)證

4 結(jié)論

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

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