Phase-field material point method for dynamic brittle fracture with isotropic and anisotropic surface energy

TL;DR

This paper introduces a new phase-field material point method for simulating dynamic brittle fracture in elastic media, explicitly accounting for anisotropic surface energy and impact-driven fracture, verified through analytical and experimental comparisons.

Contribution

A novel phase-field material point method that explicitly incorporates anisotropic surface energy for robust dynamic fracture simulation in elastic media.

Findings

Method accurately predicts crack paths influenced by anisotropy.

Verifies robustness through analytical and experimental benchmarks.

Outperforms standard phase field finite element methods in key scenarios.

Abstract

A novel phase field material point method is introduced for robust simulation of dynamic fracture in elastic media considering the most general case of anisotropic surface energy. Anisotropy is explicitly introduced through a properly defined crack density functional. The particular case of impact driven fracture is treated by employing a discrete field approach within the material point method setting. In this, the equations of motion and phase field governing equations are solved independently for each discrete field using a predictor-corrector algorithm. Contact at the interface is resolved through frictional contact conditions. The proposed method is verified using analytical predictions. The influence of surface energy anisotropy and loading conditions on the resulting crack paths is assessed through a set of benchmark problems. Comparisons are made with the standard Phase Field Finite Element Method and experimental observations.

Tables

	HTA fibre	6376 epoxy	laminate ($\varphi$ = 0o)	laminate ($\varphi$ = 90o)
Young’s modulus [MPa]	235000	3600	136000	8750
Tensile strength [MPa]	3920	105	1670	60
Density [kg/m3]	1770	1310	1586	1586

Equations

$${\Psi }_s={\Psi }_{el}+{\Psi }_f=\underset{\Omega }{\int }{\psi }_{el}\left(\mathbf{\epsilon }\right)dV+\underset{\Gamma }{\int }{\mathcal{G}}_c\left(\theta \right)d\Gamma$$

$$\bm{\varepsilon}=\frac{1}{2}\Big{(}\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\Big{)}$$

$${\Psi }_f=\underset{\Gamma }{\int }{\mathcal{G}}_c\left(\theta \right)d\Gamma \approx \underset{\Omega }{\int }{\bar{\mathcal{G}}}_c{\mathcal{Z}}_{c,Anis}dV$$

$${\mathcal{Z}}_{c,Anis}=\left[\frac{{\left(c-1\right)}^2}{4l_0}+l_0|\nabla c{|}^2+l_0^3\sum _{\begin{array}{c}ijkl\end{array}}{\gamma }_{ijkl}\frac{{\partial }^{2}c}{\partial x_i\partial x_j}\frac{{\partial }^{2}c}{\partial x_k\partial x_l}\right]$$

$$\mathbf{\gamma }=\left[\begin{array}{cccccc}{\gamma }_{1111}& {\gamma }_{1122}& {\gamma }_{1133}& {\gamma }_{1112}& {\gamma }_{1123}& {\gamma }_{1113}\\ {\gamma }_{2211}& {\gamma }_{2222}& {\gamma }_{2233}& {\gamma }_{2212}& {\gamma }_{2223}& {\gamma }_{2213}\\ {\gamma }_{3311}& {\gamma }_{3322}& {\gamma }_{3333}& {\gamma }_{3312}& {\gamma }_{3323}& {\gamma }_{3313}\\ {\gamma }_{1211}& {\gamma }_{1222}& {\gamma }_{1233}& {\gamma }_{1212}& {\gamma }_{1223}& {\gamma }_{1213}\\ {\gamma }_{2311}& {\gamma }_{2322}& {\gamma }_{2333}& {\gamma }_{2312}& {\gamma }_{2323}& {\gamma }_{2313}\\ {\gamma }_{1311}& {\gamma }_{1322}& {\gamma }_{1333}& {\gamma }_{1312}& {\gamma }_{1323}& {\gamma }_{1313}\end{array}\right].$$

$$\mathbf{\gamma }=\left[\begin{array}{ccc}{\gamma }_{1111}& {\gamma }_{1122}& {\gamma }_{1112}\\ {\gamma }_{2211}& {\gamma }_{2222}& {\gamma }_{2212}\\ {\gamma }_{1211}& {\gamma }_{1222}& {\gamma }_{1212}\end{array}\right].$$

$${\mathbf{x}}_{\theta }={\mathbf{R}}_{\theta }\mathbf{x}$$

$${\mathcal{G}}_c\left(\theta \right)=\underset{\Gamma }{\int }{\mathcal{G}}_c\left(\theta \right)d\Gamma \approx {\int }_{-\infty }^{+\infty }{\bar{\mathcal{G}}}_c{\mathcal{Z}}_{c,Anis}d{x}_{2_{\theta }}\approx {\int }_{-x_{lb}}^{+x_{lb}}{\bar{\mathcal{G}}}_c{\mathcal{Z}}_{c,Anis}d{x}_{2_{\theta }}$$

$${\bar{\mathbf{f}}}_1^{cont}={\bar{\mathbf{f}}}_1^{nor}+{\bar{\mathbf{f}}}_1^{tan}={\bar{f}}_1^{nor}\cdot {\mathbf{n}}_1^{cont}+{\bar{f}}_1^{tan}\cdot {\mathbf{s}}_1^{cont}$$

$$\dot{\mathcal{K}}\left(\dot{\mathbf{u}}\right)+{\dot{\mathcal{W}}}^{int}\left(\dot{\mathbf{u}},\dot{c},\nabla \dot{c}\right)-{\dot{\mathcal{W}}}^{ext}\left(\dot{\mathbf{u}}\right)-{\dot{\mathcal{W}}}^{cont}\left(\dot{\mathbf{u}}\right)=0$$

$$\nabla \dot{c}=\frac{d}{dt}\left(\frac{\partial c}{\partial x_i}\right)$$

$$\dot{\mathcal{K}}\left(\dot{\mathbf{u}}\right)=\frac{d}{dt}\underset{\Omega }{\int }\frac{1}{2}\rho |\dot{\mathbf{u}}{|}^{2}dV$$

$${\dot{\mathcal{W}}}^{int}\left(\dot{\mathbf{u}},\dot{c},\nabla \dot{c}\right)=\frac{d{\Psi }_s}{dt}=\frac{d}{dt}\underset{\Omega }{\int }\left({\psi }_{el}+{\bar{\mathcal{G}}}_c{\mathcal{Z}}_{c,Anis}\right)dV$$

$${\psi }_{el}=g(c){\psi }_{el}^{+}+{\psi }_{el}^{-}$$

$$g=(1-k_f)c^2+k_f$$

$$\mathbf{\sigma }=\frac{\partial {\psi }_{el}}{\partial \mathbf{\epsilon }}.$$

$${\dot{\mathcal{W}}}^{ext}\left(\dot{\mathbf{u}}\right)={\int }_{\partial {\Omega }_{\bar{t}}}(\bar{\mathbf{t}}\cdot \dot{\mathbf{u}})d\partial {\Omega }_{\bar{t}}+{\int }_{\Omega }(\mathbf{b}\cdot \dot{\mathbf{u}})dV.$$

$${\dot{\mathcal{W}}}^{cont}\left(\dot{\mathbf{u}}\right)={\int }_{\partial {\Omega }_{\bar{f}}}({\bar{\mathbf{f}}}^{cont}\cdot \dot{\mathbf{u}})d\partial {\Omega }_{\bar{f}}={\int }_{\partial {\Omega }_{\bar{f}}}\left(({\bar{\mathbf{f}}}^{nor}+{\bar{\mathbf{f}}}^{tan})\cdot \dot{\mathbf{u}}\right)d\partial {\Omega }_{\bar{f}}.$$

$${\dot{\mathcal{W}}}^{cont}\left(\dot{\mathbf{u}}\right)={\dot{\mathcal{W}}}_1^{cont}\left({\dot{\mathbf{u}}}_1\right)+{\dot{\mathcal{W}}}_2^{cont}\left({\dot{\mathbf{u}}}_2\right)=$$

$${\int }_{\partial {\Omega }_{1\bar{f}}}({\bar{\mathbf{f}}}_1^{cont}\cdot {\dot{\mathbf{u}}}_1)d\partial {\Omega }_{1\bar{f}}+{\int }_{\partial {\Omega }_{2\bar{f}}}({\bar{\mathbf{f}}}_2^{cont}\cdot {\dot{\mathbf{u}}}_2)d\partial {\Omega }_{2\bar{f}}=0$$

$$\{\begin{array}{ll}\nabla \cdot \mathbf{\sigma }+\mathbf{b}=\rho \ddot{\mathbf{u}}& \text{on }\Omega \\ & \\ \begin{array}{rl}\left(\frac{4l_0\left(1-k_f\right)\mathcal{H}}{{\bar{\mathcal{G}}}_c}+1\right)c& -4l_0^2\Delta c\\ & +4l_0^4\sum _{\begin{array}{c}ijkl\end{array}}{\gamma }_{ijkl}\frac{{\partial }^{4}c}{\partial x_i\partial x_j\partial x_k\partial x_l}=1\end{array}& \text{on }\Omega \end{array}$$

$$\begin{array}[]{*{20}{c}}{\psi_{e}^{+}-\mathscr{H}\leq 0}&{\dot{\mathscr{H}}\geq 0}&{\mathscr{\dot{H}}\left({\psi_{e}^{+}-\mathscr{H}}\right)=0}.\end{array}$$

$$\{\begin{array}{ll}\mathbf{\sigma }\cdot \mathbf{n}=\bar{\mathbf{t}},& \text{on }\partial {\Omega }_{\bar{t}}\\ \mathbf{u}=\bar{\mathbf{u}},& \text{on }\partial {\Omega }_{\bar{u}}\\ \mathbf{u}=\text{prescript}(0)\mathbf{u},& \text{on }\text{prescript}(0)\Omega \\ \dot{\mathbf{u}}=\text{prescript}(0)\dot{\mathbf{u}},& \text{on }\text{prescript}(0)\Omega \\ \ddot{\mathbf{u}}=\text{prescript}(0)\ddot{\mathbf{u}},& \text{on }\text{prescript}(0)\Omega \\ \left[4l_0^2\nabla c-2l_0^4{\sum }_{\begin{array}{c}ijkl\end{array}}{\gamma }_{ijkl}\left(\frac{{\partial }^{3}c}{\partial x_j\partial x_k\partial x_l}\right)-2l_0^4{\sum }_{\begin{array}{c}ijkl\end{array}}{\gamma }_{ijkl}\left(\frac{{\partial }^{3}c}{\partial x_i\partial x_j\partial x_k}\right)\right]\cdot \mathbf{n}=0,& \text{on }\partial \Omega \\ 2l_0^4{\sum }_{\begin{array}{c}ijkl\end{array}}{\gamma }_{ijkl}\left(\frac{{\partial }^{2}c}{\partial x_k\partial x_l}\right)+2l_0^4{\sum }_{\begin{array}{c}ijkl\end{array}}{\gamma }_{ijkl}\left(\frac{{\partial }^{2}c}{\partial x_i\partial x_j}\right)=0,& \text{on }\partial \Omega \\ c=\text{prescript}(0)c,& \text{on }\text{prescript}(0)\Omega \end{array}$$

$${\mathbf{n}}_1^{cont}=-{\mathbf{n}}_2^{cont},$$

$${\bar{\mathbf{f}}}_1^{nor}=-{\bar{\mathbf{f}}}_2^{nor},$$

$${\bar{f}}^{nor}\le 0,$$

$${\gamma }_n\le 0,$$

$${\gamma }_n{\bar{f}}^{nor}=0,$$

$${\mathbf{s}}_1^{cont}=-{\mathbf{s}}_2^{cont},$$

$${\bar{\mathbf{f}}}_1^{tan}=-{\bar{\mathbf{f}}}_2^{tan},$$

$$|{\bar{f}}^{tan}|\le {\mu }_f|{\bar{f}}^{nor}|,$$

$$|{\gamma }_s|\ge 0,$$

$$|{\gamma }_s|\left(|{\bar{f}}^{tan}|-{\mu }_f|{\bar{f}}^{nor}|\right)=0,$$

$${\bar{f}}^{nor}={\bar{\mathbf{f}}}_1^{cont}\cdot {\mathbf{n}}_1^{cont}={\bar{\mathbf{f}}}_2^{cont}\cdot {\mathbf{n}}_2^{cont}$$

$${\bar{f}}^{tan}={\bar{\mathbf{f}}}_1^{cont}\cdot {\mathbf{s}}_1^{cont}={\bar{\mathbf{f}}}_2^{cont}\cdot {\mathbf{s}}_2^{cont}.$$

Full text

††journal: Journal of Computer Methods in Applied Mechanics and Engineering

1 Introduction

Failure of materials subjected to dynamic loading is commonly associated with complex yet intriguing phenomena, i.e., crack merging, branching and arrest [1, 2]. These phenomena become even more pronounced in the case of anisotropy and high-rate loading conditions, e.g., impact. Anisotropy governs the fracture response of both natural and manufactured materials as in the case of granitic rocks [3], biological tissues [4], single crystals [5], and composite sheets [6]. Furthermore, the response of such materials under impact loading is being receiving considerable attention as it pertains to numerous industrial applications particularly within the automotive and aerospace sector, see, e.g., [7]. Numerical simulation of fracture propagation under such conditions can provide valuable insight into the underlying mechanical processes while also providing a framework for optimum design of materials considering their post-fracture response under impact loading. However, robust and accurate simulation of impact driven dynamic fracture is a challenging task as it requires the fusion of robust fracture propagation modelling with contact induced non-linearities and large displacement kinematics.

Within the framework of Computational Fracture Mechanics, a variety of methods has been introduced to address the problem of crack propagation. Among the most commonly used mesh based methods are the element deletion method [8], the Cohesive Zone Method [9, 10] , the eXtended Finite Element Method (XFEM) [11, 12] and crack-driving configurational force approaches [13, 14]. In these methods, algorithmic tracking of individual cracks is required as these evolve, merge, or branch. This results in considerable increase of the underlying computational complexity especially in the three dimensional case. Furthermore, an ad-hoc crack growth criterion is required for crack evolution.

Francfort and Marigo [15] introduced a framework for avoiding these issues by establishing brittle fracture as an energy minimization problem within a robust variational structure. More recently, Bourdin et al. [16] provided a regularization of the variational formulation which is more suitable for numerical solution schemes using as point of departure the phase field approximation of the Mumford-Shah potential presented in [17]. Within this variational setting, brittle fracture is formulated as a coupled, i.e., displacement and phase field problem, and the crack path naturally emerges from the solution of corresponding field equations.

As a result, standard re-meshing or enrichment strategies near the crack tip and the requirement for algorithmic tracking of the crack front are avoided. Complex crack topologies, e.g., crack merging and/or branching as well as applications to three dimensional domains are efficiently resolved in the same manner, see, e.g., [18]. Finite element based phase field formulations have been introduced to treat brittle [19, 20], ductile fracture [21, 22] and hydraulic fracture [23, 24]. Phase field models for anisotropic fracture have been presented [25, 26, 27, 28] although not within a dynamic setting. Very recently, Hesch et al. [29] have developed a method to resolve contact problems involving isotropic phase field fracture. In this formulation, a finite element based mortar contact algorithm in conjunction with a hierarchical refinement scheme is employed that reduces computational costs although relying on the predefinition of contact areas. Therefore, an adaptive hierarchical refinement is required for arbitrary impact fracture problems to resolve the local contact features.

The Material Point Method (MPM) [30] has been introduced as a promising alternative to computationally expensive particle based methods that can efficiently deal with contact and large displacement problems. MPM is an extension of Particle-In-Cell (PIC) methods that efficiently treats history-dependent variables. In MPM, the continuum is represented by a set of Lagrangian particles, i.e., the material points, that are mapped onto a non-deforming Eulerian mesh (computational grid) where the governing equations are solved. This combined Eulerian-Lagrangian approach has been proven particularly advantageous in problems pertinent to high material and geometric nonlinearities since the distortion error is minimized [31, 32, 33]. Within this context, MPM has already been used to simulate very challenging engineering problems e.g. penetration [34], cutting process simulations [35] and solid-fluid interaction problems [36, 37].

To this point, few research has been conducted in damage simulation utilizing MPM using either discrete [38, 39], cohesive [40, 41], or continuum damage models [42, 43]. Taking advantage of the good qualities of phase field modelling in naturally resolving complex crack paths, a Phase Field Material Point Method (PF-MPM) has been successfully introduced by the authors in [44] for quasi-static brittle fracture problems while a variant accounting for anisotropy in the quasi-static regime has been developed in [45].

Moving beyond the state-of-the-art, we present a phase field MPM method for the solution of dynamic fracture considering materials with anisotropic fracture energy; isotropy emerges as a special case of the proposed formulation. Following, the method is extended to also account for frictional contact fracture problems. We use as point of departure the MPM contact algorithm introduced in Bardenhagen et al. [46] where multiple fields, termed discrete fields, are introduced in the non-deforming Eulerian mesh so that each contact body corresponds to a different field. We define the variational structure of our phase field implementation of impact driven fracture at each discrete field from which the coupled weak form of the contact problem naturally emerges. Finally, we develop a predictor corrector solution algorithm for the solution of the governing equations over time.

This paper is organized as follows. In section 2 phase-field modelling is briefly described in both isotropic and anisotropic brittle fracture. The discrete field formulation for phase field fracture due to impact is presented in section 3. The Material Point Method implementation for frictional contact fracture is presented in section 4. Finally, in section 5, a set of benchmark problems are examined to demonstrate the accuracy and robustness of the proposed method.

2 Preliminaries

2.1 Phase field modelling

In the following, the case of an arbitrary deformable domain $\Omega$ is considered, with an external boundary $\partial\Omega$ and a crack path $\Gamma$ as shown in Fig. (1a). The deformable domain $\Omega$ with domain volume $V$, is subjected to body forces $\bm{b}=\left[\begin{matrix}b_{1}&b_{2}&b_{3}\end{matrix}\right]^{T}$. Furthermore, a set of traction/pressure loads $\bm{\bar{t}}$ is applied on the boundary $\partial\Omega_{\bar{t}}\subseteq\partial\Omega$. A prescribed displacement field, denoted as $\mathbf{\bar{u}}$, is imposed on the boundary $\partial\Omega_{\bar{u}}\subseteq\partial\Omega$.

According to Griffith’s theory [47] the stored energy ${\Psi_{s}}$ of the body $\Omega$ can be expressed as

[TABLE]

where ${\Psi_{el}}$ and ${\Psi_{f}}$ are the elastic strain energy and the fracture energy (surface energy), respectively. Moreover, $\psi_{el}\left(\bm{\varepsilon}\right)$ corresponds to the elastic energy density and $\bm{\varepsilon}$ is the symmetric strain tensor which under the small strain assumption is defined as

[TABLE]

The $\left(\nabla\right)$ symbol in Eq. (2) stands for the gradient operator and $\mathbf{u}(\mathbf{x},t)$ for the displacement field of a point $\mathbf{x}=\left[\begin{matrix}x_{1}&x_{2}&x_{3}\end{matrix}\right]^{T}$ at time $t$.

Due to material anisotropy, the critical fracture energy density $\mathscr{G}_{c}\left(\theta\right)$ in equation (1) explicitly depends on the orientation angle of the crack $\theta\left(s\right)$, $s\in\Gamma$. In the 2D case, the orientation angle is defined as the angle between the tangent vector at any point to the crack path $\Gamma$ and the horizontal. In the 3D case, the orientation can be defined by considering the direction cosines of the normal to the tangent plane of the fracture surface with respect to the global coordinate system.

In principle, the stored energy $\Psi_{s}$ is known provided that both $\Gamma$ and $\theta$ at the current configuration are known. Hence, the computational treatment of elastic fracture mechanics gives rise to a nonlinear problem whereby standard procedures revert to path tracking and optimization algorithms to predict and resolve the crack path as this evolves. In the phase field approximation, the path dependent fracture energy surface integral is transformed into a volume integral defined over the entire domain $\Omega$ (Bourdin et al. [16]) - see also, Fig. 1b. Hence, the phase field approximation gives rise to equation (3)

[TABLE]

where the functional $\mathscr{Z}_{c,Anis}=\mathscr{Z}_{c,Anis}\left(c,\theta\right)$ and $c$ is the phase field. Parameter $\bar{\mathscr{G}}_{c}$ in Eq. (3) corresponds to the energy required to create a unit area of fracture surface $\mathscr{A}_{c}\left(\theta\right)$. This is assumed to be constant for all directions even though the actual surface energy of an anisotropic material is direction dependent. To resolve this inconsistency, directionality of the fracture toughness is accounted for in the definition of the functional $\mathscr{Z}_{c,Anis}$ by introducing a fourth-order anisotropic sensor $\gamma$ as discussed in section 2.2. This formulation enables for a direct comparison to be drawn between isotropic and anisotropic models by controlling the value of a single parameter as further discussed in section 5. Although strain rate dependence of the fracture toughness has been reported in the literature [48], such concepts are beyond the scope of this work.

2.2 Anisotropic crack density functional

To account for the general case of anisotropic material behaviour, $\mathscr{Z}_{c,Anis}$ is defined as the fourth-order functional utilized in [25, 45] according to Eq. (4)

[TABLE]

where $c(\mathbf{x},t)\in\left[0,1\right]$ is the phase field defined over the domain $\Omega$, $l_{0}\in\mathbb{R}^{+}$ is a length scale parameter and $\gamma_{ijkl}$, $i,j,k,l=1\dots 3$ are the components of the fourth-order anisotropic tensor corresponding to the anisotropic constitutive behaviour of the material. Phase field values of $c=1$ correspond to uncracked regions of the domain $\Omega$. Conversely, values of $c=0$ correspond to cracked regions. The length scale parameter $l_{0}$ controls the width of the regularized crack topology.

Remark 1

A second-order functional can be employed to model anisotropy on the fracture properties. However, on the modelling side, fourth-order phase field functionals have been shown to successfully and robustly account for strong anisotropies [49, 26] while avoiding ill-posedness associated with second-order anisotropic models [50]. On the simulation side, they improve the convergence rate of the underlying Newton procedure [51].

Remark 2

A mathematical proof on the $\Gamma-$ convergence of the fourth-order anisotropic theory presented in this work has not been yet established. The extensive numerical studies performed in this work, within the bounds of the anisotropic tensors consider, hint that solutions provided by the fourth-order functional used in this work indeed converge. A relevant discussion on the aspect of $\Gamma-$ convergence for the fourth-order isotropic functional can be found in [51].

The anisotropic tensor $\bm{\gamma}$ is conveniently defined in the three dimensional space utilizing Voigt notation as

[TABLE]

To demonstrate the versatility of the anisotropic functional in describing different material symmetries we focus on the 2D case for brevity in which case, $\bm{\gamma}$ reduces to

[TABLE]

The direction angle of the crack path $\theta$ can be explicitly introduced in the expression of $\mathscr{Z}_{c,Anis}$ through a coordinate transformation,i.e., by transforming the Cartesian coordinate system $\mathbf{x}=\left[\begin{matrix}x_{1}&x_{2}\end{matrix}\right]^{T}$ to $\mathbf{x}_{\theta}=\left[\begin{matrix}x_{\theta_{1}}&x_{\theta_{2}}\end{matrix}\right]^{T}$ where axis ${x}_{1_{\theta}}$ is defined along the crack path $\Gamma$ and axis ${x}_{2_{\theta}}$ is normal to the crack interface (see Fig. 1). Thus, the transformation relation (7) holds

[TABLE]

where $\theta$ is the counter-clockwise angle between the ${x}_{1}$-axis and ${x}_{1_{\theta}}$ and $\mathbf{R}_{\theta}$ is the standard 2D rotation matrix. Detailed derivations on the transformation are provided by the authors in Appendix B of [45]. Eventually, the surface energy density $\mathscr{G}_{c}\left(\theta\right)$ for each angle $\theta$ is cast in the following form

[TABLE]

where $x_{lb}$ is the distance from the crack. Eq. (8) enables numerical evaluation and visualization of the anisotropic surface energy density in polar form. Integration is performed along the normal to the crack path where the phase field variations are significant; assuming a value $x_{lb}=50l_{0}$ yields a reasonable approximation.

Figs. 2a and 2b illustrate the surface energy densities $\mathscr{G}_{c}\left(\theta\right)$ and their reciprocals $1/\mathscr{G}_{c}\left(\theta\right)$, respectively in polar coordinates for isotropic symmetry with the second and fourth order phase field models. The cases of cubic and orthotropic symmetry are also shown in Fig. 2.

To derive these polar plots, the parameter ${{{\bar{\mathscr{G}}_{c}}}}$ is chosen to be ${{{\bar{\mathscr{G}}_{c}}}}=0.70710$ kN/m for the fourth order isotropic, cubic and orthotropic symmetry whereas ${{{\bar{\mathscr{G}}_{c}}}}=1$ kN/m for second order isotropic symmetry. The parameter ${{{\bar{\mathscr{G}}_{c}}}}$ is chosen so that all previously mentioned models have the same minimum value of surface energy density $\mathscr{G}_{c_{min}}=1$ kN/m.

Remark 3

Similar polar plots can be derived for the 3D case by considering the transformation of the global coordinate system to the coordinate system defined by the tangent plane at the fracture surface and its normal. It is useful to note that the 3D equivalent of Eq. 8 is still a line integral as integration is performed along the normal to the fracture surface. Such aspects are beyond the scope of this work; an intuitive approach on the rotation of anisotropic tensors is provided in [52].

3 Governing equations for phase field fracture due to impact

3.1 Derivation of the coupled strong form for impact-fracture problems

In this section, the governing equations for contact induced brittle fracture are introduced. For brevity, the case of two bodies is presented herein. In Fig. 3a, two deformable domains, namely $\Omega_{1}$ and $\Omega_{2}$ are considered such that $\Omega_{1}\cup\Omega_{2}=\Omega$. Their external boundaries are defined as $\partial\Omega_{1}$ and $\partial\Omega_{2}$, respectively.

At time $t$ the two bodies are in contact along the surface $\partial\Omega_{1\bar{f}}=\partial\Omega_{2\bar{f}}=\partial\Omega_{\bar{f}}\subseteq\partial\Omega$. Furthermore, two crack paths are defined, i.e., $\Gamma_{1}$ and $\Gamma_{2}$ at $\Omega_{1}$ and $\Omega_{2}$, respectively, under the action of a set of tractions $\bm{\bar{t}}$ and body forces $\bm{b}=\left[\begin{matrix}b_{1}&b_{2}&b_{3}\end{matrix}\right]^{T}$.

When the two bodies are in contact, a contact force $\bm{\bar{f}}_{1}^{cont}$ is applied to body $\Omega_{1}$ from body $\Omega_{2}$. This is defined in component form according to Eq. (9)

[TABLE]

where $\bm{\bar{f}}_{1}^{nor}$ and $\bm{\bar{f}}_{1}^{tan}$ stand for the normal and tangential contact force vectors whereas $\bar{f}_{1}^{nor}$ and $\bar{f}_{1}^{tan}$ are their corresponding components. The normal and tangential surface unit vectors on contact surface $\partial\Omega_{1\bar{f}}$ are denoted as $\mathbf{n}_{1}^{cont}$ and $\mathbf{s}_{1}^{cont}$, respectively. Similarly, a contact force $\bm{\bar{f}}_{2}^{cont}$ is applied from $\Omega_{1}$ to $\Omega_{2}$ with components $\bm{\bar{f}}_{2}^{nor}$, $\bm{\bar{f}}_{2}^{tan}$ and normal and tangential surface unit vectors $\mathbf{n}_{2}^{cont}$ and $\mathbf{s}_{2}^{cont}$ being defined accordingly.

Employing a phase field representation of fracture allows for a robust derivation of the impact-fracture strong form by considering the energy balance equation (10)

[TABLE]

where $\dot{\mathscr{K}}\left(\dot{\mathbf{u}}\right)$ is the rate of the kinetic energy, $\dot{\mathscr{W}}^{int}\left(\dot{\mathbf{u}},\dot{c},\nabla\dot{c}\right)$ is the rate of internal work, $\dot{\mathscr{W}}^{ext}\left(\dot{\mathbf{u}}\right)$ is the rate of the work done by external forces, and $\dot{\mathscr{W}}^{cont}\left(\dot{\mathbf{u}}\right)$ is the rate of work done by contact forces. Furthermore, $\dot{\mathbf{u}}={d\mathbf{u}}/{dt}$ corresponds to the velocity field, $\dot{c}={dc}/{dt}$ is the phase field time derivative, and $\nabla\dot{c}$ corresponds to the rate of the phase field spatial derivative, i.e.,

[TABLE]

for $i=1,...,3$.

The kinetic energy rate functional $\dot{\mathscr{K}}\left(\dot{\mathbf{u}}\right)$ is expressed as

[TABLE]

where $\rho$ corresponds to the mass density.

The rate of internal work is expressed as

[TABLE]

where the phase field approximation introduced in Eq. (3) is employed to define the fracture energy corresponding to the crack paths $\Gamma_{1}$ and $\Gamma_{2}$ - see also Fig. 3b).

In this work, $\psi_{el}$ in Eq. (13) is decomposed into a purely tensile $\psi_{el}^{+}$ and a purely compressive $\psi_{el}^{-}$ parts according to the strain energy density decomposition introduced by Miehe et al. [53] that is based on a spectral decomposition of the strain tensor. In this, the strain energy density is defined as

[TABLE]

where $g(c)\in\left[0,1\right]$ is a degradation function expressed as

[TABLE]

and $0\leq k_{f}\ll 1$ is a model parameter to treat potential ill-conditioning. In this work the model parameter is considered to be $k_{f}=0$ with no impact on the results as also highlighted by Braides [54]. The stress field is derived from Eq. (16) as

[TABLE]

Eq. (14) is adopted herein for the purpose of verification however different schemes, also with significant computational advantages, can be found in the literature (see, e.g., [55] and [49] for a review of existing models).

The rate of the external work functional $\dot{\mathscr{W}}^{ext}\left(\dot{\mathbf{u}}\right)$ is defined as

[TABLE]

Finally, the rate of work done by contact forces $\dot{\mathscr{W}}^{cont}\left(\dot{\mathbf{u}}\right)$ is expressed as

[TABLE]

Clearly, $\dot{\mathscr{W}}^{cont}\left(\dot{\mathbf{u}}\right)$ must vanish as the contributing forces are always opposite. However, $\dot{\mathscr{W}}^{cont}\left(\dot{\mathbf{u}}\right)$ is retained in the energy balance equation and is further decomposed into discrete field components; this greatly facilitates numerical approximation as will be highlighted in section 4. Therefore, Eq. (18) is expressed as

[TABLE]

where $\dot{\mathbf{u}}_{1}$ and $\dot{\mathbf{u}}_{2}$ are the velocity fields at body $\Omega_{1}$ and $\Omega_{2}$, respectively.

Applying the divergence theorem in Eq. (10), performing the necessary algebraic manipulations, and finally considering that the resulting expression must hold for arbitrary values of $\dot{\mathbf{u}}$ and $\dot{c}$, the strong form of the problem is derived as (see [45] for details)

[TABLE]

where $\mathscr{H}$ is a history field defined as the maximum value of the tensile part of the elastic energy density $\psi_{el}^{+}$ obtained in time domain $[0,t]$ and $\ddot{\mathbf{u}}={d\dot{\mathbf{u}}}/{dt}$ corresponds to the acceleration field.

In our implementation, we use the history field (see, e.g., Miehe et al. [56]) to enforce the irreversibility condition pertinent to the crack propagation problem, i.e., $\prescript{(t)}{}{\Gamma}\subseteq\prescript{(t+\Delta t)}{}{\Gamma}$ according to the following Kuhn-Tucker conditions for loading and unloading, i.e.,

[TABLE]

The coupled field equations (20) are subject to the set of boundary and initial conditions defined in Eq. (21)

[TABLE]

where $\mathbf{n}$ stands for the outward unit normal vector on the boundary.

Furthermore, the coupled field equations (20) are subjected to the kinematic constraints presented in Eqs. (22) to (26) and (27) to (31) at contact surface $\partial\Omega_{\bar{f}}$ [57]. The kinematic constraints of Eqs. (22) to (26) correspond to the normal contact laws

[TABLE]

whereas Eq. (27) to (31) correspond to the tangential contact and friction laws, where the Coulomb friction model is adopted.

[TABLE]

Kinematic constraints (22), (27) and (23) and (28) are imposed to satisfy Newton’s third law at the contact surface $\partial\Omega_{\bar{f}}$. Condition (24) is imposed on the normal component of the contact force that is defined according to Eq. (32)

[TABLE]

and implies a non-tension, i.e., non-stick, condition at the contact surface $\partial\Omega_{{}\bar{f}}$.

The tangential component $\bar{f}^{tan}$ is defined accordingly as

[TABLE]

The impenetrability condition (25) is imposed to ensure no penetration between the contact surfaces $\partial\Omega_{{1}\bar{f}}$ and $\partial\Omega_{{2}\bar{f}}$ when the two bodies are in contact.

3.1.1 Discrete field formulation for the coupled governing equations

In this work, a discrete field approach is adopted for the robust and efficient numerical treatment of contact dynamics between deformable bodies whereby each body is treated independently as discrete field. In the general case, it is assumed that the entire domain consists of a set of independent discrete fields $\{\mathscr{D}\mid\mathscr{D}=1,2,\dots,N_{\mathscr{D}}\}$, where $N_{\mathscr{D}}\in\mathbb{Z}^{+}$ stands for the total number of discrete fields and $\mathscr{D}$ indexes the $\mathscr{D}^{th}$ discrete field. Furthermore, all corresponding quantities that belong to discrete field $\mathscr{D}$, i.e., body $\Omega_{\mathscr{D}}$, are denoted with the subscript $\mathscr{D}$. Hence, in the two body case considered in this section $\{\mathscr{D}\mid\mathscr{D}=1,2\}$.

Within the discrete field setting, the contact forces arising from the interaction of the discrete fields are treated as additional external forces. Hence, the energy balance equation (10), is re-defined for each discrete field $\mathscr{D}$ as

[TABLE]

where $\dot{\mathscr{W}}^{cont}_{\mathscr{D}}\left(\dot{\mathbf{u}}_{\mathscr{D}}\right)$ is the rate of work done by contact forces and is expressed as

[TABLE]

Thus, the coupled strong form introduced in Eqs. (20) is now defined for the discrete field $\mathscr{D}$ as

[TABLE]

The set of boundary and initial conditions introduced in Eqs. (21) are modified for each discrete field $\mathscr{D}$ accordingly as

[TABLE]

where the last boundary condition is due to the contact forces that in this implementation are considered as forces applied externally to the discrete field $\mathscr{D}$.

4 Material Point Method for dynamic anisotropic fracture

Dynamic fracture under impact naturally involves large displacement kinematics especially in the pre- and post-fracture regime, e.g., in the case of high velocity projectile impact problems. To accurately resolve the pre and post fracture kinematics, the Material Point Method [30] is used in this work to solve the system of coupled governing Eqs. (36).

In the Material Point Method framework employed herein, the entire domain $\Omega=\Omega_{1}\cup\Omega_{2}$ is discretized into a set of material points $\mathscr{P}=\{p\mid p=1,2,\dots,N_{p}\}$, where $N_{p}\in\mathbb{Z}^{+}$ is the total number of material points whereas $p$ indexes the $p^{th}$ material point. It is assumed herein that $N_{\mathscr{D}p}$ material points belong to discrete field $\mathscr{D}$, i.e. body $\Omega_{\mathscr{D}}$ (see Fig. 3c).

According to the MPM approximation, the mass density $\rho_{\mathscr{D}}$ and domain volume $V_{\mathscr{D}}$ corresponding to the discrete field $\mathscr{D}$ are additively decomposed into the corresponding material point contributions according to Eqs. (38) and (39), respectively, i.e.,

[TABLE]

and

[TABLE]

where $\mathbf{x}_{\mathscr{D}}$ is the position vector of discrete field $\mathscr{D}$ and $\delta$ is the Dirac delta function. The material point mass density is defined as $\rho_{\mathscr{D}p}=M_{\mathscr{D}p}/V_{\mathscr{D}p}$ where $M_{\mathscr{D}p}$ and $V_{\mathscr{D}p}$ are the material point mass and volume, respectively. Furthermore, $\mathbf{x}_{\mathscr{D}p}$ corresponds to the position vector of material point $p$ at discrete field $\mathscr{D}$.

These material points are moving within a fixed computational grid, i.e. an Eulerian grid. The Eulerian grid is a non-deforming mesh that consists a set of $N_{n}\in\mathbb{Z}^{+}$ grid nodes and $N_{cells}\in\mathbb{Z}^{+}$ grid cells (see Fig. 3c). The material points are mapped onto the Eulerian grid where the governing equations are solved. The updated solution is mapped back from grid nodes to the material points. Finally, the background grid is reset and the computational cycle proceeds. The steps of the MPM are shown in Fig. 4. In this work, mapping from material points to grid nodes and vice versa is implemented by utilizing higher order B-splines interpolation functions.

4.1 Discrete equilibrium equations for contact dynamics

Defining appropriate trial solution and weighting function spaces for the displacement field, i.e.,

[TABLE]

and

[TABLE]

respectively, the discrete form of the equations of motion introduced in the first of Eqs. (36) is expressed for each discrete field $\mathscr{D}$ as

[TABLE]

where $\mathbf{w}$ are weighting functions that satisfy the homogeneous essential boundary conditions of the problem [58].

Substituting the material point approximation introduced in Eqs. (38) and (39) into Eq. (40), Eq. (41) is established

[TABLE]

Next, the weighting functions $\mathbf{w}_{\mathscr{D}p}$ and their spatial derivatives $\mathbf{\nabla w}_{\mathscr{D}p}$ are interpolated in the Galerkin sense according to relations (42)

[TABLE]

and (43), respectively,

[TABLE]

where $N_{I}(\mathbf{x}_{p})$ are the higher-order B-spline interpolation functions evaluated at the material point positions $\mathbf{x}_{\mathscr{D}p}$ [45, 59]. Furthermore, $\mathbf{w}_{\mathscr{D}I}$ are the weighting function, values evaluated at the background grid nodes; $I$ refers to the $I^{th}$ grid node.

Similar expressions are established for the displacement, velocity, and acceleration field, i.e.,

[TABLE]

[TABLE]

and

[TABLE]

respectively, where $\mathbf{u}_{\mathscr{D}I}$, $\dot{\mathbf{u}}_{\mathscr{D}I}$, and $\ddot{\mathbf{u}}_{\mathscr{D}I}$ are the components of the nodal displacement, velocity and acceleration vectors, respectively, evaluated at node $I$.

Substituting Eqs. (42) and (43) in relation (41) and performing the necessary algebraic manipulations, the following expression is established

[TABLE]

where $\bm{F}^{irt}_{\mathscr{D}I}$ are the nodal components of the inertia forces evaluated as

[TABLE]

whereas $\bm{F}^{int}_{\mathscr{D}I}$ are the nodal components of the internal forces

[TABLE]

Similarly, the nodal components of the external force vector $\bm{F}^{ext}_{\mathscr{D}I}$ assume the following form

[TABLE]

Finally, $\bm{F}^{cont}_{\mathscr{D}I}$ corresponds to the contact force nodal vector defined as

[TABLE]

As the weighting functions in Eq. (41) are chosen arbitrarily, Eq. (47) should hold for every set of nodal values $\mathbf{w}_{\mathscr{D}I}$. Hence, the following equilibrium equation is finally established

[TABLE]

where $\bm{R}^{u}_{\mathscr{D}I}$ is the nodal residual force vector at grid node $I$.

Finally, substituting Eq. (46) in Eq. (48), Eq. (52) is rewritten in the following form

[TABLE]

where $\bm{M}^{u}_{\mathscr{D}}$ is the global lumped mass matrix of the structure whose $M_{\mathscr{D}I}^{u}$ component is expressed as

[TABLE]

Eq. (53) lends itself conveniently into an explicit predictor-corrector time integration scheme as will be further discussed in Section 4.3.

4.2 Discrete phase field equations

The discrete form of the anisotropic phase field governing equations introduced in the second of Eqs. (36) can be also derived on the basis of the Material Point setting and Galerkin approximation. Similarly to the case of the displacement field, the phase field $c$ and the corresponding weighting functions $q$ are defined with respect to the following spaces, i.e.,

[TABLE]

and

[TABLE]

respectively.

Hence, the weak form of the phase field governing equations is expressed for each discrete field $\mathscr{D}$ as

[TABLE]

Introducing the MPM approximation (Eq. (39)) into Eq. (55), the following expression is obtained

[TABLE]

where $c_{\mathscr{D}p}$, $q_{\mathscr{D}p}$ and $\gamma_{ijkl_{\mathscr{D}p}}$ are the phase field, weighting functions and anisotropic tensor components evaluated at the material point $p$. Parameter $\mathscr{F}_{\mathscr{D}p}$ in Eq. (56) is expressed as

[TABLE]

where $l_{0_{\mathscr{D}p}}$, $k_{\mathscr{D}p}$, $\mathscr{H}_{\mathscr{D}p}$ and $\bar{\mathscr{G}}_{c_{\mathscr{D}p}}$ are the length scale parameter, model parameter, history field and critical fracture energy density of material point $\mathbf{x}_{p}$. Both $c_{\mathscr{D}p}$ and $q_{\mathscr{D}p}$ are interpolated at the nodal points of the background mesh, similarly to the case of the displacement field. Following the same procedure as in Section 4.1 the nodal residual vector for the phase field is established as

[TABLE]

where

[TABLE]

and

[TABLE]

Finally, applying the phase field interpolation and its spatial derivatives as

[TABLE]

[TABLE]

and

[TABLE]

Eq. (58) can be rewritten in the following convenient form as

[TABLE]

where $\bm{K}^{c}_{\mathscr{D}}$ is an $(N_{n}\times N_{n})$ coefficient matrix whose ${K}^{c}_{I,J}$ component is expressed as

[TABLE]

The $(N_{n}\times 1)$ vector $\bm{c}_{\mathscr{D}}$ holds the nodal values of the phase field, defined at the background grid nodes, and $\bm{F}^{c}_{\mathscr{D}}$ is the $(N_{n}\times 1)$ vector whose $F^{c}_{\mathscr{D}I}$ component is defined from Eq. (60).

4.3 Solution procedure

A staggered solution procedure [56] is employed to numerically solve the coupled Eqs. (53) and (64). In this, the two sets of equations are treated independently, by allowing the equation of motion to be solved either implicitly or explicitly [18]. Although an explicit time integration scheme is utilized herein to integrate Eqs. (53) in the time domain [34], an implicit time integration scheme can also be employed in a straightforward manner [60].

4.3.1 Explicit time integration scheme

To numerically solve the equation of motion (53), we employ the momentum formulation of the Material Point Method algorithm [61]. Eq. (53), is rewritten at the grid node $I$ at time $t$ as

[TABLE]

and considering, a forward Euler integration scheme, the acceleration field is expressed as

[TABLE]

where $\Delta t$ stands for the corresponding time step.

In view of Eq. (67), Eq. (66) is rewritten as

[TABLE]

where $\prescript{(t+\Delta t)}{}{\mathbf{p}}_{\mathscr{D}I}$ and $\prescript{(t)}{}{\mathbf{p}}_{\mathscr{D}I}$ are the nodal momentum at time $t+\Delta t$ and $t$, respectively.

At time $t$ the nodal momentums $\prescript{(t)}{}{\mathbf{p}}_{\mathscr{D}I}$ are unknown; hence, these are mapped from material points to grid node $I$ using Eq. (69)

[TABLE]

Similarly, the nodal internal forces $\prescript{(t)}{}{\bm{F}}^{int}_{\mathscr{D}I}$ are evaluated as

[TABLE]

Eq. (68) is numerically solved by extending a predictor-corrector algorithm introduced by [46] for granular media and further improved by [34] for the case of impact induced plasticity. In this, the trial momentums are initially evaluated for each discrete field $\mathscr{D}$, neglecting the contact forces $\prescript{(t)}{}{\bm{F}}^{cont}_{\mathscr{D}I}$, as

[TABLE]

The corresponding trial nodal velocities $\prescript{(t+\Delta t)}{}{\dot{\mathbf{u}}}^{trl}_{\mathscr{D}I}$ are then computed accordingly as

[TABLE]

The trial velocities correspond to the velocities of each discrete field $\mathscr{D}$ when no contact force is exerted between them.

The predicted trial velocities $\prescript{(t+\Delta t)}{}{\dot{\mathbf{u}}}^{trl}_{\mathscr{D}I}$ (evaluated from Eq. (72)) are then corrected according to Eq. (73)

[TABLE]

where $\prescript{(t+\Delta t)}{}{\dot{\mathbf{u}}}_{\mathscr{D}I}$ is the vector of corrected nodal velocities at time $t+\Delta t$. To evaluate the corrected nodal velocities using Eq. (73), the contact forces $\prescript{(t)}{}{\bm{F}}^{cont}_{\mathscr{D}I}$ must be evaluated first. The procedure for evaluating the contact forces between two discrete fields is presented in section 4.3.2.

4.3.2 Contact force evaluation

The contact force vector $\prescript{(t)}{}{\bm{F}}^{cont}_{\mathscr{D}I}$ is the sum of a normal $\prescript{(t)}{}{\bm{F}}^{nor}_{\mathscr{D}I}$ and a tangential $\prescript{(t)}{}{\bm{F}}^{tan}_{\mathscr{D}I}$ force vector. Hence, the corresponding components of these vectors, i.e. $\prescript{(t)}{}{F}^{nor}_{\mathscr{D}I}$ and $\prescript{(t)}{}{F}^{tan}_{\mathscr{D}I}$, should be initially computed taking into account the kinematic contact constraints presented in Eqs. (22) to (26) and (27) to (31). Their evaluation is performed through the following procedure.

The nodal centre of mass velocities are calculated using Eq. (74) below

[TABLE]

These correspond to the velocities that each discrete field $\mathscr{D}$ would have if these were to move as a single field (non-slip contact). The normal component of the contact force $\prescript{(t)}{}{F}^{nor,s}_{\mathscr{D}I}$ is evaluated considering the impenetrability condition defined in Eq. (25), at contact grid node $I$ as

[TABLE]

As aforementioned in section 3.1, when two bodies come into contact at contact grid node $I$ it holds that $\gamma_{nI}=0$. Substituting relation (73) into (75), considering the equilibrium of contact forces on the contact surface, i.e.,

[TABLE]

and also Eq. (74), the normal component of contact force is expressed as

[TABLE]

The surface unit normal vector is computed by using the mass gradients [43, see, e.g.,] as

[TABLE]

However, as also mentioned in [34], Eq. (78) should be modified to satisfy the collinearity conditions (22) and (27) at the contact surface $\partial\Omega_{\bar{f}}$ as

[TABLE]

to insure conservation of momentum.

To satisfy the non-tensional constraint (Eq. (24)) during contact, the normal component should be modified as

[TABLE]

Similarly, the tangential component of the contact force is evaluated considering the non-slip condition introduced in Eq. (30) as

[TABLE]

Substituting relation (73) into (81) and then making use of Eqs. (76) and (74), the tangential component of contact force is expressed as

[TABLE]

where the surface unit tangential vector $\mathbf{s}^{cont}_{\mathscr{D}I}$ can be derived as the unit vector that forms an orthogonal basis with $\mathbf{n}^{cont}_{\mathscr{D}I}$. The tangential component can be further modified to account for sliding at contact grid node $I$, considering the Coulomb friction model, as

[TABLE]

Therefore, the contact force is eventually evaluated as

[TABLE]

when the impenetrability condition

[TABLE]

is satisfied at contact grid node $I$.

Finally, once the contact force vector is computed from Eq. (84), the initially predicted nodal velocities $\prescript{(t+\Delta t)}{}{\dot{\mathbf{u}}}^{trl}_{\mathscr{D}I}$ should be corrected according to Eq. (73).

4.3.3 Material point properties update

The corrected nodal velocities $\prescript{(t+\Delta t)}{}{\dot{\mathbf{u}}}_{\mathscr{D}I}$ are utilized to update the material point properties. Hence, the total strains at the $p^{th}$ material point are evaluated as

[TABLE]

The total stresses are evaluated from Eq. (16). Finally, the displacement, velocity and acceleration of all material points are updated as

[TABLE]

[TABLE]

and

[TABLE]

respectively. The material point positions are also updated as

[TABLE]

4.3.4 Staggered solution algorithm

The solution procedure is summarized in Algorithm 1 where $E_{\mathscr{D}p}$ and $\nu_{\mathscr{D}p}$ are the Young’s modulus and Poisson ratio at the material points. At each time increment $\left(m=0,...,N_{steps}-1\right)$, the active part of the Eulerian grid is determined and the total number of active grid nodes $N_{n}$, unconstrained degrees of freedom $N_{dofs}$ and active cells $N_{cells}$ are evaluated (see [44] for details).

The background grid basis functions and their derivatives, i.e., $\bm{N}\left(\prescript{(m)}{}{\mathbf{x}}_{p}\right)$, $\nabla\bm{N}\left(\prescript{(m)}{}{\mathbf{x}}_{p}\right)$ and $\Delta\bm{N}\left(\prescript{(m)}{}{\mathbf{x}}_{p}\right)$ are evaluated at the material points with respect to the global coordinate system. To account for the arbitrary material orientation, the first and second spatial derivatives of the basis functions, i.e., $\nabla\bm{N}_{\phi_{p}}\left(\prescript{(m)}{}{\mathbf{x}}_{p}\right)$, and $\Delta\bm{N}_{\phi_{p}}\left(\prescript{(m)}{}{\mathbf{x}}_{p}\right)$, respectively are also evaluated in the local material coordinate system. In the 2D cases examined herein, the principal material orientation at the $p^{th}$ material point is defined with regards to the angle $\phi_{\mathscr{D}p}$ between the global axis $x_{1}$ and the principal material axis.

Following, the contact grid nodes are detected among the discrete fields according to Remark 4.

Remark 4

Two discrete fields are in contact at grid node $I$ when at least one material point from both discrete fields is projected into grid node $I$. In this case, the grid node $I$ is a contact grid node for this pair of discrete fields.

Next, the outward normal $\prescript{(m)}{}{\mathbf{n}}^{cont}_{\mathscr{D}I}$ and tangential $\prescript{(m)}{}{\mathbf{s}}^{cont}_{\mathscr{D}I}$ unit vectors are computed at the contact grid nodes. Mass, momentum and internal forces are projected from material points to grid nodes; thus, the quantities $\prescript{(m)}{}{M}^{u}_{\mathscr{D}I}$, $\prescript{(m)}{}{\mathbf{p}}_{\mathscr{D}I}$ and $\prescript{(m)}{}{\bm{F}}^{int}_{\mathscr{D}I}$ are obtained. Finally, the solution of the coupled Eqs. (68) and (64) is obtained within a set of $N_{staggs}$ staggered iterations $\left(k=1,...,N_{staggs}\right)$.

In the employed staggered scheme, the phase field Eq. (64) is initially solved for a specific value of the history field $\prescript{(m)}{}{\mathscr{H}}^{(k)}_{\mathscr{D}p}$. The basis functions $\bm{N}\left(\prescript{(m)}{}{\mathbf{x}}_{p}\right)$ and their spatial derivatives evaluated with respect to the material principal axes are utilized to compute the phase field coefficient matrix $\bm{K}^{c}_{\mathscr{D}}$ from relation (65). Thus, the phase field nodal values $\prescript{(m)}{}{c}^{(k)}_{\mathscr{D}I}$ are obtained for each discrete field $\mathscr{D}$. Next, the phase field nodal values are mapped back onto the material points and the degradation function $\prescript{(m)}{}{g}^{(k)}_{\mathscr{D}p}$ is computed at each material point. Next, the equation of motion (68) is integrated in time employing the predictor-corrector algorithm described in Section 4.3.1 and updated values for the history field $\prescript{(m)}{}{\mathscr{H}}^{(k)}_{\mathscr{D}p}$ are obtained.

Finally, the phase field nodal residual vector $\prescript{(m)}{}{R}^{c(k)}_{I}$ is evaluated according to the updated value of the history field $\prescript{(m)}{}{\mathscr{H}}^{(k)}_{\mathscr{D}p}$ and convergence is checked as $\|\prescript{(m)}{}{\bm{R}}^{c(k)}\|\leq tol_{c}$ or $k\geq N_{staggs}$ where $tol_{c}$ and $\|\cdot\|$ stand for the phase field tolerance value and the Euclidean norm, respectively. After convergence, the material point properties are updated and the algorithm proceeds to the next increment $m$.

Four conditions, namely C.1 to C.4, are also included in Algorithm 1. These are employed to verify that the kinematic constraints introduced in Eqs. (22)-(26) and (27)-(31) are satisfied at contact grid nodes.

5 Numerical examples

In this section, a set of two-dimensional numerical examples is presented. The numerical examples demonstrate the accuracy of the proposed PF-MPM against the standard PF-FEM as well as its computational efficiency in impact-fracture problems. Both the isotropic and anisotropic phase field models are examined within both single and multi discrete-field examples. Quadratic B-splines ( $C^{1}$ ) are utilized for the background grid. The initial cell density is chosen to be at least $3x3=9$ material points per cell element. Extensive numerical experiments performed in this work have demonstrated that this cell density results in accurate estimates when quadratic basis functions are utilized. Higher order B-splines are employed not only to compute the anisotropic phase field matrix in Eq. (65), but also to treat the “cell-crossing error” of the Material Point Method [62]. In all cases examined in this section, stability of the explicit integration scheme is established on the basis of the following upper bound for the time increment $\Delta t$

[TABLE]

where

[TABLE]

and $\Delta t_{cr}$ corresponds to the critical time step prescribed by the Courant – Friedrichs – Lewy (CFL) condition. Parameter $\alpha_{c}\in\left[0.8,0.98\right]$ in Eq. (91) depends on the nonlinearities of the system [31]. In all the numerical experiments, we consider $\alpha_{c}=0.80$.

In all numerical experiments presented the phase field residual tolerance was set to $tol_{c}=10^{-6}$ and a single stagger iteration was required for solution convergence. This is due to the small time step, imposed by the stability requirements (Eq. (91)).

5.1 Plate under impact loading

A plate under impact loading is examined. The same problem has been previously analysed by Borden et al. [18] with a finite element phase field implementation, considering a second order isotropic phase field formulation. The geometry and boundary conditions are presented in Fig. 5a. Herein, three cases are considered, i.e., (i) isotropic symmetry, (ii) cubic symmetry, and (iii) orthotropic symmetry. The material orientation is considered to be $\phi=+30^{o}$ with respect to $x$ axis (clockwise) as shown in Fig. 6a.

The cell (patch) spacing is $h=0.125$ mm and plane strain conditions are assumed. The grid is formed by two knot vectors $\Xi=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.001240\mathchar 44\relax\penalty 00.002481\mathchar 44\relax\penalty 0...\mathchar 44\relax\penalty 00.997518\mathchar 44\relax\penalty 00.998759\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$ and $H=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.003067\mathchar 44\relax\penalty 00.006134\mathchar 44\relax\penalty 0...\mathchar 44\relax\penalty 00.993865\mathchar 44\relax\penalty 00.996932\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$, $265024$ control points and $806x326=262756$ cells. The total number of material points is $2304000$ and the elastic material properties are $E=32000$ N/mm2, $\nu=0.20$ and $\rho=2450$ kg/m3.

The length scale parameter is chosen to be $l_{0}=0.25$ mm and $k_{f}=0.00$. In case (i), all the anisotropic material parameters are chosen such that $\gamma_{ijkl}=0$. Hence, the anisotropic phase field model reduces to the second order isotropic case. The maximum and minimum surface energy densities are equal to $\mathscr{G}_{c}\left(\theta\right)=\bar{\mathscr{G}}_{c}=\mathscr{G}_{c_{max}}=\mathscr{G}_{c_{min}}=0.003$ N/mm. In case (ii) cubic symmetry of the surface energy density is considered with $\bar{\mathscr{G}}_{c}=0.002121$ N/mm and anisotropic parameters $\gamma_{1111}=\gamma_{2222}=1.00$, $\gamma_{1122}=0.00$ and $\gamma_{1212}=74.00$. These parameters result into maximum and minimum surface energy densities $\mathscr{G}_{c_{max}}=0.0049$ N/mm and $\mathscr{G}_{c_{min}}=0.003$ N/mm, respectively. In case (iii) the anisotropic parameter $\gamma_{2222}$ is increased to $\gamma_{2222}=80.00$ giving rise to orthotropic symmetry with maximum and minimum surface energy densities $\mathscr{G}_{c_{max}}=0.0067$ N/mm and $\mathscr{G}_{c_{min}}=0.003$ N/mm, respectively.

The surface energy densities and their reciprocals for material orientation $\phi=+30^{o}$ are shown in Fig. 5b and Fig. 5c, respectively. In cases (ii) and (iii), the parameter $\bar{\mathscr{G}}_{c}$ is chosen so that $\mathscr{G}_{c_{min}}=0.003$ N/mm and to facilitate comparisons between all cases (see also, Figs. 5b and 5c).

A single discrete field is considered in this example. The solution procedure is implemented with a time step $\Delta t=0.025$ $\mu$s for $N_{steps}=3200$ steps. The critical time step is $\tilde{\Delta t_{cr}}=0.026$ $\mu$s. The traction is considered to be constant $\sigma=1$ N/mm2 during the analysis. The initial crack is modelled by introducing an initial history field at the corresponding material points as in [18]. The Rayleigh wave speed is $\dot{u}_{R}=2125$ m/s for the material parameters of that specimen [63].

5.1.1 Case (i): Isotropy

Initially, the PF-MPM is compared against the Phase Field Finite Element Method (PF-FEM) with the results obtained in [18] for the same cell (patch) spacing $h=0.125$ mm.

The total energy time-histories for the two solutions are shown in Figs. 6a and 6b where the 2 methods demonstrate a very good agreement. The total fracture energy results (see Fig. 6b) are in perfect agreement with the results reported in [18]. Minor differences are observed, especially for time $t>50$ $\mu$s. The total elastic strain energies (see Fig. 10a) also demonstrate very good agreement with some differences after $t>30$ $\mu$s.

The evolution of the phase field is presented in Fig. 7 for specific time steps. In these, the occurrence of a branched crack is observed at approximately $t=35$ $\mu$s. The evolution of the hydrostatic stress for the same timesteps is shown in Fig. 8. To demonstrate the influence of the surface energy density into the resulting crack paths, the reciprocal of the surface energy density is also plotted (black circle); the hydrostatic stresses are also shown in the same figure. Since the surface energy density is isotropic, hence independent of the material orientation, the crack naturally propagates along the vertical axis (see Figs. 7b and 8b) until branching occurs. Crack branching is perfectly symmetrical due to structure, load symmetry and the isotropic phase field model.

The crack tip velocities for the two methods are presented in Fig. 6c. As already mentioned in numerous works (see [18], [64], [65]), the crack tip and the exact location of crack branching cannot be identified uniquely due to the smooth description of the crack. Therefore, the crack tip velocity is measured with the methodology employed in [18] to facilitate verification. The results of both methods illustrate very good agreement. The crack widening and branching regions are almost the same for the two solutions and they are also shown in Fig. 6c. Crack widening here refers to the broadening of the damage zone prior to branching in accordance with the definition introduced in [18].

In Fig. 6c, the crack tip velocity is clearly below the Rayleigh wave speed which stands for the crack speed limit as elaborated by [63] [2]. However, experimental studies have shown that cracks rarely propagate at speeds close to the Rayleigh wave speed. In fact, they propagate at a fraction of the Rayleigh wave speed, i.e. 60% $\dot{u}_{R}$ [1]. As shown in Fig. 6c, the resulting crack tip velocities are below this limit.

5.1.2 Case (ii): Cubic Symmetry

The evolution of the phase field and the hydrostatic stress for specific time steps is presented in Fig. 9 and 10, respectively. From Figs. 9b and 10b, it is observed that due to the anisotropic cubic model and the material orientation $\phi=+30^{o}$ the crack does not initiate along the vertical axis. The crack propagates until the crack branches at approximately $t=50$ $\mu$s. In this case, the branched crack is not symmetrical and it branches along its two preferential weak directions (see Figs. 10c and 10d). To further illustrate this, the reciprocal of the surface energy density is also plotted in Fig. 10.

5.1.3 Case (iii): Orthotropic Symmetry

The evolution of the phase field and the hydrostatic stresses for several time steps in case (iii) are shown in Figs. 11 and 12, respectively. Contrary to the isotropic case, crack initiation does not occur along the the vertical axis as shown in Figs. 11b and 12b. Indeed, the anisotropic orthotropic model and the material orientation trigger the crack to propagate along a weak direction that is not aligned with vertical axis, similar to case (ii). However, contrary to case (ii) no crack branching is observed in case (iii). This can be justified by the fact that there is only one preferential weak direction. As a result, the crack path continues to propagate at one half of the plate. To further illustrate the effect of anisotropy on the resulting crack path, the reciprocal of the surface energy density is also plotted in Fig. 12.

The energy time-histories for all cases are shown in Figs. 13a and 13b. In all cases the crack initiates when the total elastic energy becomes approximately equal to $0.12$ J/m. It should be stressed that the total elastic strain energy evolves in an almost identical fashion in cases (ii) and (iii) until approximately $t=50$ $\mu$s. After that, the two models diverge as a result of the crack branching in the case of cubic symmetry only.

The crack tip velocities for the three cases are shown in Fig. 13c. The branched regions of both case (i) (isotropic) and (ii) (cubic) are also highlighted in Fig. 13c. The crack branching in isotropic symmetry is observed earlier than in cubic symmetry with a corresponding decrease in crack speed. In all cases, the crack initiates at approximately $t=10$ $\mu$s; the resulting crack tip velocities in all cases are comparable. This can be justified by the fact that the impact energy imposed as well as $\mathscr{G}_{c_{min}}$ are identical in all cases. In the orthotropic case, all the results are below the 60% $\dot{u}_{R}$ limit.

5.1.4 Crack branching and merging

The efficiency of phase field models to deal with complex crack paths, i.e. including crack branching and crack merging is demonstrated herein. The traction is increased to $\sigma=2.3$ N/mm2. All other model parameters are kept constant. The second order isotropic phase field model is used. The total duration of the analysis is $t_{tot}=130$ $\mu$s within $N_{steps}=5200$ steps. The evolution of the phase field for several time steps is presented in Fig. 14. In particular, Fig. 14b illustrates 5 branched cracks, i.e., 1 main, 2 secondary and 2 tertiary branches. More branched cracks are observed in Fig. 14c while in Fig. 14d 4 merged cracks are presented. The total strain energy together with the total fracture energy is shown in Fig. 16 whereas The evolution of the hydrostatic stress for that case is shown in Fig. 15 for several time steps

The crack tip velocity for that case is shown in Fig. 16c. The crack tip velocity is measured along the paths $\mathbf{C}_{3,1}$, $\mathbf{C}_{3,2}$, $\mathbf{C}_{3,3}$ and $\mathbf{C}_{3,4}$ that are marked in Fig. 14b. The increased impact loading, i.e. $\sigma=2.3$ N/mm2, leads to a crack initiation at approximately $t=5$ $\mu$s. This occurs earlier than in case (i), where the first crack initiates at approximately $t=10$ $\mu$s. The crack tip rapidly accelerates to the 60% $\dot{u}_{R}$ limit. Although, some points exceed the 60% $\dot{u}_{R}$ limit, the majority of measured points satisfy this condition while all points are clearly below the Rayleigh wave speed. In particular, the crack tip propagates with an average speed close to the 60% $\dot{u}_{R}$ limit. The three branched regions are also illustrated in Fig. 16c where a decrease in crack speed is observed.

5.2 Collision of two rings

Herein, the collision of two rings is analysed with the geometry and boundary conditions of the problem shown in Fig. 17a. The aim of this example is to demonstrate the robustness of the proposed method into resolving fragmentation problems and the interactions occurring between fragments.

The cell (patch) spacing is chosen to be $h=0.50$ mm for the numerical implementation and plane stress conditions are assumed with thickness $2$ mm. The grid is formed by two knot vectors $\Xi=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.00263\mathchar 44\relax\penalty 00.0052\mathchar 44\relax\penalty 0...\mathchar 44\relax\penalty 00.9947\mathchar 44\relax\penalty 00.9973\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$ and $H=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.00357\mathchar 44\relax\penalty 00.0071\mathchar 44\relax\penalty 0...\mathchar 44\relax\penalty 00.9928\mathchar 44\relax\penalty 00.9964\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$, $107724$ control points and $380x280=106400$ cells. Two discrete fields are considered, i.e., field A (left ring) and B (right ring). The corresponding friction coefficient is $\mu_{f}=0.65$. The total number of material points is $325620$.

The elastic material parameters are chosen to be $E=190000$ N/mm2, $\nu=0.30$ and $\rho=8000$ kg/m3 for both bodies. A time step $\Delta t=0.0125$ $\mu$s for $N_{steps}=50000$ is considered. The critical time step is $\tilde{\Delta t_{cr}}=0.071$ $\mu$s. The initial distance between the two rings is assumed to be $2h=1.00$ mm. An initial velocity is applied to the material points of the two rings as $\dot{\mathbf{u}}_{Ap_{(0)}}=\dot{\mathbf{u}}_{{(0)}}$ and $\dot{\mathbf{u}}_{Bp_{(0)}}=-\dot{\mathbf{u}}_{{(0)}}$. To examine the influence of the initial velocity into the resulting crack paths, two cases are considered, namely (i) $\dot{\mathbf{u}}_{{(0)}}=0.01$ mm/$\mu$s and (ii) $\dot{\mathbf{u}}_{{(0)}}=0.02$ mm/$\mu$s, respectively. The second order isotropic model is utilized for that problem, therefore $\gamma_{ijkl}=0$ with length scale parameter $l_{0}=1.00$ mm, $k_{f}=0.00$ and $\mathscr{G}_{c}\left(\theta\right)=\bar{\mathscr{G}}_{c}=\mathscr{G}_{c_{max}}=\mathscr{G}_{c_{min}}=6.00$ N/mm.

The total fracture energy time-history for both cases is shown in Fig. 17b. The evolution of the phase field and the hydrostatic stress for points ((1)-(6)) labelled in Fig. 17b is shown in Figs. 18 and 19 for case (i) and in Figs. 20 and 21 for case (ii), respectively.

5.2.1 Case (i): $\dot{\mathbf{u}}_{{(0)}}=0.01$ mm/$\mu$s

In case (i), a crack initiates at the contact surface of the two rings due to their initial impact (see Fig. 18b) followed by a second crack that initiates and fully propagates on the opposite side of each ring (see Fig. 18c and Fig. 18d, respectively). Material degradation also occurs on the top and bottom surfaces of each ring. Crack opening gradually increases (see Fig. 18e) and eventually both rings are fully separated in two fragments (see Fig. 18f). As also shown in Fig. 17b the fracture process has been fully developed by point (4), hence the fracture energy remains constant along the path (4)-(6).

5.2.2 Case (ii): $\dot{\mathbf{u}}_{{(0)}}=0.02$ mm/$\mu$s

Similar to case (i), in case (ii) a crack initiates at the contact surface of the two rings due to their initial impact (see Fig. 20b). Next, and contrary to case (i), two additional cracks simultaneously propagate of the top right (left) and bottom right (left) of each ring (see Fig. 20c). This is due to the increased impact velocity compared to case (i) where the corresponding points underwent material degradation only. Two more cracks are observed on the top left (right) and bottom left (right) at each ring (see Fig. 20d). The complete crack paths are presented in Fig. 20e. After that point, the fracture energy remains constant; existing cracks do not propagate and new cracks are not initiated.

The final deformed configuration of the problem is shown in Fig. 20f where each ring is split into five fragments. The PF-MPM method naturally resolves the large displacement motion of the fragments, accounting also for the non-stationarity of the contact surfaces (see, also, Figs. 21d, 21e and 21f). Using a phase field driven fracture approximation allows both, the crack paths and the contact surfaces to not be tracked algorithmically during the simulation process. Furthermore, this is accomplished with no mesh distortion induced errors, contrary to a FEM based approach.

5.3 Sphere-beam impact fracture problem

In this case, a sphere to beam impact fracture problem is examined. The geometry and boundary conditions of the problem are presented in Fig. 22. To examine the dependence of the resulting crack patterns and overall response of the beam on the level of assumed material anisotropy, three cases are examined, namely case (i), case (ii) and case (iii) with different surface energy densities. The material orientation is considered to be $\phi=+45^{o}$ with respect to $x$ axis (clockwise) for all cases.

The cell (patch) spacing is $h=0.125$ mm and plane stress conditions are assumed with thickness $2$ mm. The grid is formed by two knot vectors $\Xi=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.0025\mathchar 44\relax\penalty 00.0050\mathchar 44\relax\penalty 0...\mathchar 44\relax\penalty 00.9950\mathchar 44\relax\penalty 00.9975\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$ and $H=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.00167\mathchar 44\relax\penalty 00.00333\mathchar 44\relax\penalty 0...\mathchar 44\relax\penalty 00.99667\mathchar 44\relax\penalty 00.99833\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$, $242004$ control points and $400x600=240000$ cells. Two discrete fields are considered in this example, namely A for the sphere and B for the beam with the corresponding friction coefficient being $\mu_{f}=0.65$. The total number of material points is $536796$.

An initial velocity is applied to all material points in the sphere $\dot{\mathbf{u}}_{Ap_{(0)}}=0.02$ mm/$\mu$s while the beam is at rest at this stage. The initial distance between the sphere and the beam is considered to be $h=0.125$ mm. The elastic material parameters are chosen to be $E=190000$ N/mm2, $\nu=0.30$ and $\rho=8000$ kg/m3 for both bodies. The solution procedure is implemented with a time step $\Delta t=0.0125$ $\mu$s. The critical time step is $\tilde{\Delta t_{cr}}=0.018$ $\mu$s.

In case (i), the second order isotropic phase field model (i.e. $\gamma_{ijkl}=0$) is chosen with surface energy density $\mathscr{G}_{c}\left(\theta\right)=\bar{\mathscr{G}}_{c}=\mathscr{G}_{c_{max}}=\mathscr{G}_{c_{min}}=10.6066$ N/mm for the beam. In case (ii) , the second order isotropic phase field model is chosen again, but with a reduced surface energy density $\mathscr{G}_{c}\left(\theta\right)=\bar{\mathscr{G}}_{c}=\mathscr{G}_{c_{max}}=\mathscr{G}_{c_{min}}=9.75$ N/mm for the beam. Finally, in case (iii) the fourth order orthotropic model is utilized with anisotropic parameters $\bar{\mathscr{G}}_{c}=7.50$ N/mm, $\gamma_{1111}=80.00$, $\gamma_{2222}=1.00$, $\gamma_{1122}=0.00$ and $\gamma_{1212}=74.00$. These parameters result in maximum and minimum surface energy densities $\mathscr{G}_{c_{max}}=23.6892$ N/mm and $\mathscr{G}_{c_{min}}=10.6066$ N/mm, respectively for the beam. The surface energy density of the sphere is taken us sufficiently large (i.e. $\mathscr{G}_{c_{A}}\left(\theta\right)=100\mathscr{G}_{c_{B}}\left(\theta\right)$) so that the sphere remains undamaged in all cases. The length scale parameter is $l_{0}=0.25$ mm and $k_{f}=0.00$ in all cases. In all cases reported in this section, the projectile does not penetrate the beam, rather it bounces back and the beam undergoes free vibrations.

5.3.1 Case (i): Isotropy - $\mathscr{G}_{c}\left(\theta\right)=10.6066$ N/mm

The time history of the total fracture energy is shown in Fig. 23a. In Fig. 23a, the path points (1-7) are labelled to facilitate discussion on the material response. Phase field and hydrostatic stress snapshots corresponding to points (1-6) are shown in Figs. 24 and 25, respectively.

The sphere initially comes into contact with the beam and fracture initiates at the contact surface (see. Fig. 24b and point (2) in Fig. 23a). Next, the right edge of the beam gradually degrades (see. 24c and point (3) in Fig. 23a) just before a crack initiates at the middle right-edge point. However, as the beam vibrates, the degradation continues at the left edge of the beam (see Fig. 24d and point (4) in Fig. 23a) a median crack develops and propagates just below the crack nucleation region (see Fig. 24e and point (5) in Fig. 23a). The complete crack path is shown in Fig. 24f.

The results of Fig. 23a can be further examined in view of the total fracture energy evolution. The evolution of the total fracture energy from point (1) to (2) corresponds to damage initiating between the sphere and the beam at their contact surface. Degradation at the right edge of the beam results in a further increase of the fracture energy corresponding to the path (2-3). Finally, the crack rapidly propagates from point (4) to (6). Hence, the total fracture energy corresponding to crack propagation is

[TABLE]

This is in very good agreement with the analytical prediction as $A_{f}\cdot\mathscr{G}_{c}\left(\theta\right)=10\cdot 2\cdot 10.6066=212.13$ mJ, where $A_{f}$ stands for the fracture surface. The slight increase of the total fracture energy from point (6) to (7) corresponds to the marginal degradation of the beam material during the free vibration regime of its response.

5.3.2 Case (ii): Isotropy - $\mathscr{G}_{c}\left(\theta\right)=9.75$ N/mm

Even though the variation in $\mathscr{G}_{c}$ is small compared to case (i), it results in a significantly different material response. The total fracture energy time-history for case (ii) is shown in Fig. 23b. The evolution of the phase field and the hydrostatic stress for points ((1)-(6)) labelled in Fig. 23b is shown in Figs. 26 and 27, respectively.

Similar to case (i), the sphere initially comes into contact with the beam and causes damage at their contact surface. As a result, material degradation is observed at the right edge of the beam (see. Fig. 26b and point (2) in Fig. 23b) as in case (i). Contrary to case (i) however, a flexural crack initiates at the middle right-edge point of the beam due to maximum principal tensile stresses developing at the tensile fibre of the beam.

As the beam oscillates the maximum tension region alternates between the two edges and the crack arrests (see Fig. 26c and point (3) in Fig. 23b). A second crack then initiates at the left edge (see Fig. 26d and point (4) in Fig. 23b) and propagates (see Fig. 26e and point (5) in Fig. 23b) until the two cracks finally merge as shown in Fig. 26f.

As in case (i), the evolution of the fracture energy (shown shown in Fig. 23b) is consistent with the observed response. The first crack (right crack) initiates at point (2) and stops at point (3). The second crack (left crack) propagates from point (4) to (6). Therefore, the total fracture energy is

[TABLE]

This is again in very good agreement with the analytical prediction as in this case $A_{f}\cdot\mathscr{G}_{c}\left(\theta\right)=10\cdot 2\cdot 9.75=195$ mJ.

5.3.3 Case (iii): Orthotropy

Orthotropic anisotropy with a material orientation $\phi=+45^{o}$ results in two cracks at each beam edge (right and left) that do not coincide with the horizontal axis as in case (ii). The evolution of the phase field and the hydrostatic stress are represented for several time steps in Figs. 28 and 29, respectively. The characteristic points ((1)-(7)) of that analysis are shown in Fig. 23c.

Similarly to the previous cases, damage initiation is observed at the contact surface (see Fig. 28b and point (2) in Fig. 23c). Next, degradation occurs at the left edge of the beam (see Fig. 28c and point (3) in Fig. 23c). The first crack (right crack) initiates at middle right-edge point of the beam and propagates along the material’s week direction until it arrests in the vicinity of the beam’s neutral axis (see Fig. 28d and point (4) in Fig. 23c). After impact, further degradation occurs due to the beam’s free vibration resulting in degradation to its left edge (see Fig. 28e and point (5) in Fig. 23c). Finally, a second crack (left crack) initiates at the middle left-edge of the beam and propagates along the material’s week direction (see Fig. 28f and point (6) in Fig. 23c). Similar to the first crack, the second crack arrests in the vicinity of the beam’s neutral axis. The final crack paths are shown in Fig. 29f where the two cracks do not merge as in case (ii).

5.3.4 Discussion on observed fracture patterns

Fracture patterns emerging from static indentation of a practically rigid sphere against a deformable solid as well as from low and high speed impact tests have been the focus of extensive experimental investigations, see, e.g., [66, 67]. The failure modes observed vary considerably with the velocity of the projectile, the flexibility of the impacted beam, and the interface properties [67].

We focus here in case (i) with isotropic fracture energy equal to $\mathscr{G}_{c}\left(\theta\right)=10.6066$ N/mm. The fracture patterns shown in Fig. 24 correspond to a median type of crack, with the crack at the left end nucleating due to impact and then propagating towards the right edge driven by the principal tensile stresses at the mid-span. To investigate the effect of the projectile velocity on the induced fracture pattern, a total of 33 analysis cases is performed keeping the geometry, the elastic and the fracture properties of the beam similar to those reported in section 5.3. In each case the projectile velocity is varied from 0.02 mm/$\mu$s to 0.18 mm/$\mu$s at a step size of 0.005 mm/$\mu$s. All analysis parameters, the background cell size and the cell density are similar to section 5.3.

We define the Hertzian cone index $\beta$ with a value $\beta=0$ corresponding to a cone not developing and $\beta=1$ when a cone develops. This is plotted versus the projectile velocity in Fig. 30; a cone fracture pattern occurs for velocities larger than 0.10 mm/$\mu$s. Conversely, for velocities smaller than 0.10 mm/$\mu$s the crack pattern is consistent with the flexure failure mode described in section 5.3. In Fig. 31, phase field snapshots are shown for six particular cases of projectile velocity. All snapshots correspond to time $t=16\mu$ when the maximum value of fracture energy for all case has been attained.

For the case of the lowest impact velocity considered in Fig. 31a, the crack pattern involves a median crack propagating from the left and towards the right edge of the beam. Secondary, flexural cracks appear at the right edge of the beam. Such a response is consistent with experimental observations on brittle materials at low impact loads where a plastic band initiates at the impact zone prior to crack formation, see, e.g., [68]. In the framework presented herein, the material degradation prior to fracture assumes this role.

Increasing impact velocities result in a Hertz cone formation at the vicinity of the impact zone. Secondary cracks also propagate from the left edge. Of interest is also the evidently smooth transition from a median to a Hertz cone fracture pattern from Fig. 31c to Fig. 31d.

5.4 Anisotropic plate with centred crack

The case of the anisotropic rectangular plate shown in Fig. 32 is examined, made from a unidirectional HTA/ 6376 composite laminate and subjected to an initial velocity field $\dot{u}\left(\mathbf{x}\right)_{(0)}=0.0002y$ mm/$\mu$s. The plate contains a pre-existing crack at its centre with length $25$ mm. The material properties of the composite are summarized in Table 5.4. In addition, the following values hold, i.e., $\nu_{xy}$=0.3 and $G_{xy}=5500$ N/mm2.

In their experimental campaign, Cahil et al. [69] have shown that cracks grow parallel to the fibre direction hence indicating that the damage originates only through matrix failure. We consider herein the case of fibre orientation at $\phi=+45^{o}$ as shown in Fig. 32. This allows us to use the same Young’s modulus and Poisson’s ratio along $x$ and $y$. The elastic material properties considered are $E=14980$ N/mm2 and $\nu=0.36$.The mass density is $\rho=1586$ kg/m3. The length scale parameter is considered $l_{0}=1$ mm and the anisotropic parameters are taken to be $\gamma_{1111}=1.00$, $\gamma_{2222}=2900$, $\gamma_{1122}=0.00$ and $\gamma_{1222}=74.00$ and $\bar{\mathscr{G}}_{c}=4.175$ N/mm. These parameters correspond to an orthotropic surface energy with $\mathscr{G}_{c_{min}}=5.9067$ N/mm along the fibre orientation and $\mathscr{G}_{c_{max}}=30.9044$ N/mm normal to the fibre.

The grid is formed by two knot vectors $\Xi=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.00666\mathchar 44\relax\penalty 00.01333\mathchar 44\relax\penalty 0...\mathchar 44\relax\penalty 00.98666\mathchar 44\relax\penalty 00.99333\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$ and $H=\{0\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00\mathchar 44\relax\penalty 00.0033\mathchar 44\relax\penalty 00.0066\mathchar 44\relax\penalty 00.01\mathchar 44\relax\penalty 0...0.9899\mathchar 44\relax\penalty 00.9933\mathchar 44\relax\penalty 00.9966\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\mathchar 44\relax\penalty 01\}$, $45904$ control points and $150x300=45000$ cells. The cell (patch) spacing is $h=1.00$ mm and plane stress conditions are assumed. The total number of material points is $281250$. The solution procedure is implemented with a time step $\Delta t=0.0125$ $\mu$s for a total time of $25$ $\mu$s. The critical time step is $\tilde{\Delta t_{cr}}=0.201$ $\mu$s.

The results for the numerical simulations together with the experimental observations are shown in Fig. 33. The reciprocal of the surface energy density (black eclipse) is also plotted on these snapshots. In Figs. 33b and 33c, the phase field evolves along the material orientation $\phi=+45^{o}$. The crack paths derived from our simulation agree well with the experimental crack paths presented in Fig. 33d.

6 Conclusions

In this work, a novel numerical method is introduced for the treatment of dynamic brittle fracture in both isotropic and anisotropic media. The evolution of crack paths is represented by means of phase field models within a Material Point Method setting. Anisotropy is explicitly introduced in the fracture energy through a crack density functional. The method is further extended to account for frictional contact problems involving phase field fracture adopting a discrete field approach. A notable advantage of the proposed formulation is that both the equilibrium and phase field governing equations are solved independently for each discrete field rendering the method suitable for parallel implementation.

The method is rather appealing for the case of phase field modelling where very fine meshes are commonly required due to the regularized crack topology. Rather than employing a uniform background mesh and material point density, multiple small scale problems can be solved separately for each discrete field at their corresponding background domain. In terms of contact driven fracture, contrary to standard FEM implementations that necessitate the algorithmic treatment of local contact features, these now naturally emerge from the interaction of material points within a fixed Eulerian mesh. Indeed, the fixed Eulerian grid is utilized to identify the contact surfaces using the material points’ projection on the grid.

A set of representative numerical examples is presented where the computational advantages of PF-MPM are demonstrated. The method is verified against the standard Phase Field Finite Element Method; the two methods are in good agreement. The influence of anisotropy is examined in terms of crack path, time history energy results and crack tip velocities. Benchmark problems with complex crack path i.e. crack branching and merging are considered and the robustness of the method is established. It is shown that different loading velocities and fracture material parameters strongly influence the dynamic failure response of the structure and the resulting crack paths. Fracture energies computed from the proposed method are compared and indeed verified against the corresponding analytical predictions. Finally, crack paths derived from the method are validated against experimental observations.

Acknowledgement

The research described in this paper has been financed by the University of Nottingham through the Dean of Engineering Prize, a scheme for pump priming support for early career academic staff. The authors are grateful to the University of Nottingham for access to its high performance computing facility.