Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 55, 1, pp. 331-341, Warsaw 2017 DOI: 10.15632/jtam-pl.55.1.331 NUMERICAL CREEP ANALYSIS OF FGM ROTATING DISC WITH GDQ METHOD Hodais Zharfi, Hamid Ekhteraei Toussi Ferdowsi University of Mashhad, Faculty of Engineering, Mashhad, I.R. Iran e-mail: ekhteraee@um.ac.ir Rotating discs are the vital part of many kinds of machineries. Usually, they are operating at relatively high angular velocity and temperature conditions. Accordingly, in practice, the creep analysis is an essential necessity in the study of rotating discs. In this paper, the ti- me dependent creep analysis of a thin Functionally Graded Material (FGM) rotating disc investigated using the GeneralizedDifferential Quadrature (GDQ)method. Creep is descri- bed with Sherby’s constitutivemodel. Secondary creep governing equations are derived and solved for a disc with two various boundary conditions and with linear distribution of SiC particles in pure Aluminum matrix. Since the creep rates are a function of stresses, time and temperature, there is not a closed form solution to these equations. Using a solution algorithm and the GDQ method, a solution procedure for these nonlinear equations is pre- sented. Comparison of the results with other existing creep studies in literature reveals the robustness, precision and high efficiency beside rapid convergence of the present approach. Keywords: creep, FGM,GDQ method, rotating disc, boundary conditions 1. Introduction A wide area research has assigned to rotating discs because of their numerous utilizations in various industries and rotating machines such as: gas and steam turbines, pumps, turbo gene- rators, compressors, flywheels, ship propellers and automotive braking systems (Gupta et al., 2005; Hojjati andHassani, 2008; Singh andRay, 2002). Inmost of these applications, discs have to operate at high rotational speed leading to large centrifugal forces, and in the presence of high temperatures, thematerials strength can be reduced explicitly. These operating conditions lead to a vast continuous time dependent deformations, so the creep phenomenon finds the prio- rity in the research. Creep deformation can affect the performance of systems entirely. As an example, in a turbine rotor there is always a possibility that the heat from the external surface is transmitted to the shaft and then to the bearings, which has adverse effects on the functio- ning and efficiency of the rotor (Bayat et al., 2008). Under high thermo-mechanical loading, the monolithic materials cannot do well. Augmenting the second strong and stable phase to the based phase causes a significant reduction in creep deformations. In other words, increasing the reinforcements leads to creep strength rising. For this reason and because of more other great benefits of FGMmaterials, thesematerials have attracted the interest ofmany researchers (Loghman et al., 2011; Ghorbani, 2012). In recent years, creep analysis in rotating discs made of functionally gradedmaterials has been addressed bymany researchers. Whal et al. (1954) studied steady state creep behavior of a turbine rotating disc using a power function creep law and Huber-Mises-Hencky and Tresca yield criteria. They compared results of their work with the experimental results. Arya and Bhatnagar (1979) investigated creep responses of an orthotropic rotating disc. They illustrated that tangential stress in all over the disc and the radial creep rate at inner radius increased with intensification of mate- rial anisotropy. Nieh (1984) demonstrated that a disc with SiC whiskers in Aluminum matrix 332 H. Zharfi, H.E. Toussi has a better creep responses than a disc made of pure Aluminum. Białkiewicz (1986) presen- ted a theoretical finite strain analysis in rotating discs based on creep rupture using Kachanov theory and Norton’s creep law. Bhatnagar et al. (1986) studied steady state creep behavior in three different discs with anisotropic materials. The discs had constant, linear and hyperbolic decreasing thicknesses. The results illustrated that discs with decreasing thickness in the radial direction had better creep responses. Pandey et al. (1992) studied steady state creep behavior of an Al-SiC composite rotating disc in a uniaxial load state and temperature between 623-723K. They investigated various compositions of particle size (45.9µm, 14.5µm, 1.7µm)with the total volume content of 10%, 20%, and 30%. They observed that the disc with finer particles shows the better creep responses. Gupta et al. (2004) investigated the secondary creep response in a rotating disc made of pure Aluminum with SiC particle distribution in the radial direction. The creep law was Sherby’s model. Their results showed that the tangential and radial stres- ses were not affect by governing temperature and particle distributions. Radial and tangential creep rate reduced with particle size reduction, additional of particle volume content and de- creased governing temperature. Singh and Ray (2001) studied secondary creep behavior in a thin anisotropic rotating disc made of Aluminum with SiC whickers. The effects of anisotropy and non-homogenous distribution of reinforcements on creep behavior of the rotating disc were investigated. TheGeneralized Differential Quadraturemethodwas proposed byBellman et al. (1972) and extended by Shu and coworkers (Shu and Richards, 1992). In this numerical method, function derivatives are approximated in terms of linear summation andweighting coefficients in various grid points. Computation of the weighting coefficient is done using high order polynomial simu- lation and linear vector space analysis. The weighting coefficient for the first order derivatives is computed by a simple algebraic formulation and higher order coefficients are obtained using a recursive relationship based on the first order weighting coefficient. The major advantage of the GDQmethod over the DQmethodwhich was extended by Bellman et al. (1972), is its ease in computing the weighting coefficients without any restriction to the grid point selection. The GDQ method has been successfully applied to some fluid flow problems (Shu et al., 1995) and structural vibration analysis (Shu, 1996). In all cases, the method valuable characteristic such as high precision and fast convergence was confirmed. Fig. 1. Studied disc and its parameter In this paper, the generalized differential quadrature method is used to obtain numerical solution of creep governing equations. For this purpose, the thermo-elastic creep equations of an FGMrotating disc shown inFig. 1, under free-free andfixed-free boundary conditions is derived. Using a solution algorithm and the GDQ method, creep response of the considered structure is achieved. There is no modeling involved, and the convergence to the final response is very fast. Comparison of the results also illustrates that themethod is very efficient and accurate. It avoids simplifications and restrictions which other creep solution methods in the literature are involved, therefore it is a great merit for creep studies of structures. Numerical creep analysis of FGM rotating disc with GDQ method 333 2. Mathematical formulation Before starting to develop a closed form solution for the distribution of creep stress and defor- mation, there are some assumptions that must be asserted. Amongst them presumptions are made that the material of the disc is isotropic and its yield criterion complies with the Huber- -Mises-Henckymodel. In addition, it is supposed that compared to other dimensions of the disc its thickness is too small so that axial stresses can be ignored and plain stress conditions can be assumed. The total strain in an FGM rotating disc is an ensemble of elastic, thermal and creep strains in the following form εtotal = εelastic +εthermal +εcreep (2.1) For derivation of the primary and secondary creep equation, at firstwe use a combination of the strain-displacement relation and Eq. (2.1) εr = du dr = 1 E (σr−νσθ)+αT +εr,c εθ = u r = 1 E (σθ−νσr)+αT +εθ,c (2.2) In the above two equations, εr and εθ are the radial and tangential strains. σr and σθ are the radial and tangential stresses. u is the radial displacement and E, α and T are the modulus of elasticity, thermal expansion coefficient and governing temperature in Kelvin, respectively. εr,c and εθ,c are the radial and tangential creep strains. The stress relations can be concluded from Eqs. (2.2) as below σr = E 1−ν2 [du dr +ν u r − (1+ν)αT − (εr,c+νεθ,c) ] σθ = E 1−ν2 [u r +ν du dr − (1+ν)αT − (νεr,c+εθ,c) ] (2.3) Using the equilibrium condition for an element of the rotating disc in form of Fig. 2, we obtain the equilibrium equation as dσr dr + σr−σθ r +ρrω2 =0 (2.4) With replacing the radial and tangential stresses fromEqs. (2.3) into Eq. (2.4), the simplest form of the creep equation is obtained as below A d2u dr2 +B du dr +Cu+D=0 (2.5) Inwhich the coefficientsA,B,C andD inEq. (2.5) presented in relations given below. It should be noted that all FGM disc parameters are a function of the radius A= E 1−ν2 B= E r(1+ν) + 1 1−ν2 dE dr + 2νE (1−ν2)2 dν dr + νE r(1−ν2) C = −E r2(1+ν) + ν r(1−ν2) dE dr + 2Eν2 r(1−ν2)2 dν dr + E r(1−ν2) dν dr − νE r2(1−ν2) D= E r(1−ν2) [(ν−1)εr,c+(1−ν)εθ,c] + 1 1−ν2 [ (1−ν2) dE dr +2νE dν dr ] [−(1+ν)αT −εr,c−νεθ,c] + E 1−ν2 [ − dν dr αT − (1+ν) dα dr T − (1+ν)α dT dr − dεr,c dr − dν dr εθ,c−ν dεθ,c dr ] (2.6) 334 H. Zharfi, H.E. Toussi Fig. 2. Schematic of an element of the rotating disc Creep governing equation for the rotating disc in terms of the displacement rate can be repre- sented as A d2u̇ dr2 +B u̇ dr +Cu̇=−Ḋ (2.7) where u̇ is the radial displacement rate, A, B and C are presented in Eqs. (2.6)1−3 and Ḋ is a function of the creep rates as below Ḋ= E r(1−ν2) [(ν−1)ε̇r,c+(1−ν)ε̇θ,c]+ 1 1−ν2 [ (1−ν2) dE dr +2νE dν dr ] (−ε̇r,c−νε̇θ,c) + E 1−ν2 [ − ε̇r,c dr − dν dr ε̇θ,c−ν ε̇θ,c dr ] (2.8) where ε̇r,c and ε̇θ,c are the radial and tangential creep rate components, respectively, which can be computed using Sherby’s creep law as ε̇r,c = [M(σe−σ0)] n 2σe (2σr −σθ) ε̇θ,c = [M(σe−σ0)] n 2σe (2σθ −σr) (2.9) where in the above two equations, n is the creep exponent and, in this study, is considered to be 8. σe is the effective stress which is obtained from the Huber-Mises-Hencky isotropic yield criterion for the plane stress problem as Eq. (2.10).M and σ0 are creep parameters and depend on the particle size P which is considered to be 1.7µm, particle distribution function v(r), and the level of prevailing temperatureT . Based on the experimental creep data reported byPandey et al. (1992) and using a regression technique, Gupta et al. (2004) representedM and σ0 as the following functions Eqs. (2.11) σe = √ σ2r +σ 2 θ −σrσθ (2.10) and lnM =−35.38+0.2077lnP +4.98lnT −0.622lnv(r) σ0 =−2.11916−0.03507P +0.01057T +1.00536v(r) (2.11) As it is mentioned above, no exact solution to Eq. (2.7) can be obtained because its right hand side contains the functions of creep rates which all are time, stress and temperature dependent. However, there are numerous numerical methods for solving this problem. In this study, the GeneralizedDifferential Quadraturemethod is usedwith a solution techniquewhich is explained in the next section. Numerical creep analysis of FGM rotating disc with GDQ method 335 3. Solution algorithm The procedure of creep analysis in anFGM rotating disc is initiatedwith aGDQ thermo-elastic solution of Eq. (2.5) in which the creep strains are ignored. Afterwards, using this displacement field, the radial and tangential stresses are calculated using Eqs. (2.3). Within this stress field, thedistributionof radial and tangential strain rates is obtainedusingEqs. (2.9). ThenwithGDQ analysis of Eq. (2.7), the displacement rate can be computed so the distribution of the radial and tangential stress rates are obtained. Afterwards, by selecting a suitable time interval ∆t, the next radial and tangential stress at the next time step are calculated using equations σt+1 = σ̇t∆t+σt εt+1 = ε̇t∆t+εt (3.1) This procedure is continued until the stress distributions converge and reach the steady state condition. Asmentioned above, for the numerical creep analysis of the FGM rotating disc, a numerical procedure is inevitable. In this study the Generalized differential Quadrature method is used because of high precision and fast convergence of this method. The Differential Quadrature idea was proposed by Bellman et al. (1972), and extended by many other researchers later on. In this method, the partial differential of a function with respect to a coordinate direction is expressed as a linear weighted sumof all the functional values at all grid points in that direction (Tornabene et al., 2016). For a smooth function f(r), GDQ discretizes its n-th order derivative with respect to r at the grid point (ri) as f(n)r (ri)= N ∑ k=1 a (n) ik f(rk) n=1,2, . . . ,N−1 i=1,2, . . . ,N (3.2) where N is the number of grid points in the r direction. a (n) ik is the weighting coefficient for the second and higher order derivative which can be completely determined from the first order derivatives as a (n) ij =            n ( a (n−1) ii a (1) ij − a (n−1) ij ri−rj j 6= i N ∑ k=1,k 6=i a (n) ik j= i for i,j=1,2, . . . ,N, n=2,3, . . . ,N −1 (3.3) In the above equation a (1) ij is the first weighting coefficient to be obtained using a (1) ij =          A(1)(ri) (ri−rj)A(1)(rj) j 6= i N ∑ k=1,k 6=i a (1) ik j= i for i,j=1,2, . . . ,N (3.4) where A(1)(ri)= N ∏ j=1,j 6=i (ri−rj) (3.5) When the coordinates of grid points are known, the weighting coefficient and so the discretized derivative canbeeasily calculated fromEqs. (3.2)-(3.5) (ShuandChew,1999).Various gridpoint 336 H. Zharfi, H.E. Toussi distributions are investigated in the literature. In this study roots of the Chebyshev polynomial of the first kind are used for grid point generation as (Tornabene et al., 2012) ri = gi−g1 gN −g1 gi =cos (2i−1 2N π ) i=1,2, . . . ,N (3.6) whereN is the total number of grid points along the radial direction. For the problem studied herein, the GDQmethodwhich is presented by the above equations demonstrates its numerical accuracy, high speed and computation amount reduction. 4. Numerical application and discussion In this Section, the above formulation for creep analysis is applied to a FGM rotating disc with two different boundary conditions. For this purpose, a FGMdisc, as can be seen in Fig. 1, with inner radius a= 0.05m and outer radius b= 0.2m is considered. It is subjected to T = 623K temperature and rotates withω=15000rpm.Thematerial of the disc is silicon carbide particles distributed with a linear volume fraction in pure Aluminummatrix v(r)= vmax− (vmax−vmin) r−a b−a (4.1) In the above equation vmax = 0.4 and vmin = 0.1 are the maximum and minimum volume fraction of silicon carbide at the inner and outer radius, respectively. The boundary conditions of the disc are are established as: — free-free boundary condition σr(r= a)= 0 σr(r= b)= 0 (4.2) — fixed-free boundary condition ur(r= a)= 0 σr(r= b)= 0 (4.3) In the distribution patterns of the non-homogeneous disc, thematerial heterogeneity ismild and the change in thermo-mechanical properties is not very high. Therefore, similarly to the work done by Loghman et al. (2011), the disc parameter such as density, Poisson’s ratio, thermal expansion coefficient and elasticity modulus can be calculated using a linearmixture rule which is presentedbyEq. (4.4). In otherwords, based onEq. (4.4), thematerial and thermal properties depend on the reinforcement particles volume fraction at each point P(r)=PAl+(PSiC −PAl)v(r) (4.4) In the above equation,P(r) is the intended property in the composite disc,PAl andPSiC are the value of the same property in pure Aluminum and silicon carbide, respectively. The following data for pure material properties are used in this paper: EAl = 69GPa, ESiC = 410GPa, ρAl =2700Kgm −3, ρSiC =3200kgm −3, αAl =23.1 ·10 −6K−1, αSiC =4 ·10 −6K−1. Based on the creep analysis algorithm and the above specified parameters, a computer code has been developed to find the creep response of the FGM rotating disc. Before providing the results of the FGM rotating disc, to validate the analysis and the developed computer code, the results of tangential strain rate have been computed for a disc which was used by Loghman et al. (2011) and Ghorbani (2012). A comparison between the results is shown in Fig. 3. It can be seen that there is a good agreement between the results. However, the time consumed and computation amount by the other method is much higher. Numerical creep analysis of FGM rotating disc with GDQ method 337 Fig. 3. Results validation The radial stress distribution versus the radial direction is presented in Fig. 4. Through the free-free boundary condition, the value of radial stress in the inner and the outer radius of the disc becomes zero and the maximum value of radial stress occurs near the middle of the disc. The disc with the fixed-free boundary condition behaves in a different manner. As it can be seen, for this condition the radial stress has its maximum value at the inner radius and then a decreasing trend along the radial direction appears, so theminimum values are observed at the outer radius. Fig. 4. Radial stress distribution for two various boundary conditions Figure 5 confirms that the tangential stress for both boundary conditions show a decreasing trend in the radial direction. It means that tangential stress experiences a maximum near the inner radii and gradually decreases to reach a minimum at the outer rim. It seems that the reason for this type of behavior is the existence of a relatively higher content of particles at the inner radius in comparisonwith other points in the disc. But it is obviously seen that the values of tangential stress in the disc with fixed-free boundary conditions is lower than in the other disc. Thedistributionof the radial componentof creep strain rate along thedisc radius is presented in Fig. 6. As it can be seen, for both boundary conditions the radial creep rate has itsmaximum value at the inner rim and a decreasing trend along the radial direction. But their behaviour is different, the fixed-free disc has a positive radial creep rate and the free-free disc has a 338 H. Zharfi, H.E. Toussi Fig. 5. Tangential stress distribution for two various boundary conditions Fig. 6. Radial strain rate distribution for two various boundary conditions negative one, as well as the absolute value of the radial creep rate in the disc with the fixed-free boundary condition is higher than the disc with the free-free boundary condition. Therefore, it is important to choose appropriate boundary conditions during the disc installation because of high dependency of the FGMdisc behavior on boundary conditions. As depicted, the strain rate values in fixed-free boundary conditions are much higher than in the other ones, however this situation is the most practical condition. Figure 7 shows the disc with the fixed-free boundary condition. It has an increasing- -decreasing trend via the radial direction, but in the disc with the free-free boundary condi- tion, the tangential creep rate is maximal at the inner side of the disc and minimal at the outer side. It has amonotonic decreasing behavior at the intermediate radii. It is observed that similarly to the radial creep rate, the maximum rate of the tangential creep rate takes place near the inner radius of discs, also themaximum value of the effective stress occurs at the inner rims, so it is necessary to reinforce the inner side of disc with high strength materials. The- refore, the FGM materials with the higher values of particle reinforcements at the inner radii, in other words, discs with the decreasing type of particle distribution along the radial direc- tion have a better creep responses in comparison with the increasing type ones or homogenous materials. Numerical creep analysis of FGM rotating disc with GDQ method 339 Fig. 7. Tangential strain rate distribution for two various boundary conditions Fig. 8. Radial displacement rate distribution for two various boundary conditions The variation of the radial displacement rate for two boundary conditions is presented in Fig. 8.As it canbeseen, thediscwith the free-freeboundaryconditionhasadecreasing trendand thediscwith thefixed-free boundaryconditionhas an increasing trendalong the radial direction, as expected. Obviously, becouse of the integrated nature of the displacement (compared to the differential nature of strain) this picture does not reveal many details appearing in the previous stress or strain curves. 5. Conclusion In thispaper, thegeneralizeddifferential quadraturemethod isused toobtain the timedependent creep responses of anFGMrotating disc. The creep equations are derived for a discwith a linear distribution of the particle content and are solved for two various boundary conditions. It can be demonstrated that the GDQ method is an efficient and accurate method in creep analysis of rotating discs. As it is shown, comparisons of the results with other available approaches show a good agreement. However, the time consumed by the other methods for solving the same problem is much higher. The results also illustrate that the material inhomogeneity has a considerable effect on the creep behavior of FGM rotating discs. Creep rates and stress fields at the inner side of the disc is much greater, so the disc must to be reinforced at this internal 340 H. Zharfi, H.E. Toussi side. Making use of the decreasing type of particle distribution leads to better creep responses as well as higher creep life time. References 1. Arya V.K., Bhatnagar N.S., 1979, Creep analysis of rotating orthotropic discs, International Journal of Nuclear Engineering and Design, 55, 323-330 2. Bayat M., Sleem M., Sahari B.B., Hamouda A.M.S., Mahdi E., 2008, Analysis of func- tionally graded rotating disks with variable thickness, Mechanics Research Communications, 35, 283-309 3. Bellman R., Kashef B.G., Casti J., 1972, Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations, Journal of Computational Physics, 10, 40-52 4. Bhatnagar N.S., Kulkarni P.S., Arya V.K., 1986, Steady state creep of orthotropic rotating discs of variable thickness, International Journal of Nuclear Engineering and Design, 91, 121-141 5. Białkiewicz J., 1986,Dynamic creep rupture of a rotatingdisc of variable thickness, International Journal of Mechanical Science, 28, 671-681 6. GhorbaniM.T., 2012,A semi-analytical solution for time-variant thermoelastic creep analysis of functionally graded rotating disks with variable thickness and properties, International Journal of Advanced Design and Manufacturing Technology, 5, 2, 41-50 7. Gupta V.K., Singh S.B., Chandrawat H.N., Ray S., 2004, Steady state creep and material parameters in a rotating disc of Al-SiC composites,European Journal of Mechanics A Solids, 23, 335-344 8. Gupta V.K., Singh S.B., Chandrawat H.N., Ray S., 2005, Modeling of creep behavior of a rotating disc in the presence of both composition and thermal gradients, Journal of Engineering Materials and Technology, 127, 97-105 9. Hojjati M.H., Hassani A., 2008, Theoretical and numerical analysis of rotating discs of non- uniformthicknessanddensity, International Journal ofPressureVessels andPiping,85, 10, 694-700 10. Loghman A., Gorbanpour Arani A., Shajari A.R., Amir S., 2011, Time dependent ther- moelastic creep analysis of rotating diskmade of Al-SiC composite,Archive of Applied Mechanics, 81, 1853-1864 11. Nieh T.G., 1984, Creep rupture of a silicon carbide reinforced aluminum composite,Metallurgical Transactions, 15A, 139-146 12. PandeyA.B.,MishraR.S.,MahajanY.R., 1992,Steady state creepbehavior of silicon carbide particulate reinforced aluminum composites,Acta Metallurgica et Materialia, 40, 2045-2082 13. Shu C., 1996, Free vibration analysis of composite laminated conical shells by generalized diffe- rential quadrature, Journal of Sound and Vibration, 194, 587-604 14. Shu C., Chew Y.T., 1999, Application of multi-domain GDQ method to analysis of waveguides with rectangular boundaries,Progress in Electromagnetics Research, 21, 1-19 15. Shu C., Chew Y.T., Richards B.E., 1995, Generalized differential integral quadrature and their application to solve boundary layer equations, International Journal of Numerical Methods in Fluids, 21, 723-733 16. Shu C., Richards B.E., 1992, Application of generalized differential quadrature to solve two- dimensional incompressible Navier-Stokes equations, International Journal of Numerical Methods in Fluids, 15, 791-798 17. SinghS.B.,Ray S., 2001,Steady state creepbehavior in an isotropic functionally gradedmaterial rotating disc of Al-SiC composite,Metallurgical and Materials Transactions, 32A, 1679-1685 18. Singh S.B., Ray S., 2002, Modeling the anisotropy and creep in orthotropic aluminum-silicon carbide composite rotating disc,Mechanics of Materials, 34, 363-372 Numerical creep analysis of FGM rotating disc with GDQ method 341 19. Tornabene F., Fantuzzi N., Bacciocchi M., 2016, The GDQ method for the free vibration analysis of arbitrarily shaped laminated composites shells using a NURBS-based isogeometric ap- proach,Composite Structures, 153, 190-218 20. Tornabene F., Liverani A., Caligiana G., 2012, Laminated composite rectangular and an- nular plates: AGDQ solution for static analysis with a posteriori shear and normal stress recovery, Composite: Part B, 43, 1847-1872 21. WhalA.M., SankeyG.O.,ManjoineM.J., ShoemakerE., 1954,Creep test of rotating discs at elevated temperature and comparisons with theory, Journal of Applied Mechanics, 21, 225-235 Manuscript received May 28, 2016; accepted for print September 12, 2016