A mortar formulation for frictionless line-to-line beam contact

This paper describes the quasi-static 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 cross-sections is proposed for the computation of contact regions. It is combined with the geometrically exact beam formalism on the Lie group SE(3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$SE(3)$\end{document}. 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, multi-wire cable structures for the automotive industry or endoscopes for medical applications 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 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 [24][25][26][27], so that the models are formulated on the Lie group R 3 × 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 [37][38][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 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 [9][10][11][12][13][40][41][42][43][44][45][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 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 lineto-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.

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 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 H OC ∈ SE (3). It is composed of a translation x OC ∈ R 3 which represents the Cartesian coordinates of the frame center C in the frame {O} and a rotation 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 δπ C C ∈ 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 se(3), which is isomorphic to R 6 through the map R 6 → se(3) : δπ C C → δπ C C : where the translation part δx C C = R T OC δx O C is an element of R 3 . The rotation part δ θ C C = R T OC δR OC belongs the Lie algebra of SO(3), denoted so (3). This Lie algebra is the set of 3 Fig. 1 Kinematic description of the beam configuration using frame transformations as unknowns by 3 skew-symmetric matrices and is isomorphic to Whether the (·) 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 δx C C and δθ C C or in the inertial frame {O} as δx O C and δθ O C . 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 ∈ [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 R → SE(3) : s → 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 elementf C (s) of the Lie algebra 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 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 S denotes the space of admissible configurations q and V the space of admissible infinitesimal motions, the weak form of the equilibrium is given by: Find q ∈ S such that The virtual work done by the external forces is where f ext C (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 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 S h ⊂ S and V h ⊂ 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 = (H OA , H OB ). 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 d 0 AB 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, 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), d AB (and as a consequence 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 δπ e = δπ A δπ B 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 Q and P must be adapted. Here, since strains are constant elementwise 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 centerlinex OF is associated to a point on the slave centerline x OC by means of a projection, which will be described in Sect. 3.5. Similarly, rotation matrixR OF is evaluated at the same location asx 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 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), the local constraint gradient is obtained as where 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)R T OF R OC = R F C is the relative orientation between the local slave and the local master beam. Thus n F =R T OF R OC n C may be seen as the same normal vector as n C , but expressed in the local frame of the master. The matrix M is constant and defined as Remarks 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, M has a simple expression.

Continuous formulation
We suppose that contact interactions develop along a beam segment. The normal contact pressure, denoted by λ, 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 S and the space of admissible Lagrange multipliers by M. Finally, the space of admissible variations is written as V. With these notations the weak form of the equilibrium is given by the following saddle point problem: Find q ∈ S and λ ∈ M such that where we define -δW int and δW 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. c = [s a , s b ] 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 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 S, V and M 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 × 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 M h ⊂ 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: 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 δW con e (q e , λ e , δπ e ) = (δπ e ) T G T e λ e = (δπ e ) T f con e , (28) 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 f con e represents the weighted nodal contact forces and moments in the local frame. The contribution of one element to Eq. (22b) becomes δW lag e (q e , λ e , δλ e ) = (δλ e − λ e ) T g con e , where g con e 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 quasi-static equilibrium of the entire two beam system is then summarized as where q is the configuration of the entire discretized system, λ contains all the nodal values of the Lagrange multipliers, f int and f ext are the assembled weighted nodal internal and external force vectors, B is the assembled matrix of constraint gradients and g con is the vector of assembled weighted constraints. The ⊥ sign denotes componentwise orthogonality, meaning that λ i g i = 0, ∀i = 1, 2, . . . , 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 ξ i = kλ i − pg i is introduced and the contact contribution in (28) is replaced by the variation of the functional where C is the set of discrete constraints, k is a scaling factor, p is a penalty coefficient and dist(a, A) denotes the distance of a point a ∈ 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 A = {i ∈ C : ξ i ≥ 0} defines the set of active constraints andĀ is its complement and thus is the set of inactive constraints.
Then, the semi discrete problem is stated as where ξ A is a vector that collects the active augmented multipliers, g con A are the active assembled constraints, λĀ are the inactive standard Lagrange multipliers and B 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.

Projection
The gap function introduced in (17) requires the determination of coordinates 2 of the point located atx 2 . It is obtained by projection of the point located at x 1 on the slave beam onto the master beam. This projection is given by the following condition: where m = R OC e x is the normal of the cross-section on the slave beam and e x = 1 0 0 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.

Integration boundaries
In Eqs. (29) and (31), the potential contact region for the contact element [s e a , s e b ] ∈ [0, L e 1 ] is defined by the set of values s 1 such that Eq. (37) has a solutions 2 ∈ [0, L e 2 ]. Therefore, as represented on Fig. 3, an integration boundary s e k 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 x Oβ 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 r is the residual, r i k 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 · is the L 2 norm, tol r and 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 r f and 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 x ref OC 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 cm and length L = 1 m are placed on top of each other and separated by a distance 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 λ = p. If 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 = 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.

Cantilever
In this example, the behavior of a cantilever beam of radius r = 1 mm and length L = 0.3 m entering contact with a rigid substrate separated by 0.5 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 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 cm, length L = 1 m, separated by 2.5 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 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 m and radius r = 1 mm, and they are separated by a distance of 0.5 mm. The material parameters of the beams are the same as those for the thin beam in the previous  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.

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.

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. and the associated jacobian is given by Integral 31 is then computed as The integration boundaries s e a and s e b 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.

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:

Linearization of the active constraints
Starting from here, the linearizations are given for one contact element to keep notations simple. Thus in the following 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 g was computed in Sect. 3.1. For the shape functions of the active Lagrangian multipliers one has ζ 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 s e a and s e b come from the varying domain of the integral. They can be expressed in terms of the configuration increments π . As detailed in Sect. 3.5.1, the integration boundaries may be slave nodes, in which case we have s e k = 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 s e k = S k π , where k = a, b. The projection condition is of the form P = P (q e , s k (q e )) = 0.
By the implicit function theorem its linearization is given by P = P q π + ∂P ∂s 1 s 1 =s e k s e k = 0, where P q is the derivative of P with respect to the configuration variables. and where x C CF = R T OC (x OF − x OC ) is the relative configuration vector expressed in the local frame of the slave beam.

Tangent stiffness
The active constraint variation may be linearized as follows: N A and J nat have already been derived in the previous section. A closed form expression for 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 Q e and G T in the iteration matrix. However, these terms probably cannot be neglected for more advanced numerical tests.