# Modeling of Thermal Distributions around a Barrier at the Interface of Coating and Substrate

Authors: Ali Sahin
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Due to constant heat flux, the thermal distribution around an insulated barrier at the interface of substrate and functionally graded material (FGM) which are essentially two-phase particulate composites is examined  in such a way that the volume fractions of the constituents vary continuously in the thickness direction. Using integral transform method, two-dimensional steady-state diffusion equation with variable conductivity is turned into constant coefficient differential equation. Reducing that equation to a singular integral equation with Cauchy type, the temperature distribution around the barrier is obtained by defining an unknown function, which is called density function, as a series expansion of orthogonal polynomials. Results are shown for different thickness and nonhomogeneity parameters of FGM.


---

## Body

## 1. Introduction

There are many engineering applications under severe thermal loading that require high temperature resistant materials in various forms of composites and bonded materials such as power generation, transportation, aerospace, and thermal barrier coatings. New developments in science and technology rely on the developments of new materials. Composites appear to provide the necessary flexibility in the design of these new materials, which are essentials that every part of the material in use exhibits uniform properties. In the mid-1980s, a new composite material, which was initially designed as a kind of thermal barrier coating used in aerospace structural applications and fusion reactors, was found by a group of scientists in Japan. Because of the material’s structure, it is named as functionally graded material (FGM). FGMs were used in modern technologies as advanced structures where the composition or the microstructure is locally varied so that a certain variation of the local material properties is achieved [1]. FGMs are also developed for general use as structural components in extremely high-temperature environments. The concept is to make a composite material by varying the microstructure from one material to another material with a specific gradient. The transition between the two materials can usually be approximated by means of a power series or an exponential function [2–5].The aircraft and aerospace industry and the computer circuit industry are very interested in the possibility of materials that can withstand very high thermal gradients. This is normally achieved by using a ceramic layer connected with a metallic layer. The composition profile which varies from 0% ceramic at the interface to 100% ceramic near the surface, in turn, is selected in such a way that the resulting nonhomogeneous material exhibits the desired thermomechanical properties. The concept of FGMs could provide great flexibility in material design by controlling both composition profile and microstructure [6].A number of reviews dealing with various aspects of FGMs have been published in the past few decades. They show that most of early research studies in FGMs had focused more on thermal stress analysis and fracture mechanics. Fracture mechanics of FGMs have been studied analytically by Erdogan and coworkers [7–9]. Erdogan identified a number of typical problem areas relating to the fracture of FGMs by considering mainly the investigation of the nature of stress singularity near the tip of a crack in a different geometry [10]. Erdogan also investigated the nature of the crack-tip stress field in a nonhomogeneous medium having a shear modulus with a discontinuous derivative. The problem was considered for the simplest possible loading and geometry, namely, the antiplane shear loading of two bonded half spaces, in which the crack is perpendicular to the interface. It was shown that the square-root singularity of the crack-tip stress field is unaffected by the discontinuity in the derivative of the shear modulus [11]. Related to the fracture problems in composite materials, the solution of integral equations with strongly singular kernels is examined by Kaya and Erdogan [12, 13]. In an axisymmetric coordinate system, an embedded axisymmetric crack in a nonhomogeneous infinite medium was studied by Ozturk and Erdogan [14]. They showed the effect of the material nonhomogeneity on the stress intensity factors under constant Poisson’s ratio.Due to the material mismatch at the interface of substrate and coating, the thermal distributions and thermal stresses on the crack or insulated barrier at the interface are examined by researchers. A general analysis of one-dimensional steady-state thermal stresses in a hollow thick cylinder made of functionally graded material is developed by Jabbari et al. [15]. The material properties were assumed to be nonlinear with a power law distribution. The mechanical and thermal stresses were obtained through the direct method of solution of the Navier equation. The problem of general solution for the mechanical and thermal stresses in a short length functionally graded hollow cylinder due to the two-dimensional axisymmetric steady-state loads was solved using the Bessel functions by Jabbari et al. [16]. A standard method was used to solve a nonhomogeneous system of partial differential Navier equations with nonconstants coefficients, using Fourier series.Jin and Noda [17] examined an internal crack problem in nonhomogeneous half-plane under thermal loading using the airy stress function method and Fourier transform. They reduced the problem to a system of singular integral equations and solved it by numerical methods. They used superposition method by defining the problem in two cases. One is the linear one-dimensional heat conduction problem under constant heat flux without crack, and the other one is the two-dimensional heat conduction problem with an insulated crack subject to constant heat flux which is acted in the opposite direction. It was also considered the problem of an axisymmetric penny-shaped crack embedded in an isotropic graded coating bonded to a semi-infinite homogeneous medium by Rekik et al. [18]. The coating’s material gradient is parallel to the axisymmetric direction and is orthogonal to the crack plane. They used Hankel transform to convert the equations into coupled singular integral equations along with the density function, and they solved it numerically.In this study, the Hankel integral transform method will be used to solve the heat equation in axisymmetric coordinate system. Problem will be examined as a one-dimensional and a two-dimensional heat conduction problem that is a mixed boundary value problem over the real line. Using mixed boundary conditions a Fredholm integral equation will be obtained with Cauchy type singularity and then it will be solved by using some known numerical techniques [19, 20].


## 2. Definition of the Problem

The thermal distribution around a penny-shaped barrier at the interface of graded composite coating and a substrate is given by the following steady-state heat equation in axisymmetric coordinate system:(1)1r∂∂r(rk(z)∂T∂r)+∂∂z(k(z)∂T∂z)=0,
and the conductivities of the substrate and the graded composite coating are, respectively, given as
(2)k(z)=k0,z≤0,k(z)=k0eδz,z≥0,
where k0 is a constant and δ is the nonhomogeneity parameter related to the graded coating. Note that the conductivity is continuous at the interface of substrate and the graded composite coating. As shown in Figure 1, it is considered a penny-shaped barrier with radius a, centered at the origin of the axisymmetric coordinate system. It is assumed that a uniform heat flux is applied over the stress free boundary, and the barrier faces remain insulated.Figure 1
Geometry of the heat conduction problem.The solution can be obtained using superposition method which is an addition of one- and two-dimensional heat conduction problems,T1(z) and T2(r,z), respectively, as shown in Figures 2(a) and 2(b). As in Figure 2(a), it will be assumed that there will be no barrier and in flux causes thermal distribution only z direction. On the other hand, in Figure 2(b), flux will be assumed in an opposite direction on the barrier that causes thermal distribution in (r,z) plane.(a) One-dimensional heat conduction without an insulated barrier, (b) Two-dimensional heat conduction with an insulated barrier.
(a)
(b)Rewriting (1) assuming that no changing in r direction:
(3)δdT1dz+d2T1dz2=0,0≤z≤h;d2T1dz2=0,z≤0,
along with suitable boundary conditions:
(4)k(z)ddzT1(z)⟶finite,z⟶-∞,k(z)ddzT1(z)=-Q0,z=h,T1(0+)=T1(0-),
the solution of the one-dimensional heat conduction is obtained straightforward as
(5)T1(z)={Q0k0δe-δz,0≤z≤h,Q0k0δ,z≤0,
where the continuity can be seen as z→0. For two-dimensional heat conduction problem, (1) can be simplified for 0<z≤h as
(6)∂2∂r2T2(r,z)+1r∂∂rT2(r,z)+∂2∂z2T2(r,z)+δ∂∂zT2(r,z)=0,
and for z<0,
(7)∂2∂r2T2(r,z)+1r∂∂rT2(r,z)+∂2∂z2T2(r,z)=0
with standard boundary conditions:
(8)k(z)∂∂zT2(r,z)⟶finite,z⟶-∞,k(z)∂∂zT2(r,z)=0,z=h,0≤r<∞,∂∂zT2(r,0+)=∂∂zT2(r,0-),a<r<∞,
and mixed boundary conditions:
(9)k(z)∂∂zT2(r,z)=Q0,z⟶0,0≤r<a,(10)T2(r,0+)=T2(r,0-),a<r<∞.
Equation (6) will be solved using Hankel integral transform such that T2(r,z) denote the Hankel transform of zero order and τ(ρ,z) denote inverse Hankel transform of zero order [21] shown as below, respectively,
(11)T2(r,z)=∫0∞τ(ρ,z)J0(rρ)ρdρ,τ(ρ,z)=∫0∞T2(r,z)J0(rρ)rdr.
Using Hankel transform, the solution of (6) along with boundary conditions (8) can be obtained as
(12)τ(ρ,z)={A(ρ)e-(m2+m1)z+A(ρ)m2+m1m2-m1×e-2m2h+(m2-m1)z,0<z≤h,A(ρ)1ρ(e-2m2h-1)(m2-m1)×eρz,-∞<z<0,
where m1=δ/2, m2=ρ(δ/2ρ)2+1 and observing that m2+m1>0 and m2-m1>0.


## 3. Evaluation of Integral Equation

The unknown valueA(ρ) can be obtained by defining a new function, which is called density function [11], such as
(13)ψ(r)=∂∂r(T2(r,0+)-T2(r,0-)),
where ψ(r) satisfies the following conditions:
(14)∫0aψ(r)dr=0,ψ(r)=0,ata<r<∞,ψ(r)=-ψ(-r).
Substituting (12) into (13) along with the conditions given in (14), the unknown function A(ρ) can be obtained:
(15)A(ρ)=1F(ρ)∫0aψ(s)J1(sρ)sds,
where
(16)F(ρ)=-ρ-ρ(m2+m1m2-m1)e-2m2h+(e-2m2h-1)(m2-m1).
Using the boundary condition (9) in transformed domain as z→0-, it can be obtained as
(17)∫0∞(e-2m2h-1)(m2-m1)A(ρ)J0(rρ)ρdρ=Q0k0,
and substituting the value of A(ρ) into (17), the integral equation to be solved for unknown ψ(s) can be obtained as
(18)∫0a(∫0∞(e-2m2h-1)η(ρ)J0(rρ)J1(sρ)ρdρ)ψ(s)sds=Q0k0,
where
(19)η(ρ)=m2-m1F(ρ).
Defining new normalized variables and parameters such as
(20)s′=sa,r′=ra,ρ′=aρ,h′=ha,δ′=aδ,ψ(s)=ψ(as′)=ψ(s′),
the integral equation in (18) can be expressed in the form of
(21)∫01ψ(s)sds∫0∞J0(rρ)J1(sρ)ρdρ+∫01ψ(s)sds∫0∞2e-2m2hη(ρ)J0(rρ)J1(sρ)ρdρ-∫01ψ(s)sds∫0∞(2η(ρ)+1)J0(rρ)J1(sρ)ρdρ=q0,
where q0=2Q0/k0 and the prime sign is removed for the sake of simplicity. The first double integral in (21) can be expressed in terms of the first and second kinds of elliptic integrals, K and E, respectively. As k→1, the second kind of elliptic integral E(k) has a finite value while the first kind of elliptic integral K(k) has a logarithmic singularity such as
(22)K(k)=log(41-k2),E(k)=1.
Defining a new function M(s,r), given in the appendix, and using some algebraic manipulations, a Cauchy type singularity can obtained by using the properties of ψ(s) in (14). Hence, the first integral of (21) becomes
(23)∫01ψ(s)sds∫0∞J0(rρ)J1(sρ)ρdρ=1π∫-11ψ(s)s-rds+1π∫01(M(s,r)-1s-r+M(s,r)-1s+r)ψ(s)ds.
The second double integral in (21) has an exponential integrand with negative exponent so that as ρ→∞, the integrand asymptotically approaches to zero. The asymptotic expansion of the integrand can be expressed as
(24)e-2hρρ(-1+hδ2-δ22ρ-h2δ4-2hδ325ρ2+(1/3)h3δ6-h2δ5-2hδ4+2δ327ρ3-(1/3)h4δ8+(4/3)h3δ7-8h2δ6-40hδ5211ρ4+O(ρ-5)h2δ4-2hδ325ρ2).
Since there is no singularity and any discontinuity over the interval, and due to the asymptotic expansion as ρ→∞, the infinite integral can be approximated over the interval [0,B]. Depending on parameters h and δ, the value of B can be chosen large enough to obtain small enough value of the integral (less than 10-25) over [B,∞). Hence, the second double integral in (21) may be expressed as
(25)∫01ψ(s)sds∫0∞2e-2m2hη(ρ)J0(rρ)J1(sρ)ρdρ=∫01ψ(s)sds∫0B(h,δ)2e-2m2hη(ρ)J0(rρ)J1(sρ)ρdρ.
Finally, the last integral in (21) can be evaluated by using a series
(26)GA(ρ)=-δ22+δ3261ρ2-δ5291ρ4+5δ72141ρ6-⋯+429δ152301ρ14-O(ρ-16),
an asymptotic expansion of the integrand G(ρ)=ρ(2η(ρ)+1) as ρ→∞. Let us define a value C that depends on parameters h and δ and superposes the infinite integral as follows:
(27)∫0∞G(ρ)J0(rρ)J1(sρ)dρ=∫0C(h,δ)G(ρ)J0(rρ)J1(sρ)dρ×∫C(h,δ)∞GA(ρ)J0(rρ)J1(sρ)dρ+∫C(h,δ)∞[G(ρ)-GA(ρ)]J0(rρ)J1(sρ)dρ.
For a particular value of C, G(C)-GA(C)≈0 can be obtained; then the last infinite integral can be expressed as
(28)∫0∞G(ρ)J0(rρ)J1(sρ)dρ=∫0C(h,δ)G(ρ)J0(rρ)J1(sρ)dρ+∫C(h,δ)∞GA(ρ)J0(rρ)J1(sρ)dρ.
Since the function G(ρ) is a smooth function over the interval [0,C], it can be evaluated numerically using Gaussian quadrature along with Chebyshev orthogonal polynomials. On the other hand, the evaluation of the integral over the interval [C,∞) can be evaluated separately for each term in GA(ρ) like
(29)∫C(h,δ)∞GA(ρ)J0(rρ)J1(sρ)dρ=-δ22∫C(h,δ)∞J0(rρ)J1(sρ)dρ+∫C(h,δ)∞(δ3641ρ2-δ55121ρ4+5δ7163841ρ6-⋯)×J0(rρ)J1(sρ)dρ,
where
(30)H0(r,s)=∫C(h,δ)∞J0(rρ)J1(sρ)dρ={1s-∫0C(h,δ)J0(rρ)J1(sρ)dρ,s>r,12s-∫0C(h,δ)J0(rρ)J1(sρ)dρ,s=r,-∫0C(h,δ)J0(rρ)J1(sρ)dρ,s<r,
and the other integrals in (29) can be evaluated iteratively for k=1,2,3,… as follows:
(31)Hk(r,s)=∫C(h,δ)∞J0(rρ)J1(sρ)ρkdρ=J0(rC)J1(sC)kCk-1+skLk-1(r,s)+rkNk-1(r,s),
where
(32)Lk(r,s)=∫C(h,δ)∞J0(rρ)J0(sρ)ρkdρ=J0(rC)J0(sC)(k-1)Ck-1-rHk-1(s,r)k-1-sHk-1(r,s)1-k,(33)Nk(r,s)=∫C(h,δ)∞J1(rρ)J1(sρ)ρkdρ=J1(rC)J1(sC)(k+1)Ck-1+rHk-1(r,s)k+1+sHk-1(s,r)k+1.
The initial values of the integral H1(r,s), L2(r,s), and N3(r,s) are shown in the appendix.


## 4. Numerical Evaluation of Integral Equation

The integral equation in (21) can be solved using Gaussian quadrature method. Using the condition given in (14), the unknown function ψ(s) can be defined in terms of the truncated Nth term series expansion of Chebyshev orthogonal polynomials of the first kind, Tn(s), as
(34)ψ(s)=∑n=1NA2n-1T2n-1(s)1-s2,-1<s<1,
where A2n-1 are coefficients. Substituting the series expansion of the unknown function ψ(s) into the first integral in (23), it can be found that
(35)1π∫-11ψ(s)s-rds=∑n=1NA2n-11π∫-11T2n-1(s)(s-r)1-s2ds=∑n=1NA2n-1U2n-2(r),
where U2n-2(r) is the Chebyshev polynomials of the second kind, and eliminating the logarithmic singularities, shown in appendix, the second integral in (23) can be written as
(36)1π∫01(M(s,r)-1s-r+M(s,r)-1s+r)ψ(s)ds=∑n=1NA2n-11π×∫01(M(r,s)-1s-r+M(r,s)-1s+r+12rlog|s-rs+r|)×T2n-1(s)1-s2ds.
Finally, substituting the truncated series representation of ψ(s) into (25) and (29) system of algebraic equation can be obtained the to be solved for A2n-1,
(37)∑n=1NA2n-1(U2n-2(ri)+12riT2n-1(ri)2n-1+Z1(ri)+Z2(ri)-Z3(ri)T2n-1(ri)2n-1)=q0,
where ri, i=1,2,3,…, are collocation points and
(38)Z1(ri)=1π∫01(M(ri,s)-1s-ri+M(ri,s)-1s+ri+12rilog|s-ris+ri|)T2n-1(s)1-s2ds,Z2(ri)=∫01(∫0B(h,δ)2e-2m2hη(ρ)J0(riρ)J1(sρ)ρdρ)×sT2n-1(s)1-s2ds,Z3(ri)=∫01(∫0C(h,δ)G(ρ)J0(riρ)J1(sρ)dρ+∫C(h,δ)∞GA(ρ)J0(riρ)J1(sρ)dρ)×sT2n-1(s)1-s2ds.


## 5. Results

Because of the nature of the problem it is necessary to increase the density of the collocation points near the endsr=±1. Thus, these points may be selected as follows:
(39)Tn(ri)=0,ri=cos((2i-1)π2N),i=1,2,…,N.
Then, we get a (N×N) system of equations whose solution gives the coefficients A2n+1. With known coefficient values, the temperature distribution around the insulated barrier may be obtained by integrating (13) such as
(40)T2(r,0+)-T2(r,0-)=∑n=1NA2n-1∫-1r/aaT2n-1(s)1-s2ds.
Defining new variable like
(41)s=cosθ,π≤θ≤arccos(ra),
the integral in (40) can be evaluated using the relation
(42)Un(t)=sin{(n+1)arccost}sin(arccost),
and the difference in temperature distribution on the plane of the insulated barrier can be  obtained as
(43)T*(r)=T2(r,0+)-T2(r,0-)a=-1-(ra)2∑n=1NA2n-1U2n-2(r/a)2n-1.



---

*Source: 968464-2013-09-09.xml*