On Boussinesq’s Problem for a Power-Law Graded Elastic Half-Space on Elliptical and General Contact Domains

The indentation of a power-law graded elastic half-space by a rigid counter body is considered in the framework of linear elasticity. Poisson’s ratio is assumed to be constant over the half-space. For indenters with an ellipsoidal power-law shape, an exact contact solution is derived, based on the generalizations of Galin’s theorem and Barber’s extremal principle for the inhomogeneous half-space. As a special case, the elliptical Hertzian contact is revisited. Generally, elastic grading with a positive grading exponent reduces the contact eccentricity. Fabrikant’s approximation for the pressure distribution under a flat punch of arbitrary planform is generalized for power-law graded elastic media and compared with rigorous numerical calculations based on the boundary element method (BEM). Very good agreement between the analytical asymptotic solution and the numerical simulation is obtained for the contact stiffness and the contact pressure distribution. A recently published approximate analytic solution for the indentation of a homogeneous half-space by a counter body, whose shape slightly deviates from axial symmetry but is otherwise arbitrary, is generalized for the power-law graded half-space. The approximate procedure for the elliptical Hertzian contact exhibits the same asymptotic behavior as the exact solution. The approximate analytic solution for the indentation by a pyramid with square planform is in very good agreement with a BEM-based numerical solution of the same problem.


Introduction
Functionally graded materials (FGMs) are widespread amongst tribological systems in biology [1], biotechnology [2], and technics [3]. The main reason for that is their enhanced ability to withstand various forms of mechanical damage associated with contact and deformation [4]. The analysis of mechanical contacts of FGMs started in the framework of geo-and soil mechanics [5], as early as in the 1940s. Only decades later, the tribological perspectives of FGMs became a focus of research [6], which sparked a broad series of works on their contact mechanical properties (for a comprehensive, up-to-date overview, see the handbook [7], p. 251 ff.).
Meanwhile, tribology and materials science of FGMs and functionally graded coatings have developed into large scientific branches within their respective fields, with a very broad corresponding literature (see, e.g., [8][9][10]). However, in the present manuscript, we will concern ourselves only with the contact mechanics of FGMs.
In that regard, several forms of functional material inhomogeneity have been considered in the literature [11]. Nevertheless, only a few cases allow for a closed-form analytical treatment. One of them is constituted by materials, which exhibit power-law grading of the elastic properties in the normal direction, i.e., into the material.
While axisymmetric contacts of power-law graded elastic materials have been studied in depth, and the corresponding state of the start is basically on par with the one for homogeneous media-including general solutions for non-adhesive normal contacts [12,13], with some constants E 0 and z 0 (mathematicians might treat E 0 /(z 0 ) k as one constant, but in the above formulation E 0 and z 0 have immediate physical meaning). Poisson's ratio ν shall be constant. Effects of microscopic surface roughness, inelasticity, or surface tension are neglected. The gap between the contacting bodies at the moment of first contact constitutes the contact profile f (x,y), with the cartesian coordinates x and y in the contact plane. The generalization for the contact of two elastic bodes is straightforward if both materials have the same exponent k of elastic grading. The case of different values for k is mathematically very complicated, although, very recently also for this problem, analytic contact solutions have been provided [26].

Fundamental Solution for Point Loading
If a normal point load F acts on the surface of an elastic half-space with elastic grading in the form (1), the half-space surface will experience elastic normal displacements [27]: with the Gamma function Γ and the polar radius r from the force application point. Hence, for a pressure distribution p(x,y), the resulting surface normal displacement is provided by the convolution with the contact domain Ω.

Generalization of Galin's Theorem
Rostovtsev [20,21] has generalized Galin's theorem for the homogeneous elastic halfspace, as he proved that a pressure distribution of the form with some polynomial S m (x,y) of the order m, corresponds to a polynomial normal surface displacement of the same order in the contact domain, and vice versa. Notably, Galin's theorem describes a logical equivalence, i.e., it is valid in both directions. Rostovtsev proved both directions of the generalization separately in the publications [20,21], respectively.

Generalization of Barber's Extremal Normal Force Principle
Another important theorem, which will be used in the following, is Barber's [28] principle for the homogeneous elastic half-space, that the contact domain in the corresponding indentation problem (if the domain is not fixed) can be determined from the condition with the indentation depth d, and the normal force integral [29] (p. 52) with the pressure distribution p* due to the unit indentation of the elastic half-space by a flat punch with the planform Ω.
Barber's proof is based on a harmonic function representation of the elastic displacements. Let us briefly consider an alternative proof, which is far easier to generalize mathematically for an inhomogeneous half-space. Equation (6) can be written in the form with the contact stiffness K of the flat punch contact and a functional G of the indenter profile, As was first pointed out by Mossakovskii [30] in the context of axisymmetric normal contacts without a slip of homogeneous elastic materials, the incremental difference between two contact configurations with the contact domains Ω and Ω + dΩ (dΩ is the change in the contact domain due to the change in indentation depth from d to d + dd) can be thought of as an infinitesimal indentation dd by a rigid flat punch with the planform Ω. Therefore, the (incremental) contact stiffness K(Ω) applies to any (frictionless) normal contact with the same contact domain, and we have from Equation (7) dF from which Equation (5) can be immediately deduced. Moreover, the flat punch superposition idea is in no way bound to the homogeneity of the elastic half-space. Hence, Barber's principle can be also applied straightforwardly for the power-law graded half-space.

Elliptical Flat Punch
The pressure distribution under an elliptical rigid flat punch with the half-axes a and b is an indispensable basis for the analysis in the following sections. It follows immediately as a special case of Rostovtsev's (or Galin's) theorem, because, according to Equation (4), a pressure distribution will produce a constant surface normal displacement inside the elliptical contact domain. For convenience, we may rewrite the pressure distribution in polar coordinates {r,ϕ} as with the contour of the contact domain (without loss of generality, we put a > b) The total normal force is defined by

Contact with a Power-Law Ellipsoidal Indenter
An interesting contact problem, which always exhibits an elliptical contact domain because of Galin's theorem, is the indentation by a power-law ellipsoidal indenter, with a constant eccentricity of the elliptical horizontal indenter cross sections. The contact profile reads with the indenter eccentricity ε and some constant B n . A very common special case of such a profile is the elliptical Hertzian contact with n = 1. For the contact with a homogeneous elastic half-space, very recently a general exact solution has been provided by the author for the Boussinesq problem with this class of counter bodies [31]-note that the case n = 2 was already considered for homogeneous materials by Shtayerman [32] (pp. 257 ff). Based on the fundamental relations laid out in the previous section, the general homogeneous solution can be generalized for the contact with a power-law graded elastic half-space without severe difficulties.
The following procedure is also applicable to more general elliptical contact profiles if the profile can be expanded as a series of the power-law functions (15). However, if the series is infinite, the procedure will only be approximate, because, generally, the contact domain will not be strictly elliptical in that case.

Macroscopic Contact Solution
The contact solution is based on the integral expression (6) for the force as a function of the indentation depth and the contact domain. For profiles in the form of (15), the contact domain is elliptical, and so p* and the contour of the contact domain are provided by Equations (11) and (12). Hence, the normal force integral reads which can be evaluated as Note that P n are polynomials of the order n in e 2 and ε 2 , which have no dependence on k.
In other words, they are already known from the corresponding homogeneous solution.
Let us characterize the contact domain by the smaller half-axis b and the eccentricity e. According to Barber's theorem, the values of b and e, which solve the non-adhesive normal contact problem, will maximize the force in Equation (17) for a given indentation depth. Maximizing with respect to b, we obtain Hence, and maximizing with respect to e leads to the condition (20) Equations (18)- (20) exactly solve the normal contact problem. For e = ε = 0, it is easily verified that the known axisymmetric results [7] (p. 265) are exactly recovered. Note that the trivial solution of Equation (20), e = 0, corresponds to a minimum of the force. For k = 0, the results simplify to the known homogeneous solution.

Pressure Distribution
For the determination of the pressure distribution, we return to the idea of understanding the indentation process as a series of incremental flat punch indentations.
As Equation (20) has no dependence on the normal load (i.e., either the indentation depth or the normal force), the contact eccentricity will remain constant during the indentation process. Moreover, the pressure distribution at the end of the indentation procedure directly follows from the integral of pressure distributions for the incremental flat punch indentations. Hence, in polar coordinates, using Equations (11), (12), and (18), we obtain It is interesting that, except for the constant factor before the integral, this corresponds to the known axisymmetric solution [7] (p. 265) scaled to the elliptical contact domain.

Example: Elliptical Hertzian Contact
Based on the general solution for arbitrary (but natural) n, let us recover the known solution for the Hertzian case with n = 1. We have P 1 (e 2 , ε 2 ) = (2 − e 2 − ε 2 )/2 and, therefore, The contact eccentricity is defined as the nontrivial solution of which can be shown to coincide with the known result [23].
In the limiting case k → 1, the solution is always e = ε. This is to be expected, as the Gibson solid (with k = 1 and ν = 0.5) reacts like a three-dimensional foundation of independent linear springs and, therefore, the contact domain always corresponds to horizontal indenter cross sections. For small eccentricities, Equation (23) has the asymptotic solution In Figure 1, the relation between the contact eccentricity e and the profile eccentricity ε is shown for three values of k, based on the exact solution (23). The thin drawn-through lines correspond to the respective asymptotic solution (24), which apparently is valid up to ε ≈ 0.2. Obviously, for increasing values of k, the contact eccentricity becomes smaller, i.e., closer to the profile eccentricity. This is in line with the above comment that, for k → 1, the solution is always e = ε. Finally, the pressure distribution equals 4 . 3 e k ε ≈ + (24) In Figure 1, the relation between the contact eccentricity e and the profile eccentricity ε is shown for three values of k, based on the exact solution (23). The thin drawn-through lines correspond to the respective asymptotic solution (24), which apparently is valid up to ε ≈ 0.2.
Obviously, for increasing values of k, the contact eccentricity becomes smaller, i.e., closer to the profile eccentricity. This is in line with the above comment that, for k → 1, the solution is always e = ε. Finally, the pressure distribution equals

A Fabrikant-Type Approximation for General Contact Domains
To constructively apply the extremal principle (5) to the force integral (6), one requires an expression for the pressure distribution p* under a flat punch with arbitrary cross section. In that regard, the form (11) for the pressure distribution under a rigid flat elliptical punch suggests that Fabrikant's approximation [24] for the pressure distribution under a flat punch of arbitrary planform, indenting a homogeneous elastic half-space, can be generalized straightforwardly to account for the type of functional elastic inhomogeneity considered in the present manuscript.
Noteworthily-as was pointed out very recently in [33]-what generally is referred to as "Fabrikant's approximation" is, in fact, a special case of an asymptotic expansion, which was developed earlier by the group of Mossakovskii [34]. Nevertheless, to avoid confusion, we will keep to the usual term.
Let us suppose a flat punch of arbitrary shape with the contour R(φ) in polar coordinates is pressed into a power-law graded elastic half-space. We postulate that the pressure

A Fabrikant-Type Approximation for General Contact Domains
To constructively apply the extremal principle (5) to the force integral (6), one requires an expression for the pressure distribution p* under a flat punch with arbitrary cross section. In that regard, the form (11) for the pressure distribution under a rigid flat elliptical punch suggests that Fabrikant's approximation [24] for the pressure distribution under a flat punch of arbitrary planform, indenting a homogeneous elastic half-space, can be generalized straightforwardly to account for the type of functional elastic inhomogeneity considered in the present manuscript.
Noteworthily-as was pointed out very recently in [33]-what generally is referred to as "Fabrikant's approximation" is, in fact, a special case of an asymptotic expansion, which was developed earlier by the group of Mossakovskii [34]. Nevertheless, to avoid confusion, we will keep to the usual term.
Let us suppose a flat punch of arbitrary shape with the contour R(ϕ) in polar coordinates is pressed into a power-law graded elastic half-space. We postulate that the pressure distribution will be defined approximately by the expression (11) (which is, of course, exact for the elliptical flat punch). Accordingly, the total normal force will be and, similarly, for the indentation depth, Hence, the contact stiffness approximately equals where the brackets denote averaging over the polar angle.
Example: Flat Punch with Square Planform As an illustrative example, let us consider a square with side length 2a as the face of the flat punch. For k = 0, we obtain for the contact stiffness with the incomplete Beta function B. Hence, if we write K = ca 1+k /(1 + k)/S, we obtain for the shape factor S It is interesting to compare this value with the result of Boyer's [35] approximate solution procedure for the same problem. Partitioning the square into four sub-squares with side length a, and calculating the displacement of one sub-square under the assumption that the loads on the other three sub-squares can be treated as point loads, applied at the sub-square center, one obtains with the fundamental solution (2) and above definition for the shape factor that It should be noted that the term "shape factor" in this framework can be slightly misleading, as the power-law graded elastic half-space has an intrinsic length scale z 0 , and therefore any contact problem of power-law graded materials also has size effects; this is blurred in the above definition (which is chosen such that the shape factor for the cylindrical flat punch always equals 1/2).
In Figure 2, the shape factors according to the approximations (30) and (31) are shown in comparison with a rigorous numerical solution obtained with the boundary element method (BEM) [36]. While Fabrikant's approximation for all k agrees very closely with the numerical solution, Boyer's approximation only in the homogeneous case provides satisfactory results.
Nevertheless, the "corner singularity" of the pressure distribution at the corners of the square is, in fact, more severe than predicted by Fabrikant's approximation, which only obtains an "edge stress singularity" for every point of the contact boundary. In other words, it treats all boundary points as parts of an edge, which is smooth along the contact contour. This effect is demonstrated in Figure 3. There, the BEM results for the pressure distributions along the x-axis and along the diagonal of the contact square are shown in normalized variables for a graded elastic half-space with k = 0.5 and compared with the respective results from the generalization of Fabrikant's approximation. While very good agreement between the approximation and the rigorous numerical solution is achieved along the x-axis, the approximation underestimates the pressure distribution along the diagonal and, especially, the stress singularity at the corner. Thus, a more detailed asymptotic analysis of Nevertheless, the "corner singularity" of the pressure distribution at the corners o the square is, in fact, more severe than predicted by Fabrikant's approximation, which only obtains an "edge stress singularity" for every point of the contact boundary. In othe words, it treats all boundary points as parts of an edge, which is smooth along the contac contour. This effect is demonstrated in Figure 3. There, the BEM results for the pressur distributions along the x-axis and along the diagonal of the contact square are shown in normalized variables for a graded elastic half-space with k = 0.5 and compared with th respective results from the generalization of Fabrikant's approximation. While very good agreement between the approximation and the rigorous numerical solution is achieved along the x-axis, the approximation underestimates the pressure distribution along th diagonal and, especially, the stress singularity at the corner. Thus, a more detailed asymp totic analysis of the corner singularity behavior will probably result in a more refined ap proximation for the pressure distribution.  Nevertheless, the "corner singularity" of the pressure distribution at the corners o the square is, in fact, more severe than predicted by Fabrikant's approximation, whic only obtains an "edge stress singularity" for every point of the contact boundary. In othe words, it treats all boundary points as parts of an edge, which is smooth along the contac contour. This effect is demonstrated in Figure 3. There, the BEM results for the pressur distributions along the x-axis and along the diagonal of the contact square are shown i normalized variables for a graded elastic half-space with k = 0.5 and compared with th respective results from the generalization of Fabrikant's approximation. While very goo agreement between the approximation and the rigorous numerical solution is achieve along the x-axis, the approximation underestimates the pressure distribution along th diagonal and, especially, the stress singularity at the corner. Thus, a more detailed asymp totic analysis of the corner singularity behavior will probably result in a more refined ap proximation for the pressure distribution.

General Approximate Solution for Almost Axisymmetric Contacts
Nevertheless, not only contact problems of flat punches with arbitrary contact domains can be solved approximately for homogeneous materials, but also general smooth profiles have been considered. Recently, based on Fabrikant's approximation and Barber's extremal principle, Popov [25] presented an analytic (first order) approximate contact solution for the indentation of a homogeneous elastic half-space by a rigid indenter, whose shape slightly deviates from axis-symmetry but is otherwise arbitrary. As now all constituent parts of that solution have provided for the case of the power-law graded elastic half-space, we can easily generalize Popov's solution for the elastically inhomogeneous problem.

General Approximate Solution Procedure for Non-Symmetric Profiles
The solution is based on a general approximate procedure, which, for the homogeneous case, was first suggested by Barber and Billings [37] and can straightforwardly be generalized for the power-law graded problem.
Let us suppose a rigid indenter with the profile f = f (r,ϕ) is pressed into the powerlaw graded elastic half-space. There shall be no tilting moments around the origin of the coordinate system, defined such that f (r = 0) = 0. The normal force is defined by Equation (6). If, for the pressure distribution p*, we use the generalization of Fabrikant's approximation, discussed in the previous section, the force integral approximately simplifies to with the generalized Abel transform According to Barber's principle, applied to the power-law graded case, the correct contour R(ϕ) of the contact domain for a provided indentation depth results from the solution of the variational problem (5). Let us first find the correct contour for fixed values of the non-linear moment J k . According to Lagrange, this is carried out by finding the unconditional extremum of the functional with the Lagrange multiplier λ. Applying the extremal principle results in In general, this will be a transcendent algebraic equation for R(ϕ), defying a closedform solution. However, in analogy to the homogeneous case, a closed-form asymptotic solution can be derived if the profile is only slightly deviating from axial symmetry.

Closed-Form Approximate Solution for Almost Axisymmetric Profiles
Considering the almost axially symmetric problem, we write f (r, ϕ) = f 0 (r) + f 1 (r, ϕ), ∀{r, ϕ} : and also introduce the corresponding generalized Abel transforms 1−k , g 0 (R(ϕ)) = dG 0 (R(ϕ)) dR(ϕ) , The contact domain will be almost circular with the contour Substituting Equation (38) into Equation (35), expanding in powers of R 1 /R 0 , and neglecting all terms of higher than linear order, we obtain where the prime denotes the derivative with respect to the first functional argument. The determination of the Lagrange multiplier is completely analogous to the homogeneous derivation presented by Popov [25] and shall therefore not be detailed here for reasons of space. We obtain, after some elementary calculations, Hence, the force integral (32) further simplifies to Here, the second and third integrals vanish in perfect analogy to the respective homogeneous solution (see [25] for the respective details). Finding the extremum of the remaining trivial first integral with respect to R 0 results in the axisymmetric solutions and Finally, with Equation (42), Equation (40) can be further simplified to provide Equations (42)-(44) define the closed-form analytic, approximate solution to the almost axially symmetric Boussinesq problem. Of course, for k = 0, the known homogeneous approximate solution is exactly recovered. Note that either d, F, or R 0 can be used to completely characterize the contact configuration.

Pressure Distribution for Almost Axisymmetric Profiles
As was pointed out before, the indentation process from the indentation depth d* = 0 to d* = d can be understood as a series of incremental flat punch indentations dd*, with the planform of the punch coinciding with the contact domain at the current indentation depth d*. Once the macroscopic contact problem has been solved, i.e., the contact domain as a function of the indentation depth has been determined, the resulting pressure distribution can be obtained as the integral of the pressure distributions from the corresponding series of incremental flat punch indentations. In an approximate sense, this procedure can always be executed analytically based on the approximation for the pressure distribution under a flat punch of arbitrary shape. Hence, we obtain for the pressure distribution where R 0 must be understood as a function of R(ϕ), by inverting Equation (44).

Example I: Elliptical Hertzian Contact
As an illustrative example, let us consider the elliptical Hertzian contact with The separation into axisymmetric and non-axisymmetric parts results in Hence, and Equations (42)-(44) provide the approximate contact solution For small eccentricities, the last relation describes an ellipse with half-axes and the eccentricity of the contact domain is approximately provided by which is in perfect agreement with Equation (24). It is also checked easily that, for small eccentricities, the obtained solutions for d and F agree with the asymptotes of the respective exact solutions in Equation (22).

Example II: Indentation by a Pyramid with Square Planform
As a second example for the approximate procedure, let us consider the normal contact with a rigid pyramid of square planform. The application of the analytic approximate solutions (42)-(44) for this profile geometry is straightforward and shall not be detailed here for reasons of space. In Figure 4, a comparison is shown between a BEM-based numerical solution and the approximate solution for the contour of the contact domain in normalized variables. As expected, the analytic approximate solution works very well, except for the edges of the pyramid, for which the contact domain is "smoothed" in the approximate solution.

Example II: Indentation by a Pyramid with Square Planform
As a second example for the approximate procedure, let us consider the normal contact with a rigid pyramid of square planform. The application of the analytic approximate solutions (42)-(44) for this profile geometry is straightforward and shall not be detailed here for reasons of space. In Figure 4, a comparison is shown between a BEM-based numerical solution and the approximate solution for the contour of the contact domain in normalized variables. As expected, the analytic approximate solution works very well, except for the edges of the pyramid, for which the contact domain is "smoothed" in the approximate solution.

Discussion and Conclusions
Several non-adhesive, frictionless indentation problems for the power-law graded elastic half-space have been considered in the case of elliptical and non-standard contact domains. Different powerful techniques, which were originally developed for the closedform analysis of the Boussinesq problem for a homogenous half-space, and which already proved their usefulness, were generalized for the elastically inhomogeneous problem. Nonetheless, some restrictions have to be kept in mind when applying the obtained results to real engineering contacts.
These restrictions can be classified into two groups. The first stems from the physical modelling and mainly comprises the assumptions of linearly elastic material behavior and the very specific form of functional grading (i.e., the power law). Real contacts will usually exhibit different forms of inelasticity; moreover, the precise form of the functional elastic

Discussion and Conclusions
Several non-adhesive, frictionless indentation problems for the power-law graded elastic half-space have been considered in the case of elliptical and non-standard contact domains. Different powerful techniques, which were originally developed for the closedform analysis of the Boussinesq problem for a homogenous half-space, and which already proved their usefulness, were generalized for the elastically inhomogeneous problem. Nonetheless, some restrictions have to be kept in mind when applying the obtained results to real engineering contacts.
These restrictions can be classified into two groups. The first stems from the physical modelling and mainly comprises the assumptions of linearly elastic material behavior and the very specific form of functional grading (i.e., the power law). Real contacts will usually exhibit different forms of inelasticity; moreover, the precise form of the functional elastic grading might deviate from the power-law description, especially if it originates from a hard-to-control technical grading procedure. However, the power law is a very convenient way of approximating more complicatedly graded properties, and, due to the applicability of all results for positive and negative values of the grading exponent, allows for the analysis of both hard surfaces with a softer material core as well as soft surfaces with a harder material core.
The second group of restrictions stems only from the mathematical analysis. As the Boussinesq problem for a general contact domain, despite its linearity, cannot be solved analytically in closed form, several approximate procedures were facilitated to allow for a fast, analytic solution of the contact problem. Although these procedures are asymptotically exact, i.e., they provide exact first-order solutions of a general problem around a specific expansion "point" (usually, the axisymmetric problem), their applicability should, if possible, be checked by more rigorous numerical simulations, e.g., based on the BEM or the finite element method (FEM), at least for some random problem samples.
On the other hand, the obtained normal contact solutions can be used to solve other classes of contact problems that can, exactly or approximately, be reduced to the elastic Boussinesq problem. This is true for viscoelastic normal contacts, via the elastic-viscoelastic correspondence principle [38], or tangential contacts with friction, via the generalization for power-law graded materials [19] of the principle by Jäger [39] and Ciavarella [40].
Finally, the approximate solution presented in Section 5 can be considered an extension of the method of dimensionality reduction (MDR) for power-law graded elastic materials [18] to non-axisymmetric contacts. It can be used for different problems in tribology and engineering, e.g., within the framework of indentation testing or as a tribological tool for the analysis of local features of a random rough surface.