Simulation of 3D fatigue crack propagation using an implicit time integration scheme

International Journal of Fatigue 50 (2013) 2–9

Contents lists available at SciVerse ScienceDirect

International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue

Simulation of 3D fatigue crack propagation using an implicit time integration scheme W. Weber ⇑, P. Steinmann, G. Kuhn Chair of Applied Mechanics, University of Erlangen-Nuremberg, Egerlandstr. 5, 91058 Erlangen, Germany

a r t i c l e

i n f o

Article history: Received 15 May 2011 Received in revised form 18 June 2012 Accepted 10 July 2012 Available online 28 July 2012 Keywords: 3D crack growth simulation 3D crack growth criterion Predictor–corrector scheme

a b s t r a c t Due to the non-linear behavior of fatigue crack growth an incremental numerical procedure has to be applied for the simulation. In the present work a new implicit time integration scheme for arbitrary loading conditions is presented ensuring accurate simulation results. Within this scheme special focus is on the evolution of the stress field between two successive discrete crack fronts of the incremental procedure. By the consideration of the changing stress field the numerically obtained crack fronts are optimized with respect to their shape and location. The whole procedure is implemented in the framework of a predictor–corrector scheme. Two examples are shown to demonstrate the benefit of the presented scheme. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Cracks, which trace back to damaging during the manufacturing process, are often the origin of failure of machinery parts. The collapse of safety-relevant parts results in perilous situations for human beings. Therefore, the fracture mechanics assessment of these structures becomes more important in the design process. In presence of cyclic loading conditions fatigue crack propagation is very critical, because crack growth occurs for lower stresses compared to static loadings. Due to the non-linear nature of crack growth an incremental procedure has to be applied for the simulation of crack propagation. Each increment starts with a complete stress analysis including the determination of the fracture mechanics parameters along the crack front. Then, the 3D crack growth criterion is evaluated for the calculation of the crack extension and the kink angle, cf. Fig. 1a and b. Finally, the discretization is adjusted to the new crack geometry for the next incremental loop, see Fig. 1c. For the stress analysis the dual boundary element method (dual-BEM) in terms of the dual discontinuity formulation is applied [1–3] which is a special formulation for crack problems. Due to its nature the BEM captures stress concentration problems better than the finite element method [4]. Moreover, the modification of the mesh during the simulation of crack propagation is easier by using boundary elements compared to volume orientated methods. By the application of the adaptive cross approximation

⇑ Corresponding author. Tel.: +49 9131 8528502. E-mail address: [email protected] (W. Weber). 0142-1123/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijfatigue.2012.07.003

the numerical complexity of the stress analysis is reduced significantly [5–8]. Since the 3D crack growth criterion is mainly based on the fracture mechanics parameters the calculation of these values follows the stress analysis. In the present context of linear-elastic fracture mechanics these parameters are at first the stress intensity factors (SIFs). Additionally, the non-singular T-stresses are taken into account. The SIFs and the T-stresses are determined by an extrapolation method from the singular stress field in the crack front near field [9]. In case of surface breaking cracks corner singularities have to be analyzed at the surface intersection points, cf. Section 2.2. The main topic of this work is on the determination of the crack growth. Starting from a given initial crack front the new crack front is defined relative to the initial one by the crack deflection and the crack extension. These data have to be calculated for the relevant points of the crack – namely the crack front nodes of the discretization. For the calculation of the crack deflection the maximum tangential stress (MTS) criterion [10] has been established. In the present context an extended formulation of this criterion considering the non-singular T-stresses is utilized which shows a stable behavior in the simulation [11]. The crack extension is obtained by the evaluation of a crack propagation rate formulation. At first only the fracture mechanics parameters of the initial crack front are known. Therefore, the crack growth is predicted in a linear way. In the next incremental loop the fracture mechanics parameters of the new crack front are known additionally. Therewith, the evolution of the changing stress field between successive crack fronts is approximated and the error of the prediction is estimated. By minimizing the linearization error the non-linear nature of crack growth is considered and the predicted crack front is

3

W. Weber et al. / International Journal of Fatigue 50 (2013) 2–9

Fig. 1. Three steps of an increment.

corrected. As a consequence this predictor–corrector procedure yields an optimization of the new crack front with respect to its location and shape. Finally, the numerical model has to be updated at the end of each increment for the next loop. In case of a predictor step large crack extensions are present at the whole crack front. Therefore, the gap between the two crack fronts is closed by inserting a new row of elements, see Fig. 1. If the crack front is corrected the new crack front is located near the previous one. In this case the crack front nodes are moved to their new position. The accuracy of the presented predictor–corrector-scheme and its applicability on the transition from crack arrest to crack growth is shown on examples. 2. Singularities along a 3D crack front Along a 3D crack front two different types of singularities can occur. If the crack front is smooth at the crack front point under consideration the well-known wedge singularity characterized by pffiffiffi the 1= r-stress-singularity is present. At kinks of the crack front and especially at the intersection points with the free surface the type of singularity is not a priori known. Here, 3D corner singularities have to be analyzed. 2.1. Wedge singularity pffiffiffi At smooth parts of the crack front the classical 1= r -stress-singularity is present. For each point P along such a part of the crack front the stresses in the near field are given in terms of the local coordinate system, cf. Fig. 2, by [12]:

rij ðr; u; PÞ ¼

III X pffiffiffi K M ðPÞ M pffiffiffiffiffiffiffiffiffi fij ðuÞ þ T ij ðPÞ þ O r : 2 p r M¼I

ð1Þ

The SIFs KM(P) denote the intensity of the typical stress singularity for each crack opening mode M and fijM ðuÞ are the corresponding angular functions. Tij(P) are the constant T-stresses. The SIFs as well as the T-stresses are calculated by an error controlled local extrapolation technique based on the numerically obtained stresses in the crack near field [9]. In case of mixed mode conditions an equivalent stress intensity factor is required for the assessment against brittle rupture or for the determination of the crack growth under cyclic loading conditions. A reasonable approach is the utilization of the energy release rate (ERR). The classical formulation of Irwin [13]

GðPÞ ¼

1m 1 ½K III ðPÞ2 ½½K I ðPÞ2 þ ½K II ðPÞ2  þ 2G G

ð2Þ

with the Poisson ratio m and the shear modulus G is based on the assumption of co-planar crack propagation. However, if mode II is active a crack deflection is present and Irwins formula gives incorrect results. To overcome this problem Nuismers formulation [14] via the maximum ERR including the mode II kinking is extended by mode III in this work. For the consideration of mode III the type of mode III crack propagation has to be defined at first. Here, it is assumed that the mode II kinking occurs immediately whereas the mode III rotation is formed continuously during crack growth and the crack front rotates in one piece [15], see Fig. 3a. Therefore, the formation of mode III facets (cf. Fig. 3b) is excluded in this formulation. Since the ERR is defined for an infinitesimal crack extension no mode III rotation is observable. As a consequence, only the mode II kinking (kink angle u) has to be taken into account and the ERR is written as:

Gðu; PÞ ¼

2 1m 1  K III ðu; PÞ : ½½K I ðu; PÞ2 þ ½K II ðu; PÞ2  þ 2G 2G

ð3Þ

Here, K M denote the SIFs of the kinked branch crack which are defined as follows:

cos u2 ½K I ðPÞ½1 þ cos u  3K II ðPÞ sin u; 2 u cos 2 K II ðu; PÞ ¼ ½K I ðPÞ sin u þ K II ðPÞ½3 cos u  1 and 2 K I ðu; PÞ ¼

K III ðu; PÞ ¼ K III ðPÞ cos Fig. 2. Crack front coordinate system.

u 2

:

ð4Þ ð5Þ ð6Þ

4

W. Weber et al. / International Journal of Fatigue 50 (2013) 2–9

Fig. 3. Basic propagation possibilities of mode III.

The derivation of Eqs. (4) and (5) is shown in [14] and Eq. (6) is obtained accordingly. In detail, K M are determined by identification from the current stress field in the coordinate system of the branch crack. Since the branch crack is virtual K M are independent from the length of the branch crack. In order to obtain the maximum ERR Eq. (3) has to be maximized with respect to u. The direction uG max ðPÞ of the maximum ERR Gmax ðPÞ is determined via the condition @Gðu; PÞ=@ u ¼ 0 under the restriction @ 2 Gðu; PÞ=@ u2 < 0 and Gmax ðPÞ is defined as Gmax ðPÞ ¼ GðuG max ; PÞ. Finally, the equivalent SIF is calculated under the assumption of pure mode I conditions of crack growth by

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G Gmax ðPÞ K eq ðPÞ ¼ : 1m

2.2. Corner singularity At special points like corners and especially at the intersection points of the crack front with the outer boundary corner singularities have to be considered. In contrast to points at a smooth crack front part the type of singularity is not a priori known and depends on the elastic material parameters as well as on the geometrical situation around the singular point, e.g. point Q in Fig. 5. The stresses around such a point are expressed with respect to a spherical coordinate system, c.f. Fig. 5, in asymptotic form [16]:

rij ðq; u; h; Q Þ ¼

N0 X K L ðQÞqaL 1 g Lij ðu; h; Q Þ þ T ij ðQ Þ L

ð7Þ

The direction of the maximum energy release rate and the resulting fracture limit surface Keq = KIc are illustrated in Fig. 4. In case of cyclic loading conditions the cyclic equivalent SIF is required to decide whether crack arrest, stable or unstable crack growth is present. For mixed mode conditions the equivalent cyclic SIF is determined via min max DK eq ðPÞ ¼ K max eq ðPÞ  K eq ðPÞ ¼ ½1  RK eq ðPÞ:

ð8Þ

K max and K min eq eq denote the maximum and the minimum equivalent SIF during a load cycle. R is the ratio of K min to K max eq eq ðR ¼ min max K eq =K eq Þ and for proportional loading conditions it is equal to the applied stress ratio.

  þ O qaN0 þ1 1 :

ð9Þ

The stress field is mainly characterized by the asymptotic exponents aL, which are determined by the solution of a quadratic eigenvalue problem [17,18]. K L ðQ Þ are the generalized stress intensity factors of the singularities [aL  1] while g Lij ðu; h; Q Þ are the corresponding eigenvectors. Reasonable values of aL are between 0.5 and 1.0, [17,18]. If aL pffiffiffi is greater than 0.5 the singularity is weaker compared to the 1= r stress-singularity at smooth parts of the crack front. This can be expressed within the classical SIF-concept by setting KM(Q) = 0. If aL is less than 0.5 the opposite situation is present. Here, the SIF tends to infinity. For the simulation of fatigue crack growth the results of the comparison of numerical analyses and experimental observations

Fig. 4. Direction of maximum ERR and fracture limit surface for m = 0.3.

5

W. Weber et al. / International Journal of Fatigue 50 (2013) 2–9

criterion it is assumed that the crack grows in the (x1, x2)-plane of the local crack front coordinate system (cf. Fig. 2) perpendicular to the maximum tangential stress rmax u ðPÞ. The kink angle is determined by @ ru(P)/@ u = 0 under the restriction @ 2 ru(P)/@ u2 < 0. This results in [10]

Fig. 5. Coordinate system at a corner of the crack front.

are taken into account. It has been shown that a crack propagates in the vicinity of a corner within only few load cycles in such a way pffiffiffi that the 1= r -stress-singularity is present even at these points [19]. This fact corresponds to a special geometric situation around the singular point. It means that inside the volume a smooth crack front is present. Moreover, a special crack front intersection angle c is observed at the surface breaking points. This angle can be determined by a singularity analysis [18]. Within the crack growth algorithm which is presented in Sections 3 and 4 the knowledge on corner singularity is considered as follows: The crack front intersection angle c that corresponds to aL = 0.5 is determined by singularity analysis. After calculation of the crack propagation for crack front points that are characterized by wedge singularity the geometry of the crack front close to the intersection points is defined by a pure geometric approach in that way that the angle c between crack front and outer surface is present. 3. 3D crack growth criterion The 3D crack growth criterion is set up by formulas that describe the mechanism of crack growth under the present state of stress. It is assumed that at each point of the crack front the crack starts growing in radial direction perpendicular to the maximum tangential stress if the equivalent cyclic SIF exceeds the lower threshold value of crack arrest. In the presence of mode III it is assumed that the crack front rotates gradually in one piece without the formation of facets.

0

1

uMTS ðPÞ ¼ 2 arctan B @

2K II ðPÞ ffiC qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A: 2 2 K I ðPÞ þ ½K I ðPÞ þ 8½K II ðPÞ

This angle describes the tangent of a deflecting crack relative to the current crack path. The utilization of this angle within the simulation for a curved crack path in association with a finite straight incremental length Da(P) leads to a new crack front point P MTS nþ1 , which differs from the real curved crack path, cf. Fig. 6. Moreover, a ‘‘zigzag’’ path around the smooth curved crack path may occur. In order to obtain a better approximation of the real curved crack path and to avoid a ‘‘zigzag’’ path the T-stresses are considered additionally. Then, based on the MTS-criterion the curved crack path is characterized by the path of the maximum tangential stress ru with respect to the moving crack front coordinate system, ðx1 ; x2 Þ and (r⁄,u⁄), respectively. Depending on the local crack extension Da(P) and the T-stress T11(P), the kink angle uMTS⁄(P) is obtained by the evaluation of

@ ru ðPÞ ¼0 @u

ð11Þ

under the restriction @ 2 ru ðPÞ=@ u2 < 0 with

ru ðu; Da; PÞ ¼

rr ðu; Da; PÞ þ ru ðu; Da; PÞ 1

þ 2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ½rr ðu; Da; PÞ  ru ðu; Da; PÞ2 þ 4sru ðu; Da; PÞ ð12Þ

and [12]





K I ðPÞ u 3u rr ðu; Da; PÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5 cos  cos 2 2 4 2 p DaðPÞ   K II ðPÞ u 3u þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5 sin þ 3 sin 2 2 4 2 p DaðPÞ þ T 11 ðPÞ cos2 u; 

ð13Þ



K I ðPÞ u 3u 3 cos þ cos ru ðu; Da; PÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 4 2 p DaðPÞ   K II ðPÞ u 3u þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 sin  3 sin 2 2 4 2 p DaðPÞ 2 þ T 11 ðPÞ sin u;

3.1. Crack deflection For the direction of crack growth the direction of the maximum energy release rate according to Section 2.1 would be an appropriate choice. However, the application of the kink angle according to this criterion leads to a ‘‘zigzag’’ path, which is typical for criteria for the crack deflection that utilize only the SIFs [11]. Due to the definition of the energy release rate based on an infinitesimal crack extension the consideration of a finite crack extension within the simulation is hardly possible. For mixed mode I/II problems the direction of the maximum energy release rate is identical to the direction of the maximum tangential stress [14]. Therefore, the maximum tangential stress (MTS) criterion [10] is taken into account for the determination of the kink angle. Moreover, the influence of the crack extension on the kink angle can be considered within this criterion. The MTS criterion traces back to Erdogan and Sih [10] and has been verified successfully by experiments [20]. Within this

ð10Þ

Fig. 6. Determination of the kink angle.

ð14Þ

6

W. Weber et al. / International Journal of Fatigue 50 (2013) 2–9





K I ðPÞ u 3u sin þ sin sru ðu; Da; PÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 4 2 p DaðPÞ   K II ðPÞ u 3u þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos þ 3 cos  T 11 ðPÞ 2 2 4 2 p DaðPÞ  cos u sin u: ð15Þ The solution of Eq. (11) cannot be given explicitly. It has to be determined indirect by an iterative method, e.g. the Newtons method. The MTS-criterion has been developed for plane problems. Therefore, it can be utilized for plane problems and crack configurations with a small mode III portion. Usually, these conditions are present in crack problems from industrial practice. However, if a relevant mode III portion has to be considered an extended criterion is required. In this case, the direction of the maximum energy release rate which has been introduced in Section 2.1 can be utilized for general 3D mixed mode problems hazarding a ‘‘zigzag’’ path. 3.2. Crack extension The crack growth speed in the present case of fatigue crack growth is described by the crack propagation rate. It is defined as the crack extension per load cycle N. The crack growth rate depends mainly on the equivalent cyclic SIF. This fact has been realized by Paris and his co-workers [21] who gave a first formula for this behavior:

da ðPÞ ¼ C ½DK eq ðPÞm : dN

ð16Þ

The crack propagation rate da/dN depends on the equivalent cyclic SIF DKeq(P) as well as on the material parameters C and m. It has to be considered that DKeq(P) itself depends on the crack length a. This formulation is valid if the cyclic equivalent SIF far exceeds the lower threshold of crack growth and is far lower than the critical SIF. The formula of Erdogan and Ratwani [22] considers the threshold value DKth of fatigue crack growth as well as the fracture toughness KIc:

da ½DK eq ðPÞ  DK th m ðPÞ ¼ C : dN ½1  RK Ic  DK eq ðPÞ

ð17Þ

This formulation includes all phases of fatigue crack growth which occur at a crack front – crack arrest, stable crack propagation and rupture. Based on this definition the crack growth rate has been described by even more extensive formulations, e.g. the modified Forman/Mettu-equation [23]. 4. Time integration scheme The basic equations of crack growth in the section above are written in a time derivative form. Therefore, these equations have to be integrated in time in order to obtain the crack growth. This time integration is done in the present paper in the framework of a predictor–corrector-scheme in order to ensure a high accuracy. 4.1. Predictor step Starting from an initial crack front a new one is predicted. In order to define the new crack front relative to the initial one the crack deflection u(P) and the crack extension Da(P) are required for each point of the crack front. At first, the equivalent cyclic SIF DKeq(P) based on the energy release rate is determined, see Section 2.1. Since only the state of stress of the current crack front is known, it is assumed that the cyclic equivalent SIF DKeq(P) is

constant between the initial and the predicted crack front. Therefore, the crack extension length Da(P) is calculated by the evaluation of a crack propagation rate formulation with the loading parameters of the current crack front:



DaðPÞ ¼

 da ðDK eq ðPÞÞ DN lc : dN

ð18Þ

DNlc denotes a user-specified incremental number of load cycles. For the crack propagation rate da/dN(DKeq(P)), e.g. Eq. (16) or Eq. (17) can be applied. By the knowledge of the local crack extension Da(P) the crack  deflection uMTS ðPÞ is determined by the extended MTS-criterion (11) considering the non-singular T-stresses. Finally, at the surface breaking points the crack front is defined by ensuring the crack front intersection angle determined by the analysis of corner singularities. 4.2. Corrector step Within the next incremental loop of the crack growth simulation the boundary value problem for the predicted crack front configuration is solved and the cyclic equivalent SIFs DK Beq ðPÞ of the predicted crack front are known. These values are utilized to estimate the linearization error of the predictor step. In order to prevent confusion between the fracture mechanics parameters of the initial crack front and the predicted one in the following the values of the initial crack front are marked with the superior A and the values of the predicted crack front with the superior B, respectively. For all kinds of crack configurations the SIF in 2D is given by

DK ¼ Drn

pffiffiffiffiffiffiffiffi p a YðaÞ;

ð19Þ

where Drn denotes the cyclic normal stress, a represents the crack length and Y(a) is a geometric function which depends on the crack length a. In order to transfer this condition to the present 3D problem a virtual initial crack length a0(P) is introduced for each point P of the initial crack front. Therewith the virtual crack length of the predicted crack front reads as a0(P) + Da(P). For the determination of the virtual initial crack length, Eq. (19) is evaluated at the initial (superscript A) and the predicted (superscript B) crack front. This leads to the equations:

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi DK Aeq ðPÞ ¼ Drn p a0 ðPÞ Yða0 ðPÞÞ; p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi DK Beq ðPÞDrn p ½a0 ðPÞ þ DaðPÞ Yða0 ðPÞ þ DaðPÞÞ:

ð20Þ ð21Þ

The calculation of the ratio of both equations eliminates the cyclic normal stress from the system of equations. In the next step, a statement on the geometric function has to be made. Usually, this function Y(a) is a polynomial of the ratio of the crack length a to a characteristic length b of the structure,

hai2 h a in a YðaÞ ¼ c0 þ c1 þ c2 þ . . . þ cn : b b b

ð22Þ

The constants of the polynomial are denoted by ci. At the beginning of a crack growth simulation the crack length a is very small compared to the structure which means a  b (comparable to a crack in an infinite structure). Therefore, the geometric function is approximately constant:

Yða0 ðPÞÞ  Yða0 ðPÞ þ DaðPÞÞ:

ð23Þ

Otherwise, if the crack length a has increased, the local crack extension is much smaller than the crack length, Da  a and the value of the geometrical function is again approximately constant. By the utilization of this approximation the virtual initial crack length a0(P) is obtained by

7

W. Weber et al. / International Journal of Fatigue 50 (2013) 2–9

a0 ðPÞ ¼ DaðPÞ

½DK Aeq ðPÞ2

ð24Þ

½DK Beq ðPÞ2  ½DK Aeq ðPÞ2

and consequently the cyclic equivalent SIF between the successive crack fronts is approximated via A DK app eq ða; PÞ ¼ DK eq ðPÞ

rffiffiffiffiffiffiffiffiffiffiffiffi a ; a0 ðPÞ

a0 ðPÞ 6 a 6 a0 ðPÞ þ DaðPÞ:

ð25Þ

This approximation of the SIF is utilized for the re-calculation of a more accurate number of load cycles DNacc(P) that considers the changing stress field and which is therefore a good approximation to the real number of load cycles DNreal(P). It is calculated by the evaluation of the crack propagation rate:

DNacc ðPÞ ¼

Z

a0 ðPÞþDaðPÞ

1 da ðDK app eq ða; PÞÞ dN

a0 ðPÞ

da:

ð26Þ

Since the stress field has changed during crack growth the obtained incremental number of load cycles differs from the userspecified one DNlc. Moreover, the number of load cycles DNacc(P) is not constant for all points P of the crack front. In order to correct this mismatch the predicted crack extension Da(P) is replaced by

Dacorr ðPÞ ¼ DaðPÞ þ



 da B ðK eq ðPÞÞ  ½DN lc  DNacc ðPÞ: dN

ð27Þ

The modification of the crack front by this procedure leads within a few iterations to the desired crack front that is characterized by the prescribed number of load cycles. This means that the relative error is neglectable with respect to a user-specified accuracy e:

jDNlc  DNacc ðPÞj=DNlc < e:

ð28Þ

As a consequence of the changed crack extension the crack deflection has to be adapted. Here, the extended MTS-criterion Eq. (11) is evaluated with the corrected crack extension Dacorr(P). Since the stress intensity factors for all three crack opening pffiffiffiffiffiffiffiffi modes depend on p a the presented method is in principle applicable for arbitrary mixed mode problems. It requires only that the method for calculation of equivalent stress intensity factor ensures pffiffiffiffiffiffiffiffi the p a behavior even if the crack kinks. As an example the maximum energy release rate (3) which is applied within this work shows this behavior. In contrast, the energy release rate according to Irwin (2) fails this requirement. 5. Examples Two examples are chosen to demonstrate the efficiency of the presented 3D crack growth criterion. Based on the model of a penny-shaped crack the accuracy of the proposed predictor–correctorscheme is shown. Within this example mixed mode conditions are present which vary along the crack front. Second, a complex crack in a four-point-bending-specimen is investigated. Here, crack arrest is present at parts of the crack front.

Fig. 8. Example for a grown penny-shaped crack.

5.1. Penny-shaped crack The model of a penny-shaped crack as sketched in Fig. 7 is under investigation. At half of the length of a circular bar (length L = 250 mm, diameter D = 40 mm) a circular initial crack (diameter d = 2 mm) is located. In order to obtain a variable state of stress with changing mixed mode conditions along the crack front the crack surface is rotated relative to the mid cross section by the angle b = 45°. The bar consists of the material steel with a Young’s modulus of E = 210 GPa and a Poisson’s ratio of m = 0.3. It is symmetrically loaded by the cyclic tensile traction t0 with a maximum value of t0 = 100 MPa under the stress ratio R = 0. For the simulation of crack growth a Paris formulation (cf. Eq. (16)) with the material parameters m = 2.36 and C = 1.11  pffiffiffiffiffiffiffiffiffiffiffi 1011 mm (DKeq in MPa mm) is utilized for the description of the crack propagation rate. In Fig. 8 a grown crack is exemplarily illustrated. In order to verify the accuracy of the predictor–correctorscheme crack growth is simulated up to overall one million load cycles with different step sizes. Within the analysis the step sizes of DN0 = 50,000, DN0 = 100,000, DN0 = 200,000, and DN0 = 500,000 are used. Furthermore, the resulting crack growth is compared to the simulation results that are obtained with the common proce-

(a) Fig. 7. Model of penny-shaped crack.

(b)

8

W. Weber et al. / International Journal of Fatigue 50 (2013) 2–9

1.2

Table 1 Number of increments within the simulations.

1.0

a mm

0.8

Number of predictor steps

Number of corrector steps

correctors predictor

Overall

50 100 200 500

20 10 5 2

52 34 19 20

2.6 3.4 3.8 10.0

72 44 24 22

0.6

0.4

0.2

0

50

100

200

500

Δ N0 1000 Fig. 9. Crack length results explicit time integration.

1.2

1.0

0.8

a mm

DN 0 1000

0.6

0.4

5.2. Four point bending specimen

0.2

0

is mainly caused by the worse mesh quality that influences the accuracy of the SIFs calculation. In order to obtain the converged crack length of approximately 1.1 mm with the explicit time integration scheme a very small incremental step size has to be chosen. This leads directly to a high number of simulation increments which corresponds to an enormous numerical effort. In contrast, a larger step size can be specified by applying the predictor–corrector scheme. Here, corrector steps follow each prediction of a new crack front. Table 1 shows the number of predictor and corrector steps. It can be seen that the number of corrector steps depends strongly on the incremental step size. Since a larger step size causes a higher linearization error in the predictor step more corrector steps are required to ensure the prescribed accuracy of 1.25% which is a very strict specification. Moreover, a high number of corrector steps indicate a too large step size. As an optimum with respect to the quality of the results as well as the numerical complexity the simulation with the step size DN0 = 200,000 can be identified. Here, overall 24 simulation increments are required. For the simulation with the explicit time integration procedure far more than 20 increments are needed to obtain the same accuracy in the resulting crack length. Therefore, beside the guarantee of a high accuracy also a speed up in the simulation time is obtained by the predictor–corrector-scheme.

50

100

200

500

Δ N0 1000 Fig. 10. Crack length results implicit time integration.

dure of applying predictor steps only. Within all simulations the extended MTS-criterion is utilized for calculating the crack deflection in order to obtain a smooth crack path. In order to compare the simulation results the crack length due to crack growth of the upper middle point that is subjected to mode I/II-conditions only is under consideration as scetched in Fig. 8. Although there are only Modes I and II present at this point it is also influenced by mode III. Since the SIFs are obtained from a stress analysis of the whole 3D crack the stresses in front of this point are strongly influenced by the surrounding crack geometry which is additionally subjected to mode III. Fig. 9 shows the results of the simulations with predictor steps only and in Fig. 10 the results by the application of the predictor– corrector–procedure are illustrated. In Fig. 9 a strong influence of the incremental step size on the crack length is observed if an explicit time integration procedure is applied. In contrast, almost the same crack length is obtained by the implicit scheme as it can be seen in Fig. 10. The difference in the crack length between the largest step size and the other ones

As a second example the model of a four point bending specimen (cf. Fig. 11) of the material PMMA with a complex initial crack front is analyzed. It has been investigated experimentally in [24]. The initial crack front has been formed by simple mode I loading conditions. Hereby, the strong curvature of the crack front is caused by a state of residual stress in the specimen. This residual stresses with high tensile stress in the middle and the compressive stress near the free surfaces lead to an advanced crack growth in the middle of the specimen. Then, the residual stresses have been reduced by tempering the specimen just below the glass transition temperature. Afterwards, the crack propagates only in the vicinity of the free surfaces, which is a challenge for the simulation. For the simulation the crack propagation rate according Eq. (17) is utilized which has been identified experimentally for this material as [24]

da ½DK eq  DK th m ½DK eq  130:75 ¼ 0:00024 ¼C dN ½1  RK Ic  DK eq 55  DK eq

ð29Þ

pffiffiffiffiffiffiffiffiffi with da/dN in mm per loadcycle and DKeq in MPa mm. The step size is set to 10,000 and the predictor–corrector scheme is applied.

Fig. 11. Model of four point bending specimen.

W. Weber et al. / International Journal of Fatigue 50 (2013) 2–9

9

the prescribed step size within the incremental simulation procedure. Moreover, a speed up of the simulation is reached. In addition, reliable simulation results are obtained even if crack arrest is present at parts of the crack front. Here, unphysical stress concentrations are eliminated by detecting and smoothing the area of the crack front which starts growing during the simulation step. References

Fig. 12. Crack growth of four point bending specimen.

By the utilization of the predictor–corrector scheme parts of the crack front that start growing between the discrete increments can be identified. Furthermore, this area is smoothed by a cubic spline-line in order to prevent unphysical stress concentrations. Fig. 12 shows the numerically obtained crack fronts in front of the experimentally observed crack path [24]. Overall a good agreement between the numerical and the experimental obtained crack fronts is observed. Special attention within this example is also on the surface breaking points. Here, the numerical crack front is defined by the knowledge on corner singularities. An intersection angle of c  14° obtained by analysis of corner singularities has been applied in simulation. This angle is also observed for the crack fronts of the test specimen. 6. Conclusion The simulation of 3D fatigue crack propagation using an implicit time integration scheme in terms of a predictor–corrector procedure has been presented. For mixed mode problems the equivalent stress intensity factor is calculated by the presented criterion of the maximum energy release rate. The maximum tangential stress criterion under the consideration of the non-singular T-stresses is utilized for the determination of the crack deflection and a crack propagation law is evaluated in order to obtain the crack extension. Numerical examples have confirmed the necessity of the corrector steps in order to ensure a high accuracy independent from

[1] Mi Y, Aliabadi MH. Dual boundary element method for three-dimensional fracture-mechanics analysis. Eng Anal Bound Elem 1992;10:161–71. [2] Portella A, Aliabadi MH, Rooke DP. The dual boundary element method: effective implementation for crack problems. Int J Num Meth Eng 1992;33:1269–87. [3] Partheymüller P, Haas M, Kuhn G. Comparison of the basic and the discontinuity formulation of the 3D-dual boundary element method. Eng Anal Bound Elem 2000;24:777–88. [4] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques – theory and applications in engineering. Springer-Verlag; 1984. [5] Bebendorf M. Approximation of boundary element matrices. Num Math 2000;86:565–89. [6] Bebendorf M, Rjasanow S. Adaptive low-rank approximation of collocation matrices. Computing 2003;70:1–24. [7] Kolk K, Weber W, Kuhn G. Investigation of 3D crack propagation problems via fast BEM formulations. Comp Mech 2005;37:32–40. [8] Weber W, Kolk K, Kuhn G. Acceleration of 3D crack propagation simulation by the utilization of fast BEM-techniques. Eng Anal Bound Elem 2009;33:1005–15. [9] Kolk K, Kuhn G. 3D crack growth simulation. In: Proceedings of 2nd ECCM; 2001. [10] Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J Bas Eng 1963;85:519–25. [11] Weber W, Steinmann P, Kuhn G. Precise 3D crack growth simulations. Int J Fract 2008;149:175–92. [12] Leblond JB, Torlai O. The stress field near the front of an arbitrarily shaped crack in a three-dimensional elastic body. J Elast 1992;29:97–131. [13] Irwin GR. Analysis of stresses and strains near the end of a crack traversing a plate. J App Mech 1957;24:361–4. [14] Nuismer RJ. An energy release rate criterion for mixed mode fracture. Int J Fract 1975;11:245–50. [15] Lazarus V, Buchholz FG, Fulland M, Wiebesiek J. Comparison of predictions by mode II of mode III criteria on crack front twisting in three of four point bending experiments. Int J Fract 2008;153:141–51. [16] Kondratiev VA. Boundary value problems for elliptic equations in domains with conical or angular points. Trans Moscow Math Soc 1967;16:227–313. [17] Dimitrov A, Andrä A, Schnack E. Efficient computation of order and mode of corner singularities in 3D-elasticity. Int J Num Meth Eng 2001;52:805–27. [18] Kolk K. Automatische 3D-Rissfortschrittssimulation unter Berücksichtigung von 3D-Effekten und Anwendung schneller Randelementformulierungen. VDIVerlag; 2005. [19] Heyder M, Kolk K, Kuhn G. Numerical and experimental investigations of the influence of corner singularities on 3D fatigue crack propagation. Eng Fract Mech 2005;72:2095–105. [20] Plank R, Kuhn G. Fatigue crack propagation under non-proportional mixed mode loading. Eng Fract Mech 1999;62:203–29. [21] Paris PC, Erdogan F. A critical analysis of crack propagation laws. J Bas Eng 1963;85:528–34. [22] Erdogan F, Ratwani M. Fatigue and fracture of cylindrical shells containing a circumferential crack. Int J Fract Mech 1970;6:379–92. [23] Forman RG, Mettu SR. Behavior of surface and corner cracks subjected to tensile and bending loads in Ti–6Al–4V alloy. ASTM STP 1992;1131:948–53. [24] Heyder M, Kuhn G. 3D fatigue crack propagation: experimental studies. Int J Fatigue 2006;28:627–34.