丁 雪

（揚(yáng)州大學(xué) 水利與能源動(dòng)力工程學(xué)院，江蘇 揚(yáng)州 225009）

1 引言

河網(wǎng)非恒定流的數(shù)值計(jì)算，是生產(chǎn)中經(jīng)常要進(jìn)行的工作，具有重要的實(shí)際應(yīng)用價(jià)值，其計(jì)算問題可歸結(jié)為一維圣維南方程組的求解問題。在目前的技術(shù)水平下，得到圣維南方程組的解析解比較困難，只能采用數(shù)值方法求解其近似解。用的較多的數(shù)值計(jì)算方法有特征線法、有限分析法、有限差分法和有限體積法。其中有限差分法又可以分為直接解法、分級(jí)解法、汊點(diǎn)分組解法、矩陣標(biāo)識(shí)法、非線性方法以及汊點(diǎn)水位預(yù)測(cè)修正法。在有限差分法中，Preissmann隱式格式因?yàn)榫哂蟹€(wěn)定性好、計(jì)算速度快的優(yōu)點(diǎn)，所以得到了廣泛應(yīng)用。但是在河網(wǎng)非恒定流水力計(jì)算中，由于天然河道的形狀不規(guī)則，而且河道之間的連接是交叉的，所以形成的系數(shù)矩陣是一個(gè)大型的、不規(guī)則、不對(duì)稱的稀疏矩陣。當(dāng)河網(wǎng)中汊點(diǎn)比較多的時(shí)候，矩陣運(yùn)算會(huì)占用巨大的存儲(chǔ)單元，同時(shí)由于大量的零元素參加運(yùn)算，使得消元速度大為降低，計(jì)算次數(shù)猛增，計(jì)算時(shí)間較長(zhǎng)，精度難以保證。因此河網(wǎng)非恒定流計(jì)算的核心問題是：如何壓縮系數(shù)矩陣的尺度，減少對(duì)計(jì)算機(jī)內(nèi)存的占用量。

文獻(xiàn)[1]對(duì)平原河網(wǎng)水動(dòng)力模型及求解方法進(jìn)行過介紹。本文僅對(duì)有限差分法中的各解法基本思路進(jìn)行討論，比較其優(yōu)缺點(diǎn)。

2 河網(wǎng)方程組的基本形式

2.1 微段方程基本形式

描述一維明渠非恒定流的Saint-Venant方程組是由連續(xù)與運(yùn)動(dòng)方程組構(gòu)成：

式中：A為過水?dāng)嗝婷娣e；Q為流量；q為單位河長(zhǎng)的旁側(cè)入流量，入流為正，出流為負(fù)；y為水位；k為流量模數(shù)；α為動(dòng)量校正系數(shù)；g為重力加速度；t為時(shí)間；x為沿水流方向的距離。

對(duì)于某一河段的任意兩個(gè)相鄰計(jì)算斷面 i，i+1之間的微段，用Preissmann四點(diǎn)隱式差分格式差分逼近Saint-Venant方程組，可得微段差分方程組：

式中：Δy、ΔQ為待求水位、流量的增量；a，b，c，d，e 及， a′， b′， c′， d′， e′是已知水位、流量（包括迭代初值）的函數(shù)。

2.2 邊點(diǎn)方程

邊點(diǎn)方程形式是由邊界控制條件決定的。邊界控制條件有流量控制、水位控制、水位流量關(guān)系控制等多種情況[2]，但相應(yīng)的控制方程可以統(tǒng)一用以下形式表示：

其中a、b、c為按實(shí)際條件導(dǎo)得的系數(shù)和右端項(xiàng)。這一類方程的非零系數(shù)最多只有兩個(gè)，分別位于主對(duì)角線及緊鄰該線（前或后）的位置上，排列也是比較規(guī)整的。

2.3 汊點(diǎn)連接方程

2.3.1 流量銜接條件

進(jìn)出每一汊點(diǎn)的流量必須與該汊點(diǎn)內(nèi)實(shí)際水量的增減率相平衡：

式中：n，m分別表示與某一汊點(diǎn)相連的河段數(shù)和汊點(diǎn)號(hào)，Vm表示m汊點(diǎn)蓄水量，在假設(shè)匯合區(qū)很小，水位變化引起的匯合區(qū)水體積的變化不計(jì)的情況下，可以簡(jiǎn)化成如下形式：

2.3.2 能量銜接條件

在汊點(diǎn)k的第i斷面和第j斷面間建立能量銜接方程時(shí)，如果考慮流速水頭的影響和斷面間能量損失，可以得到的方程(8)，其中， Z （ k, i） ，V（ k, i）分別表示k汊點(diǎn)第i斷面水位和平均流速。

3 求解方法

在選擇差分格式將圣維南方程組進(jìn)行差分離散時(shí)，有顯格式和隱格式兩種不同的方法。顯格式差分法的優(yōu)點(diǎn)是容易理解，方便編制計(jì)算程序，但由于顯式差分的穩(wěn)定需要具備一定的條件，對(duì)計(jì)算時(shí)間步長(zhǎng)要求苛刻，同時(shí)需要采取一定的技巧處理流量和水位之間的不協(xié)調(diào)性[3]，因而很少利用它來求解河網(wǎng)非恒定流。隱式差分法具有穩(wěn)定性好，精度高的優(yōu)點(diǎn)，但是在計(jì)算河網(wǎng)非恒定流時(shí)，如果使用隱式差分法，每進(jìn)行一次迭代就要求解一個(gè)高階代數(shù)方程組，使得隱格式的實(shí)際應(yīng)用價(jià)值降低了。這篇文章中主要介紹隱格式的有限差分法，顯格式方法可參考文獻(xiàn)[4-5]。

隱格式差分法又可分為直接解法和分級(jí)解法兩大類。

3.1 直接解法

直接解法的基本思想是把河道斷面的水力要素作為基本未知量，去直接求解由內(nèi)斷面方程和邊界方程構(gòu)成的方程組。中山大學(xué)數(shù)學(xué)力學(xué)系從河網(wǎng)矩陣的特點(diǎn)出發(fā)，提出了“網(wǎng)河不恒定流隱式方程組的稀疏矩陣解法”，這種方法能夠有效節(jié)省計(jì)算機(jī)內(nèi)存，提高計(jì)算速度[6-7]。

由于直接解法系數(shù)矩陣階數(shù)是河網(wǎng)中選取的計(jì)算斷面數(shù)的兩倍，所以對(duì)于大型河網(wǎng)，幾乎是不可能用直接解法去進(jìn)行求解。

3.2 分級(jí)解法

荷蘭學(xué)者Dronkers首先提出了分級(jí)解法，在這之后又有許多學(xué)者對(duì)此法進(jìn)行了進(jìn)一步的完善。分級(jí)解法的基本思想[8]是先將未知水力要素集中到汊點(diǎn)上，待汊點(diǎn)未知數(shù)求出后，再將各河段當(dāng)作單一河道進(jìn)行求解。實(shí)踐證明，分級(jí)解法與直接解法相比更能節(jié)約計(jì)算機(jī)存儲(chǔ)量，更加實(shí)用。按照汊點(diǎn)上所保留未知數(shù)的個(gè)數(shù)對(duì)分級(jí)解法進(jìn)行分類，如果在汊點(diǎn)上保留水位和流量?jī)蓚€(gè)未知數(shù)，稱為二級(jí)解算法；僅保留水位或流量一個(gè)未知數(shù)，稱為三級(jí)解算法。

3.2.1 二級(jí)解法

二級(jí)解法的基本思想是將所有的邊界方程和河段方程構(gòu)成的二級(jí)連接方程組，可得到河段內(nèi)部斷面的水力要素，然后利用微段方程，即可求出全部?jī)?nèi)部計(jì)算斷面的未知數(shù)[8]。

3.2.2 三級(jí)解法

三級(jí)解法其思想是將問題歸結(jié)于節(jié)點(diǎn)水位（或水位增量）的方程組，再求解節(jié)點(diǎn)間斷面的水位、流量。此法相比于直接解法所需求解的代數(shù)方程組的階數(shù)降得多，并且更加準(zhǔn)確便捷[3]。

3.2.3 四級(jí)解法

四級(jí)解法是從三級(jí)解法的基礎(chǔ)上提出來的。進(jìn)一步從三級(jí)連接方程組中分離出外邊界方程和汊點(diǎn)能量銜接方程，就可以得到以河網(wǎng)內(nèi)部汊點(diǎn)水力要素為基本未知量的四級(jí)連接方程組[8]。四級(jí)解法所用機(jī)時(shí)要比三級(jí)解法直接消元所用的機(jī)時(shí)少，而且能夠更有效的節(jié)省計(jì)算機(jī)存儲(chǔ)量。

文獻(xiàn)[11]對(duì)各種分級(jí)解法的推導(dǎo)和求解做了詳細(xì)介紹。

3.3 矩陣標(biāo)識(shí)法

矩陣標(biāo)識(shí)法是在三級(jí)分解法的基礎(chǔ)上，根據(jù)節(jié)點(diǎn)水位方程系數(shù)矩陣的高稀疏性，對(duì)矩陣非零元素進(jìn)行標(biāo)識(shí)[9]。按照代碼指示，將非零元素用一維數(shù)組存儲(chǔ)，排除零元素，求解時(shí)只對(duì)非零元素運(yùn)算。該方法避免了零元素的處理，并以一維數(shù)組存儲(chǔ)系數(shù)矩陣，提高了檢索效率，整體上使方程組求解的效率大為提高[10]。但該方法未闡明系數(shù)矩陣與河網(wǎng)拓?fù)浣Y(jié)構(gòu)的關(guān)系，矩陣標(biāo)識(shí)與求解之間具有強(qiáng)耦合關(guān)系，與面向?qū)ο蟪绦蛟O(shè)計(jì)提倡對(duì)象之間的低耦合性相左，增加了面向?qū)ο缶幊痰碾y度，也降低了程序的靈活性和可維護(hù)性，限制了方法的應(yīng)用和推廣。

3.4 汊點(diǎn)分組解法

汊點(diǎn)分組解法的特點(diǎn)是能根據(jù)實(shí)際需要，將河網(wǎng)中的汊點(diǎn)分片分組編碼后進(jìn)行計(jì)算，對(duì)汊點(diǎn)方程組的系數(shù)矩陣進(jìn)行壓縮，使得壓縮后的矩陣與分組后每組中的汊點(diǎn)數(shù)的階數(shù)相同[12]。將河網(wǎng)汊點(diǎn)分組編碼后進(jìn)行計(jì)算，使得計(jì)算的工作量大幅度減少，節(jié)約了內(nèi)存，提高了工作效率。

實(shí)際求解時(shí)，識(shí)別各河段與汊點(diǎn)之間的關(guān)系形成汊點(diǎn)方程組是比較困難的，如果將汊點(diǎn)分組，將會(huì)使得問題的復(fù)雜性進(jìn)一步增加。另外汊點(diǎn)分組解法的遞推計(jì)算相當(dāng)復(fù)雜，程序編制計(jì)算難度大。

侯玉為解決上述問題提出了汊點(diǎn)分組新解法，其基本思想是利用矩陣的分塊計(jì)算技術(shù)將一般分級(jí)解法形成的原汊點(diǎn)水位關(guān)系應(yīng)用于汊點(diǎn)分組，從而使遞推過程得到簡(jiǎn)化，這種解法具有較好的通用性[13]。

3.5 非線性方法

徐小明將樹狀河網(wǎng)的松弛迭代法應(yīng)用到環(huán)狀河網(wǎng)，將任意復(fù)雜的河網(wǎng)視為一系列的單一河道，而每一單一河道所形成的線性方程組系數(shù)都是五對(duì)角矩陣[14]。由于松弛迭代法是一條一條河道進(jìn)行求解，所以可不對(duì)河道匯流節(jié)點(diǎn)進(jìn)行編碼、河道斷面也不必全局編碼，而且河道數(shù)目可以靈活增減，矩陣帶寬減少從而提高計(jì)算精度。采用 Newton-Raphson方法直接求解非線性代數(shù)方程組，一般情況下收斂速度快，并且自動(dòng)校正，也就是說前面迭代產(chǎn)生的誤差并不會(huì)一步一步傳遞下去。這種方法的缺點(diǎn)是：對(duì)初始值的要求非常嚴(yán)格，只有迭代初值選取合適的時(shí)候，才能得到滿足一定精度要求的解。

3.6 牛頓迭代解法

追趕法是經(jīng)常用來求解差分格式的方法，在河網(wǎng)計(jì)算中，追趕法在很多情況下不通用。牛頓迭代解法在河網(wǎng)計(jì)算應(yīng)用方面相比于追趕法有其獨(dú)到的優(yōu)勢(shì)，并且具有一定的通用性[15]。河道斷面的編碼可以任意，使得處理一些較復(fù)雜的自然情況變得方便，另外這種方法的計(jì)算速度比較快，較好的滿足了實(shí)踐的要求。

3.7 汊點(diǎn)水位預(yù)測(cè)修正法

利用汊點(diǎn)水位預(yù)測(cè)修正法處理汊點(diǎn)處的回流效應(yīng)，無需建立和求解整體分支河道方程組，則不需要求解大型稀疏矩陣所對(duì)應(yīng)的代數(shù)方程組，因此節(jié)省內(nèi)存空間；每個(gè)分支河道都能獨(dú)立地進(jìn)行計(jì)算，使得河道數(shù)目可以靈活增減；河道及計(jì)算斷面編號(hào)簡(jiǎn)單方便，不需要進(jìn)行節(jié)點(diǎn)優(yōu)化編碼；此法還適用于環(huán)狀和樹狀河網(wǎng)并且對(duì)流動(dòng)方向沒有要求[16]。

文獻(xiàn)[17]給出了汊點(diǎn)水位預(yù)測(cè)修正法（又稱雙松弛迭代法）求解河網(wǎng)水流的具體步驟。

上述這些求解方法在緩流狀態(tài)下一般應(yīng)用效果較好，但是當(dāng)河段中出現(xiàn)急流或激波現(xiàn)象時(shí)，往往無法準(zhǔn)確模擬[18]。后來又有學(xué)者提出，采用具有良好激波捕捉能力的 Godunov格式建立新的一維河網(wǎng)水動(dòng)力數(shù)學(xué)模型。

4 結(jié)語

無論是哪一種數(shù)值解法求解偏微分方程，最終原問題都會(huì)被歸結(jié)為線性方程組的求解。而且求解的精確度受到線性方程組系數(shù)矩陣特性的直接影響。河網(wǎng)非恒定流計(jì)算是非常復(fù)雜的，特別是考慮到實(shí)際應(yīng)用的情況時(shí)，數(shù)值求解計(jì)算方法會(huì)變得復(fù)雜很多，仍然有許多問題需要學(xué)者們進(jìn)一步研究，比如在建立汊點(diǎn)能量銜接方程時(shí)考慮水流運(yùn)動(dòng)及水環(huán)境變化規(guī)律。

[1] 盧士強(qiáng),徐祖信.平原河網(wǎng)水動(dòng)力模型及求解方法探討[J].水資源保護(hù),2003,(3)：5-9.

[2] 焦?jié)櫦t.河口圍墾工程與一二維銜接水沙數(shù)學(xué)模型的研究[D].天津：天津大學(xué)建筑工程院,2005.

[3] 閆秀平,姚慧敏.一維河網(wǎng)非恒定流數(shù)值計(jì)算研究及其應(yīng)用[J].河北水利,2014,(4)：33.

[4] 郭新蕾.河網(wǎng)的一維水動(dòng)力及水質(zhì)分析研究[D].武漢：武漢大學(xué)環(huán)境工程,2005.

[5] 顧元掞,王尚毅,郭傳鎮(zhèn).具有岔道支流的明渠不恒定計(jì)算[J].天津大學(xué)學(xué)報(bào),1995,(1)：54-64.

[6] 李岳生.網(wǎng)河不恒定流隱式方程組系數(shù)矩陣解法[J].中山大學(xué)(自然科學(xué)版),1977,(3)：27-37.

[7] 楊世孝,肖子良.河網(wǎng)不恒定流隱式方程組的稀疏矩陣解法[J].中山大學(xué)學(xué)報(bào)(自然科學(xué)學(xué)報(bào)),1994,(1)：61-68.

[8] 于麗麗.天津市北系河網(wǎng)洪水資源化蓄滯洪區(qū)蓄水研究[D].天津：天津大學(xué)建筑工程學(xué)院,2005.

[9] 李光枳,王船海.大型河網(wǎng)水流模擬的矩陣標(biāo)識(shí)法[J].河海大學(xué)學(xué)報(bào),1995,23(1)：36-43.

[10] 錢真,賈衛(wèi)紅,李世陽.河網(wǎng)水流模擬的矩陣標(biāo)識(shí)法研究[J].人民長(zhǎng)江,2014,45(14)：85-88.

[11] 白玉川,萬春艷,黃本勝,等.河網(wǎng)非恒定流數(shù)值模擬的研究進(jìn)展[J].水利學(xué)報(bào),2000,(12)：43-47.

[12] 李義天.河網(wǎng)非恒定流隱式方程組的汊點(diǎn)分組解法[J].水利學(xué)報(bào),1997,(3)：49-57.

[13] 侯玉,卓建明,鄭國權(quán).河網(wǎng)非恒定流汊點(diǎn)分組解法[J].水科學(xué)進(jìn)展,1999,10(1)：48-52.

[14] 徐小明,何建京,汪徳,等.求解大型河網(wǎng)非恒定流的非線性方法[J].水動(dòng)力學(xué)研究與進(jìn)展,2001,16(1)：18-24.

[15] 張大偉,董增川,李少華.河網(wǎng)非恒定流計(jì)算的牛頓迭代解法[J].西北水利發(fā)電,2004,20(3)：31-33.

[16] 陳永燦,王智勇,朱德軍,等.復(fù)雜河網(wǎng)水動(dòng)力數(shù)值模型[J].水力發(fā)電學(xué)報(bào),2011,22(2)：203-207.

[17] 姜偉,吉慶豐,宋振華,等.南方河網(wǎng)區(qū)引排水過程數(shù)值模擬研究[J].灌溉排水學(xué)報(bào),2013,32(6)：91-95.

[18] 張大偉,權(quán)錦,馬建明.應(yīng)用 Godunov格式模擬復(fù)雜河網(wǎng)明渠水流運(yùn)動(dòng)[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào),2015,23(6)：1088-1096.