Plane Thermoelastic Waves in Infinite Half-Space Caused
FACTA UNIVERSITATIS
Series: Mechanical Engineering Vol. 14, N
o
3, 2016, pp. 293 - 300
DOI: 10.22190/FUME1603293K
Original scientific paper
THE BOUNDARY ELEMENT METHOD FOR VISCOELASTIC
MATERIAL APPLIED TO THE OBLIQUE IMPACT OF SPHERES
UDC 539.3
Stephan Kusche
Department of System Dynamics and the Physics of Friction, TU Berlin, Germany
Abstract. The Boundary Element Method (BEM) for elastic materials is extended to deal
with viscoelastic media. This is obtained by making use of a similar form of the
fundamental solution for both the materials. Some considerations are attributed to the
difference of the normal and the tangential contact problem. Both normal and tangential
problems are furthermore assumed to be decoupled. Then the oblique impact of hard
spheres with an incompressible viscoelastic half-space (linear standard-model) is studied.
By assuming stick conditions during impact, one obtains the dependence of the two
coefficients of restitution as functions of two input parameters. This result is expressed in
an elegant and compact form of the fitting function.
Key Words: Contact Mechanics, Boundary Element Method, Linear Viscoelastic
Material, Rebound Test, Oblique Impact, Tangential Problem
1. INTRODUCTION
The history of analytical research of the spheres’ field of the impact started with Hertz
[1] in 1882. His theory of normal contacts has been extended to tangential contact theories,
considering partial slip [2-9]. Like in other fields of contact mechanics, a solution in the
case of viscoelastic time dependent material behavior is difficult and its closed form has
not been found yet. For the pure normal impact of spheres Hunter [10], Sabin [11] and
Calvit [12] discovered some analytical solutions. These solutions are given in an implicit
form; some kind of numerical treatment is needed in order to obtain results.
In the field of elastic contact mechanics, the Boundary Element Method (BEM) is well
established. This method is based on the fundamental solution for elastic half-space. Since
the fundamental solution for the viscoelastic half-space is known as well, it should be possible
to adapt many procedures from elastic to viscoelastic problems. This can be achieved by using
Received September 22, 2016 / Accepted November 18, 2016
Corresponding author: Stephan Kusche
Institute of Mechanics, Berlin Institute of Technology, Str. d. 17. Juni 135, 10623 Berlin, Germany
E-mail: s.kusche@tu-berlin.de
294 S. KUSCHE
the Prony series for modeling the material behavior. In this paper the given approach is
used to study the viscoelastic impact problem of spheres.
2. PROBLEM DESCRIPTION
The problem under consideration is the following: A sphere of mass m and radius R has
initial velocities vx0
s
, vz0
s
and an initial angular velocity ω0
s
. After rebound from viscoelastic
ground the velocities change to vx1
s
, vz1
s
and ω1
s
(see Fig. 1). These velocities have to be
calculated.
To simplify this problem approach it will be assumed, that the normal- and the
tangential-contact problems are decoupled. This assumption becomes true under two
circumstances. Firstly, the sphere is a rigid body and the viscoelastic half-space is
incompressible. Then the normal load will not lead to any tangential displacement and vice
versa. Secondly, the displacement of the centre of the sphere in tangential direction ux
s
is
assumed to be much smaller than radius R of the sphere. Then only tangential contact force
Fx has to be taken into account for the calculation of the resulting moment acting on the
sphere. Furthermore, the contact area will conserve a rotational symmetric shape and the
action line of resulting normal force Fz points toward the centre of the sphere.
Fig. 1 Velocities before and after impact (left hand side). Displacement of the center
of the sphere and resulting contact forces during impact (right hand side)
With these considerations tangential velocity vx
s
and angular velocity ω
s
are not
independent of each other. By using the principle of impulse and angular momentum one
obtains:
22
,
5
s s
x x x
mu F mR F R , (1)
and this results in:
0 0
2
( )
5
s s s s
x x
v v R . (2)
This means that the after-impact velocities can be expressed by two numbers:
1 1 11
0 0 0 0
,
s s ss
x kz
z xs s s s
z x k
v R vv
e e
v v R v
. (3)
The Boundary Element Method for Viscoelastic Material Applied to the Oblique Impact of Spheres 295
For the normal contact problem the coefficient of restitution ez is used. In tangential
direction the velocity of the contact point is the major variable. The ratio of the contact
point’s velocity after and before impact ex is the second number used in this paper. The
tangential and angular velocity after the impact can be calculated from this number.
The viscoelastic half-space is modeled by the linear standard-model. The Maxwell-Element
and the Kelvin-Voigt-Element are included as special cases. For the linear standard-model the
time dependent relaxation modulus to shear has the following form:
( ) exp
t
G t G G
. (4)
It turns out that the numerical solution for ex and ez can be expressed by the following
two dimensionless numbers:
1/5 1/5
0 0
1 22 3 2 3
, ,
s s
z z
Rv Rv
G
m G m G
. (5)
At first glance it seems to be unreasonable that neither the tangential nor the angular
velocity appears therein. But, under consideration that normal and tangential contact
problems are decoupled, at least the coefficient of restitution ez cannot be dependent on any
tangential velocity. Then again, the tangential problem is, under the assumption of a full
stick regime, only dependent on the contact point’s velocity and the contact area’s shape
during the impact. Since the system is linear, the magnitude of the contact point’s velocity
has no significance. It is finally used for normalization of ex, as this has been done in the
definition stated before.
Both numbers, δ1 and δ2, can be adapted to the special cases of the material model. The
Maxwell-Element is obtained if G∞ tends towards zero, which means that δ2 becomes
infinite. The Kelvin-Voigt-Element is obtained if G tends toward infinity but η is held
constant. Since δ1~ηG
-3/5
and δ2~η, the parameter δ1 becomes zero. Finally in the elastic
case η can be obtained to be either zero or infinity, which correspondents to an infinite soft
or an infinite stiff damper. This means, expressed by the parameters, that either δ2 becomes
zeros (and δ1 is finite, the elastic constant is G∞) or δ1 tends towards infinite (and δ2 is finite,
the elastic constant is then equal to G∞+G). Summarized this means:
2
1
1 2
Maxwell-Element:
Kelvin-Voigt-Element: 0
Elastic Spring: or 0
(6)
3. THE BOUNDARY ELEMENT METHOD
The used BEM is based on the fundamental solution for a point force acting on a
viscoelastic, incompressible half-space [13-15]. Let x and y be Cartesian coordinates on the
surface of a half-space. A point force with components Fx, Fy and Fz, is applied at the centre
of the mentioned coordinate system and causes the following deflection of the surface:
296 S. KUSCHE
2
2 2
3 3
2 1
4 4,
( ),
,z
x x x y z
i i
x xy
u u
r r rr
J t r xF y
r
. (7)
These equations can only be applied, if the acting force is added at time zero and the
force is held to be constant. From this solution, the elastic fundamental solution can be obtained
if J(t)→1/G is substituted. Herein J(t) is the time dependent creep function for shear:
/
0
( ) ( ) 1
tt
J t J J e
. (8)
This creep function (8) is related to the time dependent relaxation modulus G(t) given
in Eq. (4). In the case of the standard-material:
0 0
0
1
, , ,
G
J J J
G G G J G
, (9)
and in case of the Maxwell-element:
0
1
, 0, J J G
G
. (10)
The fundamental solution for the viscoelastic material looks similar to the one from
elastic material. Especially the dependency in space is the same. With this conclusion many
schemes from elastic contact mechanics can be used. Firstly, a standard procedure is to
integrate the fundamental solution (Boussinesq) over a rectangle, assuming constant pressure
[16]. This solution (Love) can be superimposed to find the deflection for an arbitrary, but
piecewise constant, pressure distribution [17]. To speed this up, it is common practice to
interpret this procedure as a convolution, which can be carried out on massive parallel
processing units [18-20]. The next, more complicated problem arising in the field of elastic
contact mechanic is the opposite task: Solve the indentation problem of a rigid indenter
pressed in an elastic half-space with an a priori unknown contact area. Then conjugate
gradient algorithms are used, like the modified form by Polonsky and Keer [21]. The
procedures named above are supposed to be known and are used in the following calculations.
The method described above can be applied to the viscoelastic case, if the acting forces
are kept constant during one time step. The overall solution will be superimposed by
adding the elastic solution multiplied with the time dependent creep function delayed by
the time span since the force has been applied. At this point, a sum appears whose numerical
costs grow linear in simulation time. By utilizing the special form of the creep function used
here, an iterative algorithm can be derived. By applying this algorithm, only the previous time
step is needed, and the contacting indenter shape is modified by a viscoelastic term. By using
this method, the normal contact can be solved. A description and the application to the
tangential sliding of an indenter can be found in [22]. The algorithm is shown below:
The Boundary Element Method for Viscoelastic Material Applied to the Oblique Impact of Spheres 297
1 1
1
, 1 1
0 00 0 0
, 1 1 1
0 0 0
1
1
( ) exp
ex
1
)p
( )
(
t
z z
n n
z z
n n
h
n n
n i
z n n i i i i i
i i
a b
z zt
z n n n t n n n n
a b
n n
t tJ J
f t t f f f
J J J
hJ J
f a h f b f f
J J J
f a
u
u
1
0
exp
z
z zt
n n
D
h J
b f
J
. (11)
Herein uz,n is the normal deflection of the surface and fn is the deflection due to a pressure
distribution pn - each at time tn. Furthermore, a and b are variables from the previous time step,
initialized with value zero at the beginning of the simulation, and it is J∞=J0+J. As can be seen
in the last line of Eq. (11), deformation D
z
has to been taken into account to include viscoelastic
behavior (in the elastic case D
z
is equal to zero). The only unknown term in Eq. (11) is fn+1,
which means that pressure distribution pn+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 for instant deformation (G(0)<∞ and J(0)>0). That excludes the
Kelvin-Element because it does not fit into the scope of this algorithm.
The tangential contact is very similar to the normal one. Only the calculation of the
deflection in tangential direction ux has to be adopted. In all simulations carried out here,
the deformation - perpendicular to the plane of the motion - is neglected. It has been shown
for the case of parabolic bodies that this assumption causes only a negligible error [23,24].
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 step size ht the solution of the elastic impact
problem has been used:
1
2 5
4
, 2
0
10 2.87
(4 ( ))
t c c elastic s
z
m
h T T
R G t v
. (12)
Since the shear modulus is time dependent, one has to start with G(0) and iteratively set
t→Tc. With the elastic solution also an approximation for contact radius a can be found:
, 0
,
2.94
s
c elastic z
elastic
T v
a a R . (13)
The geometric discretisation has been chosen in a way that 256256 points are used. A
comparison with a finer discretisation shows only a slight improvement of the error. 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 the previous time steps. This
can be handled by adding only the current increment of tangential movement at the boundary
of the sphere in each time step.
298 S. KUSCHE
4. RESULTS
The simulation carried out shows that the assumptions, condensed in formula (3) and
(5), are correct. This has been proven by choosing random simulation parameters, that all
coincide in the same functional relation. This relation has been fitted to the numerical data
with two functions, one for each number ex and ez. In the following contour plots (Figs. 2
and 3) the obtained data and the fitting are shown. The fitting function and the parameter
are given below.
The fitting function for the coefficient of restitution in normal direction ez is (parameters are
given in Table 1):
1 2
2 2
10 1 1
10 2
2
1 1
tanh( )
2 2
max 0, log ( )
log ( )
max 0,
2
z
p
p
c c
x y
x
y
e
(14)
Fig. 2 Contour lines of the coefficient of restitution ez (left hand side). Special case of the
Kelvin-Voigt-Element (δ1→0) and the Maxwell-Element (δ2→∞) (right hand side)
Table 1 List of the fitting parameters for the coefficient of restitution in normal direction ez
corresponding to the fitting function given in (14)
Parameter δ1
p
δ2
p
c1 c2
Value -1.342 0.9583 2.079 -2.892
The Boundary Element Method for Viscoelastic Material Applied to the Oblique Impact of Spheres 299
The fitting function for the coefficient of restitution in tangential direction ex is (parameters
are given in Table 2):
2
1
3 4
2
2
5 6
7
8 9
10 1 1
10 2 2
exp
( )
( ) ( ) ,
[1 tanh( )] ,
2
1
1 tanh( ) 2 exp 2
4
1 t
( )
anh( ) ,
log
log
el min el max el
x x x x x x x x
x
x
p
p
e e e e e e e e
c x
e c c y
c
y
e c c y
c
c c x
x
y
(15)
Fig. 3 Contour lines of the coefficient of restitution ex (left hand side). Special case of the
Kelvin-Voigt-Element (δ1→0) and the Maxwell-Element (δ2→∞) (right hand side)
Table 2 List of the fitting parameters for the coefficient of restitution in normal direction
ex corresponding to the fitting function given in Eq. (15)
Parameter ex
min
ex
el
ex
max
δ1
p
δ2
p
Value -0.0399 0.1038 0.2028 -0.4 -0.8
Parameter c1 c2 c3 c4 c5
Value 2.335 0.5406 -2.29 2.253 1.718
Parameter c6 c7 c8 c9
Value 2.335 0.7999 -2.105 4.000
300 S. KUSCHE
5. CONCLUSIONS
The BEM has been extended to the normal and tangential contact problem of viscoelastic
media. The approach has different advantages: It implements a relatively simple adaptation
to the existing BEM procedures developed for elastic media and its applicability to any other
transient problem dealing with time dependent viscoelastic behavior modeled by the Prony
series. Exemplary the impact of a sphere with an elastomer, modeled by a Maxwell-element
and a standard-model, has been investigated. The obtained compact and elegant solution is
given in form of a fitting function for further use.
REFERENCES
1. Hertz, H., 1882, Über die Berührung fester elastischer Körper, Journal für die reine und angewandte
Mathematik, 92, pp. 156-171.
2. Deresiewicz, H., 1968, A Note on Hertz’s Theory of Impact, Acta Mechanica, 6, pp. 110–112.
3. Mindlin, R.D., 1949, Compliance of Elastic Bodies in Contact, Journal of Applied Mechanics, 16, pp. 259–268.
4. Maw, N., Barber, J.R., Fawcett, J.N., 1976, The Oblique Impact of Elastic Spheres, Wear, 38(1), pp. 101–114.
5. Maw, N., Barber, J.R., Fawcett, J.N., 1977, The Rebound of Elastic Bodies in Oblique Impact, Mechanical
Research Communications, 4(1), pp. 17–22.
6. Maw, N., Barber, J.R., Fawcett, J.N., 1981, The Role of Elastic Tangential Compliance in Oblique Impact,
Journal of Lubrication Technology, 103(1), pp. 74–80.
7. Jäger, J., 1994, Analytical Solutions of Contact Impact Problems, Applied Mechanics Review, 47(2), pp. 35–54.
8. 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, Sci. Rep., 5, 8479.
9. 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, Journal of
Applied Mathematics and Mechanics (ZAMM), 96(9), pp. 1089-1095.
10. 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.
11. Sabin, G. C. W., 1987, The impact of a rigid axisymmetric indenter on a viscoelastic half-space, International
Journal of Engineering Science, 25(2), pp. 235–251.
12. Calvit, H. H., 1967, Numerical solution of the problem of impact of a rigid sphere onto a linear viscoelastic
half-space and comparison with experiment, International Journal of Solids and Structures, 3(6), pp. 951–960.
13. 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.
14. 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.
15. 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.
16. Johnson, K. L., 1985, Contact Mechanics, Cambridge University Press, Cambridge.
17. Pohrt, R., Li, Q., 2014, Complete Boundary Element Formulation for Normal and Tangential Contact
Problems, Physical Mesomechanics, 17(4), pp. 334-340.
18. 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.
19. 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.
20. 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, 145–154.
21. Polonsky, I., Keer, L., 1999, A numerical method for solving rough contact problems based on the multi-level
multi-summation and conjugate gradient techniques, Wear, 231(2), 206–219.
22. Kusche, S., 2016, Frictional force between a rotationally symmetric indenter and a viscoelastic half-space.
ZAMM - Journal of Applied Mathematics and Mechanics, pp. 1-14.
23. Johnson, K.L., 1955, Surface interaction between elastically loaded bodies under tangential forces, Proc. R.
Soc. A., 230, 531-548.
24. Munisamy, R.L., Hills, D.A., Nowell, D., 1994, Static axisymmetric Hertzian contacts subject to shearing
forces, ASME J. Appl. Mech., 61(2), pp. 278–283.