Modeling of combined slip and finite sliding at spherical grain contacts

The present paper is aimed at developing the analytical description of the interaction of two contacting spheres for several classes of slip and sliding trajectories, typical in the experimental testing. The analysis accounts for memory effects in the slip regime and configurational effects in the sliding regime, expressed in terms of an active loading surface and memory surfaces within the space of contact forces. Analytical relations for contact response are derived for linear and piecewise-linear motion trajectories of the sphere. The problem of multiple contact interaction of the sphere moving over the regularly packed granular bed is also considered analytically. It is demonstrated that the dual contact activation-separation processes occur within the combined slip–sliding modes, essentially affecting the distribution of contact tractions. The results obtained are relevant for the class of contact problems requiring analysis of interaction of slip and sliding displacements, in particular in testing grain contact interaction aimed at specification of elastic, frictional and wear parameters.


Introduction
The deformational response of soils, rocks or granular materials has usually been described by constitutive relations within the formalism and rules of the theories of elasticity and plasticity. However, the state variables of inelastic response are related to micromechanical contact interaction. The models of contact slip and sliding then provide an input for the structure of constitutive models and evolution rules of state variables. Applying homogenization within a representative volume with account for micromechanical contact parameters, the incremental relations can be generated at the macro-level. The static and dynamic response of contact interfaces therefore becomes a fundamental issue in development of micromechanical approach.
A landmark study on contact mechanics was presented by Hertz [1], who provided the analytical formulae for normal static contact interaction and for impact of elastic bodies. Hertz model is based on the assumptions of frictionless, flat contact surface, linear elasticity, small strain theory and material isotropy. The FEM-based verification has matched the Hertzian solution exceptionally well under these assumptions [2]. Moreover, for the various realistic contact scenarios, which are theoretically not covered by the Hertz' model (frictional contact, non-flat contact surface), no significant effect on the Hertzian force-deflection relationship has also been observed, only the large strains have been found to cause important prediction errors [2].
The contact response of two elastic spheres under normal or oblique static loads has been later extended by Mindlin and Deresiewicz [3]. Relying on these fundamental works, Cundall and Strack [4] pioneered for introducing the numerical discrete element method (DEM), where the granular material has been treated as a conglomeration of solid particles undergoing a time-dependent dynamic motion subjected to the resultant action of forces, including the contact tractions and contact torques, according to second and third Newton's laws within the rigid boundary domain. The method has been recognized as a powerful modeling approach specifying direct introspection of the particle micro-properties at steady or transient deformation states without any a priori global assumptions, usually required for the continuum-based modeling [5][6][7][8].
The grain-grain contact interaction requires definition of tangential compliances. There are two fundamental deformation regimes, namely, the partial slip and sliding, developed consecutively in the progressive motion. For the reciprocal motion, these two regimes develop consecutively during each grain motion reversal. A third regime can occur for the particles bonded by the adhesive forces which suffer the damage and destruction, induced by progressive slip and sliding, thus, leading to a purely frictional response. It should be noted that the problem of multiple contact interaction with a simultaneous contact separation-activation of a sliding sphere over the regularly packed granular bed was not treated analytically earlier in the literature. In this case, after the first contact activation accompanied by slip, new contact activation with another sphere from the granular bed can occur. The following slip and sliding or their mixed modes induce the contact state evolution requiring the detailed analysis.
A well known Mindlin-Deresiewicz (M-D) theory can only be applied for the slip regimes occurring for a very small range of tangential displacement between the contacting sphere centers. For identical spheres under oblique elastic contact, it has been found that initially circular contact zone remains circular after being tangentially loaded [1]. However, the deviation from a circle to an oval shape was demonstrated by Larson and Storåkers [9], in a creeping oblique contact problem. In particular, this deviation was modest, while in the most extreme case (perfect plasticity, full adhesion, obliquity of π/3) it reached six percent. The contact interaction of spheres including the cohesive/adhesive forces has been also analyzed in [10,11].
For an oblique contact of elastic identical spheres under the proportional loading, Walton [12] determined that the normal and tangential directions are reciprocally independent. However, for a non-proportional loading, Elata and Berryman [13] found that the slip regime has to be taken into account, even for identical spheres, and solution becomes path-dependent. The semi-analytical model for oblique impacts of elasto-plastic spheres was proposed and verified by FEM and experiments in [14]. In [15,16], analyzing the relative sphere motion along a specified linear trajectory without the initially activated contact zone, it was demonstrated that both normal and tangential forces vary proportionally, generating the continuing evolution of the sliding regime without the initial slip mode.
For a cyclic variation of tangential force, the progressive and reverse slip zones are generated. In modeling such complex phenomena, the slip memory rules expressed in terms of the loading surfaces were firstly proposed by Dobry et al. [17] and next developed in a more detail by Jarzębowski and Mróz [18]. The alternative slip memory diagrams for the oblique loading were also presented in [19].
Thus, in the solution of these problems, the proper specification of contact characteristics of two grain-grain interaction is very important. The static testing of grain contact in combined slip and sliding regimes, aimed at specifying the friction coefficient, was performed by Lukaszuk et al. [20]. Two grain frictional contact response under monotonic and cyclic loading in the slip regime was also experimentally investigated by Cavarretta et al. [21] and Cole et al. [22,23]. The reference can also be made to fretting and sliding wear testing for the reciprocal sphere motion of amplitudes corresponding to the slip or combined slip-sliding regimes, cf. Fouvry et al. [24], Dini et al. [25].
The frictional dissipation, generated during the relative slip or sliding motions on the highly stressed asperity contacts, induces the flash heating at the elevated temperature affecting the values of stiffness moduli and friction coefficient. The induced contact instability along rock fault plane due to highly localized slip by such thermally driven weakening mechanism has recently been analyzed by Rice [26] and Platt et al. [27]. Also an attempt was made by Ciavarella [28] to extend Mindlin's solution by introducing Griffith's fracture model for the inception of slip regime. It modifies the area of adherence zone and induces increase of friction force before sliding can be developed.
The present paper is aimed at developing the analytical description of interaction of two contacting spheres for several classes of slip and sliding trajectories, typical in the experimental testing. In Sect. 2, the slip regime and its transition to sliding regime is briefly discussed. In Sect. 3, the combined slip and sliding modes are analyzed for a linear motion trajectory of sphere, entering in contact in the sliding mode or starting its motion from a formed contact zone. Next, in Sect. 4, the sphere slip and sliding along a regularly packed granular bed is considered.

Slip regime analysis
The contact interaction of two elastic spheres under normal and oblique loading has been studied by Mindlin and Deresiewicz (M-D) in [1], where two identical spheres were placed symmetrically in a mutual contact and subjected to normal compressive force N with subsequently applied tangential load T (Fig. 1).
A classical distinction between the slip and slide terms should first be mentioned. In particular, when the tangential force applied is primarily small, there is no slip on the entire contact surface. With increase in the magnitude of this force, the slip starts developing at the edge of the contact zone, since the concentrated shear stress has been distributed here in the absence of slip. Thus, the slip ring of radius c ≤ r ≤ a H develops at the contact zone perimeter ( Fig. 1a) with the slip displacement discontinuity oriented parallel to this stress ( Fig. 1a), whereas at the center, the adherence zone of radius 0 ≤ r ≤ c is created. Due to presence of this zone, the center point of contact surface remains fixed under the evolution of the displacement δs at the center of the sphere. For T = μN (where μ is the coefficient of friction between the spheres), the slip zone expands to a whole contact zone and the adherence zone radius vanishes, c = 0, inducing a dissipative sliding process. Thus, the slip regime is characterized by T < μN and 0 ≤ r ≤ c. Next, in a continuing deformation process, the sliding regime develops, since the sphere starts to move until contact separation. If both N and T increase proportionally and the limit friction condition occurs, T = μN , then the sliding regime is generated instantaneously within the whole contact domain. For periodically varying force T , the consecutive slip zones develop from the contact zone boundary and the hysteretic slip occurs (Fig. 1b).
Theoretically, the loading and unloading processes induced by varying force T can be represented by the respective loading surfaces in the N , T -plane (Fig. 1c, d). The normal loading surface F 0 = 0, the slip loading and unloading surfaces F 1 = 0 and F 2 = 0 correspond to the initial normal loading and the subsequent tangential loading, and unloading occurring within the domain bounded by the limit sliding surface F L = 0. Thus, the tangential force during cyclic loading inducing progressive and reverse slip depends on the accumulated slip zones in the contact domain, can be expressed by a function where c 1 , c 2 , . . . , c n are the radii of internal slip zones generated in consecutive loading events. They constitute the set of slip memory parameters of contact loading history and are erased in a consecutive sliding regime.
From Hertz contact law, the force N induced overlap h 0 between two spheres of radii R 1 = R 2 = R ( Fig. 1) can be determined by the following relationships: where k n is the sphere stiffness modulus in normal contact direction, K = 3 1 − ν 2 /(4E), in which ν and E are the sphere material Poisson's ratio and Young modulus. The radius of contact zone is then specified by the Hertz formula: For the created contact zone with N = const, the ultimate displacement between two identical spheres centers δs = δsu, at which the sphere starts to slide can specified by the following formula where G = E/2/(1 + ν) is the shear modulus. Consider now these relations for the spheres of different radii R 1 and R 2 and different linear elastic materials with Poisson ratios ν 1 and ν 2 , and Young moduli E 1 and E 2 . In this case, the effective moduli and the representative sphere radius are used [29]: For the identical spheres, the effective characteristics (5) can be rewritten as follows: Substituting these relations into (2) and (3), the normal force and the overlap formulae are as follows: For two different spheres, in view of (4) and (6), the ultimate slip displacement δsu is as follows For the spheres of the same material, but of different radii the following expressions can written For identical spheres, when the contact zone is created first, the δsu is extremely low, i.e., 3-16 times lower than the overlap h 0 , for 0 ≤ ν ≤ 0.5 and μ ∈ [0.1, 0.6]. Thus, to obtain the slip regime, the displacement condition δs ≤ δsu should be valid, otherwise the sliding regime occurs and the M-D theory is not appropriate, requiring a separate treatment which will be presented below.
To summarize our discussion, the ultimate slip displacement can be defined as follows. When T ≤ μN and N ≥ 0, the contact zone center is fixed and the slip displacement develops at the sphere center, oriented tangentially to the contact plane. Its ultimate value δsu is reached in the limit state T = μN and is expressed uniquely by (7) and (8) in terms of h 0 or N . It does not depend on the slip history occurring before reaching the limit state.
For the various impact problems, the contact velocity vector is not horizontally pointed as it was assumed in M-D theory. In this case, the sliding regime can be also specified by a limit impact angle φ u , while the angle φ represents the angle of the sphere center velocity inclined to the normal of contact (Fig. 2a). Following the analysis of [15], the limit condition for the pure sliding regime (when c/a H = 0) for the spheres of different radii in contact can be expressed as follows For the case of the same material of the spheres, relying on (8), we can get [15], thus: In Fig. 2a, the plotted circle of radius R s = (R 1 + R 2 ) sin φ u indicates the sliding regime limit for all paths of the upper sphere, emanating from the point O 1 and not intersecting the circle and otherwise. Thus, a conical domain with its vertex at O 1 and tangential to the sphere of radius R s specifies the slip regime and the external domain specifies the sliding regime, bounded by the plane normal to OO 1 . The contact separation occurs for the paths emanating from O 1 in the exterior of the slip and sliding domains.
In Fig. 2b, a bilinear motion path of the upper sphere is represented. First, the contact zone is created by the overlap  (Fig. 2b). This case is typical for the sphere motion initiation from the static equilibrium configuration, where the primary slip regime always precedes the sliding regime up to the contact separation. The plotted circle of radius R s demonstrates the slip (inside the circle) and sliding (outside the circle) limits for the bilinear motion paths. These cases are considered in detail in the next section. Now, referring to the evolution of loading surfaces plotted in Fig. 1c, it can be demonstrated that, if the stress point A reaches the limit friction locus, the active loading surface F 1 = 0 moves toward the origin and F L = F 1 = F 2 = 0, and the adherence zone is erased, while the whole contact zone reaches the sliding state with the slip displacement distributed parallel to the tangential force T . This evolution occurs, when the contact zone is created first by normal load N and next the slip zone starts to progress due to the tangential load T . In this case, when T = T A < μN is applied as shown in Fig. 1b, then no sliding is generated, (i.e., surface F 1 = 0 lies within the domain bounded by the limit surface F L = 0, in Fig. 1c) and only the slip zones propagate from the external contact perimeter. For linear sliding trajectories, satisfying the condition φ > φ u , the contact zone is generated simultaneously with the distributed slip displacement over the sliding path.
Thus, the limit slip displacement δsu, acting on the sphere center parallel to the contact plane, evolves simultaneously with the sliding displacement s r (Fig. 1b), travels along the sliding path and vanishes at contact separation. Both slip and sliding displacement components are superposed: where s s is the slip component of displacement δsu, when the contact plane rotates during sliding. For the slip or sliding reversal, s s affects the shape of hysteresis loop [15,16].

The main parameters of contact geometry
Let us consider translation of the sphere of radius R 2 along the imposed linear path OO 2 ( Fig. 3a) with velocity v s (t 0 ). Suppose the bottom sphere to be fixed, and the upper sphere center to be attached to a moving tool imposing this path. The rotation of this sphere is supposed to be restricted by the moving tool and translational motion is being considered. The restriction in rotational motion is a typical case in the experimental testing of the spherical grains contact interaction [3,20,22]. For example, the grains were attached to stainless mounting pins of the apparatus with epoxy in the experimental work [22]. The tangential viscosity, artificially damping the oscillatory behavior of the sliding sphere during its slow down, is not considered. To this end, the application of the force-slip displacement law for the dynamical contact interaction has been demonstrated in the preceding work [16]. The contact condition is R − y 0 > 0 (where R = R 1 + R 2 , y 0 is the vertical coordinate of the sphere center), and it is supposed that the maximal overlap h 0 value is given; then, y 0 = R − h 0 = const. For the rightward motion, the contact engagement is at the point A, and separation is at the point D (Fig. 3).
In view of Eq. (11), the sliding regime develops instantaneously, when y 0 ≥ y lim = (2−ν)μ 4(1−ν) R = R tan φ u and Coulomb friction rule is valid. For high overlaps, the slip regime can develop from the contact engagement and the N − T path diagram (cf. Fig. 1) is initially located inside the sliding friction locus. Such case is discussed in Sect. 4, when the contact engagement state is generated by vertical force.
In Hertz theory, the force N is expressed via the total overlap h t by Eq. (7). The maximal individual h 1 and h 2 (Fig. 3a) could be assumed to be proportional to their compliances, so that and In the present work, the contact point is specified geometrically via intersection of sphere profiles (Fig. 3). In this case, by denoting R 2 /R 1 = ξ, we have (considering the second order geometric approximation): These two specifications of h 1 and h 2 are equivalent, when η = ξ, so then the ratios of radii and the stiffness moduli are the same. However, it turns out that there is no need to select a specific contact trajectory for each case of different radii spheres in contact, since the contact interaction depends on the varying overlap h t and the orientation of con-tact plane angle α t , and the evolution of the driving force can be expressed in terms of these two variables. The plotted contact path (Fig. 3a) should be regarded as the reference path, which, in turn, represents the actual path for the spheres of equal elastic moduli. Using formulae (13,14), any path can be reconstructed from the specified overlap and used for calculation of frictional dissipation referred to the contact path.
The contact path is constructed geometrically by tracking the center of the line obtained by intersection of two spheres. It represents by functions ds 1,2 = r 1,2 d α t , defined by the contact path radii r 1 (α t ) = const, r 2 (α t ) = const, for α t ∈ [± α u , ∓ α u ] (Fig. 3a). The upper signs are referred to the contact engagement at +α u and separation at −α u , for the leftward motion and lower signs are used for the rightward motion (cf., Fig. 3).
Referring to Fig. 3a, it can be seen that h r is composed of two different overlaps h 1,2 = R 1,2 − r 1,2 , and , and the contact center point is symmetrically located relative to a sphere-sphere intersection points. Note, in the limit case, if the upper rigid sphere slides over the bottom elastic sphere, the overlap develops on the bottom sphere only, and then: In the opposite case, when the bottom sphere is rigid, while the upper sphere is elastic, then Primarily, by assuming that the slip displacement is s s = 0, the sphere center coordinate and displacement can be expressed as s r (α t ) = y 0 tan α t and s r = s r − s 0 = y 0 tan α t − tan s g |α u | (where s 0 is the horizontal coordinate of the upper sphere contact engagement, s g |α u | specifies the contact separation angle, s g = sign (v s ) defines the direction of motion according to the axis x (Fig. 3), v s is the sphere center horizontal velocity). The h r corresponding to path coordinate s r can be specified as where cos , and α u = atan (s u /y 0 ). Moreover, the sign ± implies the symmetry conditions, i.e., s t ∈ [±s u , ∓s u ]. The radius of the contact zone is a r = ± R 2 1 − r 2 1 . Also, for the contact radius a H , calculated from Hertz formula, is smaller than the sphere-sphere intercept value a u , i.e., a u /a H → √ 2, for h 0 /R 1 → 0. The following set of equations allows for finding the sphere radii at the contact point C (Fig. 3), thus, Solution of this set yields the following formulae for the contact path radii where C * = R 2 1 − R 2 2 /y 0 . In (Fig. 3a), the path defines as:

Basic rate relationships for sliding motion
Let us first specify the rate equations by neglecting the slip displacement, s s = 0. Hence, in view of the above formulae, the overlap and displacement rates (or increments) can be expressed as follows (Fig. 3b): Thus, the sphere center velocity v s =ṡ t can be decomposed into the tangential to contact plane sliding velocity ṡ c and the normal to this plane velocityḣ r , expressing the overlap change during sliding ( Fig. 3b): where ṡ 1 and ṡ 2 are the velocities of contact point C motion along the curved contact paths at the bottom and upper spheres (Fig. 3a). Using Eq. (7), the normal N (h t ) and tangential ±μN (h t ) contact forces (Fig. 4a) can be expressed as Hence, the equilibrium conditions in the sliding regime are defined in terms of the forces oriented in the normal and parallel directions to the translation path ( Fig. 4a): where ϕ = arctan μ is the angle of friction. Note, that the vertical force N * is equilibrated by the vertical reaction force As the sliding sphere is constrained against rotation, the reaction moments M 1 = F R 1 and M 2 = −F R 2 are equilibrated by the constraints. It is assumed that the condition | α t | ≤ φ is satisfied along the entire path. Accounting for the angle of friction, the active loading surface, expressed in terms of forces N R , F and configuration parameter α t , can be now defined as follows Figure 5 presents the diagrams of N * and T * forces specified by Eqs. (23) and (24) as the functions of the sphere center coordinate s t = s r . The overlap h r is predicted by Eq. (16). The values of h 0 = 0.01R 2 , R 1 = 1 m were assumed. The effective Young's modulus was set E * = 1Pa, providing It allows for a simple rescaling of the given results, e.g., Fig. 5, T * and N * diagrams are not symmetrical relative to α t = 0, since the contact plane rotates during sliding process. In particular, for s t ∈ [−s u , 0], the horizontal component of normal force N is pointed in the same direction as the friction force, while, for s t ∈ [0, s u ], it is oriented oppositely to this force (cf., Fig. 4) and due to this asymmetry develops towards the direction of motion.
The contact orientations corresponding to the maximal values of forces T * max and N * max can be determined from Eqs. (22)- (24), requiring d T * /d α t = 0 and d N * /d α t = 0, and the following equations are derived: where α T and α N are the contact angles corresponding to extremal values of forces T * and N * . In Eqs. (26,27), the angles α T and α N depend only on the friction angle φ and initial position cos α u = y 0 /R, and are independent of the contact stiffness. This condition is very important for the experimental testing, i.e., when the values of relative position y 0 /R and angles α T (or α N ) are measured, the experimental friction coefficient can be established from the formulae (26,27). Moreover, T * max and N * max are also important for selecting the loading ring capacity in developing the inter-particle friction apparatuses. The results of numerical solution of Eqs. (26,27) are shown in Fig. 4b, for the rightward motion of upper sphere of different radii and initial overlap of h 0 = 0.1 R 2 (R 1 = 1 m).
As can be seen in Fig. 4b the maximal values of resultant horizontal force correspond to orientation angles α T < 0 and depart from the vertical line of contact symmetry with decrease in the coefficient of friction and increase in the spheres radii ratio. The angles α N > 0 are also mainly shifted from the symmetry line with the same trend. The values of T * max and N * max are also depicted in Fig. 5. Meanwhile, Fig. 6 presents the force path F − N R for the case of monotonic rightward motion along the linear trajectory and evolution of the limit surface F L = 0 represented by two Coulomb lines rotating about the origin O with the symmetry axis oriented at the angle α t to the N R -axis.
In Fig. 6, during the sphere sliding process the limit locus rotates as α t varies from −α u to α u . The normal load effect surface F 0 = 0 translates with the loading point C.
In the sliding process, the path reaches the maximal value of the driving force F at point P F , for α T ≤ 0 (Fig. 6a,  Fig. 4b). Next, the extremal point P N is reached at α N ≥ 0 at which N R attains its maximum (Fig. 4b). The graphs in Fig. 6b-d present the evolution of force path until the contact separation at point O. The tangent moduli values at O are different for activation and separation states, namely, Fig. 6d), for the contact activation and separation points, respectively (for the rightward motion). When the reverse translation is imposed from the state represented by the loading point C, then the unloading path follows the curve CC 1 in the slip mode and finally reaches the sliding regime at point C r . In this case, the active loading surface F L = 0 and the normal load surface F 0 = 0 evolve and specify the reverse slip displacement. The reverse slip and sliding modes were discussed in detail in [15].
Thus, by combining the diagrams of Figs. 1 and 6, the combined slip and sliding memory rules are now represented by the N R −F limit and loading surfaces. The memory effects in slip regime and the configurational effects in sliding regime are coupled and defined using the active loading surfaces which evolve in the space of contact forces.
In work [16], the main relationships for dynamic sphere motion have been derived for parameters related to two equal spheres. In these relationships, by applying substitutions k n = k * n and R 5/2 0 = R 5/2 / 4 √ 2 , the derived formulae can be adopted in the case of different radii spheres in contact. The coefficient of restitution in tangential direction is the basic parameter for defining the frictional energy dissipation. It can be predicted by the following formula: where m is the sphere and moving device mass, K (k) and E (k) are the first and second kind complete elliptical integrals of the elliptic modulus k = 1 − cos 2 α u /(1 + cos α u ). For predicting the time of contact, t c , etc., the formulae, given in [16], can also be adopted.

Basic rate relations for the combined slip and sliding motion
Let us consider the imposed translation of the sphere with account for both the slip and sliding velocities (Fig. 3a, c). The slip displacement is referred to the sphere center motion paralelly to the contact plane, but at a fixed contact zone center position at the point C. Meanwhile, the sliding is manifested by the contact zone motion relative to both contacting spheres, and by the variation of the overlap in the normal direction. Thus, the translation velocity at the sphere center can be now decomposed as v s = ṡ t = ṡ s + ṡ r (Fig. 3c).
In this case, the sliding regime develops simultaneously with the upper sphere contact engagement, but the ultimate slip displacement evolves along the whole translation path, since the overlap h t is subjected by the change in contact plane size. For the left-or rightward motions it is expressed by Eq. (9), as follows: where χ = s g (2 + ν) μ/4 andḣ t =ḣ s +ḣ r are the normal contact velocities, resulted from the slip and sliding motions due to change in the overlaps h s and h r , respectively (Fig. 3c). Using these formulae, the slip displacement can be defined as ṡ s =u 0 /cos α t = χ ḣ s +ḣ r /cos α t , wherė h s = −u 0 tan α t (cf. Fig. 3c). Using Eqs. (19) and (29), we can geṫ then, the sphere center translation and total overlap velocities are expressed as follows where the approximate forms are also provided. Now, the exact rate Eqs. (31, 32) are integrated to provide the relationships, specifying the effect of slip displacement on the total displacement s t and the total overlap h t , thus: where α 0 = −s g |α u | is the initial angle of contact engagement.
At contact engagement, when s r = 0, we have s t = 0 and h t = 0. At contact separation, for s r = 2y 0 tan s g |α u | and α t = s g |α u | = s g atan (|s u |/y 0 ), we have s t > s r , h t > h r and h t = 0. In this case, the contact separation angle α * u ≥ α u is found numerically using the condition h t α t = α * u = 0 in Eq. (33b). The approximate forms of the above formulae are as follows: The coupled effect of slip induces the overlap growth and affects the sliding displacement values, but these interactive effects are rather small. As the slip displacement is related to the sphere center motion, then it results from the sphere deformation. The obtained prediction by Eq. (33a, b) could be further rigorously verified using a FEM-based elastic 3D solution of the sliding contact problem. Eq. (33a, b), we obtain: while the approximate forms of these equations are obtained from Eqs. (33c, d): From Eq. (34c), for α t = α u , the coupled sliding displacement s t (α u ) is the same as the pure sliding displacement s r (α u ), while the contact separation position is specified from h t α t = α * u = 0 in Eq. (33b). Thus, in contact interaction analysis, either exact or approximate formulae can be used. It should also be noted that neglecting the effect of slip induced overlap component h s and assuming h t = h r , we directly obtain the condition ṡ t = ṡ r (1 − χ tan α t ). Then slip displacement affects only the tangential sliding velocity, while the contact separation is at the angle α t = α u . In summary, formulae (33, 34) express analytically a coupling of slip and sliding modes during sphere translation along the linear path.
In view of our analysis the slip displacement effect can be treated by the following schemes: (a) exact prediction (Eq. 33a, b); it means a solution of the case ofḣ t =ḣ r +ḣ s accounting for the term 1/(1 + χ tan α t ); (b) approximate prediction; it means that the mathematical approximation 1/(1 + χ tan α t ) ≈ 1 − χ tan α t is applied and formulae (33c, d) are derived; (c) simplified prediction; it means that the assumptionḣ s = 0,ḣ t =ḣ r is applied, yielding directly the slip effect in the form of 1 − χ tan α t with formulae (16) and (33c).
In Fig. 7a, the effect ofḣ t =ḣ r +ḣ s (regular lines) againsṫ h t =ḣ r (dashed lines) is shown on the total overlaps h t and h r . The evolution of slip displacement s s versus sliding  (Fig. 7b), its rate ṡ s (Fig. 7d) and the overlap rateḣ s due to traveling slip (Fig. 7c) are also given. From Fig. 7 a, it can be seen that the account for slip component rate,ḣ s = 0 by Eq. (33b) generates the additional portion of the overlap which, in turn, results in the increased values of contact separation angle, i.e., α * u > α u . It occurs due to the sphere center distortion towards the direction of motion, accumulated from the changes of overlap (cf., h s ( s r ), in Fig. 7b), since the positive rate ofḣ s (α t ) is held over the whole sliding path (Fig. 7c). It results in supporting the deformed sphere before contact separation. Moreover, the slip displacement rate ṡ s is distributed asymmetrically, for α t ∈ [−α u , 0] and α t ∈ 0, α * u , which finally results in the negative value of s s α * u (Fig. 7d).

The combined slip and sliding from the initially activated contact zone
Consider the most frequent case in the experimental contact testing, when the sphere motion proceeds from the activated contact at α δ t0 = 0 (Fig. 2b). The transverse slip and sliding from the formed contact zone was discussed by Fouvry et al. [23] and Dini et al. [25] in relation to the analysis of fretting wear.
The initially applied load N 0 activates the contact zone of radius a δ 0 and overlap h δ 0 (Fig. 8a). The subsequent slip and sliding of the upper sphere relative to a fixed bottom sphere is then induced along parallel or inclined paths relative to the contact zone. After activation of contact zone, the coordinate system (x 1 , y 1 ) is introduced with the x 1 -axis following the sliding path (Fig. 8a). For δ > 0, the trajectory path is inclined clockwise to the transverse direction and the initial slip process is related to growth of contact zone. For δ < 0 the slip is combined with the reduction of contact zone size. Note that the leftward motion for δ > 0 along the path O 1 O 3 generates the same response as for the rightward motion for δ < 0 (Fig. 8a).

Basic relations for the combined slip and sliding motion
In general, the resultant forces acting onto x 1 , y 1 axes either for the slip or sliding modes are given as: where N h δ t = k * n h δ t 3/2 .
Slip mode After contact activation and next imposed sphere translation the slip displacement δs of the sphere center is generated along the initial contact plane remaining horizontally fixed at the point A (Fig. 8a). The slip displacement proceeds along the transverse direction and induces the sphere center displacement and the penetration depth normal to the contact plane s δ s = s δ t = δs/cos δ and h δ s = δs tan δ, respectively. Thus, the overlap is composed of the initial value and the additional component due to sphere motion, i.e., h δ t = h δ 0 + h δ s , and normal contact force equals: The slip regime terminates, when δs = δsu = χ h δ t = χ N su /k * Thus, at the slip-sliding transition state the following relations occur: Using M-D theory, the tangential contact force can be expressed by following formula: Finally, for the selected slip displacement range δs ∈ [0, δsu], the forces T (δ s ) and N h δ t are calculated from Eqs. (40) and (37). Setting α δ t = 0 in Eqs. (35, 36), the evolution of resultant forces along the slip path s δ s = δs/cos δ is determined (cf., Fig. 9).
Sliding mode After reaching the slip displacements s su and δs = δsu. (Fig. 8c), the sliding of the sphere develops at any α δ t . The subsequent displacement is s δ t = s δ s + s δ r and the friction force is T (δ s ) = μN h δ t . It yields now that the problem becomes statically determinate. So, putting T (δ s ) into equilibrium equations (35, 36), the driving force equals where N * 1 α δ t = N h δ t cos α δ t − δ − s g φ /cos φ is the resultant force along the y 1 -axis representing the reaction of the moving device. Hence, the limit surface (41) evolves with the variation of the contact plane inclination angle α δ t (cf, Fig. 9c).
During sphere sliding (Fig. 8a), the overlap and sliding displacement can be expressed as follows: Referring to velocity decomposition in the sliding regime (Fig. 8c), the sphere center translation velocity v s =ṡ δ t = s δ s +ṡ δ r involves the slipṡ δ s =δ su /cos α δ t − δ and slidinġ s δ r = y 0 cos δα δ t /cos 2 α δ t − δ velocity components. In this case, by depending on the overlap values, the limit slip displacement, δsu = χ h δ t , travels to the contact separation with the velocity ofδ su = χḣ δ t (Fig. 8c). Meanwhile, the normal contact velocity is also decomposed by the slip and sliding components,ḣ δ t =ḣ δ s +ḣ δ r , represented as the overlap rates induced during slip and sliding modes (Fig. 8c). Sincė h δ s = −δ su tan α δ t − δ , combining the these three equations we haveδ su = χ ḣδ r −δ su tan α δ t − δ , which finally yields the following rate equation: Referring to Fig. 8c, we haveḣ δ r = −ṡ δ r sin α δ t − δ , and we can find that: Finally, the sphere center translation velocity and overlap rate are expressed by following formulae: wherė Next, the integration of these equations couples the slip and sliding modes, thus that finally yields the following formulae: and for the approximate model, we have In Eq. (50a, b), by setting h δ su = 0 and δ = 0, we arrive at formulae (33b) and (33d), for α 0 = 0. It yields the correctness of the obtained equations. In the coupled slip and sliding mode, the contact separation angle α * u is found numerically by setting h t α t = α * u = 0 in in Eq. (50a, b). In Fig. 9a, b, the resultant force-path diagrams are plotted by assuming h 0 = 0.01R 1 . The same motion direction is supposed, but the positive and negative trajectory angles are applied.
From Fig. 9a, b, it can be seen that the resultant force-path diagrams indicate the different deformation processes after contact activation. In both cases, first the slip regime develops, while the sliding mode evolves out of the slip mode and proceeds up to the sphere separation. The effect of slip displacement evolution during the sliding mode is shown by red lines (Fig. 9a). This effect generates α * u > α u , due to sphere center distortion towards the direction of motion, accumulated from the changes of overlap during the slip mode. For the softening behavior this effect is very small (Fig. 9a) and invisible. However, for the class of the problems involving contact activation and separation or slip and sliding reversals, the slip effect is important and cannot be omitted. For δ > 0, the resultant force levels are up to 4 times higher than for δ < 0, where the maximal normal force is generated after contact activation. The evolution of limit surface F − N * 1 in the slip and sliding domains is presented in Fig. 9c. The slip regime initiation is denoted by bold points in the diagram. Then the slip mode forces evolve almost instantaneously and, when the sphere starts to slide, the limit friction curves suddenly change the directions and next develop along quasi-elliptical loops, reflecting the variation of contact plane inclination angle α δ t for the case of leftward motion. Increase in the overlap values generates the elongations of the loops, and the contact separation occurs when F = N * 1 = 0. The effect of slip displacement evolving during sphere sliding is given by the dashed lines. It increases, when the contacting sphere radii ratio grows. For the leftward motion, this effect is practically invisible due to softening behavior.
The level of driving force, which can be achieved by the adjustment of the trajectory angle imposed, is very important in the experimental testing of contact parameters of two spherical grains. The driving force evolution in the slip-sliding modes for different angles δ is schematically represented in Fig. 10.
It can be seen (Fig. 10) that different maximal values of driving forces can be reached at point positions P m under different slip-sliding displacement limits. The slip to sliding transition point P s can be reached via positive or negative values of angle δ. For δ < 0, the driving force, slip displacement limit and whole path length exhibit the small values due to softening behavior (cf., Fig. 9b). In this case, the slip to sliding transition value of the driving force P s can occur on the descending branch of the force-displacement curve. Meanwhile, for the cases of δ = 0, the maximal value of driving force is attained at the termination of slip regime (P m = P s ). Hence, when the driving force is monotonically increased during experimental cyclic loading near the force peak value P m , the small contact surface imperfection of the grain can lead a spurious jump from the slip displacement to excessive sphere sliding with the subsequent effect of ratcheting, cf., [22]. The limit angles at which the driving force takes peak values are given in Fig. 4b. Moreover, it can be seen (Fig. 9d) that s δ su << s δ ru are valid; however, for δ < 0, the s δ su is smaller than s δ ru by about 3.3 times and cannot be neglected, especially in the cyclic loading analysis.

Sphere contact interaction in a mixed control process
Let us consider now a two-point contact interaction of an elastic sphere of radius R 2 with a bed of fixed regularly packed spheres of radius R 1 resting on a rigid substrate (Fig. 11a). For simplicity, only the plane problem is considered by assuming the spheres to enter contact and slide along their diameters. In this case, the upper sphere out-of-plane motion is prevented by the pooling device, i.e., the sphere is attached to the actuator, oriented to produce a plane shearing motion [22]. The actuator can be equipped with a load-cell capable to induce a specified load level at a specified rate [23]. Consider first the contact engagement process by placing the upper sphere center at the point O 0 (Fig. 11a) to make the sphere diameter tangential to two bottom spheres and next applying the vertical (gravity) load Q to generate the equilibrated contact configuration. Then the slip or slip-sliding regimes can develop. Next, the sphere sliding process is activated by a pulling device driven by the horizontal force T 0 under fixed load Q. In this case, the upper sphere horizontal displacement is controlled and the evolution of T 0 depends on the contact stiffness, coefficient of friction, ratio of contacting spheres radii and contact separation mode.
Generally, there are four stages of contact interaction: (a) symmetric dual contact engagement of the upper sphere with two bottom spheres under action of force Q (Fig. 11a); (b) interactive contact separation induced by the driving force T 0 (starting from the point C 1 ) and right contact sliding mode activation (at the point C 2 ) (Fig. 14a); (c) sliding motion of the sphere along the right bottom sphere in a single contact mode; (d) non-symmetric interactive dual contact separation-(a) (b) Fig. 11 Contact interaction process: a symmetric contact engagement by the force Q; b interactive contact separation-engagement by the horizontal force T 0 engagement of the sphere with the subsequent bottom sphere from the bed (Fig. 11b).

Symmetric contact interaction due to application of force Q
For Q = 0, the initial angle α 0 u specifies the initial contact configuration for y 0 u = R cos α 0 u . It increases to value α * u as the load Q is introduced and the sphere center point O 0 translates vertically along the symmetry line O 0 O by the displacement s t (Fig. 11a). The contact response is the same as that discussed in the Sect. 3.1, for the linear trajectory of sphere center motion. Thus, for the initial and loaded configurations, respectively, we have: where h 0 = h 1,2 is the overlap between the upper and bottom spheres caused by load Q (Fig. 11a).

Contact slip mode
The slip regime arises first after application of the vertical force Q. The planes of contact are at the initial angular configuration α 0 u . The slip condition (10) requires tan α 0 u ≤ tan φ u , so using Eq. (51a), we have: where φ u = arctan χ .
From this condition, it follows that large upper spheres, for which R 2 ≥ ξ c R 1 , after loading by the force Q engage contact in a slip mode. For μ = 0.5 and ν = 0.3, the vertical slip response can be obtained only when R 2 ≥ 2.45R 1 . The smaller upper spheres of size R 2 < ξ c R 1 enter contact in the sliding mode.
The applied load Q deforms the spheres and equal overlaps h 1 and h 2 with the equal reaction forces N 1 and N 2 on the left and right bottom spheres are generated. It requires the following equilibrium condition: In the slip response, the sphere center vertical displacement equals s t = s s generating the normal and tangential components of h 0 = s s cos α 0 u and δs = s s sin α 0 u , respectively (Fig. 11a). Then, the contact tractions equal Note, the positive values of Q and α 0 u should used, while s s should be treated as negative following the orientation of y-axis (cf., Fig. 11a). For a frictionless contact, it yields, s s = Q 2/3 2k * n −2/3 cos −5/3 α 0 u .

Contact sliding response
For ξ < ξ c , the coupled slip and sliding displacements develop simultaneously. So, the upper sphere vertical displacement is now composed of s t = s s + s r (Fig. 11a) and the obtained rate Eqs. (29)-(32) are valid with the substitution of α t = π/2 − α * u . In this case, in Fig. 3a denoted distance y 0 now is equal to R 1 (cf., Fig. 11a). Hence, the overlap and the sliding displacement are: The above rate relations remain valid and in view of Eqs. (31) and (32) we have: Integrating these approximate forms and substituting Eq. (55), one obtains: The overlap h 1,2 can also be expressed in terms of s r . So, using Eqs. (55) and (57b) it can be written Finally, the displacement s t can be calculated by Eq. (57a) by substituting angle of loaded configuration, α * u , which is obtainable from the solution of equilibrium equation: Note, the solution for α * u can be achieved numerically (e.g., Newton-step iterative technique can be embraced). The slip displacement components are calculated as, δs = δsu = χ h 1,2 α * u and s s = s t − s r . The vertical displacements and overlap are plotted in Fig. 12. The value of μ = 0.5 is used, while the load Q is applied gradually. The value of E * = 30000 Pa is selected as for a soft material.
As it is seen in Fig. 12, with the increase in the upper sphere size the stiffer response of granular bed is generated due to increase in the contact stiffness k * n . It can be noted that, in the slip mode (ξ > ξ c ), the deformation of granular bed is highly reduced. The obtained final values of the overlap h 1,2 , the tangential displacement δs (Fig. 12a) are applied as the initial values, for the subsequent simulation, when the horizontal force T 0 is applied. They are used for predicting the initial contact tractions N 1,2 and T 1,2 . Table 1 also presents the predicted contact engagement parameters of different size polypropylene (PP) spheres with regularly packed bed. The vertical force applied for i-th upper sphere is assumed to be equal Q i = Q 1 (R i /R 1 ) 3 (Q 1 = 8.71 N, R 1 = 0.03 m are the force and radius of a smallest particle considered). Other material parameters are: E 1,2 = 1.82 GPa, ν 1,2 = 0.3, μ = 0.3 yielding E * = 1 GPa. The radius of bottom sphere is held constant and supposed to be equal R 1 = 0.1 m. Table 1 shows that the increase of the upper sphere radius above R 2 = 4.88R 1 generates the contact engagement within the slip response. Otherwise, the combined slip and sliding displacement occurs. Generally, when ξ < 1, the slip effect is small and s s << s r . It becomes more pronounced, when ξ → ξ c . Also, it is worth noting that for R 2 > 4.88R 1 , the overlap and vertical displacement values are approximately equal, but tendency of δs << h 1,2 is hold independently of the value ξ.

Contact interaction induced by application of force T 0
Consider now the contact interaction under action of the driving force T 0 . To this end, the contact state evolution in the force plane N − T is analyzed in Fig. 13. In Fig. 13a, for the left contact zone (cf., Fig. 11a, point C 1 ), at the initial configuration, the Coulomb cone is inclined at the angle α 0 u > 0 to x-axis (Fig. 13a). After application of vertical force Q, either slip (ξ ≥ ξ c ) or sliding (ξ < ξ c ) regimes can develop. In the sliding regime, the Coulomb cone rotates due to variation of contact plane orientation by the angle increment d α > 0, and the point A moves along the limit friction locus due to rotation generating the curve OA (Fig. 13a). For the slip regime, there is no variation of the contact plane orientation and force T 1 (N 1 ) response, illustrated by the curve OA 1 , evolves inside this locus (Fig. 13a). When the driving force T 0 is gradually applied, the reverse slip path AB is generated with the subse-quent sliding path BO terminated by the contact separation at the left contact zone.
In Fig. 13b, the contact force evolution is illustrated for the right contact zone (cf., Fig. 11b, point C 2 ). Now, at the initial configuration, the Coulomb cone is oriented by the angle α 0 u < 0 (Fig. 13b), and load Q induces the force T 2 (N 2 ) path OA generated after rotation of contact plane by the increment d α > 0 to the angle α * u . The slip response path OA 1 can be also generated for larger upper sphere. The subsequent application of T 0 induces the reversed slip displacement with the gradual increasing in the magnitude of the normal contact force N 2 following the path AC (cf., Fig. 13b). Subsequently, it reaches the reverse sliding mode at the point C and the single contact sliding of the upper sphere develops along the right bottom sphere. Finally, it terminates, when the sphere touches the next bottom sphere. In this case, the contact planes orientations on the left and right contact zones evolve asymmetrically (Fig. 11b). During new contact activa-(a) ( c ) (b) (d) Fig. 13 Contact force paths: a, b contact engagement-separation after application of force Q and subsequent gradual loading by driving force T 0 ; c, d contact separation-engagement with the new sphere after single contact sliding mode tion the normal contact force N 1 gradually decreases to zero in the sliding mode OB 1 (Fig. 13c), while the driving force increases along the slip path OA and turns out into the sliding path at the point A (Fig. 13d). The incremental numerical analysis supports the presented theoretical diagrams. Referring to Fig. 14a, consider the equilibrium conditions required for loading by the driving force T 0 in the final contact state generated after loading by the force Q. These conditions can be written as follows: while the incremental conditions are expressed in the form: where α lt u , α rt u denote the contact plane orientation angles at the left and right contact zones (Fig. 14a) = α * u , for ξ < ξ c , and α lt u = α rt u = α 0 u , for ξ ≥ ξ c . In the primary state after application of load T 0 , the slip regime develops at both contact zones. The slip displacement follows from M-D theory, while Masing's rule is applied for the reverse slip, thus: where A = χ k * n −2/3 /μ = 3/(32G * ) (2K * R * ) −1/3 , K * = 3/(8E * ), δ A s and T A are the slip displacement and force values at force reversal.
In Eq. (63), the contact forces N 1 , T 1 and N 2 , T 2 , acting on the left and right contact zones (cf., Fig. 14a At the initial state of loading by force T 0 , the slip modes develop at both contact zones (Fig. 13a, b) with d α lt u = d α rt u = 0 and α lt u = α rt u = |α u | (where either α u = α 0 u or α u = α * u are dependent on the contact mode generated by the force Q). So, the solution of Eq. (62) refers to two unknowns d s and β, while the compatibility equations between the overlap and displacement increments are as follows: while the sphere center trajectory increments are: d x = d s sin (|β + |α u ||), d y = d s cos (|β + |α u ||).
The sliding mode at the left contact zone (cf. lines BO or B 1 O in Fig. 13a), accompanied by the slip regime at the right zone (cf. curves AC or A 1 C 1 , in Fig. 13b), requires the solution of Eq. (62) for the unknowns d s, β and d α lt u . However, the effect of contact configuration change d α lt u on the contact forces evolution is quite small up to contact separation and can be neglected and the same compatibility equations are adopted. In this case, the effect of slip displacement on the sliding regime is included by stating Eq. (64) in the form d s cos β = −d h 1 /(1 + χ tan |α u |), d s cos (β + 2 |α u |) = −d h 2 /(1 + χ tan α u ). Finally, it results in the following solution of Eq. (62) for the trajectory angle β and displacement increment d s, for the right and leftward motion of the upper sphere: where ω and B (β) are given explicitly in "Appendix A", for the states at the left and right contact zones. The contact activation by after application of force T 0 terminates, when N 1 = 0, and sphere starts sliding over the bottom sphere in a single contact mode. In this case, the last incrementally calculated load value T 0 is the limit friction force T 0 μ . Its subsequent values continuously decrease up to new contact interaction with the bottom sphere. During the single contact sliding the solution is continued by the imposing increments of horizontal sphere center displacements, following the equations given in [15].
After the single contact sliding mode, the upper sphere comes into new contact with the bottom sphere (Fig. 11b). However, during this contact activation, the initial overlaps h 1 and h 2 are no longer equal and develop non-symmetrically. Moreover, the reaction forces N 1 , N 2 are inclined at different angles α lt u , α rt u (where α lt u + α rt u = 2 |α u |) (Fig. 11b). It requires new geometrical compatibility conditions to be imposed. For simplified analysis, it is supposed again that configurational change at the left contact zone during new contact activation can be non-significant, d α lt u = d α rt u = 0. Hence, from Fig. 11b, two angle relations can be written, namely, (R − h 1 ) cos α lt u = R cos α rt u and R sin α rt u + (R − h 1 ) sin α lt u = 2R 1 . A solution of these equations yields the following formulae: During simulation, as the upper sphere touches point C 2− , we replace the coordinate system into the center of new bottom sphere (Fig. 11b). Therefore, the positive angle α lt u found at the sliding regime termination takes a negative value and serves as the initial value for new contact activation analysis. In the compatibility conditions (64, 65), the angle 2 |α u | is replaced by the sum of angles α lt u + α rt u . The formulae for prediction of the sphere center trajectory angle and path increment are calculated from Eq. (66), using parameters ω and B (β), specified in "Appendix A". Also, the trajectory increments are calculated as d During new contact activation, a mixed sliding/slip mode generally occurs. The unloading sliding regime lasts unchanged at the left contact zone without any switch to slip mode (Fig. 11c, curve AB), while at the right contact zone (point C 2 ) the slip mode progresses in the loading state. Finally, the new contact activation process terminates, when N 1 = 0 and the contact separation from the left bottom sphere occurs. All relationships can be applied in the same manner.

Computational implementation
Let us discuss computational implementation of the above equations. Firstly, the upper sphere is placed in the contact with two bottom spheres at the position defined by angle α 0 u . Next, the vertical force Q is applied and the initial angle α 0 u or α * u is determined either for slip or sliding regimes. Next, the sign of this angle is specified in terms of the coordinate system for the left bottom sphere. The value of this angle is kept unchanged during the contact activation analysis in the slip regime.
Then, the driving force is applied and increased gradually using small increment value d T 0 > 0. The trajectory angle β and slip displacement d s are determined from Eq. (66) along with checking the appropriate contact mode conditions. After prediction d s, the overlap and tangential displacement increments are specified from Eqs. (64, 65). The overlaps, h 1,2t = h 1,2 + d h 1,2 , forces N 1,2t = N 1,2 + d N 1,2 , T 1,2t = T 1,2 + d T 1,2 and sphere center path s t = s u + d s are determined for any time t of simulation.
After termination of contact activation, when N 1 = 0, a single contact sliding of the upper sphere proceeds. The last value of the applied horizontal force and the last coordinate s 0 t value, determined in contact activation analysis, are treated as the initial values for the single contact sliding mode. Then, the sliding motion is analyzed for the increasing horizontal component of the sphere displacement, at which, the driving force is gradually reduced (cf., Fig. 15), next finding the angle increment dα t , overlap and path length s t = s 0 t + d s [15].
During this mode, the determined absolute value of α t is compared with the limit value of angle α lt u from Eq. (67) in order to control find new contact activation at the point C 2 (Fig. 11b). For |α t | ≥ α lt u , the upper sphere comes to new contact activation process with the next bottom sphere. Then, the last value of angle α t is selected to be initial for the new contact activation analysis. The sign of this angle is changed to be appropriate with new coordinate system (Fig. 11b). The last applied horizontal force and trajectory coordinate values, determined in sliding mode, are then also selected to be the initial values for this stage of analysis. Next, the horizontal force increment d T 0 > 0 is imposed again. The next calculation steps are similar and rely on prediction of β and d s.

Results and discussion
In Fig. 15, the slip-sliding and sliding-slip mode results are presented for the rightward motion of the upper sphere for ξ = R 2 /R 1 = 0.4 and 0.8, μ = 0.5, Q = 100 N and E * = 30000 Pa (R 1 = 1 m). Thus, after application of the force Q, the upper sphere deforms moving vertically down (Fig. 15c) under the sliding mode for both cases. Due to the higher contact stiffness, the generated overlap is smaller for the case of R 2 /R 1 = 0.8 (Fig. 15a). The positive values of tangential forces T 1,2 and displacement δ1,2 are also generated after application of force Q (cf., Fig. 15b). When the driving force T 0 is gradually applied, the contact activation in the slip mode develops at both contact zones. The sphere center gradually shifts to the right, leading to formation of the trajectory path with the rising and descending portions (Fig. 15c). The deeper rising-descending effect is seen for the case of smaller sphere (R 2 = 0.4m) having the lower value of contact stiffness.
The applied force T 0 gradually reduces the magnitudes of tangential contact forces and eventually changes their initial directions (cf., Fig. 15b). It occurs gradually. So, the unloading regimes proceed along with the increase in the slip displacements δ1,2 (cf., Fig. 15b). In this case, the normal contact force for the left sphere N 1 monotonically decreases, while the force N 2 increases for the right bottom sphere (Fig. 15a). Finally, when h 1 decreases to zero the contact separation at the left bottom sphere occurs. In this case, the tangential force T 1 also tends to zero along the softening curve with a non-zero value of the final displacement δ1 (Fig. 15b). Meanwhile, the tangential force T 2 and displacement δ2 monotonically increase. When the upper sphere separates from the left bottom sphere, it moves in a single  Fig. 16a) over the right bottom sphere. However, as can be seen in Fig. 16a, for the smaller spheres (R 2 /R 1 = 0.4), the single contact mode can proceed in the slip regime, and the sliding regime develops thereafter. In Fig. 15d, the driving force evolution along the sliding path is demonstrated. First, the value of force T 0 monotonically increases due to contact activation. At the maxima of T 0 the slip mode terminates and sliding regime develops at a single contact mode and then the limit friction force gradually decreases along the sliding path. Next, the sphere comes into contact with another bottom sphere and new contact activation begins. Such contact response was experimentally analyzed by Lukaszuk et al. [16]. As it can be seen in Fig. 15d, the contact engagement occurs at the negative value of driving force, in order to maintain the force equilibrium. This means that the upper sphere should be even supported by reversing T 0 direction until it reaches the new contact with the bottom sphere.
The local N-T force evolution during initial contact activation after application of driving force is shown in Fig. 16a. It is seen that the deformational behavior evolves within the slip mode, and contact forces develop inside the limit friction locus (Fig. 16a). In particular, the normal force N 1 decreases non-linearly with the decreasing tangential force T 1 tending to the limit locus. When the force path T 1 (N 1 ) touches the limit friction locus, the sliding mode evolves on the bottom left sphere and the unloading process continues within this mode until the contact separation reaching N 1 = 0 (cf., Fig. 16a). Meanwhile, the contact force N 2 monotonically grows up with the decrease in the values of tangential force T 2 .
In Fig. 16b, new non-symmetrical contact activation is illustrated. As it is seen, the sliding regime at the single contact interaction persists during the new contact activation. The normal contact force N 1 (regular lines) on the left contact zone decreases to zero along the force path remaining on the limit friction locus. At the new right contact zone the slip mode first develops. The force path T 2 (N 2 ) is located in the slip domain and evolves nonlinearly with progressing contact activation, finally reaching the limit friction locus. The sliding mode at the right contact zone is accompanied by full contact separation at the left zone.
In Fig. 16b, the path T 2 (N 2 ) starts at small value of N 2 , since the small overlap value (0.005-0.01R 2 ) is assumed at the instant when the upper sphere touches new bottom sphere at point C 2 (cf., Fig. 11b).

The practical aspects of application
Any direct measurement of the parameters of the spheresphere contact interaction within a granular material is not possible. The available data on contact friction, elasticity, viscosity and wear parameters, resulting from two spheres contact testing, are limited. As a rule, the developed interparticle friction test rigs (cf., [10,20,21]) are mainly based on the force controlled slip and sliding process. First, the contact zone is formed by a normal force, next a transverse driving force induces slip regime with its transition to the sliding regime. Such transition is required in order to measure the magnitude of the coefficient of friction. The measurement is complex because it is difficult to perform contact slip activation due to surface shape imperfections leading the spurious oscillations of the magnitude of the driving force [20]. More-over, the force oscillations tracked during experiment may lead near the transition point to the uncontrolled switch from slip to sliding mode due to a configuration-softening-based contact response [21].
On the contrary, for an apparatus based on the displacement control technique described here, these difficulties could be avoided. In particular, a simple application of the moving device generating linear trajectory of the upper sphere center does not require the contact zone activation by the normal force. It generates the sliding mode instantaneously with no preceding slip regime. So, in prediction of the coefficient of friction between two spherical grains using Eq. (26), the measurement of angle α T , at which the driving force magnitude reaches the maximal value, is only required. It should be emphasized, that this prediction is independent of the contact stiffness modulus. Next, by application of Eq. (28), the same test allows for prediction of the coefficient of restitution in tangential direction. This coefficient is always required in the DEM-based simulations of granular flow and suffers from lack of a proper experimental measurement.
The mixed (Q − v 0 ) control of finite sliding of the sphere along the bed of fixed spheres assures continuity in the evolution of contact parameters despite the complex interaction at the transition of sliding mode from one to the other grain and softening response during the single-contact sliding stage. The dual contact response model (or multiple contact in 3D case), discussed in Sect. 4, can be applied in identification of elastic, viscous and frictional contact moduli at both macro and nano-scales.
Moreover, the reference can be made to wear modeling with the account for the asperity interaction at the contact interface in the relative (monotonic or reciprocal) sliding of two bodies. The other reference is related to silo discharge flow of a granular material, where the moving particles in a central funnel zone slide down towards the outlet along the particles fixed in a stagnant zone.

Conclusions
For the sliding or combined slip-sliding regimes, a well known Mindlin and Deresiewicz theory cannot be applied and requires extension. In the present analysis, the combined slip and sliding rules for two spheres of different radii in contact have been proposed for simulation of the displacement and mixed controlled processes. The analytical relations of contact traction evolution along the sliding paths have been derived in the static analysis for the sphere motion along the linear and combined trajectory paths. The account of the memory effects in the slip regime and of the configurational effects in the sliding regime was also presented in terms of the active loading surface and memory surfaces in the space of contact forces.
In particular, the analytical coupling of slip and sliding modes during sphere displacement controlled processes is of sufficient complexity. Therefore, the exact, approximate and simplified solutions have been proposed. The simplified solution, when the overlap rate in the slip mode is neglected, allows for prediction of contact separation angle analytically and affects only the coupled slip and sliding displacement evolution. It greatly simplifies the coupling of slip and sliding modes providing the required analytical predictions. The coupled effect of slip mode induces growth of the overlap and of total sphere center displacement until the contact separation. The verification of this effect is further required using FEM 3D solution.
The problem of a multiple contact interaction, particularly the simultaneous contact separation and activation of the sphere with the regularly packed granular bed has been considered analytically for the case of plane interaction. This class of problems has so far not been considered in the literature. Therefore, the incremental expressions specifying the sphere center trajectory angle and displacements have been proposed. It was demonstrated that the contact separation process occurs within the combined slip-sliding mode, essentially affecting the contact traction distribution, while the contact activation proceeds in the slip mode.
The analysis of mechanical response of a granular aggregate in its equilibrium configuration under the static or dynamic perturbations-induced small slip or sliding amplitudes should be carried out with the account for the multiple contact slip-sliding interaction effects discussed in the paper. On the other hand, a finite sliding mode is only slightly affected by the slip displacement, so some simplified models could be applied, or the slip effect even neglected.
The analytical treatment and derived formulae generate a class of benchmark problems, which can be used in assessing accuracy of the discrete element incremental schemes, especially, for the contacting sphere translatory motions. Meanwhile, the experimental applicability of the proposed approach was also discussed.

Compliance with ethical standards
Conflict of interest The material has been developed without any actual or potential conflict of interest including any financial, personal or other relationships with other people or organizations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.