Investigation of time-dependent mechanical loading on the start-stop process of functionally graded rotating disks under thermal conditions

Abstract

Rotating disks are used in industries, such as aerospace, automotive, and power plants. These disks are under time-dependent mechanical loading when the machine starts or stops working conditions. Thermal stresses caused by temperature changes with this time-dependent mechanical load create dangerous conditions. Therefore, the estimation of the elastic limit angular velocity and acceleration as a criterion for the initiation of plastic deformations has special importance in the start-stop process. An analytical Homotopy perturbation method is used to solve the heat transfer equation and the governing Navier equations in both radial and tangential directions of functionally graded rotating disks with non-uniform thickness. The Tamura-Tomoto-Ozawa model is implemented to calculate the yield stress at a different radius of the disk. The obtained results will be verified with the higher-order finite difference method. The effect of thickness parameter, type of thermal and boundary conditions on the angular velocity and acceleration, and also the starting radius of the plastic deformations are investigated by numerical examples. It is shown that by defining the appropriate temperature gradient on the outer surface of the disk and the parameters in the time-dependent angular velocity function, the level of stresses can be controlled and optimized in all working processes.

Keywords

Rotating disk start-stop process functionally graded material homotopy perturbation method time-dependent loading shear stress

Introduction

Rotating disks are mechanical parts that are used in a wide range of mechanical devices and equipment. The angular velocity of the rotating disks is usually constant throughout normal work condition. However, during the start or stop process of the machine, the rotating disks experience angular acceleration due to the change in angular velocity over time. For example, during the starting process of centrifugal compressors, the disks and other rotating parts are under angular acceleration, producing tangential displacement and shear stress, as well as radial and tangential stresses. If the shear stress exceeds permissible limits, it can have damaging effects on the disk and related components. Therefore, it is necessary to study rotating disks with time-dependent angular velocities and accelerations as mechanical loading to achieve optimal design specifications.

Materials of disks are exposed to high temperatures under these conditions. Therefore, materials that have a high ability to withstand both mechanical and thermal stresses should be selected. In a rotating disk with uniform thickness, there is a high difference of stress along the radius. It is better to use disks with non-uniform thickness to improve the stress distribution. For these reasons, it is highly recommended to construct rotating disks with non-uniform thickness profiles, made of functionally graded materials (FGMs).

Functionally graded materials were first used in 1984 in Japan for the construction of space equipment. These materials are non-homogeneous composite materials designed to enhance the efficiency of various components of a structure and control destructive deformations and stresses.

The goal is to achieve desired properties such as high strength, low weight, good conductivity, high resistance to corrosion, and high temperature resistance. In FGMs, the properties at each point are defined by a suitable mixing law based on the properties of the constituent components, typically metal and ceramic. Over the past two decades, extensive research has been conducted on FGM structures, particularly their application in the design of rotating disks. Estimating stress in FGM rotating disks is crucial due to its numerous engineering applications.

Real engineering problems are typically modeled based on nonlinear differential equations. However, these equations are difficult to solve exactly, and attempts to do so often fail. Research in the field of rotating disks has focused on finite element simulations and numerical solutions for disks with uniform or variable thickness and density under constant angular velocity.^1–18 However, previous research has shown that numerical methods for solving nonlinear problems are not always highly accurate and can have high error rates. Additionally, finite element solutions require commercial software and can be costly. To address the need for high-precision solutions, this paper uses analytical methods to analyze the mechanical behavior of FGM rotating disks under time-dependent mechanical and thermal loading. Specifically, the Homotopy perturbation method (HPM)^19–23 is used to obtain displacement-stress distributions in both radial and tangential directions of the FGM rotating disk under time-dependent mechanical and thermal loading conditions.

Among the first research carried out on rotating disks, the analytical solution of homogeneous elastic-plastic disks by Gamer using the Tersca yield criterion is more referable.¹ Subsequent research has investigated rotating disks from different aspects, assuming homogeneous materials with constant mechanical properties in the radial direction. For example, Eraslan presented a computational model to study the plastic deformation of annular rotating disks with variable thicknesses placed on a rigid shaft as a boundary condition. The model used the von Mises yield criterion and Swift’s law to simulate the nonlinear strain-hardening behavior of the material and calculated the angular velocity of the plastic limit for this selected model at different values of geometric parameters.² Eraslan and Orcan showed that depending on the shape of the disk thickness profile, the radial stress in the central area of the disk can exceed tangential stress. They also presented the expansion of plastic areas resulting from increasing angular velocity with displacement and stress distribution.³ In another study, Eraslan and Orcan investigated the elastic-plastic deformation of rotating disks with variable thickness under three types of boundary conditions: free, under pressure, and constraint on the inner surface and free on the outer surface. They found that depending on the type of boundary condition, the plastic zone can include one, two, or three regions, which are presented by different forms of the yield criterion.⁴

In the field of functionally graded rotating disks, several researchers have proposed different approaches to study their mechanical and thermal properties. You et al. were the first to propose the use of composite fibers to make rotating disks, assuming that the elastic modulus, thermal expansion coefficient, and density change based on a power function along the radial direction of the disk under constant angular velocity and uniform temperature distribution.⁵ Kordkheili and Naghdabadi presented a semi-analytical thermo-elastic solution for solid and hollow symmetric rotating disks made of graded materials under plane stress conditions.⁶ Bayat et al. studied the elastic deformation of FGM disks using the theory of first-order shear deformation.⁷ Hojjat et al. performed several works on rotating disks,^8–11 studying their elastic behavior using variational iteration method (VIM),⁸ variable material property method (VMP),⁹ Adomian’s decomposition method (ADM),¹⁰ and homotopy perturbation method (HPM).¹¹ The thickness and density function are considered non-uniform and parametric studies for different values of thickness and density parameters are carried out after the verification process.^8–11 Nayak et al. used a VIM to evaluate the elastic-plastic behavior of thermo-mechanically loaded FGM disks.¹²

In the field of optimization, Jafari et al. introduced classical and modern optimization methods in the minimal weight configuration of an elastic rotating disk with variable thickness and density. They used Karush-Kuhn-Tucker (KKT), simulated annealing (SA), and particle swarm methods (PSO) and found that the performance of PSO and SA methods is simpler and provides greater flexibility. The constraint of optimization is defined in such a way that the von Mises equivalent stress in the disk remains lower than the yield strength of the disk material.¹³ Recently, Alashti and Jafari investigated the effect of ductile damage models on the plastic deformation behavior of rotating disks with variable thickness. They showed that by considering these models in the simulations, more realistic predictions of the full plastic limit angular velocities in rotating disks can be predicted.¹⁴

In the case of time-dependent angular velocity, most research is done using numerical methods. Zheng et al. provided stress analysis for non-uniform thickness and time-dependent angular velocity in functionally graded rotating disks using the ordinary finite difference method (FDM).¹⁵ Salehian et al. performed a thermo-elastic analysis of functionally graded rotating hollow circular disks with variable thickness and angular velocity using the Galerkin method.¹⁶ Recently, Jafari used the homotopy perturbation method (HPM) and the finite difference method (FDM) for the elastic analysis of a rotating annular disk with non-uniform thickness and density under time-dependent angular velocity. The results show that HPM is easier, more practical, with no restrictions in use, and more convenient than FDM. Boundary conditions have a significant effect on reducing the displacement-stress levels in non-uniform disks compared to uniform disks.¹⁷ In the next research, Jafari used HPM and FDM for the elastic analysis of a rotating annular disk with non-uniform thickness and density under time-dependent angular velocity, considering tangential displacement.¹⁸ This research shows that HPM successfully handled the solutions of the equilibrium equation of a rotating disk in both radial and tangential directions simultaneously. Their results show that shear stress has an important influence on the distribution of equivalent von Mises stress in the elastic region and must be considered in general rotating disk analysis.¹⁸ However, no complete analytical solution has been provided for functionally graded rotating disks with time-dependent angular velocities and accelerations.

At the end of the introduction, it is often assumed in the literature on rotating disks that the angular velocity as mechanical loading is constant over time. However, considering the importance of time-dependent mechanical load, change of thickness profile along the radial direction, the advantage of FGM in the manufacture of a rotating disk, temperature changes in the radial direction and the need for high accuracy of analytical, it is essential to conduct a realistic analysis of rotating disks. In this paper, the displacement-stress distributions in both radial and tangential directions of FGM rotating disks under time-dependent mechanical and thermal loading conditions are studied using HPM and higher-order FDM for different boundary conditions that occur in real work environments. The elastic-linear hardening behavior is considered for the FGM stress-strain curve. The values of equivalent von Mises stress are implemented as a criterion for estimating the elastic limit of angular velocity and acceleration. The effect of shear stress on the values of equivalent von Mises stress is also considered. The results are compared to the higher-order FDM. After verifying the results, the values of critical parameters in the start-stop process of a rotating disk are determined in such a way that plastic deformations do not occur. A parametric study is carried out in the form of numerical examples for different values of the thickness parameter, the type of thermal loading, the values of parameters in angular velocity function, and the disk boundary conditions. The elastic limits and the radius of the beginning of plastic deformation in each case are computed. It shows that the presented model can calculate the displacement-stress in disks under time-dependent mechanical-thermal loads for all types of boundary conditions corresponding to real working environments. In addition, this model based on the HPM solution can be used for FGM rotating disks with any form of geometrical and mechanical properties function and combined loading conditions.

In this paper, we analyze the actual working conditions of rotating disks and present the results using analytical modeling. The study focuses on two different types of boundary conditions applied to the internal and external surfaces of the rotating disk, as shown in Figure 1. In each of these boundary conditions, the angular acceleration is transmitted from the rigid shaft to the annular disk.

Figure 1.

(a) constrained-guided and (b) constrained-free.

Geometrical-mechanical-thermal properties of rotating disk

Geometrical properties

An annular rotating disk with a variable thickness profile in the radial direction is being considered. The radius of the disk on the inner and outer surfaces is indicated by $r_{i}$ and $r_{o}$ , respectively. Rotating disks with variable thickness are modeled in a cylindrical coordinate system ( $r, θ, z$ ) as shown in Figure 2. The disk is symmetrical to the mid-plane and its profile is supposed to vary as a function of radius ( $r$ ) as shown in Figure 3¹⁰:

h (r) = h_{o} {(\frac{r}{r_{o}})}^{m_{5}}

(1)

where $m_{5}$ is geometric parameters ( $- 1 \leq m_{5} \leq 0$ ) and $h_{o}$ is the thickness of the disk at $r = r_{o}$ . The geometrical parameters of a disk are listed in Table 1.

Figure 2.

An example of a disk profile with variable thickness in a cylindrical coordinate system.

Figure 3

Disk profile for different thickness parameter $m_{5}$ .

Table 1.

Dimensional parameters of annular rotating disk.

Disk thickness: $h_{o}$	Outer radius: $r_{o}$	Inner radius: $r_{i}$
0.1 m	1 m	0.2 m

Mechanical properties

In FGM rotating disks, the inner surface is typically made of a metal material, such as aluminum, while the outer surface is made of a ceramic material. This choice is made because the outer surface is exposed to high working temperatures, and ceramics have properties that make them suitable for such conditions. By using an appropriate mixing rule, the mechanical properties of the disk can be varied based on the volume fraction of metal to ceramic, transitioning from the inner to the outer surfaces. In this particular study, the inner surface of the disk is made of aluminum, and the outer surface is made of zirconia. The mechanical properties of these materials are provided in Table 2. It is worth noting that zirconia, being a ceramic material, does not have a specific yield stress value due to the inherent brittleness of ceramics.

Table 2.

Mechanical properties of metal and ceramic in the FGM rotating disk.¹⁸

Mechanical properties	Zirconia	Aluminum
Young’s modulus: $E$ (GPa)	151	70
Poisson coefficient: $v$	0.3	0.3
Density: $ρ (\frac{k g}{m^{3}})$	5700	2700
Thermal expansion coefficient: $α (\frac{10^{- 6}}{° C})$	10	23.1
Thermal conductivity coefficient: $k (\frac{W}{m^{°} C})$	2	209
Yield stress: $σ_{°}$ (MPa)	–	300
Tangent modulus: $E_{t}$ (GPa)	–	35
Constant: $q$ (GPa)	7.5
Angular velocity constant: $ω_{0}$ ( $rad / s$ )	100
Angular velocity constant: $λ$ ( $1 / s$ )	$[- 1, + 1]$

The mechanical property $P (r)$ is defined using the law of mixing with the volume fraction function in power form¹⁸:

P (r) = P_{o} {(r / r_{o})}^{m}

(2)

The volume fraction function $V (r)$ and the property grading index $m$ can be defined as equation (3), based on equation (2). In this volume fraction function, the ceramic properties increase in the radial direction of the disk. The temperature of the fluid around the disk is generally higher than its initial temperature, and the disk will experience an increase in temperature. The use of this form of the volume function is justified, considering that the temperature resistance of ceramics is higher than that of metals.

V (r) = {(r / r_{o})}^{m}, m = \frac{L n (P_{i} / P_{o})}{L n (r_{i} / r_{o})}

(3)

The mechanical properties of the inner and outer surface of the disk are denoted by $P_{i}$ and $P_{o}$ respectively. The grading index $m$ is obtained based on the FGM properties on these surfaces. Elastic modulus, density, thermal expansion, thermal conductivity, and the disk thickness profile are varied based on Table 3, as equations (2) and (3). In Table 3, $E_{e}, ρ_{e}, α_{e}, k_{e}$ are reference values for elastic modulus, density, thermal expansion coefficient, and conductivity coefficient respectively, and $h_{e}$ is the thickness of the disk on the outer surface ( $r = r_{o}$ ). Additionally, $m_{1}, m_{2}, m_{3}, m_{4}$ are known as the grading index of these mechanical properties in the structure of FGM and $m_{5}$ is the thickness parameter. An elastic-linear hardening² model is used to model the stress-strain relationship of the disk material:

{\begin{matrix} ε = \frac{σ}{E} σ \leq σ_{°} \\ ε = \frac{σ_{°}}{E} + \frac{1}{E_{t}} (σ - σ_{°}) σ 〉 σ_{°} \end{matrix}

(4)

where $σ_{°}$ and $E_{t}$ are the material yield stress and tangent modulus, respectively.

Table 3.

Mechanical properties of functionally graded material with the definition of the parameters.¹⁷

Thickness profile	Thermal conductivity coefficient	Thermal expansion coefficient	Density	Elastic modulus
$\begin{array}{l} h (r) = h_{e} {(\frac{r}{r_{o}})}^{m_{5}} \\ h_{e} = h_{o} \\ - 1 \leq m_{5} \leq 0 \end{array}$	$\begin{array}{l} k (r) = k_{e} {(\frac{r}{r_{o}})}^{m_{4}} \\ k_{e} = k_{o} \\ m_{4} = \frac{L n (\frac{k_{i}}{k_{o}})}{L n (\frac{r_{i}}{r_{o}})} \end{array}$	$\begin{array}{l} α (r) = α_{e} {(\frac{r}{r_{o}})}^{m_{3}} \\ α_{e} = α_{o} \\ m_{3} = \frac{L n (\frac{α_{i}}{α_{o}})}{L n (\frac{r_{i}}{r_{o}})} \end{array}$	$\begin{array}{l} ρ (r) = ρ_{e} {(\frac{r}{r_{o}})}^{m_{2}} \\ ρ_{e} = ρ_{o} \\ m_{2} = \frac{L n (\frac{ρ_{i}}{ρ_{o}})}{L n (\frac{r_{i}}{r_{o}})} \end{array}$	$\begin{array}{l} E (r) = E_{e} {(\frac{r}{r_{o}})}^{m 1} \\ E_{e} = E_{o} \\ m_{1} = \frac{L n (\frac{E_{i}}{E_{o}})}{L n (\frac{r_{i}}{r_{o}})} \end{array}$

Thermal properties

In this paper, the rotating disk is subjected to a thermal field along the radial direction. The temperature distribution is approximated as one-dimensional, and the governing equation for heat transfer in the rotating disk, without considering the heat source, is expressed as follows¹⁷:

\frac{d}{dr} (rh (r) k (r) \frac{dT (r)}{dr}) = 0

(5)

In this relation $T (r)$ , $k (r)$ , and $h (r)$ are the temperature distribution function, the thermal conductivity function, and the thickness profile of the disk, respectively. The thermal boundary conditions on the inner and outer surfaces can be defined in a general form with the following expressions¹⁷:

C_{11} T (r_{i}) + C_{12} {\frac{dT (r)}{dr} |}_{r = r_{i}} = ξ_{i}

(6)

C_{21} T (r_{o}) + C_{22} {\frac{dT (r)}{dr} |}_{r = r_{o}} = ξ_{o}

(7)

where $C_{ij} (i, j = 1, 2)$ and $ξ_{ij} (i, j = 1, 2)$ are constants related to the heat transfer conditions on these surfaces. This paper applies two types of special thermal boundary conditions:

I. Type A: the temperature on the inner surface $T_{i}$ and on the outer surface $T_{o}$ of the disk is known and constant. In this case, the constants in equations (6) and (7) are summarized as follows¹³:

ξ_{i} = T_{i}, ξ_{o} = T_{o}

C_{11} = C_{21} = 1, C_{12} = C_{22} = 0

(8)

II. Type B: the temperature on the inner surface $T_{i}$ is constant and the temperature gradient on the outer surface of the disk is equal to $T_{o}$ . In this case, the constants in equations (6) and (7) are summarized as follows¹⁷:

ξ_{i} = T_{i}, ξ_{o} = T_{o}

C_{11} = C_{22} = 1, C_{12} = C_{21} = 0

(9)

Equivalent stress in FGM disk

One of the applications of studying rotating disks under time-dependent mechanical loading is to correctly predict the equivalent von Mises stress. The analysis of yield conditions in graded materials differs from that of homogeneous materials. In functionally graded materials (FGMs), the yield stress in different parts of the material has various values. According to Figure 4, the stress-strain diagram for FGMs is placed between the diagrams related to metal and ceramics. In this case, the yield stress in the rotating disk can be calculated using an experimental parameter called $q$ . This parameter is calculated using the ratio of stress difference in metal and ceramic to the strain difference in them as follows^21,22:

q = \frac{σ_{c} - σ_{m}}{ε_{c} - ε_{m}}

(10)

Figure 4.

Stress-strain diagram for metals, ceramics, and functionally graded materials.¹⁷

The value of $q$ depends on various factors and is calculated experimentally. The Tamura-Tomoto-Ozawa model presents the yield stress function ( $σ_{y} (r)$ ) based on this parameter from the following equation¹⁷:

σ_{y} (r) = σ_{ym} ((1 - V_{c} (r)) + \frac{q + E_{m}}{q + E_{c}} \frac{E_{c}}{E_{m}} V_{c} (r))

(11)

In this equation, $σ_{ym}$ is the yield stress in the metal and $V_{c} (r)$ is the volume fraction function of the ceramic. When the equivalent von Mises stress in the rotating disk made of FGM reaches the yield stress in equation (11), the plastic deformations are initiated, and the related angular velocity is named the elastic limit angular velocity of the disk ( $ω_{e}$ ). For rotating disks that are in plane stress conditions, von Mises stress is calculated based on the following relationship:

σ_{eq} (r) = \sqrt{σ_{r}^{2} + σ_{θ}^{2} - σ_{r} σ_{θ} + 3 τ_{r θ}^{2}}

(12)

In this relation, $σ_{r}$ and $σ_{θ}$ are radial and tangential components of the normal stresses, $τ_{r θ}$ is the in-plane shear stress, and $σ_{eq}$ is equivalent von Mises stress.

Time-dependent mechanical load

In this paper, we assume that the angular velocity and acceleration as mechanical loading are a function of time ( $t$ ), rather than a constant state as discussed in the previous literature.^1–14 For the convenience of a modeling rotating disk in general form, angular velocity and acceleration are defined as an exponential function¹⁸:

ω (t) = ω_{0} e^{- λ t}

(13)

α (t) = \frac{d ω (t)}{dt} = - λ ω_{0} e^{- λ t}

(14)

These relations focus on the accelerating or decelerating process on the rotating disk. In these relations, $ω_{0}$ represents the initial angular velocity, and $λ$ is a constant that influences the law of transition. If $λ 〉 0$ , it indicates a decelerating process that occurs when the rotating disk is stopping. On the other hand, an accelerating process occurs when $λ 〈 0$ and the disk starts the work condition. If $λ = 0$ , the disk is rotating with a constant angular velocity and zero angular acceleration, which corresponds to normal work. The elastic limit angular velocity ( $ω_{e}$ ) and elastic limit angular acceleration ( $α_{e}$ ) occur when the von Mises stress, as equivalent stress, exceeds the yield strength of the material of disk.

In Figure 5, the angular velocity and acceleration variations based on the parameters in Table 2 are plotted for various time values ( $t$ ). As seen in these figures, for $λ 〉 0$ , the changes in angular velocity and acceleration are completely different over time. The angular velocity is positive; it has $ω_{0}$ value at $t = 0$ and tends toward zero as time increases. However, the angular acceleration is negative; it has $- ω_{0} λ$ value at $t = 0$ and tends toward zero as time increases. For $λ = 0$ , the angular velocity remains constant over time ( $ω (t) = ω_{0}$ ) and the angular acceleration is zero. For $λ 〈 0$ , the changes in angular velocity and acceleration over time are similar.

Figure 5.

(a) Angular velocity and (b) angular acceleration, for different value of $λ$ .

Theoretical background

Governing equations of rotating disk

According to Figure 2, the FGM annular rotating disk is under symmetrical mechanical and thermal loads ( $T = T (r)$ ). Since the ratio of radius to thickness is high, it is possible to consider the disk as thin with acceptable accuracy. The displacement-stress analysis in the rotating disk is done by writing the governing equilibrium equations in Navier form. The advantage of Navier’s method is that the unknown displacements obtained from solving this form of equilibrium equation automatically satisfy the governing consistency equation of the disk (equation (15))^17,18:

ε_{r} (r) = \frac{d}{dr} (r ε_{θ} (r))

(15)

In this relation, $ε_{r} (r)$ and $ε_{θ} (r)$ are radial strain and tangential strain in the rotating disk. A disk element with all in-plane forces in the radial ( $r$ ) and tangential ( $θ$ ) directions is displaced in Figure 6. The equilibrium equations based on the time-dependent angular velocity and acceleration can be given as equations (16) and (17):

\frac{d}{dr} [σ_{r} h (r) r] + \frac{d τ_{r θ}}{d θ} h (r) - σ_{θ} h (r) = - ρ (r) ω (t)^{2} r^{2} h (r)

(16)

\frac{d}{dr} [τ_{r θ} h (r) r] + τ_{r θ} h (r) = - ρ (r) α (t) r^{2} h (r)

(17)

where $f (r) = - ρ (r) ω^{2} r^{2} h (r)$ is a radial force and $f (θ) = - ρ (r) α (t) r^{2} h (r)$ is a tangential force. Due to the plane stress and axisymmetric conditions, the stress component along the thickness ( $σ_{z}$ ) is zero and displacement fields are functions of $r$ only. Therefore, equations (16) and (17) can be simplified further to:

\frac{d}{dr} [σ_{r} h (r) r] - σ_{θ} h (r) = - ρ (r) ω^{2} (t) r^{2} h (r)

(18)

\frac{d}{dr} [τ_{r θ} (r) h (r) r] + τ_{r θ} (r) h (r) = - ρ (r) α (t) r^{2} h (r)

(19)

Figure 6.

Acting forces at the disk element.^12,13

The stress-strain relationship for this analysis are²:

σ_{r} (r) = \frac{E (r)}{1 - v^{2}} (ε_{r} (r) + v ε_{θ} (r) - (1 + v) ε^{T} (r))

(20)

σ_{θ} (r) = \frac{E (r)}{1 - v^{2}} (ε_{θ} (r) + v ε_{r} (r) - (1 + v) ε^{T} (r))

(21)

τ_{r θ} (r) = \frac{E (r)}{2 (1 + v)} γ_{r θ}

(22)

In this relations $ε^{T} (r)$ is the thermal strain can be calculated as follows¹⁷:

ε^{T} (r) = α (r) T (r)

(23)

Strain-displacement relationships can be defined as:

ε_{r} (r) = \frac{d u_{r} (r)}{dr}

(24)

ε_{θ} (r) = \frac{u_{r} (r)}{r} + \frac{1}{r} \frac{d u_{θ} (r)}{d θ} \overset{Rotationalsymmetry}{\to} ε_{θ} (r) = \frac{u_{r} (r)}{r}

(25)

\begin{matrix} γ_{r θ} (r) = \frac{1}{r} \frac{du (r)}{d θ} + \frac{d u_{θ} (r)}{dr} - \frac{u_{θ} (r)}{r} \overset{Rotationalsymmetry}{\to} γ_{r θ} (r) \\ = \frac{d u_{θ} (r)}{dr} - \frac{u_{θ} (r)}{r} \end{matrix}

(26)

In the Navier method, the governing equilibrium equations in the term of radial ( $u_{r} (r)$ ) and tangential ( $u_{θ} (r)$ ) displacements in relative direction are obtained by substituting equations (20)–(22) and (23)–(26) in the relations (18) and (19):

\begin{matrix} \frac{d^{2} u_{r} (r)}{d r^{2}} + (\frac{1}{r} + \frac{1}{h (r)} \frac{dh (r)}{dr} + \frac{1}{E (r)} \frac{dE (r)}{dr}) \frac{d u_{r} (r)}{dr} \\ + (- \frac{1}{r^{2}} + \frac{v}{rh (r)} \frac{dh (r)}{dr} + \frac{v}{rE (r)} \frac{dE (r)}{dr}) u_{r} (r) \\ - α (r) (1 + v) \frac{dT (r)}{dr} - (\frac{1}{h (r)} \frac{dh (r)}{dr} + \frac{1}{E (r)} \frac{dE (r)}{dr} + \frac{1}{α (r)} \frac{d α (r)}{dr}) \\ α (r) (1 + v) T (r) = - \frac{(1 - v^{2}) ρ (r) ω^{2} (t) r}{E (r)} \end{matrix}

(27)

\begin{matrix} \frac{d^{2} u_{θ} (r)}{d r^{2}} + [\frac{1}{r} + \frac{1}{h (r)} \frac{dh (r)}{dr} + \frac{1}{E (r)} \frac{dE (r)}{dr}] \frac{d u_{θ} (r)}{dr} \\ - [\frac{1}{r^{2}} + \frac{1}{h (r) r} \frac{dh (r)}{dr} + \frac{1}{E (r) r} \frac{dE (r)}{dr}] u_{θ} (r) \\ = \frac{2 (1 + v) ρ (r) α (t) r}{E (r)} \end{matrix}

(28)

These are non-linear and in-homogeneous second-order differential equations with different degrees of derivatives of unknown displacements. The analytical solution of equations (27) and (28) for constant angular velocity and only density change along the radius of disk without thermal effect is investigated in Hojjati and Jafari.^10,11 However, this study assumes a graded material with thermal effect, and the Homotopy Perturbation Method (HPM) is implemented to solve these equilibrium equations for time-dependent angular velocity and acceleration.

Homotopy perturbation method (HPM)

To present HPM,^19–23 we consider the following differential equation and boundary condition:

A (u) - f (r) = 0, r \in Θ,

(29)

B (u, \partial u / \partial n) = 0, r \in Γ,

(30)

In this relation, $u$ is the unknown function, $A$ is a general differential operator, $B$ is a boundary operator, $f (r)$ is a known analytic functional, and $Γ$ is the boundary of domain $Θ$ .

One of the most important steps in solving the differential equations by HPM is to find the linear and nonlinear parts of the function $A$ . Hence, function $A$ in equation (29) can be divided into two parts $L$ and $N$ . In this definition, $L$ is the linear part whereas $N$ is the nonlinear part. So, equation (29) is rewritable as:

L (u) + N (u) - f (r) = 0 .

(31)

The Homotopy perturbation structure is established as the following equation:

H (v, p) = L (v) - L (u_{°}) + pL (u_{°}) + p [N (v) - f (r)] = 0

(32)

where $v$ in HPM is a character variable defined as:

v (r, p) : Θ \times [0, 1] \to R

(33)

At equation (32), $p \in [0, 1]$ is the embedding parameter and $u_{°}$ is the first approximation that satisfies the boundary condition in equation (30). For $p = 0$ and $p = 1$ , equation (32) reduces to the following equations, respectively:

H (v, 0) = L (v) - L (u_{°}) = 0, H (v, 1) = A (v) - f (r) = 0

(34)

This relation shows that $L (v) = L (u_{°})$ and $A (v) = f (r)$ . This information is used to write equation (32). To begin the process of solving, we consider $v$ as:

v = v_{0} + p v_{1} + p^{2} v_{2} + p^{3} v_{3} + . . . .

(35)

Following, equation (32) is arranged according to the various powers of $p$ -terms ( $p^{0}, p^{1}, p^{2}, \dots . .$ ). Each coefficient of $p$ -terms is a differential equation in terms of $v$ that must be solved. After this stage, the process of changes $p$ from zero to unity is equal to changing from $u_{°}$ to $u (r)$ in $v (r, p)$ function. Finally, the best approximation for the solution is:

u = {lim}_{p \to 1} v = v_{°} + v_{1} + v_{2} + . . . . . .

(36)

The above convergence is discussed in (Refs. 19–23).

Higher-order finite difference method (HFDM)

One of the effective methods for solving difficult or impossible differential equations is HFDM (26, 27) as a numerical method. In HFDM, the approximation of derivatives plays a key role in the numerical solutions. To solve differential equations with specified boundary conditions, a set of grid points must be identified within the variable interval. The precision of the solution depends on the number of grid points. Increasing the number of these points can improve the precision of the solution to the desired degree. Central difference equations are often used to approximate the derivatives of a function $f (r)$ . Derivatives for grid points in boundary conditions are written based on the forward and backward finite difference form. Using Lagrange’s interpolation formula, the higher-order first derivative in the central-forward-backward difference form for equally grid point can be written as^24,25:

\begin{matrix} {\frac{df}{dr} |}^{r = r_{i}} \overset{Central}{\to} \frac{- f^{i + 2} + 8 f^{i + 1} - 8 f^{i - 1} + f^{i - 2}}{12 δ_{r}} + O (h^{4}), \\ {\frac{df}{dr} |}^{r = r_{i}} \overset{Forward}{\to} \frac{- 11 f^{i} + 18 f^{i + 1} - 9 f^{i + 2} + 2 f^{i + 3}}{6 δ_{r}} + O (h^{3}), \\ {\frac{df}{dr} |}^{r = r_{i}} \overset{Backward}{\to} \frac{11 f^{i} - 18 f^{i - 1} + 9 f^{i - 2} - 2 f^{i - 3}}{6 δ_{r}} + O (h^{3}) \end{matrix}

(37)

The second derivative of a function $f (r)$ for a given variable $r$ is:

\begin{matrix} {\frac{d^{2} f}{d r^{2}} |}^{r = r_{i}} \overset{Central}{\to} = \frac{- f^{i + 2} + 16 f^{i + 1} - 30 f^{i} + 16 f^{i - 1} - f^{i - 2}}{{δ_{r}}^{2}} + O (h^{4}), \\ {\frac{d^{2} f}{d r^{2}} |}^{r = r_{i}} \overset{Forward}{\to} = \frac{- f^{i + 3} + 4 f^{i + 2} - 5 f^{i + 1} + 2 f^{i}}{{δ_{r}}^{2}} + O (h^{3}), \\ {\frac{d^{2} f}{d r^{2}} |}^{r = r_{i}} \overset{Backward}{\to} = \frac{- f^{i - 3} + 4 f^{i - 2} - 5 f^{i - 1} + 2 f^{i}}{{δ_{r}}^{2}} + O (h^{3}) \end{matrix}

(38)

To obtain a correct and effective solution, this higher-order formulation in one dimension is used to solve the governing differential equations of the rotating disk. The Navier equations in both radial and tangential direction are written in the finite difference form for each grid point in the interval $(r_{i} . . r_{o})$ , as shown Figure 7. Additionally, the boundary points are also rewritten in the finite difference form. There are $N - 2$ equations for the inner points and two equations for the inner and outer surface boundary conditions of the disk. At the end of the finite difference process, a set of linear algebraic equations can be obtained in the following matrix expression, which must be solved using an appropriate method.²⁵

A \cdot U = B

(39)

Figure 7.

Gridding point of finite differences for disk in radial direction.

Analytical solution: Application of HPM

In this section, the governing Navier equations of FGM rotating disks with variable thickness in both the radial and tangential directions under time-dependent mechanical and thermal loading conditions are simultaneously solved by HPM as an analytical method to achieve the goals of this paper.

Heat transfer equation

The final form of the governing differential equation of heat transfer in an FGM disk is obtained by substituting the thermal conductivity coefficient and the disk thickness profiles from Table 3 in equation (5) as follows:

\frac{d^{2}}{d r^{2}} T (r) + \frac{(1 + m_{4} + m_{5})}{r} \frac{d}{dr} T (r) = 0

(40)

The linear and non-linear parts of equation (40) are defined in a way that ensures convergence of the solution:

L (T (r)) = \frac{d^{2}}{d r^{2}} T (r) + \frac{(1 + m_{4} + m_{5})}{r} \frac{d}{dr} T (r)

(41)

N (T (r)) = 0

(42)

f (r) = 0

(43)

The homotopy function is then constructed as:

\begin{matrix} H (T, p) = L (T (r)) - L (T_{⊙} (r)) + pL (T_{⊙} (r)) \\ + p [N (T (r)) - f (r)] = 0 \end{matrix}

(44)

Where, an unknown function $T (r)$ and an unknown function of the initial conditions $T_{⊙} (r)$ are considered. These unknown functions are defined in equations (45) and (46).

L (T (r)) = \frac{d^{2}}{d r^{2}} T (r) + \frac{(1 + m_{4} + m_{5})}{r} \frac{d}{dr} T (r)

(45)

L (T_{⊙} (r)) = \frac{d^{2}}{d r^{2}} T_{⊙} (r) + \frac{(1 + m_{4} + m_{5})}{r} \frac{d}{dr} T_{⊙} (r)

(46)

Equation (47) defines the next step in the process $T (r)$ :

T (r) = T_{0} (r) + p T_{1} (r) + p^{2} T_{2} (r) + p^{3} T_{3} (r) + . . .

(47)

By substituting equations (45)–(47) in the homotopy equation (44) and rearranging the resulting relationship according to $p$ -terms, equations (48)–(50) are obtained:

\begin{matrix} p^{0} \to - \frac{m_{4}}{r} \frac{d}{dr} T_{⊙} (r) + \frac{m_{4}}{r} \frac{d}{dr} T_{0} (r) + \frac{m_{5}}{r} \frac{d}{dr} T_{0} (r) \\ - \frac{m_{5}}{r} \frac{d}{dr} T_{⊙} (r) - \frac{d^{2}}{d r^{2}} T_{⊙} (r) + \frac{d^{2}}{d r^{2}} T_{0} (r) \\ - \frac{1}{r} \frac{d}{dr} T_{⊙} (r) + \frac{1}{r} \frac{d}{dr} T_{0} (r) \end{matrix}

(48)

\begin{matrix} p^{1} \to \frac{m_{4}}{r} \frac{d}{dr} T_{1} (r) + \frac{m_{5}}{r} \frac{d}{dr} T_{1} (r) + \frac{m_{4}}{r} \frac{d}{dr} T_{⊙} (r) \\ + \frac{m_{5}}{r} \frac{d}{dr} T_{⊙} (r) + \frac{d^{2}}{d r^{2}} T_{⊙} (r) + \frac{d^{2}}{d r^{2}} T_{1} (r) \\ + \frac{1}{r} \frac{d}{dr} T_{⊙} (r) + \frac{1}{r} \frac{d}{dr} T_{1} (r) \end{matrix}

(49)

p^{2} \to \frac{m_{4}}{r} \frac{d}{dr} T_{2} (r) + \frac{m_{5}}{r} \frac{d}{dr} T_{2} (r) + \frac{d^{2}}{d r^{2}} T_{2} (r) + \frac{1}{r} \frac{d}{dr} T_{2} (r)

(50)

To calculate the unknown functions, these differential equations must be solved. Equation (48) is solved by assuming a certain value for $T_{⊙} (r) = T_{0} (r)$ , as shown in equation (51):

T_{0} (r) = c_{1} + c_{2} r^{- m_{4} - m_{5}}

(51)

Substituting equation (51) into equation (48) and solving for yields $T_{1} (r)$ as equation (52):

T_{1} (r) = c_{3} + c_{4} r^{- m_{4} - m_{5}}

(52)

Similarly, equation (53) is obtained for equation (50):

T_{2} (r) = c_{5} + c_{6} r^{- m_{4} - m_{5}}

(53)

Since the disk is not under an external heat source, $f (r)$ function in the heat transfer differential equation is equal to zero, and resulting the same response for different coefficients of $p$ . To find the unknown function $T (r)$ , equations (51)–(53) are replaced in equation (47) and simultaneously $p \to 1$ . The optimum number of HPM iterations considering the converged result is two. Therefore, the solution of this equation is given by equation (54):

T (r) = C_{1} + C_{2} r^{- m_{4} - m_{5}}

(54)

In this relation, $C_{1}$ and $C_{2}$ are unknown parameters that are determined based on the thermal boundary conditions.

Navier equation: Radial direction

In this section, the equilibrium equation in a radial direction is solved using the HPM. The material properties function as Table 3, angular velocity function as equation (13) and temperature distribution function along the radius of the disk as equation (54) are substituted into equation (27), resulting in the Navier equation in terms of the unknown radial displacement ( $u_{r} (r)$ ):

\begin{matrix} \frac{d^{2} u_{r} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{r} (r)}{dr} + (\frac{(v (m_{1} + m_{5}) - 1)}{r^{2}}) u_{r} (r) \\ = (1 + v) α_{o} {(\frac{r}{r_{o}})}^{m_{3}} (C_{1} + C_{2} r^{- m_{4} - m_{5}}) (\frac{m_{1} + m_{3} + m_{5}}{r}) \\ + (1 + v) α_{o} {(\frac{r}{r_{o}})}^{m_{3}} \frac{d (C_{1} + C_{2} r^{- m_{4} - m_{5}})}{dr} - \frac{(1 - v^{2}) ρ_{o} {(\frac{r}{r_{o}})}^{m_{2}} ω_{o}^{2} {(e^{- λ t})}^{2} r}{E_{o} {(\frac{r}{r_{o}})}^{m_{1}}} \end{matrix}

(55)

The linear and non-linear parts of the equation (55) are defined separately:

\begin{matrix} L (u_{r} (r)) = \frac{d^{2} u_{r} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{r} (r)}{dr} \\ + (\frac{(v (m_{1} + m_{5}) - 1)}{r^{2}}) u_{r} (r) \end{matrix}

(56)

N (u_{r} (r)) = 0

(57)

\begin{matrix} f (r) = - (1 + v) α_{o} {(\frac{r}{r_{o}})}^{m_{3}} (C_{1} + C_{2} r^{- m_{4} - m_{5}}) (\frac{m_{1} + m_{3} + m_{5}}{r}) \\ - (1 + v) α_{o} {(\frac{r}{r_{o}})}^{m_{3}} \frac{d (C_{1} + C_{2} r^{- m_{4} - m_{5}})}{dr} \\ + \frac{(1 - v^{2}) ρ_{o} {(\frac{r}{r_{o}})}^{m_{2}} ω_{o}^{2} {(e^{- λ t})}^{2} r}{E_{o} {(\frac{r}{r_{o}})}^{m_{1}}} \end{matrix}

(58)

The homotopy function is constructed as follows:

\begin{matrix} H (u_{r}, p) = L (u_{r} (r)) - L (u_{⊙} (r)) + pL (u_{⊙} (r)) \\ + p [N (u_{r} (r)) - f (r)] = 0 \end{matrix}

(59)

In equation (59), an unknown function $u_{r} (r)$ and a function of the initial conditions $u_{⊙} (r)$ are defined:

\begin{matrix} L (u_{r} (r)) = \frac{d^{2} u_{r} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{r} (r)}{dr} \\ + (\frac{(v (m_{1} + m_{5}) - 1)}{r^{2}}) u_{r} (r) \end{matrix}

(60)

\begin{matrix} L (u_{⊙} (r)) = \frac{d^{2} u_{⊙} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{⊙} (r)}{dr} \\ + (\frac{(v (m_{1} + m_{5}) - 1)}{r^{2}}) u_{⊙} (r) \end{matrix}

(61)

We considered $u_{r} (r)$ as:

u_{r} (r) = {u_{r}}_{0} (r) + p {u_{r}}_{1} (r) + p^{2} {u_{r}}_{2} (r) + p^{3} {u_{r}}_{3} (r) + . . .

(62)

Next, equations (60)–(62) are replaced into equation (59) and the results are rearranged based on the powers of $p$ -terms. To calculate the unknown function $u_{r} (r)$ , the differential equations that are the coefficients of $p$ -terms must be solved. Assuming $u_{r 0} (r) = u_{⊙} (r)$ , the answer of these differential equations is substituted into equation (62) to calculate the unknown function for $u_{r} (r)$ as $p \to^{yields} 1$ .The optimum number of HPM iterations, considering the converged result, is two. Therefore, the analytical solution for $u_{r} (r)$ is given by equation (63):

\begin{matrix} u_{r} (r) = C_{3} r^{- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} + \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4}} + C_{4} r^{- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} - \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4}} - \frac{(1 + v)}{E_{0} \sqrt{m_{1}^{2} + (2 m_{5} - 4 v) m_{1} - 4 m_{5} v + 4}} \\ ((- \frac{2 r^{(- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} + \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4})} ρ_{0} ω_{0}^{2} e^{- 2 λ t} (v - 1) r e^{(- m_{1} + m_{2}) \ln (\frac{r}{r_{o}})} e^{(\frac{1}{2} m_{1} + \frac{1}{2} m_{5} - \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} + 2) \ln (r)}}{- \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} - m_{1} + 2 m_{2} + m_{5} + 6} \\ + \frac{2 r^{(- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} - \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4})} ρ_{0} ω_{0}^{2} e^{- 2 λ t} (v - 1) r e^{(- m_{1} + m_{2}) \ln (\frac{r}{r_{o}})} e^{(\frac{1}{2} m_{1} + \frac{1}{2} m_{5} + \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} + 2) \ln (r)}}{\sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} - m_{1} + 2 m_{2} + m_{5} + 6} \\ + α_{0} E_{0} (C_{2} (m_{1} + m_{3} - m_{4}) (\frac{- 2 r^{(- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} + \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4})} r e^{m_{3} \ln (\frac{r}{r_{o}})} e^{(\frac{1}{2} m_{1} - \frac{1}{2} m_{5} - \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} - 4) \ln (r)}}{- \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} + m_{1} + 2 m_{3} - 2 m_{4} - m_{5} + 2}) \\ + \frac{2 r^{(- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} - \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4})} r e^{m_{3} \ln (\frac{r}{r_{o}})} e^{(\frac{1}{2} m_{1} - \frac{1}{2} m_{5} + \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} - 4) \ln (r)}}{\sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} + m_{1} + 2 m_{3} - 2 m_{4} - m_{5} + 2}) \\ + C_{1} (m_{5} + m_{1} + m_{3}) (\frac{2 r^{(- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} - \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4})} r e^{m_{3} \ln (\frac{r}{r_{o}})} e^{(\frac{1}{2} m_{1} + \frac{1}{2} m_{5} + \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} - 4) \ln (r)}}{\sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} + m_{1} + 2 m_{3} + m_{5} + 2} \\ - \frac{2 r^{(- \frac{1}{2} m_{1} - \frac{1}{2} m_{5} + \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4})} r e^{m_{3} \ln (\frac{r}{r_{o}})} e^{(\frac{1}{2} m_{1} + \frac{1}{2} m_{5} - \frac{1}{2} \sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} - 4) \ln (r)}}{\sqrt{m_{1}^{2} + 2 m_{1} m_{5} - 4 m_{1} v + m_{5}^{2} - 4 m_{5} v + 4} + m_{1} + 2 m_{3} + m_{5} + 2}))) \end{matrix}

(63)

In this relation, $C_{3}$ and $C_{4}$ are constant that can be calculated by applying boundary conditions on both inner and outer surfaces of the disk in a radial direction. The radial displacement is a function of the thickness parameter ( $m_{5}$ ), the Poisson’s ratio ( $v$ ), the grading index of material properties ( $m_{1}, m_{2}, m_{3}, m_{4}$ ), the initial angular velocity ( $ω_{0}$ ), the constant ( $λ$ ), the variable time ( $t$ ) and $C_{1}$ and $C_{2}$ as unknown parameters in the thermal equation.

Navier equation: Tangential direction

In this section, the governing equilibrium equation in the tangential direction is solved by HPM. This equation is obtained with by substituting the material properties function as Table 3 and the angular acceleration function as equation (14) in to equation (28) as shown in equation (64):

\begin{matrix} \frac{d^{2} u_{θ} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{θ} (r)}{dr} - (\frac{(1 + m_{1} + m_{5})}{r^{2}}) u_{θ} (r) \\ = \frac{2 (1 + v) ρ_{0} {(\frac{r}{r_{o}})}^{m_{2}} (λ ω_{0} e^{- λ t}) r}{E_{0} {(\frac{r}{r_{o}})}^{m_{1}}} \end{matrix}

(64)

The linear and non-linear parts of equation (64) are defined as equations (65)–(67):

\begin{matrix} L (u_{θ} (r)) = \frac{d^{2} u_{θ} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{θ} (r)}{dr} \\ - (\frac{(1 + m_{1} + m_{5})}{r^{2}}) u_{θ} (r) \end{matrix}

(65)

N (u_{θ} (r)) = 0

(66)

f (r) = \frac{2 (1 + v) ρ_{0} {(\frac{r}{r_{o}})}^{m_{2}} (λ ω_{0} e^{- λ t}) r}{E_{0} {(\frac{r}{r_{o}})}^{m_{1}}}

(67)

The homotopy function is constructed as follows:

\begin{matrix} H (u_{θ}, p) = L (u_{θ} (r)) - L (u_{⊙} (r)) + pL (u_{⊙} (r)) \\ + p [N (u_{θ} (r)) - f (r)] = 0 \end{matrix}

(68)

Where, $u_{θ} (r)$ as an unknown function and $u_{⊙} (r)$ as a function of the initial conditions are defined as equations (69) and (70), respectively:

\begin{matrix} L (u_{θ} (r)) = \frac{d^{2} u_{θ} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{θ} (r)}{dr} \\ - (\frac{(1 + m_{1} + m_{5})}{r^{2}}) u_{θ} (r) \end{matrix}

(69)

\begin{matrix} L (u_{⊙} (r)) = \frac{d^{2} u_{⊙} (r)}{d r^{2}} + (\frac{(1 + m_{1} + m_{5})}{r}) \frac{d u_{⊙} (r)}{dr} \\ - (\frac{(1 + m_{1} + m_{5})}{r^{2}}) u_{⊙} (r) \end{matrix}

(70)

We considered $u_{θ} (r)$ as:

u_{θ} (r) = {u_{θ}}_{0} (r) + p {u_{θ}}_{1} (r) + p^{2} {u_{θ}}_{2} (r)

(71)

Next, equations (69)–(71) be replaced into equation (68) and the results rearrangement based on the powers of $p$ -terms. To calculate the unknown function $u_{θ} (r)$ , the differential equations that are coefficient of $p$ -terms must be solved. To solve these differential equations, it is assumed that $u_{θ 0} (r) = u_{⊙} (r)$ . Substituting the answer of these differential equations in equation (71), the unknown function for $u_{θ} (r)$ as $p \to^{yields} 1$ are calculated. The optimum number of HPM iterations considering the converged result is two. Therefore, the analytical solution of $u_{θ} (r)$ is given by equation (72):

\begin{matrix} u_{θ} (r) = C_{5} r + C_{6} r^{(- 1 - m_{1} - m_{5})} \\ - \frac{2 (λ ω_{0} e^{- λ t}) ρ_{0} r^{3 + m_{1} + m_{2}} {(r_{o})}^{m_{1} - m_{2}}}{(- 2 - m_{2} + m_{1}) E_{0} (4 + m_{5} + m_{2}) {(r^{m_{1}})}^{2}} \end{matrix}

(72)

Tangential displacement in equation (72) is a function of the thickness parameter ( $m_{5}$ ), the constant parameters in the elastic modulus and density functions ( $m_{1}, m_{2}$ ), the initial angular velocity ( $ω_{0}$ ), the constant ( $λ$ ), and the variable time ( $t$ ). In this relation, $C_{5}$ and $C_{6}$ are constant to be determined by applying boundary conditions of a rotating disk at inner and outer surfaces in the tangential direction.,

Numerical solution: Application of HFDM

Heat transfer equation

In this section, the governing equation for heat transfer of the FGM disk is solved by the HFDM. The process involves several steps:

I. Rewriting Derivatives: the derivatives in equation (5) are rewritten using difference relations (equations (37) and (38)):

\begin{matrix} \frac{- T^{i + 2} + 16 T^{i + 1} - 30 T^{i} + 16 T^{i - 1} - T^{i - 2}}{{δ_{r}}^{2}} \\ + \frac{(1 + m_{4} + m_{5})}{r} \frac{- T^{i + 2} + 8 T^{i + 1} - 8 T^{i - 1} + T^{i - 2}}{12 δ_{r}} = 0 \end{matrix}

(73)

II. Grid Points: the radial direction of the disk ( $r ϵ [r_{i}, r_{o}]$ ) is divided into $N$ grid points ( $i = [1, N]$ ). For the internal points ( $i = [2, N - 1]$ ), the finite difference of equation (73) is written.

III. Boundary Conditions: the boundary conditions of the heat transfer equation are also written in the finite difference form. There are two thermal boundary conditions, A and B, represented by equations (74) and (75) respectively.

Thermal boundary conditions A:

{\begin{matrix} i = 1, T (r) = T_{i} \to T^{1} = T_{i} \\ i = N, T (r) = T_{o} \to T^{N} = T_{o} \end{matrix}

(74)

Thermal boundary conditions B:

{\begin{matrix} i = 1, T (r) = T_{i} \to T^{1} = T_{i} \\ i = N, \frac{d}{dr} T (r) = T_{o} \to \frac{11 T^{N} - 18 T^{N - 1} + 9 T^{N - 2} - 2 T^{N - 3}}{6 δ_{r}} = T_{o} \end{matrix}

(75)

IV. Solving Linear Algebraic Equations: By solving a final system of linear algebraic equations, the temperature value is calculated for each points in the disk grid.

Navier equations: Radial and tangential direction

Next, the Navier equations for the radial and tangential directions are considered. In equations (27) and (28), the derivatives of displacement are replaced by equations (37) and (38). This results in the differential equations being expressed in finite difference form as equations (76) and (77):

\begin{matrix} \frac{- {u_{r}}^{i + 2} + 16 {u_{r}}^{i + 1} - 30 {u_{r}}^{i} + 16 {u_{r}}^{i - 1} - {u_{r}}^{i - 2}}{{δ_{r}}^{2}} \\ + (\frac{(1 + m_{1} + m_{5})}{r^{i}}) \frac{- {u_{r}}^{i + 2} + 8 {u_{r}}^{i + 1} - 8 {u_{r}}^{i - 1} + {u_{r}}^{i - 2}}{12 δ_{r}} \\ + (\frac{(v (m_{1} + m_{5}) - 1)}{{(r^{i})}^{2}}) {u_{r}}^{i} = (1 + v) α_{o} {(\frac{r^{i}}{r_{o}})}^{m_{3}} T^{i} (\frac{m_{1} + m_{3} + m_{5}}{r^{i}}) \\ + (1 + v) α_{o} {(\frac{r^{i}}{r_{o}})}^{m_{3}} \frac{- T^{i + 2} + 8 T^{i + 1} - 8 T^{i - 1} + T^{i - 2}}{12 δ_{r}} \\ - \frac{(1 - v^{2}) ρ_{o} {(\frac{r^{i}}{r_{o}})}^{m_{2}} ω_{o}^{2} {(e^{- λ t})}^{2} r^{i}}{E_{o} {(\frac{r^{i}}{r_{o}})}^{m_{1}}} \end{matrix}

(76)

\begin{matrix} \frac{- {u_{θ}}^{i + 2} + 16 {u_{θ}}^{i + 1} - 30 {u_{θ}}^{i} + 16 {u_{θ}}^{i - 1} - {u_{θ}}^{i - 2}}{{δ_{r}}^{2}} \\ + (\frac{(1 + m_{1} + m_{5})}{r^{i}}) \frac{- {u_{θ}}^{i + 2} + 8 {u_{θ}}^{i + 1} - 8 {u_{θ}}^{i - 1} + {u_{θ}}^{i - 2}}{12 δ_{r}} \\ - (\frac{(1 + m_{1} + m_{5})}{{(r^{i})}^{2}}) {u_{θ}}^{i} = \frac{2 (1 + v) ρ_{0} {(\frac{r^{i}}{r_{o}})}^{m_{2}} (λ ω_{0} e^{- λ t}) r^{i}}{E_{0} {(\frac{r^{i}}{r_{o}})}^{m_{1}}} \end{matrix}

(77)

The finite difference form of these differential equations is written for internal grid points. The boundary condition for rotating disk at the inner and outer surfaces in the radial and tangential direction are represented in the finite difference form as shown in Table 4.

Table 4.

Boundary conditions of the disk in finite difference form.

Cylindrical direction	FD form	Inner surface: under constraint outer surface: Free	Inner surface: under constraint outer surface: guided
Radial direction	Inner surface: $i = 1$	$\begin{matrix} u_{r} (r) = 0 \\ {u_{r}}^{1} = 0 \end{matrix}$	$\begin{matrix} u_{r} (r) = 0 \\ {u_{r}}^{1} = 0 \end{matrix}$
	Outer surface: $i = N$	$\begin{matrix} σ_{r} (r) = 0 \\ \frac{11 {u_{r}}^{N} - 18 {u_{r}}^{N - 1} + 9 {u_{r}}^{N - 2} - 2 {u_{r}}^{N - 3}}{6 δ_{r}} \\ + \frac{v}{r^{N}} {u_{r}}^{N} - (1 + v) α (r^{N}) T (r^{N}) = 0 \end{matrix}$	$\begin{matrix} u_{r} (r) = 0 \\ {u_{r}}^{N} = 0 \end{matrix}$
Tangential direction	Inner surface: $i = 1$	$\begin{matrix} u_{θ} (r) = 0 \\ u_{θ}^{1} = 0 \end{matrix}$	$\begin{matrix} u_{θ} (r) = 0 \\ u_{θ}^{1} = 0 \end{matrix}$
	Outer surface: $i = N$	$\begin{matrix} τ_{r θ} (r) = 0 \\ \frac{11 {u_{θ}}^{N} - 18 {u_{θ}}^{N - 1} + 9 {u_{θ}}^{N - 2} - 2 {u_{θ}}^{N - 3}}{6 δ_{r}} - \frac{u_{θ}^{N}}{r^{N}} = 0 \end{matrix}$	$\begin{matrix} u_{θ} (r) = 0 \\ u_{θ}^{N} = 0 \end{matrix}$

The radial and tangential displacement of the rotating disk for each grid point is determined by solving the system of equations. Using these displacements, the strain and stress values in each of these points can be calculated.

Numerical analysis

In this section, we verify the results of both heat transfer and Navier equations in radial and tangential directions using analytical HPM. We then provide numerical analysis examples based on the parameters of the problem. The values of the parameters in Table 3 are listed in Table 5, based on the geometric and mechanical properties of aluminum and zirconia corresponding to the inner and outer surfaces of the disk.

Table 5.

The grading index of FGM rotating disk made of aluminum-zirconia.

$m_{1}$	$m_{2}$	$m_{3}$	$m_{4}$
0.4776	0.4662	−0.5175	−2.8887

Rotating disks are mechanical parts that are used in a wide range of mechanical devices and equipment such as gears, turbine rotors, flywheels, shrink fits, etc. In this equipment, the inner surface of the disk is at environment temperature ( $T_{i} = 25^{°} C$ ) and the outer surface in contact with the boundary conditions such as fluids has a higher temperature ( $T_{o} = 125^{°} C$ ). This higher temperature on the outer surface can be constant (thermal condition A: $T (r) = T_{i}$ ) or include a temperature gradient (thermal condition B: $\frac{d}{dr} T (r) = T_{o}$ ).

Verification of HPM

The main objective of this section is to demonstrate the ability of the HPM to manage non-uniform FGM rotating disks under time-dependent mechanical and thermal loading conditions. The results obtained from solving the heat transfer equation (equation (54)) and the Navier equations (equations (63) and (72)) by HPM are verified with the help of HFDM. In HFDM, a mesh sensitivity analysis was performed to ensure that the results are independent of the meshing size. The number of points in the finite difference gridding of the disk was determined $N = 401$ for the aluminum disk with uniform thickness and material properties and free-free boundary conditions, as shown in Table 6. This type of disk has the highest level of stress compared to FGM disks with variable thickness and material properties. Therefore, this converged number can be used for FGM disks with all types of boundary conditions in this paper.

Table 6.

Convergence test of for HFDM along the radial direction of the rotating disk.

The number of points in the finite difference grid	Von Mises equivalent stress in $r = r_{i}$ (MPa)
$N = 401$	300.06
$N = 301$	300.05
$N = 201$	300.004
$N = 101$	299.53
$N = 51$	298.86

The results of solving the steady-state heat transfer equation at different radius of the disk for the thickness parameter $m_{5} = - 0.5$ are given in Table 7, based on this convergence test. As shown, the results of HPM and HFDM have converged to the same values with acceptable accuracy. For the thermal boundary condition type $A$ (equation (8)), the results of the two methods are very close. However, for the boundary condition of the type $B$ with the temperature gradient, the analytical solution provides more accurate results. This is because in HFDM, existing derivatives are approximated with difference functions, which shows the advantage of analytical solutions over numerical solutions. The temperature distribution of the disk for the two thermal boundary conditions is shown in Figure 8, verified by HPM. As seen in this figure, the disk with the temperature gradient on the outer surface (thermal condition $B$ ) experiences more uniform temperature changes than the thermal condition $A$ . Therefore, it is expected that the displacement and stress level in the disk with the thermal condition $B$ will be lower.

Table 7.

Comparison of the results for the heat transfer equation by HPM and HFDM: $m_{5} = - 0.5$ .

Disk radius (m)	Temperature based HFDM $(^{0} C$ )		Temperature based HPM $(^{0} C$ )		Percentage error HPM to HFDM	Percentage error HPM to HFDM
	Thermal conditions
	A	B	A	B	A	B
0.2	25	25	25	25	0	0
0.4	29.0716	26.4990	29.0717	26.4955	0.00034	0.013
0.6	42.3562	31.3900	42.3564	31.3749	0.00047	0.048
0.8	71.7181	42.1703	71.7182	42.1593	0.00013	0.026
1	125	61.7371	125	61.7294	0	0.012

Figure 8.

Temperature distribution for two types of thermal boundary conditions and different thickness parameter $m_{5}$ .

Table 8 presents the results obtained from two methods, HPM and HFDM, for Navier equations in radial and tangential directions at different radius of the disk. It is important to note that there is no exact solution for rotating disks with variable specifications. The disk is subjected to two different boundary conditions on the inner and outer surfaces, and two thermal boundary conditions can be considered separately for each type. For verification, it is assumed that the disk has constrained-free boundary conditions on the inner and outer surfaces. The angular velocity and acceleration are considered time-dependent, and the rotating disk is in the start process. Hence, parameters in equations (13) and (14) is considered as: $ω = ω_{0} e^{- λ t} |_{t = 1 s, λ = - 0.5} \to ω = 164.87 \frac{r a d}{s}$ and $α = - λ ω_{0} e^{- λ t} |_{t = 1 s, λ = - 0.5} = 82.43 \frac{r a d}{s^{2}}$ . As known, the results demonstrate the ability of these techniques to solve the governing equation of the rotating disk at radial and tangential directions in its general form. In HPM, two iterations lead to high precision of the solution for the governing equation of the problem. This solution is in good agreement with the results of the HFDM. The HPM is much simpler and more straightforward than HFDM or other analytical methods. It doesn’t require specific algorithms, complex calculations, or large computational capacity.

Table 8.

Comparison of the results for the Navier equations in the rotating disk by two methods of HPM and HFDM for the thickness parameter ( $m_{5} = - 0.5$ ) and ( $λ = - 0.5, t = 1 s$ ).

Disk radius (m)	Cylindrical direction	Displacement based FDM (m)		Displacement based HPM (m)		Percentage error HPM to HFDM	Percentage error HPM to HFDM
		Thermal conditions
		A	B	A	B	A	B
0.2	Radial direction	0	0	0	0	0	0
	Tangential direction	0		0		0
0.4	Radial direction	0.0002486	0.0002193	0.0002496	0.0002196	0.371	0.123
	Tangential direction	$7.38 * 10^{(- 6)}$		$7.3780 * 10^{(- 6)}$		0.0361
0.6	Radial direction	0.0004313	0.0003640	0.0004324	0.0003638	0.258	0.0617
	Tangential direction	$1.3039 * 10^{(- 5)}$		$1.3034 * 10^{(- 5)}$		0.0363
0.8	Radial direction	0.0006220	0.0004890	0.0006230	0.0004875	0.172	0.290
	Tangential direction	$1.8140 * 10^{(- 5)}$		$1.8133 * 10^{(- 5)}$		0.0366
1	Radial direction	0.0008594	0.0006120	0.0008603	0.0006085	0.108	0.581
	Tangential direction	$2.2879 * 10^{(- 5)}$		$2.2871 * 10^{(- 5)}$		0.0370

Example 1

By confirming the results of solving the governing equations using HPM, a parametric analysis was conducted in Figure 9 to study the distributions of radial, tangential, and shear stress for various thickness parameters $m_{5}$ , constant angular velocity function $λ = - 0.5$ , time $t = 2 s$ , and thermal boundary conditions $A$ and $B$ . According to equations (3) and (4), the angular velocity and acceleration were estimated as ( $ω = 271.82 \frac{r a d}{s}$ ) and ( $α = 135.91 \frac{r a d}{s^{2}}$ ), respectively. This position of the disk is equivalent to 2 s after the start process. The boundary condition of the disk is considered constrained-free for the inner and outer surfaces of the disk. As shown in Figure 9, this boundary condition is well satisfied.

Figure 9.

Distribution of stresses for two types of thermal boundary conditions and different thickness parameter $m_{5}$ , $λ = - 0.5, t = 2 s$ : (a) radial and tangential stress and (b) shear stress.

The maximum location of radial stress is at the inner surface of the disk, and as the thickness parameter increases, it moves toward the middle radius of the disk. For tangential stress, the middle radius of the disk has the maximum values of stress. In the thermal condition $B$ , the tangential stress in the outer surface of the disk has larger values than in the thermal condition $A$ . The levels of shear stress compared to other stresses are low, and thermal boundary conditions do not affect its values. As the thickness parameter increases, the level of these stresses decreases for both thermal boundary conditions. In general, the stresses for the thermal condition $B$ are lower than for the condition $A$ . This issue follows the results stated in Figure 8, where the disk in the gradient thermal condition $B$ experiences uniform temperature changes along the radius of the disk compared to the thermal condition $A$ .

Example 2

In this example, the distribution of equivalent von Mises stress for a rotating FGM disk is demonstrated for various thickness parameters $m_{5}$ , constant angular velocity function $λ = - 0.5$ , time $t = 2 s$ , and thermal boundary condition $A$ and $B$ . The angular velocity, acceleration, and position of the disk are estimated in Example 1. The von Mises stress distribution along the disk radius is important since plastic deformations in rotating disks are not desirable and start from the radius where the yield stress reaches its equivalent von Mises stress. The findings are illustrated based on the constrained-free and constrained-guided boundary conditions in Figure 10.

Figure 10.

Distribution of von Mises stress for two types of thermal boundary conditions and different thickness parameter $m_{5}$ , $λ = - 0.5, t = 2 s$ : (a) constrained-free and (b) constrained-guided.

Constrained-free boundary condition: The maximum values of von Mises stress for different thickness parameters $m_{5}$ occur at the inner radius of the disk ( $r = r_{i}$ ), and as this parameter increases, it moves toward the middle radius of the disk. For radius close to the inner surface of the disk, the thermal condition $A$ produces a higher equivalent stress in the disk, while for radius close to the outer surface of the disk, a thermal condition $B$ exerts more equivalent stress on the disk.

Constrained-guided boundary condition: The maximum values of equivalent stress for all thickness parameters $m_{5}$ occur at the outer radius of the disk ( $r = r_{o}$ ). In this case, the combination of radial, tangential, and shear stress on the disk’s outer surface is such that as the thickness parameter increases, the value of maximum von Mises stress increases.

Approximately, for all disk radius, the equivalent von Mises stress for the thermal condition $A$ is higher than the thermal condition $B$ .

Example 3

In this example, the yield stress function for the FGM rotating disk is calculated. The elastic limit occurs when the yield stress in the disk reaches the equivalent von Mises stress according to equation (11). Considering the ceramic volume fraction function as ( $V_{c} (r) = {(\frac{r}{r_{o}})}^{m_{1}}$ ), equation (11) is rewritten in the following form:

σ_{y} (r) = σ_{y m} ((1 - {(r / r_{o})}^{m_{1}}) + \frac{q + E_{m}}{q + E_{c}} \frac{E_{c}}{E_{m}} {(r / r_{o})}^{m_{1}})

(78)

In this relation, the elastic modulus of aluminum $E_{m}$ , the elastic modulus of zirconia $E_{c}$ , and the coefficient $m_{1}$ according to Table 3 are placed. Considering relation ( $σ_{y} (r) = σ_{eq} (r)$ ) as a criterion for the initiation of plastic deformations, it is possible to find the angular velocity and acceleration of the elastic limit in the FGM disk according to equation (79). There are two unknowns in this relation, the radius and the angular velocity.

\begin{matrix} \sqrt{σ_{r}^{2} + σ_{θ}^{2} - σ_{r} σ_{θ} + 3 τ_{r θ}^{2}} \\ = σ_{ym} ((1 - {(\frac{r}{r_{o}})}^{m_{1}}) + \frac{q + E_{m}}{q + E_{c}} \frac{E_{c}}{E_{m}} {(\frac{r}{r_{o}})}^{m_{1}}) \end{matrix}

(79)

For better understanding, the distribution of the ratio ( $\frac{σ_{e q} (r)}{σ_{y} (r)}$ ) for various thickness parameters $m_{5}$ , constant in angular velocity function $λ = - 0.5$ , time $t = 3.5 s$ , and thermal boundary condition $A$ and $B$ , are plotted in Figure 11. The angular velocity and acceleration are estimated as ( $ω = 575.46 \frac{r a d}{s}$ and $α = 287.73 \frac{r a d}{s^{2}}$ ). This position of the disk is equivalent to 3.5 s after the start process. As it is clear, for some radius of the disk, plastic deformation is initiated. It means that in these radius, the von Mises equivalent stress has exceeded the yield stress of the metal component of the disk in Table 2 (300 MPa). For constrained-free boundary condition and both thermal conditions, for different values of the thickness parameter, plastic deformation is started and most of the disk radius has yielded. For constrained-guided boundary condition and both thermal conditions, some radius near the outer surface of the disk are yielded, but most disk radius still have elastic deformation.

Figure 11.

Distribution of von Mises stress to yeild stress for two types of thermal boundary conditions and different thickness parameter $m_{5}$ , $λ = - 0.5, t = 3.5 s$ : (a) constrained-free and (b) constrained-guided.

Example 4

The example investigates the effect of the thickness parameter $m_{5}$ , boundary thermal condition, and boundary conditions of the rotating disk on the values of elastic limit angular velocity ( $ω_{e}$ ) and acceleration ( $α_{e}$ ). Tables 9 to 11 record the radius of plastic initiation, the time to start these deformations, and elastic limit angular velocity and acceleration for these parameters. It is assumed that the disk is in the start process, and the constant in angular velocity is negative as, $λ = - 0.2, - 0.35, - 0.5$ . At the time $t$ , when the von Mises equivalent stress exceeds the yield stress, the disks are yielded. These tables show that the time at which the disks are yielded is increased by decreasing the constant $λ$ in all thermal boundary and boundary conditions of the disk. Additionally, the influence of the constant $λ$ on the values of elastic limit angular acceleration is greater than its influence on the elastic limit angular velocity.

Table 9.

Results for elastic limit angular velocity and acceleration for different values of thickness parameter $m_{5}$ and $λ = - 0.2$ .

Thickness parameter ( $m_{5}$ )	Boundary condition	Radius of plastic initiation		Time $t (s)$		Elastic angular velocity $ω_{e} (\frac{r a d}{s})$		Elastic angular acceleration $α_{e} (\frac{r a d}{s^{2}})$
		Thermal boundarycondition
		$A$	$B$	$A$	$B$	$A$	$B$	$A$	$B$
0	Constrained-free	0.2	0.2	5.97	6.15	330.01	342.12	66.00	68.42
	Constrained-guided	1	1	6.80	7.60	389.61	457.22	77.92	91.44
−0.5	Constrained-free	0.2	0.2	7.36	7.47	435.79	445.48	87.15	89.09
	Constrained-guided	1	1	6.38	7.36	358.22	435.79	71.64	87.15
−1	Constrained-free	0.569	0.57	8.27	8.32	522.78	528.03	104.55	105.60
	Constrained-guided	1	1	5.91	7.12	326.08	415.37	65.21	83.07

Table 10.

Results for elastic limit angular velocity and acceleration for different values of thickness parameter $m_{5}$ and $λ = - 0.35$ .

Thickness parameter ( $m_{5}$ )	Boundary condition	Radius of plastic initiation		Time $t (s)$		Elastic angular velocity $ω_{e} (\frac{r a d}{s})$		Elastic angular acceleration $α_{e} (\frac{r a d}{s^{2}})$
		Thermal boundary condition
		$A$	$B$	$A$	$B$	$A$	$B$	$A$	$B$
0	Constrained-free	0.2	0.2	3.41	3.52	329.86	342.80	115.45	119.98
	Constrained-guided	1	1	3.89	4.34	390.20	456.76	136.57	159.86
−0.5	Constrained-free	0.2	0.2	4.21	4.27	436.44	445.71	152.75	155.99
	Constrained-guided	1	1	3.65	4.21	358.76	436.44	125.56	152.75
−1	Constrained-free	0.569	0.57	4.72	4.76	521.74	529.09	182.60	185.18
	Constrained-guided	1	1	3.38	4.07	326.41	415.57	114.24	145.45

Table 11.

Results for elastic limit angular velocity and acceleration for different values of thickness parameter $m_{5}$ and $λ = - 0.5$ .

Thickness parameter ( $m_{5}$ )	Boundary condition	Radius of plastic initiation		Time $t (s)$		Elastic angular velocity $ω_{e} (\frac{r a d}{s})$		Elastic angular acceleration $α_{e} (\frac{r a d}{s^{2}})$
		Thermal boundary condition
		$A$	$B$	$A$	$B$	$A$	$B$	$A$	$B$
0	Constrained-free	0.2	0.2	2.38	2.46	328.70	342.12	115.45	171.06
	Constrained-guided	1	1	2.73	3.04	391.57	457.22	195.78	228.61
−0.5	Constrained-free	0.2	0.2	2.95	2.99	437.10	445.93	218.55	222.96
	Constrained-guided	1	1	2.55	2.95	357.87	437.10	178.93	218.55
−1	Constrained-free	0.569	0.57	3.31	3.33	523.30	523.30	261.65	261.65
	Constrained-guided	1	1	2.37	2.85	327.06	415.78	163.53	207.89

The values of elastic limit angular velocity and acceleration increase by increasing the thickness parameter $m_{5}$ in constrained-free conditions. However, as this parameter increases, the values of elastic limit angular velocity and acceleration decrease for the constrained-guided conditions. According to Figure 10, in this boundary condition, the maximum value of von Mises stress occurs on the outer surface of the disk, and a disk with a higher thickness parameter has a more equivalent von Mises stress. In terms of boundary conditions, the disk with the constrained-free condition has more elastic limit angular velocity and acceleration compared to the disk with the constrained-guided conditions. On the effect of the thermal boundary condition, a thermal condition $B$ due to a temperature gradient on the outer surface and uniform changes of temperature along the radial direction of the disk has higher values for the elastic limit angular velocity and acceleration. The presented points can be used to control the plastic deformations in the rotating disk during the accelerating process until the normal work condition with constant angular velocity is reached.

Example 5

The example investigates the effect of boundary thermal conditions of a rotating disk on the values of angular velocity ( $ω$ ) and acceleration ( $α$ ) during the deceleration process. Tables 12 to 14 present the results of the study, where the constrained-free disk boundary conditions are considered and the disk thickness parameter is assumed as $m_{5} = - 0.5$ . The study presents the values of maximum radial displacement and radial stress and their locations for different times in the disk at the stopping process. The constant in angular velocity is assumed to be positive as, $λ = 0.2, 0.35, 0.5$ . The values of angular velocity and acceleration at any time of the disk process are independent of the thermal boundary condition, according to equations (13) and (14).

Table 12.

Results for angular velocity and acceleration over time for thickness parameter $m_{5} = - 0.5$ and $λ = 0.2$ .

Time $t (s)$	angular velocity $ω (\frac{r a d}{s})$	angular acceleration $α (\frac{r a d}{s^{2}})$	Maximum value of $u_{r}$ ( $\times 10^{- 4} m$ )		Radial coordinate of maximum $u_{r}$ ( $m$ )		Maximum value of $σ_{r}$ ( $\times 10^{7} Pa$ )		Radial coordinate of maximum $σ_{r}$ ( $m$ )
			Thermal boundary condition
			$A$	$B$	$A$	$B$	$A$	$B$	$A$	$B$
0	100	−20	7.573	5.056	1	1	5.267	3.942	0.2	0.2
5	36	−7.35	7.060	4.543	1	1	3.916	2.590	0.2	0.2
10	13.53	−2.70	6.991	4.474	1	1	3.733	2.408	0.2	0.2
15	4.97	−0.99	6.981	4.465	1	1	3.708	2.383	0.2	0.2

Table 13.

Results for angular velocity and acceleration over time for thickness parameter $m_{5} = - 0.5$ and $λ = 0.35$ .

Time $t (s)$	angular velocity $ω (\frac{r a d}{s})$	angular acceleration $α (\frac{r a d}{s^{2}})$	Maximum value of $u_{r}$ ( $\times 10^{- 4} m$ )		Radial coordinate of maximum $u_{r}$ ( $m$ )		Maximum value of $σ_{r}$ ( $\times 10^{7} Pa$ )		Radial coordinate of maximum $σ_{r}$ ( $m$ )
			Thermal boundary condition
			$A$	$B$	$A$	$B$	$A$	$B$	$A$	$B$
0	100	−35	7.573	5.056	1	1	5.267	3.942	0.2	0.2
5	17.37	−6.08	6.998	4.481	1	1	3.751	2.426	0.2	0.2
10	3.01	−1.05	6.980	4.464	1	1	3.706	2.380	0.2	0.2
15	0.524	−0.183	6.980	4.463	1	1	3.704	2.379	0.2	0.2

Table 14.

Results for angular velocity and acceleration over time for thickness parameter $m_{5} = - 0.5$ and $λ = 0.5$ .

Time $t (s)$	angular velocity $ω (\frac{r a d}{s})$	angular acceleration $α (\frac{r a d}{s^{2}})$	Maximum value of $u_{r}$ ( $\times 10^{- 4} m$ )		Radial coordinate of maximum $u_{r}$ ( $m$ )		Maximum value of $σ_{r}$ ( $\times 10^{7} Pa$ )		Radial coordinate of maximum $σ_{r}$ ( $m$ )
			Thermal boundary condition
			$A$	$B$	$A$	$B$	$A$	$B$	$A$	$B$
0	100	−50	7.573	5.056	1	1	5.267	3.942	0.2	0.2
5	8.20	−4.10	6.984	4.467	1	1	3.715	2.389	0.2	0.2
10	0.673	−0.33	6.980	4.463	1	1	3.704	2.379	0.2	0.2
15	0.055	−0.027	6.980	4.463	1	1	3.704	2.379	0.2	0.2

The tables show that more time is required to stop the rotating disk as the constant $λ$ decreases. Therefore, the disk with the specifications of Table 14 will stop earlier than the other two disks. For all these disks, the maximum radial displacement occurs on the outer surface, and the maximum radial stress occurs on the inner surface of the disk. The values of these parameters for the thermal condition $B$ are lower than for the thermal condition $A$ . Based on the results of examples 4 and 5, the optimal value of the constant $λ$ can be determined according to the conditions of the working environment of the disk. This will ensure that plastic deformation does not occur at the start process, and the stop process occurs at the optimal time.

Conclusion

In this paper, an analytical model is presented to investigate the deformation-stress distribution in FGM rotating disks under time-dependent mechanical and thermal loading conditions. The model considers all the geometrical and mechanical properties of the rotating disk as variable in the radial direction. The heat transfer equation and the governing Navier equations in the radial and tangential direction of the rotating disk are solved using the HPM as the analytical method and verified using the HFDM as the numerical method. The results provided by HPM have good accuracy for both heat and stress analysis, and successfully handle the solutions of these equations simultaneously, as shown in Tables 7 and 8. Since the results of the two analytical and numerical methods are in good agreement, it can be predicted that the results are in good agreement with the experimental data regardless of laboratory errors.

The numerical analysis showed that the disk with a thermal boundary condition $B$ experiences less stress compared to the disk with a thermal boundary condition $A$ for various boundary conditions of the disk and different values of the thickness parameter. Therefore, by defining the appropriate temperature gradient on the outer surface of the disk as a thermal boundary condition $B$ , the level of thermal stresses can be controlled and optimized. In this paper, the level of thermal stresses by thermal condition $B$ has been reduced by 15%–20% compared to the thermal boundary condition $A$ for different thickness parameters $m_{5}$ .

The type of boundary conditions affects the distribution of equivalent von Mises stress. It was shown that the yielding along the radial direction of the rotating disk is not a fixed function, and the ratio of $\frac{σ_{e q}}{σ_{y}}$ along the radial direction of the disk is variable for different available parameters. In the disk with constrained-free boundary conditions, with the increase of the thickness parameter, the stress level decreased, and the angular velocity and acceleration increased. However, for the disk with constrained-guided boundary conditions, the behavior of the disk was different, and with the increase of the thickness parameter, the stress level increased, and then the angular velocity and acceleration of the disk decreased.

In the start process of the disk, the time in which plastic deformation is introduced can be improved by considering the smaller value for the constant $λ$ . In the stop process of the disk, as the constant $λ$ is smaller, more time is required to stop the disk. Finally, it can be concluded that the proposed analytical model can successfully handle the FGM rotating disk problem in the general form under time-dependent mechanical and thermal loading. The presented model leads to time and cost savings in handling complicated cases without the need to use finite element analysis software, which is commercially available.

As the future works, this model could be expanded to examine more complicated problems of combined loading cases of rotating disk in gears, turbine rotors, flywheels. In this equipment, any form of functions for variable parameters along with creep analysis can be considered. Laboratory tests and the use of experimental data are also the future part of research in this topic.

Footnotes

Appendix

Notation

r	Radius of disk
r _i	Inner radius of disk
r _o	Outer radius of disk
h ₀	Thickness of disk at r=r_o
E	Young’s modulus (GPa)
ν	Poisson coefficient
ρ	Density, $ρ (\frac{k g}{m^{3}})$
α	Thermal expansion coefficient, $α (\frac{10^{- 6}}{° C})$
k	Thermal conductivity coefficient, $k (\frac{W}{m^{°} C})$
σ₀	Yield stress (MPa)
E _t	Tangent modulus (GPa)
q	Constant (GPa)
ω ₀	Angular velocity constant (rad/s)
λ	Angular velocity constant (1/s)
E _e	Reference values for elastic modulus
ρ _e	Reference values for density
α _e	Reference values for thermal expansion coefficient
k _e	Reference values for thermal conductivity coefficient
m ₁	Grading index of elastic modulus
m ₂	Grading index of density
m ₃	Grading index of thermal expansion coefficient
m ₄	Grading index of thermal conductivity coefficient
m ₅	Grading index of thickness function
C _ij _(i,j=1,2)	Constant related to heat transfer function
ξ _ij _(i,j=1,2)	Constant related to heat transfer function
T _i	Temperature on the inner surface
T ₀	Temperature on the outer surface
σ _ym	Yield stress in the metal
V _c _(r)	Volume fraction function
ω	Angular velocity
ω _e	Elastic limit angular velocity
α	Angular acceleration
α _e	Elastic limit angular acceleration
σ_r	Radial stress
σ_θ	Tangential stress
τ _rθ	Shear stress
σ_eq	Equivalent von Mises stress
t	Time (s)
ε _r	Radial strain
ε _θ	Tangential strain
γ_rθ	Shear strain
ε^T	Thermal strain
u _r	Radial displacement
u _θ	Tangential displacement
u	Unknown function in HPM

Handling Editor: Sharmili Pandian

Declaration of conflicting interests

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Sanaz Jafari

References

Gamer

. Elastic-plastic deformation of the rotating solid disk. Arch Appl Mech 1984; 54: 345–354.

Eraslan

Orcan

. Elastic–plastic deformation of a rotating solid disk of exponentially varying thickness. Mech Mater 2002; 34: 423–432.

Eraslan

Orcan

. On the rotating elastic–plastic solid disks of variable thickness having concave profiles. Int J Mech Sci 2002; 44: 1445–1466.

Eraslan

Orcan

. Elastic–plastic deformations of rotating variable thickness annular disks with free, pressurized and radially restricted boundary conditions. Int J Mech Sci 2003; 45: 643–667.

You

Zhang

, et al. On rotating circular disks with varying material properties. Z Angew Math Phys 2007; 58: 1068–1084.

Kordkheili

SAH

Naghdabadi

. Thermoelastic analysis of a functionally graded rotating disk. Compos Struct 2007; 79: 508–516.

Bayat

Saleem

Sahari

, et al. Mechanical and thermal stresses in a functionally graded rotating disk with variable thickness due to radially symmetry loads. Int J Press Vessels Piping 2009; 86: 357–372.

Hojjati

Jafari

. Variational iteration solution of elastic non uniform thickness and density rotating disks. Far East J Appl Math 2007; 29: 185–200.

Hojjati

Hassani

. Theoretical and numerical analyses of rotating discs of non-uniform thickness and density. Int J Press Vessels Piping 2008; 85: 694–700.

10.

Hojjati

Jafari

. Semi exact solution of elastic non uniform thickness and density rotating disks by homotopy perturbation and Adomian’s decomposition methods. Part I: elastic solution. Int J Press Vessels Piping 2008; 85: 871–878.

11.

Hojjati

Jafari

. Semi-exact solution of non-uniform thickness and density rotating disks. Part II: elastic strain hardening solution. Int J Press Vessels Piping 2009; 86: 307–318.

12.

Nayak

Bhowmick

NathSaha

. Elasto-plastic analysis of thermo-mechanically loaded functionally graded disks by an iterative variational method. Eng Sci Technol 2020; 23: 42–64.

13.

Jafari

Hojjati

Fathi

. Classical and modern optimization methods in minimum weight design of elastic rotating disk with variable thickness and density. Int J Press Vessels Piping 2012; 92: 41–47.

14.

Alashti

Jafari

. The effect of ductile damage on plastic behavior of a rotating disk with variable thickness subjected to mechanical loading. Sci Iran B 2016; 23: 174–193.

15.

Zheng

Bahaloo

Mousanezhad

, et al. Stress analysis in functionally graded rotating disks with non-uniform thickness and variable angular velocity. Int J Mech Sci 2016; 119: 283–293.

16.

Salehian

Shahriari

Yousefi

. Thermo-elastic analysis of a functionally graded rotating hollow circular disk with variable thickness and angular speed. J Braz Soc Mech Sci Eng 2018; 41: 33–45.

17.

Jafari

. Analytical modeling of elastic limit angular velocity in a rotating disk of a functionally graded material under mechanical-thermal loading conditions. AmirKabir J Mech Eng 2021; 53: 2589–2612.

18.

Jafari

. Elastic limit angular velocity and acceleration investigation in non-uniform rotating disk under time-dependent mechanical loading. J Appl Comput Mech 2022; 8: 791–808.

19.

. Homotopy perturbation technique. Comput Methods Appl Mech Eng 1992; 178: 257–262.

20.

. Homotopy perturbation method: a new nonlinear analytical technique. Appl Math Comput 2003; 135: 73–80.

21.

. Asymptotology by homotopy perturbation method. Appl Math Comput 2004; 6: 156–591.

22.

. Limit cycle and bifurcation of nonlinear problems. Chaos Solitons Fractals 2005; 26: 827–833.

23.

. Homotopy perturbation method for bifurcation of nonlinear problems. Int J Nonlinear Sci Numer Simul 2005; 6: 207–208.

24.

Nakmura

. Applied numerical methods with software. Englewood Cliffs, NJ: Prentice-Hall International Inc, 1991.

25.

Ashok

Singh

Bhadauria

. Finite difference formulae for unequal sub-intervals using Lagrange’s interpolation formula. Int J Math Anal 2009; 3: 815–827.