0 Introduction

Ceramic thermal barrier coatings (TBCs) have been applied to metal substrates to make them resistant to thermal shock in the aerospace and electronic devices, due to its many well-known advantageous properties, such as insulation and corrosion resistance. The accurate prediction of heat propagation in these devices may be important in ensuring their efficient work and design. Its heat transfer and thermal stress properties have been analyzed by the Fourier law, which is traditionally used to analyze the behavior of macroscopic heat transport. Recently, as nanotechnology and ultrafast processes continue to develop, TBCs with thicknesses between 10 nm and 100 nm are being used extensively in devices (Zeng and Chen^[1]). The average free path of the energy carrier has been reached (microscale effect in size), resulting in the discontinuity of heat conduction because of the insufficient number of collisions of the energy-required heat conduction. In addition to microscale effect in time, the high-speed thermal shock reaches femtosecond or picosecond order^[2], which is shortened into the range of the heat carrier mean free time, resulting in the premise of Fourier heat conduction theorem, in which the heat conduction velocity is infinite, thus not established. Therefore, there are obvious differences between the heat conduction theory at the microscopic scale and those at the macroscopic scale. New research indicates that the heat transfer mechanics at the micro-and nano-scales may deviate significantly from those at the macro-scale. For example, by the experiment of short-pulse laser heating on metals, Qiu and Tien^[3] have verified that the actual heat flux is not as fast as the Fourier law prediction.

Many scholars have studied and published non-Fourier heat transport models that are based on different physical theories. Maxwell^[4] first reported that the heat propagation speed is limited and modified the Fourier law by adding a relaxation time. Landau^[5] proved heat waves in liquid helium Ⅱ propagate with a finite speed. The CV model by modifying the Fourier law to interpret the thermal wave phenomenon was recently applied by some researchers (Cattaneo^[6] and Vernotte^[7]). Guyer and Krumhansl^[8] brought forward the phonon scattering model by directly solving the linearized Boltzmann equation for heat transport by phonon collision and scattering. This model was improved later by Joseph and Preziosi^[9], who emphasized the interrelation with a Jeffrey's type of heat flux equation. Tzou^[10] proposed the dual-phase-lag model with an emphasis on the lagging behavior at the micro-scale. Joshi and Majumdar^[11-12]solved the linearized Boltzmann transport equation to derive the EPRT and analyzed heat transfer in dielectric thin films. The EPRT is valid for microscale in both time and size. Qiu et al.^[13-14]proposed a two-step model with a quantum mechanical and statistical basis for a short-pulsed laser heating a metal film. Choudhuri^[15] established the three-phase delay hyperbolic heat conduction equation (TPL) by joining the phase delay of heat flow, temperature gradient and thermal displacement gradient.

At present, lots of studies researched the objects in the fields of biology and semiconductor. Li^[16] compared the analytical and numerical solutions of non-Fourier temperature field of laminated plate and laminated cylinder, and proved that the numerical solution obtained by backward difference separation is absolutely stable. Atfi and Talaee^[17-18]studied the non-Fourier thermal response of finite hollow cylinder subjected to periodic thermal shock based on C-V model. Cao^[19]researched the dynamic response of piezoelectric pole under thermal shock of heat source based on fractional order theory. Yu^[20]solved the temperature field and stress of the laminated material. On this basis, he studied the stress distribution at the crack tip under thermal shock. Li^[21] analyzed the fracture behavior of solid plates under thermal shock by the non-Fourier. Zhang and Li^[22-23] analyzed transient thermal stress and intensity factors for a circumferential crack in a hollow cylinder by generalized fractional heat conduction. Guo^[24]analyzed the fracture mechanics of plate, sandwich plate, coating, cylinder and other structures with the two-phase delayed model. Zhang and Chen^[25] studied the temperature field, displacement field and stress intensity factor of cracked hollow cylinder under thermal shock by using single-phase delayed non-Fourier heat conduction mode. Mao et al.^[26] studied silicon crystal with point defect scatter by molecular dynamics. Those researches are mainly aimed at only on the TBCs rather than the TBCs structure (thermal barrier coating and substrate), and kingdom is mainly in heat conduction.

In this paper, the equations combining the EPRT (for the TBCs) and the Fourier law (for the substrate) are numerically solved by the temperature field under thermal shock. The thermal stresses in the ceramic TBCs and J-integral of crack at the ceramic TBCs surface are studied by the finite element method, and the results are compared with those obtained by the Fourier law on both the TBCs and the substrate. The influence of the physical heat properties of the TBCs (such as the relaxation time and phonon speed) on the temperature field and the thermal stress as well as J-integral at an edge crack tip is analyzed.

1 Calculated Temperature Field and Stress and J-Integral 1.1 Heat Transportation Theoretical Model

Thermal shock in the direction of the x-axis has been studied (Fig. 1).

Because the thickness size of the ceramic TBCs is comparable with the mean free path(MFP) of typical phonons, which are considered to be the major heat carriers in the dielectric material, the phonon intensity for the ceramic TBCs ((0≤x≤b)) is numerically predicted by the EPRT^[11]:

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{v}\frac{{\partial I_\omega ^ + }}{{\partial t}} + \mu \frac{{\partial I_\omega ^ + }}{{\partial x}} = }&{}\\ {\frac{{\frac{1}{2}\int_{ - 1}^1 {{I_\omega }} {\rm{d}}\mu - I_\omega ^ + }}{{\tau v}}}&{, \;\;\;{\kern 1pt} \mu > 0}\\ {\frac{1}{v}\frac{{\partial I_\omega ^ - }}{{\partial t}} + \mu \frac{{\partial I_\omega ^ - }}{{\partial x}} = }&{}\\ {\frac{{\frac{1}{2}\int_{ - 1}^1 {{I_\omega }} {\rm{d}}\mu - I_\omega ^ - }}{{\tau v}}, }&{\mu < 0} \end{array}} \right. $

(1)

where I_ω is the phonon intensity in heat transport, v is the phonon speed, and τ is the phonon relaxation time. Furthermore, μ=cosθ, where θ is the angle between the phonon propagation and the x-axis, ω is the phonon frequency (ω_D=κ_BT_D/ ħ, where ω_D is the Debye cutoff phonon frequency, κ_B=1.381×10^-23 J/K is Boltzmann's constant, T_D= 647 K is the Debye temperature, and ħ =6.626×10^-34/2π J·s is Plank's constant). One could assume that the phonon intensity is in energy equilibrium under per frequency. The ''-'' and ''+'' superscripts represent the backward and forward transport phonons in the correspondence with the directional cosine.

Once, the phonon intensity is numerically solved by Eqs. (1), the temperature can be calculated by the Bose-Einstein distribution function at an equilibrium state^[11]:

$ \begin{array}{l} I_\omega ^0(T) = \frac{1}{2}\int_{ - 1}^1 {{I_\omega }} {\rm{d}}\mu = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_p {{v_p}} \frac{{\hbar \omega D(\omega )}}{{\exp \left[ {\frac{{\hbar \omega }}{{{\kappa _B}T(x)}}} \right] - 1}} \end{array} $

(2)

where I_ω⁰ is the equilibrium phonon intensity obtained by the equilibrium temperature, D(ω) is the phonon density of states per unit volume (D(ω)=ω²/2π²v³), and p is the index of the phonon polarization direction^[11].

Because the thickness of the substrate (b≤x≤b+ h) is much larger than the MFP, the temperature field of the substrate is predicted by the Fourier law^[25]:

$ \frac{{{\partial ^2}T}}{{\partial {x^2}}} = \alpha \frac{{\partial T}}{{\partial t}} $

(3)

where α is the thermal diffusivity.

1.2 Initial Conditions and Boundary Conditions

At t=0, the ceramic TBCs and the substrate are assumed to be at a uniform temperature T₀ and remain adiathermal.

At t>0, the temperature at x=0 is instantaneously raised to T₁. For EPRT, the intensity function for incoming phonon remains constant, and for Fourier's law, the boundary condition is constant temperature T₁. An idealized thermal shock on the surface of the TBCs is assumed, i.e., the heat transfer coefficient is infinitely large. Under this condition, the thermal stress is the greatest.

Because the substrate material is a metal or an alloy, the phonons are absolutely reflected at the bottom of the TBCs.

One assumes that there is no heat loss, and the two materials are in perfect thermal contact at the interface of the ceramic TBCs and substrate, and at the bottom of the substrate the temperature is a constant T₀.

The above conditions of the phonon and heat transport are expressed as follows:

$ \left\{ {\begin{array}{*{20}{l}} {t = 0:I_\omega ^ + (0, x, \mu ) = }&{}\\ {I_\omega ^ - (0, x, \mu ) = I_\omega ^0\left( {{T_0}} \right), }&{0 \le x \le b}\\ {t = 0:T = {T_0}, }&{b \le x \le b + h}\\ {t = 0:\frac{{\partial T}}{{\partial t}} = 0, }&{b \le x \le b + h}\\ {x = 0:I_\omega ^ + (t, 0, \mu ) = I_\omega ^0\left( {{T_1}} \right), }&{t > 0}\\ {x = b:I_\omega ^ - (t, b, \mu ) = I_\omega ^ + (t, b, \mu ), }&{t > 0}\\ {x = b + h:\frac{{\partial T}}{{\partial x}} = 0, }&{t \ge 0} \end{array}} \right. $

(4)

1.3 Temperature Field Numerical Calculation

Eq. (1) are two wave equations in the integrodifferential form. Eq. (1) are numerically solved by an explicit upstream differencing method, i.e., by backward differencing when the phonon intensity wave is moving in the positive x-direction (the first equation in Eq.(1)) and by forward differencing when it is moving in the negative x-direction (the second equation in Eq. (1)). To guarantee the stability of the numerical computation, the forward difference in time is applied while ensuring that Δt≤Δx/v μ in Eq. (1) and Eq. (3). The range of the cosine of the angle μ(-1≤μ≤1) is divided into K divisions. The integral in Eq. (1) can be approximated by a Gaussian quadrature^[12]as

$ \int_{ - 1}^1 I (\mu ){\rm{d}}\mu = \sum\limits_{i = 1}^k {{w_i}} I\left( {{\mu _i}} \right) $

(5)

where μ_i shows discrete directions and w_i is the weighting factor. The integrodifferential equations can be converted into the differential equation of the unknown I_i. In Eq. (3), a complete implicit difference method is used to calculate temperature field. Eq. (1) and Eq. (3) can be represented in the finite difference expressions as follows:

$ \left\{ {\begin{array}{*{20}{l}} {I_{j, \omega i}^{n + 1, + } = \left( {1 - {\mu _i}A - B} \right)I_{j, \omega i}^{n, + } + {\mu _i}AI_{j - 1, \omega i}^{n, + } + }&{}\\ {(B/2)\sum\limits_{i = 1}^k {{w_i}} I_{j, \omega i}^n, }&{\mu > 0, 0 \le x \le b}\\ {I_{j, \omega i}^{n + 1, - } = \left( {1 + {\mu _i}A - B} \right)I_{j, \omega i}^{n, - } - {\mu _i}AI_{j + 1, \omega i}^{n, - } + }&{}\\ {(B/2)\sum\limits_{i = 1}^k {{w_i}} I_{j, \omega i}^n, }&{\mu < 0, 0 \le x \le b}\\ { - \alpha \lambda T_{j - 1}^{n + 1} + (1 + 2\alpha \lambda )T_j^{n + 1} - }&{}\\ {\alpha \lambda T_{j + 1}^{n + 1} = T_j^n, }&{b \le x \le b + h}\\ {I_{j, \omega i}^{0, + } = I_{j, \omega i}^{0, - } = I_\omega ^0\left( {{T_0}} \right), }&{0 \le x \le b}\\ {T_j^0 = {T_0}, }&{b \le x \le b + h}\\ {T_j^{ - 1} = T_j^0, }&{b \le x \le b + h}\\ {I_{0, \omega i}^{n, + } = I_\omega ^0\left( {{T_1}} \right), }&{t > 0}\\ {I_{D, \omega i}^{n, - } = I_{D, \omega i}^{n, + }, }&{t > 0}\\ {T_{L + 1}^n = T_{L - 1}^n, }&{t \ge 0} \end{array}} \right. $

(6)

where the superscript n represents nth time step, subscript i is the index of the unit direction vector, subscript j represents the jth node, A=vΔt/Δx, B=Δt/τ, D is the last node of TBCs, and L is the last node of the substrate. Material parameters^[27] in the TBCs and substrate are shown in Table 1.

1.4 Thermal Stress and Fracture Calculation by Finite Element

In this work, the thermal stress and J-integral at the crack tip in the TBCs are analyzed using the finite element method after numerical calculation of the temperature field. The numerical calculations make use of the element software MSC.MARC, which can be applied in the thermal stress and J-integral^[28]. The stress and temperature analyses are uncoupled. The TBCs and the substrate are modeled by a four-node plane strain and quadratic elements in the stress analysis, because the model is a good elastic body without plastic deformation. In fracture mechanics analysis, it is assumed that the fracture is an edge crack that opens to the surface and its length is 2×10^-8 m. Quarter point 8-node isoparametric degenerate singular elements are employed in modeling the region near the crack tip, as shown in Fig. 2.

The assumed boundary conditions put constraints on the displacement field in the y-axis. The displacements only take place along the x-axis, but at the bottom of substrate, the x-displacement is assumed to be zero. The initial temperature should be T₀ at all nodes. The thermal loads are derived from temperature fields. The parameters and the material properties used in the heat transport, thermal stress, and J-integral numerical calculations are listed in Tables 1 and 2^[29-32]. The parameters and properties remain constant unless otherwise specified in the numerical calculations and subsequent analysis.

2 Results and Discussion 2.1 Temperature Distribution and Comparison

Figs. 3 and 4 show temperature field from EPRT and the Fourier law, respectively. It is clear that the temperature jumps on the TBCs surface are calculated by the EPRT, which are not shown by the Fourier law. In the TBCs and substrate, the amplitude of temperature attenuation obtained from the EPRT is lower than the attenuation from the Fourier law. Because the heat transfer velocity is finite according to the EPRT, where the heat transport velocity is represented by income and out-coming of phonon, it is infinite for the Fourier law.

2.2 Results of Thermal Stress and J-integral

The comparison between the thermal stress in the y-direction obtained from the EPRT and that obtained from the Fourier law is presented in Fig. 5. Both of the thermal stress profiles along the x-direction are similar, i.e., the thermal stresses gradually change from tensile stress to compressive stress along the x-direction in the TBCs and are only compressive stresses along the x-direction within the substrate. But the thermal stress, special tensile stress in the TBCs calculated with the EPRT, is lower than that calculated from the Fourier law, because temperature contrast from the surface of the TBCs to the bottom of the substrate, which is obtained from the EPRT, is lower than that obtained from the Fourier law. At the interface, a thermal stress discontinuity occurs for the difference in the characters between the two materials.

For reasons given above, the amplitude of temperature attenuation from the EPRT is lower than the attenuation from the Fourier law, and the J-integral from the EPRT is lower than the results from the Fourier law in the TBCs and changes slowlier with time, and J-integral from the Fourier law decreases with time in Fig. 6, because temperature difference between the top of TCBs and bottom of substrate from EPRT changes slowlier with time than that from the Fourier law.

3 Influence of the TBCs Heat Physical Properties on Temperature, Thermal Stress and J-Integral 3.1 Influence of Phonon Speed

The temperature field, thermal stress profiles and J-Integral at the crack tip predicted by the EPRT with different phonon speeds (v₁=8×10² m/s, v₂=1.3804×10³ m/s, v₃=5×10³ m/s) are shown in Figs. 3, 7, 8, 9 and 10.

At the surface of the TBCs, the amplitude of temperature attenuation decreases with increasing phonon speed. However, in the TBCs, the amplitude of temperature attenuation increases with increasing phonon speed, because the intensity function of incoming phonon remains constant at the surface of the TBCs for the EPRT, and the intensity of out-coming phonon increases with increasing phonon speed, so the temperature jump increases, whereas in the TBCs, the phonon transport increases with increasing phonon speed, so the heat transport increases.

At the surface of the TBCs, the magnitude of the thermal stress decreases as the phonon speed increases. But in the TBCs, the thermal stress rises with the rise of phonon speed. The J-integral at the crack tip decreases with increasing phonon speed.

3.2 Influence of Relaxation Time

The temperature field, thermal stress profile and J-Integral at the crack tip predicted by the EPRT with different relaxation times (τ₁=3×10^-12 s, τ₂=6×10^-12 s, τ₃=1.2×10^-11 s) are shown in Figs. 3, 11, 12, 13 and 14.

At the surface of the TBCs, the amplitude of temperature attenuation descends with the increase of relaxation time. But the amplitude of temperature attenuation increases with increasing relaxation time in the TBCs.

On the surface of the TBCs, the thermal stress decreases with the increase of relaxation time. Nevertheless, in the TBCs, the thermal stress and compressive stress increase for the same reason as stated above. The J-integral at the crack tip decreases with increasing relaxation time.

4 Conclusions

This study numerically analyzed temperature field, thermal stress and J-integral obtained from the EPRT and from the Fourier law. The following conclusions can be drawn:

1) For nano-ceramic TBCs under thermal shock, the temperature field, thermal stress and J-integral derived from the EPRT are lower than those derived from the Fourier law.

2) The temperature attenuation and thermal stress decrease on the surface of the TBCs with increasing phonon speed and relaxation time, but the amplitude of temperature attenuation and thermal stress in the TBCs increase with increasing phonon speed and relaxation time.

3) J-integral at the crack tip decreases with increasing phonon speed and relaxation time.