Numerical stress intensity factor calculation in flawed round bars validated by crack propagation tests

Engineering Fracture Mechanics 108 (2013) 37–49

Contents lists available at SciVerse ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Numerical stress intensity factor calculation in flawed round bars validated by crack propagation tests J. Lebahn ⇑, H. Heyer, M. Sander University of Rostock, Institute of Structural Mechanics, Albert-Einstein-Str. 2, 18059 Rostock, Germany

a r t i c l e

i n f o

Keywords: Surface crack Stress intensity factor Vertex singularity Boundary layer Cyclic crack propagation

a b s t r a c t In this paper the results of crack propagation tests on round bars under cyclic tension and bending loading are presented. The tests were performed to investigate the crack geometry and its dependencies e.g. on the stress ratio and overloads. Furthermore, the crack propagation tests were used to validate the numerical stress intensity factor (SIF) calculation. This validation is necessary, due to the uncertainties in the numerical calculation of pffiffiffi three-dimensional crack problems. As shown in literature, the 1= r-singularity of the stress field is not fulfilled in general for surface cracks in the near surface domain. Due to this boundary layer effect, the classical SIF is not valid. In order to calculate SIFs at the surface point an extrapolation is used. The method’s parameters are validated by the crack propagation tests and ensured by numerical investigations concerning the singularity of the stress field and the boundary layer thickness. These investigations were performed on a CT-specimen with a straight crack front. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Efficient residual lifetime calculations of cracked structures can be performed by analytical crack propagation simulations [1]. The knowledge of the crack path and the related solution of the stress intensity factor (SIF) is a basic requirement for the simulation. Especially for semi-elliptical surface cracks it is often assumed that the crack front exhibits an elliptical geometry during propagation. Several authors [2–4] have attested the development of elliptical crack fronts by crack propagation tests. However, in the numerical crack propagation simulations performed by Hou [5], deviations from the elliptical shape are reported and explained by the crack closure effect. In order to achieve SIF solutions, numerical simulations are usually performed by use of the Finite Element Method (FEM) [6–12]. For three-dimensional crack problems, uncertainties of the calculated SIFs arise, especially in the region where the crack front intersects the structure’s surface. pffiffiffi In literature [13–18] it is well known that the 1= r-singularity of the stress field in the surface domain (boundary layer) is in general not fulfilled, whereby the conventional SIF is not valid. According to the experimental and numerical investigations of several authors [14,15,17,18] the stress field’s singularity kr depends on the Poisson’s ratio m and the angle b between the crack front and the surface of the structure. Furthermore, the thickness of the boundary layer was investigated experimentally [19] and numerically [16,18] leading to different results. From all these investigations no general conclusion on the calculation of SIFs were deduced. For semi-elliptical surface cracks in round bars many SIF solutions are available in literature [6–11]. In order to consider the boundary layer effect on the SIF calculation for the intersection point, some authors make use of an extrapolation [6] or take the SIFs from nodes behind the surface [7]. Due to the fact that all methods are not validated by experiments, in the ⇑ Corresponding author. E-mail address: [email protected] (J. Lebahn). 0013-7944/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2013.04.013

38

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

Nomenclature a crack depth, half axis of ellipse A area a/c, a/b aspect ratio Aij, Bij coordinate functions of the stress field A1, A2 coefficients of the displacement field b half axis of ellipse B1, B2 coefficients of the stress field c crack length Ci, Di coordinate functions of the displacement field D diameter E Young’s modulus ET tangent modulus F force e error value CB correction factor to calculate the SIF at point B G energy release rate J J-integral K stress intensity factor (SIF) Kmin, Kmax minimum and maximum stress intensity factor DK cyclic stress intensity factor Le element length n, N number R stress ratio: Kmin/Kmax r, u, z cylindrical coordinates R, u, h spherical coordinates rmin, rmax lower and upper bounds for regression analysis t thickness ui components of displacement vector v crack opening displacement in y-direction w width of the CT-specimen x, y, z Cartesian coordinates b angle between crack front and structure’s surface bs natural angle between crack front and structure’s surface d distance between y-axis and point B k eigenvalue, vertex singularity ku exponent of the displacement field kr singularity exponent of the stress field m Poisson’s ratio rF yield stress rij coordinates of the stress tensor xpl size of the plastic zone

present paper crack propagation tests on round bars are performed and used to investigate the development of the crack geometry and to validate the SIF calculation. Furthermore, the boundary layer effect is numerically investigated. These information are additionally used for the validation of the SIF calculation. 2. Numerical investigations of the boundary layer effect First of all, the numerical investigations concerning the singularity of the intersection point’s stress fields, the boundary layer effect and the calculation of SIFs are discussed and the related theoretical background is presented. 2.1. Special feature of surface cracks At first Benthem [13] stated the existence of a vertex singularity at the intersection point, which differs from the classical pffiffiffi 1= r -singularity. Shivakumar and Raju [15] show the existence of two singular stress fields characterized by different singularities:

rij ¼ Aij ðu; zÞ  rkr þ Bij ðu; hÞ  Rk ði; j ¼ 1; 2; 3Þ

ð1Þ

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

39

Fig. 1. Definition of the crack front geometry for investigating the boundary layer effect.

with the theoretical solution kr = 0.5. Both singularities are dominant in the surface domain of a crack front, which is pffiffiffi also called the boundary layer. The first term of Eq. (1) exhibits the classical 1= r-singularity in cylindrical coordinates, known from two-dimensional crack problems. The second term of Eq. (1) is based on the vertex singularity and can be defined in spherical coordinates, Fig. 1. Outside the boundary layer the influence of the vertex singularity is insignificant and the three-dimensional crack problem can be considered as a series of two-dimensional crack problems in planes normal to the crack front. In this domain the classical SIF KI is valid and related to the energy release rate GI and J-integral

K 2I ¼ GI  E0 ¼ J  E0 with E0 ¼



E plane stress E=ð1  m2 Þ plane strain

ð2Þ

in the case of Mode-I problems, within the LEFM. Shivakumar and Raju [15] ascertained that only the vertex singularity is dominant at the intersection point. Furthermore, they have shown that the energy release rate tends to zero for kr > 0.5 and to infinity for kr < 0.5 at the intersection point. Pook [14] and De Matos and Nowell [16] have shown that the singularity kr of the stress field depends on the Poisson’s ratio m and the angle b between the crack front and the surface of the structure. For instance for an angle b = 90° and a Poisson’s ratio m = 0.3 Benthem [13] calculated analytically a singularity of kr = 0.452, which is valid for a crack in a semi infinite plate. De Matos and Nowell [16] and Hutar and Náhlík [18] estimated the thickness of the boundary layer of middle tension specimens by numerical simulations. They have shown that the domain depends on the Poisson’s ratio. For m = 0.3 Hutar and Náhlík [18] indicate a boundary layer thickness of about 18% of the half specimen thickness.

2.2. Numerical SIF calculation of semi-elliptical surface cracks in round bars Within the scope of the validation, SIFs for semi-elliptical surface cracks in round bars were calculated from energy release rates, which was calculated numerically by the Virtual Crack Closure Technique (VCCT) [12]. According to Fig. 2 the Mode-I part of the energy release rate

GI ¼

Fz  v 2A

ð3Þ

is obtained from the z-component of the nodal force and the crack opening displacement normal to the crack plane. Herein A is the average element area. To calculate the crack opening displacement by VCCT two calculation steps are required. In the first step the nodal force is calculated. In the second step the crack opening displacement is calculated, after the crack has been extended by an incremental step Da. To reduce the computational effort, only one calculation is performed by the Modified Virtual Crack Closure Integral (MVCCI)-method [20]. Here, the crack opening displacement is taken from a neighbored node of the crack plane. However, equal element sizes ahead and behind the crack front are required for this method [21]. According to the convergence studies performed by Müller et al. [22] the element size should fulfill Da < a/20. Fig. 3a shows the quarter FE model of a round bar with a semi-elliptical surface crack. To obtain the required SIF solutions for the load cases tension and bending, six crack depths reaching from 0.05 6 a/D 6 0.7 and six aspect ratios reaching from 0.1 6 a/ b 6 1.25, Fig. 3b, were modeled. For all crack geometries the Mode-I part of the energy release rates are calculated along the crack front using the MVCCI-method. Therewith, the SIFs were determined from Eq. (2) by assuming the plain-strain condition.

40

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

(a)

(b)

Fig. 2. Definition of crack geometry for the MVCCI-method: (a) 3D mesh, definition of average element area A and (b) 2D mesh, definition of crack depth and crack length.

Fig. 3. (a) Quarter FE model of a round bar with semi-elliptical surface crack and (b) definition of the crack geometry.

2.3. Two-dimensional FE analysis on a CT-specimen The thickness of the boundary layer is determined exemplarily for a standard compact tensile specimen (CT) from the three-dimensional displacement and stress fields, Section 2.4. To find a suitable methodology, at first extensive simulations with different fine discretizations were carried out for the plane-strain problem. The specimen with a width w = 72 mm and a thickness t = 10 mm, Fig. 4, is loaded by a force F = 10 kN. A crack with a depth a = 22.5 mm and a/w = 0.31 is assumed, respectively. The FE simulations are carried out under the assumption of a linear-elastic material behavior with a Young’s

Fig. 4. CT-specimen and local mesh in the vicinity of the crack tip.

41

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

Modulus E = 210 GPa and a Poisson’s ratio m = 0.3. Due to the symmetry, only one-half of the specimen is meshed with eightnode isoparametric elements. In Fig. 4 the element arrangement around the crack tip is shown in the range 0 6 r 6 4 mm. First, the SIF was determined from energy release rate by use of the MVCCI-method assuming plane-strain conditions. Within the scope of a convergence study a value of GI = 2.047 N/mm was obtained with an element length Le = 0.1 mm around the pffiffiffiffiffiffiffiffiffi crack tip. This corresponds to a SIF KI = 687 MPa mm. From literature it is known that boundary value problems in plane elasticity theory can be solved with the method of complex stress functions [23,30,31]. For crack problems using polar coordinates, the boundary condition of traction-free crack surfaces at u ¼ p leads to a homogeneous equation system with complex coefficients and a positive, real exponent k. An infinite number of eigenvalues k = n/2 (with n = 1, 2, 3, . . .) arises from the condition that the coefficient determinant must disappear for a non-trivial solution. With it, the complex solution of the boundary value problem can be formed from the superposition of the according eigenfunctions. Therefore, with use of the Kolosov formulas, the displacement field of a plane crack under Mode-I loading is described by:

ui ðr; uÞ ¼

1 X n ðnÞ r2  ui ðuÞ ði ¼ 1; 2Þ

ð4Þ

n¼1

The first term, according to the eigenvalue k = 1/2, dominates for r ! 0 and is linked with the stress intensity factor KI. The crack opening displacements u2(r) at u = 180° can be described approximately by

u2 ðr; u ¼ 180 Þ ¼ A1  r ku þ A2  rku þ1 þ . . .

ð5Þ

with the theoretical solution ku = 0.5 for a plane crack problem. Herefrom, the known relation between the crack opening displacement and the SIF

4ð1  m2 Þ 1 u2 ðrÞ ¼ pffiffiffiffiffiffiffi  K I  r2 2p  E

ð6Þ

follows for the near-field. From the numerically approximated displacement field, the exponent ku was first determined by use of a one-term ansatz of Eq. (5). Its expression in a double-logarithmic representation becomes:

logðu2 Þ ¼ logðA1 Þ þ ku  logðrÞ

ð7Þ

The near-field was defined as the range, in which the relative error between the FE solution and the theoretical solution according to Eq. (6) is less than 1%. This was achieved for rmin 6 r 6 rmax, Table 1. To perform a convergence study the domain r 6 0.2 mm around the crack tip was meshed more finely in four steps. Table 1 shows the element sizes Le and the number n of elements in radial direction. According to the above definition of the near-field, the corresponding ranges rmin and rmax are listed in Table 1, with the calculated values ku. Taking into account the uncertainty in the definition of an upper bound for the near-field and in order to improve convergence, the second term was included in the evaluation according to Eq. (5). The lower bound rmin was held constant for each mesh and the upper bound was fixed at rmax = 1.5 mm, Table 1. This corresponds to a relative error of less than 2% between the FE solution and the theoretical solution according to the near-field, Eq. (6). The parameters A1, A2 and ku were determined by a non-linear regression with the program MATLAB. The results of ku are shown in Table 1. Concerning the exponent ku a monotonous convergence against the theoretical solution (ku = 0.5) is ascertained. The convergence of the exponent ku was significantly improved by including the second term in the regression, Fig. 5. In this manner it is possible to reach a satisfactory approximation for ku with a relatively coarse discretization surrounding the crack tip. Furthermore, it was examined in what way the asymptotic behavior of the stress field 1

1

ð2Þ ð3Þ 2 rij ðr; uÞ ¼ r2  rð1Þ ði; j ¼ 1; 2Þ ij ðuÞ þ rij ðuÞ þ r  rij ðuÞ þ . . .

ð8Þ

surrounding the crack-tip, is described by the FE solution. The stress distribution r22(r) at u = 0° can be described approximately by

r22 ðr; u ¼ 0 Þ ¼ B1  rkr þ B2  rkr þ1 þ . . .

ð9Þ

Table 1 Convergence and analysis method of the exponent ku. u2 ðrÞ ¼ A1  r ku þ A2  r ku þ1

u2 ðrÞ ¼ A1  r ku Le (mm)

n

rmin (mm)

rmax (mm)

ku

rmax (mm)

ku

0.1 0.05 0.025 0.0125 0.00625

2 4 8 16 32

0.3000 0.1500 0.0875 0.0437 0.0219

1.000 0.850 0.850 0.813 0.794

0.5137 0.5101 0.5075 0.5058 0.5045

1.5 1.5 1.5 1.5 1.5

0.5088 0.5063 0.5050 0.5030 0.5018

42

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

0.516 0.514 displacements (one-term) displacements (two-terms) stresses (two-terms)

0.512 0.51 0.508 0.506 0.504 0.502 0.5 0

5

10

15

20

25

30

35

number of elements n Fig. 5. Convergence behavior of the exponent ku depending on the number of elements.

Table 2 Convergence and analysis method of the exponent kr.

r22 ðrÞ ¼ B1  rkr

r22 ðrÞ ¼ B1  rkr þ B2  rkr þ1

Le (mm)

n

rmin (mm)

rmax (mm)

kr

rmax (mm)

kr

0.1 0.05 0.025 0.0125 0.00625

2 4 8 16 32

0.3000 0.1500 0.0875 0.0437 0.0219

1.000 0.850 0.850 0.813 0.794

0.5264 0.5118 0.5107 0.5083 0.5073

1.5 1.5 1.5 1.5 1.5

0.5357 0.4923 0.4951 0.4963 0.4982

with the theoretical solution kr = 1/2 for a planar problem. Herewith, the condition kr = ku  1 is strictly fulfilled. Thus, the known relation between the stress distribution and the SIF,

K

1

I ffi  r 2 r22 ðrÞ ¼ pffiffiffiffiffiffi 2p

ð10Þ

corresponding to the first term of Eq. (9), arises for the near-field. The evaluation was again performed in a double-logarithmic representation

logðr22 Þ ¼ logðB1 Þ þ kr  log

ð11Þ

Also that area was defined as the near field, in which the relative error between the FE solution and the theoretical solution, Eq. (10), is smaller than 1%. Contrary to the approximation of the displacement field, this requirement for the quality of the FE solution could not be realized for all investigated element sizes, due to the singularity of the stress field. That is the reason why the ascertained borders from the approximation of the displacement field, Table 1, were used for the evaluation by Eq. (11), Table 2. The convergence behavior of the exponent kr is much worse than for the exponent ku. The values ascertained within the scope of the convergence study for ku fulfill the condition ku > 0.5, Table 1. According to the requirement kr = ku  1 the condition kr > 0.5 must be fulfilled therefore. If only the singular term of Eq. (9) is used, this condition is not kept for all element sizes, Table 2. To improve the convergence behavior and the approximation of the stress field, the second term of Eq. (9) was included in the regression. The bounds rmin and rmax were both adopted for the evaluation from the analysis of the displacement field, Tables 1 and 2. The parameters B1, B2 and kr were determined by a non-linear regression with the program MATLAB. The obtained results of kr are shown in Table 2. With the exception of the coarse mesh with Le = 0.1 mm, a monotonous convergence against the theoretical solution is ascertained. A comparison of the convergence behavior of the exponent ku, calculated by the different methods under consideration of the relation kr = ku = 1, is shown in Fig. 5. 2.4. Three-dimensional FE analysis on a CT-specimen The three-dimensional investigations of the exponent ku and the size of the boundary layer resulting from the vertex singularity, are carried out by use of Eq. (5). The investigations are performed on a CT-specimen with a straight crack front,

43

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

G1 [N/mm]

2.5

2.0 plane strain 3D

1.5

1.0 0.0

1.0

2.0

3.0

4.0

5.0

z [mm] Fig. 6. Energy release rate along the crack front.

Fig. 1. Due to the symmetry, only one quarter of the specimen is meshed with 20-node isoparametric elements. The vicinity of the crack front was meshed, in every plane z = constant, with a radial element pattern characterized by a length Le = 0.0125 mm. Near the free surface a fine mesh with an element size in through-thickness direction of 0.0125 mm is used. Towards the middle plane of the specimen (z/t = 0) the element size is gradually increased in through-thickness direction up to a size of 0.05 mm. A plane-strain state was generated at first by preventing the displacements normal to the surface with z/t = 0.5 in order to make it comparable to the results to the 2D analysis. Therefore, the resulting energy release rate G, Fig. 6, and the exponent ku are independent of the z-coordinate. The change of the energy release rate along the crack front for the 3D problem is also shown in Fig. 6. It is influenced by the test geometry, the load and the vertex singularity at the free surface. Based on elasticity-theoretical considerations, for Poisson’s ratio m > 0 and ku > 1/2, the energy release rate at the free surface must be zero [15]. The distribution of the energy release rate shows a strong gradient for z ? 5 mm, Fig. 6, according to the theoretical trend. Shivakumar and Raju [15] suggested that the displacement field in the vicinity of the crack front can be described by a superposition of two terms

ui ¼ C i ðu; zÞ  r 0:5 þ Di ðu; hÞ  Rk

ði ¼ 1; 2; 3Þ

ð12Þ

outgoing from the known cylindrical singularity of the stress field inside the body and an unknown vertex singularity in spherical coordinates with a constant exponent k. The displacement field of a crack’s near-field is affected by the vertex singularity dependent on the distance R, Eq. (12). The vertex singularity itself depends on the Poisson’s ratio. The FE analyses are evaluated concerning the crack opening displacements in y-direction according to Eq. (5) in a cylindrical coordinate system. In the middle plane of the specimen the cylindrical singularity dominates with ku  0.5 and at the free surface the vertex singularity dominates with ku = k. There exists a boundary layer, influenced by the vertex singularity, in which the cylindrical singularity changes significantly and the SIF is not valid. The exponent ku is determined for different planes z/t, Fig. 7. In each case the computed displacements in a range 0.04 mm 6 r 6 1 mm, are described by a two-term ansatz, Eq. (5), using the program MATLAB. In the middle plane of the model ku = 0.5033 is achieved with an adequate match of the theoretical value. At the free surface a value of 0.5503 arises for the exponent ku. The size of the boundary layer, identified as the domain, in which ku deviates more than 1% of the cylindrical singularity, corresponds to about 19% of the half crack front length and agrees well with the investigations of Hutar and Náhlík [18].

0.56

Boundary Layer

0.55 0.54 0.53 0.52 0.51 0.5 0.3

0.32 0.34 0.36 0.38

0.4

0.42 0.44 0.46 0.48

0.5

Fig. 7. Influence of the vertex singularity on the exponent ku along the crack front.

44

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

3. Crack propagation experiments on round bars The crack propagation tests were performed on round bars with a diameter of 20 mm. All bars consist of the high strength steel 34CrNiMo6. 3.1. Test schedule The specimens were loaded by cyclic tension or bending. Different stress ratios were investigated, Table 3. Due to the 4point-bending device only positive stress ratios were tested for the bending load case. For both load cases constant amplitude sequences were adopted. The stress amplitudes decrease from block to block so that the crack propagation occurs in the Paris-regime (105 mm/cycle–104 mm/cycle). Between the blocks, overloads with an overload ratio OL = 2 or a red dyestuff were inserted to generate beach marks. After the tests, the beach marks were measured under a traveling microscope to investigate the development of the crack geometry. To perform the test procedure, a continuous crack growth measurement was needed. Therefore, the DC potential drop method was used. The required calibration curve was obtained from pretests. To initiate a crack, the specimens were notched by bore holes, Fig. 8a, or milled ellipses, Fig. 8b. The analysis of the fracture surfaces contains the measurement of the beach marks, the screening of their elliptical shape and the allocation to load cycles. Measurement values are the crack depth a, the crack length 2c, the distance D and the angle b between the crack front and the surface at point B, Fig. 3b. A visual screening of the elliptical shape was used. Therefore, each beach mark was fitted by an ellipse, Fig. 8. For the first crack propagation tests some fracture surfaces exhibit beach marks, which deviate from the elliptical shape, Fig. 8b. To investigate this, the beach marks were separated into elliptical and non-elliptical ones and plotted over the maximum SIFs of the overloads. For the calculation of the SIFs, the geometry function from Raju and Newman [10] was used. It could be observed, not shown here, that the beach marks differ from the elliptical shape, if the maximum SIF at point A exceed about 1.900 N/mm3/2. To investigate the reasons of this limit, elastic–plastic FE simulations were performed on a natural crack front with a = 5 mm and a/c = 0.89, resulting from the crack propagation tests for the bending load case, Fig. 9. Isotropic hardening with an initial yield strength rF = 825 MPa and a linear flow curve with a tangent modulus of ET = 825 MPa was used as material model. The nominal stresses were increased step wise till a maximum value of 800 MPa was reached. For each step the diameter xpl of the plastic zone in front of the crack at the points A and B were measured. Additionally, the geometry factors for both points were calculated in a linear-elastic FE simulation from energy release rate. It could be observed that the plastic zone xpl at the surface is bigger than in the interior, due to the different stress states and the higher SIF. Thus, the retardation of the crack propagation is larger at the surface of the crack front than in the interior

Table 3 Overview of the performed crack propagation experiments. Stress ratio

1 0.1 0.5

Overload ratio = 2

Overload ratio = 1

Bending

Tension

Bending

Tension

– 19 5

4 5 –

– 2 5

– 4 –

Fig. 8. Fracture surfaces with beach marks and ellipse. (a) Elliptical shaped beach marks and (b) beach marks differ from the elliptical shape.

45

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

Aspect ratio a/c

2.0 Bending Tension

1.5

1.0

0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Normalized crack depth a/D Fig. 9. Development of the aspect ratio over the crack depth of all beach marks.

after an overload. If the plastic zone and the overload are too big, only crack propagation in the interior is possible, leading to non-elliptical crack geometries, Fig. 8b, which is in accordance to the numerical investigations from Hou [5]. According to the theory of LEFM the plastic zone size has to be small compared to the crack size and the ligament. Related to the investigated crack geometry it is demanded that xpl < 0.5 mm. Therefore, the maximum SIF at point A of an overload is limited to Kmax < 1.700 N/mm3/2. For all further tests the limit of Kmax was applied and only elliptical crack geometries were observed. If the assumptions of LEFM are fulfilled, the cracks propagate with an elliptical geometry, in agreement to [2–4]. 3.2. Development of the crack geometry The investigations of the crack geometry were performed on 290 elliptical beach marks. For tension and bending, Fig. 9 shows the development of the aspect ratio a/c over the crack depth. The regression functions for bending and tension are also plotted. As can be seen, there is approximately a quadratic dependence between the aspect ratio and the crack depth for both load cases. The deviations between tension and bending decrease with increasing crack depth. To identify the dependencies of the crack geometry a sensitivity analysis was performed using the software Visual-XSel 12.0 [24]. The input values are the stress ratio R, the overload ratio OL, the load case and the crack depth a/D of all beach marks, while the output value is the aspect ratio a/c. To derive the sensitivities a regression analysis was performed by use of a quadratic model with interaction effects. The reduced model, containing only the significant factors, explains 93% of the test data. The relative effects of the significant input values on the aspect ratio are listed in Table 4. As can be seen, the crack depth in terms of a/D dominates the

Table 4 Sensitivities of the crack geometry from the investigated factors. Output value

Relative effect (%) a/c

b

Input value

a/D

Load case

a/D

Load case

OL

Main effect (linear) Main effect (quadratic) Interaction effect

38.2 45.5 5.3

5.7 0.0 5.3

52.7 29.9 0.0

8.5 0.0 0.0

8.9 0.0 0.0

140

Intersection angle β [°]

130 120 110 100 90 Bending Tension Mean value bending Mean value tension

80 70 60 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Normalized crack depth a/D Fig. 10. Intersection angles of all beach marks and their mean value.

46

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

aspect ratio. There is also a significant dependence of the load case and an interaction effect between the crack depth and the load case, Fig. 9. The interaction effect becomes obvious by the decreasing distance of the regression functions for bending and tension with increasing crack depth. All other input values are insignificant for the aspect ratio. Furthermore, the intersection angles b between the crack front and the surface at point B are investigated, Fig. 10. As can be seen by the regression functions, the intersection angles depend on the crack depth and the load case. In a regression analysis, according to those from the aspect ratio, the significant input factors were identified, Table 4. The intersection angle mainly depends on the crack depth. The load case and the overload ratio are also significant. No interaction effect between the crack depth and the load case on the intersection angle was attested, as it is also obvious from the parallel regression functions, Fig. 10. pffiffiffi According to the investigations of Heyder et al. [17] the crack geometry of natural cracks fulfills the 1= r-singularity of the stress field along the whole crack front. If this yields, in cracked round bars the intersection angle bs which fulfills the pffiffiffi 1= r -singularity of the stress field at point B depends on the load case and the crack depth. This would imply that Pook’s [14] empirical equation

bs ¼ arctan



v 2

 ð13Þ

m

is not valid for cracks in structures with curved surfaces. In this equation bs only depends on Poisson’s ratio. For the considered material bs = 100°, which equals the mean value of the intersection angles for the bending load case but differs significantly for the tension load case, Fig. 10. 4. Validation of the SIF calculation The validation of the SIF calculation is based on the experimental development of the crack geometry, quantified by the aspect ratio, which is characterized by the regression functions of the beach marks, Fig. 9. Due to a larger data base, only the bending load case was used. As mentioned, an extrapolation is used to avoid the boundary layer effect for the calculation of SIFs at point B. To validate the parameters of the extrapolation, the numerical investigations of the boundary layer thickness, Section 2.4, are considered. 4.1. Extrapolation and correction of the SIFs at point B To extrapolate the SIFs at point B a polynomial regression function is fitted through the SIFs in the regression domain of the crack front, Fig. 11. This polynomial function is extrapolated to the surface in order to calculate the SIF at point B. The

Boundary layer

from

KB,ex

to Regression domain

640

3/2

K in [N/mm ] I

620

·

600 580

Data points used for the extrapolation Data points not used for the extrapolation Regression function for the extrapolation Extrapolated SIF at point B

560 540 520

0

0.2

0.4

0.6

KA

0.8

1

s/s

max

Fig. 11. SIF calculation at the surface point by use of an extrapolation.

Table 5 Investigated parameter combinations of the validation of the SIF calculation. Parameter

Start value

End value

Step width

Number

CB Order From To Combinations

0.9 1 0 0.3 –

1.0 4 0.3 1.0 –

0.01 1 0.025 0.1 –

11 4 13 8 4576

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

47

parameters of the extrapolation are the order of the regression function as well as the limits of the regression domain. The range of the parameters and the step widths are listed in Table 5. It could be observed that all parameter combinations lead to aspect ratios in the analytical crack propagation simulation, which are smaller than those from the crack propagation tests. This trend was also attested for surface cracks in literature [2,28] and was explained by different crack closure effects between the inner and the outer surface domain of the crack front. In order to consider this in the analytical crack propagation simulation Newman and Raju [28] introduced a correction function

K B ¼ K B 

8 2 4 > < 0:9 þ 0:2R  0:1R > :

R>0 ð14Þ

for

0:9

R60

which depends on R. This function is based on crack propagation tests of surface cracks in aluminum plates. Müller et al. [2] also performed crack propagation tests on surface cracks in aluminum and steel plates. They could not identify a significant dependence between the stress ratio and the crack geometry as it was also attested by the sensitivity analysis in the present study. Therefore, a constant correction factor CB is used to calculate the SIF at point B. This factor has been validated as shown in Table 5. 4.2. Validation concept The validation concept consists in a comparison of the experimental crack geometry with those from the analytical crack propagation simulation. The crack propagation simulations are performed by the software NASGRO 6.02 [25], which contains a data base of SIF solutions of different crack problems. It is also possible to enter additional SIF solutions in NASGRO. In the case of semi-elliptical surface cracks the SIFs of the deepest and the surface point have to be provided to NASGRO as 2D tables. Herewith, the SIFs for arbitrary crack geometries are calculated by 2D interpolation routines [25]. In NASGRO the Forman/Mettu-equation [26] is used to describe the crack propagation curve, while the crack closure effect is considered by Newman’s crack closure function [27]. Because the same crack closure effect is calculated for point A and B, the development of the crack geometry in the analytical crack propagation simulation is only dependent on the SIF solution. Thus, it follows a characteristic a/D–a/c-curve for each SIF solution, independent of the stress ratio. If the initial aspect ratio differs from this curve, it converges to it with increasing crack depth. The SIFs of the different semi-elliptical surface cracks in a round bar were calculated along the crack fronts, Section 2.2. Herewith, the SIFs at point B are determined by an extrapolation and correction. Each combination of the parameters, used for the extrapolation, leads to different SIFs at point B and therefore to different SIF solutions and different characteristic a/ D–a/c-curves. By comparison these curves with the a/D–a/c-curve of the crack propagation tests (regression function, Fig. 9), the optimum parameter combination of the extrapolation is found. The validation schedule is as follows: for a parameter combination i the SIFs of all crack geometries of the round bar at point A and B are inserted in NASGRO; the SIFs at point B are calculated by an extrapolation and a correction depending on the parameter combination i; an analytical crack propagation simulation is performed with the initial crack geometry a/ D0 = 0.05 and a/b0 = 0.96 (The initial aspect ratio was calculated with the regression function fbending for the bending load case, Fig. 9.); finally from this simulation the characteristic a/D–a/c-curve is obtained as N discrete pairs of value ½Da ðiÞ ; acjðiÞ . j Thereby, an error

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u  ðiÞ 2 N  ðiÞ u1 X a a ei ¼ t   f bending N j¼1 c j Dj

ð15Þ

is calculated for the parameter combination i by a comparison of the characteristic a/D–a/c-curve with the regression function f bending of the crack propagation tests. In this way, an error ei is calculated for every parameter combination i. 4.3. Results of the validation According to the numerical investigations of the boundary layer, Section 2.4, it is demanded that the distance between the surface and the regression domain is at least 19% of the half crack front length. With respect to this condition and a minimum error ei between simulation and experiment, the optimal parameter combination is as follows: CB = 0.97, Order = 2, From = 0.2 s/smax, To = 1.0 s/smax. A comparison of the development of the aspect ratios from simulation and experiment is shown in Fig. 12. The parameters of the extrapolation are in agreement with those from Shin and Cai [6]. It should be noted that the correction factor CB = 0.97 is bigger than the correction factor from Newman and Raju for |R| < 0.6. For instance, if R = 0.1 the correction function is nearly equal to 0.9. One reason for the different values of the correction factors could be the different materials used in the crack propagation tests. According to Silva [29] the crack closure effect is influenced by the cyclic hardening behavior of the material. The material 34CrNiMo6 of the round bars shows cyclic softening whereas the aluminum alloy used by Newman and Raju is cyclically neutral [29].

48

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

Aspect ratio a/c

1.6 Beach marks

1.4

Simulation 1.2

Data fit beach marks

1.0 0.8 0.6 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Normalized crack depth a/D Fig. 12. Development of the crack geometry of crack propagation tests and simulation.

As mentioned before, many SIF solutions are available in literature for semi-elliptical surface cracks in round bars. For the bending load case and an aspect ratio a/c = 1, some SIF solutions are plotted over the crack depth a/D at point A, Fig. 13a and at point B, Fig. 13b. The SIFs at point A from different researchers are in very good agreement for a/D < 0.5, which is the significant domain for cyclic crack propagation simulations. Contrary to this, there is a high scatter between the SIFs at point B due to the uncertainties mentioned before. Although Shin and Cai [6] use a similar extrapolation their SIFs at point B differ to those from the present study.

(a)

(b)

Fig. 13. Comparison of geometry functions, a/c = 1, bending load case: (a) at point A and (b) at point B.

J. Lebahn et al. / Engineering Fracture Mechanics 108 (2013) 37–49

49

5. Conclusions The numerical investigations of the boundary layer effect allow to get the following conclusions: 1. The convergence of ku followed from the stress field is much worse as for the exponent ku obtained from the displacement field. 2. In both cases the convergence behavior is improved by the two-term ansatz. 3. It is possible to deduce the thickness of the boundary layer from the development of ku over the crack front. For the investigated straight crack front the thickness of the boundary layer is 19% of the half crack front length. 4. For a straight crack front the energy release rate tends to zero at the surface. The crack propagation tests and the validation can be summarized as follows: 1. If the assumptions of the LEFM are fulfilled, the cracks propagate with an elliptical crack geometry. 2. The naturally arising aspect ratio as well as the interaction angle are almost independent on the stress ratio and the overload ratio, but dependent on the crack depth and the load case. 3. For the calculation of the SIFs at point B an extrapolation and correction are used. The parameters were validated by the crack propagation tests under consideration of the numerically calculated boundary layer thickness. 4. A R-independent correction factor is used for the determination of the SIFs at point B to improve the accordance between simulation and test.

References [1] Sander M, Richard HA, Lebahn J, Wirxel M. Fracture mechanical investigations on wheelset axles. In: Proceedings of 16th international wheelset congress, Cape Town, South Africa; 2010. [2] Müller HM, Müller S, Munz D, Newman J. Extension of surface cracks during cyclic loading. Frac. Mech.: 17th vol., ASTM STP 905, Philadelphia; 1986. p. 625–43. [3] Lin XB, Smith RA. Shape growth simulation of surface cracks in tension fatigued round bars. Int J Fatigue 1997;19:461–9. [4] Carpinteri A, Brighenti R. Part-through cracks in round bars under cyclic combined axial and bending loading. Int J Fatigue 1996;18:33. [5] Hou C-Y. Simultaneous simulation of closure behavior and shape development of fatigue surface cracks. Int J Fatigue 2008;30:1036–46. [6] Shin CS, Cai CQ. Experimental and finite element analysis on stress intensity factors of an elliptical surface crack in a circular shaft under tension and bending. Int J Frac 2004;129:239–64. [7] Madia M, Beretta S, Schödel M, Zerbst U, Luke M, Varfolomeev I. Stress intensity factor solutions for cracks in railway axles. Engng Frac Mech 2011;78:764–92. [8] Levan A, Royer J. Part-circular surface cracks in round bars under tension, bending and twisting. Int J Frac 1993;61:71–99. [9] Carpinteri A. Elliptical-arc surface cracks in round bars. Fatigue Frac Engng Mater Struc 1992;15:1141–53. [10] Raju IS, Newman JC. Stress-intensity factors for circumferential surface cracks in pipes and rods under tension and bending loads. Frac. Mech. 17th vol., ASTM STP 905, Philadelphia; 1986. p. 789–805. [11] Shiratori M, Miyoshi T, Sakai T, Zhang GR. Analysis of stress intensity factors for surface cracks subjected to arbitrarily distributed surface stresses. Trans Jpn Soc Mech Engng 1986:25. [12] Shivakumar KN, Tan PW, Newman JC. A virtual crack-closure technique for calculating stress intensity factors for cracked three dimensional bodies. Int J Frac 1988;36:43–50. [13] Benthem JP. State of stress at the vertex of a quarter-infinite crack in a half-space. Int J Solids 1977;13:479–92. [14] Pook LP. Some implications of corner point singularities. Engng Frac Mech 1994;48:367–78. [15] Shivakumar KN, Raju IS. Treatment of singularities in cracked bodies. Int J Frac 1990;45:159–78. [16] De Matos PFP, Nowell D. The influence of the poisson’s ratio and corner point singularities in three-dimensional plasticity-induced fatigue crack closure. Int J Fatigue 2008;30:1930–43. [17] Heyder M, Kolk K, Kuhn G. Numerical and experimental investigations of the influence of corner singularities on 3D fatigue crack propagation. Engng Frac Mech 2005;72:2095–105. [18] Hutar P, Náhlík L. Fatigue crack shape prediction based on vertex singularity. Appl Comput Mech 2008;2:45–52. [19] Smith CW, Epstein JS, Olaosebikan O. Boundary layer effects in cracked bodies: an engineering assessment. Frac. Mech., 17th vol., ASTM STP 905, Philadelphia; 1986. p. 775–88. [20] Buchholz F-G. Finite element analysis of a 3D mixed-mode fracture problem by virtual crack closure integral methods. In: Krishna Murthy AV, Buchholz F-G-, editors. Fracture mechanics. Bangalore: Interline Publishing; 1994. p. 7–12. [21] MSC Marc – Theory and user information, reference manual, vol. A, MSC Software, R 2007. [22] Müller T, Heyer H, Sander M. Numerische Bestimmung von Spannungsintensitätsfaktoren mittels automatisierter Netzgenerierung an rissbehafteten Wellen unter Berücksichtigung technologisch bedingter Eigenspannungen. DVM-Bericht 2011;243:55–64. [23] Gross D, Seelig T. Fracture mechanics: with an introduction to micromechanics. Berlin Heidelberg: Springer-Verlag; 2011. [24] Ronniger C. Design of experiments & statistics. Visual-XSel 12.0, CRGRAPH, Munich; 2012. [25] NASGRO – Fracture mechanics and fatigue crack growth analysis software, reference manual, Version 6.0; 2009. [26] Forman RG, Mettu SR. Behavior of surface and corner cracks subjected to tensile and bending loads in Ti-6Al-4V alloy. Frac. Mechanics 22, vol. I, ASTM STP 1131; 1992. p. 519–46. [27] Newman JC. A crack opening stress equation for fatigue crack growth. Int J Frac 1984;24(3):131–5. [28] Newman Jr JC, Raju IS. Prediction of fatigue crack-growth patterns and lives in three-dimensional cracked bodies. Adv Frac Res, Sixth Int Conf Frac 1984;3:1597–608. [29] Silva FS. Crack closure inadequacy at negative stress ratios. Int J Fatigue 2004;26:241–52. [30] Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity. Leyden: Noordhoff International Publishers; 1975. [31] Williams ML. The complex-variable approach to stress singularities – II. ASME J Appl Mech 1956;23:477–8.