Abstract
A corotational finite element formulation for large displacement analysis of planar functionally graded sandwich (FGSW) beam and frame structures is presented. The beams and frames are assumed to be formed from a metallic soft core and two symmetric functionally graded skin layers. The EulerBernoulli beam theory and von Kármán nonlinear straindisplacement relationship are adopted for the local strain. Exact solution of nonlinear equilibrium equations for a beam segment is employed to interpolate the displacement field for avoiding the membrane locking. An incrementaliterative procedure is used in combination with the arclength control method to compute the equilibrium paths. Numerical examples show that the proposed formulation is capable of evaluating accurately the large displacement response with just several elements. A parametric study is carried out to highlight the effect of the material distribution, the core thickness to height ratio on the large displacement behaviour of the FGSW beam, and frame structures.
1. Introduction
Large displacement analysis of structures is an important topic in the field of structural mechanics. This topic grows in importance due to the development of new materials which enable structures to undergo large deformation. Many investigations on the large displacement analysis of structures using both analytical and numerical methods are reported in the literature. Numerical methods, especially the finite element method with its versatility in spatial discretization, are an effective tool for the large displacement analysis of structures.
In order to analyse beam and frame structures undergoing large displacements by the finite element method, a nonlinear beam element which can model accurately the nonlinear behaviour of a structure is necessary to formulate. Various nonlinear beam elements for analysis of planar beam and frame structures are available in the literature. Depending on the choice of reference configuration, the nonlinear beam elements can be classified into three types: the total Lagrangian formulation, the updated Lagrangian formulation, and the corotational formulation. In the corotational formulation, which will be discussed herein, the kinematics are described in an element attached local coordinate system. The finite element formulation is firstly formulated in the local system and then transformed into a global system with the aid of transformation matrices. The elements proposed by Hsiao et al. [1, 2], Meek and Xue [3], Pacoste and Eriksson [4], and Nguyen [5] are some amongst the other corotational beam elements for large displacement analysis of planar beams and frames.
Analyses of structures made of functionally graded materials (FGMs) have been extensively carried out since these materials were created by Japanese scientists in 1984 [6]. The finite element analysis is often employed to handle the complexities arising from the material inhomogeneity. Chakraborty et al. [7] proposed a finite element formulation for studying the thermoelastic behaviour of shear deformable FGM beams. The formulation using the exact solution of an equilibrium FGM Timoshenko beam to interpolate the displacement field is free of shear locking. Based on the higherorder shear deformation beam theory, Kadoli et al. [8] investigated the static behaviour of FGM beams under ambient temperature. The Ritz method was employed by Lee et al. [9] in studying the postbuckling response of FGM plates subject to compressive and thermal loads. Almeida et al. [10] extended the total Lagrange formulation proposed by Pacoste and Eriksson in [4] to investigate the geometrically nonlinear behaviour of FGM beams. A finite element formulation based on the linear exact shape functions was employed by Taeprasartsit [11] for buckling analysis of perfect and imperfect FGM columns. Recently, the first author and his coworkers [12–15] derived the corotational finite element formulations for large displacement analysis of FGM beam and frame structures. In these works, in order to prevent the formulations from the membrane locking, the average strain has been introduced to replace the membrane strain.
Sandwich structures with height strengthtoweight ratio are widely used in aerospace application such as skin of wings, aileron, and spoilers. In order to improve the performance of the structures in high thermal environment, FGMs could be incorporated in the sandwich construction. Several analyses, mainly the vibration and buckling, of functionally graded sandwich (FGSW) beams have been carried out in recent years [16–18]. To maintain minimum weight for a given mechanical loading condition, FGSW beams and frames are often designed to be slender and they might undergo large displacements during service. Investigation on the large displacement behaviour of FGSW beams and frames, to the authors’ best knowledge, has not been reported in the literature, and it is studied in this paper for the first time. To this end, a finite element formulation based on EulerBernoulli beam theory is derived in the context of the corotational approach and used in the investigation. Different from previous works, the interpolation functions for the displacement field are obtained by solving the equilibrium nonlinear equations of a beam segment. The exact interpolation functions lead to a balance between the axial and transverse displacements, and the finite element formulation based on these functions is free of the membrane locking [19]. The concept of the average membrane strain mentioned above is therefore not necessary to use in the present work. Using the derived formulation, the equilibrium equations are constructed and they are solved by an incrementaliterative procedure in combination with the arclength control method. Numerical examples are given to show the accuracy of the derived formulation and to illustrate the effect of the material distribution and the core thickness ratio on the large displacement behaviour of the FGSW beam and frame structures.
2. Finite Element Formulation
2.1. FGSW Beam
Figure 1 shows a FGSW beam formed from a metallic soft core and two symmetric FGM skin layers in a Cartesian coordinate system. The axis is chosen to be on the midplane of the beam and the axis is directed upwards. Denoting the beam height and the soft core thickness by and , respectively, the FGM is assumed to be made of metal and ceramic with the volume fraction of ceramic, , varying in the thickness direction by the following powerlaw distribution:where is the nonnegative grading index. The volume fraction of metal is . As seen from (1), the top and bottom surfaces contain only ceramic, whereas the core is pure metal. The effective Young’s modulus of the beam, , evaluated by Voigt model readswhere and are Young’s moduli of ceramic and metal, respectively.
2.2. Corotational Approach
The corotational approach is a convenient way to derive geometrically nonlinear finite element formulations. The central idea of the approach is to introduce a local coordinate system that continuously moves and rotates with the element during its deformation process. The element formulation is firstly derived in the local coordinate system and then transferred into the global one by using transformation matrices. Depending on the definition of local coordinate system, various corotational beam formulations are available in the literature. The formulation adopted in the present work is closely related to the ones described in [4, 22], and its main points are summarized in the following.
Figure 2 shows a planar twonode beam element, , and its kinematics in two coordinate systems, a local and a global . The element is initially inclined to the axis an angle . The global system is fixed, while the local one continuously moves and rotates with the element during its deformation. The local system is chosen with its origin at node and with the axis directed towards node . In such a local system, the axial displacement at node and the transverse displacements at the two nodes and are always zero: . The element vector of local nodal displacements thus contains only three components:where is the local axial displacement at node and and are the local rotations at nodes and , respectively. In (3) and hereafter, a superscript “” denotes the transpose of a vector or a matrix, and the bar suffix is used to indicate a variable with respect to the local system.
The global nodal displacements in general are nonzero, and the global element vector of nodal displacements has six components aswhere , , are, respectively, the global axial and transverse displacements and rotation at node and , , are the corresponding quantities at node . The local displacement and rotations in (3) are related to the global ones byThe initial and current lengths, and , in (5) can be calculated aswith and being, respectively, the coordinates of the nodes and , and in (5) is the rigid rotation of the element which can be computed from the element coordinates and the current global nodal displacements [22].
Differentiating the local displacements in (5) with respect to the global displacements leads to the relation between the virtual vector of local displacements and that of global ones:where the transformation matrix has the following form:withEquating the internal virtual work for the element in both two systems leads to the relation between the global and local vectors of nodal forces aswhere and denote the local and global vectors of internal forces for the element, respectively.
The global element tangent stiffness matrix is obtained by differentiation of the global force vector with respect to the global nodal displacements, having the formIn the above equation, , , and are components of the vector , is the local tangent stiffness matrix, andDetail derivations of the matrix and vectors and are given in [22]. Equations (10) and (11) completely define the global nodal force vector and tangent stiffness matrix, provided that and are known.
2.3. Local Formulation
Based on EulerBernoulli beam theory, the displacements in the and directions at any point with respect to the local system, and , respectively, are given bywhere and are, respectively, the axial and transverse displacements of the point on the midplane and a suffix coma is used to denote the derivative with respect to the followed variable.
The von Kámán nonlinear relationship can be used for the local axial strain in large displacement analysis [13, 15]:Assuming linearly elastic behaviour for the beam material, the axial stress is related to the axial strain bywhere is the effective Young modulus, defined by (2) with replaced by .
The strain energy for an element with initial length of readswhere and are, respectively, the extensional and bending rigidities, and they are defined aswith being the crosssectional area. Substituting (2) into (17), one gets the explicit forms for and aswhere is the beam width. The interpolation functions for the displacements and in the present work are found by solving the homogeneous equilibrium equations of a beam element, which resulted from the strain energy (16) asThe element boundary conditions for the element are as follows:The solution in series form of (19) can be easily derived by using the command dsolve of Maple, and by keeping the thirdorder terms for , the solution has the following forms:where are constants which can be determined from element end condition (20). Finally, one can express the displacements and in terms of the local nodal displacements asThe finite element formulation based on the quintic function for the axial displacement and cubic polynomials for the transverse displacement does not encounter any membrane locking problem [19, 22]. The average strain, used by the first author and his coworkers [13, 15] for avoiding the membrane locking when a linear function is adopted for , is not necessary to use herein.
Substituting (22) into (16), one can express the strain energy in terms of the local nodal displacements , , and . The local nodal force vector and tangent stiffness matrix for the element are obtained by differentiating the strain energy with respect to the local nodal force vector as [4]Equation (23) together with (10) and (11) completely defines the finite element formulation of the element. A Maple code for derivation of the interpolation functions and the local element formulation as described above is given in the Appendix. The code also generates a Fortran code for the local internal force vector and tangent stiffness matrix.
3. Numerical Procedures
The derived element internal force vector and tangent stiffness matrix are assembled to construct nonlinear equilibrium equations for the structures, which can be written in the following form [22]:where and are the structural vectors of nodal displacements and nodal internal forces, respectively, is the fixed external loading vector, and the scalar is a load parameter. Vector in (24) is known as the residual force vector.
The system of (24) can be solved by an incrementaliterative procedure. The procedure results in a predictorcorrector algorithm, in which a new solution is firstly predicted from a previous converged solution and then successive corrections are added until some chosen convergence criterion is satisfied. In the present work, a convergence criterion based on Euclidean norm of the residual force vector is used for the iterative procedure aswhere is the tolerance, chosen by 10^{−4} for all numerical examples given in Section 4.
In order to deal with the limit point, the snapthrough, and snapback situations, at which the structure tangent stiffness matrix ceases to be positive definite, the arclength constraint method developed by Crisfield [23, 24] is adopted herewith. Numerical procedure in the present work is implemented by using the spherical arclength constraint method with the details described in [22].
4. Numerical Examples
Numerical examples are given in this section to show the accuracy and efficiency of the proposed formulation as well as illustrate the effect of the material distribution and the core thicknesstoheight ratio, , on the large displacement behaviour of the FGSW beams and frames.
4.1. Cantilever Beam under Tip Load
A FGSW beam composed of aluminum (Al) and zirconia (ZrO_{2}) subjected to a transverse load at its free end is considered. Young’s modulus of Al is 70 GPa and that of ZrO_{2} is 151 GPa [25]. An aspect ratio is assumed for the beam.
The validation of the derived formulation firstly needs to be confirmed. From the literature review, it is clear that there is no result available on the large displacement of FGSW beams and frames; the validation therefore is carried out on a pure metal cantilever beam. In Table 1, the normalized tip deflection of the isotropic cantilever beam obtained in the present work is compared to the analytical solution of Mattiasson [20] and the finite element results of Nguyen et al. [15] and Nanakorn and Vu [21]. In the table (and in the following also), denotes Young’s modulus of the metal. As seen from the table, the present formulation is more accurate than the two finite element formulations of [15, 21], which have been derived by using the corotational and total Lagrangian approaches, respectively. Thus, in addition to avoiding using the average strain, the exact interpolation adopted in the present work is capable of improving the accuracy also. The convergence of the proposed formulation is illustrated in Table 2, where the normalized tip displacements at various load amplitudes are given for different number of the elements. Irrespective of the load amplitude, the convergence is achieved by using just ten elements, which is very fast. In this regard, ten elements are used to discretize the cantilever beam in the computation reported in the following.
Figures 3 and 4, respectively, illustrate the effect of the grading index and the core thickness to the beam height, , on the large displacement response of the FGSW cantilever beam. At a given value of the applied load, as seen from Figure 3, the tip displacements increase as the grading index increases. The increase of the displacements can be explained by the fact that, as seen from (1), the beam associated with a higher index contains less ceramic. Since Young’s modulus of the ceramic is considerable higher than that of the metal, the rigidities of the beam with less ceramic percentage are smaller, and this leads to the larger displacements. The influence of the core thickness to the beam height ratio on the large displacement response of the beam, as seen from Figure 4, is similar to that of the grading index , and the tip displacements of the beam are increased by the increase of the ratio. This phenomenon can also be explained by the decrease of the rigidities of the beam with a higher ratio. Furthermore, the large displacement curves in Figures 3 and 4 gradually approach the curves of a homogeneous beam obtained in [15, 21] when the index grows to infinity or the core thickness approaches the beam thickness. This is reasonable since the rigidities of the FGSW beam gradually decrease when raising the index and the ratio, and as seen from (2), the beam is fully aluminum when or . It is worth mentioning that a load increment has been used in the numerical procedure in this example, and the maximum number of iterations is 8.
In Figure 5, the thickness distribution of the axial stress at the clamped end of the FGSW cantilever beam corresponding to a transverse load is illustrated for various values of the grading index and the core thicknesstoheight ratio. The effect of the material inhomogeneity and the core thicknesstoheight ratio on the stress distribution is clearly seen from the figure again. As seen from the figure, different from homogeneous and functionally graded beams, the curves for stress distribution of the FGSW beam are composed of three distinct parts, in which the stress in the two functionally graded layers is not straight due to the powerlaw variation of the effective modulus. The stress distribution is influenced by both the material inhomogeneity and the core thicknesstoheight ratio, and the maximum stress is increased by the increase of the index and the ratio. As the initial yield stress of aluminum is just N/m^{2}, the plastic deformation may be involved when the beam undergoes the large deformation. In order to take the effect of plastic deformation into account, an elastoplastic analysis should be employed instead of the elastic analysis used herein.
(a)
(b)
4.2. Asymmetric Frame
An asymmetric frame under a downward load as depicted in the righthand side of Figure 6 is analysed. The frame, as also known as Lee frame in the literature, is widely used by researchers to test the nonlinear finite element formulations and numerical algorithms because of its snapthrough and snapback behaviour [1, 4, 15]. The frame is formed from two FGSW beams with length cm, width cm, and height cm. The FGSW material is also assumed to be formed from aluminum and zirconia as the above cantilever beam. The frame is discretized by ten elements, five for each beam.
(a)
(b)
In Figure 6, the loaddisplacement curves for the frame are depicted for various values of the grading index . The transverse displacement and axial displacement in the figures have been computed at the loaded point. The loaddisplacement curves for the frame with different values of the core thicknesstoheight ratio are shown in Figure 7. The snapthrough and snapback behaviour of the FGSW frame considered herein is similar to that of the isotropic frame [1, 4]. The grading index remarkably affects the limit load of the frame, and the limit load decreases as the index increases. The core thickness to the beam height ratio, as seen from Figure 7, also changes the limit load of the frame, and the limit load decreases when the ratio increases. The effect of the grading index and the core thicknesstoheight ratio on the limit load of the frame can also be explained by the decrease of the frame rigidities as in case of the cantilever beam.
(a)
(b)
The relation between the applied load and axial stress at intersection point of the loaded line and interface surface of the soft core with lower functionally graded layer of the asymmetric frame is depicted in Figure 8 for various values of the index and the ratio. The stress at the point gradually grows when raising the applied load, and it then reaches a limit point. Again, the plastic deformation may have occurred during the frame undergoing the large deformation as the cases and in the figure, and as mentioned above, an elastoplastic analysis should be adopted to handle the effect of plastic deformation.
(a)
(b)
4.3. Portal Frame
In this last example, a FGSW portal frame made of steel and alumina under a download as depicted in Figure 9 is considered. Young’s modulus of steel is 210 GPa, and that of alumina is 390 GPa [25]. The beam is formed from three beams with cm, cm, and cm.
(a)
(b)
In Figure 9, the loaddisplacement curves of the frame are shown for various values of the index and for . The loaddisplacement curves of the frame with different values of the core to thickness ratio are depicted in Figure 10 for an index . Six elements, two for each beam, have been used in the analysis. The frame undergoes very large displacements before it reaches a limit point. The effect of the grading index and the core thickness to beam height on the large displacement response of the portal frame is similar to that of the asymmetric frame.
(a)
(b)
5. Conclusions
A nonlinear corotational finite element formulation for large displacement analysis of FGSW beam and frame structures has been derived. The FGSW beams and frames are assumed to be composed of a metallic soft core and two symmetric functionally grade layers. Based on EulerBernoulli beam theory, the nonlinear equilibrium equations for a beam element have been solved, and the obtained solution was employed to interpolate the displacement field. An incrementaliterative method was used in combination with the arclength control method to compute the large displacement response of the structures. Numerical results show that the convergence of the proposed formulation is fast, and the large displacement response of the structures can be accurately evaluated by just several elements. A parametric study has been carried out to highlight the effect of the material distribution and the core thickness to the beam height ratio on the large displacement behaviour of the FGSW beam and frame structures. As demonstrated in the numerical examples, the stress at some parts of the structures may exceed the yield stress, and it is important to take the effect of plastic deformation into consideration in the large displacement analysis of FGSW beams and frames. To this end, a nonlinear beam element, which is capable of modeling both the large displacement and elastoplastic behaviour of the FGSW beams and frames, is necessary to develop, and this work requires more efforts.
Appendix
This Appendix lists the Maple code necessary to derive the exact interpolation functions and to generate a Fortran code for the local internal force vector and tangent stiffness matrix described in Section 2.3 (see Algorithm 1).

Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under Grant no. 107.022015.02.