Dynamic Characteristic Analysis of the Multi-Stage Centrifugal Pump Rotor System with Uncertain Sliding Bearing Structural Parameters

Abstract

The traditional dynamic characteristic analysis of the multi-stage centrifugal pump rotor system is developed assuming the bearing structural parameters with constant values. However, the manufacturing errors will cause the structural parameters to vary around their nominal values and then affect the dynamic characteristics of the bearing-rotor system. Thus, this paper proposes a method for analyzing the dynamic characteristics of the bearing-rotor system with uncertain bearing structural parameters. First, dynamic characteristic coefficients of the sliding bearing are identified to establish the dynamic model of the rotor system, and its dynamic characteristics are analyzed through finite element simulations. Next, the sliding bearing structural parameters are taken as the variables to establish an optimization model, which is solved by the improved particle swarming optimization algorithm to obtain the extreme critical speed of the rotor system. A case study was carried out on a multi-stage centrifugal pump. The obtained extreme values of the critical speeds were close to those calculated using the multiple samples generated by the Monte Carlo method, indicating that the proposed method can provide accurate variation ranges of critical speeds efficiently and lay a theoretical basis for selecting robust operational speeds and designing the rotor system of the multi-stage centrifugal pump.

1. Introduction

Multi-stage centrifugal pump is a key piece of equipment for the fluid transportation, which has been widely used in various industries such as the metallurgy, chemical industry, electric power and so on [1,2,3]. Rapid developments of these industries push the multi-stage centrifugal pump to increase its speed, resulting in safety and stability failures of the rotor system under operational conditions. Once the working speed approximates the critical speed of the rotor system, a resonance will occur and induce severe vibrations, causing deformations of rotor components and affecting safeties of operators [4]. Therefore, several studies have analyzed the dynamic characteristics of the multi-stage centrifugal pump rotor system. Generally, the dynamic model of the rotor system is constructed to predict the corresponding critical speeds, providing an important approach to select the working speeds and ensure the reliable operations of the rotor system [5,6,7].

Since sliding bearing is an important component to support the rotor system of a multi-stage centrifugal pump, its dynamic characteristics can directly affect the dynamic characteristics of the rotor system [8]. Thus, the prerequisite for establishing an accurate dynamic model of the rotor system is studying the dynamic characteristics of the sliding bearing. Many methods have been provided to identify the dynamic characteristic coefficients of the sliding bearing. The equation proposed by O. Reynolds in 1886 lays a theoretical basis of fluid lubricated bearing research, and solving the Reynolds equation is the key problem of analyzing the hydrodynamic characteristics of the sliding bearing [9]. Researchers at home and abroad have proposed a variety of algorithms to obtain the analytical solution of the Reynolds equation. The numerical analysis and the fluid simulation are the mainly used two approaches. In the numerical analysis, the Reynolds equation is generally solved with the aid of the MATLAB software to obtain the oil film pressure and the dynamic characteristic coefficients of the sliding bearing. Li et al. [10] analyzed the transient response of the rotor system and the pressure distribution of oil film flow field by combining the finite element method and finite difference method; the partial differential equation about the disturbance of dynamic characteristic coefficients of the sliding bearing were solved by the numerical integration. Zhou [11] used the bisection and secant method to iteratively calculate the static equilibrium position of the sliding bearing. On that basis, the small parameter method was adopted to solve the dynamic characteristic coefficients of the sliding bearing. In the fluid simulation, the dynamic model of the bearing fluid domain was established with the aid of the fluid simulation software, and the dynamic characteristic coefficients of the sliding bearing were solved by setting the boundary conditions. Hanoca et al. [12] analyzed the influence of oil film thickness on the dynamic characteristics of the sliding bearing using ANSYS Workbench software. The numerical solution of the computational fluid dynamics was compared with the numerical analysis solution of the Reynolds equation. Manshoor et al. [13] investigated the three-dimensional computational fluid dynamics (CFD) analysis on the performance characteristics of a thin film lubricated journal bearing and solved the characteristics of the sliding bearing using three turbulence models with the ANSYS Fluent software.

The aforementioned studies focused on identifying the dynamic characteristics of the sliding bearing with constant structural parameters, and the disturbances of the structural parameters have not been further discussed. Recently, few researchers have investigated the effects of the manufacturing errors of the sliding bearing structural parameters on the dynamic characteristics of the sliding bearing-rotor system. Wei et al. [14] established the dynamic model of the rotor system and analyzed the influences of the dimensional errors of the dynamic viscosity, bearing width, bearing diameter and journal diameter on the friction power loss. Wang et al. [15] numerically studied the influence of the bearing surface waviness on the performance of the aerostatic sliding bearing, and used the finite difference method to solve the Reynolds equation. The obtained dynamic characteristic coefficients were utilized to analyze the influences of bearing surface waviness amplitude on the bearing capacity, stiffness and friction coefficient. Shin et al. [16] studied the influence of the bearing seal surface roughness on the lubrication characteristics and concluded that the small machining error and surface waviness would affect the stability of the rotor system. Song et al. [17] proposed a new algorithm for dealing with the non-triangular matrix based on the periodic boundary conditions of the journal bearing, and then investigated the half frequency whirling motion and critical speed caused by the shape error of the bearing surface. Later, Zoupas et al. [18] analyzed the influences of the manufacturing errors on the minimum oil film thickness, friction torque and maximum pressure of the large tilting pad thrust bearing.

Although effects of the manufacturing errors on the dynamic characteristics of the sliding bearing have already been investigated, how the manufacturing errors further affect the dynamic characteristics of the sliding bearing-rotor system have not been studied. Moreover, many studies have already found that the parameter disturbances of the rotor system can affect its stability and provided theoretical methods to perform the dynamic analysis of rotor system with uncertainty [19,20]. Therefore, this paper presents a method to study the dynamic characteristics of the multi-stage centrifugal pump rotor system, where the structural parameters of the sliding bearing are uncertain because of the manufacture errors. First, the mathematic relationship between the structural parameters and dynamic characteristic coefficients of the sliding bearing is established according to the Reynolds equation. These structural parameters contain the length, clearance ratio, ratio of length to diameter and oil film viscosity coefficient. On that basis, the dynamic model of the multi-stage centrifugal pump rotor system can thus be established to perform the dynamic characteristic analysis. Next, an optimization model is further established to obtain the extreme value of the critical speed of the bearing-rotor system, where the sliding bearing structural parameters are taken as the independent variables. The finite element simulations and the improved particle swarm optimization (PSO) algorithm are combined to solve the optimization model, and the extreme value of the critical speed can be obtained in the variation ranges of the sliding bearing structural parameters.

Henceforth, this paper is organized as follows. The principles of identifying the dynamic characteristic coefficients of the sliding bearing and constructing the dynamic model of the bearing-rotor system are provided in Section 2. In Section 3, an optimization model for obtaining the extreme value of the critical speed of the rotor system is established and solved by the improved PSO algorithm. In Section 4, a case study for illustrating the proposed method is performed on a real multi-stage centrifugal pump rotor system, and variation ranges of the critical speeds are obtained considering uncertain sliding bearing structural parameters. Finally, conclusions from the current research are summarized in Section 5.

2. Dynamic Modeling of the Bearing-Rotor System of a Multi-Stage Centrifugal Pump

The bearing-rotor system of a multi-stage centrifugal pump is described in Figure 1, which is mainly composed of the impeller, rotating shaft and sliding bearing. As an important component supporting the rotor system, the structural parameters and dynamic characteristics of the sliding bearings can directly affect the dynamic characteristics of the rotor system. Thus, the dynamic modeling of the rotor system considering the dynamic characteristics of the sliding bearing is discussed in this section.

2.1. Calculation of the Dynamic Characteristic Coefficients of the Sliding Bearing

The sliding bearing plays the role of supporting the rotating shaft and reducing the friction during the operation of a multi-stage centrifugal pump [21], a typical structure of which is shown in Figure 2.

According to the bearing fluid lubrication theory, the radial oil film characteristics of the sliding bearing obey the generalized Reynolds equation of fluid mechanics expressed in Equation (1):

$$\frac{1}{R_b^2}\frac{\partial }{\partial \phi }\left(\frac{\rho h^3}{12\eta }\frac{\partial p}{\partial \phi }\right)+\frac{\partial }{\partial z}\left(\frac{\rho h^3}{12\eta }\frac{\partial p}{\partial z}\right)=\frac{\left(v_j+v_b\right)}{2R_b}\frac{\partial \left(\rho h\right)}{\partial \phi }+\frac{\partial \left(\rho h\right)}{\partial t},$$

(1)

where P is the oil film pressure, H is the thickness of oil film, $R_b$ is the bearing radius, η is the oil film viscosity, ρ is the oil film density, φ and z are the circumferential and axial coordinates of the sliding bearing respectively, and v_j and v_b are the linear velocities of the journal and bearing respectively.

Considering the actual working condition, the bearing bush is assumed to be fixed with the assembling unit as v_j =0. Then, Equation (1) can be transformed into the following Equation (2) if the influence of the lubricating fluid medium density is ignored.

$$\frac{1}{R_b^2}\frac{\partial }{\partial \phi }\left(\frac{h^3}{12\mu }\frac{\partial p}{\partial \phi }\right)+\frac{\partial }{\partial z}\left(\frac{h^3}{12\mu }\frac{\partial p}{\partial z}\right)=\frac{{\mathsf{\Omega }}_j}{2}\frac{\partial h}{\partial \phi }+\frac{\partial h}{\partial t},$$

(2)

Subsequently, $\partial h/\partial t=0$ is defined when the sliding bearing is in a stable state, and the following dimensionless parameters are introduced into the equation:

$$\overline{z}=\frac{2z}{L_b},\epsilon =\frac{e}{c},\varphi =\frac{c}{R_b},\overline{p}=\frac{{\varphi }^2}{2\eta {\mathsf{\Omega }}_j}p,\lambda =\frac{L_b}{2R_b},\overline{h}=\frac{h}{c}=1+\epsilon \cos \phi ,$$

(3)

where $\overline{z}$ is the dimensionless axial coordinate, ε is the eccentricity, ϕ is the dimensionless clearance ratio, and λ is the dimensionless ratio of length to diameter. Thus, Equation (2) can be further simplified into a dimensionless form as shown in Equation (4):

$$\frac{\partial }{\partial \phi }\left(\frac{{\overline{h}}^{{}^3}\partial \overline{p}}{\partial \phi }\right)+\frac{\partial }{\partial \overline{z}}\left(\frac{{\overline{h}}^{{}^3}}{{\lambda }^2}\frac{\partial \overline{p}}{\partial \overline{z}}\right)=-3\epsilon \sin \phi ,$$

(4)

A five-point difference method is adopted to discretize the partial differential equation in Equation (4), and then the over-relaxation iterative method is used to calculate the oil film pressure. After the oil film pressure distribution of the sliding bearing is obtained, the dynamic characteristic coefficients of the sliding bearing can be solved by the small parameter method. First, the small perturbations of the displacement and velocity are applied at the equilibrium position, and then the dimensionless thickness and pressure of the oil film under the disturbance are written as:

$$\begin{array}{c}\overline{h}={\overline{h}}_0+\overline{x}\sin \theta +\overline{y}\cos \theta ,\\ \overline{p}={\overline{p}}_0+{\overline{p}}_{\overline{x}}\overline{x}+{\overline{p}}_{\overline{y}}\overline{y}+{\overline{p}}_{{\overline{x}}^{\prime }}{\overline{x}}^{\prime }+{\overline{p}}_{{\overline{y}}^{\prime }}{\overline{y}}^{\prime },\end{array}$$

(5)

where

$${\overline{p}}_{\overline{x}}=\frac{\partial \overline{p}}{\partial \overline{x}},{\overline{p}}_{\overline{y}}=\frac{\partial \overline{p}}{\partial \overline{y}},{\overline{p}}_{{\overline{x}}^{\prime }}=\frac{\partial \overline{p}}{\partial {\overline{x}}^{\prime }},{\overline{p}}_{{\overline{y}}^{\prime }}=\frac{\partial \overline{p}}{\partial {\overline{y}}^{\prime }},$$

(6)

where ${\overline{p}}_{\overline{x}}$ and ${\overline{p}}_{\overline{y}}$ are the displacements of the oil film pressure in x and y directions, ${\overline{p}}_{{\overline{x}}^{\prime }}$ and ${\overline{p}}_{{\overline{y}}^{\prime }}$ are the velocities of the oil film pressure in x and y directions, ${\overline{h}}_0$ is the thickness of the oil film in a static equilibrium condition, and ${\overline{p}}_0$ is the pressure of the oil film pressure in a static equilibrium condition.

The above Equations (5) and (6) are utilized to redescribe Equation (4), and the following dimensionless disturbance equation is obtained after ignoring the high-order term:

$$\left[\frac{\partial }{\partial \theta }({\overline{h}}_0^3\frac{\partial }{\partial \theta })+\frac{1}{{\lambda }^2}\frac{\partial }{\partial \overline{z}}({\overline{h}}_0^3\frac{\partial }{\partial \overline{z}})\right]\left\{\begin{array}{c}{\overline{p}}_{\overline{x}}\\ {\overline{p}}_{\overline{y}}\\ {\overline{p}}_{{\overline{x}}^{\prime }}\\ {\overline{p}}_{{\overline{y}}^{\prime }}\end{array}\right\}=\left\{\begin{array}{c}3(\cos \theta -3\frac{\sin \theta }{{\overline{h}}_0}\frac{\partial {\overline{h}}_0}{\partial \theta })-3{\overline{h}}_0^3\frac{\partial }{\partial \theta }(\frac{\sin \theta }{{\overline{h}}_0})\frac{\partial {\overline{p}}_0}{\partial \theta }\\ -3(\sin \theta +3\frac{\cos \theta }{{\overline{h}}_0}\frac{\partial {\overline{h}}_0}{\partial \theta })-3{\overline{h}}_0^3\frac{\partial }{\partial \theta }(\frac{\cos \theta }{{\overline{h}}_0})\frac{\partial {\overline{p}}_0}{\partial \theta }\\ 6\sin \theta \\ 6\cos \theta \end{array}\right\},$$

(7)

Thus, after the ${\overline{p}}_{\overline{x}}$, ${\overline{p}}_{\overline{y}}$, ${\overline{p}}_{{\overline{x}}^{\prime }}$ and ${\overline{p}}_{{\overline{y}}^{\prime }}$ are calculated by the integral operation, the dimensionless stiffness and damping coefficients of the sliding bearing can be solved according to Equations (8) and (9), respectively.

$$\left[\begin{array}{cc}{\overline{k}}_{xx}& {\overline{k}}_{xy}\\ {\overline{k}}_{yx}& {\overline{k}}_{yy}\end{array}\right]=-{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_0^{2\pi }\left[\begin{array}{c}\sin \theta \\ \cos \theta \end{array}\right]}}\left[\begin{array}{cc}{\overline{p}}_{\overline{x}}& {\overline{p}}_{\overline{y}}\end{array}\right]d\theta d\overline{z},$$

(8)

$$\left[\begin{array}{cc}{\overline{c}}_{xx}& {\overline{c}}_{xy}\\ {\overline{c}}_{yx}& {\overline{c}}_{yy}\end{array}\right]=-{\displaystyle {\int }_{-1}^1{\displaystyle {\int }_0^{2\pi }\left[\begin{array}{c}\sin \theta \\ \cos \theta \end{array}\right]}}\left[\begin{array}{cc}{\overline{p}}_{{\overline{x}}^{\prime }}& {\overline{p}}_{{\overline{y}}^{\prime }}\end{array}\right]d\theta d\overline{z},$$

(9)

where k and c represent the stiffness and damping coefficients, and the subscripts x and y mean the corresponding x and y directions.

According to Equations (1)–(9), the length L_b, clearance ratio φ, ratio of length to diameter λ and oil film viscosity coefficient η are the main factors affecting the stiffness and damping coefficients of the sliding bearing. Then, the L_b, φ, λ and η are taken as the uncertain structural parameters of the sliding bearing in this paper, and the effects of their variations on the dynamic characteristics of the sliding bearing-rotor system of a multi-stage centrifugal pump are studied.

2.2. The Finite Element Modeling of the Multi-Stage Centrifugal Pump Rotor System

The prerequisite to precisely predict the dynamic characteristics of the multi-stage centrifugal pump rotor system is constructing an accurate dynamic model, which can well reflect the complex structure and actual working conditions. In this paper, the ANSYS software is adopted to construct the finite element model of the multi-stage centrifugal pump in the virtual environment. When constructing the finite element model, the small holes, keyways, chamfers on the rotating shaft and impeller are ignored, and the following simplification principles are also considered:

• The centrosymmetric pump shaft is simplified as a circular section beam element, and the parameters such as the element density, elastic modulus and Poisson’s ratio are defined.
• The multiple impellers are equivalent to the centralized mass points by applying the Mass 21 units at their positions. The translational mass and moment of inertia for a Mass 21 unit are set through its six real constants as MASSX, MASSY, MASSZ, IXX, IYY and IZZ.
• The spring-damping element COMBI214 is adopted to simulate the dynamic characteristics of the sliding bearing. The stiffness and damping coefficients identified through the method presented in Section 2.1 are set in the COMBI214 to simulate the support of the rotating pump shaft. The COMBI214 unit has two nodes, and the four stiffness coefficients and four damping coefficients of the sliding bearing are set through K11, K22, K12, K21, C11, C22, C12 and C2. K11 and K22 are the main stiffness values, K12 and K21 are the cross stiffness values, C11 and C22 are the main damping values, and C12 and C21 are the cross damping values.

On the basis, a schematic diagram is provided in Figure 3 to describe the equivalent approach for the multi-stage centrifugal pump rotor system. In Figure 3, the red points mean the Mass 21 units which are used to represent the impellers at these positions, and the COMBI214 spring-damping elements represent the dynamic characteristics of the sliding bearings.

3. Dynamic Characteristic Analysis of the Multi-Stage Centrifugal Pump Rotor System with Uncertain Sliding Bearing Structural Parameters

The dynamic characteristic analysis of the multi-stage centrifugal pump rotor system mainly focuses on predicting the critical speeds, which can benefit the selections of the operational speeds to avoid the resonance and guarantee the stable operation. However, the manufacturing errors can cause the variations of the sliding bearing structural parameters around their nominal values, and lead to changes of the stiffness and damping coefficients of the sliding bearing. As a consequence, the critical speeds of the multi-stage centrifugal pump rotor system will vary within a certain range. Considering this situation, an optimization model is established in this section to obtain the extreme value of the critical speed within the variation ranges of the sliding bearing structural parameters. The improved particle swarm optimization algorithm is introduced to solve this optimization model based on the co-simulations of the MATLAB and ANSYS. The Monte Carlo algorithm is used to verify the accuracy of the obtained extreme value.

3.1. Establishing an Optimization Model Based on the Critical Speed of the Rotor System

3.1.1. The Design Variables

The length L_b, clearance ratio φ, ratio of length to diameter λ and oil film viscosity coefficient η are the maim structural parameters of the sliding bearing. The manufacturing errors cause the variations of these structural parameters, and then result in uncertain dynamic characteristics of the multi-stage centrifugal pump rotor system. Therefore, these structural parameters of the sliding bearing V_S = (L_b, φ, λ, η) are taken as the design variables to obtain robust critical speeds of the rotor system.

3.1.2. The Optimization Objective

The extreme value of the critical speed of the multi-stage centrifugal pump rotor system is defined as the optimization objective. At present, the critical speeds of the multi-stage centrifugal pump rotor system are mainly obtained by establishing the finite element model of the rotor system in simulation software and obtaining the Campbell diagram through the modal analysis. According to Section 2.2, the prerequisite of establishing the finite element model is to obtain the 8 stiffness and damping coefficients, which represents the dynamic characteristics of the sliding bearing. For each combination of sliding bearing structural parameters, its corresponding 8 dynamic coefficients can be calculated by writing and solving the programs of the Reynolds equation in the MATLAB software.

Therefore, in order to establish the direct relationship between the sliding bearing structural parameters and the critical speeds of the rotor system, the ANSYS and MATLAB are combined to compute the critical speeds automatically when the sliding bearing structural parameters change. The co-computations are realized based on the re-development language of ANSYS and the data exchange interface of the ANSYS and MATLAB. The main control program to calculate the 8 dynamic coefficients of the sliding bearing is written in the MATLAB software. The calculated dynamic characteristic coefficients are transferred into the ANSYS software through the ANSYS parametric design language (APDL). Then, the corresponding finite element model of the multi-stage centrifugal pump rotor system can be established. Subsequently, the modal analysis considering the gyroscopic effect is carried out on the finite element model to yield the Campbell diagram shown in Figure 4, where the abscissa represents the speed and the ordinate represents the frequency. The ray from the origin intersects with the other lines, and the intersections stand for the natural frequencies [22]. The natural frequency f and the critical speed n have the following mathematic relationship:

$$n=60\cdot f,$$

(10)

Then, the critical speeds output by the ANSYS software are transferred to the MATLAB software for the subsequent analysis.

In summary, the optimization model for obtaining the extreme value of the critical speed of the multi-stage centrifugal pump rotor system is expressed in Equation (11):

$$\begin{array}{l}{f}_{\min speed}(V_s)=\min (n)or{f}_{\max speed}(V_s)=\max (n),\\ V_s=(L_b\text{, }f\text{, }\lambda \text{, }\eta {)}^T,\\ s.t.\{\begin{array}{c}{L}_{b_{\min }}\le L_b\le L_{b_{\max }}\\ f_{\min }\le f\le f_{\max }\\ {\lambda }_{\min }\le \lambda \le {\lambda }_{\max }\\ {\eta }_{\min }\le \eta \le {\eta }_{\max }\end{array},\end{array}$$

(11)

3.2. The Improved PSO Algorithm for Solving the Optimization Model

As a numerical optimization algorithm based on the swarm intelligence, PSO algorithm can effectively solve the optimization problem of complex functions [23,24]. Thus, an improved PSO algorithm is introduced to solve the optimization model for obtaining the extreme value of the critical speed in this research. For a PSO algorithm, each particle has two attributes as the position vector and velocity vector. The position represents the direction of flight and the particle flies to a new position according to the velocity vector. Assuming that there are N D-dimensional particles forming a community, the i_th particle with the position vector X_i = ($x_{i1}$, $x_{i2}$, ···, $x_{id}$) flies to a better position with the velocity vector V_i = ($v_{i1}$, $v_{i2}$, ···, $v_{id}$). During the optimization process, the position vector and velocity vector are updated iteratively based on the individual optimal position and the global optimal position described in Equation (12).

$$\{\begin{array}{l}{v}_{id}^{k+1}=\omega v_{id}^k+c_1{r}_1\left({{\displaystyle P}}_{bes{t}_{id}}^k-x_{id}^k\right)+c_2{r}_2\left({{\displaystyle G}}_{bes{t}_{id}}^k-x_{id}^k\right)\hfill \\ x_{id}^{k+1}=x_{id}^k+v_{id}^{k+1}\hfill \end{array},$$

(12)

where k is the current iteration number, ω is the inertial weight, $c_1$ and $c_2$ are the learning factors, $r_1$ and $r_2$ are the random parameters varying in the range of 0 to 1, $v_{id}$ means the d_th element in the velocity vector of the i_th particle, and $x_{id}$ means the d_th element in the position vector of the i_th particle, P_{best_id} is the d_th element in the individual optimal position vector of the i_th particle determined by comparing the fitness values, and G_{best_id} is the d_th element in the global optimal position vector determined by comparing the fitness values of the N particles. The fitness value indicates the superior or inferior of the current particle position, and it can be calculated from the fitness function. Moreover, the inertia weight ω is dynamically modified as expressed in Equation (13), for improving the global and local search ability at different stages of the iterative process.

$$\omega ={\omega }_{\max }-k\left({\omega }_{\max }-{\omega }_{\min }\right)/k_{\max },$$

(13)

where k_max is the maximum iteration number, and ω_min and ω_max stand for the defined minimum and maximum inertial weights respectively.

According to the aforementioned theory of the PSO algorithm and the specific optimization problem described in Section 3.1, the sliding bearing structural parameters V_s = (L_b, φ, λ, η) are defined as the position vector of one particle, and the fitness value is the critical speed of the rotor system calculated using the related position vector X. Then, the optimal critical speed of the rotor system can be obtained through the following steps:

• The basic parameters for carrying out the optimization are initialized, including the number of particles, the maximum iteration number, the position and velocity vectors for each particle, the learning factors and the extreme values of the inertial weights.
• In each iteration, the critical speeds corresponding to the position vector of each particle are calculated with the co-simulations of the MATLAB and ANSYS software, and then the individual optimal position for each particle and the global optimal position among all the particles can be determined by comparing the critical speeds.
• After each iteration, the position and velocity vectors for each particle are updated by Equation (12), and these updated vectors are utilized to find better individual and globe optimal positions by repeating the step 2.
• The iteration terminates when the convergence condition is met or the iteration number reaches its maximum value. Then, the fitness value corresponding to the global optimal position is the extreme critical speed within the variation ranges of the sliding bearing structural parameters.

3.3. The Monte Carlo Algorithm for the Validation

The extreme values of the critical speed can be obtained by solving the optimization models according to Section 3.1 and Section 3.2. To verify whether the obtained values belong to local optimums and guarantee their accuracies, the Monte Carlo algorithm is adopted in this paper to generate a large number of sample data within the variation ranges of the independent variables. The Monte Carlo algorithm is a numerical sampling process, belonging to a branch of stochastic mathematics. The algorithm has the characteristic of strong universality, and its calculation accuracy will be improved with the increase of simulation times. Since the repeatedly calculations under a large number of random samples are time-consuming, the Monte Carlo algorithm is mainly used for a verification in practical engineering. A brief introduction of the random variable sampling process in the Monte Carlo method is described as follows [25,26].

The random variable β within the interval [a, b] is assumed to meet the distribution of the probability density function p(β), which has the following normalization condition:

$${\displaystyle {\int }_a^{b}p}(\beta )\mathrm{d}\beta =1,$$

(14)

In order to estimate the random variable, the pseudo-random number evenly distributed between 0 and 1 is generated, and then the cumulative distribution function of is expressed in Equation (15).

$$F_{\xi }(\xi )=\{\begin{array}{ll}1\hfill & ,\xi >1\hfill \\ \xi \hfill & ,0<\xi \le 1\hfill \\ 0\hfill & ,\xi \le 0\hfill \end{array},$$

(15)

If a non-decreasing function can be found to describe the mapping relationship between the ξ and β as β = f(ξ), the following probability relationship can be obtained.

$$p\left(a<\beta \le {\beta }_1\right)=p\left(0<\xi \le {\xi }_1\right),$$

(16)

Then, the cumulative distribution function in Equation (15) and the probability density function in Equation (16) are combined to yield the following equation.

$$F_{\beta }\left({\beta }_1\right)=F_{\xi }\left({\xi }_1\right),$$

(17)

For Equation (17), its left side is further expressed by the integral of the probability density function, and its right side is substituted by the equation in Equation (15). Then, a following relationship between the β and $\xi$ can be obtained as described in Equation (18).

$${\displaystyle {\int }_a^{{\beta }_1}p}(\beta )\mathrm{d}\beta ={\xi }_1,{\xi }_1\in (0,1),$$

(18)

Thus, a β₁ can be obtained by Equation (18) once a random number ${\xi }_1$ is generated.

Based on the probability and statistics theory of the Monte Carlo method, its application in validating the accuracy of the optimization model in this research can be concluded as the following steps:

• Taking the sliding bearing structural parameters L_b, φ, λ and η as the random variables, and M random combinations of the structural parameters are generated to calculate the dynamic characteristic coefficients of the sliding bearing.
• These calculated stiffness and damping coefficients are written into the finite element model of the multi-stage centrifugal pump rotor system to predict the corresponding critical speeds.
• The obtained critical speeds of the multiple samples are compared to determine the maximum and minimum values, which are further compared with those calculated by solving the optimization models with the improved PSO algorithm.

4. Case Study and Results Discussion

In this section, the previously presented method was applied to a real multi-stage centrifugal pump to validate its feasibility. Effects of the variations of the sliding bearing structural parameters on the dynamic characteristics of the multi-stage centrifugal pump rotor system were investigated through the co-simulations of the MATLAB and ANSYS. Then, optimization models were established and solved by the improved PSO algorithm to obtain the optimal critical speeds within the variation ranges of the sliding bearing structural parameters. The accuracies of the obtained extreme values of the critical speeds were validated through the Monte Carlo algorithm and the rated speed of the multi-stage centrifugal pump.

4.1. Effects of the Structural Parameters on the Dynamic Characteristics of the Sliding Bearing

The studied sliding bearing is described in Figure 5, and the nominal values of its structural parameters are L_b = 0.45 m, φ = 0.002, λ = 2.04 and η = 0.0016 Pa·s. The nominal values were first used to calculate the stiffness and damping coefficients of the sliding bearing according to the method provided in Section 2.1, and the dynamic characteristic coefficients computed through the programs in the MATLAB virtual environment are provided in

Table 1. Dynamic characteristic coefficients of the sliding bearing.

Stiffness Coefficients (×107 N/m)				Damping Coefficients (×105 N·s/m)
Kxx	Kxy	Kyx	Kyy	Cxx	Cxy	Cyx	Cyy
1.09	−2.18	2.30	1.23	7.13	3.38	−3.75	8.87

.

Furthermore, the structural parameters containing the length L_b, clearance ratio φ, ratio of length to diameter λ and oil film viscosity coefficient η were taken as the variables, and the effects of the variations of the structural parameters on the dynamic characteristic coefficients of the sliding bearing were studied following the single variable principle. In Figure 6, the stiffness and damping coefficients of the sliding bearing increase with the increase of the length L_b. In Figure 7, the stiffness and damping coefficients of the sliding bearing decrease with the increase of the clearance ratio φ. In Figure 8, the stiffness and damping coefficients of the sliding bearing increase with the increase of the ratio of the length to diameter λ. In Figure 9, the stiffness and damping coefficients of the sliding bearing increase with the increase of the viscosity coefficient η. Therefore, it is significant to take the uncertain sliding bearing structural parameters into consideration, since their variations can lead to the changes of the dynamic characteristic coefficients.

4.2. Dynamic Characteristic Analysis of the Multi-Stage Centrifugal Pump Rotor System

The specific parameters of the impeller and shaft section of the multi-stage centrifugal pump rotor system are described in

Table 2. The specific parameters of the impellers.

Impeller Series	Axial Coordinate(m)	Mass(kg)	Polar Moment of Inertia (kg·m2)	The Moment of Inertia about a Diameter (kg·m2)
1	0.37	13	0.1	0.19
2	0.47	12	0.09	0.18
3	0.57	14	0.1	0.19
4	0.97	13	0.1	0.19
5	1.07	12	0.09	0.18
6	1.17	13	0.09	0.19

and

Table 3. Dimension parameters of the rotor shaft section.

Order Number	Start-Point Coordinate (m)	End-Point Coordinate (m)	Outer Diameter (m)
1	0	0.25	0.016
2	0.25	0.37	0.028
3	0.37	1.17	0.036
4	1.17	1.29	0.028
5	1.29	2.04	0.016

, respectively. One the basis, the beam elements, spring-damping elements and concentrated mass elements were integrated to establish the dynamic model of the multi-stage centrifugal pump rotor system in the ANSYS virtual environment as shown in Figure 10. The 8 stiffness and damping coefficients of each sliding bearing are solved as listed in Table 1, which were written into the corresponding spring-damping elements. Then, the modal analysis was carried out to obtain the Campbell diagram shown in Figure 4 for predicting the critical speeds of the multi-stage centrifugal pump rotor system.

Obeying the single variable principle, effects of the variations of the siding bearing length L_b, clearance ratio φ, ratio of length to diameter λ and oil film viscosity coefficient η on the dynamic characteristics of the multi-stage centrifugal pump rotor system were also investigated. The first three critical speeds were selected to represent the dynamic characteristics of the rotor system. In Figure 11, the first three critical speeds increase with the increase of the bearing length L_b, and the first critical speed is more sensitive to the variation of L_b. In Figure 12, the first critical speed decreases with the increase of the clearance ratio φ, and the second and third critical speeds increase with the increase of the clearance ratio φ. In Figure 13, the first two critical speeds increase with the increase of the ratio of length to diameter λ, and the third critical speed decreases with the increase of the ratio of length to diameter λ. However, the variation of λ has slight effects on all the three critical speeds. In Figure 14, all the first three critical speeds increase with the increase of the oil film viscosity coefficient η, and the first critical speed is more vulnerable to the variation of η since it has a greater rate of change.

4.3. Robust Dynamic Characteristic Analysis of the Rotor System Considering Uncertain Sliding Bearing Structural Parameters

The simulated results in Section 4.2 show that the sliding bearing structural parameters have important impacts on the dynamic characteristics of the multi-stage centrifugal pump rotor system. Thus, it is necessary to take the variations of the sliding bearing structural parameters caused by the manufacturing errors into account when predicting the critical speeds. Therefore, the length L_b, clearance ratio φ, ratio of length to diameter λ and oil film viscosity coefficient η were taken as the independent variable, and the 5% and −5% of the nominal value of each structural parameter were calculated to define the boundary values. Then, the corresponding variation ranges of the sliding bearing structural parameters L_b, φ, λ and η are [0.43 m, 0.47 m], [0.0019, 0.0021], [1.94, 2.14], [0.0015 Pa.s, 0.0017 Pa.s] respectively. Considering the first three critical speeds, 6 optimization models were established for obtaining the minimum and maximum values of each critical speed within the variation ranges of these sliding bearing structural parameters.

To carry out the improved particle swarm optimization algorithm for solving the established optimization models, the number of particles was 80, the maximum number of iterations was 150, the learning factors C₁ and C₂ were defined as 1.5, respectively, and the maximum and minimum inertia weights were 1.2 and 0.6, respectively. During the optimization process, the MATLAB and ANSYS were combined to calculate the critical speed of each particle. The iterative process for obtaining the maximum value of the first-order critical speed is shown in Figure 15, where the optimal critical speed trends to a stable value after 58 iterative calculations. The finally obtained extreme values of the first three critical speeds are listed in

Table 4. The ranges of the critical speeds calculated by different methods.

Different Methods	Critical Speeds(rpm)
	The First Order	The Second Order	The Third Order
The improved PSO algorithm	[2451.3,2505.7]	[4137.6,4177.3]	[11,527.8,11,541.6]
The Monte Carlo method	[2449.9,2509.5]	[4132.1,4183.4]	[11,525.1,11,544.8]

, and obvious differences between the maximum and minimum values can be observed. Comparing the extreme values in Table 4, the deviations of the first three critical speeds are 54.4 rpm, 39.7 rpm and 13.8 rpm respectively, and the corresponding variation rates are 2.22 %, 1.0% and 0.12%. Thus, the first-order critical speed shows more vulnerable to the variations of the sliding bearing structural parameters.

In all, 10⁶ random combinations of the structural parameters L_b, φ, λ and η were defined by the Monte Carlo method, and the corresponding maximum and minimum values of the first three critical speeds are also listed in Table 4. It can be seen from Table 4 that the maximum error between the extreme critical speeds calculated by the Monte Carlo method and those obtained by solving the optimization model is within 1%, indicating that the improved particle swarm optimization algorithm for obtaining the optimal critical speed can not only improve the calculation efficiency by reducing the number of samples but also have a desired accuracy. Moreover, in the actual operating conditions, the rated working speed of the rotor system of the multi-stage centrifugal pump is 1980 rpm. It is far from the first-order critical speed range [2451.3 rpm, 2505.7 rpm] calculated in this paper, avoiding the resonance and guaranteeing a safe operation.

5. Conclusions

In this paper, a new method was proposed to predict the critical speeds of the multi-stage centrifugal pump with uncertain sliding bearing structural parameters. The structural parameters are taken as the variables to establish an optimization model for obtaining the extreme value of the critical speed, which can be solved by the improved particle swarm optimization algorithm based on the co-simulations of the MATLAB and ANSYS. Therefore, variation ranges of the critical speeds of the focused orders can be obtained, laying a basis for selecting robust speeds that are away from the variation ranges to guarantee the stable operation of the multi-stage centrifugal pump under actual working conditions. The main contributions of this paper include:

• The oil film pressure distribution of the sliding bearing was first solved by combing the finite difference method and the successive over-relaxation method. Then, the small parameter method was used to compute the dynamic characteristic coefficients of the sliding bearing in the virtual environment of the MATLAB. On that basis, effects of the variations of the sliding bearing structural parameters L_b, φ, λ and η on the stiffness and damping coefficients were investigated.
• The beam elements, spring-damping elements and concentrated mass elements were combined to construct the finite element model of the multi-stage centrifugal pump rotor system in the virtual environment of ANSYS. The computed dynamic characteristic coefficients of the sliding bearing were assigned to the spring-damping elements. Moreover, the MATLAB and ANSYS were combined to investigate the effects of the variations of each structural parameter on the first three critical speeds of the rotor system.
• Six optimization models for obtaining the extreme values of the first three critical speeds were established, where the sliding bearing structural parameters L_b, φ, λ and η were taken as the independent variables. These optimization models were solved by the improved PSO algorithm within the variation ranges of the structural parameters, and the maximum and minimum values of each critical speed were obtained. These extreme values were close to those calculated with the 10⁶ random samples generated by the Monte Carlo method. The deviation between the maximum and minimum values of the first-order critical speed is 54.5 rpm, showing the necessity for considering the effects of the manufacturing errors on the dynamic characteristics of the rotor system. Moreover, the actual operational speed is away from the variation range of the first-order critical speed, indicating that the proposed method is promising in selecting robust operational speed of a multi-stage centrifugal pump.