一管双机式混流式水轮机压力脉动预测与振动分析

杨梦娇，陈　端，张昌兵，李金遥

(四川大学水利水电学院， 成都 610065)

0　引　言

水力发电机组安全、稳定运行一直是广大专家和学者重点关注的问题，尤其是机组运行过程中的水力振动问题[1-3]。引起水力振动的原因是多方面的，包括叶轮机械固有的转动部件与静止部件间的动静干扰，叶片绕流产生的卡门涡列，偏离最优工况下尾水管中发展出的旋转涡带[4]等。由于水轮机流道的复杂性，采用实验方法测量其内部流场比较困难，且成本高、周期长。随着计算机和计算流体力学(CFD)的不断发展与完善，应用数值模拟方法计算水轮机内部的复杂流场，预测压力脉动成为一种重要工具，并取得了一系成果[5-8]。杨建明在文献[9]中指出尾水管涡带是造成水力机械在部分负荷时引起机组振动和出力摆动的主要根源。文献[10]中说明水电站过渡过程也可以引起振动。基于振动带来的危害，苏华山[11]等人研究了振动条件下液力耦合器叶轮内部两相流动特性，探索了振动对液力耦合器内外特性影响。曾立飞[12]等人针对调节阀振动问题，考虑阀杆系统振动对阀内流场的影响，进一步揭示由流体诱发振动的原因。庞立军[13]等人针对三峡右岸部分机组振动、产生异常噪音的问题，模拟了固定导叶开展出水边处的涡街振动，确定了产生异常噪声的激振源，并提出有效优化方案。

对于引水管路较长的水电站，为了节约成本，往往采用一根引水管带两台、甚至多台机组的方式。但机组之间不可避免地会出现水力干扰的问题，严重时可能影响机组的稳定运行。文献[14]报道了某电站一根总管带四台水轮机组，机组的振动问题十分突出，大修周期约为2～3年，平均每两年就要更换一个转轮，这严重影响了电站的正常运行及其综合效率。文献[15]对其进行了数值分析，得出了振动主要是由于转轮引起的结论。文献[16]对水轮机进行了尾水管补气数值模拟，有效的补气方式能有效减轻水轮机的水力振动。

由于一根总管带多台水轮发电机组更容易造成水力振动。但到目前为止，水轮机的数值模拟主要是针对单台机组，没有考虑机组之间的水力干扰导致的水力因素，因此，有必要对其水力振动原因作进一步分析。研究表明[17]，标准k-ε模型用于应变率特别大时可能导致负的正应力，为了模拟的流动更加符合湍流的物理定律，Realizablek-ε模型将原来标准k-ε模型中不变的黏度系数Cμ与应变率联系起来，引入与旋转和曲率有关的内容，因而Realizablek-ε模型能更好地模拟射流撞击、分离流、二次流、旋流和圆口射流等复杂的流动问题，且计算结果较准确，精确较高[18]。因此，本文以一管两机混流式水轮机为例，采用时均N-S控制方程组结合Realizablek-ε模型对其水力振动进行数值分析，为机组的稳定运行提供科学依据。

1　数学模型

水轮机的水流按不可压缩流处理。其中岔管、蜗壳、尾水管的水流采用绝对坐标系下控制方程；转轮内的水流采用的旋转坐标系下的控制方程。

1.1　在绝对坐标系统下的控制方程

水轮机的水流视为不可压缩流动，通过岔管、蜗壳、尾水管的水的连续性方程和动量方程如下：

▽V=0

(1)

(2)

式中：V是水流的绝对速度，m/s；t是时间，s；ρ0是水的密度，kg/m3；p是压力Pa；ν和λ分别是第一和第二流体黏度系数。

1.2　在旋转坐标系统下的控制方程

在旋转坐标系下，通过转轮的水的连续性方程和动量方程如下：

▽W=0

(3)

(4)

式中：W是水流的相对速度，m/s；ω是旋转角速度，rad/s；r是半径，m；p是压力Pa。

1.3　Realizable k-ε模型

Realizablek-ε模型的方程如下[19]：

湍动能k输运方程：

湍流耗散率ε输运方程：

(6)

式中：ρ为流体密度；t为时间；k为湍动强度；ε为耗散率；μ为动力黏度系数；v为运动黏度系数；μt为湍流黏性系数；ui、uj分别为流体在xi、xj方向的瞬时流速；Gk为平均速度梯度引起的湍动能k的产生项；S为平均应变率张量系数；G2为常数；σk、σε分别为k、ε的普朗特数。

湍流黏性系数计算：

(7)

上述公式中，取C2=1.9，σk=1.0，σε=1.2。

2　计算方法

2.1　计算域与网格

计算域包括一根岔管和两台相同型号的混流式水轮机。水轮机包括蜗壳、转轮和尾水管，如图1所示。水轮机型号为HLD126-LJ145，固定导叶数8个，活动导叶数Z1=16个，转轮叶片数Z2=14，尾水管形状为弯肘型，其参数如表1所示。

表1　水轮机参数Tab.1 Turbine parameters

本文对水轮机从岔管入口至尾水管出口的所有部件进行全流道三维计算，运用基于压力法的Pressure-Based求解器，动量方程中压力与速度的耦合采用广泛应用的SIMPLE[19]算法。采用非结构化网格，在空间离散上，动量方程采用二阶迎风分格式，湍流量输运方程采用一阶迎风分格式， 压力项离散具有二阶精度。在时间上，采用一阶隐式格式进行离散。

由于结构比较复杂，采用适用性强的非结构化网格进行划分，如图1。转动部件和非转动部件之间采用的交界面(interface)模式为：动静部件耦合，采用瞬态计算的方法。

图1　计算区域及网格Fig.1 Computational domain and mesh

为了保证计算模型能真实模拟水流内部的流动过程，且满足精度要求，在比较复杂的部件如转轮处采用细网格，在相邻部件连接处进行了均匀细化处理，使得网格变化均匀过渡。网格划分总数为250万。

2.2　边界条件

流体状态：流体选用液态水，状态为湍流，采用不可压缩黏性模型，密度为998.2 kg/m3，黏度为0.001 kg/(m·s)。

入口条件：入口采用速度进口边界条件，采用平均流速。

出口条件：出口采用压力出口条件。根据水轮机的吸出高度Hs=-8.0 m，设定尾水管出口压力为8.0 m。

壁面条件：固壁面采用无滑移边界条件(wall)，近壁区采用标准壁面函数来确定其附近区域的流动。

“壁面函数”适用于实体边界：

(8)

绝对和旋转坐标系之间的界面条件：

V=W+ωr

(9)

式中：ωr是圆周速度, m/s.

2.3　计算工况

为了研究不同开度和水头下机组振动的情况，对表2所示工况进行计算。

表2　计算工况Tab.2 Computational operation conditions

为了更好地研究内部水流的变化情况，通过在各过流部件内设置监测点，计算其在不同工况下的压力脉动情况。监测点的位置如图2所示。图2(a)中P1～P5为Z=0平面上的点；图2(b)中P6～P10为1号机组上Z=-0.6 m面上分布的点，其中P6为该截面上的中点，P7～P10为绕P6周向分布、距离中心0.6 m的点；P16为Z=-1.2 m平面上的中点。(P11～P15、P17为2号机组上与1号机组对应位置上的监测点)

图2　监测点位置Fig.2 The location of monitoring points

3　结果与分析

先进行定常计算，当计算收敛时，将定常计算的结果作为非定常计算的初始条件。非定常计算时，依据收敛条件，时间步长应满足Courant-Friedrichs-Lewy(CFL)条件：CCFL=uΔt/Δx, 对于水轮机模拟计算，取CCFL≤1，并由最小网格尺寸和最大流速，时间步长应小于5×10-4s，计算时取1×10-4s。计算过程中，对选定测点位置的压力脉动监测，并进行频谱分析，分析水力振动的原因。

3.1　转轮叶片背面压力分析

以1号转轮为例，其叶片背面压力分布如图3所示。

由图3可知，从叶片进水边至出水边，压力呈减小趋势，转轮叶片靠近上冠处的压力分布比靠近下环位置的压力分布更均匀。叶片背面进水边靠近下环位置处压力梯度较大，由图3(f)可见，该处产生了负压，表明在大开度、高水头工况下，该处水流流速变化大，导致较大的压力梯度，对叶片造成冲击与振动，引起疲劳破坏。下环靠近出水边存在大范围的负压区，该区域容易产生气蚀破坏，影响机组的稳定运行。计算结果表明，在同一水头下，开度越大，下环出水边负压区的范围越大，其表面最大压力值越大。在同一开度下，水头越高，表面的最大负压值越大。2号转轮的叶片压力分布规律与1号转轮叶片压力分布规律大致相同，但2号机组负压区的范围更大，表明两台水轮机的工作状态存在差异。

图3　转轮叶片背面压力分布图Fig.3 Pressure contour on blade back surface of runner

3.2　压力脉动

为了进一步分析压力脉动的变化情况，对各个监测点压力随时间的变化进行快速傅里叶变换(FFT)。通过比较，选择蜗壳进口处的P4、P5，截面Z=-0.6 m 的中心P6、P11和在该截面上离中心0.6 m的点P7、P12几个典型监测点进行分析。图4为工况1下的压力脉动计算结果和对应的频谱分布情况，脉动压力幅值及对应的功率谱见表3。

由图4中压力脉动变化过程和表3中压力脉动幅值可见。蜗壳中的和尾水管中心线的压力脉动周期较短，而尾水管涡带掠过的区域，压力脉动周期相对较长。由表3可见，尾水管中心线处的压力脉动幅值较小，而涡带掠过的区域压力脉动幅值很大，其值是蜗壳里压力脉动幅值的5倍左右，是尾水管中心线处的20多倍。表中监测点4、6、7(属于1号机组)的幅值和功率密度值分别大于对应监测点5、11、12(属于2号机组)。说明1号机组的振动强于2号机组，且振动产生的能量(功率密度的大小代表能量高低)也高于2号机组，这与实际相符合。并且1号机组P6的功率密度的值是对应2号机组P11的值的数倍，在转轮出口处两台机组产生了较大差距。这可能是由于在转轮内，两台机组的水流的流态不同，产生的压力脉动不同，导致了进入尾水管的水流流态不同，也就造成了二者在尾水管能量耗散的差异性。压力脉动的不稳定，引起水流的冲击与振动，故而导致厂房出现振动，产生异常噪声。这正是引起两台机组振动不同的原因。

图4　压力脉动时域图和功率密度频域图Fig.4 Time process of the pressure pulsation and power spectrum density (PSD)

表3　压力脉动的幅值及最大功率谱密度Tab.3 Amplitude of pressure pulsation and maximum power spectral density

注：压力脉动幅值表示峰-峰值。

由频谱分析结果可知，发现两台机组对应监测点的主频基本一致，说明相同位置受影响的主要因素相同。蜗壳压力脉动的主频为140.0 Hz，这与水轮机转轮叶片的固有频率f2=nZ2/60=140 Hz相同，表明蜗壳里压力脉动主要受水轮机转轮叶片动静干涉的影响，另外也包含一些低频成分；尾水管中心线处压力脉动频率主要为10.0 Hz和140.0 Hz，与水轮机固有频率f1=n/60=10 Hz和转轮叶片固有频率相同， 主要受到转轮的影响。在尾水管涡带掠过的区域，压力脉动的主频主要为10.0 Hz及其倍频。

由表3中的最大功率谱密度值来看，尾水管中心线位置的功率谱密度值很小，表明该处的压力脉动所含能量较小。蜗壳里的最大功率谱密度值在10.0 MPa2·s左右，而尾水管涡带区域的最大功率谱密度值很大，达到180.0 MPa2·s，是蜗壳里最大功率谱密度值的18倍，表明尾水管涡带单位时间所含振动能量较大，这是造成机组振动的主要因素。

3.3　尾水管涡带

几种典型工况下的尾水管涡带如图5所示，图5中左边为1号机组，右边为2号机组。

图5　稳定时刻尾水管内三维涡带图Fig.5 Three-dimensional Vorex rope in the draft tube at a certain stabilization time

当机组偏离最优工况区运行时，由于转轮叶片出口产生圆周方向的速度分量，会引起尾水管内水流的旋转，产生涡带，进而产生能量损失，影响机组的稳定工作[20]。由图5可知，工况不同，尾水管涡带的形状有差别。同一台机组，在同一开度下，随着水头的增加，流量增大，涡带逐渐减小，涡带强度减弱，能量耗散减小，振动减小。在设计工况(即工况7)下无漩涡。比较两台水轮机，在相同水头、相同开度下，尾水管进口里的涡带是相似的，对应的功率谱密度值也基本相等。

因此可以发现，两台不同的机组，在相同水头、相同开度下，两台联合运行机组的尾水管进口处的涡带差别并不明显，可以推测尾水管处涡带产生的振动并不是两台机组产生差异的症结所在。

4　结论与展望

本文通过对一根引水管带两台水轮机的非恒定流场进行了三维数值预测，可以得到如下结论：

(1)采用不可压缩流N-S方程结合real Realizablek-ε模型较好地模拟了一管带双水轮机的流场，得到了水轮机的压力分布和涡带分布。

(2)转轮叶片背面靠近下环出水边出现负压区，同一水头下，开度越大，最大负压值越大；同一开度下，水头越高，最大负压值越大，越容易产生空蚀破坏，也是引起机组产生振动的重要原因之一。计算结果表明，两台机组转轮叶片负压区的面积和最大负压值有差异，导致两台机的空蚀破坏和振动程度不相同。

(3)由压力脉动的频谱分析可知，蜗壳里压力脉动主频为140.0 Hz，与水轮机和转轮叶片固有频率相同，主要受转轮叶片的动静干扰的影响；尾水管里压力脉动的影响因素较多，除了和转轮叶片频率影响外，还包括水轮机固有频率10.0 Hz及其倍频、水轮机导叶频率160.0 Hz等其他因素的影响。影响最大的还是尾水管涡带掠过的区域，其最大功率谱密度值是蜗壳里的18倍，表明尾水管蜗带单位时间所含振动能量较大，这是造成机组振动的主要因素。

(4)在设计开度下，没有尾水涡带出现。随着活动导叶开度减小, 尾水管涡带由转轮泄水锥底部产生，并逐渐发展，在25°开度下发展为相对于尾水管中心线对称的三根涡带，并且涡带强度也逐渐增大。在活动导叶开度从25°减小到15°的过程中，尾水涡带又逐渐过渡到一根涡带，形成偏心涡带，涡带强度进一步增大。模拟结果与实际情况吻合。

(5)由于两台机组之间尾水管进口振动因子存在差异化，说明两台机组之间振动原因的不同，振动结果不同，导致1号、2号机组的振动情况不同；同一台机组，在相同开度下，随着水头增加，尾水管处涡带减小，振动减小，但在相同工况下，两台联合运行机组的尾水管进口处的涡带差别却不大。说明引起两台机组振动的差异的位置主要发生在转轮处而非尾水管的涡带，但振动的主要因素却是由尾水管涡带引发的。

(6)尽管采用Real Realizablek-ε模型比较成功地模拟了一管多机式混流式水轮机的涡带和压力脉动特性，但由于该模型本身的时均性，两机振动的差异并不明显。在本文的基础上，有必要进一步对模型进行研究和对比，以及开展流-固耦合的多物理场的研究，以便能更好地模拟出机组之间振动的差异和规律，为机组稳定运行和维护提供科学依据。

□

参考文献：

[1]陈婧,马震岳,刘志明,等.水轮机压力脉动诱发厂房振动分析[J].水力发电,2004,30(5):24-27.

[2]肖若富,韦彩新,韩风琴,等.混流式水轮机转轮的动力学研究[J].大电机技术,2001:41-43.

[3]王正伟,周凌九,黄源芳.尾水管涡带引起的不稳定流动计算与分析[J].清华大学学报 (自然科学版), 2002,42(12):1 647-1 650.

[4]刘宇,杨建明,戴江,等.混流式水轮机三维非定常湍流计算[J].水力发电学报,2004,23(4):102-105.

[5]ZUO Zhi-gang, LIU Shu-hong, LIU De-min et al.Numerical analyses of pressure fluctuations induced by interblade vortices in a model Francis turbine[J].Journal of Hydrodynamics, 2015,27(4):513-521.

[6]XIAO Ye-xiang, WANG Zheng-wei, ZHANG Jin, et al.Numerical predictions of pressure pulses in a Francis pump turbine with misaligned guide vanes[J].Journal of Hydrodynamics, 2014,26(2):250-256.

[7]QIAN Zhong-dong, YANG Jian-dong, HUAI Wen-xin.Numerical simulation and analysis of pressure pulsation in Francis Hydraulic Turbine with air admission[J].Journal of Hydrodynamics, 2007,19(4):467-472.

[8]桂中,唐澍,潘罗平.混流式水轮机尾水管非定常流动模拟及不规则压力脉动预测[J].中国水利水电科学研究院学报, 2006, 4(1):68-73.

[9]杨建明.水轮机尾水管和转轮中湍流计算研究[D].北京：清华大学,1999.

[10]李文锋,冯建军,罗兴锜,等.基于动网格技术的混流式水轮机转轮内部瞬态流动数值模拟[J].水力发电学报,2015,34(7):64-73.

[11]苏华山,陈从平,赵美云,等.泵轮轴向振动条件下高速液力耦合器特性[J].2017,33(7):51-57.

[12]曾立飞,刘观伟,毛靖儒,等.调节阀振动对阀内流场影响的数值模拟[J].中国电机工程学报, 2015,35(8):1 977-1 982.

[13]庞立军,吕桂萍,钟苏,等.水轮机固定导叶的涡街模拟与振动分析[J].机械工程学报,2011,4(22):159-166.

[14]陈世云.大桥电厂水轮机过流部件汽蚀分析及处理[J].四川水力发电,2006,25(3):88-90.

[15]唐巍,张昌兵,李雷.中低比转速混流式水轮机空化与水力振动分析[J].水电能源科学,2014,32(8):145-147.

[16]李雷,张昌兵,唐巍.混流式水轮机水力振动及补气减振研究[J].中国农村水利水电,2015,(2):132-135.

[17]Shih T H, Lion W W, Shabbir A, et al.A newk-εeddy viscosity model for high Reynolds number turbulent flow[J].Computer Fluids, 1995,24(3):227-238.

[18]张师帅.CFD技术原理与应用[M].武汉:华中科技大学出版社,2016:84-85.

[19]王福军.计算流体力学分析——CFD软件原理与应用[M].北京:清华大学出版社,2004:114-141.

[20]符杰,曾永忠,宋文武,等.贯流混流式水轮机尾水管压力脉动分析[J].设计与研究,2005,6(32):5-6,24.