Stress intensity factors of cracks in finite plates subjected to thermal loads

&&e&ng fiaeture hiec&&c~ Vol. 43, No. 4, pp. 641-650, 1992 Printed in Great Britain.

oats-?944/92 g5.00+ 0.00 0 1992Pergamon Press Ltd.

STRESS INTENSITY FACTORS OF CRACKS IN FINITE PLATES SUBJECTED TO THERMAL LOADS K. Y. LAM, T. E. TAY and W. G. YUAN Department of Me&an&al and Production Enp;ineering, National University of Singapore, IO Kent Ridge Crescent, Singapore 0511, Republic of Singapore Abstract-In the framework of cracked plate tbermoelastic problems, the perturbation effects caused by the presence of a crack on thermal stresses, displacements and stress intensity factors in an isotropic linear elastic medium with varying crack surface heat conductivity under uniform heat flow are studied. The disturbance due to the crack results in a non-smooth temperature modihcation of the original field and induces singular character of the thermal stresses near the crack tip. This ~sturbing ante distribution is solved numericafjy by using the distributed dipole me&d, Finite clement solutions under the above temperature distribution are presented. Numerical results of mixed mode stress intensity factors are given for various geometric cases of a plate. It is shown that, for the case of a crack situated in the finite plate center subjected to a uniform heat flow, the geometry does not a&% the mode I thermal stress intensity factor Ki significantly, and Kg is negligible compared with the mode II thermal stress intensity factor Ka. Only the mode II thermal stress intensity factor Ka is significaut and shows a very Iwxdized phenomenon in the nea&eid region of a crack. For all cases studied here, the vaiues of the mode II thermal stress intensity factor in a smaller plate are dependent on the geometric size of the plate. As the geometric size of the plate increases, the value of Kl, is close to the solution of an i&&e plate.

I. INTRODUCTION Trm STUDY of th~~~lasti~ fracture mechanics, which deals with the events of catastrophic spreading of existing cracks subjected to thermal loading, is considered to be of great importance in the design of structures such as aerospace components, turbines and nuclear pressure vessels. A steady heat flow disturbed by the presence of a crack gives rise to a local intensification of the temperature gradient, resulting in large thermal stresses in the neighborhood of the crack tips. Sih and Liebowitz [l] have shown that singular behavior can be found in both heat flow and stress field. Sib [2] considered the singularities of two~~ensional thermal stresses at the crack tips in an infinite medium and showed that the I/.,/% stress sin~l~ty is preserved in the thermal stress problem. Static and transient thermal stress intensity factors have dso been discussed by many investigators [3-S]. Conventionally, crack surfaces are assumed to be insulated [2,6] when solving for thermal stress intensity factors. Such an assumption often leads to the question of how much the stress intensity factors will change when a partial thermal co~du~ti~ty is used for crack surfaces. Barber [7,8] obtained a steady-state thermal stress solution to the problem of a circular crack in an elastic solid with finite crack surface heat conductivity, Kuo [9] also treated finite heat conductivity at crack surfaces by imposing a set of distributed dipoles at the crack surfaces. He derived more general formulas for mixed mode stress intensity factors with varying crack surface heat conducti~ty by using the method of an influence function. In this work, the disturbing temperature field for the central cracked plate subjected to a uniform temperature gradient is solved n~e~~lly, and then the finite element method is used to solve for the thermal stress and displacement fields using the above disturbing temperature dist~b~tion as the thermal loads. Several examples with various geometries of the plates are given for calculating the stresses and displacements using MSCINASTRAN. Once a finite element solution has been obtained, the values of stress intensity factors must be extracted from it by using the relationship between the stress intensity factor and ~spla~ments in the crack region. The results of mixed mode thermal stress intensity factors are compared with the analytical solutions obtained in refs [2,9] to show the effects of the geometries of a plate on the mixed mode thermal stress intensity factors. 2. THEORY As illustrated in Fig. 1, B homogeneous isotropic linear elastic body t/ with a boundary surface S and a crack surface S, with no body force or distributed heat source is considered, It is assumed 641

K. Y. LAM et al.

642

Fig. 1. A crack in a linear elastic solid and iocal coordinate at the crack border.

that the crack surface SC lies in a plane. It is also assumed that the quasi-static thermoelasticity theory is valid, The thermoelastic stress-strain relations for homogeneous and isotropic media are given by: V

V

-&&-l-2v

1-2v

aT6,

(1)

where bti is the stress tensor, Q is the strain tensor, G is shear modulus, v is Poisson’s ratio, a is the thermal expansion coefficient and 6, is the Kronecker symbol. The differential equations of equilibrium with no body forces are: ir#,j = 0.

The str~n~spla~ent

(2)

relations are: (3)

&ij=(U,j+Uj,i)/2*

The mechanical boundary conditions are: triinj = c ui=up crvnj= 0

P on S,

(4)

Pon

(51

S,

P on SC

(6)

where fl is the prescribed traction on the traction boundary S,, uy is the prescribed displacement on the displacement boundary S,,, SCis the crack surface boundary and 5 are the unit vectors in the directions shown in Fig. 1. Fourier’s law of heat conduction is given by: 5 V2T = T,, where T is temperature,

(7)

t is time and V* is the Laplacian operator.

We define

where K is heat conductivity, p is mass density and c is specific heat. The thermal initial and boundary conditions including crack surfaces are: T(P, t) = F(P, g

t)

P on S,

(9)

(P, t) = q”(P, t)

P on S,

(10)

where To(P, t) is the prescribed temperature prescribed heat flux on the boundary S,.

on the temperature

boundary

Sr and q” is the

Stress intensity factors of cracks in finite plates

The crack surface thermal conductivity

conduction

643

equation is given by

f$(P,t)=i[T(P+,t)-T(P-,t)] Pon SC

(11)

where 1 is the dimensionless Biot number [7’J,L is the characteristic length of the problem [lo] and P+, P- represent the points at the upper and lower crack surfaces respectively. In eq. (1 l), the crack surfaces are insulated when 1 is zero and are perfectly heat conductive when R is infinity. Since this is a linear problem, we can assume that: t)

iW

o,(P , t) = d!‘(P rl ’ t) Y ’ t) + u’T’(P

(13)

T(P, t) = T”‘(P, t) + TyP,

where functions with superscript (1) are solutions to the problem for Z = co and functions with superscript (2) are solutions which account for the proper crack thermal conductivity. In the case of perfect heat conductivity of the crack surface, the temperature is continuous across the surfaces. T(i), 08) and ui’) satisfy eqs (l)-(10). The crack surface thermal conductivity condition, eq. (1 l), is replaced by FyP+,

t) = T@‘(P-, t).

(14)

In order to satisfy the actual crack surface thermal conductivity condition, eq. (1 l), the solution for imperfect heat conductivity of the crack surface, T’*), of) and UP) must be added to the solution of perfect crack surface thermal conductivity. The solution for imperfect heat conductivity of the crack surface has to satisfy eqs (l)-(3), (7), (8) and the following boundary conditions: cr@n.= rl 3 0 #(P

?t) = 0

T(*‘(P f t)=O

aF2)

an(P,t)=O

‘[““alf ‘*”(Q, t) = i

P on S D and S E

(19

P on S ”

(16)

P on ST

(17)

Pon

(18)

S,

[T(*)(Q+, t) - T”)(Q-, t)]

Q on SC,

(1%

From eq. (19) it is seen that, for finite R, the temperature across the crack surfaces is not continuous except at the crack tips; there is a temperature jump from one side of the crack to the other. Such a temperature jump across the crack surface resembles the temperature jump caused by distribution dipoles situated at a surface in a solid [7]. Therefore, eq. (19) can be satisfied if we can find a set of dipoles properly distributed at the crack surface SC, It has been shown by Carslaw and Jaeger [lo] that, with a set of dis~but~ dipoles of strength M(Q, t) at crack surface SC, the temperature distribution in Y is given by T”(P,

t) =

’ TG(P, Q, t - z)M(Q, t - r) dz dS(Q) IsSC 0

(20)

where T’(P, Q, T) is the temperature solution at point P due to an instan~n~us dipole of unit strength applied at point Q at t = 0. In other words, TG has to satisfy eqs (7), (17) and (18). With a method similar to that given by ref. [lo], it can be shown that f12) given by eq. (20) will create a temperature jump of it4(Q, T)/pc{ at point Q on the crack surface SC. Therefore, eq. (19) becomes:

* ss SC 0

WQt 0

aTG(P, Q, t - r) an(Q)

dr

dS(Q) + &Qt 0 = &

M(Qt f>

(21)

644

K. Y. LAM et al.

where

q(Q, t) =

aF’)(Q, t)

an

(22)

is the heat flux across the crack surfaces for the case of perfect heat conductivity of the crack surface. Equation (21) is a Fredholm integral equation. It is often difficult to obtain TG and T”) analytically except for some special cases.

3. EXAMPLE

PROBLEM

Let us consider a finite crack of length 2a in an infinite plate subjected to a steady-state temperature gradient of R at infinity. As shown in Fig. 2, the thickness of the plate is assumed to be D and the temperature gradient n is assumed to have an angle of inclination of /I with the x-axis. The plate is assumed to be stress-free at the boundary. The solution consists of two parts. For the case of Iz = co, T(l) is a linear function of x and y. According to Boley and Weiner [ll], such a linear temperature field with perfect heat conductivity of the crack surface will not generate any stresses, i.e. UC)= 0, K’;‘)= K{i) = K$f{= 0. Therefore, attention is focused on the second part of the solution, involving the Fredholm equation defined by eq. (21) for the dipole distribution function M(Q), Q on SC. It is known [lo] that for a dipole of unit strength applied at point (0,O) referred to in Fig. 2, the steady-state temperature solution is

’

TG=

2&(X2 + y2) *

(23)

Equation (23) is substituted into eq. (21) and x is replaced by x - CCfor the consideration of TG distributed along crack surface SC. Thus, the thermal boundary condition at the crack surfaces is given by DRsin/.? +

+a M(Z)ti AM(x) 2xK(x n)2 = 2pcra s -a

(24)

where Z is the integral variable at the crack surface and M(x) is in the units of energy per time. Temperature jumps across crack surfaces due to these distributed dipoles are Tc2’(x,0+) - T'2'(x, O-) =

I

a

I

a

Wx)

o.

I

Fig. 2. Heat flux through the crack border.

(25)

Stress intensity factors of cracks in finite plates

645

It is known [7,8, 121 that &=2KD

Jr@ +~)(a

-x)]Qsin/I

(26)

is the solution of Dosin/

+

+a M#)~ =O * _-(I 27&(x n)2 s

(27)

Thus, it is reasonable [9] to express M(x) as

M(x)=&u

-x)(u

-x,&0 f j

0

where d, are coefficients to be determined. Both eqs (26) and (28) implicitly satisfy the crack tip temperature conditions, T( f a, Of) = T( f a, O-) or M( +a) = 0. The coefficients 4 in eq. (28) can be determined by the following Fredhohn equation [13]: xDRsin/I

+A’=

+I?M(Z) - 1 -d.z 2nK s -,x-z

where A ’ is a constant to be determined.

Fig. 3. The disturbing temperature distribution (insulated).

(2%

K. Y. LAM et al.

646

-0.75 ml

-ml

-0.25 &on 0.25 0.50 Cl.75 toll

X/A Fig. 4. Temperature jumps at the crack surface.

e, = 1 e2 = f

ezj =

$ . .fJ(2i -

3) (j 2 2).

I-

Temperature ~st~bution TC2)(P,t) and tem~rature jumps at the crack surface for varying 1 can be calculated by substituting eq. (28) into eqs (20) and (25). Figure 3 shows the temperature distribution T”) caused by the insulated crack surface and Fig. 4 shows the temperature jumps at the crack surface. Physically, F2’ represents the disturbing effect on the linear temperature field due to the presence of a crack. This kind of disturbing temperature distribution is non-smooth and valid in the near-field region of the crack. As shown in Fig. 3, it can be seen that the temperature Tc2)is antisymmetric with respect to the line y = 0. In other words, the temperature F2) in the lower portion has the negative values of the corresponding quantity in the upper portion. Furthermore, the temperature TC2)must he zero outside the crack surface, Ix 12 a at y = 0. Inside the crack, the temperature is discontinuous and changes sign across the crack surfaces.

4. FINITE

ELEMENT

MODEL

The case of a plate with a central crack with the crack surfaces perfectly insulated (A = 0) is studied. All material properties and geometries of the plate are defined in Tables 1 and 2 and Fig. 5. Due to geometrical symmetry, it is suiIicient to solve only the right side of the plate, whose finite element mesh is displayed in Fig. 6. The far-field regular elements used are eight-noded isoparametric elements. The elements near the crack tips are singular elements, as shown in Fig. 7. The crack surfaces are prevented from overlapping by appropriate restraints. Table 1. Material properties Young’s modulus Material properties

E=2x

105MPa

Poisson’s ratio

Coefficient of thermal expansion

v = 0.3

0: = 6.5 x 10-6/oC

Stress intensity factors of cracks in finite plates Table 2. Stress intensity factors for the case of uniform tensile stress Geometrical cases fJ/W l/10

115

Stress intensity factors WPa(mm)‘“] for half-crack length 4 = 1 mm

HIW 0.5

iu; (ref. [ 141) 221.5

216.2

:I 3.0

212.6 212.6

208.6 210.0 209.1

0.5

265.8

262.4

:.II 3:o

248.1 239.2 235.6

244.2 235.1 230.8

&,(=I 2.4 1.9 1.2 1.4 1.3 1.6 1.7 2.1

In order to verify the adequacy of the mesh used, a uniform tensile stress a0 normal to the crack plane is applied to the top edge of the model, The bottom edge is restrained from moving in the direction of cr,. The mode I stress intensity factors Ki, calculated from the finite element solution of the &ack opening displacements for various geometries, are listed in Table 2. The values of Ki are compared with equivalent stress intensity factors K; obtained using the handbook [14]. For the cases considered here, the percentage error in the mode I stress intensity factor is less than 2.4%. Having established the adequacy of the mesh, it is then employed in the calculation of the thermal stress intensity factors where the appropriate temperatures are applied to the nodes. The strength of the temperature gradient 0 sin fi ~C/~) is given a nominal value of one. The deformed shape of the original mesh is displayed in Fig. 8 and the displacements of the crack surfaces are shown in Fig. 9. The opening and sliding deformations of the crack surfaces are easily seen in Fig. 9. The stress intensity factors can be evaluated using the above displacement IMd. The well-known relations~ps [I51 used to evaluate the stress intensity factors from nodal displacements are given by

(31)

Fig. 5. Geometry of the plate.

Fig. 6. Finite element mesh.

648

K. Y. LAM et al.

Fig. 7. Details of the finite element mesh near the crack tip.

Fig. 8. Deformed

mesh for the case of a/W = l/IO, H/W= 1.

and K,

8 38 ( - 2~ + 3) sin - - sin 2 2

(32) (2x -3)cos;+cos;

where K = 3 - 4v is for plane strain, r and B are the local coordinates as shown in Fig. 5, u and v are respectively the x and y components of the ~spla~ments vector in the region of the crack tip. The values of stress intensity factors for the case of fI sin #I = l”C/mm are shown in Table 3. The analytical solutions for g, er for an infinite plate obtained in ref. [2] for an insulated crack are listed below: KP=O

(33)

Kfr = $Ga(l + v)aJ(rca)R

Crack

center

R

sin /?

(34)

crad!tip

Fig. 9. Deformed shape of the crack surfaces for the case of a/W - l/IO, H/W = 1.

Stress intensity factors of cracks in finite plates

649

Table 3. Thermal stress intensity factors for each geometric case Geometrical cases

O/W l/10

115

HP

Stxess intensity factors ~Pa(~)‘~] for half-crack length o = I mm

R; (x10-2)

0.5 I.0 2.0 3.0

X-E

0.5

0.05 0.11 0.13 0.12

0:16 0.14

q

41

(0.0)

0.526 0.554 0.567 0.571 0.387 0.513 0.540 0.549

a

(0.576)

m= -4, I4 I 0.913 0.963 0.985 0.991 0.673 0.892 0.937 0.946

where G is the shear modulus, v is Poisson’s ratio and Q is the half-crack length. Using the values of material properties listed in Table 1 and the geometrical p~amete~ listed in Table 2, e, for the case of s2 sin /I = l”C/mm assumes the value of 0.576 MPa(mm)‘tz. Hence the analytical solutions can be compared with the values from the finite element solutions in order to determine the effect of finite geometry on the thermal stress intensity factors. 5. DISCUSSION A non-smoothly varying temperature gradient along the crack, especially around the crack tips, is induced by the perturbation of the temperature distribution in the neighborhood of a crack. This kind of temperature distribution is obtained numerically by solving a Fredholm integral equation. In the numerical analysis using the finite element method, the sin~a~ty elements are adopted around the crack tip. In the case of insulated crack surface heat conductivity, eight examples with different plate geometries are calculated. Based on a displacement field computed by MSC/NASTRAN, the mixed mode stress intensity factors are obtained (see Table 3). As shown in Table 3, Ki is negligible compared with K,, in all eight cases considered. This shows that the geometry of the plate does not significantly affect the mode I component. Table 3 also shows the variations of mode II stress intensity factor with increasing values of H/W for the cases of a/W = l/10 and l/5 respectively. For each case of the same a/W, K,, increases with increasing H/W. For each case of the same H/W, the value of Kit in the case of a/W = l/l0 is larger than that in the case of a/W = l/5. For the case of a/W = l/l0 and H/W = 3, the solution for 4, approaches that of an infinite plate since Ku/e, is 0.991. This result differs from the infinite solution by less than 1%. This shows that for the case of a/W > l/IO and H/W > 3, the ratio m( = lu,/J$ is essentially equal to one. From the above process of solutions, it can be seen that the disturbance due to the crack does not induce the mode I thermal stress intensity factor. For smaller plates, the mode II thermal stress intensity factor shows a very localized phenomenon. The value of the thermal stress intensity factor KiI depends on the geometric size of the plate. The value of Ku in a smaller plate is lower than that in a bigger plate. As the size of the plate increases, the value of thermal stress intensity factor E;,, is close to the solution of an infinite plate. The ratio m ( = K,, /fir) can be used as the geometric factor in thermal fracture analysis. From eqs (28) and (30), we can see that, for the case of partial heat conductivity of the crack surface, the nodal temperature loading will be reduced proportionally. According to the procedure obtained for the case of the insulated crack surface, the thermal stress and displacement fields for partial heat ~nducti~ty of the crack surface can be calculated easily. The above conclusions should be valid for the case of partial heat ~nducti~ty of the crack surface. This is more widely applicable in engineering structure design. 6. CONCLUSIONS The finite element method has been used successfully to calculate thermal stress intensity factors. The numerical analysis procedures for calculating the thermal stress intensity factors for a

650

K. Y. LAM et al.

finite plate containing a central crack are presented in this paper. The method essentially involves the calculation of the temperature distribution caused by the presence of the crack. This temperature distribution is used in the finite element analysis to obtain the stresses, displacements and thermal stress intensity factors. The results show that, for the cases studied here, the mode I thermal stress intensity factor K, is negligible compared with the mode II thermal stress intensity factor K,, and the geometry of a plate does not affect the mode I stress intensity factor significantly. This conclusion is consistent with refs [2,9]. Only the mode II thermal stress intensity factor is significant and shows a very localized phenomenon. For all cases studied here, the value of the thermal stress intensity factor K,, in smaller plates depends on the geometric size of the plate. As the size of the plate increases, the value of thermal stress intensity factor K,, approaches that of the infinite plate. The ratio m ( = K,/Kfi) can be used as the geometric factor in the realistic design of structures. REFERENCES G. C. Sih and H. Liebowitz, Mathematical Fundamentals, Fracture (Edited by H. Liebowitz), Vol. 2, pp. 67496. Academic Press, New York (1968). G. C. Sih, On the singular character of thermal stresses near a crack tip. J. appl. Mech. 29, 587-589 (1962). S. Davidson, G. Melrose and J. Tweed, The linear thermoelastic problem for a strip with a line crack parallel to its edges, in Basticity (Edited by G. Eason and R. W. Ogden), pp. 49-57. Ellis Honvood, England (1990). H. Sckine, Thermal stress singula~ties at tips of a crack in a ~mi-infinite medium under uniform heat flow. Engng Fracture Mech. 7, 713-729 (1975). W. H. Chen and K. Ting, Hybrid finite element analysis of transient thermoelastic fracture problems subjected to general heat transfer conditions. Comput. Mech. 4, I-10 (1989). P. C. Paris and G. C. Sih, Stress analysis of crack. Fracture Toughness Testing and its Applications, ASTM STP 381, 30-83 (1964). J. R. Barber, The disturbance of a uniform steady-state heat flux by a partially conducting plane crack. Znt. J. Mass Heat Transfer 19, 956-958 (1976). J. R. Barber, Steady-state stress caused by an imperfectly conducting ~nny-shah crack in an elastic solid. J. Therm. Stress 3, 77-83 (1979). A. Y. Kuo, Effects of crack surface heat conductance on stress intensity factors. J. appl. Mech. 57, 354-358 (1990). H. S. Carslaw and J. G. Jaeger, Conduction of Heat in Solids. Oxford University Press, London (1959). B. A. Boley and J. H. Weiner, Theory of Thermal Stresses. John Wiley, New York (1960). W. D. Collins, On the solution of some axisymmetric boundary value problems by means of integral equations--II. Mathematics 6, 120-133 (1959). N. I. M~khelishi~, Singtdar Integral .E+ations (Translated by J. R. Radok). Noordhoff, Groningen (1946). D. P. Rooke and D. J. Cartwright, &ompendium of Stress Intensity Factors. Her Majesty’s Stationary Office, London (1976). D. R. J. Owen and A. J. Fawkes, Engtneering Fracture Mechanics-Numerical Methodr and Applications. Pineridge Press, Swansea, U.K. (1983). (Received 28 August 1991)