抽水蓄能电站水力瞬变的有限体积法建模模拟

周领¹,吴金远¹,王丰²,刘静¹,卢坤铭³

(1.河海大学 水利水电学院,南京 210098;2.扬州大学 水利科学与工程学院,江苏 扬州 225009; 3.中国三峡建工(集团)有限公司,成都 610041)



摘要:

针对抽水蓄能电站管道内水力瞬变问题,采用二阶有限体积法（Finite Volume Method, FVM）Godunov格式进行数值建模和模拟。首先根据有限体积法将数学模型的控制方程进行离散,采用Riemann求解器进行通量计算。为了得到二阶精度的Godunov格式,避免数据重构时产生虚假振荡,引入并分析比较了3种不同的斜率限制器。机组全特性曲线采用改进的Suter变换,在边界处加入虚拟边界,并与机组边界控制方程相结合,实现了管道内部控制体与边界处计算的统一性。将本研究所建数学模型计算结果与传统特征线法（Method of Characteristics, MOC）计算结果、实测值进行对比,并分析了相关参数的敏感性。结果表明:在库朗数等于1时,二阶FVM与MOC计算结果相同,当库朗数小于1时,FVM计算结果比MOC更准确、更稳定;针对抽蓄机组甩负荷过程,二阶FVM的机组转速计算值与实测值基本吻合,且比MOC计算结果更接近实测值;在抽水蓄能电站水力瞬变计算中,MOC需要对不同特性的管道内水锤波速进行调整,不仅增加计算复杂性,也会产生计算误差；而FVM只需降低库朗数,计算简单方便,且精度更高。因此,本研究提供了一种稳定、精确的水力瞬变计算方法。

关键词:  抽水蓄能水电站  有限体积法  Godunov格式  特征线法  水力瞬变

DOI：10.11918/202101057

分类号:TV143.1

文献标识码:A

基金项目:国家自然科学基金(6,8)；霍英东教育基金会青年教师基金(161068)



Numerical simulation of hydraulic transients in pumped storage power station with finite volume method

ZHOU Ling¹,WU Jinyuan¹,WANG Feng²,LIU Jing¹,LU Kunming³

(1.College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098,China; 2.College of Hydraulic Science and Engineering, Yangzhou University, Yangzhou 225009, Jiangsu, China; 3.China Three Gorges Construction Engineering Corporation, Chengdu 610041, China)

Abstract:

The second-order finite volume method (FVM) Godunov scheme was adopted to investigate the hydraulic transients in pumped storage power station. Firstly, the governing equations of mathematical model were discretized according to FVM, and the flux was calculated by Riemann solver. Three slope limiters were introduced and compared to avoid spurious oscillations during data reconstruction. The whole characteristic curves of the unit were transformed by using improved Suter approach. The virtual-boundary approach combined with the governing equations of the unit was presented to achieve a unified computation scheme for all the control volumes at internal domain and boundaries. The results calculated by the proposed scheme were compared with those of method of characteristics (MOC) scheme and measured data, and parameter sensitivity was discussed. Results show that when Courant number was equal to 1, the second-order FVM and MOC had the same accuracy; when Courant number was less than 1, the second-order FVM was more accurate and stable than MOC. During load rejection process of pumped storage station, the unit speed calculated by second-order FVM was basically consistent with the measured value and more accurate than that predicted by MOC. For the calculation of hydraulic transients in pumped storage power station, the values of wave speed in pipe sections with different properties should be adjusted to realize MOC simulations, resulting in computation complexity and numerical errors. While FVM is simpler in computation with higher accuracy, which only requires reducing the Courant number. Therefore, the proposed numerical approach is stable and accurate for hydraulic transients.

Key words:  pumped storage power station  finite volume method (FVM)  Godunov scheme  method of characteristics (MOC)  hydraulic transients


周领, 吴金远, 王丰, 刘静, 卢坤铭. 抽水蓄能电站水力瞬变的有限体积法建模模拟[J]. 哈尔滨工业大学学报, 2022, 54(6): 79-86. DOI: 10.11918/202101057.
ZHOU Ling, WU Jinyuan, WANG Feng, LIU Jing, LU Kunming. Numerical simulation of hydraulic transients in pumped storage power station with finite volume method[J]. Journal of Harbin Institute of Technology, 2022, 54(6): 79-86. DOI: 10.11918/202101057.
基金项目 国家自然科学基金(51679066，51839008)；霍英东教育基金会青年教师基金(161068) 作者简介 周领(1985-)，男，教授，博士生导师 通信作者 周领，zlhhu@163.com 文章历史 收稿日期: 2021-01-15



Abstract            Full text            Figures/Tables            PDF


抽水蓄能电站水力瞬变的有限体积法建模模拟
周领¹, 吴金远¹, 王丰², 刘静¹, 卢坤铭³    
1. 河海大学 水利水电学院，南京 210098;
2. 扬州大学 水利科学与工程学院，江苏 扬州 225009;
3. 中国三峡建工(集团)有限公司，成都 610041

收稿日期: 2021-01-15
基金项目: 国家自然科学基金(51679066，51839008)；霍英东教育基金会青年教师基金(161068)
作者简介: 周领(1985-)，男，教授，博士生导师
通信作者: 周领，zlhhu@163.com


摘要: 针对抽水蓄能电站管道内水力瞬变问题，采用二阶有限体积法(Finite Volume Method, FVM)Godunov格式进行数值建模和模拟。首先根据有限体积法将数学模型的控制方程进行离散，采用Riemann求解器进行通量计算。为了得到二阶精度的Godunov格式，避免数据重构时产生虚假振荡，引入并分析比较了3种不同的斜率限制器。机组全特性曲线采用改进的Suter变换，在边界处加入虚拟边界，并与机组边界控制方程相结合，实现了管道内部控制体与边界处计算的统一性。将本研究所建数学模型计算结果与传统特征线法(Method of Characteristics, MOC)计算结果、实测值进行对比，并分析了相关参数的敏感性。结果表明: 在库朗数等于1时，二阶FVM与MOC计算结果相同, 当库朗数小于1时，FVM计算结果比MOC更准确、更稳定; 针对抽蓄机组甩负荷过程，二阶FVM的机组转速计算值与实测值基本吻合，且比MOC计算结果更接近实测值; 在抽水蓄能电站水力瞬变计算中，MOC需要对不同特性的管道内水锤波速进行调整，不仅增加计算复杂性，也会产生计算误差；而FVM只需降低库朗数，计算简单方便，且精度更高。因此，本研究提供了一种稳定、精确的水力瞬变计算方法。
关键词: 抽水蓄能水电站    有限体积法    Godunov格式    特征线法    水力瞬变    
Numerical simulation of hydraulic transients in pumped storage power station with finite volume method
ZHOU Ling¹, WU Jinyuan¹, WANG Feng², LIU Jing¹, LU Kunming³    
1. College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China;
2. College of Hydraulic Science and Engineering, Yangzhou University, Yangzhou 225009, Jiangsu, China;
3. China Three Gorges Construction Engineering Corporation, Chengdu 610041, China



Abstract: The second-order finite volume method (FVM) Godunov scheme was adopted to investigate the hydraulic transients in pumped storage power station. Firstly, the governing equations of mathematical model were discretized according to FVM, and the flux was calculated by Riemann solver. Three slope limiters were introduced and compared to avoid spurious oscillations during data reconstruction. The whole characteristic curves of the unit were transformed by using improved Suter approach. The virtual-boundary approach combined with the governing equations of the unit was presented to achieve a unified computation scheme for all the control volumes at internal domain and boundaries. The results calculated by the proposed scheme were compared with those of method of characteristics (MOC) scheme and measured data, and parameter sensitivity was discussed. Results show that when Courant number was equal to 1, the second-order FVM and MOC had the same accuracy; when Courant number was less than 1, the second-order FVM was more accurate and stable than MOC. During load rejection process of pumped storage station, the unit speed calculated by second-order FVM was basically consistent with the measured value and more accurate than that predicted by MOC. For the calculation of hydraulic transients in pumped storage power station, the values of wave speed in pipe sections with different properties should be adjusted to realize MOC simulations, resulting in computation complexity and numerical errors. While FVM is simpler in computation with higher accuracy, which only requires reducing the Courant number. Therefore, the proposed numerical approach is stable and accurate for hydraulic transients.
Keywords: pumped storage power station    finite volume method (FVM)    Godunov scheme    method of characteristics (MOC)    hydraulic transients    
抽水蓄能电站作为现代电力中不可或缺的能源调节结构^[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]：

$g H_{x}+V V_{x}+V_{t}+\left(\frac{f V|V|}{2 D}-g S_{0}\right)=0$ (1)

$V H_{x}+H_{t}+\frac{a^{2}}{g} V_{x}=0$ (2)

式(1)~(2)可改写成矩阵形式：

$\frac{\partial \boldsymbol{U}}{\partial t}+\boldsymbol{A} \frac{\partial \boldsymbol{U}}{\partial x}=\boldsymbol{S}$ (3)

其中，$\mathit{\boldsymbol{U}} = \left( {\begin{array}{*{20}{l}}H\\V\end{array}} \right), \mathit{\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}}V&{{a^2}/g}\\g&V\end{array}} \right), \mathit{\boldsymbol{S}} = \left( {\begin{array}{*{20}{c}}0\\{g{S_0} - \frac{{fV|V|}}{{2D}}}\end{array}} \right)$

式中：H为测压管水头，m；V为流速，m/s；g为重力加速度，m/s²；a为波速，m/s；D为管道直径，m；f为管道恒定摩阻系数；S₀为管道坡度；x为沿管轴线距离，m；t为时间，s。

若采用Riemann求解格式进行求解，且不考虑对流项，则可将式(3)改为经典水锤方程：

$\frac{\partial \boldsymbol{U}}{\partial t}+\frac{\partial \boldsymbol{F}}{\partial x}=\boldsymbol{S}$ (4)

式中：$\mathit{\boldsymbol{F}} = \mathit{\boldsymbol{\overline A}} \mathit{\boldsymbol{U}}, \mathit{\boldsymbol{\overline A}} = \left( {\begin{array}{*{20}{c}}0&{{a^2}/g}\\g&0\end{array}} \right)$。

有限体积法是将计算区域离散为多个单元体，并对各单元体单独积分求解，如图 1所示。

Fig. 1
图 1 计算区域网格 Fig. 1 Grid of computational region


鉴于控制变量U在各时间和空间内均为连续分布，且在各控制体内分布平均，由此可得到相应的积分公式：

$\boldsymbol{U}_{i}^{n+1}=\boldsymbol{U}_{i}^{n}-\frac{\Delta t}{\Delta x}\left(\boldsymbol{F}_{i+\frac{1}{2}}-\boldsymbol{F}_{i-\frac{1}{2}}\right)+\frac{\Delta t}{\Delta x} \int_{i-\frac{1}{2}}^{i+\frac{1}{2}} \boldsymbol{S} \mathrm{d} x$ (5)

式中：i为第i个控制体；F_i+1/2为控制体右边界处的通量；F_i-1/2为控制体左边界处的通量；Δt为时间步长，s；Δx为空间步长，m；上标n表示t时刻；上标n+1表示t+Δt时刻。

1.2 水泵水轮机控制方程 1.2.1 全特性曲线水泵水轮机的全特性曲线用以求解水泵水轮机在各瞬变时刻下各项特征参数的值。由于全特性曲线在各象限内均存在着开度线交叉、聚集、多值性等特点，因此在具体计算时需要对全特性曲线进行曲线转换。Suter变换^[11]是在水泵水轮机中最常用的变换方式，其具体形式如下：

$W H(x, y)=\frac{1}{\left(\frac{Q_{11}}{Q_{11 r}}\right)^{2}+\left(\frac{N_{11}}{N_{11 r}}\right)^{2}}=\frac{h}{q^{2}+n^{2}}$ (6)

$W M(x, y)=\frac{\frac{M_{11}}{M_{11 r}}}{\left(\frac{Q_{11}}{Q_{11 r}}\right)^{2}+\left(\frac{N_{11}}{N_{11 r}}\right)^{2}}=\frac{m}{q^{2}+n^{2}}$ (7)

$\left\{\begin{array}{cc}x=\arctan (q / n), & n \geqslant 0 \\x=\arctan (q / n)+{\rm{ \mathsf{ π} }}, & n<0\end{array}\right.$ (8)

式中：WH、WM分别为水头特性函数和力矩特性函数, x为相对流量角, y为相对导叶开度, q=Q₁₁/Q_11r为相对单位流量, n=N₁₁/N_11r为相对单位转速, h=H/H_r为相对水头, m=M₁₁/M_11r为相对单位力矩, 下标11表示单位值，下标r表示额定值。

但是常规的Suter变换无法表示0开度线下转速、流量、力矩间的关系，且常规的Suter变换后的曲线在小开度情况下曲线远稀疏于大开度情况下，导致计算结果在小开度情况下不够精确^[12]，而抽水蓄能电站在小开度情况下极易产生不稳定的情况。介于上述原因，本文采用改进的Suter变换^[13]，不仅可以表示0开度线下各参数的关系，对于小开度下曲线的疏密问题也有所改善，具体形式如下：

$\begin{aligned}W H(x, y)=& \frac{1}{\left(\frac{Q_{11}}{Q_{11 r}}+c\right)^{2}+\left(\frac{N_{11}}{N_{11 r}}\right)^{2}}=\\& \frac{h}{(q+c \sqrt{h})^{2}+n^{2}}\end{aligned}$ (9)

$W M(x, y)=\frac{M_{11}}{M_{11 r}}=\frac{m}{h}$ (10)

$\begin{cases}x=\arctan [(q+c \sqrt{h}) / n], & n \geqslant 0 ; \\ x=\arctan [(q+c \sqrt{h}) / n]+{\rm{ \mathsf{ π} }}, & n<0 。\end{cases}$ (11)

式中c为常数，一般取1.0~1.5，在本文计算中，c取1.2。其余各参数含义同式(6)~(8)。

1.2.2 机组边界条件1) 转速平衡方程。在机组全甩荷工况下，转速平衡方程如下^[13-14]：

$n=n_{0}+\frac{\Delta t}{2 T_{\mathrm{a}}}\left(m+m_{0}\right)$ (12)

式中：T_a为机组惯性时间常数，${T_{\rm{a}}} = \frac{{G{D^2}N_r^2}}{{365{P_r}}}$；GD²为机组转动惯量，t · m²；N_r为机组额定转速，r/min；P_r为机组额定功率，kW；下标0表示上一时刻的值，其他符号含义同上。

2) 水头平衡方程。设蜗壳前和尾水管后压力钢管分别为节点1和2，对这两个节点分别用特征线方程，并带入水轮机水头计算方程，则可得到水头平衡方程^[13-14]：

$h=\left[C_{\mathrm{p} 1}-C_{\mathrm{m} 2}-\left(B_{\mathrm{p} 1}+B_{\mathrm{m} 2}\right) Q_{r} q+C_{2}|q| q\right] / H_{r}$ (13)

式中：系数${C_2} = Q_r^2\left( {\frac{1}{{2gA_1^2}} - \frac{1}{{2gA_2^2}}} \right)$；A为压力钢管面积，m²；其他符号含义同上。

联立式(9)、(10)、(12)和式(13)，即可求出各瞬变时刻机组的水头、流量、转速、力矩等参数。

1.3 二阶Godunov求解格式 1.3.1 通量计算由于各物理变量在各单位体内均是连续的，而在通量边界上是间断的，因此为求出Godunov格式下通量边界处的值，可采用Riemann问题的求解方式^[15]：

$\boldsymbol{U}_{x}^{n}= \begin{cases}\boldsymbol{U}_{\mathrm{L}}^{n}, & x<x_{i+1 / 2} \\ U_{\mathrm{R}}^{n}, & x>x_{i+1 / 2}\end{cases}$ (14)

式中：U_Lⁿ为单元体i在边界i+1/2左侧界面的平均值，U_Rⁿ为单元体i在边界i+1/2右侧界面的平均值。

由此，可求出各单元体在各边界处的通量值：

$\begin{aligned}\boldsymbol{F}_{i+\frac{1}{2}}=& \boldsymbol{A}_{i+\frac{1}{2}} \boldsymbol{U}_{i+\frac{1}{2}}(t)=\\& \frac{1}{2} \boldsymbol{A}_{i+\frac{1}{2}}\left\{\left(\begin{array}{cc}1 & \frac{a}{g} \\\frac{g}{a} & 1\end{array}\right) \boldsymbol{U}_{\mathrm{L}}^{n}+\left(\begin{array}{cc}1 & \frac{-a}{g} \\\frac{-g}{a} & 1\end{array}\right) \boldsymbol{U}_{\mathrm{R}}^{n}\right\}\end{aligned}$ (15)

对于二阶精度的Godunov通量计算格式，需要得到二阶精度的U_Lⁿ和U_Rⁿ，因此需要对原数据进行线性重构。

Step 1??数据重组。为避免重组数据时产生虚假振荡，分别引入MINMOD、MUSCL和SUPERBEE 3种斜率限制器函数，并在后文中分析比较其异同，则可得到：

$\boldsymbol{U}_{i}^{\mathrm{L}}=\boldsymbol{U}_{i}^{n}-\frac{\Delta x}{2} \Delta i$ (16)

$\boldsymbol{U}_{i}^{\mathrm{R}}=\boldsymbol{U}_{i}^{n}+\frac{\Delta x}{2} \Delta i$ (17)

式中，Δi由斜率限制器函数计算得到，对于MINMOD函数：

$\varDelta i= \begin{cases}\sigma_{i}^{n}, & \left|\sigma_{i}^{n}\right|<\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ \sigma_{i-1}^{n}, & \left|\sigma_{i}^{n}\right|>\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ 0, & \sigma_{i}^{n} \sigma_{i-1}^{n} \leqslant 0\end{cases}$ (18)

对于MUSCL函数:

$\varDelta i= \begin{cases}2 \sigma_{i}^{n}, & \left|\sigma_{i}^{n}\right| \leqslant \frac{1}{3}\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ \frac{\sigma_{i}^{n}+\sigma_{i-1}^{n}}{2}, & \frac{1}{3}\left|\sigma_{i-1}^{n}\right|<\left|\sigma_{i}^{n}\right|<3\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ 2 \sigma_{i-1}^{n}, & \left|\sigma_{i}^{n}\right| \geqslant 3\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ 0, & \sigma_{i}^{n} \sigma_{i-1}^{n} \leqslant 0\end{cases}$ (19)

对于SUPERBEE函数:

$\varDelta i= \begin{cases}2 \sigma_{i}^{n}, & 2\left|\sigma_{i}^{n}\right|<\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ \sigma_{i-1}^{n}, & \left|\sigma_{i-1}^{n}\right| \leqslant 2\left|\sigma_{i}^{n}\right| \leqslant 2\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ \sigma_{i}^{n}, & \left|\sigma_{i-1}^{n}\right|<\left|\sigma_{i}^{n}\right|<2\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ 2 \sigma_{i-1}^{n}, & \left|\sigma_{i}^{n}\right| \geqslant 2\left|\sigma_{i-1}^{n}\right| \text { 且 } \sigma_{i}^{n} \sigma_{i-1}^{n}>0 \\ 0, & \sigma_{i}^{n} \sigma_{i-1}^{n} \leqslant 0\end{cases}$ (20)

$\sigma_{i}^{n}=\left(\boldsymbol{U}_{i+1}^{n}-\boldsymbol{U}_{i}^{n}\right) / \Delta x$ (21)

$\sigma_{i-1}^{n}=\left(\boldsymbol{U}_{i}^{n}-\boldsymbol{U}_{i-1}^{n}\right) / \Delta x$ (22)

Step 2??推进时间计算。

$\overline{\boldsymbol{U}}{}_{i}^{\mathrm{R}}=\boldsymbol{U}_{i}^{\mathrm{R}}+\frac{\Delta t}{2 \Delta x}\left(\overline{\boldsymbol{A}} \boldsymbol{U}_{i}^{\mathrm{L}}-\overline{\boldsymbol{A}} \boldsymbol{U}_{i}^{\mathrm{R}}\right)$ (23)

$\overline{\boldsymbol{U}}{}_{i}^{\mathrm{L}}=\boldsymbol{U}_{i}^{\mathrm{L}}+\frac{\Delta t}{2 \Delta x}\left(\overline{\boldsymbol{A}} \boldsymbol{U}_{i}^{\mathrm{L}}-\overline{\boldsymbol{A}} \boldsymbol{U}_{i}^{\mathrm{R}}\right)$ (24)

Step 3??Riemann问题求解。

$\boldsymbol{U}_{\mathrm{L}}^{n}=\overline{\boldsymbol{U}}{}_{i}^{\mathrm{R}}$ (25)

$\boldsymbol{U}_{\mathrm{R}}^{n}=\overline{\boldsymbol{U}}{}_{i+1}^{\mathrm{L}}$ (26)

将求解出的式(25)、(26)带入式(15)，则可求出各单元体边界处的通量。

1.3.2 时间积分在得到二阶精度的通量计算值后要得到n时刻到n+1时刻的解，需要对式(5)进行积分求解，采用二阶显式Runge-Kutta，以得到二阶计算精度，计算过程如下：

$\boldsymbol{U}_{i}^{n+1}=\overline{\boldsymbol{U}}{}_{i}^{n+1}+\Delta t \boldsymbol{S}\left(\overline{\overline{\boldsymbol{U}}}{}_{i}^{n+1}\right)$ (27)

$\overline{\overline{\boldsymbol{U}}}{}_{i}^{n+1}=\overline{\boldsymbol{U}}{}_{i}^{n+1}+\frac{\Delta t}{2} \boldsymbol{S}\left(\overline{\boldsymbol{U}}{}_{i}^{n+1}\right)$ (28)

$\overline{\boldsymbol{U}}{}_{i}^{n+1}=\boldsymbol{U}_{i}^{n}-\frac{\Delta t}{\Delta x}\left(\boldsymbol{F}_{i+\frac{1}{2}}-\boldsymbol{F}_{i-\frac{1}{2}}\right)$ (29)

若采用二阶求解格式的通量计算值，则式(27)~(29)所得到的格式在空间与时间上均为二阶精度。

1.3.3 虚拟边界根据上述求解方式可知，在对任一单元体i进行二阶求解时，需要单元体i左右各两个单元体的物理变量值，因此对于边界处的单元体需要进行特定处理如图 2所示。在本文中，采取在边界两边分别添加-1，0和N+1，N+2的虚拟单元的处理方法。

Fig. 2
图 2 计算区域与计算网格 Fig. 2 Computational region and grid


对于添加的虚拟单元满足以下条件：

$\boldsymbol{U}_{-1}=\boldsymbol{U}_{0}=\boldsymbol{U}_{1 / 2}$ (30)

$\boldsymbol{U}_{N+1}=\boldsymbol{U}_{N+2}=\boldsymbol{U}_{N+1 / 2}$ (31)

且各虚拟单元的物理变量分别满足边界处的Riemann不变量方程。

1) 水库处。

对于上游水库，方程如下：

$V_{1 / 2}-\frac{g}{a} H_{1 / 2}=V_{1}^{n}-\frac{g}{a} H_{1}^{n}$ (32)

对于下游水库，方程如下：

$H_{N}^{n}+\frac{a}{g} V_{N}^{n}=H_{N+1 / 2}+\frac{a}{g} V_{N+1 / 2}$ (33)

式中：V₁ⁿ、H₁ⁿ分别为与上游水库相连的压力管道的第一个单元体的流速和水头, V_Nⁿ、H_Nⁿ分别为与下游水库相连的压力管道的最后一个单元体的流速和水头, H_1/2、H_N+1/2分别为上、下游水库的水位高。因此在任意时刻可直接求出水库处虚拟单元的水头和流速值。

2) 机组处。

根据机组控制方程，只需求出蜗壳处与尾水管处的虚拟单元的物理变量值，便可求出机组在各瞬变时刻下的物理变量值。因此，结合Riemann不变量方程，可得：

$C_{\mathrm{p} 1}=H_{N}^{n}+\frac{a}{g} V_{N}^{n}$ (34)

$B_{\mathrm{p} 1}=\frac{a}{g A_{1}}$ (35)

$C_{\mathrm{m} 2}=H_{1}^{n}+\frac{a}{g} V_{1}^{n}$ (36)

$B_{\mathrm{m} 2}=\frac{a}{g A_{2}}$ (37)

式中：V_Nⁿ、H_Nⁿ分别为与蜗壳相连的压力管道的最后一个单元体的流速和水头, V₁ⁿ、H₁ⁿ分别为与尾水管相连的压力管道的第1个单元体的流速和水头。将得到的式(34)~(37)带入式(13)即可求解虚拟边界条件下的水头平衡方程。

2 计算分析 2.1 简单管道系统算例设置一上游为水库，下游为阀门的简单管道，管道长500 m，计算区域共10个，波速为1 000 m/s，上游水库为10 m，初始流速为0.1m/s，重力加速度为9.8 m/s²，总的计算时间取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斜率限制器。

Fig. 3
图 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可以有效抑制数值耗散，计算更稳定，结果更精确。

Fig. 4
图 4 MOC计算结果 Fig. 4 MOC calculation results


Fig. 5
图 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。

Fig. 6
图 6 FVM与不同网格数的MOC对比图 Fig. 6 Comparison of FVM and MOC with different grid numbers


表 1
表 1 各模型计算时间 Tab. 1 Computational time of models ? s

网格数 MOC计算时间 FVM计算时间

64 0.072 0.310

128 0.159 0.992

256 0.386 3.747

512 1.080 14.118



表 1 各模型计算时间 Tab. 1 Computational time of models ?


2.2 某抽水蓄能电站水力瞬变计算分析 2.2.1 电站参数电站设有两台150 MW的可逆式机组，上游引水系统采用“一管双机”、下游尾水系统采用“一管一机”的布置方式，布置如图 7所示，相关设计参数见表 2、3。

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


表 2
表 2 管道参数 Tab. 2 Parameters of water pipe system 管道 直径/m 长度/m 糙率 倾角/(°) 波速/(m·s^-1)

L1 8.00 15.39 0.014 0 976.4

L2 8.00 169.26 0.014 -50.71 976.4

L3 8.00 20.77 0.014 0 976.4

L4 8.00 56.40 0.015 0 976.4

L5 8.00 26.60 0.014 0 976.4

L6 4.20 100.33 0.013 0 1 202.3

L7 4.55 5.40 0.013 0 1 210.8

L8 6.09 14.00 0.014 0 1 045.1

L9 6.20 70.94 0.014 0 1 045.1

L10 6.20 25.52 0.014 60.00 1 152.75

L11 6.20 13.60 0.014 0 1 152.75



表 2 管道参数 Tab. 2 Parameters of water pipe system


表 3
表 3 电站机组参数 Tab. 3 Unit parameters of power station 机组参数 额定水头(扬程)/m 额定流量/(m³·s^-1) 额定转速/(r·min^-1) 额定功率/kW 转动惯量/(t·m²)

常规机组 105.8 148.8 200 139 000 10 920

水泵机组 126.7 154.8 200 157 900 10 920



表 3 电站机组参数 Tab. 3 Unit parameters of power station


2.2.2 计算结果在机组甩荷过程中，机组转速的波动是造成水轮机无法稳定运行最直接原因。因此将MOC和二阶Godunov格式的FVM计算出的转速与甩荷实验值进行对比。但是由于MOC在Cr < 1.0时会产生较严重的数值衰减，因此需要对管道内波速进行调整，来满足Cr=1.0的条件。调整后的波速与管道分段数见表 4。

表 4
表 4 管道系统各部分管段波速及其分段数 Tab. 4 Section number and wave speeds of pipe system 管道 MOC FVM

C/(m·s^-1) N Cr C/(m·s^-1) N Cr

L1 961.875 4 1.000 976.400 3 0.761

L2 984.070 43 1.000 976.400 43 0.992

L3 1 038.500 5 1.000 976.400 5 0.940

L4 1 007.143 14 1.000 976.400 14 0.969

L5 950.000 7 1.000 976.400 6 0.881

L6 1 194.405 21 1.000 1 202.300 20 0.959

L7 1 350.000 1 1.000 1 210.800 1 0.897

L8 1 166.667 3 1.000 1 045.100 3 0.896

L9 1 043.235 17 1.000 1 045.100 16 0.943

L10 1 063.333 6 1.000 1 152.750 5 0.903

L11 1 133.333 3 1.000 1 152.750 2 0.678



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


本文选取了3组实验计算工况，已知3组实验工况的输入功率、初始流量与初始转速均为额定值，其他相关参数见表 5。

表 5
表 5 机组100%甩荷计算工况参数 Tab. 5 Calculation parameters of 100% load rejection 实验工况 上库水位/m 下库水位/m 初始导叶开度/% 第1段关闭时长/s 第2段关闭时长/s 拐点导叶开度/%

1 412.40 290.97 73.80 3.62 32.53 60.00

2 406.08 290.60 74.30 5.64 30.53 61.13

3 404.70 290.38 78.93 4.80 33.00 59.96



表 5 机组100%甩荷计算工况参数 Tab. 5 Calculation parameters of 100% load rejection


两种方法在整个计算时间段内计算的转速以及蜗壳处水头的变化情况见表 6以及如图 8~10所示。

表 6
表 6 机组100%甩荷工况结果 Tab. 6 Calculation results of 100% load rejection ? (r·min^-1)

实验工况 实验最大转速 MOC计算转速 FVM计算转速

1 283.84 281.417 281.642

2 279.20 278.401 278.612

3 283.08 279.006 279.005



表 6 机组100%甩荷工况结果 Tab. 6 Calculation results of 100% load rejection ?


Fig. 8
图 8 甩荷实验1 Fig. 8 Load rejection experiment 1


Fig. 9
图 9 甩荷实验2 Fig. 9 Load rejection experiment 2


Fig. 10
图 10 甩荷实验3 Fig. 10 Load rejection experiment 3


由图 8(a)、图 9(a)、图 10(a)所示，MOC与FVM在对该抽水蓄能电站模型进行甩荷计算时，计算结果与实验值基本吻合，表明二阶Godunov格式的FVM在处理较复杂管网模型时的适用性和准确性。由表 6中数据与图 8(a)、图 9(a)、图 10(a)所示，MOC与FVM计算的最大转速均略小于实验值，但FVM计算结果更加接近于实验值，这是因为MOC在计算过程中，需要对管道波速进行调整来满足库朗数条件，由此带来了一定的误差，且调整波速间接地增加了计算时间，导致其计算效率的降低。而FVM仅适当降低库朗数条件，无需进行波速调整，简化了计算过程，具有较高的计算效率与精度。

由图 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仅适当降低库朗数条件，无需进行波速调整，简化了计算过程，具有较高的计算效率和精度。


参考文献
[1] GUITTET M, GAPEZZALI M, GAUDARD L, et al. Study of the drivers and asset management of pumped-storage power plants historical and geographical perspective[J]. Energy, 2016, 111: 560. DOI:10.1016/j.energy.2016.04.052


[2] STEFFEN B. Prospects for pumped-hydro storage in Germany[J]. Energy Policy, 2012, 45: 420. DOI:10.1016/j.enpol.2012.02.052


[3] 黄伟, 杨开林, 郭新蕾, 等. 抽水蓄能电站极端甩负荷工况球阀协同调节[J]. 清华大学学报(自然科学版), 2019, 59(8): 635.
Huang Wei, Yang Kailin, Guo Xinlei, et al. Coordinated regulation of ball valves in pumped storage stations for extreme conditions[J]. Tsinghua Univ(Sci & Technol), 2019, 59(8): 635. DOI:10.16511/j.cnki.qhdxxb.2019.22.013


[4] TRIVEDI C, CERVANTES M J, GANDHI B K, et al. Transient pressure measurements on a high head model Francis turbine during emergency shutdown, total load rejection, and runaway[J]. Journal of Fluids Engineering, 2014, 136(12): 121107. DOI:10.1115/1.4027794


[5] GUINOT V. Riemann solvers for water hammer simulations by Godunov method[J]. International Journal for Numerical Methods in Engineering, 2000, 49(7): 851. DOI:10.1002/1097-0207(20001110)49:73.0.CO;2-#


[6] ZHAO Ming, GHIDAOUI M S. Godunov-type solutions for water hammer flows[J]. Journal of Hydraulic Engineering-ASCE, 2004, 130(4): 341. DOI:10.1061/(ASCE)0733-9429(2004)130:4(341)


[7] 郑劼恒, 蒋明, 郭芮, 等. 顺序输送管道水力瞬变模拟的有限体积法[J]. 计算力学学报, 2015, 32(3): 418.
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. DOI:10.7511/jslx201503019


[8] ZHOU Ling, WANG Huan, LIU Deyou, et al. A second-order Finite Volume Method for pipe flow with water column separation[J]. Journal of Hydro-environment Research, 2017, 17: 47. DOI:10.1016/j.jher.2016.11.004


[9] 赵越, 周领, 刘德有, 等. 基于有限体积法Godunov格式的水锤计算模型[J]. 水利水电科技进展, 2019, 39(1): 76.
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. DOI:10.3880/j.issn.1006-7647.2019.01.013


[10] WYLIE E B, STREETER V L. Fluid transients in systems[M]. New York: Prentice Hall, 1993.


[11] SUTER P. Representation of pump characteristics for calculation of water hammer[J]. SulzerTechnical Review, 1966(11): 45.


[12] 张健, 郑源. 抽水蓄能技术[M]. 南京: 河海大学出版社, 2011.
ZHANG Jian, ZHENG Yuan. Pumped storage technology[M]. Nanjing: Hohai University Press, 2011.


[13] 杨开林. 电站与泵站中的水力瞬变及调节[M]. 北京: 中国水利水电出版社, 2000.
YANG Kailin. Hydraulic transient and regulation in power station and pumping station[M]. Beijing: China Water & Power Press, 2000.


[14] 陈乃祥. 水利水电工程的水力瞬变仿真与控制[M]. 北京: 中国水利水电出版社, 2005.
CHEN Naixiang. Hydraulic transient simulation and control of water conservancy and hydropower Engineering[M]. Beijing: China Water & Power Press, 2005.


[15] TORO E F. Riemann solvers and numberical methods for fluid dynamics[M]. 3rd ed. Berlin: Springer-Verlag, 2009: 115.