Numerical Modelling of Crack Initiation, Propagation and Branching under Dynamic Loading

TL;DR

This paper presents a novel pseudo-spring SPH method to simulate crack initiation, propagation, and branching under dynamic loading in 2D and 3D, showing good agreement with experimental results.

Contribution

The paper introduces a pseudo-spring analogy in SPH for damage modeling, enabling accurate simulation of complex crack behaviors in dynamic conditions.

Findings

Crack speeds and paths match experimental data

Crack branching and surface evolution are effectively modeled

Framework successfully simulates 3D crack phenomena

Abstract

In this paper crack initiation, propagation and branching phenomena are simulated using the Pseudo-Spring Smooth Particle Hydrodynamics (SPH) in two and three-dimensional domains. The pseudo-spring analogy is used to model material damage. Here, the interaction of particles is limited to its initial immediate neighbours. The particles are connected via springs. These springs do not provide any extra stiffness in the system but only define the level of interaction between the connecting pairs. It is assumed that a crack has passed through a spring connecting a particle pair if the damage indicator of that spring becomes more than a predefined value. The crack branching of a pre-notched plate under dynamic loading and the effect of loading amplitude are studied here. The computed crack speeds, crack paths and surfaces are compared with experimental and numerical results available in the literature and are found to be in good agreement. The effect of notch location for a plate with a circular hole is studied here. The ability of the framework to model arbitrary crack paths and surfaces are demonstrated via three-dimensional simulations of chalk under torsion, Kalthoff-Winkler experiment and Taylor bullet impact.

Tables

Table 1. Table 1 : Parameters for pre-notched plate under tensile loading

	Material Proerties				Discretization		Artificial Viscosity
Parameter	$\rho$	$E$	$\nu$	${ϵ}_{max}$	$\Delta p$	$h$	${\beta }_1$	${\beta }_2$
	($kg/m^3$)	(GPa)			(mm)	(mm)
Value	2450	32	0.2	0.000509	0.125	0.250	1.0	1.0

Table 2. Table 2 : Parameters for pre-notched plate with a circular hole under tensile stress

	Material Proerties				Discretization		Artificial Viscosity
Parameter	$\rho$	$E$	$\nu$	${ϵ}_{max}$	$\Delta p$	$h$	${\beta }_1$	${\beta }_2$
	($kg/m^3$)	(GPa)			(mm)	(mm)
Value	2700	71.4	0.25	0.00483	0.1	0.15	1.0	1.0

Table 3. Table 3 : Parameters for the chalk under torsional loading

	Material Proerties				Discretization		Artificial Viscosity
Parameter	$\rho$	$E$	$\nu$	${\sigma }_{max}^p$	$\Delta p$	$h$	${\beta }_1$	${\beta }_2$
	($kg/m^3$)	(GPa)		(MPa)	(mm)	(mm)
Value	1150	2	0.18	15	0.4	0.52	1.0	1.0

Table 4. Table 4 : Parameters for Klathoff-Winkler experiment simulation

	Material Proerties				Discretization		Artificial Viscosity
Parameter	$\rho$	$E$	$\nu$	${\delta }_{tc}$	$\Delta p$	$h$	${\beta }_1$	${\beta }_2$
	($kg/m^3$)	(GPa)			(mm)	(mm)
Value	8000	190	0.3	0.0044	0.8	1.04	1.0	1.0

Table 5. Table 5 : Computational parameters used for Taylor bullet impact

Inter-Particle	Smoothing	Time-step	Artificial Viscosity
Spacing (p)	Length (h)	($\Delta$t)	Parameters (${\beta }_1,{\beta }_2$)
0.3 mm	0.6 mm	$5\times {10}^{-9}$s	(1.0,1.0)

Equations

$$\frac{d\rho }{dt}=-\rho \frac{\partial u^{\beta }}{\partial x^{\beta }}$$

$$\frac{d{\rho }_i}{dt}=\sum _{j\in {\bar{N}}^i}{m}_j{u}_{ij}^{\beta }{\bar{W}}_{ij,\beta }$$

$$\frac{d{u}^{\alpha }}{dt}=\frac{1}{\rho }\frac{\partial {\sigma }^{\alpha \beta }}{\partial x^{\beta }}$$

$$\frac{d{u}_i^{\alpha }}{dt}=\sum _{j\in {\bar{N}}^i}{m}_j\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}-{\pi }_{ij}{\delta }^{\alpha \beta }-P_{ij}^a{\delta }^{\alpha \beta }\right){\bar{W}}_{ij,\beta }$$

$$\frac{de}{dt}=\frac{{\sigma }^{\alpha \beta }}{\rho }\frac{\partial u^{\alpha }}{\partial x^{\beta }}$$

$$\frac{d{e}_i}{dt}=-\frac{1}{2}\sum _{j\in {\bar{N}}^i}{m}_j{u}_{ij}^{\alpha }\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}-{\pi }_{ij}{\delta }^{\alpha \beta }-P_{ij}^a{\delta }^{\alpha \beta }\right){\bar{W}}_{ij,\beta }$$

$${\hat{W}}_{ij,\beta }=B^{\beta \alpha }{W}_{ij,\alpha }$$

$$\mathbf{B}={\mathbf{A}}^{-1},and{A}^{\beta \alpha }=-\sum _j\frac{m_j}{{\rho }_j}{x}_{ij}^{\beta }{W}_{ij,\alpha }$$

$$W(q,h)={\alpha }_d\{\begin{array}{ll}1-\frac{3}{2}{q}^2+\frac{3}{4}{q}^3,& \text{if }0\le q\le 1\\ \frac{1}{4}(2-q{)}^3,& \text{if }1\le q\le 2\\ 0& \text{otherwise}\end{array}$$

$${\pi }_{ij}=\{\begin{array}{ll}\frac{-{\beta }_1{\bar{C}}_{ij}{\mu }_{ij}+{\beta }_2{\mu }_{ij}^2}{{\bar{\rho }}_{ij}},& \text{if }{\mathbf{U}}_{ij}.{\mathbf{X}}_{ij}\le 0\\ 0,& \text{otherwise}\end{array}$$

$$P_{ij}^a=\gamma \left(\frac{|P_i|}{{\rho }_i^2}+\frac{|P_j|}{{\rho }_j^2}\right){\left[\frac{W(\Delta x,h)}{W(\Delta p,h)}\right]}^l$$

$$P=K\left(\frac{\rho }{{\rho }_0}-1\right)$$

$${\dot{S}}_{\alpha \beta }=2\mu \left({\dot{ϵ}}_{\alpha \beta }-\frac{1}{3}{\delta }^{\alpha \beta }{\dot{ϵ}}^{\gamma \gamma }\right)+S^{\alpha \gamma }{\omega }^{\beta \gamma }+S^{\gamma \beta }{\omega }^{\alpha \gamma }$$

$${\dot{ϵ}}^{\alpha \beta }=\frac{1}{2}\left(\frac{\partial u^{\alpha }}{\partial x^{\beta }}+\frac{\partial u^{\beta }}{\partial x^{\alpha }}\right)$$

$${\omega }^{\alpha \beta }=\frac{1}{2}\left(\frac{\partial u^{\alpha }}{\partial x^{\beta }}-\frac{\partial u^{\beta }}{\partial x^{\alpha }}\right)$$

$$\Delta {ϵ}_{pl}^{\alpha \beta }=\frac{1-c_f}{2\mu }{S}^{\alpha \beta }$$

$$\Delta {ϵ}_{pl}=\sqrt{\frac{2}{3}\Delta {ϵ}_{pl}^{\alpha \beta }\Delta {ϵ}_{pl}^{\alpha \beta }}=\frac{1-c_f}{3\mu }\sqrt{\frac{3}{2}{S}^{\alpha \beta }{S}^{\alpha \beta }}$$

$$\Delta W_p=\Delta {ϵ}_{pl}^{\alpha \beta }{S}_n^{\alpha \beta }$$

$$\frac{d{\rho }_i}{dt}=\sum _{j\in N_U^i}{m}_j{u}_{ij}^{\beta }{\bar{W}}_{ij,\beta }+\sum _{j\in N_D^i}{m}_j{u}_{ij}^{\beta }\left(f_{ij}{\bar{W}}_{ij,\beta }\right)$$

$$\begin{array}{r}\frac{d{u}_i^{\alpha }}{dt}=\sum _{j\in N_U^i}{m}_j\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}-{\pi }_{ij}{\delta }^{r\beta }-P_{ij}^a{\delta }^{r\beta }\right){\bar{W}}_{ij,\beta }\\ +\sum _{j\in N_D^i}{m}_j\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}-{\pi }_{ij}{\delta }^{r\beta }-P_{ij}^a{\delta }^{r\beta }\right)\left(f_{ij}{\bar{W}}_{ij,\beta }\right)\end{array}$$

$$\begin{array}{r}\frac{d{e}_i}{dt}=-\frac{1}{2}\sum _{j\in N_U^i}{m}_j{u}_{ij}^{\alpha }\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}-{\pi }_{ij}{\delta }^{\alpha \beta }-P_{ij}^a{\delta }^{\alpha \beta }\right){\bar{W}}_{ij,\beta }\\ -\frac{1}{2}\sum _{j\in N_D^i}{m}_j{u}_{ij}^{\alpha }\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}-{\pi }_{ij}{\delta }^{\alpha \beta }-P_{ij}^a{\delta }^{\alpha \beta }\right)\left(f_{ij}{\bar{W}}_{ij,\beta }\right)\end{array}$$

$$\frac{d{u}_i^{\alpha }}{dt}={\left(\frac{d{u}_i^{\alpha }}{dt}\right)}_{Equation\text{ref}eq6df}-\sum _{j\in N^i}{m}_j\left(\frac{{\sigma }_0^{\alpha \beta }}{{\rho }_i{\rho }_j}\right)W_{ij,\beta }$$

$$D=\{\begin{array}{ll}1,& ϵ\ge {ϵ}_{max}\\ 0,& Otherwise\end{array}$$

$$u=u_{max}\sin (\frac{\theta +\frac{\pi }{2}}{2})$$

$${\delta }_t=\frac{r_{ij}{|}_t-\Delta p}{\Delta p}$$

Full text

††journal: Engineering Fracture Mechanics

1 Introduction

Numerical modelling of crack initiation, propagation and branching under dynamic loading is still an open problem. The tracking of arbitrary paths and their propagation are modelled by different grid-based or mesh-free numerical techniques. The grid-based method such as Finite Element Method (FEM) or mesh-free methods such as Smooth Particle Hydrodynamics (SPH) cannot model crack discontinuity in their standard form. Each of these methods needs some special treatment to model the crack discontinuity.

FEM is used to model crack path and its propagation by element deletion or kill method [1], [2] or by adaptive re-meshing [3], [4], [5], [6]. In the element deletion or kill technique [1] [2], an element is assumed to fail once the stress, strain or damage values exceed a predefined value to model arbitrary cracks. The edges of the failed meshes act as the new crack surfaces. The main problem in element deletion or kill algorithms is the numerical instabilities arise due to excessive mesh distortion. This may lead to unreal element volume and error in the numerical computation [7]. Though FEM coupled with automatic and adaptive re-meshing is used for modelling of arbitrary crack propagation [3], [4], [5], [6], this technique is not quite useful for tracking multiple cracks. The method of re-meshing is also complicated as multiple field variables require internal mapping between different meshes during re-meshing.

The crack growth is also simulated through the element edges using the inter-element separation techniques [8], [9], [10], [11]. The crack path always propagates through the element edges or boundaries. Hence, the element size and shape are important parameters for tracking the crack path correctly. The effect of mesh size in computation is reduced in Embedded Finite Element Method (EFEM) [12] [13]. The EFEM is capable of tracking arbitrary crack propagation without any re-meshing. In this method an enrichment is introduced at the element level [14] [15]. Belytschko and his group [16], [17], [18] developed the extended finite element method (XFEM) for modelling of arbitrary crack propagation. They utilise the partition of unity [19] for local enrichment. In XFEM, only the surrounding elements of a crack and the crack propagation path (known a priori) are enriched and the crack propagated through these enriched elements (not only through the element boundaries). Though these element based methods can model evolving cracks and track their paths, the requirement of re-meshing or local enrichments makes these numerically and computationally intensive. Moreover, these methods may not be so effective for multiple cracks, crack interaction or excessive deformations. In this context, mesh-free methods are more suitable.

Different types of mesh-free/particle methods are introduced over the years. Belytschko and his colleagues [20], [21], [22], [23] developed the Element-Free Galerkin (FFG) method for static and dynamic crack propagation. As the method is mesh-free, there is no re-meshing. In their works, they used moving least square interpolants as test and trial functions. Rabczuk and Belytschko [24] [25] proposed another cracking particle method for modelling of crack path and its evolution. They introduce discontinuous enrichments, and the crack is modelled by particle splitting. As the cohesive traction at a particle reduces to zero, the particle is assumed to be cracked and split into two particles. Later, another version of cracking particle is proposed [26] without any enrichment. However, this still requires particle splitting. The need for discontinuous enrichments or particle splitting or visibility criteria [27] for crack path tracking makes these methods computationally intensive.

Another particle method, SPH is truly a mesh-free method [28], introduced for astronomical problems. This is extended to fluid problems and later on solids [29]. In SPH, the computational domain is discretised into a set of particles. At any time instant, a particle interacts with its neighbour particles within its influence domain (known as the neighbourhood) through a kernel function. The bell-shaped kernel function ensures that the level of interaction with the nearest neighbour particles is higher than the particles close to the neighbourhood boundary. The particles have no interaction with particles outside its influence domain which is fixed by the smoothing length of the kernel function. If a particle moves out of the neighbourhood of a particle, their interaction ceases. The kernel function is independent of material properties. Naturally, it is unable to compute the damage state of different particles. But, this does not pose any problem in astronomical or fluid problems as there is no concept of damage. However, for solid continua, it poses significant difficulty as large deformation may often lead to damage/failure. If a particle which moves out of the influence domain of another particle is considered to be damaged, the damage variable will be a function of the smoothing length/radius of the kernel function. As the smoothing length/radius of the kernel function increases, the influence domain of a particle increases. A particle has to move further to get out of the influence domain of another particle. Hence, the damage state of a particle will change if the kernel radius increases or decreases for material under the same loading conditions. This may lead to non-physical behaviour in the numerical computation.

In this context, a ‘pseudo-spring’ strategy in the SPH framework is developed [30] for crack modelling. Springs connect all the particles in the computational domain. The springs are termed as ‘pseudo-spring’ as these do not provide any extra stiffness to the system but only define the level of interaction between particles based on the damage state. The interaction of particles is not the same as the standard SPH. Here, particles only interact through standard kernel function with its immediate neighbours (from initial configuration). This framework does not require any crack enrichment, particle splitting or tracking of crack paths. When a certain parameter of a pseudo-spring is more than a predefined critical value, the spring is considered to be damaged, and the interaction ceases between the connected particles. Moreover, it is assumed that the crack has passed through that particle pair. For partially damage pseudo-springs, the level of interaction is decreased accordingly. In this framework, the crack path can be tracked automatically by the damaged pseudo-springs. The simplicity of the SPH is preserved in the computation. This method is successfully used to simulated impact induced damages [31], [32], [33], [34], [35].

In this work, this pseudo-spring analogy in SPH framework is used to simulate the crack initiation, propagation and branching in two and three dimensions. The propagation of cracks, the deviation of crack paths under different conditions and generation of multiple crack branches are simulated. The effects of tensile loading amplitude in crack branching are discussed. The generation of arbitrary crack surfaces under dynamic loading in three dimensions is modelled. The propagation of cracks under impact loading in three dimensions is also demonstrated. The ability of the pseudo-spring SPH to simulate the crack initiation, propagation and branching under dynamic loading in two and three-dimensional space is demonstrated via numerical examples. The computed crack speeds, their paths and crack surfaces are compared with experimental and numerical results available in the literature.

The paper is organised as follows: The necessary steps of pseudo-spring SPH framework and the boundary condition are discussed in Section 2. Numerical examples of two and three-dimensional crack initiation, propagation and branching under dynamic loading are discussed in Section 3. Some conclusions are drawn in Section 4.

2 Computational framework : Pseudo-Spring SPH

A set of particles represents the computational domain in SPH. Here, the basics of SPH and the ‘pseudo-spring’ analogy incorporated in the framework are discussed briefly. A more detailed description of standard SPH and the ‘pseudo-spring’ framework can be found in [29], [36], [37], [38], [30], [34], [35], and references therein. In SPH, a particle interacts with its surrounding neighbourhood particles through a kernel function. The mass, momentum and energy equations and their approximated particle forms are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, $i$ denotes the particle of interest, $j$ denotes its neighbour particles within its influence radius and $\mathbb{\bar{N}}^{i}\in\mathbb{Z^{+}}$ is the set of total number of neighbour particles within its influence domain. The density of the material is represented by $\rho$, $u^{\alpha}$ is the velocity component. $\sigma^{\alpha\beta}$ represents the Cauchy stress components. $e$ is the specific internal energy. $\alpha$, $\beta$ are the indices for coordinate axes. The time derivative is in moving Lagrangian frame and $u^{\alpha}_{ij}=u^{\alpha}_{i}-u^{\alpha}_{j}$. $\bar{W}_{ij,\beta}$ is the symmetric form of the kernel gradient correction computed as $\bar{W}_{ij,\beta}=0.5\left(\hat{W}_{ij,\beta}+\hat{W}_{ji,\beta}\right)$ [39]. The kernel gradient correction $\hat{W}_{ij,\beta}$ is calculated as [40]

[TABLE]

where,

[TABLE]

The kernel gradient is represented as $W_{ij,\beta}=\frac{\partial W(\bm{x}_{i}-\bm{x}_{j},h)}{\partial x^{\beta}_{i}}$. $h$ is the smoothing length of the kernel function. Here, the following cubic B-spline kernel is used.

[TABLE]

where, $\alpha_{d}=\frac{10}{7\pi h^{2}}$ and $\alpha_{d}=\frac{1}{\pi h^{3}}$ in 2D and 3D respectively. $q$ represents the normalised position vector as $q=\frac{|\bm{x}_{i}-\bm{x}_{j}|}{h}$.

The artificial viscosity $\pi_{ij}$ shown in Equation 10 is utilised in the current computation to remove any unphysical oscillation present near any jump or shock [28].

[TABLE]

The intensity of the artificial viscosity is determined using $\beta_{1}$, $\beta_{2}$. $\mu_{ij}=\frac{h(\bm{U}_{ij}.\bm{X}_{ij})}{\bm{X}^{2}_{ij}+\epsilon h^{2}}$ and $\epsilon=0.01$. The sound speed through the medium is computed as $C=\sqrt{\frac{E}{\rho}}$, $E$ denotes the Young’s modulus, $\bm{X}_{ij}=\bm{X}_{i}-\bm{X}_{j}$, $\bm{U}_{ij}=\bm{U}_{i}-\bm{U}_{j}$, $\bar{C}_{ij}=0.5(C_{i}+C_{j})$ and $\bar{\rho_{ij}}=0.5(\rho_{i}+\rho_{j})$. The artificial/Monaghan pressure [41] $P^{a}_{ij}$ as shown in Equation 11 is used to remove the tensile instability.

[TABLE]

where, $\gamma$ is a controlling parameter taken as $0.3$ for tension and $0.01$ for compression. The hydrostatic pressure is denoted by $P$. $\Delta x$ and $\Delta p$ are the current particle spacing and the average particle spacing in the reference configuration. $l$ is 4 [41].

The pressure $P$ is calculated from the linear equation of compressibility (Equation of State).

[TABLE]

where, $K$, $\nu$ and $\rho_{0}$ are the bulk modulus, Poisson’s ratio and initial mass density of the material respectively. The Cauchy stress components $\sigma^{\alpha\beta}$ are calculated from the deviatoric stress components $S^{\alpha\beta}$ and hydrostatic pressure $P$ as $\sigma^{\alpha\beta}=S^{\alpha\beta}-P\delta^{\alpha\beta}$. The rate of change in deviatoric stress components is calculated by Jaumann rate as,

[TABLE]

where, $\mu$ represents the shear modulus and the components of strain rate ($\dot{\bm{\epsilon}}$) and spin tensor ($\bm{\omega}$) are computed as,

[TABLE]

[TABLE]

Von-Mises yield function $y_{f}=\sqrt{J_{2}}-\frac{\sigma_{y}}{\sqrt{3}}$ is used to account for the effects of metal plasticity. $\sigma_{y}$ and $J_{2}=\frac{1}{2}S^{\alpha\beta}S^{\alpha\beta}$ are the yield stress of the material and the second invariant of deviatoric stress tensor respectively. The Wilkins criterion $S_{n}^{\alpha\beta}=c_{f}S^{\alpha\beta}$ is used for return mapping of the yield surfaces and $c_{f}=\min\left(\frac{\sigma_{y}}{\sqrt{3J_{2}}},1\right)$. The plastic strain increment, effective plastic strain increment and accumulation of plastic work density are estimated as

[TABLE]

[TABLE]

[TABLE]

The standard predictor-corrector method [42] is used to update the discretised Equations 2, 4, 6 and the stable time step for integration is determined from CFL (Courant-Fredrich-Levy) condition as, $\bigtriangleup t=\min_{i}\left(c_{s}\frac{h_{i}}{C_{i}+|v_{i}|}\right)$ where, CFL number $c_{s}$ is 0.3, $C_{i}$ is elastic wave velocity, $h_{i}$ is smoothing length and $v_{i}$ is velocity of the respective particles.

2.1 Pseudo-Spring analogy for crack

The pseudo-spring model for crack propagation is discussed in details in [30], [34], [35]. Here, a quick review of the algorithm is provided. In SPH, the particles interact with each other through the kernel function. However, the kernel function does not have any material property as shown in Equation 10. Hence, the initiation of crack and its growth cannot be modelled through it directly. A particle influences all the particles in its neighbourhood and the level of interaction between any particle pair decreases as their intermediate distance increases. As the material deforms, a particle may move out of another particle’s influence/neighbour domain, and the interaction between the particles stops. However, this cannot be considered as material damage/failure as different values of smoothing length would lead to different types of material response as discussed in [34]. Therefore, the ‘pseudo-spring’ analogy is used to incorporate the effects of material damage and fracture in the kernel function.

Here, a particle does not interact with all the particles found in its neighbourhood. The interaction of a particle is restricted to its immediate neighbourhood only as shown in Figure 1. All the particles are connected to its modified neighbourhood through a set of springs. These springs define the level of interaction between particles and most importantly these do not provide any extra stiffness to the system. Hence, the term ‘pseudo’ is used. At the start of the computation, the material is undamaged and the springs ensure complete interaction between connecting particles. However, as the material undergoes deformations, the damage accumulates in the springs. When the accumulated damage is more than a predefined value, it is assumed that the springs break. This leads to no interaction between connecting particles and propagation of crack path. This analogy does not require any particle splitting [43]. The concept is illustrated in Figure 1 and 2.

The analogy is incorporated in the framework in the following way. A interaction factor $f_{ij}$ between a connecting pair of particles is defined as $f_{ij}=1-D_{ij}$ and $D_{ij}$ is the accumulated damage in the spring between the $i-j$ particle pair. At the start of the computation, there is no material damage. Hence, the damage index $D_{ij}=0$ and this leads to the interaction factor $f_{ij}=1$ implying complete interaction between particle pairs. However, as the damage accumulates in the springs, the interaction factor $f_{ij}$ decreases i.e. $D_{ij}\in(0,1)\Rightarrow f_{ij}\in(0,1)$ or partial interaction. When the damage index $D_{ij}=1$ i.e. the interaction factor $f_{ij}=0$ or no interaction (leading to formation of crack surface). $D_{ij}$ is computed as $D_{ij}=0.5(D_{i}+D_{j})$ , where $D_{i}$ and $D_{j}$ are the damage parameters at the $i-th$ and $j-th$ particles respectively.

$\mathbb{N}^{i}\in\mathbb{Z^{+}}$, $\mathbb{N}_{D}^{i}=\{j\in\ {\mathbb{N}^{i}}\mid 0\leq f_{ij}<1\}\subseteq\mathbb{N}^{i}$ and $\mathbb{N}_{U}^{i}=\mathbb{N}^{i}\setminus\mathbb{N}_{D}^{i}$ are the set of immediate neighbouring particles that interact with the $i-th$ particle, the set of particles which are connected to the $i-th$ particle through damaged and undamaged springs at a given time instant respectively. The interaction between particles $i$ and $j\in\mathbb{N}^{i}$ in a damaged configuration is reduced by modifying the corresponding derivative of the kernel function as $\bar{W}_{ij,\beta}=f_{ij}\bar{W}_{ij,\beta}\forall j\in\mathbb{N}_{D}^{i}\subseteq\mathbb{N}^{i}$. Figure 3 shows the typical modified kernel in 2D for different damage states i.e.,- $f_{ij}=1$ (no damage), $0<f_{ij}<1$ (partial damage) and $f_{ij}=0$ (open crack).

Finally the modified forms of the discrete conservation equations 2, 4 and 6 become,

[TABLE]

[TABLE]

[TABLE]

2.2 Stress boundary condition

A stress boundary condition for SPH is developed in the context of solid continua in [44]. The boundary condition can be applied with or without ghost/dummy particles. The modified form of the momentum equation to account for the applied stress is shown in Equation 22.

[TABLE]

Here, $\sigma^{\alpha\beta}_{0}$ is the stress component applied to the boundary particles.

3 Numerical simulation

3.1 Two dimensional numerical examples

In this section, the crack initiation, propagation and branching in solids under dynamic loading are simulated for two-dimensional geometry. These phenomena are simulated using the Pseudo-Spring SPH. First, the crack branching of a pre-notched plate under dynamic axial loading is studied, and the effect of the loading amplitude is also investigated here. The computed crack speeds, their paths, and crack surfaces are compared with experimental and numerical results available in the literature. The crack path from a notch in a plate with a circular hole is also modelled.

3.1.1 Crack branching

A major advantage of the pseudo spring analogy compared to other crack tracking approaches is that it can also capture complex crack path evolutions without requiring any special treatment. In this section, the crack branching phenomenon in a brittle material is modelled. Towards this, a rectangular glass plate of 100 mm length and 40 mm width is considered. A pre-existing notch of 50 mm length is located symmetrically as shown in Figure 4. The plate is under tensile loading ($\sigma$ = 1 MPa) at the boundaries along the length of the plate. The material parameters used for the present computation are from [45, 46]. The material and SPH parameters used for the current simulations are given in Table 1. The damage parameter $D$ in the pseudo-springs is obtained as below.

[TABLE]

In the experiment [47], it is observed that initially the crack propagates in a self-propagating mode and after few microseconds, the branching takes place as shown in Figure 5. The predicted crack paths (crack initiation, propagation and branching) as obtained via the present simulations are shown in Figure 6. The crack propagation initiates around 10 $\mu$s and propagates in a straight line up to 27 $\mu$s. The initiation of branching of crack is observed at 28 $\mu$s. The crack propagates in two symmetrical branches and reaches the plate boundary at 62 $\mu$s.

The crack branching problem is studied in [9, 18, 24, 46, 49]. The final crack path obtained from the present simulation is compared with other numerical predictions available in the literature as shown in Figure 7. The crack propagation speed over time is compared with different techniques. The present results compare reasonably well with the results obtained in [46] using peridynamics, [45] using XFEM-CZM, and [50] using XFEM-CNM. The crack propagation speed increases after the crack initiation and decreases after the crack branching. All the predicted phenomena of the crack branching problem can be captured using the present pseudo-spring SPH framework, and the present predictions are comparable with the experimental and other numerical observations.

3.1.1.1 Effect of load intensity

Next, the effect of different tensile loading ($\sigma$) intensities on branching behaviour is explored. $\sigma$ values of 0.3, 1.1, 1.2, 2.0 and 4.0 MPa are considered for the simulations. All the other parameters are kept unchanged. It can be observed that the branching phenomenon depends on the loading ($\sigma$) intensity. For, $\sigma$ = 0.3 MPa loading, no branching is observed, and the crack propagates in a straight line as shown in Figure 9.

But the branching can be observed for all the other load intensities. The crack branching starts around 15 $\mu$s for $\sigma$ = 1.1 MPa (Figure 10), 14 $\mu$s for $\sigma$ = 1.2 MPa (Figure 11), 9 $\mu$s for $\sigma$ = 2.0 MPa and 6 $\mu$s for $\sigma$ = 4.0 MPa (Figure 12). As the loading amplitude ($\sigma$) is increased, the point where crack branching takes place, moves closer to the initial pre-notch position. The increase in loading amplitudes also results in secondary crack branches in the glass plates. No secondary crack branch is observed for $\sigma$ = 1.1 MPa but the secondary branches are found for $\sigma$ = 1.2, 2.0 and 4.0 MPa. The number of secondary branches also increases with the increase in loading amplitudes. These observations are consistent with [51].

Crack branching or splitting of cracks into multiple paths are one of the features of brittle fracture. The brittle materials undergoing fracture have a process zone [52] and the microcracks in this zone undergo nucleation, growth, and coalescence. These lead to crack path instabilities [53]. The crack branching occurs possibly due to the interaction between the reflected wave from the boundaries and the crack tip (waves generated due to crack propagation) or due to the pile-up of waves in the process zone [54]. This leads to the multiple branches of crack paths for the higher amplitude of boundary stress.

3.1.2 Pre-notched plate with a circular hole

In this section, the present framework is applied to model the crack path for a plate with a $10$ mm notch and an off-centre circular hole under tensile loading as shown in Figure 13. One end of the plate is clamped, and another end is under uniform tensile stress $\sigma$. This numerical study aims to simulate the crack propagation originating from the pre-existing notch and to track the changes in crack path due to the presence of the off-centre circular hole. The problem is studied in [55] by the arbitrary local mesh replacement method (a finite element based strategy). Tabiei and Wu [56] used DYNA3D and Kosteski et al. [57] used truss-like discrete element method to model the crack paths. Dipasquale et al. [58] used adaptive grid refinement in 2D peridynamics to track the crack paths. Here, the present framework is used to track the crack paths. The material and SPH parameters for the simulation are given in Table 2. The damage parameter $D$ is obtained at any particle from Equation 23.

Three different notch locations from the fixed support, named as case I ($h=5$ mm), case II ($h=10$ mm) and case III ($h=15$ mm) [59] are considered. Figure 14 shows the crack evolution for case I. The crack path does not interact with the circular hole, and the crack propagates from the notch to the other end without going near the circular hole.

The crack propagation for case II (notch distance 10 mm) and case III (notch distance 15 mm) are given in Figure 15 and 16 respectively. The crack path deviation from the off-centre circular hole depends on the notch location. As the distance from the support increases, the deviation of the crack path from the circular hole decreases. This deviation is largest for the case I (notch distance 5 mm). The crack path interacts with the circular hole for case III. The predicted crack paths for all the three cases are consistent with the observations in [57, 59] as shown in the Figure 14, 15 and 16.

3.2 Three dimensional numerical examples

Numerical modelling of arbitrary propagation and interaction of cracks in three dimension (3-D) is a challenging exercise. Crack propagation in 3-D is much more difficult to capture than the same in 2-D, and it poses a real test for a computational framework to demonstrate its efficacy. This motivated the extension of the present computational framework in 3-D. Herein, the pseudo spring SPH is explored to model arbitrary crack propagation in the 3-D solid continuum. Towards this end, a few example problems are considered. The first example is the formation of a helical crack in cylindrical chalk subjected to torsional loading at both ends. The Kalthoff-Winkler experiment and the failure of Taylor bullet are also modelled. The predicted fracture surfaces are found to be in good agreement with the experimental results available in the literature.

3.2.1 Cylindrical chalk under torsion

In this section, the failure of a cylindrical chalk bar under torsional loading is modelled. The diagram of the set-up and loading on the chalk specimen are shown in Figure 17. The specimen has a length of 100 mm and a diameter of 8 mm. The torsional loading is generated by applying a velocity field of $u_{max}$ = 28.3 m/s in the tangential direction to the chalk bar as shown in Figure 17. The distribution of the velocity field is given as

[TABLE]

The velocity field is applied for a length of 20 mm each at both ends of the specimen. The application of the velocity field is such that it generates a twisting behaviour in the chalk. This has the similar effect of the traction field applied in the tangential direction of the chalk creating torsional loading [60, 61]. The material properties and other computational parameters are given in Table 3.

The maximum principal stress criterion is used for modelling of the failure. The springs are considered to fail completely when the maximum principal stress exceeds $\sigma^{p}_{max}$ value. The failure can happen anywhere and results in an arbitrary crack surface as shown in Figure 18. The present fracture surfaces are found to be helical. These surfaces are consistent with the crack surface obtained in [60, 61] for the similar problem using the extended element-free Galerkin (XEFG) method.

3.2.2 Kalthoff-Winkler experiment

As the second example, the Kalthoff-Winkler experiment [62] is considered. A cylindrical projectile impacts a plate with two edge notches as shown in Figure 19. Numerical simulations of similar problems are performed in [48, 26, 49, 57, 51, 30]. At a low-velocity impact, the plate undergoes brittle failure, and crack propagates at an angle 70 ∘ to the notch direction. The material and SPH parameters are given in Table 4. The impact velocity is taken as 16.5 m/s.

In the present numerical simulation, half of the projectile-target system is considered due to the symmetry and accordingly a symmetric boundary condition is applied. The Rankine’s criteria or critical stretch between particles is used as the failure [30, 49] criterion. The critical stretch at any time instant is calculated as

[TABLE]

When the calculated stretch ($\delta_{t}$) at any time instant exceeds the critical stretch ($\delta_{tc}$) value, the material is assumed to be fully damaged. The crack pattern using the current framework is shown in Figure 20 and 21. The crack path from the present simulation is found to be approximately 75∘ which is close to the observed crack path through experiment. A second crack path is found in the numerical results. This is not found in the experimental observation. However, a similar type of second crack path was reported in [18, 51, 57, 58]. This crack path originates possibly due to the reflection of the wave from the boundary or due to improper estimation of material parameters.

3.2.3 Taylor bullet impact

A flat-ended cylindrical projectile of length 30 mm and diameter 6 mm is moving towards a rigid wall as shown in Figure 22. The failure of the projectile is modelled in this section. The bullet is made of Weldox 460E steel. The Johnson-Cook material and damage models are used. The material and damage parameters may be found in [35]. The SPH parameters are given in Table 5. Initial impact velocity of the projectile is 600 m/s, and initial particle distribution is shown in Figure 23. The petalling fracture mode for the projectile is simulated in [63]. The present framework also predicts petalling fracture mode as shown in Figure 24.

4 Closure

Numerical modelling of crack initiation, propagation and branching is always a challenging task. Towards this, several attempts are made, but not all are efficient. The element-based methods can simulate the arbitrary cracks but have their own set of disadvantages namely intensive computation for local enrichment, re-meshing, element deletion etc. The mesh-free methods also require enrichments or particle splitting.

In this work, the SPH along with the ‘pseudo-spring’ analogy is used to model arbitrary crack propagation. The ‘pseudo-spring’ analogy provides a way to model discontinuity in solid continua without any enrichment, particle splitting or crack path tracking. The particles are connected by springs which do not provide any extra stiffness to the system but reduce the level of interaction between particles based on the damage state of the springs. The interaction of particles is only with the immediate neighbours. If a spring fails completely, the interaction between connecting particles ceases. This implies the passing of crack paths through the spring and generation of a new surface. This framework also preserves the simplicity of SPH.

The model is used to study two and three-dimensional problems. First, A pre-notched plate under tensile loading is studied. Initially, the crack propagates in a straight line and then splits into two branches. Both the crack branches propagate until the end of the plate. The crack speed and its path are in good agreement with the experimental and other numerical results from the literature. It is observed that the crack propagates in a straight line without any branches for low-intensity load. However, the crack splits into branches for higher load intensity. If loading intensity is further increased, multiple branches and sub-branches can be observed. Next, a notched plate with a circular hole under tensile load is considered. The position of the notch is changed and the crack path propagation is studied. The propagation of crack path and their interaction with the circular hole are found to be in good agreement with the results from the literature. The computational framework is also used to model arbitrary crack paths and surfaces in three dimensions. The fracture of a chalk bar under torsion is considered. The numerically obtained crack surface is similar to the experimental and numerical results. The Kalthoff-Winkler experiment and the Taylor impact test are also simulated for brittle fracture and the crack paths are found to be comparable to the experimental and numerical results. The ‘pseudo-spring’ SPH is found to be very efficient for the present simulations of crack initiation, propagation and branching without any special treatment. In three dimensional problem, the arbitrary crack surface can also be traced without much difficulty. The methods seem to be very good in tacking discontinuity and easy to implement. It is anticipated that the pseudo-spring analogy has the potential to constitute an efficient computational tool for modelling more complex fracture problems in both 2D and 3D continua.

Acknowledgement

The authors would like to acknowledge the Defence Research and Development Organization (DRDO), India for their support of this work.