Theories and Analysis of Functionally Graded Beams

: This is a review paper containing the governing equations and analytical solutions of the classical and shear deformation theories of functionally graded straight beams. The classical, ﬁrst-order, and third-order shear deformation theories account for through-thickness variation of two-constituent functionally graded material, modiﬁed couple stress (i.e., strain gradient), and the von Kármán nonlinearity. Analytical solutions for bending of the linear theories, some of which are not readily available in the literature, are included to show the inﬂuence of the material variation, boundary conditions, and loads.


Preliminary Comments
Beams are structural members that have a ratio of length-to-height that is very large (say, a/h > 10) and are subjected to forces both in-plane and transverse to the plane that tend to bend about an axis perpendicular to their length. Such members are known as structural elements and their study constitutes structural mechanics, which is a subset of solid mechanics. Due to the particular structural configuration (i.e., one dimension being much bigger in comparison to the cross-sectional dimensions), the deformation and stress fields can be predicted, for most practical engineering problems, with structural theories known as the beam theories. Beam theories are derived from the three-dimensional elasticity theory by making certain simplifying assumptions concerning the deformation (kinematics) and stress states. The development of such theories dates back to Galileo Galilei, Leonardo da Vinci, and Jacob Bernoulli and Euler. The first one is the Euler-Bernoulli beam theory, in which the transverse shear strain is neglected, making the beam infinitely rigid in the transverse direction. The second one that accounts for the transverse shear strain (γ xz ) is popularly known as the Timoshenko beam theory [1,2].
All modern developments are refinements to the above stated two theories, where the displacement fields are expanded in terms of powers of the thickness [3] and accounting for other non-classical continuum mechanics aspects (e.g., stress and strain gradient effects and material length scales). For example, a general higher-order theory is of the form u = u xêx + u yêy + u zêz (1) where x (x), u y = 0, u z (x, z) = p ∑ i=0 z i ψ (i) z (x) (2) where φ (0) x = u and ψ (0) z = w are the midplane displacements along the x and z directions, respectively, and φ (i) x and ψ (i) x can be mathematically interpreted as higher-order generalized displacements with the meaning For a general third-order beam theory, we have m = 3 and p = 2 in Equation (2). The thirdorder beam theory of Reddy, derived from the third-order plate theory (see Reddy [4][5][6][7] and Heyliger and Reddy [8]), adopts a displacement field that is a special case of Equation (1) and imposes zero transverse shear stress conditions on the bounding planes (i.e., top and bottom faces) of the beam to express the variables introduced with the higher order terms in terms of the variables appearing in the Euler-Bernoulli and Timoshenko beam theories.

Functionally Graded Structures 1.2.1. Background
Functionally gradient materials (FGM) are a class of composites that have a gradual variation of material properties from one surface to another. These novel materials were proposed as thermal barrier materials for applications in space structures, nuclear reactors, turbine rotors, flywheels, and gears, to name only a few. In general, all the multi-phase materials, in which the material properties are varied gradually in a predetermined manner, fall into the category of functionally gradient materials. Such property enhancements endow FGMs with material properties such as the resilience to fracture.
In the last two decades, a large number journal papers dealing with functionally graded beams and plates have appeared in the literature and a critical review of these papers is not a focus of this introduction to FGM structures (see, e.g., Birman [9] and Klusemann [10] for a review). The works of Praveen and Reddy [11] also considered von Kármán nonlinearity in functionally graded plates.

FGM Material Models
A typical FGM represents a particulate composite with a prescribed distribution of volume fractions of constituent phases. In the case of beams, the material properties are assumed to vary continuously through the beam height. Several models are available in the literature, but the Voigt (or power-law) and Mori-Tanaka schemes [12] have been generally used for the study of FGM structures.
The advantage of the Voigt scheme is the simplicity of implementation and the ease of computation. According to the Voigt scheme, the effective properties are the arithmetic average of constituent property (P) and are given by (see, e.g., [11,13])

Modified Couple Stress Effects
Theories that account for microstructural length scales are the modified couple stress theory of Mindlin [14], Koiter [15], and Toupin [16] and the strain gradient theory of [17][18][19]. A more complete review of the early developments can be found in the work of Srinivasa and Reddy [20]. The strain gradient theory is a more general form of the modified couple stress theory and the relationship between the modified couple stress theory and the strain gradient theory can be found in the recent work of Reddy and Srinivasa [21]. In recent years a number of attempts have been made to bring microstructural length scales into the continuum description of beams and plates. Such models are useful in determining the structural response of micro and nano devices made of a variety of new materials that require the consideration of small material length scales over which the neighboring secondary constituents interact, especially when the spatial resolution is comparable to the size of the secondary constituents.
Microstructure-dependent theories are developed for the Bernoulli-Euler beam by Park and Gao [22], for the shear deformable beams and plates by Ma, Gao, and Reddy [23,24], and for vibration and buckling of shear deformable beams by Araujo dos Santos and Reddy [25][26][27]. In the last two decades, Reddy and his colleagues [23][24][25][26][27][28][29][30] have published a large number of papers dealing with linear and nonlinear bending of classical and firstand third-order shear deformable beams using the modified couple stress theory. Some of these works have accounted for the von Kármán nonlinearity and functionally graded materials. Of course, there are many papers by other colleagues on the same topics, which are not cited here and references to them can be found in the works already cited here. The von Kármán nonlinearity may have significant contribution to the response of beamlike elements used in micro-and nano-scale devices such as biosensors and atomic force microscopes (see, e.g., Li et al. [31] and Pei et al. [32]).

The Strain Energy Functional
The modified couple stress theory is based on the hypothesis that the rate of change of macrorotations cause additional stresses, called couple stresses, in the continuum. The rate of change of rotation is represented by the curvature tensor χ, which is defined by where ω is the rotation vector and u is the displacement vector of an arbitrary point in the beam. Physically, ω denotes the macrorotation at a point of the continuum. According to the modified couple stress theory [18], the strain energy potential of an elastic beam can be expressed as where L is the length of the beam, σ is the Cauchy stress tensor, ε is the simplified Green-Lagrange strain tensor, and m is the deviatoric part of the symmetric couple stress tensor. In the coming sections, these relations will be specialized to various beam theories. The couple stress tensor m is related to the curvature tensor χ through the constitutive relation [14]: where is the length scale parameterand G is the shear modulus. As pointed out in [33] the material length scale parameter of the modified couple stress theory is not constant for an especial material and changes as the size of a structure changes. To determine this value, experimental data for all different sizes are required. The present paper outlines the displacement fields of the three theories (classical, first-order, and third-order), the governing equations, and analytical solutions of straight beams for the linear case. To keep the size of the paper within reasonable limits, many details are not included here, and interested readers may consult the forthcoming book by the senior author on beams and circular plates [34], which is very comprehensive in its treatment of the theories, analytical solutions by exact means, the Navier solution approach, and numerical solutions by variational and finite element methods.

Kinematics
The displacement field of the classical beam theory (CBT) is constructed assuming that transverse lines perpendicular to the beam axis (x) remain: (1) straight, (2) inextensible, and (3) perpendicular to the tangents of the deflected x-axis. These assumptions, known as the Euler-Bernoulli hypothesis, result in the following displacement field: where (ê x ,ê z ) are the unit basis vectors along the (x, y) coordinates and (u, w) denote the axial and transverse displacements, respectively, of a point on the midplane of the beam. Based on the displacement field in Equation (9), the only nonzero strain in the present case is (see Reddy [35]) ε xx : The only nonzero components of the rotation and curvature are

Equations of Equilibrium
First we introduce the stress resultants N xx , M xx , and P xy where M xy is the couple stress induced by the difference between rates of rotations. Then using the principle minimum total potential energy, we obtain the Euler equations of equilibrium as whereM xx = M xx + P xy takes into account the modified couple stress effects in both the governing equations and the boundary conditions. The duality pairs of the CBT are (the first element of each pair is the primary variable and the second element is the secondary variable) where V eff is the effective shear force

Governing Equations in Terms of Displacements
For an isotropic material the one-dimensional stress-strain relation We assume that the beam is graded with two material combination through the beam height according to the relation where E 1 and E 2 are Young's moduli of the two materials used, and n is the index that dictates the relative dominance of volume fractions V 1 (z) and V 2 (z) = 1 − V 1 (z). We assume that Poisson's ratio ν is a constant for the FGM material. The stress resultants can be expressed in terms of the displacements as where A xx , B xx , D xx , and A xy are the extensional, extensional-bending, bending, and inplane shear stiffness coefficients. For beams with E = E(x) (i.e., n = 0 or E 1 = E 2 = E), we have A xx = EA 0 , B xx = 0, D xx = EI 0 , and A xy = GA 0 2 . For n = 0 (i.e., FGM beams), we have Figure 1a contains the variation of the non-dimensional axial stiffness (Ā xx = A xx /E 2 A 0 , A 0 = bh) and bending stiffnessD xx = D xx /E 2 I 0 , I 0 = bh 3 /12) as functions of the volume fraction index n for various values of the modulus ratio M = E 1 /E 2 ≥ 1 and Figure 1b contains similar plots for the non-dimensional extensional-bending coupling It is clear that bothĀ xx andD xx are the maximum at n = 0 and decrease with increasing values of n. However,B xx is zero at n = 0, increases to a maximum at n = √ 2, and then decreases with increasing value of n. Thus, beams with nonzero B xx will have a response that is not monotonic with respect to n. The equations of equilibrium can be expressed in terms of the displacements u and w using the beam constitutive relations from Equation (21). We obtain

General Solution
The exact solution to the linearized equations of equilibrium (i.e., static case) under distributed load q(x) is given by where K 1 through K 6 are constants of integration and ξ, η, ζ, and µ are dummy coordinates introduced to indicate the order of integration. The six constants of integration are determined using six boundary conditions, three at each end of the beam (i.e., one element of the each of the three duality pairs at each point: (u, N xx ), (w, V eff = dM xx /dx), and (dw/dx,M xx ). The stress resultants N xx andM xx can be computed using Equations (20a) and (20b). We can determine the constants of integration for various boundary conditions (pinned: u = w =M xx = 0; hinged: N xx = w =M xx = 0; and clamped u = w = dw/dx = 0). In the following we present the exact solutions for beams with various boundary conditions at x = 0 and x = L, L being the length of the beam.

Pinned-Hinged Beams
The exact solution for this case, with FGM and modified couple stress (MCS) effect, is We note that the effect of B xx on the mechanical deflection is zero, while the bending moment is independent of B xx .

Pinned-Pinned Beams
The solution for the pinned-pinned FGM beam is (ξ = x/L) When B xx = 0, the solutions for the pinned-hinged beams and pinned-pinned beams coincide.

Numerical Results
To present numerical results, we consider pinned-pinned functionally graded beams of length L = 100 in (254 cm), height h = 1 in (2.54 cm), and width b = 1 in (2.54 cm) and subjected to uniformly distributed load of intensity q 0 lb/in (1 lb/in = 175 N/m). The FGM beam is made of two materials with the following values of the moduli, Poisson's ratio, and shear correction coefficient: We shall investigate the parametric effects of the volume fraction index, n, and boundary conditions on the transverse deflections and bending moment.  Figure 1b) is reflected in the variation of u(x) with n. In particular, u(x) increases with increasing values of n for n < 5 but the magnitude of u(x) decreases with increasing values of n. On the other hand, the deflection and slope increase their magnitudes with the increasing values of n as the bending stiffness D xx dominates bending. Figure 4 contains variations of the maximum displacements u(0.25L) and w(0.5L) with the volume fraction index n. It is clear from the plots that the displacement u increases with n for n ≤ 5 and then decreases with the increasing values of n, as dictated by the variation of B xx with n (see Figure 1b). Both w(0.5L) and −(dw/dx)(L) (not shown here) increase with n but there are two parts, the first part exhibits rapid increase with n followed by slow increase due to the interplay between B xx and D xx in the bending response.

Clamped Beams
For beams clamped at both ends and subjected to uniformly distributed transverse load of intensity q(x) = q 0 , we havê The variation of u(x) for the clamped-clamped beam is the same as that of the pinnedpinned beam. Figure 5 contains plots of w(x) vs. x/L and −(dw/dx)(x) vs. x/L while Figure 6 contains the bending moment M xx (x) vs. x/L. The data used here is the same as that used for pinned-pinned beams. The results obtained show the same trends as those discussed for the pinned-pinned beams.

Preliminary Comments
In this section we consider the first-order shear deformation theory (TBT), most commonly known as the Timoshenko beam theory. The TBT brings the transverse shear strain γ xz = 2ε xz and shear stress σ xz into the calculations. However, in the TBT the transverse shear stress through the beam thickness is only represented as a constant, whereas the elasticity equations (as discussed in mechanics of materials books) show that the variation should be quadratic. To account for the inaccuracy in predicting the transverse shear force magnitude (not the shear stress distribution itself), shear correction factor (SCF) has been introduced (see [2,36,37], among many others). According to Timoshenko, the shear correction factor is the ratio of the average shear strain on a section to the shear strain at the geometric centroid of the cross section. The SCF, in general, is a function of the cross-sectional shape, Poisson's ratio, material properties, boundary conditions, and so on. For rectangular sections, Timoshenko [2] proposed a SCF K s = 5(1+ν) (6+5ν) , which takes the range values (5/6) ≤ K s ≤ (15/17) for 0 ≤ ν ≤ 0.5. Following these preliminary comments, we proceed, in somewhat parallel fashion to the developments presented for the classical beam theory (CBT), to develop the equations of equilibrium and exact solutions for the linear case.

Displacements and Strains
The TBT is based on the displacement field where φ x denotes the rotation (independent of the slope, θ x = −dw/dx) of the crosssectional plane about the y-axis. In the TBT, the normality assumption of the classical beam theory (CBT) is relaxed and a constant state of transverse shear strain (and thus constant shear stress computed from the constitutive equation) with respect to the thickness coordinate z is included. As stated earlier, the TBT requires a shear correction factor to compensate for the error due to this constant shear stress assumption. The von Kármán nonlinear strains of the TBT are where G the shear modulus [G(z) = E(z)/2(1 + ν)] and ν is Poisson's ratio, which is assumed to be a constant.

Equations of Equilibrium
The principle of minimum total potential energy for the TBT has the same form as that for CBT, except that one must add the strain energy terms associated with the transverse stress σ xz . The curvature in the TBT is given by The stress resultants of the TBT are defined as Here K s denotes the shear correction factor. The Euler equations of the TBT are The three duality pairs for the TBT are where the effective shear force and bending moments are We note that the effective shear force in the TBT has the modified couple stress term.

Governing Equations in Terms of Displacements Beam Constitutive Equations
The stress resultants (N xx , M xx , N xz , P xy ) in terms of the strains are where A xx , B xx , D xx , and A xy are as defined in Equation (21), and K s denotes the shear correction factor and S xz is the shear stiffness The equations of equilibrium in Equations (43)-(45) now can be expressed in terms of the displacements u, w, and φ x with the help of the beam constitutive relations in Equations (48a)-(48c) as

General Solution
In this section we present exact solutions to the linear equations of equilibrium of FGM beams without the modified couple stress effect. By setting f = 0 in Equations (50)-(52), we obtain Again, we further assume that the beam stiffness coefficients are all constant and f = 0. Equations (53)-(55), when expressed in terms of the stress resultants N xx , N xz , and M xx (see Equations (43)-(45) with P xy = 0) take the following form: Substituting for N xz from the third equation into the second equation, we obtain Integrating the above equations, Integrating the second equation, we obtain Here c i (i = 1, 2, 3) denote the constants of integration.
The left sides of Equations (58) and (59) can be expressed in terms of the displacements (u, φ x ) using Equations (48a) and (48b); we have Solving for du/dx and dφ x /dx, we obtain where Integrating Equations (61) and (62), we obtain From Equations (56) and (58), we arrive at and using Equation (48d) we obtain The six constants of integration are determined using six boundary conditions, three at each end of the beam. One (and only one) element of the each of three duality pairs at each boundary point must be known (see Equation (43)): (u, N xx ), (w, N xz ), and (φ x , M xx ). We note that, in TBT, φ x has replaced −dw/dx as the primary variable, and it is dual to the bending moment M xx . One should not specify dw/dx in place of φ x in the TBT. The stress resultants (N xx , M xx , P xy , N xz ) can be computed with the help of Equations (48a)-(48d).

Pinned-Pinned Beams
The exact solution of a beam pinned at both ends (u(0) = 0, w(0) = 0, M xx (0) = 0, u(L) = 0, w(L) = 0, and M xx (L) = 0) is It is clear that all functions except the transverse deflection are the same as those predicted by the classical beam theory (CBT). The transverse deflection has an additional positive term that adds to the value predicted by the CBT. Thus, the first-order shear deformation theory (TBT) deflection w is larger than those predicted by the CBT (i.e., the CBT underpredicts w). The numerical results obtained by the TBT are either the same or very close to those obtained using the CBT. Because of the fact that the beam considered is a thin beam with length-to-height ratio of L/h = 100, the effect of the shear deformation is not seen. Table 1 shows the numerical results obtained with the two theories for the data: Transverse deflections obtained with the CBT and TBT for three different length-toheight ratios, L/h = 50, L/h = 20, and L/h = 10 are presented in Table 2. It is clear that when the beam is moderately thick (L/h = 20) to very thick (L/h = 10), the CBT under predicts the deflections, although the difference between the two solutions may not be significant. Table 1. Numerical results obtained with the classical (CBT) and the shear deformation (TBT) beam theories for the displacementsū = u(0.25L) × 10 2 and w(0.5L) and slopes θ x (L) = −(dw/dx)(L) and φ x (L) of a pinned-pinned FGM beams under a uniformly distributed load (all results are normalized with the load).

Preliminary Comments
From the discussions presented in the previous sections, it is clear that the transverse shear stress distribution through the beam height, computed using the stress-strain relation σ xz = 2Gε xz , is either zero (in CBT) or constant (in TBT), although the actual variation of σ xz (x, z) with z, determined using the 3-D equations of equilibrium of linearized elasticity, is cubic. Therefore, it is necessary to have the displacement field (especially u 1 ) to be a cubic function of z. In this section, a third-order beam theory which accounts for the von Kármán geometric nonlinearity, through thickness variation of the material and modified couple stress effect, is presented. The theory presented herein accounts for the vanishing of transverse shear stress on the bottom and top surfaces of the beam (see [4,8,38,39]).

Kinematics
The displacement field of the Reddy-Bickford third-order beam theory (RBT) is where α = 4/3h 2 . The nonzero strain and curvature components are where (omitting the higher-order terms in the thickness strain)

Equations of Equilibrium
The equations of equilibrium of the RBT are derived using the principle of minimum total potential energy. We introduce the following stress resultants: The equations of equilibrium of the RBT are: The duality pairs (the first element of the pair denotes a generalized displacement while the second element denotes a generalized force): where Thus, there are four boundary conditions at each boundary point. Requiring ∂w/∂x as well as φ x to vanish at a support necessarily implies that the shear force, when shear stress is computed using the constitutive relation σ xz = Gγ xz , is zero. However, the effective shear force V eff is not zero. One of the challenges of higher-order theories is the ability to specify boundary conditions that involve higher-order stress resultants. In most cases, one does not know the known values of the higher-order stress resultants. Therefore, whenever the lower-order stress resultant is specified, we assume that the corresponding higher-order stress resultant is known to be zero. For example, if M xx is specified at a point, we assume thatM xx = M xx (implying that P xx = 0 there).

Beam Constitutive Relations
In the RBT, as in in the case of CBT and TBT, we invoked the inextensibility of the transverse normal lines, which amounts to setting ε zz = 0. Therefore, we can use one-dimensional constitutive relations. In particular, the one-dimensional constitutive relations are Substituting the constitutive relations from Equations (82) and (83) into the definition of the stress resultants in Equations (75a)-(75c), we obtain where (see Equation (21))

(86b)
If the higher-order terms are neglected in the governing equations of motion but not in the constitutive relations, we obtain the third-order theory developed by Levinson [38].
If one neglects the higher-order terms selectively in the constitutive relations, one obtains the so-called simplified Reddy-Bickford beam theory. These ideas will be discussed shortly.

Equilibrium Equations in Terms of the Displacements
With the help of Equations (84a) and (84b), the equations of equilibrium, Equations (76)-(78), can be expressed in terms of the generalized displacements (u, w, φ x ). We obtain

Exact Solutions for Bending
In this section we present exact solutions to the linear equations of equilibrium of the RBT for FGM beams without the effect of the modified couple stress. First, the equations of equilibrium in terms of the stress resultants can be obtained from Equations (87)-(89) by omitting the nonlinear terms and time-depedent terms: Integrating the equations with respect to x, we obtain where K 1 and K 2 are constants of integration. Expressing N xx and M xx in Equations (94) and (96) in terms of the generalized displacements, and solving for du/dx and dφ x /dx in terms of d 2 w/dx 2 , we obtain where (see Equations (87) and (88)) Integrating the two equations in (97a) and (97b), we obtain where K 4 and K 5 are the constants of integration. We return to Equation (95) and write it in terms of the generalized displacements (the differential of the constant part involving P 1 , P 2 , and K 1 is set to zero): where Integrating Equation (100) once and collecting the like terms together, we obtain where It is clear from Equation (102) that the analytical solution to the RBT is not algebraic but hyperbolic (because c 1 > 0 and c 2 > 0). The homogeneous solution of Equation (102) is w h (x) = K 7 cosh µx + K 8 sinh µx, µ = c 1 c 2 .
The total solution w(x) is obtained by adding the particular solution, w p (x) due to g(x): w(x) = w h (x) + w p (x). In addition, there are eight constants of integration, including the two constants of integration introduced in Equation (105). In the RBT, one is required specify −dw/dx in addition to φ x (or their dual variables, P xx andM xx , respectively), providing the required eight boundary conditions. The second-order derivative appearing in Equation (102) comes from P xx of Equation (92). If we neglect the second-order derivative of w in Equation (102), by reasoning that they are higher-order terms (i.e., very small compared to the other terms in the equation), we obtain Use of these boundary conditions yield u(x) = B xx D * xx q 0 L 3 12 w(x) = A xx D * xx q 0 L 4 24

Summary
In this paper three different beam theories, namely, the classical, first-order, and thirdorder beam theories are presented for beams, accounting for the through-thickness variation of the material, modified couple stress effect, and the von Kármán nonlinearity. Exact solutions for bending of the three theories are presented for several boundary conditions.
Numerical examples are also presented to illustrate the accuracy of various models and bring out certain salient features of functionally graded beams. A study of the FGM beams also revealed that the dimensionless bending deflections (w = wD xx /q 0 L 4 ) are not monotonic functions of the power-law index/exponent n because the coupling stiffness B xx is not a monotonically increasing or decreasing function of the modulus ratio.
Finite element models of the nonlinear theories presented herein can be found in the monograph by Reddy [34], which also contains detailed discussions of obtaining analytical and numerical solutions. A companion paper on FGM circular plates will appear following the publication of this paper (see Reddy, et al. [40]). Extensions of the theories presented herein to buckling and vibration [34,[41][42][43], and to account for nonlocal effects [44], are also awaiting. Extension of the theories and solutions to curved beams is another major topic for interested researchers.