A Corotational Formulation for Large Displacement Analysis of Functionally Graded Sandwich Beam and Frame Structures

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 Euler-Bernoulli beam theory and von Kármán nonlinear strain-displacement 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 incremental-iterative 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.


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 higher-order 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][13][14][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 strength-to-weight 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][17][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 Euler-Bernoulli 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 arc-length 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.

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 ℎ 0 , 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 power-law distribution: where  is the nonnegative grading index.The volume fraction of metal is   = 1 −   .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 reads where   and   are Young's moduli of ceramic and metal, respectively.

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 two-node beam element, (, ), and its kinematics in two coordinate systems, a local (, ) and a global (, ).The element is initially inclined to the axis an angle  0 .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:   =   =   = 0.The element vector of local nodal displacements (d) 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 (d) has six components as where   ,   ,   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 by The initial and current lengths,  and   , in (5) can be calculated as with (  ,   ) 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 (T) has the following form: with Equating the internal virtual work for the element in both two systems leads to the relation between the global and local vectors of nodal forces as where f in and f in denote the local and global vectors of internal forces for the element, respectively.The global element tangent stiffness matrix k  is obtained by differentiation of the global force vector f in with respect to the global nodal displacements, having the form In the above equation,   ,   , and   are components of the vector f in , k  is the local tangent stiffness matrix, and Detail derivations of the matrix T and vectors r and z are given in [22].Equations ( 10) and ( 11) completely define the global nodal force vector and tangent stiffness matrix, provided that f in and k  are known.

Local Formulation.
Based on Euler-Bernoulli beam theory, the displacements in the  and  directions at any point with respect to the local system,  1 and  3 , respectively, are given by where () 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  by where () is the effective Young modulus, defined by (2) with  replaced by .
The strain energy for an element with initial length of  reads where  11 and  22 are, respectively, the extensional and bending rigidities, and they are defined as with  being the cross-sectional area.Substituting ( 2) into (17), one gets the explicit forms for  11 and  22 as where  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) as The 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 third-order terms for (), the solution has the following forms: where  1 ,  2 , . . .,  6 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 as The 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 f in and tangent stiffness matrix k  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.

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 p and q in are the structural vectors of nodal displacements and nodal internal forces, respectively, q ef is the fixed external loading vector, and the scalar  is a load parameter.
Vector g 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 as     g     ≤      q ef     , where  is the tolerance, chosen by 10 −4 for all numerical examples given in Section 4.
In order to deal with the limit point, the snap-through, and snap-back situations, at which the structure tangent stiffness matrix ceases to be positive definite, the arc-length constraint method developed by Crisfield [23,24] is adopted herewith.Numerical procedure in the present work is implemented by using the spherical arc-length constraint method with the details described in [22].

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 thickness-to-height ratio, ℎ 0 /ℎ, on the large displacement behaviour of the FGSW beams and frames.

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 /ℎ = 50 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, ℎ 0 /ℎ, 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 ℎ 0 /ℎ ratio.This phenomenon can also  be explained by the decrease of the rigidities of the beam with a higher ℎ 0 /ℎ 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 ℎ 0 approaches the beam thickness.This is reasonable since the rigidities of the FGSW beam gradually decrease when raising the index  and the ℎ 0 /ℎ ratio, and as seen from ( 2), the beam is fully aluminum when  = ∞ or ℎ 0 = ℎ.It is worth mentioning that a load increment Δ =   /2 2 has been used in the numerical procedure in this example, and the maximum number of iterations is 8.

Pure Al
In Figure 5, the thickness distribution of the axial stress at the clamped end of the FGSW cantilever beam corresponding to a transverse load  = 10  / 2 is illustrated for various values of the grading index and the core thickness-to-height ratio.The effect of the material inhomogeneity and the core thickness-to-height 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 power-law variation of the effective modulus.The stress distribution is influenced by both the material inhomogeneity and the core thickness-to-height ratio, and the maximum stress is increased by the increase of the index  and the ℎ 0 /ℎ ratio.As the initial yield stress of aluminum is just 2 × 10 9 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.

Asymmetric
Frame.An asymmetric frame under a downward load  as depicted in the right-hand 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 snap-through and snap-back behaviour [1,4,15].The frame is formed from two FGSW beams with length  = 120 cm, width  = 3 cm, and height ℎ = 2 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.
In Figure 6, the load-displacement 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 thickness-to-height ratio are shown in Figure 7.The snap-through and snap-back 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 ℎ 0 /ℎ ratio increases.The effect of the grading index  and the core thickness-to-height 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.
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 ℎ 0 /ℎ 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  = 3 and ℎ 0 /ℎ = 0.8 in the figure, and as mentioned above, an elastoplastic analysis should be adopted to handle the effect of plastic deformation.

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  = 120 cm,  = 3 cm, and ℎ = 2 cm.
In Figure 9, the load-displacement curves of the frame are shown for various values of the index  and for ℎ 0 /ℎ = 0.5.The load-displacement curves of the frame with different values of the core to thickness ratio are depicted in Figure 10 for an index  = 2. 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.

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  # local axial and traverse displacements u(x) and w(x) # local displacement and rotations ujL, tiL, tjL # local element vector feL, tangent stiffness matrix keL with(linalg); deq1flA11*(diff(u(x),x,x) + diff(w(x),x)*diff(w(x),x,x))=0; deq2flA11*(diff(u(x),x,x)*diff(w(x),x) + diff(u(x),x)*diff(w(x),x,x) + 3/2*(diff(w(x),x)) ∧ 2*diff(w(x),x,x)) + A22*diff(w(x),x$4)=0; dsolve(eq1,u(x)); deq2flA11*(diff(u(x),x,x)*diff(w(x),x) + diff(u(x),x)*diff(w(x),x,x) + 3/2*(diff((x),x)) ∧ 2*diff(w(x),x,x)) + A22*diff(w(x),x$4)=0; dsolve(eq2,w(x),series); w two symmetric functionally grade layers.Based on Euler-Bernoulli 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 incremental-iterative method was used in combination with the arc-length 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.

Figure 1 :
Figure 1: A FGSW beam in a Cartesian coordinate system.

Figure 2 :
Figure 2: A planar corotational beam element and its kinematics.

Figure 4 :
Figure 4: Effect of core thickness-to-height ratio on large displacement response of FGSW cantilever beam ( = 5).

Figure 8 :Figure 9 :
Figure 8: Relation between applied load and axial stress at intersection point of loaded line and interface surface of core and lower layer of asymmetric frame.

Table 2 :
Convergence of formulation in evaluating tip displacements of isotropic beam.
* Note: nELE is the number of elements.