Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

56, 1, pp. 147-156, Warsaw 2018
DOI: 10.15632/jtam-pl.56.1.147

ON AXISYMMETRIC HEAT CONDUCTION PROBLEM FOR MULTILAYER

GRADED COATED HALF-SPACE

Dariusz M. Perkowski, Roman Kulchytsky-Zhyhailo,Waldemar Kołodziejczyk

Bialystok University of Technology, Faculty of Mechanical Engineering, Białystok, Poland

e-mail: perkowski.d.m@gmail.com; r.kulczycki@pb.edu.pl; w.kolodziejczyk@pb.edu.pl

The technique of integral Hankel transform to find the solution of heat conduction in half-
-space coated by a multilayered package of homogenous laminae is applied. The half-space
is heated by the given heat flux on the boundary surface. The temperature and heat flux
distribution in the radial direction is analyzed for two types of coatings: 1) when the heat
conductivity coefficient is described by a power or exponential function of the distance to
the boundary surface; 2) multilayered coating has a periodic structure.

Keywords: temperature, heat flux, FGM, composite, homogenizedmodel

1. Introduction

Modern engineering construction require the use ofmaterials with appropriate thermal andme-
chanical properties whichmake themmore andmore complicated structures. To suchmaterials
belongmedia with functionally changing gradation properties (i.e. gradientmaterials). Further-
more, the surface coverage of gradient materials significantly affect the behaviour of the bodies
under the influence of mechanical and thermal loads.

The formulations of the problems lead to boundary value problems of partial differential
equations with varying coefficients. For the power (or exponential) law of variation of the heat
conduction coefficient (or Young’s modulus), the analytical methods of solutions are known
(Guler andErdogan, 2004, 2006; LiuandWang, 2008;Matysiak et al., 2011;Kulchytsky-Zhyhailo
andBajkowski, 2015). If the thermal andmechanical properties are describedbyother functions,
obtaininganalytical solutions encounters considerablemathematical difficulties.Parallelwith the
applicationof analyticalmethods for the solutionof partial differential equations, inhomogeneous
layers are also modeled by using an approach according to which the coating is replaced by
a package of homogeneous or inhomogeneous layers (Ke and Wang, 2006, 2007; Kulchytsky-
-Zhyhailo and Bajkowski, 2012, 2015; Liu T.-J. et al., 2008; Liu and Wang, 2009; Liu J. et al.,
2011, 2012).Most studies are basedon two-dimensional problemsof elasticity or thermoelasticity
(Barik et al., 2008; Choi and Paulino, 2008; Diao, 1999; Diao et al., 1994; Liu J. et al., 2011,
2012). However, the axisymmetric and three-dimensional problems are dealt with in a much
lesser degree.

A special type of graded coating is themultilayered coatingwithperiodic structure (Farhat et
al., 1997;Vevodin et al., 2001). In themodelingof laminatedhalf-spaces or coatingswithperiodic
structures, it is customary to use two different approaches. In the first of these approaches the
layers are considered as separate continuousmedia.The secondapproach is based on the analysis
of a homogenized uniform coating whose properties are determined on the basis of thematerial
properties and geometric characteristics of the strip of periodicity (Matysiak and Woźniak,
1987; Woźniak, 1987). The solution obtained for the laminated half-space is compared with
Kulchytsky-Zhyhailo andMatysiak (2005).



148 D. Perkowski et al.

In the presentwork,we consider an axisymmetric problemof heat conductivity for half-space
withamultilayered coating.Theboundarysurface isheatedby thegivenheatflux. Inparallel,we
obtain: 1) analytical solutions of the problems for the coatingwhoseheat conductivity coefficient
is described by a continuous function of the distance to the boundary surface; 2) the solution for
themultilayered coating with a periodic structure which is described by a homogenized uniform
layer (Matysiak and Woźniak, 1987; Woźniak, 1987). We analyze the difference between the
temperature and heat flux in the non-homogeneous half-space caused by the use of two different
models of nonuniform coatings. The obtaining of the smallest deviations in consideration of the
heat conduction problem will be a strong argument in favor for application of the proposed
methods for solving axisymmetric and three-dimensional problems of thermoelasticity for the
functionally graded coated half-space.

2. Formulation of the problem

Assume that the surface z=h of the non-homogeneous half-space is heated by a heat flux q(r)
on the circle of radius a, where r, z are dimensionless cylindrical coordinates referred to the
linear size a, h=H/a,H is thickness of the coating (Fig. 1).

Fig. 1. The scheme of the body

The non-homogeneous half space is formed by the homogeneous half-space with the heat
conductivity coefficientK0 and a systemof non-homogeneous layers with thicknessesHi and the
heat conductivity coefficientsKi, i=1,2, . . . ,n, respectively, where the value of the parametern
corresponds to the number of the layer in the package. Assume that the conditions of perfect
thermal contact are realized between the layers of the coating and between the coating and the
base.
The analyzed problemof theory of heat conduction is reduced to the solution of the following

partial differential equations

∂2Ti
∂r2
+
1

r

∂Ti
∂r
+
1

Ki

∂

∂z

(
Ki
∂Ti
∂z

)
=0 i=0,1, . . . ,n (2.1)

with boundary conditions imposed on the surface of the non-homogeneous half-space

∂Tn
∂z
=
q(r)a

Kn
H(1−r) z=h (2.2)

conditions of perfect thermal contact between the components of the considered half-space

Ti =Ti+1 Ki(h
∗

i)
∂Ti
∂z
=Ki+1(h

∗

i)
∂Ti+1
∂z

z=h∗i i=0, . . . ,n−1
(2.3)

and conditions imposed at infinity

Ti → 0 r2+z2 →∞ i=0,1, . . . ,n (2.4)



On axisymmetric heat conduction problem for multilayer... 149

where Ti is the temperature in the i-th component of the non-homogenous medium, the index
i = 0 describes the parameters and functions of state in the homogeneous half-space, h∗i is
the coordinate z of the upper surface of the i-th component of the non-homogenous half-space,
h∗0 =0, h

∗
i =h

∗
i−1+hi, hi =Hi/a, i=1, . . . ,n, h

∗
n =h,H(r) – Heaviside step function.

3. Method of solution

The solution of the boundary value problem is sought by applying theHankel integral transfor-
mation (see Sneddon, 1966)

T̃i(s,z)=

∞∫

0

Ti(r,z)rJ0(sr) dr (3.1)

where J0(sr) is the Bessel function. The solution to equation (2.1) was determined by Hankel
integral transformation technique (3.1). The temperature for the homogeneous half-space in the
Hankel transform space which satisfies the regularity conditions at infinity (2.4) can be written
in the form

T̃0(s,z)= t0(s)exp(sz) (3.2)

where t0(s) is the unknown function.
We considered the following cases.

Case A

Let n=1. The dependence of the heat conductivity coefficient on the coordinate z is described
by the formula

K1(z)=K0exp(αz) α=
1

h
ln
(KS
K0

)
0¬ z¬h (3.3)

where KS is the heat conductivity coefficient on the surface of the inhomogeneous half-space.
The general solution to differential equation (2.1) in theHankel transform space specified in the
coating can be written in the form

T̃1(s,z)= t1(s)exp(α
(−)z)+ t2(s)exp(α

(+)z) (3.4)

where 2α(±) =−α±
√
α2+4s2, t1(s) and t2(s) are the unknown functions.

Case A’

Let n=1. The dependence of the heat conductivity coefficient on the coordinate z is described
by the power function

K1(z)=K
∗(c±z)α c=±h

((KS
K0

)1/α
−1

)−1

K∗ =
K0
cα

0¬ z¬h

(3.5)

In equation (3.5) for the caseK0 <KS, we take sign “+”, whenK0 >KS – sign “−”.
The Hankel transform of temperature in the coating can be written in the form

T̃1(s,z)= t1(s)ζ
pIp(sζ)+ t2(s)ζ

pKp(sζ) (3.6)

where 2p+α=1, ζ = c±z, Ip(sζ),Kp(sζ) are modified Bessel functions.



150 D. Perkowski et al.

Moreover, if the analytical solution of partial differential equation with variable coefficient
(2.1) is not known, the non-homogeneous coating can be replaced by a multilayered system of
homogeneous layers. Their thermal properties are described by their heat conductivity coeffi-
cients

Ki =
1

hi

h∗
i∫

h∗
i−1

K1(z) dz (3.7)

Case B

Thecoating composedofnhomogeneous layers.Thegeneral solution toequations (2.1) expressed
in the Hankel transform domain takes the form

T̃i = t2i−1(s)sinh[s(h
∗

i −z)]+ t2i(s)cosh[s(h
∗

i −z)] i=1,2, . . . ,n (3.8)

where ti(s), i=1, . . . ,2n are the unknown functions.

Case C

Themultilayered coatingwith amicroperiodical structure.Assume that the repeated fundamen-
tal layer comprises twohomogeneous elastic sublayerswithdifferent thicknesses (hI andhII) and
thermal conductivities (KI andKII). A large number of equations and boundary conditions on
the interfaces complicates the solution of the problem.Another approach is using a homogenized
model (Choi andPaulino, 2008; Diao, 1999) in which properties of the homogenized coating are
determined on the base of properties of the components.
Applying the homogenized model to the coating, we solve the boundary value problem de-

scribed by the equation (Matysiak andWoźniak, 1987; Woźniak, 1987)

∂2T1
∂r2
+
1

r

∂T1
∂r
+
1

p21

∂2T1
∂z2
=0 (3.9)

where T1 is the temperature in the homogenized coating

p21 = K̃K
−1
c Kc =

KIKII
(1−η)KI +ηKII

K̃ = ηKI +(1−η)KII η=
hI

hI +hII

(3.10)

the boundary condition imposed on the surface of the non-homogeneous half-space

∂T1
∂z
=
q(r)a

Kc
H(1−r) z=h (3.11)

and boundary conditions of perfect thermal contact between the homogenized coating and the
substrate

T0 =T1 K0
∂T0
∂z
=Kc

∂T1
∂z

z=0 (3.12)

Boundary conditions (2.4) stay without change. The general solution to equation (3.9) in the
Hankel transform takes the form

T̃1 = t1(s)sinh[(h−z)sp1]+ t2(s)cosh[(h−z)sp1] (3.13)



On axisymmetric heat conduction problem for multilayer... 151

Equations (3.2), (3.4), (3.6), (3.8) and (3.13) contain the unknown functions ti(s). These
functions are obtained satisfying boundary conditions (2.2) and (2.3) (or (3.11)) and (3.12) in
Case C. Satisfying the boundary conditions, the functions ti(s) may be written as

ti(s)=
q̃(s)a

K̂s
t̂i(s) (3.14)

where the functions t̂i(s) are obtained from the solution to linear equations (see Appendix A),
q̃(s) is the Hankel transform of heat flux, K̂ =Kn in Case A, A’, B and K̂ =Kcp1 in Case C.

Applying the inverse Hankel transform to equations (3.2), (3.4), (3.6), (3.8) and (3.13),
temperature can be find at the desired location (see Sneddon, 1966)

Ti(r,z) =

∞∫

0

sT̃i(s,z)J0(sr) ds (3.15)

At internal points of the non-homogeneous half space (z <h) the integrals are evaluated nume-
rically using the Gaussian quadrature (Abramowitz and Stegun, 1964). On the surface z = h,
we take into account the asymptotic behaviour of the functions t2n−1(s) and t2n(s) as s→∞.
The continuity of the results when z→h has been also verified.

4. Numerical results and discussion

Assume that the heat flux is elliptically distributed as follows

q(r)=

{
Q0
√
1−r2 for r < 1

0 for r­ 1
(4.1)

where r = r∗/a = 1 is radius of the circle heat zone and the heat flux (4.1) in the Hankel
transform space takes the form (Gradshteyn and Ryzhik, 2015)

q̃(s)=

√
π

2
Q0
J3/2(s)

s3/2
(4.2)

where J3/2(s) is the Bessel function.

Theanalysis of the original relations inCaseA (orA’) shows that the solution of the problem
of modeling of the inhomogeneous coating by the package of homogeneous layers depends on
three (Case A) or four (Case A’) dimensionless parameters: thickness of the coating h, ratio of
heat conductivity coefficients on the surfaces of thenon-homogeneous half space and the substra-
teKS/K0, parameterα (only for Case A’) and the number of layers in the package n. A similar
solution obtained for an inhomogeneous coating with regard to the continuous dependence of
thermal properties on the coordinate is independent of the parameter n. In what follows, we
assume that h=0.4,K0/KS =2, 4, or 8, α=1 (only for Case A’), and n=10, 20, 40, or 80.

The temperature and the heat flux in the radial direction r on the considered nonhomo-
geneous surface are shown in Figs. 2 and 3 ((a) – Case A, (b) – Case A’). In this figures, the
continuous lines correspond to the solution of the problem with continuous variation of the
thermal properties. The rhombi correspond to the results obtained for the package formed by
40 homogeneous layers. The results of calculations presented in Fig. 2 show good agreement
between the solutions obtained using the analyzed two models of the coating. As follows from
Fig. 3, the maximum absolute error of calculation of the radial heat flux component can be
observed on the end of the heated zone for r=1.



152 D. Perkowski et al.

Fig. 2. The dimensionless temperature distribution on the surface z=h: (a) – Case A, (b) – Case A’;
1 –K0/KS =2, 2 –K0/KS =4, 3 –K0/KS =8; h=0.4

Fig. 3. The dimensionless radial heat flux on the surface z=h: (a) – Case A, (b) – Case A’;
1 –K0/KS =2, 2 –K0/KS =4, 3 –K0/KS =8; h=0.4

The values of the radial heat flux in this point are presented in Table 1. The analytical
solution is presented in the rows of Table 1 with n→∞. The relative error of their evaluation
with the help ofmodeling of the inhomogeneous coating by the package ofn homogeneous layers
is given in columns with n= 10, 20, 40, and 80. It is easy to see that as the number of layers
becomes twice larger, the corresponding error becomes almost twice lower. In the case where
there are 80 layers in the package andK0/KS ¬ 8, the error of finding the heat flux at the point
r=1, z=h does not exceed 2.2%.

Estimating the original relations, we conclude that the distributions of temperature andheat
flux in the problem of homogenized coating (Case C) depend on four dimensionless parameters:
thickness of the coating h, ratios of heat conductivity coefficients KI/K0 and KII/K0 and the
ratio of the thicknesses of layers in the strip of periodicity hI/hII. Similar distributions for the
non-uniform coating additionally depend on the number of layers in the coating n. To decrease
the number of input parameters, we assume that the thermal properties of one layer in the strip
of periodicity coincide with the thermal properties of the base (KI/K0 = 1 or KII/K0 = 1)
and the thicknesses of all layers in the stack are identical (hI/hII = 1). We also assume that
K0/KI(or K0/KII)= 4, h=0.2, 0.4, or 0.8, and n=10, 20, 40, or 80.



On axisymmetric heat conduction problem for multilayer... 153

Table 1. The dimensionless radial heat flux at the point r = 1, z = h (n → ∞) and the
errors of their evaluation as a result of modeling of the inhomogeneous coating by the package
of n homogeneous layers n=10, 20, 40, and 80

K0/KS n=10 n=20 n=40 n=80 n→∞
Case A

2 2.79% 1.43% 0.68% 0.29% 0.6169

4 5.40% 2.73% 1.29% 0.54% 0.5021

8 7.79% 3.90% 1.82% 0.71% 0.4236

Case A’

2 3.95% 2.03% 0.98% 0.43% 0.6022

4 11.00% 5.68% 2.72% 1.14% 0.4562

8 23.00% 11.97% 5.65% 2.18% 0.3435

The dimensionless temperature at the centre of the heating area for different thicknesses
of the coating and different numbers of layers is presented in Table 2. The temperature cal-
culated for the homogenized model is in the last column. The relative differences between the
temperature in the non-homogeneous coating and the temperature in the homogenized coating
are presented in columns with n = 10, 20, 40, and 80. The upper numbers in cells have been
calculated for the case K0/KI =4, K0/KII =1, the lower numbers were obtained for the case
K0/KI = 1, K0/KII = 4. It can be seen that as the number of layers becomes twice larger,
the errors become twice lower then. For the same value of the parameter δ=h/n (for example:
h = 0.2, n= 10; h = 0.4, n= 20; and h = 0.8, n = 40), these errors are in the same order of
magnitude.

Table 2. The dimensionless temperature at the centre of the heating zone at the point r = 0,
z=h (n→∞) and the errors of their evaluation as a result of modeling of the inhomogeneous
coating by the package of n homogeneous layers n=10, 20, 40, and 80

h n=10 n=20 n=40 n=80 n→∞

0.2
1.12% 0.56% 0.28% 0.14%

1.0841
−1.13% −0.56% −0.28% −0.14%

0.4
3.35% 1.67% 0.83% 0.41%

1.2481
−3.26% −1.64% −0.83% −0.42%

0.8
7.99% 3.92% 1.94% 0.96%

1.3887
−7.31% −3.76% −1.90% −0.96%

Figures 4a and 4b show the dimensionless radial heat flux as functions of z for r = 1 and
two numbers of layers (n = 20 and n = 40). In Fig. 4, the rhombi mark the numerical results
obtained for the non-homogeneous laminated coating, whereas the solid lines correspond to the
homogenized coating. It should be emphasized that in the case of homogenized coating, we do
not know which layer of the slip of periodicity is located at the analyzed point of the coating.
Hence, the radial heat flux at every point of the coating is described by the two curves. Curves 1
and 2 correspond to the heat flux acting in the layers with smaller and larger heat conductivity
coefficients, respectively. In the homogeneous substrate, curves 1 and 2 coincide.
Comparing the heat flux obtained in both analyzed problems, we conclude that only in the

case of the heat flux acting in the homogeneous substrate we get deviations comparable with
the deviations of temperature. In the layers of the coating, the deviations of heat flux vary from
1%-5% (K0/KI ¬ 4, n = 20) up to 10%-20% for the heat flux acting on the boundary of the
region of heating. The indicated deviations strongly depend on the gradient of the analyzed
parameter in the investigated layer of the slip of periodicity, which explains the following obse-



154 D. Perkowski et al.

Fig. 4. The dimensionless radial heat flux on the line r=1: (a) – n=20, (b) – n=40;K0/KI =4,
h=0.4, hI/hII =1

rvations: in the layers with lower heat conductivity coefficients, the deviations aremuch smaller
(Fig. 4) and the maximum deviations are observed at the point (1,h). As could be expected,
the agreement between the solutions improves as the number of layers in the coating increases.

5. Conclusions

This paper provides the solution to the problem of the inhomogeneous half-space heated by the
heat flux. It is shown that the solution to the problem for a package of 20-80 homogeneous layers
is in good agreement with the analytical solution to the problem for the coating whose depen-
dence of the heat conductivity coefficient on the coordinate z is described by an exponential (or
power) function.This is a strong argument for the possibility ofmodeling of the gradient coating
with continuous variation of thermal and mechanical properties by a package of homogeneous
layers.
It is shown that the solution of the axisymmetric problem of heat conductivity for the half-

-spacewitha laminated coating of theperiodic structureheatedbyheatflux is ingoodagreement
with the solution of the problem inwhich the coating ismodeled by a homogenized coating. The
smallest deviations are obtainedwhile finding the temperature andheat flux in the homogeneous
substrate.

A. Appendix

A system of linear equations for determination of the functions t̂i(s), i=0,1,2:
—Case A

− t̂0(s)+ t̂1(s)+ t̂2(s)= 0

− t̂0(s)+α(−)s−1t̂1(s)+α(+)s−1t̂2(s)= 0

α(−)s−1t̂1(s)exp(α
(−)h)+α(+)s−1t̂2(s)exp(α

(+)h)= 1

(A.1)

—Case A’

− t̂0(s)+ t̂1(s)cpIp(sc)+ t̂2(s)cpKp(sc)= 0
t̂0(s)∓ t̂1(s)cpIp−1(sc)± t̂2(s)cpKp−1(sc)= 0
t̂1(s)(c±h)pIp−1

(
s(c±h)

)
− t̂2(s)(c±h)pKp−1

(
s(c±h)

)
=±1

(A.2)



On axisymmetric heat conduction problem for multilayer... 155

—Case C

− t̂0(s)+ t̂1(s)sinh(sp1h)+ t̂2(s)cosh(sp1h)= 0

K0(Kcp1)
−1t̂0(s)+ t̂1(s)cosh(sp1h

)
+ t̂2(s)sinh(sp1h)= 0

t̂1(s)=−1

(A.3)

A system of linear equations for determination of the functions t̂i(s), i=0,1, . . . ,2n in Case B

− t̂0(s)+ t̂1(s)sinh(sh1)+ t̂2(s)cosh(sh1)= 0
K0K

−1
1 t̂0(s)+ t̂1(s)cosh(sh1)+ t̂2(s)sinh(sh1)= 0

− t̂2i(s)+ t̂2i+1(s)sinh(shi+1)+ t̂2i+2(s)cosh(shi+1)= 0 i=1, . . . ,n−1
−KiK−1i+1t̂2i−1(s)+ t̂2i+1(s)cosh(shi+1)+ t̂2i+2(s)sin(shi+1)= 0 i=1, . . . ,n−1
t̂2n−1(s)=−1

(A.4)

Acknowledgments

Thisworkwascarriedoutwithin theproject“selectedproblemsof thermomechanics formaterialswith

temperature dependent properties”. The project was financed by the National Science Centre awarded

based on the decision number DEC-2013/11/D/ST8/03428.

References

1. Abramowitz M., Stegun I.A., 1964, Handbook of Mathematical Functions: with Formulas,
Graphs, and Mathematical Tables, United States Department of Commerce, National Bureau of
Standards

2. Barik S.P., Kanoria M., Chaudhuri P.K., 2008, Steady state thermoelastic contact problem
in a functionally gradedmaterial, International Journal of Engineering Science, 46, 775-789

3. Choi H.J., Paulino G.H., 2008, Thermoelastic contact mechanics for a flat punch sliding over
a graded coating/substrate system with frictional heat generation, Journal of the Mechanics and
Physics of Solids, 56, 4, 1673-1692

4. DiaoD.F., 1999,Finite element analysis on local yieldmapandcriticalmaximumcontactpressure
for yielding in hard coating with an interlayer under sliding contact, Tribology International, 32,
25-32

5. Diao D.F., Kato K., Hayashi K., 1994, The maximum tensile stress on a hard coating under
sliding friction,Tribology International, 27, 4, 267-272

6. Farhat Z.N., Ding Y., Northwood D.O., Aplas A.T., 1997, Nanoindentation and friction
studies on Ti-based nanolaminated films, Surface and Coatings Technology, 89, 24-30

7. Gradshteyn I.S., Ryzhik I.M., 2015,Table of Integrals, Series and Products, Academic Press

8. Guler M.A., Erdogan F., 2004, Contact mechanics of graded coatings, International Journal
of Solids and Structures, 41, 3865-3889

9. GulerM.A., ErdoganF., 2006,Contactmechanics of two deformable elastic solidswith graded
coatings,Mechanics of Materials, 38, 633-647

10. Ke L.-L., Wang Y.-S., 2006, Two-dimensional contact mechanics of functionally graded mate-
rials with arbitrary spatial variations of material properties, International Journal of Solids and
Structures, 43, 18/19, 5779-5798

11. Ke L.-L., Wang Y.-S., 2007, Two-dimensional sliding frictional contact of functionally graded
materials,European Journal of Mechanics – A/Solids, 26, 171-188



156 D. Perkowski et al.

12. Kulchytsky-ZhyhailoR., Bajkowski A., 2012,Analytical and numericalmethods of solution
of three-dimensional problem of elasticity for functionally graded coated half-space, International
Journal of Mechanical Sciences, 54, 105-112

13. Kulchytsky-Zhyhailo R., Bajkowski A., 2015, Three-dimensional analytical elasticity solu-
tion for loaded functionally graded coated half-space, Mechanics Research Communications, 65,
43-50

14. Kulchytsky-ZhyhailoR.,Matysiak S.J., 2005,On some heat conduction problem in a perio-
dically two-layered body. Comparative results, International Communications in Heat and Mass
Transfer, 32, 3-4, 332-340

15. Liu J., Ke L.-L.,WangY.-S., 2011, Two-dimensional thermoelastic contact problem of functio-
nally gradedmaterials involving frictional heating, International Journal of Solids and Structures,
48, 18, 2536-2548

16. Liu J., Ke L.-L., Wang Y.-S., Yang J., Alam F., 2012, Thermoelastic frictional contact of
functionally gradedmaterials with arbitrarily varying properties, International Journal of Mecha-
nical Sciences, 63, 1, 86-98

17. Liu T.-J., Wang Y.-S., 2008,Axisymmetric frictionless contact problem of a functionally graded
coating with exponentially varyingmodulus,Acta Mechanica, 199, 151-165

18. Liu T.-J., Wang Y.-S., 2009, Reissner-Sagoci problem for functionally graded materials with
arbitrary spatial variation of material properties, Mechanics Research Communications, 36,
322-329

19. LiuT.-J.,WangY.-S., ZhangC., 2008,Axisymmetric frictionless contact of functionally graded
materials,Archive of Applied Mechanics, 78, 267-282

20. Matysiak S., Kulchytsky-Zhyhailo R., Perkowski D., 2011, Reissner-Sagoci problem for
a homogeneous coating on a functionally graded half-space,Mechanics Research Communications,
38, 4, 320-325

21. Matysiak S.J., Woźniak C., 1987, On the modelling of heat conduction problem in laminated
bodies,Acta Mechanica, 65, 223-238

22. Sneddon I.N., 1966 The Reissner-Sagoci problem,Proceedings of the Glasgow Mathematical As-
sociation, 7, 136-144

23. Voevodin A.A., Iarve E.V., Ragland W., Zabinski J.S., Donaldson S., 2001, Stress ana-
lyses and in-situ fracture observation of wear protective multilayer coatings in contact loading,
Surface and Coatings Technology, 148, 38-45

24. Woźniak C., 1987, A nonstandard method of modelling of thermoelastic periodic composites,
International Journal of Engineering Science, 25, 483-499

Manuscript received January 31, 2017; accepted for print July 31, 2017