- Email: [email protected]
Probabilistic theory for mixed mode fatigue crack growth in brittle plates with random cracks
Engineering Fracture Mechanics 66 (2000) 305±320
www.elsevier.com/locate/engfracmech
Probabilistic theory for mixed mode fatigue crack growth in brittle plates with random cracks D.W. Nicholson*, P. Ni, Y. Ahn Department of Mechanical, Materials and Aerospace Engineering, University of Central Florida, Orlando, FL 32816, USA Received 8 April 1999; received in revised form 3 January 2000; accepted 11 January 2000
Abstract A probabilistic fracture mechanics theory is proposed for fatigue life prediction in which a mixed-mode fatigue crack growth model is combined with extreme value probabilistic methods. The Sih±Barthelemy fatigue crack growth criterion based on strain energy intensity factor under uniaxial loading has been slightly modi®ed and extended to multiaxial loading including biaxial loading and torsion. A probability density function is computed for the number of cycles to failure, and scaling relations are also computed. Crack number, length, and orientation at the onset of fatigue are treated as random varixsables whose probability distributions are known. Numerical results are presented to demonstrate dependence on plate size, on the variance of the crack length distribution, on the multiaxial loading factors, and on the parameters of the fatigue model. 7 2000 Elsevier Science Ltd. All rights reserved.
1. Introduction High strength steels, alloys, and ceramics have been used increasingly in components of aerospace structures to achieve enhanced durability under severe operational stresses and temperatures. However, these materials often behave in a brittle manner and their fracture and fatigue behavior often exhibits randomness. Prediction of their fatigue life using probabilistic models is a subject of great interest. Acknowledging the statistical nature of fatigue crack growth in brittle materials, a number of probabilistic theories have been proposed [1±4]. Several well-known theories appear to introduce ad hoc models for fatigue crack growth, typically con®ned to Mode I. A well-known example is the Birnbaum± Saunders distribution [4,5]. It regards the incremental crack extension under a load as a random variable, with failure occurring when the crack attains a critical length; it does not use an established crack growth model from fracture mechanics, and it takes no account of the random number, length or orientation of * Corresponding author. Tel.: +1-407-823-2416; fax: +1-407-823-0208. E-mail address: [email protected] (D.W. Nicholson). 0013-7944/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 0 0 ) 0 0 0 0 7 - 2
306
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
pre-existing cracks. Issues of fracture mode and multiaxial loading are not addressed. A number of more recent but similar models are discussed in Ref. [4]. Many treatments emphasize the randomness of the applied load. Several treatments rigorously combine fatigue crack models from fracture mechanics with probabilistic methods, starting from the randomness of the initial defect [1,2]. However, to the best knowledge of the current authors, the fracture criteria used in existing treatments are based on Mode I relations. Furthermore, the number of cracks in the structure appears not to be treated as a random variable. It appears that few, if any, treatments are available which rigorously apply fatigue crack growth models from fracture mechanics to cracks under mixed mode and which accommodate random number, random length and random orientation of the cracks at the onset of the fatigue. In contrast, the present study introduces a probabilistic theory based on established criteria from fracture mechanics under mixed mode, accommodating the phenomena of crack kinking and curving. The crack number, length and orientation at the onset of fatigue are treated as random variables. The treatment of crack number is based on the cell concept introduced by Nicholson and Ni [6]. For the sake of illustration, a model due to Sih and Barthelemy [7], hereafter called the SB model, is used, with some further development, for mixed-mode fatigue crack growth under multiaxial loading, and eventual failure. An important assumption in the SB model, also used here, is that the kinked and curved crack can be replaced by an equivalent straight crack to be de®ned. Biaxial loading and torsion are accommodated. A discussion and comparison is given of several other mixed mode criteria, making it clear that they can be accommodated in the probabilistic theory by simply substituting the relations for the corresponding crack growth criteria. Recently, for static fracture under ®xed load, the authors introduced a conditional probability model based on a cell concept and the Sih [8] mixed-mode static fracture criterion. A probability density function (PDF) was derived for the strength of brittle plates with random cracks [6,9]. The strength distribution depended on plate size, the parameters of the crack length distribution, the parameters of the fracture criterion and load ratios. Unlike related earlier models not using the cell concept [10±12], the model is applicable without restriction on crack number or length. The present study extends the earlier probabilistic model proposed above to fatigue life prediction under ¯uctuating multiaxial load. For reasons to be explained, at the onset of fatigue crack number follows a binomial distribution, crack length follows a Gamma distribution, and crack orientation follows a uniform distribution. The SB model for mixed mode fatigue crack growth is slightly modi®ed and is extended to biaxial load and torsion using the solution of Eftis and Subramonian [13]. Under mixed mode loading, the crack kinks. In the SB model, the kinked crack is modeled using an equivalent straight crack connecting the two endpoints. Fracture occurs when at least one equivalent crack elongates and rotates to the point of satisfying the mixed-mode static fracture criterion, based on Sin [8], for multiaxial loads. The use of a cell notion together with conditional probability relations allows the current model to avoid limitation to large numbers of long cracks, unlike previous models such as in Ref. [10]. An extreme value probabilistic model is formulated, leading to an probability density function for Nf , the number of cycles to plate failure. The mean and variance of Nf are computed to show dependence on plate size, on the variance of the crack length distribution, on the load ratios, and on the parameters of the fatigue model. Scaling relations relating laboratory specimens and service plates are derived and computed. 2. Mixed mode fatigue model under multiaxial loadings 2.1. General Mixed mode fatigue is very complex, especially since the crack initially kinks and then follows a
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
307
curved path. For small kink lengths, the stress state at the crack tip has been analyzed by Hussain et al. [14] in Modes I and II, and Sih [15] in Mode III. Despite the complexity of the stress state, three scalar parameters have been proposed which are reasonably successful both in predicting the initial kink angle and the subsequent path under fatigue. A detailed comparison of the three parameters is given in Appendix A. The models are: (1) the Sih±Barthelemy [7] model based on the critical strain energy intensity factor, denoted here by Smin ; (2) the maximum circumferential stress criterion, denoted here by tmax [16]; and (3) the maximum strain energy release rate criterion, denoted as Gmax [14,15,17]. In the subsequent development a slight modi®cation of the Sih±Barthelemy [7] model is used. However, Gmax and tmax can be used just as easily by simply substituting the corresponding relations into the probabilistic model. In the SB model, crack growth is controlled by the ¯uctuation of the strain energy intensity factor S, initially expressed for uniaxial loading by Sih [8]. Here S is extended to biaxial load using the Eftis± Subramonian [13] solution, and also accommodates torsion. Incremental growth of an inclined crack according to the SB model is illustrated in Fig. 1. In the SB model a criterion was introduced to predict the growth of a fatigue crack. The criterion has mostly been applied to angled cracks under uniaxial tension. It is expressed in terms of the ¯uctuation of the strain energy intensity factor DS: This contrasts with the classical Paris relation in which the cracks grow along their axes, in Mode I. Also the Paris relation is expressed in terms of the ¯uctuation of the stress intensity factor. The main diculty posed by mixed mode is that the crack does not grow in a self-similar fashion, but kinks and curves. Accordingly, the SB model uses an equivalent straight crack illustrated in Fig. 2 and explained below. 2.2. Equivalent straight crack To resolve the relations between a and b, a detailed geometric representation is shown in Fig. 2. Let a' denote the half length of a line joining the endpoints of the crack, and let b be the angle of this line relative to the direction associated with s: The SB crack growth relation is now da 0 n a DS , dN
1
where a and n are material parameters, N is number of cycles, and S is the strain energy intensity factor (explained in Section 2.4) calculated using the equivalent crack length a ' and orientation b: Eq. (1) is easily generalized to incorporate thresholds. Suppose that the half-crack length is a and the inclined
Fig. 1. Straight line approximation of kinked crack.
308
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
angle of the crack is b at N cycles. The minimum value of S will correspond to an angle y0 which depends on b (and kB under biaxial load as well as kT under torsional load). The crack will grow by Da along y0 : The following relations have been derived in this study [9]: da
da 0 cos y0 b
2
and the angle b is changed incrementally according to db
da 0 tany0 : a0
3
By way of proof, the equivalent crack after incremental elongation and rotation has length a da 0 : From elementary geometry, using a 0 a, a da 0 a cos db da cos y0 ÿ db
4
furnishing Eq. (2) to ®rst order in db: Also, a 0 sin db da siny0 da 0 tany0 ,
5
furnishing Eq. (3) to ®rst order in db: 2.3. Fatigue crack growth law Numerical integration furnishes the relations a 0 N and b N , which de®ne the path of the crack endpoints. Static failure occurs at Nf , a 0 N and b N for which S Scr : Clearly, Nf depends on the applied stress states, as well as on the initial crack length and orientation, a0 and b0 : If there are M precracks in a plate in which a0 and b0 have known probability density functions, Nf will likewise have a probability distribution to be derived using extreme value methods explained below. For present purposes we introduce a convenient modi®cation of the SB model (Eq. (1)) as follows: dS n0 a 0 DS : dN
6
Fig. 2. Geometry for crack growth model, right-half only.
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
309
In the following, for the sake of illustration we will consider the case in which n 0 > 1: Hereafter the primes in a 0 and n ' will not be displayed. (However, dropping the primes does not imply that the coecients in Eqs. (1) and (6) have the same values!) For suitable choices of a and n, this model implies crack growth characteristics very similar to the SB model. Now DS Smax ÿ Smin , and we assume that Smin 0: The subscripts, max and min, and the D need not be displayed further, so that the crack growth law may be written as dS aS n : dN
7
In fact Eq. (7) is just a convenient example of the more general relation S=Scr L N , either in closed form or in numerical form. Eqs. (2) and (3) likewise imply such relations. The method described in the subsequent sections in fact applies to the general relation. To integrate Eq. (7), ®rst note that Nf sc dS a dN 8 n S0 S 0 We integrate both sides of Eq. (8) and rearrange it, so that initial value of strain energy intensity factor S0 is obtained as follows: 1 1ÿn S0 S1ÿn cr ÿ a 1 ÿ n Nf
9
where, Nf is the number of cycles to failure and Scr is the critical value of strain energy intensity factor, which is a material property. The above equation will be used to derive the fracture criterion based on the number of cycles to failure.
2.4. Strain energy intensity factor for biaxial loads and torsion The strain energy intensity factor S suited for biaxial loading and torsion was presented in detail in Ref. [9]. S satis®es a relation of the form S s 2 aC b, kB , kT , y0 b, kB , kT
10
where the biaxial and torsional factors kB and kT are illustrated in Fig. 3. Also, Fig. 4 demonstrates the crack propagation angles y0 for a range of biaxial parameters, which are very important for the probabilistic model presented here. In particular, if KI , KII and KIII are the stress intensity factors in Model I, II and III, respectively, then [8] 2 S a11 KI2 2a12 KI KII a22 KII2 a33 KIII
in which a11
a12
3 ÿ 4n cos y 1 cos y , 16pm
2sin ycos y ÿ 1 ÿ 2n 16p m,
11
310
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
Fig. 3. Biaxially±torsionally loaded center-cracked plate.
Fig. 4. Angle of crack propagation for a range of biaxial parameters.
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
a22
311
4 1 ÿ n 1 ÿ cos y 1 cos y 3 cos y ÿ 1 , 16pm
a33
4 : 16p m
12
Also n are m are the Poisson's ratio and the shear modulus, respectively. For mixed mode fatigue crack growth, the direction of crack propagation from the axis of the equivalent crack coincides with the angle y0 , which minimizes the strain energy intensity factor. In general, y0 must be evaluated numerically as a function of b, kB and kT : The growth and rotation of the equivalent crack leads to catastrophic crack propagation when the static fracture criterion is satis®ed, i.e., when S Scr where Scr is a material property called the critical strain p energy intensity p factor. For biaxial load and torsion, S is obtained as follows [9]. Let K s a k , K s a kII , and I I II p KIII s a kIII : Eq. (11) therefore becomes S s 2a C
13
where C a11 kI2 2a12 kI kII a22 k and
2 II
2 a33 kIII
14
p 1 kB ÿ 1 ÿ kB cos 2b , kI p 2 kII
p 1 ÿ kB sin 2b p , 2
kIII
p pkT sin b:
15
3. Extreme value probabilistic model 3.1. General: the cell concept Several earlier probabilistic fracture theories, for example Jayatilaka [10], obtained extreme value distributions by assuming that the number of cracks in the plate is large (and by using the right tail of the crack length distribution). We consider this assumption very restictive. In contrast, Ni [6] and Nicholson and Ni [9] introduced the cell concept. As illustrated in Fig. 5, we consider the plate to consist of Nc cells each of which has a probability of M=Nc of containing a center crack. Now asymptotic relations are achieved by letting the number of cells become large. However, most importantly, the number and size of the cracks can remain modest. The cell concept still has limitations. Currently each cell is assumed to contain at most one (initial) crack, and no cracks cross the cell borders, initially or during the fatigue process. Clearly, the characteristic cell dimensions must realistically be several times the critical crack length for static fracture under the assumed stress. Thus, strictly speaking, the asymptotic argument applies to plates which are very large compared to the critical crack length. The asymptotic argument has the
312
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
mathematical bene®t of furnishing `elegant' expressions appearing the extreme value distribution. However, if the conditions for a valid asymptotic argument are not satis®ed, a probabilistic model can still be formulated, using and evaluating comparatively cumbersome expressions. Under the current cell assumption, the crack number in the plate is random and satis®es a binomial distribution in which M is the expected number of cracks. The crack initial orientation follows a uniform probability distribution (i.e. isotropic). Now as a working hypothesis we assume that the crack initial length follows a normalized two-parameter probability distribution denoted as f a, and in particular we have selected a Gamma distribution for this purpose. Restrictions on admissible crack length distributions are derived by Ni [9] and Nicholson and Ni [6]. Suitable alternative distributions include the logarithmic normal and the Weibull distribution. All three distributions are commonly used in reliability theory [5] and have a very similar overall appearance. However, if the mean and variance are ®xed, they dier somewhat in the important tail region. At this time the best crack length distribution appears to be an open question. In any event, the logarithmic normal and Weibull distributions can be implemented in the current probabilistic formulation simply by substituting the corresponding probability density function. A joint probability density function, in which the crack length a and the crack angle b are statistically independent and all crack orientations are equally likely between 0 and p=2, is expressed as fab a, b
2 f a p
16
It should be noted that many other important random variables are not currently being accommodated, such as the critical fracture parameter (e.g. KIC) and the applied load.
3.2. Normalized crack length distribution The Gamma distribution with parameters p > 0 and q > 0 is expressed as 1 q p pÿ1 x pÿ1 exp ÿ xdx a exp ÿ qa , G p f a G p 0 The mean and variance of a are
Fig. 5. A plate with N cells and M inclined center-cracks M < N).
17
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
a p=q,
Sa p=q 2
313
18
To study the eect of crack length dispersion (variance) it is convenient to introduce the normalized The mean of c is unity, and the variance of c is denoted as Sc Sa =a 2 ÿ 1: crack length c a=a:
3.3. Failure probability for one crack The probability of failure at number of cycles N for one crack is given by FN N 1 ÿ
p af 2
0
0
2 f a da db p
19
where af is obtained as a function of b and the number of cycles N using both the fracture criterion and the mixed mode fatigue crack growth model. In particular, " # 1 1ÿn S1ÿn a 1 ÿ n cr N 1ÿ af 2 f s C S1ÿn cr
20
where C has been stated in Eqs. (14) and (15).
3.4. Extreme value distribution The following extreme value argument is adapted from the static fracture case in Ref. [6]. The probability that one crack will fail in N cycles is FN N , and that it will survive is 1 ÿ FN N : The probability that a cell will survive is 1 ÿ FN N M=Nc : The probability that at least one cell will fail in N cycles is the complement of the probability that all cells will survive. Now assuming that the number of cells is large compared to the expected number of cracks, we obtain the probability of failure for the whole plate after N cycles as p 1 ÿ exp ÿ MFN GN , N large 21 Pr amax > af , 0 < b < 2 1 ÿ exp ÿ M provided that at least one cell contains a crack [6]). Now, if Nfk is the number of cycles to failure of the kth cell, Eq. (21) is interpreted after a change of variables as N gN Nf dNf 22 Pr Nf min < N GN N 0
in which Nf min minfNf1 , Nf2 , . . . ,Nki , . . . ,NfN g:
23
The corresponding probability density function, which is the desired extreme value distribution for the number of cycles to failure, is given by gN
dGN : dNf
24
314
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
4. Scaling relations It is of great practical interest to compute scaling relations between service and laboratory structures. Consider two plate populations with areas A1 > A2 and M1 > M2 , and assume that A1 =A2 M1 =M2 : The subscripts 1 and 2 denote the service and laboratory plates, respectively. It is assumed that both experience the same value of s, kB and kT , and contain pre-existing cracks coming from the same population. The ratios of the means and variances of the number of cycles to failure are to be calculated as a function of M, i.e.: mN1 E M1 , Sc , kB , kT , mN2 E M2 , Sc , kB , kT
SN1 Var M1 , Sc , kB , kT : SN2 Var M2 , Sc , kB , kT
25
(The mean crack lengths cancel owing to the use of the normalized crack length distribution.) The mean E Nf and variance Var Nf of the number of cycles to failure for the whole plate are given (under suitable conditions [6]) by 1 1 1 dGN 1 ÿ GN dNf Nf gN dNf Nf dNf 26 E Nf dNf 0 0 0 ÿ E Nf2
1 0
Nf2 gN dNf
1 0
Nf2
ÿ Var Nf E Nf2 ÿ E 2 Nf 2
dGN dNf 2 dNf 1 0
1 0
Nf 1 ÿ GN dNf
Nf 1 ÿ GN dNf ÿ
1 0
1 ÿ GN dNf
27 2
:
28
For purposes of integration we introduce a variable change Nf
S1ÿn cr 1 ÿ x1ÿn a 1 ÿ n
29
S1ÿn cr x ÿn dx a
30
and dNf ÿ
where x is expressed as "
a 1 ÿ n Nf x 1ÿ S1ÿn cr
#
1 1ÿn
:
31
Then x
2C cas : S1ÿn cr
32
Note that Nf 4 0 as x4 1,while Nf 4 1 as x 4 0: Substituting Eqs. (28) and (29) into Eqs. (26) and (27) yields
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
E Nf
S1ÿn cr a
1 0
1 ÿ GN x ÿ x ÿn dx
ÿ S 2 1 ÿ n E Nf2 cr2 a 1 ÿ n
1 0
ÿ 1 ÿ GN x 1 ÿ x1ÿn ÿ x ÿn dx:
315
33
34
5. Numerical results A computer program has been developed for calculating the foregoing ratios. For the sake of illustration, we select AISI 4130 steel. Sih and Macdonald [18] gave the material properties of this material: the elastic modulus is E 30 106 psi, Poisson's ratio is n 0:25 and the critical strain energy intensity factor is Scr 32:2 lb/in. Also let the applied stress be s 170 ksi, a 0:01, and n 2: The foregoing expressions involve triple integration. They have been implemented using 20 node Gaussian quadrature for the ®nite intervals and 41 node Kronrod quadrature for the in®nite intervals [19]. Three dierent stress states are considered in the numerical calculations. For a uniaxial stress state, the mean and variance of number of cycles to failure only depend on the crack number M and the variance of the crack length distribution Sc : Therefore, kB 0, and kT 0 in Eq. (8). We have taken the values of the parameter Sc to be 2, 2.5, 3, and 4. With varying Sc and M1 =M2 , the ratios of the mean of the number of cycles to failure under uniaxial stress have calculated. As can be seen in Figs 6 and 7 along the kb 0 axis, the results show that as the ratio of expected crack
Fig. 6. The ratio of mean number of cycles-to-failure under biaxial-torsional loading.: M1 =M2 = service to laboratory expected crack number; kb = biaxiality factor; mN1 =mN2 = service to laboratory strength mean strengths.
316
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
Fig. 7. The ratio of variance of number of cycles-to-failure under biaxial-torsional loading:M1 =M2 = service to laboratory expected crack number; kb = biaxiality factor; SN1 =SN2 = service to laboratory strength variance.
number increases, the ratios of mean and variance of number of cycles to failure decrease. The model thus predicts the following: if the mean number of cycles to failure mN2 is measured in laboratory plates, and if the crack length distributions of service and laboratory plates are assumed to be the same, then the mean and variance of number of cycles to failure of a service plate decreases as its size (or number of crack M1 gets larger. Regarding the mean strength ratio, the results are as expected. In intuitive terms, for the same crack density and crack probability density function, the probability is greater in larger plates that a crack is present which is long compared to the critical static crack length. Thus, service plates will have lower mean strengths, as predicted. It likewise can be seen that the decreasing variance ratio should be expected. For very large plates, the probability is high that each plate contains a crack whose length approaches the critical value for static fracture, and hence has a small value of Nf : (This is consistent with the fact that the mean values of Nf are lower in the service plates.) The distribution of Nf will thus tend to become clustered near a small values. On the other hand, for small plates, there is a very low probability of a crack near the static fracture value, and there is no reason for the values of Nf to cluster in any particular range. For biaxial stress kT 0), the mean and variance of the number of cycles-to-failure depend not only on the crack number M and the variance of crack length distribution Sc , but also on the biaxiality parameter kB as shown in Fig. 4. By taking the variance of crack length distribution Sc to be 3 and varying kB and M1 =M2 , the ratios of the mean and variance of the number of cycles to failure under
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
317
biaxial stress state are obtained. For biaxial-torsional stress, the mean and variance of the number of cycles-to-failure depend on the crack number M, the variance of crack length distribution Sc , the biaxiality pparameter kB , and the torsional parameter kT : Taking the value of parameters Sc 3 and kT 1= 3, and varying kB and M1 =M2 , the ratios of the mean and variance of the number of cycles to failure under biaxial-torsional stress state are found as illustrated in Figs. 6 and 7.
6. Conclusions In the current study, an extreme value probabilistic theory has been formulated for fatigue life prediction. Previous studies have used adhoc models for crack growth, while the current study uses models established in fracture mechanics. In particular, the Sih±Barthelemy model for mixedmode fatigue crack growth has been further developed and extended to multiaxial loads. The extreme value distribution of interest has been derived based on the cell concept. The asymptotic distribution is obtained by allowing the number of cells to become large, rather than the number of cracks as in earlier studies. The model is thus suitable for a ®nite number of cracks. The crack number, crack length and crack orientation are considered random variables. The crack number is given by a binomial distribution. The crack length is given by a normalized Gamma distribution Scaling relations for number of cycles to failure in laboratory against service plates are formulated and computed.
Appendix A. Comparison of models for mixed mode fracture and fatigue In mixed mode fracture the crack typically kinks when fracture occurs. Several criteria accommodating kinking have been proposed for mixed mode. All three criteria are consistent with limited available data, consisting of kink angles measured in angled cracks under uniaxial tension. In addition, the crack propagation path of a fatigue problem under uniaxial cyclic load is computed by adopting the eective straight crack assumption. All three criteria show reasonable agreement with measured paths. The authors have found little well-documented experimental data on fatigue under multiaxial loads. For the crack under general loading as shown in Fig. 2(a), the three parameters are based on the following expressions: 1 3 y 1 3y 3 y 3 3y ÿ KII sin sin A1 syy p KI cos cos 4 2 4 2 4 2 4 2 2pr S
1 k 16pm
ÿcos y 1 cos y KI2
1 1 k 1 1 ÿ cos y 1 sin y 2cos y ÿ k ÿ 1 KI KII 8pm 16pm
1 2 K cos y 3cos y ÿ 1 KII2 4pm III
A2
318
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
y 2 k1 1 pÿy p ÿ 1 3cos 2 y KI2 8sin y cos yKI KII 9 ÿ 5cos 2 y KII2 G 2 2m 3 cos y py y 1 pÿy p 2 KIII 2m p y
A3
where r, y denote the polar coordinates as shown in Fig. 1(b), m is the shear modulus, and k 3 ÿ 4n or k 3 ÿ n= 1 n for plane strain and plane stress, respectively n is the Poisson ratio). The stress intensity factors KI , KII and KIII are given by KI
s p pa 1 kb ÿ 1 ÿ kb cos 2b , 2
A4
KII
s p p a 1 ÿ kb sin 2b, 2
A5
p KIII kt s pa sin j:
A6
Mixed mode criteria based on the three parameters are, respectively, called the maximum circumferential stress criterion tmax criterion) [13], the minimum strain energy density criterion Smin criterion) [4,10] and the maximum strain energy release rate criterion Gmax criterion) [2,17]. The tmax and Gmax criteria predict that the crack will propagate along the direction that maximizes syy and G, respectively, whereas the Smin criterion predicts that the crack will propagate along the direction that minimizes S. The predicted crack propagation angle y0 for various values of b is computed w.r.t. y and shown in Table 1 along with experimental data. All three models give reasonable and comparable agreement. A direct application of a mixed mode fracture criterion is the mixed mode fatigue crack propagation problem. The crack propagation path for the fatigue problem under a mixed mode loading condition can be computed by taking the following steps. (1) Given the initial crack shape (we assume the initial crack to be a straight crack), ®nd the angle of crack growth using a mixedmode fracture criterion. (2) Along the direction of crack growth, extend the crack by increment da; (3) For the new crack shape (generally, a kinked crack), ®nd the new angle of crack growth. (4) Repeat Steps 2 and 3 until the desired path is obtained. Although the algorithm described above is a straightforward one, Step 3 poses a particularly dicult problem since ®nding the crack growth direction for a kinked crack is extremely hard both theoretically and computationally. To overcome this diculty the assumption of the eective straight crack is used again, in which the kinked crack is replaced by a straight crack joining the two crack tips of the kinked crack. By employing this assumption, we can carry out Step 3 by using a mixed-mode fracture criterion on the eective straight crack. The computation of the fatigue crack growth path is now reduced to a series of eective straight crack problems where the crack orientation Table 1 Computed and measured values of the crack propagation angle y0 b
308
408
508
608
708
808
tmax [3] Smin [3] n 1=3) Gmax Experiment [3]
ÿ60.28 ÿ63.58 ÿ64.78 ÿ62.48
ÿ55.78 ÿ56.78 ÿ60.08 ÿ55.68
ÿ50.28 ÿ49.58 ÿ54.38 ÿ51.68
ÿ43.28 ÿ41.58 ÿ46.68 ÿ43.18
ÿ33.28 ÿ31.88 ÿ35.98 ÿ30.78
ÿ19.38 ÿ18.58 ÿ20.78 ÿ17.38
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
319
Fig. 8. Computed and experimental fatigue crack propagation paths.
angle for the new eective straight crack is obtained by use of Eqs. (A1)±(A3). In the present investigation, we compute the crack propagation path for the fatigue problem under a uniaxial cyclic loading using the tmax , Smin and Gmax mixed mode fracture criteria. The crack orientation angle of the initial crack is assumed to be 158, 308, 458, and 608 and the results are shown in Fig. 8 against data reported by Mageed and Pandey [20]. All three mixed mode criteria seem to produce reasonably good, and very comparable results.
References [1] Wirshung PH. Probabilistic fatigue analysis. In: Sundararajan C, editor. Probabilistic structural mechanics handbook. New York: Chapman and Hall, 1995. [2] Harris DO. Probabilistic fracture mechanics''. In: Sundararajan C, editor. Probabilistic structural mechanics handbook. New York: Chapman and Hall, 1995.
320
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
[3] Bloom JM, Ekvall JC. Probabilistic fracture mechanics and fatigue methods: application for structural design and maintenance. Baltimore, MD:ASTM STP 798, 1983. [4] Castillo E. Extreme value distributions in engineering. NY: Academic Press, 1988. [5] Mann NR, Schafer RE, Singpurwalla ND. Methods for statistical analysis of reliability and life data. New York: Wiley, 1974. [6] Nicholson DW, Ni P. Extreme value probabilistic theory for mixed mode brittle fracture. Engineering fracture mechanics 1997;58:121±32. [7] Sih G, Barthelemy BM. Mixed mode fatigue crack growth predictions. Engineering Fracture Mechanics 1980;13:439±51. [8] Sih G. Strain energy density factor applied mixed-mode crack problems. International Journal of Fracture 1974;10(3):305±21. [9] Ni P. Probabilistic theory for mixed-mode brittle fracture and fatigue with random cracks. Doctoral dissertation, University of Central Florida, Orlando, FL, 1996. [10] Jayatilaka A, De S. Fracture of engineering brittle materials. London: Applied Science, 1979. [11] Ichikawa M. Probability evaluation of brittle solids under polyaxial stress states based on maximum energy release rate criterion. Engineering Fracture Mechanics 1995;51(4):629±35. [12] She S, Landes JD. Statistical analysis of fracture in graphite. International Journal of Fracture 1993;63:189±200. [13] Eftis J, Subramonian N. The inclined crack under biaxial load. Engineering Fracture Mechanics 1978;10:43±67. [14] Hussain MA, Pu SL, Underwood J. Strain energy release rate for a crack under combined Mode I and Mode II. In: Fracture analysis. ASTM STP 560, New York: American Society or Testing and Materials, 1974. p. 2±28. [15] Sih GC. Stress distribution near internal crack tips for longitudinal shear problems. Journal of Applied Mechanics 1965;32:51±8. [16] Erdogan F, Sih GC. On the crack extension in plates under loading and transverse shear. ASME Journal of Basic Engineering 1963;85:519±27. [17] Ichikawa M, Tanaka S. A critical analysis of the relationship between the energy release rate and the stress intensity factors for non-coplanar crack extension under combined mode loading. International Journal of Fracture 1982;18:19±28. [18] Sih G, Macdonald B. Fracture mechanics applied to engineering problemsÐstrain energy density fracture criterion. Engineering Fracture Mechanics 1974;6:361±86. [19] Squire WB. Integration for engineers and scientists. New York: American Elsevier, 1970. [20] Mageed AMA, Pandey RK. Fatigue crack closure in kinked cracks and path of crack propagation. International Journal of Fracture 1990;44:R39±R42.
www.elsevier.com/locate/engfracmech
Probabilistic theory for mixed mode fatigue crack growth in brittle plates with random cracks D.W. Nicholson*, P. Ni, Y. Ahn Department of Mechanical, Materials and Aerospace Engineering, University of Central Florida, Orlando, FL 32816, USA Received 8 April 1999; received in revised form 3 January 2000; accepted 11 January 2000
Abstract A probabilistic fracture mechanics theory is proposed for fatigue life prediction in which a mixed-mode fatigue crack growth model is combined with extreme value probabilistic methods. The Sih±Barthelemy fatigue crack growth criterion based on strain energy intensity factor under uniaxial loading has been slightly modi®ed and extended to multiaxial loading including biaxial loading and torsion. A probability density function is computed for the number of cycles to failure, and scaling relations are also computed. Crack number, length, and orientation at the onset of fatigue are treated as random varixsables whose probability distributions are known. Numerical results are presented to demonstrate dependence on plate size, on the variance of the crack length distribution, on the multiaxial loading factors, and on the parameters of the fatigue model. 7 2000 Elsevier Science Ltd. All rights reserved.
1. Introduction High strength steels, alloys, and ceramics have been used increasingly in components of aerospace structures to achieve enhanced durability under severe operational stresses and temperatures. However, these materials often behave in a brittle manner and their fracture and fatigue behavior often exhibits randomness. Prediction of their fatigue life using probabilistic models is a subject of great interest. Acknowledging the statistical nature of fatigue crack growth in brittle materials, a number of probabilistic theories have been proposed [1±4]. Several well-known theories appear to introduce ad hoc models for fatigue crack growth, typically con®ned to Mode I. A well-known example is the Birnbaum± Saunders distribution [4,5]. It regards the incremental crack extension under a load as a random variable, with failure occurring when the crack attains a critical length; it does not use an established crack growth model from fracture mechanics, and it takes no account of the random number, length or orientation of * Corresponding author. Tel.: +1-407-823-2416; fax: +1-407-823-0208. E-mail address: [email protected] (D.W. Nicholson). 0013-7944/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 0 0 ) 0 0 0 0 7 - 2
306
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
pre-existing cracks. Issues of fracture mode and multiaxial loading are not addressed. A number of more recent but similar models are discussed in Ref. [4]. Many treatments emphasize the randomness of the applied load. Several treatments rigorously combine fatigue crack models from fracture mechanics with probabilistic methods, starting from the randomness of the initial defect [1,2]. However, to the best knowledge of the current authors, the fracture criteria used in existing treatments are based on Mode I relations. Furthermore, the number of cracks in the structure appears not to be treated as a random variable. It appears that few, if any, treatments are available which rigorously apply fatigue crack growth models from fracture mechanics to cracks under mixed mode and which accommodate random number, random length and random orientation of the cracks at the onset of the fatigue. In contrast, the present study introduces a probabilistic theory based on established criteria from fracture mechanics under mixed mode, accommodating the phenomena of crack kinking and curving. The crack number, length and orientation at the onset of fatigue are treated as random variables. The treatment of crack number is based on the cell concept introduced by Nicholson and Ni [6]. For the sake of illustration, a model due to Sih and Barthelemy [7], hereafter called the SB model, is used, with some further development, for mixed-mode fatigue crack growth under multiaxial loading, and eventual failure. An important assumption in the SB model, also used here, is that the kinked and curved crack can be replaced by an equivalent straight crack to be de®ned. Biaxial loading and torsion are accommodated. A discussion and comparison is given of several other mixed mode criteria, making it clear that they can be accommodated in the probabilistic theory by simply substituting the relations for the corresponding crack growth criteria. Recently, for static fracture under ®xed load, the authors introduced a conditional probability model based on a cell concept and the Sih [8] mixed-mode static fracture criterion. A probability density function (PDF) was derived for the strength of brittle plates with random cracks [6,9]. The strength distribution depended on plate size, the parameters of the crack length distribution, the parameters of the fracture criterion and load ratios. Unlike related earlier models not using the cell concept [10±12], the model is applicable without restriction on crack number or length. The present study extends the earlier probabilistic model proposed above to fatigue life prediction under ¯uctuating multiaxial load. For reasons to be explained, at the onset of fatigue crack number follows a binomial distribution, crack length follows a Gamma distribution, and crack orientation follows a uniform distribution. The SB model for mixed mode fatigue crack growth is slightly modi®ed and is extended to biaxial load and torsion using the solution of Eftis and Subramonian [13]. Under mixed mode loading, the crack kinks. In the SB model, the kinked crack is modeled using an equivalent straight crack connecting the two endpoints. Fracture occurs when at least one equivalent crack elongates and rotates to the point of satisfying the mixed-mode static fracture criterion, based on Sin [8], for multiaxial loads. The use of a cell notion together with conditional probability relations allows the current model to avoid limitation to large numbers of long cracks, unlike previous models such as in Ref. [10]. An extreme value probabilistic model is formulated, leading to an probability density function for Nf , the number of cycles to plate failure. The mean and variance of Nf are computed to show dependence on plate size, on the variance of the crack length distribution, on the load ratios, and on the parameters of the fatigue model. Scaling relations relating laboratory specimens and service plates are derived and computed. 2. Mixed mode fatigue model under multiaxial loadings 2.1. General Mixed mode fatigue is very complex, especially since the crack initially kinks and then follows a
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
307
curved path. For small kink lengths, the stress state at the crack tip has been analyzed by Hussain et al. [14] in Modes I and II, and Sih [15] in Mode III. Despite the complexity of the stress state, three scalar parameters have been proposed which are reasonably successful both in predicting the initial kink angle and the subsequent path under fatigue. A detailed comparison of the three parameters is given in Appendix A. The models are: (1) the Sih±Barthelemy [7] model based on the critical strain energy intensity factor, denoted here by Smin ; (2) the maximum circumferential stress criterion, denoted here by tmax [16]; and (3) the maximum strain energy release rate criterion, denoted as Gmax [14,15,17]. In the subsequent development a slight modi®cation of the Sih±Barthelemy [7] model is used. However, Gmax and tmax can be used just as easily by simply substituting the corresponding relations into the probabilistic model. In the SB model, crack growth is controlled by the ¯uctuation of the strain energy intensity factor S, initially expressed for uniaxial loading by Sih [8]. Here S is extended to biaxial load using the Eftis± Subramonian [13] solution, and also accommodates torsion. Incremental growth of an inclined crack according to the SB model is illustrated in Fig. 1. In the SB model a criterion was introduced to predict the growth of a fatigue crack. The criterion has mostly been applied to angled cracks under uniaxial tension. It is expressed in terms of the ¯uctuation of the strain energy intensity factor DS: This contrasts with the classical Paris relation in which the cracks grow along their axes, in Mode I. Also the Paris relation is expressed in terms of the ¯uctuation of the stress intensity factor. The main diculty posed by mixed mode is that the crack does not grow in a self-similar fashion, but kinks and curves. Accordingly, the SB model uses an equivalent straight crack illustrated in Fig. 2 and explained below. 2.2. Equivalent straight crack To resolve the relations between a and b, a detailed geometric representation is shown in Fig. 2. Let a' denote the half length of a line joining the endpoints of the crack, and let b be the angle of this line relative to the direction associated with s: The SB crack growth relation is now da 0 n a DS , dN
1
where a and n are material parameters, N is number of cycles, and S is the strain energy intensity factor (explained in Section 2.4) calculated using the equivalent crack length a ' and orientation b: Eq. (1) is easily generalized to incorporate thresholds. Suppose that the half-crack length is a and the inclined
Fig. 1. Straight line approximation of kinked crack.
308
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
angle of the crack is b at N cycles. The minimum value of S will correspond to an angle y0 which depends on b (and kB under biaxial load as well as kT under torsional load). The crack will grow by Da along y0 : The following relations have been derived in this study [9]: da
da 0 cos y0 b
2
and the angle b is changed incrementally according to db
da 0 tany0 : a0
3
By way of proof, the equivalent crack after incremental elongation and rotation has length a da 0 : From elementary geometry, using a 0 a, a da 0 a cos db da cos y0 ÿ db
4
furnishing Eq. (2) to ®rst order in db: Also, a 0 sin db da siny0 da 0 tany0 ,
5
furnishing Eq. (3) to ®rst order in db: 2.3. Fatigue crack growth law Numerical integration furnishes the relations a 0 N and b N , which de®ne the path of the crack endpoints. Static failure occurs at Nf , a 0 N and b N for which S Scr : Clearly, Nf depends on the applied stress states, as well as on the initial crack length and orientation, a0 and b0 : If there are M precracks in a plate in which a0 and b0 have known probability density functions, Nf will likewise have a probability distribution to be derived using extreme value methods explained below. For present purposes we introduce a convenient modi®cation of the SB model (Eq. (1)) as follows: dS n0 a 0 DS : dN
6
Fig. 2. Geometry for crack growth model, right-half only.
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
309
In the following, for the sake of illustration we will consider the case in which n 0 > 1: Hereafter the primes in a 0 and n ' will not be displayed. (However, dropping the primes does not imply that the coecients in Eqs. (1) and (6) have the same values!) For suitable choices of a and n, this model implies crack growth characteristics very similar to the SB model. Now DS Smax ÿ Smin , and we assume that Smin 0: The subscripts, max and min, and the D need not be displayed further, so that the crack growth law may be written as dS aS n : dN
7
In fact Eq. (7) is just a convenient example of the more general relation S=Scr L N , either in closed form or in numerical form. Eqs. (2) and (3) likewise imply such relations. The method described in the subsequent sections in fact applies to the general relation. To integrate Eq. (7), ®rst note that Nf sc dS a dN 8 n S0 S 0 We integrate both sides of Eq. (8) and rearrange it, so that initial value of strain energy intensity factor S0 is obtained as follows: 1 1ÿn S0 S1ÿn cr ÿ a 1 ÿ n Nf
9
where, Nf is the number of cycles to failure and Scr is the critical value of strain energy intensity factor, which is a material property. The above equation will be used to derive the fracture criterion based on the number of cycles to failure.
2.4. Strain energy intensity factor for biaxial loads and torsion The strain energy intensity factor S suited for biaxial loading and torsion was presented in detail in Ref. [9]. S satis®es a relation of the form S s 2 aC b, kB , kT , y0 b, kB , kT
10
where the biaxial and torsional factors kB and kT are illustrated in Fig. 3. Also, Fig. 4 demonstrates the crack propagation angles y0 for a range of biaxial parameters, which are very important for the probabilistic model presented here. In particular, if KI , KII and KIII are the stress intensity factors in Model I, II and III, respectively, then [8] 2 S a11 KI2 2a12 KI KII a22 KII2 a33 KIII
in which a11
a12
3 ÿ 4n cos y 1 cos y , 16pm
2sin ycos y ÿ 1 ÿ 2n 16p m,
11
310
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
Fig. 3. Biaxially±torsionally loaded center-cracked plate.
Fig. 4. Angle of crack propagation for a range of biaxial parameters.
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
a22
311
4 1 ÿ n 1 ÿ cos y 1 cos y 3 cos y ÿ 1 , 16pm
a33
4 : 16p m
12
Also n are m are the Poisson's ratio and the shear modulus, respectively. For mixed mode fatigue crack growth, the direction of crack propagation from the axis of the equivalent crack coincides with the angle y0 , which minimizes the strain energy intensity factor. In general, y0 must be evaluated numerically as a function of b, kB and kT : The growth and rotation of the equivalent crack leads to catastrophic crack propagation when the static fracture criterion is satis®ed, i.e., when S Scr where Scr is a material property called the critical strain p energy intensity p factor. For biaxial load and torsion, S is obtained as follows [9]. Let K s a k , K s a kII , and I I II p KIII s a kIII : Eq. (11) therefore becomes S s 2a C
13
where C a11 kI2 2a12 kI kII a22 k and
2 II
2 a33 kIII
14
p 1 kB ÿ 1 ÿ kB cos 2b , kI p 2 kII
p 1 ÿ kB sin 2b p , 2
kIII
p pkT sin b:
15
3. Extreme value probabilistic model 3.1. General: the cell concept Several earlier probabilistic fracture theories, for example Jayatilaka [10], obtained extreme value distributions by assuming that the number of cracks in the plate is large (and by using the right tail of the crack length distribution). We consider this assumption very restictive. In contrast, Ni [6] and Nicholson and Ni [9] introduced the cell concept. As illustrated in Fig. 5, we consider the plate to consist of Nc cells each of which has a probability of M=Nc of containing a center crack. Now asymptotic relations are achieved by letting the number of cells become large. However, most importantly, the number and size of the cracks can remain modest. The cell concept still has limitations. Currently each cell is assumed to contain at most one (initial) crack, and no cracks cross the cell borders, initially or during the fatigue process. Clearly, the characteristic cell dimensions must realistically be several times the critical crack length for static fracture under the assumed stress. Thus, strictly speaking, the asymptotic argument applies to plates which are very large compared to the critical crack length. The asymptotic argument has the
312
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
mathematical bene®t of furnishing `elegant' expressions appearing the extreme value distribution. However, if the conditions for a valid asymptotic argument are not satis®ed, a probabilistic model can still be formulated, using and evaluating comparatively cumbersome expressions. Under the current cell assumption, the crack number in the plate is random and satis®es a binomial distribution in which M is the expected number of cracks. The crack initial orientation follows a uniform probability distribution (i.e. isotropic). Now as a working hypothesis we assume that the crack initial length follows a normalized two-parameter probability distribution denoted as f a, and in particular we have selected a Gamma distribution for this purpose. Restrictions on admissible crack length distributions are derived by Ni [9] and Nicholson and Ni [6]. Suitable alternative distributions include the logarithmic normal and the Weibull distribution. All three distributions are commonly used in reliability theory [5] and have a very similar overall appearance. However, if the mean and variance are ®xed, they dier somewhat in the important tail region. At this time the best crack length distribution appears to be an open question. In any event, the logarithmic normal and Weibull distributions can be implemented in the current probabilistic formulation simply by substituting the corresponding probability density function. A joint probability density function, in which the crack length a and the crack angle b are statistically independent and all crack orientations are equally likely between 0 and p=2, is expressed as fab a, b
2 f a p
16
It should be noted that many other important random variables are not currently being accommodated, such as the critical fracture parameter (e.g. KIC) and the applied load.
3.2. Normalized crack length distribution The Gamma distribution with parameters p > 0 and q > 0 is expressed as 1 q p pÿ1 x pÿ1 exp ÿ xdx a exp ÿ qa , G p f a G p 0 The mean and variance of a are
Fig. 5. A plate with N cells and M inclined center-cracks M < N).
17
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
a p=q,
Sa p=q 2
313
18
To study the eect of crack length dispersion (variance) it is convenient to introduce the normalized The mean of c is unity, and the variance of c is denoted as Sc Sa =a 2 ÿ 1: crack length c a=a:
3.3. Failure probability for one crack The probability of failure at number of cycles N for one crack is given by FN N 1 ÿ
p af 2
0
0
2 f a da db p
19
where af is obtained as a function of b and the number of cycles N using both the fracture criterion and the mixed mode fatigue crack growth model. In particular, " # 1 1ÿn S1ÿn a 1 ÿ n cr N 1ÿ af 2 f s C S1ÿn cr
20
where C has been stated in Eqs. (14) and (15).
3.4. Extreme value distribution The following extreme value argument is adapted from the static fracture case in Ref. [6]. The probability that one crack will fail in N cycles is FN N , and that it will survive is 1 ÿ FN N : The probability that a cell will survive is 1 ÿ FN N M=Nc : The probability that at least one cell will fail in N cycles is the complement of the probability that all cells will survive. Now assuming that the number of cells is large compared to the expected number of cracks, we obtain the probability of failure for the whole plate after N cycles as p 1 ÿ exp ÿ MFN GN , N large 21 Pr amax > af , 0 < b < 2 1 ÿ exp ÿ M provided that at least one cell contains a crack [6]). Now, if Nfk is the number of cycles to failure of the kth cell, Eq. (21) is interpreted after a change of variables as N gN Nf dNf 22 Pr Nf min < N GN N 0
in which Nf min minfNf1 , Nf2 , . . . ,Nki , . . . ,NfN g:
23
The corresponding probability density function, which is the desired extreme value distribution for the number of cycles to failure, is given by gN
dGN : dNf
24
314
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
4. Scaling relations It is of great practical interest to compute scaling relations between service and laboratory structures. Consider two plate populations with areas A1 > A2 and M1 > M2 , and assume that A1 =A2 M1 =M2 : The subscripts 1 and 2 denote the service and laboratory plates, respectively. It is assumed that both experience the same value of s, kB and kT , and contain pre-existing cracks coming from the same population. The ratios of the means and variances of the number of cycles to failure are to be calculated as a function of M, i.e.: mN1 E M1 , Sc , kB , kT , mN2 E M2 , Sc , kB , kT
SN1 Var M1 , Sc , kB , kT : SN2 Var M2 , Sc , kB , kT
25
(The mean crack lengths cancel owing to the use of the normalized crack length distribution.) The mean E Nf and variance Var Nf of the number of cycles to failure for the whole plate are given (under suitable conditions [6]) by 1 1 1 dGN 1 ÿ GN dNf Nf gN dNf Nf dNf 26 E Nf dNf 0 0 0 ÿ E Nf2
1 0
Nf2 gN dNf
1 0
Nf2
ÿ Var Nf E Nf2 ÿ E 2 Nf 2
dGN dNf 2 dNf 1 0
1 0
Nf 1 ÿ GN dNf
Nf 1 ÿ GN dNf ÿ
1 0
1 ÿ GN dNf
27 2
:
28
For purposes of integration we introduce a variable change Nf
S1ÿn cr 1 ÿ x1ÿn a 1 ÿ n
29
S1ÿn cr x ÿn dx a
30
and dNf ÿ
where x is expressed as "
a 1 ÿ n Nf x 1ÿ S1ÿn cr
#
1 1ÿn
:
31
Then x
2C cas : S1ÿn cr
32
Note that Nf 4 0 as x4 1,while Nf 4 1 as x 4 0: Substituting Eqs. (28) and (29) into Eqs. (26) and (27) yields
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
E Nf
S1ÿn cr a
1 0
1 ÿ GN x ÿ x ÿn dx
ÿ S 2 1 ÿ n E Nf2 cr2 a 1 ÿ n
1 0
ÿ 1 ÿ GN x 1 ÿ x1ÿn ÿ x ÿn dx:
315
33
34
5. Numerical results A computer program has been developed for calculating the foregoing ratios. For the sake of illustration, we select AISI 4130 steel. Sih and Macdonald [18] gave the material properties of this material: the elastic modulus is E 30 106 psi, Poisson's ratio is n 0:25 and the critical strain energy intensity factor is Scr 32:2 lb/in. Also let the applied stress be s 170 ksi, a 0:01, and n 2: The foregoing expressions involve triple integration. They have been implemented using 20 node Gaussian quadrature for the ®nite intervals and 41 node Kronrod quadrature for the in®nite intervals [19]. Three dierent stress states are considered in the numerical calculations. For a uniaxial stress state, the mean and variance of number of cycles to failure only depend on the crack number M and the variance of the crack length distribution Sc : Therefore, kB 0, and kT 0 in Eq. (8). We have taken the values of the parameter Sc to be 2, 2.5, 3, and 4. With varying Sc and M1 =M2 , the ratios of the mean of the number of cycles to failure under uniaxial stress have calculated. As can be seen in Figs 6 and 7 along the kb 0 axis, the results show that as the ratio of expected crack
Fig. 6. The ratio of mean number of cycles-to-failure under biaxial-torsional loading.: M1 =M2 = service to laboratory expected crack number; kb = biaxiality factor; mN1 =mN2 = service to laboratory strength mean strengths.
316
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
Fig. 7. The ratio of variance of number of cycles-to-failure under biaxial-torsional loading:M1 =M2 = service to laboratory expected crack number; kb = biaxiality factor; SN1 =SN2 = service to laboratory strength variance.
number increases, the ratios of mean and variance of number of cycles to failure decrease. The model thus predicts the following: if the mean number of cycles to failure mN2 is measured in laboratory plates, and if the crack length distributions of service and laboratory plates are assumed to be the same, then the mean and variance of number of cycles to failure of a service plate decreases as its size (or number of crack M1 gets larger. Regarding the mean strength ratio, the results are as expected. In intuitive terms, for the same crack density and crack probability density function, the probability is greater in larger plates that a crack is present which is long compared to the critical static crack length. Thus, service plates will have lower mean strengths, as predicted. It likewise can be seen that the decreasing variance ratio should be expected. For very large plates, the probability is high that each plate contains a crack whose length approaches the critical value for static fracture, and hence has a small value of Nf : (This is consistent with the fact that the mean values of Nf are lower in the service plates.) The distribution of Nf will thus tend to become clustered near a small values. On the other hand, for small plates, there is a very low probability of a crack near the static fracture value, and there is no reason for the values of Nf to cluster in any particular range. For biaxial stress kT 0), the mean and variance of the number of cycles-to-failure depend not only on the crack number M and the variance of crack length distribution Sc , but also on the biaxiality parameter kB as shown in Fig. 4. By taking the variance of crack length distribution Sc to be 3 and varying kB and M1 =M2 , the ratios of the mean and variance of the number of cycles to failure under
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
317
biaxial stress state are obtained. For biaxial-torsional stress, the mean and variance of the number of cycles-to-failure depend on the crack number M, the variance of crack length distribution Sc , the biaxiality pparameter kB , and the torsional parameter kT : Taking the value of parameters Sc 3 and kT 1= 3, and varying kB and M1 =M2 , the ratios of the mean and variance of the number of cycles to failure under biaxial-torsional stress state are found as illustrated in Figs. 6 and 7.
6. Conclusions In the current study, an extreme value probabilistic theory has been formulated for fatigue life prediction. Previous studies have used adhoc models for crack growth, while the current study uses models established in fracture mechanics. In particular, the Sih±Barthelemy model for mixedmode fatigue crack growth has been further developed and extended to multiaxial loads. The extreme value distribution of interest has been derived based on the cell concept. The asymptotic distribution is obtained by allowing the number of cells to become large, rather than the number of cracks as in earlier studies. The model is thus suitable for a ®nite number of cracks. The crack number, crack length and crack orientation are considered random variables. The crack number is given by a binomial distribution. The crack length is given by a normalized Gamma distribution Scaling relations for number of cycles to failure in laboratory against service plates are formulated and computed.
Appendix A. Comparison of models for mixed mode fracture and fatigue In mixed mode fracture the crack typically kinks when fracture occurs. Several criteria accommodating kinking have been proposed for mixed mode. All three criteria are consistent with limited available data, consisting of kink angles measured in angled cracks under uniaxial tension. In addition, the crack propagation path of a fatigue problem under uniaxial cyclic load is computed by adopting the eective straight crack assumption. All three criteria show reasonable agreement with measured paths. The authors have found little well-documented experimental data on fatigue under multiaxial loads. For the crack under general loading as shown in Fig. 2(a), the three parameters are based on the following expressions: 1 3 y 1 3y 3 y 3 3y ÿ KII sin sin A1 syy p KI cos cos 4 2 4 2 4 2 4 2 2pr S
1 k 16pm
ÿcos y 1 cos y KI2
1 1 k 1 1 ÿ cos y 1 sin y 2cos y ÿ k ÿ 1 KI KII 8pm 16pm
1 2 K cos y 3cos y ÿ 1 KII2 4pm III
A2
318
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
y 2 k1 1 pÿy p ÿ 1 3cos 2 y KI2 8sin y cos yKI KII 9 ÿ 5cos 2 y KII2 G 2 2m 3 cos y py y 1 pÿy p 2 KIII 2m p y
A3
where r, y denote the polar coordinates as shown in Fig. 1(b), m is the shear modulus, and k 3 ÿ 4n or k 3 ÿ n= 1 n for plane strain and plane stress, respectively n is the Poisson ratio). The stress intensity factors KI , KII and KIII are given by KI
s p pa 1 kb ÿ 1 ÿ kb cos 2b , 2
A4
KII
s p p a 1 ÿ kb sin 2b, 2
A5
p KIII kt s pa sin j:
A6
Mixed mode criteria based on the three parameters are, respectively, called the maximum circumferential stress criterion tmax criterion) [13], the minimum strain energy density criterion Smin criterion) [4,10] and the maximum strain energy release rate criterion Gmax criterion) [2,17]. The tmax and Gmax criteria predict that the crack will propagate along the direction that maximizes syy and G, respectively, whereas the Smin criterion predicts that the crack will propagate along the direction that minimizes S. The predicted crack propagation angle y0 for various values of b is computed w.r.t. y and shown in Table 1 along with experimental data. All three models give reasonable and comparable agreement. A direct application of a mixed mode fracture criterion is the mixed mode fatigue crack propagation problem. The crack propagation path for the fatigue problem under a mixed mode loading condition can be computed by taking the following steps. (1) Given the initial crack shape (we assume the initial crack to be a straight crack), ®nd the angle of crack growth using a mixedmode fracture criterion. (2) Along the direction of crack growth, extend the crack by increment da; (3) For the new crack shape (generally, a kinked crack), ®nd the new angle of crack growth. (4) Repeat Steps 2 and 3 until the desired path is obtained. Although the algorithm described above is a straightforward one, Step 3 poses a particularly dicult problem since ®nding the crack growth direction for a kinked crack is extremely hard both theoretically and computationally. To overcome this diculty the assumption of the eective straight crack is used again, in which the kinked crack is replaced by a straight crack joining the two crack tips of the kinked crack. By employing this assumption, we can carry out Step 3 by using a mixed-mode fracture criterion on the eective straight crack. The computation of the fatigue crack growth path is now reduced to a series of eective straight crack problems where the crack orientation Table 1 Computed and measured values of the crack propagation angle y0 b
308
408
508
608
708
808
tmax [3] Smin [3] n 1=3) Gmax Experiment [3]
ÿ60.28 ÿ63.58 ÿ64.78 ÿ62.48
ÿ55.78 ÿ56.78 ÿ60.08 ÿ55.68
ÿ50.28 ÿ49.58 ÿ54.38 ÿ51.68
ÿ43.28 ÿ41.58 ÿ46.68 ÿ43.18
ÿ33.28 ÿ31.88 ÿ35.98 ÿ30.78
ÿ19.38 ÿ18.58 ÿ20.78 ÿ17.38
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
319
Fig. 8. Computed and experimental fatigue crack propagation paths.
angle for the new eective straight crack is obtained by use of Eqs. (A1)±(A3). In the present investigation, we compute the crack propagation path for the fatigue problem under a uniaxial cyclic loading using the tmax , Smin and Gmax mixed mode fracture criteria. The crack orientation angle of the initial crack is assumed to be 158, 308, 458, and 608 and the results are shown in Fig. 8 against data reported by Mageed and Pandey [20]. All three mixed mode criteria seem to produce reasonably good, and very comparable results.
References [1] Wirshung PH. Probabilistic fatigue analysis. In: Sundararajan C, editor. Probabilistic structural mechanics handbook. New York: Chapman and Hall, 1995. [2] Harris DO. Probabilistic fracture mechanics''. In: Sundararajan C, editor. Probabilistic structural mechanics handbook. New York: Chapman and Hall, 1995.
320
D.W. Nicholson et al. / Engineering Fracture Mechanics 66 (2000) 305±320
[3] Bloom JM, Ekvall JC. Probabilistic fracture mechanics and fatigue methods: application for structural design and maintenance. Baltimore, MD:ASTM STP 798, 1983. [4] Castillo E. Extreme value distributions in engineering. NY: Academic Press, 1988. [5] Mann NR, Schafer RE, Singpurwalla ND. Methods for statistical analysis of reliability and life data. New York: Wiley, 1974. [6] Nicholson DW, Ni P. Extreme value probabilistic theory for mixed mode brittle fracture. Engineering fracture mechanics 1997;58:121±32. [7] Sih G, Barthelemy BM. Mixed mode fatigue crack growth predictions. Engineering Fracture Mechanics 1980;13:439±51. [8] Sih G. Strain energy density factor applied mixed-mode crack problems. International Journal of Fracture 1974;10(3):305±21. [9] Ni P. Probabilistic theory for mixed-mode brittle fracture and fatigue with random cracks. Doctoral dissertation, University of Central Florida, Orlando, FL, 1996. [10] Jayatilaka A, De S. Fracture of engineering brittle materials. London: Applied Science, 1979. [11] Ichikawa M. Probability evaluation of brittle solids under polyaxial stress states based on maximum energy release rate criterion. Engineering Fracture Mechanics 1995;51(4):629±35. [12] She S, Landes JD. Statistical analysis of fracture in graphite. International Journal of Fracture 1993;63:189±200. [13] Eftis J, Subramonian N. The inclined crack under biaxial load. Engineering Fracture Mechanics 1978;10:43±67. [14] Hussain MA, Pu SL, Underwood J. Strain energy release rate for a crack under combined Mode I and Mode II. In: Fracture analysis. ASTM STP 560, New York: American Society or Testing and Materials, 1974. p. 2±28. [15] Sih GC. Stress distribution near internal crack tips for longitudinal shear problems. Journal of Applied Mechanics 1965;32:51±8. [16] Erdogan F, Sih GC. On the crack extension in plates under loading and transverse shear. ASME Journal of Basic Engineering 1963;85:519±27. [17] Ichikawa M, Tanaka S. A critical analysis of the relationship between the energy release rate and the stress intensity factors for non-coplanar crack extension under combined mode loading. International Journal of Fracture 1982;18:19±28. [18] Sih G, Macdonald B. Fracture mechanics applied to engineering problemsÐstrain energy density fracture criterion. Engineering Fracture Mechanics 1974;6:361±86. [19] Squire WB. Integration for engineers and scientists. New York: American Elsevier, 1970. [20] Mageed AMA, Pandey RK. Fatigue crack closure in kinked cracks and path of crack propagation. International Journal of Fracture 1990;44:R39±R42.