Abstract
This paper describes the quasistatic formulation of frictionless line contact between flexible beams by employing the mortar finite element approach. Contact constraints are enforced in a weak sense along the contact region using Lagrange multipliers. A simple projection appropriate for thin beams with circular crosssections is proposed for the computation of contact regions. It is combined with the geometrically exact beam formalism on the Lie group \(SE(3)\). Interestingly, this framework leads to a constraint gradient and a tangent stiffness invariant under rigid body transformations. The formulation is tested in some numerical examples.
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, multiwire cable structures for the automotive industry or endoscopes for medical applications^{Footnote 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, Kirchhofftype 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 crosssection [7, 8]. This type of formulation involves the handling of finite rotations, which belong to a nonlinear space, \(SO(3)\). Often they are treated as decoupled from the translation field [24–27], 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 beamtobeam contact are relatively scarce. Due to their particular geometric features, beams of circular crosssection 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 beamtobeam contact were dedicated to pointtopoint 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 [37–39] 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 degreesoffreedom (DOFs) should be adjusted to the density of collocation points to avoid overconstrainment. 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 pointtopoint and linetoline contact are covered by switching between the two contact models.
Mortar methods for \(3D\) continuum contact mechanics have been extensively studied in the literature [9–13, 40–46]. 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 overconstrained 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 linetoline beam contact. As a perspective, the ultimate goal is a formulation that treats both pointtopoint and linetoline 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.
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.
Frames
Reference frames are everpresent 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:
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
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}\):
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 skewsymmetric matrices and is isomorphic to \(\mathbb{R}^{3}\) through the map \(\mathbb{R}^{3}\rightarrow \mathfrak{so}(3):\mathbf{a}\rightarrow \widetilde{\mathbf{a}}\)
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.
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:
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
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].
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
The virtual work done by the external forces is
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
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).
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
The interpolation formula follows as
It is based on the exponential and logarithm maps, for more details see [1, 3, 4]. The discretized strain is then given by
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].
In general the discretizations of infinitesimal motions and deformation gradients can be written
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:
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.
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
and the variations are collected in the 12 dimensional vector
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:
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 crosssections, 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),
the local constraint gradient is obtained as
where \(\mathbf{n}^{C}\) denotes the unit outward normal to the slave beam expressed in its local frame. It is defined as:
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
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 crosssections is rather simple, \(\mathbf{M}\) has a simple expression.
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
where we define

\(\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 crosssection 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.
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
The discretizations of the infinitesimal motions are
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
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 infsup 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:
\(\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
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
where the constraint gradient of the contact element was identified as
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
where \(\mathbf{g}_{e}^{\text{con}}\) is the vector of weighted normal constraints for element \(e\). Its components are computed as
The contributions in (28) and (30) are assembled in the usual finite element manner. The discrete quasistatic equilibrium of the entire two beam system is then summarized as
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.
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
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
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
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 quasistatic Lie group solver [16]. At each load step a semismooth Newton method is applied to obtain the solution of the resulting nonlinear system. The necessary linearizations are given in the Appendix. They are not restricted to the augmented Lagrangian method and might be reused for other Newtontype algorithms.
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:
where \(\mathbf{m}=\mathbf{R}_{OC}\mathbf{e}_{x}\) is the normal of the crosssection 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 crosssection on the slave side and the centerline curve of the master beam. The solution to (36) is found using a Newton procedure.
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
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.
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
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
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
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.
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\) Nm^{2} and bending stiffnesses \(EI_{22}=EI_{33}=24.54\) Nm^{2}. 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 nonmatching discretization are shown on Fig. 4b. The contact pressure distribution is uniform, which indicates that the methodology passes the patch test.
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\) Nm^{2} and bending stiffnesses \(EI_{22}=EI_{33}=0.16\) Nm^{2}. A schematic description of the test is given in Fig. 5a.
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 nocontact, 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\) Nm^{2}, bending stiffnesses \(EI_{22}=EI_{33}=24.54\) Nm^{2} 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.
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 nocontact 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 wellsuited 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 illconditioning [21], regardless of the approach selected for solving contact.
Practical applicability
The test cases presented in this section showed the potential of mortar for solving linetoline 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 pointtopoint type model should be applied to complement the mortar’s shortcomings when contact is highly localized.
Conclusion and perspectives
In this paper the mortar formulation of frictionless linetoline contact among beams with circular crosssections 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 overconstrainment 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 pointtopoint 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.
Notes
 1.
Application examples treated in the THREAD project: https://threadetn.eu
References
 1.
Sonneville, V., Cardona, A., Brüls, O.: Geometrically exact beam finite element formulated on the Special Euclidean group \(SE(3)\). Comput. Methods Appl. Mech. Eng. 268, 451–474 (2014)
 2.
Sonneville, V.: A geometric local frame approach for flexible multibody systems. PhD Thesis, University of Liège (2015)
 3.
Sonneville, V., Cardona, A., Brüls, O.: Geometric interpretation of a nonlinear finite element on the Lie group \(SE(3)\). Arch. Mech. Eng. 61, 305–329 (2014)
 4.
Borri, M., Botasso, C.: An intrinsic beam model based on a helicoidal approximation. Int. J. Numer. Methods Eng. 37, 2267–2289 (1994)
 5.
Kirchhoff, G.: Über das Gleichgewicht und die Bewegung eines unendlich dünnen elastischen Stabes. J. Reine Angew. Math. 56, 285–313 (1859)
 6.
Love, A.E.H.: A Treatise on the Mathematical Theory of Elasticity, 4th edn. Cambridge University Press, UK (1927)
 7.
Reissner, E.: On onedimensional largedisplacement finitestrain beam theory. Stud. Appl. Math. 2, 87–95 (1973)
 8.
Simo, J.C.: A finite strain beam formulation. The threedimensional dynamic problem. Part I. Comput. Methods Appl. Mech. Eng. 49, 55–70 (1985)
 9.
Belgacem, F.B.: The Mortar finite element method with Lagrange multipliers. Numer. Math. 84, 173–197 (1999)
 10.
Wohlmuth, B.: Variationally consistent discretization schemes and numerical algorithms for contact problems. Acta Numer. 20, 569–734 (2011)
 11.
Belgacem, F.B., Renard, Y.: Hybrid finite element methods for the Signorini problem. Math. Comput. 72, 1117–1145 (2003)
 12.
Popp, A.: Mortar methods for computational contact mechanics and general interface problems. PhD Thesis, Teschnische Universität München (2012)
 13.
Cavalieri, F.J., Cardona, A.: An augmented Lagrangian technique combined with a mortar algorithm for modeling mechanical contact problems. Int. J. Numer. Methods Eng. 93(4), 420–442 (2012)
 14.
Alart, P., Curnier, A.: A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput. Methods Appl. Mech. Eng. 92, 353–375 (1990)
 15.
Lolić, D., Zupan, D., Brojan, M.: A consistent strainbased beam element with quaternion representation of rotations. Comput. Mech. 65, 1397–1412 (2020)
 16.
Brüls, O., Cardona, A., Arnold, M.: Lie group generalized\(\alpha \) time integration of constrained flexible multibody systems. Mech. Mach. Theory 48, 121–138 (2012)
 17.
Meier, C., Popp, A., Wall, W.A.: A finite element approach for the linetoline contact interaction of thin beams with arbitrary orientation. Comput. Methods Appl. Mech. Eng. 308, 377–413 (2016)
 18.
Denoël, V.: Advantages of a semianalytical approach for the analysis of an evolving structure with contacts. Commun. Numer. Methods Eng. 0, 1–6 (2002)
 19.
Meier, C., Popp, A., Wall, W.A.: A lockingfree finite element formulation and reduced models for geometrically exact Kirchhoff rods. Comput. Methods Appl. Mech. Eng. 290, 314–341 (2015)
 20.
Meier, C., Popp, A., Wall, W.A.: An objective 3D large deformation finite element formulation for geometrically exact curved Kirchhoff rods. Comput. Methods Appl. Mech. Eng. 278, 445–478 (2014)
 21.
Meier, C., Wall, W.A., Popp, A.: Geometrically exact finite element formulations for curved slender beams: KirchhoffLove theory vs. Simo–Reissner theory. Arch. Comput. Methods Eng. 26, 163–243 (2019)
 22.
Audoly, B., Pomeau, Y.: Elasticity and Geometry—From Hair Curls to the Nonlinear Response of Shells. Oxford University Press, Oxford, UK (2006)
 23.
Bertails, F., Audoly, B., Cani, M., Querleux, B., Leroy, F., Lévêque, J.: SuperHelices for Predicting the Dynamics of Natural Hair. In: Proceedings of the ACM SIGGRAPH (2006)
 24.
Jelenić, G., Crisfield, M.A.: Geometrically exact 3D beam theory: implementation of a straininvariant finite element for statics and dynamics. Comput. Methods Appl. Mech. Eng. 171, 141–171 (1999)
 25.
Crisfield, M.A., Jelenić, G.: Objectivity of strain measures in the geometrically exact threedimensional beam theory and its finiteelement implementation. Proc. R. Soc. Lond. A 455, 1125–1147 (1999)
 26.
Cardona, A., Géradin, M.: A beam finite element nonlinear theory with finite rotations. Int. J. Numer. Methods Eng. 26, 2403–2438 (1988)
 27.
Cottanceau, E., Thomas, O., Véron, P., Alochet, M., Deligny, R.: A finite element/quaternion/asymptotic numerical method for the \(3D\) simulation of flexible cables. Finite Elem. Anal. Des. 139, 14–34 (2018)
 28.
Grazioso, S., Di Gironimo, G., Siciliano, B.: A geometrically exact model for soft continuum robots: the finite element deformation space formulation. Soft Robot. 6, 790–911 (2019)
 29.
Jödicke, R., Jungnickel, U., Müller, A.: Lie group modeling of nonlinear helical beam elements. In: Proceedings of the ASME IDETC/CIE (2014)
 30.
Sonneville, V., Brüls, O., Bauchau, O.: Interpolation schemes for geometrically exact beams: a motion approach. Int. J. Numer. Methods Eng. 112, 1129–1153 (2017)
 31.
Kikuchi, N., Oden, J.T.: Contact Problems in Elasticity. SIAM Studies in Applied and Numerical Mathematics (1988)
 32.
Laursen, T.A.: Computational Contact and Impact Mechanics. Springer, Berlin (2003)
 33.
Wriggers, P., Zavarise, G.: On contact between threedimensional beams undergoing large deflections. Commun. Numer. Methods Eng. 13, 42–438 (1997)
 34.
Zavarise, G., Wriggers, P.: Contact with friction between beams in 3D space. Int. J. Numer. Methods Eng. 49, 977–1006 (2000)
 35.
Chamekh, M., ManiAouadi, S., Moakher, M.: Modeling and numerical treatment of elastic rods with frictionless selfcontact. Comput. Methods Appl. Mech. Eng. 198, 3751–3764 (2009)
 36.
Meier, C., Grill, M.J., Wall, W.A., Popp, A.: Geometrically exact beam elements and smooth contact schemes for the modeling of fiberbased materials and structures. Int. J. Solids Struct. 154, 124–146 (2018)
 37.
Durville, D.: Contactfriction modeling within elastic beam assemblies: an application to knot tightening. Comput. Mech. 49, 687–707 (2012)
 38.
Durville, D.: Numerical simulation of entangled materials mechanical properties. J. Mater. Sci. 40, 5941–5948 (2005)
 39.
Bertails, F., Cadoux, F., Daviet, G., Acary, V.: A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies. ACM Trans. Graph. 30, 1–14 (2011)
 40.
Hueber, S., Wohlmuth, B.: A primal–dual active set strategy for nonlinear multibody contact problems. Comput. Methods Appl. Mech. Eng. 194, 3147–3166 (2005)
 41.
Farah, P., Wall, W.A.: A mortar finite element approach for point, line and surface contact. Int. J. Numer. Methods Eng. 00, 1–44 (2017)
 42.
Gitterle, M., Popp, A., Gee, M.W., Wall, W.A.: Finite deformation frictional mortar contact using a semismooth Newton method with consistent linearization. Int. J. Numer. Methods Eng. 84, 543–571 (2010)
 43.
Popp, A., Wall, W.A.: Dual mortar methods for computational contact mechanics – overview and recent developments. GAMMMitt. 37(1), 66–84 (2014)
 44.
Puso, M.A., Laursen, T.A.: A mortar segmenttosegment contact method for large deformation solid mechanics. Comput. Methods Appl. Mech. Eng. 193, 601–629 (2004)
 45.
Raviart, P.A., Thomas, J.M.: Primal hybrid finite element methods for 2nd order elliptic equations. Math. Comput. 31, 391–413 (1977)
 46.
Yang, B., Laursen, T.A., Meng, X.: Two dimensional mortar contact methods for large deformation frictional sliding. Int. J. Numer. Methods Eng. 62, 1183–1225 (2005)
 47.
Han, S., Bauchau, O.: On the global interpolation of motion. Comput. Methods Appl. Mech. Eng. 337, 352–386 (2018)
 48.
Bauchau, O., Han, S.: Interpolation of rotation and motion. Multibody Syst. Dyn. 31, 339–370 (2014)
 49.
Chen, JS.: On the contact behavior of a buckled Timoshenko beam constrained laterally by a plane wall. Acta Mech. 222, 225–232 (2011)
 50.
Odin: a research code for the simulation of nonsmooth flexible multibody systems. University of Liège, Department of Aerospace and Mechanical Engineering. To be released as opensource under the Apache v2 license
 51.
Arnold, M., Brüls, O., Linn, J.: THREADnumerical modeling of highly flexible structures for industrial applications. In: Celledoni, E., Münch, A. (eds.) Mathematics with industry: driving innovationECMI Annual Report 2019, pp. 20–25 (2019)
Acknowledgements
This work is partly funded by the Robotix Academy project of the Greater Region. The first author would like to thank the Fraunhofer Institute for Industrial Mathematics Kaiserslautern for its financial support.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Affiliations
Corresponding author
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
6.1 Integration of contact contributions
Integrals (29) and (31) involved in formulation (35a)–(35c) have moving boundaries. Indeed \(s_{a}\) and \(s_{b}\), or rather \(s_{a}^{e}\) and \(s_{b}^{e}\) for the discretized element, depend on the configuration, which needs to be taken into account during the linearization procedure. In this section we show how to compute these integrals using Gauss points and make the necessary preparations for the linearization procedure. For that purpose, the case of the element constraint vector is detailed here. First we introduce the mapping \(\zeta (\eta ):[1,1]\rightarrow [s_{a}^{e},s_{b}^{e}]:\eta \rightarrow s\) to rescale the integrals. Its expression is
and the associated jacobian is given by
Integral 31 is then computed as
The integration boundaries \(s_{a}^{e}\) and \(s_{b}^{e}\) are found through the projection procedure described in Sect. 3.5.1. Similarly, in Eq. (44), the projection explained in Sect. 3.5 is used to associate points on the master beam to Gauss points on the slave beam, where the integration is performed.
6.2 Linearizations for the contact element
The linearization of the Lagrange functional variation is given by
The first term hides the linearization of the active constraints as can be seen when written explicitly:
6.2.1 Linearization of the active constraints
Starting from here, the linearizations are given for one contact element to keep notations simple. Thus in the following \(\mathcal{A}\) will refer to the set of element nodes associated to active constraints. These element contributions are then assembled to a global iteration matrix. The linearization of the constraint is given by
\(\Delta g\) was computed in Sect. 3.1. For the shape functions of the active Lagrangian multipliers one has
\(\zeta \) is a quantity related to the slave beam and its linearization can immediately be obtained from Eq. (41) as
Thereby, it is linked to the linearizations of the integration boundaries. The same holds for the jacobian of the natural transformation. Its linearization is given by
\(\Delta s_{a}^{e}\) and \(\Delta s_{b}^{e}\) come from the varying domain of the integral. They can be expressed in terms of the configuration increments \(\Delta \boldsymbol{\pi }\). As detailed in Sect. 3.5.1, the integration boundaries may be slave nodes, in which case we have \(\Delta s_{k}^{e}=0\), where \(k=a,b\). They may also be implicitly defined by the projection of a master node onto a slave beam. In that case Eq. (37) is linearized. The aim is to find the operator that links the linearization of the boundary to the linearization of the configuration as \(\Delta s_{k}^{e} =\mathbf{S}_{k}\Delta \boldsymbol{\pi }\), where \(k=a,b\). The projection condition is of the form
By the implicit function theorem its linearization is given by
where \(\mathcal{P}_{q}\) is the derivative of \(\mathcal{P}\) with respect to the configuration variables. Solving for \(\Delta s_{k}^{e}\) gives
and thus \(\mathbf{S}_{k}=\left (\dfrac{\partial \mathcal{P}}{\partial s_{1}} \biggr _{s_{1}=s_{k}^{e}}\right )^{1}\mathcal{P}_{q}\). Finally, the derivatives are computed as
and
where \(\mathbf{x}_{CF}^{C}=\mathbf{R}_{OC}^{T}(\bar{\mathbf{x}}_{OF} \mathbf{x}_{OC})\) is the relative configuration vector expressed in the local frame of the slave beam.
6.2.2 Tangent stiffness
The active constraint variation may be linearized as follows:
\(\Delta \mathbf{N}_{\mathcal{A}}\) and \(\Delta J_{\mathit{nat}}\) have already been derived in the previous section. A closed form expression for \(\Delta \mathbf{Q}_{e}\) may be found in [1] and the linearization of the constraint gradient is given by
where the same notations as in Sect. 3.1 apply. In the present contribution all the numerical results were obtained without including both \(\Delta \mathbf{Q}_{e}\) and \(\Delta \mathbf{G}^{T}\) in the iteration matrix. However, these terms probably cannot be neglected for more advanced numerical tests.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Bosten, A., Cosimo, A., Linn, J. et al. A mortar formulation for frictionless linetoline beam contact. Multibody Syst Dyn (2021). https://doi.org/10.1007/s11044021097995
Received:
Accepted:
Published:
Keywords
 Beam contact
 Mortar
 Special Euclidean group