Theories and Analysis of Functionally Graded Beams

Abstract

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, first-order, and third-order shear deformation theories account for through-thickness variation of two-constituent functionally graded material, modified 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 influence of the material variation, boundary conditions, and loads.

1. Introduction

1.1. 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 (${\gamma }_{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

$$\mathbf{u}=u_x{\widehat{\mathbf{e}}}_x+u_y{\widehat{\mathbf{e}}}_y+u_z{\widehat{\mathbf{e}}}_z$$

(1)

where

$$u_x(x,z)=\sum _{i=0}^m{z}^i{\varphi }_x^{\left(i\right)}\left(x\right),u_y=0,u_z(x,z)=\sum _{i=0}^p{z}^i{\psi }_z^{\left(i\right)}\left(x\right)$$

(2)

where ${\varphi }_x^{\left(0\right)}=u$ and ${\psi }_z^{\left(0\right)}=w$ are the midplane displacements along the x and z directions, respectively, and ${\varphi }_x^{\left(i\right)}$ and ${\psi }_x^{\left(i\right)}$ can be mathematically interpreted as higher-order generalized displacements with the meaning

$${\varphi }_x^{\left(i\right)}=\frac{1}{\left(i\right)!}{\left(\frac{{\partial }^i{u}_x}{\partial z^i}\right)}_{z=0},{\psi }_z^{\left(i\right)}=\frac{1}{\left(i\right)!}{\left(\frac{{\partial }^i{u}_z}{\partial z^i}\right)}_{z=0}$$

(3)

For a general third-order beam theory, we have $m=3$ and $p=2$ in Equation (2). The third-order 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.

1.2. 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.

1.2.2. 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])

$$P\left(z\right)=\left(P_1-P_2\right)f\left(z\right)+P_2,f\left(z\right)={\left(\frac{1}{2}+\frac{z}{h}\right)}^n$$

(4)

1.3. Modified Couple Stress Effects

1.3.1. Background

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 first- and 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 beam-like 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]).

1.3.2. 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 $\chi$, which is defined by

$$\chi =\frac{1}{2}\left[\mathbf{\nabla }\omega +{\left(\mathbf{\nabla }\omega \right)}^{\mathrm{T}}\right]$$

(5)

where $\omega$ is the rotation vector

$$\omega =\frac{1}{2}\mathbf{\nabla }\times \mathbf{u}$$

(6)

and u is the displacement vector of an arbitrary point in the beam. Physically, $\omega$ 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

$$U=\frac{1}{2}{\int }_0^L\left[{\int }_A\left(\sigma :\epsilon +\mathbf{m}:\chi \right)dA\right]dx$$

(7)

where L is the length of the beam, $\sigma$ is the Cauchy stress tensor, $\epsilon$ is the simplified Green–Lagrange strain tensor, and $\mathbf{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 $\chi$ through the constitutive relation [14]:

$$\mathbf{m}=2G{\ell }^2\chi$$

(8)

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.

2. Classical Theory of Beams (CBT)

2.1. 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:

$$\mathbf{u}(x,z)=\left[u\left(x\right)+z{\theta }_x\left(x\right)\right]{\widehat{\mathbf{e}}}_x+w\left(x\right){\widehat{\mathbf{e}}}_z,{\theta }_x\equiv -\frac{dw}{dx},$$

(9)

where $({\widehat{\mathbf{e}}}_x,{\widehat{\mathbf{e}}}_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]) ${\epsilon }_{xx}$:

$$\begin{array}{cc}\hfill {\epsilon }_{xx}(x,z)& =\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2+z\left(-\frac{d^{2}w}{d{x}^2}\right)\equiv {\epsilon }_{xx}^{\left(0\right)}+z{\epsilon }_{xx}^{\left(0\right)},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\epsilon }_{xx}^{\left(0\right)}\left(x\right)& =\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2,{\epsilon }_{xx}^{\left(1\right)}\left(x\right)=-\frac{d^{2}w}{d{x}^2}.\hfill \end{array}$$

(11)

The only nonzero components of the rotation and curvature are

$${\omega }_y=\frac{1}{2}\left(\frac{\partial u_1}{\partial z}-\frac{\partial u_3}{\partial x}\right)=-\frac{dw}{dx},{\chi }_{xy}=\frac{1}{2}\frac{\partial {\omega }_y}{\partial x}=-\frac{1}{2}\frac{d^{2}w}{d{x}^2}$$

(12)

2.2. Equations of Equilibrium

First we introduce the stress resultants $N_{xx}$, $M_{xx}$, and $P_{xy}$

$$N_{xx}={\int }_A{\sigma }_{xx}dA,M_{xx}={\int }_{A}z{\sigma }_{xx}dA,P_{xy}={\int }_A{\mathcal{M}}_{xy}dA,$$

(13)

where ${\mathcal{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

$$\begin{array}{cc}\hfill \frac{d{N}_{xx}}{dx}+f& =0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{d^2{\overline{M}}_{xx}}{d{x}^2}+\frac{d}{dx}\left(\frac{dw}{dx}{N}_{xx}\right)+q& =0,\hfill \end{array}$$

(15)

where ${\overline{M}}_{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)

$$(u,N_{xx}),(w,V_{\mathrm{eff}}),\left(-\frac{dw}{dx},{\overline{M}}_{xx}\right).$$

(16)

where $V_{\mathrm{eff}}$ is the effective shear force

$$V_{\mathrm{eff}}\equiv \frac{d{\overline{M}}_{xx}}{dx}+N_{xx}\frac{dw}{dx}$$

(17)

2.3. Governing Equations in Terms of Displacements

For an isotropic material the one-dimensional stress–strain relation

$${\sigma }_{xx}=E\left(z\right){\epsilon }_{xx}(x,z)$$

(18)

We assume that the beam is graded with two material combination through the beam height according to the relation

$$E\left(z\right)=\left(E_1-E_2\right)V_1\left(z\right)+E_2,V_1\left(z\right)={\left(\frac{1}{2}+\frac{z}{h}\right)}^n$$

(19)

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\left(z\right)$ and $V_2\left(z\right)=1-V_1\left(z\right)$. We assume that Poisson’s ratio $\nu$ is a constant for the FGM material.

The stress resultants can be expressed in terms of the displacements as

$$\begin{array}{cc}\hfill N_{xx}& ={\int }_A{\sigma }_{xx}dA=A_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]-B_{xx}\frac{d^{2}w}{d{x}^2}\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill M_{xx}& ={\int }_{A}z{\sigma }_{xx}dA=B_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]-D_{xx}\frac{d^{2}w}{d{x}^2}\hfill \end{array}$$

(20b)

$$\begin{array}{cc}\hfill P_{xy}& ={\int }_A{\mathcal{M}}_{xy}dA=-A_{xy}\frac{d^{2}w}{d{x}^2}\hfill \end{array}$$

(20c)

where $A_{xx}$, $B_{xx}$, $D_{xx}$, and $A_{xy}$ are the extensional, extensional-bending, bending, and in-plane shear stiffness coefficients. For beams with $E=E\left(x\right)$ (i.e., $n=0$ or $E_1=E_2=E$), we have $A_{xx}=E{A}_0$, $B_{xx}=0$, $D_{xx}=E{I}_0$, and $A_{xy}=G{A}_0{\ell }^2$. For $n\ne 0$ (i.e., FGM beams), we have

$$\begin{array}{cc}\hfill A_{xx}& ={\int }_{A}E\left(z\right)dA=E_2{A}_0\frac{M+n}{1+n}\hfill \\ \hfill B_{xx}& ={\int }_{A}E\left(z\right)zdA=E_2{B}_0\frac{n(M-1)}{2(1+n)(2+n)}\hfill \\ \hfill D_{xx}& ={\int }_{A}E\left(z\right)z^{2}dA=E_2{I}_0\left[\frac{(6+3n+3n^2)M+(8n+3n^2+n^3)}{6+11n+6n^2+n^3}\right]\hfill \\ \hfill A_{xy}& =\frac{{\ell }^2}{2(1+\nu )}{\int }_{A}E\left(z\right)dA={\ell }^2\frac{E_2{A}_0}{2(1+\nu )}\frac{M+n}{1+n}\hfill \\ \hfill A_0& =bh,B_0=b{h}^2,I_0=\frac{b{h}^3}{12},M=\frac{E_1}{E_2}.\hfill \end{array}$$

(21)

Figure 1a contains the variation of the non-dimensional axial stiffness (${\overline{A}}_{xx}=A_{xx}/E_2{A}_0$, $A_0=bh$) and bending stiffness ${\overline{D}}_{xx}=D_{xx}/E_2{I}_0$, $I_0=b{h}^3/12$) as functions of the volume fraction index n for various values of the modulus ratio $M=E_1/E_2\ge 1$ and Figure 1b contains similar plots for the non-dimensional extensional-bending coupling stiffness ${\overline{B}}_{xx}=B_{xx}/E_2{B}_0$, $B_0=bh$). It is clear that both ${\overline{A}}_{xx}$ and ${\overline{D}}_{xx}$ are the maximum at $n=0$ and decrease with increasing values of n. However, ${\overline{B}}_{xx}$ is zero at $n=0$, increases to a maximum at $n=\sqrt{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

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left\{A_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+B_{xx}\left(-\frac{d^{2}w}{d{x}^2}\right)-N_{xx}^T\right\}-f& =0\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill -\frac{d^2}{d{x}^2}\left\{B_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+\left(D_{xx}+A_{xy}\right)\left(-\frac{d^{2}w}{d{x}^2}\right)-M_{xx}^T\right\}& \hfill \\ \hfill -\frac{d}{dx}\left\{A_{xx}\frac{dw}{dx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]-B_{xx}\frac{dw}{dx}\frac{d^{2}w}{d{x}^2}\right\}-q& =0\hfill \end{array}$$

(23)

2.4. Exact Solutions

2.4.1. General Solution

The exact solution to the linearized equations of equilibrium (i.e., static case) under distributed load $q\left(x\right)$ is given by

$$\begin{array}{cc}\hfill u\left(x\right)& =\frac{{\widehat{D}}_{xx}}{{\widehat{D}}_{xx}^{*}}{K}_{1}x+\frac{B_{xx}}{{\widehat{D}}_{xx}^{*}}{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\eta d\xi \hfill \\ \hfill & -\frac{B_{xx}}{{\widehat{D}}_{xx}^{*}}\left(K_2\frac{x^2}{2}+K_{3}x+K_4\right)\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill w\left(x\right)& =\left(\frac{B_{xx}}{{\widehat{D}}_{xx}^{*}}{K}_1\right)\frac{x^2}{2}-\frac{A_{xx}}{{\widehat{D}}_{xx}^{*}}\left(K_2\frac{x^3}{6}+K_3\frac{x^2}{2}+K_{5}x+K_6\right)\hfill \\ \hfill & +\frac{A_{xx}}{{\widehat{D}}_{xx}^{*}}{\int }^x{\int }^{\xi }{\int }^{\eta }{\int }^{\zeta }q\left(\mu \right)d\mu d\zeta d\eta d\xi \hfill \end{array}$$

(25)

where $K_1$ through $K_6$ are constants of integration and $\xi$, $\eta$, $\zeta$, and $\mu$ 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_{\mathrm{eff}}=d{\overline{M}}_{xx}/dx$), and ($dw/dx,{\overline{M}}_{xx}$). The stress resultants $N_{xx}$ and ${\overline{M}}_{xx}$ can be computed using Equations (20a) and (20b).

We can determine the constants of integration for various boundary conditions (pinned: $u=w={\overline{M}}_{xx}=0$; hinged: $N_{xx}=w={\overline{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.

2.4.2. Pinned-Hinged Beams

The exact solution for this case, with FGM and modified couple stress (MCS) effect, is ($\xi =x/L$):

$$\begin{array}{cc}\hfill {\widehat{D}}_{xx}^{*}u\left(\xi \right)& =B_{xx}\frac{q_0{L}^3}{12}\left(-3{\xi }^2+2{\xi }^3\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\widehat{D}}_{xx}^{*}w\left(\xi \right)& =A_{xx}\frac{q_0{L}^4}{24}\left(\xi -2{\xi }^3+{\xi }^3\right),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\overline{M}}_{xx}& =\frac{q_0{L}^2}{2}\left(\xi -{\xi }^2\right),N_{xx}=0,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill N_{xz}& =\frac{d{M}_{xx}}{dx}=\frac{q_{0}L}{2}\left(1-2\xi \right).\hfill \end{array}$$

(29)

We note that the effect of $B_{xx}$ on the mechanical deflection is zero, while the bending moment is independent of $B_{xx}$.

2.4.3. Pinned-Pinned Beams

The solution for the pinned-pinned FGM beam is ($\xi =x/L$)

$$\begin{array}{cc}\hfill {\widehat{D}}_{xx}^{*}u\left(\xi \right)& =B_{xx}\frac{q_0{L}^3}{12}\left(\xi -3{\xi }^2+2{\xi }^3\right),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\widehat{D}}_{xx}^{*}w\left(\xi \right)& =-\frac{B_{xx}^2}{{\widehat{D}}_{xx}}\frac{q_0{L}^4}{24}\left(\xi -{\xi }^2\right)+A_{xx}\frac{q_0{L}^4}{24}\left(\xi -2{\xi }^3+{\xi }^4\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\overline{M}}_{xx}& =\frac{q_0{L}^2}{2}\left(\xi -{\xi }^2\right),N_{xx}=0,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill N_{xz}& =\frac{d{M}_{xx}}{dx}=\frac{q_{0}L}{2}\left(1-2\xi \right).\hfill \end{array}$$

(33)

When $B_{xx}=0$, the solutions for the pinned-hinged beams and pinned-pinned beams coincide.

2.4.4. 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:

$$E_1=30\times {10}^6\mathrm{psi}(210\mathrm{GPa}),E_2=3\times {10}^6\mathrm{psi}(21\mathrm{GPa}),\nu =0.3$$

We shall investigate the parametric effects of the volume fraction index, n, and boundary conditions on the transverse deflections and bending moment.

Figure 2 contain plots of $u\left(x\right)$ vs. $x/L$ and $w\left(x\right)$ vs. $x/L$, while Figure 3 slope $-(dw/dx)\left(x\right)$ vs. $x/L$ and bending moment $M_{xx}\left(x\right)$ vs. $x/L$. The axial displacement $u\left(x\right)$ exists only because of the coupling coefficient $B_{xx}$ (because $f=0$), and the way $B_{xx}$ varies with n (see Figure 1b) is reflected in the variation of $u\left(x\right)$ with n. In particular, $u\left(x\right)$ increases with increasing values of n for $n<5$ but the magnitude of $u\left(x\right)$ 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\left(0.25L\right)$ and $w\left(0.5L\right)$ with the volume fraction index n. It is clear from the plots that the displacement u increases with n for $n\le 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\left(0.5L\right)$ and $-(dw/dx)\left(L\right)$ (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.

2.4.5. Clamped Beams

For beams clamped at both ends and subjected to uniformly distributed transverse load of intensity $q\left(x\right)=q_0$, we have

$$\begin{array}{cc}\hfill {\widehat{D}}_{xx}^{*}u\left(\xi \right)& =B_{xx}\frac{q_0{L}^3}{12}\left(\xi -3{\xi }^2+2{\xi }^3\right),\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill -{\widehat{D}}_{xx}^{*}\frac{dw}{dx}& =A_{xx}\frac{q_0{L}^3}{12}\left(-\xi +3{\xi }^2-2{\xi }^3\right),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\widehat{D}}_{xx}^{*}w\left(\xi \right)& =A_{xx}\frac{q_0{L}^4}{24}\left({\xi }^2-2{\xi }^3+{\xi }^4\right),\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill M_{xx}\left(x\right)& =-\frac{q_0{L}^2}{12}\left(1-6\xi +6{\xi }^2\right),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill N_{xz}=\frac{d{M}_{xx}}{dx}& =\frac{q_{0}L}{2}\left(1-2\xi \right),N_{xx}=0.\hfill \end{array}$$

(38)

The variation of $u\left(x\right)$ for the clamped-clamped beam is the same as that of the pinned-pinned beam. Figure 5 contains plots of $w\left(x\right)$ vs. $x/L$ and $-(dw/dx)\left(x\right)$ vs. $x/L$ while Figure 6 contains the bending moment $M_{xx}\left(x\right)$ 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.

3. First-Order Theory of Beams (TBT)

3.1. 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 ${\gamma }_{xz}=2{\epsilon }_{xz}$ and shear stress ${\sigma }_{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=\frac{5(1+\nu )}{(6+5\nu )}$, which takes the range values $(5/6)\le K_s\le (15/17)$ for $0\le \nu \le 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.

3.2. Displacements and Strains

The TBT is based on the displacement field

$$\mathbf{u}(x,z)=\left[u\left(x\right)+z{\varphi }_x\left(x\right)\right]{\widehat{\mathbf{e}}}_x+w\left(x\right){\widehat{\mathbf{e}}}_z$$

(39)

where ${\varphi }_x$ denotes the rotation (independent of the slope, ${\theta }_x=-dw/dx$) of the cross-sectional 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

$$\begin{array}{cc}\hfill {\epsilon }_{xx}(x,z)& =\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2+z\frac{d{\varphi }_x}{dx}\equiv {\epsilon }_{xx}^{\left(0\right)}+z{\epsilon }_{xx}^{\left(1\right)}\hfill \end{array}$$

(40a)

$$\begin{array}{cc}\hfill {\gamma }_{xz}\left(x\right)={\varphi }_x& +\frac{dw}{dx},{\epsilon }_{xx}^{\left(0\right)}=\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2,{\epsilon }_{xx}^{\left(1\right)}=\frac{d{\varphi }_x}{dx}\hfill \end{array}$$

(40b)

where G the shear modulus [$G\left(z\right)=E\left(z\right)/2(1+\nu )$] and $\nu$ is Poisson’s ratio, which is assumed to be a constant.

3.3. 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 ${\sigma }_{xz}$. The curvature in the TBT is given by

$${\omega }_y=\frac{1}{2}\left(\frac{\partial u_1}{\partial z}-\frac{\partial u_3}{\partial x}\right)=\frac{1}{2}\left({\varphi }_x-\frac{dw}{dx}\right),{\chi }_{xy}=\frac{1}{2}\frac{\partial {\omega }_y}{\partial x}=\frac{1}{4}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)$$

(41)

The stress resultants of the TBT are defined as

$$\begin{array}{cc}\hfill N_{xx}& ={\int }_A{\sigma }_{xx}dA,N_{xz}=K_s{\int }_A{\sigma }_{xz}dA,\hfill \\ \hfill M_{xx}& ={\int }_{A}z{\sigma }_{xx}dA,P_{xy}={\int }_A{\mathcal{M}}_{xy}dA.\hfill \end{array}$$

(42)

Here $K_s$ denotes the shear correction factor. The Euler equations of the TBT are

$$\begin{array}{cc}\hfill -\frac{d{N}_{xx}}{dx}-f& =0\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill -\frac{d{N}_{xz}}{dx}-\frac{1}{2}\frac{d^2{P}_{xy}}{d{x}^2}-\frac{d}{dx}\left(N_{xx}\frac{dw}{dx}\right)-q& =0\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill -\frac{d{M}_{xx}}{dx}-\frac{1}{2}\frac{d{P}_{xy}}{dx}+N_{xz}& =0\hfill \end{array}$$

(45)

The three duality pairs for the TBT are

$$(u,N_{xx}),(w,V_{\mathrm{eff}}),\left({\varphi }_x,{\overline{M}}_{xx}\right).$$

(46)

where the effective shear force and bending moments are

$$\begin{array}{cc}\hfill V_{\mathrm{eff}}& \equiv N_{xz}+N_{xx}\frac{dw}{dx}+\frac{1}{2}\frac{d{P}_{xy}}{dx}\hfill \end{array}$$

(47a)

$$\begin{array}{cc}\hfill {\overline{M}}_{xx}& =M_{xx}+\frac{1}{2}{P}_{xy}\hfill \end{array}$$

(47b)

We note that the effective shear force in the TBT has the modified couple stress term.

3.4. 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

$$\begin{array}{cc}\hfill N_{xx}& ={\int }_A{\sigma }_{xx}dA=A_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+B_{xx}\frac{d{\varphi }_x}{dx}\hfill \end{array}$$

(48a)

$$\begin{array}{cc}\hfill M_{xx}& ={\int }_{A}z{\sigma }_{xx}dA=B_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+D_{xx}\frac{d{\varphi }_x}{dx}\hfill \end{array}$$

(48b)

$$\begin{array}{cc}\hfill P_{xy}& ={\int }_A{\mathcal{M}}_{xy}dA=\frac{1}{2}{A}_{xy}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)\hfill \end{array}$$

(48c)

$$\begin{array}{cc}\hfill N_{xz}& =K_s{\int }_A{\sigma }_{xz}dA=S_{xz}\left({\varphi }_x+\frac{dw}{dx}\right),\hfill \end{array}$$

(48d)

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

$$S_{xz}=\frac{K_s}{2(1+\nu )}{\int }_{A}E\left(z\right)dA$$

(49)

The equations of equilibrium in Equations (43)–(45) now can be expressed in terms of the displacements u, w, and ${\varphi }_x$ with the help of the beam constitutive relations in Equations (48a)–(48c) as

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left\{A_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+B_{xx}\frac{d{\varphi }_x}{dx}\right\}& =f,\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left[S_{xz}\left({\varphi }_x+\frac{dw}{dx}\right)\right]-\frac{1}{4}\frac{d^2}{d{x}^2}\left[A_{xy}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)\right]& \hfill \\ \hfill -\frac{d}{dx}\left\{A_{xx}\frac{dw}{dx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+B_{xx}\frac{dw}{dx}\left(\frac{d{\varphi }_x}{dx}\right)\right\}-q& =0\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left\{B_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+D_{xx}\frac{d{\varphi }_x}{dx}\right\}+S_{xz}\left({\varphi }_x+\frac{dw}{dx}\right)& \hfill \\ \hfill -\frac{1}{4}\frac{d}{dx}\left[A_{xy}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)\right]& =0\hfill \end{array}$$

(52)

3.5. Exact Solutions

3.5.1. 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

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left[A_{xx}\left(\frac{du}{dx}\right)+B_{xx}\frac{d{\varphi }_x}{dx}\right]& =f,\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left[S_{xz}\left({\varphi }_x+\frac{dw}{dx}\right)\right]-q& =0,\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left(B_{xx}\frac{du}{dx}+D_{xx}\frac{d{\varphi }_x}{dx}\right)+S_{xz}\left({\varphi }_x+\frac{dw}{dx}\right)& =0.\hfill \end{array}$$

(55)

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:

$$\frac{d{N}_{xx}}{dx}=0,-\frac{d{N}_{xz}}{dx}-q=0,-\frac{d{M}_{xx}}{dx}+N_{xz}=0$$

(56)

Substituting for $N_{xz}$ from the third equation into the second equation, we obtain

$$\frac{d{N}_{xx}}{dx}=0,-\frac{d^2{M}_{xx}}{d{x}^2}-q=0.$$

(57)

Integrating the above equations,

$$N_{xx}=c_1,\frac{d{M}_{xx}}{dx}=-{\int }^{x}q\left(\xi \right)d\xi +c_2$$

(58)

Integrating the second equation, we obtain

$$M_{xx}=-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +c_{2}x+c_3\equiv F\left(x\right),$$

(59)

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,{\varphi }_x$) using Equations (48a) and (48b); we have

$$\left[\begin{array}{cc}\hfill A_{xx}& \hfill B_{xx}\\ \hfill B_{xx}& \hfill D_{xx}\end{array}\right]\left\{\begin{array}{c}\frac{du}{dx}\\ \frac{d{\varphi }_x}{dx}\end{array}\right\}=\left\{\begin{array}{c}{c}_1+N_{xx}^T\\ F\left(x\right)+M_{xx}^T\end{array}\right\}$$

(60)

Solving for $du/dx$ and $d{\varphi }_x/dx$, we obtain

$$\begin{array}{cc}\hfill \frac{du}{dx}& =\frac{1}{D_{xx}^{*}}\left[D_{xx}{c}_1-B_{xx}F\left(x\right)\right]\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill \frac{d{\varphi }_x}{dx}& =\frac{1}{D_{xx}^{*}}\left[-B_{xx}{c}_1+A_{xx}F\left(x\right)\right]\hfill \end{array}$$

(62)

where

$$D_{xx}^{*}=A_{xx}{D}_{xx}-B_{xx}^2,F\left(x\right)=-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +c_{2}x+c_3$$

(63)

Integrating Equations (61) and (62), we obtain

$$\begin{array}{cc}\hfill u\left(x\right)& =\frac{D_{xx}}{D_{xx}^{*}}{c}_{1}x+\frac{B_{xx}}{D_{xx}^{*}}{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\eta d\xi \hfill \\ \hfill & -\frac{B_{xx}}{D_{xx}^{*}}\left(c_2\frac{x^2}{2}+c_{3}x+c_4\right),\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill {\varphi }_x\left(x\right)& =-\frac{B_{xx}}{D_{xx}^{*}}{c}_{1}x-\frac{A_{xx}}{D_{xx}^{*}}{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\eta d\xi \hfill \\ \hfill & +\frac{A_{xx}}{D_{xx}^{*}}\left(c_2\frac{x^2}{2}+c_{3}x+c_5\right).\hfill \end{array}$$

(65)

From Equations (56) and (58), we arrive at

$$N_{xz}=\frac{d{M}_{xx}}{dx}=-{\int }^{x}q\left(\xi \right)d\xi +c_2$$

(66)

and using Equation (48d) we obtain

$$\begin{array}{cc}\hfill \frac{dw}{dx}& =\frac{1}{S_{xz}}\left(-{\int }^{x}q\left(\xi \right)d\xi +c_2\right)-{\varphi }_x\hfill \\ \hfill & =\frac{1}{S_{xz}}\left(-{\int }^{x}q\left(\xi \right)d\xi +c_2\right)+\frac{B_{xx}}{D_{xx}^{*}}{c}_{1}x\hfill \\ \hfill & -\frac{A_{xx}}{D_{xx}^{*}}\left(-{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\eta d\xi +c_2\frac{x^2}{2}+c_{3}x+c_5\right)\hfill \end{array}$$

(67)

or

$$\begin{array}{cc}\hfill w\left(x\right)& =\frac{1}{S_{xz}}\left(-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +c_{2}x\right)+\frac{B_{xx}}{D_{xx}^{*}}{c}_1\frac{x^2}{2}\hfill \\ \hfill & -\frac{A_{xx}}{D_{xx}^{*}}\left(-{\int }^x{\int }^{\xi }{\int }^{\eta }{\int }^{\zeta }q\left(\mu \right)d\mu d\zeta d\eta d\xi +c_2\frac{x^3}{6}+c_3\frac{x^2}{2}+c_{5}x+c_6\right)\hfill \end{array}$$

(68)

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 (${\varphi }_x,M_{xx}$). We note that, in TBT, ${\varphi }_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 ${\varphi }_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).

3.5.2. Pinned-Pinned Beams

The exact solution of a beam pinned at both ends ($u\left(0\right)=0$, $w\left(0\right)=0$, $M_{xx}\left(0\right)=0$, $u\left(L\right)=0$, $w\left(L\right)=0$, and $M_{xx}\left(L\right)=0$) is

$$\begin{array}{cc}\hfill D_{xx}^{*}u\left(\xi \right)& =B_{xx}\frac{q_0{L}^3}{12}\left(\xi -3{\xi }^2+2{\xi }^3\right),\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill D_{xx}^{*}{\varphi }_x\left(\xi \right)& =\frac{B_{xx}^2}{D_{xx}}\frac{q_0{L}^3}{24}\left(1-2\xi \right)\hfill \\ \hfill & -A_{xx}\frac{q_0{L}^3}{24}\left(1-6\xi +4{\xi }^3\right),\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill D_{xx}^{*}w\left(\xi \right)& =-\frac{B_{xx}^2}{D_{xx}}\frac{q_0{L}^4}{24}\left(\xi -{\xi }^2\right)\hfill \\ \hfill & +A_{xx}\frac{q_0{L}^4}{24}\left(\xi -2{\xi }^3+{\xi }^4\right)+\frac{D_{xx}^{*}}{S_{xz}}\frac{q_0{L}^2}{2}\left(\xi -{\xi }^2\right),\hfill \end{array}$$

(71)

$$\begin{array}{cc}\hfill {\overline{M}}_{xx}& =\frac{q_0{L}^2}{2}\left(\xi -{\xi }^2\right),N_{xz}=\frac{q_{0}L}{2}\left(1-2\xi \right),N_{xx}=0.\hfill \end{array}$$

(72)

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. Numerical results obtained with the classical (CBT) and the shear deformation (TBT) beam theories for the displacements $\overline{u}=u\left(0.25L\right)\times {10}^2$ and $w\left(0.5L\right)$ and slopes ${\theta }_x\left(L\right)=-(dw/dx)\left(L\right)$ and ${\varphi }_x\left(L\right)$ of a pinned–pinned FGM beams under a uniformly distributed load (all results are normalized with the load).

n	$\overline{\mathit{u}}$-CBT	$\overline{\mathit{u}}$-TBT	w-CBT	w-TBT	${\mathit{\theta }}_{\mathit{x}}$	${\mathit{\varphi }}_{\mathit{x}}$
0.0	0.00000	0.00000	0.5208	0.5210	0.01657	0.01657
1.0	0.09973	0.09973	1.0014	1.0016	0.03062	0.03062
2.0	0.20118	0.20118	1.2635	1.2638	0.03722	0.03722
3.0	0.26301	0.26301	1.4261	1.4265	0.04148	0.04148
4.0	0.29297	0.29297	1.5440	1.5445	0.04495	0.04495
5.0	0.30523	0.30523	1.6415	1.6420	0.04806	0.04806
6.0	0.30850	0.30850	1.7286	1.7292	0.05097	0.05097
10.0	0.29367	0.29367	2.0305	2.0312	0.06131	0.06131
20.0	0.24302	0.24302	2.6047	2.6056	0.08080	0.08080

shows the numerical results obtained with the two theories for the data:

$$E_1=30\times {10}^6\mathrm{psi}(210\mathrm{GPa}),E_2=3\times {10}^6\mathrm{psi}(21\mathrm{GPa}),\nu =0.3,K=\frac{5}{6}$$

Transverse deflections obtained with the CBT and TBT for three different length-to-height ratios, $L/h=50$, $L/h=20$, and $L/h=10$ are presented in

Table 2. Numerical results obtained with the classical (CBT) and the Timoshenko (TBT) beam theories for the transverse deflections of a pinned–pinned FGM beams under a uniformly distributed load (the results are normalized with the load, $q_0$); $\overline{w}=w\left(0.5L\right)10$ and $\widehat{w}=w\left(0.5L\right)\times {10}^2$.

	$\mathit{L}/\mathit{h}=50$		$\mathit{L}/\mathit{h}=20$		$\mathit{L}/\mathit{h}=10$
n	w-CBT	w-TBT	$\overline{w}$-CBT	$\overline{w}$-TBT	$\widehat{w}$-CBT	$\widehat{w}$-TBT
0.0	0.06510	0.06517	0.04167	0.04193	0.05208	0.05338
1.0	0.12517	0.12529	0.08011	0.08058	0.10014	0.10250
2.0	0.15793	0.15809	0.10108	0.10173	0.12635	0.12960
3.0	0.17827	0.17847	0.11409	0.11409	0.14261	0.14661
4.0	0.19300	0.19324	0.12352	0.12445	0.15440	0.15905
5.0	0.20518	0.20544	0.13132	0.13236	0.16415	0.16935
6.0	0.21608	0.21636	0.13829	0.13943	0.17286	0.17855
10.0	0.25382	0.25417	0.16244	0.16387	0.20305	0.21020
20.0	0.32559	0.32605	0.20838	0.21020	0.26047	0.26957

. 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.

4. The Third-Order Beam Theory

4.1. 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 ${\sigma }_{xz}=2G{\epsilon }_{xz}$, is either zero (in CBT) or constant (in TBT), although the actual variation of ${\sigma }_{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]).

4.2. Kinematics

The displacement field of the Reddy–Bickford third-order beam theory (RBT) is

$$\mathbf{u}(x,z)=\left[u\left(x\right)+z{\varphi }_x\left(x\right)-\alpha z^3\left({\varphi }_x+\frac{dw}{dx}\right)\right]{\widehat{\mathbf{e}}}_x+w\left(x\right){\widehat{\mathbf{e}}}_z,$$

(73)

where $\alpha =4/3h^2$. The nonzero strain and curvature components are

$$\begin{array}{c}{\epsilon }_{xx}{\epsilon }_{xx}^{\left(0\right)}+z{\epsilon }_{xx}^{\left(1\right)}+z^3{\epsilon }_{xx}^{\left(3\right)},{\gamma }_{xz}={\gamma }_{xz}^{\left(0\right)}+z^2{\gamma }_{xz}^{\left(2\right)},\hfill \end{array}$$

(74a)

$$\begin{array}{cc}\hfill {\chi }_{xy}& ={\chi }_{xy}^{\left(0\right)}+z^2{\chi }_{xy}^{\left(2\right)},{\chi }_{yz}=z{\chi }_{yz}^{\left(1\right)}\hfill \end{array}$$

(74b)

where (omitting the higher-order terms in the thickness strain)

$$\begin{array}{cc}\hfill {\epsilon }_{xx}^{\left(0\right)}& =\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2,{\epsilon }_{xx}^{\left(1\right)}=\frac{d{\varphi }_x}{dx},{\epsilon }_{xx}^{\left(3\right)}=-\alpha \left(\frac{d{\varphi }_x}{dx}+\frac{d^{2}w}{d{x}^2}\right),\hfill \\ \hfill & {\gamma }_{xz}^{\left(0\right)}={\varphi }_x+\frac{dw}{dx},{\gamma }_{xz}^{\left(2\right)}=-\beta \left({\varphi }_x+\frac{dw}{dx}\right),\beta =3\alpha =\frac{4}{h^2}\hfill \\ \hfill & {\chi }_{xy}^{\left(0\right)}=\frac{1}{4}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right),{\chi }_{xy}^{\left(2\right)}=-\frac{1}{4}\beta \left(\frac{d{\varphi }_x}{dx}+\frac{d^{2}w}{d{x}^2}\right),\hfill \\ \hfill & {\chi }_{yz}^{\left(1\right)}=-\frac{1}{2}\beta \left({\varphi }_x+\frac{dw}{dx}\right)\hfill \end{array}$$

(74c)

4.3. 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:

$$\begin{array}{cc}\hfill (N_{xx},M_{xx},P_{xx})& ={\int }_A(1,z,z^3){\sigma }_{xx}dA,(N_{xz},P_{xz})={\int }_A(1,z^2){\sigma }_{xz}dA,\hfill \end{array}$$

(75a)

$$\begin{array}{cc}\hfill (P_{xy},R_{xy})& ={\int }_A(1,z^2){\mathcal{M}}_{xy}dA,Q_{yz}={\int }_{A}z{\mathcal{M}}_{yz}dA,\hfill \end{array}$$

(75b)

$$\begin{array}{cc}\hfill {\overline{M}}_{xx}=M_{xx}& -\alpha P_{xx},{\overline{N}}_{xz}=N_{xz}-\beta P_{xz},\beta =3\alpha =\frac{4}{h^2}.\hfill \end{array}$$

(75c)

The equations of equilibrium of the RBT are:

$$\begin{array}{cc}\hfill -\frac{d{N}_{xx}}{dx}& =f\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left(N_{xx}\frac{dw}{dx}\right)-\frac{d{\overline{N}}_{xz}}{dx}-\alpha \frac{d^2{P}_{xx}}{d{x}^2}& \hfill \\ \hfill -\frac{1}{2}\frac{d^2{P}_{xy}}{d{x}^2}-\frac{1}{2}\beta \frac{d^2{R}_{xy}}{d{x}^2}+\beta \frac{d{Q}_{yz}}{dx}& =q\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill -\frac{d{\overline{M}}_{xx}}{dx}+{\overline{N}}_{xz}-\frac{1}{2}\frac{d{P}_{xy}}{dx}+\frac{1}{2}\beta \frac{d{R}_{xy}}{dx}-\beta Q_{yz}& =0\hfill \end{array}$$

(78)

The duality pairs (the first element of the pair denotes a generalized displacement while the second element denotes a generalized force):

$$\left(u,N_{xx}\right),\left(w,V_{\mathrm{eff}}\right),\left(-\frac{dw}{dx},M_{\mathrm{eff}}\right),\left({\varphi }_x,{\tilde{M}}_{xx}\right),$$

(79)

where

$$\begin{array}{cc}\hfill V_{\mathrm{eff}}& ={\overline{N}}_{xz}+\alpha \frac{d{P}_{xx}}{dx}+N_{xx}\frac{dw}{dx}+\frac{1}{2}\frac{d{P}_{xy}}{dx}+\frac{1}{2}\beta \frac{d{R}_{xy}}{dx}-\beta Q_{yz},\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill M_{\mathrm{eff}}& =\frac{1}{2}{P}_{xy}+\frac{1}{2}\beta R_{xy},{\tilde{M}}_{xx}={\overline{M}}_{xx}+\frac{1}{2}{P}_{xy}-\beta \frac{1}{2}{R}_{xy}.\hfill \end{array}$$

(81)

Thus, there are four boundary conditions at each boundary point. Requiring $\partial w/\partial x$ as well as ${\varphi }_x$ to vanish at a support necessarily implies that the shear force, when shear stress is computed using the constitutive relation ${\sigma }_{xz}=G{\gamma }_{xz}$, is zero. However, the effective shear force $V_{\mathrm{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 that ${\overline{M}}_{xx}=M_{xx}$ (implying that $P_{xx}=0$ there).

4.4. 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 ${\epsilon }_{zz}=0$. Therefore, we can use one-dimensional constitutive relations. In particular, the one-dimensional constitutive relations are

$$\begin{array}{cc}\hfill {\sigma }_{xx}& =E{\epsilon }_{xx},{\sigma }_{xz}=G{\gamma }_{xz},\hfill \end{array}$$

(82)

$$\begin{array}{cc}\hfill {\mathcal{M}}_{xy}& =2{\ell }^{2}G{\chi }_{xy},{\mathcal{M}}_{yz}=2{\ell }^{2}G{\chi }_{yz}.\hfill \end{array}$$

(83)

Substituting the constitutive relations from Equations (82) and (83) into the definition of the stress resultants in Equations (75a)–(75c), we obtain

$$\begin{array}{cc}\hfill \left\{\begin{array}{c}{N}_{xx}\\ M_{xx}\\ P_{xx}\end{array}\right\}& =\left[\begin{array}{ccc}{A}_{xx}& B_{xx}& E_{xx}\\ B_{xx}& D_{xx}& F_{xx}\\ E_{xx}& F_{xx}& H_{xx}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{xx}^{\left(3\right)}\end{array}\right\},\hfill \end{array}$$

(84a)

$$\begin{array}{cc}\hfill \left\{\begin{array}{c}{N}_{xz}\\ P_{xz}\\ P_{xy}\\ R_{xy}\\ Q_{yz}\end{array}\right\}& =\left[\begin{array}{ccccc}{A}_{xz}& D_{xz}& 0& 0& 0\\ D_{xz}& F_{xz}& 0& 0& 0\\ 0& 0& A_{xy}& D_{xy}& 0\\ 0& 0& D_{xy}& F_{xy}& 0\\ 0& 0& 0& 0& D_{xy}\end{array}\right]\left\{\begin{array}{c}{\gamma }_{xz}^{\left(0\right)}\\ {\gamma }_{xz}^{\left(2\right)}\\ 2{\chi }_{xy}^{\left(0\right)}\\ 2{\chi }_{xy}^{\left(2\right)}\\ 2{\chi }_{yz}^{\left(1\right)}\end{array}\right\}.\hfill \end{array}$$

(84b)

where (see Equation (21))

$$\begin{array}{c}\hfill (A_{xx},B_{xx},D_{xx},E_{xx},F_{xx},H_{xx}) ={\int }_A(1,z,z^2,z^3,z^4,z^6)E\left(z\right)dA,\end{array}$$

(85a)

$$\begin{array}{cc}\hfill (A_{xz},D_{xz},F_{xz})& =\frac{1}{2(1+\nu )}{\int }_A(1,z^2,z^4)E\left(z\right)dA,\hfill \end{array}$$

(85b)

$$\begin{array}{cc}\hfill (A_{xy},D_{xy},F_{xy})& =\frac{{\ell }^2}{2(1+\nu )}{\int }_A(1,z^2,z^4)E\left(z\right)dA,\hfill \end{array}$$

(85c)

4.5. Beam Stiffness Coefficients for FGM Beams

For the FGM beams, the integrals in Equations (85a)–(85c) can be evaluated as:

$$\begin{array}{cc}\hfill A_{xx}& =E_{2}bh\frac{M+n}{1+n},B_{xx}=E_2\frac{b{h}^2}{2}\frac{n(M-1)}{(1+n)(2+n)},\hfill \\ \hfill D_{xx}& =E_2\frac{b{h}^3}{12}\left[\frac{(6+3n+3n^2)M+(8n+3n^2+n^3)}{(1+n)(2+n)(3+n)}\right],\hfill \\ \hfill E_{xx}& =E_2\frac{b{h}^4}{8}(M-1)\left[\frac{n(8+3n+n^2)}{(1+n)(2+n)(3+n)(4+n)}\right],\hfill \\ \hfill F_{xx}& =E_2\frac{b{h}^5}{80}f\left(n\right),H_{xx}=E_2\frac{b{h}^7}{448}g\left(n\right),\hfill \\ \hfill A_{xy}& =E_2{\ell }^2\frac{bh}{2(1+\nu )}\frac{M+n}{1+n},F_{xy}=\frac{E_2{\ell }^{2}b{h}^5}{120(1+\nu )}f\left(n\right),\hfill \\ \hfill D_{xy}& =E_2{\ell }^2\frac{b{h}^3}{24(1+\nu )}\left[\frac{(6+3n+3n^2)M+(8n+3n^2+n^3)}{(1+n)(2+n)(3+n)}\right],\hfill \end{array}$$

(86a)

where

$$\begin{array}{cc}\hfill f\left(n\right)& =\frac{f_{1}M+n{f}_2}{f_3},g\left(n\right)=\frac{g_{1}M+g_2}{g_3},M=\frac{E_1}{E_2},\hfill \\ \hfill f_1& =(24+18n+23n^2+6n^3+n^4),\hfill \\ \hfill f_2& =(184+110n+55n^2+10n^3+n^4),\hfill \\ \hfill f_3& =(1+n)(2+n)(3+n)(4+n)(5+n),\hfill \\ \hfill g_1& =(720+660n+964n^2+405n^3+115n^4+15n^5+n^6),\hfill \\ \hfill g_2& =(720+1764n+1624n^2+735n^3+175n^4+21n^5+n^6),\hfill \\ \hfill g_3& =(1+n)(2+n)(3+n)(4+n)(5+n)(6+n)(7+n).\hfill \end{array}$$

(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.

4.6. 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,{\varphi }_x$). We obtain

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left\{A_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+{\overline{B}}_{xx}\frac{d{\varphi }_x}{dx}-\alpha E_{xx}\frac{d^{2}w}{d{x}^2}\right\}& =f,\hfill \end{array}$$

(87)

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left[\frac{dw}{dx}\left\{A_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+{\overline{B}}_{xx}\frac{d{\varphi }_x}{dx}-\alpha E_{xx}\frac{d^{2}w}{d{x}^2}\right\}\right]& \hfill \\ \hfill -\frac{1}{4}\frac{d^2}{d{x}^2}\left[A_{xy}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)-\beta D_{xy}\left(\frac{d{\varphi }_x}{dx}+\frac{d^{2}w}{d{x}^2}\right)\right]& \hfill \\ \hfill -\frac{1}{4}\beta \frac{d^2}{d{x}^2}\left[D_{xy}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)-\beta F_{xy}\left(\frac{d{\varphi }_x}{dx}+\frac{d^{2}w}{d{x}^2}\right)\right]& \hfill \\ \hfill -\alpha \frac{d^2}{d{x}^2}\left\{E_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+{\overline{F}}_{xx}\frac{d{\varphi }_x}{dx}-\alpha H_{xx}\frac{d^{2}w}{d{x}^2}\right\}& \hfill \\ \hfill -{\widehat{A}}_{xz}\frac{d}{dx}\left({\varphi }_x+\frac{dw}{dx}\right)-{\beta }^2{D}_{xy}\frac{d}{dx}\left({\varphi }_x+\frac{dw}{dx}\right)& =q,\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill -\frac{d}{dx}\left\{{\overline{B}}_{xx}\left[\frac{du}{dx}+\frac{1}{2}{\left(\frac{dw}{dx}\right)}^2\right]+{\widehat{D}}_{xx}\frac{d{\varphi }_x}{dx}-\alpha {\overline{F}}_{xx}\frac{d^{2}w}{d{x}^2}\right\}& \hfill \\ \hfill -\frac{1}{4}\frac{d}{dx}\left[A_{xy}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)-\beta D_{xy}\left(\frac{d{\varphi }_x}{dx}+\frac{d^{2}w}{d{x}^2}\right)\right]& \hfill \\ \hfill +\frac{1}{4}\beta \frac{d}{x}\left[D_{xy}\left(\frac{d{\varphi }_x}{dx}-\frac{d^{2}w}{d{x}^2}\right)-\beta F_{xy}\left(\frac{d{\varphi }_x}{dx}+\frac{d^{2}w}{d{x}^2}\right)\right]& \hfill \\ \hfill +{\widehat{A}}_{xz}\left({\varphi }_x+\frac{dw}{dx}\right)+{\beta }^2{D}_{xy}\left({\varphi }_x+\frac{dw}{dx}\right)& =0,\hfill \end{array}$$

(89)

where

$$\begin{array}{cc}\hfill & {\overline{B}}_{xx}=B_{xx}-\alpha E_{xx},{\overline{D}}_{xx}=D_{xx}-\alpha F_{xx},{\overline{F}}_{xx}=F_{xx}-\alpha H_{xx},\hfill \\ \hfill & {\widehat{D}}_{xx}={\overline{D}}_{xx}-\alpha {\overline{F}}_{xx},{\overline{A}}_{xz}=A_{xz}-\beta D_{xz},{\overline{D}}_{xz}=D_{xz}-\beta F_{xz},\hfill \\ & {\widehat{A}}_{xz}={\overline{A}}_{xz}-\beta {\overline{D}}_{xz},\alpha =\frac{4}{3h^2},\beta =\frac{4}{h^2}=3\alpha .\hfill \end{array}$$

(90)

4.7. 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:

$$\begin{array}{cc}\hfill -\frac{d{N}_{xx}}{dx}& =0\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill -\frac{d{\overline{N}}_{xz}}{dx}-\alpha \frac{d^2{P}_{xx}}{d{x}^2}& =q\hfill \end{array}$$

(92)

$$\begin{array}{cc}\hfill -\frac{d{\overline{M}}_{xx}}{dx}+{\overline{N}}_{xz}& =0\hfill \end{array}$$

(93)

Integrating the equations with respect to x, we obtain

$$\begin{array}{cc}\hfill N_{xx}& =K_1\hfill \end{array}$$

(94)

$$\begin{array}{cc}\hfill {\overline{N}}_{xz}& +\alpha \frac{d{P}_{xx}}{dx}=-{\int }^{x}q\left(\xi \right)d\xi +K_2\hfill \end{array}$$

(95)

$$\begin{array}{cc}\hfill M_{xx}& =-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +K_{2}x+K_3\equiv F\left(x\right).\hfill \end{array}$$

(96)

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{\varphi }_x/dx$ in terms of $d^{2}w/d{x}^2$, we obtain

$$\begin{array}{cc}\hfill \frac{du}{dx}& =\frac{P_1}{{\overline{D}}_{xx}^{*}}+\frac{{\overline{D}}_{xx}}{{\overline{D}}_{xx}^{*}}{K}_1-\frac{{\overline{B}}_{xx}}{{\overline{D}}_{xx}^{*}}F\left(x\right)+\frac{J_1}{{\overline{D}}_{xx}^{*}}\frac{d^{2}w}{d{x}^2},\hfill \end{array}$$

(97a)

$$\begin{array}{cc}\hfill \frac{d{\varphi }_x}{dx}& =\frac{P_2}{{\overline{D}}_{xx}^{*}}-\frac{B_{xx}}{{\overline{D}}_{xx}^{*}}{K}_1+\frac{A_{xx}}{{\overline{D}}_{xx}^{*}}F\left(x\right)+\frac{J_2}{{\overline{D}}_{xx}^{*}}\frac{d^{2}w}{d{x}^2},\hfill \end{array}$$

(97b)

where (see Equations (87) and (88))

$$\begin{array}{cc}\hfill P_1& ={\overline{D}}_{xx}{M}_T^{\left(0\right)}-{\overline{B}}_{xx}{M}_T^{\left(1\right)},P_2=A_{xx}{M}_T^{\left(1\right)}-B_{xx}{M}_T^{\left(0\right)},\hfill \\ \hfill J_1& =\alpha \left({\overline{D}}_{xx}{E}_{xx}-{\overline{B}}_{xx}{F}_{xx}\right),J_2=\alpha \left(A_{xx}{F}_{xx}-B_{xx}{E}_{xx}\right),\hfill \\ \hfill {\overline{D}}_{xx}^{*}& =A_{xx}{\overline{D}}_{xx}-B_{xx}{\overline{B}}_{xx},D_{xx}^{*}=A_{xx}{D}_{xx}-B_{xx}{B}_{xx}.\hfill \end{array}$$

(97c)

Integrating the two equations in (97a) and (97b), we obtain

$$\begin{array}{cc}\hfill {\overline{D}}_{xx}^{*}u\left(x\right)& =-{\overline{B}}_{xx}\left(-{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\xi d\eta +K_2\frac{x^2}{2}+K_{3}x+K_4\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+J_1\frac{dw}{dx}+\left(P_1+{\overline{D}}_{xx}{K}_1\right)x,\hfill \end{array}$$

(98)

$$\begin{array}{cc}\hfill {\overline{D}}_{xx}^{*}{\varphi }_x\left(x\right)& =A_{xx}\left(-{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\xi d\eta +K_2\frac{x^2}{2}+K_{3}x+K_5\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+J_2\frac{dw}{dx}+\left(P_2-B_{xx}{K}_1\right)x,\hfill \end{array}$$

(99)

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):

$$\begin{array}{c}0=-{\widehat{A}}_{xz}\left({\varphi }_x+\frac{dw}{dx}\right)+\left(-{\int }^{x}q\left(\xi \right)d\xi +K_2\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\alpha \frac{d}{dx}\left(E_{xx}\frac{du}{dx}+{\overline{F}}_{xx}\frac{d{\varphi }_x}{dx}-\alpha H_{xx}\frac{d^{2}w}{d{x}^2}\right)\hfill \\ \hspace{1em}\hspace{1em}=-\frac{{\widehat{A}}_{xz}}{{\overline{D}}_{xx}^{*}}[A_{xx}\left(-{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\xi d\eta +K_2\frac{x^2}{2}+K_{3}x+K_5\right)\hfill \\ +J_2\frac{dw}{dx}+\left(P_2-B_{xx}{K}_1\right)x+{\overline{D}}_{xx}^{*}\frac{dw}{dx}]\hfill \\ -\frac{1}{{\overline{D}}_{xx}^{*}}\frac{d}{dx}\left[P_{3}F\left(x\right)+\alpha \left(E_{xx}{J}_1+{\overline{F}}_{xx}{J}_2-\alpha H_{xx}{\overline{D}}_{xx}^{*}\right)\frac{d^{2}w}{d{x}^2}\right]\hfill \\ +\left(-{\int }^{x}q\left(\xi \right)d\xi +K_2\right).\hfill \end{array}$$

(100)

where

$$P_3=\alpha \left(A_{xx}{\overline{F}}_{xx}-{\overline{B}}_{xx}{E}_{xx}\right),{\widehat{D}}_{xx}^{*}={\overline{D}}_{xx}^{*}-P_3=A_{xx}{\widehat{D}}_{xx}-{\overline{B}}_{xx}{\overline{B}}_{xx}.$$

(101)

Integrating Equation (100) once and collecting the like terms together, we obtain

$$\begin{array}{cc}\hfill 0& =-\frac{{\widehat{A}}_{xz}}{{\overline{D}}_{xx}^{*}}\left({\overline{D}}_{xx}^{*}+J_2\right)w+\frac{\alpha }{{\overline{D}}_{xx}^{*}}\left(\alpha {\overline{D}}_{xx}^{*}{H}_{xx}-E_{xx}{J}_1-{\overline{F}}_{xx}{J}_2\right)\frac{d^{2}w}{d{x}^2}\hfill \\ \hfill & -\frac{{\widehat{A}}_{xz}}{{\overline{D}}_{xx}^{*}}\left[A_{xx}\right(-{\int }^x{\int }^{\xi }{\int }^{\eta }{\int }^{\zeta }q\left(\mu \right)d\mu d\zeta d\eta d\xi \hfill \\ \hfill & +K_2\frac{x^3}{6}+K_3\frac{x^2}{2}+K_{5}x+K_6)+\left(P_2-B_{xx}{K}_1\right)\frac{x^2}{2}]\hfill \\ \hfill & +\frac{{\widehat{D}}_{xx}^{*}}{{\overline{D}}_{xx}^{*}}\left(-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +K_{2}x\right)-\frac{P_3}{{\overline{D}}_{xx}^{*}}{K}_3\hfill \\ \hfill & =-c_{1}w+c_2\frac{d^{2}w}{d{x}^2}+g\left(x\right),\hfill \end{array}$$

(102)

where

$$\begin{array}{cc}\hfill c_1& =\frac{{\widehat{A}}_{xz}}{{\overline{D}}_{xx}^{*}}({\overline{D}}_{xx}^{*}+J_2)=\frac{{\widehat{A}}_{xz}{D}_{xx}^{*}}{{\overline{D}}_{xx}^{*}},\hfill \\ \hfill c_2& =\frac{\alpha }{{\overline{D}}_{xx}^{*}}\left(\alpha {\overline{D}}_{xx}^{*}{H}_{xx}-E_{xx}{J}_1-{\overline{F}}_{xx}{J}_2\right),\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill g\left(x\right)& =-\frac{{\widehat{A}}_{xz}}{{\overline{D}}_{xx}^{*}}\left[A_{xx}\right(-{\int }^x{\int }^{\xi }{\int }^{\eta }{\int }^{\zeta }q\left(\mu \right)d\mu d\zeta d\eta d\xi +K_2\frac{x^3}{6}+K_3\frac{x^2}{2}\hfill \\ \hfill & +K_{5}x+K_6)+\left(P_2-B_{xx}{K}_1\right)\frac{x^2}{2}]\hfill \\ \hfill & +\frac{{\widehat{D}}_{xx}^{*}}{{\overline{D}}_{xx}^{*}}\left(-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +K_{2}x\right)-\frac{P_3}{{\overline{D}}_{xx}^{*}}{K}_3.\hfill \end{array}$$

(104)

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\left(x\right)=K_{7}cosh\mu x+K_{8}sinh\mu x,\mu =\sqrt{\frac{c_1}{c_2}}.$$

(105)

The total solution $w\left(x\right)$ is obtained by adding the particular solution, $w_p\left(x\right)$ due to $g\left(x\right)$: $w\left(x\right)=w_h\left(x\right)+w_p\left(x\right)$. 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 ${\varphi }_x$ (or their dual variables, $P_{xx}$ and ${\overline{M}}_{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

$$\begin{array}{cc}\hfill 0& =-\frac{{\widehat{A}}_{xz}{D}_{xx}^{*}}{{\overline{D}}_{xx}^{*}}w-\frac{{\widehat{A}}_{xz}}{{\overline{D}}_{xx}^{*}}\left[A_{xx}\right(-{\int }^x{\int }^{\xi }{\int }^{\eta }{\int }^{\zeta }q\left(\mu \right)d\mu d\zeta d\eta d\xi \hfill \\ \hfill & +K_2\frac{x^3}{6}+K_3\frac{x^2}{2}+K_{5}x+K_6)\hfill \\ \hfill & +\left(A_{xx}{M}_T^{\left(1\right)}-B_{xx}{M}_T^{\left(0\right)}-B_{xx}{K}_1\right)\frac{x^2}{2}]\hfill \\ \hfill & +\frac{{\widehat{D}}_{xx}^{*}}{{\overline{D}}_{xx}^{*}}\left(-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +K_{2}x\right)-\frac{P_3}{{\overline{D}}_{xx}^{*}}{K}_3,\hfill \end{array}$$

(106)

The solution to these equations is (after a lengthy algebra)

$$\begin{array}{cc}\hfill u\left(x\right)& =\frac{1}{A_{xx}}\left(N_{xx}^T+K_1\right)x,\hfill \end{array}$$

(107)

$$\begin{array}{cc}\hfill {\varphi }_x\left(x\right)& =\frac{1}{D_{xx}}\left(-{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\xi d\eta +K_2\frac{x^2}{2}+K_{3}x\right)\hfill \\ \hfill & +\frac{17}{56}\frac{1}{A_{xz}}\left(-{\int }^{x}q\left(\xi \right)d\xi +K_2\right)+\frac{1}{D_{xx}}{M}_{xx}^{T}x+\frac{1}{D_{xx}}{K}_5,\hfill \end{array}$$

(108)

$$\begin{array}{cc}\hfill w\left(x\right)& =-\frac{1}{D_{xx}}(-{\int }^x{\int }^{\xi }{\int }^{\eta }{\int }^{\zeta }q\left(\mu \right)d\mu d\zeta d\eta d\xi +K_2\frac{x^3}{6}+K_3\frac{x^2}{2}\hfill \\ \hfill & +K_{5}x)-\frac{1}{D_{xx}}{M}_{xx}^T\frac{x^2}{2}-\frac{1}{D_{xx}}{K}_6-\frac{2}{7}\frac{1}{A_{xz}}{K}_3\hfill \\ \hfill & +\frac{17}{14}\frac{1}{A_{xz}}\left(-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +K_{2}x\right),\hfill \end{array}$$

(109)

The corresponding solutions obtained using the first-order shear deformation theory (TBT) of beams are:

$$\begin{array}{cc}\hfill u\left(x\right)& =\frac{1}{A_{xx}}\left(N_{xx}^T+K_1\right)x,\hfill \end{array}$$

(110)

$$\begin{array}{cc}\hfill {\varphi }_x\left(x\right)& =\frac{1}{D_{xx}}\left(-{\int }^x{\int }^{\xi }{\int }^{\eta }q\left(\zeta \right)d\zeta d\eta d\xi +K_2\frac{x^2}{2}+K_{3}x\right)\hfill \\ \hfill & +\frac{1}{D_{xx}}{M}_{xx}^{T}x+\frac{1}{D_{xx}}{K}_5,\hfill \end{array}$$

(111)

$$\begin{array}{cc}\hfill w\left(x\right)& =-\frac{1}{D_{xx}}(-{\int }^x{\int }^{\xi }{\int }^{\eta }{\int }^{\zeta }q\left(\mu \right)d\mu d\zeta d\eta d\xi +K_2\frac{x^3}{6}\hfill \\ \hfill & +K_3\frac{x^2}{2}+K_{5}x)-\frac{1}{D_{xx}}{M}_{xx}^T\frac{x^2}{2}-\frac{1}{D_{xx}}{K}_6\hfill \\ \hfill & +\frac{1}{A_{xz}}\left(-{\int }^x{\int }^{\xi }q\left(\eta \right)d\eta d\xi +K_{2}x\right),\hfill \end{array}$$

(112)

A comparison of Equation (112) with Equation (109) shows that the shear correction factor predicted by the “simplified" RBT is $K_s=14/17$, which is slightly smaller than the value suggested for rectangular cross-section beams, which is $K_s=5/6$. Of course, the TBT solution is different from the RBT solution. Particularly, the expression for ${\varphi }_x$ in the RBT contains a term due to transverse shear coefficient.

As an example, we present here the exact solutions using the simplified RBT of a functionally graded beam with both ends pinned and subjected to a uniformly distributed load of magnitude $q_0$. We take the origin of the x-coordinate at the left end of the beam (i.e., $0\le x\le L$). For this case, the boundary conditions at both ends, $x=0$ and $x=L$, are:

$$u=w=0,{\overline{M}}_{xx}=0\to M_{xx}=0.$$

(113)

Use of these boundary conditions yield

$$\begin{array}{cc}\hfill u\left(x\right)& =\frac{B_{xx}}{D_{xx}^{*}}\frac{q_0{L}^3}{12}\frac{x}{L}\left(1-3\frac{x}{L}+2\frac{x^2}{L^2}\right)+\frac{J_1{P}_3}{{\widehat{A}}_{xz}{D}_{xx}^{*}}\frac{q_{0}L}{2}\left(1-\frac{x}{L}\right),\hfill \\ \hfill {\varphi }_x\left(x\right)& =-\frac{A_{xx}}{D_{xx}^{*}}\frac{q_0{L}^3}{12}\frac{x}{L}\left(1-3\frac{x}{L}+2\frac{x^2}{L^2}\right)-\frac{1}{D_{xx}}\frac{q_0{L}^3}{24}\left(1-2\frac{x}{L}\right)\hfill \end{array}$$

(114)

$$\begin{array}{cc}\hfill & -\frac{M_{xx}^{T}L}{2D_{xx}}\left(1-2\frac{x}{L}\right)+\frac{J_2{P}_3{q}_0}{2D_{xx}^{*}{\widehat{A}}_{xz}}\left(2J_2+J_1\frac{B_{xx}}{D_{xx}}\right)\left(1-2\frac{x}{L}\right),\hfill \\ \hfill w\left(x\right)& =\frac{A_{xx}}{D_{xx}^{*}}\frac{q_0{L}^4}{24}\frac{x}{L}\left(1-2\frac{x^2}{L^2}+\frac{x^3}{L^3}\right)-\frac{B_{xx}^2}{D_{xx}{D}_{xx}^{*}}\frac{q_0{L}^4}{24}\frac{x}{L}\left(1-\frac{x}{L}\right)\hfill \end{array}$$

(115)

$$\begin{array}{cc}\hfill & +\frac{M_T^{\left(1\right)}{L}^2}{2D_{xx}}\frac{x}{L}\left(2-\frac{x}{L}\right)+\frac{P_3{q}_0{L}^2}{2{\widehat{A}}_{xz}{D}_{xx}^{*}}\left(2-J_1\frac{B_{xx}}{D_{xx}}\right)\frac{x}{L}\left(1-\frac{x}{L}\right).\hfill \end{array}$$

(116)

5. Summary

In this paper three different beam theories, namely, the classical, first-order, and third-order 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 ($\overline{w}=w{\widehat{D}}_{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.