Composite Structures 91 (2009) 385–390
Contents lists available at ScienceDirect
Composite Structures journal homepage: www.elsevier.com/locate/compstruct
A novel approach to stress analysis of pressurized FGM cylinders, disks and spheres Naki Tutuncu a,*, Beytullah Temel b a b
Cukurova University, Department of Mechanical Engineering, 01330 Balcali-Adana, Turkey Cukurova University, Department of Civil Engineering, 01330 Balcali-Adana, Turkey
a r t i c l e
i n f o
Article history: Available online 30 June 2009 Keywords: Functionally graded material Pressure vessel Disk Cylinder Sphere Complementary functions method
a b s t r a c t Axisymmetric displacements and stresses in functionally-graded hollow cylinders, disks and spheres subjected to uniform internal pressure are determined using plane elasticity theory and Complementary Functions method. The material is assumed to be functionally graded in the radial direction. Variations in the material properties such as Young’s modulus and Poisson’s ratio may be arbitrary functions of the radial coordinate. This assumption yields a two-point boundary value problem with a governing differential equation of variable coefficients. General analytical solutions of such equations are not available. Infusion of Complementary Functions method into the analysis of stresses in pressure vessels is a novel approach. Complementary Functions method reduces the boundary value problem to an initial-value problem which can be solved accurately by one of many efficient methods such as Runge–Kutta method. Various material models from the literature are used and corresponding radial displacement and stresses are computed. Verification of the proposed method is done using benchmark solutions available in the literature for some special cases and virtually exact results are obtained using the fifth-order Runge– Kutta method (RK5). Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Recent developments in materials engineering have allowed the use of inhomogeneous materials called functionally graded materials (FGMs) in structural elements. Functionally graded materials are heterogeneous composite materials with gradient compositional variation of the constituents (e.g. metallic and ceramic) from one surface of the material to the other which results in continuously varying material properties. The materials are intentionally designed in such a way that they possess desirable properties for specific applications. The concept of FGM, initially developed for super heat resistant materials to be used in space planes or nuclear fusion reactors, is now of interest to designers of functional materials for energy conversion, dental and orthopedic implants, sensors and thermo generators [1]. Additional potential applications of FGMs include their use as interfacial zones to improve the bonding strength and to reduce residual stresses in bonded dissimilar materials and as wear resistant layers such as gears, cams, ball and roller bearings and machine tools [2]. The variation of material properties in an FGM is usually obtained by changing the volume fractions of two or more compatible constituents. FGMs with material properties varying in the thickness direction can be manufactured by high speed centrifugal casting or by depositing ceramic layers on a metallic substrate. Functionally graded fiber * Corresponding author. E-mail address:
[email protected] (N. Tutuncu). 0263-8223/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2009.06.009
reinforced composites can be manufactured by varying the volume fraction and/or the orientation of fibers in the preform prior to infusing resin into it [3]. Increasing use of functionally graded pressure vessels has necessitated comprehensive stress analyses. Horgan and Chan [4] have presented exact static solutions for the pressurized hollow cylinder or disk problem for functionally graded isotropic linearly elastic materials. A similar study for FGM pressure vessels in a uniform magnetic field is carried out by Jabbari et al. [5,6], Shao [7] and Dai et al. [8]. Jabbari et al. [9] have also obtained mechanical and thermal stress solutions in functionally-graded hollow cylinders. Stresses in disks involve similar analyses and works by Reddy and Srinath [10] and Gurushankar [11] may be cited among the earlier studies conducted on the subject. Concerning the stress analysis of cylindrical and spherical structural elements Fukui et al. [12] and Salzar [13] have presented analyses involving finite elements and other numerical techniques due to the nature of functions chosen to describe the inhomogeneous material properties. Analytical solutions to benchmark problems provide invaluable checks on the accuracy of numerical or approximate analytical schemes and allow for widely applicable parametric studies [4]. Assuming constant Poisson’s ratio and Young’s modulus obeying a simple power law of the radial coordinate renders the standard Euler–Cauchy equation whose solution is readily available. These solutions are presented in the Appendix for verification purposes. Employment of sophisticated micromechanical methods such as the Mori–Tanaka and the self-consistent methods to estimate the
386
N. Tutuncu, B. Temel / Composite Structures 91 (2009) 385–390
effective material properties yields more complex functions for Poisson’s ratio and Young’s modulus [14]. The resulting governing differential equation then possesses variable coefficients. General closed-form solutions of such equations are not available. In such situations solution methods include integral transformations, development of finite element models and, in some special cases, series solutions are attempted. Assuming the member is composed of many homogeneous layers of different properties emulating the FGM behavior is another way of tackling the problem on hand. All of these approaches require heavy mathematical manipulations and, in the case of having to discretise the domain into many elements, high amount of computational time. Linear elastic analysis of isotropic pressure vessels with axisymmetric conditions essentially results in a two-point boundary value problem. The governing differential equation will have variable coefficients which are functions of material properties. A novel approach will be attempted to obtain displacements, strains and stresses in a simple and efficient manner: the Complementary Functions method will be infused into the analysis to convert the problem to an initial-value problem which can efficiently be solved by, for example, the fifth-order Runge–Kutta method (RK5) with great accuracy. The theoretical background for the method is available in the literature [15–17]. The method is also successfully applied in other structural mechanics problems such as those involving curved bars [18] and beams [19]. Material properties for typical FGMs made from a mixture of ceramic and metal will be used in the analysis. The ceramic constituent is able to withstand high temperature environments due to its better thermal resistance, while the metal constituent provides stronger mechanical performance. Three material models will be used: (a) simple power law with constant Poisson’s ratio [4] for which case analytical benchmark solutions are available, (b) exponentially-varying properties [20], and (c) volume fractions of the constituents varying according to a power law [21,22]. It should be emphasized once again that the solution procedure is not confined to any particular choice of material model; it is equally suitable for arbitrary functions defining the gradient variation of material properties.
for the cylinder in plane strain. An annular disk is essentially a thin circular plate and plane stress assumptions will be valid. The stiffness terms will then take the following form:
C 11 ¼
E ; 1 m2
C 12 ¼
Em 1 m2
Note that for functionally graded materials E, m and, thus, the stiffness terms will be functions of the non-dimensional radial coordinate x. The only non-trivial equilibrium equation is
drr rr rh þ ¼0 dx x
ð3Þ
Using Eqs. (1)–(3) the governing equation of radial displacement becomes 2
d V 2
dx
þ MðxÞ
dV þ NðxÞV ¼ 0 dx
ð4Þ
where
MðxÞ ¼
1 dC 11 1 þ C 11 dx x
and
1 1 dC 12 1 x C 11 dx x
NðxÞ ¼
Considering the fact that stiffness terms C11 and C12 will be functions of x, a general closed-form solution of Eq. (4) is not possible. For a uniform internal pressure p the boundary conditions are given as
rr jx¼1 ¼ p and rr jx¼k ¼ 0
ð5Þ
2.2. FGM spheres The strain–displacement relations for a sphere are
er ¼
dV dx
and
eh ¼ e/ ¼
V x
ð6Þ
2. Governing equations
Under axisymmetric conditions the thick-walled sphere will have two independent stress components. The stress–strain relations are
Consider a thick-walled hollow cylinder (or a disk) and sphere of inner radius a and outer radius ka where k is a constant. The material is isotropic and functionally graded. Modulus of elasticity E and Poisson’s ratio m vary in the radial direction.
rr ¼ C 11 er þ 2C 12 eh rh ¼ r/ ¼ C 12 er þ ðC 11 þ C 12 Þeh
2.1. FGM cylinders and disks The axisymmetric governing equations to be presented in what follows are applicable to both thick-walled hollow cylinders and annular disks. The strain–displacement relations in term of the normalized radial displacement V are
er ¼
dV dx
and
eh ¼
V x
ð1Þ
2
d V 2
þ Q ðxÞ
dV þ RðxÞV ¼ 0 dx
where
1 dC 11 2 þ C 11 dx x
and
C 12 ¼
Em ð1 þ mÞð1 2mÞ
ð8Þ
Using Eqs. (6)–(8) the displacement equation is derived as
ð2Þ
where C11 and C12 are stiffness terms given by
Eð1 mÞ ; ð1 þ mÞð1 2mÞ
drr 2ðrr rh Þ ¼0 þ x dx
QðxÞ ¼
rr ¼ C 11 er þ C 12 eh rh ¼ C 12 er þ C 11 eh
C 11 ¼
Application of general Hooke’s law in terms of Lame constants results in stiffness terms C11 and C12 being the same as in the cylinder in plane strain. Both problems describe pure radial expansion [23]. The equilibrium equation is
dx
where V and x are non-dimensional variables defined as V = u/a and x = r/a with u being the radial displacement and r the radial coordinate. The stress–strain relations are
ð7Þ
RðxÞ ¼
2 1 dC 12 1 x C 11 dx x
subjected to the same boundary conditions given by Eq. (5).
ð9Þ
387
N. Tutuncu, B. Temel / Composite Structures 91 (2009) 385–390
Here, too, a general closed-form solution is not possible. Benchmark solutions for cylinders (or disks) for special cases of material models (homogenous and FGM obeying a single power law) are derived and given in the Appendix. 3. Solutions by the Complementary Functions method The method of Complementary Functions (CFM) transforms the solution of two-point boundary value problems (BVP) to initial-value problems (IVP). It reduces to a particularly simple solution scheme when applied to the present class of problems. 3.1. Cylindrical vessels (or disks) The governing equation of displacement is rewritten as 00
V þ MV 0 þ NV ¼ 0
ð10Þ
0
where () denotes the derivative with respect to x, and M and N are functions of x. The boundary conditions given by Eq. (5) written in terms of displacement take the form
V þ C 12 ¼ p x¼1 x x¼1 V C 11 V 0 x¼k þ C 12 ¼0 x
C 11 V 0
ðiÞ
ð14aÞ
and
¼
ðiÞ MZ 2
ðiÞ NZ 1
ð14bÞ
The above system of equations can be solved numerically for the homogeneous solutions. Kronecker delta initial conditions given below will be used to assure the linear independence of the solutions [16] ðjÞ
Z i ¼ dij ;
i; j ¼ 1; 2
ð15Þ
The numerical scheme chosen is the fifth-order Runge–Kutta method (RK5). The solutions are calculated at 11 collocation points through the thickness (10 divisions) in the interval 1 6 x 6 k. Having determined the homogeneous solutions Vi and their first derivatives V 0i , the constants bi can now be found by imposing the boundary conditions given by Eqs. (11) and (12). This process results in the following simple system of equations for the constants b1 and b2:
A11
A12
A21
A22
#(
b1 b2
(
)
¼
p
)
ð16Þ
0
where
C 12 V 1 C 11 V 01 þ x x¼1 C 12 ¼ C 11 V 01 þ V 1 x x¼k
A11 ¼ A21
Eð1 mÞ ; ð1 þ mÞð1 2mÞ
C 12 ¼
Em ð1 þ mÞð1 2mÞ
and the coefficients in Eq. (14b) take the form
1 ; x
N¼
1 x2
m ¼ lnðEc =Em Þ= ln k
ðiÞ ðZ 2 Þ0
"
C 11 ¼
Using the properties E = 200 GPa, m = 0.3 and k = 2, the radial displacement V, radial stress rr and hoop stress rh are calculated at 11 collocation points through the thickness (10 divisions) and these results are compared with those obtained from analytical solutions as shown in Table 1. Another verification attempt has been made regarding the FGM cylinder with constant Poisson’s ratio and Young’s modulus varying as E(x) = Emxm, where now Em is the metal phase on the inner surface (x = 1). On the outside (x = k), ceramic phase is present with modulus Ec. The exponent m is then expressed as
where Vi are the linearly independent homogeneous solutions, and bi are constants to be determined via the boundary conditions. ðiÞ ðiÞ Letting V i ¼ Z 1 and V 0i ¼ Z 2 yields ðiÞ
Analytical solutions are readily available for homogeneous cylinders and FGM cylinders with constant Poisson’s ratio and elastic modulus obeying a simple power law as E(x) = Emxm, where Em is the modulus of the constituent on the inner surface [4]. These benchmark solutions are given in the Appendix expressed using the notation of the present study. For a homogeneous cylinder with constant modulus E and Poisson’s ratio m the stiffness terms are
ð12Þ
ð13Þ
ðZ 1 Þ0 ¼ Z 2
4. Verification of the proposed solution scheme
M¼
Solution of Eq. (10) is
i ¼ 1; 2
The only difference from the cylindrical case is the variable coefficients in the governing differential equation. In the spherical case the variable coefficient Q and R will take the place of N and M, respectively, in Eqs. (10) and (14b). The complexity of these functions will entirely depend upon the material models used for FGM behavior.
ð11Þ
x¼k
V ¼ bi V i ;
3.2. Spherical vessels
C 12 V 2 A12 ¼ C 11 V 02 þ x x¼1 C 12 A22 ¼ C 11 V 02 þ V 2 x x¼k
Note that since the first derivatives of the radial displacement are readily available calculation of the stress values will not require additional computations.
ð17Þ
The stiffness terms take the form
C 11 ¼
Em ð1 mÞ xm ; ð1 þ mÞð1 2mÞ
C 12 ¼
Em m xm ð1 þ mÞð1 2mÞ
The variable coefficients now become
1 M ¼ ðm þ 1Þ ; x
N¼
mm 1 1 2 1m x
For numerical results the values Em = 200 GPa, Ec = 360 GPa, m = 0.3 and k = 2 have been used and these results are compared with the exact ones as given in Table 2. Table 3 gives a summary of the convergence study performed on the results of this case. Relative error and is calculated at the midpoint through is defined as AnalyticalCFM Analytical the thickness (x = 1.5) for radial displacement, radial stress and hoop stress values. For number of divisions greater than two the relative error decreases very fast. In fact, at 10 divisions the relative error is in the order of 107 for all results. Note that the analytical (exact) results listed in Tables 1 and 2 are results obtained from the analytical benchmark solutions presented in the Appendix for special cases. Examining these results reveals the great accuracy and efficiency achieved by CFM; calculations performed at only 11 points through the thickness yielded exact numerical results. 5. Analysis of other FGM models Another FGM model is taken from Ref. [21]. The region between the inner and outer surfaces is made of the combined
388
N. Tutuncu, B. Temel / Composite Structures 91 (2009) 385–390
Table 1 Comparison of CFM results with analytical results for homogeneous cylinder. x
rr (GPa)
V
1. 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.
rh (GPa)
CFM
Analytical
CFM
Analytical
CFM
Analytical
0.009533 0.008832 0.008262 0.007793 0.007404 0.007078 0.006803 0.006571 0.006375 0.006208 0.006067
0.009533 0.008832 0.008262 0.007793 0.007404 0.007078 0.006803 0.006571 0.006375 0.006208 0.006067
1. 0.768595 0.592593 0.455621 0.346939 0.259259 0.187499 0.128028 0.078189 0.036011 0.
1. 0.768595 0.592593 0.455621 0.346939 0.259259 0.1875 0.128028 0.078189 0.036011 0.
1.666667 1.435262 1.259259 1.122288 1.013605 0.925926 0.854167 0.794694 0.744856 0.702678 0.666667
1.666667 1.435262 1.259259 1.122288 1.013605 0.925926 0.854167 0.794694 0.744856 0.702678 0.666667
Table 2 Comparison of CFM results with analytical results for FGM cylinder using Model 1. x
rr (GPa)
V
1. 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.
Analytical
CFM
Analytical
CFM
Analytical
0.007422 0.006811 0.006328 0.005940 0.005624 0.005365 0.005150 0.004971 0.004821 0.004694 0.004587
0.007422 0.006811 0.006328 0.005940 0.005624 0.005365 0.005150 0.004971 0.004821 0.004694 0.004587
1. 0.803165 0.644350 0.513670 0.404335 0.311538 0.231791 0.162507 0.101729 0.047951 0.
1. 0.803165 0.644350 0.513670 0.404335 0.311538 0.231791 0.162507 0.101729 0.047951 0.
1.202558 1.131098 1.076507 1.034230 1.001160 0.975110 0.954510 0.938199 0.925311 0.915184 0.907305
1.202558 1.131098 1.076507 1.034230 1.001160 0.975110 0.954510 0.9382 0.925311 0.915184 0.907305
Table 3 Relative error values at x = 1.5 for Model 1 (n: number of divisions). n
V
rr
rh
2 4 6 8 10
2.60E04 1.20E05 1.70E06 4.00E07 1.30E07
2.10E04 1.00E05 1.40E06 3.40E07 1.10E07
2.70E04 1.20E05 1.70E06 4.00E07 1.30E07
ceramic–metal material with different mixing ratios of the ceramic and metal. Defining the volume fraction of the ceramic constituent as
x1 Vc ¼ k1
n ð18Þ
The overall Poisson’s ratio m and Young’s modulus E of the FGM can be given by
m ¼ mc V c þ mm ð1 V c Þ E ¼ Ec V c þ Em ð1 V c Þ
ð19Þ ð20Þ
An additional model to be used of a metal–ceramic FGM is given in [20]. Overall Poisson’s ratio and modulus of elasticity are assumed to vary as
ð1 mc Þecx 1 þ bx 2ð1 þ bxÞ ð1 mc Þecx
m¼1 E ¼ Ec
ð1 þ mc Þð1 þ bxÞ2
ð21Þ ð22Þ
where the constants b and c are given by
b¼
rh (GPa)
CFM
lc 1 lm
c ¼ lnð1 þ bÞ þ ln
ð23Þ 1 þ mm 1 mc
ð24Þ
In the above equations Ec, mc, lc are Young’s modulus, Poisson’s ratio and the shear modulus of the ceramic constituent and Em, mm, lm are those of the metal constituent. The material properties given by these models are substituted first in the stiffness terms C11 and C12. M and N coefficients for cylinders and disks and Q and R for spheres are subsequently calculated. These coefficients would yield extremely complex expressions and they will not be explicitly given in the present paper. 6. Results and conclusions For numerical results the properties used are Em = 200 GPa, Ec = 360 GPa, mm = 0.333, mc = 0.2, lm = 75 GPa, lc = 150 GPa. The exponent in Eq. (18) is taken to be n = 2 and the ratio of outer to inner radius for all cases is k = 2. The results for non-dimensional radial displacement V, radial stress rr and hoop stress rh are all calculated for a unit inner pressure and presented in the form of graphs. The FGM model with constant Poisson’s ratio and Young’s modulus varying as a single power law is referred to as Model 1. The properties varying according to the rule of mixtures of volume fractions are used in Model 2. Model 3 includes overall Poisson’s ratio and modulus of elasticity varying as exponential functions of the radial coordinate as given by Eqs. (21) and (22). Different models result in substantial differences in the radial displacement distribution through the thickness. Model 2 yields the highest and Model 3 the lowest displacements in all cases as Figs. 1, 4 and 7 show. Hoop stress values converge to about the same magnitude at around midpoint of thickness. Hoop stress and radial displacement magnitudes are about the same in cylinder and disk but much less in the sphere. The order does not change with respect to the models used. Fig. 2 shows that, in the case of the cylinder, the radial stress due to Model 3 lie between those of Models 1 and 2. As seen in Figs. 5 and 8 Models 1 and 2 yield radial stresses of about the same magnitude in disk and sphere. Models 1 and 2
389
N. Tutuncu, B. Temel / Composite Structures 91 (2009) 385–390
Model 1
0.008
Model 2 Model 3
0.006
σr (GPa)
V
0.007
0.005 0.004 0.003 1
1.2
1.4
1.6
1.8
2
0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
Model 1 Model 2 Model 3
x
1
1.2
1.4
x
1.6
1.8
2
Fig. 1. Radial displacements in cylinder for various material models. Fig. 5. Radial stresses in disk for various material models.
Model 1
1.60
Model 2
1.40
-0.4 -0.5 -0.6 -0.7 -0.8
σθ (GPa)
σr (GPa)
0 -0.1 -0.2 -0.3
Model 1 Model 3
1
1.2
1.4
x
1.6
1.8
1.20 1.00 0.80
Model 2
-0.9 -1
Model 3
0.60 1
2
Fig. 2. Radial stresses in cylinder for various material models.
1.2
1.4
1.6
x
1.8
2
Fig. 6. Hoop stresses in disk for various material models.
0.0037
1.60
0.0032
Model 2 Model 3
Model 2 Model 3
0.0027 V
1.40
σθ (GPa)
Model 1
Model 1
1.20
0.0022
1.00
0.0017
0.80
0.0012 0.0007
0.60 1
1.2
1.4
x
1.6
1.8
2
1
Fig. 3. Hoop stresses in cylinder for various material models.
0.0083
Model 2 Model 3
V
σr (GPa)
0.0063 0.0053 0.0043 0.0033 1
1.2
1.4
x
1.6
1.8
Fig. 4. Radial displacements in disk for various material models.
1.4
x
1.6
1.8
2
Fig. 7. Radial displacements in sphere for various material models.
Model 1
0.0073
1.2
2
0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1
Model 1 Model 2 Model 3
1
1.2
1.4
x
1.6
1.8
Fig. 8. Radial stresses in sphere for various material models.
2
390
N. Tutuncu, B. Temel / Composite Structures 91 (2009) 385–390
k
Model 2
V¼
Model 3
0.59
σθ (GPa)
Radial displacement
Model 1
0.69
k
k 2 xk 1 k1
ðk k
k2
0.49
Radial stress
0.39
rr ¼
ÞðC 012
þ
k
0.29
C 011 k1 Þ
k 1 xk 2 ðk k ÞðC 012 þ C 011 k2 Þ
1
1.2
1.4
x
1.6
1.8
2
k
xm1 ðk 2 xk1 k 1 xk2 Þ k
k 1 k
ðA:5Þ
k2
rh ¼
xm1 ðS1 þ S2 þ S3 Þ ðk k ÞðC 012 þ C 011 k1 ÞðC 012 þ C 011 k2 Þ k1
k2
Appendix A. Analytical benchmark solutions The following expressions are derived for a unit inside pressure. A.1. Homogeneous cylinder Radial displacement
x 2
ðk 1ÞðC 12 þ C 11 Þ
2
k x1 2
ðk 1ÞðC 12 C 11 Þ
ðA:1Þ
Radial stress 2
1 k x2 2
k 1
ðA:2Þ
Hoop stress 2
k x2 þ 1 2
k 1
ðA:3Þ
A.2. FGM cylinder using Model 1 Poisson’s ratio is assumed constant and Young’s modulus varies as E = Emxm. Define C 011 ; C 012 ; m ; k1 and k2 as
Em ð1 mÞ ð1 þ mÞð1 2mÞ Em m ¼ ð1 þ mÞð1 2mÞ
C 011 ¼ C 012
m ¼
m
1m pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 k1 ¼ ðm 4 þ m2 4mm Þ 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 k2 ¼ ðm þ 4 þ m2 4mm Þ 2
ðA:6Þ
S1 ¼ C 011 C 012 ð1 þ k1 k2 Þðk 2 xk1 k 1 xk2 Þ k
yield similar results because they essentially include the material properties varying according to a power law. Model 3, however, gives noticeably different values due to exponential variation of properties.(see Figs. 3, 6 and 9). It should be pointed out once again that the purpose of the present work has been the introduction of CFM to the solution procedure of the class of problems on hand. Converting the two-point boundary value problem to an initial-value problem gave way to the implementation of well-established numerical schemes. Runge–Kutta method of fifth-order (RK5) has been used to solve the system of equations [24]. The procedure is simple, efficiently implemented and accurate to the extent that exact numerical results have been obtained to six-digit accuracy by picking only 11 collocation points in RK5.
rh ¼
ðA:4Þ
where
Fig. 9. Hoop stresses in sphere for various material models.
rr ¼
k2
Hoop stress
0.19
V¼
k1
S2 ¼ S3 ¼
k ðC 011 Þ2 ðk 2 k2 xk1 k ðC 012 Þ2 ðk 2 k1 xk1
k1
k
k2
k k1 x Þ k
k 1 k 2 xk 2 Þ
References [1] Kashtalyan M. Three-dimensional elasticity solution for bending of functionally graded rectangular plates. Eur J Mech A – Solid 2004;23:853–64. [2] Erdogan F. Fracture mechanics of functionally graded materials. Compos Eng 1995;5(7):753–70. [3] Qian LF, Batra RC, Chen LM. Static and dynamic deformations of thick functionally graded elastic plates by using higher-order shear and normal deformable plate theory and meshless local Petrov–Galerkin method. Compos Part B – Eng 2004;35:685–97. [4] Horgan CO, Chan AM. The pressurized hollow cylinder or disk problem for functionally graded isotropic linearly elastic materials. J Elasticity 1999;55(1):43–59. [5] Jabbari M, Sokrabpour S, Eslami MR. Mechanical and thermal stresses in a functionally graded hollow cylinder due to radially symmetric loads. Int J Pres Ves Pip 2002;79(7):493–7. [6] Jabbari M, Sokrabpour S, Eslami MR. General solution for mechanical and thermal stresses in a functionally graded hollow cylinder due to nonaxisymmetric steady-state loads. J Appl Mech – Trans ASME 2003;70(1):111–8. [7] Shao ZS. Mechanical and thermal stresses of a functionally graded circular hollow cylinder with finite length. Int J Pres Ves Pip 2005;82(3):155–63. [8] Dai HL, Fu YM, Dong ZM. Exact solutions for functionally graded pressure vessels in a uniform magnetic field. Int J Solids Struct 2006;41(18– 19):5570–80. [9] Jabbari M, Bahtui A, Eslami MR. Axisymmetric mechanical and thermal stresses in thick long FGM cylinders. J Therm Stresses 2006;29(7):643–63. [10] Reddy TY, Srinath H. Elastic stresses in rotating anisotropic annular disk of variable thickness and variable density. Int J Mech Sci 1974;16:85–9. [11] Gurushankar GV. Thermal stresses in a rotating, nonhomogeneous anisotropic disk of varying thickness and density. J Strain Anal 1975;10(3):137–42. [12] Fukui Y, Yamanaka N. Elastic analysis for thick-walled tubes of functionally graded materials. JME Int J Ser I: Solid Mech, Strength Mater 1992;35(4):379–85. [13] Salzar RS. Functionally graded metal matrix composite tubes. Compos Eng 1995;5(7):891–900. [14] Vel SS, Batra RC. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. J Sound Vib 2004;272:703–30. [15] Aktas Z. Numerical solutions of two-point boundary value problems. Ankara, Turkey: METU, Dept of Computer Eng; 1972. [16] Roberts SM, Shipman JS. Fundamental matrix and two-point boundary-value problems. J Optimiz Theory Appl 1979;28(1):77–8. [17] Agarwal RP. On the method of complementary functions for nonlinear boundary-value problems. J Optimiz Theory Appl 1982;36(1):139–44. [18] Yildirim V. Free vibration analysis of non-cylindrical coil springs by combined use of the transfer matrix and the complementary functions methods. Commun Numer Methods Eng 1997;13:487–94. [19] Calim FF. Free and forced vibrations of non-uniform composite beams. Compos Struct 2009;88:413–23. [20] Jin ZH, Batra RC. Some basic fracture mechanics concepts in functionally graded materials. J Mech Phys Solids 1996;44:1221–35. [21] Liew KM, He XQ, Ng TY, Kitipornchai S. Finite element piezothermoelasticity analysis and the active control of FGM plates with integrated piezoelectric sensors and actuators. Comput Mech 2003;31(3-4):350–8. [22] Ferreira AJM, Batra RC, Roque CMC, Qian JF, Jorge RMN. Natural frequencies of functionally graded plates by a meshless method. Compos Struct 2006;75:593–600. [23] Love AEH. A treatise on the mathematical theory of elasticity. New York: Dover Publications; 1944. p. 142. [24] Chapra SC, Canale RP. Numerical methods for engineers. New York: McGraw-Hill; 1998. p. 703.