Acta Polytechnica CTU Proceedings doi:10.14311/APP.2020.26.0024 Acta Polytechnica CTU Proceedings 26:24–29, 2020 © Czech Technical University in Prague, 2020 available online at https://ojs.cvut.cz/ojs/index.php/app ON EVALUATION OF THE THREE-DIMENSIONAL ISOGEOMETRIC BEAM ELEMENT 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. The exact description of the arbitrarily curved geometries, including conic sections, is an undeniable advantage of isogeometric analysis (IGA) over standard finite element method (FEM). With B-spline/NURBS approximation functions used for both geometry and unknown approximations, IGA is able to exactly describe beams of various shapes and thus eliminate the geometry approximation errors. Moreover, naturally higher continuity than standard C0 can be provided along the entire computational domain. This paper evaluates the performance of the nonlinear spatial Bernoulli beam adapted from formulation of Bauer et al. [1]. The element formulation is presented and the comparison with standard FEM straight beam element and fully three-dimensional analysis is provided. Although the element is capable of geometrically nonlinear analysis, only geometrically linear cases are evaluated for the purposes of this study. Keywords: Bernoulli theory, curved beam element, isogeometric analysis, NURBS. 1. Introduction The desire for the automatic connection between Computer-aided design (CAD) and Finite element analysis (FEA) has been the crucial impulse for the development of isogeometric analysis [2]. The idea of IGA is to use the basis functions used for the geom- etry description in CAD also as the approximation functions for the analysis. This results in the possibil- ity of only one model shared between the design and analysis. Isogeometric approach introduces into the analysis some other very important features. One of them is a possibility to exactly model arbitrarily curved geometries, including a conic sections which can be only approximated by standard polynomial basis func- tions. This makes it very convenient for the use in the analysis of curved beams. The focus of this paper is placed on the geometri- cally nonlinear three-dimensional Bernoulli beam [1, 3]. The beam formulation is briefly introduced and the performance over standard FEM approaches is evalu- ated. 2. NURBS-based analysis The most wide-spread technology in CAD indus- try are NURBS (Non-Uniform Rational B-Splines). 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) Each NURBS function is generated by weighting B- spline functions Npi with a given weight wi associated with control point Pi R p i (ξ) = Ni,p(ξ)wi∑n j=1 Nj,p(ξ)wj . (2) B-spline functions form a special subset of NURBS functions (corresponding to constant wi) and are de- rived recursively starting with a piecewise constant functions Ni,0(ξ) = { 1 if ξi ≤ ξ < ξi+1 0 otherwise , (3) Ni,p(ξ) = ξ − ξi ξi+p − ξi Ni,p−1(ξ) + ξi+p+1 − ξ ξi+p+1 − ξi+1 Ni+1,p−1(ξ), (4) where p is the degree of the B-spline function, ξi is the coordinate of the ith-knot and ξ ∈ 〈0, 1〉 is a parametric coordinate. 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). Understanding of knotspans is similar to elements in standard FEM, nevertheless higher conti- nuity up to Cp−1 between knotspans can be achieved naturally using NURBS, unlike C0 in standard FEM. Moreover, the NURBS basis functions are generally non-interpolatory. In the knots (points which are di- viding patch into knotspans), the continuity can be locally decreased up to C0 by knot multiplication, 24 https://doi.org/10.14311/APP.2020.26.0024 https://ojs.cvut.cz/ojs/index.php/app vol. 26/2020 On evaluation of the three-dimensional isogeometric beam element N4,2 ξ Knotspan Patch Continuity: C0 C0C1C∞ C∞ 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. which is a standard technique in NURBS technology. This procedure imposes the interpolatory behavior of the functions at the particular point. For better understanding of NURBS see Fig. 1, where mapping between real and parametric geometry is illustrated. For more information the reader could refer to [4]. 3. Beam element formulation The formulation of the presented three-dimensional beam element [1] is based on Bernoulli beam the- ory and accounts for the geometrically nonlinear be- haviour. Bernoulli kinematics assume that the orthog- onality of the cross-section to the center line is pre- served after deformation and that the cross-sectional dimensions remain unchanged. The element has four degrees of freedom in each control point: three global displacements u,v, and w and additional degree of freedom corresponding to rotation around the cen- ter line ψ. The rotational degree of freedom enables the element to develop a torsion (warping effects are neglected here) and also to model initially twisted beams. 3.1. Geometric description In the following, the standard notation using upper- case and lower-case letters for the undeformed and deformed configuration, respectively, is adapted. The beam formulation (see Figure 2) starts from three- dimensional approximation reduced using Bernoulli kinematic assumptions to the displacements of the cen- ter line, which is given by the position vector Xc (xc). The position vector of the center line is described as a linear combination of the control points coordinates X̂i and the corresponding basis functions Rpi Xc = ∑ i R p i X̂ i, (5) xc = ∑ i R p i x̂ i. (6) Deformed position vector is given as x̂i = X̂i + ûi, (7) where ûi is a vector of global degrees of freedom (u,v,w). Additionally to the center line, the cross-section orientation is described by a moving trihedral given by the base vectors Ai (ai) with i ∈{1, 2, 3}. A position vector of an arbitrary point of a beam continuum given by coordinates X (x) can be expressed as X(θ1,θ2,θ3) = Xc(θ1) + θ2A2(θ1) + θ3A3(θ1), (8) x(θ1,θ2,θ3) = xc(θ1) + θ2a2(θ1) + θ3a3(θ1), (9) where θi are the convective contravariant coordinates. The first components of a moving trihedral A1 and a1 are aligned with a normalized tangents T and t, respectively, and are computed as A1 = ∑ i R p i,1X̂ i, (10) a1 = ∑ i R p i,1x̂ i, (11) where (·),i denotes the derivative with respect to the coordinate θi. The remaining components Aα with α ∈{2, 3} are orthogonal to the tangent of the center line (resp. A1) and are described by the reference trihedral, given by A0α, T0, as Aα = R̄T(Ψ)Λ(T0, T)A0α, (12) where Ψ is a initial rotation about a tangent of a beam, and Λ and R̄T are the two key operations used 25 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 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 cross-section at the current position. for the alignment of a moving trihedral illustrated in Figure 3. While the mapping matrix Λ(T0, T) aligns the reference trihedral given by the tangent T0 to the tangent T at the current position, the rotation matrix R̄T(Ψ) rotates the aligned reference trihedral about T with given rotation Ψ. The rotation Ψ = Ψ(θ1) is calculated using basis functions and values Ψ̂i assigned to the control points. The Euler-Rodriguez formula [5] is used to define both the mapping matrix Λ and the rotation matrix RT. Analogical procedures denoted as Λ(T, t) and R̄t(ψ) are used to align the moving trihedral of the undeformed to the deformed configuration. The base vectors aα of the deformed configuration are given as aα = R̄t(ψ)Λ(T, t)Aα, (13) where ψ = ψ(θi) is a rotational degree of freedom evaluated at the current position θi. The base vectors of the undeformed continuum Gi are defined as Gi = X,i leading to G1(θ1,θ2,θ3) = A1(θ1) + θ2A2,1(θ1) + θ3A3,1(θ1), G2(θ1) = A2(θ1), (14) G3(θ1) = A3(θ1) with analogical equations for the base vectors gi = x,i of the undeformed configuration. In the sequel, dot products Gi · Gj ( gi · gj) and Ai · Aj ( ai · aj) are denoted as Gij (gij) and Aij (aij), respectively. 3.2. Green-Lagrange strain tensor The Green-Lagrange strain tensor E calculated for the curvilinear coordinate system is defined as Eij = 1 2 (gij −Gij) . (15) and for the orthogonal base vectors, the transformation to the Cartesian coordinate system denoted with (̃·) is given by Ẽij = Eij ‖Gi‖2‖Gj‖2 . (16) By substituting geometric relations into the strain ten- sor, the individual components of the tensor can be derived. During the derivation additional simplifica- tions are made assuming that only a slender beam is considered (b,h << L where b and h are cross-sectional dimensions and L is a length of the beam) resulting in the following equations for the diagonal term E11 E11 = 1 2 (a11 −A11)︸ ︷︷ ︸ � +θ2 a2,1 · a1 − A2,1 · A1︸ ︷︷ ︸ κ21 +θ3 a3,1 · a1 − A3,1 · A1︸ ︷︷ ︸ κ31 . (17) 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α)︸ ︷︷ ︸ κβα , (18) where (α,β) ∈{(2, 3), (3, 2)}. 26 vol. 26/2020 On evaluation of the three-dimensional isogeometric beam element 3.3. Constitutive equations The energetically conjugated stress tensor to the Green-Lagrange strain tensor is the Second Piola- Kirchhoff tensor S. Elastic isotropic material is con- sidered within this work and St. Venant–Kirchhoff material is applied. As the beam formulation is re- duced to the center line and Bernoulli theory is as- sumed, the shear forces S̃23 and S̃32 and normal forces S̃22 S̃33 vanish and the full constitutive relation can be reduced to  S̃11S̃12 S̃13   =   E 0 00 G 0 0 0 G   ︸ ︷︷ ︸ D   Ẽ11Ẽ12 Ẽ13   , (19) where E is Young’s modulus, G is the shear modulus and D is a reduced elasticity matrix. 3.4. Principle of Virtual Work The system is in equilibrium when the overall virtual work of the system is equal to zero δW = −δWint + δWext = 0. (20) The internal end external virtual work of the system is given by δWint = ∫ Ω S : δEdx, (21) δWext = ∫ ΓN t : δudx + ∫ Ω ρ0B : δudx, (22) where Ω and ΓN are the domain and the Neumann boundary, respectively. The external virtual work depends on the boundary forces t, body forces B and material density ρ0. By substitution of expressions for the virtual Green- Lagrange strain tensor obtained from (15), constitu- tive relations (19) and performing linearization the following incremental relation can be obtained ∑ s Krs∆us = −Rr, (23) where ∆us are the displacement increments and Rr = −(Fintr + F ext r ) (24) = ∂Wint ∂ur + ∂Wext ∂ur , (25) Krs = ∂Rr ∂us = ∂2Wint ∂ur∂us . (26) 4. Numerical example The described beam element has been implemented and tested on the example of helicoidal spring (see Figure 4). This problem represents fully three- dimensional structure where the all membrane, bend- ing and torsion effects are present. The geometry is F = 1 t t E = 108 ν = 0.0 Figure 4. Helicoidal spring cantilever subjected to the unit force tip load and its cross-section. Figure 5. Helicoidal spring geometry for different thickness/length ratios: 0.1, 0.05, 0.025, 0.0125. given by Xc(θ1) =   10 sin ( θ12π ) 10 cos ( θ12π ) 20 θ1   . (27) The structure is subjected to the unit tip force load. Different dimensions of the square cross-section have been used in the analysis, in order to evaluate the performance and the applicability of the underlying assumptions. The studied geometries can be observed in Figure 5, where helicoidal springs with different thickness/length (t/L) ratios (0.1, 0.05, 0.025 and 0.0125) are illustrated. The results have been compared with results ob- tained using standard FEM straight beam element with cubic approximation [6]. This element allows for the analysis of both thick and thin beams and is naturally locking-free. The beam has been used in two different configurations: i) the element satisfying Bernoulli assumptions; ii) the element accounting also for the shear effects. The results can be observed in the Figure 6. As expected, for the thick beam configuration the results obtained using classical FEM or IGA elements with Bernoulli assumptions differ from the standard FEM beam element with account for the shear deforma- 27 Edita Dvořáková, Bořek Patzák Acta Polytechnica CTU Proceedings -1.01e-06 -1.00e-06 -9.90e-07 -9.80e-07 -9.70e-07 -9.60e-07 -9.50e-07 -9.40e-07 -9.30e-07 -9.20e-07 -9.10e-07 10 15 20 25 30 35 40 45 t/l = 0.1 -1.60e-05 -1.58e-05 -1.56e-05 -1.54e-05 -1.52e-05 -1.50e-05 -1.48e-05 -1.46e-05 10 15 20 25 30 35 40 45 t/l = 0.05 -2.55e-04 -2.50e-04 -2.45e-04 -2.40e-04 -2.35e-04 -2.30e-04 10 15 20 25 30 35 40 45 t/l = 0.025 -4.05e-03 -4.00e-03 -3.95e-03 -3.90e-03 -3.85e-03 -3.80e-03 -3.75e-03 -3.70e-03 10 15 20 25 30 35 40 45 t/l = 0.0125 IGA p=3 IGA p=4 IGA p=6 FEM - thin FEM - shear -1.01e-06 -1.00e-06 -9.90e-07 -9.80e-07 -9.70e-07 -9.60e-07 -9.50e-07 -9.40e-07 -9.30e-07 -9.20e-07 -9.10e-07 10 15 20 25 30 35 40 45 t/l = 0.1 -1.60e-05 -1.58e-05 -1.56e-05 -1.54e-05 -1.52e-05 -1.50e-05 -1.48e-05 -1.46e-05 10 15 20 25 30 35 40 45 t/l = 0.05 -2.55e-04 -2.50e-04 -2.45e-04 -2.40e-04 -2.35e-04 -2.30e-04 10 15 20 25 30 35 40 45 t/l = 0.025 -4.05e-03 -4.00e-03 -3.95e-03 -3.90e-03 -3.85e-03 -3.80e-03 -3.75e-03 -3.70e-03 10 15 20 25 30 35 40 45 t/l = 0.0125 IGA p=3 IGA p=4 IGA p=6 FEM - thin FEM - shear -1.01e-06 -1.00e-06 -9.90e-07 -9.80e-07 -9.70e-07 -9.60e-07 -9.50e-07 -9.40e-07 -9.30e-07 -9.20e-07 -9.10e-07 10 15 20 25 30 35 40 45 t/l = 0.1 -1.60e-05 -1.58e-05 -1.56e-05 -1.54e-05 -1.52e-05 -1.50e-05 -1.48e-05 -1.46e-05 10 15 20 25 30 35 40 45 t/l = 0.05 -2.55e-04 -2.50e-04 -2.45e-04 -2.40e-04 -2.35e-04 -2.30e-04 10 15 20 25 30 35 40 45 t/l = 0.025 -4.05e-03 -4.00e-03 -3.95e-03 -3.90e-03 -3.85e-03 -3.80e-03 -3.75e-03 -3.70e-03 10 15 20 25 30 35 40 45 t/l = 0.0125 IGA p=3 IGA p=4 IGA p=6 FEM - thin FEM - shear Figure 6. Comparison of isogeometric beam element with cubic (IGA p = 3), quartic (IGA p = 4) and hexic (IGA p = 6) approximation with standard FEM straight beam element with Bernoulli assumptions (FEM - thin) and with account for shear effects (FEM - shear). The vertical axis corresponds to the vertical tip displacement w, the horizontal axis corresponds to the number of nodes (control points). tion. This difference is diminishing with the smaller thickness/length ratio. In addition to the beam elements, the fully three- dimensional analysis has been assessed. The enormous number of nodes, in comparison with beam elements, have to be used in order to obtain sufficiently accu- rate results. Finally, the mesh consisting of 102800 nodes (20×20×257) has been used. For the thin beam with t/L = 0.0125, the obtained deflection is w = −3.877e−3. This result is slightly higher than the results obtained using beam assumptions. The difference between the full 3D and beam models is getting bigger as the height to length ratio increases. This difference is indicating, that the beam theory as- sumptions are no longer adequate. Overall, the results demonstrate the superior convergence properties of the IGA element. This is due to its ability to capture curved geometry and natural high-order interpolation. 5. Conclusions The spatial isogeometric element based on the work of A. M. Bauer et al. [1] has been implemented and its convergence properties have been evaluated in this contribution. The element formulation is geometri- cally nonlinear, however only the geometrically linear analysis has been considered in this paper. The per- formance of the presented element has been compared with standard FEM beam element using benchmark problem of helicoidal spring. For the higher thickness/length ratios t/L = 0.1, 0.5 the element does not provide good results due to the limitations of Bernoulli assumptions. For the thin beams analysed, the excellent results are obtained. Acknowledgements The financial support of this research by the Grant Agency of the Czech Technical University in Prague (SGS project No. SGS19/032/OHK1/1T/11) is gratefully acknowledged. 28 vol. 26/2020 On evaluation of the three-dimensional isogeometric beam element References [1] A. M. Bauer, M. Breitenberger, B. Philipp, et al. Nonlinear isogeometric spatial bernoulli beam. Computer Methods in Applied Mechanics and Engineering 303:101–127, 2016. [2] 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. [3] L. Greco, M. Cuomo. B-spline interpolation of Kirchhoff-Love space rods. Computer Methods in Applied Mechanics and Engineering 256:251 – 269, 2013. doi:10.1016/j.cma.2012.11.017. [4] L. Piegl, W. Tiller. The NURBS Book. Springer-Verlag Berlin Heidelberg, New York, 1997. [5] J. S. Dai. Euler–rodrigues formula variations, quaternion conjugation and intrinsic connections. Mechanism and Machine Theory 92:144–152, 2015. [6] Z. Bittnar, J. Šejnoha. Numerické metody mechaniky 1. České vysoké učení technické, Praha, 1992. 29 https://doi.org/10.1016/j.cma.2012.11.017 Acta Polytechnica CTU Proceedings 26:24–29, 2020 1 Introduction 2 NURBS-based analysis 3 Beam element formulation 3.1 Geometric description 3.2 Green-Lagrange strain tensor 3.3 Constitutive equations 3.4 Principle of Virtual Work 4 Numerical example 5 Conclusions Acknowledgements References