水泵水轮机启动过程转轮动应力分析

桂中华 夏 翔 王 薇 周凌九 刘殿海 王正伟

(1.国网新源控股有限公司抽水蓄能技术经济研究院， 北京 100761; 2.中国农业大学水利与土木工程学院， 北京 100083;3.清华大学能源与动力工程系， 北京 100084)

0 引言

在运行过程(特别是过渡过程)中，水轮机转轮承受剧烈的水压力脉动，可能产生严重振动和高水平动应力，甚至可能产生共振和疲劳[1-3]。与常规水轮机组相比，水泵水轮机长期处于启动、停机和负荷变化等过渡过程中，面临的稳定性问题更加突出。准确预测水泵水轮机转轮在过渡过程(如启动过程)中的动应力水平具有重要意义。

试验测量是一种有效分析转轮动应力的手段[3-4]。但是，试验方法的成本和难度很高，并且无法应用在设计阶段。随着计算机技术的进步，数值模拟为转轮动应力的研究提供了便利。水轮机转轮在水压力作用下产生动力响应是典型的流固耦合(Fluid-structure interaction, FSI)问题。而流固耦合问题的求解方法可分为强耦合法和弱耦合法。其中弱耦合(即分离式求解)方法的实施相对容易，因此得到了广泛的应用。一些学者采用分离式求解的双向流固耦合方法对水轮机和水泵的内部流动及结构响应进行了研究[5-9]。上述工作都对研究对象作了一定简化，因为对于复杂结构而言，双向耦合方法需要消耗大量的计算资源，在实际工程中并不适用。事实上，在水力机械正常运行时，转轮的弹性变形量远小于流体流动的特征长度，忽略结构振动对流体流动的影响可满足精度要求。因此，也有学者采用分离式求解的单向流固耦合方法进行研究[10-14]。然而，上述单向流固耦合瞬态计算方法都没有考虑水体的附加质量，理论上对于动应力幅值特别是共振工况动应力的预测存在一定误差。文献[15]基于声固耦合提出了改进的单向流固耦合(Acoustic-structure based one-way FSI, ASOW-FSI)方法，并以此预测了模型水泵水轮机转轮在不同转速工况下的动应力，结果表明计算得到的共振曲线与试验数据高度吻合。文献[16]借助圆柱涡激振动研究了ASOW-FSI方法在流固耦合问题中的适用性，结果表明这一方法不能预测频率锁定等自激振动现象，但对于受迫振动计算具有较高的精度。更多关于流固耦合问题的研究现状可参见文献[17-19]。

对于复杂的水力机械转轮，动应力的获取需要消耗大量的计算资源。而一个完整的过渡过程通常长达几十秒，通过数值计算准确预测其动应力特性十分困难。因此，目前关于瞬态过程转轮动应力数值计算的报道较少。本文采用基于声固耦合构建的单向流固耦合方法研究某原型混流式水泵水轮机转轮在启动过程中的振动响应，重点分析水力激励力的瞬态特性对动应力的影响。

1 数值模型

1.1 动态响应分析方法

线性系统的结构动力学方程可以表示为

(1)

式中Ms——结构系统的质量矩阵

Cs——结构系统的阻尼矩阵

Ks——结构系统的刚度矩阵

u(t)——节点位移向量

F(t)——外部激励力

在流固耦合问题中，F(t)即为水压力，可通过计算流体动力学(CFD)分析得到。

在进行双向流固耦合计算时，流体流动产生的非恒定水压力使得结构产生振动。同时，结构的运动也会反作用于流场。而在单向耦合中，忽略了结构运动对流场带来的影响。事实上这种简化不仅引起流场计算不准确，还会导致对结构动力特性的误判。当结构在流体中变速运动时，将受到额外的与加速度方向相反的流体力。这种现象被称为附加质量效应，现已被广泛研究[20-23]。其中，文献[20]通过试验研究了静水中淹没平板的附加质量效应，结果表明附加质量是结构模态振型的函数，附加质量在低阶模态尤为显著。文献[21]利用试验和数值方法研究了水体附加质量对水泵水轮机转轮动力特性的影响，并从振型分析出发讨论了这种现象产生的机理。文献[22-23]研究了壁面对附加质量的影响，结果表明壁面的接近程度可以极大地改变附加质量，结构离壁面越近，附加质量效应越明显。其中，文献[22]采用声固耦合方法计算了淹没圆盘的模态，并与试验数据进行了比较。结果表明，声固耦合方法可以准确预测考虑壁面影响的水下结构固有模态。

ASOW-FSI方法在动态响应分析中采用声学单元对结构周边的水体进行建模，可以较充分地考虑水体对结构动力特性的影响。此时，系统的结构动力学控制方程可以写为

(2)

式中Ma——附加质量矩阵

Ca——附加阻尼矩阵

Ka——附加刚度矩阵

文献[15]运用这种方法预测了模型水泵水轮机转轮在不同相对转速(即转速与额定转速之比N/Nrated)下的动应力，结果如图1所示。可以看出，ASOW-FSI方法对共振转速和共振幅值的预测都与试验数据吻合较好。与常规的单向耦合(OW-FSI)方法相比，ASOW-FSI方法在共振点附近的计算精度得到了很大的提升。

1.2 流动分析模型

本文研究的混流式水泵水轮机组的基本参数如下：转轮直径Dr=5.1 m，叶片数Zb=9，额定转速Nrated=300 r/min，固定导叶及活动导叶数Zs=Zg=20，导叶最大开度Gmax=36°，额定水头Hrated=295 m，额定流量Qrated=118 m3/s。

为了分析启动过程中水力激励力的瞬态特性对动应力的影响，本文针对某个特定启动过程选取不同转速工况点分别进行稳态和局部瞬态动应力计算。其中，稳态动应力是指假定机组在某工况点稳定运行时的动应力水平；而局部瞬态动应力是指在过渡过程中选取某一时间段进行独立求解得到的结果。机组启动过程中转速和导叶开度的变化规律如图2(图中G表示导叶开度)所示。在转速上升段(图中虚线框所示时间段)选取了4个工况点，对应的运行参数展示在表1(表中H表示运行水头)中。

为了给动应力计算提供相对准确的边界条件，建立了全流道三维非定常流动数值分析模型，其计算域包括蜗壳、固定导叶、活动导叶、转轮流体域、间隙、均压管及尾水管等部分。使用六面体和四面体网格对计算域进行离散，其中，六面体网格主要分布于间隙、均压管、尾水管及活动导叶等区域。计算域网格如图3所示，蜗壳及固定导叶区域的网格单元数8.6×105，活动导叶2.78×105，转轮流体域1.499×106，间隙4.43×105，均压管1.896×106，尾水管8.24×105。

表1 选定工况点的运行参数Tab.1 Operating parameters of selected working conditions

在蜗壳进口设置总压边界，尾水管出口为开放边界，进出口压力均由实测数据换算得到。设置转轮流体域和上、下间隙为转动域，其中靠近转轮室一侧的壁面保持静止，而其它区域均设为静止域，在各计算子域之间设置交界面(包括动-静和静-静交界面)。为了保证交界面上数据传递的准确性，采用CFX提供的GGI网格连接方式进行插值。

选用SSTk-ω湍流模型求解URANS方程。采用有限体积法对控制方程在空间上进行离散，而时间上的离散采用二阶全隐式格式。非定常计算的时间步长Δt设为转动周期的1/200，待测点压力呈现出明显周期性后，至少再计算5个转动周期。

针对该水泵水轮机组在额定工况下稳定运行时的水压力进行了现场测试，重点监测了球阀前后、蜗壳以及尾水管进口处的压力脉动(采样频率为100 Hz)。其中，蜗壳处测点的位置、现场实测得到的压力脉动以及本文CFD计算的结果如图4所示。

可以看出，本文通过CFD计算得到的压力脉动的均值和幅值都与实测数据吻合较好。其中，均值误差为0.08%，幅值误差为7.55%，测点位置的偏离可能是造成压力脉动幅值产生误差的原因之一。本文所述“幅值”均指某一变量时域信号中的峰峰值。由于现场测试时的采样频率较低，其结果无法真实反映压力脉动的频率，因而不能直接用来判定数值模拟对于脉动频率预测的准确性。对数值计算的结果重新进行了采样(采样频率同样为100 Hz)，结果也展示在图4中。结果与实测数据基本吻合，证明了该数值分析模型对压力脉动频率的预测是准确的。

1.3 有限元模型

采用基于声固耦合的单向流固耦合方法计算转轮在运行过程中的动态响应，求解过程在ANSYS APDL中进行。计算域包括转轮结构以及由声学单元模拟的间隙流体和转轮流体域。其中流体域沿用了流动分析中的网格，并且流固耦合交界面两侧的节点一一对应。重点针对叶片及其倒圆处进行了网格加密。有限元模型如图5所示。

作用在转轮结构上的荷载包括非定常水压力、重力及离心力。其中，水压力由流体动力学计算得到，加载的时间步长为转动周期的1/200。在主轴螺栓中心线上施加固定约束，将转轮进出口和间隙进出口的流体域边界设为全吸收边界，而顶盖和座环壁面设为全反射边界。转轮材料为结构钢，密度为7 850 kg/m3，泊松比0.3，弹性模量200 GPa。而水体的密度为1 000 kg/m3，水中声速取1 482 m/s。

2 结果与讨论

2.1 转轮模态

为了解转轮的动力特性，采用声固耦合法计算了转轮在流道中的固有模态。模态分析中的网格、边界条件和材料属性均与动应力分析保持一致。本文主要关注转轮的节径模态(节径是指结构在自由振动过程中位移保持为零的直线)。图6展示了前4阶节径模态的振型，其中1ND、2ND和3ND模态以上冠和下环的振动为主，而4ND模态的最大位移在叶片上。前4阶节径模态的固有频率分别为24.84、57.76、119.04、161.90 Hz。

2.2 动静干涉

由于高水头水泵水轮机转轮的进水边(发电工况)设计较矮，且活动导叶的出水边较厚，在运行过程中转轮叶片和导叶间的干涉作用较强。文献[24]研究表明，高水头水轮机转轮在最佳效率运行点附近时，动静干涉贡献的动应力约占总应力幅值的80%。因此，开展动静干涉理论计算有助于理解和分析非定常水压力以及结构的动态响应。

由转轮叶片和活动导叶联合作用产生的最终压力场可以表示为旋转压力场和静止压力场的累加，经调制，最终压力场在静止参考系中可表示为[25]

(3)

其中

ω=2πN

(4)

式中ω——转轮转动角速度

N——转轮转速

m——转动系统的谐波阶数

n——固定系统的谐波阶数

Amn——联合幅值

θs——静止坐标系中的角位置

φm——转动系统对应阶次谐波的相位

φn——固定系统对应阶次谐波的相位

这个压力场具有两种径向压力模态，节径数k分别为

k1=mZb-nZg

(5)

k2=mZb+nZg

(6)

一般情况下，谐波阶次越高，振幅越小，因此可以只考虑k=k1的情况。其在旋转坐标系和静止坐标系下角速度和激励频率分别为[1]

ωr=nZgω/k

(7)

ωs=mZbω/k

(8)

fr=nZgfN

(9)

fs=mZbfN

(10)

式中ωr、ωs——旋转坐标系、静止坐标系下角速度

fr、fs——旋转坐标系、静止坐标系下激励频率

fN——转动频率

对于本文研究的对象，动静干涉模式(部分结果)如表2所示。以2节径振型为例：在静止坐标系下，机组各部件受到的动静干涉激励频率为18fN，角速度为-9ω；针对转动部件进行分析时，需将其置于旋转坐标系下来考虑，其2ND水力激振力的频率为20fN，角速度为-10ω。当转轮结构的2ND模态频率等于20fN时，可能引发共振。

2.3 稳态动应力

首先针对前文选定的4个工况点进行稳态计算。在转轮振动和动应力研究过程中，需分析其周边的水压力特性，因此设置了如图7所示的7个压力脉动测点。以OP1工况为例，不同测点处的压力脉动时程曲线如图8所示。由于p2和p3处的压力脉动特性基本一致，在分析中略去p3的结果。

表2 动静干涉模式(部分结果)Tab.2 Modes of rotor-stator interaction

从图8可以看出，OP1工况下无叶区和叶片进水边处的压力脉动最显著，其幅值分别为0.12 MPa和0.10 MPa。此外，叶片出水边处也出现了明显的压力脉动，约为0.07 MPa。而间隙内的压力脉动相对较弱，幅值约为0.04 MPa。受篇幅限制，本文没有将OP2～OP4工况测点处的压力-时间曲线列出。但是图9展示了不同工况下各测点处的压力脉动频域信号。各工况下压力脉动的频率特征基本一致：无叶区的压力脉动频率为叶片通过频率9fN及其倍频；转轮叶片进水边及间隙处为导叶通过频率20fN；出水边的频率成分相对复杂，有20fN、18fN以及由漩涡运动引起的压力脉动频率(如OP1中的2.33fN和OP4中的0.667fN)。其中，0.667fN对应尾水管涡带的转动频率0.333fN，由于压力测点处于旋转坐标系下，监测到的脉动频率与涡带转频相差fN。

OP1至OP4工况下的转频fN分别为2.87、3.2、3.8、4.5 Hz，对应的20fN分别为57.4、64、76、90 Hz。其中，OP1工况的20fN约等于转轮结构在流道中的2ND模态频率57.76 Hz，结合动静干涉理论分析可知，此时转轮可能会发生共振。

完成CFD计算后，将转轮表面的水压力转换为节点荷载，即可进行结构响应分析。图10为OP1工况下某一时刻转轮的等效应力云图。可以看出，叶片与上冠、下环相接处出现了明显的应力集中，其中叶片进水边靠近上冠处最为明显，最大等效应力为317 MPa。在后续分析过程中，分别在叶片进水边靠近上冠处和靠近下环处设置动应力监测点pc和pb。

用测点处每一时刻的等效应力减去时均值，即可得到动应力时程曲线。不同工况下pc和pb处动应力如图11所示。可以看到，在OP1工况，转轮产生了明显的共振，pc和pb处的动应力幅值分别达到了47.9 MPa和47.1 MPa，并且还有进一步增大的趋势。此时，转轮结构的2ND模态被激发，振动频率为20fN。而在转速较大的OP2、OP3、OP4工况，转轮的振动较弱，pc处的动应力幅值分别为7.5、6.2、6.6 MPa。其中OP3和OP4工况下动应力的频率成分相对简单，主频为20fN。而在OP2工况，转轮的振动形式较为复杂，上冠和下环处的动应力主频分别为18fN和20fN。

事实上，OP2工况的转频fN为3.2 Hz，而静止坐标系下机组的动静干涉频率18fN=57.6 Hz，恰好对应转轮结构的2ND模态频率，但是其激励振型与结构模态振型存在一定偏差。因此，转轮的固有频率在结构响应中表现突出，但不会引发共振。根据动静干涉分析结果可知，在旋转坐标系下，节径数为2的水力激振力频率为20fN，转速为-10N，每隔时间Δt=1/(20fN)激励一次，并且在此过程中激振力旋转了半个周期，而转轮结构在旋转坐标系下保持静止，因此恰好可以激发出结构的2ND振型；同理，在静止坐标系下，节径数为2的激振力频率为18fN，转速为-9N，每经过时间Δt=1/(18fN)，激振力旋转半个周期，但是由于转轮本身的转动，每一次激励都会产生20°的偏差。从另一个角度来说，即18fN对应的激振力在旋转坐标系下并不是严格对应2ND振型。上述分析可以用图12进行描述。

2.4 瞬态动应力

实际上在启动过程中，转轮所受的水力激励力是非周期性的。总体上，随着转速逐渐增大，激励力的脉动频率越来越高。针对上述4个工况点进行局部瞬态计算，分析激励力的瞬变特性对转轮动应力的影响。首先，转速在OP1工况(即共振点)附近时危险点pc处的瞬态应力计算结果如图13所示。从应力水平(即时均应力)的角度来看，随着转速不断增大，等效应力水平逐渐上升。从动应力的角度看，当转速增至OP1工况时，并没有出现明显异常的应力波动。但是随着时间推移，应力的波动逐渐增强，达到峰值后逐渐减弱。最大应力波动出现的时间较共振点滞后约0.25 s。这一现象是由转轮结构的惯性引起的。另外，应力波动的幅值约为24 MPa，比稳态计算的结果(47.9 MPa，见图11a)小得多。这是因为结构的振动需要一定时间才能激励出来，而在机组启动过程中激振力的频率一直在增大，与转轮固有频率重合的时间较短，所以并未产生非常强烈的振动。由此可以推断，启动过程中转速变化得越慢，在共振点附近能够激励出的振动越强烈。上述结果表明，在共振点附近采用稳态方法估算过渡过程中的动应力水平是不可行的。

接下来，以OP4工况为例，分析非共振区域的瞬态计算结果(如图14a所示)。同样，可以看出随着转速不断增大，应力水平逐渐上升。实际上，这是由于转轮离心力和流场中的压力梯度逐渐增大所致。为了剔除应力水平的变化对动应力幅值评估的干扰，假定应力水平与转速在局部呈线性关系，然后采用最小二乘法进行线性回归分析，结果如图中的虚线所示。由此，可以得到转速在OP4工况附近时pc处的动应力，如图14b所示。结果表明，动应力幅值为6.7 MPa，与稳态计算的结果(6.6 MPa，见图11d)非常接近。说明在非共振区域可以采用稳态方法估算过渡过程中的动应力水平。

根据上述分析可知，可以采取稳态和局部瞬态相结合的方法来研究过渡过程中结构的动应力。具体的实施过程如下：首先，进行非稳态流动分析；第2步，进行流固耦合结构模态计算；第3步，找到共振工况点以及由于导叶开度、转速等运行参数瞬变引起高水平压力脉动的工况点；第4步，针对上述工况点进行局部瞬态动应力计算，其它区域可选择性地取点进行稳态分析。

3 结论

(1)采用基于声固耦合构建的单向流固耦合方法研究了某原型混流式水泵水轮机组在水轮机模式启动过程中转轮的振动响应和动应力。结果表明，在非共振区域，激振力的瞬变特性对转轮的振动几乎没有影响，具体表现为稳态计算和瞬态计算得到的动应力水平基本一致；而在共振区域，激振力的瞬变特性对转轮的振动有显著影响，一方面会使启动过程中的动应力幅值远比采取稳态方法估算的值低，另一方面，还会使得高水平振动滞后于共振点。上述结果可为过渡过程转轮动应力的评估提供参考。

(2)借助模态分析和动静干涉理论对转轮的共振特性进行了分析。研究表明，仅当某一水力激励力与转轮固有模态的频率和振型均对应时才会产生共振；当频率相等而振型之间存在相对转动时，这种激振力会在一定程度上主导转轮的振动，但振动幅值只是较常规情况略微增大。