抽水蓄能电站球阀自激振系统参数辨识与特性研究

李佰霖, 方子铃, 贵树鑫, 王云峰

(1.三峡大学 电气与新能源学院,湖北 宜昌 443002;2.三峡大学 梯级水电站运行与控制湖北省重点实验室,湖北 宜昌 443002;3. 福建仙游抽水蓄能有限公司,福建 莆田 351267)

水力振动是抽水蓄能电站中一种常见的危害机组安全的现象,水力振动产生的压力远超管道正常运行的安全范围,会造成爆管等事故,严重威胁抽水蓄能电站的安全运行。索丽生等[1]曾通过实例对水电站水力振动的发生进行了分析与预测;周建旭等[2-3]基于改进的水体弹性模型揭示了压力管道的水力振动特性;段炼等[4-6]将矩阵法应用于水电站的水力振动研究,并通过实例分析了该方法的可靠性。除了以上方法,胡光华等[7]基于特征线法对不同波动的耦合效应进行了分析,刘冬等[8]在此基础上提出了改进的特征线法;楼勇等[9]则通过CFD(computational fluid dynamics)建立三维模型的方式得到了自激振动时的压力与流量规律;程加堂等[10]将粒子群BP(back propagation)神经网络算法应用于水电机组振动的故障诊断。

球阀在抽水蓄能电站中起到控制流量与稳定压力的作用,当球阀密封劣化漏水会导致球阀变成柔性阀而发生自激振动[11]。张绍春[12]对柔性阀的自激振动进行分析得到在有分叉管时,自激振动的最大压力可以达到静水压力的3倍;朱渊岳等[13-14]通过特征线法(method of characteristics,MOC)和阻抗法计算的仿真结果得出了自激振动的原因,并提出了消除自激振动的措施;李红辉等[15]建立了自激振动的数值仿真模型[16],并通过经验小波变换对模拟自激振动信号和实际信号的频谱进行分析,以验证模型的可行性;斯静等[17]设计了两个对比方案,得出了管道与阀门漏水量对自激振动振幅和周期的影响。除此之外,有些学者则从机理角度对自激振动进行了分析, Bouzidi等[18]研究了弹簧加载阀的自激振动机制,发现振动幅度与管道长度和弹簧刚度呈正相关。Awad等[19]建立了一个简单的球阀模型来研究环形密封的自激振动机理,并分析了各种因素对模型稳定性的影响。

基于以上研究现状,过往学者对自激振动原理进行了详细研究缺乏与现场实际相结合,为此本文通过特征线法建立了自激振动系统的综合仿真模型,研究框架如图1所示。期望通过物理系统和数字模型的完全结合来分析球阀的自激振动特性,从而对现场工作进行指导,这有利于减少自激振动对水力设备的破坏,维护抽水蓄能电站的正常运行。

图1 球阀自激振动研究框架Fig.1 Research framework for self-excited vibration of ball valve

1 建 模

某抽水蓄能电站的结构如图2所示。本节选择图2中的虚线部分进行建模。该结构主要包括上游水库、主管道L1以及分叉管道L2,L3、阀门1和阀门2。在该模型中,自激振动时由于阀门2处密封漏水引起的,因此阀门2被视为柔性阀门。阀门1是关闭的,因此阀门1被视为封闭端。

图2 某抽水蓄能电站的引水结构Fig.2 The water diversion structure of a pumped storage power station

1.1 柔性阀门建模

抽水蓄能电站球阀根据漏水特性可分为刚性阀门和柔性阀门,通过分析活动环的受力情况可以判断它能否正常与球阀表面贴合来确定阀门的漏水特性。活动环受到水压力的合力如下

F=P2S2+P3S3-P1S1

(1)

式中:F为活动环受到水压力的合力;P1为输入腔的水压;P2为阀体内部水压;P3为输出腔的水压;S1为输入腔的水压对活动环的压力面积;S2为阀体内的水压对活动环的压力面积;S3为输入腔的水压对活动环的压力面积;S1=S2+S3。

工作密封的结构如图3所示。刚性阀门漏水特性分析:正常情况下密封不会发生漏水。由于输入腔的水是从球阀的上游侧管道引入的,因此输入腔的水压等于阀体内部水压。输出腔内的水是常压水,其压力远小于输入腔的水压,并且输入腔的水压对活动环的压力面积大于阀体内部水压对活动环的压力面积。活动环受到的水压力将其向左推动,使其与球阀表面贴合。

图3 工作密封结构Fig.3 Structure of working seal

柔性阀门漏水特性分析:当输入腔与输出腔之间或阀体与输出腔之间的密封泄漏时,输出腔的水压也等于球阀的上游侧管道压力,活动环受到的水压力处于平衡状态。活动环不会向阀体移动,也无法完成工作密封的投入。如果此时出现微小的压力扰动,使阀门漏水面积减小,球阀上游侧的管道压力会随之增大。阀体、输入腔和出输出腔中的水压也会随着上游水压的增加而增加。由于输出腔的水压变化是由密封漏水引起的,因此它的水压变化是一个缓慢的过程,具有一定的滞后性。因此,输出腔的水压小于输入腔的水压。此时,活动环受到水压力向左,阀门的漏水面积进一步减小。球阀上游水压的持续增加导致活动环继续受到向左的力的影响。这形成了一个正反馈过程。

半个周期后,从上游水库反射的水锤波导致阀门上游侧压力减小,阀体、输入腔和输出腔中的水压也随之减小。由于输出腔的水压变化具有一定的滞后性,此时输出腔的水压大于输入腔的水压。活动环受到的水压力向右,导致阀门漏水面积增大,这也是一个正反馈过程。水锤波在管道中来回传播使这一过程不断重复,这会导致压力振幅不断增大最终发生自激振动。

柔性阀门的漏水面积变化如图4所示。柔性阀的漏水面积与水压之间的关系为

图4 柔性阀门漏水面积Fig.4 Leakage area of flexible valve

(2)

式中:H0为初始压力;A1为静压下的漏水面积,H1为密封处水压。

柔性阀门与刚性阀门的漏水特性如图5所示。刚性阀门的漏水量和阀门前后的压差的关系式如下

图5 刚性阀门和柔性阀门的漏水特性Fig.5 Leakage characteristics of rigid valve and flexible valve

(3)

式中:k为常数;Q为阀门的漏水量;H为阀门前后的压差。

柔性阀门的漏水量和阀门前后的压差的关系式如下

(4)

式中:D1为阀门直径;x0为阀门前后压差为0时的密封间隙;Af为阀门密封盖的面积;K为阀门密封的弹性系数;ρ为水的密度。

由于密封的压力水是从球阀上游侧的管道引入的,因此H1=H。因此,柔性阀门漏水量可以表示为

(5)

1.2 特征线法

以水锤基本方程为基础,通过特征线法将水锤基本方程转化为有限差分方程进行计算。特征线法的矩形计算网络如图6所示。将管道划分为每段长为Δx,将时间被划分为每个时间段Δt=Δx/a。C+为正特征线,C-为负特征线,i为管道节点。计算从t=0开始。除边界点外的其他节点的参数可以通过特征线方程计算获得,边界点可以通过特征线方程结合边界条件计算获得。

图6 特征线法的矩形计算网络Fig.6 Rectangular calculation network of MOC

特征线方程

C+∶HPi=CP-BQPi

(6)

C-∶HPi=CM+BQPi

(7)

CP=Hi-1+BQi-1-RQi-1|Qi-1|

(8)

CM=Hi+1-BQi+1+RQi+1|Qi+1|

(9)

上游边界条件

HP=Hu

(10)

式中,Hu为上游水库水头。

下游边界条件

(11)

2 参数优化

由于阀门漏水量无法准确测量,因此很难获得柔性阀门的特征参数k1和k2,无法精确模拟该抽水蓄能电站球阀自激振动过程。为了使模拟的自激振动过程符合实际情况,需要对柔性阀门的特性参数进行优化。

2.1 适应度函数

为了找到k1和k2的最优解本文将仿真结果与监测数据相结合得到了反映特征参数优良程度的适应度函数。由于监测数据的采样率较低,很难充分反映压力随时间的变化过程。如果通过比较每个时间点的数据来评估仿真结果和监测数据之间的一致性就会产生很大的误差。为了使比较结果更加准确,本文对监测数据与仿真结果的上下包络线进行了比较,以确定实际和模拟自激振动发展过程的一致性。形成监测数据包络线方程的点如下

(12)

式中:HP(t) 为压力数据;Hm为上包络线的点;Hn为下包络线的点;t为自激振动时间。

为了使曲线更平滑,对曲线进行拟合。由于曲线的函数类型尚未确定,因此可以通过泰勒级数展开将曲线的表达式转换为多项式函数,最小二乘法通过计算最小误差的平方和来找到拟合多项式函数的最佳函数匹配。拟合后上下包络线的函数表达式为

(13)

式中:Hm1(t) 为监测数据优化后的上包络线;Hn2(t) 为监测数据优化后的下包络线。

仿真结果的包络线是通过相同的方法求得,将监测数据与仿真结果的包络线进行了比较,并采用时间乘绝对误差积分准则作为性能指标,得到反映仿真结果和监测数据一致性程度的函数如下

(14)

式中:Hm2(t) 为仿真结果的上包络线;Hn2(t) 为仿真结果的下包络线。

由于不同的柔性阀门特性参数会导致不同的阀门漏水量,需要限定漏水量的范围,以排除不符合要求的特征参数。从式(2)中可以得出,初始漏水面积为A1,漏水面积在0和2A1之间来回变化。因此,0

(15)

式中,Q0为阀门的初始漏水量。

对于不满足式(15)的柔性阀特征参数,在性能指标中加入惩罚函数。惩罚函数如下

(16)

获得适应度函数如下

y(j)=F1(t)+ε(t)

(17)

2.2 粒子群优化

采用粒子群算法进行优化,粒子的速度和位置更新公式如下

(18)

式中:xj为粒子位置;vj为粒子的更新速度;w为惯性因子;c1和c2为加速度常数;pbestj为局部最优位置;gbestj为全局最优位置。

将基于MOC的自激振动模型代入PSO进行参数优化的流程如图7所示。主要步骤如下:

图7 参数优化流程图Fig.7 Parameters optimization flow chart

步骤1通过设置适当的系统参数并限定粒子的范围随机生成粒子。

步骤2将每组粒子代入MOC程序得到仿真结果,通过现场状态监测装置获得监测数据,通过适应度函数计算每组粒子的适应度。

步骤3将初始适应度作为每个粒子的局部最优解,并找到全局最优解。

步骤4通过速度和位置更新公式,更新粒子的速度和位置。超出范围的粒子将根据边界条件进行处理。

步骤5计算每个新粒子的适应度与局部最优解比较,找到新的局部最优解,并从所有新的局部最优解中找到全局最优解。

步骤6重复步骤4和步骤5,直到找到最优解。

抽水蓄能电站初始水头为516 m,具体管道参数如表1所示。监测数据的拟合曲线函数为4阶,参数为a1=52.38,a2=-3.504×10-3,a3=6.581×10-5,a4=-2.323×10-7,a5=3.776×10-10,b1=50.9,b2=4.508×10-3,b3=-6.935×10-5,b4=2.321×10-9,b5=-3.7×10-10。粒子群优化算法设置50个粒子并迭代30次;w=0.8;c1=1.5;c2=1.5。

表1 管道参数Tab.1 Pipeline Parameters

代入数据计算后,适应度函数曲线如图8所示。在20次迭代之后,曲线趋于平缓,结果达到稳定。求得柔性阀的特征参数为k1=0.039 8,k2=3.19×10-5。

图8 适应度函数曲线Fig.8 Fitness function curve

优化前后仿真结果与监测数据的对比如图9所示。可以看到优化前的压力数据与实际情况存在较大差距;优化后仿真结果在400 s前与实际自激振动过程存在较小误差,在400 s后的发展过程基本相同。通过粒子群算法优化的柔性阀门特征参数符合实际情况。

图9 优化前后球阀上游侧压力对比Fig.9 Comparison of upstream side pressure of ball valve before and after optimization

3 分 析

在获得与实际情况一致的模型后,本章通过控制变量法分析了初始扰动,管道长度,管道直径,双柔性阀门系统对自激振动的影响,确定自激振动特性为抽水蓄能电站的优化设计提供了依据。

3.1 初始压力扰动对自激振动的影响

仅改变阀门1处受到的初始压力扰动,不改变其他参数,设计对比方案1,其球阀上游侧压力如图10所示,统计不同初始扰动其振幅时间如表2所示。

表2 不同压力扰动的振幅时间Tab.2 Amplitude time of different pressure disturbances

图10 方案1对比结果Fig.10 Comparison results of option 1

如表2所示:从1 m水头到2.5 m水头在两种不同初始扰动情况下所需时间相同,均为105 s;从2.5 m水头到400 m水头在3种不同初始扰动情况下所需时间相同,均为600 s。在不同的初始压力扰动下,自激振动从一个振幅到另一个振幅的增长速度是相同的。对比3种不同初始扰动的发展过程,可以发现自激振动在初始振幅较小的阶段缓慢增加,在后期振幅较大时增加迅速,因此初始扰动越大,自激振动增长速度也就越快。

3.2 管道长度对自激振动的影响

为探究自激振动与不同管道长度的关系,设计对比方案2,仅改变主管道长度与分叉管道长度,其他条件不变,共设计了9种管道长度组合,其球阀上游侧压力如图11所示,统计不同管道长度组合的振幅与频率如表3所示。

表3 不同管道长度组合的振幅与频率Tab.3 Amplitude and frequency of different pipeline length combinations

图11 方案2对比结果Fig.11 Comparison results of option 2

在前6种管道长度组合中仅改变主管道长度,球阀上游侧压力信号的主频率仍为1.67 Hz,而后3种管道长度组合中仅改变分叉管道长度,自激振动压力信号的主频率随着分叉管道长度的增加而减小,这是由于水锤波在异径点处发生了反射,因此自激振动频率主要与分叉管道长度有关。

仅增大主管道长度,自激振动幅值并没有随之增加,由幅值结果可以判断自激振动与不同管道长度的比例有关。在复杂管道的情况下,假设每条管道的水锤波速相等,当主管道长度是分叉管道长度的偶数倍时的自激振动的增长速度远大于当主管道长度不是分叉管道长度偶数倍时的增长速度。这是由于水锤波在管道中来回流动时,会在不同管道的异径点处发生透射与反射。当主管道长度是分叉管道长度的偶数倍时,不同频率的水锤波之间的叠加效应更明显;当主管道长度不是分叉管道长度的偶数倍时,不同频率的水锤波之间的抵消效应更明显。仅增大分叉管道长度,自激振动的频率减小,这就导致了自激振动的增长速度也随之减小。

3.3 管道直径对自激振动的影响

为探究管道直径对自激振动的影响,设计对比方案3,仅改变主管道直径与分叉管道直径,其他条件不变,其球阀上游侧压力如图12所示,不同管道直径的自激振动幅值如表4所示。

表4 不同管道直径的自激振动幅值Tab.4 Self excited vibration amplitudes of different pipeline diameters

图12 方案3对比结果Fig.12 Comparison results of option 3

当仅改变主管道直径时,与初始状态相比自激振动幅值基本未发生改变,因此主管道直径对自激振动的增长速度影响不大;当仅改变主管道直径时,分叉管道与阀门直接相连,分叉管道直径变化会导致阀门附近的水锤波变化从而直接影响自激振动的发展过程,因此从仿真结果可以看到自激振动的增长速度随着分叉管道直径的增大而减小;同时减小或增大主管道和分叉管道直径的自激振动结果与以上分析结果相同。

3.4 双柔性阀门系统对自激振动的影响

当两处阀门均发生密封漏水变为柔性阀门时,设计对比方案4,将阀门1由刚性阀门变为柔性阀门,其他管道等参数不变,设计了3种情况下的自激振动仿真,阀门2上游侧压力如图13所示,3种情况下的自激振动幅值如表5所示。

表5 双柔性阀门系统不同情况下的自激振动幅值Tab.5 Amplitude of self-excited vibration of double flexible valve system under different conditions

图13 方案4对比结果Fig.13 Comparison results of option 4

当系统中的两个阀门均为柔性阀门时,其他条件均不发生改变,其仿真结果如情况1所示,300 s时的自激振动幅值明显增大,这是由于有两个阀门同时发生振动加速了自激振动的发展。在柔性阀门1处增加1个与柔性阀门2处相同的扰动,其仿真结果如情况2所示,300 s时的自激振动幅值仅有14 m,这是由于两处阀门上的扰动产生的水锤波在管道中传播时相向运动发生了抵消导致自激振动增长缓慢。在柔性阀门1处增加一个与柔性阀门2处相反的扰动,其仿真结果如情况3所示,300 s时的自激振动幅值与情况1相比进一步增大,这是由于两处阀门扰动产生的水锤波在管道中传播时同向运动发生了叠加导致管道中的扰动进一步增大,使得自激振动幅值也快速增大。当柔性阀门扰动发生的时间不同时,其仿真结果如情况4、情况5所示,由于扰动有时差,水锤波在传播时出现抵消或者叠加,最终导致自激振动结果不同。

4 结 论

为研究抽水蓄能电站球阀自激振动的发展过程本文提出了一种球阀自激振动的研究框架,通过物理系统的监测数据来优化数字模型的建立,并分析数字模型的自激振动特性,以指导物理系统的工作。在优化过程中,本文基于MOC和PSO对模型参数进行优化。通过将优化后仿真结果与监测数据进行比较,验证了优化参数的可靠性。最后通过两个对比方案发现了自激振动特性,结论如下:

(1)初始扰动越大,自激振动增长速度也就越快。

(2)由于水锤波在异径点处的反射,分叉管道长度会影响自激振动频率,同时又会影响自激振动的发展过程。

(3)不同管道的水锤波速相同时,由于水锤波在异径点处的透射与反射,当主管道长度是分叉管道长度的偶数倍时自激振动的增长速度远大于当主管道长度不是分叉管道长度偶数倍时的增长速度。

(4)分叉管道直径会影响阀门附近的水锤波,从而导致自激振动的增长速度随着分叉管道直径变化而变化。

(5)当系统中的两个阀门均为柔性阀门时,两处阀门处的扰动产生的水锤波会在管道中发生抵消或叠加进而导致自激振动增长速度的增大或减小。