A finite sliding model of two identical spheres under displacement and force control. Part II: dynamic analysis

In Part I of the present study, the static analysis for the sliding of two identical spheres under displacement and force control was carried out. For linear and circular sliding trajectories, the contact traction evolution was analytically specified for both monotonic and reciprocal sliding regimes. Similarly, for the specified gravity loading, the driving force evolution and the sliding path were also determined. In the present Part II of the analysis, the dynamic response for the same sliding modes is presented. The contact traction and velocity evolutions are considered in detail. The analytical formulae are proposed for prediction of the tangential restitution coefficient, critical velocity, and time of contact for the displacement and load-controlled motions. The effects of loading and reloading in reciprocal sliding are also considered with account for the slip and sliding regimes. The generated results have practical aspects and can be implemented in modeling of the asperity and rough surface interaction, wear analysis and also in the development of the numerical discrete element method.


Introduction
Dynamic interaction of contacting structural elements is of fundamental importance in numerous engineering applications, including friction and flash temperature of rough surfaces, wear, erosion, or particle deposition. Any modeling of these interactions requires a detailed understanding of the mechanics between two individual bodies under various external actions. A general analytical Mindlin's solution [1] for the slip regimes of two contacting spheres under normal and tangential forces is still incomplete as the sliding regimes occurring for larger displacement values require a separate treatment.
Alternatively, to be numerically applicable in the dynamical simulation of multiparticle systems, the simplified tangential and normal force-displacement models describing the impact of two spheres have been adopted. In most cases, the simplified models enable us to do the analytical solution for the range of experimentally accessible parameters, such as the coefficients of restitution and friction or the period of contact, and are of sufficient accuracy relative to the measurements [2].
The pioneering analysis within the simplified model assumption was developed by Hertz. Following the theory of elastic contact, an approximate solution was derived for the normal impact of elastic bodies in terms of the nonlinear relationship between the normal force and overlap between the contacting bodies. Following the Hertzian contact law, the duration of contact of two identical spheres, as presented in [3], is proportional to the radius of the sphere and inversely proportional to the initial velocity of order 1/5. The accuracy of the Hertz model was demonstrated by the FEM analysis, cf. [4].
During contact interaction, a loss of kinetic energy occurs. In the case of normal impact of elastic sphere with elastic substrate, the investigations were conducted by Hunter [5] and Wu et al. [6]. They found the amount of kinetic energy of the ball that is converted into elastic waves. In the case of normal impact for elastoplastic spheres, the dissipation of kinetic energy occurs mainly due to plastic deformation, while the energy loss due to stress wave motion is very small [7]. The coefficient of tangential restitution for the dissipative collision of frictional viscoelastic spheres was semi-analytically investigated by Schwager et al. [8]. They numerically simulated the tangential force evolution by using a linear spring (the model by Cundall and Strack [9]) in the direction of the contact plane, while the cut-off force was described by Coulomb's friction law. It was shown that the dynamics of the system depends on three independent parameters, including material constants and initial velocities of the collision. However, the detailed dynamics of the collision may be rather complex, revealing multiple transitions between (dissipative) Coulomb friction and (conservative) loading and unloading of the spring.
Thornton [10] proposed a theoretical model for the collinear impact of two elastic-perfectly plastic spheres. In view of the transition from elastic to fully plastic impacts, the analytical solution for the coefficient of restitution was also given in this work. Recently, the energy absorption in compression and impact of dry elastic-plastic spherical granules has been theoretically and experimentally analyzed by Antoniuk et al. [11]. From a single-particle compression and impact tests, they measured Young's modulus, contact stiffness, yield stress, restitution coefficient, and strength of spherical granules ( γ -Al 2 O 3 , zeolites 4A and 13X, sodium benzoate). In the elastic range of compressive loading, the contact Hertz model was tested and the analytical contact model was proposed to describe the elastic-plastic nonadhesive deformation beyond the yield point of a force-displacement curve. In the case of repeated loading at constant load amplitude, the analyzed granules exhibited cyclic hardening with an increase in the restitution coefficient up to a certain saturation level during the plastic deformation process. A model describing the increase in the contact stiffness with the number of cycles was then proposed.
During an oblique impact, when the tangential forces dominate, the contact response depends not only on the initial values of the normal and tangential forces, but also on the loading history applied. To simplify the matter, instead of a deformable, the rigid-body approach was introduced by Goldsmith [12] to predict the dynamical behavior of two objects during impact. A theoretical model was developed and based on the impulse-momentum law, not accounting for the influence of the contact deformation and developed tractions. In particular, Maw et al. [13] demonstrated that this technique does not accurately predict the oblique rebound response, particularly at small impact angles, where the tangential contact velocity reverses the direction during the impact. In Maw's model, an elastic oblique impact was considered describing the normal forces by Hertz contact theory without any energy loss in this direction, resulting in the coefficient of restitution equal to one. By examination of the FEM-based results from the rebound kinematics of elastic and elastoplastic spheres, a simple semi-analytical model accounting for impact velocity, impact angle and degree of plastic deformation, was proposed by Wu et al. [14].
In general, for the tangential direction, the variety of simplified models of different complexity can be found in the literature. They mainly incorporate three basic elements, i.e., spring (linear or nonlinear deformation), dashpot (viscous damping), and slider (frictional damping). Probably the earliest approach was proposed by Cundall and Strack [14] by implementing a linear-damped tangential spring in combination with a sliding friction. The tangential spring stiffness k s was selected proportionally to the stiffness k n of the spring in normal direction by varying the ratio within the range of 2/3 ≤ k s /k n ≤ 1, following from Mindlin's model bounds [1]. A damping term should also be included to dissipate the energy and eliminate the perpetual oscillations in tangential direction due to tangential spring usage alone. The introduction of a damping term can also assist in eliminating the spurious oscillations occurring in tangential direction when modeling the sliding regimes. Instead of a constant tangential stiffness k s , Kohring [15,16] proposed a most general formula for the tangential spring stiffness, assuming k s to depend on the overlap evolution between the particles during contact. In this case, k s was derived from Mindlin's theory by assuming that no slip on the contact surface occurs. In general, these two approaches match sufficiently accurately the experimental measurements by Foerster et al. [17], for the oblique rebound response of the sphere. However, the models are sensitive to the tangential spring stiffness, when the friction coefficient is large. In this case, there is a limit of the tangential spring model, when the discrete particle simulation cannot adequately reproduce the macroscopically representative wall pressure near the hopper bottom after the end of filling [18]. In addition, the simplest models accounting for the single sliding friction law or a single viscous damping term only, or a viscous damping combined with the sliding friction rule, have been implemented in various simulations predicting the granular material behavior [19][20][21][22][23].
Along with the linear spring models, variable-stiffness spring models with a sliding friction rule were proposed. In particular, Langston et al. [24] again adopted the theory of Mindlin deriving the simplified equation for evaluation of the friction force relying on the normal force described by the continuous potential interaction model with inverse power potential choosing to be n = 36 to obtain a compromise between excessive overlap without the sphere being too stiff. In this model, the friction force was described by the same equation, which produced the friction force resulted from the sliding or slip mode in terms of the maximal tangential displacement dependent on the overlap. As a result, the macroscopically representative wall pressure was reproduced by discrete element simulation of the filling process in hoppers. The similar nonlinear model, additionally accounting for the hysteretic behavior, was proposed by Walton and Braun [25]. Following Mindlin's approach for a constant normal force, Tsuji et al. [26] proposed a model based on the tangential spring of nonlinear stiffness related to the normal displacement. In the model by Maw et al. [13], a functional expression was adapted to describe the tangential displacements and allowed a stepwise solution without memorizing the whole loading history in Mindlin and Deresiewicz [27] theory. The accuracy and differences of the above tangential force schemes were precisely demonstrated by Schafer et al. [22] in comparison with the experimental measurements for the oblique rebound response of the sphere [17].
In Part 1 of the analysis [28], a finite sliding model of two identical spheres under the displacement and force control has been proposed. The model was based on the approximate prediction of a real deformation of the spheres by the overlap approach along with the Hertz formula for normal contact force and Coulomb's law for sliding friction. The analytical relationships for the contact force and sliding path were derived, and the hysteretic behavior was also taken into account assuming an instantaneous switch of the contact force onto a reversal sliding path or a gradual unloading/reloading process relying on the slip displacement with transition to the sliding mode.
In the present part of the analysis, the proposed model is implemented to predict the dynamical behavior of two identical spheres under displacement and force control. The simulation is mainly devoted to an explicit analysis of the main dynamical parameters as well as providing the relations for the coefficient of restitution, critical velocities, time period of contact interaction, etc.
The paper is organized as follows. Section 2 is devoted to the displacement control approach, with the numerical and explicit analysis performed to study the dynamics of two spheres in contact for the linear and circular motion trajectories. In Sect. 3, the dynamical actions during the sphere sliding under the mixed control are presented, while the concluding remarks are stated in Sect. 4.

Governing relations
Consider the dynamics of two identical spheres in contact in a displacement controlled motion. Let the upper sphere center be attached to a moving tool imposing the translation along the linear or circular path with a constraint set on rotational motion. Let the sphere, moving with the initial velocity v s (t 0 ), come into an elastic, isolated contact with a bottom stationary sphere. Due to contact, a force N normal to the contact plane and the friction contact force T = ±μ |N | arise. These forces, acting at the contact interfaces, produce the resultant force T * , acting on the moving sphere center, which reduces the imposed velocity v s (t 0 ) due to dissipation of the kinetic energy generated by this velocity. The resultant normal force component N * is not analyzed as it does not affect the sphere motion and becomes the tool reaction.
If the kinetic energy is not fully dissipated during the contact, the contacting spheres can separate at a certain time instant t c . The sliding paths during the sphere motion were schematically illustrated in Figs. 4 and 9 of Part 1 for the linear and eccentric circular trajectories.
Thus, the equation of motion of the upper sphere center path specified by the coordinate s t ≡ s (t) at any time t can be written as follows: where T * (s t ) ≡ T * (s (t)) is the resultant contact force component induced from the imposed trajectory of the moving sphere, m denotes its mass (the mass of moving tool elements can also be included).

Case 1: a linear path of the sliding sphere
For the imposed linear trajectory, the upper sphere center moves with the specified initial velocity v s (t 0 ) relatively to the fixed sphere center (Fig. 4, in Part 1). The imposed trajectory induces the sliding contact path, which was defined in the Cartesian coordinate system (see, Fig. 4, in Part 1). For the sphere center, this path is described by the evolution of the horizontal coordinate bounded over the range of s t ∈ [±s u , ∓s u ] (where s u = ± 4Rh 0 − h 2 0 are the ultimate horizontal coordinates of the upper sphere center path). The sign in s t bounds depends on the direction of motion. In particular, for the leftward motion of the upper sphere center, the upper signs refer to the contact engagement at the point x = +s u and separation at the point x = −s u . For the rightward motion, the sliding path starts at the position x = −s u and ends at x = +s u . In this case, the lower signs in the bounds description are accounted for. A sliding time period (or time of contact) will be denoted by t c .
The resultant horizontal force (see, Fig. 8, in Part 1), induced by the linear path of the upper sphere at the time instant t, was provided by the relation where s g = sign (v s (t)) is the function evaluating the friction force direction according to the moving sphere is the normal force to the contact plane, determined in Part 1 (see, Fig. 8 there), y c = 2y 0 = 2R − h 0 = const. is the fixed vertical coordinate of the upper sphere center for the linear trajectory (see, Fig. 4, in Part 1), h 0 represents the maximal overlap between two spheres of radii R resulted from the normal force N (0), while k n = 1 is the stiffness of the nonlinear Hertzian spring depending on the Poisson's ratio ν, the coefficient of friction between the spheres is denoted by μ and elastic modulus of sphere by E (see, Part 1). At the contact engagement or separation points there is T * (±s u ) = 0, and T * (0) = −s g μ |N (0)| for the nominal contact plane.
The resultant of vertical force (i.e., the moving tool reaction force) versus the sphere center path, following Part 1 analysis, can be defined by formula First, we solve Eq. (1) numerically. A 5-th order Gear predictor-corrector scheme [29] was implemented for the time integration of this equation. Following this scheme, all dynamical parameters are predicted in a predictor procedure and subsequently corrected by the acceleration increment in the corrector procedure, cf. [30]. In Fig. 1, the evolution of the main dynamical parameters for different values of the initial velocity v s (t 0 ) is illustrated for the rightward motion. The initial velocity was related to its critical value v s,crit for which the sphere stops exactly at the contact separation point. The critical velocity determination will be given in the next section. Here, the unit values k n = 1 N/m 3/2 , R = 1 m and μ = 0.6 were simply selected for the current illustration. This parameter selection corresponds to the nondimensional force and length characterization by definingN = N /k n R 3/2 ,ĥ 0 = h 0 /R,ŝ t = s t /R, etc. The hysteretic response, presented in the next section, was also included into the analysis.
As can be seen in Fig. 1a, the path coordinate s(t) grows in time depending on the value of the initial velocity. For increasing initial velocity, the character of function s(t) becomes linear versus time. In this case, along the path, the contacting sphere looses its kinetic energy due to the frictional dissipation and sphere center velocity v s (t) decreases in time (Fig. 1b). For v s (t 0 ) < v s,crit , the moving sphere stops and the contact separation does not occur. When the sphere completely stops, the resultant of horizontal contact tractions acting at the sphere center vanishes, T * = 0, but the vertical resultant force does not vanish, N * = 0 (cf. the case v s (t 0 ) = 0.8v s,crit , in Fig. 1c) because the sphere cannot separate from contact. During the sphere slowdown, the hysteretic response (cf. the case of v s (t 0 ) = 0.8v s,crit in Fig. 1d) takes place in numerical computation. The details of the numerical implementation of hysteretic behavior will be considered later in the paper. The  Fig. 1c). As can be seen, the contact time decreases with the increase in the initial sphere velocity.

Case 2: a circular eccentric path of the contacting sphere
Consider the upper sphere sliding along a circular eccentric trajectory relatively to the fixed bottom sphere. Let the upper sphere center be attached to a moving tool imposing the circular path of radiusR with the tangential velocity v s (t 0 ) applied (Fig. 2). Due to contact, the force N normal to the contact plane (denoted by AB in Fig. 2) and the friction contact force T = ±μ |N | arise. By acting at the contact center point C, these forces produce the resultant tangential force T * , acting at the moving sphere center along its trajectory and reducing the imposed velocity. The resultant force N * acting along the vectorñ, normal to the trajectory, represents the reaction of the moving tool and does not affect the motion. The upper sphere center motion is described with respect to the axisτ in terms of the curved path coordinates t , and the contact forces N , T = ±μ |N | are projected onto the trajectory axesñ andτ (Fig. 2). Let us determine the angleβ t =α t −γ t generating the inclination of contact forces relative to the eccentric circular path (Fig. 2). Thus, due to the eccentricityẽ, the curved sliding path coordinates 1t , generated at the contact center C during the upper sphere sliding, is governed by the variable radiusr t = const. (Fig. 2). On the other hand, the curved sliding path coordinates t at the upper sphere center is controlled by the radius R 1 = 2r t and can be defined ass t = α t 0 R 2 1 + dR 1 dα 2 dα = α t 0 4r 2 t + 4 dr t dα 2 dα = 2s 1t . Hence, the curved sliding path coordinate developed at the moving sphere center is a twice of that at the contact interface. Geometrically, it can also be proven that rectangular coordinates describing these paths are interrelated

Circular path with eccentricity
Curved sliding path at contact interface

Fig. 2
Motion the contacting sphere along the imposed circular trajectory with eccentricity asx t = 2x 1t andỹ t = 2ỹ 1t (Fig. 2). The angleγ t between the tangential axisτ for the circular path of the sphere center and the horizontal axis can be expressed asγ t = ± arctan x t y t +ẽ = ± arctan Similarly, the angleγ 1t between the tangential axisτ 1 for the curved sliding path at the contact interface and the horizontal axis was given in Part 1 of the analysis and was expressed asγ 1t = ± arctan t =γ t and finally,β t =α t −γ t =β 1t . Hence, the contact forces, acting at the contact interface at point C, are inclined relative to the sliding path at the same angle as those acting at the upper sphere center path (Fig. 2). Hence, for definition of the resultant tangential force acting on the upper sphere center, we can apply Eqs. (41-42) and (48-49), given in Part 1, as follows: is the normal force at the contact plane referred to the sphere center path coordinates t , given in Part 1. In this case, the resultant normal force acting on the moving tool was defined by the following expression Here,s t ∈ [−s u ,s u ] describes the generalized coordinate of the curved sliding path as the arc length evolution at the instant t, for the path radius 2r t and the angleα t ∈ [−α u ,α u ] (see, Fig. 10 and Sect. 3.3, in Part 1). Following relations in Part 1, the ultimate values of this coordinate at the contact engagement/separation points were defined ass u =R(α u − arcsin(k sinα u )) (whereα u = ± arccos(( 1 is the angle at the sphere engagement/separation points). The eccentricity ratiok =ẽ/R describes the shape of the circular path. In particular, there are two limit cases: a perfectly circular trajectory, whenk = 0, and the linear path, fork = 1. In the latter 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 formulaẽ Also, fork → 1, we getα t ≈ α t and the eccentric circular path transforms into the linear sliding trajectory.
For brevity, in Fig. 3, the resultant tangential force evolution in time for the linear and circular paths specified by different eccentricities was calculated numerically and plotted. In these plots, the initial velocity was selected within the range of v s (t 0 ) /v s,crit = −[1, 1.5, 3], for the leftward motion (where v s,crit is used from the linear path of the contacting sphere).
As can be seen in Fig. 3, with increasing the eccentricity ratiok, the circular path tends to the linear path (cf. curve fork = 0.8 and curve for the linear path). In the case ofk = 0.3 and v s (t 0 ) = v s,crit , the sphere moving along the circular path cannot reach the contact separation point because of too small value of the initial velocity. As the kinetic energy tends to zero, the unloading/reloading force-slip displacement curves were taken into account resulting in the slowdown of the moving sphere, and the tangential force drops to zero.

Hysteretic response for reciprocal sliding
The static force-path diagrams accounting for the hysteretic response were presented in Part 1 [28] of the analysis. These diagrams allow us to evaluate descending/ascending branches of the unloading/reloading forces-path functions in terms of the slip displacement. The hysteretic response curves mainly generate the kinetic energy dissipation during the unloading/reloading loops and can be used instead of velocity-dependent damping, when the sphere-sphere contact is treated as viscoelastic.
In numerical prediction of the hysteretic response during the contact dynamics, it is important to recognize the loading, unloading, or reloading regimes. When all dynamical parameters are tracked over the time integration of Eq. (1), any change in direction of the tangential (or horizontal, in the case of linear path) velocity can be treated as a trigger for the unloading or reloading regimes. Thereby, these regimes can be recognized by the condition of change in the velocity orientation expressed by the sign conversion, i.e., sign (v s (t)) = sign (v s (t + t)) (where t is the time step). If this condition is first satisfied, the unloading regime starts, and if the sign conversion is detected again, the reloading process takes place.
For the dynamic analysis, the unloading/reloading force-path curves may be implemented in terms of the analytical formulae (58) and (59) presented in Part 1. The unloading or reloading force-path curves at any time t + t can be generally expressed as follows: where are the values of friction force and path coordinate, at the onset of the unloading or reloading regime, )| is the final value of friction force on the unloading or reloading curve at point ) is the ultimate slip displacement at the unloading initiation point s A determined from M-D theory in Part 1.
An incremental form of the unloading/reloading force-path curve applied in Eq. (2) yields the following: This equation can be also reexpressed in a velocity-dependent fashion by substituting ( In view of unloading or reloading regimes, the resultant contact force tractions are predicted by Eqs. (1a-1b) by using the substitution −s In Fig. 4, the hysteretic response of two spheres in contact is demonstrated for the linear path. In the calculations, the unloading/reloading regimes are modeled by the nonincremental form (Eq. 2) and those defined by the incremental form (Eq. 3). To get the unloading/reloading cycle, an initial velocity v s (t 0 ) was selected to be 0.7v s,crit allowing for the slowdown of the moving sphere. As can be seen in Fig. 4a, the prediction of the unloading/reloading regime according to the nonincremental equation results in spurious oscillations during the sphere slowdown. In particular, at time instant t = 30 s, instead of the tangential force vanishing to zero, the spurious oscillations suddenly appear. They originate from the instantaneous jump from the unloading to the reloading curve. To attain the gradual transition of these curves, a very small time step is required. This step should be even less than the error of double-precision variables in most cases, when the sphere motion tends to slowdown. In contrast (see, Fig. 4b), application of the incremental form for the unloading/reloading curves effectively damps the oscillations resulting in vanishing tangential force and complete stop of the sphere.
The hysteretic behavior with neglect of the slip displacement was also considered (see Fig. 4c). In this case, the hysteretic behavior occurs due to a velocity sign change in Coulomb's friction law. As can be seen, the application of the pure Coulomb's friction law leads to stepwise oscillations in T * (t) and N * (t) functions during the unloading/reloading cycle (Fig. 4c).
The effect of the power factor p is plotted in Fig. 4d. Here, a decrease in p values from 1/2 up to 1/4 produces a more rapid damping on the unloading/reloding cycle than at p = 1/2. For example, in M-D theory, the function power of p = 1/3 was derived analytically.
Finally, the damping by the hysteresis loop may be demonstrated by assuming that the upper sphere has not a single-contact sliding direction, i.e., at a certain time instant during contact, the upper sphere is instantly forced to move in the opposite direction (Fig. 5a). In this case, the reversal of relative contact velocity was considered in simulation. In fact, a real sphere sliding process requires application of an impulsive loading reducing the sphere velocity and gradually accelerating to the reverse motion in a very short time period. In the present analysis, this period is neglected and a numerical procedure is applied to the case of instantaneous In Fig 5a, it was supposed that the sphere moving with the initial velocity of v s (t 0 ) = 1.05v s,crit comes into contact with a stationary sphere at time t 0 ; then, at a certain time instant, the sphere motion is changed instantly due to the applied impact forcing the sphere to move in the opposite direction. Finally, after a small time period, the sphere is subjected to another impact resulting in the instantaneous change of velocity direction (Fig. 5a). Due to these impacts, the unloading/reloading cycles can be effectively demonstrated (Fig. 5c, d). As can be seen, the changes in velocity direction result in the hysteretic loops on T * and N * evolution functions (Fig. 5c). These loops are determined by the incremental form of Eq. (3). During the unloading-reloading cycles, the sphere looses its kinetic energy and finally stops, remaining in contact for vanishing tangential force and fixed normal force (Fig. 5d).
Integrating Eq. (4) along the sliding path, we obtain where W (s t ) = s(t) s(t 0 ) T * (s t ) ds t is the work done by the sliding sphere from point s (t 0 ) to point s (t) under the action of the resultant tangential contact force T * (s t ).
Using the T * (s t ) relation for the linear path imposed, we can combine the work W n (s t ) done by the normal contact force (not perpendicular to the sliding path, see, Fig. 8, in Part 1) with the work W μ (s t ) lost due to the friction force. Thus, Inserting N (s t ) into Eq. (6), the following expressions are obtained: Due to the presence of the integrand of the form 2R − √ • 3/2 , Eqs. (7-8) cannot be expressed in terms of elementary functions. In order to eliminate the square root within the power function, we rearrange the polynomials (7) and (8) into the trigonometric formulae. From Fig. 4 in Part 1, the function of overlap can be expressed by the formula where Now, the sphere center path coordinate and its differential are related to the angle α t by the following expressions: Combining Eq. (7) with Eqs. (9) and (10), the work of the normal force can be expressed in the following form: Integration of expression (12) is simple, providing the following formula for the work done by the normal force at any s t ≡ s (t):  9) and (10), the integrand can be rearranged into the following trigonometric expression: Expression (14) is an elliptic integral of the form R x, √ ax 3 + bx 2 + cx + e dx (where R denotes a rational function of its arguments). Integrals of this type are not integrable in elementary terms, but can be reduced into relation involving the integrals of rational functions and Legendre's canonical form integrals of the first, second, or third kind [31]. The forms of these integrals and results are well studied (e.g., [31,32]).
Thus, after the sequence of transformations, the work of friction force can be expressed by the following relation:  Finally, the horizontal velocity for the linear path is expressed from Eq. (5) as follows: where The sign "±" specifies the horizontal velocity direction according to the coordinate system used. In Eq. (16), it is also convenient to introduce the function s g = sign (v s (t 0 )) instead of "±" to analytically fulfill the sign conversion due to velocity direction change.
The graphs of the analytical tangential velocity function and those calculated numerically are plotted in Fig. 6.
As can be seen in Fig. 6, the analytically predicted velocities according Eq. (16) along the sphere center path are consistent with those computed numerically.

Specification of the coefficient of restitution and critical velocity
For two spheres in contact, the imposed initial dynamic parameters, such as position, velocity, can vary during the contact interaction of two spheres. Measurements of these parameters are conducted mainly on the macroscopic level and are very limited for the grain-grain interaction.
The coefficient of restitution is a very useful parameter characterizing the quantity of the kinetic energy dissipation in the contact sliding motion. The coefficient of restitution is generally given by the formula where are the initial and final values of the velocity (the upper signs refer to the contact engagement point at the sphere center x = +s u and its separation at the point x = −s u , for the leftward motion; for the rightward motion, the path starts at x = −s u ending at the point x = +s u , and the lower signs account for this case). Using expression (16), we can get Generally, e s > 0 corresponds to a monotonic motion along the sliding path, and e s < 0 corresponds to a case of velocity reversal. The value e s = 0 specifies the complete loss of kinetic energy before the contact separation.
When the contacting sphere moves from the contact engagement to separation points, i.e., s (t) : ±s u → ∓s u , the normal force work expressed by Eq. (13) turns to W n (±s u ) = 0 (i.e., no potential energy accumulation) and only the work dissipated by the frictional force remains. This allows for simplification of the analytical expression for the coefficient of restitution. In particular, for the frictional force work specified by Eq. (14), we can impose the symmetry conditions, because the frictional force evolves symmetrically along the path. Hence, in Eq. (15), we can adopt a range of integration over α t : 0 → α u and use the multiplication factor 2, for W μ (s t ). In this case, the incomplete elliptic integrals of the first and second kinds simply turn to the complete elliptic integrals, such that, K (k) = F (1, k) and E (k) = E (1, k), for s u > 0.
Applying these conditions in Eq. (15), we finally obtain Here, the term −s g was simply replaced by "−", since the work of the frictional force is negative, and the complete elliptic integrals are used only. Now, inserting Eq. (19) into Eq. (18), we obtain formula for the restitution coefficient for the linear path of the contacting sphere center, and thus, Thus, the restitution coefficient due to the frictional damping during sphere sliding mainly depends on the initial velocity, mass of the moving sphere, its radius, friction coefficient, contact stiffness, and the sphere path factor k r = y 0 /R = 1 − h 0 2R specified by the maximal overlap h 0 . Finally, we can analytically define the critical horizontal velocity v s,crit , when the moving sphere dissipates all kinetic energy exactly at the contact separation point. Thus, when v s (t 0 ) < v s,crit , the sphere stops without the contact separation. Thereby, holding for v s (t c ) = v s (±s u ) = 0 in Eq. (16) and using Eq. (19) we get the following: In Fig. 7, the function for the restitution coefficient predicted by Eq. (20) is plotted versus the initial velocity, for μ = 0.6 and different values of the initial velocity of the contacting sphere. As can be seen in Fig. 7, the restitution coefficient is strongly dependent on the initial velocity parameter. With increasing initial velocity, the restitution coefficient increases as well. For the imposed high value of initial velocity, we have e s → 1. Also, the values of restitution coefficient calculated analytically and those are numerically coincident, thus proving the correctness of the formulae derived.

Specification of the time of contact
The duration of the contact interaction t c (time of contact) during the upper sphere sliding with respect to the bottom stationary sphere is a very important parameter in modeling of the particulate systems. For example, in a correct DEM simulation, the assumed time step for integration of the motion equations should be much lower than the duration of contact for the normal and tangential directions.
To find the time of contact interaction of two spheres, double integration of Eq. (4) is required; thus, It is obvious that Eq. (22) is analytically nonintegrable because the function of W (s t ) contains the elliptical integrals occurring in Eq. (15). In this case, the numerical solution is feasible by integration over the small intervals (e.g., via the trapezoidal sum). A small increment for the argument is required, and the numerical integration is quite inconvenient. Thereby, the quadrature formulas of Gauss can be adopted resulting in very accurate approximations, provided the interpolation nodes are chosen in an appropriate way. In the current analysis, we use a quadrature operating with five interpolation nodes with known approximation error [33]: = 0.703 h 10 9 1 f (ξ ) (8) , a < ξ < b, (24) where h = b − a, c = (b + a)/2. In Eq. (24), the 8-th derivative is denoted as (8). Thus, the approximate integration of Eq. (22) yields the general formula for the contact time: where v s (•) is the velocity function specified in Eq.
are the initial and final values of this function. Results for the contact time determined by the numerical analysis and those predicted by the approximate formula (25) are plotted in Fig. 8.
As can be seen in Fig. 8, the quadrature-based formula (25) yields a very accurate prediction. A mild deviation from the numerical prediction is observed for the initial velocity near the critical value.
2.5 Analytical solution for an eccentric circular path of the contacting sphere

Determination for the tangential velocity
As in the preceding section, the tangential velocity of the sphere, moving along the imposed circular path with eccentricity, may be found by integrating the total work expended by force T * (s t ) over the curved paths t at the sliding sphere center. Thus, the work of force N (s t ) (normal to the contact plane, but not perpendicular to the curved sliding path) can be expressed as follows: A substitution ofk (27) allows for a simple integration giving the following formula: The obtained formula (28) providesW n (s t ) = 0, but when the trajectory proceeds along the finite path from the contact engagement to separation points, i.e.,s t : ±s u → ∓s u , thenW n (∓s t ) = 0 yielding the complete loading/unloading cycle durings t : ±s u → ∓s u . The work dissipated due to the frictional force is Again, Eq. (29) represents a nonintegrable integral of the form R x, √ ax 3 + bx 2 + cx + e dx. It can be expressed in terms of the elliptic integrals by the sequence of suitable technique of transformations. However, in this case, the performed transformations into elementary functions were unsuccessful for the integrand. Hence, we adopted the approximate integration of Eq. (29) via quadrature (23) and get the following relationship: where h =s t − (±s u ) , c = (s t + (±s u ))/2, andT μ (s t ) = −s g k n μ h t (s t ) is the frictional force developed at anys t , such thats t ∈ ±s t , ∓s t . The signs are accounted for in the same manner as in the case of linear path with respect to the leftward or rightward motions. Finally, the function of tangential velocity againsts t can be determined by applying Eqs. (16), (28), and where v s (t 0 ) ≡ v s (±s u ) is the initial tangential velocity applied at the center of sliding sphere. In Fig. 9, the graphs plotted for the tangential velocity function againsts t compare the results obtained by Eq. (31) with those determined numerically. The analysis was conducted for different eccentricity ratios, specified by valuesk = 0.3,k = 0.8 and different initial velocities. The initial data used here were the same as for the above numerical analysis.
As can be seen in Fig. 9, the approximate prediction by formulae (30) and (31) for the tangential velocity against the curved sliding path of the moving sphere center along a circular trajectory with eccentricity is perfectly consistent with that computed numerically.

Specification of the critical velocity, coefficient of restitution, and contact time
The critical velocity and coefficient of restitution may be predicted from the total work expended by the moving sphere from the contact engagement to separation points, i.e.,s t : ±s u → ∓s u . Thus, for the finite path s t : ±s u → ∓s u ,W n (±s u ) = 0, and the work dissipated by the frictional force contributes to the critical velocity only. It allows us for the approximate prediction of the total work in the following form: Explicitly, Eq. (32) can also be expressed by rearranging Eq. (29) into a more simple form for the numerical integration:W where cot β t = arccos r 2 +R 2 1−k 2 2r tRt ,r t (0) ≡r t (α t = 0), andr u ≡r t (α u ) are the radii specified by equationr t (α t ) =R 1 2 =R 2 −k cosα t ± 1 −k 2 sin 2 (α t ) , derived in Part 1. Note that for the curved sliding path at point ofα t = 0, the angleβ t = 0, and the function of cot β t is singular. In this case, when the numerical integration is applied, this point should be slightly shifted to omit the singularity. Now, using Eqs. (21) and (32), we can obtain an approximate formula for the critical tangential velocity, when the moving sphere along the eccentric circular trajectory loses the kinetic energy exactly at the contact separation point; thus, In Fig. 10, the graphs, plotted for the tangential restitution coefficient versus initial tangential velocity, provide the comparison of the results obtained by formula (35) with those determined numerically. The analysis was conducted for different ratios of the eccentricity (k = 0.3,k = 0.8) and different values of the initial velocities. The function for the restitution coefficient determined for the linear path of the upper sphere center was also included.
As can be seen in Fig. 10, the approximate prediction by formula (35) perfectly coincides with that computed numerically. For the eccentricity ratiok = 0.8, the coefficient of restitution tends to that predicted for the linear path, because fork → 1, the circular path turns into the linear path. Also, as can be seen, the increase in the initial velocity reduces the energy loss due to the frictional forceT μ (s t ), and in turn, the restitution coefficient increases.
The time of contact, when the trajectory of the moving sphere is defined by a circular path with eccentricity, can also easily be predicted by expression (25). In this case, five values of the tangential velocity function predicted by Eq. (31) should be substituted into Eq. (25). These velocity function values are as follows: , v s (s t = 0.6546 |s u |), and v s (s t = −0.6546 |s u |).
The results of the contact time predicted by the numerical analysis and those determined by the approximate formula (25), when the trajectory of the moving sphere is defined by an eccentric circular path, are now plotted in Fig. 11. The effect of eccentricity is demonstrated accounting fork = 0.8 andk = 0.3 (Note that the tilde above a symbol was omitted in Fig. 11 as well in all above plots). The contact time predicted for the linear paths of the contacting sphere was also included.
As can be seen in Fig. 11, the approximate formula (32) results in very accurate prediction, when the trajectory of the moving sphere is defined by a circular trajectory with eccentricity. For eccentricity ratiok → 1, the contact time tends to the value corresponding to the linear path of the moving sphere (cf.k = 0.8, in Fig. 11).

Finite sliding of rigid spheres for the force controlled motion
In this section, the vertical load (e.g., gravity) is assumed as fixed and next the horizontal driving force control is applied in the dynamic analysis of finite sliding of two identical spheres. Due to the presence of nonintegrable functions describing the deformation of spheres in terms of the Hertz overlap formula, the main contact dynamics parameters for a finite sliding of the deformable spheres cannot be derived explicitly (see, Part 1). In this case, the prediction of the main contact dynamics parameters is limited here to the case of finite sliding of the rigid spheres. Let the vertical load Q be applied first and kept constant, next the horizontal displacement u 0 is induced by the pulling device with the conjugate horizontal force T 0 at time t 0 (Fig. 12). Suppose that the application of vertical force Q and horizontal displacement is accompanied by the constraint preventing the rotational motion of the sphere. When driving force T 0 overcomes the static friction force, i.e., |T 0 | > T μ , the upper sphere starts to move relative to the lower fixed sphere. From the static analysis (see, Sect. 4, in Part 1), the static friction force had its maximal value at the contact engagement, i.e., T μ = T μ (±s u ), next decreasing in the course of progressive sliding along the contact surface. Thus, the upper sphere can accelerate, and its tangential velocity v s (t) varies with time.
Consider the contact dynamics in the n−τ coordinate system. The direction of axes n and τ is specified to be coincident with the Cartesian x, y axes when the angle equals to α t = 0 (Fig. 12).
The unbalanced tangential and normal forces are defined as a sum of their projections onto τ and n axes. As can be seen (Fig. 12), the projection of vertical force Q onto τ -axis and the projection of driving force T 0 onto n-axis change their orientation along the sliding path. Accounting for the sign conversion, the sum of force projections can be expressed by the following formulae: where s g = sign (T 0 )and α t = s t /(2R) is the angle at time t (Fig. 12).
During the upper sphere center motion relative to the fixed sphere along a perfect circular path of radius R, the orientation of tangential velocity vector v s varies in time yielding the change in velocity v s (Fig. 12), which induces the centrifugal acceleration a n pointed toward the center of the circular sliding path. In view of Eqs. (36) and (37), the application of D'Alambert principle yields the following system of equations: The second equation of the system (38) can be adopted to define the contact reaction force in terms of the static and inertial contact forces as follows: where N stat (α t ) = −T 0 sin α t + Q cos α t is the static reaction force and Φ n = −a n m = −m v 2 s (t) 2R is the inertial reaction force of contact.
The inertial reaction Φ n is pointed oppositely to the centrifugal acceleration a n and acts as a repulsive contact force growing in time due to increasing in the tangential velocity. When |Φ n | = |N stat (α t )|, the contact separation occurs and the tangential contact velocity at the separation point can be defined as follows: where the sign ± can be replaced by sign (T 0 ). The contact separation coordinate, defined at the sphere center, s c and contact time t c can be obtained by solving the system of Eq. (38).

Numerical analysis
As previously, the solution of the system (38) will be primarily performed numerically via the 5th order Gear predictor-corrector scheme [29]. Plots of the main dynamical parameters are depicted in Fig. 13. The analysis was conducted for different ratios of the driving force T 0 and the static friction force T μ . Additionally, the effect of centrifugal acceleration a n on the distribution of dynamical parameters is also shown. The contact engagement was supposed to be at the point defined by the vertical coordinate y 0 /R = 0.86. The values for Q = 1 N, R = 1 m, μ = 0.6 were also simply selected.
As can be seen in Fig. 13a, the sliding sphere center coordinate s (t) progresses with time until the separation occurs. Since the driving force is applied, the sphere accelerates with increasing velocity v s (t) in time  Fig. 13b), because the static friction force T μ diminishes with progressing along the sliding path (cf. Fig. 19, in Part 1) allowing the unbalanced force T * (s t ) to grow (Fig. 13d).
The contact reaction static force N stat (t) diminishes (cf. Fig. 13c) after the contact engagement due to decrease in the contact plane orientation, when the sliding path runs around the bottom sphere center. However, the inertial contact reaction Φ n , caused by centrifugal acceleration, after contact engagement starts to grow in time. At the time instant t c , when |Φ n | = |N stat (α t )|, the resultant contact reaction becomes N (t c ) = 0 and the contact separation of the interacting spheres occurs (Fig. 13c). In fact, the inertial reaction force is repulsive and yields a shorter contact period in comparison with the case of a n = 0 (cf. Fig. 13a-c).

Analytical solution for the tangential velocity
Analytical solution of the system of Eq. (38) for a sliding path s(t) can be performed by specifying the function of tangential velocity versus path length and its integrals with respect to t. Substitution of formula (39) into the first relation of (38) yields the nonlinear differential equation of the second order: Introducing the velocity functions v s = ds t dt and d 2 s t dt 2 = v s dv s ds t , Eq. (41) reduces to Bernoulli first-order nonlinear differential equation of the following form: where F τ (s t ) = 3s g μ (Q sin (s t /(2R)) + T 0 cos (s t /(2R))) + 1 − 2μ 2 N stat (s t /(2R)) , s 0 ≡ s (t 0 ) = ±s u is the initial coordinate of the curved sliding path at the center of upper sphere. To match the direction of sphere motion, a velocity function ± √ • can simply be replaced by sign (T 0 ) √ •.
This coordinate at the separation point is specified by finding a root of Eq. (43) in combination with Eq. (40); thus, However, the root (44) cannot be expressed by radicals, and in this case, the numerical solution is only feasible. The graphs of the tangential velocity determined by formula (43) and those predicted numerically are plotted in Fig. 14, for the different ratios of the driving force T 0 to the static friction force T μ .
As can be seen in Fig. 14, at the contact engagement (s (t 0 ) = −s u ), for the rightward motion, the tangential velocity is zero, and the upper sphere starts to move under the applied driving force T 0 , while the tangential velocity increases with progressing the sliding path up to the separation point. As mentioned above, this increase results from the decrease in the static friction force T μ along the upper sphere center path. The calculated velocities according to Eq. (43) are consistent with those computed numerically showing the correctness of the above-obtained relationships.

Analytical prediction of contact time
Another integration of the tangential velocity function, , is required to find the contact time t c . However, function (43) is analytically nonintegrable. An approximate integration of the function v (s t ) −1 using the previously applied Gauss quadrature (23) is also not suitable in view of an improper integral with an infinite discontinuity at the point approaching the initial coordinate of the path, where v s (±s u ) = 0.
The tangential velocity graphs plotted in Fig. 13b show that the neglect of the centrifugal force effect, that is, the assumption of a n = 0 induces the higher velocity value at the contact separation, the increased separation distance and the contact time, in comparison with the case of a n = 0. Approximately calculating the contact time, let us supposes that the higher area of function v (s t ) −1 , in the case of a n = 0, could be adequately reduced to the area of v (s t ) −1 , in the case of a n = 0, if we use the separation coordinate obtained from (44). This allows us for avoiding the solution of the nonlinear differential equation produced by the inertial contact reaction Φ n and calculating the tangential velocity function by splitting variables in the total work integral.
Graphs of the contact time versus the applied driving force predicted numerically and those determined by formula (47) are plotted in Fig. 15. The initial data used in the calculation were described in Sect. 3.2.
In Fig. 15, the values of contact time elapsed from the contact engagement to its separation increase, when the applied driving force decreases. From these graphs, it can be seen that the relative error in approximate prediction of t c by Eq. (47) does not exceed (+1 − 2) % for the driving force close to T μ . The performed analysis has also shown that this error is dependent only on the coefficient of friction. In particular, with a decrease in μ, the error of t c prediction by Eq. (47) tends to zero.

Conclusions
The main dynamical parameters for the finite sliding of two identical spheres subjected to displacement and force controlled motions have been considered in the presented analysis. For the linear sliding path of the contacting sphere under displacement control, the formulae relating the horizontal velocity to the sliding path coordinate of the upper sphere center and next the coefficient of restitution were derived analytically in terms of elliptical integrals of the first and second kinds. For the contact time prediction, an approximation via five point Gauss quadrature was also proposed allowing for a high accuracy of the analysis.
For the curvilinear sliding path, the analytical prediction of these parameters cannot be presented in elementary terms or reduced into elliptical integrals of rational functions. Hence, all these parameters have been specified by using the Gauss quadrature resulting in a low error of approximation as well.
For the case of force control, applied in the dynamic analysis of a finite sliding of rigid spheres, the formula for the tangential velocity was derived analytically. The contact time was resolved approximately obtaining a high accuracy of prediction. Meanwhile, in the case of a finite sliding of the deformable sphere, the nonintegrable functions describing the deformation of the sphere in terms of the overlap occur and the numerical technique is only feasible for the dynamical analysis.
The effect of unloading/reloading was also taken into account relying on the proposed simple power law relationships specifying the force-slip displacement curves. In particular, the nonincremental form of this relation may result in spurious oscillations, when the contacting sphere slows down or has multiple contacts. The incremental form proposed allows for elimination of these drawbacks.
The practical aspects in application of the derived dynamical parameters are suitable in modeling of the following: (a) asperity interaction and rough surface contact mechanics; (b) wear analysis; (c) flash temperature specification; (d) proper time step selection in DEM analysis and improvement of the incremental numerical procedure.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.