A finite sliding model of two identical spheres under displacement and force control—part I: static analysis

A simplified analytical model of tangential contact engagement, sliding and separation of two elastic, identical spheres is developed assuming the kinematically induced sphere motion trajectory or load controlled sliding motion. The evaluation of driving force during contact sliding motion is determined for both monotonic and reciprocal sliding motion. The analytical formulae and diagrams of driving force versus sliding path are specified for linear and circular paths. The sliding trajectories are also determined for the load controlled programs. The results presented can be applied in the experimental testing of frictional response of contacting bodies, in a wear study of rough surfaces or in the contact interaction analysis of granular material during flow. The results can also be relevant for the development of the discrete element method widely applied in simulation of granular material flow, where the sliding regime conditions prevail in grain contact interaction.


Introduction
The contact interaction of two elastic spheres under normal or oblique loading and torsional couple has been studied by Mindlin and Deresiewicz [1], Lubkin [2], Deresiewicz [3], Walton [4], Vu-Quoc et al. [5], Segalman et al. [6]. For an oblique loading, only the contact slip regime was studied with a sticking (adherence) zone, in the central part and the slip zone, in the outer part of contact area. The sliding regime then occurs when the slip zone develops within the whole contact area. For a cyclic variation of the tangential force, the progressive and reverse slip zones are generated, and the hysteretic deformation response is accompanied by the consecutive evolution of loading, unloading and reloading slip zones from the boundary of the contact zone. The discrete memory of contact response with consecutive creation and erasure of loading events can be geometrically presented by means of consecutive loading surfaces in the T of two spheres under the fixed normal load and monotonically or periodically varying torsion couple can be analytically treated for both slip and sliding regimes as the contact zone remains fixed and the torsional couple reaches a limiting value, cf. [2,3,6]. In order to simulate oblique impacts with adhesion, Mindlin and Deresiewicz theory was extended by Thornton and Yin [10]. When the sliding regime develops under the fixed normal load N and the increasing tangential load T , the central sticking zone is erased, and sliding occurs along the whole area of contact. For a specified trajectory of the sphere center, both N and T vary, and the contact zone changes its size and orientation when moving along the boundary of the fixed sphere. This type of contact response will be analyzed in the present paper. In fact, the sphere displacement reached during the slip regime is very small relative to the radius of the contact zone, and contact interaction during the flow of granular material is associated for a large part with the sliding regime. Such factors, as the contact force evolution, length of sliding path and time period of the contact interaction and frictional dissipation during monotonic or reciprocal sliding, become important parameters in the deformation and flow analysis of granular matter. Both static and dynamic sliding motion will be considered in Parts I and II of this paper.
The paper is organized as follows. In Sect. 2, some basic formulae for the slip regime are presented by following Mindlin's and Deresiewicz's (M-D) analysis [1]. In Sect. 3, the sliding regime is analyzed for the imposed linear and circular trajectories of the moving sphere. In Sect. 4, the sphere sliding is considered for the specified load control process, while the concluding remarks are stated in Sect. 5.

Slip regime analysis
Consider two identical spheres of radius R made of a linear elastic material characterized by Young's modulus E and Poisson's ratio ν. The spheres are placed symmetrically with respect to each other in a mutual contact (Fig. 1). Applying the normal compressive force N to each sphere, the circular contact zone of radius r = a H is generated. When the tangential force T is subsequently applied and monotonically increases, the annular slip zone with radius c ≤ r ≤ a H develops at the contact perimeter with the slip displacement discontinuity oriented parallel to the force T and, in the central sticking (adherence) zone of radius 0 ≤ r ≤ c, the slip vanishes. When T reaches its limit value T = μN (where μ is the coefficient of friction between the spheres), the slip zone expands to the whole contact area, and the adherence zone vanishes, c = 0. Next, for the continuing deformation process, the sliding regime develops and is characterized by the motion of the contact zone, while its size evolution is terminated by the contact separation. If both N and T increase proportionally and T = μN , then there is no slip regime and the sliding regime takes place instantaneously.
The radius of the contact zone is specified from the Hertz solution: where K = 3 1 − ν 2 /(4E). The normal force induces the penetration depth (overlap) h 0 between two spheres in contact ( Fig. 1), specified by the relations: where k n = 1 K represents the coefficient of normal contact stiffness. In fact, the tangent and secant stiffness moduli are derived from (2).
The tangential displacement δ relative to the centers of contacting spheres ( Fig. 1) is expressed as follows from the M-D solution [1]: where G = E 2(1+ν) denotes the shear modulus. An ultimate displacement δ u at which the sphere starts to slide is specified by setting c = 0, that is, for the vanishing sticking zone. Thus, in view of (1) and (3), there is It is seen that the ultimate tangential displacement, at which the slip regime turns out into the sliding regime, is proportional to the overlap value h 0 or to the square of contact radius a H , depending also on the Poisson's ratio ν and the friction coefficient μ. Its value is significantly lower than the value of overlap h 0 .
In Fig. 2a, the dependence of the ratio δ u /h 0 on μ and ν is graphically illustrated, and the function of δ u /R versus a H /R, μ and ν is presented in Fig. 2b.
These diagrams (Fig. 2) provide the bounds on the ultimate tangential displacement resulted from the slip of spheres first loaded by normal force N , next by increasing force T . In fact, for most common materials Poisson's ratio is, 0 ≤ ν ≤ 0.5, the value of δ u is much smaller than the overlap value, and the ratio δ u /h 0 tends to zero for the vanishing friction coefficient. The slip regime of the contacting spheres takes place, when the ultimate tangential displacement is up to 3-16 times smaller than the overlap between the spheres, for μ ∈ [0.1, 0.6].
Consider now the contact engagement, typical for the flow or impact problems, when the contacting sphere moves with the velocity inclined at the angle ϕ to the normal direction of the contact. Assuming the slip displacement to be collinear with the velocity vector, we can write in view of (2) and (3) The slip zone specified by c is now instantaneously generated and varies with the angle ϕ, thus and a pure sliding regime occurs, when c/a H = 0, thus Figure 3 presents the impact angle ϕ u specifying the sliding and slip domains. It can be seen that the slip regime develops only for small values of ϕ u . Let us refer to the analysis of Di Renzo and Di Maio [11] who made an attempt to apply M-D theory to predict the evolution of the tangential force, displacement and velocity during the collision of an aluminum oxide sphere with a soda-lime glass anvil for varying impact angles. For fairly inclined impacts, the sliding mode dominates during the whole sphere motion period. However, for small impact angles, the sliding regime was observed only in the final short period of the sphere contact and was preceded by the longer period of slip regime for varying normal and tangential displacements.
The following conclusions can be stated in view of the presented analysis: • The M-D theory can be applied for the slip regime occurring in the small range of tangential displacement 0 ≤ δ ≤ δ u , where δ u is significantly smaller than the sphere overlap h 0 . Similarly, for the impact loading, the slip regime occurs only for small impact angles. • When δ > δ u or ϕ > ϕ u , the sliding regime occurs, for which the M-D theory cannot be applied and new relations should be derived for the contact force sliding trajectory evolution.  (Fig. 4) with the velocity v s (t 0 ). This sliding mode can be referred to specific cases, for instance, the spherical asperities interaction in relative translation of a punch over a rough substrate or to translation controlled testing of friction and wear of two interacting bodies. The linear path O 1 O 2 of the sphere corresponds to the contact path A-B-C, where the point A represents the contact engagement at time t 0 , B is the center of the sliding line for the contacting spheres corresponding to the maximal overlap h 0 , and C is the contact separation point occurring at time t c , provided the kinetic energy has not been earlier lost due to the frictional dissipation process (Fig. 4). Now, the sliding regime can easily be specified by requiring ϕ > ϕ u , where ϕ is the angle between the initial velocity of the contacting sphere to the linear path O 1 O 2 (Fig. 5). An ultimate angle ϕ u is specified v s (t 0 )>0 Linear path of the contacting sphere  The contact separation occurs for the paths emanating from O 1 in the exterior of slip and sliding domains. During sliding, the contact zone moves with respect to both spheres, changing its orientation and size (Fig. 4). In particular, at the contact engagement, the overlap is zero, while with progress of sliding, the overlap h t starts to grow up reaching the maximal value at the symmetry line, and subsequently, it diminishes tending to zero again at the contact separation.
The analytical solution within the continuum mechanics formulation is not available and only the numerical incremental procedure can be applied. Therefore, in the present paper, a simplified model is developed by assuming the normal contact traction N to be specified by the Hertz formula (2) in terms of the contacting spheres overlap geometry and material characteristics, while the tangential traction is defined by the sliding friction rule, T = μN .
The sliding path s t ∈ [±s u , ∓s u ] will be specified in an appropriate coordinate system x, y with the reference configuration of the system assumed to be a right-handed Cartesian with an origin point coincident with the center of the lower sphere, (Fig. 4). The sign conversion in s t bounds is resulted from the leftward or rightward motions. In particular, the upper signs are referred to the contact engagement at the sliding path coordinate x = +s u and its separation at x = −s u during the leftward motion of the contacting sphere. For the rightward motion (Fig. 4), the sliding path starts at position x = −s u , and ends at x = +s u , while the lower signs are accounted for. Instead of using the perpendicular projections onto x and y axis, the sliding path may be also analyzed in terms of angle α t ∈ [±α u , ∓α u ].
On the contact surface center, the slip displacement is zero, and the sliding displacement is defined in terms of the sliding path length s t ∈ [±s u , ∓s u ]. At the contacting sphere center, the contact evolution process is accompanied by both slip and sliding displacements, where the slip displacement varies according to formula (4) and is proportional to the contact overlap. In fact, when contact separation occurs, the slip displacement is erased. However, during the sliding reversal motion under preserved contact pressure, the slip process precedes the sliding mode and generates continuous reorientation of contact friction stress. The slip displacement is accounted for in the description of transition effects during the sliding reorientation process in Sect. 3.4. The problem of progressive and reciprocal sliding motion will be discussed in detail in this section.
3.2 Case 1: a linear path of the contacting sphere

Contact geometry
Let the upper sphere move along a linear path from point O 1 to point O 2 relative to the lower sphere, with its center fixed at point O (Fig. 4). The spheres come into contact at point A, and after a sliding time period t c , the separation of the spheres occurs at point C. During this period, the sphere center translation length is 4s u , while the contact sliding path is of length 2s u . Referring to Fig. 4, it can be seen that the overlap h t and contact zone radius a t developed at any s t are symmetrical relative to the spheres intersection points.
The contact condition between two identical spheres may be expressed as 2R − y > 0, where y = y (t) = const (for the linear path). Suppose the maximal overlap h 0 between the contacting spheres is given. The fixed vertical coordinate of the upper sphere may be set as y = 2y 0 = 2R − h 0 , while the sphere center position may be denoted by x = 2s t (Fig. 4).
The overlap h t corresponding to the sliding path coordinate s t , developed at time t, is specified from the formulae and For s t = 0, there is h t = h 0 , and, for s t = s u , there is h t = 0. The ultimate coordinate of the linear sliding path is expressed using (9) as follows: The ultimate coordinate of the upper sphere center specifies as 2s u (h 0 ) (Fig. 4), while the overlap dependency versus sphere center trajectory is simply determined by Eq. (9) substituting s 2 t for 4s 2 t . Mathematically, the sign ± implies the symmetry conditions, i.e., the linear sliding path developed at any t can be bounded over the range of s t ∈ [±s u , ∓s u ], depending on the direction of motion (as explained above).
The radius of the contact zone a t is specified from the equation where a t (h t ) > 0 should be specified in analysis.
Comparing (10) and (11), it can be seen that the contact zone radius a t changes in the same way as the ultimate coordinate of the linear sliding path. In fact, when h t = h 0 , then, The evolution of the contact overlap h t (s t ) and contact zone radius a t (s t ) over the linear sliding path is plotted in Fig. 6.
As can be seen from the functions plotted in Fig. 6, the increase in the maximal overlap h 0 with respect to the sphere radius produces the increase in the linear sliding path. With progress of the sliding path coordinate s t , the overlap h t grows up reaching the maximal value at the symmetry line, and subsequently, it diminishes tending to zero again at the separation of the spheres. The contact zone radius a t has the maximum at the maximal overlap h 0 and zero at the contact engagement and separation.
It is easy to note that the value of contact radius specified by (11) is larger than that defined by the Hertz formula (2). In fact, for the small value of the ratio h 0 /R 1, we have Figure 7 illustrates the main difference between the Hertzian contact zone radius a H and that determined by the overlap approach a u .
The Hertz formula results from the elastic solution accounting for the contact deformation between two bodies and a H are smaller than the intercept value resulting from the intersection of undeformed profiles in the overlap approach. In the limit, for h 0 /R → 0, the ratio a u /a H → √ 2.

Evolution of contact tractions
The imposed linear trajectory induces the contact forces acting at the center of the contacting surfaces between the spheres. Now, the normal and tangential forces acting on the contact surface ( Fig. 8) are related to the Undeformed profiles (Overlap approach) linear sliding path in terms of the overlap geometry by using the Hertz formula (2) and the sliding friction rule, thus

Fig. 7 Contact geometry interpretation using Hertz contact model
or in non-dimensional form Accounting for the contact normal vector rotation under the linear trajectory of the sphere along O 1 O 2 (cf. force directions at time instants t 1 and t, in Fig. 8), the normal and tangential contact forces can be referred to the nominal contact plane at s t = 0, corresponding to the maximal overlap h 0 . For a rightward or leftward motion, the resultant of the vertical and horizontal forces denoted by N * and T * can be expressed as follows: Finally, in view of formulae sin α t = s t s 2 t +y 2 0 and cos α t = y 0 s 2 t +y 2 0 , the above expressions can be expressed for the linear sliding path as follows: In Fig. 9, the normal and tangential force diagrams are presented for both sliding directions and different ratios of h 0 /R. Here, the unit values for k n = 1 N/m 3/2 , R = 1 m and μ = 0.6 were simply selected. They correspond to the non-dimensional characterization form specified by Eq. (15), so the same symbols are retained.
As can be seen in Fig. 9, the force diagrams are not symmetric with respect to the nominal plane s t = 0. The asymmetry is produced by the contact normal vector rotation along the linear sliding path. In particular, for s t ∈ [−s u , 0] and rightward motion, the horizontal component of normal force is pointing in the same direction as the friction force, but, for s t ∈ [0, s u ], this component is oppositely oriented to the friction force, (Fig. 8). The direction of asymmetry develops toward the direction of sliding velocity. In Fig. 9a, the direction of velocity is pointed to the right (i.e., positive); hence, the resultant force T * (s t ) is negative. The slip displacement range from M-D theory corresponding to the slip regime, when the spheres first enter into contact along the y-axis and next the slip and sliding motion occurs along BC or BA lines (Fig. 4), is marked in Fig. 9. It is seen that the slip displacement is much smaller than the sliding displacement.
3.3 Case 2: a circular path of the contacting sphere

Contact geometry
Consider now the contact response for an imposed circular trajectory of radiusR and eccentricityẽ with respect to the fixed sphere center O (Fig. 10). Suppose, the upper sphere moves with the initial velocity v s (t 0 ) along 10 Two spheres in the relative sliding motion along an eccentric circular path this path and comes into contact at point A, and after a finite time t c the spheres separate at point C (Fig. 10). As the upper sphere moves along the circle with a non-zero eccentricity, it slides with respect to the fixed sphere along a curved sliding path which can be specified in terms of the generalized coordinatess t ,α t andr t . In this case, the curved sliding path at any t deviates from a perfect circular trajectory. All variables referred to the circular trajectory will now be marked by the tilde superscript.
The evolution of the curved sliding path can be defined in the polar coordinates (r t ,α t ), wherer t = const and varies with the angleα t ∈ [±α u , ∓α u ] at any t. The contact condition between two identical spheres may be expressed as 2R −R 1 > 0, whereR 1 =R 1 (t) = const.
In view of the relationR the radius of the curved sliding path can be expressed as follows: wherek =ẽR is defined as the eccentricity ratio.
Forα t = 0, we obtainR 1 = −ẽ ±R. Hence, following Fig. 10, the positive sign should be taken into account. The derivative ofr t (α t ) with respect toα t has the form: Due to the symmetry, the curved sliding path may be described by integration ofr t (α) overα ∈ [0,α t ], thuss Assumingk < 1 and integrating equation (23), it is obtained Thus, the evolution of the overlap versus angle for two identical spheres in contact undergoing relative motion along a circular path with eccentricity may be defined as Solving Eq. (24) forα t , we arrive at the following formula: where Inserting Eq. (26) into relation (25) and using the positive sign in (26), we can specify the distribution of the overlap over the curved sliding paths t ∈ [±s u , ∓s u ] as follows: To rearrange (28) into a closer form, we refer the curved sliding path to the rectangular coordinates (x t ,ỹ t ), (Fig. 10). In this case, defining the path radius by the relatioñ the contact overlap can be expressed by the formulã In view of the relatioñ and Eq. (20), we can obtain:ỹ Inserting this relation into expression (28) (using the positive sign), we can rewrite the formula for the contact overlap in the following form: Now, the arc length of the curved sliding path is expressed by the integral, thus Integration of (34) is simple, thus yielding the following formula: The inverse form of relation (35) is as follows: Inserting (36) into the relationship forh t (x t ), we finally express the contact overlap in the form The following identities result from formulae (25) and (37): The contact zone radiusã t , (Fig. 10), is specified by the relation (11), thus Finally, the formulae for the ultimate values of the above functions should be derived. The contact engagement condition ish t (±α u ) = 0, and for the separation of spheres, it turns out toh t (∓α u ) = 0. Thus, substituting one of these conditions into relation (25), we obtain a limit value of the angle specifying the curved sliding path:α The ultimate value specifyings u is simply calculated from Eq. (24), for theα t =α u . Thus, the presented formulae specify the curved sliding path geometry for the circular trajectory with eccentricity of the contacting sphere, when the boundss t ∈ [±s u , ∓s u ],α t ∈ [±α u , ∓α u ] are fulfilled.

Interpretation for the eccentricity ratio
The limit values of eccentricity ratiok =ẽR may be interpreted as an index for transition from the perfect circular trajectory to the linear path. Thus, for the perfect circular trajectory,h t (α t ) = const, and this condition is fulfilled whenk = 0. In this case, from Eq. (39) and (21), we getα u = ∞ andR 1 =R = const, meaning that the contacting sphere slides along the perfect circular trajectory.
For the linear path, the vertical coordinate stands forỹ t = y t = const. In this case, the ultimate coordinate of the curved sliding path should be approached as close as to the ultimate coordinate of the linear sliding path, i.e.,s u ≈ s u ≈x u . To prove this condition, Eq. (35) may be rewritten in the following form: Here, the solution fors u may be obtained only when the angles are small. This yields thats u ≈ s u , for 2s u /R 1, andk → 1, fork =R/ẽ. Hence, the transition from the circular path with eccentricity to the linear trajectory of the contacting sphere will occur whenk → 1. In this case, it is convenient to select a certain value of the maximal overlap h 0 , assume the ratio ofk to be close unity, and finally predict the eccentricity and radius of the trajectory by the formulae e = (2R − h 0 )k 1−k ,R =ke. In Fig. 11, the graphs of transition from the linear to circular paths are plotted in terms of different eccentricity ratios.
As can be seen in Fig. 11, whenk ≈ 0, the contacting sphere slides along the perfect circular trajectory over the other sphere. For 0 <k < 1, the motion of the contacting sphere occurs with different eccentricity ratios resulting in the different lengths of the curved sliding paths. Finally, fork → 1, the contacting sphere trajectory approaches a linear path.

Evolution of the contact tractions
Similar to the case of a linear path, the contact force vectors acting on the curved sliding path at time t appear to be non-coincident with those acting at the overlap center. In particular, the force N (s t ) specified on the contact plane is not normal to the curved sliding path but has the projections N τ (s t ) , T τ (s t ) on this path (Fig. 12). As can be seen here, the forces acting on the curved sliding path and those on the contact plane can be related by the angleβ t .
As follows from Fig. 12, the sum of force projections onto theτ andñ axes provides the resultant tangential and normal forces for a rightward or leftward motion. Thus

Fig. 12 Contact force vectors acting on the contact surface at time t under a circular path with eccentricity
It should be noted that these expressions are the same as relations (16), (17), forβ t = α t . Now, referring to Fig. 12, the relation between anglesβ t andα t may be determined asβ t =α t −γ t .
In order to get the formula for the angleγ t , the following relationship can be used: An inverse form of the integral equation (34) yields The derivative for the arc length of the curved sliding path specified by Eq. (35) can be expressed as Thus, using Eqs. (43-45), we obtaiñ In view of formulae (36) and (46), the angleγ t may be related tos t as: Now, using Eq. (26), the sine and cosine functions in Eqs. (41), (42) can be expressed in the following form: Finally, the contact forces versus trajectory angleα t are expressed in the form: The correctness of the obtained relationships can be demonstrated. Thus, substitutingk → 1, we get α t ≈ α t , and formulae (50), (51) turn out into expressions (16), (17) for the linear sliding path. This transition was already demonstrated in Fig. 11.
Also, it can be demonstrated that Eq. (51) specifies the contact traction, corresponding to the perfect circular trajectory of the contacting sphere. From Fig. 11, it is clear that for the perfectly circular path only the frictional force should arise, and the overlap remains constant. Substitutingk = 0 into (51), we simply obtain T * (s t ) = T (s t ) = −s g μ |N (s t )|. Hence, for the perfectly circular path, only the sliding friction force contributes to the response of two identical spheres in contact.
Finally, for the curved sliding path expressed in terms of s t ∈ [±s u , ∓s u ], the relations (50), (51) can be specified as follows: where is the normal force acting on the sliding contact zone, expressed in terms of the overlap geometry and Hertz formula (2). In Fig. 13, the evolutions of the contact tractions for the rightward or leftward motion, illustrating the effect of different eccentricity ratiosk, are plotted. Here, the unit values of k n = 1 N/m 3/2 , R = 1 m and μ = 0.6 were also simply selected for the illustration.
As can be seen in Fig. 13, for the eccentricity ratiok → 1, the circular and linear paths tend to each other. The force-path functions are non-symmetrical, and the asymmetry grows along the direction of the sphere motion. The normal force projection also contributes to the slow down of the sphere motion along the induced trajectory, similar to the friction force action. Both friction and contact configuration effects interact during sliding, thus affecting the required driving force. A rigorous theory evaluating the slip regime of two elastic spheres under the fixed normal load and monotonically or periodically varying tangential force has been presented by Mindlin and Deresiewicz [1]. In this case, the slip mode occurs for very small sphere centers displacement relative to its radius with zero slip at the contact surface center. In the discrete element method (DEM), the slip regime is modeled in a simplified manner. In this case, the tangential force acting at the contact center is usually calculated accounting for the linear dependence between the slip displacement and tangential force of a spring whose stiffness is related to the overlap by M-D theory. The slip displacement vector, in turn, is determined by integration of the contact tangential velocity versus time yielding a small motion of the sphere and resulting in a non-zero displacement relative to the contact center. This approach treats the spheres rather more as rigid than deformable neglecting the existence of an adherence zone, which follows from M-D theory, and allows us for the simplified analysis.
For a specified sliding trajectory of the contacting sphere, both contact tractions, N and T , vary, and the contact zone changes its size and orientation. In this case, an explicit account for the unloading/reloading regimes within the continuum mechanics formulation is not available. Therefore, two simplified models will be discussed here.
First, consider the unloading/reloading regime simply neglecting the slip mode. In this case, the loading curve instantly drops onto a reversal sliding path, while the reloading curve turns out instantaneously onto the loading path. In this case, the unloading curve, bounded for s A ≤ s unl ≤ s C , and the reloading function, for s E ≤ s rld ≤ s C , are discontinuous at the unloading and the reloading starting points T A and T C (Fig. 14).
The discontinuous unloading/reloading functions can be generated in terms of the Heaviside step function for a discrete variable, thus  Secondly, let us account for the slip displacement at the start of unloading or reloading regimes. In this case, the unloading curve attains the reversal sliding path without the discontinuity jump, i.e., throughout the force-slip displacement curve (Fig. 14).
As demonstrated above by M-D theory, the ultimate slip displacement, at which the sphere starts to slide, is small relative to the overlap. Hence, the force-slip displacement curve can be described approximately by a simple power law function. Thus, assume that the unloading force-slip displacement curve can be defined by the function where f ( s) = s p , s is the increment of the slip displacement, p is the power exponent, (0 < p < 1), and k H is the slope of the curve. The unloading curve slope k H = | is the value of the load at which the slip curve intersects the reverse sliding curve, and δ u (s A ) is the ultimate slip displacement at s A (Fig. 14).
Finally, the unloading/reloading force-path curves can be generated as follows: In general, the ultimate slip displacement, when the slip regime turns out into sliding regime, should be dependent on the radius of the contact zone and the point at which the unloading or reloading regime starts. In view of Eqs. (4), (9) and (11), we can obtain In Fig. 15, the relations of the ultimate slip displacement versus sliding path for different values of the overlap are demonstrated. As can be seen from these plots, the maximal values of δ u (s t ) predicted by (60) have  Fig. 16 for the linear trajectory of the contacting sphere.
As can be seen here, for a full sliding path, i.e., s t ∈ [±s u , ∓s u ] along the linear trajectory, the relation T * = f (N * ) is represented by the (elliptical like) force loops originating from the contact engagement point and terminated at the contact separation point, where the contact tractions vanish. For the leftward motion, the loop develops in a positive (N * , T * ) plane, while for the rightward motion the function values are negative. With increasing the overlap between the contacting spheres, the loops become more elongated. The force loops represent the configurational effect associated with the contact plane rotation during sliding.
The force-path diagrams for a partial sliding combined with the reverse motion inducing loadingunloading-reloading regimes with/without account for the slip are plotted in Fig. 17. Here, the contact tractions were computed by Eqs. (55-62) for the linear trajectory of the contacting sphere. The power exponent value p = 0.5 was applied in Eqs. (58), (59).
In Fig. 17, the contacting sphere undergoes a leftward motion (v s (t 0 ) < 0) and comes into contact at the point position s u . Later, at the point s t1 < 0, the direction of velocity changes, inducing the reverse motion. In this case, the unloading regime develops and the tangential force curve suddenly (denoted by dashed lines) drops down or following the slip curve gradually moves along the reversal sliding path. The unloading regime ends at the sliding path coordinate s t2 > 0, and the reloading regime develops under the negative velocity. According to Eq. (60), the reloading slip curve has a smaller ultimate slip displacement value than the unloading curve, because of the inequality |s t2 | < |s t1 |. The main features of the application of the unloading/reloading curves are also illustrated in detail during the dynamic analysis (Part II of the current study).

Analysis of finite sliding of rigid spheres
Referring to Fig. 18, consider the case of relative sliding of two rigid spheres under the applied driving force T 0 . The lower sphere is fixed in its position, and the upper sphere initially rests on the horizontal plane at position y 1 . The normal vector of contact is located at point −x u and is inclined to the y-axis at the angle −α u . The gravity force Q is initially equilibrated by the supporting plane reaction. When the horizontal driving force T 0 is applied, the sphere separates from the plane and the force equilibrium is reached for the developed contact normal and tangential forces N and T . The sliding regime develops instantaneously, when the driving force T 0 reaches its limit value |T 0 | = T μ (where T μ is the resultant horizontal force). When the spheres are rigid, the sliding path lies on the fixed sphere surface, and for plane motion, it is a perfect circular path of radius R (Fig. 18). The initial orientation angle α u of the contact normal vector is specified as α u = arccos (y 0 /R) (where y 0 is the initial vertical coordinate of the contact engagement), while the support elevation is defined as During the sphere sliding, the orientation angle α t of the contact normal vector varies. The positive driving force T 0 > 0 induces the rightward motion with the positive velocity v s > 0, and the sphere path starts at the angle α t = −α u (Fig. 18). For T 0 < 0, a leftward motion (v s < 0) is obtained with the sliding path starting point located at α t = α u . It represents the configurational variable evolving during the sliding process, and the equilibrium for the sliding sphere path (Fig. 18) combined with the limit friction condition μN provides the following relations: where s g = sign (v s ) , Q = |Q|. The solution of (63) simply yields where ϕ = arctan μ, N (α t ) is the positive contact reaction (18), while T μ (α t ) is the limit friction force which should be reached to induce the motion of the upper sphere. The relation (64) can be simply illustrated in the plane T 0 , Q (Fig. 19). The limit friction for the fixed angle ϕ and varying angle α t is represented by the Coulomb cone with its axis oriented at the angle α t to the Q-axis. The initial contact is presented in Fig. 19a, for the contact orientation angle α u . During sliding, the orientation angle decreases and becomes negative; so, the contact separation point P s is reached, for α f = −φ (Fig. 19b). From relation (64) and Fig. 19, it also follows that the sliding regime can occur only when α u ≤ π 2 − ϕ. For α u = π 2 − φ, Eq. (64) yields T 0 → ∞ and the sphere cannot be moved by the applied force. Figure 20 presents the functions of driving and normal forces calculated by Eqs. (64) and (65) versus the curved sliding path, s t = α t R, for the different values of friction coefficients. In fact, for the reverse sliding, the driving force at the onset of sliding is represented by the point P r in Fig. 19b, and next the force varies depending on the evolution of the orientation angle α t . Here, the values of Q = 1 N, R = 1 m and y 0 /R = 0.7 were simply selected as the input parameters.
As can be seen in Fig. 20, before the contact engagement, the frictional force acts between the upper sphere and support, and at the contact engagement point (i.e., s t = s u ), the static friction force instantly grows up to the maximal value. This is due to the assumed rigid sphere contact response. For continuing sliding, the static friction force decreases until the contact separation occurs. As demonstrated in Fig. 20, the static friction force and the contact reaction forces increase for increasing friction coefficient.

Analysis of the finite plane sliding of deformable spheres
Consider the contact of two elastic deformable spheres. Suppose the upper sphere to be initially placed at any point x u with contact vector inclined to the y-axis at the angle α u and the constant vertical load Q is applied first, and next the horizontal force T 0 is increased monotonically. When |T 0 | ≤ T μ , it can be stated that the slip regime occurs and the resultant tangential force induces a small tangential displacement δ between the contacting sphere centers during contact engagement. When the horizontal force reaches the friction force T μ , the slip regime passes into the sliding regime and the sphere starts to slide over the contact surface, when |T 0 | > T μ . For both regimes, the force equilibrium along the τ -axis, (cf., Fig. 11, in Part 2), can be expressed as follows: where is the contact reaction force, and T is the resultant tangential force. For the slip regime, substituting T and N forces from Eqs. (66) and (67) into Mindlin's Eq. (3), the tangential displacement is defined as: where s κ = sign (Q sin α t + T 0 cos α t ) is the sign specifying the tangential displacement direction due to the resultant tangential force.
Using α t = 0 in Eq. (68), we simply arrive at the Mindlin equation (3). For T 0 = 0, the upper sphere center slides down along the contact plane under the action of vertical force Q, when the support is removed. Also, Eq. (68) yields that the upper sphere should be placed in a position specified by the initial angle |α u | ≤ φ.
Actually, besides the variable force T 0 , δ is also dependent on α t in the right hand side of Eq. (68), because the change in tangential displacement yields the change in the resultant tangential force. However, the above analysis has been shown that the limit values of the tangential displacements δ u are very small (see Sect. 4.3) and the constant value of angle α t ≈ α u can be reasonably substituted in Eq. (68), for modeling real materials.
When |T 0 | = T μ , substituting Eq. (65) into Eq. (68), the approximate equation specifying the ultimate tangential displacement is obtained. The root of this equation cannot be derived explicitly and is expressed in terms of the arc length at the centers of contacting spheres as follows: The slip regime terminates on the curved sliding path at the point s δ,u ±s u , when the friction force, specified by Eq. (64), takes the value T μ s δ,u ± s u /(2R) .
For the sliding regime, substitution of μN in Eq. (66) results in the solution given by Eq. (64). For the case of deformable spheres, this equation can be resolved accounting for the curved sliding path specified by the variable radius r t (α t ) depending on the developed overlap (Fig. 18). Thus, according to the Hertz solution and Eq. (65), the overlap can be expressed as follows: The radius r (α t ) of the curved sliding path may be expressed as the function of the overlap: where c = 2 is taken for the curved sliding path specified at the centers of contacting spheres and c = 1 stands for the path taken at the contact surface.
Starting from the expressions (70), (71), the relation of r (α t ) to the applied force Q is obtained, thus In this case, the vertical coordinate of the sliding paths is expressed as follows: The curved sliding path for deformable spheres can be defined by the integral Inserting Eq. (72) into Eq. (74), the solution for s t cannot be expressed explicitly. Hence, the function of the friction force against the sliding path can be defined in incremental form. Thus, Eq. (74) may be simply rearranged as follows: In this case, the angle increment can be expressed by the following relation: Combining Eq. (76) with (77), we obtain the tangential compliance in the following form: For the i-th iteration, the frictional force against the curved sliding path can be expressed as follows: In Fig. 21, the horizontal force versus curved slip/sliding path predicted for the sphere centers is plotted. For the slip mode (Fig. 21a), the constant vertical load Q is applied first, and the horizontal force T is increased subsequently up to the limit value T μ s δ,u ± s u /(2R) . In this case, the tangential displacement was predicted according to Eq. (68) for each step of T . The sliding regime (Fig. 21b) developed after the slip mode is calculated incrementally according to Eq. (79). The value of the vertical load Q = 100 N is selected, the contact stiffness k n values used are: 1,000, 2,000 Nm −3/2 , while the other parameters are the same as in the previous analysis. Also, the support at the contact engagement is not considered because the upper sphere is placed in a position specified by the initial angle of |α u | ≤ ϕ. This allowed to increase the horizontal force from zero value. In the case of support application, the initial value of the applied force should be used as Qμ.
As can be seen in Fig. 21a, the monotonic increase in the horizontal force produces the nonlinear increase in tangential displacement developed at the center of the upper sphere. This process terminates at the intersection point between the slip and sliding modes, and next the sliding friction force gradually decreases along the progressing sliding path. For increasing contact stiffness, the deformable sphere response tends to the rigid   Fig. 21 are presented in the diagram of evolution of the contact force N and driving force T 0 depicted in Fig. 22.
During the monotonic increase in the horizontal force T 0 , the slip regime develops from the point P 0 = Q cos α u (where T 0 = 0), resulted from Eq. (67) (Fig. 22). Due to relation between the tangential displacement δ and angle α t , the N − T 0 curve is slightly nonlinear. When the driving force reaches the friction force limit, i.e., |T 0 | = T μ , the slip regime terminates at the point P s and the sliding regime starts. During this process, the friction force decreases with progressing the sliding path according to Eq. (79), and the N − T 0 curve reaches the point P 0 again, where T 0 = 0. The unloading and reloading processes emanating from the slip regime are also shown in Fig. 22 by dashed lines.
When |T 0 | > T μ , the force equilibrium at a certain point of the sliding path is violated yielding N = 0, and the contact separation occurs. This case will be discussed in Part 2 devoted to the dynamic analysis.

Verification of the model prediction by experimental measurements
A theoretical prediction of the friction force for two rigid contacting spheres according to Eq. (64) will be compared with the experimental measurements. Łukaszuk et al. [12] measured the coefficient of friction between pairs of the metallic spheres and pea grains. Spheres and grains of 5 mm diameter were tested in a specially developed apparatus to measure the components of small forces and displacements. In the measurement, the simultaneous measurement of the friction force was induced to maintain slow steady motion, when the force Q is applied. A low value of the vertical force Q equal to 0.98 N was applied. The measured coefficient of friction μ was approximately equal to 0.68 (steel spheres), 0.22 (lubricated steel spheres) and 0.4 (pea grains). The initial position of the contact engagement y 0 was derived from the experimental graphs [12], for the vertical and horizontal displacements. In Fig. 23, the measured and predicted friction force versus horizontal displacement of the sliding sphere center are depicted.
As can be seen in Fig. 23, the theoretical and experimental friction force evolutions versus the horizontal sliding path component are sufficiently adequate. However, at the contact engagement, the experimental and theoretical predictions differ. The small value of force Q applied to a stiff material of the spheres resulted in small values of the theoretical ultimate tangential slip displacement of s δ,u = 0.4 × 10 −7 m and s δ,u = 0.7 × 10 −5 m for the steel and pea grains, respectively. In view of this fact (Fig. 23), the slip regime curves are almost vertical, and their values are smaller than the accuracy level of measured displacements. The difference in the maximal driving force prediction during slip regime may be attributed to the initial material and contact surface shape imperfections, combined with the development of localized inelastic deformations. On the other hand, during sliding regime, the model prediction is fairly close to the experimental data. Let us note that the precise measurement of the static force is quite complicated, depending on the technique of mounting the spheres in the initial contact position and application of vertical force.

Conclusions
A finite sliding model of two identical spheres under the displacement and force control has been proposed. The model is based on the approximate analysis of the contact interaction process of two spheres by the overlap approach along with the Hertzian nonlinear elastic spring analogy and Coulomb law for the sliding friction. In particular, the analytical and incremental (in the case of force control) relationships between the contact force and the sliding path length were derived. The effect of unloading/reloading was also analyzed for the case of reciprocal sliding motion.
The well-known Mindlin and Deresiewicz theory may be applied only for the partial slip regimes occurring for a very small range of tangential displacements induced by the motion of contacting sphere centers. For the sliding regimes, the M-D theory is not appropriate, and the present analysis provides the complete set of relations specifying the sphere motion.
Further applications and extensions of the analysis can be envisaged to other problems, such as (a) rough contact surface interaction, (b) wear analysis and (c) specification of the flash temperature effect in the contact zone.