The Boussinesq flat-punch indentation problem within the context of linearized viscoelasticity

The Boussinesq problem, namely the indentation of a flat-ended cylindrical punch into a viscoelastic half-space is studied. We assume a linear viscoelastic model wherein the linearized strain is expressed as a function of the stress. However, this expression is not invertible, which makes the problem very interesting. Based on the Papkovich–Neuber representation in potential theory and using the Fourier–Bessel transform for axisymmetric bodies, an analytical solution of the resulting time-dependent integral equation is constructed. Consequently, distribution of the displacement and the stress fields in the half space with respect to time is obtained in the closed form. © 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Indentation is probably one of the most often used experimental procedures to determine the material properties of a body. Unlike other experimental procedures that are global and provide an average value of the material moduli, indentation can be used to determine the material moduli locally. Thus, in inhomogeneous materials, indentation can be used to determine how the material properties vary over the body. Nano-indentation can also be used to determine the properties of very small bodies, and studying the problem of indentation is essential when one is concerned with determining the material properties of inhomogeneous bodies.
In this study, we are interested in analyzing the problem of the response of a linear viscoelastic body when it is subject to indentation. With regard to linear viscoelasticity one can express the stress in terms of the history of strain or vice-versa. However, the constitutive relations that are used are invariably invertible. The constitutive relation that we consider expresses the linearized strain in terms of the history of the stress, an interesting feature of the constitutive relation being that it cannot be inverted to express the stress as a function of the strain. We have purposely chosen such a model. Providing an expression of the strain in terms of the history of the stress is in keeping with causality as force/stress is the cause and the deformation/strain is the effect. Moreover, such an approach can be viewed as a simplification of more general implicit constitutive relations wherein the one has an implicit relationship between the stress and the strain (see Rajagopal, 2003 ). Boltzmann (1874) was the first to introduce an integral model to describe linear viscoelastic response. Based on the linear viscoelastic model Lockett (1972) introduced nonlinear models wherein stress is expressed as a nonlinear function of strain or the strain as a nonlinear function of the stress. Earlier, Green and Rivlin (1957) and Green, Rivlin, and Spencer (1959) had introduced nonlinear models to describe the response of viscoelasticity, but these models are not amenable to analysis. In view of this, Fung (1981) introduced a quasilinear viscoelastic model, however in his model the stress was given as a function of the history of the strain. Muliana, Rajagopal, and Wineman (2013) developed quasilinear viscoelastic models wherein the strain is function of the history of the stress. The model that we study can be viewed as a linearization of a model belonging to the class of viscoelastic models introduced by Muliana et al. (2013) wherein the linearized strain is expressed as a nonlinear function of the history of the stress.
It is well known that solution of problems in linear viscoelasticity can be obtained by appealing to a correspondence principle between the solutions in linearized elasticity and linear viscoelasticity. However, for the correspondence principle to be applicable certain conditions with regard to the deformation and the loading have to be satisfied. Wineman and Rajagopal (20 0 0) provide several examples where the correspondence principle fails. The problem under consideration is an example where one cannot appeal to the correspondence principle. Thus, we have to study the governing equations for the viscoelastic solid that we are indenting without appealing to the correspondence principle. Finally it is necessary to remark that while the solution that we seek are time dependent, we are not considering inertial effects. Thus, we are only interested in solving the balance of linear momentum in the absence of inertial effects. The time dependence of the solutions is a consequence of the material under consideration being viscoelastic.
In this paper we consider the response for the linearized viscoelastic model which takes the following integral form for time t > 0 where F is the compliance tensor. The positive kernel K that appears in the Volterra convolution operator in (1) is usually assumed to be given by the exponential sum with parameters K(0) , K 1 , . . . , K N , τ 1 , . . . τ N ≥ 0 , that characterize the generalized creep. If K ( t ) ≡ 1, then (1) reduces to the classical Hooke's law In general, the representation (1) cannot be inverted even being linear in T . If N = 1 and K(0) = 0 in (2) , then K (t) = −K (t) /τ 1 and K (0) = K 1 /τ 1 , on differentiating (1) with respect to t we derive the linearized Kelvin-Voigt model where dot stands for the time derivative. For relevant studies of the Kelvin-Voigt model (see Bulíček, Málek, and Rajagopal (2012) , Erbay and Ş engül (2015) and Itou, Kovtunenko, and Rajagopal (2018) ). The constitutive model (1) is a specific case of the more general class of nonlinear viscoelastic models when we substitute FT with a nonlinear relationship F (T ) . In the context of quasi-linear viscoelasticity, such models were investigated with respect to well-posedness in Itou, Kovtunenko, and Rajagopal (2019a) . For related issues of nonlinear relationships F (T ) within the context of limiting small strain ε ≤ M , we refer the reader to Itou, Kovtunenko, and Rajagopal (2017) , Itou, Kovtunenko, and Rajagopal (2019b) and Kulvait, Málek, and Rajagopal (2019) . Within the context of contact mechanics, an important subclass of the problems adjacent to cracks subject to non-penetration was developed by Khludnev and Kovtunenko (20 0 0) , and Khludnev and Sokołowski (1997) .  studied the problem of cracks in a linear viscoelastic body, andHintermüller, Kovtunenko, andKunisch (2009) , Itou, Kovtunenko, and , Khludnev and Shcherbakov (2018) , Lazarev and Rudoy (2014) studied the problem of a non-penetrating crack within the context of the linearized elastic relationship.
In the current paper, we study the linearized viscoelastic model (1) with respect to time dependent contact indentation by a rigid punch. The indentation test on viscoelastic materials has technological importance because of the ubiquitous use of inhomogeneous materials such as composites and inhomogeneous polymeric bodies. For instance, the identification of viscoelastic properties of bitumen is determined by indentation tests (see in Jäger, Lackner, & Eberhardsteiner, 2007 ). For the theory of the surface loading of a half-space we refer the reader to Argatov and Mishuris (2018) , Galin (2008) , Popov and Heß (2015) , Selvadurai (2001) and others. This is closely related to the Boussinesq problem (see Boussinesq, 1885 ) h r 0 consisting in finding the deformation of a half-space subject to a concentrated load applied at the origin on its plane boundary and perpendicular to it. We study the deformation of a viscoelastic half-space under normal indentation of a flat-ended cylindrical punch, then pass to the limit as the punch radius tends to zero thus solving the Boussinesq problem.
The indentation problem is formulated as follows. For spatial points i, j=1 which solve the equilibrium equations neglecting body forces and inertial effects (where the index j after comma means the partial derivative with respect to x j ) and the constitutive relation for the viscoelastic body (1) takes the form At infinity the following asymptotic condition is assumed, At the plane boundary { x 3 = 0 } the cylindrical punch with the flat end of radius h > 0 is indented perpendicularly with a given normal force F ( t ) ≥ 0 such that Neglecting sliding and the effects of friction, the contact conditions at x 3 = 0 consist in the tangential stress-free conditions with the unknown indentation depth D ( t ) ≥ 0 to be found. Due to the axisymmetric structure of the problem, relations (5) -(12) can be rewritten in the cylindrical coordinates , and these relations depend on t, r and z , but not on θ . The geometric configuration determined by u z ( t, r , 0) and T zz ( t, r , 0) at the boundary z = 0 for fixed t is portrayed in Fig. 1 . It is worth recalling that we cannot invert the implicit relationship (6) and substitute the stress T into the balance Eq. (5) as it is usually done when solving indentation problems expressed by explicit linear relations (i.e., when T is expressed in the terms of ε [ u ]). Therefore, in order to construct a solution to the indentation problem (5) -(12) in the closed form, we rely on the following arguments.
We introduce an auxiliary vector v (t, x ) = (v 1 , v 2 , v 3 ) and rewrite the response (6) equivalently as where the tensor ε Indeed, if the Eq. (14) holds, then applying ε [ ·] to both sides we can interchange ε and the integral Volterra operator I t and with the help of (13) to derive (6) . The converse holds true up to a rigid motion, which vanishes due to the condition at infinity (8) . We also note that the variable v is redundant and it can be eliminated after solving the problem.
After setting ψ r = ψ θ ≡ 0 in (22) and from (20) the formula for the stresses is derived (compare with Korsunsky, 1995;Zhou & Gao, 2013 ) T rr = 2 μ(−zψ z,rr + 2 νψ z,z − φ ,rr ) , T θθ = 2 μ − z r ψ z,r + 2 νψ z,z − 1 r φ ,r , We apply to Eq. (24) the Fourier-Bessel (also called Hanckel) transform using the Bessel function J m of order m = 0 , 1 , to obtain (see Harding & Sneddon, 1945 ) which obeys the general solution decaying at the infinity according to (8) For further use we note the following differentiation rules similar results hold for φ. To reduce the number of unknown coefficients, we insert (29) and (26) into the tangential stressfree condition (10) and use (30) to calculate Inserting (29) and (31) into (25) and (26) and for the nonzero components of the stress To treat the normal stress-free condition T zz (t, r, 0) = 0 for r > h in (12) , we represent the free factor C with the help of the cosine Fourier transform as by means of unknown functions p ( t ) and q ( l ) (see Sneddon, 1960 ). Indeed, substituting (34) into (33) , interchanging the order of integration, using integration by parts for ξ cos (lξ ) = d dl sin (lξ ) and the discontinuous Weber-Schafheitlin integral (see Abramowitz and Stegun, 1972 , 11.4.39) Inserting (34) into (32) and using the other discontinuous Weber-Schafheitlin integral and from (14) obtain the representation for the normal displacement u z ( t, r, z ) at the boundary z = 0 as Finally, applying (38) to the normal displacement in (11) , we derive the resulting integral equation for determining the unknowns from the given contact force due to (9) −2 π h 0 T zz (t, r, 0) r dr = F (t) . (40)

The main result
Gathering together representations of Section 2 we restate the indentation problem (5) -(12) for the flat-ended axisymmetric punch. Given the punch radius h > 0 and the time-dependent normal force F ( t ) ≥ 0 with rate ˙ F (t) , for all radii r ≥ 0, depths z ≥ 0, and times t ≥ 0 find the punch indentation depth D ( t, h ), displacement u r ( t, r, z ), u z ( t, r, z ), and stress T rr ( t, r, z ), T θθ ( t, r, z ), T rz ( t, r, z ), T zz ( t, r, z ) in the viscoelastic body which satisfy: the balance equations T rr,r + T rz,z + 1 r the constitutive equations the linearized strain equations ε rr = u r,r , ε θθ = 1 r u r , ε zz = u z,z , ε rz = 1 2 (u r,z + u z,r ) ; and the asymptotic condition at infinity u r , u z → 0 as r 2 + z 2 → ∞ .
As z = 0 , the normal stress T zz fulfills (40) , and the following boundary conditions hold: the tangential stress-free condition T rz (t, r, 0) = 0 (45) and the complementary contact conditions Theorem 1 (Flat-punch indentation) . The flat-punch indentation problem (40) -(47) has the solution expressed in the analytical form as follows: the indentation depth is given by the displacement is given by and the stress is given by The integrals J m n (h ) for integer m, n ≥ 0 are defined by the imaginary part of the mth-order Hanckel transform J m which is calculated analytically for (50) as The proof of Theorem 1 is given in Appendix A .
The computed solutions (49) and (50) for an example indentation problem are depicted versus ( r, z )-coordinates in Fig. 2 .
Here the material parameters are ν = 0 . 5 and rather small E = 1 MPa that is typical for soft materials. The flat-ended punch of radius h = 1 mm is indented with the constant loading rate ˙ F = 10 N/s subject to the relations such that after t = 1 s the indentation depth D ( t, h ) ≈ 1.38 mm. In Fig. 2 we can observe the boundary behavior of the solution at z = 0 , and the discontinuous stress at the punch end-point (r, z) = (h, 0) .
On passing the punch radius h → 0 we get the following Boussinesq problem subject to a concentrated load: find the displacement u r ( t, r, z ), u z ( t, r, z ) and the stress T rr ( t, r, z ), T θθ ( t, r, z ), T rz ( t, r, z ), T zz ( t, r, z ) in the viscoelastic body which satisfy: the balance Eq. (41) ; the constitutive equations (42) ; the linearized strain (43) ; the asymptotic condition at infinity (44) ; the tangential stress-free condition (45) and the Neumann-type boundary condition for the normal stress at z = 0 The Dirac delta-function δ (Dirac measure) in (54) implies that such that in virtue of (54) and (55) −2 π ∞ 0 T zz (t, r, 0) r dr = F (t) .
We note that, for the indentation problem, the normal stress (40) and the relation T zz (t, r, 0) = 0 for r > h in complementary conditions (47) together lead to the equality (56) . From Theorem 1 we derive the next result.
Theorem 2 (Boussinesq problem and the stress The integrals here are defined for integer m, n ≥ 0 The proof of Theorem 2 is given in Appendix B . The computed solutions (57) and (58) of the corresponding Boussinesq problem are depicted in Fig. 3 with the parameter set from the previous example (see (53) ). Here we clearly observe the singular behavior of displacement and stress at the origin r = z = 0 causing evidently numerical instabilities. We recall that Hooke's law (3) is provided by the identity where (53) is a particular case. Thus, we obtain the following result. (57) , (58) in Theorem 2 hold true in the following cases: for an elastic body given by Hooke's constitutive relation when inserting (61) in these relations, and for the Kelvin-Voigt model by substituting (62) in these relations .

Corollary 1. The analytical formulas (48) -(50) in Theorem 1 and
The corollary suggests that though the constitutive relation we started with is not invertible and even the domain and boundary conditions are not appropriate for the correspondence principle to hold, we see a connection between the elastic and viscoelastic solutions.

Appendix A. Proof of Theorem 1
We start with the separation of variables in (39)