余热排出泵叶轮流固耦合特性分析

摘 要：余热排出泵长期在高温高压环境下运行，其结构的可靠性对整个机组的安全运行有着重要的影响。该文采用单向流固耦合方法研究了几种常用材料及增加叶轮盖板厚度对转子结构动力学特性的影响，同时对比了常温和高温下叶轮的振动特性。结果表明，余热排出泵高温运行时，叶轮的第1阶固有频率为394.17 Hz，较常温下升高了2.28%，远高于叶轮叶片通过频率。在各阶模态振型下，1Cr13MoS叶轮固有频率最大，ZG225-450叶轮固有频率最小。设计工况和偏工况下叶轮的应力及变形变化趋势基本一致，叶轮的变形随半径的增加而不断增大，最大变形量出现在后盖板叶轮出口处。叶轮的应力分布不均，最大应力均出现在叶片尾缘与后盖板接合处。沿着前盖板和叶片接合线，应力的峰值随前盖板厚度增加而减小。沿着后盖板和叶片接合线，增加后盖板厚度明显减小了应力峰值。高温和常温应力变化趋势基本一致，但是高温下叶轮应力明显高于常温。尤其在叶轮进口处附近，高温应力相对于常温应力增加超过300%。该研究结果为下一步进行更加复杂的动力学分析、疲劳分析以及结构优化提供参考。

关键词：泵；叶轮；动力学；余热排出泵；盖板厚度

0 引 言

随着中国对核电的高度重视，国内高校和泵企业等正加紧研制开发核电用泵，实现核电用泵的国产化。在核电站内的泵类设备中，余热排出泵是重要的设备之一，可保护核电站的安全运行：当反应堆进行正常停堆时，排出堆芯余热和核主泵在一回路中产生的热量；当发生严重事故时，减小核电站受损程度[1-2]。由于余热排出泵高温高压运行环境的特殊性，其结构的可靠性对整个机组的安全运行有着重要的影响，因此，准确地对余热排出泵转子部件进行强度分析可为其结构设计及寿命延长提供依据。

近年来，流固耦合有限元分析已成为分析流体机械可靠性的重要方法。国内外不少学者尝试采用流固耦合的方法对转子系统进行振动特性分析，取得了许多有价值的成果。刘厚林等[3]分别运用单、双向流固耦合方法对导叶式离心泵叶轮进行静应力强度分析，指出单向耦合即可满足叶轮静力学分析。Kan等[4]采用单向流固耦合计算得到了水轮机转子在多工况下应力和变形变化规律。牟介刚等[5]采用单向流固耦合方法，分析了悬臂式离心泵在不同悬臂比、不同工况下叶轮应力和变形的变化规律，并对不同悬臂比离心泵的转子动力学特性进行研究。Schneider等[6]对多级离心泵进行了单向流固耦合计算，分析了叶轮结构设计参数和温度对叶轮应力和变形的影响。刘厚林等[7]对比分析了无预应力下、流固耦合作用力和旋转离心力作用下的余热排出泵转子的模态。刘建瑞等[8]考虑了水的附加质量对余热排出泵转子振动性能的影响，对转子进行水中模态分析。Zhou等[9-10]采用双向流固耦合计算方法，研究了考虑流固耦合作用的余热排出泵在不同工况下内部流场及转子应力和变形分布。考虑到余热排出泵实际高温高压的运行环境，目前很少有学者考虑对叶轮盖板进行加厚处理。同时，材料属性对转子动力学特性的影响研究较少，尤其是对于有特殊用途的余热排出泵而言，合理选择叶轮材料至关重要。

因此，本文应用单向流固耦合计算方法研究了余热排出泵实际高温运行环境下叶轮模态分布，详细分析了叶轮几种常用材料和不同盖板厚度下叶轮变形量和应力分布，以期为研究整个系统的动态特性提供参考。

1 流固耦合动力学方程

1.1 基于模态分析的静流体整体耦合法

对于具有n个自由度的系统，其结构动力方程为

 （1）

式中Ms是系统的质量矩阵；Cs是系统的阻尼矩阵；Ks是系统的刚度矩阵；Fs是系统承受的载荷，Pa；是节点的加速度矢量，m/s2；是节点的速度矢量，m/s；是节点的位移矢量，m[11]。

在静流体中，考虑到流固耦合作用对方程做了一些假设：流体是可压缩的，压力波动将会引起密度变化；流体是非黏性的（即无黏性耗散）以及非旋转的；流体没有运动。因此，流体将被视为声学流体，控制方程可用如下的声波方程表示

 （2）

 （3）

式中c是声音在水中的传播速度，m/s；p是由声波引起的流体压力以及作用于流体的力，Pa；Δ是拉普拉斯算子；k是流体体积模量；t是时间，s；ρ0是流体平均密度，kg/m3[11]。

在流固耦合交界面，用如下的方程描述了流体与结构的相互作用

 （4）

式中是交界面处单位法向量；是沿着法线向量的压力梯度[11]。因此，可将方程（2）离散化为

 （5）

式中，，和分别表示流体等效质量矩阵、流体等效阻尼矩阵、流体等效刚度矩阵以及流体压力载荷。此外，作用于流体的等效耦合质量矩阵为，其中为在交界面处由结构运动产生的流体荷载[11]。

考虑到在交界面处压力对结构的影响，方程（1）改写为

 （6）

其中，Ksf = -R表示等效耦合刚度[11]。

因此，可以耦合流体和结构方程，基于单一求解器通过整体求解方法同时求解[11]，方程表达如下

（7）

在对叶轮进行模态分析时，自由振动方程系统视为非阻尼情况[11]，方程（7）简化如下

 （8）

1.2 基于应力和变形分析的动流体分区耦合法

在分区耦合法中，有2种耦合求解流体动力学和结构响应策略：单向耦合和双向耦合[12]。双向耦合方法主要是针对具有强物理耦合效应的问题，即除了考虑流动对结构的影响，结构的变形或运动对流动的反作用也需要考虑[13]。在泵机械中，大多数情况下由于泵内不稳定流引起的叶轮的变形很小，而叶轮变形和振动对流场影响可以忽略不计[10,14]。因此，本文采用单向耦合方法，如图1所示，在每一个时间步长流体载荷通过交界面单向传递到结构。

注：ti表示第i步的时间，s；ti+1表示第i+1步的时间，s；Δt表示时间步长，s；p(t)表示流场压力，Pa；v(t)表示流场速度，m·s-1；F(t)表示流固耦合面的流体力，Pa。

Note: ti is time at the i step, s; ti+1 is time at the i+1 step, s; Δt is time step, s; p(t) is flow field pressure, Pa; v(t) is flow field velocity, m·s-1; F(t) is fluid force along fluid-structure interfaces, Pa.

图1 单向耦合求解流程图

Fig.1 Flowchart of one-way coupled solution

2 模型基本参数及边界条件

2.1 计算模型

余热排出泵过流部件包括后掠扭曲叶轮、径向导叶和环形蜗壳3部分，其三维结构如图2所示。余热排出泵的主要性能和几何参数分别为设计流量Qd=910 m3/h，扬程Hd=77 m，转速n=1 490 r/min，比转数ns=104.5，叶轮进口直径Dj=270 mm，叶轮出口直径D2=511 mm，叶轮叶片数Zi=5，导叶进口直径D3=515 mm，导叶出口直径D4=718 mm，导叶叶片数Zg=7，蜗壳基圆直径D5=840 mm。为探讨增加盖板厚度对转子结构动力学特性的影响，设计了4种不同的叶轮结构：1）原始盖板；2）前盖板厚度增加50%；3）后盖板厚度增加50%；4）前后盖板厚度均增加50%，详细参数见表1。余热排出泵长期在高温高压环境恶劣环境下运行，叶轮材料的选取至关重要。因此本文选取5种常用材料[15]如表2所示，分析了不同材料叶轮的动力学特性。

注：a表示前盖板厚度，b表示后盖板厚度，mm。

Note: a is front shroud thickness, b is rear shroud thickness, mm.

图2 泵及叶轮模型

Fig.2 Pump and impeller model

2.2 网格划分

2.2.1 流体域网格

网格划分是数值模拟中的重要步骤之一，网格的好坏直接影响数值模拟的准确性。对余热排出泵计算域采用ICEM软件进行结构化网格划分，对壁面进行网格加密[16-17]，如图3所示。对数值模拟进行网格无关性检查，如图4所示，确定余热排出泵计算域的最终网格总数为5.20×106。

表1 叶轮结构参数

Table 1 Geometrical parameters of impeller

表2 叶轮材料参数

Table 2 Material properties of impeller

图3 计算域网格

Fig.3 Mesh of computational domain

图4 网格无关性分析

Fig.4 Mesh independent analysis

2.2.2 结构域网格

为研究湿模态下叶轮的动力学特性，将叶轮完全淹没在1.12 m×0.36 m×1.12 m的水体中[18]。如图5所示，运用ANSYS Workbench中自带的网格划分软件，对叶轮及水体部分进行网格划分[19]。分别基于Method和Hex Dominant网格类型通过Sizing命令对叶轮和水体域进行网格划分。经网格无关性检查，确定最终网格节点总数为129 220，网格单元总数为49 556。

图5 数据插值映射过程

Fig.5 Data interpolation mapping procedure

2.3 边界条件

采用整体耦合法对叶轮进行模态分析时，为了获得更高的求解精度，用Fluid220和221高阶单元作为水体域单元，默认使用Solid186和187作为结构域单元。由于水体域单元不支持完全积分法，所以本文对水体域积分控制使用缩减积分法，而对流固耦合面的固体域积分控制使用完全积分法。考虑到湿模态分析时，构成结构动力学有限元方程的质量矩阵是非对称的，因此在求解中需选用非对称模态算法[20]。

假设流体在泵内部是不可压缩三维湍流流动，满足连续方程和动量方程，采用SST k-ω湍流模型求解雷诺时均方程[21-24]。采用分区耦合法对叶轮进行应力和变形量分析时，先对流体域进行了三维定常和非定常数值模拟：在定常计算过程中，边界条件设置为总压进口和质量流量出口，动静域之间交界面采用Frozen rotor，计算域壁面采用无滑移网格函数[25-26]，采用高阶求解精度，收敛残差RMS设置为10-5，计算迭代步数最大500；在非定常数值计算过程中，以定常计算结果作为初始值，边界条件设置不变，而动静域之间交界面改成Transient frozen rotor[27]。计算周期为6圈，叶轮旋转1度为1个时间步长，即5.747 13×10-5 s。以非定常计算最后一步的计算结果作为初值对叶轮进行了6圈单向流固耦合计算：时间步长设置和非定常计算一致；流体场和固体场叶轮叶片、前后盖板设置为流固耦合交界面（fluid-structure interface）；有限元时间积分采用采用Hilber-Hughes- Taylor（HTT）离散方法，使用瑞利阻尼（Rayleigh damping）系数来定义结构瞬态阻尼效应[28]。如图1所示，设置叶轮后盖板的孔面为Fixed Support支撑面。本文分别模拟了设计工况下25 ℃常温水以及180 ℃高温水余热排出泵内部流动以对比实际运行环境下叶轮的动力学特性。

2.4 试验验证

为了验证数值模拟的准确性，对余热排出泵按比例缩小30%的模型泵进行了外特性试验验证。试验在江苏大学流体机械质量技术检验中心的开式试验台上进行，试验现场如图6所示。采用上海自仪九仪表有限公司生产的型号规格为LWGY-200A涡轮流量计进行测量流量，量程为0～1.6 MPa，精确度为0.5级；采用上海威尔泰仪器仪表有限公司生产的WT200智能压力变送器测量进出口压力，进出口处压力变送器的测量量程分别为−0.1～0.1 MPa和0～1.6 MPa，测量精度均为0.1级。如图7所示，根据相似换算原理[29-30]，得到了真实尺寸泵的试验和数值模拟性能曲线。由图7可知，试验扬程和效率略高于计算值，预测扬程的最大偏差为6.78%，预测效率的最大偏差为9%。但是试验结果趋势一致，计算模型较准确地预测了余热排出泵的外特性，数值计算结果较为可信。

图6 开式试验台

Fig.6 Open test rig

图7 试验结果与模拟结果对比

Fig.7 Comparison of experimental and numerical results

3 结果及分析

3.1 模态分析

为了了解叶轮模态振型情况，首先对原始盖板的合金结构钢叶轮的计算结果提取分析。如图8所示，叶轮第1、2阶振型相似，都是在叶轮某2个对称位置的变形量最大，即叶轮是摆动变形。第3、4阶模态振型相似，都是叶轮沿着轴向的扭动变形。第5阶和第6阶模态振型相似，都是扭曲变形，可以看到2条明显的节径，且分布角度不同。第7、8、9、10阶振型为叶轮前后盖板沿着轴向的扭曲变形。通过对其他几种形式的叶轮前10阶模态振型分析发现，叶轮振型基本相同，但是固有频率有一定差别。

表3列出了不同材料、盖板厚度和温度下叶轮前10阶固有频率。由表3可知，改变叶轮盖板厚度对前6阶固有频率的影响可以忽略。从第7阶开始，叶轮后盖板的厚度对固有频率有很大影响，增加后盖板厚度可以提高叶轮的固有频率。5种材料的叶轮固有频率变化规律趋于一致。在常温各阶模态振型下，1Cr13MoS(B_3)叶轮固有频率最大，ZG225-450(B_2)叶轮固有频率最小。以合金结构钢(B_0)原始盖板叶轮为例，余热排出泵在高温运行时，叶轮的第1阶固有频率为394.17 Hz，较常温下升高了2.28%，远高于高于叶轮叶片通过频率，在安全运行范围内。

注：上面一排视图是叶轮正视图，下面一排视图是叶轮右视图。

Note: Top row figures are front view of impeller, and bottom row figures are right view of impeller.

图8 原始叶轮前10阶模态振型图

Fig.8 First ten-order mode of vibration on original impeller

表3 叶轮前10阶固有频率

Table 3 First-ten order natural frequency of impeller

3.2 应力和变形分析

如图9、10所示分别为叶轮在3种工况下变形及应力分布，其中对变形量和应力进行了无量纲处理[7]，从图中可以看出在设计工况和偏工况下叶轮的应力及变形变化趋势基本一致。叶轮的变形随半径的增大而不断增大，最大总变形出现在后盖板叶轮出口处。叶轮的应力分布不均，最大应力均出现在叶片尾缘与后盖板接合处。

注：上面一排视图是叶轮正视图，下面一排视图是叶轮后视图。设计流量Qd=910 m3·h-1。

Note: Top row figures are front view of impeller, and bottom row figures are back view of impeller. Design flow rate Qd is equal to 910 m3·h-1.

图9 不同工况叶轮变形量分布图

Fig.9 Dimensionless displacement of impeller at different flow rates

注：上面一排视图是叶轮后视图，下面一排视图是叶轮右视图。

Note: Top row figures are back view of impeller, and bottom row figures are right view of impeller.

图10 不同工况叶轮应力分布图

Fig.10 Dimensionless stress of impeller at different flow rates

为了详细分析叶轮关键部位的应力变化，如图11定义了前后盖板和叶片接合处的4条路径，这4条路径在叶轮表面从进口前缘延伸到出口尾缘。

图11 盖板与叶片接合处路径

Fig.11 Investigated paths in junction between shroud and blade

图12为不同工况下叶轮4条路径应力分布图。从图12中可以看出3种工况下4条路径应力变化趋势基本一致：在前后盖板与叶片压力面接合处，路径中间附近出现应力峰值，在前缘附近可以明显的看出大流量工况下应力大，小流量工况下应力小；在前后盖板与叶片吸力面接合处，最大应力均出现在尾缘处。其中，叶片吸力面与盖板接合处的应力峰值远高于与压力面接合处的应力峰值。

注：横坐标s表示路径上的点距离叶轮进口的无量化大小。其中，横坐标0表示进口处，1.0表示出口处，下同。

Note: Path length s is dimensionless distance away from leading edge of impeller. Path length 0 represents leading edge, while 1.0 represents trailing edge, the same as below.

图12 不同工况路径应力分布图

Fig.12 Distribution of equivalent stress in defined paths for different flow rates

余热排出泵的实际运行环境对叶轮材料的选择提出了苛刻的要求，如图13所示为设计工况下5种叶轮材料的4条路径应力分布图。从图13中可以看出5种材料的4条路径应力变化趋势基本一致，但是在前缘附近合金结构钢(B_0)的应力明显高于其他4种材料。

图13 不同材料路径应力分布图

Fig.13 Distribution of stress in defined paths using different materials

如图14所示为设计工况下4种叶轮盖板厚度对应的4条路径应力分布图。从图14b和14d中可以看出，沿着前盖板和叶片接合线，应力的峰值随前盖板厚度增加而减小。从图14a和14c中可以看出，沿着后盖板和叶片接合线，增加后盖板厚度明显减小了应力峰值。

图14 不同盖板厚度路径应力分布图

Fig.14 Distribution of stress in defined paths with different shroud thicknesses

图15所示对比分析了设计工况下原始模型叶轮在常温和高温下4条路径应力。由图15可知，两种温度下四条路径应力变化趋势基本一致，但是高温下叶轮应力明显高于常温。尤其在叶轮进口处附近，高温应力相对于常温应力增加超过300%。在前盖板和叶片压力面接合处，常温时应力峰值在路径中间附近，而高温时峰值在尾缘。其余3种路径峰值位置没有改变。

图15 不同温度路径应力分布图

Fig.15 Distribution of stress in defined paths under different temperatures

4 结论与讨论

本文应用单向流固耦合计算方法分析了余热排出泵实际高温运行环境下叶轮动力学特性，并定量研究了叶轮不同结构设计参数和不同材料下的变形量和应力，结论如下：

1）余热排出泵高温运行时，叶轮的第1阶固有频率为394.17 Hz，较常温下升高了2.28%，远高于叶轮叶片通过频率。在各阶模态振型下，1Cr13MoS的叶轮固有频率最大，ZG225-450叶轮固有频率最小。改变叶轮盖板厚度对前6阶固有频率影响可以忽略。但是从第7阶开始，增加后盖板厚度可以提高叶轮的固有频率。

2）设计工况和偏工况下叶轮的应力及变形变化趋势基本一致。叶轮的变形随半径的增大而不断增大，最大总变形出现在后盖板叶轮出口处。叶轮的应力分布不均，最大应力均出现在叶片与后盖板接合的进口处。其中，沿着前盖板和叶片接合线，应力的峰值随前盖板厚度增加而减小。沿着后盖板和叶片接合线，增加后盖板厚度明显减小了应力峰值。

3）高温和常温应力变化趋势基本一致，但是高温下叶轮应力明显高于常温。尤其在叶轮进口处附近，高温应力相对于常温应力增加超过300%。在前盖板和叶片压力面接合处，常温时应力峰值在路径中间附近，而高温时峰值在尾缘。其余3种路径峰值位置没有改变。后文将在本文基础上深入分析高温高压条件下泵的转子动力学特性及叶轮使用寿命评估。