- Email: [email protected]
Three-dimensional crack growth simulation using BEM
Pergamon
00457949(94)Eo268-7
Compufers & Swucrures Vol. 52. No. 5. pp. W-878, 1994 Copyright 0 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 004%7949/94 $7.00 + 0.00
THREE-DIMENSIONAL CRACK GROWTH SIMULATION USING BEM Y. Mif and M. H. Aliabadi Wessex Institute of Technology, Portsmouth University, Ashurst Lodge, Southampton SO4 2AA, U.K. (Received 25 May 1993) Abstract-In this paper an application of dual boundary element method to the analysis of three-dimensional mixed-mode crack growth is presented. The crack growth processes are simulated numerically with an incremental crack-extension analysis based on the minimum strain energy density criterion and a fatigue crack growth law. For each crack extension, the size of the increment (crack growth) and direction of crack growth is evaluated along the crack front. The dual boundary element analysis is applied to perform a single-region analysis and the stress intensity factors are evaluated using a displacement-based method. The crack extension is modelled with the introduction of new boundary elements, removing the requirement of interior remeshing, an intrinsic feature of the dual boundary element method. Results of analysis are presented for several geometries, including planar and non-planar crack growth.
INTRODUCTION
Numerical methods for the reliability analysis of structural components have accomplished a level of development which makes them viable tools for design engineers. In particular, through advances in the application of finite and boundary element methods to fracture mechanics, it is possible to address critical problems such as life expectancy of flawed or cracked structures. The simulation of general mixed-mode crack growth using numerical techniques requires the capability of predicting the direction and the amount of crack growth for each given load increment as well as the robustness to update the numerical model to account for the chang-
ing crack geometry. Finite element applications to crack growth analysis can be found in Refs [l-5] using techniques such as node grafting, local remeshing and automatic remeshing. However, an intrinsic feature of the finite element method common to all of these formulations is the need for continuous remeshing to follow the crack growth. Application of the boundary element method to fracture mechanics is now well established (see Aliabadi and Rooke [6]) and was first applied to mode I crack growth analysis by Cruse and Besuner [7] in 1975. It is well known that the solution of general crack problem is not possible with direct application of the BEM in a single region analysis. This is due to the coincidence of the crack surfaces which gives rise to a singular system of algebraic equations. Blandford et al. [8] used the subregion method to overcome this difficulty. The subregion method was subsequently
ton
leave from Taiyuan Heavy Machinery Institute,
Taiyuan, People’s Republic of China.
used to simulate two- and three-dimensional crack growth [9, IO]. The subregion method introduces artificial boundaries into the body which connect the crack front to the boundary in such a way that the domain is divided into subregions without cracks. In an incremental crack extension analysis, these artificial boundaries must be repeatedly re-introduced for each increment of the crack extension. Hence the main drawback of this method is that the introduction of artificial boundaries is not unique and thus cannot be easily implemented in an automatic procedure for incremental analysis of crack extension. Furthermore, the method generates a larger system of algebraic equations (particularly in 3D) than is required. These drawbacks have led to the development of alternative formulations [ll-161 to remove the requirement of subregions. Among these alternative formulations the dual boundary element method developed by Portela et al. [ 151 for two-dimensional and Mi and Aliabadi [16] for three-dimensional crack problems is the most general and has been applied to problems involving embedded cracks, edged cracks and kinked cracks. The dual boundary element method incorporates two independent boundary integral equations, with the displacement equation applied for collocation on one of the crack surfaces and the traction on the other. As a consequence, general mixed-mode crack problems can be solved in a single region. Application of the dual boundary element method to two-dimensional crack growth analysis has been presented in Ref. [17] and the extension of the method to dynamic and elastoplastic analysis are reported in references (181 and [19], respectively. In this paper the application of the DBEM [16] to analysis of three-dimensional crack growth in linear elastic fracture
871
Y. Mi and M. H. Aliabadi
872
mechanics is presented. An incremental crack extension analysis is performed to determine the crack path. The minimum strain energy density criterion [20] is applied to obtain both the size of the crack extension increment and the crack growth direction. Also presented is a modified fatigue criterion based on the Paris formula to obtain the size of the crack extension. Results of incremental crack extension analysis are presented for several geometries including planar and non-planar crack growth. THE DUAL BOUNDARY
ELEMENT
METHOD
The dual boundary integral equations, on which the DBEM is based, are the displacement and the traction integral equations. The displacement boundary integral equation relating the boundary displacement components uj and traction components $ can be written as [21]
smooth boundary, the stress components gii are given by
fcgx’) +
T&‘, X)UL (x) dr 6)
=
U&‘, f Jr
f tj(x’) + n,(x’)
T&‘,
f
U&‘,
x)tk(x) dr(x),
(3)
)-IFr
I=
xh (4 dr (x)
f r
T&x’, x)u,(x) dT(x)
+
(2)
where Tkij(x', x) and Ukii(x’, x) contain derivatives of Tj(x’, x) and U,(x’, x) respectively. The symbol f stands for the Hadamard principal value integral. On a smooth boundary, the traction components, tj are given from eqn (2) by
= n,(x’ C,(X’)llj(X’)
xh(x) dr(x),
Uo(X’,X)?,(X)dT(x),
(1)
where i, j denote Cartesian components, Tu(x', x)and UJx‘, x) represent the Kelvin traction and displacement fundamental solutions, respectively, at a boundary point x. The symbol f stands for the Cauchy principal value integral. In absence of body force and assuming continuity of both strains and traction at x’ on a
where n,(x’) denotes the component of outward unit normal to the boundary at x’. The boundary integral equations (1) and (3) constitute the dual boundary element method. The integrals in eqns (1) and (3) are regular, provided r # 0 (r = lx’ - xl) in the fundamental solutions. As the distance r tends to zero, the fundamental solutions exhibit singularities, they are of order l/r in U,, l/r2 in To and Ukij and l/r’ in Tkv.Accurate integration procedures for these singular integrals have been reported in Ref. [16].
discontinuous element
Crack
edge
discontinuous
continuous
Fig. 1. Crack modelling strategy.
Crack growth simulation using BEM
O-Geometry x -Collocation
(4
nodes w
nodes
Fig. 2. Discontinuous and edge discontinuous elements.
CRACK MODELLINC
STRATEGY
for nodes at 5, = 0, r~.= f 1; and
The use of dual boundary integral equations imposes certain restrictions on the choice of elements required for the discretization of crack surfaces. These restrictions are due to the continuity requirements of the field variables for the existence of Cauchy and Hadamard principal value integrals. For this reason and to maintain efficiency and simplicity of the boundary elements the present formulation uses discontinuous quadratic elements for the crack modelling, as shown in Fig. 1. The geometry parameter xi is approximated using quadratic quadrilateral serendipity elements. The shape functions are given as
N’(C,?)=~(d+55.)(~*-~*).
(5)
for nodes at {, = f 1, qU= 0, where 1 is the parametric position of collocation points and 0 < 1 d 1 (see Fig. 2). On parts of the surface intersecting the crack surface, the displacement and traction fields are approximated using an ‘edge discontinuous’ element. The shape functions for an element, depicted in Fig. 2(b), are given as N,(5 rl) =
(1- 00
- 7)(-T - rl - 1) ’ 2(A + 1)
MY59 rl) = %l + X.)(1 + ?rl.)(rL + VIZ- l), N2(r
q)
=
1
for nodes at &, q. = f 1; MY<, rl) = $1 - <*)(I + VI.),
N3({
,
(1- 5*)Q- rl) ’
(n+l) q)
=
t1 + t)(’ - q)(t - q - 1) ’
2(1 + I)
for nodes at 5, = 0, r~.= f 1; and MY57 rt) = %l + er,,<1 - tt2),
(4)
N4(L4) =
(1 +q)Q -rl)(l 21
+o ’
for nodes at 5, = + 1, q. = 0. The displacement and traction fields over the crack surface are approximated using quadratic discontinuous elements (see Fig. 2a). The shape functions are given as
for nodes at {,, qDl= f 1;
Fig. 3. Crack front coordinate.
874
Y. Mi and M. H. Aliabadi
h a=Zb R=lOa h=6R
Fig. 4. An embedded elliptical crack in a cylindrical bar.
N5(5
rl)
=
(1 + 5)U + rlwt +
rl -
N6(5,rl) = N7(5
(1 -
1) ’
21(1 +A)
t21u+ II) ’
(1+1)
tl)=(1-5)(1+11)(--lr+‘1-1) ’
21(1+X?)
Ns(Lv I=
(1 +?I@ -rl)U 21
-5) 1
(6)
where the parameter i is the same as in eqn (5). On all other surfaces, shape functions in eqns (5) or (6) are used with 1 = 1 (i.e. N”(& 9) = M”(& n)). The general modelling strategy developed here is similar to that reported in Ref. [16] and can be summarized as follows: (i) crack surfaces are modelled with discontinuous quadratic elements;
iz Fig. 5. Propagated crack fronts for an embedded elliptical crack (4 = 0).
(ii) surfaces intersecting a crack surface are modelled with edge-discontinuous elements; (iii) all other surfaces are modelled with continuous elements; (iv) the displacement integral equation (1) is applied for collocation on one of the crack surfaces (say the upper surface); (v) the traction integral equation (3) is applied on the opposite crack surface (say the lower surface); (iv) the displacement integral equation (1) is used for collocation on all other surfaces. CRACK EXTENSIONANALYSIS
Several criteria have been proposed to describe the local direction of mixed-mode crack growth. Among them, the most popular are the maximum principal stress, the maximum energy release rate and the minimum strain energy density. Here, the minimum strain energy density criterion is adopted, which was
o.700’1”~“~~~~‘~r”““~~““‘~“~~~~’ so eo 30
120
Isa
110
Front Fig. 6. Normalized stress intensity factors for the embedded elliptical crack (4 = 0). a:
Angle
Along
the
Crack
875
Crack growth simulation using BEM
-
Hypothesis (2): crack extension occurs when the strain energy density factor in the region determined by hypothesis (l), S = Smin,reaches a critical value, say S,,. Hypothesis (3): the length, r,,, of the initial crack extension is assumed to be proportional to Smiosuch that S,,/r,, remains constant along the new crack front.
2h
Y i I
X a
!!!ifP-
I I
,(----___ ,’
-
:‘i
m
Fig. 7. A quarter penny-shaped edge crack in a square bar.
formulated by Sih [20]. The theory is based on the hypotheses. Hypothesis (1): the direction of crack propagation at any point along the crack front is toward the region with the minimum value of strain energy density factor S as compared with other regions on the same spherical surface surrounding the point.
The incremental crack extension analysis assumes a piece-wise linear discretization of the unknown crack path. For each increment of the crack extension, the DBEM is applied to carry out a stress analysis of cracked struture and displacement extrapolation technique [16] is used for the evaluation of the mixed-mode stress intensity factors. The steps of the computational cycle are summarized as follows: (i) carry out a DBEM stress analysis of cracked structure; (ii) compute the stress intensity factors along the crack front; (iii) compute the direction of crack extension increment using Sih’s strain energy criterion; (iv) compute the amount of increment along the crack front using Sib’s criterion or Paris fatigue growth formula; (v) repeat all the above steps sequentially until a specified number of crack extension increments is reached. CRACK EXTENSION
Three applications of the incremental crack extension analysis will be studied in present section. The first and third applications characterized the behaviour of embedded cracks deformed under pure
Y
4th
incr.
front
ti
Al;= 0.1497 Lig= a.1 003
3rd
incr.
front
4
2nd
incr.
front
+
AS= 0.1811
1st
incr.
front
*
front
I,
Ax,= 0.2215 AXZ=0.2002 Axs= 0.2049 Aa = 0.2032
Original
APPLICATION
4%5= 0.1800
Fig. 8. Propagated crack fronts for an edge crack configuration.
Y. Mi and M. H. Aliabadi
816
mode I and mixed-mode. The second application analyses a corner crack. For the first and second application, a quasi-static crack propagation is assumed. Under this assumption, the strain energy density factors S, after each analysis, along the crack front, are evaluated by
S=all~~+2a,~KlK,,+a,,K~,+a3,~~,,,
(7)
where, K,, K,, and K,,, are mode I, II and III stress intensity factors and n,,=&(3-4v-
cos @(l + cos (3)
a,,=&sin8(cosB-1+2v)
C$between the y-axis and the crack surface is taken as 0 so that only the incremental size of the crack extension need be considered. The ratios b/a = 0.5, R/a = 10 and h/R = 6 are chosen and the Young’s modulus is taken to be 1000 units and the Poisson’s ratio is 0.3. A total of five steps of crack extensions were performed and the maximum incremental size of the crack at each state was fixed at 0.2 times the semi-major axis a. Figure 5 shows the propagated crack fronts together with the original one. The normalized stress intensity factors for each crack extension are plotted in Fig. 6. As it can be seen that the crack grows into the shape (i.e. circular) which gives constant stress intensity factors along the crack front. The second example is a square bar specimen of width w and height 2h containing a quarter pennyshaped corner crack of radius a as shown in Fig. 7. The ratios of a/w and a/h are chosen as l/4. This corner crack configuration is a more challenging
1 a33 = G
t
I
t
t
in which p stands for the shear modulus of elasticity, v is Poisson’s ratio and 6 is angle in the crack front coordinate system, as shown in Fig. 3. The crack incremental direction (i.e. 8 in Fig. 3), at each point along the front, is obtained through minimizing the strain energy density factor S of eqn (7) with respect to 0. The crack incremental size r, is decided by keeping Smin/r,,= constant, as stated in the previous section. For the first application, an elliptical crack of the major semi-axis a and minor semi-axis b in a cylindrical bar of radius R and height h subjected to a tensile stress u, as shown in Fig. 4, is considered. The angle t
t
t
t
t
-
Initial
crack
-
Propagated
mesh
crack
t
(b)
I.20
* 0
10
20
w : Angle
54
40
Along
50
the
60
Crack
70
10
90
Front
Fig. 9. Normalized stress intensity factors for the edge crack problem.
Fig. 10. (a) Propagated crack fronts. (b) Deformed crack mesh after fourth crack extension.
Crack growth simulation using BEM problem for assessing the accuracy of the proposed procedure as the crack boundary interactions influence the crack front profile. The same material properties are used as in the first example. The crack fronts of five stages are projected on the original crack plane and shown in Fig. 8; the stress intensity factors along the crack front, corresponding to each stage (excluding the last one), are plotted in Fig. 9, from which similar behaviour to the first example can be seen. The embedded elliptical crack problem of the first example is considered again as an example, but the crack surface makes an angle 4 = 45” with the y-axis, as shown in Fig. 4. This is a mixed problem and the fatigue crack growth is considered instead of quasi-static one. The Paris model is used to model fatigue crack propagation, which states that da = c(AKeff)“,
(12) emin is chosen as zero and a,,,,, constant in this example. After each DBEM analysis, the directions of the next increment of the crack growth are evaluated by minimizing Sih’s strain energy density factor. The maximum amount Au, of the next increment is fixed (0.2u), corresponding to the crack front point where maximum AKF occurs. From eqn (9)
AK,,=K,“,-K$=‘=KF(l-R),
Aa x c(AK~)~ AN
(13)
z c(AKz)”
(14)
and Aa,
AN.
From eqns (13) and (14), the following relationship is obtained
where da/dN is the rate of change of crack length with respect to the loading cycles, c = 1.5463 x lo-l6 and m = 3.88 are material constants; (10)
Therefore, the incremental size at any front point can be evaluated by
in which X;, stands for the effective stress intensity factor, adopted here as
f& = WI + Ku, 1)’+ 2K:,.
(11)
V-
LY N
Since linear elasticity is considered, R in eqn (10) can be written as
(9)
dN
811
B
_r
-0.2
L
!-0.6 -
FW
/6
I 0 Original
p/
-0.8 0
1
0 lncrimint I
I
I
20
40
60
w
,
80
: ANGLE
I
100
I
120
ALONG
4
I
,
140
THE
160
I
180
CRACK
I
,
ZOO
220
,
240
I
260
,
280
,
3DO
,
320
FRONT
Fig. 11. Normalized stress intensity factors for inclined embedded elliptical crack.
,
2.40
360
878
Y. Mi and M. H. Aliabadi
A total of five steps of increment are mrformed. The crack surface after the last analysis, including five crack fronts, is shown in Fig. 10(a) and the deformed
crack mesh after the last analysis in Fig. 10(b). Model I, II and III stress intensity factors corresponding to each analysis are plotted in Fig. 11.
CONCLUSIONS
An automatic procedure for simulation of threedimensional crack growth using the dual boundary element method has been presented in this paper. The method has the advantage of eliminating interior remeshing as the crack growths take place, hence resulting in a substantial amount of computational savings as well as simple modelling. The test examples demonstrate that, through using the minimum strain energy criterion and a modified fatigue formula, cracks grow into a shape with constant stress intensity factors. REFERENCES I.
2.
3.
4.
5.
6. 7.
A. R. Ingraffea, Nodal grafting for crack propagation studies. &t. J. numer. Eigng 11, 1185-1187 (i977). R. Perucchio and A. R. Ineraffea. An integrated boundary element analysis system with interaczve computer graphics for three-dimensional linear-elastic fracture mechanics. Compur. Struct. 20, 157-171 (1985). S. Valliappan and V. Murt, Automatic remeshing technique in quasi-static and dynamic crack propagation. Proc. NCJMETA, ‘85 ConJ, pp. 107-116, Swansea (1985). M. S. Shepherd, N. A. S. Yehia, G. S. Burd and T. J. Weinder, Automatic crack propagation tracking. Comput. Struct. 20, 211-223 (1985). E. M. Remzi and W. S. Blackburn, Automatic crack propagation studies in T-junctions and cross bars. Engng Compur. 7, 116-124 (1990). M. H. Aliabadi and D. P. Rooke, Numerical Fracture Mechanics. Computational Mechanics Publications, Southampton and Kluwer Academic, Dordrecht (1991). T. A. Cruse and P. M. Besuner, Residual life prediction for surface cracks in complex structural details. J. Aircraff 12, 369-375 (1975). I
I
8. G. Blandford, A. R. Ingraffea and J. A. Liggett, Two dimensional stress intensity factor computation using the boundary element method. Inr. J. numer. Engng 17, 387-404 (1981). 9 A. R. Ingraffea, G. Blandford and J. A. Liggett, Automatic modelling of mixed-mode fatigue and quasi-static crack propagation using the boundary element method. Proc. 14th Nat. Symp. on Fracture Mech. ASTM STP 791, I 407-I 426 (1987). 10. P. A. Wawrzynek, L. F. Martha and A. R. Ingraffea, A computational environment for the simulation of fracture process in three-dimensions. In Annul. Numer. Exper. Aspects of Three-dimensional Fracture Process
(Edited by A. J. Rosakis et al.), ASME, AMD-91, pp. 321-327 (1988). 11. J. Weaver, Three-dimensional crack analysis. Znr. J. Solids Struct. 13. 321-330 (1977). 12. J. 0. Watson, Hermitian cubic anh singular elements for plain strain. In Development in Boundary Element Merhod. Vol. 4 (Edited bv P. K. Baneriiee and T. 0. Watsonj. Elsevie; Applieci Science, Barcing (1986). 13. E. Z. Polch, T. A. Cruse and C. J. Huang, Traction BIE solution for flat cracks. Comput. Mech. 2, 253-267 (1987). 14. L. J. Gray,
L. F. Martha and A. R. Integraffea, Hypersingular integrals in boundary element fracture analysis. Znt. J. numer. Merh. Engng 29, 1135-1158, (1990). 15. A. Portela, M. H. Aliabadi and D. P. Rooke, The dual boundary element method: effective implementation for crack problems. Znr. J. numer. Meth. Engng 33, 1269-1287
(1992).
16. Y. Mi and M. H. Aliabadi, Dual boundary element method for three-dimensional fracture mechanics analysis. Engng Anal. Boundary Elements 10, 161-171 (1992). 17. A. Portela, M. H. Aliabadi and D. P. Rooke, Dual boundary element incremental analysis of crack propagation. Comput. Struct. 46, 237-247 (1993). 18. P. Fedelinski, M. H. Aliabadi and D. P. Rooke, The dual boundary element method in dynamic fracture mechanics. Engng. Anal. Boundary Elements 12, 203-210 (1993).
19. V. Leitlo, M. H. Aliabadi and D. P. Rooke, The dual boundary element formulation for elastoplastic fracture mechanics, in press. 20. G. C. Sih, Mechanics of Fracture Initiation and Propagation. Kluwer Academic, Dordrecht (1991). 21. T. A. Cruse, Mathematical formulation of the boundary integral equation methods in solid mechanics. Report No. AFOSR-TR-77-1002, Pratt and Whitney Aircraft Group (1977).
00457949(94)Eo268-7
Compufers & Swucrures Vol. 52. No. 5. pp. W-878, 1994 Copyright 0 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 004%7949/94 $7.00 + 0.00
THREE-DIMENSIONAL CRACK GROWTH SIMULATION USING BEM Y. Mif and M. H. Aliabadi Wessex Institute of Technology, Portsmouth University, Ashurst Lodge, Southampton SO4 2AA, U.K. (Received 25 May 1993) Abstract-In this paper an application of dual boundary element method to the analysis of three-dimensional mixed-mode crack growth is presented. The crack growth processes are simulated numerically with an incremental crack-extension analysis based on the minimum strain energy density criterion and a fatigue crack growth law. For each crack extension, the size of the increment (crack growth) and direction of crack growth is evaluated along the crack front. The dual boundary element analysis is applied to perform a single-region analysis and the stress intensity factors are evaluated using a displacement-based method. The crack extension is modelled with the introduction of new boundary elements, removing the requirement of interior remeshing, an intrinsic feature of the dual boundary element method. Results of analysis are presented for several geometries, including planar and non-planar crack growth.
INTRODUCTION
Numerical methods for the reliability analysis of structural components have accomplished a level of development which makes them viable tools for design engineers. In particular, through advances in the application of finite and boundary element methods to fracture mechanics, it is possible to address critical problems such as life expectancy of flawed or cracked structures. The simulation of general mixed-mode crack growth using numerical techniques requires the capability of predicting the direction and the amount of crack growth for each given load increment as well as the robustness to update the numerical model to account for the chang-
ing crack geometry. Finite element applications to crack growth analysis can be found in Refs [l-5] using techniques such as node grafting, local remeshing and automatic remeshing. However, an intrinsic feature of the finite element method common to all of these formulations is the need for continuous remeshing to follow the crack growth. Application of the boundary element method to fracture mechanics is now well established (see Aliabadi and Rooke [6]) and was first applied to mode I crack growth analysis by Cruse and Besuner [7] in 1975. It is well known that the solution of general crack problem is not possible with direct application of the BEM in a single region analysis. This is due to the coincidence of the crack surfaces which gives rise to a singular system of algebraic equations. Blandford et al. [8] used the subregion method to overcome this difficulty. The subregion method was subsequently
ton
leave from Taiyuan Heavy Machinery Institute,
Taiyuan, People’s Republic of China.
used to simulate two- and three-dimensional crack growth [9, IO]. The subregion method introduces artificial boundaries into the body which connect the crack front to the boundary in such a way that the domain is divided into subregions without cracks. In an incremental crack extension analysis, these artificial boundaries must be repeatedly re-introduced for each increment of the crack extension. Hence the main drawback of this method is that the introduction of artificial boundaries is not unique and thus cannot be easily implemented in an automatic procedure for incremental analysis of crack extension. Furthermore, the method generates a larger system of algebraic equations (particularly in 3D) than is required. These drawbacks have led to the development of alternative formulations [ll-161 to remove the requirement of subregions. Among these alternative formulations the dual boundary element method developed by Portela et al. [ 151 for two-dimensional and Mi and Aliabadi [16] for three-dimensional crack problems is the most general and has been applied to problems involving embedded cracks, edged cracks and kinked cracks. The dual boundary element method incorporates two independent boundary integral equations, with the displacement equation applied for collocation on one of the crack surfaces and the traction on the other. As a consequence, general mixed-mode crack problems can be solved in a single region. Application of the dual boundary element method to two-dimensional crack growth analysis has been presented in Ref. [17] and the extension of the method to dynamic and elastoplastic analysis are reported in references (181 and [19], respectively. In this paper the application of the DBEM [16] to analysis of three-dimensional crack growth in linear elastic fracture
871
Y. Mi and M. H. Aliabadi
872
mechanics is presented. An incremental crack extension analysis is performed to determine the crack path. The minimum strain energy density criterion [20] is applied to obtain both the size of the crack extension increment and the crack growth direction. Also presented is a modified fatigue criterion based on the Paris formula to obtain the size of the crack extension. Results of incremental crack extension analysis are presented for several geometries including planar and non-planar crack growth. THE DUAL BOUNDARY
ELEMENT
METHOD
The dual boundary integral equations, on which the DBEM is based, are the displacement and the traction integral equations. The displacement boundary integral equation relating the boundary displacement components uj and traction components $ can be written as [21]
smooth boundary, the stress components gii are given by
fcgx’) +
T&‘, X)UL (x) dr 6)
=
U&‘, f Jr
f tj(x’) + n,(x’)
T&‘,
f
U&‘,
x)tk(x) dr(x),
(3)
)-IFr
I=
xh (4 dr (x)
f r
T&x’, x)u,(x) dT(x)
+
(2)
where Tkij(x', x) and Ukii(x’, x) contain derivatives of Tj(x’, x) and U,(x’, x) respectively. The symbol f stands for the Hadamard principal value integral. On a smooth boundary, the traction components, tj are given from eqn (2) by
= n,(x’ C,(X’)llj(X’)
xh(x) dr(x),
Uo(X’,X)?,(X)dT(x),
(1)
where i, j denote Cartesian components, Tu(x', x)and UJx‘, x) represent the Kelvin traction and displacement fundamental solutions, respectively, at a boundary point x. The symbol f stands for the Cauchy principal value integral. In absence of body force and assuming continuity of both strains and traction at x’ on a
where n,(x’) denotes the component of outward unit normal to the boundary at x’. The boundary integral equations (1) and (3) constitute the dual boundary element method. The integrals in eqns (1) and (3) are regular, provided r # 0 (r = lx’ - xl) in the fundamental solutions. As the distance r tends to zero, the fundamental solutions exhibit singularities, they are of order l/r in U,, l/r2 in To and Ukij and l/r’ in Tkv.Accurate integration procedures for these singular integrals have been reported in Ref. [16].
discontinuous element
Crack
edge
discontinuous
continuous
Fig. 1. Crack modelling strategy.
Crack growth simulation using BEM
O-Geometry x -Collocation
(4
nodes w
nodes
Fig. 2. Discontinuous and edge discontinuous elements.
CRACK MODELLINC
STRATEGY
for nodes at 5, = 0, r~.= f 1; and
The use of dual boundary integral equations imposes certain restrictions on the choice of elements required for the discretization of crack surfaces. These restrictions are due to the continuity requirements of the field variables for the existence of Cauchy and Hadamard principal value integrals. For this reason and to maintain efficiency and simplicity of the boundary elements the present formulation uses discontinuous quadratic elements for the crack modelling, as shown in Fig. 1. The geometry parameter xi is approximated using quadratic quadrilateral serendipity elements. The shape functions are given as
N’(C,?)=~(d+55.)(~*-~*).
(5)
for nodes at {, = f 1, qU= 0, where 1 is the parametric position of collocation points and 0 < 1 d 1 (see Fig. 2). On parts of the surface intersecting the crack surface, the displacement and traction fields are approximated using an ‘edge discontinuous’ element. The shape functions for an element, depicted in Fig. 2(b), are given as N,(5 rl) =
(1- 00
- 7)(-T - rl - 1) ’ 2(A + 1)
MY59 rl) = %l + X.)(1 + ?rl.)(rL + VIZ- l), N2(r
q)
=
1
for nodes at &, q. = f 1; MY<, rl) = $1 - <*)(I + VI.),
N3({
,
(1- 5*)Q- rl) ’
(n+l) q)
=
t1 + t)(’ - q)(t - q - 1) ’
2(1 + I)
for nodes at 5, = 0, r~.= f 1; and MY57 rt) = %l + er,,<1 - tt2),
(4)
N4(L4) =
(1 +q)Q -rl)(l 21
+o ’
for nodes at 5, = + 1, q. = 0. The displacement and traction fields over the crack surface are approximated using quadratic discontinuous elements (see Fig. 2a). The shape functions are given as
for nodes at {,, qDl= f 1;
Fig. 3. Crack front coordinate.
874
Y. Mi and M. H. Aliabadi
h a=Zb R=lOa h=6R
Fig. 4. An embedded elliptical crack in a cylindrical bar.
N5(5
rl)
=
(1 + 5)U + rlwt +
rl -
N6(5,rl) = N7(5
(1 -
1) ’
21(1 +A)
t21u+ II) ’
(1+1)
tl)=(1-5)(1+11)(--lr+‘1-1) ’
21(1+X?)
Ns(Lv I=
(1 +?I@ -rl)U 21
-5) 1
(6)
where the parameter i is the same as in eqn (5). On all other surfaces, shape functions in eqns (5) or (6) are used with 1 = 1 (i.e. N”(& 9) = M”(& n)). The general modelling strategy developed here is similar to that reported in Ref. [16] and can be summarized as follows: (i) crack surfaces are modelled with discontinuous quadratic elements;
iz Fig. 5. Propagated crack fronts for an embedded elliptical crack (4 = 0).
(ii) surfaces intersecting a crack surface are modelled with edge-discontinuous elements; (iii) all other surfaces are modelled with continuous elements; (iv) the displacement integral equation (1) is applied for collocation on one of the crack surfaces (say the upper surface); (v) the traction integral equation (3) is applied on the opposite crack surface (say the lower surface); (iv) the displacement integral equation (1) is used for collocation on all other surfaces. CRACK EXTENSIONANALYSIS
Several criteria have been proposed to describe the local direction of mixed-mode crack growth. Among them, the most popular are the maximum principal stress, the maximum energy release rate and the minimum strain energy density. Here, the minimum strain energy density criterion is adopted, which was
o.700’1”~“~~~~‘~r”““~~““‘~“~~~~’ so eo 30
120
Isa
110
Front Fig. 6. Normalized stress intensity factors for the embedded elliptical crack (4 = 0). a:
Angle
Along
the
Crack
875
Crack growth simulation using BEM
-
Hypothesis (2): crack extension occurs when the strain energy density factor in the region determined by hypothesis (l), S = Smin,reaches a critical value, say S,,. Hypothesis (3): the length, r,,, of the initial crack extension is assumed to be proportional to Smiosuch that S,,/r,, remains constant along the new crack front.
2h
Y i I
X a
!!!ifP-
I I
,(----___ ,’
-
:‘i
m
Fig. 7. A quarter penny-shaped edge crack in a square bar.
formulated by Sih [20]. The theory is based on the hypotheses. Hypothesis (1): the direction of crack propagation at any point along the crack front is toward the region with the minimum value of strain energy density factor S as compared with other regions on the same spherical surface surrounding the point.
The incremental crack extension analysis assumes a piece-wise linear discretization of the unknown crack path. For each increment of the crack extension, the DBEM is applied to carry out a stress analysis of cracked struture and displacement extrapolation technique [16] is used for the evaluation of the mixed-mode stress intensity factors. The steps of the computational cycle are summarized as follows: (i) carry out a DBEM stress analysis of cracked structure; (ii) compute the stress intensity factors along the crack front; (iii) compute the direction of crack extension increment using Sih’s strain energy criterion; (iv) compute the amount of increment along the crack front using Sib’s criterion or Paris fatigue growth formula; (v) repeat all the above steps sequentially until a specified number of crack extension increments is reached. CRACK EXTENSION
Three applications of the incremental crack extension analysis will be studied in present section. The first and third applications characterized the behaviour of embedded cracks deformed under pure
Y
4th
incr.
front
ti
Al;= 0.1497 Lig= a.1 003
3rd
incr.
front
4
2nd
incr.
front
+
AS= 0.1811
1st
incr.
front
*
front
I,
Ax,= 0.2215 AXZ=0.2002 Axs= 0.2049 Aa = 0.2032
Original
APPLICATION
4%5= 0.1800
Fig. 8. Propagated crack fronts for an edge crack configuration.
Y. Mi and M. H. Aliabadi
816
mode I and mixed-mode. The second application analyses a corner crack. For the first and second application, a quasi-static crack propagation is assumed. Under this assumption, the strain energy density factors S, after each analysis, along the crack front, are evaluated by
S=all~~+2a,~KlK,,+a,,K~,+a3,~~,,,
(7)
where, K,, K,, and K,,, are mode I, II and III stress intensity factors and n,,=&(3-4v-
cos @(l + cos (3)
a,,=&sin8(cosB-1+2v)
C$between the y-axis and the crack surface is taken as 0 so that only the incremental size of the crack extension need be considered. The ratios b/a = 0.5, R/a = 10 and h/R = 6 are chosen and the Young’s modulus is taken to be 1000 units and the Poisson’s ratio is 0.3. A total of five steps of crack extensions were performed and the maximum incremental size of the crack at each state was fixed at 0.2 times the semi-major axis a. Figure 5 shows the propagated crack fronts together with the original one. The normalized stress intensity factors for each crack extension are plotted in Fig. 6. As it can be seen that the crack grows into the shape (i.e. circular) which gives constant stress intensity factors along the crack front. The second example is a square bar specimen of width w and height 2h containing a quarter pennyshaped corner crack of radius a as shown in Fig. 7. The ratios of a/w and a/h are chosen as l/4. This corner crack configuration is a more challenging
1 a33 = G
t
I
t
t
in which p stands for the shear modulus of elasticity, v is Poisson’s ratio and 6 is angle in the crack front coordinate system, as shown in Fig. 3. The crack incremental direction (i.e. 8 in Fig. 3), at each point along the front, is obtained through minimizing the strain energy density factor S of eqn (7) with respect to 0. The crack incremental size r, is decided by keeping Smin/r,,= constant, as stated in the previous section. For the first application, an elliptical crack of the major semi-axis a and minor semi-axis b in a cylindrical bar of radius R and height h subjected to a tensile stress u, as shown in Fig. 4, is considered. The angle t
t
t
t
t
-
Initial
crack
-
Propagated
mesh
crack
t
(b)
I.20
* 0
10
20
w : Angle
54
40
Along
50
the
60
Crack
70
10
90
Front
Fig. 9. Normalized stress intensity factors for the edge crack problem.
Fig. 10. (a) Propagated crack fronts. (b) Deformed crack mesh after fourth crack extension.
Crack growth simulation using BEM problem for assessing the accuracy of the proposed procedure as the crack boundary interactions influence the crack front profile. The same material properties are used as in the first example. The crack fronts of five stages are projected on the original crack plane and shown in Fig. 8; the stress intensity factors along the crack front, corresponding to each stage (excluding the last one), are plotted in Fig. 9, from which similar behaviour to the first example can be seen. The embedded elliptical crack problem of the first example is considered again as an example, but the crack surface makes an angle 4 = 45” with the y-axis, as shown in Fig. 4. This is a mixed problem and the fatigue crack growth is considered instead of quasi-static one. The Paris model is used to model fatigue crack propagation, which states that da = c(AKeff)“,
(12) emin is chosen as zero and a,,,,, constant in this example. After each DBEM analysis, the directions of the next increment of the crack growth are evaluated by minimizing Sih’s strain energy density factor. The maximum amount Au, of the next increment is fixed (0.2u), corresponding to the crack front point where maximum AKF occurs. From eqn (9)
AK,,=K,“,-K$=‘=KF(l-R),
Aa x c(AK~)~ AN
(13)
z c(AKz)”
(14)
and Aa,
AN.
From eqns (13) and (14), the following relationship is obtained
where da/dN is the rate of change of crack length with respect to the loading cycles, c = 1.5463 x lo-l6 and m = 3.88 are material constants; (10)
Therefore, the incremental size at any front point can be evaluated by
in which X;, stands for the effective stress intensity factor, adopted here as
f& = WI + Ku, 1)’+ 2K:,.
(11)
V-
LY N
Since linear elasticity is considered, R in eqn (10) can be written as
(9)
dN
811
B
_r
-0.2
L
!-0.6 -
FW
/6
I 0 Original
p/
-0.8 0
1
0 lncrimint I
I
I
20
40
60
w
,
80
: ANGLE
I
100
I
120
ALONG
4
I
,
140
THE
160
I
180
CRACK
I
,
ZOO
220
,
240
I
260
,
280
,
3DO
,
320
FRONT
Fig. 11. Normalized stress intensity factors for inclined embedded elliptical crack.
,
2.40
360
878
Y. Mi and M. H. Aliabadi
A total of five steps of increment are mrformed. The crack surface after the last analysis, including five crack fronts, is shown in Fig. 10(a) and the deformed
crack mesh after the last analysis in Fig. 10(b). Model I, II and III stress intensity factors corresponding to each analysis are plotted in Fig. 11.
CONCLUSIONS
An automatic procedure for simulation of threedimensional crack growth using the dual boundary element method has been presented in this paper. The method has the advantage of eliminating interior remeshing as the crack growths take place, hence resulting in a substantial amount of computational savings as well as simple modelling. The test examples demonstrate that, through using the minimum strain energy criterion and a modified fatigue formula, cracks grow into a shape with constant stress intensity factors. REFERENCES I.
2.
3.
4.
5.
6. 7.
A. R. Ingraffea, Nodal grafting for crack propagation studies. &t. J. numer. Eigng 11, 1185-1187 (i977). R. Perucchio and A. R. Ineraffea. An integrated boundary element analysis system with interaczve computer graphics for three-dimensional linear-elastic fracture mechanics. Compur. Struct. 20, 157-171 (1985). S. Valliappan and V. Murt, Automatic remeshing technique in quasi-static and dynamic crack propagation. Proc. NCJMETA, ‘85 ConJ, pp. 107-116, Swansea (1985). M. S. Shepherd, N. A. S. Yehia, G. S. Burd and T. J. Weinder, Automatic crack propagation tracking. Comput. Struct. 20, 211-223 (1985). E. M. Remzi and W. S. Blackburn, Automatic crack propagation studies in T-junctions and cross bars. Engng Compur. 7, 116-124 (1990). M. H. Aliabadi and D. P. Rooke, Numerical Fracture Mechanics. Computational Mechanics Publications, Southampton and Kluwer Academic, Dordrecht (1991). T. A. Cruse and P. M. Besuner, Residual life prediction for surface cracks in complex structural details. J. Aircraff 12, 369-375 (1975). I
I
8. G. Blandford, A. R. Ingraffea and J. A. Liggett, Two dimensional stress intensity factor computation using the boundary element method. Inr. J. numer. Engng 17, 387-404 (1981). 9 A. R. Ingraffea, G. Blandford and J. A. Liggett, Automatic modelling of mixed-mode fatigue and quasi-static crack propagation using the boundary element method. Proc. 14th Nat. Symp. on Fracture Mech. ASTM STP 791, I 407-I 426 (1987). 10. P. A. Wawrzynek, L. F. Martha and A. R. Ingraffea, A computational environment for the simulation of fracture process in three-dimensions. In Annul. Numer. Exper. Aspects of Three-dimensional Fracture Process
(Edited by A. J. Rosakis et al.), ASME, AMD-91, pp. 321-327 (1988). 11. J. Weaver, Three-dimensional crack analysis. Znr. J. Solids Struct. 13. 321-330 (1977). 12. J. 0. Watson, Hermitian cubic anh singular elements for plain strain. In Development in Boundary Element Merhod. Vol. 4 (Edited bv P. K. Baneriiee and T. 0. Watsonj. Elsevie; Applieci Science, Barcing (1986). 13. E. Z. Polch, T. A. Cruse and C. J. Huang, Traction BIE solution for flat cracks. Comput. Mech. 2, 253-267 (1987). 14. L. J. Gray,
L. F. Martha and A. R. Integraffea, Hypersingular integrals in boundary element fracture analysis. Znt. J. numer. Merh. Engng 29, 1135-1158, (1990). 15. A. Portela, M. H. Aliabadi and D. P. Rooke, The dual boundary element method: effective implementation for crack problems. Znr. J. numer. Meth. Engng 33, 1269-1287
(1992).
16. Y. Mi and M. H. Aliabadi, Dual boundary element method for three-dimensional fracture mechanics analysis. Engng Anal. Boundary Elements 10, 161-171 (1992). 17. A. Portela, M. H. Aliabadi and D. P. Rooke, Dual boundary element incremental analysis of crack propagation. Comput. Struct. 46, 237-247 (1993). 18. P. Fedelinski, M. H. Aliabadi and D. P. Rooke, The dual boundary element method in dynamic fracture mechanics. Engng. Anal. Boundary Elements 12, 203-210 (1993).
19. V. Leitlo, M. H. Aliabadi and D. P. Rooke, The dual boundary element formulation for elastoplastic fracture mechanics, in press. 20. G. C. Sih, Mechanics of Fracture Initiation and Propagation. Kluwer Academic, Dordrecht (1991). 21. T. A. Cruse, Mathematical formulation of the boundary integral equation methods in solid mechanics. Report No. AFOSR-TR-77-1002, Pratt and Whitney Aircraft Group (1977).