# Stable Cracking Particles Method Based on Stabilized Nodal Integration and Updated Lagrangian Kernel

Authors: S. Xu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A stable cracking particles method (CPM) based on updated Lagrangian kernels is proposed. The idea of CPM is to model the crack topology by a set of cracked particles. Hence no representation of the crack surface is needed making the method useful for problems involving complex fracture patterns as they occur in dynamics and under fast loading conditions. For computational efficiency, nodal integration is exploited in the present paper. In order to avoid instabilities, a scheme is presented to stabilized the integration. Moreover, a set of simple cracking rules are proposed in order to prevent numerical fracture. The method is applied to two benchmark problems and shows good accuracy.


---

## Body

## 1. Introduction

Meshless methods have been a competitor to finite element method due to their ability to add particles and model large deformations, dynamic fracture, and fragmentation with ease [1–14]. In contrast to finite element methods that require the deletion of elements for complex dynamic fracture problems involving penetration and perforation, meshless methods [14–20] treat such problems quite naturally. However, meshless methods are computationally expensive. Moreover, numerical fracture is reported when fracture and material modelling is not accounted for carefully [21–23]. Recently developed efficient methods such as finite element methods with edge rotations [24–27], phase-field models [28, 29], and partition-of-unity enriched finite element [30–44] and meshless methods have been widely applied to fracture problems with a few number of cracks [23, 45–52]. However, their application to complex dynamic fracture and fragmentation remains a major challenge. A powerful method for complex fracture is the cracking particles method (CPM) [53, 54] that is based on enriching nodes with step-enrichment function once a fracture criterion is met. The crack’s topology is modeled by discrete plane crack segments and the well-posedness of the initial boundary value problem is restored by means of cohesive zone models. The advantage of the CPM is its robustness and efficiency. Complex fracture is modeled naturally with relatively coarse discretizations as the process zone does not need to be smeared over several particles. The CPM has been applied to numerous challenging problems [55–59]. However, due to the “discontinuous” representation of the crack surface, stresses might be transferred over the opening crack leading to spurious cracking [60–63].In this paper, we present a stabilized nodal integrated CPM based on updated Lagrangian kernels that alleviates several shortcomings of the original CPM. In particular, the following occurs.(i)
We propose a set of simple cracking rules according to [64] in order to avoid spurious cracking in the CPM. The original CPM shows suffering from spurious cracking adjacent to the crack surface that might artificially increase the dissipated energy.(ii)
We combine for the first time the stabilized nodal integration and the CPM. CPM was only used in combination with computationally expensive Gauss quadrature, stress-point integration and nodal integration. Nodal integration is computationally very efficient but suffers from instabilities. Stress-point integration eliminates these instabilities but cracks are restricted to cross nodes—not stress points. Moreover, there is still the need to update the position of the stress points. Therefore, we employ the stabilized conforming nodal integration [65] also used in the smoothed finite element method [66–70].(iii)
While the original CPM employs a Lagrangian kernel that is only applicable to moderate deformations, an updated Lagrangian kernel formulation is proposed here for the first time. It guarantees the applicability of the method to extremely large deformation.The paper is organized as follows. First, we present the CPM. Then, the weak form is stated and the discrete equations are derived. Subsequently, the fracture model and the cohesive zone models are discussed before the paper ends with examples and conclusions.


## 2. CPM

The key idea of the CPM is to decompose the displacement field into two parts: the continuous partuC, sometimes also referred to as the “usual” part, and the discontinuous or “enriched” part uD:
(1)u(X,t)=uC(X,t)+uD(X,t).
The discretization of the continuous part of the displacement field is based on moving least squares (MLS) shape functions [71] of linear completeness:
(2)uC(X,t)=∑I∈𝒩‍NI(X)uI(t),
where NI(X) are the meshless shape functions of node I at position X and uI(t) are nodal parameters at time t. Note that (1) the nodal parameters are not the true physical displacement values at node I and (2) the shape functions are expressed in terms of material coordinates X. It can be shown that the shape functions are given by
(3)NI(X)=pT(X)A-1(X)DI(X)
with
(4)DI(X)=w(X-XI,h)pT(XI),AI(X)=∑I∈𝒩‍w(X-XI,h)p(XI)pT(XI),
where w(X-XI,h) denotes the kernel function and h its support size; p(X)=(1,X,Y) is a linear polynomial basis. We point out again that the kernel function is expressed in terms of material coordinates and therefore is called Lagrangian kernel [22]. It was shown by [22] that a Lagrangian kernel avoids numerical fracture often observed for simulations based on Eulerian kernels. However, the Lagrangian kernel formulation limits the amount of large deformations. Therefore, we update the kernel functions every nth time step referring to a new reference configuration and hence call this kernel an updated Lagrangian kernel. A similar approach was mentioned in [54].The discontinuous part of the displacement field is obtained by simply multiplying the shape function with enrichment functions accounting for the jump in the displacement field (see Figure2):
(5)u(X,t)=∑I∈𝒩c‍NI(X)ℋΓc(X)qI(t).𝒩c being the set of cracking particles qI(t) are additional degrees of freedom (DOF) and the enrichment function is the step function
(6)ℋ(X)=X·nΓc∥X·nΓc∥.
The key strength of the CPM is that it does not require any representation of the crack's topology. In the CPM, the crack topology is described as a set of crack segments as illustrated in Figure 1 allowing the simulation of very complex crack patterns with ease.(a) Continuous crack and (b) representation of the crack with discrete cohesive crack segments.
(a)
Continuous crack
(b)
Discretization with crack segmentsFigure 2
Normal of the crack segment.It is advantageous to write the approximation of the displacement field in vector matrix notation:(7)d=ND
and store the “standard” and enriched shape functions in the matrix N; the DOFs of the “standard” and enriched part of the displacement field are stored in the vector D. Moreover, we can define the matrix [[N]] and the matrix containing the derivatives of the meshless shape functions B needed later.Though the CPM can handle complex crack patterns, it might lead to spurious cracking adjacent to the crack that can be avoided by a set of simple rules as shown in [64]. We propose similar rules and firstly distinguish between propagation and initiating cracks by a simple criterion based on a circular support domain. When a cracking particle does not find another cracking particle within the circular domain of size αh, then a new crack is initiated. The factor α is chosen to be 1.1 in all simulations. Secondly, we define a zone adjacent to a crack surface where no new cracks can evolve; see Figure 3.Figure 3
(a) Spurious cracking during crack propagation and (b) crack propagation without spurious cracking.As such a criterion might prevent crack branching, these rules need to be adjusted at the crack tip. Since the crack surface is not continuous in the CPM, we employ the following simple algorithm to detect the crack tip for branching cracks: crack branching occurs when the angle of existing and newly created particles at the crack tip exceeds a prescribed tolerance. According to Figure4, the deviation in this angle is calculated by
(8)γmax⁡=max⁡︸I∈𝒩~c(nct·nI)∀I∈𝒩~c,γmin⁡=min⁡︸I∈𝒩~c(nct·nI)∀I∈𝒩~c𝒩~c⊂𝒩c being the set of new cracking particles. Crack branching is assumed when γmax⁡-γmin⁡≥TOL.Figure 4
Crack with crack tip node (solid line) at timen and nodes that meet cracking criterion at time n+1 (dashed line).


## 3. Weak Form and Discretization

The linear momentum equation is(9)∇·P+ϱb=ϱu¨,X∈Ω,
where P is the nominal stress tensor, ϱ is the density, b are body forces, and the superimposed dots denote material time derivatives. The displacement and traction boundary conditions are
(10)u=u¯,X∈Γunt·P=t¯,X∈Γtnc·P=tc([[u]]),X∈Γc,
where the index c refers to the crack, the index t refers to traction boundaries, and the index u refers to displacement boundaries.We solve the equation of motion in weak form that can be stated in variational form. Find the displacement fieldu∈𝒰∀δu∈𝒱 such that the first variation in the energy is zero:
(11)δW=δWint⁡-δWext+δWinertia-δWcoh=0∀δu∈𝒰0
with
(12)δWint⁡=∫Ω‍∇δu:PdΩ,δWext=∫Γt‍δu·t¯dΓ+∫Ω‍ϱδu·bdΩ,δWinertia=∫Ω‍ϱδu·u¨dΩ,δWcoh=∫Γc‍δ[[u]]·tcdΓ.
The approximation spaces 𝒰 and 𝒰0 are given by
(13)𝒰={u(X,t)∣u∈H1,u=u¯on⁡Γu,udiscontinuouson⁡Γcu(X,t)∣u∈H1},𝒰0={δu∣δu∈H1,δu=0on⁡Γu,δudiscontinuouson⁡Γcu(X,t)∣u∈H1}.The essential and natural boundary conditions are given by(14)u=u¯,X∈Γunt·P=t¯,X∈Γtnc·P=tc([[u]]),X∈Γc,
the index c referring to crack boundaries, the index t to natural boundaries, and the index u to essential boundaries.By substituting the approximation of the displacement fieldu and the virtual displacement field δu into the weak form (11), the discrete system of equations is obtained:
(15)∑J=1n‍∫ΩJ‍∇δuJ:PdΩ-∑J=1n‍∫Γt,J‍δu·t¯dΓ-∑J=1n‍∫ΩJ‍ϱδu·bdΩ+∫Γc,J‍δ[[u]]·tcdΓ+∑J=1n‍∫ΩJ‍ϱδu·u¨dΩ=0
or in vector-matrix form:
(16)MIJDJ=FIext-FIint⁡
with
(17)FIext=∫Γt‍(NI)TtdΓ+∫Ω‍(NI)TbdΩ+∫Γc‍[[(NI)T]]tcdΓ,FIint⁡=∫Ω‍(BI)TPdΩ,MIJ=∫Ω‍ϱNINJTdΩ.For computational efficiency, nodal integration is employed. Nodal integration does not drastically reduce the number of quadrature points; it was shown, for example, by [72] that the critical time step is increased by orders of magnitude; care has to be taken in the existence of cracks [73]. However, nodal integration leads to instabilities due to rank deficiency. Therefore, we employ the stabilized conforming nodal integration technique as mentioned previously. Stabilized conforming nodal integration has also shown great performance for fracture problems in the context of the finite element method; see [74–80]. The explicit central difference time integration scheme is used (for the time integration).


## 4. Fracture Criterion and Cohesive Zone Model

Fracture is governed by the maximum principal stress. The crack is introduced perpendicular to the direction of the maximum principal stress. Complex crack patterns are obtained automatically as the method does not require any representation of the crack's topology. There is no need to distinguish between crack nucleation and crack propagation and complex crack patterns including crack branching are natural outcome of the simulation.The cohesive zone model related the jump in the displacement field[[u]] to the cohesive traction at the crack surface. The discontinuous displacement field can be split into a part acting perpendicular and another tangential to the crack surface:
(18)δn=[[u]]·nδt=|[[u]]-δnn|.
Defining an effective crack opening displacement as suggested by [81]
(19)δ=β2δt2-δn2
an effective cohesive zone model can be derived by
(20)t=tmax⁡δmax⁡δifδ≤δmax⁡orδ˙<0.
The traction vector is then obtained by
(21)t=tδ(β2δt+δnn),
where β determines the amount of tangential and normal tractions.


## 5. Results

### 5.1. Crack Branching

Let us consider a specimen with initial crack as illustrated in Figure5. It is subjected to uniaxial tensile loading of σ=1MPa at the top and bottom. This classical benchmark problem of dynamic fracture has been studied by several people to test the accuracy and robustness of their computational method [39, 53, 82, 83]. Also, experimental data can be found for such type of problems; see, for example, [84–86]. The material properties for this example are modulus of elasticity E=32,000MPa and Poisson's ratio ν=0.20. The maximum crack speed is restricted by the Rayleigh wave speed cR=2119.0m/s. Discretizations ranging from only 4000 nodes up to more than 16000 nodes are investigated. Also the influence of the cracking rules will be demonstrated in this section.Figure 5
Plate with a horizontal initial notch under tensile tractions.The fracture patterns at various time steps are depicted in Figure6. Spurious cracks adjacent to the “main” crack are observed when no cracking rules are applied. The velocity of the propagating crack is illustrated in Figure 7. We note the following.(i)
The cracking rules do not influence the crack speed indicating that much less energy is dissipated in the spurious cracks.(ii)
The crack propagates faster just before it branches.(iii)
The maximum crack speed is far below the Raleigh wave speed. This agrees well with experimental observations of the microbranch instability problem as reported by [84–86]. The ability to capture this physical behaviour naturally is one advantage of the CPM over other methods with continuous crack surface [23, 39, 51, 87].Crack pattern for the crack branching problem.
(a)
4141 nodes; with spurious cracking
(b)
16281 nodes; with spurious cracking
(c)
4141 nodes; without spurious cracking
(d)
16281 nodes; without spurious crackingCrack tip speed for the crack branching problem.
(a)
Without spurious cracking
(b)
Comparison with and without spurious cracking


### 5.2.  The Kalthoff Experiment

The second classical benchmark example studied here is the Kalthoff problem [88]. Therefore consider the double-notched specimen under impact loading (v=20m/s) as illustrated in Figure 8. The impact leads to mode I dominated fracture with a crack propagating almost orthogonal to the impact loading.Figure 8
The Kalthoff problem.We exploit the symmetry of the model and carry out simulations of discretizations ranging from 10,000 to 40,000 particles. According to [89], the modulus of elasticity in these experiments is E=190GPa, the initial density is ϱ0=8000kg/m3, and Poisson’s ratio is ν=0.3. Figure 9 shows the final fracture pattern that matches the experimental data. Figure 10 reveals that the crack speed here is also not influenced by the cracking rules. However, the crack looks more erratic without using the cracking rules.Final crack path of the Kalthoff problem.
(a)
10201 nodes; with spurious cracking
(b)
10201 nodes; without spurious cracking
(c)
40804 nodes; with spurious cracking
(d)
40804 nodes; without spurious crackingFigure 10
Crack tip speed and cohesive fracture energy of the Kalthoff problem.


## 5.1. Crack Branching

Let us consider a specimen with initial crack as illustrated in Figure5. It is subjected to uniaxial tensile loading of σ=1MPa at the top and bottom. This classical benchmark problem of dynamic fracture has been studied by several people to test the accuracy and robustness of their computational method [39, 53, 82, 83]. Also, experimental data can be found for such type of problems; see, for example, [84–86]. The material properties for this example are modulus of elasticity E=32,000MPa and Poisson's ratio ν=0.20. The maximum crack speed is restricted by the Rayleigh wave speed cR=2119.0m/s. Discretizations ranging from only 4000 nodes up to more than 16000 nodes are investigated. Also the influence of the cracking rules will be demonstrated in this section.Figure 5
Plate with a horizontal initial notch under tensile tractions.The fracture patterns at various time steps are depicted in Figure6. Spurious cracks adjacent to the “main” crack are observed when no cracking rules are applied. The velocity of the propagating crack is illustrated in Figure 7. We note the following.(i)
The cracking rules do not influence the crack speed indicating that much less energy is dissipated in the spurious cracks.(ii)
The crack propagates faster just before it branches.(iii)
The maximum crack speed is far below the Raleigh wave speed. This agrees well with experimental observations of the microbranch instability problem as reported by [84–86]. The ability to capture this physical behaviour naturally is one advantage of the CPM over other methods with continuous crack surface [23, 39, 51, 87].Crack pattern for the crack branching problem.
(a)
4141 nodes; with spurious cracking
(b)
16281 nodes; with spurious cracking
(c)
4141 nodes; without spurious cracking
(d)
16281 nodes; without spurious crackingCrack tip speed for the crack branching problem.
(a)
Without spurious cracking
(b)
Comparison with and without spurious cracking


## 5.2.  The Kalthoff Experiment

The second classical benchmark example studied here is the Kalthoff problem [88]. Therefore consider the double-notched specimen under impact loading (v=20m/s) as illustrated in Figure 8. The impact leads to mode I dominated fracture with a crack propagating almost orthogonal to the impact loading.Figure 8
The Kalthoff problem.We exploit the symmetry of the model and carry out simulations of discretizations ranging from 10,000 to 40,000 particles. According to [89], the modulus of elasticity in these experiments is E=190GPa, the initial density is ϱ0=8000kg/m3, and Poisson’s ratio is ν=0.3. Figure 9 shows the final fracture pattern that matches the experimental data. Figure 10 reveals that the crack speed here is also not influenced by the cracking rules. However, the crack looks more erratic without using the cracking rules.Final crack path of the Kalthoff problem.
(a)
10201 nodes; with spurious cracking
(b)
10201 nodes; without spurious cracking
(c)
40804 nodes; with spurious cracking
(d)
40804 nodes; without spurious crackingFigure 10
Crack tip speed and cohesive fracture energy of the Kalthoff problem.


## 6. Conclusions

In conclusion, a cracking particles method based on nodal integration and updated Lagrangian kernels was proposed where a set of simple cracking rules was suggested in order to avoid spurious cracking. The stabilized nodal integration guarantees the computational efficiency of the method while maintaining the stability. Moreover, the updated Lagrangian kernel ensures that the method remains applicable also for extremely large deformations.The method was applied to two benchmark examples where the performance of the method was demonstrated.



---

*Source: 646514-2014-05-07.xml*