Enhanced multiple-point beam-to-beam frictionless contact finite element

In this paper a new contact finite element for beams with circular cross-sections is presented. The element is an enhancement of the previously formulated point-wise beam-to-beam contact finite elements to be used in cases when beams get in contact at very acute angles. In such situation, if beam deformations in the vicinity of the contact zone are taken into account, the contact is not point-wise but it extends to a certain area. To cover such a case in a more realistic way, two additional pairs of contact points are introduced to accompany the original single pair of contact points. The central pair is determined using the orthogonality conditions for the beam axes and the positions of two extra points are defined on one beam axis by a shift of the local co-ordinate. This shift depends on beams geometry and the current angle between tangent vectors at the central contact point. The appropriate kinematic variables for normal contact together with their finite element approximation are derived. Basing on the weak form for normal contact and its linearisation, the tangent stiffness matrix and the residual vector are derived. The new element is tested using author’s computer programs and comparisons with the point-wise contact elements are made.

in [11] and continued in [9,10,12] where contact without and with Coulomb friction for beams of circular and rectangular cross-sections was covered. The further development concerned inclusion of thermal and electric coupling, see [1]. Some subsequent research was also devoted to smoothing procedures for 3D curves representing axes of beams in contact, e.g. [7]. A rigorous approach to the question of point-wise contact was also suggested by Konyukhov and Schweizerhof in [5]. There the authors focused their interest on the closest-point projection procedure, which for the beam-to-beam contact leads to the orthogonality conditions, see [11]. The same authors used their approach to analyse the problem of rope wound around a cylinder and the question of knot-tightening [6]. The latter issue was also considered and solved by Durville in [4].
One of the key assumptions in the majority of those analyses was, that the closest points of two curves can always be uniquely determined, at least locally. In other words, the point-wise contact between beams was assumed. However, such an approach may fail in some special cases, for instance when irregular assemblies of fibre-like objects are considered [3]. Hence, the general approach must include possibilities when the contacting beam-like objects form acute angles and are parallel or conforming, see Fig. 1. The issues related to the existence and uniqueness of the closest points location were discussed in detail in [5].
The new beam-to-beam contact finite element presented in this paper is developed to model in a more precise way the situations, when contact between beams cannot be considered as point-wise. The main assumption for this new element is that, locally, there still exists the unique solution to the closest-point projection equations. Therefore, the presented approach is not suitable for the case of conforming/parallel beams. Such a case requires a completely different approach. In the present proposal, besides the main central pair of contact points found using the orthogonality conditions, two additional pairs of points are involved. Thus, the line along which the contact may occur between two almost parallel or almost conforming beams, is discretised by three contact point pairs, see Fig. 1a. The problem of finding the location of two additional contacts is discussed in Sect. 2. Section 3 presents some important details of the contact formulation for the additional points. The weak form of the contact contribution to the virtual work formulation and its consistent linearisation are presented. All the kinematic variables involved are expressed in terms of position vectors and displacements of contact points. The finite element discretisation of this formulation is introduced in Sect. 4. As a result the explicit form of residual vector and consistent tangent stiffness matrix are derived. They are ready to use in any computer program for static analysis of beams by the FEM. Numerical examples presented in Sect. 5 solved using the author's computer program confirm the effectiveness and advantages of the new beam-to-beam contact element with respect to the previously derived one-point contact element in the cases of beams forming acute angles. Section 6 contains final remarks, conclusions and outlook. Let is be stressed once again, that the presented approach with two additional points is based on the existence and uniqueness of the closest point location procedure. If this fails to determine uniquely the closest points-the central points in the new approach, then the additional points cannot be found, either. That is why the introduced algorithm may be seen as an alternative for the standard point-wise contact in the cases of beams forming acute angles, but not in the cases of parallel or conforming beams or when one beam wraps around the second one.

The additional pairs of contact points
Let us consider contact between two beams denoted by m and s with their axes presented in Fig. 2. The beams have circular cross-sections with radii r m and r s . The presented multiple-point beam-to-beam contact element involves one pair of central contact points C mn − C sn and two additional ones. The location of the central points is found using the orthogonality conditions. This issue was exactly described in [6,11] and as result of the minimisation of the axes distance d N one gets the local co-ordinates of the closest points on  Fig. 3 Additional contact points on beam axes and additional normal gaps each beam-ξ mn and ξ sn (−1 ≤ ξ in ≤ 1). The following formula yields the basic kinematic variable for the normal contact-the gap function The additional contact pairs are found in a two-step way. First, the contact candidate points C snb , C sn f are found on the beam s using the local co-ordinates shifted back (subscript b) and forward (subscript f ) with respect to the local co-ordinate of the central contact point ξ sn as shown in Fig. 3. The definition of the shift ξ s of the local co-ordinate is somewhat arbitrary. Here a result of some geometric considerations with respect to the layout of beam contact zones is adopted. Let us consider a pair of contacting beams. The plan view of the contact zone is presented in Fig. 4, while the cross-section A-A along the beam s is shown in Fig. 5  distance between the centre contact point C sn on the axis of the beam s and the edge of the overlap region depicted in the plan view is Now the distance between a shifted contact point, C snb or C sn f , and the central one C sn is assumed as a fraction 1/k of a s . This yields the following shift value for the local coordinate where l s is the length of the beam element within which the contact point lies. The value of the fraction denominator k can be tuned to a given application. In the numerical tests it was found that k = 3 yields the best results, so finally With the additional points located on the beam s, the orthogonality condition can be applied to find their counterparts C mnb and C mn f on the beam m. The appropriate relations read where dot · denotes the scalar product. Generally, the conditions (7) and (8) represent non-linear equations, which can be solved iteratively using the Newton-Raphson method. With a finite element discretisation adopted, the co-ordinate of element mid-point ξ m = 0 can be used as the starting value. Following the lines of the solution for the central points with a set of orthogonality conditions (see [11]) with the use of the following expressions for the position vectors and their derivatives at the additional points x mnb,m = u mnb,m + x mnb,mm ξ mnb x mn f,m = u mn f,m + x mn f,mm ξ mn f the iterative updates for the local co-ordinates of the points C mnb and C mn f are obtained in the form Having located the additional contact points at both beams m and s two additional normal gap functions can be defined in a way analogous to the normal gap at the central point (1). This yields where the distances between the contact points are With three gap functions given by (2) and (13) the normal contact formulation for three separate contact point pairs can be formulated. The active set strategy can be adopted and the required contact search algorithm, described in a detailed way in [8] can be effectively applied here, too. This allows to express the weak form of the contact contribution in an analogous way as in [10]. Finally, it should be pointed out that the presented procedure changes the character of the beam-to-beam contact. The standard point-wise contact between beams features a specific symmetry-a parity of both contacting beams, making the distinction between the beam m and the beam s fully reversible. The results of analysis do not depend on this distinction. Contrary to that situation, in the present approach, due to the nature of the shifted additional points on the beam s, this parity is no longer present and computation results will depend on the choice of the beam s. In this sense, the new beam-to-beam element will be similar to the nodeto-segment or node-to-surface contact elements used in the modelling of contact between regular solid bodies, where the distinction between slave and master bodies plays a vital role.
Due to the adopted method to find the shifted contact points, as presented in Fig. 4, it is advisable to treat the beam with a smaller cross-section as the beam s. With this provision, the resulting contact points will be more effectively distributed along the real contact zone which is longer along the beam with a smaller cross-section and narrower along the beam with a larger cross-section.

The weak formulation of normal contact with additional points and the kinematic variables
For the contact problem modelled by a set of three contact point pairs the appropriate weak formulation can be used. In this paper the penalty method is applied with the same value of the penalty parameter for all three contacts. Hence, following the derivation from [9,11], the contact term in the weak formulation for a pair of contacting beams can be expressed in the form with three separate contributions for the central and two additional contact points. The subsequent linearisation of the variation (15), required in the general Newton-Raphson iterative scheme used to solve the non-linear contact problem, consists of the respective three contributions, too. One gets It should be pointed out, that a contact check must be applied separately for each from three pairs of contact points: C mnb and C snb , C mn and C sn , C mn f and C sn f . If any of these pairs is inactive, i.e. the gap function for the contact pair exhibits separation, then the appropriate terms in (15) and (16) should be excluded from the formulation. The expressions (15) and (16), which after discretisation yield the residual vector and the tangent stiffness matrix for the contact element, include the kinematic variables related to all three pairs of contact points. The variables for the central contact point have been expressed in terms of displacements in [9,11] and the appropriate formulae will not be repeated in this paper. Here the focus is locked on the derivations corresponding to the variables defined at the additional points.
Due to a different method used to localise them, including the shift of local co-ordinate along the beam s and the orthogonality conditions on the beam m only, there will be some serious differences between the formulation for the pairs of the shifted points and the central pair points.
Variations and linearisations of the gap functions (13) and (14) can be expressed as where the vectors normal to the beam m at the contact points C mnb and C mn f were introduced It is worth to note, that the final expressions (17) and (18) contain the terms-the last ones in the parentheses, which are non-zero due to lack of orthogonality of the normal vectors n b and n f to the tangent to the beam s axis. This is obviously not the case for the central pair of contact points, where the orthogonality holds for both beams. The linearisations of the variations (17) 1 and (18) 1 will consequently involve additional terms, when compared to the corresponding expressions for the central contact pair. After the basic computations, analogous to those presented in [9,11] one gets The kinematic variables (17), (18), (21) and (22) involve variations and linearisations of local co-ordinates on both contacting beams. Additionally, unlike the analogue expressions for the central contact pair [8], the formulae (21) and (22) contain the components (the fifth parentheses) resulting from the lack of orthogonality at the beam s, which include δξ sn f and δξ snb . The required kinematic variables for the beam m can be obtained carrying out the appropriate mathematical operations on the orthogonality conditions (7) and (8). The variations and linearisations take the form In order to evaluate the variables given by (23)-(26), related to the local co-ordinate ξ m , one has to find variations and linearisations and variations of linearisations of local coordinates ξ s at the contact points C sn f and C snb on the beam s. These kinematic variables can be obtained from the relations (3) and (5). From (3) we get The expressions (27)-(29) include the kinematic variables δξ sn , ξ sn and δξ sn related to the central contact point C sn . They were derived in [8] and the appropriate formulae are not given here. The additional computations are required to get the variations and linearisations of the local co-ordinate shift, i.e.: δξ s , ξ s and δξ s . These variables are obtained from the respective operations on (5). To facilitate these derivations, let us first express the variable angle ϕ between the axes of the beams m and s at the central contact points C mn and C sn in terms of the scalar product of respective tangent vectors. The reciprocal of the function sine is given by where c = cosϕ. The scalar product of the tangent vectors takes the form With (30) taken into account the expression (5) is reformulated as Using the chain rule we obtain from (33) The variation and the linearisation of the cosine c is readily obtained from (32) in the following way: Including the auxiliary results one can finally get where Similarly for the linearisations we get where Variations and linearisations of the position vector derivatives with respect to local co-ordinates are taken from [8] δx sn,s = x sn,ss δξ sn + δu sn,s δx mn,m = x mn,mm δξ mn + δu mn,m x sn,s = x sn,ss ξ sn + u sn,s Now the linearisation of (34) has to be found. We get The last expression to be calculated is δc which can be found from the linearisation of (38) as and its components are found separately by the linearisation of (39)-(41). Quite tedious but basic calculations lead to the following results The relation (49) involves linearisations of variations of position vector derivatives with respect to ξ s , which are given in [10] by δx sn,s = δu sn,s + δu sn,ss ξ sn + u sn,ss δξ sn + x sn,sss δξ sn ξ sn + x sn,ss δξ sn δx mn,m = δu mn,m + δu mn,mm ξ mn + u mn,mm δξ mn + x mn,mmm δξ mn ξ mn + x mn,mm δξ mn In this way all the terms in the weak form (15) and its linearisation (16) have been expressed in terms of position vectors and displacements at the contact points C mnb , C snb , C mn , C sn , C mnb , C snb and can now be subjected to the finite element discretisation.

Finite element discretisation -the residual vector and the tangent stiffness matrix
The finite element discretisation of the contact contributions (15) and (16) yields the matrix quantities characterising the contact finite element. In the case of contact existing between two bodies, this element couples the degrees of freedom related to each of these bodies, which in the case of separation would remain uncoupled. In general, the idea of the finite element approach is to express the displacement vector of an arbitrary point on the beam, in particular of the contact points C mnb , C snb , C mn , C sn , C mn f , C sn f , in terms of the nodal displacements assembled in the vectors u M and u S , related to the finite elements discretising the beam m and the beam s, respectively. The specific form of the nodal displacement vectors depends on the method of representation of the beam axis curve. The author's previous work summarised in [8] contains a detailed discussion of this issue, which is also closely related to the all-important problem of curve smoothing. This procedure is necessary to ensure an efficient treatment of large movement of contact points between adjacent finite elements representing beams. In order not to repeat this discussion, let us just assume, that the variations and linearisations of displacements of contact points are given by: The explicit form of all the matrices G, H and M in the above definitions depends on the details of curve representation, which were given in [8].
In the expressions (53)-(61) provision has been made for the general case when the central contact points C mn , C sn do not lie within the same element where the additional contact points C mnb , C snb , C mn f , C sn f are. Thus, the appropriate displacement vectors for the forward and back elements, u M f , u S f , u Mb , u Sb , respectively, were introduced. An exemplary case for the beam m, with different elements M, M b and M f is shown in Fig. 6a. On the contrary, if any of the additional points is located within the same element as the central point, as shown in Fig. 6b, then the appropriate displacement vectors are identical.
In the following only the discretisation of the quantities related to additional points will be given in detail. The dis- cretised variables regarding the central points can be taken directly from [8].
Discretisations of variables including the local coordinates of the additional contact points on the beam s, given by (27)-(29), can be expressed by where the matrices F s and R s are related to the central contact point C sn and were given in [8], while the matrices F s and R s must be evaluated by the discretisation of (34), (35) and (47). Let us introduce the following matrices and parameters to simplify the notation and Now one can define the matrix what allows to evaluate the discretisations (62) and (63). Let us also introduce further matrices The linearisations of variations of u mn,m and u sn,s are taken as in [8] The matrices R m and R s in (64), (69) and (70) are related to the discretisation of δξ mn and δξ sn for the central contact points C mn and C sn which were given in [8]. These expressions allow to determine the following matrices which finally define the matrix necessary to evaluate the discretisations (64). It is worth to note, that those discretisations, even though related to the additional contact points C snb and C sn f moved back and forward with respect to the central contact point C sn , are expressed in terms of nodal displacements of beam elements where central contact points C mn and C sn lie. It is so because of the adopted definition of the shifted points, presented in Sect. 2, where all the variables involved are estimated at the central points.
Discretisations of variables including the local coordinates of the additional contact points on the beam m can be expressed by where the involved matrices yield from the discretisation of components included in (23)-(26). They take the following form where Now the discretisations of variables related to the gap functions (17) and (18) can be determined. For the variation and linearisations we get where The linearisations of the variations of the gap functions can be discretised in the following way In the evaluation of the matrices involved in (85)-(89) discretisations of δu mnb , δu mn f , δu snb and δu sn f are necessary. As discussed in [8] they are required due to the expression of displacements at an arbitrary point of a smoothed curve which is non-linear in terms of nodal displacements. These variables for the additional points can be derived in an analogous way as in [8] for the central points to yield It is also useful to express the normal vectors (19) and (20) by their components: and to introduce further four matrices in order to simplify the notation Now the matrices in (88) and (89) can be given in the following form Having discretised all the kinematic variables one can finally express the residual vector and tangent stiffness matrix for the newly developed multiple point beam-to-beam contact finite element. The weak form (15) after the discretisation yields the residual vector: It contains the component R c for the central pair of contact points, which was given in [8], and the following components related to the additional pairs evaluated from the second and the third component of (15) concerning the additional pairs of contact points. The linearisation of the weak form (16) after the discretisation yields the tangent stiffness matrix: The component matrix K c for the central pair of contact points was given in [8]. The remaining ones related to the additional pairs are evaluated from the components in the second and the third parentheses in (16) and take the form The characteristic vectors and matrices (101) and (103) are ready to embed in any computer program of beam-to-beam contact analysis. An inherent contact search routine, like the one described in [8], determines if a contact contribution has to be switched on. In a case of an active contact pair the appropriate components of (101) and (103) have to be added to the global residual vector or stiffness matrix used within the Newton-Raphson iterative solution scheme.

Numerical examples
In this section two examples of frictionless contact between beams are considered to check the performance of the newly developed three-point contact element. The axes of beams are approximated using Hermite polynomials combined with inscribed curve segments, as discussed in [7,8]. The contact search algorithm from [8] is applied. In the analysed cases beams axes form an acute angle in the plan view to illustrate better the merits of the new approach. Data and results in the examples are given without any specified units but in fact any consistent unit set could be adopted. The computer times for both considered approaches differ. However, the introduction of two additional points, which means the 200 % increase in the number of active constraints, increases the time of computations only by about 15-30 %, depending on the beam layout. The effect is so favourable because the computations of the additional points involve many matrices and vectors, which are identical as the ones used in the computations for the central contact points. Also the majority of new matrices are computed using the same subroutines, which allows for some optimisation of the computer code.

Example 1
In this example two clamped-clamped beams 1 and 2 with the lengths l 1 = 6.000, l 2 = 6.325 and the angle between their axes ϕ = 16.26 • in the plan view, as presented in Fig. 7, are considered. Both beams have circular cross-sections with the radius r = 0.1 and are made of a linearly elastic material characterised with the Young's modulus E = 250 × 10 5 and the Poisson's ratio ν = 0.3. The initial gap separating the beams is g N 0 = 0.001. The beams are forced into contact by the nodal displacements 1 = 2 = 1.0 applied in 60 increments at the ends of beam 2. Each beam is discretised by 10 co-rotational finite elements developed by Crisfield [2].
This example is solved using the currently developed three-point contact element and the standard point-wise contact element. The comparison of the available results is done for various values of the penalty parameter.
The graphs in Fig. 8 present the values of normal contact forces in the case of new element applied. The dotted curves represent the total value of the force while the solid and dashed lines show its components related to the central and additional contact pair, respectively. Due to the symmetry of the beam layout both additional pairs, forward and back, feature the same force. Comparing these graphs one can clearly notice the influence of increasing penalty parameter what makes the contact more stiff. The penetration decreases and the moment when the additional contact pairs become active is more and more delayed.
There is also one more interesting effect related to the increase of the penalty parameter-for its extreme values the central contact point is switched off at the end of the deformation process and the entire contact force is localised in the additional points.
The newly derived element was also compared with the standard point-wise contact case. Figure 9 presents the values of the final total contact force computed using both elements with the solid line corresponding to the three-point contact and the dotted line-to the point-wise contact. It can be clearly seen that the distribution of the force between three contact pairs makes the approximation less dependent on the value of the penalty parameter. The analysis with smaller values of this crucial coefficient yields a better approximation of the contact force. Moreover, the new element allows for larger values of this parameter to be taken into accountthe solution can be found for ε N reaching 10 9 , while the point-wise contact element becomes ill-conditioned already for 3 × 10 6 .
These results indicate that the three-point contact element is an attractive replacement for the point-wise contact element in the cases when beam axes form small acute angles.
The 3D view of the final deformed configuration of beam axes with equal scale for displacements and beam dimensions, obtained for the case of the penalty parameter equal to 10 9 is presented in Fig. 10. Solid lines depict the deformed axes, while the dotted lines-the initial undeformed configuration.

Example 2
In this example two cantilever beams 1 and 2 with the lengths l 1 = l 2 = 6.083 and with the angle between their axes ϕ = 18.92 • in the plan view, as presented in Fig. 11, are considered. Both beams have circular cross-sections with the radius r = 0.1 and are made of a linearly elastic material characterised with the Young's modulus E = 250 × 10 5 and the Poisson's ratio ν = 0.3. The initial gap separating the beams is g N 0 = 0.001. The beams are forced into contact by the nodal displacements 1 = 2 = 0.2 applied in 60 increments simultaneously at the beams tips. Each beam is discretised by 10 co-rotational finite elements developed by Crisfield [2]. The problem was solved using various values of the penalty parameter to check the performance of the new element and compare it to the point-wise contact element.
The graph showing the evolution of the penetration at contact points, computed for the penalty parameter 10 4 , is presented in Fig. 12. The solid line represents the values for the central contact pair, while the dashed line-the additional contact pairs. The contact zone layout is not symmetric in this example but the differences between the contact behaviour in the additional points shifted back and forward are undistinguishable. It can be observed, that initially just the The comparison of results obtained using both contact elements is given in Fig. 13, where the total contact force (104) versus the penalty parameter is shown. The solid line corresponds to the three-point contact and the dotted lineto the point-wise contact. Similarly as in Example 1, the case with the contact force distributed at three points yields better results in terms of convergence of the contact force to the exact value. Its better approximation can be obtained for smaller penalty parameter values what is advantageous due to a better conditioning of the penalty parameter approach.
The 3D view of the final deformed configuration of beam axes with equal scale for displacements and beam dimensions, obtained for the case of the penalty parameter equal to 10 5 , is presented in Fig. 14. Solid lines depict the deformed axes, while the dotted lines-the initial undeformed configuration.

Example 3
In this example four identical beams with the lengths l = 15.31 forming a symmetric assembly where their axes cross in the plan view at the angles ϕ = 21.80 • , as presented in  3. The initial gaps separating the beams are g N 0 = 0.033. The beams are forced into contact by a set of 8 nodal displacements = 1.7 applied in 30 increments simultaneously at all the beams tips. Each beam is discretised by 10 co-rotational finite elements developed by Crisfield [2]. The problem was solved using the penalty parameter ε N = 10 4 .
There are four contact points in this case but because of the symmetry it is enough to give the results for two of them, denoted by A and B in Fig. 15. The total contact force at the point A from the point-wise contact analysis is 146.0, while that from the three-point contact analysis is 147.8. The respective values of the total contact force at the point B are 148.9 and 150.8, respectively. Hence, just like in the previous examples, the new beam-to-beam contact model allows to compute the contact force with a better accuracy. It is also worth to mention, that in the three-point contact model the better result yields from a much smaller violation of the geometry, i.e. the maximal penetration is 0.005415, which represents just 2.7 % of the cross-section depth of the beam. On the other hand, the point-wise contact model yields the penetration of 0.01489, which is equal to 7.4 % of the crosssection depth.
The 3D view of the final deformed configuration of beam axes with equal scale for displacements and beam dimensions is presented in Fig. 16. Solid lines depict the deformed axes, while the dotted lines-the initial undeformed configuration.

Concluding remarks
This paper reports the derivation and numerical tests to assess the performance of a new beam-to-beam contact finite element for the frictionless case. The novelty with respect to the previously developed element is the enhancement of the contact zone modelling in the situation of acute angles between the beam axes. The doubtful point-wise model is replaced by a more realistic distribution of the contact force at three contacting points discretising a real line contact zone. The result is a better representation of the contact phenomenon. The calculated examples point out that this more realistic approach has serious advantages when the most popular penalty parameter method of contact constraint enforcement is used. The new element reduces the localised contact forces, thus yields a better numerical model which is less prone to ill-posedness, provides a wider margin of acceptable penalty parameters   It may also be worth noting that the presence of three pairs of contact points enhances the simple beam-to-beam point-wise contact model, as it introduces a possibility to include the resultant bending moment at the contact zone. This is clearly not the case for the point-wise contact. Hence, the presented approach may be viewed as a model which describes better the 3D nature of the beam-to-beam contact.
The next step in the development will be the inclusion of Coulomb friction at the contact points. This will require definition and discretisation of some extra kinematic variables, analogously to what was done in Sects. 3 and 4 in the present frictionless formulation.