Frictional interactions for non‐localized beam‐to‐beam and beam‐inside‐beam contact

This contribution presents the extensions of beam‐to‐beam and beam‐inside‐beam contact schemes of the same authors towards frictional interactions. Since the schemes are based on the beams' true surfaces (instead of surfaces implicitly deduced from the beams' centroid lines), the presented enhancements are not only able to account for frictional sliding in the beams' axial directions, but also in the circumferential directions. Both the frictional beam‐to‐beam approach as well as the frictional beam‐inside‐beam approach are applicable to shear‐deformable and shear‐undeformable beams, as well as to beams with both circular and elliptical cross‐sections (although the cross‐sections must be rigid). A penalty formulation is used to treat unilateral and frictional contact constraints. FE implementation details are discussed, where automatic differentiation techniques are used to derive the implementations. Simulations involving large sliding displacements and large deformations are presented for both beam‐to‐beam and beam‐inside‐beam schemes. All simulation results are compared with those of the frictionless schemes.


INTRODUCTION
Beams are widely used to describe the mechanical behavior of slender structures. [1][2][3][4][5][6][7] The contact description between the beams (with or without friction 8,9 ) is often an essential part of these mechanical models.
Large-deformation frameworks to treat contact between beams with rigid cross-sections can be classified in different ways. Most of the schemes are formulated for shear-undeformable beams, [10][11][12][13] whereas others can also treat shear-deformable beams. 14 Most schemes can only be used for beams with circular cross-sections, [10][11][12] whereas a few can also be used for beams with elliptical cross-sections. 8,13,[14][15][16] Most schemes can only be used for small contact areas 8,15 (assuming point-wise contact interactions), whereas others are able to treat finite contact areas. 9,10,14,16 Of all the frameworks able to treat frictional contact between beams, only a few are capable to not only account for frictional sliding in the beams' axial directions, but also in the circumferential directions. 8,17 Finally, all frameworks are beam-to-beam contact schemes, that is, they repel beams if they touch each other. Only recently, we have proposed a beam-inside-beam contact approach that enforces one beam to remain inside another beam. 14 Its goal is the exact opposite of beam-to-beam contact schemes: to keep beams embedded inside each other.
To date, no frictional beam-to-beam contact scheme exists for the general case of shear-deformable beams, with both circular and elliptical cross-sections, that experience finite contact areas (non-localized contacts). A frictional beam-inside-beam contact scheme is also missing.
The aim of this contribution is to fill this gap by extending our previously developed contact frameworks 14,16 towards frictional interactions. Because the frameworks use the true beam surfaces to quantify penetration -instead of surfaces implicitly deduced from the beams' centroid lines -frictional sliding is not only quantified in the beams' axial direction, but also in the circumferential direction. To this end, the local problem to quantify penetration (or exclusion) does not only involve the unknown surface parameters of the master beam. It also involves one unknown surface parameter of the slave beam. This (circumferential) surface parameter describes the location of the (slave) contact point around the cross-section's perimeter for which penetration is quantified.
Although the aims of the beam-to-beam and beam-inside-beam contact frameworks are the exact opposite, their methodology is similar. Both schemes follow a master-slave approach, and select cross-sections along one of the two beams in contact (the slave). For each selected cross-section, the penetration between two beams is quantified (although arguably for the beam-inside-beam scheme it may be a called a measure of exclusion, instead of a measure of penetration). This measure is then used to establish a contact virtual work in order to repel beams from each other (for the beam-to-beam scheme) or to keep them embedded (for the beam-inside-beam scheme). The similarities between the beam-to-beam and the beam-inside-beam schemes also manifest when quantifying the amount of frictional sliding.
The structure of the article is the following. The penalty formulation for large deformation frictional contact in the continuum setting is provided in Section 2. Section 3 details the spatial discretization of the finite-element framework and provides implementation details. Section 4 presents a set of multibody numerical examples for beam-to-beam contact and beam-inside-beam contact and compares the results of the frictional frameworks with those of the frictionless frameworks. Conclusions are provided in Section 5.

Parametrization of the beams' surface
We consider beams consisting of rigid cross-sections which are attached to the beam's centroid lines at the cross-sectional centers of gravity. The beams of interest can deform in two ways. First, their centroid lines can elongate, bend, and revolve. Second, due to shear deformation, the normal vector to the cross-section's plane is not necessarily aligned with the tangent to the centroid line, see Figure 1.
In the following, we consider beam  and its surface denoted by . We also consider a two-parameter vector function, x = x(h), that maps any surface point on , with surface coordinates h = [h 1 , h 2 ] T , to its location in the global F I G U R E 1 A typical cross-section (in gray) in the undeformed (left) and deformed (right) configurations for (A) a plain cross-section, (B) a hollow cross-section. The centroid-line is presented with a dotted line in both configurations. A surface point is represented with a red dot. Surface basis vectors are presented with orange arrows. The normal vector to the cross-sections' plane is presented with a blue arrow (A) (B) coordinate system. h 1 ∈ [0, L] denotes the arc-length parameter of the beam's undeformed centroid line x 0c ∶ (0, L) → R 3 and L denotes the length of the undeformed centroid line. h 2 ∈ [0, 2 ] denotes a circumferential parameter of the perimeter of the cross-section attached to x 0c (h 1 ) (see References 8,[14][15][16].
The location of a surface point in the undeformed configuration can be obtained from: where v 0 denotes a vector in the plane of the cross-section attached to x 0c (h 1 ).
In the deformed configuration, the centroid line deforms according to: where u ∶ (0, L) → R 3 denotes the centroid-line's displacement. The location of any surface point in the deformed configuration with surface coordinates h is given by: where denotes the deformation mapping relating the location of a surface point in the deformed configuration to its location in its undeformed configuration. As the beams' cross-sections are assumed to be rigid, there exists a rotation tensor where denotes the field of variables used to parametrize SO(3), for example, quaternions 1,2 or the smallest rotation with respect to a reference triad 4 or rotation vectors 3 as used in this work. For further use, we define the local basis vectors at surface point x. For ∈ {1, 2}, the covariant tangent surface vectors are given by: Assuming proper orientation of 1 and 2 , the outward-pointing normal unit vector to the surface is computed from: where = −1 for the inner surface of the master body for the beam-inside-beam contact scheme (such that n points towards the center of the hollow beam), and 1 otherwise. Note that local basis { 1 , 2 , n} is not necessarily orthonormal in the deformed configuration, meaning that in general 1 ⋅ 2 ≠ 0 and || || ≠ 1. Remark: If a geometrically exact beam formulation is used, in Equation (4) is an element of the solution. Different FE formulations exist to suitably treat the parametrization of (see References 7,18), but are not discussed here in order to not distract from the contact formulation. In our approach, we use the vector-like parametrization of three-dimensional finite rotations provided in Reference 3.

Contact kinematics
In the present contribution, we focus on non-localized contact between beams, meaning that the contact area is finite. This stands in contrast to localized contact, which occurs on a narrow area of the surface. In general, localized contact is modeled with a point-wise contact interaction, see References 8,11,12,15, which is not suitable for non-localized contact cases.
Although the present framework is applicable for multibeam systems, for the sake of clarity we only discuss systems of two beams, denoted by  I and  J . To repel these two beams in case of penetration (in the beam-to-beam scheme) or to keep them embedded (in the beam-inside-beam scheme), we integrate the contact virtual work along the centroid line of  I .  I is therefore called the slave and  J the master. 19

Local problem and normal gap
In the following, we focus on a cross-section of  I denoted by , that is attached to centroid point x I c (h  ), where h  denotes the (fixed) arc-length parameter locating  along the centroid line. The perimeter of  is denoted . We first quantify penetration of  into  J , where  J denotes the surface of  J . For that purpose, we determine two surface points, one on  and one on  J . These points are denoted by , respectively, and will be used to quantify penetration and relative sliding of contacting surfaces. x I and x J map surface parameters to the location of the corresponding surface point in the global coordinate system.
To determine the surface coordinates of these two points, denoted by h I and h J , we solve a set of equations. This is usually referred to as the local problem or the projection problem. 19 The local problem presented here involves three surface parameters: the circumferential parameter of , h I2 , and the two surface parameters of  J , h J1 and h J2 . h I1 does not need to be determined as  is fixed at Three of the equations we solve for are expressed as: where: denotes the array of unknowns. q denotes the array solution of Equation (7). In the following, a bar over a quantity indicates that this quantity is evaluated at the solution of the local problem. Variable g denotes an unknown scalar, for which we can write at the solution of Equation (7): where g denotes the so-called "gap vector." Hence, g = g N denotes the amount of penetration.
As the system of equations of Equation (7) is underdetermined, an additional equation is needed. Here, we impose that at the solution of the local problem, n I and n Jp must be orthogonal to J 2 (see Figure 2). n Jp denotes the (normalized) projection of n J on the plane spanned by vectors I 2 and n I . This plane has the following normal unit vector:ñ n Jp is then obtained from: The additional equation to be added to the system reads: where has the dimension of length and is used to ensure that components of f 1 and f 2 have the same units (if  is elliptical, we can for example set = a  where a  denotes the dimension of  along its largest semiaxis). The final set of (B) (A) F I G U R E 2 Solution of the local problem for (A) beam-inside-beam contact; (B) beam-to-beam contact. The slave cross-section attached to x I c (h  ) is shown in red. The master surface is shown in gray equations of the local problem is abbreviated as follows: Equation (13) is nonlinear and can be solved using Newton's method for which the following Jacobian is required:

Sliding increment
The true novelty of this contribution lies in the treatment of frictional contact interactions. We introduce the formulation directly in the time-discretized setting, that is, we assume that the simulation time is divided in numerous time increments. From here onward, subscript n refers to the previous time increment. If no subscript is present, the quantity refers to the current time increment. The particularity of the contact kinematics employed in this contribution prevents the use of conventional frictional frameworks, as for instance developed for the node-to-surface (NTS) approach. 20,21 In the NTS approach, it is sufficient to measure the relative sliding of the slave node (that has fixed surface coordinates) over the master surface. In the present contribution, however, also the circumferential coordinate of the slave contact point (h I2 ) generally varies between two time increments.
In more detail, as h I2 is not fixed, x I can be located at a different surface point on  in the previous and current configuration. Thus, both slidings of x I and x J must be incorporated. The change in the location of x I between two increments can be approximated as (see Figure 3):  Similarly, we introduce Δx J as: The fact that the position of the contact points are mapped to the previous configuration ensures that the sliding are located on different finite elements.
We express Δg T , the increment of tangential sliding between two increments, as the following frame-indifferent measure (see also Figure 3): where the summation on repeated indices holds. Components Δg T are approximated by the following projection: where I n denotes the contravariant basis vector in the configuration from the previous time increment: M I n are contravariant components of the metric tensor of the surface of  I in the previous configuration,  I n : where covariant components of the metric tensor of  I n are given by: and in Equation (21), I n denotes a covariant surface tangent vector in the previous configuration at surface coordinates h I n defined as: The measure of relative sliding proposed in Equation (17)  deformation-dependent via the coupling between q and p IJ , see Section 2.3.4. Thus, the first solution is adopted in this work.
Remark: If one assumes exact normal contact in the previous time increment such that: then the right-hand side of Equation (18) reduces to: In this case, no internal variables would need to be stored from the previous increment, while are needed in Equation (18). In the present contribution, however, we use Equation (18) because the penalty method is used (see below). Thus, condition (23) is not exactly fulfilled.

Normal contact
For a given cross-section , the impenetrability of  and  J is enforced via unilateral contact conditions: where: where T denotes the nominal traction vector, that is, the traction acting in the current configuration, yet integrated over the contact area in the reference configuration.

Penalty regularization of the constraint
The penalty method is employed to enforce the constraints of Equation (25). Typical for the method is that after enforcing the constraint, some residual penetration remains. This can be interpreted as some compliance of the otherwise rigid cross-sections (see References 14,16). Contact traction T N , acting between surface points x J and x I , is given by: where N > 0 denotes the penalty stiffness and ⟨⟩ denote the Macaulay brackets.

Coulomb's friction law
The Coulomb friction law 19,21,22 is used in this work to describe the magnitude of the tangential forces. For two bodies in contact, the corresponding constraints are specified by the limit friction condition, the slip rule and the complementarity condition, respectively: 20 where T T denotes the nominal frictional contact traction vector, the static friction coefficient, and Φ the yield function. Simply stated, Equation (28) implies that if ||T T || < − T N , surfaces do not slide with respect to each other. If they slide, the tangential traction, ||T T ||, is equal to − T N .

Penalty regularization of Coulomb friction
The penalty regularization is applied to regularize tangential contact constraints. In practice, it means that we allow for a small elastic relative displacement of the surfaces. The tangential gap vector reads: where g Tel n updated denotes the previous elastic tangent gap transferred to the current configuration (see Equations (30)-(36)). The tangential gap vector can also be split as follows: where g Tel denotes the elastic part of g T and Δg Tsl the increment in irreversible tangential sliding. Both quantities can be obtained from a return mapping procedure, see Equation (41) and thereof.

Transfer of g Tel n and update of the tangential gap
The (stored) elastic gap of the previous configuration, g Tel n must be transferred to the current configuration to properly treat rigid body motions. 20 In the current configuration, the contravariant components of the (projected) elastic tangential gap of the previous increment read: These components are associated to covariant basis vectors (see Equation (5)) in the current configuration such that: To preserve the elastic gap's norm, the following scaling of the components is performed: with: otherwise. (34) Contravariant components of g T are computed according to: such that: Once g T is computed, a return mapping procedure is employed to compute the split into g Tel and Δg Tsl , see Equation (30).

Return mapping procedure
First, the trial state is computed, assuming g T is entirely elastic: where T denotes the (user-defined) tangential penalty stiffness. Inserting T T tr into Φ in Equation (28).1 gives: Depending on the sign of Φ (T T tr ), frictional sliding occurs or not, which is referred to the sliding and sticking case, respectively. The two cases are treated as shown below.
In this case, ||T T tr || does not exceed the threshold value − T N in Equation (28). The trial state fulfills the Coulomb friction conditions Equation (28) and: We can note here that the constraint in Equation (28).3 is not respected in general as the penalty regularization allows small elastic sliding displacements, but is respected for Δg Tsl as Φ||Δg Tsl || = 0. The higher T , the lower the magnitude of g Tel for a given T N . However, large values of T increase the condition number of the global stiffness matrix. 23 In this case ||T T tr || exceeds T N and Equation (28).1 is violated. The trial frictional traction is corrected to take the limit frictional traction given by the Coulomb's law: In both cases (Φ (T T tr ) ≤ 0 and Φ (T T tr ) > 0), the elastic tangential vector g Tel is given by: where the tangential traction, T T , is defined by Equations (39) or (40) for the sticking and sliding case, respectively. Note that frictional sliding increment, Δg Tsl , can be explicitly retrieved from Equation (30).

Contact virtual work
In case of contact, a contact virtual work, Π c , is added to the virtual work equation of the system and the space of admissible variations is modified. 19 In the quasi-static settings as considered in this contribution, the virtual work including contact reads: where Π  i denotes the internal and external virtual work of beam  i (excluding contact interactions). Kinematic variables associated with  i are stored in p i and the associated test functions in p i . p i gathers the field of displacement of the centroid line, u  i , as well as the field of variables used to parametrize SO 3 ,  i . p i is only admissible Following Reference 21, by using the action-reaction principle, the contact virtual work, Π c , can be suitably expressed with respect to the slave part only. The infinitesimal virtual work produced by dL  I , an infinitesimal part of  I 's centroid line, denoted by d Π c , can be written as in Reference 14: where: dL  I is related to the differential of h I1 , dh I1 , by: Integration over all penetrated slave sections gives: h 1I L and h 1I U denote the lower and upper bounds of the integral, respectively, corresponding to the first and last penetrated slave sections if we assume a unique contact area.

Variation of the local parameters, q
For the purpose of linearization procedure introduced below, variations of local parameters q with respect to variations of the kinematic variables p IJ are needed. Again, the corresponding equations were previously provided in References 14, 16, but are repeated here to make the contribution self-contained. To express q, in terms of p IJ , we start from the stationarity of local residual f in (13) with respect to p IJ as follows: where we recognize the Jacobian of the local problem, H, defined in Equation (14). After rearrangement, we obtain: where:

Interpolation of the beams' surfaces
Beams  I and  J are now discretized with a series of consecutive beam finite elements (BFEs). The nodal variables of all BFEs are gathered in arrayp, and the associated variations in arrayp. The finite-dimensional trial and test functions, p h and p h , obtained by combiningp andp to properly chosen interpolation functions, replace their infinite dimensional counterparts p IJ and p IJ , respectively. A variety of BFEs exists. 6 They usually differ by whether or not shear deformations can occur, by the type of rotational variables used, the interpolation schemes employed for the different types of variables and the treatment of locking. In order to remain as general as possible, we do not restrict ourselves to one specific type of BFE in this section. We assume that the discretized surface is sufficiently smooth. In the numerical examples in Section 4, however, we have used two-nodes beam elements with rotation vectors as rotational variables. The associated discretized surface is discontinuous with gaps and overlappings at the nodes. To overcome this problem, we work with an alternative surface which possesses the desired continuity, which was previously introduced in Reference 14. Thus, each contact element constructed for every integration point involves two BFEs on the slave side and two on the master side. The contact element's nodal variables are denoted asp e (see below). In the following, we detail the procedure to obtain the contact contributions to the global force vector and the global stiffness matrix, denoted by r g and K g , respectively.

Contact residual and stiffness
The discretized form of the virtual work in Equation (42) leads to a set of nonlinear equations. Newton's method is generally used to iteratively determine global solutionp sol of the virtual work statement. This requires the linearization of Equation (42) around an estimate ofp sol , denotedp, which yields: where △p denotes an increment of the nodal variables. The global residual force column, r g , reads: where f int denotes the internal force column stemming from the contributions of all BFEs, and f ext the external force column. r c contains all the (assembled) contact contributions from all contact elements, where a contact element refers here to all cross-sections attached to an integration point (see below) along  I 's centroid line and their projection on discretized surface  J .
Since f ext generally does not depend onp, the global stiffness obtained after the linearization of r g , can be decomposed as follows: where K int denotes the stiffness matrix associated with the BFEs, and K c denotes the stiffness matrix associated with all contact elements.
The contact virtual work consists of the contributions of all contact elements: where S denotes the set of active contact elements (i.e., those for which g N < 0), and Π c,e denotes the contact virtual work associated with contact element e.p e denotes the column of nodal variables involved in this contact element.
If no smoothing procedure of the surface is required, each contact element involves two BFEs: the first one is part of the discretization of  I , and the second one is part of the discretization of  J . As mentioned above, however, if a smoothing of the beam's surface is necessary to improve its surface continuity, each contact element depends on several BFEs of  I and of  J . 14,16,20 The elements of  I and  J required to construct contact element e are denoted by  and  , respectively. The associated nodal variables are denoted byp  andp  such thatp .
Next, we discuss how to construct the contributions of a single contact element to the total force column and stiffness matrix. The force column and stiffness matrix associated with element e are denoted by r ce and K ce . We distinguish the contributions stemming from normal and tangential contact. Hence, r ce is decomposed as: where r cNe and r cTe denote the normal and tangential contact element contact contribution. Similarly, the contact element stiffness is decomposed as follows:

Force vector and stiffness of a single contact element
To numerically integrate Π c,N in Equation (46), n  IP integration points (to which a cross-section is attached where penetration is to be quantified) are placed along 's centroid-line. Whether we integrate along the centroid line of a beam element or an artificial smoothed centroid line constructed from several BFEs, the variable over which we integrate is denoted by (and hence, a mapping is constructed between h M 1 and if necessary). The contact virtual work of the normal contact interactions reads: Jacobian  = h  1 maps differential length dh  1 to differential increment d . The weight of the k th integration point is denoted by w k and its coordinate in the parameter space is denoted by k . For a given cross-section of  attached to x c ( k ), r cek in Equation (59) is expressed as: wherêp denotes the total derivative with respect to variables p performed by the automatic differentiation (AD) algorithm. 20,24 Here, the exception in AD indicates that the variation of the gap vector, g ( k ), is not influenced by the local variables as discussed in our previous work. 14 K cNek , stemming from the contribution of the k th integration point, can be obtained using AD. It allows to include the implicit dependency of local variables q on global variables p IJ (see Equation (48)): Similarly to Equation (59), the numerical integration of the frictional term Π cT in Equation (46) governing friction reads: where: The corresponding contribution to the tangent stiffness can once again be obtained with AD:

Contribution of all contact elements
r c and K c , which contain the contributions of all contact elements in set S, are assembled as follows: where A denotes the finite-element assembly operator.

NUMERICAL EXAMPLES
In this section, the proposed frictional framework is applied to three beam-to-beam examples and one beam-inside-beam example. So far, we have presented our contact framework without specifically referring to a certain type of BFEs. In all presented simulations, however, Simo-Reissner geometrically exact beam elements are used. The associated BFEs have a linearly interpolated centroid line. [1][2][3] This entails that the surface of a series of consecutive BFEs is C 0 -continuous only if the undeformed configuration is straight. If the string of BFEs is not straight in the reference configuration, the surface associated with such strings of BFEs is not C 0 -continuous and hence, contact constraints are hard (if not impossible) to apply.
An approximated but C 1 -continuous surface description was proposed in Reference 16 to alleviate this issue, see also Figure 4. This surface is used in all the numerical examples below and is only used for the contact treatment; it has no influence on the employed beam theory, nor on the BFEs.

Beam-to-beam contact: A 1+6 strand in tension
Strands are assemblies of wires used in tire reinforcement or as components of wire ropes. In this example, we focus on a 1+6 (a central wire surrounded by six helical wires) strand subjected to tension. The geometry of the strand is reported in Table 1. A gap of 0.05 mm is inserted between the central wire and its six surrounding wires in order to prevent numerical issues at the strand's ends. The Young modulus is set to E = 188 GPa and its Poisson ratio to 0.3. Only one pitch of the strand is modeled. Twenty BFEs are used to discretize each wire. One strand's end is clamped while the nodes at the other end are moved in the z-direction by a final displacement u end such that the strand's axial engineering strain reaches strand = 0.015, see Figure 5. Only one Gauss integration point is used per smoothed patch for the integration of the contact virtual work.

F I G U R E 4
Two connected beams j and j + 1 whose centroid lines is shown with a green and orange dashed line, respectively. Their nodes are shown with gray circles. The smoothed centroid line constructed for these two beams is shown with a plain blue curve. The associated control points are shown with red circles. As an example, the quadrature points of the three-point Gauss-Lobatto rule are indicated by black lines and the single quadrature point of the one-point Gauss-Legendre rule is indicated by a green cross TA B L E 1 Geometrical parameter for the 1+6 strand (a straight wire surrounded by six helical wires) 25 Central wire diameter 3.94 mm Helical wire diameter 3.73 mm Pitch length 115 mm F I G U R E 5 Beam-to-beam contact 4.1: 1+6 strand subjected to an axial displacement u end at one of this end while the other end is clamped F I G U R E 6 Beam-to-beam contact 4.1:

Component of the reaction forces in the z-direction
The initial penalty stiffness is estimated using the Hertz theory for parallel cylinders, 26 resulting in N = E 2 8(1− ) 2 = 82 × 10 9 N/m. As the tensile force gradually increases with the loading, the contact forces between the wires increase. Thus, the residual penetration due to the use of the penalty method increases. As soon as the penetration measured at an integration point according to Equation (9) falls below −5% of the smaller slave beam cross-section dimension at the end of a time step, the penalty stiffness of all quadrature points of the slave BFE/smoothed curve is increased by 10% and the corresponding time step is repeated.
The tangential penalty stiffness is set to T = 8.2 × 10 9 N/m and allows for reasonably small elastic relative displacements of the contacting surfaces without causing convergence issues. The friction coefficient is set to = 0.115. The loading is applied in 150 equally spaced time steps. Reaction forces in the axial direction are shown in Figure 6. It can be observed that for small strains, the response predicted by our model has the same slope as the one predicted by Costello's theory. Figure 6 also shows the reaction force predicted for two equivalent strand models in commercial FE software Abaqus © . One uses C3D8R elements (hexahedra with reduced integration and hourglass control) and the other B31 elements (linearly interpolated beam elements). Both simulations give a similar response and capture the first experimental points well. As the material behavior is elastic, the reduction of the reaction force's slope due to the wires' plastification is not captured. Figure 7 shows the evolution of the number of active contact points. Initially, the wires are not in contact due to the (small) initial gap between the core wire and the helical ones. Then, as the strand elongates, wires come in contact in the center of the strand. The contact propagates towards the strand's ends. Despite a relatively large number of contact interactions, only a few global Newton-Raphson iterations are necessary to converge, see Figure 7

The choice of the quadrature rule and penalty stiffness
The influence of the choice of the quadrature rule and the penalty parameter is studied. The goal is to investigate if contact locking appears when more quadrature points are used and/or for high values of the penalty stiffness. The numerical example of Figure 5 is repeated with: To study the influence of the quadrature rule on contact traction T N , the contact traction of the different contact elements between the central wire and a chosen peripheral wires have been reported for: • k = 1 but a varying number of Gauss-Legendre quadrature points, see Figure 8, • 1 and 2 Gauss-Legendre quadrature points but k ∈ {0.5, 1, 2, 10}, see Figure 9, • k = 1 but a varying number of Gauss-Lobatto quadrature points, see Figure 10, • 3 Gauss-Lobatto quadrature points but k ∈ {0.5, 1, 2}, see Figure 11. Figure 8 shows that when the number of Gauss-Legendre quadrature points nQP is more than two, T N oscillates along the contact line. Figure 9(A) shows that for a 1-point Gauss-Legendre quadrature rule, T N does not oscillate, even for k = 10. Figure 9(B) shows that for a 2-point Gauss-Legendre quadrature rule, T N does not oscillate for k ≤ 1, but higher values of k (and thus N ) induce oscillations. For a Gauss-Lobatto integration rule, Figure 10 shows that even with a number of quadrature points as low as 2, oscillations of T N are present. Figure 11 shows that increasing k increases the amplitude of these oscillations. None of the simulations using k = 10 converges with a Gauss-Lobatto quadrature. For all converging simulations, the evolution of the reaction force and the  26 The tangential penalty stiffness is set to T = N ∕10 = 4.3 × 10 9 N/m.
As the rotation increases, beams wrap around each other, which causes the contact area and the number of penetrated sections to increase (bottom diagram in Figure 13). The deformed configuration is similar for all simulations and is shown on the right in Figure 12. Despite the substantial number of penetrated sections, the number of global iterations required to converge according to ||f int + r c − f ext || < 10 −8 , remains low (top diagram on Figure 13). This is thanks to the proper linearization of r c , see Equations (61) and (64), with the AD tool. Figure 14 reports the evolution of the reaction forces and the total torques around z axis at the support. The influence of friction is relatively small for this example, although friction does have a substantial influence on the reaction force in the axial direction of the beams (top-right diagram in Figure 14).

Beam-to-beam contact: Twisting and pulling of a fiber
In this example, four beams of length L = 70 × 10 −3 m are aligned along the z direction ( Figure 15). The beams' cross-sections are circular with a radius of 3.6 × 10 −3 m. During the first part of the simulation, the displacements and rotations of all beam nodes at one end are fully restrained. The sections at the other end are rotated around the z axis with an angle of 180 • . This loading is applied in 1800 increments. Each beam is discretized with 40 BFEs. The nodes at the rotated end of the beams are only free to move along z while the other kinematic variables are prescribed. During the second part of the loading, one of the beam is extracted from the deformed structure by pulling it (at the end node) in the z direction. The initial penalty stiffness is once again estimated using Hertz theory for the case of perfectly parallel cylinders in contact, resulting in N = 4.4 × 10 10 N/m. The tangential penalty stiffness is set to T = N ∕10 = 4.4 × 10 9 N/m. A one-point Gauss-Legendre quadrature rule is used to integrate the contact virtual work. The simulation is performed without friction and with a friction coefficient of 0.25. It is also performed for circular (with a radius of 3.6 × 10 −3 m) and elliptical cross-sections (with the same cross-sectional area as the circular cross-sections, with a = 2.16 × 10 −3 m).
The undeformed configuration as well as the deformed configuration at the end of the two part of the simulation is shown in Figure 15 for circular and elliptical cross-sections. The number of global iterations required to attain the desired accuracy, ||f int + r c − f ext || < 10 −8 is shown in Figure 16. Figure 17 shows the reaction force, revealing that both friction as well as cross-sectional shape have a substantial influence. The force-displacement curves of the cases with friction are less smooth than their frictionless counterparts. This is because of the change in the sticking/slipping status of contact points.

Beam-inside-beam contact: Insertion
This example involves an initially straight thin beam that is pushed in a hollow beam (see Figure 18(A)). Both beams have elliptical cross-sections. Initially, only a small part of the inner beam is present in the hollow one. The curved part of the hollow beam's centroid-line is a half-circle of radius 150 × 10 −3 m. The kinematic variables of the outer beam's end node near the thin beam's insertion location are restrained. Its wall thickness is 10 −3 m and the lengths of its cross-sectional semi-axes are a = 20 × 10 −3 m and b = 16 × 10 −3 m. The inner beam has a length of 54 × 10 −2 m, and a Young's modulus of 100 GPa. The cross-sectional shape is given by a = 5.4 × 10 −3 m and b = 4.3 × 10 −3 m. The outer hollow beam is more compliant with E = 10 GPa. The z-displacement of the inner beam's end node furthest away from the outer beam is prescribed to reach 54 × 10 −2 m in 300 increments, whilst the other kinematic variables at this end node are restrained. The thin beam is discretized with 20 BFEs and the hollow beam with 45 BFEs. A one-point Gauss-Legendre quadrature rule is used to integrate the contact virtual work.
Once again the simulation is performed with different static friction coefficients: 0, 0.5, and 1. Both beams have a Poisson's ratio of 0.33. The initial penalty stiffness is set to N = 1 × 10 3 N/m which is several order of magnitude less than for the examples in Sections 4.2 and 4.3. A higher penalty stiffness causes convergence issues with oscillations of the contact status, meaning that some contact elements penetrate the thick beam wall and then detach (g N > 0) from one iteration to the next. The tangential penalty stiffness is set to T = 1 × 10 2 N/m and allows acceptably small elastic tangential gaps, while allowing the (global) Newton-Raphson scheme to converge.
Both structures deform due to contact, see Figure 18. Figure 19 shows that numerous sections of the inner beam penetrate the wall of the outer beam, which indicates that the contact is non-localized. Only a few iterations are necessary to reach convergence for the "beam-inside-beam" framework as Figure 19 shows. Figure 20 shows the component of the reaction force in the z direction. The presence of friction clearly has a substantial influence on the force-displacement response, indicating that friction not only influences the results of beam-to-beam contact schemes, but also those of beam-inside-beam contact schemes.
The contact of the tip of the inner beam plays a crucial role during the entire simulation. It is enforced with a contact at the closest pair of surface points between the surface of the last section of the inner beam and the inner surface of the tube. A similar contact element was used in Reference 8. This contact interaction needs a specific treatment in order to avoid a complete loss of contact between the tip and the inner surface when the tip slides out. More details can be found in Reference 14.

CONCLUSION
This contribution presents the extension of beam-to-beam and beam-inside-beam contact frameworks towards friction. It is applicable to shear-deformable and shear-undeformable beams with circular and elliptical cross-sections. The framework is not only able to account for frictional sliding in the beams' axial direction, but also in the circumferential direction. It is suitable for non-localized contact interactions, occurring for instance when beams are parallel to each other or wrapped around each other. The contact kinematics are formulated in terms of surface parameters of the master and slave beams. Thus, the introduced framework can be exploited for a variety of BFE formulations, provided that their cross-sections are rigid and their discretized surface is C 1 -continuous. An important specificity of the introduced framework is that both slave and master contact points change their location at the beams' surfaces during the relative tangential motion of beams. This is unlike common node-to-segment approaches, in which the location of the slave contact points are fixed. We propose a measure of relative tangential displacement that is frame indifferent and does not involve higher-order dependencies on global kinematic variables. Thanks to that, the measure is suitable for finite-deformation and finite-sliding problems, and also leads to computationally efficient linearization.
The presented formulation is shown to efficiently regularize contact constraints in a series of numerical examples for beam-to-beam and beam-inside-beam contact interactions, with and without friction for beams with circular and elliptical cross-sections. Even if numerous contact interactions are present and the beams' deformations, rotations and curvatures are substantial, only a few global iterations are necessary to converge. This is thanks to the consistent linearizations achieved using AD, which automatically incorporates the dependencies between global and local variables.
While the present formulation is suitable for non-localized beam-to-beam contact, it is not truly adapted to enforce localized contact. In such cases, frameworks as presented by Gay Neto et al., 8,15 which enforce contact at a single pair of points, seem more accurate. A scheme that automatically decides whether to enforce point-wise contact or non-localized contact remains for future work. Although Reference 10 provides a geometrically based choice of the contact formulation which depends on the spatial arrangement of the two beams' centroid-lines, it cannot be used here because it is limited to beams with circular cross-sections. Another geometrical criterion would be necessary to decide which formulation to employ for beams with elliptical cross-sections (point-wise surface-to-surface contact 8,15,28 or surface-to-surface non-localized contact).