Thermoelastic Analysis of a Cracked Substrate Bonded to a Coating Using the Hyperbolic Heat Conduction Theory

In this article, the dynamic temperature and thermal stresses around a crack in a substrate bonded to a coating are obtained using the hyperbolic heat conduction theory. Fourier and Laplace transforms are applied and the hyperbolic heat conduction and thermoelastic crack problems are reduced to solving singular integral equations. The overshooting phenomenon is observed and the crack kinking phenomenon under thermal loading is investigated by applying the criterion of maximum hoop stress. Numerical results show that the hyperbolic heat conduction parameters and the geometric size of the composite have significant influence on the dynamic stress field. It seems that high temperature loading on the surface may lead to crack kinking away from the surface and low temperature loading may cause crack kinking toward the coating. Moreover, the hyperbolic heat conduction theory may give more conservative results than that the Fourier heat conduction theory.


INTRODUCTION
High-rate heat transfer has become a major concern in modern industries especially in material processing, such as the pulsed laser heat and ultrasonic waves, and accurate heat conduction analysis is of great importance for the material and structural integrity. Investigation of the temperature and stress fields is essential to the safety design of the composite structures under severe temperature loading. The Fourier heat conduction model, although given sufficient accuracy for many engineering applications, implies infinite thermal wave propagation speed, and renders ineffective at the very small length and time scales associated with smallscale systems [1]. Examples include the transient temperature field caused by pulsed laser heating of thin structures, and the measured surface temperature of a slab immediately after an intense thermal shock is 300 C higher than that obtained from the parabolic heat conduction model [2]. Consideration of the hyperbolic heat conduction model becomes important if irreversible physical processes, such as THERMOELASTIC ANALYSIS OF BONDED CRACKED SUBSTRATE 271 crack or void initiation in a solid, are involved in the process of heat transport. In these cases, the hyperbolic heat conduction model should be used [3].
Inherent defects in materials such as dislocations and cracks may disturb the temperature distribution when thermal loading is applied, and singularities may be developed in the neighborhood of discontinuities. Some researchers have studied the heat conduction problems of cracked materials using the classical Fourier heat conduction model [4][5][6][7]. Extensive study of the distribution of thermal stress in the vicinity of a crack in an elastic body has been performed since 1950s [8,9]. General coupled thermoelastic theories have been adopted to study the fracture problems in thermoelasticity [10,11]. If the effects of both the inertial term and the thermoelastic coupling term are neglected, uncoupled thermoelastic theory can be applied to solve the thermal-elastic fracture problem [12,13]. In most of the usual engineering applications it is proper to use the uncoupled thermoelastic theory without significant error [14][15][16]. Numerous studies have been devoted to solving thermoelastic crack problems in advanced materials under thermal loading based on the classical Fourier heat conduction theory [17][18][19][20][21][22][23]. The influence of transient thermal gradients and substrate constraint on delamination of thermal barrier coating has been recently studied [24].
Some investigations on crack problems in thermoelastic materials have been made using the hyperbolic heat conduction model. The solution based on the hyperbolic heat conduction equation for a travelling point heat source around a propagating crack tip has been derived and the temperature distribution at the tip of a propagating crack was measured experimentally [25,26]. A cracked layer under thermal shock has been analyzed and the effect of second sound of Lord-Shulman theory has been studied [27]. The problem of a finite crack in a material layer under transient non-Fourier heat conduction was investigated in [28] and the problem of an interface crack in layered composite media under applied thermal flux was studied using the hyperbolic heat conduction theory [29]. A thermoelastic analysis of a cracked substrate under a thermal shock was given in [30] based on the hyperbolic heat conduction theory; and based on the same theory the transient temperature and thermal stresses around a partially insulated crack in a thermoelastic strip under a temperature impact were obtained [31]. The transient temperature field around a thermally insulated crack in a substrate bonded to a coating has been obtained by using the hyperbolic heat conduction model [32].
In this article, the thermoelastic problem of a cracked substrate bonded by a coating under transient thermal loading is studied using the hyperbolic heat conduction model. Fourier and Laplace transforms are employed to reduce the problem to solving singular integral equations. The asymptotic temperature and stress fields around the crack tip are obtained in an explicit form and Laplace inversion is applied to get the dynamic temperature and stress fields. The overshooting phenomenon is observed when applying the hyperbolic heat conduction model and the crack kinking phenomenon is investigated by applying the criterion of maximum hoop stress.

STATEMENT OF PROBLEM AND BASIC EQUATIONS
Consider a thermoelastic substrate containing a crack of length 2c parallel to the interface between the substrate and the coating, as shown in Figure 1 composite is initially at certain temperature and the free surface of the coating y = −h = − a + b is suddenly heated to by a temperature change T 0 . The crack surfaces are assumed to be thermally insulated.

Heat Conduction Equations
In the heat-transfer situations include extremely high temperature gradients, extremely large heat fluxes or extremely short transient durations, Fourier's law may be modified to a relation of the type [33], where q is the heat flux, T is the temperature, k is the thermal conductivity of the material, is the spatial gradient operator, and is the so-called relaxation time. The local energy balance equation with vanishing heat source follows: where and C are the mass density and the specific heat capacity, respectively. Incorporating Eq. (1) with Eq. (2) leads to the hyperbolic heat conduction equation for the substrate and the coating, where 2 is Laplace's differential operator, i , d i = k i / i C i and k i (i = 1 2) are the relaxation times, the thermal diffusivities and the thermal conductivities of the substrate and the coating, respectively. It is noted that here and afterwards the subscript "i = 1 2" denote the quantities of the substrate and the coating, respectively. By introducing the dimensionless variables where T 0 is the reference temperature and d 0 is the reference thermal diffusivity, the governing equations (3) have the following dimensionless forms: It is noted that here and afterwards, the hat "¯" of the dimensionless variables is omitted for simplicity. The hyperbolic heat equation (5) is subjected to the following boundary and initial conditions in dimensionless forms It should be noted that the relaxation time ( ) for most engineering materials is of the order of 10 −14 − 10 −6 s and the thermal diffusivity 10 −8 − 10 −3 m 2 /m 2 s, and therefore the Fourier parabolic heat conduction model can give good results. Experiments have shown that some special materials may have thermal relaxation time up to the order of 10s, which are important materials to work as thermal insulators [34]; in this case it is necessary to take into account the effect of the relaxation time and use the hyperbolic heat conduction model. Other applications of hyperbolic heat conduction model are in the area of transient thermal disturbance with small length and time scales.

Thermal-Elastic Field Equations
The equilibrium equations, the strain-displacement relations, the compatibility equation and the constitutive law of plane stress thermoelasticity can be expressed as follows: where u x and u y are components of displacement, E, v and are Young's modulus, Poisson's ratio and the coefficient of linear thermal expansion, respectively.
Let U x y be the Airy stress function, then the stresses can be expressed in terms of U as Substitution of Eq. (12) into (10) and (11) leads to By introducing the following dimensionless quantities the governing equations (13) can be rewritten in the following dimensionless forms: The mechanical boundary conditions can be expressed as

Temperature Field
Applying the Laplace transform to Eqs. (5): where Br stands for the Bromwich path of integration. By assuming that the composite is at rest in the beginning, we can have . Following the procedure in [30], the temperature field in the Laplace domain satisfying the boundary condition and regularity condition can be expressed as and D j = D j p (j = 1 2) are unknowns to be determined, and where constants C 0 1 , C 0 2 and f can be obtained from the boundary conditions as detailed in Appendix A. As the unknown functions are dependent, all other unknown functions can be expressed by only one independent unknown function, for example D 1 , as follows Introduce the temperature density function as It is clear from the boundary conditions (27) Substituting (21) and (22) into (27) and using Fourier inverse transform, we can get Substituting (21) and (22) into (7) and applying the relation (29), we get the singular integral equation for t = ct as follows where the kernel function K x t is given as and R = 1 − 2 / 1 . The singular integral equation (30) under the singlevaluedness condition in (28) has the following form of solution [35]: where x is a bounded and continuous function. From the condition (28), it is clear that x is an odd function of x, i.e., −x = − x . Function x can be solved numerically, as detailed in [32], and function D 1 can be calculated by using the Chebyshev quadrature for integration as By substituting Eq. (33) into (21)-(23) the temperature in the Laplace domain can be obtained, and the temperature in the time domain can be determined by applying the Laplace inverse transform.

Thermal Stresses
Considering the temperature expressions (23), the general solution of the Airy function in the Laplace domain satisfying the regular condition at infinity can be expressed as where B i , C i (i = 1 − 4), A 1 and A 2 are unknowns to be determined, and ij (i = 1 − 3 j = 1 2) are known functions determined from the boundary conditions as shown in Appendix A.
The stress components can be obtained by substituting Eqs. (34)(35)(36) into Eq. (12), the components of displacement can be determined by integration of strain, and the satisfaction of the boundary conditions (16)(17)(18)(19) leads to the following relations: i.e., the unknowns B i , C i (i = 1 − 4), can be expressed as functions of the two independent unknown functions A 1 and A 2 , and the coefficients X mi and Z mi (i = 0 − 4) are given in Appendix A.
Denote the jumps of displacement across the plane y = 0 in the Laplace domain by u x and u y , Following the procedure in [36] and introducing two dislocation density functions f j x (j = 1 2) as By applying the boundary conditions (20), the singular integral equations for f i x (i = 1 2) are obtained as: where functions U i x , M ij x t (i j = 1 2) are known functions defined in Appendix A. Functions f j x (j = 1 2) also need to satisfy the single-valuedness condition The solution of the integral Eqs. (40), f j x (j = 1 2 , can be expressed as and the singular integral equations can be reduced to solving a system of algebraic equations in terms of F 1 t p and F 2 t p as [37,38]: where,

Asymptotic Stress Field Near Crack Tip
By considering the asymptotic nature of the integrands for large values of the integration variable and following the procedure in [36,37] and using the asymptotic formula [39]: The singular stresses near the crack tip can be obtained as where the dynamic SIFs K 1 t and K 2 t are given by the stress intensity factors (SIFs) in the Laplace domain K * 1 p and K * 2 p are defined as and r are the polar coordinates measured from the crack tip defined by The hoop and shear stresses at an angle near the crack tip can be obtained in terms of the polar coordinates r as [40] r t = y r t cos 2 + x r t sin 2 − xy r t sin 2 r r t = sin 2 y r t − x r t /2 + xy r t cos 2 (51) Upon substituting Eqs. (45)(46)(47) in to Eqs. (51), the hoop stresses and the shear stresses can be expressed as These expressions are in the same form for the stress state in the neighborhood of the crack tip in the case of plane strain or generalized plane stress [41]. The hoop stress intensity factor and shear stress intensity factor associated with the hoop and shear stresses at an arbitrary angle can be defined, respectively, as If we are interested in the initiation of the slow crack extension in a given material, the criterion of maximum hoop stress can be applied to predict the crack extension direction. This criterion imply that the crack will start to grow from the tip in the direction along which the hoop stress (intensity factor) is maximum and the shear stress (intensity factor) is zero, providing that the fracture toughness of the material is homogeneous in all directions [41]. It should be noted that the condition K r = 0 is a generalization of the crack propagation condition K 2 = 0 for mixed mode crack extension [42].
To find the angle of crack kinking, from Eq. (52), we can write the derivative of with respect to equal to zero and let the second derivative of be negative. It is evident that the crack kinking angle depends on the values of K 1 and K 2 , and the analytical solution would be very complex if not impossible. In this article, we show the solution of crack kinking angle in graphical form.

NUMERICAL RESULTS AND DISCUSSIONS
The temperature field in the time domain can be obtained from Eqs. (23) and (33) by applying numerical inversion of Laplace transform, as detailed in [43] and Appendix B. In the computation of the temperature field the values of the parameters used in the algorithm for the numerical inversion of Laplace Transform are chosen as N = 10, = 0, 0 2 ≤ ≤ 0 3, such that the function f t can be best described within a particular period of time.
Assuming the temperature boundary conditions on the coating to be and the geometric size of the substrate and coating as a/c = 1, b/c = 0 3, the material properties k 1 = 1, k 2 = 0 5, = 0 4, 1 = 2, 2 = 1 and the boundary temperature condition T 0 = +1, the temperature field and the stress perturbation around the crack can be obtained consecutively. The dynamic temperature at the mid-points of the crack faces for hyperbolic heat conduction model and Fourier model are displayed in Figure 2, in which "UF" and "LF" denote the temperatures of the midpoints on the upper and lower crack faces, respectively. The results based on the hyperbolic heat conduction model show that the temperature oscillates at the very beginning of temperature impact loading, and then increases to reach a peak point and fluctuate and converge to a steady-state value as time becomes long enough. It is shown that the temperature on the mid-point of the upper face may exceed Figure 2 The variation of temperature with time at the crack face mid-points.

281
the temperature on the boundary. The results based on the Fourier heat conduction is the special case when 1 = 2 = 0, and the comparison to the hyperbolic heat conduction model shows the effect of the relaxation time on the thermal response.
The temperature distribution near the crack at the time instant (t = 3 5 is shown in Figures 3a and 3b for the hyperbolic heat conduction model and Fourier model, respectively. The disturbance of the thermal-insulate crack on the temperature field can be observed from the iso-temperature lines, and there is a temperature jump across the crack faces. For hyperbolic heat conduction model, due to the interference of the insulated crack, temperature waves may interact and reinforce each other and result in the higher temperature in the inner region of the heat conduction medium than that on the boundary, leading to the "overshooting phenomenon", which is of great importance in thermal engineering applications such as safety design of the electronic or mechanical devices under severe thermal loadings [44]. There is no overshooting phenomenon found based on the Fourier heat conduction model.
The effect of the thickness of coating b/c on the dynamic SIFs is shown in Figure 4 when the geometry of the cracked substrate is a/c = 1, the material properties are k 1 = k 2 = 1, = 0 4 and 1 = 1, 2 = 0 5 and the boundary temperature condition T 0 = +1. The SIFs oscillate and increase in magnitude till they reach the peak values and oscillate for some time before approaching their static values. The magnitude of the peak value decreases as the size ratio b/c increases, which indicates that a thicker coating may decrease the magnitude of SIF. Figure 5 exhibits the variation of the dynamic SIFs with different thermal relaxation times 1 and 2 when k 1 = k 2 = 1, = 0 4, a/c = 1, b/c = 0 5 and T 0 = +1. The results based on the Fourier heat conduction model shows that the dynamic SIF increases smoothly to reach the peak values and then decrease to the static values. The oscillating feature of the dynamic SIFs is observed when the hyperbolic heat conduction theory is applied, and as the relaxation time of the coating decreases the magnitude of the dynamic SIF decreases. This conclusion indicates that the structure safety design based on the hyperbolic heat conduction theory may lead to more conservative design than Fourier's heat conduction theory.
The angular functions of the normalized dynamic hoop stresses 2 r/c and shear stresses r 2 r/c at a specific time instant t = 3 0 are displayed in Figure 6 when a/c = 1, b/c = 0 1, k 1 = k 2 = 1 and = 0 4. The maximum hoop stress criterion is used to study crack kinking in this article. It is shown that when the hoop stress reaches the maximum value, the shear stress is zero. When the temperature loading on the surface of the coating is positive, the maximum hoop  stress appears at the angle = +60 degrees, which indicates that the crack may kink in this direction; when a negative temperature is applied on the coating surface, the hoop stress reaches the maximum value at the angle = −70 degrees, which means that the crack may kink toward the coating in this direction. It seems that high temperature loading on the coating surface may lead to crack kinking away from the surface and low temperature loading may cause the crack to kink toward the coating. This is in agreement with the conclusion of [45] that a tensile residual stress in the coating increases the likelihood of kinking toward the coating while a compressive residual stress discourages kinking toward the coating.

CONCLUSIONS
The transient thermal stresses around a crack in a substrate bonded to a coating have been obtained using the hyperbolic heat conduction theory. Integral transform method is applied and the hyperbolic heat conduction and thermoelastic crack problem are reduced to solving singular integral equations. The overshooting phenomenon may be observed when the hyperbolic heat conduction is applied. The crack kinking phenomenon is investigated by applying the criterion of maximum hoop stress. Numerical results show that the hyperbolic heat conduction parameters, the material properties and the geometric size of the composite have significant influence on the dynamic stress field. Hot or cold temperature loadings may lead to different crack kinking directions. The results predicted by the hyperbolic heat conduction model are more conservative than that by the Fourier's heat conduction model.
The functions U i x , M ij x t , (i j = 1 2) are functions defined as: where, M 60 = L 61 X 10 + L 63 X 30 (A.13-2) X m0 m = 1 − 4 is defined as: (A.14-1) (A.14-2) (A.14-3)   where, Then X jk (j k = 1 − 4) are defined as  and it can be evaluated at N discrete points along the real positive p-axis in the Laplace domain: The function f t can be approximated by where P 0 n is a Jacobi polynomial and the coefficients C n can be determined by solving the following equations k n=0 k k − 1 · · · k − n − 1 + k + 1 + k + 2 · · · + k + n + 1 C n = · F * p k The parameters N , and are selected such that the function f t can be best described within a particular period of time.