Extension of Boley’s method to functionally graded beams

The objective of this contribution is the computation of the Airy stress function for functionally graded beam-type structures subjected to transverse and shear loads. For simplification, the material parameters are kept constant in the axial direction and vary only in the thickness direction. The proposed method can be easily extended to material varying in the axial and thickness direction. In the first part an iterative procedure is applied for the determination of the stress function by means of Boley’s method. This method was successfully applied by Boley for two-dimensional (2D) isotropic plates under plane stress conditions in order to compute the stress distribution and the displacement field. In the second part, a shear loaded cantilever made of isotropic, functionally graded material is studied in order to verify our theory with finite element results. It is assumed that the Young’s modulus varies exponentially in the transverse direction and the Poisson ratio is constant. Stresses and displacements are analytically determined by applying our derived theory. Results are compared to a 2D finite element analysis performed with the commercial software ABAQUS. It is found that the analytical and numerical results are in perfect agreement.


Introduction
Functionally graded materials (FGMs) are composite materials where the material properties vary throughout the body. The simplest way of designing a graded material is by taking isotropic homogenous layers with different material parameters, which are perfectly bonded together. This special case is a laminate, where the material properties vary discontinuously. The first FGMs were compositions of metallic and ceramic materials where the thermal expansion coefficients varied continuously. A prominent example of a previous FGM application was the space shuttle where a material made of ceramic and metal provided thermal protection and increased the load-carrying capability during re-entering the earth's atmosphere. One of the most significant aspects of designing FGMs is to close the gap in the properties of different materials continuously. For an introduction the reader is referred to Mahamood [1] and Udupa [2] who give a compact overview of the topic. A mechanical advantage of FGMs is the stress continuity (where the stress distribution is nonlinear in general) and the lower stress gradient inside the material which reduce peak stresses, compared to rigidly bonded plies. The literature on FGMs, where the focus is laid on the derivation of a mechanical model which is valid for arbitrary material property variations, is limited. Here, beam types structures are investigated and a very general modeling strategy is presented to find accurate beam models for a large class of FGMs. Mian [3] focuses on a three-dimensional solution with transversal-dependent elasticity parameters by solving two-dimensional equations for an equivalent homogeneous plate with averaged elastic parameters. The work of Sankar [4] contains the study of a sinusoidally-loaded, simply supported beam with orthotropic material parameters that vary exponentially. The Bernoulli-Euler theory solution is obtained by starting with the assumed axial and transversal displacements. Strains and stresses are computed afterwards. The Airy stress function for an anisotropic FGM beam subjected to normal and axial loads is computed by Ding [5] who assumes the stress function as a sum of products of the independent variables. The suggested method of [5] is applied similarly by Zong [6], Hu [7] and Yang [8]. A cantilever subjected to a linearly distributed normal force is analyzed in [7], whereas in [8] the Airy stress function for a two-layer FGM cantilever with concentrated end loads is computed. A number of various different loading and boundary value problems is solved in [5], and a finite element analysis is performed to compare analytical and numerical results in [5] and [7]. Sakurai [9] developed the Airy stress function for exponentially graded orthotropic beams by an infinite series of products of independent variables. The stress distributions are presented for a simply supported beam subjected to a sinusoidally distributed normal force and a cantilever under uniform pressure. A two-dimensional solution for pure bending and tension is derived in an analogous way by [5][6][7][8] and by Chu [10] for transversal and axial isotropic exponentially graded material. The elastic buckling and bending of a functionally graded, Timoshenko theory-based porous beam, with transversal cosine-dependent material parameters is studied by Chen [11]. Hadji [12] analyzes the static bending and free vibrations of a shear-deformable FGM beam where the strains are assumed to be exponentially distributed in the transverse direction and the material parameters are powerlaw distributed. A similar analysis of a shear-deformable FGM beam is found in Guenfod [13]. In the work of Xia [14] a comparison of an FGM Reddy-Bickford beam versus a homogenous Bernoulli-Euler beam is presented: The investigation starts with strain assumptions where the transversal strain is disregarded, and as a result, the solutions of determinate and indeterminate beams are shown. A higher-order shear-deformable finite beam element for an FGM sandwich structure is derived by Li [15] where higher-order strain assumptions are taken into account, but the transversal strain is neglected as in [14]. Our motivation in this contribution is to provide an analytical method (under classical plane-stress assumption) to determine a beam model without any a priori assumptions on the displacement or the stress field. This strategy dates back to works of von Karman [16] and Sewald [17] who showed that the curvature of an isotropic beam does depend not only on the bending moment, but also on its even-numbered spatial derivatives. In our contribution, we adopt an iterative solution method for the bi-potential equation, which arises by evaluating the compatibility equations when an isotropic, homogeneous material is considered, see Boley and Tolins [18]. The method of Boley and Tolins [18] consists of a step-by-step procedure to compute the Airy stress function as the solution of the bi-potential equation. The result of the first iteration step yields the classical Bernoulli-Euler (BE) solution, whereas the higher-order iterations are correction terms which take into account the thickness-to-length ratio and become more relevant for thicker beams. Furthermore, one can derive more accurate distributions for the stress components, and cross-sectional warping effects are automatically accounted for. Boley's method is extended up to the fourth iteration by Irschik [19], who calculates the bending moment for different thicknessto-length ratios to demonstrate the influence of different boundary conditions for a statically indeterminate beam. In a further contribution, the thermoelastic effect is considered in Boley [20], see also the book on thermal stresses by Boley and Weiner [21]. Krommer and Irschik [22] use this method to approximately solve the charge equation of electrostatics for the electric potential if a certain displacement field is assumed. For piezoelectric materials, Schoeftner and Benjeddou [23] give an iterative solution for the electric potential and the Airy stress function in case of piezoelectricity. It is found that the compatibility and the charge equation, which both depend on the stress and the electric potential, are coupled. For a simply supported piezoelectric bimorph with sinusoidal mechanical and electrical actuation, the result is compared to (2D) analytical results. Motivated by the latter contributions, the iterative method from Boley is extended to FGMs. Instead of the bi-potential equation, we derive a fourth-order partial differential equation for the Airy stress function from the compatibility equations. Hereby, no a priori assumptions are made, neither on the displacement field nor on the stress field. The Airy stress function is computed iteratively, and the external loads are replaced by the stress resultants such as normal force, shear force and bending moment. As previously mentioned the first iteration solution is denoted as an elementary solution, and the higher-order solutions can be interpreted as an extension of the Bernoulli-Euler result for a graded material. The results of the further iterations contain the thickness-to-length ratio and further components of the compliance matrix as well. Numerical results are given for an isotropic FGM cantilever with exponentially varying Young's modulus subjected to a shear load. The results are presented to demonstrate the influence of the graded material on stresses and deformations. The beam model is verified by a two-dimensional (2D) finite element calculation performed with ABAQUS where analytical and numerical results are in excellent agreement.

Problem statement
The following presented solution, based on Boley's method [18], allows the determination of the Airy stress function of a functionally graded beam. The tractions acting on the upper and on the lower surface are polynomials in x. The cross section of the beam is rectangular and constant along the beam length L, the thickness is 2c, and the width is b is set to unit one. It is assumed that the thickness of the beam is sufficiently larger than b, i.e., that the rectangular beam is narrow, such that the state of plane stress can be taken into account. A global Cartesian (x,y)-coordinate system is attached to the mid-span axis, where the axial coordinate is denoted by x, and the transverse coordinate by y (Fig. 1).
The traction vectors at the upper surface at y = c and at the lower one at y = −c read In Eq. (1), the subscripts u and l denote the upper and lower surface. The material characteristics of the beam are described by the compliance matrix S, where we restrict ourselves to transversal dependency only in this contribution: with Poisson's ratio ν and the y-dependent Young's modulus E(y).

Method of solution
Under the prescribed constitutions, we have to find solutions of the plane linear theory of elasticity, which satisfies exactly the equilibrium equations, the strain-displacement relations, the stress boundary conditions at the upper and lower surface of the beam at y = ±c and the stress-strain relations, see Eq. (7). The normal stresses are σ x and σ y , τ xy is the shear stress, ε x and ε y are normal strains, γ xy is the shear angle, u is the axial, and v is the transversal displacement, respectively. Note that we take advantage of the following definitions for the partial derivatives with respect to x and y: If the stresses in Eq. (3) are expressed as derivatives of the introduced Airy stress function ψ(x, y), see Timoshenko and Goodier [24], the equilibrium equations in Eq. (3) are automatically satisfied. The stress-strain relation expressed in Voigt's notation reads For all further iterations, (i.e., ψ i>1 ) the stress boundary conditions are trivial: for i ≥ 2.

First iteration
For the computation of the first potential ψ 1 we consider Eqs. (12) and (13), and we find for the first iteration Equation (16) is a fourth-order ordinary differential equation (ODE) with variable coefficients. A solution of Eq. (16) is see Lang and Pucker [28], with four remaining constants k j and functions g j (y) which are results when solving Eq. (16). To ensure the x, y dependency of ψ 1 the constants k j are replaced by functions f j (x) which leads to Integration of the first row of Eq. (14) once, and the second of Eq. (14) twice with respect to x leads to the following vector-matrix notation: The vector q contains the integrated surface loads, A denotes a 4×4 coefficient matrix, f is the vector with the unknowns f j (x), and g denotes the vector of fundamental solutions g j (y). Taking into account Eqs.

Stress function in terms of stress resultants
The remaining constants C 1 , C 2 , C 3 and K 1 , K 2 , K 3 in Eq. (23) are computed by introducing the stress resultants, such as bending moment M(x), shear force Q(x) and the normal force N (x), Substituting Eqs. (6) and (24) into Eq. (25), it follows that where The results in Eq. (29) correspond to the sign convention for a right-handed coordinate system x yz as that shown in Fig. 2 for a differential beam element (i.e., at x + dx the positive normal force N (x) and the positive shear force Q(x) are parallel to the x-and y-directions, respectively, but the positive bending moment M(x) is in opposite direction to the positive z-direction). It is noted that if the positive directions of the coordinates do not match with those in Fig. 2, one observes that the resulting equilibrium equations on beam level and those in Eq. (29) will differ in signs.

Second iteration
To compute the second stress function ψ 2 , we set i = 2 in Eq. (12), and by considering Eq. (13) we find with the trivial boundary conditions from Eq. (15). The solution of ψ 2 is a long expression in general, but can be determined by symbolic computation software.

Further iterations
We repeat the prescribed procedure and obtain for the third iteration step for i = 3 in Eq. (12) The solution of ψ 3 depends on the previous solutions ψ 1 , ψ 2 . As for the second iteration, trivial boundary conditions are assumed, see Eq. (15).

Shear-loaded beam
The developed solution procedure is verified by the following example. An FGM beam is investigated which is subjected to a shear force at the upper surface at y = c, see Fig. 3. The boundary tractions are as follows: We assume an isotropic graded material with an exponential variation of the Young's modulus in transverse direction where the compliance matrix of Eq. (2) becomes with an introduced nondimensional gradient α and the Young's modulus E 0 at the midspan axis at y = 0. The stress-strain relations become Reinserting the relations of Eq. (34) into Eq. (8) leads to the iteration rule for computation of the stress function ψ with the conventions of Eq. (13) to be considered. If we set α = 0, which corresponds to isotropic material, we obtain the iteration rule as stated in the paper of Boley [18]: The solution of Eq. (37) follows as The unknown functions f j (x) can be determined from the tractions at y = ±c, see Eq. (32); with the coefficient matrix A and the load vector q we can write: ⎡ The result of ψ 1 is lengthy, and to simplify the expression, we introduce entities corresponding to the elementary theory: the so-called effective axial rigidity [E A] eff , the flexural rigidity [E I ] eff and the location of the neutral axis c 0 . For computing these elementary entities see e.g. Ziegler [26]: The normal stress σ 1x = ψ 1,yy as a result of the first iteration in terms of the introduced entities of Eqs. (40) and (41) reads and the shear stress τ 1xy = −ψ 1,xy For the normal stress σ 1x there remains a linear distribution due to the bending moment M(x) and a constant part due to the normal force N (x):

Second iteration
Inserting ψ 1 into Eq. (35) leads to the computation rule for ψ 2 : following the boundary conditions (15). If a linear distribution of the shear force is assumed, the iteration procedure ends at i = 2, which corresponds to the second-order polynomial of the load vector q, see Eq. (39).

Displacements
The displacements u(x, y) and v(x, y) are computed by integration of the normal strains ε x and ε y and then adjusting the arbitrary functions F(x) and G(y) so that u ,y + v ,x − γ xy = 0 is satisfied. During this procedure, three constants preventing rigid body motion, u 0 , v 0 and r 0 , remain, which have to be determined by setting up a boundary value problem, see Sect. 5. For more details of this standardized procedure see, e.g., Timoshenko and Goodier [24]. The expressions of the displacements are very lengthy, so we show the displacements of the first iteration u 1 and v 1 only: higher-order terms  Table 2.

Numerical example
The previously determined results are applied to a cantilever beam with length L, thickness 2c and a thicknessto-length ratio λ = 2c/L = 1/10, see the scaled sketch in Fig. 4. The beam is loaded by a linearly distributed shear force at the upper surface, and the stress distributions and the displacements are determined due to variations of the Young's modulus in the y-direction. The graded material is assumed for a constant Poisson's ratio ν and an exponentially distributed Young's modulus E(y) = E 0 e αy/c , see Eq. (33). The derived results are compared to a two-dimensional (2D) finite element analysis (FEA) performed with the commercial software ABAQUS. The example beam (length L = 50 mm and thickness 2c = 5 mm) is meshed with CPS8, an 8-node biquadratic plane stress quadrilateral element. To model the FGM characteristics of the beam, the cross section is separated into 50 isotropic linear elastic layers with constant Poisson's ratio ν and constant Young's modulus E n . The layer-wise Young's modulus E n is computed from  with β n = (2n − 51)/50, 1 ≤ n ≤ 50. For the FEA model, we only investigate the case α = 1/2 and one finds the parameters for the FEA collected in Table 1.
The three kinematical constants u 0 , v 0 and r 0 (see Eqs. (47) and (48)) are computed by three kinematical boundary conditions, in order to prevent two rigid body translations and one rigid body rotation. It holds that and The latter constraint for the cross-sectional rotation was suggested by Szabo [27] and applied in the work of Gahleitner and Schoeftner [25] in order to get similar results as for a rigid clamped end with u(0, y) = v(0, y) = 0. In our simulation the FE model is subjected to the same tractions as the analytical model (see Fig. 3). The upper surface traction t u at y = c contains the imposed shear load n u (x), while the lower layer traction t l at y = −c vanishes. It is important to note that at the clamped end at x = 0 the traction components of t 0 read σ x (0) = ψ 1,yy (0, y) + ψ 2,yy (0, y) and τ xy (0) = −ψ 1,xy (0, y), see Eqs. (42) and (43) for the first iterations and the solution of (46) for the second iteration (lengthy expression). For the free end at x = L the x-and the y-components of the traction t L read σ x (L) = ψ 2,yy (L , y) (see Eq. (46), and the bending moment vanishes, consequently ψ 1,yy (L , y) = 0 and −ψ 1,xy (L , y) = 0) (Fig. 5).

Stress field
The stress resultants can be determined by integrating the beam equilibrium relations, see Eq. (29): Considering 3 boundary conditions at x = L, one finds for M, Q and N : where the nondimensional axial coordinate ξ = x/L is introduced. In Figs. 7a, b, normal stresses are plotted as functions of Young's modulus E(η) = E 0 e αη with α = {0, 1/10, 1/5, 1/2}. Figure 7a shows the stress as a function of the normal force N (x), and Fig. 7b shows the stress as a function of the bending moment M(x). It is noted that the normalized axial stress is plotted: the reference stress σ * = |M|c/I = 75 N/mm 2 denotes the absolute axial stress at the upper fiber at η = 1 and ξ = 1/2 due to the bending moment M(1/2) of a homogeneous configuration with α = 0 and the moment of inertia I = 2bc 3 /3, see the linear portion of Eq. (45). It can be seen that the nonlinearity due to the normal force in Fig. 7a is negligible also for highly graded materials (i.e., even for α ≥ 1/2 because the nonlinear part in Fig. 7a is relatively small compared to the mean value 0.33). On the other hand the nonlinearity due to the bending moment can be observed for α > 1/5. In Fig. 7c the total normalized axial stress σ x /σ * = (σ x N + σ x M )/σ * is shown and compared to the FEA (black points for α = 1/2) showing an excellent agreement of analytical and numerical results. It has to be noted that for α = 1/2, the maximum normal stress σ x max is 27% larger than σ * . The axial stress of the higher iterations σ 0 x (derived from ψ i with i > 1) is plotted in Fig. 7d. It is noted that this is a self-equilibrated stress, which means that stress resultants vanish by integration over the cross-sectional area. One observes that the quantity of σ 0 x is very low (i.e., below 0.75% of σ * ) and the influence of the residual stress σ 0 x is practically negligible.
In Fig. 8a the normal stress σ y is plotted. This stress is often denoted as compressive stress in the literature, and it vanishes at both surfaces, and also its derivative vanishes at η = −1: This follows directly from the boundary conditions because σ y ∝ ψ = 0 at η = ±1 and σ y,y ∝ ψ ,y ∝ τ xy = 0 at η = −1 hold. In relation to the stress σ 0 x presented in Fig. 7d, the quantity of σ y is of very low influence: in general it holds that O(σ y ) = O(σ x )λ 2 , i.e., the order of magnitude of the compressive stress is lower by a factor λ 2 than the axial stress. The shear stress distribution is shown in Fig. 8b: It vanishes at the lower surface at η = −1 and fulfils the normalized traction load n 0 /σ * x at the upper surface at η = 1. The FEA results for σ y and τ xy are given for α = 1/2 and those match perfectly with our analytical solution.
Some final remarks concerning the stresses are made: The nondimensional residual stress σ 0 x /σ * , in Fig. 7d, and the compressive stress σ y /σ * in Fig. 8a are lower-order terms (quadratic order with respect to thicknessto-length ratio λ) in comparison with the nondimensional axial stress σ x /σ * (e.g., σ ymax /σ * = 0.0033). Independent of the thickness-to-length ratio self-equilibrated stresses are of low impact and the compressive stress is low. In addition to the beam considered in this study (λ = 1/10) this influence remains negligible even for a thicker beam (λ = 1/5: the maximum value is σ y /σ * ≈ 0.014).

Displacements
The displacements are presented for different Young's moduli in Fig. 9.  1 and v 1 (red curve for α = 1/2). It has to be noted that for our example beam, where a thickness-to-length ratio of λ = 1/10 is considered which refers to a moderate thick beam, the higher-order solution is almost equal to the first-order solution. The obtained analytical results are in a very good accordance with the FEA solution. One observes that the higher the material is graded, the lower is the horizontal deflection at η = 0 (Fig. 9a). This also holds for the vertical deflection (Fig. 9b). Furthermore, considering only the first iterations (red curves for α = 1/2), one finds a very good agreement with FEA results (black dots): The difference is less than < 1%. Taking into account also the higher-order terms (black), the analytical solutions (black) practically coincide with FE results. In Figs. 10a and b the normalized horizontal and vertical deflections are shown over the cross section at ξ = 1/2. Again, the vertical tip deflection of an isotropic homogeneous beam v * = 1.07 mm is used for normalization. The FEA results are performed for the material gradient α = 1/2.
The influence of the higher-order terms, besides the fundamental solution from Eqs. (47) and (48), is marginal. At the lower fiber at η = −1, the normalized axial deflection is 0.038 for all material configurations, see Fig. 10a. The distribution remains almost linearly, but it decreases for higher graded materials. Again, the FEA results agree with the analytical results, regardless of whether the higher-order terms are considered or not (c.f. red, black and pointed curves). The vertical cross-sectional deflection is almost constant with respect to the thickness (Fig. 10b). The first-order solution only differs slightly from the higher order and the FEA result (error < 0.1%).
As a final statement the authors suggest that for thin beams (λ < 1/20) the first-order results yield accurate solutions for industrial use, see Eqs. (44) and (45) for the normal stress distributions and column 1 of Table 2 in the Appendix for the displacements.

Comparison to literature
An additional verification of our theory is performed in this section, where we compare our results to an example given in Chu [10]. In order to obtain a stress function in the case of vanishing tractions with t u = t l = {0, 0} T , (1), the computation procedure for the stress function is similar to Sect. 4, but the load vector is adjusted A cantilever beam is investigated which is loaded by the end moment B, the shear force V and the normal force H at x = L (see Fig. 11).
The constants C 1 , C 2 , C 3 and K 1 , K 2 , K 3 are related to three constants

Conclusion
In the present contribution, an iterative solution procedure has been derived for the computation of the Airy stress function of a functionally graded beam with rectangular cross section and transversely varying material parameters. For this purpose, Boley's method, which was initially developed for deriving analytical beam solutions for isotropic thick beams, is here extended. We assume a functionally graded material (FGM) and do not take advantage of a priori axiomatic assumptions, neither for the stress field nor for the displacement field. Explicit results are given for a cantilever beam subjected to a linearly distributed shear load. For other kinematic (no matter if statically determinate or indeterminate beams) and load boundary conditions the proposed method can be easily adjusted. The analytical solution is obtained by solving algebraic equations for the first iteration and simple ordinary differential equations for the higher iteration. Finally, a numerical verification of our functionally graded beam model is presented. The solution is compared to two-dimensional (2D) finite element results performed in ABAQUS, showing that the analytical results are in excellent agreement with the ABAQUS outcome. The obtained solution is an extension of the Bernoulli-Euler (BE) theory for FGM which takes into account the thickness-to-length ratio of the structure; hence, it is valid for thicker structures. Table 2 Midspan axis displacements for an isotropic FGM beam with Young's modulus E(y) = E 0 e y/2c , constant Poisson's ratio ν, length L, thickness 2c subjected to a linear distributed shear force with n u (x) = n 0 (x/L − 1) Six constants, u 0 , v 0 , r 0 , M 0 , Q 0 and N 0 , are to be determined by setting up three appropriate boundary conditions at each lateral beams ends