Dynamics and Vibration Analysis of a Rotor-Bearing System in a Turbofan Engine with Emphasis on Bearings Modeling

Hossein Rahmani^*, Mohammadreza Elhami and Amin Moslemi Petrudi

Department of Mechanical Engineering, IHU University, Tehran, Iran
E-mail: hsnrahmani@ihu.ac.ir
*Corresponding Author

Received 17 August 2021; Accepted 25 July 2022; Publication 20 August 2022

Abstract

In this research, while investigating the vibration analysis of rotary axes, we specifically investigate the rotor of a turbofan engine used in the industry. The features of this rotor range are high-performance, lightweight, and low-vibration range. These three factors are in contradiction with each other, resulting in a thorough examination of the total vibration of the complete turbofan rotor. To achieve this, various parameters such as the concentration of properties, rotational inertia, gyroscopic torque, rotational loading, effects of unbalanced mass, crevillous effects, bearing flexibility, etc. have been studied in modeling. The rotor’s natural frequencies, along with the critical velocity, are plotted as well as the shape of its modes. The software is required to perform the computations written by Ansys software and after ensuring its accuracy.

Keywords: Vibration, shaft modeling, turbofan engine, natural frequency, mood shape.

1 Introduction

Rotordynamics is a field under mechanics, that mainly deals with the vibration of rotating structures. In recent days, the study of rotor dynamics has gained more importance within Jet engine industries. The main reason is Jet engine consists of many rotating parts that constitutes a complex dynamic system. While designing rotors of high-speed turbomachinery, it is of prime importance to consider rotor dynamics characteristics to account. Considering these characteristics at the design phase may prevent the jet engine from severe catastrophic failures. These rotordynamic characteristics can be determined with the help of the much relied Finite elements method. Traditionally, Rotordynamic analyses were performed with specialized commercial tools. On the other hand capabilities of more general FEA software have gradually been developed and the most commonly used software is Ansys. The aim of this thesis work is to build an RM12 Jet engine rotor model in Ansys and evaluate its rotor dynamics capabilities with the specialized rotor dynamics tools. These works help in understanding, modeling, simulation, and post-process techniques for the rotor dynamics analysis of RM12 Jet engine rotor using Ansys. The science of vibrations studies the oscillating motion of objects and the forces associated with them. The subject of vibrations is one of the most important and practical topics in the field of engineering. All objects that have the properties of elasticity and mass are capable of vibration, which causes an error and wear in all machines and structures. Rotating shafts and their vibrations have played a very important role in the industry because the shafts are the main means of power generation and power transmission in the industry. Axes as the main part of internal combustion engines, types of steam, water, and gas turbines, generators, electric motors, pumps, fans, compressors, and all kinds of machine tools such as lathes, drills, and milling. Also have many applications in the automotive, rail, shipping, and aviation industries. Creating severe vibrations and applying large forces in a short time or creating small and light vibrations and applying small forces in a long time can cause rotation systems to fail. This failure can manifest itself in any of the components of the rotating system such as shafts, and bearings [1]. Investigated the dynamic behavior of rotors, they obtained equations for the elements of a rotor such as disks, the main axis of the rotor, and bearings using the Rayleigh-Ritz method for a simple rotor [2]. Conducted research on the critical speed of shafts and developed his famous method for determining natural frequency [3]. Analyzed the flexible rotor system with magnetic bearings and obtained the specific values, the optimal dynamic properties for the rotor, which make it stable at different speeds, as well as the critical speed of the rotor, using Campbell diagrams. They examined the shape of the modes and the frequency and time of the responses [4]. Performed various effects such as shear deformation, rotational inertia, gyroscopic torque, axial force, torsional torque, axial force, gyroscopic torque, and nonlinear factors in the shaft model. Other axial components such as vanes, seals, dampers, foundations, and enclosures were also gradually incorporated into the axle models. Thus, today, quite complex models are created for the axes that make it possible to accurately analyze the vibrations of the axis. The issue at hand was initially limited to finding the critical axis and finding the amplitude of the fluctuations caused by the nomenclature. The vibration of the shaft due to the force of gravity and the misalignment of the bearings are other issues related to the shafts. Various factors of instability in the axes such as internal shaft damping, bearing oil layer, shaft asymmetry, alternating axial force, and the coupling of transverse and torsional vibrations were identified and investigated. In the analysis of transverse axis vibrations [5], first used the finite element method. They considered the shaft to be composed of shafts, bearings, and thin rigid discs. The shaft element is based on the Euler-Bernoulli beam element previously acquired by Archer in 1963 and is intersected by two transverse vibrations considering the transverse vibrations. The bearings were modeled as a parallel spring and damper with 16 coefficients corresponding to different directions, and the gyroscopic effect was also considered in rigid disk modeling. In this article, they calculated the threshold of instability and the response to engagement [6]. Studied the nonlinear dynamic behavior of the rotor with roller bearings numerically and then compared it with experimental data [7]. Used the finite element method in numerical analysis. Examining the position, he modeled the rotor using beam elements [8]. Compared the dynamic behavior of a real model of a rotor with a disk and the blades with a simplified model of a rotor with a disk. They examined the complete diagrams and frequencies of both models using Ansys software Gave [9]. The effect of bearings on both rotors on the dynamic behavior of the GTF transmission system has been explored and quantified, especially for the influence of the location of bearings on both rotors on the load sharing behavior of the star gearing system in the GTF gearbox. The results show that the position of bearings on both rotors could affect the load sharing performance and maximum floating amounts of members in the star gearing system of the GTF gearbox, therefore, the location of bearings on rotors can be rearranged to achieve better performance in vibration control and load distribution of the system [10]. Numerical study on the rotor dynamic analysis of a dual-spool turbofan engine in the context of blade defect events. The simulation results indicate that the high-vibration status of the rotor commences beyond 10,000 rpm, which is identified as the first critical speed of the lower speed rotor. Moreover, we monitored the unbalanced stages near the inter-rotor bearing, which prominently influences the overall rotor dynamic status, and the corrosion of the HPC to prevent further instability [11]. A second bearing is positioned aft of the gear reduction and supports the gear reduction. The second bearing is a thrust bearing. A fan drive turbine drive shaft drives the gear reduction. The fan drive turbine drive shaft has a weakened link which is aft of the second bearing such that the fan drive turbine drive shaft will tend to fail at the weakened link and a location aft of the second bearing [12]. In this study, for a turbofan engine-powered unmanned aerial vehicle, Global Warming Potential (GWP) calculations have been performed covering a range of flight Mach number and altitude. Results were presented to illustrate the GWP magnitudes for different flight speeds and altitudes in a comparative manner for the turbofan-powered UAV. Optimization calculations were done by using a genetic algorithm. a genetic algorithm (GA) is a metaheuristic inspired by the process of natural selection that belongs to the larger class of evolutionary algorithms (EA). Genetic algorithms are commonly used to generate high-quality solutions to optimization and search problems by relying on biologically inspired operators such as mutation, crossover and selection.

2 Statement of the Problem

2.1 Vibrations of Rotating Axes

Although the complete set of the rotor consists of many components, which play a key role in the vibration analysis of a rotor, there are 4 main components, which are: shaft, components mounted on the shaft, bearing, base. Figure 1 shows the complete set of the rotor.

Figure 1 Schematics of rotor assembly.

In this research, using finite element analysis, each of these four members is modeled and finally, the matrix related to each component is obtained, which is used to analyze the whole rotor set, the total mass matrix of the system, and the stiffness matrix of the whole system. And the damping matrix of the whole system is assembled and placed in the following vibration equation, which will give the last equation after solving the system frequency response [13]. In the following equation, $\mathrm{\Omega }$ is the rotational speed of the rotor.

$$[M]\{\ddot{\delta }\}+(\mathrm{\Omega }[G]+[C])\{\dot{\delta }\}+[K]\{\delta \}=\{f\}$$

(1)

2.2 Kinematics of Rotating Elements

A rotating element is considered as a rigid body and its rotational kinetic energy is calculated as follows:

$$T_{rot}=\frac{1}{2}(I_{xx}{\omega }_x^2+I_{yy}{\omega }_y^2+I_{zz}{\omega }_z^2)$$

(2)

In general, the position of a rigid object in space by three angles $\alpha$, $\beta$, $\mathrm{\Phi }$ to Euler angles and a fixed coordinate system (bare) with fixed unit vectors and another with rotating unit vectors is a coordinate system on the center of a fixed rotating element Has been. When all of the Euler angles are zero, the two coordinate systems overlap. For a general position for the element in space, the connection between the two devices can be successfully and completely calculated in three steps with the mentioned Euler angles. Figure 2 shows the kinematics of rotating elements [14].

Figure 2 Kinematics of rotating elements (Cao, 2018).

2.3 Rail Shaft Element

Strain energy for a small disk located at a distance s along the element, the flexural strain energy relationship is as follows [15]:

$$d{u}_B^e=\frac{1}{2}{\left\{\begin{array}{c}\hfill v^{\prime }\hfill \\ \hfill w^{\prime \prime }\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill EI\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill EI\hfill \end{array}\right]\left\{\begin{array}{c}\hfill v^{\prime }\hfill \\ \hfill w^{\prime \prime }\hfill \end{array}\right\}ds$$

(3)

Which is EI flexural axis. If the element is also under a constant axial force, the strain energy stored in the disk will be called the strain energy under a constant axial load as follows:

$$d{u}_A^e=-\frac{1}{2}{\left\{\begin{array}{c}\hfill v^{\prime }\hfill \\ \hfill w^{\prime \prime }\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill P\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill P\hfill \end{array}\right]\left\{\begin{array}{c}\hfill v^{\prime }\hfill \\ \hfill w^{\prime \prime }\hfill \end{array}\right\}ds$$

(4)

Strain energies for the whole element are obtained by integrating differential expressions over the entire length of the element:

$U_A^e$	$={\displaystyle \frac{EI}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right]\mathit{d}y$	(5)
$U_A^e$	$={\displaystyle \frac{F}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2\right]\mathit{d}y$	(6)

The total energy of the axis is obtained from Equation (7):

$U_A^e$	$={\displaystyle \frac{EI}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right]\mathit{d}y$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{F}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2\right]\mathit{d}y$	(7)

The kinetic energy formula of an element uses the expansion of the disk kinetic energy equation located at a distance it is obtained from the length of the element [15].

$d{T}^e$	$={\displaystyle \frac{1}{2}}{\left\{\begin{array}{c}\hfill \dot{v}\hfill \\ \hfill \ddot{w}\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill \mu \mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mu \hfill \end{array}\right]\left\{\begin{array}{c}\hfill v^{\prime }\hfill \\ \hfill w^{\prime \prime }\hfill \end{array}\right\}ds+{\displaystyle \frac{1}{2}}\dot{\mathrm{\Phi }}{\tau }_{p}ds$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{2}}{\left\{\begin{array}{c}\hfill \dot{\beta }\hfill \\ \hfill \dot{\gamma }\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill I_d\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_d\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \dot{\beta }\hfill \\ \hfill \dot{\gamma }\hfill \end{array}\right\}ds-\dot{\mathrm{\Phi }}\dot{\gamma }\beta {\tau }_{p}ds$	(8)
$T_s$	$={\displaystyle \frac{\rho s}{2}}{\displaystyle {\int }_0^L}({\dot{u}}^2+{\dot{w}}^2)\mathit{d}y+{\displaystyle \frac{\rho l}{2}}{\displaystyle {\int }_0^L}({\dot{\beta }}^2+{\dot{\gamma }}^2)\mathit{d}y$
	$\mathrm{\hspace{1em}}+\rho IL{\mathrm{\Omega }}^2-2\rho I\mathrm{\Omega }{\displaystyle {\int }_0^L}\gamma \beta \mathit{d}y$	(9)

2.4 Element Shaft Euler-Bernoulli

If the effect of rotational inertia on the rail shaft element is ignored, the resulting model becomes the Euler-Bernoulli model. Rotational mass matrices and gyroscope matrices are removed from the model [16]. The relationship between the strain energy of an axis and the symmetric surface is in the form of Equation (15).

$$U_1=\frac{EI}{2}{\int }_0^L\left[{\left(\frac{{\partial }^{2}u}{\partial y^2}\right)}^2+{\left(\frac{{\partial }^{2}w}{\partial y^2}\right)}^2\right]\mathit{d}y$$

(10)

Which E is Young’s modulus. If the axis is under a constant axial force F, a second expression is added to the axis strain energy, which is obtained as follows:

$$U_2=\frac{F}{2}{\int }_0^L\left[{\left(\frac{\partial u}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial y}\right)}^2\right]\mathit{d}y$$

(11)

Therefore, the strain energy of the whole axis is obtained from the following equation:

$U_T$	$={\displaystyle \frac{F}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2\right]\mathit{d}y$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{EI}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right]\mathit{d}y$	(12)

The general formula for the axial kinetic energy is obtained by extending the disk kinetic energy equation. For an element, the length of the L element, the expression energy is as follows.

$$T_S=\frac{\rho s}{2}{\int }_0^L({\dot{u}}^2+{\dot{w}}^2)\mathit{d}y$$

(13)

2.5 Element Shaft Timoshenko

The general formula for the kinetic energy of an element is obtained by extending the kinetic energy equation of a disk located at a distance S from the length of the element [15].

$d{T}^e$	$={\displaystyle \frac{1}{2}}{\left\{\begin{array}{c}\hfill \dot{v}\hfill \\ \hfill \dot{w}\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill \mu \mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mu \hfill \end{array}\right]\left\{\begin{array}{c}\hfill \dot{v}\hfill \\ \hfill \dot{w}\hfill \end{array}\right\}ds+{\displaystyle \frac{1}{2}}\dot{\mathrm{\Phi }}{\tau }_{p}ds$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{2}}{\left\{\begin{array}{c}\hfill \dot{\beta }\hfill \\ \hfill \dot{w}\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill I_D\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_D\hfill \end{array}\right]\left\{\begin{array}{c}\hfill \dot{\beta }\hfill \\ \hfill \dot{w}\hfill \end{array}\right\}ds-\dot{\varphi }\dot{\gamma }\beta {\tau }_{p}ds$	(14)

For an element of length L, the expression of the total kinetic energy is as follows:

$T_s$	$={\displaystyle \frac{\rho s}{2}}{\displaystyle {\int }_0^L}({\dot{u}}^2+{\dot{w}}^2)\mathit{d}y+{\displaystyle \frac{\rho I}{2}}{\displaystyle {\int }_0^L}({\dot{\beta }}^2+{\dot{\gamma }}^2)\mathit{d}y$
	$\mathrm{\hspace{1em}}+\rho IL{\mathrm{\Omega }}^2-2\rho I\mathrm{\Omega }{\displaystyle {\int }_0^L}\dot{\gamma }\beta \mathit{d}y$	(15)

For a small disk located at a distance s along with the element, the flexural strain energy relationship is as follows:

$d{u}_B^e={\displaystyle \frac{1}{2}}{\left\{\begin{array}{c}\hfill v^{\prime \prime }\hfill \\ \hfill w^{\prime \prime }\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill EI\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill EI\hfill \end{array}\right]\left\{\begin{array}{c}\hfill v^{\prime \prime }\hfill \\ \hfill w^{\prime \prime }\hfill \end{array}\right\}ds$

(16)

Which is axis-bending stiffness. If the element is also subject to a constant axial force P, the strain energy stored in the disk will be called the strain energy under a constant axial load as follows:

$d{u}_A^e=-{\displaystyle \frac{1}{2}}{\left\{\begin{array}{c}\hfill v^{\prime }\hfill \\ \hfill w^{\prime }\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill P\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill P\hfill \end{array}\right]\left\{\begin{array}{c}\hfill v^{\prime }\hfill \\ \hfill w^{\prime }\hfill \end{array}\right\}ds$

(17)

Due to the shear deformation, the strain energy resulting from the shear must also be added to the equations:

$d{u}_s^e={\displaystyle \frac{1}{2}}{\left\{\begin{array}{c}\hfill v_{\mathit{shear}}^{\prime \prime }\hfill \\ \hfill w_{\mathit{shear}}^{\prime \prime }\hfill \end{array}\right\}}^T\left[\begin{array}{cc}\hfill \mathit{kAG}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathit{kAG}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill v_{\mathit{shear}}^{\prime \prime }\hfill \\ \hfill w_{\mathit{shear}}^{\prime \prime }\hfill \end{array}\right\}ds$

(18)

Strain energies for the whole element are obtained by integrating differential expressions over the entire length of the element:

$U_s^e$	$={\displaystyle \frac{kAG}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{{\partial }^2{u}_s}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right]\mathit{d}y$
$U_B^e$	$={\displaystyle \frac{EI}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right]\mathit{d}y$
$U_A^e$	$={\displaystyle \frac{F}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2\right]\mathit{d}y$	(19)

Therefore, the strain energy of the whole element is obtained from the following equation:

$U_T^e$	$={\displaystyle \frac{F}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^2+{\left({\displaystyle \frac{\partial w}{\partial y}}\right)}^2\right]\mathit{d}y$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{EI}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right]\mathit{d}y$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{kAG}{2}}{\displaystyle {\int }_0^L}\left[{\left({\displaystyle \frac{{\partial }^2{u}_s}{\partial y^2}}\right)}^2+{\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)}^2\right]\mathit{d}y$	(20)

For the element of Timoshenko’s beam the functions of the Fig are as follows [15].

The disk is structured in such a way that the disk is assumed to be rigid, so there will be no hard matrix for the disk. We consider the disk to be a rigid body and define kinetic energy for it [20]. The term kinetic energy is defined for a disk that has 4 degrees of freedom of rotation and rotation as follows:

$T_D={\displaystyle \frac{1}{2}}{M}_D({\dot{u}}^2+{\dot{w}}^2)+{\displaystyle \frac{1}{2}}(I_{\mathit{DX}}{w}_X^2+I_{\mathit{DY}}{w}_Y^2+I_{\mathit{DZ}}{w}_Z^2)$

(21)

That u and w, are the coordinates of the center of the disk at R and $I_{\mathit{DX}}$, $I_{\mathit{DY}}$, $I_{\mathit{DZ}}$, The moment of inertia of the disk concerning the x, y, and z axes. Assuming the disk is symmetrical ($I_{\mathit{DZ}}=I_{\mathit{DZ}}$) And small angles $\theta$ and $\mathrm{\Psi }$ constant angular velocity $\dot{\varphi }=\mathrm{\Omega }$ the kinetic energy equation will be:

$T_D$	$={\displaystyle \frac{1}{2}}{M}_D({\dot{u}}^2+{\dot{w}}^2)+{\displaystyle \frac{1}{2}}{I}_{DX}({\dot{\theta }}^2+{\dot{\mathrm{\Psi }}}^2)$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{2}}{I}_{DY}({\mathrm{\Omega }}^2+2\mathrm{\Omega }\dot{\mathrm{\Psi }}\theta )$	(22)

If the displacement vector is considered as follows:

$$\delta ={[u,w,\theta ,\mathrm{\Psi }]}^T$$

(23)

Using Lagrangian equations:

The first matrix is the classical mass matrix and the second matrix is the gyroscopic matrix.

${\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial T}{\partial \dot{\delta }}}\right)-{\displaystyle \frac{\partial T}{\partial \delta }}$	$=\left[\begin{array}{cccc}\hfill M_D\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill M_D\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_{DX}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_{DX}\hfill \end{array}\right]\left[\begin{array}{c}\hfill \ddot{u}\hfill \\ \hfill \ddot{w}\hfill \\ \hfill \ddot{\theta }\hfill \\ \hfill \ddot{\psi }\hfill \end{array}\right]$
	$\mathrm{\hspace{1em}}-\mathrm{\Omega }\left[\begin{array}{cccc}\hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_{Dy}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -I_{Dy}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill \dot{u}\hfill \\ \hfill \dot{w}\hfill \\ \hfill \dot{\theta }\hfill \\ \hfill \dot{\psi }\hfill \end{array}\right]$	(24)

Where $\rho$ is the amount of mass per unit volume. The cross-sectional area of the beam is assumed to be constant and the moment of inertia of the cross-sectional area of the axis relative to its diameter is assumed to be constant. Therefore, the mass and hard disk element matrices are as follows:

$[M_D^e]=\left[\begin{array}{cccc}\hfill M_D\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill M_D\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_{DX}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_{DX}\hfill \end{array}\right]\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}[G_d^e]=\left[\begin{array}{cccc}\hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill I_{Dy}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -I_{Dy}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]$

(25)

In the above equations, $M_D$ is the mass of the disk and $I_{\mathit{DX}}$ is the moment of diagonal inertia and $I_{\mathit{Dy}}$ is the moment of polar inertia.

In the case of bearings, the element is such that the bearing is on one side of the node and the other node of the bearing element is connected to the base element because each node has four degrees of freedom, so the element has eight degrees of freedom and therefore the matrix Its stiffness and damping will be 8 $\times$ 8, and using the short bearing assumption, we will have a stiffness and damping matrix [21].

${[k_b]}_e=\left[\begin{array}{cccccccc}\hfill K_{xx}\mathit{\hspace{1em}}\hfill & \hfill K_{xy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -K_{xx}\mathit{\hspace{1em}}\hfill & \hfill -K_{xy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill K_{yx}\mathit{\hspace{1em}}\hfill & \hfill K_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -K_{yx}\mathit{\hspace{1em}}\hfill & \hfill -K_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill -K_{xx}\mathit{\hspace{1em}}\hfill & \hfill -K_{xy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill K_{xx}\mathit{\hspace{1em}}\hfill & \hfill K_{yx}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill -K_{yx}\mathit{\hspace{1em}}\hfill & \hfill -K_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill K_{yx}\mathit{\hspace{1em}}\hfill & \hfill K_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]$

(26)

The matrix elements are the hardness coefficients of the oil film. The first index indicates the direction of the force and the second index indicates the direction of displacement.

${[C_b]}_e=\left[\begin{array}{cccccccc}\hfill C_{xx}\mathit{\hspace{1em}}\hfill & \hfill C_{xy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -C_{xx}\mathit{\hspace{1em}}\hfill & \hfill -C_{xy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill C_{yx}\mathit{\hspace{1em}}\hfill & \hfill C_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -C_{yx}\mathit{\hspace{1em}}\hfill & \hfill -C_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill -C_{xx}\mathit{\hspace{1em}}\hfill & \hfill -C_{xy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill C_{xx}\mathit{\hspace{1em}}\hfill & \hfill C_{yx}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill -C_{yx}\mathit{\hspace{1em}}\hfill & \hfill -C_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill C_{yx}\mathit{\hspace{1em}}\hfill & \hfill C_{yy}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]$

(27)

The matrix elements are the damping coefficients of the oil film.

Assume that the terms viscous damping and stiffness are known and the effects of bending can be ignored. The virtual work of the forces entering the axis is written as follows:

$\delta w$	$=-k_{xx}u\delta u-k_{xz}w\delta u-k_{zz}w\delta w-k_{zx}u\delta w-c_{xx}\dot{u}\delta u-c_{xz}\dot{w}\delta u$
	$\mathrm{\hspace{1em}}-c_{zz}w\delta w-c_{zx}\dot{u}\delta w$
$\delta w$	$=F_u\delta u+F_w\delta w$	(28)

Where $F_u$ and $F_w$ are the general components of force. If we want to represent the equations in the form of a matrix, we will have the assumption $F_{\theta }=F_{\mathrm{\Psi }}=0$.

$\left[\begin{array}{c}\hfill F_u\hfill \\ \hfill F_w\hfill \\ \hfill F_{\theta }\hfill \\ \hfill F_{\mathrm{\Psi }}\hfill \end{array}\right]=-\left[\begin{array}{cccc}\hfill k_{xx}\mathit{\hspace{1em}}\hfill & \hfill k_{yz}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill k_{zz}\mathit{\hspace{1em}}\hfill & \hfill k_{zx}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill u\hfill \\ \hfill w\hfill \\ \hfill \theta \hfill \\ \hfill \psi \hfill \end{array}\right]-\left[\begin{array}{cccc}\hfill C\mathit{\hspace{1em}}\hfill & \hfill C\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill C\mathit{\hspace{1em}}\hfill & \hfill C\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill \dot{u}\hfill \\ \hfill \dot{w}\hfill \\ \hfill \dot{\theta }\hfill \\ \hfill \dot{\psi }\hfill \end{array}\right]$

(29)

The first matrix is the hardness matrix and the second matrix is the viscous damping matrix. Matrices are generally asymmetric ($K_{zz}\ne K_{zz}$ and $C_{zz}\ne C_{xz}$). The components of the matrices can be a function of the rotation speed of the rotor.

By combining the matrices for a single-rotor system, the general form of the system equation will be as follows:

$$[M]\{\ddot{\delta }\}+([C]+\mathrm{\Omega }[G])\{\dot{\delta }\}+[K]\{\delta \}=\{F(t)\}$$

(30)

Where $\delta$ contains the node displacement, M is the symmetric matrix of mass, C is an asymmetric matrix containing the antisymmetric matrix of the gyroscope as a function of $\mathrm{\Omega }$, and an asymmetric matrix is obtained from the bearing profile; It is asymmetric due to the nature of the bearings. And $F(t)$ is the force corresponding to the mass nomination, asynchronous forces, or harmonic forces fixed in space.

Calculate the natural frequency, the shape of modes.

Considering the attenuation term, solving the equation does not simply lead to the discussion of eigenvalues and eigenvectors. Rather, it requires the use of more complex methods that can be examined by examining the orthogonality of fashions [20].

$$[M]\{\ddot{\delta }\}+([C]+\mathrm{\Omega }[G])\{\dot{\delta }\}+[K]\{\delta \}=\{0\}$$

(31)

To calculate the frequency and natural modes, the system vibration equation:

$\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [M]\hfill \\ \hfill [M]\mathit{\hspace{1em}}\hfill & \hfill (\mathrm{\Omega }[G]+[C])\hfill \end{array}\right]\dot{\eta }+\left[\begin{array}{cc}\hfill -[M]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [K]\hfill \end{array}\right]\eta =\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \end{array}\right]$

(32)

Defined in the above relation:

$M^{*}$	$=\left[\begin{array}{cc}\hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [M]\hfill \\ \hfill [M]\mathit{\hspace{1em}}\hfill & \hfill ([C]+\mathrm{\Omega }[G])\hfill \end{array}\right],K^{*}=\left[\begin{array}{cc}\hfill -[M]\mathit{\hspace{1em}}\hfill & \hfill [0]\hfill \\ \hfill [0]\mathit{\hspace{1em}}\hfill & \hfill [K]\hfill \end{array}\right]$	(33)
$\overrightarrow{\eta }$	$=\left[\begin{array}{c}\hfill \dot{\delta }\hfill \\ \hfill \delta \hfill \end{array}\right]$	(34)

By placing in the above relationships:

$${[M]}^{*}\{\dot{\delta }\}+{[K]}^{*}\{\delta \}=\{0\}$$

(35)

The answer to the above equation can be considered as follows:

$$\overrightarrow{\eta }=\{\mathrm{\Psi }\}{e}^{\sigma t}$$

(36)

And to solve the nonlinearity of the above equation:

$$|\sigma M^{*}+K^{*}|=0$$

(37)

Response to stimulus forces

Here the steady-state response to the constant force is calculated as follows.

$$K\delta =F_1$$

(38)

The equation is easily solved using methods such as Gaussian.

Mass nomination force, asynchronous force, constant harmonic force in space.

According to the stated equation, the equations to be solved are as follows [20]:

$[M]\{\ddot{\delta }\}+([C]+\mathrm{\Omega }[G])\{\dot{\delta }\}+[K]\{\delta \}$
$\mathrm{\hspace{1em}}=f_s\sin (w^{*}t)+f_c\cos (w^{*}t)$	(39)

Which $f_c,f_s$ are synchronous forces and mass nominee forces. By placing and specifying the expressions multiplied by $\sin (\mathrm{\Omega }t)$ and $\cos (\mathrm{\Omega }t)$, the following equation is obtained:

$\left[\begin{array}{cc}\hfill K-M{\omega }^{{*}^2}\mathit{\hspace{1em}}\hfill & \hfill -{\omega }^{*}(\mathrm{\Omega }[G]+[C])\hfill \\ \hfill {\omega }^{*}(\mathrm{\Omega }[G]+[C])\mathit{\hspace{1em}}\hfill & \hfill K-M{\omega }^{{*}^2}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\delta }_s\hfill \\ \hfill {\delta }_s\hfill \end{array}\right]=\left[\begin{array}{c}\hfill f_s\hfill \\ \hfill f_c\hfill \end{array}\right]$

(40)

3 Results and Discussion

3.1 Validation of Numerical Modeling

To ensure the correctness of the software used and check the accuracy of the results, examples were solved with Ansys and MATLAB software that the model used in this example is designed in Katia software. These examples include the ones mentioned incredible sources and compare the answers with each other. After correcting the mistakes and ensuring the correct operation of the software, it has been used in this project. Here are some suggestions on how to look or get an appointment for Ansys Workbench software.

3.2 Simple Rotor

In this example, a simple steel rotor is shown in Figure 3. The rotor has a shaft and a disk that is rigidly attached to the end of the shaft, and the center of mass of the disk is due to an imbalanced mass with the axis of rotation. It is not true. This rotor is under rotation with a natural frequency of 700 RPM. The first six frequencies, the Campbell diagram, the harmonic response, and the shape of the modes are given below.

Figure 3 Simple rotor with the disk.

Table 1 The specifications of the simple rotor

Properties	Unit
Disk diameter	0.6 m
Disk thickness	0.5 m
Shaft length	1m
Shaft diameter	0.03 m
Density	7850 kg/m
Young’s modulus	211 GPa
Poisson’s ratio	0.3
Disk mass	110 kg
Distance from the axis	0.001 m
Unbalanced mass	110 kg

3.3 Simple Rotor Analytical Modeling

Analytical modeling of the rotor is considered with four degrees of freedom, which includes two transverse states and two torsional states $\{u\}={\{yz\mathrm{\Psi }\theta \}}^T$

$[C_{\mathit{gyro}}]$	$=\left[\begin{array}{cccc}\hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\Omega }{J}_p\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -\mathrm{\Omega }{J}_p\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right],[M]=\left[\begin{array}{cccc}\hfill m\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill m\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill J_T\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill J_T\hfill \end{array}\right],$
$[k]$	$=\left[\begin{array}{cccc}\hfill k_{11}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{12}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{11}\mathit{\hspace{1em}}\hfill & \hfill -k_{12}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -k_{12}\mathit{\hspace{1em}}\hfill & \hfill k_{22}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill k_{12}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{22}\hfill \end{array}\right]$	(41)
$I$	$={\displaystyle \frac{\pi d^4}{64}},J_p={\displaystyle \frac{m{D}^2}{8}},J_T={\displaystyle \frac{m{D}^2}{16}},$
$k_{11}$	$={\displaystyle \frac{12EI}{L^3}},k_{12}={\displaystyle \frac{6EI}{L^2}},K_{22}={\displaystyle \frac{4EI}{L}}$	(42)

The effects of the gyroscopic matrix are very important during rotor analysis.

Table 2 shows the comparison of the first five frequencies, which are calculated by both analytical and numerical methods.

Table 2 Analytical and numerical response of a simple rotor under rotational loading

	Type	First	Second	Three	Four	Five
Row	of Analysis	Mode	Mode	Mode	Mode	Mode
1	Analytical response (Ansys)	2.4128	2.4131	5.851	19.554	19.559
2	Numerical response (Matlab)	2.3117	2.391	5.671	19.1	19.451

To evaluate the validity of the modeling, a simple Campbell Rotor diagram in analytical and numerical methods is presented in Figures 4 and 5.

Figure 4 Campbell diagram (Matlab analytical modeling).

Figure 5 Campbell diagram (Ansys numerical modeling).

After reviewing the results, a Campbell diagram is drawn for a simple rotor. The modes in the complete diagram are the coefficients $1*\omega$ and $2*\omega$.

3.4 Simple Rotor Harmonic Response

To calculate the harmonic response, the mass of the unbalance is equal to the mass of the disk at point e on the axis of rotation on the disk. To calculate the balancing force:

$\mathrm{\{}F(t)\}=\{\begin{array}{c}\hfill m_{u}e{\mathrm{\Omega }}^2\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\}\cos (\mathrm{\Omega }t)+\{\begin{array}{c}\hfill 0\hfill \\ \hfill m_{u}e{\mathrm{\Omega }}^2\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\}\sin (\mathrm{\Omega }t)$

(43)

Figure 6 shows the location of the force applied to the unbalanced disk in Ansys software. This force is also applied to the mentioned point. The rotor harmonic response is also compared by analytical and numerical methods. Figure 7 shows the simple rotor harmonic response (Matlab analytical model) and, Figure 8 shows the harmonic response of the rotor (Ansys numerical model). Figures 9–11 shows the modes obtained by Ansys software. Which indicates different states of motion, including torsion and transition.

Figure 6 The location of the force applied in the software.

Figure 7 Simple rotor harmonic response (Matlab analytical model).

Figure 8 The harmonic response of the rotor (Ansys numerical model).

Figure 9 The first mode of the rotor (bending) in Ansys software.

Figure 10 The second (torsional) mode of the rotor in Ansys.

Figure 11 The third (bending) mode of the rotor in Ansys.

3.5 Rotor with Roller Bearing

In this step, a rotor with the dimensions and specifications used in the first part is shown with the difference that a roller bearing is attached to the end of the disc. As in the previous example, the rotor is made of steel. This rotor is under rotation with a natural frequency of 700 RPM. In the following, the first five frequencies, Campbell diagram, harmonic response, and shape of modes are examined. The specifications of the roller bearing are as follows, which are connected to the disc on the rotation plate.

$K_y$	$=3*{10}^5\mathrm{N}/{\mathrm{M}}^2$
$K_z$	$=6*{10}^5\mathrm{N}/{\mathrm{M}}^2$	(44)
$[k]$	$=\left[\begin{array}{cccc}\hfill k_{11}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{12}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{11}\mathit{\hspace{1em}}\hfill & \hfill -k_{12}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -k_{12}\mathit{\hspace{1em}}\hfill & \hfill k_{22}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill k_{12}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{22}\hfill \end{array}\right]+\left[\begin{array}{cccc}\hfill k_y\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_z\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]$	(45)

The amount of stiffness in the vertical direction is twice the stiffness in the transverse direction. In analytical modeling, the amount of bearing stiffness is added to the total stiffness of the matrix and the final shape of the stiffness matrix is calculated. Table 3 shows the frequency calculated from both analytical and numerical methods. In analytical modeling, as in the previous example, a rotor with four degrees of freedom is considered, which includes two transverse states and two torsional states. To evaluate the validity of the modeling, the Campbell-rotor diagram with a roller bearing is shown in analytical and numerical methods in Figures 12 and 13.

Table 3 Comparison of the first five analytical and numerical frequencies of the rotor with a roller bearing

	Type	First	Second	Three	Four	Five
Row	of Analysis	Mode	Mode	Mode	Mode	Mode
1	Analytical response (Matlab)	1.96	2.1	6.12	19.02	19.1
2	Numerical response (Ansys)	2.2374	2.25	5.9	18.93	18.97

Figure 12 Campbell diagram (Matlab analytical modeling).

Figure 13 Campbell diagram (Ansys numerical modeling).

To calculate the harmonic response of a rotor with a roller bearing, as in the previous example, an unbalanced mass equal to the mass of the disk is placed on the disk at point end of the axis of rotation. The harmonic response with the rotor shown in Figures 14 and 15 has been compared to analytically and numerically. The presence of two peaks in both diagrams indicates the existence of two critical velocities during the period and the maximum value in the frequency response diagram changes. The location is the end of the shaft. Figures 16–18 shows the modes obtained for the rotor with roller bearings obtained by Ansys software.

Figure 14 The harmonic response of the rotor with a roller bearing (Matlab analytical model).

Figure 15 The harmonic response of the rotor with a roller bearing (Ansys numerical model).

Figure 16 The first (torsional) mode of the rotor in Ansys software.

Figure 17 The second mode (bending) of the rotor in Ansys software.

Figure 18 The third (bending) mode of the rotor in Ansys software.

Figure 19 Rotor with hydrodynamic bearing.

3.6 Rotor with Hydrodynamic Bearing

In this rotor, with the dimensions and specifications used in the rotor, the same as before is used, which is shown in Figure 19. The difference is that a hydrodynamic bearing is attached to the end of the disc. As in the previous example, the rotor is made of steel. This rotor is under rotation with a natural frequency of 700 RPM. In the following, the first five frequencies, Campbell diagram, harmonic response, and shape of modes are examined.

The specifications of the selective disc-connected hydrodynamic bearing that rotates on the plate are as follows:

$$\begin{array}{c}{K}_y=3*{10}^5\mathrm{N}/{\mathrm{M}}^2\hfill \\ K_y=6*{10}^5\mathrm{N}/{\mathrm{M}}^2\hfill \end{array}\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}\begin{array}{c}{C}_y=500\mathrm{Ns}/\mathrm{m}\hfill \\ C_z=500\mathrm{Ns}/\mathrm{m}\hfill \end{array}$$

(46)

In this bearing, the amount of damping in the vertical direction is twice the horizontal direction, and as a result, the number of dampers used in the equation will be as follows.

$[C]$

$=[C_{\mathit{gyro}}]+[C]=\left[\begin{array}{cccc}\hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\Omega }{J}_p\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -\mathrm{\Omega }{J}_p\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]+\left[\begin{array}{cccc}\hfill C_y\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill C_z\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]$

The stiffness matrix is used in the same way as in the second example

$[k]$

$=\left[\begin{array}{cccc}\hfill k_{11}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{12}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{11}\mathit{\hspace{1em}}\hfill & \hfill -k_{12}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -k_{12}\mathit{\hspace{1em}}\hfill & \hfill k_{22}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill k_{12}\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_{22}\hfill \end{array}\right]+\left[\begin{array}{cccc}\hfill k_y\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill k_z\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \end{array}\right]$

(48)

In analytical modeling, as in the previous example, a rotor with four degrees of freedom is considered, which includes two transverse states and two torsional states. Table 4 shows a comparison of the first five frequencies calculated by both analytical and numerical methods. For modeling validation, the Campbell-rotor diagram with a hydrodynamic bearing is shown in analytical and numerical methods in Figures 20 and 21.

Table 4 Comparison of analytical and numerical frequencies

	Type	First	Second	Three	Four	Five
Row	of Analysis	Mode	Mode	Mode	Mode	Mode
1	Analytical response (Matlab)	6.14	0	13.97	18.96	19.1
2	Numerical response (Ansys)	6.51	0	14.57	18.93	18.99

Figure 20 Campbell diagram (analytical model).

Figure 21 Campbell diagram (Ansys model).

Figure 22 The harmonic response of the rotor with hydrodynamic bearing (Matlab analytical model).

To calculate the harmonic response of a rotor with a hydrodynamic bearing, as in the second example, an unbalanced mass equal to the mass of the disk is placed on the disk at point e of the axis of rotation. Figure 22 shows the harmonic response of the rotor with hydrodynamic bearing (Matlab analytical model) and, Figure 23 shows the harmonic response of the rotor with hydrodynamic bearing (Ansys analytical model). The rotor harmonic response is also compared by analytical and numerical methods. Figures 24–26 shows the modes obtained for the rotor with hydrodynamic bearings obtained by Ansys software.

Figure 23 The harmonic response of the rotor with hydrodynamic bearing (Ansys analytical model).

Figure 24 The first (torsional) mode of the rotor in Ansys software.

Figure 25 The second mode (bending) of the rotor in Ansys software.

Figure 26 The third mode (bending) of the rotor in Ansys software.

3.7 Simulation of RM12 Turbofan Rotor System

The turbofan engine selected in this project will be used in the Saab JAS 39 Gripe fighter. The working speed of this rotor is 30000 RPM (3140 radians per second). The general view of the turbine and the parts that are important in vibrations are shown below. This turbofan engine consists of three fan stages, seven axial compressor stages, one low-pressure turbine stage, and one high-pressure turbine stage with a shaft connected between them. The rotor is mounted on four bearings whose positions are marked. The components of the rotor are turbines and shafts, respectively, which are made of steel alloys, the compressor is made of aluminum alloys, and the fan used is made of titanium. Table 5 shows the specifications of materials used in the rotor and, Table 6 shows the weight of parts used in the gas turbine rotor.

Table 5 Specifications of materials used in the rotor

Mechanical Properties	Steel	Aluminum	Titanium
Young modulus (GPa)	207	71	116
Shear modulus (GPa)	80	26	44
Poisson’s ratio	0.3	0.35	0.35
Volumetric mass (kg/m3)	7800	2700	4506
Thermal expansion coefficient (1/C)	14 $\times$ 10${}^6$	23 $\times$ 10${}^6$	8.6 $\times$10 ${}^6$

Table 6 Weight of parts used in the gas turbine rotor

		Unit		The Total
Name	Material	Weight (kg)	Number	Weight (kg)
Low-pressure steel turbine	Steel	25.11	1	25.11
High-pressure steel turbine	Steel	11.78	1	11.78
High-pressure compressor	Steel + Aluminum	18.31	1	18.31
Fan set	Steel + Titanium	29.13	1	29.13
The whole set of rotors	–	–	–	139.68

All turbofan rotor parts have been modeled in Catia software after measuring dimensions in millimeters and finally, the assembly has been disassembled. The following are the rotor parts. Figure 27 shows the first stage spool (low-pressure rotor) and, Figure 28 shows the second-stage spool (high-pressure rotor) and, Figure 29 shows the complete model of the RM12 rotor. The main model of turbofan with blades is shown in Figure 30. The turbofan rotor finite element model is shown in Figure 31.

Figure 27 First stage spool (low-pressure rotor).

Figure 28 Second stage spool (high-pressure rotor).

Figure 29 The complete model of the RM12 rotor.

Figure 30 Turbofan RM12 model.

Figure 31 Turbofan rotor finite element model.

Table 7 Results of the round floor bearing Ansys program

Ansys Results (Hz Frequency)	Critical Velocity (Hz)	Mode Number
473	951.51	1BW
1293	1290.76	2FW
1353	1342.03	3FW
1387	1372.45	4BW
1449	1404.93	5FW
1581	1470.06	6FW
1599	1584.55	7FW
1763	1611.78	8FW
1778	1769.42	9FW
1936	1783.6	10FW
2018	1944.74	11FW
3737	2027.22	12FW

Figure 32 Campbell rotor diagram with rigid bearing.

For the case of using the floor bearing, which of course is the closest model to the real state of the rotor in question, the results of the Ansys program are shown in Table 7. In the model of simulation of a rotor with a rigid bearing, the rotor in question does not have any rotation in the direction of the axis of rotation or movement in the directions at the place of application of the bearings and is completely restricted. The range of rotational speed of the rotor with rigid bearing (RPM) is 0–30000. The damping frequencies associated with the critical velocity for a rigid bearing rotor are shown in Table 7. Figure 32 shows a Campbell diagram of the RM12 rotor with a round bottom bearing in the applied frequency range. The rigid bearing in the RM12 rotor is shown in the sections shown in Figure 45. Figure 33–36 shows the RM12 rotor modes under rotational loading.

Figure 33 Location of round floor bearings.

Figure 34 The first mode (torsion) of the free rotor.

Figure 35 The second mode (bending) of the rotor.

Figure 36 The third mode (bending) of the rotor.

3.8 RM12 Rotor with Roller Bearing

To use the roller bearing in the RM12 rotor, the stiffness parameters of this bearing can be defined as follows: The location of the bearings is as shown in Figure 37.

$$\begin{array}{c}{K}_{11}=1\times {10}^7\mathrm{N}/{\mathrm{M}}^2\hfill \\ K_{22}=1\times {10}^7\mathrm{N}/{\mathrm{M}}^2\hfill \\ K_{12}=2\times {10}^6\mathrm{N}/{\mathrm{M}}^2\hfill \\ K_{21}=2\times {10}^6\mathrm{N}/{\mathrm{M}}^2\hfill \end{array}$$

(49)

Figure 37 RM12 rotor with roller bearing.

The bearings are shown in the locations shown in the figure. In this bearing, the rotational movement in the direction of the axis is restricted. The lame frequencies along with the critical speed for a rotor with rigid bearings are shown in Table 8. Figure 38 shows a Campbell RM12 rotor diagram with a roller bearing in the applied frequency range. Figures 39–41 shows the RM12 rotor modes of a roller bearing.

Table 8 Results from the Ansys rotor program with roller bearings

Ansys Results (Hz Frequency)	Critical Velocity (Hz)	Mode Number
1288	1419.42	1FW
1348	1396.88	2BW
1392	1410.83	3BW
1455	1487.9	4BW
1573	1599.56	5BW
1605	1678.09	6BW
1758	1766.91	7FW
1778	1723.47	8BW
1945	1942.19	9BW
2002	2109.1	10FW
2382	2313.69	11FW

Figure 38 Campbell rotor diagram with roller bearing.

Figure 39 The first mode (torsion) of the rotor.

Figure 40 The second mode (bending) of the rotor.

Figure 41 The third mode (bending) of the rotor.

Figure 42 COMBI element 214.

Figure 43 RM12 rotor schematic.

Figure 44 RM12 rotor cross-section.

Table 9 Results of Ansys Rotor software with hydrodynamic bearing (asymmetric)

Ansys Results (Hz Frequency)	Critical Velocity (Hz)	Mode Number
1287	1298.09	1BW
1348	1367.47	2BW
1381	1392.83	3FW
1442	1465.85	4BW
1573	1586.88	5BW
1591	1605.25	6BW
1751	1769.91	7BW
1765	1785.51	8FW
1918	1942.03	9BW
1999	2051.56	10BW
2039	2313.21	11BW

Figure 45 Campbell rotor diagram with hydrodynamic bearing (asymmetric).

Table 10 Results of Ansys Rotor program with hydrodynamic bearing (symmetrical)

Ansys Results (Hz Frequency)	Critical Velocity (Hz)	Mode Number
1299	1305.12	1BW
1362	1378.02	2BW
1401	1409.17	3FW
1465	1481.17	4BW
1591	1616.01	5BW
1605	1621.90	6BW
1774	1789.14	7BW
1786	1799.02	8FW
1929	1941.11	9BW
2024	2037.02	10BW
2046	2311.84	11BW

Figure 46 Campbell rotor diagram with hydrodynamic bearing (symmetrical).

Figure 47 Place of application of harmonic force.

Figure 48 The first (torsional) mode of the rotor.

Figure 49 The second mode (bending) of the rotor.

Figure 50 Three mode (bending) of the rotor.

Figure 51 The harmonic response of free RM12 rotor to unbalanced mass.

Figure 52 The harmonic response of RM12 rotor with rigid bearing to unbalanced mass.

Figure 53 The harmonic response of RM12 rotor with roller bearing to unbalance mass.

Figure 54 The harmonic response of RM12 rotor with hydrodynamic bearing to unbalanced mass.

3.9 RM12 Rotor with Hydrodynamic Bearing

To use the hydrodynamic bearing in the RM12 rotor, the stiffness parameters, and the damper properties of this bearing can be defined as follows (Srikrishnanivas, 2012). This bearing is inserted in the software in the form of COMIN214 in the simulation, which is a two-dimensional asymmetric element of spring and damper. The COMIN214 element and the shape of the rotor in question along with the bearing element are shown in Figure 42. Figure 43. RM12 shows the rotor schematic and, Figure 44 shows the location of the bearings. In this bearing, the rotational motion is limited in the direction of the axis. The Campbell rotor diagram in question shows in the range of 0–30.000 RPM. Crank frequencies and critical velocities for a rotor with a hydrodynamic (asymmetric) bearing are shown in Table 9. Figure 45 shows a Campbell diagram of an RM12 rotor with a hydrodynamic (asymmetric) bearing in the applied frequency range. Table 10 shows the critical frequency and velocity values for the bearing symmetrical element. Figure 46 shows a Campbell RM12 rotor diagram with a hydrodynamic (symmetrical) bearing in the applied frequency range. Figure 47 shows the location of the harmonic force. Figures 48–50 shows the modes of the RM12 rotor under rotational loading and the use of a hydrodynamic (asymmetric) bearing. Figure 51 shows the harmonic response of free RM12 rotor to unbalanced mass and, Figure 52 shows the harmonic response of RM12 rotor with rigid bearing to unbalanced mass and, Figure 53 shows the harmonic response of RM12 rotor with roller bearing to unbalance mass and, Figure 54 shows the harmonic response of RM12 rotor with hydrodynamic bearing to unbalanced mass.

$\begin{array}{c}{K}_{11}=1\times {10}^7\mathrm{N}/{\mathrm{M}}^2\hfill \\ K_{22}=1\times {10}^7\mathrm{N}/{\mathrm{M}}^2\hfill \\ K_{12}=2\times {10}^6\mathrm{N}/{\mathrm{M}}^2\hfill \\ K_{21}=2\times {10}^6\mathrm{N}/{\mathrm{M}}^2\hfill \\ C_{11}=100\mathrm{Ns}/\mathrm{m}\hfill \\ C_{22}=100\mathrm{Ns}/\mathrm{m}\hfill \\ C_{12}=200\mathrm{Ns}/\mathrm{m}\hfill \\ C_{21}=200\mathrm{Ns}/\mathrm{m}\hfill \end{array}$

(50)

To calculate the unbalance response in the RM12 rotor, an unbalanced mass $1\times {10}^{-5}$ is placed on the second-row fan blade and a rotational speed in the range (0–18000 RPM) is applied over some time. The frequency response has been calculated for all four types of rotors along with the bearings considered in the previous section. To calculate the frequency response:

$$\{F(t)\}=\left\{\begin{array}{c}\hfill me{\mathrm{\Omega }}^2\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right\}\cos (\mathrm{\Omega }t)+\left\{\begin{array}{c}\hfill 0\hfill \\ \hfill me{\mathrm{\Omega }}^2\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right\}\sin (\mathrm{\Omega }t)$$

(51)

4 Conclusions

The results of numerical and analytical simulations had a very good agreement with each other.

The main frequencies of the rotor were calculated at the specified rotational speed, the results of which are fully specified.

The main values of critical velocities are obtained for different rotor states, which are reduced when using the hydrodynamic journal bearing. Which itself is effective in system stability.

According to the results, as expected, in the free rotor mode, the highest lame frequencies are observed in the system.

Due to the shape of the modes obtained in each case, the most displacement has occurred in the RM12 rotor fan, so reducing the vibrations of the input fan to prevent failure, should be considered as an important parameter.

In the bearing modeling mode, when the bearings are rigidly modeled, the crank frequencies are higher than in the case of the hydrodynamic and rolling journal bearing mode, due to the higher system stiffness ratio. It is the same as before. Also, the use of hydrodynamic journal bearings reduces the system frequency, which in turn increases the efficiency and reduces the system vibrations.

In calculating the frequency response of the RM12 rotor using an unbalanced mass that was applied to the input fan, in the case of using a hydrodynamic bearing, the amplitude of the system vibrations is greatly reduced.

References

[1] Lalanne M, Ferraris G. Rotordynamics prediction in engineering. Wiley; 1998 Feb 4.

[2] Dimarogonas AD, Paipetis SA, Chondros TG. Analytical methods in rotor dynamics. Springer Science & Business Media; 2013 Feb 19.

[3] Lei S, Palazzolo A. Control of flexible rotor systems with active magnetic bearings. Journal of Sound and Vibration. 2008 Jul 8;314(1-2):19–38.

[4] Dimentberg FM. Flexural vibrations of rotating shafts. Butterworths’; 1961.

[5] Ruhl, Roland L., and J. F. Booker. “A finite element model for distributed parameter turborotor systems.” (1972): 126–132.

[6] Bai C, Zhang H, Xu Q. Experimental and numerical studies on nonlinear dynamic behavior of rotor system supported by ball bearings. Journal of engineering for gas turbines and power. 2010 Aug 1;132(8).

[7] Creci G, Menezes JC, Barbosa JR, Corra JA. Rotordynamic analysis of a 5-kilonewton thrust gas turbine by considering bearing dynamics. Journal of Propulsion and Power. 2011 Mar;27(2):330–336.

[8] Taplak H, Parlak M. Evaluation of gas turbine rotor dynamic analysis using the finite element method. Measurement. 2012 Jun 1;45(5):1089–1097.

[9] Wang S, Zhu R. Study on load sharing behavior of coupling gear-rotor-bearing system of GTF aero-engine based on multi-support of rotors. Mechanism and Machine Theory. 2020 May 1;147:103764.

[10] Kim S, Jung K, Lee J, Park K, Yang K. Rotordynamic Analysis of a Dual-Spool Turbofan Engine with Focus on Blade Defect Events. Tribology and Lubricants. 2020;36(2):105–115.

[11] Sheridan, W. G., & McCune, M. E. U.S. Patent No. 10,612,555. Washington, DC: U.S. Patent and Trademark Office. 2020.

[12] Dinc A, Elbadawy I. Global warming potential optimization of a turbofan powered unmanned aerial vehicle during surveillance mission. Transportation Research Part D: Transport and Environment. 2020 Aug 1;85:102472.

[13] Ehrich FF, editor. Handbook of rotordynamics. Krieger Publishing Company; 1999.

[14] Cao H, Niu L, Xi S, Chen X. Mechanical model development of rolling bearing-rotor systems: A review. Mechanical Systems and Signal Processing. 2018 Mar 1;102:37–58.

[15] Xiang J, Jiang Z, Chen X. A class of wavelet-based Rayleigh-Euler beam element for analyzing rotating shafts. Shock and Vibration. 2011 Jan 1;18(3):447–458.

[16] Nelson, H. D. “A finite rotating shaft element using Timoshenko beam theory.” (1980): 793–803.

[17] Kaghazian A, Foruzande HR, Hajnayeb A, Mohammad Sedighi H. nonlinear free vibrations analysis of a piezoelectric bimorph nanoactuator using nonlocal elasticity theory. Modares Mechanical Engineering. 2016 Jun 10;16(4):55–66.

[18] Palaninathan R, Chandrasekharan PS. Curved beam element stiffness matrix formulation. Computers & structures. 1985 Jan 1;21(4):663–669.

[19] Weaver Jr W, Timoshenko SP, Young DH. Vibration problems in engineering. John Wiley & Sons; 1990 Feb 14.

[20] Elhami, M R, Mazdgir, A, Vibration Modeling and Analysis of Complete Rotor Set by Finite Element Method, 2nd International Conference on Acoustics and Vibrations, Tehran, 2012.

[21] Jamshidi H, Gholi H, Vibration analysis of rotor-disk assembly with nonlinear bearing Journal of Technical and Engineering Modares, 10(20), 2005. 35–44.

[22] Srikrishnanivas, D. Rotor dynamic analysis of RM12 jet engine rotor using ANSYS. 2012.

Biographies

Hossein Rahmani received his B.Sc. in the Shahid Bahonar University of Kerman and M.Sc degrees in Mechanical Engineering from IHU University, Tehran, Iran. He has written numerous educational, research, and administrative records in his career. He is also capable and interested in research, Rotor dynamic, numerical simulation, and vibrational analysis.

Mohammadreza Elhami received a Ph.D. degree in Mechanical Engineering from the University of Liverpool. Now he works as an Associate Professor at the school of mechanical engineering, IHU University, tehran, Iran. His research interests include vibrational and dynamic analysis.

Amin Moslemi Petrudi received his B.Sc. in Automobile Engineering from Allameh Amini University and M.Sc degrees in Mechanical Engineering from IHU University, Tehran, Iran. He has written numerous educational, research, and administrative records in his career. He is also capable and interested in research, design, modeling and simulation, impact and penetration mechanics, and stress analysis.

European Journal of Computational Mechanics, Vol. 31_2, 275–318.
doi: 10.13052/ejcm2642-2085.3125
© 2022 River Publishers

Abstract

1 Introduction

2 Statement of the Problem

2.1 Vibrations of Rotating Axes

2.2 Kinematics of Rotating Elements

2.3 Rail Shaft Element

2.4 Element Shaft Euler-Bernoulli

2.5 Element Shaft Timoshenko

3 Results and Discussion

3.1 Validation of Numerical Modeling

3.2 Simple Rotor

3.3 Simple Rotor Analytical Modeling

3.4 Simple Rotor Harmonic Response

3.5 Rotor with Roller Bearing

3.6 Rotor with Hydrodynamic Bearing

3.7 Simulation of RM12 Turbofan Rotor System

3.8 RM12 Rotor with Roller Bearing

3.9 RM12 Rotor with Hydrodynamic Bearing

4 Conclusions

References

Biographies