Research Article  Open Access
Geometrically Nonlinear Analysis of Beam Structures via Hierarchical OneDimensional Finite Elements
Abstract
The formulation of a family of advanced onedimensional finite elements for the geometrically nonlinear static analysis of beamlike structures is presented in this paper. The kinematic field is axiomatically assumed along the thickness direction via a Unified Formulation (UF). The approximation order of the displacement field along the thickness is a free parameter that leads to several higherorder beam elements accounting for shear deformation and local crosssectional warping. The number of nodes per element is also a free parameter. The tangent stiffness matrix of the elements is obtained via the Principle of Virtual Displacements. A total Lagrangian approach is used and NewtonRaphson method is employed in order to solve the nonlinear governing equations. Locking phenomena are tackled by means of a Mixed Interpolation of Tensorial Components (MITC), which can also significantly enhance the convergence performance of the proposed elements. Numerical investigations for large displacements, large rotations, and small strains analysis of beamlike structures for different boundary conditions and slenderness ratios are carried out, showing that UFbased higherorder beam theories can lead to a more efficient prediction of the displacement and stress fields, when compared to twodimensional finite element solutions.
1. Introduction
Many structural elements, such as aircraft wings, rotor blades, robot arms, or structures in civil construction, can be idealised as beams. Furthermore, the hypothesis that the unstrained and deformed configurations are coincident at equilibrium is often not true and geometrical nonlinearities cannot be neglected. Engineering fields such as aeronautics, space, and automotive need more and more accurate models since an accurate prediction of the mechanics of beams plays a paramount role in their optimal design. Therefore, geometrically nonlinear modelling of beams represents an important and uptodate research topic.
A general overview on linear and nonlinear structural mechanics can be found in Nayfeh and Pai [2]. Nonlinear structural analysis via finite elements was thoroughly discussed in Crisfield [3] and Bathe [4]. Hodges et al. [5] provided a variationalasymptotical method that allowed obtaining an asymptotically correct strain energy for the approximation of stiffness coefficients for the prediction of geometrically nonlinear behaviour of composite beams. A paralinear isoparametric element for the geometrically nonlinear analysis of elastic twodimensional bodies was presented by Wood and Zienkiewicz [6]. NewtonRaphson method was used in order to solve the nonlinear equilibrium equations. Surana [7] provided a geometrically nonlinear formulation for twodimensional curved beams. A total Lagrangian approach was used and the beam element was derived using linear, paralinear, and cubiclinear plane stress elements. Dufva et al. [8] presented a twodimensional sheardeformable beam element for large deformation analyses. Cubic interpolation was used for the rotation angles caused by bending and linear interpolation polynomials were used for the shear deformations. The absolute nodal coordinate formulation was used for the finite element discretization along the beam axis. Further works on large rotations and large displacements analysis of sheardeformable beams by using an absolute nodal coordinate formulation accounting for a nonrigid crosssectional kinematics were carried out by Dufva et al. [9] and Omar and Sharana [10]. A geometric and material nonlinear analysis was carried out by Chan [11] for beamcolumns and frames. An optimum nonlinear solution technique within the NewtonRaphson scheme was obtained by minimizing the residual displacements. The evaluation of geometrically exact beam theories and models based on a second order approximation of finite rotations for the buckling and postbuckling analysis of beam structures was carried out by Ibrahimbegovic et al. [12]. Yu et al. [13] developed a generalised Vlasov theory for composite beams by means of the variationalasymptotic method. The geometrically nonlinear, threedimensional elasticity problem was split into a linear, twodimensional crosssectional analysis and a nonlinear onedimensional beam analysis. A novel twodimensional finite element solution for the nonlinear buckling and wrinkling of sandwich plates has been developed by Yu et al. [14]. Kirchoff’s theory was adopted for the kinematics of the skins, whereas a higherorder displacement field was considered for the core mechanics. The nonlinear governing equations were derived by the Principle of Virtual Displacements and solved via the Asymptotic Numerical Method (ANM). Based on this work, Huang et al. [15] proposed a more effective twodimensional model by using the technique of slowly variable Fourier coefficients. GarciaVallejo et al. [16] introduced a new absolute nodal coordinate finite element together with a reduced integration procedure in order to mitigate the locking phenomena in dynamic structural problems.
A static analysis of beamlike structures via a hierarchical onedimensional approach is addressed in this paper. The kinematics along the thickness is axiomatically assumed via a Unified Formulation (UF). UF had been previously proposed for plate and shell structures; see Carrera [17]; and it has been applied to the study of beams by Carrera et al. [18] and Carrera and Giunta [19]. Thanks to this approach, the derivation of a family of higherorder beam theories is made formally general regardless of the throughthethickness approximating functions and the approximation order. UF has been extended to large deflections and postbuckling analysis of beam structures in recent works by Pagani and Carrera [20, 21], where the capability of such approach to investigate global and local deformations in solid and thinwalled beam structures was demonstrated by using Lagrange polynomials as approximating functions for the crosssection kinematics within a layerwise approach. An elastoplastic analysis via UFbased onedimensional finite elements has been carried out by Carrera et al. [22], showing that such formulation can lead to a 3Dlike accuracy in terms of displacements and stresses for compact and thinwalled structures subjected to localized loadings. Applications of the UF approach for the investigation of multiphysics problems, vibration analyses, and static structural analyses of composite structures can be found in Carrera et al. [23], Biscani et al. [24], Koutsawa et al. [25], and Giunta et al. [26–29]. In this study, a large displacements, large rotations, and small strains analysis of beamlike structures are carried out using UFbased finite elements and Taylor expansion of the displacement field along the thickness using an equivalent singlelayer formulation. The approximation order is a free parameter and, therefore, several kinematic models can be straightforwardly obtained, accounting for nonclassical effects, such as transverse shear and local crosssectional warping. The number of nodes per element is also not a priori fixed. Linear, quadratic, and cubic elements are formulated. The tangent stiffness matrix of the element is derived from the weak form of the governing equations obtained via the Principle of Virtual Displacements (PVD). A total Lagrangian (TL) formulation is used and the global problem is solved by classical NewtonRaphson prediction/correction method. As far as the stress prediction is concerned, once the second PiolaKirchoff stresses have been obtained by the TL formulation, they are transformed into the true Cauchy stresses to have a direct comparison with results coming from updated Lagrangian formulations implemented in the commercial software ANSYS. As opposed to the Lagrange layerwise approximation [20, 21], in both linear and nonlinear analyses, the use of Taylor polynomials allows an enrichment of the crosssectional kinematics by simply increasing the order N and with no need for additional crosssectional nodes. This feature makes Taylorbased refined models particularly suitable for the analysis of multilayered structures in the framework of an equivalent singlelayer approach that will be presented in a future work. As further novelties with respect to [20, 21], the correction of shear and membrane locking phenomena in the nonlinear regime via the MITC method has been introduced, allowing an improved convergence performance of the proposed finite elements in slender structures. Finally, a detailed description of the stress analysis under large displacements, not often encountered in the literature, has been also provided, together with the discussion of some limitations of the proposed formulation with respect to finite elements with large strains capabilities.
2. Preliminaries
A beam is a structure whose axial extension () is predominant with respect to any other dimension orthogonal to it. The crosssection () is defined by intersecting the beam with planes orthogonal to its axis. Figure 1 presents the beam geometry and the reference system, being and , respectively, the beam’s thickness and width. A fixed Cartesian reference system is adopted. The coordinate is coincident with the axis of the beam and it is bounded such that , whereas the  and axis are two orthogonal directions laying on . The displacement field in a twodimensional approach iswhere and are the displacement components along the  and axis, respectively. Superscript ‘’ represents the transposition operator. For the sake of convenience, the displacements gradient vector is introduced:Subscripts ‘’ and ‘’, when preceded by comma, represent derivation versus the corresponding spatial coordinate.
Geometrical nonlinearity is accounted for in a GreenLagrange sense. Large displacements and rotations are, therefore, considered. The strain vector isEquation (3) can be written in the following matrix form:whereA virtual variation of the strain vector cam be written as (see Crisfield [3])where stands for the virtual variation operator.
The vectorial form of second PiolaKirchhoff’s stress tensor isThe material is supposed to withstand small strains. Hooke’s law is, therefore, considered:In the case of an anisotropic material, the reduced material stiffness matrix readsCoefficients are not reported here for the sake of brevity. They can be found in Reddy [30]. The Cauchy stress tensor can be derived from the deformation tensor and the PiolaKirchoff stress tensor through the following relation:whereand J = det() is the determinant of . The weak form of the governing equations is obtained by means of the Principle of Virtual Displacement:where is the total work and the internal one: is the volume of the reference undeformed configuration. is the work done by the external forces. An infinitesimal variation of the total work readsAfter few manipulations (see Crisfield [3]), (16) can be rewritten in the following form:where isThe variation of the virtual work is finally written in terms of the actual and virtual variation of the gradient vector:
3. Hierarchical Beam Elements
3.1. Kinematic and Finite Element Approximations
Within the assumed Unified Formulation and the finite element framework, the displacement components are approximated along the beam thickness via the base functions and along the axis by the shape functions :where are the nodal unknowns. Einstein’s compact notation is used in (20): a repeated index implicitly implies summation over its variation range. is the number of terms accounted in the throughthethickness expansion and it is arbitrary. Index varies over the element number of nodes and it is also a free parameter of the formulation. Linear, quadratic, and cubic elements along the beam axis are considered. These elements are addressed by “B2”, “B3”, and “B4”, respectively. The finite element shape functions approximate the displacements along the beam axis in a sense up to an order . For the sake of brevity, these functions are not reported here, but they can be found in Bathe [4]. Taylor polynomials are used as the class of expansion functions . The generic explicit form of the displacement field expanded via order Taylor polynomials is given by is the order of the approximating polynomials along the thickness and it is a free parameter of the formulation. By this approach, several displacementbased theories and finite elements accounting for nonclassical effects are straightforwardly derived. By replacing (20) within (2), the kinematic and finite element approximation of the displacements gradient vector readswhere and areand
3.2. Tangent Stiffness Matrix
Once (22) is replaced within (19), the variation of the total virtual work readswhere the superscript ‘e’ refers to the considered element, is the element volume at the reference unstrained configuration, and are the “fundamental nuclei” of the linear, initialdisplacement, and geometric contributions to the tangent stiffness matrix:The nuclei are very general regardless of the approximation order over the thickness, the class of approximating functions , and the number of nodes per element along the beam axis; see Carrera et al. [18]. Their explicit form can be found in the Appendix. Once the approximation order and the number of nodes per element are fixed, the element tangent stiffness matrix is obtained straightforwardly via summation of the previous nucleus corresponding to each term of the expansion. Finally, the nonlinear system is solved via the classical NewtonRaphson prediction/correction method.
3.3. Shear and Membrane Locking: MITC Beam Elements
In the geometrically nonlinear analysis of straight beams, the displacement components are coupled by the quadratic terms in the geometric relations; see (3). Therefore, membrane as well as shear locking phenomena will degrade the element performance and need to be mitigated, especially when slender structures and loworder shape functions are considered (see Reddy [31] and Malkus and Hughes [32] for more details). In this study, locking phenomena are overcome via the MITC method (see Bathe et al. [33–35]), consisting in the following interpolation of all the strain components along the beam element axis:where denotes an implicit summation and varies from 1 to . , , and are the strain components coming from the geometrical relations in (3) evaluated at the th tying point and are the assumed interpolating functions. Their expressions as functions of the natural beam element coordinate can be found in Carrera et al. [36] and they are reported below. For linear elements, the interpolation is reduced to a point evaluation, sinceFor quadratic elements, the assumed interpolating functions and tying points areAnd for cubic elements
4. Numerical Results
The beam support is . The crosssection is square with m. Slender () and short beams () are investigated. Cantilever, doubly clamped, and simply supported (hingedhinged) beams made of aluminium ( GPa and are considered. A concentrated load is applied at for the cantilever case and at for doubly clamped and simply supported boundary conditions. A dimensionless load factor is used, being the moment of inertia of the beam crosssection. Both displacement and stress values are given with respect to the initial fixed coordinate system.
Results for a plane stress, large displacements, large rotations, and small strains analysis provided by the proposed family of onedimensional finite elements are compared with twodimensional finite elements based on a total Lagrangian formulation and small strains hypothesis, referred to as “FEM 2D TL” (see Hu et al. [37]). Reference solutions from the available literature as well as classical onedimensional corotational ANSYS finite elements “Beam3”, with both EulerBernoulli (EBT) and Timoshenko (TBT) kinematics, are considered for comparison and validation purposes. Results given by twodimensional large strains ANSYS finite elements “Plane183” based on an updated Lagrangian formulation are also provided as a further assessment. About the computational costs, in order to be able to predict an accurate stress field in both short and slender beams, the most refined model used in the following numerical investigations for the proposed onedimensional formulation is given by a mesh of nodes and beam theory , corresponding to degrees of freedom (DOFs), whereas a mesh of elements was used for 2D FEM solutions, corresponding to DOFs. It should be noticed that the computational advantage coming from the UF approach is even more significant in nonlinear analyses when compared to linear analyses, since an iterative solution procedure is required and, therefore, a computational gain is obtained at every solution step.
4.1. Locking Assessment
In order to correct the shear and membrane locking phenomena affecting nonlinear onedimensional elements, MITC method was adopted. Figure 2 shows the effectiveness of the MITC B2 elements in predicting the normalised displacement for increasing slenderness ratios , where is the converged solution obtained with 40 B4 elements. It can be noticed that the locking correction strategy is effective regardless the beam theory order and the considered boundary conditions. Due to the presence of membrane locking, simply supported beams are the most critical case among those investigated, as far as element performance is concerned. Figure 3 shows that MITC correction in slender beams can significantly reduce the number of nodes needed for convergence. The converged reference solution is here obtained with 140 B4 elements. The improvement is even more significant when lowerorder shape functions are used, such as in linear and quadratic beam elements. It should be also noticed that, unlike a classical displacementbased finite element, an element adopting MITC correction strategy no longer assures a monotonic convergence “from below”, as shown in Laulusa et al. [38]. Following the previous convergence analyses, 121 nodes and MITC B4 elements are used for the plot results, whereas 121 nodes and MITC B2, B3, and B4 elements are compared in the table results at a fixed load parameter.
4.2. Cantilever Beam
A slender cantilever beam () is first considered in order to validate the model towards classical reference beam solutions. Dimensionless displacements , Cauchy stresses , and thickness are considered. Figures 4 and 5 show the evolution of the displacement components with the load parameter. The displacement is evaluated at , whereas at . For the case of slender beam, the reference solution based on EBT kinematics can accurately predict the nonlinear deformation and it matches the higherorder beam theories as well as the twodimensional FEM results. On the other hand, if the slenderness ratio is reduced, the shear deformation effects as well as local crosssectional warping become relevant and at least a beam theory with order equal to should be used for an accurate displacement prediction, as shown in Figures 6 and 7. A more detailed numerical comparison is given in Table 1, showing that beam theories with can reduce the error given by TBT from to , when compared to 2D FEM solutions. As far as stress prediction is concerned, Figures 8–10 show the throughthethickness profile of the stress components for the short beam case. Stress components are given in the initial fixed reference system. A small difference can be noticed between the large strains 2D FEM solution “Plane183” and the small strains “FEM 2D TL”. Furthermore, the higherorder beam theories match the 2D small strains solution, with relative differences being lower than for and B4 elements, as shown in Table 2. Unlike classical theories, the proposed higherorder models can predict global as well as localized displacements and stresses, even in the proximity of boundary conditions and throughout the nonlinear regime, as shown in Figures 26–30.


4.3. Doubly Clamped Beam
Clampedclamped boundary conditions are considered. The evaluation point for is , whereas is evaluated at . Figures 11 and 12 show the loaddisplacement curves for a slender beam, whereas the short beam case is shown in Figures 13 and 14. A significant difference between large and small strains hypothesis can be noticed in this latter case. Similarly to the cantilever case, Table 3 shows that higherorder beam theories can improve the accuracy with respect to the 2D FEM small strains solution from (TBT) to (), at worse. Figures 15–17 as well as Table 4 show that higherorder beam theories are required in order to accurately predict the stress profile in a thick beam, being the relative error given by a nd order beam theory about in the worst case and the one given by theories with lower than .


4.4. Simply Supported Beam
Simply supported beams are considered. and are evaluated at and , respectively. The loaddisplacement curves are presented in Figures 18 and 19, whereas Figures 20–22 show the throughthethickness profile of the stress components at the fixed load parameter . The same considerations of the previous sections also apply to this case: the displacement prediction can be improved from an error of for a nd order theory to for a th order theory, as shown in Table 5. Similarly for the stresses given in Table 6, and B4 elements lead to a relative difference with respect to “FEM 2D TL” of about for , for , and for , whereas the errors given by a theory and B4 finite elements reduce to about . The capability of the proposed formulation for an accurate stress prediction is preserved along the full deformation path, as shown in Figures 23–25, where the stress profile given by the model and the results obtained by the 2D FEM solution are presented for different load factors .

