A peridynamics-based finite element method for quasi-static fracture analysis

The bond-based peridynamic model is proposed by Silling in [8], which assumes that a point ${\boldsymbol{x}}$ interacts with all the points in its neighborhood, $H_{\delta({\boldsymbol{x}})}$, where $\delta$ is a horizon that denotes the cut-off radius of the interaction. Based on this, the equation of motion at point ${\boldsymbol{x}}$ yields

where ${\boldsymbol{b}}({\boldsymbol{x}})$ is the external body force field, ${\boldsymbol{\xi}}={\boldsymbol{x^{\prime}}}-{\boldsymbol{x}}$ is a relative position vector referred to as a bond, and ${\boldsymbol{f}}({\boldsymbol{\xi}})$ is the pairwise force field of the bond ${\boldsymbol{\xi}}$. Particularly, ${\boldsymbol{f}}({\boldsymbol{\xi}})$ denotes the nonlocal force vector that point ${\boldsymbol{x}}^{\prime}$ exerts on point ${\boldsymbol{x}}$. For quasi-static problems, the equilibrium equation at point ${\boldsymbol{x}}$ can be obtained by setting $\ddot{{\boldsymbol{u}}}\left({\boldsymbol{x}}\right)=0$ in Eq. (1).

To ensure the balance of the linear momentum, the pairwise force vector ${\boldsymbol{f}}({\boldsymbol{\xi}})$ must be anti-symmetric [8], i.e., ${\boldsymbol{f}}({\boldsymbol{\xi}})=-{\boldsymbol{f}}(-{\boldsymbol{\xi}}^{\prime})$. Thus, the pairwise force is assumed to be central [14], i.e.

where $\hat{{\boldsymbol{f}}}[{\boldsymbol{x}}]\left\langle{\boldsymbol{\xi}}\right\rangle$ (respectively $\hat{{\boldsymbol{f}}}[{\boldsymbol{x}}^{\prime}]\left\langle-{\boldsymbol{\xi}}\right\rangle$) is the partial interaction due to the action of point ${\boldsymbol{x^{\prime}}}$ over point ${\boldsymbol{x}}$ (with respect to point ${\boldsymbol{x}}$ over point ${\boldsymbol{x^{\prime}}}$). For homogeneous materials with linear elasticity and small deformations, a possible constitutive model [8, 15] is

where ${\boldsymbol{u}}$ is the displacement field. ${\boldsymbol{C}}({\boldsymbol{\xi}})$ is the micromodulus tensor of the bond ${\boldsymbol{\xi}}$, which is defined as [8]

where $c(\|{\boldsymbol{\xi}}\|)$ is the micromodulus coefficient.

To sum up, for quasi-static problems, the governing equations of the bond-based peridynamics for a configuration $\Omega$, $\Omega\subset\mathbb{R}^{d}(d=1,2,3)$, can be summarized as

where ${{\boldsymbol{u}}}^{*}$ is the prescribed displacement on $\partial\Omega_{{\boldsymbol{u}}}$, ${\boldsymbol{\eta}}$ is a measure of deformation of bond ${\boldsymbol{\xi}}.$ Eq. (5a), Eq. (5b) and Eq. (5c) are the static admissibility, constitutive equation and kinematic admissibility and compatibility, respectively.

In the peridynamics, once the failure of the structure is considered, the bond break is needed. Different kinds of bond break criteria have been introduced in the literature [10, 16, 17], and here we adopt the stretch-based criterion proposed by Silling and Askari [10]. The stretch of a bond ${\boldsymbol{\xi}}$ is defined as

Then the bond break can be recorded with a history-dependent scalar-valued function, $\mu$, which is defined as

where $t$ and $\tau$ denote computational steps. $s_{crit}$ is the critical bond stretch. Note that Eq. (7) means the brittle fracture mode. Multiply $\mu({\boldsymbol{\xi}},t)$ to the right-hand side of Eq. (5b), and then the constitutive equation including bond break is obtained.

Based on $s_{crit}$, the critical energy dissipation of the broken bond yields [18]

Note that the critical energy dissipation $\omega_{crit}$ is different for bonds with different lengths. With $\mu$ and $\omega_{crit}$ at hand, the effective damage for each point ${\boldsymbol{x}}$ is defined as [18]

which then indicates the failure of the structure.

In 2016, Han et al.[19] derived the principle of minimum potential energy for a hybrid classical continuum mechanics and state-based peridynamic model. As a special case of the hybrid model, the principle of minimum potential energy for the bond-based peridynamic model can be directly obtained as: the solution of Eq. (5) is also the solution of

where

is the kinematically admissible displacement space [13], and the total potential energy $\Pi({\boldsymbol{u}})$ is defined as

where the first and second terms on the right-hand side of Eq. (12) are the deformation energy and external work, respectively.

Besides, it has been proved in [19] that the sufficient and necessary condition for $\Pi({\boldsymbol{u}})$ to reach the minimum yields $\delta\Pi({\boldsymbol{u}})={\boldsymbol{0}}$, i.e.,

where

is the trial space.

The classical FEM is based on the classical continuum mechanics, in which the expression of the total potential energy is a single integral on the configuration $\Omega$. However, in the peridynamics, the expression of the total potential energy contains a double integral term, i.e., the deformation energy of the structure, see Eq. (12), which results in the inconsistency of the finite element framework for the peridynamics with that of the classical FEM.

To avoid the above issue, we reconstruct the formulation of the deformation energy. Firstly, note that ${\boldsymbol{f}}({\boldsymbol{\xi}})=0,\forall{\boldsymbol{x^{\prime}}}\notin H_{\delta({\boldsymbol{x}})}$, thus the inner integral defined on $H_{\delta({\boldsymbol{x}})}$ can be extended to the entire configuration $\Omega$. Therefore, the deformation energy of $\Omega$ can be rewritten as

Further, we define a new type of integral operation as

where $\bar{\Omega}$ is an integral domain generated by two $\Omega$s, $\bar{{\boldsymbol{g}}}({\boldsymbol{x^{\prime}}},{\boldsymbol{x}})$ is a double-parameter function related to ${\boldsymbol{g}}({\boldsymbol{\xi}})$ and defined on $\bar{\Omega}$. Then the total potential energy $\Pi({\boldsymbol{u}})$ can be presented in a single integral form, i.e.

In this section, we introduce the peridynamic elements (PEs) in the domain $\bar{\Omega}$ to preserve the single integral characteristic when discretizing the total potential energy.

Before introducing PEs, we briefly review the elements in classical FEM because the PEs are generated based on them. In classical FEM, the configuration $\Omega$ is divided by a set of elements, $\left\{e_{i}\right\}_{i=1}^{m}$, where $m$ is the total number of the elements. These elements are not overlapped but share edges and vertices, called nodes, between adjacent elements. Therefore, we name the element in classical FEM as continuous element (CE) to distinguish it from peridynamic element (PE). The set of mesh nodes is denoted as $\left\{{\mathrm{P}}_{l}\right\}_{l=1}^{n}$, with $n$ the total number of the nodes.

Fig. 1 shows the relation between PEs and CEs in the 2D case. Suppose $\Omega$ is divided into a set of CEs, $\left\{e_{i}\right\}_{i=1}^{m}$. For each $e_{i}$, its neighborhood $H_{i}$ is defined as the minimal coverage of the union of the neighborhood of all points in $e_{i}$, as shown in Fig. 1a. Then for each CE $e_{j}\in H_{i}$, including $e_{i}$, there will be a new PE, denoted as $\bar{e}_{k}$, combined by $e_{j}$ and $e_{i}$, as shown in Fig. 1b. Fig. 1c displays some PEs in detail, one can find that PEs may have different geometries, but they are all composed of two CEs.

In general, if we assume that CE $e_{i}$ has $n_{i}$ nodes and the nodes labels are $\left[{\begin{array}[]{*{20}{c}}{{{\mathrm{P}}_{{i_{1}}}}}&{{{\mathrm{P}}_{{i_{2}}}}}&\cdots&{{{\mathrm{P}}_{{i_{{n_{i}}}}}}}\end{array}}\right]$, CE $e_{j}$ has $n_{j}$ nodes and the nodes labels are $\left[{\begin{array}[]{*{20}{c}}{{{\mathrm{P}}_{{j_{1}}}}}&{{{\mathrm{P}}_{{j_{2}}}}}&\cdots&{{{\mathrm{P}}_{{j_{{n_{j}}}}}}}\end{array}}\right]$, where ${\mathrm{P}}_{{\alpha}_{\beta}}\in\left\{{\mathrm{P}}_{l}\right\}_{l=1}^{n}(\alpha=i,j;\beta=1,2,\cdots,n_{\alpha})$, then PE $\bar{e}_{k}$ is an element with $\bar{n}_{k}=n_{i}+n_{j}$ nodes and the nodes labels are $\left[{\begin{array}[]{*{20}{c}}{{{\mathrm{P}}_{{k_{1}}}}}&{{{\mathrm{P}}_{{k_{2}}}}}&\cdots&{{{\mathrm{P}}_{{k_{{\bar{n}_{k}}}}}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{{\mathrm{P}}_{{j_{1}}}}}&{{{\mathrm{P}}_{{j_{2}}}}}&\cdots&{{{\mathrm{P}}_{{j_{{n_{j}}}}}}}&{{{\mathrm{P}}_{{i_{1}}}}}&{{{\mathrm{P}}_{{i_{2}}}}}&\cdots&{{{\mathrm{P}}_{{i_{{n_{i}}}}}}}\end{array}}\right]$. Based on the CEs, $\left\{e_{i}\right\}_{i=1}^{m}$, and the above method of generating PE, a set of PEs, $\left\{\bar{e}_{k}\right\}_{k=1}^{\bar{m}}$, will finally be obtained. For example, $n_{i}=n_{j}=4$ and $\bar{n}_{k}=8$ in Fig. 1c.

Now we have two sets of elements, CEs and PEs. The former is used to calculate local quantities, and the latter is used to calculate nonlocal quantities.

According to the interpolation technique of the classical FEM, the displacement field ${\boldsymbol{u}}_{i}({\boldsymbol{x}})$ on any CE $e_{i}$ can be approximately expressed as

where

are the shape function matrix of CE and the nodal displacement vector of CE for $e_{i}$, respectively. Here, $\left\{N_{i_{l}}({\boldsymbol{x}})\right\}_{l=1}^{n_{i}}$ is the shape function for node ${\mathrm{P}}_{i_{l}}$, $\left\{u_{i_{l}}\right\}_{l=1}^{n_{i}}$, $\left\{v_{i_{l}}\right\}_{l=1}^{n_{i}}$ and $\left\{w_{i_{l}}\right\}_{l=1}^{n_{i}}$ are the displacement components of node ${\mathrm{P}}_{i_{l}}$ along the $X$, $Y$, and $Z$ directions, respectively.

The displacement field $\bar{{\boldsymbol{u}}}_{k}({\boldsymbol{x^{\prime}}},{\boldsymbol{x}})$ on any PE $\bar{e}_{k}$ can be approximately expressed as

where

are the shape function matrix of PE and the nodal displacement vector of PE for $\bar{e}_{k}$, respectively.

Then, based on Eq. (5c) and Eq. (25), for any PE $\bar{e}_{k}$ we have

where

is the difference matrix for shape function for $\bar{e}_{k}$. In Eq. (33),

is the difference operator matrix and ${\boldsymbol{I}}$ is an identity matrix of dimension $d\ (d=1,2,3)$.

Next, base on Eq. (5b) and Eq. (32), for any PE $\bar{e}_{k}$ we have

where

is the partial interaction force matrix for $\bar{e}_{k}$. In Eq. (36), ${\boldsymbol{D}}({\boldsymbol{\xi}})$ is the matrix form of the micromodulus tensor ${\boldsymbol{C}}({\boldsymbol{\xi}})$. For $d=3$ as an example, we have

With Eq. (18), Eq. (25), Eq. (32) and Eq. (35) at hand, the total potential energy $\Pi({\boldsymbol{u}})$ can be approximated as

where ${\boldsymbol{d}}$ is the total nodal displacement vector and

are the total stiffness matrix and total load vector, respectively. In Eq. (39), $\bar{{\boldsymbol{G}}}_{k}$ and ${\boldsymbol{G}}_{i}$ are the transform matrix of the degree of freedom for the nodes of $\bar{e}_{k}$ and the transform matrix of the degree of freedom for the nodes of $e_{i}$, respectively, which satisfies

respectively. Further,

are the element stiffness matrix and the element load vector, respectively.

Finally, a linear system including the solution of the nodal displacement vector ${\boldsymbol{d}}$ can be derived from Eq. (38) using the condition Eq. (13). That is

We first conduct a four-point bending test of a single-edge notched beam. The geometry and the loading setup are shown in Fig. 3a and the mesh is shown in Fig. 3b. The Young’s modulus is $E=30\text{GPa}$ and the critical stretch is set to be $s_{crit}=0.02$. The average size of the CEs is $h\approx 1.0\text{mm}$. The simulation is implemented through 20 equably progressive increments steps, i.e., $\text{N}=20$.

Fig. 4 shows the evolution of the effective damage contours during the loading process of the four-point bending for the notched beam. First, the damage initiates, i.e., the bond breaks for the first time, at step 13. Then, the damage develops slowly in the next several steps. After that, the damage propagates suddenly and violently at step 18, which is in line with brittle fracture characteristics.

Fig. 5 displays the curves of the average force on the loading points on the top of the beam versus the displacement. Before step 17, the effective damage is small, so the curve approximates linear. Then between step 17 and step 18, the curve drops drastically, which is related to the sudden propagation of the crack between these two steps.

We now investigate the compression test of a central notched Brazilian disk. The geometry and the loading setup of the disk are shown in Fig. 6a and the mesh is shown in Fig. 6b. The Young’s modulus is $E=3.1\text{GPa}$ and the critical stretch is set to be $s_{crit}=0.02$. The average size of the CEs is $h\approx 0.7\text{mm}$. The simulation is implemented through 64 equably progressive increments steps, i.e., $\text{N}=64$.

Fig. 7 displays the effective damage contours at given loading steps. The simulation results indicate that the cracks initiate at the notched corners at step 38. Then the cracks grow slowly at both corners until step 63 and suddenly penetrate the disk at step 64. The crack geometries of the disk are very similar to the experimental results reported in [21, Fig. 17].

Fig. 8 shows the curves of the average force on the loading points on the top of the disk versus the displacement. Before the effective damage occurs, the curve rises with a constant slope, and then the slope gradually decreases with the increase of breaking bonds. Then, along with the sudden propagation of the cracks between step 63 and step 64, the curve drops drastically.

In this example, we consider the mixed mode fracture of a double-edge-notched plate. The geometry and the loading setup of the plate are shown in Fig. 9a and the mesh is shown in Fig. 9b. The Young’s modulus is $E=30\text{GPa}$ and the critical stretch is set to be $s_{crit}=0.02$. The average size of the CEs is $h\approx 1.25\text{mm}$. The simulation is performed through 25 equably progressive increments steps, i.e., $\text{N}=25$.

The effective damage evolution contours of this test are shown in Fig. 10. The damage appears first at the notched corners at step 9. At step 11, the damage is more obvious than at step 9, and it could be found that the damage of the left notched corner propagates down, and the damage of the right notched corner propagates up. Such asymmetric crack path stems from the asymmetry of the boundary conditions, and similar results were reported in [22, Fig. 17] and [23, Fig. 21] with the hybrid model. Then at step 12, the cracks propagate destructively.

Fig. 11 shows the curves of the average force on the loading points on the left and upper edge of the plate versus the displacement. It can be seen that before step 12, the tensile loads on the upper and lower edges of the board play a major role. However, with the destructive propagation of the cracks at step 12, the plate can almost no longer withstand tensile loads. Therefore, the reaction force on the upper edge keeps a very small level after step 12. On the other hand, the force curve related to the left edge keeps rising after a slight drop, which means that the plate can still withstand the shear loads after step 12, and then it is mainly damaged under the shear loads.