DOI:10.14311/APP.2021.30.0012 Acta Polytechnica CTU Proceedings 30:12–17, 2021 © Czech Technical University in Prague, 2021 available online at http://ojs.cvut.cz/ojs/index.php/app ON COMPARISON OF 3D ISOGEOMETRIC TIMOSHENKO AND BERNOULLI BEAM FORMULATIONS Edita Dvořáková∗, Bořek Patzák Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: edita.dvorakova@fsv.cvut.cz 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 three- dimensional example of helicoidal spring. Keywords: Bernoulli theory, curved beam element, isogeometric analysis, NURBS, Timoshenko theory. 1. Introduction Isogeometric analysis is an alternative of stan- dard finite element analysis introduced in 2005 by T. Hughes [3]. The crucial impulse for its develop- ment was given by the desire for the automatic con- nection between Computer-aided design (CAD) and Finite element analysis. This goal has been achieved by the application of the same NURBS-based geome- try description in CAD also for the unknown approx- imation in the analysis. The use of CAD basis functions enables the exact description of arbitrarily curved geometries includ- ing 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 prob- lem. The first beam formulation is based on Bernoulli beam theory, while the Timoshenko kinematics is as- sumed in the second beam element presented. 2. 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, in- cluding the conic sections widely used in architectural design. There are many effective algorithms devel- oped 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 Rpi C(ξ) = n! i=1 R p i (ξ)Pi. (1) 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 de- rived 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. 3. 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 orthogo- nality of the cross-section to the center line after de- formation is assumed to be preserved and the cross- sectional dimensions remain unchanged. The formu- lation accounts for the geometrically nonlinear behav- ior, nevertheless only the linear behavior is studied in the presented examples. The element has four de- grees 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 ψ. 12 http://dx.doi.org/10.14311/APP.2021.30.0012 http://ojs.cvut.cz/ojs/index.php/app vol. 30/2021 Comparison of Timoshenko and Bernoulli beams N4,2 ⇠ Knotspan Patch Continuity: C0 C0C1C1 C1 Knotspan p = 2, ⌅ = {0, 0, 0, 0.5, 1, 1, 1} N1,2 N2,2 N3,2 Knots Control points P1 P2 P4P3 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 vec- tor Xc (xc) Xc = ! i R p i X̂ i, (2) xc = ! i R p i x̂ i, (3) where X̂i are the control points coordinates and Rpi are the corresponding basis functions. Deformed po- sition vector is given as x̂i = X̂i + ûi, (4) where ûi is a vector of global degrees of freedom (u, v, w). The cross-section orientation is described by a mov- ing trihedral given by the base vectors Ai (ai) with i ∈ {1, 2, 3}. With assumption of Bernoulli kinemat- ics, a position vector of an arbitrary point of a beam continuum given by coordinates X (x) can be ex- pressed as X(θ1, θ2, θ3) = Xc(θ1) + θ2A2(θ1) + θ3A3(θ1), (5) x(θ1, θ2, θ3) = xc(θ1) + θ2a2(θ1) + θ3a3(θ1), (6) where θi are the convective contravariant coordi- nates. The first components of a moving trihedral A1 and a1 are aligned with a normalized tangents T and t, re- spectively, 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. A1) and are derived from the reference trihedral, given by A0α, T0, using two key operations used for the alignment at the desired po- sition. 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 defini- tion of Green-Lagrange strain tensor, the individual components of the tensor can be derived as E11 = 1 2 (a11 − A11) " #$ % " +θ2 a2,1 · a1 − A2,1 · A1" #$ % κ21 +θ3 a3,1 · a1 − A3,1 · A1" #$ % κ31 . (7) The diagonal terms Eαα as well as the shear term E23 are equal to zero. This yields from the Bernoulli assumptions (undeformable cross-section). For the off-diagonal terms E12, E13 E1α = 1 2 θβ (aβ,1 · aα − (Aβ,1 · Aα)" #$ % κβα , (8) where (α, β) ∈ {(2, 3), (3, 2)}. The relations between strains and displacements are not linear and to the purposes of the linear analy- sis 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 stan- dard procedure [5]. The fully derived equations for the presented beam formulation can be found in [1]. 4. Timoshenko beam element formulation The second beam element presented here has been proposed by G. Zhang et al. [2]. The beam formu- 13 Edita Dvořáková, Bořek Patzák Acta Polytechnica CTU Proceedings A1(✓ 1) A2(✓ 1) A3(✓ 1) a1(✓ 1) a3(✓ 1) a2(✓ 1) u(✓1, ✓2, ✓3) X(✓1, ✓2, ✓3) Xc(✓ 1) xc(✓ 1) x y z Figure 2. Illustration of the Bernoulli beam in its undeformed and deformed configurations. T0 T0 A02 A03 T ⇤(T0, T) ⇤(T0, T)⇤(T0, T)A 0 3 A03 ⇤(T0, T)A 0 2 T0 A02 A03 T ⇤(T0, T)A 0 3 ⇤(T0, T)A 0 2 R̄T ( ) R̄T ( ) A3 A2 Figure 3. The Λ(T0, T) and R̄T(Ψ) operations used for the alignment of the Bernoulli beam cross-section at the current position. t(s) n(s) b(s) ut ✓t ub ✓b un ✓n r(s) Figure 4. Local coordination system and degrees of freedom of Timoshenko beam element. 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 dis- placement ut, normal displacement un, binormal dis- placement ub and rotations θt, θn, θb (see Fig. 4). A strain-displacement matrix B is defined as ε = Br, (9) where ε = {εm, γn, γb, χt, χn, χb}T and r = {ut, un, ub, θt, θn, θb}T . B is derived from the follow- ing relations for membrane strain εm, shear strains γn and γb, torsional strain χt, and bending strains χn and χb εm = u′t − κun, χt = θ′t − κθn, γn = κut + u′n − τ ub − θb, χn = κθt + θ′n − τ θb, (10) γb = τ un + u′b + θn, χb = τ θn + θ′b. Curvature κ is given as κ = &&&& &&&& d2r(s) ds2 &&&& &&&& (11) and torsion τ τ = dn(s) ds b(s), (12) where r(s) is the position vector and s is the curvilin- ear coordinate measured along the centerline of the beam. Conjugated generalized stress vector contains nor- mal and shear forces, and torsional and bending mo- ments. Similarly to the previous formulation, by fol- lowing standard procedure the stiffness matrix is ob- tained. 4.1. Numerical locking The well-known drawback of the Timoshenko beam element formulations is numerical locking. Locking phenomena rises from the use of the same approxima- tion order for displacements and for rotations, which leads to the shear strains of a higher order than bend- ing moments. Several methods for unlocking beam elements are available, B̄-method [6] is used for the purposes of 14 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. 5. Numerical examples Both presented beam formulations have been imple- mented 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. 5.1. 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 ex- ample of the thin helix depicted in Figure 5. The convergence of the tip vertical displacement has been observed. Additionally, the results have been com- pared 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. When B̄-locking treatment is used a superior conver- gence is obtained. The slow convergence of the tra- ditional FEM beam is due to the inability to exactly represent curved geometry. See Figure 6 for the re- sults. On the other hand, isogeometric formulation with the locking treatment provides superior conver- gency. 5.2. Convergence comparison In order to compare the Timoshenko and the Bernoulli beam formulations the error of the verti- cal tip displacement has been calculated as e = wref − w wref . (13) 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 con- vergence is provided by the presented Bernoulli beam formulation, nevertheless, the Timoshenko beam ele- ment performs also very well. 5.3. Applicability of Bernoulli formulation Finally, the applicability of the Bernoulli beam for- mulation has been assessed. The reference solution of the Timoshenko beam formulation has been used r1 r2 Helicoidal spring: R = 2 m H = 3 m Load: f = 3 kN/m Cross-section: r1 = 0.11 m r2 = 0.15 m Material: E = 200 GPa Figure 5. Helicoidal spring cantilever subjected to the constant force load and its cross-section. 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 ra- tios. The results can be observed in Table 1. As ex- pected, the error rises significantly with the increas- ing radius of the cross-section. r1 0.11 0.31 0.51 0.71 0.91 r2 0.15 0.35 0.55 0.75 0.95 r2/L 0.012 0.027 0.042 0.058 0.074 Error 0.7% 4.5% 10.9% 18.7% 27.1% Table 1. The error of Bernoulli beam formulation to the Timoshenko reference solution for different cross- sections. 6. Conclusions Two isogeometric beam elements have been imple- mented: 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] with B̄- method [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 suit- ability of the locking treatment used has been ver- ified. Secondly, the convergence properties of both 15 Edita Dvořáková, Bořek Patzák Acta Polytechnica CTU Proceedings -0.023 -0.0229 -0.0228 -0.0227 -0.0226 -0.0225 -0.0224 -0.0223 20 40 60 80 100 120 140 160 Ve rti ca l t ip d is pl ac em en t Number of control points (nodes) REF Timoshenko Timoshenko BBAR FEM -0.02293 -0.02292 -0.02291 -0.02290 -0.02289 -0.02288 -0.02287 20 40 60 80 100 120 140 160 Ve rti ca l t ip d is pl ac em en t Number of control points (nodes) REF Timoshenko BBAR FEM Figure 6. Comparison of Timoshenko beam element with no locking treatment with Timoshenko beam element using B̄-method and with naturally locking-free standard FEM beam element. The vertical axis corresponds to the vertical tip displacement w, the horizontal axis corresponds to the number of nodes (control points). 16 vol. 30/2021 Comparison of Timoshenko and Bernoulli beams 0 1x10-6 2x10-6 3x10-6 4x10-6 5x10-6 6x10-6 7x10-6 8x10-6 9x10-6 30 35 40 45 50 55 60 65 70 75 E rr or o f t he v er tic al ti p d is pl ac em en t Number of control points (nodes) Bernoulli Timoshenko BBAR Figure 7. Comparison of covergence of Timoshenko beam element and Bernoulli beam element. The vertical axis corresponds to the erro of vertical tip displacement w, the horizontal axis corresponds to the number of nodes (control points). The error is calculated with respect to the reference solution of particular beam theory. 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 su- perior convergence compared to the traditional FEM formulation can be observed for both formulations. Acknowledgements The financial support of this research by the Grant Agency of the Czech Technical University in Prague (SGS project No. SGS20/038/OHK1/1T/11) is gratefully ac- knowledged. References [1] A. Bauer, M. Breitenberger, B. Philipp, et al. Nonlinear isogeometric spatial bernoulli beam. Computer Methods in Applied Mechanics and Engineering 303:101–127, 2016. [2] G. Zhang, R. Alberdi, K. Khandelwal. Analysis of three-dimensional curved beams using isogeometric approach. Engineering Structures 117:560–574, 2016. [3] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Engrg 194:4135–4195, 2005. [4] L. Piegl, W. Tiller. The NURBS Book. Springer-Verlag Berlin Heidelberg, New York, 1997. [5] Z. Bittnar, J. Šejnoha. Numerické metody mechaniky 1. České vysoké učení technické, Praha, 1992. [6] R. Bouclier, T. Elguedj. Locking free isogeometric formulations of curved thick beams. Comput Methods Appl Mech Engrg 245-246:144–162, 2012. 17