An analysis of stress intensity factor for thermal transient problems based on green's function

EngineerbigFtac&re &fechaaicsVol. 49, No. 3, pp. 393403, 1994

Pergamon

~1~~~~1~3

Elscvier Scienca Ltd Printedin Gnat Britain.

OOIE7944/94 n&l + 0.00

AN ANALYSIS OF STRESS I~ENSI~ FACTOR FOR THERMAL TRANSIENT PROBLEMS BASED ON GREEN’S FUNCTION YONG-WAN KIM,t JAE-HAN LEE and BONG YOO Mechanical Structure Development Department, Korea Atomic Energy Research Institute, P.O. Box 7, Daeduk, Taejon, 305-606, Korea Abetract-An

ethcient numerical method is presented to compute the stress intensity factor for the cracked body subjected to a thermal transient loading. It was developed on the basis of the Green’s function concept and Duhamel’s theorem. The stress intensity factor due to the thermal transient loading was calculated by an integration over the elapsed time interval for the product of the boundary temperature and the &en’s function. The Green’s function of the present method represents the change in stress intensity factor for the cracked body subjected to a Heaviside step function of the boundary ornate, and it can be determined by any method proven so far. The cracked wall and cracked pipe subjected to thermal transient loading were analyzed, for practical applications. The stress intensity factors obtained from the present method have showngood agnxmentwiththose evaluated by the standard finite element analysis. Compared to the transient finite element method, the present method is known to be more effective in view of computational aspect without sacrificing the solution accuracy.

1. ~RODU~ON ILMDS induced from thermal transients have long been recognized as one of the major causes of the mechanical failure. The stress intensity factors due to such loads are the function not only of the instantaneous boundary temperature loading but also of the past history of the boundary t~~ratu~, The stress intensity factor data of the thermal transient problem have been reported only for several typical cases. Although the general purpose numerical techniques such as finite element method can he applied to the analysis of more general cases, they require laborious time-consuming analysis. Finite element method has been used for the analysis of the cracked body subjected to a thermal loading. Wilson [1], Bahr [2] and Hellen [3] have used the tinite element method to compute the stress intensity factor for the thermal loading problem. For the thermal transient problem, however, the finite element analysis requires a step-by-step computation for the entire time range due to the transient nature of the problem. The Green’s function method has been applied to the crack tip loading problem for the computation of the stress intensisty factor [4,5]. The crack-face Green’s function obtained from the solution of point forces acting on the crack faces was then used to obtain the stress intensity factors for any problem of the same geometry under an arbitrary loading condition. The weight fiction technique first introduced by Bueckner 161has been used for the determination of the stress intensity factor. He has shown that the stress intensity factor due to an arbitrary set of applied loads can be obtained by integrating the product of these loads with the weight function of the cracked body over crack length. These two methods are especially suited when, for a given geometry, a large number of the stress intensity factor solutions for complex loading are desired. Also, it has been proved that these methods have remarkable ~mpu~tional efficiency. Wu [7,8], Xuejun [9] and Tsai [lo] employed the weight function method for the computation of the stress intensity factor in the thermal shock problem, while Shiratori [l 1J used the influence function method which has the basis of the superposition principle. These studies have been confined to the analysis of the simple thermal shock or time independent problem. In order to analyze the thermal transient problem by the use of the above methods, time-consuming analysis for the evaluation of temperature and/or stress ~st~bution for the untracked body is a task yet to be worked out. Therefore, more efficient methodology is necessary to compute the stress intensity factor in the practical thermal transient problem, tAuthor

to whom correspondence

should be addressed. 393

394

YONG-WAN KIM ef al.

In this study, a new numerical technique for the computation of the stress intensity factor was developed based upon the Green’s function concept and the Duhamel’s integral of the structural dynamics. It was derived for the two-dimensional (ZD) problem of homogeneous isotropic body for simplicity, and it can be extended directly to more practical problems such as the computation of stress intensity factor for three-dimensional (3D) cracked body, In contrast to the conventional use of the Green’s function as the response for a pair of point forces acting at the crack face, the Green’s function represents the time history solution of the cracked body subjected to the unit step change of the boundary temperature in the present method. For the pipe and the wall subjected to thermal transient loading the stress intensity factor was computed as a function from the Duhamel’s integration for the boundary temperature and the Green’s function. The restriction of the present method was discussed and also, the result of the present method was compared with that of the standard finite element analysis in order to investigate the computational efficiency and accuracy. 2. METHOD OF ANALYSIS 2.1. Green function approach The Green’s function is usually defined as the response of a system to a standard step or impulse input. The important property of the Green’s function is that, when suitably defined, it contains all the essential information of the system. There, it can be utilized to compute the response of the system to any input by considering it as being composed of a large number of small impulses. For the valid representation of the Green’s function, the system must satisfy three properties which are causality, invariance and linearity. In the linear elastic fracture mechanics, the Green’s function concept has already been used to obtain the stress intensity factor for cracks subjected to boundary pressure acting on the crack faces from the solution of the same crack subjected to point forces acting on the crack faces. It has been discussed that the Green’s function method is valid for the uncoupled thermoelastic problem and the stress distribution is the unique function of the temperature distribution under the quasi-static thermoelastic theory [12]. It has been assumed in the present analysis that the thermo-mechanical properties do not change during the thermal transient analysis and the strain rates due to thermal transient loading are much less than the strain rates caused by stress wave propagation in the elastic body. When the boundary input loadings vary continuously with time, the Duhamel’s theorem used in the structural dynamics is quite useful for the extension of known solutions in terms of that for the corresponding case with time-independent boundary conditions. Based upon these backgrounds, the stress intensity factor of the cracked body subjected to thermal transient loading can be evaluated from the intention for the product of boundary temperature and the Green’s function as explained in Fig. 1. The Green’s function means the fundamental solution obtained from a step change of the boundary temperature. Only the stress intensity factor for mode I will be considered in the present study and it is formally defined as follows

where au,, is the stress acting at distance r along the line ahead of crack. The input temperature at the boundary for the completion of the Green’s function is expressed as follows W)

= @oW),

(2)

where H(t) means the Heaviside step function. 8 represents the change of boundary temperature and Q. is the magnitude of the step function. From the definition of the stress intensity factor, the Green’s function of the present method is expressed as follows G,,(t) = lii’?%~,.~

(t),

395

SIF for thermal transient loading

‘8(t)-

#@-

tfme

Boundory

temperoture

Duhamel integration

--c 4

I

time

8

Stress intensity factor during thermal tmnsients

k

tlma

Step

time

function input

Fundamental

solution

Fig. 1. Procedure for determining the stress intensity factor of the thermal transient problem.

where aYv(t) means the stress induced from the step change of boundary temperature and G, stands for the Green’s function of mode I. The Green’s function needs to be computed once only for a given set of boundary conditions. 2.2. Stress intensity factor during thermal transients Consider a cracked body in which the boundary temperature changes as below G(t) = T,(t) - T,, where T,, represents the function of time for an On the basis of the stress intensity factor at be expressed as follows

(4)

initial temperature and T, is the boundary temperature expressed as a arbitrary thermal transient. Green’s function concept and the Duhamel’s theorem, the change of the time t due to a small change of the temperature boundary at time T can AK, (t, 2) = G, (t - z)AO (r) =

AQ@)

G,(t - ~1~

AT

,

where AK, means the change of mode I stress intensity factor induced from the change of boundary temperature A8 during the time segment AZ. The Duhamel’s integration is schematically illustrated in Fig. 2 where the stress intensity factor at time t can be obtained from the summation of all the individual responses acting separately. K,(t) = lim i AK,(t, q) . n-m i=, By replacing AK, with eq. (5), eq. (6) becomes K,(t)=lim iGKI(t-rJTAr. n-03,-I

AQ (ri)

Equation (7) can be written as follows K,(t)=

‘G,,(t-r)Tdr. s

EFM 490-F

0

d@(r)

(8)

YONG-WAN KIM et al.

396

t-td

0

7

t

Transient temperature

t-r

0

time

at boundary

time

Green’s function

Fig. 2. Duhamel’s integration for Green’s function and boundary temperature.

For the convenience of the numerical manipulation, separated into two parts as below:

s

I- t&q

I

G,U I- ,d

wo=

the integration

-

7)

dr +

y

G&

range of eq. (8) is

d@(~)dz --T+-

(9)

9

s0

where t,, is the decay period. After a decay period, the Green’s, function GK,(f) converges to a constant value for most of the elastic solids and it is not the function of time any more. In the right hand side of eq. (9), the first term means the integration during the decay period and the second term means the integration when the Green’s function is expressed as a constant value. The gomputation time can be significantly reduced due to the presence of decay time since the second inBgration can be done explicitly. By incorporating the concept of decay time into eq. (9), it can bo simply expressed as follows K,(t)=

s’ t-Q

G,(t - 2)q

dt + G&J{@(t

- td) - S(0)) .

The integration over the time range td is only necessary for the calculation of the stress intensity factor. The integration range td is divided into n steps for the numerical integration. 4(t)=

i

G,,(t~-ri)Q(Ti)-~~-AT)

AT +G,,(t,){Q(t

-4-@(o))

9

(11)

i-l

where 7iEti-l

+AZ

(12)

Finally, eq. (11) is expressed as below: 4(t)=

i G~~(f-~i){~(Ti)-~(~i-~)}+G~~(t~){Q(t-td)-~(O)).

(14)

i-l

The Green’s function and the Duhamel’s theorem can be applied to the calculation of the stress intensity factor for mode I without an explicit calculation for the distribution of the temperature and stress. Also, these procedures can be applied to the computation of stress intensity factor for other modes of fracture and 3D problems. 3. NUMERICAL

ANALYSIS

As a numerical example of the present method, the cracked pipe and the cracked wall are analyzed. It is assumed for the computation of Green’s function that stress analysis and temperature analysis are not coupled as specified in the previous section. In order to show that

SIF for thermal transient loading

397

the Green’s function can be calculated by any proven technique, the finite element method and the Bueckner’s weight function method are used in the calculation of the Green’s function for the analysis of the cracked pipe and the cracked wall, respectively. For the methodological comparison, the full transient finite element analysis was done with the commercial finite element code ANSYS by applying the temperature change directly at the boundary of the finite element model and performing a standard thermal and linear elastic fracture analysis. The results of the present method are compared with those obtained by the standard finite element analysis. 3.1. Pipe subjected to thermal transients The stress intensity factor during thermal transient loading is computed for a pipe with two axial cracks of the crack length to pipe thickness ratio 0.5. The inner wall of the pipe contacts with the fluid of temperature T, and the outer wall of the pipe is insulated. It is assumed that the presence of crack does not have an effect on the heat flow. The Green’s function is computed by the finite element analysis where the unit step change of fluid temperature is imposed on the inner wall of the pipe. As shown in Fig. 3, the quarter of the pipe is modeled in the finite element analysis. The finite element analysis includes the heat conduction analysis and the linear elastic fracture analysis. The temperature distribution for each time step is obtained from the heat conduction analysis and the nodal temperature is then used for the input load of the static finite element analysis. The stress intensity factors, which are the Green’s’ function of the present analysis, are computed from the results of static analysis by employing the quarter point element introduced by Barsoum [13]. The Green’s function computed from finite element analysis is plotted in Fig. 4 with nondimensionalized parameters. (151

(17) where p, c, u and h stand for density, specific heat capacity, heat conductivity and the convection coefficient, respectively. The Biot number j3 represents the thermo-physical property of the pipe. The thickness of the pipe and the crack length are expressed as Wand a, respectively. The Green’s function decreases to zero after the certain decay period. As the Biot number increases, the maximum stress intensity factor has increased. By using the Green’s function calculated, the stress intensity factor is computed for several typical thermal transients which frequently occur in the real piping system.

crack

flnite element

model

Fig. 3. pipe subjected to thermal transient loading and finite element model.

398

YONG-WAN KIM et al.

I

a/w = 0.5

0.6-

_o

00’ . 10-J -

L

A loo

.

lo-’

lo-’ t*

Fig. 4. Green’s function for a cracked pipe calculated by finite element analysis,

For the case when the fluid temperature

decreases as below

s(t) = -6&t*

(18)

the stress intensity factors calculated by the present method and standard finite element analysis are compared. As shown in Table 1, the difference of the stress intensity factor obtained by these two methods is less than 1%. As for the other two cases of more complex thermal transient shown in Fig. 5, in which the Biot number is 5, the stress intensity factors computed by the present method have shown good agreement with those by the transient finite element analysis. However, the present method has required a negligible computation time comparing with the standard finite element analysis. For the evaluation of the stress intensity factor during thermal transients in the present analysis, the computation of the Green’s function and the Duhamel’s integration for each time step are necessary. Comparing to the transient finite element analysis, the finite element analysis for each time step is replaced by the simple time integration which consumes much less computational time and efforts. As the number of thermal transient events increases, the present approach becomes more cost effective. 3.2. Wall subjected to thermal transients A cracked wall shown in Fig. 6 is analyzed by the present method. One side of the wall is exposed to the fluid and the other side is insulated. For this problem, the temperature and stress distribution of untracked body can be expressed explicitly. The Green’s function of the cracked wall is computed by the Bueckner’s weight function method which has been used in the analysis of the thermal shock problem [2, 141. The Green’s functions are computed from the integration over the crack length for the product of the stress of Table 1. Comparison of stress intensity factors computed from the present method and the standard finite eleqent method (B = 5) Time (r*) 0.1 0.3 0.5 0.7 0.9

Standard FEM (4/&J 0.02160 0.08591 0.13775 0.17690 0.20635

where t*=-

lcr

pew*’

Resent method (K&J 0.02175 0.08638 0.13859 0.17791 0.20731

40= -J;;;;. 1 _-v

Difference (%) 0.69 0.55 0.61 0.57 0.46

SIF for thermal transient loading

399

0.6 -

present

method

-

0.8

0.4

0.6

0.2

0.4

0.0

-0.2

0.2

0.0 -0.0

present method

-0.4 0.2

0.4

0.6

0.8

1.0

I 0.0

0.2

0.4

t*

0.6

0.8

1.0

t*

(a)

09

Fig. 5. Comparison of stress intensity factor obtained from Green’s function method and those computed from transient finite element analysis. (a) Case 1 thermal transient loading (J = 5); (b) case 2 thermal transient loading (B = 5).

untracked body subjected to a step change of the boundary temperature and the Bueckner’s weight function. It is expressed as follows

1

m(a*, X) dX,

where A, tan An= /I

(19)

(20)

(21)

fluid temperature t 0)

insulated

Fig. 6. Cracked wall subjected to thermal transients.

400

YONG-WAN KIM er al.

The length of the crack and the thickness of the wall are expressed as a and W, respectively. The A,,is the solution of the characteristic eq. (20). The Bueckner’s weight function is given as below m(a*, X) =

&4+m,~~~+m,(q5fy].

(23)

where m, = 0.6147 + 17.1844@*)‘+ 8.7822(~*)~

(24)

m, = 0.2502 + 3.2899@ *)’ + 70.0444(u *)’ .

(25)

The solution of eq. (18), that is the Green’s function of this problem, is plotted in Fig. 7 for cases where Biot numbers are 1, 5, 10, 20, 100, co and the crack length to wall thickness ratio is 0.1. Unlike the previous example of the cracked pipe in the previous section, the Green’s function is calculated without a finite element analysis. From eqs (8) and (19), the stress intensity factor for an arbitrary thermal transient loading is expressed as follows

-2cos&(l

X 1

X

-X)+4(3X-

sin 3, l)Y+ n

12(1-2X)

A, sin A, + cos 12,- 1 no }m(u*, X)dX]Fdr.

(26)

The stress intensity factor computed by the present method is compared with that evaluated by the transient finite element analysis in Fig. 8. The results obtained from two completely different approaches have shown good consistency. It can be deduced that any proved numerical method can be used for the determination of the Green’s function in the present analysis. By the use of the present method, the response of the stress intensity factor for various Biot numbers and the crack length are plotted in Figs 9 and 10, respectively. The full history of the stress intensity factor during thermal transient loading is efficiently obtained from the present method. This method can be usefully applied to the computation of the stress intensity factor for various design transients of the real structure such as pressure vessel and pipes.

Fig. 7.

function for a cracked wall calculated by Bueckner’s weight function method.

401

SIF for thermal transient loading

0.2

-0.2

-0.4 1.0

0.6

0.6

0.4

0.2

0.0

t* Fig. 8. Comparison of stress intensity factors by Green’s function method and those by transient finite element analysis for cracked wall (~9= 10).

3.3. Crack closure during thermal transients It has been specified in the derivation of the present method that the present method is developed on the basis of superposition principle. The restriction of the superposition method is discussed by Aamodt [IS]. When the applied load tends 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. 1.0 a3 & 8

0.0 -1.0 0.4

0.2 0 q

0.0

\ -0.4

0.6

0.6

1.0

t* Fig. 9. Time history of stress intensity factors during thermal transient for various Biot numbers.

1

0.0

a*= a/W

I

I

I

I

I

0.2

0.4

0.6

0.6

1.0

t* Fig. 10. Time history of stress intensity factors during thermal transient for various crack lengths (B = IO).

YONG-WAN KIM er al.

402

t

0

present method FEM(interface element)

-0.2

-0.4

0.0

I

I

I

0.2

0.4

0.6

a/W=O.l

I

0.8

1.0

t* Fig. 11. Comparison with result of finite element analysis utilizing interface elements for cracked wall (B = 10).

Because of the presence of the crack closure, the crack surfaces in the final configuration must always be separated along their entire length for the valid application of the superposition principle. Also, it should be noted that the stress intensity factor may be negative but it is meaningful only in the case of superposition where the total stress intensity factors are positive. In order to model the crack closure in the standard finite element analysis, the interface element is used at the crack face. The interface element has zero stiffness when the tensile loading is applied to the interface element but it has large stiffness when the compressive loading is applied to the interface element. Finite element analysis using the interface element is performed for the cracked wall subjected to cyclic thermal transient loading. As shown in Fig. 11, the results of the present method and the finite element method with interface element have shown good agreement when the total load tends to open the crack. However, results of the finite element analysis become zero when the total load tends to close crack. Therefore, the present method can be applied to the open crack problem. This is true for other methods which have the basis of superposition. However, for the application of the present method to the practical engineering problem, it is not so severe a problem since the negative stress intensity factor can be replaced by a zero value of stress intensity factor and the maximum value of the stress intensity factor is generally important. 4. CONCLUSIONS A numerical technique is presented for the efficient computation of stress intensity factor for the cracked body subjected to the thermal transient loading. The present method was developed on the basis of the Green’s function and the Duhamel’s theorem. For the verification of the present method, numerical analysis for the cracked wall and the cracked pipe subjected to thermal transients was carried out. The stress intensity factor computed by the present method has shown good agreement with those obtained from the standard transient finite element analysis of thermal transient problem. Compared to the standard finite element analysis, the present method requires much less computational time and effort since the transient finite element analysis for each time step can be replaced by the simple numerical integration. The present method is more effective as the thermal transient events increase. On the other hand, it has limitations similar to that of the superposition principle.

SIF for thermal transient loading

403

REFERENCES [l] W. K. Wilson and I.-W. Yu, The use of the J-integral in thermal stress crack problems. Int. J. Fracture 15, 377-387 [2] g??Bahr, H. Balke, M. Kuna and H. Liesk, Fracture analysis of a single edge cracked strip under thermal shock. Theor. appl. Fracture Mech. 8, 33-39 (1987). [3] T. K. Hellen, F. Cesari and A. Maitan, The application of fracture mechanics in thermally stressed structures. Int. J. Press. Vess. Piping 10, 181-204 (1982). [4] M. H. Aliabadi and D. P. Rooke, in Numerical Fracture Mechanics, pp. 193-202. Computational Mechanics Publications, Boston (1991). [S] V. Shivakumar and R. G. Forman, Green’s function for a crack emanating from a circular hole in an infinite sheet. Inr. J. Fracture 16, 305-316 (1980). [6] H. F. Bueckner, A novel principle for the computation of stress intensity factors. Z. angew. Murk. Mech. So, 529-546 t197m

171 R:&eira

and X. R. Wu, Stress intensity factors for axial cracks in hollow cylinders subjected to thermal shock. Engng

Fracture Med. 27, 185-197 (1987). X. R. Wu and A. J. Carlsson, Weighf Functions and Stress Intensity Factor Solutions. Per8amon Press, Oxford (1991). 6; F. Xuejun and Y. Shouwen, Thermal shock fracture in a surface cracked plate. Engng Fructure Mech. 41, 223-228

m

(1992). C. H. Tsai and C. C. Ma, Thermal weight function of cracked bodies subject to thermal loading. Engag Fracture Med. 41, 27-40 (1992).

VII M. Shiratori, Analysis of stress intensity factor for surface cracks in pipes by an influence function method. Proc. ASME ‘89 Conf. on P. V. Tech. (Pp. 45-50), Toyko (1989). B. A. Boley and J. H. Weiner, Theory of Thermal Stress. John Wiley, New York (1960). t::; R. S. Barsoum, Gn the use of isoparametric finite elements in linear fracture mechanics. Int. J. numer. Meth. Engng 10, 25-37 (1976). 1141K. Y. Leeand K. B. Sim, Thermal shock stress intensity factor by Bueckner’s weight function methods. Engng Fracture Mech. 37, 799-804 (1990). VI B. Aamodt and P. G. Bergan, On the principle of superposition for stress intensity factor. Engng Fracture Mech. 8, 437-440

(1976).

(Received 26 August 1993)