New benchmark problems for verification of the curve‐to‐surface contact algorithm based on the generalized Euler–Eytelwein problem

Development of the numerical contact algorithms for finite element method usually concerns convergence, mesh dependency, etc. Verification of the numerical contact algorithm usually includes only a few cases due to a limited number of available analytic solutions (e.g., the Hertz solution for cylindrical surfaces). The solution of the generalized Euler–Eytelwein, or the belt friction problem is a stand alone task, recently formulated for a rope laying in sliding equilibrium on an arbitrary surface, opens up to a new set of benchmark problems for the verification of rope/beam to surface/solid contact algorithms. Not only a pulling forces ratio TT0 , but also the position of a curve on a arbitrary rigid surface withstanding the motion in dragging direction should be verified. Particular situations possessing a closed form solution for ropes and rigid surfaces are analyzed. The verification study is performed employing the specially developed Solid‐Beam finite element with both linear and C1 ‐continuous approximations together with the Curve‐to‐Solid Beam (CTSB) contact algorithm and exemplary employing commercial finite element software. A crucial problem of "contact locking" in contact elements showing stiff behavior despite the good convergence is identified. This problem is resolved within the developed CTSB contact element.


INTRODUCTION
• a rigid punch problem, configured as pressing of a rigid block into elastic half space, see References 1,3. . The next set of models is very popular in tribology for both verification and validation. They are known as generalization of the Hertz problem into adhesion region: • the Johnson-Kendall-Roberts (JKR) model; 4 • the Derjaguin-Muller-Toporov (DMT) model. 5 The Hertz problem, however, including friction becomes even more specific especially for rolling applications, see References 1,6. These solutions are widely used for verification and corresponding to two-or three-dimensional numerical models resp. to 2D or 3D finite elements.
A standalone problem for objects of different nature is the frictional contact between a rope and a cylinder initially formulated by Leonhard Euler as the 2D problem and published in "Remarks on the effect of friction on equilibrium". 7 This problem after popularization by Eytelwein in his text book 8 nowadays is known as the Euler-Eytelwein problem or the belt friction problem as well as the Euler formula for rope friction. This problem, even formulated as 2D, is widely used in wire rope mechanics. 9 In Reference 10 this problem has been generalized into the sliding contact of a rope on a orthotropic surface with arbitrary geometry.
Among computational contact problems a problem dealing with rope/beam to shell/solid is requiring special algorithms, which is related to beam-to-beam contact algorithm, see for example selected publications. [11][12][13][14] The current contribution is aimed to formulate a set of particular problems based on the generalized Euler-Eytelwein formula for the verification of computational contact algorithms for rope/beam to shell, solid type. This set is formulated as frictional contact between a rigid surface of various geometry and a rope with particular geometry laying on this surface in equilibrium.

GENERALIZED EULER-EYTELWEIN PROBLEM: PARTICULAR SOLUTIONS
The famous Euler-Eytelwein problem describes the equilibrium condition of an unextensible rope being pulled on both ends and laying on a rigid cylinder represented by a circle in 2D plane. This problem gives a formula which describes the ratio of pulling forces T and T 0 on both ends of the rope for a given arch contact angle and friction coefficient between the rope and cylinder. The rope starts to slide if the ratio reaches the following critical value: For a given friction coefficient and if the rope is wrapped several times around the cylinder a huge difference between the pulling forces on the ends of a rope can be balanced. This problem has been known in practice since the 18th century when 1786 Euler has published the famous solution (1) in Reference 7. However, if a rope or a belt is wrapped more than once around the cylinder then the rope possesses no longer the circular in plane configuration due to its thickness and three dimensional description becomes necessary. In Konyukhov 10 this problem has been generalized for a rope laying on an orthotropic rough arbitrary surface. We here represent the necessary results in order to concentrate further on particular surfaces. The understanding of differential geometry for both a free curve, and for a curve laying on a surface are essential in order to describe the equilibrium condition. The Serret-Frenet frame for a free curve r( ) is constructed from orthogonal vectors: • a unit tangent vector evaluated as a derivative of the curve vector r with respect to its arc length s: = dr ds , where ds = √ r ′ ⋅ r ′ d ; • a unit normal vector evaluated as a normalized derivative of with respect to the arc length s: = . This vector is pointing into the center of curvature; • a unit binormal vector evaluated as the cross product = × . F I G U R E 1 Gaussian frame 1 , 2 , n for a surface and Darboux frame , n, g for a curve on the surface Equation (2) represents the three Serret-Frenet formulas, which describe the relation between the basis vectors and their derivatives with respect to the arc-length parameter s: where is the torsion and k is the curvature of the curve.
In case of a curve laying on a surface, the geometrical characteristics are described using the Darboux frame , n, g. In order to evaluate these vectors the Gauss surface frame ( 1 , 2 , n) for a surface ( 1 , 2 ) first should be defined: The Darboux frame is formed with the unit tangent vector to a curve , the unit normal vector to a surface n, and the unit tangent normal vector g = × n.
The Gauss, the Serret-Frenet and the Darboux frames are depicted with their geometrical relationships in Figure 1. It is easy to see, that the the Serret-Frenet frame is obtained via rotation of the Darboux basis over the vector at the angle between a unit normal vector of the surface n and a unit normal vector of the curve . This relation is written as the following transformation, Equation (4): The Darboux frame defines three essential directions for kinematics of a rope on a surface: • vector defines the pulling direction, in which the tangential pulling force T for a rope is acting; • vector g defines the dragging direction, in which the dragging force G for a rope is acting; • vector n defines the normal direction, in which the normal force N for a rope is acting;.
The curvature vector k of a free curve can be decomposed into two components in the Darboux frame: where k n is the normal curvature and k g is the geodesic curvature.
Remark 1 (Geodesic curvature for a geodesics). A geodesic line, defined as a curve representing the shortest part between two points on a surface, has, at the each point of this curve, zero geodesic curvature k g = 0.
The generalization of the famous 2D Euler-Eytelwein problem for arbitrary surfaces with orthotropic properties given in the Darboux frame is derived in Reference 10 as the following theorem: Theorem 1. If a rope is laying in equilibrium under tangential forces on a rough orthotropic arbitrary surface then three following conditions (all of them) are satisfied: For further verification we establish a set of benchmark problems with particular geometry of a surface and a curve (representing mid-line of a rope) laying on it: 1. A helix (representing geodesics) on a cylinder; 2. A circle on a cone; 3. A geodesic line on a cone; 4. A circle on a sphere; 5. A rope in a tubular sinusoidal profile; 6. A rope on an ellipse as a plane curve.
In these particular cases, except 3 the integral in Equation (10) is taken in a closed form. The case 3, however, allows a remarkable closed form solution for a geodesics and only numerical integration of the integral equation (11) is required. Cases 1-4 will be used further for the verification of a particular finite element model.

Specifications of the critical case for the finite element verification
Various finite element models will be employed further for the verification: • combinations of shell finite elements for a surface and beam finite elements for a rope together with the beam-to-surface contact type; • combinations of Solid-Beam finite elements with various approximation for both a surface and a rope together with curve-to-solid beam (CTSB) contact type.
A surface, on which a rope is laying, in both models is assumed to be rigid and is modeled with fixed degrees of freedoms for corresponding finite elements.
In order to represent the critical state -the whole rope is sliding condition for verification of Equation (10)-a rope is attached at both ends to spring finite elements, see for example Figure 2 for the cylinder. Displacements are applied at the selected end of a rope until all contact elements will be detected in sliding condition while another spring node is fixed.
We assume isotropic friction in all further considered cases: g = n = .
F I G U R E 2 Specification of a finite element model. A rope as a helix over a cylinder

A rope as a helix on a cylinder
The case of a rope wrapping a cylinder is the most common case for the application of the 2D Euler-Eytelwein formula representing circular geometry, however, in reality a rope is laying along 3D helix curve on a cylinder, see Figure 3. A cylinder is given by the following vector equation: where h is height of the cylinder. Vector equation of the helix is given as: where H is a pitch of the helix. All necessary for further computation geometrical parameters of a helix are: • the arc length • the curvature k F I G U R E 3 Representation of a helix on a cylinder Remark 3 (Curve with constant curvature and torsion.). It should be noticed that the curvature k and torsion of a helix are both constant. Inverse statement is also valid: if the curvature k and torsion are both constant and not zero for a curve then this curve is a helix. A curve with zero curvature k = 0 and zero torsion = 0 is a straight line, and a curve with nonzero constant curvature and zero = 0 torsion is a circle.
Computing a normal to a cylindrical surface in Equation (12) as n = (cos , sin , 0) and a normal of a helix in Equation (13) as = (− cos , − sin , 0), we obtain for the angle its cosine cos = (n ⋅ ) = −1. This leads to the zero geodesic curvature k g = k sin = 0, see Equation (6). Thus, with regards to Remark 1, a helix is a geodesic line of a cylinder.
The critical force ratio for the helix is obtained from Equation (11) for a geodesics: Taking into account equations for curvature (15) and for the arc length (14), we obtain: Remark 4 (Maximal pulling forces at both ends are leading to geodesics). If in an imaginary experiment we are trying to balance the maximal forces at both ends on of a rope (i.e., both forces T 0 → ∞, T → ∞), laying in a certain domain  on a rough arbitrary surface outside the equilibrium domain in the dragging direction assuming that the condition opposite to Equation (8), namely − g ≥ tan ≥ g is satisfied, thus the rope is allowed to slide in the dragging direction in the domain , then relative motion in the dragging direction will be eliminated and the curve in the limit will follow the geodesics in this domain  on a surface, see proof in Reference 10 and representation in Figure 4.

A rope as a circle on a cone
Our next particular example for the verification is a rope following a conic section as intersection of a right circular cone with a plane. The cone surface is given by a generatrix which makes an angle * with the vertical axis and can be written either using and z as coordinates: where h 1 and h 2 are z-coordinates of the lower and upper base surfaces of the cone segment with the H = h 2 − h 1 as height of the segment; or using and r as coordinates: where R 1 and R 2 are radii of the lower and upper base surfaces of the cone segment, respectively. For the numerical example we consider a cone segment with the height H and angle * produced by intersection of a cone with planes z = h 1 and z = h 2 ), respectively. A half of a circle with radius R follows the intersection of the cone segment with a plane z = h c ∈ (h 1 , h 2 ), see Figure 5.
According to the second condition in the theorem in Equation (8), an equilibrium is possible if the friction coefficient in dragging direction g does not increase the critical value of friction crit = | tan | depending in the current case on the cone angle * , namely = − * . All necessary equations for the circle geometrical parameters are evaluated as: • the arc-length • the curvature • the normal curvature F I G U R E 5 Specification of a finite element model. A rope with a half-circle configuration over a cone • the geodesic curvature It should be noted that the equilibrium in this case is reached for a rope which is not forming a geodesic line. This is the simplest benchmark test for both the second and third conditions of the theorem, namely: • equilibrium is possible if > tan * , sliding in the dragging direction starts at crit = tan * ; • the critical ratio of pulling forces in Equation (10) is derived, taking into account all necessary geometrical parameters in Equations (20)-(23), as and finally the compact formula is obtained:

A rope as a geodesic line on a cone
A geodesic line can be defined as the shortest distance curve on a surface. Then the main property is the zero geodesic curvature on the surface, see Remark 1. Geodesics on a cone, see Figure 6, does not follow the intersection of a cone by a plane in general. A remarkable property of the geodesics on a cone: if the cone would be unrolled into a plane (hemidisk), then the geodesic line is represented by a straight line on this plane. The problem of finding a geodesic line on a cone has a closed form solution, see Reference 15. This solution is written as a curve on the cone with regards to Equation (18). The equation for geodesic line between points with coordinates ( 1 , z 1 ) and ( 2 , z 2 ) has the following form: with z( ) = z 1 cos(a ) − sin(a ) ⋅ z 2 cos(a 2 ) − z 1 z 2 sin(a 2 ) and a = sin( * ).
We consider a geodesic line laying on the outer cone surface, therefore, the angle between normal vectors = and curvatures are computed as F I G U R E 6 A rope as a geodesic line on a cone. Specification of a finite element model Further geometrical parameters, such as the arc-length s and the full curvature k should be computed by their definition for a spatial curve: • the differential of an arc-length • the curvature The ratio of pulling forces for the critical state is obtained from the Equation (11) for geodesics: Remark 5. Computation of the critical ratio of forces for the geodesic line on a cone in Equation (31) cannot be fulfilled in the closed form. For further verification an appropriate numerical method for computation of the integral with a given tolerance should be employed.

A rope as a circle on a sphere
The case of a circular rope on a sphere is kinematically more reachable in comparison to that in previous cases, which has included only developed surfaces with zero Gaussian curvature-cylinder and cone, because a sphere is a double curved surface. Let r be the radius of a circle, R be the radius of a sphere, and H the distance between the principal circle of the sphere and a plane giving this circle as intersection, see Figure 7. From the right triangle, see Figure 7, we can directly derive trigonometrical parameters for the angle between a circle and a sphere normal vectors n: Circle as intersection of a sphere with a plane-geometrical parameters The equilibrium conditions from the second and the third conditions of the Theorem 1 are written as follows: • equilibrium is possible if > tan * , the sliding starts in the dragging direction at • for the case of isotropic friction = g = n , the critical state of the Euler-Eytelwein formula with regards to Equation (10) is derived as and taking into account Equations (32) and (33) These conditions are forming an effective benchmark test for verification of the rope/beam to surface contact algorithm because only in the region near the great circle defined by Equation (34) equilibrium is possible, otherwise the rope/beam starts to slide in the dragging direction. This will be observable, however, only for the transient numerical solution. For the quasistatic numerical solution a priori disconvergence will be observed due to disability to find an equilibrium state.

Application to a plane curve with an inflection point. Plane sinus profile: a rope in a tubular profile
The curvature k for the plane 2D curve allows simple geometrical interpenetration as a derivative of the rotation angle for the unit tangent vector with respect to the arc length s, see Figure 8: The 2D case can be interpreted as a 3D case, in which a cylindrical surface with a generatrix parallel to the third axis OZ is generated by a curve in XOY -plane. All active forces are parallel to the plane XOY , respectively. Theorem 1 for arbitrary 3D case was especially represented for 2D case, see Reference 16, as follows: If a rope is laying in equilibrium under tangential pulling forces on a rough curve in 2D, see Figure 8; then the following conditions (all of them) are satisfied:

No separation -normal reaction N is positive for all points of the curve:
this condition can be interpreted as • curvature k is keeping its sign; or • angle is a monotone function.

Limit values of the tangential forces:
The forces at both ends T and T 0 are satisfying the following inequality where the is coefficient of friction between a rope and a curve in the tangential direction (the dragging coefficient of friction g ≠ 0 just should be zero); and 0 − is the final angle of rotation for the unit tangent , see Figure 8.
It is remarkable that the form of famous Euler-Eytelwein formula T T 0 = e ( 0 − ) is completely preserved. As an interesting application see comparison of the results for a circle and an ellipse in Appendix (presented in book 16 ). It is possible to apply the theorem to a nonmonotone (s) curve assuming that "no separation" condition (7) is fulfilled by the construction. This is an example of a tubular profile and a rope inside it: in this case "no separation" condition if fulfilled.
Consider a sinusoidal tubular profile with amplitude A and periodical on the length l as shown in Figure 9. This profile is described in its longitudinal coordinate x by the function: Nonmonotone angle curve has inflection points at x n = n. A rope laying inside the sinusoidal tubular profile The function has inflection points at periodical coordinates x n = n, n ∈ N. Tubular profile provides the double-side contact constructively because in the region x ∈ contact appears with a lower surface of the tubular profile, or with an upper part of the sinusoidal curve, and in the next region x ∈ contact appears with an upper surface of the tubular profile, or with lower part of the sinusoidal curve. Taking derivative of the function (39) and evaluating it for x = 0 we obtain tan 1 , see Figure 9: The slope for tangent line gives an initial angle 1 : The finite angle of rotation − 0 for the unit tangent should be taken for the region in which the corresponding function is monotone (monotonicity domain), that is, in the region between inflection points, or for the region with the same sign of curvature. This can be easily computed as an increment of the slope for tangent line to the graph 2 − 1 , see Figure 9. The critical ratio of tangential forces is computed in each part of monotonicity domain via increments of the slope Δ : So far we are free to split the function in various ways as soon as it is geometrically easier to capture the monotonicity domains. Below we reproduce solutions for two and four monotonicity domains. a) Solution for 1 4 of the period. In this case it is very easy to see from the graph in Figure 9, that the increment of the slope is Δ = 1 − 0 = 1 .

OVERVIEW OF THE CTSB CONTACT ALGORITHM
A special CTSB contact algorithm was developed in References 17,18 based on the Solid-Beam finite element with elliptical cross section, see Reference 19. In order to construct a basic solid-beam finite element with linear approximation of the mid-line, six nodes are used, as can be seen in Figure 10. Only three nodes are sufficient to define an ellipse with principal axes 1 − 3 and 1 − 4 in the local polar coordinate system. Nodes 1 and 2 define the mid-line of the element and nodes 1, 3, 4 define the cross-sectional ellipse on the left hand side, and nodes 2, 5, 6 the one on the right-hand side, respectively. The CTSB contact algorithm consists of the master side, which is a deformable surface taken from the one solid-beam element and the slave side, which consists of the mid-line from another solid-beam element. Deformation of finite elements on the slave side in the transverse direction is not taken into account for the contact, because for the computation of penetration the slave surface is assumed to be a tube with radius R. From a geometrical point of view, the master side for the linear element is a surface spanned by vectors between two arbitrary elliptical bases. Thus, it includes the possibility to represent the following surfaces: linear and elliptic cylinders, right and elliptic conical surfaces away from the apex. The linear solid beam finite element is easily extended into modeling of C1-continuous iso-geometric finite element describing tubular surfaces with varying elliptical cross-sections. These surfaces are classified as sweep surfaces obtained by sweeping of the cross-section along a curve, see more about surfaces in Reference 20. The slave side is represented as a surface of the cylinder with circular base for linear elements, and a tubular surface enveloping the slave curve with constant radius for C1-continuous set of elements, respectively. In order to describe the contact kinematics of CTSB element two coordinate systems are employed: the Gaussian surface frame 1 , 2 , n for external surface of a master solid-beam element and the Serret-Frenet frame , , for a mid-line of a slave solid-beam element, see Figure 11.
In the case of linear approximations for CTSB contact element the following shape functions are involved:

F I G U R E 11
Curve-to-solid beam contact algorithm. Gaussian and Serret-Frenet frames • six shape functions for the master surface arising from the shape functions for the corresponding solid-beam finite element with radial coordinate set up to r = 1: where coordinate ∈ [0, 2 ] corresponds to a circumferential direction in a local polar coordinate system and M i are standard Lagrange linear approximation shape functions for the mid-line: where ∈ [−1, 1] is an axial coordinate; • two linear shape functions for the mid-line of the slave element: It should be noticed that "shear locking free" options for bending should be taken into account for better performance of solid-beam elements in full similarity to techniques widely used for shell and solid-shell finite elements. The solid-beam finite element used in the current contribution is the "shear locking free" element.
For the unified description in a covariant form we reintroduce variables on the master surface 1 ≡ and 2 ≡ . Following the geometrically exact theory of contact, see in Reference 21, it is sufficient for the forthcoming finite element discretization of all necessary contact parameters to introduce approximations of the relative displacement vector where the approximation matrix A( 1 , 2 , ) is defined as: (51) and x is a nodal vector of the CTSB contact element defined as direct sum of master and slave nodal coordinates: The generalization into arbitrary C1-continuous approximation is straightforward: • arbitrary C1-continuous approximations functions M i ( ), i = 1, 2, … , k are chosen employing, for example, the corresponding set of NURB, iso-geometric etc. functions with k parameters as a number of nodes or knots instead of two in Equation (48). This approximation is applied for both mid-lines of master and slave solid-beam elements. The derived solid-beam finite element has 3k nodes/knots.
• the set M i ( ), i = 1, 2, … , k constitutes directly the shape functions for a slave part in the matrix A as a block with dimension 3k , see Equation (51) • the set of shape functions for the C1-continuous finite element is derived similar to the linear shape functions set in Equation (47)  The weak form for the CTSB contact algorithm is evaluated along the slave curve with its differential dl s = |r ′ |dl s as: Within the computational algorithm, this contact integral is computed via a numerical quadrature formula of Gauss or Lobatto types. Thus, the contact points are the integration points j corresponding to the slave line r( j ). They are projected onto the master surface employing the Closest-Point-Projection (CPP) procedure defining the shortest distance 3 between a contact slave point r(eta j ) and a master surface ( 1 , 2 ), see Figure 11. Convective coordinates of this projection point on the master surface are 1 j,c , 2 j,c . The measure of normal contact interactions is a penetration and is computed by subtracting the radius R of a slave tube from the shortest distance 3 at each contact point r( j ): The measure of tangential interactions, following the geometrically exact theory in Reference 21, first, is formulated in the rate form aṡi where (v m − v s ) is a vector of relative velocity corresponding to the vector of relative displacement in Equation (50), a ij are contravariant metrics tensor components and i are contravariant surface vectors conjugated with covariant ones i . For the numerical implementation an incremental measure Δ i is then computed. The residual vector for nonfrictional contact is computed by substitution a relative virtual displacement vector (50) into a normal part of the contact integral (53): Integral in Equation (56) in the current contribution is computed via the complex numerical integration formula with subdomains: where • m -number of subdomains in which the slave element l s is split; • N ip -number of integration points for the chosen quadrature formula (Gauss or Lobatto) in each subdomains; • J j = |r ′ ( j )| -Jacobian for the slave length element; W j are weights of the quadrature formula.
The normal force N is computed with regards to the penalty method at each integration points as where N is a normal penalty parameter. Thus, only penetrating points with p < 0 are contributing into the integral in Equation (56). The residual vector for a tangential part is computed in a similar fashion, however, taking into account conditions in Equation (58), namely, only penetrating points with p < are contributing into the residual vector R T : The real contravariant components of the tangential force T i,real are computed via the standard return-mapping scheme for the Coulomb friction law at the current load step (n): where absolute value of the tangential force T is computed via contravariant components as follows: The trial contravariant components of the tangential forces T i,tr (n) are computed via the following finite difference scheme (incremental analog of the rate Equation 55): where T is the tangential penalty parameter and i 0 , i = 1, 2 are initial values of the convective coordinates from the CPP procedure, assuming sticking and zero initial tangential force at the initial state of contact interactions.
Remark 6. The updated scheme in Equation (62) is correct for all cases of unidirectional loading (i.e., without the reverse loading). For more general case the initial sliding point with coordinates i 0 should be updated within the return-mapping scheme for the "sliding" case also, see more details in Reference 22. This procedure, however, for the forthcoming examples, including only unidirectional loading, is not necessary.
The stiffness matrix K in the case of contact with friction consists of normal K N and tangential parts K T , which are derived from the discretization and corresponding transformations of the contact integral, respectively, from the normal and tangential parts of the contact integral, see more about their derivations in References 22 and 23: The normal part is responsible for geometrical nonfrictional contact where h ij are contravariant components of the curvature tensor for the master surface and A ,j = A j , j = 1, 2 are partial derivatives of the approximation matrix A.
The tangential part K T depends on the real tangential force T real and is computed for the sticking case as and for the sliding case with regards to Equation (60) as In all parts of the tangential part for both sticking and sliding cases T is a trial force with corresponding covariant T i components. The structure of the introduced stiffness matrices is explained in References 21,23 in more details.
The contact is implemented as a separated finite element. These are following terms: the residual vector for normal forces R N in Equation (57) with the corresponding tangent matrix K N in Equation (64) and the residual vector for tangent forces R T in Equation (59) with the corresponding tangent matrix K T in Equation (65) for sticking or in Equation (66) for sliding depending on the result of the return-mapping scheme in Equation (60). Further study on variation of numerical integration is concerning both the number of subdomains m and the number of integration points for the chosen quadrature formula N ip exclusively for the contact finite element.
Another interesting approach for modeling a beam with elliptical cross-section together with the contact has been published in References 24 and 25. The Simo-Reissner geometrically exact beam theory has been employed for the beam finite element and contact has been detected directly as the shortest distance between two elliptical cross-sections.

VERIFICATION OF SPECIAL CASES
In this section particular cases with the analytical closed form solutions, shown in Section 2, are used here for the verification, namely, for comparisons with numerical solutions and further discussions. Numerical solutions are performed with the following finite element approaches: • Computation in commercial software package ANSYS; • Computation via the specially developed CTSB contact algorithm in the academical university software package FEAP.
In order to ensure the quality of the results the standard mesh convergence procedure on sequences of different meshes for each new geometry has been employed. Visual match of two resulting curves computed with different meshes for the corresponding diagram derived in the range of friction coefficients ∈ [0, 0.3] was the decision criteria to select the final mesh. In ANSYS computation it was found out that starting from three beam finite elements being in contact with a single shell finite element is the sufficient discretization. All following figures are showing final meshes in ANSYS model. Each FEAP model includes a rigid body as a single solid-beam finite element representing the exact geometry (cylinder and cone) and corresponding meshes of solid-beam finite elements for a rope.

Computation in commercial software package ANSYS
The rope is modeled using the three-node BEAM189 finite elements with quadratic approximation along the mid-curve and six degree of freedoms per node. Element BEAM189 according to ANSYS Manual has been implemented with regard to the Timoshenko beam theory. These elements are providing approximations of the curvature, which is necessary for the analyzed problems. The rope is discretized with 100 BEAM189 finite elements with radius r rope = 0.001 m. A large longitudinal together with low shear stiffness's are taken in order to simulate unextensible rope kinematic behavior via beam finite elements. In the follow-up simulations the elasticity modulus of taken is E = 2.1 ⋅ 10 11 N m 2 and the shear modulus is taken 4000 times smaller providing shear soft rope behavior.
All rigid surfaces-cylinder, cone and sphere-are modeled using the SHELL281 shell finite elements with quadratic serendipity approximations. The elasticity modulus of the shells is E = 2.1 ⋅ 10 11 N m 2 and all degrees of freedom for all nodes are fixed, which provides the rigid body condition of the contact surfaces for the Euler-Eytelwein problem.
Though the elasticity modulus does not influence the behavior of surfaces with all fixed nodes, this value is required by the conventional ANSYS contact algorithm, namely for the computation of penalty parameters N , T for the contact elements.
Only isotropic friction is involved in forthcoming computations, therefore, we will not distinguish between friction coefficients in "pulling" and in "dragging" directions: = g = . The rope is supported on a soft linear spring at one end and the tangential force at another end is applied, as discussed in Section 2.1. The quality of numerical results is estimated by the comparison of a force ratio at both ends with the analytic solution. The pulling force or corresponding displacement in numerical simulations is applied until the "sliding" status is reached at all contacting elements.

Computation employing the CTSB contact algorithm in FEAP
The rigid surfaces such as cylinder and cone are modeled geometrically with a single solid-beam element and the rope is modeled using N el Solid-Beam elements. Rope nodes/knots are positioned such that the zero initial penetration relatively the rigid surface is satisfied. The specially developed CTSB contact algorithm is used. The model is implemented into the academic finite element code FEAP-MEKA. It is important for rope solid-beam finite elements, that the axes for the cross-section 1 − 3 and 1 − 4, see Figure 10, are orthogonal in the initial configuration. In order to satisfy this condition for an arbitrary mid-line curve, the node generation is performed in the corresponding Serret-Frenet frame , , . Figure 12 shows a curve as a position vector r representing a mid-line of the rope. A circular cross-section of the solid beam finite element is modeled with three nodes. The approximation functions for nodes 1, 3, 4 of the solid beam finite element are defined in the local cylindrical coordinate system associated with the Serret-Frenet frame similar to follows Equation (47), therefore, orthogonality conditions are fulfilled.

Node generations for the solid-beam finite element
We show the node generation procedure step-wisely for the case of a helical rope with radius r rope laying on a cylinder with radius R, which will be used in the follow-up numerical example.
1. Coordinates of mid-line nodes for N el linear slave solid beam elements are offset by the value r rope from the cylinder surface in the radial direction. Afterwards, mid-line nodal coordinates are computed according to the helix Equation (13) with increments in the circumferential direction N el within the domain of contact angle ∈ [0, ]: 2. The Serret-Frenet basis vectors , , , see Section 2, are computed via slightly different from their definition procedure, which is computationally less expensive Setup of nodes for the solid-beam finite element: a mid-line for the rope and nodes for a cross-section in the Serret-Frenet frame • first, a unit vector is computed as normalized tangent vector via parameter : Both original and generalized Euler-Eytelwein equations have been derived for an unextensible rope or a thread. That is why the very important issue for further numerical investigations is the modeling of a rope using beam or solid-beam finite elements.
A "bending of a flexible beam by a rigid beam" test has been introduced in Konyukhov and Schweizerhof 26 in order to study a rope-like kinematic behavior compared with a beam-like kinematic behavior. In other words this test is served to answer unequivocally what is the difference between a rope and a beam in terms of their mechanical properties.
The test is configured as follows: a flexible weightless beam modeled with a selected set of finite elements-beam, shell, solid-shell or solid-beam finite elements-is going to be deformed by another rigid beam, modeled with a single finite element. The rigid beam is translated along the vector as shown in Figure 13 by step-by-step applied displacements. Corresponding contact elements should be specified between flexible and rigid beams.
During deformation process we distinguish between the pure Kirchhoff beam-like kinematic behavior, see Figure 14 left, and and the rope-like kinematic behavior, see Figure 14 right. It has been shown in Reference 26, that in order to observe the rope-like kinematics the beam theory allowing rotation of the normal (so far different from the Kirchhoff's beam) should be employed and the corresponding shear modulus G should be scaled as at least 1000 times less than the corresponding shear modulus for elastic beam G = E 2(1 + ) . A beam finite element based on the Timoshenko theory has been implemented in Reference 26. Further on it has been proved in Reference 19 for the "bending of a flexible beam by a rigid beam" test using the solid-beam finite element that the rope-like kinematic behavior can be reached for 3D continuum only by using an orthotropic material and by setting corresponding shear modulus G and G to much lower value in comparison with the isotropic material. It has been also found out by numerical experiments that either setting bending stiffness to much lower value or another combinations within elastic coefficients for the orthotropic material does not deliver "a rope-like" kinematic behavior.
Since a value of the scaling factor for shear modulus is not specified exactly in order to reach "a rope-like" kinematics, further on in the current contribution several computations with different scale factors will be provided for the verification.

A rope as a helix on a cylinder
The main goal of the current contribution-the analysis of the 3D generalized Euler-Eytelwein problem as a set of benchmark problems, which can be used for the verification of different Curve-To-Surface contact algorithms. We start here with the first benchmark problem-a rope as a helix on a cylinder, discussed in Section 2.2-providing comparison with various numerical algorithms varying different parameters. Figure 15 shows a finite element model for this case as a helix on a half cylinder which is modeled in ANSYS as specified in Section 4.1. In the current example radius of the cylinder R = 1.00 m and the radius of a rope r rope = 0.01 m, pitch for the helix H = 1.00 m and contact angle for a rope in the cylindrical coordinate system − 0 = .

F I G U R E 15
A rope as a helix wrapped over a cylinder. ANSYS model Figure 16 shows comparison between the numerical results of the force ratio T 0 T in the modeled problem using the penalty type regularization and the analytic solution (17). The friction coefficient is varied from = 0 (no friction) to = 0.7. As it can be seen the numerical solution is very close to the analytic for relatively low friction coefficient ∈ [0, 0.3]. The higher the friction coefficient is chosen, the larger the deviation of the numerical solution from the analytical one.
Remark 7 (Locking inside the contact element). The possible explanation for the undisclosed ANSYS numerical solution is the "contact locking" inside the contact element leading to more stiff behavior. This can be better understood later during the analysis using the CTSB contact algorithm.

CTSB contact algorithm. Variations of subintegration domains and approximations for the contact integral
Since the solid-beam finite element is providing a 3D continuum model for beams and ropes, first of all, the shear softening is modeled to reach a soft rope type behavior. The material model is Saint Venant-Kirchhoff linear elastic and orthotropic with regards to the corresponding Serret-Frenet frame. In order to achieve a rope like kinematics of deformation in comparison with the Kirchhoff beam, the shear modules G and G are taken to be consequently 1000, 10,000 and 100,000 times less than a corresponding isotropic module G , = E∕(1 + 2 ) with the Poisson coefficient as = 0.3. Elasticity modules in all direction are taken to be equal E , = E , = E , = E = 2.1E11 N m 2 . The angle of contact here is Δ = radians.
The following approximation types together with corresponding CTSB contact algorithms are implemented for further numerical analysis: • six-nodes Solid-Beam with linear approximations for the mid-line with the corresponding 8-nodes linear CTSB contact element as discussed in Section 3. Setting all parameters for verification, the critical force ratio T T 0 in Equation (17), Section 2.2, becomes: Variations of shear modulus to verify a rope like kinematics. Figure 17 shows the numerical solutions of the force ratio compared to the theoretical solution in Equation (73) by varying the helix H height from 0, which represents a circular configuration, up to twice the cylinder radius H = 2 ⋅ R = 2. The first three curves (marked as "linear" after the "theor sol" one) are given for the case of linear approximation functions and the scaled modules G subsequently as mentioned before. The Gauss integration formula with one subdomain m = 1 and two points N gp = 2 is taken (the simplest case without subdivisions). As it is observed from Figure 17, the accuracy of the solution increases with a reduction of the shear modulus G, therefore a soft shear rope kinematics can be effectively modeled.
Variations of numerical integration.
The next three computations in Figure 17 shows geometrical error due to selection of various integration scheme for the linear middle line of a finite element: 1. Lobatto formula with 2 integration points without subdivisions, or nominally only with a single subdomain m = 1; 2. Gauss formula with 6 integration points without subdivisions; 3. Gauss formula with N gp = 2 integration points with m = 3 subdomains; The Lobatto quadrature formula, compared to the Gauss one, always contains 2 integration points located at boundaries. In this case the slave nodes are exactly representing a helix geometry and included into computation with zero initial penetration. It can be observed that all combinations of the Gauss quadrature deliver more error and higher values force ratio, than the simplest Lobatto formula with two integration points. This is explained as "initial penetration effect" because the integration slave point for the linear element of the rope does not lay on the curved rigid solid beam master element-obviously the linear approximation functions are not able to represent the curvature of the rope sufficiently. Figure 18 shows a schematic representation of this problem. This leads to stiff behavior and obviously explains "contact locking" inside a contact element mentioned in Remark 7.
By using higher-order approximation functions such as Hermite-splines, the curvature of the rope can be approximated better than with linear approximation functions. Figure 17 shows also the numerical solution using Hermite-spline approximations and the Gauss formula with N gp = 2 integration points with m = 3 subdomains for the corresponding contact integrals. It is remarkable that the numerical solution using the 2-point Lobatto quadrature delivers less numerical error than using the 6-point Gauss quadrature. The reason for this is that the C1-approximations of Hermite type do not exactly describe the curvature of the rope and, therefore, the slave points do not lay on the cylinder, in contrast to the case of two points Lobatto quadrature formula. Again the "contact locking" inside a contact element is observed. Computational efforts for the cases 2 and 3 are equivalent-6 total computation for the segment, however the case with the maximal number of subdivisions 3 delivers slightly better results due to well-known effect of integration of nonsmooth function, see more study in Reference 21.
The very important conclusion of this paragraph is formulated as a following remark.
Remark 8 (Remedy for the "contact locking"). For the best result in verification all contact points must lay on the initial geometry. This can be reached exclusively with the Lobatto quadrature formula with two points. Further variations either with the number of subdomains or the number of integration points as well as employing the Gauss quadrature formula are not delivering better results. Therefore, the Lobatto quadrature formula with two points for both residual and tangent matrices within the contact element is the remedy for the "contact locking." All further computations employing the CTSB contact will be performed with the Lobatto quadrature formula with 2 points only.
The conventional ANSYS contact algorithm does not contain any further parameters such as integration type etc. Those options are represented only within the CTSB contact element.

4.5
A rope as the circle on a cone Figure 19 shows a schematic representation of a rope with circular configuration on a cone which is modeled in ANSYS.
In this computation, first, the angle * for the cone is selected based on the critical value for for a coefficient of friction crit = 0.1. According to the Section 2.3, the rope can be in equilibrium only for a coefficient of friction above the critical value, otherwise it slides in the dragging direction.
In due course, in all computations the coefficient of friction must be large than this critical value > 0.1 in order to satisfy equilibrium condition in dragging direction. Inverse situation < 0.1 has been correctly detected as the divergence in quasi-static numerical examples.

F I G U R E 19
Schematic representation of a rope as a circle-line wrapped over a cone, modeled in ANSYS (75) Figure 20 shows the numerical results of the force ratio T 0 T from both ANSYS and FEAP computations for 8-and 16-node CTSB contact element, compared to the analytic solution. Only Lobatto quadrature formula with two integration points has been used. Both penalty and Augmented Lagrange methods delivered the converged solution for the analyzed problem in ANSYS. As it can be seen, the friction coefficient is varied from = 0.11 (due to crit = 0.1) until = 0.7. In addition, it can be seen from the Figure 20 that the numerical results using the Solid Beam element with the corresponding CTSB contact element deliver small numerical error for the all range of coefficient of friction. The following material parameters have been used: elasticity modulus is E = 2.1E11 N m 2 , Poisson's ratio is = 0.3, shear modulus G is reduced as shown in figure than the elasticity modulus E, penalty parameter is N = T = 2.1E8 N m 2 .

A rope as a geodesic line on a cone
and the following equation for the verification: Figure 22 shows the numerical results for the force ratio For larger values the same "contact locking" effect for the contact element in ANSYS is observed, see Remark 7, and can be easily explained and eliminated withing CTSB contact element. Furthermore, Figure 22 shows the numerical results of the presented benchmark problem computed in FEAP using Solid Beam elements with linear and Hermite-splines (corresponding 8-and 16-node CTSB contact elements) approximation functions in the longitudinal direction. An elasticity modulus of E = 2.1E11 N m 2 and a Poisson coefficient of = 0.3 of the rope were applied in both cases. The shear modulus was reduced as (G = 10 −5 ⋅ E). Penalty contact was used with normal and tangential penalty parameter of N = T = 2.1E6 N m 2 . The Lobatto quadrature with 2 integration points has been used for the computation. The relative error of these numerical results with respect to the theoretical ones is laying below 0.25%.

A rope as a circle on a sphere
The last numerical example in the current benchmark set is a circular arch on a sphere discussed in Section 2.5. Figure 23 presents a schematic representation of a rope on a half sphere modeled in ANSYS. The circle is placed at an angle 10 • with the principal circle plane of the sphere with respect to the center. It is obvious, that this circle is not a geodesic line and the rope can be in equilibrium only for a coefficient of friction which is above a critical value crit , otherwise it slides in dragging direction. This critical value is calculated as: crit = | tan | = | tan(180 • − 10 • )| = | tan 170 • | = 0.176 .
The value of > 0.176 has been applied in order to satisfy this equilibrium condition. Figure 24 shows a comparison between the numerical results for the force ratio T 0 T using the penalty contact type and the analytic solution (35). The friction coefficient is varied from = 0.2 (due to crit = 0.176) to = 0.7. The "contact locking" effect in contact element is more pronounced for the double curved surface here.
The current example has been analyzed exclusively with ANSYS software. Nevertheless, several modeling strategies described in monograph 21 would be possible within the academical university software package FEAP: • A rigid sphere is meshed with solid-shell finite elements of quadratic order; a rope can be modeled either with beam or solid-beam finite elements of the various approximation order. A contact algorithm depending on the combination of the order of approximation for both sphere and beam/solid-beam might be chosen as the Node-To-Segment (NTS) or the Segment-To-Segment (STS) contact algorithm.
• It is possible to avoid the meshing of a rigid sphere if the Curve-To-Rigid Surface (CTRS) contact as a sub-set of the Segment-To-Analytical Surface (STAS) contact algorithm is employed. In this case one is free to choose either beam or solid-beam finite elements of the various approximation order in order to model a rope.
Due to a relatively large number of possible combinations we preferred to stay in line with our strategy in the current example, namely solid-beam finite elements for the exact geometry together with the CTSB contact (but not possible with this example) as well as ANSYS computations.

CONCLUSION
A new set of particular benchmark problems possessing closed form analytic solutions has been studied for the verification of the computational curve-to-surface/beam-to-surface/beam-to-solid contact algorithm based on the generalized Euler-Eytelwein formula, or the belt friction problem. This problem is formulated as equilibrium of a rope/string/belt laying under friction conditions on a rough arbitrary surface. The set includes the following geometrical combinations: a helix on a cylinder, a circle on a cone, a geodesics on a cone, a circle on a sphere. First important thing for the verification of the quasi-static solution is the position of a whole curve strictly within the equilibrium domain in the dragging direction on a surface in order to make the equilibrium this direction possible. Even partial positioning outside the equilibrium domain will cause disconvergence for a quasi-static solution as a result of impossibility to find an equilibrium position numerically. No separation condition-a curve is laying on a convex part of a surface-can be, however, constructively fulfilled, for example, by placing a string/rope/beam within a tube profile. In order to verify the force ratio T T 0 in the pulling direction applying beam/solid-beam finite elements material properties should be selected in order to satisfy shear softness: by setting a transverse shear cross-section parameter to low value for conventional beam finite elements; by taking into account orthotropic material with low value for transverse shear modules G , , G , with regards to the local Serret-Frenet frame. The benchmark problems have solved numerically using the commercial finite element software ANSYS and the developed CTSB contact algorithm for both linear and C1-continuous approximations. During verification process for the pulling force ratio T T 0 , a problem of "contact locking" in contact element showing stiff behavior in comparison with the analytical solution has been detected and studied using various numerical integration techniques of both Lobatto and Gauss types with subdomain subdivisions. Convergence in these cases, however, has not been affected. It has been found that application of the Lobatto quadrature formula with lowest order 2 for the integration of Curve-To-Surface contact integral gives the best results if the initial geometry of a curve is not represented exactly as a curve on a surface. The set of benchmark problems can be widely employed during developments of other novels algorithms for contact between beams, ropes, string and surfaces.

DATA AVAILABILITY STATEMENT
Data openly available in a public repository that issues datasets with DOIs.