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.
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 stress and of deviatoric strain by means of the shear modulus (i.e., ). 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 and the creep compliance :
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 , chosen so that initially (i.e., at ) 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 -axis and -axis laying in the common plane of contact, a computational domain is chosen as a rectangular and uniformly spaced grid with elementary cells, whose side lengths are denoted by and , with . A set of integers , 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 , . 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 , with time steps, is thus imposed, adding a third argument to problem parameters; for example, is the discrete counterpart of the continuous pressure , with , and denotes the elementary pressure in the cell of the spatial mesh, achieved during the contact process after time steps in the simulation window. It is required that, at (i.e., the beginning of the simulation window), the viscoelastic body was undisturbed. This elementary pressure is uniform in the cell and constant during the kth 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:
Here, 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), the initial (i.e., at time ) composite contact geometry (or the gap between the undeformed surfaces), the gap between the deformed surfaces, 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 the viscoelastic influence coefficient, expressing the displacement induced in the patch of the spatial mesh at time , by a uniform pressure of Pa that acted in the patch at the time . 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 , all pressure distributions in the previous 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., with k 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 distribution , with fixed but otherwise arbitrary, under the assumption that all historical pressures , , 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: (2)Initialize all nodal pressures with the mean pressure on the computational domain (i.e., ) in the first iteration, with being the contact area, that is, the set of cells . 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 field from (7) by introducing the pressure history , , as well as the current pressure , , 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): 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 cells : (6)Compute the descent direction in the CGM: 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: .(8)Assess the length of the step to be made in the descent direction computed in step : where denotes the convolution between the influence coefficients matrix and the descent direction, adjusted by its mean value computed on the set : (9)Memorize the current contact tractions in a new variable for subsequent error computation: , .(10)Update the system solution by making a step of length along the descent direction : It should be noted that the CGM loop consisting in the algorithm steps to 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 : (12)Enforce the cut-off value for contact pressure to account for plastic effects in the set of cells : (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 set are than recalculated by adding the positive term , , as discussed in [16]: Steps and 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 of the next iteration.(15)Adjust the nodal pressures to comply with the static force equilibrium (2): (16)Verify solution convergence to the imposed precision :
If pressure distribution is not stabilized, a new iteration is executed starting from step . 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 history , (where denotes the upper limit of the simulation window), the initial contact geometry , and the contact compliance (i.e., the creep compliance function for the viscoelastic material).(2)Adopt the initial guess for the computational domain P and the discretization parameters: , , , , , and .(3)Digitize all problems inputs. Obtain , , , , , .(4)Compute the array of influence coefficients , , , as discussed in [19].(5)Solve an initial contact state (), as a purely elastic contact problem, but with the contact compliance , and obtain the initial pressure distribution , .(6)Increase (i.e., apply a new time increment) and solve a new instantaneous contact state to get , .(7)If the end of the simulation window is reached (i.e., ), stop program execution and export data; otherwise, go to step .
The computational complexity of a convolution product between pressure and the influence coefficients matrix for a specified time step is of the order with direct multisummation but can be conveniently reduced to by the DCFFT technique [20]. When summation over time dimension is applied to simulate the memory effect, the computational complexity becomes of the order . Practically, a viscoelastic contact scenario with 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 radius . The simulated loading program is a step loading , where denotes the Heaviside step function and is fixed but otherwise arbitrarily chosen. The Standard Linear Solid Model, also known as the Zener model, consists in a spring of shear modulus 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: where and are the shear modulus and the dashpot viscosity of the Kelvin component, respectively, and is the relaxation time, namely, the time that stress takes in the Kelvin unit to decay by a factor of (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 after the first contact, at the radial coordinate : where the contact radius at the time is given by and 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 and contact radius , computed for the loading level by setting 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.
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. Figure 2 depicts the predicted radial pressure distributions for various material hardness values, defined as ratios to Hertz central pressure . 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.
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 stiffness connected in parallel. Its relaxation modulus function can be expressed with the aid of the spring stiffness and the dashpot viscosity of each Maxwell component , as a Prony series:
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:
Unlike the case of purely elastic material, the viscoelastic creep and relaxation functions are not reciprocal in the time domain (i.e., ), but a corresponding relation can be established in the Laplace transform domain (i.e., ), with 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.
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 in a time interval and subsequently ramped down in the time range (i.e., a triangular loading profile):
The variation of the indentation depth with the imposed load, traced in Figure 4 for and , matches well the curve obtained experimentally in [33].
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 equidistant points laying in a rectangular patch of side lengths mm are superimposed to a smooth sphere of radius mm, having its center at coordinates , . 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.
The contact is driven by a step loading with 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 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.
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 Figures 6 and 7, showing additional contact spots in the end of the simulation window at s. The evolution of the contact area and that of the indentation depth with time are depicted in Figure 8.
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.
Competing Interests
The authors declare that they have no competing interests regarding the publication of this paper.
Acknowledgments
This work was partially supported by the project Integrated Center for Research, Development and Innovation in Advanced Materials, Nanotechnologies, and Distributed Systems for Fabrication and Control, Contract no. 671/09.04.2015, Sectoral Operational Program for Increase of the Economic Competitiveness cofunded by the European Regional Development Fund.