An explicit integral expression for the stress intensity factor of a semi-elliptic surface crack subjected to thermal transient loading

International Journal of Pressure Vessels and Piping 76 (1999) 631–639

An explicit integral expression for the stress intensity factor of a semi-elliptic surface crack subjected to thermal transient loading Y.W. Kim*, J.I. Kim, M.H. Chang Korea Atomic Energy Research Institute, P.O. Box 105, Yu-Sung, Taejon, 305-600, South Korea Received 26 February 1999; accepted 14 March 1999

Abstract An explicit integral expression for the stress intensity factor of a semi-elliptic surface crack in a plate subjected to thermal transient loading was developed. The stress intensity factor of a semi-elliptic surface crack in a plate, which is exposed to a step change of fluid temperature, was calculated on the basis of the weight function method. The change of the stress intensity factor for a semi-elliptic surface crack subjected to an arbitrary change of the boundary fluid temperature was obtained by Duhamel integration for the product of the step function result and the time varying fluid temperature. The result obtained by the present method has shown good agreement with those obtained by the influence function method. As a practical application, a parametric analysis was performed for the crack behavior during the emergency cool down of reactor coolant in the reactor pressure vessel. Also, the present expression can be effectively applied to the simulation of fatigue crack growth of a semi-elliptic surface crack subjected to various thermal transient loading. q 1999 Elsevier Science Ltd. All rights reserved. Keywords: Stress intensity factor; Semi-elliptic surface crack; Thermal transient loading

1. Introduction The structural integrity of the pressure vessel is dominated by the pressure loading and the thermal transient loading. The load induced from these thermal transients has been recognized as one of the major causes of a mechanical failure. In the nuclear reactor, the emergency cool down of the coolant temperature induces large tensile stress on the inner wall of the pressure vessel and may cause unstable crack growth. Therefore, the monitoring of the semi-elliptic surface crack on the inner wall of the reactor is very important to the safety of the nuclear pressure vessel. As for another practical problem, the numerical simulation for the fatigue crack growth of a semi-elliptic surface crack subjected to thermal transient loading needs a huge amount of calculation for the stress intensity factor as a function of the crack length and the thermal transient loading. The stress intensity factor of this semi-elliptic surface crack has been calculated with continuous improvement to the methodology [1–4]. Several numerical procedures were presented for the evaluation of the stress intensity factor subjected to thermal shock loading [5–11]. Due to its complex geometry and transient nature, the stress intensity * Corresponding author. Tel.: 1 82-42-868-8981; fax: 1 82-42-8688990. E-mail address: [email protected] (Y.W. Kim)

factor of a semi-elliptic surface crack subjected to thermal transient loading is evaluated by the transient finite element analysis which needs laborious meshing and time-consuming analysis. The influence function method, which requires the polynomial fitting of the stress distribution at each time step, has been used as an alternative approach. A numerical computation procedure for the evaluation of the stress intensity factor of the two-dimensional (2D) through crack was presented based upon the Green’s function concept and the Duhamel’s integral and the results of the through crack have been presented [12]. The present study is directed to develop a simple explicit expression for the evaluation of the stress intensity factor of a semielliptic surface crack subjected to various thermal transient loading. Results obtained by the present method were compared with those by the influence function method. As for practical problem, the stress intensity factors were calculated for the surface crack in a reactor pressure vessel subjected to the emergency cool down of the coolant.

2. Method of approach Consider a semi-elliptic surface crack in a plate in which the crack size and the thickness is much smaller than the other dimension as shown in Fig. 1. One side of the plate is exposed to fluid of time-varying temperature and the other

0308-0161/99/$ - see front matter q 1999 Elsevier Science Ltd. All rights reserved. PII: S0308-016 1(99)00026-5

632

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

for weight function and loading KI ˆ

Za 0

syy …x†m…x†dx

…2†

where syy …x† is the crack face stress distribution and m…x† is the weight function related to the crack opening displacement. From the definition of the weight function and the stress field induced from the step change of the boundary fluid temperature, Green’s function for the stress intensity factor of a crack in a plate subjected to a change of boundary fluid temperature can be expressed as follows: GKI …t† ˆ Fig. 1. Semi-elliptic surface crack in a large plate.

side is insulated. There exists a semi-elliptic surface crack of depth a, and width 2c in the surface facing the fluid. It has been assumed that the thermo-mechanical properties do not change during the thermal transient analysis and the strain rates due to the thermal transient loading are much less than the strain rates caused by the stress wave propagation in the elastic body. Also, there is no convection through the crack face and the presence of the semi-elliptic crack does not have an effect on the heat conduction since the crack lies parallel to the heat conduction direction. When the boundary input loading varies continuously with time, the Duhamel’s theorem used in the structural dynamics, is quite useful for the extension of known solutions in terms of the corresponding case with time-independent boundary conditions. On the basis of this background, the stress intensity factor of the cracked body subjected to thermal transient loading can be evaluated by the integration of the product of boundary temperature and the fundamental solution obtained from a step change of boundary temperature. 2.1. Stress intensity factor after a step change of boundary temperature Imagine a step change of the boundary fluid temperature

u…t† ˆ u0 H…t†

…1†

where H…t† means the Heaviside step function; u represents the change of boundary fluid temperature and u 0 is the magnitude of the step function. Let us define the variation of mode I stress intensity factor due to a unit step change of the boundary fluid temperature as the Green function, GKI …t†. The stress intensity factor can be effectively evaluated utilizing the weight function method [13–15]. The characteristic of the weight function is independent of the loading, which enables the stress intensity factor to be calculated using a simple integral expression of a work-like product

Za 0

‰syy …x; t†ŠuˆH…t† m…x†dx

…3†

where ‰syy …x; t†ŠuˆH…t† is the stress distribution as a function of time after the step change of boundary temperature [16]. When a step temperature change was applied, the stress intensity factor of the semi-elliptic crack can be obtained from Eq. (3) and the weight function of the semi-elliptic surface crack " ∞ Za Ea X sin li A exp…2tl2i † GKI …t† ˆ 1 2 n iˆ1 li 1 sin li cos li 0 !! ! ( x x sin li 14 3 21 × 2 2cos li 1 2 li W W ! )# x li sin li 1 cos li 2 1 1 12 1 2 2 W l2i " ×

1 M2A

( !1=2 2 x p 1 1 M1A 1 2 a 2p…a 2 x†

x 12 a

! 1 M3A

x 12 a

!3=2 )# dx …4†

GBKI …t†

ˆ

Za

"

∞ Ea X sin li exp…2tl2i † 1 2 n iˆ1 li 1 sin li cos li 0 !! ! ( x x sin li 14 3 21 × 2 2cos li 1 2 li W W ! )# x li sin li 1 cos li 2 1 1 12 1 2 2 W l2i

" ×

( !1=2 ! 2 x x p 1 1 M1B 1M2B px a a

1 M3B

x a

!3=2 )# dx …5†

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

633

from Eqs. (4) and (9). " ∞ Zt Za Ea X sin li A KI …t† ˆ 1 2 n iˆ1 li 1 sin li cos li 0 0 !! ( x 2 × exp…2…t 2 t†li † × 2 2cos li 1 2 W ! x sin li 21 14 3 li W ! )# x li sin li 1 cos li 2 1 1 12 1 2 2 W l2i " Fig. 2. Calculation of stress intensity factor for time-dependent boundary temperature.

x 12 a

1 M3A

where

bˆ

×

( !1=2 ! 2 x x p 1 1 M1A 1 2 1M2A 1 2 a a 2p…a 2 x†

hW k

dx

du… t † dt : dt …10†

…6†

Superscripts A and B, respectively, stand for the deepest point …f ˆ p=2† and the surface point …f ˆ 0†. h, k , and a are convection coefficient, conductivity, and the thermal expansion coefficient, respectively. E is the Young’s modulus and n represents Poisson’s ratio. All of the material properties were assumed to be constant in order to maintain linearity of the governing equations. l i is the solution of the following characteristic equation.

li tan li ˆ b:

…7†

The weight function developed by Shen and Glinka [4] was used in the present analysis. 2.2. Stress intensity factor expression for arbitrary thermal transients

In a similar manner, the stress intensity factor at the surface …f ˆ 08 † can be expressed from Eqs. (5) and (9). " ∞ Zt Z a Ea X sin li B KI …t† ˆ exp…2…t 2 t†l2i † 1 2 n iˆ1 li 1 sin li cos li 0 0 !! ! ( x x sin li 14 3 21 × 2 2cos li 1 2 li W W ! )# x li sin li 1 cos li 2 1 1 12 1 2 2 W l2i " ×

( !1=2 ! 2 x x p 1 1 M1B 1M2B px a a x a

1 M3B

!3=2 )#

An arbitrary change of fluid temperature can be expressed as follows:

u…t† ˆ Ta …t† 2 T0

…8†

where T0 represents the initial temperature where t ˆ 0 and Ta is the boundary temperature expressed as a function of time for an arbitrary thermal transient. Incorporating the definition of Green’s function to Duhamel’s theorem, the stress intensity factor at time t can be represented as follows: KI …t† ˆ

!3=2 )#

Zt 0

GKI …t 2 t†

du… t † dt dt

…9†

The calculation procedure for the stress intensity factor is explained in a schematic diagram in Fig. 2. The stress intensity factor for a semi-elliptic crack in a plate subjected to arbitrary thermal transient can be derived

dx

du… t † dt : dt …11†

The stress intensity factor for the crack having other shapes, like a corner crack subjected to thermal transient can be derived by a similar procedure. When the Green’s function and the fluid temperature can be expressed as an explicit function the stress intensity factor during the thermal transient can be expressed as a simple integral function. However, the numerical integration should be used when the fluid temperature cannot be expressed as a simple explicit function. The right-hand side of Eq. (9) can be divided into two terms by introducing the concept of a decay time: KI …t† ˆ

Zt 2 td 0

GKI …t 2 t†

Zt du du dt 1 GKI …t 2 t† dt dt dt t 2 td …12†

634

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

Fig. 3. Stress intensity factor of semi-elliptic surface crack for step change of boundary temperature.

where td is the decay time. After the decay period, the Green’s function converges to a constant value for most elastic solids and it is no longer the function of time. For numerical integration, the decay time td is divided into n equal time steps. Therefore Eq. (12) can be rewritten as follows: KI …t† ˆ

n X

GKI …td 2 ti †

iˆ1

{u…ti † 2 u…ti 2 Dt†} Dt Dt

1 GKI …td †{u…t 2 td † 2 u…0†}

…13†

where

ti ˆ ti21 1 Dt n X

D ti ˆ t d

…14† …15†

iˆ1

Dt is a time increment for numerical integration and it should be small enough to minimize the numerical integration error. The numerical value of decay time td can be obtained from the analysis result of the step change of the boundary fluid temperature. When the convection coefficient and other material properties change with the temperature, the present calculation procedure cannot be used because of the violation of the linearity of the whole system. 3. Results of numerical calculation

the crack aspect ratio of 0.2 and 0.4. As non-dimensionalized time t p, which is defined below, approaches unity, the stress intensity factor becomes zero. tp ˆ

kt : rcW 2

…16†

The fluid temperature which varies according to the trigonometric function shape was used for a simple application of the present calculation to eliminate other complexities. The stress intensity factor for the simple trigonometric thermal transient has been calculated by the present method to be compared with those of the influence function method. In the influence function method, the stress intensity factor was evaluated from the product of the polynomial coefficient of the stress distribution and the influence function of Raju [1–2].  i x syy …x† ˆ Ai a iˆ0 3 X

…17†

where Ai is the coefficient calculated from the polynomial curve-fitting of the stress distribution at every time step. The stress intensity factor can be obtained from the product of the polynomial coefficients of the stress field and the influence function Fi. KI ˆ

r X pa 3 AF Q iˆ0 i i

…18†

where  1:65 a c

3.1. Stress intensity factor for linear change of boundary fluid temperature

Q ˆ 1 1 1:464

The stress intensity factor for the crack in a plate subjected to a step change of boundary temperature, i.e. Green’s function in this analysis, is shown in Fig. 3 for

Result of the present method has shown a good agreement with that of the influence function method as shown in Fig. 4(a) and (b). The present method is developed on the basis

…19†

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

635

Fig. 4. Stress intensity factor subject to linear change of boundary temperature. (a) a=c ˆ 0:2; (b) a=c ˆ 0:4.

of the superposition principle. When the applied loads tend to close the crack the singular behavior does not occur in the vicinity of the crack tip. That is, the stress intensity factor becomes zero. The restriction of the superposition method is discussed by Aamodt [17]. It should be noted that the stress intensity factor may be negative but it is meaningful only in

the case of superposition when the total stress intensity factors are positive. 3.2. Stress intensity factor after direct vessel injection As for practical problems, the thermal shock transient

636

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

Fig. 5. Stress intensity factor for exponential decay of coolant temperature. (a) f ˆ p=2; (b) f ˆ 0.

Fig. 6. Temperature change of reactor coolant after direct vessel injection.

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

637

Fig. 7. Change of stress intensity factor after direct vessel injection. (a) a=c ˆ 0:2; (b) a=c ˆ 0:4.

after the emergency injection of the coolant into the reactor vessel was solved. Since the other dimension are much larger than the thickness and the crack, the reactor vessel can be simply modeled as an infinite plate. The emergency cool down temperature of the reactor coolant is frequently modeled as a simple exponential function for the convenience of calculation in practical problem of the pressurized thermal shock problem.

u…t† ˆ uf 1 …u0 2 uf †e…2gt†

…20†

where g is determined from the shape of the coolant temperature distribution as a function of time. For this case, the time integration of Eq. (10) can be performed symbolically. " ∞ Za Ea‰2…u 2 u †Š X sin li g 0 f A KI …t† ˆ 2 l 1 sin li cos li li 2 g 12n 0 nˆ1 i

n

 e

2lt

2e

2l2i t

o

    x × 22cos li 1 2 W

  x sin li 21 14 3 li W   x li sin li 1 cos li 2 1 112 1 2 2 gŠ W l2i "

   2 x 1=2 × p 1 1 M1A 1 2 a 2p…a 2 x†      x x 3=2 1 M3A 1 2 dx: 1M2A 1 2 a a Results of the exponentially decaying transient calculated by using Eq. (21) are shown in Fig. 5(a) for the stress

638

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

Fig. 8. Change of stress intensity factor for various aspect ratio of semi-elliptic surface crack. f ˆ p=2; (b) f ˆ 0.

intensity factor of the deepest point …f ˆ p=2†. Similarly, the stress intensity factor at the surface …f ˆ 0† is derived from Eq. (11) and the numerical result is plotted in Fig. 5(b). The parametric study for this model can be easily performed utilizing Eq. (21) without time history finite element analysis. In general, an explicit function cannot describe the coolant temperature obtained by the direct vessel injection experiment [18]. As shown in Fig. 6, the thermal transient of the coolant calculated by the numerical simulation program of the direct vessel injection, also, cannot be expressed as a simple functional form. The transient solution is obtained from the numerical integration for the step function result and the transient boundary fluid temperature. The variation of the stress intensity factor for the crack tip …f ˆ 908† and surface …f ˆ 08† during this thermal transient is shown in Fig. 7(a) and (b), respectively. Similar results are shown in Fig. 8(a) and (b) for the semi-elliptic crack of a

different aspect ratio and compared with those of the influence function method. Results obtained from the present method have shown good agreement with those of the influence function method, which needs a laborious curve-fitting for the stress distribution along the cracked distance at each time step. For the elliptic crack having diverse aspect ratios, the variation of the stress intensity factor is shown in Fig. 8(a) and (b), respectively. At the deepest point the aspect ratio of the elliptic axes is dominant to the stress intensity factor but at the surface it is less.

4. Conclusion A simple explicit integral expression was developed for the evaluation of the behavior of the surface crack subjected to thermal transient loading. The stress intensity factor for the semi-elliptic surface crack was calculated by Duhamel’s

Y.W. Kim et al. / International Journal of Pressure Vessels and Piping 76 (1999) 631–639

integration for the product of step function results obtained from the weight function method and the time-varying fluid temperature. The results were compared with those obtained by the influence function method and have shown good agreement. The present calculation method is powerful to the calculation of the stress intensity factor for the semielliptic crack in the reactor vessel subjected to thermal transient. However, the temperature-dependent material property can not be considered in the present method since it was developed on the basis of linearity.

Acknowledgements This work has been performed under the nuclear research and development program sponsored by MOST.

References [1] Raju IS, Newman JC. Stress intensity factors for a wide range of semielliptical surface cracks in finite thickness plates. Engng Fract Mech 1979;11:817–829. [2] Newman JC, Raju IS. An empirical stress intensity factor equation for the surface crack. Engng Fract Mech 1981;15(1/2):185–192. [3] Shiratori M, Niyoshi T, Tanikawa K. Analysis of stress intensity factors for surface cracks subjected to arbitrarily distributed surface stress. In: Murakami Y, editor. Stress intensity factor handbook, 2. Oxford: Pergamon Press, 1987. [4] Shen G, Glinka G. Weight function for a surface semi-elliptic crack in a finite thickness plate. Theor Appl Fract Mech 1991;15:247–255.

639

[5] Xuejun F, Shouwen Y. Thermal shock fracture in a surface-cracked plate. Engng Fract Mech 1992;41(2):223–228. [6] Kim YW, Kim JH, Park JS, Kim JI. Evaluation on the stability of a crack in a reactor pressure vessel subjected to pressurized thermal shock. In: Proceeding of the 5th NPP Component Integrity Workshop by KSME, Korea, 1998, p. 473–78. [7] Kim JH, Kim YW, Kim TW, Park JS, Kim JI. Probabilistic analysis for the reactor vessel integrity subjected to pressurized thermal shock. In: Proceeding of the 5th NPP Component Integrity Workshop by KSME, Korea, 1998, p. 459–72. [8] Oliveira R, Wu XR. stress intensity factors for axial cracks in hollow cylinders subjected to thermal shock. Engng Fract Mech 1987;27(2):185–197. [9] Kim YW, Lee HY, Yoo B. Numerical evaluation of stress intensity factor for vessel and pipe subjected to thermal shock. Int J Press Ves Piping 1994;58:215–222. [10] Xuejun F, Shouwen Y. Thermal shock fracture in a surface-cracked plate. Engng Fract Mech 1992;41(2):223–228. [11] Bahr HA, Balke H, Kuna M, Liesk H. Fracture analysis of a single edge cracked strip under thermal shock. Theor Appl Fract Mech 1987;8:33–39. [12] Kim YW, Lee JH, Yoo B. An analysis of stress intensity factor for thermal transient problems based on Green’s function. Engng Fract Mech 1994;49(3):393–403. [13] Rice JR. Some remarks on elastic crack-tip stress fields. Int J Solids Struct 1972;8:751–758. [14] Bueckner HF. novel principle for the computation of stress intensity factors. Z Angew Math Mech 1970;50:129–146. [15] Xu XR, Carlsson AJ. Weight functions and stress intensity factor solutions, Oxford: Pergamon Press, 1991. [16] Boley BA, Weiner JH. Theory of thermal stress, New York, NY: Wiley, 1960. [17] Aamodt B, Bergan PG. On the principle of superposition for stress intensity factor. Engng Fract Mech 1976;8:437–440. [18] Cha JH, Jun HG. An investigation of fluid mixing with direct vessel injection. J Korean Nucl Soc 1994;26(1):63–77.