Thermo-viscoelastic response of a cracked, functionally graded half-plane under a thermal shock

Accepted Manuscript Thermo-viscoelastic response of a cracked, functionally graded half-plane under a thermal shock Wenzhi Yang, Zengtao Chen PII: DOI: Reference:

S0013-7944(18)30976-7 https://doi.org/10.1016/j.engfracmech.2018.11.042 EFM 6255

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

11 September 2018 15 November 2018 21 November 2018

Please cite this article as: Yang, W., Chen, Z., Thermo-viscoelastic response of a cracked, functionally graded halfplane under a thermal shock, Engineering Fracture Mechanics (2018), doi: https://doi.org/10.1016/j.engfracmech. 2018.11.042

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Thermo-viscoelastic response of a cracked, functionally graded half-plane under a thermal shock Wenzhi Yang, Zengtao Chen Department of Mechanical Engineering, University of Alberta, Edmonton, AB, Canada

Highlights:   

A thermo-viscoelastic model is developed for the crack problem in an FGM SIFs of elastic/viscoelastic FGMs are compared to show the effect of viscosity Parametric studies are conducted to show the effect of gradient parameters

Abstract: In this paper, a thermo-viscoelastic model is developed for the crack problem in an infinite, functionally graded half plane under a thermal shock. The moduli are assumed to be separable forms of both spatial coordinates and time, and two types of relaxation function of viscoelasticity are considered in this paper. By employing the Fourier transform and Laplace transform, coupled with the singular integral equations, the governing partial differential equations under mixed thermo-mechanical boundary conditions are solved numerically. The results show that the variations of stress intensity factors (SIFs) in viscoelastic functionally graded materials (FGMs) are significantly different than those in the elastic ones, which cannot be neglected in designing FGMs. key words: Functionally graded materials, Viscoelastic, Thermal shock, Crack problem, Stress intensity factors Nomenclature

c eij E G1 G2 h k K I , K II p sij

half-length of crack deviator components of strain tensor Young’s modulus shear relaxation function bulk relaxation function distance between crack and free surface heat conductivity stress intensity factors Laplace transform variable deviator components of stress tensor

t

T U

  , , ,  ij     ij  mi (t )

time temperature Airy stress function thermal expansion coefficient material constants of FGM strain thermal diffusivity Poisson’s ratio Fourier transform variable stress thermal expansion coefficient viscoelastic relaxation functions

1. Introduction Over the last decades, with the aid of developed manufacturing techniques like 3D printing, various soft materials such as hydrogels, elastomers, polymers and other composites have been invented and put into use in a broad range of applications such as biomedical, aerospace, tissue, automotive, mechanical, drug delivery, civil, and nuclear engineering [1-5]. Among them, polymer-based functionally graded materials (FGMs) are currently receiving much attention [69]. As a kind of multi-constituent material, to avoid severe thermal stresses and improve the structural strength, FGMs always possess gradual change in composition and microstructure so as to achieve continuous, spatial variations in physical and mechanical properties. For polymer-based FGMs, materials may exhibit significant creep and stress relaxation even under room temperature. Unlike elastic materials, deformation in viscoelastic FGMs is a timedependent process, and both stress and strain depend on their time variation. Moreover, in the framework of classical, linear viscoelasticity, materials will show more significant creep and stress relaxation under elevated temperature, and therefore, viscoelasticity cannot be neglected in designing polymer-based FGMs. FGMs can help reduce mechanical and thermal stresses to avoid catastrophic failure of structures, especially for the stresses around the cracks. There has been extensive research about

the crack problems in FGMs under mechanical or thermal loading [10-17]. Noda and Jin [18] pointed out that the stress intensity factors are still applicable in fracture problems of FGMs and their singularity and angular distribution in crack tip fields are the same with the homogenous materials. Following these works, a number of analytical or semi-analytical models are introduced to solve the crack problems, such as the piecewise linear or exponential models [1920]. However, compared with the elastic problems, there are only a few investigations focusing on the crack problem in viscoelastic FGMs. By solving the governing equations directly, Schovanec et al. [20-25] studied the stationary and dynamic mode I and mode III crack problems in nonhomogeneous, viscoelastic materials. Paulino and Jin [26-29] extended the elasticviscoelastic correspondence principle to FGMs under the assumption that the relaxation functions can be written in separable forms of space and time, which greatly simplified the analytical solutions. Using this correspondence, they analyzed a crack problem in viscoelastic FGMs strip under in-plane loading and anti-plane loading; Wang et al. [30] solved the viscoelastic

crack

problem

using

their

piecewise-exponential

model.

However,

the

correspondence principle can only be used in simple problems, in which one can turn the final elastic solutions to viscoelastic ones directly by corresponding substitutions. When the problem become complex, such as the thermal-viscoelastic problems, or when dimensionless variables exist, it is difficult to get the correct solutions by this correspondence. In this article, we aim to establish a thermo-viscoelastic analytical model for a cracked halfplane under a thermal shock in FGMs. By employing Laplace transform, the convolution in viscoelastic constitutive equation is avoided to simplify the derivation process. The partial differential equations are solved by Fourier transform, and the boundary conditions on crack faces are reduced to a set of singular integral equations. The Lobatto-Chebyshev method [31] is employed to solve the singular integral equations numerically, and the Laplace inversion [32] is used to get the results in the time domain.

2. Statement of the problem and basic equations As illustrated in Fig. 1, a semi-infinite, nonhomogeneous, viscoelastic plane containing a crack of length 2c under plane stress condition is considered. The crack is parallel to the free boundary, and is assumed to remain completely, thermally insulated, so the temperature field will be disturbed by the presence of this crack. The whole temperature is set to be zero initially,

and a thermal shock is suddenly applied on the free boundary at time t  0 . Inertia and body forces are neglected.

Fig. 1. Crack geometry and coordinates

In viscoelastic, nonhomogeneous materials, the moduli are functions of both spatial coordinates and time. To simplify the theoretical derivation, separable forms of the properties [12, 26-29, 33] are assumed as:

E  E0 exp(  y )m1 (t )

  0 (1   y ) exp(  y )m2 (t )    0 exp( y )m3 (t ) k  k0 exp( y )   0

(1)

where  ,  ,  ,  are material constants; E ,  ,  , k and  are the Young’s modulus, Poisson’s ratio, thermal expansion coefficient, heat conductivity and thermal diffusivity; and

mi (t ) (i  1, 2,3) are the relaxation functions in viscoelasticity. In the present work, two classical relaxation functions [34] are employed as: t

E E  m1 (t )  m2 (t )  m3 (t )  m(t )  (   (1   )e t0 ) E0 E0 t m1 (t )  m2 (t )  m3 (t )  m(t )  ( 0 ) q , 0  q  1 t

(2)

2.1 Heat conduction

The heat conduction when the heat generation is negligible is: 1 1 T  2T  k T  k  t

(3)

By introducing the following dimensionless variables,

T  T / T0 , t  t / (c 2 /  ), ( x, y, h)  ( x, y, h) / c,     c

(4)

Considering the material properties in equation (1), the governing equation for heat conduction can be reduced to:  2T  

T T  y t

(5)

Here and after, the hats of the dimensionless variables have been omitted for simplicity. The initial and boundary conditions for heat conduction in dimensionless forms are:

T  0, (t  0) t T ( x, h)  1, (t  0, x  ) T  0,

T  0, ( y  ) T  0, ( y  0, x  1) y

(6)

T ( x, 0 )  T ( x, 0 ), ( x  1) T ( x, 0 ) T ( x, 0 ) , ( x  1)  y y 2.2 Thermal-viscoelastic field equations The equilibrium equations for plane stress problems without considering the inertia effect are:  xy  y  x  xy   0,  0 x y x y

(7)

The strain-displacement relations are:

x 

u v 1 u v ,  y  ,  xy  (  ) x y 2 y x

(8)

The compatibility equation for two-dimensional problems is: 2  2 xy  2 x   y   2 y 2 x 2 xy

The constitutive law of viscoelasticity can be expressed as:

(9)

t

sij   G1 ( x, y, t   ) 0

 kk

deij d

d

d dT   G2 ( x, y, t   ) kk d  3  ( x, y, t   ) d   d d 0 0 t

t

(10)

with 1 1 sij   ij   kk  ij , eij   ij   kk  ij 3 3

where sij , eij are deviatoric components of the stress and strain tensors, T is temperature,

G1 ( x, y, t ) and G2 ( x, y, t ) are the shear and the bulk relaxation functions,  ( x, y, t ) is the thermal relaxation function. The relations between these functions and the relaxation function of Young’s modulus E * and Poisson’s ratio  * in the Laplace domain are expressed as [33]: G1* 

E* E* * , G  ,  *  pG2* * 2 * * 1  p 1  2 p

(11)

where p is the Laplace transform variable. Here and after, the superscript “*” denotes the variables in the Laplace domain. By using the Laplace transform, the above equations can be reduced to: 1 ( x*  p  * y* )  p *T * pE * 1 ( y*  p  * x* )  p *T *  * pE

 x*   y*

 xy* 

1  p * *  xy pE *

Combined with equation (1), the following equations are obtained:

(12)

u * 1  [ x*  p0 m2 ( p)(1   y ) exp(  y ) y* ]  p 0 m3 ( p) exp( y )T * x pE0 m1 ( p) exp(  y ) v* 1  [ y*  p0 m2 ( p)(1   y ) exp(  y ) x* ]  p 0 m3 ( p) exp( y )T * y pE0 m1 ( p) exp(  y )

(13)

u * v* 2(1  p0 (1   y ) exp(  y )m2 ( p)) *    xy y x pE0 exp(  y )m1 ( p) where mi ( p ) (i  1, 2,3) are the Laplace transform of mi (t ) (i  1, 2,3) . Let U * be the Airy stress function in the Laplace domain, in terms of which the stress can be expressed as:

 x* 

 2U *  2U *  2U * * *      , , xy y xy x 2 y 2

(14)

Substitute the above equations into equation (9) and (13), the governing equation is obtained:  2 2U *  2 

T *  2U *  2 * 2       2T * )  0 E  p m p m p   y T  ( 2U * )   2 ( ) ( ) exp(( ) )( 2 0 0 1 3 y y 2 y (15)

Besides (4), we introduce the following dimensionless variables,

 ij   ij / ( E0 0T0 ), U  U / ( E0 0T0 c 2 ) (u , v)  (u , v) / (c 0T0 ),  ij   ij / ( 0T0 )

(16)

( x, y , h )  ( x, y , h ) / c, (  ,  ,  )  (  ,  ,  )  c

The governing equation can be reduced to:  2 2U *  2 

  2U * T * 2 2 * ( 2U * )   2  p m ( p ) m ( p ) exp((    ) y )(  T  2    2T * )  0 1 3 2 y y y (17)

And the dimensionless constitutive equations are reduced to:

1 u * [exp(  y ) x*  p0 m2 ( p)(1   y ) y* ]  pm3 ( p) exp( y )T *  x pm1 ( p) 1 v* [exp(  y ) y*  p0 m2 ( p)(1   y ) x* ]  pm3 ( p) exp( y )T *  y pm1 ( p)

(18)

u * v* 2(exp(  y )  p0 (1   y )m2 ( p)) *  xy   pm1 ( p) y x Similar to Equation (5), the hats of the dimensionless variables have been omitted for simplicity. The mechanical boundary conditions can be expressed as:

 xy ( x, h)   y ( x, h)  0, ( x  )  xy ( x, 0)   y ( x, 0)  0, ( x  1)  xy ( x, 0 )   xy ( x, 0 ), ( x  1)

(19)

 y ( x, 0 )   y ( x, 0 ), ( x  1) u ( x, 0 )  u ( x, 0 ),

( x  1)

v( x, 0 )  v( x, 0 ),

( x  1)

3. Thermal stress We assume the thermal-viscoelastic coupling effect can be neglected, and the stress field will not influence the temperature field. By employing Laplace transform, the solution of (5) subjected to conditions (6) are given as [12]: 

T * ( x, y , p ) 

 D( , p) exp(

2

y  ix )d  

 

T ( x, y , p )  *

y0

2 D( , p) 1  exp(2 (h  y )) exp( 1 y  ix )d 1   2 exp( 2  h)







1 exp( ( y  h)); p

1 exp( ( y  h)); p

(20)

y0

where

1 

 2

  , 2 

 2

 ,  

p   2

2 4

, 

 2

 p

2 4

 is the Fourier transform variable and D( , p ) is associated with the density function of temperature, which can be calculated from the singular integral equation in [12]. In the Laplace

domain, substituting the solution of temperature distribution (20) into the following governing equation for thermal stresses:

 2 2U *  2 

T *  2U *  2 2 *        2T * ) (21) p m p m p   y T  ( 2U * )   2 ( ) ( ) exp(( ) )( 2 1 3 2 y y y

the general solution of (21) can be expressed as: U * ( x, y , p ) 









 ( B1  B2 y) exp(s2 y  ix )d 

 C exp[(     ) y  ix ]d , 1

2

y0 

U * ( x, y , p ) 

(22)

 {( A  A y)  ( A  A y) exp(2sy)}exp(s y  ix )d 1

2

3

4

1

 

  {C21  C22 exp(2  y ) exp[(     1 ) y  ix ]d  ,

y0



where A1 , A2, A3 , A4, B1 , B2 can be derived from the boundary conditions (19), and

s1  

 2

 s, s1  

 2

 s, s    2

2 4

and the following equations can be obtained by the particular solution of (21), C1 ( , p )  [(     2 )(  2 )   2 ]2 [ 2  p  (2   ) 2 ]D( , p ) p 2 m1 ( p )m3 ( p )

2 D( , p ) p 2 m1 ( p )m3 ( p ) 1  2 exp(2  a )  D( , p ) exp(2  a ) 2 p m1 ( p )m3 ( p ) C22 ( , p )  [(     2 )(  2 )   2 ]2 [(2   ) 2   2  p ] 2 1  2 exp(2  a ) C21 ( , p )  [(     1 )(  1 )   2 ]2 [ 2  p  (2   ) 1 ]

(23) From constitutive equation (18) and the corresponding boundary conditions, the jumps of the displacement on the line y  0 are:

[u * ] 1  [ x* ]  pm3 ( p)[T * ] x pm1 ( p)  2 [v* ] 1   [ exp(  y ) y*]  pm3 ( p) [T * ] 2 x pm1 ( p) y

(24)

The stress components can be obtained from solution (22) combined with equation (14). Substituting the stresses and temperature jump into (24), we have:

2 2  1 2 s2 ( B2  A4 )  s2 ( B1  A3 )  (2 s1 A2  s1 A1 )   [u * ]    exp(ix )d  2 2 ( ) x pm p ( ) ( ) ( )   C    C  C         1  1 21 2 22 1  



1  2 exp(ix )d  2  a  e  1 2

 pm ( p) D( )  3



(   s2 )[ s2 2 ( B1  A3 )  2 s2 ( B2  A4 )]  s2 2 ( B2  A4 )   [v ] i 1  2 2 2   ( ) (   s1 )(2 s1 A2  s1 A1 )  s1 A2  (     1 ) (  1 )C21  exp(ix )d  x  pm1 ( p)    2 (     2 ) (  2 )(C22  C1 )  

*





1  2 exp(ix )d  2  a 1  2e

i

 (  ) pm ( p) D( )  3



(25) Following the procedure in [31], introducing two dislocation density functions as:  1*  x, p  

[u * ( x, p )] , x

 2*  x , p  

[v* ( x, p )] x

(26)

By applying the second mechanical boundary condition on crack faces in equation (19), the following singular integral equations can be obtained as: 1

2

 ij

  [  x  k

1 j 1

ij

( x, )] j * ( , p )d  4 pm3 ( p )Wi* ( x, p ), i  1, 2,

1  x  1

(27)

with 1



* i

( x, p )dx  0,

i  1, 2

(28)

1

The kernels are given by: 

k11 ( x, )   [1  4 f11 ( )]sin[( x   ) ]d  0



k22 ( x, )   [1  4 2 f 22 ( )]sin[( x   ) ]d  0



k12 ( x, )   4 f12 ( ) cos[( x   ) ]d  0



k21 ( x, )   4 2 f 21 ( ) cos[( x   ) ]d  0

and

(29)



W ( x, p )  2   w1 ( , p ) sin( x )d  * 1

0



W2* ( x, p )  2   2 w2 ( , p ) cos( x )d 

(30)

0

h11 (  g1  2 g 2 )  2 sh12 ( s2 g1  g 2 )  g3 8s 3 h (  g1  2 g 2 )  2 sh22 ( s2 g1  g 2 ) w2* ( , p )   21  g4 8s 3 w1* ( , p )  

Where the expressions of fij , g1 , hij can be found in the Appendix. Using the Lobatto-Chebyshev method [31], the above singular integral equations can be transformed to algebraic equations: n 1  k11 ( xk , i )]F1* ( i , p )   Ai k12 ( xk , i )F2* ( i , p )  4 pm3 ( p )W1* ( xk , p )  i  xk i 1

n

 Ai [ i 1 n

 AF i 1

* i 1

( i , p )  0

n

n

1  k22 ( xk , i )]F2* ( i , p )  4 pm3 ( p )W2* ( xk , p ) Ai k21 ( xk , i )F ( i , p )   Ai [    x i 1 i 1 i k * 1

n

 AF i 1

i

2

*

(31)

( i , p )  0

where

 i ( , p)  *

Fi* ( , p ) 1 2

, (i  1, 2),   1

(i  1) , i  1, 2,....n; n 1 (2k  1) , k  1, 2,....n  1; xk  cos 2(n  1)

(32)

 i  cos

Ai 

 2(n  1)

, i  1, n; Ai 

 n 1

(33)

, i  2,3,.....n  1.

The stress intensity factors (SIFs) are defined as: K I * ( p)  

 4

F2* (1, p ), K II * ( p )  

 4

F1* (1, p )

(34)

To get the SIFs in the time domain, the inverse Laplace transformation has to be employed for above solutions. Since analytical inversion is difficult, the numerical inversion by Miller and Guy [32] will be used.

4. Numerical results and discussions The SIFs in the time domain can be obtained after the inverse Laplace transform. To begin with, the following standard, linear viscoelastic, relaxation function: t

m(t )  (

E E  E  (1   )e t0 ) , with   0.5, t0  1 E0 E0 E0

(35)

is considered [34], and the results of SIFs are shown in comparison with the corresponding elastic solutions. The elastic solution for the same thermal shock and same boundary conditions can be obtained when we assume the relaxation function to be m(t )  1 . In the following figures, the SIFs for comparison of elastic and viscoelastic materials and a series of parametric investigations are plotted versus dimensionless time. In all the cases, the distance between the boundary and crack faces, h , is assumed to be 1. It is worth noting the nonhomogeneous parameter  has no influence on the SIFs since it has been eliminated in the governing equations. Besides, from the singular integral equations (27), it is easy to find that only the relaxation function m3 (t ) in the thermal expansion coefficient will affect the results of SIFs. From a physical viewpoint, this can be explained by the space independence of relaxation functions. All three relaxation functions are assumed to be independent of space coordinates. However,  T (change of temperature) is a variable depending on the space coordinates in the present problem, which leads to the fact that only the relaxation function of thermal expansion, m3 (t ) , will enhance the nonhomogeneity of thermal stresses, thus resulting the variation of SIFs.

Fig. 2. Comparison of SIFs in viscoelastic and elastic materials, when (  ,  ,  )  (1, 1, 0.1) . Gradient parameters were selected to avoid possible contact of crack faces caused by negative values of K I .

As shown in Fig. 2, K I is positive while K II is negative in the select values of the nonhomogeneous material parameter [12, 35]. Compared to K I , the magnitudes of K II are more than five times of K I in both elastic and viscoelastic materials, which implies Mode II fracture plays a predominant role under thermal shock in both elastic and viscoelastic FGMs. Significant discrepancies are shown between elastic and viscoelastic results: the elastic SIFs increase from 0 at t=0 to a steady state around t=1.5, while the viscoelastic SIFs increase from 0 at t=0 to a peak value around t=0.3~0.5, after which the SIFs start to decrease dramatically, even approaching 0 at t=3.0, which indicates that although there is no much difference in the maximum of SIFs, fracture may occur at a much earlier stage in viscoelastic FGMs under thermal shock if the stress intensity factors exceed the fracture toughness.

Fig. 3. Variation of the cleavage stress versus angle  at time t  0.5, t  2.0 in elastic and viscoelastic nonhomogeneous materials when (  ,  ,  )  (1, 1, 0.1) .

Fig. 4. Variation of the maximum of cleavage stress versus time in elastic and viscoelastic nonhomogeneous materials when (  ,  ,  )  (1, 1, 0.1) .

To give a clearer comparison and illustrate the dynamic stress field around the crack tip, the variation of cleavage stress,   :

 

1 2 r

3  1 3 3  3 3    K I [ cos( )  cos( )]  K II [ sin( )  sin( )] 4 2 4 2 4 2 4 2  

(36)

versus angle  are shown in Fig. 3 when the nonhomogeneous material parameters (  ,  ,  )  (1, 1, 0.1) . Two different times t  0.5, t  2.0 are considered for both elastic and

viscoelastic materials. All curves show the cleavage stresses reach their maximum at the same angle, no matter whether the materials are elastic or viscoelastic, which indicates the crack propagation direction will keep the same at different times. Compared with elastic materials, the cleavage stresses in viscoelastic materials in later stage will be significant decreased. In fracture mechanics, the maximum of cleavage stresses can characterize the crack initiation under mixed mode facture, and their variations with time are depicted in Fig. 4. The results show that the main difference between elastic and viscoelastic materials are the maximum of cleavage stress reach its peak value at a much earlier time for the viscoelastic materials.

Fig. 5. The effect of the gradient parameter  on the SIFs when   1,   1 .

Fig. 6. The effect of the gradient parameter  on the SIFs when   1,   0.1 .

Fig. 7. The effect of the gradient parameter  on the SIFs when   1,   0.1 .

In Figs. 5-7, parametric investigations are conducted to illustrate the influence of nonhomogeneous parameters on the SIFs. Fig. 5 shows the effect of the gradient parameter  on

the SIFs when   1,   1 . It’s obvious that the variation of  , the gradient parameter of thermal expansion coefficient, has almost no influence on either K I or K II . The effect of the gradient parameter  on the SIFs is illustrated in Fig. 6. Compared with  ,  has a great influence on both K I and K II . When  is decreased from 1 to -1, the peak value of K I increases fourfold, while the peak value of K II is nearly doubled, which greatly enhance the fracture risk. In Fig. 7, the SIFs K I and K II are plotted versus dimensionless time considering the variation of the gradient parameter  when   1,   0.1 . Obviously,  has very little influence on K II but a great impact on K I . When   0.5 and   1 , K I takes a positive value, which means the crack faces are under tension. When   0.5 and   1 , K I becomes negative, which indicates the contact of crack faces would occur. As mentioned above, the Mode II fracture would play a predominant role in the present problem. It seems the variation of K I will have limited effect on the fracture behaviors. However, the contact of crack faces will increase the friction and make the heat transfer easier [12]. Besides, the crack contact also influences Mode II fracture as reported by reference [35-37]. Therefore, a detailed discussion considering crack closure problem in nonhomogeneous, viscoelastic material will be analyzed in a future work. Another power law relaxation function is also considered as t 1 m(t )  ( 0 ) q , q  , t0  1. t 2

(37)

Variations of SIFs with dimensionless time obtained from the power law relaxation function are compared with those of the standard linear function under the same nonhomogeneous parameters, as shown in Fig. 8. For K I , the magnitude of peak values for the power law function is a little higher than that for the standard linear one, and the rate of decay of K I is slower. Compared with K I , a significant difference is observed for K II , where the magnitude of peak values for the power law function is much higher than that for the standard linear one, and the decay rate is much slower. At time t  3.0 , the value of K II for the standard linear, viscoelastic material has approached zero, while it is still much higher than zero for the power law, viscoelastic material. All these results indicate that the actual type of relaxation function is critical for designing the viscoelastic FGMs.

Fig. 8. Comparison of SIFs for two different viscoelastic relaxation functions when   1,   1,   0.1 .

5. Conclusions In this paper, the transient SIFs and the thermal stresses under a thermal shock in FGMs are investigated when the crack is parallel to the free boundary of infinite half plane. Two classical types of relaxation function of linear viscoelasticity are considered and the moduli of FGMs are assumed to be forms of both spatial coordinates and time. By employing the Fourier transform and Laplace transform, coupled with the singular integral equations, the governing partial differential equations under mixed thermo-mechanical boundary conditions are solved numerically. The results show that significant discrepancies of transient SIFs exist between nonhomogeneous elastic and viscoelastic materials. The fracture risk occurs at a much earlier stage in viscoelastic FGMs under thermal shock than the elastic ones. An increase of the gradient of thermal conductivity in viscoelastic FGMs can lower the thermal stresses and fracture risks. The actual type of relaxation function is critical for designing the viscoelastic FGMs. All the results show that the viscoelastic properties should be concerned in designing FGMs, especially the “soft” polymer-based FGMs.

Acknowledgement The authors thank the Natural Sciences and Engineering Research Council of Canada (NSERC) and China Scholarship Council for the financial support to the present work.

APPENDIX

h11 ( )   s1  exp[( s1  s2 )h][ s1  hs2 ( s2  s1 )] h12 ( )  1  exp[( s1  s2 )h][1  h( s2  s1 )(1  hs2 )] h21 ( )  1  exp[( s1  s2 )h][1  h( s2  s1 )] h22 ( )  h 2 ( s2  s1 ) exp[( s1  s2 )h] f11 ( )  [  h11  s2 ( s1  s2 )h12 ]( s1  s2 ) 3 f12 ( )  [2 h11   ( s1  s2 )h12 ]( s1  s2 ) 3 f 21 ( )  [  h21  s2 ( s1  s2 )h22 ]( s1  s2 ) 3 f 22 ( )  [2 h21   ( s1  s2 )h22 ]( s1  s2 ) 3 g1 ( )   s2 2 f3  2 s2 f 4  f5 g 2 ( )  (2 s2   ) s2 2 f3  (3s2  2  ) s2 f 4  f 6 g3 ( )  exp( s2 h)[(1  hs2 ) f 2  hs2 2 f1 ]  (     1 ) I 21  (     2 ) I 22 g 4 ( )  exp( s2 h)[(1  hs2 ) f1  hf 2 ]  I 21  I 22

I1 ( , p ) 

1 C1 ( , p ) p m1 ( p )m3 ( p ) 2

I 21 ( , p ) 

1 C21 ( , p ) p m1 ( p )m3 ( p )

I 22 ( , p ) 

1 C22 ( , p ) p m1 ( p )m3 ( p )

2

2

f1 ( )  I 21 exp[(     1 )h]  I 22 exp[(     2 )h] f 2 ( )  (     1 ) I 21 exp[(     1 )h]  (     2 ) I 22 exp[(     2 )h] f3 ( )  I1  I 21  I 22 f 4 ( )  (     1 ) I 21  (     2 )( I 22  I1 ) f5 ( )  (     1 ) 2 I 21  (     2 ) 2 ( I 22  I1 ) 

1  2 D( ) 1  2 e 2  h

f 6 ( )  (     1 ) 2 (  1 ) I 21  (     2 ) 2 (  2 )( I 22  I1 ) 

1  2 D( ) 1  2 e 2  h

Reference [1] Drury JL, Mooney DJ. Hydrogels for tissue engineering: scaffold design variables and applications. Biomaterials. 2003;24:4337-51. [2] Guo M, Pitet LM, Wyss HM, Vos M, Dankers PY, Meijer EW. Tough stimuli-responsive supramolecular hydrogels with hydrogen-bonding network junctions. Journal of the American Chemical Society. 2014;136:6969-77. [3] Luo F, Sun TL, Nakajima T, Kurokawa T, Zhao Y, Sato K, et al. Oppositely charged polyelectrolytes form tough, self-healing, and rebuildable hydrogels. Advanced materials. 2015;27:2722-7. [4] Sun JY, Zhao X, Illeperuma WR, Chaudhuri O, Oh KH, Mooney DJ, et al. Highly stretchable and tough hydrogels. Nature. 2012;489:133-6. [5] Haag SL, Bernards MT. Polyampholyte Hydrogels in Biomedical Applications. Gels. 2017;3:41. [6] Singh AK, Siddhartha. A Novel Technique for Manufacturing Polypropylene Based Functionally Graded Materials. International Polymer Processing. 2018;33:197-205. [7] Singh AK, Vashishtha S. Mechanical and Tribological Peculiarity of Nano-TiO2-Augmented, Polyester-Based Homogeneous Nanocomposites and Their Functionally Graded Materials. Advances in Polymer Technology. 2018;37:679-96. [8] Parameswaran V, Shukla A. Dynamic fracture of a functionally gradient material having discrete property variation. Journal of Materials Science. 1998;33:3303-11. [9] Lambros J, Santare M, Li H, Sapna G. A novel technique for the fabrication of laboratory scale model functionally graded materials. Experimental Mechanics. 1999;39:184-90. [10] Bao G, Wang L. Multiple cracking in functionally graded ceramic/metal coatings. International Journal of Solids and Structures. 1995;32:2853-71. [11] Choi HJ, Lee KY, Jin TE. Collinear cracks in a layered half-plane with a graded nonhomogeneous interfacial zone–Part I: Mechanical response. International Journal of Fracture. 1998;94:103-22. [12] ZhiHe J, Naotake N. Transient thermal stress intensity factors for a crack in a semi-infinite plate of a functionally gradient material. International Journal of Solids and Structures. 1994;31:203-18. [13] Erdogan F, Wu B. The surface crack problem for a plate with functionally graded properties. Journal of applied Mechanics. 1997;64:449-56. [14] Gu P, Asaro RJ. Cracks in functionally graded materials. International Journal of Solids and Structures. 1997;34:1-17. [15] Gu P, Asaro RJ. Crack deflection in functionally graded materials. International Journal of Solids and Structures. 1997;34:3085-98. [16] Jin ZH, Batra R. Some basic fracture mechanics concepts in functionally graded materials. Journal of the Mechanics and Physics of Solids. 1996;44:1221-35.

[17] Jin ZH, Batra RC. Stress Intensity Relaxation at the Tip Of an Edge Crack In a Functionally Graded Material Subjected To a Thermal Shock. Journal of Thermal Stresses. 1996;19:317-39. [18] Noda N, Jin ZH. Crack-tip singularity fields in nonhomogeneous body under thermal stress fields. JSME international journal Ser A, Mechanics and material engineering. 1995;38:364-9. [19] Guo LC, Noda N. Modeling method for a crack problem of functionally graded materials with arbitrary properties—piecewise-exponential model. International Journal of Solids and Structures. 2007;44:6768-90. [20] Wang B, Mai YW, Noda N. Fracture mechanics analysis model for functionally graded materials with arbitrarily distributed properties. International Journal of Fracture. 2002;116:161-77. [21] Alex R, Schovanec L. An anti-plane crack in a nonhomogeneous viscoelastic body. Engineering fracture mechanics. 1996;55:727-35. [22] Herrmann J, Schovanec L. Quasi-static mode III fracture in a nonhomogeneous viscoelastic body. Acta mechanica. 1990;85:235-49. [23] Herrmann J, Schovanec L. Dynamic steady-state mode III fracture in a nonhomogeneous viscoelastic body. Acta mechanica. 1994;106:41-54. [24] Schovanec L, Walton JR. The quasi-static propagation of a plane strain crack in a power-law inhomogeneous linearly viscoelastic body. Acta Mechanica. 1987;67:61-77. [25] Schovanec L, Walton J. The energy release rate for a quasi-static mode I crack in a nonhomogeneous linearly viscoelastic body. Engineering fracture mechanics. 1987;28:445-54. [26] Paulino GH, Jin ZH. Correspondence Principle in Viscoelastic Functionally Graded Materials. Journal of Applied Mechanics. 2001;68:129. [27] Paulino GH, Jin ZH. Viscoelastic Functionally Graded Materials Subjected to Antiplane Shear Fracture. Journal of Applied Mechanics. 2001;68:284. [28] Paulino GH, Jin ZH. A crack in a viscoelastic functionally graded material layer embedded between two dissimilar homogeneous viscoelastic layers–antiplane shear analysis. International journal of fracture. 2001;111:283-303. [29] Jin ZH, Paulino GH. A viscoelastic functionally graded strip containing a crack subjected to in-plane loading. Engineering Fracture Mechanics. 2002;69:1769-90. [30] Wang ZH, Zhang L, Guo LC. A viscoelastic fracture mechanics model for a functionally graded materials strip with general mechanical properties. European Journal of Mechanics - A/Solids. 2014;44:75-81. [31] Chen ZT, Hu KQ. Thermo-Elastic Analysis of a Cracked Half-Plane Under a Thermal Shock Impact Using the Hyperbolic Heat Conduction Theory. Journal of Thermal Stresses. 2012;35:342-62. [32] Miller MK, Guy J, WT. Numerical inversion of the Laplace transform by use of Jacobi polynomials. SIAM Journal on Numerical Analysis. 1966;3:624-35. [33] Cheng ZQ, Meguid SA, Zhong Z. Thermo-mechanical behavior of a viscoelastic FGMs coating containing an interface crack. International Journal of Fracture. 2010;164:15-29.

[34] Sladek J, Sladek V, Zhang C, Schanz M. Meshless local Petrov-Galerkin method for continuously nonhomogeneous linear viscoelastic solids. Computational Mechanics. 2005;37:279-89. [35] El-Borgi S, Erdogan F, Hidri L. A partially insulated embedded crack in an infinite functionally graded medium under thermo-mechanical loading. International Journal of Engineering Science. 2004;42:371-93. [36] Tamuzs V, Petrova V, Romalis N. Plane problem of macro-microcrack interaction taking account of crack closure. Engineering Fracture Mechanics 1996;55:957–67. [37] Petrova V, Schmauder S. Crack closure effects in thermal fracture of functionally graded/ homogeneous bimaterials with systems of cracks. ZAMM‐Journal of Applied Mathematics and Mechanics/Zeitschrift Für Angewandte Mathematik und Mechanik. 2015;95:1027-36.