*Article* **Effect of Impeller Diameter on Dynamic Response of a Centrifugal Pump Rotor**

**Alireza Shooshtari 1,\*, Mahdi Karimi 1, Mehrdad Shemshadi <sup>1</sup> and Sareh Seraj <sup>2</sup>**


**\*** Correspondence: shooshta@basu.ac.ir

**Abstract:** This paper investigates the effect of impeller diameter on the dynamic response of a centrifugal pump using an inverse dynamic method. For this purpose, the equations of motion of the shaft and the impeller are derived based on Timoshenko beam theory considering the impeller as a concentrated mass disk. For practical modeling, the model of Jones and Harris is added to the equation to include the effect of bearings. As a case study, the model is applied to a process pump used in an oil refinery. Computing the eigenvalues of the model and comparing them with the natural frequencies of the structure, the model updating of the problem is performed through an indirect method. Three impellers with different diameters are applied to the updated model. The results show that increasing the diameter of the pump impeller can increase the amplitude of vibration up to 52% at critical speeds of the rotor. It is found that in addition to the hydraulic condition and efficiency, the impeller diameter should be considered as an important factor in the selection of centrifugal pumps.

**Keywords:** rotor dynamic; bearing; centrifugal pump; impeller diameter; Lagrangian equations

## **1. Introduction**

Rotating machineries such as pumps, compressors, turbines, etc., play an important role in many different industries. Accurate predictions of rotor system dynamic characteristics are very important in the design of any type of machines. There have been many studies relating to the field of rotor dynamic systems during the recent years. Engineers have developed several new techniques about the dynamics of rotating machines.

The first recorded theory of rotor dynamics was in a classic paper of Jeffcott [1]. The Jeffcott rotor model has been used to explain the whirling effect. It consists of a simply supported flexible massless shaft with a rigid disc mounted at the mid-span. In the Jeffcott model, the moments of inertia *Ip* and *It* are not considered. This is because there are no gyroscopic moments exerted on the shaft. The disc is assumed to move in a plane that is perpendicular to the shaft spin axis. By developments in the technology of rotating machines, the rotational speed of rotors became higher, and so, the non-conservative forces generated through the bearings of the rotor become considerable. To determine the critical speeds in which resonance has occurred, it is necessary to know the natural frequencies, mode shapes, and forced responses, which are caused by unbalancing in rotor systems. Prohl studied the critical speed evaluation of a turbine shaft, and he suggested the transfer matrix method [2]. The first application of the finite element method for a rotor system was made by Ruhl and Booker [3]. In their study, the influences of the rotary inertia, gyroscopic moment, bending, shear deformation, axial load, and internal damping were neglected to simplify the model. The theory has been developed by considering the rotary inertia, gyroscopic moment, and axial forces. Nelson and McVaugh extended this to include gyroscopic effects. They derived the equations of motion for the shaft and the effects of translational and rotary inertia and gyroscopic moments on it [4]. Erturka et al.

**Citation:** Shooshtari, A.; Karimi, M.; Shemshadi, M.; Seraj, S. Effect of Impeller Diameter on Dynamic Response of a Centrifugal Pump Rotor. *Vibration* **2021**, *4*, 117–129. https://doi.org/10.3390/ vibration4010010

Academic Editors: Hamed Kalhori and Aleksandar Pavic Received: 26 November 2020 Accepted: 26 January 2021 Published: 9 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

presented an analytical method based on Timoshenko beam theory for calculating the frequency response function (FRF) of a spindle–holder–tool combination. They proposed a mathematical model and obtained the point FRF for a tool [5]. Subbiah et al. showed that a rotor has certain speed ranges in which a large amplitude of vibration could occur. These speed ranges are known as critical speed, which results in excessive rotor deflection [6]. Phadatare and Barun described a step forward in calculating the nonlinear frequencies and resultant dynamic behavior of a high-speed rotor bearing system with unbalanced mass. In this study, Fast Fourier Transform (FFT) analysis was established for finding the fundamental frequencies of the rotor according to the variation of shaft diameter and location of unbalanced mass [7]. Metsebo et al. have focused on the influence of the rotating shaft on the dynamic of a rotor ball bearing system. They carried out a mathematical modeling for the system considering the shaft as a Timoshenko beam [8].

Ball bearings are the essential elements of rotating machineries. So, the influence of bearings on the performance of rotor-bearing systems is very important. El-Sayed derived a set of equations for the stiffness of bearings using the Hertz theory and determined the total deflections of inner and outer races caused by an applied load [9]. Tamura and Tsuda performed a theoretical study about fluctuations of radial spring characteristics of a ball bearing due to ball revolutions [10]. Many researchers also estimated bearing stiffness by carrying out some experiments. Stone and Walford developed a rotor-bearing test rig to estimate the bearing's radial stiffness and damping by measuring the response of the rotor [11]. Jairo et al. presented an experimental validation for a mathematical modeling of ball bearing. In this model, the bearing was considered as a mass-springdamper system based on Hertz equations for contact deformation [12]. Xia sheng et al. proposed a mathematical model for the stiffness of bearings, which is varying by speed. They explained that the speed of a rolling bearing varies the stiffness of the bearing [13]. Zhang et al. investigated the effect of ring misalignment on the service characteristics of ball bearing and rotor systems [14]. They improved a quasi-static model of ball bearing considering the misalignment in its ring. Neisi et al. worked on the effect of off-sized balls on contact stresses in a touchdown bearing [15]. Yi Qin et al. developed a dynamics model for deep groove ball bearings with local faults based on coupled and segmented displacement excitation [16]. Chandrasekaran et al. used computational fluid dynamic (CFD) methods and mathematical modeling to investigate the impeller design parameters on the effect of fluid follows in the centrifugal pumps [17].

The objective of this paper is investigation on the effect of impeller diameter on the amplitude of vibration at critical speeds in an overhung centrifugal pump. Modeling of the shaft and impeller is based on Timoshenko theory, and modeling of the bearing is based on the work of Jones and Harris, which is added to the model of shaft simultaneously. Then, numerical analysis for a real centrifugal pomp in the oil refinery has been done based on the proposed model. It has shown that in addition of the effect of the geometrical parameter of the shaft, the effect of the diameter of the impeller on the dynamical behavior of the pump is very important.

The novelty of this research is the development of the mathematical model of the shaft, impeller, and bearing simultaneously. Using this developed model, the effect of impeller size diameter on the dynamic behavior of a centrifugal pump has been investigated, and it has been shown that it is very considerable.

This paper includes four sections. In the first section, the literature has been reviewed and the necessity of research has been recognized. In the second section, the model of the system has been introduced and the equation of motion based on the energy method has been derived by calculating the potential and kinetic energy of the system. In the third section, the model of the bearing has been presented, and the governing equations of motion for the ball bearing have been introduced. Finally, in the last section, the derived equations have been used and solved.

For an actual centrifugal pump, the effect of impeller diameter on the amplitude of vibration has been investigated.

#### **2. System Modeling**

The rotor of a centrifugal pump is approximated by a simple model as shown in Figure 1a. The model is composed of a shaft of length *L* and supported by two bearings located at *L*<sup>1</sup> and *L*<sup>2</sup> along the shaft. *u*, *v*, and *w* reflect the displacement in the *ex*, *ey*, and *ez* directions, respectively.

**Figure 1.** Model of centrifugal pump (**a**) definitions; (**b**) force of the bearing on the rotor system.

The shaft is modeled as a Timoshenko beam. In this model, first-order shear deformation theory with rotary inertia and gyroscopic effect has been considered. The shaft rotates at a constant speed around its longitudinal axis. In addition, it has a uniform, circular cross-section.

The equations of motion of the system are obtained using the Lagrangian equations as:

$$\frac{d}{dt}\frac{\partial T}{\partial \dot{q}\_n} - \frac{\partial T}{\partial q\_n} + \frac{\partial I}{\partial q\_n} = (Q\_{nc})^T \tag{1}$$

where *T* is the total kinetic energy of the system, *U* is the total potential energy of the system, *qn* is the generalized coordinate and *Qnc* represents non-conservative forces that are not directly related to the potential energy of the system. So, to derive the governing equations of motion using Equation (1), one must calculate the kinetic and potential energy of the system.

The total kinetic energy of a rotor system is estimated by the dynamic motion of the shaft and disk [14].

$$T = T\_{disk} + T\_{shaft} \tag{2}$$

The kinetic energy due to the rotation of the disk is difficult to calculate. Therefore, we assume that the disk is symmetric and rigid and has been fixed at the end of the shaft. The motion of the disk can be defined as two superimposed rotations *θx*, *θ<sup>y</sup>* and two translational deflections *u*, *v* in directions *ex*, *ey*. So, the kinetic energy of the disk can be expressed as:

$$T\_{disk} = \frac{1}{2}M\left(\dot{\boldsymbol{u}}^2 + \dot{\boldsymbol{v}}^2\right) + \frac{1}{2}I\_l\left(\dot{\boldsymbol{\theta}}\_x^2 + \dot{\boldsymbol{\theta}}\_y^2\right) - I\_p\Omega\dot{\boldsymbol{\theta}}\_x\boldsymbol{\theta}\_y + \frac{1}{2}\left.I\_pL\Omega^2\right.\tag{3}$$

where *M* is the mass of the rigid disk, *It* is the diametral mass moment of inertia, and *Ip* is the polar mass moment of inertia. The term *Ip*<sup>Ω</sup> . *θxθ<sup>y</sup>* represents the gyroscopic effect, and finally, the last term defines the kinetic energy due to the rotation of the disk.

The kinetic energy of the shaft involves the kinetic energy due to bending of the shaft, effect of rotatory inertia, and gyroscopic effect. The kinetic energy of the shaft can be derived as:

$$T\_{shaft} = \int\_0^l \frac{1}{2} \rho\_{shaft} A(\dot{\hat{u}}^2 + \dot{\hat{v}}^2) dz + \int\_0^l \frac{1}{2} l\_l (\dot{\hat{\theta}}\_x^2 + \dot{\hat{\theta}}\_y^2) dz + \frac{1}{2} l\_l L \Omega^2 - \int\_0^l l\_p \Omega \dot{\hat{\theta}}\_x \theta\_y \, dz \tag{4}$$

where *A* is the cross-sectional area, *ρshaf t* is the density of the shaft, and *Jt* and *Jp* are the diametric and polar inertia of the shaft, respectively. So, substituting Equations (3) and (4) in Equation (2), the total kinetic energy of system has been obtained.

The potential energy of the system includes the strain energy due to the deformation of the shaft (*U*1) and the strain potential energy due to the deflection of the bearing installed on the shaft (*U*2) [12]:

$$
\mathcal{U}I = \mathcal{U}I\_1 + \mathcal{U}I\_2.\tag{5}
$$

The strain energy of the shaft can be expressed as:

$$dI\_1 = \frac{1}{2} \int\_0^l EI((\frac{\partial^2 u}{\partial^2 z})^2 + (\frac{\partial^2 v}{\partial^2 z})^2) dz + \frac{1}{2} \int\_0^l kGA\left(\Phi\_x^2 + \Phi\_y^2\right) dz\tag{6}$$

where *E* represents the modulus of elasticity, *G* is the shear modulus, and *k* is the shear coefficient. In addition, *A* and *I* are the cross-section area and moment of inertia of the shaft receptivity. In the above equation, the first term is related to the shaft bending, and the second term is due to shear deformation. In addition, the potential energy caused by bearing forces can be expressed as:

$$\mathcal{U}\_2 = \left(-F\_{1x}\mu - F\_{1y}\upsilon\right)\_{l\_1} + \left(-F\_{2x}\mu - F\_{2y}\upsilon\right)\_{l\_2}.\tag{7}$$

Different forces work on the impeller as a consequence of the fluid. These forces are non-conservative forces and are unknown, and they have been ignored for simplicity. However, mass unbalance generates an additional centrifugal force that makes it possible to calculate this non-conservative force. If the impeller is out of balance, the resulting centrifugal force will induce the rotor to vibrate. When the shaft rotates at a speed equal to the natural frequency, this vibration becomes large. Unbalance is defined by a small mass *mu* situated at a distance *du* from the geometric center of the impeller, as shown in Figure 2 [13].

The out of balance force at the end of the rotor is:

$$F\_{\mathfrak{U}} = m\_{\mathfrak{U}} \, d\_{\mathfrak{U}} \, \Omega^2. \tag{8}$$

This rotating force can be resolved into two components in the *x* and *y* direction as:

$$\begin{cases} F\_{\mu-x} = m\_{\mu} \, d\_{\mu} \, \Omega^2 \cos \Omega t \\ F\_{\mu-x} = m\_{\mu} \, d\_{\mu} \, \Omega^2 \sin \Omega t \end{cases} \tag{9}$$

**Figure 2.** Mass unbalance.

It is assumed that the unbalance mass is in the *x* direction in the initial state. The deflections in the *x* and *y* directions are expressed as:

$$\begin{array}{l} u(z.t) = f(z)q\_1(t) = f(z)q\_1\\ v(z.t) = f(z)q\_2(t) = f(z)q\_2 \end{array} \tag{10}$$

where *q*<sup>1</sup> and *q*<sup>2</sup> are generalized independent coordinates and *f*(*z*) is the displacement function that satisfies the boundary conditions of the system. As the rotor of the centrifugal pump has simply supported at both ends, *f*(*z*) has been selected as [12]:

$$f(z) = \sin\left(\frac{n\,\pi}{L\_2 - L\_1}z - L\_1\right). \tag{11}$$

The angular displacements *θ<sup>x</sup>* and *θ<sup>y</sup>* are assumed to be small. So, they are approximated by the derivative of *u* and *v* with respect to the *z*-direction. Therefore, using Equations (10) and (11), *θ<sup>x</sup>* and *θ<sup>y</sup>* have been expressed as:

$$\begin{array}{l} \theta\_x(z.t) = -\frac{\partial u}{\partial z} = -\frac{df(z)}{dz}q\_1 = -\operatorname{g}(z)q\_1\\ \theta\_y(z.t) = \frac{\partial v}{\partial z} = \frac{df(z)}{dz}q\_2 = \operatorname{g}(z)q\_2. \end{array} \tag{12}$$

To derive the equations of motion, the kinetic energy and the potential energy are specified by generalized coordinates. So, using Equation (3), the kinetic energy of the disk in generalized coordinates is expressed as:

$$T\_{\rm disk} = \frac{1}{2}Mf^2(L)\left(\dot{\vec{q}}\_1^2 + \dot{\vec{q}}\_2^2\right) + \frac{1}{2}I\_t\chi^2(L)\left(\dot{\vec{q}}\_1^2 + \dot{\vec{q}}\_2^2\right) + \frac{1}{2}I\_p\Omega^2 + I\_p\Omega\chi^2(L)\left(\dot{\vec{q}}\_1\eta\_2\right) \tag{13}$$

where

$$g(z) = \frac{\partial f(z)}{\partial z} = \frac{n\,\,\pi}{(L\_2 - L\_1)} \cos\left(\frac{n\,\,\pi}{(L\_2 - L\_1)} z - L\_1\right). \tag{14}$$

In addition, Equation (4) in the generalized equation becomes:

$$T\_{shaft} = \frac{1}{2}\rho A \int\_0^L f^2(z)dz \left(\dot{\vec{q}}\_1^2 + \dot{\vec{q}}\_2^2\right) + \frac{1}{2}I\_l \int\_0^L \mathbf{g}^2(z)dz \left(\dot{\vec{q}}\_1^2 + \dot{\vec{q}}\_2^2\right) + \frac{1}{2}I\_p\Omega^2 + I\_p\Omega \int\_0^L \mathbf{g}^2(z)dz \left(\dot{\vec{q}}\_1\eta\_2\right). \tag{15}$$

Likewise, using the generalized coordinates, the potential energy in the generalized coordinate will be:

$$\begin{array}{ll} \mathrm{LI} = \frac{1}{2} \mathrm{EI} \int\_{0}^{L} h^{2}(z) dz \left( \dot{q}\_{1}^{2} + \dot{q}\_{2}^{2} \right) + \frac{1}{2} \mathrm{K} \mathrm{GA} \int\_{0}^{L} \left( \beta\_{x}^{2} + \beta\_{y}^{2} \right) dz \, + \left( -\mathrm{F}\_{1x} \, q\_{1} - \mathrm{F}\_{1y} \, q\_{2} \right) \delta(z - L\_{1}) \\ \qquad + \left( -\mathrm{F}\_{2x} \, q\_{1} - \mathrm{F}\_{2y} \, q\_{2} \right) \delta(z - L\_{2}) \end{array} \tag{16}$$

where *δ* represnts the dirac delta function and

$$h(z) = \frac{\partial f(z)}{\partial z} = -\left(\frac{n\ \pi}{(L\_2 - L\_1)}\right)^2 \sin\left(\frac{n\ \pi}{(L\_2 - L\_1)} z - L\_1\right). \tag{17}$$

In addition, substituting the displacement function into the kinetic energy of the mass unbalance expression of Equation (11) gives:

$$T\_u \cong m\_u d\Omega \; f(L) \left( q\_1 \cos \Omega t - q\_2 \sin \Omega t \right). \tag{18}$$

#### *Ball Bearing Model*

In this section, a mathematical model for calculating bearing stiffness is proposed by analyzing the equations in the bearing dynamic model, which is based on Jones and Harris's efforts [14]. This mathematical model aims to give a comprehensive consideration of the nonlinear stiffness of the ball bearing, and it can be seen that the bearing stiffness is critically dependent on the preloading, *g*, of the rolling elements. In case of rolling element bearings, the elastic deformation takes place at both the inner raceway and the outer raceway with the rolling element. Based on the Hertzian contact theory, the relation between load *F* and deflection is [15]:

$$F = k\_p \delta\_p^{3/2} \tag{19}$$

where *δ<sup>p</sup>* is the contact deformation and *kp* is a load–deformation constant for a single point contact (either at the inner or outer raceway). If the ball and raceway are made of steel, then [16,17]:

$$k\_p = 2.15 \times 10^5 (\sum \rho)^{-\frac{1}{2}} (\delta^\*)^{-3/2} \qquad \frac{N}{mm^{1.5}} \tag{20}$$

where *δ*∗ is the dimensionless contact deformation and ∑ *ρ* is the curvature sum [14]. So, one can write the stiffness of the inner and outer ring contact as:

$$\begin{array}{llll}k\_{\rm pi} = 2.15 \times 10^5 (\sum \rho\_i)^{-\frac{1}{2}} (\delta^\*\_{\rm i})^{-3/2} & \frac{N}{mm^{1.5}}\\k\_{\rm po} = 2.15 \times 10^5 (\sum \rho\_o)^{-\frac{1}{2}} (\delta^\*\_{\rm o})^{-3/2} & \frac{N}{mm^{1.5}}\end{array} \tag{21}$$

where

$$\begin{array}{c} \Sigma \rho\_i = \frac{1}{D\_b} \left( 4 - \frac{1}{f\_i} + \frac{2\gamma}{1 - \gamma} \right) \\ \Sigma \rho\_o = \frac{1}{D\_b} \left( 4 - \frac{1}{f\_o} + \frac{2\gamma}{1 - \gamma} \right) \end{array}$$

and

$$f\_i = \frac{r\_i}{D\_b} \quad f\_o = \frac{r\_o}{D\_b} \quad \gamma = \frac{D}{D\_m} \quad D\_m = \frac{1}{2}(d\_i + d\_o) \cong \frac{1}{2}(D+d).$$

In the above equations, subscripts *i* and *o* represent inner and outer raceways and *Db* is the ball diameter. In addition, *Dm* is the pitch diameter, *di* and *do* are the inner and outer ring raceway contact diameter, and *ri* and *ro* are the inner and outer raceway groove radius, respectively.

The total deformation, *δ*, at a single rolling element location is given by:

$$
\delta = \delta\_{pi} + \delta\_{po} = (\frac{1}{k\_{pi}} + \frac{1}{k\_{po}})F^{\frac{2}{3}}.\tag{22}
$$

Equation (19) can be rewritten as:

$$F = k\_{pio} \delta^{\frac{3}{2}} \tag{23}$$

where

$$k\_{pio} = \frac{k\_{pi}k\_{po}}{(k\_{pi}^{\frac{2}{3}} + k\_{po}^{\frac{2}{3}})} - $$

*kpio* is the load–deformation constant for two point contacts of a ball with raceways.

To calculate the deformations, the conventions shown in Figure 3 are used. Considering the presence of radial clearance, the force produced by every spring in the horizontal and vertical direction can be obtained [18]:

$$\begin{aligned} F\_{\mathbf{x}} &= k\_{pio} (\mathcal{g} + \iota \cos \psi\_i + \upsilon \sin \psi\_i)^{\frac{3}{2}} \cos \psi\_i \\ F\_{\mathbf{y}} &= k\_{pio} (\mathcal{g} + \upsilon \cos \psi\_i + \upsilon \sin \psi\_i)^{\frac{3}{2}} \cos \psi\_i \end{aligned} \tag{24}$$

where *g* is the radial preload between the ball and races, *u* and *v* are the displacements of the moving ring in the *x* and *y* directions, respectively, and *ψ<sup>i</sup>* is the angular position of the *i* th element.

$$\mathbf{u}\cos\psi\_{\ell} + \mathbf{v}\sin\psi\_{\ell}$$

**Figure 3.** Bearing model.

Finally, the bearing stiffness can be simplified to [15]:

$$\begin{aligned} k(u) = k\_{pio} \sum\_{i=1}^{z} \left\{ g + \iota \cos \psi\_i \right\}^{\frac{1}{2}} \{ \cos \psi\_i - \frac{B}{15A} \sin \psi\_i \} \cos \psi\_i \\\ A = \sum\_{i=1}^{z} \left( g + \iota \cos \psi\_i \right)^{\frac{1}{2}} \sin^2 \psi\_i \\\ B = \sum\_{i=1}^{z} \left( g + \iota \cos \psi\_i \right)^{1/2} \sin \psi\_i \cos \psi\_i \\\ \psi\_i = \frac{\pi}{Z} (2i - 1) \qquad i = 1, 2, 3, \dots, z \end{aligned} \tag{25}$$

where *Z* is the number of rolling elements.

It can be seen that the bearing stiffness is critically dependent on value of the preloading *g* of the rolling elements.

#### **3. Numerical Solution**

The model is applied for an installed pump (P-502B) in a Kermanshah oil refinery in Iran. This centrifugal pump and its impeller are shown in Figure 4.

**Figure 4.** Rotor pump P-502B in a Kermanshah oil refinery.

The schematic of the shaft and impeller of the above pump is shown in Figure 5. The flowing fluid is light petroleum gases including propane and butane.

**Figure 5.** The schematic diagram of the pump rotor (P-502B).

The shaft of the mentioned pump is supported by two bearings SKF7308 and SKF6307, which have been installed in distance *L*<sup>1</sup> = 0.125 m and *L*<sup>2</sup> = 0.521 m respectively, see Figure 1a. In addition, three impellers with different diameters have been installed at the end of the shaft.

The characteristics of shaft and impeller are shown in Tables 1 and 2, respectively.

**Table 1.** Characteristics of shaft P-502B.



**Table 2.** Characteristics of impeller P-502B.

#### *Bearing Stiffness of P-502B*

To obtain the stiffness of the bearing, knowledge of the properties is necessary. The general properties of these ball bearings are shown in Table 3.


**Table 3.** Dimensional specifications and extracted parameters in bearings.

It can be seen that the bearing stiffness is critically dependent on the value of preloading force *g* for the rolling elements.

Substituting the above values in Equation (25), the values of stiffness in the *x* and *y* directions have been obtained as:

Bearing stiffness for SKF6307:

$$k(\mathbf{x}) = k(y) = \left(1.101 - 94.6 \mathbf{x}^2\right) \times 10^5 \quad \frac{N}{m}.\tag{26}$$

Bearing stiffness for SKF7308:

$$k(\mathbf{x}) = k(y) = \left(1.13 - 97.8\mathbf{x}^2\right) \times 10^5 \quad \frac{N}{m}.\tag{27}$$

Finally, using the Equations (3), (4), (6), (7) and Lagrangian Equation (1), the mathematical equation for the case study can be expressed as:

$$
\begin{bmatrix} 10.86 & 0 \\ 0 & 10.86 \end{bmatrix} \begin{bmatrix} \dot{q}\_1 \\ \dot{q}\_2 \end{bmatrix} + \begin{bmatrix} 0 & 0.902\Omega \\ -0.902\Omega & 0 \end{bmatrix} \begin{bmatrix} \dot{q}\_1 \\ \dot{q}\_2 \end{bmatrix} + $$

$$
\begin{bmatrix} 2.44 \times 10^7 + 125 \times 10^5 q\_1^2 & 0 \\ 0 & 2.44 \times 10^7 + 125 \times 10^5 q\_2^2 \end{bmatrix} \begin{bmatrix} q\_1 \\ q\_2 \end{bmatrix} = \tag{28}
$$

$$
\begin{bmatrix} 28.5 \times 10^{-5} \sin \Omega t & 0 \\ 0 & 28.5 \times 10^{-5} \cos \Omega t \end{bmatrix} . \tag{29}
$$

The above equation is a lumped parameter model that can be described as:

$$
\tilde{M}\ddot{q} - \Omega \tilde{D}\dot{q} + \tilde{K}q = \tilde{\mathbb{Q}}^{nc}.\tag{29}
$$

Looking at Equation (28), one can see that the damping matrix is dependent on the angular velocity of the shaft. The stiffness matrix is implicitly dependent on the amount of amplitude of vibration and so, the equations of motion are nonlinear equations. This will lead to a complicated model. So, for simplicity, normal clearance has been considered to make the value of the stiffness coefficient independent from the deflection, and so, the equations of motion will become linear ordinary differential equations. According to the SKF General catalog, hence *g* = 0.051 and *g* = 0.053 used for SKF6307 and SKF7308 respectively, the equivalent stiffness for bearings in a normal clearance case have been derived as follows:

$$k(r) = \textbf{1.040} \times \textbf{10^5} \,\frac{\text{N}}{\text{m}} \text{ for SKF6307}$$

$$k(r) = \textbf{1.066} \times \textbf{10^5} \,\frac{\text{N}}{\text{m}} \text{ for SKF7308}$$

To obtain the critical speed of the rotor, the eigenvalues and thus the natural frequencies of the system must be determined. To do this, a slightly different but equivalent approach can be used to combine the pair in Equation (29). By letting *α* = *q*<sup>2</sup> − *jq*<sup>1</sup> and subtracting *j*× in the second Equation (29) from the first Equation (29), we have [19]:

$$
\tilde{M}\ddot{\alpha} - j\Omega \tilde{D}\dot{\alpha} + \tilde{K}\alpha = 0. \tag{30}
$$

Now, we can solve this equation by letting *α* = *α*0*ej*Ω*<sup>t</sup>* and taking the positive root to obtain backward critical speed. If we let *α* = *α*0*e*−*j*Ω*<sup>t</sup>* in which Ω is the forward critical speed, one must add the out of balance moments to the system to determine the response. Thus:

$$
\delta \check{M}\ddot{\alpha} - j\Omega \check{D}\dot{\alpha} + \check{K}\alpha = F\_{\mu}.\tag{31}
$$

Likewise, by *α* = *q*<sup>2</sup> − *jq*<sup>1</sup> and subtracting *j*× in the second equation from the first, the rotor amplitude is calculated [19]. In addition, to determine the effect of impeller diameter, the various values of this parameter have been considered. To do this procedure, a code in MATLAB software has been written, and the following outputs shown in Figures 6–8 have been derived.

The amplitudes of the rotor response plots for all pump impellers that are recommended by pump manufacture are shown in Figures 6–8. The summary of the results is shown in Table 4.

**Figure 6.** Calculated response at impeller diameter = 0.254 m for API unbalance, normal bearing clearance.

**Figure 7.** Calculated response at impeller diameter = 0.281 m for API unbalance, normal bearing clearance.

**Figure 8.** Calculated response at impeller diameter = 0.31 m for API unbalance, normal bearing clearance.


**Table 4.** Effect of impeller diameter on the amplitude of vibration and critical speed of the rotor.

So, one can see that by increasing the impeller diameter, the amplitude of vibration at critical speed will be increased up to 52%, while the value of critical speed will decrease. The value of amplitude of vibration shown in this table has been obtained without considering the structural damping of the shaft, bearings, and impeller and also the effects of fluid flow in the pump. So, in reality, these values are much smaller than the tabulated values in Table 4.

## **4. Conclusions**

In this paper, the equations of motion for the shaft of a centrifugal pump have been derived using Lagrangian equations and based on Timoshenko beam theory and Jones and Harris [20,21] modeling for bearings based on Hertzian theory for rolling element [22,23]. By these equations, it is possible to define a lumped parameter model for the system and determine the amplitude of rotor in critical speeds [12,24]. The solutions of these equations for an applied case study show that as the impeller diameter increases, the amplitude of vibration at critical speeds increases, too. In this case study, by replacing the impeller diameter in the equations with a larger size (impeller diameter 0.256 m replaced with 0.32 m), it is shown that the amplitude is approximately increased by 52% in the critical speed. So, as the impeller trimming or impeller exchange is a general method for the maintenance in plants, it must be considered that the bigger diameter will yield more amplitude of vibration at critical speed.

It must be considered that modeling of the shaft, bearings, and impellers with details and assembling them and then conducting finite element analysis is very complex, and as the case study is a process pump in an oil refinery that works 24 h a day, stopping the production line and then measuring, modeling, and conducting FEM analysis is very expensive. On the other hand, the purpose of this manuscript is to introduce an analytical model for all the essential rotary parts of a centrifugal pump in detail and then conduct a sensitive analysis for the effect of impeller diameter on the critical speed of it inversely. For this purpose, the dynamic inverse solution is very suitable when there are practical limitations for FEM analysis, such as it being very complex and expensive.

**Author Contributions:** Investigation, A.S., M.K. and M.S.; Methodology, A.S. and S.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

