抽水蓄能電站水力瞬變的有限體積法建模模擬

周 領，吳金遠，王 豐，劉 靜，盧坤銘

(1.河海大學 水利水電學院，南京 210098；2.揚州大學 水利科學與工程學院，江蘇 揚州 225009；3.中國三峽建工(集團)有限公司，成都 610041)

抽水蓄能電站作為現代電力中不可或缺的能源調節結構[1-2]，有著調峰填谷、黑啟動等功能[3]。然而，為了滿足電力系統的動態服務要求，抽水蓄能電站有著一機多用、工況轉換迅速、啟停頻繁等特點[4]，常會導致整個機組產生較大的水力振動、共振等異常運行問題。因此，準確模擬各種水力瞬變對保證抽水蓄能電站安全穩定運行和精準控制極為重要。

目前，常采用特征線法(Method of characteristics,MOC)進行抽水蓄能電站水力瞬變的模擬計算。但是，抽水蓄能電站過流系統會涉及較多的短管、岔管、支管等情況，在特征線法求解中一般采取以下兩種方法解決，一是插值計算，但會導致計算效率與精度的降低；二是調整波速或者網格長度，或者直接忽略短管，簡化管網模型，相比于前者計算效率有所升高，但是會引入新的計算誤差。

近年來，有限體積法(Finite volume method,FVM)逐漸被運用于有壓管道水力瞬變計算。在保證系統質量與能量守恒前提下，FVM可以有效地處理非連續問題，避免虛假數值振蕩。Guinot[5]最早將有限體積法運用于水錘問題，得到了和特征線法相類似的格式。隨后，Zhao等[6]基于Godunov格式與Riemann求解器，得到了一階和二階的水錘求解格式。鄭劼恒等[7]在順序輸送管道的水力瞬變中，采用有限體積法的Godunov格式，并通過Riemann求解器以及MUSCL方法對通量進行計算，從而有效避免了特征線法求解中的虛假振蕩。Zhou等[8]將控制體等分，假定僅在控制體中心出現空穴，從而實現了Godunov格式對水柱分離-彌合水力現象的模擬。趙越等[9]研究了Godunov格式下庫朗數、對流項等參數的敏感性。

為解決MOC在處理抽蓄管網系統工序復雜、精度較低的問題，本文采用二階Godunov格式的FVM，并將虛擬邊界與機組控制方程相結合，實現了對某抽水蓄能電站的水力瞬變模擬，并將結果與MOC計算值、實驗值進行對比研究。

1 數學模型及其求解

1.1 水錘控制方程

對于管道內的非定常流，其連續方程和動量方程可寫成如下的微分方程形式[10]：

(1)

(2)

式(1)～(2)可改寫成矩陣形式：

(3)

式中：H為測壓管水頭，m；V為流速，m/s；g為重力加速度，m/s2；a為波速，m/s；D為管道直徑，m；f為管道恒定摩阻系數；S0為管道坡度；x為沿管軸線距離，m；t為時間，s。

若采用Riemann求解格式進行求解，且不考慮對流項，則可將式(3)改為經典水錘方程：

(4)

有限體積法是將計算區域離散為多個單元體，并對各單元體單獨積分求解，如圖1所示。

圖1 計算區域網格Fig.1 Grid of computational region

鑒于控制變量U在各時間和空間內均為連續分布，且在各控制體內分布平均，由此可得到相應的積分公式：

(5)

式中：i為第i個控制體；Fi+1/2為控制體右邊界處的通量；Fi-1/2為控制體左邊界處的通量；Δt為時間步長，s；Δx為空間步長，m；上標n表示t時刻；上標n+1表示t+Δt時刻。

1.2 水泵水輪機控制方程

1.2.1 全特性曲線

水泵水輪機的全特性曲線用以求解水泵水輪機在各瞬變時刻下各項特征參數的值。由于全特性曲線在各象限內均存在著開度線交叉、聚集、多值性等特點，因此在具體計算時需要對全特性曲線進行曲線轉換。Suter變換[11]是在水泵水輪機中最常用的變換方式，其具體形式如下：

(6)

(7)

(8)

式中：WH、WM分別為水頭特性函數和力矩特性函數,x為相對流量角,y為相對導葉開度,q=Q11/Q11r為相對單位流量,n=N11/N11r為相對單位轉速,h=H/Hr為相對水頭,m=M11/M11r為相對單位力矩,下標11表示單位值，下標r表示額定值。

但是常規的Suter變換無法表示0開度線下轉速、流量、力矩間的關系，且常規的Suter變換后的曲線在小開度情況下曲線遠稀疏于大開度情況下，導致計算結果在小開度情況下不夠精確[12]，而抽水蓄能電站在小開度情況下極易產生不穩定的情況。介于上述原因，本文采用改進的Suter變換[13]，不僅可以表示0開度線下各參數的關系，對于小開度下曲線的疏密問題也有所改善，具體形式如下：

(9)

(10)

(11)

式中c為常數，一般取1.0～1.5，在本文計算中，c取1.2。其余各參數含義同式(6)～(8)。

1.2.2 機組邊界條件

1)轉速平衡方程。在機組全甩荷工況下，轉速平衡方程如下[13-14]：

(12)

2)水頭平衡方程。設蝸殼前和尾水管后壓力鋼管分別為節點1和2，對這兩個節點分別用特征線方程，并帶入水輪機水頭計算方程，則可得到水頭平衡方程[13-14]：

h=[Cp1-Cm2-(Bp1+Bm2)Qrq+C2|q|q]/Hr

(13)

聯立式(9)、(10)、(12)和式(13)，即可求出各瞬變時刻機組的水頭、流量、轉速、力矩等參數。

1.3 二階Godunov求解格式

1.3.1 通量計算

由于各物理變量在各單位體內均是連續的，而在通量邊界上是間斷的，因此為求出Godunov格式下通量邊界處的值，可采用Riemann問題的求解方式[15]：

(14)

由此，可求出各單元體在各邊界處的通量值：

(15)

Step1數據重組。為避免重組數據時產生虛假振蕩，分別引入MINMOD、MUSCL和SUPERBEE 3種斜率限制器函數，并在后文中分析比較其異同，則可得到：

(16)

(17)

式中，Δi由斜率限制器函數計算得到，對于MINMOD函數：

(18)

對于MUSCL函數:

(19)

對于SUPERBEE函數:

(20)

(21)

(22)

Step2推進時間計算。

(23)

(24)

Step3Riemann問題求解。

(25)

(26)

將求解出的式(25)、(26)帶入式(15)，則可求出各單元體邊界處的通量。

1.3.2 時間積分

在得到二階精度的通量計算值后要得到n時刻到n+1時刻的解，需要對式(5)進行積分求解，采用二階顯式Runge-Kutta，以得到二階計算精度，計算過程如下：

(27)

(28)

(29)

若采用二階求解格式的通量計算值，則式(27)～(29)所得到的格式在空間與時間上均為二階精度。

1.3.3 虛擬邊界

根據上述求解方式可知，在對任一單元體i進行二階求解時，需要單元體i左右各兩個單元體的物理變量值，因此對于邊界處的單元體需要進行特定處理如圖2所示。在本文中，采取在邊界兩邊分別添加-1，0和N+1，N+2的虛擬單元的處理方法。

圖2 計算區域與計算網格Fig.2 Computational region and grid

對于添加的虛擬單元滿足以下條件：

U-1=U0=U1/2

(30)

UN+1=UN+2=UN+1/2

(31)

且各虛擬單元的物理變量分別滿足邊界處的Riemann不變量方程。

1)水庫處。

對于上游水庫，方程如下：

(32)

對于下游水庫，方程如下：

(33)

2)機組處。

根據機組控制方程，只需求出蝸殼處與尾水管處的虛擬單元的物理變量值，便可求出機組在各瞬變時刻下的物理變量值。因此，結合Riemann不變量方程，可得：

(34)

(35)

(36)

(37)

2 計算分析

2.1 簡單管道系統算例

設置一上游為水庫，下游為閥門的簡單管道，管道長500 m，計算區域共10個，波速為1 000 m/s，上游水庫為10 m，初始流速為0.1m/s，重力加速度為9.8 m/s2，總的計算時間取10 s，下游閥門設置為瞬時關閉，管道無模阻，則結果中的所有壓力衰減均是由于數值耗散引起的。

分別用MOC和二階Godunov格式的FVM對上述簡單系統進行水錘求解，主要研究內容包括：1)比較分析MINMOD、MUSCL與SUPERBEE 3種斜率限制器函數對水錘計算的影響；2)研究庫朗數Cr(1.0、0.5、0.2和0.1)對兩種求解格式計算結果的影響；3)比較分析FVM與MOC水錘計算的精度和效率。

如圖3(a)所示，當Cr=1.0時，MINMOD、MUSCL和SUPERBEE 3種斜率限制器函數計算完全相同；而如圖3(b)所示，當Cr=0.5時，3種斜率限制器函數均會導致一定程度的數值耗散。然而，MUSCL與SUPERBEE會產生一定的虛假的數值振蕩，而MINMOD計算更加穩定。因此，本文選擇MINMOD斜率限制器。

圖3 不同斜率限制器在不同庫朗數條件下的比對圖Fig.3 Comparison of different slope limiters with different Courant numbers

如圖4、5所示，當Cr=1.0時，兩種方法計算結果與精確解完全相同，證明了二階Godunov格式的FVM水錘計算的準確性。但是當Cr<1.0時，MOC和FVM均會出現不同程度的數值耗散，且Cr越小，耗散程度越嚴重。如圖5所示，在相同庫朗數條件下，二階Godunov格式的FVM數值耗散遠小于MOC。說明在Cr<1.0的條件下，二階Godunov格式的FVM可以有效抑制數值耗散，計算更穩定，結果更精確。

圖4 MOC計算結果Fig.4 MOC calculation results

圖5 FVM計算結果Fig.5 FVM calculation results

如圖6所示，當Cr=0.5時，網格數為512個的MOC格式與網格數為64個的FVM格式，計算精度基本一致。由表1可知，網格數為64個的FVM格式CPU計算時間為0.310 s，而網格數為512個的MOC格式CPU計算時間為1.080 s，約為FVM計算時間的3倍。說明在相同計算精度的情況下，FVM的計算效率要遠遠大于MOC。

圖6 FVM與不同網格數的MOC對比圖Fig.6 Comparison of FVM and MOC with different grid numbers

表1 各模型計算時間Tab.1 Computational time of models s

2.2 某抽水蓄能電站水力瞬變計算分析

2.2.1 電站參數

電站設有兩臺150 MW的可逆式機組，上游引水系統采用“一管雙機”、下游尾水系統采用“一管一機”的布置方式，布置如圖7所示，相關設計參數見表2、3。

圖7 某抽水蓄能電站布置圖Fig.7 Layout of a pumped storage power station

表2 管道參數Tab.2 Parameters of water pipe system

表3 電站機組參數Tab.3 Unit parameters of power station

2.2.2 計算結果

在機組甩荷過程中，機組轉速的波動是造成水輪機無法穩定運行最直接原因。因此將MOC和二階Godunov格式的FVM計算出的轉速與甩荷實驗值進行對比。但是由于MOC在Cr<1.0時會產生較嚴重的數值衰減，因此需要對管道內波速進行調整，來滿足Cr=1.0的條件。調整后的波速與管道分段數見表4。

表4 管道系統各部分管段波速及其分段數Tab.4 Section number and wave speeds of pipe system

本文選取了3組實驗計算工況，已知3組實驗工況的輸入功率、初始流量與初始轉速均為額定值，其他相關參數見表5。

表5 機組100%甩荷計算工況參數Tab.5 Calculation parameters of 100% load rejection

兩種方法在整個計算時間段內計算的轉速以及蝸殼處水頭的變化情況見表6以及如圖8～10所示。

表6 機組100%甩荷工況結果Tab.6 Calculation results of 100% load rejection (r·min-1)

由圖8(a)、圖9(a)、圖10(a)所示，MOC與FVM在對該抽水蓄能電站模型進行甩荷計算時，計算結果與實驗值基本吻合，表明二階Godunov格式的FVM在處理較復雜管網模型時的適用性和準確性。由表6中數據與圖8(a)、圖9(a)、圖10(a)所示，MOC與FVM計算的最大轉速均略小于實驗值，但FVM計算結果更加接近于實驗值，這是因為MOC在計算過程中，需要對管道波速進行調整來滿足庫朗數條件，由此帶來了一定的誤差，且調整波速間接地增加了計算時間，導致其計算效率的降低。而FVM僅適當降低庫朗數條件，無需進行波速調整，簡化了計算過程，具有較高的計算效率與精度。

圖8 甩荷實驗1Fig.8 Load rejection experiment 1

圖9 甩荷實驗2Fig.9 Load rejection experiment 2

圖10 甩荷實驗3Fig.10 Load rejection experiment 3

由圖8(b)、圖9(b)、圖10(b)所示，MOC與FVM地瞬變壓力計算結果基本一致，但是很明顯兩者在轉速變化較快時，均會產生一定的數值波動，但FVM計算結果相對更穩定。

本文進行該抽水蓄能電站模擬時進行了一定的簡化處理，如忽略了平水建筑物的影響，從而導致機組轉速達到最大值后的計算結果與實測數據存在一定的差異。不過，這對于本文結論的影響不大。

3 結 論

1)本文所建的二階Godunov格式的FVM可準確、高效地實現對某抽水蓄能電站的水力瞬變計算模擬，且比MOC計算結果更接近實測值。

2)在二階Godunov格式中，MINMOD、MUSCL與SUPERBEE 3種不同斜率限制器在Cr=1.0時均能得到精確的水錘結果；但在Cr<1.0時，3種斜率限制器函數均會導致一定的數值耗散，但MUSCL與SUPERBEE會產生虛假的數值振蕩，而MINMOD計算更加穩定。

3)在Cr=1.0時，二階Godunov格式的FVM和MOC兩者計算結果一致。但是在Cr<1.0時，兩種方法均會出現一定程度的數值耗散，但二階Godunov的FVM可以有效地抑制數值耗散，計算更加穩定。

4)為了得到相同精度的計算結果，MOC需要更加密集的網格數，計算耗時長，而二階Godunov格式的FVM在只需較稀疏的網格即可實現同樣的計算精度，計算效率更加高效。

5)在處理抽水蓄能電站的水力瞬變問題時，MOC為了滿足庫朗數條件，需要進行波速調整，或者直接忽略短管，因而產生了計算誤差；二階Godunov的FVM僅適當降低庫朗數條件，無需進行波速調整，簡化了計算過程，具有較高的計算效率和精度。