Linear static isogeometric analysis of an arbitrarily curved spatial Bernoulli-Euler beam

: The present research is focused on the linear analysis of a spatial Bernoulli-Euler beam. Metrics of the reference and deformed configurations are rigorously defined with respect to the convective coordinate frame of reference. No higher order terms are neglected which makes the formulation ideally suited for analysis of arbitrarily curved spatial beams in the frame of finite (but small) strain theory. The well-known issue of nonorthogonality of local coordinate system at an arbitrary point of a spatial beam is solved by the introduction of a new coordinate line that is orthogonal to the normal plane of the beam axis at each point. Generalized coordinates of the present model are translations of the beam axis and the angle of twist of a cross section. Two different parameterizations of this angle are discussed and implemented. Both geometry and kinematics are described with the same set of NURBS functions, in line with the isogeometric approach. Numerical analysis proved that the theoretical considerations are correct. An in-depth analysis of convergence properties has confirmed the fact that models with the highest interelement continuity have an improved accuracy per DOF, for problems that result in a smooth structural response.


Introduction
Efficient application of curved spatial beams represents a challenge for contemporary engineering.Undeniable aesthetical characteristics of these structures make them appealing for designers while, on the other hand, a correct understanding of their mechanical behavior is not a trivial task.Additionally, development of new materials enables construction of complex structural shapes, not feasible just a few decades ago.These trends are followed by improvements of numerical models for analysis of these structures.Some of the most notable contributions since the 1980s are from Reissner [1], Simo [2] and [3], Iura and Atluri [4] and [5], Ibrahimbegović [6], Crisfield [7], Kapania [8] and Pi [9].Compact reviews of the development of computational models for analysis of curved beams, with historical perspective, can be found in [10] and [11], where a special attention is paid to the nonlinear analysis.
One of the most prominent influences on the development of computational mechanics in this century came from the introduction of isogeometric analysis (IGA) by Hughes and co-workers, [12] and [13].The motivation for the development of IGA stems from the problem of integration of design and analysis.Traditionally, different approaches are applied for these two inseparable segments of structural modeling which result in the inefficiency of the whole design process.The solution is accomplished by a groundbreaking approach where the reference geometry and kinematics are described with the same set of basis functions, mostly some type of spline functions.Although the same logic is valid for the standard isoparametric finite element (FE) analysis [14], the main difference is the fact that the isogeometric approach requires a primary description of the reference geometry.This results in many advantages over the classic reverse approach of the isoparametric FE analysis, the exact representation of the reference geometry being the most obvious one.Since its origination, application of IGA found its place in almost all areas of mechanics, [15].
Considerations of spatial Bernoulli-Euler beams (BE) in the frame of IGA are not frequent and works of Raknes et al. [16], Bauer et al. [17], as well as series of papers by Greco and Cuomo [18], [19] and [20], stand out.Concretely, a spatial cable formulation without torsion is developed purely in terms of displacement DOFs in [16] with attractive nonlinear dynamic applications.However, cables with thin circular cross sections are considered solely.The bending strip method is introduced for modeling of multiple NURBS patches, which is not straightforward in the rotation-free analysis.Authors in [18] introduced a spatial Kirchhoff-Love rod modeled as a ribbon using B-splines and Bezier interpolations.The model has four DOFs, three translations and one rotation around the tangent, which is only independent rotation component of the BE beam.The same authors extended their approach to the assemblies of Kirchhoff-Love space rod elements, introducing a multi-patch implicit G 1 formulation, [19].Furthermore, the authors suggested an usage of mixed approach in [20] for the analysis of spatial non-polar elements.Based on [18], the authors in [17] successfully introduced a geometrically nonlinear spatial curved BE beam including torsion.Furthermore, a different methodology for the so-called cylindrical element was developed in [21], which enabled the analysis of a rod described as a full 3D continuum, while only surface control points are parameterized.
Regarding the shear deformable spatial beams, a collocation technique is successfully applied for linear IGA in [22] where the simple constitutive matrix is used and coupling between axial and bending effects are not considered, which is valid only for beams with a small value of curviness.A similar approach is presented in [23].
The main drawback of almost all mentioned references, regarding IGA of arbitrarily curved spatial beams, is the fact that they do not consider complete beam metric, as explicitly stated in [17].Namely, beams are considered very slender and square order terms with respect to the principal axes of inertia are neglected.Additionally, the tangential base vector at an arbitrary point of a cross section is approximated with the tangential base vector at the centroid, for the purposes of calculating values of metric coefficients, [17].The latter assumption is necessary in order to decouple torsion from axial straining as well as to exclude warping.However, this is here performed strictly by the introduction of a new coordinate line.Although beam theories are mostly valid for slender beams, it is of substantial theoretical interest to study geometrically exact beam models, regardless of the level of slenderness.
The standard measure of slenderness for straight beams is not straightforwardly applicable for curved beams since the effect of initial curvature must be also considered.Therefore, for planar circular beams, product Kh is often introduced where K is the curvature of the centroid axis in the reference configuration, while h is the height of a cross section of a uniform beam, [24].Moreover, for beams with a variable curvature, the maximum value of curvature of the beam axis, Kmax, should be used.It follows that Kmaxh is one valid measure of curviness for beams with a smoothly distributed curvature.For spatial beams, it is more complicated to define this parameter since the height of a cross section h must be replaced with the maximum dimension of the cross section in the plane parallel to the osculating plane, hmax.Evidently, this measure of slenderness lacks more explicit introduction of the influence of beam length and it requires amelioration that will be considered in further research.It should be pointed out that the term arbitrarily curved does not refer just to the variable curvature of the beam axis, but, more importantly, to its curviness defined as Kmaxhmax.More discussion and results on this topic, regarding the plane case, can be found in [25], [26] and [27].
In the present research, the isogeometric finite element is formulated utilizing the exact metric of an arbitrarily curved spatial BE beam.Only assumptions are those from the BE hypothesis itself, as well as those regarding finite (but small) strain theory.The complete derivation is presented with sufficient details in order to enable reproducibility.The model has four DOFs: three translations of the beam axis and one rotation of a cross section around the tangential basis of the beam axis.Two formulations are postulated, depending on the rotational term that is used as DOF.All quantities, except this rotational term, are expressed with respect to the global Cartesian frame of reference.One of the prominent contributions of the present research is the introduction of the new orthogonal coordinate system at each point of the beam, which enabled an elegant composition of equations.The present study is based on the theory given in [28] and [29] with considerable meliorations and novel examples.
The paper is organized into six sections.The next section deals with some basic concepts of the NURBS-based IGA modeling.It is followed by a detailed derivation of the reference and deformed metric of the spatial BE beam.The fourth section gives a derivation of virtual work via energetically conjugated pairs of section forces and strains, which resulted with appropriate stiffness matrix and vector of external load.Afterward, two numerical examples are considered in detail with a thorough convergence check and comparison with relevant results.At the end section, brief conclusions are drawn.

NURBS-based isogeometric analysis
Following the isogeometric concept, both geometry and kinematics are discretized with the same Non-Uniform Rational B-spline (NURBS) functions.Literature devoted to the NURBS-based IGA is growing exponentially and a detailed discussion on this topic can be found in [30] and [13].The aim of this section is to deliver a brief introduction of the NURBS-based isogeometric modeling as a background to the present work.

Non-uniform rational B-spline
B-spline, or basis spline, is a piecewise polynomial function expressed as a linear combination of B-spline basis functions BI,p(ξ) and coefficients PI which represent coordinates of control points: The piecewise linear function constructed with control points is called control polygon.In general, control points are not interpolatory to a B-spline curve.Furthermore, the B-spline function is defined by a set of non-decreasing real numbers called knot vector ξ={ξ1, ξ2,…, ξN+p+1}, where N is a number of basis functions while p is a polynomial order.Elements ξi are called knots and they represent coordinates of control points in the parametric space, Fig. 1.Geometries defined by a single patch, i.e. interval [ξ1, ξN+p+1], are considered here solely.Additionally, a knot span is defined by an interval (ξi, ξi+1) and it exists in both parent and physical domain, Fig. 1.In the present research, each knot span, for which ξi≠ξi+1, corresponds to an isogeometric finite element.Multiplication of knots in a knot vector generates patches with different interelement continuities.If the same knot is repeated m times, B-spline functions are C p-m continuous at that knot.It is interesting to note that a knot becomes interpolatory for a case of C 0 interelement continuity.However, the present considerations of the spatial BE beam require utilization of, at least, C 1 continuity.Moreover, if the first and the last knots are repeated p+1 times, the so-called open knot vector is obtained.Since this type of knot vector is utilized here exclusively, it should be stressed out that basis functions are interpolatory at the ends of observed function.An important property that follows from this fact is that the control polygon is tangential to a curve at these points.Additionally, unequal knot spans constitute a non-uniform knot vector, which is more versatile than its uniform counterpart.
The B-spline basis functions are defined by the Cox-de Boor recursion formula: .
Three important properties of B-splines are emphasized from the perspective of their utilization in the numerical analysis: partition of unity, pointwise non-negativeness, and local support property, [12].Naturally, B-spline curves inherit their properties directly from their basis functions.For a more detailed discussion on these matters, see [13].
Generalization of the B-splines generates a NURBS, a significantly more versatile type of basis functions.This is due to fact that the NURBS curve in ℝ d is actually a projection of the corresponding B-spline curve in ℝ d+1 , [30].If we introduce rational basis functions RI,p and weights wI, associated with the I th control point, a NURBS curve can be defined as: , . ( It follows that rational basis functions share the same main properties with the B-spline basis functions, which are then inherited by NURBS curves.For clarity, two spatial cubic NURBS curves, which slightly differ in the weight associated with the control point P5, are displayed in Fig. 2a.This weight affects the curve only at the interval (1,3).The previously mentioned fact that the control polygon is tangent to a curve at its ends is observed.Additionally, the curve is C 0 continuous at ξ=3 and control point P6 is interpolatory since this knot value is repeated p=3 times in the knot vector.Matching basis functions are given in Fig. 2b where the effect of different weights at the interval (1,3), as well as C 0 interelement continuity at ξ=3 are evident.Also, all three emphasized important properties of NURBS, inherited by B-splines, are evident from Fig. 2b.

Mesh refinement strategies
Besides the classic hp-refinement strategies, well-known from the isoparametric FE approach, [14], a novel type of mesh refinement is introduced in IGA, the k-refinement.Detailed considerations of the hpk-refinement procedures can be found in [31] while only basic ideas are presented here.
Principal mesh refinement strategy used in FE analysis is the h-refinement.It refers to a decrease of characteristic length h of elements utilized for the mesh.For IGA, it implies an insertion of new knots in a knot vector.The main constraint, which must be satisfied during this knot insertion, is that a curve stays the same.This constraint leads to a definition of new basis functions and new control points, which are formed as linear combinations of the original ones.Consequently, the h-refinement results with a splitting of existing elements into new ones.Repeating of existing knots does not create new elements but decreases interelement continuity.In this research, the h-refinement implies the insertion of new elements while continuity and polynomial order remain fixed.
The second standard refinement strategy is based on the elevation of polynomial order used for approximation while number and size of elements stay the same.It is referred to as the p-refinement and, similar to the h-refinement, shape, and dimension of the original curve must not change.An increase of polynomial order p is followed with the multiplication of knots, which preserves interelement continuity at each knot.
Finally, the IGA specific procedure is the k-refinement, which is enabled by fact that the knot insertion and order elevation do not commute.This property allows us to apply the p-refinement on the coarsest mesh, followed with the h-refinement.In this way, interelement continuity can be exactly prescribed.The k-refinement introduced many interesting features into the NURBS-based isogeometric structural modeling and concrete conclusions are yet to be drawn.

Metric of the spatial Bernoulli-Euler beam
A detailed and rigorous definition of beam metric is presented in this section.The classic BE assumption that a beam's cross section is rigid and perpendicular to the beam axis in the deformed configuration is applied.This presumption leads to the degeneration of a 3D continuum model of the beam into an arbitrarily shaped line, Fig. 3.The analysis is performed with respect to the convective frame of reference (ξ, η, ζ) utilizing its respective base vectors (g1, g2, g3).At the centroid, the η and ζ coordinate lines correspond to the principal axes of the second moment of area of a cross section, while the ξ coordinate axis coincides with the centroidal axis of the beam.Due to the BE hypothesis, the complete kinematics of beam is defined from translations of the beam axis and rotation of a cross section around the g1 basis.Since the metric of a line is naturally described with the Frenet-Serret coordinate system (s, τ, υ) and its base vectors (t, n, b), we will simultaneously follow both coordinate systems while noticing the initial angle θ between (n, b) and (g2, g3) which can be arbitrary, Fig. 3.

Metric of the beam axis
The position vector of the beam axis is r=r(ξ) while its global rectangular Cartesian components are x m =x m (ξ), where m=1,2,3 (x 1 =x, x 2 =y, x 3 =z).The derivative of this vector with respect to the ξ coordinate returns the base vector g1: where ( .),k designates differentiation with respect to the k th coordinate of the system (ξ, τ, υ).im are the unit base vectors of the global rectangular Cartesian coordinate system while t is the tangent of the centroid curve.According to the NURBS-based isogeometric approach, m I x are the Cartesian components of the position vector of the I th control point.Notice the relation between the convective and arc-length coordinate: Metric and reciprocal metric tensors of the beam axis, as well as their respective determinants, are: Furthermore, the well-known relation for the vector of curvature states: where K represents the modulus of curvature vector with respect to the Frenet-Serret frame of reference.
The relation for the normal as a function of ξ coordinate stems from this expression: ,1 where 1   11   Γ is the Christoffel symbol of the second kind defined as: If we introduce the following designation: , the normal of the beam axis can be compactly written as: ,2 .
Positions of the indices are irrelevant for the Cartesian coordinate system since there is no distinction between the basis and its dual counterpart.Furthermore, the binormal is: g n (12) and its components with respect to the Cartesian frame of reference are: , where emnk is the permutation symbol.Now, the binormal can be represented as: , 3 .
Change of the tangential base vector g1 of the beam axis with respect to the ξ coordinate is: while changes of the normal and binormal are: , , n n (16) where τ is the modulus of torsion: The relations ( 8) and ( 16) actually represent the Frenet-Serret formulae that, with respect to the arclength coordinate, state: This relation can be written as: . , , where ti are tangent, normal and binormal for i=1,2,3, respectively.Having in mind the vector property of a skew-symmetric system of the second order (spinor), we can define a pseudovector of the curvature with components defined with respect to the Frenet-Serret triad as: where K is known as the Darboux vector [32].Now, the equation ( 19) can be represented as: It is mentioned that the coordinates η and ζ coincide with the principal axes of inertia of a cross section at the centroid.They form an arbitrary angle θ with the normal and binormal of the beam axis, Fig. 3, which results in simple relations: cos sin cos sin .sin cos sin cos In order to distinguish between the Frenet-Serret triad and the convective frame of reference, the derivatives and displacement components with respect to the η and ζ axes are designated with an overbar.Accordingly, components of the curvature vector with respect to the convective frame of reference are: and they satisfy the relation ( .
Notice that the moduli of these base vectors are: Consequently, derivatives of the base vectors gi with respect to the ξ coordinate, as functions of the same set of base vectors are: where it can be observed that this tensor is not skew-symmetric, as in (18), due to fact that g1 basis is not a unit vector.
Finally, the metric of the beam axis is completed with the Christoffel symbols of the second kind: , , , Observe that the components of the curvature vector can be obtained from the relation ( 26):

Metric of an arbitrary point
Since the position vector of an arbitrary point is: the base vector with respect to the ξ coordinate is: Using the equation ( 26), the base vectors at an arbitrary point are: where the standard designation [4], is introduced: Since the base vectors, g2 and g3 are functions of one variable ξ, the metric tensor at an arbitrary point of a cross section is: while its reciprocal counterpart is: A brief analysis of the derived expressions reveals the well-known fact that base vector 1 g is not orthogonal with respect to the other two, g2 and g3, which implies that the coordinate system (ξ, η, ζ) is not orthogonal, except at the centroid.Usage of these coordinates makes equations of beam element complex and theory itself inefficient.Therefore, we will introduce a new coordinate, ξλ, as a function of ξ, η and ζ: ( ) , , , which, by itself, can represent differential of the position vector of an arbitrary point: Furthermore, it is necessary to constrain its base vector to be collinear with the tangential basis at the centroid: where the scaling factor of unit value is chosen for convenience.Notice that ξ=ξ(ξλ), η=η(ξλ) and ζ=ζ(ξλ).In order for this constraint to be satisfied, the following relations must hold: which is easily checked if (38) is inserted into (37).In this way, the new coordinate system (ξλ, η, ζ ) which is orthogonal at each point of the beam is constructed.Additionally, it is evident from (38) that ξλ and ξ coincide at the centroid.Note that (ξλ, η, ζ ) coordinates are dependent since (35) holds.
However, they represent a valid coordinate system since they uniquely represent every point of the beam, due to the fact that (ξ, η, ζ ) are independent.
It should be remarked that the convective property can also be attributed to the newly introduced ξλ coordinate.However, if the convective ξλ coordinate satisfies the conditions of orthogonality (38) in the reference configuration, it cannot satisfy them in the deformed configuration.Analogously, if the conditions of orthogonality (38) are satisfied for both reference and deformed configuration, the ξλ coordinate is not convective.These considerations are not as crucial for linear analysis as they are for nonlinear one.
Remark.It is interesting to note that the constraints (38) can be transformed into the following system of differential equations: Finally, the Christoffel symbols of the second kind, with respect to the coordinate system (ξλ, η, ζ ), follow from (38): 3.2 Deformed configuration

Metric of the beam axis
If we designate the displacement vector of the beam axis with u(ξ), the position vector of the beam axis in the deformed configuration is: An asterisk in superscript is used throughout this paper as a designation for quantities in the deformed configuration.Next, the tangential base vector of the deformed beam axis is: while the determinant of the metric tensor of the deformed beam axis is: The base vectors g2 and g3 in the deformed configuration keep the unit intensities, and they are defined as: where R is an infinitesimal rotation matrix while ϕ is the axial vector of spinor Φ: Components of this axial vector are: where the hat symbol designates local components of displacement, contrary to the ones with the respect to the Cartesian coordinates.Additionally, designation ( ) k ⋅ stands for covariant derivative with respect to the k th coordinate of the observed coordinate system.It should be emphasized that for the BE theory, the rotations ϕ 2 and ϕ 3 depend on the displacement field of the beam axis: while the rotation ϕ 1 is independent of it.Therefore, this rotation must be introduced as a generalized coordinate of the spatial BE beam: Here, ω is the physical component of total rotation of a cross section around tangential basis.According to the isogeometric concept, the adopted generalized coordinates are described with the same set of NURBS functions: where m I u is the m th displacement component at the I th control point, while ωI is its associated rotation.
Finally, from the equations ( 44) and (45) follow displacement gradients with respect to the η and ζ coordinates:

Reference strains of the beam axis
Normal strain of the beam axis along the tangential direction is: which, for linear analysis, reduces to: and it represents the first reference strain of the beam axis.Moreover, in order to define strain at an arbitrary point of a cross section, it is necessary to introduce the differences of curvatures in the deformed and reference configurations, κ1, κ2 and κ3, with respect to the convective coordinate frame.
We will designate the first of these reference strains as the torsional strain and the other two as the flexural strains of the beam axis.From the equation (28) follows: , g g g g g u g u g g u g g u g g g g g u g u g g u g g u g g g g g u g u g g u g 2 1,1 ,2 , + ⋅ g u (54) where products of displacement gradients are neglected.Using the equations ( 26) and (51), we obtain: ,11 ) The relations ( 53) and (56) represent the reference strains of the beam axis as functions of reference geometry, displacement gradients of the beam axis, the total rotation of a cross section around g1 basis and gradient of this rotation.

Strain at an arbitrary point
In order to distinguish stress and strain components with respect to the two considered coordinate systems, components with respect to the (ξ, η, ζ ) coordinate system are designated with tilde while the ones with respect to the (ξλ, η, ζ ) system are designated with an overbar.The normal strain at an arbitrary point, with respect to the (ξ, η, ζ) coordinate system, is: where the component of the metric tensor in the reference configuration is defined with (33) while its counterpart in the deformed configuration is, analogously: Introduction of these two terms into (57) returns: Likewise, the shear strains due to torsion, at an arbitrary point of a cross section, are: Remark.The well-known problem of coupling of torsional and normal strain in (59) exists due to initial torsion of the beam.This problem is treated differently in literature.One of the most consistent approaches is given in [4], where the authors utilized the so-called conjugated stress and strain.
Components of these tensors are defined with respect to the base vectors of the centroid and an arbitrary point of a cross section.Although rather unconventional for beam theories, this approach allowed authors to derive a correct value of virtual work, [4].Here, the newly introduced coordinate system (ξλ, η, ζ) returns the same value of virtual work, but in (art.a) more elegant manner.It should be pointed out that the conditions of orthogonality (38) are actually contained in the transformation of base vectors from the centroid to an arbitrary point in [4], but they are not explicitly stated.
The normal strain at an arbitrary point with respect to the coordinate system (ξλ, η, ζ) is obtained using the tensor property: If the relations (59) and ( 60) are introduced into (61), we obtain the final expression: ( ) Contrary to the strain ε11(ξ, η, ζ), the strain ε11(ξλ, η, ζ) is expressed with respect to the same base vectors over the whole cross section (gλ=g1).The intensity of this base vector is arbitrary, but its orthogonality with respect to the cross-sectional plane is quintessential.As already mentioned, the most significant result of the introduction of the new coordinate system is that the torsional strain κ1 is rigorously excluded from the equation ( 62) which means that the normal strain is independent of torsion with respect to the coordinate system (ξλ, η, ζ).
Remark.Let us introduce curvature changes with respect to the Frenet-Serret frame of reference: and recall that the flexural strains can be represented as (see [26]): Inserting (65) into (62) results with the classic expression for the normal strain with respect to the Frenet-Serret coordinates: ( ) ( ) where the indices in parentheses indicate physical components of strain.This discussion is analogous to the one given in [26] for the plane case.However, for the spatial case, it is important to emphasize that the equation ( 66) is only valid with respect to some sλ coordinate, which can be introduced analogously to (38).This is due to fact that the standard arc-length coordinate s is orthogonal with respect to the cross-sectional plane only at the centroid.

Alternative formulation
The total angle of rotation of a cross section around tangential basis g1, ω, is composed of two parts, Fig. 4. The first one is the rotation of normal plane, ϖ, while the second one is the change of angle between the normal and binormal and principal axes of inertia, ∆θ: where the rotation ∆θ is described analogously to (50): Since the rotation of normal plane: is a function of displacement field of the beam axis, it follows that we can adopt the change of angle between the Frenet-Serret triad and principal axes of inertia as an alternative generalized coordinate.A detailed derivation of displacement gradients with respect to the normal and the binormal is given in the Appendix.
From (56), ( 22) and ( 23) it follows that: ( ) ( ) which can be written as: where the following designations are introduced: while the derivatives of the Christoffel symbol and the modulus of curvature are: .
Remark.It is interesting to show that an identical result can be obtained directly from the equations ( 23) and ( 54): where the angle between (n * , b * ) and (g2 * , g3 * ) is: It is now necessary to find the difference of torsions of the beam axis in the deformed and reference configuration in the equation (74).We will perform this by approximating torsion of the deformed axis with the first two terms of the Taylor series: ( ) After some algebraic manipulation, we can find: ,111 , which, after inserting into (74), returns (71).
Consequently, the expressions for flexural strains (56), as functions of alternative generalized coordinate, reduce to: .
This alternative formulation is more complex than the standard one because it requires at least C 2 interelement continuity, due to the presence of displacement gradients of the third order.However, it has profound theoretical significance since it demonstrates potential to introduce the change of angle between the Frenet-Serret frame and principal axes of inertia as a generalized coordinate.It should be emphasized that this formulation actually represents a rotation-free model of the spatial BE beam if the nature of a problem is such that influences of terms containing ∆θ and ∆θ,1 in (71) and (78) are negligible.

Finite element formulation
In this section, a detailed derivation of stress resultants energetically conjugated with the reference strains of the beam axis is presented.It allowed us to rigorously define internal and external virtual work, which ultimately lead to a formulation of appropriate stiffness matrix and external load vector of an isogeometric finite element of the spatial BE beam.Rotation boundary conditions for both formulations are briefly discussed at the end.

Stress-strain relation
For linear elastic material considered here, the generalized Hooke's law is applied: where µ and λ are the Lame's parameters: while E is the modulus of elasticity and ν is Poisson's ratio.i j σ are mixed components of the second Piola-Kirchhoff stress tensor with respect to the reference configuration.Using the condition that components 2 2 σ and 3 3 σ are zero we obtain: ( ) Accordingly, non-zero contravariant components of stress are: ( ) This relation is valid for general beam element and it is independent of coordinate system.With respect to the (ξλ, η, ζ ) coordinates, three non-zero stress components for the BE beam are: If we recall the reciprocal metric tensor (34), it follows that, with respect to the (ξ, η, ζ ) coordinates, 2 2 σ and 3 3 σ are non-zero at an arbitrary point of a cross section.In addition, it is evident that usage of (ξλ, η, ζ ) coordinates decouples normal and torsional strains, not only for the BE beam but also for the Timoshenko beam, contrary to the (ξ, η, ζ ) coordinate system.
It is worthwhile to show relations of these three stress components with respect to both local coordinate systems:

Stress resultants
In order to perform the integration over a cross section, it is required to recall that the differential volume element, for both coordinate systems, is: , where εijk is the Levi-Civita tensor.Now, the only section force that explicitly follows from the BE theory is: It is evident from (84) and (86) that this section force is identical with respect to both coordinate systems, which is expected since the ξ and ξλ coordinates coincide at the centroid.Using the equations (62), ( 83) and (86), we obtain: which can be written as: ( ) ( ) where the following geometric properties are introduced: Resultant moment of normal stress with respect to the centroid is: where: M1 is the torsional moment, while M2 and M3 are bending moments with respect to g2 and g3, respectively: ( ) ( ) Again, the section moments are equal for both coordinate systems, analogously to the normal force.
Geometric property labeled as I11 reduces to the second polar moment of area for the straight beam.It is well-known fact that it is equal to the torsion constant for circular cross sections only, [33].Therefore, a correction factor is applied here for examples regarding beams with rectangular cross sections.The property I11 is scaled according to the relation of the polar moment of area and torsion constant, [33].A similar scaling is required for all shapes of a cross section, except for circular one.
If vectors of the generalized section forces and reference strains are introduced: then we can write: where the constitutive matrix is: It should be pointed out that this matrix is not symmetric because, for curved beams, the section forces do not produce virtual work on the reference strains but their energetically conjugated counterparts.

Internal virtual work -stiffness matrix
Virtual work produced by internal forces, with respect to both coordinate systems, is: , , , , , .
It can be easily shown that utilization of stresses and strains, with respect to both coordinate systems, gives the same scalar value of internal virtual work, which validates the introduction of ξλ coordinate.
Therefore, after the integration over cross section area we obtain: ( ) where stress resultants, energetically conjugated with the reference strains, are introduced: ( ) ( ) Final expressions for these stress resultants are obtained by inserting (83) into (100): Furthermore, we will define vector: which allows us to write: where the constitutive matrix is: while the additional geometric property is introduced as: Observe that the constitutive matrix D, contrary to D , is symmetric because it couples energetically conjugated variables.Now, the stiffness matrix can be defined using the following vectors and matrices: These designations enable us to represent the internal component of virtual work as: where: is the linear stiffness matrix of a spatial BE isogeometric finite element for which the rotation around tangential basis is parameterized with the total angle of twist of a cross section, ω.The alternative formulation, using ∆θ as a generalized coordinate, returns the following vectors and matrices: where the submatrix W is the same as in (106).
Remark.The present derivations of metric of the spatial beam and its finite element formulation represent the general case of the model developed in [26].By introduction of restrictions ω=θ=θ * =0, the present spatial formulation reduces to the plane case given in [26].

External virtual work -vector of external load
Virtual work of external forces is equal to negative work of external loads on virtual displacements.Volume and surface loads are reduced to the beam axis: where the vectors: , , consist of the global rectangular Cartesian components of a distributed load and distributed moment.
Influence of their concentrated counterparts can be considered as a special case.
If we introduce a vector of components of distributed loads and moments: and utilize the expressions for infinitesimal rotations, (47), the transformation matrix is obtained: , Now, the external part of virtual work reduces to: where, excluding the effect of non-homogenous boundary conditions, the external load vector is: For the alternative formulation, the fourth row of the submatrix NI (113) must be modified in order to obtain the complete angle of twist of a cross section:

Boundary conditions
Boundary conditions are here expressed via the global rectangular Cartesian coordinates, just as all quantities, except the angle of twist.A detailed discussion on displacement and rotation boundary conditions is given in [26].However, the imposition of rotation boundary conditions is somewhat different for the spatial case.Concretely, the clamped boundary condition is prescribed by the elimination of rotations: 2 For the fixed value ξ=ξ0, the relations (117) represent two constraint equations which must be satisfied and therefore lead to the reduction of a number of DOFs by two.Since the clamped boundary conditions are frequently applied at the ξ=0 or ξ=1, only the first or the last two control points have an influence on rotation at these coordinates.Seeing that a control polygon is tangent to a curve at these coordinates, it follows from the condition (117) that the displacement components of these control points, perpendicular to the tangent polygon and collinear with the principal axes of inertia at the clamped section, must be equal.To the best of our knowledge, these conditions are applied for the first time to an arbitrary directed spatial beam axis in this paper, concretely in the second example in the following section.
Restraining the angle of twist is trivial for the standard formulation where this angle is a generalized coordinate.However, for the alternative formulation, it requires an additional effort, since: Concretely, for the boundary condition ω(0)=0, the constraint equation is: where α, β, γ and ψ are functions of first and second derivatives of basis functions at ξ=0, while indices 1,2 and 3 refer to serial numbers of control points.

Numerical analysis
In order to verify the presented formulation, two numerical examples are studied.The first one is taken from literature while the other one is conceived by the authors.In addition, all presented results are calculated using a code built in the software package Wolfram Mathematica.This package provides an efficient programming environment and allows implementation of attractive data visualization capabilities, [34].
Although it is well-known fact that the Gauss quadrature rule is not optimal for IGA, due to increased continuity, [35], numerical integration is here performed with the classic full Gauss quadrature rule with p+1 integration points.Additionally, no explicit analysis of membrane and torsional locking phenomena is performed here.However, similar effects of membrane locking as those for plane beams are expected, [26].The torsional locking deserves a special attention and it will be tackled in forthcoming research.
Regarding the two developed formulations with respect to the parameterization of the angle of twist, it should be stressed out that both of them return virtually indistinguishable results for the selected examples.For clarity, the additive decomposition of the total angle of twist (67) is represented as a diagram along the beam axis in both examples.

Convergence analysis
The convergence properties are the ones of the most prominent features of IGA and they are examined thoroughly.Since the order of convergence should be independent of interelement continuity, [13], an improved accuracy per DOF is expected for high continuity IGA, compared to the standard FEA.There are many researches devoted to this problem.For example, posteriori error estimate analysis performed in [36] suggested that the higher order continuities are promising for smooth problems.However, authors in [37] analyzed the error estimates dependent on the hpk-refinements and they observed that the convergence with respect to the NURBS order is better for C 1 continuity, compared to the higher ones.Nevertheless, they used maximum continuity of C (p-1)/2 and left the question of higher regularities opened.Therefore, the presented examples are utilized for a test of the optimal order of convergence, which is p+1 for displacements, p for normal force and torsional moment and p-1 for bending moments, [38] and [18].Since these orders vary nonlinearly under the h-refinement, a better insight is enabled by representing the order of convergence as a function of number of elements nel, similar as in [27].A uniform mesh in the parametric space is utilized and the order of convergence, for a mesh with nel elements, is calculated as: where X ex is the reference solution while h m X is an approximate solution for a mesh with m elements.
For examples where there is no analytical solution available, results obtained with very dense meshes are utilized as the reference values.During the analysis, it is found that the order of convergence oscillates sharply for ranges where the convergence of the observed relative error is nearly reached or the error itself gets very small.In order to obtain readable graphics, some oscillating parts of these graphs are truncated.A detailed discussion regarding this behavior is given in [27].

A pre-twisted circular beam
The first example deals with a pre-twisted circular beam given in Fig. 5.It is loaded with a uniformly distributed load pz=-1 kN/m.All translations and the total angle of twist are restrained at the end sections.The most interesting property of this example is the presence of initial geometric torsion of the beam while its axis is restricted to the plane.The beam is taken from the paper by Greco and Cuomo [18] where three cases are considered: positive, zero and negative twist.However, here we study only the case that they designated as the negative twist, Fig. 5.It should be noted that the authors in the paper [18] stated different dimensions of the cross section than those used here, but through a private communication, they confirmed that the results reported in [18] are in accordance with the properties given in Fig. 5.
Moreover, the simple straight spatial B33 element is utilized and an appropriate FE model is made, which is not a trivial task, [18].A special attention is paid to the correct definition of section normals and a dense mesh of 314 B33 elements is created in Abaqus, [39].
Convergence properties for the displacement magnitude at the middle section, as well as for stress resultants at the start section are presented in Fig. 6, Fig. 7, Fig. 8, Fig. 9 and Fig. 10.It is evident from Fig. 6b that the order of convergence for p=2 is lower than the expected value of p+1.For the cubic NURBS, these orders are converging to the analytically predicted value.It is interesting to observe that the order of convergence for the model with C 1 continuity converges from below, while for the C 2 model it converges from above.This behavior is generally found, and it is valid for most variables considered in this research.It returns an improved accuracy per DOF for the models with the highest continuity.For the quartic NURBS, the behavior of the order of convergence is oscillatory, but it can be remarked that the model with C 1 continuity returns higher values, compared to the C 2 ones.However, for the model with the highest, C 3 continuity, the order of convergence is extremely high.All these observations are appropriately reflected at the convergence of relative error in Fig. 6a.
Next, convergence properties are somewhat different for the normal force, Fig. 7. Firstly, the quadratic NURBS shows an unexpectedly high order of convergence, which is degrading for denser meshes.However, its accuracy is very low, which is expected.The orders of convergence for all the other models are close to the predicted value of p, except for the quartic NURBS with the highest continuity, which has an extremely good convergence.For the torsional moment in Fig. 8, quadratics and cubics return the orders of convergence that are close to their respective degrees p.However, the convergence is very fast for the quartics and orders could not be detected.
Regarding the bending moment M2, the orders of convergence for all models converge to the analytical values, Fig. 9b.Convergence analysis of the bending moment M3 in Fig. 10 reveals that the quadratics behave according to predictions, while the cubics return somewhat a better order of convergence than the analytical value of p-1=2.Regarding the quartic NURBS, the relative error for the models with C 1 and C 2 continuity converges with an extremely high order.Nevertheless, the orders of convergence are close to the theoretical prediction for dense meshes for the model with C 3 continuity.
Comparison of kinematic and static variables in Fig. 11, Fig. 12, Fig. 13 and Fig. 14 shows an excellent compliance with the results from the reference [18], the B33 model and the present approach.Still, small discrepancies exist and the most prominent ones are detected for the distribution of normal force in Fig. 13a, where the results from the reference [18] slightly deviate from the results obtained with the present approach and B33 element.It should be noted that the authors in [18] used an approximate coefficient for torsion constant, 1/3, while in this research a more accurate value of 0.312 for b/h=1/10 is used, [33].In the end, the scaled deformed configuration is displayed in Fig. 15.
Interestingly, the additive decomposition of twist angle in Fig. 12b reveals the fact that influence of rotation of the normal plane is dominant along the nearly complete length of the beam.Nevertheless, the change of angle between the Frenet-Serret frame and the principal axes of inertia is of the same order of magnitude.It is noticed that the twist angle of the cross section has an extremum at ξ=0.55, while ∆θ changes sign in the vicinity of ξ=0.43.

A cubic beam
The second example considers a spatial cubic beam represented in Fig. 16.The beam is clamped at the lower end and loaded with a vertical concentrated force at the upper end, F=1000 kN.In order to make a detailed comparison of results, an appropriate model is made in Abaqus using 618 B33 elements.Since the beam is statically determined, stress resultants at the clamped end are calculated analytically.
A detailed analysis of convergence properties is given in Fig. 17, Fig. 18, Fig. 19, Fig. 20 and Fig. 21.Contrary to the other graphics, the orders of convergence are estimated in Fig. 17 and Fig. 19, since their behavior is too oscillatory for the presentation.Fig. 17 shows that the cubic elements behave according to the predictions.However, an extremely high of convergence is noticed for the quartics and quintics with the highest continuity.The C 1 and C 2 quartics are close to the expectations, while the quintics, lower than C 4 , show a much higher order of convergence than the predicted value of p+1=6.
Regarding the convergence of the normal force, Fig. 18b reveals that the orders of convergence for the cubic NURBS are close to the expected value, while for the C 1 and C 2 quartics this order does not converge.On the other hand, the order of convergence for the C 3 quartic element is close to the theoretical value of p=4.Considering the quintics, it is evident that higher continuities return higher orders of convergence.It is interesting to note that the C 2 cubic element yields more accurate solutions than C 1 , for the same number of elements.
A similar behavior is observed in Fig. 19 for the torsional moment.Estimated orders of convergence are close to the theoretical predictions while the relative error shows a non-monotonical behavior for the quintic NURBS.
For the bending moment M2 in Fig. 20b, all polynomials return the orders of convergence close to the predictions.Contrary to some previous findings, here it is evident that lower continuities give higher values of the order of convergence.Nevertheless, Fig. 20a reveals that the accuracy of models with a higher interelement continuity is improved.In Fig. 21b an extremely high order of convergence for M3 is detected for the cubic NURBS.Moreover, the quartics behave according to the predictions, while the quintics deserve a special attention.Namely, models with all continuities converge to the expected value of p-1=4.However, it is evident that the order of convergence for the model with C 3 continuity converges from below, while for the model with C 4 continuity converges from above, which results with an increased accuracy per DOF for the model with the highest continuity.
Distributions of kinematic and static variables in Fig. 22, Fig. 23 and Fig. 24 show that the results obtained with the present formulation and those from B33 element are virtually indistinguishable.Regarding the decomposition of twist angle in Fig. 24a, it is interesting to note that the rotation of normal plane prevails, similar as for the pre-twisted beam.The extremum is at the free end where the angle ∆θ is zero and the rotation of the cross section equals that of the normal plane, ω=ϖ.The scaled deformed configuration is presented in Fig. 25 from two points of view.
In the end, comparison of displacement components at the free end for nine different values of curviness is performed.Besides the classic models, composed of one-dimensional B33 and B31 elements, the C3D10 tetrahedral solid element is utilized and full three-dimensional models are created in Abaqus, [39].The aspect ratio of the cross section remained fixed and the whole section is scaled for all nine models.Since there is no initial twist in this example, the maximum dimension of the cross section parallel to the osculating plane, hmax, reduces to b, Fig. 16.Displacement components at the free end for B33, B31, C3D10 and the present models are given in Table 1 for nine different values of Kmaxb.It should be remarked that the maximum curvature of the considered beam is Kmax≈1.48.The graphical representation of relative differences of displacement components with respect to the results from the C3D10 model is given in Fig. 26.
A good agreement of the results for all models is noticed for lower values of curviness.However, as Kmaxb increases, the relative differences increase, as well.It is evident from Table 1 that the present formulation returns the stiffest response of all considered models.This is due to two important features of the present model.The first one is the rigorous inclusion of the complete metric of beam continuum, contrary to the straight elements B33 and B31.The second one is the strict application of the BE assumptions.Since they are inevitably erroneous for strongly curved beams, the stiffer structural response is obtained, compared to the general solid model without kinematic restrictions.Additionally, the effect of torsional locking could be influential, and it will be examined in the upcoming paper.
There are several interesting observations regarding the results displayed in Fig. 26.Firstly, the relative difference between the beam and solid models is detected for displacement component along the y-axis.On the other hand, this component has the lowest magnitude of all three, which proportionally affects the relative difference of displacement vector.Furthermore, each considered beam model gives the best compliance with the C3D10 element for exactly one of the displacement components.These are ux, uy and uz for the present, B31, and B33 models, respectively.Finally, the most unexpected result from this analysis is an extremely good compliance of uz displacement component obtained with the simple straight B33 element, which is even better than that of the shear deformable B31 element.This brief analysis suggests that the present example is complex and a further comparison of theories and elements is required in order to draw firm conclusions.
A satisfactory level of accuracy, as well as the limits of applicability of the present IGA BE beam formulation, are problem dependent.Nevertheless, the results in Fig. 26 suggest that we can expect a reasonable compliance with a general solid model for strongly curved beams with Kmaxb<0.3.In order to depict considered geometries more clearly, models with three high values of curviness are displayed in Fig. 27 where it is evident that these strongly curved beams hardly resemble the classic beam geometry.

Conclusions
Linear analysis of an arbitrarily curved spatial BE beam is successfully performed.Strains at an arbitrary point are rigorously defined via reference strains of the beam axis.The introduction of the new orthogonal coordinate system enabled decoupling of the normal strain at an arbitrary point and the torsional strain of the beam axis.This led to a simple and exact definition of stress resultants that are energetically conjugated with the reference strains.Elegant expressions are obtained, being a function of Cartesian components only.
However, the angle of twist is an exception.It is the only quantity of the spatial BE beam that must be expressed with respect to the local coordinates.The introduction of the novel parameterization of this variable proved valid and an additional insight into its nature is enabled.
The convergence analysis confirmed previous findings that models with the highest interelement continuity return an improved accuracy per DOF for the considered problems.However, due to the highly nonlinear behavior of the order of convergence with respect to the h-refinement, a more detailed analysis is needed.
The numerical study proved that the present approach is applicable for strongly curved beams.However, in order to define firm limits of applicability of the developed theory, it is required to make more comparisons with different analytical, numerical and physical models.Further research will be focused on nonlinear analysis, for which the present formulation is ideally suited.

Fig. 3 .
Fig. 3. Reference configuration of an arbitrarily curved spatial BE beam; Coordinate axes and respective base vectors.

Fig. 4 .
Fig.4.Rotation of a cross section around tangential basis g1 -twist.Reference and deformed configurations of a rectangular cross section are designated with green and red lines, respectively.

Fig. 6 .
Fig.6.The pre-twisted beam.Convergence properties for the displacement magnitude at the middle section using the h-refinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 7 .
Fig.7.The pre-twisted beam.Convergence properties for the normal force at the start section using the hrefinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 8 .
Fig.8.The pre-twisted beam.Convergence properties for the torsional moment at the start section using the hrefinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 9 .
Fig.9.The pre-twisted beam.Convergence properties for the bending moment M2 at the start section using the hrefinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .
Fig. 10.The pre-twisted beam.Convergence properties the bending moment M3 at the start section using the h-refinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 15 .
Fig. 15.The pre-twisted beam.Reference (transparent) and deformed configurations.Displacements are scaled with a factor of 500.

Fig. 17
Fig. 17.The cubic beam.Relative error for the displacement magnitude at the free end the h-refinement with different NURBS orders and continuities.The orders of convergence are estimated.
Fig. 17.The cubic beam.Relative error for the displacement magnitude at the free end the h-refinement with different NURBS orders and continuities.The orders of convergence are estimated.

Fig. 18 .
Fig.18.The cubic beam.Convergence properties for the normal force at the fixed section using the h-refinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 19 .
Fig.19.The cubic beam.Relative error for the torsional moment at the fixed section using the h-refinement with different NURBS orders and continuities.The orders of convergence are estimated.

Fig. 20 .
Fig.20.The cubic beam.Convergence properties for the bending moment M2 at the fixed section using the hrefinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 21 .Fig. 22 .
Fig. 21.The cubic beam.Convergence properties for the bending moment M3 at the fixed section using the hrefinement with different NURBS orders and continuities: a) relative error; b) order of convergence.

Fig. 24 .
Fig. 24.The cubic beam.Distributions of: a) rotations around basis g1; b) displacement components and displacement magnitude u.

Fig. 25 .
Fig. 25.The cubic beam.Reference (transparent) and deformed configurations from two different points of view.Displacements are scaled with a factor of 2.

Fig. 26 .
Fig. 26.Relative differences of the displacement components at the free end vs. curviness Kmaxb.

Fig. 27 .
Fig. 27.Reference configurations of the cubic beam for three high values of curviness.