朱德軍，陳永燦，劉昭偉，王智勇

（清華大學(xué)水利水電工程系，北京 100084）

河網(wǎng)模擬JPWSPC方法和分級(jí)解法的對(duì)比

朱德軍，陳永燦，劉昭偉，王智勇

（清華大學(xué)水利水電工程系，北京 100084）

比較了汊點(diǎn)水位預(yù)測-校正法（junction－point water stage prediction and correction，JPWSPC）和經(jīng)典的分級(jí)解法，在處理緩流河網(wǎng)汊點(diǎn)處回流效應(yīng)時(shí)的不同，并對(duì)比了它們的計(jì)算效率。應(yīng)用這2類方法時(shí)，Saint-Venant方程組都采用Preissmann格式離散，生成的非線性離散方程用Newton-Raphson方法求解。比較表明：這2類方法都能處理普適河網(wǎng)，JPWSPC法無需求解整體連接矩陣，同時(shí)不會(huì)增加每一時(shí)間步的迭代次數(shù)，因而節(jié)約了系統(tǒng)內(nèi)存，提高了計(jì)算效率；河網(wǎng)中河段數(shù)目越多，JPWSPC法的效率優(yōu)勢越明顯。

河網(wǎng)水動(dòng)力數(shù)值模擬；汊點(diǎn)水位預(yù)測-校正法（JPWSPC）；分級(jí)解法；效率

1 概 述

隱式差分法穩(wěn)定性好，能允許較大的時(shí)間步長，在一維水動(dòng)力數(shù)值模擬中應(yīng)用最為廣泛，應(yīng)用隱式差分法離散控制方程后，將形成離散方程組，在每一時(shí)間步都需要求解此方程組。單一河道離散方程組的系數(shù)矩陣有明顯的帶狀特征，大大節(jié)省了計(jì)算量。當(dāng)眾多河道相互連接形成河網(wǎng)時(shí)，如果流動(dòng)為緩流，由于回水效應(yīng)，河網(wǎng)中各分支河道互相影響，必須對(duì)河網(wǎng)中的所有支流同時(shí)進(jìn)行模擬［1］。與此同時(shí)，河網(wǎng)連接條件破壞了離散矩陣的帶狀特征，使得求解大型河網(wǎng)系數(shù)矩陣所需的內(nèi)存和計(jì)算量都非常巨大。

針對(duì)這一問題，誕生了很多處理方法（詳見文獻(xiàn)［1］和［2］中的綜述），其中分級(jí)解法［1，3］既能求解樹狀河網(wǎng)，又能求解環(huán)狀河網(wǎng)，被大量學(xué)者和工程人員采用，尤其在我國，圍繞此類方法開展了大量研究，并提出了一系列的改進(jìn)形式，如三級(jí)、四級(jí)解法等［4－7］。分級(jí)解法主要分為3步：第一步，利用各河網(wǎng)系數(shù)矩陣的帶狀特性，通過消元，建立各河段首末斷面變量之間的聯(lián)系；第二步，將第一步建立的關(guān)系與汊點(diǎn)處的連接條件以及外邊界條件聯(lián)合組成整體矩陣，求解此整體矩陣，得到各河段首末斷面變量值；第三步：將各河段首末斷面的值，代入第一步變換后的系數(shù)矩陣，得到所有內(nèi)節(jié)點(diǎn)的變量值。分級(jí)解法的第一和第三步利用了系數(shù)矩陣的帶狀特性，第二步中的整體矩陣與整個(gè)河網(wǎng)的離散矩陣相比，尺寸大大降低，因而提高了計(jì)算效率。

Islam等人［1］曾比較了2種分級(jí)解法，我國學(xué)者提出的三級(jí)和四級(jí)等改進(jìn)解法，都可以建立在這類分級(jí)解法之上。但是，呂滿英等人［8］發(fā)現(xiàn)，當(dāng)初始誤差比較大時(shí)分級(jí)解法容易失穩(wěn)。另外，當(dāng)河網(wǎng)中河段數(shù)目很多時(shí)，第二步中的矩陣尺寸依然很大，由于該矩陣是不規(guī)則的稀疏矩陣，求解相當(dāng)耗時(shí)。文獻(xiàn)［2］中，陳永燦和朱德軍等人首次提出汊點(diǎn)水位預(yù)測校正（JPWSPC）法，利用緩流的特征處理汊點(diǎn)處的回水效應(yīng)。該方法具有分級(jí)解法普適性強(qiáng)的優(yōu)點(diǎn)，無需建立和求解總體矩陣。本文將JPWSPC法和Islam等人選取的2個(gè)分級(jí)解法進(jìn)行對(duì)比，重點(diǎn)分析了他們的效率差別。

2 河網(wǎng)控制方程組

河網(wǎng)水動(dòng)力過程由一維圣維南方程組描述，如方程（1）和（2）所示：

式中：A為過水面積；Q為流量；q為側(cè)向入流量；Z為水位；g為重力加速度；Sf為摩阻坡度，Sf＝n2Q｜Q｜／（A2R4／3），R為水力半徑，n為糙率系數(shù)。汊點(diǎn)連接條件如（3）、（4）所示：式中下標(biāo)i和o分別代表流入或流出汊點(diǎn)的河道斷面變量值。

3 河網(wǎng)控制方程組數(shù)值離散

采用Preissmann隱式差分格式離散方程（1）、（2），并利用Newton-Raphson方法求解離散形成的非線性方程組，得到如下方程：

其中，F(xiàn)C，F(xiàn)M分別表示連續(xù)和動(dòng)量方程，RFC，RFM分別表示其余量，Δ表示連續(xù)兩個(gè)Newton-Raphson迭代步的變量增量。同樣，采用Newton-Raphson法求解（3）、（4）得：

式中：B表示渠道水面寬度；f和g分別表示方程（3）和（4）的左邊項(xiàng)的余量。若汊點(diǎn)處有n個(gè)分支河道交匯，則有一個(gè)形如（7）式的方程和n－1個(gè)形如（8）式的方程。由這2個(gè)方程可見，因?yàn)榛厮?yīng)，不同分支河道的變量互相聯(lián)系，河網(wǎng)離散方程組的系數(shù)矩陣不再是單一河道的五對(duì)角矩陣，這就是用隱式差分法求解緩流河網(wǎng)的難點(diǎn)所在。

單一河道的兩端都是外邊界，其求解方法已很成熟，而且非本文的重點(diǎn)，故不再贅述。不失一般性，取河段一端是汊點(diǎn)，另一端是外邊界的情況進(jìn)行分析。為了便于比較，本文將各種方法都用系數(shù)矩陣的形式表達(dá)，方程（5）至（8）以及外邊界條件表示成矩陣形式如圖1所示。其中EBC和IBC分別表示外和內(nèi)邊界條件，零元素沒有標(biāo)出，＃代表非零元素，既可能是零元素也可能是非零元素，可見，在內(nèi)邊界處，系數(shù)矩陣的帶狀特征被破壞。

4 河網(wǎng)處理方法

4．1 Schaffranek等人的方法

該方法由Schaffranek等人于1981年提出（見文獻(xiàn)［1］），首先將方程（5）和（6）寫成矩陣形式如（9）所示：

圖1 原始系數(shù)矩陣示意圖Fig．1 O riginal coefficientm atrix

式（9）是一個(gè)遞歸表達(dá)式，對(duì)各微段（相鄰節(jié)點(diǎn)間的單元）應(yīng)用此方程，可得方程

式中：m為河段中的節(jié)點(diǎn)數(shù)；E1＝M1；F1＝N1；Ei＋1＝Mi＋1Ei；Fi＋1＝Mi＋1Fi＋Ni＋1。實(shí)際上，由方程（5）和（6）推導(dǎo)方程（10）的過程，等價(jià)于將圖1所示系數(shù)矩陣經(jīng)過行變換得到圖2所示矩陣的過程。圖2所示矩陣的第一行，以及最后三行中只含有首尾節(jié)點(diǎn)的變量，由此可得整體河段方程組。

圖2 Schaffranek等人方法系數(shù)矩陣示意圖Fig．2 Coefficientmatrix of themethod by Schaffranek et al．

若河網(wǎng)由M個(gè)河段組成，則整體方程組的系數(shù)矩陣是一個(gè)4 M×4 M的稀疏矩陣，常采用高斯消去法求解。如引言中所述，我國學(xué)者對(duì)河網(wǎng)模擬方法作出了重要貢獻(xiàn)，在該類分級(jí)解法的基礎(chǔ)上，又提出了三級(jí)和四級(jí)解法，進(jìn)一步減小了整體矩陣的大小（譬如三級(jí)解法通過變換，只保留各河段首尾斷面的水位或者流量，將整體矩陣縮小為2 M×2 M）。

4．2 Sen和Garg的方法

Sen和Garg的方法［1，3］，其思想與Schaffranek等人的方法類似，也是首先從各河段中，提取出首尾節(jié)點(diǎn)的變量，由這些變量組成整體河段方程組，求解河段方程組，然后代入求解各內(nèi)節(jié)點(diǎn)的變量值。不同的是，他們首先由圖1所示矩陣，變形為圖3所示矩陣，顯然前三行和最后一行中只含有首尾節(jié)點(diǎn)的變量，由此可得整體河段方程組，同樣是一個(gè)4 M× 4 M的稀疏矩陣。

4．3 JPWSPC法

上述2種方法處理汊點(diǎn)連接條件時(shí)，采用的都是純粹的數(shù)學(xué)方法。陳永燦和朱德軍等人［2］根據(jù)緩流的特征，提出了JPWSPC法，用以處理汊點(diǎn)處的回水效應(yīng)。根據(jù)非恒定漸變緩流的特點(diǎn)，流入和流出汊點(diǎn)斷面的流量受汊點(diǎn)水位的影響，若規(guī)定流入為正，流出為負(fù)，當(dāng)汊點(diǎn)水位過高時(shí)，汊點(diǎn)處凈流量為負(fù)，反之汊點(diǎn)處凈流量為正。根據(jù)這一特點(diǎn)，在一次時(shí)間步進(jìn)過程中，首先采用一個(gè)預(yù)測步，預(yù)測各汊點(diǎn)水位，再用若干校正步，使汊點(diǎn)處的條件滿足方程（3）和（4）的要求，這就是JPWSPC法。

如圖4所示，A點(diǎn)代表一汊點(diǎn)，其坐標(biāo)為x0，UA和AD分別代表匯于汊點(diǎn)A的兩分支河道，水流方向如圖中箭頭所示，λ＋和λ－分別為流經(jīng)點(diǎn)（x0，t0＋Δt）的正負(fù)特征線，根據(jù)圣維南方程組的性質(zhì)，在分支河道UA和AD中，水深和流量分別近似滿足式（11）和（12）所示關(guān)系。

圖4 汊點(diǎn)處特征線Fig．4 Characteristic curves at a junction point

根據(jù)這2個(gè)關(guān)系，代入方程（7），可以構(gòu)造汊點(diǎn)水位的迭代關(guān)系，如式（13）所示，

為了簡單起見，引入變量AC如式（15）所示，

其中：α為可調(diào)整的常數(shù)，反映式（11）和（12）推導(dǎo)過程中所作假設(shè)的影響，根據(jù)經(jīng)驗(yàn)，α可以取為1．0～2．0，較大的α值有利于計(jì)算穩(wěn)定，較小的α值有利于提高收斂速度。式（14）進(jìn)一步變形為

將式（16）代入汊點(diǎn)處的內(nèi)邊界條件，圖1所示矩陣將變形為圖5所示矩陣。

圖5 JPWSPC法系數(shù)矩陣示意圖Fig．5 Coefficientmatrix of the JPWSPC method

這樣，通過JPWSPC法，實(shí)現(xiàn)了汊點(diǎn)處的解耦，各河段的變量形式上不再互相聯(lián)系。在每一New-ton-Raphson迭代步，河網(wǎng)的離散矩陣都由彼此獨(dú)立的五對(duì)角矩陣組成，各五對(duì)角矩陣可以獨(dú)立求解。顯然，應(yīng)用JPWSPC方法，求解過程非常簡潔，除了將方程（7）和（8）用方程（16）代替外，沒有任何額外的步驟，易于程序?qū)崿F(xiàn)，而且不需要求解不規(guī)則的稀疏整體連接矩陣。

5 計(jì)算結(jié)果及討論

假想渠網(wǎng)如圖6所示，渠道的參數(shù)見表1。該渠網(wǎng)既包含樹狀又包含環(huán)狀組成渠道，是一個(gè)理想的檢查河網(wǎng)模擬方法的算例，該算例也曾被Sen和Garg采用［3］。在渠道8-1的外端給定流量條件，其他渠道外邊界固定水深為5 m。本研究采用分級(jí)解法和JPWSPC方法取得的結(jié)果幾乎完全一致，所以只畫出一條曲線，同時(shí)，畫出了文獻(xiàn)［3］的結(jié)果，與本文結(jié)果進(jìn)行比較，由圖7可見，結(jié)果幾乎完全一致。

表1 假想渠網(wǎng)參數(shù)［3］Table 1 Parameters of the hypothetic channel network［3］

采用分級(jí)解法時(shí)，連接矩陣采用選主元高斯消去法求解；所有方法中，采用同樣的收斂判據(jù)，即｜ΔQ／Q｜≤0．001和｜ΔA／A｜≤0．001，同時(shí)，采用同樣的時(shí)間步長Δt＝180 s。由圖7結(jié)果可見，渠道7中產(chǎn)生反方向流動(dòng)，這一相當(dāng)復(fù)雜的流態(tài)，所有的方法都能準(zhǔn)確模擬。

圖6 假想渠網(wǎng)示意圖［3］Fig．6 The hypothetic channel network［3］

表2 各種方法效率對(duì)比（應(yīng)用于圖6河網(wǎng)）Table 2 Comparison of the efficiencies of every method used in the network depicted in Fig．6

圖7 假想渠網(wǎng)邊界條件和計(jì)算結(jié)果Fig．7 Boundary conditions and calculated results for the hypothetic channel network

為了比較各方法的計(jì)算效率和對(duì)內(nèi)存的要求，表2比較了各方法產(chǎn)生的整體連接矩陣尺寸、迭代次數(shù)以及模擬時(shí)間。結(jié)果顯示，因?yàn)椴恍枰⒄w連接矩陣，JPWSPC法節(jié)約了內(nèi)存，而且減少了每一迭代步的計(jì)算時(shí)間。表2同時(shí)顯示，JPWSPC法并不明顯增加迭代次數(shù)，所以該方法比分級(jí)解法快大概28%。

顯然，整體連接矩陣的大小，與河網(wǎng)中河段數(shù)目有關(guān)，河網(wǎng)中河段數(shù)目越多，JPWSPC法所節(jié)省的內(nèi)存就越大。為了進(jìn)一步比較，將河網(wǎng)的節(jié)點(diǎn)加密，使微元長度是表1中的一半，相應(yīng)地，時(shí)間步長也縮小一半，Δt＝90 s，比較結(jié)果如表3所示。表3中的結(jié)論與表2類似，JPWSPC法節(jié)省了系統(tǒng)內(nèi)存，減少了模擬時(shí)間，同時(shí)迭代步數(shù)沒有明顯的增加。表3還顯示，因?yàn)楣?jié)點(diǎn)加密，而河段數(shù)目不變，JPWSPC法相對(duì)分級(jí)解法的效率優(yōu)勢有所減弱，比分級(jí)解法只快了26%。

表3 各種方法效率對(duì)比（應(yīng)用于加密的圖6河網(wǎng)）Table 3 Com parison of the efficiencies of every method used in the refined Fig．6 network

6 結(jié) 論

本文比較了最近提出的JPWSPC法和分級(jí)解法，在處理緩流河網(wǎng)汊點(diǎn)處回流效應(yīng)時(shí)的表現(xiàn)。比較表明：這2類方法都能處理普適河網(wǎng)，JPWSPC法無需求解整體連接矩陣，同時(shí)不會(huì)增加每一時(shí)間步的迭代次數(shù)，因而節(jié)省了系統(tǒng)內(nèi)存，并且提高了計(jì)算效率；河網(wǎng)中河段數(shù)目越多，JPWSPC法的效率優(yōu)勢越明顯。

［1］ ISLAM A，RAGHUWANSHI N S，SINGH R，et al．Comparison of gradually varied flow computation algo-rithms for open－channel network［J］．J．Irrig．Drain．Eng．，2005，131（5）：457－465．

［2］ ZHU De-jun．，CHEN Yong-can，WANG Z Y．A novel method for gradually varied subcritical flow simulation in general channel networks［C］∥33rd IAHR Congress．Vancouver，Canada，2009：6327－6335．

［3］ SEN D J，GARG N K．Efficient algorithm for gradually varied flows in channel networks［J］．J．Irrig．Drain．Eng．，2002，128（6）：351－357．

［4］ 張二駿，張東生，李 挺．河網(wǎng)非恒定流三級(jí)聯(lián)合算法［J］．華東水利學(xué)院學(xué)報(bào)，1982，（1）：1－13．（ZHANG Er-jun，ZHANG Dong-sheng，LITing．Three Steps Method to compute unsteady flow for river networks［J］．Journal of East China Institute of Water Re-sources，1982，（1）：1－13．（in Chinese））

［5］ 吳壽紅．河網(wǎng)非恒定流四級(jí)解算法［J］，水利學(xué)報(bào)，1985，（8）：42－50．（WU Shou-hong．Four Steps Method to compute unsteady flow for river networks［J］．Journal of Hydraulic Engineering，1985，（8）：42－50．（in Chinese））

［6］ 侯 玉，卓建民，鄭國權(quán)．河網(wǎng)非恒定流汊點(diǎn)分組解法［J］．水科學(xué)進(jìn)展，1999，10（1）：48－52．（HOU Yu，ZHUO Jian-min，ZHENGGuo-quan．Unsteady flow computation using junctions group method in river net-work［J］．Advances in Water Science，1999，10（1）：48－52．（in Chinese））

［7］ 楊開林，白正裕．調(diào)水渠網(wǎng)非恒定流的線性變換求解方法［J］．水利學(xué)報(bào)，2004，（3）：35－41．（YANG Kai-lin，BAI Zheng-yu．Linear transformation method for calculation of unsteady flow in water diversion net-works［J］．Journal of Hydraulic Engineering，2004，（3）：35－41．（in Chinese））

［8］ 呂滿英，江 洧，詹杰民．河網(wǎng)節(jié)點(diǎn)水流連接條件處理方法研究［J］．人民黃河，2007，29（3）：31－32．（LV Man-ying，JIANG Wei，ZHAN Jie-min．Study of treatment for linking conditions at junction points of river networks［J］．Yellow River，2007，29（3）：31－32．（in Chinese））

（編輯：周曉雁）

Comparison of JPWSPC M ethod and Three-phase Algorithm s for Hydrodynam ic Simulation in General Open-channel Networks

ZHU De-jun，CHEN Yong-can，LIU Zhao-wei，WANG Zhi-yong
（Department of Hydraulic Engineering，Tsinghua University，Beijing 100084，China）

This paper compares a newly proposed method，namely the junction-pointwater stage prediction and cor-rection（JPWSPC）method，for gradually varied subcritical flow simulation in general open-channel networks with conventional three-phase algorithms．Emphasis is on their difference in backwater effect treatment and efficiency．When the two classes of algorithms are implemented，the Saint-Venant equations are discretized using the Preiss-mann scheme and the resulting nonlinear equation system is solved using the Newton-Raphsonmethod．The compar-ison shows that the two classes of methods are applicable in general networks，that，furthermore，the JPWSPC method saves computermemory and computation time by avoiding to establish and solve the globalmatrix without increasing the required number of iterations per time step，and that the advantage of the JPWSPCmethod over the three-phase algorithms increases with increase of the number of junction points．

hydrodynamic simulation in channel networks；junction-point water stage prediction and correction（JPWSPC）method；three-phase algorithms；efficiency

TV131

A

1001－5485（2010）10－0039－05

2010-08-18

國家水體污染控制與治理科技重大專項(xiàng)（2008ZX07207－010－05）；國家自然科學(xué)基金資助項(xiàng)目（50779026）

朱德軍（1980-），男，江蘇鹽城人，博士后，主要從事水環(huán)境和水生態(tài)問題研究，（電話）010-62772255（電子信箱）zhudj02＠m(xù)ails．tsinghua．edu．cn。