# Geometrically Nonlinear Analysis of Beam Structures via Hierarchical One-Dimensional Finite Elements

Authors: Y. Hui; G. De Pietro; G. Giunta; S. Belouettar; H. Hu; E. Carrera; A. Pagani
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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 (l) is predominant with respect to any other dimension orthogonal to it. The cross-section (Ω=h×b) is defined by intersecting the beam with planes orthogonal to its axis. Figure 1 presents the beam geometry and the reference system, being h and b, respectively, the beam’s thickness and width. A fixed Cartesian reference system is adopted. The x coordinate is coincident with the axis of the beam and it is bounded such that 0≤x≤l, whereas the y- and z-axis are two orthogonal directions laying on Ω. The displacement field in a two-dimensional approach is(1)uTx,z=ux,zwx,zwhere u and w are the displacement components along the x- and z-axis, respectively. Superscript ‘T’ represents the transposition operator. For the sake of convenience, the displacements gradient vector θ is introduced:(2)θ=u,xu,zw,xw,zSubscripts ‘x’ and ‘z’, when preceded by comma, represent derivation versus the corresponding spatial coordinate.Figure 1
Beam geometry and reference system.Geometrical nonlinearity is accounted for in a Green-Lagrange sense. Large displacements and rotations are, therefore, considered. The strain vectorE is(3)Exx=u,x+12u,x2+w,x2Ezz=w,z+12u,z2+w,z2Exz=u,z+w,x+u,xu,z+w,xw,zEquation (3) can be written in the following matrix form:(4)E=H+12Aθθwhere(5)ET=ExxEzzExz(6)H=100000010110(7)Aθ=u,x0w,x00u,z0w,zu,zu,xw,zw,xA virtual variation of the strain vector cam be written as (see Crisfield [3])(8)δE=δH+12Aθθ=H+Aθδθwhere δ stands for the virtual variation operator.The vectorial form of second Piola-Kirchhoff’s stress tensorS is(9)ST=SxxSzzSxzThe material is supposed to withstand small strains. Hooke’s law is, therefore, considered:(10)S=QEIn the case of an anisotropic material, the reduced material stiffness matrix Q reads(11)Q=Q11Q13Q15Q13Q33Q35Q15Q35Q55Coefficients Qij 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 F and the Piola-Kirchoff stress tensor through the following relation:(12)σ=1JFSxxSxzSxzSzzFTwhere(13)F=1+u,xu,zw,x1+w,zand J = det(F) is the determinant of F. The weak form of the governing equations is obtained by means of the Principle of Virtual Displacement:(14)δL=δLint-δLext=0where L is the total work and Lint the internal one:(15)δLint=∫V0δETSdVV0 is the volume of the reference undeformed configuration. Lext is the work done by the external forces. An infinitesimal variation of the total work reads(16)dδL=∫V0δETdS+dδETSdVAfter few manipulations (see Crisfield [3]), (16) can be rewritten in the following form:(17)dδL=∫V0δETQdE+δθTS^dθdVwhere S^∈R4×4 is(18)S^=SxxSxz00SxzSzz0000SxxSxz00SxzSzzThe variation of the virtual work is finally written in terms of the actual and virtual variation of the gradient vector:(19)dδL=∫V0δθTHT+ATθQH+Aθ+S^dθdV


## 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 functionsFσz and along the axis by the shape functions Njx:(20)ux,z=FσzNjxqσjuwx,z=FσzNjxqσjwwithσ=1,2,…,Nu,j=1,2,…,Nnewhere qσjn:n=u,w are the nodal unknowns. Einstein’s compact notation is used in (20): a repeated index implicitly implies summation over its variation range. Nu is the number of terms accounted in the through-the-thickness expansion and it is arbitrary. Index j varies over the element number of nodes Nne 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 C0 sense up to an order Nn-1. 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 Fσz. The generic explicit form of the displacement field expanded via N-order Taylor polynomials is given by(21)ux=ux1+ux2z+ux3z2+⋯+uxN+1zNuz=uz1+uz2z+uz3z2+⋯+uzN+1zNN 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 reads(22)θ=FσNj,xqσjuFσ,zNjqσjuFσNj,xqσjwFσ,zNjqσjw=Gσjqσjwhere Gσj∈R4×2 and qσj∈R2×1 are(23)Gσj=FσNj,x0Fσ,zNj00FσNj,x0Fσ,zNjand(24)qσjT=qσjuqσjw


### 3.2. Tangent Stiffness Matrix

Once (22) is replaced within (19), the variation of the total virtual work reads(25)dδLe=δqτiT∫V0eGτiTHT+ATθQH+Aθ+S^GσjdVdqσj=δqτiTKτσijel+Kτσijet1+Kτσijet2dqσjwhere the superscript ‘e’ refers to the considered element, V0e=le·be·he is the element volume at the reference unstrained configuration, and KτσijelKτσijet1Kτσijet2∈R2×2 are the “fundamental nuclei” of the linear, initial-displacement, and geometric contributions to the tangent stiffness matrix:(26)Kτσijel=∫V0eGτiTHTQHGσjdVKτσijet1=∫V0eGτiTHTQA+ATQH+AGσjdVKτσijet2=∫V0eGτiTS^GσjdVThe nuclei are very general regardless of the approximation order N over the thickness, the class of approximating functions Fτ, and the number of nodes per element Nne 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:(27)E¯xx=N¯pExxpE¯zz=N¯pEzzpE¯xz=N¯pExzpwhere p denotes an implicit summation and varies from 1 to Nne-1. Exxp, Ezzp, and Exzp are the strain components coming from the geometrical relations in (3) evaluated at the p-th tying point rTp and N¯p are the assumed interpolating functions. Their expressions as functions of the natural beam element coordinate r∈-1,1 can be found in Carrera et al. [36] and they are reported below. For linear elements, the interpolation is reduced to a point evaluation, since(28)N¯1=1rT1=0For quadratic elements, the assumed interpolating functions and tying points are(29)N¯1=-123r-13N¯2=123r+13rT1=-13rT2=13And for cubic elements(30)N¯1=56rr-35N¯2=-53r-35r+35N¯3=56rr+35rT1=-35rT2=0rT3=35


## 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 functionsFσz and along the axis by the shape functions Njx:(20)ux,z=FσzNjxqσjuwx,z=FσzNjxqσjwwithσ=1,2,…,Nu,j=1,2,…,Nnewhere qσjn:n=u,w are the nodal unknowns. Einstein’s compact notation is used in (20): a repeated index implicitly implies summation over its variation range. Nu is the number of terms accounted in the through-the-thickness expansion and it is arbitrary. Index j varies over the element number of nodes Nne 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 C0 sense up to an order Nn-1. 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 Fσz. The generic explicit form of the displacement field expanded via N-order Taylor polynomials is given by(21)ux=ux1+ux2z+ux3z2+⋯+uxN+1zNuz=uz1+uz2z+uz3z2+⋯+uzN+1zNN 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 reads(22)θ=FσNj,xqσjuFσ,zNjqσjuFσNj,xqσjwFσ,zNjqσjw=Gσjqσjwhere Gσj∈R4×2 and qσj∈R2×1 are(23)Gσj=FσNj,x0Fσ,zNj00FσNj,x0Fσ,zNjand(24)qσjT=qσjuqσjw


## 3.2. Tangent Stiffness Matrix

Once (22) is replaced within (19), the variation of the total virtual work reads(25)dδLe=δqτiT∫V0eGτiTHT+ATθQH+Aθ+S^GσjdVdqσj=δqτiTKτσijel+Kτσijet1+Kτσijet2dqσjwhere the superscript ‘e’ refers to the considered element, V0e=le·be·he is the element volume at the reference unstrained configuration, and KτσijelKτσijet1Kτσijet2∈R2×2 are the “fundamental nuclei” of the linear, initial-displacement, and geometric contributions to the tangent stiffness matrix:(26)Kτσijel=∫V0eGτiTHTQHGσjdVKτσijet1=∫V0eGτiTHTQA+ATQH+AGσjdVKτσijet2=∫V0eGτiTS^GσjdVThe nuclei are very general regardless of the approximation order N over the thickness, the class of approximating functions Fτ, and the number of nodes per element Nne 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:(27)E¯xx=N¯pExxpE¯zz=N¯pEzzpE¯xz=N¯pExzpwhere p denotes an implicit summation and varies from 1 to Nne-1. Exxp, Ezzp, and Exzp are the strain components coming from the geometrical relations in (3) evaluated at the p-th tying point rTp and N¯p are the assumed interpolating functions. Their expressions as functions of the natural beam element coordinate r∈-1,1 can be found in Carrera et al. [36] and they are reported below. For linear elements, the interpolation is reduced to a point evaluation, since(28)N¯1=1rT1=0For quadratic elements, the assumed interpolating functions and tying points are(29)N¯1=-123r-13N¯2=123r+13rT1=-13rT2=13And for cubic elements(30)N¯1=56rr-35N¯2=-53r-35r+35N¯3=56rr+35rT1=-35rT2=0rT3=35


## 4. Numerical Results

The beam support is0,l×-h/2,h/2×-b/2,b/2. The cross-section is square with h=b=1 m. Slender (l/h=100) and short beams (l/h=10) are investigated. Cantilever, doubly clamped, and simply supported (hinged-hinged) beams made of aluminium (E=75 GPa and ν=0.33) are considered. A concentrated load Pz is applied at x/l=1,z/h=0 for the cantilever case and at x/l=1/2,z/h=0 for doubly clamped and simply supported boundary conditions. A dimensionless load factor λ=Pzl2/EI is used, I 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 121 nodes and beam theory N=5, corresponding to 1.5·103 degrees of freedom (DOFs), whereas a mesh of 240×24 elements was used for 2D FEM solutions, corresponding to 3.6·104 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. Figure2 shows the effectiveness of the MITC B2 elements in predicting the normalised displacement u^z=uz/uzCubic for increasing slenderness ratios l/h, where uzCubic 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 N 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 uzCubic 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.Figure 2
Locking correction via MITC method for linear elements and different beam theories, doubly clamped beam,λ=2, u^z evaluated at (l/2,-h/2).Figure 3
Convergence analysis of the normalised displacementu^z versus the total number of nodes Nn via cubic finite elements B4 and MITC B4, simply supported beam (l/h=100 and λ=3), u^z evaluated at (l/2,h/2).


### 4.2. Cantilever Beam

A slender cantilever beam (l/h=100) is first considered in order to validate the model towards classical reference beam solutions. Dimensionless displacements ui~=ui/l, Cauchy stresses σij~=σij2I/Pzlh, and thickness z~=z/h are considered. Figures 4 and 5 show the evolution of the displacement components with the load parameter. The displacement u~x is evaluated at (l,h/2), whereas u~z at (l,-h/2). 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 N equal to 2 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 N≥2 can reduce the error given by TBT from 5.5% to 0.7%, 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 0.6% for N≥3 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.Table 1
Displacements for a short cantilever beam,λ=10.
10 × - u ~ x 10 × u ~ z  Plane183 6.1163 8.7165  Beam3 TBT 6.4377 8.2409 Beam3 EBT 6.3579 8.1661 EBT [1] 6.2652 8.1061  N B2 B3 B4 B2 B3 B4 5 6.1579 6.1585 6.1585 8.6817 8.6820 8.6820 4 6.1575 6.1580 6.1581 8.6809 8.6812 8.6812 3 6.1596 6.1601 6.1601 8.6760 8.6764 8.6764 2 6.1521 6.1527 6.1527 8.6671 8.6674 8.6674Table 2
Cauchy stresses evaluated atx=l/4 and z=-h/2 for a short cantilever beam, λ=5.20.
10 × σ ~ x x 10 × σ ~ z z 10 × σ ~ x z  Plane183 2.4166 1.3372 1.7978 FEM 2D TL 2.7403 1.5302 2.0480  N B2 B3 B4 B2 B3 B4 B2 B3 B4 5 2.5622 2.7426 2.7414 1.6556 1.5336 1.5291 2.0730 2.0496 2.0474 4 2.5632 2.7435 2.7424 1.6552 1.5332 1.5287 2.0730 2.0494 2.0472 3 2.5554 2.7355 2.7344 1.6468 1.5248 1.5203 2.0779 2.0546 2.0522 2 2.4216 2.5995 2.5984 1.7133 1.5896 1.5851 2.1354 2.1114 2.1090Figure 4
Load factorλ versus dimensionless displacement u~x for a slender cantilever beam.Figure 5
Load factorλ versus dimensionless displacement u~z for a slender cantilever beam.Figure 6
Load factorλ versus dimensionless displacement u~x for a short cantilever beam.Figure 7
Load factorλ versus dimensionless displacement u~z for a short cantilever beam.Figure 8
Dimensionless Cauchy stressσ~xx along the thickness coordinate at x=l/4 and λ=5.20 for a short cantilever beam.Figure 9
Dimensionless Cauchy stressσ~zz along the thickness coordinate at x=l/4 and λ=5.20 for a short cantilever beam.Figure 10
Dimensionless Cauchy stressσ~xz along the thickness coordinate at x=l/4 and λ=5.20 for a short cantilever beam.


### 4.3. Doubly Clamped Beam

Clamped-clamped boundary conditions are considered. The evaluation point foru~x is (l/4,h/2), whereas u~z is evaluated at (l/2,-h/2). 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 4.6% (TBT) to 0.3% (N≥3), 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 2-nd order beam theory about 27.8% in the worst case and the one given by theories with N≥3 lower than 0.8%.Table 3
Displacements for a short doubly-clamped beam,λ=15.89.
10 3 × - u ~ x 10 2 × u ~ z  Plane183 9.4375 7.1822 FEM 2D TL 8.8841 6.7892  Beam3 TBT 9.2102 7.1003 Beam3 EBT 9.5180 6.4598  N B2 B3 B4 B2 B3 B4 5 8.8863 8.8848 8.8848 6.7911 6.7918 6.7918 4 8.8857 8.8841 8.8841 6.7879 6.7887 6.7887 3 8.8829 8.8814 8.8814 6.7725 6.7733 6.7733 2 8.8476 8.8461 8.8460 6.7109 6.7117 6.7117Table 4
Cauchy stresses evaluated atx=l/4 and z=0 for a short doubly-clamped beam, λ=11.04.
1 0 3 × σ ~ x x 1 0 3 × σ ~ z z 1 0 2 × σ ~ x z  Plane183 8.6194 3.2658 1.1290 FEM 2D TL 9.8958 3.1957 1.1281  N B2 B3 B4 B2 B3 B4 B2 B3 B4 5 9.9050 9.8829 9.8987 3.1738 3.2083 3.1858 1.1210 1.1315 1.1250 4 9.9077 9.8858 9.9014 3.1859 3.2205 3.1979 1.1208 1.1312 1.1248 3 9.9050 9.8830 9.8987 3.1858 3.2200 3.1976 1.1210 1.1314 1.1249 2 10.4620 10.4390 10.4550 2.3066 2.3405 2.3182 0.8477 0.8581 0.8516Figure 11
Load factorλ versus dimensionless displacement u~x for a slender doubly clamped beam.Figure 12
Load factorλ versus dimensionless displacement u~z for a slender doubly clamped beam.Figure 13
Load factorλ versus dimensionless displacement u~x for a short doubly clamped beam.Figure 14
Load factorλ versus dimensionless displacement u~z for a short doubly clamped beam.Figure 15
Dimensionless Cauchy stressσ~xx along the thickness coordinate at x=l/4 and λ=11.04 for a short doubly clamped beam.Figure 16
Dimensionless Cauchy stressσ~zz along the thickness coordinate at x=l/4 and λ=11.04 for a short doubly clamped beam.Figure 17
Dimensionless Cauchy stressσ~xz along the thickness coordinate at x=l/4 and λ=11.04 for a short doubly clamped beam.


### 4.4. Simply Supported Beam

Simply supported beams are considered.u~x and u~z are evaluated at (0,-h/2) and (l/2,h/2), 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 λ=6.03. The same considerations of the previous sections also apply to this case: the displacement prediction can be improved from an error of 3.7% for a 2nd order theory to 0.8% for a 5th order theory, as shown in Table 5. Similarly for the stresses given in Table 6, N=2 and B4 elements lead to a relative difference with respect to “FEM 2D TL” of about 2.3% for σ~xx, 60.5% for σ~zz, and 35.2% for σ~xz, whereas the errors given by a theory N=5 and B4 finite elements reduce to about 0.3%. 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 N=5 model and the results obtained by the 2D FEM solution are presented for different load factors λ.Table 5
Displacements for a short simply supported beam,λ=8.36.
1 0 2 × u ~ x 10 × u ~ z  Plane183 3.9585 1.3032 FEM 2D TL 3.8909 1.2805  N B2 B3 B4 B2 B3 B4 5 3.9239 3.9239 3.9192 1.2843 1.2843 1.2839 4 3.9011 3.9000 3.8960 1.2820 1.2819 1.2814 3 3.8267 3.8251 3.8232 1.2739 1.2738 1.2735 2 3.7501 3.7482 3.7474 1.2583 1.2580 1.2579Table 6
Cauchy stresses evaluated atx=l/4 and z=-h/2 for a short simply supported beam, λ=6.03.
1 0 2 × - σ ~ x x 1 0 3 × - σ ~ z z 1 0 2 × - σ ~ x z  Plane183 8.6470 5.2555 2.1314 FEM 2D TL 8.2064 4.8990 2.0056  N B2 B3 B4 B2 B3 B4 B2 B3 B4 5 8.2878 8.2030 8.2002 6.9946 4.8699 4.9120 2.2776 2.0010 2.0067 4 8.2696 8.1840 8.1811 6.9402 4.8190 4.8614 2.2707 1.9942 2.0002 3 8.2018 8.1142 8.1125 6.7966 4.6917 4.7352 2.2464 1.9717 1.9779 2 8.4861 8.3928 8.3917 4.0106 1.8893 1.9344 1.5683 1.2941 1.3006Figure 18
Load factorλ versus dimensionless displacement u~x for a short simply supported beam.Figure 19
Load factorλ versus dimensionless displacement u~z for a short simply supported beam.Figure 20
Dimensionless Cauchy stressσ~xx along the thickness coordinate at x=l/4 and λ=6.03 for a short simply supported beam.Figure 21
Dimensionless Cauchy stressσ~zz along the thickness coordinate at x=l/4 and λ=6.03 for a short simply supported beam.Figure 22
Dimensionless Cauchy stressσ~xz along the thickness coordinate at x=l/4 and λ=6.03 for a short simply supported beam.Figure 23
Through-the-thickness profile ofσ~xx for different load factors.Figure 24
Through-the-thickness profile ofσ~zz for different load factors.Figure 25
Through-the-thickness profile ofσ~xz for different load factors.Figure 26
Axial displacementu~x, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 27
Transverse displacementu~z, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 28
Piola-Kirchoff stressS~xx, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 29
Piola-Kirchoff stressS~zz, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 30
Piola-Kirchoff stressS~xz, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5


## 4.1. Locking Assessment

In order to correct the shear and membrane locking phenomena affecting nonlinear one-dimensional elements, MITC method was adopted. Figure2 shows the effectiveness of the MITC B2 elements in predicting the normalised displacement u^z=uz/uzCubic for increasing slenderness ratios l/h, where uzCubic 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 N 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 uzCubic 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.Figure 2
Locking correction via MITC method for linear elements and different beam theories, doubly clamped beam,λ=2, u^z evaluated at (l/2,-h/2).Figure 3
Convergence analysis of the normalised displacementu^z versus the total number of nodes Nn via cubic finite elements B4 and MITC B4, simply supported beam (l/h=100 and λ=3), u^z evaluated at (l/2,h/2).


## 4.2. Cantilever Beam

A slender cantilever beam (l/h=100) is first considered in order to validate the model towards classical reference beam solutions. Dimensionless displacements ui~=ui/l, Cauchy stresses σij~=σij2I/Pzlh, and thickness z~=z/h are considered. Figures 4 and 5 show the evolution of the displacement components with the load parameter. The displacement u~x is evaluated at (l,h/2), whereas u~z at (l,-h/2). 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 N equal to 2 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 N≥2 can reduce the error given by TBT from 5.5% to 0.7%, 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 0.6% for N≥3 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.Table 1
Displacements for a short cantilever beam,λ=10.
10 × - u ~ x 10 × u ~ z  Plane183 6.1163 8.7165  Beam3 TBT 6.4377 8.2409 Beam3 EBT 6.3579 8.1661 EBT [1] 6.2652 8.1061  N B2 B3 B4 B2 B3 B4 5 6.1579 6.1585 6.1585 8.6817 8.6820 8.6820 4 6.1575 6.1580 6.1581 8.6809 8.6812 8.6812 3 6.1596 6.1601 6.1601 8.6760 8.6764 8.6764 2 6.1521 6.1527 6.1527 8.6671 8.6674 8.6674Table 2
Cauchy stresses evaluated atx=l/4 and z=-h/2 for a short cantilever beam, λ=5.20.
10 × σ ~ x x 10 × σ ~ z z 10 × σ ~ x z  Plane183 2.4166 1.3372 1.7978 FEM 2D TL 2.7403 1.5302 2.0480  N B2 B3 B4 B2 B3 B4 B2 B3 B4 5 2.5622 2.7426 2.7414 1.6556 1.5336 1.5291 2.0730 2.0496 2.0474 4 2.5632 2.7435 2.7424 1.6552 1.5332 1.5287 2.0730 2.0494 2.0472 3 2.5554 2.7355 2.7344 1.6468 1.5248 1.5203 2.0779 2.0546 2.0522 2 2.4216 2.5995 2.5984 1.7133 1.5896 1.5851 2.1354 2.1114 2.1090Figure 4
Load factorλ versus dimensionless displacement u~x for a slender cantilever beam.Figure 5
Load factorλ versus dimensionless displacement u~z for a slender cantilever beam.Figure 6
Load factorλ versus dimensionless displacement u~x for a short cantilever beam.Figure 7
Load factorλ versus dimensionless displacement u~z for a short cantilever beam.Figure 8
Dimensionless Cauchy stressσ~xx along the thickness coordinate at x=l/4 and λ=5.20 for a short cantilever beam.Figure 9
Dimensionless Cauchy stressσ~zz along the thickness coordinate at x=l/4 and λ=5.20 for a short cantilever beam.Figure 10
Dimensionless Cauchy stressσ~xz along the thickness coordinate at x=l/4 and λ=5.20 for a short cantilever beam.


## 4.3. Doubly Clamped Beam

Clamped-clamped boundary conditions are considered. The evaluation point foru~x is (l/4,h/2), whereas u~z is evaluated at (l/2,-h/2). 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 4.6% (TBT) to 0.3% (N≥3), 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 2-nd order beam theory about 27.8% in the worst case and the one given by theories with N≥3 lower than 0.8%.Table 3
Displacements for a short doubly-clamped beam,λ=15.89.
10 3 × - u ~ x 10 2 × u ~ z  Plane183 9.4375 7.1822 FEM 2D TL 8.8841 6.7892  Beam3 TBT 9.2102 7.1003 Beam3 EBT 9.5180 6.4598  N B2 B3 B4 B2 B3 B4 5 8.8863 8.8848 8.8848 6.7911 6.7918 6.7918 4 8.8857 8.8841 8.8841 6.7879 6.7887 6.7887 3 8.8829 8.8814 8.8814 6.7725 6.7733 6.7733 2 8.8476 8.8461 8.8460 6.7109 6.7117 6.7117Table 4
Cauchy stresses evaluated atx=l/4 and z=0 for a short doubly-clamped beam, λ=11.04.
1 0 3 × σ ~ x x 1 0 3 × σ ~ z z 1 0 2 × σ ~ x z  Plane183 8.6194 3.2658 1.1290 FEM 2D TL 9.8958 3.1957 1.1281  N B2 B3 B4 B2 B3 B4 B2 B3 B4 5 9.9050 9.8829 9.8987 3.1738 3.2083 3.1858 1.1210 1.1315 1.1250 4 9.9077 9.8858 9.9014 3.1859 3.2205 3.1979 1.1208 1.1312 1.1248 3 9.9050 9.8830 9.8987 3.1858 3.2200 3.1976 1.1210 1.1314 1.1249 2 10.4620 10.4390 10.4550 2.3066 2.3405 2.3182 0.8477 0.8581 0.8516Figure 11
Load factorλ versus dimensionless displacement u~x for a slender doubly clamped beam.Figure 12
Load factorλ versus dimensionless displacement u~z for a slender doubly clamped beam.Figure 13
Load factorλ versus dimensionless displacement u~x for a short doubly clamped beam.Figure 14
Load factorλ versus dimensionless displacement u~z for a short doubly clamped beam.Figure 15
Dimensionless Cauchy stressσ~xx along the thickness coordinate at x=l/4 and λ=11.04 for a short doubly clamped beam.Figure 16
Dimensionless Cauchy stressσ~zz along the thickness coordinate at x=l/4 and λ=11.04 for a short doubly clamped beam.Figure 17
Dimensionless Cauchy stressσ~xz along the thickness coordinate at x=l/4 and λ=11.04 for a short doubly clamped beam.


## 4.4. Simply Supported Beam

Simply supported beams are considered.u~x and u~z are evaluated at (0,-h/2) and (l/2,h/2), 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 λ=6.03. The same considerations of the previous sections also apply to this case: the displacement prediction can be improved from an error of 3.7% for a 2nd order theory to 0.8% for a 5th order theory, as shown in Table 5. Similarly for the stresses given in Table 6, N=2 and B4 elements lead to a relative difference with respect to “FEM 2D TL” of about 2.3% for σ~xx, 60.5% for σ~zz, and 35.2% for σ~xz, whereas the errors given by a theory N=5 and B4 finite elements reduce to about 0.3%. 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 N=5 model and the results obtained by the 2D FEM solution are presented for different load factors λ.Table 5
Displacements for a short simply supported beam,λ=8.36.
1 0 2 × u ~ x 10 × u ~ z  Plane183 3.9585 1.3032 FEM 2D TL 3.8909 1.2805  N B2 B3 B4 B2 B3 B4 5 3.9239 3.9239 3.9192 1.2843 1.2843 1.2839 4 3.9011 3.9000 3.8960 1.2820 1.2819 1.2814 3 3.8267 3.8251 3.8232 1.2739 1.2738 1.2735 2 3.7501 3.7482 3.7474 1.2583 1.2580 1.2579Table 6
Cauchy stresses evaluated atx=l/4 and z=-h/2 for a short simply supported beam, λ=6.03.
1 0 2 × - σ ~ x x 1 0 3 × - σ ~ z z 1 0 2 × - σ ~ x z  Plane183 8.6470 5.2555 2.1314 FEM 2D TL 8.2064 4.8990 2.0056  N B2 B3 B4 B2 B3 B4 B2 B3 B4 5 8.2878 8.2030 8.2002 6.9946 4.8699 4.9120 2.2776 2.0010 2.0067 4 8.2696 8.1840 8.1811 6.9402 4.8190 4.8614 2.2707 1.9942 2.0002 3 8.2018 8.1142 8.1125 6.7966 4.6917 4.7352 2.2464 1.9717 1.9779 2 8.4861 8.3928 8.3917 4.0106 1.8893 1.9344 1.5683 1.2941 1.3006Figure 18
Load factorλ versus dimensionless displacement u~x for a short simply supported beam.Figure 19
Load factorλ versus dimensionless displacement u~z for a short simply supported beam.Figure 20
Dimensionless Cauchy stressσ~xx along the thickness coordinate at x=l/4 and λ=6.03 for a short simply supported beam.Figure 21
Dimensionless Cauchy stressσ~zz along the thickness coordinate at x=l/4 and λ=6.03 for a short simply supported beam.Figure 22
Dimensionless Cauchy stressσ~xz along the thickness coordinate at x=l/4 and λ=6.03 for a short simply supported beam.Figure 23
Through-the-thickness profile ofσ~xx for different load factors.Figure 24
Through-the-thickness profile ofσ~zz for different load factors.Figure 25
Through-the-thickness profile ofσ~xz for different load factors.Figure 26
Axial displacementu~x, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 27
Transverse displacementu~z, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 28
Piola-Kirchoff stressS~xx, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 29
Piola-Kirchoff stressS~zz, cantilever beam, λ=3.79, and l/b=10.
(a)
FEM 2D TL (b)
N=5Figure 30
Piola-Kirchoff stressS~xz, cantilever beam, λ=3.79, and l/b=10.
(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.



---

*Source: 4821385-2018-11-27.xml*