Analytical calculation of the stress intensity factor in a surface cracked plate submitted to thermal fatigue loading

Engineering Fracture Mechanics 77 (2010) 2354–2369

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Analytical calculation of the stress intensity factor in a surface cracked plate submitted to thermal fatigue loading Hoai Nam Le, Catherine Gardin * Institut Pprime, CNRS-ENSMA UPR 3346, Universite de Poitiers, 86961 Futuroscope Chasseneuil Cedex, France Departement Physique et Mecanique des Materiaux, ENSMA, Téléport 2, 1, Avenue Clement Ader, BP 40109, 86961 Futuroscope Chasseneuil Cedex, France

a r t i c l e

i n f o

Article history: Received 5 November 2009 Received in revised form 26 March 2010 Accepted 10 May 2010 Available online 31 May 2010 Keywords: Stress intensity factor Thermo-elasticity Thermal cycling Crack growth Cooling/heating time

a b s t r a c t A stainless steel specimen with a pre-existing surface notch is exposed to a convective medium of cyclic temperature. The history of stress intensity factor of the cracked body for different crack lengths is obtained by a closed-form integration of the stress field, using Duhamel’s theory with principle of superposition and appropriate weight functions. The obtained results are compared with numerical simulations performed with ABAQUS and they appear to be in very good agreement. The stress intensity factor history shows that fatigue behavior does not depend only on temperature amplitude DT ¼ T max  T min , quenching rate, and duration of thermal shock but also on heating rate and duration. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction A large number of structural components in aerospace, nuclear and electronic engineerings are subjected to a wide variety of loadings of mechanical and thermal origins [1,2]. The phenomenon of fatigue may be induced by a cyclic thermal gradient from the surface to the core of the structures. The prediction of the fatigue durability of cracked components involves an analysis of fatigue crack growth and requires accurate stress intensity factor (SIF) solutions corresponding to the transient temperature distribution [3,4]. In 1998, a leak occurred in the reactor heat removal system of the Civaux plant. The cause of cracking was identified as high cycle thermal fatigue [5,6]. This incident has inspired several studies [7–12] on thermo-mechanical fatigue of stainless steel material used in various components of the pressurized water reactors (PWR). Experimental results obtained using a facility called SPLASH [13] show that crack initiation networks can be obtained by cyclic thermal loading. This facility allows submitting a parallelepipedic specimen to cyclic thermal down-shocks by thermal sprays on two opposite sides, implying large thermal gradients in the thickness (Fig. 1). The duration of the thermal cycle is 7.75 s, with only 0.25 s for cooling, with a temperature range around 200 °C. It appears that these crack networks are limited on the specimen surface and the deepest crack does not exceed 2.5 mm even with the most severe loading. Another experimental facility, called FAT3D, has been developed in CEA Saclay (France) [8], on tubes cooled cyclically in an angular portion, implying a 3D thermal fatigue loading. The cooling time is around 15 s, and the duration of the cycle ranges between 90 and 190 s. These different tests have demonstrated that a network of cracks can be initiated and can propagate under pure cyclic thermal loading. These networks, located in the border of the cooling zone, generally present a dominant crack which propagates quickly and may, according to the level of the loading, propagate through the thickness * Corresponding author. Tel.: +33 05 49 49 82 35; fax: +33 05 49 49 82 38. E-mail address: [email protected] (C. Gardin). 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.05.012

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

2355

Nomenclature a crack length Bi Biot number Cp specific heat D thermal diffusivity E Young’s modulus H convective heat transfer coefficient K stress intensity factor (SIF) Kmax maximum value of SIF during one cycle minimum value of SIF during one cycle Kmin Kth fatigue crack propagation threshold 2L specimen thickness Rt cooling–heating time ratio minimum/maximum SIF ratio RK Ti initial temperature of specimen Tmax/Tmin maximum/minimum temperature during cycle tC duration of temperature cycle duration of heating/cooling tup/tdo a coefficient of thermal expansion k thermal conductivity t Poisson ratio q density

of the test-tube. Thus it appears that the characteristics of the thermal cycling play an important role on crack propagation. This remark is also reported through the work of Lejeail et al. [10]. In practical applications, finite element method (FEM) is used to study the stress and strain fields as well as stress intensity factor induced by fatigue thermal shock [10,14–17]. The main disadvantage of this method is that complicated numerical calculation must be repeated for the same cracked body, subjected to different mechanical or thermal loadings. To simulate cracking, the computation of stress intensity factors for a cracked body subjected to a thermal transient loading by the finite element method requires a step-by-step computation in order to build the entire crack growth. This is expensive in computational resources and time-consuming because of the number of the parameters studied and of the highly transient nature of thermal shock problem. Another approach in design and analysis of thermal shock fatigue facilities is then welcomed. In the area of the thermo-elastic stress field, Lu and Fleck [18] solved the thermal shock problem for an orthotropic plate of finite thickness suddenly exposed to a convective medium of different temperature. Closed-form expression of the transient thermal stress is obtained over the full range of Biot number. It appears that the maximum tensile stress obtained at the specimen surface increases with increasing Biot number. It is important for thermal fatigue analysis to obtain the transient temperature, thermal stress as well as stress intensity factor induced by cyclic variation of temperature. Marie [19] suggests that any variation of fluid temperature can be divided in succession of linear shocks. The analytical solution for any variation of temperature can then be obtained as the analytical solution in case of linear shock is known. Then, by using a thermo–elastic stress strain behavior, transient thermal stress can be deduced. Finally, the stress intensity factor at the surface and deepest point of a semi-elliptical crack can be achieved using the shape function method [20,21].

Fig. 1. Example of experimental temperature evolution in SPLASH test [31].

2356

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

Other authors [22,23] propose using Green’s function method (GFM) and Duhamel’s theorem. Kim et al. [22] model a semi-elliptical surface crack subjected to a thermal transient loading, studying the influence of the aspect ratio of the crack. Lee and colleagues [23] study an axisymmetric crack, submitted to triangular and sinusoidal temperature variations, putting in light that the Green’s function determined at a fixed temperature can be applied to a relatively wide range of temperatures. Green’s function separates the influences of geometry and loading. Once Green’s function for the geometry is known, the SIFs for the other loading cases can be easily determined. The efficiency of GFM has been demonstrated because it is less time-consuming compared to FEM. In the present paper, Duhamel theorem [24] and principle of superposition are used to obtain the closed-from expression of the SIF of a stainless steel plate of finite thickness with faces exposed to a convective medium whose temperature is cyclically varying. Two types of thermal cycling are considered: sinusoidal variation of temperature and saw-toothed fluctuation of temperature which approaches the real one in main thermal fatigue test facilities. Parametric studies are carried out to bring to light the influence of thermal loading history on SIF. 2. Analytical approach In this part, the aim is to obtain an analytical expression of the stress intensity factor (SIF) of a cracked plate submitted to cyclic thermal loading at its edges, in the temperature range [Tmax – Tmin]. A preliminary study of thermal cold shock is carried out, as the further investigation of cyclic boundary conditions will use the Duhamel’s principle [24]. A classical sinusoidal fluctuation of temperature is then studied. Sinusoidal signals allow expressing any type of signal using Fourier’ transform, especially for cyclic experiments. Nevertheless, in considered experiments [8,25,26], where the specimen is cooled by water, it appears that the cooling rate is far much higher than the heating one (see Fig. 1). So, a saw-toothed temperature fluctuation will also be considered in this paper, in order to approach experimental conditions. As the calculation of the SIF is based here on weight functions, we have in a first step determined the cyclic temperature and stress histories in a crack-free specimen. 2.1. Cyclic temperature and stress histories of a crack-free specimen 2.1.1. Preliminary study of thermal shock A crack-free infinite plate of thickness 2L is considered (see Fig. 2) with its Cartesian coordinates embedded at the centre of the plate. The employed material has a density q, a thermal conductivity k, a specific heat Cp. Initially, this plate is at a uniform temperature Ti, and at time t = 0, its two edges (x = ±L) are suddenly exposed to a convective medium at temperature Tf. The heat transfer coefficient H is assumed to be constant. Heat flow within the solid induces a transient temperature distribution T(x, t) and a stress state rðx; tÞ. The expression of the transient temperature distribution in this classical conduction problem is the following [24]: 1 Tðx; tÞ  T f X ¼ Ti  Tf n¼0



    x 2 sinðkn Þ 2 Dt exp kn 2 cos kn kn þ sinðkn Þ cosðkn Þ L L

ð1Þ

where kn are the roots of kntan kn = Bi. The non-dimensional heat transfer coefficient Bi ¼ HL is the Biot number. k D ¼ C pkq is the thermal diffusivity.

λ, ρ, Cp Ti

H, T f

H, T f

y x

z -L

2L

L

Fig. 2. Crack-free specimen under thermal shock conditions.

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

2357

Then, we assume an uncoupled linear thermo-mechanical behavior:

eij ðx; tÞ ¼

1þm m rij ðx; tÞ  trrðx; tÞdij þ aðTðxÞ  T i Þdij E E

i ¼ x; y; z

ð2Þ

where E is the Young’s modulus, m the Poisson ratio, and a is the coefficient of thermal expansion. We will restrict this study to the determination of ryy stresses, responsible of mode I opening for a crack in the (x, z) plane. So, we only have to calculate the axial eii strains and rii stresses which are independent from shear behavior. Then, Eq. (2) leads to three equations with six unknown quantities exx ; eyy ; ezz ; rxx ; ryy ; rzz . Note that the following determination of the stresses and strains is done without crack in the specimen! The analytical solving requires additional assumptions which have been verified through an ABAQUS numerical simulation: – The normal stress is vanishing in through-thickness direction because the plate is free to expand in that direction and no external mechanical loading is applied:

rxx ¼ 0

ð3Þ

– The plate is free to expand with vanishing axial force in the y direction:

Z

L

ryy dx ¼ 0

ð4Þ

L

Nevertheless, due to edge effect near x = ±L, we cannot write that

ryy ¼ 0 8x.

– It is verified through numerical approach that the plate stretches uniformly along the y direction and does not bend, which implies that eyy is independent of all spatial dimensions including x:

eyy ðx; tÞ ¼ eyy ðtÞ

ð5Þ

Under bi-dimensional analysis, an additional plane strain assumption is introduced:

ezz ¼ 0 leading to

Z

L

rzz dx – 0

ð6Þ

L

Then, we deduce the analytical expression of the ryy(x, t) stress that will be later responsible for the mode I propagation of a crack in the xz plane, propagating in the x direction:

ryy ðx; tÞ ¼ EaðT i  T f Þ Wn ðxÞ ¼

1 1v

1 X

!

Wn ðxÞ expðcn tÞ

n¼0



2 sinðkn Þ kn þ sinðkn Þ cosðkn Þ

with

 2  x sinðkn Þ k D and cn ¼ 2n  cos kn kn L L

ð7Þ

2.1.2. 2.1.2.Time-dependent temperature conditions 2.1.2.1. The use of the Duhamel’s principle [24]. We wish in this part to calculate in a first step the temperature field versus time in a wall submitted to a time-dependent temperature at its two edges Te(t), as shown in Fig. 3. All the other parameters remain unchanged compared to paragraph 2.1.1. In order to solve analytically this problem, it is necessary to use the Duhamel’s principle [24]. The considered system with assumption of linear elastic behavior, named system I, has then to be considered as the superposition of system I–A and system I–B. In A system, the initial temperature of the specimen is set to zero, and the boundary temperature conditions at the edges are time-dependent. In B one, the initial temperature is Ti, but the two edges are placed in a medium at a 0 °C temperature. It is important to note that by applying the Duhamel’s principle and principle of superposition, analytical expression of temperature can be obtained everywhere in the plate. Because the approach to obtain the expression of temperature and the stress is very similar, our attention is focused only on the normal stress component ryy. The expression of the normal stress in the system I–B can be easily deduced from (7) by replacing Tf by 0 I—B yy ðx; tÞ

r

¼ EaT i

1 X

!

Wn ðxÞ expðcn tÞ

ð8Þ

n¼0

Duhamel’s principle can be applied to system I–A. The normal stress of system I–A with time-dependent condition boundary condition can be deduced with:  rIA yy ðx; tÞ ¼ T e ðtÞ  r ðx; 0ÞÞ þ

Z

0

t

T e ðx; t  sÞ 

@ r ðx; tÞ ds @t

ð9Þ

where r ðx; tÞ is the solution of normal stress for system I–A with a unit step temperature obtained by replacing Ti by 0 and Tf by 1 in (7):

2358

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

Te(t)

Te(t)

Tf = Te(t)

Tf = Te(t)

Ti = Ti

Ti = 0

Ti

Tf = 0

Tf = 0

2L

2L

2L

I

I-A

I-B

Fig. 3. Principle of superposition applied to the considered system.

r ðx; tÞ ¼ rIA yystep ¼ ð1ÞEa

1 X

!

Wn ðxÞ expðcn tÞ

ð10Þ

n¼0

Leading to, by superposition: I—B rIyy ðx; tÞ ¼ rI—A yy ðx; tÞ þ ryy ðx; tÞ

ð11Þ

These general expressions will now be applied to the two specific cases of sinusoidal and saw-toothed fluctuations of the temperature at the edges of the coupon. 2.1.2.2. Sinusoidal evolution of the boundary temperature. The boundary temperature conditions are the followings:

T e ðtÞ ¼ B þ A cosðxtÞ

ð12Þ

where B is the mean temperature, A the amplitude, and x the angular frequency. min min As the specimen edge temperature is inside the range [Tmax  Tmin], it can be deduced that A ¼ T max T ; B ¼ T max þT . 2 2 Using Eqs. (8)–(11) it follows, after simplification, that: I yy ðx; tÞ

r

¼ Ea

1 X

2

Wn ðxÞðT i  B  A sin un Þe

cn t



n¼0

1 X

!

Wn ðxÞA cos un ðcosðxt  un ÞÞ

ð13Þ

n¼0

where the phase-shift un is defined as:

c

n ; sin un ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2n þ x2

x

cos un ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2n þ x2

ð14Þ

In Eq. (13), the first term corresponds to a transient response, whereas the second term is the steady one. 2.1.2.3. Saw-toothed evolution of the boundary temperature. We will now consider the saw-toothed fluctuation of temperature of the convective medium. This fluctuation has different rates during cooling and heating (Fig. 4). The temperature of the conductive medium decreases linearly from Tmax to Tmin, tdo being the cooling duration. During the heating tup time, temperature then increases linearly up to Tmax. The corresponding period of the temperature cycle is tc = tup + tdo. The saw-toothed fluctuation of temperature Te(t) and the corresponding derivative T 0e ðtÞ with respect to time can then be analytically expressed as follows:

8 > < 0 6 t 6 tdo : T e ðtÞ ¼ It þ J tdo 6 t 6 t c : T e ðtÞ ¼ Mt þ N > : tc 6 t : T e ðtÞ ¼ Tðt  Nc tc Þ

and T 0e ðtÞ ¼ I and T 0e ðtÞ ¼ M and

T 0e ðtÞ

ð15Þ

0

¼ T ðt  Nc tc Þ

where Nc is the current number of cycles.

with I ¼

T min  T max ; t do

J ¼ T max ;

M¼

T max  T min ; tup

N ¼ J  M  tc

With the same approach than previously, the final solution of the normal stress can be obtained. It is given in Appendix A.

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

2359

temperature Te(t) Tmax

Tmin

0

tdo

tup

tc

time t

Fig. 4. Saw-toothed fluctuation of temperature.

crack y 0

x a

z Tf , H

Tf , H

T

2L Fig. 5. Plate with single edge crack under thermal boundary conditions.

2.2. Stress intensity factor of plate with single edge crack under thermal loading Let us now consider a plate with a single edge crack of depth a as show in Fig. 5. This crack is crossing the entire thickness in the z direction and it will propagate along x in the xz plane. We assume that the change in temperature, if any, due to the flaw, can be neglected. In the previous part, we obtained in (7), (13) and Appendix A the expressions of the ryy(x, t) stress of the crack-free body. The stress intensity factor of the cracked plate under thermal shock or time-dependent temperature conditions can be effectively evaluated by utilizing the weight function method. The characteristic of the weight function is independent of the loading, which enables the SIF to be calculated using a simple integral expression along the crack lips of a word-like product for weight function and loading:

K I ðtÞ ¼

Z

L

mðx; aÞ:ryy ðx; tÞ dx

ð16Þ

La

where m(x, a) is the weight function proposed by Tada et al. [27]. 3. Results and discussion 3.1. Problem data The material considered in this article, the stainless steel AISI 304L, is supposed to have an isotropic mechanical and thermal behavior. The effect of varying temperatures should be considered in performing the analysis as realistically as possible. However, it was shown that the effect of material properties is small due to compensation effect of thermal and mechanical properties over the temperature range [5]. The general properties used in this work are summarized in Table 1. The considered wall has a thickness 2L = 10 mm. The choice of the convective heat transfer coefficient H plays an important role. For forced convection by water and water vapor, recent calculations carried out on thermal fatigue tests [8,28] show that Hwater can vary between 6000 and 17,000 W/ m2 K1 corresponding respectively to Bi = 1.57 and Bi = 4.47. For conservative estimation under water cooling, Bi will be set to 1.57 in the following of this article. For cyclic temperature, the temperature Tmax is set at 350 °C and Tmin at 100 °C, representative of literature thermal fatigue test conditions [8,13,25].

2360

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

Table 1 Mechanical and thermal properties of AISI 304L used in calculation [13]. Density q (kg/m3)

Young’s modulus E (GPa)

Poisson’s ratio t

Threshold SIF Kth pffiffiffiffiffi (MPa m)

Specific heat Cp (J kg1 K1)

Thermal conductivity k (W m1 K1)

Coefficient of thermal expansion a (K1)

7800

196

0.3

5.2

500

19

18  106

3.2. Numerical approach ABAQUS is used to validate the analytical approach. The simulation is carried out in transient regime with a bi-dimensional coupled temperature–displacement CPE8T plain strain element type. The represented plate has here a finite dimension along y, but it was chosen high in order not to influence the vicinity of the crack (refer to Fig. 6). The mesh is refined around the crack tip so that the stress intensity factor can be correctly evaluated [29]. The plate is free to expand and contract. Another simplified model has also been studied, without crack, in order to obtain the temperature and stress evolutions in the crack-free body.

Fig. 6. Mechanical and thermal boundary conditions of numerical simulation.

450 Analysis-Heart Analysis-Edge FEM-Heart FEM-Edge

Temperature (°C)

400 350 300 250 200 150 100

0

25

50

75

100

time (s) Fig. 7. Temperature evolution on the edge (x = L) and in the heart (x = 0) of the specimen for Ti = 350 °C with Bi = 1.57.

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

2361

The mechanical boundary conditions indicated in Fig. 6 are as follows – displacement blocked along y on the axis of symmetry, – a point of this axis is embedded to avoid the overall translation. Moreover, an adiabatic thermal condition is imposed along the top surface of the specimen and on the line of symmetry including the crack surface, and the left and right sides of the specimen (x = ± L) are exposed to a convective medium of temperature T f ðtÞ with a coefficient of heat transfer H. 3.3. Thermal fatigue by sinusoidal evolution of the boundary temperature 3.3.1. Evolution versus time of temperature, stress, and stress intensity factor In this part, we consider the coupon placed in a convective medium (Bi = 1.57) with a sinusoidal variation of temperature. We assume that this temperature varies between Tmax = 350 °C and Tmin = 100 °C. In Figs. 7–9, the circular frequency is set at x = 0.2512, corresponding to a 0.04 Hz test frequency.

300 Analysis-Heart Analysis-Edge FEM-Heart FEM-Edge

σyy (MPa)

200

100

0

-100

-200

0

25

50

75

100

time (s) Fig. 8. Stress evolution induced by cyclic temperature evolution (Ti = 350 °C) on the edge (x = L) and in the heart (x = 0) of the specimen.

6

Analysis FEM

4

max

KI

K I (MPa.m1/2)

2

0

-2 min

KI

-4

-6

0

25

50

75

100

temps (s) Fig. 9. SIF of cracked plate with a = 0.2 mm crack submitted to sinusoidal fluctuation of temperature.

2362

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

The evolution of the temperature on the edge of the specimen is compared with that in the heart (located in the middle of the thickness along x of the coupon (x = 0)) in Fig. 7. The agreement between analysis and ABAQUS simulation on the crackfree body is fairly good. An initial small transient evolution is observed, limited here to the first cycle as the calculations were done for a low frequency, followed by a steady one. For higher values of the frequency, the number of cycles necessary to reach the steady state is higher. As could be expected, the plotted curves show that the amplitude of fluctuation is less important in the heart than on the surface of the plate. In addition, a delay on the cyclic temperature in the heart (x = 0) is observed, due to the conduction phenomenon in the specimen thickness. This cyclic gradient in temperature induced is therefore responsible of cyclic mechanical stresses in the plate, as demonstrated earlier through Eq. (7). On the surface (x = L), the steady evolution of the temperature belongs in the range [140–309 °C] (amplitude of 169 °C). The amplitude of the steady evolution decreases when the test frequency increases. For example, for a 1 Hz frequency, the edge temperature fluctuates between 197 °C and 252 °C (amplitude of 55 °C). The evolution of normal stress on the edge and in the heart of the coupon is plotted in Fig. 8. As a consequence of the temperature evolution, the normal stress also has a transient and a steady evolution. Edge and heart normal stresses appear to be almost out-of-phase: when the edge is under tension, the heart is under compression, and inversely. Under steady regime, the maximum and minimum stress values are opposite. The stress amplitude is lower in the heart (x = 0) than on the edge of the plate (x = L). The very good agreement with numerical results on the crack-free body is also noted. The stress intensity factor of the plate with a 0.2 mm surface crack is now plotted versus time in Fig. 9. This evolution also presents a short transient and a steady regime. When the surface of the plate is under compression, the SIF has a negative value, illustrating the occurrence of crack closure during part of the cycle. Formin fatigue study, the only steady response will be considered. It is noteworthy that during the steady regime, the SIF ratio KI RK ¼ K max is equal to 1 because of the observed symmetry in the stresses. Thus, our attention will particularly focus on the I K max value during the steady regime. I As widely accepted under fatigue loading, we assume now that, in the range [DKth–DKu], cracks obey a classical Paris law

da ¼ C DK n dN

ð17Þ

da where dN is the crack growth rate, DK the stress intensity factor range, C and n are material parameters. As here K is sinusoidal, we consider that DK = Kmax. Below the threshold value DKth, cracks are assumed to stop, and above the ultimate value DKu they propagate instantaneously.

3.3.2. Influence of different parameters on the SIF We wish in this part to investigate the influence of several parameters on the calculated stress intensity factor. As the agreement is good between analysis and ABAQUS simulation for temperature, stress and SIF, we may restrict to analytical calculations. As the SIF is simply related to the stresses through Eq. (16), when possible, we will justify these influences using the analytical expression of the stress obtained in Eq. (13) that can clearly be divided into a transient:

Ea

1 X

! 2

cn t

Wn ðxÞðT i  B  A sin ðun ÞÞe

ð18Þ

n¼0

and a steady part

Ea

1 X

!

Wn ðxÞ  A cosðun Þ cosðxt  un Þ

ð19Þ

n¼0

The initial temperature Ti of the coupon only influences the transient regime, and not the steady one. This is illustrated in Fig. 10, with two values of Ti (Ti = 350 °C and Ti = 500 °C) for Bi = 1.57: the initial value of the SIF increases with Ti. The steady evolution depends only on the characteristics of the convective medium. The influences of the mean temperature B and of the amplitude A of the convective medium on the SIF evolution are illustrated in Fig. 11. The mean temperature B influences only the transient component of the SIF as demonstrated in Eqs. (18) and (19). The amplitude A influences both the transient and the steady responses. In particular, the steady SIF values are proportional to A. The frequency f of the test influences the transient and steady regimes through the phase-shift un defined in Eq. (14). The Biot number Bi also influences the two regimes through un and kn. As the corresponding influences are complex, and also because we are interested in fatigue study, we will focus on the steady regime, and especially on the K max maximum SIF vaI lue (Fig. 12, plotted for a = 0.2 mm crack length). In the studied range of Biot number, it is found that Bi does not change significantly the shape of the K max curve. For a given frequency, the K max value increases with the Bi value: a crack will be more I I likely to propagate under severe thermal shock. For a given Biot number, K max has a very sharp increase for the very small values of the frequency. Then, K max slowly tends I I towards zero for higher frequencies, as the temperature of the coupon undergoes very limited fluctuations. The peek is attained respectively at 0.1, 0.18, and 0.2 Hz for a Biot number of 1.57, 3.15, and 4.47. For Bi = 1.57 and a = 0.2 mm, the curve K max is below the threshold value Kth for every frequencies. For Bi = 3.15 and Bi = 4.47, below a given frequency, crack I

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

7

2363

Analysis-Ti=350°C Analysis-Ti=500°C FEM-Ti=350°C FEM-Ti=500°C

6 5 4

KI(MPa.m1/2)

3 2 1 0 -1 -2 -3 -4 -5

0

25

50

75

100

time (s) Fig. 10. SIF evolution on the edge of specimen for Ti = 350 °C and Ti = 500 °C, with Bi = 1.57, a = 0.2 mm.

x

6 5

xxx x x x x x x x

4

1/2

KI(MPa.m )

3 2 1

x

0 xx

xxx x x x x x x

x

x

-1

x

-2

x

-3

x

x x xxx

-4 -5

x

0

25

A=125-B=225 A=75-B=225 A=125-B=125

xxx x x x x x x

x

x

x x

x x

xxx x x x x x x

x

x

x x

x x

x x xxx

x x xxx

50

75

x x x x x xxx 100

time (s) Fig. 11. Influence of the mean temperature B and of the amplitude A of the convective medium, for Bi = 1.57, a = 0.2 mm.

propagation can occur as K max is higher than Kth. It is also noted that the shapes of these curves are analogous for other crack I lengths.

3.3.3. Evolution of SIF with crack length It is interesting now to investigate the evolution of K max versus crack length for prediction of fatigue propagation. The I maximum stress intensity factor in the steady regime is plotted versus crack length in Fig. 13 for different values of the frequency in the range of [0.04–10 Hz], Bi being equal to 1.57. A peek of K max is observed for a crack length value decreasing I with increase in the frequency: longer cracks are more sensitive to small frequencies. A comparison with the results of Taheri [30] shows similarities in the tendencies. With the assumptions presented earlier concerning the propagation law of the material, cracks cannot propagate with frequencies higher than 1 Hz because the curve K max lies under the threshold curve. When f = 0.5 Hz, a crack arrest should I be observed at 2 mm. For lower frequencies (f = 0.04 and 0.1 Hz), the K max curve is highly above the threshold curve: through I predictions, the crack with these conditions should propagate under thermal cycling until it reaches a crack length of 4 mm. The crack growth rate would be maximal for a crack length around 1.5 mm.

2364

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

9 8

Bi=4.47

KI

max

1/2

(MPa.m )

7 th

K

6 5 4

Bi=3.15

3

Bi=1.57

2 1 0

2

4

6

8

10

Frequency (Hz) Fig. 12. Maximum SIF K max versus the loading frequency for different Biot numbers, crack length a = 0.2 mm. I

Fig. 13. Maximum SIF K max versus crack length for different loading frequencies, with Bi = 1.57. I

3.4. Thermal fatigue by saw-toothed evolution of the boundary temperature conditions 3.4.1. Evolution versus time of temperature, stress and SIF We consider in this paragraph the specimen presented in Fig. 5 submitted to a convective medium with a saw-toothed variation of temperature. It is supposed that this temperature varies between Tmax = 350 °C and Tmin = 100 °C. The cooling– do heating time ratio Rt is defined as: Rt ¼ ttup . Rt = 1 corresponds to the saw-toothed fluctuation with the same cooling and heating rates (triangular fluctuation) and Rt = 0 to an ideal repeated cold thermal shock. In order to approach cold shock experimental conditions [8,13,25], we will restrict to 0 < Rt < 1. The evolution of the temperature in the heart of the specimen for current analysis and FEM approach is plotted in Fig. 14. The curves have been determined with Bi = 1.57 as in sinusoidal case. The other parameters for calculation are tdo ¼ 3 s; t up ¼ 197 s implying Rt = 0.0151. A very good agreement between the two curves allows validating the analytical solution for saw-toothed evolution of boundary condition. With the retained parameters, no transient evolution can be noted: the steady response is obtained from the first cycle. The difference of temperature between the edge and the heart of the specimen is quite small (below 84 °C at the end of cooling t = 4 s). During the heating period, a homogeneous temperature is observed in the specimen. This can be explained by the high value of tc and the related slow temperature rates, enabling conduction phenomena. The maximum and minimum temperatures are then very close to the imposed ones. A small delay (about 2 s) in the heart of the coupon is also observed.

2365

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

350 Analysis-Heart Analysis-Edge Difference Heart-Edge analysis FEM-Heart FEM-Edge Difference Heart-Edge FEM

Temperature (°C)

300 250 200 150 100 50 0 0

25

50

75

100

time (s) Fig. 14. Evolution of temperature under saw-toothed fluctuation, tc = 200 s, Rt = 0.015.

350

350

300

300 Analysis-Stress in heart

250

250

Analysis-Stress on edge Temperature difference heart-edge FEM-Stress in heart

σyy (MPa)

150

150

FEM-Stress on edge

100

100

σ yy Edge

50

50

0

0

-50

-50

Temperature difference

-100

-100

σ yy Heart

-150 -200

200

Temperature (°C)

200

-150 0

25

50

75

-200 100

time (s) Fig. 15. Normal stress evolution induced by saw-toothed fluctuation for Ti = 350 °C.

Then, the ryy stresses are plotted on the edge and in the heart of the coupon in Fig. 15. As for temperature, no transient regime appears. As under sinusoidal loading, we note out-of-phase evolution between edge and heart. The peeks of stresses are observed at the end of the cooling period. During the heating period, the stresses approach zero as the temperature is nearly homogeneous in the specimen. In the same figure, the temperature difference between heart and edge of the coupon is also plotted. Despite the small temperature gradient in the specimen thickness, the ryy stresses reach important values: 300 MPa on the surface and -148 MPa in the heart. It is also remarked that the peeks of stress are obtained when the difference in temperature is maximum. Fig. 16 presents the corresponding evolution of the SIF versus time. A good agreement in tendency between analytical approach and simulation is noted. However, a discrepancy (13%) at the peek of SIF is recorded. This discrepancy decreases with increasing crack length. For example the discrepancy is inferior to 1% for a = 1 mm. The maximum value of the SIF is obtained at the end of the cooling period. As the heating period is long here, a saturation of the SIF value is observed. The SIK ratio RK is no longer here equal to 1: with the conditions of the test, RK = 0.036. It appears to be close to Rt = 0.015.

2366

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

10

KI(MPa.m1/2)

8

max

KI

6 Analysis 4

FEM

2 Kmin I

0 0

25

50

75

100

time (s) Fig. 16. SIF of cracked plate with a = 0.2 mm subjected to saw-toothed fluctuation of temperature (Ti = 350 °C).

14 12

T i=500°C

10

10

8

1/2

KI(MPa.m )

1/2

KI(MPa.m )

12

T i=350°C

6 4 2

8

6

4

0 0

100

time (s)

200

300

Fig. 17. Influence of the initial temperature on the evolution of the SIF subjected to saw-toothed fluctuation of temperature (a = 0.2 mm, Bi = 1.57).

pffiffiffiffiffi For saw-tooted fluctuation, the maximum value of SIF is K max ¼ 7:75 MPa m. This value is much higher than the maximum I pffiffiffiffiffi pffiffiffiffiffi SIF for sinusoidal fluctuation of the same frequency (0:86 MPa m for f = 0.005 Hz) or of the same cooling rate (4:9 MPa m for f = 0.166 Hz). Hence, the shape of the temperature fluctuation plays a crucial role for fatigue analysis.

3.4.2. Influence of different parameters on the SIF The influence of the initial temperature Ti is illustrated in Fig. 17. The comments are analogous to those previously developed under sinusoidal thermal loading: the initial temperature only plays a role in the transient regime (cf. Eqs. (A4) and (A7)), with a higher SIF value because of the initial cold shock when Ti = 500 °C. The maximum stress intensity factor K max is plotted versus the frequency f for different values of the cooling–heating ratio I Rt in Fig. 18a, with a zoom for smallest frequencies in Fig. 18b. For the small values of f, an initial sharp increase of K max is I obtained, up to a peek, followed by a slow decrease down to 0. The peek value is higher for small values of Rt, corresponding to more severe cooling conditions. The frequency corresponding to the peek increases with Rt, but remains lower than 0.1 Hz. Above 0.1 Hz, K max increases with Rt. I With the considered parameters of the saw-toothed fluctuation, crack propagation should occur only for Rt = 0.015 and at low frequency (K max > K th ). For high frequency, even for the most severe saw-toothed fluctuation, K max always remains beI I low the threshold value. The curve corresponding to the sinusoidal fluctuation of temperature appears to be close to that of Rt = 1, demonstrating that for the same heating and cooling rate, the shape of the temperature (triangle or sinusoidal) has a limited influence on K max . I

2367

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

Rt=0.015 Rt=0.2 Rt=0.5 Rt=1.0 Sinusoidal K th

(a) 8 7

8

6 1/2

5

max

4 3

5 4 3

2

2

1

1 0.5

1

Rt=0.015 Rt=0.2 Rt=0.5 Rt=1.0 Sinusoidal th K

7

K I (MPa.m )

K max (MPa.m1/2) I

6

(b)

1.5

2

0.1

0.2

Frequency (Hz)

0.3

0.4

0.5

Frequency (Hz)

Fig. 18. a and b – Maximum SIF K max versus frequency for different values of the cooling–heating ratio Rt, a = 0.2 mm, Bi = 1.57.

16 f=0.005Hz

14 12

f=0.04Hz

KI

max

10 f=0.1Hz

8

th

K

6

f=0.4Hz

4

f=1Hz

2 f=2Hz

0

0

1

2

3

4

5

a (mm) Fig. 19. K max versus crack lengths for different loading frequencies and crack lengths, Bi = 1.57, Rt = 0.015.

3.4.3. Evolution of SIF with crack length The evolution of K max versus crack length for different frequencies with Rt = 0.015 is presented on Fig. 19. Qualitatively, I this curve presents similarities with the sinusoidal case (cf. Fig. 13): longer cracks are more sensitive to small frequencies. Quantitatively, above 0.4 Hz, no crack propagation is expected with Rt = 0.015 as K max remains below the threshold Kth. BeI low this frequency, crack propagation will be possible, until the crack length reaches approximately 4 mm. The lowest is the frequency, the easiest is the crack propagation. 4. Conclusions An efficient analytical approach based on the weight function method and on the Duhamel’s principle has been developed in order to calculate the time-dependant stress intensity factor of a plate with an edge crack submitted to thermal cyclic loading. This analysis agrees very well with ABAQUS numerical simulation. The main advantage of the present approach is that parametric study can be easily performed. Two types of thermal cycling have been investigated:  Sinusoidal temperature fluctuation.  Saw-toothed temperature fluctuation.

2368

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

The following main conclusions can be drawn:  Cyclic temperature evolution on the surface of the plate causes a cyclic gradient of temperature between the surface and the heart of the plate. This gradient is the origin of out-of-phase cyclic normal stress evolution between the surface and the heart of the considered plate.  The SIF value is proportional to the amplitude Tmax  Tmin of the applied temperature of the convective medium. Whatever the type of temperature loading, an increase in the Biot number implies a raise of the SIF value. Moreover, as the initial temperature of the plate only influences the transient evolution of the SIF, it can have only a very slight role on the fatigue behavior, in particular when the frequency of the thermal cycling is low.  A smaller thermal loading frequency leads to easier crack propagation. Nevertheless, for any frequency, the crack growth rate first increases before undergoing a sharp decrease that may lead to crack arrest, through considerations based on a SIF threshold value for the material.  Besides frequency, Biot number, crack length. . ., it is necessary to take into account the shape of thermal loading in fatigue analysis. In particular, the cooling–heating time ratio Rt plays a major role on kinetics of propagation. Please note that a three-dimensional extension of the analysis developed in this paper is presented and compared with experiments and finite element simulation in [32], but without the possibility of studying the influence of several parameters has it has been done here due to computational cost. A fairly good agreement in kinetics prediction is obtained between three-dimensional extension of the analysis and finite element simulation as well as experimental results. Appendix A

1. For 0 6 t 6 tc (Nc = 0):  if t 6 tdo

rIyy ðx; tÞ ¼ Ea

1 X



Wn ðxÞ  ecn t ðT i  T max þ

n¼0

 if tdo 6 t 6 tc

rIyy ðx; tÞ ¼ Ea

1 X

1

cn



Wn ðxÞ  ecn t ðT i  T max Þ þ

n¼0

2. For tc 6 t (Nc P 1):  If Nctc 6 t 6 Nctc + tdo: I yy ðx; tÞ

r

¼ Ea

1 X

" cn t

Wn ðxÞ e ðT i  T max Þ þ

n¼0

 ð1  ecn t Þ

1

cn

I

cn

ðA1Þ

ðecn ðttdo Þ  ecn t Þ þ

ðe

cn ðttdo Nc tc Þ

M

cn

ð1  ecn ðttdo Þ Þ

cn ðtNc tc Þ

e

Þþ

I cn tup cn ðe

 ðA2Þ

 ecn tc Þ þ cM ð1  ecn tup Þ n

ecn tc  1

# ðe

cn t

e

cn ðtNc tc Þ

Þ

ðA3Þ

Hence, the transient and the steady components of opening stress noted respectively Nctc 6 t 6 Nctc + tdo: I yytrans ðx; tÞ

r

¼ Ea

1 X

"

Wn ðxÞe

n¼0

rIyystab ðx; tÞ ¼ Ea

1 X

"

Wn ðxÞ

n¼0

 if Nctc + tdo 6 t:

rIyy ðx; tÞ ¼ Ea

1 X

cn t

I

cn

ðT i  T max Þ þ



n¼0

þ

cn tup

 ecn tc Þ þ cM ð1  ecn tup Þ

I

cn

 ecn tc Þ þ cM ð1  ecn tup Þ n

ecn tc  1

#

n

ðA4Þ

ecn tc  1

ðecn ðttdo Nc tc Þ  ecn ðtNc tc Þ Þ 

Wn ðxÞ ecn t ðT i  T max Þ þ

I cn tup cn ðe

I

cn ðe

rIyytrans ðx; tÞ and rIyystab ðx; tÞ for

I cn tup cn ðe

 ecn tc Þ þ cM ð1  ecn tup Þ n

ecn tc  1

ðecn ðtNc tc tdo Þ  ecn ðtNc tc Þ Þ þ #

M

cn

# ecn ðtNc tc Þ

ðA5Þ

ð1  ecn ðt  Nc t c  tdo ÞÞ

ðecn t  ecn ðtNc tc Þ Þ

ðA6Þ

Hence, the transient and the steady components of opening stress are the following for Nctc + tdo 6 t:

rIyytrans ðx; tÞ ¼ Ea

1 X n¼0

"

Wn ðxÞecn t ðT i  T max Þ þ

I

cn ðe

cn tup

 ecn tc Þ þ cM ð1  ecn tup Þ n

ecn tc  1

#

ðA7Þ

H.N. Le, C. Gardin / Engineering Fracture Mechanics 77 (2010) 2354–2369

rIyystab ðx; tÞ ¼ Ea

1 X



Wn ðxÞ

n¼0



I cn tup cn ðe

I

cn

ðecn ðtNc tc tdo Þ  ecn ðtNc tc Þ Þ þ

 ecn tc Þ þ cM ð1  ecn tup Þ n

ecn tc  1

# cn ðtNc tc Þ

e

M

cn

2369

ð1  ecn ðt  Nc t c  t do ÞÞ ðA8Þ

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]

Guédé Z, Sudret B, Lemaire M. Life-time reliability based assessment of structures submitted to thermal fatigue. Int J Fatigue 2007;29(7):1359–73. Velasco E, Colás R, Valtierra S, Mojica JF. A model for thermal fatigue in an aluminium casting alloy. Int J Fatigue 1995;17(6):399–406. Kikuchi K, Ue K, Kudo Y, Saito M. Crack propagation in first wall by thermal fatigue and creep. Fusion Engng Des 2000;49–50:229–34. Kudo Y, Kikuchi K, Saito M. Thermal fatigue crack propagation behaviour of F82H ferritic steel. J Nucl Mater 2002;307–311:471–4. Haddar N, Fissolo A. 2D simulation of the initiation and propagation of crack array under thermal fatigue. Nucl Engng Des 2005;235(9):945–64. Kimura S, Kogawa H, Teramoto T, Saito M. Crack propagation in first wall of fusion reactor by cyclic thermal stress. Fusion Engng Des 1998;39– 40:569–74. Ancelet O, Chapuliot S, Hénaff G. Experimental and numerical study of crack initiation and propagation under 3D thermal fatigue loading in a welded structure. Int J Fatigue 2008;30(6):953–66. Ancelet O, Chapuliot S, Henaff G, Marie S. Development of a test for the analysis of the harmfulness of a 3D thermal fatigue loading in tubes. Int J Fatigue 2007;29(3):549–64. Fissolo A, Amiable S, Ancelet O, Mermaz F, Stelmaszyk JM, Constantinescu A, et al. Crack initiation under thermal fatigue: an overview of CEA experience. Part I. Thermal fatigue appears to be more damaging than uniaxial isothermal fatigue. Int J Fatigue 2009;31(3):587–600. Lejeail Y, Kasahara N. Thermal fatigue evaluation of cylinders and plates subjected to fluid temperature fluctuations. Int J Fatigue 2005;27:768–72. Kerezsi B, Price JWH, Ibrahim R. A two-stage model for predicting crack growth due to repeated thermal shock. Engng Fract Mech 2003;70(6):721–30. Gardin C, Le HN, Benoit G, Bertheau D. Crack growth under thermal cyclic loading in a 304L stainless steel – Experimental investigation and numerical prediction. Int J Fatigue 2010;32:1650–7. Haddar N, Fissolo A, Maillot V. Thermal fatigue crack networks: an computational study. Int J Solids Struct 2005;42(2):771–88. Seyedi M, Taheri S, Hild F. Numerical modeling of crack propagation and shielding effects in a striping network. Nucl Engng Des 2006;236(9):954–64. Shahani AR, Nabavi SM. Closed form stress intensity factors for a semi-elliptical crack in a thick-walled cylinder under thermal stress. Int J Fatigue 2006;28(8):926–33. Svivastava A, Joshi V, Shivpuri R. Computer modeling and prediction of thermal fatigue cracking in die-casting tooling. Wear 2004;256:38–43. Miranda ACO, Meggiolaro MA, Castro JTP, Martha LF, Bittencourt TN. Fatigue life and crack path predictions in generic 2D structural components. Engng Fract Mech 2003;70:1259–79. Lu TJ, Fleck NA. The thermal shock resistance of solids. Acta Mater 1998;46(13):4755–68. Marie S. Analytical expression of the thermal stresses in a vessel or pipe with cladding submitted to any thermal transient. Int J Press Vessels Pip 2004;81(4):303–12. Marie S, Chapuliot S. Improvement of the calculation of the stress intensity factors for underclad and through-clad defects in a reactor pressure vessel subjected to a pressurised thermal shock. Int J Press Vessels Pip 2008;85(8):517–31. Marie S, Ménager Y, Chapuliot S. Stress intensity factors for underclad and through clad defects in a reactor pressure vessel submitted to a pressurised thermal shock. Int J Press Pip 2005;82:746–60. Kim YW, Kim JI, Chang MH. An explicit integral expression for the stress intensity factor of a semi-elliptic surface crack subjected to thermal transient loading. Int J Press Vessels Pip 1999;76(9):631–9. Lee HY, Kim JB, Yoo B. Green’s function approach for crack propagation problem subjected to high cycle thermal fatigue loading. Int J Press Vessels Pip 1999;76(8):487–94. Arpaci VS. Conduction heat transfer. Addison-Wesley Publishing Company; 1966. Fissolo A, Marini B, Nais G, Wident P. Thermal fatigue for a 316L type steel. J Nucl Mater 1996;233–237:156–61. Le HN, Gardin C, Benoit G, Bertheau D. Crack growth estimation in a parallelepipedic specimen subjected to a cyclic temperature gradient. In: Portella PD, Beck T, Okazaki M, editors. 6th International conference on low cycle fatigue. Berlin [Germany]: DVM; 2008. p. 183–8. Tada H, Paris PC, Irwin GR. The stress analysis of cracks handbook. Paris Productions Incorporated; 1985. Holman JP. Heat transfer. Tokyo: McGraw-Hill Kogakusha; 1976. Courtin S, Gardin C, Bézine G, Ben Hadj Hamouda H. Advantages of the J-integral approach for calculating stress intensity factors when using the commercial finite element software ABAQUS. Engng Fract Mech 2005;72(14):2174–85. Taheri S. Some advances on understanding of high cycle thermal fatigue crazing. J Press Vessels Technol 2007;129(3):400–10. Maillot V, Fissolo A, Degallaix G, Degallaix S. Thermal fatigue crack networks parameters and stability: an experimental study. Int J Solids Struct 2005;42(2):759–69. Le HN, Gardin C. Analytical prediction of crack propagation under thermal cycling inducing a thermal gradient in the specimen thickness – comparison with experiments and numerical approach. Engng Fract Mech, submitted for publication.