A damage and failure implementation for the simulation of ductile solids with Total-Lagrangian Smooth Particle Hydrodynamics

TL;DR

This paper introduces a new damage and failure model for Total-Lagrangian SPH to better simulate ductile solids under large deformations, incorporating damage criteria and validated against experiments and FEM results.

Contribution

It develops a novel damage and failure model for TLSPH that accounts for material property changes and discontinuities, improving ductile solid simulations.

Findings

The new model accurately predicts damage initiation and growth in ductile materials.

Validation shows improved agreement with experimental data and FEM simulations.

The implementation enhances TLSPH's stability and predictive capabilities for ductile failure.

Abstract

Smooth-Particle-Hydrodynamics is gaining popularity for the simulation of solids subjected to machining, wear, and impacts. Its attractiveness is due to its abilities to simulate problems involving large deformations resulting from the absence of mesh, recent improvements in stability conferred by the development of Total-Lagrangian version of SPH (TLSPH), but also its availability in the open-source software LAMMPS. This implementation features a damage model similar to the `pseudo-spring' method which creates instabilities when used for the simulation of ductile materials. In this contribution, we present a new damage and failure model for TLSPH suitable for ductile materials. In this implementation, not only the constitutive equations but also the TLSPH approximation are modified in order to take into account the change in material properties as well as the presence of discontinuities due to the initiation and growth of damage. This new approach is accompanied by the implementation of the Cockroft-Latham, Johnson-Cook, and Gurson-Tvergaard-Needleman damage criteria. The predictive capabilities of the implementation of this new damage model are then tested and compared against both experimental and results of Finite Element simulations.

Tables

Table 1. Table 1: Material parameters for Weldox steels. ρ 0 subscript 𝜌 0 \rho_{0} is the reference bulk density, E 𝐸 E the Young modulus, ν 𝜈 \nu the Poisson’s ratio, c 0 = E / ( 2 ​ ( 1 − 2 ​ ν ) ​ ρ 0 ) subscript 𝑐 0 𝐸 2 1 2 𝜈 subscript 𝜌 0 c_{0}=\sqrt{E/(2(1-2\nu)\rho_{0})} is the bulk speed of sound, Γ 0 subscript Γ 0 \Gamma_{0} the Grüneisen Gamma in the reference state.

${\rho }_0$ [m/s]	$E$ (GPa)	$\nu$	$c_0$ [m/s]	$S_{\alpha }$	${\Gamma }_0$
7750	211	0.33	5166	1.5	0

Table 2. Table 2: Material constants for the Johnson-Cook constitutive model. [ 36 ]

Material	Yield stress	Strain hardening
	A (MPa)	B (MPa)	$n$
Weldox 460E	499	382	0.458
Weldox 700E	859	329	0.579
Weldox 900E	992	364	0.568

Table 3. Table 3: Material constants for three damage and failure models used here. The parameters for the Cockroft-Latham and Johnson-Cook failure criteria were proposed by Dey et al.   [ 36 ] , while that for the GTN model were proposed by Øien and Schønberg  [ 42 ] .

Material	Cockroft-Latham	Johnson-Cook					GTN
	$W_{cr}$ (MPa)	$D_1$	$D_2$	$D_3$	$D_4$	$D_5$	$f_0$	$f_{cr}$	$f_F$	$k_{\omega }$
Weldox 460E	1219	0.636	1.936	-2.969	-0.0140	1.014	0.0001	0.01	0.2	5.5
Weldox 700E	1424	0.361	4.768	-5.107	-0.0013	1.333	0.005	0.25	0.3	3.5
Weldox 900E	1510	0.294	5.149	-5.583	0.0023	0.951	0.005	0.18	0.2	3.5

Equations

$$\mathbf{u}=\mathbf{x}-\mathbf{X},$$

$$\rho J$$

$$\ddot{\mathbf{u}}$$

$$\dot{e}$$

$$\mathbf{F}=\frac{d\mathbf{x}}{d\mathbf{X}}=\frac{d\mathbf{u}}{d\mathbf{X}}+\mathbf{I},$$

$$⟨\mathrm{f}(\mathbf{x})⟩={\int }_{\Omega }\mathrm{f}({\mathbf{x}}^{\prime })W({\mathbf{x}}^{\prime }-\mathbf{x},h)d{\mathbf{x}}^{\prime }$$

$${\int }_{\Omega }W({\mathbf{x}}^{\prime }-\mathbf{x},h)d{\mathbf{x}}^{\prime }=1$$

$${\int }_{\Omega }({\mathbf{x}}^{\prime }-\mathbf{x})W({\mathbf{x}}^{\prime }-\mathbf{x},h)d{\mathbf{x}}^{\prime }=0$$

$$⟨\nabla \mathrm{f}(\mathbf{x})⟩={\int }_{\Omega }\mathrm{f}({\mathbf{x}}^{\prime })\nabla W({\mathbf{x}}^{\prime }-\mathbf{x},h)d{\mathbf{x}}^{\prime }.$$

$$⟨{\nabla }_0\mathrm{f}(\mathbf{X})⟩={\int }_{\Omega }\mathrm{f}({\mathbf{X}}^{\prime })\nabla W({\mathbf{X}}^{\prime }-\mathbf{X},h)d{\mathbf{X}}^{\prime }.$$

$${\nabla }_0\mathrm{f}({\mathbf{X}}_i)=\sum _{j\in S_i}{V}_j^0\mathrm{f}({\mathbf{X}}_j)\nabla W_i(X_{ij}),$$

$$\nabla W_i(X_{ij})=\frac{d{W}_i(X_{ij})}{d{X}_{ij}}\frac{{\mathbf{X}}_j-{\mathbf{X}}_i}{X_{ij}}.$$

$$\sum _{j\in S_i}{V}_j^0{W}_i(X_{ij})$$

$$\sum _{j\in S_i}{V}_j^0\nabla W_i(X_{ij})$$

$${\nabla }_0\mathrm{f}({\mathbf{X}}_i)=\sum _{j\in S_i}{V}_j^0\left(\mathrm{f}({\mathbf{X}}_j)-\mathrm{f}({\mathbf{X}}_i)\right)\nabla W_i(X_{ij}).$$

$${\nabla }_0\mathrm{f}({\mathbf{X}}_i)=\sum _{j\in S_i}{V}_j^0\left(\mathrm{f}({\mathbf{X}}_j)-\mathrm{f}({\mathbf{X}}_i)\right)\nabla W_i(X_{ij}){\mathbf{L}}_i^{-1},$$

$${\mathbf{L}}_i=\sum _{j\in S_i^i}{V}_j^0({\mathbf{X}}_j-{\mathbf{X}}_i)\otimes \nabla W_i(X_{ij})$$

$${\int }_a^b\mathrm{f}(x)W_i(x)dx={\int }_a^d\mathrm{f}(x)W_i(x)dx+{\int }_d^b\mathrm{f}(x)W_i(x)dx.$$

$$\begin{array}{rl}{\int }_a^b\mathrm{f}(x)W_i(x)dx& =\mathrm{f}(x_i){\int }_a^d{W}_i(x)dx+\mathrm{f}(x_k){\int }_d^b{W}_i(x)dx\\ & +{\mathrm{f}}^{\prime }(x_i){\int }_a^d(x-x_i)W_i(x)dx\\ & +{\mathrm{f}}^{\prime }(x_k){\int }_d^b(x-x_k)W_i(x)dx+o(h^2).\end{array}$$

$$\begin{array}{rl}{\int }_a^b\mathrm{f}(x)W_i(x)dx& =\mathrm{f}(x_i){\int }_a^b{W}_i(x)dx+(\mathrm{f}(x_k)-\mathrm{f}(x_i)){\int }_d^b{W}_i(x)dx\\ & +{\mathrm{f}}^{\prime }(x_i){\int }_a^b(x-x_i)W_i(x)dx\\ & +{\int }_d^b[(x-x_k){\mathrm{f}}^{\prime }(x_k)-(x-x_i){\mathrm{f}}^{\prime }(x_i)]W_i(x)dx+o(h^2),\end{array}$$

$$\begin{array}{rl}{\int }_a^b\mathrm{f}(x)& \nabla W_i(x)dx=\mathrm{f}(x_i){\int }_a^b\nabla W_i(x)dx+(\mathrm{f}(x_k)-\mathrm{f}(x_i)){\int }_d^b\nabla W_i(x)dx\\ & +{\mathrm{f}}^{\prime }(x_i){\int }_a^b(x-x_i)\nabla W_i(x)dx\\ & +{\int }_d^b[(x-x_k){\mathrm{f}}^{\prime }(x_k)-(x-x_i){\mathrm{f}}^{\prime }(x_i)]\nabla W_i(x)dx+o(h^2),\end{array}$$

$$\begin{array}{rl}\sum _{j=1}^N{V}_j^0{\mathrm{f}}_j\nabla W_i(X_{ij})& =f_i\sum _{j=1}^N{V}_j^0\nabla W_i(X_{ij})+(f_k-f_i)\sum _{j=k}^N{V}_j^0\nabla W_i(X_{ij})\\ & +{\mathrm{f}}_i^{\prime }\sum _{j=1}^N{V}_j^0(X_j-X_i)\nabla W_i(X_{ij})\\ & +\sum _{j=k}^N{V}_j^0[(X_j-X_k){\mathrm{f}}_k^{\prime }-(X_j-X_i){\mathrm{f}}_i^{\prime }]\nabla W_i(X_{ij}),\end{array}$$

$$\begin{array}{rl}\sum _{j=1}^N{V}_j^0{\mathrm{f}}_j& \nabla W_i(X_{ij})={\mathrm{f}}_i\sum _{j=1}^N{V}_j^0\nabla W_i(X_{ij})+({\mathrm{f}}_k-{\mathrm{f}}_i)\sum _{j=k}^N{V}_j^0\nabla W_i(X_{ij})\\ & +\nabla {\mathrm{f}}_i\sum _{j=1}^N{V}_j^0({\mathbf{X}}_j-{\mathbf{X}}_i)\otimes \nabla W_i(X_{ij})\\ & +\nabla {\mathrm{f}}_k\sum _{j=k}^N{V}_j^0({\mathbf{X}}_j-{\mathbf{X}}_k)\otimes \nabla W_i(X_{ij})-\nabla {\mathrm{f}}_i\sum _{j=k}^N{V}_j^0({\mathbf{X}}_j-{\mathbf{X}}_i)\otimes \nabla W_i(X_{ij}).\end{array}$$

$$\begin{array}{rl}\nabla {\mathrm{f}}_i& =\sum _{j\in S_i}{V}_j^0({\mathrm{f}}_j-{\mathrm{f}}_i)\nabla W_i(X_{ij}){\mathbf{L}}_i^{-1}-({\mathrm{f}}_k-{\mathrm{f}}_i)\sum _{j\in S_i^k}{V}_j^0\nabla W_i(X_{ij}){\mathbf{L}}_i^{-1}\\ & -\nabla {\mathrm{f}}_k\sum _{j\in S_i^k}{V}_j^0({\mathbf{X}}_j-{\mathbf{X}}_k)\otimes \nabla W_i(X_{ij}){\mathbf{L}}_i^{-1},\end{array}$$

$$\nabla {\mathrm{f}}_i=\sum _{j\in S_i^i}{V}_j^0({\mathrm{f}}_j-{\mathrm{f}}_i)\nabla W_i(X_{ij}){\mathbf{L}}_i^{-1}.$$

$${\mathbf{F}}_i$$

$${\dot{\mathbf{F}}}_i$$

$$\dot{e}=\sum _{i\in \Omega }\frac{1}{m_i}{\mathbf{f}}_i{\mathbf{v}}_i,$$

$$\dot{e}=\sum _{i\in \Omega }\sum _{j\in S_i^i}\frac{1}{m_i}{V}_i^0{V}_j^0({\mathbf{v}}_j-{\mathbf{v}}_i)\otimes \nabla W_i(X_{ij}){\mathbf{L}}_i^{-1}:{\mathbf{P}}_i.$$

$$\dot{e}=\sum _{i\in \Omega }\sum _{j\in S_i^i}\frac{1}{m_i}{V}_i^0{V}_j^0\left[{\mathbf{P}}_i{\mathbf{L}}_i^{-T}+{\mathbf{P}}_j{\mathbf{L}}_j^{-T}\right]:{\mathbf{v}}_i\otimes \nabla W_i(X_{ij})$$

$${\mathbf{f}}_i=\sum _{j\in S_i^i}{V}_i^0{V}_j^0\left[{\mathbf{P}}_i{\mathbf{L}}_i^{-T}+{\mathbf{P}}_j{\mathbf{L}}_j^{-T}\right]\nabla W_i(X_{ij}).$$

$${\mathbf{f}}_{ij}=V_i^0{V}_j^0{\mathbf{P}}_i{\mathbf{L}}_i^{-T}\nabla W_i(X_{ij})\ne 0,$$

Full text

††journal: Computer Methods in Applied Mechanics and Engineering

1 Introduction

Smooth Particle Hydrodynamics (SPH) is a meshfree particle method originally developed in the seventies for modeling three-dimensional problems in astrophysics [1, 2]. Recent extensions now allow SPH to be used to simulate solid mechanics problems [3, 4]. Its principle is the use of a discrete kernel function centered about particles to determine, at their location point, an approximation of continuous fields such as density, velocity, stresses, and internal energy. This principle, when applied to solid body deformations transforms the underlying set of partial differential equations into simple algebraic equations.

Due to its abilities to simulate problems involving large deformations, SPH is gaining popularity for the study of the deformations of solids in the fields of machining [3, 4], wear [5, 6, 7], and resistance to impacts [8, 9, 10, 11]. Such popularity has also been fueled by its recent implementation in the open-source particle code LAMMPS [12]. This implementation is shipped in the user package called Smooth Mach Dynamics (SMD) created by Ganzenmüller [13] and is based on the numerically stable and efficient Total-Lagrangian version of SPH (TLSPH).

In cases involving large deformations and damage, true meshfree particle methods such as SPH have several advantages over mesh-based methods such as Finite Element Method (FEM). One advantage is the absence of need to mesh, or remesh when large distortions occur. Another advantage is its the ability to realistically simulate damage and fracture due to the absence of immutable connectivity between neighbouring particles [14]. In fact, the decrease or loss of interaction between two neighbouring particles could be seen respectively as the presence of damage or as the formation of a crack between them. This concept is exploited in the so called “pseudo-spring” damage and fracture method [11]. In this method, the kernel function in the direction of every virtual segment linking two neighbouring particles, called a “pseudo-spring”, is scaled according to the amount of damage present in that “pseudo-spring”. When the amount of damage reaches the critical limit, the kernel function in the “pseudo-spring’s” direction becomes null, and the connectivity between particles is lost. This method is similar to the use of cohesive elements [15] in mesh-based methods which are particularly adapted to the simulation of debonding, and failure of brittle solids.

However, the “pseudo-spring” fracture model has three limitations.

It imposes the calculation of a quantity in a non-integration point (i.e. between particles). 2. 2.

It breaks the mathematical approximation upon which SPH is built by scaling the kernel functions. 3. 3.

It imposes particles to be linked to their immediate neighbours only, as done by Chakraborty and Shaw [11], or the interaction with non-intermediate neighbours separated from that particle by a crack to be scaled appropriately [16].

The other currently available damage and fracture implementation for SPH is the particle cracking method [17, 18, 19]. Contrary to the “pseudo-spring” method, damage is calculated at the center of particles (i.e. the integration points). When damage reaches a critical limit, the particle is split into two. In this case, damage is not localized as its value corresponds to an average over each particle’s corresponding volume. Therefore, the precise location of cracks is not defined and the way particles need to be split is non-trivial.

Currently, the only damage models implemented in TLSPH LAMMPS package for solids were models for brittle failure using a version of the “pseudo-spring” failure method [13]. No damage law for ductile material has been fully implemented. Varga, Leroch and coworkers [5, 20, 6] have used the TLSPH LAMMPS package to simulate the mechanical response of a ductile material to scratch, without taking into account damage. By recalculating the particles’ neighbour lists after a few timesteps, they have relied on numerical failure, more specifically the loss of connectivity between neighbouring particles as their distance increase, to account for damage. This so-called “natural ability” for SPH to account for damage has also been used by other researchers [3, 7, 9, 10, 11]. However, this is unsatisfactory when it comes to predict the mechanical response of ductile materials experiencing damage.

No damage law for ductile material is currently available in LAMMPS. Only an undocumented pre-implementation of Johnson-Cook damage law is present but is flagged as currently unsupported. Moreover, when used with the “pseudo-spring” failure method, this implementation suffers instability issues triggered by the start of the propagation of failure.

In this paper, a new implementation for damage and failure for the LAMMPS TLSPH package is introduced. In this approach, no alteration of the kernel functions is used, and the particles’ neighbour list remains identical throughout each simulation. Thus, this approach does not rely on numerical fracture to account for damage. Instead, the inherent displacement discontinuity generated by the presence of cracks or voids is taken into account by scaling the velocity difference between neighbouring damaged particles. Damage is computed at the center of particles where cracks are assumed to be located. This new method is accompanied by the implementation of three damage laws: Cockroft-Latham [21], Johnson-Cook [22], and Gurson-Tvergaard-Needleman [23]. The predictive capabilities of the new implementation are then tested and compared against both experimental and Finite Element simulations.

2 Total-Lagrangian SPH theory for solids

Smooth Particle Hydrodynamics is a particle based method in which a continuous region of solid is discretized into point particles to which a certain volume is attached. Each particle $i$ interacts with its neighbours ($S_{i}$) by means of a smoothing kernel function $W_{i}$ centered around each particle (Fig. 1).

In the original formulation, the kernel function travels with the particle, making it Eulerian. Unfortunately, this formulation suffers from a series of instabilities and low accuracy issues when applied to solid mechanics. The first issue SPH suffered from was the so-called “tensile instability” [24] which leads to particle clunging and numerical fracture. This was found to be caused by the use of Eulerian kernel functions [25] and was addressed by the development of Total-Lagrangian SPH (TLSPH) [26] in which kernel functions are fixed to the particles in their reference configuration, making the kernel functions Lagrangian. However, TLSPH was impaired by a rank-deficiency problem, which necessitates artificial particle velocity damping, and also by yet another tensile instability arising when the difference between the current and the reference configurations gets too large [27, 28, 29]. Rank-deficiency is solved via the use of artificial particle velocity damping, while the new tensile instability was addressed by the development of the hourglass control scheme similar to that used in Finite Elements Methods (FEM) [30].

Another problem crippling the use of Eulerian kernel functions for the simulation of solids is the risk of numerical fracture under large deformations. This occurs when the distance between two neighbouring particles in the undeformed configuration becomes larger than the kernel cutoff radius and these particles are no longer neighbours in the deformed state. If total damage has not been reached, this would lead to numerical failure. The use of Lagrangian kernel functions also solves this problem. In fact, in TLSPH, the SPH approximation is always determined in the reference (undeformed) configuration. Thus, particles which are neighbours in the underformed configurations will always remain as such.

All these improvements featured in the previous implementation of TLSPH for solids in LAMMPS are detailed in the work by Leroch et al. [5].

2.1 Governing equations

The solid mechanics problem to solve here is described by a set of conservation equations in Lagrangian. This methodology uses a reference configuration of the simulated domain for the computation of stresses and accelerations. Alongside the unchangeable reference configuration, the current set (i.e. deformed domain) is adjusted according to the accelerations computed in the reference configuration. Therefore, all the constitutive and conservation equations are expressed in terms of the reference coordinates $\bm{X}$, which in these simulations are taken as the initial configuration, which is supposed undeformed.

The displacement $\bm{u}$ being given by:

[TABLE]

the respective conservation equations for mass, impulse, and energy in the total-Lagrangian framework are [30]:

[TABLE]

where $\rho$ is the mass density, $\bm{P}$ the first Piola-Kirchhoff stress tensor, $e$ the internal energy per unit mass, and $\bm{\nabla}$ the gradient or divergence operator. The subscript [math] designate that the quantity or operator is evaluated in the reference configuration, and when no subscript is indicated the evaluation is in the current configuration. Finally, $J$ designates the determinant of the deformation gradient $\bm{F}$ given by:

[TABLE]

which can be interpreted as the transformation matrix that describes the rotation and stretch of a line element from the reference configuration to the current (deformed) configuration.

2.2 SPH approximation in the bulk in the absence of damage

The basis of the original SPH formulation is that the value of any continuous field function $\mathrm{f}$ at a point $\bm{x}$ can be approximated as $\langle\mathrm{f}(x)\rangle$ by the following integral interpolant [31]:

[TABLE]

where $h$ is the cutoff radius, and the function $W$ is the kernel which is even, normalized, and has a compact support (i.e. $W(\bm{x^{\prime}}-\bm{x},h)=0$ for $||\bm{x^{\prime}}-\bm{x}||\geq h$). Therefore, $W$ fulfills the following conditions:

[TABLE]

[TABLE]

From these equations, it can be seen that the Dirac function is a valid kernel function. If used, the integral interpolant reproduces exactly the function $\mathrm{f}$. It can also be seen that if the sphere of radius $h$ centered on $\bm{x}$ is astride a boundary, the kernel function being truncated by the boundary, the conditions given by Equations 7 and 8 are no longer satisfied [32].

Since the kernel has a compact support, in the absence of damage, the approximation of the spatial gradient of $f$ is [33]:

[TABLE]

This can be expressed in the reference configuration and yield the Lagrangian approximation of the spatial gradient:

[TABLE]

In the second step of the SPH approximation, the continuous simulated domain $\Omega$ is discretized into particles. Integrals are replaced by their discretized form written as a summation over the neighboring particles. The particle approximation of Equation 10 evaluated at the particle $i$ yields:

[TABLE]

where $S_{i}$ is the set of neighbouring particles of $i$, $X_{ij}=||\bm{X}_{j}-\bm{X}_{i}||$, and the definition of the gradient of the kernel function is:

[TABLE]

This formulation, however, does not fulfill zeroth- and first-order completeness conditions, i.e. the approximation of zeroth- and first-order polynomials are not exactly approximated. This is due to the non respect, in general, of the following conditions:

[TABLE]

In order to restore zeroth-order completeness, Monaghan [34] introduced the following ad-hoc improvement by adding Eq. 14 to Eq. 11:

[TABLE]

On the other hand, the correction for first-order completeness can be done using a corrected kernel gradient as detailed in Bonet and Lok [35]:

[TABLE]

where $\bm{L}_{i}$ is a correction matrix given by:

[TABLE]

3 SPH approximation at the boundaries

The original SPH formulation has been derived in the bulk for continuous field functions. At the boundaries, however, the field functions are discontinuous. Thus, there, the TLSPH approximation given by Eq. 16 might not be valid. Derivation of an SPH approximation at the boundaries can be done following the work of Liu et al. [33] who have introduced a new formulation which applies for problems with or without discontinuities called DSPH (for Discontinuous SPH).

The original 1D derivation of the DSPH approximation by Liu et al. [33] is as follows. Assuming that the kernel support for the particle $i$ is bounded by $a$, and $b$ , and that the function $f$ had an integrable discontinuity in at $d$, as shown in Fig. 2, the integral interpolant of $f$ can be divided into two parts:

[TABLE]

Performing a second order Taylor series expansion of $f$ around $x_{i}$, and an arbitrary point $x_{k}$ located on the opposite side of the discontinuity from $x_{i}$ (i.e. $d\leq x_{k}\leq b$) gives:

[TABLE]

Rearranging by combining similar terms yields:

[TABLE]

from which Liu et al. [33] extracted $\mathrm{f}^{\prime}(x_{i})$ and obtained its DSPH approximation by replacing the integrals by their particle approximations.

Here, we apply the same idea but to the integral interpolant of the gradient of $\mathrm{f}$ which yields:

[TABLE]

and its Lagrangian particle approximation version:

[TABLE]

which interpolated into 3D space gives:

[TABLE]

From this equation, the kernel approximation of the gradient of $\mathrm{f}$ in $i$ is extracted and yields

[TABLE]

where $\bm{L}_{i}=\sum_{j\in S_{i}^{i}}V_{j}^{0}(\bm{X}_{j}-\bm{X}_{i})\otimes\bm{\nabla}W_{i}(X_{ij})$, $S_{i}^{i}$ and $S_{i}^{k}$ are respectively the part of the support of $W_{i}$ located on the same side as and the opposite side of the discontinuity from the particle $i$.

Following the same methodology, but by splitting the interval $[d,b]$ into as many segments as there are particles, one can show that the formulation for the Lagrangian DSPH approximation of the kernel gradient is:

[TABLE]

This equation also applies in the case where no discontinuities are present and also in the case of a boundary. When no discontinuity is present $S_{i}^{i}$ is identical to $S_{i}$, and this formulation is identical to the Total-Lagrangian SPH approximation in Eq. 16. When a boundary is present, such as due to a crack or a void, $S_{i}^{i}$ is the portion of $S_{i}$ included within the domain $\Omega$ (Fig. 3).

The deformation gradient and its time derivative are obtained by calculating the derivative of the displacement and velocity field, respectively (see Equation 5). Using the Lagrangian DSPH approximation (Equation 25), they are calculated as:

[TABLE]

where $\otimes$ is the dyadic product of two vectors, and $\bm{v}$ is the particle velocity.

The nodal forces are derived from the conservation of energy following Bonnet and Lok [35]. Noting that the variation of the internal energy is

[TABLE]

and recalling Equations 4 and 33 for the energy conservation and $\dot{\bm{F}}$, respectively, the following equality is obtained:

[TABLE]

Rearranging the summations involved and noting the anti-symmetry property of the kernel function $\bm{\nabla}W_{i}(X_{ij})=-\bm{\nabla}W_{j}(X_{ji})$:

[TABLE]

Comparing this expression with Equation 28 gives the internal forces applied onto the particle $i$ as

[TABLE]

Owing to the anti-symmetry property of the kernel function mentioned above, this force expression will conserve linear momentum exactly, as $\bm{f}_{i}=-\bm{f}_{j}$.

3.1 TLSPH approximation with damage

Implementing damage into TLSPH begins by appropriately scaling the constitutive laws similarly to what is done in FEM. If the TLSPH approximation is kept unchanged, once a particle is fully damaged, it will continue applying a force on its neighbours. In fact, if $i$ is a neighbour of a fully damaged particle $j$, the force applied onto $i$ by $j$ would be:

[TABLE]

similarly, the deformation matrix and its time derivative at $i$ would still be dependant of $j$.

In order to solve this problem, two methods have been used. First, the “pseudo-spring” method consists of the scaling of the kernel function by a function of $D_{i}$ and $D_{j}$, damages at the particle $i$ and $j$, respectively. This function is chosen such that it is equal to $1$ when both $D_{i}=1$ and $D_{j}=1$, and [math] when either $D_{i}=1$ or $D_{j}=0$. In pure elasticity, this does not cause any problem as it is physically equivalent to a total unloading of the virtual $i$-$j$ “link”, i.e. the contribution of $j$ to the deformation and velocity gradient matrices at $i$ are respectively $\bm{F}_{ij}=\bm{I}$ and $\dot{\bm{F}}_{ij}=\bm{0}$. However, when plastic deformation has occurred, permanent deformation of this virtual link remains and $\bm{F}_{ij}$ cannot be equal to $\bm{I}$. Therefore, this explains why using the “pseudo-spring” method generates instability when used for the simulation of ductile materials. The second method consists of splitting any fully damaged particles. This raises the difficult question as to in how many particles should the broken ones be spit into. It would also generate implementation problems within LAMMPS. This is why it was necessary to develop a new approach.

Here the chosen approach consist of scaling the difference of velocity between $i$ and $j$, i.e. $v_{j}-v_{i}$, by a factor $1-D_{j}$. This is equivalent to having a crack located at the center of the particle $j$, and having the part of the volume corresponding to $j$ adjacent to $i$ moving at the same velocity as $i$. This is indeed the case in the limit of particles with infinitely low mass.

The amount of damage $D_{j}$ represents here the average damage in the volume $V_{j}$ attached to the particle $j$. $D_{j}$ is evaluated at its center point and varies between 0 (undamaged) and 1 (fully damaged). When a particle is fully damaged, it has lost all stress carrying capacities but is not deleted. For practical reasons only, its velocity is set to [math]. When $j$ is undamaged, its contribution to the deformation matrix and its time derivative are respectively $\bm{F}_{ij}=V_{j}^{0}(\bm{u}_{j}-\bm{u}_{i})\bm{\nabla}W_{i}(X_{ij})\bm{L}_{i}^{-1}+\bm{I}$ and $\dot{\bm{F}}_{ij}=V_{j}^{0}(\bm{v}_{j}-\bm{v}_{i})\bm{\nabla}W_{i}(X_{ij})\bm{L}_{i}^{-1}$. When $j$ is fully damaged, they become respectively $\bm{F}_{ij}=\textrm{const.}$ and $\dot{\bm{F}}_{ij}=0$. Assuming that $\dot{\bm{F}}_{ij}$ is linearly depending on $D_{j}$ yields:

[TABLE]

and after time integration:

[TABLE]

Using Equation 33 and recalling the energy conservation equations, the nodal forces when particles are damaged are expressed as

[TABLE]

One can see that if $D_{j}=1$, $\bm{f}_{ij}=0$ since $\bm{P}_{j}=\bm{0}$, and that similarly to Equation 31, this expression conserves linear momentum exactly.

4 Constitutive model

In order to close the set of Equations 2 to 4, we modified the constitutive relationships used by Leroch et al. [5] to take into account the effects of damage. The symmetric stress tensor is expressed as the sum of its isotropic part, i.e. the hydrostatic pressure ($\sigma_{m}$), and the traceless symmetric deviatoric stress $\bm{\sigma_{d}}$:

[TABLE]

The hydrostatic pressure, on the one hand, is estimated using an equation of state (EOS). Multiple EOS are implemented in the SMD package (see the SMD documentation for more details [13]), but as Leroch et al. [5], we use here the Mie-Grüneisen equation:

[TABLE]

where $c_{0}$ is the bulk speed of sound, $\Gamma_{0}$ the Grüneisen Gamma in the reference state.

The deviatoric stress, on the other hand is determined using empirical models for the von Mises flow stress. Two elasto-plastic models were already implemented [13], namely the linear plastic and the Johnson-Cook models. Additionally, the Swift and Voce models were added to this list. Here, however, similarly to Leroch et al. [5], we use the empirical Johnson-Cook model for the von Mises flow stress scaled with damage[22]:

[TABLE]

where $\varepsilon_{p}$ is the equivalent plastic strain which is calculated jointly with the deviatoric stress tensor according to the algorithm developed by Leroch et al. [5] unless the Gurson-Tvergaard-Needleman (GTN) model is used (see Algorithm 1), $\varepsilon_{p}^{*}$ is the normalized plastic strain rate, $T^{*}$ the normalized temperature, $A$ the yield stress, $B$ and $n$ the strain hardening parameters, $C$ the strain rate parameters, and $m$ the temperature coefficient. The normalized plastic strain rate and temperature are respectively given by:

[TABLE]

where $\dot{\varepsilon}_{p}$ and $\dot{\varepsilon}_{0}$ are the plastic strain rate, and the user-defined reference plastic strain rate, $T_{0}$ is the reference temperature, and $T_{m}$ the melting temperature.

5 Damage and failure models

In order to determine the amount of damage in each particle, we have implemented in the LAMMPS SMD package[13] three different phenomenological models: Cockcroft-Latham (CL) for its simplicity [21], Johnson-Cook (JC) [22] and Gurson-Tvergaard-Needleman (GTN) for their wide application to damage and ductile solids [23]. CL and JC are models that only predict the damage onset point and need to be coupled with a damage evolution model, while the GTN model describes the nucleation, growth of voids until fracture occurs. Here two variables related to damage are used: the damage initiation variable $D_{init}$, and the damage variable $D$.

5.1 Cockcroft-Latham

Cockcroft-Latham [21] proposed that the onset of fracture be simply governed by the amount of “plastic work” per unit volume expressed as:

[TABLE]

where $\sigma_{1}$ is the maximum principal stress, and $\langle\sigma_{1}\rangle=\sigma_{1}$ if $\sigma_{1}>0$ and $\langle\sigma_{1}\rangle=0$ otherwise. Failure occurs when $W_{p}$ reaches a critical value $W_{cr}$ which means that the damage initiation variable is simply derived as:

[TABLE]

and indicates, when equal to 1, when damage starts.

5.2 Johnson-Cook

The Johnson-Cook failure model [22] is a widely used model by both scientists and in industry to simulate the damage of ductile materials. It is available in most of the commercial and open-source finite element solvers. It is a strain rate and temperature dependant phenomenological model based on the local accumulation of plastic strain. According to this model, damage starts when [36]:

[TABLE]

with

[TABLE]

and where $D_{1}$, …, $D_{5}$ are five material constants, $\sigma^{*}$ is the stress triaxiality.

5.3 Damage evolution model

Cockcroft-Latham and Johnson-Cook models describe only the point at which damage initiates. In order to have a complete model of the fracture phenomenon, a damage evolution model is required.

Here, for the sake of simplicity, it was arbitrarily assumed that the damage variable is given by:

[TABLE]

Even if this assumption affects the details of the damage propagation, it does not change the fundamentals of the implementation on which we focus here.

5.4 Gurson-Tvergaard-Needleman

The Gurson-Tvergaard-Needleman (GTN) model is one of the most referenced void-nucleation-growth-coalescence model. First, Gurson [37] developed a yield function for porous materials. Tvergaard and Needleman then extended it to fully quantitatively describe void nucleation and fracture [38]. This model supposes that the yield function $\Phi$ of a porous material with a void fraction $f$ is expressed as

[TABLE]

where $\sigma_{eq}$ the equivalent von Mises stress, $\sigma_{M}=\sigma_{f}(\varepsilon_{M}^{p},\dot{\varepsilon}_{p},T)$ is the flow stress of the undamaged matrix (the part of material free of void) calculated according to Equation 38 , and $f^{*}$ the void fraction function defined as[23]:

[TABLE]

with $f_{c}$ the critical void fraction, and $f_{f}$ the void fraction at failure. Failure is supposed to occur when the void fraction function $f^{*}$ equals $1/q_{1}$, thus the damage variable is defined as $D=q_{1}f^{*}$, and the yield function can be rewritten as:

[TABLE]

The evolution of the void fraction is governed by both the nucleation of new voids and the growth of existing ones such that the total void fraction growth rate is:

[TABLE]

where $\dot{f}_{nucleation}$ and $\dot{f}_{growth}$ represent the nucleation rate and the growth rate of voids, respectively.

Chu and Needleman have suggested that the nucleation rate follow a normal distribution [39]:

[TABLE]

where $\varepsilon_{M}^{p}$ is the undamaged matrix plastic strain, $\varepsilon_{N}$, $s_{N}$ are the mean and standard deviation of the distribution of plastic strain respectively, and $\mathrm{F}_{N}$ is the total void volume fraction that can be nucleated.

The growth rate of voids is comprised of two terms respectively representing that the matrix is plastically incompressible and the mechanism of void softening in shear [40]:

[TABLE]

where $\bm{\dot{\varepsilon}^{p}}$ is the macroscopic strain rate, and $\omega(\bm{\sigma})$ a non-dimensional metric defined as:

[TABLE]

which lies between [math] and $1$, with $\omega(\bm{\sigma})=0$ for all axisymmetric cases.

Owing to the use of a non-linear yield function, the plasticity algorithm proposed by Leroch et al. [5] needs to be adapted to the GTN model. The new algorithm developed to compute the plastic strain and the deviatoric part of the stress tensor is as follows (Algorithm 2). At a given time step $n+1$, as it is unknown if yielding has occurred, it is supposed, as a trial, that the total strain increment $\Delta\bm{\varepsilon}$ is purely elastic. The corresponding stress is therefore:

[TABLE]

where $\bm{C}$ is the stiffness matrix. Since the strain increment is small, the trial von Mises equivalent stress and the trial hydrostatic pressure can respectively be approximated as:

[TABLE]

where $G$ and $K$ are respectively the shear and the bulk moduli of the undamaged material.

If yield does not occur, $\bm{\sigma}_{trial}$ is exactly the stress tensor at the time step $n+1$, and there is no plastic strain increment. If yielding occurs, however, the yield point is exceeded and the material has deformed plastically. In order to calculate the increment of plastic strain in the undamaged matrix, $\Delta\varepsilon_{M}^{p}$, first, it is supposed that the void fraction $f$ and the stress triaxiality $\sigma^{*}$ are constant. Under this assumption, the yield function is only function of the ratio $x=\sigma_{n+1}^{eq}/\sigma_{M}^{p}$ and the value of $x$ at which yield occurs is solved numerically using the Newton-Raphson method. Second, since $\Delta\bm{\varepsilon}=\Delta\bm{\varepsilon}_{el}+\Delta\bm{\varepsilon}_{p}$. The stress tensor at the time step $n+1$ is:

[TABLE]

and its equivalent von Mises stress and hydrostatic pressure can be approximated as:

[TABLE]

where $\sigma_{y}=x\sigma_{M}^{p}$ is the yield stress. By taking the difference between 54 and 57 the equivalent plastic strain increment is obtained as:

[TABLE]

Similarly the difference between 55 and 58 yields

[TABLE]

which is linked to $\Delta\varepsilon_{p}^{eq}$ using the flow rule

[TABLE]

and gives:

[TABLE]

In order to be able to solve Equation 59, it is assumed that the rate of equivalent plastic work in the matrix material equals the macroscopic rate of plastic work such that [23]:

[TABLE]

Thus, solving for $\Delta\varepsilon_{M}^{p}$ and plugging it into Equation 59 yields

[TABLE]

Finally, this non-linear equation is numerically solved for $\Delta\varepsilon_{M}^{p}$ also using the Newton-Raphson method.

The final hydrostatic pressure $\sigma_{n+1}^{m}$ cannot be calculated directly from Equation 58, since, owing to the non-linear relationship between pressure and deformation, the bulk modulus $K$ is not explicitly known. Locally, a linear approximation of $K$ is obtained by taking the first order approximation of the ratio between Equations 58 and 57 and recalling that triaxiality is supposed constant:

[TABLE]

Therefore, the new hydrostatic pressure is:

[TABLE]

6 Case study: simulation of a series of tensile tests of smooth and notched specimens made out of three different Weldox steel alloys, and comparison with FEM results

In order to have confidence in the capabilities of the new damage implementation presented here, it was used to simulate a series of tensile tests of three different specimens each made out of three different Weldox steel alloys. The three different specimen geometries considered are: one smooth and two pre-notched cylindrical specimens with a notch radius of 2 mm and 0.4 mm, respectively (see Fig. 4). Three different Weldox steel alloys are considered: W460E, W700E, and W900E. This case study was chosen because experimental results and material parameters for the three different damage models implemented have been published in the literature [36, 41, 42].

For each combination of geometry and alloy, simulations are performed using the CL, the JC, and the GTN damage models as detailed in the previous sections. The results of the TLSPH simulations are then compared with not only the original experimental data published by Dey and coworkers [36, 41], but also with the results of the same simulations performed using FEM, taken here as a reference. However, the FEM simulations are performed only using the JC damage model as it is the only one of the three damage models to be implemented in the FEM software used here, i.e Abaqus Explicit. All the necessary material parameters are taken from the literature. The elasto-plastic parameters as well as the parameters for the CL and JC criteria are those proposed by Dey and coworkers [41, 36] (see Tables 1, 2 and 3), while the parameters for the GTN model are taken from Øien and Schønberg [42] (see Table 3).

In all simulations, the full gauge length of the smooth specimens is simulated, while, to minimize computational times, for the pre-notched specimens, only a portion of respectively 7 mm and 20 mm in height around the notch is considered. In SPH, all domains are discretized using a simple cubic array of particles which constant spacing is of 0.15 mm and 0.3 mm for the smooth and $R=2~{}mm$ notched specimen, respectively, while for the $R=0.4~{}mm$ notched specimen it varies from 0.8 mm for furthest point from the notch, to 0.1 mm in the notch. On the last top and bottom rows of particles, a velocity $\pm v$ of the form $v_{max}(1-e^{-t})$ is applied in the direction colinear to the each specimen’s axis (Fig. 4). As no strain rate effect is considered here, to speed-up simulations, a high maximum velocity is chosen: $v_{max}=10~{}mm/ms$. Such a high velocity was chosen in order to decrease the total required time to be simulated and has no impact on the results as no strain rate effect was taken into account.

6.1 Results

Fig. 5 presents the results of the SPH simulations for all combination of specimen geometries and alloys tested plotted alongside both the experimental data and FEM results. More specifically, Fig. 5a shows the evolution of the engineering stress as a function of the engineering strain for the smooth specimens, while Fig. 5b and 5c show the variation of principal stress vs. true plastic strain averaged over the cross-section at the middle of the notch of both pre-notched specimens.

First, it can be seen that, the flow stresses predicted by the SPH simulation is in good agreement overall with the FEM. The agreement is very good in the case of the smooth and $R=2~{}mm$ notched specimens, while for the $R=0.4~{}mm$ notched specimens not as good agreement is obtained. In that case, the results of the SPH simulations fit better the experimental data than FEM. In fact, FEM consistently overpredicts the flow stresses for all three alloys. This results is, however, consistent with the results of the FEM simulations performed by Dey and coworkers [36, 41]. It can be also noticed that for all these simulation cases, the flow stress given by SPH using the GTN model is slightly lower than that obtained with the other damage models. This is due to the fact that GTN uses a different yield function and that presence of voids in the underformed material is considered, i.e. $f_{0}\neq 0$.

Second, as expected, after damage initiates, a steady decay in stress is observed, and the simulations are numerically stable. The strains at failure predicted by SPH using the JC model are in very good agreement with those predicted using FEM in the case of the smooth and $R=2~{}mm$ notched specimens. In the case of the $R=0.4~{}mm$ notched specimens, however, SPH simulations predict strain at failure that are higher than those predicted by FEM, but they are a little closer to experimental results. Such discrepancies can be explained by the lowest flow stress observed in the SPH simulations. The strains at failure predicted by SPH using the CL model, however, are always higher than that using JC. This is consistent with conclusions from Dey and coworkers’ work [36, 41].

Since the results obtained using the Johnson-Cook damage model using SPH and FEM are in good agreement with each other, we have confidence that the damage implementation in SPH is successful. Also, we could see that, none of the models are quantitatively able to predict the experimental results. This was expected as all the parameters used were taken directly from the literature and efforts were not made to optimise them from the point of view of agreement with experiment. Of course, this is an important aspect of using SPH to simulate real experimental data but that is not the focus of this contribution, which is on presenting this new implementation of damage and failure in SPH.

A typical example of the result of a TLSPH simulation is given by Fig. 6. It shows, alongside the stress-strain curve, the evolution of distribution of both the equivalent von Mises stress and damage accumulated in a W700E notched cylindrical specimen which notch radius is of $R=0.4~{}mm$ using the JC damage criterion. At the onset of plasticity (snapshot (a)), it can be seen that, as expected, the stress is concentrated around the notch tip. Thus, it is natural to see that when the strain at failure is reached (snapshots (b)) that cracks initiate from this same point. One can also clearly distinguish from snapshot (c) that the two forming cracks are oriented at $\pm 45^{\circ}$ with the horizontal, i.e. in the direction of the principal stresses, and that their presence lead to a decrease of maximum principal stress. The stress relaxation in the outer parts of the specimen that follow cause a change in the crack growth direction (snapshot (d)) as they start propagating horizontally. This further growth lead to stress at the center of the specimen along a horizontal plane. There, another crack initiates (snapshot (e)) and quickly propagates outwards to connect with the upper crack and break the specimen into two (snapshot (f)).

7 Discussion

For the simulations of mechanical problems involving large deformations, particle based methods present the advantage over Finite Elements Modeling to not require re-meshing which is very CPU intensive, and is not always guaranteed to be successful. When using particle based methods for such applications, however, one needs to be aware of the risk of numerical failure that would occur when the distance between two neighbouring particles becomes eventually higher than the the kernel cutoff radius. This can lead to failure strains and energy releases that are not physical. This is a reason why TLSPH is more favourable than classical SPH for simulations where their is a probability of failure. In fact, in TLSPH, if the reference configuration is never updated, two neighbouring particles would remain as such whatever their separation distance in the deformed state, thus preventing numerical failure from occurring.

Using TLSPH with the “pseudo-spring” method for the simulation of ductile materials generates numerical instabilities as soon as damage initiates. This triggered the need for a new way of taking damage into account in the SPH approximation. The resulting new implementation was tested on a series of tensile tests for different damage criteria. The results of these tests show that, due to the good agreement with FEM and experimental results, the new TLSPH approximation can be trusted. They also show that, depending on the type of specimen and elasto-plastic parameters used, some damage criteria give better results than others. The Johnson-Cook, Cockcroft-Latham, and Gurson-Tvergaard-Needleman damage criteria have been implemented as they are one of the most popular criteria used in academia and industry. These changes have all been implemented in the user package SMD of LAMMPS for the community to test and use. The code is available in the first author’s Github account: https://github.com/adevaucorbeil/lammps/ .

8 Conclusion

A new method for damage and failure was implemented in the LAMMPS user package Smooth Mach Dynamics for the simulation of ductile solids using Total-Lagrangian SPH. Three popular damage criteria that are Cockroft-Latham, Johnson-Cook, and Gurson-Tvergaard-Needleman have also been implemented in this package. This new method was then tested on the study of tensile tests of smooth and notched specimens made out of three different Weldox steel alloys. The simulation results obtained using TLSPH were then compared to both experimental and FEM results. These results shows that TLSPH can be used successfully to simulate the mechanical response of materials experiencing damage and failure without relying on or exploiting numerical failure. The simulations presented here serve as validation of the implementation and are the first necessary step before being able to use this model for the simulation of more complicated phenomena such as machining, wear, and impacts of ductile materials.

Declaration of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors gratefully acknowledge the financial support of the Australian Research Council (ARC) Training Centre in Alloy Innovation for Mining Efficiency (IC160100036). Fruitful discussions with Prof. Matthew Barnett (Deakin University and Director of the ARC Training Centre) are acknowledged.