1 Introduction

Flexible slender structures like cables, ropes and hoses have a variety of applications in engineering systems. These are, to name a few, textile manufacturing, multi-wire cable structures for the automotive industry or endoscopes for medical applicationsFootnote 1 [51]. Often the behavior of such systems is determined by contact interactions among those slender components. These flexible components can be modeled using a finite element approach with \(1D\) elements based on the non linear theory of beams. For very slender beams, Kirchhoff-type beam models disregarding transverse shearing of the cross section w.r.t. the centerline tangent [5, 6, 19, 20], or additionally the extensibility of the centerline [22, 23], are considered as suitable. A more general class of Tymoshenko type beam models able to capture shear deformation in a simplified manner, describe the kinematics in terms of translations and orientations of reference frames attached to each beam cross-section [7, 8]. This type of formulation involves the handling of finite rotations, which belong to a non-linear space, \(SO(3)\). Often they are treated as decoupled from the translation field [2427], so that the models are formulated on the Lie group \(\mathbb{R}^{3}\times SO(3)\). Alternatively, beam models formulated on the special Euclidean group \(SE(3)\) have been developed [1, 4, 28, 29]. In this formalism rotation and translation variables are intrinsically coupled [1]. Moreover, deformations, arbitrary motion increments as well as internal and external forces are expressed in the local frame attached to the cross section. As a consequence the equilibrium equations and the tangent stiffness matrix are invariant under rigid body motions [3]. This framework is also used in this paper.

Even though \(3D\) continuum contact mechanics has been studied extensively in the literature [31, 32], publications on beam-to-beam contact are relatively scarce. Due to their particular geometric features, beams of circular cross-section may experience contact interactions of different natures. These are: point forces in the range of large contact angles and continuous distributed line forces for small contact angles. The earliest publications about beam-to-beam contact were dedicated to point-to-point contact models. The contact constraint is in general defined from a closest point projection and then enforced via the penalty or the Lagrange multiplier method [33, 34]. In [3739] the case of distributed forces is treated by collocation. A high number of collocation points is then needed to obtain accurate results and the number of kinematic degrees-of-freedom (DOFs) should be adjusted to the density of collocation points to avoid over-constrainment. Indeed, these methods may suffer from convergence issues when the number of contacts is high compared to the number of kinematic DOFs. More recent formulations describe distributed line forces as a continuous field and use a unilateral minimal distance criterion to determine the kinematics of contact [17, 35]. The constraints are enforced via the penalty approach, which yields a smooth representation of the distributed contact force. In [36] both point-to-point and line-to-line contact are covered by switching between the two contact models.

Mortar methods for \(3D\) continuum contact mechanics have been extensively studied in the literature [913, 4046]. These methods are characterized by contact constraints enforced in a weak sense and a saddle point formulation of the problem, where contact pressures play the role of Lagrange multipliers. In this setting, the issues of over-constrained formulations are avoided and optimal convergence rates without the need for smooth contact kinematics may be proven, provided appropriate finite element spaces are chosen. The paper proposes a mortar finite element formulation for frictionless line-to-line beam contact. As a perspective, the ultimate goal is a formulation that treats both point-to-point and line-to-line contact in a consistent manner. The contact model is combined with the \(SE(3)\) local frame formulation for the beam finite element (FE). As a consequence, the constraint gradient and the associated contact forces are expressed in the local frame and their expressions are invariant under rigid body motion of the contact element.

The paper is structured as follows: First the geometrically exact beam model formulated on \(SE(3)\) as developed in [1] is described in Sect. 2. It starts by an introduction to the concept of frames and frame transformations in 2.1, which will be used extensively throughout the paper. The formulation of the weak form and the discretization process follow in Sects. 2.2 and 2.3. Section 3 is dedicated to the proposed contact model. Contact kinematics and invariance properties of the related constraint gradient as they emerge from the local frame approach are discussed in 3.1. The continuous and discretized contact formulations are given in Sects. 3.2 and 3.3, followed by a section about the augmented Lagrangian approach that is employed in this paper to solve the equilibrium equations. In Sect. 3.5 the projection procedure is detailed. Finally, numerical tests are presented in Sect. 3.3. The first one is the traditional patch test. The second test illustrates the behavior of the formulation in the presence of a jump in the distributed contact pressures on a simple \(2D\) example. The last test is the \(3D\) twisting of two beams, which is used to study the properties of the proposed formulation for problems of a higher complexity.

2 Beam model

In this section, the geometrically exact beam model formulated on the special Euclidean group \(SE(3)\) [1] is presented. In this framework translation and rotation variables are naturally coupled through the group operation. This property is conserved after spatial discretization. As a consequence, the finite element is free of any locking effect.

2.1 Frames

Reference frames are ever-present in the kinematic description of multibody systems. Following Sonneville [2], frame transformations are taken as elements of the Lie group \(SE(3)\). The Lie group theory offers a rigorous definition of frame derivatives.

The frame transformation between an inertial frame \(\{O\}\) and a local frame \(\{C\}\) is represented by an element \(\mathbf{H}_{OC}\in SE(3)\). It is composed of a translation \(\mathbf{x}_{OC}\in \mathbb{R}^{3}\) which represents the Cartesian coordinates of the frame center \(C\) in the frame \(\{O\}\) and a rotation \(\mathbf{R}_{OC}\) which belongs to the special orthogonal group \(SO(3)\) and represents the orientation of the base vectors of \(\{C\}\) in the frame \(\{O\}\). It may be written as a 4 by 4 matrix:

$$ \mathbf{H}_{OC}= \begin{bmatrix} \mathbf{R}_{OC} & \mathbf{x}_{OC} \\ \mathbf{0} & 1\end{bmatrix} . $$
(1)

The group operation is simply the matrix product and can be interpreted as a sequence of frame transformations. Variations of the local frame \(\{C\}\) are represented by introducing a left invariant vector field as

$$ \delta \mathbf{H}_{OC}=\mathbf{H}_{OC} \widetilde{\delta \boldsymbol{\pi }}_{C}^{C}, $$
(2)

where \(\delta \boldsymbol{\pi }_{C}^{C}\in \mathbb{R}^{6}\) is interpreted as an arbitrary infinitesimal motion of frame \(\{C\}\) expressed in frame \(\{C\}\). It is invariant under a superimposed frame transformation. It belongs to the Lie algebra of the special Euclidean group, denoted \(\mathfrak{se}(3)\), which is isomorphic to \(\mathbb{R}^{6}\) through the map \(\mathbb{R}^{6}\rightarrow \mathfrak{se}(3):\delta \boldsymbol{\pi }_{C}^{C} \rightarrow \widetilde{\delta \boldsymbol{\pi }}_{C}^{C}\):

$$ \begin{bmatrix} \delta \mathbf{x}_{C}^{C} \\ \delta \boldsymbol{\theta }_{C}^{C}\end{bmatrix} \rightarrow \begin{bmatrix} \widetilde{\delta \boldsymbol{\theta }}_{C}^{C} & \delta \mathbf{x}_{C}^{C} \\ \mathbf{0} & 0\end{bmatrix} , $$
(3)

where the translation part \(\delta \mathbf{x}_{C}^{C}=\mathbf{R}_{OC}^{T}\delta \mathbf{x}_{C}^{O}\) is an element of \(\mathbb{R}^{3}\). The rotation part \(\delta \widetilde{\boldsymbol{\theta }}_{C}^{C}=\mathbf{R}_{OC}^{T} \delta \mathbf{R}_{OC}\) belongs the Lie algebra of \(SO(3)\), denoted \(\mathfrak{so}(3)\). This Lie algebra is the set of 3 by 3 skew-symmetric matrices and is isomorphic to \(\mathbb{R}^{3}\) through the map \(\mathbb{R}^{3}\rightarrow \mathfrak{so}(3):\mathbf{a}\rightarrow \widetilde{\mathbf{a}}\)

$$ \mathbf{a} = \begin{bmatrix} a_{1} \\ a_{2} \\ a_{3} \end{bmatrix} \rightarrow \begin{bmatrix} 0 & -a_{3} & a_{2} \\ a_{3} & 0 & -a_{1} \\ -a_{2} & a_{1} & 0 \end{bmatrix} = \widetilde{\mathbf{a}}. $$
(4)

Whether the \(\widetilde{(\cdot )}\) operator denotes the map in Eq. (3) or in Eq. (4) is generally clear from the context. The variations can be represented in the local frame \(\{C\}\) as \(\delta \mathbf{x}_{C}^{C}\) and \(\delta \boldsymbol{\theta }_{C}^{C}\) or in the inertial frame \(\{O\}\) as \(\delta \mathbf{x}_{C}^{O}\) and \(\delta \boldsymbol{\theta }_{C}^{O}\). In the following, the local frame representation will generally be adopted. For the sake of notational simplicity, the frame superscript will be dropped when referring to local frame quantities.

2.2 Continuous formulation

The choice of a left invariant representation of frame derivatives yields a formulation of static equilibrium equations in the local frame. In other words, strain measures, stress resultants and infinitesimal motions are expressed in the frame attached to the beam centerline. The equations are frame invariant and insensitive to rigid body motions of the beam.

Suppose that the centerline of a beam is parametrized by \(s\in [0,L]\), where \(s\) is the centerline arclength coordinate in the initial undeformed configuration. We attach a frame to each of its points, as schematically represented in Fig. 1. The configuration \(q(s)\) is given by the map \(\mathbb{R}\rightarrow SE(3):s\rightarrow \mathbf{H}_{OC}(s)\), which corresponds to the transformation from the inertial frame \(\{O\}\) to the local frame \(\{C(s)\}\). Strain measures are constructed in accordance with the Lie group representation of the frame derivatives:

$$ \dfrac{\mathrm{d} \mathbf{H}_{OC}(s)}{\mathrm{d} s}=\mathbf{H}_{OC}(s) \tilde{\boldsymbol{\mathit{f}}}_{C}. $$
(5)

The element \(\tilde{\boldsymbol{\mathit{f}}}_{C}(s)\) of the Lie algebra \(\mathfrak{se}(3)\) is interpreted as an objective deformation gradient expressed in the local frame. Then the sectional strains are obtained by the linear relation

$$ \boldsymbol{\epsilon }_{C}=\boldsymbol{\mathit{f}}_{C}- \boldsymbol{\mathit{f}}_{R}, $$
(6)

where \(\boldsymbol{\mathit{f}}_{R}(s)\) is the deformation gradient in the reference configuration. In this paper the reference configuration is identified with the initial undeformed configuration for simplicity. The deformation measures obtained in that way are the same as for classical geometrically exact beam formulations [8]. This procedure can be adapted to other types of structural components such as shells and superelements to define suitable objective deformation measures [2].

Fig. 1
figure 1

Kinematic description of the beam configuration using frame transformations as unknowns

If \(\mathcal{S}\) denotes the space of admissible configurations \(q\) and \(\mathcal{V}\) the space of admissible infinitesimal motions, the weak form of the equilibrium is given by: Find \(q\in \mathcal{S}\) such that

$$ \delta \mathcal{W}^{\text{int}}\left (q,\delta \boldsymbol{\pi }_{C} \right )=\delta \mathcal{W}^{\text{ext}}\left (q,\delta \boldsymbol{\pi }_{C}\right ) \quad \forall \delta \boldsymbol{\pi }_{C}(s) \in \mathcal{V}. $$
(7)

The virtual work done by the external forces is

$$ \delta \mathcal{W}^{\text{ext}}\left (q,\delta \boldsymbol{\pi }_{C} \right )=\int _{0}^{L}(\delta \boldsymbol{\pi }_{C})^{T}\mathbf{f}_{C}^{\text{ext}} \,\mathrm{d}s, $$
(8)

where \(\mathbf{f}_{C}^{\text{ext}}(s)\) contains the resultant forces and moments expressed in the local frame. They may depend on the configuration. The variation of the internal elastic strain energy is given by

$$ \delta \mathcal{W}^{\text{int}}(q,\delta \boldsymbol{\pi }_{C})=\int _{0}^{L}( \delta \boldsymbol{\epsilon }_{C})^{T}\mathbf{K}_{C} \boldsymbol{\epsilon }_{C}\,\mathrm{d}s, $$
(9)

where \(\mathbf{K}_{C}(s)\) is the usual sectional stiffness matrix for a linear elastic material. It will be assumed constant in the remainder of this paper. Equation (5) is a kinematic relation that needs to be verified by the solution together with Eq. (7).

2.3 Spatial discretization

The finite element method is selected to discretize Eqs. (5) and (7), i.e., in order to define the finite dimensional subspaces \(\mathcal{S}_{h}\subset \mathcal{S}\) and \(\mathcal{V}_{h}\subset \mathcal{V}\). For that purpose, consider a two noded beam element of length \(L_{e}\) as shown in Fig. 2. A frame \(\{A\}\) is attached to the node at \(s=0\) and a frame \(\{B\}\) is attached to the node at \(s=L_{e}\). Thus the configuration of the discrete beam element is given by \(q_{e}=\left (\mathbf{H}_{OA},\mathbf{H}_{OB}\right )\). In this work, a geometric interpolation of the two frames is chosen, such that the kinematic relation in (5) is satisfied by construction. To that end the discrete relative configuration vector is computed as

$$ \mathbf{d}_{AB}=\log _{SE(3)}\left (\mathbf{H}_{OA}^{-1}\mathbf{H}_{OB} \right ). $$
(10)

The interpolation formula follows as

$$ \mathbf{H}_{OC}(s)=\mathbf{H}_{OA}\exp _{SE(3)}\left ( \dfrac{s}{L_{e}} \mathbf{d}_{AB}\right ). $$
(11)

It is based on the exponential and logarithm maps, for more details see [1, 3, 4]. The discretized strain is then given by

$$ \boldsymbol{\epsilon }_{C} = \dfrac{\mathbf{d}_{AB}-\mathbf{d}_{AB}^{0}}{L_{e}}, $$
(12)

where \(\mathbf{d}_{AB}^{0}\) is the relative configuration vector in the reference configuration. A few distinctive characteristics of this discretization scheme should be highlighted. First, for the two noded beam element as presented here, \(\boldsymbol{\epsilon }_{C}\) is constant along the span of the element. As a consequence the discretization scheme defines piecewise constant strain elements and helicoidal shape functions [4, 15]. Moreover, as indicated by Eq. (10), \(\mathbf{d}_{AB}\) (and as a consequence \(\boldsymbol{\epsilon }_{C}\)) only depends on the relative configuration between \(\{A\}\) and \(\{B\}\). Therefore, the interpolation scheme is frame invariant. Finally, the interpolation couples translation and rotation variables. Combining it with an appropriate Lie group time integration scheme [16] yields a beam formulation without any global parametrization of rotation and which is inherently locking free [1].

Fig. 2
figure 2

Geodesic interpolation of frame transformations

In general the discretizations of infinitesimal motions and deformation gradients can be written

$$ \delta \boldsymbol{\pi }_{C}(s)=\mathbf{Q}\left (s,\mathbf{d}_{AB} \right )\delta \boldsymbol{\pi }_{e} $$
(13a)
$$ \delta \boldsymbol{\epsilon }_{C}(s)= \mathbf{P}\left (s,\mathbf{d}_{AB} \right )\delta \boldsymbol{\pi }_{e}, $$
(13b)

where \(\delta \boldsymbol{\pi }_{e}= \begin{bmatrix} \delta \boldsymbol{\pi }_{A} \\ \delta \boldsymbol{\pi }_{B} \end{bmatrix} \) collects the nodal infinitesimal motion variations. The methodology may be extended to higher order interpolation [30] and other choices of local parametrization to represent the configuration around the nodal frames [47, 48]. In that case matrices \(\mathbf{Q}\) and \(\mathbf{P}\) must be adapted. Here, since strains are constant elementwise \(\mathbf{P}\) is independent of \(s\).

Then, the discretized weak form of the variational formulation (7) can be obtained from the developed expressions, which leads to the following equations for the element internal and external forces, respectively:

$$ \mathbf{f}_{e}^{\text{int}}=\mathbf{P}^{T}\mathbf{K}_{C} \boldsymbol{\epsilon }_{C} L_{e}, $$
(14a)
$$ \mathbf{f}_{e}^{\text{ext}}=\int _{0}^{L_{e}}\mathbf{Q}^{T}\mathbf{f}_{C}^{ \text{ext}}\,\mathrm{d}s. $$
(14b)

3 Contact model

In this section, a model for frictionless beam contact based on the mortar finite element method is formulated. To keep notations simple, the derivations are made for a two beam system, but they may be generalized to an arbitrary number of beams.

3.1 Contact kinematics and constraint gradient

The state of the system is described by the configuration variable \(q\). Any beam model, with any parametrization could be considered. In this contribution, the beam formulated on the special Euclidean group as presented in Sect. 2 is selected. To each point of the centerline of beam 1 a frame \(\{C(s_{1})\}\) is attached and a frame \(\{F(s_{2})\}\) for beam 2. The configuration of the two beams is then given by two fields of frame transformations as

$$ q(s_{1},s_{2})=\left (\mathbf{H}_{OC}(s_{1}),\mathbf{H}_{OF}(s_{2}) \right ), \qquad \mathbf{H}_{OC},\mathbf{H}_{OF}\in SE(3) ,$$
(15)

and the variations are collected in the 12 dimensional vector

$$ \delta \boldsymbol{\pi }(s_{1},s_{2})= \begin{bmatrix} \delta \boldsymbol{\pi }_{C}(s_{1}) \\ \delta \boldsymbol{\pi }_{F}(s_{2})\end{bmatrix} , \qquad \delta \boldsymbol{\pi }_{C},\delta \boldsymbol{\pi }_{F}\in \mathbf{\mathbb{R}^{6}} .$$
(16)

According to the common nomenclature in contact mechanics, beam 1 is labeled the slave and beam 2 the master. A point on the master centerline \(\bar{\mathbf{x}}_{OF}\) is associated to a point on the slave centerline \(\mathbf{x}_{OC}\) by means of a projection, which will be described in Sect. 3.5. Similarly, rotation matrix \(\bar{\mathbf{R}}_{OF}\) is evaluated at the same location as \(\bar{\mathbf{x}}_{OF}\). Then, a scalar gap function is introduced:

$$ g(q)=\|\bar{\mathbf{x}}_{OF}-\mathbf{x}_{OC}\|-r_{1}-r_{2}, $$
(17)

where \(r_{1}\) and \(r_{2}\) denote the cross section radii of the slave and master beam respectively. It measures the relative position between points on the two beam surfaces under the assumption that the effect of shear deformation on the geometry of the beam’s external skin may be neglected. In case of frictionless contact between slender beams with circular cross-sections, this expression is sufficient to formulate contact conditions, such that the contact model can be entirely written over the centerlines.

By computing the variation of Eq. (17),

$$ \delta g(q) = \begin{bmatrix} \mathbf{G}_{C} & \mathbf{G}_{F} \end{bmatrix} \begin{bmatrix} \delta \boldsymbol{\pi }_{C} \\ \delta \boldsymbol{\pi }_{F} \end{bmatrix} = \mathbf{G}\delta \boldsymbol{\pi }, $$
(18)

the local constraint gradient is obtained as

$$ \mathbf{G}^{T}= \begin{bmatrix} -\mathbf{M}\mathbf{n}^{C} \\ \mathbf{M}\bar{\mathbf{R}}_{OF}^{T}\mathbf{R}_{OC}\mathbf{n}^{C}\end{bmatrix} , $$
(19)

where \(\mathbf{n}^{C}\) denotes the unit outward normal to the slave beam expressed in its local frame. It is defined as:

$$ \mathbf{n}^{C}=\mathbf{R}_{OC}^{T} \dfrac{(\bar{\mathbf{x}}_{OF}-\mathbf{x}_{OC})}{\|(\bar{\mathbf{x}}_{OF}-\mathbf{x}_{OC})\|}= \mathbf{R}_{OC}^{T}\mathbf{n}^{O}. $$
(20)

Note that in Eq. (19) \(\bar{\mathbf{R}}_{OF}^{T}\mathbf{R}_{OC}=\mathbf{R}_{FC}\) is the relative orientation between the local slave and the local master beam. Thus \(\mathbf{n}^{F}=\bar{\mathbf{R}}_{OF}^{T}\mathbf{R}_{OC}\mathbf{n}^{C}\) may be seen as the same normal vector as \(\mathbf{n}^{C}\), but expressed in the local frame of the master. The matrix \(\mathbf{M}\) is constant and defined as

$$ \mathbf{M}= \begin{bmatrix} \mathbb{I}_{3\times 3} \\ \mathbf{0}_{3\times 3}\end{bmatrix} . $$
(21)

Remarks 1

As a consequence of the \(SE(3)\) Lie group derivative the constraint gradient depends only on local relative quantities. It is invariant under a superimposed rigid body transformation. Since the contact kinematics of thin beams with circular cross-sections is rather simple, \(\mathbf{M}\) has a simple expression.

3.2 Continuous formulation

We suppose that contact interactions develop along a beam segment. The normal contact pressure, denoted by \(\lambda \), has the dimensions of a force per unit length. It plays the role of a scalar Lagrange multiplier field. The space of admissible configurations is denoted by \(\mathcal{S}\) and the space of admissible Lagrange multipliers by ℳ. Finally, the space of admissible variations is written as \(\mathcal{V}\). With these notations the weak form of the equilibrium is given by the following saddle point problem: Find \(q\in \mathcal{S}\) and \(\lambda \in \mathcal{M}\) such that

$$ \delta \mathcal{W}^{\text{int}}\left (q,\delta \boldsymbol{\pi }\right )+ \delta \mathcal{W}^{\text{con}}\left (q,\lambda ,\delta \boldsymbol{\pi }\right )=\delta \mathcal{W}^{\text{ext}}\left (q, \delta \boldsymbol{\pi }\right ) \quad \forall \delta \boldsymbol{\pi } \in \mathcal{V}, $$
(22a)
$$ \delta \mathcal{W}^{\text{lag}}\left (q,\delta \lambda \right ) \geq 0 \quad \forall \delta \lambda \in \mathcal{M}, $$
(22b)

where we define

$$ \delta \mathcal{W}^{\text{con}}\left (q,\lambda ,\delta \boldsymbol{\pi }\right )=\int _{s_{a}}^{s_{b}}(\delta \boldsymbol{\pi })^{T}\mathbf{G}^{T}\lambda \,\mathrm{d}s_{1}, $$
(23a)
$$ \delta \mathcal{W}^{\text{lag}}\left (q,\delta \lambda \right )=\int _{s_{a}}^{s_{b}}( \delta \lambda - \lambda ) g\, \mathrm{d}s_{1}. $$
(23b)
  • \(\delta W^{\text{int}}\) and \(\delta W^{\text{ext}}\) are the usual contributions stemming from the internal and external work.

  • The Lagrange multiplier is defined on the slave beam.

  • Here, the potential contact region is a portion of the centerline i.e. \(\Gamma _{c}=\left [s_{a},s_{b}\right ]\) and the contact contributions in (23a) and (23b) are line integrals along segments of the slave beam. The integration boundaries \(s_{a}\) and \(s_{b}\) are also defined through the projection in Sect. 3.5.

  • The variational inequality (22b) expresses the weak version of the normal contact constraint [12, 31].

  • In Eq. (23a) the constraint forces and moments conjugated to the local arbitrary motions appear as quantities expressed in the local frame. Moreover, as a consequence of Eq. (21), contact does not induce any torque on the beams. If this effect is taken into account, the second entry of \(\mathbf{M}\) must be replaced by the lever arm between the cross-section center and the application point of the contact force, which slightly complicates the formulation.

  • All contributions only depend the relative configuration of the two beams and thus, as for the contact free formulation, the equations are insensitive to large amplitude rigid motions of the two beam system.

  • The problem is formulated on a Lie group and consistent Lie group discretization and solution schemes need to be applied. Moreover, the usual functional spaces \(\mathcal{S}\), \(\mathcal{V}\) and ℳ involved in mortar formulations have to be adapted to the present setting.

3.3 Discrete formulation

In this section, the contact finite element that describes the interaction of two beam elements as introduced in Sect. 2.3 is developed. An illustration of the discrete model is given in Fig. 3. The configuration of the element is given by

$$ q_{e}=\left (\mathbf{H}_{OA}, \mathbf{H}_{OB}, \mathbf{H}_{OD}, \mathbf{H}_{OE}\right ). $$
(24)

The discretizations of the infinitesimal motions are

$$ \delta \boldsymbol{\pi }_{C}(s_{1})= \underbrace{\mathbf{Q}\left (s_{1},\mathbf{d}_{AB}\right )}_{= \mathbf{Q}_{1}}\delta \boldsymbol{\pi }_{e{_{1}}} $$
(25a)
$$ \delta \boldsymbol{\pi }_{F}(s_{2})= \underbrace{\mathbf{Q}\left (s_{2},\mathbf{d}_{DE}\right )}_{= \mathbf{Q}_{2}}\delta \boldsymbol{\pi }_{e{_{2}}}. $$
(25b)

Let us collect the nodal variations for the contact element in a \(24\times 1\) vector, such that the field of infinitesimal motions is given by

$$ \delta \boldsymbol{\pi }(s_{1},s_{2}) = \underbrace{\begin{bmatrix} \mathbf{Q}_{1} & \mathbf{0}\\ \mathbf{0} & \mathbf{Q}_{2} \end{bmatrix}}_{=\mathbf{Q}_{e}}\delta \boldsymbol{\pi }_{e}, \qquad \delta \boldsymbol{\pi }_{e}= \begin{bmatrix} \delta \boldsymbol{\pi }_{e{_{1}}} \\ \delta \boldsymbol{\pi }_{e{_{2}}}\end{bmatrix} = \begin{bmatrix} \delta \boldsymbol{\pi }_{A} \\ \delta \boldsymbol{\pi }_{B} \\ \delta \boldsymbol{\pi }_{D} \\ \delta \boldsymbol{\pi }_{E} \end{bmatrix} . $$
(25c)
Fig. 3
figure 3

Notations for the contact element

A key question now is the choice of interpolation functions for the Lagrange multipliers and the associated finite dimensional space \(\mathcal{M}_{h}\subset \mathcal{M}\). The main criterion is the verification of the inf-sup condition to guarantee the stability of the saddle point formulation [9, 10]. It has been shown that for usual displacement based formulations, combining linear finite elements with linear interpolation functions for the Lagrange multiplier field yields a stable mortar method [11]. In this paper, the discretization of the nodal frames is based on a first order interpolation on the Lie algebra. Therefore, the Lagrange multiplier is discretized using standard linear shape functions. These are defined over the slave element:

$$ \lambda (s_{1})=\mathbf{N}(s_{1})\boldsymbol{\lambda }_{e}, $$
(26a)
$$ \delta \lambda (s_{1})=\mathbf{N}(s_{1})\delta \boldsymbol{\lambda }_{e}. $$
(26b)

\(\boldsymbol{\lambda }_{e}= \begin{bmatrix} \lambda _{A} \\ \lambda _{B}\end{bmatrix} \) and \(\delta \boldsymbol{\lambda }_{e}= \begin{bmatrix} \delta \lambda _{A} \\ \delta \lambda _{B}\end{bmatrix} \) are vectors that collect the nodal values of the Lagrange multipliers and their variations respectively. The shape functions simply are

$$ \mathbf{N}(s_{1})= \begin{bmatrix} \dfrac{s_{1}}{L_{e{_{1}}}}, & 1-\dfrac{s_{1}}{L_{e{_{1}}}}\end{bmatrix} = \begin{bmatrix} N_{A}, & N_{B}\end{bmatrix} . $$
(27)

In this paper, the nodes of the slave beam and the nodes of the Lagrange multiplier coincide and the notation is simplified accordingly for simplicity.

The discretized version of the weak form is obtained by replacing the continuous fields involved in formulation (22a) and (22b) by their finite dimensional counterparts. The contact contribution of one element to Eq. (22a) becomes

$$ \delta \mathcal{W}_{e}^{\text{con}}\left (q_{e},\boldsymbol{\lambda }_{e}, \delta \boldsymbol{\pi }_{e}\right )=\left (\delta \boldsymbol{\pi }_{e} \right )^{T}\mathbf{G}_{e}^{T}\boldsymbol{\lambda }_{e}=\left (\delta \boldsymbol{\pi }_{e}\right )^{T}\mathbf{f}_{e}^{\text{con}}, $$
(28)

where the constraint gradient of the contact element was identified as

$$ \mathbf{G}_{e}^{T}(q_{e})=\int _{s_{a}^{e}}^{s_{b}^{e}} \begin{bmatrix} \mathbf{Q}_{1}^{T} & \mathbf{0} \\ \mathbf{0} & \mathbf{Q}_{2}^{T} \end{bmatrix} \begin{bmatrix} \mathbf{G}_{1}^{T} \\ \mathbf{G}_{2}^{T} \end{bmatrix} \begin{bmatrix} N_{A} & N_{B}\end{bmatrix} \,\mathrm{d}s_{1}. $$
(29)

It expresses the orientation of the weighted nodal contact forces and moments in the local frames. Therefore, the vector \(\mathbf{f}^{\text{con}}_{e}\) represents the weighted nodal contact forces and moments in the local frame. The contribution of one element to Eq. (22b) becomes

$$ \delta \mathcal{W}_{e}^{\text{lag}}\left (q_{e},\boldsymbol{\lambda }_{e}, \delta \boldsymbol{\lambda }_{e}\right )=\left (\delta \boldsymbol{\lambda }_{e}-\boldsymbol{\lambda }_{e}\right )^{T} \mathbf{g}_{e}^{\text{con}}, $$
(30)

where \(\mathbf{g}_{e}^{\text{con}}\) is the vector of weighted normal constraints for element \(e\). Its components are computed as

$$ \mathbf{g}_{e}^{\text{con}}(q_{e})=\int _{s_{a}^{e}}^{s_{b}^{e}} \mathbf{N}^{T} g \,\mathrm{d}s_{1}. $$
(31)

The contributions in (28) and (30) are assembled in the usual finite element manner. The discrete quasi-static equilibrium of the entire two beam system is then summarized as

$$ \mathbf{f}^{\text{int}}(q)+\mathbf{B}^{T}(q)\boldsymbol{\lambda }= \mathbf{f}^{\text{ext}}(q), $$
(32a)
$$ \mathbf{0}\leq \boldsymbol{\lambda } \perp \mathbf{g}^{\text{con}}\geq \mathbf{0}, $$
(32b)

where \(q\) is the configuration of the entire discretized system, \(\boldsymbol{\lambda }\) contains all the nodal values of the Lagrange multipliers, \(\mathbf{f}^{\text{int}}\) and \(\mathbf{f}^{\text{ext}}\) are the assembled weighted nodal internal and external force vectors, \(\mathbf{B}\) is the assembled matrix of constraint gradients and \(\mathbf{g}^{\text{con}}\) is the vector of assembled weighted constraints. The ⊥ sign denotes componentwise orthogonality, meaning that \(\lambda _{i} g_{i}=0\), \(\forall i=1,2,\ldots,m\), where the subscript \(i\) refers to the components of the vectors and \(m\) is the total number of discrete constraints. In the present framework it is equal to the total number of Lagrange multiplier nodes, as opposed to collocation based formulations, where the number of contact points is arbitrary. Hence the mortar method does not suffer from potential overconstraining.

The formulation of Eqs. (32a), (32b) requires the resolution of a linear complementarity problem (LCP). In order to use a Newton type solver, problem (32a), (32b) needs to be reformulated without inequalities. This can be done by applying nonsmooth equations [12] or by using an augmented Lagrangian method [14], which is the approach selected in this work.

3.4 Augmented Lagrangian approach

Following [13], an augmented Lagrangian method is applied to solve the LCP (32a), (32b). For that purpose, the augmented multiplier \(\xi _{i}=k\lambda _{i}-p g_{i}\) is introduced and the contact contribution in (28) is replaced by the variation of the functional

$$ \mathcal{L}^{\text{con}}\left (q,\boldsymbol{\lambda }\right )=\sum _{i \in \mathcal{C}}\left (\dfrac{p}{2}g_{i}^{2}-k\lambda _{i} g_{i}- \dfrac{1}{2p}\text{dist}^{2}(\xi _{i},\mathbb{R}^{+})\right ), $$
(33)

where \(\mathcal{C}\) is the set of discrete constraints, \(k\) is a scaling factor, \(p\) is a penalty coefficient and \(\text{dist}(\mathbf{a},A)\) denotes the distance of a point \(\mathbf{a}\in \mathbb{R}^{n}\) to the convex set \(A\). Numerical parameters \(k\) and \(p\) have no influence on the solution, but their presence improves convergence. The variation of Eq. (33) is given by

$$\begin{aligned} \delta \mathcal{L}^{\text{con}}=\sum _{i\in \mathcal{C}} \textstyle\begin{cases} -\xi _{i}\delta g_{i}-k g_{i} \delta \lambda _{i} \quad \text{if } i \in \mathcal{A,} \\ \dfrac{-k^{2}}{2p}\lambda _{i}\delta \lambda _{i} \qquad \quad \, \, \, \text{if } i\in \bar{\mathcal{A}}, \end{cases}\displaystyle \end{aligned}$$
(34)

where \(\mathcal{A}=\left \{ i\in \mathcal{C}:\xi _{i}\geq 0\right \} \) defines the set of active constraints and \(\bar{\mathcal{A}}\) is its complement and thus is the set of inactive constraints.

Then, the semi discrete problem is stated as

$$ \mathbf{f}^{\text{int}}(q)-\mathbf{B}_{\mathcal{A}}^{T}(q) \boldsymbol{\xi }_{\mathcal{A}}=\mathbf{f}^{\text{ext}}(q), $$
(35a)
$$ -k\mathbf{g}_{\mathcal{A}}^{\text{con}}=\mathbf{0}, $$
(35b)
$$ -\dfrac{k^{2}}{p}\boldsymbol{\lambda }_{\bar{\mathcal{A}}}=\mathbf{0}, $$
(35c)

where \(\boldsymbol{\xi }_{\mathcal{A}}\) is a vector that collects the active augmented multipliers, \(\mathbf{g}^{\text{con}}_{\mathcal{A}}\) are the active assembled constraints, \(\boldsymbol{\lambda }_{\bar{\mathcal{A}}}\) are the inactive standard Lagrange multipliers and \(\mathbf{B}_{\mathcal{A}}\) is the matrix of active assembled constraint gradient.

The system of equations in (35a)–(35c) may now be solved using a quasi-static Lie group solver [16]. At each load step a semi-smooth Newton method is applied to obtain the solution of the resulting non-linear system. The necessary linearizations are given in the Appendix. They are not restricted to the augmented Lagrangian method and might be reused for other Newton-type algorithms.

3.5 Projection

The gap function introduced in (17) requires the determination of coordinate \(\bar{s_{2}}\) of the point located at \(\bar{\mathbf{x}}_{2}\). It is obtained by projection of the point located at \(\mathbf{x}_{1}\) on the slave beam onto the master beam. This projection is given by the following condition:

$$ \mathcal{P}_{m}=(\bar{\mathbf{x}}_{OF}(\bar{s_{2}})-\mathbf{x}_{OC}(s_{1})) \cdot \mathbf{m}(s_{1})=0, $$
(36)

where \(\mathbf{m}=\mathbf{R}_{OC}\mathbf{e}_{x}\) is the normal of the cross-section on the slave beam and \(\mathbf{e}_{x}= \begin{bmatrix} 1 & 0 & 0\end{bmatrix} ^{T}\). Thus, given a point on the slave beam, the associated point on the master side will be the intersection between the plane defined by the cross-section on the slave side and the centerline curve of the master beam. The solution to (36) is found using a Newton procedure.

3.5.1 Integration boundaries

In Eqs. (29) and (31), the potential contact region for the contact element \([s_{a}^{e},s_{b}^{e}]\in [0,L_{e{_{1}}}]\) is defined by the set of values \(s_{1}\) such that Eq. (37) has a solution \(\bar{s}_{2}\in [0,L_{e{_{2}}}]\). Therefore, as represented on Fig. 3, an integration boundary \(s_{k}^{e}\) with \(k=a,b\), is either given by a node of the slave beam or by the projection of a master node onto the slave. This projection is given by

$$ \mathcal{P}=(\mathbf{x}_{O\beta }-\mathbf{x}_{OC}(s_{k}^{e}))\cdot \mathbf{m}(s_{k}^{e})=0, \quad \beta =D,E ,$$
(37)

where \(\mathbf{x}_{O\beta }\) denotes the position of a master node.

Note that this way of computing the potential contact region fails when the two beam elements are perfectly orthogonal. This special case is not addressed in this paper.

4 Numerical tests

In this section the formulation is tested in three numerical examples: a patch test, a planar cantilever beam test and a \(3D\) twisting test. These are implemented in the research code Odin [50]. The convergence criterion for the Newton algorithm is

$$ \mathrm{Convergence}(\mathbf{r}) = \frac{\|\mathbf{r}\|}{\dfrac{ \sum _{k}\|\mathbf{r}_{k}^{1}\|}{N_{1}} + \cdots + \dfrac{\sum _{k}\|\mathbf{r}_{k}^{n}\|}{N_{n}} + 10^{-12}}< \text{tol}_{r} \;\;\lor \;\; \|\mathbf{r}\| < \text{tol}_{a}, $$
(38)

where \(\mathbf{r}\) is the residual, \(\mathbf{r}_{k}^{i}\) denotes the unassembled contribution to the total residual of element \(k\) of type \(i\) (i.e. beam element or contact element), \(N_{i}\) is the number of elements of type \(i\), and \(\|\cdot \|\) is the \(L^{2}\) norm, \(\text{tol}_{r}\) and \(\text{tol}_{a}\) are relative and absolute tolerances. Convergence is evaluated separately for the equilibrium of forces in Eq. (35a) and the constraints in Eq. (35b) and in Eq. (35c). These components are denoted by \(\mathbf{r}_{f}\) and \(\mathbf{r}_{c}\), respectively, and the actual convergence criterion is given by

$$ \mathrm{Convergence}(\mathbf{r}_{f}) \;\;\land \;\; \mathrm{Convergence}(\mathbf{r}_{c}). $$
(39)

The following values for the tolerances were used in the numerical examples: \(10^{-4}\) for the relative tolerance of the equilibrium of forces, \(10^{-2}\) for the relative tolerance of the constraints, \(10^{-7}\) for the absolute tolerance of the equilibrium of forces and \(10^{-5}\) for the absolute tolerance of the constraints.

In two of the following numerical examples, the spatial convergence of the proposed algorithm is analyzed. For that purpose the error on the beam centerline is studied. It is computed as the sum of the errors for each individual beam. For one beam the error is given by

$$ \mathrm{error} = \sqrt{ \frac{\int _{0}^{L} \| \mathbf{x}_{OC} - \mathbf{x}^{\text{ref}}_{OC} \|^{2} \; \mathrm{d} s}{\int _{0}^{L} \| \mathbf{x}^{\text{ref}}_{OC} \|^{2} \; \mathrm{d} s}}, $$
(40)

where \(\mathbf{x}_{OC}^{\text{ref}}\) is a reference solution for the beam centerline. In the investigated examples, it is taken as the numerical solution obtained from a very fine discretization.

4.1 Patch test

First a simple patch test is carried out. Two cantilever beams of radius \(r=5\text{ cm}\) and length \(L=1\text{ m}\) are placed on top of each other and separated by a distance \(\epsilon r\). The mechanical properties of both beams are: axial stiffness \(EA=39.27\) KN, shear stiffnesses \(GA_{22}=GA_{33}=13.09\) KN, torsional stiffness \(GJ=16.36\) Nm2 and bending stiffnesses \(EI_{22}=EI_{33}=24.54\) Nm2. The beams are clamped as shown in Fig. 4a. A constant distributed load of \(p=100\) N/m is applied on both beams in opposite directions, such that they are compressed together. In the exact solution, the beams remain straight and the compressive contact pressure takes the value \(\lambda =p\). If \(\epsilon \) is chosen equal to 0 the problem becomes trivial from a numerical point of view and the solution converges after 1 iteration independently of any other parameter. For \(\epsilon =10^{-10}\), the results obtained for the non-matching discretization are shown on Fig. 4b. The contact pressure distribution is uniform, which indicates that the methodology passes the patch test.

Fig. 4
figure 4

(a) Schematic description of the patch test. (b) Resulting Lagrange multiplier field

4.2 Cantilever

In this example, the behavior of a cantilever beam of radius \(r=1\text{ mm}\) and length \(L=0.3\text{ m}\) entering contact with a rigid substrate separated by \(0.5\text{ mm}\) from the beam is studied. The material parameters of the beam are given by: axial stiffness \(EA=0.628\) MN, shear stiffnesses \(GA_{22}=GA_{33}=0.242\) MN, torsional stiffness \(GJ=0.12\) Nm2 and bending stiffnesses \(EI_{22}=EI_{33}=0.16\) Nm2. A schematic description of the test is given in Fig. 5a.

Fig. 5
figure 5

(a) Schematic description of the cantilever test. (b) Final configuration obtained for a discretization of 16 elements. (c) Lagrange multipliers for the thin beam at different load steps with a discretization of 256 elements. (d) Lagrange multipliers for a beam with a bigger cross-section radius as a comparison. (e) Resultant contact force for different number of elements. (f) Spatial convergence analysis

The distributed load \(p\) is gradually increased until it reaches a maximum of 10 N/m. The cantilever beam first enters contact with its tip. At this stage, it is subject to a point contact force, which is observed in Fig. 5c. As the load is increased it switches to line contact. At the transition point between contact and no-contact, a discontinuity in the contact pressure followed by a rapid exponential decay is observed, see Fig. 5c. That steep decay tends towards a point force when the importance of shearing deformation becomes negligible. As a comparison, Fig. 5d shows the Lagrange multipliers for a case, where the effect of shearing is more pronounced. For this case the following parameters where adopted: radius \(r=5\text{ cm}\), length \(L=1\text{ m}\), separated by \(2.5\text{ cm}\), axial stiffness \(EA=39.27\) KN, shear stiffnesses \(GA_{22}=GA_{33}=13.09\) KN, torsional stiffness \(GJ=16.36\) Nm2, bending stiffnesses \(EI_{22}=EI_{33}=24.54\) Nm2 and maximum load \(p=500\text{ Pa}\). Far from the transition region the contact pressure takes the value of the applied distributed load. As shown by analytic solutions of similar problems, in the contact zone, the curvature is imposed by the geometry of the rigid element, Fig. 5b, and such effects are to be expected [18, 49]. Figure 5e shows the spatial convergence of the resultant contact force, obtained by integrating the Lagrange multiplier along the contact region, whereas Fig. 5f shows the discretization error as defined in Eq. (40). A convergence rate of around 2, characteristic for constant strain elements, is obtained. The numerical results indicate that the mortar finite element method is able to deal well with discontinuities present in the contact pressure. As one would expect from a linear finite element description, this discontinuity induces oscillations in the Lagrange multipliers. These do not affect the optimal spatial convergence of the mortar method.

4.3 Twisting

In this example the twisting of two beams is studied. The beams are both initially straight, of length \(L=1\text{ m}\) and radius \(r=1\text{ mm}\), and they are separated by a distance of \(0.5\text{ mm}\). The material parameters of the beams are the same as those for the thin beam in the previous section. They are fully clamped on one end and fixed to a common spherical joint on the other end, which is actuated by a hinge. In that manner the end nodes of the beam travel on an imposed circular path without being constrained in their translation which could potentially force the beams to buckle away from each other. The two beams take a shape that approaches a double helix and contact occurs along a continuous line.

In Figs. 7 and 8, the results for the beams being submitted to a total of four turns are shown. In this case both beams are discretized using 32 beam elements. As an indication, for a total number of 2400 load steps the average number of global Newton iterations is 1.8 for a maximum of 5 iterations. As expected from the method, the weighted contact constraints in Fig. 8b are satisfied exactly. In Fig. 8a some oscillations are observed near the transitions between contact and no-contact states, which can be attributed to the presence of a discontinuity in the contact pressure as in the previous example. The proposed method is able to handle the nonsmooth distributed contact force at the cost of oscillations. These can be reduced by choosing appropriate numerical parameters. Improvements in the treatment of these oscillations will be addressed in future works.

In this particular example a spatial convergence rate above 2 is obtained; see Fig. 6a. This superconvergence could be due to the \(SE(3)\) interpolation scheme, which approximates the beam centerline by piecewise helices [1, 4]. These are particularly well-suited for this kind of problem and an accurate representation of the centerline is obtained using a small number of elements. The influence of the discretization and the size of the load step on the convergence of the global Newton scheme is shown on Fig. 6b. For this study the terms related to the linearization of the shape functions and the normal in equation (56) were neglected in the tangent stiffness matrix. Note that for certain cases a high number of maximum iterations is observed. However, it occurs only once during the simulation when the two beams enter contact and the Lagrangian multipliers are initialized to zero. Subsequently, the multiplier computed at the previous load step is taken as initial value. Furthermore, in this paper, a general Tymoshenko type model, not tailored to very slender beams, is used which might have a negative impact on convergence due to ill-conditioning [21], regardless of the approach selected for solving contact.

Fig. 6
figure 6

(a) Spatial convergence analysis for the \(3D\) twisting example. (b) Average and maximum number of Newton iterations. For the influence of the discretization the number of load steps is set to 60. For the influence of the load step the number of elements is set to 8. For both analyses the beams are submitted to one turn

Fig. 7
figure 7

Results obtained for the twisting of two beams. The beams have been scaled for visualization purposes

Fig. 8
figure 8

(a) Final distributed contact force as a function of the arc length parameter for the twisting example. (b) Final integrated gap

4.4 Practical applicability

The test cases presented in this section showed the potential of mortar for solving line-to-line contact problems with relatively coarse meshes. However, the imposition of the weighted constraint (31), which is needed for variational consistency, imposes limitations. Indeed, for the method to detect penetration, two beam elements need to be in contact over a sufficiently large fraction of the computed contact patch. Otherwise the weighted gap will yield a positive value despite the presence of local penetration. The size of the patch depends on the relative configuration of the contacting beam elements and the discretization. As a consequence, an extremely fine discretization would be needed to deal with situations that involve pointwise interactions or contact over very small regions, shorter than the length of the beam diameter. Reduced models, such as the beam, are usually used in applications with several components interacting with each other and the interest mostly lies in the system behavior of assemblies. In such a context relatively coarse discretizations are commonly preferred and thus a point-to-point type model should be applied to complement the mortar’s shortcomings when contact is highly localized.

5 Conclusion and perspectives

In this paper the mortar formulation of frictionless line-to-line contact among beams with circular cross-sections is addressed and combined with the \(SE(3)\) formalism for flexible multibody systems.

First, the contact kinematics in the local frame approach is discussed. Then the equilibrium equations are formulated as a saddle point problem, where the contact constraints are enforced in a weak sense. As a consequence of the \(SE(3)\) Lie group formalism the contact forces are expressed in the local frames of the beams and the equilibrium equations enjoy similar invariance properties as the contact free case.

A linear discretization scheme is applied for the Lagrange multipliers. The discretized equations are solved using the Augmented Lagrangian method. The behavior of the proposed method is tested in numerical examples and optimal spatial convergence rates and the absence of over-constrainment are shown for these, as expected for the mortar approach. The ability of the method to capture nonsmooth phenomena involved in beam contact mechanics is studied.

Future efforts will concentrate on the combination of this mortar formulation with a point-to-point contact formulation, the upscaling of the methodology towards larger sized problems involving several beams and the extension of the formulation to friction. Moreover, the invariance properties of the constraint gradient and the iteration matrix deserve more thorough numerical investigations in terms of numerical advantages they could provide.