- Email: [email protected]
Simulation of surface crack shape evolution using the finite element technique and considering the crack closure effects
International Journal of Fatigue 33 (2011) 719–726
Contents lists available at ScienceDirect
International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue
Simulation of surface crack shape evolution using the finite element technique and considering the crack closure effects Chien-Yuan Hou Department of Civil Engineering, National Chi Nan University, Puli, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 27 August 2010 Received in revised form 16 November 2010 Accepted 30 November 2010 Available online 4 December 2010 Keywords: Surface cracks Finite elements Crack shape Crack closure
a b s t r a c t Shape development and closure analysis of surface cracks using the finite element method have been employed by researchers without considering interactions of the two phenomena. This study uses the free-front finite element technique to account for such interactions. Mesh refinement effects on the crack shape development are studied by concluding that the element size ahead of the crack front is not affecting the simulated crack shape. Newman and Raju assumed different Paris equation coefficients for examining the surface crack growth. This assumption is here analysed by applying the crack closure concept using the results of the free-front technique. The results are close to those determined by employing the above assumption for the equation coefficients. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction One of the most important issues related to surface cracks subjected to fatigue loads is the shape development. Newman and Raju [1] were the pioneers with regards to this problem. They assumed that an initial semi-elliptical surface crack maintains a semi-elliptical shape during the whole propagation, and estimated the growth increments using the Paris law simply at the free surface and the deepest penetration point of the crack, i.e., points C and A in Fig. 1, respectively. Connolly and Collins [2] used the alternative current field measurement technique to analyse the shape change of elliptical surface cracks. McFadyen et al. [3] measured the shape evolution of surface cracks with various aspect ratio values and used the Newman and Raju’s approach to predict the crack aspect ratio change. Shin [4] observed the crack growth rate and shape development of corner cracks and calibrated the stress intensity factor of such cracks. The surface crack growth in round bars has been studied by some researchers [5–8]. Song and Shieh measured the crack closure behavior of surface cracks and found that the closure slowed down the crack growth rate at the surface [9]. Heyder and Kuhn observed the 3-D crack propagation in transparent PMMA specimens [10]. Since the shape of semi-elliptical surface cracks during propagation unnecessarily remains elliptical, Smith et al. [11] and Lin and Smith [12,13] developed a multiple point scheme which monitored not only the crack front growth at points A and C, but the growth at several front points between A and C was also moniE-mail address: [email protected] 0142-1123/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijfatigue.2010.11.022
tored. In the scheme, the shape assumption, i.e., the mathematically-defined ellipse, was no longer required, and different finite element models were constructed when a new crack front was given. Hence the finite element model was continuously changing due to the crack front evolution, and a computer code to generate the new model was required. Since the phenomenon of plasticity-induced crack closure was observed, it has been used to explain many metal fatigue phenomena such as mean stress effects, crack growth retardation, anomalous growth behavior of short cracks, etc. The 2-D simulation of crack closure by the finite element technique has achieved interesting results during the last twenty years of the 20th century. The 3-D crack closure analysis has been rare in that period probably because of computation facilities not sophisticated enough to perform the required complex computations. In recent years, the difficulty has been overcome and some results have been reported. Chermahini et al. [14] extended the 2-D closure analysis to semicircular and semi-elliptical cracks, but the mesh density of the finite element model was coarse due to the limitation of computation facilities. Zhang and Bowen [15] constructed a simple semi-circular crack finite element model to calculate the closure behavior when the crack was subjected to a single overload. Skinner and Daniewicz [16] performed finite element calculations for the closure level of semi-elliptical cracks. Hou [17] carried out calculations related on both semi-circular and -elliptical cracks at various loading levels and R ratios. Simandjuntak et al. studied the closure behavior of corner cracks emanating from a hole, and compared the calculated results to the observed data [18]. Kim et al. [19] investigated the stabilization behavior of the crack
720
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
opening level, the effect of mesh size, and the effect of initial crack length in the 3-D closure analysis. de Matos and Nowell [20] studied the crack front geometry influenced by the stress singularity near the free surface, and the crack opening stresses were calculated using the modified crack shape. Obviously, the crack shape development and the crack closure are two factors interacting each other when a 3-D crack propagates. For a surface crack, the crack closure effect is more pronounced at the free surface due to the plane stress state. As a result, the crack growth rate at the free surface is reduced and a front-turning-inward phenomenon is usually observed (see Fig. 1). The crack front turning inward near the free surface increases the stress intensity factor at this region and subsequently causes different crack closure levels. Unfortunately, none of the above mentioned 3-D finite element crack shape [1,11–13] or crack closure [14–20] studies, or other similar studies reported in literature take the interactions of the two phenomena into consideration. In a conventional 3-D finite element crack closure analysis, the crack shape is primarily considered in the stage of the finite element model construction. As the model is determined, the initial and the subsequent crack shapes cannot be altered due to the crack front nodes are released simultaneously, that is, the evolving crack shape must follow the prescribed mesh arcs. If the technique is to be used for simultaneous simulation of the crack closure and crack shape development, it is obvious that crack front must be free to move without restraining to the prescribed mesh arcs in the crack closure analysis. To solve the difficulty, Hou [21] developed a finite element technique without releasing the crack front nodes simultaneously. This strategy combines the multiple point scheme [11–13] and the conventional crack closure analysis and, hence, allows the crack front free to move without restraining so that the evolving crack shape can be simulated and the crack opening stresses can also be evaluated concurrently. The simulated crack shapes and crack opening levels determined by this new technique seemed to be promising. As a subsequent study, this paper uses the above free-front technique further addressing the effects of the mesh refinement and employing different meshes by considering various species of arcs on the simulated crack shape. In addition, the Newman and Raju’s assumption related to the surface crack shape development is assessed using the results of the above technique. 2. The free-front finite element technique for the crack closure simulation Fig. 1 schematically shows the crack-plane meshes of a 3-D finite element crack closure analysis. The initial crack front in the
conventional crack closure analysis is coincident to a mesh arc, and the front nodes are released simultaneously during a loading cycle. As a consequence, the subsequent crack front is also coincident to one of the mesh arcs: see fronts C0 to Cj in Fig. 1. Hence, once the mesh plan is determined, the crack shape is restrained for the following crack closure analysis, and the crack shape change cannot be captured. Hou developed a finite element technique so that both the crack closure and shape can be simultaneously simulated [21]. In such a technique, two finite element models are constructed: one is used to calculate the stress intensity factor at each front node and the other is used for the crack closure analysis. Since the first model is used for the stress intensity factors, the special elements of quartered point are deployed along the front so that the stress singularity can be captured and hence the model is analysed elastically. Elastic–plastic analysis is performed by using the second model, in a way similar to the conventional closure analysis, and the opening stresses (Sop) at the front nodes are determined. As the stress intensity factors and opening stresses are obtained, the effective stress intensity factor can be calculated at all front nodes. Along the crack front, the node with the maximum effective stress intensity factor is determined and the crack increment at this node is set to one radial element size for the closure analysis of the next loading cycle. The ratios of other front nodal increments to the increment of the maximum-effective-stress-intensity node are then calculated using the modified Paris law:
ðDK eff Þni Dai ¼ Damax ðDK eff Þnmax
ð1Þ
where Dai and Damax are the crack increments at the front node i and at the node with the maximum increment, respectively; DKeff is the effective stress intensity range; and n is the power exponent of the Paris law. The increment of other front nodes must be lower than one radial element size, and is accumulated with the increment of the previous loading cycle. Once the accumulated increment at a front node is larger than one element size, the node is released for the closure analysis of the next loading cycle. Hence, in the proposed technique, the front nodes are not released simultaneously in one loading cycle and the crack front is not restrained to the prescribed mesh arcs during the analysis. Thus the crack shape evolution can be simulated. For example, using such a technique, a crack front can be formed by the solid circular symbols, as is shown in Fig. 1. By employing the proposed technique, the calculated crack face profile is rugged as compared to that of the conventional closure analysis, but the simulated opening stresses and crack shape are promising. 3. Effects of mesh plan on the simulated crack shape
C0: initial crack front Cj: crack front of the j th load cycle
Cj C0
free surface
a
A: point of deepest penetration
c
C front turning inward
Fig. 1. Schematic illustration of the restrained-front and free-front techniques for the finite element crack closure analysis.
3.1. The crack closure model and loading conditions The issue of the mesh refinement has received numerous attentions in the crack closure analysis since different crack opening levels can be obtained from finite element models with various mesh sizes. It is also interesting to assess whether the simulated crack shape is influenced by various models with different mesh sizes. A semi-circular surface crack with a radius of 5.25 mm is here examined. The crack-plane geometry of the analysed specimen and the finite element meshes of the specimen are shown in Fig. 2. Because of the symmetrical conditions of the specimen, the model is constructed simply for a quarter of the specimen. Eight-node hexahedral elements are used in the model for the crack closure analysis. Region A (bordered by thick lines in Fig. 2) is a fine mesh area, and element sizes of 0.005, 0.01, and
721
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
25 mm
A
initial crack front
60 mm
(a)
(b)
Fig. 2. (a) The dimensions of the crack plane and the mesh plan, and (b) the finite element meshes of the analysed specimen.
0.02 mm are respectively placed along the radial direction of three different models. Along the circumferential direction of such three models, 32 elements are deployed. The material is assume to be elastic-perfectly plastic with a plastic flow stress r0 = 360 MPa. The Young’s modulus and Poisson’ ratio are assumed to be 200 GPa and 0.3, respectively. Constant amplitude loads perpendicular to the crack plane with R = 0, Smax/r0 = 0.33 and 0.5 are applied (Smax = maximum stress). Skinner and Daniewicz [16] suggested that five elements embedded in the crack front plastic zone are likely to be required for the convergence of the stabilized crack opening stresses. The above radial element sizes are chosen in this way because the model with 0.02mm elements produces a plastic zone close to the deepest penetration point containing approximately five elements when Smax/ r0 = 0.33. 3.2. Crack face node opening and closing criteria The loading-direction displacements of crack face nodes are monitored in each unloading step. When any crack face node displacement is found to be negative indicating that the node is closed, the node displacement along the loading direction is fixed to zero for the following unloading steps. During the loading phase, the loading-direction reaction forces of the fixed crack face nodes are monitored. If a face node is still in contact, the direction of the reaction force should be the same as that of the external load. Once the direction of the reaction force becomes opposite to that of the external load, the fixed node is released for the following loading steps, that is, the crack face at the node is opened. 3.3. The simulated results of the conventional crack closure analysis Fig. 3 shows the calculated relationship between Sop/Smax and crack increment using the conventional restrained front approach. Note that, for both loading conditions (Smax/r0 = 0.33 and 0.5), the calculated ratio Sop/Smax for the 0.01- and 0.005-mm models shows an increasing and stabilized trend at the free surface, whereas such a ratio at the deepest penetration point shows an increasing, then decreasing and finally stable trend. Skinner and Daniewicz [16] performed the crack closure analysis of semi-circular cracks with the restrained crack front technique. The applied loads were Smax/r0 = 0.27 and R = 0.1. Since these conditions are close to Smax/r0 = 0.33 and R = 0, their results are used in Fig. 3a for comparison. In their study, three different
initial crack dimensions with the ratio of the crack depth to the specimen thickness a/t = 0.34, 0.53, 0.71 were considered and ten loading cycles were carried out. This number of loading cycles allows the cracks to travel to a maximum crack increment approximately 0.05 mm, which is much smaller than the crack increment range considered in the present study. Their results show that along the direction of the deepest penetration point the values of Sop/Smax become stable in ten loading cycles. Hence, it is assumed that the calculated Sop/Smax values remain stable to a larger crack increment and these values are illustrated as the hollow symbols in Fig. 3a. Note that these symbols are plotted at an arbitrarily chosen crack increment. It is obvious that the present results are slight smaller but close to those reported values. Kim et al. [19] performed the crack closure analysis for Smax/ r0 = 0.54 and R = 0. The reported Sop/Smax values are illustrated in Fig. 3b for comparison. At the free surface, results of Kim et al. reach the stable value earlier than those of the present study but both results show close stable values. Along the direction of the deepest penetration point, both studies show the similar trend with different stable values. These differences are mainly caused by (1) different initial crack dimensions, (2) different stress–strain relationships, and (3) different node release scheme used in the calculations. The plastic zone sizes at the free surface and the deepest penetration point were estimated using Irwin’s approximation, thus the number of elements embedded in these two plastic zones can be calculated. The relationship between the calculated stabilized Sop/Smax and the number of elements embedded in the plastic zones is plotted in Fig. 4. Solanki et al. [22] studied the crack closure behavior in compact specimens under 2-D plane stress and plane strain conditions, and their mesh refinement results are also shown in Fig. 4. Obviously the 2-D plane stress results show that the stabilized Sop/Smax values converge to a constant as the number of elements in the plastic zone increases. Despite the number of elements in the free-surface plastic zone of the present study is larger than those of the illustrated 2-D plane stress studies, the calculated stabilized Sop/Smax values of the present 3-D cracks show more fluctuations, and the converged Sop/Smax value cannot be acquired. However, it can be expected that, when the 2-D results are used for extrapolation to the region of larger number of elements, the results of the present study will be close to the extrapolated line (see Fig. 4). On the other hand, the stabilized Sop/Smax values of the 2-D plane strain studies become smaller and no converged trend is found as the number of elements
722
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
1 0.8
Sop /Smax
(b)
Skinner et al. [16] present study S max / σ 0 = 0.27 S max / σ 0 = 0.33 radial element size a/t = 0.34 a/t = 0.53 0.005 mm a/t = 0.71 0.01 mm
1
Kim et al. [19] Smax /σ 0 = 0.54
0.8
0.6
Sop /Smax
(a)
free surface 0.4
0.005 mm 0.01 mm
0.6
free surface deepest penetration point
0.4
deepest penetration point 0.2
present study Smax /σ 0 = 0.5 element size
0.2
0
0 0
0.3
0.6
0.9
1.2
0
0.3
Crack Increment (mm)
0.6
0.9
1.2
Crack Increment (mm)
Fig. 3. Comparison of the present restrained-front results with other finite element results: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
tio n
1.0
present 3-D analysis
ec g d ir a d in
deepest point, Smax /σ 0 = 0.33 free surface, Smax /σ0 = 0.5
0.8
0.4
lo
Sop /Smax Stablized
extrapolated line 0.6
2-D analysis [22] radial direction
plane stress plane strain
Fig. 5. Illustration of the element enlargement along the loading direction in the fine mesh region.
0.2
0.0
1
10
100
1000
Number of Elements in Plastic Zone Fig. 4. Effects of mesh refinement on the stabilized opening levels obtained from the restrained-front technique.
increases. The present 3-D results at the deepest penetration point show close stabilized Sop/Smax values to the 2-D results, but they are still with larger fluctuations. These fluctuations cause the ambiguous trend of the calculated stabilized Sop/Smax and it is uneasy to determine whether the 3-D stabilized Sop/Smax values exhibit the converged values. These results retard the mesh refinement study on the crack closure behavior of 3-D cracks. Since the mesh arrangement of this study follows the recommendation of Skinner and Daniewicz [16] (using of at least five elements in the plastic zone), there should be other reasons to cause these undesired fluctuations. Fig. 5 shows a small portion of the mesh arrangement in the fine mesh region, i.e., region A in Fig. 2. For the initial crack front, the mesh arrangement causes approximately five elements embedded in the deepestpenetration-point plastic zone for the 0.02-mm model when Smax/r0 = 0.33. However, along the loading direction, simply one layer of fine meshes is deployed and element enlargement is arranged for the upper layer. This mesh arrangement reduces the elements embedded in the plastic zone along the loading direction and is probably the cause of stabilized Sop/Smax fluctuations. The fine mesh layers along the loading direction are needed to be increased when the crack closure behavior is the main issue to be addressed. However, when the crack shape evolution is also
a concerned topic, this mesh refinement becomes impractical. For the crack closure analysis, computations are usually terminated when the Sop/Smax values become stable, and this phenomenon usually takes no more than 20–30 loading cycles when the applied load is not large. Within this number of loading cycles, the crack advances just a little and the observation of shape change for a further crack growth increment becomes impossible. To overcome this difficulty, computations for more loading cycles are needed. Hence, more elements must be deployed along the radial direction. For example, in the present study, there are 64, 128, and 256 elements deployed along the radial direction of the three models with the element sizes equal to 0.02, 0.01, and 0.005 mm, respectively. These mesh arrangements allowed the simulated cracks to advance approximately 1/10th of the initial crack dimension (5.25 mm), and the required loading cycles were nearly 190 loading cycles for the 0.005-mm model. For the computation facilities that the author can access, the required CPU time for the case of the 0.005-mm model is approximately two months. Hence, increasing the layer of fine meshes causes intolerable CPU time and becomes impractical at this stage. It is obvious that the mesh plan used in the present study is inadequate for the crack closure study. However, when the crack shape is the main issue to be addressed, as is discussed in the next section, a useful conclusion can be still reached based on the fluctuated results. 3.4. Simulated results of the free-front technique For the 0.005-mm model, Fig. 6 compares the calculated Sop/ Smax results of the free-front technique with those of the conventional restrained front. It is obvious that both techniques simulate
723
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
(b)
1
S max /σ 0 = 0.33
Sop /Smax
0.8
Smax /σ 0 = 0.5
free front restrained front
free front restrained front
0.6
free surface 0.4
0.2
0
1
0.8
Sop /Smax
(a)
0.3
0.6
0.9
free surface
0.4
0.2
deepest penetration point 0
0.6
deepest penetration point 0
1.2
0
0.3
0.6
0.9
1.2
Crack Increment (mm)
Crack Increment (mm)
Fig. 6. Comparison of the calculated crack opening levels between the restrained- and free-front techniques: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
(a)
(b)
1.0
1.0
Smax /σ 0 = 0.5 radial element size
S max /σ 0 = 0.33 radial element size
0.6
free surface 0.4
0.2
0
0.3
0.6
0.6
free surface 0.4
0.2
deepest penetration point
0
0.005 mm 0.01 mm 0.02 mm
0.8
0.005 mm 0.01 mm 0.02 mm
Sop /Smax
Sop /Smax
0.8
0.9
1.2
Crack Increment (mm)
0
deepest penetration point 0
0.3
0.6
0.9
1.2
Crack Increment (mm)
Fig. 7. Mesh refinement effects on the calculated opening levels using the free-front technique: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
close results except at the deepest penetration point where the Sop/ Smax values of the free-front technique deviate from the stable value when Smax/r0 = 0.5 (see Fig. 6b). Fig. 7 shows the mesh size effects on the calculated Sop/Smax values for the free-front technique and, for both loading conditions, the different radial element sizes and the mentioned insufficient layer of fine meshes lead to the Sop/Smax differences. The simulated crack shapes are plotted in Fig. 8. It is observed that the continuously higher opening level near the free surface (see Fig. 7) leads the front to deviate from the initial shape and to turn inward. When the simulated crack shapes of the three models are compared, no significant difference is found when the cracks advance a distance of approximately 1/10th of the initial crack dimension. In other words, even there are the mesh size and the insufficient fine mesh layer problems as discussed above to cause the closure level differences, these differences have little effect on the simulated crack shape. Based on these results, attention should be paid for the mesh refinement effects in the crack closure analysis. However, in the crack shape analysis, the plan of fine meshes seems to not be an important issue. This postulation is further verified in the next section by the models using different species of mesh arcs. 3.5. Simulated crack shape by the models with different species of mesh arcs The finite element analyses discussed in the previous section have been performed based on meshes of circular arcs, and the
shape evolution has been solely studied to a short distance from the initial state. It is also interesting to check whether the calculated crack shape development is influenced by different finite element models based on different species of mesh arcs and to check whether the crack shape change is significant when the crack propagates to a larger extent. Fig. 9 shows two different mesh plans on the crack plane based on different species of mesh arcs. Fig. 9a is the mesh plan using the concentric circular arcs. There are three fine mesh regions and are denoted as A–C. In all these three regions, 64 and 16 eight-node hexahedral elements are deployed along the radial and circumferential directions, respectively. Because these three regions are bordered by circular arcs with radii of 5, 8.2, 14.6, and 24.2 mm, the element sizes along the radial direction in regions A, B, and C are 0.05, 0.1, and 0.15 mm, respectively. The initial crack front is located at region A and is semi-circular with a radius of 5.25 mm. Fig. 9b is an alternative mesh plan for the same crack problem using elliptical arc meshes. The alternative model has three fine mesh regions, too. There are 96 elements deployed along the radial direction in each of the three fine mesh regions, and the element sizes along the radial and deepest-penetration directions are shown in Fig. 9b. Because the initial crack front is semi-circular, the major and minor axes of the two elliptical borders of region A are carefully chosen so that in the region one of the mesh arcs is circular with a radius of 5.25 mm. Fig. 10 shows the calculated Sop/Smax values obtained from the circle-based and ellipse-based models. Putra and Schijve [23]
724
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
(a)
(b) initial front
initial front
Smax/σ0 = 0.5
Smax /σ0 = 0.33
radial element size 0.005 mm 0.01 mm 0.02 mm
radial element size 0.005 mm 0.01 mm 0.02 mm
free surface
free surface
Fig. 8. Mesh refinement effects on the simulated crack shape using the free-front technique: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
(a)
(b)
60 mm
60 mm
24.25 0.10
96
B
0.06 0.04
B
96
8.89 5.05
64
A
A
96
64
0.12
0.08
5 8.2 14.6
4.85
24.2
25 mm
C
14.65
64
25 mm
C
12.53
0.20
24.05
43.25
Fig. 9. The circular- and elliptical-arc based mesh plans to study the effects of different species of finite element models on the simulated crack shape.
1
measured data [23]
Sop/Smax
0.8
Smax /σ0 = 0.5 ellipse-based meshes circle-based meshes
0.6
free surface
0.4
0.2
deepest penetration point 0
0
5
10
15
20
used in the tests and in the present simulations are different. Unfortunately, the Sop/Smax data along the free surface were not reported in Putra and Schijve’s paper, the data comparison is not made along the direction. Although the calculated crack opening levels of the two models have major differences along the direction of the deepest penetration point, it is interesting that the simulated crack shapes are similar without significant deviation even when the cracks advance to a deep depth of the specimen (see Fig. 11). In addition to the measured opening levels, in the same article Putra and Schijve also reported the crack shapes at various crack depths and these shapes are shown in Fig. 11. The simulated results successfully capture the front turning-inward phenomenon, however, the measured shape tends to have a more pronounced decreasing trend for the aspect ratio (a/c). The aspect ratio discrepancy is probably caused by the inaccurate estimations of the stress intensity factor in the simulations and this issue will be discussed more detail in the next section.
Crack Increment (mm) Fig. 10. The calculated Sop/Smax of the circle-based and the ellipse-based models.
reported the measured data along the direction of the deepest penetration point and these data are also shown in the figure. It is found that the simulated results are close to the measured data despite the specimen thickness and the loading conditions
ellipse-based model
4. Verification of the Newman and Raju’s assumption for the surface crack shape Newman and Raju [1] gave equations for the stress intensity factor of semi-elliptical cracks. In addition, they assumed that an initial semi-elliptical surface crack maintains a mathematicallydefined semi-elliptical shape during the crack propagation, and
circle-based model
specimen edge
t = 25 mm Putra and Schijve [23] t = 10 mm
Fig. 11. The simulated crack shape development of different models based on various species of mesh arcs and the measured crack shape development.
725
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
C C =C A ¼ 0:9
n
ð2Þ
where CC and CA are the empirical constants of the Paris equation at the free surface and the deepest penetration point of the semi-elliptical crack, respectively. By using such an assumption, the calculated crack aspect ratio evolution is forced to be consistent with the test data, but it seems to be illogical that the material constant of the Paris law depends on the location in an isotropic and homogeneous material such as aluminum and steel. Hence, an explanation is required for the above assumption. The following paragraphs illustrate how to explain the Newman and Raju’s assumption using the crack closure concept. When the closure effect is considered in a surface crack, the modified Paris equation for the growth rates at the free surface and the deepest penetration point can be written in this way:
dc=dN ¼ C 0 ðDK eff ;C Þn ¼ C 0 ðU C DK C Þn ¼ C 0 U nC ðDK C Þn ¼ C C ðDK C Þn
ð3Þ
da=dN ¼ C 0 ðDK eff ;A Þn ¼ C 0 ðU A DK A Þn ¼ C 0 U nA ðDK A Þn ¼ C A ðDK A Þn
ð4Þ
where c and a are the crack dimensions along the free surface and the deepest-penetration directions, respectively; C0 is the material constant for the modified Paris equation; DKeff,C, DKeff,A, UC and UA are the effective stress intensity factor and effective stress intensity ratio at the free surface and the deepest penetration point, respectively. From these two equations we can obtain CC = C0 U nC and CA = C0 U nA , that is, the material constant of the Paris equation may exhibit different values at the free surface and the deepest penetration point due to the crack closure effect. If the value of UC/UA is 0.9, the relation between CC and CA is identical to the Newman and Raju’s assumption. Pang [24] further modified Newman and Raju’s assumption for welded joints and used a value of 0.85 instead of 0.9 in Eq. (2). To verify whether the crack closure effects can be used to explain Newman and Raju’s or Pang’s assumption, Hou [17] tried to use the Newman and Raju’s two-point approach, combining the calculated crack opening levels obtained from the conventional finite element technique, in order to calculate the aspect ratio evolution of a semi-circular crack. The results showed that the predicted aspect ratio evolution was erroneous because the simulated UC/UA of the conventional crack closure analysis was approximately 0.63, which is far below 0.9. Therefore, we can conclude that crack closure concept cannot be used to explain the above assumption. However, the aspect ratio evolution was calculated based on the crack growth rate using the modified Paris equation considering both the crack closure and the stress intensity factor, the erroneous aspect ratio predictions could be attribute to the inaccurate stress intensity factor calculations which were based on the assumption that the front was mathematically-defined ellipses. This is true because the front turning-inward phenomenon near the free surface was not considered in Hou’s calculations (see Ref. [17]). This difficulty can be overcome by employing the free-front technique since the stress intensity factors near the free surface including the effects of front turning inward are always monitored during the calculations. Since the aspect ratio evolution can be obtained from the crack growth rate at points C and A, the ratio (dc/dN)/(da/dN) can be calculated instead of the aspect ratio:
ðdc=dNÞ=ðda=dNÞ ¼ ðDK eff ;C Þn =ðDK eff ;A Þn ¼ ððU C DK C Þ=ðU A DK A ÞÞn
ð5Þ
The model shown in Fig. 9a is used for the calculations. In addition, the crack growth rate ratio is also calculated using the Newman and Raju’s assumption without considering the crack closure effects:
ðdc=dNÞ=ðda=dNÞ ¼ 0:9n ðDK C Þn =ðDK A Þn
ð6Þ
Note that DKC and DKA in Eq. (6) are obtained using the equations proposed by Newman and Raju [1], that is, the crack front is assumed to be an ellipse. The crack growth rate ratios obtained from Eqs. (5) and (6) are plotted against the crack depth in Fig. 12. In stead of 0.9, the values 0.8 and 0.85 are also used in Eq. (6). Newman and Raju’s results show that, when the constant is 0.9, the growth rate ratio is close to 1 during the early stage of the crack growth, that is, the crack front growth rate at the free surface is close to that at the deepest penetration point. As the front moves forward, the ratio is larger than 1, indicating that the growth rate at the free surface becomes larger than that at the deepest penetration point, even if the assumption reduces the coefficient of the Paris equation at the free surface. The Smax/r0 = 0.5 and 0.67 results obtained from the freefront technique show that Smax has insignificant effects on the growth rate ratio except for the crack depth larger than 15 mm. Obviously, when the crack closure is considered and the stress intensity factor is not calculated based on the mathematical ellipses, the value 0.9 in the Newman and Raju’s assumption cannot still be reached by using the free-front technique. However, improvement has been made since the figure shows that the results of the free-front technique fit the value 0.8 in the best way. Such a value is much better than 0.63 which is obtained from the results when the crack front shape is restricted to the mathematical functions. The difference between 0.8 and 0.9 is probably caused by the incorrect stress singularity used at the intersection point between the front and the free surface. In the free-front calculations, the stress intensity factors are obtained using the wellknown r1/2 stress singularity. It has been shown by many researchers that the stress state at the intersection of the front and the free surface no longer presents the r1/2 stress singularity [20,25,26]. When the correct stress singularity (different from the r1/2) is used, the stress intensity factor near the free surface can be calculated more accurately, and it is believed that the aspect
2
Newman and Raju's assumption: n
C C /C Α = 0.9 n 0.85n 0.8
1.5
(d c/dN)/(da /dN)
estimated the growth increments at the free surface and the deepest penetration point of the crack using the Paris equation. In this approach, the crack front variation is only monitored at the above two points, and the crack dimensions along these two directions are used to estimate the aspect ratio evolution of the semi-elliptical cracks. When the calculated aspect ratios are compared to the test data, it seems that the Paris equation overestimates the crack growth rate at the free surface. Therefore, an assumption was made to reduce the growth rate at the free surface [1]:
1
Results by free-front technique Smax /σ0 = 0.5 Smax /σ0 = 0.67
0.5
0
0
5
10
15
20
25
Crack depth, a (mm) Fig. 12. Comparison of the crack growth rate ratio between the Newman and Raju’s assumption and the results of the free-front technique.
726
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
ratio evolution of a surface crack can be captured by taking into account the crack closure concept. 5. Conclusions (1) The mesh refinement is an important issue for the finite element crack closure analysis. However, the simulated crack shape is neither significantly influenced by different finite element models of various radial mesh sizes nor by the models based on different species of mesh arcs. The crack front turning inward near the free surface can be captured by the free-front technique. However, the simulated shapes show discrepant aspect ratio as compared to the measured data. (2) Based on the observed crack shape development, Newman and Raju correlated the crack growth rates at the free surface and the deepest penetration point by assuming different Paris constants at the above two points, i.e., CC/CA = 0.9n. The conventional finite element crack closure analysis shows that the CC/CA value is approximately 0.63n. When the crack shape is free to change in the crack closure analysis and the front turning inward caused stress intensity factor variation near the free surface is considered, the estimated CC/CA value is approximately 0.8, which is nearer to the Newman and Raju’s assumption. The 0.8 and 0.9 difference is probably caused by the incorrect r1/2 stress singularity used in the stress intensity factor estimations. It is believed that when the correct stress singularity is used in the analysis, the estimated CC/CA value and the aspect ratio evolution can be further improved. Acknowledgements The author would like to express thanks to National Center for High-Performance Computing of Taiwan for providing the computation facilities. References [1] Newman Jr JC, Raju IS. An empirical stress-intensity factor equation for the surface crack. Eng Fract Mech 1981;15:185–92. [2] Connolly MP, Collins R. The measurement and analysis of semi-elliptical surface fatigue crack growth. Eng Fract Mech 1987;26:897–911. [3] McFadyen NB, Bell R, Vosikovsky O. Fatigue crack growth of semi-elliptical surface cracks. Int J Fatigue 1990;12:43–50.
[4] Shin CS. Some aspect of corner fatigue crack growth from holes. Int J Fatigue 1991;13:233–40. [5] Carpinteri A. Shape change of surface cracks in round bars under cyclic axial loading. Int J Fatigue 1993;15:21–6. [6] Carpinteri A, Vantadori S. Sickle-shaped surface crack in a notched round bar under cyclic tension and bending. Fatigue Fract Eng Mater Struct 2009;32:223–32. [7] Lin XB, Smith RA. Shape growth simulation of surface cracks in tension fatigued round bar. Int J Fatigue 1997;19:461–9. [8] Couroneau N, Royer J. Simplified model for the fatigue growth analysis of surface cracks in round bars under mode I. Int J Fatigue 1998;20:711–8. [9] Song PS, Shieh YL. Crack growth and closure behaviour of surface cracks. Int J Fatigue 2004;26:429–36. [10] Heyder M, Kuhn G. 3D fatigue crack propagation: experimental studies. Int J Fatigue 2006;28:627–34. [11] Chipalo MI, Gilchrist MD, Smith RA. A finite element technique for the investigation of the shape development of planar cracks with initially irregular profiles. Int J Mech Sci 1990;32:243–51. [12] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates, part I: the numerical technique. Eng Fract Mech 1999;63:503–22. [13] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates, part II: crack shape change. Eng Fract Mech 1999;63:523–40. [14] Chermahini RG, Palmberg B, Blom AF. Fatigue crack growth and closure behavior of semicircular and semi-elliptical surface flaws. Int J Fatigue 1993;15:259–63. [15] Zhang JZ, Bowen P. On the finite element simulation of three-dimensional semi-circular fatigue crack growth and closure. Eng Fract Mech 1998;60:341–60. [16] Skinner JD, Daniewicz SR. Simulation of plasticity-induced fatigue crack closure in part-through cracked geometries using finite element analysis. Eng Fract Mech 2002;69:1–11. [17] Hou C-Y. Three-dimensional finite element analysis of fatigue crack closure behavior in surface flaws. Int J Fatigue 2004;26:1225–39. [18] Simandjuntak S, Alizadeh H, Pavier MJ, Smith DJ. Fatigue crack closure of a corner crack: a comparison of experimental results with finite element predictions. Int J Fatigue 2005;27:914–9. [19] Kim JS, Kang JY, Song JH. Elucidation of fatigue crack closure behavior in surface crack by 3D finite element analysis. Int J Fatigue 2007;29:168–80. [20] de Matos PFP, Nowell D. The influence of the Poisson’s ratio and corner point singularities in the three-dimensional plasticity-induced fatigue crack closure: a numerical study. Int J Fatigue 2008;30:1930–43. [21] Hou C-Y. Simultaneous simulation of closure behavior and shape development of fatigue surface cracks. Int J Fatigue 2008;30:1036–46. [22] Solanki K, Daniewicz SR, Newman Jr JC. Finite element modeling plasticityinduced crack closure with emphasis on geometry and mesh refinement effects. Eng Fract Mech 2003;70:1475–89. [23] Putra IS, Schijve J. Crack opening stress measurements of surface cracks in 7075-T6 aluminum alloy plate specimen through electron fractography. Fatigue Fract Eng Mater Struct 1992;15:323–38. [24] Pang HLJ. Fatigue crack growth and coalescence of surface cracks. In: Proceedings of the 12th conference offshore mechanics arctic engineering; 1993. p. 485–91. [25] Benthem JP. State of stress at the vertex of a quarter-infinite in a half-space. Int J Solid Struct 1977;13:479–92. [26] Bazant ZP, Estenssoro L. Surface singularity and crack propagation. Int J Solid Struct 1979;15:405–26.
Contents lists available at ScienceDirect
International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue
Simulation of surface crack shape evolution using the finite element technique and considering the crack closure effects Chien-Yuan Hou Department of Civil Engineering, National Chi Nan University, Puli, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 27 August 2010 Received in revised form 16 November 2010 Accepted 30 November 2010 Available online 4 December 2010 Keywords: Surface cracks Finite elements Crack shape Crack closure
a b s t r a c t Shape development and closure analysis of surface cracks using the finite element method have been employed by researchers without considering interactions of the two phenomena. This study uses the free-front finite element technique to account for such interactions. Mesh refinement effects on the crack shape development are studied by concluding that the element size ahead of the crack front is not affecting the simulated crack shape. Newman and Raju assumed different Paris equation coefficients for examining the surface crack growth. This assumption is here analysed by applying the crack closure concept using the results of the free-front technique. The results are close to those determined by employing the above assumption for the equation coefficients. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction One of the most important issues related to surface cracks subjected to fatigue loads is the shape development. Newman and Raju [1] were the pioneers with regards to this problem. They assumed that an initial semi-elliptical surface crack maintains a semi-elliptical shape during the whole propagation, and estimated the growth increments using the Paris law simply at the free surface and the deepest penetration point of the crack, i.e., points C and A in Fig. 1, respectively. Connolly and Collins [2] used the alternative current field measurement technique to analyse the shape change of elliptical surface cracks. McFadyen et al. [3] measured the shape evolution of surface cracks with various aspect ratio values and used the Newman and Raju’s approach to predict the crack aspect ratio change. Shin [4] observed the crack growth rate and shape development of corner cracks and calibrated the stress intensity factor of such cracks. The surface crack growth in round bars has been studied by some researchers [5–8]. Song and Shieh measured the crack closure behavior of surface cracks and found that the closure slowed down the crack growth rate at the surface [9]. Heyder and Kuhn observed the 3-D crack propagation in transparent PMMA specimens [10]. Since the shape of semi-elliptical surface cracks during propagation unnecessarily remains elliptical, Smith et al. [11] and Lin and Smith [12,13] developed a multiple point scheme which monitored not only the crack front growth at points A and C, but the growth at several front points between A and C was also moniE-mail address: [email protected] 0142-1123/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijfatigue.2010.11.022
tored. In the scheme, the shape assumption, i.e., the mathematically-defined ellipse, was no longer required, and different finite element models were constructed when a new crack front was given. Hence the finite element model was continuously changing due to the crack front evolution, and a computer code to generate the new model was required. Since the phenomenon of plasticity-induced crack closure was observed, it has been used to explain many metal fatigue phenomena such as mean stress effects, crack growth retardation, anomalous growth behavior of short cracks, etc. The 2-D simulation of crack closure by the finite element technique has achieved interesting results during the last twenty years of the 20th century. The 3-D crack closure analysis has been rare in that period probably because of computation facilities not sophisticated enough to perform the required complex computations. In recent years, the difficulty has been overcome and some results have been reported. Chermahini et al. [14] extended the 2-D closure analysis to semicircular and semi-elliptical cracks, but the mesh density of the finite element model was coarse due to the limitation of computation facilities. Zhang and Bowen [15] constructed a simple semi-circular crack finite element model to calculate the closure behavior when the crack was subjected to a single overload. Skinner and Daniewicz [16] performed finite element calculations for the closure level of semi-elliptical cracks. Hou [17] carried out calculations related on both semi-circular and -elliptical cracks at various loading levels and R ratios. Simandjuntak et al. studied the closure behavior of corner cracks emanating from a hole, and compared the calculated results to the observed data [18]. Kim et al. [19] investigated the stabilization behavior of the crack
720
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
opening level, the effect of mesh size, and the effect of initial crack length in the 3-D closure analysis. de Matos and Nowell [20] studied the crack front geometry influenced by the stress singularity near the free surface, and the crack opening stresses were calculated using the modified crack shape. Obviously, the crack shape development and the crack closure are two factors interacting each other when a 3-D crack propagates. For a surface crack, the crack closure effect is more pronounced at the free surface due to the plane stress state. As a result, the crack growth rate at the free surface is reduced and a front-turning-inward phenomenon is usually observed (see Fig. 1). The crack front turning inward near the free surface increases the stress intensity factor at this region and subsequently causes different crack closure levels. Unfortunately, none of the above mentioned 3-D finite element crack shape [1,11–13] or crack closure [14–20] studies, or other similar studies reported in literature take the interactions of the two phenomena into consideration. In a conventional 3-D finite element crack closure analysis, the crack shape is primarily considered in the stage of the finite element model construction. As the model is determined, the initial and the subsequent crack shapes cannot be altered due to the crack front nodes are released simultaneously, that is, the evolving crack shape must follow the prescribed mesh arcs. If the technique is to be used for simultaneous simulation of the crack closure and crack shape development, it is obvious that crack front must be free to move without restraining to the prescribed mesh arcs in the crack closure analysis. To solve the difficulty, Hou [21] developed a finite element technique without releasing the crack front nodes simultaneously. This strategy combines the multiple point scheme [11–13] and the conventional crack closure analysis and, hence, allows the crack front free to move without restraining so that the evolving crack shape can be simulated and the crack opening stresses can also be evaluated concurrently. The simulated crack shapes and crack opening levels determined by this new technique seemed to be promising. As a subsequent study, this paper uses the above free-front technique further addressing the effects of the mesh refinement and employing different meshes by considering various species of arcs on the simulated crack shape. In addition, the Newman and Raju’s assumption related to the surface crack shape development is assessed using the results of the above technique. 2. The free-front finite element technique for the crack closure simulation Fig. 1 schematically shows the crack-plane meshes of a 3-D finite element crack closure analysis. The initial crack front in the
conventional crack closure analysis is coincident to a mesh arc, and the front nodes are released simultaneously during a loading cycle. As a consequence, the subsequent crack front is also coincident to one of the mesh arcs: see fronts C0 to Cj in Fig. 1. Hence, once the mesh plan is determined, the crack shape is restrained for the following crack closure analysis, and the crack shape change cannot be captured. Hou developed a finite element technique so that both the crack closure and shape can be simultaneously simulated [21]. In such a technique, two finite element models are constructed: one is used to calculate the stress intensity factor at each front node and the other is used for the crack closure analysis. Since the first model is used for the stress intensity factors, the special elements of quartered point are deployed along the front so that the stress singularity can be captured and hence the model is analysed elastically. Elastic–plastic analysis is performed by using the second model, in a way similar to the conventional closure analysis, and the opening stresses (Sop) at the front nodes are determined. As the stress intensity factors and opening stresses are obtained, the effective stress intensity factor can be calculated at all front nodes. Along the crack front, the node with the maximum effective stress intensity factor is determined and the crack increment at this node is set to one radial element size for the closure analysis of the next loading cycle. The ratios of other front nodal increments to the increment of the maximum-effective-stress-intensity node are then calculated using the modified Paris law:
ðDK eff Þni Dai ¼ Damax ðDK eff Þnmax
ð1Þ
where Dai and Damax are the crack increments at the front node i and at the node with the maximum increment, respectively; DKeff is the effective stress intensity range; and n is the power exponent of the Paris law. The increment of other front nodes must be lower than one radial element size, and is accumulated with the increment of the previous loading cycle. Once the accumulated increment at a front node is larger than one element size, the node is released for the closure analysis of the next loading cycle. Hence, in the proposed technique, the front nodes are not released simultaneously in one loading cycle and the crack front is not restrained to the prescribed mesh arcs during the analysis. Thus the crack shape evolution can be simulated. For example, using such a technique, a crack front can be formed by the solid circular symbols, as is shown in Fig. 1. By employing the proposed technique, the calculated crack face profile is rugged as compared to that of the conventional closure analysis, but the simulated opening stresses and crack shape are promising. 3. Effects of mesh plan on the simulated crack shape
C0: initial crack front Cj: crack front of the j th load cycle
Cj C0
free surface
a
A: point of deepest penetration
c
C front turning inward
Fig. 1. Schematic illustration of the restrained-front and free-front techniques for the finite element crack closure analysis.
3.1. The crack closure model and loading conditions The issue of the mesh refinement has received numerous attentions in the crack closure analysis since different crack opening levels can be obtained from finite element models with various mesh sizes. It is also interesting to assess whether the simulated crack shape is influenced by various models with different mesh sizes. A semi-circular surface crack with a radius of 5.25 mm is here examined. The crack-plane geometry of the analysed specimen and the finite element meshes of the specimen are shown in Fig. 2. Because of the symmetrical conditions of the specimen, the model is constructed simply for a quarter of the specimen. Eight-node hexahedral elements are used in the model for the crack closure analysis. Region A (bordered by thick lines in Fig. 2) is a fine mesh area, and element sizes of 0.005, 0.01, and
721
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
25 mm
A
initial crack front
60 mm
(a)
(b)
Fig. 2. (a) The dimensions of the crack plane and the mesh plan, and (b) the finite element meshes of the analysed specimen.
0.02 mm are respectively placed along the radial direction of three different models. Along the circumferential direction of such three models, 32 elements are deployed. The material is assume to be elastic-perfectly plastic with a plastic flow stress r0 = 360 MPa. The Young’s modulus and Poisson’ ratio are assumed to be 200 GPa and 0.3, respectively. Constant amplitude loads perpendicular to the crack plane with R = 0, Smax/r0 = 0.33 and 0.5 are applied (Smax = maximum stress). Skinner and Daniewicz [16] suggested that five elements embedded in the crack front plastic zone are likely to be required for the convergence of the stabilized crack opening stresses. The above radial element sizes are chosen in this way because the model with 0.02mm elements produces a plastic zone close to the deepest penetration point containing approximately five elements when Smax/ r0 = 0.33. 3.2. Crack face node opening and closing criteria The loading-direction displacements of crack face nodes are monitored in each unloading step. When any crack face node displacement is found to be negative indicating that the node is closed, the node displacement along the loading direction is fixed to zero for the following unloading steps. During the loading phase, the loading-direction reaction forces of the fixed crack face nodes are monitored. If a face node is still in contact, the direction of the reaction force should be the same as that of the external load. Once the direction of the reaction force becomes opposite to that of the external load, the fixed node is released for the following loading steps, that is, the crack face at the node is opened. 3.3. The simulated results of the conventional crack closure analysis Fig. 3 shows the calculated relationship between Sop/Smax and crack increment using the conventional restrained front approach. Note that, for both loading conditions (Smax/r0 = 0.33 and 0.5), the calculated ratio Sop/Smax for the 0.01- and 0.005-mm models shows an increasing and stabilized trend at the free surface, whereas such a ratio at the deepest penetration point shows an increasing, then decreasing and finally stable trend. Skinner and Daniewicz [16] performed the crack closure analysis of semi-circular cracks with the restrained crack front technique. The applied loads were Smax/r0 = 0.27 and R = 0.1. Since these conditions are close to Smax/r0 = 0.33 and R = 0, their results are used in Fig. 3a for comparison. In their study, three different
initial crack dimensions with the ratio of the crack depth to the specimen thickness a/t = 0.34, 0.53, 0.71 were considered and ten loading cycles were carried out. This number of loading cycles allows the cracks to travel to a maximum crack increment approximately 0.05 mm, which is much smaller than the crack increment range considered in the present study. Their results show that along the direction of the deepest penetration point the values of Sop/Smax become stable in ten loading cycles. Hence, it is assumed that the calculated Sop/Smax values remain stable to a larger crack increment and these values are illustrated as the hollow symbols in Fig. 3a. Note that these symbols are plotted at an arbitrarily chosen crack increment. It is obvious that the present results are slight smaller but close to those reported values. Kim et al. [19] performed the crack closure analysis for Smax/ r0 = 0.54 and R = 0. The reported Sop/Smax values are illustrated in Fig. 3b for comparison. At the free surface, results of Kim et al. reach the stable value earlier than those of the present study but both results show close stable values. Along the direction of the deepest penetration point, both studies show the similar trend with different stable values. These differences are mainly caused by (1) different initial crack dimensions, (2) different stress–strain relationships, and (3) different node release scheme used in the calculations. The plastic zone sizes at the free surface and the deepest penetration point were estimated using Irwin’s approximation, thus the number of elements embedded in these two plastic zones can be calculated. The relationship between the calculated stabilized Sop/Smax and the number of elements embedded in the plastic zones is plotted in Fig. 4. Solanki et al. [22] studied the crack closure behavior in compact specimens under 2-D plane stress and plane strain conditions, and their mesh refinement results are also shown in Fig. 4. Obviously the 2-D plane stress results show that the stabilized Sop/Smax values converge to a constant as the number of elements in the plastic zone increases. Despite the number of elements in the free-surface plastic zone of the present study is larger than those of the illustrated 2-D plane stress studies, the calculated stabilized Sop/Smax values of the present 3-D cracks show more fluctuations, and the converged Sop/Smax value cannot be acquired. However, it can be expected that, when the 2-D results are used for extrapolation to the region of larger number of elements, the results of the present study will be close to the extrapolated line (see Fig. 4). On the other hand, the stabilized Sop/Smax values of the 2-D plane strain studies become smaller and no converged trend is found as the number of elements
722
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
1 0.8
Sop /Smax
(b)
Skinner et al. [16] present study S max / σ 0 = 0.27 S max / σ 0 = 0.33 radial element size a/t = 0.34 a/t = 0.53 0.005 mm a/t = 0.71 0.01 mm
1
Kim et al. [19] Smax /σ 0 = 0.54
0.8
0.6
Sop /Smax
(a)
free surface 0.4
0.005 mm 0.01 mm
0.6
free surface deepest penetration point
0.4
deepest penetration point 0.2
present study Smax /σ 0 = 0.5 element size
0.2
0
0 0
0.3
0.6
0.9
1.2
0
0.3
Crack Increment (mm)
0.6
0.9
1.2
Crack Increment (mm)
Fig. 3. Comparison of the present restrained-front results with other finite element results: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
tio n
1.0
present 3-D analysis
ec g d ir a d in
deepest point, Smax /σ 0 = 0.33 free surface, Smax /σ0 = 0.5
0.8
0.4
lo
Sop /Smax Stablized
extrapolated line 0.6
2-D analysis [22] radial direction
plane stress plane strain
Fig. 5. Illustration of the element enlargement along the loading direction in the fine mesh region.
0.2
0.0
1
10
100
1000
Number of Elements in Plastic Zone Fig. 4. Effects of mesh refinement on the stabilized opening levels obtained from the restrained-front technique.
increases. The present 3-D results at the deepest penetration point show close stabilized Sop/Smax values to the 2-D results, but they are still with larger fluctuations. These fluctuations cause the ambiguous trend of the calculated stabilized Sop/Smax and it is uneasy to determine whether the 3-D stabilized Sop/Smax values exhibit the converged values. These results retard the mesh refinement study on the crack closure behavior of 3-D cracks. Since the mesh arrangement of this study follows the recommendation of Skinner and Daniewicz [16] (using of at least five elements in the plastic zone), there should be other reasons to cause these undesired fluctuations. Fig. 5 shows a small portion of the mesh arrangement in the fine mesh region, i.e., region A in Fig. 2. For the initial crack front, the mesh arrangement causes approximately five elements embedded in the deepestpenetration-point plastic zone for the 0.02-mm model when Smax/r0 = 0.33. However, along the loading direction, simply one layer of fine meshes is deployed and element enlargement is arranged for the upper layer. This mesh arrangement reduces the elements embedded in the plastic zone along the loading direction and is probably the cause of stabilized Sop/Smax fluctuations. The fine mesh layers along the loading direction are needed to be increased when the crack closure behavior is the main issue to be addressed. However, when the crack shape evolution is also
a concerned topic, this mesh refinement becomes impractical. For the crack closure analysis, computations are usually terminated when the Sop/Smax values become stable, and this phenomenon usually takes no more than 20–30 loading cycles when the applied load is not large. Within this number of loading cycles, the crack advances just a little and the observation of shape change for a further crack growth increment becomes impossible. To overcome this difficulty, computations for more loading cycles are needed. Hence, more elements must be deployed along the radial direction. For example, in the present study, there are 64, 128, and 256 elements deployed along the radial direction of the three models with the element sizes equal to 0.02, 0.01, and 0.005 mm, respectively. These mesh arrangements allowed the simulated cracks to advance approximately 1/10th of the initial crack dimension (5.25 mm), and the required loading cycles were nearly 190 loading cycles for the 0.005-mm model. For the computation facilities that the author can access, the required CPU time for the case of the 0.005-mm model is approximately two months. Hence, increasing the layer of fine meshes causes intolerable CPU time and becomes impractical at this stage. It is obvious that the mesh plan used in the present study is inadequate for the crack closure study. However, when the crack shape is the main issue to be addressed, as is discussed in the next section, a useful conclusion can be still reached based on the fluctuated results. 3.4. Simulated results of the free-front technique For the 0.005-mm model, Fig. 6 compares the calculated Sop/ Smax results of the free-front technique with those of the conventional restrained front. It is obvious that both techniques simulate
723
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
(b)
1
S max /σ 0 = 0.33
Sop /Smax
0.8
Smax /σ 0 = 0.5
free front restrained front
free front restrained front
0.6
free surface 0.4
0.2
0
1
0.8
Sop /Smax
(a)
0.3
0.6
0.9
free surface
0.4
0.2
deepest penetration point 0
0.6
deepest penetration point 0
1.2
0
0.3
0.6
0.9
1.2
Crack Increment (mm)
Crack Increment (mm)
Fig. 6. Comparison of the calculated crack opening levels between the restrained- and free-front techniques: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
(a)
(b)
1.0
1.0
Smax /σ 0 = 0.5 radial element size
S max /σ 0 = 0.33 radial element size
0.6
free surface 0.4
0.2
0
0.3
0.6
0.6
free surface 0.4
0.2
deepest penetration point
0
0.005 mm 0.01 mm 0.02 mm
0.8
0.005 mm 0.01 mm 0.02 mm
Sop /Smax
Sop /Smax
0.8
0.9
1.2
Crack Increment (mm)
0
deepest penetration point 0
0.3
0.6
0.9
1.2
Crack Increment (mm)
Fig. 7. Mesh refinement effects on the calculated opening levels using the free-front technique: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
close results except at the deepest penetration point where the Sop/ Smax values of the free-front technique deviate from the stable value when Smax/r0 = 0.5 (see Fig. 6b). Fig. 7 shows the mesh size effects on the calculated Sop/Smax values for the free-front technique and, for both loading conditions, the different radial element sizes and the mentioned insufficient layer of fine meshes lead to the Sop/Smax differences. The simulated crack shapes are plotted in Fig. 8. It is observed that the continuously higher opening level near the free surface (see Fig. 7) leads the front to deviate from the initial shape and to turn inward. When the simulated crack shapes of the three models are compared, no significant difference is found when the cracks advance a distance of approximately 1/10th of the initial crack dimension. In other words, even there are the mesh size and the insufficient fine mesh layer problems as discussed above to cause the closure level differences, these differences have little effect on the simulated crack shape. Based on these results, attention should be paid for the mesh refinement effects in the crack closure analysis. However, in the crack shape analysis, the plan of fine meshes seems to not be an important issue. This postulation is further verified in the next section by the models using different species of mesh arcs. 3.5. Simulated crack shape by the models with different species of mesh arcs The finite element analyses discussed in the previous section have been performed based on meshes of circular arcs, and the
shape evolution has been solely studied to a short distance from the initial state. It is also interesting to check whether the calculated crack shape development is influenced by different finite element models based on different species of mesh arcs and to check whether the crack shape change is significant when the crack propagates to a larger extent. Fig. 9 shows two different mesh plans on the crack plane based on different species of mesh arcs. Fig. 9a is the mesh plan using the concentric circular arcs. There are three fine mesh regions and are denoted as A–C. In all these three regions, 64 and 16 eight-node hexahedral elements are deployed along the radial and circumferential directions, respectively. Because these three regions are bordered by circular arcs with radii of 5, 8.2, 14.6, and 24.2 mm, the element sizes along the radial direction in regions A, B, and C are 0.05, 0.1, and 0.15 mm, respectively. The initial crack front is located at region A and is semi-circular with a radius of 5.25 mm. Fig. 9b is an alternative mesh plan for the same crack problem using elliptical arc meshes. The alternative model has three fine mesh regions, too. There are 96 elements deployed along the radial direction in each of the three fine mesh regions, and the element sizes along the radial and deepest-penetration directions are shown in Fig. 9b. Because the initial crack front is semi-circular, the major and minor axes of the two elliptical borders of region A are carefully chosen so that in the region one of the mesh arcs is circular with a radius of 5.25 mm. Fig. 10 shows the calculated Sop/Smax values obtained from the circle-based and ellipse-based models. Putra and Schijve [23]
724
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
(a)
(b) initial front
initial front
Smax/σ0 = 0.5
Smax /σ0 = 0.33
radial element size 0.005 mm 0.01 mm 0.02 mm
radial element size 0.005 mm 0.01 mm 0.02 mm
free surface
free surface
Fig. 8. Mesh refinement effects on the simulated crack shape using the free-front technique: (a) Smax/r0 = 0.33, and (b) Smax/r0 = 0.5.
(a)
(b)
60 mm
60 mm
24.25 0.10
96
B
0.06 0.04
B
96
8.89 5.05
64
A
A
96
64
0.12
0.08
5 8.2 14.6
4.85
24.2
25 mm
C
14.65
64
25 mm
C
12.53
0.20
24.05
43.25
Fig. 9. The circular- and elliptical-arc based mesh plans to study the effects of different species of finite element models on the simulated crack shape.
1
measured data [23]
Sop/Smax
0.8
Smax /σ0 = 0.5 ellipse-based meshes circle-based meshes
0.6
free surface
0.4
0.2
deepest penetration point 0
0
5
10
15
20
used in the tests and in the present simulations are different. Unfortunately, the Sop/Smax data along the free surface were not reported in Putra and Schijve’s paper, the data comparison is not made along the direction. Although the calculated crack opening levels of the two models have major differences along the direction of the deepest penetration point, it is interesting that the simulated crack shapes are similar without significant deviation even when the cracks advance to a deep depth of the specimen (see Fig. 11). In addition to the measured opening levels, in the same article Putra and Schijve also reported the crack shapes at various crack depths and these shapes are shown in Fig. 11. The simulated results successfully capture the front turning-inward phenomenon, however, the measured shape tends to have a more pronounced decreasing trend for the aspect ratio (a/c). The aspect ratio discrepancy is probably caused by the inaccurate estimations of the stress intensity factor in the simulations and this issue will be discussed more detail in the next section.
Crack Increment (mm) Fig. 10. The calculated Sop/Smax of the circle-based and the ellipse-based models.
reported the measured data along the direction of the deepest penetration point and these data are also shown in the figure. It is found that the simulated results are close to the measured data despite the specimen thickness and the loading conditions
ellipse-based model
4. Verification of the Newman and Raju’s assumption for the surface crack shape Newman and Raju [1] gave equations for the stress intensity factor of semi-elliptical cracks. In addition, they assumed that an initial semi-elliptical surface crack maintains a mathematicallydefined semi-elliptical shape during the crack propagation, and
circle-based model
specimen edge
t = 25 mm Putra and Schijve [23] t = 10 mm
Fig. 11. The simulated crack shape development of different models based on various species of mesh arcs and the measured crack shape development.
725
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
C C =C A ¼ 0:9
n
ð2Þ
where CC and CA are the empirical constants of the Paris equation at the free surface and the deepest penetration point of the semi-elliptical crack, respectively. By using such an assumption, the calculated crack aspect ratio evolution is forced to be consistent with the test data, but it seems to be illogical that the material constant of the Paris law depends on the location in an isotropic and homogeneous material such as aluminum and steel. Hence, an explanation is required for the above assumption. The following paragraphs illustrate how to explain the Newman and Raju’s assumption using the crack closure concept. When the closure effect is considered in a surface crack, the modified Paris equation for the growth rates at the free surface and the deepest penetration point can be written in this way:
dc=dN ¼ C 0 ðDK eff ;C Þn ¼ C 0 ðU C DK C Þn ¼ C 0 U nC ðDK C Þn ¼ C C ðDK C Þn
ð3Þ
da=dN ¼ C 0 ðDK eff ;A Þn ¼ C 0 ðU A DK A Þn ¼ C 0 U nA ðDK A Þn ¼ C A ðDK A Þn
ð4Þ
where c and a are the crack dimensions along the free surface and the deepest-penetration directions, respectively; C0 is the material constant for the modified Paris equation; DKeff,C, DKeff,A, UC and UA are the effective stress intensity factor and effective stress intensity ratio at the free surface and the deepest penetration point, respectively. From these two equations we can obtain CC = C0 U nC and CA = C0 U nA , that is, the material constant of the Paris equation may exhibit different values at the free surface and the deepest penetration point due to the crack closure effect. If the value of UC/UA is 0.9, the relation between CC and CA is identical to the Newman and Raju’s assumption. Pang [24] further modified Newman and Raju’s assumption for welded joints and used a value of 0.85 instead of 0.9 in Eq. (2). To verify whether the crack closure effects can be used to explain Newman and Raju’s or Pang’s assumption, Hou [17] tried to use the Newman and Raju’s two-point approach, combining the calculated crack opening levels obtained from the conventional finite element technique, in order to calculate the aspect ratio evolution of a semi-circular crack. The results showed that the predicted aspect ratio evolution was erroneous because the simulated UC/UA of the conventional crack closure analysis was approximately 0.63, which is far below 0.9. Therefore, we can conclude that crack closure concept cannot be used to explain the above assumption. However, the aspect ratio evolution was calculated based on the crack growth rate using the modified Paris equation considering both the crack closure and the stress intensity factor, the erroneous aspect ratio predictions could be attribute to the inaccurate stress intensity factor calculations which were based on the assumption that the front was mathematically-defined ellipses. This is true because the front turning-inward phenomenon near the free surface was not considered in Hou’s calculations (see Ref. [17]). This difficulty can be overcome by employing the free-front technique since the stress intensity factors near the free surface including the effects of front turning inward are always monitored during the calculations. Since the aspect ratio evolution can be obtained from the crack growth rate at points C and A, the ratio (dc/dN)/(da/dN) can be calculated instead of the aspect ratio:
ðdc=dNÞ=ðda=dNÞ ¼ ðDK eff ;C Þn =ðDK eff ;A Þn ¼ ððU C DK C Þ=ðU A DK A ÞÞn
ð5Þ
The model shown in Fig. 9a is used for the calculations. In addition, the crack growth rate ratio is also calculated using the Newman and Raju’s assumption without considering the crack closure effects:
ðdc=dNÞ=ðda=dNÞ ¼ 0:9n ðDK C Þn =ðDK A Þn
ð6Þ
Note that DKC and DKA in Eq. (6) are obtained using the equations proposed by Newman and Raju [1], that is, the crack front is assumed to be an ellipse. The crack growth rate ratios obtained from Eqs. (5) and (6) are plotted against the crack depth in Fig. 12. In stead of 0.9, the values 0.8 and 0.85 are also used in Eq. (6). Newman and Raju’s results show that, when the constant is 0.9, the growth rate ratio is close to 1 during the early stage of the crack growth, that is, the crack front growth rate at the free surface is close to that at the deepest penetration point. As the front moves forward, the ratio is larger than 1, indicating that the growth rate at the free surface becomes larger than that at the deepest penetration point, even if the assumption reduces the coefficient of the Paris equation at the free surface. The Smax/r0 = 0.5 and 0.67 results obtained from the freefront technique show that Smax has insignificant effects on the growth rate ratio except for the crack depth larger than 15 mm. Obviously, when the crack closure is considered and the stress intensity factor is not calculated based on the mathematical ellipses, the value 0.9 in the Newman and Raju’s assumption cannot still be reached by using the free-front technique. However, improvement has been made since the figure shows that the results of the free-front technique fit the value 0.8 in the best way. Such a value is much better than 0.63 which is obtained from the results when the crack front shape is restricted to the mathematical functions. The difference between 0.8 and 0.9 is probably caused by the incorrect stress singularity used at the intersection point between the front and the free surface. In the free-front calculations, the stress intensity factors are obtained using the wellknown r1/2 stress singularity. It has been shown by many researchers that the stress state at the intersection of the front and the free surface no longer presents the r1/2 stress singularity [20,25,26]. When the correct stress singularity (different from the r1/2) is used, the stress intensity factor near the free surface can be calculated more accurately, and it is believed that the aspect
2
Newman and Raju's assumption: n
C C /C Α = 0.9 n 0.85n 0.8
1.5
(d c/dN)/(da /dN)
estimated the growth increments at the free surface and the deepest penetration point of the crack using the Paris equation. In this approach, the crack front variation is only monitored at the above two points, and the crack dimensions along these two directions are used to estimate the aspect ratio evolution of the semi-elliptical cracks. When the calculated aspect ratios are compared to the test data, it seems that the Paris equation overestimates the crack growth rate at the free surface. Therefore, an assumption was made to reduce the growth rate at the free surface [1]:
1
Results by free-front technique Smax /σ0 = 0.5 Smax /σ0 = 0.67
0.5
0
0
5
10
15
20
25
Crack depth, a (mm) Fig. 12. Comparison of the crack growth rate ratio between the Newman and Raju’s assumption and the results of the free-front technique.
726
C.-Y. Hou / International Journal of Fatigue 33 (2011) 719–726
ratio evolution of a surface crack can be captured by taking into account the crack closure concept. 5. Conclusions (1) The mesh refinement is an important issue for the finite element crack closure analysis. However, the simulated crack shape is neither significantly influenced by different finite element models of various radial mesh sizes nor by the models based on different species of mesh arcs. The crack front turning inward near the free surface can be captured by the free-front technique. However, the simulated shapes show discrepant aspect ratio as compared to the measured data. (2) Based on the observed crack shape development, Newman and Raju correlated the crack growth rates at the free surface and the deepest penetration point by assuming different Paris constants at the above two points, i.e., CC/CA = 0.9n. The conventional finite element crack closure analysis shows that the CC/CA value is approximately 0.63n. When the crack shape is free to change in the crack closure analysis and the front turning inward caused stress intensity factor variation near the free surface is considered, the estimated CC/CA value is approximately 0.8, which is nearer to the Newman and Raju’s assumption. The 0.8 and 0.9 difference is probably caused by the incorrect r1/2 stress singularity used in the stress intensity factor estimations. It is believed that when the correct stress singularity is used in the analysis, the estimated CC/CA value and the aspect ratio evolution can be further improved. Acknowledgements The author would like to express thanks to National Center for High-Performance Computing of Taiwan for providing the computation facilities. References [1] Newman Jr JC, Raju IS. An empirical stress-intensity factor equation for the surface crack. Eng Fract Mech 1981;15:185–92. [2] Connolly MP, Collins R. The measurement and analysis of semi-elliptical surface fatigue crack growth. Eng Fract Mech 1987;26:897–911. [3] McFadyen NB, Bell R, Vosikovsky O. Fatigue crack growth of semi-elliptical surface cracks. Int J Fatigue 1990;12:43–50.
[4] Shin CS. Some aspect of corner fatigue crack growth from holes. Int J Fatigue 1991;13:233–40. [5] Carpinteri A. Shape change of surface cracks in round bars under cyclic axial loading. Int J Fatigue 1993;15:21–6. [6] Carpinteri A, Vantadori S. Sickle-shaped surface crack in a notched round bar under cyclic tension and bending. Fatigue Fract Eng Mater Struct 2009;32:223–32. [7] Lin XB, Smith RA. Shape growth simulation of surface cracks in tension fatigued round bar. Int J Fatigue 1997;19:461–9. [8] Couroneau N, Royer J. Simplified model for the fatigue growth analysis of surface cracks in round bars under mode I. Int J Fatigue 1998;20:711–8. [9] Song PS, Shieh YL. Crack growth and closure behaviour of surface cracks. Int J Fatigue 2004;26:429–36. [10] Heyder M, Kuhn G. 3D fatigue crack propagation: experimental studies. Int J Fatigue 2006;28:627–34. [11] Chipalo MI, Gilchrist MD, Smith RA. A finite element technique for the investigation of the shape development of planar cracks with initially irregular profiles. Int J Mech Sci 1990;32:243–51. [12] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates, part I: the numerical technique. Eng Fract Mech 1999;63:503–22. [13] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates, part II: crack shape change. Eng Fract Mech 1999;63:523–40. [14] Chermahini RG, Palmberg B, Blom AF. Fatigue crack growth and closure behavior of semicircular and semi-elliptical surface flaws. Int J Fatigue 1993;15:259–63. [15] Zhang JZ, Bowen P. On the finite element simulation of three-dimensional semi-circular fatigue crack growth and closure. Eng Fract Mech 1998;60:341–60. [16] Skinner JD, Daniewicz SR. Simulation of plasticity-induced fatigue crack closure in part-through cracked geometries using finite element analysis. Eng Fract Mech 2002;69:1–11. [17] Hou C-Y. Three-dimensional finite element analysis of fatigue crack closure behavior in surface flaws. Int J Fatigue 2004;26:1225–39. [18] Simandjuntak S, Alizadeh H, Pavier MJ, Smith DJ. Fatigue crack closure of a corner crack: a comparison of experimental results with finite element predictions. Int J Fatigue 2005;27:914–9. [19] Kim JS, Kang JY, Song JH. Elucidation of fatigue crack closure behavior in surface crack by 3D finite element analysis. Int J Fatigue 2007;29:168–80. [20] de Matos PFP, Nowell D. The influence of the Poisson’s ratio and corner point singularities in the three-dimensional plasticity-induced fatigue crack closure: a numerical study. Int J Fatigue 2008;30:1930–43. [21] Hou C-Y. Simultaneous simulation of closure behavior and shape development of fatigue surface cracks. Int J Fatigue 2008;30:1036–46. [22] Solanki K, Daniewicz SR, Newman Jr JC. Finite element modeling plasticityinduced crack closure with emphasis on geometry and mesh refinement effects. Eng Fract Mech 2003;70:1475–89. [23] Putra IS, Schijve J. Crack opening stress measurements of surface cracks in 7075-T6 aluminum alloy plate specimen through electron fractography. Fatigue Fract Eng Mater Struct 1992;15:323–38. [24] Pang HLJ. Fatigue crack growth and coalescence of surface cracks. In: Proceedings of the 12th conference offshore mechanics arctic engineering; 1993. p. 485–91. [25] Benthem JP. State of stress at the vertex of a quarter-infinite in a half-space. Int J Solid Struct 1977;13:479–92. [26] Bazant ZP, Estenssoro L. Surface singularity and crack propagation. Int J Solid Struct 1979;15:405–26.