Dynamic fracture with continuum-kinematics-based peridynamics

1.   Introduction

Predicting crack propagation and material damage is still very challenging in computational mechanics. Fracture problems have been addressed by various computational methods such as damage models or discontinuous finite element discretizations ^[1,2,3] and phase-field fracture simulations ^[4,5,6]. All these approaches are based on the classical continuum mechanics assumption of a homogeneous bulk material. Peridynamics allows an alternative approach to fracture because it models the material in a non-local form. Initially introduced by Silling ^[7,8], peridynamics uses integral equations to describe the relative displacements and forces between material points. Concepts like stress and strain are absent, and the behavior of a material point is described solely by its interactions with other material points.

The original peridynamic concepts were restricted to the interaction of bonds, which has limited the ability to account for a material's Poisson ratio other than ${}^{1}\!\!\diagup\!\!{}_{4}\;$ for 3D problems. Other formulations, like ordinary state-based peridynamics and non-ordinary state-based peridynamics, address this problem, cf. ^[9]. A new approach was recently introduced by Javili, McBride & Steinmann, who propose a continuum-kinematics-based reformulation of peridynamics (CPD), ^[10,11,12]. This geometrically exact formulation relies on an analogy to the classical continuum mechanics and is intrinsically designed to capture the lateral contraction of the material correctly. Three types of material point interactions are introduced, namely bond-, surface- and volume-based interactions, which correspond to the invariants of a general deformation. The relationships between the material parameters of CPD and isotropic linear elasticity were recently elaborated for two- and three-dimensional problems, ^[13,14].

The new kinematics of CPD require a new concept of damage and fracture. Because of the three different types of interactions, it is no longer sufficient to understand material damage as a bond-based phenomenon. Thus we enrich the material description by kinematic variables that account for the loss of load-carrying capacity in the material's force density, which is also extended by a density related to contact. In that way, crack nucleation and propagation, as well as impact damage, can be modeled. To the knowledge of the authors, this is the first concept of damage within the CPD framework.

This manuscript is organized as follows. Section 2 introduces the notation and the necessary theory of CPD. Here the one-, two-, and three-neighbor interactions are defined, and the model is extended to contact of two or more bodies. Section 3 is the paper's core; here, we introduce our damage model. Sections 4–6 present numerical examples. We start with a model validation by a simple 2d crack growth for a mode I tension test in Section 4. Next, the crack initiation due to reflected impact waves is presented in Section 5 using the example of a curved bar ^[15]. Finally, in Section 6, CPD is used to compute damage induced by the impact of a sphere on a brittle plate.

2.   Theory of continuum-kinematics-based peridynamics

In peridynamics, a body is considered as a set of $N$ points in Euclidean space $\mathbb{R}^3$, and the dynamics is described by the movement of these points. Specifically, the bijective mapping Eq 1,

describes the transformation (deformation) of the body from reference configuration $\mathcal{B}_0$ to current configuration $\mathcal{B}_t$ at the time $t$. The movement of a body can therefore be described as a parametrical (temporal) sequence of deformations $\boldsymbol{\Phi}\colon \mathcal{B}_0 \times [0, T] \mapsto \mathbb{R}^3$.

The point position in material configuration is described by $\boldsymbol{X}^i$ and in current configuration as Eq 2,

with the displacement vector $\boldsymbol{u}(\boldsymbol{X}^i, t)$ and $i = 1, \dots, N$. Points interact only with other points inside of their specified neighborhood $\mathcal{H}_1^i$, which is defined as the set of points inside the spherical space with the radius $\delta \in \mathbb{R}^+$, also called the horizon $\delta$ (see Figure 1). Accordingly, $\mathcal{H}_1^i$ includes all points $\boldsymbol{X}^j$ inside the horizon of point $\boldsymbol{X}^i$ in the reference configuration of the body $\mathcal{B}_0$, Eq 3:

The equation of motion for point $i$ reads Eq 4:

with the density $\rho$, the point acceleration vector $\boldsymbol{\ddot{u}}$, and the point force density vectors $\boldsymbol{b}_0^{\mathrm{int}}$ and $\boldsymbol{b}_0^{\mathrm{ext}}$, which denote force per unit undeformed volume. The external force density $\boldsymbol{b}_0^{\mathrm{ext}}$ results from the external forces that are acting on the body and the internal force density $\boldsymbol{b}_0^{\mathrm{int}}$ from the interactions between the individual material points. Peridynamics can be understood as a Lagrangian particle method, because all equations are mapped to the reference configuration. In the following, the notation $\boldsymbol{u}^{i} = \boldsymbol{u}(\boldsymbol{X}^i, t)$ and $\boldsymbol{b}_0^{\mathrm{int}, \, i} = \boldsymbol{b}_0^{\mathrm{int}}(\boldsymbol{X}^i, t)$ will be used for improved readability.

All our simulations consider a short period of time, therefore an explicit time integration scheme is used. We employ the Velocity-Verlet algorithm of Littlewood ^[16], where for each time $t$ and material point $i$, the acceleration, velocity and displacement are calculated as Eq 5,

with the time step ${\vartriangle} t$.

Various formulations of peridynamics exist for the calculation of the internal force density, and all of them are based on the non-local interactions between the material points. The unique point of CPD is the use of three different types of interactions, also called one-, two- and three-neighbor interactions (see Figure 2). Correspondingly, for CPD, $\boldsymbol{b}_0^{\mathrm{int}, \, i}$ is the sum of the internal force densities of these interactions, thus Eq 6:

2.1.   One-neighbor interactions

The one-neighbor interaction of point $i$ and $j$, in standard peridynamics also called the bond, is defined in material and current configuration as Eq 7:

One-neighbor interactions can be interpreted as line elements with the initial length $L^{ij}$ in material notation and the deformed length $l^{ij}$ in current configuration. These so called relative length measures of the one-neighbor interaction are defined as Eq 8:

It is assumed, that all one-neighbor interactions of point $i$ contribute equally. Therefore, an effective one-neighbor volume is defined as Eq 9,

with $N_1^i$ being the number of one-neighbor interactions for point $i$ and the neighborhood volume by Eq 10:

with the factor $\beta^i \in [0, 1]$ that takes the fullness of the neighborhood into account. As an example, it applies $\beta^i = 1$ if the neighborhood of point $i$ is completely inside the body $\mathcal{B}_0$. In contrast, if the neighborhood of point $i$ is partially outside the body $\mathcal{B}_0$, the factor $\beta^i < 1$ works as a correction factor to the volume $V_{\mathcal{H}}^i$.

The force density due to one-neighbor interactions is defined as Eq 11,

with the harmonic potential Eq 12,

and the one-neighbor interaction constant $C_1$. The constant $C_1$ can be interpreted as a resistance against the length change of one-neighbor interactions. With Eq 13:

the internal force density $\boldsymbol{b}_{1}^{\mathrm{int}, \, i}$ of point $i$ can be formulated as Eq 14,

2.2.   Two-neighbor interactions

Two-neighbor interactions are area elements, respectively triangles, spanned by the points $\boldsymbol{X}^i$, $\boldsymbol{X}^j$ and $\boldsymbol{X}^k$. They are constructed by two corresponding one-neighbor interactions $\boldsymbol{\Delta X}^{ij}$ and $\boldsymbol{\Delta X}^{ik}$ of point $i$. One important condition is that the distance between the points $\boldsymbol{X}^j$ and $\boldsymbol{X}^k$ needs to be smaller than the horizon $\delta$, as displayed in Figure 3. Therefore, the set of all corresponding point-sets for two-neighbor interactions of point $i$ is defined as Eq 15,

The deformation of two-neighbor interactions is mainly described by the relative area measure, in material and current notation defined as Eqs 16, 17:

and as scalar quantities the areas

The force density due to two-neighbor interactions is defined as Eqs 18–20:

with the harmonic potential

and the effective two-neighbor volume

The number of two-neighbor interactions of point $i$ is $N_2^i$. The two-neighbor interaction constant $C_2$ can be interpreted as a resistance against the area change. With Eqs 21, 22:

the internal force density $\boldsymbol{b}_{2}^{\mathrm{int}, \, i}$ of point $i$ can be formulated as

2.3.   Three-neighbor interactions

Three-neighbor interactions are volume elements, precisely tetrahedrons, spanned by the points $\boldsymbol{X}^i$, $\boldsymbol{X}^j$, $\boldsymbol{X}^k$ and $\boldsymbol{X}^l$. They are constructed by the three corresponding one-neighbor interactions $\boldsymbol{\Delta X}^{ij}$, $\boldsymbol{\Delta X}^{ik}$ and $\boldsymbol{\Delta X}^{il}$ of point $i$. For a valid three-neighbor interaction, the conditions Eq 23,

must be met. Consequently, the set of all corresponding point-sets for three-neighbor interactions of point $i$ is defined as Eq 24,

The deformation of three-neighbor interactions is mainly described by the relative volume measure, in material and current notation defined as Eq 25,

The force density due to three-neighbor interactions is defined as Eq 26,

with the harmonic potential Eq 27,

and the effective three-neighbor volume Eq 28,

The number of three-neighbor interactions of point $i$ is $N_3^i$. The three-neighbor interaction constant $C_3$ can be interpreted as a resistance against the volume change. With Eq 29,

the internal force density $\boldsymbol{b}_{3}^{\mathrm{int}, \, i}$ of point $i$ can be formulated as Eq 30,

2.4.   Contact

For the modeling of contact between peridynamic bodies, we employ the algorithm by Silling and Askari ^[8], which is in detail described in ^[16]. The approach of this algorithm is mainly based on short-range forces that are activated at a certain threshold of the point distance.

For the incorporation of contact, Eq 4 needs to be extended to Eq 31,

with the contact force density $\boldsymbol{b}_0^{\mathrm{con}, \, i}$. We consider two different peridynamic bodies $\mathcal{B}^I$ and $\mathcal{B}^{II}$, discretized with the point spacings $\boldsymbol{\Delta x_I}$ and $\boldsymbol{\Delta x_{II}}$. Then we define the contact point sets $\mathcal{H}_{t}^{\mathrm{con}, \, i}$ and $\mathcal{H}_{t}^{\mathrm{con}, \, j}$ as Eqs 32, 33

with the critical contact distance $l_c \approx \mathrm{max}\left(\frac{\boldsymbol{\Delta x_I}}{2}, \, \frac{\boldsymbol{\Delta x_{II}}}{2}\right)$. Now, the contact force densities yield to Eqs 34, 35:

with the contact spring constant $C^{\mathrm{con}}$ and the point volumes $V^i$ and $V^j$.

3.   Modeling damage with continuum-kinematics-based peridynamics

In classical peridynamics, damage is modeled by the failure of one-neighbor interactions. The failure quantity for the strain-based damage model reads Eq 36,

with the one-neighbor interaction stretch Eq 37,

and the critical one-neighbor interaction stretch $\varepsilon_c$. This stretch can assumed to be identical to the critical bond stretch in classical peridynamics, which has been estimated by Madenci and Oterkus ^[9] as Eq 38,

with the critical energy release rate $G_c$, the shear modulus $G$ and the bulk modulus $K$. The pointwise damage quantity $D^i$ incorporates the whole neighborhood, and is defined as Eq 39,

These equations cannot directly be used to model damage within the continuum-kinematics-based framework, because they do not take two- or three-neighbor interactions into consideration. Applying this damage model alone will not lead to crack paths but to diffuse failure zones, because two- or three-neighbor interactions are still active and lead to forces between failed points.

To address this problem, failure quantities for two- and three-neighbor interactions, $d_2^{ijk}$ and $d_3^{ijkl}$, are introduced. Here we propose that two- and three-neighbor interactions fail, if one or more corresponding one-neighbor interactions fail. Therefore, the failure quantity for two-neighbor interactions can be defined as Eq 40,

and for three-neighbor interactions as Eq 41,

With these failure quantities, we re-define the internal force density for one-neighbor interactions Eq 14 as:

for two-neighbor interactions Eq 22 as:

and for three-neighbor interactions Eq 30 as:

In such a manner, the failed point interactions do not contribute to the internal material response and their damaging effect is considered.

4.   Numerical results for the mode I tension test

In the following section, it is shown, that CPD is able to model crack growth for two- and three-dimensional problems. Therefore, a square with edge length $l$ and a predefined crack of length $a = \frac{1}{2}l$ is subjected to tension due to the expansion of the upper and lower region of the model with a constant velocity $v_0 = {0.005}\text{ m s}^{-1}$ (see Figure 4).

For the 2D setup, a uniformly distributed point cloud with $200 \times 200$ points, and for the 3D setup $60 \times 60 \times 3$ points are used. Both setups use the material parameters of steel, as listed in Table 1. As derived by Ekiz, Javili, and Steinmann ^[13], the interaction constants of the two-dimensional setup are Eq 45,

with the Young's modulus $E$ and the Poisson's ratio $\nu$. Furthermore, the constants for one- and three-neighbor interactions of the three-dimensional setup are defined as Eq 46,

for $C_2 = 0$, and with the first and second Lamé parameter $\lambda = \frac{E \, \nu}{(1 + \nu)(1-2\nu)}$ and $\mu = \frac{E}{2(1+\nu)}$ ^[14].

In Figure 5, the damage $D^i$ for the 2D and the 3D setup is shown. The crack propagates and grows as expected until the square is broken into two pieces for both setups. Without the additional failure quantities $d_2^{ijk}$ and $d_3^{ijkl}$, a diffuse damage field and no clear crack path would be the result of this simulations. Consequently, CPD can be used to model crack-growth with our proposed damage model.

5.   Curved bar under pressure

In the following section, crack initiation due to impact is investigated for a two- and three-dimensional discretizations. For this purpose, a model of a curved bar is subjected to pressure waves, which are supposed to superimpose inside the material and eventually lead to crack initiation. The material points are spatially distributed along the curve $f(x) = \cos(\frac{\pi}{2} \, x)$ (see ^[15] for more details). As shown in Figure 6, for each root point $x_i$ on $f$, $N_n$ points occur along the curve Eq 47,

with the derivative $f'(x) = \frac{\mathrm{d}f}{\mathrm{dx}} = - \frac{\pi}{2} \, \sin(\frac{\pi}{2} \, x)$. The bar has the width $W_B = {62.5}\text{ mm}$ along the curve $n_i(x)$ and is defined for the root points $x_i \in [-\frac{L_B}{2}, \frac{L_B}{2}]$ with the bar length $L_B = {1}\text{ m}$. The number of points $N_n$ on $n_i(x)$ is a measure to describe the density of the point cloud, since it is used to specify the point spacing $\Delta x = \frac{W_B}{N_n}$. For the three-dimensional implementation, $N_n$ layers of material points are equally distributed with distance $\Delta x$ along $z \in [-\frac{W_B}{2}, \frac{W_B}{2}]$. A coarse point cloud with $N_n = 5$ is shown in Figure 7 for the purpose of illustrating the discretization.

On each side of the curved bar, a pressure impulse Eq 48,

with the pressure peak $p_0$ and the impulse duration $t_1$ is applied for one layer of material points in the left and right boundary (see Figure 8). The pressure is applied via the external body force density Eqs 49–51,

with the normal vector

for the left side and

for the right side of the bar, and the angle $\alpha = \arctan{\left(-\frac{1}{f'(L_B/2)}\right)}$.

The material parameters used for the calculations are shown in Table 2. For the two-dimensional setup, a pressure impulse with the peak $p_0 = {4 × 10^5}\text{ N/m}$ and for the three-dimensional setup, $p_0 = {1× 10^6}\text{ N/m}^2$ is used. For both setup's, the pulse has the duration $t_1 = {300}$ µs. Remark that for 2D, the body force density $\boldsymbol{b}^{\mathrm{ext}}_0$ has the unit [N/m²]. The interaction constants are calculated as before (see Eqs 45, 46).

In Figure 9, the damage $D^i$ of the two-dimensional setup is shown. After $t = {2.9}\text{ ms}$, a crack in the middle of the bar is visible. The pressure waves propagate through the bar and then get reflected, which consequently transforms them into tensile waves. These tensile waves then lead to the initiation of a crack. The waves continue to propagate in the model and when superimposed again, the same effect occurs and more cracks are formed. The two-dimensional model reproduces this behavior very well, since exactly these further cracks occur for time $t = {5.8}\text{ ms}$.

The same behavior can also be observed with the 3D model (see Figure 10). As an important remark, for the visualization of the waves in the model, a stress tensor was artificially calculated. The calculation of the stresses is not part of the peridynamics and is only used to illustrate the wave reflection. After the first reflection of the pressure wave, a single crack is initiated in the center of the model. Also the cracking due to the further superposition of the waves can be detected, as seen in Figure 11. The position differs from that of the 2D model, but this could be explained by the versatile influencing factors of CPD, such as material parameters and different discretizations. Here further studies are necessary. In summary it can be stated, that CPD can be used to map cracking due to the material response to pressure waves.

6.   Impact damage

In the following section, damage due to contact between two peridynamic bodies is investigated for our proposed damage model. Here, the shot of a sphere with initial velocity $v_0 = {50}\text{ m s}^{−1}$ through a circular disc that is free in space is computed numerically. As it is displayed in Figure 12, the sphere has the radius $r = {10}\text{ mm}$ and the disc the radius $R = {250}\text{ mm}$ and height $h = {10}\text{ mm}$. The material parameters for both bodies are listed in Table 3. For the modeling of the contact, the search radius $l_c = {2.5}\text{ mm}$ and the contact spring constant $C^{\mathrm{con}} = {1000}\text{ GPa}$ are used.

In Figure 13, the results of the numerical experiment are displayed for $t = {0.26}\text{ ms}$ (Figure 13a) and $t = {0.86}\text{ ms}$ (Figure 13b). The sphere punches through the circular disc and stamps a hole in it. Slight damage occurs in the impact area, but no further cracks propagate as a result of the impact. Similar material response can be seen for variations of material or contact parameters.

7.   Summary

This study presents a new approach to dynamic fracture and impact damage within the CPD. We extend the classical uni-axial damage model and introduce failure quantities for two- and three-neighbor interactions. With these quantities, two- and three-neighbor interactions are deactivated if the corresponding one-neighbor interactions fail, which allows cracks to be modeled within the continuum-based framework.

Three different effects were investigated numerically: crack growth, crack initiation, and damage due to contact. At first, we found that for two- and three-dimensional simulations, our approach handles crack growth of a mode I tension test very well. At second, to investigate crack initiation, a curved bar is subjected to pressure waves. The pressure waves superimpose and lead to tensile waves that initiate cracks in the model. We observed both a first crack in the middle of the model as well as secondary cracks due to the further superposition of the waves. At third, additional impact simulations show that CPD is even able to model a punchthrough impact without further cracking. Clearly, the proposed damage model can capture effects due to contact between two peridynamic bodies.

To sum it all up, we have introduced a new damage model that extends the novel CPD framework and now opens up versatile possibilities for the simulation of fracture and damage. We could show that our damage model can represent the growth of existing cracks, the initiation of cracks, and the damage due to the contact of two bodies.

Acknowledgments

The authors gratefully acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG) in the projects WE 2525/15-1 and WI 1430/9-1.

Conflict of interest

The authors declare that they have no conflict of interest.