An h-adaptive meshfree-enriched finite element method based on convex approximations for the three-dimensional ductile crack propagation simulation

Journal Pre-proof An h-adaptive meshfree-enriched finite element method based on convex approximations for the three-dimensional ductile crack propagation simulation

Bo Ren, C.T. Wu, Dandan Lyu

PII:

S0167-8396(19)30104-9

DOI:

https://doi.org/10.1016/j.cagd.2019.101795

Reference:

COMAID 101795

To appear in:

Computer Aided Geometric Design

Please cite this article as: Ren, B., et al. An h-adaptive meshfree-enriched finite element method based on convex approximations for the three-dimensional ductile crack propagation simulation. Comput. Aided Geom. Des. (2019), 101795, doi: https://doi.org/10.1016/j.cagd.2019.101795.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Highlights • • • •

An interactive meshfree-enriched finite element formulation is introduced for ductile crack propagation. Convex approximation is utilized to guarantee the compatibility across elements. A crack front tracking and a surface updating algorithm are presented to model crack propagation in three-dimensional space. A meshfree particle integration method is employed for the weak integration.

An h-adaptive meshfree-enriched finite element method based on convex approximations for the three-dimensional ductile crack propagation simulation Bo Ren1, C. T. Wu1*, Dandan Lyu1 1 Computational and Multiscale Mechanics Group Livermore Software Technology Corporation, 7374 Las Positas Road, Livermore, CA 94551, USA

Abstract A new numerical approach for modeling the three-dimensional ductile crack propagation problem is presented. The approach is obtained by the introduction of an h-adaptive finite element refinement scheme using the element-wise meshfree enrichments and meshfree visibility criteria to achieve the desired accuracy due to the material nonlinearity and topological changes in the presence of displacement discontinuity. To ensure the element compatibility condition, meshfree convex approximations are employed to the finite element refinement scheme. A crack propagation algorithm combining the meshfree discretization and the crack surface tracking scheme is proposed. To further simplify the freeform surface generation in this moving discontinuity problem, a meshfree particle integration method is considered. Finally, two numerical benchmarks are studied to examine the effectiveness of present method in the three-dimensional ductile crack propagation analysis. Keywords: Meshfree; finite element method; ductile fracture; convex approximation

1. Introduction The study of automotive component failures is critical to the car structure and safety design. Failure in automotive components is an occurrence which affects the life of almost every person at one stage or another (Heyes 1998). Among various automotive component failures, failure of metallic components such as torsion bars, connection rods, driveshafts, crankshafts and some welded parts can cause crash accidents that leaves passengers to a catastrophic injury from which one may never fully recover. Thus, it is important for automotive industry to use an advanced computational tool that enables for varying metal failure analyses. The finite elements methods incorporating the classical physics of continuum and the linear elasticity theory have been the basis for most computational tools in engineering analyses of brittle fracture. Unlike standard finite element methods which utilize the element deletion strategy to model the surface decohesion and the subsequent crack propagation, the adaptive meshing technology (Bouchard et al. 2000) has been employed. The numerical treatment of crack propagation using the adaptive meshing technology G

∗Corresponding author, C. T. Wu, E-mail: [email protected]  XG G

in three-dimensional problems is not new. Early developments on adaptive remeshing for crack propagation simulation can be traced back to the works of Terzopoulos and Fleischer (1988) and Norton et al. (1991) in computer graphics and animation. Later developments in computer graphics using physically-based animation techniques (O’Brien 2002; Federl and Prusinkiewicz 2004) in the context of continuum mechanics and finite element methods have proven to be convenient to generate special effects for brittle and ductile fracture in film, television and games. Despite the quality of the results produced by those methods is sufficient for computer graphics and animation, they are not suitable for standard crack propagation analysis in the engineering practice. In essence, physically-based animation techniques require the finite element discretization of objects into tetrahedral meshes to compute the internal forces for driving the deformation. Nevertheless, most element cutting algorithms (Sifakis et al. 2007) in discretization impose many restrictions on the geometry of the cuts. As a result, they could easily produce ill-conditioned tetrahedral meshes including jagged or over-smoothed crack surfaces and sliver elements that lead to unphysical crack growth and deteriorate the critical time step in the dynamic fracture simulation. In order to resolve those numerical issues of adaptive remeshing, novel finite element methods such as embedded discontinuity finite element methods (Oliver et al. 2002), extended finite element methods (XFEM) (Moes, Dolbow, Belytschko 1999; Sukumar et al. 2000; Wu et al. 2018b) and generalized finite element methods (GFEM) (Duarte and Kim 2008; Fries and Belytschko 2010) permit a local enrichment of approximation spaces through the partition of unity to model the singularity and displacement discontinuity. To model the crack propagation, XFEM and GFEM have been combined with the level set (Osher and Sethian 1988) or fast marching (Sukumar et al. 2008) methods such that the finite element mesh does not need to be updated or re-meshed upon crack growth. The counterpart of the discrete crack approach is to employ a continuous field variable to represent cracks by the phase-field methods (Miehe et al. 2010). In contrast to discrete descriptions of fracture, phase-field descriptions do not require sophisticated bookkeeping of degrees of freedom or jump conditions in tracking the moving discontinuities. Since phase-field methods utilize a diffuse representation of cracks as a substitute for virtual discontinuities, the amount of crack resolution is controlled by a prescribed length scale. Because of that, phase-field methods usually demand a very fine mesh placed by the accuracy requirement and thus are computationally expensive. An effective alternative for modeling the crack propagation in brittle materials is the bond-based peridynamics proposed by Silling (2000). The bond-based peridynamics is a new continuum theory that describes the material deformation in the brittle fracture process by a nonlocal approach. Different from the weakly or strictly nonlocal models (Bazant and Jirasek 2002; Chen et al. 2000) in the classical continuum mechanics where their spatial derivatives required by the partial differential equations do not exist along the crack surfaces, the peridynamic theory offers the ability to include displacement discontinuities in a continuum body without explicitly modeling the crack surfaces. Bond-based peridynamic models can be implemented within the finite element framework (Ren, Wu and Askari 2017) to improve the accuracy of fracture simulations in a non-uniform discretization as well as to simplify the enforcement of boundary conditions.G Despite those advantages offered by bond-based peridynamic models, the bond failure is linked to fracture energy released rate with a fixed Poisson’s ratio which has limited their applications in nonlinear problems. Those novel numerical methods and their variants share a common feature aiming to minimize the mesh orientation and mesh size sensitivity issues in the material failure analysis. It has been shown (Bourdin et al. 2008; Silling and Lehoucq 2010) that different underlying fracture energy expressions in those methods can be interpreted as a regularization of Griffith’s fracture formulation. Nevertheless, an YG G

extension of those novel finite element methods to the modeling of three-dimensional ductile fracture is not straightforward. The ductile fracture models differ from brittle fracture models in incorporating crack front plastic nonlinear effects for crack growth analysis. As the singular enrichment functions for cracktip displacement fields are no longer applicable for ductile fracture simulation, mesh refinement near the crack front is usually recommended in order to accurately model the nonlinear stress/strain fields and curved geometry. Most of the above reviewed numerical methods are either restricted to the brittle fracture simulation or suffering from numerical difficulties when complex crack topologies are involved in three-dimensional crack propagation problems (Beese et al. 2018; Wu et al. 2018a; Aldakheel et al. 2018). Because the ductile fracture process cannot be described by the means of linear elastic fracture mechanics (LEFM) using the above approaches, developments of ductile fracture criteria along with new freeform surface generation and adaptive mesh refinement schemes for effective crack propagation simulation are needed. In comparison to finite element methods, considerable success of fracture and crack propagation analyses also has been achieved in meshfree methods (Li and Liu 2004). The meshfree fracture analyses based on discrete crack approach can be accomplished by the techniques of modification of kernel (weight) function or the enrichment of discontinuous function. The first technique includes the visibility (Belytschko, Lu, Gu 1994), diffraction (Organ et al. 1996), transparency criteria (Organ et al. 1996) and other methods based on geometric information and different screening effects (Duflot and Nguyen-Dang 2004; Barbiere and Petrinic 2014). The second technique for meshfree fracture analyses incorporates the discontinuous enrichment functions with meshfree approximations either intrinsically or extrinsically (Fleming et al. 1997; Rabczuk and Belytschko 2004; Bordas, Rabczuk, Zi 2008) similar to that in XFEM and GFEM. Among them, the visibility criterion is considered the earliest and easiest meshfree technique to model the propagation of cracks (Krysl and Belytschko 1999). An initial work of applying visibility criterion to meshfree ductile crack propagation analyses was shown by Simonsen and Li (2004). Similar approaches have been developed to the analysis of ductile fracture in thermal-mechanical (Simkins and Li 2006), high-speed (Ren et al. 2011) and structural (Ren and Li 2012) problems. Regardless the ability to avoid complex remeshing operations, those meshfree methods utilize the background cells (Chen et al. 1996) for the domain integration which makes them difficult to model the crack propagation in threedimensional problems. In principle, this demands the development of a pure nodal-integrated numerical scheme based on either strong form or weak form collocation formulations for meshfree methods. On the other hand, meshfree methods based on a direct nodal integration scheme exhibit a poor performance (Chen et al. 2001) known as the spurious energy modes in solid mechanic applications. This unstable mode is caused by rank instability of the meshfree discrete system (Belytschko et al. 2000) and needs stabilization. A pioneering approach to circumvent this numerical instability was demonstrated by Beissel and Belytschko (1996). This stabilization method reconstructs the Galerkin weak form consisting of the residual of equilibrium equation to stabilize the solution. However, the main drawback of this residual stabilization approach is the contradictory demands on the stabilization control parameter due to the accuracy requirement. Several new stabilization methods have been proposed recently that either perform a velocity smoothing based on the linear momentum equation (Pan et al. 2019; Wu et al. 2019) or augment the standard quadratic energy functional by several stabilization terms containing first-order strain gradients (Wu et al. 2017, 2018a) for stabilization. Those novel meshfree methods bypass the need of stabilization control parameters and open new opportunities to model many challenging engineering problems including the material failure and crack propagation simulations (Wu et al. 2016, 2018a, 2019). ZG G

The objective of this study is to introduce a stabilized meshfree method, the Smoothed Particle Galerkin (SPG) method (Wu et al. 2015, 2017, 2018a), to the study of crack propagation problems in the three-dimensional ductile fracture process. Different from our previous work in using the meshfree method for crack propagation simulations (Wu et al. 2016), the present method starts from the finite element method. The finite elements are enriched with meshfree particles via an h-adaptive mesh refinement scheme when certain accumulation of damage is reached. Upon crack initiation, a freeform surface generation algorithm is introduced to the meshfree-enriched finite elements for crack propagation analysis. The reminder of the paper is organized as follows: The field approximation for SPG method is introduced in Section 2. Basic equations of SPG method for the dynamic problem are provided in Section 3. Section 4 describes the h-adaptive mesh refinement scheme. The algorithm for the crack front tracking and surface updating in three-dimensional cases is introduced. Two numerical benchmarks are given in Section 5 to examine the effectiveness of present method. Conclusions are made in Section 6.

2. Overview on GMF approximation In computational mechanics, the application of meshfree-enriched finite element approach was introduced by Wu and Hu (2011a) to reinforce the performance of low-order finite element formulation in the near-incompressible analysis. The meshfree-enriched finite element formulation has been shown to pass the numerical inf-sup tests (Wu and Hu 2011a) and behave very robust in the nonlinear large strain analysis (Wu and Koishi 2012). The central idea of meshfree-enriched finite element formulation is the enrichment of a finite element with a set of meshfree particles using convex approximations (Sukumar 2004; Wang and Chen 2014). The earliest convex approximation in computational mechanics was introduced by Sukumar (2004). Hormann and Sukumar (2017) have demonstrated that the convex approximation generated from maximum entropy coordinates can be related to positive barycentric coordinates with the Lagrangian property. Following the previous works in (Wu and Hu 2011b; Wu and Koishi 2012), the generalized meshfree (GMF) approximation (Wu and Hu 2011b) is used in this study to construct the convex approximation. The foundation of GMF approximation is to revise the Shepard method (Shepard 1968) by an enrichment of basis function to exhibit the linear exactness. The choice of the enriched basis function determines if the GMF approximation has the convexity and boundary Kronecker-delta property. The GMF approximation can reproduce several well-known meshfree approximations. If the enriched basis function is chosen to be the polynomials, it degenerates to the moving least squares (MLS) approximation (Belytschko 1994 et al.; Lancaster and Salkauskas 1981) and the reproducing kernel (RK) approximation (Liu et al. 1995). The GMF approximation will be equivalent to the maximum entropy (ME) (Sukumar 2004; Arroyo and Ortiz 2006) approximation if an exponential basis function is adopted. The GMF approximation can also be considered as a new meshfree approximation if it is constructed by other basis function. The GMF approximation can naturally bear the weak Kronecker-delta property (Sukumar 2004) at boundaries despite it is convex or non-convex. This property makes the imposition of element compatibility in the h-adaptive meshfree-enrichment finite element method and the generation of freeform surfaces for crack propagation problems much easier. The first-order GMF approximation based on exponential basis functions which is equivalent to Sukumar’s maximum entropy approximation is summarized as follows (Wu and Hu 2011b; Wu and Koishi 2012):

[G G

Ψ I ( x, λr ) =

ψI = ψ

φa ( x; X I )Γ I ( X I , λr )

¦ J =1φa (x; X J )Γ J ( X J , λr ) n

for fixed

x

(1)

subjected to n

Rr ( x, λr ) = ¦ Ψ I x I = 0

(2)

ψ I = φa (x; X I )Γ I ( X I , λr ) ,

(3)

ψ = ¦ I =1ψ I = ¦ I =1φa (x; X I )Γ I ( X I , λr ) ,

(4)

I =1

where

n

XI =

n

x − xI is the normalized distance by node I support size a(x) , a(x)

(5)

x is the coordinate of interior point, x I is the coordinate of node I,

φa (x; X I ) is the kernel function of node I with support size a(x) , Γ I ( X I , λr ) is the exponential function of the GMF approximation, n is the number of nodes within the support size a(x) at fixed x , λr (x) are constraint parameters which have to be decided,

r is the number of constraints.

Subsequently, the spatial derivatives of the GMF approximation can be derived and are given by

∇Ψ I = Ψ I ,x + Ψ I ,λr λr ,x , where

(6)

n ψ ψ I ,x − Ψ I ¦ J ,x , ψ J =1 ψ ψ I ,x = φa ,x Γ I + φa Γ I ,x ,

Ψ I ,x =

Ψ I ,λr =

(7) (8)

n φa Γ I ,λ § φa Γ J ,λ · − ΨI ¦¨ ¸, ψ ψ ¹ J =1 © r

(9)

r

λr ,x = − J −1 Rr ,x , n § φa Γ I ,λr J = ¦¨ ψ I =1 ©

n · § φa Γ J ,λr X R ⊗ − ⊗ ¸ ¨ ¦ I r ψ J =1 © ¹

n

n

I =1

I =1

Rr ,x = ¦ Ψ I ,x X I + ¦ Ψ I X I ,x .

(10)

· ¸, ¹

(11) (12)

It is worthwhile to note that the property of the partition of unity in the above GMF approximation is automatically satisfied by the definition of normalization form in Eq. (1). The completion of the GMF approximation is achieved by finding λ to satisfy Eq. (2). Eq. (2) is a set of constraint equations which ensure the linear reproducing condition in GMF approximation. To determine λ at any fixed x in Eq. \G G

(1), a root-finding algorithm is required for the non-linear basis functions.

3. Basic equations of SPG Method 3.1 The SPG weak formulation This section provides an overview on the SPG method which has been utilized to analyze the dynamics crack propagation problems in this study. The SPG weak form is formulated based on the stabilization method (Wu et al. 2017, 2018a). The admissible space for the displacement fields in SPG method is defined by

{

H g1 ( Ω) = v : v Ω∈ H1 ( Ω) , v = vg on Γg

}

(13)

where v g is the applied displacements on Dirichlet boundary Γg . The total particles in discretization are denoted by an index set Z I = {x I }INP=1 ⊂ ℜ3 . The displacement field using the convex approximation from Eq. (1) is defined by NP

u h ( x ,t ) = ¦ Ψ I ( x ) u ( x I , t ) ≡uˆ ( x ,t ) ∀x ∈ Ω GGGGGGGGGGGGGGGGGGGGGGGGGGGGG(14) I =1

where NP is the total number of particles in discretization of domain Ω . ΨI ( x) , I=1,…, NP are defined in Eq. (1) which can be considered as the Lagrangian shape functions of meshfree approximation for the h

displacement field u . Now the weak form of the dynamic problem is written as (Wu et al. 2017, 2018a)

³

Ω

T δ uˆ ⋅ ρ 0 uˆ d Ω = ³ δ F : P ( F ) d Ω + ³ δ F T : P ( F )d Ω − l ext ( uˆ ) , Ω Ω   

  

standard term

{

∀δ uˆ ∈ H = v : v 1 0

stabilization term

∈ H ( Ω ) , v = 0 on Γ g 1

Ω

}

(15)

where ȡ0 is the initial material density. The variation of deformation gradient F and enhanced deformation gradient F for stabilization are defined as:

§ ∂x · ∂δ uˆ δF =δ ¨ ¸ = © ∂x ¹ ∂x

δ F = ∇ (δ F ) ⋅ β ( x )

∀x ∈ Ω

 ( x; x - ȟ )( ȟ - x ) dȍ β ( x) = ³ Ȍ Ω

(16) (17) (18)

where x is the coordinates in the current configuration, and x is the coordinates in the reference configuration.

 in Ȍ

Eq. (18) is the strain smoothing function (Wu et al. 2017, 2018a) introduced for

]G G

stabilization, and β denotes the stabilization coefficient matrices. The external force is expressed in a standard way given by:

l ext ( uˆ ) = ³ δ uˆ ⋅ bd Ω + ³ δ uˆ ⋅ hd Γ h Ω

(19)

Γh

where b is the body force, and h is the surface traction applied on the Neumann boundary Γ h . The nominal stress P can be related to the Cauchy stress ı in metal plasticity by

P = J 0 F −1 ⋅ ı

(20)

and J0 is the determinant of the deformation gradient F. The second term on the right-hand side (RHS) of

~

Eq. (15) represents the stabilization functional which contains an enhanced stress field P [see (Wu et al. 2017, 2018a) for details] to provide the stabilization effect. G 3.2 Semi-discrete formulations Using the first-order meshfree convex approximation (Wu et al. 2011b) for Ψ ( x ) and the Shepard

 ( x ) in stabilization (Wu et al. 2017, 2018a) leads to the following discrete form of the function for Ȍ momentum equation to be solved in explicit dynamics analysis:

~  M lumpUˆ = f ext − f int − f stab

(21)

NP

M Ilump = ¦ ρ 0 Ψ I ( x N ) VN0 I [3×3]

(22)

N =1

where

ext ) is the standard external force matrix, Uˆ = §¨ Uˆ1 ,",UˆNP ·¸ is the matrix f ext = ( f1ext ,", f NP T

T

©

containing nodal accelerations. V

0 N

is the initial nodal volume of node N, and M

¹

lump I

is the lumped

nodal mass matrix of node I. The computation of internal force vectors involves variation of the deformation gradient tensor. The discrete form of Eq. (16) using the convex approximation gives

∂Ψ I δ uiI I =1 ∂x j NP

δ Fij = ¦

(23)

Using the direct integration (DNI) scheme, the internal force and stabilization force are computed by: DNI NP

f Iint =

¦ B ( x ) P ( x )V T I

N

N

0 N

∀x N ∈ Ω

(24)

N =1

DNI NP

fIstab =

¦ B ( x ) P ( x )V T

I

N

N =1

^G G

N

0 N

(25)

~

The components of strain-gradient matrices B I and B I can be found in (Wu et al. 2017, 2018a), thus they are omitted in this paper. The nominal stress P defined with respect to the reference configuration also can be expressed by the Cauchy stress ı for convenience in metal plasticity computations. Using the Voigt rule, we obtain DNI NP

f Iint =

¦ B ( x ) P ( x )V T I

N

N

0 N

N =1

NP

= ¦ BIT ( x N )ı ( x N ) J 0VN0

(26)

N =1

Likewise, we have DNI NP

fIstab =

¦ B ( x ) P ( x )V T I

N

N

N =1

0 N

NP

= ¦ B IT ( x N ) ı ( x N ) J 0VN0

(27)

N =1

~ where the first-order strain-gradient matrices BI and B I can be found in (Wu et al. 2017, 2018a) as well. Finally, the semi-discrete equations in Eq. (21) can be integrated by the Newmark method in the temporal space. The critical time step for the central difference time integration is governed by the CourantFriedrichs-Lewy (CFL) condition for explicit dynamics analyses and is computed following the previous works in (Park et al. 2011; Wu et al. 2017, 2018a).

4. The h-adaptive meshfree-enriched finite element approach for crack propagation The h-adaptivity is a key ingredient in many non-linear solid and structural simulations involving high gradients and moving discontinuity solutions. The purpose of introducing the h-adaptivity in this study is to enhance the local accuracy in the coarse and low-order finite elements due to the presence of plasticity and strong discontinuity in which the high accuracy usually is difficult to achieve with fixed number of degrees of freedom. 4.1 The h-adaptive scheme In the proposed h-adaptive scheme, the coarse element will firstly be enriched by a group of SPG particles (called the SPG group in this paper) when the damage indicator D reaches a certain value. This value is a numerical parameter that is smaller than the damage threshold for crack propagation. In this study, the damage indicator is defined as an averaged value in the FEM element and is given by: Ng

D=¦ I =1

DI Ng

(28)

Here DI is the internal damage variable at Gaussian points. DI can be a damage valuable from continuum damage mechanics or simply defined as the effective plastic strain. ܰ௚ =8 is the total number of Gaussian points in the tri-linear FEM element. Once the element-wise D meets the adaptivity criterion, SPG particles are enriched into that element, and the finite element formulation is switched to SPG formulation for the crack initiation and propagation analysis. The position of an enriched SPG particle in the SPG group is calculated using standard finite element shape functions by: _G G

N

x p ,i = ¦ x I ,i Ψ FE (ξ ,i ), i = 1, 2, 3

(29)

I =1

where Ψ

FE

is the regular FEM shape function. N=8 denotes the total number of nodes in the FEM

element. x I ,i , i = 1, 2,3 are the FEM nodal coordinates. ߦ is the local coordinate in isoparametric space. Here ߦ is equally distributed along each dimension that leads to uniformly distributed SPG particles (2*2*2, 3*3*3, or 4*4*4 …in each element) inside the FEM element space as shown in Fig. 1 where the red dots denote SPG particles. The total mass and volume of that FEM element are equally distributed to all SPG particles to keep the consistence of mass. The kinematic field valuables like displacements and velocities of SPG particles also can be interpolated by Eq. (29). Other internal valuables such as stress, strain, effective plastic strain and damage value are not computed at finite element nodal points, and therefore cannot be interpolated directly at SPG particles using Eq. (29). To use Eq. (29) for interpolation, the internal valuables at each finite element nodal point are assigned to the values from the closest Gaussian point associated with that nodal point. Using this assignment, the internal valuables at each enriched SPG particle can be approximated by N

Θ p = ¦ Θ I Ψ FE ( x p )

(30)

I =1

where Θ p and ΘI are internal valuables of SPG particle and FEM Gaussian point, respectively.

Fig. 1 The sketch of an h-adaptive meshfree-enriched finite element where 4*4*4 SPG particles are used to enrich the element 4.2 Crack-front tracking and surface updating algorithm The goal of present h-adaptive meshfree-enriched finite element method is to effectively trace the crack front in a three-dimensional ductile fracture problem. To reach this goal, a crack surface that represents the strong discontinuity in the three-dimensional space has to be defined. In this paper, the crack surfaces are always defined in the reference configuration. 4.2.1 Crack initiation `G G

The crack initiates when the accumulation of damage reaches the critical value. Theoretically the location of crack initiation is a natural outcome of damage evolution. In other words, the crack initiation could start anywhere inside the computational domain. However, in practice, cracks usually initiate from the domain boundary. Therefore, we assume that cracks only can initiate from the enriched elements on the solid surface in this study. Let’s recall that each SPG group is defined as a set of SPG particles within each enriched FEM element in the h-adaptive procedure. For each enriched FEM element on the solid surface, a damage center x d is determined in the reference configuration defined by (Wu et al. 2016):

( )

NS

x I exp DI

I =1

¦ exp ( DJ )

xd = ¦

(31)

J

where NS = 4 × 4 × 4 is the total number of SPG particles per meshfree-enriched finite element used in this study. On the other hand, it is known that material damage in ductile fracture is a nonlocal phenomenon. Since the internal damage indicator DI defined in Eq. (28) is a local value, it can only be used as the indicator for h-adaptivity and cannot be used for crack initiation. To represent the nonlocality in predicting the ductile crack initiation, an integral type of non-local damage indicator D (Chen et al. 2000; Wu et al. 2016) is defined in this following for the crack initiation:

D ( xd ) = ³ e Ψ c ( xd ) D ( x )d Ωe

(32)

Ω

e

where Ψ c is shape function for regularization. Ω denotes an element-wise regularization process zone c

for the computation of non-local damage indicator D . In general, the size of Ψ is much smaller than the element size. Once the non-local damage variable D in the meshfree-enriched finite element reaches the threshold for crack initiation, the material is considered fully degraded. The transition of displacement field from continuous damage to discrete crack will then be triggered. This initiated crack surface will be formed numerically to pass the damage center x d . To completely define the crack surface, the surface normal is needed. In this work, the criterion of maximum principle stress is utilized to define the crack surface normal n which means only Mode I cracks are monitored in the crack propagation process. 4.2.2 Crack surface discretization Upon crack initiation, the crack surface cuts each meshfree-enriched finite element into two parts at one explicit time step. Based on this cutting procedure, the geometric relation between each edge of FEM element and crack surface needs to be tested. As shown in Fig. 2, let an edge segment of FEM element, ሬሬሬሬሬറ ൌ ࢞௡ െ ࢞ௗ and crack surface normal n have ሬሬሬሬሬറ, be cut by a crack surface at point x , then the vector݊݀  n

the following relation: ሬሬሬሬሬറ ݊݀ ή ࢔ ൌ Ͳ XWG G

(33)

ሬሬሬሬሬറ as follows: Meanwhile, x n can be defined by a factor ߣ and a vector ‫ܤܣ‬

JJJG x n = x A + λ AB

.

(34)

ሬሬሬሬሬറ will intersect the Substituting Eq. (34) to Eq. (33), ߣ can be solved. If െͳ ൑ ߣ ൑ ͳ, the line segment  crack surface. Then all intersection points and the damage center ࢞ௗ can be used to construct the triangularized crack surface mesh in 3D space as shown in Fig 3.

Fig. 2 The interaction between element segment and crack surface

Fig. 3 The initial crack surface in one FEM element The crack surface mesh represents a physical crack existing in the three-dimensional domain. To introduce this discontinuity in the displacement field, the connectivity between two SPG particles locating on different sides of the crack surface mesh must be cut off. In this paper, an algorithm that introduces the 3D parametric visibility condition is adopted to introduce the discontinuity. The detail of this algorithm is referred to the work by Ren and Li (2012). 4.2.3 Crack propagation

XXG G

Fig. 4 Crack surface mesh and crack front line in 3D solid For an evolving crack surface, the crack front lining with many crack front segments, shown as red lines in Fig. 4, is maintained when the crack surface mesh is updated. During the calculation, the crack front condition is monitored to determine the propagation of cracks. For each crack front segment, the following equation is used to calculate the regularized damage indicator.

D ( x c ) = ³ Ψ c ( x c ) D ( x )d Ω

(35)

Ωx

c

As shown in Fig. 5, xc is the center point of a front line segment. Ψ is regularization shape function defined in Eq. (32). The searching zone Ω x is a cone shape zone with the viewing point xc . The utilization of cone shape instead of spherical shape in detecting the next damage center for crack propagation is to eliminate the non-physical crack reversing phenomena. The angle of this cone shape 0 zone is chosen to be 150 in this study. Once the damage indicator D ( x c ) meets the crack propagation

criterion, the position of damage center ( x d ) for each meshfree-enriched finite element inside the cone shape searching zone can be obtained by Eq. (31). Subsequently, the two vectors, one comes from xc to

x d and another goes along with this front line segment (shown as ࢂ૚ and ࢂ૛ in Fig. 5), define the crack propagation surface. The defined crack surface will cut the SPG group located before this crack front segment into two parts. The same algorithm repeats to update the crack surface mesh and to define the meshfree connectivity by the newly generated surface triangles.

XYG G

Fig. 5 Illustration of searching zone in the crack propagation algorithm

5. Numerical examples Two numerical examples are studied to examine the effectiveness of the present method for crack propagation analysis. In both examples, the elasto-plastic material model is used and the effective plastic strain is adopted as the damage indicator. The computer simulations are carried out on Lenovo T40 laptop with i7 core. It is noted that since the crack surface will cut through the meshfree-enriched finite element in one time step, increasing the total number of SPG particles within the element will have no significant change in the computational results. 5.1. Two crack surfaces propagation in a ductile solid In the first example, a steel plate with two off-set notches is stretched at two ends with a constant velocity of 5.0 m/s as shown in Fig. 6. The plate has a dimension of 135m × 200m × 15m. The problem is initially discretized by 1564 hexahedral elements. The Young’s modulus of the steel plate is 190GPa and a linear hardening modulus is 100 MPa. The Poisson’s ratio is 0.3. The purpose of this example is to demonstrate the effectiveness of the present method in modeling propagating and merging cracks. The value of effective plastic strain for h-adaptivity in Eq. (28) is set as 0.35%. The critical effective plastic strain for crack initiation and crack propagation is 0.4%.

XZG G

Fig. 6 The configuration of the steel plate The result of dynamic crack propagation is shown in Fig. 7. It is shown that the effective plastic strain accumulates near two notches initially. Once the damage indicator meets the adaptive criterion, FEM elements are replaced by a number of SPG groups as shown in Fig. 7 (a). The effective plastic strain in SPG groups keeps increasing as steel plate deforms. When the crack initiation criterion is satisfied, the crack surface is nucleated as shown in Fig. 7(b). The crack propagates with the increasing effective plastic strain in front of crack front as shown in Fig. 7(c). Because of the specific geometrical configuration in this example, a high plastic strain zone was built up between two crack fronts as shown in Fig. 7(d). This leads to a change of crack propagation directions for two cracks as shown in Fig. 7(e). Finally, these two cracks merge to form a major crack connecting by a closed-loop crack as shown in Fig. 7(f).

(a) t=38 ms

(b) t=42 ms

X[G G

(c) t=54 ms

(d) t=61ms

(e) t=66 ms

(f) t= 68 ms

Fig. 7 The dynamic crack propagation history with the effective plastic strain contour Fig. 8 illustrates the crack surfaces constructed by the present crack-front tracking and surface updating algorithm. No crack is locked from propagation. The result suggests the present method is able to produce the stable crack propagation in the three-dimensional ductile fracture simulation.

Fig. 8 The illustration of crack surfaces 5.2. The failure process of two jointed bar and convergence study In this numerical example, two rectangular bars are joined together at the location inside the red square in Fig. 9. Each bar has a 0.5m × 0.5m cross section. This example is used to mimic the ductile X\G G

fracture of welded parts in automotive applications. The top and bottom bars are stretched along the x direction at a constant velocity of 5.0 m/s. Since the imposed velocity boundaries are offset in the y direction, they introduce material failure outside the joint zone. Unlike the first numerical example, there is no preset crack existing in this example. In other words, cracks can be initiated in any surface element of the jointed bars. The numerical model is originally discretized by 1920 hexahedral elements. For simplicity, the material is assumed to be isotropic and homogenous. The Young’s modulus of the material is 190GPa with a linear hardening modulus of 200MPa. The Poisson’s ratio is 0.3. The effective plastic strain for hadaptive is 0.35%. The critical effective plastic strains for crack initiation and crack propagation are 0.4%.

Fig. 9 The configuration and discretization of two joined bars (length unit: m)

The final deformed shape with the crack pattern is shown in Fig. 10 where two high effective plastic zones can be found at top and bottom bars close to the joint position. Even though the h-adaptive refinement is triggered in both bars, the crack is only observed at low bar because of the non-symmetric geometric configuration.

Fig. 10 The final deformation (effective plastic strain) of two joined bars Fig. 11 reports the crack propagation history. As shown in Figs. 11(a) and (b), the crack is initialized at the inner surface of the bottom bar. At the same time, the h-adaptive refinement is also triggered at the upper bar but without any crack initiation. Then the crack propagates very fast in the bottom bar as shown from Fig. 11(b) to Fig. 11(e). Fig. 11(f) shows the final crack propagation pattern. X]G G

In the following, three numerical simulations with different discretization (coarse mesh: 1920 hexahedral elements, medium mesh: 3750 hexahedral elements and fine mesh: 5940 hexahedral elements) are carried out to study the discretization effect of the proposed method. The results of reaction force response are presented in Fig. 12. Three force curves match well before cracks initiate. However, there is a noticeable offset of force curves between the coarse mesh and two others when the crack starts to initiate. As the mesh is continuously refined, the reaction force curve converges as shown in Fig. 12. This observation indicates the element-wise regulation process zone defined in Eq. (32) still produces certain mesh-dependent results in the coarse mesh. Therefore, a full regularization process zone across neighbor elements will need to be considered in the future. On the other hand, three discretization models predict same slope of force curves in the crack propagation phase. This suggests crack propagation speeds are approximately the same in three discretization models. Fig. 13 compares the final crack patterns in coarse and medium discretization models. The comparable results are obtained which suggest the effectiveness of the present method for the three-dimensional ductile crack propagation simulation.

(a) t=6.3 ms

(b) t=6.7 ms

(c) t=7.0 ms

(d) t=7.2 ms

(e) t=7.3 ms

(f) t=10 ms X^G

G

Fig. 11 The dynamic crack propagation history with the effective plastic strain contour

Fig. 12 Comparison of reaction force curves with different meshes

(a) 1920 elements

GG

(b) 3750 elements

G

Fig. 13 Comparison of final crack patterns in two discretization models

6. Conclusions In this paper, we have presented a novel numerical approach to model the three-dimensional ductile crack propagation problem. It adopts a local meshfree discretization and incorporates visibility criteria for efficiently updating shape functions upon topological changes when the crack is present. Motived by the meshfree-enriched finite element method for near-incompressible analysis, the element-wise meshfree enrichments based on meshfree convex approximations has been adopted in this study to ensure the X_G G

element compatibility condition across element edges. The crack-front tracking and surface updating algorithm using the triangle mesh has been developed to propagate the cracks. The system of equations is integrated using the particle-based SPG method to explore the possibility in the three-dimensional ductile crack propagation analysis. The numerical results in this study reveal the robustness and stability of present approach in the threedimensional ductile crack propagation analysis. In particular, the cracks can be propagated smoothly without smearing or experiencing numerical locking from propagation. Enhancements of present method including a full regularization process zone for Eq. (32) will be needed. Further study for more complicated crack patterns and geometries in automotive applications will be considered in the near future.

Acknowledgements The authors wish to thank Dr. John O. Hallquist of LSTC for his support to this research. The support from Honda R&D Co., Ltd., Automobile R&D Center, Japan and JSOL Corporation, Japan is also gratefully acknowledged.

References Aldakheel F., Wrigger P., Miehe C. 2018. A modified Gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling. Comput. Mech. 62, 815-833. Arroyo M., Ortiz M. 2006. Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods. Int. J Numer. Methods Eng. 65, 2167-2202. Barbiere E., Petrinic N. 2014. Three-dimensional crack propagation with distance-based discontinuous kernels in meshfree method. Comput. Mech. 53, 325–342. Bazant Z., Jirasek M. 2002. Nonlocal integration formulations of plasticity and damage: survey and progress. J Eng. Mech. 128, 1119-1149. Beese S., Loehnert S., Wrigger P. 2018. 3D ductile crack propagation within a polycrystalline microstructure using XFEM. Comput. Mech. 61, 71-88. Belytschko T., Lu Y.Y., Gu L. 1994. Element-free Galerkin methods. Int. J Numer. Method Eng. 37, 229-256. Belytschko T., Guo Y., Liu W.K., Xiao S.P. 2000. A unified stability analysis of meshless particle methods. Int. J Numer. Methods Eng. 48, 1359–1400.͒ Beissel S., Belytschko T. 1996. Nodal integration of the element-free Galerkin method. Comput. Methods Appl. Mech. Eng. 139. 49-74. Bordas S., Rabczuk T., Zi G. 2008. Three-dimension crack initiation, propagation, branching and junction in nonlinear materials by an extended meshfree method without asymptotic enrichment. Eng. Fract. Mech. 75, 943–960. Bouchard P.O., Bay F., Chastel Y. 2000. Crack propagation modelling using an advanced remeshing technique. Comput. Methods Appl. Mech. Eng. 189, 723-742. Bourdin B., Francfort G., Marigo J.J. 2008. The variational approach to fracture. J Elastic. 91, 1-3. Chen J.S., Pan C., Wu C.T., Liu W.K. 1996. Reproducing kernel particle methods for large deformation analysis of non-linear structures. Comput. Methods Appl. Mech. Eng. 139, 195–227. ͒ Chen J.S., Wu C.T., Belytschko T. 2000. Regularization of material instabilities by meshfree approximations with intrinsic length scales. Int. J Numer. Methods Eng. 47, 1303-1322. Chen J.S., Wu C.T., Yoon S., You Y. 2001. A stabilized conforming nodal integration for Galerkin meshfree methods. Int. J Numer. Methods Eng. 50, 435–466. Duarte C.A., Kim, D.J. 2008. Analysis and applications of a generalized finite element method with global-local enrichment. Comput. Methods Appl. Mech. Eng. 197, 487-504. Duflot M., Nguyen-Dang H. 2004. A meshless method with enriched weight functions for fatigue crack growth. Int. J Numer. Methods Eng. 59,1945–1961. Fleming M., Chu Y.A., Moran B., Belytschko T. 1997. Enriched element-free Galerkin methods for crack tip fields.

X`G G

Int. J Numer. Methods Eng. 40,1483–11504. Heyes, A.M. 1998. Automotive component failures. Eng. Fail. Anal. 5, 129-144. Hormann K., Sukumar N. 2017. Generalized Barycentric Coordinates in Computer Graphics and Computational Mechanics. CRC Press, Boca Raton, FL. Federl P, Prusinkiewicz P. 2004. Finite element model of fracture formulation on growing surfaces. Lect. Notes Comput. Sc. 3037, 138-145. Fries, T.P., Belytschko T. 2010. The extended/generalized finite element method: an overview of the method and its applications. Int. J Numer. Methods Eng. 84, 253-304. Krysl P., Belytschko T. 1999. The element-free Galerkin method for dynamic propagation of arbitrary 3D cracks. Int. J Numer. Methods Eng. 44, 767–800. Lancaster P., Salkauskas K. 1981. Surfaces generated by moving least squares methods. Math. Comput. 37, 141-158. Li S., Liu W.K. 2004. Meshfree Particle Method. Springer Verlag, Berlin. Liu W.K., Jun S., Zhang Y.F., 1995. Reproducing kernel particle methods. Int. J Numer. Methods Fl. 20, 1081-1106. Miehe C., Hofacker M., Welschinger F. 2010. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput. Methods Appl. Mech. Eng. 199, 2765-2778. Moes N., Dolbow J., Belytschko T. 1999. A finite element method for crack growth without remeshing. Int. J Num. Methods Eng. 46, 131-150. Norton A., Turk G., Bacon B., Gerth J., Sweeney P. 1991. Animation of fracture by physical modeling. Vis. Comp. 7, 201-219. O’Brien J.F., Bargteil A.W., Hodgins J.K. 2002. Graphical modeling and animation of ductile fracture. ACM Trans. Graph. 21, 291-294. Oliver J., Huespe A., Pulido M.D.G., Chaves E. 2002. From continuum mechanics to fracture mechanics: the strong discontinuity approach. Eng. Fract. Mech. 69, 113-136. Organ D., Fleming M., Terry T., Belytschko T. 1996. Continuous meshless approximations for nonconvex bodies by diffraction and ͒transparency. Comput. Mech. 18, 1–11. ͒ Osher S., Sethian J.A. 1988. Fronts propagation with curvature-dependent speed: Algorithms based on HamiltonJacobi formulations. J Comput. Phys. 79, 12-49. Pan X., Wu C.T., Hu W., Wu Y.C. 2019. A momentum-consistent stabilization algorithm for Lagrangian particle methods in the thermo-mechanical friction drilling analysis. Comput. Mech. DOI: https://doi.org/10.1007/s00466019-01673-8. Park C.K., Wu C.T., Kan C.D. 2011. On the analysis of dispersion property and stable time step in meshfree method using the generalized meshfree approximation. Finite Elem. Anal. Des. 47, 683-697. Rabczuk T., Belytschko T. 2004. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int. J Numer. Methods Eng. 61, 2316–2343. Ren B., Li S., Qian J., Zeng X. 2011. Meshfree simulations of spall fracture. Comput. Methods Appl. Mech. Eng. 200, 797–811. Ren B., Li S. 2012. Modeling and simulation of large-scale ductile fracture in plates and shells. Int. J Numer. Methods Eng. 49, 2373-2393. ͒ Ren B., Wu C.T., Askari E. 2017. A 3D discontinuous Galerkin finite element method with the bond-based peridynamics model for dynamic brittle failure analysis. Int. J Impact Eng. 99, 14-25. Shepard D. 1968. A two-dimensional interpolation function for irregularly-spaced data. ACM Proceedings of the 1968 23rd ACM National Conference, New York, 517-524, DOI: 10.1145/800186.810616. Sifakis E., Der K.G., Fedkiw, R. 2007. Arbitrary cutting of deformable tetrahedralized objects. ACM Proceedings of the Symposium on Computer Animation’07, 73-80. Silling S.A. 2000. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech. Phys. Solids 48, 175-209. Silling S.A., Lehoucq R. 2010. Peridynamic theory of solid mechanics. Adv. Appl. Mech. 44, 73-168. Simkins D.C., Li S. 2006. Meshfree simulations of thermos-mechanical ductile fracture. Comput. Mech. 38, 235–249. ͒ Sukumar N., Moes N., Moran B., Belytschko, T. 2000. Extended finite element method for three-dimensional crack modelling. Int. J Numer. Methods Eng. 48, 1549-1570. Sukumar N., 2004. Construction of polygonal interpolants: a maximum entropy approach. Int. J Numer. Methods Eng. 61, 2159-2181. Sukumar N., Chopp D., Bechet E., Moes N. 2008. Three-dimensional non-planar crack growth by a coupled extended finite element and fast marching method. Int. J Numer. Eng. 76, 727-748.

YWG G

Terzopoulos D., Fleischer K. 1988. Modeling inelastic deformation: viscoelasticity, plasticity, fracture. Comput. Graph. 22, 269-278. Wang D., Chen P. 2014. Quasi-convex reproducing kernel meshfree method. Comput. Mech. 54, 689-709. Wu C.T., Hu W. 2011a. Meshfree-enriched simplex elements with strain smoothing for the finite element analysis of compressible and nearly incompressible solids. Comput. Methods Appl. Mech. Eng. 200, 2991-3010. Wu C.T., Park C.K., Chen J.S. 2011b. A generalized approximation for the meshfree analysis of solids. Int. J Numer. Methods Eng. 85, 693-722. Wu C.T., Koishi M. 2012. Three-dimensional meshfree-enriched finite element formulation for micromechanical hyperelastic modeling of particulate rubber composites. Int. J Numer. Methods Eng. 91, 1137-1157. Wu C.T., Koishi M., Hu W. 2015. A displacement smoothing induced strain gradient stabilization for the meshfree Galerkin nodal integration method. Comp. Mech. 56, 19-37. Wu C.T., Ma N., Takada K., Okada H. 2016. A meshfree continuous-discontinuous approach for the ductile fracture modeling in explicit dynamics analysis. Comput. Mech. 58, 391-409. Wu C.T., Wu Y.C., Crawford J.E., Magallanes J.M. 2017. Three-dimensional concrete impact and penetration simulations using the smoothed particle Galerkin method. Int. J Impact Eng. 106, 1-17. Wu C.T., Bui T.Q., Wu Y.C., Luo T.L., Wang M., Liao C.C., Chen P.Y., Lai Y.S. 2018a. Numerical and Experimental validation of a Particle Gakerkin Method for Metal Grinding Simulation. Comput. Mech. 61, 365-383. Wu C.T., Ma N., Guo Y., Hu W., Takada K., Okada H., Saito K. 2018b. A dynamic ductile failure analysis of shell structures using a nonlocal XFEM method with experimental validation. Adv. Eng. Softw. 123, 1-12. Wu C.T., Wu Y.C., Lyu D.D., Pan X., Hu W. 2019. The momentum-consistent smoothed particle Galerkin (MCSPG) method for the extreme thread forming simulation in flow screwing process. Comput. Part. Mech. DOI: 10.1007/s40571-019-00235-2.

YXG G