Thermo-mechanical stress intensity factors evaluation

Nuclear Engineering and Design 240 (2010) 2579–2588

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Thermo-mechanical stress intensity factors evaluation Thomas Menouillard ∗ Theoretical and Applied Mechanics, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208-3111, USA

a r t i c l e

i n f o

Article history: Received 16 March 2010 Received in revised form 7 June 2010 Accepted 21 June 2010

a b s t r a c t This paper deals with the evaluation of the static mixed-mode stress intensity factors in case of thermal and mechanical loadings. A method and a formulation based on the J-integral and the energy release rate formulation is developed. The numerical results emphasize the effects of the temperature field on the stress intensity factors; thus, we show that an increasing temperature in the direction of the crack increases the mode 1 stress intensity factor. In addition, the formulation is also applied to the computation of an experiment consisting in a precracked glass specimen into a temperature gradient induced by a water bath and an oven. © 2010 Elsevier B.V. All rights reserved.

1. Introduction The cooling process and the fission in nuclear power plants induce stress state in many components, such as the tank and pipelines. Thermal loading (shock and fatigue) (James, 1976; McGowan, 1979; Lindau and Möslang, 1994; Chapuliot et al., 2005; Reytier et al., 2006; Hasselman, 2006) has become a wide research area according to the literature (Noda and Jin, 1993; Erdogan and Wu, 1996; Janssens et al., 2009) because the temperature gradient affects the behavior of existing microcracks and cracks, and thus the stress state may become severe, and because some nuclear power plants are already becoming 30–40 years old. The computational methods has to give tools to be able to predict whether there are some risks in using these power plants by computing fracture mechanics. Thus the evaluation of the stress intensity factors is major point in the different materials with varying properties in space and with thermal effects. The computation of the stress intensity factors (Rice, 1968a) for functionally graded materials (FGM), materials whose properties vary spatially, was developed by Kim and Paulino (2002, 2003a,b) and Rao and Rahman (2003). FGM represent a wide research field linked to metals and metal-ceramic composites, as described by Markworth et al. (1995), Gu and Asaro (1997) and Suresh and Mortensen (1997). Noda and Jin (1993) dealt with the stress intensity factors in this particular case; their analysis was restricted to the assumption that all material properties depend only on the coordinate y (perpendicular to the crack faces) in such a way that the properties are exponential functions of y. Erdogan and Wu (1996) developed formulas for the thermal stress intensity factors for Functionally Graded Material in the thickness of the struc-

∗ Tel.: +1 847 491 3046. E-mail address: [email protected]. 0029-5493/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2010.06.034

ture. Jin and Paulino (2001) studied an edge crack in a strip of Functionally Graded Material (FGM) under transient thermal loading; the material is assumed to have constant Young’s modulus and Poisson’s ratio, but the thermal properties varied along the thickness direction of the strip, which corresponded to the crack direction. This paper presents a development for evaluating the mixedmode stress intensity factors when there is a temperature field in the structure and space variable material properties. This variation may be due to temperature effects on material properties such as Young’s modulus. On contrast to Menouillard et al. (2006), which considered mixed mode for graded materials, the effect of the temperature on the fracture parameters, the energy release rate and the stress intensity factors, is developed. In our analytical development, the temperature field is assumed to linearly vary in the vicinity of the crack tip. This paper is organized as follows. Section 2 presents the static thermo-mechanical problem and the governing equations corresponding to thermal and mechanical fields, such as temperature, stress and displacement. Section 3 presents the form of these three different fields near the crack tip. Section 4 presents the energy release rate formulation taking into account both thermal and mechanical effects and also develops a formulation and a strategy to determine the two stress intensity factors with such thermo-mechanical effects. Section 5 presents the numerical examples: first, a validation example compares the theoretical mode 1 stress intensity factor to the computation. Second, a mixed-mode specimen under a fixed mechanical loading and different thermal loadings; this allows to emphasize the effect of the temperature on the value of the energy release rate and the stress intensity factors with the same mechanical loading. Third, an experiment deals with a quasi-static crack propagation induced by thermal loading. Section 6 presents the concluding remarks.

2580

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

Nomenclature

The mechanical boundary conditions on the displacement and tracand ∂M ) of the domain are: tion boundaries (∂M 1 2

C0

u = ud

C E, E0 ex , ey er , e fd , F d G Id K1 , K2 n r,  T T,x T,y u u0 us x x, y

Linear elastic behavior law for constant mechanical properties in space Linear elastic behavior law for varying mechanical properties in space Young’s modulus Cartesian basis near crack tip Polar basis near crack tip External forces Energy release rate Identity tensor Stress intensity factors Outward normal Polar coordinates Temperature Temperature gradient in ex direction Temperature gradient in ey direction Displacement field Asymptotic displacement for constant mechanical properties in space Asymptotic displacement for varying mechanical properties in space Position vector Cartesian coordinates

Greek letters ˛ Coefficient of thermal expansion ε Strain tensor  Thermo-mechanical coefficient Crack c  Kolosov’s coefficient ∇ Gradient operator  Poisson’s ratio  Domain Stress tensor Direction of propagation c  Virtual extension field

Let us consider a body  with a crack c shown in Fig. 1. The body is assumed to be a linear isotropic elastic material whose Young’s modulus varies continuously in space. Bimaterials are not considered in this paper. Static thermal and mechanical loads are applied on the body; these lead to a equilibrium defined by the displacement field u(x), the temperature T (x), and the stress (x). x denotes the location of a material point. Vectors are denoted with bold characters, and tensors are underlined two times. A prescribed displacement ud is applied on the boundary ∂M and 1 M and ∂M are such external forces Fd on the boundary ∂M . ∂ 2 1 2 ∪ ∂M = ∂M and ∂M ∩ ∂M = ∅. We do not consider that ∂M 1 2 1 2 body forces fd in . In addition, the body is subjected to thermal loading; a temperature Td is applied on the boundary ∂T1 , and a heat flux d on the boundary ∂T2 . ∂T1 and ∂T2 are such that ∂T1 ∪ ∂T2 = ∂T , and ∂T1 ∩ ∂T2 = ∅. ∂M is used for the mechanical boundary condition, whereas ∂T is for the thermal boundary conditions. The mechanical and thermal equilibrium equations in the domain  are: in 

∇ 2 T = 0 in 

· n = Fd ·n=0

on

(3)

∂M 2

(4)

on crack surfaces

(5)

where n is the outward normal to the material boundary. The thermal boundary conditions on the different parts ∂T1 and ∂T2 are: T = Td

on ∂T1

(6)

∇ T · n = d on ∂T2

(7)

Fig. 1 presents a schematic of the mechanical and thermal boundary conditions. The behavior law, denoted by C, describes the relationship between the stress , the displacement u and the temperature T as: = C · (ε(u) − ˛ T Id )

in 

(8)

where ˛ denotes the coefficient of thermal expansion, and Id the identity matrix. The behavior law C depends on the space or on the temperature as C (T ) or C (x), or even on both. This is due to the Young’s modulus dependency with the space or/and the temperature. Thus the general behavior law form is C (T (x), x), which can be recast in C (x), and be finally considered for computation. The method to solve the problem is developed as follows. First of all, the temperature field can be extracted and evaluated using a thermal analysis. Temperature fields encountered in a lot of structures under thermal boundary conditions are linear (or quasilinear). Even though, the continuous temperature field can be linearized at least near the crack tip vicinity. Second, the shape of the stress can be determined from Eq. (1). The displacement field is solution of a differential equation with a nonzero right term due to the thermal effect; this means that a particular solution has to be found. The shape of the displacement u and stress solution fields are written as a linear combination of the asymtotic fields. The constants appearing in the combination characterize the fractures modes, called the stress intensity factors. 3. Solution fields near crack tip

2. Governing equations

div = 0

on ∂M 1

(1) (2)

3.1. Temperature The assumption that the temperature field has the property to be linearly dependent of the material position near the crack tip vicinity, is made. Thus, the temperature field near crack tip is then defined by: T (x, y) = T,x x + T,y y + T0

(9)

where T,x , T,y and T0 are constants whose values are determined from the thermal problem described by Eqs. (2), (6) and (7). Fig. 2 presents the different basis and coordinates system; the origin is chosen to be the crack tip, and the direction ex is given by the direction of the crack. (x, y) denotes the coordinates of a point M located at position x as shown in Fig. 2. Solving the thermal problem alone leads to the determination of the temperature field in the structure, which allows to determine the constants T,x , T,y and T0 expressed in Eq. (9). 3.2. Stress The stress is expressed in Eq. (1), which is a system of firstorder differential equations over the continuous domain , whose solution involves basis functions, denoted by g1 and g2 defined over

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

2581

Fig. 1. Domain  and (a) mechanical and (b) thermal boundary conditions.

the domain , and some constants denoted by K1 and K2 , which are given by Williams (1957) and Williams and Ewing (1972). The stress field in the vicinity of the crack tip (r is small) is written as: lim (x) = K1 g1 (x) + K2 g2 (x)

r→0

(10)

where the tensors g1 and g2 can be viewed as the asymptotic stress fields near the crack tip, corresponding respectively to modes 1 and 2 (Bui, 1978). Detailed expressions of g1 and g2 are given in Appendix A.1, in the polar coordinates (see Fig. 2 for nomenclature). The scalars K1 and K2 denote the stress intensity factors corresponding respectively to modes 1 and 2. So far, no reference has been made to the elastic linear constitutive relation, neither the space dependence of the Young’s modulus. Therefore, these analytical results are still true for any spatial distribution of the behavior law, provided that it is continuous over the domain  and provided that r is sufficiently small. 3.3. Displacement As Eq. (10) remains true for any linear elastic law, one can write for a linear elastic material with constant mechanical properties: (x) = C0 · ε(u0 (x))

C

· ε(u0i (x))

= gi (x) ∀i ∈ {1, 2}

(12)

However, the constitutive relation C of our problem varies in space. By combining the behavior law (Eq. 8) and the stress shape near crack tip (Eq. 10), leads to the following equation verified by the displacement u: C (x) · (ε(u(x)) − ˛ T (x)Id ) = K1 g1 (x) + K2 g2 (x)

u(x) = K1 us1 (x) + K2 us2 (x) + uT (x)

(14)

where us1 and us2 are the asymptotic displacement fields near the crack tip corresponding to the material model C(x), and uT a particular displacement solution of Eq. (13). Menouillard et al. (2006) showed that the asymptotic displacement fields us1 and us2 for the material with varying characteristics in space are related to the displacement fields u01 and u02 and to the constitutive relations C (x) and C0 , as: ∀i ∈ {1, 2} ε(usi ) = (C (x))−1 · (C0 · ε(u0i ))

(15)

Eq. (15) established in (Menouillard et al., 2006), allows to evaluate the asymptotic strain field ε(usi ) for any continuously varying graded material using the known asymptotic strain field of a constant mechanical property. It then allows to evaluate the asymptotic displacement field usi from the knowledge of u0i and C (x) through an analytical integration of the strain field. The displacement field uT denoted by (ur , u ) in the basis (er , e ) due to the thermal term ˛ T (x), can be extracted from Eqs. (12)–(15) as:

(11)

where u0 is the displacement field near the crack tip. The 0 exponent corresponds to a linear elastic material with constant properties. According to Eqs. (10) and (11), the displacement field u0 (x) can be written as K1 u01 (x) + K2 u02 (x) near the crack tip, where u01 and u02 are defined by: 0

Thus, the displacement field u, solution of Eq. (13), takes the following form:

(13)

ε(uT (x)) = ˛ T (x)Id

(16)

The particular displacement field uT needs to be determined. Solving Eq. (14) with the form T (x, y) = T0 + T,x x + T,y y for the temperature field, leads to the solution (see Appendix B):

⎛ uT = ⎝

⎞

˛

r2 (T,x cos() + T,y sin()) + ˛T0 r 2 ⎠ r2 ˛ (T,x sin() − T,y cos()) 2 (er ,e )

(17)

Thereby the associated displacement uT is completely determined by the thermal gradient. 4. Mixed-mode analysis 4.1. Energy release rate formulation

Fig. 2. Cartesian and polar coordinates used near crack tip.

The J-integral developed by Rice (1968b) is the original way to define the energy release rate. A virtual extension field  which comes form the path independent formulation of the J-integral (Rice, 1968a) is used in the formulation of the energy release rate. Fig. 3 presents the continuous virtual extension field  parallel to the crack face, whose norm decreases from 1 at the crack tip to zero far enough, as shown in Fig. 3. The general formulation of the energy release rate G for space variable materials and thermal effects comes from both studies of Suo and Combescure for nonhomogeneous materials (Suo and Combescure, 1992a) and thermal

2582

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

Fig. 3. Description of the field  near the crack’s tip.

• and m is defined by the following scalar:

effects (Suo and Combescure, 1992b), as:



G





tr( ∇ u ∇ ) d −

=





w div() d 

tr( ) ∇ (˛ T ) ·  d −

+





1 2

m

=

(˛ T )2 div()tr(C)

− 

+2∇ (˛ T ) · tr(C) +(˛ T )2 tr(∇ C · ) d

(18)

(22)

tr[∇ C (ε(u) − ˛ T Id)(ε(u) − ˛ T Id)] d







+ 

Note that all the forms (a, l and m) own the gradient of the behavior law, denoted by ∇ C, and the form a does not include any thermal effect due the temperature T.

∇ fd  · u d

fd · u div() d + 

where the work w is defined by: w=

1 tr[ · (ε(u) − ˛ T Id)] 2

(19)

and div() = (∂ x /∂x) + (∂ y /∂y). Here, the general formulation of the energy release rate G is shown for both variable material and thermal effects, and can be evaluated from the knowledge of the temperature, stress and displacement fields. 4.2. Stress intensity factors formulation

: R4 → R (u, v) → J(u, v) = a(u, v) + l(u + v) + m

• a(·, ·) introduced by Menouillard et al. (2006), is defined by the following bilinear symmetric form:



tr[C ε(u) · ∇ v · ∇  + C ε(v) · ∇ u · ∇ ]

= 



−tr(C ∇ v ∇ ) ˛ T 

+∇ fd  · v + tr[C ε(v)]∇ (˛ T ) ·  +fd · v div() + ˛ T tr[∇ C ·  ε(v)] +˛ T div()tr[C ε(v)] d

=

1 1 a(us1 , us1 )K12 + a(us2 , us2 )K22 2 2 +a(us1 , us2 )K1 K2 +(a(uT , us1 ) + l(us1 ))K1 +(a(uT , us2 ) + l(us2 ))K2 1 + a(uT , uT ) + l(uT ) + m 2

(24)

4.4. Determination of the stress intensity factors

• l(·) is defined by the following linear form: =

(23)

(20)

−tr[(div() C + ∇ C ) ε(u) · ε(v)] d

l(v)

J(u, u) = 2 G = a(u, u) + 2 l(u) + m

G

where:

a(u, v)

These different functions (i.e. a, l and m) have the following property: if u is the real displacement field, then the form J is related to the energy release rate as:

According to Eq. (14) and the properties of the linearity and bilinearity of the different forms a(·, ·) and l(·), the following relation between the energy release rate and the stress intensity factors is obtained:

First, let us introduce the following function: J

4.3. Link to the energy release rate

(21)

In order to determine the two stress intensity factors K1 and K2 , two scalar equations involving the two quantities must be written. The bilinear form a(·, ·) is applied to two couples of displacement fields as follows: a(u − uT , us1 ) = K1 a(us1 , us1 ) + K2 a(us2 , us1 ) a(u − uT , us2 ) = K1 a(us1 , us2 ) + K2 a(us2 , us2 )

(25)

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

2583

The solutions of this linear system of 2 equations and 2 unknowns are the following: K1 =

K2 =

a(us2 , us2 ) a(u − uT , us1 ) − a(us1 , us2 ) a(u − uT , us2 ) a(us1 , us1 ) a(us2 , us2 ) − a(us1 , us2 )2 a(us1 , us1 ) a(u − uT , us2 ) − a(us1 , us2 ) a(u − uT , us1 ) a(us1 , us1 ) a(us2 , us2 ) − a(us1 , us2 )2

(26)

(27)

4.5. Remarks First, if the linear elastic material has constant properties (the related characteristics are denoted with a zero exponent) and if there is no thermal effect, one obtains the following relations: a(u01 , u01 ) =

2(1 − 2 ) E0

(28)

a(u02 , u02 ) =

2(1 − 2 ) E0

(29)

a(u01 , u02 ) = 0

(30)

where E0 defines the Young’s modulus, and  the Poisson’s ratio. Thus, the stress intensity factors are determined by: K1 = K2 =

E0 a(u, u01 )

(31)

2(1 − 2 ) E0 a(u, u02 )

(32)

2(1 − 2 )

This formulation is the particular case studied by Attigui and Petit (1997). Second, Eqs. (26) and (27) are new because of the term uT corresponding to the displacement field due to the thermal loading. However Menouillard et al. (2006) showed that when the Young’s modulus depends on space, the use of the auxiliary displacement fields u01 and u02 instead of us1 and us2 is correct whether a small radius of  field presented in Fig. 3 is chosen. This is due to the continuity of the mechanical properties in space. Practically, the asymptotic displacement field u01 and u02 corresponding to the constant material, are chosen for computation, together with a small radius for the parameter . Third, the fact that the temperature may affect the stress intensity factors is taken into account by the term uT only. Beside this, the temperature only affects the energy release rate, as shown by the additional terms in Eq. (24) compared to the particular linear elastic case without thermal loading known as: G=

1 − 2 E

(K12 + K22 )

in plane strain

(33)

These terms due to uT , l(·) and m, are due to the temperature field T, and are zero when there is no temperature (i.e. T = 0). Finally, the theory developed allows the computation of all quantities of the problem: the temperature, the shape of the stress and the displacement, and the latter part gives the way to compute the constants characterizing the formulation of the stress and displacement fields (i.e. stress intensity factors). Thus, we have defined the entire method to evaluate mechanical, thermal and fracture characteristics. 5. Numerical examples 5.1. Example 1: mode 1 case Here, we compare our method for computing stress intensity factors to a reference solution. The simple example made by Wilson and Yu (1979) has the particularity to have both mechanical and

Fig. 4. Geometry and loadings of the experiment.

thermal loading, and their numerical result is taken as a reference. Fig. 4 presents the geometry, boundary conditions and the loading of the cracked structure. The structure is rectangular, whose width is W and height 2W . The dimension W is taken as 0.1 m, and a denotes the crack length, whose value is defined by W/2. The vertical displacement of the top and bottom are constrained to be zero. The crack tip is located at the center of the structure, and the crack surface is horizontal until the left edge. The temperature field in the structure is imposed to linearly vary from −T0 on the left edge to +T0 on the right edge (T0 = 50 K). Thus, the temperature at the crack tip is zero and is negative where the crack remains, and positive on the right part of the structure. This configuration coupled to the displacement boundary conditions makes the crack open in pure mode 1. So, only the stress intensity factor K1 is interesting to compute. The material properties are: Young’s modulus E = 75 GPa, Poisson’s ratio  = 0.3, coefficient of thermal expansion ˛ = 26 ␮m/(m K−1 ). Table 1 presents the results of the normalized stress intensity factor K1 as a function of the  field radius. This radius is the only parameter in this simulation. The reference value for the stress intensity factor K1 is given by Wilson and Yu (1979) and can be written as: 0.4978 ≤

K1 ≤ 0.5144 √ (E ˛ T0 /(1 − )) a

(34)

Table 1 shows that the results of the thermal stress intensity factor K1 agree with the literature (Wilson and Yu, 1979), quite independently of the virtual extension field  in the formulation. The accuracy of the stress intensity computation decreases and can become unacceptable when the radius value becomes as great as the crack length. Nevertheless, when the radius tends to zero, the computed stress intensity factor definitely tends to be in the range of the numerical value obtained by Wilson and Yu (1979). One can Table 1 Normalized stress intensity factors calculated for thermal stress problem to be compared to the range [0.4978, 0.5144] Wilson and Yu, 1979. radius on a

K1 √ (E ˛ T0 /(1−)) a

0.02 0.05 0.1 0.2 0.3 0.4 0.5 0.6

0.5023 0.5069 0.4926 0.4985 0.4821 0.4633 0.4398 0.4512

2584

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588 Table 2 Difference in the formulation between Dolbow and Gosz (2002) and Menouillard et al. (2006). Fields

Fig. 5. Geometry and loading of the numerical experiment.

notice that a maximum relative error of 4% is ensured with a radius ratio of 30%, which defines a quite great contour size. 5.2. Example 2: illustrative mixed-mode case Here we study the influence of the temperature field near the crack tip through the energy release rate and stress intensity factors, with thermal and mechanical loadings. The temperature field is given as T (x, y) = T0 + T,x x + T,y y where the origin (x = 0, y = 0) is the crack tip as shown in Fig. 2. The example’s geometry is described in Fig. 5. The length L is 2 m. On the left and right side of the specimen, traction free is assumed: Fx = Fy = 0. On the bottom of the specimen, the displacement is imposed to be zero: ux = uy = 0. And on top of the specimen, a force with an angle of 25◦ with the vertical axis, is applied: Fx = F cos(65◦ ) and Fy = F sin(65◦ ), where F = 5 kN. The loading creates a mixed mode at the crack tip. The material model is isotropic linear elastic. The main characteristic is a linear dependence of the Young’s modulus with the temperature as follows: E = ETref −  (T − Tref ) = E0 −  T

(35)

where  is the material coefficient of this dependence (in Pa/m), Eref and Tref are reference Young’s modulus and temperature at a reference point. For this numerical example, the material parameters used are the following:  = 0.3 Pa/m, E0 = 70 GPa. Despite the space dependence of the Young modulus, the computation of the stress intensity is performed with the displacement auxiliary field u0i corresponding to the constant material parameters whose values are the one of the variable material taken at the crack tip position, i.e. E0 = E(xtip ). The behavior law is still depending on the space as C(x). Dolbow and Gosz (2002) use this auxiliary field u0i to have a computation as simple as possible. However Dolbow and Gosz (2002) and Menouillard et al. (2006) showed that the condition to keep having good results with the simple constant asymptotic auxiliary field, is to use a small radius for the virtual field in the formulation a(·, ·). Another method would be to use the actual asymptotic displacement field usi , but the computation then becomes slightly more costly, because we need to evaluate them. Thereby, Dolbow and Gosz (2002) showed that the fact that the material properties depend on the space does not affect the complexity of the computation. Table 2 presents the difference in the formulation between Dolbow and Gosz (2002) and Menouillard et al. (2006). ε(·) denotes the deformation gradient operator defined by 1/2(∇ · +∇ ·T ). The difference between the two references is

Dolbow and Gosz (2002)

Menouillard et al. (2006) Option 1

Option 2

Asymptotic displacement Asymptotic strain

u0 −1 C(x) · 0

u0 ε(u0 )

us ε(us )

Asymptotic stress

0

C(x) · ε(u0 )

C(x) · ε(us )

the fact that Menouillard et al. (2006) defines all fields (strain and stress) from the knowledge of the asymptotic displacement fields (u0i and usi ), so that the compatibility equations are satisfied, whereas they are not in Dolbow and Gosz (2002). However the equilibrium of the asymptotic fields are not ensured in Menouillard et al. (2006). Under these mechanical loadings and boundaries, and without any thermal effect, the energy release rate is evaluated √ 101.6 Pa m, and the stress intensity factors K1 = 15 MPa m and √ K2 = 1.5 MPa m. The direction of the crack propagation comes from the maximum hoop stress criterion: c = −10.6◦ . Under the same mechanical loading, different thermal loadings are applied. This example aims at studying the effects of the temperature field on the different fracture variables (i.e. energy release rate and stress intensity factors). As written previously, the temperature field is taken as T0 + T,x x + T,y y; so the study deals with the three parameters T0 , T,x and T,y . x indicates the crack direction, and y is perpendicular to the crack and oriented to the top of the structure. The origin (0, 0) is located at the crack tip position. 5.2.1. Influence of the temperature T0 A constant temperature field T0 applied in the specimen under the same mechanical boundary conditions. The influence of the temperature T0 on the energy release rate and stress intensity factors is studied. First, Fig. 6(a) shows the two stress intensity factors as a function of the temperature parameter T0 . It shows that the two values are constant and strictly independent of T0 . Second, as shown in Fig. 6(a), the crack direction is independent too. Third, Fig. 6(b) presents the evolution of the energy release rate as a function of the temperature parameter T0 , and it shows that the energy release rate significantly increases with the temperature parameter T0 because of the terms contained in the first Eq. (24). 5.2.2. Influence of temperature gradient T,x A temperature field T,x x is applied in the specimen under the same mechanical boundary conditions. The influence of the temperature gradient T,x in the temperature field is studied. Fig. 7(a) shows the two stress intensity factors as a function of the gradient temperature parameter T,x . It shows that the two values increase as the temperature gradient T,x increases. However, as shown in Fig. 7(a), the crack direction is almost independent of the temperature gradient, because the ratio between the stress intensity factors stays almost constant. Then, Fig. 7(b) presents the evolution of the energy release rate as a function of the temperature gradient T,x , and it shows that the energy release rate increases with the temperature gradient T,x . 5.2.3. Influence of temperature gradient T,y A temperature field T,y y is applied in the specimen under the same mechanical boundary conditions. Fig. 8(b) presents the evolution of the energy release rate as a function of the temperature gradient T,y , and it shows that the energy release rate does not vary so much with the temperature gradient T,y . Fig. 8(a) shows that the stress intensity factors slightly vary with the temperature gradient

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

Fig. 6. (a) Stress intensity factors and crack propagation direction, and (b) energy release rate as a function of the temperature field T (x, y) = T0 for the material model: E = E0 − T .

T,y . Nevertheless, Fig. 8(a) presents the direction c as a function the temperature gradient T,y , and a significant variation is shown because of K2 variation. The crack direction mainly depends on the temperature gradient T,y , which is the gradient perpendicular to the crack direction. 5.2.4. Conclusion The influence of the static temperature field in the domain can modify the fracture parameters. A positive temperature gradient in the direction of the crack encourages the crack to propagate by increasing the energy release rate, but does not affect the direction of propagation, neither the stress intensity factors. Furthermore, a temperature gradient perpendicular to the crack makes the crack turn, but does not significantly change the energy release rate. 5.3. Example 3: crack propagation in a heated strip 5.3.1. Experiment This example deals with a commonly used experiment (Marder, 1994; Bahr et al., 1995; Ronsin et al., 1995; Ronsin and Perrin, 1998; Yang and Ravi-Chandar, 2001). A long glass specimen whose width w and thickness e = 1.1 mm is moved at velocity V = 1 mm/s from an oven at temperature Thot ∈ [135; 180 ◦ C] to a cold water at ambient temperature Tcold = 20 ◦ C. The temperature difference T between the oven and the water bath is defined by T = Thot − Tcold . V is called the dipping velocity and is small enough to consider quasi-static thermal phenomenon. The thermal conductivity of the glass is 1.022 J/m/s/K. The distance between the oven and the water bath is denoted by h. Fig. 9 shows the geometry of the experiment and the four parameters of the experiment:

2585

Fig. 7. (a) Stress intensity factors and crack propagation direction, and (b) energy release rate as a function of the temperature field T (x, y) = T,x x for the material model: E = E0 − T .

h, w and T . Their values define the setting of the experiment, and are given in Fig. 10(a) and (b). The material properties of the glass are: Young’s modulus E = 75 GPa, Poisson’s ratio  = 0.23, the coefficient of thermal expansion ˛ = 7.7E −6 /K and the fracture √ toughness KIc = 0.33 MPa m. As a small dipping velocity is considered, the temperature field is linear between the water and the oven. The temperature field given by the boundary conditions in the plate creates a tensile stress at the region near the entrance of the oven (Ronsin and Perrin, 1998). This is the reason why the crack may propagates when its tip is located in this region. Thus, the crack tip remains in the region of the entrance of the oven. 5.3.2. Numerical results One observes that the crack propagates straight for the three simulations. Fig. 10(a) presents the vertical position of the crack tip as a function of time; one can notice that the crack tip moves down at velocity V = 1 mm/s when it is in the oven, then stays at the same vertical position during the remaining time of the simulation. Fig. 10(b) presents the mode 1 stress intensity factor as a function of time. This is not surprising that once the fracture toughness is reached, the relative position of the crack remains the same as the temperature field is constant. Thus a straight crack propagation occurs; the results of the experiment (Ronsin and Perrin, 1998) is shown in Fig. 11 whether or not the fracture toughness is reached. One can notice that the stress intensity factor remains constant during the quasi-static propagation of the crack, and its respective value is close to the fracture toughness. To conclude, the positive temperature gradient in the direction of the crack make it propagate.

2586

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

Fig. 8. (a) Stress intensity factors and crack propagation direction, and (b) energy release rate as a function of the temperature field T (x, y) = T,y y for the material model E = E0 − T .

Fig. 10. (a) Vertical position of the crack tip and (b) mode 1 stress intensity factor as a function of time for the straight crack propagation obtained with a dipping velocity of V = 1 mm/s.

Fig. 11. Experimental crack path obtained by Ronsin and Perrin (1998): (a) no propagation, and (b) straight propagation.

Fig. 9. Geometry and loading of the experiment.

6. Conclusions This paper presents a theoretical development and a formulation to compute the stress intensity factors for the particular static case with thermal and mechanical loadings and the restriction that the temperature field is linear in space near the crack tip vicinity. The proposed method is based on the general formulation of the energy release rate which leads to the construction of a bilinear form, a linear one and a scalar. These mathematical functions are then used to write two equations in order to capture the two stress intensity factors for the knowledge of the stress, displacement fields and the asymptotic displacement fields near crack tip.

Numerically, the results are first validated through the comparison with the literature. Methodology developed leads to small error if the -radius is sufficiently small. The second example aims at studying the effect of the temperature field on the energy release rate and the stress intensity factors. It shows the effect of the temperature field on the fracture parameters; e.g. a positive temperature gradient in the direction of the crack encourages the propagation, because this gradient makes the energy release rate increase. Moreover, a temperature gradient tends to make the crack turn in the opposite direction of the temperature gradient. Third, a numerical simulation of an existing experiment shows that a positive temperature gradient in the direction of the crack can make it propagate if the gradient is great enough.

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

Appendix A. Stress and displacement fields near crack tip

In order to solve the problem, one can notice that exponent n in Eq. (42) has to be 2 to be able to find out a solution. The particular displacement field uT , where the temperature is T (x, y) = T0 + T,x x + T,y y, becomes:

A.1. Asymptotic stress fields Mode 1:

⎡

 3 cos(/2) 1 − sin 2 sin 2 ⎣ g1 = √ 2 r sym Mode 2:



−1 g2 = √ 2 r

sin

 2

2 + cos

3  cos 2 2

 3 sin cos 2 2 3  1 + sin sin 2 2

cos

sym

 2

⎛

⎤ ⎦

(36)

(ex ,ey )



 3 sin 2 2  3  sin cos cos 2 2 2 1 − sin



u01 (r, )

1+ = E0



Mode 2: u02 (r, ) =

1+ E0



⎛

 ⎞

⎜ cos

⎟  ⎟  ⎠

r ( − cos )⎜ ⎝ 2

⎛ r ⎝ 2

sin

 2

2

(37)

(38) (er ,e )



 3 ( + cos ) + 2 sin 2 2 3  ( + cos ) − 2 cos cos 2 2

− sin

⎞ ⎠

(39) (er ,e )

where  is the Kolosov’s coefficient defined as a function of the Poisson’s ratio  by 3 − /(1 + ) in plane stress, and 3 − 4 in plane strain. Appendix B. Particular thermal displacement field uT The displacement field uT denoted (ur , u ) in the basis (er , e ) due to the thermal term ˛ T (x), can be extracted from Eqs. (12)–(15) as: ε(uT (x)) = ˛ T (x)Id

(40)

Writing out this equation leads to:

⎛ ⎜

˛T Id = ⎝

∂ur ∂r

1 2



sym

∂u 1 ∂ur u −  + r ∂ r ∂r 1 ∂u ur + r ∂ r

⎞ ⎟ ⎠

(41)

where ur , u define the displacement field uT in the polar basis at the crack tip vicinity (see Fig. 2). Here we make the assumption that the solution has the following dependence on the local polar coordinates near the crack tip: −−→ uT = r n f () = r n



f1 () f2 ()





= rn (er ,e )

F1 () F2 ()



(42) (ex ,ey )

The relation between the two basis (polar and cartesian), shown in Fig. 2, is written as: F1 () = f1 () cos() − f2 () sin()

(43)

F2 () = f1 () sin() + f2 () cos()

(44)

The temperature gradient in the polar coordinates near crack tip can be written as:



∇T =

T,x cos() + T,y sin() −T,x sin() + T,y cos()



(45) (er ,e )

uT = ⎝

⎞

˛

r2 (T,x cos() + T,y sin()) + ˛T0 r 2 ⎠ r2 ˛ (T,x sin() − T,y cos()) 2 (er ,e )

(46)

References

A.2. Asymptotic displacement fields: case of a material with constant material properties Mode 1:

2587

Attigui, M., Petit, C., 1997. Mixed-mode separation in dynamic fracture mechanics: new path independent integrals. International Journal of Fracture 84 (1), 19– 36. Bahr, H., Gerbatsch, A., Bahr, U., Weiss, H., 1995. Oscillatory instability in thermal cracking: a first-order phase-transition phenomenon. Physical Review E 52 (1), 240–243. Bui, H., 1978. Mecanique de la rupture fragile. Masson. Chapuliot, S., Lacire, M., Marie, S., Nédélec, M., 2005. Thermomechanical analysis of thermal shock fracture in the brittle/ductile transition zone. Part I: Description of tests. Engineering Fracture Mechanics 72 (5), 661–673. Dolbow, J., Gosz, M., 2002. On the computation of mixed-mode stress intensity factors in functionally graded materials. International Journal of Solids and Structures 39 (9), 2557–2574. Erdogan, F., Wu, B., 1996. Crack problems in FGM layers under thermal stresses. Journal of Thermal Stresses 19 (3), 237–265. Gu, P., Asaro, R., 1997. Cracks in functionally graded materials. International Journal of Solids and Structures 34 (1), 1–17. Hasselman, D., 2006. Unified theory of thermal shock fracture initiation and crack propagation in brittle ceramics. Journal of the American Ceramic Society 52 (11), 600–604. James, L., 1976. Fatigue-crack propagation in austenitic stainless steels (effects of nuclear reactor operating conditions). Atomic Energy Review 14, 37–86. Janssens, K., Niffenegger, M., Reichlin, K., 2009. A computational fatigue analysis of cyclic thermal shock in notched specimens. Nuclear Engineering and Design 239 (1), 36–44. Jin, Z., Paulino, G., 2001. Transient thermal stress analysis of an edge crack in a functionally graded material. International Journal of Fracture 107, 73– 98. Kim, J., Paulino, G., 2002. Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. International Journal for Numerical Methods in Engineering 53 (8), 1903–1935. Kim, J.-H., Paulino, H., 2003a. An accurate scheme for mixed-mode fracture analysis of functionally graded materials using the interaction integral and micromechanics models. International Journal for Numerical Methods in Engineering 58, 1457–1497. Kim, J.-H., Paulino, H., 2003b. Mixed-mode J-integral formulation and implementation using graded elements for fracture analysis of nonhomogeneous orthotropic materials. Mechanics of Materials 35 (1–2), 107–128. Lindau, R., Möslang, A., 1994. Fatigue tests on a ferritic-martensitic steel at 420 ◦ C: comparison between in-situ and postirradiation properties. Journal of Nuclear Materials 212, 599–603. Marder, M., 1994. Instability of a crack in a heated strip. Physical Review E 49 (1), 51–54. Markworth, A., Ramesh, K., Parks, W., 1995. Modelling studies applied to functionally graded materials. Journal of Materials Science 30 (9), 2183–2193. McGowan, J., 1979. Application of warm prestressing effects to fracture mechanics analyses of nuclear reactor vessels during severe thermal shock. Nuclear Engineering and Design 51 (3), 431–444. Menouillard, T., Elguedj, T., Combescure, A., 2006. Mixed mode stress intensity factors for graded materials. International Journal of Solids and Structures 43, 1946–1959. Noda, N., Jin, Z., 1993. Thermal stress intensity factors for a crack in a strip of a functionally gradient material. International Journal of Solids and Structures 30 (8), 1039–1056. Rao, B., Rahman, S., 2003. Mesh-free analysis of cracks in isotropic functionally graded materials. Engineering Fracture Mechanics 70, 1–27. Reytier, M., Chapuliot, S., Marie, S., Nédélec, M., 2006. Thermomechanical analysis of thermal shock fracture in the brittle/ductile transition zone—Part II: Numerical calculations and interpretation of the test results. Engineering Fracture Mechanics 73 (3), 283–295. Rice, J., 1968a. A path independent integral and the approximate analysis of strain concentration by cracks and notches. Journal of Applied Mechanics (Trans ASME) 35, 379–386. Rice, J., 1968b. Mathematical analysis in the mechanics of fracture. Fracture: An Advanced Treatise 2, 191–311. Ronsin, O., Perrin, B., 1998. Dynamics of quasistatic directional crack growth. Physical Review E 58 (6), 7878–7886. Ronsin, O., Heslot, F., Perrin, B., 1995. Experimental study of quasistatic brittle crack propagation. Physical Review Letters 75 (12), 2352–2355. Suo, X., Combescure, A., 1992a. Energy release rate and integral for any nonhomogeneous material. The American Society of Mechanical Engineers 233, 173–179.

2588

T. Menouillard / Nuclear Engineering and Design 240 (2010) 2579–2588

Suo, X., Combescure, A., 1992b. On the application of the G method and its comparison with de Lorenzi’s approach. Nuclear Engineering and Design 135, 207–224. Suresh, S., Mortensen, A., 1997. Functionally graded metals and metal-ceramic composites: Part 2. Thermomechanical behaviour. International Materials Reviews 42 (3), 85–116. Williams, M., 1957. On the stress at the base of a stationary crack. Journal of Applied Mechanics 24, 111–114.

Williams, J., Ewing, P., 1972. Fracture under complex stress—the angled crack problem. International Journal of Fracture Mechanics 8 (4), 441–446. Wilson, W., Yu, I., 1979. The use of the J-integral in thermal stress crack problems. International Journal of Fracture 15 (4), 377–387. Yang, B., Ravi-Chandar, K., 2001. Crack path instabilities in a quenched glass plate. Journal of the Mechanics and Physics of Solids 49 (1), 91–130.