ON COMPARISON OF 3D ISOGEOMETRIC TIMOSHENKO AND BERNOULLI BEAM FORMULATIONS

Abstract. Application of isogeometric analysis (IGA) for curved beams is very convenient for its ability of exact representation of curved geometries. Several beam formulation has been presented since the introduction of IGA. In this paper, two different beam formulations are presented: Bernoulli beam formulation of A. M. Bauer et al. [1], and Timoshenko beam element introduced by G. Zhang et al. [2]. Both beam elements are implemented and their performance is documented on the fully threedimensional example of helicoidal spring.


Introduction
Isogeometric analysis is an alternative of standard finite element analysis introduced in 2005 by T. Hughes [3]. The crucial impulse for its development was given by the desire for the automatic connection between Computer-aided design (CAD) and Finite element analysis. This goal has been achieved by the application of the same NURBS-based geometry description in CAD also for the unknown approximation in the analysis.
The use of CAD basis functions enables the exact description of arbitrarily curved geometries including curved beams. The focus of this paper is placed on the isogeometric analysis of spatial curved beams. Two different beam formulations are presented and their performance is compared on a benchmark problem. The first beam formulation is based on Bernoulli beam theory, while the Timoshenko kinematics is assumed in the second beam element presented.

NURBS-based analysis
Non-Uniform Rational B-Splines (NURBS) [4] are the most widespread technology used in CAD industry. Their fitness for CAD is rooted in the ability of the exact description of a vast range of geometries, including the conic sections widely used in architectural design. There are many effective algorithms developed for handling NURBS objects which can be also very useful for isogeometric analysis.
The computational domain in isogeometric analysis (IGA) is firstly divided into patches, which are further divided into knotspans (often referred to as elements in IGA), see Fig. 1. The NURBS curve is obtained by linear combination of Cartesian coordinates of the control points P and NURBS functions R p NURBS functions are rational polynomials which are generated by weighting B-spline functions with a weights associated with a control points. B-spline functions are piecewise polynomials and can be derived recursively starting with a piecewise constant functions, see [4] for the detailed description. The use of NURBS basis introduces some special aspects to the analysis. The NURBS basis functions are generally non-iterpolatory, thus the degrees of freedom in the control points lack the direct physical meaning. Additionally, the higher continuity between knotspans can be achieved naturally using NURBS, unlike C 0 in standard FEM. One of the main features of IGA is the fact, that the basis can be enriched without altering geometry representation using knot insertion and degree elevation algorithms.

Bernoulli beam element formulation
The first presented beam formulation is adopted from the work of A. M. Bauer et al. [1]. The element is based on Bernoulli beam theory, thus the orthogonality of the cross-section to the center line after deformation is assumed to be preserved and the crosssectional dimensions remain unchanged. The formulation accounts for the geometrically nonlinear behavior, nevertheless only the linear behavior is studied in the presented examples. The element has four degrees of freedom in each control point: three global displacements u, v, and w and an additional degree of freedom corresponding to the rotation around the centerline ψ.

Parametric space
Physical space x y Figure 1. Description of B-spline (NURBS) finite element geometry. The geometry is modeled by one patch consisting of two knotspans with quadratic NURBS approximation.
In the following, the standard notation using upper-case and lower-case letters for the undeformed and deformed configuration, respectively, is adopted. The geometry description is illustrated in Figure 2. The three-dimensional beam approximation can be reduced to the center line given by the position vector X c (x c ) whereX i are the control points coordinates and R p i are the corresponding basis functions. Deformed position vector is given aŝ whereû i is a vector of global degrees of freedom (u, v, w).

The cross-section orientation is described by a moving trihedral given by the base vectors
With assumption of Bernoulli kinematics, a position vector of an arbitrary point of a beam continuum given by coordinates X (x) can be expressed as (6) where θ i are the convective contravariant coordinates.
The first components of a moving trihedral A 1 and a 1 are aligned with a normalized tangents T and t, respectively, and are computed as a linear combination of the control points coordinates and the derivatives of the basis functions. The remaining components A α with α ∈ {2, 3} are orthogonal to the tangent of the center line (resp. A 1 ) and are derived from the reference trihedral, given by A 0 α , T 0 , using two key operations used for the alignment at the desired position. These operations are illustrated in Figure 3.
Analogical procedures denoted as Λ(T, t) and R t (ψ) are used to align the moving trihedral of the undeformed to the deformed configuration.
By substituting geometric relations into the definition of Green-Lagrange strain tensor, the individual components of the tensor can be derived as The diagonal terms E αα as well as the shear term E 23 are equal to zero. This yields from the Bernoulli assumptions (undeformable cross-section). For the off-diagonal terms E 12 , E 13 The relations between strains and displacements are not linear and to the purposes of the linear analysis considered in this paper need to be linearized. The energetically conjugated stress tensor to the Green-Lagrange strain tensor is the Second Piola-Kirchhoff tensor. The stiffness matrix is obtained by a standard procedure [5]. The fully derived equations for the presented beam formulation can be found in [1].

Timoshenko beam element formulation
The second beam element presented here has been proposed by G. Zhang et al. [2]. The beam formu- x y z Figure 2. Illustration of the Bernoulli beam in its undeformed and deformed configurations.  lation is based o the Timoshenko beam theory and is derived in local coordinate system (t, n, b), where t, n, b are tangent, normal and binormal vectors. There are six independent unknowns: tangential displacement u t , normal displacement u n , binormal displacement u b and rotations θ t , θ n , θ b (see Fig. 4). A strain-displacement matrix B is defined as where ε = {ε m , γ n , γ b , χ t , χ n , χ b } T and r = {u t , u n , u b , θ t , θ n , θ b } T . B is derived from the following relations for membrane strain ε m , shear strains γ n and γ b , torsional strain χ t , and bending strains χ n and χ b Curvature κ is given as and torsion τ where r(s) is the position vector and s is the curvilinear coordinate measured along the centerline of the beam. Conjugated generalized stress vector contains normal and shear forces, and torsional and bending moments. Similarly to the previous formulation, by following standard procedure the stiffness matrix is obtained.

Numerical locking
The well-known drawback of the Timoshenko beam element formulations is numerical locking. Locking phenomena rises from the use of the same approximation order for displacements and for rotations, which leads to the shear strains of a higher order than bending moments.
Several methods for unlocking beam elements are available,B-method [6] is used for the purposes of vol. 30/2021 Comparison of Timoshenko and Bernoulli beams this study. This method is based on the projection of the membrane and shear strains onto a basis of a lower order, leading to the removal of the locking. See [6] for the detailed description of the method.

Numerical examples
Both presented beam formulations have been implemented and their performance has been evaluated on the fully three-dimensional example of helicoidal spring cantilever subjected to the continuous constant force load [2] (see Figure 5). Three main aspects have been studied: numerical locking of Timoshenko beam, comparison of convergence of the Bernoulli and Timoshenko beam formulations, and error of the Bernoulli formulation caused by the neglect of the shear effects.

Numerical locking
Firstly, the removal of the numerical locking of the Timoshenko beam formulation has been verified. The presented Timoshenko formulation without and with locking treatment have been compared using the example of the thin helix depicted in Figure 5. The convergence of the tip vertical displacement has been observed. Additionally, the results have been compared to the standard FEM beam element with cubic approximation, which is naturally locking-free [5].
As expected, the numerical locking can be observed when no locking treatment is used. The convergence of the isogeometric beam is significantly deteriorated in comparison with the standard FEM beam element. WhenB-locking treatment is used a superior convergence is obtained. The slow convergence of the traditional FEM beam is due to the inability to exactly represent curved geometry. See Figure 6 for the results. On the other hand, isogeometric formulation with the locking treatment provides superior convergency.

Convergence comparison
In order to compare the Timoshenko and the Bernoulli beam formulations the error of the vertical tip displacement has been calculated as The reference solution have been calculated using fine mesh consisting of 320 control points individually for each formulation.
The results are shown in Figure 7. The better convergence is provided by the presented Bernoulli beam formulation, nevertheless, the Timoshenko beam element performs also very well.

Applicability of Bernoulli formulation
Finally, the applicability of the Bernoulli beam formulation has been assessed. The reference solution of the Timoshenko beam formulation has been used in order to verify the error caused by the neglect of the shear effects. The analysis has been performed for five different cross-sections in order to assess the behavior of beams with various thickness/length ratios.
The results can be observed in Table 1. As expected, the error rises significantly with the increasing radius of the cross-section.

Conclusions
Two isogeometric beam elements have been implemented: the Bernoulli beam element based on the work of A. M. Bauer et al. [1] and the Timoshenko beam element presented by G. Zhang et al. [2] withBmethod [6] used as a locking treatment. Three main aspects have been observed in order to compare the performance of the presented formulations.
Firstly, the numerical locking of the Timoshenko beam element has been documented and the suitability of the locking treatment used has been verified. Secondly, the convergence properties of both formulations have been analyzed, proving the better behavior of the Bernoulli formulation. Finally, the applicability of the Bernoulli beam element only to the thin structures has been demonstrated. The superior convergence compared to the traditional FEM formulation can be observed for both formulations.