{'fraction_non_alphanumeric': 0.03590584470104804, 'fraction_numerical': 0.012287997471280263, 'mean_word_length': 4.792310332990045, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Finite Element Approach for the Line-to-Line Contact Interaction of Thin Beams with Arbitrary Orientation

The objective of this work is the development of a novel finite element formulation describing the contact interaction of slender beams in complex 3D configurations involving arbitrary beam-to-beam orientations. It is shown in a mathematically concise manner that standard beam contact models based on a point-wise contact force fail to describe a considerable range of configurations, which are, however, likely to occur in practical applications. On the contrary, the formulation proposed here models beam-to-beam contact by means of distributed line forces, a procedure that is shown to be applicable for arbitrary geometrical configurations. The proposed formulation is based on a Gauss-point-to-segment type contact discretization and a penalty regularization of the contact constraint. By means of detailed theoretical and numerical investigations, it is shown that this approach is more suitable for beam contact than possible alternatives based on mortar type contact discretizations or constraint enforcement by means of Lagrange multipliers. The proposed formulation is enhanced by a consistently linearized integration interval segmentation avoiding numerical integration across strong discontinuities. In combination with a smoothed contact force law and the employed C1-continuous beam elements, this procedure drastically reduces the numerical integration error, an essential prerequisite for optimal spatial convergence rates. The resulting line-to-line contact algorithm is supplemented by contact contributions of the beam endpoints, which represent boundary minima of the underlying minimal distance problem. Finally, a series of numerical test cases is analyzed in order to investigate the accuracy and consistency of the proposed formulation regarding integration error, spatial convergence behavior and resulting contact force distributions.


Introduction
There exist many fields of application, where mechanical system behavior is crucially determined by slender fiberor rod-like components. In technical applications, such fibers occur for example in industrial webbings, high-tensile ropes and cables, fiber-reinforced composite materials or synthetic polymer materials. Furthermore, also the fibers in biological systems such as muscles and biological tissue or the filaments in biopolymer networks [4] can be identified as slender mechanical components of this type. In most cases, these fibers can be modeled with sufficient accuracy by applying a 1D beam theory. In the last three decades, many different types of beam element formulations have been proposed in order to discretize such beam models by means of the finite element method. In his recent contribution [29], Romero points out the excellent performance of one specific category of beam elements denoted as geometrically exact beam formulations. While most of the geometrically exact beam formulations available in the literature are of Simo-Reissner type (see e.g. [3,9,11,28,30,31,32,33,39]), thus incorporating the modes of axial tension, shear, torsion and bending, in our recent contributions [22] and [23], a shear-free formulation based on the Kirchhoff theory of thin rods, thus being tailored for the modeling of slender fibers, has been proposed. The latter formulation will also be applied within this contribution since it provides a C 1 -continuous beam centerline representation, a desirable property enabling smooth contact kinematics in the context of beam-to-beam contact interaction.
Most of the applications mentioned above are characterized by mechanical contact interactions that significantly influence the overall system performance and by geometrically quite complex contact configurations, which allow for arbitrary fiber-to-fiber orientations. Despite the large number of publications concerning beam element formulations in general and despite the obvious need for robust and accurate beam contact formulations in many fields of application, there exists only a comparatively limited amount of literature focusing on beam-to-beam contact interaction. On the contrary, intensive research work has been done in the solid contact modeling of 3D continua within the last two decades. Important aspects in this field of research are for example the investigation of different constraint enforcement strategies (Lagrange multiplier method, penalty method, augmented Lagrange method etc.), types of contact discretization (node-to-segment/collocation-point-to-segment, Gauss-point-to-segment, mortar-like formulations), efficient contact search and active set strategies, procedures for Lagrange multiplier condensation and accurate integration schemes. Exemplarily, the reader is referred to the monographs [15,36] and to the review articles [35,27].
In contrast to contact formulations for 3D continua, which are typically based on a 2D contact traction field acting on the contact surfaces, the arguably most popular beam contact formulation [37] known in the literature models mechanical beam-to-beam contact interaction by means of a discrete contact force acting at the closest point between the two space curves representing the contacting beams (with circular cross-sections). This model, in the following denoted as point-to-point contact formulation, results in an elegant and efficient numerical formulation, which subsequently has been extended to frictional problems considering friction forces [38] and friction torques [13], rectangular beam cross-sections [20,21], smoothed centerline geometries [17], constraint enforcement via Lagrange multipliers [16] and adhesion effects [14]. Quite recently, it has been applied to self-contact problems [10]. However, one of the limitations of these point-to-point contact formulations can be attributed to the question whether contact between beams enclosing small contact angles, i.e. nearly parallel or entangled beams, should rather be modeled by means of a distributed line force instead of a discrete point force from a mechanical point of view. This question has been addressed by the recent publications [18] and [19], which propose additional contact points located in the neighborhood of the closest point in order to somewhat distribute the contact force in such configurations. Nevertheless, this formulation still relies on the existence of a locally unique closest point projection between the two contacting beams.
It is precisely this requirement that represents the second and essential limitation of point-to-point type beam contact formulations. In very general scenarios, such as in the applications mentioned in the beginning, where arbitrary beamto-beam orientations can occur, a unique closest point projection cannot be guaranteed for all potential contact regions. Consequently, some mechanically relevant contact points might be missed leading to large nonphysical penetrations or even to an entirely undetected crossing of the considered beams. There exist only a few alternative beam contact formulations available in the literature today that can overcome this limitation. One of these alternatives is the contact formulation developed by Durville [5], [6], [7], [8], which is based on a collocation-point-to-segment type formulation and the definition of proximity zones on an intermediate geometry. A second alternative proposed by Chamekh et al. [1], [2] is based on a Gauss-point-to-segment type formulation and primarily investigates self-contact problems of beams. What these two formulations have in common is that the contact forces are distributed along the two beams. Consequently, these types of formulations will be denoted as line-to-line contact formulations in the following.
The mentioned limitations of the point-to-point contact formulations were our motivation to perform mathematically concise and rigorous investigations concerning the existence of the corresponding closest-point-projection. In contrary to Konjukhov et al. [12], [13] who have already treated this question by means of geometrical criteria, we derive a very general analytical criterion that is valid for arbitrary contact configurations and that is based on proper and easyto-determine control quantities. Based on this analytical criterion, we can conclude that the standard point-to-point contact formulation is not applicable in a considerable range of practically relevant contact configurations. This result encouraged us to develop new improved line-to-line contact formulations, which are applicable beyond the standard point contact model, and which are inspired by some well-known and successful techniques known from contact mechanics for 3D solids. Specifically, we propose a novel beam contact formulation based on a Gauss-point-to-segment type contact discretization and a penalty regularization of the contact constraint. Additionally, we make theoretical considerations concerning alternative constraint enforcement strategies by means of Lagrange multipliers and alternative contact discretizations based on mortar methods. However, detailed theoretical and numerical investigations of these different approaches suggest the penalty-based Gauss-point-to-segment formulation as the variant that is most suitable for beam-to-beam contact and as the method of choice for the applications considered within this work.
In contrast to existing line-to-line beam contact formulations, our approach is extended by a consistently linearized integration interval segmentation that avoids numerical integration across strong discontinuities. It is verified by means of suitable numerical examples that precisely this component in combination with a smoothed contact force law and the applied C 1 -continuous beam element formulation leads to a drastic reduction of the numerical integration error. This, in turn, improves spatial convergence rates and in many cases only enables optimal convergence behavior under uniform mesh refinement. Furthermore, the resulting line-to-line contact algorithm is supplemented by contact contributions of the beam endpoints, whereas all existing contact formulations, no matter if point-to-point or line-to-line, typically search for minimal distance solutions only within the beams interior, but not for possible boundary minima. On the basis of a suitable numerical example, it is shown that in many applications these endpoint contact scenarios can appear with considerable frequency. Although, the influence of these endpoint forces on the overall solution quality might be of secondary interest, it is shown that neglecting these contributions will drastically reduce the robustness of the nonlinear solution scheme in many cases and may even prohibit convergence at all. Finally, a numerical test case suitable for line-to-line contact scenarios has been designed and a corresponding analytical solution based on the Kirchhoff theory of thin rods has been derived. This test case and the associated analytical solution can serve as valuable benchmark for the proposed formulation but also for future beam-to-beam contact approaches.
The remainder of this paper is organized as follows. In Section 2, we briefly repeat the main constituents of the applied beam element formulation initially proposed in [22,23] and extend the formulation to elastodynamics. In Section 3, the theory of standard point-to-point contact formulations is presented, followed by an analytical investigation of the existence and uniqueness of the required closest point projection and the derivation of a simple but mathematically concise criterion for solvability. In Section 4, the proposed line-to-line contact formulation is introduced and intensively compared to alternative methods known from the field of solid contact mechanics. The overall contact algorithm is completed by the contact contributions arising from the beam endpoints in Section 5 before a detailed numerical verification is performed in Section 6. While the first four examples presented in Sections 6.1 -6.4 aim at investigating the accuracy and consistency of the new formulation regarding integration error, spatial convergence behavior and contact force evolutions, the final two examples in Sections 6.5 and 6.6 represent possible real-world applications in order to verify the robustness of the proposed formulation when applied to practically relevant problems.

Applied beam formulation
In our recent contributions [22,23], a geometrically exact beam element formulation according to the geometrically nonlinear Kirchhoff theory of thin rods incorporating the modes of axial tension, torsion and non-isotropic bending has been proposed. The underlying beam theory and the resulting finite element formulation are tailored to deal with problems involving highly slender fibers and the numerical challenges (e.g. membrane locking) resulting from such high beam slenderness ratios. In addition to the general element formulation, a reduced element formulation neglecting the mode of torsion has been proposed in [23] and has been shown to deliver identical results as the general beam element formulation when restricting the considered structures to initially straight beams with circular crosssections and excluding axial/torsional moments from the set of external loads. Since these restrictions are easily fulfilled for the numerical examples considered in this contribution, we will exclusively resort to this simple and efficient torsion-free variant. However, the transfer of the following derivations from the torsion-free to the general element formulation is trivial, since the beam centerline representations of both are identical. The constituents of the static, torsion-free formulation presented in [23] will be summarized and extended to elastodynamics in the following.

Continuum formulation
The current configuration of the torsion-free beam is completely described by the beam centerline represented via a parametrized space curve (s, t) → r(s, t) ∈ 3 . Here, s ∈ [0, l] ⊂ and l ∈ represent an arc-length parametrization of the curve and the beam length in the initial configuration, respectively, and (.) = ∂ ∂s (.) = (.) ,s denotes the derivative with respect to this arc-length coordinate. Furthermore, t ∈ [0, t end ] ⊂ represents the time and( .) = ∂ ∂t (.) = (.) ,t denotes the corresponding time derivative. If we neglect rotational inertia contributions, which is common practice and mechanically sensible when considering highly slender beams, the extension of the torsion-free formulation according to [23] to dynamic problems is straightforward. In this case, the kinetic and hyper-elastic stored energies are: Here, ρ is the mass density, A the cross-section area, I the moment of inertia and E the Youngs modulus. Furthermore, v = v(s) represents the material velocity field, while = (s) and κ = κ(s) are the fields of axial tension and bending curvature. The corresponding weak form of the dynamic balance equations of the considered beam reads Here,f andm ⊥ denote distributed forces and moments, whereasf andm ⊥ denote discrete point forces and moments on the Neumann boundary Γ σ of the beam. Furthermore, we have applied the following additional abbreviations: ||r || 4 and δθ ⊥ = r × δr ||r || 2 .
As indicated by the subscript (.) ⊥ , the torsion-free beam theory is only applicable if the external moment vectors contain no components parallel to the centerline tangent vector, i.e. r T (s)m ⊥ (s) ≡ 0 ∀ s ∈ [0, l] and r Tm ⊥ ≡ 0 on Γ σ .

Spatial discretization
After having defined the weak form of the dynamic equilibrium equations, corresponding boundary and initial conditions and proper spaces of trial and test functions, i.e. r ∈ S ⊂ 3 satisfying the essential boundary conditions on the Dirichlet boundary Γ u and δr ∈ V ⊂ 3 with δr = 0 on Γ u , the space-and time-continuous problem setting is completed. Spatial discretization is performed by replacing the test and trial spaces by finite-dimensional subsets, i.e. r ≈ r h ∈ S h ⊂ S and δr ≈ δr h ∈ V h ⊂ V. Here and in the following, the index h denotes the spatially discretized version of a quantity. However, in the following, this index will often be omitted when there is no danger of confusion. Concretely, we follow a Bubnov-Galerkin approach leading to the following discretized beam centerline: whered i ,t i ∈ 3 are positions and tangent vectors at the two element nodes (i = 1, 2), δd i , δt i ∈ 3 represent their variations, l ele is the initial length of the initially straight beam element and ξ ∈ [−1; 1] is an element parameter coordinate. If considering initially straight beams, the latter can explicitly be related to the arc-length coordinate according to s(ξ) = s 0 + (ξ + 1)l ele /2 and (.) ,s = (.) ,ξ · J −1 ele (ξ). Here, s 0 represents the arc-length coordinate of the first node and J ele (ξ) = l ele /2 the element Jacobian. Similar to the abbreviation (.) = (.) ,s for the arc-length derivative, we will use the notation (.) = (.) ,ξ for the derivative with respect to the element parameter coordinate. The third order Hermite shape functions N i d (ξ) and N i t (ξ) (see [22] for further information concerning their properties) are defined as They provide a C 1 -continuous beam centerline representation, thus enabling smooth contact kinematics. This property will be very beneficial for the derivation of the contact formulation in the following sections. The abbreviations d, δd and N(ξ) appearing in (4) represent proper element-wise vector-and matrix-valued assemblies of the nodal variables and shape functions. In order to avoid membrane locking in the range of very high slenderness ratios as considered in this contribution, we additionally apply the so-called MCS method introduced in [23], where the original axial strain field and its variation δ occurring in (2) are replaced by the following re-interpolations which are based on second-order Lagrange polynomials L k (ξ). The resulting element residual contributions r int , r kin and r ext of the internal, inertia and external forces and their linearizations are summarized in Appendix A. An assembly of these quantities and the corresponding element-wise contact contributions r con presented in the next sections leads to the following global system of equations representing the spatially discretized version of (2), viz.
where D is the assembled global vector of primary variables containing the nodal degrees of freedomd k ,t k of all n node nodes with k = 1, ..., n node . It is worth to mention that the presented torsion-free beam formulation is a geometrically exact representation of a real "cross-section-reduced" structural model based on a 1D continuum theory, but it neither requires the application of any rotational primary degrees of freedom nor the enforcement of additional director constraints. Furthermore, the global inertia forces are composed of a constant symmetric mass matrix M and the global acceleration vectorD, while general geometrically exact beam formulations usually lead to nonlinear inertia force contributions. Finally, as long as no external moments are acting, i.e.m ⊥ (s) ≡ 0 ∀ s ∈ [0, l] andm ⊥ ≡ 0 on Γ σ , the global stiffness matrix K int = dR int /dD is symmetric (see Appendix A). Due to the absence of rotational degrees of freedom, any time discretization scheme suitable for second-order ODEs can be applied to (7).

Point-to-point contact formulation and limitations
Within this section, we briefly repeat the main constituents of a standard point-to-point beam contact formulation as introduced in [37]. Thereto, we consider two arbitrarily curved beams with cross-section radii R 1 and R 2 , respectively. The beam centerlines are represented by two parametrized curves r 1 (ξ) and r 2 (η) with curve parameters ξ and η. Furthermore, r 1,ξ (ξ) = r 1 (ξ) and r 2,η (η) = r 2 (η) denote the tangents to these curves at positions ξ and η, respectively.
In what follows, we assume that the considered space curves are at least C 1 −continuous, thus providing a unique tangent vector at every position ξ and η. The kinematic quantities introduced above are illustrated in Figure 1.

Contact formulation and contribution to weak form
The point-to-point beam contact formulation enforces the contact constraint by prohibiting penetration of the two beams at the closest point positions ξ c and η c . Here and in the following, the subscript c indicates, that a quantity is evaluated at the closest point coordinate ξ c or η c , respectively. These closest point coordinates are determined as solution of the bilateral ("bl") minimal distance problem, also denoted as bilateral closest point projection, with This leads to two orthogonality conditions that have to be solved for the unknown closest point coordinates ξ c and η c : The contact condition of non-penetration at the closest point is formulated by means of the inequality constraint where g is the gap function. This constraint can be included into our variational problem setting via a penalty potential or alternatively via a contact contribution in terms of a corresponding Lagrange multiplier potential Throughout this work, we solely apply constraint enforcement via penalty regularization according to (11) (see also our remarks in Section 4.4.3). Variation of (11) leads to the contribution of one contact point to the weak form: In (13), we can identify the contact force vector f cε as well as the normal vector n. The two are defined as: According to (14), the point-to-point beam contact formulation models the contact force f cε that is transferred between the two beams as a discrete point force acting at the respective closest points of the beam centerlines.
Here, the first two contributions denote the change in X(ξ, η) due to a change in the parameter coordinates ξ and η, whereas the contributions δX/∆X represent the variation/linearization of X(ξ, η) at fixed parameter coordinates. As already mentioned in [37], the total variation of the gap simplifies according to δg = n T (δ (r 1c ) − δ (r 2c )) = n T δr 1c + r 1c δξ − δr 2c − r 2c δη = n T (δr 1c − δr 2c ) since n T r 1c = n T r 2c = 0, which is a consequence of the orthogonality conditions (9) satisfied at the closest points ξ c and η c .
For later use, we also define the so-called contact angle as the angle between the tangent vectors at the contact point: In a next step, spatial discretization has to be performed. Since, for simplicity, we only consider the contact contribution of one contact point, the indices 1 and 2 are directly transferred to the two finite elements where the point contact takes place. Inserting the spatial discretization (4) into the orthogonality conditions (9) allows to solve the latter for the unknown closest point parameter coordinates ξ c and η c . Since, in general, the system of equations provided by (9) is nonlinear in ξ and η, a local Newton-Raphson scheme is applied for its solution. The corresponding linearizations of (9) can for example be found in [37]. Inserting equations (4) into equation (13) leads to the following contact residual contributions r con,1 and r con,2 of the two considered elements:

Limitations of point-to-point contact formulation
The point-to-point contact formulation provides an elegant and efficient contact model as long as sufficiently large contact angles are considered. However, its limitation lies in the requirement of a unique closest point solution according to (9), which cannot be guaranteed for arbitrary geometrical configurations. In [13], the authors have already treated the question of uniqueness and existence of the closest point projection by means of geometrical criteria based on so-called projection domains. Within this section, we want to analyze this question from a different perspective: This procedure will allow us to define easy-to-evaluate control quantities and to derive proper upper and lower bounds of these control quantities within which a unique closest point solution can be guaranteed in a mathematically rigorous manner. In the following, it will be derived that the contact angle α defined in (15), the closest point distance d bl as well as the geometrical (or mathematical) curvatureκ of the beam centerline according tō are such suitable control quantities. We have introduced the parameter coordinates ∈ [0;l] representing the arc-length of the current, deformed beam centerline andl denoting the corresponding current length. For the following analytical derivations, which are based on the space-continuous problem setting, we use the current arc-length parameterss 1 and s 2 instead of the initial arc-length parameters s 1 and s 2 (required for the space-continuous problem setting of the beam element formulation) or the normalized element parameters ξ and η (required for the spatially discretized problem setting). This choice simplifies many steps due to the essential property ||r 1,s 1 || = ||r 2,s 2 || = 1. Moreover, we define the maximal cross-section to curvature radius ratio µ max according to i.e. as the quotient of the cross-section radius R and the minimal radius of curvaturer occurring in the deformed geometry. The application of beam theories in general, particularly the application of the Kirchhoff beam theory, is only justified for problems exhibiting small values of this ratio, i.e. µ max 1. This property will be useful later on in this section. In order to simplify the following derivations, we anticipate the definition of the unilateral ("ul") distance function field d ul (s 1 ) presented in Section 4, which assigns a closest partner points 2c of the second beam (in this context also denoted as master beam) for every given points 1 on the first beam (in this context also denoted as slave beam) by means of the following unilateral closest point projection (see Figure 3(a) for an illustration): Next, one has to realize that the bilateral closest point projection (8) represents a special case of the unilateral closest point projection (19). Concretely, the closest point coordinates (8) are found through minimization of the minimal distance function d ul (s 1 ) according to (19) with respect to the slave beam parameters 1 , viz.: Now, in a first step, we want to examine the requirements for the existence of a unique solution of the unilateral closest point projection. As soon as we can guarantee a unique distance function d ul (s 1 ), the investigation of the existence and uniqueness of the bilateral closest point projection simplifies from the analysis of a function with 2D support occurring in (8) to the analysis of a function with 1D support according to (20). For a given point with coordinate vector r 1 (s 1 ), the unilateral closest point projection according to (19) searches for the corresponding closest point coordinates 2c on the space curve r 2 (s 2 ). In case of C 1 -continuous curves, which is guaranteed by the applied Hermite shape functions and which leads to a uniquely defined tangent vector field, a necessary condition for the existence of the minimal distance solution (19) is satisfied in case the requirement of a vanishing first derivative is fulfilled, i.e.
In case the two beams are close enough so that the sought-after closest points 2c is relevant in terms of active contact forces (g(s 1 ) = 0 → d ul (s 1 ) = 2R) and under consideration of the worst case cos(β 2 )(s 1 ) = 1, we obtain the following final requirement for a unique solution of the unilateral closest point projection according to (19): As a consequence of the maximal cross-section to curvature radius ratio µ max 1, a uniquely defined unilateral distance function d ul (s 1 ) can be guaranteed as long as the beams are sufficiently close. A corresponding criterion for arbitrary distances defined via d ul (s 1 ) =: k · R can be derived by replacing the factor 2 by k in (26). In a second step, we want to investigate the requirements for a unique bilateral closest point solution according to (20), based on a uniquely defined distance function d ul (s 1 ) (which is provided as consequence of (26)). Again, the first derivative ∂s 2c ∂s 1 = r T 1,s 1 (s 1 ) (r 1 (s 1 ) − r 2 (s 2c )) ||r 1 (s 1 ) − r 2 (s 2c )|| → r T 1,s 1 (s 1c ) (r 1 (s 1c ) − r 2 (s 2c ))=0 (27) has to vanish. This is satisfied at the closest points 1c by the first line of (9). Furthermore, the additional identity ∂d/∂s 2c ≡ 0 ∀s 1 ∈ [0;l 1 ] is fulfilled as consequence of the second line of (9). Again, a locally unique solution of the minimal distance problem (20) additionally requires a positive second derivative. Differentiation of (27) yields: The derivative ∂s 2c /∂s 1 appearing in (28) can be derived by consistently linearizing the orthogonality condition (21): After making use of this result and calculating the derivatives of (27) with respect tos 1 /s 2c , requirement (28) yields: Using the quantities defined in (24), the contact angle α according to (15) and the additional definitions condition (30) can be reformulated. Due to the strictly positive denominator, we only have to consider the numerator: where we have assumed sufficiently close beams d bl = k · R satisfyingκ 1 d bl < 1 andκ 2 d bl < 1 as consequence of (18). In case the two beams are close enough so that the sought-after closest point pair (s 1c ,s 2c ) is relevant in terms of active contact forces (g = 0 → d bl = 2R), the inequality (32) can be reformulated by means of worst case estimates: Since we solely consider positive contact angles α ∈ [0; 90 • ], only the positive branch of the quadratic inequality (33) has to be considered. Consequently, we end up with the following lower bound for the contact angle: The importance of the final requirement in (34) is quite obvious: As long as we can provide an upper bound µ max for the admissible ratio of cross-section to curvature radius, we will directly obtain from (34) a lower bound for the admissible contact angles above which the closest point solution is unique. Again, condition (34) can be expanded to general, but still sufficiently small (κd bl < 1!), distances d bl = k · R by replacing the factor 2 by k.
The three examples (a) Two parallel beams (b) Straight + circular beam (c) Straight + helical beam  Figure 2(a)), we obtain the trivial requirement α > 0, which reflects the well-known singularity of the closest point projection for parallel beams. If we consider a straight beam and a circular beam, both being oriented in a centrical manner as depicted in Figure 2(b), we observe a constant gap g(s 1 ) = const. between both beams, thus leading to a non-unique bilateral closest point solution, but this time at a contact angle of α = 90 • . However, this case is not practically relevant, since contact in such a scenario can only occur ifr ≈ 2R, therefore leading to a cross-section to curvature radius ratio µ ≈ 0.5, which is not supported by the considered beam theory, anyway. The third situation ( Figure 2(c)) is similar to the example that will later be numerically investigated in Section 6.2. The contact interaction between a straight beam and a helical beam again leads to a constant gap function g(s 1 ) = const. and consequently to a non-unique bilateral closest point solution. With decreasing slope h, the ratio of cross-section to curvature radius as well as the contact angle at which this non-unique solution appears increase. This is in perfect agreement with (34). In this context, the helix represents an intermediate configuration between the case of two straight parallel beams according to Figure 2(a) (slope h → ∞) and the case of a straight and a circular beam according to Figure  2(b) (slope h = 0). In Section 6.2 it will be shown that for such geometries a comparatively large scope of contact angles α ∈ [0 • ; α min ] can not be modeled by means of the standard point-to-point contact formulation. In practical simulations, the lower bound (34) has to be supplemented by a proper safety factor in order to guarantee for a unique closest point solution not only when contact actually occurs (g = 0) but already for a sufficient range of small positive gaps g > 0. Furthermore, too small angles α marginally above the lower bound (34) might lead to an ill-conditioned system of equations in (9) even if a unique analytical solution exists. Thus, the important result of this section is that the standard point-to-point contact formulation is not only unfeasible for examples including strictly parallel beams, but rather for a considerable range of small contact angles, since no locally unique closest point solution is existent in this range. According to (34), the size of this range depends on the ratio of the maximal bending curvature amplitude expected for the considered mechanical problem and the cross-section radius.
So far, we have only used mathematical arguments to show why the point-to-point beam contact formulation cannot be applied in the range of small contact angles. However, it is also questionable from a physical or mechanical point of view if the model of "point-to-point contact" itself is suitable to describe the contact interaction of beams enclosing small angles at all. On the one hand, it is clear that configurations providing a strictly constant distance function, i.e. d ul,s 1 (s 1 ) ≡ 0, are best modeled by a line-to-line and not by a point-to-point contact formulation. On the other hand, if an exact constraint enforcement of beams with rigid cross-sections is assumed, a pure point-to-point contact situation would already occur for non-constant distance functions with very small slopes, i.e. 0 < ||d ul,s 1 (s 1 )|| 1. However, this is a pure consequence of the rigid cross-section assumption inherent to the employed beam model, while a 3D continuum approach would naturally lead to distributed contact tractions. Consequently, also in the context of 1D continuum theories, such scenarios should better be modeled by a line-to-line rather than a point-to-point contact formulation. In the next section, a novel line-to-line contact formulation, which is capable of modeling arbitrary beam contact scenarios spanning the entire range of possible contact angles α ∈ [0 • ; 90 • ] and which is particularly beneficial for small contact angles and nearly constant distance functions d ul (s 1 ), will be proposed.

Line-to-line contact formulation
In the following, we present a novel line-to-line contact formulation that does not formulate the contact condition in form of a point-constraint at the closest points anymore, but rather as a line constraint enforced along the entire beam length. Consequently, we do not search for one closest point pair, but rather for a closest point field η c (ξ) on the second beam (master) assigned to the parameter coordinate field ξ on the first beam (slave). The relevant kinematic quantities of this approach are illustrated in Figure 3(a). The closest master point η c to a given slave point ξ is determined as solution of the following unilateral ("ul") minimal distance problem: It has already been shown in Section 3.2 (see (26)), that a unique unilateral closest point solution according to (35) can be guaranteed in case the considered beams are close enough so that contact can occur (i.e. d ul (ξ) ≈ 2R). Condition (35) leads to one orthogonality condition that has to be solved for the unknown parameter coordinate η c : Thus, in contrary to the procedure of the last section, the normal vector is still perpendicular to the second beam but not to the first beam anymore. Furthermore, in the context of line contact the subscript c indicates that a quantity is is formulated by means of an inequality-constraint for the gap function field g(ξ) along the entire slave beam.

Constraint enforcement and contact residual contribution
In the following, we apply a constraint enforcement strategy based on the space-continuous penalty potential: In Section 4.4, it will be shown that this strategy is preferable in beam-to-beam contact applications as compared to alternative methods known from contact modeling of 3D continua. The space-continuous penalty potential in (38) does not only serve as purely mathematical tool for constraint enforcement, but also has a physical interpretation: It can be regarded as a mechanical model for the flexibility of the surfaces and/or cross-sections of the contacting beams. Variation of the penalty potential defined in (38) leads to the following contact contribution to the weak form: In the virtual work expression (39), we can identify the contact force vector f cε (ξ) and the normal vector n(ξ): According to (40), the line-to-line beam contact formulation models the contact force f cε (ξ) that is transferred between the beams as a distributed line force. For comparison reasons, we can again define the contact angle field as: Next, spatial discretization has to be performed. For simplicity, we only consider the contact contribution stemming from one finite element on the slave beam and one finite element on the master beam being assigned to the former via projection (36). Therefore, the indices 1 of the slave beam and 2 of the master beam will in the following also be used in order to denote the two considered finite elements lying on these beams. Inserting the spatial discretization in (4) into the orthogonality condition (36) allows to solve the latter for the unknown closest point parameter coordinate η c (ξ) for any given slave coordinate ξ. The linearizations of (36) required for an iterative solution procedure can be found in Appendix C. Inserting the discretization (4) into equation (39) and replacing the analytical integral by a Gauss quadrature finally leads to the following contributions of element 1 and 2 to the discretized weak form: Here, n GP is the number of Gauss points per slave element, w k are the corresponding Gauss weights, ξ k are the Gauss point coordinates in the parameter space ξ ∈ [−1; 1] and finally η c,k is the closest master point coordinate assigned to the Gauss point coordinate ξ k on the slave beam (see also Figure 3(b)). The Jacobian J(ξ k ) maps between the slave beam arc-length increment ds 1 and an increment in the parameter space used for numerical integration (see also Section 4.2). Furthermore, r con,1 and r con,2 are the residual contributions of the slave (1) and master (2) element.

Remark:
In (39), we derived a similar expression for the variation of the gap as for the point-to-point contact case. This time, the variation δξ is zero since ξ remains fixed, and, again, the contribution due to the variation of η vanishes as a consequence of the orthogonality condition on the slave side: Remark: The gap function in (37) describes the exact value of the minimal beam surface-to-surface distance at a given coordinate ξ, only if the contact normal vector is perpendicular to both beam centerlines: While both conditions in (43) are exactly satisfied at the closest point of the point-to-point contact formulation per definition, only the second condition is fulfilled for an arbitrary contact point ξ within an active line-to-line contact segment. However, on the one hand, when considering non-constant evolutions of the centerline distance field along the considered beams, i.e. d ul (ξ) const., the region of active lineto-line contact contributions characterized by g(ξ) < 0, decreases with increasing penalty parameter. In the limit → ∞, the line-to-line contact formulation converges towards the point-to-point contact formulation, where both conditions (43) are fulfilled exactly. Thus, for a sensibly chosen penalty parameter, the gap function definition (37) provides also a good approximation for the line-to-line contact formulation. On the other hand, in configurations with constant centerline distance field d ul (ξ) = const., i.e. a range where no unique bilateral closest point solution exists and the point-to-point contact formulation cannot be applied, the two orthogonality conditions (43) are exactly fulfilled for the entire beam anyway.

Integration segments
From a pratical point of view, it is desirable to decouple the beam discretization and the contact discretization. This can be achieved by allowing for n II ≥ 1 contact integration intervals per slave beam element with n GR integration points defining a Gauss rule of order p = 2n GR − 1 on each of these integration intervals, thus leading to n GP = n II · n GR integration points per slave element. In order to realize such a procedure, one has to introduce n II further parameter spacesξ i ∈ [−1; 1] with i = 1, ..., n II on each slave element: In the simplest case, the parameter coordinates ξ 1,i and ξ 2,i confining the i th integration interval are chosen equidistantly within the slave element. Further information on the general determination of ξ 1,i and ξ 2,i is provided later on in this section. The total Jacobian J(ξ(ξ i )) = ds 1 /dξ i follows directly from (44) and reads where the mapping J ele (ξ(ξ i )) from the arc-length space s 1 to the element parameter space ξ results from the applied beam element formulation. Additionally, the sum over the number of Gauss points appearing in (42) has to be split: Here, the terms r i j con,1 and r i j con,2 denote the residual contributions of one individual Gauss point j in the integration interval i and the element parameter coordinates ξ i j are evaluated according to (44) at the Gauss point coordinatesξ j : Similar to the Gauss weights w j , these Gauss point coordinatesξ j are constant, i.e. not deformation-dependent, and identical for all integration intervals in case the same Gauss rule is applied in each of these intervals. The Gauss quadrature applied for integration of (46) guarantees for exact integration of polynomials up to order p = 2n GR − 1 when using an integration rule with n GR quadrature points per integration interval. However, by simply integrating across the element boundaries of two successive master elements associated with the considered integration interval via the closest point projection (36), the integrand would not have a closed-form polynomial representation anymore and the mentioned polynomial order of exact integration can not be guaranteed. On the one hand, the integrands occurring in (46) are not of purely polynomial nature, a fact, that precludes exact integration anyway. On the other hand, strong discontinuities in the integrand, such as e.g. jumps in the contact force from a finite value to zero at the master beam endpoints, might increase the integration error drastically. In the following, we try to find a compromise between integration accuracy and computational efficiency. Thereto, we subdivide the integration intervals introduced above into sub-segments whenever the projections of master beam endpoints lie within the considered integration interval. With this integration interval segmentation, we avoid integration across strong discontinuities at the master beam endpoints (see Figure 4(b)). However, we do not create integration segments at all master element boundaries (a) Subsegments at all master element boundaries (b) Subsegments only at master beam endpoints Figure 4: Creation of integration sub-segments on the slave beam in order to avoid discontinuities of the integrand (see Figure 4(a)), where weak discontinuities in the integrand might occur. A further example for locations showing weak discontinuities in the integrand are the boundaries of active contact zones, i.e. locations where the contact line force decreases from a positive value to zero. As we will see later, the integration across this kind of discontinuities is rather unproblematic due to the applied beam formulation being C 1 -continuous at the element boundaries (see Section 2) and an applied quadratic penalty law regularization (see Section 4.3) that leads to a smoother transition between contact and non-contact zones along the beam length. In order to find the boundary coordinate ξ B of an integration sub-segment created at a given master beam endpoint η EP , the latter has to be projected onto the slave beam according to the following rule (with p 2 according to (36)): where the given parameter coordinate η EP can take on the values −1.0 and 1.0 and ξ B is in general found via an iterative solution of (48). The derivative p 2,ξ needed for such an iterative solution procedure can be found in Appendix C.
In the worst, yet very unlikely, case that two master beam endpoints have valid projections according to (48) within one integration interval, this interval has to be subdivided into three sub-segments. In this case, for one of these three sub-segments, both boundary coordinates ξ 1,i and ξ 2,i are determined via (48) and are consequently deformationdependent. Thus, in general, the boundary coordinates ξ 1,i and ξ 2,i introduced in (44) can be determined by: if a valid master beam endpoint projection exists for i = 1, ..., n II , if a valid master beam endpoint projection exists . Thus, in the standard case, these boundary coordinates are equidistantly distributed and constant. In case a valid projection of a master beam endpoint onto an integration interval exists, ξ B1 (η EP , d 12 ) denotes the resulting deformation-dependent lower boundary of an created sub-segment, whereas ξ B2 (η EP , d 12 ) denotes the corresponding upper boundary. Equation (49) together with equations (47) and (45) provide all information necessary in order to evaluate the element residual contributions according to (46). The linearization of the contributions r i j con,1 and r i j con,2 of one individual Gauss point on element 1 can be formulated by means of the following total differential: The combination of a line-to-line type contact model with a consistently linearized integration interval segmentation at the beam end points as presented in this section, a quadratically regularized smooth penalty law and a C 1 -continuous smooth beam centerline representation is a distinctive feature of the proposed contact formulation. The benefits of these additional means are a drastical reduction of the integration error which enables a consistent spatial convergence behavior for a low number of Gauss points (see Section 6 for verification), an increase of the algorithmic robustness as well as a reduction of possible contact force/energy jumps without significantly increasing the computational effort.

Penalty Laws
Up to now, we have considered the following linear penalty law as introduced in (40) and illustrated in Figure 5(a): In practical simulations, one often applies regularized penalty laws that allow for a smooth contact force transition as illustrated in Figure 5(b). This second variant is favorable from a numerical point of view: First of all, it may improve the performance of tangent-based iterative solution schemes applied to the nonlinear system of equations stemming from the considered discretized problem, since a unique tangent exists at the transition pointḡ between the states of "contact" and "non-contact". Secondly, the time integration scheme applied in dynamic simulations benefits from such a smooth contact force law. And thirdly, also numerical integration of the line-to-line contact forces along the beam length (see Section 4.2) becomes more accurate if a smooth force law is used. The quadratically regularized (a) Standard linear penalty law (b) Linear penalty law with quadratic regularization Figure 5: Graphical visualization of standard and quadratically regularized penalty law penalty law applied within this contribution has the following analytical representation: For simplicity, all theoretical derivations within this work are still based on a linear penalty law according to (52). However, a more general form of these equations that is valid for arbitrary penalty laws can easily be derived by simply replacing all linear force-like expressions of the form −ε g by the generic expression f cε (g).

Alternative constraint enforcement strategies
Similar to point-to-point contact formulations, the constraint equation resulting from the line-to-line contact formulation can be considered within a variational framework by means of a Lagrange multiplier potential or by means of a penalty potential. In contrast to the point-to-point case, however, the constraint in (37) is not only defined at a single point but rather on a parameter interval ξ ∈ [ξ a , ξ b ]. According to Section 4.1, the penalty method, which introduces no additional degrees of freedom, can be directly applied in terms of a space-continuous penalty potential, see (38), that can alternatively be interpreted as a simple hyper-elastic stored-energy function representing the accumulated cross-section stiffness of the contacting beams. The final contact formulation resulting from such a procedure after spatial discretization and numerical integration is often denoted as Gauss-point-to-segment type formulation.
In contrary, the Lagrange multiplier method applied to the constraint in (37) introduces an additional primary variable field λ(ξ), which is typically discretized in a manner consistent to the spatial discretization of the displacement variables (discrete inf-sup stable pairing). Eventually, the nodal primary variables resulting from the discretization of the Lagrange multiplier field can be considered as additional unknowns or be eliminated by means of a penalty regularization (applied to a spatially discretized version of (37)). Both variants are typically denoted as mortar-type formulations (see e.g. [25], [26]). In Section 4.4.1, the main steps of applying a mortar formulation to beam contact problems, thus representing an alternative to the formulation of Section 4.1, are provided. Finally, in Sections 4.4.2 and 4.4.3, a detailed comparison and evaluation of the variants "Gauss-point-to-segment" versus "mortar" and "penalty method" versus "Lagrange multiplier method", respectively, is performed in the context of beam-to-beam contact.

Constraint enforcement based on consistent Lagrange multiplier discretization
As an alternative to Section 4.1, we now consider constraint enforcement via a Lagrange multiplier potential: Variation of the Lagrange multiplier potential leads to the following contact contribution to the weak form: In (55), the contact force f cλ (ξ) = −λ(ξ)n(ξ) =: f cλ (ξ)n(ξ) transferred between the two beams can again be interpreted as a distributed line force. This time, the Lagrange multiplier field represents the magnitude of this line force. Next, spatial discretization has to be performed. Again, we consider the contribution of one slave beam element 1 and one master beam element 2. In addition to the spatial discretization (4), a trial space λ ≈ λ h ∈ S λh ⊂ S λ ⊂ and a weighting space δλ ≈ δλ h ∈ V λh ⊂ V λ ⊂ have to be defined for the field of Lagrange multipliers, too: Here, n λ represents the number of nodes of the Lagrange multiplier discretization per slave element, the vector N λ,1 (ξ) collects the corresponding test and trial functions with support on slave beam 1, andλ 1 as well as δλ 1 contain the corresponding discrete nodal Lagrange multipliers and their variations, respectively (see e.g. [34] concerning a proper choice of the spaces S λh and V λh ). Inserting (4) and (56) into (55) and replacing the analytical integral by a Gauss quadrature finally leads to the following contribution of elements 1 and 2 to the discretized weak form: Again, r con,1 and r con,2 represent the contact force residual contributions of slave element 1 and master element 2, whereas r λ,1,2 denotes the corresponding residual contribution stemming from constraint equation (37). Based on (57), different strategies of constraint enforcement are possible: Considering the nodal Lagrange multipliersλ 1 as additional unknowns would lead to an exact satisfaction of the discrete version of the constraints (37). Alternatively, these discrete constraint equations can be regularized by means of a penalty approach. Let n ele,s denote the total number of slave elements. Then, one typically defines so-called nodal gapsĝ j according tô In (58), a summation over all slave elements with support of the shape function N j λ,1 (ξ) assigned to the nodal gapĝ j is sufficient. Consequently, each nodal gap according to (58) represents one line of the total residual contribution R λ resulting from constraint equation (37). Now, one can replace the nodal Lagrange multipliers by nodal penalty forces: Inserting the nodal penalty forces instead of the unknown nodal Lagrange multipliers into (57) finally results in: This procedure eliminates the additional nodal unknownsλ. However, the constraint of vanishing nodal gapsĝ j will not be exactly fulfilled anymore. The only difference of the discretized weak form (42), i.e. the one resulting from a space-continuous penalty potential, and (60), i.e the one resulting from a discretized Lagrange multiplier potential and a subsequent penalty regularization, lies in the definition of the scalar contact forces λ ε (ξ) and ε g(ξ) , respectively.

Comparison of the two penalty approaches
The main advantage of the formulation presented in Section 4.4.1 is that it results from a consistent Lagrange multiplier discretization. As long as the trial and weighting spaces S h , V h , S λh and V λh are chosen such that a proper discrete inf-sup-stability condition is satisfied, no contact-related locking effects have to be expected, even for large values of the penalty parameter. This does in general not hold for the formulation presented in Section 4.1, where contact-related locking might occur for very high penalty parameters. When considering highly slender beams, moderate values of the penalty parameter are often sufficient in order to satisfy the contact constraint with the desired accuracy. In Section 6, it will be verified numerically that within this range of penalty parameters the spatial convergence behavior is not deteriorated by contact-related locking effects when applying the contact formulation according to Section 4.1. A crucial advantage of the latter formulation lies in its efficiency and its straight-forward implementation. On the one hand, the numerical implementation of the variant presented in Section 4.4.1 requires an additional element evaluation loop in order to determine the nodal gaps according to (58), or, in other words, the penalty-based elimination of the Lagrange multipliers cannot exclusively be conducted on element level. On the other hand, in combination with the standard gap function definition according to (37), this variant requires a very fine finite element discretization when applied to contact problems involving highly slender beams. This fact will be illustrated in the fol-  Figure 6(a), two straight beam elements with cross-section radii R 1 = R 2 = R characterized by a comparatively small contact angle and a large penetration of almost g(ξ c ) ≈ 2R are depicted. The resulting contact line force vector field according to (40) is illustrated in green color. Furthermore, in Figure 6(b), the evolution of the gap function is plotted over the length of the slave beam element. With increasing penalty parameter, the formulation according to Section 4.4.1 forces the nodal gaps in (58) to vanish. Roughly speaking, this means that the areas enclosed by positive gaps and the areas enclosed by negative gaps, as indicated with red and blue color in Figure 6(b), must balance each other. For small contact angles and reasonable spatial discretizations, this is possible. However, when looking at the gap-function evolution resulting from two almost perpendicular beams as illustrated in Figure 7(b), such a balancing can only be achieved if the beam element length is reduced drastically. This need for a sufficiently fine spatial discretization increases the numerical effort of this method. Alternatively, one might modify the definition of the gap function g(ξ), such that negative/positive gap contributions are weighted stronger/weaker. Since such an extra effort is not necessary for the procedure proposed in Section 4.1, we want to focus on this variant.

Penalty method vs. Lagrange multiplier method
Constraint enforcement by means of Lagrange multipliers is common practice in the field of computational contact mechanics for solids, especially in combination with mortar methods (see Section 4.4.1), due to some advantageous properties, for example concerning the accuracy of contact resolution. Even though the application of the Lagrange multiplier method for constraint enforcement in beam-to-beam contact scenarios has already been investigated in [16], the vast majority of publications in this field is based on regularized constraint enforcement via the penalty method. This fact can be justified by a couple of reasons: When considering discretizations based on structural models the ratio of surface degrees of freedom to all degrees of freedom (=1 for beams) is much larger than for solid discretizations based on a 3D continuum theory. Consequently, also the ratio of additional Lagrange multiplier degrees of freedom to displacement degrees of freedom would be comparatively high when enforcing, e.g., beam-to-beam line contact constraints (see Section 4.4.1) by means of Lagrange multipliers. Furthermore, when modeling slender structures by means of mechanical beam models, which are often based on the assumption of rigid cross-sections, computational efficiency is one of the key aspects whereas the resolution of exact contact pressure distributions and other mechanical effects on the length scale of the cross-section, which is typically by orders of magnitude smaller than the length dimension of the beam, is not of primary interest. If one is primarily interested in the global system behavior, even penetrations on the order of magnitude of the cross-section radius are often tolerable. Typically, penalty parameters required to limit the penetrations to such values decrease with the beam thickness. Often, the required values are proportional to the beam bending stiffness and therefore the penalty contributions do not significantly deteriorate the conditioning of the system matrix which is usually dominated by high axial and shear stiffness terms.
Besides the arguments above, there is one further crucial point, which makes the penalty method not only preferable to constraint enforcement via Lagrange multipliers, but which even prohibits the use of the latter method. Many of the perhaps most efficient and elegant beam models available in the literature (see e.g. the comparison of ANS beams and geometrically exact beams in [29]), are based on the assumption of rigid cross-sections. Especially when considering very thin beams, this assumption is well-justified and the properties of the resulting beam formulations are desirable from a numerical point of view. However, combining the assumption of rigid cross-sections and contact constraint enforcement via Lagrange multipliers leads to the following dilemma when considering, e.g., the dynamic collision of two beams: In the range of large contact angles, the initial kinetic energy will be transformed into elastic bending energy and back to kinetic energy during the impact. However, with decreasing contact angle the elastic bending deformation decreases and in the limit of two matching, exactly parallel beams the amount of elastic deformation during the collision drops to zero, since the cross-sections are rigid. The accelerations and contact forces resulting from such a scenario are unbounded and the resulting numerical problem become singular. Thus, undoubtedly, a certain amount of cross-section flexibility is indispensable when modeling such a scenario. This cross-section flexibility can be provided by a penalty force law such as the one in Section 4.1, which already has the structure of a typical hyper-elastic strain energy function and models the accumulated stiffness of the cross-sections of the two contacting beams. Of course, this idea can be refined by deriving more sophisticated penalty laws in form of reduced models based on a continuum mechanical analysis of the cross-section deformation and stiffness. However, since our primary intention is the regularization of parallel-impact scenarios and not the resolution of local deformations on the cross-section scale, we will keep the simple and convenient force law according to (53) in the following. Nevertheless, the adaption of the presented theory to more general penalty laws is straightforward. Furthermore, with these considerations in mind, the penalty parameter in the context of rigid-cross-section beam contact is no longer a pure mathematical tool of constraint enforcement, but it rather has a physical meaning: it serves as mechanical model of the beam cross-section stiffness. This interpretation simplifies the determination of a proper penalty parameter.

Endpoint-to-line and endpoint-to-endpoint contact contributions
(a) Element parameter space (b) Endpoint-to-line (c) Endpoint-to-endpoint The contact formulations presented in the last two sections have only considered solutions of the minimal distance problem within the element parameter domain ξ, η ∈ [−1; 1] as represented by the blue area in Figure 8(a). Due to the C 1 -continuity of our discrete centerline representation, also solutions coinciding with the element nodes are found by this procedure. However, a minimal distance solution can also occur in form of a boundary minimum at the physical endpoints of the contacting beams. The boundary solutions indicated by the four red lines in Figure 8(a) represent solutions with one parameter taking on the value −1 or 1 and the other parameter being arbitary. Mechanically, these solutions can be interpreted as the minimal distance appearing between a physical beam endpoint and an arbitrary beam segment as indicated in Figure 8(b). Additionally, a minimum can also occur in form of the distance between the physical endpoints of both beams (see Figure 8(c)), which corresponds to the four green corner points in Figure 8(a). Neglecting these boundary minima can lead to impermissibly large penetrations and even to an entirely undetected crossing of the beams. At first view, these contact configurations seem to be comparatively rare for thin beams and the mechanical influence of these contact contributions seems to be limited. However, practical simulations have shown that neglecting these contributions does not only lead to a slight inconsistency of the mechanical model itself but also to a drastically reduced robustness of the nonlinear solution scheme, since initially undetected large penetrations can lead to considerable jumps in the contact forces during the iterations of a nonlinear solution scheme. While for the endpoint-to-endpoint case, the contact point coordinates are already given, the endpoint-to-line case requires a unilateral closest-point-projection similar to the one in (35). Depending on which beams endpoint is given, this unilateral closest-point-projection either searches for the closest point η c to a given point ξ ∈ {−1, 1} or for the closest point ξ c to a given point η ∈ {−1, 1}. As soon as the contact point coordinates are known, one can directly apply the residual contribution of the point-to-point contact formulation according to (16). From a geometrical point of view, applying this model means that the beam endpoints are approximated by hemispherical surfaces. Again, it is justified to only consider the variation contribution with fixed ξ and fixed η for δg according to (13), since either the considered parameter coordinate is indeed fixed (if representing a physical endpoint) or the corresponding tangent vector is perpendicular to the contact normal (if representing the projection onto a segment). Nevertheless, one has to distinguish between the cases endpoint-to-endpoint and endpoint-to-line contact in order to correctly include the increments ∆ξ and ∆η in the linearizations of the contact residuals (see Appendix B for details).

Numerical examples
In this section, we want to verify the robustness and accuracy of our new line-to-line contact formulation presented in Section 4. For all examples, a standard Newton-Raphson scheme is applied in order to solve the nonlinear system of equations R tot resulting from the discretized weak form (7). As convergence criteria we check the Euclidean norms of the displacement increment vector ∆D k and of the residual vector R k tot at Newton iteration k. For convergence, these norms have to fall below prescribed tolerances δ R and δ D , i.e. ||R k tot || < δ R and ||∆D k || < δ D . If nothing to the contrary is mentioned, these tolerances are chosen according to the following standard values δ R = δ D = 1.0 · 10 −7 . The first example is a simple static patch test that should verify the effectiveness of the integration interval segmentation introduced in Section 4.2. As illustrated in Figure 9, the example consists of one completely fixed, rigid beam discretized with three beam elements (different element lengths) and a second, deformable beam discretized by two beam elements with cross-section radii R 1 = R 2 = 0.005, Youngs moduli E 1 = E 2 = 1.0 · 10 9 , length of the first beam l 1 = 2.0 and length of the second beam l 2 = 0.8. The second beam is loaded by a constant transverse line load p = 1.0 and its left endpoint is exposed to a Dirichlet-displacement of ∆u = 1.001 within 100 equidistant load steps. Furthermore, contact interaction between the two beams is modeled by the linear penalty law according to (52) with a penalty parameter ε = 500. As a consequence of the constant transverse line load and the chosen penalty parameter, there exists a trivial analytical solution with a constant gap g re f = −p/ε = −0.002 along the entire upper beam. In order to verify the working principle of the integration interval segmentation in the presence of strong discontinuities, we have chosen the first (rigid) beam as slave beam. In Figure 10, the average relative error e rel = n GP,tot i=1 g i − g re f n GP,tot · g re f of the gaps g i at the active Gauss points is plotted over the number of load steps for the formulations with and without integration interval segmentation at the beam endpoints in combination with different numbers of Gauss points n GP,tot . In all cases, three integration intervals per slave element have been applied. From Figure 10(a), one observes that the strong discontinuity of the contact force ε g(ξ i j ) occurring in the integrand of (46) leads to a considerable integration error that only gradually decreases when increasing the number of Gauss points. As expected, the formulation with integration interval segmentation (see Figure 10(b)) yields a significantly lower integration error level and a faster decline in the error with increasing number of Gauss points. Yet, even this formulation does not allow for an exact integration, in general, since the test functions N 1 and N 2 in (46) have no closed-form polynomial representation across the element boundaries. However, it will be shown in the next examples that the corresponding integration error is typically lower than the overall discretization error and therefore of no practical relevance. Furthermore, compared to a formulation with integration interval segmentation at all master beam element nodes, which would then allow for exact numerical integration, the proposed segmentation strategy is considerably less computationally expensive.

Example 2: Twisting of two beams
The second example aims at verifying the accuracy and consistency of the line-to-line contact formulation by investigating the spatial convergence behavior. Thereto, we consider two initially straight and parallel beams with circular cross-sections and radii R 1 = R 2 = R = 0.01, initial lengths l 1 = l 2 = l = 5 and Youngs moduli E 1 = E 2 = E = 1.0·10 9 as illustrated in Figure 11(a). The initial geometries of the two beams k = 1, 2 are given by the analytical expression: The distance of the two beams is chosen such that they exhibit an initial gap of g 0 = −0.1R. We clamp the beams at one end and move the beam cross-sections in a Dirichlet-controlled manner at the other end such that the corresponding cross-section center points move on a circular path. By this procedure, the two beams get twisted into a double-helical shape as illustrated in Figure 11(a). We try to adapt the system parameters in a way such that the analytical solution for the deformed beams is exactly represented by a helix with constant slope according to In the following, we only present the corresponding results, while the derivation based on the projected ODEs representing the strong form of the Kirchhoff theory (see [23]) is summarized in Appendix D. Before the actual twisting process starts, the two beams are pre-stressed by an axial displacement at the left endpoints (superscript "l") within one load step. Then, these points are moved on a circular path with radius r = R − |g 0 |/2 = 0.0095, i.e. within n l = 8 further load steps in order to end up with one full twist rotation. The translational displacements at the right endpoints of the right beams (superscript "r") are set to zero, i.e.
Furthermore, the x−components of all tangential degrees of freedom (see also Section 2) are set to zero, i.e.
whereas the y− and z−components of these nodal tangents are not prescribed but part of the numerical solution. As shown in Appendix D, these boundary conditions for the tangential degrees of freedom are sufficient in order to impose the necessary boundary moments at the endpoints. If, finally, the penalty parameter is chosen according to the resulting analytical solution obeys the analytical representation of (62), thus showing a gap of g 0 between the two beams that is constant along the beam lengths. As already mentioned earlier, the penalty parameter and the resulting gap between the two beams occurring in the analytical solution (62) can be interpreted as a mechanical model for the contact-surface/cross-section flexibility of the considered beams. Furthermore, the derived analytical solution corresponds to a mechanical state consisting of constant axial tension , constant bending curvature κ = (1+ )r r 2 +h 2 and vanishing torsion along both beams. In Figure 11(b), the relative L 2 -error of the FE solution for beam 1 is plotted with respect to the analytical solution over the element length for discretizations with 4, 8, 16, 32, 64, 128 and 256 elements per beam. For all convergence plots in this work, the following definition of the relative L 2 -error has been applied: Herein, r h denotes the numerical solution of the beam centerline position for a certain discretization. For all examples without analytical solution, the standard choice for the reference solution r re f is a numerical solution using a spatial discretization that is by a factor of four finer than the finest discretization shown in the corresponding convergence plot. The normalization with the element length l makes the error independent of the length of the considered beam. The second normalization leads to a more convenient relative error measure, which relates the L 2 -error to the maximal displacement u max occurring for the investigated load case.
In order to investigate the influence of the applied Gauss rule, we compare the cases of a 5-point and a 2-point Gauss rule with one integration interval per element in both cases. According to Figure 11(b), the 5-point-variant converges towards the analytical solution up to machine precision with the optimal order O(h 4 ) as expected for the applied thirdorder beam elements. Throughout this work, this 5-point-rule will be the default value if nothing to the contrary is mentioned. Reducing the number of Gauss integration points to a value of 2 leads to slight increase of the L 2 -error in the range of comparatively rough spatial discretizations. However, for finer discretizations the 2-point curve converges towards the 5-point curve. When looking at the upper right data point in Figure 11(b), one observes the remarkable result that a total of 8 contact evaluation points per beam (4 elements per beam with 2 Gauss points per element) is sufficient in order to end up with a relative error that is far below 1%.
In Section 3.2, we have derived a lower bound α min for the contact angle, above which a unique bilateral closest point projection exists. In the following, we briefly want to verify the corresponding result (34) by means of a slightly modified version of the considered twisting example. Thereto, we assume that the maximal admissible ratio of crosssection to curvature radius supported by the beam theory is 1%, i.e. µ max = 0.01. For simplicity, we additionally assume that the helix radius given in (62) equals the beam cross-section radius, i.e. r = R and consequently g 0 = 0. With µ max = 0.01, the minimal admissible slope for a helix with constant slope similar to (62) can be calculated as: Furthermore, after some geometrical considerations, one can calculate for the case h = h min the actual contact angle enclosed by two corresponding tangents, which is a constant angle in case of helical beams similar to (62): r T 1,ϕ (ϕ)r 2,ϕ (ϕ) ||r 1,ϕ (ϕ)|| · ||r 2,ϕ (ϕ)|| This is exactly the same result that we would obtain for the lower bound α min by inserting µ max = 0.01 into (34). This means that the helix geometry according to example 2 represents an extreme case, where all worst-case assumptions made in the derivation (33) become true and where, for a given admissible radius ratio µ max = 0.01, a non-unique closest point solution appears exactly at the contact angle α min predicted as lower bound by equation (34). On the other hand, this example shows that (34) provides the best possible lower bound, since it actually occurs in a practical example. Furthermore, it can be concluded that the considered twisting example, leading to a constant gap function along both beams, can of course not be modeled by means of the standard point-to-point contact formulation.

Example 3: General contact of two beams
So far, we have only considered scenarios with a constant gap function along the beam length. By means of the following examples, the more general case of non-constant gaps, and especially the case of a change in sign in the gap evolution along the beam, will be investigated. At positions with a change in sign in the gap function, the contact force according to the standard law in (52) drops to zero. As illustrated in Figure 6(c), this leads to a kink in the force evolution at this point, which becomes more and more pronounced with increasing contact angle (see Figure 7(c)). This weak discontinuity in the integrand may in general increase the numerical integration error and can be avoided by replacing the standard linear force law by the smoothed force law in (53)  influence of these two different force laws on the integration error and eventually on the spatial convergence behavior will be investigated by means of the following example: We consider beam geometries and material parameters identical to the last example. The penalty parameter is decreased to ε = 1000. Also, the initial configuration is similar to the one illustrated in Figure 11(a) of the last example. However, this time the initial distance between the beams is increased to a value of 2r = 4R = 0.04. The Dirichlet boundary conditions of the tangential degrees of freedom are slightly changed in order to completely avoid any cross-section rotation at the boundaries. Correspondingly we have: ∆t l 1,x = ∆t l 2,x = ∆t r 1,x = ∆t r 2,x = ∆t l 1,y = ∆t l 2,y = ∆t r 1,y = ∆t r 2,y = 0.
Thus, this time the tangents are completely clamped at both ends. Furthermore, no axial pre-stressing is applied, i.e.
The remaining Dirichlet conditions are similar to the last example, see (64) and (65). The resulting deformed configuration is illustrated in Figure 12(a). Due to the larger separation of the beams, the gap function increases from negative values to positive values when approaching the beam endpoints. The corresponding contact force evolutions resulting from different spatial discretizations are illustrated in Figure 13(a). In Figures 12(b) and 12(c), the relative L 2 -error with respect to a numerical reference solution is plotted for the formulation based on a linear penalty law and the formulation based on the quadratically regularized force law (regularization parameterḡ = 0.1R = 0.001). In case of the simple linear penalty law, the number of Gauss points has to be enhanced by a factor of 10 as compared to the standard 5-point rule in order to ensure O(h 4 ) convergence within the considered range of spatial discretizations (see Figure 12(b)). Thus, obviously, the increased integration error resulting from the kink in the penalty force law dominates the spatial discretization error if the standard 5-point Gauss rule is applied. Only an increase in the number of Gauss points, and therefore an increase in the numerical effort, reduces this integration error. An elimination of this kink by means of a smoothed penalty law enables the same accuracy and the optimal convergence order O(h 4 ) already with the standard 5-point Gauss rule (see Figure 12(c)) and consequently reduces the numerical effort drastically.

Example 4: Influence of integration interval segmentation on convergence behavior
In the first example, we have already illustrated how the integration error can be reduced by means of an integration interval segmentation at the beam endpoints. Now, we want to investigate the influence of this method on the spatial convergence behavior. Again, we consider beam geometry, material parameters as well as the penalty parameter to be identical to the last example. In order to enforce an integration across the beam endpoints, the initial geometry of one of the beams is shifted by a value of r = 2R = 0.02 along the positive z-axis leading to the representation: For this example, we apply the following Dirichlet boundary conditions at the endpoints of the two considered beams: ∆d r 1,x = −0.12, ∆d r 2,x = 0.12, ∆d r 1,y = ∆d r 2,y = ∆d r 1,z = ∆d r 2,z = 0, ∆d l 1,x = −r 1−cos k2π n l , ∆d l 1,y = r sin k2π n l , ∆d l 2,x = r 1−cos k2π n l , ∆d l 2,y = −r sin k2π n l , k = 1, ..., n l , ∆d l 1,z = ∆d l 2,z = 0, ∆t l 1,x = ∆t l 2,x = ∆t r 1,x = ∆t r 2,x = ∆t l 1,y = ∆t l 2,y = ∆t r 1,y = ∆t r 2,y = 0.
The two additional conditions in the first line of (74) enforce a negative gap and consequently active contact forces at the (non-matching) right endpoints of the beams. By this means, we enforce an integration across a contact force jump at these endpoints which is sensible in order to investigate the effectiveness of the integration interval segmentation. All the remaining Dirichlet conditions appearing in (74) are similar to Section 6.3. The deformed geometry resulting from these boundary conditions is illustrated in Figure 14(a). Furthermore, the contact force evolutions corresponding to different finite element meshes are presented in Figure 13(b). The contact force evolution shows the expected jump from f c (s = 0.02 − ) = 0 to f c (s = 0.02 + ) ≈ 4.5 at position s = 0.02 (see also the detail view in Figure 13 interval segmentation (Figure 14(b)) and once without a corresponding segmentation (Figure 14(c)).
According to Figure 14(b), the remaining integration error of the formulation with interval segmentation and a 5-point Gauss rule slightly deteriorates the spatial convergence behavior. However, by applying two instead of one 5-point Gauss integration intervals per element, this influence of the integration error vanishes and we observe the optimal convergence order O(h 4 ). On the contrary, the convergence behavior of the formulation without integration interval segmentation (see Figure 14(c)) is still deteriorated by the integration error for a 5-point Gauss rule even with two intervals per element. Even if the number of intervals is increased to 10, i.e. an increase of the number of Gauss points by a factor of 5, this negative influence is still visible in the range of fine discretizations. Furthermore, it is worth mentioning that this effect is expected to become even more pronounced in practical applications, where the displacements are not Dirichlet-controlled in the direct neighborhood of the strong discontinuity. All in all, it seems that the integration interval segmentation solely applied at the beam endpoints represents a sensible compromise of integration accuracy and computational efficiency. Additionally, in dynamic simulations, this strategy prevents from force and energy jumps in scenarios where active Gauss points of standard integration schemes based on fixed, nonsegmented integration intervals would slide across master beam endpoints.

Example 5: Simulation of a biopolymer network
In a first practically relevant example, we apply the presented simulation framework in order to investigate the influence of mechanical contact interaction on the three-dimensional Brownian motion of filaments in biopolymer networks. Biopolymer networks are tight meshes of highly slender polymer filaments (e.g. Actin filaments) embedded in a liquid phase, often interconnected by means of a second molecule species (so-called cross-linkers). These networks can for example be found in biological cells. There, they crucially determine the mechanical properties of cells and biologically highly relevant processes such as cell-migration or cell-division. In a recent contribution [4], a finite element model for the Brownian motion of these filaments in the absence of mechanical contact interaction has been proposed. Accordingly, these slender filaments are described by means of a geometrically nonlinear beam theory. The mechanical interaction of the filaments with the surrounding fluid is modeled by means of external distributed line loads consisting of velocity-proportional viscous drag forces and thermal excitation forces. The latter are modeled as stochastic forces characterized by a mean value of zero, a variance determined by the absolute temperature and the properties of the surrounding fluid and finally by spatial and temporal correlation lengths which are assumed to be zero (see [4] for details). Due to the physical length scales relevant for such systems, inertia forces can often be neglected, thus leading to a system of first-order stochastic partial differential equations (SPDEs). While spatial discretization is performed via the finite element method in combination with a geometrically exact beam formulation (Section 2), a Backward Euler scheme is applied in order to discretize the resulting semi-discrete problem in time.
Here, we combine the Brownian dynamics simulation framework presented in [4] with the contact algorithm proposed in the sections before in order to simulate the free diffusion of Actin filaments (without consideration of cross-linker . Thereto, we consider a system of 37 initially straight and randomly distributed Actin filaments with circular cross-section of radius R = 2.45 · 10 −3 , length l = 2 and Youngs modulus E = 1.3 · 10 9 (all quantities given in the units mg, µm, and s) as illustrated in Figure 15(a). All further physical system parameters describing the viscous and stochastic forces are identical to those applied in [24] and can be found therein.
The simulation was performed based on a spatial discretization with 8 beam elements per filament, a time step size of ∆t = 1.0·10 −4 and a total simulation time of t end = 1.0·10 −1 . Furthermore, the contact parameters have been chosen as ε = 5.0 · 10 4 andḡ = 2.0 · 10 −3 in combination with 50 integration intervals per slave element based on a 5-point Gauss rule, respectively. The spatial configurations at times t = 0.0 and t end = 0.1 as well as a corresponding detail view at t end = 0.1 are illustrated in Figure 15. Due to the stochastic forces, the velocity field of these filaments is strongly fluctuating in space and in time, thus leading to drastic and frequent changes in the active contact sets. This property in combination with the high filament slenderness ratio of approximately 800 makes this example very demanding concerning the robustness of the proposed contact algorithm. The Newton-Raphson convergence tolerances are set to δ R = δ D = 10 −6 . In this example, where dynamic collisions at all possible filament-to-filament orientations can occur, the significance of the endpoint contact contributions introduced in Section 5 becomes apparent. In order to underpin this statement, the corresponding total numbers of active line contact Gauss points and active beam endpoint contacts have been plotted over the simulation time in Figures 16(a) and 16(b). Accordingly, even for this comparatively small example, the endpoint contact contributions occur with significant frequency. Neglecting these endpoint contact forces would not only allow for nonphysically large penetrations, it would also lead to non-convergence of the Newton-Raphson scheme in many time steps. The contact angles measured for this example during the simulation time lie within the range α ∈ [4 • ; 90 • ], thus covering almost the entire possible scope. Considering the maximal curvaturē κ max ≈ 2.0 measured during the simulation time and the beam cross-section radius R = 2.45 · 10 −3 , the lower bound for the contact angle that would allow for a point-to-point contact formulation can be calculated as α min ≈ 8 • . This means that even this example, which is dominated by rather large contact angle configurations, cannot be completely covered by a standard point-to-point contact formulation. Finally, we conclude this section by an exemplary statistical analysis of a physically relevant quantity. Concretely, the influence of the mechanical contact interaction on the filament diffusion measured by the mean square displacement per time step is evaluated. This quantity is defined as: In (75) Statistical evaluation of the measured modeling errors finally yield a mean value of 8.5% and a variance of 2.5% (see Figure 16(c)). In other words, for the considered example, the mean square displacement per time step is overestimated by 8.5% in average when neglecting mechanical contact interaction. Of course, this analysis only has an exemplary character, since system parameters such as fluid and filament properties, magnitude of stochastic forces, considered simulation time and/or type of chosen (periodic) boundary conditions (not considered here) might drastically change the influence of mechanical contact interaction on the filament diffusion behavior. Nevertheless, this result represents a valuable first indication that mechanical contact may decrease diffusivity noticeably. Beyond this example, there are many questions of interest in this field of application, e.g. the influence of mechanical contact interaction on the development of thermodynamically stable or unstable phases in cross-linked biopolymer networks (see e.g. [24]), where a robust contact simulation framework such as the one proposed in this contribution is of essential importance.
6.6. Example 6: Simulation of the static twisting process of a rope In this last example, the static twisting process of a rope will be investigated. The considered rope is built from 7 × 7 individual fibers with length l = 5, circular cross-section of radius R = 0.01 and Youngs modulus E = 10 9 . The arrangement of the initially straight fibers in seven sub-bundles with seven fibers per sub-bundle is illustrated in Figure 17(a). For spatial discretization, we use 10 beam elements per fiber. The contact parameters have been chosen as ε = 5.0 · 10 5 ,ḡ = 0.1R = 0.001 in combination with seven 5-point integration intervals per element. In the first stage of the twisting process, each of the seven sub-bundles is twisted by four full rotations within 80 static load steps. The twisting process is performed in a Dirichlet-controlled manner, such that the cross-section center points at one end of the sub-bundles (front side in Figure 17) are moving on a circular path (see also Example 2 of Section 6.2) with respect to the individual sub-bundle center points, while the corresponding points at the other end of the sub-bundles (back side in Figure 17) remain fixed. The deformed configurations at characteristic load steps after one, two, three and four full rotations are illustrated in Figures 17(b)-17(e). In the second stage of the twisting process, all seven subbundles together are twisted by one further rotation within 20 additional static load steps. This time, the cross-section center points are moving on a circular path with respect to the center point of the entire 7 × 7-rope. The deformed configuration at the end of this twisting process is illustrated in Figure 17(f). While the cross-section center points of all fiber endpoints at one end of the rope (front side in Figure 17) are fixed in axial direction, the cross-section center  Figure 17) are free to move in axial direction.
Additionally, a constant axial tensile forcef ax = 1000 acting on each of these axially freely movable fiber endpoints provides axial pre-stressing during the entire twisting process. In contrary to Section 6.2, the fiber endpoints are simply supported but not clamped. Consequently, Dirichlet conditions are only applied to the positional degrees of freedomd i at the endpoints but not to the tangential degrees of freedomt i . As already mentioned in Section 6.2, each individual fiber is free of mechanical torsion at the end of this twisting process, since only external (contact and reaction) forces, but no external torsional moments are acting on the fibers. Nevertheless, of course, an overall external axial torque resulting from the moment contributions of the reaction forces at the beam endpoints with respect to the centerline of the rope is necessary in order to guarantee for static equilibrium of the twisted rope at different load steps. The corresponding evolution of this external axial torque during the deformation process normalized by the maximal torque occurring at load step 100 is plotted in Figure 18(a). Interestingly, the evolution of the twisting torque over the twisting angle is almost linear within the two stages of deformation, i.e. the behavior of the rope is similar to the twisting response of a slender continuum. The higher slope in the second twisting stage, where all sub-bundles are twisted with respect to the centerline of the rope, results from the increased overall elastic resistance. The external work which is required in order to perform the considered twisting process in a quasi-static manner is proportional to the area enclosed by the graph of the twisting torque evolution and the horizontal axis of Figure 18(a). From a purely mechanical point of view, it is quite obvious that the presented example, which is dominated by line-to-line contact interaction along the entire length of the rope, should better be modeled by a line-to-line than a point-to-point type contact formulation. However, we also want to motivate this choice from a mathematical point of view. Thereto, the minimal contact angle min (α) occurring within the entire rope at a specific load step has been plotted over the first 80 load steps in Figure 18(b). As expected, the value of this minimal angle increases with increasing twisting angle/load step. In order to mathematically evaluate the possibility of applying the point-to-point contact formulation to this example, we will compare this minimal contact angle with the lower bound α min according to (34), above which a unique bilateral closest point solution can be guaranteed. To be able to do so, we have additionally plotted the evolution α min = arccos (1 − 2Rκ max ) over the first 80 loads steps. Here, R = 0.01 represents the cross-section radius andκ max is the maximal curvature value occurring in the entire rope for the considered load step. As can be seen from Figure 18(b), the curve representing the actual minimal angle min (α) lies entirely below the curve representing the minimal admissible angle α min . Thus, a unique bilateral closest point solution cannot be guaranteed and hence, as expected, the point-to-point contact formulations is not suitable for this example.

Conclusion
Within this contribution, a new finite element formulation describing the contact behavior of slender beams in complex 3D contact configurations involving arbitrary beam-to-beam orientations has been proposed. It has been shown by means of a mathematically concise investigation of standard point-to-point beam contact models that these formulations fail to describe a considerable range of practically relevant contact configurations as consequence of a non-unique bilateral closest point projection. In contrary, the proposed line-to-line formulation models contact interaction of slender continua by means of distributed line forces. It has been shown analytically that the corresponding unilateral closest point projection relevant for this line contact formulation always possesses a unique solution and thus is applicable for any geometrical contact configuration. By means of theoretical and numerical investigations, different contact discretizations, i.e. Gauss-point-to-segment or mortar type formulations, as well as different constraint enforcement strategies, based on penalty regularization or Lagrange multipliers, have been evaluated.
On the basis of these comparisons, it has been concluded that a penalty-based Gauss-point-to-segment formulation is most suitable for the considered range of beam-to-beam contact applications. On the one hand, the penalty regularization of the contact constraint, which can be interpreted as mechanical model of the cross-section stiffness, has been shown to be indispensable when employing beam models based on the assumption of rigid cross-sections. On the other hand, as compared to mortar-type approaches, the proposed Gauss-point-to-segment contact discretization is favorable in terms of computational efficiency and implementation effort. Additionally, it has been argued and verified numerically that the required range of penalty parameters does typically not induce contact-related locking phenomena when considering thin beams. The proposed formulation is supplemented by a consistently linearized integration interval segmentation that avoids numerical integration across strong discontinuities at beam endpoints.
In combination with a smoothed contact force law and the applied C 1 -continuous beam element formulation, this procedure drastically reduces the numerical integration error. It has been verified numerically that this reduction of the integration error is an essential prerequisite in order to obtain optimal spatial convergence rates. Moreover, the resulting line-to-line contact algorithm has been supplemented by the contact contributions of the beam endpoints. It has been shown, that these contributions are likely to occur in systems with arbitrarily distributed and oriented fibers.
Finally, it has been verified by means of several numerical examples that all the described individual model components are necessary in order to obtain an accurate, consistent and robust contact algorithm that can model complex systems of slender fibers with arbitrary contact configurations. Amongst others, a new numerical test case suitable for line-to-line contact scenarios has been designed and a corresponding analytical solution based on the Kirchhoff theory of thin rods has been derived. This test case in combination with the associated analytical solution can serve as valuable benchmark for the proposed model but also for future beam-to-beam contact formulations.

Appendix A. Residual contributions and linearization of applied beam element formulation
After some reformulations of (2), the residual contributions r int , r kin and r ext of one torsion-free beam element according to Section 2 due to internal, inertia and external forces can be derived: Here, we have introduced the following additional abbreviations: Consequently, we obtain the following expressions for the corresponding linearizations k int , k kin and k ext : whered ,d is typically a constant factor depending on the applied time integration scheme and S(.) is a skew-symmetric matrix that represents the cross-product, i.e. S(a)b = a × b ∀ a, b ∈ 3 . Additionally, we have: It can easily be shown that in the absence of external moments, i.e.m =m = 0, the overall stiffness matrix is symmetric. If the MCS method is applied, it is sensible to slightly reformulate the element residual contribution due to axial tension. Eventually, the original contribution r int,EA and the alternative MCS-contribution read: Accordingly, the corresponding contributions of the axial tension terms to the element stiffness matrix yield: In equations (A.5) and (A.6) applies the summation convention over the repeated indices i and j.

(B.4)
Depending on the case (point-, line-or endpoint-contact), (B.4) can be simplified due to n T r 1,ξ = 0 and/or n T r 2,η = 0. In case of endpoint contact, only the partial derivatives dξ c /dd 12 and dη c /dd 12 have to be adapted, while all other terms remain unchanged. In case of contact between an endpoint of beam 1, i.e. ξ c = −1 or ξ c = 1, with a segment η c ∈ [−1; 1] on beam 2, we consider the second line of (B.2) in order to determine dη c /dd 12 , while dξ c /dd 12 vanishes: Appendix C. Linearization of the line-to-line contact formulation We briefly want to repeat the linearization (50) of the contributions r i j con,1 and r i j con,2 of one individual Gauss point: We focus on the most general case with an integration interval segmentation being applied on both sides of the slave element. In the line contact case, the orthogonality condition p 2 on beam 2 is relevant. Its linearization reads: (ξ i j ,η c (ξ i j )) .
Inserting equations (D.2)-(D.7) into the equilibrium equations (D.1) leads to only one remaining relation − τ 1 + κm − τm b + κ f +f n = 0, (D.8) that has to be satisfied by the system parameters, while the other three equilibrium equations of (D.1) are satisfied automatically. From the family of solutions in (D.8), we restrict ourselves to one with vanishing mechanical torsion: (D.9) Altogether, equations (D.8) and (D.9) postulate the following requirement for the penalty parameter: (1 + )r (r 2 + h 2 )g 0 EA + EI(1 + )h 2 (r 2 + h 2 ) 2 . (D.10) In a next step, the Dirichlet boundary conditions have to be determined. The relation r = R − |g 0 |/2 for the helix radius appearing in (62) stems from the simple observation that the distance between the two helix centerlines has to satisfy 2r = 2R − |g 0 | in order to generate the required gap g 0 . With r being defined this way, the derivation of the conditions (64) and (65) is trivial in order to end up with a helix with radius r. However, the condition (63) for the axial displacement requires some further calculations. Thereto, we have to express the required constant axial tension = 0.01 of the helix as a function of the total length l c of the deformed helix in order to determine the helix slope h: Equation (D.11) yields the required helix slope in case of a given helix radius r and a prescribed axial tension . The required axial displacement u of the right endpoint follows from (62) and its initial position according to (61) as ∆d l 1,z = ∆d l 2,z = u = 2πh − l. (D.12) Finally, the corresponding Dirichlet-conditions for the tangential degrees of freedom have to be determined. The constant bending moment m b along the beam has to be considered by means of proper moment boundary conditions at the beam endpoints. According to [23] (see e.g. Table 1), the virtual work contribution of an external moment vector m j b b j FS in b-direction at the left/right boundary node j = l, r leads to a residual entry into the n-component of the corresponding nodal tangential degrees of freedom. Since the local n-directions coincide with the global x-directions at the beam endpoints here, it is sufficient to prescribe the x-components of the nodal tangents via Dirichlet constraints in order to enable proper reaction moments. According to the analytical solution in (62), these x-components have to vanish: ∆t l 1,x = ∆t l 2,x = ∆t r 1,x = ∆t r 2,x = 0. (D.14) Since we have chosen the system parameters in a way that leads to vanishing mechanical torsion, no additional torsional external moments have to be applied at the beam endpoints, i.e. m l t l FS = m r t r FS = 0. This is the reason why the application of the torsion-free Kirchhoff beam element presented in Section 2 is justified for this example and leads to the correct mechanical solution. As explained in [22] and [23], the axial components of the nodal tangents represent the axial tension ( j = t j − 1) at the nodes and cannot be prescribed, but are a part of the FE solution.