# Explicit Dynamic Finite Element Method for Failure with Smooth Fracture Energy Dissipations

Authors: Jeong-Hoon Song; Thomas Menouillard; Alireza Tabarraei
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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.


---

## Body

## 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 Section2. 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(1)∂Pji∂Xj+ρ0bi-ρ0u¨i=0inΩ0,
where P is the nominal stress tensor, ρ0 is the initial mass density, and b is the body force vector. The boundary conditions are
(2)nj0Pji=t¯i0onΓt0(3)ui=u¯ionΓu0(4)nj0Pji-=-nj0Pji+=τi0c(⟦ui⟧)onΓc0,
where n0 is the unit normal vector to the boundary, τ0c is the cohesive traction across the crack surfaces, t¯0 is the applied traction on the Neumann boundary Γt, and u¯ is the applied displacement on the Dirichlet boundary Γu0; Γu0∪Γt0=Γ0, Γu∩Γt=∅. Superscript plus and minus signs refer to the two sides of the discontinuity.The spaces of admissible function are(5)𝒰={u(X,t)∣u(X,t)∈C0,u(X,t)=u¯(t)onΓu0,udiscontinuousonΓc},(6)𝒰0={δu(X,t)∣δu(X,t)∈C0,δu(X,t)=0onΓu0,δudiscontinuousonΓc}.The weak form of the momentum equation is given by foru(X,t)∈𝒰(7)δWkin=δWint-δWext+δWcoh∀δu(X)∈𝒰0,
where δWint is the internal work, δWext is the external work performed by the applied loads, δWkin is the kinetic work performed by inertia forces, and δWcoh is the work performed by the cohesive traction on the crack surface Γc. These quantities are defined as
(8)δWkin=∫Ω0δu·ρ0u¨dΩ0,δWint=∫Ω0∂δu∂X:PdΩ0,δWext=∫Ω0δu·ρ0bdΩ0+∫Γt0δu·t¯0dΓt0,(9)δWcoh=-∫Γcδ⟦u⟧·τcdΓc,
where t¯ is the normalized traction prescribed on Γt0 and τc 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 atX=a, as shown in Figure 1(a).Representation of a discontinuity in a one-dimensional finite element: (a) problem description of a one-dimensional rod and (b) representation of a discontinuity with the phantom node method; solid circles denote real nodes and hollow circles denote phantom nodes.
(a)
(b)The displacement field in cracked element 1 can be seen to consist of two separated displacement fields as shown in Figure1(b): the displacement field of element 1a for X<a, and element 1b for X>a. 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:(10)Iisaphantomnodein{element1aifXI-a>0element1bifXI-a<0.We can now rewrite the displacement field of element 1 as a set of two superimposed elements with phantom nodes(11)u1(X,t)=u1a(X,t)+u1b(X,t)=[u11aN1+p21aN2]H(-(X-a))︸u1a+[p11bN1+u21bN2]H(X-a)︸u1b,
where a superscript and subscript denote the element and node number, respectively, NI are the shape functions, p is the nodal unknown of the phantom node, and H(x) is the Heaviside step function defined by
(12)H(x)={1x≥00x<0.The displacement jump across the crack is given by(13)⟦u⟧X=a=limɛ→0[u1a(X+ɛ)-u1b(X-ɛ)]X=a={u11aN1(a)+p21aN2(a)}-{p11bN1(a)+u21bN2(a)}=(uI1a-uI1b)NI(a).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Ω0 as shown in Figure 2. The motion is described by x=Φ(X,t), where X and x denote material and spatial coordinates, respectively. In the current domain, the image of the initial domain Ω0 is denoted by Ω. We allow this domain to contain internal discontinuities Γc which is enveloped by a region Ωc.Figure 2
Representation of a two-dimensional body with an internal discontinuity in the initial and the current domains.Inside of the regionΩc, we define two local level set functions f(X) and g(X,t), where f(X) and g(X,t) are signed distance functions which describe the crack surface and tip geometry, respectively. The isozero line of the function f(X), that is, f(X)=0, corresponds to the crack surface Γc, and the function g(X,t) is defined so that g(X,t)>0 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
(14)X∈Γc0iff(X)=0g(X,t)>0,X∈Ωc0.
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
(15)∑IfINI(X)=0,
where fI=f(X) and f(X)=minX¯∈Γc∥X-X¯∥. As a consequence of (15), the surface of discontinuity can be represented by f(XI) 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 g(X,t) by a list of cracked elements.Figure 3
Two-dimensional discontinuity representation by two implicit functionsf(X) and g(X,t).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 Figure4.Figure 4
Representation of a cracked element in two dimensions; solid circles denote real nodes and hollow circles denote phantom nodes.The displacement field of this superimposed element is(16)u(X,t)=∑I=1nNuI1a(t)NI(X)H(-f(X))︸u1a(X,t)+uI1b(t)NI(X)H(f(X))︸u1b(X,t).The explicit value of the displacement jump is given by(17)⟦u⟧=limɛ→0[u1a(X+ɛ∇f)-u1b(X-ɛ∇f)]=[NIH(-f(X))uI1a-NIH(f(X))uI1b]X∈Γc=(uI1a-uI1b)NI(a).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 Figure5, the element can be replaced with three overlapping elements and the discontinuous displacement fields can be represented by
(18)u(X,t)=u1a(X,t)+u1b(X,t)+u1c(X,t)=∑I∈S1uINIH(-f1(X))H(-f2(X))+∑I∈S2uINIH(-f1(X))H(f2(X))+∑I∈S3uINIH(f1(X))H(-f2(X)),
where f1(X) and f2(X) are level set functions for crack 1 and 2, respectively.Figure 5
Decomposition of an element into three elements1a, 1b, and 1c in order to model crack branching; solid and hollow circles denote the original nodes and the added phantom nodes, 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 Figure5, the phantom node method cannot resolve this case and in this case, minimal remeshing is required to properly model crack branching.


## 3.1. Phantom Node Method in Two Dimensions

Consider an initial domainΩ0 as shown in Figure 2. The motion is described by x=Φ(X,t), where X and x denote material and spatial coordinates, respectively. In the current domain, the image of the initial domain Ω0 is denoted by Ω. We allow this domain to contain internal discontinuities Γc which is enveloped by a region Ωc.Figure 2
Representation of a two-dimensional body with an internal discontinuity in the initial and the current domains.Inside of the regionΩc, we define two local level set functions f(X) and g(X,t), where f(X) and g(X,t) are signed distance functions which describe the crack surface and tip geometry, respectively. The isozero line of the function f(X), that is, f(X)=0, corresponds to the crack surface Γc, and the function g(X,t) is defined so that g(X,t)>0 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
(14)X∈Γc0iff(X)=0g(X,t)>0,X∈Ωc0.
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
(15)∑IfINI(X)=0,
where fI=f(X) and f(X)=minX¯∈Γc∥X-X¯∥. As a consequence of (15), the surface of discontinuity can be represented by f(XI) 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 g(X,t) by a list of cracked elements.Figure 3
Two-dimensional discontinuity representation by two implicit functionsf(X) and g(X,t).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 Figure4.Figure 4
Representation of a cracked element in two dimensions; solid circles denote real nodes and hollow circles denote phantom nodes.The displacement field of this superimposed element is(16)u(X,t)=∑I=1nNuI1a(t)NI(X)H(-f(X))︸u1a(X,t)+uI1b(t)NI(X)H(f(X))︸u1b(X,t).The explicit value of the displacement jump is given by(17)⟦u⟧=limɛ→0[u1a(X+ɛ∇f)-u1b(X-ɛ∇f)]=[NIH(-f(X))uI1a-NIH(f(X))uI1b]X∈Γc=(uI1a-uI1b)NI(a).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 Figure5, the element can be replaced with three overlapping elements and the discontinuous displacement fields can be represented by
(18)u(X,t)=u1a(X,t)+u1b(X,t)+u1c(X,t)=∑I∈S1uINIH(-f1(X))H(-f2(X))+∑I∈S2uINIH(-f1(X))H(f2(X))+∑I∈S3uINIH(f1(X))H(-f2(X)),
where f1(X) and f2(X) are level set functions for crack 1 and 2, respectively.Figure 5
Decomposition of an element into three elements1a, 1b, and 1c in order to model crack branching; solid and hollow circles denote the original nodes and the added phantom nodes, 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 Figure5, 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.Figure6(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, faint 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
(19)fkin=fint-fext+fcoh-fcorrect,
where ∥fcorrect∥ 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 fcoh 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 a in Figure 6. The flowchart for numerical computation procedures is described in Table 1.Table 1
Flowchart for numerical computation procedures.
Steps
Numerical computation procedures
1
Apply initial conditions including initial boundary conditions and prestress
2
Compute initial forces:f0=f0ext-f0int
3
Compute acceleration:an=Mn-1fn
4
Time update:tn+1=tn+Δtn+(1/2)
5
First partial update for nodal velocities:vn+(1/2)=vn+(tn+(1/2)-tn)an
6
Apply velocity boundary conditions attn+(1/2)
7
Update nodal displacements:dn+1=dn+Δtn+(1/2)×vn+(1/2)
8
Compute strain and stress
9
Check fracture criterion and insert phantom nodes
10
Compute forces at time:fn+1=fn+1ext-fn+1int+fn+1coh+fn+1correct
11
Compute acceleration:an+1=Mn+1-1fn+1
12
Second partial update for nodal velocities:vn+1=vn+(1/2)+(tn+1-tn+(1/2))an+1
13
Update countern
14
Output: if simulation is not completed, go to Step 4Evolution from the state (a) to (d) through the intermediate states (b) and (c) where the correction force decreases.Mu denotes the mass matrix related to the standard degrees of freedom u, Mb to the existing discontinuous degree of freedoms b, and Ma to the new discontinuous degrees of freedom a. fext denotes the external force, fcoh denotes the cohesive force, and fint denotes the internal forces related to the different kinds of degrees of freedom (i.e. u, b, and a, respectively, denoted by a cross, a continuous circle, and a dotted circle).
(a)
(b)
(c)
(d)The evolution of the correction force is shown in Figure7 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
(20)fecorrect(t)=(1-a(t)le)feint(tinj),
where tinj is the time the corresponding degrees of freedom are injected, a(t) is the crack propagation speed, and le is the characteristic length of finite elements.Figure 7
Law of the correction force as a function of the fracture ratio of the crack tip element (i.e., the position of the crack tip in the element).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 t≤3tc=3h/c1 (where c1 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 L=10 m, the initial crack length a=5 m, and the vertical position of the crack is h=2 m. Two regular meshes are used: 78×39 and 120×59 4-node elements. The material properties are Young’s modulus E=210 GPa, Poisson’s ratio ν=0.3, and density ρ=8,000 kg/m3. The tensile stress applied on the top surface is σ0=500 MPa. The crack velocity is imposed to be zero until 1.5tc and 1,500m/s after. The mode 1 stress intensity factor is normalized by the factor σ0h.Figure 8
Geometry and loading of the semi-infinite plate example.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 factorK1 and the velocity a˙ of the crack is given by [23]
(21)K1(a˙,t)={0ift<tc2σ01-νc1(t-tc)(1-2ν)πiftc≤t<1.5tc2σ01-νc1(t-tc)(1-2ν)π×1-(a˙/cr)1-(a˙/2cr)if1.5tc≤t,
where the Rayleigh wave speed is cr=2,947 m/s and the dilatational wave speed is c1=5,944 m/s.Figure9 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.Stress field in the fine mesh: (a) without correction and (b) with correction force.
(a)
(b)Figures10(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.Normalized stress intensity factor as a function of time computed on the (a) coarse mesh and (b) fine mesh; relative error on stress intensity factor as a function of time on the (c) coarse mesh and (d) fine mesh.
(a)
(b)
(c)
(d)


### 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 J2 plasticity theory is used to model the behavior of the specimen. The material properties are Young’s modulus E=69 GPa, Poisson’s ratio ν=0.3, density ρ=2,780 kg/m3, yield stress σ0=275 MPa, and hardening slope h=640 MPa. A constant velocity of 1.8 m/s is applied at the center of two steel bars located in each hole of the specimen.Figure 11
Compact tension specimen geometry (W=0.05m) and the thickness is W/2.(a) Deformed mesh and (b) experimental postmortem specimen [24].
(a)
(b)Figure12(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.Figure 13
Comparison of the load deflection curve between the experimental data [24] and computational results.


### 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 [26–30]. 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 a0.1 m by 0.04 m prenotched specimen as shown in Figure 14. Tensile traction, σ=1 MPa, is applied on both the top and bottom edges as a step function in time. We discretized the domain with 100×51 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.Figure 14
Plate with an initial notch under tensile stress.The pattern of multiple crack branchings with damage evolution is shown in Figure15. The numerical simulation is executed until one of the crack tips reaches the boundary of the specimen: t=56.02 μ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].The pattern of crack branching and damage evolution at different time steps: (a)t=30.14 μs, (b) t=45.24 μs, and (c) t=56.02 μs.
(a)
(b)
(c)
(d)


### 5.4. Thick Cylinder under Internal Pressure

We consider a thick cylinder under high internal pressure. The inner and outer radii of the cylinder are80mm and 150 mm, respectively. An internal pressure, p=p0e-t/t0, is applied with p0=10 GPa and t0=0.1 ms. The material properties are ρ=7800 kg/m3, E=210 GPa, and ν=0.30. Also, we induced a ±5% 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 first5.0 μ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.Deformed shape of the thick cylinder under internal pressure at different time steps: (a)t=1.78 μs, (b) t=3.42 μs, (c) t=6.40 μs, and (d) t=13.50 μs; the long elements are cracked as indicated in (d).
(a)
(b)
(c)
(d)Shape of the fragments of the thick cylinder under internal pressure at different time steps: (a)t=3.72 μs and (b) t=6.71 μs.
(a)
(b)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.


## 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 t≤3tc=3h/c1 (where c1 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 L=10 m, the initial crack length a=5 m, and the vertical position of the crack is h=2 m. Two regular meshes are used: 78×39 and 120×59 4-node elements. The material properties are Young’s modulus E=210 GPa, Poisson’s ratio ν=0.3, and density ρ=8,000 kg/m3. The tensile stress applied on the top surface is σ0=500 MPa. The crack velocity is imposed to be zero until 1.5tc and 1,500m/s after. The mode 1 stress intensity factor is normalized by the factor σ0h.Figure 8
Geometry and loading of the semi-infinite plate example.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 factorK1 and the velocity a˙ of the crack is given by [23]
(21)K1(a˙,t)={0ift<tc2σ01-νc1(t-tc)(1-2ν)πiftc≤t<1.5tc2σ01-νc1(t-tc)(1-2ν)π×1-(a˙/cr)1-(a˙/2cr)if1.5tc≤t,
where the Rayleigh wave speed is cr=2,947 m/s and the dilatational wave speed is c1=5,944 m/s.Figure9 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.Stress field in the fine mesh: (a) without correction and (b) with correction force.
(a)
(b)Figures10(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.Normalized stress intensity factor as a function of time computed on the (a) coarse mesh and (b) fine mesh; relative error on stress intensity factor as a function of time on the (c) coarse mesh and (d) fine mesh.
(a)
(b)
(c)
(d)


## 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 J2 plasticity theory is used to model the behavior of the specimen. The material properties are Young’s modulus E=69 GPa, Poisson’s ratio ν=0.3, density ρ=2,780 kg/m3, yield stress σ0=275 MPa, and hardening slope h=640 MPa. A constant velocity of 1.8 m/s is applied at the center of two steel bars located in each hole of the specimen.Figure 11
Compact tension specimen geometry (W=0.05m) and the thickness is W/2.(a) Deformed mesh and (b) experimental postmortem specimen [24].
(a)
(b)Figure12(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.Figure 13
Comparison of the load deflection curve between the experimental data [24] and computational results.


## 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 [26–30]. 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 a0.1 m by 0.04 m prenotched specimen as shown in Figure 14. Tensile traction, σ=1 MPa, is applied on both the top and bottom edges as a step function in time. We discretized the domain with 100×51 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.Figure 14
Plate with an initial notch under tensile stress.The pattern of multiple crack branchings with damage evolution is shown in Figure15. The numerical simulation is executed until one of the crack tips reaches the boundary of the specimen: t=56.02 μ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].The pattern of crack branching and damage evolution at different time steps: (a)t=30.14 μs, (b) t=45.24 μs, and (c) t=56.02 μs.
(a)
(b)
(c)
(d)


## 5.4. Thick Cylinder under Internal Pressure

We consider a thick cylinder under high internal pressure. The inner and outer radii of the cylinder are80mm and 150 mm, respectively. An internal pressure, p=p0e-t/t0, is applied with p0=10 GPa and t0=0.1 ms. The material properties are ρ=7800 kg/m3, E=210 GPa, and ν=0.30. Also, we induced a ±5% 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 first5.0 μ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.Deformed shape of the thick cylinder under internal pressure at different time steps: (a)t=1.78 μs, (b) t=3.42 μs, (c) t=6.40 μs, and (d) t=13.50 μs; the long elements are cracked as indicated in (d).
(a)
(b)
(c)
(d)Shape of the fragments of the thick cylinder under internal pressure at different time steps: (a)t=3.72 μs and (b) t=6.71 μs.
(a)
(b)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.



---

*Source: 293861-2013-11-14.xml*