Advances in Tribology

Volume 2018, Article ID 6498503, 13 pages

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

## Numerical Analysis of Elastic Contact between Coated Bodies

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

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

Received 31 May 2018; Revised 28 September 2018; Accepted 11 October 2018; Published 1 November 2018

Academic Editor: Patrick De Baets

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

#### Abstract

Substrate protection by means of a hard coating is an efficient way of extending the service life of various mechanical, electrical, or biomedical elements. The assessment of stresses induced in a layered body under contact load may advance the understanding of the mechanisms underlying coating performance and improve the design of coated systems. The iterative derivation of contact area and contact tractions requires repeated displacement evaluation; therefore the robustness of a contact solver relies on the efficiency of the algorithm for displacement calculation. The fast Fourier transform coupled with the discrete convolution theorem has been widely used in the contact modelling of homogenous bodies, as an efficient computational tool for the rapid evaluation of convolution products that appear in displacements and stresses calculation. The extension of this technique to layered solids is tantalizing given that the closed-form analytical functions describing the response of layered solids to load are only available in the frequency domain. Whereas the false problem periodization can be treated as in the case of homogenous solids, the aliasing phenomenon and the handling of the frequency response function in origin require adapted techniques. The proposed algorithm for displacement calculation is coupled with a state-of-the-art contact solver based on the conjugate gradient method. The predictions of the newly advanced computer program are validated against existing results derived by a different method. Multiple contact cases are simulated aiming to assess the influence of coating thickness and of its elastic properties on the contact parameters and the strass state. The performed simulations prove that the advanced algorithm is an efficient tool for the contact analysis of coated bodies, which can be used to further understand the mechanical behavior of the coated system and to optimize its design.

#### 1. Introduction

The service life of machine elements undergoing contact load is expected to benefit from substrate protection by means of a protective hard coating, providing low friction, high wear resistance, and long fatigue life. The competent design of protective coatings requires the knowledge of the stresses generated in the coated system under contact load, as well as the material response to these stresses. The complex mathematical models arising in the analysis of layered systems, admitting closed-form solutions only under strong limiting assumptions, suggest the use of numerical methods, capable of treating deterministic contacting surfaces resulting from various manufacturing processes. The classical Hertz model is usually employed as a starting point in any nonconforming contact calculation. However, in case of coated bodies, due to the elastic mismatch between the coating and the substrate, the contact behavior may deviate significantly from the framework of the contact between homogenous bodies.

For a multilayered body, the closed-form relationship between the excitation (surface tractions) and the response (displacements and stresses) may be conveniently derived by the use of the integral transform theory. The latter substitutes a complicated problem in the spatial domain with its equivalent form in the transformed frequency domain, which may result much simpler and its solution easier to obtain. The mathematical modelling of stress and displacement in layered systems pioneered with the work of Burmister [1], involving Bessel type stress functions in the contact analysis of fully bonded or frictionless interfaces under prescribed axisymmetric normal load. Chen and Engel [2] extended the results to nonaxisymmetric normal loading, whereas O’Süllivian and King [3] obtained the frequency response functions (FRF) for the contact of bilayered materials under normal and shear loads. Closed-form expressions were explicitly derived in the transform domain, but solutions in the spatial domain required numerical inverse transformation. The fast Fourier transform (FFT) was first applied to Contact Mechanics by Ju and Farris [4], resulting in increased computational efficiency compared to the continuous integral transform method. The contact of rough surfaces with coatings was subsequently studied by Nogi and Kato [5] using the FFT technique and the FRFs derived by O’Süllivian and King [3]. The conjugate gradient method (CGM) was used for the solution of the contact problem for both homogenous and layered solids. Polonsky and Keer [6] discussed the periodicity error associated with the application of the FFT to nonperiodic problems and proposed a correction based on the multilevel multisummation (MLMS) technique. A method for analysis of the stress state developed in layered elastic solids with cracks was subsequently proposed [7]. The source of the periodicity error in the FFT technique was further studied by Liu et al. [8], and a method of discrete convolution and FFT (DCFFT) was advanced, which completely avoids any additional error except for the discretization error. Liu and Wang [9] applied the latter method to the study of contact stress fields caused by surface tractions. A similar technique was employed by Wang et al. [10] in the study of the partial slip contact of coated bodies. For layered materials, the explicit expressions of the FRFs stand at the core of the FFT-based techniques. Recently, Yu et al. [11] obtained the analytical FRFs of multilayered materials in a recurrence format and studied the elastic contact of layered bodies with various configurations.

This paper advances a numerical simulation technique for the contact of bodies with a single layer coating, by coupling a state-of-the-art contact solver [12] with a displacement computation technique based on the FRFs derived in the literature [5]. As the contact solver [12] was originally designed for the purely elastic response, the MLMS method employed for displacement computation requires the Green’s functions for the elastic half-space in the space domain and therefore cannot be applied to coated bodies. In this paper, the MLMS is substituted by a novel FFT-based technique that can take advantage of the existing FRFs [5]. It should be noted that Nogi and Kato [5] have also obtained the displacement in a coated body based on the FRFs, but their technique was found lacking [6] in precision due to a poorly controlled periodicity error. In this paper, the convolution products arising in the stress and displacement calculation as a result of superposition of effects of fundamental solutions are evaluated in the frequency domain with improved computational efficiency and high precision. A comprehensive set of contact simulations is performed, aiming to provide insight on the combined influence of the coating thickness and of the elastic mismatch between the coating and the substrate, on the contact parameters and on the stress state in the coated body.

#### 2. Spectral Response of Elastic Coated Half-Spaces

In the frame of linear theory of elasticity, stresses and displacements arising in a half-space due to general loadings are obtained by superposition of solutions of point forces, also referred to as the Green’s function method. The Green’s function describes the half-space response in the time/space domain, and its counterpart in the frequency domain, describing the spectral response of the half-space to a Dirac delta function, is referred to as the FRF. This duality is of particular interest in case of coated bodies, where only the FRFs are known, while the Green’s functions derivation is difficult, if not impossible. For example, the spectral normal displacement induced by a unit point force acting normally on the boundary of a half-space of Young modulus and Poisson’s ratio is where and are the angular frequencies in the and directions, respectively, both contained in bounding plane. Equation (1) is the spectral counterpart of the well-known Boussinesq solution (i.e., the Green’s function) that is usually employed to compute the displacement induced by a prescribed pressure in a homogenous, isotropic and linear elastic half-space:

The integral in (2) can be performed analytically only in a few cases of simple or axisymmetric contact geometry, such as the Hertz contact. In the general contact problem, the contact area and the pressure distribution are a priori unknown and therefore must be determined by trial-and-error. An iterative approach is usually employed, in which the initial guess is a domain expected to include the contact area. The search for the contact solution involves repeated computation of displacement for arbitrary pressure distributions such as the ones resulting during the iterative process. Such robustness can only be achieved by employing a numerical method, in which the initial guess domain is split up into elementary patches over which a prescribed pressure distribution is assumed. In this manner, the continuous pressure distribution is substituted by a reunion of pressure elements that are adjusted during the iterative process so that, at convergence, all boundary conditions are satisfied. This digitization implies that the contact problem is reported to discrete spatial coordinates instead of the continuous coordinates . The method employed in this paper assumes uniform pressure elements, i.e., a piecewise constant distribution, in which pressure is assumed uniform for each set . The most important part of the contact solution is the displacement computation, which concentrates most of the required computational resources. Once the displacement can be calculated for arbitrarily prescribed discretized pressure, using either the Green’s function or the corresponding FRF, the solution of the contact between homogenous or in homogenous bodies can be obtained in an iterative manner, as detailed in the next section. Equation (2) is in fact a two-dimensional convolution product that can also be assessed in the frequency domain.

The study of stresses and displacements in layered materials with application to layered soil deposits was pioneered by Burminster [1]. O’Sullivan and King [3] derived the stress fields in a layered half-space due to normal and shear tractions using Papkovich-Neuber elastic potentials and double Fourier transform. Based on these results, Nogi and Kato [5] later obtained the explicit form of the FRF for displacement in the three-dimensional normal contact problem involving a coated material: where is the coating thickness, , , , , , and . The shear moduli and the Poisson’s ratios for the layer and the substrate are denoted by and , respectively, with for the coating and for the substrate. It should be noted that the counterpart of (3) in the time/space domain, i.e., the Green’s function, is not available in closed-form expression.

The FRF (3) can be used to derive displacement from a superposition of effects of type (2), by employing the discrete convolution theorem [13]. The spectral computation of displacement is convenient for coated contact analysis for two reasons: (1) the displacement and stress response of multilayered materials is only available in the frequency domain as the FRFs, and (2) the convolution theorem states that, for series with terms, the computational complexity decreases from , when the convolution is evaluated in the time/space domain, to an improved in the frequency domain. The reduction comes from the fact that a convolution operation in the time/space domain can be computed as an element-wise product between the Fourier transforms of the convolution members. The application of this technique implies not only direct but also invers Fourier transform, the latter required to retrieve the convolution result in the time/space domain.

Additional difficulty arise due to application of the convolution theorem to a discretized system as opposed to the continuous one. Expression in (2) is a linear continuous convolution, and its discrete counterpart is a discrete cyclic convolution. The sampling of the convolution members creates series with a finite numbers of terms, and the transfer of these series to the frequency domain is achieved via the FFT algorithm. However, application of FFT to a series implies that the series is periodical with a period equal to the considered number of terms. In case of nonconforming contact problems, this implicit periodization infers that we are considering as excitation not the actual pressure, but a periodical load having as period the pressure distribution in the initial guess domain. This problem periodization may be desirable in case of contact simulation involving nominally flat surfaces, when the inherent microtopography can be replicated by the periodic repetition of a characteristic asperity. In case of nonperiodic contact problems, an error is thus introduced in the displacement computation. The latter periodicity error implies that we are computing the displacement induced not only by the actual pressure distribution, but by the spurious neighboring periods as well. It is expected that the periodicity error will be more significant at the edges of the computational domain, where the influence of the neighboring periods is more prominent. This assertion also suggest a way of reducing the periodicity error, by moving away the neighboring pressure periods so that their contribution to the target domain is minimized. A larger period, i.e., a larger than expected computational domain may be considered, and the original pressure must be zero-padded in the extended domain. This extension may result in series with more terms, thus increasing the computational effort. Numerical experimentation may be required to assess the optimum magnitude of the domain extension.

Nogi and Kato [5] advanced a robust algorithm for the contact of homogenous or bilayered bodies, with consideration of microtopography, in which the displacement solution was protected against false periodicity by doubling the number of grids along each direction. Liu et al. [8] discussed the source of the periodicity error in the convolution theorem and advanced the DCFFT method that can efficiently handle the conversion of linear convolution into cyclic convolution. Their strategy implies the doubling of the target domain in each physical dimension, thus obtaining a virtual domain in which the techniques of zero padding and wrap-around order are employed. A more in-depth study of the error sources, such as the aliasing and the Gibbs phenomena, is performed in [9]. The DCFFT technique requires the knowledge of the Green’s function, and therefore application to coated bodies may not be straightforward. A conversion process of the needed influence coefficients from the FRFs is indicated in [9], and a correction factor, i.e., the level of aliasing control, is employed.

In this paper, the elastic response of the coated body is calculated using an algorithm based on the computation of convolution into the frequency domain. The periodicity error is minimized by employing a target domain extension, in which the excitation is zero-padded. As described in [9], the accuracy of such a technique depends on the intensity of the aliasing phenomena, which can only be reduced by using a small sampling interval in the frequency domain. The latter can be achieved by considering a larger computational domain in the time/space domain while maintaining the sample interval (in the space domain) constant. The algorithm steps are described below. Without losing generality, the equations are presented in one dimension only (i.e., line contact) for brevity, but extension to the three-dimensional point contact is straightforward.

The algorithm requires as input the target computational domain in the spatial dimension, , and the number of pressure elements , . It is convenient to choose as a power of 2, considering that FFT is computed on a vector which is otherwise zero-padded to the next power of 2. A domain extension ratio is next chosen to minimize the periodicity error. In this paper, a value (in each spatial dimension) was imposed, the available RAM memory being the main limiting factor. The resulting extended domain should be discretized by considering grids, resulting in uniform pressure patches having the same size as the initial ones. By increasing the number of pressure elements, the sample interval in the frequency domain is decreased, aiming for the reduction of the aliasing phenomena. Indeed, by imposing a domain extension of ratio and by keeping the data interval in the spatial domain constant, the data interval in the Fourier transform domain is decreased from to .

The next step consists in the computation of the discrete frequency samples , which are the counterpart of the grids from the spatial domain. The set of discrete samples should be centered on the origin of the frequency domain:

A discrete counterpart of the FRF is computed next by evaluating in each frequency sample:

The terms of the latter series are subsequently rearranged in wrap-around order; i.e., the terms corresponding to the negative frequencies are transferred after the ones corresponding to positive frequencies. A new series of discrete spectral samples is obtained:

A different treatment is required for the second member of the convolution product, i.e., the excitation. Pressure is assumed known in the target domain as a set of pressure elements , but otherwise can be arbitrary. The expanded pressure series is obtained by zero padding the original pressure series in both directions: the original pressure samples are transferred to the middle terms of the extended series so that everywhere, except for

Then, the fast Fourier transform of the pressure series is computed:

According to the discrete convolution theorem, the element-wise product of and is computed next, producing the displacement output series in the frequency domain:

The convolution result in the spatial domain is then obtained by inverse FFT: from which only the middle terms, corresponding to the target domain, are kept as algorithm output:

The aforementioned algorithm can be used to compute the displacement or stress response of a coated half-space to a prescribed, but otherwise arbitrary, set of pressure elements acting on the boundary. Special care is required with the evaluation of the FRF at the origin of the frequency domain. Most FRFs describing the elastic response of a coated body are singular in origin; e.g., in relation (3) cannot be evaluated for . The discretization in the frequency domain implies that the FRF is assumed piecewise constant, and a representative point from the series (4) is used for each small interval. The characteristic function value on the latter interval is derived by a simple evaluation of the FRF in the chosen point. Steep gradients of the FRF, such as those appearing in the vicinity of the origin, where the function is usually singular, may challenge the choice of the representative point. The FRFs may be singular in origin; however, they are all integrable everywhere [5], including on a vicinity of the origin. Consequently, an average value of the FRF can be computed over the patch centered on origin:

The latter average value is the best available indicator of the local behavior of the FRF in a vicinity of origin and therefore can be used instead of the noncalculable . As discussed in the following section, the value of the FRF in origin may not even be needed for the solution of the contact problem involving coated bodies.

#### 3. The Contact of Coated Bodies

In the framework of Contact Mechanics, the nonconforming contact between idealized geometries is preponderantly analyzed as its nature is well-defined, easily controlled, repeatable and relatively insensitive to manufacturing imperfections. In this case, the contact stresses arise as a local stress concentration that can be assumed independent of stresses in the bulk of the bodies. Consequently, bodies of arbitrary surface may be considered as linear elastic half-spaces, and the results from the previous section can be directly applied to computation of stresses and displacements generated in the contacting bodies.

The significance of the term is discussed first, for the case when is the FRF describing the spectral displacement due to a normal load described by a Dirac delta function (a unit point normal force). The spectral series is implicated in displacement computation after rearrangement in wrap-around order according to relation (6). After this rearrangement, becomes the first term in the new series . When inverse FFT is applied, as in relation (10), to a series whose first term is perturbed by a constant value , each term of the resulting series is perturbed by the quantity , where is the number of terms in the series. Moreover, if the term is misevaluated by a constant, only the first term of the displacement series, , will be compromised, as is calculated by element-wise product. Consequently, after application of inverse FFT to , the displacement series in the space domain, , will be known except for a constant, related to the misevaluation of .

The situation in which only relative displacements of the loaded half-space boundary can be computed is well known in the Contact Mechanics theory from the framework of the plane contact problem, i.e., the line contact. In the latter case, this model inability is a consequence of the way the plane problem solution is derived [14]. For a three-dimensional problem (i.e., point contact), the solution is unique, as shown by the displacement solution calculated by Love [15]. However, the strategy designed [14] for the resolution of the plane contact problem, making use of relative displacements only, can be extended to the three-dimensional contact case.

The model for the frictionless contact used in this paper is based on a well-known formulation employed in contact scenarios with various constitutive laws: elastic [12, 16, 17], elastic-plastic [18], or viscoelastic [19–21]. For nonconforming contacts, the initial point of contact evolves into a contact surface as normal load is applied, and the elastic bodies deform as a result of compression. The equation of the surface of separation between the deformed contacting bodies results [14] from comparison of contact geometry before and after the deformation, along the normal contact axis (i.e., the -axis, passing through the initial point of contact and normal to the plane that separates best the contacting surfaces in the latter point of contact): where denotes the gap between the deformed surfaces, the initial (i.e., prior to deformation) gap, and is the approach between points in the contacting bodies distant to the contact region (i.e., the rigid-body motion). A contact process schematic is shown in Figure 1. The boundaries of the two contacting bodies (a) and (b) in underfomed state are displayed using dashed lines. The normal approach is indicated as the distance between the centers of the positions occupied by the spherical indenter before and after the half-space deformation. Without loosing generality, body (a) is assumed rigid and body (b) is an elastic coated half-space, for simplicity. In this manner, the superposition of initial contact geometry and of normal displacement, and , takes a simpler form, as and . The contact model is also valid when the indenter is elastic (coated or not). In this case, the corresponding displacements must be computed separately and added to to obtain the composite displacement needed in (13). The displacement computation described in the previous section assumes that the elastic contacting bodies can be approximated with elastic half-spaces.