Viscoelastic Hertzian Impact — Feature Article

1. Introduction

In his ground-laying paper [1] on the frictionless local contact of elastic solids, Heinrich Hertz estimated the half-duration of the normal dissipationless impact as

$$t_{\mathrm{m}}^0={\tau }_{\mathrm{m}}^0\frac{w_{\mathrm{m}}^0}{V_0},{\tau }_{\mathrm{m}}^0=\underset{0}{\overset{1}{\int }}\frac{\mathrm{d}\xi }{\sqrt{1-{\xi }^{5/2}}},$$

(1)

where $V_0$ is the initial relative velocity of approach, and $w_{\mathrm{m}}^0$ is the maximum value of the contact approach (evaluated from the initial contact moment).

Hertz’s theory of elastic impact can be represented as the initial-value problem for the second-order nonlinear differential equation

$$m\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{t}^2}=-k{w}^{3/2},w{|}_{t=0}=0,\frac{\mathrm{d}w}{\mathrm{d}t}{|}_{t=0}=V_0,$$

(2)

where w is the contact approach measured from the time, t, of initial contact, m is the equivalent mass, k is the stiffness coefficient in Hertz’s contact law $F=k{w}^{3/2}$, and F is the contact force (reaction).

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^{*}\sqrt{R}$ is determined in terms of the equivalent radius, R, and the effective elastic modulus, $E^{*}$, defined as

$$\frac{1}{R}=\frac{1}{R_1}+\frac{1}{R_2},\frac{1}{E^{*}}=\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}.$$

(3)

From Eqs. (2), the equation of energy conservation follows in the form

$$\frac{m}{2}\left(V_0^2-{\left(\frac{\mathrm{d}w}{\mathrm{d}t}\right)}^2\right)=\frac{2}{5}k{w}^{5/2},$$

(4)

which determines the maximum contact approach

$$w_{\mathrm{m}}^0={\left(\frac{5m}{k}\right)}^{2/5}{\left(\frac{V_0}{2}\right)}^{4/5},$$

(5)

achieved at the time moment $t=t_{\mathrm{m}}$ when $\mathrm{d}w/\mathrm{d}t=0$.

The same Hertzian equations (1)–(5) also apply in the case of impact of a rigid sphere of mass m onto an elastic half-space (see Figure 1b), as Eqs. (6) allow a passage to the limit as $E_1\to \infty$ and $R_2\to \infty$. 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 [2].

It should be noted that some of the Hertz model’s restrictions can be relaxed. For instance, in a routine manner Willis [3] extended Hertz’s theory of impact to anisotropic bodies, using his solution for the problem of local frictionless contact. By utilizing Bondareva’s solution [4] for a heavy elastic sphere on a rigid plane, Villaggio [5] 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 estimated by Hunter [6] (see also [7,8]) based on the analytical solution obtained by Miller and Pursey [9] 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 [10].

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_{\mathrm{c}}^0$, to be equal twice the time $t_{\mathrm{m}}^0$ to the maximum of the contact approach $w_{\mathrm{m}}^0$. It was shown by Hunter [6] and Deresiewicz [11] that the variation of the contact approach as a function of the time of contact can be well approximated by the simple formula

$$w\approx w_{\mathrm{m}}^{0}sin\frac{\pi t}{t_{\mathrm{c}}^0},t_{\mathrm{c}}^0=2t_{\mathrm{m}}^0.$$

(6)

The classical Hertz impact theory has been given substantial experimental verification [12] and, in particular, was extended to the frictional impact of anisotropic nonlinear elastic solids [13] and power-graded viscoelastic solids [14] as well as to tangential (oblique) impact [15,16,17] and impact with adhesion [18,19].

It is known [20,21] that the coefficient of restitution in collision of two perfectly elastic bodies equals 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 [6]) 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 [21]. A phenomenological approach (see, e.g., [22,23]) leads to the dissipative contact model

$$F=k{w}^{3/2}+\chi w^{\beta }\dot{w},$$

(7)

where $\dot{w}=\mathrm{d}w/\mathrm{d}t$ is the relative contact velocity, $\beta$ is a dimensionless constant, and $\chi$ denotes the hysteresis damping factor. As a rule, the coefficient $\chi$ is interpreted in terms of the coefficient of restitution [24]. By adopting a constitutive law similar to the Kelvin–Voigt model ${\sigma }_{zr}=2{\mu }_0{\epsilon }_{zr}+2\eta {\dot{\epsilon }}_{zr}$, where $\eta$ is the viscosity coefficient, Goldobin et al. [25] arrived at Eq. (7) with $\beta =1/2$ and $\chi$ being proportional to a linear combination of the shear and bulk viscosity coefficients (see also [26,27]).

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 [28], who complemented the analytical solution by Lee and Radok [29] 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 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_{\mathrm{c}}$. Later, Forney [30] questioned Hunter’s result about the impact duration. To date, this issue remains unresolved.

The case of a Kelvin–Voigt solid was considered by Khusid [31] (see also [14]) 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 el al. [32].

In what follows, we consider the normal impact of a rigid sphere on an isotropic viscoelastic half-space with a constant Poisson’s ratio, $\nu$, and a hereditary constitutive law

$${\sigma }_{zr}\left(t\right)=2\underset{0^{-}}{\overset{t}{\int }}\mu (t-t^{\prime })\frac{\mathrm{d}{ϵ}_{zr}\left(t^{\prime }\right)}{\mathrm{d}{t}^{\prime }}\mathrm{d}{t}^{\prime },$$

(8)

where $\mu \left(t\right)$ is the shear relaxation modulus, and $0^{-}$ indicates the instant immediately before the initial point of contact.

In his first-order perturbation analysis of the viscoelastic Hertzian impact, Aksel [33] applied the viscoelastic constitutive law in the form

$${\sigma }_{zr}\left(t\right)=2{\mu }_0{ϵ}_{zr}\left(t\right)+2\underset{0}{\overset{t}{\int }}\widehat{\mu }\left(t^{\prime }\right){ϵ}_{zr}(t-t^{\prime })\mathrm{d}{t}^{\prime },$$

(9)

where ${\mu }_0=\mu \left(0\right)$ is the instantaneous shear modulus, $\widehat{\mu }\left(t\right)=\mathrm{d}\mu \left(t\right)/\mathrm{d}t$ is the viscoelastic relaxation kernel, and it is tentatively assumed that ${ϵ}_{zr}\left(t\right)=0$ for $t<0$.

The problem of material parameters identification by means of impact tests was considered in a number of experimental studies [34,35,36,37]. Kren and Naumov [38] formulated the problem of determining the relaxation function $\mu \left(t\right)$ from the spherical impact loading history (impactor velocity, $V\left(t\right)$, contact force, $F\left(t\right)$, and contact approach, $w\left(t\right)$) without a priori adopting any material’s model. However, the problem with the Kren–Naumov method is that the Lee–Radok solution [29], which is valid only for the loading contact stage, was incorrectly applied for the unloading stage as well. It is still to note here that the approximation that utilizes for the restitution phase the same form of the equation of motion derived for compression is sometimes used for the sake of simplicity [39,40].

2. Hunter’s Model of the Viscoelastic Impact

2.1. Viscoelastic Hertzian Impact

For the sake of simplicity, we consider the single impactor configuration (Figure 1b) and start with Newton’s second law and the initial conditions

$$m\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{t}^2}=-F,w{|}_{t=0}=0,\frac{\mathrm{d}w}{\mathrm{d}t}{|}_{t=0}=V_0,$$

(10)

where F is the contact reaction.

Assuming a constant Poisson’s ratio $\nu$ and neglecting inertia effects in the viscoelastic half-space target, the Hertzian elastic solution can be generalized as follows [28,29]:

$$F\left(t\right)=\frac{8}{3(1-\nu )R}\underset{0}{\overset{t}{\int }}\mu (t-t^{\prime })\frac{\mathrm{d}}{\mathrm{d}{t}^{\prime }}\left(a^3\left(t^{\prime }\right)\right)\mathrm{d}{t}^{\prime },$$

(11)

$$w\left(t\right)=\frac{a^2\left(t\right)}{R},t\le t_{\mathrm{m}}.$$

(12)

Here, $\mu \left(t\right)$ is the shear relaxation modulus that describes the material’s response to an instantaneous unit shear deformation, $a\left(t\right)$ is the contact radius as a function of time t. It is clear that Eq. (12) is the Hertzian relation between the contact radius a and the contact approach w.

From Eqs. (10)–(12), it follows that

$$m\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{t}^2}=-\frac{8\sqrt{R}}{3(1-\nu )}\underset{0}{\overset{t}{\int }}\mu (t-t^{\prime })\frac{\mathrm{d}}{\mathrm{d}{t}^{\prime }}\left(w^{3/2}\left(t^{\prime }\right)\right)\mathrm{d}{t}^{\prime }.$$

(13)

It should be remembered that in the loading stage, the Lee–Radok solution (11), (12) as well as Eq. (13) are valid until the time moment, $t_{\mathrm{m}}$, of maximum penetration, $w_{\mathrm{m}}$, which, in view of Eq. (12), coincides with the moment of maximum contact radius, $a_{\mathrm{m}}$, and therefore, with the moment when the impactor velocity vanishes, i.e.,

$$a_{\mathrm{m}}=a\left(t_{\mathrm{m}}\right),w_{\mathrm{m}}=w\left(t_{\mathrm{m}}\right),\frac{\mathrm{d}w}{\mathrm{d}t}{|}_{t=t_{\mathrm{m}}}=0.$$

(14)

In the unloading stage ($t_{\mathrm{m}}<t<t_{\mathrm{c}}$), Hunter’s solution is given in terms of the function $t_1\left(t\right)$ that solves the equation $a\left(t\right)=a\left(t_1\right)$ for $t>t_{\mathrm{m}}$ and $t_1<t_{\mathrm{m}}$. Namely, the contact reaction is given by

$$F\left(t\right)=\frac{8}{3(1-\nu )R}\underset{0}{\overset{t_1\left(t\right)}{\int }}\mu (t-t^{\prime })\frac{\mathrm{d}}{\mathrm{d}{t}^{\prime }}\left(a^3\left(t^{\prime }\right)\right)\mathrm{d}{t}^{\prime },$$

(15)

whereas Eq. (12) is replaced with the following relation [28]:

$$Rw\left(t\right)=a^2\left(t\right)-\underset{t_{\mathrm{m}}}{\overset{t}{\int }}{\mu }^{-1}(t-t^{\prime })\frac{\mathrm{d}}{\mathrm{d}{t}^{\prime }}\left(\underset{\phantom{.}{t}_1\left(t^{\prime }\right)}{\overset{t^{\prime }}{\int }}\mu (t^{\prime }-t^{\prime \prime })\frac{\mathrm{d}}{\mathrm{d}{t}^{\prime \prime }}\left(a^2\left(t^{\prime \prime }\right)\right)\mathrm{d}{t}^{\prime \prime }\right)\mathrm{d}{t}^{\prime }.$$

(16)

Here, ${\mu }^{-1}\left(t\right)$ is the shear creep compliance (we keep the notation from [28]).

2.2. Impact for a Maxwell Solid

The reciprocal relaxation and creep functions for a Maxwell solid are

$$\mu \left(t\right)={\mu }_{0}exp(-\eta t),{\mu }^{-1}\left(t\right)=\frac{1}{{\mu }_0}(1+\eta t),$$

(17)

where ${\mu }_0$ is the instantaneous shear modulus, and $\eta$ is an inverse relaxation time.

The substitution of ${\left(17\right)}_1$ into Eq. (13), in view of the initial conditions (10), leads to the differential equation

$$\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{t}^2}+\eta \left(\frac{\mathrm{d}w}{\mathrm{d}t}-V_0\right)=-\frac{8\sqrt{R}{\mu }_0}{3(1-\nu )m}{w}^{3/2},t<t_{\mathrm{m}},$$

(18)

which is an exact result.

However, in the unloading (or withdrawal) stage, the application of Eqs. (17) allows to reduce Eq. (16) to a complicated nonlinear differential equation for $t_1\left(t\right)$. Nevertheless, by utilizing the simple approximation $t_1\approx 2t_{\mathrm{m}}-t$, Hunter derived the resulting governing differential equation

$$\frac{{\mathrm{d}}^2\mathcal{w}}{\mathrm{d}{t}^2}-\eta \left(3\frac{\mathrm{d}\mathcal{w}}{\mathrm{d}t}+V_0\right)=-\frac{8\sqrt{R}{\mu }_0}{3(1-\nu )m}{\mathcal{w}}^{3/2},t>t_{\mathrm{m}},$$

(19)

which is amendable to analytical treatment, but constitutes an approximate result.

Let us introduce the auxiliary function

$$\mathcal{w}=\frac{a^2\left(t\right)}{R},t\ge t_{\mathrm{m}},$$

(20)

which is related to the contact approach by the equation

$$\frac{\mathrm{d}\mathcal{w}}{\mathrm{d}t}=exp\left(-2\eta (t-t_{\mathrm{m}})\right)\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{a^2\left(t\right)}{R}\right).$$

(21)

We note that, in view of (12) and (20), the solution of Eq. (19) is subject to the boundary conditions (or the initial conditions for the withdrawal stage)

$$\mathcal{w}{|}_{t=t_{\mathrm{m}}}=w{|}_{t=t_{\mathrm{m}}},\frac{\mathrm{d}\mathcal{w}}{\mathrm{d}t}{|}_{t=t_{\mathrm{m}}}=0,$$

(22)

The contact duration $t_{\mathrm{c}}$ is determined by the condition $\mathcal{w}\left(t_{\mathrm{c}}\right)=0$, when the contact shrinks down.

In the dimensionless variables

$$\omega =\frac{a^2\left(t\right)}{R{w}_{\mathrm{m}}^0},\tau =\frac{V_{0}t}{w_{\mathrm{m}}^0},{\tau }_{\mathrm{m}}=\frac{V_0{t}_{\mathrm{m}}}{w_{\mathrm{m}}^0},\epsilon =\frac{\eta w_{\mathrm{m}}^0}{V_0},$$

(23)

where $w_{\mathrm{m}}^0$ is given by (5) with k being replaced by $8\sqrt{R}{\mu }_0/\left[3(1-\nu )\right]$, Eqs. (18) and (19) become

$$\frac{{\mathrm{d}}^2\omega }{\mathrm{d}{\tau }^2}+\epsilon \left(\frac{\mathrm{d}\omega }{\mathrm{d}\tau }-1\right)=-\frac{5}{4}{\omega }^{3/2},\tau <{\tau }_{\mathrm{m}},\omega {|}_{\tau =0}=0,\frac{\mathrm{d}\omega }{\mathrm{d}\tau }{|}_{\tau =0}=1,$$

(24)

$$\frac{{\mathrm{d}}^2\omega }{\mathrm{d}{\tau }^2}-\epsilon \left(3\frac{\mathrm{d}\omega }{\mathrm{d}\tau }+1\right)=-\frac{5}{4}{\omega }^{3/2},{\tau }_{\mathrm{m}}<\tau ,\omega {|}_{\tau ={\tau }_{\mathrm{m}}}={\omega }_{\mathrm{m}},\frac{\mathrm{d}\omega }{\mathrm{d}\tau }{|}_{\tau ={\tau }_{\mathrm{m}}}=0,$$

(25)

where ${\tau }_{\mathrm{m}}$ is the time-like point for which the solution to Eq. (24) first yields $\mathrm{d}\omega /\mathrm{d}\tau =0$, and the corresponding value of $\omega \left({\tau }_{\mathrm{m}}\right)$ will be denoted by ${\omega }_{\mathrm{m}}$.

2.3. Asymptotic Solution for the Loading Stage

The problem with Forney’s critique [30] of Hunter’s approximate solution is in the integral decomposition representation

$${\tau }_{\mathrm{m}}=\underset{0}{\overset{{\omega }_{\mathrm{m}}}{\int }}\frac{\mathrm{d}\omega }{\dot{\omega }}=\underset{0}{\overset{1}{\int }}\frac{\mathrm{d}\omega }{\dot{\omega }}+\underset{1}{\overset{{\omega }_{\mathrm{m}}}{\int }}\frac{\mathrm{d}\omega }{\dot{\omega }},$$

(26)

where $\dot{\omega }=\mathrm{d}\omega /\mathrm{d}\tau$, because the function $\omega$ cannot be regarded as a perturbation of the limit ($\epsilon =0$) elastic solution ${\omega }_0\left(\tau \right)$ that is defined only on the interval $[0,{\tau }_{\mathrm{m}}^0]$, where ${\tau }_{\mathrm{m}}^0$ is the dimensionless Hertzian half duration of impact. At the same time, the second integral on the right-hand side of Eq. (26) refers to the interval $\omega \in (1,{\omega }_{\mathrm{m}})$, which corresponds to $\tau \in ({\tau }_{\mathrm{m}}^0,{\tau }_{\mathrm{m}})$, that is outside the interval of validity of the limit solution.

That is why, we will make use of the formula

$${\tau }_{\mathrm{m}}=\underset{1}{\overset{0}{\int }}\frac{\mathrm{d}\omega }{\upsilon }=-\underset{0}{\overset{1}{\int }}{\omega }^{\prime }\left(\upsilon \right)\frac{\mathrm{d}\upsilon }{\upsilon },$$

(27)

where $\omega$ is regarded as a function of $\upsilon \in [0,1]$, ${\omega }^{\prime }\left(\upsilon \right)$ is its derivative, and $\upsilon$ decreases from 1 to 0 as $\tau$ increases from 0 to ${\tau }_{\mathrm{m}}$.

By setting $\upsilon =\mathrm{d}\omega /\mathrm{d}\tau$ and ${\mathrm{d}}^2\omega /\mathrm{d}{\tau }^2=\mathrm{d}\upsilon /\mathrm{d}\tau =\upsilon \mathrm{d}\upsilon /\mathrm{d}\omega$, we reduce Eq. (24) to the first-order differential equation

$$\frac{\mathrm{d}\omega }{\mathrm{d}\upsilon }=-\frac{\upsilon }{(5/4){\omega }^{3/2}+\epsilon (\upsilon -1)},\upsilon \in (0,1).$$

(28)

The first-order approximate solution to Eq. (28) subject to the boundary condition $\omega {|}_{\upsilon =1}=0$ is given by

$$\omega \simeq {\omega }_0+\epsilon {\omega }_1,{\omega }_0={(1-{\upsilon }^2)}^{2/5},{\omega }_1=\frac{{(4/5)}^2}{{(1-{\upsilon }^2)}^{3/5}}\underset{\upsilon }{\overset{1}{\int }}\frac{\xi (1-\xi )\mathrm{d}\xi }{{(1-{\xi }^2)}^{3/5}}.$$

(29)

The substitution of (29) into Eq. (27) yields

$${\tau }_{\mathrm{m}}\simeq {\tau }_{\mathrm{m}}^0+\epsilon c_{\mathrm{m}}^{\tau },{\tau }_{\mathrm{m}}^0=\frac{4}{5}\underset{0}{\overset{1}{\int }}\frac{\mathrm{d}\xi }{{(1-{\xi }^2)}^{3/5}},$$

(30)

$$c_{\mathrm{m}}^{\tau }=\frac{2^4}{5^2}\underset{0}{\overset{1}{\int }}\frac{(1-\xi )\mathrm{d}\xi }{{(1-{\xi }^2)}^{6/5}}-\frac{2^5·3}{5^3}\underset{0}{\overset{1}{\int }}\frac{1}{{(1-{\upsilon }^2)}^{8/5}}\underset{\upsilon }{\overset{1}{\int }}\frac{\xi (1-\xi )}{{(1-{\xi }^2)}^{3/5}}\mathrm{d}\xi \mathrm{d}\upsilon .$$

(31)

Moreover, in view of (29), we obtain

$${\omega }_{\mathrm{m}}\simeq 1+\epsilon c_{\mathrm{m}}^{\omega },c_{\mathrm{m}}^{\omega }=\frac{2^4}{5^2}\underset{0}{\overset{1}{\int }}\frac{\xi (1-\xi )}{{(1-{\xi }^2)}^{3/5}}\mathrm{d}\xi .$$

(32)

It can be easily verified numerically that the first-order approximations ${\left(30\right)}_1$ and ${\left(32\right)}_1$ completely agree with the corresponding results obtained by Hunter. At the same, the asymptotic formula ${\left(30\right)}_1$ disagrees with Forney’s result ${\tau }_{\mathrm{m}}-{\tau }_{\mathrm{m}}^0=O\left({\epsilon }^{1/2}\right)$ as $\epsilon \to 0$.

From a practical point of view, it is of interest to evaluate the maximum contact force, $F_{\mathrm{M}}$, and the corresponding time moment, $t_{\mathrm{M}}$, such that $F_{\mathrm{M}}=F\left(t_{\mathrm{M}}\right)$. In the case of a Maxwell solid, according to Eqs. (11), (17), (18), and (23), for $\tau \in (0,{\tau }_{\mathrm{m}})$ we have

$$\begin{array}{ccc}\hfill \frac{w_{\mathrm{m}}^0}{m{V}_0^2}F& =& \frac{5}{4}{\omega }^{3/2}+\epsilon \left(\frac{\mathrm{d}\omega }{\mathrm{d}\tau }-1\right)\hfill \\ \hfill & \simeq & \frac{5}{4}{\omega }_0^{3/2}+\epsilon \left(\frac{15}{8}{\omega }_0^{1/2}{\omega }_1-1+\upsilon \right),\hfill \end{array}$$

(33)

where ${\omega }_0$ and ${\omega }_1$ are given by ${\left(29\right)}_2$ and ${\left(29\right)}_3$, respectively.

It can be shown that

$$\frac{w_{\mathrm{m}}^0}{m{V}_0^2}{F}_{\mathrm{M}}\simeq \frac{5}{4}-\epsilon \left(1-\frac{15}{8}{c}_{\mathrm{m}}^{\omega }\right),\frac{V_0}{w_{\mathrm{m}}^0}{t}_{\mathrm{M}}\simeq {\tau }_{\mathrm{m}}^0-\epsilon \left(\frac{8}{15}-c_{\mathrm{m}}^{\tau }\right),$$

(34)

where $c_{\mathrm{m}}^{\tau }$ and $c_{\mathrm{m}}^{\omega }$ are given by (31) and ${\left(32\right)}_2$, respectively.

2.4. Asymptotic Solution for the Unloading Stage

Now, we transform Eq. (25) as

$$\upsilon \frac{\mathrm{d}\upsilon }{\mathrm{d}\omega }=-\frac{5}{4}{\omega }^{3/2}+\epsilon (3\upsilon +1),\omega \in (0,{\omega }_{\mathrm{m}}),\upsilon {|}_{\omega ={\omega }_{\mathrm{m}}}=0.$$

(35)

The first-order approximate solution $\upsilon \simeq {\upsilon }_0+\epsilon {\upsilon }_1$ to Eq. (35) is given by

$${\upsilon }_0=-\sqrt{{\omega }_{\mathrm{m}}^{5/2}-{\omega }^{5/2}},{\upsilon }_1=-\frac{1}{\sqrt{{\omega }_{\mathrm{m}}^{5/2}-{\omega }^{5/2}}}\underset{\omega }{\overset{{\omega }_{\mathrm{m}}}{\int }}\left[3\sqrt{{\omega }_{\mathrm{m}}^{5/2}-{\varphi }^{5/2}}-1\right]\mathrm{d}\varphi .$$

(36)

The end of the unloading stage is determined by the value

$${\upsilon }_{\mathrm{c}}=\upsilon {|}_{\omega =0}\simeq -1+\epsilon \left(1-\frac{5}{4}{c}_{\mathrm{m}}^{\omega }-3\underset{0}{\overset{1}{\int }}\sqrt{1-{\xi }^{5/2}}\mathrm{d}\xi \right),$$

(37)

where $c_{\mathrm{m}}^{\omega }$ is defined by formula ${\left(32\right)}_1$.

The duration of the unloading stage can be evaluated as

$${\tau }_{\mathrm{c}}-{\tau }_{\mathrm{m}}=\underset{1}{\overset{{\upsilon }_{\mathrm{c}}}{\int }}\frac{\mathrm{d}\omega }{\upsilon }=-\underset{0}{\overset{{\omega }_{\mathrm{m}}}{\int }}\frac{\mathrm{d}\omega }{\upsilon \left(\omega \right)},$$

and, in view of (36), we find that

$${\tau }_{\mathrm{c}}-{\tau }_{\mathrm{m}}\simeq {\tau }_{\mathrm{m}}^0-\epsilon c_{\mathrm{c}}^{\tau },$$

(38)

where we have introduced the notation

$$c_{\mathrm{c}}^{\tau }=\frac{1}{4}{\tau }_{\mathrm{m}}^0{c}_{\mathrm{m}}^{\omega }+\underset{0}{\overset{1}{\int }}\frac{1}{{(1-{\xi }^{5/2})}^{3/2}}\underset{\xi }{\overset{1}{\int }}\left[3\sqrt{1-{\zeta }^{5/2}}-1\right]\mathrm{d}\zeta \mathrm{d}\xi .$$

(39)

Hence, the dimensional impact duration $t_{\mathrm{c}}^0=(w_{\mathrm{m}}^0/V_0){\tau }_{\mathrm{c}}^0$ can be represented as

$$t_{\mathrm{c}}\simeq t_{\mathrm{c}}^0\left(1+\eta t_{\mathrm{c}}^0\frac{(c_{\mathrm{m}}^{\tau }-c_{\mathrm{c}}^{\tau })}{{\left({\tau }_{\mathrm{c}}^0\right)}^2}\right),$$

(40)

where ${\tau }_{\mathrm{c}}^0=2{\tau }_{\mathrm{m}}^0$, and $c_{\mathrm{c}}^{\tau }$ is given by (39).

Again, by numerical check for the involved integrals it can be easily established that formula (40) completely agrees with the corresponding Hunter’s result.

Finally, according to Eq. (21), the rebound velocity is found to be

$$V_{\mathrm{c}}=V_{0}exp\left(-2\eta (t_{\mathrm{c}}-t_{\mathrm{m}})\right)\frac{\mathrm{d}\omega }{\mathrm{d}\tau }{|}_{\tau ={\tau }_{\mathrm{c}}},$$

that is $V_{\mathrm{c}}/V_0\simeq \left(1-\epsilon ({\tau }_{\mathrm{c}}-{\tau }_{\mathrm{m}})\right){\upsilon }_{\mathrm{c}}$, and, in view of (37) and (38), we arrive at the following Hunter’s result for the coefficient of restitution:

$$e\simeq 1-\frac{4}{9}\eta t_{\mathrm{c}}^0.$$

(41)

Thus, our calculations for $w_{\mathrm{m}}$, $t_{\mathrm{c}}$, and e have conformed those by Hunter, obtained by different method.

3. Comparison with the FEM Solution

3.1. Impact for a Maxwell Solid

It is interesting that though numerical approaches for the viscoelastic Hertzian impact were developed in a number of papers [40,41,42], quite general numerical results were obtained not so long ago by Herrenbrück et al. [43] 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_{\mathrm{m}}$ scaled by the Hertzian elastic solution $w_{\mathrm{m}}^0$, for the coefficient of restitution e as well as for the maximum acceleration, which coincides with the relative maximum contact force $F_{\mathrm{M}}/m$, also scaled by the Hertzian solution $F_{\mathrm{M}}^0/m$, where

$$F_{\mathrm{M}}^0=\frac{5}{4}\frac{m{V}_0^2}{w_{\mathrm{m}}^0}.$$

(42)

In the case of a Maxwell solid, the creep function ${\left(\right)}_2$ is now represented as

$${\mu }^{-1}\left(t\right)=\frac{1}{{\mu }_0}\left(1+\frac{t}{{\tau }_{\mathrm{R}}}\right),{\tau }_{\mathrm{R}}=\frac{1}{\eta },$$

(43)

where ${\tau }_{\mathrm{R}}$ is the characteristic relaxation time, which is introduced instead of the inverse characteristic time $\eta$ used before.

In view of ${\left(1\right)}_1$, ${\left(23\right)}_4$, and ${\left(43\right)}_2$, we have

$$\frac{{\tau }_{\mathrm{R}}}{t_{\mathrm{m}}^0}=\frac{1}{\epsilon {\tau }_{\mathrm{m}}^0},\epsilon =\frac{t_{\mathrm{m}}^0}{{\tau }_{\mathrm{m}}^0{\tau }_{\mathrm{R}}},$$

(44)

where ${\tau }_{\mathrm{m}}^0$ is given by ${\left(1\right)}_2$.

According to relations ${\left(32\right)}_1$, ${\left(34\right)}_1$, and (41), the mentioned above impact parameters can be approximated as

$$\frac{w_{\mathrm{m}}}{w_{\mathrm{m}}^0}=1+\epsilon c_{\mathrm{m}}^{\omega },e=1-\epsilon \frac{8}{9}{\tau }_{\mathrm{m}}^0,\frac{F_{\mathrm{M}}}{F_{\mathrm{M}}^0}=1-\epsilon c_{\mathrm{M}}^F,$$

(45)

where $F_{\mathrm{M}}^0$ is given by (42), and we have introduced the notation

$$c_{\mathrm{M}}^F=\frac{4}{5}-\frac{3}{2}{c}_{\mathrm{m}}^{\omega }.$$

(46)

As it is seen from Figure 2 and Figure 3 (see the inserts), the analytical approximations (45) can be used with less than 5% relative error for $\epsilon \le 0.2$.

We note that for a Maxwell solid, in view of (8), (9), and ${\left(17\right)}_1$, we have $\widehat{\mu }\left(t\right)=-\eta {\mu }_{0}exp(-\eta t)$. It can be verified that Aksel [33] obtained the approximate formula $e=1-2.69\eta w_{\mathrm{m}}^0/V_0$, which can be transformed to the form (41) with the coefficient $0.914$ instead of $4/9\approx 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.

3.2. Impact for a Standard Linear Solid

By adopting the notation used in [43], the shear creep function in the case of a standard linear solid will be

$${\mu }^{-1}\left(t\right)=\frac{1}{(1-\alpha ){\mu }_0}\left\{1-\alpha exp\left(-(1-\alpha )\frac{t}{\tau }\right)\right\},$$

(47)

where $\alpha =({\mu }_0-{\mu }_{\infty })/{\mu }_0$, and ${\mu }_{\infty }$ is the relaxed shear modulus.

The short-time approximation (as t tends to zero) follows from (47) in the form

$${\mu }^{-1}\left(t\right)\simeq \frac{1}{{\mu }_0}\left(1+\frac{t}{\tau /\alpha }\right).$$

(48)

By comparing formulas ${\left(43\right)}_1$ and (48), we conclude that the first-order approximations (45) still can be employed provided

$$\epsilon =\frac{\alpha t_{\mathrm{m}}^0}{{\tau }_{\mathrm{m}}^0\tau },$$

(49)

where $\tau$ and $\alpha$ are the model parameters of a standard solid subjected to 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 $\epsilon$. We note also that the Maxwell model (43) can be recovered from the standard linear solid model in the limit as $\alpha$ tends to zero (see also [44]).

4. Material Parameters Identification via Impact Testing

Let us recall that the analytical approximations for the characteristics of impact derived in the framework of the Maxwell viscoelastic model exploit the following expansions about $t=0$:

$$\mu \left(t\right)={\mu }_0(1-\eta t+\dots ),{\mu }^{-1}\left(t\right)=\frac{1}{{\mu }_0}(1+\eta t-\dots ).$$

(50)

Moreover, the characteristic relaxation time $1/\eta$ is assumed to be much larger than the impact duration, that is $\eta t_{\mathrm{c}}\ll 1$ or, which is asymptotically the same, $\eta t_{\mathrm{c}}^0\ll 1$.

The analytical approximations are given in terms of the small dimensionless parameter $\epsilon$, which, in view of (1) and (23), is proportional to $\eta t_{\mathrm{c}}^0$, where $t_{\mathrm{c}}^0$ is the Hertzian impact duration (see Eqs. (1) and ${\left(6\right)}_2$). Namely, Hertz’s theory of impact yields the characteristic time

$$\frac{w_{\mathrm{m}}^0}{V_0}={\mathcal{C}}_0{\left(\frac{m^2}{R{V}_0}\right)}^{1/5},$$

(51)

where we have introduced the notation for the compliance coefficient

$${\mathcal{C}}_0={\left(\frac{15(1-\nu )}{32{\mu }_0}\right)}^{2/5}.$$

(52)

Now, we consider the approximate formulas (34), which can be recast as

$$\frac{F_{\mathrm{m}}}{m{V}_0}=\frac{V_0}{w_{\mathrm{m}}^0}\left(\frac{5}{4}-\epsilon c_{\mathrm{M}}^F\right),$$

(53)

$$t_{\mathrm{M}}=\frac{w_{\mathrm{m}}^0}{V_0}\left({\tau }_{\mathrm{m}}^0-\epsilon c_{\mathrm{M}}^t\right),$$

(54)

where we have introduced the short-hand notation

$$c_{\mathrm{M}}^F=1-\frac{15}{8}{c}_{\mathrm{m}}^{\omega },c_{\mathrm{M}}^t=\frac{8}{15}-c_{\mathrm{m}}^{\tau }.$$

(55)

We recall that $F_{\mathrm{M}}$ and $t_{\mathrm{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 ${\mathcal{C}}_0$ and $\eta$ from the impact data collected from several tests characterized by the governing parameter

$$\mathcal{V}={\left(\frac{R{V}_0}{m^2}\right)}^{1/5}.$$

(56)

In view of ${\left(23\right)}_4$, (51), and (56), formulas (53) and (54) can be represented as follows:

$$\frac{F_{\mathrm{M}}}{m{V}_0}=\frac{5}{4}{\mathcal{C}}_0^{-1}\mathcal{V}-c_{\mathrm{M}}^F\eta ,$$

(57)

$$t_{\mathrm{M}}={\tau }_{\mathrm{m}}^0{\mathcal{C}}_0{\mathcal{V}}^{-1}-c_{\mathrm{M}}^t\eta {\mathcal{C}}_0^2{\mathcal{V}}^{-2}.$$

(58)

First, we consider Eq. (57) and note that this formula represents a linear relation between the relative maximum contact force $F_{\mathrm{m}}/\left(m{V}_0\right)$ and the variable impact parameter $\mathcal{V}$. Provided that both of them are measured in experiment, the material parameters ${\mathcal{C}}_0^{-1}$ and $\eta$ can be evaluated via linear regression by means of fitting the linear formula (57) to the scaled experimental data. After that the instantaneous shear elastic modulus will be given by

$${\mu }_0=\frac{15(1-\nu )}{32}{\mathcal{C}}_0^{-5/2},$$

(59)

where the material Poisson’s ratio $\nu$, as usual, is supposed to be known.

Second, in order to exhibit the method of linear regression, we transform Eq. (58) to the form

$$\mathcal{V}{t}_{\mathrm{M}}={\tau }_{\mathrm{m}}^0{\mathcal{C}}_0-c_{\mathrm{M}}^{\tau }\eta {\mathcal{C}}_0^2{\mathcal{V}}^{-1}.$$

(60)

By fitting Eq. (60) to the experimental data ${(R{V}_0/m^2)}^{1/5}{t}_{\mathrm{M}}$ versus ${(m^2/R{V}_0)}^{1/5}$, we can evaluate ${\mathcal{C}}_0$ and $\eta {\mathcal{C}}_0^2$, from where, in view of (59), we readily get ${\mu }_0$. Meanwhile, the inverse characteristic relaxation time $\eta$ is simply determined from the ratio of the linear regression coefficients.

In the same way, formula (11) for the impact duration can be rewritten as

$$t_{\mathrm{c}}=\frac{w_{\mathrm{m}}^0}{V_0}\left(2{\tau }_{\mathrm{m}}^0+\epsilon c_{\mathrm{c}}^t\right)$$

(61)

and eventually transformed to the form

$$\mathcal{V}{t}_{\mathrm{c}}=2{\tau }_{\mathrm{m}}^0{\mathcal{C}}_0+c_{\mathrm{c}}^t\eta {\mathcal{C}}_0^2{\mathcal{V}}^{-1}.$$

(62)

By comparing Eqs. (60) and (62), we readily see that a similar linear regression method can be designed for evaluating the parameters ${\mu }_0$ and $\eta$ from the impact duration data.

5. Discussion

The Hertzian impact assumes a paraboloidal approximation $\phi \left(r\right)=r^2/2R$ for the initial gap between two colliding elastic solids, which eventually leads to Hertz’s contact law $F=k{w}^{3/2}$. Shtaerman [45] and Galin [46] obtained the force-displacement relations $F=k_n{w}^{(2n+1)/\left(2n\right)}$ and $F=k_{\lambda }{w}^{(\lambda +1)/\lambda }$ for the gap functions $\phi \left(r\right)=A_n{r}^{2n}$ (n is an integer) and $\phi \left(r\right)=A_{\lambda }{r}^{\lambda }$ ($\lambda \ge 1$ is a real number). The special case of a conical gap, $\phi \left(r\right)=A_{1}r$ and $F=k_1{w}^2$, was earlier considered by Love [47]. The corresponding generalizations of the elastic Hertzian impact model was given by Kilchevsky [48]. 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 $F=k_{\lambda }{w}^{(\lambda +1)/\lambda }$.

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. [49]. It would be of undoubted interest to incorporate this effect even into the elastic Hertzian impact model (see, e.g., [35,50]).

Still a puzzle remains to be solved, and this concerns the inconsistency of Hunter’s prediction for the duration of impact ($t_{\mathrm{c}}<t_{\mathrm{c}}^0$) with the experimental observations [35] that the effect of viscoelasticity increases the impact duration as compared to the elastic case (that is, $t_{\mathrm{c}}>t_{\mathrm{c}}^0$). For the linear viscoelastic Maxwell impact model (see, e.g., [44,51]), though the viscoelasticity effect increases the duration of impact, we have $t_{\mathrm{m}}=t_{\mathrm{m}}^0(1+O\left(\zeta \right))$ and $t_{\mathrm{c}}=t_{\mathrm{c}}^0(1+O\left({\zeta }^2\right))$, where $\zeta$ is the loss factor, that is the effect on the duration of the loading stage is much stronger than the effect of the overall contact duration. Apparently, the answer to the raised question might be sought in the fact that Hunter’s approximation $t_1\simeq 2t_{\mathrm{m}}-t$ is not asymptotically exact (strictly speaking, the sign ≃ should be replaced with ≈). At the same time, the approximate model (18) 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. [37] for a generalized Maxwell model result in the impact duration smaller than that predicted by Hertz’s theory.

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 (15) and (16) with the function $t_1\left(t\right)$ being determined by the equation $a\left(t\right)=a\left(t_1\right)$, where $a\left(t\right)$ is the current contact radius in the unloading stage ($t>t_{\mathrm{m}}$), and $a\left(t_1\right)$ is the same value of contact radius in the loading stage ($t_1<t_{\mathrm{m}}$). A number of numerical schemes have been designed in the literature [40,41,42], but still no in-depth numerical study of the impact problem (e.g., for a Maxwell solid) has been published. It is to note that the FEM simulations, while being very useful for the overall analysis [43], do not suit well for verifying asymptotic solutions, if the second order smallness effects should be spotlighted.

To conclude, in view of its practical significance, the problem of viscoelastic impact requires a further investigation.