Analytical Solutions for Functionally Graded Beams under Arbitrary Distributed Loads via the Symplectic Approach

Abstract

A novel symplectic approach is employed in the analysis of homogenous and functionally graded beams subjected to arbitrary tractions on the lateral surfaces. Two models of functionally graded beams are heterogeneous in the sense that the material properties are exponential functions of the length and thickness, respectively. Within the symplectic framework, the method of separation of variables alone with the eigenfunction expansion technique is adopted to obtain the exact analysis of displacement and stress fields. The complete solutions include homogenous solutions with coefficients to be determined by two end boundary conditions and particular solutions satisfying the lateral boundary conditions. Two examples are presented for functionally graded beams to demonstrate the effects of material inhomogeneity. The efficiency and accuracy of the symplectic analysis are shown through numerical results.

1. Introduction

Functionally graded materials (FGMs) as a new kind of composite materials with material properties vary continuously within the solid and have gained great interests for decades [1]. There are plenty of literatures on the mechanical behaviors of the functionally graded structures in the last two decades [2–8]. Of particular interest is the special case that the material properties vary in the thickness direction only. For instance, Ding et al. [9] and Huang et al. [10] presented the analytical solutions for FGM beams and functionally graded magneto-electro-elastic beams by using the stress function approach. Nie et al. [11] employed the displacement function approach to analyze the functionally graded beams with arbitrary material properties along the thickness. As far as I know, exact analytical solutions for FGM beams with material constants varying along the length direction have not been reported in the literatures. Recently, a few researchers have investigated the static and dynamic analysis of bidirection FGM structures with material properties varying along two directions and obtained numerical or semianalytical solutions [12–14]. It should be pointed that ample skills and experiences are often required on the exact analysis of FGM structures to seek the trial shape function.

To overcome the above difficulty in the conventional analytical methods, the symplectic method developed by Zhong [15] for elasticity problems of homogeneous materials can be extended to obtain the analytical solutions of FGM beams. The approach based on Hamiltonian system has shown great superiority in revealing the structure of solutions and their physical essence as well as predicting the accurate local behavior which is usually covered up by the Saint-Venant principle in the traditional elasticity analysis [16]. Furthermore, the symplectic approach as a rational and a unified manner has been widely applied to various problems in many various branches of applied mechanics [17–23]. In the numerical calculation, Zhao and Chen [24] reconstructed the symplectic expansion formula in the state space formalism to avoid matrix singularity and achieved the stability of numerical calculation. For the symplectic analysis of inhomogeneous materials, Zhao et al. [25–28] systematically developed a symplectic framework to plane problems of functionally graded elastic materials, functionally graded piezoelectric materials, and functionally graded magneto-electro-elastic materials whose material properties change continually along the length direction. More recently, Zhao et al. [29] further presented the analytical solutions for bidirection functionally graded cantilever beams subjected to different loads at the free end.

It is noted that, the symplectic analyses for the functionally graded beams mentioned above were considered for the lateral traction-free boundary conditions. The objective of this work is to develop a unified mathematical treatment for functionally graded beams subjected to arbitrary distributed loads on the lateral surfaces. The present work involves the well-structured analytical procedure of elasticity beams and functionally graded beams with material properties varying exponentially in the length and thickness direction, respectively. It is interesting that, though the symplectic conjugate orthogonality for functionally graded beams does not hold any more, the eigenfunction expansion technique can still be used to give accurate predictions, as clarified in Example 1 through the convergence analysis as well as comparison with corresponding results of homogeneous beam.

Different from existing studies, this study emphasizes on the particular solutions of generally supported beams subjected to arbitrary form tractions along length direction. By incorporating the particular solution with eigensolutions that satisfy the state equation and free lateral boundary conditions, exact elasticity solutions for plane beams with any combinations of end conditions are determined in a systematic manner. Therefore, symplectic approach is effective in predicting accurate local stress distributions which are in practice very important in evaluating the possible failure or fracture of materials and structures which usually originate from corners, discontinuities, and other singularities. Also, the study allows us to evaluate and compare the influence of material inhomogeneities on the mechanical behavior of structures along the length and thickness direction, respectively. The exponential form may also be used to approximate any material variation through proper subdivision (for which Taylor expansion may be used, indicating that the exponential form is closely related to the polynomial one). Furthermore, with the development of material fabrication technology, it may be possible to control the spatial distributions of material properties according to an optimized design.

2. Basic Formulations

We consider the plane stress problem of a rectangular domain (Figure 1). For isotropic material, the constitutive relations are

$\begin{array}{c}{\sigma }_x=\frac{E}{\mathrm{1}-{\nu }^{\mathrm{2}}}\left(\frac{\partial u_x}{\partial x}+\nu \frac{\partial u_z}{\partial z}\right),\\ {\sigma }_z=\frac{E}{\mathrm{1}-{\nu }^{\mathrm{2}}}\left(\nu \frac{\partial u_x}{\partial x}+\frac{\partial u_z}{\partial z}\right),\\ {\tau }_{xz}=\frac{E}{\mathrm{2}(\mathrm{1}+\nu )}\left(\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}\right),\end{array}$

(1)

where u _x and u _z are displacements in the x- and z-directions, respectively, and ν is the constant Poisson's ratio. For Young's modulus E, we will consider three models of materials as that in Section 3.

Figure 1:

The plane beam problem.

Absence of body force, the governing equations of equilibrium are

$\frac{\partial {\sigma }_x}{\partial x}+\frac{\partial {\tau }_{xz}}{\partial z}=\mathrm{0},\hspace{1em}\hspace{1em}\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\sigma }_z}{\partial z}=\mathrm{0}.$

(2)

The state equation can be obtained directly from (1) and (2) as

$\frac{\partial }{\partial z}\left\{\begin{array}{c}{u}_z\\ u_x\\ {\sigma }_z\\ {\tau }_{xz}\end{array}\right\}=\left[\begin{array}{cccc}\mathrm{0}& -\nu \frac{\partial }{\partial x}& \frac{\mathrm{1}-{\nu }^{\mathrm{2}}}{E}& \mathrm{0}\\ -\frac{\partial }{\partial x}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{2}(\mathrm{1}+\nu )}{E}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -\frac{\partial }{\partial x}\\ \mathrm{0}& -E\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}& -\nu \frac{\partial }{\partial x}& \mathrm{0}\end{array}\right]\left\{\begin{array}{c}{u}_z\\ u_x\\ {\sigma }_z\\ {\tau }_{xz}\end{array}\right\}$

(3a)

or

$\dot{\mathbf{v}}=\mathbf{H}\mathbf{v},$

(3b)

in which $\left(\dot{\hspace{0.17em}}\right)=\partial /\partial z$ and H is the operator matrix. For homogeneous materials, v = [u _z , u _x , σ _z , τ _xz ]^T is the state vector and the operate matrix H is the Hamiltonian matrix [16]. For functionally graded materials with Young's modulus, which varies exponentially along length or thickness direction, the state vector is reformed as $\mathbf{v}={[u_z,u_x,{\widehat{\sigma }}_z,{\widehat{\tau }}_{xz}]}^{\text{T}}$ , in which ${\widehat{\sigma }}_z$ and ${\widehat{\tau }}_{xz}$ are new forms of stress variables σ _z and τ _xz , respectively. And the operator matrix for inhomogeneous beams exhibits different properties compared with the Hamiltonian matrix [25].

It should be mentioned that the previous literatures for homogenous materials and functionally graded materials in the symplectic framework are emphasized on the analytical solutions for plane beams with the free boundary conditions on the lateral surfaces [25–29]. The present work attempts to obtain the complete stress distributions of functionally graded beams subjected to arbitrary tractions on the lateral surfaces. In the next context, we call the functionally graded beams whose material constants vary along the length direction as axial FGM beams, while those vary along the thickness direction as transverse FGM beams.

3. Symplectic Approach

In this section, a brief review of the symplectic approach for elasticity problems and functionally graded materials problems is given firstly. Then, the particular solutions are presented for generally supported beams subjected to arbitrary form tractions on the lateral surfaces. Thus, the exact analytical solutions are obtained by superimposing the particular solution and all eigensolutions of particular eigenvalues (zero or − α) and general eigenvalues.

3.1. Solutions for Homogenous Materials

First, the symplectic analysis of homogeneous beams subjected to arbitrary lateral tractions is presented. To obtain the eigensolutions of zero eigenvalues and general eigensolutions, we consider the following homogenous boundary conditions on the lateral surfaces as

$x=\pm h:{\sigma }_x=E\frac{\partial u_x}{\partial x}+\nu {\sigma }_z=\mathrm{0},\hspace{1em}\hspace{1em}{\tau }_{xz}=\mathrm{0}.$

(4)

The method of separation of variables is possible for the so-called state equation (3a) and (3b)

$\mathbf{v}(x,z)=\xi \left(z\right)\mathbf{\Phi }\left(x\right)={\text{e}}^{\mu z}{\left[w\left(x\right),u\left(x\right),\sigma \left(x\right),\tau \left(x\right)\right]}^{\text{T}},$

(5)

where μ is the eigenvalue of the operator matrix H and Φ(x) is the corresponding eigenvector.

Zero eigenvalue is multiple and its eigensolutions actually correspond to the classical solutions of Saint-Venant problem. The six eigenvectors of the zero eigenvalues are

$\begin{array}{c}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}={\left[\mathrm{1,0},\mathrm{0,0}\right]}^{\text{T}},\hspace{1em}\hspace{1em}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}={\left[\mathrm{0},-\nu x,E,\mathrm{0}\right]}^{\text{T}}\\ {\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}={\left[\mathrm{0,1},\mathrm{0,0}\right]}^{\text{T}},\hspace{1em}\hspace{1em}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}={\left[-x,\mathrm{0,0},\mathrm{0}\right]}^{\text{T}}\\ {\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}={\left[\mathrm{0},\frac{\nu }{\mathrm{2}}{x}^{\mathrm{2}}-\left(\frac{\mathrm{2}}{\mathrm{5}}+\frac{\nu }{\mathrm{2}}\right)h^{\mathrm{2}},-Ex,\mathrm{0}\right]}^{\text{T}},\\ {\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{3}\right)}={\left[-\frac{\mathrm{6}+\mathrm{5}\nu }{\mathrm{10}}{h}^{\mathrm{2}}x+\frac{\mathrm{2}+\nu }{\mathrm{6}}{x}^{\mathrm{3}},\mathrm{0,0},\frac{E}{\mathrm{2}}\left(x^{\mathrm{2}}-h^{\mathrm{2}}\right)\right]}^{\text{T}}.\end{array}$

(6)

The corresponding eigensolutions can be constructed as

$\begin{array}{c}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}={\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)},\hspace{1em}\hspace{1em}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}={\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}\\ {\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}={\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)},\hspace{1em}\hspace{1em}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}={\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}\\ {\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}={\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)},\\ {\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{3}\right)}={\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{3}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{6}}{z}^{\mathrm{3}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)},\end{array}$

(7)

in which the first subscript 0 indicates eigenvalue zero, the second subscripts 1 and 2 indicate the first and second Jordan chain, the superscripts (0) and (i) (i = 1, 2, 3) indicate, respectively, the fundamental eigenvectors and the ith Jordan normal form eigenvector.

The particular solutions for external forces loaded on the lateral surfaces can be determined from the Jordan normal form solution. We take the plane beam subjected to arbitrary external normal force along z-axis as example, such as uniform loads, linear distribution loads, and cosinusoidal loads applied on the upper and lower surfaces (seen in Figure 2).

Figure 2:

The plane beam with transverse loads at lateral surfaces.

The conditions at the upper and lower surfaces are

$\begin{array}{c}x=-h:({\sigma }_x=)E\frac{\partial u}{\partial x}+\nu \sigma =q_{\mathrm{1}}\left(z\right),\hspace{1em}\hspace{1em}\tau =\mathrm{0},\\ x=h:({\sigma }_x=)E\frac{\partial u}{\partial x}+\nu \sigma =q_{\mathrm{2}}\left(z\right),\hspace{1em}\hspace{1em}\tau =\mathrm{0}.\end{array}$

(8)

A particular solution of $\mathbf{H}\stackrel{\sim }{\mathbf{\Phi }}=k_{\mathrm{1}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{(\mathrm{3})}$ can be solved as follows:

$\tilde{\mathbf{\Phi }}=\left\{\begin{array}{c}\tilde{w}=\mathrm{0}\\ \tilde{u}=-\frac{k_{\mathrm{1}}(\mathrm{1}+\mathrm{2}\nu )}{\mathrm{24}}{x}^{\mathrm{4}}+\frac{k_{\mathrm{1}}(\mathrm{5}+\mathrm{6}\nu )h^{\mathrm{2}}}{\mathrm{20}}{x}^{\mathrm{2}}+\frac{\mathrm{1}-{\nu }^{\mathrm{2}}}{E_{\mathrm{0}}}cx\\ \tilde{\sigma }=k_{\mathrm{1}}{E}_{\mathrm{0}}\left(\frac{\mathrm{1}}{\mathrm{3}}{x}^{\mathrm{3}}-\frac{\mathrm{3}}{\mathrm{5}}{h}^{\mathrm{2}}x\right)+c\nu \\ \tilde{\tau }=\mathrm{0}\end{array}\right\},$

(9)

in which k₁ = 3[q₂(z) − q₁(z)]/(2Eh³) and c = [q₁(z) + q₂(z)]/2.

The particular solution of the plane problem can be constructed as

$\tilde{\mathbf{v}}=\tilde{\mathbf{\Phi }}+k_{\mathrm{1}}\left[z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{3}\right)}+\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+\frac{\mathrm{1}}{\mathrm{6}}{z}^{\mathrm{3}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{24}}{z}^{\mathrm{4}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}\right].$

(10)

It should be pointed that a particular solution can also be deduced for arbitrary distributed shear traction along z-axis. But for concentrated force and other types of discontinuous loads, particular solutions can not be obtained directly in the present symplectic framework. Maybe further endeavors will be made to seek numerical solutions using numerical technique [30].

For the general eigenvalue μ ≠ 0, the homogeneous boundary conditions on lateral surfaces lead to the following equations:

$\mathrm{2}\mu h+\sin \left(\mathrm{2}\mu h\right)=\mathrm{0},\hspace{1em}\hspace{1em}\mathrm{2}\mu h-\sin \left(\mathrm{2}\mu h\right)=\mathrm{0}.$

(11)

It can be seen that, − μ _i must be a root if μ _i is a root which is in accordance with the characteristics of Hamiltonian operator matrix [16]. Equation (11) may be solved by numerical methods, such as Newton's method in the monograph of Yao et al. [16], which is able to provide fast convergent numerical solutions.

It should be noted that, μ _i with the negative sign of real part is classified as α-set eigenvalues, while − μ _i is β-set eigenvalues. The eigensolutions of the two classes eigenvalues are investigated in detail which are signed as v _i and v_{− i} (i = 1, 2, 3, …) in Yao et al. [16] and Zhao and Chen [24].

The exact analytical solution of a loaded beam is formed by zero eigensolutions, general eigensolutions, and the particular solution mentioned above. The complete solution can be described as

$\mathbf{v}=\tilde{\mathbf{v}}+m_{\mathrm{1}}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}+m_{\mathrm{2}}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}+m_{\mathrm{3}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}+m_{\mathrm{4}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+m_{\mathrm{5}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+m_{\mathrm{6}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{3}\right)}+{\displaystyle \sum }_{i=\mathrm{1}}^N\left(a_i{\mathbf{v}}_i+b_i{\mathbf{v}}_{-i}\right).$

(12)

It should be mentioned that N is a truncated number, which should be large enough to ascertain the accuracy of the symplectic expansion.

By substituting the symplectic expansion series of (12) into the Hamiltonian mixed energy variational principle, algebraic linear equations can be constructed and all unknown constants m _i , a _i , and b _i can be solved. The solving procedure of expanding coefficients for homogeneous beams is given in detail in Leung and Zheng [17].

According to [24], numerical instability appears when N is larger than a certain number. This phenomenon is actually associated with the eigensolution ${\mathbf{v}}_{-i}={\text{e}}^{-{\mu }_{i}L}{\mathbf{\Phi }}_{-i}$ in (8) because ${\text{e}}^{-{\mu }_{i}L}$ becomes very large when i is large (Re[μ _i ] < 0). In order to overcome the numerical difficulty, (12) is rewritten as

$\mathbf{v}=\tilde{\mathbf{v}}+m_{\mathrm{1}}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}+m_{\mathrm{2}}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}+m_{\mathrm{3}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}+m_{\mathrm{4}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+m_{\mathrm{5}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+m_{\mathrm{6}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{3}\right)}+{\displaystyle \sum }_{i=\mathrm{1}}^N\left[a_i{\mathbf{v}}_i+\left(b_i{\text{e}}^{-{\mu }_{i}L}\right)\left({\text{e}}^{-{\mu }_i\left(x-L\right)}{\mathbf{\Phi }}_{-i}\right)\right]=\tilde{\mathbf{v}}+m_{\mathrm{1}}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}+m_{\mathrm{2}}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}+m_{\mathrm{3}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}+m_{\mathrm{4}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+m_{\mathrm{5}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+m_{\mathrm{6}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{3}\right)}+{\displaystyle \sum }_{i=\mathrm{1}}^N\left(a_i{\mathbf{v}}_i+b_i^{\prime }{\mathbf{v}}_{-i}^{\prime }\right).$

(13)

Then no big number will appear in the new expansion formula and the numerical stability can be achieved.

3.2. Solutions for Axial FGM Beams

In this section, the axial FGM beams which are functionally graded in the length direction are investigated. The Young modulus is the exponential function of the longitudinal direction in the form of E = E₀e^αz while Poisson's ratio keeps unaltered. Here, E₀ is a constant and α is the gradient index of the material.

The altered stresses ${\widehat{\sigma }}_{ij}={\sigma }_{ij}{\text{e}}^{-\alpha z}\hspace{0.17em}\hspace{0.17em}(i,j=x,z)$ are introduced to ensure the governing differential equations with constant coefficients. Different from the symplectic analysis of homogenous material, the operator matrix H exhibits different characteristics when compared with the Hamiltonian matrix. The so-called shift-Hamiltonian system in Chen and Zhao [25, 26] still allowed the symplectic analysis to be used. For eigen-problem analysis, the particular eigenvalues of the axial FGM are zero and − α. The eigenvectors of the particular eigenvalues are

$\begin{array}{c}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}={\left[\mathrm{1,0},\mathrm{0,0}\right]}^{\text{T}},\hspace{1em}\hspace{1em}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}={\left[\mathrm{0,1},\mathrm{0,0}\right]}^{\text{T}},\\ {\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}={\left[-x,\mathrm{0},\mathrm{0},\mathrm{0}\right]}^{\text{T}}\\ {\mathbf{\Phi }}_{-\alpha ,\mathrm{1}}^{\left(\mathrm{0}\right)}={\left[-\frac{\mathrm{1}}{\alpha }\cosh \left(\lambda x\right),-\frac{\lambda }{{\alpha }^{\mathrm{2}}}\sinh \left(\lambda x\right),E_{\mathrm{0}}\cosh \left(\lambda x\right),\mathrm{0}\right]}^{\text{T}}\\ {\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{\left(\mathrm{0}\right)}={\left[\frac{\mathrm{1}}{\lambda \alpha }\sinh \left(\lambda x\right),\frac{\mathrm{1}}{{\alpha }^{\mathrm{2}}}\cosh \left(\lambda x\right),-\frac{E_{\mathrm{0}}}{\lambda }\sinh \left(\lambda x\right),\mathrm{0}\right]}^{\text{T}}\\ {\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{\left(\mathrm{1}\right)}=\left\{\begin{array}{c}\frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}x\cosh \left(\lambda x\right)+\frac{\nu -\mathrm{1}}{{\lambda }^{\mathrm{3}}}\sinh \left(\lambda x\right)-\frac{C}{\lambda }\sinh \left(\lambda x\right)\\ -\frac{\mathrm{2}}{\alpha {\lambda }^{\mathrm{2}}}\cosh \left(\lambda x\right)+\frac{\mathrm{1}}{\alpha \lambda }x\sinh \left(\lambda x\right)+\frac{\mathrm{2}(\mathrm{1}+\nu )}{\alpha {\lambda }^{\mathrm{2}}}\cosh \left(\lambda h\right)-\frac{C}{\alpha }\cosh \left(\lambda x\right)\\ -\frac{E_{\mathrm{0}}\alpha }{{\lambda }^{\mathrm{2}}}x\cosh \left(\lambda x\right)+\frac{E_{\mathrm{0}}}{\alpha \lambda \nu }\sinh \left(\lambda x\right)+\frac{E_{\mathrm{0}}\alpha C}{\lambda }\sinh \left(\lambda x\right)\\ \frac{E_{\mathrm{0}}}{{\lambda }^{\mathrm{2}}}[\cosh \left(\lambda x\right)-\cosh \left(\lambda h\right)]\end{array}\right\},\end{array}$

(14)

where $\lambda =\alpha \sqrt{\nu }$ and C = h²[3sinh ⁡ (λh)/(λh)³ − 3cosh ⁡ (λh)/(λh)² + sinh ⁡ (λh)/(λh)]/[cosh ⁡ (λh) − sin(λh)/].

In the expressions above, the first subscripts 0 and − α indicate the corresponding eigenvalues, the second subscripts 1 and 2 indicate the first and second Jordan chain, the superscripts (0) and (1) indicate, respectively, the fundamental eigenvectors and the first-order Jordan normal form eigenvector.

Also, the origin solutions associated with the above eigenvectors are

$\begin{array}{c}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}={\mathbf{M}}_t{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)},\hspace{1em}\hspace{1em}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}={\mathbf{M}}_t{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)},\\ {\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}={\mathbf{M}}_t\left({\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}\right)\\ {\mathbf{v}}_{-\alpha ,\mathrm{1}}^{\prime \left(\mathrm{0}\right)}={\mathbf{M}}_t{\text{e}}^{-\alpha z}{\mathbf{\Phi }}_{-\alpha ,\mathrm{1}}^{\left(\mathrm{0}\right)},\hspace{1em}\hspace{1em}{\mathbf{v}}_{-\alpha ,\mathrm{2}}^{\prime \left(\mathrm{0}\right)}={\mathbf{M}}_t{\text{e}}^{-\alpha z}{\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{\left(\mathrm{0}\right)},\\ {\mathbf{v}}_{-\alpha ,\mathrm{2}}^{\prime \left(\mathrm{1}\right)}={\mathbf{M}}_t{\text{e}}^{-\alpha z}\left({\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{\left(\mathrm{1}\right)}+z{\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{\left(\mathrm{0}\right)}\right),\end{array}$

(15)

in which M _t = diag ⁡ [1, 1, e^αz, e^αz] is the transform matrix between eigensolutions and the origin solutions of the plane problem.

To obtain the eigensolutions of the general eigenvalues (μ ≠ 0, − α), we can derive the following characteristic polynomial from the eigenequation HΦ(x) = μΦ(x):

${\eta }^{\mathrm{4}}+\left(\mathrm{2}\alpha \mu +\mathrm{2}{\mu }^{\mathrm{2}}-\nu {\alpha }^{\mathrm{2}}\right){\eta }^{\mathrm{2}}+{\mu }^{\mathrm{2}}{\left(\mu +\alpha \right)}^{\mathrm{2}}=\mathrm{0}.$

(16)

In the analysis above, we assume Φ(x) = e^ηxV as the solution of eigenequation and V is the undermined constant vector. The eigen-root η can be obtained as

$\begin{array}{c}{\eta }_{\mathrm{1}}=-{\eta }_{\mathrm{3}}=\frac{\mathrm{1}}{\mathrm{2}}\alpha \sqrt{\nu }+\frac{\mathrm{1}}{\mathrm{2}}\sqrt{-\mathrm{4}{\mu }^{\mathrm{2}}-\mathrm{4}\alpha \mu +\nu {\alpha }^{\mathrm{2}}},\\ {\eta }_{\mathrm{2}}=-{\eta }_{\mathrm{4}}=\frac{\mathrm{1}}{\mathrm{2}}\alpha \sqrt{\nu }-\frac{\mathrm{1}}{\mathrm{2}}\sqrt{-\mathrm{4}{\mu }^{\mathrm{2}}-\mathrm{4}\alpha \mu +\nu {\alpha }^{\mathrm{2}}}.\end{array}$

(17)

The general solution is then

$\mathbf{\Phi }\left(x\right)=\left[\begin{array}{cccc}{A}_{\mathrm{1}}& B_{\mathrm{1}}& C_{\mathrm{1}}& D_{\mathrm{1}}\\ B_{\mathrm{2}}& A_{\mathrm{2}}& D_{\mathrm{2}}& C_{\mathrm{2}}\\ A_{\mathrm{3}}& B_{\mathrm{3}}& C_{\mathrm{3}}& D_{\mathrm{3}}\\ B_{\mathrm{4}}& A_{\mathrm{4}}& D_{\mathrm{4}}& C_{\mathrm{4}}\end{array}\right]\left\{\begin{array}{c}\cosh \left({\eta }_{\mathrm{1}}x\right)\\ \sinh \left({\eta }_{\mathrm{1}}x\right)\\ \cosh \left({\eta }_{\mathrm{2}}x\right)\\ \sinh \left({\eta }_{\mathrm{2}}x\right)\end{array}\right\}.$

(18)

The general solutions can be divided into the symmetric part (coefficients with the denotations A _i and C _i ) and the antisymmetric part (coefficients with the denotations B _i and D _i ). For the symmetric solution, the following algebraic equation for the eigenvalue μ can be deduced:

${\eta }_{\mathrm{1}}h\tanh \left({\eta }_{\mathrm{1}}h\right)-{\eta }_{\mathrm{2}}h\tanh \left({\eta }_{\mathrm{2}}h\right)=\mathrm{0}.$

(19)

For the antisymmetric solution, there is

${\eta }_{\mathrm{2}}h\tanh \left({\eta }_{\mathrm{1}}h\right)-{\eta }_{\mathrm{1}}h\tanh \left({\eta }_{\mathrm{2}}h\right)=\mathrm{0}.$

(20)

The transcendent equations above are in terms of μ and their roots are solved using the Newton-Rphson method in the present work. Also, the relationships between the constants A _i and C _i or B _i and D _i can be obtained. Hence, we get the eigenvectors Φ _i ^s and Φ _i ^a . The origin solutions of the plane problems are

${\mathbf{v}}_i^{\prime s}{\mathbf{M}}_t{\text{e}}^{{\mu }_{n}z}{\mathbf{\Phi }}_i^s,\hspace{1em}\hspace{1em}{\mathbf{v}}_i^{\prime a}={\mathbf{M}}_t{\text{e}}^{{\mu }_{n}z}{\mathbf{\Phi }}_i^a.$

(21)

Next, we investigate the particular solutions of the axial functionally graded beams subjected to lateral loads (seen in Figure 2). The conditions at the upper and lower surfaces are assumed as

$\begin{array}{c}x=-h:E_{\mathrm{0}}\frac{\partial \tilde{u}}{\partial x}+\nu \tilde{\sigma }=q_{\mathrm{1}}\left(z\right){\text{e}}^{-\alpha z},\hspace{1em}\hspace{1em}\tilde{\tau }=\mathrm{0}\\ x=h:E_{\mathrm{0}}\frac{\partial \tilde{u}}{\partial x}+\nu \tilde{\sigma }=q_{\mathrm{2}}\left(z\right){\text{e}}^{-\alpha z},\hspace{1em}\hspace{1em}\tilde{\tau }=\mathrm{0}.\end{array}$

(22)

The solution of $\mathbf{H}\stackrel{\sim }{\mathbf{\Phi }}=-\alpha \stackrel{\sim }{\mathbf{\Phi }}+k_{\mathrm{2}}{\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{(\mathrm{1})}$ can be solved as follows:

$\left(\tilde{\mathbf{\Phi }}=\right)\left\{\begin{array}{c}\tilde{w}=k_{\mathrm{2}}\left\{\frac{\alpha }{\mathrm{2}{\lambda }^{\mathrm{3}}}{x}^{\mathrm{2}}\sinh \left(\lambda x\right)-\left(\frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}+C\right)\frac{\alpha }{{\lambda }^{\mathrm{2}}}\left[x\cosh \left(\lambda x\right)+\frac{\mathrm{1}+\nu }{\lambda }\sinh \left(\lambda x\right)\right]\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(-\frac{\alpha }{{\lambda }^{\mathrm{4}}}\left(\mathrm{1}-{\nu }^{\mathrm{2}}\right)\cosh \left(\lambda h\right)x\phantom{\frac{\alpha }{\mathrm{2}{\lambda }^{\mathrm{3}}}}\right\}-\frac{\mathrm{1}-{\nu }^{\mathrm{2}}}{E_{\mathrm{0}}}\frac{D}{\alpha \nu }\\ \tilde{u}=k_{\mathrm{2}}\left\{\frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}{x}^{\mathrm{2}}\cosh \left(\lambda x\right)-\left(\frac{\mathrm{3}+\nu }{{\lambda }^{\mathrm{2}}}+C\right)\frac{\mathrm{1}}{\lambda }x\sinh \left(\lambda x\right)\right)+\mathrm{2}\left(\mathrm{1}+\nu \right)\frac{\mathrm{1}}{{\lambda }^{\mathrm{3}}}h\sinh \left(\lambda h\right)\\ \left(\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\left[\left(\mathrm{5}-\mathrm{3}\nu \right)\left(\mathrm{1}+\nu \right)\frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}+\mathrm{2}C\left(\mathrm{1}+\nu \right)\right]\frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}\cosh \left(\lambda h\right)\right\}\\ \tilde{\sigma }=\frac{k_{\mathrm{2}}{E}_{\mathrm{0}}}{\nu }\left[-\frac{x^{\mathrm{2}}}{\mathrm{2}\lambda }\sinh \left(\lambda x\right)+\left(\frac{\mathrm{2}+\nu }{{\lambda }^{\mathrm{2}}}+C\right)\left[x\cosh \left(\lambda x\right)+\frac{\sinh \left(\lambda x\right)}{\lambda }\right]+\frac{x}{{\lambda }^{\mathrm{2}}}\cosh \left(\lambda h\right)\right]+\frac{D}{\nu }\\ \tilde{\tau }=k_{\mathrm{2}}\frac{E_{\mathrm{0}}\alpha }{{\lambda }^{\mathrm{3}}}\left\{x\sinh \left(\lambda x\right)-h\sinh \left(\lambda h\right)-\left(\frac{\mathrm{2}}{\lambda }+C\lambda \right)\left[\cosh \left(\lambda x\right)-\cosh \left(\lambda h\right)\right]\right\}\end{array}\right\},$

(23)

in which k₂ = [q₁(z) − q₂(z)]λ³e^−αz/2E₀[sinh ⁡ (λh) − λhcosh ⁡ (λh)] and D = ((q₁(z) + q₂(z))/2)e^−αz.

The particular solution of the plane problem can be constructed as

$\tilde{\mathbf{v}}={\mathbf{M}}_t{\text{e}}^{-\alpha z}\left[\tilde{\mathbf{\Phi }}+k_{\mathrm{2}}\left(z{\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{-\alpha ,\mathrm{2}}^{\left(\mathrm{0}\right)}\right)\right].$

(24)

The exact analytical solution of the FGM beam is a combination of all eigensolutions and particular solution which satisfies the loading conditions on the lateral surfaces in the form of

${\mathbf{v}}^{\prime }=\tilde{\mathbf{v}}+m_{\mathrm{1}}{\mathbf{v}}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}+m_{\mathrm{2}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}+m_{\mathrm{3}}{\mathbf{v}}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+m_{\mathrm{4}}{\mathbf{v}}_{-\alpha ,\mathrm{1}}^{\prime \left(\mathrm{0}\right)}+m_{\mathrm{5}}{\mathbf{v}}_{-\alpha ,\mathrm{2}}^{\prime \left(\mathrm{0}\right)}+m_{\mathrm{6}}{\mathbf{v}}_{-\alpha ,\mathrm{2}}^{\prime \left(\mathrm{1}\right)}+{\displaystyle \sum }_{i=\mathrm{1}}^{\infty }(a_i^s{\mathbf{v}}_i^{\prime s}+a_i^a{\mathbf{v}}_i^{\prime a}+b_i^s{\mathbf{v}}_{-i}^{\prime s}+b_i^a{\mathbf{v}}_{-i}^{\prime a}).$

(25)

Similar to procedure of homogeneous materials, we present a novel form of the β-class eigensolutions (v_{− i}^′s and v_{− i}^′a) in the symplectic expansion series to overcome the numerical instability.

3.3. Solutions for Transverse FGM Beams

For comparison, the symplectic analysis is also included for transverse functionally graded materials whose Young's modulus varies exponentially along the thickness direction in the form of E = E₀e^βx and Poisson's ratio keeps unaltered. Also, E₀ is a constant and β is the gradient index of the material.

In the symplectic framework, we introduce the new stress variables ${\widehat{\sigma }}_{ij}={\sigma }_{ij}{\text{e}}^{-\beta x}\mathrm{\hspace{0.17em}\hspace{0.17em}}(i,j=x,z)$ into the state vector as $\mathit{\nu }={\left[u_z,u_x,{\widehat{\sigma }}_z,{\widehat{\tau }}_{xz}\right]}^{\text{T}}$ . Because of the material inhomogeneity, the operator matrix H performs different properties from Hamiltonian matrix. But it can degenerate to the conventional Hamiltonian matrix when β = 0.

Similar to the analysis of homogenous materials, explicit expressions of the eigenvectors and eigensolutions for zero eigenvalue are obtainable. Consider

$\begin{array}{c}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}={\left[\mathrm{1,0},\mathrm{0,0}\right]}^{\text{T}},\hspace{1em}\hspace{1em}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}={\left[\mathrm{0},-\nu x,E_{\mathrm{0}},\mathrm{0}\right]}^{\text{T}}\\ {\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}={\left[\mathrm{0,1},\mathrm{0,0}\right]}^{\text{T}},\hspace{1em}\hspace{1em}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}={\left[-x,\mathrm{0,0},\mathrm{0}\right]}^{\text{T}},\\ {\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}={\left[0,\frac{1}{2}v{x}^2,-E_{0}x,\mathrm{0}\right]}^{\text{T}}\\ {\mathbf{\Phi }}_{\mathrm{0}}^{\left(\mathrm{3}\right)}=\left\{\begin{array}{c}-\frac{\mathrm{1}}{\mathrm{6}}\nu x^{\mathrm{3}}+\mathrm{2}(\mathrm{1}+\nu )\left[\frac{\mathrm{1}}{\mathrm{2}{\beta }^{\mathrm{2}}}{x}^{\mathrm{2}}-\frac{\mathrm{1}}{{\beta }^{\mathrm{2}}}x-\frac{c_{\mathrm{1}}}{E_{\mathrm{0}}\beta }{e}^{-\beta x}-\frac{A_{\mathrm{0}}}{\beta }x\right]+\frac{\mathrm{1}}{\mathrm{2}}\nu A_{\mathrm{0}}{x}^{\mathrm{2}}+c_{\mathrm{2}}\\ \mathrm{0}\\ \mathrm{0}\\ \frac{E_{\mathrm{0}}}{\beta }x-\frac{E_{\mathrm{0}}}{{\beta }^{\mathrm{2}}}+c_{\mathrm{1}}{e}^{-\beta x}-\frac{A_{\mathrm{0}}{E}_{\mathrm{0}}}{\beta }\end{array}\right\}.\end{array}$

(26)

It should be noted that Φ₀⁽³⁾ has been determined by means of Jordan chain Φ₀⁽³⁾ = HΦ₀⁽²⁾, in which Φ₀⁽²⁾ was reconstructed as

${\mathbf{\Phi }}_{\mathrm{0}}^{\left(\mathrm{2}\right)}={\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+A_{\mathrm{0}}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)},$

(27)

where A₀ = [βhcosh ⁡ (βh)/sinh ⁡ (βh) − 1]/β, c₁ = E₀h/βsinh ⁡ (βh), and c₂ = 2(1 + ν)h/β²sinh ⁡ (βh).

The eigensolutions corresponding to the eigenvectors above can be constructed as

$\begin{array}{c}{\mathbf{v}}_{\mathrm{0,1}}^{\prime \left(\mathrm{0}\right)}={\mathbf{M}}_t^{\prime }{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)},\hspace{1em}\hspace{1em}{\mathbf{v}}_{\mathrm{0,1}}^{\prime \left(\mathrm{1}\right)}={\mathbf{M}}_t^{\prime }\left({\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}\right)\\ {\mathbf{v}}_{\mathrm{0,2}}^{\prime \left(\mathrm{0}\right)}={\mathbf{M}}_t^{\prime }{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)},\hspace{1em}\hspace{1em}{\mathbf{v}}_{\mathrm{0,2}}^{\prime \left(\mathrm{1}\right)}={\mathbf{M}}_t^{\prime }\left({\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}\right),\\ {\mathbf{v}}_{\mathrm{0,2}}^{\prime \left(\mathrm{2}\right)}={\mathbf{M}}_t^{\prime }\left({\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}\right)\\ {\mathbf{v}}_{\mathrm{0}}^{\prime \left(\mathrm{3}\right)}={\mathbf{M}}_t^{\prime }\left[{\mathbf{\Phi }}_{\mathrm{0}}^{\left(\mathrm{3}\right)}+\left(z{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{6}}{z}^{\mathrm{3}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}\right)\right)\hspace{1em}\left(+A_{\mathrm{0}}\left(z{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}\right)\right],\end{array}$

(28)

in which M _t ′ = diag ⁡ [1, 1, e^βx, e^βx] is the transform matrix between the eigenvectors and original solutions.

For the general eigenvalues μ ≠ 0, the characteristic polynomial of the eigenequation can be

${\eta }^{\mathrm{4}}+\mathrm{2}\beta {\eta }^{\mathrm{3}}+\left({\beta }^{\mathrm{2}}+\mathrm{2}{\mu }^{\mathrm{2}}\right){\eta }^{\mathrm{2}}+\mathrm{2}\beta {\mu }^{\mathrm{2}}\eta -{\beta }^{\mathrm{2}}\nu {\mu }^{\mathrm{2}}+{\mu }^{\mathrm{4}}=\mathrm{0}.$

(29)

The eigen-roots are

$\begin{array}{c}{\eta }_{\mathrm{1}}=-\frac{\mathrm{1}}{\mathrm{2}}\beta +\frac{\mathrm{1}}{\mathrm{2}}\sqrt{{\beta }^{\mathrm{2}}-\mathrm{4}{\mu }^{\mathrm{2}}-\mathrm{4}\mu \beta \sqrt{\nu }},\\ {\eta }_{\mathrm{2}}=-\frac{\mathrm{1}}{\mathrm{2}}\beta -\frac{\mathrm{1}}{\mathrm{2}}\sqrt{{\beta }^{\mathrm{2}}-\mathrm{4}{\mu }^{\mathrm{2}}-\mathrm{4}\mu \beta \sqrt{\nu }},\\ {\eta }_{\mathrm{3}}=-\frac{\mathrm{1}}{\mathrm{2}}\beta +\frac{\mathrm{1}}{\mathrm{2}}\sqrt{{\beta }^{\mathrm{2}}-\mathrm{4}{\mu }^{\mathrm{2}}+\mathrm{4}\mu \beta \sqrt{\nu }},\\ {\eta }_{\mathrm{4}}=-\frac{\mathrm{1}}{\mathrm{2}}\beta -\frac{\mathrm{1}}{\mathrm{2}}\sqrt{{\beta }^{\mathrm{2}}-\mathrm{4}{\mu }^{\mathrm{2}}+\mathrm{4}\mu \beta \sqrt{\nu }},\end{array}$

(30)

which results in the general solution as follows:

$\mathbf{\Phi }\left(x\right)=\left[\begin{array}{cccc}{A}_{\mathrm{1}}& A_{\mathrm{2}}& A_{\mathrm{3}}& A_{\mathrm{4}}\\ B_{\mathrm{1}}& B_{\mathrm{2}}& B_{\mathrm{3}}& B_{\mathrm{4}}\\ C_{\mathrm{1}}& C_{\mathrm{2}}& C_{\mathrm{3}}& C_{\mathrm{4}}\\ D_{\mathrm{1}}& D_{\mathrm{2}}& D_{\mathrm{3}}& D_{\mathrm{4}}\end{array}\right]\left\{\begin{array}{c}{\text{e}}^{{\eta }_{\mathrm{1}}x}\\ {\text{e}}^{{\eta }_{\mathrm{2}}x}\\ {\text{e}}^{{\eta }_{\mathrm{3}}x}\\ {\text{e}}^{{\eta }_{\mathrm{4}}x}\end{array}\right\}.$

(31)

The relations between coefficients A _i , B _i , C _i , and D _i can be deduced. Then, we can obtain the eigenvalue μ from the following algebraic equation:

$\left({\eta }_{\mathrm{1}}-{\eta }_{\mathrm{2}}\right)\left({\eta }_{\mathrm{3}}-{\eta }_{\mathrm{4}}\right)+\left({\eta }_{\mathrm{1}}-{\eta }_{\mathrm{4}}\right)\left({\eta }_{\mathrm{2}}-{\eta }_{\mathrm{3}}\right)\cosh \left[\left({\eta }_{\mathrm{1}}-{\eta }_{\mathrm{2}}-{\eta }_{\mathrm{3}}+{\eta }_{\mathrm{4}}\right)h\right]+({\eta }_{\mathrm{1}}-{\eta }_{\mathrm{3}})({\eta }_{\mathrm{4}}-{\eta }_{\mathrm{2}})\cosh \left[({\eta }_{\mathrm{1}}-{\eta }_{\mathrm{2}}+{\eta }_{\mathrm{3}}-{\eta }_{\mathrm{4}})h\right]=\mathrm{0}.$

(32)

The transcendent equation is solved by the Newton-Rphson method. The origin solution corresponding to each eigenvalue μ _i is

${\mathbf{v}}_i^{\prime }={\mathbf{M}}_t^{\prime }{\text{e}}^{{\mu }_{i}z}{\mathbf{\Phi }}_i.$

(33)

We seek the particular solutions for the FGM beams subjected to arbitrary tractions on the lateral surfaces. The conditions at the upper and lower surfaces are

$\begin{array}{c}x=-h:({\tilde{\sigma }}_x=)E_{\mathrm{0}}\frac{\partial \tilde{u}}{\partial x}+\nu \tilde{\sigma }=q_{\mathrm{1}}\left(z\right){\text{e}}^{\beta h},\hspace{1em}\hspace{1em}\tilde{\tau }=\mathrm{0},\\ x=h:({\tilde{\sigma }}_x=)E_{\mathrm{0}}\frac{\partial \tilde{u}}{\partial x}+\nu \tilde{\sigma }=q_{\mathrm{2}}\left(z\right){\text{e}}^{-\beta h},\hspace{1em}\hspace{1em}\tilde{\tau }=\mathrm{0}.\end{array}$

(34)

The solution of $\mathbf{H}\stackrel{\sim }{\mathbf{\Phi }}=k_{\mathrm{3}}{\mathbf{\Phi }}_{\mathrm{0}}^{(\mathrm{3})}$ can be solved as follows:

$\left(\tilde{\mathbf{\Phi }}=\right)\left\{\begin{array}{c}\tilde{w}=\mathrm{0}\\ \tilde{u}=k_{\mathrm{3}}\left\{\frac{\mathrm{1}}{\mathrm{24}}{\nu }^{\mathrm{2}}{x}^{\mathrm{4}}-\frac{\mathrm{1}}{\mathrm{6}}{\nu }^{\mathrm{2}}{A}_{\mathrm{0}}{x}^{\mathrm{3}}+\left(\mathrm{1}-{\nu }^{\mathrm{2}}\right)\left[\left(c_{\mathrm{1}}x+\frac{c_{\mathrm{1}}}{\beta }-c_{\mathrm{3}}\right)\frac{{\text{e}}^{-\beta x}}{E_{\mathrm{0}}\beta }-\frac{\mathrm{1}}{\mathrm{2}{\beta }^{\mathrm{2}}}{x}^{\mathrm{2}}+\left(\frac{\mathrm{2}}{{\beta }^{\mathrm{3}}}+\frac{A_{\mathrm{0}}}{{\beta }^{\mathrm{2}}}\right)x\right]\right)\\ \left(\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{\mathrm{2}\nu \left(\mathrm{1}+\nu \right)}{\beta }\left[\frac{\mathrm{1}}{\mathrm{6}}{x}^{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{2}\beta }{x}^{\mathrm{2}}-\frac{A_{\mathrm{0}}}{\mathrm{2}}{x}^{\mathrm{2}}+\frac{c_{\mathrm{1}}}{E_{\mathrm{0}}\beta }{\text{e}}^{-\beta x}\right]-c_{\mathrm{2}}\nu \right\}\\ \tilde{\sigma }=k_{\mathrm{3}}{E}_{\mathrm{0}}\left\{\left(c_{\mathrm{3}}-c_{\mathrm{1}}x\right)\frac{\nu }{E_{\mathrm{0}}}{\text{e}}^{-\beta x}+c_{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\nu A_{\mathrm{0}}{x}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{6}}\nu x^{\mathrm{3}}+\frac{\mathrm{2}\nu }{{\beta }^{\mathrm{3}}}+\frac{A_{\mathrm{0}}\nu }{{\beta }^{\mathrm{2}}}\right)-\frac{\nu }{{\beta }^{\mathrm{2}}}x\\ \left(\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{\mathrm{2}\left(\mathrm{1}+\nu \right)}{\beta }\left[\frac{\mathrm{1}}{\mathrm{2}}{x}^{\mathrm{2}}-\frac{\mathrm{1}}{\beta }x-A_{\mathrm{0}}x-\frac{c_{\mathrm{1}}}{E_{\mathrm{0}}}{\text{e}}^{-\beta x}\right]\right\}\\ \tilde{\tau }=\mathrm{0}\end{array}\right\},$

(35)

in which k₃ = (q₁(z) − q₂(z))/2[c₁h + E₀hcosh ⁡ (βh)/β² − (2E₀ + A₀E₀β)sinh ⁡ (βh)/β³] and c₃ = [(q₁(z)e^βh − q₂(z)e^−βh)/2 − k₃h(E₀/β² + c₁cosh ⁡ (βh))]/[k₃sinh ⁡ (βh)].

The particular solution of the plane problem can be constructed as

$\tilde{\mathbf{v}}={\mathbf{M}}_t^{\prime }\left\{\tilde{\mathbf{\Phi }}+k_{\mathrm{3}}\left[z{\mathbf{\Phi }}_{\mathrm{0}}^{\left(\mathrm{3}\right)}+\left(\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{2}\right)}+\frac{\mathrm{1}}{\mathrm{6}}{z}^{\mathrm{3}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{24}}{z}^{\mathrm{4}}{\mathbf{\Phi }}_{\mathrm{0,2}}^{\left(\mathrm{0}\right)}\right)\right)\right)\left(\left(+A_{\mathrm{0}}\left(\frac{\mathrm{1}}{\mathrm{2}}{z}^{\mathrm{2}}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{1}\right)}+\frac{\mathrm{1}}{\mathrm{6}}{z}^{\mathrm{3}}{\mathbf{\Phi }}_{\mathrm{0,1}}^{\left(\mathrm{0}\right)}\right)\right]\right\}.$

(36)

The complete analytical solution is also a combination of the particular solution and the solutions corresponding to zero and general eigenvalues as follows:

${\mathbf{v}}^{\prime }=\tilde{\mathbf{v}}+m_{\mathrm{1}}{\mathbf{v}}_{\mathrm{0,1}}^{\prime \left(\mathrm{0}\right)}+m_{\mathrm{2}}{\mathbf{v}}_{\mathrm{0,1}}^{\prime \left(\mathrm{1}\right)}+m_{\mathrm{3}}{\mathbf{v}}_{\mathrm{0,2}}^{\prime \left(\mathrm{0}\right)}+m_{\mathrm{4}}{\mathbf{v}}_{\mathrm{0,2}}^{\prime \left(\mathrm{1}\right)}+m_{\mathrm{5}}{\mathbf{v}}_{\mathrm{0,2}}^{\prime \left(\mathrm{2}\right)}+m_{\mathrm{6}}{\mathbf{v}}_{\mathrm{0}}^{\prime \left(\mathrm{3}\right)}+{\displaystyle \sum }_{i=\mathrm{1}}^{\infty }\left({\zeta }_i{\mathbf{v}}_i^{\prime }+{\zeta }_{-i}{\mathbf{v}}_{-i}^{\prime }\right).$

(37)

When a large number of eigensolutions corresponding to nonzero eigenvalues are involved to obtain accurate results, we reformulate the form of β-class eigensolutions to achieve the stability of numerical calculation.

4. Numerical Examples

4.1. Example 1

The simple supported isotropic beam of thickness 2h = 1 m and length-to-thickness ratio l/(2h) = 5 is subjected to prescribed normal tractions q₁ = − 1 N/m² at its upper surface. Three cases of plane beams are discussed in the analysis: (1) elasticity materials: Young's modulus is E₀ = 2.0 × 10¹¹ N/m², (2) axial FGMs: Young's modulus varies exponentially along z with its value at z = 0 being E₀ = 2.0 × 10¹¹ N/m². The material gradient index takes two values as αh = 0.01 and αh = 0.1, and (3) transverse FGMs: Young's modulus varies exponentially along x with its value at x = 0 being E₀ = 2.0 × 10¹¹ N/m². The material gradient index takes two values as βh = 0.01. and βh = 0.1. For all the three materials models, the Poisson's ratio takes the value ν = 0.29.

The present results for homogenous materials (seen in Figure 3) agree well with the solutions in Timoshenko and Goodier [31]. Meanwhile, the symplectic method improves the efficiency and accuracy of plane problem with arbitrary lateral loads. Figures 4, 5, 6, and 7 depict the stresses and displacement contours of plane beams which is functionally-graded in the length or thickness direction. It can be seen from Figures 4 and 6 that the well matchs are obtained for the smaller inhomogenous parameters αh = 0.01 and βh = 0.01, respectively. We clearly achieved the highly accurate local stress and displacement distributions for both FGMs models. It should also be noted from Figures 4–7 that, the materials inhomogeneity do have effects on the distributions of stress and displacement. In contrast, the inhomogeneity in the length direction plays more evident effects on both stresses and displacements. Take the normal stress for example, the normal stress at the vicinity of the lateral surfaces decrease gradually with the increasing of gradient index αh (seen in Figures 4 and 6).

Figure 3:

Contours of stress and displacement for homogenous beam. (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

Figure 4:

Contours of stress and displacement for axial FGM beam (αh = 0.01). (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

Figure 5:

Contours of stress and displacement for axial FGM beam (αh = 0.1). (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

Figure 6:

Contours of stress and displacement for transverse FGM beam (βh = 0.01). (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

Figure 7:

Contours of stress and displacement for transverse FGM beam (βh = 0.1). (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

4.2. Example 2

The simple supported isotropic beam of thickness 2h = 1m and length-to-thickness ratio l/(2h) = 5 is subjected to prescribed normal tractions q₁ = − q₀(z/l) (q₀ = 0.1 MPa) at its upper surface. We also consider three materials models as that in example 1.

For comparison, Figure 8 first shows the stress and displacement distributions of the homogenous beam. We take eighty terms of eigensolutions for nonzero eigenvalue in the symplectic expansion. Figure 9 depicts the contours of the analytical stress and displacement for axial functionally graded beam, while Figure 10 illustrates the complete stress and displacement distributions of transverse functionally graded beam. It can be seen that, when αh = 0.01 or βh = 0.01, the distributions of physical quantities all tend those for the homogeneous materials. It is verified that the whole stress and displacement fields of the functionally graded beam can be obtained accurately in the symplectic framework.

Figure 8:

Contours of stress and displacement for homogenous beam. (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

Figure 9:

Contours of stress and displacement for axial FGM beam (αh = 0.01). (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

Figure 10:

Contours of stress and displacement for transverse FGM beam (βh = 0.01). (a) σ _z ; (b) τ _xz ; (c) u _z ; (d) u _x .

5. Conclusion

On the basic of symplectic analysis, the exact solutions for stress and displacement fields in homogeneous and functionally graded beams are obtained in a rational, well-structured, and efficient approach. Particular solutions for two exponential FGM models subjected to arbitrary lateral loads are determined in a systematic way. Combined with the homogeneous solutions that correspond to the eigensolutions of particular eigenvalues (0 or α) and general eigenvalues, the complete analytical solutions are deduced. Numerical results show that the symplectic approach is effective in predicting accurate local stress and displacement distributions. The approach is suited for plane beams with any combinations of simply supported, clamed, free, and sliding-contact ends. The similar analytical procedure is also valid for bidirectional functionally graded beams and functionally graded curved beams in which the material properties vary according to a power law.