Non-localised contact between beams with circular and elliptical cross-sections

The key novelty of this contribution is a dedicated technique to efficiently determine the distance (gap) function between parallel or almost parallel beams with circular and elliptical cross-sections. The technique consists of parametrizing the surfaces of the two beams in contact, fixing a point on the centroid line of one of the beams and searching for a constrained minimum distance between the surfaces (two variants are investigated). The resulting unilateral (frictionless) contact condition is then enforced with the Penalty method, which introduces compliance to the, otherwise rigid, beams’ cross-sections. Two contact integration schemes are considered: the conventional slave-master approach (which is biased as the contact virtual work is only integrated over the slave surface) and the so-called two-half-pass approach (which is unbiased as the contact virtual work is integrated over the two contacting surfaces). Details of the finite element formulation, which is suitably implemented using Automatic Differentiation techniques, are presented. A set of numerical experiments shows the overall performance of the framework and allows a quantitative comparison of the investigated variants.


Introduction
Many engineering materials such as paper, fabrics or opencell high-porosity foams consist of slender fiber-like constituents at the microscopic scale [9,25]. Various approaches to describe their mechanical behavior explicitly incorporate their discrete micro-structure [2,12,17,22,33]. In many cases beam models and beam finite elements (BFEs) are used to represent single fibers, yarns or struts [1,3,16]. It is often crucial to incorporate beam-to-beam contact in order to obtain accurate mechanical predictions. However, due to the specificity of beam kinematics, standard techniques developed to treat contact between 3D solids cannot be directly adopted.
Beam-to-beam contact schemes are in general built upon assumptions on the contacting systems, which restrict their use to specific contact scenarios. The formulation of a particular contact scheme is typically determined by three main issues: (i) whether or not contact remains localized (in other words if contact interactions are confined to a small part of the beams' surfaces), (ii) what are the beams' cross-sectional shapes and (iii) are the cross-sections rigid or deformable. This contribution is limited to rigid cross-sections.
A point-wise contact force approach is generally used to model localized contact between beams. For beams with circular sections for example, the contact conditions are enforced at the closest pair of centroid points [24,34,36]. If the beam's cross-sections are elliptical, the contact force can be applied at the closest pair of surface points where the tangent planes of the contacting surfaces are parallel. This was demonstrated by Neto et al. [7,8] for the frictionless and frictional case, respectively.
If no unique minimum distance between the beams' surfaces exists (e.g. in case of parallel beams or if one beam is wrapped around the other), the aforementioned frameworks are difficult to apply. In such cases, and if both beams in contact have circular sections, Chamekh et al. [4] (for beams undergoing self-contact) and Meier et al. [23] have demonstrated that contact can be modeled as a continuous force acting along one of the beams' centroid lines. If elliptical sections are employed however, the contact scheme of [23] is not directly applicable as the centroid lines are not sufficient to locate possible contact interactions.
An alternative is then to integrate contact forces (contact virtual work) over the surface [21]. In such an approach, fixed material points on one of the beams' surface are projected on the other beam's surface to then determine if they are penetrated. As each projection requires the solution of a (small) non-linear problem, and because many projections are required to accurately approximate the contact area, the associated computational costs are substantial.
In the approach proposed in the current contribution, the contact virtual work is integrated along the centroid line of one of the beams (slave). The contact kinematics are based on an appropriately defined projection procedure between the contacting surfaces, which assumes that the cross-sections are rigid. In this procedure, points are fixed along the centroid line of the slave body. For each point, the projection problem solves for the circumferential parameter of the slave surface and both surface parameters of the counter-surface (master).
Two different projection schemes are presented and compared. The first one is based on the minimisation of a constrained scalar function and the second one is posed as a set of equations to solve. The latter approach is expected to be more efficient and stable because surface derivatives of lower degrees are involved, which is appealing when interpolated surfaces have a reduced continuity. As the integration is performed along a curve and not over a surface (as in [21]), the computational costs are substantially reduced. Consequently, larger models can be investigated with similar computational efforts, thereby easing the framework's application to industrially and scientifically relevant problems.
In the present contact framework we adapt and study two approaches to integrate the contact virtual work. As a first choice, we use a conventional single-pass algorithm, a master-slave approach, in which we arbitrarily choose one of the beams' surface (slave beam) to integrate contact virtual work. Consequently, the associated contact framework is biased [28,29].
To avoid this issue, we adapt the idea presented in [28,29], the so-called "two-half-pass" approach, which symmetrically treats both contacting surfaces, and which has so far only been considered for conventional finite elements and not for beam-to-beam contact. Consequently, this framework is unbiased. In this approach, similarly to the well known two-pass approach, contact conditions are evaluated twice by changing the roles of the contacting surfaces (the slave body during the first half-pass becomes master during the second half-pass and inversely for the master body). The difference between two-half-pass and two-pass algorithms is that for the two-half-pass, contact tractions are only acting on the slave surface for each half-pass. This paper is organized as follows. In Sect. 2, the kinematics of the proposed contact framework are presented. This includes a discussion on how to determine whether a section penetrates another beam and if so, how to quantify the amount of penetration. Also, a length-specific contact virtual work is introduced, which is integrated over the slave beam's centroid line. We specify this for single-pass and two-half-pass algorithms. In Sect. 3, the contact kinematics are adapted to the finite element method after the beams are discretized with BFEs. Implementation details are included. The numerical examples of Sect. 4 indicate the capabilities and efficiency of the contact framework. Finally in Sect. 5, conclusions are presented and possible extensions are discussed.

Continuum contact framework
This section explains the proposed framework to treat contact between beams with elliptical and circular rigid sections. Geometrically Exact Beam theory [10,11,23,27,30,31] is utilized here, but the contact framework can be adapted to other beam formulations with rigid cross sections. We first explain how surfaces of such beams can be parametrised. Based on these beams' surface parametrizations, we then introduce the contact kinematics for the novel contact scheme in the spacecontinuous setting. This framework will be adapted to the finite element method in Sect. 3.

Surface description of a geometrically exact beam
In this work, a continuous beam (not a beam finite element) is considered as a (slender) body, of which the cross-sections (i) do not deform, (ii) remain planar, (iii) their center of gravity form its centroid line and (iv) their normal vector can rotate with respect to the tangent of its centroid line (shear deformable). We consider beam B and the parametrization of its surface where h 1 and h 2 are longitudinal and circumferential parameters, respectively. L denotes the length of the centroid line of B in the undeformed configuration. The current centroid line position is obtained by adding displacement u ∈ R 3 to the initial centroid line location, X c ∈ R 3 : x c = X c + u. (2.1) As cross-sections are rigid and remain planar, x(h) can be obtained by adding the location vector of a centroid point, , to a vector, v(h) ∈ R 3 , that lies in the plane of cross-section C that is attached to x c (see Fig. 1): For elliptical sections, v can be expressed as follows: where a and b denote the two semi-axes of the elliptical section in the ellipse's principal directions. Note that in case a = b, a circular cross-section is recovered. Triad {e 1 , e 2 , e 3 } attached to C forms an orthonormal basis. This local triad in the undeformed configuration denoted by {e 01 , e 02 , e 03 } varies as a function of h 1 if the undeformed centroid line is not straight. The plane containing C is spanned by vectors {e 1 , e 2 } while e 3 denotes its normal (unit) vector. As shear deformation is possible, C is not necessarily normal to the beam's centroid line, i.e. ∂x c ∂h 1 × e 3 = 0. As no cross-sectional deformation occurs furthermore, vector e i can be obtained by a rigid rotation of its associated vector in the undeformed configuration, e 0i , according to: is the group of orthogonal transformations [27].
For further use, we define two vectors tangent to the surface at x(h) as τ 1 = ∂x ∂h 1 and τ 2 = ∂x ∂h 2 . In general, τ 1 and τ 2 are not necessarily orthogonal to each other. The unit vector normal to the surface at the same surface point is defined as follows: (2.5)

The contact scheme
For shear undeformable beams with circular sections, simplified contact kinematics can be found [23,34,36]. Simplifications are also possible if only one of the two beams in contact has a circular cross-section [13]. However, if both beams possess non-circular cross-sections and shear deformation is accounted for, a different approach is needed as the surfaces of the beams cannot be deduced from their centroid line alone. The contact scheme presented here seeks to treat contact interactions between beams if the signed distance function between the beams' surfaces does not possess a unique minimum (e.g. if the beams are parallel to each other or wrapped around each other). This is in contrast to the schemes of [7,8] in which penetration is prevented by forces between the closest pair of surface points.

Projection problem and signed distance function
Let us consider beams B I and B J , and their respective surfaces ∂B I and ∂B J which may be colliding. In our scheme, for a given cross section along the centroid line of B I , we determine if it penetrates ∂B J and if so, by what amount. B I and B J thus have a different role and in order to distinguish them, we call B I the slave and B J the master. C denotes a cross-section of B I and ∂C its perimeter. C is attached to centroid point x 1 c (h C1 ) (see Fig. 2). The outward pointing normal (unit vector) at surface point x I (h I ) ∈ ∂C with h I = h C1 , h I2 T is expressed as: (2.6) We introduce the following gap vector: In the following, we present two possibilities to quantify the amount of penetration between ∂C and ∂B J . The first approach is based on a constraint minimisation of a scalar function, whereas the second one relies on a set of equations to be solved. Higher order surface derivatives are expected for the first approach which are not necessarily well defined on the surface of a body after discretization with beam finite elements. The first approach also yields complicated expressions that translate into more (and hence less efficient) code.
Approach I: pair of surface points determined from the minimisation of an objective function One approach to determine if ∂C penetrates ∂B J and if so, by what amount, is to determine the minimum of scalar function g·n I . In order to prevent that objective function g · n I has an infinite number of minima, we only consider points x I ∈ C and x J ∈ ∂B J that meet constraint c: Here and in the following, a bar over a quantity indicates that it is evaluated at the solution of the local problem. Unit vectorñ I in Eq. (2.9) is defined as:ñ and denotes the unit vector normal to the plane spanned by surface vectors n I and τ I 2 (see Fig. 2). At the solution of Eq. (2.8), the gap vector is colinear with n I but not with n J , which differs from conventional master-slave approaches like the node-to-surface algorithm [18]. In this way, the first and second-order derivatives of Eq. (2.8) are shorter. The reason is that n J depends on two variables, h J1 and h J2 , whereas n I only depends on h I2 . Note also that if g is not aligned with n I at the solution of Eq. (2.8), the measure of penetration [see Eq. (2.15)] is not measured in the direction of the normal to ∂B I . This yields non-physical components of the contact traction vector when contact constraints are regularized (see Sect.

2.2.2).
If only the constraint in Eq. (2.9) is present, four solutions are possible. This is graphically illustrated in Fig. 3. To prevent this, the last two constraints are added, but they can only be verified once the minimisation problem of Eq. (2.8) is solved.
We solve the minimisation problem using the interior extremum theorem, in which only the first constraint in Eq. (2.9) is incorporated via a Lagrange multiplier: where q = h I2 , h J1 , h J2 , λ denotes the variables that we solve for, which consist of three surface parameters and Lagrange multiplier λ ∈ R. To solve Eq. (2.13), we apply Newton's method for which we linearise residual f in Eq. (2.13) which requires the following Jacobian: Once the solution of Eq. (2.13) is found, the amount of penetration, which is also called the normal gap or signed distance function, is given by: Equation (2.15) entails that g N is only negative in case of penetration, as long as the constraints in Eqs. (2.10) and (2.11) are met. The problem with Jacobian H in Eq. (2.14) is the presence of third-order surface derivatives, since n I andñ I are based on first-order surface derivatives. This comes with two disadvantages: • Third-order surface derivatives are not necessarily smooth, in particular in the spatially discretized setting (see Sect. 3). This may impair the convergence of the aforementioned Newton's procedure, which for instance yields solutions for which the constraints of (2.10) and (2.11) do not hold (see also Fig. 3). • f and H are inefficient to compute due to their complicated expressions.
Quantifying the amount of penetration in terms of a residual form [i.e. a system of non-linear equations that replaces Eq. (2.13)] in which lower-order surface derivatives are present, may therefore be computationally advantageous.
Approach II: Pair of surface points determined from a residual form As an alternative to Eq. (2.13), one can formulate a problem in which no minimization is to be performed explicitly, but where we start from some residual form. The vector equations we solve for are: such that: At the solution of the local problem,ḡ = g N where g N has been defined in Eq. (2.15). The obvious problem with Eq. (2.16) is that four unknowns are present in a system of three non-linear equations. We therefore add one more equation: with: is motivated by the fact that at the solution sought, both projected normal vectors n I and n J p must be orthogonal to the vector tangent to the (master) surface, τ J 2 . Also in Eq. (2.19), a 1 is used to make equations from f 1 and f 2 consistent in terms of units. The new set of equations to solve for thus reads: A pair of surface points that abides Eq. (2.21) is shown in Fig. 4. Equation (2.21) is nonlinear and is solved using Newton's method for which the following Jacobian is required: The Jacobian however includes only second-order surface derivatives, whereas the Jacobian of Eq. (2.14) includes thirdorder surface derivatives. Thanks to the lower order of surface derivatives of Approach II, it is expected that Approach II is more robust than Approach I and that f and H are faster to compute (see Sect. 4.1 for a comparison between the two approaches).

Frictionless contact conditions and penalty regularization
The unilateral contact conditions are for both approaches given by: where T N denotes the nominal normal contact traction, and g N is given by Eq. (2.15). If a penalty formulation is used, contact traction T N , acting at the pair of surface points used to measure penetration, is given by: where N > 0 denotes a user-defined parameter, usually referred to as the penalty stiffness, and • denote the Macaulay brackets. Nominal contact traction vector T acting on the section acts in the normal direction to the slave surface at surface point x I such that: The corresponding virtual work of the contact force is then given by: 26) where δg N denotes the variation of the normal gap g N with respect to all involved kinematic variables: the displacement vector components and rotations of the two beams in contact length of X I c and dh I1 denotes the differential of h I1 .

Contact virtual work
The virtual work equation for the two-body system B I , B J including contact interactions reads: where δ B i denotes the internal and external virtual work of B i (excluding contact interactions). p and δ p are admissible functions of the trial and test spaces, respectively. δ c denotes the virtual work of the contact forces between B I and B J . As stated above, cases of interest in this contribution are configurations in which contact interactions arise over a finite length along the beams in contact. T N , introduced in Eq. (2.24), can be seen as a length-specific contact traction because it only accounts for a single (slave) section. The virtual work of all penetrated sections follows from the integration of the virtual work dδ c along the centroid line of B I as follows: where h 1I L and h 1I U denote the lower and upper bound of the integral along h J1 . Note that only penetrated sections contribute to δ c [see Eq. (2. 26)]. Now that the contact kinematics and the virtual work of the contact forces are defined for the space continuous problem, they must be spatially discretized.

Spatial discretization
The spatial discretization method employed in this contribution is the FE Method. Beam B is now discretized with beam FEs (BFEs). When needed, subscript K denotes the index of a node. A brief explanation of the interpolation of the surface of a BFE is given in the following. This will serve as the basis for the discretized contact formulation.

Interpolation of the surface
For a BFE that is part of the discretization of beam B, the position of a surface point can still be obtained from Eq. (2.2), where the position of the centroid line in the undeformed configuration is given by: X K denotes the reference location of node K and N X K denotes the associated interpolation functions. The centroid line position in the current configuration is given by: whereû K ∈ R 3 denotes the displacement of node K and N u K its displacement interpolation function. Rotation vector θ used to compute using Rodrigues' formula (see [15]) is interpolated as follows: whereθ K ∈ R 3 denotes the nodal rotation vector of node K and N θ K the associated interpolation function. For each BFE, n u nodes are used to interpolate displacements, n X nodes are used to interpolate positions in the undeformed configuration and n θ nodes are used to interpolate rotations. Depending on the beam formulation employed n u , n X and n θ can differ 1 .
Not every type of BFE yields a C 0 -continuous surface of the discretized beam [1,23]. This obviously makes contact difficult to formulate. For instance, two-node Geometrically Exact BFEs [10] are employed in Sect. 4. The surface of beams discretized with such BFEs is only C 0 -continuous if the beam is initially straight. This poor continuity of the surface comes with several disadvantages when contact is considered [35]. This is why for each beam, the smoothing technique of [21] is employed which provides a surface close to discretized beam's surface, but C 1 -continuous. Such a smoothed surface is constructed as an assembly of consecutive patches that are based on the original discretized geometry. The resulting surface continuity reduces the risk of abrupt changes in the direction of the contact force between subsequent global iterations of the Newton-Raphson scheme. This is however not elaborated here in order to keep the focus on the contact formulation. However, the numerical examples of Sect. 4 only employ smoothed surfaces of the discretized beams.

Discretized contact weak form and linearization.
Let us now consider two BFEs, M and N which are part of the discretizations of B I and B J , respectively. Surfaces of these elements are denoted by ∂S M and ∂S N . We store the nodal variables of both beams in vec- The surfaces of both elements are assumed to collide. Consequently, force vector r c and stiffness matrix K c associated with this contact must be computed and assembled in the global force residual and global stiffness matrix. Different approaches to compute these entities exist: 1. In a single-pass approach, M is assumed to be the slave and N the master. The discretized version of Eq.
Traction continuity is preserved locally, but the approach introduces a bias, i.e. the choice of which body is the slave and which one is the master influences the results. 2. In a double-pass approach, a single-pass procedure as in 1 is performed twice. First, M acts as the slave and N as the master, then their roles are inverted, i.e. N becomes the slave and M the master. This results in a unbiased approach. Besides doubling the computational costs associated with contact, the problem might also become over-constrained [29]. 3. The "two-half-pass" algorithm was first introduced in [28] and [29] for the frictionless and frictional cases, respectively. During the first half-pass where M is the slave and N the master, contact traction vector t = t M only acts on ∂S M and no traction vector affects ∂S N . During the second half-pass, N acts as the slave and M as the master. The traction vector for the second halfpass is t N (and not −t M ) and only acts on ∂S N . This leads to additional computational effort relative to the single-pass algorithm, but entails an unbiased treatment of contact. Over-constraining is less likely to occur than with a double-pass scheme, and [28,29] have shown an increased robustness for the two-half-pass algorithm relative to the other two options, if applied to standard FE models.
The procedures necessary to obtain the contact residual and contact stiffness for a single-pass (1) and a two-half-pass (3) approach are described next. The double-pass algorithm (2) can be trivially obtained by performing a single-pass procedure for a second time after inverting the roles of the slave and the master. We therefore do not detail it in the following. We now assume for simplicity that all sections attached to integration points along the centroid line of M have their projection according to Eq.

Single pass algorithm
Equation (2.28) specialized to the contact between two BFEs can be written as: where η ∈ [−1, 1] denotes the centroid point coordinate in the parameter space and J = ∂x I c ∂η . Eq. (3.5) is a quadrature where n M I P integration points are used. The weight and coordinates of the kth integration point is denoted by w k and η k , respectively. The normal gap measured at this integration point is denoted by g N k , and J k = J (h M1 (η k )).
To find contact-residual r c , one must recast Eq. (3.5) into: where δ p denotes the variation of the nodal kinematic variables. To this end, the variation of normal gap δg N k related to the kth integration point must be expressed solely in terms of variations of the kinematic variables such that we can write: b k is obtained from Eq. (A.14) by evaluating all quantities at the surface points obtained from the solution of the local problem at this integration point.
If the formalism introduced in [14,18] is used, the implicit dependency of g N k on variables in p can be included via an exception in the automatic differentiation (AD). In this case, b k can be equivalently obtained from: where A k is obtained from Eq. (A.2) and operatorδ δw denotes differentiation with respect to variables w performed by the automatic differentiation algorithm [14,18].
By summing the contribution of the n M I P integration points, r c can be rewritten as: (3.10) The associated stiffness matrix can also be obtained using the AD procedure for which we write: (3.11)

Two-half-pass algorithm
In the two-half-pass approach, the contact traction vector is computed independently for element B M and B N . To indicate to which half-pass quantities refer to, superscripts "M → N " and "N → M" are employed.
For the first half-pass in which B M is the slave, we write for the contact residual: (3.14) For the second half-pass, the roles of the two bodies are inverted. This entails that for this half-pass the measure of penetration, g N →M N k , is measured from the perimeter of cross-sections of N to ∂S M . The contact residual that corresponds to this half-pass can be expressed as: We can also note that K two-half-pass c , defined as: is not symmetric, unlike K c in Eq. (3.11).

Numerical examples
In the previous sections, we have introduced a scheme to treat non-localized contact between beams discretized with BFEs, possibly shear-deformable and with circular or elliptical sections. In the current section, we investigate its capabilities based on three numerical examples. First, a semi-circular arch is brought in contact with an initially straight beam. In the second example, aligned beams are twisted. Finally, we consider the bending of a wire rope.
For all examples, the contact force is a linear function of g N [see Eq. 2.25]. The penalty stiffness employed is given by: where E denotes the Young's modulus of the beams and ν their Poisson's ratio. The penalty stiffness in Eq. (4.1) corresponds to the (length-specific) apparent stiffness relating penetration and contact force for parallel, isotropic, linear elastic cylinders in Hertz's theory (see [26]). Note also that radius of curvature of contacting surfaces is not present in Eq. (4.1), which is convenient. Note that finite penalty stiffness N can be interpreted as the elastic compliance of, the otherwise rigid, cross-sections.
In the three presented numerical examples, only one integration point per slave patch (see [21]) is used to evaluate r c and K c . The total number of integration points employed is thus low. This reduces the computational cost and also alleviates the risk of over-constraining. -Two different settings for the single-pass approach are studied. In the first one, the top arch is used as slave, while the initially straight beam serves as the slave in the second one (denoted as "inverted singlepass" in the following). -The two-half-pass scheme is also studied.   Figure 6 reports sum of the reaction forces at the fixed nodes of the straight beam. The curves of the single-pass schemes match well. This is however not the case for the curves of the two-half-pass schemes, as oscillations are present. For the two-half-pass scheme furthermore, the sum of the reaction forces at the supports of the two beams does not vanish (see Fig. 7), which suggests violation of actionreaction principle at the contact interface. This difference can be explained by the fact that for two colliding patches, a section of the first patch penetrates the other patch during the first half-pass, whereas this is not necessarily the case for the second half-pass (when the roles of master and slave are inverted). This is illustrated in Fig. 8. Thus, this effect violates Newton's third law. In an ideal situation in the continuum formulation in which there is no penetration, the two-half-pass does not suffer from this lack of balance of forces. This is because in this case, the tangent planes at the contacting points are parallel. In this case and if no quadra- Fig. 8 Example 1: contact tractions applied on the perimeter of sections attached to integration points at t ≈ 1.3 for the two-half-pass approach for the coarsest mesh (mesh A). On the far left, penetration is only detected for one of the two half-passes ture is employed, when a section is found to be penetrating during the first half-pass, a penetration will necessarily be found during the second half-pass (Fig. 9).     This phenomenon is exacerbated because the highest contact tractions are located at the ends of the contact zone while in the center of the arch, the contact tractions are considerably smaller.

Example 1: contact between a semi-circular arch and a straight beam
Similar effects occur for the simple truss structure of Fig. 10 for which a node-to-segment contact scheme is employed. Figure 10a presents the single-pass algorithm. Both bodies are deformed due to the effect of the contact applied to the penetrated node and at its projection on the master truss. Figure 10b presents the results of the twohalf-pass algorithm. Only the bottom structure is deformed because the contact traction is only applied to the slave body at the penetrated node. As no nodes of the top structure are penetrated, no contact traction is applied to it. Hence, it does not deform. This explains the non-vanishing sum of the reaction forces at the supports for the two-half-pass scheme in Fig. 10c. In "Appendix C", the evolution of the global residual is reported for the finest mesh (mesh C) for the three studied integration schemes.

Comparison of computational costs for the local problems
In Sect. 2.2.1, we have introduced two possible sets of equations to quantify penetration. As discussed, if Eq. (2.13) (Approach I) is employed, higher order derivatives are involved yielding a longer code and potentially a longer execution time than if Eq. (2.21) (Approach II) is used. The lengths of the generated codes to compute r c and K c is reported in Table 1 for the single-pass approach.
In order to compare execution times, the numerical example in Fig. 5 is employed (with mesh A). The loading is applied in 1000 increments such that several thousands of projection problems and evaluations of r c and K c are performed. For both approaches, the number of penetrated sections and the number of global iterations to converge is similar (see Fig. 11). This is also the case for the reaction forces (see Fig. 12).
However, the average number of local iterations necessary to converge [such that f < 10 −10 with f given by Eqs. (2.13)

Example 2: twisting of parallel beams
In the second example, we consider three parallel beams (see Figs. 13,14). The orientation of the cross-sections of each beam and their cross-sectional dimensions differ (see Table 2 and Fig. 13). The end sections on both sides of the three beams are rotated by 180 • around the [0, 0, 1] E i axis such that the beams wrap around each other. This loading is applied in 540 increments. The evolution of the number of active (i.e. penetrated) slave sections as well as the number of (global) iterations to converge is shown in Fig. 15. The amount of penetration, g N , along the center line of the slave body is shown in Fig. 16. g N is measured twice for each pair of beams in the two-half-pass algorithm. One can observe the good agreement between the penetration measured for the two-half-pass, albeit the single pass performs just as well at roughly half the computational costs.
Finally, the reaction force and torque computed at the end nodes on one side of the beams are reported in Fig. 17. One can note the good match between the curves corresponding to the single and two-half-pass algorithms. This is in contrast with Example 1 and can be explained by the fact that the penetration measured for the two half-passes is similar for a given pair of patches in contact (see Fig. 16). Thus the erroneous effect in results of Example 1 does not occur in this test case. The deformed structure in the final configuration in Fig. 14b is practically the same for the single or two-half-pass algorithms.

Example 3: bending of a rope
The last example considers a rope-like structure (see Figs. 18, 19a). It consists of seven wires (beams) in total; six are wrapped around a central beam that is initially straight. The centroid lines of the outer beams are parametrised helices [32] in the undeformed configuration. In the undeformed configuration, the beams are slightly detached. The applied BCs are as follows: 1. The sections at the end of the rope undergo a rotation of ± π 4 around the [0, 1, 0] T E i axis. This rotation is applied in 1000 equally spaced increments, 2. the nodes in the center of the rope do not displace but are free to rotate.   All beams are given a Young's modulus of 100 GPa, a Poisson's ratio of 0.3, they are discretized using 89 BFEs and their cross-sections are elliptical with semi-axes of 0.3 m and 0.23 m. The rope has a length of 3 m. A single-pass algorithm is employed. Despite the large number of contact interactions and the substantial deformations, only a few (global) iterations are necessary to converge (see Fig. 20). This is thanks to the proper linearization of the contact residual obtained using the Automatic Differentiation technique. Three-dimensional views of the initial and final configurations of the structure are shown in Fig. 19. The reaction torque and forces are reported in Fig. 21.

Conclusion
In this contribution, we have introduced an efficient methodology to treat non-localized contact between shear deformable beams with circular and elliptical sections. It treats contact  between parallel or almost-parallel beams for which contact interactions cannot be applied at a single pair of surface points. As the presented framework only quantifies penetration once per sampled cross-section, it is considerably more efficient than our previous scheme [21] in which the surfaces need to be sampled.
To numerically approximate the contact area, pairs of surface points are determined along the axial direction of the beams in contact. To this end, several cross-sections are considered along one of the contacting beams. For each section, a measure of penetration is computed. The surface of the beams is explicitly used to formulate the contact kinematics. This makes the framework not only applicable to (sheardeformable and shear-undeformable) beams with circular cross sections, but also to beams with elliptical cross sections. The proposed framework may therefore be considered as an attempt to generalize the "line-to-line" contact scheme of Meier et al. [23] that also treats non-localized contact, but is limited to shear-rigid (Kirchhoff) beams with circular sections.
We have also introduced two approaches to quantify the normal gap. In the first approach it is computed using the minimisation of an objective function. In the second approach it is obtained by solving a set of equations. The first approach requires third order surface derivatives and a lengthier code and is consequently 10-20% slower than the second approach.    The introduced contact framework is a master-slave approach. To overcome the bias introduced in the treatment of contact, we have investigated the "two-half-pass" algorithm and compared it to the classical "single-pass approach". We have observed that the measure of penetration for the two half-passes can be significantly different, which yields oscillations of the reaction forces. The sum of the reaction forces do not vanish anymore in these cases.
The introduced scheme is well suited to detect nonlocalised contact (for beams whose centroid line is parallel or almost parallel), but it is not appropriate to detect localised contact. A scheme that combines our approach for nonlocalized contact and the scheme of Neto et al. [7,8] for localized contact may thus be able to treat more general scenarios. The formulation of such a combined scheme in a variationally consistent manner, as the ABC formulation of Meier et al. [24] for beams with circular cross sections, is however out of the scope of this work and remains for future work.