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

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

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

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

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

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

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

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

1 概 述

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

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

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

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

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

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

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

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

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

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

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

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

4．1 Schaffranek等人的方法

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

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

式（9）是一個遞歸表達式，對各微段（相鄰節(jié)點間的單元）應用此方程，可得方程

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

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

若河網(wǎng)由M個河段組成，則整體方程組的系數(shù)矩陣是一個4 M×4 M的稀疏矩陣，常采用高斯消去法求解。如引言中所述，我國學者對河網(wǎng)模擬方法作出了重要貢獻，在該類分級解法的基礎上，又提出了三級和四級解法，進一步減小了整體矩陣的大?。ㄆ┤缛壗夥ㄍㄟ^變換，只保留各河段首尾斷面的水位或者流量，將整體矩陣縮小為2 M×2 M）。

4．2 Sen和Garg的方法

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

4．3 JPWSPC法

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

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

圖4 汊點處特征線Fig．4 Characteristic curves at a junction point

根據(jù)這2個關系，代入方程（7），可以構造汊點水位的迭代關系，如式（13）所示，

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

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

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

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

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

5 計算結果及討論

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

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

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

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

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

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

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

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

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

6 結 論

本文比較了最近提出的JPWSPC法和分級解法，在處理緩流河網(wǎng)汊點處回流效應時的表現(xiàn)。比較表明：這2類方法都能處理普適河網(wǎng)，JPWSPC法無需求解整體連接矩陣，同時不會增加每一時間步的迭代次數(shù)，因而節(jié)省了系統(tǒng)內(nèi)存，并且提高了計算效率；河網(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)非恒定流三級聯(lián)合算法［J］．華東水利學院學報，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］，水利學報，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］ 侯 玉，卓建民，鄭國權．河網(wǎng)非恒定流汊點分組解法［J］．水科學進展，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］．水利學報，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é)點水流連接條件處理方法研究［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

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

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