抽水蓄能电站尾水管水柱分离数值模拟

肖鹏飞，赵志高，杨桀彬，刘程鹏，何相慧，杨建东

［1．中国长江三峡集团有限公司科学技术研究院，北京市 101100；2．水资源工程与调度全国重点实验室（武汉大学），湖北省武汉市 430072］

0 引言

随着风光新能源的发展，抽水蓄能成为新型电力系统中的基础性调节电源。根据《抽水蓄能中长期发展纲要》，中国抽水蓄能装机将翻两番，达到1.2 亿kW。在资源储备站点日益紧张的背景下，抽水蓄能向着高水头、大容量、长输水管道方向发展。抽水蓄能机组复杂的流道布置、频繁的工况切换[1]可能导致水力系统出现较大的水锤振动、液柱分离甚至空化空蚀现象[2]。如何准确模拟尾水管液柱分离对判别抽水蓄能运行安全性、确定机组安装高程具有重要意义。

国内外研究学者在20 世纪80 年代就对伴有空穴流、液柱分离的气液两相瞬态过程展开了研究。杨建东[3]从空化和气液两相流理论两个方面，采用理论与实验相结合的方法，系统地研究了空化流动与液核分离的瞬态过程，为今后进一步研究其瞬态过程的机理、数值分析和满足实际需要奠定了基础。蒋劲等[4]给出了一种考虑两相间物质质量交换、管壁流质相互作用的一维计算模型。周领等[5]将Wylie[6]提出的离散空腔模型（DGCM）和有限体积法Godunov 格式结合起来计算水柱分离。杨桀彬[7]在集中空穴模型和气泡均匀分布模型的基础上提出了以临界空穴率作为水柱分离的充分条件的气液两相流计算模型。何相慧[8]在抽水蓄能机组动态特性实验平台上开展了水柱分离模型实验，给出了尾水管发生水柱分离的判别条件。

本文基于有限体积法（FVM），引入空穴率α描述气体体积，采用二阶Godunov 格式离散双曲型偏微分方程，对含气型气液两相流以及相变下的瞬态过程进行了数值模拟，尤其是摒弃了Wylie 离散空腔模型，直接采用空化模型封闭方程组，并通过工程实例与传统计算方法进行对比。

1 气液两相流及空化模型

1.1 基于有限体积法的气液两相流模型

为了衡量管道中含气量变化，引入量度气液两相流中气泡含量大小的参数，容积空穴率α，其表达式为：

式中：W——体积；

g和l——气相和液相；

L——管道长度；

A——段面横截面面积。

气液两相瞬变流的基本方程仍然是动量方程和能量方程，应对气相和液相分别列出。对于含气型气液两相瞬变流中气相体积较小（α<0.15）的情况，气泡会跟随液体一起运动，于是忽略动量交换，将动量方程合并为一个。而连续方程仍保留分相的形式，以便计算气液两相之间的质量交换[9]。

类似单相流体动量方程推导方法，可得到含气型气液两相瞬变流动量方程：

式中：p——管道中心点压强，Pa；

V——流速，m/s；

x——从任意点开始的沿管道轴线的坐标距离，m；

t——时间，s；

α——空穴率；

ρg——气相密度，kg/m3；

ρl——液相密度，kg/m3；

D——管道直径，m；

f——达西-威斯巴哈（Darcy-Weisbach）系数；

g——重力加速度，m/s2。

连续性方程的本质是质量守恒方程，为考虑两相之间的质量交换，将两相间的质量变换率=dm/dt引入单相流的连续方程，可以分别写出气相和液相的连续方程：

对上述方程进行变形再通过有限体积法进行离散，得到计算表达式：

式中：Fi+12(UL,i+12,UR,i+12)，Fi-12(UL,i-12,UR,i-12)——界面i+1/2，i-1/2的数值通量；

UL和UR——界面左右的变量值。

1.2 变波速模型与相变模型

在参考文献[9]中，给出了波速方程：

式中：a——水锤波波速，m/s；

C1——管道固定方式系数；

E——管道弹性模拟量；Kg和Kl——气体容变弹性模拟量和液体容变弹性模量。

对于空化现象的数值模拟，需要定量描述气液之间的物理相变过程，关键在于建立或选取合适的空化模型。空化模型是用于描述气相和液相之间质量交换的数学模型，目前应用较为广泛的空化模型主要分为两类：

（1）基于正压流体状态方程的空化模型；

（2）基于质量输运的空化方程[10]。

本文采用的计算模型需要定量考虑气液两相之间的相变关系，故选取基于质量输运的空化模型。有压管道的空化现象，往往伴随着气泡的破碎和融合，在极短的时间内就会发展为空化空腔，因而基于单个气泡的空化模型也不适用于方程的计算和求解。最后基于相变物理本质，选取了以Singhal 模型为基础的空化模型进行数值模拟分析，并根据气液两相管流实际情况对公式参数进行了优化。

式中：Vch——特征速度，它反映了液体和蒸汽之间局部相对速度的影响，通常认为大空隙不移动，所以Vch=V；

σ——表面张力系数。

在该模型中，特征速度、表面张力系数、不可凝结气体的影响也被考虑在内，因而被称为 “全空化模型”。

2 有限体积法和特征线法的耦合

2.1 计算模型

在抽水蓄能管道系统中，可能发生空化的区域较少，实际上一般只发生在尾水管进口。其中尾水管发生的空化和液柱分离会对水轮机和水力系统造成极大影响。在实际工程中，往往降低安装高程以防止尾水管内的压力降至蒸汽压力以下[11]。而需要重点关注的也是尾水管部分的空化现象。

因此，同时考虑尾水管部分空化模拟的准确性和计算效率，提出了两个解决方案。

（1）将FVM 和MOC 方法进行耦合，采用FVM 计算尾水管，MOC 方法计算抽水蓄能电站其余水力单元。

（2）尾水管的FVM 模拟计算和其他部分的MOC 模拟计算采用不同的时间步长进行耦合。

基于Riemann 不变量，有限体积法在边界上有以下方程：

结合特征线法基本边界方程：

如图1 所示，对于MOC 下游边界和FVM 上游边界，也就是边界1 耦合，通过式（9）、式（10）以及边界方程联立可求解。类似地，对于FVM 下游边界和MOC 上游边界，也就是边界2 耦合，通过式（8）、式（11）以及边界方程联立可求解。

图1 边界耦合示意图Figure 1 Schematic diagram of boundary coupling

2.2 方法验证

如图2 所示，采用串联点作为边界耦合点，对于上库-阀门简单的管道系统，压力管道L1 和压力管道L2 的长度均为1000m，糙率n=0.014。

图2 水库-阀门系统Figure 2 Reservoir-valve system

上游库水位20m，阀门流量由4m3/s 在0 时刻减为0m3/s。

对于管道L1 采用FVM 进行数值模拟计算，对于管道L2采用MOC 进行数值模拟计算。方案一FVM 和MOC 时间步长均为0.04s；方案二FVM 时间步长为0.04s，MOC 时间步长为0.2s。

阀门处测压管水头随时间变化值如图3 所示。显然，在不同时间步长下耦合计算结果较为吻合，说明该耦合方法可用于管道瞬态过程模拟。

图3 阀门测压管水头Figure 3 Valve gauge head

3 工程应用

3.1 工程计算模型

以某抽水蓄能电站为例，该电站为一洞四机方案，设有下游调压室，管道示意图如图4 所示。机组号为1 号、2 号、3 号和4 号，机组及管道其他参数如表1 所示。

表1 某抽水蓄能电站运行参数Table 1 Operating parameters of a pumped storage power plant

表2 某抽水蓄能电站运行参数Table 2 Operating parameters of a pumped storage power plant

图4 某抽水蓄能电站管线布置示意图Figure 4 Schematic diagram of pipeline layout of a pumped storage power plant

本文研究抽水蓄能电站尾水管空化和液柱分离受安装高程的影响情况。尾水管部分采用FVM 方法建模计算，其他部分采用MOC 方法计算。同时，为了兼顾计算效率和准确性，FVM 部分取0.001s 时间步长计算，MOC 部分取0.01s 时间步长计算。

计算工况为该工程过渡过程大波动尾水管进口压强的最不利工况CT4-5：上库正常蓄水位1262.5m，下库死水位873m，四台发电机组均以额定出力运行。在0时刻，同一水力单元3号、4 号机组甩符荷，1 号、2 号机组5s 后甩负荷，导叶紧急关闭。

3.3 安装高程敏感性分析

为研究空化和液柱分离对于抽水蓄能电站尾水管进口压强的影响，本次数值模拟，选取了6 组不同机组安装高程进行计算，并对1 号机组尾水管进口压强进行分析。计算结果以相对压强对应的压强水头为准。

6 组数值仿真模拟试验计算结果如图5 所示。显然，在不考虑压力脉动和计算误差的情况下，安装高程从初始的805m 上升到825m 时，尾水管进口压强随着安装高程上升而不断下降，波形完全一致。当机组安装高程上升到830m 时，尾水管进口压强会下降到汽化压强以下，进而发生空化。

图5 不同安装高程尾水管进口压强水头Figure 5 Tailpipe inlet pressure head at different installation elevations

3.3 水柱分离模拟对比

图6 给出了初始安装高程为830m 条件下，运用FVM 采用空化模型和不采用空化模型两种方法对1 号机组尾水管进行数值模拟的计算结果。采用空化模型计算，尾水管进口压强达到汽化压强时，压强不再下降，出现空化和液柱分离，之后的最大升压相比初始压强上升62.85m。而不发生空化的数学模型计算结果最大升压相比初始压强上升38.15m，远小于发生液柱分离后的压强升高值。并且发生液柱分离后压强上升速率更快。同时，发生液柱分离后会使得波速降低，进而水锤波动周期向后推移。在导叶完全关闭后，尾水管进口压强趋于稳定，与无空化模型计算结果相似。

图6 尾水管进口压强水头Figure 6 Tail pipe inlet pressure head

图7 单独给出了通过FVM 空化模型计算得到的尾水管进口压强随时间变化曲线图。在0 时刻，3 号、4 号机组导叶瞬时关闭，受此影响，机组尾水管压强缓慢下降。5s 时刻，该机组导叶开始关闭，尾水管压强开始迅速下降，并在9.35s 时刻达到汽化压强，压强停止下降，尾水管进口发生空化，进而出现液柱分离。10.93s 时刻，如图7 中1 点所示，下游增压波到达尾水管进口，液柱弥合，同时液柱分离处出现两个反射波，增压波到达液柱分离处，受空穴反射为减压波向下游传播，而在这同时，增压波继续向机组传播，在机组处反射为增压波向下游传播。受到增压波影响，尾水管压强迅速上升，如图7 中2 点所示。在13.48s 时刻，二次反射增压波到达尾水管进口处，形成液柱分离水锤的二次升压，升压极值为95.21m，远大于不发生空化的升压极值70.52m。

图7 尾水管进口压强水头Figure 7 Tail pipe inlet pressure head

4 结论

本文基于特征线法和有限体积法，引入相变空化模型并耦合边界条件方程，提出了一种抽水蓄能机组水柱分离计算新方法，在不考虑压力脉动和计算误差情况下，实现了对该电站尾水管空化过程的一维数值模拟，并阐明了水柱分离过程中尾水管进口压强的动态演变规律。该方法克服了传统特征线法不能模拟抽水蓄能水柱分离现象的局限性，为抽水蓄能电站安装高程设计及安全稳定运行分析提供了有效工具。