An Efficient Beam Element Based on Quasi-3D Theory for Static Bending Analysis of Functionally Graded Beams

In this paper, a 2-node beam element is developed based on Quasi-3D beam theory and mixed formulation for static bending of functionally graded (FG) beams. The transverse shear strains and stresses of the proposed beam element are parabolic distributions through the thickness of the beam and the transverse shear stresses on the top and bottom surfaces of the beam vanish. The proposed beam element is free of shear-looking without selective or reduced integration. The material properties of the functionally graded beam are assumed to vary according to the power-law index of the volume fraction of the constituents through the thickness of the beam. The numerical results of this study are compared with published results to illustrate the accuracy and convenience rate of the new beam element. The influence of some parametrics on the bending behavior of FGM beams is investigated.


Introduction
Functionally graded (FG) materials (FGM) are one of the advanced composites. In which the material properties of FGM vary continuously through one or more directions. A typical, the material properties of an FGM beam, plate and shell varies continuously through the thickness direction. Due to their advantages, the FGMs have used widely in many fields such as civil, aerospace, automobile, engineering, nuclear power plants and so on [1,2]. Since then, many scientists have been focused on the mechanical analysis of FG beams, plates and shells. In which, they used several theories and methods, for instance, analytical and numerical methods based on Euler-Bernoulli theory, Timoshenko theory or first-order shear deformation theory (FSDT), higher-order shear deformation theory (HSDT), Quasi-3D theory and Carrera Unified Formulation (CUF).
Sankar [3] developed an elasticity solution to analyze a simply supported FG beam subjected to transverse distribution loading. In his work, Sankar developed new beam theory which was similar to the Euler-Bernoulli beam theory. Zenkour [4] analyzed an exponentially graded thick rectangular plate using both 2-D and 3-D elasticity solutions. Zhong et al. [5] analyzed a cantilever FG beam using a general two-dimensional solution. The free vibration and buckling analysis of FG beams under mechanical and thermal loads were investigated by Trinh et al. [6] using the analytical method based on the state space approach and higher-order beam theory.
However, in HSDT and Quasi-3D theory, the displacement field is considered by the existence of the higher order derivative of the deflection of transverse. So that it involves the development of a C 1 continuity element, which can cause difficulty to originate the second derivative of deformation in FEM. To overcome these continuity issues, Hermite interpolation functions with C 1 elements and some C 0 approximations have been adopted. Chakraborty et al. [37] developed a new beam element based on the FSDT and an exact solution of the static part of the governing differential equations for analysis of FGM structures. Nguyen et al. [38] applied the Timoshenko beam model and FEM for dynamic response of bi-directional FG beams subjected to moving load. Khan et al. [39] investigated the static bending and free vibration of FG beams using FEM. Heyliger [40] developed a higher order beam finite element for bending and vibration of beams. Kapuria et al. [41] studied bending behavior and free vibration response of layered FG beams using a third order zigzag theory and FEM. Based on refined shear deformation theory, Vo et al. [42] developed a finite element model to analyze free vibration and buckling of FG sandwich beams. Moallemi-Oreh et al. [43] used FEM for stability and free vibration analysis of the Timoshenko beam. Pascon [44] analyzed FG beams with variable Poisson's ratio using FEM. Yarasca et al. [45] studied FG sandwich beams using Hermite-Lagrangian finite element formulation. The use of higher-order shape function will cost much computation effort in comparison with linear shape function. Furthermore, the linear shape function is simpler in formulation and transformation than higher-order shape function. However, plate and beam element using linear shape function are mainly developed based on FSDT and HSDT. To author's knowledge, there is currently no beam element using linear shape function based on a Quasi-3D theory. Therefore, the development of a beam element using linear shape function based on a Quasi-3D theory is necessary.
This paper presents a new beam element based on Quasi-3D theory, which only requires C 0 shape functions. The organization of this study is as follows. Firstly, Section 2 defines the model and material properties of FG beams. In Section 3, the governing equations of FG beams based on Quasi-3D theory are given. The finite element formulations of the proposed beam element are presented in Section 4. In Section 5, some example problems are carried out to show the convergence and accurate rate of new beam element in comparison with published data. Then, the static bending behaviors of FG beams using the proposed beam element are studied. The influences of the distribution of materials properties, length-to-thickness ratio, boundary conditions and effect of normal strain are investigated. Finally, in the conclusion section, some remarks on the proposed beam element are given.

Functionally Graded Material
Consider an FG beam as shown in Figure 1, the length of the beam is L, the width of the beam is b, and the thickness of the beam is h. The Young's modulus varies continuously through the thickness of the beams with a power law distribution [19,20]: in which subscript m denotes the metallic component and c denotes the ceramic component, E m , E c are respectively Young's modulus of the metal and ceramic, p is the power-law index. In this study, the Poisson's ratios ν of both components are assumed to be constant and equal.

Governing Equations
The displacements of a point of the beam are expressed by

Governing Equations
The displacements of a point of the beam are expressed by The functions f 1 (z), f 2 (z) are given by Shi [13] The strain field is obtained as follows in which the symbol (,) means the derivatives with respect to the quantity following it and the symbol ( ) means the derivatives with respect to z direction. Rewrite the strain components in the short form as follows where Rewrite the transverse shear strain γ xz as follows The transverse shear strain is assumed to have a quadratic distribution across the thickness of the beam. In addition, the transverse shear strain equals to zero at the top and bottom surfaces of the beam. These conditions lead to The constitutive relations between the stress field and the strain field are expressed as follows In this study, Young's modulus E of FGM is a function of the coordinate, whereas, the Poisson's ratio is assumed to be constant and equal, the coefficients C ij vary with the position according to the following formulas [17] Equation (9) may be rewritten in the short form as

Finite Element Formulation
The expression of the strain energy of the beam is The expression of the variation of strain energy can be calculated as follows After integrating Equation (14) over the beam section and rewriting it in the matrix form, the variation of the strain energy can be computed as where R and T 01 are given by Gγ 01 dz (16) and where As consequence, Equation (16) can be rewritten as where In the current work, a two-node beam element is considered, each node includes five degrees of freedom. The vector of displacement of node i-th is The nodal displacement vector of the proposed beam element, U, is defined by The isoparametric geometry on two nodes of new beam element and the nodal variables are given by For the mixed finite element formulation, the authors approve a quadratic interpolation for α, β with parameters α m , β m and a constant shear resultant T 01 , which will be eliminated later.
in which, the shape functions N 1 , N 2 and N m are defined as follows where The first equation in Equation (8) is imposed in integral form as following Substitute α and w from Equations (24) and (25) into Equation (29), the parameter α m can be deduced as where Substitute Equations (24) and (30), into Equation (17), the strain variation vector, δω, can be obtained as Using Equation (30), the variable α becomes Then Using Equations (33) and (35), the strain variation Equation (32) can be expressed as follows where The shear strain variation δγ 01 can be obtained as Rewrite Equation (39) into the matrix form as below where Using Equation (25) the strain fields δε, δγ 01 become The variation of strain energy Equation (15) becomes Substituting Rewritten Equation (46) in the matrix form where K uu , K uT , K uβ , K TT , K Tβ , K ββ , f u , f T and f β are calculated as follows The parameter T 0 and rotation β m are then removed in two steps as below.
Step 1: eliminate rotation parameter β m where K 11 , K 12 , K 22 , f 1 and f 2 are expressed as Step 2: eliminate the stress resultant constant T 0 Substitute Equation (65) into Equation (59), the expression of the variation of strain energy is manifested as The stiffness matrix of new beam element K is now calculated as The nodal load vector of the beam element is expressed as

Convergence Study
In this subsection, some examples are given to determine the convergence rate of the proposed beam element. A homogeneous beam with the width b = 1 m, Young's modulus E = 29000 Pa, Poisson's ratio ν = 0.3 as in the work of Heyliger [40] is considered here.
Firstly, a cantilever beam subject to a concentrated load P = 100 N at the free end side is considered. Secondly, a simply supported beam under a uniform load q = 10 N/m is investigated. The numerical results for some cases of L/h ratios and number of elements are shown in Tables 1 and 2. The numerical results are compared with results of Heyliger [40] using two-node beam C 1 continuous formulation element based on HSDT. The comparison shows that the new beam element has an excellent convergence rate. Although the new proposed beam element uses linear shape functions, it has a better convergence rate than that of the beam element using higher-order shape function, thus it costs less effort and time of computation.

Validation Study
Continuously, to confirm the accuracy of the proposed beam element, the static bending of an FG beam subjected to a uniform load q is investigated. The FG beams made of two components, which are Aluminum and Alumina. The material properties of Aluminum and Alumina are E m = 70 GPa, E c = 380 GPa, ν m = 0.3, ν c = 0.3. Two cases of slenderness ratios L/h = 5 and L/h = 20 of the FG beams are considered. The displacements and stresses are calculated in the normalized form as.
For simply-simply (SS) and clamped-clamped (CC) supported beams For clamped-free (CF) supported beams For axial, normal and shear stresses The numerical results of dimensionless vertical displacement, normal stress and axial stress of FG beam using proposed beam are compared with the results of Li et al. [16] and Vo et al. [20] using analytical and finite element methods, which are given in Tables 3-7. The results in these tables show that the solutions from the present theory are very close with the results from HSDT of Li et al. [16] and Quasi-3D solutions of Vo [20] for different values of the power-law index, aspect ratio, and boundary conditions.   [16] and Vo [20].    The distributions of shear stress and axial stress along the depth of FG beam are compared with results of Li et al. [16] and Vo et al. [20] as in Figures 3 and 4. According to Figure 3, the shear stress distribution is parabolic along with the thickness and asymmetric for the FG beams. From Figure 4, the axial stress variation is not linear across the thickness of the FG beam, and its variation is linear across the thickness for isotropic (full ceramic or full metal) beams only. In general, the values of the axial stress do not equal to zeros at the mid-plane of the FG beams. Both shear stress variation and axial stress variation through the thickness of the FG beam using the proposed beam element are in remarkable agreement with those of Vo [20] using Quasi-3D theory.
According to the comparison, the results of the proposed beam element are very close actual adjacent to the Li et al. [16] and Vo et al. [20] solutions. Therefore, the new beam element can be applied to analyze FG beams.

Static Behaviour of FG Beams
In this subsection, an FG beam which is produced of Aluminum and Zirconium dioxide (Al/ZrO 2 ) under uniform distribution load q is investigated using proposed beam element. Various power-law indexes, slenderness ratios and boundary conditions are considered. The material properties of Al are E m = 70 GPa, ν m = 0.3, and the material properties of ZrO 2 are E c = 200 GPa, ν c = 0.3. The non-dimensional formulas are applied as follows In this study, some cases of boundary conditions are considered. For the simply-simply supported (SS) beam: u = w = 0 at x = 0, L; For the clamped-clamped supported (CC) beam: u = w = α = β = 0 at x = 0, L; For the clamped-simple supported (CS) beam: u = w = α = β = 0 at x = 0, and u = w = 0 at x = L; For the clamped-free supported (CF) beam: u = w = α = β = 0 at x = 0. The numerical results for bending behaviors of FG beams under uniform load are shown in Tables 8-11 and Figures 5 and 6. Table 8 and Figure 5 shows the nondimensional maximum vertical displacement of FG beams for some cases of boundary conditions, power-law index and different values of the length-to-height ratio. To show more clearly the effect of the power-law index and slenderness ratio, Figure 6 shows the dependence of nondimensional maximum vertical displacement of FG beams on the continuous transformation of the power-law index and slenderness ratio. It shows that the deflection of FG beams increases when increasing the power-law index. This can be explained that increasing the value of the power-law index leads to an increase in the component of metal in FGM, so that the FG beam becomes more flexible. Furthermore, it can be observed that the nondimensional maximum vertical displacement depends not only on the power law index, slenderness ratio but also boundary conditions, which is more pronounced for CC and CS beams than SS and CF beams Furthermore, it can be observed that the nondimensional maximum vertical displacement depends on power-law distribution index, boundary conditions and the length-to-thickness ratio. In addition, boundary conditions have more strongly effects on the deflection of CC and CS beams than those of SS and CF beams.   Table 9 shows the nondimensional axial stress  of FG beams subjected to a uniform load depends on some parameters and boundary conditions. Table 7 and Figure 7 present the distributions of nondimensional axial stress along with the depth of the SS and CC FG beams for different values of the power-law distribution index. The most significant aspect of this figure is that the axial stress distribution of FG beams is much more different from those of isotropic beams. As seen from Table 7 and Figure 7, the axial stress variation is not linear along with the thickness of the FG beams, the tensile stresses at the top are maximum. The values of the axial stresses do not equal to zeros at the mid-plane of the FG beams. This indicates that the neutral plane of the FG beams does not appear at the mid-plane, it is near the top face of the FG beams. This can be explained by the variation of the modulus of elasticity across the depth of the FG beams.   Table 9 shows the nondimensional axial stress σ * x (L/2, h/2) of FG beams subjected to a uniform load depends on some parameters and boundary conditions. Table 7 and Figure 7 present the distributions of nondimensional axial stress along with the depth of the SS and CC FG beams for different values of the power-law distribution index. The most significant aspect of this figure is that the axial stress distribution of FG beams is much more different from those of isotropic beams. As seen from Table 7 and Figure 7, the axial stress variation is not linear along with the thickness of the FG beams, the tensile stresses at the top are maximum. The values of the axial stresses do not equal to zeros at the mid-plane of the FG beams. This indicates that the neutral plane of the FG beams does not appear at the mid-plane, it is near the top face of the FG beams. This can be explained by the variation of the modulus of elasticity across the depth of the FG beams.   The non-dimensional shear stress distributions across the thickness of the beams made of FGM with different values of the power-law distribution index and some cases of boundary conditions are presented in Figure 8. The shear stresses of the full ceramic (isotropic) beams are symmetric about the mid-plane of the beams. In addition, the shear stress distributions are greatly influenced by the power-law index. In addition, Figure 8 shows the great dependence of the shear stress distribution on the power-law index. The non-dimensional normal stresses of the FG beams under uniform distribution load are shown in Table 11, which highlight the effect of thickness stretching on bending behaviors of FG beam from Quasi-3D theory. Due to the thickness stretching effect, the vertical displacement obtained from present Quasi-3D theory is smaller than those of HSDT and FSDT. The variation of the vertical displacement through the thickness of the FG beam for SS and CC boundary conditions are shown in Figure 9. According to Figure 9, the difference among the present Quasi-3D theory and other HSDT or FSDT is meaningful for thickness stretching. In this present The non-dimensional normal stresses of the FG beams under uniform distribution load are shown in Table 11, which highlight the effect of thickness stretching on bending behaviors of FG beam from Quasi-3D theory. Due to the thickness stretching effect, the vertical displacement obtained from present Quasi-3D theory is smaller than those of HSDT and FSDT.
The variation of the vertical displacement through the thickness of the FG beam for SS and CC boundary conditions are shown in Figure 9. According to Figure 9, the difference among the present Quasi-3D theory and other HSDT or FSDT is meaningful for thickness stretching. In this present Quasi-3D theory, the vertical displacement is not constant through the thickness of the beams as in HSDT and FSDT. Finally, to show more obviously the influence of normal deformation on the deflection of FG beams, we suggest the deflection ratio which is well-defined as the fraction of transverse displacement obtained by present Quasi-3D beam theory to that calculated by HSDT. The effect of normal deformation on the deflection of SS and CC supported FG beams is exhibited in Figure 10 for different values of power-law distribution index and slenderness ratio. Figure 10 shows that the deflection ratio is almost smaller than unity. It shows that the deflection will be decreased when the normal deformation effect is included. In the case of CC beams, there is a range of power-law index and slenderness ratio that causes the deflection ratio to be greater than unity, which represents that the normal deformation has more affect strongly than bending and shear deformation in this case. Finally, to show more obviously the influence of normal deformation on the deflection of FG beams, we suggest the deflection ratio which is well-defined as the fraction of transverse displacement obtained by present Quasi-3D beam theory to that calculated by HSDT. The effect of normal deformation on the deflection of SS and CC supported FG beams is exhibited in Figure 10 for different values of power-law distribution index and slenderness ratio. Figure 10 shows that the deflection ratio is almost smaller than unity. It shows that the deflection will be decreased when the normal deformation effect is included. In the case of CC beams, there is a range of power-law index and slenderness ratio that causes the deflection ratio to be greater than unity, which represents that the normal deformation has more affect strongly than bending and shear deformation in this case.
beams, we suggest the deflection ratio which is well-defined as the fraction of transverse displacement obtained by present Quasi-3D beam theory to that calculated by HSDT. The effect of normal deformation on the deflection of SS and CC supported FG beams is exhibited in Figure 10 for different values of power-law distribution index and slenderness ratio. Figure 10 shows that the deflection ratio is almost smaller than unity. It shows that the deflection will be decreased when the normal deformation effect is included. In the case of CC beams, there is a range of power-law index and slenderness ratio that causes the deflection ratio to be greater than unity, which represents that the normal deformation has more affect strongly than bending and shear deformation in this case.
(a) (b) Figure 10. The deflection ratio of FG beams subjected to uniform load, (a) SS beams and (b) CC beams.

Conclusions
In this paper, a new efficient Quasi-3D beam element was developed for static bending analysis of FG beams. Using mixed formulation, only 0 C continuous shape functions are required for finite element formulation of the new beam element. In addition, the new beam element presents the excellent results of displacement and stress even for a coarse mesh. Therefore, the proposed beam

Conclusions
In this paper, a new efficient Quasi-3D beam element was developed for static bending analysis of FG beams. Using mixed formulation, only C 0 continuous shape functions are required for finite element formulation of the new beam element. In addition, the new beam element presents the excellent results of displacement and stress even for a coarse mesh. Therefore, the proposed beam element costs less effort and time of computation than those using higher order shape functions, consequently, it can be widely applied for complex structural analysis. The shear stresses vary parabolically across the thickness of the FG beam, and equal to zeros at two free surfaces of beams, so it does not need any shear correction factors. The new beam element includes shear deformation and normal deformation. Effect of normal deformation is significant, and it should be considered in the static bending analysis of FG beams, especially for medium and very thick FG beams. The numerical results of the FG beams using the proposed beam element are in good agreement with other published results. The new beam element is accurate and efficient for bending behavior of FG beams. The influences of some parameters such as the power-law distribution index and length-to-thickness ratio are investigated.