Free Vibration Analysis of a Thermally Loaded Porous Functionally Graded Rotor–Bearing System

Abstract

The present work deals with natural and whirl frequency analysis of a porous functionally graded (FG) rotor–bearing system using the finite element method (FEM). Stiffness, mass and gyroscopic matrices are derived for porous and non-porous FG shafts by developing a novel two-noded porous FG shaft element using Timoshenko beam theory (TBT), considering the effects of translational inertia, rotatory inertia, gyroscopic moments and shear deformation. A functionally graded shaft whose inner core is comprised of stainless steel (SS) and an outer layer made of ceramic (ZrO₂) is considered. The effects of porosity on the volume fractions and the material properties are modelled using a porosity index. The non-linear temperature distribution (NLTD) method based on the Fourier law of heat conduction is used for the temperature distribution in the radial direction. The natural and whirl frequencies of the porous and non-porous FG rotor systems have been computed for different power law indices, volume fractions of porosity and thermal gradients to investigate the influence of porosity on fundamental frequencies. It has been found that the power law index, volume fraction of porosity and thermal gradient have a significant influence on the natural and whirl frequencies of the FG rotor–bearing system.

1. Introduction

An important aspect to take into consideration for a superior structural performance is the material strength of the system. Traditional composite materials are impotent when they are subjected to thermo-mechanical loading, due to inter-laminar stresses which cause the de-lamination of layers. Generally, metals are preferred because of their high strength and toughness [1]. However, at high temperatures, the strength of the metal drastically deteriorates. On the other hand, ceramics are heat resistant; however, their applications are restricted due to their low toughness. The development of a new class of composites, functionally graded materials (FGMs), mitigated the problems of de-bonding, de-lamination and residual stresses in fibre-reinforced composites at elevated temperature, while making use of the advantages of both metal and ceramic material properties. FGM is an inhomogeneous micromechanical composite typically made from different phases of metal and ceramic material constituents. The volume fraction of constituent materials is arranged in the desired direction based on material laws for smooth and continuous change from one layer to another.

The term FGM was coined by a group of Japanese scientists in the mid-1980s. There are numerous types of methods based on material applications to fabricate FGMs—they include plasma spraying, chemical vapour deposition, physical vapour deposition, slurry dipping, electroplating, sintering, combustion synthesis, centrifugal casting, tape casting, hot pressing, powder metallurgy, etc. Suresh and Mortensen [2] and Miyamoto et al. [3] discussed some of the manufacturing processes of functionally graded materials. The initial functionally graded materials were made for the aerospace industry, as the highest temperature of the surface of spacecrafts is estimated to reach as high as 2100 K; hence, the material used on the surface should withstand that high temperature [4]. Moreover, on account of wide range of applications in many other industries, such as automotive, electronics, biomedical, marine, defence and construction, functionally graded materials have become a topic of interest for many researchers. Most of the previous works have reported on the dynamic characteristics of beams. Since the modelling of the FG shaft in the present work is based on Timoshenko beam theory with rotating effects, the literature review is arranged around the dynamic analysis of FG beams, FG rotors and then porous FG beams and rotors.

Many investigations were executed to analyse FG structures, and excellent review papers are reported in the literature [5,6,7]. Liu et al. [8] developed a shell element to analyse the nonlinear response of FG structures. Liu and Jeffers [9] performed an isogeometric analysis of FG sandwich plates based on layer wise displacement theory. Simsek [10] investigated a static analysis of a functionally graded simply supported beam under a uniformly distributed load by the Ritz method, using Timoshenko and higher order shear deformation beam theories. The fundamental frequency of FG beams using classical, first order and different higher order shear deformation theories for different boundary conditions was analysed [11]. Aydogdu and Taskin [12] performed a free vibration analysis of simply supported FG beams using parabolic shear deformation beam theory, classical beam theory and exponential shear deformation beam theory. Pydah and Sabole [13] performed a static analysis of an FG bi-directional circular beam based on Euler–Bernoulli theory. The free vibration analysis of Euler and Timoshenko FG beams using the Rayleigh–Ritz method subjected to different boundary conditions is investigated by Pradhan and Chakraverty [14]. Li et al. [15] developed a beam model for the accurate analysis of FG beams based on a material-based higher-order shear theory.

Rotor–bearing systems play a crucial role in the mechanical industry, and several researchers analysed rotor–bearing systems based on various models. Dimentberg [16] examined viscous and hysteretic internal damping of a rotating shaft. Ruhl and Booker [17] developed a finite element model of a turbo-rotor system by including the effects of translational kinetic energy and elastic bending energy. Nelson and McVaugh [18] generalised the Ruhl’s element by including the effects of rotatory inertia, gyroscopic moments and axial load. Extending this work, Nelson [19] developed a finite shaft element using Timoshenko beam theory by including the effects of transverse shear deformation. The works reported above are related to steel or uniform shafts. Very few works reported on the dynamic analysis of an FG rotor–bearing system. Gayen and Roy [20], based on Timoshenko beam theory (TBT), presented a work on the vibration and stability analysis of a functionally graded spinning shaft using three node beam elements. Bose and Sathujoda [21] performed natural frequency analysis of a functionally graded rotor system using a three-dimensional element, modelled in ANSYS (ANSYS 18.0, ANSYS, Canonsburg, PA, USA). As FGMs have an extensive variety of applications, it is important to consider the effect of temperature on FG structures and rotors. Recently, Mahi et al. [22] investigated the temperature-dependent free vibration analysis of FG beams with general boundary conditions. The thermo-mechanical buckling and non-vibration analysis of FG beams on a nonlinear elastic foundation are examined [23]. Zhang [24] analysed the thermal post buckling and nonlinear vibration behaviour of functionally graded beams based on the concept of physical neutral surfaces, von Karman strain–displacement relationships and higher order shear deformation theory. Kiani and Eslami [25] analysed the thermal buckling of FG beams. Very few works are reported to investigate the effect of temperature on FG rotors. The effect of thermal gradients on the natural frequencies of functionally graded rotor–bearing systems is investigated by Bose and Sathujoda [26]. The natural frequency of the FG rotor–bearing system for different temperature gradients is analysed by Gayen et al. [27]. Gayen et al. [28] also analysed the influence of temperature gradients on the whirl frequencies of a cracked rotor.

Amid all of the fabrication techniques, sintering is one of the ways to fabricate an FGM, but due to huge differences in the solidification temperatures of material constituents, micro voids and porosities are formed within the inter-layers of FGM while fabricating with this technique [29]. Wattanasakulpong et al. [30] deduced that porosities are also formed when a functionally graded material is fabricated using a multi-step sequential infiltration technique. The density and strength of an FGM are reduced due to the presence of porosities in the structure of FGM. Since the formation of porosities is inevitable while fabricating, the effect of porosity has become extremely prominent for analysing an FG system. Akbas [31] explored the forced vibration analysis of FG porous deep beams under dynamic loads. Ebrahimi and Jafari [32] investigated the vibration analysis of porous FG beams subjected to various thermal loadings based on higher order shear deformation theory carried out by utilising the Navier solution method. Ebrahimi and Mokhtari [33] presented a transverse vibration analysis of rotating porous functionally graded beams using the differential transform method. Atmane et al. [34] presented a free vibration analysis of thick, porous, functionally graded beams resting on elastic foundations using efficient shear and normal deformation beam theory. Wattanasakulpong and Chaikittiratana [35] investigated the flexural vibration analysis of a porous functionally graded beam using Timoshenko beam theory. Wattanasakulpong and Ungbhakorn [36] investigated the linear and nonlinear vibration analysis of elastically end-restrained FG beams with different porosities. Researchers have been analysing FG rotor systems with crack defects. A vibration analysis of an FG rotor–bearing system whose FG shaft had a transverse crack and a surface crack was performed using finite element formulation [27,28]. Sathujoda et al. [37] investigated the natural and whirl frequencies of a corroded FG rotor–bearing system using the finite element method (FEM). Even though a few works on the natural frequency analysis of cracked and corroded FG rotor systems are available [27,28,37] in the literature, to the best of the authors’ knowledge, research on the dynamic analysis of a porous FG rotor system is rarely found. Recently, an attempt was made by Bose and Sathujoda [38] to analyse the natural frequencies of a porous FG rotor–bearing system without thermal effects using ANSYS three-dimensional finite elements. However, this is an approximate analysis without any thermal effects and the model is not suitable for time dependent steady-state and transient dynamic vibration analyses. Since FG rotors are generally subjected to elevated temperatures to make use of FG material properties, a comprehensive and accurate study of the natural and whirl frequencies of a porous FG rotor system using efficient modelling is required, which is addressed in the present work.

This literature review reveals that the research on the vibration analysis of porous FG rotor–bearing systems is extremely limited. Since the porosity affects the dynamic characteristics of an FG rotor system, it is important to analyse the dynamic properties, such as natural and whirl frequencies, to avoid rotor failures. Our main interest in the present work is to investigate the natural and whirl frequencies of a porous FG rotor–bearing system using a novel two-node porous FG shaft element. Dimensionless natural frequencies and whirl frequencies are computed for different power law indices (k) and volume fractions of porosity (α) at different temperature gradients (∆T).

2. Materials and Methods

Material modelling of a porous FG shaft is reported in this section. Detailed descriptions of the non-linear temperature distribution (NLTD), as well as the power law gradation of nonporous and porous FG shafts, are presented in Section 2.1, Section 2.2 and Section 2.3, respectively.

Material gradation is divided into three categories: continuously graded microstructure, discretely graded microstructure and multi-phase graded microstructure [39]. Generally, in FGM, material properties (thermal conductivity K, Young’s modulus E, density and Poisson’s ratio) are diversified along the gradation direction. The radial direction is considered the gradation direction in the case of circular FG shafts; as the precise details of the graded microstructure of FGM are obscured, the volume fraction of different phases is used for material gradation. The distribution of metal and ceramic phases along the radial direction of the FG shaft is represented in Figure 1. The properties of FGM must be position and temperature dependent for a precise model of the FG shaft. The position dependency of material properties can be achieved by the Voigt model [40], which is a simple rule for mixtures of composites. The material properties of a specific layer P_l are expressed as in Equation (1).

$$P_l=P_m{V}_m+P_c{V}_c$$

(1)

where $P_m$ and $P_c$ are the material properties of metal and ceramic; $V_m$ and $V_c$ are the volume fractions of metal and ceramic of the respective layers. The temperature dependency of material properties given by Touloukian [41] can be expressed as a non-linear function of temperature, as shown in Equation (2). Although there are many other non-linear variations in the rotating systems [42,43], the non-linear gradation of material properties and non-linear temperature distribution are considered in the present work.

$$P\left(T\right) = P_0\left(P_{-1}{T}^{-1}+1+ P_{1}T+ P_2{T}^2+ P_3{T}^3\right)$$

(2)

where T is temperature in Kelvin, and $P_{-1}$, $P_0$, $P_1$, $P_2$ and $P_3$ are coefficients of temperature that are different for different material constituents. Reddy and Chin [44] listed these coefficients and they are represented in

Table 1. Thermal coefficients of the constituent materials of functionally graded materials (FGMs).

Material Properties		P0	P−1	P1	P2	P3
Steel	E	201.04 × 109	0	3.08 × 10−4	−6.534 × 10−7	0
	K	15.379	0	−1.26 × 10−3	2.09 × 10−6	−3.7 × 10−10
	v	0.326	0	−2 × 10−4	3.797 × 10−7	0
ZrO2	E	244.27 × 109	0	−1.37 × 10−3	1.21 × 10−6	−3.7 × 10−10
	K	1.7	0	1.276 × 10−4	6.648 × 10−8	0
	v	0.2882	0	1.133 × 10−4	0	0
Si3N4	E	348.4 × 109	0	−3.07 × 10−4	2.160 × 10−7	−8.9 × 10−11
	K		0	9.095 × 10−4	0	0
	v	0.24	0	0	0	0
SUS304	E	201.04 × 109	0	3.07 × 10−3	−6.534 × 10−7	0
	K		0	1.276 × 10−4	0	0
	v	0.3262	0	−2 × 10−4	3.797 × 10−7	0

.

2.1. Non-Linear Temperature Distribution (NLTD)

It is presumed that the temperature is varied only along the radial direction of the shaft. The temperature variation in a hollow cylinder at any radial distance r is given by a 1-D Fourier heat conduction equation, as shown in Equation (3). $T_m$ and $T_c$ are inner and outer temperatures of a cylinder, while $R_i\mathrm{and}{R}_o$ are the innermost and outermost radii of the cylinder, respectively, $K$ is thermal conductivity at radius r.

$$\frac{d}{dr}\left[rK\frac{dT}{dr}\right]=0$$

(3)

Considering that there is no heat generation and, for boundary conditions, $T=T_m when r= R_{i }and T=T_c when r=R_{o }$ gives the non-linear temperature distribution $T\left(r\right)$, Equation (4) is the solution of Equation (3). Only the first seven terms of a polynomial series are considered as the solution [45]. $K_{cm} = K_c- K_m$ and $R_i\le r\le R_o$

$$T\left(r\right) =T_m + \left(T_c-T_m\right)\frac{\left[{{\displaystyle \sum }}_{j=0}^5\left\{\frac{{\left(-1\right)}^j}{jk+1}{\left(\frac{K_{cm}}{K_m}\right)}^j{\left(\frac{r-R_i}{R_o-R_i}\right)}^{jk+1}\right\}\right]}{\left[{{\displaystyle \sum }}_{j=0}^5\left\{\frac{{\left(-1\right)}^j}{jk+1}{\left(\frac{K_{cm}}{K_m}\right)}^j\right\}\right]}$$

(4)

2.2. Power Law Gradation of Non-Porous FG Shaft

There are various gradation laws such as the power law, sigmoidal law and exponential law for material gradation. Amongst these gradation laws, the power law is the most widely used and is the simplest for controlling the gradation of material properties. As FGM is especially made for high-temperature conditions, it is considered that material properties vary nonlinearly with temperature (T). Thus, NLTD is always used alongside the power law. For an FG shaft with a circular cross-section, the power law is expressed as in Equation (5) [28].

$$P\left(r,T\right) = P_m+ \left(P_c-P_m\right){\left[\frac{r-R_i}{R_o-R_i}\right]}^k$$

(5)

where P(r,T) represents the varying material properties along with the radius. R_i and R_o represent the inner and outer radius of the shaft, r is the radial coordinate of the shaft and k is the power law index. The value of k can range from 0 to ∞.

2.3. Material Gradation of Porous FG Shaft

The porosity of an FG system can be classified into even porosity and uneven porosity. The even distribution of the porosities is represented in Figure 2. The distribution of porosity in an FG rectangular beam is available in the literature [36]; however, in the present work, porosities are distributed in FG circular shafts. An FG shaft whose inner core is comprised of stainless steel (SS) and an outer layer made of ceramic (ZrO₂) is considered to have porosities that are evenly distributed within the cross-section of the shaft. Typically, these porosities are formed in the shaft during the process of fabrication. Due to these porosities, the material properties of layer $P_l$ are modified as shown in Equation (6).

$$P_l=P_m\left(V_m-\frac{\alpha }{2}\right) + P_c\left(V_c-\frac{\alpha }{2}\right)$$

(6)

where α ($\ll 1$) is the porosity volume fraction for the porous FG shafts and $\alpha =0$ for the FG shaft without porosity. Here, in Equation (6), $P_m$, $P_c$ and $V_m$, $V_c$ represent material properties and volume fraction of metal and ceramic. The sum of the volume fractions of metal and ceramic is shown in Equation (7).

$$V_m+V_c=1$$

(7)

Since the FG shaft consists of porosities, the material properties of the FG shaft also depend on the volume fraction of the porosity. Hence, the modified power law for even porosity is expressed as shown in Equation (8).

$$P\left(r,T,\alpha \right) = P_m+ \left(P_c-P_m\right){\left[\frac{r-R_i}{R_o-R_i}\right]}^k- \left(P_c+P_m\right)\frac{\alpha }{2}$$

(8)

The material properties such as Young’s modulus, Poisson’s ratio and the material density of imperfect FGM with an even porosity are expressed as shown in Equations (9a)–(9c).

$$E\left(r,T,\alpha \right) =E_m\left(T\right) +\left(E_c\left(T\right) -E_m\left(T\right)\right){\left[\frac{r-R_i}{R_o-R_i}\right]}^k- \left(E_c\left(T\right)+E_m\left(T\right)\right)\frac{\alpha }{2}$$

(9a)

$$\rho \left(r,T,\alpha \right) ={\rho }_m\left(T\right) +\left({\rho }_c\left(T\right) -{\rho }_m\left(T\right)\right){\left[\frac{r-R_i}{R_o-R_i}\right]}^k - \left({\rho }_c\left(T\right)+{\rho }_m\left(T\right)\right)\frac{\alpha }{2}$$

(9b)

$$v\left(r,T,\alpha \right) =v_m\left(T\right) +\left(v_c\left(T\right) -v_m\left(T\right)\right){\left[\frac{r-R_i}{R_o-R_i}\right]}^k- \left(v_c\left(T\right)+v_m\left(T\right)\right)\frac{\alpha }{2}$$

(9c)

where $R_i$ represents the inner radius, $R_o$ represents the outer radius of the shaft, r is the radial coordinate of the shaft and k is the power law index.

A porous, functionally graded rotor–bearing system is considered in the present work. The rotor shaft is made of an FGM; it is assumed that porosity is evenly distributed in the FG shaft. The disc, made of uniform steel, is precisely located at the mid-span of the FG shaft. The ends of the shaft are mounted on linear bearings. Python (IDLE Python 3.6.1, 64-bit, Python Software Foundation, Wilmington, DE, USA) code is developed to assign the material properties using the power law and NLTD. Finite element modelling of the porous FG shaft element provides the novelty of the current research work, and is described in Section 3.

3. Finite Element Modelling of Porous FG Rotor–Bearing System

The finite element formulations have been developed using Timoshenko beam theory (TBT) for natural frequency analysis of porous and nonporous FG rotor systems. A two-noded porous FG shaft element stiffness and mass matrices are derived in this work, and is, to the best of the authors’ knowledge, not available in the literature. The effects of rotatory inertia, gyroscopic moments, elastic bending energy, translational kinetic energy, and transverse shear deformation are included in shaft element modelling. A typical FG rotor–bearing system is shown in Figure 3a. A uniform steel disc is located precisely at the mid-span of the shaft and the ends of the shaft are supported on linear bearings. The component equations of motion of the FG rotor system are discussed in the following Sections (Section 3.1, Section 3.2 and Section 3.3).

3.1. Finite Porous FG Shaft Element

The FG shaft was divided into ten finite elements. A two-node FG shaft element was developed based on Timoshenko beam theory. Each node of the element was considered to have four degrees of freedom (two translational and two rotational) as shown in Figure 3b. As the analysis was performed for natural bending frequencies, which are appropriate in the case of a Jeffcott rotor supported on linear bearings in transverse directions, the axial translational degree of freedom was not considered. The equation of motion of the uniform shaft element is mentioned in the work of Nelson and McVaugh [16]. The equation of motion of the porous FG shaft element (Equation (10)) can be developed by the application of Hamilton’s extended principle and by using the equation of energy and work functions discussed in [19]. However, in the present work, the individual matrices are also the functions of the volume fraction of porosity, the radial variation in the material properties and the temperature gradients, unlike the uniform shaft element.

$$\left(\left[{\mathrm{M}}^{\mathrm{e}}\right] + \left[{\mathrm{N}}^{\mathrm{e}}\right]\right)\left\{\ddot{q}\right\} - \mathsf{\Omega }\left[{\mathrm{G}}^{\mathrm{e}}\right]\left\{\dot{q}\right\} + \left[{\mathrm{K}}^{\mathrm{e}}\right]\left\{q\right\} = \left\{{\mathrm{Q}}_c^e\right\}\cos \mathsf{\Omega }\mathrm{t}+ \left\{Q_s^e\right\}\sin \mathsf{\Omega }\mathrm{t}$$

(10)

The coefficient matrices are

$$\left[{\mathrm{K}}^{\mathrm{e}}\right] = {\left[\mathrm{K}\right]}_0 +{\Phi }_R{\left[\mathrm{K}\right]}_1$$

(11a)

$$\left[{\mathrm{M}}^{\mathrm{e}}\right] = {\left[\mathrm{M}\right]}_0+{\Phi }_R{\left[\mathrm{M}\right]}_1+{\Phi }_R{}^2{\left[\mathrm{M}\right]}_2$$

(11b)

$$\left[{\mathrm{N}}^{\mathrm{e}}\right] = {\left[\mathrm{N}\right]}_0+{\Phi }_R{\left[\mathrm{N}\right]}_1+{\Phi }_R{}^2{\left[\mathrm{N}\right]}_2$$

(11c)

$$\left[{\mathrm{G}}^{\mathrm{e}}\right] = {\left[\mathrm{G}\right]}_0+{\Phi }_R{\left[\mathrm{G}\right]}_1+{\Phi }_R{}^2{\left[\mathrm{G}\right]}_2$$

(11d)

Equation (10) represents the elemental equation of motion of a porous FG shaft. $\mathsf{\Omega }$ is the spin speed of the shaft in $\mathrm{rad}/\sec$. $\left[{\mathrm{K}}^{\mathrm{e}}\right]$ is the elemental stiffness matrix, $\left[{\mathrm{M}}^{\mathrm{e}}\right]$ is the elemental translational mass matrix, $\left[{\mathrm{N}}^{\mathrm{e}}\right]$ is the elemental rotatory matrix and $\left[{\mathrm{G}}^{\mathrm{e}}\right]$ is the elemental gyroscopic matrix of the porous FG shaft. $\left[{\mathrm{K}}^{\mathrm{e}}\right]$, $\left[{\mathrm{M}}^{\mathrm{e}}\right]$, $\left[{\mathrm{N}}^{\mathrm{e}}\right]$ are symmetric matrices, whereas $\left[{\mathrm{G}}^{\mathrm{e}}\right]$ is a skew symmetric matrix. ${\Phi }_R$ is the transverse shear effect and is shown in Equation (12a). As the shaft element is composed of FGM, it also depends on radial distance, as well as the volume fraction of porosity and temperature.

$$\Phi \left(r,T,\alpha \right) =\frac{12E\left(r, T,\alpha \right) I\left(r\right) dr}{\kappa \left(r, T ,\alpha \right) A G\left(r, T,\alpha \right)l_e{}^2}$$

(12a)

where, $E\left(r, T ,\alpha \right)$ is Young’s Modulus, $I\left(r\right)$ is the moment of inertia, $\kappa \left(r, T,\alpha \right)$ is the transverse shear form factor, $G\left(r, T,\alpha \right)$ is the Shear Modulus of the shaft element, l_e is the length of the element and $dr$ is the thickness of the shaft element whose inner radius is $r$ and outer radius is $r+dr$. A solid FG shaft whose inner radius $R_i=0$ and outer radius $R_o=R$ is shown in Equation (12b).

$${\Phi }_R=\underset{0}{\overset{R}{\int }}\Phi \left(r,T,\alpha \right) = \underset{0}{\overset{R}{\int }}\frac{12E\left(r, T,\alpha \right) I\left(r\right) dr}{\kappa \left(r, T ,\alpha \right) A G\left(r, T ,\alpha \right)l_e{}^2}$$

(12b)

The spatial constraint matrix associated with translational shape functions is shown in Equation (13).

$$\left[\mathsf{\psi }\left(\mathrm{s}\right)\right] = \left[\begin{array}{cccccccc}{\mathsf{\psi }}_1& 0& 0& {\mathsf{\psi }}_2& {\mathsf{\psi }}_3& 0& 0& {\mathsf{\psi }}_4\\ 0& {\mathsf{\psi }}_1& -{\mathsf{\psi }}_2& 0& 0& {\mathsf{\psi }}_3& -{\mathsf{\psi }}_4& 0\end{array}\right]$$

(13)

$${\mathsf{\psi }}_{\mathrm{j}}\left(\mathrm{s}\right) = \frac{\left[{\mathsf{\alpha }}_{\mathrm{j}}\left(\mathrm{s}\right) + {\mathsf{\Phi }}_{\mathrm{R}}{\mathsf{\beta }}_{\mathrm{j}}\left(\mathrm{s}\right)\right]}{1+{\mathsf{\Phi }}_{\mathrm{R}}}$$

$$\mathrm{j}=1,2,3,4;\mu =\frac{s}{l_e}$$

$${\mathsf{\alpha }}_1=1-3{\mathsf{\mu }}^2+2{\mathsf{\mu }}^3 {\mathsf{\beta }}_1=1-\mathsf{\mu }$$

$${\mathsf{\alpha }}_2={\mathrm{l}}_{\mathrm{e}}\mathsf{\mu }\left(1-2\mathsf{\mu }+{\mathsf{\mu }}^3\right) {\mathsf{\beta }}_2={\mathrm{l}}_{\mathrm{e}}\mathsf{\mu }\left(1-\mathsf{\mu }\right)/2$$

$${\mathsf{\alpha }}_3={\mathsf{\mu }}^2\left(3-2\mathsf{\mu }\right) {\mathsf{\beta }}_3=\mathsf{\mu }$$

$${\mathsf{\alpha }}_4={\mathrm{l}}_{\mathrm{e}}{\mathsf{\mu }}^2\left(-1+\mathsf{\mu }\right) {\mathsf{\beta }}_4={\mathrm{l}}_{\mathrm{e}}\mathsf{\mu }\left(-1+\mathsf{\mu }\right)/2$$

The spatial constraint matrix associated with rotational shape functions is shown in Equation (14).

$$\left[\mathsf{\phi }\left(\mathrm{s}\right)\right] = \left[\begin{array}{c}{\mathsf{\phi }}_{\mathrm{v}}\\ {\mathsf{\phi }}_{\mathrm{w}}\end{array}\right] = \left[\begin{array}{cccccccc}0& -{\mathsf{\phi }}_1& {\mathsf{\phi }}_2& 0& 0& -{\mathsf{\phi }}_3& {\mathsf{\phi }}_4& 0\\ {\mathsf{\phi }}_1& 0& 0& {\mathsf{\phi }}_2& {\mathsf{\phi }}_3& 0& 0& {\mathsf{\phi }}_4\end{array}\right]$$

(14)

$${\mathsf{\phi }}_{\mathrm{j}}\left(\mathrm{s}\right) = \frac{\left[{\mathsf{\gamma }}_{\mathrm{j}}\left(\mathrm{s}\right)+{\mathsf{\Phi }}_{\mathrm{R}}{\mathsf{\delta }}_{\mathrm{j}}\left(\mathrm{s}\right)\right]}{1+{\mathsf{\Phi }}_{\mathrm{R}}}$$

$${\mathsf{\gamma }}_1=\mathsf{\mu }\left(6\mathsf{\mu }-6\right)/{\mathrm{l}}_{\mathrm{e}} {\mathsf{\delta }}_1=0$$

$${\mathsf{\gamma }}_2=1-4\mathsf{\mu }+3{\mathsf{\mu }}^2 {\mathsf{\delta }}_2=1-\mathsf{\mu }$$

$${\mathsf{\gamma }}_3=\mathsf{\mu }\left(-6\mathsf{\mu }+6\right)/{\mathrm{l}}_{\mathrm{e}} {\mathsf{\delta }}_3=0$$

$${\mathsf{\gamma }}_4=\mathsf{\mu }\left(3\mathsf{\mu }-2\right) {\mathsf{\delta }}_4=\mathsf{\mu }$$

These spatial constraint matrices are taken from the lines of Nelson [19]. The shape functions are established using Timoshenko beam theory. The spatial constraint matrices, which are associated with translational and rotational shape functions, are used to derive the elemental stiffness matrix, translational mass matrix, rotatory mass matrix and gyroscopic matrix of a porous FG shaft element. These elemental matrices are derived using the equations below.

Elemental stiffness matrix:

$$\left[{\mathrm{K}}^{\mathrm{e}}\right]={\mathrm{C}}_1\left(\underset{0}{\overset{l_e}{\int }}E\left(r,T,\alpha \right)I\left(r\right){\left[{\psi }^{″}\right]}^T\left[{\psi }^{″}\right]ds+\underset{0}{\overset{l_e}{\int }}\kappa \left(\mathrm{r},\mathrm{T},\alpha \right)G\left(\mathrm{r},\mathrm{T},\alpha \right)A\left(r\right){\left[{\psi }^{\prime }\right]}^T\left[{\psi }^{\prime }\right]ds\right)$$

(15a)

Elemental translatory mass matrix:

$$\left[{\mathrm{M}}^{\mathrm{e}}\right] = {\mathrm{C}}_2\underset{0}{\overset{l_e}{\int }}m{\left[\psi \right]}^T\left[\psi \right]ds$$

(15b)

Elemental rotatory mass matrix:

$$\left[{\mathrm{N}}^{\mathrm{e}}\right] = {\mathrm{C}}_3\underset{0}{\overset{l_e}{\int }}{I}_D{\left[\phi \right]}^T\left[\phi \right]ds$$

(15c)

Elemental gyroscopic matrix:

$$\left[{\mathrm{G}}^{\mathrm{e}}\right] = \left[{\mathrm{H}}^{\mathrm{e}}\right] - {\left[{\mathrm{H}}^{\mathrm{e}}\right]}^{\mathrm{T}}where \left[{\mathrm{H}}^{\mathrm{e}}\right] = {\mathrm{C}}_4\underset{0}{\overset{l_e}{\int }}{I}_P{\left[{\phi }_w\right]}^T\left[{\phi }_v\right]ds$$

(15d)

While solving the elemental matrices of a porous FG shaft element, coefficients arise for every elemental matrix. C₁, C₂, C₃ and C₄ are the coefficients of the elemental stiffness matrix, elemental translatory mass matrix, elemental rotatory mass matrix and elemental gyroscopic matrix, respectively, and these coefficients depend on the radial distance, temperature and volume fraction of porosity.

The coefficient of stiffness matrix [${\mathrm{K}}^{\mathrm{e}}$] of a solid porous FG shaft element for radius R is

$$C_1= \frac{1}{l_e{}^3\left(1+ {\mathsf{\Phi }}_{\mathrm{R}}\right)}\underset{0}{\overset{R}{\int }}E\left(r,T,\alpha \right)I\left(r\right) dr$$

(16a)

The coefficient of the translatory mass matrix $\left[{\mathrm{M}}^{\mathrm{e}}\right]$ of a solid porous FG shaft element for radius R is

$$C_2= \frac{l_e}{420{\left(1+ {\mathsf{\Phi }}_{\mathrm{R}}\right)}^2}\underset{0}{\overset{R}{\int }}\rho \left(r,T,\alpha \right)A\left(r\right) dr$$

(16b)

The coefficient of the rotatory mass matrix [${\mathrm{N}}^{\mathrm{e}}$] of a solid porous FG shaft element for radius R is

$$C_3= \frac{1}{120 l_e{\left(1+{\Phi }_R\right)}^2}\underset{0}{\overset{R}{\int }}\rho \left(r,T,\alpha \right)A\left(r\right)r^2 dr$$

(16c)

The coefficient of gyroscopic matrix [${\mathrm{G}}^{\mathrm{e}}$] of a solid porous FG shaft element for radius R is

$$C_4= \frac{1}{60l_e{\left(1+ {\Phi }_R\right)}^2}\underset{0}{\overset{R}{\int }}\rho \left(r,T,\alpha \right)A\left(r\right)r^2 dr$$

(16d)

where $I\left(r\right) = \mathsf{\pi }{r}^3$ is the moment of inertia and $A\left(r\right) = 2\mathsf{\pi }r$ is the area of a thin layer of the shaft whose inner radius is $r$, while the outer radius of the layer is $r+dr$.

3.2. Uniform Steel Disc

The equation of motion of a uniform steel disc [18] is

$$\left(\left[{\mathrm{M}}^{\mathrm{d}}\right] + \left[{\mathrm{N}}^{\mathrm{d}}\right]\right)\left\{{\ddot{q}}^d\right\} - \mathsf{\Omega }\left[{\mathrm{G}}^{\mathrm{d}}\right]\left\{{\dot{q}}^d\right\} = \left\{Q_c^d\right\}\cos \mathsf{\Omega }\mathrm{t}+\left\{Q_s^d\right\}\sin \mathsf{\Omega }\mathrm{t}$$

(17)

where $\left[{\mathrm{M}}^{\mathrm{d}}\right]$ is the translational mass matrix of the disc, $\left[{\mathrm{N}}^{\mathrm{d}}\right]$ is the rotatory mass matrix of the disc, $\left[{\mathrm{G}}^{\mathrm{d}}\right]$ is the gyroscopic matrix of the disc, $\left\{q^d\right\}$ is the nodal displacement vector of the disc and $\left\{Q^d\right\}$ is the external loading vector of the disc. Here, v and w are translational displacements, while β and τ are rotational displacements along the y and z axes.

$$\begin{gathered} \left\{q^d\right\}=\left\{\begin{array}{c}v\\ w\\ \beta \\ \tau \end{array}\right\},\left[{\mathrm{M}}^{\mathrm{d}}\right]=\left[\begin{array}{cccc}{\mathrm{m}}_{\mathrm{d}}& 0& 0& 0\\ 0& {\mathrm{m}}_{\mathrm{d}}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\left[{\mathrm{N}}^{\mathrm{d}}\right]=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& I_d& 0\\ 0& 0& 0& I_d\end{array}\right], \\ \left[{\mathrm{G}}^{\mathrm{d}}\right]=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& -{\mathrm{I}}_{\mathrm{p}}\\ 0& 0& {\mathrm{I}}_{\mathrm{p}}& 0\end{array}\right] \end{gathered}$$

(18)

3.3. Linear Bearings

The bearings used in this work obey the equation of motion [18] shown in Equation (19).

$$\left[{\mathrm{C}}^{\mathrm{b}}\right]\left\{{\dot{q}}^b\right\}+\left[{\mathrm{K}}^{\mathrm{b}}\right]\left\{q^b\right\}=\left\{Q^b\right\}$$

(19)

where $\left[{\mathrm{C}}^{\mathrm{b}}\right]$ is the damping matrix of the bearings, $\left[{\mathrm{K}}^{\mathrm{b}}\right]$ is the stiffness matrix of the bearings, $\left\{q^b\right\}$ is the nodal displacement vector and $\left\{Q^b\right\}$ is the external loading vector.

$$\left\{q^b\right\}=\left\{\begin{array}{c}v\\ w\end{array}\right\}, \begin{array}{c}\left[{\mathrm{C}}^{\mathrm{b}}\right]\end{array}=\left[\begin{array}{cc}{\mathrm{c}}^{\mathrm{b}}{}_{\mathrm{vv}}& {\mathrm{c}}^{\mathrm{b}}{}_{\mathrm{vw}}\\ {\mathrm{c}}^{\mathrm{b}}{}_{\mathrm{wv}}& {\mathrm{c}}^{\mathrm{b}}{}_{\mathrm{ww}}\end{array}\right], \begin{array}{c}\left[{\mathrm{K}}^{\mathrm{b}}\right]\end{array}=\left[\begin{array}{cc}{\mathrm{k}}^{\mathrm{b}}{}_{\mathrm{vv}}& {\mathrm{k}}^{\mathrm{b}}{}_{\mathrm{vw}}\\ {\mathrm{k}}^{\mathrm{b}}{}_{\mathrm{wv}}& {\mathrm{k}}^{\mathrm{b}}{}_{\mathrm{ww}}\end{array}\right]$$

(20)

The bearings used in the system are isotropic in nature ${\mathrm{k}}^{\mathrm{b}}{}_{\mathrm{vw}}={\mathrm{k}}^{\mathrm{b}}{}_{\mathrm{wv}}=0$ and ${\mathrm{c}}^{\mathrm{b}}{}_{\mathrm{ww}}={\mathrm{c}}^{\mathrm{b}}{}_{\mathrm{vv}}=0$.

4. System Equation of Motion and Solution Procedure

The assembled equation of motion of the FG rotor–bearing system with and without porosity is

$$\left[{\mathrm{M}}_{\mathrm{f}}\right]\left\{\ddot{q}\right\}-\mathsf{\Omega }\left[{\mathrm{G}}_{\mathrm{f}}\right]\left\{\dot{q}\right\}+\left[{\mathrm{K}}_{\mathrm{f}}\right]\left\{q\right\}=\left\{Q\right\}$$

(21)

where $\left[{\mathrm{M}}_{\mathrm{f}}\right]$ is the global mass matrix, which is incorporated into the elemental translational and rotary mass matrices of the porous FG shaft and disc mass, while $\left[{\mathrm{G}}_{\mathrm{f}}\right]$ is the global gyroscopic matrix, which includes the gyroscopic moments of the FG shaft and disc. $\left[{\mathrm{K}}_{\mathrm{f}}\right]$ is the global stiffness matrix, which includes all the elemental stiffness matrices of the porous FG shaft and the stiffness matrix $\left[{\mathrm{K}}^{\mathrm{b}}\right]$ of the bearings.

Equation (21) can be expressed as

$$\left[\mathrm{A}\right]\dot{h}+\left[\mathrm{B}\right]h=0$$

(22)

where

$$\left[\mathrm{A}\right]=\left[\begin{array}{cc}0& \left[{\mathrm{M}}_{\mathrm{f}}\right]\\ \left[{\mathrm{M}}_{\mathrm{f}}\right]& -\mathsf{\Omega }\left[{\mathrm{G}}_{\mathrm{f}}\right]\end{array}\right], \left[\mathrm{B}\right]=\left[\begin{array}{cc}-\left[{\mathrm{M}}_{\mathrm{f}}\right]& 0\\ 0& \left[{\mathrm{K}}_{\mathrm{f}}\right]\end{array}\right], h=\left\{\begin{array}{c}\dot{q}\\ q\end{array}\right\}$$

Equation (23) is assumed as the solution of Equation (22)

$$h={\mathrm{h}}_{\mathrm{o}}{\mathrm{e}}^{\mathsf{\lambda }\mathrm{t}}$$

(23)

Equation (24) can be obtained by substituting Equation (23) into Equation (22),

$$\mathsf{\lambda }\left[\mathrm{A}\right]{\mathrm{h}}_{\mathrm{o}}{\mathrm{e}}^{\mathsf{\lambda }\mathrm{t}}+\left[\mathrm{B}\right]{\mathrm{h}}_{\mathrm{o}}{\mathrm{e}}^{\mathsf{\lambda }\mathrm{t}}=0$$

(24)

Eigen values can be calculated by solving the Eigen value problem; see Equation (25).

$$\left(\left[{\mathrm{B}}^{-1}\right]\left[\mathrm{A}\right]+\mathsf{\lambda }\mathrm{I}\right){\mathrm{h}}_{\mathrm{o}}=0$$

(25)

The Eigen values are of the form

$${\mathsf{\lambda }}_{\mathrm{n}}\left(\mathsf{\Omega }\right) ={\mathsf{\xi }}_{\mathrm{n}}\left(\mathsf{\Omega }\right)\pm {\mathrm{i}\mathsf{\omega }}_{\mathrm{n}}\left(\mathsf{\Omega }\right)$$

(26)

where $\xi$ is the damping constant and $\omega$ is the natural frequency. The inverse of the imaginary part at Ω (spin speed) = 0 gives the value of the natural frequency in rad/s.

5. Validations

A finite element code was developed in Python using the formulation discussed in Section 4 to obtain the natural and whirl frequencies of the FG rotor–bearing system. The Python code was validated using the data available in the literature. A step-by-step code validation of the power law gradation, nonlinear temperature distribution, natural and whirl frequencies was carried out and is presented in the following subsections.

5.1. Non-Linear Temperature Distribution and Power Law Gradation in Circular FG Shaft

The non-linear temperature distribution as a function of radial distance is expressed as in Equation (4). A solid FG shaft with a radius of 0.05 m, an inner core (SS) temperature of 300 K and an outer ceramic layer temperature (ZrO₂) of 900 K is considered. The computed NLTD of the circular FG shafts for different power law indices, using Python code, is shown in Figure 4a. For k = 0, the temperature distribution is linear and for k = 0.5, 1, 3 and 5, the temperature distribution is nonlinear. This is because the thermal conductivity is a function of temperature. The computed values of the plot are in good accordance with the literature [28].

Power law gradation and NLTD method are used to obtain the material properties as a function of the radius and temperature of an FG shaft. The variation in the computed Young modulus of the FG (SS-ZrO₂) shaft due to the radial distance is shown in Figure 4b. The computed values are in good accordance with the literature [28].

5.2. Natural Frequencies of Simply Supported Non-Rotating FG Shaft

The dimensionless natural frequencies (DNFs) of a simply supported non-rotating FG shaft whose inner core is comprised of stainless steel (SS) with an outer layer made of zirconia (ZrO₂) are considered for validation. The power law is used for material gradation, DNFs are computed for the power law indices, k = 0.5, 1, 5, and for the slenderness ratio (R/2L) of 0.02. The length of the shaft is taken as 1.25 m and the radius of the cross-section is taken as 0.05 m. The shaft is assumed to be at 300 K. The computed values are tabulated in

Table 2. Dimensionless natural frequencies of simply supported non-rotating FG shaft.

Modes	k = 0.5		Error %	k = 1		Error %	k = 5		Error %
	[26]	Present		[26]	Present		[26]	Present
1st	3.2059	3.2031	0.087	3.1859	3.1780	0.248	3.1505	3.1328	0.562
2nd	3.2059	3.2031	0.087	3.1859	3.1780	0.248	3.1505	3.1328	0.562
3rd	6.3455	6.3492	0.058	6.3059	6.2995	0.101	6.2356	6.2098	0.414
4th	6.3455	6.3492	0.058	6.3059	6.2995	0.101	6.2356	6.2098	0.414
5th	9.3684	9.3941	0.274	9.3096	9.3206	0.118	9.2054	9.1872	0.198
6th	9.3684	9.3941	0.274	9.3096	9.3206	0.118	9.2054	9.1872	0.198

and it can be seen that the percentage of error is almost negligible. The computed values are in good accordance with the natural frequencies in the literature [28].

5.3. Natural Frequencies of Simply Supported Porous Functionally Graded Square Cross-Sectional Beam

The dimensionless fundamental natural frequency $\overline{\mathsf{\omega }}({\overline{\mathsf{\omega }}}^4={\omega }^2{\mathrm{L}}^4{\mathsf{\rho }}_{\mathrm{ss}}/{\mathrm{E}}_{\mathrm{ss}}{\mathrm{h}}^2)$ of the simply supported functionally graded beam with porosity is calculated using Python code and is validated with the published results of Ebrahimi et al. [32] in this section. The validation was performed to verify the porosity formulation and the Python code. The materials used in the porous functionally graded beam are Si₃N₄ and SUS304 and the temperature-dependent coefficients of Young’s modulus, the mass density and Poisson’s ratio of these materials are tabulated in

Table 3. Material properties of Stainless Steel and Zirconia.

Material Properties	Stainless Steel (SS)	Zirconia (ZrO2)
Young’s Modulus (GPa)	207.8	168
Density (kg/m3)	8166	5700
Poisson’s ratio	0.3	0.24

. The slenderness ratio of the beam L/h is taken as 20. The natural frequencies are validated for the different material distributions of the metal–ceramic (k = 0.5, 1, 2, 5) and porosity models (α = 0, 0.1, 0.2). The validations are tabulated in

Table 4. Dimensionless natural frequencies of a simply supported, functionally graded square cross-sectional beam with porosity.

α	k = 0.5			k = 1
	Present	[30]	Error %	Present	[30]	Error %
0	4.5283	4.5158	0.28	3.9695	3.9583	0.28
0.1	4.5947	4.5821	0.27	3.9621	3.9509	0.28
0.2	4.6806	4.6678	0.27	3.9516	3.9406	0.28
α	k = 2			k = 5
	Present	[30]	Error %	Present	[30]	Error %
0	3.5659	3.5553	0.3	3.2436	3.2332	0.32
0.1	3.5188	3.5082	0.3	3.1736	3.1634	0.33
0.2	3.4599	3.4495	0.3	3.0888	3.0787	0.33

. The computed values are in good agreement with the values found in the literature. Hence, based on this validation, we ensured that the porosity formulation and FE code used to compute the DNFs were accurate.

5.4. Whirl Frequencies of a Double-Disc Steel Rotor System

Since there is no research available to validate the whirl frequencies of the FG rotor, the whirl frequencies of the steel rotor–bearing system are considered. The computed whirl frequencies, using the developed FE code, are given in

Table 5. Backward and forward whirl (BW and FW) frequencies of a steel double-disc rotor.

Modes	0 RPM			4000 RPM
	[46]	Present	Error %	[46]	Present	Error %
1BW	13.79	13.76	0.217	13.59	13.56	0.221
1FW	13.79	13.76	0.217	13.97	13.94	0.215
2BW	43.66	43.75	0.206	40.07	40.19	0.299
2FW	43.66	43.75	0.206	46.90	46.96	0.127
3BW	114.08	115.32	1.087	95.52	96.76	1.298
3FW	114.08	115.32	1.087	131.63	132.67	0.790

. The rotor system is kept at room temperature. The length and diameter of the shaft are 1.5 m and 0.5 m, respectively. The shaft is discretised into 15 elements. The thickness of the two discs is 0.07 m; the radius of the left (first) disc is 0.14 m and the radius of the right (second) disc is 0.175 m. The stiffness of the bearings is 1 MN/m. It can be seen from Table 5 that the computed values are in good accordance with the whirl frequencies in the literature [46].

6. Results and Discussions

Jeffcott FG rotor systems with FG porous and non-porous shafts, whose inner cores are comprised of stainless steel with outer layers made of ZrO₂, as shown in Figure 3a, are considered in the present work. The shaft is divided into ten finite elements based on convergence. The power law and the non-linear temperature distribution method are used for the gradation of material properties in the radial direction. Rotor–bearing data used for this analysis are given in

Table 6. Data used in the present work.

Shaft
Length (L)	0.5 m
Diameter (D)	0.02 m
Disc
$\mathrm{Location}$	Mid-span
Mass (m)	2 kg
Polar moment of inertia (Ip)	$0.0024{\mathrm{kg}\mathrm{m}}^2$
Diametral moment of inertia (Id)	$0.0012 {\mathrm{kg}\mathrm{m}}^2$
Bearing Stiffness	105 N/m

and material properties are defined in Table 3. The effect of different parameters such as the volume fraction of porosity, power law index, and the thermal gradients of the natural and whirl frequencies are investigated and discussed in detail in this section.

6.1. Effect of Porosity on Material Properties

The material properties of the FG shaft are the functions of the power law parameters, temperature and volume fraction of porosity. The temperature of the inner core is considered to be 300 K, whereas the outer layer is considered to be 900 K. The variation in the Young modulus of the porous FG shaft through the radius is shown in Figure 5a,b for volume fractions with porosities of 0.1 and 0.2. Young’s modulus of the FG shaft reduces due to the presence of porosities in the FG shaft. Similarly, other material properties such as density and Poisson’s ratio are affected due to the presence of porosities in the FG shaft. The effect on material properties is more prevalent with the increase in the volume fraction of porosity.

6.2. Mode Shapes of Rotor–Bearing System

The first, second and third mode shapes of porous and non-porous FG rotor–bearing systems are plotted in Figure 6a–c for volume fractions of porosity, α = 0.3 and a power law index of 0.5. The dotted line represents the modal displacement curve of a non-porous FG rotor–bearing system, whereas the solid line represents the modal displacement curve of a porous FG rotor–bearing system. The plotted graphs are as expected and can be compared with the mode shapes of rotors with intermediate bearings described in the literature [47]. The natural frequencies corresponding to each mode shape are given in Figure 6a–c.

6.3. The Effect of Porosity on DNFs for Different Power Law Indices (k)

The dimensionless natural frequencies ($\overline{\omega }$) of the FG rotor–bearing system for different volume fractions of porosity are calculated using ${\overline{\omega }}^4={\mathsf{\rho }}_{\mathrm{ss}}{\mathrm{AL}\mathsf{\omega }}^2/{\mathrm{E}}_{\mathrm{ss}}\mathrm{I}$, where ${\mathsf{\rho }}_{\mathrm{ss}}$ is the density of stainless steel, $\mathrm{A}$ is the area of the cross-section of the shaft, $\mathrm{L}$ is the length of the shaft, $\mathsf{\omega }$ is the natural frequency of the shaft in rad/s, ${\mathrm{E}}_{\mathrm{ss}}$ is Young’s modulus of stainless steel at room temperature and I is a diametrical moment of inertia. Dimensionless natural frequencies (DNFs) of FG rotor–bearing systems at different volume fractions of porosity are plotted for different power law indices in Figure 7. It has been observed that the DNFs decrease with an increase in the volume fraction of porosity. The DNFs decrease due to a reduction in the material properties such as Young’s modulus and density, while there is an increase in the volume fraction of porosity. Young’s modulus and density for various volume fractions of porosity are tabulated in

Table 7. Young’s modulus and density of ZrO₂ and stainless steel (SS) at different α.

α	E (ZrO2)	E (SS)	Density (ZrO2)	Density (SS)
0	1.68 × 1011	2.078 × 1011	5700	8166
0.1	1.492 × 1011	1.8901 × 1011	5006.7	7472.7
0.2	1.304 × 1011	1.7022 × 1011	4313.4	6779.4
0.3	1.116 × 1011	1.5143 × 1011	3620.1	6086.1
0.4	9.284 × 1010	1.3264 × 1011	2926.8	5392.8

. However, the DNF curves of different power law indices converge at α = 0.4 due to a proportional drop in the stiffness and mass of the porous FG shaft for all power law indices. DNFs also decrease with the increase in k. The reason for this is explained in Section 6.4.

6.4. The Effect of Power Law Index on DNFs for Different Volume Fractions of Porosity

Figure 8 shows the variation in DNFs against the power law index for different volume fractions of porosity. It has been observed that the DNFs reach their maximum at k = 0. Since the FG shaft is completely composed of ceramic at k = 0, the ratio of Young’s modulus to the density of a ceramic is high. However, the DNFs of the system decrease with the increase in the power law index. As the power law index increases from zero to five, the composition of the FG shaft alters from an absolute ceramic to an FG shaft that is comprised of a combination of metal and ceramic. The percentage of metal increases in the FG shaft and hence the ratio of Young’s modulus to the density decreases. Eventually, the DNFs of the system are reduced. DNFs with volume fractions that have porosities (α) = 0.1, 0.2 are slightly higher than the DNFs without porosity. This is due to the fact that the reductions in the material properties lead to an increase in the ratio of stiffness and mass as compared to the ratio of the FG shaft without porosity.

6.5. The Effect of Temperature Gradient on DNFs for Different Values of k and α

The variation in DNFs with the power law index of an FG rotor system with and without porosity for different thermal gradients is plotted in Figure 9a,b. Dimensionless natural frequencies decrease with the increase in the temperature gradient, as expected. This is due to the reduction in material properties, especially Young’s modulus, at higher temperatures. However, the presence of porosities in the FG shaft result in a reduction in Young’s modulus and other material properties, which would lead to the abrupt drop in the DNFs at high temperatures (T_oi = T_o − T_i = 300 or 600). Later, the DNFs slightly increase initially with the power law index (k), but the surge in DNFs eventually becomes barely detectable.

6.6. The Effect of Power Law Index, Volume Fraction of Porosity, Thermal Gradient on Whirl Frequencies

The whirl frequencies of the porous FG rotor–bearing system at a shaft speed of 4000 rpm are calculated for different power law indices, volume fractions of porosity and temperature gradients in this work. The results are tabulated in

Table 8. Whirl frequencies of a porous FG rotor–bearing system for different values of k, α and ∆T.

Modes ∆T = 0	α = 0				α = 0.1				α = 0.2
	k = 0.5	k = 1	k = 5	k = 10	k = 0.5	k = 1	k = 5	k = 10	k = 0.5	k = 1	k = 5	k = 10
1BW	36.166	35.979	35.530	35.423	36.233	36.061	35.667	35.581	36.176	36.028	35.718	35.663
1FW	36.169	35.981	35.533	35.425	36.236	36.064	35.670	35.584	36.179	36.031	35.721	35.666
2BW	115.421	112.794	106.213	104.396	121.456	118.403	110.859	108.800	128.509	124.906	116.150	113.796
2FW	123.974	120.924	113.333	111.249	131.128	127.536	118.729	116.343	139.655	135.344	124.963	122.194
3BW	298.593	294.259	286.243	285.250	301.059	295.974	285.567	285.339	304.723	298.622	287.369	285.826
3FW	299.121	294.785	286.784	285.805	301.566	296.478	287.088	285.874	305.208	299.104	287.868	286.341
Modes ∆T = 300	α = 0				α = 0.1				α = 0.2
	k = 0.5	k = 1	k = 5	k = 10	k = 0.5	k = 1	k = 5	k = 10	k = 0.5	k = 1	k = 5	k = 10
1BW	35.384	35.307	35.128	35.102	35.326	35.288	35.225	35.240	35.109	35.128	35.231	35.301
1FW	35.387	35.310	35.131	35.105	35.330	35.292	35.228	35.243	35.113	35.132	35.234	35.304
2BW	115.009	112.453	106.028	104.253	120.921	117.968	110.638	108.637	127.783	124.330	115.882	113.605
2FW	123.727	120.717	113.217	111.159	130.816	127.279	118.595	116.242	139.247	135.014	124.803	122.079
3BW	279.577	277.400	274.872	275.731	281.708	278.952	275.531	276.345	284.924	281.370	276.713	277.439
3FW	280.103	277.924	275.413	276.286	282.213	279.454	276.051	276.880	285.405	281.849	277.211	277.953
Modes ∆T = 600	α = 0				α = 0.1				α = 0.2
	k = 0.5	k = 1	k = 5	k = 10	k = 0.5	k = 1	k = 5	k = 10	k = 0.5	k = 1	k = 5	k = 10
1BW	34.868	34.834	34.723	34.673	34.766	34.776	34.789	34.775	34.502	34.573	34.761	34.798
1FW	34.872	34.838	34.726	34.676	34.770	34.780	34.792	34.778	34.507	34.578	34.765	34.801
2BW	114.722	112.201	105.836	104.056	120.569	117.663	110.413	108.406	127.340	123.951	115.613	113.331
2FW	123.555	120.564	113.097	111.035	130.612	127.100	118.457	116.100	139.000	134.799	124.643	121.914
3BW	269.020	267.286	264.812	264.470	271.652	269.309	265.917	265.481	275.452	272.277	267.610	267.022
3FW	269.545	267.809	265.352	265.026	272.155	269.810	266.436	266.015	275.932	272.755	268.107	267.536

. As the power law index increases, the whirl frequencies of the FG rotor system decrease. The whirl frequencies decrease with the increase in the temperature gradient as expected due to the reduction in stiffness at higher thermal gradients. The split between 1BW-1FW and 3BW-3FW is minimal compared to the split between 2BW and 2FW. Campbell diagrams for k = 0.5, α = 0, 0.1 and Ω = 4000 rpm, 10,000 rpm at different temperature gradients are shown in Figure 10a–d. The split between 1FW-1BW and 3FW-3BW is barely detectable in the figures.

7. Conclusions

A dynamic analysis of a porous functionally graded rotor–bearing system has been carried out using the finite element method. The stiffness and mass matrices of a two-noded porous FG shaft are derived in this work. The FG rotor shaft’s inner core is comprised of stainless steel and its outer layer is made of ceramic ZrO₂. The effects of porosity, power law index and thermal gradients on the natural and whirl frequencies of an FG rotor system have been analysed and the following conclusions can be drawn from the analysis:

• DNFs decrease with an increase in the volume fraction of porosity (α) due to the reductions in the material properties, which lead to a drop in the stiffness and mass of the porous FG shaft.
• DNFs of the FG shaft with low volume fraction porosity are slightly higher than the DNFs of the FG shaft without porosity. The drop in stiffness and mass escalated the DNFs of the porous FG shaft.
• DNFs of the FG rotor system initially decrease with the increase in k, but later the rate of decrease in DNFs is comparatively low. This is because the percentage of metal is significantly higher at higher power law indices.
• DNFs decrease with an increase in thermal gradients—the reason being that Young’s Modulus reduces at elevated temperatures, which, in turn, reduces the stiffness of the system. However, due to the presence of porosities in the FG shaft, the DNFs drop abruptly at low values of k; however, later, the growth is hardly detectable.
• Whirl frequencies decrease with an increase in the power law index and thermal gradients. The split between the second forward and backward whirl frequencies is significant compared to the split between the first and third whirls.

As the natural and whirl frequencies are affected by various factors of a porous FG rotor system, such as the power law index, volume fraction of porosity and thermal gradients, it is important to analyse the effects of these factors to accurately predict the critical speeds of an FG rotor–bearing system in order to avoid any rotor failures due to resonance. The present analysis provides an important insight into the natural frequency analysis of a porous FG rotor system, which can be used to analyse the dynamic behaviour of a porous FG rotor system at its critical speeds. Further work in this direction could include the modelling of uneven porosity distributions in an FG rotor system in order to investigate their effects on the dynamic properties of the system.