Numerical modelling of crack initiation, propagation and branching under dynamic loading

Engineering Fracture Mechanics xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Numerical modelling of crack initiation, propagation and branching under dynamic loading ⁎

Md Rushdie Ibne Islama,b, , Amit Shawb a b

ESS Engineering Software Steyr GmbH, Berggasse 35, 4400 Steyr, Austria Department of Civil Engineering, Indian Institute of Technology, Kharagpur 721302, West Bengal, India

A R T IC LE I N F O

ABS TRA CT

Keywords: Pseudo-spring SPH Crack initiation Propagation Branching Dynamic fracture

In this paper crack initiation, propagation and branching phenomena are simulated using the Pseudo-Spring Smoothed 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. 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. Next, the effect of notch location for a plate with a circular hole is studied. The ability of the framework to model arbitrary crack paths and surfaces are also demonstrated via three-dimensional simulations of chalk under torsion, Kalthoff-Winkler experiment, Taylor bullet impact and crack branching.

1. Introduction Numerical modelling of crack initiation, propagation and branching under dynamic loading is still a challenging problem. Over the last few decades, several mesh-based and mesh-free methods are explored to simulate the propagation of arbitrary oriented cracks. Each of these methods needs some special treatment to model the crack discontinuity. Finite Element Method (FEM) is used to model crack path and its propagation by element deletion or kill method [1,2] or by adaptive re-meshing [3–6]. In the element deletion technique [1,2], an element is assumed to have failed once the stress, strain or damage values in the element exceed a predefined threshold. The edges of the failed elements act as the new crack surfaces. The difficulty in such approaches is the numerical instability arises due to excessive element 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–6], this technique is not quite useful for tracking multiple cracks. The method of remeshing 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–11]. Therein, 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)

⁎

Corresponding author at: ESS Engineering Software Steyr GmbH, Berggasse 35, 4400 Steyr, Austria. E-mail addresses: [email protected] (M.R.I. Islam), [email protected] (A. Shaw).

https://doi.org/10.1016/j.engfracmech.2019.106760 Received 3 May 2019; Received in revised form 21 October 2019; Accepted 1 November 2019 0013-7944/ © 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Md Rushdie Ibne Islam and Amit Shaw, Engineering Fracture Mechanics, https://doi.org/10.1016/j.engfracmech.2019.106760

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Nomenclature K cs cf T D S E f rij A D1 − D5 (A¯ , B¯ , n¯, W umax Tm q j m i Wp x P T0 J2 iD

Ui h C e B d dt

¯ 

u

yf T∗ ∊pl π (β1, β2) Δp σ δt δtc ∊∗̇ pl (α, β ) ∊f Δ ∊pl Δt σy σm θ ρ Δx ∊̇pl ν ∊̇0 μ ω ∊̇ σ∗ CFL FFG XEFG EFEM XFEM FEM LSM MLS SPH 2D 3D

Bulk modulus CFL number Coefficient for return mapping Current temperature Damage index Deviatoric component of Cauchy stress tensor Elastic modulus Interaction coefficient Inter particle distance between particles i and j Inverse of symmetric re-normalization matrix Johnson-Cook damage parameters C¯ , m ¯ ) Johnson-Cook material constants Kernel function Maximum velocity Melting temperature Normalised position vector Particles inside the support domain of particle of interest Particle mass Particle of interest Plastic work Position vector Pressure Room temperature Second invariant of deviatoric stress tensor Set of particles connected by partially damaged springs Set of particles connected by undamaged springs Smoothing length Sound speed Specific internal energy Symmetric re-normalization matrix Time derivative in moving Lagrangian frame Total number of particles in the computational domain Velocity vector

Von-Mises yield function Homologous temperature Accumulated effective plastic strain Artificial viscosity Artificial viscosity controlling parameters Average particle spacing Cauchy stress tensor Computed stretch Critical stretch Dimensionless plastic strain rate Dummy indices Fracture strain Incremental effective plastic strain Incremental time step Material flow stress Mean stress Measure of angle Particle density Particle spacing in current configuration Plastic strain rate Poison’s ratio Reference strain rate Shear modulus Spin tensor Strain rate Stress triaxiality Courant-Fredrich-Levy Element-Free Galerkin (FFG) Extended Element-Free Galerkin (XEFG) Embedded Finite Element Method (EFEM) Extended Finite Element Method (XFEM) Finite Element Method (FEM) Lattice Spring Models (LSM) Moving Least Square (MLS) Smoothed Particle Hydrodynamics (SPH) Two Dimension Three Dimension

[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–18] developed the extended finite element method (XFEM) where the elements are locally enriched through a partition of unity [19]. In XFEM, only the surrounding elements of a crack and the crack propagation path (known a priori) are enriched and the crack propagates through these enriched elements (not only through the element boundaries). Though the element based methods can model the evolving cracks and track their paths, the requirement of re-meshing or local enrichments makes these numerically 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 or particle methods are introduced over the years. Belytschko and his colleagues [20–23] developed the Element-Free Galerkin (FFG) method for static and dynamic crack propagation. As the method is mesh-free, there is no remeshing required. In their work, they used moving least square (MLS) 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 introduced 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. Alternatively, Peridynamics, a non-local approach to continuum mechanics, is proposed in [28–30] to model crack propagation and material fracture. Peridynamics uses an integral-differential equation where the partial derivatives of the displacement are replaced by integrations. As no explicit derivatives of field variables are involved, the Peridynamics formulation is suitable to model both continuous and discontinuous problems. Several attempts are made to simulate crack propagation by Peridynamics [31–37]. The Lattice Spring Models (LSM) is an interesting concept to model solid deformation [38,39]. The LSM is a discontinuous-based approach in which lattice points connected by springs represent the materials. As LSM is a discontinuous based approach it does not 2

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

require any special treatment for material damage. The springs can be removed once certain field variables exceed their predefined values to represent fracture. However, the standard LSM can account for only a certain range of Poisson’s ratio and elastic problems. Further improvements such as the Born spring model, the beam element model, the multibody shear spring, and the non local potential are made in LSM [40–43] to overcome these limitations. However, these developments are mathematically complex and computationally intensive. The smoothed particle hydrodynamics (SPH) [44] is relatively simpler and truely particle-based formulation. SPH is initially introduced for astronomical problems. This is extended to fluid problems and later to solids [45]. 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 chosen 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. To overcome this, a Total Lagrangian SPH shell formulation is proposed in [46]. Therein the neighbourhood of a particle does not change through the computation and all the particles inside the neighbourhood interacts through links. These links between interacting particles is assumed to be broken once the fracture criteria are met based on the initial and current distances between particles. Once a bond is broken, the connecting particles through that bond do not interact with each other, and the generation of crack plane/surface occurs. Subsequently, the kernel functions are updated. However, for the problems with large material deformation, the Lagrangian kernel may suffer from difficulties associated with large deformation and discontinuities in the computational domain due to material fracture. The damage of materials and formation of multiple crack paths are modelled in SPH by the use of moving least square (MLS) functions [47]. The damage evolution in the materials is modelled through elasto-plastic damage models. The computation of the crack path is quite similar to the ‘cracking particle method’ [24]. The cracks are represented by tracking the ruptured representative elementary volumes (REVs). This framework requires tracking of the cracked REVs and a visibility criterion to remove the interaction between cracked REVs. Once the crack paths are formed, the MLS functions are also modified. Though this method is able to capture crack paths and their propagation with reasonable accuracy, the requirement of tracking cracked REVs and the use of visibility criterion make it computationally intensive. In order to facilitate SPH to model crack propagation without any numerical artifact, the ‘pseudo-spring’ strategy is developed in [48]. Herein all the interacting particles in the computational domain are assumed to be virtually connected through a network of springs. 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. Herein, the 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. The broken spring indicates 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 broken pseudosprings. The simplicity of the SPH is preserved in the computation. The pseudo-spring SPH method is successfully used to simulated impact induced damages [49–54]. Though LSM has certain similarities with the pseudo-spring SPH model (both use spring network), their formulation and implementation are quite different. LSM is a bottom-to-top (micro-scale) model in which the springs represent the material properties, and the springs of various types (lateral and angular) are used to mimic the solid properties. However, the pseudo-spring SPH is a continuum based framework where the springs are used to represent the discontinuity developed during computation only. In this work, the 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, deviation of crack paths under different conditions and generation of multiple crack branches are simulated. The effects of tensile loading amplitude on crack branching are discussed. The generation of arbitrary crack surfaces under dynamic loading in three dimensions is modelled. The efficacy of the method is demonstrated via few benchmark 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 3

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

can be found in [45,55–57,48,52,53], and the 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

dρ ∂uβ = −ρ β dt ∂x

(1)

dρi = dt

(2)

¯ ij, β mj uijβ W

∑ j ∈ N¯ i

duα 1 ∂σ αβ = dt ρ ∂x β duiα = dt

∑ j ∈ N¯ i

(3)

αβ σ jαβ ⎛σ ⎞ ¯ mj ⎜ i 2 + 2 − πij δ αβ − Pija δ αβ ⎟ W ij, β ρi ρj ⎝ ⎠

(4)

de σ αβ ∂uα = dt ρ ∂x β dei 1 =− dt 2

(5)

∑ j ∈ N¯ i

αβ σ jαβ ⎛σ ⎞ ¯ mj uijα ⎜ i 2 + 2 − πij δ αβ − Pija δ αβ ⎟ W ij, β ρi ρj ⎝ ⎠

(6)

where, i denotes the particle of interest, j denotes its neighbour particles within the influence radius and i ∈ + is the set of neighbour particles within the influence domain. The density of the material is represented by ρ;uα is the velocity component; σ αβ represents the Cauchy stress components; e is the specific internal energy and α, β are the indices for coordinate axes. The time ¯ ij, β = 0.5(Wiĵ , β + Wjî , β ) is the symmetric form [58] of the gradient derivative is in moving Lagrangian frame and uijα = uiα − ujα . W corrected kernel Wiĵ , β . Wiĵ , β is calculated as [59]

Wiĵ , β = B βαWij, α

(7)

where

B = A−1, and Aβα = − ∑ j

mj ρj

x ijβ Wij, α

The kernel gradient is represented as Wij, β = following cubic B-spline kernel is used. 3

(8) ∂W (xi − xj, h) ∂x iβ

with h being the smoothing length of the kernel function. Here, the

3

⎧1 − 2 q2 + 4 q3, if 0 ⩽ q ⩽ 1 ⎪ W (q, h) = αd 1 (2 − q)3, if 1 ⩽ q ⩽ 2 ⎨4 ⎪0 otherwise ⎩ 10

(9) |xi − xj |

1

where, αd = 7πh2 and αd = 3 in 2D and 3D respectively, and q represents the normalised position vector as q = h . πh The artificial viscosity πij shown in Eq. (10) is utilised in the current computation to remove any unphysical oscillation present near any jump or shock [44]. 2

¯

⎧ −β1 Cij μij + β2 μij , ifu . x ⩽ 0 ij ij ρ¯ij πij = ⎨ 0, otherwise ⎩

(10)

The intensity of the artificial viscosity is determined using β1, β2 . μij = computed as C =

E ρ

h (uij . xij ) xij2+ ∊ h2

and ∊ = 0.01. The sound speed through the medium is

, E denotes the Young’s modulus, xij = xi − xj , uij = ui − uj, C¯ij = 0.5(Ci + Cj ) and ρ¯ij = 0.5(ρi + ρj ) . For pro-

blems where the structure do not experience any fracture, the optimum values of the artificial viscosity parameters may be obtained through the approach given in [60–62]. In the present study as the material undergoes fracture, the parameters β1, β2 are determined through numerical experiments such that πij removes the unphysical oscillations across the jump in the field variables but do not overdamp the system. The artificial/Monaghan pressure [63] Pija as shown in Eq. (11) is used to remove the tensile instability.

|Pj | ⎞ W (Δx , h) ⎤l ⎛ |P | Pija = γ ⎜ i2 + 2 ⎟ ⎡ W (Δp , h) ⎥ ρ ρj ⎢ ⎦ ⎠⎣ ⎝ i

(11)

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

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

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

ρ P = K ⎛⎜ − 1⎞⎟ ⎝ ρ0 ⎠

(12)

where, K , ν and ρ0 are the bulk modulus, Poisson’s ratio and initial mass density of the material respectively. The Cauchy stress components σ αβ are calculated from the deviatoric stress components S αβ and hydrostatic pressure P as σ αβ = S αβ − Pδ αβ . The rate of change in deviatoric stress components is calculated by Jaumann rate as,

Fig. 1. (a) Pseudo-springs network on immediate neighbours and crack propagation through damaged springs (b) Crack propagation through Pseudo-spring networks in three dimension. 5

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

̇ = 2μ ⎛∊̇αβ − 1 δ αβ ∊̇γγ ⎞ + S αγω βγ + S γβωαγ Sαβ 3 ⎝ ⎠

(13)

where, μ represents the shear modulus. The components of strain rate (∊̇ ) and spin tensor (ω ) are computed as,

∊̇αβ =

1 ⎛ ∂uα ∂uβ ⎞ + 2 ⎝ ∂x β ∂x α ⎠

(14)

ωαβ =

1 ⎛ ∂uα ∂uβ ⎞ − 2 ⎝ ∂x β ∂x α ⎠

(15)

⎜

⎟

⎜

⎟

Von-Mises yield function yf =

J2 −

σy 3

1

is used to account for the effects of metal plasticity. σy and J2 = 2 S αβS αβ are the yield

stress of the material and the second invariant of deviatoric stress tensor respectively. The Wilkins criterion return mapping of the yield surfaces and cf = min

(

accumulation of plastic work density are estimated as

Δ ∊αβ pl =

1 − cf 2μ

σy 3J2

Snαβ

= cf S αβ is used for

)

, 1 . The plastic strain increment, effective plastic strain increment and

S αβ

(16)

Fig. 2. Modification of kernel in 2D based on damage state [52,53]. 6

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Δ ∊pl =

1 − cf 2 αβ Δ ∊αβ pl Δ ∊ pl = 3 3μ

3 αβ αβ S S 2

(17)

αβ ΔWp = Δ ∊αβ pl Sn

(18)

The standard predictor-corrector method [64] is used to update the discretised Eqs. (2), (4), (6) and the stable time step for

(

h

)

integration is determined from CFL (Courant-Fredrich-Levy) condition as, △t = mini cs C + i| v | where, CFL number cs is 0.3, Ci is i

i

elastic wave velocity, hi is smoothing length and vi 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 [48,52,53]. 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 Eq. (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 [52]. 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 Fig. 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 they 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 [65]. The concept is illustrated in Fig. 1. The analogy is incorporated in the framework in the following way. An interaction factor fij between a connecting pair of particles is defined as fij = 1 − Dij where Dij 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 Dij = 0 and this leads to the interaction factor fij = 1 implying complete interaction between particle pairs. As the damage accumulates in the springs, the interaction factor fij decreases i.e. Dij ∈ (0, 1) ⇒ fij ∈ (0, 1) which implies partial interaction. When the damage index Dij = 1 i.e. the interaction factor fij = 0 , interaction stops leading to formation of crack surface. Dij is computed as Dij = 0.5(Di + Dj ) , where Di and Dj are the damage parameters at the i − th and j − th particles respectively. Let iD = {j ∈ i |0 ⩽ fij < 1} ⊆ i and Ui = i ⧹iD be the set of immediate neighbours with which the i-th particle pairs through damaged and undamaged springs respectively at a given time instant. The interaction between particles i and j ∈ iD in a damaged ¯ ij, β = fij W ¯ ij, β ∀ j ∈ iD ⊆ i . Fig. 2 configuration is reduced by modifying the corresponding derivative of the kernel function as W shows the typical modified kernel in 2D for different damage states i.e., fij = 1 (no damage), 0 < fij < 1 (partial damage) and fij = 0 (open crack). Finally the modified forms of the discrete conservation Eqs. (2), (4), (6) become,

dρi = dt

∑ j ∈ NUi

¯ ij, β + mj uijβ W

∑

¯ ij, β ) mj uijβ (fij W

i j ∈ ND

(19)

Fig. 3. Setup of pre-notched glass plate under tensile loading. 7

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

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

Artificial Viscosity

E (GPa)

ν

∊max

Δp (mm)

h (mm)

β1

β2

(kg/m3 ) 2450

32

0.2

0.000509

0.125

0.250

1.0

1.0

Parameter

ρ

Value

Discretization

Fig. 4. Crack propagation and branching: Experimental observation [69].

Fig. 5. Crack path evolution over time in pre-notched plate under tensile loading of σ = 1 MPa.

duiα dt

αβ

∑

=

j ∈ NUi

σ mj ⎛⎜ i 2 + ρ ⎝ i αβ

+

∑ i j ∈ ND

σ mj ⎜⎛ i 2 + ρ ⎝ i

σ αβ j ρj2

σ αβ j ρj2

¯ ij, β − πij δ αβ − Pija δ αβ ⎞⎟ W ⎠

¯ ij, β ) − πij δ αβ − Pija δ αβ ⎟⎞ (fij W ⎠

(20)

8

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 6. Comparison of final crack path under tensile loading of σ = 1 MPa: (a) [70] (b) [68] (c) [71] and (d) Present. 2500

Crack velocity (m/s)

2000

Rayleigh speed Song et al. 2008 (XFEM) Song et al. 2008 (CZM) Song & Belytschko 2009 Agwai et al. 2011 Present study

1500

1000

500

0 0

10

20

30 Time (μs)

40

50

60

Fig. 7. Comparison of crack speed over time.

Fig. 8. Crack path evolution over time in pre-notched plate under tensile loading of σ = 0.3 MPa.

9

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 9. Crack path evolution over time in pre-notched plate under tensile loading of σ = 1.1 MPa.

Fig. 10. Crack path evolution over time in pre-notched plate under tensile loading of σ = 1.2 MPa.

Fig. 11. Crack path evolution over time in pre-notched plate under different magnitude of tensile loading.

dei dt

−

1

= −2 1 2

∑ i j ∈ ND

∑ j ∈ NUi

mj uijα ⎛⎜ ⎝

mj uijα ⎛⎜ ⎝

σiαβ ρi2

σiαβ

+

ρi2

+

σ αβ j ρj2

σ αβ j ρj2

¯ ij, β − πij δ αβ − Pija δ αβ ⎞⎟ W ⎠

¯ ij, β ) − πij δ αβ − Pija δ αβ ⎞⎟ (fij W ⎠

(21)

In the present formulation, the particles remain intact. Only the pseudo-springs undergo continuous damage (Fig. 1) resulting in continuous change in particle interaction. If the interaction between two particles ceases to exist due to the formation of a crack between them (Dij = 1), the contact and sliding between crack surfaces need to be explicitly established. Towards this, a simple yet efficient particle-particle contact approach [66] is used in this study. In the computation, particles with broken pseudo-springs are identified as particles on crack surface. During the computation of the interaction between two particles, it is first checked if the pseudo-spring between these two particles exists. If there is still a pseudo-spring (Dij < 1), it means that there is no crack surface between these two particles. Hence, the particle pair interact through the kernel function. If the pseudo-spring between these two particles does not exist any more (Dij = 1) due to damage, then the interaction between these two particles (when they come inside each other’s influence domain due to crack closer) or with other particles on the crack surface (in case of sliding) are handled through the particle-particle contact approach [53]. 2.2. Stress boundary condition A stress boundary condition for SPH is developed in the context of solid continua in [67]. The boundary condition can be applied 10

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 12. Setup of pre-notched plate with a circular hole under tensile loading. Table 2 Parameters for pre-notched plate with a circular hole under tensile stress. Material properties

Artificial Viscosity

E (GPa)

ν

∊max

Δp (mm)

h (mm)

β1

β2

(kg/m3 ) 2700

71.4

0.25

0.00483

0.1

0.15

1.0

1.0

Parameter

ρ

Value

Discretization

with or without ghost/dummy particles. The modified form of the momentum equation to account for the applied stress is shown in Eq. (22).

duiα du α =⎛ i ⎞ − dt ⎝ dt ⎠Equation 20 ⎜

⎟

αβ

∑ j∈Ni

⎛σ ⎞ ¯ mj ⎜ 0 ⎟ W ij, β ρρ ⎝ i j⎠

(22)

Here, σ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. 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. The computed crack speeds, their paths, and crack surfaces are compared with experimental and numerical results available in the literature. Then crack path from a notch in a plate with a circular hole is also modelled. It is to be noted that in SPH the h/Δp ratio (h is the smoothing length and Δp is the inter-particle spacing) is an important parameter. The effect of particles spacing and the convergence of the predictions via the pseudo-spring SPH framework is extensively studied in [53] in the context of adiabatic shear plugging failure in metal plate. It was observed that the predictions converge after a certain particle resolution with a fixed h/Δp ratio. With the converged particle spacing, the failure plane was found to be consistent with the experimental studies. Therefore, a convergence study (like FEM analysis) with a constant h/Δp ratio is important to find the inter-particle spacing for any simulation. This particle spacing then can be used for any predictive simulation. In the present work, similar type of approach is adopted. The results shown in the examples are obtained with the converged interparticle spacing found through the numerical experiments. 11

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 13. Comparison of crack path for notch distance 5 mm from support (case I).

12

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 14. Comparison of crack path for notch distance 10 mm from support (case II).

13

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 15. Comparison of crack path for notch distance 15 mm from support (case III).

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 14

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 16. Diagram of chalk under torsion. Table 3 Parameters for the chalk under torsional loading. Material properties Parameter

Value

ρ

E

(kg/m3 )

(GPa)

1150

2

ν

0.18

Discretization

Artificial Viscosity

p σmax (MPa)

Δp

h

(mm)

(mm)

15

0.4

0.52

β1

β2

1.0

1.0

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 Fig. 3. The plate is under tensile loading (σ = 1 MPa) at the boundaries along the length of the plate. The material parameters used for the present computation are from [68,33]. 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.

1, ∊ ⩾ ∊max D=⎧ ⎨ ⎩ 0, Otherwise

(23)

In the experiment [69], it is observed that initially the crack propagates in a self-propagating mode and after few microseconds, the branching takes place as shown in Fig. 4. The predicted crack paths (crack initiation, propagation and branching) as obtained via the present simulations are shown in Fig. 5. The crack propagation initiates around 10 μs and propagates in a straight line up to 27 μs. The initiation of branching of crack is observed at 28 μs. The crack propagates in two symmetrical branches and reaches the plate boundary at 62 μs. The crack branching problem is also studied in [9,18,24,33,71]. The final crack path obtained from the present simulation is compared with other numerical predictions available in the literature in Fig. 6. The crack propagation speed over time is compared with different techniques (see Fig. 7). The present results compare reasonably well with the results obtained in [33] using peridynamics, in [68] using XFEM-CZM, and in [72] 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 (σ ) intensities on branching behaviour is explored. σ 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 (σ ) intensity. For, σ = 0.3 MPa loading, no branching is observed, and the crack propagates in a straight line as shown in Fig. 8. But the branching can be observed for all the other load intensities. The crack branching starts around 15 μs for σ = 1.1 MPa (Fig. 9), 14 μs for σ = 1.2 MPa (Fig. 10), 9 μs for σ = 2.0 MPa and 6 μs for σ = 4.0 MPa (Fig. 11). As the loading amplitude (σ ) 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 σ = 1.1 MPa but the secondary branches are found for σ = 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 [37]. 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 [73] and the microcracks in this zone undergo nucleation, growth, and coalescence. These lead to crack path instabilities [74]. 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 [75]. This leads to the multiple branches of crack paths for the higher amplitude of boundary stress. 15

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

(a) Crack surface by [80]

(b) Predicted crack surface by [80]

(c) Present prediction

(d) Present prediction

Fig. 17. A three dimensional chalk under torsion with helical crack surface.

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 Fig. 12. One end of the plate is clamped, and the another end is under uniform tensile stress σ . This numerical study aims to simulate the crack propagation originating from the pre-existing notch and to track the changes in crack 16

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 18. Set up for Kalthoff-Winkler experiment.

Table 4 Parameters for Klathoff-Winkler experiment simulation. Material properties Parameter

Value

Discretization

Artificial Viscosity

E (GPa)

ν

δtc

Δp (mm)

h (mm)

β1

β2

(kg/m3 ) 8000

190

0.3

0.0044

0.8

1.04

1.0

1.0

ρ

path due to the presence of the off-centre circular hole. The problem is studied in [76] by the arbitrary local mesh replacement method (a finite element based strategy). Tabiei and Wu [77] used DYNA3D and Kosteski et al. [78] used truss-like discrete element method to model the crack paths. Dipasquale et al. [34] used adaptive grid refinement in 2D peridynamics to track the crack paths. Here, the same example is simulated using the pseudo-spring SPH. The material and SPH parameters for the simulation are given in Table 2. The damage parameter D is obtained at any particle via Eq. (23). Three different notch locations which referred to as case I (h = 5 mm), case II (h = 10 mm) and case III (h = 15 mm) [79] are considered. Fig. 13 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 Figs. 14 and 15 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 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 [78,79] as shown in the Figs. 13–15. 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. 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 a 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. 17

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

(a)

(b)

(c)

(d) Damage:

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 19. Damage contour in the crack path for Kalthoff-Winkler experiment.

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 Fig. 16. 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 umax = 28.3 m/s in the tangential direction to the chalk bar as shown in Fig. 16. The distribution of the velocity field is given as

u = umax sin(

θ+ 2

Π 2

)

(24)

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 [80,81]. The material properties and other computational parameters are given in Table 3. The maximum principal stress criterion is used for modelling failure. For a given spring i − j the maximum principal stress σijp is determined as σijp = (σip + σjp)/2, where σip and σjp are the maximum principal stresses at i-th and j-th particle respectively. The springs p are considered to fail completely when σijp exceeds σmax value (Table 3). 18

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 20. Crack propagation for Kalthoff-Winkler experiment.

The failure can happen anywhere and results in an arbitrary crack surface as shown in Fig. 17. The present fracture surfaces are found to be helical. These surfaces are consistent with the crack surface obtained in [80,81] 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 [82] is considered. A cylindrical projectile impacts a plate with two edge notches as shown in Fig. 18. Numerical simulations of similar problems are performed in [70,26,71,78,37,48]. 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 [48,71] criterion. The critical stretch at any time instant is calculated as

δt =

rij ∣t − Δp Δp

(25) 19

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 21. Cylindrical projectile impacting rigid surface. Table 5 Computational parameters used for Taylor bullet impact. Inter-Particle Spacing (Δp )

Smoothing Length (h)

Time-step (Δ t)

Artificial Viscosity Parameters (β1, β2 )

0.3 mm

0.6 mm

5 × 10−9 s

(1.0, 1.0)

Fig. 22. Initial particle distribution of the projectile.

When the calculated stretch (δt ) at any time instant exceeds the critical stretch (δtc ) value, the material is assumed to be fully damaged. The crack pattern using the current framework is shown in Figs. 19 and 20. 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,37,78,34]. 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 Fig. 21. The failure of the projectile is modelled in this section. The bullet is made of Weldox 460E steel. The Johnson-Cook material [83] and damage models [84] are used. In the Johnson-Cook damage model the fracture strain is determined through five parameters as

∊f = [D1 + D2exp(D3 σ ∗)][1 + ∊pl̇ ∗]D4 [1 + D5 T ∗]

(26) 20

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 23. Petalling fracture mode for cylindrical projectile.

Here, D1-D5 are the material constants, σ ∗ is the stress tri-axiality ratio defined as plastic strain (∊pl̇ ∗ = as,

D=

∑

∊pl̇ ∊̇0

σm σy

1

and σm = 3 σii . ∊pl̇ ∗ is the dimensionless rate of

; where, ∊pl̇ and ∊0̇ are the current and reference strain rate respectively). The damage index (D) is computed

Δ ∊pl ∊f

(27)

where Δ ∊pl represents the incremental plastic strain. The Johnson-Cook material model is used to compute the yield stress σy 21

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 24. Petalling failure for cylindrical projectile at different particle spacing. ¯ σy = [A¯ + B¯ ∊npl ][1 + ∊pl̇ ∗]C¯ [1 − T ∗m¯ ]

(28)

T−T ¯ ) are the material constant. T ∗ is the homologous temperature (T ∗ = T − T0 ); where T , T0 and Tm are the current, where, (A¯ , B¯ , n¯, C¯ , m m 0 room and melting temperature respectively. The material and damage parameters may be found in [53]. 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 Fig. 22. The petalling fracture mode for the projectile is simulated in [85]. The present framework also predicts petalling fracture mode as shown in Fig. 23. The effect of particle spacing (Δp) is studied for this case. In SPH, the convergence depends on inter-particle spacing (Δp) when the h ratio is kept constant [86]. The Taylor bullet impact is simulated with two other particle spacing (Δp = 0.4 mm and 0.6 mm)

Δp

with same

h Δp

ratio. The petalling failure modes of the cylinder at different Δp is shown in Fig. 24.

3.2.4. Crack branching in 3-D As the last example, crack branching in 3-D is shown here. The predicted crack path and branching is shown in Fig. 25. Initially, the crack propagates in a straight line. After the self-propagation, the crack branches into two crack paths and these branches propagate towards the end of the plate. The crack branching behaviour obtained via the 3D formulation is found to be in good agreement with the 2-D results discussed in Section 3.1.1. 4. Closure Numerical modelling of crack initiation, propagation and branching is a challenging problem. Towards this, though several attempts are made, not all are proved to be very efficient. Most of the existing methods, either mesh-based or mesh-free, require computationally intensive task such as local enrichment, re-meshing, element deletion, etc. 22

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Fig. 25. Crack branching in 3D.

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. The criterion or measure to be adopted for breaking the springs depends on the material behaviour. The ‘pseudo spring’ strategy, as a concept, is not restricted to any specific failure measure. 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 predicted crack surface is found to be consistent with the experimental and other 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 results reported in the literature. The crack branching in 3-D is also found to be consistent with the 2-D results. The ‘pseudo-spring’ strategy is an efficient tool for modelling 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. However, the ‘pseudo-spring’ analogy is an add on to SPH method. Hence, the outcomes of the simulations depend on the inter-particle spacing for a constant h ratio to some extent. Therefore, the Δp

present examples are simulated at different particle spacing and the results at converged Δp is reported here. 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.

23

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

Acknowledgement The authors would like to acknowledge the Defence Research and Development Organization (DRDO), India for their support of this work. References [1] Beissel S, Johnson G, Popelar C. An element-failure algorithm for dynamic crack propagation in general directions. Eng Fract Mech 1998;61(3):407–25. [2] Børvik T, Hopperstad O, Berstad T, Langseth M. Numerical simulation of plugging failure in ballistic penetration. Int J Solids Struct 2001;38(34):6241–64. [3] De-Andrés A, Pérez J, Ortiz M. Elastoplastic finite element analysis of three-dimensional fatigue crack growth in aluminum shafts subjected to axial loading. Int J Solids Struct 1999;36(15):2231–58. [4] Martha LF, Wawrzynek PA, Ingraffea AR. Arbitrary crack representation using solid modeling. Eng Comput 1993;9(2):63–82. [5] Marusich T, Ortiz M. Modelling and simulation of high-speed machining. Int J Numer Meth Eng 1995;38(21):3675–94. [6] Schöllmann M, Fulland M, Richard H. Development of a new software for adaptive crack growth simulations in 3d structures. Eng Fract Mech 2003;70(2):249–68. [7] Børvik T, Hopperstad O, Berstad T, Langseth M. Perforation of 12mm thick steel plates by 20mm diameter projectiles with flat, hemispherical and conical noses: part ii: numerical simulations. Int J Impact Eng 2002;27(1):37–64. [8] Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33(20–22):2899–938. [9] Xu X-P, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42(9):1397–434. [10] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Numer Meth Eng 1999;44(9):1267–82. [11] Zhou F, Molinari J-F. Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency. Int J Numer Meth Eng 2004;59(1):1–24. [12] Dvorkin EN, Cuitiño AM, Gioia G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. Int J Numer Meth Eng 1990;30(3):541–64. [13] Oliver J. On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations. Int J Solids Struct 2000;37(48):7207–29. [14] Mosler J, Meschke G. Embedded crack vs. smeared crack models: a comparison of elementwise discontinuous crack path approaches with emphasis on mesh bias. Comput Methods Appl Mech Eng 2004;193(30):3351–75. [15] Linder C, Armero F. Finite elements with embedded strong discontinuities for the modeling of failure in solids. Int J Numer Meth Eng 2007;72(12):1391–433. [16] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45(5):601–20. [17] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Eng Fract Mech 2002;69(7):813–33. [18] Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. Int J Numer Meth Eng 2003;58(12):1873–905. [19] Melenk JM, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996;139(1–4):289–314. [20] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin methods. Int J Numer Meth Eng 1996;39(6):923–38. [21] Lu Y, Belytschko T, Tabbara M. Element-free Galerkin method for wave propagation and dynamic fracture. Comput Methods Appl Mech Eng 1995;126(1–2):131–53. [22] Belytschko T, Lu Y, Gu L, Tabbara M. Element-free Galerkin methods for static and dynamic fracture. Int J Solids Struct 1995;32(17–18):2547–70. [23] Belytschko T, Lu Y, Gu L. Crack propagation by element-free Galerkin methods. Eng Fract Mech 1995;51(2):295–315. [24] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Meth Eng 2004;61(13):2316–43. [25] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 2007;196(29):2777–99. [26] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Comput Methods Appl Mech Eng 2010;199(37):2437–55. [27] Organ D, Fleming M, Terry T, Belytschko T. Continuous meshless approximations for nonconvex bodies by diffraction and transparency. Comput Mech 1996;18(3):225–35. [28] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48(1):175–209. [29] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83(17–18):1526–35. [30] Silling SA, Bobaru F. Peridynamic modeling of membranes and fibers. Int J Non-Linear Mech 2005;40(2–3):395–409. [31] Ha YD, Bobaru F. Studies of dynamic crack propagation and crack branching with peridynamics. Int J Fract 2010;162(1–2):229–44. [32] Ha YD, Bobaru F. Characteristics of dynamic brittle fracture captured with peridynamics. Eng Fract Mech 2011;78(6):1156–68. [33] Agwai A, Guven I, Madenci E. Predicting crack propagation with peridynamics: a comparative study. Int J Fract 2011;171(1):65–78. [34] Dipasquale D, Zaccariotto M, Galvanetto U. Crack propagation with adaptive grid refinement in 2d peridynamics. Int J Fract 2014;190(1–2):1–22. [35] Wu C, Ren B. A stabilized non-ordinary state-based peridynamics for the nonlocal ductile material failure analysis in metal machining process. Comput Methods Appl Mech Eng 2015;291:197–215. [36] Ni T, Zhu Q-Z, Zhao L-Y, Li P-F. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Eng Fract Mech 2018;188:320–43. [37] Zhou X, Wang Y, Qian Q. Numerical simulation of crack curving and branching in brittle materials under dynamic loads using the extended non-ordinary statebased peridynamics. Eur J Mech A/Solids 2016;60:277–99. [38] Ladd AJ, Kinney JH, Breunig TM. Deformation and failure in cellular materials. Phys Rev E 1997;55(3):3271. [39] Ladd AJ, Kinney JH. Elastic constants of cellular structures. Phys A: Stat Mech Appl 1997;240(1–2):349–60. [40] Caldarelli G, Castellano C, Petri A. Criticality in models for fracture in disordered media. Phys A: Stat Mech Appl 1999;270(1–2):15–20. [41] Lilliu G, van Mier JG. 3d lattice type fracture model for concrete. Eng Fract Mech 2003;70(7–8):927–41. [42] Zhao G-F, Fang J, Zhao J. A 3d distinct lattice spring model for elasticity and dynamic failure. Int J Numer Anal Meth Geomech 2011;35(8):859–85. [43] Chen H, Lin E, Jiao Y, Liu Y. A generalized 2d non-local lattice spring model for fracture simulation. Comput Mech 2014;54(6):1541–58. [44] Monaghan JJ, Gingold RA. Shock simulation by the particle method SPH. J Comput Phys 1983;52(2):374–89. [45] Libersky LD, Petschek AG, Carney TC, Hipp JR, Allahdadi FA. High strain lagrangian hydrodynamics: a three-dimensional SPH code for dynamic material response. J Comput Phys 1993;109:67–75. [46] Maurel B, Combescure A. An sph shell formulation for plasticity and fracture analysis in explicit dynamics. Int J Numer Meth Eng 2008;76(7):949–71. [47] Caleyron F, Combescure A, Faucher V, Potapov S. Dynamic simulation of damage-fracture transition in smoothed particles hydrodynamics shells. Int J Numer Meth Eng 2012;90(6):707–38. [48] Chakraborty S, Shaw A. A pseudo-spring based fracture model for SPH simulation of impact dynamics. Int J Impact Eng 2013;58:84–95. [49] Shaw A, Reid S, Roy D, Chakraborty S. Beyond classical dynamic structural plasticity using mesh-free modelling techniques. Int J Impact Eng 2015;75:268–78. [50] Chakraborty S, Shaw A. Crack propagation in bi-material system via pseudo-spring smoothed particle hydrodynamics. Int J Comput Methods Eng Sci Mech 2014;15(3):294–301. [51] Chakraborty S, Shaw A. Prognosis for ballistic sensitivity of pre-notch in metallic beam through mesh-less computation reflecting material damage. Int J Solids Struct 2015;67:192–204. [52] Chakraborty S, Islam MRI, Shaw A, Ramachandra L, Reid S. A computational framework for modelling impact induced damage in ceramic and ceramic-metal

24

Engineering Fracture Mechanics xxx (xxxx) xxxx

M.R.I. Islam and A. Shaw

composite structures. Compos Struct 2017;164:263–76. [53] Islam MRI, Chakraborty S, Shaw A, Reid S. A computational model for failure of ductile material under impact. Int J Impact Eng 2017;108:334–47. [54] Lahiri SK, Shaw A, Ramachandra L. On performance of different material models in predicting response of ceramics under high velocity impact. Int J Solids Struct 2019;176–177:96–107. [55] Liu G-R, Liu MB. Smoothed particle hydrodynamics: a meshfree particle method. World Scientific; 2003. [56] Liu MB, Liu GR. Restoring particle consistency in smoothed particle hydrodynamics. Appl Numer Math 2006;56(1):19–36. [57] Liu MB, Liu GR. Smoothed particle hydrodynamics SPH: an overview and recent developments. Arch Comput Methods Eng 2010;17(1):25–76. [58] Islam MRI, Chakraborty S, Shaw A. On consistency and energy conservation in smoothed particle hydrodynamics. Int J Numer Meth Eng 2018;116(9):601–32. [59] Chen JK, Beraun JE, Carney TC. A corrective smoothed particle method for boundary value problems in heat conduction. Int J Numer Meth Eng 1999;46(2):231–52. [60] Shaw A, Reid SR. Heuristic acceleration correction algorithm for use in SPH computations in impact mechanics. Comput Methods Appl Mech Eng 2009;198(49):3962–74. [61] Shaw A, Reid SR. Applications of sph with the acceleration correction algorithm in structural impact computations. Curr Sci 2009:1177–86. [62] Shaw A, Roy D. Stabilized sph-based simulations of impact dynamics using acceleration-corrected artificial viscosity. Int J Impact Eng 2012;48:98–106. [63] Monaghan JJ. SPH without a tensile instability. J Comput Phys 2000;159(2):290–311. [64] Monaghan JJ. On the problem of penetration in particle methods. J Comput Phys 1989;82(1):1–15. [65] Ren B, Li S. Meshfree simulations of plugging failures in high-speed impacts. Comput Struct 2010;88(15):909–23. [66] Campbell J, Vignjevic R, Libersky L. A contact algorithm for smoothed particle hydrodynamics. Comput Methods Appl Mech Eng 2000;184(1):49–65. [67] Douillet-Grellier T, Pramanik R, Pan K, Albaiz A, Jones BD, Williams JR. Development of stress boundary conditions in smoothed particle hydrodynamics (sph) for the modeling of solids deformation. Comput Particle Mech 2016:1–21. [68] Song J-H, Wang H, Belytschko T. A comparative study on finite element methods for dynamic fracture. Comput Mech 2008;42(2):239–50. [69] Ramulu M, Kobayashi A. Mechanics of crack curving and branching-a dynamic fracture analysis. Int J Fract 1985;27(3):187–201. [70] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Comput Mech 2007;39(6):743–60. [71] Braun M, Fernández-Sáez J. A new 2d discrete model applied to dynamic crack propagation in brittle materials. Int J Solids Struct 2014;51(21):3787–97. [72] Song J-H, Belytschko T. Cracking node method for dynamic fracture with finite elements. Int J Numer Meth Eng 2009;77(3):360–85. [73] Ravi-Chandar K, Knauss W. An experimental investigation into dynamic fracture: Iii. On steady-state crack propagation and crack branching. Int J Fract 1984;26(2):141–54. [74] Ravi-Chandar K, Yang B. On the role of microcracks in the dynamic fracture of brittle materials. J Mech Phys Solids 1997;45(4):535–63. [75] Bobaru F, Zhang G. Why do cracks branch? a peridynamic investigation of dynamic brittle fracture. Int J Fract 2015;196(1–2):59–98. [76] Rashid M. The arbitrary local mesh replacement method: an alternative to remeshing for crack propagation analysis. Comput Methods Appl Mech Eng 1998;154(1–2):133–50. [77] Tabiei A, Wu J. Development of the dyna3d simulation code with automated fracture procedure for brick elements. Int J Numer Meth Eng 2003;57(14):1979–2006. [78] Kosteski L, D’Ambra RB, Iturrioz I. Crack propagation in elastic solids using the truss-like discrete element method. Int J Fract 2012;174(2):139–61. [79] Ni T, Zhu Q-Z, Zhao L-Y, Li P-F. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Eng Fract Mech 2018;188:320–43. [80] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Eng Fract Mech 2008;75(5):943–60. [81] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Comput Struct 2010;88(23):1391–411. [82] Kalthoff J, Winkler S. Failure mode transition at high rates of shear loading. DGM Informationsgesellschaft mbH, Impact Loading and Dynamic Behavior of Materials 1988;1:185–95. [83] Johnson GR, Cook WH. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In: Proceedings of the 7th international symposium on ballistics, vol. 21, The Hague, The Netherlands, 1983. p. 541–7. [84] Johnson GR, Cook WH. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Eng Fract Mech 1985;21(1):31–48. [85] Teng X, Wierzbicki T, Hiermaier S, Rohr I. Numerical prediction of fracture in the Taylor test. Int J Solids Struct 2005;42(9–10):2929–48. [86] Graham DI, Hughes JP. Accuracy of sph viscous flow models. Int J Numer Methods Fluids 2008;56(8):1261–9.

25