Dual-phase-lagging heat conduction and associated thermal shock fracture of sandwich composite plates

International Journal of Heat and Mass Transfer 139 (2019) 317–329

Contents lists available at ScienceDirect

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

Dual-phase-lagging heat conduction and associated thermal shock fracture of sandwich composite plates S.L. Guo a, K.F. Wang b, B.L. Wang b,c,⇑ a

School of Civil Engineering, Qingdao University of Technology, Qingdao 266033, PR China Graduate School at Shenzhen, Harbin Institute of Technology, Harbin 150001, PR China c Centre for Infrastructure Engineering, School of Computing, Engineering and Mathematics, Western Sydney University, Penrith, NSW 2751, Australia b

a r t i c l e

i n f o

Article history: Received 23 November 2018 Received in revised form 17 March 2019 Accepted 16 April 2019 Available online 14 May 2019 Keywords: Thermal shock fracture Dual-phase-lag heat conduction Sandwich composite plate Penny-shaped crack

a b s t r a c t This paper studies a penny-shaped crack at the middle of a symmetric sandwich structure and is under a sudden temperature gradient on the crack faces. The temperature field of the structure and the stress intensity factor (SIF) at the crack tip are derived based on the dual-phase-lag non-Fourier heat conduction theory. Two specific factors are identified to predict the coupling effects of the thermal property differences between the skin and the core (including the differences in thermal conductivity, thermal diffusivity, heat flux lag and temperature gradient lag) on the thermoelastic behavior. If the coefficient of thermal expansion (CoTE) of the skins is greater than that of the core, the peak value of transient SIF increases. If the CoTE of the skins is smaller than that of the core, the transient SIF is scarcely influenced. The results are useful for the designs and applications of sandwich structures work under thermal shock. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction The thermal wave behavior was observed in heat propagations at ultra-low temperature [1], in small scale and under ultra-high heat flux [2]. This is a challenge for the classical Fourier heat conduction law which describes heat conductions as diffusion behaviors. Thus non-Fourier heat conduction theories were promoted to figure the wave behaviors in these severe heating processes [3–5]. Non-Fourier heat conduction theories improve the accuracy of temperature prediction in these severe heating processes. The thermoelasticity analysis based on non-Fourier heat conduction theory may be useful for the designs and applications of thermal protection materials and structures. The widely studied non-Fourier heat conduction model is the single-phase-lag model. It was first individually formulated by Cattaneo [3] and Vernottee [4], respectively. The single-phase-lag model is also called C-V model for short. In the C-V model, the heat flux lag is introduced to approach the wave nature of heat transmission. The heat conduction law is

qðx; t þ sq Þ ¼ kc rTðx; tÞ

ð1Þ

where T is the temperature, q is the heat flux vector, x is the position vector, t is the time, kc is the thermal conductivity, sq is ⇑ Corresponding author at: Graduate School at Shenzhen, Harbin Institute of Technology, Harbin 150001, PR China. E-mail address: [email protected] (B.L. Wang). https://doi.org/10.1016/j.ijheatmasstransfer.2019.04.081 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

the thermal relaxation time (heat flux lag), and r is the differential operator. The key parameter sq is related to the collision frequency of the molecules within the energy carrier. This model has been widely applied in strong transient thermoelastic problems [6] and in many fundamental thermal shock fracture analysis of single bodies [7–10], and structures [11,12] and piezoelectric materials [13,14]. It is found that the thermal wave nature significantly intensifies the incipient transient thermoelasticity response of materials to the thermal shock. However, Tzou [5] point out that C-V model may lose too many microscope features of materials in small-scale and high-rate heating. Thus Tzou [5] presented the dual-phase-lag model (the DPL model) based on the hyperbolic two-step model [15]. The DPL model introduces the heat flux lag and temperature gradient lag that addresses more microscope features than the C-V model. The heat conduction law is

qðt þ sq Þ ¼ kc rTðt þ sT Þ

ð2Þ

where sT is the temperature gradient lag. The work [16] proved the accuracy of the DPL model in small scale and under ultra-high heat flux. Although the thermal shock fracture studies based on DPL model is scarce so far, it interests researchers in recent years. Wang and his cooperators studied the thermal shock fractures of a semiinfinite medium [17] and a thin piezoelectric layer [18] separately. Chen et al. [19] investigated the non-Fourier thermal fracture of an edge-cracked cylindrical bar. Zhang et al. [20] discussed the thermal

318

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

Nomenclature a c h1 h2 Jn kc kd l

radius of the crack [m]; specific heat capacity [Jm3K1]; half thickness of the core [m]; thickness of the upper (or bottom) skins [m]; nth order Bessel function of first kind; thermal conductivity [Wm1K1]; thermal diffusivity [m2s1]; characteristic length parameter of the material 1#, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ ð1Þ l ¼ kd sq [m];

p q q0

Laplace transform parameter heat flux vector [Wm2] constant, parameter of thermal boundary condition [Wm2] space coordinates [m] time [s] temperature [K] displacement component in r direction [m] displacement component in z direction [m] position vector [m] coefficient of thermal expansion [K1] Lame constants mass density [kgm3] thermal relaxation time (heat flux lag) [s]

r, z t T u w x

a k, G

q sq

fracture in a multilayer structure. It is demonstrated that the temperature gradient lag always enhance the thermal wave nature. Besides, a fractional heat conduction model presented by Mainardi [21] has been also used in thermal shock fracture analysis [22,23] recently. However, the physical interpretation and value of the fractional order is not clear. Above all, the DPL model is essential to understand the thermoelasticity response of materials under severe thermal shock so that it is applied in this paper. Sandwich structures are widely used in engineering [24,25], especially the thermal protection system used in space industry [26–28]. Thermomechanical analysis is critical to the designs and applications of sandwich structures. Zhao et al. [29] completed the analysis of transient thermal stress in sandwich plate with functionally graded skins. Ding et al. [30] solved the interface crack problem in layered orthotropic materials under thermomechanical loading. The previous researches about sandwich structures were based on the classic Fourier model. Under ultrahigh heat flux and in small scale, non-Fourier effect may be significant [5,17,18]. To extend the application of sandwich structures, the thermomechanical analysis based on dual-phase-lag heat conduction is necessary. This paper investigates the dual-phase-lag heat conduction and associated fracture of a symmetric sandwich structure under a sudden temperature gradient. The penny-shaped crack in the middle plane is considered. The solutions of the temperature field and displacement field are completed. The stress intensity factors at the crack tip are evaluated. The influence of the property differences between the components of the sandwich structure on the thermoelasticity behavior is discussed.

2. Formulation of the problem A symmetric sandwich structure with a penny-shaped crack in its middle plane is investigated. The half thickness of the core is h1 and the radius of the crack is a. The thicknesses of the upper and bottom skins are h2 . The geometry of the structure is depicted in

sT r

temperature gradient lag [s] stress [Pa] r differential operator superscript * Laplace transform superscript (j) material number (j = 1,2) a dimensionless radius of the crack, a ¼ a=l h1 dimensionless half thickness of the core, h1 ¼ h1 =l h2 dimensionless thickness of the upper (or bottom) skins, h2 ¼ h2 =l kc ratio of thermal conductivities of materials 1# and 2#, ð2Þ ð1Þ kc ¼ kc =kc kd ratio of thermal diffusivities of materials 1# and 2#, ð2Þ ð1Þ kd ¼ kd =kd r, z dimensionless space coordinates, r ¼ r=l, z ¼ z=l ð1Þ t dimensionless time, t ¼ t=sq ðjÞ t dimensionless phase lag of temperature gradient, ðjÞ tðjÞ ¼ sðjÞ T =sq sq ratio of heat flux lags of materials 1# and 2#, ð1Þ sq ¼ sð2Þ q =sq k, G ratio of Lame constants of materials 1# and 2#, k ¼ kð2Þ =kð1Þ , G ¼ Gð2Þ =Gð1Þ a ratio of coefficients of thermal expansion of materials 1# and 2#, a ¼ að2Þ =að1Þ

Fig. 1. Employing the Taylor series expansion, the DPL heat flux equation given in Eq. (2) is equal to [5]

!   @ 1 2 @2 @ 1 þ sq þ sq 2 qðr; z; tÞ ¼ kc 1 þ sT rTðr; z; tÞ @t 2 @t @t

ð3Þ

The second-order terms in sT and the third-order terms in sq and subsequent terms are negligible in above equation. Ignoring the inner heat source, conservation of energy thus applies

@qr @qz qr @T ¼0 þ þ þ qc @t @r @z r

ð4Þ

where q is the mass density and c is the specific heat capacity. Substituting Eq. (3) into Eq. (4), it gives the governing equation of the temperature:

! 1 2 @3 @2 @ Tðr; z; tÞ sq 3 þ sq 2 þ 2 @t @t @t !   @ @2 1 @ @2 ¼ kd 1 þ sT þ þ Tðr; z; tÞ @t @r 2 r @r @z2

ð5Þ

where kd is the thermal diffusivity. It has been proved that the DPL model with first-order Taylor series expansion in sq and sT has better agreements with the experiment data than C-V model and classic Fourier model [16,31]. The model applied in this paper is more precise than the model adopted in literature [16] due to retaining higher order (second-order) terms in sq . In addition, the model applied in this paper is identical to the hyperbolic two-step model which accurately predicted the temperatures of the electron gas and the metal lattice during short-pulse laser heating of metals [15]. For common materials like ceramics, the governing Eq. (5) based on the DPL model covers a wide scale of space and time for physical observations. Denote the displacement components in r and z directions as u and w, respectively. The normal stresses in r, h and z directions are

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

319

Fig. 1. Geometry of the cracked sandwich structure.

indicated by rrr , rhh and rzz , respectively. The shear stress is denoted as rrz . The elasticity constitutive equations are

8 rrr > > > <

9 > > > =

2

6 rhh 6 ¼6 > > 4 r zz > > > > : ; rrz

k þ 2G

k

k

k

k þ 2G

k

k

k

k þ 2G

0

0

0

0

3

07 7 7 05

G 8 9 8 9 1 @u=@r > > > > > > > > > > > > <1= < = u=r T  ð3k þ 2GÞa > > > 1> @w=@z > > > > > > > : > : ; ; 0 @u=@z þ @w=@r

Tðr; z; 0Þ ¼ 0; @Tðr; z; 0Þ=@t ¼ 0 and @ 2 Tðr; z; 0Þ=@t 2 ¼ 0 ð6Þ

where k and G are Lame constants, a is the coefficient of thermal expansion. Inertia effect can be ignored when thermal wave speed is one order of magnitude smaller than sound speed [18]. In spite of that inertia effect may increase the amplitude of the thermoelastic response, it does not change the laws of effects of non-Fourier and property differences between the structure components. For convenience, inertia effect is ignored in this paper. The balance equations of stresses are

( @r

rr

@r @ rrz @r

rrz þ @@z þ rrr r rhh ¼ 0

þ

@ rzz @z

þ

rrz r

ð7Þ

¼0

For a systematic study, we introduce the dimensionless paramð1Þ

eters according to: t ¼ t=sq , r ¼ r=l, z ¼ z=l, a ¼ a=l, h1 ¼ h1 =l, ð2Þ ð1Þ ð2Þ ð1Þ ð2Þ tð1Þ ¼ s s t ¼ sð2Þ T =sq , kc ¼ kc =kc , kd ¼ kd =kd , ð2Þ ð1Þ s k ¼ k =k , G ¼ G =G and a ¼ að2Þ =að1Þ , where

h2 ¼ h2 =l,

ð2Þ ð1Þ q = q ,

sq ¼ s l¼

ð1Þ ð1Þ ð2Þ T = q , ð2Þ ð1Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ ð1Þ kd sq is a characteristic length parameter of the material

1#, and the supper script (j) (j = 1,2) indicates the material number. Then the dimensionless governing equations of the temperature are

! !   1 @3 @2 @ @2 1 @ @2 ð1Þ ð1Þ @ þ þ þ þ T ðr; z; tÞ ¼ 1 þ t 2 @t3 @t2 @t @r 2 r @r @z2 @t T ð1Þ ðr; z; tÞ;

for material 1#;

the temperature. Therefore the initial temperature of the structure can be set as zero. However Eq. (5) contains third order derivative with respect to time which requires more initial conditions. Assume the sandwich structure is initially static. The initial thermal condition is [5,16–20]

ð8aÞ

!

1 2 @3 @2 @ T ð2Þ ðr;z;tÞ s þ sq 2 þ 2 q @t 3 @t @t !   @ @2 1 @ @2 T ð2Þ ðr;z;tÞ; formaterial2#; þ ¼ kd 1 þ tð2Þ sq þ @r 2 r @r @z2 @t ð8bÞ

ð9Þ

The temperature change is simplified as DT ¼ Tðr; z; tÞ. Applying Laplace transformation, the governing equations (8a) and (8b) and the initial conditions (9) yield the general solution of the temperature

T ðjÞ ðr; z; pÞ ¼

Z

1

  ðjÞ ðjÞ AðjÞ ðsÞec sz þ BðjÞ ðsÞec sz J 0 ðsrÞds

ð10aÞ

0

where ð1Þ

c

cð2Þ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0:5p3 þ p2 þ p ¼ 1þ ; ð1 þ tð1Þ pÞs2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u 0:5s2q p3 þ sq p2 þ p ¼ t1 þ kd ð1 þ tð2Þ sq pÞs2

ð10bÞ

and where Jn is nth order Bessel function of first kind, the superscript * indicates Laplace transformation, p is Laplace transform parameter. The unknowns AðjÞ ðsÞ and BðjÞ ðsÞ will be determined by boundary conditions. A symmetric and an anti-symmetric thermal shocks located on the crack faces will be considered separately in the following. Through an appropriate superposition, the solutions of the two cases can be used in arbitrary thermal shock problem. 3.1. Symmetric thermal shock Image the inner crack suddenly release heat. The normal derivatives of the temperatures on the crack faces are equal and opposite constants. Due to the symmetry, only the upper part of the structure is considered. The surface boundary condition and the continuity conditions are given as

@ Tðr; h1 þ h2 ; tÞ ¼ 0 @z

ð11Þ



þ

ð12Þ



þ

ð13Þ

Tðr; h1 ; tÞ ¼ Tðr; h1 ; tÞ qðr; h1 ; tÞ ¼ qðr; h1 ; tÞ The conditions on the middle plane are given as

3. Heat conduction

@ lq Tðr; 0; tÞ ¼ ð1Þ0 HðtÞ; @z kc

The solution of Eq. (5) requires initial and boundary conditions. The properties of the materials are assumed to be independent of

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

r>a

r
ð14aÞ

ð14bÞ

320

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

where q0 is a known constant, and HðtÞ is Heaviside function. Applying Laplace transformation, the conditions (14a) and (14b) are

As sac increasing, H1 ðsac tÞ approaches to Heaviside function closer and closer. The Laplace transform of the heating function is

ð1Þ

rewritten as @T  =@z ¼ lq0 =ðpkc Þ, r < a; and @T  =@z ¼ 0, r > a. Substituting the general solution of the temperature given in Eq. (10a), the conditions (14a) and (14b) result

Z

1





cð1Þ s Að1Þ ðsÞ þ Bð1Þ ðsÞ J0 ðsrÞds ¼

0

Z 0

1



lq0 ð1Þ pkc

;

r
ð15aÞ



cð1Þ s Að1Þ ðsÞ þ Bð1Þ ðsÞ J0 ðsrÞds ¼ 0; r > a

ð15bÞ

Applying inverse Hankel transform, it gives

Z a   lq0 Að1Þ ðsÞ  Bð1Þ ðsÞ cð1Þ ¼  ð1Þ nJ 0 ðsnÞdn 0 pkc

ð16Þ

Here Eq. (16) is a linear algebraic equation of the unknowns AðjÞ ðsÞ and BðjÞ ðsÞ. In addition, there are the other three linear algebraic equations of AðjÞ ðsÞ and BðjÞ ðsÞ can be deduced from the surface boundary condition and the continuity conditions. The unknowns AðjÞ ðsÞ and BðjÞ ðsÞ are solved as

fCðsÞg ¼  where



lq0  ð1Þ

pkc

va ðsÞ

the



Z

nJ 0 ðsnÞdn

Suppose the normal derivatives of the temperatures on the crack faces are the same constant. The surface boundary condition and the continuity conditions are the same as Eqs. (11)–(13). The conditions on the middle plane are given as

@ lq Tðr; 0; tÞ ¼ ð1Þ0 HðtÞ; @z kc Tðr; 0; tÞ ¼ 0;

T

is

defined

as

þ sac t þ 1Þe

sac t

ð18Þ

where sac is a dimensionless parameter associated with the heating rate. Fig. 2 draws the curves of heating function H1 ðsac tÞ for different values of sac . In the beginning of thermal loading that t 2 ½0; 10=sac , the value of H1 ðsac tÞ quickly climbs from 0 to 0.997. That simulates a general thermal shock. When t > 10=sac , it has H1 ðsac tÞ 2 ð0:997; 1Þ, which means quasi-static thermal load.

r
r>a

ð20aÞ ð20bÞ 



To satisfy Eq. (20b), the factor Að1Þ ðsÞ þ Bð1Þ ðsÞ is in the form

0

fCðsÞg

ð19Þ

3.2. Anti-symmetric thermal shock

ð17Þ

fCðsÞg ¼ Að1Þ ðsÞ Bð1Þ ðsÞ Að2Þ ðsÞ Bð2Þ ðsÞ , and the detail is given in Appendix A. Substituting the solution of the vector fCðsÞg into Eq. (10a), the temperature of the structure under symmetric thermal shock is determined. The thermal shock formulated in Eq (14a) employs common mathematical idealizations. To satisfy the required initial conditions in general case, Heaviside function in Eq (14a) may be substituted by function H1 ðsac tÞ

H1 ðsac tÞ ¼ 1  ð0:5s

1 s2ac sac 1    p ðp þ sac Þ3 ðp þ sac Þ2 p þ sac

The expression of the temperature field under the ideal thermal shock transforms to that under a general thermal shock by replac  ð1Þ ð1Þ ing the factor q0 =ðpkc Þ with q0 = H1 ðsac ; pÞkc .

a

vector

2 2 ac t

H1 ðsac ; pÞ ¼

Að1Þ ðsÞ þ Bð1Þ ðsÞ ¼

Z

a

UðnÞsinðsnÞdn

ð21Þ

0

where UðnÞ is a new unknown function associated with the crack. From Eq. (21) and Eqs. (A1)–(A3) in Appendix A, the unknowns of the temperature solution in Eqs. (10a) and (10b) are obtained

fCðsÞg ¼ where



vb ðsÞ





Z

a

UðnÞsinðsnÞdn

ð22Þ

0



vb ðsÞ is given in Appendix A.

 

By the matrix operation Að1Þ ðsÞ  Bð1Þ ðsÞ ¼ v0 fCðsÞg, where

v0 ¼ ½ 1 1 0 0 , condition (20a) yields the equation of UðnÞ

Z

1

0

Z bðsÞJ 0 ðsrÞ

a

sUðnÞsinðsnÞdnds ¼ 

0

lq0 ð1Þ

pkc

;

r
ð23Þ

  where bðsÞ ¼ v0 vb ðsÞ cð1Þ . It is found that lims!þ1 bðsÞ ¼ 1. Employ the skill bðsÞ ¼ 1 þ ðbðsÞ  1Þ, the double integral in Eq. (23) can be decomposed into two terms. In addition, we deal the Ra integral X 1 ¼ 0 sUðnÞsinðsnÞdnds by the technique of integration by part. Finally, Eq. (23) is transformed into the Abel equation

Z 0

r

U0 ðnÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn þ r 2  n2

Z

a

UðgÞ

0

Z 0

1

ðbðsÞ  1ÞsJ 0 ðsrÞsinðsgÞdsdg ¼ 

lq0 ð1Þ

pkc

ð24Þ R1

where the known result [32] of the integral X 2 ¼ 0 J 0 ðrsÞcosðsnÞds is applied and it is given in Appendix D. The solution of Abel equation is well known. Using the known result [32] of the integral Rn 1=2 J 0 ðrsÞdr (which is given in Appendix D), the X 3 ¼ 0 rðn2  r 2 Þ solution of Eq. (24) is simplified to the form of Fredholm equation

UðnÞ þ

2

p

Z 0

a

UðgÞ

Z 0

1

ðbðsÞ  1ÞsinðsnÞsinðsgÞdsdg ¼ 

2 lq0

p pkð1Þ c

n ð25Þ

Fig. 2. Images of general heating curve H1 ðsac tÞ for different values of

sac .

This is a second kind of Fredholm equation which can be solved by discretization method. Once Eq. (25) is solved, the vector fCðsÞg can be calculated by Eq. (22) and then the temperature of the structure under anti-symmetric thermal shock is obtained.

321

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

4. Thermal shock fracture Substituting the elasticity constitutive equations (6) into the stresses balance equations (7), it has the general solution of displacements

in o 8 R 1 h ðjÞ > uðjÞ ¼ l 0 N1 ðz; sÞ DðjÞ ðsÞ J 1 ðsrÞds > > > >   > R 1 ðjÞ > ðjÞ ðjÞ > < þlaðjÞ 0 f 1 ðsÞ 1s AðjÞ ðsÞesc z þ BðjÞ ðsÞesc z J 1 ðsrÞds in o R 1 h ðjÞ > > wðjÞ ¼ l 0 N2 ðz; sÞ DðjÞ ðsÞ J 0 ðsrÞds > > > >  > R 1 ðjÞ 1  ðiÞ ðjÞ ðjÞ > : þlaðjÞ f ðsÞ A ðsÞesc z  BðjÞ ðsÞesc z J ðsrÞds 0

2

ð26Þ

and the general solution of stresses

Z a Ka ðsÞJ0 ðsrÞ sua ðnÞsinðsnÞdnds 0 0 Z 1 ½Ca ðsÞfCðsÞgJ 0 ðsrÞds ¼ 0; þ að1Þ

Z

h in o 8 ðjÞ ðjÞ R 1 > rðjÞ s N3 ðz; sÞ DðjÞ ðsÞ J 0 ðsrÞds > zz ¼ 2G 0 > > >  > R 1 ðjÞ  > ðjÞ ðjÞ > < þ2GðjÞ aðjÞ 0 f 3 ðsÞ AðjÞ ðsÞec sz þ BðjÞ ðsÞec sz J 0 ðsrÞds h in o ðjÞ ðjÞ R 1 > > rðjÞ s N4 ðz; sÞ DðjÞ ðsÞ J 1 ðsrÞds > rz ¼ 2G 0 > > >  > R 1 ðjÞ  ðjÞ ðjÞ ðjÞ > : þ2GðjÞ aðjÞ f ðsÞ A ðsÞesc z  BðjÞ ðsÞesc z J ðsrÞds 0

1

ð27Þ

be determined by boundary conditions.

Z

a

sua ðnÞsinðsnÞdn ¼ cosðasÞua ðaÞ þ

0

Under symmetric thermal shock with respect to the middle plane, the displacement u is symmetric while w is antisymmetric. Thus the stress rzz is symmetric and the stress rrz is anti-symmetric. The conditions on the middle plane are

rzz ðr; 0; tÞ ¼ 0; r < a

ð28aÞ

wðr; 0; tÞ ¼ 0;

ð28bÞ ð29Þ

The continuity conditions are given as

¼

þ uðr; h1 ; tÞ



ð30Þ

r

¼r

ð31Þ

þ zz ðr; h1 ; tÞ

ð32Þ

rrz ðr; h1 ; tÞ ¼ rrz ðr; hþ1 ; tÞ rzz ðr; h1 þ h2 ; tÞ ¼ 0

ð34Þ

rrz ðr; h1 þ h2 ; tÞ ¼ 0

ð35Þ R1 0

F a ðsÞJ 0 ðsrÞds ¼ 0 is sat-

h in o  1 ð1Þ  ð1Þ F a ðsÞ ¼ N2 ð0; sÞ Dð1Þ ðsÞ þ að1Þ f 2 ðsÞ Að1Þ ðsÞ  Bð1Þ ðsÞ s

ð36Þ

Reasonably, the function F a ðsÞ is in the form

Z F a ðsÞ ¼ 0

a

ua ðnÞsinðsnÞdn

cosðsnÞu0a ðnÞdn

ð41Þ

R1 Recall the integral X 2 ¼ 0 J 0 ðrsÞcosðsnÞds listed in Appendix D. As a result, the Abel equation emerges

Z a Z 1 u0a ðnÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi dn þ u ð g Þ sðKa ðsÞ  Ka0 ÞJ 0 ðsrÞsinðsgÞdsdg a 0 0 0 r 2  n2 Z 1 ¼ að1Þ ½Ca ðsÞfCðsÞgJ 0 ðsrÞds ð42Þ Z

Ka0

r

0

Its solution is simplified to

Ka0 ua ðnÞ þ

2

Z

p

Z

a

0

ua ðgÞ

1

ðKa ðsÞ  Ka0 ÞsinðsnÞsinðsgÞdsdg ¼ Q a ðnÞ

0

where

Q a ðnÞ ¼

ð33Þ

The surface boundary conditions are given as

According to Eq. (28b), the equation isfied for r > a, where

0

a

ð43Þ

þ

wðr; h1 ; tÞ ¼ wðr; h1 ; tÞ  zz ðr; h1 ; tÞ

ð40Þ

It is found that lims!þ1 Ka ðsÞ ¼ Ka0 , where Ka0 is a constant. We decompose the double integral in Eq. (39) into two terms through Ka ðsÞ ¼ Ka0 þ ðKa ðsÞ  Ka0 Þ. Applying the technique of integration by parts, we have

Z

 uðr; h1 ; tÞ

ð39Þ

h i Ka ðsÞ ¼ N3ð1Þ ð0; sÞ fg^a0 ðsÞg; h i ð1Þ ð1Þ ½Ca ðsÞ ¼ N3 ð0; sÞ ½g^amn ðsÞ þ f 3 ðsÞ½ 1 1 0 0 

4.1. Symmetric thermal shock

rrz ðr; 0; tÞ ¼ 0

r
where

n o where DðjÞ ðsÞ are four dimensional column vectors and unknown h i ðjÞ ðjÞ functions, the row vectors Nk ðz; sÞ and functions f k ðsÞ n o (k = 1,2,3,4) are listed in Appendix B. The unknown DðjÞ ðsÞ will

r>a

ð38Þ

0

1

4

n o 1 Dð1Þ ðsÞ ¼ fg^a0 ðsÞgF a ðsÞ þ að1Þ ½g^amn ðsÞfCðsÞg s

where the detail is given in Appendix C. Thanks to Eqs. (37) and (38), condition (28a) results the equation of ua ðnÞ

0

s

where the function ua ðnÞ is unknown. Substituting the general solutions of displacements and stresses into conditions (29)–(35), it n o gives the other seven equations of DðjÞ ðsÞ . Thus the eight n o unknown functions DðjÞ ðsÞ are related to the function F a ðsÞ. The n o relations between Dð1Þ ðsÞ and F a ðsÞ are denoted as

ð37Þ

2 lq0 að1Þ

p pkð1Þ c

Pa ðn; fÞ ¼ and

Z 0

where

Rn

2

1

Z

a

0

fPa ðn; fÞdf

ð44Þ

1 s

^ a ðsÞ J ðsfÞsinðsnÞds C 0 



^ a ðsÞ ¼ ½Ca ðsÞ v ðsÞ . C a

ð45Þ Note

that

the

integral

2 1=2

J0 ðrsÞdr (which has been given in Appendix D) X 3 ¼ 0 rðn  r Þ has been employed in above procedure. Eq. (41) is a second kind of Fredholm equation. It is easy to compute the numerical results of ua ðnÞ by discretization method. In order to improve the convergence rate of numerical solution of Eq. (43), the function Pa ðn; fÞ need pre-computation. It is found ^ a ðsÞ ¼ C ^ a0 , where C ^ a0 is a constant. We decompose that lims!þ1 C the function Pa ðn; fÞ into two integral through   ^ a ðsÞ ¼ C ^ a0 þ C ^ a ðsÞ  C ^ a0 . Substituting the result [32] of the inteC

322

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

R1 gral X 4 ¼ 0 s1 J 0 ðsfÞsinðsnÞds (it is copied in Appendix D), the function Pa ðn; fÞ is rewritten as

 Hðn  fÞ þ arcsinðn=fÞHðf  nÞ 2 Z 1  ^ a ðsÞ  C ^ a0 s1 J ðsfÞsinðsnÞds C þ 0

^ a0 Pa ðn; fÞ ¼ C

p

ð46Þ

0

ð1Þ

2G Ka0 ua ðaÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ OðrÞ; r 2  a2

r ! aþ

ð47Þ

Define the model I stress intensity factor (SIF) according to pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ðr  aÞrzz ðr; 0Þ. The SIF is evaluated as

K I ¼ limþ r!a

K I

Z

a

s3=2

The second term, an infinite integral, has a high convergence rate. Finally the solution of the elasticity field is completed. To investigate the fracture behavior of the crack, the stress on the crack plane is checked. Substituting Eqs. (37) and (38) into Eq. (27), and applying the integral [32] listed in Appendix D, the stress very close to the crack tip is given as

rzz ðr; 0Þ ¼ 

The function Kb ðsÞ approaches to a constant as s increasing, where the constant is denoted as Kb0 . We decomposed the double integral in Eq. (54) into two terms through Kb ðsÞ ¼ Kb0 þ ðKb ðsÞ  Kb0 Þ. Calculating the integral of ub ðnÞ by parts, it yields

pffiffiffi ¼ 2G lKa0 uðaÞ= a ð1Þ

Z

a

þ s1=2 0

n1=2 J 1=2 ðsnÞd½nub ðnÞ

ð56Þ

R1 Fortunately, the integration 0 J 1 ðsrÞJ 1=2 ðsnÞs1=2 ds is given in literature [32]. Thus Eq. (54) is transformed into the Abel equation

rffiffiffiffi Z 2

½nub ðnÞ0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dn p 0 r2  n2 Z a Z 1 pffiffiffi ðKb ðsÞ  Kb0 Þs3=2 J 3=2 ðsgÞJ 1 ðsrÞdsdg þr gub ðgÞ 0 0 Z 1 ¼ að1Þ r ½Cb ðsÞfCðsÞgJ 1 ðsrÞds

Kb0

r

ð57Þ

0

ð48Þ

4.2. Anti-symmetric thermal shock Under anti-symmetric thermal shock with respect to the middle plane, the displacement u is anti-symmetric while w is symmetric. Thus the stress rzz is anti-symmetric and the stress rrz is symmetric. The conditions on the middle plane are

rrz ðr; 0; tÞ ¼ 0; r < a

0

pffiffiffi nub ðnÞJ 3=2 ðsnÞdn ¼ ðsaÞ1=2 ub ðaÞJ 1=2 ðsaÞ

ð49aÞ

R1 1=2 Applying the integration 0 x2 ð1  x2 Þ J 1 ðyxÞdx given in literature [32], the solution of the Abel equation (57) is simplified to Z a pffiffiffiffiffiffi Z 1 ngub ðgÞ Kb0 ub ðnÞ þ ðKb ðsÞ  Kb0 ÞsJ 3=2 ðnsÞJ3=2 ðsgÞdsdg ¼ Q b ðnÞ 0

0

ð58Þ

where

Q b ðnÞ ¼ að1Þ n1=2

Z

a

UðfÞPb ðn; fÞdf

ð59Þ

0

uðr; 0; tÞ ¼ 0;

r>a

rzz ðr; 0; tÞ ¼ 0

ð49bÞ ð50Þ

The continuity conditions and surface boundary conditions are the same as Eqs. (30)–(35). R1 From Eq. (49b), the equation 0 F b ðsÞJ 1 ðsrÞds ¼ 0 is satisfied for r > a, where

h in o  1 ð1Þ  ð1Þ F b ðsÞ ¼ N1 ð0; sÞ F ð1Þ ðsÞ þ að1Þ f 1 ðsÞ Að1Þ ðsÞ þ Bð1Þ ðsÞ s

ð51Þ

Consequently, the function F b ðsÞ can be expressed in the form

F b ðsÞ ¼

pffiffi s

Z 0

a

pffiffiffi nub ðnÞJ 3=2 ðsnÞdn

ð52Þ

Z

Pb ðn; fÞ ¼ 0

ð53Þ

0

ð54Þ where

h

^ b0 Pb ðn; fÞ ¼ C

rffiffiffiffi

p

n3=2 gHðn  gÞ 2 Z 1  ^ b ðsÞ  C ^ b0 s1=2 J ðnsÞsinðsfÞds þ C 3=2

ð61Þ

The integral in Eq. (61) converges fast. Under anti-symmetric thermal shock, model II cracking may occur in the core. On this account, the shear stress on the crack plane is investigated. With the substitution of Eqs. (52) and (53) and integral listed in Appendix D, the stress very close to the crack tip is simplified to

rffiffiffiffi ð1Þ 2 2G Kb0 aub ðaÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ OðrÞ; rrz ðr; 0Þ ¼  p r r2  a2

r ! aþ

ð62Þ

Define the model II stress intensity factor according to pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K II ¼ limþ 2ðr  aÞrrz ðr; 0Þ. It follows

i

Kb ðsÞ ¼ N4ð1Þ ð0; sÞ fg^b0 ðsÞg; h i ð1Þ ð1Þ ½Cb ðsÞ ¼ N4 ð0; sÞ ½g^bmn ðsÞ þ f 4 ðsÞ½ 1 1 0 0 

ð60Þ

0

where the detail is given in Appendix C. With the substitution of Eqs. (52) and (53), condition (49a) gives the equation of ub ðnÞ ! Z Z 1 pffiffi a pffiffiffi Kb ðsÞs s nub ðnÞJ 3=2 ðsnÞdn þ að1Þ ½Cb ðsÞfCðsÞg J 1 ðsrÞds ¼ 0 0

^ b ðsÞs1=2 J ðnsÞsinðsfÞds C 3=2

  ^ b ðsÞ ¼ ½Cb ðsÞ v ðsÞ . The above Eq. (58) is a and where C b Fredholm equation which can be easily solved by numerical method. To improve the convergence rate of the numerical solution, we conduct further computation of function Pb ðn; fÞ. It is found that ^ b0 , where C ^ b0 is a constant. Note that the integral ^ ðsÞ ¼ C lims!þ1 C R 1 1=2 b s sinðs g ÞJ ðsnÞds has been exactly calculated [32]. We 3=2 0 decompose the function Pb ðn; fÞ into two integral through   ^ b0 þ C ^ b ðsÞ  C ^ b0 . Then function Pb ðn; fÞ is rewritten as ^ b ðsÞ ¼ C C

where the function ub ðnÞ is unknown. Using the definition (51) and conditions (50) and (30)–(35) we have the equation

n o 1 Dð1Þ ðsÞ ¼ fg^b0 ðsÞgF b ðsÞ þ að1Þ ½g^bmn ðsÞfCðsÞg s

1

r!a

ð55Þ

K II

rffiffiffiffi ð1Þ 2 2G lKb0 ub ðaÞ pffiffiffi ¼ p a

ð63Þ

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

323

5. Numerical results and discussion The previous sections demonstrate the solution of thermoelasticity fields in Laplace domain. In order to obtain the corresponding solution in time domain, which is of more interest, inverse Laplace transformation is applied. However, the Laplace transforms of the thermoelasticity fields are so complicated that the analytical inversion is almost impossible. Fortunately, there are some ordinary numerical methods [33] for the inversion of Laplace transforms that be used in this work. Stehfest method [34] is efficient and reasonable accurate for a fairly wide range of Laplace transforms [35]. Therefore we complete inverse Laplace transformation by this method. The algorithm is described as:

FðtÞ ¼

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

ð64Þ Fig. 5. Effect of heating rate on the peak of the temperature.

C n ¼ ð1ÞnþN

minðn;NÞ X m¼int½ðnþ1Þ=2

mN ð2mÞ! ðN  mÞ!m!ðm  1Þ!ðn  mÞ!ð2m  nÞ! ð65Þ

Usually, it is adequate for N = 8.

Fig. 6. Temperature histories at investigated point for different values of w1 and w2 for symmetric thermal shock.

Fig. 3. Verification of stress intensity factor.

Fig. 4. Temperatures of the plate under general heating and ideal heating.

Fig. 7. Temperature histories at investigated point for different values of w1 and w2 for anti-symmetric thermal shock.

324

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

Table 1 The detail thermal properties for different values of w1 and w2 . Parameters

kc

kd

sq

tð1Þ

tð2Þ

w1

w2

Case 1 Case 2 Case 3 Homogeneous

2/3 1 1/2 1

4/3 1/2 1/2 1

2/3 1 4 1

100 100 100 100

50 300 200 100

2 1/2 2 1

3 1/2 1/2 1

Fig. 8. Model I SIFs histories at the beginning of symmetric thermal shock for different values of w1 and w2 .

Fig. 10. Model II SIFs histories at the beginning of anti-symmetric thermal shock for different values of w1 and w2 .

Fig. 9. Model I SIFs histories for different values of w1 and w2 for symmetric thermal shock.

Fig. 11. Model II SIFs histories for different values of w1 and w2 for anti-symmetric thermal shock.

To verify the accuracy of the solution, the results given in previous literature [36] are introduced to be the reference. Matsunaga and Noda [36] researched the transient thermoelastic problem for an infinite plate containing a penny-shaped crack at an arbitrary position based on classic Fourier heat conduction model. Let material 2# be the same as material 1#, and reducing the DPL model into classic Fourier model, the problem studied in this paper is just the special case that the plate containing a center crack which has been discussed by Matsunaga and Noda [36]. Mathematically, let pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kc ¼ 1, kd ¼ 1, k ¼ 1, G ¼ 1, a ¼ 1 and cð1Þ ¼ cð2Þ ¼ 1 þ p=s2 , then the results in this work reduce to that for an homogeneous plate containing a center crack based on classic Fourier model. Introduc-

ing the reference value K s ¼ 2aGð3k þ 2GÞðk þ GÞ1 q0 kc a3=2 to normalize the thermal stress intensity factor, the results for the case that under symmetric thermal shock are calculated, as shown as in Fig. 3. It can be seen that the results in this work coincide with that given by Matsunaga and Noda [36]. Above all, the solution given in this paper is accurate. The dimensionless phase lag of temperature gradient t ¼ sT =sq always enhance the thermal wave nature, which has been investigated in the range of t that from 5 to 200, in previous studies [5,18]. According to Tzou’s calculation [5], the values of t for metals such as copper, gold and lead are 152, 114 and 70, respectively. Besides, Xu et al. [37] obtained the value of t for the casting sand

1

325

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

Fig. 12. Model I SIFs histories for different values of a for symmetric thermal shock.

Fig. 15. Model II SIFs histories for different values of G for anti-symmetric thermal shock.

on the thermoelastic behavior. In consequence, the value of t(1) is assumed to be 100 for the following numerical illustrations. Fig. 4 shows the temperature of the plate under symmetric heating. Both the results of general heating and ideal heating are exhibited. The temperature converges as the heating rate factor sac increases. Fig. 5 shows the peaks of the temperatures Tmax at the center of interface in the range of sac 2 ½1; 100. For sac ¼ 100, the results under general heating (T max ¼ ð0:0235; 0:00348Þ) are very close to the results under ideal heating (T max ¼ ð0:0236; 0:00361Þ). This means that the ideal heating conditions can be used if sac P 100. For illustration, the ideal heating conditions are considered in the following discussion. The coupling effects of the heat conduction properties of Material 2# on the temperature field are studied first. For illustration, the temperature at point (0, h1) is investigated. Two thermal propqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi erty factors are introduced as w1 ¼ kd sq tð1Þ =tð2Þ =kc and Fig. 13. Model II SIFs histories for different values of a for anti-symmetric thermal shock.

w2 ¼ w1 =sq [38]. Fig. 6 depicts the temperature histories of investigated point in the cases of symmetric thermal shock. And the results under anti-symmetric thermal shock are shown in Fig. 7. The temperature histories for different values of w1 and w2 are plotted. The detail thermal properties are listed in Table 1. In Figs. 6 and 7, the dimensionless temperature is normalized by ð1Þ

T 0 ¼ lq0 =kc . The temperature climbs to the peak value and then keeps in a quasi static state. It is found that the transient peak value and quasi-static value of the temperature increase with the thermal property factors w1 and w2 , respectively. Another object of particular interest is the thermally induced stress intensity factor (SIF). The effects of the properties of Material 2# on the SIFs are studied. We introduce K I0 ¼ ^ a0 a3=2 q =kð1Þ and K II0 ¼ ð4=3pÞGð1Þ að1Þ C ^ b0 a3=2 q =kð1Þ to Gð1Þ að1Þ C 0

Fig. 14. Model I SIFs histories for different values of G for symmetric thermal shock.

pile, which is 107. This paper tries to focus on the coupling effects of the thermal property differences between the skin and the core

c

0

c

normalize the model I and Model II SIFs, respectively. Therefore the effects of the thermal load and properties of material 1# are eliminated. The reference SIFs K I0 and K II0 are deduced from the solutions demonstrated in Sections 3 and 4 with the substitution of static thermal load and h1 ! 1. In the following numerical examples, Poisson’s ratios are fixed at 0.3. Therefore, Lame constants satisfy k ¼ 1:5G. The coupling effects of the heat conduction properties of Material 2# on the SIFs are studied. Figs. 8 and 9 draw the model I SIF histories for different values of w1 and w2 . Figs. 10 and 11 show the model II SIF histories for the selected values of w1 and w2 . It shows that the SIF shakes vigorously in the beginning of the ther-

326

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

mal shock and then keeps a quasi-static state. In the stage of the SIF oscillation, it appears a peak may be higher than final static level. From Fig. 8, the transient peak value of the model I SIF increase with the thermal property factor w1 . From Fig. 10, however, the effect of the thermal properties of Material 2# on the transient peak value of the model II SIF is not obvious. For comparing, the SIFs based on classical Fourier model are also exhibited in Figs. 8–11. It can be seen that the SIF based on the DPL model is considerable greater than it based on classical Fourier model. As expected, the results based on the two models converge as loading time increasing. In addition, the curves in Figs. 9 and 11 show that the static level of the SIF depends on the ratio kc of the thermal conductivity coefficients. Fig. 12 depicts the model I SIFs for different ratios of thermal expansion coefficients a. Fig. 13 depicts the model II SIFs for the selected ratios of a. The result for a ¼ 1 is also given for reference. It is found that the transient peak values of the SIFs are increased by material 2# in the cases of a > 1. And the transient peak values of the SIFs are almost not affected by material 2# in the cases of a < 1. As shown as Fig. 12, the static level of the model I SIF decreases with the ratio a. However, Fig. 13 proves that the static level of the model II SIF is independent of the ratio a. Fig. 14 plots the model I SIFs for different ratios of shear moduli G. And Fig. 15 plots the model II SIFs for the selected ratios G. The result for G ¼ 1 is also given for reference. It is found that both the transient peak value and static level of the SIF increase with the ratio G. The effect of the ratio G on the model I SIF is more significant than it on the model II SIF. 6. Conclusions A penny-shaped crack at the middle of a symmetric sandwich structure is considered in this paper. The thermoelasticity fields of the structure subjected to a sudden temperature gradient are solved based on the dual-phase-lag heat conduction. The stress intensity factors (SIFs) are evaluated to investigate the fracture behavior of the crack. The effects of the thermal and mechanical property differences between the skins and core on the temperature field and SIF are studied. The temperature field and SIF shake vigorously at the beginning of thermal shock (about t < 2). Late, they approach the static level slowly. It is found that the transient peak value and the static value of the temperature increase with the thermal property factors qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi w1 ¼ kd sq t1 =t2 =kc and w2 ¼ w1 =sq . The transient peak value of

Acknowledgment This research was supported by the Research Innovation Fund of Shenzhen City of China (project nos. JCYJ20170413104256729, JCYJ20170811160538023), the Natural Science Foundation of Guangdong Province of China (project nos. 2016A030311006, 2016A030310367), and the National Natural Science Foundation of China (project nos. 11672084, 11602072, 11372086). Appendix A Substituting the general solution of the temperature into the surface boundary conditions and the continuity conditions yields ð2Þ sðh þh Þ 1 2

ð2Þ sðh

Að2Þ ðsÞec

ð1Þ sh

Að1Þ ðsÞec

1

 Bð2Þ ðsÞec

þ Bð1Þ ðsÞec

ð1Þ sh 1

1 þh2 Þ

¼0

ðA1Þ ð2Þ sh

¼ Að2Þ ðsÞec

1

ð2Þ sh

þ Bð2Þ ðsÞec

1

ðA2Þ

    ð1Þ ð1Þ ð2Þ ð2Þ Að1Þ ðsÞec sh1  Bð1Þ ðsÞec sh1 ¼ w Að2Þ ðsÞec sh1  Bð2Þ ðsÞec sh1 ðA3Þ where

w ¼ kc 1ðpÞcð2Þ =cð1Þ ;

1ðpÞ ¼

1 þ tð2Þ sq p þ sq p þ 1 1 þ tð1Þ p

0:5p2 þ p þ 1 0:5s

2 2 qp

ðA4Þ

Under symmetric thermal shock, the condition on the middle plane gives

Z a   lq0 Að1Þ ðsÞ  Bð1Þ ðsÞ cð1Þ ¼  ð1Þ nJ 0 ðsnÞdn pkc 0

ð16aÞ

where the detail has been given in Section 3.1. Eqs. (A1)–(A3) and (16) consist of a linear algebraic equations of AðjÞ ðsÞ and BðjÞ ðsÞ. It gives

fCðsÞg ¼ 

 lq0  v ðsÞ kp a

Z

a

nJ 0 ðsnÞdn

ð17aÞ

0

where

2 

cð1Þ

6 cð1Þ sh 1 6e



cð1Þ

0

ð1Þ sh 1

ec

ec

va ðsÞ ¼ 6 6 ecð1Þ sh1 ecð1Þ sh1 4 0

0

ð2Þ sh

31 8 9 >1> > > 7 > ð2Þ sh <0> = c 1 7 e 7 ð2Þ sh 7 c 1 > > we 5 > >0> > : ; ð2Þ 0 ec sðh1 þh2 Þ 0

1

cð2Þ sh1

we

ð2Þ sðh þh Þ 1 2

ec

ðA5Þ

the model I SIF increases with w1 . But the effect of w1 on the transient peak value of the model II SIF is not considerably. Both the transient peak values of the model I SIF and model II SIF increase with the ratio of thermal expansion coefficients a when a > 1, but they are not affected when a < 1. In addition, both the model

Under anti-symmetric thermal shock, the condition on the middle plane results

I SIF and model II SIF increase with the ratio of shear moduli G. This paper only considered a sandwich structure consisting of homogeneous materials. In engineering practice, the commodiously designed fiber-matrix composite material and functionally graded material is prevalent and have been a major research interest among mechanics and materials communities [39–44]. Future research is clearly needed in order to create the knowledge of the combing effects of material nonlinearity and dual-phase-lag heat conduction.

In this case the unknowns AðjÞ ðsÞ and BðjÞ ðsÞ should be deduced from Eqs. (A1)–(A3) and (21).

Conflict of interest The authors declared that there is no conflict of interest.

Að1Þ ðsÞ þ Bð1Þ ðsÞ ¼

Z

a

UðnÞsinðsnÞdn

ð21aÞ

0

fCðsÞg ¼



vb ðsÞ



Z

a

UðnÞsinðsnÞdn

ð22aÞ

0

where

2 



1

1

cð1Þ sh

cð1Þ sh1

6e

vb ðsÞ ¼ 6 6

1

cð1Þ sh1

4e

0

e

ð1Þ ec sh1

0

0 cð2Þ sh

e

1

cð2Þ sh1

we

ð2Þ sðh þh Þ 1 2

ec

31 8 9 1 > > > > <0> = 7 > 1 e 7 7 ð2Þ c sh1 >0> > 5 > we > : > ; ð2Þ 0 ec sðh1 þh2 Þ 0

cð2Þ sh

ðA6Þ

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

h

Appendix B

h

i



h

i

h

N1ðjÞ ðz; sÞ ¼ esz zesz esz zesz N2ðjÞ ðz;sÞ ¼ esz

h



i     ð2dðjÞ þ 1Þ=s  z esz esz ð2dðjÞ þ 1Þ=s þ z esz

h

i

    h i N3ðjÞ ðz;sÞ ¼  esz z  ðdðjÞ þ 1Þ=s esz esz z þ ðdðjÞ þ 1Þ=s esz

h

i h N4ðjÞ ðz; sÞ ¼ esz

  z  dðjÞ =s esz

  i  dðjÞ =s þ z esz

esz

ðjÞ

h

h ðjÞ

o

in

Nð2Þ 4 ðh1 þ h2 ; sÞ

Dð2Þ ðsÞ

ðC5Þ

o ðC6Þ

Nð1Þ 4 ð0; sÞ

in

o  1 ð1Þ  Dð1Þ ðsÞ þ að1Þ f 4 ðsÞ Að1Þ ðsÞ  Bð1Þ ðsÞ ¼ 0 s

ðC7Þ

solution of

ðjÞ

8n o9 > < Dð1Þ ðsÞ > = 1 n o ¼ fg a0 ðsÞgF a ðsÞ þ að1Þ ½g amn ðsÞfCðsÞg ð2Þ > > s : D ðsÞ ;

f 4 ðsÞ ¼ f 2 ðsÞ Appendix C

h

Dð2Þ ðsÞ

Together the definition (36) and Eqs. (C1)–(C7), we have the n o DðjÞ ðsÞ

ðjÞ

f 3 ðsÞ ¼ f 1 ðsÞ ðjÞ

in

Nð2Þ 3 ðh1 þ h2 ; sÞ

Under symmetric thermal shock, the condition (29) on the middle plane gives

f 2 ðsÞ ¼ cðjÞ f 1 ðsÞ ðjÞ

ðC4Þ

The surface boundary conditions (34) and (35) give

 1 ð2Þ  ð2Þ ð2Þ þ að2Þ f 4 ðsÞ Að2Þ ðsÞesc ðh1 þh2 Þ  Bð2Þ ðsÞesc ðh1 þh2 Þ ¼ 0 s

3kðjÞ þ 2GðjÞ   2 ðkðjÞ þ 2GðjÞ Þ 1  ðcðjÞ Þ

ðjÞ

h in o ð2Þ Nð1Þ Dð1Þ ðsÞ  G N4 ðh1 ; sÞ Dð2Þ ðsÞ 4 ðh1 ; sÞ  1 ð1Þ  ð1Þ ð1Þ ¼ að1Þ f 4 ðsÞ Að1Þ ðsÞesc h1  Bð1Þ ðsÞesc h1 s  1 ð2Þ  ð2Þ ð2Þ þ að2Þ G f 4 ðsÞ Að2Þ ðsÞesc h1  Bð2Þ ðsÞesc h1 s

327

o

 1 ð2Þ  ð2Þ ð2Þ þ að2Þ f 3 ðsÞ Að2Þ ðsÞec sðh1 þh2 Þ þ Bð2Þ ðsÞec sðh1 þh2 Þ ¼ 0 s

dðjÞ ¼ GðjÞ =ðkðjÞ þ GðjÞ Þ f 1 ðsÞ ¼

in

The continuity conditions (30)–(33) give

where

in o h in o ð2Þ N Dð1Þ ðsÞ  N1 ðh1 ; sÞ Dð2Þ ðsÞ  1 ð1Þ  ð1Þ ð1Þ ¼ að1Þ f 1 ðsÞ Að1Þ ðsÞesc h1 þ Bð1Þ ðsÞesc h1 s   ð2Þ ð2Þ ð2Þ 1 ð2Þ f 1 ðsÞ Að2Þ ðsÞesc h1 þ Bð2Þ ðsÞesc h1 þa s h in o h in o N2ð1Þ ðh1 ; sÞ Dð1Þ ðsÞ  N2ð2Þ ðh1 ; sÞ Dð2Þ ðsÞ  1 ð1Þ  ð1Þ ð1Þ ¼ að1Þ f 2 ðsÞ Að1Þ ðsÞesc h1  Bð1Þ ðsÞesc h1 s  1 ð2Þ  ð2Þ ð2Þ þ að2Þ f 2 ðsÞ Að2Þ ðsÞesc h1  Bð2Þ ðsÞesc h1 s h in o h in o N3ð1Þ ðh1 ; sÞ Dð1Þ ðsÞ  G Nð2Þ Dð2Þ ðsÞ 3 ðh1 ; sÞ  1 ð1Þ  ð1Þ ð1Þ ¼ að1Þ f 3 ðsÞ Að1Þ ðsÞesc h1 þ Bð1Þ ðsÞesc h1 s   ð2Þ ð2Þ ð2Þ 1 ð2Þ þ a G f 3 ðsÞ Að2Þ ðsÞesc h1 þ Bð2Þ ðsÞesc h1 s ð1Þ 1 ðh1 ; sÞ

2

ð1Þ

f 2 ðsÞ

6 ð1Þ 6 f 4 ðsÞ 6 6 ð1Þ 6 f ðsÞescð1Þ h1 6 1 6 ð1Þ scð1Þ h1 6 1 6 f 2 ðsÞe ½g amn ðsÞ ¼ ½Da ðsÞ 6 ð1Þ 6 f ðsÞescð1Þ h1 6 3 6 ð1Þ 6 f ðsÞescð1Þ h1 6 4 6 6 0 4 0

ðC1Þ

ðC2Þ

ðC3Þ

ð1Þ

f 2 ðsÞ ð1Þ

ð1Þ

ð1Þ h

ð1Þ

ð1Þ h

f 2 ðsÞesc

1

1

ð1Þ ð1Þ f 3 ðsÞesc h1 ð1Þ ð1Þ f 4 ðsÞesc h1

0 0

Nð1Þ 2 ð0; sÞ

1 T fg a0 ðsÞg ¼ ½Da ðsÞ ½ 1 0 0 0 0 0 0 0 

0

f 1 ðsÞesc

i

3 0 7 6h i 7 6 7 6 Nð1Þ ð0; sÞ 0 7 6 4 7 6 i h i 7 6h ð2Þ 6 Nð1Þ ðh ; sÞ  N1 ðh1 ; sÞ 7 7 6 1 1 7 6 i h i 7 6h 7 6 ð1Þ ð2Þ  N2 ðh1 ; sÞ 7 6 N2 ðh1 ; sÞ 7 6 ½Da ðsÞ ¼ 6 h i h i 7 7 6 ð1Þ ð2Þ G N3 ðh1 ; sÞ 7 6 N3 ðh1 ; sÞ 7 6 6h i h i 7 7 6 ð1Þ ð2Þ 6 N4 ðh1 ; sÞ G N4 ðh1 ; sÞ 7 7 6 6 h i7 7 6 ð2Þ 6 N3 ðh1 þ h2 ; sÞ 7 0 7 6 4 h i5 ð2Þ 0 N4 ðh1 þ h2 ; sÞ

0

f 4 ðsÞ

2h

af 1ð2Þ ðsÞescð2Þ h1 af 2ð2Þ ðsÞescð2Þ h1 scð2Þ h1 aGf ð2Þ 3 ðsÞe ð2Þ aGf 4 ðsÞescð2Þ h1 ð2Þ ð2Þ af 3 ðsÞec sðh1 þh2 Þ ð2Þ ð2Þ af 4 ðsÞesc ðh1 þh2 Þ

0

ðC8Þ

ðC9Þ

ðC10Þ

3

7 7 7 7 ð2Þ ð2Þ 7 sc h1 af 1 ðsÞe 7 7 ð2Þ scð2Þ h1 7 af 2 ðsÞe 7 7 ð2Þ scð2Þ h1 7 aGf 3 ðsÞe 7 7 ð2Þ ð2Þ aGf 4 ðsÞesc h1 7 7 7 ð2Þ ð2Þ af 3 ðsÞec sðh1 þh2 Þ 7 5 scð2Þ ðh1 þh2 Þ af ð2Þ ðsÞe 4 0

ðC11Þ

328

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329

n o Employing matrix operation, Dð1Þ ðsÞ can be calculated by

n o 1 Dð1Þ ðsÞ ¼ fg^a0 ðsÞgF a ðsÞ þ að1Þ ½g^amn ðsÞfCðsÞg s

ðC12Þ

where

fg^a0 ðsÞg ¼ ½ ½I4  ½04  fg a0 ðsÞg

ðC13Þ

fg^amn ðsÞg ¼ ½ ½I4  ½04  fg amn ðsÞg

ðC14Þ

and where ½I4  is the indentify matrix of size 4 by 4, and ½04  is a zero matrix of size 4 by 4. Under anti-symmetric thermal shock, the condition (50) on the middle plane gives

h

in o  1 ð1Þ  Nð1Þ Dð1Þ ðsÞ þ að1Þ f 3 ðsÞ Að1Þ ðsÞ þ Bð1Þ ðsÞ ¼ 0 3 ð0; sÞ s

ðC15Þ

n o In this case, the unknowns DðjÞ ðsÞ of the general solutions are

determined by definition (51) and Eqs. (C1)–(C6) and (C15). That is

8n o9 > < Dð1Þ ðsÞ > = 1 n o ¼ fg b0 ðsÞgF b ðsÞ þ að1Þ ½g bmn ðsÞfCðsÞg > s : Dð2Þ ðsÞ > ;

where

fg^b0 ðsÞg ¼ ½ ½I4  ½04  fg b0 ðsÞg

ðC21Þ

fg^bmn ðsÞg ¼ ½ ½I4  ½04  fg bmn ðsÞg

ðC22Þ

Appendix D The integrals [32] applied in this paper are copied in the following:

Z

1

0

Z

rJ0 ðrsÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dr n2  r 2

n

0

Z

Z

1

J 0 ðsgÞsinðsnÞ

0

where

2h

ð1Þ 1 ð0; sÞ

ðC17Þ

1 T fg b0 ðsÞg ¼ ½Db ðsÞ ½ 1 0 0 0 0 0 0 0 

ðC18Þ

2

ð1Þ

f 1 ðsÞ

6 ð1Þ 6 f 3 ðsÞ 6 6 ð1Þ 6 f ðsÞescð1Þ h1 6 1 6 ð1Þ scð1Þ h1 6 1 6 f 2 ðsÞe ½g bmn ðsÞ ¼ ½Db ðsÞ 6 ð1Þ 6 f ðsÞescð1Þ h1 6 3 6 ð1Þ 6 f ðsÞescð1Þ h1 6 4 6 6 0 4 0

Z

1

0

Z

1

; r
Z

1

n 0

xJ 0 ðsnxÞ sinðsnÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi ffi dx ¼ s 1  x2

rffiffiffiffi

p 2

s1=2 J 3=2 ðsÞ

ds ¼ s ¼

N 0 7 6h i 7 6 7 6 Nð1Þ ð0; sÞ 0 7 6 7 6h 3 i h i 7 6 ð1Þ ð2Þ 7 6 N ðh1 ; sÞ  N ðh ; sÞ 1 1 7 6 1 7 6h i h i 7 6 ð1Þ ð2Þ 6 N2 ðh1 ; sÞ  N2 ðh1 ; sÞ 7 7 6 i h i 7 ½Db ðsÞ ¼ 6 h ð2Þ 6 Nð1Þ ðh ; sÞ G N3 ðh1 ; sÞ 7 7 6 3 1 6h i h i 7 7 6 ð1Þ ð2Þ 6 N ðh1 ; sÞ G N4 ðh1 ; sÞ 7 7 6 4 6 h i7 7 6 ð2Þ 6 0 N3 ðh1 þ h2 ; sÞ 7 7 6 h i5 4 Nð2Þ ðh þ h ; sÞ 0 1 2 4

n
0

!

3

i

1 pffiffiffiffiffiffiffiffiffi ; r 2 n2

x¼r=n

x2 J 1 ðsxÞ pffiffiffiffiffiffiffiffiffiffiffiffiffi ffi dx ¼ 1  x2

1

0

ðC16Þ

(

J 0 ðrsÞcosðsnÞds ¼

J 1 ðsrÞJ 1=2 ðsnÞs

1=2

Z

1

J 0 ðxÞsinðxn=gÞ

0

p 2

dx x

Hðn  gÞ þ arcsinðn=gÞHðg  nÞ

ds ¼

8 <0 1=2

2 n2 Þ n1=2 : ðr21=2 r Cð1=2Þ

¼

qffiffiffi 2

pr

;r < n 1=2 pnffiffiffiffiffiffiffiffiffi

r 2 n2

;n 6 r

rffiffiffiffi s1=2 sinðsgÞJ 3=2 ðsnÞds ¼

0

p 2

n3=2 gHðn  gÞ

References

ð1Þ

f 1 ðsÞ

0

ð1Þ

f 3 ðsÞ ð1Þ

ð1Þ h

ð1Þ

ð1Þ h

1

1

ð1Þ ð1Þ f 3 ðsÞesc h1 ð1Þ ð1Þ f 4 ðsÞesc h1

0 0

ðC20Þ

3

7 7 7 7 ð2Þ ð2Þ af 1 ðsÞesc h1 7 7 7 ð2Þ ð2Þ af 2 ðsÞesc h1 7 7 7 aGf 3ð2Þ ðsÞescð2Þ h1 7 7 7 ð2Þ scð2Þ h1 7 aGf 4 ðsÞe 7 7 ð2Þ ð2Þ af 3 ðsÞec sðh1 þh2 Þ 7 5 scð2Þ ðh1 þh2 Þ af ð2Þ 4 ðsÞe 0

scð2Þ h1 af ð2Þ 1 ðsÞe scð2Þ h1 af ð2Þ 2 ðsÞe scð2Þ h1 aGf ð2Þ 3 ðsÞe ð2Þ aGf 4 ðsÞescð2Þ h1 ð2Þ ð2Þ af 3 ðsÞec sðh1 þh2 Þ ð2Þ ð2Þ af 4 ðsÞesc ðh1 þh2 Þ

n o Employing matrix operation, Dð1Þ ðsÞ can be calculated by

n o 1 Dð1Þ ðsÞ ¼ fg^b0 ðsÞgF b ðsÞ þ að1Þ ½g^bmn ðsÞfCðsÞg s

0

0

f 1 ðsÞesc f 2 ðsÞesc

[1] L. Landau, The theory of superfluidity of helium II, J. Phys. 5 (1941) 71. [2] T.Q. Qiu, C.L. Tien, Heat transfer mechanisms during short-pulse laser heating of metals, J. Heat Transf. 115 (1993) 835. [3] 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. [4] P. Vernotte, Les paradoxes de la theorie continue de l’equation de la chaleur, C. R. Acad. Sci. 246 (1958) 3154–3155. [5] D.Y. Tzou, The generalized lagging response in small-scale and high-rate heating, Int. J. Heat Mass Transf. 38 (1995) 3231–3240.

ðC19Þ

[6] H.T. Chen, H.J. Lin, Study of transient coupled thermoelastic problems with relaxation time, ASME J. Appl. Mech. 62 (1) (1995) 208–215. [7] B.L. Wang, J.E. Li, Thermal shock resistance of solids associated with hyperbolic heat conduction theory, Proc. R. Soc. A 469 (2013) 2153.

S.L. Guo et al. / International Journal of Heat and Mass Transfer 139 (2019) 317–329 [8] S.L. Guo, B.L. Wang, Thermal shock fracture of a cylinder with a penny-shaped crack based on hyperbolic heat conduction, Int. J. Heat Mass Transf. 91 (2015) 235–245. [9] W. Li, F. Song, J. Li, R. Abdelmoula, C. Jiang, Non-Fourier effect and inertia effect analysis of a strip with an induced crack under thermal shock loading, Eng. Fract. Mech. 162 (2016) 309–323. [10] H.M. Huang, F. Su, Y. Sun, Thermal shock of semi-infinite body with multipulsed intense laser radiation, Acta Mech. Solida Sin. 23 (2010) 175–180. [11] 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. [12] Z.T. Chen, K.Q. Hu, Thermoelastic analysis of a cracked substrate bonded to a coating using the hyperbolic heat conduction theory, J. Therm. Stress 37 (2014) 270–291. [13] B.L. Wang, J.E. Li, Hyperbolic heat conduction and associated transient thermal fracture for a piezoelectric material layer, Int. J. Solids Struct. 50 (2013) 1415– 1424. [14] R. Zhang, X.Q. Fang, Y. Pang, On the dissipative transient waves in a piezoelectric microplate under strong thermal shock, Wave Random Complex 23 (2013) 1–10. [15] S.I. Anisimov, B.L. Kapeliovich, T.L. Perel’man, Electron emission from metal surfaces exposed to ultra-short laser pulses, Sov. Phys. - JETP 39 (1974) 375– 377. [16] D.Y. Tzou, Experimental support for the lagging response in heat propagation, J. Thermophys. Heat Transf. 9 (1995) 686–693. [17] B.L. Wang, J.E. Li, C. Yang, Thermal shock fracture mechanics analysis of a semiinfinite medium based on the dual-phase-lag heat conduction model, Proc. R. Soc. A 471 (2015) 20140595. [18] S.L. Guo, B.L. Wang, H.S. Nan, Coupling effects of inertia and dual-phase-lag heat conduction on thermal shock fracture of a cracked piezoelectric layer, Eng. Fract. Mech. 179 (2017) 278–293. [19] L.M. Chen, J.W. Fu, L.F. Qian, On the non-Fourier thermal fracture of an edgecracked cylindrical bar, Theo. Appl. Fract. Mech. 80 (2015) 218–225. [20] L. Zhang, X.M. Zhang, S. Peng, Z.M. Yan, Y. Liang, B. Yan, Q.B. Li, Thermalmechanical coupling propagation and transient thermal fracture in multilayer coatings, Heat Transf. Res. 48 (2017) 935–954. [21] F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena, Chaos, Solitons & Fractals 7 (9) (1996) 1461–1477. [22] X.Y. Zhang, X.F. Li, Thermal shock fracture of a cracked thermoelastic plate based on time–fractional heat conduction, Eng. Fract. Mech. 171 (2017) 22–34. [23] X.Y. Zhang, Z.T. Chen, X.F. Li, Thermal shock fracture of an elastic half-space with a subsurface penny-shaped crack via fractional thermoelasticity, Acta Mech. 229 (12) (2018) 4875–4893. [24] L.H. Wu, Z.Q. Jiang, W.J. Feng, An analytical solution for static analysis of a simply supported moderately thick sandwich piezoelectric plate, Struct. Eng. Mech. 17 (2004) 641–654. [25] K.Q. Hu, Q.H. Qin, Y.L. Kang, Anti-plane shear crack in a magnetoelectro elastic layer sandwiched between dissimilar half spaces, Eng. Fract. Mech. 74 (2007) 1139–1147. [26] T.B. Cheng, W.G. Li, D.N. Fang, Thermal shock resistance of ultra-hightemperature ceramics under aerodynamic thermal environments, AIAA J. 51 (2013) 840–848.

329

[27] 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. Stress 37 (2014) 14–33. [28] T.B. Cheng, W.G. Li, W. Lu, Y.S. Shi, D.N. Fang, Thermal shock resistance of ultra-high-temperature ceramic thermal protection system, J. Spacecraft Rockets 51 (2014) 986–990. [29] J. Zhao, Y. Li, X. Ai, Analysis of transient thermal stress in sandwich plate with functionally graded coatings, Thin Solid Films 516 (2008) 7581–7587. [30] S.H. Ding, Y.T. Zhou, X. Li, Interface crack problem in layered orthotropic materials under thermo-mechanical loading, Int. J. Solids Struct. 51 (2014) 4221–4229. [31] D.Y. Tzou, Macro- to Microscale Heat Transfer: the Lagging Behavior, John Wiley & Sons, 2014. [32] I.S. Gradshteyn, I.M. Ryzhik, Tables of Integrals, Series and Products, Acadeimic Press, San Diego, CA, 1965. [33] G. Honig, U. Hirdes, A method for the numerical inversion of Laplace transforms, J. Comp. Appl. Math. 10 (1) (1984) 113–132. [34] H. Stehfest, Algorithm 368: numerical inversion of Laplace transforms, Commun. ACM 13 (1970) 47–49. [35] B. Davies, B. Martin, Numerical inversion of the Laplace transform: a survey and comparison of methods, J. Comput. Phys. 33 (1979) 1–32. [36] Y. Matsunaga, N. Noda, Transient thermoelastic problem for an infinite plate containing a penny-shaped crack at an arbitrary position, J. Therm. Stresses 19 (4) (1996) 379–394. [37] Y.S. Xu, Y.K. Guo, Z.Y. Guo, Experimental research on transient heat transfer in sand, Acta Mech. Sin. 12 (1) (1996) 39–46. [38] S.L. Guo, B.L. Wang, K.F. Wang, J.E. Li, Coupling effects of dual-phase-lag heat conduction and property difference on thermal shock fracture of coating/substrate structures, Int. J. Solids Struct. 152 (2018) 238–247. [39] H. Wang, Q.H. Qin, Meshless approach for thermo-mechanical analysis of functionally graded materials, Eng. Anal. Bound. Elem. 32 (2008) 704–712. [40] W.J. Feng, E. Pan, X. Wang, Stress analysis of a penny-shaped crack in a magneto-electro-thermo-elastic layer under uniform heat flow and shear loads, J. Therm. Stresses 31 (2008) 497–514. [41] W.J. Feng, P. Ma, R.K.L. Su, An electrically impermeable and magnetically permeable interface crack with a contact zone in magnetoelectroelastic bimaterials under a thermal flux and magnetoelectromechanical loads, Int. J. Solids Struct. 49 (2012) 3472–3483. [42] H.P. Wu, L. Li, G.Z. Chai, F. Song, T. Kitamura, Three-dimensional thermal weight function method for the interface crack problems in bimaterial structures under a transient thermal loading, J. Therm. Stresses 39 (2016) 371–385. [43] W.J. Li, H.M. Huang, Z.M. Zhang, X.L. Xu, Effects of gradient density on thermal protection performance of AVCOAT composites under varied heat flux, Polym. Composite 37 (2016) 1034–1041. [44] Y. Qiu, H.P. Wu, J. Wang, J. Lou, Z. Zhang, A. Liu, G. Chai, The enhanced piezoelectricity in compositionally graded ferroelectric thin films under electric field: a role of flexoelectric effect, J. Appl. Phys. 123 (2018).