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.
The discretization of pressure is extended to all problem parameters, and two indexes are used for spatial localization inside a rectangular computational domain divided into small nonintersecting rectangles. The domain is the initial guess in an iterative process aiming to find, by trial-and error, the contact area , defined as the reunion of patches having associated nonvanishing pressure elements. The contact area must be contained in the computational domain at any iteration, otherwise the searching sequence must be restarted with a larger guess. An additional assumption is required for the convergence of this iterative process: where the strict inequality applies to the contact area. Relation (14) implies adhesion is not accounted for. Whereas the latter may play an important role in the contact or rubber-like materials, the real contact area between metallic bodies, established between the higher asperities, is much smaller than the apparent contact surface, so that adhesion may be assumed trivial.
As bodies are assumed impenetrable in the frame of linear theory of elasticity, where the strict inequality applies to the noncontact region of . The boundary conditions (14) and (15) can be streamlined as the latter equation being used to assess, in the iterative process, the contact or noncontact status of every patch in the computational domain . The contact model is completed with the static force equation, relating the normal force to the pressure distribution:
The model also assumes that the resultant of the pressure distribution is aligned with the normal force, so that no tilting moment or rolling is expected to occur. The contact model (13)–(17) describes any instantaneous contact state and therefore can be used beyond the framework of linear elasticity, to simulate the contact process of bodies with complicated response (i.e., elastic-plastic, viscoelastic, inhomogeneous materials), provided the displacement is properly evaluated. The difficulty in solving the latter system of equations and inequalities stems from the fact that neither the contact area, nor the pressure distribution, are known in advance. Relation (13) written for the contacting patches yields a system of equations having the pressure elements as unknowns. But the size of this system is also unknown, because the contact area is unknown. The solution is to determine the contact area (i.e., the size of the system), and the pressure elements (i.e., the system solution) simultaneously. The existence and uniqueness of the solution of system (13)–(17) was discussed in [22], and convergence is guaranteed for strictly compressive contact tractions, i.e., is not allowed.
The system resulting from (13) is essentially linear in the pressure elements and can be solved using appropriate numerical techniques, such as the conjugate gradient method (CGM). The advantage of the latter, beside the superlinear rate of convergence, consists in the possibility of adding new constraints during the iterative process. Indeed, it is convenient not to include the static equilibrium equation (17) in the system, but to enforce outside the CGM core, thus conserving the symmetrical nature of the system matrix, as required by the CGM. It should be noted that the numerical solution of the CGM sequence is achieved by minimization of the residual , . The size of the system is iterated based on the boundary conditions (14)–(16). The patches for which, during system resolution, pressure results negative, , are excluded from the contact area, and the size of the system decreases. On the other hand, patches , for which the gap results negative, , are included in the contact area, and the size of the system increases. The stabilization of this apparently opposing processes leads to convergence, which is established based on the ratio of the change in two subsequent pressure iterations to the imposed load.
Apparently, (13) may be compromised by a misevaluation of displacement . Let be the constant introduced by setting the FRF in origin to zero, i.e., . The true value of displacement and the computed one, , are then related by
The rigid-body approach may be expressed in terms of displacement and initial contact geometry from (13), considering that , . Consequently,where denotes the numbers of contacting patches, i.e., the cardinal of the set . By combining relations (18) and (19), the relation between the true and the calculated rigid-body approach results as
Considering (13), (18), and (20), it follows that the system residual remains unchanged when is used instead of . This algorithm feature suggest that pressure distribution and the contact area in the three-dimensional contact problem can be calculated by making use of relative displacements only. For computational purposes, one may prefer to take a shortcut by setting the value of to zero, thus avoiding potential convergence issues with the numerical integration in (12). Whereas the literature [5, 9] suggests 64-point Gaussian quadrature, we encountered no convergence problem with the Matlab function “quad2d”. However, the related computational effort often exceeds that for the evaluation of all the other terms in the series. This shortcut is nonexistent when the rigid-body approach is needed (requiring absolute displacement evaluation), or in the computation of stresses developing in the contacting bodies due to contact load.
The flowchart of the contact solver employed to derive contact area and pressure distribution is presented in Figure 2. A complete description of the algorithm is beyond the point of this paper, and an interested reader is directed to the literature [12, 18, 21]. By substituting in (13) the displacement calculated with the numerical technique (4)–(11), a numerical model for the contact of coated bodies, bounded by arbitrary rough surfaces, becomes readily available.
4. Results and Discussions
The frictionless, quasistatic contact of a rigid spherical indenter pressed against an elastic coated half-space is simulated using the proposed numerical technique. The parameters of the elastic medium, that is, the Young modulus and the Poisson’s ratio, are denoted by and , with for the semi-infinite substrate and for the dissimilar elastic coating that is perfectly bonded to the substrate. The reference scenario is the Hertz contact for the homogenous elastic medium, i.e., when and , yielding the contact area , the central pressure and the rigid-body approach , which are used as normalizers for the matching contact parameters. The magnitude of the normal load, the initial gap, the Young modulus of the substrate, , the Poisson’s ratios of the elastic materials, , are kept constant in all simulations, whereas the coating thickness and the ratio are varied. Contact parameters and stresses are derived numerically based on a rectangular target computational domain of side lengths . The displacement and stress computation was performed, however, in a domain 8 times larger in each tangential direction, according to the proposed value of for the minimization of the periodicity error. The calculation of the normal approach and of the stress state in the elastic medium requires FRF evaluation at the origin of the frequency domain, and the Matlab function “quad2d” was employed with this task. No convergence errors were encountered by using the Matlab default absolute and relative tolerances. The CGM based contact solver also converged in all simulations, with an imposed precision chosen with respect to the number of grids in the target domain, as suggested in [12]. It should be noted that the derivation of contact parameters requires a two-dimensional surface grid only, whereas the stress computation involves grid extension in the depth direction. Bases on a convergence analysis, a mesh was found to provide a reasonable compromise, assuring a computational time of less than 10 minutes (for each contact case) on a 4 core 3.2GHz CPU, with a peak RAM memory utilization of 8 GB. Contact parameters only can be derived on a much finer two-dimensional mesh, e.g., with grids in the target computational domain, in a matter of seconds. The stress state computation requires the derivation of the nodal values of each of the six stress tensor components, hence the important computational effort. The propensity of plastic yielding can be assessed based on the relation between the magnitude of the von Mises equivalent stress, , and the yield strength of the elastic material, :
The contact simulation program is first validated against existing results obtained by a different method. The algorithm for the computation of elastic displacement in a layered half-space is authorized by the pressure profiles depicted in Figure 3, matching the literature results [3]. The Hertz semielliptic pressure distribution, corresponding to the case , matches well the prediction of the numerical program. The pressure profiles in the coated contact can deviate significantly from semielliptic for specific coating configurations, as shown in Figure 4. Thin coatings harder than the substrate lead to a more evenly distributed pressure on a receding contact area.
The understanding and prediction of the indentation of a layered half-space was the subject of numerous research efforts [23–25], both analytical and numerical. Yang and Ke [23] investigated the two-dimensional frictionless contact problem of a coating structure consisting of a surface coating, a functionally graded layer and a substrate under a rigid cylindrical punch. The contact problems for thin films and cover plates, in which the loading consists of any one or combination of stresses caused by uniform temperature changes and temperature excursions, far field mechanical loading, and residual stresses resulting from film processing or welding, was considered by these authors [24]. Volkov et al. [25] obtained the fields of displacements, stresses and strains arising in a linearly elastic half-space with a functionally graded coating subjected to indentation by a punch with a spherical tip.
The methodology applied in this paper differs from the aforementioned works, and the goals of the paper are different. The newly proposed numerical algorithm can address the layered contact problem for arbitrary contact geometry, and the improved computational efficiency allows for high density meshes that can tackle rough contact scenarios. To our best knowledge, 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 field, has not yet been analyzed. Extrapolation of design recommendations based on the limited number of specific cases presented in the literature may be difficult, especially considering the complexity of the involved dependencies. A parametric study is therefore performed using the newly proposed method, aiming to advance coating thickness and compliance recommendations, based on the numerically predicted contact parameters, as well as on the calculated depth and intensity of the von Mises equivalent stress. From a mechanical point of view, these parameters are of paramount importance for the improvement of the contact resistance.
The combined influence of the coating thickness and of the elastic mismatch between the coating and the substrate , on the contact parameters, is depicted in Figures 5–10. To this end, the dimensionless coating thickness was incremented in the interval with a step of 0.025, and the ratio was varied between 0.2 and 4 with an increment of 0.1, resulting in contact simulations. The predicted iso-contours of central maximum pressure, contact radius and rigid-body approach, are depicted in Figures 5, 7, and 9, respectively. Particular curves showing the influence of coating thickness on the contact parameters for specific ratios are presented in Figures 6, 8, and 10.
Figures 5–8 suggest that the more compliant the coating with respect to the substrate, the larger the contact radius and the shallower the maximum central pressure. For , the contact radius increases with the coating thickness, whereas the maximum pressure decreases. The opposite is true for . For coatings thick and stiff, i.e., and , the pressure distribution differs significantly from the Hertz case, the central pressure being at least twice its reference counterpart. Figure 8 exposes a limitation of the numerical approach: the contact area can only vary in increments equal to the sampling interval, hence the steps in the plotted data.
Figures 9 and 10 suggest that the rigid-body approach increases with the coatings thickness for , whereas the opposite is true for . In the latter case, the stiffer the coating, the more rapid the decrease of with the increase of .
The stress state analysis aims to assess the influence of the input parameters on the dimensionless magnitude and depth of the maximum von Mises equivalent stress. The use of relation (12) is first benchmarked against results from the literature [3] and a good agreement is found. Position of maximum equivalent stress is denoted by the symbol “X” in Figures 11–13. For comparison purposes, it should be noted that in [3], the stress tensor second invariant is plotted, with . Stress iso-contours in Figure 11 match well the analytical Hertz solution, giving confidence in the stress calculation method.
Figures 12 and 13 show that the von Mises equivalent stress is discontinuous across the interface between the coating and the substrate, despite the continuity of displacements and of the stress tensor components , , and . The maximum von Mises stress is always located on the contact axis. The FRFs used in the stress state computation are given in Appendix.
The effect of varying and on the intensity and depth of the maximum von Mises stress is also assessed by numerical simulation of the aforementioned contact cases. Relevant magnitude and depth iso-contours are depicted in Figures 14 and 15, whereas in Figure 16 the depth of the maximum stress is normalized by the layer thickness instead of . In the latter plot, a smaller than unity value indicates that the maximum is located in the coating. For these cases, the contact performance could greatly benefit from the coating process, considering that the coating material can be chosen with a yield strength much greater than that of the substrate, thus impeding plastic yielding and crack nucleation which eventually lead to contact failure by fretting fatigue. The iso-contours in Figure 16 show that in a domain roughly overlapping the rectangle of sides and , the maximum von Mises stress is located in the coating, but very close to the interface, i.e., .
5. Conclusions
The competent design of protective coatings, requiring detailed knowledge of stresses developing in a coated contact of complex initial geometry, may benefit from the use of numerical techniques. The latters are expected to provide reliable simulation data for real engineering applications, by solving models with fewer limiting assumptions.
The numerical method advanced in this paper combines a state-of-the-art contact solver, based on the conjugate gradient method, with a novel technique for computation of displacement in coated half-spaces undergoing prescribed normal load. In the newly proposed method, the displacement computation is performed in the Fourier transform domain for two reasons: (1) at the time of the analysis, the response of the coated half-space to a unit normal force is only available in the frequency domain, as the frequency response functions, and (2) calculation of convolution products in a discrete system can be performed more rapidly in the frequency domain than in the time/space domain.
Application of the convolution theorem to contact problems, which are essentially nonperiodic unless explicitly considered so, may be challenged by the periodicity error. In this paper, a domain extension in the space domain is employed, seconded by an increase the resolution in the frequency domain to minimize the aliasing phenomenon. Additional difficulty arises in the evaluation of the frequency response function at the origin of the frequency domain. In the newly proposed method, this computation was performed numerically, given that the frequency response function is integrable in a neighborhood of the origin, and consequently an average value can be calculated for the latter vicinity. When required, this mean value was considered instead of the missing sample in the discrete series resulted from digitization of the frequency response function. The advantages of the proposed methodology consist in (1) the ability to treat prescribed, yet arbitrary, contact geometry, including the inherent microtopography of real contacting surfaces; (2) increased computational efficiency due to the calculation of convolution products in the frequency domain; and (3) efficient control of periodicity error and of the aliasing phenomenon.
The predictions of the newly advanced numerical program agree well with results from the literature obtained by other methods. A parametric study is performed, aiming to assess 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 developed in the coated body. These results may provide basis for the design of highly efficient coatings that improve the contact performance, thus extending the service life of various machine elements undergoing contact loading.
Appendix
The stresses and displacements in the coated system can be expressed in Cartezian coordinates , , and , with the origin of the -axis for each layer located on its top. In this paper, is used for the coating and for the substrate. The FRFs of the six stress tensor components are functions of , , and in the coating and of , , and in the substrate, where and are the coordinates in the frequency domain corresponding to and , respectively, where the needed coefficients are
Data Availability
The expressions of the frequency response functions 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 cofunded from the European Regional Development Fund.