Study of energies and oscillator strengths of Fe XXI including plasma shielding effects

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Contact mechanics of two elastic spheres reinforced by functionally graded materials (FGM) thin coatings X.W. Chen, Z.Q. Yue∗ Department of Civil Engineering, The University of Hong Kong, Hong Kong, PR China

a r t i c l e

i n f o

Keywords: Contact mechanics Coated sphere Functionally graded materials (FGM) Multi-layer half space Nonlinear Fredholm integral equation

a b s t r a c t This paper examines the normal frictionless point contact between two dissimilar elastic spheres reinforced by FGM coatings with arbitrarily varied shear modulus and Poisson’s ratio. Both the coating thicknesses and contact area are small compared to the radiuses of the two spheres. The material inhomogeneity of the FGM coating is only along the radial direction. The contact problem is formulated and reduced to a nonlinear Fredholm integral equation of the second kind. A multi-layer elastic half space model is used to develop the relationship between the contact stress and the vertical displacements of the FGM coated spheres. In this model, the FGM coating is divided into large number of perfectly bonded dissimilar thin sub-layers. Each sub-layer is homogeneous and has constant values of shear modulus and Poisson’s ratio. Analytical solutions are derived for the contact force and stress. New numerical results are presented to show the effects of graded shear modulus and graded Poisson’s ratio on the contact responses. It is shown that the contact force, stress and area can be significantly influenced by the graded shear modulus, and the effects of the graded Poisson’s ratio are not negligible.

1. Introduction Nonlinear point contact between two dissimilar elastic spheres is one of the fundamental models in contact mechanics and the theoretical bases of many practical applications. They include the design of ball bearing, the statistical modeling of rough surface contact [13], the analysis of nonlinear wave propagation in particle chain [25], as well as the discretized element modeling of granular material [8]. Since the seminal contribution by Hertz [18] for the normal frictionless point contact of two homogenous elastic spheres, the extension of Hertz’s theory to more complicated cases has been a subject of long-standing interests with both practical importance and scientific fascination. Surface reinforcement with thin coating has well been used as an effective technique to improve contact performances of structure components. Recent developments in the contact of coated elastic spheres can be found in the investigations of the contact stress [9], the plasticity evolution [6], the strengthening effect [12], the contact force and area laws [11,12], as well as the electrical resistance of coated elastic sphere [21]. The application of traditional coating has revealed that the material mismatch and discontinuity at the interface can lead to stress concentration, crack initiation and delamination. These disadvantages inspire the exploration of functionally graded materials (FGM) coating as an alternative.

∗

The FGM is a type of composite material characterized by the gradual variation of composition and microstructure along a specific direction. The resulting changes in material properties enable the FGM possess many unique features and free of the abovementioned weaknesses in traditional materials [28]. The accurate contact analysis of FGM coated materials with full consideration of actual material inhomogeneity is essential for the fundamental understanding of FGM reinforced contacting surface. To make the contact problem mathematically treatable, the shear modulus in FGM coating is usually assumed to vary according to exponential function while the Poisson’s ratio is assumed to be constant along the depth direction. This special type of the material inhomogeneity has been considered in many contact problems. They include the spherical indentation of FGM coated half space [22,24], the plane sliding contact [2,15], the plane rolling contact [1,16], the frictional plane spherical contact [14] and the receding contact [26,29]. This special model was also incorporated into boundary element modeling for the analysis of the frictional contact problem [17]. It is noted that the actual material gradations can be more complicated than that can be described by the exponentially graded shear modulus and constant Poisson’s ratio. To release this restriction, Ke and Wang [20] developed a linear multilayered model to simulate the arbitrarily graded shear modulus for plane frictionless indentation, which was extended by Liu and Xing [23] to axisymmetric problem. However, the Poisson’ ratio was limited to 1/3 in this improved model.

Corresponding author. E-mail address: [email protected] (Z.Q. Yue).

https://doi.org/10.1016/j.enganabound.2019.09.009 Received 30 March 2019; Received in revised form 11 July 2019; Accepted 11 September 2019 0955-7997/© 2019 Published by Elsevier Ltd.

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 1. Point contact between two elastic spheres with FGM coating: (a) Overview; (b) Geometric relation within the contact region.

𝑟2 ( ) . 2 𝑅2 + 𝑡2

According to the above literature survey, most of the previous studies are based on the assumptions that the modulus varies with depth according to special functions and the Poisson ratio keeps constant in depth. Such treatments need to be examined when addressing the contact problem with general forms of inhomogeneity in both modulus and Poisson’s ratio. Furthermore, the contact mechanics between two elastic spheres with FGM coatings receives less attention. The above two needs motivate the present study. This paper is organized as follows. In Section 2, the contact problem is formulated and the corresponding mixed boundary values problem is converted into a nonlinear Fredholm integral equation of the second kind. Section 3 develops a multi-layer half space model to account for arbitrarily graded FGM coatings and derives the kernel function in the integral equation. Section 4 introduces the numerical method for the solution of the nonlinear integral equation. Section 5 verifies the proposed model by comparing the Hertz’ solution for two uncoated spheres and the FEM results available in the literature for two identical spheres with homogeneous coatings, respectively. Section 6 presents selective numerical results to investigate the effects of graded modulus and Poisson’s ratio on the contact force, contact stress and contact area. Section 7 makes the concluding remarks on this study.

𝑓 2 (𝑟 ) =

2. Formulation of the contact problem

(i) 𝑢𝑧 (𝑟, 0) = 𝛿 − 𝑓 (𝑟), (ii) 𝜎𝑧𝑧 (𝑟, 0) = 0, (iii) 𝜎𝑟𝑧1 (𝑟, 0) = 𝜎𝑟𝑧1 (𝑟, 0) = 0,

The problem is examined under the following three assumptions. (1) The FGM coatings are very thin as compared with the curvature radiuses of two elastic spheres, i.e. t1 , t2 ≪ R1 ,R2 . (2) The two elastic spheres are in point contact under the small strain condition. (3) The material inhomogeneities are only along the radial thickness or depth direction according to two functions 𝜇 ci (z) and vci (z). Thus, the two contacting spheres can be approximated as two FGM coated half space [14] as illustrated in Fig. 2. Consequently, the surface deformations uz1 (r,0) and uz2 (r,0) (Fig. 1(b)) are equivalent to the surface settlements of FGM coated half spaces, as shown in Fig. 2(b) and (c), respectively. These assumptions enable the use of surface Green’s function of the half spaces to calculate the uz1 (r,0) and uz2 (r,0). 2.2. The mixed boundary values problem Referring to the compatibility relation (1) and the frictionless condition, the mixed boundary values problem for the contact problem under consideration can be stated as

2.1. Problem statement

(0 < 𝑟 ≤ 𝑎)

(1)

where, as illustrated in Fig. 1(b), uz1 (r,0) and uz2 (r,0) are the normal surface deformations of sphere 1 and 2 within the contact region (0 ≤ r ≤ a), respectively. 𝛿 1 and 𝛿 2 are the contact interferences for the two spheres, respectively. f1 (r) and f2 (r) are the two initial gaps given by 𝑓 1 (𝑟 ) =

𝑟≤𝑎; 𝑟>𝑎; 𝑟≥0.

(3)

where, uz (r,0), 𝛿, f(r) and 𝜎 zz (r,0) are the total contact deformation within the contact region, total contact interference, total initial gap and contact stress, respectively. They are given by

The problem of point contact between two elastic spheres with FGM coating is illustrated in Fig. 1. The Ri and ti denote the radius and coating thickness of the sphere, respectively; The 𝜇ci , vci and 𝜇 si , vsi are the shear modulus, Poisson’s ratio at coating surface and of the sphere, respectively. The subscript ‘c’ and ‘s’ refer to the ‘coating’ and ‘sphere’, respectively. The subscript i = 1, 2 denote the sphere 1 and the sphere 2, respectively. For the contact between the two deformable spheres, the following compatibility relation within the contact region (0 ≤ r ≤ a) should be satisfied 𝑢𝑧1 (𝑟, 0) + 𝑢𝑧2 (𝑟, 0) = 𝛿1 + 𝛿2 − 𝑓1 (𝑟) − 𝑓2 (𝑟),

(2)

𝑢𝑧 (𝑟, 0) = 𝑢𝑧1 (𝑟, 0) + 𝑢𝑧2 (𝑟, 0)

(4)

𝑓 (𝑟 ) = 𝑓 1 (𝑟 ) + 𝑓 2 (𝑟 ) ,

(5)

𝛿 = 𝛿 1 + 𝛿2 ,

(6)

𝜎𝑧𝑧 (𝑟, 0) = 𝜎𝑧𝑧1 (𝑟, 0) = 𝜎𝑧𝑧2 (𝑟, 0)

(7)

For mixed boundary values problem stated by (3), the continuous normal stresses 𝜎 zz (r,0) within the contact zone (0 < r ≤ a) are the unknowns to be determined for a given contact interference 𝛿. The resultant force can thus be calculated by

𝑟2 ( ), 2 𝑅1 + 𝑡1

𝑃 = 2𝜋 58

𝑎

∫0

𝜎𝑧𝑧 (𝑟, 0)𝑟𝑑𝑟

(8)

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 2. Approximation of the two contacting FGM coated spheres by FGM coated half spaces: (a) The two contacting FGM coated spheres; (b) FGM coated half space 1 for sphere 1; (c) FGM coated half space 2 for sphere 2.

2.3. Relations between the contact stress and surface settlement

be re-written in terms of the following dual integral equations [ ] +∞ ⎧ ̂ 𝑘𝑧 1 + ̂ 𝑘33 (𝜌) 𝜏3 (𝜌, 0)𝐽0 (𝜌𝑟)𝑑𝜌 = 𝛿 − 𝑓 (𝑟) , ⎪∫0 ⎨ +∞ ⎪ 𝜏3 (𝜌, 0)𝐽0 (𝜌𝑟)𝜌𝑑𝜌 = 0 , ⎩∫0

It is evident that the key issue for solving the contact problem under consideration are the surface Green’s functions relating the contact stresses 𝜎 zz1 (r,0), 𝜎 zz2 (r,0) to the surface settlements uz1 (r,0), uz2 (r,0) of the FGM coated half spaces. These relations can be established by solving the governing for vertically inhomogeneous half space, as shown in Fig. 2(b) and (c). This issue will be addressed in Section 3. For the convenience of demonstration, here we just give the following general expressions 𝑢𝑧1 (𝑟, 0)=𝑘𝑧1 𝑢𝑧2 (𝑟, 0)=𝑘𝑧2

∫0

+∞

[ ] 1 + 𝑘331 (𝜌) 𝜏3 (𝜌, 0)𝐽0 (𝜌𝑟)𝑑𝜌 ,

+∞

[ ] 1 + 𝑘332 (𝜌) 𝜏3 (𝜌, 0)𝐽0 (𝜌𝑟)𝑑𝜌 .

∫0

1 − 𝑣𝑐1 , 𝜇𝑐1

𝑘𝑧2 = −

1 − 𝑣𝑐2 . 𝜇𝑐2

where

̂ 𝑘33 (𝜌) =

+∞

∫0

𝜎𝑧𝑧 (𝑟, 0)𝐽0 (𝜌𝑟)𝑟𝑑𝑟

𝑘𝑧1 𝑘𝑧2 𝑘 (𝜌) + 𝑘 (𝜌) . 𝑘𝑧1 + 𝑘𝑧2 331 𝑘𝑧1 + 𝑘𝑧2 332

(13)

The dual integral Eq. (12) govern the unknown transformed normal stress 𝜏 3 (𝜌,0) and can be solved by introducing the following finite cosine transform of an auxiliary function 𝜙(x) (Sneddon, [27])

(9)

𝑎

𝜏3 (𝜌, 0) =

∫0

𝜙(𝑥) cos (𝜌𝑥)𝑑𝑥.

(14)

Using the expression (14), the last set of Eq. (12) is automatically satisfied and the dual integral Eq. (12) reduce to the Fredholm integral equation of the second kind for the unknown auxiliary function 𝜙(r) as below

(10)

𝜏 3 (𝜌,0) is the Hankel transform of the contact stresses 𝜎 zz (r,0), which is defined as 𝜏3 (𝜌, 0) =

(12) 𝑟 > 𝑎.

̂ 𝑘𝑧 = 𝑘𝑧1 + 𝑘𝑧2 ,

where, kz1 and kz2 are two material constants given by 𝑘𝑧1 = −

𝑟 ≤ 𝑎,

𝜙(𝑟) +

(11)

𝑎

∫0

𝑀 (𝑟, 𝑥)𝜙(𝑥)𝑑𝑥 = 𝜔(𝑟),

0 < 𝑟 < 𝑎,

(15)

where, Jm (𝜌r) is mth order Bessel function of the second kind. The integral variable 𝜌 is the transformed counterpart of radius r. The k331 (𝜌) and k332 (𝜌) are two dimensionless exponentially decreasing functions related to the material properties of the half space 1 and 2, respectively. They uniformly decay to zero as 𝜌 → +∞.

where, the inhomogeneous item 𝜔(r) is a known function related to the prescribed contact approach 𝛿 and the shape function f(r) (5) Fig. 2, it is given by [ ] 𝑟 𝑟𝑓 ′ (𝑠) 2 1 𝜔 (𝑟 ) = 𝛿− 𝑑𝑠 , (16) √ ∫0 𝜋̂ 𝑘𝑧 𝑟2 − 𝑠2

2.4. The integral equation

and the kernel M(r, x) is expressed in terms of the following improper integral

Substituting the integral representations (9) into the first equation of Eq. (3), the mixed boundary values problem defined by Eq. (3) can

𝑀 (𝑟, 𝑥) = 59

2 𝜋 ∫0

+∞

̂ 𝑘33 (𝜌) cos (𝜌𝑟) cos (𝜌𝑥)𝑑𝜌.

(17)

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

is denoted as the layer (n + 1)th layer with shear modulus 𝜇1(n + 1) and Poisson’s ratio v1(n + 1) . We further assume that all the sub layers are perfectly bonded to the adjacent layers. Thus, any arbitrarily graded 𝜇1 (z) and v1 (z) can be accurately approximated if the number of the sub-layers is sufficiently large. For this kind of discretized system, the constitutive Eq. (21) are written piecewise as

Substituting Eq. (14) into the Eq. (11), we can derive the solution of contact stress as 𝜎𝑧𝑧 (𝑟, 0) = √

𝑎

𝜙(𝑎) 𝑎2

−

𝑟2

−

∫𝑟

𝜙′ (𝑡) 𝑑𝑡 , √ 𝑡2 − 𝑟2

(𝑟 < 𝑎).

(18)

The resultant force can then be obtained by integrating the contact stress over the contact zone [0, a] 𝑃 = −2𝜋

𝑎

∫0

𝜎𝑧𝑧 (𝑟, 0)𝑟𝑑𝑟 = −2𝜋

𝑎

∫0

𝜙(𝑟)𝑑𝑟 .

( ) ⎧𝜎 = 2𝑣1(𝑗 ) 𝜇1(𝑗 ) 𝜕 𝑢𝑟1 + 𝑢𝑟1 + 𝜕 𝑢𝑧1 + 2𝜇 𝜕 𝑢𝑟1 , 1(𝑗 ) ⎪ 𝑟𝑟1 1 − 2𝑣1(𝑗 ) 𝜕𝑟 𝑟 𝜕𝑧 𝜕𝑟 ( ) ⎪ 2𝑣1(𝑗 ) 𝜇1(𝑗 ) 𝜕 𝑢𝑟1 𝑢𝑟1 𝜕 𝑢𝑧1 𝑢𝑟1 ⎪𝜎 = + + + 2𝜇1(𝑗 ) , ⎪ 𝜃𝜃1 1 − 2𝑣1(𝑗 ) 𝜕𝑟 𝑟 𝜕𝑧 𝑟 ( ) ⎨ 2𝑣1(𝑗 ) 𝜇1(𝑗 ) 𝜕 𝑢𝑟1 𝑢𝑟1 𝜕 𝑢𝑧1 𝜕 𝑢 ⎪𝜎 = + + + 2𝜇1(𝑗 ) 𝑧1 , ⎪ 𝑧𝑧1 1 −( 2𝑣1(𝑗 ) 𝜕𝑟 ) 𝑟 𝜕𝑧 𝜕𝑧 ⎪ 𝜕 𝑢𝑧1 𝜕 𝑢𝑟1 ⎪𝜎𝑟𝑧1 = 𝜇1(𝑗 ) + , ⎩ 𝜕𝑟 𝜕𝑧

(19)

Clearly, once the k331 (𝜌) and k332 (𝜌) are given, the contact problem can be readily solved form the Fredholm integral Eq. (15). 3. Solution of FGM coated half space subjected to surface loading In this section, the problem of FGM coated half space subjected to surface loading, as discussed in Section 2.3, is addressed. For convenience, only the half space 1 (Fig. 2(b)) is analyzed. In the ensuing, the last subscript ‘1′ without the bracket ‘()’ refers to the half space 1.

(24) According to the perfectly bonded assumption, the following interface conditions are considered ( ) ( ) ( ) 𝑢𝑧1 𝑟, 𝐻1−(𝑗 ) = 𝑢𝑧1 𝑟, 𝐻1+(𝑗 ) = 𝑢𝑧1 𝑟, 𝐻1(𝑗 ) ( ) ( ) ( ) 𝑢𝑟1 𝑟, 𝐻1−(𝑗 ) = 𝑢𝑟1 𝑟, 𝐻1+(𝑗 ) = 𝑢𝑟1 𝑟, 𝐻1(𝑗 ) ( ) ( ) (25) ( ) 𝜎𝑧𝑧1 𝑟, 𝐻1−(𝑗 ) = 𝜎𝑧𝑧1 𝑟, 𝐻1+(𝑗 ) = 𝜎𝑧𝑧1 𝑟, 𝐻1(𝑗 ) ( ) ( ) ( ) 𝜎𝑟𝑧1 𝑟, 𝐻1−(𝑗 ) = 𝜎𝑟𝑧1 𝑟, 𝐻1+(𝑗 ) = 𝜎𝑟𝑧1 𝑟, 𝐻1(𝑗 )

3.1. Governing equations The equilibrium equations for the axisymmetric contact problem under consideration are ⎧ 𝜕 𝜎𝑟𝑟1 + 𝜕 𝜎𝑟𝑧1 + 𝜎𝑟𝑟1 − 𝜎𝜃𝜃1 = 0 , ⎪ 𝜕𝑟 𝜕𝑧 𝑟 ⎨ 𝜕𝜎 ⎪ 𝑟𝑧1 + 𝜎𝑟𝑧1 + 𝜕 𝜎𝑧𝑧1 = 0 . ⎩ 𝜕𝑟 𝑟 𝜕𝑧

(20)

3.3. Solution representations

In this study, we assume that both the shear modulus and Poisson’s ratio vary arbitrarily along the thickness direction of the FGM coating. Therefore, the constitutive equations for such inhomogeneous half space can be written as ( ) ⎧𝜎 = 2𝑣1 (𝑧)𝜇1 (𝑧) 𝜕 𝑢𝑟1 + 𝑢𝑟1 + 𝜕 𝑢𝑧1 + 2𝜇 (𝑧) 𝜕 𝑢𝑟1 , 1 ⎪ 𝑟𝑟1 1 − 2 𝑣 1 (𝑧 ) 𝜕𝑟 𝑟 𝜕𝑧 𝜕𝑟 ( ) ⎪ 𝜕 𝑢 2 𝑣 𝑧 𝜇 𝑧 𝜕 𝑢 𝑢 𝑢𝑟1 ( ) ( ) 𝑧1 1 1 𝑟1 𝑟1 ⎪𝜎 = + + + 2 𝜇 1 (𝑧 ) , ⎪ 𝜃𝜃1 1 − 2 𝑣 1 (𝑧 ) 𝜕𝑟 𝑟 𝜕𝑧 𝑟 ( ) (21) ⎨ 𝜕 𝑢𝑧1 2𝑣1 (𝑧)𝜇1 (𝑧) 𝜕 𝑢𝑟1 𝑢𝑟1 𝜕 𝑢𝑧1 ⎪𝜎𝑧𝑧1 = + + + 2𝜇1 (𝑧) , 1 − 2 𝑣 1 (𝑧 ) 𝜕𝑟 𝑟 𝜕𝑧 𝜕𝑧 ⎪ ( ) ⎪ 𝜕 𝑢𝑧1 𝜕 𝑢𝑟1 ⎪𝜎𝑟𝑧1 = 𝜇1 (𝑧) + , ⎩ 𝜕𝑟 𝜕𝑧

The multi-layer half space governed by the Eqs. (20), (24) and subjected to the boundary conditions (23) and interface conditions (25) can be solved following the General Kelvin’s Solution for the ThreeDimensional N-Layers Elasticity [30–32]. Firstly, the following Hankel’s transforms are employed to represent the field components ⎡𝑤11 (𝜌, 𝑧)⎤ ⎡ −𝑢𝑟1 (𝑟, 𝑧)𝐽1 (𝜌𝑟)𝜌⎤ ⎢ ⎥ ⎥ +∞ ⎢ 𝑤31 (𝜌, 𝑧)⎥ ⎢ ⎢ 𝑢𝑧1 (𝑟, 𝑧)𝐽0 (𝜌𝑟)𝜌 ⎥𝑟𝑑𝑟 , 𝐔 1 (𝑧 ) = = ⎢𝜏31 (𝜌, 𝑧) ⎥ ∫0 ⎢ 𝜎𝑧𝑧1 (𝑟, 𝑧)𝐽0 (𝜌𝑟) ⎥ ⎢𝜏 (𝜌, 𝑧) ⎥ ⎢−𝜎 (𝑟, 𝑧)𝐽 (𝜌𝑟) ⎥ 1 ⎣ 11 ⎦ ⎣ 𝑟𝑧1 ⎦

𝑧 = 0. 0 < 𝑧 < 𝑡1 . 𝑡1 ≤ 𝑧 .

(22)

𝑢𝑟1 (𝑟, 𝑧)= −

The surface settlements uz1 (r,0) due to the 𝜎 zz1 (r,0) can be determined by solving the Eq. (20)–(21) under the following boundary conditions 𝜎𝑧𝑧1 (𝑟, 0) = 𝜎𝑧𝑧 (𝑟, 0), 𝜎𝑟𝑧1 (𝑟, 0) = 0,

0 ≤ 𝑟 ≤ 𝑎; 0 ≤ 𝑟;

(26)

where, w11 (𝜌,z), w31 (𝜌,z), 𝜏 31 (𝜌,z) and 𝜏 11 (𝜌,z) are the transformed counterparts of the displacements ur1 (r,z), uz1 (r,z), the vertical stress 𝜎 zz1 (r,z) and plane shear stress 𝜎 rz1 (r,z), respectively. Jm (𝜌r) is mth order Bessel function of the second kind. The integral variable 𝜌 is the transformed counterpart of radius r. Thus, the elastic fields in physical domain can be expressed by using the reciprocal property of Eq. (26)

where, the 𝜇 1 (z) and v1 (z) are the graded shear modulus and Poison’s ratio within the entire range z ∈ [0, +∞]. Specifically, they can be expressed separately as ⎧{𝜇 ; 𝑣 }, 𝑐1 { } ⎪{ 𝑐1 } 𝜇1 (𝑧) ; 𝑣1 (𝑧) = ⎨ 𝜇𝑐1 (𝑧) ; 𝑣𝑐1 (𝑧) , { } ⎪ 𝜇𝑠 1 ; 𝑣 𝑠 1 , ⎩

𝐻𝑗 < 𝑧 < 𝐻𝑗+1

𝑢𝑧1 (𝑟, 𝑧)=

∫0

+∞

∫0

𝜎𝑧𝑧1 (𝑟, 𝑧) =

(23)

+∞

3.2. A multi-layer half space model for contact stress and displacements

𝑤31 (𝜌, 𝑧)𝐽0 (𝜌𝑟)𝑑𝜌,

+∞

∫0

𝜎𝑟𝑧1 (𝑟, 𝑧) = −

𝑤11 (𝜌, 𝑧)𝐽1 (𝜌𝑟)𝑑𝜌,

∫0

(27) 𝜏31 (𝜌, 𝑧)𝐽0 (𝜌𝑟)𝜌𝑑𝜌,

+∞

𝜏11 (𝜌, 𝑧)𝐽1 (𝜌𝑟)𝜌𝑑𝜌.

Using the Hankel’s transforms (26), the governing Eqs. (20) and (24) can be converted into the following ODE for the layer j

In general, it is difficult to solve the governing Eqs. (20) and (21) directly under arbitrary 𝜇 1 (z) and v1 (z). To tackle this problem, we propose a multi-layer half space model to approximate the material inhomogeneity. Take the FGM coated half space 1 for instance, the multi-layer half space model is illustrated in Fig. 3. As it is shown in Fig. 3(b), in the multi-layer half space model the FGM coating is divided into large number of homogeneous thin sublayers. The 𝜇 1(j) , v1(j) , h1(j) and H1(j) denote the shear modulus, Poisson’s ratio, layer thickness and vertical coordinate of the bottom surface of layer No. j (j = 1, 2, …, n + 1), respectively The elastic substrate

𝑑 𝐔 (𝑧) = 𝜌𝐂1(𝑗 ) 𝐔1 (𝑧) , 𝑑𝑧 1

(28)

where, the C1(j) is the matrix of elastic constants defined as

𝐂1(𝑗 )

60

⎡ 0 ⎢ ⎢ = ⎢ 1 − 2𝛼1(𝑗 ) ⎢ ⎢ 0 ( ) ⎢ ⎣4𝜇1(𝑗 ) 1 − 𝛼1(𝑗 )

−1

0

0

𝜇1(𝑗 ) 0 2𝛼1(𝑗 ) − 1

0 0

𝛼1(𝑗 )

1 ⎤ 𝜇1(𝑗 ) ⎥ ⎥ 0 ⎥, ⎥ 1 ⎥ ⎥ 0 ⎦

(29)

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 3. A multi-layer model to simulate the arbitrarily graded FGM coating: (a) Material properties and geometry parameters; (b) Approximated material ditribution profile in the multi-layer model. Fig. 4. Dimensional contact stresses given by presnet analytical solutios and Hertz’s soltuion for the contact of two dismmialr spheres under different radius ratio R2 /R1 .

and 𝛼 1(j) = (1 − 2v1(j) )/2(1 − v1(j) ) is generalized Poisson’s ratio.

In Eq. (31), Qp1(j) , Qq1(j) , Rp1(j) and Rq1(j) are the matrixes of elastic constants for layer j, which are given by

3.4. Solutions in transformed domain By solving the ODE (28) analytically and considering the perfectly − ) = 𝐔 (𝐻 + ) = 𝐔 (𝐻 bonded conditions 𝐔1 (𝐻1( 1 1 1(𝑗) ) between the layers, 𝑗) 1(𝑗) one can obtained the solution of U1 (z) as below [30–32] { 𝐔 1 (𝑧 ) =

( ) ( ) 𝑒 𝐐1(𝑗 ) 𝐻1(𝑗 ) − 𝑧 𝐔 𝐻1(𝑗 ) , )[ ( ( ) ] ( ) 𝑒𝜌 𝐻1(𝑛) −𝑧 𝐐𝑞1(𝑛+1) + 𝜌 𝑧 − 𝐻𝑛(𝑗 ) 𝐑𝑞1(𝑛+1) 𝐔1 𝐻1(𝑛) , ) ( 𝜌 𝐻1(𝑗 ) −𝑧

𝐐𝑞1(𝑗 )

𝑗 < 𝑛+1 , 𝑗 = 𝑛+1 . (30)

where, the Q1(j) (𝜒) is the transfer matrix for layer j defined as ( ) 𝐐1(𝑗 ) (𝜒) = 𝐐𝑞1(𝑗 ) − 𝜌𝜒𝐑𝑞1(𝑗 ) + 𝑒−2𝜌𝜒 𝐐𝑝1(𝑗 ) − 𝜌𝑧𝐑𝑝1(𝑗 ) .

𝐐 𝑝 1 (𝑗 ) (31) 61

⎡ 1 ⎢ 1⎢ 𝛼1(𝑗 ) = ⎢ 2⎢ 0 ( ) ⎢ ⎣−2𝜇1(𝑗 ) 1 − 𝛼1(𝑗 ) ⎡ 1 ⎢ 1⎢ −𝛼1(𝑗 ) = ⎢ 2⎢ 0 ( ) ⎢ ⎣2𝜇1(𝑗 ) 1 − 𝛼1(𝑗 )

𝛼1 ( 𝑗 ) 1 ( ) −2𝜇1(𝑗 ) 1 − 𝛼1(𝑗 ) 0 −𝛼1(𝑗 ) 1 ( ) 2𝜇1(𝑗 ) 1 − 𝛼1(𝑗 ) 0

0 −

−

1+𝛼1(𝑗 )

−𝛼1(𝑗 ) 1

1 −𝛼1(𝑗 ) 1+𝛼1(𝑗 ) 2𝜇1(𝑗 )

1 𝛼1(𝑗 )

2𝜇1(𝑗 )

0

2𝜇1(𝑗 )

0

1+𝛼1(𝑗 )

1+𝛼1(𝑗 ) 2𝜇1(𝑗 )

0 𝛼1(𝑗 ) 1

⎤ ⎥ ⎥ ⎥, ⎥ ⎥ ⎦

⎤ ⎥ ⎥ ⎥, ⎥ ⎥ ⎦

(32)

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 5. The frictionless point contact between two identical spheres with homogeneous coatings (a) and a sphere with homogeneous coating compressed by a rigid smooth plane (b). Fig. 6. Comparison between the dimenssional P ∼ 𝛿 relations given by FEM analysis [10] and present analytical solution for the normal contact of two identical coated spheres with radius R = 10mm, Young’s modulus ratio Ec /Es = 4, under different thickness ratio.

governed the two sets of boundary values U1 (0) and U1 (Hn )

and ⎡ −1 )⎢ 1 − 𝛼1 ( 𝑗 ) ⎢ 1 = ⎢ 2 ⎢−2𝜇1(𝑗 ) ⎢ 2𝜇 ⎣ 1(𝑗 ) (

𝐑𝑞1(𝑗 )

⎡ 1 )⎢ 1 − 𝛼1 ( 𝑗 ) ⎢ 1 = ⎢ 2 ⎢2𝜇1(𝑗 ) ⎢2𝜇 ⎣ 1(𝑗 ) (

𝐑 𝑝 1 (𝑗 )

−1 1 −2𝜇1(𝑗 ) 2𝜇1(𝑗 )

1 2𝜇1(𝑗 ) − 2𝜇1 1(𝑗 )

1 −1

−1

− 2𝜇1

−1

− 2𝜇1 1(𝑗 )

−2𝜇1(𝑗 ) −2𝜇1(𝑗 )

1(𝑗 )

−1 −1

( ) 𝐔1 (0) = 𝑒𝜌𝑡1 𝛀1 𝐔1 𝐻1(𝑛) , ( ) ( ) 𝐔1 𝐻1(𝑛) = 𝐐𝑞1(𝑛+1) 𝐔1 𝐻1(𝑛) .

1 ⎤ 2𝜇1(𝑗 ) ⎥ 1 ⎥ − 2𝜇 1(𝑗 ) ⎥,

1 ⎥ −1 ⎥⎦

1 ⎤ 2𝜇1(𝑗 ) ⎥ 1 ⎥ 2𝜇1(𝑗 ) ⎥.

1 ⎥ 1 ⎥⎦

where, t1 = h1(1) + h1(2) + ⋅⋅⋅ + h1(n) is the thickness of the FGM layer 1 and 𝛀1 is the production of transfer matrixes from the top layer to the last layer as follow ( ) ( ) ( ) ( ) 𝛀1 = 𝐐1(1) ℎ1(1) 𝐐1(2) ℎ1(2) ⋯ 𝐐1(𝑛−1) ℎ1(𝑛−1) 𝐐1(𝑛) ℎ1(𝑛) .

(33)

(36)

Eliminating the 4 linear dependent equations in Eq. (35), the boundary values U1 (H1(n) ) can be solved in terms of transformed contact stress 𝜏 31 (𝜌,0) as below

Thus, using the Eq. (30), the U1 (z) within the layer j (Hj − 1 ≤ z ≤ Hj ) can be related to the U1 (Hn ) at the last interface by ) ] [( ( ) 𝐔1 (𝑧) = 𝑒𝜌 𝐻1(𝑗 ) −𝑧 +ℎ1(𝑗+1) +ℎ1(𝑗+2) +⋯+ℎ1(𝑛) 𝐐1(𝑗 ) 𝐻1(𝑗 ) − 𝑧 𝐐1(𝑗+1) ( ) ( ) ( ) ( ) × ℎ1(𝑗+1) 𝐐1(𝑗+2) ℎ1(𝑗+2) ⋯ 𝐐1(𝑛) ℎ1(𝑛) 𝐔 𝐻1(𝑛) .

(35)

( ) 𝐌1 𝐔1 𝐻1(𝑛) = 𝑒−𝜌𝑡1 𝐈𝑢 𝜏31 (𝜌, 0) ,

(34)

(37)

where, [

Substituting z = 0 into the Eq. (34) and z = H1(n) into the second equations of Eq. (30) we can have the following 8 linear algebra equations

𝐌1 = 62

] 𝐏𝑝1 , 𝐏𝑞1 𝛀1

[ 𝐈𝑢 = 0

0

1

0

]𝑇

,

(38)

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 7. Variation of material properties along the radius direction in the coating: (a) Hard coating; (b) Soft coating. The z is the localized coordinate to decrivbe the graded shear modulus and Poisson’s ratio in each sphere. The subscript “i” refers to the sphere number. The 𝜇𝑐(𝑖) (𝑧) and 𝑣𝑐(𝑖) (𝑧) in the coating (𝑅𝑖 ≤ 𝑧 ≤ 𝑅𝑖 + 𝑡𝑖 ) under specific gradation pattern, e.g. Step, Linear, Exp or Sine can be determined by using the variation function 𝐹 (𝑧) defined in Table 1.

Fig. 8. Dimensional load-interference (P ∼ 𝛿) relation for the contact of two disimilar FGM coated spheres for Case (a)∼(d) in Table 2.

and Pq1 , Pp1 are given by [ 𝐏𝑞1 =

𝐏𝑝1

0 0

0 0

1 0

⎡ 1 ⎢ =⎢ ⎢1 − 𝛼1(𝑛+1) ⎣

Substituting the U1 (H1(n) ) given by Eq. (37) into Eq. (34) we can obtain the solutions of U1 (z) as

] 0 , 1 −1 1 − 𝛼1(𝑛+1)

−

1

2𝜇1(𝑛+1) 1 + 𝛼1(𝑛+1) 2𝜇1(𝑛+1)

[

⎤ 2𝜇1(𝑛+1) ⎥ 1 + 𝛼1(𝑛+1) ⎥. ⎥ 2𝜇1(𝑛+1) ⎦

(

)]

𝐔1 (𝑧) = 𝑒−𝜌 ℎ1(1) +ℎ1(2) +⋯+ℎ1(𝑗−1) + 𝑧−𝐻1(𝑗−1) ( ) ( ) ( ) ×𝐐1(𝑗 ) 𝐻1(𝑗 ) − 𝑧 𝐐1(𝑗+1) ℎ1(𝑗+1) 𝐐1(𝑗+2) ℎ1(𝑗+2) ⋯ ( ) −1 ×𝐐1(𝑛) ℎ1(𝑛) 𝐌 𝐈𝑢 𝜏31 (𝜌, 0).

1

(39)

63

(40)

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 9. Dimensional load-interference (P ∼ 𝛿) relation for the contact of two disimilar FGM coated spheres for Case (e)∼(h) in Table 3.

The U1 (z) can also be expressed by

Table 1 Variation function 𝜇𝑐(𝑖) (𝑧) and 𝑣𝑐(𝑖) (𝑧) for different types of material gradation. Gradation pattern Step Linear Exp Sine

Variation function 𝜇𝑐(𝑖) (𝑧) and 𝑣𝑐(𝑖) (𝑧) 𝜇𝑐(𝑖) (𝑧) = 𝜇𝑐(𝑖) 𝑣𝑐(𝑖) (𝑧) = 𝑣𝑐(𝑖) 𝜇𝑐(𝑖) (𝑧) = 𝜇𝑠(𝑖) [1+𝛼(𝜇𝑐(𝑖) )( 𝑅𝑧 − 1)] 𝑣𝑐(𝑖) (𝑧) =

𝜇𝑐(𝑖) (𝑧) = 𝑣𝑐(𝑖) (𝑧) =

𝑖

𝑣𝑠(𝑖) [1+𝛼(𝑣𝑐(𝑖) )( 𝑅𝑧 𝑖

− 1)]

𝜇𝑠(𝑖) Exp[𝛼(𝜇𝑐(𝑖) )( 𝑅𝑧 𝑖 𝑣𝑠(𝑖) Exp[𝛼(𝑣𝑐(𝑖) )( 𝑅𝑧 𝑖

Inhomogeneity index 𝛼(𝜇𝑐(𝑖) ) and 𝛼(𝑣𝑐(𝑖) )

𝛼(𝜇𝑐(𝑖) ) = 𝛼(𝑣𝑐(𝑖) ) =

− 1)]

𝛼(𝜇𝑐(𝑖) ) =

− 1)]

𝛼(𝑣𝑐(𝑖) ) =

𝜇𝑐(𝑖) (𝑧) = 𝑣𝑐(𝑖) (𝑧) =

− 1)+1]

𝛼(𝑣𝑐(𝑖) ) =

𝑅𝑖 𝑡𝑖 𝑅𝑖 𝑡𝑖 𝑅𝑖 𝑡𝑖 𝑅𝑖 𝑡𝑖 𝑅𝑖 𝑡𝑖 𝑅𝑖 𝑡𝑖

𝜇

( 𝜇𝑐(𝑖) − 1) 𝑠(𝑖)

𝑣𝑐(𝑖)

(𝑣

− 1)

𝑠(𝑖) 𝜇𝑐(𝑖)

In 𝜇 𝑣

𝑠(𝑖)

In 𝑣𝑐(𝑖) 𝑠(𝑖)

𝜇

( 𝜋2 sin−1 𝜇𝑐(𝑖) − 1) 𝑣

3.5. Solutions in physical domain

𝑠(𝑖)

( 𝜋2 sin−1 𝑣𝑐(𝑖) − 1) 𝑠(𝑖)

Substituting Eq. (41) into Eq. (27), we have the elastic fields within the layered half space expressed by

− 1)+1]

(Remark: As illustrated in Fig. 7, the 𝜇𝑐(𝑖) , 𝑣𝑐(𝑖) , 𝜇𝑠(𝑖) , 𝑣𝑠(𝑖) are shear modulus and Poisson’s ratio at the coating surface and of the sphere, respectively)

𝑢𝑟1 (𝑟, 𝑧) = −

Table 2 Parametric combinations used to investigate the effect of modulus gradation. Case 𝜇 c 1 /𝜇 s 1 𝜇 c 2 /𝜇 s 2

(a) 4 5

(b) 6 10

(c) 0.1 0.2

(41)

where, Φ131 (𝜌,z), Φ331 (𝜌,z), Ψ131 (𝜌,z) and Ψ331 (𝜌,z) are kernel functions connecting the elastic fields and contact stress in transformed domain for sphere 1. They can be determined by the simple matrix production as shown in Eq. (40).

/

𝛼(𝜇𝑐(𝑖) ) = 𝜇𝑠(𝑖) sin 𝜋2 [𝛼(𝜇𝑐(𝑖) )( 𝑅𝑧 𝑖 𝑣𝑠(𝑖) sin 𝜋2 [𝛼(𝑣𝑐(𝑖) )( 𝑅𝑧 𝑖

[ ]𝑇 𝐔1 (𝑧) = Φ131 (𝜌, 𝑧), Φ331 (𝜌, 𝑧), Ψ331 (𝜌, 𝑧), Ψ131 (𝜌, 𝑧) 𝜏31 (𝜌, 0) .

𝑢𝑧1 (𝑟, 𝑧) =

(d)

∫0

+∞

∫0

𝜎𝑧𝑧1 (𝑟, 𝑧) =

0.25 0.375

+∞

𝜎𝑟𝑧1 (𝑟, 𝑧) = −

(Remarks: The other parameters are fixed as R1 = 10mm, R2 = 15mm, t1 = 0.05mm, t2 = 0.2mm, vc1 = vs1 = vc2 = vs2 = 0.3 and 𝜇s1 = 𝜇s2 = 200GPa). 64

Φ331 (𝜌, 𝑧)𝜏31 (𝜌, 0)𝐽0 (𝜌𝑟)𝑑𝜌 ,

+∞

∫0 ∫0

Φ131 (𝜌, 𝑧)𝜏31 (𝜌, 0)𝐽1 (𝜌𝑟)𝑑𝜌 ,

(42) Ψ331 (𝜌, 𝑧)𝜏31 (𝜌, 0)𝐽0 (𝜌𝑟)𝜌𝑑𝜌 ,

+∞

Ψ131 (𝜌, 𝑧)𝜏31 (𝜌, 0)𝐽1 (𝜌𝑟)𝜌𝑑𝜌 .

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 10. Dimensional interference-contact area (𝛿 ∼ 𝜋a2 ) relation for the contact of two disimilar FGM coated spheres for Case (a)∼(d) in Table 2.

where, k131 (𝜌), k331 (𝜌) are the convergent parts of Φ131 (𝜌,0), Φ331 (𝜌,0). They reduce to zeros for homogeneous half space and have the following asymptotic values as 𝜌 → ∞

Table 3 Parametric combinations used to investigate the effect of Poisson’s ratio gradation. Case

(e)

(f)

(g)

(h)

vc1 vs1 vc2 vs2

0.5 0.05 0.5 0.1

0.5 0.1 0.5 0.2

0.05 0.5 0.1 0.5

0.1 0.5 0.2 0.5

lim 𝑘 (𝜌) 𝜌→∞ 131

𝑘𝑟1 = lim Φ131 (𝜌, 0) = − 𝜌→∞

𝑢𝑧1 (𝑟, 0) =

+∞

+∞

∫0

𝜎𝑧𝑧1 (𝑟, 0) =

𝜎𝑟𝑧1 (𝑟, 0) = −

Φ131 (𝜌, 0)𝜏31 (𝜌, 0)𝐽1 (𝜌𝑟)𝑑𝜌,

Ψ331 (𝜌, 0)𝜏31 (𝜌, 0)𝐽0 (𝜌𝑟)𝜌𝑑𝜌,

+∞

∫0

Ψ131 (𝜌, 0)𝜏31 (𝜌, 0)𝐽1 (𝜌𝑟)𝜌𝑑𝜌 .

Ψ131 (𝜌, 0) = 0 ,

[ ] Φ131 (𝜌, 0) = 𝑘𝑧1 1 + 𝑘331 (𝜌) ,

1 − 2𝑣1(1) 2𝜇1(1)

,

𝑘𝑧1 = lim Φ331 (𝜌, 0) = − 𝜌→∞

1 − 𝑣1(1) 𝜇1(1)

4. Numerical implement As demonstrated, the 𝑘̂ 33 (𝜌) in the kernel function M(r, x) defined by Eq. (17) is given implicitly via the matrix manipulations (40) and Eqs. (13), (44). Such implicit expression makes the solution of the Fredholm integral Eq. (15) inevitably numerical. There are several techniques available for the purpose, for instance, the standard discretization approach [3,5] and the meshless method [4,34–36]. In the present analysis, the former approach is adopted and the Fredholm integral

(43)

It is noted that the Ψ331 (𝜌,0), Ψ131 (𝜌,0), Φ131 (𝜌,0) and Φ331 (𝜌,0) have the following mathematical properties [30–32] Ψ331 (𝜌, 0) = −1, [ ] Φ131 (𝜌, 0) = 𝑘𝑟1 1 + 𝑘131 (𝜌) ,

(45)

where, 𝜇 1(1) = 𝜇 1c and v1(1) = v1c are the shear modulus and Poisson’s ratio of the top sub-layer of sphere 1, respectively. In this section, we have outlined the procedure to determine the uz1 (r,0) due to the 𝜎 zz1 (r,0) for sphere 1. Similarly, the uz2 (r,0) due to equal contact stress 𝜎 zz (r,0) for sphere 2 can be determined by using the corresponding parameters 𝜇 2(j) , v2(j) , h2(j) , H2(j) instead.

Φ331 (𝜌, 0)𝜏31 (𝜌, 0)𝐽0 (𝜌𝑟)𝑑𝜌,

+∞

∫0

=0.

(46)

Thus, substituting z = 0 into Eq. (42), the surface responses can be obtained as ∫0

lim 𝑘 (𝜌) 𝜌→∞ 331

The kr1 , kz1 are the constant limits given by

(Remarks: The other parameters are fixed as R1 = 10mm, R2 = 15mm , t1 = 0.05mm, t2 = 0.2mm, and 𝜇c1 = 𝜇c2 = 𝜇s1 = 𝜇s2 = 200GPa).

𝑢𝑟1 (𝑟, 0) = −

= 0,

(44) 65

.

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 11. Dimensional interference-contact area (𝛿 ∼ 𝜋a2 ) relation for the contact of two disimilar FGM coated spheres for Case (e)∼(h) in Table 3.

value 𝜌0 can be used as the up limit of the improper integral (17). Based on the numerical test, we find that if the variable change 𝜉 = a𝜌 (a is the contact radius) is introduced in Eq. (17), the truncated value 𝜉 0 = 10 is sufficient to generate satisfactory results by using global adaptive integral technique, e.g. recursive adaptive Simpson quadrature ‘quadv’ available in MATLAB. Since contact stress should vanish at the contact edge r = a, the closure condition 𝜙(a) = 0 is required to eliminate the square root singularity in Eq. (18). Therefore, the Fredholm integral Eq. (15) is essentially nonlinear in that the contact radius a is an unknown in advance and needs to be searched iteratively. It is found that the 𝜙(a) is a monotonic function √ of a. Hence, the classical dichotomy with the starting point 𝑎 = 𝑅𝛿 as predicted by the Hertz’s solution can be used for the iteration [7].

Eq. (15) is converted to a system of linear algebraic equation as below. (𝐈 + 𝐌)𝛟 = 𝛚 , where [ ( ) 𝜙 = 𝜙 𝑥1 [ ( ) 𝜔 = 𝜔 𝑟1

( ) 𝜙 𝑥2 ( ) 𝜔 𝑟2

(47)

⋯ ⋯

( )]𝑇 , 𝜙 𝑥𝑚 ( )]𝑇 , 𝜔 𝑟𝑚

(48)

and ( ) ⎡ 𝑀 𝑟1 , 𝑥1 ( ) ⎢ 𝑎 𝑀 𝑟2 , 𝑥1 𝐌= ⎢ 𝑚⎢ ( ⋮ ⎢𝑀 𝑟 , 𝑥 ) 𝑚 1 ⎣

( ) 𝑀 𝑟1 , 𝑥2 ( ) 𝑀 𝑟2 , 𝑥2 ⋮ ( ) 𝑀 𝑟𝑚 , 𝑥2

⋯ ⋯ ⋱ ⋯

( ) 𝑀 𝑟1 , 𝑥𝑚 ⎤ ( )⎥ 𝑀 𝑟2 , 𝑥𝑚 ⎥ . ⎥ ⋮ ( )⎥ 𝑀 𝑟𝑚 , 𝑥𝑚 ⎦

(49)

In the discretized forms (47)–(49), the solution interval r ∈ [0, a) is divided into m segments of equal length .The collocation points are selected as the mid-point of each segment as ri = xi = a(2i − 1)/2m. It was shown that accurate results can be obtained by increasing the number of the segments [3,5]. The m can be chosen such that further increase of m has negligible influence on the results. The determination of the M(ri ,xj ) requires the numerical evaluation of the improper integral (17). Since the integral kernel 𝑘̂ 33 (𝜌) is an exponentially decreasing function governed by |k33 (𝜌)| < B(1 + 𝜌H)e−𝜌b [30,31,33], where B, b and H are constants related to the material properties and coating thickness of the coated half space, the convergence of improper integral (17) is ensured. Hence, a sufficiently large truncate

5. Verifications 5.1. Analytical degeneration to Hertz’s solution for two uncoated spheres If the FGM coatings become homogeneous and possess the same material properties as the bulk spheres, as discussed in Section 3.5, the 𝑘̂ 33 (𝜌) = 0. Then, the kernel function M(r, x) of the Fredholm integral Eq. (15) vanishes and the auxiliary function 𝜙(r) can be solved in exact closed form as follow ( ) 2 1 𝑟2 𝜙(𝑟) = 𝜔(𝑟) = 𝛿− , 0<𝑟<𝑎. (50) 𝜋̂ 𝑅 𝑘 𝑧

66

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 12. Dimensional contact stress 𝜎 zz (r,0) for the contact of two disimilar FGM coated spheres for Case (a)∼(d) in Table 2.

Considering the closure condition 𝜙(a) = 0 and substituting 𝜙(r) into Eqs. (18) and (19) leads to 𝜎𝑧𝑧 (𝑟, 0) = −

𝑃 =𝐾

𝑎3 , 𝑅

3 1 1√ 2 𝑎 − 𝑟2 , 2𝜋 𝐾 𝑅

5.3. Numerical comparison to the FEM results for two identical spheres with homogeneous coating The correctness of the present analytical solutions is further confirmed by comparing to the existing results of coated materials. The model of a coated elastic sphere compressed by a smooth rigid flat considered by Goltsberg and Etsion [10] as shown in Fig. 5(b) is used as example. It is noted that, this selected example is equivalent to the normal frictionless contact of two identical coated spheres, which is a special case of present contact model (Fig. 5(a)). In this example, the homogenous coating is discretized into 100 identical perfectly bonded elastic layers. The comparison is made for the P ∼ 𝛿 relation under different thickness ratio t/R. The results are plotted in Fig. 6 for the case of Es = 200GPa, Ec = 800GPa and R = 10mm. As it is observed the results predicted by present analytical solution are in exact agreement with those from FEM simulation [10] for a wide range of thickness ratio.

(51)

(52)

where, 𝐾 = (3∕8)𝑘𝑧 = (3∕4)[(1 − 𝑣21 )∕𝐸1 + (1 − 𝑣22 )∕𝐸2 ] is the equivalent modulus in classical Hertz’s theory. Here, v1 , E1 are the Poisson’s ratio and Young’s modulus for the sphere 1, and v2 , E2 are those for sphere 2. The solution (51) and (52) are consistent with Hertz’s solution for the contact of two homogeneous spheres [19]. 5.2. Numerical comparison to Hertz’s solution for two uncoated spheres

6. Numerical results and discussions

A numerical comparison to Hertz’s solution is also presented herein to show the precision of our code in calculating contact stress. In the numerical example, the contact of two dissimilar elastic sphere is considered. To model homogeneous sphere, the coating of each sphere is discretized into 100 identical layers with the same material properties of the sphere. The dimensional contact stress due to the interference D = 0.0002mm for different radius ratio R2 /R1 are plotted in Fig. 4. It is shown that, the numerical results given by present degenerated solution are in exact agreement with those of Hertz’s solution.

6.1. General In this section, selective numerical results of resultant force, contact stress and contact area. for the contact between two dissimilar FGM coated spheres are presented. The attention is focused on the effects of shear modulus and Poisson’s ratio gradations. To this end, four variation patterns as shown in Fig. 7, namely Step change (homogenous coating), linear, exponential and sine gradation, are considered. The variation 67

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

Fig. 13. Dimensional contact stress 𝜎 zz (r,0) for the contact of two disimilar FGM coated spheres for Case (e)∼(h) in Table 3.

functions 𝜇𝑐(𝑖) (𝑧) and 𝑣𝑐(𝑖) (𝑧) to describe the shearmodulus or Poisson’s ratio in the coating (𝑅𝑖 ≤ 𝑧 ≤ 𝑅𝑖 + 𝑡𝑖 ) for the four gradation patterns are given in Table 1. For simplicity, we assume that the two FGM coatings have the same gradation pattern.

and (f), the sine types gradation can lead to larger contact force. These observations are consistent with the material distributions as shown in Fig. 7, where the step change and sine type gradation occupy the largest area for hard coating and soft coating, respectively. Although both the graded shear modulus and Poisson’s ratio can influence the P ∼ 𝛿 relation, the effect of graded shear modulus is found to be more remarkable than that of graded Poisson’s ratio. These numerical results indicate that the information of material gradation may be extracted from the P ∼ 𝛿 curve.

6.2. Input parameters Since there are 14 parameters for the contact problem considered, it is practically impossible to conduct full parametric study. Thus, when investigating the effect of modulus gradation, the variables are selected as the modulus ratios 𝜇 c1 /𝜇 s1 and 𝜇 c2 /𝜇 s2 while the other parameters are fixed. The parametric combinations used to investigate the effect of modulus gradation are listed in Table 2, in which the Cases (a) and (b) are for the hard coating and Cases (c) and (d) are for soft coating. Similarly, the effect of Poisson’s ratio gradation will be explored by using the parametric combination Cases (e)–(h) given by Table 3, where both the hard and soft coating are considered. In the present analysis, the FGM coating is discretized into 100 layers to accurately approximate the continuous material gradation.

6.4. Results of contact area Figs. 10 and 11 show the dimensional interference-contact area (𝛿 ∼ 𝜋a2 ) relation for the Cases (a)–(d) and Cases (e)–(h), respectively. It is seen that the gradation of shear modulus can alter the contact area to a certain extent. For the case of hard coating Cases (a) and (b), the step change of shear modulus can generate smallest contact area 𝜋a2 under a given interference 𝛿, following by linear, exponential and sine gradation. While for the soft coating, the sine type shear modulus gradation can reduce the contact area. Unlike the effect of shear modulus gradation, the influence of graded Poisson’s ratio seems to be insignificant on the contact area. The deviations between different gradation patterns are negligible, as shown in Fig. 11.

6.3. Results of P ∼ 𝛿 relation The dimensional force-interference (P ∼ 𝛿) relations for graded shear modulus and graded Poisson’s ratio are illustrated in Figs. 8 and 9. As it is observed that for the case of hard coating Cases (a) and (b) and Cases (e) and (f), the step change of shear modulus and Poisson’s ratio can result in higher contact stiffness comparing to other types of gradation pattern. While for the case of soft coating Cases (c) and (d) and Cases (e)

6.5. Results of contact stress To assess the effect of material gradation on the contact stress, the distribution of 𝜎 zz (r,0) under the contact interference 𝛿 = 6 × 10−4 mm 68

X.W. Chen and Z.Q. Yue

Engineering Analysis with Boundary Elements 109 (2019) 57–69

for graded shear modulus and graded Poisson’s ratio are depicted in Figs. 12 and 13, respectively. As it is clearly shown that, the gradation of shear modulus can significantly change the profile of the contact stress. For the hard coating Cases (a) and (b), the continuous gradation in shear modulus can considerably reduce the peak stress at the contact center, and results in smoother stress profile comparing to the case of step change. Particularly, the location of peak stress can shift from the center to somewhere near the contact edge if the shear modulus is of exponential or sine gradation. Nonetheless, for the case of soft coating Cases (c) and (d), the continuously graded shear modulus can lead to stress concentration at the center, with the sine type gradation being the most adverse case. In contrast, the effect of graded Poisson’s ratio is less remarkable. The Poisson’s ratio gradation can only cause slight change on the magnitude of contact stress, but it cannot alter the shape of stress profile.

[2] Alinia Y, Beheshti A, Guler MA, El-Borgi S, Polycarpou AA. Sliding contact analysis of functionally graded coating/substrate system. Mechanics of Materials 2016;94:142–55. [3] Atkinson KE. The numerical solution of integral equations of the second kind. Cambridge monographs on applied and computational mathematics; 1997. [4] Aziz I. Meshless methods for multivariate highly oscillatory Fredholm integral equations. Eng Anal Bound Elem 2015;53:100–12. [5] Baker CT, 1977. The numerical treatment of integral equations. [6] Chen Z, Goltsberg R, Etsion I. Plasticity evolution in a coated sphere compressed by a rigid flat. Tribol Int 2016;98:116–24. [7] Constantinescu A, Korsunsky A, Pison O, Oueslati A. Symbolic and numerical solution of the axisymmetric indentation problem for a multilayered elastic coating. Int J Solids Struct 2013;50(18):2798–807. [8] Cundall PA, Strack OD. A discrete numerical model for granular assemblies. Geotechnique 1979;29(1):47–65. [9] Garjonis J, Kačianauskas R, Stupak E, Vadlūga V. Investigation of contact behaviour of elastic layered spheres by Fem. Mechanika 2009;3:5–12. [10] Goltsberg R, Etsion I. A universal model for the load–displacement relation in an elastic coated spherical contact. Wear 2015;322:126–32. [11] Goltsberg R, Etsion I. Contact area and maximum equivalent stress in elastic spherical contact with thin hard coating. Tribol Int 2016;93:289–96. [12] Goltsberg R, Levy S, Kligerman Y, Etsion I. Strengthening effect of soft thin coatings. J Tribol 2018;140(6):064501. [13] Greenwood J, Williamson JP. Contact of nominally flat surfaces. Proc R Soc Lond Ser A Math Phys Sci 1966;295:300–19. [14] Guler M, Erdogan F. Contact mechanics of two deformable elastic solids with graded coatings. Mech Mater 2006;38(7):633–47. [15] Guler M, Erdogan F. The frictional sliding contact problems of rigid parabolic and cylindrical stamps on graded coatings. Int J Mech Sci 2007;49(2):161–82. [16] Guler MA, Adibnazari S, Alinia Y. Tractive rolling contact mechanics of graded coatings. Int J Solids Struct 2012;49(6):929–45. [17] Gun H, Gao X-W. Analysis of frictional contact problems for functionally graded materials using BEM. Eng Anal Bound Elem 2014;38:1–7. [18] Hertz, H.R. Uber die Beruhrung fester elastischer Korper und Uber die Harte. Verhandlung des Vereins zur Beforderung des GewerbefleiBes, Berlin (1882): 449. [19] Johnson KL. Contact mechanics. Cambridge University Press; 1987. [20] Ke L-L, Wang Y-S. Two-dimensional contact mechanics of functionally graded materials with arbitrary spatial variations of material properties. Int J Solids Struct 2006;43(18–19):5779–98. [21] Korchevnik O, Goltsberg R, Kligerman Y, Etsion I. Electrical resistance model of a bilayer-coated spherical contact. IEEE Trans Comp Pack Manuf Technol 2018;8(9):1614–20. [22] Liu T-J, Wang Y-S. Axisymmetric frictionless contact problem of a functionally graded coating with exponentially varying modulus. Acta Mech 2008;199(1–4):151–65. [23] Liu T-J, Xing Y-M. Analysis of graded coatings for resistance to contact deformation and damage based on a new multi-layer model. Int J Mech Sci 2014;81:158–64. [24] Liu T-J, Xing Y-M, Wang Y-S. The axisymmetric contact problem of a coating/substrate system with a graded interfacial layer under a rigid spherical punch. Math Mech Solids 2016;21(3):383–99. [25] Nesterenko V. Dynamics of heterogeneous materials. Springer Science & Business Media; 2013. [26] Rhimi M, El-Borgi S, Lajnef N. A double receding contact axisymmetric problem between a functionally graded layer and a homogeneous substrate. Mech Mater 2011;43(12):787–98. [27] Sneddon IN. The Use of Integral Transform. New York: McGram-Hill Company; 1972. [28] Suresh S. Graded materials for resistance to contact deformation and damage. Science 2001;292(5526):2447–51. [29] Yan J, Mi C. On the receding contact between an inhomogeneously coated elastic layer and a homogeneous half-plane. Mech Mater 2017;112:18–27. [30] Yue ZQ. On generalized Kelvin solutions in a multilayered elastic medium. J Elast 1995;40(1):1–43. [31] Yue ZQ. On elastostatics of multilayered solids subjected to general surface traction. Quart J Mech Appl Math 1996;49(3):471–99. [32] Yue ZQ. Yue’s solution of classical elasticity in n-layered solids: part 1, mathematical formulation. Front Struct Civil Eng 2015;9:215–49. [33] Yue ZQ. Yue’s solution of classical elasticity in n-layered solids: part 2, mathematical verification. Frontiers of Structural and Civil Engineering 2015;9:250–85. [34] Zaheer-ud-Din. Meshless methods for one-dimensional oscillatory Fredholm integral equations. Appl Math Comput 2018;324:156–73. [35] Zaheer-ud-Din. Meshless methods for two-dimensional oscillatory Fredholm integral equations. J Comput Appl Math 2018;335:33–50. [36] Zaman S. On numerical evaluation of integrals involving oscillatory Bessel and Hankel functions. Numer Algorithms 2019:1–19.

7. Concluding remarks Based on the contact model developed and the numerical results presented in this study, the following conclusions can be made regarding the nonlinear contact mechanics for two elastically dissimilar spheres with FGM coatings: (1) Any possible gradation of material properties, including shear modulus and Poisson’s ratio, in the FGM coating can be accurately and conveniently approximated by using the present multilayer half space model. (2) The present model is general in that it can analytically degenerated to classical Hertz’s solution for uncoated spheres, and it also covers the simple case of spheres with homogeneous coatings. (3) The graded shear modulus in FGM coating has significant effect on the contact force-interference (P ∼ 𝛿) relation, particularly for the case of soft coating. While the effect of graded Poisson’s ratio is relatively insignificant. (4) All other parameters being equal, the contact area 𝜋a2 can be affected to a certain degree by the graded shear modulus, but it is nearly immune to the graded Poisson’s ratio. (5) In comparison to the conventional homogeneous coating, the FGM coating with continuously graded shear modulus can effectively redistribute the contact stress and the role of graded Poisson’s ratio is also important. The results presented in this paper may be valuable references to the design of structural components such as roller bearing and connecting rods reinforced with FGM coating. With the present contact model, one can also easily explore the complicated relations between the contact quantities, e.g. contact force, area and stress, for other types of material gradation, and subsequently apply these understandings to the wear analysis of coated material. Acknowledgment The work described in this paper was supported by the Research Grants Council of Hong Kong SAR Government (GRF Nos. 17204415 and 17207518). The authors thank the 2 reviewers for their valuable comments and suggestions. The first author Mr. Chen XW also thanks The University of Hong Kong for scholarship. References [1] Alinia Y, Asiaee A, Hosseini-nasab M. Stress analysis in rolling contact problem of a finite thickness FGM layer. Meccanica 2019;54(1–2):183–203.

69