Viscoelastic Hertzian Impact

: The problem of normal impact of a rigid sphere on a Maxwell viscoelastic solid half-space is considered. The first-order asymptotic solution is constructed in the framework of Hunter’s model of viscoelastic impact. In particular, simple analytical approximations have been derived for the maximum contact force and the time to achieve it. A linear regression method is suggested for evaluating the instantaneous elastic modulus and the mean relaxation time from a set of experimental data collected for different spherical impactors and impact velocities.


Introduction
Loosely speaking, an impact between two solids is a dynamic mechanical process characterized by some exchange of their kinetic energies via intermediate transformation into the potential energy stored in elastic deformations.If a part of the initial total kinetic energy is lost during the impact, it is called inelastic, and reasons for the energy loss can be different, depending on the mechanical properties of the solids [1,2].A wide special class of materials with intrinsic dissipation of strain energy covers viscoelastic materials whose response to deformation depends on the strain rate [3].
Recent interest in viscoelastic impact arose in modeling such diverse phenomena as dynamic response of biological soft tissues [4], hopping of a rover on an asteroid [5], fruit drop tests [6], impact of nanomanipulation of nanoparticles [7], kinetic energy dissipation in granular matter [8], seismic pounding [9], erosion [10], and protective performance of flexible polymer material [11].It goes without saying that the case of linearly viscoelastic materials is a base case for more elaborate time-dependent material models such as poroviscoelasticity [12] and surface viscoelasticity [13].Additionally, a rigid sphere is one of the most frequently used types of impactor [14,15].
As a time process, impact starts at the moment of initial contact (see Figure 1), and its subsequent time evolution strongly depends on the initial contact configuration, including the initial contact geometry and the initial kinematic conditions.The so-called Hertzian contact geometry assumes a single-point initial contact and the second-degree paraboloid approximation for the initial gap between the contacting surfaces [16].For what follows, two special cases of normal contact should be distinguished among others, namely, (a) impact between two solid spheres (see Figure 1a) and (b) impact between a solid sphere and a semi-infinite solid (see Figure 1b).In both cases, we can trace up a characteristic feature of the Hertzian impact as a continuous variation of the contact area (which is circular due to the axisymmetric contact geometry and the assumption of normal impact) during the impact starting from a point of initial contact to a point of final contact, passing through a single maximum.Whereas, in fully elastic impacts, the loading and unloading stages are symmetric, there is a marked difference between the two stages in viscoelastic impacts, depending on the share of the total kinetic energy dissipated during the impact.The end of an impact is determined by the condition of vanishing the contact reaction.
V 0 In his ground-laying paper [17] on the frictionless local contact and impact of elastic solids, Heinrich Hertz estimated the half-duration, t 0 m = t 0 c /2, of the normal dissipationless impact as follows: where V 0 is the initial relative velocity of the approach, and w 0 m is the maximum value of the contact approach (evaluated from the initial contact moment).For what follows, the superscript 0 is attached to the Hertzian contact parameters, as Hertz's model will be treated as an unperturbed model.It is pertinent to note that the scaling law t 0 m ∼ w 0 m /V 0 follows from simple dimensional considerations, but the main problem is to determine the dimensionless proportionality constant τ 0 m , which is sometimes called Hertz's integral.Hertz's theory of elastic impact can be represented as the initial-value problem for the following second-order nonlinear differential equation: where w is the contact approach measured from the time, t, of initial contact; ẇ and ẅ are the impact velocity and acceleration, respectively; m is the equivalent mass; and F is the contact force (reaction), which, according to Hertz's contact law, is given by the following: with k being the stiffness coefficient.For the basic facts of Hertz's theory that we recall here, we refer to textbooks on theoretical and solid mechanics (see, e.g., [2,16]).
In the case of collision between two elastic spheres (see Figure 1a), the equivalent mass is given by m = m 1 m 2 /(m 1 + m 2 ), whereas the stiffness coefficient k = (4/3)E * √ R is determined in terms of the equivalent radius, R, and the effective elastic modulus, E * , defined as follows: 1 From Equations ( 2) and (3), the equation of energy conservation follows in the following form: which determines the maximum contact approach, as follows: achieved at the time moment t = t 0 m when ẇ = 0.The same Hertzian Equations ( 1)-( 3), (5), and (6) also apply in the case of the impact of a rigid sphere of mass m on an elastic half-space (see Figure 1b), as Equation ( 4) allows a passage to the limit as E 1 → ∞ and R 2 → ∞.We recall that Hertz's theory of frictionless contact assumes that the contact is developed in the framework of the linear theory of elasticity, the elastic bodies are assumed to be isotropic and homogeneous, the contact is local in a sense that the initial contact occurs at a single point only, and the elastic half-space approximation applies for evaluating the contact stresses by neglecting the effect of global contact geometry [18].
It should be noted that some of the Hertz model's restrictions can be relaxed.For instance, in a routine manner, Willis [19] extended Hertz's theory of impact to anisotropic bodies, using his solution for the problem of local frictionless contact.By utilizing Bondareva's solution [20] for a heavy elastic sphere on a rigid plane, Villaggio [21] showed that the global contact geometry effect slightly increases the contact duration, compared with that predicted by the classical Hertz's theory.The inertia effect revealing itself in the impact energy loss due to the elastic wave radiation in the impact problem for an elastic half-space (Figure 1b) was first estimated by Hunter [22] (see also [23,24]) based on the analytical solution obtained by Miller and Pursey [25] for the elastic wave energy radiated by a rigid disk vibrating on the half-space surface.It is pertinent to note here that the excitation of the half-space surface by a spherical impact was considered in [26].
Hertz's theory of impact predicts the unit coefficient of restitution, e; the symmetry of loading/unloading contact process; and the duration of impact, t 0 c , to be equal twice the time t 0 m to the maximum of the contact approach w 0 m .It was shown by Hunter [22] and Deresiewicz [27] that the variation of the contact approach as a function of time during the contact can be well approximated by the following simple formula: The classical Hertz impact theory has been given substantial experimental verification [1] and, in particular, was extended to the frictional impact of anisotropic nonlinear elastic solids [28] and power-graded viscoelastic solids [29] as well as to tangential (oblique) impact [30][31][32] and impact with adhesion [33][34][35].
It is known [36,37] that the coefficient of restitution in the collision of two perfectly elastic bodies is almost equal to the unit (that is, e ≈ 1), if the time of impact well exceeds the time needed for elastic waves to traverse either body.That is why the impact configuration shown in Figure 1b primarily differs from that shown in Figure 1a by the presence of the energy dissipation (absorption [22]) mechanism due to the vibrational energy radiated into the massive substrate when elastic waves propagate to the infinity.
Energy dissipation in the Hertzian impact between two spherical solids (Figure 1a) can be associated with the effects of plastic deformation or internal friction among others [37].A phenomenological approach (see, e.g., [38,39]) leads to the following Hunt-Crossley dissipative contact model [40]: Here, β is a dimensionless constant and χ denotes the hysteresis damping factor.As a rule, the dimensional coefficient χ is interpreted in terms of the coefficient of restitution [41].By adopting a constitutive law similar to the Kelvin-Voigt model σ zr = 2µ 0 ε zr + 2η εzr , where η is the viscosity coefficient, Goldobin et al. [42] arrived at Equation (8) with β = 1/2 and χ being proportional to a linear combination of the shear and bulk viscosity coefficients (see also [43,44]).However, it should be noted that the Kelvin-Voigt model has the relaxation function in terms of Dirac's delta function (see, e.g., [45]), and therefore, it is not appropriate for modeling the material behavior under a blunt impact as it predicts a jump-like contact reaction at the impact moment (see, e.g., [46]).
It should be emphasized that, in contrast with the Hunt-Crossley model (8), where the current contact reaction depends only on the current kinematic variables w(t) and ẇ(t), in the viscoelastic Hertzian contact, a preceding history of loading is important.This allows for gaining insight into the hereditary effect in a bouncing impact [47], whereas the Hunt-Crossley reaction element forgets about the prior loading as soon as it returns into the intact state (i.e., w(t) = 0).
The impact problem becomes exceedingly hard if colliding solids are assumed to possess time-dependent mechanical properties.The Hertz impact problem for a rigid spherical indenter and a viscoelastic half-space was first considered by Hunter [48], who complemented the analytical solution by Lee and Radok [49] for the Hertzian quasi-static contact problem with monotonically increasing contact by the solution when the contact radius possesses a single maximum, which is the case in single impact problems.In the special case of a Maxwell solid, Hunter obtained the first-order perturbation approximations for the coefficient of restitution, e, and the impact duration, t c .Later, Forney [50] questioned Hunter's result about the impact duration, which indirectly casts a shadow on the entirety of Hunter's solution.To date, this issue remains unresolved.
The case of a Kelvin-Voigt solid (with accounting for the unloading stage) was considered by Khusid [51] (see also [29]), who obtained some numerical results for the impact duration and the coefficient of restitution.A systematic review of modeling linear and non-linear viscoelastic contact problems was recently given by Wang et al. [52].
In what follows, we consider the normal impact of a rigid sphere on an isotropic viscoelastic half-space with a constant Poisson's ratio, ν, and a hereditary constitutive law, as follows: where σ zr and ϵ zr are the shear stress and strain, respectively; µ(t) is the shear relaxation modulus; t ′ is the integration variable; and 0 − indicates the instant immediately before the initial point of contact.
In his first-order perturbation analysis of the viscoelastic Hertzian impact, Aksel [53] applied the viscoelastic constitutive law in the following form: where µ 0 = µ(0) is the instantaneous shear modulus, and μ(t) = dµ(t)/dt is the viscoelastic relaxation kernel, and it is tentatively assumed that ϵ zr (t) = 0 for t < 0. In what follows, the constitutive equation in the compact form ( 9) is preferred over (10).Additionally, in the spherical impact problem, the kinematic contact variables start to vary from the zeroth values, and therefore, the lower limit in the hereditary integral ( 9) may be replaced with a zero value.
The practically important problem of material parameter identification by means of impact tests was considered in a number of experimental studies [54][55][56][57].Kren and Naumov [58] formulated the inverse problem of determining the relaxation function µ(t) from the spherical impact loading history (impactor velocity, V(t); contact force, F(t); and contact approach, w(t)) without a priori adopting any material model.However, the problem with the Kren-Naumov method is that the Lee-Radok solution [49], which is valid only for the loading contact stage, was incorrectly applied for the unloading stage as well.It is noteworthy here that the approximation that is utilized for the restitution phase of the same form of the equation of motion derived for compression is sometimes used for the sake of simplicity [59,60].Apparently, this simple but erroneous approach to the viscoelastic Hertzian impact problem dates back to the seventy-year-old paper [61], which was already criticized by Lee and Radok [49].
Why do we need to study the viscoelastic Hertzian impact then?If the material behavior follows a linearly viscoelastic constitutive law, the dissipative contact model (8) still can be applied for describing a spherical impact, but the model parameters χ and β should be initialized from the impact itself (indeed, χ and β are not material parameters).On the other hand, when the impact model has been derived using the solution of the viscoelastic Hertzian contact problem, the viscoelastic material parameters that enter the contact reaction model can be determined from an independent material deformation test.However, the latter problem (that is, the problem of the viscoelastic Hertzian impact) is much more mathematically difficult than the impact problem given by Equations ( 2) and ( 8).Here, we are not looking for simple solutions but rather for the correct one.

Hunter's Model of the Spherical Viscoelastic Impact 2.1. Viscoelastic Hertzian Impact
For the sake of simplicity, we consider the single impactor configuration (see Figure 1b) and start with Newton's second law and the initial conditions: where F is the contact reaction.Assuming a constant Poisson's ratio ν and neglecting inertia effects in the viscoelastic half-space target, the Hertzian elastic solution can be generalized (for the loading stage) as follows [48,49]: Here, R is the sphere radius, µ(t) is the shear relaxation modulus that describes the material's time-dependent response to an instantaneous unit shear deformation, a(t) is the contact radius as a function of time t, and t m is the duration of the loading stage.It is clear that Equation (13) is the Hertzian relation between the contact radius a and the contact approach w.
From Equations ( 11)-( 13), the following holds: It should be remembered that, in the loading stage, the Lee-Radok Equations ( 12) and ( 13) and Equation ( 14) are valid until the time moment t m of maximum penetration, w m , which, in view of Equation ( 13), coincides with the moment of maximum contact radius, a m , and, therefore, with the moment when the impactor velocity vanishes, as follows: In the unloading stage (t m < t < t c , where t c denotes the end time of contact), Hunter's solution is given in terms of the auxiliary function t 1 (t) that solves the equation a(t) = a(t 1 ) for t > t m and t 1 < t m .Namely, the contact reaction is given by the following: whereas Equation ( 13) should be replaced with the following relation [48]: Here, µ −1 (t) is the shear creep compliance (we keep the notation from [48]).

Impact for a Maxwell Solid
The reciprocal relaxation and creep functions for a Maxwell solid are the following: where µ 0 is the instantaneous shear modulus, and η is an inverse relaxation time.
The substitution of (18) 1 into Equation ( 14), in view of the initial conditions (11), leads to the following differential equation: which is an exact result.However, in the unloading (or withdrawal) stage, the application of Equations ( 18) allows for reducing Equation ( 17) to a complicated nonlinear differential equation for t 1 (t) (see [48] for details).Nevertheless, by utilizing the simple approximation t 1 ≈ 2t m − t, Hunter derived the following resulting governing differential equation: which is amendable to analytical treatment, but constitutes an approximate result.
It should be noted that the governing Equations ( 19) and ( 20) make use of slightly different but similar symbols w and w.In the first case, w coincides with the contact approach, which is related to the relative squared contact radius by Equation (13).In the second case, the symbol w preserves the latter geometrical meaning of w, but not its kinematic meaning (i.e., w does not coincide with the contact approach).
In fact, we introduce the following auxiliary function: which is related to the contact approach by the following differential equation: We note that, in view of ( 13) and ( 21), the solution of Equation ( 20) is subject to the following boundary conditions (or the initial conditions for the withdrawal stage): The contact duration t c is determined by the condition w(t c ) = 0, when the contact shrinks down.
In the following dimensionless variables: where w 0 m is given by ( 6) with k being replaced by 8 19) and ( 20) become the following: where τ m is the time-like point for which the solution of Equation ( 25) first yields dω/dτ = 0, and the corresponding value of ω(τ m ) will be denoted by ω m .

Asymptotic Solution for the Loading Stage
The problem with Forney's critique [50] of Hunter's approximate solution is in the following integral decomposition representation: where ω = dω/dτ, because the function ω cannot be regarded as a perturbation of the limit (ε = 0) elastic solution ω 0 (τ) that is defined only on the interval [0, τ 0 m ], where τ 0 m is the dimensionless Hertzian half duration of an impact.At the same time, the second integral on the right-hand side of Equation ( 27) refers to the interval ω ∈ (1, ω m ), which corresponds to τ ∈ (τ 0 m , τ m ); that is, the second term in (27) falls outside the interval of validity of the limit solution.
That is why we will make use of the following formula: where ω is regarded as a function of υ ∈ [0, 1], ω ′ (υ) is its derivative, and υ decreases from 1 to 0 as τ increases from 0 to τ m .By setting υ = dω/dτ and d 2 ω/dτ 2 = dυ/dτ = υdυ/dω, we reduce Equation ( 25) to the following first-order differential equation: The first-order approximate solution of Equation ( 29) subject to the boundary condition ω υ=1 = 0 is given by the following: The substitution of (30) into Equation ( 28) yields the following: where we have introduced the following notation: Moreover, in view of (30), we obtain the following: It can be easily verified numerically that the first-order approximations (31) 1 and (33) 1 completely agree with the corresponding results obtained by Hunter.At the same, the asymptotic formula (31) 1 disagrees with Forney's result τ m − τ 0 m = O(ε 1/2 ) as ε → 0. From a practical point of view, it is of interest to evaluate the maximum contact force, F M , and the corresponding time moment, t M , such that F M = F(t M ).In the case of a Maxwell solid, according to Equations ( 12), ( 18), ( 19) and ( 24), for τ ∈ (0, τ m ), we have the following: where ω 0 and ω 1 are given by (30) 2 and (30) 3 , respectively.Without going into details of the first-order asymptotic analysis, the following can be shown: where c τ m and c ω m are given by ( 32) and (33) 2 , respectively.

Asymptotic Solution for the Unloading Stage
Now, we transform Equation ( 26) as follows: The first-order approximate solution υ ≃ υ 0 + ευ 1 to Equation ( 36) is given by the following: The end of the unloading stage is determined by the following value: where c ω m is defined by formula (33) 1 .The duration of the unloading stage can be evaluated as follows: and, in view of (37), we find the following: where we have introduced the following notation: Hence, the dimensional impact duration t 0 c = (w 0 m /V 0 )τ 0 c can be represented as follows: where τ 0 c = 2τ 0 m , and c τ c is given by ( 40).Again, by numerical check for the involved integrals, it can be easily established that formula (41) completely agrees with the corresponding Hunter's result.
Finally, according to Equation ( 22), the rebound velocity is found to be as follows: , and in view of ( 38) and (39), we arrive at the following Hunter's result for the coefficient of restitution: Thus, our calculations for w m , t c , and e have conformed to those by Hunter, obtained by a different method.

Impact for a Maxwell Solid
It is interesting that though numerical approaches for the viscoelastic Hertzian impact were developed in a number of papers [62][63][64], quite general numerical results were obtained not so long ago by Herrenbrück et al. [65] based on finite-element simulations for both the Maxwell model and the standard linear sold model.The numerical master curves were calculated for the maximum penetration w m scaled by the Hertzian elastic solution w 0 m , for the coefficient of restitution e, and for the maximum acceleration, which coincides with the relative maximum contact force F M /m, also scaled by the Hertzian solution F 0 M /m, where the following holds: In the case of a Maxwell solid, the creep function (18) 2 is now represented as follows: where τ R is the characteristic relaxation time, which is introduced instead of the inverse characteristic time η used before.
According to relations (33) 1 , (35) 1 , and (42), the above-mentioned impact parameters can be approximated as follows: where F 0 M is given by ( 43), and we have introduced the following notation: As it is seen from Figures 2 and 3 (see the inserts), the analytical approximations ( 46) can be used with less than 5% relative error for ε ≤ 0.2.It should be noted that different symbols in Figures 2 and 3  We note that, for a Maxwell solid, in view of ( 9), (10), and ( 18) 1 , we have μ(t) = −ηµ 0 exp(−ηt).It can be verified that Aksel [53] obtained the approximate formula e = 1 − 2.69ηw 0 m /V 0 , which can be transformed to the form (42) with the coefficient 0.914 instead of 4/9 ≈ 0.444.Thus, Aksel's asymptotic solution apparently contains a computational error, as it does not agree with the numerical solution shown in Figure 3b.Additionally, it should be noted that Aksel's approach suffers from a similar inconsistency as that of Forney's approach.It should be pointed out that Aksel's approach takes the contact radius as an independent variable, but in light of (30), the range of contact radius in the perturbed problem is larger than that in the Hertzian (unperturbed) problem.Remark 1.The accuracy of asymptotic approximations, which directly depends on the value of the small parameter ε (it is clear that the smaller ε, the smaller is the error), can be improved by the use of a Padé-like transformation [66] (see also [24] for an example).For instance, formula (46) 2 can be replaced with the following asymptotically equivalent one: Formula (48) yields the blue solid line in Figure 3b.As it is seen from the insert in Figure 3b (errors produced by the asymptotic approximations (46) 2 and (48) are negative and positive, respectively), the 5% error interval of formula (48) is 50% bigger than the interval of the same accuracy for formula (46) 2 .It is clear that the FEM numerical master curves [65] for the Maxwell viscoelastic model can be fitted with more accurate analytical approximations.For the main purpose of our study, it is important to highlight that the straight line in Figure 3b (which is produced by the first-order asymptotics (46) 2 ) coincides with the initial tangent line of the FEM master curve.This fact precisely supports the correctness of Hunter's asymptotic solution.

Impact for a Standard Linear Solid
By adopting the notation used in [65], the shear creep function in the case of a standard linear solid will be the following: where τ S is a characteristic time, α = (µ 0 − µ ∞ )/µ 0 , and µ ∞ is the relaxed shear modulus.The short-time approximation (as t tends toward zero) follows from (49) in the following form: By comparing formulas (44) 1 and (50), we conclude that the first-order approximations (46) still can be employed provided that the following holds: where τ S and α are the model parameters of a standard solid subjected to a spherical impact.As it may be seen from Figure 4, the short-time approximations can be utilized in a limited range of the dimensionless parameter ε.We note also that the Maxwell model ( 44) can be recovered from the standard linear solid model (49) in the limit as α tends toward zero, provided that τ S /α is fixed to be τ R (see also [46]).

Material Parameter Identification via Impact Testing
Let us recall that the analytical approximations for the characteristics of an impact derived in the framework of the Maxwell viscoelastic model exploit the following expansions about t = 0: Moreover, the characteristic relaxation time 1/η is assumed to be much larger than the impact duration, that is, ηt c ≪ 1 or, which is asymptotically the same, ηt 0 c ≪ 1.The analytical approximations are given in terms of the small dimensionless parameter ε, which, in view of (1) and (24), is proportional to ηt 0 c , where t 0 c is the Hertzian impact duration (see Equations ( 1) and (7) 2 ).Namely, Hertz's theory of impact yields the following characteristic time: where we have introduced the notation for the following compliance coefficient: Now, we consider the approximate formulas (35), which can be recast as follows: where we have introduced the following short-hand notation: We recall that F M and t M denote the maximum contact force and the corresponding time moment.We assume that at least one of these parameters of impact can be measured experimentally.The problem is to evaluate the material parameters C 0 and η from the impact data collected from several tests characterized by the following governing parameter: In view of (24) 4 , (53), and (58), formulas ( 55) and ( 56) can be represented as follows: First, we consider Equation ( 59) and note that this formula represents a linear relation between the relative maximum contact force F m /(mV 0 ) and the variable impact parameter V. Provided that both of them are measured in an experiment, the material parameters C −1 0 and η can be evaluated via linear regression by means of fitting the linear Formula (59) to the scaled experimental data.After that, the instantaneous shear elastic modulus will be given by the following: where the material Poisson's ratio ν, as usual, is supposed to be known.Second, in order to exhibit the method of linear regression, we transform Equation ( 60) to the following form: By fitting Equation ( 62) to the experimental data (RV 0 /m 2 ) 1/5 t M versus (m 2 /RV 0 ) 1/5 , we can evaluate C 0 and ηC 2 0 , from where, in view of (61), we readily obtain µ 0 .Meanwhile, the inverse characteristic relaxation time η is simply determined from the ratio of the two linear regression coefficients.
In the same way, formula (12) for the impact duration can be rewritten as follows: and eventually transformed to the following form: By comparing Equations ( 62) and ( 64), we readily see that a similar linear regression method can be designed for evaluating the parameters µ 0 and η from the impact duration data.It is clear that the same type of equation can be easily obtained for the duration t m of the loading stage.
It remains an open issue whether the impact test results for polymer rubber-like materials (for which the standard solid model is quite suitable) align with the asymptotic model in the appropriate range of the model variables.To apply the above-suggested linear-regression-method-based analysis to experimental data, accurate measurements of the maximum contact reaction F M or the contact duration parameters t M , t m , and t c are needed.Remark 2. In light of (13), ( 21), ( 24), ( 53), (54), and (58), Formula (46) 1 can be recast in the following form: Let us recall [67] that the ratio De = τ R /t c is called the Deborah number.In view of (45), we have τ R /t c ∼ 1/ε as ε → 0. Thus, formula (65), as well as other asymptotic formulas derived above in the framework of the Maxwell viscoelastic model, holds under the assumption De ≫ 1, that is, for the so-called fast-time impact.Further, it is clear that the first term on the right-hand side of Equation (65) represents the Hertzian solution, which is known to be self-similar [28].
Recently, using the first-and second-order asymptotic solutions obtained for the spherical impact of a Maxwell-type viscoelastic Winkler foundation, Maruoka [67] has proposed a crossover of scaling law, which may be regarded as the interference from the higher-class similarity parameters.A somewhat similar interpretation, which refers to the concept of intermediate asymptotics [68], can be extended to the case of viscoelastic Hertzian impact for a Maxwell solid.

Discussion
The Hertzian impact assumes a paraboloidal approximation φ(r) = r 2 /2R for the initial gap between two colliding elastic solids, which eventually leads to Hertz's contact law F = kw 3/2 .Shtaerman [69] and Galin [70] obtained the force-displacement relations F = k n w (2n+1)/(2n) and F = k λ w (λ+1)/λ for the monomial gap functions φ(r) = A n r 2n (n is an integer) and φ(r) = A λ r λ (λ ≥ 1 is a real number).The special case of a conical gap, φ(r) = A 1 r and F = k 1 w 2 , was earlier considered by Love [71].The corresponding generalizations of the elastic Hertzian impact model was given by Kilchevsky [72] and Graham [73].However, to the best of the author's knowledge, Hunter's model of viscoelastic impact was not extended to the case of the Shtaerman-Galin contact law Another open issue is to account for the target thickness, which is realized by a nonlinear force-displacement relation that looses the self-similarity scaling.In the case of quasi-static viscoelastic Hertzian contact, the thickness effect was considered by Argatov et al. [74].It would be of undoubted interest to incorporate this effect even into the elastic Hertzian impact model (see, e.g., [55,75]).
Still, a puzzle remains to be solved, and this concerns the inconsistency of Hunter's prediction for the duration of impact (t c < t 0 c ) with the experimental observations [55] that the effect of viscoelasticity increases the impact duration as compared with the elastic case (that is, t c > t 0 c ).For the linear viscoelastic Maxwell impact model (see, e.g., [46,76]), though the viscoelasticity effect increases the duration of impact, we have , where ζ is the loss factor.In other words, the viscoelastic dissipation effect on the duration of the loading stage is much stronger than the effect on the overall contact duration.Apparently, the answer to the raised question might be sought in the fact that Hunter's approximation t 1 ≃ 2t m − t is not asymptotically exact (strictly speaking, the sign ≃ should be replaced with ≈).At the same time, the approximate model (19) for the loading stage and the corresponding solutions (see Section 2.3) are asymptotically exact.It should be also noted that the FEM simulations performed by Diani et al. [57] for a generalized Maxwell model result in an impact duration smaller than that predicted by Hertz's theory.Hence, an accurate numerical study (with a precisely controllable accuracy) of the Hunter model for a Maxwell solid is needed to shed light on the impact duration issue.
A remark should be made about numerical solutions to the unilateral viscoelastic impact problem, where the force-displacement relation in the unloading stage is given by the two nonlinear integral equations ( 16) and ( 17) with the function t 1 (t) being determined by the equation a(t) = a(t 1 ), where a(t) is the current contact radius in the unloading stage (t > t m ), and a(t 1 ) is the same value of contact radius in the loading stage (t 1 < t m ).A number of numerical schemes have been designed in the literature [62][63][64], but still no in-depth reliable numerical study of the impact problem (e.g., for a Maxwell solid) has been published.It is noteworthy that the FEM simulations, while being very useful for the overall analysis [65], do not suit well for verifying asymptotic solutions, especially if the second-order smallness effects should be spotlighted.

Conclusions
The author's contribution to the theory of viscoelastic Hertzian impact concerns a new approach to the first-order asymptotic solution of the set of governing differential Equations ( 19) and (20) in the case of a Maxwell solid.Simple analytical approximations have been derived for the maximum contact force and the time to achieve it.It has been shown that the critique of Forney (1974) of the analytical solution obtained by Hunter (1957) is in error.By comparison with the FEM results obtained by Herrenbrück et al. (2015), Hunter's asymptotic solutions for the maximum penetration and contact reaction as well as for the coefficient of restitution have been supported.A linear-regression-based method is outlined for recovering two material parameters from the impact data collected from several tests characterized by the governing impact parameter RV 0 /m 2 .
To conclude, in view of its practical significance, the problem of viscoelastic impact requires further investigation.Acknowledgments: Support by the German Research Foundation and the Open Access Publication Fund of TU Berlin is gratefully acknowledged.

Conflicts of Interest:
The author declares no conflict of interest.

Nomenclature
The following nomenclature is used in this manuscript:

Figure 1 .
Figure 1.Initial impact configuration: (a) collision of two elastic spheres; (b) impact of a rigid sphere on an elastic half-space.

Figure 2 .
Figure 2. Master curves for the Maxwell model obtained by Herrenbrück et al. [65], using FEM simulations, and the analytical approximations (46) 1 and (46) 3 : relative maximum penetration (a) and relative maximum acceleration (b) as functions of the scaled characteristic relaxation time.

Figure 3 .
Figure 3. Master curves for the Maxwell model obtained by Herrenbrück et al. [65], using FEM simulations, and the analytical approximations (46) 2 (red curve) and (48) (blue curve): coefficient of restitution as a function of the scaled characteristic relaxation time (a) and the relative inverse relaxation time (b).

Figure 4 .
Figure 4. Rescaled master curves for the standard linear solid model obtained by Herrenbrück et al. [65] using FEM simulations and the analytical approximations (46): coefficient of restitution (a) and relative maximum acceleration (b) as functions of the scaled characteristic relaxation time.

Funding:
This research received no external funding.Data Availability Statement: Data is contained within the article.