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

周 領(lǐng)，吳金遠(yuǎn)，王 豐，劉 靜，盧坤銘

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

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

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

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

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

1 數(shù)學(xué)模型及其求解

1.1 水錘控制方程

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

(1)

(2)

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

(3)

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

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

(4)

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

圖1 計算區(qū)域網(wǎng)格Fig.1 Grid of computational region

鑒于控制變量U在各時間和空間內(nèi)均為連續(xù)分布，且在各控制體內(nèi)分布平均，由此可得到相應(yīng)的積分公式：

(5)

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

1.2 水泵水輪機控制方程

1.2.1 全特性曲線

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

(6)

(7)

(8)

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

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

(9)

(10)

(11)

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

1.2.2 機組邊界條件

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

(12)

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

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

(13)

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

1.3 二階Godunov求解格式

1.3.1 通量計算

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

(14)

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

(15)

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

(16)

(17)

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

(18)

對于MUSCL函數(shù):

(19)

對于SUPERBEE函數(shù):

(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 虛擬邊界

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

圖2 計算區(qū)域與計算網(wǎng)格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)機組處。

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

(34)

(35)

(36)

(37)

2 計算分析

2.1 簡單管道系統(tǒng)算例

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

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

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

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

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

圖4 MOC計算結(jié)果Fig.4 MOC calculation results

圖5 FVM計算結(jié)果Fig.5 FVM calculation results

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

圖6 FVM與不同網(wǎng)格數(shù)的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 電站參數(shù)

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

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

表2 管道參數(shù)Tab.2 Parameters of water pipe system

表3 電站機組參數(shù)Tab.3 Unit parameters of power station

2.2.2 計算結(jié)果

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

表4 管道系統(tǒng)各部分管段波速及其分段數(shù)Tab.4 Section number and wave speeds of pipe system

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

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

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

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

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

圖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地瞬變壓力計算結(jié)果基本一致，但是很明顯兩者在轉(zhuǎn)速變化較快時，均會產(chǎn)生一定的數(shù)值波動，但FVM計算結(jié)果相對更穩(wěn)定。

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

3 結(jié) 論

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

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

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

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

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