Receding Contact Problem of Multi-Layered Elastic Structures Involving Functionally Graded Materials

: This paper studies a receding contact problem of a functionally graded layer laminate pressed against a functionally graded coated homogeneous half-plane substrate by a rigid ﬂat indenter. The shear modulus of the functionally graded materials with a constant Poisson’s ratio is modeled by an exponential function which varies along the thickness direction. Both the governing equations and the boundary conditions of the receding contact problem are converted into a pair of singular integral equations using the Fourier integral transforms, which are numerically integrated by the Chebyshev–Gauss quadrature. The contact pressure and the contact size at both contact interfaces are eventually obtained iteratively, as developed from the steepest descent algorithm. Extensive parametric studies suggest that it is possible to regulate the contact pressure and contact size by constructing the top layer from a soft functionally graded material.


Introduction
For a multi-layered structure, according to the variation of the contact size between different layers under an external load, three contact types can be identified: advancing contact, if the contact size becomes larger; Stationary contact, if the contact size remains unchanged; and receding contact, if the contact size becomes smaller. Among these, receding contact has not been well studied, since it rarely occurs in a single-layer structure. For receding contact, however, a smaller contact size means more concentrated pressure. Hence, for multi-layered structures, receding contact may give rise to several adverse effects, including layer peeling, delamination fracture, etc.
A variety of methods have been adopted to study the receding contact problems of multi-layered homogeneous elastic structures. By using Fourier transform method, Comez et al. [1] studied the receding contact problem of double elastic layers indented by a rigid cylindrical indenter. On the basis of this work, Cakiroglu et al. [2] introduced the artificial neural networks method to reproduce the problem. By using Fourier transform method, Adibelli et al. [3] studied the receding contact problem between a coated layer and a half-plane loaded by a rigid cylindrical indenter. Yaylaci and Birinci [4] analyzed a double receding contact problem between two elastic layers and two elastic quarter-planes by using an integral transform technique, further verified by the finite element method [5] and the artificial neural network method [6].
Functionally graded materials (FGM) are generally formed from several inhomogeneous particulate composites in such a way that the volume fractions of the constituents vary along any desired spatial direction. Special preparations can yield FGM gradient changes along the thickness or in any other directions. In theoretical studies, linear functions [7,8], exponential functions [9,10], and power functions [11,12] are all used to fit the variations of the material properties approximately. Among these, the exponential function

Materials and Methods
In a multi-layered structure like that shown in Figure 1, the contact surface lying between the FG layers may be partially separate, while a rigid flat indenter with a halflength presses into the top layer. Figure 1 illustrates a typical scenario in some industry applications. Taking the punching machine as an example, the upper FG layer represents an FG gasket, and the lower FG layer represents the FG-coatings established for abating the shock effect on the substrate. The thicknesses of the FG layers are respectively h 2 and h 3 . The shear modulus of the half-plane substrate is treated as a constant µ 1 , and both the FG layers are modeled by exponential functions varying along their thickness dimensions: where β 2 and β 3 are arbitrary non-zero constants, and µ 03 represents the shear modulus of the top FG layer's bottom surface. The Poisson's ratios of the FG layers are both fixed as the same constant ν, since the value can hardly influence the contact substantially [44]. To simplify the computation, this work also ignores the influences of the gravitational and frictional forces. Under the proposed assumptions, each layer's boundary condition can be enumerated as: where H denotes the Heaviside step function, and σ xxj , σ xyj and σ yyj represent the pressure components of the j-th layer (j = 1, 2, 3) respectively, p(x) represents the distribution function of the stationary contact pressure, and q(x) represents the distribution function of the receding contact pressure, u j and v j represent the displacement components of the j-th layer (j = 1, 2, 3). For the absence of the frictional forces, the vertical displacement certainly stays continuous in both contact areas. To ensure the continuity of the vertical displacement and exclude the unwanted effects of rigid-body displacement, this work quantifies the profile  For the absence of the frictional forces, the vertical displacement certainly stays continuous in both contact areas. To ensure the continuity of the vertical displacement and exclude the unwanted effects of rigid-body displacement, this work quantifies the profile of the rigid indenter by its slope function f (x) = 0. Taken together, these conditions indicate that: According to the basic formulas of mechanics (the equilibrium equation, the straindisplacement relationships, and the linear elastic stress-strain law), the following equations concerning the displacements of layers are satisfied: where β 1 = 0, the above equations are suitable for homogeneous half-plane; Kolosov's constant κ is related to the Poisson's ratio by κ = 3 − 4ν for plane strain and κ = (3 − ν)/(1 + ν) for plane stress.
To address the formal solutions for the displacements and pressure components, this study maps the boundary conditions (3)(4)(5)(6)(7)(8)(9)(10)(11)(12) into Fourier space using Fourier integral transform. The details can be obtained in our previous work [45]. By substituting the displacement components' integral expressions into Equations (13) and (14), the following system of singular integral equations can be obtained by the inverse Fourier transformation: The expressions of K j (x, t)(j = 1, 2, 3, 4) are documented in Appendix A. The unknowns b and the distribution function of the contact pressure should satisfy the equilibrium conditions of the force expressed by the following equations: where P is the known resultant force applied to the flat indenter.

Solution to Dual Singular Integral Equations
To simplify the computation, all the physics qualities are non-dimensionalized by introducing the following transformations: Thus, the equilibrium conditions Equations (19) and (20) can be defined as: Meanwhile, the singular integral Equations (17) and (18) can be simplified as: Following the method of Erdogan and Gupta [46], Equation (33) is of the type with singularity index "1", since the singular behavior of the unknown function p 1 (s) appears while s = ±1, and Equation (34) is of the type with singularity index "−1", since the unknown function q 1 (s ) is continuous and bounded throughout the closed interval |s | ≤ 1. Thus, the flat rectangular indenter can be distinguished from the smooth indenters with differentiable shapes, including circular indenters, for the singularities existed in the tip of the flat rectangular indenter. Therefore, traditional methods for solving the singular integral equations for the receding contact problem with a circular indenter or a uniformly distributed load are inadequate here. This work suggests an innovative discretizing method for addressing the singularities as follows.
First, the integral portion that contains the top layer stationary contact pressure can be discretized using the Gauss-Chebyshev integral formulas of the first kind, and the integral portion that contains the bottom layer receding contact press can be discretized using the Gauss-Chebyshev integral formulas of the second kind. At the same time, Equations (31) and (33) are solved with Chebyshev polynomials of the second kind, and Equations (32) and (34) are solved with Chebyshev polynomials of the first kind. Now we can find: in which: Obviously, there are N + 1 equations derived from Equation (36), which correspond to N + 1 different unknown variables r j . Specifically, N is an even number. The (N/2 + 1)th equation is established automatically [47], and can be safely deleted. Thus, Equations (35)- (38) only contains 2N + 1 equations, as well as 2N + 1 unknowns, i.e., the contact pressure p 2 (s k ), q 2 (s k ), (k = 1, 2, · · · , N) and the half-length of receding contact b. It is worth noting that, among all the unknowns in Equations (35)- (38), only b is non-linear. We can initialize b with a random value and substitute it into Equations (35)- (38); the remaining unknowns p 2 (s k ) and q 2 (s k ), (k = 1, 2, · · · , N) can then be solved iteratively. The numerical results of the iteration method can be verified by Equation (38), if the absolute error δ satisfies:

Numerical Results and Discussion
In all numerical experiments, Poisson's ratios are fixed as ν = 0.25, (κ = 2). Several well-studied models can be expressed by the presented study if the key parameters are modified as specific values. For instance, this study can be reduced to a receding contact problem between a homogeneous layer and a half-plane under a concentrated force, if the parameters satisfy: β 3 h 3 = 0.001, h 2 = 0 and a/h 3 = 0.01. Our simulation results show that the corresponding normalized receding contact half-length is 1.3243, which is almost identical to the results reported by the early literature [48,49]. The accuracy of the numerical results can thus be verified.
This work first investigates the receding contact distributions for different indenter sizes. Both Figures 2 and 3 show positive relationships between the indenter size and the stationary/receding contact half-length, which means a larger indenter size invariably leads to a larger stationary/receding contact half-length. Correspondingly, the minimum stationary contact pressure and the maximum receding contact pressure are both negatively correlated with the indenter size. Comparing the position in which the extreme value appears, it can be determined that the minimum stationary contact pressure always appears at the center of the contact area, while the maximum receding contact pressure does not deterministically appear at the similar position. For several cases, a pressure concave Crystals 2022, 12, 354 7 of 14 may appear at the center of the receding contact area, and the larger the receding contact half-length, the more obvious the pressure concave. stationary/receding contact half-length, which means a larger indenter size invariably leads to a larger stationary/receding contact half-length. Correspondingly, the minimum stationary contact pressure and the maximum receding contact pressure are both negatively correlated with the indenter size. Comparing the position in which the extreme value appears, it can be determined that the minimum stationary contact pressure always appears at the center of the contact area, while the maximum receding contact pressure does not deterministically appear at the similar position. For several cases, a pressure concave may appear at the center of the receding contact area, and the larger the receding contact half-length, the more obvious the pressure concave.     Figure 4 which describes the correlation between the inhomogeneity parameter of the top layer and the non-dimensionalized receding contact half-length. Along with the increasing of the inhomogeneity parameter, the non-dimensionalized receding contact half-length first decreases, and then increases, while the nondimensionalized receding contact half-length monotonically decreases for the inhomogeneity parameter of the interlayer in Figure 5. Specially, Figure 5 shows the receding contact half-lengths between the homogeneous elastic layer and the FG coated homogeneous half-plane due to the flat indenters. Comparing Figure 5 with the 7th figure of reference [45], we can also investigate the interaction between the shape of the indenter and the inhomogeneity parameter of the FG coating. Regardless of shape (circular or flat) of the indenter, the receding contact half-length decreases with the increases in the inhomogeneity parameters of the FG coating, as well as the decrease in the indenter size. However, under the flat indenter, the regulations of the inhomogeneity parameters of the FG coating  Figure 4 which describes the correlation between the inhomogeneity parameter of the top layer and the non-dimensionalized receding contact half-length. Along with the increasing of the inhomogeneity parameter, the non-dimensionalized receding contact half-length first decreases, and then increases, while the non-dimensionalized receding contact half-length monotonically decreases for the inhomogeneity parameter of the interlayer in Figure 5. Specially, Figure 5 shows the receding contact half-lengths between the homogeneous elastic layer and the FG coated homogeneous half-plane due to the flat indenters. Comparing Figure 5 with the 7th figure of reference [45], we can also investigate the interaction between the shape of the indenter and the inhomogeneity parameter of the FG coating. Regardless of shape (circular or flat) of the indenter, the receding contact half-length decreases with the increases in the inhomogeneity parameters of the FG coating, as well as the decrease in the indenter size. However, under the flat indenter, the regulations of the inhomogeneity parameters of the FG coating on the receding contact half-length are less efficient than those under the circular indenter.
sionalized receding contact half-length first decreases, and then increases, while the nondimensionalized receding contact half-length monotonically decreases for the inhomogeneity parameter of the interlayer in Figure 5. Specially, Figure 5 shows the receding contact half-lengths between the homogeneous elastic layer and the FG coated homogeneous half-plane due to the flat indenters. Comparing Figure 5 with the 7th figure of reference [45], we can also investigate the interaction between the shape of the indenter and the inhomogeneity parameter of the FG coating. Regardless of shape (circular or flat) of the indenter, the receding contact half-length decreases with the increases in the inhomogeneity parameters of the FG coating, as well as the decrease in the indenter size. However, under the flat indenter, the regulations of the inhomogeneity parameters of the FG coating on the receding contact half-length are less efficient than those under the circular indenter.  For the monotony correlation exhibited in Figure 5, the way in which the inhomogeneity parameter of the interlayer 22 h  influences the contact pressure should be investigated further. Figure 6 illustrates the distribution of the stationary contact pressure located between the indenter and the FG layer. The stationary contact half-length depends only on the indenter size; any other parameters cannot influence it. Figure 6a shows the minimum contact pressure achieved at the center of the stationary contact area. At the edge of the area, the contact pressure increases sharply for the singularity, which occurs at the border of the indenter. It is difficult to estimate the influence of the inhomogeneity parameter on the stationary contact pressure since the contact pressure at the edge area is too large. The details are shown in Figures 6b,c. A noteworthy phenomenon is also indicated by Figures 6b,c: a positive relationship emerges between the inhomogeneity parameter of the interlayer and the contact pressure of the stationary contact center area, while the relationship near the edge of the stationary contact is the reverse. Figure 7 shows the distribution of the non-dimensionalized receding contact pressure for different inhomogeneity parameters of the interlayer. There is an obvious inverse correlation between the maximum value of the receding contact pressure and the receding contact half-length. The influence of the inhomogeneity parameter 22 h  on the receding contact pressure is similar to that on the stationary contact pressure. The maximum value For the monotony correlation exhibited in Figure 5, the way in which the inhomogeneity parameter of the interlayer β 2 h 2 influences the contact pressure should be investigated further. Figure 6 illustrates the distribution of the stationary contact pressure located between the indenter and the FG layer. The stationary contact half-length depends only on the indenter size; any other parameters cannot influence it. Figure 6a shows the minimum contact pressure achieved at the center of the stationary contact area. At the edge of the area, the contact pressure increases sharply for the singularity, which occurs at the border of the indenter. It is difficult to estimate the influence of the inhomogeneity parameter on the stationary contact pressure since the contact pressure at the edge area is too large. The details are shown in Figure 6b,c. A noteworthy phenomenon is also indicated by Figure 6b,c: a positive relationship emerges between the inhomogeneity parameter of the interlayer and the contact pressure of the stationary contact center area, while the relationship near the edge of the stationary contact is the reverse.    Figure 7 shows the distribution of the non-dimensionalized receding contact pressure for different inhomogeneity parameters of the interlayer. There is an obvious inverse correlation between the maximum value of the receding contact pressure and the receding contact half-length. The influence of the inhomogeneity parameter β 2 h 2 on the receding contact pressure is similar to that on the stationary contact pressure. The maximum value of the receding contact pressure increases with the increase in the inhomogeneity parameter β 2 h 2 .   Figure 8 illustrates how the shear moduli of both FG layers' bottom surfaces regulate the normalized receding contact half-length. For any inhomogeneity parameter value, a larger shear modulus value of the top layer's bottom surface always corresponds to a smaller receding contact size and to a larger maximum value of receding contact pressure as well. By combining Figures 4 and 5, we observe an interesting effect: while β 2 h 2 = β 3 h 3 > −1, the influence on the receding contact half-length caused by the inhomogeneity of the top FG layer is larger than that caused by the inhomogeneity of the interlayer. 1 hh  =  − , the influence on the receding contact half-length caused by the inhomogeneity of the top FG layer is larger than that caused by the inhomogeneity of the interlayer. , the proposed model is reduced to a receding contact problem between a homogeneous layer and a homogeneous half-plane. The thickness of the homogeneous interlayer does not significantly influence the contact area. This indicates that a harder and thicker FG interlayer results in a large amount of receding contact, and therefore, a reduced peak contact pressure in Figure 9. The opposite is true for a soft FG interlayer. The thicker the soft FG interlayer, the smaller the receding contact halflength. This also leads to an increase in the maximum contact pressure. The thickness of the interlayer can also regulate the receding contact half-length. Figure 9 focuses on the influence of the interlayer's thickness on the receding contact halflength, while the inhomogeneity parameters of the top layer and the interlayer are equal. For a special case β 2 h 2 = β 3 h 3 = 0.001, the proposed model is reduced to a receding contact problem between a homogeneous layer and a homogeneous half-plane. The thickness of the homogeneous interlayer does not significantly influence the contact area. This indicates that a harder and thicker FG interlayer results in a large amount of receding contact, and therefore, a reduced peak contact pressure in Figure 9. The opposite is true for a soft FG interlayer. The thicker the soft FG interlayer, the smaller the receding contact half-length. This also leads to an increase in the maximum contact pressure. 1 hh  =  − , the influence on the receding contact half-length caused by the inhomogeneity of the top FG layer is larger than that caused by the inhomogeneity of the interlayer. , the proposed model is reduced to a receding contact problem between a homogeneous layer and a homogeneous half-plane. The thickness of the homogeneous interlayer does not significantly influence the contact area. This indicates that a harder and thicker FG interlayer results in a large amount of receding contact, and therefore, a reduced peak contact pressure in Figure 9. The opposite is true for a soft FG interlayer. The thicker the soft FG interlayer, the smaller the receding contact halflength. This also leads to an increase in the maximum contact pressure. Figure 9. The dependence of non-dimensionalized receding contact half-length on the thickness ratio h 2 /h 3 for different inhomogeneity parameters β 2 h 2 and β 3 h 3 (µ 03 = µ 1 , a/h 3 = 0.5). Figure 10 focuses on the influence of the interlayer's thickness on the receding contact half-length, while the inhomogeneity parameters of the top layer and the interlayer are not equal. However, the combined effect of the thickness of the FG interlayer and the inhomogeneity parameters is difficult to determine from Figure 10. The manifested features in Figure 10 are similar to those in Figure 9, but the combined effect of the thickness of the FG interlayer and the inhomogeneity parameters is still ambiguous in this figure.  Figure 10 focuses on the influence of the interlayer's thickness on the receding contact half-length, while the inhomogeneity parameters of the top layer and the interlayer are not equal. However, the combined effect of the thickness of the FG interlayer and the inhomogeneity parameters is difficult to determine from Figure 10. The manifested features in Figure 10 are similar to those in Figure 9, but the combined effect of the thickness of the FG interlayer and the inhomogeneity parameters is still ambiguous in this figure.

Conclusions
In this study, the receding contact problem of a multi-layered elastic structure involving FGM under a flat indenter is considered. The model is transformed into dual singular integral equations that are then solved numerically. From comprehensive parametric research studying the inhomogeneity parameters of FGM, the indenter size, the interlayer's thickness, and the material properties, we can finally come to a number of noteworthy conclusions: The stationary contact size is only related to the flat indenter size. A wider flat indenter produces a larger stationary contact size. As a result, the pressure in the stationary contact center area is reduced.
The center of the receding contact area sustains the maximum pressure while the receding contact half-length is small. However, when the half-length is large enough, a special pressure concave may appear at the center of the contact area. The maximum receding contact pressure is inversely related to the receding contact half-length; the longer the half-length, the deeper the pressure concave.
The receding contact half-length is inversely correlated with the inhomogeneity parameter of the interlayer and the shear modulus ratio is

Conclusions
In this study, the receding contact problem of a multi-layered elastic structure involving FGM under a flat indenter is considered. The model is transformed into dual singular integral equations that are then solved numerically. From comprehensive parametric research studying the inhomogeneity parameters of FGM, the indenter size, the interlayer's thickness, and the material properties, we can finally come to a number of noteworthy conclusions: The stationary contact size is only related to the flat indenter size. A wider flat indenter produces a larger stationary contact size. As a result, the pressure in the stationary contact center area is reduced.
The center of the receding contact area sustains the maximum pressure while the receding contact half-length is small. However, when the half-length is large enough, a special pressure concave may appear at the center of the contact area. The maximum receding contact pressure is inversely related to the receding contact half-length; the longer the half-length, the deeper the pressure concave.
The receding contact half-length is inversely correlated with the inhomogeneity parameter of the interlayer and the shear modulus ratio is µ 03 /µ 1 ; it is positively correlated with the indenter size.
The interlayer thickness can regulate the receding contact area, but the regulation effect is also influenced by the inhomogeneity parameters β 2 h 2 and β 3 h 3 .

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.