# Modelling of Rough Contact between Linear Viscoelastic Materials

Authors: Sergiu Spinu; Delia Cerlinca
Journal: Modelling and Simulation in Engineering
 (2017)

---

## Abstract

The important gradients of stress arising in rough mechanical contacts due to interaction at the asperity level are responsible for damage mechanisms like rolling contact fatigue, wear, or crack propagation. The deterministic approach to this process requires computationally effective numerical solutions, capable of handling very fine meshes that capture the particular features of the investigated contacting surface. The spatial discretization needs to be supported by temporal sampling of the simulation window when time-dependent viscoelastic constitutive laws are considered in the description of the material response. Moreover, when real surface microtopography is considered, steep slopes inevitably lead to localized plastic deformation at the tip of the asperities that are first brought into contact. A computer model for the rough contact of linear viscoelastic materials, capable of handling deterministic contact geometry, complex viscoelastic models, and arbitrary loading histories, is advanced in this paper. Plasticity is considered in a simplified manner that preserves the information regarding the contact area and the pressure distribution without computing the residual strains and stresses. The model is expected to predict the contact behavior of deterministic rough surfaces as resulting from practical engineering applications, thus assisting the design of durable machine elements using elastomers or rubbers.


---

## Body

## 1. Introduction

The design of mechanical components requires understanding of tribological phenomena such as friction, wear, or contact fatigue. When two machine elements are brought into contact and load is transmitted, interfacial tractions occur due to direct interaction of the two surfaces. As real surfaces are not smooth at microscopic levels, interaction firstly occurs at discrete contact spots, leading to important gradients of contact pressure and of subsurface stress, often causing localized plastic deformation. The knowledge of these contact stresses provides the foundation for the investigation of surface-related phenomena, such as rolling contact fatigue, wear, or crack propagation.Due to their structural complexity and time-dependent properties, viscoelastic materials, such as elastomers or rubbers, are extensively used in engineering applications involving automotive belts and tires, seals, or biomedical devices. The closed-form mathematical description of the contact behavior of viscoelastic materials is still in its early stages, especially when considering complicated rheological models for the material response, sporadic contact regions at the asperity level, or complex loading histories. A numerical model can advance the understanding of the viscoelastic contact in the presence of asperity interaction and thus assist the design of durable machine elements using viscoelastic materials.Existing analytical efforts are based on the elastic-viscoelastic correspondence principle, capitalizing on the fact that a viscoelastic contact problem formally reduces to a so-called associated elastic problem, after removal of the time dimension via Laplace transform. Assuming the boundary conditions are time-independent, the viscoelastic solution in the frequency domain is identical in form to that of the associated elastic problem. As the latter can be derived more easily, the desired viscoelastic solution results by correspondence in the frequency domain and subsequently in the time/space domain by means of inverse Laplace transform. However, this method has limited applicability to contact problems, which exhibit transient boundary conditions difficult to treat by means of the Laplace transform. The viscoelastic contact solutions [1–3] obtained by manipulation of existing elastic solutions for the associated Hertz contact problem are subjected to limitations concerning the monotony of the contact radius. These limitations were removed by Ting [4, 5], who developed an analytical formulation for the general contact problem involving arbitrary history of the contact area and axisymmetric contact geometry. The algorithmization of this solution may require up to five conditional branches for complicated loading histories and involves differentiation and repeated integration, as well as resolution of transcendental equations. The mathematical complexity of these (partially) analytical solutions generally discourages their wide range application; therefore, great research efforts have been recently dedicated to solving the problem numerically.The numerical treatment of the viscoelastic contact problem is supported by recent developments in computational contact mechanics, based on the so-called semianalytical methods (SAM) derived from the boundary element method. According to this review [6], the computational efficiency of SAM greatly exceeds that of the finite element method (FEM), to the point where a 3D SAM contact simulation can be conducted with the same computational effort as a 2D FEM analysis. SAM models are more suitable to solve 3D contact problems because only a small contact region needs to be meshed and considered in the simulation, whereas the meshing of the bulk is required in FEM. Successful applications of SAM in the field of viscoelastic contact mechanics were reported by these authors [7–10].The contact of real surfaces occurs at discrete spots due to inherent surface roughness, and the sum of the areas of all contact spots constitutes the real contact area. As rough surfaces are essentially random processes, stochastic techniques have been extensively used in modelling rough contact problems [11–13], allowing for a rapid prediction of the mean level of contact area and contact pressure. Stochastic models assume that surface roughness consists in many tiny asperities having a specific probability distribution (e.g., Gaussian) of heights and curvatures, which deform independently of each other. This very conservative assumption is unrealistic, especially at high loads. Detailed description of real contact area and pressure requires a deterministic approach, for which the numerical simulation, capable of analyzing the contact process involving measured mechanical properties and 3D roughness maps, is best suited. Whereas FEM simulation of 3D contact problems with real or computer generated surfaces proves to be very time-consuming, the discretization level allowed in SAM [14, 15] can predict accurate contact parameters with a moderate computational effort.The SAM model employed in this paper features an updated version of the contact solver advanced by Polonsky and Keer [16]. Initially developed for the purely elastic contact, this algorithm was subsequently implemented in simulation of path-dependent (though time-independent) processes like the elastic-plastic [17] or the slip-stick [18] elastic contact. This paper extends the model generality by iteratively computing the time-dependent contact process in small time steps, thus simulating the memory effect of the viscoelastic material. To this end, a recently developed [19] semianalytical formula for computation of viscoelastic displacement induced by a known but otherwise arbitrary distribution of surface tractions is employed. Moreover, a limit for pressure is imposed on the contact area to account for the nontrivial plastic deformation occurring at the asperity level. The computational efficiency of the algorithm is increased by using a well-established method [20] for rapid calculation of convolution products, allowing the grid resolution to be increased to the level of microtopography particular features.


## 2. Assumptions and Limitations

Computational contact mechanics can incorporate the theory of viscoelastic behavior provided a directly additive (i.e., linear) viscoelastic response is assumed. Linearity in the stress-strain relationship of the viscoelastic material can be achieved in the framework of infinitesimal strains and delivers a starting point for the mathematical modelling of viscoelasticity, allowing for the computation of responses to various sequences of stress or strain.The linear viscoelastic stress-strain relationship requires two sets of functions, describing the behavior in shear and in dilatation, respectively. As Poisson’s ratio of polymers usually exceeds 0.4, additional model simplification can be achieved by considering the material as incompressible. With the latter assumption, only behavior in shear is needed, and the constitutive law is reduced to a linear relation between the tensors of deviatoric stresss and of deviatoric strain e by means of the shear modulus G (i.e., s=2Ge). The linear viscoelastic material response to arbitrary sequences of stress or strain can then be assessed, according to the Boltzmann hereditary integral, by employing two interchangeable functions of time (both in shear), namely, the relaxation modulus Ψ(t) and the creep compliance Φ(t): (1)st=∫0tΨt-t′∂et′∂t′dt′,et=∫0tΦt-t′∂st′∂t′dt′.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. Consequently, (1) can be regarded as the superposition of a series of loading histories consisting in infinitesimal changes in strain or in stress, respectively, applied in a window of observation [0,t], chosen so that initially (i.e., at t=0) the viscoelastic body is undisturbed.The contact model is also subjected to assumptions and limitations. The problem is considered quasi-static (i.e., the inertia forces due to deformation are negligible). The half-space approximation authorizes the use of the Green functions for the elastic half-space (i.e., the Boussinesq solution) in computation of displacement, together with the Boltzmann hereditary integral. This calculation assumes a linear viscoelastic response only, although nontrivial plastic deformation may occur at the asperities tips, especially until the contact area becomes sufficiently large and the elementary pressure can hold the applied load. Plasticity was successfully integrated in the semianalytical resolution of the elastic-plastic contact problem starting with the work of Jacq et al. [21]. Based on Betti’s reciprocal theorem, the elastic-plastic contact model can be divided into an elastic and a residual subproblem, whose mutually dependent solutions are obtained in an iterative manner: elastic stresses computed in the elastic subproblem are employed to assess the residual plastic strains using a universal algorithm [22] for integration of elastoplasticity equations, whereas residual displacement (i.e., displacement induced by the plastic strains) is superimposed to its elastic counterpart to assess the contact pressure in the elastic subproblem. The latter displacement requires the solution of the inclusion problem [23] for arbitrary eigenstrains in an elastic matrix, which can be obtained efficiently as shown in [24]. Despite being originally applied to rather coarse meshes, the SAM-based elastic-plastic model is conceptually functional in case of rough contact problems, but the computational complexity is generally considered prohibitive, especially for very fine meshes that intend to capture specific features of surface microtopography. A recent work [25] extended the solution of the inclusion problem to viscoelastic materials, by applying Eshelby’s formalism [26] to an ellipsoidal inhomogeneity located in a viscoelastic matrix. On the other hand, the relation between the viscoelastic stresses and the plastic flow inside the viscoelastic material, as well as the extension of Eshelby’s solution to a cuboidal inhomogeneity in a viscoelastic matrix, remains to be addressed.The development of a fully fledged contact model incorporating both viscoelastic effects and plasticity is beyond the scope of this paper. In the framework proposed here, the plastic deformation inside the asperities is accounted for by using the hardness of the viscoelastic material as a cut-off value for the contact pressure resulting from the viscoelastic model. This assumption, based on the Tabor equation [27], has been largely utilized [28, 29] in the modelling of the elastic-plastic contact of rough surfaces.The most severe model limitation consists in the neglect of adhesion, which appears virtually in all contacts scenarios. The force of adhesion is often trivial in case of metallic materials, when the real contact area, established between the asperity heights, is much smaller than the apparent (i.e., between the topographically smooth bodies) contact area. In the framework proposed in this paper, the contact solution is achieved using an optimization scheme [30] seeking the minimum of a quadratic form (i.e., the complementary energy, subjected to constraints, i.e., the boundary conditions). The convergence of this quadratic optimization is guaranteed, but the method fails when adhesion-like tensile contact tractions are assumed. It should be noted that adhesion was not considered either in the classic literature of viscoelastic contact, and its integration with SAM is still in its early stages [31].


## 3. The Contact Model

The contact parameters are assessed in an iterative manner, which requires both spatial and temporal discretization. The spatial discretization is needed as the contact area and the pressure distribution are a priori unknown and, moreover, keep changing during the contact process together with the mechanical response of the viscoelastic material. Throughout algorithm execution, integration of pressure over repeatedly changing contact area is executed in each new iteration. As the goal of this paper is to advance an algorithm for deterministic rough surfaces, pressure distribution and shape of contact area can be arbitrary during iterations, as a result of the arbitrary initial contact geometry. Therefore, displacement computation can only be achieved assisted by the spatial discretization.The spatial discretization principles are common to all SAM models and were discussed in detail in [16]. In a Cartesian coordinate system with the x→1-axis and x→2-axis laying in the common plane of contact, a computational domain P is chosen as a rectangular and uniformly spaced grid with N1×N2 elementary cells, whose side lengths are denoted by Δ1 and Δ2, with Δ=Δ1Δ2. A set of integers i=1,…,N1, j=1,…,N2 is used for indexing the grid cells. Assuming a controlling point for each elementary cell (usually the center of the cell), its coordinates can be written as x1(i)=(i-1/2)Δ1, x2(j)=(j-1/2)Δ2. All problem parameters are assumed to be piecewise constant, and the sampling is based on the discrete values computed in the control points, thus generating a discrete counterpart for each continuous distribution.The temporal discretization is required in the computation of surface viscoelastic displacement, which, according to [19], relies on the pressure history and on the time-dependent material property. The aim of this paper is to derive a versatile model, which can account for arbitrary loading history, as opposed to analytical solutions that are plagued by limitations regarding the monotony of the contact radius. To achieve this level of generality, the loading history needs to be discretized in small time intervals, sufficiently short so that contact pressure and viscoelastic mechanical properties may be assumed to be constant during each time step. A uniform temporal mesh of step Δt, with Nt time steps, is thus imposed, adding a third argument to problem parameters; for example, p(i,j,k) is the discrete counterpart of the continuous pressure p(x1(i),x2(j),t(k)), with t(k)=kΔt, and denotes the elementary pressure in the cell (i,j) of the spatial mesh, achieved during the contact process after k time steps in the simulation window. It is required that, at t=0 (i.e., the beginning of the simulation window), the viscoelastic body was undisturbed. This elementary pressure is uniform in the cell (i,j) and constant during thekth time interval. Model parameters having two arguments depend only on the spatial localization, whereas those with one argument depend only on time.The viscoelastic contact model is based on the general contact model for the elastic frictionless contact problem under normal loading developed by Johnson [32], completed with the viscoelastic displacement equation and the consideration of plasticity effect: (2)Wk=Δ∑i,j∈Akpi,j,k,k=1,…,Nt.(3)hi,j,k=hii,j+ui,j,k-ωk,i,j∈P,k=1,…,Nt.(4)pi,j,k>0,hi,j,k=0,i,j∈P,k=1,…,Nt;(5)pi,j,k=0,hi,j,k>0,i,j∈P-A,k=1,…,Nt.(6)pi,j,k<H,i,j∈P,k=1,…,Nt;(7)ui,j,k=∑n=2Nt∑l=1N1∑m=1N2Ki-l,j-m,k-npl,m,n-pl,m,n-1,i,j∈P,k=2,…,Nt.Here,W denotes the applied normal force,A the digitized contact area (i.e., a set of elementary cells best fitting the shape of the contact area), hi the initial (i.e., at time t=0) composite contact geometry (or the gap between the undeformed surfaces), h the gap between the deformed surfaces, u the relative displacement in the direction normal to the common plane of contact, ω the rigid-body approach,H the hardness of the softer material, and K(i-l,j-m,k-n) the viscoelastic influence coefficient, expressing the displacement induced in the (i,j) patch of the spatial mesh at time kΔt, by a uniform pressure of (Δ1Δ2)-1 Pa that acted in the patch (l,m) at the time nΔt. Its computation based on the purely elastic solution (i.e., the Boussinesq problem of a half-space subjected to a point load) and on the Boltzmann hereditary integral (1) is detailed in [19].Equation (2) imposes the static force equilibrium in the frictionless contact problem. Equation (3) is the geometrical condition of deformation, describing the separation between the contacting bodies. Relations (4) and (5) are the complementarity conditions, assessing the contact status of each cell in the computational domain: cells in contact have positive pressure and vanishing gap, while cells with positive gap are excluded from the contact area and consequently have vanishing pressure. It is apparent here that negative contact tractions are not allowed in this model; therefore, adhesion is not accounted for. Equation (6) limits pressure to the hardness of the material according to the Tabor equation [27], and (7) computes the viscoelastic displacement as a discrete cyclic convolution between the influence coefficients matrix and the history of pressure, as shown in [19].The contact problem can be regarded as an optimization problem having (3) as the merit function, subjected to restrictions consisting in the static force equilibrium (see (2)) and in the boundary conditions (see (4)–(6)).The difficulty in achieving the solution of model (2)–(7) stems from two facts. Firstly, the contact area is a priori unknown and keeps changing with time in the course of the contact process. Secondly, the contact simulation relies on contact evolution in the past period of time. Indeed, (7) shows that, in order to compute the viscoelastic displacement at time t=kΔt, all pressure distributions in the previous (k-1) steps are required.To overcome each of these obstacles, an individual level of iteration is implemented in the numerical treatment of this problem, resulting in a two-level nested loop. The inner loop stabilizes a solution for the instantaneous (i.e., withk fixed, but otherwise arbitrary) contact state consisting in (2)–(6), based on the assumption that the history of the pressure distribution is known, but otherwise arbitrary. The simulation of the viscoelastic contact process is achieved in the outer loop, by solving a series of sequential instantaneous contact states, related to each step from the temporal discretization.


## 4. Algorithm Description

The goal of the inner loop is to find the current pressure distributionp(i,j,k),(i,j)∈P, with k fixed but otherwise arbitrary, under the assumption that all historical pressures pi,j,1,…,p(i,j,k-1), (i,j)∈P, are known. The contact status of each cell in the computational domain is a priori unknown and, moreover, can change with each new iteration in the inner loop, as well as with the reproduction of the loading history in the outer loop. Furthermore, the number of cells in the contact area dictates the size of the linear system emerging from (3). At any iteration of the inner loop, cells with negative pressure are excluded from the contact area, as adhesion is not allowed, while cells with negative gap are reincluded, as bodies are considered impenetrable in the frame of Linear Theory of Elasticity. The associated equations are added to or removed from the linear system having the nodal pressures as unknowns. Essentially, both pressure and contact area are iterated simultaneously, and convergence is attained when pressure is stabilized in relation to the imposed normal load. This iterative approach is classical [16–18] in the contact mechanics of elastic bodies, but its application to viscoelastic materials is new. The main algorithm steps are detailed below.(1)
Initialize algorithm auxiliary variables:(8)θ=0;Rold=1;di,j=0,poldi,j,k=0,i,j∈P.(2)
Initialize all nodal pressures with the mean pressure on the computational domain (i.e.,A=P) in the first iteration, with A being the contact area, that is, the set of cells A=(i,j)∈P:p(i,j)>0. It should be noted that a vanishing pressure, leading to a vanishing contact area, does not make a suitable initial guess.(3)
Assess the displacement fieldu(i,j,k) from (7) by introducing the pressure history pi,j,1,…,p(i,j,k-1), (i,j)∈P, as well as the current pressure p(i,j,k), (i,j)∈P, adopted in the previous step. From this point forward, the third index reserved for temporal localization will be omitted for brevity in the description of the inner loop.(4)
Assess the rigid-body approach from the contact complementarity conditions (3) and (4):(9)ω=∑i,j∈A1-1∑i,j∈Aui,j+hii,j.In this manner, (3) emerges as a linear system in nodal pressures, which can be solved efficiently using the conjugate gradient method (CGM).(5)
Start the CGM loop by computing the residual and subsequently its square Euclidean norm for the set of cellsD=(i,j)∈P:p(i,j)>0,p(i,j)<H: (10)ri,j=hii,j+ui,j-ω,i,j∈P;R=∑i,j∈Dri,j2.(6)
Compute the descent directiond in the CGM: (11)di,j⟵ri,j-di,jθRRold,i,j∈D;0,i,j∈P-D,where a left arrow denotes the recalculation of a variable based on its old value.(7)
Memorize the norm of the residual for subsequent computation of descent direction:Rold=R.(8)
Assess the length of the stepα to be made in the descent direction computed in step (6): (12)α=∑i,j∈Dci,jdi,j-1∑i,j∈Dri,jdi,j,wherec(i,j) denotes the convolution between the influence coefficients matrix and the descent direction, adjusted by its mean value computed on the set D: (13)ci,j=∑l=1N1∑m=1N2Ki-l,j-mdl,m,i,j∈P;ci,j⟵ci,j-∑i,j∈A1-1∑i,j∈Dci,j,i,j∈P.(9)
Memorize the current contact tractions in a new variable for subsequent error computation:pold(i,j)=p(i,j), (i,j)∈P.(10)
Update the system solution by making a step of lengthα along the descent direction d: (14)pi,j⟵pi,j-αdi,j,i,j∈D.It should be noted that the CGM loop consisting in the algorithm steps(5) to (10) only constrains the nodal pressures to comply with (3) in the sense of residual minimization. The remaining restrictions in the contact model are enforced in the next steps.(11)
Impose the complementarity conditions in (4) and (5). As the contact model only allows positive interfacial tractions, cells with negative pressure are removed from the contact area, and the corresponding nodal pressures are set to zero in the set M=(i,j)∈P:p(i,j)<0: (15)A⟵A-M,pi,j=0,i,j∈M.(12)
Enforce the cut-off value for contact pressure to account for plastic effects in the set of cellsY=(i,j)∈P:p(i,j)>H: (16)pi,j=H,i,j∈Y.(13)
On the other hand, as bodies are considered impenetrable, cells having vanishing pressure but negative gap reenter the contact area. The nodal pressures in the setI=(i,j)∈P:p(i,j)=0,h(i,j)<0 are than recalculated by adding the positive term -αh(i,j), (i,j)∈I, as discussed in [16]: (17)A⟵A∪I,pi,j⟵pi,j-αhi,j,i,j∈I.Steps(11) and (13) lead to a reassessment of the contact or noncontact status of every cell in the computational domain. With each status change, the number of equations in the linear system resulting from (3) also changes. With the introduction of new cells in the contact area, the corresponding nodal pressures have no precedent in the residual minimization process of the CGM; therefore, a new search for the optimal path (i.e., the descent direction) must be conducted. This is accomplished in the next step.(14)
Update the auxiliary variableθ by setting it to zero if any status changes, or to unity otherwise. In the former case, it will lead to a reset of the conjugate descent directions in step (6) of the next iteration.(15)
Adjust the nodal pressures to comply with the static force equilibrium (2): (18)pi,j⟵pi,jWΔ∑i,j∈Api,j-1.(16)
Verify solution convergence to the imposed precisionε:(19)W-1Δ∑i,j∈Api,j-poldi,j≤ε.If pressure distribution is not stabilized, a new iteration is executed starting from step(3). In this manner, the inner loop derives the solution of any instantaneous contact state based on the knowledge of all previous states. Capitalizing on this result, the simulation of the viscoelastic contact process can be achieved in an outer loop, whose steps are detailed below.(1)
Acquire the input of the viscoelastic contact problem: the loading historyW(t), t∈[0,t′] (where t′ denotes the upper limit of the simulation window), the initial contact geometry hi(x1,x2), and the contact compliance (i.e., the creep compliance function Φ(t) for the viscoelastic material).(2)
Adopt the initial guess for the computational domainP and the discretization parameters: N1, Δ1, N2, Δ2, Nt, and Δt.(3)
Digitize all problems inputs. ObtainW(k), k=1,…,Nt, hi(i,j), (i,j)∈P, Φ(k), k=1,…,Nt.(4)
Compute the array of influence coefficientsKi,j,k, i,j∈P, k=1,…,Nt, as discussed in [19].(5)
Solve an initial contact state (k=1), as a purely elastic contact problem, but with the contact compliance Φ(1), and obtain the initial pressure distribution p(i,j,1), (i,j)∈P.(6)
Increasek (i.e., apply a new time increment) and solve a new instantaneous contact state to get p(i,j,k), (i,j)∈P.(7)
If the end of the simulation window is reached (i.e.,k=Nt), stop program execution and export data; otherwise, go to step (6).The computational complexity of a convolution product between pressure and the influence coefficients matrix for a specified time step is of the orderO(N12N22) with direct multisummation but can be conveniently reduced to O(N1N2log⁡N1N2) by the DCFFT technique [20]. When summation over time dimension is applied to simulate the memory effect, the computational complexity becomes of the order O(Nt2N1N2log⁡N1N2). Practically, a viscoelastic contact scenario with 1024×1024 spatial discretization and 100 time steps can be solved on a 3 GHz quad-core processor in less than an hour. Due to its achievable resolution, the algorithm is expected to tackle rough contact simulations involving real microtopography data obtained from 3D surface imaging devices.


## 5. Numerical Simulations and Discussions

Program validation is achieved by comparison with existing analytical solutions for the single asperity smooth contact of linear viscoelastic materials. A nominally flat viscoelastic half-space described by the Standard Linear Solid Model is indented by a rigid spherical punch of radiusR. The simulated loading program is a step loading W(t)=W0H(t), where H(t) denotes the Heaviside step function and W0 is fixed but otherwise arbitrarily chosen. The Standard Linear Solid Model, also known as the Zener model, consists in a spring of shear modulus G in series with a Kelvin model, or in a spring in parallel with a Maxwell model, and is described by the following creep compliance function: (20)Φt=121G+1GK1-exp-tτ,where GK and ηK are the shear modulus and the dashpot viscosity of the Kelvin component, respectively, and τ=ηK/GK is the relaxation time, namely, the time that stress takes in the Kelvin unit to decay by a factor of e (i.e., Euler’s number). All model parameters are assumed to be fixed but otherwise can be arbitrarily chosen. Ting’s solution [4, 5] for this contact scenario yields the following equation for the pressure distribution achieved at time t after the first contact, at the radial coordinate r: (21)pt,r=8GπRat2-r21/2-1τ∫0texp⁡2t′-tτRe⁡at′2-r21/2dt′,where the contact radius at the time t is given by a(t)=3RW0Φ(t)/81/3 and Re() denotes the real part of its complex argument. It is recalled that (21) only holds for monotonically increasing contact radius. The pressure history predicted by computer simulation on the newly proposed model, depicted with continuous black lines in Figure 1, matches well the data resulting from numerical computation of relation (21), displayed using dashed grey lines. Dimensionless pressure distribution and radial coordinate are defined as ratio to Hertz central pressure pH and contact radius aH, computed for the loading level W0 by setting t=0 in (20). No plasticity effect is considered in this simulation. The good validation proves that the simulation is able to capture the specifics of the viscoelastic contact process in a wide temporal window.Figure 1
Comparison of pressure history in the step loading spherical indentation of a Zener half-space.The effect of the cut-off limit on contact pressure is then simulated at the single asperity level, by neglecting the viscosity of the material. Figure2 depicts the predicted radial pressure distributions for various material hardness values, defined as ratios to Hertz central pressure pH. The limiting of pressure leads to an increase in the contact radius, similar to the one predicted by other numerical models [17] for the elastic-plastic contact. It should be noted that, in [17], the computed plastic strain generates residual displacement which increases contact conformity, resulting in a more uniform pressure distribution on an enlarged contact area, whereas, in the present model, similar pressure and contact area evolution are achieved by simply imposing additional constraints in the contact solver. The latter approach, while not being able to predict the residual state (i.e., the residual stresses and strains), has the advantage of preserving the method computational complexity. This allows concentration of computational resources on the simulation of the memory effect due to viscosity rather than on particular features of pressure at the asperity level due to plasticity. Considering that the number of samples available for each individual asperity is inevitably limited, we consider that a reasonable compromise is thus reached with this approach.Figure 2
Pressure profiles in the elastic spherical contact with material hardness as cut-off limit.The Maxwell and Kelvin rheological models for viscoelasticity might be adequate for qualitative or conceptual analyses, but quantitative representation of the behavior of real materials requires an increase in the number of parameters. The naturally occurring spectrum of relaxation times of a real viscoelastic material can be described by a Wiechert model, consisting in several Maxwell units and a free spring of stiffnessG0 connected in parallel. Its relaxation modulus function can be expressed with the aid of the spring stiffness Gi and the dashpot viscosity of each Maxwell component i, as a Prony series: (22)Ψt=G0+∑iGiexp⁡-tηiGiHt.Many indentation tests aiming to describe the time-dependent modulus of typical viscoelastic materials employ polymethyl methacrylate (PMMA), a thermoplastic polymer for which a linear viscoelastic behavior can be assumed. The viscoelastic response of PMMA up to 1000 sec under uniaxial compression was measured in standard relaxation tests by Ramesh Kumar and Narasimhan [33], and the fit of their experimental data using a two-term Prony series yields the following relaxation modulus: (23)Ψt=1973+254.776·exp⁡-t8.93+263.6628·-t117.96Ht,MPa.Unlike the case of purely elastic material, the viscoelastic creep and relaxation functions are not reciprocal in the time domain (i.e.,Ψ(t)Φ(t)≠1), but a corresponding relation can be established in the Laplace transform domain (i.e., Ψ-(s)Φ-(s)=1/s2), with s being the variable of the Laplace transform. The latter equation can be used to derive the creep compliance of the PMMA in the time domain by inverse Laplace transform. The resulting creep compliance curve is plotted in Figure 3.Figure 3
Creep compliance of the PMMA.Further validation of the numerical model can be achieved by simulating the load-displacement curve in the indentation of a PMMA half-space by a rigid spherical indenter. The load history was imposed according to (24), consisting in a normal force ramped to a peak value W0 in a time interval [0,t0] and subsequently ramped down in the time range [t0,2t0] (i.e., a triangular loading profile): (24)Wt=W0tt0Ht-2tt0-1Ht-t0.The variation of the indentation depth with the imposed load, traced in Figure4 for W0=100N and t0=10s, matches well the curve obtained experimentally in [33].Figure 4
Load versus indentation depth for triangular loading of the PMMA material.The capacity of the numerical model to simulate the viscoelastic rough contact is demonstrated by performing a contact simulation involving a deterministic roughness sample [34]. The discrete heights of 1024×1024 equidistant points laying in a rectangular patch of side lengths 0.1×0.1mm are superimposed to a smooth sphere of radius 0.5mm, having its center at coordinates x1=0, x2=0. The resulting composite initial geometry, that is, the sum of corresponding heights of the two contacting bodies, is depicted in Figure 5. The consideration of the spherical counter body helps relieve the spurious pressure concentrations at the edges of the computational domain. It should be noted that the points with the lowest height come into contact first.Figure 5
Composite initial geometry for the rough contact simulation.The contact is driven by a step loading withW0=0.5 N in a simulation covering a time window of 500 s divided into 100 time steps. Due to steep slopes in surface microtopography, the purely elastic contact solver predicts for the moment t=0 a contact area consisting in a few elementary patches only, with nodal pressures two orders of magnitude greater than the hardness of the PMMA, of 300 MPa [33]. It is clear that such levels of stress cannot be sustained by the viscoelastic material and localized plastic deformation will occur near the tip of the asperities that are first brought into contact. The limiting of pressure helps alleviate these spurious nodal pressures and creates a more realistic description of the initial contact area, as depicted in Figure 6. The black spots represent elementary cells in contact.Figure 6
Contact area map at timet=0 s.The contact simulation predicts that the real contact area increases with time. Initial contact occurs near sharp asperities and leads to plastic deformation, whereas as time goes on, additional lower asperities are brought into contact while the initial ones become flattened due to creep in the viscoelastic material. These assertions are supported by comparison of the contact area maps presented in Figures6 and 7, showing additional contact spots in the end of the simulation window at t=500 s. The evolution of the contact area and that of the indentation depth with time are depicted in Figure 8.Figure 7
Contact area map at timet=500 s.Figure 8
Evolution of contact parameters.


## 6. Conclusions

A versatile algorithm for the simulation of rough contact processes involving linear viscoelastic materials is advanced in this paper, by coupling a semianalytical method for the computation of viscoelastic displacement with a contact solver optimized for rough contact simulations involving plastic deformation at the asperity tips.The contact model features both spatial and temporal discretization. The spatial discretization ensures that the instantaneous contact state can be solved for arbitrary contact geometry, in an iterative scheme that assesses simultaneously the contact area and the pressure distribution by performing a quadratic optimization whose convergence is guaranteed. The temporal discretization handles the reproduction of the memory effect characteristic to viscoelastic materials and allows simulation of arbitrary loading programs with complex constitutive models of linear viscoelasticity in a wide temporal window, capturing both the elastic and the creep response of the material. The solution of the instantaneous contact state is achieved at every new time step, based on information from all previous states.The numerical program was validated by comparison with existing analytical and experimental results for smooth contact geometry, involving the Zener rheological model or polymethyl methacrylate, in step or ramped loading. Rough contact simulations suggest that the real contact area increases with time as additional lower asperities come into contact, while the initial ones become flattened due to creep in the viscoelastic material.The novel algorithm can tackle the linear viscoelastic rough contact with roughness samples of more than 106 points without any convergence problems, and it has a remarkably high computational efficiency compared to other numerical methods capable of handling this type of problems.Although adhesion is neglected and plastic yielding is considered only in a simplified manner, lacking the assessment of the residual state, the model is expected to predict the contact behavior of deterministic rough surfaces with a degree of accuracy relevant to practical applications.



---

*Source: 2521903-2017-02-15.xml*