1 Introduction

Many engineering materials such as paper, fabrics or open-cell 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. Thus, special formulations dedicated to beams are developed [5, 6, 13, 19, 20, 23, 24, 34, 36].

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.

2 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 space-continuous setting. This framework will be adapted to the finite element method in Sect. 3.

2.1 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 \({\mathcal {B}}\) and the parametrization of its surface \({\mathbf {x}}(\underline{h})\in {\mathbb {R}}^{3},\underline{h}=\left[ h^{1},h^{2}\right] \in \left[ 0,L\right] \times \left[ 0,2\pi \right] \in {\mathbb {R}}^{2},\) where \(h^{1}\) and \(h^{2}\) are longitudinal and circumferential parameters, respectively. L denotes the length of the centroid line of \({\mathcal {B}}\) in the undeformed configuration. The current centroid line position is obtained by adding displacement \({\mathbf {u}}\in {\mathbb {R}}^{3}\) to the initial centroid line location, \({\mathbf {X}}_{c}\in {\mathbb {R}}^{3}\):

$$\begin{aligned} {\mathbf {x}}_{c}={\mathbf {X}}_{c}+{\mathbf {u}}. \end{aligned}$$
(2.1)
Fig. 1
figure 1

An elliptical cross-section in its undeformed (grey) and current configuration (cyan). Centroid lines in both configurations are presented as dashed lines. (Color figure online)

Full size image

As cross-sections are rigid and remain planar, \({\mathbf {x}}(\underline{h})\) can be obtained by adding the location vector of a centroid point, \({\mathbf {x}}_{c}(h^{1})\in {\mathbb {R}}^{3}\), to a vector,\({\mathbf {v}}({\underline{h}})\in {\mathbb {R}}^{3}\), that lies in the plane of cross-section \({\mathcal {C}}\) that is attached to \({\mathbf {x}}_{c}\) (see Fig. 1):

$$\begin{aligned} {\mathbf {x}}(\underline{h})={\mathbf {x}}_{c}(h^{1})+{\mathbf {v}}(\underline{h}). \end{aligned}$$
(2.2)

For elliptical sections, \({\mathbf {v}}\) can be expressed as follows:

$$\begin{aligned} {\mathbf {{\mathbf {v}}}}(\underline{h})=a\;\text {cos}(h^{2})\;{\mathbf {e}}_{1}(h^{1})+b\;\text {sin}(h^{2})\;{\mathbf {e}}_{2}(h^{1}), \end{aligned}$$
(2.3)

wherea 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 \(\{{\mathbf {e}}_{1},{\mathbf {e}}_{2},{\mathbf {e}}_{3}\}\) attached to \({\mathcal {C}}\) forms an orthonormal basis. This local triad in the undeformed configuration denoted by \(\{{\mathbf {e}}_{01},{\mathbf {e}}_{02},{\mathbf {e}}_{03}\}\) varies as a function of \(h^{1}\) if the undeformed centroid line is not straight. The plane containing \({\mathcal {C}}\) is spanned by vectors \(\left\{ {\mathbf {e}}_{1},{\mathbf {e}}_{2}\right\} \) while \({\mathbf {e}}_{3}\) denotes its normal (unit) vector. As shear deformation is possible, \({\mathcal {C}}\) is not necessarily normal to the beam’s centroid line, \(\textit{i.e.}\) \(\frac{\partial {\mathbf {x}}_{c}}{\partial h^{1}}\times {\mathbf {{\mathbf {e}}}}_{3}\ne {{\mathbf {0}}}\).

As no cross-sectional deformation occurs furthermore, vector \({\mathbf {{\mathbf {e}}}}_{i}\) can be obtained by a rigid rotation of its associated vector in the undeformed configuration, \({\mathbf {{\mathbf {e}}}}_{0i}\), according to:

$$\begin{aligned} {\mathbf {{\mathbf {e}}}}_{i}={\mathbf {{\varvec{\Lambda }}}}(h^{1}){\mathbf {{\mathbf {e}}}}_{0i}. \end{aligned}$$
(2.4)

\({\mathbf {{\varvec{\Lambda }}}}(h^{1})\in SO(3)\) denotes a rotation tensor, where SO(3) is the group of orthogonal transformations [27].

For further use, we define two vectors tangent to the surface at \({\mathbf {x}}(\underline{h})\) as \({\mathbf {\varvec{\tau }}}_{1}=\frac{\partial {\mathbf {x}}}{\partial h^{1}}\) and \({\mathbf {\varvec{\tau }}}_{2}=\frac{\partial {\mathbf {x}}}{\partial h^{2}}\). In general, \(\varvec{\tau }_{1}\) and \({\mathbf {\varvec{\tau }}}_{2}\) are not necessarily orthogonal to each other. The unit vector normal to the surface at the same surface point is defined as follows:

$$\begin{aligned} {\mathbf {n}}(\underline{h})=\frac{\varvec{\tau }_{1}(\underline{h})\times \varvec{\tau }_{2}(\underline{h})}{\left\| \varvec{\tau }_{1}(\underline{h})\times \varvec{\tau }_{2}(\underline{h})\right\| }. \end{aligned}$$
(2.5)

2.2 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 (\(\textit{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.

2.2.1 Projection problem and signed distance function

Let us consider beams \({\mathcal {B}}^{\mathtt {I}}\) and \({\mathcal {B}}^{\mathtt {J}}\), and their respective surfaces \(\partial {\mathcal {B}}^{\mathtt {I}}\) and \(\partial {\mathcal {B}}^{\mathtt {J}}\) which may be colliding. In our scheme, for a given cross section along the centroid line of \({\mathcal {B}}^{\mathtt {I}}\), we determine if it penetrates \(\partial {\mathcal {B}}^{\mathtt {J}}\) and if so, by what amount. \({\mathcal {B}}^{\mathtt {I}}\) and \({\mathcal {B}}^{\mathtt {J}}\) thus have a different role and in order to distinguish them, we call \({\mathcal {B}}^{\mathtt {I}}\) the slave and \({\mathcal {B}}^{\mathtt {J}}\) the master.

Fig. 2
figure 2

Slave section \({\mathcal {C}}\) (in grey). The plane containing \({\mathbf {x}}^{\mathtt {I}}\) and spanned by vectors \(\varvec{\tau }_{2}^{\mathtt {I}}(h^{{\mathcal {C}}1})\) and \({\mathbf {n}}^{\mathtt {I}}(h^{{\mathcal {C}}1})\) is presented in translucent blue. The unit normal vector to this plane is \(\tilde{{\mathbf {n}}}^{\mathtt {I}}\). (Color figure online)

Full size image

\( {\mathcal {C}}\) denotes a cross-section of \({\mathcal {B}}^{\mathtt {I}}\) and \(\partial {\mathcal {C}}\) its perimeter. \( {\mathcal {C}}\) is attached to centroid point \({\mathbf {x}}_{c}^{1}(h^{{\mathcal {C}}1})\) (see Fig. 2). The outward pointing normal (unit vector) at surface point \({\mathbf {x}}^{\mathtt {I}}(\underline{h}^{\mathtt {I}})\in \partial {\mathcal {C}}\) with \(\underline{h}^{\mathtt {I}}=\left[ h^{{\mathcal {C}}1},h^{{\mathtt {I}}{2}}\right] ^{T}\) is expressed as:

$$\begin{aligned} {\mathbf {n}}^{\mathtt {I}}(\underline{h}^{\mathtt {I}})=\frac{\varvec{\tau }_{1}^{\mathtt {I}}(\underline{h}^{\mathtt {I}})\times \varvec{\tau }_{2}^{\mathtt {I}}(\underline{h}^{\mathtt {I}})}{\left\| \varvec{\tau }_{1}^{\mathtt {I}}(\underline{h}^{\mathtt {I}})\times \varvec{\tau }_{2}^{\mathtt {I}}(\underline{h}^{\mathtt {I}})\right\| }. \end{aligned}$$
(2.6)

We introduce the following gap vector:

$$\begin{aligned} {\mathbf {g}}={\mathbf {x}}^{\mathtt {J}}-{\mathbf {x}}^{\mathtt {I}}, \end{aligned}$$
(2.7)

pointing from a surface point on \(\partial {\mathcal {C}}\) to a surface point on \(\partial {\mathcal {B}}^{\mathtt {J}}\).

In the following, we present two possibilities to quantify the amount of penetration between \(\partial {\mathcal {C}}\) and \(\partial {\mathcal {B}}^{\mathtt {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 \(\partial {\mathcal {C}}\) penetrates \(\partial {\mathcal {B}}^{\mathtt {J}}\) and if so, by what amount, is to determine the minimum of scalar function \({\mathbf {g}}\cdot {\mathbf {n}}^{\mathtt {I}}\). In order to prevent that objective function \({\mathbf {g}}\cdot {\mathbf {n}}^{\mathtt {I}}\) has an infinite number of minima, we only consider points \({\mathbf {x}}^{\mathtt {I}}\in {\mathcal {C}}\) and \({\mathbf {x}}^{\mathtt {J}}\in \partial {\mathcal {B}}^{\mathtt {J}}\) that meet constraint c:

$$\begin{aligned} {[}{\bar{h}}^{\mathtt {I}2},{\bar{h}}^{\mathtt {J}1},{\bar{h}}^{\mathtt {J}2}]=\underset{h^{{\mathtt {I}}{2}},h^{{\mathtt {J}}{1}},h^{{\mathtt {J}}{2}}}{\text {min}}\quad {\mathbf {g}}\cdot {\mathbf {n}}^{\mathtt {I}} \end{aligned}$$
(2.8)

such that

$$\begin{aligned}&c={\mathbf {g}}\cdot \tilde{{\mathbf {n}}}^{\mathtt {I}}=0, \end{aligned}$$
(2.9)
$$\begin{aligned}&{\mathbf {x}}^{\mathtt {I}}(\underline{{\bar{h}}}^{\mathtt {I}}){\mathrm {\;is\;in\;}}{\mathcal {B}}^{\mathtt {J}}{\mathrm {\;if\;}}{\bar{{\mathbf {g}}}}\cdot \bar{{{\mathbf {n}}}}^{1}<0, \end{aligned}$$
(2.10)
$$\begin{aligned}&{\mathbf {x}}^{\mathtt {J}}(\underline{{\bar{h}}}^{\mathtt {J}}){\mathrm {\;is\;in\;}}{\mathcal {B}}^{\mathtt {I}}{\mathrm {\;if\;}}{\bar{{\mathbf {g}}}}\cdot {\bar{{\mathbf {n}}}}^{\mathtt {I}}<0. \end{aligned}$$
(2.11)

Here and in the following, a bar over a quantity indicates that it is evaluated at the solution of the local problem. Unit vector \(\tilde{{\mathbf {n}}}^{\mathtt {I}}\) in Eq. (2.9) is defined as:

$$\begin{aligned} \tilde{{\mathbf {n}}}^{\mathtt {I}}(h^{{\mathcal {C}}1},h^{{\mathtt {I}}{2}})=\frac{\varvec{\tau }_{2}^{\mathtt {I}}\times {\mathbf {n}}^{\mathtt {I}}}{\left\| \varvec{\tau }_{2}^{\mathtt {I}}\times {\mathbf {n}}^{\mathtt {I}}\right\| } \end{aligned}$$
(2.12)

and denotes the unit vector normal to the plane spanned by surface vectors \({\mathbf {n}}^{\mathtt {I}}\) and \(\varvec{\tau }_{2}^{\mathtt {I}}\) (see Fig. 2).

At the solution of Eq. (2.8), the gap vector is colinear with \({\mathbf {n}}^{\mathtt {I}}\) but not with \({\mathbf {n}}^{\mathtt {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 \({\mathbf {n}}^{\mathtt {J}}\) depends on two variables, \(h^{{\mathtt {J}}{1}}\) and \(h^{{\mathtt {J}}{2}}\), whereas \({\mathbf {n}}^{\mathtt {I}}\) only depends on \(h^{{\mathtt {I}}{2}}\). Note also that if \({\mathbf {g}}\) is not aligned with \({\mathbf {n}}^{\mathtt {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 \(\partial {\mathcal {B}}^{\mathtt {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:

$$\begin{aligned} \underline{f}(\underline{{\bar{q}}})=\frac{\partial }{\partial \underline{q}}\left( {\mathbf {g}}\cdot {\mathbf {n}}^{\mathtt {I}}-\lambda c\right) \Big |_{\underline{q}=\underline{{\bar{q}}}}={\underline{0}}, \end{aligned}$$
(2.13)

where \(\underline{q}=\left[ h^{{\mathtt {I}}{2}},h^{{\mathtt {J}}{1}},h^{{\mathtt {J}}{2}},\lambda \right] \) denotes the variables that we solve for, which consist of three surface parameters and Lagrange multiplier \(\lambda \in {\mathbb {R}}\).

Fig. 3
figure 3

Approach I: planar view of the different solutions of the projection problem in Eq. (2.13) in the case of parallel straight beams in contact. a Desired solution;b–d: Undesired solutions for which the constraints in Eqs. (2.10) and (2.11) do not hold. The undesired solutions can be obtained if a poor first guess is provided to the non-linear solver. This problem comes from the cyclic property of the cylindrical coordinate system of the cross-sections [13]

Full size image

To solve Eq. (2.13), we apply Newton’s method for which we linearise residual \(\underline{f}\) in Eq. (2.13) which requires the following Jacobian:

$$\begin{aligned} \underline{{\underline{H}}}(\underline{q})=\frac{\partial \underline{f}}{\partial \underline{q}}. \end{aligned}$$
(2.14)

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:

$$\begin{aligned} g_{N}=(\bar{{\mathbf {x}}}^{\mathtt {J}}-\bar{{\mathbf {x}}}^{\mathtt {I}} )\cdot {\bar{{\mathbf {n}}}}^{\mathtt {I}}={\bar{{\mathbf {g}}}}\cdot {\bar{{\mathbf {n}}}}^{\mathtt {I}}, \end{aligned}$$
(2.15)

where \(\bar{{\mathbf {x}}}^{\mathtt {I}} ={\mathbf {x}}^{\mathtt {I}}(h^{{\mathcal {C}}1},{\bar{h}}^{\mathtt {I}2})\) and \(\bar{{\mathbf {x}}}^{\mathtt {J}}={\mathbf {x}}^{\mathtt {J}}({\bar{h}}^{\mathtt {J}1},{\bar{h}}^{\mathtt {J}2})\). 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 \(\underline{{\underline{H}}}\) in Eq. (2.14) is the presence of third-order surface derivatives, since \({\mathbf {n}}^{\mathtt {I}}\) and \(\tilde{{\mathbf {n}}}^{\mathtt {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).

  • \(\underline{f}\) and \(\underline{{\underline{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.

Fig. 4
figure 4

Approach II: solution sought of the projection problem in Eq. (2.21) illustrated for one slave section \({\mathcal {C}}\). The plane of (unit) normal vector \(\tilde{{\mathbf {n}}}^{\mathtt {I}}\) spanned by vectors \(\varvec{\tau }_{2}^{\mathtt {I}}\) and \({\mathbf {n}}^{\mathtt {I}}\) is again shown in translucent blue. The pair of red surface points corresponding to the solution of Eq. (2.16) lies on this plane. At the solution, vectors \({\mathbf {n}}^{\mathtt {I}}\) and \({\mathbf {n}}^{\mathtt {J}p}\) point in opposite directions and are both orthogonal to \(\varvec{\tau }_{2}^{\mathtt {J}}\), thus verifying Eq. (2.19)

Full size image

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:

$$\begin{aligned} {{\mathbf {f}}}_{1}(\underline{{\bar{q}}})=\bar{{\mathbf {x}}}^{\mathtt {J}}-\bar{{\mathbf {x}}}^{\mathtt {I}} -{\bar{g}}{\bar{{\mathbf {n}}}}^{\mathtt {I}}={{\mathbf {0}}}, \end{aligned}$$
(2.16)

such that:

$$\begin{aligned}&\bar{{\mathbf {x}}}^{\mathtt {I}} (\underline{{\bar{h}}}^{\mathtt {I}}){\mathrm {\;is\;in\;}}{\mathcal {B}}^{\mathtt {J}}{\mathrm {\;if\;}}{\bar{g}}<0, \end{aligned}$$
(2.17)
$$\begin{aligned}&\bar{{\mathbf {x}}}^{\mathtt {J}}(\underline{{\bar{h}}}^{\mathtt {J}}){\mathrm {\;is\;in\;}}{\mathcal {B}}^{\mathtt {I}}{\mathrm {\;if\;}}{\bar{g}}<0. \end{aligned}$$
(2.18)

where \(\underline{q}=\left[ h^{{\mathtt {I}}{2}},h^{{\mathtt {J}}{1}},h^{{\mathtt {J}}{2}},g\right] \). At the solution of the local problem, \({\bar{g}}=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:

$$\begin{aligned} f_{2}(\underline{{\bar{q}}})=a^{1}({\mathbf {n}}^{\mathtt {I}}+{\mathbf {n}}^{\mathtt {J}p})\cdot \varvec{\tau }_{2}^{\mathtt {J}}=0, \end{aligned}$$
(2.19)

with:

$$\begin{aligned} {\mathbf {n}}^{\mathtt {J}p}=\frac{{\mathbf {n}}^{\mathtt {J}}-({\mathbf {n}}^{\mathtt {J}}{}\cdot \tilde{{\mathbf {n}}}^{\mathtt {I}})\tilde{{\mathbf {n}}}^{\mathtt {I}}}{\left\| {\mathbf {n}}^{\mathtt {J}}-({\mathbf {n}}^{\mathtt {J}}{}\cdot \tilde{{\mathbf {n}}}^{\mathtt {I}})\tilde{{\mathbf {n}}}^{\mathtt {I}}\right\| }. \end{aligned}$$
(2.20)

Equation (2.19) is motivated by the fact that at the solution sought, both projected normal vectors \({\mathbf {n}}^{\mathtt {I}}\) and \({\mathbf {n}}^{\mathtt {J}p}\) must be orthogonal to the vector tangent to the (master) surface, \(\varvec{\tau }_{2}^{\mathtt {J}}\). Also in Eq. (2.19), \(a^{1}\) is used to make equations from \({{\mathbf {f}}}_{1}\) and \(f_{2}\) consistent in terms of units. The new set of equations to solve for thus reads:

$$\begin{aligned} \underline{f}({\bar{\underline{q}}})=\left[ {{\mathbf {f}}}_{1},f_{2}\right] ^{T}={\underline{0}}. \end{aligned}$$
(2.21)

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:

$$\begin{aligned} \underline{{\underline{H}}}(\underline{q})=\frac{\partial \underline{f}}{\partial \underline{q}}. \end{aligned}$$
(2.22)

The Jacobian however includes only second-order surface derivatives, whereas the Jacobian of Eq. (2.14) includes third-order 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 \(\underline{f}\) and \(\underline{{\underline{H}}}\) are faster to compute (see Sect. 4.1 for a comparison between the two approaches).

2.2.2 Frictionless contact conditions and penalty regularization

The unilateral contact conditions are for both approaches given by:

$$\begin{aligned} g_{N}\ge 0,\quad T_{N}\le 0,\quad g_{N}T_{N}=0, \end{aligned}$$
(2.23)

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:

$$\begin{aligned} T_{N}=-\epsilon _{N}\left\langle -g_{N}\right\rangle \end{aligned}$$
(2.24)

where \(\epsilon _{N}>0\) denotes a user-defined parameter, usually referred to as the penalty stiffness, and \(\left\langle \bullet \right\rangle \) denote the Macaulay brackets.

Nominal contact traction vector \({{\mathbf {T}}}\) acting on the section acts in the normal direction to the slave surface at surface point \(\bar{{\mathbf {x}}}^{\mathtt {I}} \) such that:

$$\begin{aligned} {{\mathbf {T}}}=T_{N}{\bar{{\mathbf {n}}}}^{\mathtt {I}}. \end{aligned}$$
(2.25)

The corresponding virtual work of the contact force is then given by:

$$\begin{aligned} d\delta \Pi _{c}=T_{N}\delta g_{N}dL^{{\mathcal {B}}^{\mathtt {I}}}=-\epsilon _{N}\left\langle -g_{N}\right\rangle \delta g_{N}\left\| \frac{\partial {\mathbf {X}}_{c}^{\mathtt {I}}}{\partial h^{{\mathtt {I}}{1}}}\right\| dh^{{\mathtt {I}}{1}},\nonumber \\ \end{aligned}$$
(2.26)

where \(\delta 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 gathered in \(\underline{p}=\left[ \underline{p}^{\mathtt {I}},\underline{p}^{\mathtt {J}}\right] ^{T}\). A detailed derivation of \(\delta g_{N}\) is presented in “Appendix A”. \(dL^{{\mathcal {B}}^{\mathtt {I}}}\) denotes an infinitesimal length of \({\mathbf {X}}_{c}^{\mathtt {I}}\) and \(dh^{{\mathtt {I}}{1}}\) denotes the differential of \(h^{{\mathtt {I}}{1}}\).

2.2.3 Contact virtual work

The virtual work equation for the two-body system \({\mathcal {B}}^{\mathtt {I}},{\mathcal {B}}^{\mathtt {J}}\) including contact interactions reads:

$$\begin{aligned} \delta \Pi _{{\mathcal {B}}^{\mathtt {I}}}(\underline{p}^{\mathtt {I}},\delta \underline{p}^{\mathtt {I}})+\delta \Pi _{{\mathcal {B}}^{\mathtt {J}}}(\underline{p}^{\mathtt {J}},\delta \underline{p}^{\mathtt {J}})+\delta \Pi _{c}(\underline{p},\delta \underline{p})=0,\nonumber \\ \end{aligned}$$
(2.27)

where \(\delta \Pi _{{\mathcal {B}}^{i}}\) denotes the internal and external virtual work of \({\mathcal {B}}^{i}\) (excluding contact interactions). \(\underline{p}\) and \(\delta \underline{p}\) are admissible functions of the trial and test spaces, respectively. \(\delta \Pi _{c}\) denotes the virtual work of the contact forces between \({\mathcal {B}}^{\mathtt {I}}\) and \({\mathcal {B}}^{\mathtt {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\delta \Pi _{c}\) along the centroid line of \({\mathcal {B}}^{\mathtt {I}}\) as follows:

$$\begin{aligned} \delta \Pi _{c}=\int _{h_{L}^{1\mathtt {I}}}^{h_{U}^{1\mathtt {I}}}d\delta \Pi _{c}, \end{aligned}$$
(2.28)

where \(h_{L}^{1\mathtt {I}}\) and \(h_{U}^{1\mathtt {I}}\) denote the lower and upper bound of the integral along \(h^{{\mathtt {J}}{1}}\). Note that only penetrated sections contribute to \(\delta \Pi _{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.

3 Spatial discretization

The spatial discretization method employed in this contribution is the FE Method. Beam \({\mathcal {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.

3.1 Interpolation of the surface

For a BFE that is part of the discretization of beam \({\mathcal {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:

$$\begin{aligned} {{\mathbf {X}}}_{c}(h_{1})=\sum _{K=1}^{n_{X}}N_{K}^{X}(h_{1}){\hat{{\mathbf {X}}}}_{K}. \end{aligned}$$
(3.1)

\({\hat{{\mathbf {X}}}}_{K}\) denotes the reference location of node K and \(N_{K}^{X}\) denotes the associated interpolation functions. The centroid line position in the current configuration is given by:

$$\begin{aligned} {\mathbf {x}}_{c}(h^{1})={\mathbf {X}}_{c}+\sum _{K=1}^{n_{u}}N_{K}^{u}(h^{1}){\hat{{\mathbf {u}}}}_{K}, \end{aligned}$$
(3.2)

where \({\hat{{\mathbf {u}}}}_{K}\in {\mathbb {R}}^{3}\) denotes the displacement of node K and \(N_{K}^{u}\) its displacement interpolation function. Rotation vector \(\varvec{\theta }\) used to compute \({\varvec{\Lambda }}\) using Rodrigues’ formula (see [15]) is interpolated as follows:

$$\begin{aligned} \varvec{\theta }(h^{1})=\sum _{K=1}^{n_{\theta }}N_{K}^{\theta }(h^{1}){\hat{\varvec{\theta }}}_{K}, \end{aligned}$$
(3.3)

where \({\hat{\varvec{\theta }}}_{K}\in {\mathbb {R}}^{3}\) denotes the nodal rotation vector of node K and \(N_{K}^{\theta }\) 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_{\theta }\) nodes are used to interpolate rotations. Depending on the beam formulation employed \(n_{u}\), \(n_{X}\) and \(n_{\theta }\) can differFootnote 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.

3.2 Discretized contact weak form and linearization.

Let us now consider two BFEs, \({\mathcal {M}}\) and \({\mathcal {N}}\) which are part of the discretizations of \({\mathcal {B}}^{\mathtt {I}}\) and \({\mathcal {B}}^{\mathtt {J}}\), respectively. Surfaces of these elements are denoted by \( \partial {\mathcal {S}}^{{\mathcal {M}}}\) and \( \partial {\mathcal {S}}^{{\mathcal {N}}}\). We store the nodal variables of both beams in vectors \(\underline{p}^{{\mathcal {M}}}=\left[ {\hat{{\mathbf {u}}}}_{1}^{{\mathcal {M}}},\dots ,{\hat{{\mathbf {u}}}}_{n_{u}}^{{\mathcal {M}}},{\hat{\varvec{\theta }}}_{1}^{{\mathcal {M}}},\dots ,{\hat{\varvec{\theta }}}_{n_{\theta }}^{{\mathcal {M}}}\right] \) and \(\underline{p}^{{\mathcal {N}}}=\left[ {\hat{{\mathbf {u}}}}_{1}^{{\mathcal {N}}},\dots ,{\hat{{\mathbf {u}}}}_{n_{u}}^{{\mathcal {N}}},{\hat{\varvec{\theta }}}_{1}^{{\mathcal {N}}},\dots ,{\hat{\varvec{\theta }}}_{n_{\theta }}^{{\mathcal {N}}}\right] \).

The surfaces of both elements are assumed to collide. Consequently, force vector \(\underline{r}_{c}\) and stiffness matrix \(\underline{{\underline{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. 1.

    In a single-pass approach, \({\mathcal {M}}\) is assumed to be the slave and \({\mathcal {N}}\) the master. The discretized version of Eq. (2.28) is integrated along \(h^{{\mathcal {M}}1}\). The resulting contact traction vector acting on \(\partial {\mathcal {S}}^{{\mathcal {M}}}\) is \({{\mathbf {t}}}={{\mathbf {t}}}_{{\mathcal {M}}}\) and acts along the normal to the surface of \({\mathcal {M}}\), but not necessarily along the normal to the surface of \({\mathcal {N}}\). On \(\partial {\mathcal {S}}^{{\mathcal {N}}}\), the traction vector acting at the surface point solution of the projection problem \({\mathbf {x}}^{{\mathcal {N}}}(\underline{h}^{{\mathcal {N}}})\) is \(-{{\mathbf {t}}}_{{\mathcal {M}}}\). 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. 2.

    In a double-pass approach, a single-pass procedure as in 1 is performed twice. First, \({\mathcal {M}}\) acts as the slave and \({\mathcal {N}}\) as the master, then their roles are inverted, i.e. \({\mathcal {N}}\) becomes the slave and \({\mathcal {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. 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 \({\mathcal {M}}\) is the slave and \({\mathcal {N}}\) the master, contact traction vector \({{\mathbf {t}}}={{\mathbf {t}}}_{{\mathcal {M}}}\) only acts on \(\partial {\mathcal {S}}^{{\mathcal {M}}}\) and no traction vector affects \(\partial {\mathcal {S}}^{{\mathcal {N}}}\). During the second half-pass, \({\mathcal {N}}\) acts as the slave and \({\mathcal {M}}\) as the master. The traction vector for the second half-pass is \({{\mathbf {t}}}_{{\mathcal {N}}}\) (and not \(-{{\mathbf {t}}}_{{\mathcal {M}}}\)) and only acts on \(\partial {\mathcal {S}}^{{\mathcal {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 \({\mathcal {M}}\) have their projection according to Eq. (2.13) or Eq. (2.21) on \(\partial {\mathcal {S}}^{{\mathcal {N}}}\), and not on another element’s surface. Similarly, it is assumed that points on the perimeter of sections of \({\mathcal {N}}\) have their projections on \(\partial {\mathcal {S}}^{{\mathcal {M}}}\) and not on another element’s surface. In practice, the different sections of a slave element may have their projections on different master elements. Kinematic variables associated with contact between \({\mathcal {M}}\) and \({\mathcal {N}}\) are gathered in vector \(\underline{p}=\left[ \underline{p}^{{\mathcal {M}}},\underline{p}^{{\mathcal {N}}}\right] ^{T}\).

3.2.1 Single pass algorithm

Equation (2.28) specialized to the contact between two BFEs can be written as:

$$\begin{aligned} \delta \Pi _{c}= & {} - \epsilon _{N}\int _{-1}^{1}\left\langle -g_{N}(\eta )\right\rangle \delta g_{N}(\eta )\left\| {\mathcal {J}}(\eta )\right\| d\eta \end{aligned}$$
(3.4)
$$\begin{aligned}&\thickapprox&- \epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {M}}}}w_{k}\left\langle -g_{N}(\eta _{k})\right\rangle \delta g_{N}(\eta _{k})\left\| {\mathcal {J}}(\eta _{k})\right\| , \end{aligned}$$
(3.5)
$$\begin{aligned}&\thickapprox&-\,\epsilon _N\sum _{k}^{n_{IP}^{{\mathcal {M}}}}w_{k}\left\langle -g_{Nk}\right\rangle \delta g_{Nk}\left\| {\mathcal {J}}_{k}\right\| , \end{aligned}$$
(3.6)

where \(\eta \in [-1,1]\) denotes the centroid point coordinate in the parameter space and \( {\mathcal {J}}=\ \frac{\partial {\mathbf {x}}_{c}^{\mathtt {I}}}{\partial \eta }\). Eq. (3.5) is a quadrature where \(n_{IP}^{{\mathcal {M}}}\) integration points are used. The weight and coordinates of the kth integration point is denoted by \(w_{k}\) and \(\eta _{k}\), respectively. The normal gap measured at this integration point is denoted by \(g_{Nk}\), and \({\mathcal {J}}_{k}={\mathcal {J}}(h^{{\mathcal {M}}1}(\eta _{k}))\).

To find contact-residual \(\underline{r}_{c}\), one must recast Eq. (3.5) into:

$$\begin{aligned} \delta \Pi _{c}=\delta \underline{p}\,^{T}\underline{r}_{c}, \end{aligned}$$
(3.7)

where \(\delta \underline{p}\) denotes the variation of the nodal kinematic variables. To this end, the variation of normal gap \(\delta g_{Nk}\) related to the kth integration point must be expressed solely in terms of variations of the kinematic variables such that we can write:

$$\begin{aligned} \delta g_{Nk}=\underline{b}_{k}^{T}\delta \underline{p}. \end{aligned}$$
(3.8)

\(\underline{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_{Nk}\) on variables in \(\underline{p}\) can be included via an exception in the automatic differentiation (AD). In this case, \(\underline{b}_{k}\) can be equivalently obtained from:

$$\begin{aligned} \delta g_{Nk}=\left( \frac{{\hat{\delta }}g_{Nk}}{{\hat{\delta }}\underline{p}}\Big |_{\frac{\partial \underline{{\bar{q}}}_{k}}{\partial \underline{p}}=\underline{{\underline{A}}}_{k}}\right) ^{T}\delta \underline{p}=\underline{b}_{k}^{T}\delta \underline{p}, \end{aligned}$$
(3.9)

where \(\underline{{\underline{A}}}_{k}\) is obtained from Eq. (A.2) and operator \(\frac{{\hat{\delta }}}{{\hat{\delta }}{\underline{w}}}\) denotes differentiation with respect to variables \({\underline{w}}\) performed by the automatic differentiation algorithm [14, 18].

By summing the contribution of the \(n_{IP}^{{\mathcal {M}}}\) integration points, \(\underline{r}_{c}\) can be rewritten as:

$$\begin{aligned} \underline{r}_{c}\left( \underline{p}^{{\mathcal {M}}},{\bar{\underline{q}}}^{{\mathcal {M}}}_{1},...,{\bar{\underline{q}}}_{n_{IP}^{{\mathcal {M}}}}\right) \thickapprox -\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {M}}}}w_{k}\left\langle -g_{Nk}\right\rangle \underline{b}_{k}\left\| {\mathcal {J}}_{k}\right\| .\nonumber \\ \end{aligned}$$
(3.10)

The associated stiffness matrix can also be obtained using the AD procedure for which we write:

$$\begin{aligned} \underline{{\underline{K}}}_{c}= & {} -\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {M}}}}\frac{{\hat{\delta }}}{{\hat{\delta }}\underline{p}}(\underline{r}_{ck})\Big |_{\frac{\partial \underline{{\bar{q}}}_{k}}{\partial \underline{p}} =\underline{{\underline{A}}}_{k}}\nonumber \\= & {} -\epsilon _{N}\left( \sum _{k}^{n_{IP}^{{\mathcal {M}}}}w_{k}\frac{{\hat{\delta }}}{{\hat{\delta }}\underline{p}}\left( \left\langle -g_{Nk}\right\rangle \left\| {\mathcal {J}}_{k}\right\| \underline{b}_{k}\right) \Big |_{\frac{\partial \underline{{\bar{q}}}_{k}}{\partial \underline{p}}=\underline{{\underline{A}}}_{k}}\right) . \nonumber \\ \end{aligned}$$
(3.11)

3.2.2 Two-half-pass algorithm

In the two-half-pass approach, the contact traction vector is computed independently for element \({\mathcal {B}}^{{\mathcal {M}}}\) and \({\mathcal {B}}^{{\mathcal {N}}}\). To indicate to which half-pass quantities refer to, superscripts “\( {\mathcal {M}}\rightarrow {\mathcal {N}}\)” and “\({\mathcal {N}}\rightarrow {\mathcal {M}}\)” are employed.

For the first half-pass in which \({\mathcal {B}}^{{\mathcal {M}}}\) is the slave, we write for the contact residual:

$$\begin{aligned}&\underline{r}_{c}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\left( \underline{p}^{{\mathcal {M}}},\bar{\underline{q}_{1}},...,{\bar{\underline{q}}}_{n_{IP}^{{\mathcal {M}}}}\right) \nonumber \\&\quad \thickapprox -\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {M}}}}w_{k}\left\langle -g_{Nk}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\right\rangle \underline{b}_{k}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\left\| {\mathcal {J}}_{k}\right\| . \end{aligned}$$
(3.12)

The difference between \(\underline{b}_{k}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\) in Eq. (3.12) and \(\underline{b}\) in Eq. (3.8) is that only the kinematic variables of \({\mathcal {B}}^{{\mathcal {M}}}\) are used to construct \(\underline{b}_{k}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\). The reason is that the contact traction is considered to only act on \(\partial {\mathcal {S}}^{{\mathcal {M}}}\). However, every \(g_{Nk}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\) (and consequently \(\underline{r}_{c}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\)) in Eq. (3.12) depends on \(\underline{p}^{{\mathcal {N}}}\) because of the (implicit) dependency of the local problem with respect to variables in \(\underline{p}^{{\mathcal {N}}}\). The linearization of \(\underline{r}_{c}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\) then yields the two following sub-matrices:

$$\begin{aligned}&\underline{{\underline{K}}}_{\mathcal {MM}} \nonumber \\&\quad =-\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {M}}}}w_{k}\frac{{\hat{\delta }}}{{\hat{\delta }}\underline{p}^{{\mathcal {M}}}}\left( \left\langle -g_{Nk}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\right\rangle \underline{b}_{k}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\left\| {\mathcal {J}}_{k}\right\| \right) \Big |_{\frac{\partial \underline{{\bar{q}}}_{k}}{\partial \underline{p}}=\underline{{\underline{A}}}_{k}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}},\nonumber \\ \end{aligned}$$
(3.13)
$$\begin{aligned}&\underline{{\underline{K}}}_{\mathcal {MN}} \nonumber \\&\quad =-\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {M}}}}w_{k}\frac{{\hat{\delta }}}{{\hat{\delta }}\underline{p}^{{\mathcal {N}}}}\left( \left\langle -g_{Nk}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\right\rangle \underline{b}_{k}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}\left\| {\mathcal {J}}_{k}\right\| \right) \Big |_{\frac{\partial \underline{{\bar{q}}}_{k}}{\partial \underline{p}}=\underline{{\underline{A}}}_{k}^{{\mathcal {M}}\rightarrow {\mathcal {N}}}}.\nonumber \\ \end{aligned}$$
(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_{Nk}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\), is measured from the perimeter of cross-sections of \({\mathcal {N}}\) to \(\partial {\mathcal {S}}^{{\mathcal {M}}}\). The contact residual that corresponds to this half-pass can be expressed as:

$$\begin{aligned}&\underline{r}_{c}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}(\underline{p}^{{\mathcal {N}}},\bar{\underline{q}_{1}},...,{\bar{\underline{q}}}_{n_{IP}^{{\mathcal {N}}}}) \nonumber \\&\quad \thickapprox -\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {N}}}}w_{k}\left\langle -g_{Nk}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\right\rangle \underline{b}_{k}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\left\| {\mathcal {J}}_{k}\right\| . \end{aligned}$$
(3.15)
Fig. 5
figure 5

Example 1: schematic presentation of the structure with BCs and deformed configuration at different times of the loading. The configuration at \(t=0\) and \(t=1\) is presented in translucent in the central and right picture, respectively

Full size image

Only the kinematic variables in \(\underline{p}^{{\mathcal {N}}}\) are used to construct \(\underline{b}_{k}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\). Thus, \(\underline{r}_{c}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\) only affects the entries corresponding to the kinematic variables of \({\mathcal {N}}\). The linearization of \(\underline{r}_{c}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\) yields two new sub-matrices:

$$\begin{aligned}&\underline{{\underline{K}}}_{\mathcal {NN}} \nonumber \\&\quad =-\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {N}}}}w_{k}\frac{{\hat{\delta }}}{{\hat{\delta }}\underline{p}^{{\mathcal {N}}}}\left( \left\langle -g_{Nk}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\right\rangle \underline{b}_{k}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\left\| {\mathcal {J}}_{k}\right\| \right) \Big |_{\frac{\partial \underline{{\bar{q}}}_{k}}{\partial \underline{p}}=\underline{{\underline{A}}}_{k}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}},\nonumber \\ \end{aligned}$$
(3.16)
$$\begin{aligned}&\underline{{\underline{K}}}_{\mathcal {NM}} \nonumber \\&\quad =-\epsilon _{N}\sum _{k}^{n_{IP}^{{\mathcal {N}}}}w_{k}\frac{{\hat{\delta }}}{{\hat{\delta }}\underline{p}^{{\mathcal {M}}}}\left( \left\langle -g_{Nk}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\right\rangle \underline{b}_{k}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}\left\| {\mathcal {J}}_{k}\right\| \right) \Big |_{\frac{\partial \underline{{\bar{q}}}_{k}}{\partial \underline{p}}=\underline{{\underline{A}}}_{k}^{{\mathcal {N}}\rightarrow {\mathcal {M}}}}.\nonumber \\ \end{aligned}$$
(3.17)

We can also note that \(\underline{{\underline{K}}}_{c}^{\text {two-half-pass}}\), defined as:

$$\begin{aligned} \underline{{\underline{K}}}_{c}^{\text {two-half-pass}}=\left[ \begin{array}{cc} \underline{{\underline{K}}}_{\mathcal {MM}} &{} \underline{{\underline{K}}}_{\mathcal {MN}}\\ \underline{{\underline{K}}}_{\mathcal {{\mathcal {N}}{\mathcal {M}}}} &{} \underline{{\underline{K}}}_{\mathcal {NN}} \end{array}\right] \end{aligned}$$
(3.18)

is not symmetric, unlike \(\underline{{\underline{K}}}_{c}\) in Eq. (3.11).

4 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:

$$\begin{aligned} \epsilon _{N}=\frac{\pi E}{8(1-\nu ^{2})}, \end{aligned}$$
(4.1)

whereE denotes the Young’s modulus of the beams and \(\nu \) 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 \(\epsilon _{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 \(\underline{r}_{c}\) and \(\underline{{\underline{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.

4.1 Example 1: contact between a semi-circular arch and a straight beam

In the first example, a semi-circular arch is brought in contact with a straight beam. Young’s moduli of the beams are set to \(100\,\text {GPa}\) and their Poisson’s ratios to 0.3. The semi-circular arch has a radius of \(0.9\,\text {m}\) and the straight beam a length of \(2.7\,\text {m}\). Both beams have elliptical sections with principal axis’ lengths of \(0.1\,\text {m}\) and \(0.06\,\text {m}\). For both beams, local basis vector \({\mathbf {E}}_{01}\) points in the \([-1,0,0]_{{\mathbf {E}}_{i}}\) direction in the undeformed configuration (see Fig. 5).

The nodes at the base of the arch are moved vertically towards the bottom beam with a final vertical displacement of \(0.3\,\text {m}\) for \(0<t<1\) in 300 increments, and are then moved horizontally for \(1<t<2\) with a displacement of \(0.1\,\text {m}\) in 300 increments. The rotations of these nodes are also blocked. All kinematics variables at the ends of the straight beam are restrained during the entire simulation.

Fig. 6
figure 6

Example 1: sum of the reaction forces of the fixed nodes of the initially straight beam along the \({\mathbf {E}}_{2}\) and \({\mathbf {E}}_{1}\) directions

Full size image
Fig. 7
figure 7

Example 1: difference between the sum of the reaction forces of the fixed nodes of the straight beam and those of the curved beam

Full size image

Several numerical aspects are investigated, \( \textit{i.e.}\):

  • The number of BFEs: three (uniform) refinements for both beams with a constant ratio of the number of BFEs of both bodies. The mesh of the straight beam is coarser than the one of the arch in order to show a possible influence of the choice of the roles of master and slave. For mesh A, 90 and 60 BFEs are used; for mesh B, 120 and 80 BFEs; and for mesh C, 150 and 100 BFEs.

  • The integration scheme:

    • 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 single-pass” 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 action-reaction 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 violatesNewton’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 quadrature 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).

Fig. 8
figure 8

Example 1: contact tractions applied on the perimeter of sections attached to integration points at \(t\thickapprox 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

Full size image
Fig. 9
figure 9

Example 1: contact traction at \(t=1\) for the two-half-pass approach. Vectors are located at slave surface points with local coordinates solution of the local problem

Full size image
Fig. 10
figure 10

Effect of a single and two-half-pass variants of a contact scheme for truss networks with a node-to-segment contact scheme. The BCs applied to the truss networks are such that the nodes in red are shifted upwards, while the grey node are fixed. The undeformed configurations are presented with solid lines while the deformed configurations are presented with dashed lines. For the single-pass scheme, the bottom structure serves as slave. The trusses have a unit Young’s modulus and a unit cross-sectional area. c shows the vertical component of the sum of the reaction forces for the two contact schemes

Full size image
Table 1 Performance comparison for the codes generated byAceGen
Full size table
Fig. 11
figure 11

Evolution of the number of penetrated sections and number of global iterations necessary to converge for both local schemes. The convergence criterion for the global problem is set as \(\left\| \underline{r}_{g}^{free}\right\| < 10^{-8}\) where \( \left\| \underline{r}_{g}^{free}\right\| \) denote the components of the global residual force vectors that are not subjected to BCs

Full size image
Fig. 12
figure 12

Reaction forces for both local schemes introduced

Full size image

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.

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 two-half-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.

4.1.1 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 \(\underline{r}_{c}\) and \(\underline{{\underline{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 \(\underline{r}_{c}\) and \(\underline{{\underline{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 \(\left\| \underline{f}\right\| <10^{-10}\) with \(\underline{f}\) given by Eqs. (2.13) or (2.21)] as well as the average CPU time to determine \(\underline{r}_{c}\) and \(\underline{{\underline{K_{c}}}}\) are different. The scheme of Eq. (2.21) clearly outperforms the scheme of Eq. (2.13). Note the effect of the simplifications of \(\delta g_{N}\) (see “Appendix A”).

Fig. 13
figure 13

Example 2: view in plane \(({\mathbf {E}}_{1},{\mathbf {E}}_{2})\) of the aligned beams in the undeformed configuration. Beam A is at the bottom right, beam B the bottom left, and beam C on top

Full size image
Fig. 14
figure 14

Example 2: initial and final configuration

Full size image
Table 2 Example 2: properties of the three string of BFEs
Full size table

4.2 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^{\circ }\) around the \(\left[ 0,0,1\right] _{{\mathbf {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.

Fig. 15
figure 15

Example 2: evolution of the number of penetrated sections and number of global iterations necessary to converge. For small twist angles, the beams are not in contact due to the small initial separations between them (see Fig. 14). The convergence criterion for the global problem is set as \(\left\| \underline{r}_{g}^{free}\right\| <10^{-8}\)

Full size image
Fig. 16
figure 16

Example 2: evolution of \(g_{N}\) along the slave centroid line measured for the different couples of beams in contact for two-half-pass and single-pass approaches, for half of the loading and at the end of the loading

Full size image

4.3 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. 1.

    The sections at the end of the rope undergo a rotation of \(\pm \frac{\pi }{4}\) around the \(\left[ 0,1,0\right] ^{T}_{{\mathbf {E}}_{i}}\) axis. This rotation is applied in 1000 equally spaced increments,

  2. 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 \(\text {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\,\text {m}\) and \(0.23\,\text {m}\). The rope has a length of \(3\,\text {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.

Fig. 17
figure 17

Example 2: reaction force and torque obtained by summing the contributions of nodes of the different beams subject to BCs for a single-pass and two-half-pass algorithms. For this example, no oscillations are present for the two-half-pass

Full size image
Fig. 18
figure 18

Example 3: setup. Sections at both ends of the rope are rotated around the \({\mathbf {E}}_{2}\) axis. The displacements of the nodes in the center of the rope are restrained but not their rotations

Full size image
Fig. 19
figure 19

Example 3: reference and final configurations of the structure. The discontinuous surface of the strings of BFEs, that is improper for contact treatment, is shown. On the right, one can observe the effect of the rotations applied to the sections at the ends of the rope

Full size image

5 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 (shear-deformable 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.

Fig. 20
figure 20

Example 3: evolution of the number of penetrated sections and number of global iterations necessary to converge. The convergence criterion for the global problem is set as \(\left\| \underline{r}_{g}^{free}\right\| < 10^{-8}\)

Full size image
Fig. 21
figure 21

Example 3: reaction torque and force at one end of the rope

Full size image

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 non-localised 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 non-localized 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.