- Email: [email protected]
A simulation of fatigue crack propagation in a welded T-joint using 3D boundary element method
International Journal of Pressure Vessels and Piping 80 (2003) 111–120 www.elsevier.com/locate/ijpvp
A simulation of fatigue crack propagation in a welded T-joint using 3D boundary element method Zhihai Xianga, Seng Tjhen Lieb, Bo Wanga, Zhangzhi Cena,* b
a Department of Engineering Mechanics, Tsinghua University, Beijing 100084, People’s Republic of China School of Civil and Structural Engineering, Nanyang Technological University, Nanyang Avenue, Singapore, Singapore 639798
Received 13 December 2000; revised 21 November 2002; accepted 17 January 2003
Abstract A general procedure to investigate the fatigue propagation process of a 3D surface crack based on multi-region Boundary Element Method is detailed in this paper. The mesh can be automatically regenerated as the crack propagates. A new formula for estimating the effective stress intensity factor is used to calculate the crack extension. The maximum principal stress criterion is then employed to predict the crack growth direction. Comparison between numerical and experimental results of a welded T-joint shows that the proposed procedure is reliable. q 2003 Elsevier Science Ltd. All rights reserved. Keywords: 3D fatigue surface crack; Crack propagation; T-joint; Boundary element method
1. Introduction The simulation of fatigue crack propagation in flawed components plays an important role in structural safety assessment. Most research work [1 – 3] has been carried out using the finite element method. However, this approach is limited, in practice, to some simple problems due to the difficulty of mesh generation. With the important features of high calculation accuracy and the reduction of the model dimension, boundary element method (BEM) requires only a relatively coarse surface mesh. This substantially speeds up the initial data preparation time and makes the re-meshing strategy much easier than other numerical methods. Furthermore, with its ability to model high stress gradients accurately and efficiently, BEM has naturally become a popular approach in this field [4]. Ingraffea, Blandford and Ligget [5] first attempted to automatically simulate mixed-mode crack growth for 2D problems with multi-region BEM. Then Grestle [6] extended this method to 3D problems. Mi and Aliabadi * Corresponding author. Tel.: þ 86-10-627-85454; fax: 86-10-62770349. E-mail address: [email protected] (Z. Cen).
[7 –9], Cisilino and Aliabadi [10], Wilde and Aliabadi [11] presented the application of the Dual Boundary Element Method (DBEM) to simulate the propagation of 3D mixed-mode cracks. Irrespective of whether the multi-region BEM or DBEM is used, two essential factors should be carefully considered for the successful simulation of 3D crack growth. The first factor relates to the criteria for the rate and direction of crack growth. For crack growth rate, Paris’ Law [12] must be extended to the mixed-mode case, usually by replacing the Mode I stress intensity factor range DKI with an effective stress intensity factor range DKeff : For crack growth direction, the maximum principal stress criterion [6, 14] and the minimum strain energy density criterion [7,8,15] are commonly used. However, it was reported [1] that the difference between these two criteria is negligible. The second factor is the generation of the initial mesh and the re-meshing strategy. These will be thoroughly discussed in this paper. Another important aspect in the crack development procedure is the treatment of the corner regions near the intersection points between the crack front and the free surface [16]. This was ignored in most of the previous work. In this paper, it is treated by an extrapolation method.
0308-0161/03/$ - see front matter q 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0308-0161(03)00022-X
112
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
Finally, the strategies proposed in this paper are verified by a simulation of the fatigue crack growth in a welded T-joint. The performances of three fracture criteria are illustrated by comparison between the numerical results and the experimental data.
2. Fatigue crack growth criteria The Paris’ Law [12] is expressed as da ¼ cðDKI Þm dN
ð1Þ
where a is the crack length, N is the number of fatigue loading cycles. The two constants c and m are material parameters. DKI is the Mode I stress intensity factor range. For extension to the mixed-mode case, it is replaced by DKeff : To the authors’ knowledge, the following formulae were proposed for evaluating DKeff in the 3D case. Tanaka [13] derived: !1=4 4 8DKIII 4 4 ð2Þ DKeff ¼ DKI þ 8DKII þ 12n
Fig. 1. Local coordinate system at the crack tip.
When using Eq. (2), DKIeff in Eq. (5) is taken as DKeff ¼ DKI ; when using Eq. (3), DKIeff ¼ DKI þ lDKIII l; with Eq. (4),
KIeff
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 KIII ¼ KI2 þ 12n
Writing DKIeff ¼ DKI þ blDKIII l; Grestle [6] proposed DKeff ¼ ½ðDKI þ blDKIII lÞ2 þ 2DKII2 1=2
ð3Þ
where b is an empirical factor and is taken as 1 in this paper. It is noted [1,22] that as the crack propagates, the percentage contributions of Mode II and Mode III decrease rapidly compared to Mode I, and the total energy to drive the crack growth is influenced by the ratio of DKII and DKIII to DKI : Accordingly, Lie et al. [21] proposed that ( DKeff ¼ aDKI ð4Þ a ¼ 1 þ A½eBr ðr 2 CÞ þ C
and DKIeff ¼ Max KIeff 2 Min KIeff : The above three fracture criteria are summarized in Table 1.
3. Mesh generation and re-meshing In this paper, the multi-region BEM is used to simulate crack propagation in a welded T-joint specimen [17]. As Fig. 2(a) shows, the make-up of the specimen includes two attachments fillet-welded onto a plate.
where r¼
Table 1 List of fracture criteria
2 DKII2 þ DKIII =ð1 2 nÞ ; 2 DKI
A ¼ 676:51; B ¼ 21010:91 and C ¼ 3:31 £ 1024 : For simplicity and without loss of accuracy, the maximum principal stress criterion is used to determine the crack growth direction in this paper. Following this criterion, the crack growth direction from one crack front point is predicted in the local coordinate system (Fig. 1) by 2 3 s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u0 1 4 DKIeff DKIeff 2 tan ¼ 2 þ8 5 4 DKII 2 DKII ¼ DKIeff
22DKII pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ ðDKIeff Þ2 þ 8ðDKII Þ2
DKIeff
Criterion
DKeff
1
DKeff ¼½DKI4 þ 8DKII4 þ
2
2 nÞ
DKIeff ¼ DKI 1=4
DKeff ¼½ðDKI þ lDKIII lÞ2 þ
3
4 8DKIII =ð1
(
DKIeff ¼ DKI þ lDKIII l
2DKII2 1=2
DKeff ¼ aDKI
a ¼ 1 þ AðeBr ðr 2 CÞ þ CÞ
KIeff
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 KIII ¼ KI2 þ 12n
, DKIeff ¼MaxðKIeff Þ
ð5Þ
2 MinðKIeff Þ
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
113
Fig. 2. Dimensions and boundary element mesh of the T-joint specimen: (a) specimen dimensions; (b) boundary element mesh on the crack face; (c) details of the weld toe.
114
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
Fig. 3. Fatigue crack shape development.
3.1. Initial mesh generation The whole structure is divided into triangular and quadrilateral isoparametric elements with 4– 8 variablenumber nodes. The asymptotic behaviours of the displacements near the crack front can be correctly obtained using quarter-point elements [9,18 – 21] and transition elements [18 –21]. From these displacements, the mode I, II and III stress intensity factors KI ; KII and KIII can be evaluated accordingly [7 –9]. A difficult task is the mesh generation for the surface crack, especially when the surface crack is in mixed mode and the crack front is not regular in shape. A typical BEM mesh is shown in Fig. 2(b) for one side of the crack face, which is slightly twisted and the crack front is pseudo-semielliptical. Following Ref. [21], the assumptions and steps are taken as follows:
Step 1 Assume the shape of the two crack corners are two elliptical surfaces. Step 2 Use a linear interpolation method to obtain the intersection line L (Fig. 2(b) and (c)) between the crack face and the free surface, according to experimental points on that line. Step 3 Calculate the intersection points of the line L and the two elliptical surfaces to estimate the positions of the two crack corner points. Step 4 Use the spline to interpolate the middle part of the crack front with experimental points. Step 5 Construct the crack face with experimental points and calculated interpolation points along the crack front and the line L; by using ordinary eight-node boundary elements. Step 6 Generate quarter-point elements [9,18 –21] along the two sides of the crack front. Local coordinates
Fig. 4. Crack cross-section profiles.
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
of corner nodes of each element along the crack front are calculated at the same time, except for the two nodes located at the place where the crack front intersects the free surface. Then generate transition elements [18 – 21] just behind the two rows of quarter-point elements. Step 7 Generate the rest of the elements. Note that the z-coordinate of the nodes on the ligament is
115
determined by the average value of the z-coordinates of the nodes on the crack face. This makes the artificial boundary between two sub-regions unique. In a typical mesh for this specimen (Fig. 2), there are 2214 and 1974 degrees of freedom for sub-region I and sub-region II, respectively.
Fig. 5. Comparison of crack cross-section profiles: (a) at probe 2; (b) at probe 3; (c) at probe 4; (d) at probe 5; (e) at probe 6; (f) at probe 7.
116
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
Fig. 5 (continued )
3.2. Re-meshing strategy To simulate crack propagation, the following steps are taken: Step 1 Use the multi-region BEM program to obtain the displacements of the specimen. Calculate the stress
intensity factors [9,18 – 21] and DKeff for corner nodes of each element along the crack front, except for the two nodes located at the place where the crack front intersects the free surface. Step 2 Limit the maximum crack extension Da to a certain value (0.4 mm is adopted in the current work), and calculate the increment of fatigue loading cycles
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
117
Table 2 Numerical and experimental results of crack depth Probe no.
Cycles
Measured crack depth (mm)
Numerical results Criterion 1
Criterion 2
Criterion 3
Depth (mm)
Error (%)
Depth (mm)
Error (%)
Depth (mm)
Error (%)
2
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
1.3 1.4 1.5 1.7 1.9 2.1 2.3 2.5 2.8
1.3 1.3 1.4 1.5 1.6 1.7 1.7 1.8 1.8
0.0 21.9 26.4 28.9 213.8 219.3 225.2 228.5 235.6
1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1
0.0 21.3 24.1 24.8 28.0 212.6 217.7 219.7 226.4
1.3 1.4 1.6 1.8 2.0 2.3 2.5 2.7 3.0
0.0 2.8 3.4 7.2 9.3 8.6 7.0 9.3 5.3
3
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
4.1 4.3 4.5 4.7 4.9 5.1 5.3 5.5 5.8
4.1 4.3 4.4 4.6 4.7 4.8 4.9 5.1 5.2
0.0 21.1 22.4 22.1 24.3 25.5 26.6 28.2 210.5
4.1 4.3 4.4 4.6 4.8 4.9 5.0 5.2 5.3
0.0 20.8 21.8 21.1 23.1 24.1 25.0 26.5 28.7
4.1 4.4 4.7 5.1 5.5 5.8 6.2 6.6 7.0
0.0 2.0 4.6 9.5 12.0 14.3 17.0 19.2 20.6
4
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
5.7 6.0 6.2 6.3 6.6 6.9 6.9 7.2 7.4
5.7 5.9 6.0 6.2 6.3 6.4 6.6 6.7 6.8
0.0 22.7 22.5 22.5 24.5 26.1 25.5 26.6 27.6
5.7 5.9 6.1 6.2 6.4 6.6 6.7 6.8 7.0
0.0 22.3 21.6 21.3 22.9 24.3 23.5 24.5 25.3
5.7 6.0 6.3 6.6 7.0 7.3 7.6 7.9 8.3
0.0 20.8 2.0 4.6 5.4 6.0 9.3 10.8 12.9
5
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
6.1 6.4 6.6 6.8 6.9 7.2 7.3 7.6 7.8
6.1 6.4 6.4 6.6 6.7 6.9 7.0 7.1 7.2
0.0 21.5 22.0 23.9 22.5 24.3 24.3 26.4 27.0
6.1 6.3 6.4 6.6 6.8 6.9 7.0 7.1 7.3
0.0 21.5 21.9 23.7 22.1 23.8 23.7 25.8 26.4
6.1 6.4 6.8 7.2 7.7 8.0 8.3 8.6 8.9
0.0 1.1 3.5 5.8 11.5 11.6 13.8 13.4 15.0
6
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
5.1 5.3 5.5 5.8 6.1 6.3 6.5 6.7 7.0
5.1 5.2 5.3 5.9 5.6 5.7 5.9 6.0 6.1
0.0 21.3 23.6 24.9 27.6 29.0 29.7 211.2 212.3
5.1 5.2 5.4 5.5 5.7 5.8 5.9 6.0 6.2
0.0 21.1 23.3 24.5 27.0 28.3 28.9 210.5 211.5
5.1 5.3 5.6 5.9 6.2 6.6 6.9 7.3 7.7
0.0 0.9 1.0 2.1 2.4 3.8 6.5 8.1 10.3
7
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
2.3 2.5 2.9 3.1 3.4 3.7 3.9 4.3 4.3
2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.2
0.0 24.4 214.4 216.3 218.9 222.3 223.9 228.0 227.4
2.3 2.4 2.5 2.7 2.8 2.9 3.0 3.1 3.2
0.0 23.7 213.4 214.9 217.3 220.6 222.1 226.2 225.6
2.3 2.5 2.8 3.1 3.5 3.8 4.0 4.4 4.7
0.0 0.9 23.3 0.7 3.3 2.8 4.7 3.2 8.6
118
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
DN from Da ¼ cðDKeff Þm DN
Step 3
Step 4 Step 5
Step 6
Step 7
Step 8
Step 9
ð6Þ
for each corner node along the crack front, and determine the minimum DN: Then recalculate Da for each corner node from Eq. (6) by using the minimum DN: Calculate the crack growth direction u0 in the local coordinate system for each corner node along the crack front by using Eq. (5). Extend each corner node along the crack front to new position according to Da and u0 : Assume the shape of the two crack corners as two elliptical surfaces. Calculate the intersection points of the line L and the two elliptical surfaces to estimate the new positions of the two crack corner points. Use the spline to interpolate the middle part of the crack front with new coordinates of the corner nodes. Reconstruct the crack face according to the new crack front points and points sampled from the old crack face, using the same eight-node boundary elements. Without changing the topology of the old multiregion BEM mesh, regenerate the new mesh and recalculate the local coordinates of each corner node, except for the two nodes located at the place where the crack front intersects the free surface. Go to Step 1 for the next crack increment, or stop the simulation when all the incremental steps are completed.
the actual crack cross-section profiles (Fig. 4). Thus, the actual 3D shape of the crack surface is detailed. In this test, the plastic zone is very small, and there is no crack closure because the applied stess ratio R ¼ 0: Therefore, an elastic algorithm is sufficient to simulate the crack growth. Following the strategies described in Section 3, the crack propagation is simulated using the fracture criteria one to three (Table 1). Because the tolerance of the test measurement is 0.1 mm, we do not have high degree of confidence on the data measured during
4. Numerical results The T-joint, simulated in this paper, was made of structural steel BS4360 grade 43, for which c and m in Paris’ Law are equal to 1.834 £ 10213 and 3.0, pffiffiffirespectively; for Da=DN in mm/cycle and DK in MPa m: In the test, the specimen was subjected to an axial sinusoidal constant amplitude (200 kN) loading of 7 Hz with stress ratio R ¼ 0 until failure. Owing to the severe stress concentration, fatigue crack growth always initiated at the weld toe. An Alternating Current Potential Drop (ACPD) technique was used to monitor the development of the 3D fatigue crack in the expected crack location. By measuring the depth at fixed points along the crack front, and by choosing suitable inspection intervals, it was possible to follow the crack growth from initiation to through wall penetration. The fatigue crack development shape was recorded at every 2000 cycles, and this was obtained by using eight probes spaced equally at 10 mm apart as shown in Fig. 2(c) and Fig. 3. Finally, the plate was sliced through the thickness at each of the probe positions to get
Fig. 6. Re-meshing process of the surface crack.
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
119
Crack corner nodes are indirectly estimated to circumvent the treatment of the different singularity of the crack tip. For the simulation of a welded T-joint, numerical results from three fatigue crack growth criteria are compared. The calculated crack depths generally agree with the experimental ones. The maximum principal stress criterion works quite well in the simulation. As the crack propagates, the crack front tends to be an ellipse, and Mode I becomes dominant compared to Mode II and III.
Acknowledgements This project is partly supported by National Nature Science Foundation of China, Failure Mechanics Laboratory in Tsinghua University (Open Laboratory in Chinese Ministry of Education), and Basic Research Foundation in Tsinghua University. The experimental test is supported by the Applied Research Grant under No. RG32/95 at the Nanyang Technological University in Singapore. Fig. 7. Crack front propagation of the surface crack.
References the initial phase of the crack propagation. Therefore, the simulation starts from cycle 2,625,734 when the crack propagation is stable. The comparison between the numerical results and the experimental data is shown in Fig. 5 and Table 2. Fig. 5 reveals that no matter which criterion is used, the numerical simulation can approximately follow the crack curvature. Criterion three behaves better at probe 2 (Fig. 5(a)). For the crack growth rate, criterion three gives a more conservative evaluation than the other two criteria. This can also be observed from Table 2, which shows that after eight iteration steps, the maximum errors of crack length calculated by criterion one to criterion three are 2 35.6, 2 26.4 and 20.6%, respectively. The re-meshing process and the crack propagation results are shown in Figs. 6 and 7. From these results, it is observed that as the crack propagates, the crack front tends to become an ellipse, and Mode I becomes more dominant than Modes II and III.
5. Conclusions Based on the multi-region BEM, a general procedure to simulate the propagation of a 3D mixed-mode fatigue surface crack is presented. Steps for generating the initial BEM mesh and the automatic re-meshing strategy are detailed. This automatic re-meshing algorithm retains the topology of the original mesh. Thus, the computing time for each crack increment step is constant.
[1] Bowness D, Lee MMK. Fatigue crack curvature under the weld toe in an offshore tubular joint. Int J Fatigue 1998;20:481–90. [2] Dhondt D. Automatic 3-D Mode I crack propagation calculations with finite elements. Int J Numer Meth Engng 1998;41:739–57. [3] Newman JC, Harris CE, Piascik RS, Dawicke DS. Methodology for predicting the onset of widespread fatigue damage in lapsplice joints and structure. Proceedings of the Seventh International Fatigue Congress Beijing, P.R. China; June 1999. p. 81–92. [4] Aliabadi MH. Boundary element formulations in fracture mechanics. Appl Mech Rev 1997;50:83 –96. [5] Ingraffea AR, Blandford G, Liggett JA. Automatic modeling of mixed-mode fatigue and quasi-static crack propagation using the boundary element method. Fourteenth National Symposium on Fracture, ASTM STP 1987;791:1407–26. [6] Grestle WH. Finite and boundary element modeling of crack propagation in two- and three-dimensions using interactive computer graphics. PhD Thesis. Cornell University, Ithaca, NY; 1986. [7] Mi Y, Aliabadi MH. Three-dimensional crack growth simulation using BEM. Int J Comput Struct 1994;52:871–8. [8] Mi Y, Aliabadi MH. Automatic procedure for mixed-mode crack growth analysis. Commun Numer Meth Engng 1995;11:167–77. [9] Mi Y. Three-dimensional analysis of crack growth. Topics in Engineering, vol. 28. Southampton, UK: Computational Mechanics Publications; 1997. [10] Cisilino AP, Aliabadi MH. Three-dimensional BEM analysis for fatigue crack growth in welded components. Int J Pressure Vessels Piping 1997;70:135– 44. [11] Wilde AJ, Aliabadi MH. A 3-D dual BEM formulation for the analysis of crack growth. Comput Mech 1999;23:250–7. [12] Paris PC, Gomez RE, Anderson WE. A rational analytic theory of fatigue. Trend Engng 1961;13:9–14. [13] Tanaka K. Fatigue crack propagation from a crack inclined to the cyclic tensile axis. Engng Fracture Mech 1974;6:493–507.
120
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
[14] Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J Basic Engng Transition ASME 1963; 85:519–27. [15] Sih GC, Cha BCK. A fracture criterion for three-dimensional crack problems. Engng Fracture Mech 1974;6:699–732. [16] Zdenak P, Bazant A, Estenssoro LF. Surface singularity and crack propagation. Int J Solids Struct 1979;15:405–26. [17] Yeo WL. Monitoring the crack growth using the ACPD equipments— effect of transverse attachment. BSc Thesis. Nanyang Technological University, Singapore; 1998. [18] Zhu H. On three-dimensional crack analysis by boundary element method and some engineering applications. MEngng Thesis. Tsinghua University; 1987.
[19] Lie ST, Li G, Cen Z. Modelling tubular joints using a combination of finite and boundary elements. Theory and application of boundary element methods. Proceedings of the Eight China–Japan Symposium on BEM, Beijing, P.R. China; May 1998. p 257–66. [20] Lie ST, Li G, Cen Z. Analysis of cracked tubular joints using coupled finite and boundary elements. Engng Struct 2000;22: 272 –83. [21] Lie ST, Xiang Z, Wang B, Cen Z. Experimental and numerical simulation of 3D fatigue crack for plate-to-plate welded joints. Int J Fatigue 2000;22:411–24. [22] Pook LP. Fatigue crack paths. Proceedings of the Seventh International Fatigue Congress, Beijing, P.R. China; June 1999. p. 407–14.
A simulation of fatigue crack propagation in a welded T-joint using 3D boundary element method Zhihai Xianga, Seng Tjhen Lieb, Bo Wanga, Zhangzhi Cena,* b
a Department of Engineering Mechanics, Tsinghua University, Beijing 100084, People’s Republic of China School of Civil and Structural Engineering, Nanyang Technological University, Nanyang Avenue, Singapore, Singapore 639798
Received 13 December 2000; revised 21 November 2002; accepted 17 January 2003
Abstract A general procedure to investigate the fatigue propagation process of a 3D surface crack based on multi-region Boundary Element Method is detailed in this paper. The mesh can be automatically regenerated as the crack propagates. A new formula for estimating the effective stress intensity factor is used to calculate the crack extension. The maximum principal stress criterion is then employed to predict the crack growth direction. Comparison between numerical and experimental results of a welded T-joint shows that the proposed procedure is reliable. q 2003 Elsevier Science Ltd. All rights reserved. Keywords: 3D fatigue surface crack; Crack propagation; T-joint; Boundary element method
1. Introduction The simulation of fatigue crack propagation in flawed components plays an important role in structural safety assessment. Most research work [1 – 3] has been carried out using the finite element method. However, this approach is limited, in practice, to some simple problems due to the difficulty of mesh generation. With the important features of high calculation accuracy and the reduction of the model dimension, boundary element method (BEM) requires only a relatively coarse surface mesh. This substantially speeds up the initial data preparation time and makes the re-meshing strategy much easier than other numerical methods. Furthermore, with its ability to model high stress gradients accurately and efficiently, BEM has naturally become a popular approach in this field [4]. Ingraffea, Blandford and Ligget [5] first attempted to automatically simulate mixed-mode crack growth for 2D problems with multi-region BEM. Then Grestle [6] extended this method to 3D problems. Mi and Aliabadi * Corresponding author. Tel.: þ 86-10-627-85454; fax: 86-10-62770349. E-mail address: [email protected] (Z. Cen).
[7 –9], Cisilino and Aliabadi [10], Wilde and Aliabadi [11] presented the application of the Dual Boundary Element Method (DBEM) to simulate the propagation of 3D mixed-mode cracks. Irrespective of whether the multi-region BEM or DBEM is used, two essential factors should be carefully considered for the successful simulation of 3D crack growth. The first factor relates to the criteria for the rate and direction of crack growth. For crack growth rate, Paris’ Law [12] must be extended to the mixed-mode case, usually by replacing the Mode I stress intensity factor range DKI with an effective stress intensity factor range DKeff : For crack growth direction, the maximum principal stress criterion [6, 14] and the minimum strain energy density criterion [7,8,15] are commonly used. However, it was reported [1] that the difference between these two criteria is negligible. The second factor is the generation of the initial mesh and the re-meshing strategy. These will be thoroughly discussed in this paper. Another important aspect in the crack development procedure is the treatment of the corner regions near the intersection points between the crack front and the free surface [16]. This was ignored in most of the previous work. In this paper, it is treated by an extrapolation method.
0308-0161/03/$ - see front matter q 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0308-0161(03)00022-X
112
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
Finally, the strategies proposed in this paper are verified by a simulation of the fatigue crack growth in a welded T-joint. The performances of three fracture criteria are illustrated by comparison between the numerical results and the experimental data.
2. Fatigue crack growth criteria The Paris’ Law [12] is expressed as da ¼ cðDKI Þm dN
ð1Þ
where a is the crack length, N is the number of fatigue loading cycles. The two constants c and m are material parameters. DKI is the Mode I stress intensity factor range. For extension to the mixed-mode case, it is replaced by DKeff : To the authors’ knowledge, the following formulae were proposed for evaluating DKeff in the 3D case. Tanaka [13] derived: !1=4 4 8DKIII 4 4 ð2Þ DKeff ¼ DKI þ 8DKII þ 12n
Fig. 1. Local coordinate system at the crack tip.
When using Eq. (2), DKIeff in Eq. (5) is taken as DKeff ¼ DKI ; when using Eq. (3), DKIeff ¼ DKI þ lDKIII l; with Eq. (4),
KIeff
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 KIII ¼ KI2 þ 12n
Writing DKIeff ¼ DKI þ blDKIII l; Grestle [6] proposed DKeff ¼ ½ðDKI þ blDKIII lÞ2 þ 2DKII2 1=2
ð3Þ
where b is an empirical factor and is taken as 1 in this paper. It is noted [1,22] that as the crack propagates, the percentage contributions of Mode II and Mode III decrease rapidly compared to Mode I, and the total energy to drive the crack growth is influenced by the ratio of DKII and DKIII to DKI : Accordingly, Lie et al. [21] proposed that ( DKeff ¼ aDKI ð4Þ a ¼ 1 þ A½eBr ðr 2 CÞ þ C
and DKIeff ¼ Max KIeff 2 Min KIeff : The above three fracture criteria are summarized in Table 1.
3. Mesh generation and re-meshing In this paper, the multi-region BEM is used to simulate crack propagation in a welded T-joint specimen [17]. As Fig. 2(a) shows, the make-up of the specimen includes two attachments fillet-welded onto a plate.
where r¼
Table 1 List of fracture criteria
2 DKII2 þ DKIII =ð1 2 nÞ ; 2 DKI
A ¼ 676:51; B ¼ 21010:91 and C ¼ 3:31 £ 1024 : For simplicity and without loss of accuracy, the maximum principal stress criterion is used to determine the crack growth direction in this paper. Following this criterion, the crack growth direction from one crack front point is predicted in the local coordinate system (Fig. 1) by 2 3 s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u0 1 4 DKIeff DKIeff 2 tan ¼ 2 þ8 5 4 DKII 2 DKII ¼ DKIeff
22DKII pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ ðDKIeff Þ2 þ 8ðDKII Þ2
DKIeff
Criterion
DKeff
1
DKeff ¼½DKI4 þ 8DKII4 þ
2
2 nÞ
DKIeff ¼ DKI 1=4
DKeff ¼½ðDKI þ lDKIII lÞ2 þ
3
4 8DKIII =ð1
(
DKIeff ¼ DKI þ lDKIII l
2DKII2 1=2
DKeff ¼ aDKI
a ¼ 1 þ AðeBr ðr 2 CÞ þ CÞ
KIeff
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 KIII ¼ KI2 þ 12n
, DKIeff ¼MaxðKIeff Þ
ð5Þ
2 MinðKIeff Þ
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
113
Fig. 2. Dimensions and boundary element mesh of the T-joint specimen: (a) specimen dimensions; (b) boundary element mesh on the crack face; (c) details of the weld toe.
114
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
Fig. 3. Fatigue crack shape development.
3.1. Initial mesh generation The whole structure is divided into triangular and quadrilateral isoparametric elements with 4– 8 variablenumber nodes. The asymptotic behaviours of the displacements near the crack front can be correctly obtained using quarter-point elements [9,18 – 21] and transition elements [18 –21]. From these displacements, the mode I, II and III stress intensity factors KI ; KII and KIII can be evaluated accordingly [7 –9]. A difficult task is the mesh generation for the surface crack, especially when the surface crack is in mixed mode and the crack front is not regular in shape. A typical BEM mesh is shown in Fig. 2(b) for one side of the crack face, which is slightly twisted and the crack front is pseudo-semielliptical. Following Ref. [21], the assumptions and steps are taken as follows:
Step 1 Assume the shape of the two crack corners are two elliptical surfaces. Step 2 Use a linear interpolation method to obtain the intersection line L (Fig. 2(b) and (c)) between the crack face and the free surface, according to experimental points on that line. Step 3 Calculate the intersection points of the line L and the two elliptical surfaces to estimate the positions of the two crack corner points. Step 4 Use the spline to interpolate the middle part of the crack front with experimental points. Step 5 Construct the crack face with experimental points and calculated interpolation points along the crack front and the line L; by using ordinary eight-node boundary elements. Step 6 Generate quarter-point elements [9,18 –21] along the two sides of the crack front. Local coordinates
Fig. 4. Crack cross-section profiles.
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
of corner nodes of each element along the crack front are calculated at the same time, except for the two nodes located at the place where the crack front intersects the free surface. Then generate transition elements [18 – 21] just behind the two rows of quarter-point elements. Step 7 Generate the rest of the elements. Note that the z-coordinate of the nodes on the ligament is
115
determined by the average value of the z-coordinates of the nodes on the crack face. This makes the artificial boundary between two sub-regions unique. In a typical mesh for this specimen (Fig. 2), there are 2214 and 1974 degrees of freedom for sub-region I and sub-region II, respectively.
Fig. 5. Comparison of crack cross-section profiles: (a) at probe 2; (b) at probe 3; (c) at probe 4; (d) at probe 5; (e) at probe 6; (f) at probe 7.
116
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
Fig. 5 (continued )
3.2. Re-meshing strategy To simulate crack propagation, the following steps are taken: Step 1 Use the multi-region BEM program to obtain the displacements of the specimen. Calculate the stress
intensity factors [9,18 – 21] and DKeff for corner nodes of each element along the crack front, except for the two nodes located at the place where the crack front intersects the free surface. Step 2 Limit the maximum crack extension Da to a certain value (0.4 mm is adopted in the current work), and calculate the increment of fatigue loading cycles
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
117
Table 2 Numerical and experimental results of crack depth Probe no.
Cycles
Measured crack depth (mm)
Numerical results Criterion 1
Criterion 2
Criterion 3
Depth (mm)
Error (%)
Depth (mm)
Error (%)
Depth (mm)
Error (%)
2
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
1.3 1.4 1.5 1.7 1.9 2.1 2.3 2.5 2.8
1.3 1.3 1.4 1.5 1.6 1.7 1.7 1.8 1.8
0.0 21.9 26.4 28.9 213.8 219.3 225.2 228.5 235.6
1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1
0.0 21.3 24.1 24.8 28.0 212.6 217.7 219.7 226.4
1.3 1.4 1.6 1.8 2.0 2.3 2.5 2.7 3.0
0.0 2.8 3.4 7.2 9.3 8.6 7.0 9.3 5.3
3
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
4.1 4.3 4.5 4.7 4.9 5.1 5.3 5.5 5.8
4.1 4.3 4.4 4.6 4.7 4.8 4.9 5.1 5.2
0.0 21.1 22.4 22.1 24.3 25.5 26.6 28.2 210.5
4.1 4.3 4.4 4.6 4.8 4.9 5.0 5.2 5.3
0.0 20.8 21.8 21.1 23.1 24.1 25.0 26.5 28.7
4.1 4.4 4.7 5.1 5.5 5.8 6.2 6.6 7.0
0.0 2.0 4.6 9.5 12.0 14.3 17.0 19.2 20.6
4
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
5.7 6.0 6.2 6.3 6.6 6.9 6.9 7.2 7.4
5.7 5.9 6.0 6.2 6.3 6.4 6.6 6.7 6.8
0.0 22.7 22.5 22.5 24.5 26.1 25.5 26.6 27.6
5.7 5.9 6.1 6.2 6.4 6.6 6.7 6.8 7.0
0.0 22.3 21.6 21.3 22.9 24.3 23.5 24.5 25.3
5.7 6.0 6.3 6.6 7.0 7.3 7.6 7.9 8.3
0.0 20.8 2.0 4.6 5.4 6.0 9.3 10.8 12.9
5
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
6.1 6.4 6.6 6.8 6.9 7.2 7.3 7.6 7.8
6.1 6.4 6.4 6.6 6.7 6.9 7.0 7.1 7.2
0.0 21.5 22.0 23.9 22.5 24.3 24.3 26.4 27.0
6.1 6.3 6.4 6.6 6.8 6.9 7.0 7.1 7.3
0.0 21.5 21.9 23.7 22.1 23.8 23.7 25.8 26.4
6.1 6.4 6.8 7.2 7.7 8.0 8.3 8.6 8.9
0.0 1.1 3.5 5.8 11.5 11.6 13.8 13.4 15.0
6
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
5.1 5.3 5.5 5.8 6.1 6.3 6.5 6.7 7.0
5.1 5.2 5.3 5.9 5.6 5.7 5.9 6.0 6.1
0.0 21.3 23.6 24.9 27.6 29.0 29.7 211.2 212.3
5.1 5.2 5.4 5.5 5.7 5.8 5.9 6.0 6.2
0.0 21.1 23.3 24.5 27.0 28.3 28.9 210.5 211.5
5.1 5.3 5.6 5.9 6.2 6.6 6.9 7.3 7.7
0.0 0.9 1.0 2.1 2.4 3.8 6.5 8.1 10.3
7
2,625,734 2,631,720 2,637,726 2,643,732 2,649,738 2,653,742 2,657,746 2,661,743 2,665,747
2.3 2.5 2.9 3.1 3.4 3.7 3.9 4.3 4.3
2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.2
0.0 24.4 214.4 216.3 218.9 222.3 223.9 228.0 227.4
2.3 2.4 2.5 2.7 2.8 2.9 3.0 3.1 3.2
0.0 23.7 213.4 214.9 217.3 220.6 222.1 226.2 225.6
2.3 2.5 2.8 3.1 3.5 3.8 4.0 4.4 4.7
0.0 0.9 23.3 0.7 3.3 2.8 4.7 3.2 8.6
118
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
DN from Da ¼ cðDKeff Þm DN
Step 3
Step 4 Step 5
Step 6
Step 7
Step 8
Step 9
ð6Þ
for each corner node along the crack front, and determine the minimum DN: Then recalculate Da for each corner node from Eq. (6) by using the minimum DN: Calculate the crack growth direction u0 in the local coordinate system for each corner node along the crack front by using Eq. (5). Extend each corner node along the crack front to new position according to Da and u0 : Assume the shape of the two crack corners as two elliptical surfaces. Calculate the intersection points of the line L and the two elliptical surfaces to estimate the new positions of the two crack corner points. Use the spline to interpolate the middle part of the crack front with new coordinates of the corner nodes. Reconstruct the crack face according to the new crack front points and points sampled from the old crack face, using the same eight-node boundary elements. Without changing the topology of the old multiregion BEM mesh, regenerate the new mesh and recalculate the local coordinates of each corner node, except for the two nodes located at the place where the crack front intersects the free surface. Go to Step 1 for the next crack increment, or stop the simulation when all the incremental steps are completed.
the actual crack cross-section profiles (Fig. 4). Thus, the actual 3D shape of the crack surface is detailed. In this test, the plastic zone is very small, and there is no crack closure because the applied stess ratio R ¼ 0: Therefore, an elastic algorithm is sufficient to simulate the crack growth. Following the strategies described in Section 3, the crack propagation is simulated using the fracture criteria one to three (Table 1). Because the tolerance of the test measurement is 0.1 mm, we do not have high degree of confidence on the data measured during
4. Numerical results The T-joint, simulated in this paper, was made of structural steel BS4360 grade 43, for which c and m in Paris’ Law are equal to 1.834 £ 10213 and 3.0, pffiffiffirespectively; for Da=DN in mm/cycle and DK in MPa m: In the test, the specimen was subjected to an axial sinusoidal constant amplitude (200 kN) loading of 7 Hz with stress ratio R ¼ 0 until failure. Owing to the severe stress concentration, fatigue crack growth always initiated at the weld toe. An Alternating Current Potential Drop (ACPD) technique was used to monitor the development of the 3D fatigue crack in the expected crack location. By measuring the depth at fixed points along the crack front, and by choosing suitable inspection intervals, it was possible to follow the crack growth from initiation to through wall penetration. The fatigue crack development shape was recorded at every 2000 cycles, and this was obtained by using eight probes spaced equally at 10 mm apart as shown in Fig. 2(c) and Fig. 3. Finally, the plate was sliced through the thickness at each of the probe positions to get
Fig. 6. Re-meshing process of the surface crack.
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
119
Crack corner nodes are indirectly estimated to circumvent the treatment of the different singularity of the crack tip. For the simulation of a welded T-joint, numerical results from three fatigue crack growth criteria are compared. The calculated crack depths generally agree with the experimental ones. The maximum principal stress criterion works quite well in the simulation. As the crack propagates, the crack front tends to be an ellipse, and Mode I becomes dominant compared to Mode II and III.
Acknowledgements This project is partly supported by National Nature Science Foundation of China, Failure Mechanics Laboratory in Tsinghua University (Open Laboratory in Chinese Ministry of Education), and Basic Research Foundation in Tsinghua University. The experimental test is supported by the Applied Research Grant under No. RG32/95 at the Nanyang Technological University in Singapore. Fig. 7. Crack front propagation of the surface crack.
References the initial phase of the crack propagation. Therefore, the simulation starts from cycle 2,625,734 when the crack propagation is stable. The comparison between the numerical results and the experimental data is shown in Fig. 5 and Table 2. Fig. 5 reveals that no matter which criterion is used, the numerical simulation can approximately follow the crack curvature. Criterion three behaves better at probe 2 (Fig. 5(a)). For the crack growth rate, criterion three gives a more conservative evaluation than the other two criteria. This can also be observed from Table 2, which shows that after eight iteration steps, the maximum errors of crack length calculated by criterion one to criterion three are 2 35.6, 2 26.4 and 20.6%, respectively. The re-meshing process and the crack propagation results are shown in Figs. 6 and 7. From these results, it is observed that as the crack propagates, the crack front tends to become an ellipse, and Mode I becomes more dominant than Modes II and III.
5. Conclusions Based on the multi-region BEM, a general procedure to simulate the propagation of a 3D mixed-mode fatigue surface crack is presented. Steps for generating the initial BEM mesh and the automatic re-meshing strategy are detailed. This automatic re-meshing algorithm retains the topology of the original mesh. Thus, the computing time for each crack increment step is constant.
[1] Bowness D, Lee MMK. Fatigue crack curvature under the weld toe in an offshore tubular joint. Int J Fatigue 1998;20:481–90. [2] Dhondt D. Automatic 3-D Mode I crack propagation calculations with finite elements. Int J Numer Meth Engng 1998;41:739–57. [3] Newman JC, Harris CE, Piascik RS, Dawicke DS. Methodology for predicting the onset of widespread fatigue damage in lapsplice joints and structure. Proceedings of the Seventh International Fatigue Congress Beijing, P.R. China; June 1999. p. 81–92. [4] Aliabadi MH. Boundary element formulations in fracture mechanics. Appl Mech Rev 1997;50:83 –96. [5] Ingraffea AR, Blandford G, Liggett JA. Automatic modeling of mixed-mode fatigue and quasi-static crack propagation using the boundary element method. Fourteenth National Symposium on Fracture, ASTM STP 1987;791:1407–26. [6] Grestle WH. Finite and boundary element modeling of crack propagation in two- and three-dimensions using interactive computer graphics. PhD Thesis. Cornell University, Ithaca, NY; 1986. [7] Mi Y, Aliabadi MH. Three-dimensional crack growth simulation using BEM. Int J Comput Struct 1994;52:871–8. [8] Mi Y, Aliabadi MH. Automatic procedure for mixed-mode crack growth analysis. Commun Numer Meth Engng 1995;11:167–77. [9] Mi Y. Three-dimensional analysis of crack growth. Topics in Engineering, vol. 28. Southampton, UK: Computational Mechanics Publications; 1997. [10] Cisilino AP, Aliabadi MH. Three-dimensional BEM analysis for fatigue crack growth in welded components. Int J Pressure Vessels Piping 1997;70:135– 44. [11] Wilde AJ, Aliabadi MH. A 3-D dual BEM formulation for the analysis of crack growth. Comput Mech 1999;23:250–7. [12] Paris PC, Gomez RE, Anderson WE. A rational analytic theory of fatigue. Trend Engng 1961;13:9–14. [13] Tanaka K. Fatigue crack propagation from a crack inclined to the cyclic tensile axis. Engng Fracture Mech 1974;6:493–507.
120
Z. Xiang et al. / International Journal of Pressure Vessels and Piping 80 (2003) 111–120
[14] Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J Basic Engng Transition ASME 1963; 85:519–27. [15] Sih GC, Cha BCK. A fracture criterion for three-dimensional crack problems. Engng Fracture Mech 1974;6:699–732. [16] Zdenak P, Bazant A, Estenssoro LF. Surface singularity and crack propagation. Int J Solids Struct 1979;15:405–26. [17] Yeo WL. Monitoring the crack growth using the ACPD equipments— effect of transverse attachment. BSc Thesis. Nanyang Technological University, Singapore; 1998. [18] Zhu H. On three-dimensional crack analysis by boundary element method and some engineering applications. MEngng Thesis. Tsinghua University; 1987.
[19] Lie ST, Li G, Cen Z. Modelling tubular joints using a combination of finite and boundary elements. Theory and application of boundary element methods. Proceedings of the Eight China–Japan Symposium on BEM, Beijing, P.R. China; May 1998. p 257–66. [20] Lie ST, Li G, Cen Z. Analysis of cracked tubular joints using coupled finite and boundary elements. Engng Struct 2000;22: 272 –83. [21] Lie ST, Xiang Z, Wang B, Cen Z. Experimental and numerical simulation of 3D fatigue crack for plate-to-plate welded joints. Int J Fatigue 2000;22:411–24. [22] Pook LP. Fatigue crack paths. Proceedings of the Seventh International Fatigue Congress, Beijing, P.R. China; June 1999. p. 407–14.