Frictional beam-to-beam multiple-point contact ﬁnite element

Thepresentpaperistheextensionofauthor’sear-lier research devoted to more accurate numerical modelling of beam-to-beam contact in the cases when beam axes form acute angles in the contact zone. In such situation with beam deformations taken into account, the contact cannot be considered as 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 from the point-wise formulation. The Coulomb friction model is introduced and advantage is taken from the analogy to plasticity. The penalty method is used to enforce the contact and friction constraints. The appropriate kinematic variables for tangential contact and their ﬁnite element approximation are derived. Basing on the weak form for frictional contact and its linearisation, the tangent stiffness matrix and the residual vector are derived. The enhanced element is tested using author’s computer programs and comparisons with the point-wise contact elements are made.


Introduction
The beam-to-beam contact represents a special case in the general field of 3D contact analysis. The efforts to analyze it numerically were started by Wriggers and Zavarise in [1] B Przemysław Litewka przemyslaw.litewka@gmail.com 1 Institute of Structural Engineering, Poznań University of Technology, ul. Piotrowo 5, 60-965 Poznan, Poland and subsequently continued in [2][3][4] where contact without and with Coulomb friction for beams of circular and rectangular cross-sections was considered. Since then further developments appeared and they included inclusion of thermal and electric coupling [5], smoothing procedures for 3D curves representing axes of beams in contact, e.g., [6] as well as a rigorous approach to the question of solution existence and uniqueness in the point-wise contact formulation, which was presented by Konyukhov and Schweizerhof in [7]. These authors focused their interest on the closest-point projection procedure, which in the particular case of the beam-to-beam contact leads to the orthogonality conditions, see [1]. The same authors used their approach to analyse the problem of rope wound around a cylinder and the question of knottightening [8]. The latter issue was also considered and solved by Durville in [9].
A fine analysis of various types of contact scenarios, including the beam-to-beam case, based on the geometry consideration was given in [10].
One important assumption in the majority of those approaches was the uniqueness of the mutual closest point projection procedure for two curves, ensured at least locally. The obvious consequence was the concept of the point-wise contact between beams. Yet, it is perfectly clear that such an approach fails in some special situations, e.g., when irregular assemblies of fibre-like objects are considered [11]. Thus, the more precise approach should cater for 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 [7,10].
Also the problem of contact between beams and other objects like rigid surfaces has been recently investigated. This topic was addressed in the further paper by Konyukhov and Schweizerhof [12] and by Gay Neto et al. [13,14]. Further- Fig. 1 Non point-wise contact between beams: a beams with axes crossing at acute angles, b conforming beams (b) (a) more, the latter authors took into consideration the problem of self-contact of a rope sustaining an off-shore riser [15].
The multiple-point beam-to-beam contact finite element was developed in [16] to model in a more precise way the situations, when contact between beams cannot be considered as point-wise. The frictionless interaction between beams was considered therein. This element is now enhanced to the frictional case. The main geometrical assumptions remain valid, in particular it is necessary that, locally, there still exists the unique solution to the closest-point projection equations. Hence, the presented approach is not applicable for the case of conforming/parallel beams which requires a completely different approach. In such a case a full 3D analysis with node-to-segment approach or a mortar method could be used instead of a special beam-to-beam treatment.
In the suggested model the line along which the contact may take place between two almost parallel or almost conforming beams, is discretised by three contact point pairs, see Fig. 1a. This model with two additional pairs is the most straightforward possible extension of the point-wise approach.
The adopted interpretation of the contact zone between two cylinders as a curved line is valid, if the assumption of small strains and undeformable cross-sections is kept, as was assumed in the previous papers [3,4,16].
Moreover, the contacting beams are modelled as onedimensional curves, with circumferential frictional effects neglected, only the tangential friction force components are taken into account. Hence, the proposed approach can be practically used in the cases when the beam cross-section dimension is much smaller than the beam length. The more general case was considered successfully in [7,10] using a geometrically exact covariant approach.
For a consistence of the present paper, in the Sect. 2 a brief introduction of geometrical considerations regarding the two additional contact pairs is presented repeating the ideas from [16]. To this end also some results of the finite element discretisation for the frictionless case are summarised in Appendix 3. In the remaining sections the attention is focused mainly on the frictional components of the contact formulation, since the detailed description the normal contact part can be found in the above mentioned paper. Furthermore, the present paper is mainly related to the influences due to the additional contact pairs, which are added to the formulation of point-wise contact given earlier in [4,17]. Section 3 includes the introduction of friction contributions to the weak form in the case of multiple-point contact. The frictional interactions are defined using the penalty method, and the analogy to non-associated plasticity is used with the distinction of stick and slip friction status. The consistent linearisation required for an effective use of the Newton-Raphson method to solve the non-linear problem at hand is also presented.
The kinematic variables present in the formulation of frictional contact are introduced and expressed in terms of displacements in Sect. 4. Again, the focus is put mainly on the terms due to the additional contact pairs and a comparison with point-wise contact variables is made.
Section 5 presents the finite element discretisation of frictional terms. The components of the weak form and its linearisation are suitably related to the nodal parameters of the beam elements. In this way the components due to friction of the residual vector and the consistent tangent stiffness matrix for the contact element are determined. Combined with the normal contact components from [16] they yield the complete frictional contact finite element formulation, which was embedded in the author's computer program.
Results of some numerical examples presented in Sect. 6 confirm the effectiveness and advantages of the multiplepoint beam-to-beam contact element with respect to the previously derived one-point contact element in the cases of beams forming acute angles. Section 7 contains final remarks, conclusions and outlook.

Additional contact pairs
Let us consider a contact between two beams labelled m and s. Their axes are shown in Fig. 2. The beams are of circular  3 Additional contact pairs of points on beam axes and additional normal gaps cross-sections with radii r m and r s . In the suggested multiplepoint beam-to-beam contact element three pairs of contact points are considered-the central contact points C mn − C sn and two additional ones. The central points are found using the orthogonality conditions, as described in [1,3] where the minimisation of the axes distance d N yields the local coordinates of the closest points on each beam −ξ mn and ξ sn (−1 ≤ ξ in ≤ 1, i = s or m). Hence, the normal contact gap function is To find the additional contact pairs, first, the contact candidate points C snb , C sn f are found on the beam s using the backward (subscript b) and forward (subscript f ) shift of local co-ordinates with respect to the value ξ sn for the central contact point as shown in Fig. 3. The definition of the shift ξ s of the local co-ordinate is taken as a result of simple geometric considerations regarding the layout of beam contact zones. Let us consider the plan view of the contact zone presented in Fig. 4. The crosssection A-A along the beam s is shown in Fig. 5. The distance between the centre contact point C sn on the axis of the beam s and the edge of the overlap region shown in the plan view (Fig. 4) is where ϕ is the angle between the beams axes in the plan view.
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 . In the numerical tests (see [18]) it was found that k = 3 yields the best results, so where l s is the length of the beam element within which the contact point lies.
Having located the additional points on the beam s, the orthogonality condition can be applied to find their counterparts C mnb and C mn f on the beam m where dot · denotes the scalar product. Generally, each of the conditions (6) and (7) represents a non-linear equation, which can be solved iteratively using the Newton-Raphson method. Following the lines of the solution for the central points with a set of orthogonality conditions (see [1]) with the use of the following expressions for the position vectors and their derivatives with respect to local co-ordinates, determined at the additional points x mnb = u mnb + x mnb,m ξ mnb x mnb,m = u mnb,m + x mnb,mm ξ mnb (8) x mn f = u mn f + x mn f,m ξ mn f x sn f = u sn f + x sn f,s ξ sn f x mn f,m = u mn f,m + x mn f,mm ξ mn f (9) 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 (1) and (12) the normal contact formulation for three separate contact point pairs can be formulated. The point to note is the character change of the beam-tobeam contact which is introduced by the presented procedure. Instead of the parity of both contacting beams like in the standard point-wise approach, the present model, due to the shifted additional points on the beam s, brings the distinction between slave and master bodies. It is advisable to treat the beam with a smaller cross-section as the beam s. With this provision, the resulting contact points are 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 (see Fig. 4).

Friction contribution to the weak form
The weak form for the frictional contact problem between two bodies can be written down in the following way: where the subscripts at the components of the virtual work δ denote the influences of the deformation of the body m, the deformation of the body s, normal contact N and friction T . Here the latter term only will be discussed. It is assumed  6 Components and resultant friction forces at the backward additional contact pair that the contact model is independent of the beam model itself, as was assumed in previous research (e.g., in [17]), while the normal contact component was analyzed in detail in [16]. For the case with three contact pairs as introduced in Sect. 2, the friction contribution to (14) can be given as where the subscripts c, b and f stand for the central pair, the backward pair and forward pair, respectively, as shown in Figs. 2 and 3. The contribution from the central point with the appropriate kinematic variables and the results of the finite element discretisation in the form of tangent stiffness matrix and the residual vector were presented in the papers on point-wise contact model [2,4,17]. The present paper is in principle devoted to the derivation of components concerning friction at the additional points, described as forward and backward contributions. Thus, the appropriate friction forces and virtual displacements are introduced which allow to write down: The friction forces F T mb , F T sb and F T m f , F T s f in each of the relations (16) and (17) do not represent action-reaction forces from Newton's third principle of dynamics. Rather than this, they have to be understood as pairs of components, which result from a relative movement of a contact point along each of the contacting beams, m or s. Let it also be indicated, that special care must be taken when computing the resultant friction forces F T b and F T f , when the angle between their components acting along the beams m and s is arbitrary. Their definition must include tangent vectors and proportionality coefficients, as introduced in the continuation. These resultant forces, can then be viewed as acting simultaneously on both beams, each of them can be treated as an action or a reaction force. This is depicted in Fig. 6 for the particular case of the backward contact pair.
Similarly, the corresponding variations of tangential displacements present in (16) and (17), i.e., the tangential gaps, must be related to the contacting beams, m and s. For each contact pair two independent movements are considered.
It must be also pointed out that the assumed way to model the beams as one-dimensional curves limits the possibilities of friction modelling to the tangential direction only. The circumferential friction force is not included in the analysis. From the practical point of view this means that the present model is limited to the cases when the cross-sectional diameter of contacting beams is much smaller than the beam length, like in the case of fibres assemblies, thin contacting cables or ropes, etc.
The analysis presented in this paper follows the footsteps of the point-wise contact formulation and the analogy between the frictional contact and the plasticity [19] is utilized. This concept includes the additive split of tangential gaps at backward and forward pairs into sticking (superscript e, for the elastic analog) and sliding (superscript p, for the plastic analog) components: The Coulomb law with the constant friction parameter μ is assumed to model the interface physical relation and the sliding functions f b and f f must fulfil the conditions where F N b anf F N f are the normal contact forces at the additional points (see [16]). An important thing to note is that the formulae (20) involve the resultant friction forces. There is a theoretical possibility to check the sliding criterion separately for each sliding component as was discussed e.g., in [17]. However, such an approach seems to be non-physical. Rather than this, each contact spot is to be treated consistently and thus it should have one contact status attributed, despite the fact that the resultant is determined using two components. Now the sliding rule can be introduceḋ which is a counterpart of the non-associated flow rule. The sliding rule can be integrated with respect to pseudotime using the incremental analysis. Hence, introducing the current new (subscript n) and previous (subscript p) configurations the following relations can be written down The second terms on the right hand sides of these expressions represent the values of increments of tangential gaps. The method to find these basic kinematic variables for friction is given in Sect. 4.
In the process of the incremental solution at the current new step, the trial values of sticking tangential gaps components are assumed using the current total gaps (22, 23) and the history variables, i.e., the previous values of the sliding gap components. The penalty method is now used to determine the trial values of the friction forces components where ε T is the frictional penalty parameter. The unit tangential vectors have to be introduced in order to evaluate the trial resultant friction forces. These vectors follow from the definition of the kinematic variables given in Sect. 4. Here, let it only be stated that they can be determined using the current and previous position vectors of the contact points. Hence, the tangential vectors can be found using The magnitude of the trial resultant friction forces can now be given by and the trial values of the sliding functions can be evaluated using With (32) in hand the status of the friction at the additional contact points can be checked. If the condition is met, then the backward pair is in the stick state, and if it is not-then sliding takes place. Similar check is carried out separately for the forward pair using In the case of stick the trial value of the friction force resultant is within the sliding limit and it corresponds to the correct current friction force components for the backward pair and for the forward pair. In this situation there is no incremental increase of the plastic part of the tangential gaps, either. For the backward pair one has and for the forward pair- If the sliding status is encountered, then the trial friction force exceeds the limit imposed by the sliding rule and the Euler return procedure has to be employed. For the Coulomb friction model it leads to the closed-form solution and the friction force components are given by their limiting values- for the backward pair and for the forward pair. It should be pointed out, that there are also proportionality parameters introduced in (39) and (40). Their presence is necessary due to the fact, that the sliding status check is carried out for the resultant friction force and the formulation of the weak form (16) and (17) requires the values of the friction force components along the beams m and s. These parameters are given by The sliding state involves the update of the plastic tangential gaps. For the backward pair one gets and for the forward pair- Let it again be emphasized, that the current state of the friction is checked separately for each of three contact pairs-central, backward and forward. Hence, the status at these points may be different and correspondingly to (33) or (34) the appropriate relations for the additional points (35)-(44) have to be used. The results for sticking or sliding cases can now be introduced into the weak form relations (16) and (17) to get the following results. The weak form for the sticking case for the backward points reads and for the forward points The linearisations yield In the case of the sliding status the weak form is for the backward pair and for the forward pair. These two formulae require additional parameters s mb , s sb , s m f , s s f to control the sliding direction. They can be found from the simple relations using the local co-ordinates of the current and previous contact points if these points are located in the same contact facet.
In the contrary situation, for instance the element numbers can be used. Besides, the values of the normal contact forces present in (39) and (40) are evaluated using the penalty method with the normal contact penalty parameter ε N introduced (see [10]) The linearisations of the weak forms (49) and (50) read The weak form expressions (45) (12) and thus concern the normal contact. They were given in detail in [16] and are not discussed in this paper. The variables related to the tangential gaps and to the proportionality parameters (41) and (42) are considered in the following section.

Kinematic variables for friction
The tangential gaps representing the key kinematic variables for friction, given by (22) and (23) Fig. 7. The increment of the tangential gap is found as the distance between the current contact points C m f n , C s f n defined by the current local co-ordinates ξ m f n , ξ s f n in the current configuration and the mappings C m f p , C s f p of the previous contact points C m f p , C s f p from the previous configuration onto the current configuration, defined by the previous local co-ordinates ξ m f p , ξ s f p . With the appropriate position vectors introduced, see Fig. 7, the gap increments can be given as For the backward pair one gets similar expressions These relations include the sliding direction control parameters given by (51) and (52). Beam m previous config. The local co-ordinates ξ sbn and ξ s f n of the current points at the beam s are found from (3) and those concerning the beam m follow from the orthogonality conditions (6) and (7) applied to the beam m only, as discussed in [16]. The variations and linearisations of these variables were also given in that paper. The only exceptions were δξ mbn and δξ m f n which were absent in the frictionless contact formulation. These variables can be derived from the orthogonality conditions. The corresponding tedious but elementary derivation is given in Appendix 1.
The local co-ordinates of the previous point mappings require a comment and a special treatment. The nature of the proposed contact model is that the additional points on the beam s are defined by the shift of local co-ordinate ξ s . This value depends, among others, on the angle ϕ included between the tangents to the beam axes at the central contact points C mn and C sn . If the angle decreases, than the spacing between the points is increased to cover the longer contact zone, see Fig. 8. The angle ϕ undergoes changes during the deformation process. Therefore, it can be concluded that the relative movement of the contact point on the beam s will result from two independent factors-the frictional sliding and the change of location due to the varying angle ϕ. If the beams are modelled as 1D curves and the circumferential effects are neglected, then the latter factor would result in an unrealistic behaviour and unphysical addition to the frictional force at the contact points. Therefore, the local co-ordinates of the previous contact point mapping cannot be treated as mere history variables and taken directly from the previous increment.
To overcome this problem it is suggested, that the relative movement of the additional points is directly related to the relative movement of the central points. The idea behind this assumption is that the three contact points modelling the continuous contact line for almost parallel beams are tied together. Hence, the following relations are proposed to define the local co-ordinates required to localize the previous contact point mappings C mbp , C sbp and C m f p , C s f p In other words it is assumed that the distance between the additional contact point and the central contact point on the beam s is determined by the same difference of the local co-ordinate ξ s in the case of the current points (C s f n , C sn ) and the previous point mappings (C s f p , C sn ). And similarly, for the beam m it is assumed that the distance between the current additional point C m f n and the mapping of the previous contact point C m f p is given by the same difference in the local co-ordinates (ξ mn − ξ mp ) as in the case of the central points (C mn , C mp ).
It should be pointed out, that simultaneously the influence of the angle ϕ on the formulation of the local co-ordinate shift in (5) remains unchanged. The results of this influence on the contact element tangent stiffness matrix and the residual vector were derived in [16] and are summarised in Appendix 3.
Now the tangential gap variations present in the weak forms and their linearisations derived in Sect. 3 can be evaluated. The results are similar to the kinematic variables for friction at the central points, which were given in [17] for the point-wise contact formulation. One gets for the forward points. The difference between the formulation for the central points and the additional points is the presence of the second components in the parentheses in (61) and (62) for the latter case. These components were null for the central points because the variation of the local co-ordinates for the previous point mapping was vanishingthe co-ordinates were treated as history variables and were independent on the current displacements variations. On the contrary, for the additional points these local co-ordinates are defined by (59) for the points on the beam s. The second components at the right-hand sides of (63) are the variables related to central points, which were expressed in terms of displacements in [17]. The first components at the right-hand sides of (63) concern the additional points; they are given in [16]. Finally, the right-hand sides of (64) can also be expressed in terms of displacements of the central points, as was shown in [16], too.
The linearisations of the weak forms variations include the linearisations of the tangential gaps and of the variations of the tangential gaps. The former variables are computed in the way analogous to the variations given by (61) for the points on the beam s. Likewise with (63) and (64), the particular components of these relations were expressed in terms of displacements in [16,17]. The weak form linearisation for the slip state includes also the linearisations of the proportionality parameters (41) and (42) related to the friction force components. The similar variables were derived in [17] for the central contact points. Their detailed form is given in Appendix 2 in (139) and (140).
In this way all the terms in weak forms and their linearisations are given by displacements and position vectors and can be subjected to the finite element discretisation.

Finite element discretisation
A finite element for contact uses the same nodal parameters as the finite elements for the beams themselves. In the presented examples the beams are modelled using the corotational beam finite element derived in [20].
An important issue influencing the form of the contact finite element is the smoothing technique used to ensure the proper continuity between the adjacent contact facets. The importance of smoothing in a general contact context was discussed for instance in [21]. The particularities of the smoothing process for the beam-to-beam contact with an analysis of several possible techniques were considered in [17]. Out of the models presented therein, the one using the inscribed curve and Hermite polynomials is used in this paper. This leads to the definition of the contact finite element involving two pairs of adjacent beam elements from two contacting beams, see Fig. 9. There are three nodes m1, m2, m3 of the beam m and three nodes s1, s2, s3 of the beam s involved in the formulation. The vector of nodal displacements related to the contact element can be given by where subscripts a = 1, 2, 3 and c = 1, 2, 3 in the displacement notation u sac and u mac refer to the node numbers and the co-ordinate numbers, respectively. However, the issue can become more complex, if one considers the fact, that the friction contact formulation for the three-pair layout involves indeed not one but six different points on each beam. These are: the current contact points C mbn , C mn , C m f n (C sbn , C sn , C s f n correspondingly for the beam s) and the mappings of the previous contact points C mbp , C mp , C m f p (C sbp , C sp , C s f p for the beam s). It must be assumed that these points may not lie within a single contact facet. The situation where the mapping of the additional (forward or backward) previous contact point is separated from the current point or where the additional contact point itself is separated from the central point should be taken into account. Therefore, the appropriate vectors of nodal displacements are introduced which may or may not coincide with the vector q, given by (67), related to the current central contact points C mn and C sn . It must be pointed out, that the history variables for friction and sliding are attributed to the particular contact pointscentral and additional. They are stored and revised at each iteration. If a local co-ordinate value at the current state exceeds the limit (from -1 to 1) then the corresponding number of contact segment and contact elements is changed accordingly.
The method of smoothing determines the form of the relation between displacements of an arbitrary point on the contact facet and the nodal displacements. Let us assume the following matrix representations for displacements at central points and at additional points where the subscripts i and I are related to m-beam or s-beam (i = s or m, I = S or M), the subscript j is related to backward or forward additional pair ( j = b or f ) and the subscript k is related to current (new) or previous configuration (k = n or p). Similar relations are introduced for the derivatives of displacements with respect to the local co-ordinates The explicit form of all the matrices G and H present in (72)-(75) was given in [17].
The variations and linearisations of the current local coordinates can be discretised to yield the following relations for the points on the beam s and with j = b or f , for the points on the beam m.
The matrices F and R present in (76)-(80) were detailed in [16] and for the self-consistence of the present formulation are repeated in Appendix 3. The missing discretisations in this group are for δξ mbn and δξ m f n and they yield The matrices involved in (81) are derived in Appendix 1. Similar relations can be found for the appropriate kinematic variables related to the previous contact points. Using (63)-(66) and the matrices introduced above one can get The discretisation of the proportionality parameters (139) and (140) leads to the following matrix relations The matrices present in (88) and (89) take the form with a = 0, 1, 2 and where the auxiliary matrices were introduced. For the simplification of notation some further matrices are used: For the notation in (96) and (97) the matrices given in (128)-(129) and the split of tangent vectors (28) and (29) into components were used. The finite element discretisation of the weak forms for the stick state, (45) and (46), leads to the following expressions where j = b or f and the following residual vectors were introduced R e T jma = ε T g e T m jn t m j S m ja R e T jsa = ε T g e T s jn t s j S s ja with a = 0, 1, 2. Similarly, the finite element discretisation of the weak forms for the slip state (49) and (50) yield where j = b or f and the residual vectors are with a = 0, 1, 2.
The linearisations of the stick state weak forms (47) and (48) after the discretisation can be given in the form where i = s or m, j = b or f, a = 0, 1, 2 and c = 0, 1, 2.
The contact element residual vectors and tangent stiffness matrices (99, 103) or (101, 110) are ready to be embedded in any computer program of beam-to-beam contact analysis. An inherent contact search routine, like the one described in [17], determines for each increment and each iteration, used within the Newton-Raphson iterative solution scheme, if a contact contribution has to be switched on-separately for the central, forward and backward contact candidate pair. In a case of an active contact pair the appropriate components resulting from the normal contact (detailed in [16]) as well as the friction components derived in this paper have to be added to the global residual vector or stiffness matrix. At each contact pair the friction status must be checked using the criterion (33)-(34) and the appropriate, stick or slip versions have to be applied.

Numerical examples
This section of the paper includes two examples of frictional contact analysis for beams. The examples are meant to verify the performance of the newly developed three-point contact element in the frictional regime and compare it to the standard point-wise approach. The axes of beams in the contact facets are approximated by Hermite polynomials as inscribed curve segments, as discussed in [6,17]. The contact search algorithm presented in [17] is applied. In both analysed cases beams axes form acute angles in the plan view, i.e., the layout for which the multiple-point contact formulation is devised. In the both examples data and results are given without any specified units but they may be interpreted using any consistent unit set.
The computer times for both considered approaches differ. However, it is worth to emphasize, that the introduction of two additional points, meaning the 200 % increase in the number of contact candidates or active constraints, increases the time of computations only by about 20-30 %, depending on the beam layout. This ratio is quite favourable because the computations for the additional points involve several The values of the penalty parameters were adopted to keep the penetration below 4 % of the beam diameter.

Example 1
Let us consider a set of two cantilever beams 1 and 2 with the lengths l 1 = l 2 = 6.021 and the angle between their axes ϕ = 9.46 • in the plan view, as presented in Fig. 10. 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 Coloumb friction parameter is taken as μ = 0.1. The beams are forced into contact by the nodal displacements 1 = 2 = 0.2 applied in 60 increments at the free ends of both beams. Each beam is discretised by 5 co-rotational finite elements proposed by Crisfield [20]. The penalty parameters used in the calculations are ε N = 3 × 10 3 ε T = 1 × 10 3 .
The friction state in this example is slip. The values of the measure of the relative sliding-the plastic tangential gap g p T during the deformation process in 60 increments are presented in Figs. 11a, b for various contact points and various formulations. Due to symmetry between the beams the sliding distance is the same along each of the beams, so only one set of results for one beam is presented. Figure 11b gives a comparison between the point-wise formulation and the present multiple-point formulation. It is clear that the new approach is characterised by larger easy for sliding. It is due to the fact that the normal contact forces are distributed along the contact zone (line) instead of being concentrated at a single spot. In this situation the resistance to sliding is smaller. Figure 11a showing the accumulated sliding distance for the central and the backward contact points indicates, that the three contact points (the results for the forward and the backward points are virtually identical) move together. The difference between the distance is due to the fact that for technical reasons the computations for the additional points must start one increment later than for the central point. This difference is kept constant during the process, what means that the points indeed move together along the beam. Figure 12 presents the results of the friction force component along one beam; again due to symmetry the results are identical for beams 1 and 2. In the slip state the friction force for the Coulomb model is proportional to the normal contact force. The distribution of the contact interaction along the contact zone in the multiple-point formulation results in smaller values of penetration and the contact forces, what is evident in the presented graph.
The final deformed shape of beam axes is shown in Fig. 13.

Example 2
A symmetric assembly consisting of four identical beams is increments at all the beams tips. Each beam is discretised by 10 co-rotational finite elements developed by Crisfield [20]. The problem was solved using the penalty parameters ε N = 10 3 and ε T = 10 2 . Due to symmetry there are only two different contact points, depicted by A and B in Fig. 14. Also, at every contact spot the contact behaviour including the sliding distance, elastic tangential gaps and friction force components is the same for each of two contacting beams.
The example presented here is characterised by the stick state. Figure 15a, b presents the evolution of the elastic part  Fig. 16a, b presents the similar comparison, but related to the plastic part of the tangential gap g p T which in this case represents the total sliding distance in the resistanceless case.
In the fricional case it is clearly visible that the distribution of contact interactions in the multiple-point formulation leads to the smaller values of elastic tangential gap and the resulting friction force than in the case of typical pointwise model where the interaction is concentrated in a single point.
It can also be noted that, similarly as in Example 1, the three points in the presented contact model are moving together. The evolution of the tangential gaps in the additional points, forward and backward is virtually identical, that is why only one of them is depicted in Figs. 15 and 16. Each of these two evolutions follows the one for the central point, with a shift related to the starting value of the gap at the central point computed at the very first increment, when due to computational reasons the additional points cannot be included in the analysis.
The final deformed shape of beam axes is shown in Fig. 17.

Concluding remarks
The present paper includes an extension of the multiple-point beam-to-beam contact formulation [16] to the frictional case. The introduction of two additional contact pairs allows to model the situation when the contacting beams form acute angles and the resulting contact zone is not a single point. One could certainly think of introducing more points but the presented approach is devoted only to the simplest possible three-point approach. A further analysis of models with more additional points constitutes an interesting topic for a future research.
The paper includes the full formulation of the contact contributions to the weak form and its linearisation describing the general contact problem. The appropriate kinematic variables are introduced taking into account the possibility of independent frictional movement of the contact point along each of contacting beams. These variables are determined for the additional points using the central point coming from the point-wise formulation as a reference. The required variations of the variables with respect to the displacements are computed and then discretised using the finite element methodology with appropriate contact finite element introduced. The resulting consistent tangent stiffness matrix and the residual vectors are ready to be embedded in the finite element model of contacting beams. The author's computer programs written in Fortran language were used to solve numerical examples which confirm the effectiveness of the proposed solution in the frictional cases.
It should be emphasized once again, that the proposed beam-to-beam contact element with additional points is based on the existence and uniqueness of the closest point location procedure. If this procedure is unable to determine uniquely the closest points which are the central points in the suggested model, than the additional points cannot be found, either. Hence, the proposed approach may be viewed 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 another. It might also be used as an intermediate link between the point-wise contact approach for large angles between tangents and the node-to-segment or mortar approach for conforming beams.
It is planned to use the derived formulation in the general beam-to-beam contact model which would automatically switch between the point-wise contact for perpendicular or almost perpendicular beams and the full node-to-segment model for conforming beams, while the presented multiplepoint contact would be a transition link between the two classical contact models.
With the simplifying notation +x m f n,m x s f n,s (R s + R s ) (136) to yield the discretised form of the kinematic variables δξ m f n . The variable δξ mbn can be derived in the same way with all the subscripts f replaced by b.

Appendix 2: Details of linearisation of friction kinematic variables
This appendix includes some relations used in the linearisation of kinematic variables for friction.
where the subscript i is related to s-beam or m-beam (i = s or m) and the subscript j is related to backward or forward additional pair ( j = b or f ). Linearisations of proportionality parameters given by (41) and (42) introduced in friction contributions to the weak form (53) and (54)