Geometrically exact static isogeometric analysis of arbitrarily curved plane Bernoulli-Euler beam

We present a geometrically exact nonlinear analysis of elastic in-plane beams in the context of finite but small strain theory. The formulation utilizes the full beam metric and obtains the complete analytic elastic constitutive model by employing the exact relation between the reference and equidistant strains. Thus, we account for the nonlinear strain distribution over the thickness of a beam. In addition to the full analytical constitutive model, four simplified ones are presented. Their comparison provides a thorough examination of the influence of a beam's metric on the structural response. We show that the appropriate formulation depends on the curviness of a beam at all configurations. Furthermore, the nonlinear distribution of strain along the thickness of strongly curved beams must be considered to obtain a complete and accurate response.


Introduction
Curved beams are fundamental structural components and their analysis is a classic topic of mechanics [1]. Although some analytical solutions exist, see, e.g., [2], numerical methods are inevitable for dealing with the complex nonlinear behavior of beams. The finite element method (FEM) is the most versatile numerical procedure for solving partial differential equations. The first research on the geometrically-exact FEM analysis of curved spatial beams was presented in a seminal paper by Reissner [3]. This theory has been enhanced in a series of papers [4,5,6,7,8]. In [4], Simo introduced the term geometrically exact beam theory to designate that the governing equations of motion are valid regardless of the magnitude of kinematic quantities. The papers referred to are mainly focused on the shear-deformable beam model, and it is not until recently that the Bernoulli-Euler (BE) spatial model was scrutinized [9,10,11]. One reason for this is that the equations of the Simo-Reissner beam require only C 0 interelement continuity. On the other hand, the BE model requires C 1 interelement continuity -a requirement not fulfilled by standard Lagrange polynomials commonly used in FEM. Herein, we resolve this issue by using isogeometric analysis (IGA), [12,13]. This approach utilizes (smooth) splines as basis functions for the spatial discretization of the weak form of the equations of motion, and hence, it allows for an arbitrarily high continuity between elements.

Metric of the beam continuum
A detailed and rigorous definition of the beam metric is presented in this section. Following the classical BE assumption, a cross section is rigid and remains perpendicular to the beam axis in the deformed configuration. This assumption leads to a degeneration of a 3D continuum beam model into an arbitrarily shaped line. The present analysis is performed with respect to the convective frame of reference while the complete beam kinematics is defined by the translation of the beam axis.
In the notation, lowercase and uppercase boldface letters are used for vectors and tensors or matrices, respectively. An overbar designates quantities at the equidistant line of the beam, and the asterisk sign denotes the deformed configuration. Finally, the hat symbol specifies the local component of a vector, with respect to the curvilinear coordinates. The direct and index notations are applied simultaneously, depending on the context, and the standard summation convention is adopted. Greek index letters take values of 1 and 2. Partial and covariant derivatives with respect to the convective coordinates are designated with () ,α and () |α , respectively. The elaboration on the NURBS-based IGA modeling of curves is excluded for brevity since it is readily available in the literature. For a detailed discussion on IGA and NURBS, references [12,36,37] are recommended.

Metric of the beam axis
Let us revise some basic expressions of the metric of beam axis. A beam axis is a curve that passes through the centroids of all cross sections. It can be defined with either the arc-length coordinate s or some parametric coordinate ξ . The position vector of the beam axis in Cartesian coordinates is r = {x = x 1 , y = x 2 } and it is here defined as a linear combination: where R I are univariate NURBS basis functions, r I = {x I , y I } are the position vectors of the control points I, and N is the total number of control points [12]. Furthermore, i α = i α are the base vectors of the Cartesian coordinate system, see Fig. 1. The tangent base vector of a beam axis is: The other base vector, g 2 , is perpendicular to g 1 and has unit length. Hence, we can choose the normal of the beam axis, as in [28], or we can rotate g 1 by 90 degrees and scale it. The latter approach is utilized in present formulation and the base vector g 1 is rotated anti-clockwise: Here, g is a component of the metric tensor and also its determinant: This quantity relates differentials of arc-length and parametric coordinate as ds = √ g dξ . Therefore, it equals the square of Jacobian of the coordinate transformation from convective to arc-length coordinate.
The metric of the beam axis is completely defined by the introduction of the Christoffel symbols. They follow from the differentiation of base vectors with respect to the curvilinear coordinates [38]: where Γ γ αβ are the Christoffel symbols of the second kind. The reciprocal base vectors of the beam axis are: and the reciprocal metric tensor is: The expression for the derivatives of base vectors in matrix form is: where K is so-called signed curvature: andK = gK is the signed curvature of beam axis with respect to the convective coordinate frame. It should be noted that, since we are dealing with both the signed curvature and the modulus of curvature, the curviness Kh itself can be defined with either of them. In the remainder of the paper, Kh usually refers to the curviness obtained with the modulus of curvature. It is only in Subsection 5.3 that the so-called signed curviness is used, merely for graphical representation.

Metric of an equidistant line
The degeneration from a 3D to a 2D beam model is trivial since we are dealing with plane beams, and all quantities are constant along the ζ coordinate, as illustrated in Fig. 1. In order to reduce the beam continuum from 2D to 1D, the metric of the complete beam must be defined by some reference quantities. It is common to use the metric of the beam axis as the reference. Therefore, let us define an equidistant line which is a set of points with η = const. Its position and tangent base vectors are: Evidently, the base vector of this line is parallel to the base vector of beam axis. Their lengths differ for the initial curvature correction term g 0 [27]. The other base vector of the equidistant line is the same as the one of beam axisḡ 2 = g 2 . The metric tensor of an equidistant line is: Note that the initial curvature correction term g 0 becomes zero for ηK = 1 which is physically prohibited. Since η ∈ [−h/2, h/2], it follows that the upper bound of curviness is 2.

Bernoulli-Euler beam theory
After the metric of the beam continuum is defined, the next step is to introduce a strain measure. For the convective coordinate frame, the Lagrange strain equals the difference between the current and reference metrics: Due to the BE hypothesis, the shear strain vanishes and the position of the cross section is completely determined by the components of translation of the beam axis. In other words, these components are the only generalized coordinates of the BE beam. This fact gives rise to the so-called rotation-free beam theories [28].

Kinematics
In the deformed configuration, the position vector of an equidistant line is: where r * denotes the position vector of the deformed beam axis and the base vector g * 2 is calculated similarly as in Eq. (3). The position vector of the deformed beam axis is: where u is the displacement vector of the beam axis. It is discretized with NURBS, in the same way as the geometry, Eq. (1): u I is the vector of displacement components of a control point I with respect to the Cartesian system. For the sake of a seamless transition to the discrete equation of motion which follows in the Section 4, let us introduce the matrix of basis functions N such that Eq. (15) can be written as: where: The material derivative of the displacement field yields the velocity field, i.e.: v =ṙ * =u =u n i n = v n i n =ẋ * n i n , while its spatially discretized form follows from (15): v =u = Nq. The work conjugate pair for the Cauchy stress tensor is the strain rate tensor. It is the symmetric part of the velocity gradient and equals the material derivative of (12). The covariant components of the strain rate tensor, with respect to the Cartesian coordinates, can be written as: wherev α are the components of the velocity with respect to the local curvilinear coordinates [39]. For BE beams, the only non-zero component of the strain rate is d 11 : In order to define the kinematics of a beam continuum, the velocity of an equidistant line must be introduced. It is the material derivative of (13): where v ,2 follows from the definition of shear strain: For the purpose of further derivation, let us express the components of this gradient by index notation, [28]:

Strain rate at an equidistant line
The strain rate at an equidistant line is, analogously to Eq. (20): For the evaluation of this expression, let us first define the material derivative of mixed velocity gradient, using Eq. (8): With this relation and Eqs. (10), (21) and (24), the equidistant strain rate reduces to: At this point, we will introduce a curvature change of beam axis with respect to the convective coordinate: and its material derivative, the rate of curvature change: κ =ġ * K * + g * K * = 2d 11 K * + g * K * ⇒ g * K * =κ − 2d 11 K * .
By inserting the last expression into Eq. (26), the final form of the equidistant strain rate is obtained: This expression is analogous to the one derived in [28] for the equidistant axial strain. It is valid for finite, but small, strain analysis, which falls into the scope of the geometrical exact beam theory [6]. Finally, representing the rate of curvature change as a function of velocity gradients of the beam axis is a requirement. It follows from Eqs. (8), (22), and (27):

Finite element formulation
In line with the previous derivation, we will formulate the isogeometric BE element using the principle of virtual power. The generalized coordinates are the components of the velocities of the control points. We start from the generalized Hooke law for linear elastic material, also known as the Saint Venant-Kirchhoff material model. This material model is well-suited for the small strain and large rotation analysis [40]. The only non-zero components of the stress and strain rates are related as: where E is the Young's modulus of elasticity. It is a well-known fact that this axial state of stress is an erroneous consequence of the BE assumptions. Nevertheless, it readily used due to its simplicity.

Principle of virtual power
The principle of virtual power represents a weak form of the equilibrium. It states that at any instance of time, the total power of the external, internal and inertial forces is zero for any admissible virtual state of motion. If the inertial effects are neglected, and body and surface loads are reduced to the beam axis, it can be written for plane BE beams as: where σ σ σ is the Cauchy stress tensor, d is the strain rate tensor, and p is the vector of external line loads. All these quantities are measured with respect to the current, unknown, configuration. For problems with deformation-independent loads, the only requirement is to linearize the Cauchy stress. At the current (n + 1) configuration this stress is approximated by: where (n) (n+1) σ σ σ is the stress from the previous (n) configuration expressed with respect to the metric of the current (n + 1) configuration.
Remark. The additive decomposition of stress, as in Eq. (33), is only valid if all of the terms are expressed with respect to the same metric. Therefore, the stress from the previous configuration is here adjusted to the current metric by an adjustment term. This term is approximately equal to the ratio of determinants of the metric tensors in two configurations, c.f. Appendix A. The exact integration of the adjusted stress is impractical because the adjustment term consists of the ratio of the initial curvature correction terms at two configurations. Three approaches are considered here for dealing with this issue. These are discussed in detail in Appendix A and compared in Subsection 5.1.4.
An appropriate time derivative in Eq. (33) is designated with d t ()/ dt while ∆t (n+1) = t (n+1) −t (n) is the time increment. The rate form of stress-strain relations requires an objective time derivative. Note that the material derivative is not objective, while the corotational (Jaumann) and convective time derivatives fulfill this property [41,42]. The convective derivative of stress tensor results in the stress rate tensor. An important fact is that the components of the stress rate tensor are equal to the material derivatives of the components of the stress tensor, [39]. Having this in mind, and after the insertion of Eq. (33) into Eq. (32), the linearized form of the principle of virtual power at the current configuration is obtained: In the remainder of this paper, we neglect the time indices and asterisks for the sake of readability. Note that this simplification does not introduce any notational ambiguity since (i) the stress and strain rates are instantaneous quantities, while the known stress is calculated at the previous configuration, and (ii) all integrations are performed with respect to the metric of the current configuration, in accordance with the updated Lagrangian procedure [43]. In order to reduce the dimension from 3D to 1D, it is necessary to integrate the left-hand side of Eq. (34) over the area of the cross section. Thus, integrals along the length of the beam axis are obtained: whereÑ andM are the stress resultant and the stress couple, which are energetically conjugated with the reference strain rates of the beam axis, d 11 , andκ , whileṄ andṀ are their respective rates: These expressions are analogous to those obtained in [28]. By introducing the vectors of generalized section forces, strain rates of the beam axis, and external line loads: Eq. (34) can be written in compact matrix form as: This equation is nonlinear and it will be linearized in Section 4.4.

Relation between energetically conjugated pairs
The geometrically exact relations (35) and (36) are crucial for the accurate formulation of structural beam theories. In particular, they allow a rigorous definition of energetically conjugated pairs, and guarantee that the appropriate constitutive matrix is symmetric. By the introduction of Eqs. (29) and (31) into Eq. (36), we obtain the exact relation between energetically conjugated pairs of stress and strain rates. The resulting symmetric constitutive matrix D is derived in [28]: where: Importantly, these integrals can be analytically determined for the conventional solid cross section shapes [28]. In the following, we refer to the exact constitutive model given in Eq. (39) by D a . We introduce four reduced models to examine the various influences of the exact constitutive relation. The first and the simplest of these is designated with D 0 . It is often employed for thin beams, e.g. [25,19], since it decouples axial and bending actions: The second reduced model, D 1 , is based on the approximation: g 0 → 1. It is readily utilized for the analysis of curved beams with small curvature [9,32]. The model is specified by: where the quadratic term in the integrand of property A 1 is disregarded.
Finally, the models D 2 and D 3 are based on the Taylor approximation of the exact expressions (40): where the 4 th moment of area can be easily calculated for rectangular and circular cross sections.

Variation of strains
Since the strain rate is a function of the generalized coordinates as well as the metric, we must vary it with respect to both arguments. By noting that the variation of the tangent vectors can be expressed as: the variations of reference strains are given by: Most of the terms in the previous expression are easily computed since the variation is explicitly performed with respect to the unknown variables. However, the first two addends ofκ are exceptions since they require the variations of the velocity of the normal and the Christoffel symbol. These variations must be represented via the variations of generalized coordinates [33]. The variation of the velocity of the normal follows directly from Eq. (22): For the variation of the Christoffel symbols, we start from Eq. (5) and, after some straightforward calculation, obtain: By inserting Eqs. (47) and (48) into Eq. (46), the final expression for the variation of the rate of curvature change is found:

Discrete equation of motion
In this subsection, the virtual power is spatially discretized and linearized. We start by introducing the matrix B L , which relates the reference strain rates of the beam axis with the velocities of control points: Using Eqs. (20) and (30), the vector of the reference strain rates can be represented as: where the vector w is defined by: The submatrices B α I for an arbitrary control point I consist of the derivatives of the basis functions, i.e.: The matrix B L now follows as: Next, we need to specify the explicit matrix form of the virtual power generated by the known stress and the variation of strain rate with respect to the metric, Eqs. (38) and (46): where the matrix G β α is: with: By introducing the total matrix of the generalized section forces: expression (55) can be written as: A careful inspection of Eqs. (56) and (58) reveals that the matrix of generalized section forces, G, is symmetric [33]. Now, the integrands in the equation of the virtual power (38) reduce to: where the first term is linearized by neglecting the variation of strain rate with respect to the metric: Finally, the equation of equilibrium reduces to: which is readily written in the standard form: where: is the tangent stiffness matrix and: are the vectors of the external and internal forces, respectively. The vector ∆q in Eq. (63) contains increments of the displacements. Due to the approximation introduced in Eq. (61), the solution of Eq. (63) does not satisfy the principle of virtual power directly, but has to be enforced by additional numerical schemes. Here, both the Newton-Raphson and arc-length methods are employed. The present derivation of the geometric stiffness matrix differs from the conventional procedure for nonlinear beam formulations [25]. In particular, the full BE beam metric is incorporated in Eq. (55) and the symmetric nature of the geometric stiffness term is stressed by the elegant and compact form of Eqs. (56) and (59). The latter confirms that the formulation is valid and that the adopted force and strain quantities are energetically conjugated.

Numerical examples
The aim of the subsequent numerical experiments is to validate the proposed approach and to examine the influence of curviness on structural response.
The boundary conditions are imposed strongly and the rotations are treated with special care, see [28], since they are not utilized as DOFs. Locking issues are not considered, but the presence of membrane locking is already detected for linear analyses in [28], and it is present in nonlinear analysis as well. However, its influence is alleviated with the usage of higher order basis functions [44]. It is interesting that high interelement continuity increases the locking effect due to the presence of more constraints, which must be satisfied. This issue can be dealt with an appropriate numerical integration scheme [44]. The quest for optimal quadrature rules in IGA is an ongoing topic, see e.g. [45], but the standard Gauss quadrature with p + 1 integration points are used here.
Since no rotational DOFs are required for BE beams, the implementation of multi-patch structures receives much attention [46,25,47,48]. A straightforward algorithm for handling multi-patch structures is implemented here [36]. The rotation at the end of NURBS curve depends on the displacements of the last two control points [37]. Therefore, the rotation of two NURBS curves at a joint is a function of six displacement components. In order to prescribe rigid connection between two patches, one displacement component is constrained as a function of the other five. This approach is straightforward and does not require the introduction of bending strips or end rotational DOFs [49,25].
All the results are related to the load proportionality factor (LPF), rather than the load intensity itself. Although the implemented code allows the usage of large load increments for some examples, the results are mostly calculated with small increments in order to enable a clear comparison of the obtained results via smooth equilibrium paths. For problems that do not exhibit a snap behavior, it is possible to obtain the deformed configurations for the prescribed load increments. Thus, the relative differences between different models can be compared along the equilibrium path, c.f. Subsection 5.1. These graphs give a clearer insight into discrepancies of the models than a simple comparison of equilibrium paths.
Since the time is a fictitious quantity in the present static analysis, strain and stress rates are equal to strains and stresses, respectively.

Cantilever beam subjected to an end concentrated moment
This example is a classic benchmark test for the validation of nonlinear beam formulations [25,23,24]. The cantilever beam is loaded at its free end with a moment 2nEIπ/L, where n defines the number of circles into which the beam rolls-up, Fig. 2  A moment load can be applied in various ways for rotation-free models. For example, an additional knot can be inserted to define a couple [23], or the linear stress distribution can be imposed, [20]. Here, no additional knots are added and the force couple is defined by the external virtual power, analogous to [28]. Since we are dealing with rotation-free nonlinear analysis, this couple changes direction and intensity throughout the whole deformation. Therefore, the external load vector must be updated at each iteration to define the applied moment accurately. Similar to [23], the contribution from these non-conservative external forces to the tangent stiffness matrix is disregarded.

Convergence analysis
We examine the convergence properties of the developed formulation for two constitutive models, i.e., D 1 and D a , using n = 1 and h = 0.2 → Kh ≈ 0.13. Since these two models return different structural responses, different reference solutions must be utilized. For the D 1 model, reference analytical solution for thin beam is employed [23]. Regarding the D a model, the analytical solution is not available in literature and a mesh of 60 quintic C 1 elements (484 DOFs) is adopted as a reference solution. The L 2 -norm of the relative error of displacement for LPF = 1 is considered in Fig. 3 for three different polynomial orders with highest available interelement continuity. Since the expected order of convergence is p + 1, [9], the obtained orders are higher. Fig. 4 illustrates the convergence of the normal force and the bending moment using the D a model. The analytical solution is used as reference. The expected order is p and it is obtained in all tests. The exception is the mesh with cubic splines. A similar behavior has been observed in [14] and attributed to the higher order derivatives involved in the strong form. In the rest of this example, quartic splines with C 3 interelement continuity are exclusively used.

Comparison of constitutive models
A discretization with 24 elements is employed to calculate the equilibrium path of the tip for n = 2 and two different values of curviness. The results are compared with the analytical ones and shown in Fig.  5. For a beam with the curviness Kh ≈ 0.06, the results of all models are in almost full agreement with the analytical solution, Fig. 5a. However, when the cross section dimensions are increased such that the curviness becomes Kh ≈ 0.50, the discrepancies become significant, Fig. 5b. The D 1 model is fully aligned with the analytical solution for thin beams, while the D 0 model slightly deviates. On the other hand, the results obtained with the D a and D 3 models are virtually indistinguishable, but differ from those of the thin beam. The D 2 model deviates from the exact predictions, but with a different sign than the D 0 and D 1 models. These results show that the present approach can return accurate results for problems with large displacements and rotations. For beams with small curvature, the results are aligned with the classic predictions. As the curviness increase, the difference between constitutive models becomes clearly visible. In order to present these findings more concisely, the relative differences of the reduced models with respect to the D a model are presented in Fig. 6 for four different values of curviness. These results confirm that the reduced model D 3 returns reasonably accurate results, even for strongly curved beams. It is interesting to note the relatively large differences of the D 2 model, which suggest that the inclusion of higher order terms, as in the reduced constitutive model D 3 , is required if accurate results for strongly curved beams are sought, see Eqs. (43) and (44). Another important conclusion follows from this specific example. Since the curviness is constant along the beam, the introduced exact metric has more impact on the beam response than it is the case with examples where the maximum curviness is local. We can observe, e.g., that if the curviness along the whole beam is less than 0.15, the error of the simple decoupled D 0 model is less than 0.3 %. As the curviness increases over 0.25, the error increase over 1.1 %, etc.

Specific aspects of strongly curved beams
This example is sometimes misinterpreted due to the pure bending conditions. It is assumed that there is no change of the length of beam axis and the beam rolls-up into a perfect circle with the circumference equal to the length of the beam [19]. However, when we deal with a geometrically exact analysis that considers the effect of curviness, the strain is distributed nonlinearly across the height of the cross section with a non-zero value at the centroid. This strain is extensional and its value is negligible for very thin beams, but becomes significant for strongly curved beams. Consequently, the axis of a cantilever loaded with tip moment will reach full circle for n < 1 and its circumference will be somewhat larger than the original length. In order to depict this, the four configurations of the beam with h = 0.2 are visualized in Fig. 7. Note how the initially thin and slender beam becomes strongly curved during the deformation. As a consequence, the metric, axial strain, and stress across the height of the cross section are distributed nonlinearly. The distribution of stress is shown next to the deformed configurations and compared with the classic linear distribution. Evidently, the intrados stress is greater, and the extrados stress is lesser than the equivalent linear stress. Due to the small extension of beam axis, the ends of the beam overlap for LPF = 0.5 and LPF = 1. This overlap for LPF = 0.5 is zoomed in Fig. 7b.
Furthermore, due to the axial strain of the beam axis, it is necessary to include the influence of flexural strain when calculating the normal section force [28]. To examine this effect, 40 elements are used for a model with n = 1 and Kh = 0.25, and the influence of both axial and flexural strains is observed. All five constitutive models are employed and the results are presented in Fig. 8. Clearly, the normal force oscillates, but with negligible amplitude. We see that the axial strain gives tension while flexural strain results in compression force. For the D a model, these contributions have nearly equal absolute values and the normal force oscillates around zero. The result is similar for the D 3 model, the normal force oscillates about the value close to 0.1. For the D 2 model, the error becomes pronounced since the normal force oscillates around approximately 9.3. Due to the fact that the influence of flexural strain on normal section force is disregarded in the D 1 model, the only contribution to the normal force comes from the axial strain, which oscillates around zero for this model. This result is considered accurate, but it stems from the fact that the two errors have canceled each other out. I.e. the dilatation is erroneously obtained as near zero, and the influence of the flexural strain is disregarded. These two approximations gave the correct value of normal force. Finally, the D 0 model gives a completely inaccurate result. For this model, the axial strain is negative and there is no contribution from the flexural strain.

Test of the algorithms for the update of internal forces
The simplicity of this example is ultimately utilized in testing the algorithms for the update of internal stress. Three approaches are considered: (i) the current stress is calculated from the total strain, (ii) the current stress is calculated by Eq. (33) with the first order approximation of adjustment term, and (iii) the current stress is calculated by Eq. (33) with the zeroth order approximation of adjustment term. These approaches are designated as F1, F2 and F3, respectively, and they are elaborated in Appendix A. The physical normal force, derived in [28], is: where χ is the change of curvature with respect to the Frenet-Serret frame of reference [28,50]. From the condition of pure bending (N = 0), we can calculate the exact physical axial strain of beam axis as: since χ = K = LPF (2nπ/L) for this example. This exact value of axial strain is compared with the results obtained with three different numerical approaches for the calculation of internal forces, Fig. 9. For the first load increment, which is 0.05 here, the internal forces are calculated in a same way for all algorithms. Interestingly, for the model with Kh ≈ 0.06, all three algorithms return the same result as Eq. (67) for LPF = 0.05. The F1 approach gives the most accurate results and this is also approach utilized in all present examples. The F2 algorithm can be successfully employed, however, for those beams with moderate curvature. The algorithm F3 is not recommended since it gives large errors, even for the small values of curviness. Inevitably, both F1 and F2 algorithms give erroneous results as the curviness and strain increase.

Lee's frame
This example is also well-established in the literature in the context of nonlinear analysis of in-plane beams [51,24,52]. The frame consists of two rigidly connected beams and the displacement components at the point of force application are observed, Fig. 10a. Each beam is discretized with five quartic C 1 elements. The results obtained with the presented models are given in Fig. 10b and they are in full agreement with those from literature [25]. This is expected because the maximum curviness in this example is local and less than 0.1. In order to examine the influence of curviness on this structure, the height of the cross section is multiplied by 2.  Figs. 11 and 12, while the curviness at three characteristic sections is given in Fig. 13. It is evident that the curviness at these sections follows similar path for both observed examples. The difference in the value of curviness practically equals the difference in cross section heights. The maximum curviness for LF1 is Kh ≈ 0.23, while for LF2 it is Kh ≈ 0.46. In the context of the structural response of these frames, all constitutive relations return similar results in the case of LF1. inspecting the zoomed equilibrium path for the component v A allows us to estimate the difference between D a and D 0 models which is close to 0.2% for LPF = 1, see Fig. 11b. For the LF2 example, the discrepancies are more noticeable. The largest Figure 11: Lee's frame. Equilibrium paths for the frame LF1 using five different constitutive models: a) complete equilibrium path, b) zoomed parts of the equilibrium path.  difference between D 0 and D a models is near 0.8% at the final configuration while the results returned by the D 2 , D 3 , and D a models are virtually indistinguishable. Although the maximum curviness for the examples with increased cross section height is large, the difference between displacements obtained by different constitutive models is significantly lower in comparison with the previous example, Fig. 6. This is mainly due to the fact that the maximum curviness of Lee's frame is local, while in the previous example it was constant along the whole beam and, therefore, more impactful on the structural response. This local character of the maximum curviness is underlined in Fig. 14 where the deformed configurations for LF1 and LF2 are displayed for three LPFs. Moreover, the rotations of adjacent sections at the joint of patches match, which confirms the correct application of the constraint equation.
Finally, the stress at the section B is considered in Fig. 15. The constitutive model D a is utilized for the calculation of displacements and reference strains. Then, two expressions for the distribution of stress are utilized: (i) the exact one, based on Eqs. (29) and (31), and (ii) the linear one, based on the same equations but with K = 0. These distributions are shown below the equilibrium paths for the stress at outer fiber in Fig. 15. It is clear that, the exact stress distribution should be utilized in a post-processing phase even for structures with a maximum local curviness less than 0.1, such as the original example LF. The difference of extreme stresses between the exact and the linear stress distributions is close to 13% for the LF example. As the curviness increases, this error becomes more pronounced.

Multi-snap behavior of a parabolic arch
The final example deals with the nonlinear response of a parabolic arch depicted in Fig. 16a. The arch is simply supported and loaded with the vertical force at the apex. For shallow arches, the snap-through phenomenon occurs after which the structure reaches a stable equilibrium and shows stiffening behavior. However, for deeper arches, a multiple snap behavior can be expected [53,54,55,56,33]. The aim of this example is to show that the present formulation is capable of describing highly-complex responses of arbitrarily curved structures. Furthermore, we discuss a physical reasoning behind the observed behavior [53]. We consider four parabolic arches labeled R1, R2, R3, and R4. They differ in rise f and the magnitude of the applied load P, as detailed in Fig. 16a. Due to the symmetry, one half of the arch is modeled with 16 cubic C 2 elements. As a guiding solution, we use an Abaqus simulation, which employs about 30 straight cubic B23 elements [57]. It is worth noting that Abaqus could not converge with denser meshes. Fig. 16b shows the signed curviness at the apex for the beam with the largest curvature, i.e., R4. The graph is interesting due to the multiple limit points. The important fact to note is that the maximum local curviness of this beam is lower than 0.1.
For the R1 case, a relatively simple structural response is obtained, and the equilibrium path for v s is shown in Fig. 17a. As the rise increases to 0.5 m, one snap-back occurs, Fig. 17b. Since these two arches have curviness less than 0.025, only the results obtained with the simplest, D 0 , and the most elaborated, D a , models are displayed. They agree well for almost the complete equilibrium path while negligible  discrepancies can be observed at the load limit points. Also, the results for the R1 case are practically indistinguishable from those obtained by Abaqus. Regarding the R2 case, small differences occur at the load limit points. As the rise increases further, the arch behavior becomes more complex and multiple snaps occur, see Fig. 18. Since the results for the different constitutive models are in good agreement, it can be concluded that the influence of the curviness is not significant for these arches. However, an inspection of the equilibrium paths near the load limit points reveals that a difference between the present approach and the B23 element discretization exists, which increases with the curviness and complexity of the equilibrium path. The results obtained by the D 3 and D 2 models are practically identical to those of the D a model. Furthermore, the models D 0 and D 1 deviate slightly from these results in the vicinity of limit points, especially for the R4 case. It is noteworthy that all models return the same final deformed configuration.
In order to analyze this multi-snap phenomenon more thoroughly, the eight different deformed configurations of the R4 case at LPF = 0 are given in Fig. 19. It can be seen that these configurations are similar to the buckling/vibration eigenmodes of a simple straight beam. In particular, the configurations designated with a, b, c, and d resemble the 3 rd , 5 th , 7 th , and 9 th beam eigenshape, respectively. Each of these configurations is stiffer than the previous one and the absolute value of load limit points increases as well. After configuration d, which is seemingly straight, the arch cannot generate any more half-waves and starts to release accumulated strain energy. Hence, the point d in Fig. 19 approximately represents an inflection point of the equilibrium path. After this the load limit points decrease and the arch passes through all the previous configurations in reverse order. Finally, a stable equilibrium branch is reached after the last load limit point.
It should be noted that this multi-snap behavior depends heavily on the imposed conditions of symmetry. In reality some small perturbations would probably exist in either the force position or the geometry, which would result in a simpler response [54].
The detected behavior is further investigated by observing the normal force at the apex. The applied forces for all cases are scaled in relation to the largest one. Fig. 20 illustrates the resulting equilibrium paths and marks the values of 3 rd , 5 th , 7 th , and 9 th critical buckling loads of a simply supported, axially loaded beam with length L. Additionally, an estimate of the maximum possible axial strain of the beam axis is taken into consideration by calculating the initial length of the parabolic arch L * and by finding the Green-Lagrange strain: This is the strain of the parabola which deforms into a straight line. Multiplication of this strain with the  axial stiffness EA 0 gives an estimate of a maximum normal compression force that can be generated in these arches due to the considered applied load. These forces are also designated in Fig. 20.
The number of half-wavelength forms that an arch can attain during this symmetric snapping is limited. This limit is directly influenced by the length of the arch and the critical buckling forces of a simple straight beam with length L. To summarize, the arch R4 cannot make the 11 th eigenshape, since its maximum axial strain cannot generate the normal force that is large enough to snap the beam into this mode. The analogous conclusion is valid for all observed cases. When there are no more new stable forms to reach, the curvature of equilibrium path changes sign and the arch goes through all these configurations in reverse, until it finds the stable form. A similar discussion takes place in [53], where the rise, thrust and flexural stiffness are identified as the causes for multi-snap behavior of circular arches. However, the graphical description given in Fig. 20 provides additional insight into the nature of the phenomenon.

Conclusions
A rigorous metric of the plane BE beam is utilized consistently for the derivation of the weak form of the equilibrium. The spatial discretization of the virtual power is performed by IGA. The introduction of the full beam metric provides a higher-order accurate BE beam formulation. Four simplified models are derived in addition to the exact constitutive relation. The comparison of these models via numerical examples allows a detailed analysis of the influence the formulation components have and of the most importance, the beam's curviness. The present formulation is geometrically exact, rotation-free, and fully capable of dealing with complex responses of multi-patch structures.
The results show that, in order to correctly determine the axial strain at the centroid of a curved beam, a rigorous computational model must be employed. For beams with curviness Kh < 0.1, the simple decoupled equations return reasonably accurate results for the displacement field. As the curviness increases, its influence becomes noticeable and as a result a more involved and complex model is required. An important factor is the domain where the strong curviness exists. If this is relatively local, its influence on the global response will not be as significant as when large portions of the structure are strongly curved. Regardless of the constitutive model, utilizing the exact expressions for the equidistant strains and stresses is recommended, since these have a nonlinear distribution and this even for beams with a small curvature. The inclusion of these expressions in a post-processing phase is a simple method for improving accuracy.
A phenomenon of multi-snap behavior of a parabolic arch is revised and elaborated. The arch accumulates strain energy and deforms into the configurations that resemble the eigenmodes of a straight beam. After the limit eigenshape is reached, the arch releases the accumulated strain energy and finds a stable equilibrium branch.
An interesting direction for future research is the application of the proposed formulation to the statics and dynamics of spatial beams.
This approach for the update of internal forces is the rational one, since the integration across the cross section should be avoided at each increment. An alternative method, that proved even more accurate in this research, is to calculate the internal forces from the total strain and the current constitutive matrix. However, this approach is not theoretically sound since the linear constitutive relation is valid only for the increments of strain and stress. Nevertheless, if the strains are not large, the results in Section 5 suggest that this approach can return acceptably accurate results. As noted previously, the latter approach is designated as F1 and former one as F2. Additionally, an incremental approach for the update of internal forces without any adjustment of previously calculated forces is implemented and marked as F3. When observing the Eq. (A3), it can in fact be considered as the zeroth order approximation. The results of these three approaches are compared in Subsection 5.1.4.