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.
In the time domain, the creep compliance and relaxation modulus are not reciprocal (like in the purely elastic case); that is, . However, in the Laplace transform domain, the following relation applies [18] to their transforms:where is the variable in the Laplace transform domain. The latter equation can be used to derive the creep compliance function of PMMA by computing first its Laplace transform from (9) and by subsequently applying inverse Laplace transform to obtain in the time domain, leading to
The latter equations express the PMMA creep compliance in terms of shear modulus by using a Poisson’s ratio [21]. The semi-analytical solution of the contact problem involving linear viscoelastic materials described by complex rheological models like the Prony series is discussed in the following sections.
3. Contact Model
The contact model employed in this paper is based on the general contact model developed by Johnson [20] for the elastic domain. The contact equations, as well as the imposed assumptions and limitations, are repeated here for clarity, and the newly established dependencies, related to the viscoelastic constitutive law described in the previous section, are then discussed in detail.
The set of equations and inequalities governing the contact process are written in a Cartesian coordinate system with and axes laying in the common plane of contact (i.e., the plane passing through the initial point of contact, which separates best the two contacting surfaces). The two contacting solids are subjected to a normal force aligned with axis, compressing the two bounding surfaces together. As opposed to a time-independent purely elastic contact process, in which the final state depends only on loading level, the viscoelastic contact state depends on time, as well as on the loading history, due to the memory effect of the viscoelastic materials, thus adding a third time parameter to the elastic contact model.
The static force equilibrium relates the normal force to the pressure distribution at any time in the observation window . To keep the number of independent parameters to a minimum, the contact is assumed to be frictionless, meaning that shear tractions cannot be sustained at the contact interface. Moreover, the problem is considered as quasi-static, meaning that the inertia forces due to deformation are negligible:
The equations of the surface of separation between the two contacting bodies yield the geometrical conditions of deformation in the normal direction:where is the gap between the undeformed (i.e., initial, at time ) surfaces, is the gap between the deformed surfaces, is the relative normal displacement, and is the rigid-body approach.
The contact model is completed with the boundary conditions and constraints, also referred to as the complementarity conditions in the literature of the elastic contact. The latter equations are required as the contact area is not known a priori and consequently must be found in an iterative manner by a trial-and-error approach. The gap between the deformed contacting surfaces vanishes on the contact area, as no interpenetration of the contacting solids is allowed in the frame of elasticity. On the other hand, the gap must be positive outside the contact area, where there is clearance between the contacting bodies. In the same manner, pressure is positive on the contact area and vanishes outside the contact area. These boundary conditions and constraints that must be satisfied simultaneously can be expressed as
The assumption of non-negativity of pressure leads to neglect of contact adhesion and can be considered very conservative in the case of viscoelastic materials. Adhesion appears virtually in all contacts between real surfaces, but the force of adhesion can be often neglected in case of metallic materials, when the actual 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 that requires the non-negativity of contact tractions. The latter are obtained as the solution of a variational problem originally formulated in the field of contact mechanics by these authors [22], seeking the minimum of a quadratic form, that is, the complementary energy, subjected to constraints, that is, 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 in either the classic or modern literature [3–15] of viscoelastic contact. The contact complementarity conditions imposed in the elastic contact model used in this paper match the boundary conditions employed in [7] and lead to a surface displacement compatible with the indenter profile within the contact area. As the integration of adhesion in analytical or semi-analytical contact models is still in its early stages [23], no step back is made in the framework proposed in this paper.
The difficulty in solving the contact model (11)–(14) stems from the fact that neither the contact area nor the pressure distribution is known in advance. An iterative approach is therefore needed, involving a trial-and-error approach, in which a contact region is assumed, and the pressure distribution is then computed based on this assumption. If all constraints in the contact model are verified by the obtained solution, the contact problem solution is achieved. This solution is unique based on the theorem of uniqueness of solution of the elastostatic problem. Otherwise, the contact area is adjusted and a new pressure distribution is computed with the new guess. This iterative approach requires that the response of the contacting material, that is, the displacement induced by the surface tractions, is computed for arbitrary contact area and pressure distribution. The latter computation can only be achieved semi-analytically, and therefore a problem discretization is imposed to perform the numerical analysis of the contact process.
The contact model reviewed herein was also used extensively in the simulation of contact scenarios involving history-dependent processes like plasticity [24], wear [25], or friction [26] by adding an external loop in which the load was applied incrementally. In the latter contact scenarios, the time parameter does not need to be considered explicitly as long as the history of the contact process is properly simulated (i.e., the load is applied with sufficiently small increments). The present work attempts to link the contact model to the theory of viscoelastic behavior, in which the material properties depend explicitly on time. Manipulation of the existing semi-analytical solution in the elastic domain, aiming to achieve a viscoelastic contact algorithm, is detailed in the following sections.
4. Viscoelastic Displacement
A surface distribution of normal tractions, such as the pressure resulting from a mechanical contact process, induces a displacement field whose knowledge is essential in solving the contact problem and in performing stress analysis in the contacting bodies. Although the limiting boundary of a real solid is intrinsically rough, computational contact mechanics employ the half-space assumption, allowing for the use of fundamental solutions (i.e., the Green functions) derived in the theory of linear elasticity for a semi-infinite body bounded by a plane surface. For this approximation to remain valid, the slope of the contact geometry must remain small throughout the contact region. The normal displacement field generated in a linear elastic and isotropic solid by a distribution of normal tractions is computed by applying the superposition principle to the Green function for the elastic half-space derived by Boussinesq [27]:where is the normal displacement induced at a point of coordinates by a unit concentrated force acting in origin along direction of , and and are Poisson’s ratio and the shear modulus of the elastic half-space.
Lee and Radok [3] obtained the contact radius in the viscoelastic spherical contact problem by applying the hereditary integral operator of type to the Hertz (i.e., purely elastic) solution in which the elastic compliance was replaced by the viscoelastic creep compliance . This course of action is justified by the classic method for solving the linear viscoelastic problems of stress analysis, which is based on the concept of associated elastic problem [1, 2]. Capitalizing on the fact that basic integral equations for stress analysis in viscoelastic materials reduce in the Laplace transform domain to the type of integral equations describing stresses in elastic materials, it has been shown [1, 2] that a viscoelastic problem has an associated elastic problem, to which the former reduces after removal of time dependency by application of the Laplace transform. Consequently, if the boundary conditions are time-independent, a solution in the frequency domain is identical in form to the corresponding elastic solution. This technique of deriving viscoelastic solutions from their elastic counterpart is also referred to as the correspondence principle. The indentation of a viscoelastic half-space by a rigid indenter cannot generally be solved in this manner, as the contact problem features time-dependent boundary conditions, which impede transfer to Laplace domain. When applying this technique to the contact radius formula in the associated Hertz elastic problem, the viscoelastic solution holds true as long as the contact radius increases monotonically with time, but additional manipulations are required when the time-dependent contact area is an arbitrary function of time, as shown in [6, 7].
In this paper, the same technique of replacing the elastic contact compliance with the viscoelastic creep compliance function is applied to (15):yielding the viscoelastic displacement generated by a known history of pressure in a window of observation :
Unlike its counterpart expressing the contact radius in a Hertz-type viscoelastic contact, the displacement equation (17) does not require additional manipulations and can be used in conjunction with any history of boundary conditions. Interchanging differentiation and integration in (17) yields
In a viscoelastic contact problem, the contact area and the pressure distribution are not known in advance and, moreover, keep changing during the contact process, as the response of the viscoelastic material also changes with time. Consequently, integral (18) must be evaluated for various loading histories, implying integration of arbitrary functions over arbitrary domains. The semi-analytical treatment of these equations to attain a computationally friendly form is detailed in the following sections.
5. Problem Discretization
Neither integral (15) for the elastic case nor (18) for the viscoelastic framework can be computed analytically for general contact geometry, loading history, or material properties. Important research efforts were dedicated to obtaining the solution of these integrals in a semi-analytical manner [17, 28]. The principle of the semi-analytical method consists in considering all continuous distributions as piecewise constant functions, uniform within each surface element in a uniformly spaced rectangular mesh established in the common plane of contact. Control points must be chosen for all the elementary cells of the grid (the centers of the cells make good candidates), and all continuous parameters are evaluated in these representing points, resulting in a digitized counterpart for each continuous distribution. This discretization encourages a simplified notation taking as arguments the indexes of the cells rather than the continuous coordinates. For example, , with and integers, denotes the pressure value computed in the center of the cell , and , where and are coordinates of the center of the cell . Consequently, pressure is assumed to be uniform in any rectangular patch and therefore can be factored outside the integral operator in (15). The integral of the Green function , taken over the elementary patch of side lengths and along directions of and , respectively, yields the influence coefficient [28] for the elastic displacement :which expresses the normal displacement induced in the observation cell by a uniform pressure of magnitude Pa acting in the cell . The closed-form solution of the double integral (19) was derived by Love [29]:
Within this framework, the semi-analytical counterpart of (15) results aswhere and denote the numbers of grids along directions of and , respectively. The discrete double convolution in (21) can be performed for any imposed pressure distribution. Optimum algorithmic efficiency is achieved using the DCFFT algorithm advanced by Liu et al. [17]. The reduction of the order of computation stems from the convolution theorem, which states that the convolution operation reduces to an element-wise product in the Fourier transform domain. The semi-analytical displacement computation using (21) together with the DCFFT technique is now widely used in computational contact mechanics. In this paper, a generalization of this technique to the case of viscoelastic behavior is proposed.
As equations describing the purely elastic model are intrinsically time-independent, spatial discretization is adequate to circumvent the continuous integration in (15). The additional integration over the time span in which the body was loaded in (17) requires an additional temporal mesh capable of simulating the memory effect specific to viscoelastic materials (i.e., the property that the current state depends upon all previous states succeeded from the first loading). This temporal discretization should be chosen so that at the body was undisturbed, and the time increment should be small enough so that, during each step, the problem parameters can be assumed to be constant. A piecewise constant law is thus imposed along the temporal axis, adding a third parameter to the notation implemented in the purely elastic model. For example, is the discrete counterpart of , denoting the pressure in the elementary cell in the spatial mesh, achieved after time increments, where , with . These assumptions regarding the temporal variation of model parameters authorize the substitution of the partial derivative in (18) with the finite difference . Discretization of (18) in the time domain yields
It should be noted that, in the latter equation, the simplified notation related to problem digitization was used only for the temporal parameter, whereas continuous coordinates are employed for spatial localization. As pressure is uniform within each elementary patch, it can be factored outside the spatial integral operator as in the purely elastic model, allowing for a viscoelastic influence coefficient defined similarly to its elastic counterpart in (19):
The latter equation employs notation related to discretized parameters in both spatial and temporal dimensions. The influence coefficient (23) expresses the displacement observed after time steps in the elementary patch of the spatial mesh, due to a uniform pressure of Pa which acted in the patch in the th time step after the reference moment in which the body was undisturbed, with , , and . The semi-analytical counterpart of (18), discrete in both time and space dimensions, can thus be expressed aswhere , , and . This equation clearly shows that the memory effect is considered explicitly in the displacement computation, as pressure distributions in all previous states (i.e., in all previous time increments), together with the current pressure, are needed to evaluate the current displacement. It is noteworthy that the contribution of the historical pressures can be separated from the contribution of the current pressure, which leads to the algorithmic strategy described in the following section.
6. Algorithm Overview
The viscoelastic contact simulation is achieved by constructing a series of elastic contact problems with boundary conditions that match the ones of the viscoelastic contact problem at each new time increment in the temporal discretization. This approach is based on the fact that, provided the compatibility and internal equilibrium equations are satisfied instantaneously, any elastic solution to the contact problem befits an instantaneous viscoelastic solution. The set of equations and inequalities (11)–(15) describe in fact a purely elastic contact process, whereas substitution of (15) with relation (18) accomplishes the algorithm generalization to the viscoelastic constitutive law.
Due to the robust nature of the semi-analytical formulation for the purely elastic contact process, history-dependent processes can also be simulated using this algorithm by applying the load in small increments, assuring that the loading path is properly reproduced. In the case of viscoelastic materials, however, the time parameter appears explicitly, and the contact parameters are expected to vary even when the load level is kept constant.
Essentially, the nodal pressures result as the solution of the system of equations arising from digitization of (12). The latter is essentially non-linear due to the presence of the rigid-body approach but can be linearized by adopting an estimate for this parameter based on the current pressure. Once initial guess values for pressure and the contact area are adopted, (12) can be considered for every cell in the contact area, and the resulting set of equations, all having vanishing gap , provide an estimation of . In other words, the rigid-body approach, the pressure distribution, and the contact area are iterated simultaneously, and the true value of is found when all contact parameters converge. The advantage of this approach consists in linearization of (12), which can thus be solved using appropriate numerical methods for linear systems of equations. According to this review [30] on the subject of numerical methods applied to contact mechanics, the best candidate is the CGM, providing the fastest rate of residual decrease.
Equation (12) applied to all cells in the computational domain generates discrete equations, but only the ones related to cells in the contact area form the system to be solved. Additional difficulty arises as the contact area is not known in advance, and therefore both pressure distribution and contact area have to be iterated simultaneously. Consequently, the size of the linear system, which is directly linked to the number of cells in the contact area, may also vary during the iterative process. At any iteration, cells with negative pressure are excluded from the contact area, while cells with negative gap (i.e., cells where the contacting bodies are predicted to overlap) are (re)included. These adjustments are requested by the contact complementarity conditions in (13) and (14) and force the computed nodal pressures to comply to the boundary conditions of the contact problem. As depicted in Figure 2, these restrictions are imposed outside the core of the CGM-type iteration. The resizing of the linear system leads to reconsideration of the descent directions and the descent steps previously computed in the CGM residual minimization process.
In order to use the CGM, which is proven to converge only for systems having a symmetric and positive definite matrix, it is not convenient to include the static equilibrium equation (11) in the system. In this manner, the system matrix is in fact the influence coefficients matrix, which is symmetric and positive definite (as a diagonally dominant matrix). It should be noted that the diagonal entries in the influence coefficients matrix (for any time ) are reserved for the coefficients expressing the contribution of the pressure located in each cell to the displacement in the same cell. Moreover, (20) suggests that the influence coefficients decay rapidly with the distance between the excitation (i.e., pressure) and the effect (i.e., displacement); therefore the influence coefficients matrix is diagonally dominant for fine meshes. The static equilibrium equation (11) is imposed during each iteration of the CGM algorithm in an additional correction of the system solution outside the CGM core, as depicted in Figure 2. This correction consists in a modification of the nodal pressures proportional to the ratio of the current load to the imposed load.
The iterative process stops when such a system is found that its solution in pressure features only positive entries and the gap distribution takes only positive values on the entire computational domain. The finding of a set of equations of type (12) whose solution additionally satisfies both boundary constraints and static equilibrium marks the achievement of simultaneously converging solutions for pressure, contact area, and rigid-body approach. For each new time interval, a new set of such contact parameters has to be found, as the viscoelastic displacement is updated at each new time increment to account for the new pressure history.
At any iteration and for each newly considered time interval, the computational domain should fully enclose the contact area; otherwise, the contact process simulation should be restarted with a larger estimate. As discussed in [28], if the predicted contact area reaches the boundaries of , spurious pressure concentrators due to end effects compromise the solution precision. This is especially true for rough contact scenarios, in which the specimen roughness sample should be properly handled by extending the geometry information with smooth convex shapes. For a detailed description of the algorithm to solve the frictionless normal contact problem, the reader is referred to [15, 28].
The contact solver is adapted in this paper to allow more general boundary conditions, thus addressing an important part of indentation scenarios, in which the rigid-body approach (i.e., the displacement history), rather than the load (i.e., the surface stress history), is imposed. The contact solver advanced by Polonsky and Keer [28] is load-driven (LD) (i.e., the applied normal force is expected as input) but can be modified for displacement-driven (DD) scenarios (i.e., in which the load is not specified, but the rigid-body approach is imposed). A distinctive feature of the original solver [28] is that the normal approach is not explicitly computed. This can be achieved as the algorithm [28] minimizes the residual computed from a modified interference equation of type (12):in which the normal approach is not directly considered. As rigid-body approach is independent of spatial localization (but dependent of time), its contribution is intrinsically accounted for when the residual is corrected by its mean value [28]. In the DD formulation, on the other hand, the normal approach is known as input data, so the interference computation should be conducted as in the theoretical model. The system residual matches the gap between the deformed surfaces and vanishes on the contact area:
As a result, the residual correction by its mean value in the original algorithm [28] becomes redundant and should be removed. The lack of normal load as input, on the other hand, compromises the adoption of the guess value for the pressure distribution, as the mean pressure can no longer be computed. In our numerical simulations, no convergence problems occurred with nonvanishing uniform pressure distribution as initial guess. The algorithm sequence performing pressure correction with respect to the imposed load also becomes redundant and should be eliminated. The fact that the correct value of the rigid-body approach is known from the start coerces pressure to converge without any additional algorithm modifications. The speed of convergence for the modified DD version of the contact solver was found to be of the same order of magnitude as its LD counterpart. The flowchart for the generalized algorithm, comprising both common and specific steps pertinent to each type of boundary conditions, is depicted in Figure 2.
It was shown in the previous section that the instantaneous viscoelastic displacement can be derived based on the knowledge of all previous contact states. Capitalizing on this result, the simulation of the viscoelastic contact process can be achieved by solving a series of sequential contact states related to a temporal discretization fine enough so that the memory effect specific to viscoelastic materials is accurately captured. The main algorithm steps are detailed as follows:(1)Acquire the input of the viscoelastic contact problem: the loading history , , the initial contact geometry , and the contact compliance, that is, the creep compliance function of the viscoelastic material.(2)Adopt the computational domain and the spatial and temporal discretization parameters: , , , , , and .(3)Digitize all problems inputs. Obtain , , .(4)Compute the array of influence coefficients , , , which can be stored as a three-dimensional array, having two dimensions related to the spatial dimensions and the third dimension for the temporal discretization. In the th temporal step, with , only the first layers of this array are used.(5)Solve an initial contact state () with the contact compliance and obtain the related pressure distribution , . In this particular case, the surface displacement is simply the convolution between pressure and the influence coefficients, with no memory effect to account for:(6)Apply a new time increment (i.e., increase ) and solve the th instantaneous contact state to get . This can be achieved as all historical pressures are known at this point. The product between the influence coefficients and the historical pressure distributions can be included in (12) as an initial state (i.e., superimposed to ). In this manner, the structure of the model (11)–(14) remains unchanged, and the same type of algorithm can be applied to solve each new time step in the viscoelastic contact problem. This technique is similar to the one used in the semi-analytical resolution [24] of elastic-plastic contact problems, in which the residual displacement due to the plastic region is added to the contact geometry, forming a modified initial state.(7)If the end of the simulation window is reached (i.e., ), stop program execution and export the computed data. Otherwise, go to step 6.
Based on (24), it follows that the computational complexity of the proposed algorithm is of order when direct multi-summation is applied. The algorithmic efficiency can be increased dramatically by implementing the DCFFT technique [17] in computation of viscoelastic displacement, leading to an improved order of operations of . Practically, a viscoelastic contact scenario with a spatial discretization and 200 time steps can be solved on a 3 GHz quad-core processor in a matter of minutes.
In this paper, an additional conditional branch is added to allow for computation of displacement which persists (for a limited or an indefinite time period, according to the imposed rheological model) after contact opening. The semi-analytical equation (24) for viscoelastic displacement computation holds true even after contact opening by assuming vanishing pressure distributions for the time steps subsequent to load removal and until the load is resumed. However, an additional algorithm conditional branch should be added, as the contact solver cannot handle vanishing loads directly. The flowchart of the extended algorithm, proposed for the simulation of contact problems involving linear viscoelastic materials, with general boundary conditions and arbitrary loading programs, is presented in Figure 3. In that figure, spatial localization is omitted for brevity and integer superscripts are used to index the referred time increment.
7. Code Validation and Results
The newly proposed algorithm was benchmarked against the partially analytical results derived by Ting [6, 7] for the Maxwell and Zener viscoelastic half-spaces indented by a rigid spherical punch of radius in a step loading , where denotes the Heaviside step function and is fixed but otherwise arbitrarily chosen. It should be noted that whereas Ting’s framework requires the punch to be axisymmetric and leads to tedious algebraic manipulations due to the five possible algorithm branches when the contact area does not vary monotonically with time, the algorithm proposed in this paper can readily handle irregular contact geometry, generalized boundary conditions, and complex material properties.
The pressure distributions predicted by the newly advanced computer program at several times after the first contact, depicted with discrete symbols in Figures 4 and 5, match well the data resulting from numerical computation of relations (2) and (6), respectively, displayed using continuous lines. Dimensionless pressure distribution and radial coordinate are defined as ratio to Hertz central pressure and contact radius , both corresponding to the maximum loading level . The good agreement between the two data sets legitimizes the novel algorithm for viscoelastic contact simulation.
It should be noted that the contact model for the Maxwell material predicts a contact radius that grows indefinitely in indentation under constant load. The solution should, however, be rejected for large strain.
Three different loading programs are considered in the following simulations concerning the Maxwell half-space: case 1, instantaneous loading followed by ramp unloading ; case 2, instantaneous loading and subsequent polynomial unloading ; and case 3, sinusoidal loading . Dimensionless contact radii depicted in Figure 6 are normalized by the Hertz counterpart corresponding to the maximum load attained in the loading program (i.e., in all three cases). The time axis uses as normalizer for loading programs 1 and 2 and for the third case, with being the relaxation time of the Maxwell unit. The predicted contact radius trends agree well with Ting’s results [6], proving the ability of the computer code to handle arbitrary loading histories (i.e., contact radii whose evolution is described by increasing and/or decreasing functions of time). It is noteworthy that the maxima and minima of the contact radii do not necessarily coincide with those of the loading program.
In order to test the novel DD formulation for the linear viscoelastic domain, the simulation of a Zener viscoelastic half-space with , indented by a rigid sphere of radius , was conducted following an imposed history of rigid-body approach, as resulting from Ting’s model applied to a step loading, that is, described by the equation , with the contact radius given by (3). The predicted history of pressure distribution over a time range, similar to the one depicted in Figure 5, was compared to that generated by the numerical computation of (6). A good agreement between the two data sets was found, legitimizing the new displacement-driven contact model.
Figure 7(a) depicts the surface displacement profiles in the step loading of a contact between a Maxwell viscoelastic half-space and a rigid spherical indenter, followed by instantaneous unloading at time . The elastic component of the post-unloading displacement, that is, that corresponding to the elastic spring from the Maxwell model, is recovered instantaneously, while the component related to the dashpot unit is predicted to persist indefinitely. A more realistic scenario is presented in Figure 7(b), in which a Zener viscoelastic contact with is unloaded at time , subsequent to a step loading. In this case, the displacement is recovered gradually and is expected to vanish after a specific time interval. The post-unloading displacement profiles are displayed using red lines in Figure 7.
(a)
(b)
8. Dynamic Mechanical Analysis (DMA)
The inducement of a steady-state undulated response in a viscoelastic specimen (e.g., a uniaxial bar) during a steady-state oscillation test is often used [31] as a method for assessing the mechanical properties of the viscoelastic material. The DMA methodology employs the so-called complex modulus of the viscoelastic material, consisting in a real part, that is, the storage modulus related to the elastic energy storage, and an imaginary part, that is, the loss modulus related to the internal energy dissipation. This complex modulus results by applying the Fourier transform to the Boltzmann hereditary integral of type . The phase angle of the complex modulus is defined as the hysteresis angle of the strain response falling behind the stress excitation. In a DMA procedure, the response in amplitude and the phase angle of the viscoelastic specimen to various excitation frequencies are employed to compute the complex modulus in the frequency domain. The latter modulus can be further used [32, 33] in the characterization of the contact behavior of the viscoelastic material under harmonic load indentation.
An approximate analytical solution in terms of the viscoelastic material complex modulus for the indentation depth amplitude in a frictionless spherical indentation under sinusoidal load was derived by Huang et al. [34]. The latter solution is, however, based on the tentative viscoelastic contact solution [3] that assumes a monotonically increasing contact area during the loading history. The latter limitation can be considered too severe unless, in a loading program consisting in a sinusoidal load superposed on a carrier step load, the indentation depth amplitude due to the oscillation is negligible compared to that due to the step load. The aforementioned limiting assumptions can be released in the numerical simulation, and the analysis of the hysteresis energy loss under the cyclic load can be conducted with a precision challenged only by the spatial and temporal discretization errors.
The Zener viscoelastic contact performance under harmonic load is first analyzed in this section using the newly proposed semi-analytical model. The indentation loading history is sinusoidal, and the phase is chosen so that at the load vanishes together with its first derivative:where denotes the number of oscillations (cycles) that occur each time interval equal to . The latter parameter and the load amplitude are fixed but otherwise can be arbitrarily chosen. By varying (i.e., by varying the excitation frequency with respect to ), the effect of loading frequency on contact parameters can be assessed. Dimensionless time parameter is defined as ratio to .
The numerical simulation predicts that, apart from the first cycles, the harmonic excitation induces a steady-state sinusoidal response with the same frequency but out of phase. When is small, that is, in the range of 0.1 , the steady-state regime is attained almost immediately, as shown by the force-displacement curve depicted in Figure 8(a). The latter is similar to that of materials exhibiting elastic hysteresis. In this case, only the first loading (i.e., the ascending branch of the first cycle) is unique; subsequent loadings and unloadings overlap to an imposed precision. With increasing , a specific, but limited, number of cycles is needed for response stabilization, as shown in Figure 8(b). Figure 9 suggests contact parameters stabilization by comparison of responses in the 7th and 8th loading cycles. The presence of a hysteresis-type delay between the indentation depth and the excitation load is also demonstrated.
(a)
(b)
Further observations can be made based on analysis of contact response in the steady-state regime in an imposed frequency spectrum of the excitation load (with ranging between 0.1 and 1.5 Hz·s). The numerical simulations suggest that the amplitude of contact radius decays with increasing excitation frequency, while the amplitude of peak pressure increases, as shown in Figure 10. Another interesting feature is that, at the same loading level, the contact is more conformal on the descending branch of the sinusoid, as proven by the pressure profiles depicted in Figure 11. The gap is more apparent at smaller loading frequencies.
The force-displacement curves depicted in Figure 12 prove that the contact process is hysteretic. The out-of-phase response between load and indentation depth is the source of energy damping in the viscoelastic material. The energy loss occurring each cycle is quantifiable by the area enclosed by the hysteresis loop. As this loop is obtained in discrete form, fine temporal discretization is required to obtain its detailed description. The loops presented in Figure 12 were obtained by imposing 100 temporal steps per loading cycle, and the simulation time was extended until a steady-state regime was reached (a maximum of 10 cycles were necessary for the case Hz·s). The simulation results suggest that, for the investigated spectrum of excitation frequency and for the considered rheological model, the energy loss decays with the increasing of the excitation frequency.
A study of the variation of energy dissipation under various frequencies of the excitation load is conducted for the PMMA material, involving the spherical indentation of a PMMA half-space in a loading program consisting in a sinusoidal load with an amplitude of , fixed but otherwise arbitrarily chosen, superposed to a step carrier load of . The oscillation was applied after the contact stabilization under the step load. Considering the creep behavior of the contact condition obtained for the rough contact of a PMMA specimen in [15], 500 s were allowed for the contact to reach a steady state:
By conducting the numerical simulation under an extended frequency spectrum of the excitation load, the stabilized viscoelastic contact behavior and the energy loss per cycle were obtained without any simplifying assumptions. The predicted hysteresis loops are similar to those presented in Figure 12, but the energy loss per cycle, related to the internal viscous friction in the viscoelastic material, was found to possess a maximum at a specific harmonic load frequency, as depicted in Figure 13. The numerically computed areas of the hysteresis loops enclosed by the force-displacement curves, normalized by the maximum value, suggest that the energy loss rises to a peak value when increases from 0 to 0.0031 and then decays with further increase of the loading frequency to an asymptotic value. This behavior may be attributed to the use of two relaxation times in the PMMA relaxation curve, as opposed to the single relaxation time in the Zener rheological model. This result proves the ability of the proposed numerical model to simulate and predict the contact performance of viscoelastic materials with complicated mechanical properties, as resulting from actual experimental measurements.
9. Conclusions
A robust algorithm for the resolution of linear viscoelastic contact scenarios is advanced in this paper by generalizing an existing well-known numerical solution for the contact of linear elastic bodies. The method relies on applying the correspondence principle together with the Boltzmann hereditary integral to achieve numerical computation of viscoelastic displacement, based on a temporal discretization in which every new state is assessed considering all previous time increments, thus simulating the memory effect specific to viscoelastic materials. Compared to other existing methods, the new approach can readily handle arbitrary contact geometry, arbitrary creep compliance (which can be inputted in discrete form as resulting from experimental data), and arbitrary loading programs, leading to contact radii described by increasing or decreasing functions of time. Whereas the classical viscoelastic contact solution still requires numerical treatments that can often raise convergence issues, the newly proposed algorithm is based on the conjugate gradient method, whose convergence is guaranteed. The computational efficiency is increased using a well-established method for rapid calculation of convolution products.
The contact solver assessing the contact parameters from the geometrical condition of deformation is enhanced to allow for more general boundary conditions. The solution for the displacement-driven contact is achieved by modifying the algorithm for the load-driven contact scenario. Another model extension deals with computation of post-unloading displacement, allowing for contact openings during the loading program.
The predictions of the newly advanced numerical program are benchmarked against existing results for a linear viscoelastic half-space described by the Maxwell and Zener rheological models, indented by a spherical punch in step loading.
The contact simulations under harmonic cyclic load suggest that a steady-state, out-of-phase harmonic response is attained after a few cycles. For the Zener theoretical model, the amplitude of contact radius, as well as the energy loss per cycle, decays with increasing excitation frequency, while the amplitude of peak pressure increases. The contact is found to be more conformal on the descending branch of the loading profile. The numerical computation of the energy dissipation in polymethyl methacrylate under an extended frequency spectrum of the excitation load suggests that a maximum is attained at a specific frequency. Further increase of the harmonic load frequency leads to an asymptotic value of the area enclosed by the resulting hysteresis loop.
The performed numerical simulations prove the ability of the newly proposed algorithm to assist the dynamic mechanical analysis in derivation of the complex modulus of the viscoelastic material, providing an important theoretical tool for the study of the viscoelastic contact with fewer simplifying assumptions than other existing partially analytical solutions.
Data Availability
The data supporting the PMMA constitutive law are from previously reported studies, which have been cited. The results of the numerical simulations that support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was partially supported from 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) and Sectoral Operational Program for Increase of the Economic Competitiveness co-funded from the European Regional Development Fund.