Stress Analysis of Nonhomogeneous Rotating Disc with Arbitrarily Variable Thickness Using Finite Element Method

Rotating discs with variable thickness and nonhomogeneous material properties are frequently used in industrial applications. The nonhomogenity of material properties is often caused by temperature change throughout the disc. The governing differential equation presenting this problem contains many variable coefficients so that no possible analytical closed form solution for this problem. Many numerical approaches have been proposed to obtain the solution. However, in this study the Finite Element Method (FEM), which presents a powerful tool for solving such a problem, is used. Thus, a turbine disc modeled by using ax symmetric finite elements was analyzed. But, in order to avoid inaccuracy of the stress calculation quite fine meshing is implemented. The analysis showed that maximum displacement occurs at the boundary of the disc, either at the outer or inner boundary, depending on the loadings. The maximum radial stress occurs at an area in the middle of the disc which has the smallest thickness. In this study, rotational blade load was shown to give the largest contribution to the total displacement and stress. Also, the radial displacement and stress in a disc with variable thickness are found to be affected by the contour of the thickness variation. In general, the results obtained show excellent agreement with the published works.


INTRODUCTION
Rotating discs have many practical engineering applications such as in steam and gas turbine discs, turbo generators, internal combustion engines, casting ship propellers, turbojet engines, reciprocating and centrifugal compressors just to mention a few.Brake disk can be an example of solid rotating disk where only body force is involved.
In reality, the thickness of the disc is frequently not constant, such as turbine discs, due to economic considerations and in order to improve mechanical performance.Turbine discs are usually thick near their hub and taper down to a smaller thickness toward the periphery, since a constant thickness pronounces stress concentration near the center, even for solid disc.Furthermore, it has been shown that the stresses in variable thickness rotating annular and solid discs are much lower than those in constant thickness discs at the same angular velocity.
In many applications, the disc is working under high temperature which presents thermal loading.Beside the fact that, the temperature throughout the disc is usually not constant, i.e., there is temperature gradient present throughout the disc.This temperature gradientusually resulted in changes in material properties throughout the disc and therotating disc becomes nonhomogeneous in the radial direction.
Generally, there are two approaches for the solution of rotating discs, namely, theoretical and numerical methods.For homogeneous rotating disc with constant thickness, closed-form analytical solution is available.However, for nonhomogeneous rotating disc with variable thickness, analytical solution is not possible to obtain.Hence, many numerical attempts has been presented to solve such a problem.Timoshenko and Goodier (1970) was the first to obtain a closed form solution for homogeneous constant-thickness rotating discs without any temperature gradient.Stodola (1927) obtained an analytical solution of the problem of homogeneous rotating discs with hyperbolic profile thickness.Den Hartog (1987) reported the closed form formula for homogeneous rotating discs with constant and hyperbolic thickness under several mechanical loadings.Boresi and Richard (2003) included thermal stress in the closed form formula. Gamer (1985) achieved a good adaptive numerical solution for the constant thickness solid discs with a linear hardening material.Recently, Sharma et al. (2011) conducted analysis of stresses and strains in a rotating homogeneous thermoplastic circular disk using FEM.Guven (1997Guven ( , 1998) ) obtained an analytical solution of the problem of rotating variable thickness discs with rigid inclusion for the elastic-plastic (Guven, 1997) and fully plastic state (Guven, 1997).Leopold (1984) calculated elastic stress distributions in rotating discs withvariable thickness by using a semi-graphical method.Later, semi-analytical analysis of Functionally Graded Materials (FGM) rotating discs with variable thickness was proposed by Bayat et al. (2008).Hojjati and Jaffari (2007) proposed Variation Iteration Method (VIM) to obtain the elastic analysis of non-uniform thickness and density rotating discs subjected to only centrifugal loadings.However, You et al. (2000) proposed numerical analysis using Runge-Kutta method compared to FEM for elastic-plastic rotating discs with arbitrary variable thickness and density.Likewise, Hojjati and Hassani (2008) proposed Variablematerial Properties (VMP) method and numerical analysis using Runge Kutta's method compared to FEMfor rotating discs of non-uniform thickness and density.Recently, they also proposed semi-exact solution for thermomechanical analysis of FGM elastic-strain hardening rotating discs (Hassani et al., 2012).Furthermore, to solve nonhomogeneous rotating disc with variable thickness, Jahed et al. (2005) proposed discretization of the disc into a finite number of rings, where each ring has constant thickness and material properties, but different rings may have different thickness and material properties.
In this study, therefore, the problem of nonhomogeneous rotating disc with arbitrarily variable thicknessis addressed and other related issues are discussed and presented.It is demonstrated that the analytical solution for such a problem is not possible to be obtained.Consequently, to solve this problem the FE technique is used.The complete details of the FE formulation and analysis of the problemalong with numerical examples to verify the solution are presented in the next sections.

METHODOLOGY Analytical solution for elastic nonhomogeneous disc with variable thickness under mechanical and thermal loading:
The problem can be considered as plane stress with variable thickness.Hence, σ z = 0, as this is a static problem, the solution has to satisfy equilibrium, compatibility and constituve law of the material properties.
Equilibrium: Assuming that the weight of the disc is neglected, equilibrium of infinitesimal element of the disc, as shown in Fig. 1, in the radial direction is: Because the problem is axisymmetric, stress components in tangential direction vanish: By simplification and taking sin (dθ/2) = dθ/2 due to small angle Ө, we obtained: Compatibility: Displacaments occurs in the disc is shown in Fig. 2. The compatibility ends up with straindisplacement relation as follows: Properties of material, which are expressed by constitutive law of material: Because the problem is elastic, Hooke's law is used to get the stress-strain relation.
Due to only mechanical loading and σ z = 0: ( ) ( ) Because there is also thermal loading; then, and Hence, the total strains become: ( ) ( ) Re-arranging Eq. ( 10) and ( 11), we get: and Subtituting Eq. ( 13) into (12) to get σ r : ( ) Now subtituting Eq. ( 12) in to (13) to get σ t : ( ) Subtituting Eq. ( 4) and ( 5) into ( 14) we get: Similarly, substituting Eq. ( 4) and ( 5) into (15) we get: Now subtituting Eq. ( 29) and ( 31) into (2), which is the equilibrium equation in radial direction, we obtained a second order Ordinary Differential Equation (ODE) of the form: where, The coefficients of the ODE contain variable parameters: , ( Because the coefficients of the ODE contain many variable parameters, there are no exact solutions for this ODE.However, as stated earlier, there have been numerous numerical approaches to solve such a problem.But in this study, the FE approach is utilized.The details of the FE analysis and solution are presented below.
Fem solution for elastic nonhomogeneous disc with variable thickness under mechanical and thermal loading: FEM is one of the most successful and dominant numerical method in the last century.It is extensively used in modeling and simulation of engineering and science due to its versatility for complex geometries of solids and structures and its flexibility for many non-linear problems.
Most discs with variable thickness in the applications are ax symmetric.For such a case, ax symmetric element is the most economical but adequate to use in the finite element analysis.For any other cases in which the disc is not ax symmetric and therefore not adequate to be modeled by axisymmetric element, cyclic element and 3D solid element can be used.
The axisymmetric symmetric element has 2 translational degrees of freedom per node.Using cylindrical coordinate system where z is the axis, r is the radius and Ө is the circumference angle, the stresses in the axisymmetric problem is shown in Fig. 3.
Internal and external pressure working at the inner and outer surface of the disc are surface forces.The element surface force vector is given by:

∫∫
Centrifugal forces due to rotational speed of the disc are body forces.The element body force vector is given by:

∫∫
The element thermal forces working on the disc are given by: The total element forces {f} are sum of th surface forces, the element body forces and the element thermal forces.The global force vectors {F} are obtained by assembling all the element forces.The element stiffness matrix is given by: Assembling all of the element stiffness matrices, the global stiffness matrix [K] is obtained.The global displacement vectors are given by {F} = {K}{d}.The displacement function vectors are given by: The global strain vectors are given by: Finally, the global stress vector are given by: Numerical example: In this example, gas turbine disc was analyzed.Noncommercial Workbench was used to conduct the App.Sci. Eng. Technol., 7(15): 3114-3125, 2014 3117 Stresses in the axisymmetric problem  The steady state temperature at the outer boundary is equal to the temperature of high temperature gas.Based on the data given by Barack and Domas (1976), Liu et al. (2002), Claudio et al. (2002), Jahed (2005), Maruthi et al. (2012) and Elhefny and (2012), this temperature has a range of 550 this study, a value of 800°C is used.Throughout the disc, the temperature gradually decreases as the radius gets closer to the inner boundary.The temperature throughout the disc never reaches cooling system is applied inside the disc.Based on data published by Barack andDomas (1976) (2005) and recently by Maruthi temperature at the inner boundary has a range of 450 500°C.In this study, a value of 500 Following the aforementioned boundary temperatures along with assumption that t changethroughout the disc is linear, the following relation is used to express the temperature change: T (r) = 3 r + 324; r in mm The temperature distribution is shown in Fig. 5 Due to the temperature variation, poisson ratio and thermal coefficient vary throughout the disc as follows: A half cross sectional view of the disc with he thickness of the disc varies along the radial direction.The temperature varies along the radial direction as well.Due to the temperature variation, and thermal coefficient change as function of the temperature.Because the temperature changes as function of radius, then Young he thermal coefficient can also be expressed as function of radius.Because analytical solutions of such a problem is not available as already explained, finite element analysis was carried

Geometry and material properties:
The disc is geometry of the disc is shown in change of temperature, and thermal coefficient The density of 8220 kg/m 3 is The steady state temperature at the outer boundary is equal to the temperature of high temperature gas.Based on the data given by Barack and Domas (1976), .( 2002), Jahed et al. . (2012) and Elhefny and Guozhu (2012), this temperature has a range of 550-900°C.In C is used.Throughout the disc, the temperature gradually decreases as the radius gets closer to the inner boundary.The temperature throughout the disc never reaches the same value as a cooling system is applied inside the disc.Based on data published by Barack and Domas (1976) and Jahed et al. (2005) and recently by Maruthi et al. (2012), the temperature at the inner boundary has a range of 450a value of 500°C is used.Following the aforementioned boundary temperatures along with assumption that the temperature linear, the following relation is used to express the temperature change: (28) The temperature distribution is shown in Fig. 5.
Due to the temperature variation, Young Modulus, and thermal coefficient vary throughout As the temperature is a function of radius, Eq. ( 29) to ( 31) can be expressed as functions of radius, such that E, v and α are given by: The change of Young modulus, Poisson ratio and thermal coefficient is shown in Fig. 6 to 8.

Element type, loads and boundary conditions:
Axisymmetric model is used due to the axisymmetric geometry of the disc.Furthermore, because the model is symmetric in longitudinal direction, then half of the section area was used.
Shrink fit with the rotor shaft causes a radial force in outward direction on the inner surface of the disc.This force results in 300 MPa compressive pressures on the inner surface of the disc.The rotation of the disc causes a centrifugal body force in outward direction.The rotation velocity is 15,000 rpm = 1570 rad/s.The As boundary condition, frictionless support (i.e., rollers) is put on the disc side which functions as symmetry plane.Thus, no displacement perpendicular to this plane is allowed, yet displacement along the plane is allowed.The model is shown in Fig. 9.
Meshing: Coarse and refined meshing qualities were applied to the model as shown in Fig. 10.For both, nodes are only at the element corners.

RESULTS AND DISCUSSION
Case 1: The disc is subjected to shrink fit load only: By taking only the shrink fit load into account, the radial displacement and stress are shown in Fig. 11.The displacement and stress have outward direction due to the shrink fit load coming from the shaft.The maximum displacement as well as the maximum absolute value of the radial stress occurs at the inner surface of the disc.As the radius increases, the displacement and stress decrease until it reaches the minimum value on the outer surface of the disc.An abrupt increase of the stress distribution occurs at the smallest thickness of disc as stress concentration occurs there.The rounded-off value of the minimum stress is zero, which satisfies zero pressure boundary condition on the outer surface of the disc.
Case 2: The disc is subjected to rotational body load only: The centrifugal body force causes radial displacement and stress in outward direction, as shown in Fig. 12.It is shown that the displacement curve follows the thickness contour of the disc.Increasing thickness tends to result in increasing displacement.The maximum displacement occurs at the inner surface of the disc whereas the minimum occurs at the outer surface.The maximum radial stress occurs at the middle of the disc.
Case 3: The disc is subjected to rotational blade load only: Beside the rotational body load, the rotating attached blades cause additional centrifugal load working on the outer surface of the disc.This load causes displacement and stress in outward direction, as shown in Fig. 13.The maximum displacement occurs at the outer surface of the disc, whereas the minimum occurs near the hub of the disc.Large radial stress occurs as the radius increases.With the given geometry of the disc, the maximum stress occurs at the smallest thickness of the disc as stress concentration occurs there.

Case 4: The disc is subjected to thermal load only:
The thermal load causes radial displacement and stress as shown in Fig. 14.The displacement has negative values at most of the radius values.The negative values represent displacement in inward direction.Near the outer surface, the displacement has positive values, which means that the displacement has outward direction.Thus, due to the thermal load, the disc expands to inward and outward direction.The radial stress has outward direction.The rounded-off stress at the inner and outer surface of the disc is zero.The maximum stress occurs at the middle of the disc.
Case 5: The disc is subjected to combined loads: By combining all the loads, the radial displacament of the disc with refined meshingis shown in Fig. 15.Large radial displacement occurs at the inner and outer surface of the disc.The displacements at the inner and outer surface are 1.31 and 1.33 mm, respectively.The minimum displacement of 1.136 mm occurs at the middle of the disc.
The maximum values of the radial displacement and stress due to each load as well as combined loads are shown in Table 1.It is shown that the rotational blade load contributes the most to the total displacement and stress.To reduce its value, reducing the weight of the blades and/or the rotor speed can be proposed.The maximum displacement occurs at the outer surface, yet it is only slightly larger than that at the inner surface.Reducing the rotational blade load can be proposed to reduce the displacement at the outer surface, as excessive displacement at that part may cause rubbing against the stationary parts.Assuming that the shrink fit is fixed, reducing the rotor speed will reduce the displacement at the inner surface.The rotor speed should not exceed a certain speed which causes looseness between the shaft and the disc hub.The maximum total stress occurs at the middle, near the smallest thickness of the disc.Thickening this part can be proposed to reduce the maximum total stress as it was shown that stress concentration occurs there.

CONCLUSION
Analytical solution of rotating disc with variable thickness and nonhomogeneous material properties can not be obtained because there are many variable parameters in the coefficients of its governing differential equation.For this reason, numerous numerical approaches have been proposed to obtain the approximate solutions.One of the approaches having been widely used recently is FEM.
Most rotating discs used in applications are ax symmetric.Therefore, ax symmetric element is the most economical but adequate to use for the FE analysis.A turbine disc was analyzed as an example in this study.Several loads were applied.It was shown that maximum displacement occurs at the boundary of the disc, either at the outer boundary or the inner boundary, depending on the loadings.The maximum radial stress occurs at the area in the middle of the disc which has the smallest thickness.Furthermore, the analysis showed that each load gives different contribution to the total radial displacement and stress of the disc.The rotational blade load was shown to give the largest contribution.However, this does not apply to all turbine discs as it depends on the values of different loads for any specific turbine disc.
The radial displacement and stress in a disc with variable thickness are affected by the contour of the thickness variation.The middle area of the disc is proposed to have smaller thickness, but should not be too thin as stress concentration will occur there.
In the stress analysis, a quite fine meshing particularly in area with stress concentration is required to avoid inaccuracy of the stress calculation due to concentrated stress which cannot be spread well to surrounding area.

Fig. 3 :
Fig. 3: Stresses in the axisymmetric problem (21) Centrifugal forces due to rotational speed of the body forces.The element body force vector is (22)The element thermal forces working on the disc are(23)The total element forces {f} are sum of the element surface forces, the element body forces and the element thermal forces.The global force vectors {F} are obtained by assembling all the element forces.The element stiffness matrix is given by: (24) Assembling all of the element stiffness matrices, the global stiffness matrix [K] is obtained.The global displacement vectors are given by {F} = {K}{d}.The displacement function vectors are given by:

Fig. 4 :
Fig. 4: A half cross sectional view of the disc with dimensions in mm analysis.The thickness of the disc varies along the radial direction.The temperature varies direction as well.Due to the temperature variation, Young modulus, Poisson ratio and thermal coefficient change as function of the temperature.Because the temperature changes as function of radius, then Young modulus, Poisson ratio and the thermal coefficient also be expressed as function of radius.Because analytical solutions of such a problem is not available as already explained, finite element analysis was carried out to obtain the solutions.Geometry and material properties: axisymmetric.The geometry of the disc Fig. 4. For convenience, the change Young Modulus, Poisson ratio and thermal coefficient is assumed linear.The density of assumed to be constant.The steady state temperature at the outer boundary is equal to the temperature of high temperature gas.Based on the data given byBarack and Domas (1976),Liu et al. (2002),Claudio et al. (2002),Jahed (2005),Maruthi et al. (2012) andElhefny and (2012), this temperature has a range of 550 this study, a value of 800°C is used.Throughout the disc, the temperature gradually decreases as the radius gets closer to the inner boundary.The temperature throughout the disc never reaches cooling system is applied inside the disc.Based on data published byBarack and Domas (1976)  (2005)  and recently by Maruthi temperature at the inner boundary has a range of 450 500°C.In this study, a value of 500 Following the aforementioned boundary temperatures along with assumption that t changethroughout the disc is linear, the following relation is used to express the temperature change:

Fig. 9 :
Fig. 9: FEM model with the loads and boundary conditions Fig. 11: (a) Radial displacement (mm) and (b) stress (MPa) due to shrink fit only

Table 1 :
Maximum radial displacement and stress