- Email: [email protected]
A modified zigzag approach to approximate moving crack front with arbitrary shape
Engineering Fracture Mechanics 78 (2011) 234–251
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
A modified zigzag approach to approximate moving crack front with arbitrary shape Yan-Ping Liu, Chuan-Yao Chen, Guo-Qing Li ⇑ School of Civil Engineering and Mechanics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, PR China
a r t i c l e
i n f o
Article history: Received 30 March 2010 Received in revised form 28 July 2010 Accepted 8 August 2010 Available online 27 August 2010 Keywords: Virtual crack closure technique Strain energy release rate Finite element method Zigzag mesh pattern Fatigue crack growth
a b s t r a c t This paper proposed a modified zigzag approach to approximate crack front with arbitrary shape based on virtual crack closure-integral technique. The important features of this approach are that: (i) more detailed consideration to determine the required virtually closed area and displacement opening were presented and (ii) adaptive re-meshing technique was avoid in studying crack growth. Three cases were presented to assess the validity of this approach by comparison to the results obtained from other approaches. Reasonable agreement between the present study and the analytical solutions were obtained. Also, the shape changes during crack propagation can be tracked with ease. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction In linear elastic fracture mechanics (LEFM), the strain energy release rate (G) is an important parameter that has been used to establish various fracture criteria. However, it is very difficult to determine the closed-form G solutions for a crack with an arbitrary front in a three-dimensional finite body, and most often it is an intractable task. Therefore, the numerical methods such as finite element method (FEM) have been extensively developed. In FEM, The virtual crack closure technique (VCCT) had been widely used to evaluate the strain energy release rate (G) because of its particular convenience in FEM analysis. The virtual crack closure technique was originated from Irwin’s crack closure integral [1]. Rybicki and Kanninen [2] modified Irwin’s crack closure integral and proposed VCCT for two dimensional finite element calculations. And then, Shivakumar et al. [3] extended this VCCT to three dimensions calculations. Review on the topic can be found in [4,5]. In general, with VCCT, the expression for G is pretty straightforward and the FEM analysis for G is very simple. The required nodal forces, virtually closed area and the displacement openings to evaluate G can be easily obtained from a single finite element analysis without any special treatment in the region around crack tip. With VCCT, although no requirement on a theoretical basis, calculations are usually performed with orthogonal mesh patterns where the elements are oriented normal and parallel to the crack front. When the crack front is determined or the crack growth is self-similar, the orthogonal finite element meshes can be relatively easily achieved and retained. However, for complex crack shapes, such as the intersection of crack with free surface [6,7], the coalesced multi-cracks [8,9] or the kinking cracks [10,11], this kind of fine matching between meshes and geometry is time consuming and sometimes difficult to achieve. In order to overcome these difficulties, the VCCT calculation using non-orthogonal mesh pattern had attracted attentions from many researchers. Fawaz [6,7] assessed the validity of 3-D VCCT with non-orthogonal finite element meshes for ⇑ Corresponding author. E-mail address: [email protected] (G.-Q. Li). 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.08.007
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
235
Nomenclature Ai,1, Ai,2, Ai,3 and Ai,4 area around crack tip node Ni a, c, C1, and C2 crack length b the width of element faces at the crack front E and v Young’s modulus and Poisson’s rate crack growth parameter. Ed Fx, Fy and Fz the shear forces and opening force for the crack tip node GI, GII and GIII modes I–III strain energy release rate GIC, GIIC and GIIIC the critical values of G corresponding to modes I–III loading KI, KII and KIII modes I–III stress intensity factor m parameter to describe the location of nodes ni,x, ni,y the coordinates of normal vector for node Ni in global coordinate system Ni,x, Ni,y the coordinates of node Ni in global coordinate system R, H radius and height of solid cylinder Rbi,x, Rbi,y, Rai,x and Rai,y the coordinates of shape vector for node Ni in global coordinate system S, Si,1 and Si,2 the totally virtually closed area, the virtually closed area ahead of the crack front and the virtually closed area behind the crack front for node Ni t, L, h plate thickness, plate width and plate height ux, uy and uz the shear and displacement opening from the node behind the crack front ^ the base vector of global coordinate system ^i; ^j and k ~ ni ;~ t i and ~ ki the local coordinates system for node Ni ~ ni;a ; ~ ni;b normal to the shape vector for node Ni Da1, Da2 length of finite element across the crack front Du2, Dv2 and Dw the displacement openings in local system (~ ni ; ~ ti ; ~ ki Þ a the ratio of element lengths across the crack front x, b and c the exponents in mix-mode fracture criterion n0, g0 and f0 the coordinates of point N00i in the local coordinate system (n, g, f) rz tensile stress in ‘Z’ direction h parametric angle of an ellipse and a circle
calculating mode I stress intensity factors. For the straight-through crack model, the FEM calculations of SIF using orthogonal and non-orthogonal meshes were analyzed. Comparison on the results indicated that the difference in SIF results is only 0.14% for the linear skew mesh and 0.008% for the chevron skew mesh when compared to the results for the orthogonal mesh. Recently, some corrections on VCCT were made to investigate the influence of non-orthogonal mesh patterns on G calculation. Okada et al. [12,13] presented some developments on VCCT for three dimensional fracture analyses. Using their formulations, the non-equal sized and non-orthogonally arranged mesh pattern across the crack front can be used to calculate G with excellent correlation to analytical solutions. Smith and Raju [14] proposed new formulations of VCCT and other methods to evaluate stress intensity factors with general finite element models that lack orthogonality with the crack front. Average unit normal to the crack front was used to determine the required virtually closed area which was constructed from the slope of line segments adjacent to crack tip node. This work has features similar in many ways to the present paper. In the above work, the non-orthogonal finite element meshes were used to calculate G. The crack front was approximated in a smooth curve and crack tip nodes are arranged to precisely lie along the actual crack front. On the other hand, Ferrie and Rousseau [15] proposed a ‘‘stepped” pattern to approximate the cracked region instead of a smooth curve for a circular crack. A modified virtual crack closure technique was developed to evaluate the G at the corner crack nodes. The directions of crack propagation of each corner node were assumed to parallel to the line bisecting the element attached to the corner node and behind the crack tip, and a factor kf was proposed to determine the required virtual closed area. The FEM results indicated that the G at corner nodes had an acceptable degree of accuracy. Xie and Biggers [16,17] suggested to use an approximated zigzag line to represent the complex crack front which has the same general shape for the given crack front. The VCCT procedure based on this approximate zigzag crack front was established and the corresponding ABAQUS user subroutine was developed in their publications. Of course, this zigzag approximation does not require the use of orthogonal meshes. Even more interesting, the approach avoids re-meshing technique during the FEM analysis as crack front changes and therefore, it is particular attractive in studying crack growth with a relatively simple and effective way. Using this approximated zigzag crack front, Xie et al. [18] performed the crack growth analyses of surface cracks based discrete cohesive zone mode proposed by Xie and Waas [19]. Boscolo et al. [20] followed the same idea as Xie and Biggers to study their own problem of a debonding joint and achieved a very close shape of the debonded front compared to their test coupons. Since the nodes are not lie exactly along the actual crack front, the G results for the approximate zigzag crack front oscillate around the analytical solutions with non-negligible amplitudes. In this paper, an effort was made to investigate
236
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
a modified zigzag approach based on the approach proposed by Xie and Biggers [16,17]. In the present study, the normal vectors for each node on the crack front were determined firstly. Based on this normal vector, the approach to determine virtually closed area, nodal force and displacement opening was proposed. For comparison purpose, the differences of approaches between present study and Xie’ work [16] were discussed in details. Apart from Xie and Bigger’s work [16,17], more detailed consideration on virtually closed area calculation was given in present study. Also, the influence of different finite element size and different faces of element across the crack front on G calculations were considered. Finally, three cases were studied to illustrate the validity and accuracy of present approach. Two cases were firstly examined to assess the accuracy of G calculation using present approach. The G results are compared with these obtained by Xie and Biggers [16,17] and analytical solutions. The third case was examined to assess the ease of use in studying crack growth. The crack shape variations with crack propagation were discussed in details. 2. Methodology and formulation To evaluate strain energy release rate numerically, an eight-node element (hexahedron element) is usually used as shown in Fig. 1. The strain energy release rate (G) for the node i with an orthogonal mesh pattern can be written as [4,5]
1 GI ¼ 2S F z ðujz ujz Þ
1 F x ðujx ujx Þ GII ¼ 2S
ð1Þ
1 F y ðujy ujy Þ GIII ¼ 2S
where Fx and Fy are the shear forces and Fz is the opening force from the node i on the crack front; ux, uy and uz are the relative shear and relative displacement openings obtained from the node pair j–j behind the crack front; S is the virtually closed area where S = bDa as shown in Fig. 1. 2.1. Determination of the approximate zigzag crack front For a fixed crack with arbitrary shape, the cracked surface and uncracked surface can be detected by projecting the crack front onto the crack plane, as shown in Fig. 2. The locations for the finite element nodes embedded in the crack plane are defined by a variable m, as
m ¼ 1 for the nodes located on the uncracked surface m ¼ 0 for the nodes located on the cracked surface
u zj
j u zj*
u yj
u
j* z
u
Fz
u xj j* y
u
j*
z kˆ
ð2Þ
Fy
Fx
i
S
k
j* x
y ˆj x iˆ
a
Δa
Δa
Fig. 1. Three dimensional VCCT with eight nodes isoparametric elements.
b
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
237
Uncracked surface
mi , j +1 = 1
Cracked surface mi , j = 1 mi +1, j = 0
mi −1, j = 1
y ˆj mi , j −1 = 1
Crack tip nodes
x iˆ
: Crack tip nodes
Fig. 2. Determination of the approximate zigzag crack front.
n i = n i , a + n i ,b
N i +1 ni ,a
Ra i
Ni
ni ,b
y
ˆj
Rbi ti
x iˆ
N i −1
Fig. 3. Calculation of normal vector for the crack tip node Ni.
The crack tip nodes were determined by the variable m, as
mi;j ¼ 1 mi1;j miþ1;j ¼ 0 or mi;j1 mi;jþ1 ¼ 0
ð3Þ
where the subscripts ‘‘i” and ‘‘j” denote the locations of nodes in crack plane as shown in Fig. 2. All the crack tip nodes are physical finite element nodes. In order to ensure the geometrical continuity in the finite element analysis, the distances between the adjacent crack tip nodes in x and y directions should not exceed the element sizes. In Fig. 2, the lines which link up the adjacent crack tip nodes were considered as the approximate zigzag crack front. In this paper, the basic assumption for the G calculation is that the crack tip nodes will be propagated along the normal ! ! vector. In order to calculate ~ ni for node Ni, two shape vectors Rbi and Rai were proposed in Fig. 3, as !
Rbi ¼ Rbi;x^i þ Rbi;y^j ¼ ðNi1;x Ni;x Þ^i þ ðNi1;y Ni;y Þ^j
ð4Þ
!
Rai ¼ Rai;x^i þ Rai;y^j ¼ ðNiþ1;x Ni;x Þ^i þ ðNiþ1;y Ni;y Þ^j
ð5Þ
where Rbi,x, Rbi,y, Rai,x, Rai,y, Ni1,x, Ni1,y, Ni,x, Ni,y, Ni+1,x and N iþ1;y are the coordinates in global coordinate system, ^i and ^j are the unit vectors corresponding to the global coordinate (x, y) direction. Following the right-hand-rule for the cross product, the vector ~ ni which is normal to the approximate zigzag crack front at node Ni can be obtained by !
^ Ra ¼ k ^ ðRa ^i þ Ra ^jÞ ¼ Ra ^i þ Ra ^j ~ ni;a ¼ k i i;x i;y i;y i;x
ð6Þ
!
^ Rb ¼ k ^ ðRb ^i þ Rb ^jÞ ¼ Rb ^i Rb ^j ~ ni;b ¼ k i;x i;y i;y i;x i ^ ~ ni ¼ ~ n þ~ ni;a ¼ ni;x i þ ni;y^j ¼ ðRbi;y Rai;y Þ^i þ ðRai;x Rbi;x Þ^j i;b
ð7Þ ð8Þ
238
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
^ is the unit vector normal to ^i and ^j. where k The corresponding tangent vector for ~ ni are
^ ~ ~ ti ¼ k ni ¼ ni;y^i þ ni;x^j ~ ^ k ¼k
ð9Þ ð10Þ
i
Therefore, the local coordinate system for node Ni can be obtained from these three vectors ~ ni ;~ ti and ~ ki . And then, components of strain energy release rate (GI, GII and GIII) can be calculated within this local coordinate system. 2.2. Calculation of G for the moving approximate zigzag crack front To evaluate G, the nodal forces come from the nodes on the crack front, the virtually closed area ahead of the crack front and the displacement opening behind the crack front are required. Since the crack tip nodes were determined in Section 2.1, the nodal force can be obtained directly from the FEM analysis results. The virtually closed area and displacement opening will be calculated in sequence in this section. A component composed of crack tip nodes Ni1, Ni and Ni+1 was illustrated to calculate the virtually closed area for node Ni, as shown in Fig. 4. The virtual points on the normal to the crack front ahead of and behind the crack tip node Ni are shown as N 0i and N 00i in Fig. 4. Area around crack tip node Ni can be divided into four parts: Ai,1, Ai,2, Ai,3 and Ai,4. The corresponding four points to enclose these areas are referring to Fig. 4. For the areas Ai,1 and Ai,3, crack tip nodes Ni and Ni1 have been already known from the finite element meshes, and the other four points (N 0i ; N 0i1 ; N 00i and N 00i1 ) are to be determined. For the node Ni1, the normal vector is coincided with the edge of element. Thus, the coordinates of N 0i1 and N 00i1 are equal to the coordinates of the finite element nodes on the normal to
(a)
N i'+1
N i +1
Ai ,2 N
' i
Ai ,4
N i''+1
Ni Ai ,1 y
ˆj
N i'−1
Ai ,3
N i''
Δa2
N i''−1
N i −1 Δa1
x iˆ
a
(b)
N i'+1
N i +1
N
' i
Si ,1
Ni
N i''+1
Si ,2
N i'' y ˆj
N i'−1
N i −1
a
N i''−1
x iˆ Fig. 4. (a) Schematic of area around of nodes Ni in present work and (b) virtual area ahead of and behind the crack front.
239
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
the crack front, as shown in Fig. 4. The locations of the remaining two points N 0i and N 00i are calculated as follows. In Fig. 4, the point N 0i is the intersection of: (1) the line parallel to the vector ~ ni passing through crack tip node Ni and (2) the line parallel to ! the vector Rbi passing through point N 0i1 . Therefore, the two linear equations can be proposed as
(
ðN0i;x N i;x Þni;y ðN0i;y Ni;y Þni;x ¼ 0
ð11Þ
ðN0i;x N 0i1;x ÞRbi;y ðN0i;y N0i1;y ÞRbi;x ¼ 0 Solving the above two equations yields the coordinates of point N 0i , that is
8 ðN 0 N ÞRb n ðN0 N i;y ÞRbi;x ni;x > < N0i;x ¼ Ni;x þ i1;x i;x Rbi;y ni;x Rbi1;y i;y i;x i;x ni;y
ð12Þ
0 0 > : N0 ¼ Ni;y þ ðNi1;x Ni;x ÞRbi;y ni;y ðNi1;y Ni;y ÞRbi;x ni;y i;y Rb n Rb n i;y i;x
i;x i;y
Following the same calculating procedure, the coordinates of point N 00i can be obtained as
8 ðN 00 N ÞRb n ðN00 N i;y ÞRbi;x ni;x > < N00i;x ¼ Ni;x þ i1;x i;x Rbi;y ni;x Rbi1;y i;y i;x i;x ni;y
ð13Þ
00 00 > : N00 ¼ Ni;y þ ðNi1;x Ni;x ÞRbi;y ni;y ðNi1;y Ni;y ÞRbi;x ni;y i;y Rbi;y ni;x Rbi;x ni;y
Once the coordinates of Ni, Ni1, N 0i ; N 0i1 ; N 00i and N 00i1 are determined, the areas Ai,1 and Ai,3 can be calculated. The same procedure can be also used to calculate areas Ai,2 and Ai,4. With the areas Ai,1, Ai,2, Ai,3 and Ai,4, the totally virtual closed areas S were given by Okada [12], which is
S¼
pffiffiffi
a Si;1
1 Si;2 Si;1 2ð1 þ aÞ a
ð14Þ
where a is the ratio of element widths across the crack front, a = Da2/Da1 in Fig. 4. Si,1 is the area of face of the finite element ahead of the crack front. Si,2 is the area of face of the finite element behind the crack front. For the approximate zigzag crack front shown in Fig. 4, Si,1 and Si,2 for node Ni are obtained by
: m=0 : m=1
Fracture occured Step 2:
Step 1:
Determine the new crack tip, cracked surface and uncracked surface according to parameter m.
Calculate G for the current crack front, and determine whether fracture occurs?
Fig. 5. Crack front determination for each step.
Uncracked surface
Ni N i −1
Uncracked surface
N i +1
Ni+2
N i N i +1 N i −1
Cracked surface
Cracked surface y ˆj x iˆ
N i −2
N i −2
N i −3
N i −3
(a)
(b)
Fig. 6. Two approximate zigzag crack fronts with different shapes.
Ni+2
240
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
N i'+1
(a)
N
Si ,+11
Si ,1
Ni
N
ˆj
Ni+2
N i"+1
N i"+2
N i''
N i −1
N
y
N i +1 ' i −1
Si −1,1
N i'−2
N i'+2
' i
Si −2,1
N i'−3
" i −1
N i −2
N i"−2
N i −3
N i"−3
x iˆ
(b)
N i' Si +1,1
Si ,1
Ni N i +1
N i'−1 Si −1,1
Ni+2
N i''
N i −1
N i''+1
N i''+2
N i'−2
N i"−1 Si −2,1
N i −2
y
N i"−2
ˆj N i'−3
N i −3
N i"−3
x iˆ Fig. 7. Schematic of virtually closed crack areas using Xie’s [16] work, (a) for the crack front (a) in Fig. 6; (b) for the crack front (b) in Fig. 6.
Si;1 ¼ cb Ai;1 þ ca Ai;2 and Si;2 ¼ cb Ai;3 þ ca Ai;4
ð15Þ
where ca = 1 when Ni1 does not correspond to a physical node or 0.5 when Ni1 dose correspond to a node. The same can be said for cb but the governing node is Ni+1.The values of crack displacement opening at points N 00i could be calculated by interpolating the values from the nodes of the finite element model. This interpolation could be accomplished using the element shape functions. The shape functions Lj for the eight nodes isoparametric elements are
241
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
N i'
(a)
N i'+1 Si +1,1
Si ,1
N i +1
Ni
N i'−1
N i''+1
N i''
Si −1,2
Ni+2
Si +1,2
Si ,2
Si −1,1 N i −1
N i'+2
N i''+2
N i"−1 N i'−2
y
ˆj
Si −2,1
N i −2
Si −2,2
N i''−2 N i"−3
N i −3
N i'−3
Δa2
Δa1
x iˆ
N i'
(b)
N i'+1
N i'+2
Si +1,1
Si ,1
Ni+2
N i'−1
N i +1
Ni
Si +1,2
Si ,2
Si −1,1
N i''
N i −1
N i'−2
N i''+1
N i''+2
Si −1,2 Si −2,1
N i −2
N i"−1 Si −2,2
N i''−2 y
ˆj N i'−3
N i −3
N i"−3
x iˆ Fig. 8. Schematic of virtually closed crack areas using present work, (a) for the crack front (a) in Fig. 6; (b) for the crack front (b) in Fig. 6.
Lj ðn; g; 1Þ ¼
1 ð1 þ nj nÞð1 þ gj gÞð1 þ fj fÞ ðj ¼ 1; 2; . . . ; 8Þ 8
ð16Þ
where the subscript ‘‘j” denotes the node for the element, (n, g, f) is the local coordinate system located at the centre of element.
242
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
The crack displacement opening at point N 00i in the global system are
Du ¼
8 X
Dv ¼
Lj ðn0 ; g0 ; 10 ÞDuj ;
j¼1
8 X
Lj ðn0 ; g0 ; 10 ÞDv j ;
Dw ¼
j¼1
8 X
Lj ðn0 ; g0 ; 10 ÞDwj
ð17Þ
j¼1
where (n0, g0, f0) are the locations of point N 00i in the local coordinate system, (Duj, Dvj, Dwj) are the displacement opening of the nodes in the global system. The Dw is required to calculate GI and can be obtained directly using Eq. (17). The Du and Dv should be transformed into the local system ð~ ni ; ~ ti ; ~ ki Þ to calculate GII and GIII as
ni;y ni;x Du2 ¼ Dv qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Du qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðni;y Þ þ ðni;x Þ ðni;y Þ2 þ ðni;x Þ2 ni;x ni;y Dv 2 ¼ Dv qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ Du qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðni;y Þ þ ðni;x Þ ðni;y Þ2 þ ðni;x Þ2
ð18Þ
Since the nodal force on the crack front, the displacement opening behind the crack front and the virtually closed area are determined, the components of strain energy release rate can be calculated by Eq. (1). When the components of strain energy release rates are calculated, a mix-mode fracture criterion is then used to predict the crack growth.
Ed ¼
GI GIC
x
þ
GII GIIC
b
þ
GIII GIIIC
c
P1
ð19Þ
where Ed is the crack growth parameter. GIC, GIIC and GIIIC are the critical values corresponding to modes I–III loading, respectively. In this work, the exponents x, b and c were taken as 1 as suggested by Refs. [16,21]. Once the fracture criterion in Eq. (19) is satisfied (Ed P 1), variable m was shifted from 1 to 0. The new cracked surface, uncracked surface, crack tip nodes and the approximate zigzag crack front can be detected, as shown in Fig. 5. The above procedure for calculation of energy release rate and the application of the failure criterion were written in ANSYS Parametric Design Language (APDL) and were carried out through the ANSYS user subroutine.
2.3. Comparisons of G calculations between present study and Xie’s work In this section, the differences to calculate G between Xie’s approach and proposed approach are discussed in aspects: (1) virtually closed area correction; (2) displacement opening correction. Fig. 6 shows two approximate zigzag crack fronts with different shape. In Fig. 6, although the locations of nodes Ni, Ni1 and Ni+1 are same for each other, the nodes Ni2 and Ni3 have different locations. The corresponding virtual areas around the crack tip nodes are plotted in Fig. 7 using Xie and Bigger’s approach. For comparison purpose, the virtually areas determined by present approach are plotted in Fig. 8. In Fig. 7a and b, it can be found that the virtually closed areas around node Ni are same for each other although the approximate zigzag crack fronts have different shape. Meanwhile, the virtually areas are different for each other as shown in Fig. 8a and b. The reason cased to this result is that, the virtually area calculation in present approach are not only depending on the Ni components but also related to the crack front shape near the components.
N i +1 Uncracked surface
Ni Cracked surface
N i −1 y
ˆj
N i −2 x iˆ
Δa1
Δa2
Fig. 9. The approximate zigzag crack front with different element length.
243
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
N i'+1
(a)
N i' N i'−1
y
Si ,1
Si −1,1
N i'−2
ˆj x iˆ
Δa1
N i'+1
N i +1
N i"+1
Ni
N i''
N i'
N i −1
N i"−1
N i'−1
N i −2
N i"+2
Δa1
Δa2
(b)
y
Si ,1
Si −1,1
N i'−2
ˆj x iˆ
N i"+1
N i +1 Ni
N i −1
Si ,2
N i''
Si −1,2
N i"−1
N i −2 Δa1
N i"+2
Δa2
Fig. 10. Schematic of virtually closed area for the crack front in Fig. 9. (a) Xie and Biggers’ approach; (b) present study.
Another important feature can be also found in Figs. 7 and 8. In Fig. 7a and b, there is a blank space between Si,1 and Si1,1. The virtually area around the crack front for neighbor nodes (Ni1 and Ni) are not fully covered using Xie and Biggers’ approach. With the proposed approach, this problem is solved, where the blank space is not presented as shown in Fig. 8a and b. As to the displacement opening correction, the approaches to determine the displacement opening are different between Xie’s work and the present study. Fig. 9 shows a simple zigzag approximate crack front, in which the finite elements across the crack front have different length. The displacement openings determined by Xie and Biggers’ approach and by present study are illustrated in Fig. 10. It is observed that, in Xie’s work, the point N 00i is located on a line the direction of which is determined by the normal vector. This point N 0i has the same distance but in the opposite direction of point N 0i , as shown in Fig. 10a. In present work, the point N 00i is coincided with the finite element nodes, as shown in Fig. 10b. Because the point N 0i is located very near the crack front, to ensure the accuracy of interpolated values of displacement at this point, the displacement evaluated at the nodes of finite element model are preferable [4,14]. 3. Case studies In this section, three cases were examined to illustrate the validity and accuracy of the proposed VCCT process. The present approach was firstly applied to evaluate G of Circular crack and semi-elliptical surface crack. Three kinds of mesh patterns including orthogonal mesh pattern, non-orthogonal mesh pattern and approximate zigzag crack front were used. The third case was a through crack with an elliptical crack front, which is embedded in a plate. The crack propagation was studied using present approach in the third case. For all the analyses presented, we used eight nodes isoparametric element. Young’s modulus and Poisson’s ratio are set to be 200 GPa and 0.3, respectively. The analyses ensure the correct implementation of the 3-D VCCT into the user subroutine of ANSYS. In order to compare with the analytical solutions or the other established FEM solutions, the stress intensity factor (SIF) was used for conveniences. For the isotropic material, the mode I stress intensity factor, KI, can be calculated by the following relations:
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi GI E KI ¼ for plane strain 1 v2 pffiffiffiffiffiffiffiffi K I ¼ GI E for plane stress where E and
ð20Þ
v are the material Young’s modulus and Poisson’s ratio, respectively.
3.1. A circular crack embedded in a solid cylinder subject to tensile loading The first example is the problem of a circular crack embedded in a solid cylinder subjected to tensile loading, as depicted in Fig. 11. Due to symmetry of geometry and loading condition, one quarter of the solid cylinder was modeled as shown in Fig. 12a. The finite elements with same length are symmetrically placed across the crack front in Fig. 12b. The length of elements across the crack front is set to be 0.05a. The total numbers of elements and nodes are 17,500 and 19,526, respectively. According to the geometry of the crack, the cracked surface, uncracked surface and the crack tip nodes can be automatically identified. The required points ahead of and behind the crack tip nodes can be determined by using the presented approach. It is noticed that, when the meshes pattern in Fig. 12b are used, the points N 00i and N 0i do correspond to the physical finite element node, as shown in Fig. 12b.
244
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
σz
a
R=10a
H=2R
Z Y
σz
X
Fig. 11. A circular crack embedded in a solid cylinder.
The stress intensity factors (SIF) normalized with respect to the value K0 of analytical solution are plotted in Fig. 13. The K results qualitatively agree with the analytical solutions [22] as expected. The difference between present work and analytical solutions is just 1.2%. It can be also observed that, the K results evaluated by present study have small deviation from the results obtained by Xie’s approach [16]. Their difference is just 0.06%. The reason for this result is that, virtual areas around of the crack tip node (Si,1 and Si,2) and the displacement openings are almost same for each other. The mesh pattern in Fig. 12b is then applied to investigate the influence of the different element length across the crack front on G calculations. The detailed finite elements across the crack front with different length are illustrated in Fig. 14. The ratios a of the element length are set to be 1:5, 3:5, 5:5, 5:3 and 5:1 in Fig. 14.
(a)
(b) Crack front
UnCracked surface
N i"
Ni N i' Cracked surface
0.05a
0.05a
Fig. 12. The finite element model used to calculate G for the circular crack: (a) finite element model; (b) detailed mesh around the crack front.
245
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Fig. 13. Comparison of SIF between present study and the analytical solution.
The normalized SIF under the different ratios a calculated by presented study are illustrated in Fig. 15a. It is found that the computed SIF increased with the increasing element length ratios in Fig. 15a. When the element length ahead of the crack front is smaller, the computed stress intensity factors show slightly higher and vice versa. These trends agree well with the three dimensional solutions given in Okada [12] for through crack problem. For comparison purpose, the K results calculated by Xie and Biggers’ approach [16] are plotted in Fig. 15b. It is seen that the computed stress intensity factors decreased with the increasing element ratios in Fig. 15b. These trends are contrary to the observations in Fig. 15a and in Okada et al. [12]. The reason caused to these results maybe boiled down to the displacement opening calculation as discussed in Section 2.3. Nevertheless, the present modified approach provided more accuracy K results when dealing with the finite element meshes with different element length.
3.2. Semi-elliptical surface crack subjected to tensile loading The second example is the problem of semi-elliptical surface crack, as shown in Fig. 16. Due to the symmetry of geometry and loading condition, a quarter of the plate was modeled as well. Two mesh patterns were used for this case as shown in Fig. 17a and b. In Fig. 17a, the crack tip nodes lie precisely on the actual crack front and the mesh is non-orthogonal. The aim is to evaluate the ability of the presented VCCT process with the non-orthogonal mesh pattern. The element length across the crack front is set to be 0.05a. The total numbers of elements and nodes are 18,900 and 21,168 for the model in Fig. 15a,
UnCracked surface
Crack front N i"
Ni N
' i
Cracked surface
Δ2
Analysis cases No. 1 2 3 4 5
α = Δ1: Δ 2
1:5 3:5 5:5 5:3 5:1
Δ1
Fig. 14. Detailed mesh around the crack front with different element length.
246
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Fig. 15. Comparison of SIF between analytical solutions and FEM results with different element length: (a) using the present approach; (b) using Xie and Biggers’ approach.
σz
t
a h
2c Crack
A
A
L Section A-A c/L=0.25 a/t=1/3 h/L=1.5 a/c=0.5
Z Y X
σz Fig. 16. A semi-elliptical surface crack subject to tensile loading.
247
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
(a)
y UnCracked surface
Cracked surface
0.05a
a
x
c
0.05a
(b) Approximate crack front
y
UnCracked surface
a Cracked surface
x c Fig. 17. Detailed mesh around the surface crack front. (a) Exact elliptical crack front, non-orthogonal mesh; (b) approximate zigzag crack front, nonorthogonal mesh.
In Fig. 17b, the semi-elliptical surface crack was approximated to the zigzag line. The nodes on either side of the actual crack front are defined as crack tip nodes. The zigzag line in Fig. 17b approximates the smooth crack front curve as closely as possible. In order to diminish errors of geometrical shape simulations, the element length and width are set to be 0.025a. The total numbers of elements and nodes are 35,000 and 38,556, for the mode in Fig. 15b. Fig. 18 shows the normalized SIF for the surface crack with the mesh pattern in Fig. 17a, which are obtained by three different approaches. The results given by Newman and Raju SIF equation [23] are also plotted in Fig. 18. Because the crack tip
Fig. 18. Comparison of analytical solution and present study with the mesh (a) in Fig. 17.
248
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Fig. 19. Comparison of analytical solution, present study and Xie’s work, with the mesh (b) in Fig. 17.
nodes lie exactly on the actual crack front, the K results in Fig. 16 reveal a good correlation for this case with non-orthogonal mesh pattern as expected. The maximum deviation between present study and Newman and Raju’ solution is 2.9% with the non-orthogonal mesh pattern. It can be also observed that, the K results evaluated by present study have small deviation from the results obtained by Xie’s approach [16,17] and by Smith’s approach [14]. Their maximum difference for Xie’ approach and Smith’ approach are 0.50% and 0.31%, respectively. Fig. 19 shows the normalized SIF for this case with the mesh pattern in Fig. 15b. The K results calculated by Xie and Biggers [16,17] with the mesh pattern in Fig. 15b and by Newman and Raju SIF equation [23] are also plotted in Fig. 19. At the extremes of 2//p = 0 and 2//p = 1, the normalized K results obtained by different approaches have a very small deviation for each other. Only 1.4% and 1.1% errors between present study and Newman and Raju’ solutions are found for these two nodes, respectively. For the other nodes, the K results calculated by present study and Xie and Biggers [16,17] oscillate around the analytical solution as expected. The reason caused to this result is because of the locations of the nodes on the approximate zigzag crack front. The nodes are lie either outside or inside the curved front. However, by using present approach, the amplitude of the oscillations is obviously smaller than that for Xie’s approach as shown in Fig. 19. 3.3. Through crack with an elliptical crack front The application of the present approach to analyze moving crack front is discussed in this section, where the dynamic effects should be taken into consideration [24,25]. A through crack with an elliptical crack front was examined, in which
Fig. 20. Through crack with an elliptical crack front.
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Crack front
Top surface
249
C2
t
bottom surface
C1
Fig. 21. Detailed mesh around a through crack with an elliptical crack front.
the crack shapes will change instantaneously during crack propagation. The geometry of the plate and the initial crack shapes are shown in Fig. 20. The through crack with initial shape was approximated to the zigzag line, as shown in Fig. 21. The tensile stress rZ = 150 MPa is applied along the Z-direction on plate. This is a typical mode I dominated crack. Thus, the fracture criterion can be simplified as
Ed ¼
GI GIC
x
1
ð21Þ
where x = 1 and GIC is assumed to be 7.144 N/mm in this work. Fig. 22 shows the normalized SIF for the through crack with the initial elliptical crack front, in which the zigzag mesh pattern shown in Fig. 21 are used. The K results calculated by Xie and Biggers [16,17] with the same mesh pattern are also plotted in Fig. 21. The K results calculated by present study and Xie [16,17] oscillate around the Fawaz’ solution [6] as expected. The amplitude of the oscillations for present study is smaller than that for Xie’s approach. It can be also found that, the normalized K has the maximum value for the nodes at the top surface. The K results are much lower for the nodes at the bottom surface. The crack shapes variations with crack propagation are illustrated in Fig. 23 using the proposed approach. The crack shape varied in two stages. At the first stage, only the nodes near top surface of the plate were break, as shown in Fig. 23a–e. The reason can be boiled down to the G distributions along the crack front, which had been illustrated in Fig. 22. With the crack propagation, the distance between the nodes at top surface and at bottom surface decreased. As a result, the elliptical crack front changes into a straight line, as shown in Fig. 23a–e. At the second stage, the nodes at the mid-thickness were firstly break, as shown in Fig. 23f–h. Because, the strain energy release rate G has the largest value at mid-thickness and drops at the free surface of the plate [12] with a through-thickness crack. The crack shape seems to a convex line at the mid-thickness, as shown in Fig. 23h. The results in Fig. 23 indicated that, the present approach can enable crack growth to be tracked as the shape changes. Since the K results oscillate around the analytical solution, the accuracy of present approach, the sensitivity study of modeling details and the software convergence are presented in authors’ following paper [26].
Fig. 22. The normalized SIF for the through crack with initial elliptical crack front using zigzag mesh pattern.
250
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
(c)
(b)
(a)
(e)
(d)
Step: 14
Step: 10
Step: 5
(f)
Step: 25
Step: 20
Step: 30
(h)
(g)
Step: 40
Step: 50 Fig. 23. Crack shapes under different steps.
4. Conclusions In this paper, a modified zigzag approach to approximate moving crack front with arbitrary shape was proposed based on the VCCT. Some useful conclusions are summarized as follows. 1. A modified zigzag approach to evaluate G was presented for arbitrary crack front using non-orthogonal mesh pattern. There are two major modifications on the Xie and Biggers’ approach in present study. Firstly, the present approach provided a more detailed consideration to determine the virtually closed area. There is no blank space between the virtually closed areas for the consecutive crack tip nodes. Secondly, the present approach can be used to calculate G when the finite elements across the crack front have different sizes. 2. Three cases were examined to study the validity and accuracy of the proposed approach. The G results calculated in present study show excellent correlation with the analytical solutions or the other established FEM solutions. For the zigzag pattern, the G results calculated by present study and Xie oscillate around the analytical solution. However, by using present study, the amplitude of the oscillations is obviously smaller than that for Xie’s approach. 3. In order to study the crack growth behavior, a through crack with an elliptical crack front was examined using the zigzag mesh pattern. Although the crack shape changes during crack propagation, a simple stationary finite element mesh was used without adaptive re-meshing. And the results indicated that the changes of crack shape can be tracked by using the present approach.
References [1] Irwin GR. Analysis of stresses and strain near the end of a crack traversing a plate. J Appl Mech 1957;24:361–4. [2] Rybicki EF, Kanninen MF. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng Fract Mech 1977;9:931–8. [3] Shivakumar KN, Tan PW, Newman Jr JC. A virtual crack-closure technique for calculating stress intensity factors for cracked three dimensional bodies. Int J Fract 1988;36:R43–50. [4] Krueger R. Virtual crack closure technique: history, approach, and applications. Appl Mech Rev 2004;57(2):109–43. [5] Xie D, Biggers SB. Progressive crack growth analysis using interface element based on the virtual crack closure technique. Finite Elem Anal Des 2006;42:977–84.
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
251
[6] Fawaz SA. Application of the virtual crack closure technique to calculate stress intensity factors for through cracks with an elliptical crack front. Eng Fract Mech 1998;59(3):327–42. [7] Fawaz SA. Stress intensity factor solutions for part-elliptical through cracks. Eng Fract Mech 1999;63:209–26. [8] Lin XB, Smith RA. A numerical simulation of fatigue growth of multiple surface initially semicircular defects under tension. Int J Pres Ves Piping 1995;62:281–9. [9] Lin XB, Smith RA. Fatigue growth analysis of interacting and coalescing surface crack. Int J Fract 1997;85:283–99. [10] Xie D, Wass AM, Shahwan KW, Schroeder JA, Beoman RG. Fracture criterion for kinking cracks in a tri-material adhesively bonded joint under mixed mode loading. Eng Fract Mech 2005;72:2487–504. [11] Xie D, Wass AM, Shahwan KW, Schroeder JA, Beoman RG. Computation of energy release rates for kinking cracks based on virtually crack closure technique. CNES-Comput Model Engng Sci 2004;6:515–24. [12] Okada H, Higashi M, et al. Three dimensional virtual crack closure-integral method (VCCM) with skewed and non-symmetric mesh arrangement at the crack front. Engng Fract Mech 2005;72:1717–37. [13] Okada H, Endoh S, Kikuchi M. On fracture analysis using an element overlay technique. Engng Fract Mech 2005;72:773–89. [14] Smith SA, Raju IS. Evaluation of stress-intensity factors using general finite-element models. Fatigue Fract Mech. In: Panontin TL, Sheppard SD, editors. ASTM STP 1332, vol. 29. West Canshohockem (PA): American Society for testing and materials; 1999, p. 176–200. [15] Ferrie CH, Rousseau CQ. A method of applying VCCT to corner crack nodes. In: Proc of Am Soc for composites, 16th annual technical conference on composite materials, Blacksburg, VA, USA; 2001. [16] Xie D, Biggers SB. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique. Part I: formulation and validation. Engng Fract Mech 2006;73:771–85. [17] Xie D, Biggers SB. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique. Part II: sensitivity study on modeling details. Engng Fract Mech 2006;73:786–801. [18] Xie D, Gary M, Huang DD, Abdi F. Cohesive zone model for surface cracks using finite element analysis. AIAA-2008-106742, Chicago, Illinois; April 2008. [19] Xie D, Waas AM. Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Engng Fract Mech 2006;73:1783–96. [20] Boscolo M, Allegri G, Zhang X. Design and modeling of selective reinforcements for integral aircraft structures. AIAA J 2009;46(9):2323–31. [21] Kutlu Z, Chang FK. Modeling compressive failure of laminated composite containing multiple through-the-width delaminations. J Compos Mater 1992;26:350–87. [22] Murakami Y et al. Stress intensity factors handbook. Berlin: Pergamon Press; 1987. [23] Newman JC, Raju IS. An empirical stress-intensity factor equation for the surface crack. Engng Fract Mech 1981;15(1–2):185–92. [24] Chen XD, He XM. The effect of the recess shape on performance analysis of the gas-lubricated bearing in optical lithography. Tribol Int 2006;39(11):1336–41. [25] He XM, Chen XD. The dynamic analysis of the gas lubricated stage in optical lithography. Int J Adv Manufact Tech 2007;32(9–10):978–84. [26] Liu YP. Study on fatigue crack growth behavior in welded steel plates used for bridge. Ph.D thesis. China: Huazhong University of Science and Technology; Aug 2010.
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
A modified zigzag approach to approximate moving crack front with arbitrary shape Yan-Ping Liu, Chuan-Yao Chen, Guo-Qing Li ⇑ School of Civil Engineering and Mechanics, Huazhong University of Science and Technology, Wuhan 430074, Hubei, PR China
a r t i c l e
i n f o
Article history: Received 30 March 2010 Received in revised form 28 July 2010 Accepted 8 August 2010 Available online 27 August 2010 Keywords: Virtual crack closure technique Strain energy release rate Finite element method Zigzag mesh pattern Fatigue crack growth
a b s t r a c t This paper proposed a modified zigzag approach to approximate crack front with arbitrary shape based on virtual crack closure-integral technique. The important features of this approach are that: (i) more detailed consideration to determine the required virtually closed area and displacement opening were presented and (ii) adaptive re-meshing technique was avoid in studying crack growth. Three cases were presented to assess the validity of this approach by comparison to the results obtained from other approaches. Reasonable agreement between the present study and the analytical solutions were obtained. Also, the shape changes during crack propagation can be tracked with ease. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction In linear elastic fracture mechanics (LEFM), the strain energy release rate (G) is an important parameter that has been used to establish various fracture criteria. However, it is very difficult to determine the closed-form G solutions for a crack with an arbitrary front in a three-dimensional finite body, and most often it is an intractable task. Therefore, the numerical methods such as finite element method (FEM) have been extensively developed. In FEM, The virtual crack closure technique (VCCT) had been widely used to evaluate the strain energy release rate (G) because of its particular convenience in FEM analysis. The virtual crack closure technique was originated from Irwin’s crack closure integral [1]. Rybicki and Kanninen [2] modified Irwin’s crack closure integral and proposed VCCT for two dimensional finite element calculations. And then, Shivakumar et al. [3] extended this VCCT to three dimensions calculations. Review on the topic can be found in [4,5]. In general, with VCCT, the expression for G is pretty straightforward and the FEM analysis for G is very simple. The required nodal forces, virtually closed area and the displacement openings to evaluate G can be easily obtained from a single finite element analysis without any special treatment in the region around crack tip. With VCCT, although no requirement on a theoretical basis, calculations are usually performed with orthogonal mesh patterns where the elements are oriented normal and parallel to the crack front. When the crack front is determined or the crack growth is self-similar, the orthogonal finite element meshes can be relatively easily achieved and retained. However, for complex crack shapes, such as the intersection of crack with free surface [6,7], the coalesced multi-cracks [8,9] or the kinking cracks [10,11], this kind of fine matching between meshes and geometry is time consuming and sometimes difficult to achieve. In order to overcome these difficulties, the VCCT calculation using non-orthogonal mesh pattern had attracted attentions from many researchers. Fawaz [6,7] assessed the validity of 3-D VCCT with non-orthogonal finite element meshes for ⇑ Corresponding author. E-mail address: [email protected] (G.-Q. Li). 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.08.007
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
235
Nomenclature Ai,1, Ai,2, Ai,3 and Ai,4 area around crack tip node Ni a, c, C1, and C2 crack length b the width of element faces at the crack front E and v Young’s modulus and Poisson’s rate crack growth parameter. Ed Fx, Fy and Fz the shear forces and opening force for the crack tip node GI, GII and GIII modes I–III strain energy release rate GIC, GIIC and GIIIC the critical values of G corresponding to modes I–III loading KI, KII and KIII modes I–III stress intensity factor m parameter to describe the location of nodes ni,x, ni,y the coordinates of normal vector for node Ni in global coordinate system Ni,x, Ni,y the coordinates of node Ni in global coordinate system R, H radius and height of solid cylinder Rbi,x, Rbi,y, Rai,x and Rai,y the coordinates of shape vector for node Ni in global coordinate system S, Si,1 and Si,2 the totally virtually closed area, the virtually closed area ahead of the crack front and the virtually closed area behind the crack front for node Ni t, L, h plate thickness, plate width and plate height ux, uy and uz the shear and displacement opening from the node behind the crack front ^ the base vector of global coordinate system ^i; ^j and k ~ ni ;~ t i and ~ ki the local coordinates system for node Ni ~ ni;a ; ~ ni;b normal to the shape vector for node Ni Da1, Da2 length of finite element across the crack front Du2, Dv2 and Dw the displacement openings in local system (~ ni ; ~ ti ; ~ ki Þ a the ratio of element lengths across the crack front x, b and c the exponents in mix-mode fracture criterion n0, g0 and f0 the coordinates of point N00i in the local coordinate system (n, g, f) rz tensile stress in ‘Z’ direction h parametric angle of an ellipse and a circle
calculating mode I stress intensity factors. For the straight-through crack model, the FEM calculations of SIF using orthogonal and non-orthogonal meshes were analyzed. Comparison on the results indicated that the difference in SIF results is only 0.14% for the linear skew mesh and 0.008% for the chevron skew mesh when compared to the results for the orthogonal mesh. Recently, some corrections on VCCT were made to investigate the influence of non-orthogonal mesh patterns on G calculation. Okada et al. [12,13] presented some developments on VCCT for three dimensional fracture analyses. Using their formulations, the non-equal sized and non-orthogonally arranged mesh pattern across the crack front can be used to calculate G with excellent correlation to analytical solutions. Smith and Raju [14] proposed new formulations of VCCT and other methods to evaluate stress intensity factors with general finite element models that lack orthogonality with the crack front. Average unit normal to the crack front was used to determine the required virtually closed area which was constructed from the slope of line segments adjacent to crack tip node. This work has features similar in many ways to the present paper. In the above work, the non-orthogonal finite element meshes were used to calculate G. The crack front was approximated in a smooth curve and crack tip nodes are arranged to precisely lie along the actual crack front. On the other hand, Ferrie and Rousseau [15] proposed a ‘‘stepped” pattern to approximate the cracked region instead of a smooth curve for a circular crack. A modified virtual crack closure technique was developed to evaluate the G at the corner crack nodes. The directions of crack propagation of each corner node were assumed to parallel to the line bisecting the element attached to the corner node and behind the crack tip, and a factor kf was proposed to determine the required virtual closed area. The FEM results indicated that the G at corner nodes had an acceptable degree of accuracy. Xie and Biggers [16,17] suggested to use an approximated zigzag line to represent the complex crack front which has the same general shape for the given crack front. The VCCT procedure based on this approximate zigzag crack front was established and the corresponding ABAQUS user subroutine was developed in their publications. Of course, this zigzag approximation does not require the use of orthogonal meshes. Even more interesting, the approach avoids re-meshing technique during the FEM analysis as crack front changes and therefore, it is particular attractive in studying crack growth with a relatively simple and effective way. Using this approximated zigzag crack front, Xie et al. [18] performed the crack growth analyses of surface cracks based discrete cohesive zone mode proposed by Xie and Waas [19]. Boscolo et al. [20] followed the same idea as Xie and Biggers to study their own problem of a debonding joint and achieved a very close shape of the debonded front compared to their test coupons. Since the nodes are not lie exactly along the actual crack front, the G results for the approximate zigzag crack front oscillate around the analytical solutions with non-negligible amplitudes. In this paper, an effort was made to investigate
236
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
a modified zigzag approach based on the approach proposed by Xie and Biggers [16,17]. In the present study, the normal vectors for each node on the crack front were determined firstly. Based on this normal vector, the approach to determine virtually closed area, nodal force and displacement opening was proposed. For comparison purpose, the differences of approaches between present study and Xie’ work [16] were discussed in details. Apart from Xie and Bigger’s work [16,17], more detailed consideration on virtually closed area calculation was given in present study. Also, the influence of different finite element size and different faces of element across the crack front on G calculations were considered. Finally, three cases were studied to illustrate the validity and accuracy of present approach. Two cases were firstly examined to assess the accuracy of G calculation using present approach. The G results are compared with these obtained by Xie and Biggers [16,17] and analytical solutions. The third case was examined to assess the ease of use in studying crack growth. The crack shape variations with crack propagation were discussed in details. 2. Methodology and formulation To evaluate strain energy release rate numerically, an eight-node element (hexahedron element) is usually used as shown in Fig. 1. The strain energy release rate (G) for the node i with an orthogonal mesh pattern can be written as [4,5]
1 GI ¼ 2S F z ðujz ujz Þ
1 F x ðujx ujx Þ GII ¼ 2S
ð1Þ
1 F y ðujy ujy Þ GIII ¼ 2S
where Fx and Fy are the shear forces and Fz is the opening force from the node i on the crack front; ux, uy and uz are the relative shear and relative displacement openings obtained from the node pair j–j behind the crack front; S is the virtually closed area where S = bDa as shown in Fig. 1. 2.1. Determination of the approximate zigzag crack front For a fixed crack with arbitrary shape, the cracked surface and uncracked surface can be detected by projecting the crack front onto the crack plane, as shown in Fig. 2. The locations for the finite element nodes embedded in the crack plane are defined by a variable m, as
m ¼ 1 for the nodes located on the uncracked surface m ¼ 0 for the nodes located on the cracked surface
u zj
j u zj*
u yj
u
j* z
u
Fz
u xj j* y
u
j*
z kˆ
ð2Þ
Fy
Fx
i
S
k
j* x
y ˆj x iˆ
a
Δa
Δa
Fig. 1. Three dimensional VCCT with eight nodes isoparametric elements.
b
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
237
Uncracked surface
mi , j +1 = 1
Cracked surface mi , j = 1 mi +1, j = 0
mi −1, j = 1
y ˆj mi , j −1 = 1
Crack tip nodes
x iˆ
: Crack tip nodes
Fig. 2. Determination of the approximate zigzag crack front.
n i = n i , a + n i ,b
N i +1 ni ,a
Ra i
Ni
ni ,b
y
ˆj
Rbi ti
x iˆ
N i −1
Fig. 3. Calculation of normal vector for the crack tip node Ni.
The crack tip nodes were determined by the variable m, as
mi;j ¼ 1 mi1;j miþ1;j ¼ 0 or mi;j1 mi;jþ1 ¼ 0
ð3Þ
where the subscripts ‘‘i” and ‘‘j” denote the locations of nodes in crack plane as shown in Fig. 2. All the crack tip nodes are physical finite element nodes. In order to ensure the geometrical continuity in the finite element analysis, the distances between the adjacent crack tip nodes in x and y directions should not exceed the element sizes. In Fig. 2, the lines which link up the adjacent crack tip nodes were considered as the approximate zigzag crack front. In this paper, the basic assumption for the G calculation is that the crack tip nodes will be propagated along the normal ! ! vector. In order to calculate ~ ni for node Ni, two shape vectors Rbi and Rai were proposed in Fig. 3, as !
Rbi ¼ Rbi;x^i þ Rbi;y^j ¼ ðNi1;x Ni;x Þ^i þ ðNi1;y Ni;y Þ^j
ð4Þ
!
Rai ¼ Rai;x^i þ Rai;y^j ¼ ðNiþ1;x Ni;x Þ^i þ ðNiþ1;y Ni;y Þ^j
ð5Þ
where Rbi,x, Rbi,y, Rai,x, Rai,y, Ni1,x, Ni1,y, Ni,x, Ni,y, Ni+1,x and N iþ1;y are the coordinates in global coordinate system, ^i and ^j are the unit vectors corresponding to the global coordinate (x, y) direction. Following the right-hand-rule for the cross product, the vector ~ ni which is normal to the approximate zigzag crack front at node Ni can be obtained by !
^ Ra ¼ k ^ ðRa ^i þ Ra ^jÞ ¼ Ra ^i þ Ra ^j ~ ni;a ¼ k i i;x i;y i;y i;x
ð6Þ
!
^ Rb ¼ k ^ ðRb ^i þ Rb ^jÞ ¼ Rb ^i Rb ^j ~ ni;b ¼ k i;x i;y i;y i;x i ^ ~ ni ¼ ~ n þ~ ni;a ¼ ni;x i þ ni;y^j ¼ ðRbi;y Rai;y Þ^i þ ðRai;x Rbi;x Þ^j i;b
ð7Þ ð8Þ
238
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
^ is the unit vector normal to ^i and ^j. where k The corresponding tangent vector for ~ ni are
^ ~ ~ ti ¼ k ni ¼ ni;y^i þ ni;x^j ~ ^ k ¼k
ð9Þ ð10Þ
i
Therefore, the local coordinate system for node Ni can be obtained from these three vectors ~ ni ;~ ti and ~ ki . And then, components of strain energy release rate (GI, GII and GIII) can be calculated within this local coordinate system. 2.2. Calculation of G for the moving approximate zigzag crack front To evaluate G, the nodal forces come from the nodes on the crack front, the virtually closed area ahead of the crack front and the displacement opening behind the crack front are required. Since the crack tip nodes were determined in Section 2.1, the nodal force can be obtained directly from the FEM analysis results. The virtually closed area and displacement opening will be calculated in sequence in this section. A component composed of crack tip nodes Ni1, Ni and Ni+1 was illustrated to calculate the virtually closed area for node Ni, as shown in Fig. 4. The virtual points on the normal to the crack front ahead of and behind the crack tip node Ni are shown as N 0i and N 00i in Fig. 4. Area around crack tip node Ni can be divided into four parts: Ai,1, Ai,2, Ai,3 and Ai,4. The corresponding four points to enclose these areas are referring to Fig. 4. For the areas Ai,1 and Ai,3, crack tip nodes Ni and Ni1 have been already known from the finite element meshes, and the other four points (N 0i ; N 0i1 ; N 00i and N 00i1 ) are to be determined. For the node Ni1, the normal vector is coincided with the edge of element. Thus, the coordinates of N 0i1 and N 00i1 are equal to the coordinates of the finite element nodes on the normal to
(a)
N i'+1
N i +1
Ai ,2 N
' i
Ai ,4
N i''+1
Ni Ai ,1 y
ˆj
N i'−1
Ai ,3
N i''
Δa2
N i''−1
N i −1 Δa1
x iˆ
a
(b)
N i'+1
N i +1
N
' i
Si ,1
Ni
N i''+1
Si ,2
N i'' y ˆj
N i'−1
N i −1
a
N i''−1
x iˆ Fig. 4. (a) Schematic of area around of nodes Ni in present work and (b) virtual area ahead of and behind the crack front.
239
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
the crack front, as shown in Fig. 4. The locations of the remaining two points N 0i and N 00i are calculated as follows. In Fig. 4, the point N 0i is the intersection of: (1) the line parallel to the vector ~ ni passing through crack tip node Ni and (2) the line parallel to ! the vector Rbi passing through point N 0i1 . Therefore, the two linear equations can be proposed as
(
ðN0i;x N i;x Þni;y ðN0i;y Ni;y Þni;x ¼ 0
ð11Þ
ðN0i;x N 0i1;x ÞRbi;y ðN0i;y N0i1;y ÞRbi;x ¼ 0 Solving the above two equations yields the coordinates of point N 0i , that is
8 ðN 0 N ÞRb n ðN0 N i;y ÞRbi;x ni;x > < N0i;x ¼ Ni;x þ i1;x i;x Rbi;y ni;x Rbi1;y i;y i;x i;x ni;y
ð12Þ
0 0 > : N0 ¼ Ni;y þ ðNi1;x Ni;x ÞRbi;y ni;y ðNi1;y Ni;y ÞRbi;x ni;y i;y Rb n Rb n i;y i;x
i;x i;y
Following the same calculating procedure, the coordinates of point N 00i can be obtained as
8 ðN 00 N ÞRb n ðN00 N i;y ÞRbi;x ni;x > < N00i;x ¼ Ni;x þ i1;x i;x Rbi;y ni;x Rbi1;y i;y i;x i;x ni;y
ð13Þ
00 00 > : N00 ¼ Ni;y þ ðNi1;x Ni;x ÞRbi;y ni;y ðNi1;y Ni;y ÞRbi;x ni;y i;y Rbi;y ni;x Rbi;x ni;y
Once the coordinates of Ni, Ni1, N 0i ; N 0i1 ; N 00i and N 00i1 are determined, the areas Ai,1 and Ai,3 can be calculated. The same procedure can be also used to calculate areas Ai,2 and Ai,4. With the areas Ai,1, Ai,2, Ai,3 and Ai,4, the totally virtual closed areas S were given by Okada [12], which is
S¼
pffiffiffi
a Si;1
1 Si;2 Si;1 2ð1 þ aÞ a
ð14Þ
where a is the ratio of element widths across the crack front, a = Da2/Da1 in Fig. 4. Si,1 is the area of face of the finite element ahead of the crack front. Si,2 is the area of face of the finite element behind the crack front. For the approximate zigzag crack front shown in Fig. 4, Si,1 and Si,2 for node Ni are obtained by
: m=0 : m=1
Fracture occured Step 2:
Step 1:
Determine the new crack tip, cracked surface and uncracked surface according to parameter m.
Calculate G for the current crack front, and determine whether fracture occurs?
Fig. 5. Crack front determination for each step.
Uncracked surface
Ni N i −1
Uncracked surface
N i +1
Ni+2
N i N i +1 N i −1
Cracked surface
Cracked surface y ˆj x iˆ
N i −2
N i −2
N i −3
N i −3
(a)
(b)
Fig. 6. Two approximate zigzag crack fronts with different shapes.
Ni+2
240
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
N i'+1
(a)
N
Si ,+11
Si ,1
Ni
N
ˆj
Ni+2
N i"+1
N i"+2
N i''
N i −1
N
y
N i +1 ' i −1
Si −1,1
N i'−2
N i'+2
' i
Si −2,1
N i'−3
" i −1
N i −2
N i"−2
N i −3
N i"−3
x iˆ
(b)
N i' Si +1,1
Si ,1
Ni N i +1
N i'−1 Si −1,1
Ni+2
N i''
N i −1
N i''+1
N i''+2
N i'−2
N i"−1 Si −2,1
N i −2
y
N i"−2
ˆj N i'−3
N i −3
N i"−3
x iˆ Fig. 7. Schematic of virtually closed crack areas using Xie’s [16] work, (a) for the crack front (a) in Fig. 6; (b) for the crack front (b) in Fig. 6.
Si;1 ¼ cb Ai;1 þ ca Ai;2 and Si;2 ¼ cb Ai;3 þ ca Ai;4
ð15Þ
where ca = 1 when Ni1 does not correspond to a physical node or 0.5 when Ni1 dose correspond to a node. The same can be said for cb but the governing node is Ni+1.The values of crack displacement opening at points N 00i could be calculated by interpolating the values from the nodes of the finite element model. This interpolation could be accomplished using the element shape functions. The shape functions Lj for the eight nodes isoparametric elements are
241
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
N i'
(a)
N i'+1 Si +1,1
Si ,1
N i +1
Ni
N i'−1
N i''+1
N i''
Si −1,2
Ni+2
Si +1,2
Si ,2
Si −1,1 N i −1
N i'+2
N i''+2
N i"−1 N i'−2
y
ˆj
Si −2,1
N i −2
Si −2,2
N i''−2 N i"−3
N i −3
N i'−3
Δa2
Δa1
x iˆ
N i'
(b)
N i'+1
N i'+2
Si +1,1
Si ,1
Ni+2
N i'−1
N i +1
Ni
Si +1,2
Si ,2
Si −1,1
N i''
N i −1
N i'−2
N i''+1
N i''+2
Si −1,2 Si −2,1
N i −2
N i"−1 Si −2,2
N i''−2 y
ˆj N i'−3
N i −3
N i"−3
x iˆ Fig. 8. Schematic of virtually closed crack areas using present work, (a) for the crack front (a) in Fig. 6; (b) for the crack front (b) in Fig. 6.
Lj ðn; g; 1Þ ¼
1 ð1 þ nj nÞð1 þ gj gÞð1 þ fj fÞ ðj ¼ 1; 2; . . . ; 8Þ 8
ð16Þ
where the subscript ‘‘j” denotes the node for the element, (n, g, f) is the local coordinate system located at the centre of element.
242
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
The crack displacement opening at point N 00i in the global system are
Du ¼
8 X
Dv ¼
Lj ðn0 ; g0 ; 10 ÞDuj ;
j¼1
8 X
Lj ðn0 ; g0 ; 10 ÞDv j ;
Dw ¼
j¼1
8 X
Lj ðn0 ; g0 ; 10 ÞDwj
ð17Þ
j¼1
where (n0, g0, f0) are the locations of point N 00i in the local coordinate system, (Duj, Dvj, Dwj) are the displacement opening of the nodes in the global system. The Dw is required to calculate GI and can be obtained directly using Eq. (17). The Du and Dv should be transformed into the local system ð~ ni ; ~ ti ; ~ ki Þ to calculate GII and GIII as
ni;y ni;x Du2 ¼ Dv qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Du qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðni;y Þ þ ðni;x Þ ðni;y Þ2 þ ðni;x Þ2 ni;x ni;y Dv 2 ¼ Dv qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ Du qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðni;y Þ þ ðni;x Þ ðni;y Þ2 þ ðni;x Þ2
ð18Þ
Since the nodal force on the crack front, the displacement opening behind the crack front and the virtually closed area are determined, the components of strain energy release rate can be calculated by Eq. (1). When the components of strain energy release rates are calculated, a mix-mode fracture criterion is then used to predict the crack growth.
Ed ¼
GI GIC
x
þ
GII GIIC
b
þ
GIII GIIIC
c
P1
ð19Þ
where Ed is the crack growth parameter. GIC, GIIC and GIIIC are the critical values corresponding to modes I–III loading, respectively. In this work, the exponents x, b and c were taken as 1 as suggested by Refs. [16,21]. Once the fracture criterion in Eq. (19) is satisfied (Ed P 1), variable m was shifted from 1 to 0. The new cracked surface, uncracked surface, crack tip nodes and the approximate zigzag crack front can be detected, as shown in Fig. 5. The above procedure for calculation of energy release rate and the application of the failure criterion were written in ANSYS Parametric Design Language (APDL) and were carried out through the ANSYS user subroutine.
2.3. Comparisons of G calculations between present study and Xie’s work In this section, the differences to calculate G between Xie’s approach and proposed approach are discussed in aspects: (1) virtually closed area correction; (2) displacement opening correction. Fig. 6 shows two approximate zigzag crack fronts with different shape. In Fig. 6, although the locations of nodes Ni, Ni1 and Ni+1 are same for each other, the nodes Ni2 and Ni3 have different locations. The corresponding virtual areas around the crack tip nodes are plotted in Fig. 7 using Xie and Bigger’s approach. For comparison purpose, the virtually areas determined by present approach are plotted in Fig. 8. In Fig. 7a and b, it can be found that the virtually closed areas around node Ni are same for each other although the approximate zigzag crack fronts have different shape. Meanwhile, the virtually areas are different for each other as shown in Fig. 8a and b. The reason cased to this result is that, the virtually area calculation in present approach are not only depending on the Ni components but also related to the crack front shape near the components.
N i +1 Uncracked surface
Ni Cracked surface
N i −1 y
ˆj
N i −2 x iˆ
Δa1
Δa2
Fig. 9. The approximate zigzag crack front with different element length.
243
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
N i'+1
(a)
N i' N i'−1
y
Si ,1
Si −1,1
N i'−2
ˆj x iˆ
Δa1
N i'+1
N i +1
N i"+1
Ni
N i''
N i'
N i −1
N i"−1
N i'−1
N i −2
N i"+2
Δa1
Δa2
(b)
y
Si ,1
Si −1,1
N i'−2
ˆj x iˆ
N i"+1
N i +1 Ni
N i −1
Si ,2
N i''
Si −1,2
N i"−1
N i −2 Δa1
N i"+2
Δa2
Fig. 10. Schematic of virtually closed area for the crack front in Fig. 9. (a) Xie and Biggers’ approach; (b) present study.
Another important feature can be also found in Figs. 7 and 8. In Fig. 7a and b, there is a blank space between Si,1 and Si1,1. The virtually area around the crack front for neighbor nodes (Ni1 and Ni) are not fully covered using Xie and Biggers’ approach. With the proposed approach, this problem is solved, where the blank space is not presented as shown in Fig. 8a and b. As to the displacement opening correction, the approaches to determine the displacement opening are different between Xie’s work and the present study. Fig. 9 shows a simple zigzag approximate crack front, in which the finite elements across the crack front have different length. The displacement openings determined by Xie and Biggers’ approach and by present study are illustrated in Fig. 10. It is observed that, in Xie’s work, the point N 00i is located on a line the direction of which is determined by the normal vector. This point N 0i has the same distance but in the opposite direction of point N 0i , as shown in Fig. 10a. In present work, the point N 00i is coincided with the finite element nodes, as shown in Fig. 10b. Because the point N 0i is located very near the crack front, to ensure the accuracy of interpolated values of displacement at this point, the displacement evaluated at the nodes of finite element model are preferable [4,14]. 3. Case studies In this section, three cases were examined to illustrate the validity and accuracy of the proposed VCCT process. The present approach was firstly applied to evaluate G of Circular crack and semi-elliptical surface crack. Three kinds of mesh patterns including orthogonal mesh pattern, non-orthogonal mesh pattern and approximate zigzag crack front were used. The third case was a through crack with an elliptical crack front, which is embedded in a plate. The crack propagation was studied using present approach in the third case. For all the analyses presented, we used eight nodes isoparametric element. Young’s modulus and Poisson’s ratio are set to be 200 GPa and 0.3, respectively. The analyses ensure the correct implementation of the 3-D VCCT into the user subroutine of ANSYS. In order to compare with the analytical solutions or the other established FEM solutions, the stress intensity factor (SIF) was used for conveniences. For the isotropic material, the mode I stress intensity factor, KI, can be calculated by the following relations:
rffiffiffiffiffiffiffiffiffiffiffiffiffiffi GI E KI ¼ for plane strain 1 v2 pffiffiffiffiffiffiffiffi K I ¼ GI E for plane stress where E and
ð20Þ
v are the material Young’s modulus and Poisson’s ratio, respectively.
3.1. A circular crack embedded in a solid cylinder subject to tensile loading The first example is the problem of a circular crack embedded in a solid cylinder subjected to tensile loading, as depicted in Fig. 11. Due to symmetry of geometry and loading condition, one quarter of the solid cylinder was modeled as shown in Fig. 12a. The finite elements with same length are symmetrically placed across the crack front in Fig. 12b. The length of elements across the crack front is set to be 0.05a. The total numbers of elements and nodes are 17,500 and 19,526, respectively. According to the geometry of the crack, the cracked surface, uncracked surface and the crack tip nodes can be automatically identified. The required points ahead of and behind the crack tip nodes can be determined by using the presented approach. It is noticed that, when the meshes pattern in Fig. 12b are used, the points N 00i and N 0i do correspond to the physical finite element node, as shown in Fig. 12b.
244
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
σz
a
R=10a
H=2R
Z Y
σz
X
Fig. 11. A circular crack embedded in a solid cylinder.
The stress intensity factors (SIF) normalized with respect to the value K0 of analytical solution are plotted in Fig. 13. The K results qualitatively agree with the analytical solutions [22] as expected. The difference between present work and analytical solutions is just 1.2%. It can be also observed that, the K results evaluated by present study have small deviation from the results obtained by Xie’s approach [16]. Their difference is just 0.06%. The reason for this result is that, virtual areas around of the crack tip node (Si,1 and Si,2) and the displacement openings are almost same for each other. The mesh pattern in Fig. 12b is then applied to investigate the influence of the different element length across the crack front on G calculations. The detailed finite elements across the crack front with different length are illustrated in Fig. 14. The ratios a of the element length are set to be 1:5, 3:5, 5:5, 5:3 and 5:1 in Fig. 14.
(a)
(b) Crack front
UnCracked surface
N i"
Ni N i' Cracked surface
0.05a
0.05a
Fig. 12. The finite element model used to calculate G for the circular crack: (a) finite element model; (b) detailed mesh around the crack front.
245
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Fig. 13. Comparison of SIF between present study and the analytical solution.
The normalized SIF under the different ratios a calculated by presented study are illustrated in Fig. 15a. It is found that the computed SIF increased with the increasing element length ratios in Fig. 15a. When the element length ahead of the crack front is smaller, the computed stress intensity factors show slightly higher and vice versa. These trends agree well with the three dimensional solutions given in Okada [12] for through crack problem. For comparison purpose, the K results calculated by Xie and Biggers’ approach [16] are plotted in Fig. 15b. It is seen that the computed stress intensity factors decreased with the increasing element ratios in Fig. 15b. These trends are contrary to the observations in Fig. 15a and in Okada et al. [12]. The reason caused to these results maybe boiled down to the displacement opening calculation as discussed in Section 2.3. Nevertheless, the present modified approach provided more accuracy K results when dealing with the finite element meshes with different element length.
3.2. Semi-elliptical surface crack subjected to tensile loading The second example is the problem of semi-elliptical surface crack, as shown in Fig. 16. Due to the symmetry of geometry and loading condition, a quarter of the plate was modeled as well. Two mesh patterns were used for this case as shown in Fig. 17a and b. In Fig. 17a, the crack tip nodes lie precisely on the actual crack front and the mesh is non-orthogonal. The aim is to evaluate the ability of the presented VCCT process with the non-orthogonal mesh pattern. The element length across the crack front is set to be 0.05a. The total numbers of elements and nodes are 18,900 and 21,168 for the model in Fig. 15a,
UnCracked surface
Crack front N i"
Ni N
' i
Cracked surface
Δ2
Analysis cases No. 1 2 3 4 5
α = Δ1: Δ 2
1:5 3:5 5:5 5:3 5:1
Δ1
Fig. 14. Detailed mesh around the crack front with different element length.
246
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Fig. 15. Comparison of SIF between analytical solutions and FEM results with different element length: (a) using the present approach; (b) using Xie and Biggers’ approach.
σz
t
a h
2c Crack
A
A
L Section A-A c/L=0.25 a/t=1/3 h/L=1.5 a/c=0.5
Z Y X
σz Fig. 16. A semi-elliptical surface crack subject to tensile loading.
247
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
(a)
y UnCracked surface
Cracked surface
0.05a
a
x
c
0.05a
(b) Approximate crack front
y
UnCracked surface
a Cracked surface
x c Fig. 17. Detailed mesh around the surface crack front. (a) Exact elliptical crack front, non-orthogonal mesh; (b) approximate zigzag crack front, nonorthogonal mesh.
In Fig. 17b, the semi-elliptical surface crack was approximated to the zigzag line. The nodes on either side of the actual crack front are defined as crack tip nodes. The zigzag line in Fig. 17b approximates the smooth crack front curve as closely as possible. In order to diminish errors of geometrical shape simulations, the element length and width are set to be 0.025a. The total numbers of elements and nodes are 35,000 and 38,556, for the mode in Fig. 15b. Fig. 18 shows the normalized SIF for the surface crack with the mesh pattern in Fig. 17a, which are obtained by three different approaches. The results given by Newman and Raju SIF equation [23] are also plotted in Fig. 18. Because the crack tip
Fig. 18. Comparison of analytical solution and present study with the mesh (a) in Fig. 17.
248
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Fig. 19. Comparison of analytical solution, present study and Xie’s work, with the mesh (b) in Fig. 17.
nodes lie exactly on the actual crack front, the K results in Fig. 16 reveal a good correlation for this case with non-orthogonal mesh pattern as expected. The maximum deviation between present study and Newman and Raju’ solution is 2.9% with the non-orthogonal mesh pattern. It can be also observed that, the K results evaluated by present study have small deviation from the results obtained by Xie’s approach [16,17] and by Smith’s approach [14]. Their maximum difference for Xie’ approach and Smith’ approach are 0.50% and 0.31%, respectively. Fig. 19 shows the normalized SIF for this case with the mesh pattern in Fig. 15b. The K results calculated by Xie and Biggers [16,17] with the mesh pattern in Fig. 15b and by Newman and Raju SIF equation [23] are also plotted in Fig. 19. At the extremes of 2//p = 0 and 2//p = 1, the normalized K results obtained by different approaches have a very small deviation for each other. Only 1.4% and 1.1% errors between present study and Newman and Raju’ solutions are found for these two nodes, respectively. For the other nodes, the K results calculated by present study and Xie and Biggers [16,17] oscillate around the analytical solution as expected. The reason caused to this result is because of the locations of the nodes on the approximate zigzag crack front. The nodes are lie either outside or inside the curved front. However, by using present approach, the amplitude of the oscillations is obviously smaller than that for Xie’s approach as shown in Fig. 19. 3.3. Through crack with an elliptical crack front The application of the present approach to analyze moving crack front is discussed in this section, where the dynamic effects should be taken into consideration [24,25]. A through crack with an elliptical crack front was examined, in which
Fig. 20. Through crack with an elliptical crack front.
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
Crack front
Top surface
249
C2
t
bottom surface
C1
Fig. 21. Detailed mesh around a through crack with an elliptical crack front.
the crack shapes will change instantaneously during crack propagation. The geometry of the plate and the initial crack shapes are shown in Fig. 20. The through crack with initial shape was approximated to the zigzag line, as shown in Fig. 21. The tensile stress rZ = 150 MPa is applied along the Z-direction on plate. This is a typical mode I dominated crack. Thus, the fracture criterion can be simplified as
Ed ¼
GI GIC
x
1
ð21Þ
where x = 1 and GIC is assumed to be 7.144 N/mm in this work. Fig. 22 shows the normalized SIF for the through crack with the initial elliptical crack front, in which the zigzag mesh pattern shown in Fig. 21 are used. The K results calculated by Xie and Biggers [16,17] with the same mesh pattern are also plotted in Fig. 21. The K results calculated by present study and Xie [16,17] oscillate around the Fawaz’ solution [6] as expected. The amplitude of the oscillations for present study is smaller than that for Xie’s approach. It can be also found that, the normalized K has the maximum value for the nodes at the top surface. The K results are much lower for the nodes at the bottom surface. The crack shapes variations with crack propagation are illustrated in Fig. 23 using the proposed approach. The crack shape varied in two stages. At the first stage, only the nodes near top surface of the plate were break, as shown in Fig. 23a–e. The reason can be boiled down to the G distributions along the crack front, which had been illustrated in Fig. 22. With the crack propagation, the distance between the nodes at top surface and at bottom surface decreased. As a result, the elliptical crack front changes into a straight line, as shown in Fig. 23a–e. At the second stage, the nodes at the mid-thickness were firstly break, as shown in Fig. 23f–h. Because, the strain energy release rate G has the largest value at mid-thickness and drops at the free surface of the plate [12] with a through-thickness crack. The crack shape seems to a convex line at the mid-thickness, as shown in Fig. 23h. The results in Fig. 23 indicated that, the present approach can enable crack growth to be tracked as the shape changes. Since the K results oscillate around the analytical solution, the accuracy of present approach, the sensitivity study of modeling details and the software convergence are presented in authors’ following paper [26].
Fig. 22. The normalized SIF for the through crack with initial elliptical crack front using zigzag mesh pattern.
250
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
(c)
(b)
(a)
(e)
(d)
Step: 14
Step: 10
Step: 5
(f)
Step: 25
Step: 20
Step: 30
(h)
(g)
Step: 40
Step: 50 Fig. 23. Crack shapes under different steps.
4. Conclusions In this paper, a modified zigzag approach to approximate moving crack front with arbitrary shape was proposed based on the VCCT. Some useful conclusions are summarized as follows. 1. A modified zigzag approach to evaluate G was presented for arbitrary crack front using non-orthogonal mesh pattern. There are two major modifications on the Xie and Biggers’ approach in present study. Firstly, the present approach provided a more detailed consideration to determine the virtually closed area. There is no blank space between the virtually closed areas for the consecutive crack tip nodes. Secondly, the present approach can be used to calculate G when the finite elements across the crack front have different sizes. 2. Three cases were examined to study the validity and accuracy of the proposed approach. The G results calculated in present study show excellent correlation with the analytical solutions or the other established FEM solutions. For the zigzag pattern, the G results calculated by present study and Xie oscillate around the analytical solution. However, by using present study, the amplitude of the oscillations is obviously smaller than that for Xie’s approach. 3. In order to study the crack growth behavior, a through crack with an elliptical crack front was examined using the zigzag mesh pattern. Although the crack shape changes during crack propagation, a simple stationary finite element mesh was used without adaptive re-meshing. And the results indicated that the changes of crack shape can be tracked by using the present approach.
References [1] Irwin GR. Analysis of stresses and strain near the end of a crack traversing a plate. J Appl Mech 1957;24:361–4. [2] Rybicki EF, Kanninen MF. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng Fract Mech 1977;9:931–8. [3] Shivakumar KN, Tan PW, Newman Jr JC. A virtual crack-closure technique for calculating stress intensity factors for cracked three dimensional bodies. Int J Fract 1988;36:R43–50. [4] Krueger R. Virtual crack closure technique: history, approach, and applications. Appl Mech Rev 2004;57(2):109–43. [5] Xie D, Biggers SB. Progressive crack growth analysis using interface element based on the virtual crack closure technique. Finite Elem Anal Des 2006;42:977–84.
Y.-P. Liu et al. / Engineering Fracture Mechanics 78 (2011) 234–251
251
[6] Fawaz SA. Application of the virtual crack closure technique to calculate stress intensity factors for through cracks with an elliptical crack front. Eng Fract Mech 1998;59(3):327–42. [7] Fawaz SA. Stress intensity factor solutions for part-elliptical through cracks. Eng Fract Mech 1999;63:209–26. [8] Lin XB, Smith RA. A numerical simulation of fatigue growth of multiple surface initially semicircular defects under tension. Int J Pres Ves Piping 1995;62:281–9. [9] Lin XB, Smith RA. Fatigue growth analysis of interacting and coalescing surface crack. Int J Fract 1997;85:283–99. [10] Xie D, Wass AM, Shahwan KW, Schroeder JA, Beoman RG. Fracture criterion for kinking cracks in a tri-material adhesively bonded joint under mixed mode loading. Eng Fract Mech 2005;72:2487–504. [11] Xie D, Wass AM, Shahwan KW, Schroeder JA, Beoman RG. Computation of energy release rates for kinking cracks based on virtually crack closure technique. CNES-Comput Model Engng Sci 2004;6:515–24. [12] Okada H, Higashi M, et al. Three dimensional virtual crack closure-integral method (VCCM) with skewed and non-symmetric mesh arrangement at the crack front. Engng Fract Mech 2005;72:1717–37. [13] Okada H, Endoh S, Kikuchi M. On fracture analysis using an element overlay technique. Engng Fract Mech 2005;72:773–89. [14] Smith SA, Raju IS. Evaluation of stress-intensity factors using general finite-element models. Fatigue Fract Mech. In: Panontin TL, Sheppard SD, editors. ASTM STP 1332, vol. 29. West Canshohockem (PA): American Society for testing and materials; 1999, p. 176–200. [15] Ferrie CH, Rousseau CQ. A method of applying VCCT to corner crack nodes. In: Proc of Am Soc for composites, 16th annual technical conference on composite materials, Blacksburg, VA, USA; 2001. [16] Xie D, Biggers SB. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique. Part I: formulation and validation. Engng Fract Mech 2006;73:771–85. [17] Xie D, Biggers SB. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique. Part II: sensitivity study on modeling details. Engng Fract Mech 2006;73:786–801. [18] Xie D, Gary M, Huang DD, Abdi F. Cohesive zone model for surface cracks using finite element analysis. AIAA-2008-106742, Chicago, Illinois; April 2008. [19] Xie D, Waas AM. Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Engng Fract Mech 2006;73:1783–96. [20] Boscolo M, Allegri G, Zhang X. Design and modeling of selective reinforcements for integral aircraft structures. AIAA J 2009;46(9):2323–31. [21] Kutlu Z, Chang FK. Modeling compressive failure of laminated composite containing multiple through-the-width delaminations. J Compos Mater 1992;26:350–87. [22] Murakami Y et al. Stress intensity factors handbook. Berlin: Pergamon Press; 1987. [23] Newman JC, Raju IS. An empirical stress-intensity factor equation for the surface crack. Engng Fract Mech 1981;15(1–2):185–92. [24] Chen XD, He XM. The effect of the recess shape on performance analysis of the gas-lubricated bearing in optical lithography. Tribol Int 2006;39(11):1336–41. [25] He XM, Chen XD. The dynamic analysis of the gas lubricated stage in optical lithography. Int J Adv Manufact Tech 2007;32(9–10):978–84. [26] Liu YP. Study on fatigue crack growth behavior in welded steel plates used for bridge. Ph.D thesis. China: Huazhong University of Science and Technology; Aug 2010.