Characterization of the Functionally Graded Shear Modulus of a Half-Space

In this article, a method is proposed for determining parameters of the exponentialy varying shear modulus of a functionally graded half-space. The method is based on the analytical solution of the problem of pure shear of an elastic functionally graded half-space by a strip punch. The half-space has the depth-wise exponential variation of its shear modulus, whose parameters are to be determined. The problem is reduced to an integral equation that is then solved by asymptotic methods. The analytical relations for contact stress under the punch, displacement of the free surface outside the contact area and other characteristics of the problem are studied with respect to the shear modulus parameters. The parameters of the functionally graded half-space shear modulus are determined (a) from the coincidence of theoretical and experimental values of contact stresses under the punch and from the coincidence of forces acting on the punch, or (b) from the coincidence of theoretical and experimental values of displacement of the free surface of the half-space outside the contact and coincidence of forces acting on the punch, or (c) from other conditions. The transcendental equations for determination of the shear modulus parameters in cases (a) and (b) are given. By adjusting the parameters of the shear modulus variation, the regions of “approximate-homogeneous” state in the functionally graded half-space are developed.


Introduction
Spatial inhomogeneity of material properties widely occurs in natural and engineered materials. Research of inhomogeneous materials was initially stimulated by civil engineering due to the necessity in calculation of foundations resting on inhomogeneous soils. Also, spatial variations of material properties were observed in concrete products such as panels, foundation blocks, various slabs, etc. Increased demand in strength, thermal stability, weight reduction of products, their resistance to chemical influences, etc. led to intended fabrication of inhomogeneous materials. Initially it resulted in adoption of multilayered composites, and then to creation of functionally graded materials. Nowadays, functionally graded materials (FGMs) are used in various branches of modern industry. One of notable application areas for FGMs is semiconductor microelectronics [1,2]. Many FGMs are based on traditional materials such as aluminium, zirconium, etc. Those materials, coated with appropriate oxides of technological or natural origin, are used to manufacture semiconductor products for microelectronics. Elastic moduli of oxides, nitrides, etc. are significantly different from those of the base material, and the thickness of the coating sometimes comprises a substantial part of the total product thickness. For correct analysis of the mechanical, electrical and magnetostatic state of such products, it is necessary to take into account changes in elastic moduli, and electrical and magnetic resistance along the length or depth of the product material.
Recently, much attention has been paid to the creation of calculation schemes to determine the stress-strain state of FGMs having Young's modulus or shear modulus varying by the depth of the material. Knowing the law of elastic moduli variation within the material enables one to obtain closed-form solutions for some particular contact conditions. References to early works on this subject can be found in [3]. An exponential law of variation is most often considered due to the simplicity of analysis but also due to its wide occurrence in materials as a result of natural and technological processes. Selvadurai et al. [4] considered torsion of an exponentially graded half-space. Giannakopoulos and Suresh [5,6] considered action of a point force on a half-space with power law and exponential law inhomogeneity and indentation of the half-space by an axisymmetric indenter. The authors obtained some closed-form solutions and compared them with results of finite element analysis (FEA). Guler and Erdogan [7] analysed plane indentation of a functionally graded layer on a homogeneous substrate in the presence of friction. Mindlin's problem for an incompressible elastic half-space with an exponential variation in the linear elastic shear modulus was considered by Selvadurai and Katebi in [8]. Tokovyy and Ma [9] analysed the elastic equilibrium of an exponentially inhomogeneous layer and constructed a numerical procedure for calculation of stresses in the layer in the case of an arbitrary specified loading on the surface. Selvadurai and Katebi [10] considered joint action of a rigid circular punch and a concentrated force on a half-space with exponentially varying shear modulus. Other forms of elastic moduli variation were also considered: Altenbach and Eremeyev [11] analysed eigen-vibrations of a functionally graded material with a power law variation that corresponds to a graded distribution of the porosity in the material. Tangential loading of a power-law inhomogeneous layer bonded to a homogeneous half-space was considered by Kulchytsky-Zhyhailo and Bajkowskiy [12]. Awojobi [13] considered a hyperbolic variation of the elastic modulus, which modelled smooth transition of the sub-soil into a rigid bed. Chen and Chen [14] analysed contact between a homogeneous half-space with a linearly graded layer and a rigid punch. There are also several approaches to construct effective solutions for an arbitrary variation of elastic properties by depth of the material [15][16][17][18][19][20][21].
Nanoindentation techniques are widely used for determining Young's modulus of homogeneous materials and parameters of elastic moduli of inhomogeneous materials. Simulation of a nanoindentation procedure is reduced to solving contact problems about the indentation of functionally graded materials by an indenter (punch). Numerical approaches are being developed for determination of Young's modulus variation by depth [22][23][24][25][26]. Nakamura et al. [22] proposed the inverse analysis of spherical indentation of an FGM layer on a thick substrate. This analysis is used to evaluate two parameters: the exponent n in the power law compositional variation and the stress-strain transfer variable q. From these two parameters, one can determine Young's modulus and yield strength distribution through the layer thickness. Poisson's ratio is assumed to be constant across the FGM layer. The analysis is based on matching a finite element model of elastoplastic indentation with experimental load-displacement data using the Kalman filter technique. Practical application of this approach was illustrated in [23]. Bocciarelli et al. [24] used a similar finite element elastoplastic model of conical indentation of a power-law FGM layer and provided their own inverse analysis technique for identification of n and q. They tested the analysis on theoretical data and showed that most of the parameters can be determined, including those of material constituents. At the same time, the analysis showed poorer quality in determination of elastic moduli of metal and ceramic phases in the FGM layer. Huang et al. [25] developed a finite element model of an FGM plate subjected to tension load. It was assumed that Young's modulus and Poisson's ratio vary by linear or exponential law in direction orthogonal to the applied load. Based on theoretical data, they demonstrated the ability of the inverse analysis to determine Young's modulus of one of the constituent materials if the second one is known and Poisson's ratio is constant across the plate. Recently, a meshfree method was proposed by Chen et al. [26] for characterization of elastic moduli of a linear elastic isotropic FGM with volume fraction governed by power or exponential law. The methods of fundamental solutions and radial basis function approximation were used to solve a direct problem on uniform compression of an FGM beam. The authors proposed to fit the model to the experimental displacements data by the Neider-Mead simplex (NMS) approach. Using the theoretical data, they showed the possibility to determine Young's moduli of constituents and exponent factor in power or exponential law. Again, Poisson's ratio was assumed to be constant in the considered examples.
To fully characterize a linearly elastic isotropic material, the second elastic modulus is needed, which can be identified from analysis of the experiment on shear or torsion of a functionally graded material. An experiment using inverse analysis of torsion for characterization of a graded shear modulus of a coating was proposed in [27]. There are also some works dealing with the shear problem: Aizikovich [28] considered shear of an exponentially inhomogeneous linearly elastic half-space; anti-plane shear deformations for homogeneous and arbitrarily depth-wise inhomogeneous anisotropic linearly elastic solids were considered by Horgan and Miller [29]. Anti-plane surface waves in a functionally graded media with exponential [30,31] or arbitrary [32,33] dependence of shear modulus and density on depth are also objects of ongoing research. Particular attention is paid to shear of polymeric damping materials. Pritz [34] suggested a five-parameter fractional derivative model to describe viscoelastic materials behaviour and illustrated dynamic shear properties of some materials used in practice. Alaimo et al. [35] applied the fractional derivative model to analysis of damped free-vibration of multilayered composite plates. Unfortunately, these works do not reference approaches to identification of inhomogeneous material parameters.
In the present work we develop the theoretical and experimental methodology for determining the parameters of the functionally graded half-space shear modulus varying by depth according to the exponential law. At the initial stage the auxiliary problem about the anti-plane shear of the surface of a functionally graded half-space by the base of a strip punch is formulated and solved. The solution of the auxiliary problem is reduced to an integral equation of the first kind with a difference kernel. The analytical solutions of the integral equation are constructed by asymptotic methods with respect to the dimensionless parameter of the problem, as was proposed in [36]. The dependence of basic mechanical characteristics of the problem on two parameters of exponential shear modulus of the half-space is studied. Earlier, the method was proposed to determine exponential shear modulus parameters from values of contact stresses and by using a static condition in [37]. However, this setup is problematic to implement due to low accuracy of stress measurements. Because of this, here we propose an alternative scheme for determination of the shear modulus parameters from the equality of experimental and theoretical data on surface displacements outside the contact area, in combination with the static condition in the contact region. Also, other possible conditions are suggested for determination of the shear modulus parameters in this paper. Differences between the proposed approach and the numerical schemes discussed above consist in using the analytical solution to the contact problem, from which the two parameters of the functionally graded half-space shear modulus are unambiguously determined.

Formulation of the Research Problem
Let us consider the elastic half-space made of the functionally graded material having the shear modulus changing with depth by the exponential law where the parameter µ 0 = µ(0) represents the value of the shear modulus at the surface of the half-space, and the parameter d defines the intensity of change of the shear modulus with depth in the half-space. The aim is to determine the parameters µ 0 and d of the shear modulus that vary by depth according to Equation (1).
To achieve this goal, we propose the experimental approach of anti-plane shear of the sample of a functionally graded material by a strip punch. In the process of implementing such an experiment, some of main parameters of contact interaction need to be measured: the value of shear (displacement) of the punch, the contact force, the contact stresses under the punch, the value of displacements of the free surface of the sample outside the contact area, the slope of the displacement curve outside the contact area, etc. To use the data obtained from the experiment to determine the parameters µ 0 and d, a mathematical simulation of the shear of a functionally graded material by a strip punch should be conducted. It is reduced to the solution of the contact problem of shear of a functionally graded half-space by a strip punch, when shear modulus µ(y) of the half-space changes exponentially with y according to Equation (1). The analytical solution to this problem will provide relations for determining the abovementioned characteristics: punch displacement, contact force, contact stress, free surface displacement outside the punch, slope of the displacement curve, etc. The comparison of experimental and theoretical results will make it possible to determine parameters µ 0 and d of the shear modulus µ(y).

Mathematical Formulation of the Problem
We consider the contact problem of elasticity theory about the shear of the elastic half-space surface (y = 0) by a strip punch of width 2a (|x| ≤ a, y = 0) with magnitude ε (w(x, 0) = ε, |x| ≤ a), where w(x, y) is the shear displacement in the half-space 0 ≤ y < ∞ along the axis z. The scheme of the contact problem is presented on Figure 1.
The shear modulus of the functionally graded material in the half-space depends on the depth y according to Equation (1). The contact stresses ϕ(x) = −σ yz (x, 0) arising under the punch and other mechanical characteristics of the problem are to be determined. It is assumed that at infinity in half-space the displacements w(x, y) and stresses σ xz (x, y), σ yz (x, y) vanish as The exponential law of elastic modulus variation is chosen because of its wide occurrence in engineered systems and nature. Composition of bimaterial graded coatings often vary exponentially [38]. When studying the Oxford geological sediments [8], it was found that the shear modulus µ(y) changes with depth by the exponential law with parameters µ 0 = 3.33 MPa, 2d = 0.0879 m −1 .
In the half-space under the pure shear deformation regime the conditions σ xx = σ xy = σ yx = σ yy = σ zz = 0 are fulfilled, where σ xx , σ yy and σ zz are normal pressure components, and σ yx , σ xy are tangential pressure components. Then the equations of the theory of elasticity in stresses [39] in the case of anti-plane deformation are as follows where σ zx and σ zy , in case of functional inhomogeneity of elastic medium, are expressed through strains by formulae By substituting Equation (3) into differential Equation (2) and then dividing it by µ(y) 0, we obtain a differential equation in displacements The differential Equation (4), with a factor µ (y)/µ(y) dependent ony, in case of µ(y) = µ 0 e 2dy is transformed into a differential equation with constant coefficients Mixed boundary conditions of the formulated contact problem about the pure shear of surface of the functionally graded half-space by a strip punch at y = 0 have the form where ϕ(x) are the contact stresses under the punch to be determined.
At infinity x 2 + y 2 → ∞ displacement w and its derivatives ∂w ∂x , ∂w ∂y satisfy lim √

Integral Equation of the Problem
To construct the solution to the differential Equation (5) the Fourier integral transform is used: where the upper index F means the Fourier transform of the function. By means of the Fourier integral transform shown in Equation (9) and conditions at infinity shown in Equation (8) the differential Equation (5) is reduced to the solution of an ordinary differential equation with constant coefficients relative to the transform w F (α, y): The general solution to this equation, which vanishes at infinity ( y → ∞ ), has the form where C(α) is determined from the boundary condition of Equation (7) after application of the Fourier integral transform shown in Equation (9) and has the form The displacements w(x, y) are determined from Equation (11) by means of the inverse Fourier transform of Equation (9) and have the form where Satisfying the mixed boundary conditions of Equations (6) and (7) and using Equation (12) we obtain the integral equation of the problem where λ − is from Equation (11).
It should be noted that at d = 0 we obtain a classical integral equation for a homogeneous half-space with µ(y) = µ 0 = const.
By substituting α = d α , dα = d dα into the internal integral and making internal ξ = a ξ , dξ = a dξ and external x = a x variables dimensionless, the integral Equations (13) and (14) are reduced to the integral equation in dimensionless form (primes omitted) which is an integral equation of the first kind of the Fourier convolution type with a difference kernel with respect to unknown dimensionless contact stresses ϕ(x) = aε −1 µ −1 0 ϕ(x). The transform K(α) of the kernel k(t) of Equation (16) in the integral Equation (15) is an even function, multi-valued in the complex plane α with algebraic branch points α = ±i and expresses the following asymptotic properties: where K(0) = 0.5. After calculating the integral in Equation (16), the kernel k(t) looks like where K 0 (t) is the McDonald's function. At small t for k(t) the following behaviour takes place: The properties listed for k(t) allow us to conclude [40] that the solution to the integral Equation (15) exists and is unique in the class of functions ϕ( (15) and (16), we use asymptotic methods with respect to the dimensionless parameter λ = (da) −1 .

Effective Solution for Small λ
At small values of a dimensionless parameter λ the solution to the integral Equations (15) and (16) are constructed as a Neumann series, the zeroth term of which is set by the ratio [40] ϕ Functions ϕ ± (x) are determined from the Wiener-Hopf integral equation on a semi-axis [40,41], corresponding to Equation (15 The function ϕ ∞ (x) is determined from the integral Equation of the convolution type on an infinite interval, corresponding to Equation (15 The solution to the integral Equation (23) is constructed using the generalized Fourier integral transform [42] and has the form The solution to the integral Equation (22) on a semi-axis is based on the Wiener-Hopf method [43]. In the first step of the method, the integral equation is extended to the whole axis. Thus, for ϕ + (x) we have where v − (x) = ε −1 w − (x) are the dimensionless displacements of the free surface of the half space outside the contact area (to the left of the punch), which are given by the operator and w − (x) are the dimensional horizontal displacements of the free surface of the half space outside the contact area (to the left of the punch).
In the second step of the Wiener-Hopf method, the Fourier integral transform of Equation (9) is used to solve the integral Equation (25), by means of which the integral Equation (25) is reduced to a functional equation with respect to transforms ϕ F in which the unknown transforms are regular in the upper (Im(α) > η − ) and lower (Im(α) < η + ) half-planes of the complex plane α, respectively. Assuming that the factorization K(α) is known, i.e.,K(α) is represented as where K + (α) is regular in the upper (Im(α) > η − ), and K − (α) is regular in the lower (Im(α) < η + ) half-plane, and after substitution of Equation (27) in Equation (26), the left and right parts of the resulting equation are divided on K − (α). Let us introduce the representation where f + (α) is regular in the upper half-plane (Im(α) > η − ) and f − (α) in the lower half-plane (Im(α) < η + ). Then Equation (26) can be written as which holds in the strip η − < Im(α) < η + , the left and right part of which are continuations of some complex variable function F(α) in the upper (Im(α) > η − ) and lower (Im(α) < η + ) half-planes, correspondingly. Since the left and right parts vanish at infinity ( |α| → ∞ ), the function F(α) = 0 in the complex plane α = ξ + iη. Then two equations follow from Equation (29): which are used to determine the unknown transforms. Inverting the Fourier transform in Equations (30) and (31) with the help of Equation (9) we obtain the solution to integral Equation (25): where η − < c < d < η + , and η ± is taken as ±1 based on the properties of K(α) from Equation (16). The integral representations for ϕ − (x) and v + (x) completely coincide with Equations (32) and (33), respectively, where v + (x) = ε −1 w + (x), v + (x) and w + (x) are the dimensionless and dimensional (correspondingly) displacements of the free surface of the half-space outside the contact area (to the right of the punch).
In order to calculate integrals in Equations (32) and (33) with K ± (α), it is necessary to perform the factorization of kernel transform K(α) of Equation (27). The exact factorization K(α) of Equation (16) was previously given in [43] but it is not very effective for obtaining a foreseeable solution. To get the effective solution, we use the K(α) approximation that can be factorized by elementary means.
The punch shear force T, which is an integral characteristic of the problem, is determined by the integration of the obtained solution over the interval [−a, a] where is a dimensionless shear force. Taking ϕ(x) in the form of a superposition Equation (21) and integrating over the interval [−1, 1], we obtain the dimensionless shear force T in the form The use of an asymptotic solution to the integral Equation (15) in multiplicative form [40] ϕ where ϕ ± (1 ± x)λ −1 from Equations (40)- (42), and ϕ ∞ (xλ −1 ) from Equation (24), allows for more convenient formulae for T to analyse all considered approximations of K(α) listed in Equations (34)-(36): where Equations (50)-(52) have the same order of accuracy for small λ as in Equation (46).

Effective Solution for Large λ
When constructing the solution to the integral equation of the problem listed in Equations (15) and (16) in the case of large values of λ, it is necessary to take into account the fact that K(α) can be represented by a series in which c k at k = 0, 1, 2 , . . . , N and for N = 4 have the following values: With the help of series representation by Equation (53) and the asymptotic property shown in Equation (17), it is not difficult to obtain the decomposition of k(t) in the form of the following series [45] 1 2 k(t) = − ln|t| + a 20 |t| + a 11 t 2 ln|t| + a 21 |t| 3 + a 31 t 2 + O t 4 ln|t| at t → 0, where a 20 = − π 2 c 1 , a 11 = 1 2 c 2 , a 21 = 1 12 πc 3 , and After substitution of Equation (54) into Equation (15) and the differentiation by an external variable x, a singular integral equation is obtained with respect to ϕ(x). Inverting the singular integral equation, an equivalent integral equation of the second kind is formed [45]. The solution to the latter is constructed by the successive approximations method in the form of a double functional series with the following expressions for functions ω mn (x): where b 12 = 8 9 a 11 a 20 , b 10 (x) = 6a 21 The dimensionless shear force (47) then takes the form

Analysis and Application of the Obtained Solutions
The obtained asymptotic solutions by the parameter λ = (da) −1 of the exponentially functionally graded material shear problem cover the entire range of the parameter λ ∈ (0, ∞). Table 2 shows the numerical values of the relative contact stresses aT −1 ϕ(0) calculated by Equations (21)  The analysis of values of the relative contact stresses shows that the solutions obtained for small and large values of λ coincide on the interval λ ∈ (1 , 2).

Determination of Shear Modulus Parameters
Previously, it was noted that researchers are increasingly focused on the development of methods for determining the parameters of elastic moduli of functionally graded materials [22][23][24][25][26][27], based on numerical solutions to contact problems on indentation. The analytical solutions of the problem obtained in the previous paragraphs allow us to analyse the correlation between mechanical parameters of the problem. Table 3 shows the values of the relative stresses aT −1 ϕ(0) in the centre of the punch, where x = 0 at a = 1 m for the different values of the parameter d. At small d, for example at d = 0.001 (λ = 1000), the value of the relative stresses aT −1 ϕ(0) coincides with the value of the relative contact stress aT 0 −1 ϕ 0 (0) at x = 0 in the classical solution to the homogeneous half-space shear problem [45] aT 0 1.53 times the value aT −1 ϕ(0) at d = 0.001 (λ = 1000), which also coincides with the classical solution.
It follows that the parameter d of the shear modulus µ(y) as in Equation (1) affects the contact stresses. The meaning of the parameter h in Table 3 is explained in Section 9. Table 3. Relative contact stresses aT −1 ϕ(0) for a = 1 and different values of shear modulus parameter d. The data in Table 3 can be used to solve the problem of determining the shear modulus parameter d by the values of the relative contact stresses aT −1 ϕ(0) in the centre of the contact area in the following way.
Firstly, the parameter d (0 < d < ∞) of the functionally graded half-space shear modulus is determined from the solution to the transcendental equation with respect to d (at a = 1 m) where ω mn (x) are from Equation (57), ϕ ± (x) are from Equations (40)- (42), and ϕ ∞ (x) is from Equation (24). In Equation (59) T is taken from Equation (58), and in Equation (60) T is taken from Equation (24). T and ϕ(0) are the force and value of contact stresses in the middle of the punch, determined experimentally or by some other means.
Secondly, to determine the shear µ 0 modulus parameter, which coincides with the value of the shear modulus on the half-space surface (at y = 0), the static condition of Equation (46) is used: where T is a shear force determined from experiment. The value T of Equation (47) is determined for 0 < d < 1 and 1 < d < ∞ from Equations (59) and (60), respectively. Table 4 shows the values of the relative stresses aT −1 ϕ(x) at points x = 0.00, 0.25, 0.50, 0.75 for the same values of parameters λ and d as in Table 3 (at a = 1m). The abovementioned example of determining the parameter d of the depth-varying shear modulus µ(y) as defined in Equation (1) by values of contact stresses aT −1 ϕ(x) is difficult to implement in practice due to the fact that contact stresses are not directly measured (or almost cannot be measured) and at best are determined indirectly or approximated. In the latter case, the approximate value of contact stresses ϕ(0) is often taken as the ratio of an experimentally determined shear force T to the contact length 2a, i.e., ϕ(0) = 0.5 Ta −1 . The indirect and approximated methods for determining contact stresses lead to additional errors in determining the shear modulus parameters. In this regard, we also propose another method for determination of the parameter d of shear modulus µ(y) using an experimentally and theoretically determined value of displacement w(x) of the free material surface outside the contact area. Figure 2 shows plots of dimensionless free surface displacement v(x) to the right of the punch for d = 5, 2, 1, 0.5, 0.2, 0.1 (at a = 1).  Table 5 shows the values of the dimensionless displacements v(x) (v(x) = ε −1 w(x)) outside the contact area for 1 < |x| < ∞ at x = 1.01, 1.1, 1.2, 1.3, 1.4, 1.5, calculated for the values λ and d in Table 5 (at a = 1). Displacements w(x) at large λ and 1 < |x| < ∞ are calculated by the formula where ϕ(x) is defined by the Equation (56), and k(t) by Equation (16). At small λ the displacements of the free surface w(x) are calculated using the equation where v ∓ (x) are defined by Equations (43)- (45); H(x) is the Heaviside function. The parameter d (0 < d < ∞) of the shear modulus µ(y) from Equation (1) is determined from the solution to the following transcendental equation for |x| ∈ (1, ∞) (at a = 1): The second parameter of the shear modulus µ 0 is determined, as in the previous case, from the static condition by Equation (61) through an experimentally determined punch shear force T.
Here we should also pay attention to the possibility of the parameter d estimation through the slope of the displacements curve w(x) outside the contact area (at 1 < |x| < ∞), which is numerically equal to the value of the derivative w (x) for 1 < |x| < ∞ and which can be determined by experimental means. To determine w (x), it is sufficient to differentiate the Equations (64) and (65), that is in which ϕ(x) is from Equation (56), and the function k * (t) is given by the equation and v ∓ (x) is defined in Equations (43)- (45). As before, the second parameter µ 0 of the shear modulus µ(y) from Equation (1) is determined by the Equation (61). For practical implementation of this method, one needs to measure displacements at several points of the free surface, for example, by using laser triangulation measurement [46] or interferometry [47].
Other additional conditions can be proposed from which the parameters of the exponential shear modulus as in Equation (1) of the functionally graded material can be determined.
Thus, the analytical solutions of the contact problem obtained here allowed us to effectively determine the parameters d and µ 0 of the shear modulus µ(y) of Equation (1) of a functionally graded half space.

"Approximate Homogeneous" Area
The fifth column of Table 3 indicates the depth coordinate h in the functionally graded half-space for which the values of the shear modulus µ(y) on the segment [0, h] and the values of the shear modulus µ 0 on the half-space surface differ by no more than 1%. This means that using the shear modulus parameterd in a half-space, with an error of 1%, one can create areas of "approximate-homogeneous" state with the same error. The same effect on the creation of regions of "approximate-homogeneous" state in a functionally graded half-space was achieved in [13] by using the parameter h of hyperbolic shear modulus µ(z) = µh(h − z) −1 .
When solving the problem about a shear of a homogeneous half-space by a strip punch it is impossible to establish an unambiguous connection between the displacementsw(x, y) in the half-space {|x| < ∞,0 < y < ∞} and contact stresses ϕ 0 (x) on {|x| < 1,y = 0 }, because the displacements w(x, y) in the half-space in this case are calculated by the equations where the integral in Equation (70) is divergent for any t and y. At the same time, the correctness of Equations (69) and (70) is also confirmed by the fact that they can be derived from Equation (12) applicable to the problem about the shear of a functionally graded half-space, which was considered here, if we set d = 0. It should be noted here that the Equations (69) and (70) point to another defect in the solution to the problem on shear of a homogeneous half-space. During the formulation of the homogeneous half-space shear problem, the boundary conditions on infinity ( x 2 + y 2 → ∞ ) are accepted, which consist of vanishing of displacements w(x, y) and stresses σ xz (x, y),σ yz (x, y) (i.e., partial derivatives of displacements ∂w ∂x , ∂w ∂y ). These conditions are of the same kind as Equation (8) specified in the formulation of the functionally graded half-space shear problem. Having obtained the solution to the problem on the shear of a homogeneous half-space in the form of Equations (69) and (70) it is clear that for this solution the condition shown in Equation (8) is not satisfied for 0 < y < ∞ because the integral in Equation (70) diverges.
From the obtained solution to the problem about a functionally graded half-space shear with the shear modulus governed by Equation (1) it follows that the relationship between w(x, y) and ϕ(x) is defined unambiguously by the operator where the integral in (72) converges at d > 0. It should be noted hence that the Equation (71) can be used as an approximate solution to the homogeneous half-space shear problem. To achieve this goal, one needs to select the shear modulus µ(y) parameterd in such a way that the half-space will be "approximate homogeneous" to a certain depth h (see Table 3) with a predetermined error. Stresses σ yz (x, y) in the functionally graded half-space are determined unambiguously from the contact stresses ϕ(x) using the following formula: σ yz (x, y) = − εdµ(y) π The obtained solutions for w(x, y) and σ yz (x, y) represented by Equations (71)-(74) allow analysis of their behaviour in the half-space |x|< ∞ , y > 0 and indicate their dependence on the shear modulus µ(y) and the depth y. Representing ζ(α) as ζ(α) = 2 + α 2 ζ −1 (α) in Equations (71)-(74) we obtain that the shear stress σ yz (x, y) and displacement w(x, y) in the functionally graded half-space decrease with depth y proportional to 1 √ y and 1 √ yµ(y) , respectively. This confirms that the obtained solution satisfies the condition at infinity in Equation (8) in the functionally graded half-space shear problem formulation.
In conclusion, it should be noted that the shear modulus µ(y) in the "approximate homogeneous" region of the half-space 0 ≤ y ≤ h is determined with the error of determining the "approximate homogeneity" region by Equation (61) taking into account that d << 1 and µ(y) = µ 0 = const for 0 ≤ y ≤ h.
As an option, the shear modulus of a homogeneous half-space can be determined through Young's modulus E and Poisson's coefficient ν by the method proposed in [48]. To do that one just needs to use the formula µ 0 = 0.5 E (1 + ν) −1 .

Conclusions
We proposed a scheme for determination of the parameters of the shear modulus exponentially varying by depth in the half-space. The scheme is based on a solution of the contact problem on shear of a functionally graded half-space by a strip punch. The obtained solution was used to analyse the effect of the parameters of the depth-wise exponential shear modulus of the elastic half-space on the main contact parameters: shear stresses, shear force acting on the punch, displacements of the surface of the half-space outside the contact area and other parameters. From this solution, we derived formulas or equations for determining the parameters of the exponential shear modulus of the functionally graded half-space from experimental values of the contact parameters.
At the same time, the obtained formulas can be used to select the parameters of the exponential shear modulus of the functionally graded half-space: • to determine the boundaries of the "approximately homogeneous" region of the functionally graded half-space with a predetermined error; • to establish an unambiguous relationship between the contact stresses and the displacements in the "approximately homogeneous" region of the functionally graded half-space. This cannot be performed in solution to the contact problem on shear of the homogeneous half-space. Thus, the methodology developed in this paper expands the range of solvable problems.

Conflicts of Interest:
The authors declare no conflict of interest.