- Email: [email protected]
Finite element based simulation of fatigue crack growth with a focus on elastic–plastic material behavior
Computational Materials Science 57 (2012) 73–79
Contents lists available at SciVerse ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Finite element based simulation of fatigue crack growth with a focus on elastic–plastic material behavior P. Zerres ⇑, M. Vormwald Materials Mechanics Group, Technische Universität Darmstadt, Petersenstr. 12, D-64287 Darmstadt, Germany
a r t i c l e
i n f o
Article history: Received 26 November 2010 Received in revised form 23 December 2011 Accepted 13 January 2012
Keywords: Fatigue crack growth Crack closure Finite element method
a b s t r a c t The numerical crack growth simulation has become more and more attractive today for reasons of time and economy. The presentation of a procedure for such a simulation of fatigue crack growth in elastic– plastic materials is the aim of this paper. The remeshing of the structure and the subsequent mapping of the status variables from the old mesh to the new one after each step of crack advance are thereby the core characteristics of the proposed procedure. The crack growth life is determined by integration of a crack growth law which is based on the effective range of the stress intensity factor. The results obtained by the procedure for autofrettaged Diesel-engine injection tubes are compared with experimental results and finite element results based on a nodal release scheme. A good agreement between the two numerical analyses and between the results with the proposed procedure and the experimental results concerning the endurance limit can be stated. By looking at the fatigue lives the experimental results are underestimated by the numerical ones. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction For calculating the fatigue crack growth, which often is responsible for failure in structures under predominant cyclic loading, reliable and in research established methods differ in numerical crack growth simulations based either on linear-elastic or on elastic–plastic fracture mechanics. Choosing the first approach the assumption of small scale yielding is necessary. Research in this field focus on the description of the crack path. The range of the stress intensity factor is the fundamental parameter together with its connection to the crack growth rate, which has to be determined empirically. Two methods are successful in use today to calculate the stress intensity factors, the finite element method (FEM) and the boundary element method (BEM). The crack growth simulation based on linear-elastic fracture mechanics is implemented in several two- and three-dimensional, commercial and non-commercial software tools, such as FRANC2D [1,2], FRANC3D [3,4], ZEN CRACK [5,6], ADAPCRACK3D [7,8] and BEASY [9]. Using the FEM or the BEM crack growth is generally simulated in a stepwise manner. First the stress intensity factors for the different crack opening modes are calculated. By introduction of an appropriate mixed-mode criterion, e.g. the maximum tangential stress criterion [10], which is used in FRANC2D and FRANC3D and in its three-dimensional extension [11,12] in ADAPCRACK3D, the ⇑ Corresponding author. E-mail address: [email protected] (P. Zerres). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2012.01.018
crack growth direction is predicted. With a prescribed amount of crack advance the position of the new crack tip is determined. Knowing the new tip, the crack is enlarged by introduction of new crack faces. This model change requires a remeshing of the structure after the crack advance. The result of a crack growth analysis is the relationship between crack length and stress intensity factor. Taking these values as an input, the fatigue life is calculated by integrating an appropriate crack growth law. The described procedure is similar in all afore mentioned tools. However, they differ in the remeshing strategy, the stress intensity factor calculation and the crack growth direction determination. It should be noted that other approaches exist, which overcome the difficulties of the adaptive remeshing, e.g. the extended finite element method [13] and the meshless Galerkin-method [14]. Focusing on load sequence effects or mean stress effects or assuming the clear violation of small scale yielding, methods based on elastic-plastic fracture mechanics are favored. To gain substantial knowledge of the mechanisms of fatigue crack growth the elastic–plastic material behavior has to be taken into consideration. Thereby the crack path is commonly prescribed and it is assumed to be straight-lined. Fatigue crack growth is generally modeled by a nodal-release scheme. Beginning with the work of Newman [15], who used a nodal-release scheme to identify the influence of the mesh density on crack closure, many researches used nodal-release schemes to investigate the role of plasticity induced crack closure [16] on fatigue crack growth. Within these studies, the following four aspects were found to be essential: the finite element discretization, the material model,
74
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
the definition of a closed crack and the crack growth scheme. Concerning the mesh density, the size of the plastic zone was often used as a characterizing value [17–20], without having a uniform postulat, rather than having contradictions. Jiang et al. [21] found that by using an advanced material model, crack closure is to some extend independent of the mesh size. They also found that the results are strongly mesh dependent when using a simple material model, e.g. linear elastic-ideal plastic. Similar observations can be found in [18,22–24]. A further point, which influences the crack closure behavior is the definition of a closed crack. In most of the investigations found in the literature, a crack is defined as closed, if the first node behind the crack tip is in contact with its opposite crack face [17,18,25]. In other publications, the second node behind the crack tip [26–28] or even the whole crack front [21,29] is investigated. A summary of different crack closure definitions can be found in [30]. Within the nodal-release analysis a point has to be defined, when the nodes at the crack tip are released, i.e. the crack is growing. Again, different estimations are found in the literature: Releasing the nodes at maximum load [31,21], at minimum load [32,24], during unloading [29] and if the reaction force at the crack tip node is equal to zero [33]. However, investigations show that all estimations lead to similar results [34,25]. Recent surveys on the finite element simulation of plasticity induced crack closure were given from Solanki et al. [35], Antunes [30], Jiang et al. [21] and Herz et al. [36]. A drawback of the elastic–plastic methods is that the actual fatigue life is calculated after the simulation by integrating an empirically determined crack growth law. In general the crack path is not known a priori and has to be determined during the analysis of a crack growing in a complex structure under a high cyclic loading level. As in addition the crack growth rate is affected by plasticity effects, both above mentioned methods, used separately from each other, are not qualified solving the problem. Keeping this in mind this paper aims to introduce a procedure to simulate fatigue crack growth by combining these two methods. The identification of the crack growth is thereby including elastic–plastic material behavior. During the simulation the crack path could be determined by remeshing the structure after every increment of crack growth.
The basic inputs are the geometry of the structure, information of the material in the form of the parameters for the used materialmodel as well as the load sequence. Fig. 1 shows the steps of the numerical simulation of the fatigue crack growth. After the development of the finite element model, including a mesh refinement at the crack tip due to high stress gradients, the structure could be analyzed. Within the postprocessing the crack opening and closure level is determined and the effective range of the crack tip parameter is calculated. With this result the crack growth direction and by integrating an appropriate crack growth law the number of cycles to reach a given crack advance (or the crack advance for a prescribed number of cycles) could be determined. A new model is generated based on an update of the geometry through extending the crack in the calculated direction by the crack increment. The displacements and the status variables have to be mapped from the old to the new mesh for reestablishing the stress state prior to crack advance. The isoparametric element formulations are used to transfer the displacements. According to the isoparametric element formulation the coordinates x and the displacements u for an arbitrary point within an element (of the old mesh) can be determined by using the nE nodal coordinates xold and displacements uold j j , as well as the isoparametric shape functions N j :
x¼ u¼
nE X j¼1 nE X
N j xold j ;
ð1Þ
Nj uold j :
ð2Þ
j¼1
The isoparametric shape functions take values equal to one at node j and zero at all nodes i with i–j and are functions of the local element coordinates g 1 and g 2 :
Nj ¼ f ðgÞ;
g¼
g1 g2
In this section the developed crack growth procedure, which is named mapping analysis in the following, is presented, see Fig. 1. All the modules of the procedure are implemented as user routines to the commercial finite element program ABAQUS [37].
ð3Þ
In the example in Section 3 only quadrilateral elements with linear shape functions are used. For this type of element the local coordinates meet the restriction:
1:0 6 g 1 ; g 2 6 1:0:
2. Crack growth procedure
:
ð4Þ
To determine, if a node of the new mesh lies within an element in the old mesh, its coordinates xnew are inserted into Eq. (1). Afterwards Eq. (1) is rearranged and solved via Newtons-method for the local coordinates g 1old and g 2old :
Fig. 1. Crack growth procedure (from [50]).
75
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
F¼
nE X
N j xold xnew ¼ 0: j
ð5Þ
j¼1
If Eq. (4) is not fulfilled by g old and g old 1 2 , the node of the new mesh is not within this element and Eq. (5) has to be solved for the next element in the old mesh. But if Eq. (4) is fulfilled, the new node is within this element, and by inserting the local coordinates in Eq. (2) the displacement of this node can be calculated. By using a routine, which is implemented in ABAQUS by default all the state variables, i.e. the stresses, the strains, the components of the backstress-tensors and others, which depend on the applied material model, are mapped afterwards from the old mesh to the new one. Then the structure is analyzed again to calculate the next increment of crack advance. Due to the mapping of the status variables at the beginning of the analysis illegal conditions, e.g. that the stress tensor is outside the yield surface, are possibly prevailing. The applied material model has to be modified to take these circumstances into account. The procedure is repeated until the crack length reaches a user defined critical value. It should be noted that without the mapping of the variables, i.e. for linear-elastic material behavior, the procedure was successfully validated in [38] for the determination of the crack growth path under mixed-mode loading.
Fig. 2. Finite element models for the nodal-release analysis.
3. Example 3.1. Introduction The examined problem is a pressurized tube with closed ends made of S 460 nbk. The problem was investigated by Herz et al. [39,29], where a nodal release scheme and an advanced material model, which was originally introduced by Döring et al. [40,41] (see Section 3.3), was used to calculate the endurance limit. To increase the fatigue limit the tubes have been autofrettaged. The initial high overload of the autofrettage process introduces beneficial compressive residual stresses in the fatigue critical areas of the structure. For the calculations presented below an autofrettage pressure of paf = 4700 bar is investigated. 3.2. Geometry and mesh For pressurized tubes with closed ends, neither the plane stress nor the plane strain assumption correctly renders the stress state in the wall of the tube. But to avoid the difficulties and computational costs of a three-dimensional analysis, the straight pressurized tube with closed ends is modeled as a completely axis symmetric pressurized torus with the same cross section and a very large radius of 600 mm. This radius was found [39,29] to be sufficient, i.e. the stress state in the tube and the torus are approximately identical. The outer diameter of the cross section of the tube (see Figs. 2 and 3) is do = 6.0 mm and the inner diameter is di = 2.5 mm. Axis symmetric elements with linear shape functions and full integration are employed. For the nodal-release analysis, half of the cross section is modeled and symmetry boundary conditions are introduced for x2 = 0, see Fig. 2. To account for crack face contact a rigid surface is introduced in the x1–x3-plane. Based on convergence studies [29], the element size along the prescribed crack front is set of 2.5 lm. As shown in Fig. 3, for the mapping analysis the whole cross section area is modeled. Here, the crack flank contact is accomplished by introducing contact definitions between the two respective surfaces. The element size at the crack tip is set to 0.625 lm (convergence studies can be found in [42]). By using an advancing front algorithm and a smooth transition from fine to coarse mesh, a good overall mesh quality is ensured. The initial crack length is set to a0 = 10 lm for both types of
Fig. 3. Finite element models for the mapping analysis.
analysis. This initial crack length is in accordance with measurements of defects found in the specimens due to the manufacturing process [39,29]. Despite of the fact that crack advance is assumed to be only in x1-direction, by using the model in Fig. 3 the possibility of an arbitrary crack advance should be demonstrated. 3.3. Material model In this example a problem is investigated, where a high residual stress field is superposed to alternate stress due to cyclic loading. Hence stress relaxation effects and ratcheting will occur and have to be realistically described by the material model. For this reason also in this study the Döring-model, which basic equations are described in the following, is applied. The Mises yield function and the normality rules are the fundaments of the Döring-model:
Uðs; a; rð0Þ Þ ¼ kðs aÞk r ð0Þ ; depl ¼ dp n; dp ¼ kdepl k :
with n ¼
sa ks ak
ð6Þ ;
ð7Þ ð8Þ
with r ð0Þ beeing the yield stress, a the backstress-tensor, s the deviatoric stress-tensor and p the accumulated plastic strain.
76
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
The backstress-tensor is – following the suggestion of Chaboche [43] – decomposed into M parts:
a¼
M X
ðiÞ
a :
ð9Þ
i¼1
and for each part aðiÞ the kinematic hardening rule is expressed as
2 ðiÞ
da
ðiÞ !vðiÞ a
¼ c 4r ðiÞ n ðiÞ
jrðiÞ j
3 ðiÞ 5
a
dp þ aðiÞ
dr ðiÞ ; rðiÞ
ð10Þ
with
aðiÞ vðiÞ ¼ v0ðiÞ þ v0ðiÞ þ 0:1 cðiÞ v 1 1 n : kaðiÞ k
v0ðiÞ ¼ Q ðiÞ 1 þ
and
ð11Þ
!
av
ð12Þ
;
ð1 þ bv RMe Þ2
In Eqs. (10)–(12) the hardening behavior is controlled by cðiÞ and rðiÞ and the ratcheting behavior by vðiÞ . The differential equations for the description of isotropic hardening are formulated as follows: 8 ðiÞ > < b rTðiÞ rðiÞ dp rT rðiÞ 6 1 ðiÞ dr ¼ ; i ¼ 0; 1; . . . ; M: ðiÞ > : b rTðiÞ rðiÞ rTðiÞ rðiÞ dp rT rðiÞ > 1
ð13Þ ðiÞ rT
b is an intern variable of the model and the target-functions describe the relation between r ðiÞ and other state variables by: ! ðiÞ ðiÞ ðiÞ a1 a2 a3 ðiÞ ðiÞ rT ¼ r1 1 þ þ þ ; i ¼ 0; 1; . . . ; M; ð1 þ b1 pÞ2 ð1 þ b2 pÞ2 ð1 þ b3 pÞ2
r ð0Þ 1 r ðiÞ 1
ð14Þ
h i ð0Þ ¼ A qN ðRMe Þ þ ð1 AÞ qP ðRMe Þ ; ! ðiÞ ak ðiÞ ðiÞ ¼ r 1;0 qP ðRMe Þ 1 þ ; i ¼ 1; 2; . . . ; M; 2 ð1 þ bk vðiÞ Þ ð0Þ r1;0
ð15Þ ð16Þ
ðiÞ
aP
ðiÞ
qP ¼ 1 þ
ð1 þ bP RMe Þ2
qN ¼ qN;0 1 þ
i ¼ 0; 1; . . . ; M;
;
ð17Þ
!
aN ð1 þ bN RMe Þ2
ð18Þ
: ðiÞ
The cyclic hardening and softening, which is described by qP (for proportional hardening) and qN (for non-proportional hardening), ðiÞ
ðiÞ
ðiÞ
is controlled by the parameters a1 ; a2 ; a3 ; b1 ; b2 , and b3. According to Tanaka [44] the strength of the non-proportional hardening is determined by the non-proportional-factor:
dA ¼ cA ðAT AÞdp;
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðC : nÞ : ðC : nÞ ; AT ¼ 1 kCk2
dC ¼ cT ðn n CÞdp;
ð19Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðepl bÞ : ðepl bÞ RMe 6 0;
ð21Þ
with its normal
n ¼
epl b RMe
:
db ¼ ð1 gÞ HðU Þ hn : n in dp;
ð23Þ
with H() beeing the Heavyside function and <. . .> beeing the MacCauley-brackets. Following the proposal of Jiang and Kurath [45], the evolution of the radius RMe of the memory surface is described by:
dRMe ¼ ½g HðU Þ ð1 HðU ÞÞ ðcMe RMe Þx hn : n idp:
ð24Þ
Three more material parameters, cMe ; g and x, are introduced in Eq. (24). The material parameters for the investigated material S 460 nbk are shown in Table 1. Additionally the isotropic elastic constants are equal to E = 208.5 GPa and m = 0.3. As stated before, the mapping of the variables can result in illegal conditions in the structure at the beginning of the analysis. The required improvements of the Döring-model are described in [42,50]. 3.4. Load sequence and crack growth The crack growth for the nodal-release and the mapping analyses are akin and subdivided into two phases, see Fig. 4. The first phase is the autofrettage process. There, the specimen, including an initial crack of a0 = 10 lm, is loaded up to autofrettage pressure and then unloaded to pmin = 50 bar. Afterwards 250 load cycles with the operating pressure for the redistribution of the residual stress field due to the combined effect of ratcheting and mean stress relaxation are applied. The second phase is the crack growth phase. First the specimen is loaded up to operating pressure. For the nodal-release analysis the specimen is then unloaded to 90% of the operating pressure and in the same time the boundary condition of fixed displacements in x2-direction (according to Fig. 2) of the node representing the actual crack tip is released and by this a crack advance of one element, which is Da = 2.5 lm in this example, is realized. For the mapping analysis the specimen is also unloaded to 90% of the operating pressure. At this stage the crack is enlarged (Da = 2.5 lm for the comparison with the nodal-release analysis and the endurance limit calculation and Da = 15% of the actual crack length for the fatigue life calculation), the new model is generated and the mapping of the variables is performed. Then one step with a constant load of 90% of the operating pressure is performed in the new model in order to establish equilibrium within the structure. Afterwards—for both types of analysis—the specimen is further unloaded up to pmin. Between two crack growth steps, 10 additional load cycles are applied to stabilize the stress state at the crack tip. Within the last of these cycles the values needed for the calculation of the crack advance, which are the displacements in x2-direction for the nodal-release and the crack opening for the mapping analysis, are read out. This cycle is subdivided into at
ð20Þ
with C being the Tanaka-tensor of fourth order. Chaboche [43]suggested a strain based memory surface to improve the accuracy of modeling the ratcheting behavior. Taking this proposal into account Eq. (21) shows the definition of the memory surface
U ¼
b, the center of memory surface, and its movement is described by
ð22Þ
Table 1 Material parameters for the Döring-model for S 460 nbk [42].
g 0.1
cMe 5 ð0Þ
rinit (MPa) 408
x 1
cT 50
cA 80
qN,0 1.74
bN 100
b1
b2
b3
bP
av
bv
100
10
0.3
1.0
250
ðiÞ
0.7 Q(i) — 18 5 3 3 18
— 50 30 30 30 1
i
c(i)
r 1;0 (MPa)
a1
a2
a3
ðiÞ
aP
0 1 2 3 4 5
— 1260 408 126 40.8 12.6
201 69.4 49.6 69.9 79.2 163
1.1 1.0 0.7 0.0 0.0 0.0
0.0 0.4 1.0 0.3 0.0 0.0
0.1 0.7 1.2 0.5 0.5 0.6
0 0 0 0 0 0
ðiÞ
aN 0.26
ðiÞ
ðiÞ
ðiÞ
cv
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
Fig. 4. Crack growth calculation scheme.
least 20 increments for both loading and unloading in order to obtain enough values for an accurate calculation. In this study the crack opening point is defined as the point, where the displacements in x2-direction or the crack openings at all nodes on the crack front are greater than zero. This definition is in accordance with [39,29]. By knowing the crack opening point the effective pressure range D peff can be calculated and by this the effective range of the stress intensity factor DKeff can be determined by integrating a weight function solution presented by Ma et al. [46]. 3.5. Results and discussion 3.5.1. Comparison of the two types of analysis In this section the results obtained with the nodal-release analysis and with the mapping analysis are compared. For this purpose, the crack advance per increment for the mapping analysis is also set to D a = 2.5 lm. In Fig. 5 the calculated crack opening pressures for the two types of analysis for an operating pressure of Dp = 2690 bar are depicted. It becomes obvious that the results for both types of analysis are very similar, which leads to the conclusion that the inaccuracies due to the mapping of the variables are very small.
77
Unlike in the previous section the crack advance per increment is set to 15% of the actual crack length. A reason for this choice is given in the next section. Fig. 6 shows the results of these calculations. For an operating pressure of Dp = 2950 bar DKeff is always significantly higher than the threshold value and no crack arrest is predicted. By decreasing the operating pressure it is found that at an pressure of Dp = DpE,calc = 2620 bar DKeff falls slightly short of DKeff,th and therefore crack arrest is predicted. By comparing the numerical endurance limit DpE,calc with the experimental one DpE, exp,50% the difference is only 2.6%. From this it follows that the proposed procedure is able to predict the endurance limit with a high accuracy. 3.5.3. Fatigue life calculation In this section the fatigue lives of the tube under different operating pressures is calculated. The crack advance per increment for these calculations is also set to 15% of the actual crack length. The choice of this relative crack increment is more appropriate than a constant increment when modeling the crack growth from very small crack length of some microns up to a crack length of one millimeter (or even more). The final crack length is set to 0.75 mm, i.e. 31 crack growth increments are simulated. To calculate the fatigue life an appropriate crack growth law has to be adopted. The one used in these calculations is taken from [48]:
pffiffiffiffiffi da ¼ 8:8 109 mm=cyc ðMPa mÞ3:15 DK 3:15 eff : dn
ð25Þ
By using a Paris-law (Eq. (25)) the crack growth rate nearby the threshold value is overestimated and by this shorter fatigue lives are calculated. For this reason a second crack growth law, which was proposed by Donahue et al. [49] and involves the threshold value, is adopted:
pffiffiffiffiffi da 3:15 ¼ 8:8 109 mm=cyc ðMPa mÞ3:15 DK 3:15 : eff DK eff;th dn ð26Þ
3.5.2. Endurance limit calculation The aim of these calculations is the determination of the endurance limit of the autofrettaged tubes. The experimental endurance limit was found to be DpE,exp,50% = 2690 bar [39,29]. To determine the numerical endurance limit fatigue crack growth analyses for different operating pressures are performed. The results of these analyses are the evolution of D Keff over the crack length. If DKeff is equal or below the threshold value, which is known from previpffiffiffiffiffi ous investigations [47] as DKeff,th = 3 MPa m, crack arrest is predicted and the structure is considered to be fatigue resistant.
In Fig. 7 the evolution of the number of nodes and elements during the crack growth analysis for an operating pressure of Dp = 2950 bar is depicted. The corresponding finite element meshes for crack lengths of about 0.01 mm, 0.16 mm and 0.66 mm at 90% of the operating pressure are shown in Fig. 8. By looking at these figures one big advantage of the mapping analysis, which uses the remeshing technique, becomes obvious. It is possible to model the crack growth from very small crack lengths up to large crack lengths with a moderate number of nodes and elements.
Fig. 5. Crack opening pressure for the nodal-release and the mapping analysis.
Fig. 6. Evolution of DKeff over crack length.
78
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
Fig. 7. Evolution of the number of nodes and elements during the crack growth simulation.
Fig. 9. Calculated and experimental [39,29] fatigue lives.
numerical results is expected. Another factor, which has to be also considered is that the crack growth constants in this computations have been identified in the medium air. In this example the medium oil is existent. This might lead to different constants due to the fact of fluid-induced crack closure [16] and with that to different fatigue lives. 4. Conclusions The numerical crack growth simulation has become more and more attractive today for reasons of time and economy. The presentation of a procedure for such a simulation of fatigue crack growth in elastic–plastic materials is the aim of this paper. The remeshing of the structure and the subsequent mapping of the status variables from the old mesh to the new one after each step of crack advance are thereby the core characteristics of the proposed procedure. For the validation of the procedure results on autofrettaged Diesel-engine injection tubes are used. The results obtained by this procedure using an advanced material model, which is able to realistically simulate cyclic plasticity effects such as mean stress relaxation and ratcheting, are compared with experimental results and finite element results based on a nodal-release scheme. From these results the following conclusions can be drawn:
Fig. 8. Initial mesh (a0 = 0.01 mm) and meshed geometries for crack lengths of a 0.16 mm and a 0.66 mm (paf = 4700 bar, Dp = 2950 bar).
In Fig. 9 the calculated fatigue lives and the experimental ones are shown. As expected, the fatigue lives calculated via Eq. (26) are by a factor of about two longer than the ones calculated via Eq. (25). But they are by a factor of about three to four shorter than the experimental results. The reasons for this underestimation of the fatigue lives are not only due to the mapping of the variables, but also on other factors, which will be discussed in the following. The most important point is that because of the axis symmetric modeling of the structure a straight through the thickness crack is implied. But it is known from the experiments that semi elliptical surface cracks initiated. These cracks have a stress intensity factor, which is (at the deepest position) much lower than the one for a through the thickness crack of the same length. By taking into account these three-dimensional effects, a significant improvement of the
The results obtained by the proposed procedure are in good agreement with the results obtained by the nodal-release scheme. The results obtained by the proposed procedure are in good agreement with the experimental results concerning the endurance limit. The fatigue lives are underestimated by the proposed procedure. This is not only due to the inaccuracies resulting from the mapping of the variables but also due to the simplified axis symmetric modeling and the choice of the crack growth constants. A great advantage of this procedure is that arbitrary loading conditions, i.e. proportional as well as non-proportional loading, can be examined. But further validations of the procedure are required. Also an extension of the implementation of the procedure to account for three-dimensional structures is under progress. References [1] D.V. Swenson, A.R. Ingraffea, Computational Mechanics 3 (1988) 381–397. [2] T.N. Bittencourt, P.A. Wawrzynek, A.R. Ingraffea, J.L. Sousa, Engineering Fracture Mechanics 55 (1996) 321–334.
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79 [3] L.F. Martha, Topological and Geometrical Modeling Approach to Numerical Discretization and Arbitrary Fracture Simulation in Three Dimensions, Ph.D. Thesis, School of Civil Engineering, Cornell University, 1989. [4] B.J. Carter, P.A. Wawrzynek, A.R. Ingraffea, Gallagher Special Issue of International Journal for Numerical Methods in Engineering 47 (2000) 229– 253. [5] Zentech International Ltd., Zencrack. Tool for 3D Fracture Mechanics Simulation. [6] C. Timbrell, G. Cook, 3D FE fracture mechanics analysis for industrial applications, in: Seminar on Inelastic Finite Element Analysis, Institute of Mechanical Engineers, London, 1997. [7] M. Fulland, Risssimulation in dreidimensionalen Strukturen mit automatischer adaptiver Finite-Elemente-Netzgenerierung, Fortschrittbericht-Bericht VDI, Reihe 18: Mechanik/Bruchmechanik Nr. 280, VDI-Verlag, Dnsseldorf, 2003. [8] M. Schöllmann, M. Fulland, H.A. Richard, Engineering Fracture Mechanics 70 (2003) 249–268. [9] S. Mellings, J. Baynham, R.A. Adey, Engineering Analysis with Boundary Elements 26 (2002) 479–488. [10] F. Erdogan, G.C. Sih, Journal of Basic Engineering 85 (1963) 519–527. [11] M. Schöllmann, G. Kullmer, M. Fulland, H.A. Richard, A new criterion for 3D crack growth under mixed-mode (I + II + III) loading, in: M. de Freitas (Ed.), Proceedings of the 6th International Conference on Biaxial/Multiaxial Fatigue and Fracture, 2001, pp. 589–596. [12] M. Schöllmann, H.A. Richard, G. Kullmer, M. Fulland, International Journal of Fracture 117 (2002) 129–141. [13] N. Moes, J. Dolbow, T. Belytschko, International Journal for Numerical Methods in Engineering 46 (1999) 131–150. [14] T. Belytschko, Y.Y. Lu, L. Gu, International Journal for Numerical Methods in Engineering 37 (1994) 229–256. [15] J.C. Newman Jr., A crack-closure model for predicting fatigue crack growth under aircraft spectrum loading, in: J.B. Chang, C.M. Hudson (Eds.), Methods and Models for Predicting Fatigue Crack Growth Under Random Loading, ASTM STP 748, 1981, pp. 53–84. [16] W. Elber, Engineering Fracture Mechanics 2 (1970) 37–45. [17] R.C. McClung, H. Sehitoglu, Engineering Fracture Mechanics 33 (1989) 237– 252. [18] R.C. McClung, H. Sehitoglu, Engineering Fracture Mechanics 33 (1989) 253– 272. [19] J.D. Dougherty, J. Padovan, T.S. Srivatsan, Engineering Fracture Mechanics 56 (1997) 189–212. [20] A. Gonzales-Herrera, J. Zapatero, Engineering Fracture Mechanics 72 (2005) 337–355. [21] Y. Jiang, M. Feng, F. Ding, International Journal of Plasticity 21 (2005) 1720– 1740. [22] N.E. Ashbaugh, B. Dattaguru, M. Khobaib, T. Nicholas, R.V. Prakash, T.S. Ramamurthy, B.R. Seshadri, R. Sunder, Fatigue & Fracture of Engineering Materials & Structures 20 (1997) 963–974. [23] S. Pommier, P. Bompard, Fatigue & Fracture of Engineering Materials & Structures 23 (2000) 129–139. [24] D.M. Rodrigues, F.V. Antunes, Engineering Fracture Mechanics 76 (2009) 1215–1230. [25] K. Solanki, S.R. Daniewicz, J.C. Newman Jr., Engineering Fracture Mechanics 70 (2003) 1475–1489. [26] K. Solanki, S.R. Daniewicz, J.C. Newman Jr., Engineering Fracture Mechanics 71 (2004) 1165–1175. [27] S. Roychowdhury, R.H. Dodds Jr., Fatigue & Fracture of Engineering Materials & Structures 26 (2003) 663–673.
79
[28] S. Roychowdhury, R.H. Dodds Jr., Engineering Fracture Mechanics 70 (2003) 2363–2383. [29] E. Herz, R. Thumser, J.W. Bergmann, M. Vormwald, Engineering Fracture Mechanics 73 (2006) 3–21. [30] F.V. Antunes, D.M. Rodrigues, Engineering Fracture Mechanics 75 (2008) 3101–3120. [31] J.C. Newman, Jr., A finite-element analysis of fatigue crack closure, in: J.R. Rice, P.C. Paris (Eds.), Mechanics of Crack Growth, ASTM STP 590, 1976, pp. 281– 301. [32] S. Pommier, International Journal of Fatigue 23 (2001) 111–118. [33] K. Ogura, K. Ohji, K. Honda, Influence of mechanical factors on the fatigue crack closure, in: D.M.R. Taplin (Ed.), Fracture, 1977 (vol. 2), ICF4, Waterloo, Canada, 1977, pp. 1035–1047. [34] R.C. McClung, B.H. Thacker, S. Roy, International Journal of Fracture 50 (1991) 27–49. [35] K. Solanki, S.R. Daniewicz, J.C. Newman Jr., Engineering Fracture Mechanics 71 (2004) 149–171. [36] E. Herz, O. Hertel, M. Vormwald, Engineering Fracture Mechanics 78 (2011) 559–572. [37] Dassault Systémes Simulia Corp., Abaqus Analysis User’s Manual (6.8ef), 2008. [38] P. Zerres, M. Vormwald, Modeling of crack propagation under mixed-mode loading with advanced remeshing, in: Proceedings of the 12th International Conference on Fracture, Ottawa, Kanada, 2009. [39] E. Herz, R. Thumser, M. Vormwald, J. Bergmann, Endurance limit design of high pressure Diesel injection pipes for low failure probability, in: Proceedings of Fatigue Design, Senlis, 2005, pp. 1–8. [40] R. Döring, J. Hoffmeyer, T. Seeger, Computational Materials Science 28 (2003) 587–596. [41] R. Döring, Zum Deformations- und Schädigungsverhalten metallischer Werkstoffe unter mehrachsig nichtproportionalen zyklischen Beanspruchungen, Ph.D. Thesis, Institut für Stahlbau und Werkstoffmechanik, TU Darmstadt, Heft 78, 2006. [42] P. Zerres, Numerische Simulation des Ermüdungsrissfortschrittes in metallischen Strukturen unter Berücksichtigung zyklischer Plastizitätseffekte, Ph.D. Thesis, Institut für Stahlbau und Werkstoffmechanik, TU Darmstadt, Heft 88, 2010. [43] J.L. Chaboche, K. Dang Van, G. Cordier, Modelization of the strain memory effect on the cyclic hardening of 316 stainless steel, in: Transactions of the 5th International Conference on Structural Mechanics in Reactor Technology, Div. L, 1979, p. L11/3. [44] E. Tanaka, European Journal of Mechanics, A/Solids 13 (1994) 155–173. [45] Y. Jiang, P. Kurath, Journal of Engineering Materials and Technology 119 (1997) 161–170. [46] C.-C. Ma, J.-I. Huang, C.-H. Tsai, Journal of Pressure Vessels Technology 116 (1994) 423–430. [47] Windte, K., Experimentelle Bestimmung des Schwellwertes fnr Risswachstum unter zyklischer Beanspruchung zur Verbesserung der rechnerischen Lebensdauervorhersage bei Betriebsbelastung, 1986. [48] M. Vormwald, Anrisslebensdauervorhersage auf der Basis der Schwingbruchmechanik für kurze Risse, Ph.D. Thesis, Institut für Stahlbau und Werkstoffmechanik, TH Darmstadt, Heft 47, 1989. [49] R.J. Donahue, H. McI Clark, P. Atanmo, R. Kumble, A.J. McEvily, Crack opening displacement and the rate of fatigue crack growth, International Journal of Fracture 8 (1972) 209–219. [50] P. Zerres, R. Döring, M. Vormwald, Materials Science and Engineering Technology 42 (2011) 1093–1098.
Contents lists available at SciVerse ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Finite element based simulation of fatigue crack growth with a focus on elastic–plastic material behavior P. Zerres ⇑, M. Vormwald Materials Mechanics Group, Technische Universität Darmstadt, Petersenstr. 12, D-64287 Darmstadt, Germany
a r t i c l e
i n f o
Article history: Received 26 November 2010 Received in revised form 23 December 2011 Accepted 13 January 2012
Keywords: Fatigue crack growth Crack closure Finite element method
a b s t r a c t The numerical crack growth simulation has become more and more attractive today for reasons of time and economy. The presentation of a procedure for such a simulation of fatigue crack growth in elastic– plastic materials is the aim of this paper. The remeshing of the structure and the subsequent mapping of the status variables from the old mesh to the new one after each step of crack advance are thereby the core characteristics of the proposed procedure. The crack growth life is determined by integration of a crack growth law which is based on the effective range of the stress intensity factor. The results obtained by the procedure for autofrettaged Diesel-engine injection tubes are compared with experimental results and finite element results based on a nodal release scheme. A good agreement between the two numerical analyses and between the results with the proposed procedure and the experimental results concerning the endurance limit can be stated. By looking at the fatigue lives the experimental results are underestimated by the numerical ones. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction For calculating the fatigue crack growth, which often is responsible for failure in structures under predominant cyclic loading, reliable and in research established methods differ in numerical crack growth simulations based either on linear-elastic or on elastic–plastic fracture mechanics. Choosing the first approach the assumption of small scale yielding is necessary. Research in this field focus on the description of the crack path. The range of the stress intensity factor is the fundamental parameter together with its connection to the crack growth rate, which has to be determined empirically. Two methods are successful in use today to calculate the stress intensity factors, the finite element method (FEM) and the boundary element method (BEM). The crack growth simulation based on linear-elastic fracture mechanics is implemented in several two- and three-dimensional, commercial and non-commercial software tools, such as FRANC2D [1,2], FRANC3D [3,4], ZEN CRACK [5,6], ADAPCRACK3D [7,8] and BEASY [9]. Using the FEM or the BEM crack growth is generally simulated in a stepwise manner. First the stress intensity factors for the different crack opening modes are calculated. By introduction of an appropriate mixed-mode criterion, e.g. the maximum tangential stress criterion [10], which is used in FRANC2D and FRANC3D and in its three-dimensional extension [11,12] in ADAPCRACK3D, the ⇑ Corresponding author. E-mail address: [email protected] (P. Zerres). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2012.01.018
crack growth direction is predicted. With a prescribed amount of crack advance the position of the new crack tip is determined. Knowing the new tip, the crack is enlarged by introduction of new crack faces. This model change requires a remeshing of the structure after the crack advance. The result of a crack growth analysis is the relationship between crack length and stress intensity factor. Taking these values as an input, the fatigue life is calculated by integrating an appropriate crack growth law. The described procedure is similar in all afore mentioned tools. However, they differ in the remeshing strategy, the stress intensity factor calculation and the crack growth direction determination. It should be noted that other approaches exist, which overcome the difficulties of the adaptive remeshing, e.g. the extended finite element method [13] and the meshless Galerkin-method [14]. Focusing on load sequence effects or mean stress effects or assuming the clear violation of small scale yielding, methods based on elastic-plastic fracture mechanics are favored. To gain substantial knowledge of the mechanisms of fatigue crack growth the elastic–plastic material behavior has to be taken into consideration. Thereby the crack path is commonly prescribed and it is assumed to be straight-lined. Fatigue crack growth is generally modeled by a nodal-release scheme. Beginning with the work of Newman [15], who used a nodal-release scheme to identify the influence of the mesh density on crack closure, many researches used nodal-release schemes to investigate the role of plasticity induced crack closure [16] on fatigue crack growth. Within these studies, the following four aspects were found to be essential: the finite element discretization, the material model,
74
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
the definition of a closed crack and the crack growth scheme. Concerning the mesh density, the size of the plastic zone was often used as a characterizing value [17–20], without having a uniform postulat, rather than having contradictions. Jiang et al. [21] found that by using an advanced material model, crack closure is to some extend independent of the mesh size. They also found that the results are strongly mesh dependent when using a simple material model, e.g. linear elastic-ideal plastic. Similar observations can be found in [18,22–24]. A further point, which influences the crack closure behavior is the definition of a closed crack. In most of the investigations found in the literature, a crack is defined as closed, if the first node behind the crack tip is in contact with its opposite crack face [17,18,25]. In other publications, the second node behind the crack tip [26–28] or even the whole crack front [21,29] is investigated. A summary of different crack closure definitions can be found in [30]. Within the nodal-release analysis a point has to be defined, when the nodes at the crack tip are released, i.e. the crack is growing. Again, different estimations are found in the literature: Releasing the nodes at maximum load [31,21], at minimum load [32,24], during unloading [29] and if the reaction force at the crack tip node is equal to zero [33]. However, investigations show that all estimations lead to similar results [34,25]. Recent surveys on the finite element simulation of plasticity induced crack closure were given from Solanki et al. [35], Antunes [30], Jiang et al. [21] and Herz et al. [36]. A drawback of the elastic–plastic methods is that the actual fatigue life is calculated after the simulation by integrating an empirically determined crack growth law. In general the crack path is not known a priori and has to be determined during the analysis of a crack growing in a complex structure under a high cyclic loading level. As in addition the crack growth rate is affected by plasticity effects, both above mentioned methods, used separately from each other, are not qualified solving the problem. Keeping this in mind this paper aims to introduce a procedure to simulate fatigue crack growth by combining these two methods. The identification of the crack growth is thereby including elastic–plastic material behavior. During the simulation the crack path could be determined by remeshing the structure after every increment of crack growth.
The basic inputs are the geometry of the structure, information of the material in the form of the parameters for the used materialmodel as well as the load sequence. Fig. 1 shows the steps of the numerical simulation of the fatigue crack growth. After the development of the finite element model, including a mesh refinement at the crack tip due to high stress gradients, the structure could be analyzed. Within the postprocessing the crack opening and closure level is determined and the effective range of the crack tip parameter is calculated. With this result the crack growth direction and by integrating an appropriate crack growth law the number of cycles to reach a given crack advance (or the crack advance for a prescribed number of cycles) could be determined. A new model is generated based on an update of the geometry through extending the crack in the calculated direction by the crack increment. The displacements and the status variables have to be mapped from the old to the new mesh for reestablishing the stress state prior to crack advance. The isoparametric element formulations are used to transfer the displacements. According to the isoparametric element formulation the coordinates x and the displacements u for an arbitrary point within an element (of the old mesh) can be determined by using the nE nodal coordinates xold and displacements uold j j , as well as the isoparametric shape functions N j :
x¼ u¼
nE X j¼1 nE X
N j xold j ;
ð1Þ
Nj uold j :
ð2Þ
j¼1
The isoparametric shape functions take values equal to one at node j and zero at all nodes i with i–j and are functions of the local element coordinates g 1 and g 2 :
Nj ¼ f ðgÞ;
g¼
g1 g2
In this section the developed crack growth procedure, which is named mapping analysis in the following, is presented, see Fig. 1. All the modules of the procedure are implemented as user routines to the commercial finite element program ABAQUS [37].
ð3Þ
In the example in Section 3 only quadrilateral elements with linear shape functions are used. For this type of element the local coordinates meet the restriction:
1:0 6 g 1 ; g 2 6 1:0:
2. Crack growth procedure
:
ð4Þ
To determine, if a node of the new mesh lies within an element in the old mesh, its coordinates xnew are inserted into Eq. (1). Afterwards Eq. (1) is rearranged and solved via Newtons-method for the local coordinates g 1old and g 2old :
Fig. 1. Crack growth procedure (from [50]).
75
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
F¼
nE X
N j xold xnew ¼ 0: j
ð5Þ
j¼1
If Eq. (4) is not fulfilled by g old and g old 1 2 , the node of the new mesh is not within this element and Eq. (5) has to be solved for the next element in the old mesh. But if Eq. (4) is fulfilled, the new node is within this element, and by inserting the local coordinates in Eq. (2) the displacement of this node can be calculated. By using a routine, which is implemented in ABAQUS by default all the state variables, i.e. the stresses, the strains, the components of the backstress-tensors and others, which depend on the applied material model, are mapped afterwards from the old mesh to the new one. Then the structure is analyzed again to calculate the next increment of crack advance. Due to the mapping of the status variables at the beginning of the analysis illegal conditions, e.g. that the stress tensor is outside the yield surface, are possibly prevailing. The applied material model has to be modified to take these circumstances into account. The procedure is repeated until the crack length reaches a user defined critical value. It should be noted that without the mapping of the variables, i.e. for linear-elastic material behavior, the procedure was successfully validated in [38] for the determination of the crack growth path under mixed-mode loading.
Fig. 2. Finite element models for the nodal-release analysis.
3. Example 3.1. Introduction The examined problem is a pressurized tube with closed ends made of S 460 nbk. The problem was investigated by Herz et al. [39,29], where a nodal release scheme and an advanced material model, which was originally introduced by Döring et al. [40,41] (see Section 3.3), was used to calculate the endurance limit. To increase the fatigue limit the tubes have been autofrettaged. The initial high overload of the autofrettage process introduces beneficial compressive residual stresses in the fatigue critical areas of the structure. For the calculations presented below an autofrettage pressure of paf = 4700 bar is investigated. 3.2. Geometry and mesh For pressurized tubes with closed ends, neither the plane stress nor the plane strain assumption correctly renders the stress state in the wall of the tube. But to avoid the difficulties and computational costs of a three-dimensional analysis, the straight pressurized tube with closed ends is modeled as a completely axis symmetric pressurized torus with the same cross section and a very large radius of 600 mm. This radius was found [39,29] to be sufficient, i.e. the stress state in the tube and the torus are approximately identical. The outer diameter of the cross section of the tube (see Figs. 2 and 3) is do = 6.0 mm and the inner diameter is di = 2.5 mm. Axis symmetric elements with linear shape functions and full integration are employed. For the nodal-release analysis, half of the cross section is modeled and symmetry boundary conditions are introduced for x2 = 0, see Fig. 2. To account for crack face contact a rigid surface is introduced in the x1–x3-plane. Based on convergence studies [29], the element size along the prescribed crack front is set of 2.5 lm. As shown in Fig. 3, for the mapping analysis the whole cross section area is modeled. Here, the crack flank contact is accomplished by introducing contact definitions between the two respective surfaces. The element size at the crack tip is set to 0.625 lm (convergence studies can be found in [42]). By using an advancing front algorithm and a smooth transition from fine to coarse mesh, a good overall mesh quality is ensured. The initial crack length is set to a0 = 10 lm for both types of
Fig. 3. Finite element models for the mapping analysis.
analysis. This initial crack length is in accordance with measurements of defects found in the specimens due to the manufacturing process [39,29]. Despite of the fact that crack advance is assumed to be only in x1-direction, by using the model in Fig. 3 the possibility of an arbitrary crack advance should be demonstrated. 3.3. Material model In this example a problem is investigated, where a high residual stress field is superposed to alternate stress due to cyclic loading. Hence stress relaxation effects and ratcheting will occur and have to be realistically described by the material model. For this reason also in this study the Döring-model, which basic equations are described in the following, is applied. The Mises yield function and the normality rules are the fundaments of the Döring-model:
Uðs; a; rð0Þ Þ ¼ kðs aÞk r ð0Þ ; depl ¼ dp n; dp ¼ kdepl k :
with n ¼
sa ks ak
ð6Þ ;
ð7Þ ð8Þ
with r ð0Þ beeing the yield stress, a the backstress-tensor, s the deviatoric stress-tensor and p the accumulated plastic strain.
76
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
The backstress-tensor is – following the suggestion of Chaboche [43] – decomposed into M parts:
a¼
M X
ðiÞ
a :
ð9Þ
i¼1
and for each part aðiÞ the kinematic hardening rule is expressed as
2 ðiÞ
da
ðiÞ !vðiÞ a
¼ c 4r ðiÞ n ðiÞ
jrðiÞ j
3 ðiÞ 5
a
dp þ aðiÞ
dr ðiÞ ; rðiÞ
ð10Þ
with
aðiÞ vðiÞ ¼ v0ðiÞ þ v0ðiÞ þ 0:1 cðiÞ v 1 1 n : kaðiÞ k
v0ðiÞ ¼ Q ðiÞ 1 þ
and
ð11Þ
!
av
ð12Þ
;
ð1 þ bv RMe Þ2
In Eqs. (10)–(12) the hardening behavior is controlled by cðiÞ and rðiÞ and the ratcheting behavior by vðiÞ . The differential equations for the description of isotropic hardening are formulated as follows: 8 ðiÞ > < b rTðiÞ rðiÞ dp rT rðiÞ 6 1 ðiÞ dr ¼ ; i ¼ 0; 1; . . . ; M: ðiÞ > : b rTðiÞ rðiÞ rTðiÞ rðiÞ dp rT rðiÞ > 1
ð13Þ ðiÞ rT
b is an intern variable of the model and the target-functions describe the relation between r ðiÞ and other state variables by: ! ðiÞ ðiÞ ðiÞ a1 a2 a3 ðiÞ ðiÞ rT ¼ r1 1 þ þ þ ; i ¼ 0; 1; . . . ; M; ð1 þ b1 pÞ2 ð1 þ b2 pÞ2 ð1 þ b3 pÞ2
r ð0Þ 1 r ðiÞ 1
ð14Þ
h i ð0Þ ¼ A qN ðRMe Þ þ ð1 AÞ qP ðRMe Þ ; ! ðiÞ ak ðiÞ ðiÞ ¼ r 1;0 qP ðRMe Þ 1 þ ; i ¼ 1; 2; . . . ; M; 2 ð1 þ bk vðiÞ Þ ð0Þ r1;0
ð15Þ ð16Þ
ðiÞ
aP
ðiÞ
qP ¼ 1 þ
ð1 þ bP RMe Þ2
qN ¼ qN;0 1 þ
i ¼ 0; 1; . . . ; M;
;
ð17Þ
!
aN ð1 þ bN RMe Þ2
ð18Þ
: ðiÞ
The cyclic hardening and softening, which is described by qP (for proportional hardening) and qN (for non-proportional hardening), ðiÞ
ðiÞ
ðiÞ
is controlled by the parameters a1 ; a2 ; a3 ; b1 ; b2 , and b3. According to Tanaka [44] the strength of the non-proportional hardening is determined by the non-proportional-factor:
dA ¼ cA ðAT AÞdp;
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðC : nÞ : ðC : nÞ ; AT ¼ 1 kCk2
dC ¼ cT ðn n CÞdp;
ð19Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðepl bÞ : ðepl bÞ RMe 6 0;
ð21Þ
with its normal
n ¼
epl b RMe
:
db ¼ ð1 gÞ HðU Þ hn : n in dp;
ð23Þ
with H() beeing the Heavyside function and <. . .> beeing the MacCauley-brackets. Following the proposal of Jiang and Kurath [45], the evolution of the radius RMe of the memory surface is described by:
dRMe ¼ ½g HðU Þ ð1 HðU ÞÞ ðcMe RMe Þx hn : n idp:
ð24Þ
Three more material parameters, cMe ; g and x, are introduced in Eq. (24). The material parameters for the investigated material S 460 nbk are shown in Table 1. Additionally the isotropic elastic constants are equal to E = 208.5 GPa and m = 0.3. As stated before, the mapping of the variables can result in illegal conditions in the structure at the beginning of the analysis. The required improvements of the Döring-model are described in [42,50]. 3.4. Load sequence and crack growth The crack growth for the nodal-release and the mapping analyses are akin and subdivided into two phases, see Fig. 4. The first phase is the autofrettage process. There, the specimen, including an initial crack of a0 = 10 lm, is loaded up to autofrettage pressure and then unloaded to pmin = 50 bar. Afterwards 250 load cycles with the operating pressure for the redistribution of the residual stress field due to the combined effect of ratcheting and mean stress relaxation are applied. The second phase is the crack growth phase. First the specimen is loaded up to operating pressure. For the nodal-release analysis the specimen is then unloaded to 90% of the operating pressure and in the same time the boundary condition of fixed displacements in x2-direction (according to Fig. 2) of the node representing the actual crack tip is released and by this a crack advance of one element, which is Da = 2.5 lm in this example, is realized. For the mapping analysis the specimen is also unloaded to 90% of the operating pressure. At this stage the crack is enlarged (Da = 2.5 lm for the comparison with the nodal-release analysis and the endurance limit calculation and Da = 15% of the actual crack length for the fatigue life calculation), the new model is generated and the mapping of the variables is performed. Then one step with a constant load of 90% of the operating pressure is performed in the new model in order to establish equilibrium within the structure. Afterwards—for both types of analysis—the specimen is further unloaded up to pmin. Between two crack growth steps, 10 additional load cycles are applied to stabilize the stress state at the crack tip. Within the last of these cycles the values needed for the calculation of the crack advance, which are the displacements in x2-direction for the nodal-release and the crack opening for the mapping analysis, are read out. This cycle is subdivided into at
ð20Þ
with C being the Tanaka-tensor of fourth order. Chaboche [43]suggested a strain based memory surface to improve the accuracy of modeling the ratcheting behavior. Taking this proposal into account Eq. (21) shows the definition of the memory surface
U ¼
b, the center of memory surface, and its movement is described by
ð22Þ
Table 1 Material parameters for the Döring-model for S 460 nbk [42].
g 0.1
cMe 5 ð0Þ
rinit (MPa) 408
x 1
cT 50
cA 80
qN,0 1.74
bN 100
b1
b2
b3
bP
av
bv
100
10
0.3
1.0
250
ðiÞ
0.7 Q(i) — 18 5 3 3 18
— 50 30 30 30 1
i
c(i)
r 1;0 (MPa)
a1
a2
a3
ðiÞ
aP
0 1 2 3 4 5
— 1260 408 126 40.8 12.6
201 69.4 49.6 69.9 79.2 163
1.1 1.0 0.7 0.0 0.0 0.0
0.0 0.4 1.0 0.3 0.0 0.0
0.1 0.7 1.2 0.5 0.5 0.6
0 0 0 0 0 0
ðiÞ
aN 0.26
ðiÞ
ðiÞ
ðiÞ
cv
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
Fig. 4. Crack growth calculation scheme.
least 20 increments for both loading and unloading in order to obtain enough values for an accurate calculation. In this study the crack opening point is defined as the point, where the displacements in x2-direction or the crack openings at all nodes on the crack front are greater than zero. This definition is in accordance with [39,29]. By knowing the crack opening point the effective pressure range D peff can be calculated and by this the effective range of the stress intensity factor DKeff can be determined by integrating a weight function solution presented by Ma et al. [46]. 3.5. Results and discussion 3.5.1. Comparison of the two types of analysis In this section the results obtained with the nodal-release analysis and with the mapping analysis are compared. For this purpose, the crack advance per increment for the mapping analysis is also set to D a = 2.5 lm. In Fig. 5 the calculated crack opening pressures for the two types of analysis for an operating pressure of Dp = 2690 bar are depicted. It becomes obvious that the results for both types of analysis are very similar, which leads to the conclusion that the inaccuracies due to the mapping of the variables are very small.
77
Unlike in the previous section the crack advance per increment is set to 15% of the actual crack length. A reason for this choice is given in the next section. Fig. 6 shows the results of these calculations. For an operating pressure of Dp = 2950 bar DKeff is always significantly higher than the threshold value and no crack arrest is predicted. By decreasing the operating pressure it is found that at an pressure of Dp = DpE,calc = 2620 bar DKeff falls slightly short of DKeff,th and therefore crack arrest is predicted. By comparing the numerical endurance limit DpE,calc with the experimental one DpE, exp,50% the difference is only 2.6%. From this it follows that the proposed procedure is able to predict the endurance limit with a high accuracy. 3.5.3. Fatigue life calculation In this section the fatigue lives of the tube under different operating pressures is calculated. The crack advance per increment for these calculations is also set to 15% of the actual crack length. The choice of this relative crack increment is more appropriate than a constant increment when modeling the crack growth from very small crack length of some microns up to a crack length of one millimeter (or even more). The final crack length is set to 0.75 mm, i.e. 31 crack growth increments are simulated. To calculate the fatigue life an appropriate crack growth law has to be adopted. The one used in these calculations is taken from [48]:
pffiffiffiffiffi da ¼ 8:8 109 mm=cyc ðMPa mÞ3:15 DK 3:15 eff : dn
ð25Þ
By using a Paris-law (Eq. (25)) the crack growth rate nearby the threshold value is overestimated and by this shorter fatigue lives are calculated. For this reason a second crack growth law, which was proposed by Donahue et al. [49] and involves the threshold value, is adopted:
pffiffiffiffiffi da 3:15 ¼ 8:8 109 mm=cyc ðMPa mÞ3:15 DK 3:15 : eff DK eff;th dn ð26Þ
3.5.2. Endurance limit calculation The aim of these calculations is the determination of the endurance limit of the autofrettaged tubes. The experimental endurance limit was found to be DpE,exp,50% = 2690 bar [39,29]. To determine the numerical endurance limit fatigue crack growth analyses for different operating pressures are performed. The results of these analyses are the evolution of D Keff over the crack length. If DKeff is equal or below the threshold value, which is known from previpffiffiffiffiffi ous investigations [47] as DKeff,th = 3 MPa m, crack arrest is predicted and the structure is considered to be fatigue resistant.
In Fig. 7 the evolution of the number of nodes and elements during the crack growth analysis for an operating pressure of Dp = 2950 bar is depicted. The corresponding finite element meshes for crack lengths of about 0.01 mm, 0.16 mm and 0.66 mm at 90% of the operating pressure are shown in Fig. 8. By looking at these figures one big advantage of the mapping analysis, which uses the remeshing technique, becomes obvious. It is possible to model the crack growth from very small crack lengths up to large crack lengths with a moderate number of nodes and elements.
Fig. 5. Crack opening pressure for the nodal-release and the mapping analysis.
Fig. 6. Evolution of DKeff over crack length.
78
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79
Fig. 7. Evolution of the number of nodes and elements during the crack growth simulation.
Fig. 9. Calculated and experimental [39,29] fatigue lives.
numerical results is expected. Another factor, which has to be also considered is that the crack growth constants in this computations have been identified in the medium air. In this example the medium oil is existent. This might lead to different constants due to the fact of fluid-induced crack closure [16] and with that to different fatigue lives. 4. Conclusions The numerical crack growth simulation has become more and more attractive today for reasons of time and economy. The presentation of a procedure for such a simulation of fatigue crack growth in elastic–plastic materials is the aim of this paper. The remeshing of the structure and the subsequent mapping of the status variables from the old mesh to the new one after each step of crack advance are thereby the core characteristics of the proposed procedure. For the validation of the procedure results on autofrettaged Diesel-engine injection tubes are used. The results obtained by this procedure using an advanced material model, which is able to realistically simulate cyclic plasticity effects such as mean stress relaxation and ratcheting, are compared with experimental results and finite element results based on a nodal-release scheme. From these results the following conclusions can be drawn:
Fig. 8. Initial mesh (a0 = 0.01 mm) and meshed geometries for crack lengths of a 0.16 mm and a 0.66 mm (paf = 4700 bar, Dp = 2950 bar).
In Fig. 9 the calculated fatigue lives and the experimental ones are shown. As expected, the fatigue lives calculated via Eq. (26) are by a factor of about two longer than the ones calculated via Eq. (25). But they are by a factor of about three to four shorter than the experimental results. The reasons for this underestimation of the fatigue lives are not only due to the mapping of the variables, but also on other factors, which will be discussed in the following. The most important point is that because of the axis symmetric modeling of the structure a straight through the thickness crack is implied. But it is known from the experiments that semi elliptical surface cracks initiated. These cracks have a stress intensity factor, which is (at the deepest position) much lower than the one for a through the thickness crack of the same length. By taking into account these three-dimensional effects, a significant improvement of the
The results obtained by the proposed procedure are in good agreement with the results obtained by the nodal-release scheme. The results obtained by the proposed procedure are in good agreement with the experimental results concerning the endurance limit. The fatigue lives are underestimated by the proposed procedure. This is not only due to the inaccuracies resulting from the mapping of the variables but also due to the simplified axis symmetric modeling and the choice of the crack growth constants. A great advantage of this procedure is that arbitrary loading conditions, i.e. proportional as well as non-proportional loading, can be examined. But further validations of the procedure are required. Also an extension of the implementation of the procedure to account for three-dimensional structures is under progress. References [1] D.V. Swenson, A.R. Ingraffea, Computational Mechanics 3 (1988) 381–397. [2] T.N. Bittencourt, P.A. Wawrzynek, A.R. Ingraffea, J.L. Sousa, Engineering Fracture Mechanics 55 (1996) 321–334.
P. Zerres, M. Vormwald / Computational Materials Science 57 (2012) 73–79 [3] L.F. Martha, Topological and Geometrical Modeling Approach to Numerical Discretization and Arbitrary Fracture Simulation in Three Dimensions, Ph.D. Thesis, School of Civil Engineering, Cornell University, 1989. [4] B.J. Carter, P.A. Wawrzynek, A.R. Ingraffea, Gallagher Special Issue of International Journal for Numerical Methods in Engineering 47 (2000) 229– 253. [5] Zentech International Ltd., Zencrack. Tool for 3D Fracture Mechanics Simulation. [6] C. Timbrell, G. Cook, 3D FE fracture mechanics analysis for industrial applications, in: Seminar on Inelastic Finite Element Analysis, Institute of Mechanical Engineers, London, 1997. [7] M. Fulland, Risssimulation in dreidimensionalen Strukturen mit automatischer adaptiver Finite-Elemente-Netzgenerierung, Fortschrittbericht-Bericht VDI, Reihe 18: Mechanik/Bruchmechanik Nr. 280, VDI-Verlag, Dnsseldorf, 2003. [8] M. Schöllmann, M. Fulland, H.A. Richard, Engineering Fracture Mechanics 70 (2003) 249–268. [9] S. Mellings, J. Baynham, R.A. Adey, Engineering Analysis with Boundary Elements 26 (2002) 479–488. [10] F. Erdogan, G.C. Sih, Journal of Basic Engineering 85 (1963) 519–527. [11] M. Schöllmann, G. Kullmer, M. Fulland, H.A. Richard, A new criterion for 3D crack growth under mixed-mode (I + II + III) loading, in: M. de Freitas (Ed.), Proceedings of the 6th International Conference on Biaxial/Multiaxial Fatigue and Fracture, 2001, pp. 589–596. [12] M. Schöllmann, H.A. Richard, G. Kullmer, M. Fulland, International Journal of Fracture 117 (2002) 129–141. [13] N. Moes, J. Dolbow, T. Belytschko, International Journal for Numerical Methods in Engineering 46 (1999) 131–150. [14] T. Belytschko, Y.Y. Lu, L. Gu, International Journal for Numerical Methods in Engineering 37 (1994) 229–256. [15] J.C. Newman Jr., A crack-closure model for predicting fatigue crack growth under aircraft spectrum loading, in: J.B. Chang, C.M. Hudson (Eds.), Methods and Models for Predicting Fatigue Crack Growth Under Random Loading, ASTM STP 748, 1981, pp. 53–84. [16] W. Elber, Engineering Fracture Mechanics 2 (1970) 37–45. [17] R.C. McClung, H. Sehitoglu, Engineering Fracture Mechanics 33 (1989) 237– 252. [18] R.C. McClung, H. Sehitoglu, Engineering Fracture Mechanics 33 (1989) 253– 272. [19] J.D. Dougherty, J. Padovan, T.S. Srivatsan, Engineering Fracture Mechanics 56 (1997) 189–212. [20] A. Gonzales-Herrera, J. Zapatero, Engineering Fracture Mechanics 72 (2005) 337–355. [21] Y. Jiang, M. Feng, F. Ding, International Journal of Plasticity 21 (2005) 1720– 1740. [22] N.E. Ashbaugh, B. Dattaguru, M. Khobaib, T. Nicholas, R.V. Prakash, T.S. Ramamurthy, B.R. Seshadri, R. Sunder, Fatigue & Fracture of Engineering Materials & Structures 20 (1997) 963–974. [23] S. Pommier, P. Bompard, Fatigue & Fracture of Engineering Materials & Structures 23 (2000) 129–139. [24] D.M. Rodrigues, F.V. Antunes, Engineering Fracture Mechanics 76 (2009) 1215–1230. [25] K. Solanki, S.R. Daniewicz, J.C. Newman Jr., Engineering Fracture Mechanics 70 (2003) 1475–1489. [26] K. Solanki, S.R. Daniewicz, J.C. Newman Jr., Engineering Fracture Mechanics 71 (2004) 1165–1175. [27] S. Roychowdhury, R.H. Dodds Jr., Fatigue & Fracture of Engineering Materials & Structures 26 (2003) 663–673.
79
[28] S. Roychowdhury, R.H. Dodds Jr., Engineering Fracture Mechanics 70 (2003) 2363–2383. [29] E. Herz, R. Thumser, J.W. Bergmann, M. Vormwald, Engineering Fracture Mechanics 73 (2006) 3–21. [30] F.V. Antunes, D.M. Rodrigues, Engineering Fracture Mechanics 75 (2008) 3101–3120. [31] J.C. Newman, Jr., A finite-element analysis of fatigue crack closure, in: J.R. Rice, P.C. Paris (Eds.), Mechanics of Crack Growth, ASTM STP 590, 1976, pp. 281– 301. [32] S. Pommier, International Journal of Fatigue 23 (2001) 111–118. [33] K. Ogura, K. Ohji, K. Honda, Influence of mechanical factors on the fatigue crack closure, in: D.M.R. Taplin (Ed.), Fracture, 1977 (vol. 2), ICF4, Waterloo, Canada, 1977, pp. 1035–1047. [34] R.C. McClung, B.H. Thacker, S. Roy, International Journal of Fracture 50 (1991) 27–49. [35] K. Solanki, S.R. Daniewicz, J.C. Newman Jr., Engineering Fracture Mechanics 71 (2004) 149–171. [36] E. Herz, O. Hertel, M. Vormwald, Engineering Fracture Mechanics 78 (2011) 559–572. [37] Dassault Systémes Simulia Corp., Abaqus Analysis User’s Manual (6.8ef), 2008. [38] P. Zerres, M. Vormwald, Modeling of crack propagation under mixed-mode loading with advanced remeshing, in: Proceedings of the 12th International Conference on Fracture, Ottawa, Kanada, 2009. [39] E. Herz, R. Thumser, M. Vormwald, J. Bergmann, Endurance limit design of high pressure Diesel injection pipes for low failure probability, in: Proceedings of Fatigue Design, Senlis, 2005, pp. 1–8. [40] R. Döring, J. Hoffmeyer, T. Seeger, Computational Materials Science 28 (2003) 587–596. [41] R. Döring, Zum Deformations- und Schädigungsverhalten metallischer Werkstoffe unter mehrachsig nichtproportionalen zyklischen Beanspruchungen, Ph.D. Thesis, Institut für Stahlbau und Werkstoffmechanik, TU Darmstadt, Heft 78, 2006. [42] P. Zerres, Numerische Simulation des Ermüdungsrissfortschrittes in metallischen Strukturen unter Berücksichtigung zyklischer Plastizitätseffekte, Ph.D. Thesis, Institut für Stahlbau und Werkstoffmechanik, TU Darmstadt, Heft 88, 2010. [43] J.L. Chaboche, K. Dang Van, G. Cordier, Modelization of the strain memory effect on the cyclic hardening of 316 stainless steel, in: Transactions of the 5th International Conference on Structural Mechanics in Reactor Technology, Div. L, 1979, p. L11/3. [44] E. Tanaka, European Journal of Mechanics, A/Solids 13 (1994) 155–173. [45] Y. Jiang, P. Kurath, Journal of Engineering Materials and Technology 119 (1997) 161–170. [46] C.-C. Ma, J.-I. Huang, C.-H. Tsai, Journal of Pressure Vessels Technology 116 (1994) 423–430. [47] Windte, K., Experimentelle Bestimmung des Schwellwertes fnr Risswachstum unter zyklischer Beanspruchung zur Verbesserung der rechnerischen Lebensdauervorhersage bei Betriebsbelastung, 1986. [48] M. Vormwald, Anrisslebensdauervorhersage auf der Basis der Schwingbruchmechanik für kurze Risse, Ph.D. Thesis, Institut für Stahlbau und Werkstoffmechanik, TH Darmstadt, Heft 47, 1989. [49] R.J. Donahue, H. McI Clark, P. Atanmo, R. Kumble, A.J. McEvily, Crack opening displacement and the rate of fatigue crack growth, International Journal of Fracture 8 (1972) 209–219. [50] P. Zerres, R. Döring, M. Vormwald, Materials Science and Engineering Technology 42 (2011) 1093–1098.