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.
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 part , sometimes also referred to as the “usual” part, and the discontinuous or “enriched” part : The discretization of the continuous part of the displacement field is based on moving least squares (MLS) shape functions [71] of linear completeness: where are the meshless shape functions of node at position and are nodal parameters at time . Note that (1) the nodal parameters are not the true physical displacement values at node and (2) the shape functions are expressed in terms of material coordinates . It can be shown that the shape functions are given by with where denotes the kernel function and its support size; 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 th 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 Figure 2): being the set of cracking particles are additional degrees of freedom (DOF) and the enrichment function is the step function 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
(b) Discretization with crack segments
It is advantageous to write the approximation of the displacement field in vector matrix notation: and store the “standard” and enriched shape functions in the matrix ; the DOFs of the “standard” and enriched part of the displacement field are stored in the vector . Moreover, we can define the matrix and the matrix containing the derivatives of the meshless shape functions 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 , 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.
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 Figure 4, the deviation in this angle is calculated by being the set of new cracking particles. Crack branching is assumed when .
3. Weak Form and Discretization
The linear momentum equation is where is the nominal stress tensor, is the density, are body forces, and the superimposed dots denote material time derivatives. The displacement and traction boundary conditions are where the index refers to the crack, the index refers to traction boundaries, and the index refers to displacement boundaries.
We solve the equation of motion in weak form that can be stated in variational form. Find the displacement field such that the first variation in the energy is zero: with The approximation spaces and are given by
The essential and natural boundary conditions are given by the index referring to crack boundaries, the index to natural boundaries, and the index to essential boundaries.
By substituting the approximation of the displacement field and the virtual displacement field into the weak form (11), the discrete system of equations is obtained: or in vector-matrix form: with
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 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: Defining an effective crack opening displacement as suggested by [81] an effective cohesive zone model can be derived by The traction vector is then obtained by 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 Figure 5. It is subjected to uniaxial tensile loading of 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 and Poisson's ratio . The maximum crack speed is restricted by the Rayleigh wave speed . 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.
The fracture patterns at various time steps are depicted in Figure 6. 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].
(a) 4141 nodes; with spurious cracking
(b) 16281 nodes; with spurious cracking
(c) 4141 nodes; without spurious cracking
(d) 16281 nodes; without spurious cracking
(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 () as illustrated in Figure 8. The impact leads to mode I dominated fracture with a crack propagating almost orthogonal to the impact loading.
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 , the initial density is , and Poisson’s ratio is . 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.
(a) 10201 nodes; with spurious cracking
(b) 10201 nodes; without spurious cracking
(c) 40804 nodes; with spurious cracking
(d) 40804 nodes; without spurious cracking
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.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.