A numerical method for dynamic failure analysis through the phantom node method is further developed. A distinct feature of this method is the use of the phantom nodes with a newly developed correction force scheme. Through this improved approach, fracture energy can be smoothly dissipated during dynamic failure processes without emanating noisy artifact stress waves. This method is implemented to the standard 4-node quadrilateral finite element; a single quadrature rule is employed with an hourglass control scheme in order to decrease computational cost and circumvent difficulties associated with the subdomain integration schemes for cracked elements. The effectiveness and robustness of this method are demonstrated with several numerical examples. In these examples, we showed the effectiveness of the described correction force scheme along with the applicability of this method to an interesting class of structural dynamic failure problems.

1. Introduction

It has been shown that the extended finite element method (XFEM) [1, 2] can be successfully applied to several types of internal discontinuity problems, focusing on failure problems. For example, the XFEM has been applied to arbitrary branched and intersecting cracks [3], three dimensional crack propagation [4, 5], cohesive crack models [6] and dynamic shell [7], and 2D [8] fracture problems.

Although the standard XFEM has been successfully applied to dynamic fracture problems by Belytschko et al. [8], they mostly considered fracture problems with a single initial notch and/or simple crack geometry. This limitation arises from the difficulty in the representation of a complicated crack geometry and with the numerical integration. To circumvent this difficulty, Xu and Needleman [9], Ortiz and Pandolfi [10], Repetto et al. [11], and Cirak et al. [12] proposed the interelement crack model; but this interelement crack model can have mesh sensitivity problems as pointed out in Belytschko et al. [8] and Song et al. [13].

Our motivation in presenting this paper is to illustrate a new method which is efficient for dynamic propagation of multiple cracks and fragmentation problems but nevertheless free from mesh sensitivity by using intra-element discontinuities. In this method, we use an element superposition concept to represent cracked elements which was proposed by Song et al. [14].

In addition, we describe a method to deal with the description of the crack tip element in terms of stiffness. It aims at characterizing the release of the crack tip element when the crack propagates. An additional correction force introduced by Menouillard and Belytschko [15] takes into account this crack tip element and makes the new additional degrees of freedom release continuously. This artificial correction force will be used to smooth the stress near the crack tip because no tip enrichment is used in our XFEM discretization.

Menouillard et al. [16] developed a mass lumping strategy for the XFEM formulation and more particularly for the discontinuous enrichment part. They found that the enrichment does not significantly decrease the stable time step. Then, Menouillard et al. [17] used another decomposition of the enriched shape function developed by A. Hansbo and P. Hansbo [18] which is used in the phantom node method developed by Song et al. [14] and Song and Belytschko [19].

The outline of this paper is as follows. The governing equation and its weak form are given in Section 2. The representation of a discontinuity in a cracked element with a phantom node method is presented in Section 3. Section 4 presents a new method for smoothly releasing the newly cracked element near the tip by taking into account the position of the tracked crack tip in the tip element. Several numerical examples are given in Section 5. Section 6 presents the conclusions of this paper.

2. Governing Equations and Weak Form

For a two-dimensional dynamic problem, the strong form of the linear momentum equation in a total Lagrangian description is where is the nominal stress tensor, is the initial mass density, and is the body force vector. The boundary conditions are where is the unit normal vector to the boundary, is the cohesive traction across the crack surfaces, is the applied traction on the Neumann boundary , and is the applied displacement on the Dirichlet boundary ; , . Superscript plus and minus signs refer to the two sides of the discontinuity.

The spaces of admissible function are

The weak form of the momentum equation is given by for where is the internal work, is the external work performed by the applied loads, is the kinetic work performed by inertia forces, and is the work performed by the cohesive traction on the crack surface . These quantities are defined as where is the normalized traction prescribed on and is the cohesive traction applied on the discontinuity surface; an updated Lagrangian form is used for the cohesive work in (9).

3. Representation of a Discontinuity with Phantom Nodes

We first illustrate the modeling of a one-dimensional cracked element with phantom nodes and then give the general description for the modeling of the two-dimensional case, subsequently. Consider a one-dimensional rod and let a crack be located inside of element 1 at , as shown in Figure 1(a).

The displacement field in cracked element 1 can be seen to consist of two separated displacement fields as shown in Figure 1(b): the displacement field of element 1a for , and element 1b for . To construct new elements 1a and 1b from the element 1, we add new nodes which are replicas of the original nodes; we called these nodes phantom nodes.

We define phantom nodes by the following:

We can now rewrite the displacement field of element 1 as a set of two superimposed elements with phantom nodes where a superscript and subscript denote the element and node number, respectively, are the shape functions, is the nodal unknown of the phantom node, and is the Heaviside step function defined by

The displacement jump across the crack is given by

This procedure for cracked elements is similar to the standard XFEM nodal enriching scheme. However, the phantom nodes method simplifies the implementation of cracked elements within the context of existing finite element codes, since it is only necessary to add an extra element with phantom nodes and modify the element connectivity matrices.

3.1. Phantom Node Method in Two Dimensions

Consider an initial domain as shown in Figure 2. The motion is described by , where and denote material and spatial coordinates, respectively. In the current domain, the image of the initial domain is denoted by . We allow this domain to contain internal discontinuities which is enveloped by a region .

Inside of the region , we define two local level set functions and , where and are signed distance functions which describe the crack surface and tip geometry, respectively. The isozero line of the function , that is, , corresponds to the crack surface , and the function is defined so that along the crack surface and vanishes at the crack tip; see Figure 3. By using a set of these level set functions, we can implicitly define the crack geometry by For the numerical representation, instead of employing an implicit definition of the crack surface, generally we can approximate the path of an internal discontinuity by where and . As a consequence of (15), the surface of discontinuity can be represented by at the nodes of the cracked elements [20, 21]. Note that for the element-by-element cracking scheme which is employed in this study, we can replace the function by a list of cracked elements.

For a two-dimensional element, the superposed displacement fields in the cracked element can be developed in a similar manner to the one-dimensional case. Consider cracked element 1 and replace the element with element 1a and 1b as shown in Figure 4.

The displacement field of this superimposed element is

The explicit value of the displacement jump is given by

The concept of element overlapping method can be easily extended to modeling of an arbitrary crack junction or branching problems. When the original crack 1 branches into a new crack 2 or a crack 1 junctions with a crack 2, as shown in Figure 5, the element can be replaced with three overlapping elements and the discontinuous displacement fields can be represented by where and are level set functions for crack 1 and 2, respectively.

If the crack branching angle between crack 1 and crack 2 is acute, that is, both cracks cut the element edge of nodes 2 and 5 in Figure 5, the phantom node method cannot resolve this case and in this case, minimal remeshing is required to properly model crack branching.

4. Correction Force

In this section, we aim at developing a method to deal with the release of the crack tip element when the crack propagates through. It is to avoid sudden element release near the crack tip during propagation and thus avoid unphysical stress wave propagation due to the crack propagation.

Figure 6(a) shows the crack tip element with the associated phantom nodes. When the crack tip reaches the next element, the new crack tip element is suddenly released (see the sudden passage from Figures 6(a) to 6(d)) because the corresponding internal force takes a significant value when the phantom nodes are injected (i.e, is nonzero in Figure 6(d)). Our proposed method makes a progressive release of the crack tip element. This even happens with the cohesive force. Figures 6(b) and 6(c) show the crack tip element on a dotted line and the additional correction force acting in the momentum equation, which aims at releasing smoothly the element when the crack tip travels through from one edge to the next one. Thus the modified discrete momentum equation for the newly added degrees of freedom becomes where tends to zero when the crack tip reaches the new edge, and thus the element becomes completely cut by the discontinuity (see Figure 6(d)). Note that is not displayed in Figure 6. The initial value of the correction force (when the crack tip is on the previous edge) is such that the sum of the four forces in the equation above is zero as it is shown in Figure 6(b). At this point, the correction force is the same as the internal force, and thus no acceleration occurs yet on the new additional degrees of freedom, denoted by in Figure 6. The flowchart for numerical computation procedures is described in Table 1.

The evolution of the correction force is shown in Figure 7 which describes the magnitude of the correction force as a function of the crack tip position in the tip element. Indeed the correction force goes from the initial internal force to zero when the crack tip propagates from one edge to another. Between these two crack tip positions, the correction force is taken to be linear in our simulation. However this is not a restriction to use a linear law. In other words, the correction force is only applied to the newly added degrees of freedom, for example, the additional degrees of freedom of the crack tip element (see Figure 7). The law is where is the time the corresponding degrees of freedom are injected, is the crack propagation speed, and is the characteristic length of finite elements.

The continuity of the internal force related to the new additional degrees of freedom gives the same property to the acceleration through the momentum equation. Therefore, the velocity and displacement remains quite continuous in time when additional degrees of freedom are injected, and thus the property of continuity in time remains in the strain and stress field too, and a continuous progressive release of the tip element occurs.

5. Numerical Examples

5.1. Moving Semi-Infinite Mode I Crack

The example considered in this section is an infinite plate with a semi-infinite crack [22] loaded as shown in Figure 8. A theoretical solution of this problem for a given crack velocity is given in Freund [23]. According to the geometry described in Figure 8, the analytical solution is valid until time (where is the dilatational wave speed). Beyond that, the reflected stress wave reaches the crack tip and the analytical solution is no longer valid. The dimensions of the structure are the following: the length is  m, the initial crack length  m, and the vertical position of the crack is  m. Two regular meshes are used: and 4-node elements. The material properties are Young’s modulus  GPa, Poisson’s ratio , and density  kg/m3. The tensile stress applied on the top surface is  MPa. The crack velocity is imposed to be zero until and after. The mode 1 stress intensity factor is normalized by the factor .

We study the effect of the correction force on the accuracy of the stress intensity factor of a moving crack. The analytical relation between the stress intensity factor and the velocity of the crack is given by [23] where the Rayleigh wave speed is  m/s and the dilatational wave speed is  m/s.

Figure 9 shows the stress contour at the end of the computation for the two cases, that is, with and without the correction force. One can notice that the correction force makes the stress fields smoother in the structure when the crack propagation occurs and eliminates majority of the released stress waves due to the abrupt injection of phantom node as shown in Figure 9(a). In contrary spurious stress waves appear in Figure 9(a) due to the crack propagation and the sudden release of the crack tip elements.

Figures 10(a) and 10(b), respectively, present the normalized stress intensity factor as a function of time for the coarse mesh and fine mesh, respectively, with and without the correction force. Both figures underline that the correction force improves the result during propagation by decreasing the magnitude of the oscillations due to the released crack tip element. Indeed the number of oscillations are directly related to the number of newly cracked elements. To evaluate the improvement, Figures 10(c) and 10(d) show the relative error between the computations using the correction force and not for the coarse and fine meshes. The error is decreased from 20% to 5% by adding the correction force on the newly added degrees of freedom during the crack propagation.

5.2. Stiffened Compact Tension Specimen

The stiffened compact tension specimen is used in various experiments [24, 25]. The particularity of the stiffened test is that an additional part of material opposite to the initial crack has the effect of a stiffener. With such a configuration, the crack will not be able to propagate straight toward the stiffener, but an instability will make the crack propagate up or down as a curve. Figure 11 presents the geometry of the specimen without showing the stiffened part; the stiffened part is glued on the right edge of the specimen as shown in Figure 12(b). A plasticity theory is used to model the behavior of the specimen. The material properties are Young’s modulus  GPa, Poisson’s ratio , density  kg/m3, yield stress  MPa, and hardening slope  MPa. A constant velocity of  m/s is applied at the center of two steel bars located in each hole of the specimen.

Figure 12(a) shows the final fracture pattern of the specimen at the end of the computation. The numerical result shown in Figure 12(a) agrees well with the experimental result obtained by Galanis [24] as shown in Figure 12(b). The computed load-deflection curve is also in good agreement with the experiments as shown in Figure 13.

5.3. Dynamic Multiple Crack Branchings in a Square Plate

For dynamic fracture problems, crack branching due to a dynamic instability is a common phenomenon. Several experimental results on crack branching have been previously reported [2630]. However, because of difficulties in the representation of branched crack paths, only a few numerical results have been reported [8, 9, 31] and so forth. Note that Belytschko et al. [8] allowed the original crack to branch only once; Xu and Needleman [9] used an element edge crack model which is less complex than intraelement crack models but has a certain mesh sensitivity; for the issues regarding mesh sensitivity, refer to Belytschko et al. [8]. Also, Rabczuk and Belytschko [31] discretely modeled the crack with the meshfree cracked particle method. In the following, we examine the performance of the proposed method in a crack branching problem.

We consider a  m by  m prenotched specimen as shown in Figure 14. Tensile traction,  MPa, is applied on both the top and bottom edges as a step function in time. We discretized the domain with uniform quadrilateral elements and used explicit time integration with a Courant number of 0.1. Material softening is modeled with a Lemaitre damage law [32] and a linear cohesive law was imposed once a discontinuity developed. To capture the crack branching phenomena, we monitored the maximum principal stress criterion at several additional points around the crack tip. If the Lemaitre damage criterion is satisfied and the maximum principal directions show relatively different crack growth angles, we initiate crack branches.

The pattern of multiple crack branchings with damage evolution is shown in Figure 15. The numerical simulation is executed until one of the crack tips reaches the boundary of the specimen: μs. The maximum crack tip speed is around 1250 m/s and the Rayleigh wave speed is 2100 m/s. The overall dynamic crack branching pattern agrees with the results which were already reported by Rabczuk and Belytschko [31].

5.4. Thick Cylinder under Internal Pressure

We consider a thick cylinder under high internal pressure. The inner and outer radii of the cylinder are and  mm, respectively. An internal pressure, , is applied with  GPa and  ms. The material properties are  kg/m3,  GPa, and . Also, we induced a perturbation to the elastic modulus to introduce some asymmetry for initiating the cracks. We modeled the thick cylinder with 20,000 uniform quadrilateral elements.

Because of the high internal pressure, the fragmentation process occurs only in the first μs; then each fragments moves outward in the radial direction with no further cracks initiating. Figure 16 shows the magnified deformed mesh at different time steps. In this simulation, we obtained 6 relatively big fragments and 10 strip shape fragments; for a clear illustration of the fragments, see Figure 17.

This overall pattern of cylinder fragmentation is similar to that already reported by Rabczuk and Belytschko [31]. However, the finite element simulation shown here does not exhibit the small fragments seen in [31]. It is also found that the reduced 4-node quadrilateral finite element with the hourglass control scheme is quite sensitive to mesh distortion; it would be desirable to use the smoothed finite elements [33] or the method proposed by Areias and Rabczuk [34] for this type of simulations, that is, problems with severe mesh distortion.

6. Conclusion

A numerical method for the simulation of the dynamic propagation of multiple cracks is presented. The method employs the phantom node method with a one point integration scheme. Though the phantom node method is another form of the standard XFEM, it provides us with a simple implementation within the framework of the standard FEM. Also, by using one point integration with hourglass control, we can decrease the computational cost and circumvent the subdomain integration which is generally used for cracked elements. Moreover a correction force can handle the progressive opening of the crack tip element due to the crack propagation, and this improves the results in terms of stress intensity factors.

To evaluate the applicability of the proposed method, several numerical examples which have a certain complexity in the representation of crack geometry have been analyzed. The smooth crack propagation obtained with a correction force is first checked from an energy point of view. Then the accuracy of the method is shown by computating the stress intensity factors for a dynamic mode 1 crack propagation. The simulation of a compact tension specimen with a stiffener is in good agreement with the experiment observations. We simulated a dynamic multiple crack branching problem and found that the method is particularly successful in this type of simulations. In the simulation of a multiple crack propagation problem, the numerical example shows that the proposed model can simulate the growth, interconnection and, finally, failure of a plate containing multiple cracks. Also, the simulation of a fragment processes is quite effectively analyzed. An attractive feature of this method is its low computational cost and simplicity within the context of the conventional explicit finite element method.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The support of the Office of Naval Research under Grants N00014-13-1-0386 and N00014-11-1-0925 are gratefully acknowledged.