Calculation of transient strain energy release rates under impact loading based on the virtual crack closure technique

ARTICLE IN PRESS

International Journal of Impact Engineering 34 (2007) 1047–1060 www.elsevier.com/locate/ijimpeng

Calculation of transient strain energy release rates under impact loading based on the virtual crack closure technique De Xiea,, Sherrill B. Biggers, Jr.b a

Alpha STAR Corporation, 5199 East Pacific Coast Highway, Suite 410, Long Beach, CA 90804, USA b Department of Mechanical Engineering, Clemson University, Clemson, SC 29634, USA Received 21 August 2005; received in revised form 21 February 2006; accepted 25 February 2006 Available online 11 May 2006

Abstract This paper describes an interface element to calculate the strain energy release rates based on the virtual crack closure technique (VCCT) in conjunction with finite element analysis (FEA). A very stiff spring is placed between the node pair at the crack tip to calculate the nodal forces. Dummy nodes are introduced to extract information for displacement openings behind the crack tip and the virtual crack jump ahead of the crack tip. This interface element leads to a direct calculation of the strain energy release rate (both components GI and GII) within a finite element analysis without extra post-processing. Several examples of stationary cracks under impact loading were examined. Dynamic stress intensity factors were converted from the calculated transient strain energy release rate for comparison with the available solutions by the others from numerical and experimental methods. The accuracy of the element is validated by the excellent agreement with these solutions. No convergence difficulty has been encountered for all the cases studied. Neither special singular elements nor the collapsed element technique is used at the crack tip. Therefore, the fracture interface element for VCCT is shown to be simple, efficient and robust in analyzing crack response to the dynamic loading. This element has been implemented into commercial FEA software ABAQUSs with the user defined element (UEL) and should be very useful in performing fracture analysis at a structural level by engineers using ABAQUSs. r 2006 Elsevier Ltd. All rights reserved. Keywords: Virtual crack closure technique; Dynamic fracture; Strain energy release rate; Stress intensity factor; Interfacial element

1. Introduction In recent years, due to the growing demands for predicting the behavior of cracked structures during impact and crash events, the importance of fracture mechanics has increased considerably in engineering applications. To evaluate structural integrity under dynamic loading, it is fundamental to know the transient fracture parameters such as dynamic stress intensity factors and/or dynamic strain energy release rates. Some analytical solutions are available but are limited to idealized materials, simple geometrical configurations and specific loading histories [1–6]. These fracture parameters may also be obtained by experimental measurements Corresponding author. Tel.: +1 734 657 9805; fax: +1 562 985 0786.

E-mail address: [email protected] (D. Xie). 0734-743X/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijimpeng.2006.02.007

ARTICLE IN PRESS 1048

D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

[7–12]. However, the experiment procedure is very time consuming and expensive and is not applicable to most engineering materials. Due to rapid and extensive developments in both computational hardware and software, numerical methods are becoming economical and feasible for performing dynamic fracture analyses on practical structures with novel materials, having complex geometry and under varied loading history. Atluri and Nakagaki [13] even suggest that numerical methods are mandatory to deal with practical situations of fractures. Many numerical methods such as finite difference method (FDM) [14,15], boundary element method (BEM) [16–21] and meshless Petrov-Galerkin method (MPGM) [22,23] have been attempted. With their own finite element method (FEM) codes, Aberson et al. [24] and Aturi and Nishioka [25,26] developed singular elements to compute dynamic stress intensity factors (DSIF) for both propagating and stationary cracks. Sun and Jih [27] used the virtual crack closure technique (VCCT) to compute dynamic strain energy release rates. Their work is invaluable in providing a basic understanding for dynamic fracture mechanics. However, these codes have not been commercialized and therefore, are not available to engineers who need to perform dynamic fracture analysis in an industrial setting. On the other hand, the universal finite element analysis (FEA) software which is widely and successfully used in industry may not provide fracture parameters. Extra post-processing is commonly needed to calculate fracture parameters based on the FEA output primary variables such as nodal displacements and forces. To fill the gap between the commercial FEA and fracture mechanics analysis, this paper describes an interface element developed to accommodate fracture analysis as a user modification to FEA commercial software. The interface element uses VCCT as described by the authors [28–31]. The VCCT was initially proposed by Rybicki and Kanninen [32]. It is one of the most efficient and robust tools to compute strain energy release rates within one step analysis without using singular elements. The proposed interface element has five nodes with a very stiff spring between the node pair at the crack tip to calculate the internal forces. Three dummy nodes are introduced to extract information for displacement openings behind the crack tip and the virtual crack jump ahead of the crack tip. All the components required to calculate the strain energy release rates (GI and GII) based on VCCT are available within the interface element and therefore, the strain energy release rates can be calculated and output simultaneously as the FEA software performs conventional analysis. It has been successfully implemented within ABAQUSs through a user element UEL for crack kinking [28] and crack growth in 2D [29] and 3D [30,31] under static loading. This paper further validates the interface element for stationary cracks subjected to the dynamic loading. Five examples of 2D panels with different crack locations and orientations under impulsive loading (Heaviside function) were examined. The dynamic stress intensity factors converted from dynamic strain energy release rates output from the present interface element are compared with those obtained by others with FDM, BEM or FEM. Further, Kalthoff tests on the notched bend specimen under impact loading were fully simulated. The losts of contact (LOC), the most significant observation in the test, was captured by applying contact condition between the load and support rollers and the specimen. The dynamic stress intensity factors calculated by the present interface element with ABAQUSs agree well both quantitatively and qualitatively with those measured by Kalthoff using the optical method of caustics. Therefore, this interface element technique is a powerful tool to perform failure analysis during crash and collision for engineering structures at industrial level via commercial FEA software. 2. Interface element based on VCCT The VCCT was first proposed for 2D crack configurations by Rybicki and Kanninen [32] and was extended to three dimensions (3D-VCCT) by Shivakumar et al. [33]. Recently, Krueger [34] presented a summary of historical development and a discussion with respect to different applications. Fig. 1 shows the definition and node numbering of the VCCT interface element for 2D fracture problems. When the element is applied, it is placed in such a way that the nodes 1 and 2 are located at the crack tip, with nodes 3 and 4 behind and node 5 ahead of the crack tip. The element contains two sets of node groups: the top set (nodes 1, 3 and 5) and the bottom set (nodes 2 and 4). When the element is placed at the crack tip between conventional body elements, the top set nodes coincide with their corresponding nodes in the bottom set so that the gap between the top node set and the bottom node set, exaggerated in Fig. 1, will vanish. Since the

ARTICLE IN PRESS D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

1049

Body mesh

Δa 3

1

5 Interface element

4

2

Body mesh

Y

X Fig. 1. Definition of interface element and its node numbering when crack is oriented along the X-axis.

element has five nodes and each nodes has two degrees of freedom, the full instantaneous displacement array for the element is {U1, U2, U3, U4, U5, U6, U7, U8, U9, U10} transmitting between ABAQUSs main code and its user element subroutine UEL. However, only those components related to the crack tip, {U1, U2, U3, U4}, are updated in UEL. In order to calculate the nodal forces at the crack tip, very stiff spring is placed between nodes ‘‘1’’ and ‘‘2’’. The strain energy within the spring is 1 1 U ¼ K x ðU 1  U 3 Þ2 þ K y ðU 2  U 4 Þ2 2 2

(1)

and its variation is dU ¼ ðdU 1  dU 3 ÞK x ðU 1  U 3 Þ þ ðdU 2  dU 4 ÞK y ðU 2  U 4 Þ ¼ dUT KU,

(2)

where 2

Kx

6 6 0 K¼6 6 K x 4 0

0

K x

Ky 0

0 Kx

K y

0

0

3

7 K y 7 7; 0 7 5 Ky

9 8 U1 > > > > > > < U2 = ; U¼ > > > U3 > > > ; : U4

9 8 dU 1 > > > > > > < dU 2 = dU ¼ , > > > dU 3 > > > ; : dU 4

(3)

where ‘‘K’’ is the stiffness matrix for the element and ‘‘U’’ the displacement array related to nodes 1 and 2. Therefore, the nodal forces at the crack tip are: F x ¼ K x ðU 1  U 3 Þ;

F y ¼ K y ðU 2  U 4 Þ

(4)

The force is updated only through U1 to U4 within the element. Nodes 3, 4, 5 are dummy nodes and do not have contributions to the stiffness matrix. They are introduced only to extract information for displacement opening behind crack tip and crack jump length ahead of crack tip. Since nodes 3 and 4 are behind the crack tip, the displacement openings are: Du ¼ U 5  U 7 ;

Dv ¼ U 6  U 8 ,

(5)

The crack jump length is the distance between nodes 1 and 5 and is therefore calculated by: Da ¼ jx5  x1 j in the coordinate system shown in Fig. 1.

(6)

ARTICLE IN PRESS D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

1050

Δa 1

Y X θ

Y

5 θ

2

3 4

X Fig. 2. Interface element for inclined crack lying in (X, Y) plane.

Based on one-step-analysis 2D-VCCT, the strain energy release rates can be approximated as the product of the nodal forces at the crack tip and the nodal displacement openings behind the crack tip. That is F y Dv F x Du , (7) ; GII ¼ 2BDa 2BDa where B is the thickness of the body. In practical cases, due to the orientation of the structural geometry or large deformation during loading, the crack may lie in the (X, Y) plane while not aligning with the X-axis shown in Fig. 2. In this case, some simple modifications are necessary. The crack jump length is now qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 Da ¼ ðx5  x1 Þ2 þ y5  y1 . (8) GI ¼

The included angle is determined by x5  x1 y  y1 cos y ¼ ; sin y ¼ 5 (9) Da Da Then the nodal forces and the displacement openings are projected into local coordinate system (X¯ ; Y¯ ) as F¯ x ¼ F x cos y þ F y sin y, F¯ y ¼  F x sin y þ F y cos y

ð10Þ

D¯u ¼ Du cos y þ Dv sin y, D¯v ¼  Du sin y þ Dv cos y.

ð11Þ

and

The strain energy release rates are calculated by F¯ y D¯v F¯ x D¯u . (12) ; GII ¼ 2BDa 2BDa However, these modifications are not valid in crack kinking problems since a kinking crack is not a self-similar one. The kinking VCCT was developed in Ref. [28]. The above strain energy release rates can be directly output through ABAQUSs state variables ‘‘SDV’’. Therefore, they can be obtained simultaneously with other results as ABAQUSs performs conventional finite element analysis without extra post-processing. GI ¼

3. Validation examples To demonstrate the capability and reliability of the VCCT interface element for solving dynamic crack problems in linear fracture mechanics, 2D crack problems with existing numerical or experimental results are solved and results are compared. In order to compare with the dynamic stress intensity factors (K dyn I ðtÞ and

ARTICLE IN PRESS D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

1051

dyn dyn K dyn II ðtÞ) available in the literature, they are related to the strain energy release rate (G I ðtÞ, G II ðtÞ) for isotropic materials as follows [35]. For plane stress

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¯ K dyn ðtÞ ¼ signð F Þ EG dyn y I I ðtÞ,

(13a)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi EG dyn II ðtÞ

(13b)

¯ K dyn II ðtÞ ¼ signðF x Þ and for plane strain K dyn I ðtÞ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E ¯ ¼ signðF y Þ G dyn ðtÞ, 1  v2 I

(14a)

K dyn II ðtÞ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E ¼ signðF¯ y Þ G dyn ðtÞ, 1  v2 II

(14b)

σ (t)

where E and v are the Young’s Modulus and Poisson’s ratio, respectively. The above equations have the same form as their counterparts under static loading. The body geometry was modeled with 2D ABAQUSs standard elements such as ‘‘CPS4’’ or ‘‘CPE4’’ unless stated otherwise. Neither special singular elements nor the collapsed element technique is used at the crack tip. The interface element is placed at the location of crack tip and is used as easy as the ABAQUSs standard element. An implicit dynamic analysis procedure was employed and no convergence difficulty was encountered. No damping is included in the model.

σ0=0.4GPa

t Fig. 3. Heavisde step function type impulsive loading.

σ (t)

σ (t)

20 mm

kI

3.0

2.0

FDM BEM Present

1.0

kIsta

0.0

40 mm -1.0

Crack length 2a=4.8mm (a)

(b)

0

5

10 t/μs

Fig. 4. Central crack in a rectangular plate (a) Geometry (b) Curve of KI vs. t.

15

ARTICLE IN PRESS 1052

D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

3.1. Examples under impulsive loading In this section, five examples of 2D panels with different crack locations and orientations are examined. They are all subjected to Heaviside step-function type impulsive loading with magnitude of s0 ¼ 0:4 MPa, see Fig. 3. For all the examples in this section, the time increment is uniform with Dt ¼ 0:2 ms. One of the extensively treated dynamic crack problems is a rectangular panel with a central crack, see Fig. 4(a). The panel is assumed to be linear elastic material with the mechanical properties: shear modulus G ¼ 76:923 GPa, Poisson’ ratio v ¼ 0:3 and mass density r ¼ 5 Mg=m3 . Due to symmetry, only one quarter of the panel is modeled by 100  50 elements. Plane strain element ‘‘CPE4’’ were used. The mode I dynamic stress intensity factor K dyn using Eq. (14a) and then normalized by a factor of I ðtÞ is obtained pffiffiffiffiffiby ffi K0:kI ðtÞ ¼ K dyn pa . The normalized dynamic stress intensity factor, kI(t), is plotted ðtÞ=K , where K ¼ s 0 0 0 I against time in Fig. 4(b) using open circular dots. The results from FDM [14,15] and BEM [18] are also shown in Fig. 4(b) with solid and dashed lines. In FDM, the linear extrapolation scheme was used to compute the mode I dynamic stress intensity factor K dyn I ðtÞ. It is clear that the three agree well both quantitatively and

3.0 FDM Present 2.0

kIsta

1.0

σ (t)

σ (t)

20 mm

kI

32 mm

0.0 0

Crack length a=4.8mm

(a)

4

8

(b)

12 t/μs

16

20

Fig. 5. Edge crack in a rectangular plate (a) Geometry (b) Curve of KI vs. t.

3.0 FDM Present

kIsta

1.0

σ (t)

20 mm

kI

2.0

32 mm 0.0 0

Crack length 2a = 4.8 mm

(a)

5

10

15 t/μs

(b) Fig. 6. Central crack in a rectangular plate with end clamped (a) Geometry (b) Curve of KI vs. t.

20

25

ARTICLE IN PRESS D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

1053

qualitatively. The normalized static stress intensity factor for the same problem ( kst I ¼ 1:038) is also shown as the horizontal dashed line. The second example is a rectangular panel with a centrally located edge crack, as shown in Fig. 5(a). The panel has the same material properties as those of the previous example. Due to symmetry, only the right half of the panel is modeled by a 40  50 mesh of ‘‘CPE4’’ plane strain elements. The normalized mode I dynamic stress intensity factor kI(t) is plotted against time in Fig. 5(b). The FDM result is also shown in Fig. 5(b) by the solid line [15]. The normalized static stress intensity factor for this example is also shown in Fig. 5(b) as the horizontal dashed line. Next, consider a clamped rectangular plate with a centrally located crack as shown in Fig. 6(a), with a step function load applied to the right boundary in a plane strain condition. This plate has the same material properties as those of the previous examples. No symmetry is used here and the whole panel is modeled with 80  50 ‘‘CPE4’’ elements. The normalized mode I dynamic stress intensity factor kI(t) is plotted against time

2.0 FDM Present

1.6

kI and kII σ (t)

45°

30mm

σ (t)

1.2

kI

0.8 0.4 0.0

kII

60mm

-0.4

Crack length 2a=14.14m (a)

0 (b)

5

10 t/μs

Fig. 7. Inclined crack in a rectangular plate (a) Geometry (b) Curves of kI and kII vs. t.

Fig. 8. FEA mesh before and after deformation for the inclined crack example.

15

20

ARTICLE IN PRESS 1054

D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

in Fig. 6(b) [15]. It is qualitatively quite similar to the results for the case of a central crack in a rectangular panel except that the period of the standing wave is doubled. This is due to the asymmetric loading. The horizontal dashed line Fig. 6(b) is the normalized static stress intensity factor for this case. The fourth example in this section is a rectangular plate of length 60 mm and width 30 mm containing a central crack of length 14.14 mm inclined at an angle a ¼ 451, as shown in Fig. 7(a). The material properties are the same as in the previous examples. The whole plate is modeled with ‘‘CPE4’’ elements for the major part while ‘‘CPE3’’ are used to accommodate the inclined crack. The normalized DSIF kI(t) and kII(t) are plotted in Fig. 7(b) and compared with those of Aliabadi [18] who used the domain formulation and a sub-region technique in the BEM. The solutions obtained by the two methods are similar. Fig. 8 shows the mesh before deformation and at several values of time. Finally, consider the rectangular panel with a centrally located crack shown in Fig. 9(a). The material properties are taken to be: shear modulus G ¼ 29.4 GPa, Poisson’ ratio v ¼ 0:286 and mass density r ¼ 2:45 Mg=m3 . Due to symmetry, only a quadrant of the panel is modeled with 104  40 ‘‘CPE4’’ elements.

σ (t)

40mm

24mm

104mm

σ (t)

(a)

2.5

2.0

kIdyn (t)

1.5

1.0 Atluri Sun

0.5

Present

0 (b)

0

4

8

12 t/μs

Fig. 9. Central crack in a rectangular plate.

16

20

ARTICLE IN PRESS D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

1055

The presently computed normalized dynamic stress-intensity factor solution is shown in Fig. 9(b). The numerical solutions by Atluri and Nishiko [26] and Sun and Jih [27], also shown in Fig. 9(b), are found to be nearly identical to the results developed here. 3.2. Examples under impact loading The Charpy or drop weight test of pre-cracked specimens is used extensively for dynamic fracture toughness determination over a wide range of impact velocities [36]. It is known that the inertia effects become significant when the specimens are subjected to impact loading and must be taken into account in analysis. Aberson et al. [37] attempted a numerical analysis for a three point bending (TPB) specimen of steel subjected to a falling load with the loading history shown in Fig. 10. Fig. 11 shows the dynamic stress intensity factors (solid line) by Aberson et al. [37] and the static stress intensity factor (dash line) obtained from pffiffiffi   a 2  a 3  a 4  3PS a a K¼ þ 14:7 1:93  3:12  25:3 þ 25:9 . (15) W W W W 2BW 2 Here, S is the beam span (250 mm), B the thickness (25 mm), W the width (300 mm), and a the length of the crack (25 mm). The dynamic and static stress intensity factors obtained using the present interface element with ABAQUSs are also shown in Fig. 11 with circular and triangular dots, respectively. The excellent agreement contributes validation of the present approach. The fracture toughness was calculated by substituting the peak load Pmax into Eq. (15) for p the ffiffiffiffi static stress intensity factor. For the peak load Pmax ¼ 55.6 kN, one obtains the value KIC ¼ 43.7 MPa m. When the specimen inertia is taken into account, the peak stress intensity factor is elevated by more than 8% compared to the value obtained without considering the inertia effects, and is attained earlier than the peak of the quasistatic value. Hence, the results of dynamic analysis are in contradiction with the assumption that the crack begins to propagate at the instant when load attains its peak value. This example illustrates the need to take into account the inertia effects while analyzing the experimental results. Kalthoff tests [7,8] have yield more information on the understanding of this physical process. Some analytical models, based on beam and lumped mass-spring models, were introduced to interpret the experimental observations [38,39]. In the present study, these tests were fully modeled at a structure level.

60 Pmax=55.6kN at t=0.575ms

40 P (t), kN

P (t)

20

0 0

0.1

0.2

0.3 t/μs

Fig. 10. Loading history.

0.4

0.5

0.6

ARTICLE IN PRESS 1056

D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060 50

KIdyn (t), MPa·m1/2

40

30

20

10

FDM Present

0 0.1

0

0.2

0.3 t/μs

0.4

0.5

0.6

Fig. 11. Dynamic and quasi-static stress intensity factors.

Instrumented Tup, PH

Instrumented Anvil, PA

Incident Light Fig. 12. Experimental set-up.

Fig. 12 shows the experimental set-up used by Kalthoff. Specimens made from the epoxy resin Araldite B were impacted by a drop weight of 4.9 kg mass at a velocity of 1m/s. The specimens were 412 and 550 mm long, 100 mm wide and 10 mm thick. A fixed support span of 400 mm was used throughout the investigations. The length of the initial notches was 30 mm. The response of the specimen during the impact process was recorded in terms of the impact load, PH(t), the anvil reaction force, PA(t), and the dynamic stress intensity factor at the notch tip, K dyn I ðtÞ. The loads PH(t) and PA(t) were determined by strain gages mounted on the tup of the striking hammer and the anvil. The K dyn I ðtÞ was measured by means of the shadow optical method of caustics. Fig. 13 shows the finite element model used in present study. The body of specimen was modeled with ABAQUSs standard plane stress ‘‘CPS4’’ elements with material properties of Araldite B. Both the impact hammer and the anvil were modeled with ‘‘CPS4’’ and ‘‘CPS3’’ plane stress elements. The anvil is assumed to be steel. The hammer has the properties of steel except for the mass density. The equivalent hammer density,

ARTICLE IN PRESS D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

Fig. 13. Finite element mesh.

3 Test Present

PH/kN

2

1

0 Test Present

KIdyn/MPa·m1/2

3

2

1

0

PA/kN

Test 1.0

Present

0.5 0.0 0

0.2

0.4

0.6

0.8 1.0 Time, t/ms

1.2

1.4

Fig. 14. Results for short specimens (W ¼ 412 mm).

1.6

1.8

1057

ARTICLE IN PRESS 1058

D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

the total mass of the hammer divided by the volume the one-quarter cylinder, is used in the finite element model to simplify modeling of the hammer. Contact conditions were applied between the hammer and the specimen and also the anvil and the specimen. The contact force between the hammer and specimen in the model is comparable to PH(t) while the contact force between the anvil and specimen is PA(t). The present interface element is placed at the crack tip to calculate dynamic strain energy release rates and then the dynamic stress intensity factors K dyn I ðtÞ. Results obtained for the specimens 412 mm long and 550 mm long are shown in Figs. 14 and 15, respectively. The contact forces agree qualitatively with the test results and capture the LOC phenomenon which is one of the most significant observations in the Kalthoff’s tests. The differences between the plots of the forces are partly due to the fact that the contact forces from the model are the global responses between the specimen and the hammer and the anvil while the loads measured through use of the strain gauges are localized at a location on the hammer and the anvil. Through the more detailed modeling of both the hammer and the anvil, the differences may vanish. For both cases, the dynamic stress intensity factors have excellent

3 Test Present

PH/kN

2

1

0 Test Present

KIdyn/MPa·m1/2

3

2

1

0

PA/kN

Test 1.0

Present

0.5 0.0 0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

Time, t/ms Fig. 15. Results for elongated specimens (W ¼ 550 mm)

1.6

1.8

ARTICLE IN PRESS D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

1059

agreement with the measurements by the optical method of caustics [8]. This example further validates the present interface element for VCCT. 4. Conclusions This paper describes an interface element for 2D-VCCT for use in dynamic fracture analysis. The element has been implemented into commercial FEA software ABAQUSs via the user element subroutine UEL. The element is evaluated with five examples of impulsive loading and two examples of impact loading. The results from the present interface element agree well with the available numerical and experimental results obtained by the others. With this interface element, fracture mechanics can be directly applied to crack growth problem using commercial FEA software. The element has several significant advantages. It does not require extra post-processing to extract fracture parameters nor does it place any special burden on definition of the body mesh. The user can apply it in conjunction with conventional finite elements. Finally, with the interface element proven to be reliable in a variety of cases, it can be employed with confidence in other cases of dynamic fracture. In summary, the present 2D-VCCT interface element is simple, efficient, and robust. The procedure should be also applicable to other commercial FEA software codes with user subroutines. The stress intensity factor and the J-integral can also be calculated simultaneously with commercial FEA software if the linear extrapolation procedure and the domain integral technique are implemented through the idea of the interface element. Therefore, this interface element should prove to be a powerful means of analyzing fracture behavior of engineering structures in industry via commercial FEA software. Acknowledgements Financial support for this work was provided by NASA, the State of South Carolina, and Clemson University through the EPSCoR grant ‘‘Development and Enhancement of Research Capability for Aircraft Structures and Materials.’’ References [1] Chen EP, Sih GC. Transient response of cracks to impact loads. In: Sih GC, editor. Elastodynamic crack problems. Leider: Noordhoff Publishing; 1976. p. 1–58. [2] Freund LB. Dynamic Fracture Mechanics. Cambridge, New York: Cambridge University Press; 1990. [3] Broberg KB. Cracks and fracture. San Diego: Academic Press; 1999. [4] Achenbach JD. Wave propagation in elastic solids. Amsterdam/New York: North-Holland: American Elsevier; 1973. [5] Achenbach JD, Li ZL. Propagation of horizontally polarized transverse-waves in a solid with a periodic distribution of racks. Wave Motion 1986;8:371–9. [6] Itou S. Dynamic stress intensity factors around two rectangular cracks in an infinite elastic plate under impact load. Mech Res Commun 2002;29:225–34. [7] Kalthoff JF. On the measurement of dynamic fracture toughness—a review of recent work. Int J Fract 1985;27:277–98. [8] Bohme W, Kalthoff JF. The behavior of notched bend specimens in impact testing. Int J Fract 1982;20:R139–43. [9] Liu C, Knauss WG, Rosakis AJ. Loading rates and the dynamic initiation toughness in brittle solids. Int J Fract 1998;90:103–18. [10] Kobayashi AS, Ramulu M. Dynamic stress-intensity factors for unsymmetric dynamic isochromatics. Exp Mech 1981;21:41–8. [11] Der VK, Barker DB, Holloway DC. Split birefringent coating technique to determine dynamic stress intensity factors. Mech Res Commun 1978;5:313–8. [12] Katsamanis F, Raftopoulos D, Theocaris PS. Static and dynamic stress intensity factors by method of transmitted caustics. J Eng Mater Technol 1977;99:105–9. [13] Atluri NS, Nakagaki M. In: Alturi SN, editor. Computational methods in the mechanics of fracture. North-Holland: Amsterdam; 1986. p. 169–228. [14] Chen YM. Numerical computation of dynamic stress intensity factors by a Lagrangian finite-difference method (The HEMP code). Eng Fract Mech 1975;7:653–60. [15] Chen YM, Wilkins ML. Numerical analysis of dynamic crack problems. In: Sih GC, editor. Elastodynamic crack problems. Leiden: Noordhoff Publishing; 1976. p. 295–345. [16] Mukhopadhyay NK, Maiti SK, Kakodkar A. A review of SIF evaluation and modeling of singularities in BEM. Comput Mech 2000;25:358–75. [17] Wen PH, Aliabadi MH, Young A. A time-dependent formulation of dual boundary element method for 3D dynamic crack problems. Int J Numer Methods Eng 1999;45:1887–905.

ARTICLE IN PRESS 1060

D. Xie, S.B. Biggers Jr. / International Journal of Impact Engineering 34 (2007) 1047–1060

[18] Fedelinski P, Aliabadi MH, Rooke DP. Boundary element formulations for the dynamic analysis of cracked structures. Eng Anal Boundary Elements 1996;17:45–56. [19] Israil ASM, Dargush GF. Dynamic fracture mechanics studies by time-domain BEM. Eng Fract Mech 1991;39:315–28. [20] Sladek J, Sladek V, Dargush GF. A boundary integral equation method to dynamic crack problems. Eng Fract Mech 1985;21:269–77. [21] Mackerle J. Fracture mechanics parameters and finite element and boundary element methods—a bibliography (1992–1004). Finite Elements Analysis Design 1995;19:209–23. [22] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin methods. Int J Numer Methods Eng 1996;39:923–38. [23] Atluri SN, Shen S. The meshless local Petrov-Galerkin (MLPG) method. Forsyth, GA: Tech Science Press; 2002. [24] Aberson JA, Anderson JM, King MM. Dynamic analysis of cracked structures using singularity finite elements. In: Sih GC, editor. Elastodynamic crack problems. Leiden: Noordhoff Publishing; 1976. [25] Nishioka T. Recent developments in computational dynamic fracture mechanics. In: Aliabadi MH, editor. Dynamic fracture mechanics. Southampton: Computational Mechanics Publications; 1995. [26] Nishioka T, Atluri SN. Numerical Modeling of dynamic crack propagation in finite bodies, by moving singular elements. J Appl Mech. [27] Jih CJ, Sun CT. Evaluation of finite element based crack-closure method for calculating static and dynamic strain energy release rates. Eng Fract Mech 1990;37:313–22. [28] Xie D, Waas AM, Shahwan KW, Schroeder JA, Boeman RG. Computation of energy release rates for kinking cracks based on virtual crack closure technique. CMES-Comput Modeling Eng Sci 2004;6:515–24. [29] Xie D, Biggers SB. Progressive crack growth analysis using interface element based on the virtual crack closure technique. Finite Elements in Analysis and Design, in press. [30] Xie D, Biggers SB. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on virtual crack closure technique, Part I: formulation and validation. Eng Fract Mech 2006;73:771–85. [31] Xie D, Biggers SB. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on virtual crack closure technique. Part II: Sensitivity study on modeling details. Eng Fract Mech 2006;73:786–801. [32] 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. [33] 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. [34] Krueger R. The virtual crack closure technique: history, approach and applications. NASA/CR-2002-211628. [35] Kobayashi T. Analysis of impact properties of A533 steel for nuclear reactor pressure vessel by instrumented Charpy test. Eng Fract Mech 1984;19:49–79. [36] Kobayashi T. Progress in the instrumented Charpy impact test. Mater Sci Res Int 2002;8:141–50. [37] Aberson JA, Anderson JM, King WW. A finite element analysis of an impact test. Advance Research on Strength and Fracture Material, 4th International Conference on Fracture Waterloo, 1977, Vol. 3A, 1978, p. 85-80. [38] Williams JG. The analysis of dynamic fracture using lumped mass-spring models. Int J Fract 1987;33:47–59. [39] Marur PR. Dynamic analysis of one-point bend impact test. Eng Fract Mech 2000;67:41–53.