Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 15, N o 2, 2017, pp. 269 - 284 DOI: 10.22190/FUME170420006W © 2017 by University of Niš, Serbia | Creative Commons Licence: CC BY-NC-ND Original scientific paper THE INFLUENCE OF VISCOELASTICITY ON VELOCITY- DEPENDENT RESTITUTIONS IN THE OBLIQUE IMPACT OF SPHERES UDC 539.3 Emanuel Willert 1 , Stephan Kusche 1 , Valentin L. Popov 1,2,3 1 Berlin University of Technology, Berlin, Germany 2 National Research Tomsk State University, Tomsk, Russia 3 National Research Tomsk Polytechnic University, Tomsk, Russia Abstract. We analyse the oblique impact of linear-viscoelastic spheres by numerical models based on the Method of Dimensionality Reduction and the Boundary Element Method. Thereby we assume quasi-stationarity, the validity of the half-space hypothesis, short impact times and Amontons-Coulomb friction with a constant coefficient for both static and kinetic friction. As under these assumptions both methods are equivalent, their results differ only within the margin of a numerical error. The solution of the impact problem written in proper dimensionless variables will only depend on the two parameters necessary to describe the elastic problem and a sufficient set of variables to describe the influence of viscoelastic material behaviour; in the case of a standard solid this corresponds to two additional variables. The full solution of the impact problem is finally determined by comprehensive parameter studies and partly approximated by simple analytic expressions. Key Words: Oblique Impacts, Friction, Viscoelasticity, Standard Solid Model, Method of Dimensionality Reduction, Boundary Element Method 1. INTRODUCTION Collisions of macroscopic particles determine the dynamics of granular gases. As long as the particle density in the granular gas is small enough and hence the impact durations are small compared to the mean free time between two collisions, these will in general be binary. In many cases the difference of the particle velocities before and after the impact Received April 20, 2017 / Accepted June 21, 2017 Corresponding author: Willert, Emanuel Affiliation: Berlin University of Technology, Sekr. C8-4, Straße des 17. Juni 135, D-10623 Berlin E-mail: e.willert@tu-berlin.de 270 E. WILLERT, S. KUSCHE, V. L. POPOV can be described by two coefficients of restitution, one for each the normal and tangential direction of the impact. Due to friction, adhesion, viscoelasticity, plasticity or other effects those coefficients of restitution will in general exhibit strong and non-trivial dependencies not only of the geometric or material parameters but of the impact velocities themselves. Among the vast literature about granular media only few publication lines account for this velocity-dependence, which is mostly because of two reasons: on the one hand, the various analytical methods of statistical physics applied to deal with granular media are severely complicated by the fact that the restitution coefficients are actually velocity-dependent. On the other hand, the rigorous solution of the single contact-impact problem even in the simplest case of spherical colliding particles is a rather non-trivial task. Lun and Savage [1] and Walton and Braun [2] were the first to study the effects of the described velocity-dependence on the granular dynamics using the granular-flow kinetic theory of Lun, Savage, Jeffrey and Chepurnity. However, lacking rigorous solutions, they only used an ad-hoc model of a restitution coefficient in normal direction exponentially decreasing with the impact velocity, which can be realistic only in few cases. Besides, they did not account for inter-particle friction during the collisions and could hence achieve only rough agreement with their experimental data. Only ten years later a research group around Brilliantov and Pöschel started a series of publications to tackle this problem again. Brilliantov et al. [3] gave models for the collisions of spheres accounting for viscoelasticity and friction. However, their material model is equivalent to a Kelvin- Voigt body, which is only realistic if the time scale of interest is large compared to the relaxation time of the elastomer. As the impact times are short, this might be problematic. Moreover, their tribological friction model of broken welds and asperities leads to a stepwise linear dependence of the tangential force on the tangential displacement between the contacting bodies. For spherical profiles this cannot be true due to the profile shape. These collision models have been implemented in granular gas simulations by Schwager and Pöschel [4], Brilliantov and Pöschel [5] and Dubey et al. [6]. The history of rigorous impact solutions started with Hertz [7], who solved the frictionless and non-adhesive normal contact problem of two parabolic surfaces and the associated quasi- static impact problem. Hunter [8] studied the influence of the quasi-stationarity and found that the proportion of kinetic energy lost during the impact due to elastic wave propagation is negligible, if the impact velocities are small compared to the speed of sound in the elastic medium. Cattaneo [9] and Mindlin [10] solved the tangential contact problem of two elastically similar spheres in the case of a constant normal force and an increasing tangential force. The circular contact area will consist of an inner circular stick area and an annular region of local slip. The tangential traction distribution in the contact is a superposition of two Hertzian distributions. Their work has been extended by Mindlin and Deresiewicz [11] for various different and by Jäger [12] for arbitrary loading protocols. Based on the results of Mindlin and Deresiewicz, Maw et al. [13] and Barber [14] studied the oblique impact of elastic spheres without adhesion; they found out that the problem written in proper dimensionless variables only depends on two parameters, one describing the elastic and the other (containing a generalized angle of incidence and hence the impact velocities) the frictional properties. Moreover, the authors carried out experiments to validate their calculations. The oblique impact problem of elastic spheres with and without adhesion was also studied by Thornton and Yin [15]. A nice overview of elastic impact problems and several analytical solutions including torsional loading can be found in the paper by Jäger [16]. The Influence of Viscoelasticity on Velocity-Dependent Restitutions in the Oblique Impact of Spheres 271 In a series of publications – see for example [17, 18] and the summarizing book [19] – Popov and his co-workers have shown that the generalized Hertz-Mindlin problem for any convex axisymmetric indenter and arbitrary loading histories can be exactly mapped onto a contact between a properly chosen plain profile and a one-dimensional foundation of independent linear springs in such a way that the solution of the obtained one-dimensional model will exactly coincide with the one of the original three-dimensional problem. Due to the enormous simplification and effort reduction of analytical or numerical calculations achieved by this so-called Method of Dimensionality Reduction (MDR) Lyashenko and Popov [20] were able to give a comprehensive solution for the problem studied earlier by Maw and his co-workers in the no-slip regime, i.e. an infinite coefficient of friction. Those results have later been generalized by Willert and Popov [21] for the partial slip regime, i.e. a finite friction coefficient. The viscoelastic contact problem was first addressed by Lee and Radok [22-24]. From the close relationship between the fundamental equations of elasticity and viscosity the authors deduced a method of functional equations to obtain the solution of a viscoelastic problem if the solution of the associated elastic problem is known and the contact radius is a monotonically increasing function in time. This has been generalized to the case of any number of maxima and minima of the contact radius by Graham [25], [26] and Ting [27, 28]. An equivalent but somewhat easier formulation of Ting’s solution was given by Greenwood [29]. However, with every maximum or minimum of the contact radius the analytic calculations get more and more cumbersome. The Hertz impact problem for viscoelastic media was treated by Pao [30] and Hunter [31]. They used arbitrary viscoelastic rheologies to formulate the problem but gave only few concrete solutions. Argatov [32] found analytical solutions for the respective flat punch problem in the case of Kelvin-Voigt-, Maxwell- or standard solid model. The viscoelastic contact problem in the case of convex axisymmetric indenters and arbitrary loading protocols can also be exactly mapped within the framework of the MDR, which was proven by Kürschner and Filippov [33] and Argatov and Popov [34]. Hence, the aim of the present paper is to give a comprehensive solution of the viscoelastic oblique impact of spheres with and without slip based on the MDR. Very recently Kusche [35, 36] presented the no-slip solution of this impact problem using the Boundary Element Method (BEM). However, the BEM-calculations are numerically much more costly compared with the MDR. As the parameter space for the more general case with slip is larger by one dimension, the comprehensive solution based on BEM will be numerically very expensive. Nevertheless, the BEM-algorithm to solve the impact problem with slip has been implemented and can serve as a validation for the faster MDR-based model. We will use a standard solid for modelling viscoelastic properties because it exhibits all characteristics of general elastomers. As a limiting case the Kelvin-Voigt solid is also studied at some point. Finally, we will focus on the velocity-dependence of the coefficients of restitution as this is the main point of interest for the implementation of the obtained solutions into simulation algorithms for granular media. The paper is organized as follows: In Section 2 we will give a formulation of the studied problem. Section 3 is devoted to the description of the numerical model based on the MDR, the results of which are given in Section 5. Section 4 will present a BEM-based algorithm to solve the impact problem, which was used to validate the MDR model described before. Section 6 will give conclusions. 272 E. WILLERT, S. KUSCHE, V. L. POPOV 2. PROBLEM FORMULATION The present paper is concerned with the oblique impact of two linear-viscoelastic spheres of similar materials. This problem is equivalent to the one of a rigid sphere impacting on a viscoelastic half-space, which is why we will restrict ourselves to the latter one. During contact the frictional interaction between the two surfaces shall be assumed to obey the Amontons-Coulomb’s law with the static and the kinetic coefficients of friction being constant and equal to each other: μS = μF ≡ μ. The sphere shall have initial velocities vx0 and vz0, z pointing into the half-space, and initial angular velocity ω0. The mass, radius and moment of inertia of the sphere are m, R and J S , respectively. The point on the sphere which first comes into contact shall be denoted as K. The half-space shall possess a constant Poisson number ν and a creep function giving the response in shear. Actually a viscoelastic material may possess a second creep function for the response to hydrostatic stress, but this shall be neglected. As most elastomers can be considered incompressible (this will also fulfil the condition of elastic similarity) our assumption does not pose a considerable loss of generality. In this case we can introduce time-dependent shear modulus G(t). For the standard solid model G reads: 2 1 2 ( ) exp . G t G t G G          (1) The Kelvin-Voigt model can be recovered from this expression via the limit 2 1 ( ) lim ( ) ( ), KV G G t G t G t     (2) with the Dirac δ-distribution. A scheme of the impact with notations is shown in Fig. 1. We will make further following assumptions: Quasi-stationarity: The impact velocities shall be much smaller than the speed of sound in the viscoelastic material. We therefore neglect all inertia effects like wave propagation. Half-space hypothesis: The surface gradients shall be small. For an axisymmetric contact with parabolic indenter shapes in the vicinity of the contact point, this can be written as max max ,d a R (3) with the maximum values of indentation depth d and contact radius a. Very short impact: The displacement of the contact point due to the change of position and the rotation of the sphere shall be small compared to the contact radius. This ensures that the contact configuration stays axisymmetric and the contact problem can be treated like a tangential one. Rolling will then be accounted for only kinematically. The displacement in vertical direction is of the order of magnitude of the maximum indentation depth. The displacement in tangential direction is of the order of magnitude 0 0 , max 0 .x x K z v R u d v    (4) Hence, this assumption will be covered by the half space hypothesis if the ratio of tangential and vertical initial velocity of the contact point is of the order of 1 or smaller. The Influence of Viscoelasticity on Velocity-Dependent Restitutions in the Oblique Impact of Spheres 273 Fig. 1 Scheme of the analysed impact problem – a rigid sphere is impacting on a viscoelastic half-space 3. NUMERICAL MODEL BASED ON THE MDR Under the assumptions made, the motion of point K fully determines the motion of the sphere. The normal and tangential displacements of this point shall be uK,z and uK,x. The equations of motion for those displacements are elementary given by 2 K , K , 1 , , x x S z z F mR u m J F u m         (5) where Fx and Fz are the contact forces while the dots denote the time derivative. To determine these forces and thereby solve the axisymmetric problem described above within the framework of the MDR, two preliminary steps are necessary. First an equivalent plain profile g(x) has to be obtained from axisymmetric indenter profile f(r) via the Abel-like integral transform 2 2 0 d d ( ) . d x f r g x x r x r    (6) A spherical indenter in the vicinity of the contact can be described by the parabolic profile 2 ( ) 2 r f r R  (7) and the equivalent profile accordingly is given by the expression 2 ( ) . x g x R  (8) 274 E. WILLERT, S. KUSCHE, V. L. POPOV Fig. 2 Single element to model a standard solid Fig. 3 Single element to model a Kelvin-Voigt solid As the second step the viscoelastic properties of the half-space must be replaced by a one-dimensional foundation of independent, linear-viscoelastic elements. In case of a linear standard solid with the time-dependent shear modulus given in Eq. (1) those elements consist of a spring in series with a dashpot, the pair in parallel with a second spring (see Fig. 2). In case of a Kelvin-Voigt model (see Fig. 3) the spring in series with the dashpot is rigid. The elements are at a distance Δx of each other. This value is arbitrary if small enough. Let us first consider the standard solid and write down the necessary relations of the model and the numerical algorithm. All equations for the Kelvin-Voigt model can be derived afterwards by the limiting process. The reaction force for a single element at position xi = i Δx, with outer and inner displacement vectors i u and i u has the components , 1 , , , , 1 , , , 4 ( ) , 2 2 ( ) . 1 i x i x i x i x i z i z i z i z x f G u u u x f G u u u                   (9) The inner point must fulfil the equilibrium conditions 2 , , , 2 ,z ,z ,z ( ) 0, ( ) 0. i x i x i x i i i G u u u G u u u         (10) For the time integration we will use the least order explicit Euler integration scheme with constant time step Δt. The current time step number shall be denoted by an upper index j. In the beginning all displacements are set to zero. Then, in each time step, first the normal contact problem must be solved. For the elements in contact the normal displacement is enforced by the motion of K, 1 , , K, , for contact. j j j i z i z z u u u t     (11) The elements not in contact are free of forces, i.e. the left side of Eqs. (9) is zero, and one obtains The Influence of Viscoelasticity on Velocity-Dependent Restitutions in the Oblique Impact of Spheres 275 1 1 , , , 1 ( ), for no contact, ( ) j j j i z i z i z u u u G t           (12) where we introduce relaxation time 2 / G  . An element gets into contact if , j i z u  K, ( ) j z i u g x and leaves contact if , 0 j i z f  . To solve the tangential contact problem the tangential displacements must be calculated. The elements outside the contact area progress according to 1 1 , , , 1 ( ), for no contact. ( ) j j j i x i x i x u u u G t           (13) For the elements in contact one has to distinguish between sticking and slipping elements. For all the sticking elements, the displacement is enforced by the movement of K, 1 , , K, , for sticking contact. j j j i x i x x u u u t     (14) An element in contact is able to stick if the resulting tangential force does not exceed the maximum value given by the Amontons-Coulomb law, i.e. if , ,z , for sticking contact. j j i x i f f  (15) Any element violating this condition will slip. In this case the tangential force is known to be 1 , ,z , sgn( ), for slipping contact. j j j i x i i x f f f    (16) After the total contact forces are calculated by summation over all elements, , , , j j x i x i j j z i z i F f F f       (17) the equations of motion (5) can be solved in each time step. Note that it is impossible that contact is re-established by the viscoelastic creep. For the Kelvin-Voigt model only τ and hence the inner displacements must be set to zero in the equations above. The algorithm was implemented in MATLAB™. Only time steps j and j-1 have to be stored. That is why this algorithm requires only little memory space. Also all operations are elementary, which makes the algorithm very fast (this is also why we are able to use a least order explicit integration scheme without stability problems) and enables us to do comprehensive parameter studies on an ordinary desktop PC (the calculation of a single impact took around one or two seconds on a machine with an Intel i5 processor). 4. NUMERICAL INVESTIGATION USING BEM The results acquired with the MDR have been validated using the Boundary Element Method (BEM). The BEM-solution of the described problem is numerically exact under the assumptions stated before: the half-space approximation, quasi-static conditions and elastic similarity between the contacting surfaces. Since the BEM does not rely on axis- symmetry, this assumption is only made to have results comparable with the MDR. 276 E. WILLERT, S. KUSCHE, V. L. POPOV The application of the BEM consists of two steps. Firstly the problem of calculating the deflection field from a given pressure distribution and vice versa must be solved. This can be done by utilizing the fundamental solution for a point load acting on a viscoelastic half-space [37-39]. The material is assumed to be incompressible and components Fx, Fy, Fz of the point load are applied at time zero and are kept constant. The deflection of the surface can then be written as 2 2 2 3 3 , 2 1 ( 4 4 , . ), z x x x y z i i J t r x x xy u u r r rr r F y                        (18) In difference to the MDR, the time-dependent creep function for shear J(t) has been used. It is clear and known that J(t) and G(t) are not independent of each other. The creep function can be written by using the constants introduced in equation (1) in the following form: 2 1 2 1 2 1 1 2 1 ( ) 1 1 exp . ( ) G G G J t t G G G G G                  (19) Since the geometric dependences in the viscoelastic and elastic cases are the same, the developed algorithm can be used with only small modifications. In elastic contact mechanics it is a standard procedure to integrate the fundamental solution over a rectangle, assuming constant pressure [40]. This analytic solution is used to find the deflection field for an arbitrary but piecewise constant pressure distribution [41]. This task can be performed very fast and efficiently by using convolution techniques on a parallel computing architecture [42- 44]. The corresponding inverse problem, namely finding the pressure distribution to a given deflection field can be tackled by using the biconjugate gradient stabilized method [45]. The above described methods have been applied to the viscoelastic problem. Since the pressure distribution will change in time, a discretisation is necessary. If, for each time step, the pressure distribution is assumed to be constant, the overall solution in the deflection field can be obtained by adding two solutions in each time step: one to remove the prior load and one to add the current load. Based on the fact that the arising sum grows linearly in time, it is crucial to reduce the numerical effort. This can be achieved by applying the special form of the creep function (19) and by observing the following time step. Then an iterative algorithm can be developed: 1 1 1 , 1 1 0 00 0 0 , 1 1 1 0 0 0 1 1 1 ( ) ( )exp exp ( ) z z j j z z j j j j j i z j n i i i i i i i a b z z z j j j j j j j a b j j t z t tJ J f t t f f f J J J J J f a f b f f J J J f a u t t u                                                    0 exp z zt j j D h J b f J         (20) Herein uz,n is the normal deflection of the surface, J∞ = J(t = ∞), J0 = J(t = 0), and fj is the deflection due to a pressure distribution pj – each at the time tj. In the last line of Eq. (20) The Influence of Viscoelasticity on Velocity-Dependent Restitutions in the Oblique Impact of Spheres 277 it can be seen that an additional deformation D z has to been taken into account to include viscoelastic behaviour (in the elastic case D z is equal to zero). The unknown term in the last line of Eq. (20) is fj+1, which means that the pressure distribution pj+1 is unknown. This can be calculated with the elastic algorithms mentioned before. It should be noted that this algorithm can handle only materials with a finite modulus of instant deformation, which excludes the Kelvin-Voigt solid. The tangential contact can be solved very similarly to the normal contact so that the same scheme can be used [46]. Only the calculation of the deflection in tangential direction ux, caused by shear stress has to be adopted. If a partial slip is involved, the calculation is modified in the following way: starting with a complete stick area, the deflection is given by the increment of displacement in one time step. If this leads to shear stress that is larger than the value allowed by the Coulomb’s law, this part of the contact area will slip. In the slip areas the tangential stress is set to |τ| = μp. Then the stress in the remaining stick area is calculated again, under the consideration of the deflection caused by the shear stress in the slip area. This is done until the stick area does not change anymore. In all performed simulations, the deformation perpendicular to the plane of the motion, uy , is neglected. It turns out that this assumption, in the case of parabolic bodies, causes a negligible error [47]. At this point, the contact problem itself is solved. For the integration in time both an explicit Euler scheme and the velocity Verlet algorithm have been used. In comparison, they show no difference in the global error of the velocities at the end of the simulation and in the contact time itself. For an estimation of the step size Δt the MDR solution has been used. For the geometric discretization a matrix of 256256 points has been chosen. The comparison with a finer discretization shows only a slight error reduction. For implementation it has to be considered that the total deflection in normal direction within the contact area is known at every time step since the indentation depth of the sphere is known. Contrariwise in tangential direction: the points coming into contact have a pre-deformation through coupling to the points within the contact area from a previous time step. This can be handled by adding only the current increment of tangential movement at the boundary of the sphere in each time step. The systematic investigation of the problem has been done with the MDR. The processing time for the BEM is much higher compared to the MDR. Therefore, only a few hundred parameter sets spreading over the full range covered by the investigation done with the MDR have been calculated with the BEM. It turns out that the relative differences in the coefficients of restitution have always been smaller than 0.5%. Therefore, it is reasonable that the MDR can be applied. 5. RESULTS OF THE NUMERICAL MODEL: THE RESTITUTION COEFFICIENTS As a solution we are interested in the coefficients of restitution in normal and tangential direction , , 0 , , 0 0 ( ) , ( 0) ( ) . ( 0) K z e z z K z z K x e x x K x x u t t v e u t v u t t v R e u t v R                 (21) 278 E. WILLERT, S. KUSCHE, V. L. POPOV Maw et al. [13] have shown that in an ideally elastic case (the coefficient in normal direction being obviously unity) the coefficient in tangential direction only depends on the two dimensionless parameters 2 0 0 0 1 2 2 1 , . 2 2 x S z v RmR vJ                  (22) In the case of a sphere impacting on a viscoelastic half space modelled as a linear standard solid, two more dimensionless parameters are of interest, describing the viscoelastic material properties, namely el,1 1 1 1 2 , , (1 ) a G G M G        (23) with the maximum contact radius for the impact with an elastic half space, 1/ 5 2 2 0 el, 15 (1 ) , 1, 2 32 z i i mv R a i G        (24) Of course, any combinations of those two additional parameters would also be possible to choose as governing variables. For example, in the previous publication on the no-slip impact Kusche [35] used the parameters el,1 el,2 3 / 5 1 2 1 1 2 , (1 ) (1 ) a a G M G M              (25) to capture the influence of the material behaviour. However, as we are interested mainly in the velocity-dependence of the coefficients of restitution, it seems convenient to select δ1 and γ, because the latter one is velocity-independent and therefore the velocity- dependence due to viscoelasticity can be fully covered by parameter δ1. Moreover, the Kelvin-Voigt model can be recovered as the limiting case γ = 0. Also limit γ → ∞ corresponds to the elastic result. To reduce the number of governing parameters, we restrict ourselves mostly to χ = 7/6, which, amongst other cases, corresponds to the case of incompressible, homogenous spheres. To prove that actually 1 1 ( , ) and ( , , , ) z z x x e e e e       (26) we made comprehensive numerical studies, the results of which are shown in the upcoming figures. Thereby we first focus on the limiting case of a Kelvin-Voigt solid and afterwards look at the more general standard solid. In Fig. 4 the coefficient of restitution in normal direction is shown for a Kelvin-Voigt solid as a function of δ1. All free input parameters for the simulations, i.e. velocities, measures of inertia and so on, have been generated randomly. Nevertheless, the points create continuous curves and hence our hypothesis is proven for the normal direction. It is easy to interpret the results, as the coefficient of normal restitution shows the often-used quasi-exponentially decreasing behaviour. This, however, only remains true for this material model of a Kelvin-Voigt solid, which corresponds to an infinitely fast relaxation within The Influence of Viscoelasticity on Velocity-Dependent Restitutions in the Oblique Impact of Spheres 279 the elastomer. It was already pointed out that this is problematic as the impact times are considered to be small and the relaxation time has to be accounted for in some way. We will see the effects later in the results for the standard solid. Fig. 4 Coefficient of restitution in normal direction for the impact on a Kelvin- Voigt solid as a function of δ1 Fig. 5 Coefficient of restitution in tangential direction as a function of δ1 and ψ with χ = 7/6. Online version in colour Fig. 5 gives the tangential coefficient of restitution ex as a function of δ1 and ψ for the impact on a Kelvin-Voigt half-space. The value of χ was fixed at 7/6, all other input parameters for the impact problem have been generated randomly and yet the solutions create continuous, smooth curves. The tangential restitution has a global maximum for an impact without viscosity around ψ = 2. On the right side of the contour plot the behaviour gets quite simple and can be explained the following way: for any material model configurations are possible for which the contact will completely slip during the whole impact. In this case the total tangential force is known due to the Coulomb’s law and hence the tangential restitution coefficient for full slip (and any material model) is given by the relation 2 (1 ) 1. fs x z e e      (27) Let us now look into the results for the standard solid. We restrict ourselves again to the case χ = 7/6 to spare the generally least important parameter. Fig. 6 gives the results for the normal restitution coefficient as a function of δ1. Several logarithmically-equally-distributed values for γ have been chosen and all other input parameters for the impact problem, as always, have been generated randomly. Nevertheless, the solutions create continuous curves and it is easy to observe the influence of γ on the velocity-dependent restitution: as said before γ → ∞ corresponds to the trivial elastic case and γ = 0 to the monotonically decreasing Kelvin-Voigt solution. For intermediate values of γ the coefficient of restitution has a global minimum. After that it increases again with increasing δ1, i.e. increasing normal inbound velocities. This distinguishes the general standard solid from its limiting case with infinitely fast relaxation and has, for example, a very interesting consequence for a (driven) granular gas of viscoelastic particles: as on the increasing part of the restitution curve, the coefficient 280 E. WILLERT, S. KUSCHE, V. L. POPOV of restitution is larger for larger inbound velocities, a region of locally higher internal energy of the granular gas, i.e. higher velocities of the particles, might dissipate less energy than regions of lower energy, which might result in unstable states of the granular gas. Fig. 6 Coefficient of restitution in normal direction as a function of δ1 (logarithmic) and different values of γ for the impact with a standard solid Fig. 7 Coefficient of restitution in tangential direction without slip (ψ = 0) as a function of δ1 (logarithmic) and different values of γ for the impact with a standard solid; χ = 7/6 Fig. 8 Coefficient of restitution in tangential direction as a function of δ1 (logarithmic) and ψ for the impact with a standard solid; χ = 7/6 and γ = 0.0825 Fig. 9 Coefficient of restitution in tangential direction as a function of δ1 (logarithmic) and ψ for the impact with a standard solid; χ = 7/6 and γ = 1 Fig. 7 presents the results for the tangential restitution in the case of no slip, which have been reported by Kusche [35] with a slightly different set of governing dimensionless parameters. In Fig. 8 and 9 the results are shown for the behaviour with slip. For increasing values of γ a bulb with ex ≈ 0.5 around ψ ≈ 2 is stretching to the left, i.e. the area with less viscosity. The other areas are less affected by the material properties. Finally, we come back to the full slip solution and the different regimes for parameter ψ. The Influence of Viscoelasticity on Velocity-Dependent Restitutions in the Oblique Impact of Spheres 281 In the elastic case Maw et al. [13] distinguish three different regimes: for ψ < 1 the impact will start in a completely sticking contact and remain like this during the whole compression phase; for ψ > 4χ – 1 the contact will fully slip during the whole impact; the intermediate values correspond to a mixed regime. Now, in the viscoelastic case, the time derivatives of the contact forces in the MDR-model in the very first moment of contact are given by , , x 2 ( 0) ( 0) ( 0) 1 4 ( 0) ( 0) ( 0) 2 z K z x K x F t G t u t x F t G t u t               (28) Hence, for a finite instantaneous stiffness (this excludes the Kelvin-Voigt body), the impact will begin with sticking contact, if ( 0) ( 0) 1. x z F t F t      (29) In case of the Kelvin-Voigt body the contact forces in the first moment of contact are nonzero and the no-slip condition will be ( 0) ( 0) 1. x z F t F t      (30) Hence, this lower transition value for ψ is unaffected by viscoelasticity. For any standard solid characterised by the two parameters δ1 and γ – and probably any material behaviour – there also exists a value ψc, for which the contact will completely slip during the whole impact if ψ > ψc. For complete slip the tangential coefficient of restitution is given by Eq. (27). Fig. 10 Critical value ψc, for which the contact will completely slip during the whole impact if ψ > ψc for the impact with a standard solid Fig. 11 Relative error between the numerical result for ψc and the analytic approximation (31) In Fig. 10 the value of c  is shown for different materials. Obviously this transition value strongly correlates with the normal restitution coefficient. The global maximum is elastic case 282 E. WILLERT, S. KUSCHE, V. L. POPOV ψc = 11/3, and a very good approximation (with a relative error always smaller than 0.2%, see Fig. 11) is given by the expression 2 (1 ). c z z e e     (31) 6. CONCLUSIONS Based on the numerical models we investigated the oblique impact of linear-viscoelastic spheres under the assumptions of quasi-stationarity, the validity of the half-space hypothesis, Amontons-Coulomb friction and short impact times. Numerical models based on both the Method of Dimensionality Reduction (MDR) and the Boundary Element Method (BEM) have been implemented. As expected both methods in their results only differ within the margin of a numerical error. Due to the enormous reduction of mathematical and computational effort achieved by the MDR we were able to perform comprehensive parameter studies for the examined impact problem. It is found that the problem solution, i.e. the coefficients of normal and tangential restitution, written in proper dimensionless variables will depend on exactly four different values, at least two of which contain explicit dependencies on the inbound velocities. By accounting for the finite relaxation time within the elastomer it is possible to increase the normal restitution coefficient with increasing inbound velocities. This is in contrast with most viscoelastic collision models used in the literature about granular media and may have interesting applications in granular chains or gases. As in the elastic case, three different regimes are possible depending on the inbound velocities: the contact may fully slip during the whole impact, completely stick during the compression phase or be in a mixed regime. Viscoelasticity reduces the angle of incidence necessary to ensure complete slip but does not affect the transition between the two other regimes. The transition to full slip strongly correlates with the coefficient of normal restitution. Of course, in practice the here-given assumptions pose severe restrictions, especially the half-space hypothesis, the assumed short impact time and the assumption of perfectly linear material behaviour. Nevertheless, the proposed model and its solution to the best our knowledge are the first – from a contact-mechanical point of view – rigorous and self- consistent approach to the topic despite the extensive existing literature dealing with it. The proposed methods can without problems be applied to more general forms of the time-dependent shear modulus, for example represented in a Prony series. REFERENCES 1. Lun, C.K.K., Savage, S.B., 1986, The Effects of an Impact Velocity Dependent Coefficient of Restitution on Stresses Developed by Sheared Granular Materials, Acta Mechanica, 63(1), pp. 15–44. 2. Walton, O.R., Braun, R.L., 1986, Stress Calculations for Assemblies of Inelastic Spheres in Uniform Shear, Acta Mechanica, 63(1), pp. 73–86. 3. Brilliantov, N.V., Spahn, F., Hertzsch, J.M., Pöschel, T., 1996, Model for Collisions in Granular Gases, Physical Review E, 53(5), pp. 5382–5392. 4. Schwager, T., Pöschel, T., 1998, Coefficient of Normal Restitution of Viscous Particles and Cooling Rate of Granular Gases, Physical Review E, 57(1), pp. 650–654. 5. Brilliantov, N.V., Pöschel, T., 2000, Velocity Distribution in Granular Gases of Viscoelastic Particles, Physical Review E, 61(5B), pp. 5573–5587. The Influence of Viscoelasticity on Velocity-Dependent Restitutions in the Oblique Impact of Spheres 283 6. Dubey, A.K., Brodova, A., Puri, S., Brilliantov, N.V., 2013, Velocity Distribution Function and Effective Restitution Coefficient for a Granular Gas of Viscoelastic Particles, Physical Review E, 87, 062202. 7. Hertz, H., 1882, Über die Berührung fester elastischer Körper, Journal für die reine und angewandte Mathematik, 92, pp. 156–171. 8. Hunter, S.C., 1957, Energy Absorbed by Elastic Waves during Impact, Journal of the Mechanics and Physics of Solids, 5(3), pp. 162–171. 9. Cattaneo, C., 1938, Sul Contato di Due Corpo Elastici, Accademia dei Lincei, Rendiconti, Series 6, 27, pp. 342–348, 434–436 and 474–478. 10. Mindlin, R.D., 1949, Compliance of Elastic Bodies in Contact, Journal of Applied Mechanics, 16, pp. 259–268. 11. Mindlin, R.D., Deresiewicz, H., 1953, Elastic Spheres in Contact under Varying Oblique Forces, Journal of Applied Mechanics, 20, pp. 327–344. 12. Jäger, J., 1993, Elastic contact of Equal Spheres under Oblique Forces, Archive of Applied Mechanics, 63(6), pp. 402–412. 13. Maw, N., Barber, J.R., Fawcett, J.N., 1976, The Oblique Impact of Elastic Spheres, Wear, 38(1), pp. 101–114. 14. Barber, J.R., 1979, Adhesive Contact during the Oblique Impact of Elastic Spheres, Journal of Applied Mathematics and Physics (ZAMP), 30, pp. 468–476. 15. Thornton, C., Yin, K.K., 1991, Impact of Elastic Spheres with and without Adhesion, Powder Technology, 65(1- 3), pp. 153–166. 16. Jäger, J., 1994, Analytical Solutions of Contact Impact Problems, Applied Mechanics Review, 47(2), pp. 35–54. 17. Popov, V.L., Heß, M., 2014, Method of Dimensionality Reduction in Contact Mechanics and Friction: A Users Handbook. I. Axially Symmetric Contacts, Facta Universitatis, Series Mechanical Engineering, 12(1), pp. 1–14. 18. Popov, V.L., Pohrt, R., Heß, M., 2016, General Procedure for Solution of Contact Problems under Dynamic Normal and Tangential Loading Based on the Known Solution of Normal Contact Problem, Journal of Strain Analysis for Engineering Design, 51(4), pp. 247–255. 19. Popov, V.L., Heß, M., 2015, Method of Dimensionality Reduction in Contact Mechanics and Friction, Springer, Heidelberg, ISBN 978-3-642-53875-9. 20. Lyashenko, I.A., Popov, V.L., 2015, Impact of an Elastic Sphere with an Elastic Half Space Revisited: Numerical Analysis based on the Method of Dimensionality Reduction, Scientific Reports, 5, 8479. 21. Willert, E., Popov, V.L., 2016, Impact of an Elastic Sphere with an Elastic Half Space with a Constant Coefficient of Friction: Numerical Analysis Based on the Method of Dimensionality Reduction, ZAMM Zeitschrift für Angewandte Mathematik und Mechanik, 96(9), pp. 1089–1095. 22. Lee, E.H., 1955, Stress Analysis in Visco-Elastic Bodies, Quarterly Applied Mathematics, 13(2), pp. 183–190. 23. Radok, J.R.M., 1957, Visco-Elastic Stress Analysis, Quarterly Applied Mathematics, 15(2), pp. 198–202. 24. Lee, E.H., Radok, J.R.M., 1960, The Contact Problem for Viscoelastic Bodies, Journal of Applied Mechanics, 27(3), pp. 438–444. 25. Graham, G.A.C., 1965, The Contact Problem in the Linear Theory of Viscoelasticity, International Journal of Engineering Science, 3(1), pp. 27–46. 26. Graham, G.A.C., 1967, The Contact Problem in the Linear Theory of Viscoelasticity When the Time-Dependent Contact Area Has any Number of Maxima and Minima, International Journal of Engineering Science, 5(6), pp. 495–514. 27. Ting, T.C.T., 1966, The Contact Stresses between a Rigid Indenter and a Viscoelastic Half Space, Journal of Applied Mechanics, 33(4), pp. 845–854. 28. Ting, T.C.T., 1968, Contact Problems in the Linear Theory of Viscoelasticity, Journal of Applied Mechanics, 35(2), pp. 248–254. 29. Greenwood, J.A., 2010, Contact between an Axisymmetric Indenter and a Viscoelastic Half Space, International Journal of Mechanical Sciences, 52(6), pp. 829–835. 30. Pao, Y.H., 1955, Extension of the Hertz Theory of Impact to the Viscoelastic Case, Journal of Applied Physics, 26(9), pp. 1083–1088. 31. Hunter, S.C., 1960, The Hertz Problem for a Rigid Spherical Indenter and a Viscoelastic Half Space, Journal of the Mechanics and Physics of Solids, 8(4), pp. 219–234. 32. Argatov, I.I., 2013, Mathematical Modeling of Linear Viscoelastic Impact: Application to Drop Impact Testing of Articular Cartilage, Tribology International, 63, pp. 213–225. 33. Kürschner, S., Filippov, A.E., 2012, Normal Contact between a Rigid Surface and a Viscous Body: Verification of the Method of Reduction of Dimensionality for Viscous Media, Physical Mesomechanics, 15(4), pp. 25–30. 34. Argatov, I.I., Popov, V.L., 2015, Rebound Indentation Problem for a Viscoelastic Half-Space and Axisymmetric Indenter – Solution by the Method of Dimensionality Reduction, ZAMM Zeitschrift für Angewandte Mathematik und Mechanik, 96(8), pp. 956–967. 284 E. WILLERT, S. KUSCHE, V. L. POPOV 35. Kusche, S., 2016, The Boundary Element Method for Visco-elastic Material Applied to the Oblique Impact of Spheres, Facta Universitatis, Series Mechanical Engineering, 14(3), pp. 293–300. 36. Kusche, S., 2016, Simulation von Kontaktproblemen bei linearem viskoelastischem Materialverhalten, Doctoral Dissertation, Technische Universität Berlin 37. Talybly, L., 2010, Boussinesq’s viscoelastic problem on normal concentrated force on a half-space surface, Mechanics of Time-Dependent Materials, 14(3), pp. 253–259. 38. Gasanova, L., Gasanova, P., Talybly, L., 2011, Solution of a viscoelastic boundary-value problem on the action of a concentrated force in an infinite plane, Mechanics of Solids, 46(5), pp. 772–778. 39. Peng, Y., Zhou, D., 2012, Stress Distributions Due to a Concentrated Force on Viscoelastic Half-Space, Journal of Computation & Modeling, 2(4), pp. 51–74. 40. Johnson, K. L., 1985, Contact Mechanics, Cambridge University Press, Cambridge. 41. Pohrt, R., Li, Q., 2014, Complete Boundary Element Formulation for Normal and Tangential Contact Problems, Physical Mesomechanics, 17(4), pp. 334-340. 42. Cho, Y. J., Koo, Y. P., Kim, T. W., 2000, A new FFT technique for the analysis of contact pressure and subsurface stress in a semi-infinite solid, KSME International Journal, 14(3), pp. 331–337. 43. Liu, S., Wang, Q., Liu, G., 2000, A versatile method of discrete convolution and FFT DC-FFT for contact analyses, Wear, 243(1-2), pp. 101–111. 44. Wang, W.Z., Wang, H., Liu, Y.C., Hu, Y.Z., Zhu, D., 2003, A comparative study of the methods for calculation of surface elastic deformation, Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, 217, pp. 145–154. 45. van der Vorst, H.A., 1992, BI-CGSTAB: A Fast and Smoothly Converging Variant of BI-CG for the Solution of Nonsymmetric Linear Systems, SIAM Journal of Scientific and Statistical Computing, 13(2), pp. 631–644. 46. Kusche, S., 2016, Frictional force between a rotationally symmetric indenter and a viscoelastic half-space, ZAMM Zeitschrift für Angewandte Mathematik und Mechanik, DOI: 10.1002/zamm.201500169. 47. Munisamy, R.L., Hills, D.A., Nowell, D., 1994, Static axisymmetric Hertzian contacts subject to shearing forces, Journal of Applied Mechanics, 61(2), pp. 278–283.