Skip to main content
  • Research Article
  • Open access
  • Published:

Introducing the Logarithmic finite element method: a geometrically exact planar Bernoulli beam element

Abstract

We propose a novel finite element formulation that significantly reduces the number of degrees of freedom necessary to obtain reasonably accurate approximations of the low-frequency component of the deformation in boundary-value problems. In contrast to the standard Ritz–Galerkin approach, the shape functions are defined on a Lie algebra—the logarithmic space—of the deformation function. We construct a deformation function based on an interpolation of transformations at the nodes of the finite element. In the case of the geometrically exact planar Bernoulli beam element presented in this work, these transformation functions at the nodes are given as rotations. However, due to an intrinsic coupling between rotational and translational components of the deformation function, the formulation provides for a good approximation of the deflection of the beam, as well as of the resultant forces and moments. As both the translational and the rotational components of the deformation function are defined on the logarithmic space, we propose to refer to the novel approach as the “Logarithmic finite element method”, or “LogFE” method.

Background

We propose a novel finite element formulation, the Logarithmic finite element, or “LogFE” method, that significantly reduces the number of degrees of freedom necessary to obtain accurate approximations of boundary-value problems. The LogFE method focuses on the low-frequency part of a deformation and minimizes spurious high-frequency components in the solution.

In order to keep the exposition as simple as possible, we restrict the model presented in this paper to the case of a planar Bernoulli beam, i.e. a beam endowed with Bernoulli kinematics embedded in the Euclidean plane. In addition, we limit the degrees of freedom to coefficients related to rotations and dilatations at the nodes of the element. While we restrict the numerical examples to the evaluation of a beam consisting of one single element only, we explicitly show that degrees of freedom related to adjacent finite elements can be linked together by linear maps, based on geometrically meaningful continuity conditions. Hence, the construction of a global finite element system based on the beam elements presented in this work is possible.

Following the approach outlined in [24] and [25], we intend to present a formulation that includes, in addition to rotations and dilatations, the translations at the nodes of the configuration in a future publication.

Interpolating on the logarithmic space

Geometrically exact beam formulations generally include both translations and rotations as kinematic variables, following the continuum mechanics model of the elastica developed by Cosserat and Cosserat [9]. In most of the literature, the interpolation of the kinematic variables at the integration points is performed separately for the translational and the rotational variables. Well-known deficiencies of many of these formulations, such as lack of geometric invariance, path dependence and poor accuracy in coarse meshes, have been traced back to the characteristics of the various methods proposed for the interpolation of the rotational variables by Jelenić and Crisfield [14] and Romero [22], among others.

The modelling of beams endowed with Bernoulli kinematics within this setting presents particular challenges, especially membrane locking, as rotations and translations, although interpolated separately, must jointly satisfy the Bernoulli condition along the neutral axis. Similar difficulties characterize finite element formulations based on Lagrange functions for the axial and Hermitean functions for the transversal displacements [3, 20]. As a result, the characteristics of existing finite element formulations for geometrically exact Bernoulli beams remain inferior to those available for kinematics that include shear deformations [17]. Armero and Valverde provide an account of the historical development of geometrically exact beam formulations, including a brief description of different methods that aim to eliminate or reduce the effects of the drawbacks associated with the classical beam formulations [3].

To preserve the orthogonality of the interpolated directors, a number of existing geometrically exact models, such as [11, 14, 27], identify the kinematic variables with elements of a Lie group, the special orthogonal group SO(3). (In the work of Ericksen and Truesdell, the use of the Lie group and its associated Lie algebra is implied by the differential equation given on page 306.)

In a different approach, Armero and Valverde interpolate the director frame associated with the neutral axis of the underlying geometry of a Bernoulli beam by using its angle in the plane case [3], and in the general linear group \({GL \,}{(3, \mathbb {R}{})}\) for the three-dimensional case [4]. In the latter case, orthogonality is achieved by applying a projection operator to the resulting interpolated set of vectors that preserves the direction of the vector tangent to the neutral axis, effectively resulting in an interpolation on the special orthogonal group. As the shape functions are defined not on the global coordinate system, but with regard to the local coordinate systems induced by the director frame, membrane locking is minimized due to an appropriate coupling of axial and transversal displacements. Meier, Popp and Wall extend this approach from the linear domain to large deformations, introducing an orthogonal interpolation method based on the torsion of the neutral axis [17].

Betsch and Steinmann enforce the orthogonality of the director frame at the nodes by introducing Lagrange multipliers, thus restricting the function space of the weak form of the equations of motion [6]. The effects of a non-orthogonal interpolation of the director frame can be addressed by reformulating the weak form of the equations of motions in skew coordinates, thereby obtaining a geometrically invariant formulation [12].

Although the concept of the Lie group and the Lie algebra is central to the method we are proposing in the present work, we do not interpolate the kinematic variables themselves in a Lie group or in its associated Lie algebra. Instead, we aim to identify, in a function space generated by polynomial shape functions on the logarithmic space, i.e. the space of the Lie algebra, a deformation function that, in the case of the Bernoulli beam, transforms the neutral axis of the given initial configuration so as to obtain an equilibrium configuration.

Choosing appropriate shape functions, both with regard to their scalar part and to their respective vector in the Lie algebra of the deformation function, is of crucial importance for the performance of the proposed method. In particular, although we do not allow translations at the nodes in the present work, the vectors associated with the shape functions are not restricted to either the rotational/dilatational or the translational part of the Lie algebra associated with group of planar similarity transformations, \({sim}{(2, \mathbb {R}{})} = \left( {\mathfrak {gl}}\,(1, \mathbb {R}{}) \oplus {\mathfrak {so}}(2)\right) \oplus _{\mathrm{Id}} {\mathfrak t}({2, \mathbb {R}{}})\). (In this formula, \(\mathrm{Id}\) identifies the function \(\varphi :\mathfrak {g}\rightarrow End(\mathfrak {n})\) defining the commutator of the semidirect sum of two Lie algebras, \(\mathfrak {g}\) and \(\mathfrak {n}\), \(g_i\in \mathfrak {g}, n_i\in \mathfrak {n},\left[ (g_1,n_1),(g_2,n_2)\right] :=(g_1,n_1)(g_2,n_2)-(g_2,n_2)(g_1,n_1)\), with \((g,n)(g{'},n{'}):=(gg{'},\varphi _{g}(n{'}))\) [14, p. 84–85] [17, p. 38–40]. For matrix Lie algebras, we assume that \(g(n) = gn\) is given by the matrix product). Instead, the vectors span across the entire Lie algebra \({sim}({2, \mathbb {R}{}})\), inducing a strong coupling of rotational and translational kinematic variables, which results from the multiplicative group operation given as a semidirect product.

Accessible introductions to Lie group theory can be found in [13] and [17]. [2] develops formulations for the description of large rotations from elementary calculations, leading to a geometrically motivated approach to the theory of Lie groups and Lie algebras. A number of more advanced topics in Lie group theory, including the construction of semidirect products, are discussed in [26] and [16].

Interplay of rotations and translations

Due to the specific problems related to the interpolation on the group of rotations, which are not present with regard to translations, many geometrically exact formulations for finite beam elements rely on a complete separation of the rotations from the translations, e.g. [14, 23, 27]. This uncoupling has the obvious advantage of isolating the shape functions related to the rotations from spurious effects of the approximation of the translational part of the deformation. However, it also results in a larger number of degrees of freedom, opening up dimensions of the deformation space which are largely irrelevant in order to achieve a reasonably close approximation of the exact deformation, especially when focusing on its low-frequency component.

In contrast to the approach referred to above, which effectively is based on linear interpolations on the Lie algebra \({\mathfrak {so}}\left( 3\right) \oplus {\mathfrak t} \left( {3, \mathbb {R}{}}\right) \), the direct sum of the Lie algebras associated with the Lie groups of rotations \({ SO }(3)\) and of translations \(T ({3, \mathbb {R}{}}) \simeq \mathbb {R}^3\) in three dimensions, Sonneville, Cardona and Brüls offer a geometrically exact formulation based on the Lie algebra associated with the semidirect product of translations and rotations, i.e. the Lie algebra \({\mathfrak {se}}(3)\) associated with the special Euclidean group \({ SE }(3) := { SO }(3) \ltimes _{\mathrm{Id}} T ({3, \mathbb {R}{}})\), the semidirect product of rotations and translations in three dimensions [28]. By identifying the local material frames, i.e. the positions and orientations of the cross-sections of the beam along the neutral axis, with elements of the Lie group \({ SE }(3)\), and linearly interpolating on the associated Lie algebra, the formulation allows for a coupling of the translational and rotational components of the deformation. An earlier formulation that implicitly uses an interpolation of the neutral axis of a beam as well as of the orientation of its cross-section on a the Lie algebra \({\mathfrak {se}}(3)\), instead of polynomial shape functions, has been presented by Borri and Bottasso [7]. Selig and Ding explicitly introduce Lie groups and Lie algebras in their exposition of a screw-theoretic formulation for a planar beam [26]. The formulations proposed in these works are ultimately based on concepts from screw theory [5, 15], also referred to as motor calculus [30, 31].

To the knowledge of the authors, existing formulations for geometrically exact beam models, whether based on the Lie algebra \({\mathfrak {so}}(n) \oplus {\mathfrak t}({n, \mathbb {R}{}})\), or on the Lie algebra \({\mathfrak {se}}(n)\), for \(n = 2\) or \(n = 3\), rely on a strictly linear interpolation between the degrees of freedom given at the respective nodes (in the logarithmic space) and thus do not make use of internal degrees of freedom associated with additional shape functions, such as polynomial shape functions of higher order, that may be constructed on the Lie algebra. As a result, the image of the neutral axis of a beam element is located on an arc or a section of a helix with curvature and torsion depending exclusively on the nodal degrees of freedom. In particular, the curvature and the torsion remain constant along the neutral axis on each beam element for both the formulation based on the special orthogonal group [22, p. 125] as well as for the formulation based on the special Euclidean group [28, p. 460]. Due to the relatively low computational complexity of a single element, this restriction can often be compensated for by a finer discretization of the beam model. For example, Sonneville, Cardona, and Brüls [28] achieve good convergence characteristics and small approximation errors for a standard test case based on a bent Cantilever beam subjected to a fixed load.

However, in order to reduce the number of nodal degrees of freedom, as well as the overall number of internal and nodal degrees of freedom in the approximation algorithm, it is necessary to overcome the restriction given by a constant curvature along the neutral axis of a single beam element. In addition, by reducing the degrees of freedom, it is possible to focus on the degrees of freedom associated with the low-frequency component of the deformation while minimizing the spurious high-frequency effects. In a multigrid setting, this has the advantage of being able to separate the low-frequency component of the deformation, which is being approximated on the coarse grid, and the high-frequency part which may be approximated by a high-frequency smoothing algorithm on the fine grid.

In the various formulations referred to above, torsion and shear deformations result from the interpolation of translations related to nodal translational degrees of freedom. As the formulation of the planar Bernoulli beam element presented in this work does not include translational degrees of freedom, there is no specific mechanism that would provide for the optimization of shear deformations and shear stresses. The introduction of explicit shear deformations by including shear-related degrees of freedom at the nodes would greatly complicate the functions linking the degrees of freedom to the respective functionals of the interpolant on the border of the finite element, as shear deformations and rotations do not commute. In order to keep the subsequent exposition, which intends to highlight the basic properties of the Logarithmic finite element method, as simple as possible, we do not include degrees of freedom related to translations or shear deformations at this point. The deformation space resulting from restricting the formulation to rotational and dilatational degrees of freedom, in turn, is clearly better suited to the modelling of Bernoulli kinematics rather than Timoshenko kinematics. In particular, if the orientation of the cross-section would be taken directly from the deformation function, rather than implicitly calculated from the orientation of the neutral axis, the formulation would imply the absence of shear deformation at the nodes, also for Timoshenko kinematics, a clearly non-physical result.

The deformation function

In contrast to [11, 14, 27], the formulation presented in this paper is based on finding a deformation function acting on the given initial configuration that realizes an equilibrium configuration of the system. It is thus the deformation function that results from a function that maps elements of the parametrization of the beam to elements of the Lie group. Thus, both the initial configuration and the current configuration are fully defined by the positions of the points along the neutral axis of the beam, with rotations obtained from the derivative of the neutral axis. As a result, in order to identify individual points on the interpolant with elements of a Lie group, one would necessarily need to take the location of the neutral axis in a neighborhood of the respective points into account.

While the interpolant of the resulting finite element contains both a rotational and a translational component, the degrees of freedom at the nodes—in the simplified variant that we describe in this work—only include rotations and dilatations (i.e. radial displacements with respect to a given center). Dilatational and rotational degrees of freedom determine the local dilatation and rotation in a neighborhood of the associated node, respectively, while the position of the node itself remains unchanged. The translational component on the interior of the interpolant thus arises from the impact of rotations and dilatations associated with different nodes belonging to the same finite element, as well as their interaction.

In general, the deformation functions obtained with the formulation proposed in the present work will result in current configurations that are not characterized by a constant curvature within the single finite elements, even when restricting the shape functions on the Lie algebra to linear interpolations.

Essential characteristics

The essential characteristics of the LogFE method, which distinguish it from the different approaches referenced above, can thus be summarized as follows: In the proposed formulation, the degrees of freedom, together with the shape functions, determine a multiplicative deformation function that is defined on the Lie group of planar similarity transformations \({Sim}({2, \mathbb {R}{}})\). Associating the scalar part of the shape functions with vectors that span across the rotational and translational subalgebras of its associated Lie algebra induces a strong coupling of translational and rotational variables of the configuration. In addition, it allows for the introduction of polynomial shape functions of higher order on the logarithmic space.

Given an appropriate choice of shape functions, an additively separable linear correspondence between certain field values (obtained as functionals of the interpolant on the border of a finite element) and the function parameters of the interpolant on a given finite element can be achieved. This allows for the formulation of continuity conditions between finite elements, so that a global finite element system can be set up, following the standard procedure. LogFE elements endowed with Lagrange-type or augmented Lagrange-type shape functions can be linked by rotational values, both with each other as well as with other finite elements that accommodate continuity conditions based on rotations. Elements endowed with Hermite-type shape functions may also be linked to other finite elements through common values of curvature and of the derivative of the strain at a common node. Similarly, boundary conditions may be incorporated into the set-up of a global finite element system.

Multigrid methods

Multigrid methods rely on the interplay of a smoothing algorithm on the fine grid and a general solver on the coarse grid. To obtain good convergence characteristics, the coarse grid algorithm must operate on the low-frequency part of the approximation error. At the same time, its influence on the high frequencies must be minimized. Employing a coarse grid correction that supplies, together with an approximation of the degrees of freedom, an interpolant characterized by a small high-frequency component minimized two potential difficulties that may arise in multigrid-based calculations: it can prevent the emergence of situations in which effects related to the high-frequency component of the deformation obstruct the approximation on the coarse grid, and it minimizes the risk that changes in the approximation of the high-frequency component induced by the coarse grid correction degrade the convergence characteristics of the algorithm applied on the fine grid. For a discussion of the impact of the order of the prolongation function being used to transfer the coarse grid correction to the fine grid, see [19], as well as the literature referred to in that publication.

We envisage multigrid methods as a major application of the LogFE formulation. However, we have chosen not to present a complete multigrid formulation in the present paper. The model satisfies the conditions necessary for its incorporation in a multigrid formulation as a coarse grid solver, as it focuses on the low-frequency component of the deformation and allows for a straightforward calculation of the values of the degrees of freedom on the fine grid from the interpolant on the coarse grid. We therefore leave the practical implementation of the LogFE model in a multigrid context to future research.

Outline of the paper

The outline of the remaining sections is as follows: “Methods” section describes the application of the LogFE method for planar beams endowed with Bernoulli kinematics. “Kinematics” section describes the kinematic assumptions of the model. It provides some mathematical background on the theory of Lie groups and Lie algebras. “Continuity conditions” section presents the field values—given as functionals of the interpolant—that may be used as continuity conditions, and establishes a linear correspondence between these values and the parameters of the interpolant. In “Consistency with the linear beam theory” section, the model, in the limit of infinitesimal deformations, is related to the linear theory of the Bernoulli beam. The following part of the paper, “Quasi-static equilibrium” section, presents the equilibrium condition and formulates the evaluation of the non-linear kinematics at the element level. “Results and discussion” section examines the results obtained by the LogFE method for different boundary conditions, load characteristics and load intensities. “Conclusion” section summarizes the results and relates them to the basic characteristics of the model.

Some aspects more closely related to the implementation of the method have been included in the appendices. Appendix 1 outlines how LogFE element formulations adjusted for given boundary conditions can be constructed. Appendix 2 presents a sharp upper bound for the absolute value of the error of an approximate solution for the derivative of the matrix exponential.

Methods

Kinematics

The basic concepts of the LogFE model can be used to develop formulations for different classes of finite elements. However, in order to keep the exposition as simple as possible, and in order to focus on the essential aspects of the model, we restrict the following description of the LogFE model to the example of large deformations of a prismatic beam in the Euclidean plane with rectangular cross-section, endowed with Bernoulli kinematics. In addition, we limit the degrees of freedom to coefficients related to rotations and dilatations at the nodes of the element. The beam is assumed to consist of a homogeneous, isotropic material. The beam model can be classified as a geometrically exact beam formulation for finite rotations.

Fig. 1
figure 1

Configurations of a beam in the Euclidean plane

The deformation function

The continuous body \(\fancyscript{B}\), which is composed of the particles X of the beam, is given as a domain in a topological space (see Fig. 1). After placing the particles of the beam in the reference configuration, we discretize the beam into finite elements and endow the elements with a parametrization, such that set of particles of the beam along the codimension of the parameterization (i.e. the width) can be unambiguously identified with a tuple of parameterizing variables \(\left( \mathchoice{\xi ,\eta }{{\xi ,\eta }}{{\xi ,\eta }}{{\xi ,\eta }}\right) \in \mathbb {R}^2\). The variable \(\xi \) determines the location of the cross-section containing a given material point of the beam along the neutral axis. As a result of assuming Bernoulli kinematics, the position and orientation of the points in a given cross-section of the beam are fully determined by the location and the orientation of the neutral axis.

The location of a material point in the reference configuration is given by \(\mathbf X = \kappa _{\mathrm {ref}}({X})\), with \(\kappa _{\mathrm {ref}}\) mapping the domain \(\fancyscript{B}\) continuously into the Euclidean plane, which we will, after choosing an arbitrary origin, identify with the complex plane. This identification serves two purposes: it simplifies the notation, and it immediately clarifies, by the use of complex numbers instead of matrices representing linear maps on \(\mathbb {R}^2\), that, within this work, all linear maps operating on \(\mathbb {R}^2\) are restricted to the commutative subgroup \({GL \,}{(1,\mathbb {R}{})} \times SO{(2)} \simeq {GL \,}{(1,\mathbb {C}{})}\) of the general linear group \({GL \,}{(2,\mathbb {R}{})}\) given by the dilatations and rotations. In particular, the deformation function \(g_\xi \) (see Fig. 1), as a function of the degrees of freedom, assumes values that can be represented by matrices with complex entries.

We discretize the beam as a single finite element. The parameterization, a Lipschitz continuous embedding, maps into the reference configuration according to the function

$$\begin{aligned} \chi _{\mathrm {ref}}:\left[ 0,l\right] \times \left[ {-\tfrac{1}{2} h,\tfrac{1}{2} h}\right] \rightarrow \mathbb {C}{},\quad \left( \mathchoice{\xi ,\eta }{{\xi ,\eta }}{{\xi ,\eta }}{{\xi ,\eta }}\right) \mapsto \xi + \mathrm{i}\eta , \end{aligned}$$
(1)

in which l and h denote the length and the height of the beam, respectively.

In order to express both rotations and translations as a single, multiplicative operation, we will use the concept of homogeneous coordinates. By placing a vector space V into a larger construct \(V \times \{1\}\), we can express the translation of a vector in V as a multiplication of an element of \(V \times \{1\}\) by a matrix. Given \(\mathbf v_0, \mathbf v \in V\), \(\mathbf R \in {GL \,}(V)\), we have

$$\begin{aligned} \mathbf R \mathbf v_0 + \mathbf v \simeq \begin{pmatrix} \mathbf R \mathbf v_0 \\ 1 \end{pmatrix} + \begin{pmatrix} \mathbf v \\ 1 \end{pmatrix} = \begin{pmatrix} \mathbf R &{} \quad \mathbf v_0 \\ \mathbf 0 &{}\quad 1 \end{pmatrix} \begin{pmatrix} \mathbf v \\ 1 \end{pmatrix}. \end{aligned}$$
(2)

Expressing both linear maps and translations in this way also allows us to describe the logarithm of the function \(\mathbf v_0 \mapsto \mathbf R \mathbf v_0 + \mathbf v\), as the logarithm of the matrix \( {\bigg (\begin{array}{ll} \mathbf R &{} \quad \mathbf v \\ \mathbf 0 &{} \quad 1 \end{array}\bigg )}\). The logarithm plays an essential role in the theory of Lie groups and Lie algebras, on which the LogFE method is based. It is given as the (generally multi-valued) inverse of the exponential function, which is defined as \(\exp (\bullet ) := \sum _{k=0}^\infty \frac{1}{k!} \bullet ^k\). In this definition, the symbol “\(\bullet \)” may denote a scalar, a matrix, or a function. In the case of a function, taking the kth power is defined as consecutively applying the function for a total of k times.

We therefore embed both the initial configuration, \(\mathbf x _0\), and the current configuration, \(\mathbf {x}\), into the homogenized Euclidean plane \(\mathbb E^2 \times \{1\}\), which we identify with \(\mathbb {C}{} \times \{1\}\). Figure 2 illustrates the map from the parametrization space into the space of the initial configuration.

Fig. 2
figure 2

Parameterization of a two-dimensional beam element. The coordinate axes are related to the canonical coordinate system in the Euclidean plane by the relations \(x = \mathrm{Re}\) (x \(_0\)) and \(z = -\mathrm{Im}\) (x \(_0\))

The motion of the body from the initial configuration to the current configuration is given by a continuous map \(g(\xi )\) which depends on the parameter \(\xi \) and acts on the initial configuration, i.e. on \({\bigg (\begin{array}{ll} \mathbf x _0 \\ 1 \end{array}\bigg ) = \varphi _0(\mathbf {X}) = \varphi _0 \left( {\chi _{\mathrm {ref}}\left( {\xi ,\eta }\right) }\right) }\). A point \(\mathbf x\) in the current configuration is given by

$$\begin{aligned} \left( {\begin{array}{l} \mathbf{x } \\ 1 \\ \end{array}} \right) = g_\xi \left( {\begin{array}{l} {{\mathbf{x }_0}} \\ 1 \\ \end{array} } \right) =g_\xi ({\varphi _0}(\mathbf X )) = g_\xi \left( {\varphi _0}\left( {\chi _{{\text {ref}}}}(\xi ,\eta )\right) \right) . \end{aligned}$$
(3)

For points of the neutral axis, we have \(\eta \equiv 0\), thus their position in the current configuration only depends on the parameter \(\xi \). The map \(g(\xi )\) depends on the parameterization variable \(\xi \), but acts on the space of the initial configuration, i.e. on \(\mathbb {C}\times \{1\}\).

Lie groups and Lie algebras

All functions \(g(\xi )\) are given in the form \(\exp {({{\bar{\mathbf {Z}}}\left( \xi \right) })}\) and act on elements of \(\mathbb {C}\times \{1\}\). \({\bar{\mathbf {Z}}}(\xi )\) is an element of a Lie algebra and is therefore endowed with specific properties. We recall some basic results from the theory of Lie groups and Lie algebras, which are relevant for the subsequent formulation of the model.

Members of a set form a group, denoted G, if there exists an associative binary operation \(G \times G \rightarrow G\), there exists a neutral element \(e \in G\) with regard to this operation, and for every element \(g \in G\) there is an inverse element \(g'\) such that \(g' g = e\).

The values of the function \(g(\xi )\) are elements of the subgroup

$$\begin{aligned} \left\{ \begin{pmatrix} a &{}\quad b \\ 0 &{} \quad 1\end{pmatrix} \left| \phantom {\begin{pmatrix} a &{}\quad b \\ 0 &{} \quad 1\end{pmatrix}} a \in \mathbb {C}{} \setminus \{0\}, b \in \mathbb {C}{}\right\} \right. \end{aligned}$$
(4)

of the 2 \(\times \) 2 matrices over the complex numbers. They form a group with regard to the multiplication, which we will also denote G. The identity map, denoted \(\mathrm{Id}\), is the neutral element of G. While this group can be embedded in the general linear group \({GL \,}{(2, \mathbb {C}{})}\), as shown in Eq. (4), a more concise description characterizes this Lie group of the values \(g(\xi )\) as the group of complex similarity transformations \({Sim}({1, \mathbb {C}{}}) = {GL \,}{(1, \mathbb {C}{})} \ltimes _\mathrm{Id} T ({1, \mathbb {C}{}})\), which is isomorphic to the group of planar similarity transformations \({Sim}({2, \mathbb {R}{}}) = \left( {{GL \,}{(1, \mathbb {R}{})} \times SO(2)}\right) \ltimes _\mathrm{Id} T ({2, \mathbb {R}{}})\), with the group action defined as \(\left( \mathchoice{s,R,t}{{s,R,t}}{{s,R,t}}{{s,R,t}}\right) \circ \left( \mathchoice{s',R',t'}{{s',R',t'}}{{s',R',t'}}{{s',R',t'}}\right) = \left( \mathchoice{s s',R R',s R t' + t}{{s s',R R',s R t' + t}}{{s s',R R',s R t' + t}}{{s s',R R',s R t' + t}}\right) \), \(SO(2)\) being the special planar orthogonal group. With regard to these operations, s is a scaling parameter, R is a rotation matrix, and t is a translation vector. The action of this Lie group on a vector in the Euclidean space is defined by \(\left( \mathchoice{s,R,t\phantom {t'}}{{s,R,t\phantom {t'}}}{{s,R,t\phantom {t'}}}{{s,R,t\phantom {t'}}}\right) \circ x = s R x + t\), with \(x \in \mathbb {R}^2\). For the isomorphic complex Lie group, the group operation is given as \(\left( \mathchoice{a,b}{{a,b}}{{a,b}}{{a,b}}\right) \circ \left( \mathchoice{a',b'}{{a',b'}}{{a',b'}}{{a',b'}}\right) = \left( \mathchoice{a a',a b' + b}{{a a',a b' + b}}{{a a',a b' + b}}{{a a',a b' + b}}\right) \), and the action on a vector space as \(\left( \mathchoice{a,b}{{a,b}}{{a,b}}{{a,b}}\right) \circ z = a z + b\), with \(z \in \mathbb {C}{}\). As these group operations indicate, calculations on a Lie group generally do not involve multiplications of complete matrices, although many Lie groups, including those referred to in this work, can be embedded in the general linear group, thus facilitating the understanding of the group action and the action of the group on a vector space.

Each element of the group can be obtained by taking the exponential of an element \({\bar{\mathbf {Z}}}(\xi )\). Furthermore, \(\exp \left( {{t {\bar{\mathbf {Z}}}(\xi )}}\right) \in G\) for all \(t \in \mathbb {R}{}\). Therefore, G is a Lie group, and the set of elements \({\bar{\mathbf {Z}}}(\xi )\) forms the Lie algebra of G, denoted \(\mathfrak g\).

As a Lie algebra, \(\mathfrak g\) is a vector space together with the adjoint map, a skew-symmetric function \(\mathrm{ad}:\mathfrak g \times \mathfrak g \rightarrow \mathfrak g\), given as \(\left( \mathchoice{X, Y}{{X, Y}}{{X, Y}}{{X, Y}}\right) \mapsto \left[ X, Y\right] \), \(X, Y \in \mathfrak g\), with \(\left[ X, Y\right] := XY - YX\) for matrix Lie algebras, i.e. for Lie algebras that are subalgebras of a Lie algebra \({\mathfrak {gl}}\,\left( {n, \mathbb {R}{}}\right) \) associated with a general linear group \({GL \,}\left( {n, \mathbb {R}{}}\right) \), \(n \in {\mathbb {N}}{}\). The term \(\left[ X, Y\right] \) is also called the commutator of the Lie algebra. The specific Lie algebra that will be used in the remainder of the text, \(\mathfrak g := {sim}({1, \mathbb {C}{}}) = \mathfrak {gl} \left( 1,\mathbb C\right) \ltimes _\mathrm{Id} {\mathfrak t} ({1, \mathbb {C}{}})\), can be represented by the matrices

$$\begin{aligned} \left\{ \begin{pmatrix} a &{}\quad b \\ 0 &{}\quad 0\end{pmatrix} \left| \phantom {\begin{pmatrix} a &{}\quad b \\ 0 &{}\quad 1\end{pmatrix}}\right. a \in \mathbb {C}{}, b \in \mathbb {C}{}\right\} . \end{aligned}$$
(5)

In particular, it contains the elements \({\mathbf {Z}}_{I,r}\) given in Table 1. As will be described below in more detail, each element is associated with one or more shape functions. The index I designates the node associated with the shape functions constructed from the element \({\mathbf {Z}}_{I,r}\), while \(r\), which can assume the values 1 or 2, denotes the type of the deformation. The case \(r=1\) indicates that the shape function characterizes a dilatation with the position of the node I as its fixed point, while \(r=2\) indicates a rotation around the position of the node \(I\). By the term “dilatation”, we refer to a radial displacement of the material points with regard to a given center, the fixed point of the deformation. A dilatation thus induces volumetric strain in the body.

Table 1 Elements of the Lie algebra \(\mathfrak g\) and the associated Lie group G

Each of the elements \({\mathbf {Z}}_{I,r}\) generate a subalgebra of \(\mathfrak g\). We also note that \({\mathbf {Z}}_{I,1}\) and \({\mathbf {Z}}_{I,2}\) are independent vectors of \(\mathfrak g\), understood as a vector space. For \(\lambda _{I,1}, \lambda _{I,2} \in \mathbb {R}{}\), the element \(\exp {({{\bar{\mathbf {Z}}}_I})} \in G\), with \({\bar{\mathbf {Z}}}_I= \lambda _{I,1}{\mathbf {Z}}_{I,1} + \lambda _{I,2}{\mathbf {Z}}_{I,2}\), represents the action of simultaneously rotating and dilatating the initial configuration, with the position \(\mathbf x _0^I\) of the node \(I\) as the invariant point. The exponential of a linear combination of the elements \({\mathbf {Z}}_{I,r}\) represents a mixture of their respective actions on the initial configuration whose invariant point generally does not coincide with the position of either node.

A parameterized curve on the Lie algebra, seen as a vector space, can be given by

$$\begin{aligned} \gamma (\xi ) := \sum \limits _{\begin{array}{c} 1 \le I\le {{nel}}\\ 1 \le r\le 2 \end{array}} {\lambda _{I,r}} (\xi ) {{\mathbf {Z}}_{I,r}}, \end{aligned}$$
(6)

with \(\gamma (\xi ) \in \mathfrak g\). In this equation, \({{nel}}\) denotes the number of nodes per element. For each location \(\xi \), the curve \(\gamma (\xi )\) determines the map \(\exp \big (\gamma {(\xi )}\big ) \in G\), that will be applied to the material points of the initial configuration associated with \(\xi \). In particular, as every node \(I\) is associated with a specific value \(\xi _I\) of the parameterization variable, \(\gamma (\xi )\) provides an interpolation of the transformation applied to the respective nodes.

Shape functions

In the following, we will construct the shape functions for a finite beam element consisting of two nodes, i.e. \({{nel}}= 2\). Both nodes are located on the neutral axis. The coordinates of the locations of the node with index \(I\) are given as \(\left( \mathchoice{\xi _I,0}{{\xi _I,0}}{{\xi _I,0}}{{\xi _I,0}}\right) \). We define the function \(g(\xi )\) in Eq. (3) as

$$\begin{aligned} g_\xi := \exp \left( {{\bar{\mathbf {Z}}}\left( \xi \right) }\right) = \exp \left( {\sum _{\begin{array}{c} 1 \le I\le {{nel}}\\ 1 \le r\le 2\\ 0 \le q\le n_q-1 \end{array}} u_{I,r,q} \underbrace{\overbrace{ {p_{I,r,q}} {\left( {\psi _I} (\xi )\right) }}^{{N_{I,r,q}} (\xi )} {\mathbf {Z}}_{I,r}}_{{{\mathbf {N}}_{I,r,q}} (\xi )}}\right) , \end{aligned}$$
(7)

In this equation, \(r\) denotes the index of the respective element of the Lie algebra, and \(n_q\) denotes the number of polynomial functions \(p_{I,r,q}\) that are being used to construct the shape functions \({{\mathbf {N}}_{I,r,q}}(\xi )\) related to an element \({\mathbf {Z}}_{I,r}\) of the the Lie algebra \(\mathfrak g\), such that, for a given node \(I\), \({\mathbf {Z}}_{I,r}\) are linearly independent. \({N_{I,r,q}}(\xi )\) may be understood as the scalar part of the shape function, determining the “intensity” of the deformation, while its basic characteristics (dilatation/rotation, invariant point) are given by \({\mathbf {Z}}_{I,r}\). The term \({\mathbf {Z}}_{I,r}\) in Eq. (7) is generally not present in conventional finite element models, in which shape functions are understood as translations defined by the scalar value \(N(\xi )\) and a basis vector \(\mathbf e_i\), and are defined separately for each dimension in the vector space of the configuration. In such formulations, a coordinate-based description if often used, thus dropping the basis vector \(\mathbf e_i\). In the context of the beam model described, we have \({{nel}}= 2\), \(r \in \{1,2\}\), and \(n_q\ge 1\). \({\psi _I} (\xi )\) are the barycentric coordinates constructed on the interval \(\Omega = \left[ \xi _1,\xi _2\right] \) and associated with the node \(I\). The function \({\bar{\mathbf {Z}}}(\xi )\) thus is Lipschitz continuous, differentiable and bounded on \(\Omega \).

Continuity conditions

All numerical calculations in this work are restricted to dilatations and rotations at the nodes. It is, however, natural to ask whether the model can, in principle, be extended to a formulation including translations, and whether it is possible, based on that formulation, to perform numerical simulations not just for a single finite element, but for a complete finite element system, by applying standard procedures of the finite element on the global level. In order to be able to answer these questions in the affirmative, we show

  1. (1)

    that the formulation presented here is indeed a special case of an extended formulation that includes translations, and

  2. (2)

    that the extended formulation, given suitable shape functions, results, in a sufficiently large domain of the space of the degrees of freedom, in a completely additively separable linear isomorphism (i.e., a one-to-one correspondence with possibly different proportionality factors) between the degrees of freedom and certain functionals of the interpolant which can be used as continuity conditions. This condition corresponds to the third criterion for finite elements as given by Ciarlet in [8, p. 78–9].

Preliminary considerations

In addition to the Lie group G, we are considering the Lie group \(\tilde{G}\), and show that the application of \(\tilde{G}\), together with a suitable set of embeddings and projections, results in the same deformation function as the application of the Lie algebra G with the embeddings and projections used above. We identify the Euclidean plane with the complex numbers and consider the embeddings

$$\begin{aligned} i_1&:\mathbb {C}\rightarrow \mathbb {C}\times \{1\}, \left( \mathchoice{\mathbf {x}_0}{{\mathbf {x}_0}}{{\mathbf {x}_0}}{{\mathbf {x}_0}}\right) \mapsto \left( \mathchoice{\mathbf {x}_0, 1}{{\mathbf {x}_0, 1}}{{\mathbf {x}_0, 1}}{{\mathbf {x}_0, 1}}\right)&\text {and}&\end{aligned}$$
(8a)
$$\begin{aligned} i_2&:\mathbb {C}\rightarrow \mathbb {C}^3 , \mathbf {x}_0\mapsto \left( \mathchoice{\mathbf {x}_0, \mathbf {x}_0^1, \mathbf {x}_0^2}{{\mathbf {x}_0, \mathbf {x}_0^1, \mathbf {x}_0^2}}{{\mathbf {x}_0, \mathbf {x}_0^1, \mathbf {x}_0^2}}{{\mathbf {x}_0, \mathbf {x}_0^1, \mathbf {x}_0^2}}\right) .&\end{aligned}$$
(8b)

We also define the projections

$$\begin{aligned} \mathrm{pr}_1&:\mathbb {C}\times \{1\} \rightarrow \mathbb {C}, \left( \mathchoice{\mathbf {x}, 1}{{\mathbf {x}, 1}}{{\mathbf {x}, 1}}{{\mathbf {x}, 1}}\right) \mapsto \mathbf {x}&\text {and}&\end{aligned}$$
(9a)
$$\begin{aligned} \mathrm{pr}_2&:\mathbb {C}^3 \rightarrow \mathbb {C}, \left( \mathchoice{\mathbf {x}, \mathbf {x}^1, \mathbf {x}^2}{{\mathbf {x}, \mathbf {x}^1, \mathbf {x}^2}}{{\mathbf {x}, \mathbf {x}^1, \mathbf {x}^2}}{{\mathbf {x}, \mathbf {x}^1, \mathbf {x}^2}}\right) \mapsto \mathbf {x}. \end{aligned}$$
(9b)

Table 2 shows the generators of the Lie algebra \(\tilde{\mathfrak g}\), as well as their exponentials, which are elements of the Lie group \(\tilde{G}\).

Table 2 Elements of the Lie algebra \(\tilde{\mathfrak g}\) and the associated Lie group \(\tilde{G}\)

We note that the Lie algebra \(\tilde{\mathfrak g}\) can be extended to a larger Lie algebra \(\mathfrak s := \tilde{\mathfrak g} \oplus _\mathrm{Id} \mathfrak t(3, \mathbb {C})\), which includes translations. Table 3 shows the bases of the vector space of the Lie algebra \(\mathfrak t(3, \mathbb {C}) = \mathfrak t_1 \oplus \mathfrak t_2\oplus \mathfrak t_3\) and their respective exponentials.

Table 3 Elements of the Lie algebra \({\mathfrak t}\) and the associated Lie group T

For matrix Lie algebras, the multiplication is given by the canonical matrix multiplication, and the commutator [XY ] is defined as \( [X,Y ] := XY - YX\). A Lie algebra \(\mathfrak g\) is called abelian if its commutator vanishes identically, i.e. if \( [X,Y ] = \mathbf 0\) for all \(X, Y \in \mathfrak g\). If a Lie algebra \(\mathfrak g\) is abelian, then the elements of its associated Lie group commute, i.e. for \(X, Y \in \mathfrak g\), we have \(\mathrm{{e}}^X \mathrm{{e}}^Y = \mathrm{{e}}^Y \mathrm{{e}}^X\) and therefore \(\mathrm{{e}}^X \mathrm{{e}}^Y = \mathrm{{e}}^{X+Y}\). We note that the Lie algebra \(\tilde{\mathfrak g}\) has the abelian subalgebras \(\tilde{\mathfrak g}_1\) and \(\tilde{\mathfrak g}_2\), giving rise to their respective commutative Lie groups \(G_1\) and \(G_2\), and that the Lie algebras \(\mathfrak s_1 := \tilde{\mathfrak g}_1 \oplus _\mathrm{Id} \mathfrak t_1\) and \(\mathfrak s_2 := \tilde{\mathfrak g}_2 \oplus _\mathrm{Id} \mathfrak t_2\) are abelian and constitute subalgebras of the Lie algebra \(\mathfrak s\). Neither of the Lie algebras \(\tilde{\mathfrak g}\) and \(\mathfrak s\), however, is abelian.

Thus, at the border of a finite element, the deformation function assumes values belonging to commutative subgroups of the Lie group \(\tilde{G}\), while its values in the interior of a finite element generally do not commute. It is this property, arising from the specific embedding of the two abelian subalgebras into a larger Lie algebra, in which the interpolation takes place, that results in a strong coupling of rotational and translational components of the deformation function in the interior of a finite element, while preserving the separability of the components on its borders.

In future work, we intend to demonstrate that most of the steps in the subsequent calculations can be readily applied to deformation functions based on the larger Lie algebra \(\mathfrak s\), which includes the translational Lie algebra \(\mathfrak t(3, \mathbb {C})\). In order to keep the calculations as simple as possible, however, we restrict the following exposition to the case of rotations and dilations, i.e. to a deformation function based on the Lie algebra \(\tilde{\mathfrak g}\). In this context, we note that our current research on a model involving translations has shown that in order to obtain good approximations of solutions involving large, simultaneous translations and rotations, the application of a co-rotational approach is necessary. The restricted model presented in the subsequent exposition, however, can be formulated without that additional layer of complexity.

For a given value of \(\xi \), the deformation function \(g_\xi \) is an element of the Lie group G. Given the values \({N_{I,r,q}} (\xi )\) of the shape functions at the position \(\xi \) and the degrees of freedom \(u_{I,r,q}\), which may be assembled into a vector of d.o.f. \(\mathbf u \in U \simeq \mathbb {R}^{4n_{q}}\), the deformation function results from the application of the map

$$\begin{aligned} \psi _1 :U \rightarrow {C^1} ({\mathbb {R}{}, \mathfrak g}),\quad \mathbf u \mapsto \varphi _1 :\mathbb {R}{} \rightarrow \mathfrak g,\quad \xi \mapsto \sum \limits _{\begin{array}{c} 1 \le I\le {{nel}}\\ 1 \le r\le 2\\ 0 \le q\le n_{q}-1 \end{array}} u_{I,r,q} {N_{I,r,q}} (\xi ) {\mathbf {Z}}_{I,r}, \end{aligned}$$
(10)

followed by subsequent exponentiation. Similarly, a deformation function \(\tilde{g}_\xi \) results from the map

$$\begin{aligned} \psi _2 :U \rightarrow {C^1} ({\mathbb {R}{}, \tilde{\mathfrak g}}), \quad \mathbf u \mapsto \varphi _2 :\mathbb {R}{} \rightarrow \tilde{\mathfrak g},\quad \xi \mapsto \sum \limits _{\begin{array}{c} 1 \le I\le {{nel}}\\ 1 \le r\le 2\\ 0 \le q\le n_{q}-1 \end{array}} u_{I,r,q} {N_{I,r,q}} (\xi ) \tilde{\mathbf {Z}}_{I,r}, \end{aligned}$$
(11)

and subsequent exponentiation. By elementary calculations, it can be shown that, for identical shape functions \({N_{I,r,q}} (\xi )\) and initial configuration \({\mathbf {x}}_0({\xi , \eta })\), the results of the action of \(g_\xi := \exp \left( {{{\psi _1} ({\mathbf u})}(\xi )}\right) \) and \(\tilde{g}_\xi := \exp \left( {{{\psi _2}({\mathbf u})}(\xi )}\right) \) on the respective embeddings \(i_1\) resp. \(i_2\) of the initial configuration result in the same current configuration [25]. That is, we have

$$\begin{aligned} \mathbf {x}({\xi , \eta }) \equiv \mathrm{pr}_1 \circ \underbrace{\exp \left( {{{\psi _1}({\mathbf u})} (\xi )}\right) }_{= g_\xi } \circ \, i_1 \circ {\mathbf {x}}_0({\xi , \eta }) \equiv \mathrm{pr}_2 \circ \underbrace{\exp \left( {{{\psi _2}({\mathbf u})}(\xi )}\right) }_{= \tilde{g}_\xi } \circ \, i_2 \circ {\mathbf {x}}_0({\xi , \eta })\quad \end{aligned}$$
(12)

for every \({\mathbf {x}}_0\) and every \(\mathbf u\). The calculations are available from the authors upon request.

The Logarithmic finite element method does not depend on the application of the isoparametric concept. Therefore, the functionals on which the construction of the global finite element system is based are taken from the deformation function rather than from the interpolant of the current configuration, and we do not restrict the initial configuration other than by the condition that it be an immersion of the parametrization space into the physical space. Without loss of generality, we assume, in the remainder of this section, that the first node of the beam element is located at the origin, i.e. \({\mathbf {x}}_0({0, 0}) = \mathbf 0\).

Conditions related to the first derivative of the deformation function

Due to the underlying Bernoulli kinematics, the orientation of the cross-section depends solely on the orientation of the neutral axis of the beam. We also note that, as a result of the application of constitutive equations for beams with invariant cross-sections, the dilatational component of the deformation of the cross-section implied by the construction of the deformation function does not enter into the evaluation of the internal energy. The derivative of the current configuration \(\mathbf {x}\) with regard to the parameterization variable of the neutral axis, \(\xi \), is given by

$$\begin{aligned} {\mathinner { \frac{\partial {^{}}\mathbf {x}}{\partial {\xi ^{}}} }} ({\xi , \eta }) = \mathrm{pr}_2 \circ {\mathinner { \frac{\partial {^{}}\exp {({\psi _2}({\mathbf {u}}))}}{\partial {\xi ^{}}} }} (\xi )&\circ&i_2 \circ {\mathbf {x}}_0({\xi , \eta }) \nonumber \\&\quad + \,\mathrm{pr}_2 \circ \exp {({\psi _2}({\mathbf {u}})(\xi ))} \circ i_2 \circ {\mathinner { \frac{\partial {^{}}{\mathbf {x}}_0}{\partial {\xi ^{}}} }}({\xi ,\eta }), \end{aligned}$$
(13)

Given an appropriate choice of shape functions \({N_{I,r,q}} (\xi )\), this equation can be reduced to a much simpler expression. In particular, we will impose the following restrictions on the basis functions:

$$\begin{aligned} {N_{I,r,q}} (0)&\ne 0 \qquad \text {for nodes and indices} \left( \mathchoice{I, q}{{I, q}}{{I, q}}{{I, q}}\right) = \left( \mathchoice{1,0}{{1,0}}{{1,0}}{{1,0}}\right) , \end{aligned}$$
(14a)
$$\begin{aligned} {N_{I,r,q}} (0)&= 0 \qquad \text {for all nodes } I \text { and indices } q, \text {except for} \left( \mathchoice{I, q}{{I, q}}{{I, q}}{{I, q}}\right) = \left( \mathchoice{1,0}{{1,0}}{{1,0}}{{1,0}}\right) , \end{aligned}$$
(14b)
$$\begin{aligned} {\mathinner { \frac{\partial {^{}}N_{I,r,q}}{\partial {\xi ^{}}} }} (0)&= 0 \qquad \text {for node} \,\,I = 2. \end{aligned}$$
(14c)

With regard to the derivative of the exponential of a matrix-valued function, we note that for a differentiable function \(X(t) :\mathbb {R}{} \rightarrow \mathfrak h\), whose codomain is the Lie algebra \(\mathfrak h\), we have

$$\begin{aligned} {\mathinner { \frac{\partial {^{}}\mathrm{{e}}^{X(t)}}{\partial {t^{}}} }}({t^*}) = \mathrm{{e}}^{X({t^*})} \sum _{k=0}^\infty \frac{({-1})^k}{({k+1})!} {\mathrm{ad}_{X ({t^*})}}^k {\mathinner { \frac{\partial {^{}}X}{\partial {t^{}}} }}({t^*}). \end{aligned}$$
(15)

This equation results from transforming the fractional expression involving the exponential of the adjoint map given in [13, p. 71],

$$\begin{aligned} {\mathinner { \frac{\partial {^{}}\mathrm{{e}}^{X(t)}}{\partial {t^{}}} }}({t^*}) = \mathrm{{e}}^{X({t^*})} \left( {{{\left( {\mathrm{ad}_{X({t^*})}{}}^{-1}{} \left( {\mathbf {I}- \exp {\left( -\mathrm{ad}_{X({t^*})}\right) }}\right) \right) }} \left( {{\mathinner { \frac{\partial {^{}}X}{\partial {t^{}}} }}({t^*})}\right) }\right) , \end{aligned}$$
(16)

into a power series. In this equation, we use the adjoint operator \(\mathrm{ad}\), which, for \(X, Y \in \mathfrak h\), is given by

$$\begin{aligned} \mathrm{ad}:\mathfrak h \rightarrow {GL \,}({\mathfrak h}), \quad X \mapsto \mathrm{ad}_X :\mathfrak h \rightarrow \mathfrak h,\quad Y \mapsto \left[ X,Y\right] , \end{aligned}$$
(17)

where [XY] is the commutator as defined above. Furthermore, \({\mathrm{ad}_X}^k\) denotes the repeated application of \(\mathrm{ad}_X\). Thus, if X and Y are elements of a Lie algebra \(\mathfrak h\), the power series in (15), as well as the sum of any subset of its summands, evaluates to an element of \(\mathfrak h\). In particular, if \(\mathfrak h\) is an abelian Lie algebra, then Eq. (15) simplifies to

$$\begin{aligned} {\mathinner { \frac{\partial {^{}}\mathrm{{e}}^{X(t)}}{\partial {t^{}}} }}({t^*})= \mathrm{{e}}^{X({t^*})} {\mathinner { \frac{\partial {^{}}X}{\partial {t^{}}} }} ({t^*}). \end{aligned}$$
(18)

Given the restrictions in (), we observe that \({{\psi _2}({\mathbf u})}(0) \in \tilde{\mathfrak g}_1\), \({\mathinner { \frac{\partial {^{}}{\psi _2}({\mathbf u})}{\partial {\xi ^{}}} }} (0) \in \tilde{\mathfrak g}_1\), and therefore, as \(\tilde{\mathfrak g}_1\) is abelian,

$$\begin{aligned} {\mathinner { \frac{\partial {^{}}\exp \left( {{\psi _2}({\mathbf u})}\right) }{\partial {\xi ^{}}} }} (0) = \exp \left( {{{\psi _2}({\mathbf u})}(0)}\right) {\mathinner { \frac{\partial {^{}}{\psi _2} ({\mathbf u})}{\partial {\xi ^{}}} }} (0). \end{aligned}$$
(19)

As \({\mathbf {x}}_0({0, 0}) = \mathbf 0\) and \({\mathbf {x}}_0^1 = 0\), for \(\xi = 0\), the first summand in Eq. (13) is of the form

$$\begin{aligned} \mathrm{pr}_2 \circ \begin{pmatrix} z &{} \quad -z &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 \\ \end{pmatrix} \begin{pmatrix} 0 \\ 0 \\ {\mathbf {x}}_0^2 \end{pmatrix}, \end{aligned}$$
(20)

with \(z \in \mathbb {C}{}\), and thus vanishes. Given the restriction (14b), the remaining terms in Eq. (13) simplify to

$$\begin{aligned} {\mathinner { \frac{\partial {^{}}\mathbf {x}}{\partial {\xi ^{}}} }}({0, 0}) = \exp \left( {\sum _{1 \le r\le 2} u_{1,r,0} {N_{1,r,0}} (0) z_{1,r}}\right) \circ {\mathinner { \frac{\partial {^{}}{\mathbf {x}}_0}{\partial {\xi ^{}}} }}({0,0}), \end{aligned}$$
(21)

With \(\mathbf {x}_\xi = \mathinner { \frac{\partial {^{}}\mathbf {x}}{\partial {\xi ^{}}} }\), \({{\mathbf {x}}_0}_{{\xi }} = \mathinner { \frac{\partial {^{}}{\mathbf {x}}_0}{\partial {\xi ^{}}} }\), the main branch of the logarithm of the directional derivative of the deformation function, \({{\mathrm{\mathrm{Log}}}}\circ {\mathinner { \frac{\partial {^{}}\mathbf {x}_\xi }{\partial {{{\mathbf {x}}_0}_\xi ^{}}} }} ({0, 0}) = \sum \nolimits _{1 \le r\le 2} u_{1,r,0} {N_{1,r,0}} (0) z_{1,r}\), a (nonlinear) functional of the deformation function, is an element of the Lie algebra \(\mathfrak g_1\) and, as \({N_{1,r,0}} (0) \ne 0\) for \(r\in \{1,2\}\) due to condition (14a), depends linearly on the degrees of freedom \(u_{1, r, 0}\), \(r \in \{1,2\}\) in a neighborhood \(\mathcal U\) of the origin. We note that \(\mathcal U\) is a strip of width \(2\pi \) in the space of the degrees of freedom U, due to the ambiguous nature of the logarithm function.

Thus, in \(\mathcal U\), there is a one-to-one relationship between the parameters \(u_{1, r, 0}\) in the deformation function of a single finite element and the respective functional of the deformation function at the border of the finite element. Therefore, by defining continuity conditions based on these functionals, the respective parameters can be linked together at the global level and thus serve as global degrees of freedom of a finite element system composed of LogFE-based finite elements, as well as other compatible finite elements.

Conditions related to the second derivative of the deformation function

If one intends to use not only the first derivative of the current configuration w.r.t. the parameterization of the neutral axis, but also its second derivative, then additional restrictions must be imposed on the shape functions. We impose the following additions restrictions on the shape functions:

$$\begin{aligned} {\mathinner { \frac{\partial {^{}}N_{I,r,q}}{\partial {\xi ^{}}} }} (0)&\ne 0 \qquad \text {for nodes and indices} \left( \mathchoice{I, q}{{I, q}}{{I, q}}{{I, q}}\right) = \left( \mathchoice{1,1}{{1,1}}{{1,1}}{{1,1}}\right) , \end{aligned}$$
(22a)
$$\begin{aligned} {\mathinner { \frac{\partial {^{}}N_{I,r,q}}{\partial {\xi ^{}}} }} (0)&= 0 \qquad \text {for all nodes}~ I ~\text {and indices}~ q, \text {except for}~\left( \mathchoice{I, q}{{I, q}}{{I, q}}{{I, q}}\right) = \left( \mathchoice{1,1}{{1,1}}{{1,1}}{{1,1}}\right) , \end{aligned}$$
(22b)
$$\begin{aligned} {\mathinner { \frac{\partial {^{2}}N_{I,r,q}}{\partial {\xi ^{2}}} }} (0)&= 0 \qquad \text {for node}\, I = 2. \end{aligned}$$
(22c)

The second derivative of the current configuration w.r.t. the parameterization of the neutral axis is given by

$$\begin{aligned} {\mathinner { \frac{\partial {^{2}}\mathbf {x}}{\partial {\xi ^{2}}} }} ({\xi , \eta })&= \mathrm{pr}_2 \circ {\mathinner { \frac{\partial {^{2}}\exp ({{\psi _2}({\mathbf u})})}{\partial {\xi ^{2}}} }} (\xi ) \circ i_2 \circ {\mathbf {x}}_0({\xi , \eta }) \nonumber \\&\quad + 2 \mathrm{pr}_2 \circ {\mathinner { \frac{\partial {^{}}\exp ({{\psi _2}({\mathbf u})})}{\partial {\xi ^{}}} }} (\xi ) \circ i_2 \circ \mathinner { \frac{\partial {^{}}{\mathbf {x}}_0({\xi , \eta })}{\partial {\xi ^{}}} } \nonumber \\&\quad + \mathrm{pr}_2 \circ \exp ({{{\psi _2}({\mathbf u})}(\xi )}) \circ i_2 \circ {\mathinner { \frac{\partial {^{2}}{\mathbf {x}}_0}{\partial {\xi ^{2}}} }}({\xi ,\eta }), \end{aligned}$$
(23)

With restriction (22c), we have \({\mathinner { \frac{\partial {^{2}}{\psi _2}({\mathbf u})}{\partial {\xi ^{2}}} }} (0) \in \tilde{\mathfrak g}_1\) and thus, \({\mathinner { \frac{\partial {^{2}}\exp ({{\psi _2}({\mathbf u})})}{\partial {\xi ^{2}}} }} (0) \in \tilde{\mathfrak g}_1\). As for the first derivative, the first summand in Eq. (23) vanishes for \(\xi = 0\), i.e. \(\mathrm{pr}_2 \circ {\mathinner { \frac{\partial {^{2}}\exp ({{\psi _2}({\mathbf u})})}{\partial {\xi ^{2}}} }} (0) \circ i_2 \circ {\mathbf {x}}_0({0, 0}) = \mathbf 0\). With restriction (22b), we obtain

$$\begin{aligned} {\mathinner { \frac{\partial {^{2}}\mathbf {x}}{\partial {\xi ^{2}}} }} ({0, 0})&= 2\mathrm{pr}_2 \circ {\mathinner { \frac{\partial {^{}} \exp \left( { {\psi _2} ({\mathbf u}})\right) }{\partial {\xi ^{}}} }} (0) \circ i_2 \circ {\mathinner { \frac{\partial {^{}}{\mathbf {x}}_0}{\partial {\xi ^{}}} }}({0, 0}) \nonumber \\&\quad + \mathrm{pr}_2 \circ { \exp { \left( {\psi _2} ({\mathbf u})\right) }} (0) \circ i_2 \circ {\mathinner { \frac{\partial {^{2}}{\mathbf {x}}_0}{\partial {\xi ^{2}}} }}({0, 0}) \nonumber \\&= \mathrm{pr}_2 \circ \exp {\left( { {\psi _2} ({\mathbf u})} (0)\right) } \left( {2{\mathinner { \frac{\partial {^{}} {\psi _2} ({\mathbf u})}{\partial {\xi ^{}}} }} (0) \circ i_2 \circ {\mathinner { \frac{\partial {^{}}{\mathbf {x}}_0}{\partial {\xi ^{}}} }}({0, 0}) + i_2 \circ {\mathinner { \frac{\partial {^{2}}{\mathbf {x}}_0}{\partial {\xi ^{2}}} }}({0, 0})}\right) \nonumber \\&= {\mathinner { \frac{\partial {^{}}\mathbf {x}_\xi }{\partial {{{\mathbf {x}}_0}_\xi ^{}}} }} ({0, 0}) \left( {2\sum _{1 \le r\le 2} \left( {u_{1,r,1} {\mathinner { \frac{\partial {^{}}N_{1,r,1}}{\partial {\xi ^{}}} }} (0) z_{1,r}}\right) {\mathinner { \frac{\partial {^{}}{\mathbf {x}}_0}{\partial {\xi ^{}}} }}{(0, 0)} + {\mathinner { \frac{\partial {^{2}}{\mathbf {x}}_0}{\partial {\xi ^{2}}} }}({0, 0})}\right) . \nonumber \\ \end{aligned}$$
(24)

Based on the well-known definitions of the (geometric) curvature \(\kappa \) of the current configuration and \(\kappa _0\) of the initial configuration, we define the “material curvature” \(\kappa ^\mathrm {mat}\) of the current configuration and \(\kappa _0^\mathrm {mat}\) of the initial configuration of the neutral axis, as

$$\begin{aligned} \kappa ^\mathrm {mat}&:= \frac{\vert \partial _\xi \mathbf {x} \vert }{\vert \partial _\xi \mathbf {X} \vert } \kappa = \frac{\vert \partial _\xi \mathbf {x}\times \partial ^2_\xi \mathbf {x} \vert }{\vert \partial _\xi \mathbf {x} \vert ^2 \vert \partial _\xi \mathbf {X} \vert },&\kappa _0^\mathrm {mat}&:= \frac{\vert \partial _\xi {\mathbf {x}}_0 \vert }{\vert \partial _\xi \mathbf {X} \vert } \kappa _0 = \frac{\vert \partial _\xi {\mathbf {x}}_0\times \partial ^2_\xi {\mathbf {x}}_0 \vert }{\vert \partial _\xi {\mathbf {x}}_0 \vert ^2 \vert \partial _\xi \mathbf {X} \vert }. \end{aligned}$$
(25)

We denote \(s := {\mathinner { \frac{\partial {^{}}\mathbf {x}_\xi }{\partial {{{\mathbf {x}}_0}_\xi ^{}}} }} ({0, 0}) \in \mathbb {C}{}\). Then, by inserting the derivatives from Eqs. (13) and (24), followed by elementary calculations, we obtain

$$\begin{aligned} {\kappa _\mathrm {mat}}({0, 0})= & {} \frac{\vert s \vert ^2 {\left( 2\mathrm{Im} \left( {\sum _{1 \le r\le 2} \left( {u_{1,r,1} {\partial _\xi N_{1,r,1}} (0) z_{1,r}}\right) }\right) \vert {\partial _\xi {\mathbf {x}}_0} ({0,0}) \vert ^2 + \vert {\partial _\xi {\mathbf {x}}_0} ({0,0}) \times {\partial ^2_\xi {\mathbf {x}}_0} ({0,0}) \vert \right) }}{\vert s \vert ^2 \vert {\partial _\xi {\mathbf {x}}_0}({0,0}) \vert ^2 \vert {\partial _\xi \mathbf {X}} ({0,0}) \vert } \nonumber \\= & {} 2\frac{u_{1,1,1} {\partial _\xi N_{1,1,1}} (0) }{\vert {\partial _\xi \mathbf {X}} ({0,0}) \vert } + \kappa _0^\mathrm {mat}. \end{aligned}$$
(26)

Therefore, the change in the material curvature of the current configuration relative to the initial configuration is a linear function in the degree of freedom \(u_{1,1,1}\), and condition (22a) ensures that this function is non-trivial. Obviously, care must be taken to ensure that different finite elements, which may not have the same parameterization with regard to their material points, are linked through appropriately formulated continuity conditions w.r.t. their material curvature.

We define the “material derivative of strain” \(\dot{\varepsilon }^\mathrm {mat}\) of the current configuration and \(\dot{\varepsilon }_0^\mathrm {mat}\) of the initial configuration of the neutral axis, as

$$\begin{aligned} \dot{\varepsilon }^\mathrm {mat}&= \frac{\left\langle \partial _\xi \mathbf {x},\partial ^2_\xi \mathbf {x}\right\rangle }{\vert \partial _\xi \mathbf {x} \vert ^2 \vert \partial _\xi \mathbf {X} \vert },&\dot{\varepsilon }_0^\mathrm {mat}&= \frac{\left\langle \partial _\xi {\mathbf {x}}_0,\partial ^2_\xi {\mathbf {x}}_0\right\rangle }{\vert \partial _\xi {\mathbf {x}}_0 \vert ^2 \vert \partial _\xi \mathbf {X} \vert }. \end{aligned}$$
(27)

An elementary calculation similar to that performed in Eq. (26) yields

$$\begin{aligned} {\dot{\varepsilon }^\mathrm {mat}}({0, 0}) = 2\frac{u_{1,0,1} {\partial _\xi N_{1,0,1}} (0) }{\vert {\partial _\xi \mathbf {X}} ({0,0}) \vert } + {\dot{\varepsilon }_0^\mathrm {mat}}({0, 0}). \end{aligned}$$
(28)

The change in the material derivative of strain definition of the current configuration relative to the initial configuration is thus a linear and non-trivial function in the degree of freedom \(u_{1,0,1}\).

Table 4 Degrees of freedom and related functionals of the deformation function, for node 1 with \({\mathbf {x}}_0^1 = {\mathbf {x}}_0 \left( 0,0\right) \), and shape functions satisfying conditions () and ()

The choice of shape functions

Table 4 summarizes the degrees of freedom, their respective functionals, which may be used to impose continuity conditions at the global level, and their respective geometric meaning.

A set of shape functions that satisfies the conditions () and () and thus accommodates the application of continuity conditions both with regard to the first and the second derivative of the deformation function, generally does not generate the same deformation space as a set of shape functions that satisfies none or only some of these conditions, and may also have inferior numerical properties. Therefore, the optimal choice of shape functions for a finite element constructed on the basis of the Logarithmic finite element method will depend on the presence of continuity conditions and boundary conditions at its borders.

Fig. 3
figure 3

Polynomials of the shape functions for a two-node beam element. The scalar shape functions \(N_{I,r,q}\) associated with the dilatational degrees of freedom (\(r = 1\)) and with the rotational degrees of freedom (\(r = 2\)) are identical

For a beam element endowed with one shape function per node \(I\) and Lie algebra element \({\mathbf {Z}}_{I,r}\) (Lagrange-type shape functions), the polynomials \(p_{I,r,0} :\alpha \mapsto \alpha ^2\) therefore satisfy the conditions (). For the case of two shape functions per node \(I\) and Lie algebra element \({\mathbf {Z}}_{I,r}\) (augmented Lagrange-type shape functions), the polynomials \(p_{I,r,0} :\alpha \mapsto \alpha ^3\) and \(p_{I,r,1} :\alpha \mapsto \alpha ^4 - \alpha ^3\) satisfy these conditions. Figure 3 displays the Lagrange-type and augmented Lagrange-type shape functions for a two-node beam element. Note that both sets of Langrangian-type shape functions (see Fig. 3) do not satisfy the additional conditions (), and therefore cannot be used to model finite elements subjected to continuity conditions or boundary conditions involving the second derivative of the deformation function. However, the Hermite-type shape functions (see Fig. 3) satisfy both sets of conditions and can therefore be used to construct finite elements that can accommodate all types of continuity and boundary conditions described above (see Table 4). The choice of appropriate shape functions for different boundary conditions is discussed in more detail in Appendix 1.

Consistency with the linear beam theory

A non-linear boundary-value problem is often solved by iteratively obtaining the solutions of a sequence of linear problems, obtaining a sequence of solutions that converges toward the non-linear solution. In this section, we therefore investigate in more detail the linear steps involved in the iterative process of finding a non-linear solution. In particular, we focus on the initial linear step, for which the initial estimate for the degrees of freedom is given by the zero vector, i.e. the initial estimate for the deformation function is the identity function. We will compare the characteristics of this initial linear step with the linear theory of the Bernoulli beam. We will see that for certain standard load cases, the reaction of the beam in a neighborhood of the identity deformation predicted by the linearization of the LogFE model is identical, up to terms of higher order, to the reaction predicted by the linear theory of the Bernoulli beam.

The linearization of the LogFE formulation around the origin, i.e. \(\mathbf u= \mathbf 0\), results in a simplified model that can be regarded as a linear variant of the general approach. In this case, the degrees of freedom assume values that are proportional to the load intensity. As a result, the trajectories of the materials points of the configuration are located on the orbits of their initial locations under the action of a one-parameter Lie group. As we will see, these orbits are located on circles (and straight lines, which can be regarded as degenerate circles) if the dilatation coefficients vanish. In the case of non-zero dilatation coefficients, the orbits are located on logarithmic spirals. In order to analytically determine the linearization for different load cases, the strain and the curvature (and, as a result, the normal force and the bending moment) at a given parameter value \(\xi \) of the configuration must be derived from the deformation function.

Normal force and bending moment

The neutral axis of the beam shall be parameterized by the curve

$$\begin{aligned}\xi \mapsto \chi _0\left( f \left( \xi \right) \right) = \bigg (\begin{array}{l} \xi \\ 1 \end{array}\bigg ), \end{aligned}$$

connecting the nodal positions \({\mathbf {x}}_0^1 = 0\) and \({\mathbf {x}}_0^2 = 1\) on a straight line on the real axis. Thus, we have \({\mathbf {x}}_0\left( \xi \right) = \xi \). In the following, the scalar product and the determinant are based on the Euclidean vectors associated with the complex values, i.e. for \(w, z \in \mathbb {C}{}\), we have, in this context, the scalar product

$$\begin{aligned}\left\langle w,z\right\rangle := \left\langle \bigg (\begin{array}{l} \mathrm{Re}\left( w\right) \\ \mathrm{Im}\left( w\right) \end{array}\bigg ),\bigg (\begin{array}{l} \mathrm{Re}\left( z\right) \\ \mathrm{Im}\left( z\right) \end{array}\bigg )\right\rangle \end{aligned}$$

and the determinant

$$\begin{aligned} \left| w,z\right| := \det \bigg (\begin{array}{ll} \mathrm{Re}\left( w\right) &{} \quad \mathrm{Re}\left( z\right) \\ \mathrm{Im}\left( w\right) &{}\quad \mathrm{Im}\left( z\right) \end{array}\bigg ) . \end{aligned}$$

For the purpose of calculating the strain measures, we identify the complex plane with the Euclidean space, such that a distance of 1 equals 1 m.

With \(\dot{\mathbf {x}}:= {\partial }_{\xi } {\mathbf {x}}\), \(\ddot{\mathbf {x}}:= {\partial }^{2}_{\xi } {\mathbf {x}}\), the strain \(\varepsilon \) is given as \(\varepsilon (\xi ) = \left\| \dot{\mathbf {x}}\right\| - 1\) and the curvature \(\kappa \) is given as \(\kappa (\xi ) = \left\| \dot{\mathbf {x}}\right\| ^{-3} \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}\right| \). The derivatives with respect to a given degree of freedom \(u\) are

$$\begin{aligned} {\partial }_{u} {\varepsilon } (\xi )&= {\left\| \dot{\mathbf {x}}\right\| }^{-1} \left\langle {\partial }_{u} \dot{\mathbf {x}},\dot{\mathbf {x}}\right\rangle \quad \text {and} \end{aligned}$$
(29a)
$$\begin{aligned} {\partial }_{u} {\kappa } (\xi )&= -3 \left\| \dot{\mathbf {x}}\right\| ^{-5} \left\langle {\partial }_{u} \dot{\mathbf {x}},\dot{\mathbf {x}}\right\rangle \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}\right| + \left\| \dot{\mathbf {x}}\right\| ^{-3} \left( {\left| {\partial }_{u} \dot{\mathbf {x}}, \ddot{\mathbf {x}}\right| + \left| \dot{\mathbf {x}},{\partial }_{u} \ddot{\mathbf {x}}\right| }\right) . \end{aligned}$$
(29b)

Thus, for \(\left\| \dot{\mathbf {x}}\right\| = 1\), \(\left\| \ddot{\mathbf {x}}\right\| = 0\), we obtain the derivative of the strain as \({\partial }_{u} {\varepsilon } (\xi ) = \left\langle {\partial }_{u} \dot{\mathbf {x}},\dot{\mathbf {x}}\right\rangle \) and the derivative of the curvature as \({\partial }_{u} {\kappa }(\xi ) = \left| \dot{\mathbf {x}},{\partial }_{u} \ddot{\mathbf {x}}\right| \). For \(\mathbf u= \mathbf 0\), we have \({\bar{\mathbf {Z}}}(\xi ) \equiv \mathbf 0\), \({{e}}^{{\bar{\mathbf {Z}}}(\xi )} \equiv \mathrm{Id}\), from which follows \(\left[ \partial _u{\bar{\mathbf {Z}}}(\xi ),{\bar{\mathbf {Z}}}(\xi )\right] \equiv \mathbf 0\) and therefore \(\partial _u{e}^{{\bar{\mathbf {Z}}}(\xi )} \equiv \partial _u{\bar{\mathbf {Z}}}(\xi )\). We obtain

$$\begin{aligned} \dot{\mathbf {x}}&= \partial _\xi \mathrm{{e}}^{{\bar{\mathbf {Z}}}(\xi )} \begin{pmatrix} {\mathbf {x}}_0(\xi ) \\ 1 \end{pmatrix} + \mathrm{{e}}^{{\bar{\mathbf {Z}}}(\xi )} \begin{pmatrix} 1 \\ 0 \end{pmatrix}= \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \end{aligned}$$
(30a)
$$\begin{aligned} \partial _u\dot{\mathbf {x}}&= \partial _u\partial _\xi {\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} {\mathbf {x}}_0(\xi ) \\ 1 \end{pmatrix} + \partial _u{\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \end{aligned}$$
(30b)
$$\begin{aligned} \partial _u\ddot{\mathbf {x}}&= \partial _u\partial ^2_\xi {\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} {\mathbf {x}}_0(\xi ) \\ 1 \end{pmatrix} + 2\partial _u\partial _\xi {\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} 1 \\ 0 \end{pmatrix}. \end{aligned}$$
(30c)

Thus, the derivatives of the strain and the curvature are given by

$$\begin{aligned} \partial _u\varepsilon (\xi )&= \left\langle \partial _u\partial _\xi {\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} {\mathbf {x}}_0(\xi ) \\ 1 \end{pmatrix} + \partial _u{\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} 1 \\ 0 \end{pmatrix},\begin{pmatrix} 1 \\ 0 \end{pmatrix}\right\rangle \quad \text {and} \end{aligned}$$
(31a)
$$\begin{aligned} \partial _u\kappa (\xi )&= \left| \begin{pmatrix} 1 \\ 0 \end{pmatrix},\partial _u\partial ^2_\xi {\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} {\mathbf {x}}_0(\xi ) \\ 1 \end{pmatrix} + 2\partial _u\partial _\xi {\bar{\mathbf {Z}}}(\xi ) \begin{pmatrix} 1 \\ 0 \end{pmatrix}\right| , \end{aligned}$$
(31b)

respectively. With \(\eta _1 = 1\) (dilatation), \(\eta _2 = \mathrm{i}\) (rotation), the elements of the Lie algebra are given as \({\mathbf {Z}}_{I,r,q} = \eta _r \bigg (\begin{array}{cc} 1 &{} \quad - {\mathbf {x}}_0^I\\ 0 &{}\quad 0\end{array}\bigg )\). For the derivatives of the strain and the curvature with respect to the degrees of freedom \(\mathbf u\), evaluated at \(\mathbf u= 0\), we therefore obtain

$$\begin{aligned} \partial _{u_{I,r,q}} {\varepsilon (\xi )}&= \mathrm{Re} \left( \eta _r \left( \partial _\xi {N_{I,r,q}} (\xi ) \left( {{\mathbf {x}}_0(\xi ) - {\mathbf {x}}_0^I} \right) + {N_{I,r,q}} (\xi )\right) \right) \quad \text {and} \end{aligned}$$
(32a)
$$\begin{aligned} \partial _{u_{I,r,q}} {\kappa (\xi )}&= \mathrm{Im} \left( {\eta _r \left( {\partial ^2_\xi {N_{I,r,q}} (\xi ) \left( {{\mathbf {x}}_0(\xi ) - {\mathbf {x}}_0^I}\right) + 2\partial _\xi {N_{I,r,q}} (\xi )}\right) }\right) . \end{aligned}$$
(32b)

Equations (32a) and (32b) indicate that the linearized strain depends on the dilatation coefficients only, while the linearized curvature depends solely on the rotation coefficients. For a two-node, simply supported beam element without internal degrees of freedom, i.e. \(\bar{Q}_{I} = 0\), the differentials of the strain and of the curvature at the origin are given by

$$\begin{aligned} {{\mathrm{d \!}}}\varepsilon (\xi )&= \sum \limits _{I\in \{1,2\}} \left( {\partial _\xi {N_{I,1,0}} (\xi ) \left( {{\mathbf {x}}_0(\xi ) - {\mathbf {x}}_0^I}\right) + {N_{I,1,0}} (\xi )} \right) {{\mathrm{d \!}}}u_{I,1,0}, \end{aligned}$$
(33a)
$$\begin{aligned} {{\mathrm{d \!}}}\kappa (\xi )&= \sum \limits _{I\in \{1,2\}} \left( {\partial ^2_\xi {N_{I,2,0}} (\xi ) \left( {{\mathbf {x}}_0(\xi ) - {\mathbf {x}}_0^I}\right) + 2\partial _\xi {N_{I,2,0}} (\xi )} \right) {{\mathrm{d \!}}}u_{I,2,0}. \end{aligned}$$
(33b)

Linearization based on rotation coefficients

For a simply supported beam, subjected to an external moment \(M_0\) at node 1, we have

$$\begin{aligned} {N_{1,2,0}}(\xi )&= {p_{1,2,0}}\left( {{\psi _1}(\xi )}\right) = ({1-\xi })^3 = -\xi ^3 + 3\xi ^2 - 3\xi +1 , \end{aligned}$$
(34)
$$\begin{aligned} {N_{2,2,0}} (\xi )&= {p_{2,2,0}}\left( {{\psi _2} (\xi )}\right) = -2\xi ^3 + 3\xi ^2. \end{aligned}$$
(35)

Therefore,

$$\begin{aligned} {{\mathrm{d \!}}}\kappa (\xi ) =\varphi _{1}(\xi ) {{\mathrm{d \!}}}u_{1,2,0}+\varphi _{2}(\xi ) {{\mathrm{d \!}}}u_{2,2,0}, \end{aligned}$$
(36)

with

$$\begin{aligned} \varphi _{1}(\xi )=-12\xi ^2+18\xi -6,\end{aligned}$$
(37a)
$$\begin{aligned} \varphi _{2}(\xi )=-24\xi ^2+30\xi -6. \end{aligned}$$
(37b)

If the rotation coefficients are set to zero, i.e. \(u_{1,1,1} = 0\), \(u_{2,1,1} = 0\), the linearized strain vanishes, as \(\eta _2 = \mathrm{i}\) and therefore, as a result of Eq. 32a, \(\partial _u\varepsilon (\xi ) = 0\). A linearization of the total energy U at the initial configuration, with vanishing initial curvature \(\kappa _{0}(\xi )\equiv 0\), is given by

$$\begin{aligned} U=\frac{1}{2} EI\int _{0}^{1}\kappa (\xi )^2 \left\| \partial _{\xi }\mathbf {x}_0 \right\| {\mathrm d}\xi -M_{0}u_{1,2,0}, \end{aligned}$$
(38)

with \(\kappa (\xi )=(\varphi _{1}(\xi )u_{1,2,0}+\varphi _{2}(\xi )u_{2,2,0}) {\mathrm m}^{-1}\). The equilibrium condition is given by \(\partial _{u_{1,2,0}}U=0\) for all I, resulting in \(u_{1,2,0}=\frac{1}{3}(EI)^{-1}M_{0}\hbox {m}, u_{2,2,0}=-\frac{1}{6}(EI)^{-1}M_{0}\hbox {m}\). With these values, the curvature is given by \(\kappa (\xi )=(\xi -1)(EI)^{-1}M_{0}\). These results are consistent with the linear theory of the Bernoulli beam.

In a finite element formulation based on the LogFE method, the function

$$\begin{aligned} {g_0(\xi )} \left( {{\mathbf {x}}_0(\xi )}\right)= & {} \exp \left( {\sum \limits _{I,q} \underbrace{s \beta _{I,2,q}}_{=u_{I,2,q}} {N_{I,2,q}} (\xi ) {\mathbf {Z}}_{I,2,q}}\right) {\mathbf {x}}_0(\xi ) \nonumber \\= & {} \mathrm{{e}}^s \underbrace{\exp \left( {\sum \nolimits _{I,q} \beta _{I,2,q} {N_{I,2,q}} (\xi ) {\mathbf {Z}}_{I,2,q}}\right) }_{{\tilde{\mathbf {Z}}} (\xi )} {\mathbf {x}}_0(\xi ), \end{aligned}$$
(39)

with \(I\in \{1,2\}\), \(0 \le q\le Q_{I,r}\), \(\beta _{I,2,q} \in \mathbb R\) and \(s \in \mathbb {R}{}\) can be thought of as the analogon of the deflection curve in the conventional linear theory. In this formulation, the scalar degree of freedom s is exclusively related to the load intensity. By setting s proportional to the load intensity, based on the sensitivity of s to an infinitesimal increase of the load, we obtain a formula that is fully linearized on the space of the Lie algebra, while still resulting in a nonlinear deformation function in the physical space.

Fig. 4
figure 4

Linearized deformation (only rotational degrees of freedom): deformed configuration and orbits \(\mathbf {x} \left( s\right) \) of points on the neutral axis of the beam

For the present load case, i.e. a simply supported beam, subjected to an external moment at the left node (see also Fig. 4), we have

$$\begin{aligned} N_{1,2,0}&= -3({1 - \xi })^4 + 4({1 - \xi })^3, \end{aligned}$$
(40a)
$$\begin{aligned} N_{1,2,1}&= -3\xi ^4 + 4\xi ^3. \end{aligned}$$
(40b)

With \(\beta _{1,2,0} = 1\), \(\beta _{2,2,0} = -\frac{1}{2}\), we obtain

$$\begin{aligned} \tilde{\mathbf {Z}}:= {N_{1,2,0}} (\xi ) {\mathbf {Z}}_{1,2,0} - \frac{1}{2} {N_{2,2,0}} (\xi ) {\mathbf {Z}}_{2,2,0} = \mathrm{i} \begin{pmatrix} \frac{3}{2} \xi ^2 - 3\xi + 1 &{}\quad - \xi ^3 + \frac{3}{2} \xi ^2 \\ 0 &{}\quad 0 \end{pmatrix}. \end{aligned}$$
(41)

This formula also reveals that the deformation can be expressed as a rotation of each point of the initial configuration, with rotation center and angle depending on the parameter \(\xi \). Thus, the orbits of the material points under the given deformations are located on circles of different radius and location. For \(\frac{3}{2} \xi ^2 - 3\xi + 1 = 0\), i.e. \(\xi = 1 - \frac{1}{3} \sqrt{3}\), the deformation is a pure translation in the z direction.

$$\begin{aligned}&{g_0(\xi )}\left( {{\mathbf {x}}_0(\xi )}\right) \nonumber \\&\quad = {\left\{ \begin{array}{ll} {e}^s \exp \left( {\left( {\frac{3}{2} \xi ^2 - 3\xi + 1}\right) \begin{pmatrix} \mathrm{i} &{}\quad \frac{- \xi ^3 + \frac{3}{2} \xi ^2}{\frac{3}{2} \xi ^2 - 3\xi + 1} \mathrm{i} \\ 0 &{}\quad 0 \end{pmatrix}}\right) \begin{pmatrix} \xi \\ 1 \end{pmatrix} &{}\quad \text {if}\,\, \xi \ne 1 - \frac{1}{3} \sqrt{3}, \\ {e}^s \exp \begin{pmatrix} 0 &{}\quad \left( {- \xi ^3 + \frac{3}{2} \xi ^2} \mathrm{i} \right) \\ 0 &{}\quad 0 \end{pmatrix} \begin{pmatrix} \xi \\ 1 \end{pmatrix} = \begin{pmatrix} 1 - \frac{1}{3} \sqrt{3} \\ 1 \end{pmatrix} + \begin{pmatrix} \frac{1}{9} \sqrt{3} \mathrm{i} \\ 0 \end{pmatrix} s&\quad \text {if} \,\,\xi = 1 - \frac{1}{3} \sqrt{3}. \end{array}\right. }\nonumber \\ \end{aligned}$$
(42)

For a beam simply supported at the left side and clamped at the right side, subjected to an external moment at the left support (see Fig. 4), we omit the shape function \(N_{2,2,0}\) (see Table 8 in Appendix 1), and obtain \(N_{1,2,0} = ({1 - \xi })^2\), \(\beta _{1,2,0} = 1\), \({\partial _s \kappa } (\xi ) = 6\xi - 4\), as predicted by the linear beam theory, and

$$\begin{aligned} {g_0(\xi )} \left( {{\mathbf {x}}_0(\xi )}\right) = {e}^s \exp \left( {\left( {\xi ^2 - 2\xi + 1}\right) \begin{pmatrix} \mathrm{i} &{}\quad 0 \\ 0 &{}\quad 0 \end{pmatrix}}\right) \begin{pmatrix} \xi \\ 1 \end{pmatrix} = \begin{pmatrix} {e}^{s \left( {\xi ^2 - 2\xi + 1 }\right) \mathrm{i} } \xi \\ 1 \end{pmatrix}. \end{aligned}$$
(43)

For a simply supported beam subjected to a constant line load (see Fig. 4), we have \(N_{1,2,0} = -3({1 - \xi })^4 + 4({1 - \xi })^3\), \(N_{2,2,0} = -3\xi ^4 + 4\xi ^3\). With \(\beta _{1,2,0} = 1\), \(\beta _{2,2,0} = -1\), we obtain \({\partial _s \kappa }(\xi ) = 12\xi ^2 - 12\xi \), also consistent with the linear beam theory.

Finally, for a simply supported beam subjected to a linearly decreasing line load (see Fig. 4), we use the same shape functions \(N_{1,2,0} = -3({1 - \xi })^4 + 4({1 - \xi })^3\), \(N_{2,2,0} = -3\xi ^4 + 4\xi ^3\). With \(\beta _{1,2,0} = 1\), \(\beta _{2,2,0} = 1\), we obtain \({\partial _s \kappa }(\xi ) = -120\xi ^3 + 180\xi ^2 - 60\xi \), conforming to the results of the linear beam theory.

These cases demonstrate that for some simple load cases, the tangent space \({T_{\mathbf u}}({g_\xi })\) of the function space of the LogFE formulation at the origin, i.e. for \(\mathbf u= \mathbf 0\), contains the deflection functions obtained by the linear beam theory. Naturally, for more complicated load cases, \({T_{\mathbf u}}({g_\xi })\) will not be large enough to contain the exact linear deflection function. For larger deformations, the linearized logarithmic deformation function, just as the conventional linearized deflection curve, does not adequately represent the actual deformation of the beam.

Fig. 5
figure 5

Linearized deformation (general case): deformed configuration and orbits \(\mathbf {x} \left( s\right) \) of points on the neutral axis of the beam

The general case

In the general case, which includes, e.g., the deformation within a single step of a Newton-Raphson approximation, the orbits of the material points of the configuration are not located on circles, but on logarithmic spirals with different points of origin.

Figure 5 presents an example of the general case. For a diagonally applied constant line load f, resulting from a combination of a line load perpendicular to the neutral axis, q, and a line load parallel to the neutral axis, \(f_N\), we use \(N_{1,2,0} = -3({1 - \xi })^4 + 4({1 - \xi })^3\), \(N_{2,2,0} = -3\xi ^4 + 4\xi ^3\), as above. We set \(\beta _{1,2,0} = \gamma _2\), \(\beta _{2,2,0} = \gamma _2\), \(\gamma _2 \in \mathbb {R}{}\), such that \({\partial _s \kappa }(\xi ) = 12\gamma _1\xi ^2 - 12\gamma _1\xi \), consistent with the linear beam theory. As there is a line load in the direction of the beam, we must also determine the coefficients \(\beta _{1,1,0}\), \(\beta _{2,1,0}\), which are related to the dilatations (\(r= 1\)). Because the impact of the line force in the longitudinal direction on the curvature of the beam vanishes for \(\mathbf u= \mathbf 0\), we can drop the optional restriction \(\partial p_{I,2,0} = 0\) (see Table 8 in Appendix 1). Thus, \(\bar{Q}_I= 0\) for \(I\in \{1,2\}\), and the shape functions related to the dilatations, which result in a change in the normal force in the linear beam theory, are given as \(N_{1,1,0} = ({1 - \xi })^2\), \(N_{2,1,0} = \xi ^2\). From the linear theory, we know that \({\partial _{f_N} \varepsilon }(\xi ) = \frac{1}{EA} ({{-\xi + \frac{1}{2}}})\). Using Eq. 33a, we see that \(\beta _{1,1,0} = \frac{1}{2} \gamma _1\), \(\beta _{2,1,0} = -\frac{1}{2} \gamma _1\), \(\gamma _1 \in \mathbb {R}{}\), result in \({\partial _s \varepsilon } (\xi ) = -\gamma _1\xi + \frac{1}{2} \gamma _1\), conforming to the linear theory.

Introducing the single parameter \(s \in \mathbb {R}{}\), as above, we obtain the orbits of the material points of the beam as

$$\begin{aligned}&{g_0(\xi )} \left( {{\mathbf {x}}_0(\xi )}\right) \nonumber \\&\quad = \mathrm{{e}}^s \exp \begin{pmatrix} \frac{\gamma _1}{s} \left( {-\xi + \frac{1}{2}}\right) + \frac{\gamma _2}{s} \mathrm{i} \left( {4\xi ^3 - 6\xi ^2 + 1}\right) &{} \frac{1}{2} \frac{\gamma _1}{s} \xi ^2 + \frac{\gamma _2}{s} \mathrm{i} \left( {-3\xi ^4 + 4\xi ^3}\right) \\ 0 &{} 0 \end{pmatrix} \cdot \begin{pmatrix} \xi \\ 1 \end{pmatrix}.\nonumber \\ \end{aligned}$$
(44)

In this formula, the coefficients \(\gamma _1\), \(\gamma _2\), depend on the relative strengths of the line loads perpendicular and parallel to the neutral axis, q and \(f_N\), as well as on the axial and flexural rigidity of the beam. The results for \(\gamma _2 = 4\gamma _1\) and \(\gamma _2 = \gamma _1\) are shown in Fig. 5.

Quasi-static equilibrium

In this section, we consider the non-linear deformation of a linear elastic Bernoulli beam, modelled as a single element. The stress resultants, i.e. the axial force N and the bending moment M, are given as \(\varvec{\sigma }(\xi ) = \mathbf {C}\varvec{\epsilon }(\xi )\), with \(\varvec{\sigma }(\xi ) := \bigg (\begin{array}{l} N(\xi ) \\ M(\xi ) \end{array}\bigg )\), the matrix \(\mathbf C := \bigg (\begin{array}{ll} EA &{}\quad 0 \\ 0 &{} \quad EI \end{array}\bigg )\), and the beam strains \(\varvec{\epsilon }(\xi ) := \bigg (\begin{array}{l} \varepsilon (\xi ) \\ \kappa (\xi ) \end{array}\bigg )\). In these definitions, E is Young’s modulus, A is the area of the cross-section, I is the moment of inertia, \(\varepsilon \) denotes the strain and \(\kappa \) the curvature of the neutral axis. The neutral axis of the beam is parameterized by the variable \(\xi \in \Omega \). Partial derivatives are given in index notation. The local internal energy is given as \(\frac{1}{2} {\varvec{\sigma }(\xi )}^{\mathrm{T}}{\varvec{\epsilon }(\xi )}\). External loads, such as distributed loads, concentrated loads, and external moments, are given as \(\mathbf p_{\mathrm{ext}}\). With this, we obtain the residuum \(\mathbf {r}\) as

$$\begin{aligned} \mathbf {r}(\mathbf {u}) = \int _\Omega \varvec{\epsilon }^\mathrm{T} \mathbf {C}\varvec{\epsilon }_{\mathbf {u}}{{\mathrm{d \!}}}\xi - \mathbf p_{\mathrm{ext}}. \end{aligned}$$
(45)

It is important to recall that in the present formulation, the shape functions are linear in the degrees of freedom only on the logarithmic space. With regard to the physical space, the displacements of the material points induced by a change in the degrees of freedom are non-linear. Thus, even for conservative loads, \(\mathbf p_{\mathrm{ext}}\) generally is a non-linear function of the degrees of freedom \(\mathbf {u}\), and, as a result, its derivative \(\mathbf p_{\mathrm{ext},\mathbf {u}}\) does not vanish.

Given a suitable initial estimate for the values of the degrees of freedom, the Newton-Raphson method can be used to obtain a root of Eq. (45). An update of the Newton-Raphson root-finding algorithm is given by

$$\begin{aligned} {\mathbf {u}}_{n+1}= {\mathbf {u}}_n+ \left( {{\mathbf k_{\mathrm{int}}^\mathrm {mat}}({{\mathbf {u}}_n}) + {\mathbf k_{\mathrm{int}}^\mathrm {geo}}({{\mathbf {u}}_n}) - {\mathbf k_{\mathrm{ext}}^\mathrm {geo}}({{\mathbf {u}}_n})}\right) \left( {{\mathbf p_{\mathrm{int}}}({{\mathbf {u}}_n}) + {\mathbf p_{\mathrm{ext}}}({{\mathbf {u}}_n})}\right) , \end{aligned}$$
(46)

with the material stiffness matrix \(\mathbf k_{\mathrm{int}}^\mathrm {mat}:= \int _\Omega \left( {\varvec{\epsilon }^\mathrm{T}}\right) _{\mathbf {u}}\otimes \mathbf {C}\varvec{\epsilon }_{\mathbf {u}}{{\mathrm{d \!}}}\xi \), the geometric stiffness matrix \(\mathbf k_{\mathrm{int}}^\mathrm {geo}:= \int _\Omega \varvec{\epsilon }^\mathrm{T} \mathbf {C}\varvec{\epsilon }_{\mathbf {u}\mathbf {u}} {{\mathrm{d \!}}}\xi \), the derivative of the external loads \(\mathbf k_{\mathrm{ext}}^\mathrm {geo}:= \mathbf p_{\mathrm{ext},\mathbf {u}}\) and the internal forces \(\mathbf p_{\mathrm{int}} := - \int _\Omega \varvec{\epsilon }^\mathrm{T} \mathbf {C}\varvec{\epsilon }_{\mathbf {u}}{{\mathrm{d \!}}}\xi \).

For moderate load intensities, the geometric stiffness matrix, \(\mathbf k_{\mathrm{int}}^\mathrm {geo}\), and the derivative of the external loads, \(\mathbf k_{\mathrm{ext}}^\mathrm {geo}\), are small in relation to the material stiffness matrix. Therefore, one or both of these comparatively computationally expensive geometric components of the derivative of the residuum may be omitted in the root-finding algorithm. In this case, however, the reduction in computation effort during each iteration may be offset by an increase in the number of iterations necessary to obtain a sufficiently good solution, as the root-finding algorithm will no longer show quadratic converge near the exact solution, and the convergence characteristics is likely to be affected negatively throughout the entire iteration process. Nevertheless, if the iteration is taking place in a setting in which quadratic convergence cannot be expected as a result of other causes, such as in the context of a multigrid algorithm, omitting one or both of the terms may be advantageous with regard to the overall efficiency of the calculation.

Beam strains and external loads

In the initial configuration, as well as in the current configuration, the Euclidean plane will be identified with the homogenized complex plane, \(\mathbb {C}{} \times \{1\}\). It is assumed that the beam, in its initial configuration, is represented by a straight line without initial stresses. In the following, we will use the notations \(\bar{\mathbf {x}}:= \left( \mathchoice{\mathrm{Re}\, \mathbf {x}, \mathrm{Im}\, \mathbf {x}}{{\mathrm{Re}\, \mathbf {x}, \mathrm{Im}\, \mathbf {x}}}{{\mathrm{Re}\, \mathbf {x}, \mathrm{Im}\, \mathbf {x}}}{{\mathrm{Re}\, \mathbf {x}, \mathrm{Im}\, \mathbf {x}}}\right) \in \mathbb {R}^2\), \(\dot{\mathbf {x}}:= \bar{\mathbf {x}}_{\xi }\) and \(\ddot{\mathbf {x}}:= \bar{\mathbf {x}}_{\xi \xi }\) for brevity. We use the symbol “\(\bullet \)” in the tensor notation in order to indicate that expressions such as \(\left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}_{\bullet \mathbf {u}}\right\rangle \) involve taking the tensor product, with the scalar product as multiplicative operation. Thus, as \(\dot{\mathbf {x}}_{\mathbf {u}}\) and \(\dot{\mathbf {x}}_{\bullet \mathbf {u}}\) belong to the tensor space \(\left( {\mathbb {R}^2}\right) ^{{ndf}}\), with \({{ndf}}\) denoting the number of degrees of freedom per element, \(\left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}_{\bullet \mathbf {u}}\right\rangle \) is a tensor of order 2, belonging to the tensor space \(\mathbb {R}^{{ndf}}\otimes \mathbb {R}^{{ndf}}\).

Using the notation described above, the strain and its derivatives are given by

$$\begin{aligned} \varepsilon&= \left\| \dot{\mathbf {x}}\right\| - 1, \end{aligned}$$
(47a)
$$\begin{aligned} \varepsilon _{\mathbf {u}}&= {\left\| \dot{\mathbf {x}}\right\| {}}^{-1}{} \left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}\right\rangle , \end{aligned}$$
(47b)
$$\begin{aligned} \varepsilon _{\mathbf {u}\mathbf {u}}&= -\left\| \dot{\mathbf {x}}\right\| ^{-3} \left\langle \dot{\mathbf {x}}_{\bullet \mathbf {u}},\dot{\mathbf {x}}\right\rangle \left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}\right\rangle + {\left\| \dot{\mathbf {x}}\right\| {}}^{-1}{} \left( {\left\langle \dot{\mathbf {x}}_{\mathbf {u}\mathbf {u}},\dot{\mathbf {x}}\right\rangle + \left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}_{\bullet \mathbf {u}}\right\rangle }\right) , \end{aligned}$$
(47c)

the curvature and its derivatives by

$$\begin{aligned} \kappa&= \left\| \dot{\mathbf {x}}\right\| ^{-3} \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}\right| , \end{aligned}$$
(48a)
$$\begin{aligned} \kappa _{\mathbf {u}}&= -3 \left\| \dot{\mathbf {x}}\right\| ^{-5} \left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}\right\rangle \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}\right| + \left\| \dot{\mathbf {x}}\right\| ^{-3} \left( {\left| \dot{\mathbf {x}}_{\mathbf {u}}, \ddot{\mathbf {x}}\right| + \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}_{\mathbf {u}}\right| }\right) , \end{aligned}$$
(48b)
$$\begin{aligned} \kappa _{\mathbf {u}\mathbf {u}}&= 15 \left\| \dot{\mathbf {x}}\right\| ^{-7} \left\langle \dot{\mathbf {x}}_{\bullet \mathbf {u}},\dot{\mathbf {x}}\right\rangle \left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}\right\rangle \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}\right| \nonumber \\&\quad -3 \left\| \dot{\mathbf {x}}\right\| ^{-5} \left( {\left( \left\langle \dot{\mathbf {x}}_{\mathbf {u}\mathbf {u}},\dot{\mathbf {x}}\right\rangle + \left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}_{\bullet \mathbf {u}}\right\rangle \right) } \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}\right| + \left\langle \dot{\mathbf {x}}_{\mathbf {u}},\dot{\mathbf {x}}\right\rangle \left( {\left| \dot{\mathbf {x}}_{\bullet \mathbf {u}}, \ddot{\mathbf {x}}\right| + \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}_{\bullet \mathbf {u}}\right| }\right) \right. \nonumber \\&\quad \left. + \left\langle \dot{\mathbf {x}}_{\bullet \mathbf {u}},\dot{\mathbf {x}}\right\rangle \left( {\left| \dot{\mathbf {x}}_{\mathbf {u}}, \ddot{\mathbf {x}}\right| + \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}_{\mathbf {u}}\right| }\right) \right) \nonumber \\&\quad +\left\| \dot{\mathbf {x}}\right\| ^{-3} \left( {\left| \dot{\mathbf {x}}_{\mathbf {u}\mathbf {u}}, \ddot{\mathbf {x}}\right| + \left| \dot{\mathbf {x}}_{\mathbf {u}}, \ddot{\mathbf {x}}_{\bullet \mathbf {u}}\right| + \left| \dot{\mathbf {x}}_{\bullet \mathbf {u}}, \ddot{\mathbf {x}}_{\mathbf {u}}\right| + \left| \dot{\mathbf {x}}, \ddot{\mathbf {x}}_{\mathbf {u}\mathbf {u}}\right| }\right) . \end{aligned}$$
(48c)

In the case of a distributed load \({\mathbf q}(\xi )\) on an interval \(I \subseteq \Omega \), the external load and its derivative with regard to the degrees of freedom \(\mathbf {u}\) are given by \(\mathbf p_{\mathrm{ext}} = \int _I \left\langle {\bar{\mathbf {x}}_{\mathbf {u}}} (\xi ),{\mathbf q}(\xi )\right\rangle {{\mathrm{d \!}}}\xi \) and \(\mathbf k^\mathrm {geo}_{\mathrm{ext}} = \int _I \left\langle {\bar{\mathbf {x}}_{\mathbf {u}\mathbf {u}}} (\xi ),{\mathbf f}(\xi )\right\rangle {{\mathrm{d \!}}}\xi \). For a concentrated load \(\mathbf f\) at \(\xi _0\), the external load and its derivative with regard to the degrees of freedom \(\mathbf {u}\) are given by \(\mathbf p_{\mathrm{ext}} = \left\langle {\bar{\mathbf {x}}_{\mathbf {u}}} ({\xi _0}),\mathbf f\right\rangle \) and \(\mathbf k^\mathrm {geo}_{\mathrm{ext}} = \left\langle {\bar{\mathbf {x}}_{\mathbf {u}\mathbf {u}}} ({\xi _0}),\mathbf f\right\rangle \). For an external moment \(M_0\), applied at \((\xi _0)\), the external load and its derivative with regard to the degrees of freedom \(\mathbf {u}\) are given by \(\mathbf p_{\mathrm{ext}} = -M_0 \left\| \dot{\mathbf {x}}\right\| ^{-2} \left| \dot{\mathbf {x}},\dot{\mathbf {x}}_{\mathbf {u}}\right| \) and \(\mathbf k^\mathrm {geo}_{\mathrm{ext}} = M_0 \left( {2 \left\| \dot{\mathbf {x}}\right\| ^{-4} \left\langle \dot{\mathbf {x}}_{\bullet \mathbf {u}},\dot{\mathbf {x}}\right\rangle \left| \dot{\mathbf {x}},\dot{\mathbf {x}}_{\mathbf {u}}\right| - \left\| \dot{\mathbf {x}}\right\| ^{-2} \left( {\left| \dot{\mathbf {x}}_{\bullet \mathbf {u}},\dot{\mathbf {x}}_{\mathbf {u}}\right| + \left| \dot{\mathbf {x}},\dot{\mathbf {x}}_{\mathbf {u}\mathbf {u}}\right| }\right) }\right) \), evaluated at \(\xi _0\).

Gauss quadrature

The integral of the local energy of the beam, as well as its derivatives, over the domain \(\Omega \) is being approximated by Gauss quadrature. The minimum number of Gauss points, \(n_{\mathrm {Gauss}}\), that are required to achieve an exact integration, in the limit of \({\bar{\mathbf {Z}}}\rightarrow \mathbf 0\), results from the polynomial degree of the function space generated by the shape functions, denoted d, which in turn depends on the values \(\bar{q}_I\), \(I\in \{1,2\}\), of the nodes associated with the beam element. For the meaning of \(\bar{q}_I\), see “Full p-refinement and selective p-refinement” section in Appendix 1. With \(d = \bar{q}_1 + \bar{q}_2 + 2\), \(2n_{\mathrm {Gauss}}- 1 \ge 2({d-1})\), we obtain \(n_{\mathrm {Gauss}}\ge \frac{1}{2} ({2({d-1}) + 1}) = \bar{q}_1 + \bar{q}_2 + \frac{3}{2}\).

Thus, for a beam element without Dirichlet boundary conditions, based on one shape function per nodal degree of freedom, i.e. \(\bar{q}_1 = \bar{q}_2 = 0\), we obtain \(n_{\mathrm {Gauss}}\ge 2\), while for a beam element based on two shape functions per nodal degree of freedom, we have \(\bar{q}_1 = \bar{q}_2 = 1\) and \(n_{\mathrm {Gauss}}\ge 4\). For large deformations, the interaction of the shape functions due to the exponentiation increases and may necessitate a higher number of Gauss points in order to achieve a good approximation of the total energy and its derivatives with regard to the nodal degrees of freedom.

The spatial derivatives of the current configuration

In order to obtain the derivatives of the current configuration with regard to the degrees of freedom, \(\dot{\mathbf {x}}\) and \(\ddot{\mathbf {x}}\), we apply a calculation that can be subdivided into two steps.

In the first step, we set up formulas for the derivatives of \(\dot{\mathbf {x}}\) and \(\ddot{\mathbf {x}}\) with respect to the degrees of freedom \(\mathbf {u}\). Those formulas contain the derivatives of \(\mathrm{{e}}^{{\bar{\mathbf {Z}}}( {\xi ,\mathbf {u}})}\) with regard to the parametrization variable \(\xi \) and the degrees of freedom.

In the second step, we derive the expressions for the derivatives of the exponential \(\mathrm{{e}}^{{\bar{\mathbf {Z}}}({\xi ,\mathbf {u}})}\). Instead of the well-known analytical form for the derivative of the matrix exponential given in Eq. (16), we use the definition of the exponential as a power series (see Eq. (15)), and obtain the derivatives of that power series by taking the derivatives of each summand, resulting in an expression of \(\partial _\xi \mathrm{{e}}^{{\bar{\mathbf {Z}}}( {\xi ,\mathbf {u}})}\) as a power series. In addition to simplifying the calculation of the derivatives, the results do not contain any expressions involving the inverse of a matrix, so that they are defined at the origin and numerically stable in the vicinity of the origin. The latter property is of particular importance, as the origin in the vector space of the coefficients related to the degrees of freedom corresponds to the identity map, which is often the preferred choice of an initial estimate for a non-linear approximation algorithm. We note that there exist a number of other series expansions approximating the directional derivative of the matrix exponential that may offer increased computational efficiency [18].

In this context, we note that any closed-form expression of the exponential of a non-trivial semidirect product of two Lie algebras \(\mathfrak h \ltimes \mathfrak n\), including for the Lie algebra associated with the special Euclidean group, will include a rational expression that is not defined at the origin of \(\mathfrak h\) and induces numerical instabilities in its vicinity, and that the associated numerical difficulties increase with the order of differentiation of the exponential map as given by such a closed-form expression.

A good approximation of the power series involved in the computation, with a relatively sharp estimate of the approximation error, can be obtained by calculating the sum of a small number of its leading summands. Thus, it is possible to significantly reduce the computational effort involved in the calculation of the derivatives. The calculation of the upper bound of the approximation error is described in Appendix 2. The calculation of the derivative is the crucial step to determine the relevant derivatives of the current configuration \(\mathbf {x}\left( \xi ,\mathbf {u}\right) = {e}^{{\bar{\mathbf {Z}}}\left( \xi ,\mathbf {u}\right) } {\mathbf {x}}_0\left( \xi \right) \), which is otherwise straightforward. We also note that small proportional errors in the calculation of the derivatives do not affect the order of convergence of the Newton-Raphson root-finding algorithm [1, p. 155].

We define the complex numbers \(\zeta _1\), \(\zeta _2\) such that \(\bigg (\begin{array}{ll} \zeta _1 &{} \zeta _2 \\ 0 &{} 0\end{array}\bigg ) \equiv {\bar{\mathbf {Z}}}\). Using this notation, the matrix exponential of \({\bar{\mathbf {Z}}}\) is given by

$$\begin{aligned} \mathrm{{e}}^{\bar{\mathbf {Z}}}= & {} \sum _{k=0}^\infty \frac{1}{k!} \begin{pmatrix} \zeta _1 &{} \quad \zeta _2 \\ 0 &{} 0\end{pmatrix}^k = \left\{ \begin{array}{llll} &{}\mathbf {I}+ \sum \nolimits _{k=1}^\infty \frac{{\zeta _1}^k}{k!} \begin{pmatrix} 1 &{}\quad \frac{\zeta _2}{\zeta _1} \\ 0 &{} 0\end{pmatrix} &{}\quad &{} \text {if }\quad \zeta _1 \ne 0\\ &{}\mathbf {I}+ \begin{pmatrix} 0 &{}\quad \zeta _2 \\ 0 &{} \quad 0 \end{pmatrix}&\quad&\text {if }\quad \zeta _1 = 0 \end{array}\right\} \nonumber \\= & {} \begin{pmatrix} \mathrm{{e}}^{\zeta _1} &{} \quad \sum \nolimits _{k=0}^\infty \frac{{\zeta _1}^k}{({k+1})!} \zeta _2 \\ 0 &{}\quad 1 \end{pmatrix}. \end{aligned}$$
(49)

We obtain the first and the second derivative of \(\mathrm{{e}}^{\bar{\mathbf {Z}}}\) with regard to the parameterization variable \(\xi \) as

$$\begin{aligned} \partial _\xi \mathrm{{e}}^{\bar{\mathbf {Z}}}&= \begin{pmatrix} \mathrm{{e}}^{\zeta _1} \zeta _{1,\xi }&{}\quad \sum _{k=0}^\infty \frac{{\zeta _1}^k}{({k+2})!} \left( {({k+1})\zeta _2 \zeta _{1,\xi }+ ({k+2}) \zeta _{2,\xi }}\right) \\ 0 &{}\quad 0 \end{pmatrix}, \end{aligned}$$
(50a)
$$\begin{aligned} \partial ^2_\xi \mathrm{{e}}^{\bar{\mathbf {Z}}}&= \begin{pmatrix} \mathrm{{e}}^{\zeta _1} \left( {\zeta _{1,\xi }+ \zeta _{1,\xi \xi }}\right) &{} \quad \sum _{k=0}^\infty \frac{{\zeta _1}^k}{({k+3})!} \left( { \begin{array}{ll} &{} ({k+1}) ({k+2}) \zeta _2 {\zeta _{1,\xi }}^2 + \\ &{} ({k+1}) ({k+3}) \left( {2\zeta _{1,\xi }\zeta _{2,\xi }+ \zeta _2 \zeta _{1,\xi \xi }}\right) + \\ &{} ({k+2}) ({k+3}) \zeta _{2,\xi \xi } \end{array}}\right) \\ 0 &{}\quad 0 \end{pmatrix}. \end{aligned}$$
(50b)

These derivatives, as well as all other derivatives of \({\bar{\mathbf {Z}}}\) with respect to the parameterization variable \(\xi \) and the degrees of freedom \(\mathbf {u}\), can be written in the form of matrices with non-zeros entries of the form \(\sum _{k=0}^\infty \frac{1}{({k+m})!} {\zeta _1}^k \sum _{n=1}^N \prod _{j\in J_n}({k+j}) {\delta _n} (\xi )\), with \(m \in {\mathbb {N}}{}\), \(N \in {\mathbb {N}}{}\), \(J_n \subset {\mathbb {N}}{}\). In this equation, \( {\delta _n} (\xi )\) denotes a linear combination of products of derivatives of \(\zeta _1\) and \(\zeta _2\). In particular, \( {\delta _n} (\xi )\), which is the most computationally expensive term, does not depend on k. The calculation of the derivatives, although leading to several quite lengthy expressions, is straightforward and thus not included in this outline of the model.

Results and discussion

In order to evaluate the characteristics and the performance of the LogFE method, we have implemented a formulation for a Bernoulli beam in the plane as a program written in the numerical computing environment MatLab. This section presents the results of a number of numerical examples characterized by different external loads and boundary conditions. The prismatic beam has length \(l = 1\, \hbox {m}\) and a constant, rectangular cross-section given by the height \(h ={0.08}\,\hbox {m}\) and the width \(b ={0.1}\,\hbox {m}\). Young’s modulus is set to \(E = {3.4 \cdot 10^9}\,{\hbox {Nm}^{-2}}\). This value is consistent with the material properties of a hardened epoxy resin. As in the previous sections, \(\xi \) denotes the parameterization of the neutral axis of the beam, running from \(\xi = 0\) on the left side to \(\xi = 1\) on the right side of the beam. The deflection of the beam is given by \(w := - \mathrm{Im}({\mathbf {x}- {\mathbf {x}}_0})\,\hbox {m}\), and the axial displacement is given by \(u := \mathrm{Re}({\mathbf {x}- {\mathbf {x}}_0})\,\hbox {m}\). Both quantities are expressed in the unit metre (m). (See Fig. 2 for the identification of the complex plane with the 2-dimensional Euclidean space.) Note that the stress resultants are given as normalized values, i.e. as relative values with respect to the standard load value given for the load case.

In the presence of boundary conditions, selective \( p \)-refinement, as described in “Full p-refinement and selective p-refinement” in Appendix 1, has been applied. At simply supported nodes, the hinge condition modification given in “Hinge condition” in Appendix 1 has been used. The shape functions that have been used in the numerical examples, shown in Table 7, can be obtained from Table 8.

The focus of the following presentation and discussion of numerical examples is on characterizing the approximation properties induced by the discretization of the function space of the deformation function. Therefore, in order to eliminate, for all practical purposes, the effects of the approximation of the exponential, and of the approximate evaluation of the current configuration, the number of summands in the truncated power series approximating the exponential has been set to sixteen, and eighteen linearly spaced evaluation points have been used to evaluate the current configuration. It is obvious that, in a production-oriented implementation, these parameters should be adjusted to reduce computational cost.

Finally, the numerical algorithm can be set up with the element matrices and vectors described in “Quasi-static equilibrium” section. These element matrices are based on the geometrically exact evaluation of the strain and the curvature at the evaluation points. The application of the Newton-Raphson algorithm follows standard procedures, and the element-related data enter into the iterative calculation in the usual way.

As reference, we use the solutions given by a discretization of the beam by 96 non-linear, finite rotation beam elements endowed with linear shape functions and Reissner kinematics [21]. The formulation of this element is derived from the shell element presented in [32]. The simulation was realized in the finite element analysis program FEAP [29]. As a finite value for the shear module was required by the specific implementation of the element in the program, a very high value of \(5 \cdot 10^{11}\,{\mathrm{Nm}^{-2}}\) was chosen, thus virtually eliminating shear strains. In addition, several solutions based on the Reissner formulation and a small number of elements have been obtained from simulations performed in FEAP.

For comparison, we include the results of a non-linear simulation in Abaqus/Standard [10] using the Euler-Bernoulli beam element B23, which is endowed with a linear interpolation function for the axial and a cubic interpolation function for the transversal displacements. Due to the strong nonlinearity of the observed deformations, the simulations based on the B23 beam element are based on a large number of load steps, applying an exponential increase in the load factor. For some load cases, only models using a small number of elements succeed, while higher numbers of elements lead to a breakdown of the convergence. Therefore, we report the results for the 12-element beam model, where available. If the 12-element model fails, we report the results of the model with the highest number of elements that resulted in a converged solution. According to the Abaqus manual, the B23 beam element should be employed for “small-strain, large-rotation” analysis. A loss of accuracy for the load cases characterized by a strong load intensity is therefore to be expected. For the other load cases, the displacements and bending moments obtained by Abaqus are generally in good agreement with the results obtained by FEAP. For the normal force, larger differences between the solutions obtained by Abaqus and by FEAP occur. However, it should be noted that in most cases, the normal forces contribute only a relatively small share of the total normal stress of the beam that can be calculated from the bending moment and normal force along its neutral axis.

The number of elements, nodes, and degrees of freedom for a simply supported beam (pin/pin) and a beam with mixed boundary conditions (pin/clamped) are reported in Table 5. The LogFE formulation, in general, includes both nodal and internal degrees of freedom. Nodal degrees of freedom may be linked to the nodal degrees of freedom of other finite element elements, e.g. by continuity conditions. In contrast, internal degrees of freedom are not linked to any degrees of freedom of other elements, and thus can be eliminated from the model by static condensation.

When interpreting the results reported below, one should bear in mind that the primary goal of the LogFE method is to approximate the low-frequency component of a deformation resulting from a given load condition with a small number of degrees of freedom. It is thus to be expected that the results display some deviations from the reference solution, especially with regard to high-frequency components of the deformation. This is especially evident in load cases involving concentrated forces. We also stress that the calculation of the estimates for both the position and the stress resultants for any material point of the beam is straightforward in the case of the LogFE formulation, while similar calculations based on the conventional FE formulation require the identification of appropriate interpolation functions.

Table 5 Number of elements, nodes, and degrees of freedom of the LogFE and the conventional finite element formulation

For specific load cases, the polynomials presented in “Consistency with the linear beam theory” section result, for small deformations, in a solution of the model that is consistent with the linear deflection curve of the Bernoulli beam theory. Note that for small deformations, the degrees of freedom, which are defined on the logarithmic space, are located close to the origin of the vector space of the Lie group, which represents the application of the identity function to the initial configuration. For other load cases, the bases provided by the shape functions do not allow for a solution consistent with the linear theory in a neighborhood of the identity function. With increasing load, the components of the exact solution that do not lie within the function spaces generated by the shape functions generally increase, and, as a result, the quality of the finite element solution diminishes.

Similar to the \( hp \)-FEM concept, better solutions for higher load intensities can be obtained by increasing the number of finite elements or by increasing the number of shape functions associated with a single finite element, e.g. by adding shape functions based on polynomials of higher order. For many load cases, increasing the number of elements will be an efficient method to increase the accuracy of the results. However, a meaningful study of the advantages and disadvantages of increasing the number of elements versus an increase in the number of shape function would have to include translational degrees of freedom at each node, in addition to the rotational and dilatational degrees of freedom that have been incorporated into the finite element formulation presented in this paper. As the inclusion of nodal translations leads to a significant extension of the method at the conceptional level, we will not consider the possibility of increasing the number of elements at this time, and instead focus on the effects of an increase of the number of shape functions.

Fig. 6
figure 6

Load cases with different boundary conditions

External moments and distributed loads

External moments and distributed loads generally result in deformations that are characterized by a large relative low-frequency component for the deformation as well as for the stress resultants. While fairly good approximations of the deformation and the stress resultants can already be obtained with a single shape function per nodal degree of freedom, the results of the application of two shape functions, i.e. two polynomial functions per nodal degree of freedom, are often quite close to the exact solution, not only with regard to the deformed configuration, but also with regard to its derivatives and thus the stress resultants, which are derived from those derivatives. In order to illustrate some of the characteristics of the LogFE method, we present the results of several load cases, given in Fig. 6.

Table 6 Strain energy (\(U_{\mathrm {strain}}\)) and bending energy (\(U_{\mathrm {bending}}\)) for different load cases and load intensities
Table 7 Shape functions employed in the numerical examples

Due to the given boundary conditions, which preclude any dislocation of the nodes, the results for all load cases show a rather high share of the strain energy in the total internal energy of the beam, leading to strong non-linear deformations. The absolute values of the strain energy and the bending energy, as well as the relative share of the strain energy, are reported in Table 6.

Fig. 7
figure 7

Comparison of LogFE and conventional finite element formulations (Load case A). FEAP: Reissner kinematics, Abaqus: Cubic interpolation. 1p./2p.: One/two polynomial(s) per basis on the Lie algebra

Fig. 8
figure 8

Comparison of LogFE and conventional finite element formulations (Load case B). FEAP: Reissner kinematics, Abaqus: Cubic interpolation. 1p./2p.: One/two polynomial(s) per basis on the Lie algebra.

For load case A, a simply supported beam subjected to a moment on the left support (see Fig. 6a), Fig. 7 compares the results of the LogFE formulation with a 96-element reference solution as well as with the results of a 3-element and a 6-element finite element formulation of the same beam. For the moderate load (see Fig.  7a), the formulation results in fairly accurate solutions both for one shape function (no internal degrees of freedom) and for two shape functions (one internal degree of freedom) per nodal degree of freedom. In particular, the axial displacement, which is rather poorly approximated by the 3-element discretization, is fairly well captured by the results of the LogFE formulation.

The estimate for the normal force, interpreted in isolation, shows a significant deviation of the LogFE solutions with regard to the reference solution. However, for moderate load intensities, the normal force represents a rather small share of the total normal stress of the beam. As the respective graph in Fig. 7 shows, the impact of the error of the estimation of the normal force is fairly small, while the error in the estimation of the bending moment dominates the overall error in the normal stress at the upper and lower side of the beam. For stronger loads (see Fig. 7b), the interplay of the bending energy and the strain energy results in a general deterioration of the estimate of the normal force. In the load case presented in Fig. 7, the estimate of the normal force close to the left support of the beam becomes particularly poor, especially in the case of two shape functions per nodal degree of freedom. However, the graph of the normal stress shows that the normal force contributes only a small share of the total internal energy of the beam, as well as of the local contributions of the cross-sections along the neutral axis to the total internal energy.

Figure 8 displays the results of load case B, the application of a line load across the entire neutral axis of the beam. The beam is simply supported on the left side and clamped on the right side (see Fig. 6b). As in load case A, for moderate load intensities (see Fig. 8a), the LogFE formulation achieves a good approximation of the actual deformation for moderate load intensities. For higher load intensities (see Fig. 8b), the solution given by one shape function per nodal degree of freedom deteriorates, while the solution using two shape functions remains quite close to the reference solution. For this load case, the local relative approximation error remains stable as the load intensity increases. It is noteworthy that the rather complex axial displacement of the beam under the given load conditions is fairly well approximated by the LogFE solution, given the small number of degrees of freedom employed by the formulation.

Fig. 9
figure 9

Comparison of LogFE and conventional finite element formulations (Load case C). FEAP: Reissner kinematics, Abaqus: Cubic interpolation. 1p./2p.: One/two polynomial(s) per basis on the Lie algebra

Fig. 10
figure 10

Comparison of LogFE and conventional finite element formulations (Load case D). FEAP: Reissner kinematics, Abaqus: Cubic interpolation. 1p./2p.: One/two polynomial(s) per basis on the Lie algebra

Other load cases yield similar results. For load case C (see Fig. 6c), a beam that is simply supported at the left side, clamped at the right side, and subjected to an external moment at the left support, the solution based on one shape function per nodal degree of freedom, for higher load intensities, does not result in a valid approximation of the actual stress resultants. However, the solution employing two shape functions yields results that are very close to the reference solution (see Fig. 9b). On the other hand, for lower load intensities, the reference solution, as well as the approximations, remain quite close to the linear theory, which would predict a normalized bending moment given as \({M(\xi )}/{M_0} = -\frac{3}{2} \xi + 1\) (see Fig. 9a).

Figure 10 shows the forces and moments resulting from a line load that varies along the neutral axis of the beam (Load case D, see Fig. 6d). The LogFE solution is close to the reference solution for both the moderate and the high load intensity case. The load case also shows that the LogFE formulation can provide good approximations for load cases that induce deformations which fall outside of the deformation space generated by the shape functions in the linear limit.

Load case E (see Figs. 6e, 11) illustrates the frequency-dependence of the accuracy of the solution provided by the LogFE formulation. For the solutions based on one and on two shape functions per nodal degree of freedom, the estimates of the bending moment and of the normal force captures overall shape of the reference solution. The error consists mainly of higher frequency components, changing signs several times on the parameter interval along the neutral axis of the beam.

Load case F (see Fig. 6f) involves a clamped beam subjected to a triangular load. Due to the boundary conditions, only shape functions of higher polynomial degree are being used to model the rotations. Figure 12 shows that the LogFE formulation quite accurately predicts the deformation of the beam for both the moderate and the strong load intensity.

Fig. 11
figure 11

Comparison of LogFE and conventional finite element formulations (Load case E). FEAP: Reissner kinematics, Abaqus: Cubic interpolation. 1p./2p.: One/two polynomial(s) per basis on the Lie algebra

Fig. 12
figure 12

Comparison of LogFE and conventional finite element formulations (Load case F). FEAP: Reissner kinematics, Abaqus: Cubic interpolation. 1p./2p.: One/two polynomial(s) per basis on the Lie algebra

Fig. 13
figure 13

Load case G: A concentrated load

Fig. 14
figure 14

Comparison of LogFE and conventional finite element formulations (Load case G). FEAP: Reissner kinematics, Abaqus: Cubic interpolation. 1p./2p.: One/two polynomial(s) per basis on the Lie algebra

Concentrated loads

Load case G (see Figs. 13, 14), includes a concentrated load. It illustrates that the LogFE formulation is ill-suited for loads that induce high-frequency or highly localized changes in the internal forces and moments. The approximation of the solution based on one shape function per nodal degree of freedom (i.e. without internal degrees of freedom) is of poor quality across the entire parameter interval, and the solution based on two shape functions is still characterized by a large error near the location of the concentrated load. In these cases, it will probably be preferable to increase the number of elements and to place an additional node at the location of the concentrated load, rather than to increase the number of shape functions.

Conclusions

In this work, we have presented a novel finite element formulation which allows to formulate polynomial shape function on the space of the Lie algebra associated with the deformation function, i.e. on the logarithmic space. We therefore refer to the proposed approach as the “Logarithmic finite element method”, or “LogFE method”. The deformation function is defined on the Lie group of planar similarity transformations \( {Sim}({2, \mathbb {R}{}}) = ({{GL \,}{(1, \mathbb {R}{})} \times SO(2)}) \ltimes _\mathrm{Id} T ({2, \mathbb {R}{}})\) and acts on the initial configuration. Associating the polynomial scalar part of the shape functions with vectors that span across the rotational/dilatational and translational subalgebras of \( {sim}{(2, \mathbb {R}{})}\) induces a strong coupling of the translational and rotational variables. The numerical examples show that this coupling conforms well to the interaction of the rotational and the translational components of deformations observed in typical load cases.

Using the example of a two-node beam element endowed with Bernoulli kinematics, the concepts involved in the construction of the LogFE method have been explained in detail. In particular, by choosing appropriate basis functions and polynomial shape function on the Lie algebra \({sim}{(2, \mathbb {R}{})}\), the LogFE method results in element formulations that are, in the case of a beam element, not limited to a constant curvature along the neutral axis within a single element. Therefore, the LogFE method is able to provide a smoother interpolant than other methods that rely on the identification of rotational degrees of freedoms with elements Lie algebras associated with the respective rotational groups. The stress resultants obtained in the numerical experiments indicate that the resulting estimate of the curvature, although not constraint to a constant function within a single finite element, provides a good approximation to the reference solution, especially for those models that include additional polynomial shape functions associated with internal degrees of freedom.

The results obtained from the numerical experiments show that models based on the LogFE method can indeed, with fewer degrees of freedom, achieve approximations of the low frequency component of deformations induced by different external loads that display an accuracy comparable to existing finite element formulations. The numerical experiments also confirmed that the method indeed focuses on the approximation on the low-frequency component of the deformation function. Thus, as a standalone formulation, the LogFE method is ill-suited to solve problems characterized by spatial high-frequency deformations. In these cases, including the modelling of cracks, the use of extended finite element methods (XFEM), i.e. enriching the finite element formulation by high-frequency or discontinuous shape functions, will often be the preferred approach.

In this paper, we have restricted the model to rotational and dilatational degrees of freedom. We intend to present an extended formulation that includes translational degrees of freedom, in addition to rotations and dilatations, in a future publication.

References

  1. Absil PA, Mahony R, Sepulchre R. Optimization algorithms on matrix manifolds. Princeton: Princeton Univ. Press; 2008.

    Book  MATH  Google Scholar 

  2. Argyris J. An excursion into large rotations. Comput Methods Appl Mech Eng. 1982;32:85–155.

    Article  MathSciNet  MATH  Google Scholar 

  3. Armero F, Valverde J. Invariant Hermitian finite elements for thin Kirchhoff rods. I: The linear plane case. Comput Methods Appl Mech Eng. 2012;213–216:427–57.

  4. Armero F, Valverde J. Invariant Hermitian finite elements for thin Kirchhoff rods. II: The linear three-dimensional case. Comput Methods Appl Mech Eng. 2012;213–216:458–85.

    Article  MathSciNet  MATH  Google Scholar 

  5. Ball RS. A treatise on the theory of screws. Cambridge: Cambridge Univ. Press.

  6. Betsch P, Steinmann P. Frame-indifferent beam finite elements based upon the geometrically exact beam theory. Int J Numer Methods Eng. 2002;54(12):1775–88.

    Article  MathSciNet  MATH  Google Scholar 

  7. Borri M, Bottasso CL. An intrinsic beam model based on a helicoidal approximation. Part I: formulation. Int J Numer Methods Eng. 1994;37(13):2267–89.

    Article  MATH  Google Scholar 

  8. Ciarlet PG. The finite element method for elliptic problems, studies in mathematics and its applications. vol. 4. 2nd ed. Amsterdam; 1979.

  9. Cosserat EMP, Cosserat F. Théorie des corps déformables. Paris: A. Hermann et fils; 1909.

    MATH  Google Scholar 

  10. Dassault Systèmes. Abaqus/standard. 2014. http://www.3ds.com/products-services/simulia/portfolio/abaqus/abaqus-portfolio/abaqusstandard/.

  11. Ericksen JL, Truesdell C. Exact theory of stress and strain in rods and shells. Arch Ration Mech Anal. 1958;1(4):295–323.

    MathSciNet  MATH  Google Scholar 

  12. Eugster SR, Hesch C, Betsch P, Glocker C. Director-based beam finite elements relying on the geometrically exact beam theory formulated in skew coordinates. Int J Numer Methods Eng. 2014;97(2):111–29.

    Article  MathSciNet  Google Scholar 

  13. Hall BC. Lie groups, lie algebras, and representations: an elementary introduction, graduate texts in mathematics, vol. 222. New York: Springer-Verlag; 2010.

    Google Scholar 

  14. Jelenić G, Crisfield MA. Geometrically exact 3D beam theory: implementation of a strain-invariant element for statics and dynamics. Comput Methods Appl Mech Eng. 1999;171:141–71.

    Article  MathSciNet  MATH  Google Scholar 

  15. Klein F. Zur Schraubentheorie von Sir Robert Ball. Mathematische Annalen. 1906;62(3):419–48.

    Article  MathSciNet  MATH  Google Scholar 

  16. Knapp AW. Lie groups beyond an introduction. Boston: Birkhäuser; 2005.

    MATH  Google Scholar 

  17. Meier C, Popp A, Wall WA. An objective 3D large deformation finite element formulation for geometrically exact curved Kirchhoff rods. Comput Methods Appl Mech Eng. 2014;278:445–78.

    Article  MathSciNet  Google Scholar 

  18. Najfeld I, Havel T. Derivatives of the matrix exponential and their computation. Adv Appl Math. 1995;16(3):321–75.

    Article  MathSciNet  MATH  Google Scholar 

  19. Napov A, Notay Y. Smoothing factor, order of prolongation and actual multigrid convergence. Numerische Mathematik. 2011;118(3):457–83.

    Article  MathSciNet  MATH  Google Scholar 

  20. Prathap G. The curved beam/deep arch/finite ring element revisited. Int J Numer Methods Eng. 1985;21:389–407.

    Article  MATH  Google Scholar 

  21. Reissner E. On one-dimensional finite-strain beam theory: the plane problem. Zeitschrift für angewandte Mathematik und Physik ZAMP. 1972;23(5):795–804.

    Article  MATH  Google Scholar 

  22. Romero I. The interpolation of rotations and its application to finite element models of geometrically exact rods. Comput Mech. 2004;34(2):121–33.

    Article  MathSciNet  MATH  Google Scholar 

  23. Sander O. Geodesic finite elements for cosserat rods. Int J Numer Methods Eng. 2010;82:1645–70.

    MathSciNet  MATH  Google Scholar 

  24. Schröppel C, Wackerfuß J. Low-frequency shape functions on the logarithmic space. In: Presentation held at the \(11^{\rm th}\) World congress on computational mechanics, Barcelona, 20–25 August 2014.

  25. Schröppel C, Wackerfuß J. Polynomial shape functions on the logarithmic space: the LogFE method. PAMM. 2015;15:469–70.

    Article  Google Scholar 

  26. Selig JM, Ding X. A screw theory of static beams. In: Institute of Electrical and Electronics Engineers, Robotics Society of Japan, IROS; 2001, p. 312–317.

  27. Simo JC. A finite strain beam formulation. The three-dimensional dynamic problem. Part I. Comput Methods Appl Mech Eng. 1985;49(1):55–70.

    Article  MathSciNet  MATH  Google Scholar 

  28. Sonneville V, Cardona A, Brüls O. Geometrically exact beam finite element formulated on the special Euclidean group. Comput Methods Appl Mech Eng. 2014;268:451–74.

    Article  MathSciNet  MATH  Google Scholar 

  29. Taylor RL. FEAP. A finite element analysis program: user manual, v8.4. 2013. www.ce.berkeley.edu/projects/feap/manual84

  30. von Mises R. Motorrechnung: Ein neues Hilfsmittel der Mechanik. Zeitschrift für Angewandte Mathematik und Mechanik. 1924;4(2):155–81.

    Article  MATH  Google Scholar 

  31. von Mises R. Motor calculus: a new theoretical device for mechanics. Technische Universität Graz: Institut für Mechanik; 1996.

  32. Wagner W. A finite element model for non-linear shells of revolution with finite rotations. Int J Numer Methods Eng. 1990;29(7):1455–71.

    Article  MATH  Google Scholar 

Download references

Author's contributions

The method, initially envisaged by CS, was developed and formulated in close cooperation between both authors. Both authors have made substantive contributions to the development and implementation of the proposed method. Both authors read and approved the final manuscript.

Acknowledgements

Financial support for this research was provided by the Deutsche Forschungsgemeinschaft (DFG) via the Emmy Noether Program under Grant No. Wa 2516/3-1. This support is gratefully acknowledged.

Competing interests

The authors declare that they have no competing interests.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Jens Wackerfuß.

Appendices

Appendix 1: Boundary conditions

In the context of conventional finite element methods, boundary conditions are set by fixing the respective degrees of freedom and thus eliminating them from the calculation of the quasi-static equilibrium. Because the number of elements located at the constrained boundary is generally small in comparison to the total number of elements, the loss of accuracy resulting from a smaller set of degrees of freedom within the finite elements at the boundary is often negligible with regard to the entire configuration. In the presence of boundary conditions, the accuracy of the model may be increased by different strategies, such as refining the mesh near the boundary (h-refinement) or choosing more complex finite elements near the boundary (p-refinement).

Full p-refinement and selective p-refinement

In full p-refinement, all shape functions of an element located at the boundary are being refined, while in selective p-refinement, only those shape functions associated with degrees of freedom belonging to the node(s) of the element that are located at the constrained boundary are being refined.

In general, the refinement of any shape function in a finite element necessitates some modification of the other shape functions of that element, in order to keep the impact of the different degrees of freedom of the element as independent from each other as possible. In the following, we assume that shape functions for node I are chosen such that \( {\partial ^\mu N_{I,r,q}} ({\xi _I}) = 0 \,\forall \, \{\mu | \mu < q\}\). Then, in order that the shape for the node \(J \ne I\) do not interfere with the local impact of the shape functions for node I in the neighborhood of \(\xi _I\), the respective derivatives must vanish. If \(\bar{q}_I\) is the highest q in the set of shape functions associated with node I and Lie algebra base \(r\), then all derivatives of the shape functions associated with node J and Lie algebra base \(r\) up to the order \(\bar{q}_I + 1\) must vanish at position \(\xi _I\). Thus, a change in \(\bar{q}_I\) due to a selective refinement at node \(I\) will, in general, entail a modification of the shape functions \(N_{J,r,q}\) for \(J\ne I\).

Hinge condition

In the case of a simply supported node, the value of the degree of freedom related to the derivative of the rotation can be obtained solely from the local characteristics of the configuration and the load case. At a simply supported node, the curvature of the neutral axis is a function of the external moment, with the proportionality factor given by the flexural rigidity. In the absence of an external moment, the curvature will vanish, implying that the value of the degree of freedom related to the derivative of the rotation will vanish likewise. This, in turn, implies that the derivative of the dilatation will vanish as well. As a result, the values of the associated degrees of freedom, \(N_{I,1,1}\) and \(N_{I,2,1}\), may be set to zero. Such a situation is referred to as a hinge condition.

We note, however, that the value of the respective degree of freedom calculated in this way corresponds to the unconstrained solution in the exact (i.e. non-discretized) solution space. In a non-linear setting, it generally differs from the estimated value resulting in the unconstrained solution in the discretized space of deformation functions. Yet, although applying the hinge condition will not, in general, result in a local energy minimum in the discretized function space with regard to a variation in that degree of freedom, it will generally result in a good approximation with using fewer degrees of freedom.

Table 8 Polynomial functions for different boundary conditions and up to two degrees of freedom per Lie algebra element \({\mathbf {Z}}_{I,r}\)

Boundary-dependent shape functions

Table 8 lists possible polynomial functions \(p_{I,r,q}\) for different boundary conditions and assumptions. It contains the polynomial functions of lowest degree that fulfill the conditions given in Eq. () and, if applicable, Eq. () in “Continuity conditions” section and satisfy \(\partial ^q p_{I,r,q} = 1\). Polynomials of higher degree, and other functions, which are square integrable on the interval \(\left[ 0,1\right] \), whose derivatives \(\partial ^\mu p_{I,r,q}\), \(\mu \le q\), are locally Lipschitz continuous at \(\xi _I\) and whose derivatives \(\partial ^\mu p_{I,r,q}\), \(\mu \le \bar{q}_J + 1\), are locally Lipschitz continuous at \(\xi _J\), are also admissible, if they satisfy the conditions given in Eq. () and, if applicable, Eq. (), as well as the additional condition with regard to the shape functions associated with node J given in “Appendix: Full p-refinement and selective p-refinement”.

Figure 15 presents the scalar-valued shape functions \( {N_{I,r,q}} ({\xi })\) for the case of a beam with a simply supported node on the left side and a clamped node on the right side.

Appendix 2: Error bounds for the derivatives of the matrix exponential

The calculation of various derivatives of matrix exponentials is an essential step in the implementation of the LogFE method. In order to construct a reliable algorithm based on approximations of these derivatives, the availability of an upper limit for the absolute value of the error with regard to the exact solution is crucial. In the following, we derive a sharp and computationally inexpensive upper bound of the absolute value of the error induced by truncating the power series that converge to the exact solution in the analytical limit.

The necessity to truncate, as part of a numerical implementation, the power series in Eq. (15) after a number of summands induces an error with regard to the exact value of the derivative. Naturally, the error depends on the number of summands in the truncated series, but it also depends on the value assumed by the indeterminate. While the exact value of the error can only be computed if the exact value of the derivative is known, it is possible to compute an upper bound for the absolute error using only information obtained in the calculation of the truncated series. The upper bound on the absolute error, which decreases as the number of summands increases, allows for the calculation of the respective power series with an error bounded by an arbitrarily small upper absolute value.

Fig. 15
figure 15

Polynomials of the shape functions for a two-node beam element simply supported (without load) at node 1 and clamped at node 2. The shape functions \({N_{1,1,1}} \left( \xi \right) \), \({N_{1,2,1}} \left( \xi \right) \) and \({N_{2,2,0}} \left( \xi \right) \) are not included in the element formulation due to the given boundary conditions. “rhs”: graph refers to the scale on the right hand side

In order to obtain the upper bound for this error, we consider power series of the form \(\lim _{\nu \rightarrow \infty } s_\nu \) with

$$\begin{aligned} s_\nu = \sum _{k=0}^\nu a_k = \sum _{k=0}^\nu \frac{\prod _{j\in J}({k+j})}{({k+m})!} z^k \end{aligned}$$
(51)

\(z \in \mathbb {C}{}\), \(J \subset {\mathbb {N}}{}\), \(\nu \in {\mathbb {N}}{}\). We obtain an upper limit for the absolute value of the remainder \( {r_{n,m,J}} (z) = \vert s_n - \lim _{\nu \rightarrow \infty } s_\nu \vert \) as follows:

$$\begin{aligned} \frac{\vert a_{\mu +1} \vert }{\vert a_\mu \vert }= & {} \frac{\frac{\prod _{j\in J}({\mu +1+j})}{({\mu +1+m})!} \vert z \vert ^{\mu +1}}{\frac{\prod _{j\in J}({\mu +j})}{({\mu +m})!} \vert z \vert ^\mu } = \frac{1}{\mu + m + 1} \prod \nolimits _{j\in J} \frac{\mu + j + 1}{\mu + j} \vert z \vert \nonumber \\= & {} \frac{1}{\mu + m + 1} \prod \nolimits _{j\in J} \left( {1 + \frac{1}{\mu + j}}\right) \vert z \vert \end{aligned}$$
(52)

Thus, for \(\mu \ge n\), we have \(\frac{\vert a_{\mu +1} \vert }{\vert a_\mu \vert } \le \frac{1}{n + m + 1} \prod \nolimits _{j\in J} \left( {1 + \frac{1}{n + j}}\right) \vert z \vert =:q\) and \(\vert a_\nu \vert \le q^{\nu - n} \vert a_n \vert \). Therefore, we obtain an upper bound \({\hat{r}_{n,m,J}} (z)\) for the absolute remainder as

$$\begin{aligned} {r_{n,m,J}} (z)\le & {} \sum _{k = n+1}^\infty \vert a_k \vert \le \sum _{k = n+1}^\infty q^{k - n} \vert a_n \vert = \vert a_n \vert \sum _{k = 1}^\infty q^k = \vert a_n \vert q \sum _{k = 0}^\infty q^k \mathop {=}\limits ^{q < 1} \vert a_n \vert \frac{q}{1-q}\nonumber \\=: & {} {\hat{r}_{n,m,J}} (z) \end{aligned}$$
(53)
Table 9 Absolute error bound \(\hat{r}_n\) of the estimate for the value of the power series
Fig. 16
figure 16

Approximation of a power series derived from the exponential. The power series \(\sum _{k=0}^\infty \frac{z^k}{\left( k+1\right) !}\) is being approximated by \(\sum _{k=0}^n \frac{z^k}{\left( k+1\right) !}\) on the unit disc. For \(z \in \mathbb C \setminus \{0\}\), the power series \(\sum _{k=0}^\infty \frac{z^k}{\left( k+1\right) !}\) is equivalent to the term \(\frac{{e}^z - 1}{z}\)

Let \( {\Delta _n} ({\bullet })\) denote the absolute error of a power series for n summands and \(\left[ \bullet \right] _{1,2}\) the upper right entry of a \( 2\times 2\) matrix. Then, using Eqs. 50a and 50b, we obtain

$$\begin{aligned} \left[ {\Delta _n}\left( {\partial _\xi {e}^{\bar{\mathbf {Z}}}}\right) \right] _{1,2}&= \vert \zeta _2 \zeta _{1,\xi } \vert \hat{r}_{n,2,\{1\}} ({\zeta _1}) + \vert \zeta _{2,\xi } \vert \hat{r}_{n,2,\{2\}} ({\zeta _1}), \end{aligned}$$
(54)
$$\begin{aligned} \left[ {\Delta _n}\left( {\partial ^2_\xi {e}^{\bar{\mathbf {Z}}}}\right) \right] _{1,2}&= \vert \zeta _2 {\zeta _{1,\xi }}^2 \vert {\hat{r}_{n,3,\{1,2\}}} ({\zeta _1}) + \vert 2\zeta _{1,\xi }\zeta _{2,\xi }+ \zeta _2 \zeta _{1,\xi \xi } \vert {\hat{r}_{n,3,\{1,3\}}} ({\zeta _1}) \nonumber \\&\quad \,+ \vert \zeta _{2,\xi \xi } \vert {\hat{r}_{n,3,\{2,3\}}} ({\zeta _1}). \end{aligned}$$
(55)

The entry in the upper left corner of \(\mathrm{{e}}^{\bar{\mathbf {Z}}}\) and its derivatives depends only on the exponential of \(\zeta _1\) and is thus assumed to be given exactly. Let \( {\Delta _n} \left( {\partial ^\mu _\xi \mathbf {x}}\right) \) denote the absolute error bound of the \(\mu \)th derivative of the current configuration with regard to its parameterization for n summands. Then, for the given parameterization, we have \({\Delta _n} \left( {\partial _\xi \mathbf {x}}\right) = {\Delta _n} \left( {\partial _\xi \mathrm{{e}}^{\bar{\mathbf {Z}}}}\right) \begin{pmatrix} \xi \\ 1 \end{pmatrix} = \left[ {{\Delta _n} \left( {\partial _\xi \mathrm{{e}}^{\bar{\mathbf {Z}}}}\right) }\right] _{1,2}\) and \({\Delta _n} \left( {\partial ^2_\xi \mathbf {x}}\right) = {\Delta _n} \left( {\partial ^2_\xi \mathrm{{e}}^{\bar{\mathbf {Z}}}}\right) \begin{pmatrix} \xi \\ 1 \end{pmatrix} = \left[ {{\Delta _n} \left( {\partial ^2_\xi {e}^{\bar{\mathbf {Z}}}}\right) }\right] _{1,2}\).

Table 9 shows the error bounds \(\hat{r}_n\) for two different power series. For complex numbers z of small magnitude, the relative error becomes negligible after only a few summands, while for complex numbers of larger magnitude—typically associated with large displacement gradients—the error diminishes much slower as the number of summands rises. Figure 16 illustrates the dependence of the error on the number of summands and on the magnitude of the complex number z. It also shows that, while the actual value of the error changes as z varies along a circle given by a constant \(\vert z \vert \), the absolute error remains within a small range that depends on the magnitude of z, making it possible to obtain a relatively sharp upper bound on the magnitude of the error as a function of \(\vert z \vert \). While about three summands are sufficient to obtain very good approximations for moderate deformations, about five summands are necessary to achieve a precision that allows for the phase of quadratic convergence in the Newton-Raphson algorithm to extend up to the point at which commonly used tolerance levels in finite element calculations are being met.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Schröppel, C., Wackerfuß, J. Introducing the Logarithmic finite element method: a geometrically exact planar Bernoulli beam element. Adv. Model. and Simul. in Eng. Sci. 3, 27 (2016). https://doi.org/10.1186/s40323-016-0074-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-016-0074-8

Keywords