CFD-Rotordynamics Sequential Coupling Simulation Approach for the Flow-Induced Vibration of Rotor System in Centrifugal Pump

Featured Application

The numerical approach in this article can provide abundant information about the flow induced vibration of the pump rotor and may offer theoretical guidance to the design of reducing vibration and safety monitoring of centrifugal pumps.

Abstract

Vibration of the rotor system is closely related with the operation stability of centrifugal pump, and it is inevitably induced by the unsteady inner flow. An unsteady computational fluid dynamics model coupling with a rotordynamics model was presented, and the corresponding numerical calculation program, including a self-designed rotordynamics code, was developed on the commercial package ANSYS Workbench. The validity of the numerical calculation model was verified by a hydraulic performance and vibration test based on an industrial centrifugal pump. The hydraulic radial forces on impeller, pressure pulsation, deformation, and vibration of the main shaft under nine different flow rates were systematically investigated and explained preliminarily from the view of inner unstable flow. Results show that the blade passing frequency is the dominant frequency in the fluctuation of the above dynamical behaviors, which is closely related to the rotor–stator interaction between the impeller and volute casing. This study built the connection between internal fluid flow in the centrifugal pump and the vibration of its external rotor structure, and may provide theoretical references for the design of vibration-reducing and safety monitoring strategies design of centrifugal pumps.

1. Introduction

As an indispensable important energy conversion device, centrifugal pump is extensively applied in various agricultural and industrial fields. For the reduction of equipment investment and energy consumption, modern centrifugal pumps have developed towards the direction of higher flow capacity, rotation speed, head, and efficiency [1,2,3,4,5]. This means that the vibration problem of centrifugal pump tends to become more prominent due to the increased energy density, and it will inevitably bring about unfavorable effects on the operational safety, stability, and life. The hazard of the vibration problem is specially reflected in the rotor assembly of pump, and it may lead to the serious damage of high-speed running parts such as the shaft, impeller, and bearings [6].

There are mainly two kinds of factors influencing the vibration of centrifugal pump rotor assembly—the mechanical factors and the hydraulic factors. The mechanical factors generally relate to the unbalance mass of the rotating parts, misalignment of shaft system and radial rub between different parts, so those factors can be effectively weakened by improving the machining accuracy and assembly quality [7,8]. While for the hydraulic factors, they are predominated by the complicated fluid–structure interaction (FSI) due to inner wall pressure fluctuation caused by the unsteady flow in centrifugal pump. As a result, it is very necessary to obtain the induction mechanism of the rotor assembly vibration from the view of internal unsteady flow, so the vibration level can be restrained at the source via inner flow field optimization.

The unsteady flow behaviors in centrifugal pump has been systematically investigated in the way of numerical calculation and experiment test. The intense rotor–stator interaction due to the relative movement between impeller and volute forms the primary character of internal unsteady flow, so the dominant frequencies of unsteady pressure fluctuation are usually the blade passing frequency (f_BPF) and multiple shaft frequency (f_n) [9]. Similar pulsation tendency is also fund in the hydraulic radial force on the impeller, because the hydraulic force mainly depends on the integration of pressure on the impeller surfaces. In addition, there are also other obvious frequency components rather than f_BPF or multiple f_n, which are stimulated by a serious complex unsteady flow behaviors, such as the jet-wake pattern along the blade extension direction, the flow impact on the tongue, and the flow separation on the blade suction surface [10,11,12]. The above behaviors usually become more evident under off-design conditions and would take much greater risks for the operation of centrifugal pump [13,14].

Based on the research of internal unsteady flow in centrifugal pump, the structural dynamic response analysis has been introduced to perform the FSI study. Because of the complex structure and twisted shape of pump, currently the typical application for the FSI problem of pump is the loosely coupled (explicit) simulation between computational fluid dynamics (CFD) and computational structural dynamics (CSD) [15]. Considering the fluid-induced vibration of the pump brings very limited displacement of the fluid–solid interface, and their interaction acts as a sustained effect in the whole operation process, the FSI problem can be solved with the partitioned method, in which the equations governing the internal flow and structural deformation can be solved separately and expediently with independent field solvers. The principal advantage of the partitioned method is that the whole modeling and solving procedure for the FSI problem can be built on existing and mature CFD and CSD codes with the widely used arbitrary Eulerian–Lagrangian description for the materials, where the Eulerian (spatial) method is used to descript the fluid domain and the Largangian method for the solid component. In this way, not only the required computational resources for both CFD and CSD solvers can be largely saved, but the information exchange (interface mesh deformation) and loads transfer (pressure from fluid domain and displacement on the solid surfaces) representing the FSI effect can be applied in an efficiency and stable way. In recent years, coupling simulation based on the partitioned method has been widely applied in the FSI problems of pump and other turbomachineries, either with the strong two-way coupling strategy or the one-way sequential coupling strategy, and most of the application cases achieved satisfactory prediction results in the deformation, oscillation, and orbit curves of the impeller [16,17,18,19,20].

However, present reported numerical methods are not able to meet the requirements for the prediction of rotor vibration in centrifugal pump. The main reason is that the solving of CSD in current work is built on the general dynamic equation, thus neither the rotary inertia nor gyroscopic effects can be considered, and the stiffness and damping characteristic of bearing has been ignored. In such a case, the dynamics strain and stress distribution of the impeller can be obtained for use of strength checking and safety evaluation, but the motion features about the shaft system cannot be acquired with acceptable accuracy. In fact, the flow-induced vibration behaviors of the shaft system contain a lot of useful information, and they can be easily measured in a cost-effective manner [21]. Therefore, the monitoring and analyzing of shaft vibration data has become the mainstream approach for the safety assessment and fault diagnosis of fluid machinery [22,23].

Based on the previous work by Zhang et al. [1,24,25], this research presents a new one-way coupling simulation approach for the flow-induced vibration of rotor assembly in centrifugal pump, with special attention on the association features of unsteady internal flow and shaft vibration. Compared with former works, the first and foremost improvement of this work is the application of rotordynamics in the CSD, where the bearing damping matrix and gyroscopic matrix have been added to the general dynamic equation. With this advantage, the above new approach is able to provide a wealth shaft oscillation information including vibration velocity and shaft centerline orbit, with both the unsteady hydraulic force on the impeller and the rotor dynamic properties in consideration.

In this article, an industrial split-case double suction centrifugal pump was chosen for study, and the unsteady flow as well as rotor assembly vibration behaviors under different operation conditions were investigated. The remaining of this article was organized as follows. In Section 2, the mathematical model and its overall frame are presented. Section 3 describes the computation details and experiment setup for a research case. Section 4 and Section 5 are the analysis and discussion for the numerical and experimental results, followed by the conclusions in Section 6.

2. Sequential Coupling Model

2.1. Hypothesis and Simplification

This study constructs the CFD-rotordynamics sequential coupling numerical model of centrifugal pump mainly based on following hypotheses and simplifications:

(1) The coupling between flow field and rotordynamics belongs to one-way coupling. It means that only influences of flow induced exciting force on the vibration behavior of the rotor system were considered, but interferences of the deformation and whirl orbits of impeller to the flow field are ignored. According to the calculation and experimental results of a centrifugal pump made by Pei et al. [14,20], compared with the two-way coupling method with consideration of impeller deformation, the calculated radial force based on one-way coupling shows little differences in numerical values. In addition, the vibration displacements of impeller predicted by both one-way coupling and two-way coupling show very small errors with experimental results. However, compared with the two-way coupling method, the one-way coupling way can save as much as 75% computing time and achieves stronger convergence ability during numerical solving.

(2) In flow field analysis, internal leakage of fluid media through sealing rings was neglected.

(3) In rotordynamics analysis, mass distribution was hypothesized uniform on the impeller, namely the unbalanced mass force of the rotating impeller was treated to be zero. In addition, hole and slot structures of the main shaft are ignored. Therefore, the rotor system is considered as a rotationally symmetric structure composed of some mass points of impellers, a main shaft, and its bearings.

2.2. Overall Model Framework

The overall framework of the unsteady flow field-rotordynamics coupling model of centrifugal pump is shown in Figure 1. The coupling model is mainly composed of two parts. One is the CFD model based on finite volume method and the other is the rotordynamics model based on finite element method. The CFD model implements the steady-state CFD solving by using a multiple reference frame (RMF) model in the rotating impeller region, and its computation results are used as initial conditions for the following unsteady CFD solving, where the time varying hydraulic excitation forces acting on the impeller can be obtained. Then the exciting forces on three orthogonal spatial directions are transported to the rotordynamics model though load time series; thus, the sequential coupling process between those two models can be completed.

2.3. CFD Model

2.3.1. Governing Equations

A governing equation group of fluid mechanics in a centrifugal pump is set up according to Navier–Stokes equation and continuity equation

$$\frac{\partial \rho }{\partial t}+\nabla ·(\rho \overrightarrow{u})=0$$

(1)

$$\frac{\partial \rho }{\partial t}(\rho \overrightarrow{u})+\nabla ·(\rho \overrightarrow{u}\overrightarrow{u})=-\nabla P+\nabla ·[\mu (\nabla \overrightarrow{u}+{(\nabla \overrightarrow{u})}^T)]+\rho \overrightarrow{g}+\overrightarrow{f}$$

(2)

where t represents time (s), $\overrightarrow{u}$ is the velocity vector (m·s⁻¹), p is static pressure (Pa), ρ is density (kg·m⁻³), $\overrightarrow{f}$ is external body force term (N·m⁻³), µ is the effective viscosity (kg·m⁻¹·s⁻¹) and g is the gravity constant (m·s⁻²).

2.3.2. Turbulence Model

The governing equations are closed by introducing the SST k-ω turbulence model [26]. The SST k-ω model integrates advantages of both k-ε turbulence model and k-ω turbulence model, and it controls value of dissipation source term in a specific dissipation equation with the blending function F₁. As shown in Equation (3), coefficient φ can be changed with the variation of F₁ function

$$\phi =F_1{\phi }_1+(1-F_1){\phi }_2$$

(3)

where $\phi$, ${\phi }_1$, and ${\phi }_2$ are coefficients of the SST k-ω turbulence model, k-ω turbulence model and k-ε turbulence model, respectively. When F₁ = 0, turbulence in the region far away from the wall is solved by the k-ε turbulence model. On contrary, the k-ω turbulence model under low Reynolds number is applied in the region close to the wall when F₁ = 1. As a result, the SST k-ω turbulence model are believed to be able to provide a relatively high prediction accuracy of the performance curves of pump with acceptable computational cost [27,28].

2.3.3. Treatment of the Gird of Impeller Zone

A grid technique is needed to treat the transition and connection issues between the rotating impeller and other non-rotating flow zones of centrifugal pump. There are two different grid techniques in this work. In the steady CFD model, the MRF model The MRF model was introduced in the steady CFD model, which applies a stationary reference frame in the suction and discharge chambers, and a rotating reference frame for the impeller zone. In this way, the governing equations for fluid zones with different reference fames are solved independently, and local velocity and pressure information on their interfaces is exchanged by coordinate transformation processing. While in the unsteady CFD model, the sliding meshes technique is used in the rotating impeller zone, and a uniform stationary reference frame is set for the whole fluid domain, which can actually reflect more detailed flow behaviors caused by the rotor–stator interaction (RSI) in pump.

The sliding meshes technique is realized by modifying the velocity term of the convection diffusion reaction (CDR) equation for any fluid control volume within the computation domain based on its movement feature relative to the stationary reference frame

$$\frac{d}{dt}{\displaystyle \int }\rho \varphi dV+{\displaystyle \int }\rho \varphi (\overrightarrow{u}-\overrightarrow{u_g})·d\overrightarrow{A}={\displaystyle \int }\psi \nabla \varphi ·d\overrightarrow{A}+{\displaystyle \int }{S}_{\varphi }dV$$

(4)

where Φ represents a general flow scalar for an arbitrary control volume V (m³), $\overrightarrow{A}$ is the area vector of the control volume (m²), $\overrightarrow{u_g}$ is the mesh velocity vector of the moving mesh (m·s⁻¹), ψ and S_ψ are the diffusion coefficient and source term of scalar Φ respectively.

For a centrifugal pump, no volume overlapping occurs between dynamic and stationary meshing regions during the rotating of the impeller zone, which is attributed to the cylinder interface of dynamic and stationary meshing regions. Therefore, volume of the control volume in Equation (4) keeps constant. Hence, it can be simplified as

$$\frac{d}{dt}{\displaystyle \int }\rho \varphi dV=\frac{{(\rho \varphi V)}_{i+1}-{(\rho \varphi V)}_i}{\Delta t}$$

(5)

where i represents the respective value at the ith time level, and i+1 denotes its next one.

2.4. Rotordynamics Model

Based on the Jeffcott rotor model [29], a rotordynamics model of centrifugal pump can be constructed. It can be seen from Figure 2 that a typical rotordynamics model is composed of a principal axis, an impeller and two bearings.

The rotordynamics model is built based on the rotordynamics equation, which is transformed from the universal dynamics equation by adding the gyroscopic matrix and rotating damping matrix. It is expressed as

$$\left[M\right]\{\ddot{U}\}+(\left[C\right]+\left[G\right])\{\dot{U}\}+(\left[K\right]+\left[B\right])\left\{U\right\}=\left\{F\right\}$$

(6)

where is {U} is the displacement matrix; [M], [C] and [K] are the mass matrix, damping matrix and stiffness matrix, respectively; [G] is the gyroscopic matrix for describing the rotating inertia of the rotor system, and its value is determined by rotating speed; [B] is the rotating damping matrix closely connected to the rotor stability, and its value is also related with rotating speed; {F} is the external force vector, which mainly includes the hydraulic force caused by the inner unsteady flow and the gravity of the impeller itself in this study. The Equation (6) is applicable to calculate rotor movement in the same stationary reference coordinate system.

3. Numerical Calculation Application and Experimental Setup

3.1. Case Examined

In this study, commercial computational codes Fluent and ANSYS Mechanical were applied as the CFD and rotordynamics analysis tools, respectively. The whole process including pre-processing, post-processing, physical field coupling, and numerical solving were implemented on the ANSYS Workbench co-simulation platform [30].

A horizontal split centrifugal pump with a six-blade double suction impeller was chosen as the research object. The pump was manufactured by Hunan M&W pump Co., Ltd., and it was designed for the clean water transportation. Its whole major structure is completely symmetric, so the axial hydraulic force on the impeller is fully balanced in theory due to the double suction structures of the suction chamber and impeller. The main structural and performance parameters of the centrifugal pump are listed in

Table 1. Basic parameters of the research object

Parameter	Symbol	Value
Rated flow	Qdes	620 m3/h
Rated head	Hdes	96 m
Rated speed	ndes	1480 rev min−1
Specific speed	ns	52
NPSH required	NPSHr	3.6 m
Impeller inlet diameter	D1	0.18 m
Impeller outlet diameter	D2	0.56 m
Number of blades	Z	6
Impeller outlet width	b2	0.029 m

. The three-dimensional structure is shown in Figure 3 and the middle section of hydraulic model is represented in Figure 4, where the suction chamber, the impeller with six blade passages and the discharge chamber were exhibited with different colors. In addition, an X-Y orthogonal coordinate system with its origin located at the impeller center was established, and the rotating of impeller was set as counterclockwise.

Figure 5 shows the three-dimensional structure and parameters of the main shaft, where the installation positions of bearing were marked with upward arrows. The main shaft is made of structural steel and can be divided into five sections with different diameters. Specifically, the middle section with a diameter of 0.049 m and a length of 0.664 m is designed for carrying the impeller and its accessories. Two sections next to the impeller (diameter = 0.045 m and length = 0.194 m) are designed for the mounting of ball bearings and the end sections (diameter = 0.042 m and length = 0.135 m) are used for the connecting of coupler.

3.2. Steady State CFD Model Setup

3.2.1. Mesh Generation and Turbulence Model

Hexahedral dominated grids for CFD analysis were generated by the meshing module of ANSYS Workbench. Meanwhile, transition regions like casing tongue were processed by mesh refinement. For near-wall grids, the allowable maximum y⁺ was set as 50 by with the meshing adaptive tool of Fluent after the first CFD analysis (trial computation) [30]. Therefore, the mesh refinement function was initiated and grids in the near-wall region were re-meshed based on the flow distribution results.

In this study, the SST k-ω turbulence model was chosen, whose parameters were listed in

Table 2. Parameter settings of the SST k-ω turbulence model.

Parameter	α1	β	β′	σk1	σω1	α2	β2	σk2	σω2
Value	0.556	0.075	0.09	0.85	0.5	0.44	0.828	1	0.856

.

3.2.2. MRF Model Setup and Boundary Condition

The whole computational domain was divided into three sub-domains, namely, the suction chamber, the impeller and the discharge chamber. These sub-domains exchange the solution information through their interfaces. The MRF model was applied in the rotating impeller region, in which the reference frame was set to be clockwise rotating with a speed of 1480 rev min⁻¹, which equals to the rated rotating speed of the centrifugal pump.

Boundary conditions were set as follows: the entrance boundary was set as a velocity inlet perpendicular to the entrance plane, where the turbulent intensity and viscosity ratio were 5% and 10, respectively. An outflow condition with a 100% backflow ratio was adopted for the outlet to reflect the fully developed flow pattern. The wall condition was used for other surfaces. Specifically, the roughness was set as 2 × 10⁻⁵ m for walls of the impeller and 5 × 10⁻⁵ m for walls in the suction and discharge chambers, which match the actual processing parameters of the pump manufacturer.

3.2.3. Solution Setup

Nine flow rates from Q_0.6 to Q_1.4 at an interval of 0.1Q_des were investigated in this work. The fluid media in pump was set as 20 °C clean water, and its density and viscosity are 998.2 kg m⁻³ and 1.003 × 10⁻³ Pa s, respectively. Second order upwind scheme was applied for discretion of the momentum equation and SIMPLE algorithm was used for the iterative computation. The convergence criterion was less than 10⁻⁴ for the maximum residual values for the mass and momentums in three directions.

The case with Q_des flow rate was chosen to run the grid independence test. Calculated results including head and efficiency under six different mesh densities were presented in Figure 6, whose calculation formulas were expressed as Equations (7) and (8).

$$H=\frac{P_0-P_i}{\rho g}+\Delta {\rm Z}$$

(7)

where H is pump head, P_i and P_o are the total pressure of the inlet and outlet of the pump respectively, and $\Delta {\rm Z}$ represents the height difference between the center points of inlet and outlet faces.

$$\eta =\frac{955\rho gQH}{nT}$$

(8)

where Q is volume flow (m³/h), n is the shaft’s rotational speed (rev·min⁻¹) and T denotes the torque on the main shaft (N m).

Figure 6 shows that both head and efficiency of the pump began to keep stable when the number of girds exceeds 3.70 × 10⁶. And the fluctuation ranges of head and efficiency can be controlled within 0.7% and 0.5%, respectively. Therefore, the grid with number of 3.70 × 10⁶ was chosen for further computation. Numbers of element and node of each sub-domain were listed in

Table 3. Mesh information for calculation.

Fluid Zone	Element Number	Node Number
Suction chamber	974,659	801,071
Discharge chamber	1,261,920	1,261,920
Impeller zone	1,463,735	1,463,738
Total	3,700,314	3,526,729

, and the element quality, skewness factor, and orthogonal quality of the entire gird were 0.79, 0.35, and 0.86, respectively.

3.3. Unsteady State CFD Model Setup

The unsteady CFD model was set basically consistent with the steady-state CFD model except for solving settings and the processing of the impeller rotating regions.

In the unsteady CFD model, the frame of the impeller region was changed into a stationary reference one, and the sliding meshes were applied to grids in the impeller region. In other words, all grids in the impeller region make counterclockwise rotation around the impeller center at the rate of 1480 rev min⁻¹.

With respect to solving, the calculated results of steady-state CFD model were used as the initial conditions for unsteady CFD model under the same flow rate to assure convergence performance. In this study, the time step for solution was 1.1261 × 10⁻⁴ s, which equals to the time for the impeller to rotate by 1°. And the whole solving time is 0.2432 s, corresponding to six revolutions of the impeller. The convergence standard adopted by each solving step was that the maximum residual values reached 10⁻⁴ of the initial values of mass and momentum, but the maximum iteration number for each time step was controlled to be less than 20. Besides, various monitoring settings were added during the unsteady CFD solving process, aiming to get time series data of calculated parameters such as mean pressure and flow rate at entrance and exit of the pump, torque on impeller and hydraulic forces. The computation was carried out on the HP workstation (CPU: E5-2695V2, Memory: 48GB), and it took about 6 days to complete one case.

3.4. Rotordynamics Model Setup

3.4.1. Resonance Check

Before rotordynamics analysis, modal analysis on the rotor system of the studied centrifugal pump was carried out to eliminate possibility of structural resonance under the rated rotating speed. The modal analysis shows that the first six orders of inherent frequencies of the rotor system were 67.04 Hz, 118.61 Hz, 138.60 Hz, 235.26 Hz, 253.01 Hz, and 661.44 Hz, respectively. Since none of these were close to the rotating frequency of main shaft or its harmonics, there is no risk of resonance.

3.4.2. Element Type Selection

Finite element modeling and calculation of rotordynamics model were completed by ANSYS Parametric Design Language (APDL) language programming [30]. The codes of the whole program were provided in the Appendix A of this article. During the programming, different types of finite elements were applied for different structural parts.

Beam 188 element with a three-dimensional linear beam and two nodes was chosen for the main shaft, which is suitable for the analysis the problem with large rotation angle and large structure deformation. Mass 21 element was assigned to the impeller, and information of the mass and rotating inertia of the mass point can be added through real constants. COMBI 214 element was used for the bearings, which is a kind of spring-damper element with two nodes.

3.4.3. Mesh and Material Attribution

As shown in Figure 7, the main shaft meshed into high-quality structural hexahedral grids. The minimum grid size was set 0.002 m and each radial section was assured to have at least 20 layers of elements in order to assure solving accuracy. Local grids at ends of shaft were zoomed up and positions of the bearing were marked with ‘K0’.

The main shaft was made of structural steel, whose physical properties were listed in

Table 4. Material properties of the main shaft.

Density (kg m−3)	Elasticity Modulus (Pa)	Poisson’s Ratio
7890	2.09 × 1011	0.269

. According to the data provided by the manufacturer, the rigidity of the rolling bearings was set 5 × 10⁹ N/m, and the damping can be ignored.

3.4.4. Constraint and Load Condition

Nodes on the bearings which are far away from the main shaft (representing the outer rings of the bearings) were fixed. Loading conditions of the rotordynamics model include two parts: (1) forces changing with time was applied on the mass point of impeller, including hydraulic radial forces on the impeller along the X and Y directions (coming from the calculated results of unsteady CFD model) and the gravity of impeller along the Y direction. (2) Rotational angular velocity was applied onto the main shaft and impeller to calculate the Coriolis force. The Coriolis force comes from the rotating inertia of the rotor system, and it is used to describe the migration effect of rectilinear motion mass in the rotating system. The calculation formula of the Coriolis force is shown as

$$F_k=-2m\overrightarrow{{\omega }_r}\times \overrightarrow{v^{\prime }}$$

(9)

where $F_k$ is the geostrophic force (N), m denotes the mass (kg), $\overrightarrow{{\omega }_r}$ is the angular velocity vector of rotor (rad/s), and $\overrightarrow{v^{\prime }}$ is the velocity vector relative to the reference rotating frame (m/s).

3.4.5. Solution Setup and Damping Coefficient Selection

Unsteady state solving was also applied in the rotordynamics finite element analysis model, whose time step and solving time keep consistent with those in the unsteady CFD analysis. Transient complete method was activated in the computation process, namely the unsteady responses of the structural to external loads were calculated by complete mass, damping, and stiffness matrixes. In this way, various nonlinear features such as elasticity, large strain, and large deformation of main shaft can be obtained more accurately.

To assure solving stability during dynamics analysis of structures, a damping matrix is usually needed to consider influences of damping on structural movements in practical engineering. According to the viscous damping theory, the rotor damping is proportional to the vibration velocity and its acting direction is opposite to the velocity. Therefore, a proportional damping was applied onto the entire finite element analysis model. This implies that the system damping matrix [C] in Equation (6) is approximately treated as the linear combination of mass matrix [M] and the stiffness matrix [K], as expressed in the Equation (10)

$$\left[C\right]=\alpha ·\left[M\right]+\beta ·\left[K\right]$$

(10)

where α and β denote the mass and stiffness damping coefficients respectively.

With references to research results of Pei et al. [14,20], a stiffness damping coefficient was added into the rotordynamics model and its value was determined by trial-and-error method and engineering experiences. A trial operation was implemented under the flow rate of Q_des, in which a relatively small stiffness damping coefficient (β = 0.001) was set firstly. Figure 8 shows that the velocity curves of mass point of impeller was very unstable with β = 0.001, and no valid information can be obtained from them. Then, when the value of stiffness damping coefficient was increased to β = 0.002, the velocity curves began to be stable and appear periodic fluctuation feature after oscillation in the first revolution of impellers (corresponding to 0.041 s). Thirdly, when the stiffness damping coefficient β reached 0.003, it takes a shorter time for regulation compared with that under β = 0.002. However, it takes more system input disturbance under β = 0.003. Therefore, the stiffness damping coefficient in the computation was finally determined as β = 0.002.

3.5. Experiment Setup

3.5.1. Hydraulic Performance Test

To verify the numerical results, a hydraulic performance test of the industrial centrifugal pump was carried out in an industrial test rig, during which vibration signals were collected at the bearing end cover. The test system was shown in Figure 9, which agree with the level-1 standards in ISO 9906-2012 [31]. The test rig mainly includes two valves, two connected water tanks, a pump, an electromotor and some measuring instruments. During test, the water temperature was controlled within 19.8°C ± 0.3 °C. The random uncertainty in the process of data acquisition was within ±0.15%. The flow rate was measured by electromagnetic flowmeter with a measuring error of ±0.2%. Two pressure gauges were utilized to obtain the pump head, and the measuring error of each one was ±0.1%. A torque transducer and a laser tachometer were applied to measure the torque and rotation speed of the main shaft respectively, and their measurement errors were both ±0.1%. Thus, based on uncertainty analysis, the total measurement uncertainties for the head and efficiency are 0.21% and 0.38% respectively.

3.5.2. Vibration Measurement

Vibration behaviors of the bearing end cover of the centrifugal pump was collected during the test, which can be used to approximately represent the vibration level of the rotor system according to Chinese standards of GB/T 29531-2013 [32] and GB/T 6075.1-2012 [33].

It can be seen from Figure 10 three sensors with a common base were mounted at the center of the bearing end cove. With that test instrument, vibration accelerations on three orthogonal spatial directions (two radial directions and one axial direction) were measured with PTZ piezoelectric ceramic sensors, and the sample frequency was set as 10.24 kHz. The frequency response and measurement ranges of the sensors are 0.5 Hz–5 kHz and 0–50 g, respectively. The INV3018A data collection instrument with 24-bit resolution and 4 simultaneous-capture analog input channels was utilized for analog acquisition, digital conversion, signal filtering, and data storage, and the acquisition system accuracy was 0.02%.

4. Results Analysis

4.1. Numerical Model Validation

4.1.1. Comparison of Performance Curves

The tested and calculated performance curves of the centrifugal pump were shown in Figure 11, where the unsteady calculation results were the mean values gathered from the unsteady CFD calculation for the last impeller revolution.

In general, the predicted performance curves agree well with test results, especially in the nearby range of Q_des. For the flow of Q_des, the maximum errors between calculated head coefficients C_h and its test value are less than 0.036, with a relative deviation of about 3.3%, and the absolute difference between calculated efficiency η and its test value is controlled within 0.5%. When the flow moves far away from Q_des, the error between calculated value and test value increases with the increase of flow offset, especially in the part-load conditions. However, the maximum relative errors of C_h and η are only about 4.6% and 3.9% (around Q_0.6), respectively.

The computational errors are mainly attributed to the following two aspects. Firstly, both leakage and mechanical loss were ignored in the CFD model. Then, errors may come from the discretion and solving processes due to the complicated flow behaviors in centrifugal pump. Specifically, the solving accuracy of the SST k-ω turbulence model declines in part-load conditions, thus resulting in the higher solving error. Additionally, it can be seen from Figure 11 that the calculation curves of steady-state CFD analysis and unsteady CFD analysis differ to some extent, and their differences are much more significant under part-load conditions. Generally, the unsteady CFD results provide higher prediction accuracy. A possible explanation is that the inner flow is more unstable under part-load condition, and the unsteady CFD model with sliding meshes technique can reflect the RSI better than the steady CFD model with MRF setting for the impeller grids.

4.1.2. Comparison of Vibration Characteristics of Main Shaft

The variation curves of vibration accelerations along X and Y directions versus time were gained from numerical calculation and vibration measurement. In the numerical model, a node of the main shaft which connects with the bearing was chosen for study and comparison. It was found that spectral distribution characters of vibration acceleration under different flow rates are similar in some aspects, so the Q_des condition was chosen for the comparisons of the spectral distribution of the vibration accelerations, as shown in Figure 12.

In Figure 12, distribution laws can be seen between computation and measurement results. Specifically, the vibration acceleration magnitude reaches the peak at f_BPF (148 Hz), which is the dominant frequency. Meanwhile, the magnitudes at f_n (24.7 Hz), f_BPF and their harmonics. Moreover, the acceleration magnitude along X direction is generally higher than that along the Y direction. To sum up, the numerical calculation of vibration information of the rotor system is reliable.

However, the calculated results are different from measurement results to some extent with respect to specific vibration acceleration magnitude and distribution features. On one hand, the calculated acceleration magnitudes under dominant frequency is higher than the measurement results. For example, the calculated vibration acceleration magnitudes along X and Y directions under the f_BPF are 0.88 m s⁻² and 0.26 m s⁻² respectively, which are higher than the corresponding measurement results (0.46 m s⁻² and 0.19 m s⁻²). On the other hand, the measured vibration acceleration magnitudes at some other multiple shaft frequencies (e.g., 98.7 Hz and 296 Hz) are significantly higher than the calculated results.

Differences between calculated results and measurements are mainly caused by the following three aspects. Firstly, the unsteady CFD model applies Reynolds-average method to average the instantaneous values of turbulence phenomenon, thus ignoring the effects of high-frequency random disturbance. Secondly, vibration of the rotor system in the numerical model is only attributed to hydraulic exciting. However, due to mechanical processing and assembly, it is inevitable to have some other vibration sources in practical operation, such as external forces caused by mass imbalance of impeller, rotor asymmetry, and electrical motor. Thirdly, the computation model cannot provide vibration information of bearing directly due to the limitation of element types in the rotordynamics finite element model. Instead, bearing vibration is approximately reflected by information on the shaft node connected with the bearing. Considering the high stiffness of the bearing, the vibration acceleration measurements are smaller than the calculated values, and there are multiple harmonic components in measurement results.

In this study, a quantitative evaluation on vibration energy of a centrifugal pump was carried out by using the root-mean-square (RMS) value of vibration velocity. At the first step, curves of vibration velocity versus time were gained through time integral of vibration acceleration series. Then, RMS of vibration velocities under different flow rates were calculated and shown in Figure 13. It is easy to see that the RMS velocity distribution feature of computation results is highly similar to measurements. Generally speaking, vibration velocity level is relatively low around Q_des, and the RMS ranges are between 0.4 m s⁻¹–1.0 m s⁻¹. And RMS velocity along the X direction is generally higher than that along the Y direction. Nevertheless, RMS velocity increases significantly under par-load and large flow conditions, and it increases to the peak under Q_0.6. Moreover, the calculated RMS velocity along two directions are slightly higher than measurements, which is similar with the patterns in Figure 12.

4.2. Fluctuation of the Hydraulic Radial Force on the Impeller

In unsteady CFD analysis, variations of hydraulic radial forces along the X and Y directions of impeller of the centrifugal pump with time under different flow rates in the last revolution of impeller were shown in Figure 14a,b, respectively. Radial forces along X and Y directions present sine fluctuating characteristics with f_BPF as the dominant frequency. Under Q_0.6, the hydraulic radial force along the X direction has a high value and it points to the negative direction of the X-axis. When the flow rate increases to Q_0.7, the radial force along the X direction decreases to some extent, but it still points to the negative direction of X-axis. While under Q_0.8, the absolute value of radial force along the X direction further declines and its positive and negative direction alter with the rotation of impeller. Then under large flow conditions, the hydraulic radial force along the X direction increases gradually with the continuous increase of flow rate and it keeps positive in X-axis. To sum up, the hydraulic radial force along Y direction points to the negative direction under different flow rates and its value is positively correlated with flow rate.

4.3. Pressure Pulsation Characteristics

Fluctuation of hydraulic radial force on the impeller is closely related with pressure fluctuation of the internal unsteady flow field. On the mid-section of discharge chamber, five monitoring nodes (P1–P5) close to the volute wall were chosen to investigate the pressure fluctuation laws. The distance between each node and the volute wall is 0.02 m. As seen in Figure 15, P1 locates near the volute tongue and the rest four points are on the X-axis and Y-axis, respectively. Flow rate Q_des was taken as an example, under which the pressure fluctuations at five monitoring points in one impeller revolution were shown in Figure 16. Figure 16 shows the dominate pressure fluctuation period is 1/6 of the rotating period of the impeller, which agrees with the experimental regularity obtained by Wang et al. [12]. However, pressure values and their amplitude variations are different for each monitoring point. Similar pressure fluctuation features were observed under other flow rates.

To further compare pressure fluctuation characteristics under different flow rates, the time-varying curves $P(node,t)$ of pressure at different monitoring nodes are decomposed into the fixed component $\overline{P}(node)$ and the fluctuating component $\tilde{P}(node,t)$, as shown in Equation (11).

$$P(node,t)=\overline{P}(node)+\tilde{P}(node,t)$$

(11)

where the numerical value of $\overline{P}(node)$ is defined the arithmetic mean value of pressure in one revolution

$$\overline{P}(node)=\frac{1}{N}\sum _{i=0}^{N-1}P(node,t_0+i\Delta t)$$

(12)

where N is the number of time steps, t₀ is the initial moment and $\Delta t$ is the time step length.

Thus the value of $\tilde{P}(node,t)$ is the difference between instantaneous pressure value $P(node,t)$ and $\overline{P}(node)$

$$\tilde{P}(node,t)=P(node,t)-\overline{P}(node)$$

(13)

$\tilde{P}(node,t)$ is worth in-depth study as it indicates the instability of the inner flow of centrifugal pump. Hence, the dimensionless quantity namely the pressure fluctuation intensity coefficient $C_P^{\prime }$ is defined as the root-mean-square of $\tilde{P}\left(node,t\right)$ divided by reference pressure P_ref

$$C_P^{\prime }=\frac{{\tilde{P}}_{rms}(node)}{P_{ref}}$$

(14)

where ${\tilde{P}}_{rms}(node)$ and P_ref are shown in Equations (15) and (16), respectively.

$${\tilde{P}}_{rms}(node)=\sqrt{\frac{1}{N}\sum _{i=0}^{N-1}\tilde{P}(node,t_0+i\nabla t)}$$

(15)

$$P_{ref}=\frac{1}{2}{\rho u}_2^2$$

(16)

where u₂ is the circumferential speed at the impeller outlet.

Figure 17 shows $C_P^{\prime }$ of the five monitoring nodes under different flow conditions. $C_P^{\prime }$ at P1 is higher than that at other nodes, and it presents a V-shaped variation trend with the increase of flow rate, where it reaches the maximum value (about 0.02) under Q_0.6. $C_P^{\prime }$ decreases to the minimum (slightly higher than 0.01) under Q_0.9. Subsequently, it increases continuously with the continuous increase of flow rate, but it is still less than that in part-load conditions. $C_P^{\prime }$ values at P2–P5 are significantly lower than that at P1, and they are generally higher under part-load conditions.

4.4. Deformation and Vibration Behaviors of the Main Shaft

4.4.1. Deformation and Stress Distribution

This work focuses on the mechanics behaviors of the rotor system caused by the inner unsteady flow. At first, the displacement and equivalent stress distribution of the main shaft at the end moment of unsteady rotordynamics analysis under Q_des were shown in Figure 18 and Figure 19, respectively.

Figure 18 shows that the deformation occurs mainly in radial directions. The maximum displacements in X and Y directions reach 0.127 mm and 0.194 mm, while it is only about 0.02 mm in Z direction (axial direction). Furthermore, deformation of the main shaft is symmetric around the XY plane located on the origin, where the deformation reaches maximum. On the contrary, the minimum deformation is at the bearing node due to its constraint function.

Figure 19 the equivalent stress reaches maximum (22.3 MPa) is at the center of the main shaft, and the strength safety coefficient is about 15.2. The equivalent stress decreases gradually from the center of the main shaft to its end.

The time series of equivalent stress at the center of main shaft under different flow rates were gained and analyzed. Similarly to the data processing of pressure fluctuation curves in Section 4.3, the time-varying equivalent stress was decomposed into mean and fluctuation components, as shown in Figure 20.

In Figure 20, the mean and fluctuation components of equivalent stress reach their minimums at Q_0.8 and Q_0.9, respectively. Both of them increase significantly with flow rate under large flow conditions. Stress and its fluctuation level are important factors for the fatigue failure of main shaft [34], which is one of the most common failure modes of pumps. Therefore, the equivalent stress information predicted by the proposed model in this work is able to provide useful quantitative reference information for the fatigue strength analysis of centrifugal pumps.

4.4.2. Vibration Behaviors

The central point of main shaft was chosen for the vibration behavior analysis in this work, where the impeller is mounted. Its time-varying acceleration and velocity curves under different flow conditions are shown in Figure 21 and Figure 22, respectively. In addition, the corresponding frequency distributions were shown in Figure 23 and Figure 24.

It can be seen from Figure 23 and Figure 24 that for both acceleration and velocity, the dominant frequency of the shaft vibration is f_BPF (148 Hz), and its harmonics (296 Hz, 444 Hz, …) are also obvious. Besides, harmonic component of both the vibration acceleration and velocity is more obvious under part-load conditions, so there are some distinct small peaks for the corresponding time domain curves in Figure 21 and Figure 22. Nevertheless, harmonic component is weakened gradually with the increase of flow rate and there are little harmonic components under large flow rates, where the time-domain sine curves are relatively smooth.

In the engineering practice, the peak and RMS values in a certain period of time are counted for the vibration acceleration and velocity signals, respectively. The acceleration peak indicates the external force induced the vibration, and the RMS velocity reflects the strength of vibration energy. Figure 25 compares the above two vibration indicators in the center of the main shaft under different flow conditions.

It can be seen from Figure 25 that acceleration peak and velocity RMS are relatively low around Q_des. Given large flow (Q_1.3 or higher) or part-load conditions (Q_0.7 or lower), both acceleration peak and velocity RMS increase significantly with the increase of deviation from Q_des, and this phenomena is more obvious under part-load conditions.

4.5. Evolution Characters of the Unsteady Flow Field

4.5.1. Flow Fields under Q_des

It is the unstable flow in the centrifugal pump that induces the various radial force fluctuations, pressure pulsations, and vibration behaviors presented in the preceding sections. Inner flow field observation reveals there are obvious flow distribution changes in the discharge chamber under different flow rates and different impeller phases, while it keeps fairly uniform and stable in the suction chamber and impeller zone. Thus, the internal flow analysis mainly concentrates on the discharge chamber in this work.

Firstly, the inner flow distribution under Q_des was analyzed and used as the basis for the follow-up comparative study of flow fields under different flow rates. In this study, the rotation phase of impeller was defined zero when the impeller center, tip of aspecific blade, and tip of volute tongue meet on a straight line, as illustrated in Figure 4.

The velocity vector distributions on middle section of the discharge chamber with impeller phase of 0° and 30° were shown in Figure 26a,b, respectively. Obviously, fluid media sheds from outlet of blade passages along the tangential direction and then enters the discharge chamber with relatively regular distribution excepet for gentle impact near the tongue and the consequent secondary flow in the diffuser (marked with dotted frame).

Figure 27 shows the turbulence kinematic energy distribution on the middle plane of the impeller zone and discharge chamber. It is easy to find that there are not any obvious differences between the 0° and 30° impeller phases for both velocity and turbulence energy distribution. In Figure 28, high turbulence energy areas mainly locate at the blade wake region and the diffuser region nearby the tongue, where there are typical unstable flow structures such as jet-wake, cutting, and impact. Under Q_des, the maximum values of turbulence energy under 0° and 30° of impeller phases are close, and they are 5.31 J kg⁻¹ and 6.54 J kg⁻¹, respectively.

4.5.2. Flow Fields under Q_0.6 and Q_1.4

Due to space limitations, the minimum and maximum flow rates (Q_0.6 and Q_1.4) in the calculation range were chosen as the representative to analyze the evolution laws of flow field. The velocity vector and turbulence kinematic energy distributions under these two operation conditions were shown from Figure 28, Figure 29, Figure 30 and Figure 31.

Figure 28 indicates strong flow impingement and cutting near the tongue, followed by highly disordered secondary flow and stall flow phenomenons (marked by dotted frame) in the diffuser. In response, Figure 30 shows the turbulence kinematic energy is high in large scale, especially for the clearance between blade tips and volute wall, the blade passages near pressure surfaces, and laryngeal. Furthermore, the flow field differences between the 0° and 30° impeller phases are significant both for the velocity and turbulence kinematic energy distributions.

While under Q_1.4, as shown in Figure 29 and Figure 31, there is no visible differences between those two blade phases. The velocity vector distributions are stable and uniform in the impeller zone and tip clearance. Although there is large scale stall flow with high turbulence kinematic energy in the diffuser, but it is in the downstream and far away from impeller zone, so its impact on the impeller is very limited.

5. Discussion

The above series analysis proves that hydraulic-induced vibration and deformation of the main shaft are closely related with the unsteady flow in centrifugal pump. Around the rated flow rate, the overall flow field is uniform and the secondary flow phenomenon is controlled well. Consequently, the pressure pulsation and turbulence kinematic energy are both weak, which ensures the main shaft work in stable operation with relatively low vibration and stress levels. Given the part-load conditions, the RSI is strong when the impeller blade sweeps over the volute tongue, and it brings a series of unstable flow phenomena, such as the flow cutting and impingement in the nearby region of tongue, which has been observed by PIV measuring [10]. The inner unstable flow finally causes high fluctuation of radial forces of the impeller and significant vibration of the main shaft. While in the large flow conditions, although the radial forces on the impeller is high, their fluctuation degree is controlled effectively and the main shaft vibration level is lower than that in part-load conditions. One possible explanation is that the high flow suppresses the separation degree of fluid media, and it is the flow rate factor rather than pump structure that dominant the entire flow process. However, illustrating the flow induced mechanism of the rotor vibration is a very complex task, and a more accurate CFD computation, such as detached eddy simulation [10] and some new comparative analysis methods, may be needed in the further research.

Compared with the two-way FSI study of the centrifugal pump made by Benra [16] and Pei [14,20], this work pays more attention on the vibration behaviors of the main shaft rather than the impeller, considering the commonly occurring shaft fracture failures. Furthermore, the Coriolis force and bearing characteristics have been included in the model of this work, which were ignored in previous research [14,16,20].

6. Conclusions

The present work emphasizes on the development of a sequential coupling mathematical model for unsteady flow field-rotordynamics of the centrifugal pump as well as the corresponding numerical calculation program. A horizontal split centrifugal pump with double-entry impeller was chosen for research, and the validity and accuracy of the computation were verified through its performance test and vibration measurement. The coupling model presented in this work provided a series of interconnected computation information including the inner flow field and deformation and vibration behaviors of the external rotor system in a flexible high efficiency way, and some conclusions can be drawn as follows:

(1) f_BPF is the dominant frequency in the fluctuation of hydraulic radial force on the impeller, pressure pulsation in the blade tip clearance, and vibration of main shaft, which are closely related to the RSI of the centrifugal pump.

(2) Operating stability indexes like fluctuation of hydraulic radial forces on the impeller, vibration, and deformation levels of the main shaft are predominated by the unsteady flow process in the centrifugal pump. Those indexes are at the optimal level around the rated flow rate, but get deteriorated when it deviates from the best work condition, which is more serious in part-load conditions.

(3) Under large flow rates, although the radial forces on the impeller are larger than other conditions, the vibration levels can be controlled better than those under part-load conditions due to their relatively stable inner flow field.

In sum, this work mainly discusses the rotor vibration behaviors induces by inner flow process of centrifugal pump. A more precise unsteady CFD model should be conducted in future work to obtain more accurate time-varying laws of the hydraulic exciting forces. Moreover, more measuring points, including other non-intrusive sensors, needed to be supplemented in the experimental validation work.

7. Patents

There has been a Chinese invention patent resulted from this work: Zhang, H.H.; Zhou, T.H. Vibration Simulation Analysis Method and Device of Vane Pump Rotor System. China Patent ZL2016110994972, 2016.