周 領(lǐng)，胡安妮，吳金遠

泵站水力瞬變的二階Godunov格式模型構(gòu)建

周 領(lǐng)1,2，胡安妮1，吳金遠1

（1. 河海大學(xué)水利水電學(xué)院，南京 210098；2. 河海大學(xué)長江保護與綠色發(fā)展研究院，南京 210098）

為了提高復(fù)雜泵站系統(tǒng)水力瞬變數(shù)值模擬的高效性和穩(wěn)定性，該研究基于泵站系統(tǒng)水力瞬變問題，建立有限體積法Godunov格式的數(shù)學(xué)模型，對簡單管道系統(tǒng)和復(fù)雜泵站系統(tǒng)進行模擬研究。與常用的特征線法求解泵站水力模型方程不同，該模型引進有限體積法二階Godunov格式對模型進行離散，用Riemann求解器對離散通量進行求解。使用MUSCL-Hancock方法進行界面數(shù)值重構(gòu)，采用MINMOD斜率限制器避免虛假震蕩。提出雙虛擬單元邊界處理方法，實現(xiàn)計算區(qū)域與邊界同時達到二階精度。將所建模型計算結(jié)果與精確解、經(jīng)典算例數(shù)據(jù)進行對比，并針對庫朗數(shù)取值和計算網(wǎng)格數(shù)進行敏感性分析。結(jié)果表明：所建模型模擬結(jié)果與精確解、經(jīng)典算例數(shù)據(jù)吻合較好；與特征線法相比，二階Godunov格式更加準確、穩(wěn)定且高效。對于簡單管道系統(tǒng)，特征線法計算耗時0.227 s，二階Godunov格式計算耗時0.017 s。對于實際泵站系統(tǒng)，由于存在多特性的管道結(jié)構(gòu)，二階Godunov格式模擬時需稍微降低庫朗數(shù)。而采用特征線法進行泵站水力過渡過程計算時，若不調(diào)整管道長度或者波速，管道中庫朗數(shù)會小于1，在該文算例中，庫朗數(shù)為0.72～0.76，模擬計算結(jié)果偏差很大。所以需要調(diào)整局部管道長度或波速，以達到庫朗數(shù)為1的條件，這樣處理因改變管道特性而引入計算誤差。綜上，二階Godunov格式模擬方法可以更有效提高傳統(tǒng)泵站系統(tǒng)水力瞬變模擬的高效性、穩(wěn)定性以及準確性。

泵站；模型；有限體積法；二階Godunov格式；水力過渡過程；特征線法；一階Godunov格式

0 引 言

泵站輸水系統(tǒng)在啟、停機過程因流量變化而發(fā)生水力瞬變（或水錘），如果水力控制元件操作或水錘防護設(shè)施設(shè)置不當(dāng)，會引起異常的壓力波動，從而導(dǎo)致機組損壞、爆管等重大安全事故。準確模擬泵站系統(tǒng)水力瞬變對保證泵站系統(tǒng)安全穩(wěn)定運行、實現(xiàn)智慧精準化調(diào)度運行極為重要[1-4]。

常采用特征線法（Method of Charicteristic, MOC）對泵站系統(tǒng)水力瞬變進行建模模擬[5]。對于簡單管道系統(tǒng)，該方法具有簡單、準確的優(yōu)勢。然而，在實際工程建設(shè)中，泵站系統(tǒng)會存在復(fù)雜多特性管道結(jié)構(gòu)，水力計算時，為了滿足庫朗數(shù)為1的穩(wěn)定計算條件，MOC需要改變局部管段波速、網(wǎng)格長度（管道長度），產(chǎn)生計算誤差[5]。

近年來，有限體積法逐漸被用于有壓管道水力瞬變計算。Guinot等[6]將有限體積法一階和二階格式數(shù)值解與精確解進行比對，指出有限體積法可以準確模擬簡單管路中水錘壓力波動。Zhao等[7]基于Godunov格式（有限體積法的一種數(shù)值格式）與Riemann求解器，得到了一階和二階的水錘求解格式，并指出Godunov格式可以用于水錘計算，具有較好的準確性和穩(wěn)定性，但是在邊界處理時需要與特征線法相結(jié)合。趙越等[8]針對有壓管道內(nèi)水力瞬變問題，采用有限體積法Godunov格式進行一維數(shù)學(xué)建模，發(fā)現(xiàn)一階Godunov格式與特征線法的計算精度一致，庫朗數(shù)變化會影響計算精度。Zhou等[9]提出一種基于Godunov格式的二階顯式有限體積法，并將其應(yīng)用基于特征線法的經(jīng)典離散空穴模型，發(fā)現(xiàn)有限體積法模擬值與試驗數(shù)據(jù)吻合較好，并且可以避免虛假數(shù)值震蕩。上述研究主要針對簡單管道水力瞬變，而實際泵站系統(tǒng)含水泵、閥門、各種特性管道等，水力計算更為復(fù)雜。

本文采用有限體積法二階Godunov格式[10-16]對泵站系統(tǒng)水力過渡過程進行建模，提出將虛擬單元體與水庫、水泵等邊界元件相結(jié)合，通過通量求解模擬管道水錘變化。詳細分析了庫朗數(shù)和計算網(wǎng)格數(shù)對MOC方法和有限體積法Godunov格式的數(shù)值耗散影響，并采用Chaudhry泵站系統(tǒng)經(jīng)典算例數(shù)據(jù)[17]與本文有限體積二階Godunov格式計算值進行對比。針對某實際泵站系統(tǒng)，對比分析MOC方法和二階Godunov格式模擬結(jié)果。

1 數(shù)學(xué)模型

1.1 水錘基本控制方程

泵站輸水系統(tǒng)中，管道內(nèi)水錘的基本控制方程[18-22]可由以下連續(xù)方程和動量方程進行描述：

式中是測壓管水頭，m；為沿管軸線方向的長度，m；為時間，s；是水錘波速，m/s；是重力加速度，m/s2；是斷面平均流速，m/s；是管道直徑，m；是按照Darcy-Weisbach方程求得的摩阻系數(shù)。

將水錘基本控制方程（1）～（2）以矩陣形式表示：

采用黎曼求解格式進行求解，將式(3)可改寫為常系數(shù)雙曲方程（4）：

1.2 水泵控制方程

1）機組轉(zhuǎn)動方程

式中是機組轉(zhuǎn)動部分加上進入該部分液體的轉(zhuǎn)化慣量，kg·m2；M是電機轉(zhuǎn)矩，N·m；是泵的軸力矩，N·m。

根據(jù)機組轉(zhuǎn)動方程式，用泰勒展開作近似變換得到泵的轉(zhuǎn)速平衡方程為

式中、是本時刻機組的轉(zhuǎn)速和力矩相對值；0、0是上一時刻機組的轉(zhuǎn)速和力矩相對值；T為機組的慣性時間常數(shù)（s），表示機組在額定軸力矩作用下，額定轉(zhuǎn)速減小到0所需要的時間，表征機組慣性的大小。

2）水頭平衡方程

聯(lián)立機組轉(zhuǎn)動方程（6）和水頭平衡方程（8），利用牛頓迭代法可求得水力瞬變過程中每一時刻機組的水頭、流量、轉(zhuǎn)速、力矩等參數(shù)。

1.3 Godunov求解格式

1.3.1 通量計算

有限體積法將管道離散成若干單元體，每個控制單元內(nèi)的數(shù)值是連續(xù)的，但是在兩個單元控制單元相交的各個界面處的數(shù)值是間斷的。采用Godunov格式進行求解時，可通過黎曼問題計算出每個界面上的離散通量，具體如下：

1.3.2 一階Godunov格式

一階Godunov格式通過計算單元內(nèi)變量的平均值來表示界面左右的值：

方程（12）中邊界通量由黎曼解確定，先根據(jù)黎曼不變量方程對邊界處的離散通量進行計算，然后再計算控制單元內(nèi)的離散通量。通過虛擬邊界處理方法，對計算區(qū)域內(nèi)邊界和控制單元的離散通量統(tǒng)一模擬。

本文為統(tǒng)一計算整個計算區(qū)域的通量值，分別在管道上游邊界和管道下游邊界加入1個虛擬單元0和虛擬單元+1。

由式（12）以及負向特征線方程可推導(dǎo)出如下的黎曼不變量方程：

由式（12）以及正向特征線方程可推導(dǎo)出如下的黎曼不變量方程：

1.3.3 二階Godunov格式

1）數(shù)據(jù)重構(gòu)

將式（18）代入（17）可以得到：

2）推進時間計算

3）Riemann問題的近似求解

將式（23）～（24）代入式（12）即可得到具有二階精度的Godunov格式每個計算單元界面+1/2處的數(shù)值通量。

1.3.4 虛擬邊界

一階Godunov格式為實現(xiàn)全瞬變區(qū)域內(nèi)的統(tǒng)一模擬，在上下游邊界處各加入一個虛擬單元。由于計算精度的提高，二階Godunov格式需在上游邊界處左邊和下游邊界處的右邊分別加入2個虛擬單元0、-1和虛擬單元+1、+2。

對于上游水庫：

對于下游水庫：

對于泵站機組，將式（29）～（32）代入式（8）中即可求解虛擬邊界下的水頭平衡方程：

1.3.5 時間積分

一階Godunov格式通過歐拉法進行，二階Godunov格式采用具有二階精度龍格庫塔法進行計算，如下：

2 算例分析

2.1 模型驗證

2.1.1 簡單管道水錘問題

設(shè)置一上游為水庫，下游為閥門的簡單管道。管道長20.0 m，水錘波速為1 319 m/s，上游水庫水位為30.0 m，初始流速為0.1 m/s，重力加速度為9.8 m/s2，計算網(wǎng)格數(shù)為16，下游閥門設(shè)置為瞬時關(guān)閉，光滑管道（無摩阻）。圖1和圖2中精確解是根據(jù)水錘基本控制方程[5]求解得到，在上文式（1）～（2）中，無摩阻項時，基本方程變成線性偏微分方程，理論上解析解即精確解[23-27]。則結(jié)果中的所有壓力衰減均是由于數(shù)值耗散引起的。

圖1和圖2給出了庫朗數(shù)=0.2、0.6、1.0時，MOC、一階Godunov格式（First-order Godunov scheme, 1st-order GTS）、二階Godunov格式（Second-order Godunov scheme, 2nd-order GTS）模擬閥門處水錘計算結(jié)果。

圖1 不同庫朗數(shù)條件下水錘壓力計算結(jié)果（計算網(wǎng)格數(shù)Ns=16）

如圖1a所示，在=1條件下，MOC方法與一階Godunov格式模擬結(jié)果與精確值完全吻合。當(dāng)1時，MOC方法與一階Godunov格式模擬結(jié)果精度完全一致，且越小，數(shù)值衰減越嚴重。從數(shù)值模型方面分析，兩種方法是相似的。MOC模型可寫成如下形式：

式中下標(biāo)代表待求節(jié)點，-1代表左邊節(jié)點，+1代表左邊節(jié)點。

式（37）～（39）也可寫作式（40）：

如圖1b所示，在=1的條件下，二階Godunov格式方法模擬結(jié)果與精確值完全吻合。在<1時，與MOC和一階Godunov格式方法相比，二階Godunov格式能有效抑制數(shù)值耗散。

此外，考慮計算網(wǎng)格數(shù)對數(shù)值模擬結(jié)果的影響，如圖2所示，增加計算網(wǎng)格數(shù)量，使=256，對比不同庫朗數(shù)條件下MOC方法和二階Godunov格式數(shù)值耗散程度。

圖2 不同庫朗數(shù)條件下水錘壓力計算結(jié)果（計算網(wǎng)格數(shù)Ns=256）

如圖2所示，對于MOC方法，適當(dāng)增加計算網(wǎng)格數(shù)對計算結(jié)果影響很大，網(wǎng)格數(shù)的增加會使數(shù)值耗散減小。對于二階Godunov格式，計算網(wǎng)格數(shù)的影響較小。同樣庫朗數(shù)情況下，二階Godunov格式下稀疏網(wǎng)格（=16）的模擬結(jié)果與MOC方法網(wǎng)格加密后（=256）的模擬結(jié)果基本一致。在相同精度要求下，=256時，MOC的計算耗時0.227 s，二階Godunov格式（=16）計算耗時0.017 s。這說明二階Godunov格式效率更高。

2.1.2 泵站系統(tǒng)水力瞬變模擬驗證

為驗證水泵模型的準確性，本文模型計算值與Chaudhry[17]經(jīng)典數(shù)據(jù)（被學(xué)者們廣泛用于仿真模型驗證）進行對比。水泵系統(tǒng)包括2臺水泵，閥門，上下游水庫。水泵采用并聯(lián)，水泵直接與上游水庫相連。2臺泵的參數(shù)一致：水泵中液體的綜合轉(zhuǎn)動慣量2=16.85，額定效率η=0.84，額定轉(zhuǎn)速N=1 100 r/min，額定流量Q=0.25 m3/s，額定揚程H=60 m。管道參數(shù)如表1所示。

表1 Chaudhry算例的管道參數(shù)[17]

采用特征線法對該系統(tǒng)進行水力分析，水泵轉(zhuǎn)速、流量以及水泵出口水頭計算結(jié)果如表2所示。

如表2所示，本文提出的二階Godunov格式和已有常用的特征線法進行水泵失電工況計算。與Chaudhry[17]計算結(jié)果對比（圖3），水泵流量峰值均為0.25 m3/s；水泵轉(zhuǎn)速峰值均為1 100 r/min；水泵揚程峰值均為87 m。泵轉(zhuǎn)速、流量、壓力的最小值也基本吻合。因此，本文所建立的二階Godunov格式模型能準確模擬泵站系統(tǒng)水力瞬變。

2.2 復(fù)雜泵站系統(tǒng)瞬態(tài)計算

2.2.1 泵站參數(shù)

某泵站系統(tǒng)設(shè)置4臺設(shè)計軸功率為3 500 kW的水泵機組，采用“2用2備”運行方式，水庫1水位6.4 m，水庫2水位3.15 m。2臺泵的參數(shù)一致：設(shè)計揚程=52.7 m，設(shè)計流量=5.0 m3/s，額定轉(zhuǎn)速=500 r/min，轉(zhuǎn)動慣量2=14.5 t·m2，設(shè)計軸功率=3 500 kW。布置方式如圖4所示，管道參數(shù)如表3所示。

圖3 Chaudhry泵站算例計算結(jié)果

圖4 某泵站系統(tǒng)示意圖

表3 泵站系統(tǒng)管道參數(shù)

2.2.2 計算結(jié)果

針對泵機組失電工況進行計算分析，采用兩階段關(guān)閥規(guī)律[28-30]，快關(guān)時間30 s，總關(guān)閥時間90 s，拐點處閥門開度為0.166 7。實際泵站系統(tǒng)水力計算時，MOC一般采用調(diào)整各管道水錘波速，以保證庫朗數(shù)為1（計算方案1：Case1）；另外也可采用保持各管道水錘波速不變，庫朗數(shù)小于1（計算方案2：Case2）。

本文主要計算庫朗數(shù)變化對MOC方法、二階Godunov格式數(shù)值耗散的影響，以及2種計算方案對計算結(jié)果的影響，相應(yīng)庫朗數(shù)取值見表4。

表4 兩種計算方案Cr取值

如表5所示，方案1中MOC方法對管道波速進行調(diào)整，二階Godunov格式方法管道波速取原值。方案1管道系統(tǒng)參數(shù)見表5。

2種方法模擬泵站系統(tǒng)發(fā)生失電后120 s水力瞬變過程。整個計算時段內(nèi)水泵最大反轉(zhuǎn)轉(zhuǎn)速、最大倒流流量、水泵出口閥后最小壓力、水泵出口閥后最大壓力見表6所示。

表5 計算方案1中管道系統(tǒng)參數(shù)處理

表6 兩種計算方案下不同計算方法的計算結(jié)果

如表6所示，在實際泵站系統(tǒng)中，庫朗數(shù)條件的變化對于二階Godunov格式模擬結(jié)果影響很小，水泵最大反轉(zhuǎn)轉(zhuǎn)速穩(wěn)定在?417 r/min，水泵最大倒流流量穩(wěn)定在?3.27 m3/s，揚程峰值變化范圍很小，這說明二階Godunov格式模擬結(jié)果穩(wěn)定，這與簡單管道結(jié)論一致。然而，對于MOC方法，方案1的管道波速經(jīng)調(diào)整后，均高于實際管道波速，故水泵出口閥后壓力峰值增大，而方案2采用管道實際波速，揚程峰值減小3.13 m，倒流流量變化至?3.88 m3/s，庫朗數(shù)條件的變化使水泵最大反轉(zhuǎn)轉(zhuǎn)速、最大倒流流量、水泵出口閥后最大、最小壓力的峰值和周期均出現(xiàn)較大誤差。因此對于水泵系統(tǒng)水力瞬變模擬，本文采用的二階Godunov格式更具有穩(wěn)定性、準確性。機組轉(zhuǎn)速、流量、水泵出口閥后壓力變化過程線見圖5。

圖5 泵站系統(tǒng)計算結(jié)果

為考慮網(wǎng)格數(shù)對計算結(jié)果的影響，在上文2個算例基礎(chǔ)上設(shè)置方案3和方案 4，參數(shù)設(shè)置如表7所示。

表7 方案3、4的參數(shù)設(shè)置

如圖6a所示，對于MOC方法，庫朗數(shù)小于1時，增加網(wǎng)格數(shù)可以適當(dāng)提高計算結(jié)果精確度，但是計算耗時將成倍增加（MOC中Case2和Case4的計算耗時分別為0.063和0.094 s）。如圖6b所示，對于二階Godunov格式，庫朗數(shù)等于或小于1，網(wǎng)格數(shù)對計算結(jié)果精確度影響均很小，但計算耗時會增大（二階Godunov中Case2和Case4的計算耗時分別為0.077 和0.109 s）。因此，當(dāng)模擬水泵系統(tǒng)瞬變過程面臨庫朗數(shù)小于1情況，為達到相同精度要求，MOC需要更精細網(wǎng)格，其計算耗時將大于二階Godunov格式。

圖6 泵站系統(tǒng)水泵轉(zhuǎn)速模擬結(jié)果

3 結(jié) 論

本文采用有限體積法二階Godunov格式對泵站系統(tǒng)水力瞬變進行建模模擬，管道邊界采用虛擬邊界的處理方法，水力元件邊界采用傳統(tǒng)的特征線法進行處理。將數(shù)值模擬結(jié)果與MOC計算所得結(jié)果進行對比分析。主要結(jié)論：

1）本文建立的二階Godunov格式可以準確模擬水泵系統(tǒng)水力瞬變過程中水泵轉(zhuǎn)速、流量、壓力等參數(shù)變化過程線。通過與已有文獻計算結(jié)果對比可知，本文建立模型模擬結(jié)果精確，水泵流量峰值均為0.25 m3/s；水泵轉(zhuǎn)速峰值均為1 100 r/min；水泵揚程峰值均為87 m。泵轉(zhuǎn)速、流量、壓力的最小值也基本吻合。

2）與MOC模型相比，二階Godunov模型具有更好的穩(wěn)定性、準確性以及高效性。當(dāng)庫朗數(shù)＜1時，二階Godunov格式能保持較高的計算精度、結(jié)果穩(wěn)定性。簡單管道系統(tǒng)中，在相同精度要求下，同一臺計算機中特征線法計算耗時是0.227 s，二階Godunov格式計算耗時是0.017 s。復(fù)雜泵站系統(tǒng)，相同精度要求下，MOC計算耗時在二階Godunov格式基礎(chǔ)上成倍增加，二階Godunov格式計算效率更高。

3）在實際泵站系統(tǒng)中，為避免數(shù)值耗散，MOC模型需要對管道波速和管道分段進行前處理，例如調(diào)整波速或者網(wǎng)格長度，或者直接忽略短管，這樣處理不僅繁瑣，而且因改變管道特性而引入計算誤差。對于有限體積法二階Godunov格式，管道波速或者管道分段無需調(diào)整，因為庫朗數(shù)變化對于二階Godunov格式模擬結(jié)果影響很小。在本文實際算例中，水泵最大反轉(zhuǎn)轉(zhuǎn)速穩(wěn)定在?417 r/min，水泵最大倒流流量穩(wěn)定在?3.27 m3/s，揚程峰值變化范圍很小，這說明二階Godunov格式模擬結(jié)果穩(wěn)定。既簡單方便，又可保證計算精度。

[1] 劉金昊，吳建華. 西山供水工程事故停泵水力過渡過程計算及水錘防護[J]. 水電能源科學(xué)，2021，39(7)：113-116.

Liu Jinhao, Wu Jianhua. Calculation of hydraulic transition process and water hammer protection for accidental pump-stop in xishan water supply project [J]. Water Resources and Power, 2021, 39(7): 113-116. (in Chinese with English abstract)

[2] 黃偉，楊開林，郭新蕾，等. 抽水蓄能電站極端甩負荷工況球閥協(xié)同調(diào)節(jié)[J]. 清華大學(xué)學(xué)報（自然科學(xué)版），2019，59(8)：635-644.

Huang Wei, Yang Kailin, Guo Xinlei, et al. Coordinated regulation of ball valves in pumped storage stations for extreme conditions[J]. Journal of Tsinghua University (Science & Technology), 2019, 59(8): 635-644. (in Chinese with English abstract)

[3] 楊玉思，徐艷艷，羨巨智. 長距離高揚程多起伏輸水管道水錘防護的研究[J]. 給水排水，2009，45(4)：108-211.

Yang Yusi, Xu Yanyan, Xian Juzhi. Research on water hammer prevention in high lift, hilly and long distance water transmission pipeline[J]. Water & Wastewater Engineering, 2009, 45(4): 108-211. (in Chinese with English abstract)

[4] 方永，周宏. 管道系統(tǒng)中的水錘及其防護研究[J]. 石化技術(shù)，2016，23(2)：115.

Fang Yong, Zhou Hong. Water hammer and preventive measures in pipeline[J]. Petrochemical Industry Technology, 2016, 23(2): 115. (in Chinese with English abstract)

[5] Wylie E B, Streeter V L, Suo L. Fluid transients in systems[M]. Englewood Cliffs, NJ: Prentice Hall, 1993.

[6] Guinot V. Riemann solvers for water hammer simulations by Godunov method[J]. International Journal for Numerical Methods in Engineering, 2000, 49(7): 851-870.

[7] Zhao M, Ghidaoui M S. Godunov-type solutions for water hammer flows[J]. Journal of Hydraulic Engineering-ASCE, 2004, 130(4): 341-348.

[8] 趙越，周領(lǐng)，劉德有，等. 基于有限體積法Godunov格式的水錘計算模型[J]. 水利水電科技進展，2019，39(1)：76-81.

Zhao Yue, Zhou Ling, Liu Deyou, et al. Water hammer model based on finite volume method and Godunov-type scheme[J]. Advances in Science and Technology of Water Resources. 2019, 39(1): 76-81. (in Chinese with English abstract)

[9] Zhou L, Wang H, Liu D Y, et al. A second-order finite volume method for pipe flow with water column separation [J]. Journal of Hydro-environment Research, 2016, 17: 47-55.

[10] 向小華，吳曉玲，牛帥，等. 基于顯式有限體積法的一維河網(wǎng)模型[J]. 水利水電科技進展，2015，35(4)：6-9，43.

Xiang Xiaohua, Wu Xiaoling, Niu Shuai, et al. Construction of one-dimensional river network model based on explicit finite volume method[J]. Advances in Science and Technology of Water Resources, 2015, 35(4): 6-9, 43. (in Chinese with English abstract)

[11] Hu Y Y, Zhou L, Pan T W, et al. Godunov-type solutions for free surface transient flow in pipeline incorporating unsteady friction[J]. Journal of Water Supply: Research and Technology-Aqua, 2022, 71(4): 546-562.

[12] 鄭劼恒，蔣明，郭芮，等．順序輸送管道水力瞬變模擬的有限體積法[J]. 計算力學(xué)學(xué)報，2015，32(3)：418-422，428.

Zheng Jieheng, Jiang Ming, Guo Rui, et al. Finite volume method for hydraulic transient simulation of batching pipeline[J]. Chinese Journal of Computational Mechanics, 2015, 32(3): 418-422, 428. (in Chinese with English abstract)

[13] 楊開林. 電站與泵站中的水力瞬變及調(diào)節(jié)[M]. 北京：中國水利水電出版社，2000.

[14] Toro E F. Riemann solvers and numerical methods for fluid dynamics[M]. Berlin: Verlag Berlin Heidelberg, 2009: 115-125.

[15] 劉錦濤，章少輝，許迪，等. 灌溉輸配水系統(tǒng)明滿流的全隱式耦合模擬及驗證[J]. 農(nóng)業(yè)工程學(xué)報，2017，33(19)：124-130.

Liu Jintao, Zhang Shaohui, Xu Di, et al. Coupled simulation and validation with fully implicit time scheme for free-surface-pressurized water flow in pipe/channel[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2017, 33(19): 124-130. (in Chinese with English abstract)

[16] 王振華，馬習(xí)賀，李文昊，等．基于改進4-方程摩擦模型的輸水管道水錘壓力計算[J]. 農(nóng)業(yè)工程學(xué)報，2018，34(7)：114-120.

Wang Zhenhua, Ma Xihe, Li Wenhao, et al. Calculation of water hammer pressure of flow pipeline based on modified four-equation friction model[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2018, 34(7): 114-120. (in Chinese with English abstract)

[17] Chaudhry M H. 實用水力瞬變過程[M]. 北京：中國水利水電出版社，2015.

[18] Lyu J W, Zhang J, Wang T Y. Study on water hammer protection of the siphon breaking structure in the water supply system[J]. AQUA-Water Infrastructure, Ecosystems and Society, 2022, 71(3): 478-489.

[19] Abdel-Fatah M A, Amin A, Al Bazedi G A. Model and protected design of water piping system to minimize the water hammer effect[J]. Egyptian Journal of Chemistry, 2022, 65(1): 1-2.

[20] Lupa S I, Gagnon M, Muntean S, et al. The impact of water hammer on hydraulic power units[J]. Energies, 2022, 15(4): 1526.

[21] 富友，蔣勁，李燕輝，等. 改進雙流體模型計算有液柱分離的管路水錘壓力[J]. 農(nóng)業(yè)工程學(xué)報，2018，34(15)：58-65.

Fu You, Jiang Jing, Li Yanhui, et al. Calculation of pipe water hammer pressure with liquid column separation by improved two-fluid model[J]. Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2018, 34(15): 58-65. (in Chinese with English abstract)

[22] 周領(lǐng)，陸燕清．排水管道瞬變流的SWMM模擬能力研究[J]. 中國給水排水，2022，38(5)：108-115.

Zhou Ling, Lu Yanqing. SWMM simulation capability for transient fow in drainage pipe[J]. Water & Wastewater Engineering, 2022, 38(5): 108-115. (in Chinese with English abstract)

[23] 吳峰. 輸水管道中含氣水流瞬變特性試驗研究及數(shù)值模擬[D]. 西安：西安理工大學(xué)，2021.

Wu Feng. Experimental Research and Numerical Simulation on Transient Characteristics of Air-Bearing Water Flow in Water Pipeline[D]. Xi’an: Xi’an University of Technology, 2021. (in Chinese with English abstract)

[24] 丁澤. 有壓管道中閥門工作狀態(tài)檢測研究[D]. 哈爾濱：哈爾濱商業(yè)大學(xué)，2019.

Ding Ze. Research on Valve Working State Detection in Pressurized Pipeline[D]. Harbin: Harbin University of Commerce, 2019. (in Chinese with English abstract)

[25] 趙思茂. 輸水管道內(nèi)瞬變流特性的三維數(shù)值模擬研究[D]. 西安：西安理工大學(xué)，2020.

Zhao Simao. Three Dimensional Numerical Simulation of Transient Flow in Water Pipeline[D]. Xi’an: Xi’an University of Technology, 2020. (in Chinese with English abstract)

[26] Al-Khayat R H, Kadhim A A, Al-Baghdadi M A R S, et al. Flow parameters effect on water hammer stability in hydraulic system by using state-space method[J]. Open Engineering, 2022, 12(1): 215-226.

[27] Waqar M, Louati M, Ghidaoui M S. Time-reversal of water-hammer waves[J]. Journal of Hydraulic Research, 2022, 60(1): 25-45.

[28] Li Y J, Hu X D, Zhou F J, et al. A new comprehensive filtering model for pump shut-in water hammer pressure wave signals during hydraulic fracturing[J]. Journal of Petroleum Science and Engineering, 2022, 208: 109796.

[29] Kim S G, Lee K B, Kim K Y. Water hammer in the pump-rising pipeline system with an air chamber[J]. Journal of Hydrodynamics, 2014, 26(6): 960-964.

[30] Wan W Y, Li F Q. Sensitivity analysis of operational time differences for a pump-valve system on a water hammer response[J]. Journal of Pressure Vessel Technology, 2016, 138(1): 011303.

Construction of the second-order Godunov scheme model for hydraulic transients in pumping stations

Zhou Ling1,2, Hu Anni1, Wu Jinyuan1

(1.,,210098,; 2.,,210098,)

This study aims to implement a more efficient and stable numerical simulation of the hydraulic transient in a complex pumping station system. A finite volume method (FVM) Godunov scheme was established to simulate the simple pipeline and complex pumping station system. The FVM was then introduced to discretize the mathematical models, while the Riemann solver was selected to solve the discrete flux. The MUSCL-Hancock method was utilized to reconstruct the numerical data at the interface of control volumes. The higher numerical accuracy and stability were realized in the Godunov scheme, compared with the frequently-used method of characteristics. Meanwhile, the MINMOD slope limiter was used to avoid false oscillation. The boundary processing of the dual virtual unit was then presented for the second-order accuracy of both the computational region and the boundary, particularly for the simpler computation. The simulation of the improved model was in good agreement with the exact solution and the classical examples. The sensitivity analysis was also performed on the Courant and grid number. Furthermore, a more accurate, stable, and efficient performance was achieved in the second-order Godunov scheme, compared with the method of characteristics. More importantly, there was more outstanding attenuation with the decrease of the Courant number for a simple pipeline system. The computation time of the second-order Godunov scheme was 0.017 s at the same accuracy, compared with the method of characteristics (0.227 s). Consequently, a more stable and efficient performance was achieved in the second-order Godunov scheme. In the actual pumping system with the multiple-characteristics pipe structure, the second-order Godunov scheme required only a slight reduction in the Courant numbers, indicating the simple and convenient way for high numerical accuracy. Once the method of characteristics was used to calculate the hydraulic transition of the pumping station, the Courant number in the pipeline was less than 1 at the same length or wave velocity of the pipeline. By contrast, the Courant number was 0.72-0.76 in this case, indicating a very different simulation from the actual. Therefore, it is necessary to adjust the local pipeline length or wave velocity for the condition that the Courant number was 1. The tedious operation can lead to calculation errors, due to the change in pipeline characteristics. The accuracy can be improved but with less computational efficiency, if the wave velocity remained unchanged to increase the number of computational grids. In the method of characteristics, the number of grids can properly improve the accuracy of the calculation but with the doubled computation time, when the Courant number was less than 1. In the second-order Godunov scheme, there was little effect of grid number on the accuracy of the calculation but with the longer calculation time, whether the Courant number was equal to or less than 1. Therefore, a finer grid was preferred in the method of characteristics for the same accuracy requirements, when the Courant number was less than 1 in the transient process of the simulated pump system. Therefore, the second-order Godunov scheme can accurately simulate the process lines of rotational speed, discharge, and outlet pressure parameters during the hydraulic transient of the pump system. Anyway, the second-order Godunov scheme can be expected to effectively improve the efficiency, stability, and accuracy of hydraulic transient simulation of traditional pumping station systems.

pump station; models; finite volume method; second-order Godunov scheme; hydraulic transition process; method of characteristics; first-order Godunov scheme

10.11975/j.issn.1002-6819.2022.19.005

TV143.1

A

1002-6819(2022)-19-0042-09

周領(lǐng)，胡安妮，吳金遠. 泵站水力瞬變的二階Godunov格式模型構(gòu)建[J]. 農(nóng)業(yè)工程學(xué)報，2022，38(19)：42-50.doi：10.11975/j.issn.1002-6819.2022.19.005 http://www.tcsae.org

Zhou Ling, Hu Anni, Wu Jinyuan. Construction of the second-order Godunov scheme model for hydraulic transients in pumping stations[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(19): 42-50. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2022.19.005 http://www.tcsae.org

2022-05-12

2022-08-24

國家自然科學(xué)基金項目（51679066，51839008）；霍英東教育基金會青年教師基金（161068）

周領(lǐng)，博士，教授，研究方向為水電站、泵站、輸水工程水力學(xué)。Email：zlhhu@163.com