Abstract
The formulation of a family of advanced one-dimensional finite elements for the geometrically nonlinear static analysis of beam-like 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 higher-order beam elements accounting for shear deformation and local cross-sectional 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 Newton-Raphson 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 beam-like structures for different boundary conditions and slenderness ratios are carried out, showing that UF-based higher-order beam theories can lead to a more efficient prediction of the displacement and stress fields, when compared to two-dimensional 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 up-to-date 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 variational-asymptotical 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 two-dimensional bodies was presented by Wood and Zienkiewicz [6]. Newton-Raphson method was used in order to solve the nonlinear equilibrium equations. Surana [7] provided a geometrically nonlinear formulation for two-dimensional curved beams. A total Lagrangian approach was used and the beam element was derived using linear, paralinear, and cubic-linear plane stress elements. Dufva et al. [8] presented a two-dimensional shear-deformable 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 shear-deformable beams by using an absolute nodal coordinate formulation accounting for a nonrigid cross-sectional 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 beam-columns and frames. An optimum nonlinear solution technique within the Newton-Raphson 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 variational-asymptotic method. The geometrically nonlinear, three-dimensional elasticity problem was split into a linear, two-dimensional cross-sectional analysis and a nonlinear one-dimensional beam analysis. A novel two-dimensional 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 higher-order 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 two-dimensional model by using the technique of slowly variable Fourier coefficients. Garcia-Vallejo 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 beam-like structures via a hierarchical one-dimensional 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 higher-order beam theories is made formally general regardless of the through-the-thickness 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 thin-walled beam structures was demonstrated by using Lagrange polynomials as approximating functions for the cross-section kinematics within a layer-wise approach. An elastoplastic analysis via UF-based one-dimensional finite elements has been carried out by Carrera et al. [22], showing that such formulation can lead to a 3D-like accuracy in terms of displacements and stresses for compact and thin-walled 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 beam-like structures are carried out using UF-based finite elements and Taylor expansion of the displacement field along the thickness using an equivalent single-layer 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 cross-sectional 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 Newton-Raphson prediction/correction method. As far as the stress prediction is concerned, once the second Piola-Kirchoff 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 layer-wise approximation [20, 21], in both linear and nonlinear analyses, the use of Taylor polynomials allows an enrichment of the cross-sectional kinematics by simply increasing the order N and with no need for additional cross-sectional nodes. This feature makes Taylor-based refined models particularly suitable for the analysis of multilayered structures in the framework of an equivalent single-layer 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 cross-section () 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 two-dimensional 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 Green-Lagrange 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 Piola-Kirchhoff’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 Piola-Kirchoff 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 through-the-thickness 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 displacement-based 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, initial-displacement, 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 Newton-Raphson 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 low-order 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 cross-section is square with m. Slender () and short beams () are investigated. Cantilever, doubly clamped, and simply supported (hinged-hinged) 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 cross-section. 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 one-dimensional finite elements are compared with two-dimensional 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 one-dimensional corotational ANSYS finite elements “Beam3”, with both Euler-Bernoulli (EBT) and Timoshenko (TBT) kinematics, are considered for comparison and validation purposes. Results given by two-dimensional 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 one-dimensional 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 one-dimensional 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 lower-order shape functions are used, such as in linear and quadratic beam elements. It should be also noticed that, unlike a classical displacement-based 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 higher-order beam theories as well as the two-dimensional FEM results. On the other hand, if the slenderness ratio is reduced, the shear deformation effects as well as local cross-sectional 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 through-the-thickness 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 higher-order 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 higher-order 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
Clamped-clamped boundary conditions are considered. The evaluation point for is , whereas is evaluated at . Figures 11 and 12 show the load-displacement 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 higher-order 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 higher-order 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 load-displacement curves are presented in Figures 18 and 19, whereas Figures 20–22 show the through-the-thickness 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 .
(a) FEM 2D TL
(b) N=5
(a) FEM 2D TL
(b) N=5
(a) FEM 2D TL
(b) N=5
(a) FEM 2D TL
(b) N=5
(a) FEM 2D TL
(b) N=5
5. Conclusions
A family of refined one-dimensional finite elements derived through a Unified Formulation of the displacement field has been proposed for the geometrically nonlinear analysis of beam-like structures. Slender as well as short beams have been investigated for different boundary conditions. MITC method has been adopted in order to tackle the shear and membrane locking phenomena and improve the convergence performance of the proposed elements. The capability of UF-based structural theories to accurately yet efficiently predict both the displacement and the stress fields in the nonlinear regime via a one-dimensional model has been demonstrated in this work. As far as computational costs are concerned, the use of UF-based one-dimensional finite elements can save at least one order of magnitude in terms of DOFs, with respect to two-dimensional elements. The extension of the proposed formulation based on Taylor polynomials for an equivalent single-layer analysis of composite beam structures will be presented in a future work.
Appendix
Fundamental Nuclei of the Tangent Stiffness Matrix
The components of the linear stiffness matrix are The generic term is a cross-section moment:and it is a weighted sum (in the continuum) of each elemental cross-section area where the weight functions account for the spatial distribution of the geometry and the material. is the integral along the element axis of the product of the shape functions and/or their derivatives:These integrals are numerically evaluated through Gaussian quadrature. is the initial-displacement, or initial-slope, contribution to the tangent stiffness matrix. Its components are is the initial stress, or geometric, contribution to the tangent stiffness matrix. Its components are The integrals , , , and in (A.6) and (A.9) are given by If a MITC beam element is considered, the -integrals in (A.3), (A.11), and (A.13) are replaced, respectively, by the following integrals:
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work has been partially supported by the European Union within the Horizon 2020 Research and Innovation Program under Grant Agreement no. 642121.