Finite Element Approach of Bending and Buckling Analysis of FG Beams Based on Refined Zigzag Theory

This is the first attempt to utilize the refined zigzag theory (RZT) to study the bending and buckling behavior of a functional graded (metal/ceramic) thick beam. RZT, which has been exploited for the analysis of multilayered composite and sandwich beams does not employ shear correction factor. Furthermore, the number of kinematics variables of the RZT is not dependent to the number of layers in comparison with the layerwise theory. With regarding to the numerical solution, RZT, also requires C0 continuity interpolation, which leads to the development of this theory in finite element method. It is assumed that the mechanical properties of the beam varies through the thickness. According to the volume fraction of metal and ceramic, it is discretized across the thickness; consequently, the functionally graded beam (FGB) is modeled as a multi-layered one. The beam subjected to uniformly transverse and axial loadings. The equilibrium equations are established using the principle of virtual work. Using the shape functions of the first and second order forms, the non-isoparametric finite element consisting of three nodes and nine degrees of freedom are extracted. To confirm the excellent accuracy of the present approach, some numerical examples are provided and compared with those available in the literature that reveals that RZT is a trusty and validated theory to analyze the FG thick beams.


Introduction
Since multilayer composite structures are included of different mechanical properties through the thickness, depended to material consistent of each layer; there is a discontinuity in normal and transverse shear stresses at the interface between the layers. This leads to the separation of the layers in the borderline and consequently reducing the loading capability.
Nowadays, during development in aerospace technology, the requirement for the production of modern composite materials, which can withstand severe environmental phenomena including high thermal gradient, is perceived. Therefore, the idea of establishing the gradual variations in composition of new composites, from heat-resistant ceramics and metals with high machinery ability is formed [1][2]. This demands for new engineering materials has led to the production and development of new kinds of composite materials so called the functional grading material (FGM), which the mechanical and thermal properties change smoothly and continuously. The excellent features of the FGMs has overcome the undesirable and sudden changes of the mechanical and thermal properties in passing through a layer to the other one in laminated composite structures with gradual variations in the microscopic structure and compositional functional composition [3]. FGM structures have many advantages over multi-layer composite materials, for instance: low value of stress concentration, improved thermal resistance. Some applications of this kind of composite can be observed in nuclear reactors, membranes and catalysts in chemical industry, implanting artificial teeth, bones or artificial organs in medical engineering and worthy ability of thermal coating when structure subjected to high temperature environment.
Several studies have been devoted to provide various theories to anticipate the mechanical behaviour of FG structures. Swaminathans [4] presented a comprehensive review article to represent of the various methods employed to study the static, dynamic and stability behaviour of FG plates. Cheng and Batra [5] examined the buckling and steady state vibrations of a simply supported functionally gradient isotropic polygonal plate resting on a Winkler-Pasternak elastic foundation, using Reddy's third order shear deformation theory. Kapuria et al. [6] presented the third-order zigzag model to analyze the vibration and thermoelasticity of a FG beam, and demonstrated the 148 Finite Element Approach of Bending and Buckling Analysis of FG Beams Based on Refined Zigzag Theory benefits and confidence of proposed model on the FG beam of Al / SiC and Ni / Al2O3 using experimental tests. Sina et al. [7] presented a novel first-order shear deformation beam theory beam theory in order to analyze free vibration of FG beams. They ignored the lateral normal stresses and derived the governing differential equation of motion through Hamilton's principal. Farhatnia et al. [8] studied numerical and analytical solution to obtain the distribution of stresses in a three-layer beam consisting of metal-FGM-ceramic layers, using the Bernoulli Euler theory. The material properties is assumed to vary through the thickness following the power law distribution in terms of material constituents. In the other work, Farhatnia et al. [9] employed the first order shear deformation to analyze the thick FG beam. They discretized the differential equations and solved them using differential quadrature method. Hosseini et al. [10] studied the thermal buckling of FG beam using the theory of Timoshenko beam and Von-Karman nonlinearity hypothesis. They determined the critical buckling temperature using the Fourier series and Stokes' transformation.
As the analytical solutions are not applicable for all engineering problems, the numerical methods enable the researchers to solve many complicated problems; several studies have been devoted to the development of numerical methods, particularly the finite element method. Filippi et al. [11] presented a finite element solution based on the 1D Carrera Unified Formulation (CUF) to perform static analysis of FG structures. They derived the governing equations via weak form of Principle of Virtual Displacements. Khan et al. [12] presented one dimensional finite element model to consider static and vibration behaviour of FG beam. In order to define displacement field, they used linear interpolation and cubic Hermit interpolation functions for axial and lateral displacement, respectively. Kahya and Turan [13] suggested an element based on the first-order shear deformation theory to analyze the vibration and critical buckling load of a multi-layered FG sandwich beam. The proposed element having 3N + 7 degrees of freedom that N denotes the number of layers in the sandwich beam. Tessler et al. [14] and Gherlone et al. [15] developed the RTZ for multilayer composite and sandwich composite plates using kinematic assumptions of first order shear deformation theory. By introducing a C0-continuous Zigzag function, they provided a combination for the displacement field. The zigzag function of this theory is spontaneously zero on the surface and underneath of the multi-layer beam. Also, for a single-layer beam, it turns into beam Timoshenko's theory. Using this theory, it is not necessary to apply the shear correction factor and is possible to model the clamped edge condition Oñate et al. [16] presented a new simple linear two-noded beam element based on combination of Timoshenko theory and RTZ. The element has the capability to descend the order of integration to eliminate the shear locking effect. Moreover, using this element, they examined the delamination in a composite multi-layer beam. Discuva et al. [17] assessed a class of efficient higher-order continuous RZT-based beam elements. They adopted the interdependent/anisoparametric interpolations to have shear locking effect and to approximate the four independent kinematic variables.
Good accuracy and the convenient application of the RTZ, beside the benefit of avoiding of shear correction factor, as well as the independency of the degree of freedom of elements from the number of layers, evaluated this theory appropriate for analyzing of multi-Layer composite and sandwich structures using the finite element method. Hence, owing to the widespread usage of functionally graded materials in modern engineering structures with at least one layer of FG material, the assessment of exploiting RTZ is of great importance in the development of numerical solution for this kind of structure. Most notably, to the best of our knowledge, this is the first study, which investigates the three-layered FG beam by RZT. We modelled a FG beam as a multi-layered one by discretizing of the mechanical properties across the thickness direction and afterwards accomplished the buckling and bending analysis using RTZ.

Theory and Formulation
Consider a FGB with length of L and rectangular cross section of b×h, as shown in Fig. 1-a. The upper and lower surfaces are consisted of 100% ceramics and metal, respectively. Traditional functional graded materials, which are mixture of metal and ceramic, are non-homogeneous in microstructure and the mechanical properties are smoothly and continuously varied from one surface to another, as shown Fig. 1-b. It is assumed mechanical properties of FGB vary through the thickness direction with the power-law distribution. In this work, in order to employ the (RZT), the FGB is modeled as a N-layer with equal thicknesses, as shown in Fig. 1-c, which has the effective properties of the material in the k-dimensional layer. Young's modulus E and the shear modulus G is determined by (1) and (2) proportional to the position across the thickness. In the following, by discretizing the FGB as a multilayer beam, Young's modulus and shear modulus is determined in each layer [3].

Figure1. Geometry of a FG beam subjected to axial and transverse loadings
Wherein Vc denote volume fraction of ceramic constituents, subscripts m and c denote to metal and ceramics, respectively. The volumetric fraction of k th -layer is defined based on rule of the power law, as follows: The kinematic field in zigzag beam theory is generally written as [14]: In (4) the axial displacement ( ), the rotation ( ) and the transverse deflection ( ) are the primary kinematic variables of the underlying equivalent single-layer Timoshenko beam theory. ( ) ( ) denotes a piecewise linear zigzag function,. ( ) is a primary kinematic variable that defines the amplitude of the zigzag function along the length of the beam. In RZT, ( ) ( ), is defined as [14]: For a beam with N layers, there are N+1 zigzag displacement, ( ) , in the interfacial layers and outer surfaces of beam. For instance, as seen in Fig. 2, N equals to 3; therefore, (0) , (3) denotes for lower and upper surfaces, respectively and (1) , (2) indicates for zigzag displacements in the layer interfaces. The following notation is generally expressed for a three-layer beam in zigzag theory as shown in Fig. 2.
The displacement in interface of layers may be determined as follows [14]: By assumption small deformation and appropriate deformation of (4), the nonzero strain components are defined as: where ( ) , the derivative of zigzag function with respect to z, which is deduced from the (9c) [14]: It should be noted before, the shear correction factor is not applicable to calculate shear stress.

Principal of Virtual Work
By employing principal of virtual work, the Lagrangian governing equilibrium equations within boundary conditions can be extracted: where denotes variational operator. , indicate the strain energy and virtual work of external loading, respectively.
Upon substitution of (9) into (13) and integration over thickness, the variation of strain energy is obtained as follows: Where is stress couple, , are the normal and shear stress resultants, respectively.
As shown in Fig. 3, by considering the FG beam is subjected to axial and distributed transverse loading, the virtual work of external loadings, may be derived as follow: Where, The virtual work of axial load is determined, as follows: By substituting (14) and (18) and (19) into (12), the statement of principal of virtual work is obtained as follows:

Finite Element Solution
In order to obtain the stiffness matrix of the RZT-based elements, it is necessary to select an appropriate and consistent interpolation strategy and then applying the principle of virtual work [19]. The standard Timoshenko element ( 2 nodes and 6 degrees of freedom), which uses linear Lagrangian polynomials for u, w, θ, and the complete Guass integration to obtain the integral of the element matrix, gives a stiff solution for thin beams with curvature nearly equals to zero. This phenomenon is associated with the shear locking phenomena owing to the elimination of the curvature effect in the measurement of the shear strain when employing the linear isoparametric model in the thin beam. In order to prevent the elimination of the curvature effect , , it is necessary that the polynomial order of the lateral displacement w(x) is at least one order higher than the polynomial used for estimating θ(x). This approach is designated as non-isoparametric interpolation strategy, which implies the use of different interpolation polynomials in the approximation of θ(x), w(x). The appropriate element applicable to the aforementioned strategy is shown in Fig. 4a, c [15] As shown in Fig.4, a beam element with length e is consisted of three nodes. The coordinate of the beam with respect to the global and local coordinates are shown. In order to achieve characteristics of non-isoparametric beam element based on RZT, we use the simple interpolation functions for kinematic parameters, which according with the degrees of freedom for u, θ, ψ are the standard linear shape functions (21) and for w the standard shape functions of the order two is used (22), as follows [19]: By employing of aforementioned shape functions, the kinematics variables are defined as follows: By substituting (23) into (15), and consequently the resultants into (20) and carrying out some mathematical manipulation, the standard matrix form of the element-level equilibrium equations are obtained as follows: where [ ] , g e , { } , { } denote the stiffness and geometrical matrices, degrees of freedom vector, distributed transverse loading, respectively. The components of the aforementioned matrices are given in Appendix. Also, { } represent for forces and moments that concentrated in the nodes of the element. The critical buckling load and corresponding modes can be obtained by setting zero of { } ,{ } ; therefore for determination a non-trivial solution, an eigenvalue problem is solved by setting zero the determinant of coefficient. In this case, P and { } are eigenvalue (critical buckling load) and eigenvector (buckling modes), respectively.
In Table 1, the governing boundary conditions for a beam according with various end supports are represented as:

Numerical Results
Herein, the influence of boundary conditions, span-to-depth ratio, FG power law index on axial and 152 Finite Element Approach of Bending and Buckling Analysis of FG Beams Based on Refined Zigzag Theory lateral displacements, axial and shear stress distribution and critical buckling loads on a ceramic-metal FGB is considered. The combination of materials consists of Aluminum as metal and Alumina as ceramic. The mechanical properties are given in Table 2  In the present work, to generalize the following results, the dimensionless parameters and variables are defined as: The shear stress , ( ) , for each layer are obatined by integrating the following differential equilibrium equation (27): In equation (27), the axial stress ( ) for each layer, is obtained via the solution of finite element zigzag theory. In Fig. 5, the convergence of numerical procedure is illustrated for clamped-clamped beam with a power-law index of m = 0 and m = 5, L/h= 5 and 10, the convergence is measured by the absolute value of the relative error defined as: (28) where i denotes the number of elements to determine the lateral displacement. As depicted, the percentage of error is descended by increasing the number of elements, wherein the discretization the length of the beam to 12 elements, yields good accuracy (less than 0.4). In Table 3, the variation of the critical buckling load is listed for different boundary conditions; FG power index is taken in the range of 0-10. The present FEM results are compared with Euler Bernoulli theory (EBT) and finite element solution of refined shear deformation theory (RSDT). It is evident that the present results have good agreement and can be used as the benchmark in the design of the columns. As expected, when the FG power law index increases, that means the percentage of metal contents increases, the non-dimensional critical buckling load is reduced; whereas, variation of length-to thickness ratio does not affect the critical loads significantly. The maximum and minimum values of the critical buckling load are associated with the clamped-clamped beam and clamped-free one, respectively. As expected, the EBT overestimated the critical buckling loads compared with the aforementioned two theories. In Table 4, the displacements and stresses of simply supported FGB obtained subjected to the uniform distributed loading are presented for different values of the power law index and l/h= 5, 20. Also, in order to validate the present results, the comparison study is carried out. It is compared with the results of the two theories, namely, the Euler-Bernoulli theory [19] and the third shear deformation theory [19]. As expected, the CBT underestimates the displacements and stresses, whereas results of RZT and TBT are in close agreement.

Figure 5. A mesh convergence test for clamped -clamped FGB
The transverse displacement along the FGB length is depicted in Fig. 6 with l/h=5 for different values of the power law index, subjected to uniform transverse loading under different end supports. As expected, by increasing the ceramic volumetric fraction in FGB, the transverse displacement is reduced, due to decreasing of bending rigidity of the FG beam. The minimum amount of transverse displacement is devoted to clamped-clamped beam.
The distribution of axial displacement, normal stress and transverse strain of the FGB through the thickness of the beam are exhibited in Fig. 7 to Fig. 9 at position x =L/4, L/h=5.
In Fig. 7, the distribution of axial displacement through the thickness are sketched. In all cases, the distribution of axial displacement are smooth. In FGB, neutral axis is not coincident with the centroid of the beam cross section. The axial displacement increases nonlinearly away from the neutral axis until the maximum values at the top and bottom of the beam. The linear distribution of is devoted to pure ceramic or pure metallic beam.
As illustrated in Fig. 8, the boundary conditions have the considerable influence on the distribution of the axial stress, and except for a ceramic beam (m = 0), the maximum amount of axial stress occurs in the upper surface of the beam. Also, the axial stress variations for the FGB are not linear across the thickness; by increasing the volumetric percentage of ceramic constituents, the axial stress increases on the top of the FGB. It is be noted that for m>0, the neutral axis is not coincident with middle surface; it is placed above the middle surface of FG beam at ≅ 0.22. From bottom surface of beam toward the neutral axis, tensile stress occurs and become compressive afterwards. As shown, the effect of power law index (m=2,5,10) on location of neutral axis is not considerable.    As shown in Fig. 9, the quantity of the shear stress at the upper and lower surface of the beam equals to zero. The location of the maximum shear stress is dependent to the value of the power law index; therefore for a homogenous beam (m=0), it is located at h=0, whereas for m =2, 5, 10, the location moves toward the upper surface. As observed, by increasing the value of m, the maximum amount of shear stresses increases, though this variation is not significant for m>0. Among C-C, S-S and C-F beams, it is seen that the minimum amount of non-dimensional shear stress is devoted to C-F beam. Also, the discrepancy in values of ̅ between C-C and S-S beams is negligible.

Conclusions
In present work, the finite element approach of buckling and bending of FG (ceramic-metal) Timoshenko were considered by exploiting of RZT. For this aim, the FG beam was divided by N layers as multi-layer composite beam and employed the principal of virtual work based on three-noded non-isoparametric elements with nine degree of freedom. The influences of different boundary conditions, span-to-thickness ratio(L/h) and various values of FG power law index on bending and bucking of FGB are studied. The following conclusions are noteworthy: (1) Convergence of the FEM results was achieved rapidly with less than 12 elements of the three-noded and non-isoparametric elements. (2) The critical buckling load is increased only up to a length-to-thickness ratio of 15 and for the higher values, the critical buckling load remains constant. The amounts of deflection based on RZT is closer to the third order shear theory with L/h=5 (thick beam) and closer to those of the classical theory for L/h=20. (3) It is remarkably noting that RZT is a trustworthy procedure to examine the static behavior of FG beams, accurately. Future works would concern on investigation on dynamic and static analysis of the multi-layered FGB based on higher-order shear deformation theory.

Appendix
The submatrix can be obtained as: