Advances in Tribology

Volume 2018, Article ID 9432894, 16 pages

https://doi.org/10.1155/2018/9432894

## Viscoelastic Contact Simulation under Harmonic Cyclic Load

^{1}Department of Mechanics and Technologies, Stefan cel Mare University of Suceava, 13th University Street, 720229 Suceava, Romania^{2}Integrated Center for Research, Development and Innovation in Advanced Materials, Nanotechnologies, and Distributed Systems for Fabrication and Control (MANSiD), Stefan cel Mare University of Suceava, Suceava, Romania

Correspondence should be addressed to Sergiu Spinu; or.vsu.mif@unips.uigres

Received 29 January 2018; Accepted 16 April 2018; Published 20 May 2018

Academic Editor: Huseyin Çimenoǧlu

Copyright © 2018 Sergiu Spinu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Characterization of viscoelastic materials from a mechanical point of view is often performed via dynamic mechanical analysis (DMA), consisting in the arousal of a steady-state undulated response in a uniaxial bar specimen, allowing for the experimental measurement of the so-called complex modulus, assessing both the elastic energy storage and the internal energy dissipation in the viscoelastic material. The existing theoretical investigations of the complex modulus’ influence on the contact behavior feature severe limitations due to the employed contact solution inferring a nondecreasing contact radius during the loading program. In case of a harmonic cyclic load, this assumption is verified only if the oscillation indentation depth is negligible compared to that due to the step load. This limitation is released in the present numerical model, which is capable of contact simulation under arbitrary loading profiles, irregular contact geometry, and complicated rheological models of linear viscoelastic materials, featuring more than one relaxation time. The classical method of deriving viscoelastic solutions for the problems of stress analysis, based on the elastic-viscoelastic correspondence principle, is applied here to derive the displacement response of the viscoelastic material under an arbitrary distribution of surface tractions. The latter solution is further used to construct a sequence of contact problems with boundary conditions that match the ones of the original viscoelastic contact problem at specific time intervals, assuring accurate reproduction of the contact process history. The developed computer code is validated against classical contact solutions for universal rheological models and then employed in the simulation of a harmonic cyclic indentation of a polymethyl methacrylate half-space by a rigid sphere. The contact process stabilization after the first cycles is demonstrated and the energy loss per cycle is calculated under an extended spectrum of harmonic load frequencies, highlighting the frequency for which the internal energy dissipation reaches its maximum.

#### 1. Introduction

Important engineering applications involving products like automotive belts and tires, seals, or biomedical devices require accurate prediction of tribological processes between viscoelastic materials such as elastomers or rubbers. Considering that a closed-form description of the viscoelastic contact is difficult to achieve due to complexity of the emerging equations, numerical simulation presents itself as a worthy substitute, capable of assisting the design of tribologically competent products using viscoelastic materials.

The classic method for solving the linear viscoelastic problems of stress analysis is based on the concept of associated elastic problem [1, 2]. This approach involves removal of time dimension via Laplace transform, thus reducing the viscoelastic problem to a formally identical elastic problem whose solution is easier to achieve. Although application of this method is limited as the transient boundary conditions encountered in most contact problems cannot be directly treated by means of Laplace transform, solutions of limited viability were successfully obtained. Lee and Radok [3] derived the solution of a Hertz-type problem involving linear viscoelastic materials for the case when the contact area increases monotonically with time. Hunter [4] solved the same problem for monotonic contact radius or when the radius possesses a single maximum. The contact problem of viscoelastic bodies was extended by Yang [5] to cover general linear materials and arbitrary quadratic contact geometry. A more versatile solution, allowing for any number of loadings and unloadings, in which the contact area is a simply connected region, was later achieved by Ting [6]. A further iteration [7] involves multiple connected contact regions and contact radii described by arbitrary functions of time. The contact problem between an axisymmetric indenter and a viscoelastic half-space was recently revisited by Greenwood [8]. The mathematical complexity of these partially analytical solutions challenges their wide range application, especially in the case of the contact process under harmonic cyclic load, when up to five different cases [7] have to be considered due to the specifics of contact radius dependence on the history of the loading program. Algorithmization may also be problematic as the resulting implicit solutions require numerical integration and differentiation, as well as the resolution of transcendental equations, which may raise convergence issues.

Recent developments in numerical resolution of elastic contact problems encouraged a new approach to the problem of viscoelastic contact. Kozhevnikov et al. [9] advanced a new algorithm for the indentation of a viscoelastic half-space based on the Matrix Inversion Method (MIM). Chen et al. [10] developed a new semi-analytical method (SAM) for contact modeling of polymer-based materials with complicated properties and surface topography. These authors [11] studied the multi-indentation of a viscoelastic half-space by rigid bodies using a two-scale iterative method (TIM).

Spinu and Gradinaru [12] advanced a semi-analytical method for the computation of displacement in linear viscoelastic bodies subjected to arbitrary surface tractions. A solution for the Cattaneo-Mindlin problem involving viscoelastic materials was advanced by Spinu and Cerlinca [13], based on the algorithm for the frictionless viscoelastic contact problem reported in [14]. More recently, the same authors studied [15] the rough contact of viscoelastic materials by imposing a simplified manner in which plasticity effects at the tip of the asperities are accounted for. While dealing with friction or with surface microtopography, these contact simulations employ simple loading programs consisting usually in step or ramped loadings. The contact process under cyclic harmonic load is investigated in this paper based on an extended version of the viscoelastic contact algorithm advanced in [14]. New algorithm developments allow for more general boundary conditions, involving displacement driven contact scenarios, as well as the assessment of the post-unloading contact state.

From the point of view of algorithmic complexity, the semi-analytical method (SAM) employed herein, derived from the boundary element method (BEM), has significant computational advantages over the finite element method (FEM), as it requires a 2D spatial discretization only (i.e., the meshing of the potential contact surface), whereas a FEM simulation entails the 3D meshing of the entire contacting bodies. According to this review [16], the computational efficiency of SAM greatly exceeds that of FEM; for example, a 3D SAM contact simulation can be conducted with a computational effort comparable to that of a 2D finite element analysis. In this paper, the algorithmic computational efficiency is optimized by employing state-of-the-art numerical techniques: the conjugate gradient method (CGM), with its superlinear rate of convergence, is used for the resolution of the emerging linear system of equations, whereas the Discrete Convolution Fast Fourier Transform (DCFFT) technique [17] is engaged in the rapid computation of discrete convolution products.

#### 2. Viscoelastic Constitutive Law and Associated Contact Solutions

In the framework of linear theory of viscoelasticity [18], the material exhibits a linear stress-strain relationship; that is, an increase in stress by a constant factor leads to an increase in the strain response by the same factor. The response functions to excitations conveyed by Heaviside step functions are referred to as the material functions of the viscoelastic body, namely, the creep compliance and the relaxation modulus, which are both functions of time . The creep compliance function describes the viscoelastic strain response to a unit step change in stress, and the relaxation modulus , conversely, describes the stress response to a unit step change in strain. With these functions, the linear viscoelastic response to various sequences of stress or strain is assessed, according to the Boltzmann hereditary integral, by either of the two Volterra integral equations:where and are the tensors of deviatoric stress and deviatoric strain, respectively. Consequently, (1) can be regarded as the superposition of a series of loading histories consisting in infinitesimal changes in strain or stress, respectively, applied separately in a window of observation , chosen so that initially (i.e., at ) the viscoelastic body was undisturbed.

Analog mechanical models, constructed from linear springs and dashpots, arranged in series or in parallel, are convenient tools to model the linear viscoelastic behavior under uniaxial loading. The combination rules for these basic units state that creep compliances add in series and relaxation moduli in parallel.

The ideal spring, also referred to as the Hooke model or the ideal solid, is the elastic element in which the force is proportional to the extension. By identifying force with stress and elongation with strain and according to Hooke’s law, , with being the pertinent (i.e., longitudinal or shear) elasticity modulus. The dashpot, also referred to as the Newton model or the perfect liquid, is the viscous element in which the force is proportional to the rate of extension. According to the Newton equation, , where is the rate of strain and is the viscosity coefficient. Both Hooke and Newton models represent limiting cases of viscoelastic bodies.

A branch constituted by a spring in parallel with a dashpot is known as the Kelvin-Voigt model, whereas a branch constituted by a spring in series with a dashpot is known as the Maxwell model. The differential equation for the Maxwell model can be expressed as [18] for , under the assumption that for . The creep compliance function for the viscoelastic half-space described by a Maxwell model consisting in a spring of shear modulus in series with a dashpot of viscosity results as [19] , where denotes the relaxation time: .

Ting’s formalism [6, 7] can be employed to describe the indentation of the Maxwell viscoelastic half-space by a rigid spherical punch of radius in a step loading , where denotes the Heaviside step function and is fixed but otherwise arbitrarily chosen. The equation for the pressure distribution achieved at time after the first contact, at the radial coordinate , results aswhere denotes the contact radius at time :and denotes the real part of its complex argument. This partially analytical solution is a more computationally friendly form of the corresponding equation derived in [20].

The Maxwell model accounts well for relaxation but handles badly both creep (model creeps without bound at constant rate; therefore it is also referred to as the Maxwell fluid) and recovery (only elastic deformation is recovered, and this is done instantaneously). The Kelvin-Voigt model, on the other hand, handles creep and recovery fairly well but does not account for relaxation. Moreover, the latter model exhibits no instantaneous elastic response; consequently, the elasticity modulus is formally infinite, and contact pressure results to be infinite at the contact boundary in the beginning of the loading process, as demonstrated by Ting’s model [6, 7] implementation reported in [19]. Therefore, it can be asserted that the assumptions of the Kelvin model make it inappropriate for contact analyses.

The Maxwell and Kelvin models are adequate for qualitative or conceptual analyses, but quantitative representation of the behavior of real materials requires an increase in the number of parameters. The generalized Maxwell model is composed of a number of Maxwell models and an isolated spring in parallel, whereas the generalized Kelvin model consists in a number of Kelvin units plus an isolated spring in series. The Standard Linear Solid Model, also known as the Zener model, can be represented as a spring of shear modulus in series with a Kelvin model of parameters and or as a spring in parallel with a Maxwell model. By adopting the former representation, the differential equation for the Zener model results as [18]yielding the following creep compliance function:

The Zener model exhibits instantaneous elastic strain when stress is instantly applied; if the stress is held constant, the strain creeps towards a limit, whereas, under constant strain, the stress relaxes towards a limit. Moreover, when stress is removed, instantaneous elastic recovery occurs, followed by gradual recovery towards vanishing strain. Using Ting’s formalism [6, 7], the pressure distribution in the step loading of a Zener half-space by a rigid spherical punch of radius results aswith being the contact radius given by (3).

These basic models, having only one relaxation time, are capable of providing only qualitative description of viscoelastic behavior, whereas precise quantitative assessments require more parameters, related to the naturally occurring spectrum of relaxation times of the real viscoelastic material. Such a goal can be accomplished by using a complex model such as the generalized Wiechert model, which consists in several Maxwell units and a free spring, connected in parallel. The shear relaxation modulus function of the Wiechert model can be expressed as [18]where is the long-term modulus (longitudinal or shear) once the material is totally relaxed and and , with , are the relaxation time and the spring stiffness of each Maxwell subunit. The naturally occurring spectrum of relaxation times of a viscoelastic material can be described by including as many exponential terms as needed. Relation (7) is also referred to as a Prony series. The Prony series of a viscoelastic material is usually obtained by a one-dimensional relaxation test, in which the viscoelastic material is subjected to a sudden strain that is kept constant, while measuring the stress response over time. The initial stress is related to the purely elastic response of the material. Later on, the stress relaxes due to the viscous effects in the viscoelastic material. Mathematical description employing the Prony series can be achieved by fitting the experimental data to (7) by adjusting the model parameters , , and .

The constitutive law of the real viscoelastic material considered in this paper is that of the polymethyl methacrylate (PMMA), a thermoplastic polymer whose mechanical properties were studied extensively by Ramesh Kumar and Narasimhan [21]. These authors obtained experimentally the PMMA relaxation modulus data under uniaxial compression in a window of observation of 1000 s. Based on their results, the two-term Prony series of PMMA results as [10]which is the modulus relaxation function of the material (expressed in terms of longitudinal modulus), depicted in Figure 1.