Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 54, 1, pp. 113-125, Warsaw 2016 DOI: 10.15632/jtam-pl.54.1.113 ELASTIC-PLASTIC ANALYSIS OF PRESSURE VESSELS AND ROTATING DISKS MADE OF FUNCTIONALLY GRADED MATERIALS USING THE ISOGEOMETRIC APPROACH Amir T. Kalali, Behrooz Hassani, Saied Hadidi-Moud Mechanical Engineering Department, Ferdowsi University of Mashhad, Iran e-mail: amir.kalali.121@stu-mail.um.ac.ir; b hassani@um.ac.ir; hadidi@um.ac.ir An NURBS-based isogeometric analysis for elastic-plastic stress in a cylindrical pressure vessel is presented. The vessel is made of a ceramic/metal functionally gradedmaterial, i.e. a particle-reinforced composite. It is assumed that the material plastic deformation follows an isotropic strain-hardening rule based on the von Mises yield criterion. The mechanical properties of the graded material are modelled by the modified rule of mixtures. Selected finite element results are also presented to establish the supporting evidence for validation of the isogeometric analysis. Similar analyses areperformedand solutions for spherical pressure vessel and rotating disk made of FGMs are also provided. Keywords: isogeometric analysis, NURBS, functionally graded material, modified rule of mixtures 1. Introduction The intensity and variation of stress distributions due to large mismatch in material proper- ties can be substantially reduced if micro-structural transition behaviour, i.e. a gradedmaterial model, is used. Advances in material synthesis technologies have spurred the development of functionally gradedmaterials (FGM) with promising applications in aerospace, transportation, energy, cutting tools, electronics, and biomedical engineering (Chakraborty et al., 2003). An FGMcomprises amulti-phasematerial with volume fractions of the constituents varying gradu- ally in a predetermined profile, thus yielding a non-uniformmicrostructure in thematerial with continuously graded properties (Jin et al., 2003). Elastic and elastic-plastic analyses of thick-walled pressure vessels have always attracted a lot of research interest because of their importance in engineering applications. Figueiredo et al. (2008) proposed a numerical methodology in order to predict the elastic-plastic stress beha- viour of functionally graded cylindrical vessels subjected to internal pressure. It was assumed that the structures undergo small strain and that the material properties of the graded layer were modelled by the modified rule of mixtures approximation. Furthermore, the plastic do- main for ductile phases was defined through the von Mises yield criterion. They proposed an iterativemethod for solving the nonlinear system combining a finite element approximation and an incremental-iterative scheme. Haghpanah Jahromi et al. (2009, 2010) extended the Variable Material Property (VMP) method developed by Jahed and Dubey (1997) for materials with varying elastic and plastic properties. In the VMP method, the linear elastic solution to the boundaryvalue problemwas used as a basis to generate the inelastic solution. Through iterative analysis, the VMP method was used to obtain the distribution of material parameters which were considered as field variables. The application of the VMP method, generally applied to homogeneous elastic-plastic materials (Jahed and Shirazi, 2001; Jahed et al., 2005, 2006), was extended to materials with varying elastic-plastic properties in order to calculate the residual stresses in an autofrettaged FGM cylindrical vessel. 114 A.T. Kalali et al. Although there are several papers on the elastic analysis of FGMspherical pressurevessels in the literature (You et al., 2005;Dai et al., 2006;ChenandLin,2008), elastic-plastic stressanalysis ofFGMspherical pressurevessels is not sucha customary study. SadeghianandEkhteraei (2011) studied thermal stress field for an FGM spherical pressure vessel made of an elastic-perfectly plastic and a power law material model. Similar to the FGM cylindrical and spherical vessels, much of the studies on FGM rotating disks has been carried out in elasticity cases (Durodola and Attia, 2000; Bayat et al., 2008). HaghpanahJahromi et al. (2012) applied theVMPmethod to estimate the elasto-plastic stresses in a rotating disk with varying elastic and plastic properties in the radial direction. In this paper, isogeometric analysis is proposed for predicting stress components of a strain- -hardening cylinder based on the von Mises yield criterion under plane stress conditions. Isoge- ometric analysis was introduced byHughes et al. (2005) as a generalisation of the standard finite element analysis. In isogeometric analysis, the solution space for dependent variables is represen- ted in terms of the same functionswhich represent the geometry.Thegeometric representation is typically smooth, whereas the solution space for the standard finite element analysis is continu- ous but not smooth.Adoption of the isogeometric concept has shown computational advantages over the standard finite element analysis in terms of accuracy and analysis time in many ap- plication areas, including solid and structural mechanics. Most CAD systems use spline basis functions and oftenNon-UniformRational B-Splines (NURBS) of different polynomial orders to represent geometry. Results obtained fromfinite element analysis using the commercial software ABAQUS (v. 6.10)were used to validate the results from the isogeometric analysis. The analysis was further extended to obtain solutions for FGM spherical vessels and rotating disks. A brief review of the isogeometric analysis based on NURBS is presented in Section 2. This is followed in Section 3 by describing the details of isogeometric analysis formulation for elastic-plastic cases (functionally graded cylindrical and spherical vessels and rotating disks). In Section 4, we describe material properties of the graded layer modelled by the modified rule of mixtures, whereas in Section 5 the results of elastic-plastic analyses are presented. Finally, in Section 6, key conclusions are pointed out. 2. Fundamentals of NURBS-based isogeometric analysis 2.1. B-splines and NURBS Non-uniform rational B-splines (NURBS) are a standard tool for describing and modelling curves and surfaces in the computer aided design and computer graphics. B-splines are piecewise polynomial curves composed of linear combinations of B-spline basis functions. The piecewise definition allows approximation of a large number of control points using lower order polynomials. The coefficients are points in space, referred to as the control points. A knot vector Ξ is a set of non-decreasing real numbers representing coordinates in the parametric space of the curve Ξ= [ξ1,ξ2,ξ3, . . . ,ξi, . . . ,ξn+p+1] (2.1) where p is the order of the B-spline and n refers to the number of the basis functions (also control points). The interval [ξ1,ξn+p+1] is called a patch. The B-spline basis functions for a given degree p are defined recursively over the parametric domain by the knot vector. The piecewise constants are first defined as Ni,0(ξ)= { 1 if ξi ¬ ξ < ξi+1 0 otherwise (2.2) Elastic-plastic analysis of pressure vessels and rotating disks... 115 For p> 0, the basis functions are defined by the following recursion formula Ni,p(ξ)= ξ− ξi ξi+p− ξi Ni,p−1(ξ)+ ξi+p+1− ξ ξi+p+1− ξi+1 Ni+1,p−1(ξ) (2.3) A B-spline surface is constructed by the basis functions in two directions, Ni,p(ξ) and Mj,q(η), and a set of control points Pij, i = 1,2, . . . ,n, j = 1,2, . . . ,m. Similar to the first para- metric direction ξ, Mj,q(η) is also defined by Eqs. (2.2) and (2.3), but another knot vector H= [η1,η2,η3, . . . ,ηj, . . . ,ηm+q+1] constitutes the foundation. Often, the B-spline order is the same in both directions, i.e. p= q. The surface is to be drawn in the two-dimensional space, Pij ∈R 2. The B-spline surface is then interpolated by S(ξ,η) = n ∑ i=1 m ∑ j=1 Ni,p(ξ)Mj,q(η)Pij (2.4) The B-spline surface is the result of a tensor product. The patch for the surface is now the domain [ξ1,ξn+p+1]×[η1,ηn+p+1]. Identifying the logical coordinates (i,j) of theB-spline surface with the traditional notation of the node A and the Cartesian product of the associated basis functionswith the shape functionNA(ξ,η) =Ni,p(ξ)Mj,q(η), the familiar finite element notation is recovered, namely S(ξ,η) = nm ∑ A=1 NA(ξ,η)PA (2.5) B-splines are non-rational functions that form non-rational B-spline curves and surfaces. A rational curve or surface can represent conical sections in an exactmanner.Non-uniformrational B-splines (NURBS) are therefore introduced by including weights on the control points. The NURBS basis functions will differ from the B-spline basis functions, but the knot vectors, the tensor product nature, and refinementmechanisms are unchanged. The NURBS surface is given by S(ξ,η) = 1 w(ξ,η) nm ∑ A=1 NA(ξ,η)wAPA = nm ∑ A=1 NA(ξ,η)PA (2.6) where w(ξ,η) = nm ∑ A=1 NA(ξ,η)wA NA(ξ,η) = NA(ξ,η)wA w(ξ,η) 2.2. Fundamentals of the isogeometric analysis The isogeometric analysis was defined by Hughes et al. (2005) and means that the analysis model uses the samemathematical description as the geometry model. This notion of using the same basis for geometry and analysis is called the isoparametric concept, and it is quite common in the classical finite element analysis.The fundamental differencebetween the isogeometric ana- lysis and the finite element analysis is that, in the FEA, the basis chosen for the approximation of the unknown solution fields is used to approximate known geometrywhereas the isogeometric analysis turns this idea around and selects a basis capable of exactly representing the known geometry, and uses it as a basis for the fields we wish to approximate (Cottrell et al., 2009). The main advantages of the isogeometrical method, compared to other numerical methods, can be summarised as below: • feduction in size of the system of equations, 116 A.T. Kalali et al. • flexibility and accuracy in the definition of geometry and its boundaries, • thepossibility of keeping the originalmodel in thewhole processwithout several remeshing in problems with a varying domain of interest, • considerable ease in implementing adaptively andmesh refinement, • accuracy in satisfaction of the essential boundary conditions, • applicability of the method in problems of functionally graded materials (Hassani et al., 2011). And themain disadvantages of the isogeometric analysis can bementioned as in the following: • the control points of geometry commonly are not a part of the physical domain of the problem, • the relative difficulty of establishing a correspondence between the point in the domain and the solution. 3. Isogeometric analysis formulation for the elastic-plastic case 3.1. Elastic formulation We use the principle of virtual displacement applied to a plane elastic body 0= ∫ V (σijδεij +ρüiδui) dV − ∫ V fiδui dV − ∮ S tiδui ds (3.1) where σij is the Cauchy stress, εij is strain, ρ is density, üi is the acceleration component, ui is the displacement component (i.e. u,v,w), fi is the body force component, ti is the traction component (i.e. tx, ty, tz), V is volume and S is surface area corresponding to the volume. In the cylindrical coordinate system and axisymmetric condition, Eq. (3.1) can be rewritten as follows 0=2π ∫∫ (σijδεij +ρüiδui)r drdz−2π ∫∫ fiδuir dr− ∮ S tiδui ds (3.2) By using the NURBS basis functions, the approximated displacement functions can be written as r(ξ,η) = nm ∑ A=1 NA(ξ,η)rA = [N1(ξ,η),N2(ξ,η), . . . ,Nnm(ξ,η)]       r1 r2 ... rnm       =Nr z(ξ,η) = nm ∑ A=1 NA(ξ,η)zA =Nz u(ξ,η) = nm ∑ A=1 NA(ξ,η)uA =Nu w(ξ,η) = nm ∑ A=1 NA(ξ,η)wA =Nw δu(ξ,η) = nm ∑ A=1 NA(ξ,η)δuA =Nδu δw(ξ,η) = nm ∑ A=1 NA(ξ,η)δwA =Nδw (3.3) where rA and zAare the x- and y-coordinates of the control points of the surface, uA and wA are the control points of the displacement. Elastic-plastic analysis of pressure vessels and rotating disks... 117 The stress and strain relationships are given by ε=TN [ u w ] =B [ u w ] σ=Cε=CB [ u w ] (3.4) whereT is thematrix of differential operators,C is the constitutivematrix (constitutivematrix is calculated viaYoung’smodulus andPoisson’s ratio) andB=TN. In this study it is assumed thatPoisson’s ratio ν is amaterial constantwhileYoung’smodulusE(r) varieswith the position across the wall thickness of the vessel (disk). Substituting Eqs. (3.3) and (3.4) into Eq. (3.2), and in the absence of inertia forces, we obtain 0=2π ∫∫ ( B T CB [ u w ]) r drdz−2π ∫∫ N T [ fr fz ] r drdz− ∮ S N T [ tr tz ] ds (3.5) Note that in Eq. (3.5) all variables are written in terms of the parameters ξ and η which is similar to mapping in the standard finite element method where the base or unit elements are used 0=2π ξn+p+1 ∫ ξ1 ηm+q+1 ∫ η1 B T CB [ u w ] r(detJ) dηdξ −2π ξn+p+1 ∫ ξ1 ηm+q+1 ∫ η1 N T [ fr fz ] r(detJ) dηdξ− ∮ S N T [ tr tz ] ds (3.6) where J= [ ∂r/∂ξ ∂z/∂ξ ∂r/∂η ∂z/∂η ] and thematrix form is as follows Ku=F+T (3.7) where K=2π ξn+p+1 ∫ ξ1 ηm+q+1 ∫ η1 B T CB [ u w ] r(detJ) dηdξ F=2π ξn+p+1 ∫ ξ1 ηm+q+1 ∫ η1 N T [ fr fz ] r(detJ) dηdξ T= ∮ S N T [ tr tz ] ds (3.8) Integrals in Eq. (3.8) can be calculated using the Gauss-Legendre method of numerical integra- tion. Inorder toobtain stressdistributions for a spherical thick-walled functionally gradedpressure vessel, the isogeometric analysis formulation is rewritten in the spherical coordinate system (r,θ,ϕ). 118 A.T. Kalali et al. 3.2. Plastic formulation In order to formulate a theorywhichmodels elasto-plastic material deformation, three requ- irements have to bemet: • before the onset of plastic deformation, an explicit relationship between stress and strain must be formulated to describe material behaviour under elastic conditions, • a yield criterion must be postulated to indicate the stress level at which plastic flow commences, • a relationship between stress and strainmust be developed for post yield behaviour, when the deformation is made up of both elastic and plastic components (Owen and Hinton, 1980). According to thenormality hypothesis of plasticity, theplastic strain increment dεp is defined as: dεp = dλ ∂Q ∂σ (3.9) whereQ is the yield function and dλ is called the plastic multiplier. Assuming that the material plastic deformation follows the isotropic strain-hardening rule based on the vonMises yield criterionn≡ ∂Q/∂σ= [3/(2σe)]S and dλ= dε p e, Eq. (3.9) may be rewritten as dε p = 3 2 dεpe σe S (3.10) where dεpe is the equivalent plastic strain increment. The superscripts p and e denote plasticity and elasticity conditions respectively, also the subscript e denotes equivalent (effective) parame- ters (stress or plastic strain). The equivalent stress σe and the deviatoric stress S for the plane stress field are defined as σe = √ σ2r +σ 2 θ −σrσθ S=    Sr Sθ Sz    =        2σr−σθ 3 2σθ−σr 3 − σr+σθ 3        (3.11) where σr and σθ are the radial and hoop stresses, respectively. For a linear strain hardening material (Fig. 1), the yield stress σy and the plastic multiplier dλ are determined by σy =σy0(r)+hp(r)ε p e dλ= nTCdε nTCn+hp (3.12) where hp(r) is the plasticity modulus (i.e. the gradient of the stress-plastic strain curve) and σy0(r) is the initial yield stress of the FGM material. Both hp(r) and σy0(r) are functions dependent on the radial position r. The stress increment is given by dσ=Cdεe =C(dε−dεp)=C(dε−dλn) (3.13) by substitutingEq. (3.12)2 intoEq. (3.13), we can obtain the complete elasto-plastic incremental stress-strain relation dσ=Cepdε Cep =C− CnnTC nTCn+hp (3.14) where the superscripts ep denote the elasto-plastic behaviour. Elastic-plastic analysis of pressure vessels and rotating disks... 119 Fig. 1. Stress-strain curve for linear strain hardening If we denote all quantities at the iteration k with a superscript k, and those at the next increment by k+1 in a similar way, then wemay write dλk = nk T Cdεk nk T Cnk+hp dσ k =C ( dε k−dλknk ) (3.15) The integration to obtain the quantity at the end of the time step∆tmay be then written as σ k+1 =σk+dσk (3.16) If relatively large load increments are to be permitted, the process described can lead to an inaccurate prediction of the stresses. Two parameters R (reduction factor) and m (the excess of the yield stress is divided into m parts) can help to minimize the error (Owen and Hinton, 1980). The algorithm for plasticity isogeometric analysis is summarized in Fig. 2. Fig. 2. Flow chart of the algorithm for the elastic-plastic analysis 120 A.T. Kalali et al. 4. Mechanical behaviour of FGM It is assumed that the functionally gradedmetal-ceramic composite is locally isotropic and fol- lows the von Mises yield criterion. The three important material properties for elastic-plastic analysis are the elastic modulus E(r), the initial yield stress σy0(r), and the tangent modu- lusH(r). These properties can be calculated using themodified rule of mixtures for composites (Suresh and Mortensen, 1998). Note that the modified rule of mixtures is appropriate for mo- deling of isotropic materials E = [ (1−fc)Em q+Ec q+Em +fcEc ][ (1−fc) q+Ec q+Em +fc ] −1 σy0 =σy0m [ (1−fc)+ q+Em q+Ec Ec Em fc ] H = [ (1−fc)Hm q+Ec q+Hm +fcEc ][ (1−fc) q+Ec q+Hm +fc ] −1 hp = EH E−H (4.1) where the subscripts c andm indicate the ceramic andmetal material, respectively. The volume fraction of ceramic particles is denoted by fc, and q is the ratio of the stress to strain transfer, where σc, εc and σm, εm are the average stresses and strains in ceramic andmetal, respectively (see Fig. 3) q= σc−σm |εc−εm| 0 1, an increase in fc0 elevates vonMises stress in the outer surface, and the plastic region decreases in the inner surface. 122 A.T. Kalali et al. Fig. 4. Comparison of the VMP method (Haghpanah Jahromi et al., 2009) with the finite element analysis. The results show the residual stresses in the autofrettaged vessel. In this calculation: Em =56GPa,Ec =20GPa, σy0m =106MPa,Hm =12GPa and νm = νc =0.25. The vessel undergoes a plane-strain condition Fig. 5. VonMises stress along the thickness in the FGM cylindrical vessel subjected to internal pressure 300MPa with q=4.5GPa and n=2 for different fc0. The vessel has t/ri =1 and plane-stress condition (properties listed in Table 1) For the spherical vessel subjected to internal pressure of 600MPa, the results obtained from the isogeometric analysis method have been compared with the finite element results. Excellent agreement is observed as shown in Fig. 6. For the purpose to investigate the effect of n on the initiation of yielding,we introduce two parametersPv1 (pressures corresponding to the initiation of yielding at the inner radius) andPv2 (pressures corresponding to the initiation of yielding at the outer radius). Figure 7 shows thatPv1 andPv2 decreasewith growingn. Also, in the cylinder subjected to Pv2, by increasing n, the plastic region gradually spreads from the inner surface. Note that by increasing n, the metal properties dominate overcoming the ceramic properties and, therefore, the plastic behaviour of the material becomes more evident. Similar to the previous cases, the excellent agreement of the isogeometric analysis with the finite element predictions of elastic-plastic stresses for the rotating disc is shown in Fig. 8. The distribution of von Mises stresses across the thickness in the disk rotating at different angular velocities with n= 2 and fc0 =0.8 is presented in Fig. 9. The results clearly indicate that the growth of the plastic zone across the thickness is initiated from both sides of the disc. In this study, density of metal is ρm =4420kg/m 3 and density of ceramic is ρc =2000kg/m 3. Elastic-plastic analysis of pressure vessels and rotating disks... 123 Fig. 6. Comparison with the finite elementmethod. The results show the stress components in the spherical vessel with fc0 =1 and n=2 subjected to internal pressure 600MPa. In this calculation q=4.5GPa. The vessel has t/ri =1 (properties listed in Table 1) Fig. 7. VonMises stress along the thickness in the FGM spherical vessel subjected to Pv1 and Pv2 with fc0 =1 and different n. In this calculation q=4.5GPa. The vessel has t/ri =1 (properties listed in Table 1) Fig. 8. Comparison with the finite elementmethod. The results show the stress components in the rotating disk with fc0 =0.8 and n=2 at the angular velocity ω=230rad/s. In this calculation q=4.5GPa. The disk has t/ri =1 (properties listed in Table 1) 124 A.T. Kalali et al. Fig. 9. VonMises stress along the thickness in the FGM rotating disk at different angular velocities with fc0 =0.8 and n=2. In this calculation q=4.5GPa. The disk has t/ri =1 (properties listed in Table 1) 6. Conclusion Using the isogeometric analysis method, elastic-plastic stress distributions in a cylindrical and spherical pressurized vessels and rotating disksmade of anFGMmaterial have been determined. As expected, this approach to the plasticity problem is computationally cost effective and results in a much smaller system of equations to solve. Finite element analysis of the problem using ABAQUS commercial code has been used for verification of the isogeometric method. The nu- merical analysis within the software has been performedby the application of a “virtual thermal load” that enabled continuous variation of the material behaviour through the wall thickness. The analysis results obtained in this work also indicate the possibility of formation and growth of a plastic region within the wall thickness from the external surface of the FGM vessels or rotating diskswhereas in cylindrical (spherical) vessels and rotating disksmade of homogeneous materials, the plasticity essentially starts from the inner surface. References 1. Bayat M., Saleem M., Sahari B.B., Hamouda A.S.M., Mahdi E., 2008, Analysis of func- tionally graded rotating disks with variable thickness, Mechanics Research Communications, 35, 283-309 2. Carpenter R.D., Liang W.W., Paulino G.H., Gibeling J.C., Munir Z.A., 1999, Fracture testing and analysis of a layered functionally graded Ti/TiB beam in 3-point bending, Materials Science Forum, 308, 837-842 3. Chakraborty A., Gopalakrishnan S., Reddy J., 2003, A new beam finite element for the analysis of functionally gradedmaterials, International Journal ofMechanical Sciences,45, 519-539 4. Chen Y.Z., Lin X.Y., 2008, Elastic analysis for thick cylinders and spherical pressure vessels made of functionally gradedmaterials,Computational Materials Science, 44, 581-587 5. Cottrell J.A., HughesT.J.R., BazilevsY., 2009, Isogeometric Analysis: Toward Integration of CAD and FEA, JohnWiley, U.K. 6. Dai H.L., Fu Y.M., Dong Z.M., 2006, Exact solutions for functionally graded pressure vessels in a uniformmagnetic field, International Journal of Solids and Structures, 43, 5570-5580 7. Durodola J.F., AttiaO., 2006, Deformation and stresses in functionally graded rotating disks, Composites Science and Technology, 60, 987-995 Elastic-plastic analysis of pressure vessels and rotating disks... 125 8. Figueiredo F., Borges L., Rochinha F., 2008, Elasto-plastic stress analysis of thick-walled FGMpipes,AIP Conference Proceedings, 147-152 9. Haghpanah Jahromi B., Ajdari A., Nayeb-Hashemi H., Vaziri A., 2010, Autofrettage of layered and functionally graded metal-ceramic composite vessels, Composite Structures, 92, 1813-1822 10. Haghpanah Jahromi B., Farrahi G.H., Maleki M., Nayeb-Hashemi H., Vaziri A., 2009, Residual stresses in autofrettaged vessel made of functionally gradedmaterial,Engineering Struc- tures, 31, 2930-2935 11. Haghpanah Jahromi B., Nayeb-Hashemi H., Vaziri A., 2012, Elasto-plastic stresses in a functionally graded rotating disk, ASME Journal of Engineering Materials and Technology, 134, 21004-21015 12. Hassani B., Tavakkolin S.M., Moghadama N.Z., 2011, Application of isogeometric analysis in structural shape optimization, Scientia Iranica, 18, 864-852 13. Hughes T.J.R., Cottrell J.A., Bazilevs Y., 2005, Isogeometric analysis: CAD, finite ele- ments, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering, 194, 4135-4195 14. Jahed H., Dubey R.N., 1997, An axisymmetric method of elastic-plastic analysis capable of predicting residual stress field, Journal of Pressure Vessel Technology, 119, 264-273 15. Jahed H., Farshi B., Bidabadi J., 2005, Minimum weight design of inhomogeneous rotating discs, International Journal of Pressure Vessels and Piping, 82, 35-41 16. Jahed H., Farshi B., Karimi M., 2006, Optimum autofrettage and shrink-fit combination in multi-layer cylinders, Journal of Pressure Vessel Technology, 128, 196-201 17. Jahed H., Shirazi R., 2001, Loading and unloading behaviour of a thermoplastic disc, Interna- tional Journal of Pressure Vessels and Piping, 78, 637-645 18. Jin Z.H., Paulino G.H., Dodds Jr R.H., 2003, Cohesive fracture modelling of elastic-plastic crack growth in functionally gradedmaterials,Engineering Fracture Mechanics, 70, 1885-1912 19. Karlsson, Hibbitt, Sorensen, 2008, ABAQUS/CAE, v. 6.8-1. 20. Owen D.R.J., Hinton E. 1980, Finite Elements in Plasticity: Theory and Practice, Pineridge Press, U.K. 21. Sadeghian M., Ekhteraei H., 2011, Axisymmetric yielding of functionally graded spherical vessel under thermo-mechanical loading,Computational Materials Science, 50, 975-981 22. SetoodehA., Kalali A., Hosseini A., 2008,Numerical analysis of FGMplate by applying vir- tual temperature distribution,Proceedings of 7th Conference of Iranian Aerospace Society, Tehran 23. Suresh S., Mortensen A., 1998, Fundamentals of Functionally Graded Materials, IOM Com- munications Ltd. 24. You L.H., Zhang J.J., You X.Y., 2005, Elastic analysis of internally pressurized thick-walled sphericalpressurevessels of functionallygradedmaterials, International Journal of PressureVessels and Piping, 82, 374-345 Manuscript received December 24, 2014; accepted for print July 13, 2015