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.
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 Figure 1.
(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, is the rotation whose displacement is along the axis .
(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 and directions as where and . 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 As shown in Figure 3, the discrete Kirchhoff constraints along each side are introduced at the midpoint : where is the absciss along the side of the element.
The different functions appearing in (2) and (3) are given by Note that, in (2) and (3), and 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 , , , , , and . 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 is the same in terms of substituting (2) and (4). Furthermore, for this particular edge that has the midpoint at , we could further develop the following relation: So, the expression of is given as Similarly, we can determine the other unknowns and as However, the other three unknowns , , and in the formalism of the rotation should be determined by solving system equations given by Throughout this procedure, that is, determining three unknowns , , and , 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 which are zero at the nodes.
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 where and are the first and second halves of the pair of overlapping element, respectively, 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 , and is the step function given by 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 and can be expressed as where again an overlapping pair (elements and ) 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 where ( and , resp.) denotes the displacement (velocity and acceleration, resp.) at time . is the time integration step, is the mass matrix, and (, resp.) is the external (internal, resp.) forces at time .
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 ; 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 and stiffness matrix since where is the maximum frequency determined by solving eigenvalue of problems in (16): 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 of the mass corresponding to the enriched DOFs also depends on the enrichment function as follows: where is the total actual mass of the element, is the number of nodes of the element, and is the measure of the finite element domain . 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: where is the volume or area of the regular element and 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: where is damage parameter, is the plastic strain, is a reference strain envelope, and is a material constant. The reference strain envelope is a function of the pressure and the Lode angle: where is the initial reference strain, is the pressure, and is the Lode angle. The functions and are defined as where and are material constants and and are parameters determining the shape of the strain envelope. The weakening effect caused by the damage was also considered in this model: where is the stress of the undamaged material, 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 . For the computation of the averaged strain, , we used a pointwise weighted averaging scheme which is given by where is the cubic spline weight function, is the distance from the crack tip to the material points , is the size of the averaging domain, and is the size of the crack tip element; see Figure 5. Note that, for the computation of the averaged strain , we only included material points which show tension dominant states as shown in Figure 5(b).
(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 length with two extension parts of length at the ends; see Figure 6. Rigid plugs are perfectly bonded to the extension parts. The thickness of the cylinder is , and the diameter is . 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.
In this example, we focused on predicting final fracture pattern of two experiments. The first specimen is denoted by IMP26 experiment, where mm, mm, and mm. The length of the extension parts is 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 mm plus two 25.4 mm extension parts at the ends, mm, and 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 GPa, the density kg/m3, Poisson ratio , yield stress MPa, and the hardening modulus MPa.
The following coordinate system was used: the -axis is along the axis of the cylinder, and the - and -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 -axis, and all the other DOFs of these nodes, including the translational DOFs and and the rotational DOFs , , and , 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 to , that is, the experimental critical buckling pressure, in time and then was kept constant at until the end of simulation. The schematic of time history of the pressure is shown in Figure 7; MPa for IMP26 experiment and MPa for IMP25 experiment. We used the parameter s in both simulations.
Geometry imperfection is introduced into the radius to evoke circumferential buckling easily. The actual radius with imperfection has the form , where is the unperturbed radius, is the imperfection magnitude, and is the circumferential angle. is the number of wavelengths in the circumferential direction, for IMP26 simulation, and 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.
Figures 8(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. % was used in the simulation of IMP26, which is close to 0.107%, the maximum ovalization measured in the experimental specimen.
(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 magnitudes 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.
Figure 9 shows four snapshots of IMP25 results of % 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.
(a)
(b)
(c)
(d)
(e)
The final configurations of different imperfection magnitudes are compared in Figure 10. No large effect can be observed for imperfection magnitudes though the crack opening of % looks larger than that of %. The time histories of velocity 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.
(a)
(b)
(c)
(d)
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.
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 mm. The length of the specimen is m, the diameter m, and the shell thickness 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).
(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: where is the distance away from the initiation point along the axial direction, is time, is the peak value of the pressure wave, is the time for the wave to travel from the initiation point to the evaluation point, is the velocity of the wave, and the time parameter . In the current simulation, Mpa and m/s. The pressure was applied to all the elements from the internal side. The material properties are Young’s modulus kg/m3, Poisson ratio , initial yield MPa, and hardening modulus MPa.
The configurations at different times of numerical results are shown in Figures 14 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.
(a)
(b)
(c)
(d)
(a)
(b)
The final fracture patterns are shown in Figure 15. 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The support of the Office of Naval Research under Grants nos. N00014-13-1-0386 and N00014-11-1-0925 is gratefully acknowledged.