# Explicit Dynamic Finite Element Method for Predicting Implosion/Explosion Induced Failure of Shell Structures

Authors: Jeong-Hoon Song; Patrick Lea; Jay Oswald
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A simplified implementation of the conventional extended finite element method (XFEM) for dynamic fracture in thin shells is presented. Though this implementation uses the same linear combination of the conventional XFEM, it allows for considerable simplifications of the discontinuous displacement and velocity fields in shell finite elements. The proposed method is implemented for the discrete Kirchhoff triangular (DKT) shell element, which is one of the most popular shell elements in engineering analysis. Numerical examples for dynamic failure of shells under impulsive loads including implosion and explosion are presented to demonstrate the effectiveness and robustness of the method.


---

## Body

## 1. Introduction

In this work, we describe a method for modeling fractured discrete Kirchhoff triangular (DKT) shell elements [1] based on the extended finite element method (XFEM) [2, 3]. One of the underlying key concepts in the XFEM is the partition of unity approach [4, 5]. In the partition of unity approach, the approximation basis is spanned by the standard finite element approximation space and extended by the products of the standard finite element shape functions with special local characteristic functions which are constructed from knowledge about the solution. The XFEM has since been developed for two-dimensional dynamic fracture problems [6, 7], static shell fractures [8], and implicit dynamic shell fracture problems [9].Even though numerous references are available for continuum shell elements, the literature on dynamic crack propagation in shells is quite limited. Cirak et al. [10, 11] have developed a method for dynamic crack propagation in Kirchhoff type shells based on interelement cohesive crack methods [12–14]; in the interelement cohesive crack methods, the crack is limited to propagation along the element edges with local remeshing.Mehra and Chaturvedi [15] used the smooth particle hydrodynamics (SPH) method for simulations of tearing of the thick plates. Combescure et al. [16] and Maurel and Combescure [17] recently developed SPH shell formulations for explicit dynamic method and successfully applied the method to the prediction of dynamic fractures in shell structures.The described implementation scheme is mainly based on the XFEM, but its actual implementation follows the phantom node method [7, 18] that has been developed by the author of this paper. In this approach [7, 18], the element which contains the crack is replaced by two superposed elements with additional nodes. Though this discontinuity representation scheme uses the same linear combination of enrichment functions as the conventional XFEM, it allows for considerable simplifications in fractured thin shell element formalisms and furthermore is applicable to arbitrary large deformations.An elementwise progression of the crack is also employed; that is, the crack tip is always on an element edge. The elementwise crack propagation scheme may cause some noise during the crack propagation with coarse meshes. However, in Song et al. [7], it is shown that such noise diminishes with mesh refinements, and the crack propagation speeds converge to the progressive crack propagation results [6, 19, 20].


## 2. Discrete Kirchhoff Triangular Shell Element

The main advantage of the DKT shell element is that a mesh can easily be generated from any kind of surfaces. The geometry of the element is described by three linear shape functions in the reference coordinates. The kinematic of the DKT shell elements is described by superimposing the membrane, the bending, and the rotational (drilling) behavior of shells with different corresponding degrees of freedoms (DOFs) as shown in Figure1.Kinematic data of the DKT triangular element: (a) the corotational coordinates, (b) the bending degrees of freedom, (c) the in-plane membrane degrees of freedom, and (d) the drilling degrees of freedom.
(a)
(b)
(c)
(d)However, for further explanation on the salient features of the DKT element, henceforth, we will useβ instead of θ for the rotation as shown in Figure 2; note that this is only for a clear illustration purpose. In this notation, βx is the rotation whose displacement is along the axis x.Positive directions ofβx and βy and correspondence between the rotations θx and θy and βx and βy.
(a)
(b)The discrete Kirchhoff assumption [1] is the following: the normal rotations must be equal to the first derivative of the transverse displacement. These constraints are imposed at a discrete number of points, which leads to the relation between the normal rotations and the displacements at the element joints. Let us consider the constraints in terms of x and y directions as
(1)βx(ξ)+wx=0,βy(ξ)+wy=0,
where βx=θy and βy=-θx. This sign difference is due to the orientation of the rotation, which will generate in-plane displacements. However, with a linear discretization, these conditions cannot be verified within the entire domain of the shell finite element. One alternative approach is that they can be only verified at some discrete parts of the shell elements, such as a midpoint of each side of the DKT elements; this is the implication of the discrete Kirchhoff assumption in the DKT shell elements.To verify the discrete Kirchhoff assumption, one has to add additional shape functions which do not change the nodal values of any field but only are allowed to modify the values on the midpoint. Thus, the rotational DOFs are discretized by(2)βx=∑i=13Ni·βxi+∑i=46Pi·Ci·αi,(3)βy=∑i=13Ni·βyi+∑i=46Pi·Si·αi.
As shown in Figure 3, the discrete Kirchhoff constraints along each side i-j are introduced at the midpoint k:
(4)βs(ξ)+ws=0,
where s is the absciss along the side of the element.Figure 3
Geometry and local tangential-normal coordinate system of the DKT element.The different functionsPi appearing in (2) and (3) are given by
(5)P4=4(1-ξ-η)ξ,(6)P5=4ηξ,(7)P6=4(1-ξ-η)η.
Note that, in (2) and (3), C and S define cosinus and sinus values of the current geometry in the reference coordinates as shown in Figure 4. The particular values for the reference triangular element are C4=1, S4=0, C5=-2/2, S5=2/2, C6=0, and S6=-1. These values are the cosinus and sinus of the vector representing the side of the reference triangular element. Thus, for example, at the node 4 in Figure 3, we can show that x is the sames in terms of substituting (2) and (4). Furthermore, for this particular edge that has the midpoint at (0,1/2), we could further develop the following relation:
(8)-ws=βs(12,0)=βx(12,0)=βx1+βx22+α4=-wx=-w2-w12.
So, the expression of α4 is given as
(9)α4=w1-w2-βx12-βx22.
Similarly, we can determine the other unknowns α5 and α6 as
(10)α6=w1-w3-βy12-βy32,α5=-2w3+2w2+24(βx2+βx3-βy2-βy3).
However, the other three unknowns α4, α5, and α6 in the formalism of the rotation should be determined by solving system equations given by
(11)4·[1-2ξ-η-22η0022ξξ+2η-1-ξ224(η-ξ)η2]·[α4α5α6]=[βx,xβy,y12(βx,y+βy,x)].
Throughout this procedure, that is, determining three unknowns α4, α5, and α6, the Kirchhoff assumption is enforced in the middle of each side of the shell element. Note that the nodal value is still intact because of the shape functions Pi which are zero at the nodes.Figure 4
Representation of the crack by the XFEM and the phantom node method with DKT elements.


## 3. Representation of Fractured Shell Element

Based on the phantom node approach [7, 18], the DKT shell element which contains a crack is replaced by two superposed DKT elements with additional phantom nodes as shown in Figure 4.As with the standard approach to phantom nodes [7, 18], cracks will be inserted elementwise at propagation and the crack surface will be limited to normal to the shells midsurface. While the equations here will be for the phantom node method, the equivalence between phantom node method and XFEM is shown by Song et al. [7].The discontinuous velocity fields in the midsurface of the fractured shell elements can be described by(12)vmid=∑i=13Ni(H(-f(x,t)))vi+∑j=13Nj(H(f(x,t)))vj,
where i and j are the first and second halves of the pair of overlapping element, respectively, f(x) is an implicit function that can describe the geometry of the crack surface in the midsurface of the shell with iso-zero line, that is; crack path in the midsurface of shell can be represented by f(ξ)=0, and H(x) is the step function given by
(13)H(x)={1,x≥0,0,x<0.
The velocity fields can also be expressed in corotational coordinates as is common in shell, but we have omitted it for brevity. The velocity equation for the DKT shells does not vary from that of the other shell elements but the rotations do. Keeping with the notation from the last section the rotations in x and y can be expressed as
(14)βx=∑i=13Ni(H(-f(x,t)))βxi+∑i=46Pi(H(-f(x,t)))Ci·αi+∑j=13Nj(H(-f(x,t)))βxj+∑j=46Pj(H(-f(x,t)))Cj·αj,(15)βy=∑i=13Ni(H(-f(x,t)))βyi+∑i=46Pi(H(-f(x,t)))Si·αi+∑j=13Nj(H(-f(x,t)))βyj+∑j=46Pj(H(-f(x,t)))Sj·αj,
where again an overlapping pair (elements i and j) is used. One advantage of using the phantom node method with DKT shells is that the procedure for finding α’s does not need to be varied from the standard method. This greatly simplifies the implementation of the XFEM for DKT elements, leaving (5)–(11) unchanged.


## 4. Computation Procedures

### 4.1. Time Integration: Newmark Scheme

In this work, Newmark scheme for the explicit time integration is used. The time integration procedure is written as(16)Ut+Δt=Ut+Δt·U˙t+Δt22U¨t,(17)U˙t+Δt=U˙t+12Δt(U¨t+U¨t+Δt),(18)M·U¨t+Δt=Fext-Fint,
where Ut (U˙t and U¨t, resp.) denotes the displacement (velocity and acceleration, resp.) at time t. Δt is the time integration step, M is the mass matrix, and Fext (Fint, resp.) is the external (internal, resp.) forces at time t.A diagonal mass matrix is frequently used in this explicit time integration scheme because it allows us to avoid a matrix inversion for solving (14); that is, no matrix inversion appears in this scheme. Consequently, the main advantage of using explicit time integration scheme is to speed up the computation and use less memory by storing only vectors instead of matrices to the computer.However, this explicit integration scheme is conditionally stable, and the stability condition is defined in terms of a maximum time stepΔtc; we usually name it a critical time step since it is the largest time step that can be used. The critical time step is evaluated from the eigenvalue analysis with the mass M and stiffness matrix K since
(19)Δtc=2ωmax,
where ωmax is the maximum frequency determined by solving eigenvalue of problems in (16):
(20)det(K-ω2M)=0.
The stability of explicit time integration for the XFEM is defined by the same condition on the mass and the stiffness; the computation frequency must be greater than the greatest vibration frequency of the structure.


### 4.2. Computation of Lumped Mass Matrix for Cracked Elements

In the explicit dynamic analysis method, constructions of lumped mass are essential to ensure the computation of nodal accelerations without implicit solution procedures. However, the mass lumping scheme for cracked elements which employ the XFEM approach is not obvious. To circumvent such difficulties, several methods have been proposed: implicit (in cracked elements)-explicit (in continuum elements) time integration scheme [6] and modified mass lumping schemes [21, 22].In this study, the lumped mass for regular DOFs is diagonalized by the conventional row sum mass lumping technique, but, for the cracked elements, we used the mass lumping scheme that was proposed by Menouillard et al. [22]. Thus, the diagonal term i of the mass corresponding to the enriched DOFs also depends on the enrichment function H as follows:
(21)Mi=mnnode1∥Ωe∥∫ΩH2dΩe,
where m is the total actual mass of the element, nnode is the number of nodes of the element, and ∥Ωe∥ is the measure of the finite element domain Ωe. For the particular case of discontinuous enriched functions such as Heaviside function, the term of the mass matrix corresponding to an enriched node is in fact just a fraction of the regular finite element term; in other words, the lumped mass matrix for the enriched nodes is written as follows:
(22)Mi=mnnodeAiA,
where A is the volume or area of the regular element and Ai represents the fraction ratio of the cut element. One imperative advantage of this mass lumping scheme is that this method does not significantly decrease the critical time step of the continuum element [22].


## 4.1. Time Integration: Newmark Scheme

In this work, Newmark scheme for the explicit time integration is used. The time integration procedure is written as(16)Ut+Δt=Ut+Δt·U˙t+Δt22U¨t,(17)U˙t+Δt=U˙t+12Δt(U¨t+U¨t+Δt),(18)M·U¨t+Δt=Fext-Fint,
where Ut (U˙t and U¨t, resp.) denotes the displacement (velocity and acceleration, resp.) at time t. Δt is the time integration step, M is the mass matrix, and Fext (Fint, resp.) is the external (internal, resp.) forces at time t.A diagonal mass matrix is frequently used in this explicit time integration scheme because it allows us to avoid a matrix inversion for solving (14); that is, no matrix inversion appears in this scheme. Consequently, the main advantage of using explicit time integration scheme is to speed up the computation and use less memory by storing only vectors instead of matrices to the computer.However, this explicit integration scheme is conditionally stable, and the stability condition is defined in terms of a maximum time stepΔtc; we usually name it a critical time step since it is the largest time step that can be used. The critical time step is evaluated from the eigenvalue analysis with the mass M and stiffness matrix K since
(19)Δtc=2ωmax,
where ωmax is the maximum frequency determined by solving eigenvalue of problems in (16):
(20)det(K-ω2M)=0.
The stability of explicit time integration for the XFEM is defined by the same condition on the mass and the stiffness; the computation frequency must be greater than the greatest vibration frequency of the structure.


## 4.2. Computation of Lumped Mass Matrix for Cracked Elements

In the explicit dynamic analysis method, constructions of lumped mass are essential to ensure the computation of nodal accelerations without implicit solution procedures. However, the mass lumping scheme for cracked elements which employ the XFEM approach is not obvious. To circumvent such difficulties, several methods have been proposed: implicit (in cracked elements)-explicit (in continuum elements) time integration scheme [6] and modified mass lumping schemes [21, 22].In this study, the lumped mass for regular DOFs is diagonalized by the conventional row sum mass lumping technique, but, for the cracked elements, we used the mass lumping scheme that was proposed by Menouillard et al. [22]. Thus, the diagonal term i of the mass corresponding to the enriched DOFs also depends on the enrichment function H as follows:
(21)Mi=mnnode1∥Ωe∥∫ΩH2dΩe,
where m is the total actual mass of the element, nnode is the number of nodes of the element, and ∥Ωe∥ is the measure of the finite element domain Ωe. For the particular case of discontinuous enriched functions such as Heaviside function, the term of the mass matrix corresponding to an enriched node is in fact just a fraction of the regular finite element term; in other words, the lumped mass matrix for the enriched nodes is written as follows:
(22)Mi=mnnodeAiA,
where A is the volume or area of the regular element and Ai represents the fraction ratio of the cut element. One imperative advantage of this mass lumping scheme is that this method does not significantly decrease the critical time step of the continuum element [22].


## 5. Material Model and Modeling of Fracture

### 5.1. Damage Plasticity Model

A damage plasticity model that can account for the effects of stress triaxiality and Lode angle was proposed by Xue [23] and Xue and Wierzbicki [24]. In this constitutive model, damage of a material point is accessed by measuring the accumulation of the following damage increment:
(23)dD=m(εpεf)m-1dεpεf,
where D is damage parameter, εp is the plastic strain, εf is a reference strain envelope, and m is a material constant. The reference strain envelope is a function of the pressure and the Lode angle:
(24)εf=εf0μp(p)μθ(θL),
where εf0 is the initial reference strain, p is the pressure, and θL is the Lode angle. The functions μp and μθ are defined as
(25)μp(p)=1-qlog(1-pplim),μθ(θL)=γ+(1-γ)(6|θL|π)k,
where q and plim are material constants and γ and k are parameters determining the shape of the strain envelope. The weakening effect caused by the damage was also considered in this model:
(26)σ=w(D)σM=(1-Dβ)σM,
where σM is the stress of the undamaged material, w(D) is the weakening function, and β is another material constant.


### 5.2. Fracture Criterion

A critical strain based fracture criterion is used to determine the onset point of a poststrain localization behavior of a material, that is, fracture. When the strain at a crack tip material point reaches a fracture threshold, we inject a strong discontinuity ahead of the previous crack tip according to maximum principle tensile strain direction of an averaged strainεavg. For the computation of the averaged strain, εavg, we used a pointwise weighted averaging scheme which is given by
(27)εavg=∑αw(rα)εα∑αw(rα),
where w(r) is the cubic spline weight function, rα is the distance from the crack tip to the material points α, rc(≃3he) is the size of the averaging domain, and he is the size of the crack tip element; see Figure 5. Note that, for the computation of the averaged strain εavg, we only included material points which show tension dominant states as shown in Figure 5(b).Schematic of averaging domain for the evaluation of the fracture criterion: (a) the size of averaging domain and (b) possible principal stress states through the shell depth.
(a)
(b)


### 5.3. Dissipation of Fracture Energy

In this study, a cohesive crack model is prescribed along the newly injected strong discontinuity surfaces until the crack opening is fully developed, that is, until cohesive traction has vanished. The roles of a prescribed cohesive model can be summarized as follows.(1)
It can be a remedy to spurious mesh-dependent pathological behaviors by providing a bounded solution at the crack tip. For linear elastic fracture simulations, if the crack tip is not smoothly closed with cohesive forces, finite element solutions are unbounded at the crack tip due to the crack tip stress singularity and a crack path is determined by the surrounding mesh resolution. Also, for fracture in plastic bulk materials, the crack tip stress singularity can be slightly alleviated by plasticity. However, the finite element solutions still depend on the mesh resolution.(2)
If the crack opening displacement is not governed by a cohesive model, the normal stress component to the crack surface suddenly drops to zero due to lack of fracture energy dissipations; note that injecting a strong discontinuity without prescribing cohesive force is the same as creating two free surfaces without dissipating new surface initiation energies. In this case, the total system suffers from an excessive accumulation of elastic energy and this excessively accumulated energy accelerates the crack propagation speed; more discussions on the relationship between crack propagation speed and dissipated fracture energy can be found in Rabczuk et al. [27].In this study, we only prescribe the normal traction of a linear cohesive model since the early stage of crack initiation due to implosion or explosion is mostly due to mode 1 fracture behavior.


## 5.1. Damage Plasticity Model

A damage plasticity model that can account for the effects of stress triaxiality and Lode angle was proposed by Xue [23] and Xue and Wierzbicki [24]. In this constitutive model, damage of a material point is accessed by measuring the accumulation of the following damage increment:
(23)dD=m(εpεf)m-1dεpεf,
where D is damage parameter, εp is the plastic strain, εf is a reference strain envelope, and m is a material constant. The reference strain envelope is a function of the pressure and the Lode angle:
(24)εf=εf0μp(p)μθ(θL),
where εf0 is the initial reference strain, p is the pressure, and θL is the Lode angle. The functions μp and μθ are defined as
(25)μp(p)=1-qlog(1-pplim),μθ(θL)=γ+(1-γ)(6|θL|π)k,
where q and plim are material constants and γ and k are parameters determining the shape of the strain envelope. The weakening effect caused by the damage was also considered in this model:
(26)σ=w(D)σM=(1-Dβ)σM,
where σM is the stress of the undamaged material, w(D) is the weakening function, and β is another material constant.


## 5.2. Fracture Criterion

A critical strain based fracture criterion is used to determine the onset point of a poststrain localization behavior of a material, that is, fracture. When the strain at a crack tip material point reaches a fracture threshold, we inject a strong discontinuity ahead of the previous crack tip according to maximum principle tensile strain direction of an averaged strainεavg. For the computation of the averaged strain, εavg, we used a pointwise weighted averaging scheme which is given by
(27)εavg=∑αw(rα)εα∑αw(rα),
where w(r) is the cubic spline weight function, rα is the distance from the crack tip to the material points α, rc(≃3he) is the size of the averaging domain, and he is the size of the crack tip element; see Figure 5. Note that, for the computation of the averaged strain εavg, we only included material points which show tension dominant states as shown in Figure 5(b).Schematic of averaging domain for the evaluation of the fracture criterion: (a) the size of averaging domain and (b) possible principal stress states through the shell depth.
(a)
(b)


## 5.3. Dissipation of Fracture Energy

In this study, a cohesive crack model is prescribed along the newly injected strong discontinuity surfaces until the crack opening is fully developed, that is, until cohesive traction has vanished. The roles of a prescribed cohesive model can be summarized as follows.(1)
It can be a remedy to spurious mesh-dependent pathological behaviors by providing a bounded solution at the crack tip. For linear elastic fracture simulations, if the crack tip is not smoothly closed with cohesive forces, finite element solutions are unbounded at the crack tip due to the crack tip stress singularity and a crack path is determined by the surrounding mesh resolution. Also, for fracture in plastic bulk materials, the crack tip stress singularity can be slightly alleviated by plasticity. However, the finite element solutions still depend on the mesh resolution.(2)
If the crack opening displacement is not governed by a cohesive model, the normal stress component to the crack surface suddenly drops to zero due to lack of fracture energy dissipations; note that injecting a strong discontinuity without prescribing cohesive force is the same as creating two free surfaces without dissipating new surface initiation energies. In this case, the total system suffers from an excessive accumulation of elastic energy and this excessively accumulated energy accelerates the crack propagation speed; more discussions on the relationship between crack propagation speed and dissipated fracture energy can be found in Rabczuk et al. [27].In this study, we only prescribe the normal traction of a linear cohesive model since the early stage of crack initiation due to implosion or explosion is mostly due to mode 1 fracture behavior.


## 6. Numerical Examples

### 6.1. Implosion Induced Dynamic Failure of Cylinder

The cylinder is of lengthL with two extension parts of length Le at the ends; see Figure 6. Rigid plugs are perfectly bonded to the extension parts. The thickness of the cylinder is t, and the diameter is D. In the experiments, the specimens are loaded with hydrostatic pressure and the pressure was increased until the cylinders buckled. It was observed in the experiment that specimens will buckle in some specific modes, which have a dependence on the geometry.Figure 6
Setup for implosion induced failure of cylinder.In this example, we focused on predicting final fracture pattern of two experiments. The first specimen is denoted by IMP26 experiment, whereL=366.8 mm, D=38.085 mm, and t=0.701 mm. The length of the extension parts is Le=25.4 mm. The material of the cylinder is aluminum alloy (AL) 6061-T6. In the experiment, the result indicated a mode 2 implosion, but there was no fracture observed. The second specimen, denoted by IMP25, is L=143 mm plus two 25.4 mm extension parts at the ends, D=38.087 mm, and t=0.701 mm. Because it has a much shorter length than that of IMP26, IMP25 specimen buckled in a higher mode, that is, mode 3 buckling, and showed fracture on the interfaces of the end plugs and the main specimen.For numerical analysis, we modeled IMP26 and IMP25 specimens with 49200 and 24000 shell elements, respectively. The average element size is about 1 mm. The material behavior of AL6061-T6 is modeled with the damage plasticity [23, 24] with Young’s modulus E=69.5 GPa, the density ρ=2780 kg/m3, Poisson ratio ν=0.3, yield stress σy0=276 MPa, and the hardening modulus ET=634 MPa.The following coordinate system was used: thex-axis is along the axis of the cylinder, and the y- and z-axes are in the radial direction. A uniform surface pressure loading was applied on all the shell elements. Another set of concentrated forces along the axis direction were applied to the nodes on the rim of the two ends. These forces served to compress the cylinder along its longitudinal axis. The summation of the magnitude of these concentrated forces is equivalent to the total force exerted on the end plug due to the pressure. All the nodes in the two extension parts are only allowed to move along the x-axis, and all the other DOFs of these nodes, including the translational DOFs uy and uz and the rotational DOFs θx, θy, and θz, are constrained. The employment of the applied concentrated forces and the sliding boundary conditions is to model the effect of the plug.A bilinear load curve for the pressure was used in both simulations. The pressure started at zero and was increased topcr, that is, the experimental critical buckling pressure, in time t0 and then was kept constant at pcr until the end of simulation. The schematic of time history of the pressure is shown in Figure 7; pcr=1.15MPa for IMP26 experiment and pcr=2.83MPa for IMP25 experiment. We used the parameter t0=5.0×10-5 s in both simulations.Figure 7
Loading curve used in the simulation of Texas experiments.Geometry imperfection is introduced into the radius to evoke circumferential buckling easily. The actual radius with imperfection has the formr=r0(1-acosnθ), where r0 is the unperturbed radius, a is the imperfection magnitude, and θ is the circumferential angle. n is the number of wavelengths in the circumferential direction, n=2 for IMP26 simulation, and n=3 for IMP25 simulation. It should be noted that the buckling mode in numerical results does not depend on the imperfection mode. We also tested the IMP25 simulation with injection of mode 2 imperfection and obtained very similar results.Figures8(a)–8(d) show four snapshots of the numerical results for IMP26. The center of the specimen yielded first and then buckled in mode 2. The buckling region evolved toward the two ends, and the specimen entirely collapsed. a=0.1% was used in the simulation of IMP26, which is close to 0.107%, the maximum ovalization measured in the experimental specimen.Snapshots of numerical results of IMP26 experiment: (a) initial configuration and the mesh, (b) deformed configuration at timet=1.48 ms; the center of the cylinder begins to collapse, (c) deformed configuration at time t=1.78 ms; buckling region enlarges toward the two ends, and (d) final deformed configuration; the cylinder buckles in mode 2.
(a)
(b)
(c)
(d)In IMP25 simulation, we allowed an injection of the discontinuity near the interface of the main part and the extension part of the specimen. We also observed large plastic strain, large damaged or unstable material points along the central buckling lines. However, this may be due to the repulsive forces generated during the contact, so no crack was allowed to initiate in these regions.Several imperfection magnitudesa from 0.05% to 1% were tested in IMP25 simulation. The maximum oval imperfection measured in the experiment is 0.043%. It should be noted there should not be imperfection as large as 1% in the actual specimen, and this is only to examine the effect of the imperfection.Figure9 shows four snapshots of IMP25 results of a=0.05% at different times. The central part of the cylinder collapses first, and then the collapse region enlarges toward the ends. Some cracks were initiated at the ends and propagated along the circumferential direction.Snapshots of numerical results of IMP25 experiment: (a) deformed configuration at timet=0.92 ms; cylinder collapses at mode 3, (b) deformed configuration at time t=1.08 ms; cracks initiated at the two ends, (c) the final deformed configuration, (d) local zoom view of the crack, and (e) view from the axial direction.
(a)
(b)
(c)
(d)
(e)The final configurations of different imperfection magnitudes are compared in Figure10. No large effect can be observed for imperfection magnitudes though the crack opening of a=1.0% looks larger than that of a=0.05%. The time histories of velocity y of finite element node 11670 (at the center of cylinder) are plotted and compared in Figure 11. Let the collapse time be defined as the time to reach the peak velocity. It can be observed that smaller imperfection leads to longer collapse time, but the peak velocity of all the imperfections does not show an obvious dependence on the imperfection magnitude.Comparison of final configurations of cylinders with different imperfection magnitudes: (a) 0.05% imperfection, (b) 0.1% imperfection, (c) 0.5% imperfection, and (d) 1.0% imperfection.
(a)
(b)
(c)
(d)Figure 11
Comparisons of time history of nodal velocity at the center node with different amplitudes of imperfection.


### 6.2. Explosion Induced Failure of Cylinder

Chao [25] and Chao and Shepherd [26] performed a series of experiments with gaseous detonation loading. The schematic setting of the experiment is shown in Figure 12. The preflawed cylinder was linked to a detonation tube and an extension tube. The initial surface notch was located at the center of the cylinder. All the tubes were filled with explosive gas. The detonation source point is located inside the detonation tube, 1.52 m away from the left end of the preflawed tube.Figure 12
Preflawed cylinder with gaseous detonation loading.The pressure wave was initiated at the source point and then passed the specimen and the extension tube, causing the original surface notch to form a crack cutting through the cylinder wall and propagate.Chao [25] and Chao and Shepherd [26] launched 9 shots with the above specimen with the initial notch length Lc=25.4mm. The length of the specimen is L=0.61 m, the diameter D=0.038 m, and the shell thickness t=0.89 mm. They found two types of fracture behaviors. One was that both the forward and the backward crack tips curved after they were formed and went straight for a short distance, as shown in Figure 13(a). The other type of results also showed a backward crack curving but the forward crack tip bifurcated and finally cut the specimen into two segments. The configuration of the second fracture pattern is shown in Figure 13(b).Experimental results of Shepherd experiment [25, 26]: (a) both the forward and backward tips curved and (b) backward tip curving and forward tip bifurcation.
(a)
(b)The specimen is modeled with 40680 shell elements. The left and the right ends of the numerical models were fully clamped in the simulation. The following fitted exponential-decay curve [28] was used to represent the detonation pressure:
(28)p(x,t)={0,0<t<tcj,pcjexp(-t-tcjT),t>tcj,
where x is the distance away from the initiation point along the axial direction, t is time, pcj is the peak value of the pressure wave, tcj=x/vcj is the time for the wave to travel from the initiation point to the evaluation point, vcj is the velocity of the wave, and the time parameter T=3.0tcj. In the current simulation, pcj=6.1 Mpa and vcj=2404m/s. The pressure was applied to all the elements from the internal side. The material properties are Young’s modulus E=2780 kg/m3, Poisson ratio ν=0.3, initial yield σy0=275MPa, and hardening modulus ET=640 MPa.The configurations at different times of numerical results are shown in Figures14 and 15. The stress concentration can be seen in front of the crack tip. The crack tips went straight at the early stage of propagation. Both the forward and the backward tips curved after around 0.3 ms. The forward crack tip shows sharper curving at about 90 degrees to the axis. The backward crack tip shows a slanted path.Numerical results of Shepherd experiment at different times: deformed configurations with effective stress contour plots at (a) timet=0.255 ms and (b) time t=0.3 ms and with damage plot at (c) time t=0.42 ms and (d) time t=0.54 ms.
(a)
(b)
(c)
(d)Final configurations of Shepherd experiment. (a) Top view. (b) Side view.
(a)
(b)The final fracture patterns are shown in Figure15. The crack propagation paths are similar to the experimental results shown in Figure 13. The difference is that our results show a little shorter crack length, which may be due to the difference of the loading and boundary conditions between numerical modeling and real experiment.


## 6.1. Implosion Induced Dynamic Failure of Cylinder

The cylinder is of lengthL with two extension parts of length Le at the ends; see Figure 6. Rigid plugs are perfectly bonded to the extension parts. The thickness of the cylinder is t, and the diameter is D. In the experiments, the specimens are loaded with hydrostatic pressure and the pressure was increased until the cylinders buckled. It was observed in the experiment that specimens will buckle in some specific modes, which have a dependence on the geometry.Figure 6
Setup for implosion induced failure of cylinder.In this example, we focused on predicting final fracture pattern of two experiments. The first specimen is denoted by IMP26 experiment, whereL=366.8 mm, D=38.085 mm, and t=0.701 mm. The length of the extension parts is Le=25.4 mm. The material of the cylinder is aluminum alloy (AL) 6061-T6. In the experiment, the result indicated a mode 2 implosion, but there was no fracture observed. The second specimen, denoted by IMP25, is L=143 mm plus two 25.4 mm extension parts at the ends, D=38.087 mm, and t=0.701 mm. Because it has a much shorter length than that of IMP26, IMP25 specimen buckled in a higher mode, that is, mode 3 buckling, and showed fracture on the interfaces of the end plugs and the main specimen.For numerical analysis, we modeled IMP26 and IMP25 specimens with 49200 and 24000 shell elements, respectively. The average element size is about 1 mm. The material behavior of AL6061-T6 is modeled with the damage plasticity [23, 24] with Young’s modulus E=69.5 GPa, the density ρ=2780 kg/m3, Poisson ratio ν=0.3, yield stress σy0=276 MPa, and the hardening modulus ET=634 MPa.The following coordinate system was used: thex-axis is along the axis of the cylinder, and the y- and z-axes are in the radial direction. A uniform surface pressure loading was applied on all the shell elements. Another set of concentrated forces along the axis direction were applied to the nodes on the rim of the two ends. These forces served to compress the cylinder along its longitudinal axis. The summation of the magnitude of these concentrated forces is equivalent to the total force exerted on the end plug due to the pressure. All the nodes in the two extension parts are only allowed to move along the x-axis, and all the other DOFs of these nodes, including the translational DOFs uy and uz and the rotational DOFs θx, θy, and θz, are constrained. The employment of the applied concentrated forces and the sliding boundary conditions is to model the effect of the plug.A bilinear load curve for the pressure was used in both simulations. The pressure started at zero and was increased topcr, that is, the experimental critical buckling pressure, in time t0 and then was kept constant at pcr until the end of simulation. The schematic of time history of the pressure is shown in Figure 7; pcr=1.15MPa for IMP26 experiment and pcr=2.83MPa for IMP25 experiment. We used the parameter t0=5.0×10-5 s in both simulations.Figure 7
Loading curve used in the simulation of Texas experiments.Geometry imperfection is introduced into the radius to evoke circumferential buckling easily. The actual radius with imperfection has the formr=r0(1-acosnθ), where r0 is the unperturbed radius, a is the imperfection magnitude, and θ is the circumferential angle. n is the number of wavelengths in the circumferential direction, n=2 for IMP26 simulation, and n=3 for IMP25 simulation. It should be noted that the buckling mode in numerical results does not depend on the imperfection mode. We also tested the IMP25 simulation with injection of mode 2 imperfection and obtained very similar results.Figures8(a)–8(d) show four snapshots of the numerical results for IMP26. The center of the specimen yielded first and then buckled in mode 2. The buckling region evolved toward the two ends, and the specimen entirely collapsed. a=0.1% was used in the simulation of IMP26, which is close to 0.107%, the maximum ovalization measured in the experimental specimen.Snapshots of numerical results of IMP26 experiment: (a) initial configuration and the mesh, (b) deformed configuration at timet=1.48 ms; the center of the cylinder begins to collapse, (c) deformed configuration at time t=1.78 ms; buckling region enlarges toward the two ends, and (d) final deformed configuration; the cylinder buckles in mode 2.
(a)
(b)
(c)
(d)In IMP25 simulation, we allowed an injection of the discontinuity near the interface of the main part and the extension part of the specimen. We also observed large plastic strain, large damaged or unstable material points along the central buckling lines. However, this may be due to the repulsive forces generated during the contact, so no crack was allowed to initiate in these regions.Several imperfection magnitudesa from 0.05% to 1% were tested in IMP25 simulation. The maximum oval imperfection measured in the experiment is 0.043%. It should be noted there should not be imperfection as large as 1% in the actual specimen, and this is only to examine the effect of the imperfection.Figure9 shows four snapshots of IMP25 results of a=0.05% at different times. The central part of the cylinder collapses first, and then the collapse region enlarges toward the ends. Some cracks were initiated at the ends and propagated along the circumferential direction.Snapshots of numerical results of IMP25 experiment: (a) deformed configuration at timet=0.92 ms; cylinder collapses at mode 3, (b) deformed configuration at time t=1.08 ms; cracks initiated at the two ends, (c) the final deformed configuration, (d) local zoom view of the crack, and (e) view from the axial direction.
(a)
(b)
(c)
(d)
(e)The final configurations of different imperfection magnitudes are compared in Figure10. No large effect can be observed for imperfection magnitudes though the crack opening of a=1.0% looks larger than that of a=0.05%. The time histories of velocity y of finite element node 11670 (at the center of cylinder) are plotted and compared in Figure 11. Let the collapse time be defined as the time to reach the peak velocity. It can be observed that smaller imperfection leads to longer collapse time, but the peak velocity of all the imperfections does not show an obvious dependence on the imperfection magnitude.Comparison of final configurations of cylinders with different imperfection magnitudes: (a) 0.05% imperfection, (b) 0.1% imperfection, (c) 0.5% imperfection, and (d) 1.0% imperfection.
(a)
(b)
(c)
(d)Figure 11
Comparisons of time history of nodal velocity at the center node with different amplitudes of imperfection.


## 6.2. Explosion Induced Failure of Cylinder

Chao [25] and Chao and Shepherd [26] performed a series of experiments with gaseous detonation loading. The schematic setting of the experiment is shown in Figure 12. The preflawed cylinder was linked to a detonation tube and an extension tube. The initial surface notch was located at the center of the cylinder. All the tubes were filled with explosive gas. The detonation source point is located inside the detonation tube, 1.52 m away from the left end of the preflawed tube.Figure 12
Preflawed cylinder with gaseous detonation loading.The pressure wave was initiated at the source point and then passed the specimen and the extension tube, causing the original surface notch to form a crack cutting through the cylinder wall and propagate.Chao [25] and Chao and Shepherd [26] launched 9 shots with the above specimen with the initial notch length Lc=25.4mm. The length of the specimen is L=0.61 m, the diameter D=0.038 m, and the shell thickness t=0.89 mm. They found two types of fracture behaviors. One was that both the forward and the backward crack tips curved after they were formed and went straight for a short distance, as shown in Figure 13(a). The other type of results also showed a backward crack curving but the forward crack tip bifurcated and finally cut the specimen into two segments. The configuration of the second fracture pattern is shown in Figure 13(b).Experimental results of Shepherd experiment [25, 26]: (a) both the forward and backward tips curved and (b) backward tip curving and forward tip bifurcation.
(a)
(b)The specimen is modeled with 40680 shell elements. The left and the right ends of the numerical models were fully clamped in the simulation. The following fitted exponential-decay curve [28] was used to represent the detonation pressure:
(28)p(x,t)={0,0<t<tcj,pcjexp(-t-tcjT),t>tcj,
where x is the distance away from the initiation point along the axial direction, t is time, pcj is the peak value of the pressure wave, tcj=x/vcj is the time for the wave to travel from the initiation point to the evaluation point, vcj is the velocity of the wave, and the time parameter T=3.0tcj. In the current simulation, pcj=6.1 Mpa and vcj=2404m/s. The pressure was applied to all the elements from the internal side. The material properties are Young’s modulus E=2780 kg/m3, Poisson ratio ν=0.3, initial yield σy0=275MPa, and hardening modulus ET=640 MPa.The configurations at different times of numerical results are shown in Figures14 and 15. The stress concentration can be seen in front of the crack tip. The crack tips went straight at the early stage of propagation. Both the forward and the backward tips curved after around 0.3 ms. The forward crack tip shows sharper curving at about 90 degrees to the axis. The backward crack tip shows a slanted path.Numerical results of Shepherd experiment at different times: deformed configurations with effective stress contour plots at (a) timet=0.255 ms and (b) time t=0.3 ms and with damage plot at (c) time t=0.42 ms and (d) time t=0.54 ms.
(a)
(b)
(c)
(d)Final configurations of Shepherd experiment. (a) Top view. (b) Side view.
(a)
(b)The final fracture patterns are shown in Figure15. The crack propagation paths are similar to the experimental results shown in Figure 13. The difference is that our results show a little shorter crack length, which may be due to the difference of the loading and boundary conditions between numerical modeling and real experiment.


## 7. Conclusion

We described a new finite element method for prediction of dynamic fractures in thin shells. The method is incorporated within an explicit time integration scheme and able to represent the crack paths free from initial mesh topologies. For the representation of discontinuities due to cracks, the described method employs a simplified version of the conventional XFEM based on the phantom node method. In this approach, the cracked shell element is treated by two superimposed elements with newly added phantom nodes on the cracked portions.The method is implemented for the DKT shell element. This facilitates the implementation of the method into standard finite element programs. Another attractive feature of the method is that it provides an easy mesh generation and a relatively low computational cost and this allows large scale nonlinear dynamic fracture problems to be solved efficiently.



---

*Source: 957286-2013-10-09.xml*