Thermal shock fracture of a cylinder with a penny-shaped crack based on hyperbolic heat conduction

International Journal of Heat and Mass Transfer 91 (2015) 235–245

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Thermal shock fracture of a cylinder with a penny-shaped crack based on hyperbolic heat conduction S.L. Guo a,⇑, B.L. Wang a,b a b

Graduate School at Shenzhen, Harbin Institute of Technology, Harbin 150001, PR China Institute for Infrastructure Engineering, University of Western Sydney, Penrith, NSW 2751, Australia

a r t i c l e

i n f o

Article history: Received 17 April 2015 Received in revised form 22 July 2015 Accepted 22 July 2015

Keywords: Fracture mechanics Cylinder specimen Thermal shock Hyperbolic heat conduction

a b s t r a c t This paper studies the thermal shock fracture of a cracked cylinder based on the hyperbolic heat conduction. The crack faces are subjected to a sudden anti-symmetric thermal flux and a sudden symmetric thermal flux, respectively. By Laplace transform and dual integral equation technique, the mode II stress intensity factor and the mode I stress intensity factor are developed at the crack front for the two cases, respectively. Numerical results of stress intensity factor for selected thermal relaxation time and crack size are shown graphically. It is found that the stress intensity factor is considerably enhanced for large thermal relaxation time (which is a material constant) or small crack radius. In addition, the stress intensity factor at the crack front increases with the thermal relaxation time. For the case of anti-symmetric thermal flux, the mode II stress intensity factor increases rapidly with crack size. Whereas for the case of symmetric thermal flux, with increasing crack size, the mode I stress intensity factor increases slowly. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction

the non-Fourier heat conduction law and balance equation, the temperature governing equation is obtained as:

The classical Fourier heat conduction law assumes that a body will be affected at the instant of heating:

qðtÞ ¼ krTðtÞ;

ð1Þ

where t is time, q is the heat flux vector, k is the thermal conductivity and T is the temperature. This assumption means that the speed of heat propagation in the body is infinite. In practical situations, however, the speed of heat propagation in a body is always finite [1]. It has been proved that Fourier’s law is not accurate in some extreme cases like under highly-varying thermal loading condition, at ultra-low temperature and in nano materials. Cattaneo [2] and Vernottee [3] first independently formulated the non-Fourier heat conduction problem based on the local energy balance. The relaxation behavior is introduced to approach the wave nature of heat propagation and gives a non-Fourier heat conduction law:

qðt þ sq Þ ¼ krTðtÞ;

ð2Þ

where sq is the thermal relaxation time, which is related to the collision frequency of the molecules within the energy carrier. By ⇑ Corresponding author. Tel.: +86 075526032119. E-mail address: [email protected] (S.L. Guo). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.07.081 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.

qcsq

@2T @T @Q þ qc ¼ k r 2 T þ sq þ Q; @t @t @t 2

ð3Þ

where q is the mass density, c is the specific heat, and Q is the internal heat generation rate per unit volume. This model is named as hyperbolic heat conduction and its accuracy have been proved by some micro/macro scope experiments [4–6]. Basically, many studies have been carried out for the classic Fourier heat conduction model (e.g., Refs. [7–10] use temperaturedependent material properties and Ref. [11] studies the heat transfer in skin tissues). Based on the work of Cattaneo and Vernotte, researchers gave several significant solutions of thermal problems. The one dimensional heat conduction in semi-infinite body [12] and wall [13] were studied early. Then the two dimensional problem in cylinder was solved [14–16]. Later, Moosaie [17–19] solved one dimensional hyperbolic heat conduction in a finite medium subjected to arbitrary periodic/non-periodic surface disturbance, and a finite medium with insulated boundaries and arbitrary initial conditions. Recently, Zhao and Wu [20] analyzed a solid sphere under sudden surface temperature change, simple harmonic periodic surface temperature change, triangular surface temperature change and pulse surface temperature changes. In addition, Babaei and Chen [21] studied the hyperbolic temperature fields of a functionally graded hollow sphere and cylinder.

236

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

It is well known that material manufacture processing usually products defects or cracks in the materials, which may disturb the local temperature distribution and intensify the temperature gradient, which introduces high thermal stress intensity that may cause rapid crack growth. Hence, the research on the mechanical behaviors of materials under thermal environments is essential [22–24]. Early, Hasselman [25] established a unified theory of thermal shock fracture initiation and crack propagation in brittle ceramics. Wilson and Yu [26] present a method of using the J-line integral to extract the magnitude of crack tip stress intensity factors from displacement solutions for thermal stress crack problems. Then Sumi and Katayama [27], and Tsai [28] investigated the thermal stress in a finite rectangular plate with a Griffith crack and thermal stress in an infinite body containing a penny-shaped crack, respectively. Lee and Shul [29] calculated the thermal stress intensity factors for an insulated interface crack in an infinite two-dimensional elastic biomaterial an under uniform heat flow. Later, Lee and Park [30] evaluated the thermal stress intensity factors for a partially insulated interface crack. The problems discussed in literatures [27–30] are under steady state thermal load. However, in practice, the transient thermal problems are significant and interesting as well. Nied and Erdogan [31] studied the thermal shock problems for a circumferentially cracked hollow cylinder. Yu and Qin [32] completed a two-dimensional analysis of thermal and electric fields of a thermopiezoelectric solid damaged by cracks based on Fourier transformations. Later, they [33] developed a generalized self-consistent approximate method for determining the thermoelectroelastic properties of piezoelectric materials weakened by microcracks. Lu and Fleck [34] solved the transient thermal problem for an orthotropic plate with an edge crack. Ma et al. [35] considered A magnetoelectrically permeable interface crack between two semi-infinite magnetoelectroelastic planes under the action of a heat flow. Guo and co-workers investigated an analytical method [36] for thermal shock crack problems of a functionally graded cylindrical shell with a circumferential crack and a combined analytical–numerical method [37] for a functionally graded plate with a surface crack, respectively. The most existing researches on the transient thermal problems are based on classical Fourier heat conduction theory. Recently, Wang and co-workers [38–42], Chen and Hu [43] and Hu and Chen [44] analyzed some thermal shock problems of infinite media, semi-infinite media, strip and plate, under the framework of hyperbolic non-Fourier heat conduction. That inspires the current consideration of the thermal shock fracture of a cylinder with an embedded penny-shaped crack. This paper establishes a solution technique for the thermal shock fracture of cracked cylinder based on hyperbolic heat conduction. Laplace transform is used to deal with the time-varying behavior of the temperature and thermoelasticity fields. Solutions for the problem of thermally insulated crack (Section 3) and heated crack (Section 4) are obtained separately. Numerical results of stress intensity factor of the crack are evaluated to show the effects of thermal relaxation time and geometry size of the crack. Some significant differences between the classical heat conduction and the hyperbolic, non-classical heat conduction were observed.

2. Basic governing equations Consider the axisymmetric problem of a cracked cylinder given in Fig. 1, which will be loaded with a uniform and axisymmetric axial thermal flux in the follow sections. Thus all the field variables are functions of coordinates R (radical direction) and Z (axial direction) only. The radius of the cylinder is denoted as Rb and the crack radius is denoted as Ra. The constitutive equations for the thermal

Fig. 1. A cylinder with a penny-shaped crack.

flux and the governing equation for temperature under hyperbolic heat conduction in cylindrical coordinate system are:

qZ ðR; Z; tÞ ¼ k

@TðR; Z; tÞ @q ðR; Z; tÞ  sq Z ; @Z @t

ð4Þ

and

qcsq

@ 2 TðR; Z; tÞ @TðR; Z; tÞ þ qc @t @t 2

¼k

@ 2 TðR; Z; tÞ @R2

! 1 @TðR; Z; tÞ @ 2 TðR; Z; tÞ : þ þ R @R @Z 2

ð5Þ

It has been assumed that the thermal properties are not affected by the mechanical behaviors of the material and the heat source is neglected. In order to make the study more convenient, we introduce the following dimensionless parameters according to: t ¼ t=sq , r ¼ R=Ra , z ¼ Z=Ra , r b ¼ Rb =Ra , a ¼ Ra =l0 , where pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l0 ¼ ksq =ðqcÞ is a characteristic length parameter of the material. Accordingly, Eqs. (4) and (5) can be re-written as:

qz ðr; z; tÞ ¼ 

k @Tðr; z; tÞ @qz ðr; z; tÞ  ; @z Ra @t

ð6Þ

and 2

@ Tðr; z; tÞ 1 @Tðr; z; tÞ @ 2 Tðr; z; tÞ þ þ r @r2 @r @z2 ! 2   @ Tðr; z; tÞ @Tðr; z; tÞ ¼ a2 þ : @t 2 @t

ð7Þ

Both boundary and initial conditions are necessary for the solution of Eq. (7). Two cases of boundary conditions will be considered in Sections 3 and 4 separately. The initial conditions Tðr; z; 0Þ ¼ T 0 and qr ðr; z; 0Þ ¼ qz ðr; z; 0Þ ¼ 0 can simulate most practical cases. However, the thermal stress of interest is associated with temperature change DT ¼ Tðr; z; tÞ  T 0 . To simplify the deducing of solutions for mechanical fields, T 0 ¼ 0 is assumed, which means that DT ¼ Tðr; z; tÞ. This paper does not consider the temperature dependency of material properties. Thus, assumption of T 0 ¼ 0 does not affect the thermal stress level. Now, the initial conditions of Eq. (7) are:

Tðr; z; 0Þ ¼ 0;

@Tðr; z; 0Þ=@r ¼ @Tðr; z; 0Þ=@z ¼ 0

ð8Þ

In the following analysis, Laplace transform will be applied to Eqs. (6) and (7). Then the temperature field in Laplace transform will be obtained, and be used for solving the thermal stresses. Base on the linear thermoelastic theory, the mechanical constitutive equations are:

8 9 rrr > > > > > > < =

2

c11 c12 c13

rhh 16 6 c12 c11 c13 ¼ 6 > > Ra 4 c13 c13 c33 r zz > > > > : ; 0 0 0 rrz

9 8 9 38 @u=@r v11 > > > > > > > > > > > < = = u=r 0 7 7 11  T; 7 > > @w=@z 0 5> v33 > > > > > > > > > : ; : ; @u=@z þ @w=@r c44 0 0

ð9Þ

237

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

where rrr ; rhh and rzz are normal stresses, rrz is shear stress, cij (i; j ¼ 1; 2; 3) are elastic contents, u and w are displacements along X and Z directions respectively, and vii are temperature-stress coefficients. In the absence of body forces, the displacements governing equations deduced from the equilibrium equations are given as: ! ! @ 2 u 1 @u u @2w @2u @2 w @T c11 þ  þ c þ þ c ¼ Ra v11 13 44 @r 2 r @r r2 @r@z @z2 @r@z @r ! ! @ 2 u 1 @u @ 2 w 1 @w @ 2 u 1 @u @2w @T c44 þ þ þ þ þ c13 þ c33 2 ¼ Ra v33 @r@z r @z @r 2 r @r @r@z r @z @z @z ð10Þ

where Aa ðsÞ and Ba ðsÞ are unknown functions to be determined, and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k0 ðsÞ ¼ 1 þ a2 pðp þ 1Þ=s2 . Here and in the following a variable with a superscript * represents its Laplace transform. The thermal flux qz ðr; z; pÞ associated with Eq. (13) can be calculated from the Laplace transform of Eq. (6). The unknown functions Aa ðsÞ and Ba ðsÞ will be determined by the thermal boundary conditions of the cylinder. It follows from the thermal boundary conditions on crack surface (Eq. (11)) and on the flank of cylinder (Tðr ¼ rb ; z; tÞ ¼ 0) that

Z

0 6 r < 1;

ð11aÞ

1 < r 6 rb ;

ð11bÞ

Z

1

0

Z

Z

1

sBa ðsÞI0 ðrsk0 ðsÞÞds

0

¼

Usually, the propagation of heat is hindered by the cross crack. Then the approach models the thermal insulated crack in which the normal component of the thermal flux vector on the crack faces is zero. Assume the cylinder in Fig. 1 is subject to a couple of anti-symmetric uniform thermal load on its upper and lower surfaces. Naturally, the temperature field is anti-symmetric and the axial thermal flux is symmetric with respect to z ¼ 0 plane. By a proper superposition, the solution for the un-cracked cylinder under given thermal load has been separated and the following perturbed problem of the cracked cylinder only loaded with self-equilibrating crack surfaces thermal flow is considered (Fig. 2). Obviously, the radial displacement is anti-symmetric and the axial displacement is symmetric with respect to z ¼ 0 plane. Therefore, rzz ¼ 0 holds for the entire z = 0 plane. According to symmetry, only upper part of cylinder is taken into analysis. The mixed boundary conditions on z ¼ 0 plane are:

Tðr; z ¼ 0; tÞ ¼ 0;

sk0 ðsÞAa ðsÞJ 0 ðrsÞds 

0

3. Thermally insulated crack

qz ðr; z ¼ 0þ ; tÞ ¼ q0 HðtÞ;

1

1

Ra 1 þ p q0 ; k p

0 6 r < 1;

Aa ðsÞJ 0 ðrsÞds ¼ 0;

ð14aÞ

1 < r 6 rb ;

Ba ðsÞI0 ðr b sk0 ðsÞÞ sinðszÞds ¼ 

0

ð14bÞ

Z 0

1

Aa ðnÞJ 0 ðrb nÞenk0 ðnÞz dn:

Pursuant to Eq. (14b), Aa ðsÞ can be expressed in terms of a new unknown function Ua ðnÞ as:

Aa ðsÞ ¼

Z

1

Ua ðnÞ sinðsnÞdn:

ð16Þ

0

Substituting Eq. (16) and using some integral identities [45], the Fourier inversion of Eq. (15) results in Ba ðsÞ in terms of Ua ðnÞ:

Ba ðsÞ ¼ 

2

p

K 0 ðr b sk0 ðsÞÞ k0 ðsÞI0 ðrb sk0 ðsÞÞ

Z

1

Ua ðgÞ sinhðgsk0 ðsÞÞdg:

ð17Þ

0

Then Eq. (14a) can be re-written as an Abel integral equation:

Z 0

r

0

½Ua ðnÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn þ r2  n2

Z

1

Ua ðgÞKðr; gÞdg ¼

0

q0 Ra ð1 þ pÞ ; kp

0 6 r < 1; ð18Þ

and

rrz ðr; z ¼ 0; tÞ ¼ 0; 0 6 r < 1;

ð12aÞ

where

uðr; z ¼ 0; tÞ ¼ 0;

ð12bÞ

Kðr; gÞ ¼

1 < r 6 rb :

Z

1

s½k0 ðsÞ  1 sinðsgÞJ 0 ðrsÞds þ

0

where H(x) is the Heaviside function. The boundary conditions on the flank of the cylinder are Tðr ¼ rb ; z; tÞ ¼ 0, rrr ðr ¼ r b ; z; tÞ ¼ 0 and rzr ðr ¼ r b ; z; tÞ ¼ 0. 3.1. The temperature field

Z

1 0



Aa ðsÞJ 0 ðrsÞesk0 ðsÞz ds þ

Z

1

Ba ðsÞI0 ðrsk0 ðsÞÞ ð13Þ

Z

p

1

s

0

K 0 ðrb sk0 ðsÞÞ sinhðgsk0 ðsÞÞ I0 ðrsk0 ðsÞÞds: k0 ðsÞI0 ðrb sk0 ðsÞÞ

ð19Þ

1

Hðn; gÞUa ðgÞdg ¼

0

2 q0 Ra ð1 þ pÞ n; kp

0 6 n < 1;

p

ð20Þ

where

2

p 

Z

1

½k0 ðsÞ  1 sinðsnÞ sinðsgÞds þ

0

K 0 ðr b sk0 ðsÞÞ k20 ðsÞI0 ðr b sk0 ðsÞÞ

4

p2

Z

1 0

sinhðnsk0 ðsÞÞ sinhðgsk0 ðsÞÞds:

ð21Þ

Numerical solution for Eq. (20) can be obtained by an appropriate collocation in n. Once Ua ðnÞ is determined, the functions Aa ðsÞ and Ba ðsÞ can be evaluated to result in the temperature field in Laplace transform domain.

R

Z

2 Ra q0

Z

Hðn; gÞ ¼

0

 sinðszÞds;

2

The solution of Eq. (18) is:

Ua ðnÞ þ

By applying Hankel transform, a solution for Eq. (7) for the temperature field in Laplace transform domain can be obtained as:

T  ðr; z; pÞ ¼

ð15Þ

2 Rb

q0

Fig. 2. A cylinder with a thermally insulated crack subject to a prescribed thermal flow.

3.2. The elasticity field Thanks for the temperature field in Laplace transform, the complete solution of governing Eq. (10) in Laplace transform consists of the homogeneous solution and a particular solution. The homogeneous solution, which is independent of temperature, can be found through Hankel transform with respect to variable r. It is given as:

238



uh

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245



  a1m J 1 ðrsÞ skm z e ds F am ðsÞ a2m J 0 ðrsÞ 0 m¼1   Z 1X 2 a1m I1 ðrs=km Þ sinðszÞ þ Ra ds; Eam ðsÞ a2m I0 ðrs=km Þ cosðszÞ 0 m¼1 Z

¼ Ra

wh

1

2 X

Condition (12b) gives:

Z ð22Þ

where the subscript h means homogeneous solution, Jn and In (n = 0,1) are nth Bessel function of first kind and nth modified Bessel function of first kind, respectively, F am ðsÞ and Eam ðsÞ (m = 1, 2) are unknown functions. The characteristic value km is found to obey:

 ðc13 þ c44 Þkm   ¼ 0; c33 k2  c44 

  c  c k2 44 m  11   ðc13 þ c44 Þkm

ð23Þ

1

a1m

(

 ¼

a2m

)

1 2

c44 km  ðcc11 13 þc 44 Þkm

up

)

ð24Þ

:

wp

  b1 ðsÞJ 1 ðrsÞ sk0 ðsÞz 1 a ¼ Ra e ds A ðsÞ s b2 ðsÞJ 0 ðrsÞ 0   Z 1 b3 ðsÞI1 ðrsk0 ðsÞÞ sinðszÞ 1 a þ Ra ds; B ðsÞ s b4 ðsÞI0 ðrsk0 ðsÞÞ cosðszÞ 0 Z

1

ð25Þ

where the subscript p means particular solution, bj(s) (j = 1, 2, 3, 4) can be calculated by:

"

c11  c44 k20 ðsÞ

ðc13 þ c44 Þk0 ðsÞ

ðc13 þ c44 Þk0 ðsÞ

c44  c33 k20 ðsÞ

#

b1 ðsÞ





c11 k20 ðsÞ  c44

ðc13 þ c44 Þk0 ðsÞ

ðc13 þ c44 Þk0 ðsÞ

c44 k20 ðsÞ  c33

¼

b2 ðsÞ

v11 ; v33 k0 ðsÞ

#

b3 ðsÞ b4 ðsÞ



 ¼



v11 k0 ðsÞ : v33

2 X

rzz ¼ > :  > rrz ;

Z 0

1

½g 1m ðr; sÞeskm z 

¼

2 X

½h1m ðr; sÞ sinðszÞ

p11 ðr; sÞesk0 ðsÞz

6 4 ½g 2m ðr; sÞeskm z  ½h2m ðr; sÞ sinðszÞ p21 ðr; sÞesk0 ðsÞz ½g 3m ðr; sÞeskm z  ½h3m ðr; sÞ cosðszÞ p31 ðr; sÞesk0 ðsÞz

m¼1

F am ðsÞa1m

ð30Þ

rzz ¼ 0 on the entire z = 0 plane, one

1 C 2m F bm ðsÞ ¼  Ab ðsÞD21 ðsÞ: s

ð31Þ

1 þ Aa ðsÞb1 ðsÞ: s

1 a a F am ðsÞ ¼ F a0 ðsÞf m1  Aa ðsÞf m2 ðsÞ; s

ð28Þ

ð32Þ

where (

a

f 11

)

 ¼

a f 21

1

1

C 21

1   1 0

C 22

( ;

a

f 12 ðsÞ

)

 ¼

a f 22 ðsÞ

1 

1

1

C 21

C 22



b1 ðsÞ D21 ðsÞ

;

ð33Þ

where the functions Dim ðsÞ and the contents C im are listed in Appendix A. Secondly, the relationship between Eam ðsÞ and F am ðsÞ can be determined by the Fourier inversion of the mechanical conditions on the flank of the cylinder. That means Eam ðsÞ can be expressed in terms of Xa ðnÞ, as follows:

Eam ðsÞ ¼

2

p

Z

1

pffiffiffi a 2 xX ðxÞbam1 ðs; xÞdx þ

Z

p

0 bam3 ðsÞBa ðsÞ;

1

0

bam2 ðs; nÞAa ðnÞdn ð34Þ

where the function bami is listed in Appendix B. Finally, substituting Eqs. (30), (32) and (34) into the boundary condition (12a), an equation for Xa ðnÞ is obtained as:

l

where the functions gim(r,s), him(r,s), pi1(r,s), pi2(r,s) are listed in Appendix A. The unknown functions F am ðsÞ and Eam ðsÞ can be determined by mechanical boundary conditions. Firstly, a new unknown function Xa ðnÞ is introduced, and the relationship between F am ðsÞ and Xa ðnÞ can be determined by mechanical boundary conditions. Define:

F a0 ðsÞ

pffiffiffi a xX ðxÞJ 3=2 ðsxÞdx:

Using Eqs. (28) and (31), the relationship between F am ðsÞ and X ðnÞ can be obtained:

a

By substituting the homogeneous solution and particular solution into the constitutive relations in Eq. (9), the following thermal stresses can be yield:

2

1 0



ð26bÞ

8  9 > < rrr > =

Z

Recalling the condition obtains:



ð26aÞ "

ð29Þ

a

Obviously, Eq. (23) has four solutions for km. If we arrange the order of the roots of km such that Re(k1) < Re(k2) < Re(k3) < Re(k4), then it must have k1 = k4 and k2 = k3. Considering the convergence of displacements at z = 1, only k1 and k2 (which have negative real parts) will be used. A particular solution to Eq. (10) has been verified as:

(

pffiffi s

F a0 ðsÞ ¼

m¼1



1 < r 6 rb :

According to Eq. (29), the function F a0 ðsÞ can be expressed in terms of a new function Xa ðnÞ:

m

and the constants a1m and a2m are:

F a0 ðsÞJ1 ðrsÞds ¼ 0;

0

Z

0

½nXa ðnÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn þ 0 r 2  n2 rffiffiffiffi r

¼

p 2

rP a ;

rffiffiffiffi Z 2 r

p

1

pffiffiffi

gXa ðgÞ

Z

0

1

0

2 X

h3m ðr; sÞbam1 ðs; gÞdsdg

m¼1

0 6 r < 1;

where

8 9 3> fF am ðsÞg > > > p12 ðr; sÞ sinðszÞ > > < a = 7 fEm ðsÞg ds p12 ðr; sÞ sinðszÞ 5 a > A ðsÞ > > > > > ; p12 ðr; sÞ cosðszÞ : a B ðsÞ

Pa ðrÞ ¼

Z

! 2 X 2 a C 3m f m2 ðsÞ þ D31 ðsÞ J 1 ðrsÞAa ðsÞds 

1 0



p

m¼1

Z

1

0

þ

ð27Þ

Z

1 0

2 X m¼1 2 X m¼1

h3m ðr; sÞ

Z 0

1

bam2 ðs; nÞAa ðnÞdn ds

h3m ðr; sÞbam3 ðsÞBðsÞds 

Z 0

1

p32 ðr; sÞBa ðsÞds;

ð36Þ

239

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

P a and la ¼ 2m¼1 C 3m f m1 is a constant. Eq. (35) is an Abel integral equation and has a solution is

la Xa ðnÞ þ

Z

1

Ga ðn; gÞXa ðgÞdg ¼ Q a ðnÞ;

ð37Þ

This section solves the insulated crack problem. The case of heated crack will be analyzed in next section. 4. The heated crack

0

where

Z 2 pffiffiffiffiffiffi Ga ðn; gÞ ¼ ng

p

Q a ðnÞ ¼

1

0

2 pffiffiffi n

Z

p

1

! 2 X C 3m pffiffiffiffiffiffi s1=2 I3=2 ðns=km Þbam1 ðs; gÞ ds; km m¼1

Ua ðxÞPa ðn; xÞdx;

ð38Þ

ð39Þ

0

and where

Pa ðn; gÞ ¼ ððp=2Þ3=2 ma n3=2 gÞHðn  gÞ þ

p 2 Z

Z

1 0

ðMa ðsÞ  ma Þs1=2 J3=2 ðnsÞ sinðsgÞds

2 X C 3m pffiffiffiffiffiffi bam3 ðsÞs1=2 I3=2 ðns=km Þ  km 0 m¼1 ! D32 ðsÞ K 0 ðr b sk0 ðsÞÞ sinhðgsk0 ðsÞÞ  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I3=2 ðnsk0 ðsÞÞ ds k0 ðsÞI0 ðr b sk0 ðsÞÞ ðsk0 ðsÞÞ Z 1 Z 1X 2 C 3m pffiffiffiffiffiffi  bam2 ðs; fÞ sinðfgÞdf s1=2 I3=2 ðns=km Þds km 0 0 m¼1

qz ðr; z ¼ 0; tÞ ¼ q0 HðtÞ;

rffiffiffiffi 2

qz ðr; z ¼ 0; tÞ ¼ 0;

pffiffiffiffiffi

la Xa ð1Þ Ra : p

ð40Þ

ð41Þ

to be K II0 ¼ 2ma R3=2 a q0 =ð3pkÞ, which is consistent with the previous result [46]. Although the solutions in Laplace transform domain have been obtained, numerical Laplace inversion can not be avoided to obtain the final solutions in the time domain. The Riemann sum approximation of the Fourier integral transformed from the Laplace inversion integral is accurate enough and has been extensively used in evaluating the dynamic responses associated with non-Fourier heat conduction [1]. In this paper, a more efficient method, Stehfest method [47] is adopted. This method has reasonable accuracy for a fairly wide range of Laplace transform [48]. It gives

 2N lnð2Þ X n lnð2Þ ; CnF t n¼1 t

C n ¼ ð1ÞnþN

ð44aÞ

1 < r 6 rb ;

ð44bÞ

rzz ðr; z ¼ 0; tÞ ¼ 0; 0 6 r < 1;

ð45aÞ

wðr; z ¼ 0; tÞ ¼ 0;

ð45bÞ

and

If the diameter of the cylinder is large enough comparing to the internal crack, the cylinder can be simplified as an infinite media. Hence, the second term of Ha ðn; gÞ and the function bami varnish, and the key equations (20) and (37) become the same as in literature [38]. In addition, if the thermal load is assumed to be steady, which means the function M a ðsÞ becomes the constant ma . Then the steady stress intensity factor for the cracked infinite media is found

FðtÞ ¼

0 6 r < 1;

1

P a in which M a ðsÞ ¼ 2m¼1 C 3m f m2 ðsÞ þ D31 ðsÞ. It can be seen that the a function M ðsÞ converges to a content when s ¼ 1, which is denoted as ma , which is associated with the properties of material. Numerical solution for Eq. (37) can be obtained by an appropriate collocation in n. Once Xa ðnÞ is determined, the functions F am ðsÞ and Eam ðsÞ can be evaluated so that the elasticity fields are solved. The stress intensity factor at the crack front is of particular interest. The mechanical boundary conditions on crack surface suggest that the crack in this problem is a mode II one. Define the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi stress intensity factor as K II ¼ limR!Ra 2ðR  Ra Þrzr ðr; 0Þ, then the Laplace transform of K II is given as:

K II ¼ 

The internal thermal source is ignored in this paper. However, in some cases, the thermal generation by a heat source is strong enough that gives raise to significant stress which may products defect or crack. Therefore, the response due to the action of a heated crack in a cylinder (Fig. 3) is investigated in this section. In such case, the temperature and radial displacement are symmetric with respect to z ¼ 0 plane, whereas the axial displacement is anti-symmetric. Therefore, rrz ¼ 0 holds for the entire z ¼ 0 plane. According to symmetry, only upper part of cylinder is taken into analysis. The boundary conditions on z ¼ 0 plane are:

1 < r 6 rb

The boundary conditions on the flank of the cylinder are the same as those given in the previous section. 4.1. The temperature field Considering the symmetry of temperature field, another solution for Eq. (7) for the temperature field in Laplace transform domain is given as:

T  ðr; z; pÞ ¼

Z

1

Ab ðsÞJ 0 ðrsÞesk0 ðsÞz ds þ

Z

0

1

Bb ðsÞI0 ðrsk0 ðsÞÞ cosðszÞds;

0

ð46Þ where Ab ðsÞ and Bb ðsÞ are unknown functions to be determined. The thermal flux qz ðr; z; pÞ associated with Eq. (46) can be calculated from the Laplace transform of Eq. (6). The unknown functions Ab ðsÞ and Bb ðsÞ will be determined by the thermal boundary conditions of the cylinder. It follows from the thermal boundary conditions on crack surface (Eq. (44)) and on the flank of cylinder (Tðr ¼ rb ; z; tÞ ¼ 0) that

Z

1

sk0 ðsÞAb ðsÞJ 0 ðrsÞds ¼

0

Z

1

Ra 1 þ p q0 ; k p

sk0 ðsÞAb ðsÞJ 0 ðrsÞds ¼ 0;

0 6 r < 1;

ð47aÞ

1 < r 6 rb ;

ð47bÞ

0

Z

1

Bb ðsÞI0 ðr b sk0 ðsÞÞ cosðszÞds ¼ 

0

Z 0

1

Ab ðnÞJ 0 ðrb nÞenk0 ðnÞz dn:

ð48Þ

Pursuant to Eq. (47b), Ab ðsÞ can be expressed in terms of a new unknown function Ub ðnÞ as:

R ð42Þ

Z

2 Ra q0

minðn;NÞ X

mN ð2mÞ! : ðN  mÞ!m!ðm  1Þ!ðn  mÞ!ð2m  nÞ! m¼int½ðnþ1Þ=2 ð43Þ

2 Rb

q0

Fig. 3. A cylinder with a heated crack subject to a prescribed thermal flow.

240

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

1 sk0 ðsÞ

Ab ðsÞ ¼

Z

1

Ub ðxÞ sinðsxÞdx:

ð49Þ

0

Then Eq. (47a) can be re-written as an Abel integral equation:

Z r

1

Ub ðxÞ Ra 1 þ p pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dx ¼ q0 ; 2 2 k p x r

Ub ðxÞ ¼  ¼

2 Ra 1 þ p d p dx

p k

x

1

ð50Þ

nq0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn n2  x2

1 < r 6 rb :

ð58Þ

term of new function Xb ðnÞ:

Z

1

Xb ðxÞ sinðsxÞdx:

ð59Þ

0

Thanks for the condition have:

2 Ra 1 þ p x q0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi : p 1  x2

ð51Þ

p k

Then the function A ðsÞ can be expressed as:

2 X

C 3m F bm ðsÞ ¼

m¼1

Ra 1 þ p J ðsÞ : q0 1 sk0 ðsÞ k p

ð52Þ

Substitution of Eq. (52) into the Fourier inversion of Eq. (48) results in:

Bb ðsÞ ¼ 

F b0 ðsÞJ0 ðrsÞds ¼ 0;

0

F b0 ðsÞ ¼

b

AðsÞ ¼

Condition (45b) gives: 1

According to Eq. (58), the function F b0 ðsÞ can be expressed in the

0 6 r < 1:

The solution of Eq. (50) is:

Z

Z

2 Ra 1 þ p I1 ðsk0 ðsÞÞK 0 ðr b sk0 ðsÞÞ : q p k p 0 sk0 ðsÞI0 ðrb sk0 ðsÞÞ

ð53Þ

rrz ¼ 0 on the entire z ¼ 0 plane, we

1 b A ðsÞD31 ðsÞ: s

ð60Þ

Using Eqs. (57) and (60), the relationship between F bm ðsÞ and

Xb ðnÞ can be obtained: 1 b b F bm ðsÞ ¼ F b0 ðsÞf m1  Ab ðsÞf m2 ðsÞ; s

ð61Þ

with (

b

f 11

) ¼

b

f 21

4.2. The elasticity field



a21

a22

C 31

C 32

1   1 0

( ;

b

f 12 ðsÞ

)

 ¼

b

f 22 ðsÞ

a21

a22

C 31

C 32

1 

b2 ðsÞ D31 ðsÞ

 ;

ð62Þ

Considering the symmetry of radial and axial displacements, the homogeneous solution of governing Eq. (10) is given as:



uh



wh

  2 1X a1m J 1 ðrsÞ skm z ¼ Ra e ds F bm ðsÞ a2m J 0 ðrsÞ 0 m¼1   Z 1X 2 a1m I1 ðrs=km Þ cosðszÞ ds; Ebm ðsÞ þ Ra a2m I0 ðrs=km Þ sinðszÞ 0 m¼1 Z

F bm ðsÞ

flank of the cylinder. That means Ebm ðsÞ can be expressed in terms

ð54Þ

Ebm ðsÞ

where and are unknown functions, the characteristic value km and constants a1m and a2m have been given by Eq. (23). A particular solution has been verified as:

(

up

)

wp

  b1 ðsÞJ 1 ðrsÞ sk0 ðsÞz 1 b ¼ Ra e A ðsÞ ds s b2 ðsÞJ 0 ðrsÞ 0   Z 1 b3 ðsÞI1 ðrsk0 ðsÞÞ cosðszÞ 1 b ds; þ Ra B ðsÞ s b4 ðsÞI0 ðrsk0 ðsÞÞ sinðszÞ 0 Z

rzz ¼ > :  > rrz ;

Z 0

2 1

½g 1m ðr; sÞeskm z 

½h1m ðr; sÞ cosðszÞ

lb

p11 ðr; sÞesk0 ðsÞz

6 4 ½g 2m ðr; sÞeskm z  ½h2m ðr; sÞ cosðszÞ p21 ðr; sÞesk0 ðsÞz ½g 3m ðr; sÞeskm z  ½h3m ðr; sÞ sinðszÞ p31 ðr; sÞesk0 ðsÞz

function Xb ðnÞ is introduced, and the relationship between F bm ðsÞ and Xb ðnÞ can be determined by the mechanical boundary conditions. Define:

1 a2m F bm ðsÞ þ Ab ðsÞb2 ðsÞ: s m¼1

Z

2

p

1

Xb ðxÞbbm1 ðs; xÞdx 

0

2

p

Z

1

0

bbm2 ðs; nÞAb ðnÞdn

 bbm3 ðsÞBb ðsÞ;

ð55Þ

The unknown functions F bm ðsÞ and Ebm ðsÞ can be determined by mechanical boundary conditions. Firstly, a new unknown

2 X

Ebm ðsÞ ¼

ð63Þ

about Xb ðnÞ is obtained and equivalent to:

where the functions gim(r,s), him(r,s), pi1(r,s), pi2(r,s) are listed in Appendix A.

F b0 ðsÞ ¼

of Xb ðnÞ, as follows:

where the function bbmi is listed in Appendix C. Finally, substituting Eqs. (59), (61) and (63) into boundary condition (45a), a equation

1

where the functions bj(s) have been given by Eq. (26). By substituting the homogeneous solution and particular solution into the constitutive relations in Eq. (9), the following thermal stresses can be yield:

8  9 > < rrr > =

where the functions Dim ðsÞ and the contents C im are listed in Appendix A. Secondly, the relationship between Ebm ðsÞ and F bm ðs can be determined by the Fourier inversion of conditions on the

ð57Þ

Z 0

r

0

½Xb ðnÞ 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn þ 2 2 p r n

¼ Pb ðrÞ;

Z

1

Xb ðgÞ

0

Z 0

1

2 X C 2m m¼1

km

sI0 ðrs=km Þbbm1 ðs; gÞdsdg

0 6 r < 1;

ð64Þ

with

8 b 9 fF m ðsÞg > 3> > > > > > > < b = 7 fEm ðsÞg ds p12 ðr; sÞ cosðszÞ 5 b > > > A ðsÞ > > > > > p12 ðr; sÞ sinðszÞ : ; Bb ðsÞ p12 ðr; sÞ cosðszÞ

Pb ðrÞ ¼

Z

! 2 X 2 b C 2m f m2 ðsÞ  D21 ðsÞ J 0 ðrsÞAb ðsÞds þ

1 0





p

m¼1

Z

1

0

þ

ð56Þ

Z Z

1 0

m¼1 2 X

h2m ðr; sÞ

Z 0

1

bbm2 ðs; fÞAb ðfÞdfds

h2m ðr; sÞbbm3 ðsÞBb ðsÞds

m¼1 1

0

2 X

p22 ðr; sÞI0 ðrsk0 ðsÞÞBb ðsÞds;

ð65Þ

241

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

P b where lb ¼ 2m¼1 C 2m f m1 is a constant. Eq. (64) is an Abel integral equation and has a solution:

lb Xb ðnÞ þ

Z

1

Gb ðn; gÞXb ðgÞdg ¼ Q b ðnÞ;

ð66Þ

0

where

Gb ðn; gÞ ¼

Q b ðnÞ ¼

Z

4

p

2

1

0

2 X

C 2m bbm1 ðs; gÞ sinhðns=km Þds;

ð67Þ

m¼1

4 Ra 1 þ p q0 Pb ðnÞ; p

ð68Þ

p2 k

and where

p

Pb ðnÞ ¼

4 þ 



qffiffiffiffiffiffiffiffiffiffiffiffiffiffi

insulated crack and heated crack, respectively). To systematically investigate the effects of thermal relaxation time and geometry size on the stress intensity factor, the radius of the cylinder is given as Rb ¼ 1 m, and a reference thermal relaxation time sq0 which satpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi isfies l0ref ¼ ksq0 =ðqcÞ ¼ Rb is introduced to normalize time t. In addition, constants dI and dII which will be exactly given in the following are used to normalize the thermal stress intensity factors K I and K II respectively. The results in Figs. 4 and 5 show the satisfactory convergence of the object with N in Stehfest Laplace inversion method. It can be seen that the object becomes stable for N P 6. Well, the great efficiency of Stehfest method for Laplace inversion in the crack problem associated with hyperbolic heat conduction is manifested again [38,39]. Therefore, N = 6 in the numerical Laplace inversion has been adopted in all following calculations.

mb arcsinðnÞ þ n 1  n2

p 2 Z 0

Z

1

0 1

 M b ðsÞ  mb s2 sinðsnÞJ 1 ðsÞds

" 2 X

5.1. The insulated crack In this case, a cracked cylinder subjected a sudden uniform thermal flux q0 in the positive Z direction on the crack faces is

C 2m bbm3 ðsÞ sinhðns=km Þ

m¼1

 sinhðnsk0 ðsÞÞ I1 ðsk0 ðsÞÞK 0 ðrb sk0 ðsÞÞ ds sk0 ðsÞ sk0 ðsÞI0 ðr b sk0 ðsÞÞ Z 1 Z 1X 2 J ðfÞ þ C 2m bbm2 ðs; fÞ 1 df sinhðns=km Þds; ð69Þ fk0 ðfÞ 0 0 m¼1

 D22 ðsÞ

P b in which M b ðsÞ ¼ 2m¼1 C 2m f m2 ðsÞ  D21 ðsÞ. The function M b ðsÞ converges to a content when s ¼ 1, which is denoted as mb . The content mb is associated with the properties of material. Numerical solution for Eq. (66) can be obtained by an appropriate collocation in n. Once

0.25 Ra =0.25Rb

N=2 N=4 N=6 N=8

l0 =0.25Rb

0.20

Xb ðnÞ is determined, the functions F bm ðsÞ and Ebm ðsÞ can be evaluated so that the elasticity fields come out. The stress intensity factor at the crack front is of particular interest. The mechanical boundary conditions on crack surface suggest that the crack in this problem is a mode I one. The stress pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi intensity factor is defined as K I ¼ limR!Ra 2ðR  Ra Þrzz ðr; 0Þ. Then the Laplace transform of K I is given as:

K I

b

¼ l

pffiffiffiffiffi X ð1Þ Ra : b

0.15

0.10 0.0

0.1

0.2

ð70Þ

If the diameter of the cylinder is large enough comparing to the internal crack, the cylinder can be simplified as an infinite media. Hence, the functions Bb ðsÞ and bbmi varnish, the solution for Ab ðsÞ and the key equation (37) become the same as in literature [38]. In addition, if the thermal load is assumed to be steady, which b

0.3

0.4

0.5

Fig. 4. Convergence for the normalized stress intensity factor with numerical Laplace inversion parameter (insulated crack).

0.12

b

means the function M ðsÞ becomes the constant m . Then the steady stress intensity factor for the cracked infinite media is found to be K I0 ¼ mb R3=2 a q0 =ð2kÞ, which is consistent with the previous result [46]. With the solution in Laplace transform, the final solution can be obtained by numerical Laplace inversion.

t /τq0

Ra =0.25Rb

N=2 N=4 N=4 N=6

l0 =0.25Rb

0.11

0.10

5. Numerical results and discussions Previous sections established a general non-Fourier model for the thermomechanical coupling. The model developed in this paper applies to any materials. For illustration purpose, material PZT-5H is considered in the numerical analysis but only the thermal and mechanical properties of the material are used. The electrical properties of the materials do not influence the numerical results of this paper. The material properties cij and vii and associative contents la, lb, ma and mb are listed in Appendix D. Numerical results of thermal stress intensity factor are obtained for the cylinder with transverse isotropy in two cases (thermal

0.09

0.08 0.0

0.1

0.2

0.3

t /τq0

0.4

0.5

0.6

0.7

Fig. 5. Convergence for the normalized stress intensity factor with numerical Laplace inversion parameter (heated crack).

242

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

0.40

0.7 l 0 =0.4Rb , H

Ra =0.25Rb

l 0 =0.4Rb

l 0 =0.3Rb , H

0.35

0.6

l 0 =0.3Rb

l 0 =0.2Rb , H l 0 =0.1Rb , H

0.30

l 0 =0.2Rb

0.5

Steady State

l0 is an arbitrary value, F

0.25

0.4

0.20

0.3

0.15

0.2

0.10

0.1

0.05 0.0

0.1

0.2

t /τq0

0.3

0.4

0.5

0.0

0.10

0.20

0.30

0.40

0.50

0.60

Ra / Rb

Fig. 6. Normalized stress intensity factors for selected parameter l0 (insulated crack). Lines noted by H and noted by F are results from hyperbolic heat conduction and Fourier’s law, respectively.

Fig. 8. Amplitudes of stress intensity factor for selected parameter l0 (insulated crack).

0.18

0.35 l 0 =0.25Rb

0.30

Ra =0.4Rb , H

l 0 =0.4Rb , H

Ra =0.25Rb

0.16

l 0 =0.3Rb , H

Ra =0.4Rb , F

l 0 =0.2Rb , H

0.14

0.25 Ra =0.3Rb , H

0.20

Ra =0.3Rb , F

0.15

Ra =0.2Rb , H

l 0 =0.1Rb , H l 0 is an arbitrary value, F

0.12 0.10 0.08

Ra =0.2Rb , F

0.10

0.06 Ra =0.1Rb , H

0.05 0.00 0.0

Ra =0.1Rb , F

0.04

0.1

0.2

t /τq0

0.3

0.4

0.5

0.02 0.0

0.1

0.2

0.3

t /τq0

0.4

0.5

0.6

0.7

Fig. 7. Normalized stress intensity factor for selected crack radium (insulated crack). Lines noted by H and noted by F are results from hyperbolic heat conduction and Fourier’s law, respectively.

Fig. 9. Normalized stress intensity factors for selected parameter l0 (heated crack). Lines noted by H and noted by F are results from hyperbolic heat conduction and Fourier’s law, respectively.

The stress intensity factor is given as K II ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi l X ð1Þ 2Ra =p, where Xa ðxÞ is found to contain a factor Ra q0 =k. The steady stress intensity factor for the infinite medium

clearly higher than that for the classical Fourier law. Especially, for long thermal relaxation time or small crack radium, the relative error between the results under two theories of heat conduction is very significant. For example, when l0 ¼ 0:25Rb and Ra ¼ 0:1Rb , the amplitude of K II under hyperbolic heat conduction is 3.5 times greater than the corresponding value from Fourier’s law. Obviously, the non-Fourier effect cannot be ignored. On the other hand, the results under two theories go to uniform over time as expected. And it can be seen that the results from hyperbolic heat conduction converge slowly for long thermal relaxation time or large crack radium. In practice, the maximum of stress intensity factor which is always responsible to the fracture of the specimen is of particular interest. Fig. 8 depicts the maximal stress intensity factors during the thermal shock under hyperbolic heat conduction. Where, the steady stress intensity factor which is independent on the thermal relaxation time is also displayed for comparison. Note that the maximal stress intensity factor from Fourier’s law is just the result under steady state. Fig. 8 clearly shows the influence of parameter

analyzed. a

a

under the same thermal load is K II0 ¼ 2ma R3=2 a q0 =ð3pkÞ. Applying dII ¼ 2ma R3=2 b q0 =ð3pkÞ to normalize K II , the effect of crack size on the stress intensity factor can be investigated clearly. Figs. 6 and 7 show the variation of the thermal stress intensity factors in time domain for selected values of parameter l0 and crack radium Ra . It can be easily read that the stress intensity factor increases with increasing l0 or Ra . Note that the parameter l0 is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dependent on the thermal relaxation time sq , l0 ¼ ksq =ðqcÞ. Thus the influence of sq on the stress intensity factors can emerge from l0 . For comparison, results for the classical Fourier heat conduction are also given. The curves of stress intensity factors based on classical Fourier law climb smoothly. Whereas the results for the hyperbolic heat conduction show considerable oscillations before converge to the corresponding steady solutions. In all cases, the stress intensity factor for the hyperbolic heat conduction is

243

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

0.18 l0 =0.25Rb

Ra =0.4Rb

Ra =0.4Rb

0.16 Ra =0.3Rb

0.14

Ra =0.3Rb

0.12

Ra =0.2Rb

0.10

Ra =0.2Rb

0.08 0.06

Ra =0.1Rb Ra =0.1Rb

0.04 0.02 0.00 0.0

0.1

0.2

0.3

t /τq0

0.4

0.5

0.6

0.7

Fig. 10. Normalized stress intensity factor for selected crack radium (heated crack). The solid lines and broken lines are results from hyperbolic heat conduction and Fourier’s law, respectively.

l 0 =0.4Rb l 0 =0.3Rb l 0 =0.2Rb

0.25

6. Conclusions This paper studies the thermal shock fracture of cracked cylinder based on hyperbolic heat conduction. By Laplace transform and dual integral equation technique, solutions for insulated crack and heated crack are established respectively. The stress intensity factors are calculated to observe the effects of the thermal relaxation time and crack size on it. It is found that the results from hyperbolic non-Fourier heat conduction is higher than ones from Fourier’s law during the beginning of the thermal shock, especially for large thermal relaxation time or small crack radium relative to cylinder radium. This means that the non-Fourier effect is considerable in highly-varying thermal shock environments or at a small scale. In addition, the significant difference between the results for given crack radium suggest that the previous research for the infinite media is not sufficient in understanding the fracture behavior of materials under thermal shock.

0.35 0.30

law are displayed too. As same as that for the insulated crack, the stress intensity factor for the hyperbolic heat conduction is higher than that for the classical Fourier law, and the results increase with increasing l0 and crack size. However, the non-Fourier effect is not as strong as it in the case of the thermally insulated crack. For example, in Fig. 10, the non-Fourier effect can be almost ignored when l0 ¼ 0:25Rb and Ra ¼ 0:4Rb . Whereas in the case of the thermally insulated crack, Fig. 7 shows that the error between the results from the non-Fourier heat conduction and Fourier’s law reaches 30% under the same settings of l0 and Ra . Fig. 11 depicts the amplitude of stress intensity factors during the thermal shock. Interestingly, with increasing crack radium, the stress intensity factor increases slowly in this case of the heated crack. Such a behavior is different from the case of the thermally insulated crack. It is speculated that the different thermal boundary conditions in z = 0 plane in two cases result this phenomenon. For the case of the thermally insulated crack, the singularity of the temperature change at the crack front is stronger than it in the case of the heated crack.

Steady State

0.20 0.15 0.10 0.05

Conflict of interest 0.00

0.10

0.20

0.30

0.40

0.50

0.60

Ra / Rb

Fig. 11. Amplitudes of stress intensity factor for selected parameter l0 (heated crack).

l0 and crack radium Ra on the stress intensity factors again. The larger ratio of parameter l0 to crack radium Ra , the more prominent non-Fourier effect on results. And the stress intensity factors increase rapidly with crack radium.

None declared. Acknowledgments The authors thank the National Natural Science Foundation of China (project nos. 11172081 and 11372086) and Natural Science Foundation of Guangdong Province of China (project no. 2014A030313696) for the support of their research. Appendix A

5.2. The heated crack In this case, a cracked cylinder subjected a sudden uniform thermal flux q0 in the positive Z direction on the upper crack face and in the negative Z direction on the lower crack face is analyzed. pffiffiffiffiffi The stress intensity factor is given as K I ¼ lb Xb ð1Þ Ra , where

1 g 1m ðr; sÞ ¼ C 1m sJ 0 ðrsÞ  ðc11  c12 Þa1m J 1 ðrsÞ r h1m ðr; sÞ ¼

C 1m 1 sI0 ðrs=km Þ  ðc11  c12 Þa1m I1 ðrs=km Þ r km

Xb ðxÞ is found to contain a factor Ra q0 =k. The steady stress intensity factor for the infinite medium under the same thermal load is b

p11 ðr; sÞ ¼ D11 ðsÞJ 0 ðrsÞ  ðc11  c12 Þ

R3=2 a q0 =ð2kÞ.

K I0 ¼ m To conveniently investigating the effect of the crack size on the stress intensity factor, a constant 3=2

dI ¼ mb Rb q0 =ð2kÞ is used to normalize K I . Figs. 9 and 10 display the thermal stress intensity factors as a function of time under different settings of parameter l0 and crack radium Ra . For comparison, the results deduced from the Fourier’s

b1 ðsÞ J ðrsÞ rs 1

p12 ðr; sÞ ¼ D12 ðsÞI0 ðrsk0 ðsÞÞ  ðc11  c12 Þ g 2m ðr; sÞ ¼ C 2m sJ 0 ðrsÞ; h2m ðr; sÞ ¼

b3 ðsÞ I1 ðrsk0 ðsÞ rs

C 2m sI0 ðrs=km Þ km

244

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245

p21 ðr; sÞ ¼ D21 ðsÞJ 0 ðrsÞ; p22 ðr; sÞ ¼ D22 ðsÞI0 ðrsk0 ðsÞÞ

b 12 ðs; nÞ

a

b 2 X C 1m km f m2 ðnÞ

¼

m¼1

C 3m g 3m ðr; sÞ ¼ C 3m sJ 1 ðrsÞ; h2m ðr; sÞ ¼ sI1 ðrs=km Þ km



p21 ðr; sÞ ¼ D31 ðsÞJ 1 ðrsÞ; p22 ðr; sÞ ¼ D32 ðsÞI1 ðrsk0 ðsÞÞ C 1m ¼ c11 a1m þ c13 a2m km ;

¼ c11 b3 ðsÞk0 ðsÞ  c13 b4 ðsÞ  v11

ab21 ðs; xÞ ¼

b 2 X C 3m f m1 s m¼1

C 2m ¼ c13 a1m þ c33 a2m km ;

"

ab22 ðs; nÞ ¼

D21 ðsÞ ¼ c13 b1 ðsÞ  c33 b2 ðsÞk0 ðsÞ  v33 ; D22 ðsÞ

D31 ðsÞ ¼ c44 ½b1 ðsÞk0 ðsÞ þ b2 ðsÞ; D32 ðsÞ ¼ c44 ½b3 ðsÞ þ b4 ðsÞk0 ðsÞ

fbam1 ðs; xÞg fbam2 ðs; nÞg fbam3 ðsÞg  1  a  a11 ðs; xÞ aa12 ðs; nÞ aa13 ðsÞ ½h1m ðrb ; sÞ ¼ a a a a21 ðs; xÞ a22 ðs; nÞ a23 ðsÞ ; ½h3m ðrb ; sÞ ¼

" a 2 X C 1m f m1 s3=2 m¼1

þ "

aa12 ðs; nÞ ¼

ðkm Þ5=2

 #  a 2 ðc11  c12 Þ X a1m f m1 s1=2 s s ; K  r  x I 1 3=2 b 3=2 rb km km m¼1 ðkm Þ #

a 2 X C 1m f m2 ðnÞ

D11 ðnÞ sJ 0 ðr b nÞ  2 2 2 2 k þ s þ s2 k20 ðsÞ n n m¼1 m " # a 2 ðc11  c12 Þ X a1m f m2 ðnÞ b1 ðnÞ sJ1 ðrb nÞ  2 ;  2 2 2 rb n n þ s2 k20 ðsÞ m¼1 n km þ s

"

2 X a ¼ C 3m f m2 ðnÞ m¼1

km n n2 k2m þ s2

 D31 ðnÞ

k0 ðnÞn n2 þ s2 k20 ðsÞ

aa23 ðsÞ ¼ D32 ðsÞI1 ðrb sk0 ðsÞÞ: Appendix C



D31 ðnÞs

#

n2 þ s2 k20 ðsÞ

J 1 ðr b nÞ;

Appendix D

fbbm1 ðs; xÞg fbbm2 ðs; nÞg fbbm3 ðsÞg  1  b  ½h1m ðrb ; sÞ a11 ðs; xÞ ab12 ðs; nÞ ab13 ðsÞ ; ¼ b b b ½h3m ðrb ; sÞ a21 ðs; xÞ a22 ðs; nÞ a23 ðsÞ "  b 2 X C 1m f m1 s s ab11 ðs; xÞ ¼ K 0  rb 2 km km m¼1  #  b ðc11  c12 Þ a1m f m1 s s sinh  x ; K 1  rb  rb km km km

c11 ¼ 12:6  1010 N=m2 ;

c12 ¼ 7:95  1010 N=m2 ;

c13 ¼ 8:41  1010 N=m2 ;

c44 ¼ 2:30  1010 N=m2 ;

The temperature-stress coefficients vii can be calculated by vii ¼ cii a, where a is the thermal expansion coefficient. Referring to the property of ceramic, set a ¼ 0:75  105 K1 . Thus vii are:

v11 ¼ 0:945  106 N=ðm2 KÞ; v33 ¼ 0:878  106 N=ðm2 KÞ: Following the material properties cij and vii , the contents la , ma ,

lb and mb can be calculated by computer:

la ¼ 3:51  1010 N=ðm3 KÞ; ma ¼ 1:62  105 N=ðm3 KÞ; lb ¼ 3:38  1010 N=ðm3 KÞ; mb ¼ 1:44  105 N=ðm3 KÞ: References

ðc11  c12 Þ b3 ðsÞ I1 ðrb sk0 ðsÞÞ; rb s

" # 1=2   a a 2 X C 3m f m1 C 3m f m1 s2 s s s a  a21 ðs; xÞ ¼  K 1  r b I3=2  x ; 3 km km km km km m¼1

a

n2 k2m þ s2

þ

c33 ¼ 11:7  1010 N=m2 :

 s K 0  rb km

aa13 ðsÞ ¼ D12 ðsÞI0 ðrb sk0 ðsÞÞ 

a 22 ðs; nÞ

k2m

ðc11  c12 Þ b3 ðsÞ I1 ðr b sk0 ðsÞÞ; rb s

The elastic contents cij for PZT-5H are [49]:

Appendix B

a

nJ 0 ðrb nÞ

ab23 ðsÞ ¼ D32 ðsÞI1 ðrb sk0 ðsÞÞ:

C 3m ¼ c44 ða1m km  a2m Þ;

a 11 ðs; xÞ

n2 þ s2 k20 ðsÞ

  s s K 1  rb sinh  x ; km km

b 2 X C 3m f m2 ðnÞs

m¼1

¼ c13 b3 ðsÞk0 ðsÞ  c33 b4 ðsÞ  v33

þ

!

! b 2 ðc11  c12 Þ X a1m km f m2 ðnÞ b1 ðnÞk0 ðnÞ þ J 1 ðr b nÞ; 2 2 2 rb n2 þ s2 k20 ðsÞ m¼1 n km þ s

ab13 ðsÞ ¼ D12 ðsÞI0 ðrb sk0 ðsÞÞ 

D11 ðsÞ ¼ c11 b1 ðsÞ  c13 b2 ðsÞk0 ðsÞ  v11 ; D12 ðsÞ



n2 k2m þ s2

D11 ðnÞk0 ðnÞ

# J 1 ðr b nÞ;

[1] D.Y. Tzou, The generalized lagging response in small-scale and high-rate heating, Int. J. Heat Mass Transfer 38 (1995) 3231–3240. [2] C. Cattaneo, Sur une forme de l’equation de la chaleur eliminant le paradoxe d’ine propagation instantanee, C.R. Acad. Sci. 247 (1958) 431–433. [3] P. Vernotte, Les paradoxes de la theorie continue de l’equation de la chaleur, C.R. Acad. Sci. 246 (1958) 3154–3155. [4] T.Q. Qiu, C.L. Tien, Short-pulse laser heating on metals, Int. J. Heat Mass Transfer 35 (3) (1992) 719–726. [5] T.Q. Qiu, C.L. Tien, Heat transfer mechanisms during short-pulse laser heating of metals, J. Heat Transfer 115 (4) (1993) 835–841. [6] K. Mitra, S. Kumar, A. Vedevarz, M.K. Moallemi, Experimental evidence of hyperbolic heat conduction in processed meat, J. Heat Transfer 117 (3) (1995) 568–573. [7] W.G. Li, D.J. Li, C.Z. Zhang, D.N. Fang, Modelling the effect of temperature and damage on the fracture strength of ultra-high temperature ceramics, Int. J. Fract. 176 (2) (2012) 181–188. [8] W.G. Li, D.Y. Li, R.Z. Wang, D.N. Fang, Numerical simulation for thermal shock resistance of thermal protection materials considering different operating environments, Sci. World J. 2013 (2013) 324186. [9] T.B. Cheng, W.G. Li, C.Z. Zhang, D.N. Fang, Unified thermal shock resistance of ultra-high temperature ceramics under different thermal environments, J. Therm. Stresses 37 (1) (2014) 14–33. [10] T.B. Cheng, W.G. Li, The temperature-dependent ideal tensile strength of ZrB2, HfB2, and TiB2, J. Am. Ceram. Soc. 98 (1) (2015) 190–196. [11] Q.H. Qin, W.C. Zhang, A new meshless method for 2D heat transfer in skin tissues with moving boundary, Univ. J. Eng. Mech. 1 (2013) 50–70. [12] K.J. Baumeister, T.D. Hamill, Hyperbolic heat-conduction equation—a solution for the semi-infinite body problem, J. Heat Transfer 91 (4) (1969) 543–548. [13] Y. Taitel, On the parabolic, hyperbolic and discrete formulation of the heat conduction equation, Int. J. Heat Mass Transfer 15 (2) (1972) 369–371.

S.L. Guo, B.L. Wang / International Journal of Heat and Mass Transfer 91 (2015) 235–245 [14] W.S. Kim, J.R. Hector, M.N. Qzisik, Hyperbolic heat conduction due to axisymmetric continuous or pulsed surface sources, J. Appl. Phys. 68 (1990) 5478–5485. [15] G. Louis, J.R. Hector, W.S. Kim, M.N. Qzisik, Propagation and reflection waves in a finite medium due to axisymmetric surface sources, Int. J. Heat Mass Transfer 35 (4) (1992) 897–912. [16] M. Torabi, S. Saedodin, Analytical and numerical solutions of hyperbolic heat conduction in cylindrical coordinates, J. Thermophys. Heat Transfer 25 (2) (2011) 239–253. [17] A. Moosaie, Non-Fourier heat conduction in a finite medium subjected to arbitrary periodic surface disturbance, Int. Commun. Heat Mass Transfer 34 (8) (2007) 996–1002. [18] A. Moosaie, Non-Fourier heat conduction in a finite medium subjected to arbitrary non-periodic surface disturbance, Int. Commun. Heat Mass Transfer 35 (3) (2008) 376–383. [19] A. Moosaie, Non-Fourier heat conduction in a finite medium with insulated boundaries and arbitrary initial conditions, Int. Commun. Heat Mass Transfer 35 (1) (2008) 103–111. [20] W.T. Zhao, J.H. Wu, Z. Chen, Analysis of non-Fourier heat conduction in a solid sphere under arbitrary surface temperature change, Arch. Appl. Mech. 84 (4) (2014) 505–518. [21] M.H. Babaei, Z.T. Chen, Transient hyperbolic heat conduction in a functionally graded hollow cylinder, J. Thermophys. Heat Transfer 24 (2) (2010) 325–330. [22] S.S. Manson, Behavior of materials under conditions of thermal stress, TN 2933: NACA, 1953. [23] W.D. Kingery, Factors affecting thermal stress resistance of ceramic materials, J. Am. Ceram. Soc. 38 (1) (1955) 3–15. [24] W.G. Li, R.Z. Wang, D.Y. Li, X.L. Shen, H.B. Kou, D.N. Fang, Effect of the cooling medium temperature on the thermal shock resistance of ceramic materials, Mater. Lett. 138 (2015) 216–218. [25] D.P.H. Hasselman, Unified theory of thermal shock fracture initiation and crack propagation in brittle ceramics, J. Am. Ceram. Soc. 52 (11) (1969) 600–604. [26] W.K. Wilson, I.W. Yu, The use of the J-integral in thermal stress crack problems, Int. J. Fract. 15 (4) (1979) 377–387. [27] N. Sumi, T. Katayama, Thermal stress singularities at tips of a Griffith crack in a finite rectangular plate, Nucl. Eng. Des. 60 (3) (1980) 389–394. [28] Y.M. Tsai, Thermal stress in a transversely isotropic medium containing a penny-shaped crack, J. Appl. Mech. 50 (1) (1983) 24–28. [29] K.Y. Lee, C.W. Shul, Determination of thermal stress intensity factors for an interface crack under vertical uniform heat flow, Eng. Fract. Mech. 40 (6) (1991) 1067–1074. [30] K.Y. Lee, S.J. Park, Thermal stress intensity factors for partially insulated interface crack under uniform heat flow, Eng. Fract. Mech. 50 (4) (1995) 475– 482. [31] H.F. Nied, F. Erdogan, Transient thermal stress problem for a circumferentially cracked hollow cylinder, J. Therm. Stresses 6 (1) (1983) 1–14.

245

[32] S.W. Yu, Q.H. Qin, Damage analysis of thermopiezoelectric properties: part Icrack tip singularities, Theor. Appl. Frac. Mech. 25 (1996) 263–277. [33] S.W. Yu, Q.H. Qin, Damage analysis of thermopiezoelectric properties: part IIeffective crack model, Theor. Appl. Frac. Mech. 25 (1996) 279–288. [34] T.J. Lu, N.A. Fleck, The thermal shock resistance of solids, Acta Mater. 46 (13) (1998) 4755–4768. [35] P. Ma, W.J. Feng, R.K.L. Su, Fracture assessment of an interface crack between two dissimilar magnetoelectroelastic materials under heat flow and magnetoelectromechanical loadings, Acta Mech. Solida Sin. 24 (5) (2011) 429–438. [36] L.C. Guo, N. Noda, Investigation methods for thermal shock crack problems of functionally graded materials-part I: analytical method, J. Therm. Stresses 37 (3) (2014) 292–324. [37] Y.Y. Zhang, L.C. Guo, N. Noda, Investigation methods for thermal shock crack problems of functionally graded materials-part II: combined analyticalnumerical method, J. Therm. Stresses 37 (3) (2014) 325–339. [38] B.L. Wang, J.C. Han, Fracture mechanics associated with non-classical heat conduction in thermoelastic media, Phys. Mech. Astron. 55 (3) (2012) 493– 504. [39] B.L. Wang, J.C. Han, Non-Fourier heat conduction in layered composite materials with an interface crack, Int. J. Eng. Sci. 55 (2012) 66–75. [40] B.L. Wang, J.C. Han, A crack in a finite medium under transient non-Fourier heat conduction, Int. J. Heat Mass Transfer 55 (17) (2012) 4631–4637. [41] D.M. Chang, B.L. Wang, Transient thermal fracture and crack growth behavior in brittle media based on non-Fourier heat conduction, Eng. Fract. Mech. 94 (2012) 29–36. [42] B.L. Wang, J.E. Li, Thermal shock resistance of solids associated with hyperbolic heat conduction theory, Proc. R. Soc. A Math. Phys. Eng. Sci. 469 (2153) (2013) 20120754. [43] Z.T. Chen, K.Q. Hu, Thermo-elastic analysis of a cracked half-plane under a thermal shock impact using the hyperbolic heat conduction theory, J. Therm. Stresses 35 (4) (2012) 342–362. [44] K. Hu, Z. Chen, Thermoelastic analysis of a partially insulated crack in a strip under thermal impact loading using the hyperbolic heat conduction theory, Int. J. Eng. Sci. 51 (2012) 144–160. [45] I.S. Gradshteyn, I.M. Ryzhik, Tables of Integrals, Series and Products, Acadeimic Press, San Diego, CA, 1965. [46] B.L. Wang, Exact electroelastic solutions for penny-shaped cracks under prescribed temperature or thermal flow, Appl. Phys. Lett. 85 (2004) 2800– 2802. [47] H. Stehfest, Algorithm 368: Numerical Inversion of Laplace Transforms, Commun. ACM 13 (1970) 47–49. [48] B. Davies, B. Martin, Numerical inversion of the laplace transform: A survey and comparison of methods, J. Comput. Phys. 33 (1979) 1–32. [49] B.L. Wang, J.C. Han, An analytical model for electrode-ceramic interaction in multilayer piezoelectric actuators, Acta. Mech. Sin. 23 (3) (2007) 199–208.