Extended finite element modeling of 3D dynamic crack growth under impact loading

Finite Elements in Analysis and Design 151 (2018) 1–17

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Extended finite element modeling of 3D dynamic crack growth under impact loading Thomas Elguedj a, ∗ , Romains Pelée de Saint Maurice a,b , Alain Combescure a , Vincent Faucher c , Benoit Prabel b a

Univ Lyon, INSA-Lyon, CNRS UMR5259, LaMCoS, F-69621, France CEA, DEN, DANS, DM2S, SEMT, DYN, Gif-sur-Yvette, F-91191, France c CEA, DEN, CAD, DTN, Dir, Saint-Paul-les-Durance, F-13108, France b

A R T I C L E

I N F O

Keywords: X-FEM Dynamic crack growth Explicit dynamics 3D crack growth

A B S T R A C T

This paper is devoted to the numerical simulation of the dynamic propagation of non-planar 3D cracks under transient loadings. For that purpose, in an explicit finite element code, the X-FEM with only discontinuous enrichment is used coupled with level sets to represent the crack geometry. We show on a complex example that the commonly used Hamilton-Jacobi based level set updating algorithm lacks robustness and numerical efficiency to model such types of problems. To circumvent this, we introduce a simple geometric updating algorithm that does not involve any equation solving, and is therefore numerically cheap. The robustness of this updating algorithm is demonstrated on several virtual propagation examples. Finally common benchmark used in the literature as well as a complex experiment involving the dynamic propagation of a non-planar 3D crack under impact loading are used to demonstrate the efficiency of the approach.

1. Introduction The numerical modeling of crack propagation has drawn a lot of attention in the scientific community in the past 25 years. New numerical methods, such has Meshfree and eXtended/Generalized Finite Element Methods opened new possibilities as the discretization does not have to fit the a priori unknown crack path. The X-FEM, first introduced by Black and Belytschko [1], Moës et al. [2], and the GFEM first introduced by Strouboulis et al. [3] is based on a local partition of unity following the work of Babuška and Melenk [4]. Enrichment or handbook functions are introduced in the approximation spaces in accordance with analytical or particular solutions that improve the approximation properties of the numerical method. For cracks modeling, two types of enrichment are commonly used. Crack front enrichment based on Westergaard’s asymptotic solutions enables to capture the singularity, and discontinuous enrichment functions represent the discontinuity away from the crack tip. Several 2D cases have been studied in the literature for various problems, such as static and fatigue crack growth [2,5], fracture in porous medium [6–8], crack branching [9], fracture with contact and friction [10–12], etc. The extension to 3D cases needs specific tools to represent the crack face and front, and level sets have

been used by many authors as it allows an implicit representation of the crack [13–21]. For the case of dynamic crack propagation, in particular when transient loadings are involved, the literature is mostly restricted to 2D cases [18,22–37]. In contrast there has been a reduced number of work focusing on dynamic 3D crack growth modeling. Krysl and Belytschko [38] introduced an Element-free Galerkin method coupled with explicit time integration, but solved only quasi-static cases. Duarte et al. [39] introduced a 3D GFEM for dynamic crack growth but studied essentially 2D examples under static and transient loadings with only one element through the thickness and one 3D example with a transient loading. Ruiz et al. [40] studied 3D bending tests with dynamic loadings but only considered planar propagation, Rabczuk et al. [41,42] studied various examples but those either had a static load or a planar crack propagation, likewise for Duan et al. [43]. We are interested in this work to model the initiation and propagation of non-planar 3D cracks under impact loadings using the X-FEM. As the characteristic time of the phenomenon involved are small, we use an explicit dynamics code that is based on the central difference scheme. The use of X-FEM for crack modeling with explicit time integrators has already been studied, and the present work uses some of the results shown in Refs. [24,25,27,29,44]. The modeling of non-planar 3D

∗ Corresponding author. E-mail address: [email protected] (T. Elguedj). https://doi.org/10.1016/j.finel.2018.08.001 Received 21 September 2017; Received in revised form 3 August 2018; Accepted 6 August 2018 Available online XXX 0168-874X/© 2018 Elsevier B.V. All rights reserved.

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

cracks, coupled with X-FEM/GFEM can be done with several methods. Implicit representation of the crack using level sets is generally used in conjunction with X-FEM [13–16,45]. Explicit representations using non-planar 2D meshes have also been used with X-FEM [46] and GFEM [39,47,48]. Hybrid representations, also called implicit-explicit methods, have also been introduced so as to use the advantages of both approaches [20,21]. In the present work, we have chosen to use level sets and therefore an implicit representation of the crack surface and front using two signed distance functions. In particular, we show that standard level set evolution algorithms using Hamilton-Jacobi equations are not robust enough when applied to crack modeling to ensure an efficient representation of complex propagating cracks and can be costly from a numerical point of view. Previously cited implicit-explicit approaches are one way of circumventing the problem, as well as the use of approximate element local level sets as in Duan et al. [43]. Geometrical approaches were first introduced to circumvent this issues by Ventura et al. [49] in 2D and generalized to 3D by Agathos et al. [50]. The idea behind those methods is to continue to use level sets to represent the crack but to replace the resolution of PDEs to update them by simple geometric equations. A 2D comparison of several geometric updates (including the one of Ventura et al. [49]) and Hamilton-Jacobi based updates can be found in Duflot [17]. In addition, it is interesting to mention the work of Colombo and Massin [51] and Colombo [52]. Colombo and Massin [51] introduced a hybrid method using geometric equations to perform the velocity extension and an explicit polynomial equation to update the level sets while preserving Hamilton Jacobi equations to reinitialize and orthogonalize the level sets. Colombo [52] later simplified this approach by replacing all of the level set evolution algorithm by geometric update equations. The approach we use in this work falls into the same category of geometric updating methods. It is very similar to the 𝜙𝜓 r𝜃 method introduced in 2D by Duflot [17] and can be considered as a simplified version of the latter. As pointed out by Duflot [17], the generalization to 3D is straightforward and will be demonstrated in the examples coupled with a geometric velocity extension method. The paper is organized as follows. In section 2, the constitutive equations and space and time discretization employed are briefly introduced. Section 3 is devoted to the evaluation of the robustness of Hamilton-Jacobi based level set updating and the new geometric evolution algorithm is introduced. Section 4 presents the fracture model and its numerical adaption for fast 3D crack propagation. Section 5 introduces the numerical results of several benchmarks of dynamic crack propagation as well as the simulation results of a complex nonplanar 3D crack propagation under impact loading. Finally some concluding remarks are provided in Section 6. All the simulations are performed with EUROPLEXUS software (http://www-epx.cea.fr), further abbreviated EPX, a fast transient dynamics program jointly owned by the French Commissariat à l’Energie Atomique et aux Energies Alternatives and the Joint Research Center of the European Commission. The software development is carried out through a Consortium involving the co-owners and major industrial partners.

u(x, 0) = u0 (x)

x ∈ Ω,

𝜕u (x, 0) = u̇ 0 (x) 𝜕t

(4)

x ∈ Ω,

(5)

where n is the unit outward normal to the boundary Γ of Ω, g is the Dirichlet boundary condition on Γg and h is the Neumann boundary condition on Γh that form together the boundary Γ = Γh ∪ Γg of Ω, f is the body force. u0 and u̇ 0 are the initial displacement and velocitiy. The Cauchy stress tensor 𝝈 is given by the small strain tensor 𝜺 and Hooke’s law: 1 𝜺 = ∇s u = (∇u + ∇uT ), (6) 2

𝝈 = c ∶ 𝜺.

(7)

The trial and solution spaces used to define the variational form of the equations are given as: t = {u(·, t ) ∣ u(x, t ) = g(x, t ), x ∈ Γg , u(·, t ) ∈ H 1 (Ω)},

(8)

 = {u(·, t ) ∣ u(x, t ) = 0, x ∈ Γg , u(·, t ) ∈ H 1 (Ω)}.

(9)

The weak form can be obtained from Equations (1)–(5): Given f, g, h, u0 and u̇ 0 , find u(t ) ∈ t , t ∈ [0; T ], such that ∀w ∈ :

(w, 𝜌ü ) + a(w, u) = (w, f) + (w, h)Γh ,

(10)

(w, 𝜌u(0)) = (w, 𝜌u0 ),

(11)

(w, 𝜌u̇ (0)) = (w, 𝜌u̇ 0 ),

(12) 2

where the standard notation for the L inner product have been used:

(w, 𝜌ü ) = a(w, u) =

∫Ω ∫Ω

(w, h)Γh =

̈ Ω, w · 𝜌ud

(13)

𝜀ij (w)cijkl 𝜀kl (u)dΩ,

(14)

∫Γh

w · hdΓ.

(15)

2.2. Discrete form Galerkin’s semi-discrete form of the elastodynamic problem can be written using the previous equations: Given f, g, h, u0 and u̇ 0 , find uh = vh + gh , uh (t ) ∈ th such that ∀wh ∈  h :

(wh , 𝜌v̈ h ) + a(wh , vh ) = (wh , f) + (wh , h)Γh − (wh , 𝜌g̈ h ) − a(wh , gh ), (16) (wh , 𝜌vh (0)) = (wh , 𝜌u0 ) − (wh , 𝜌gh (0)),

(17)

2. Constitutive equations

(wh , 𝜌v̇ h (0)) = (wh , 𝜌u̇ 0 ) − (wh , 𝜌v̇ h (0)).

(18)

2.1. Continuous form

Using Equations (16)–(18) and standard finite element approximation for uh and wh , the classical matrix form is obtained: M Ü + KU = F

The strong form equations of the elastodynamic problem for a given body Ω are (cf, e.g., Hughes [53]): Let f ∶ Ω×]0; T [→ ℝ3 , g ∶ Γg ×]0; T [→ ℝ3 , and h ∶ Γh ×]0; T [→ ℝ3 , find u ∶ Ω → ℝ3 such that:

𝜕2u 𝜌 2 = div𝝈 + f 𝜕t u=g

𝝈·n=h

on Ω × [0; T [,

on Γg ×]0; T [, on Γh ×]0; T [,

t ∈]0; T [,

(19)

U (0) = U0 ,

(20)

(1)

U̇ (0) = U̇ 0 ,

(21)

(2)

where M and K are the mass and stiffness matrix, F is the vector of applied forces and Ü , U̇ , U are the acceleration, velocity and displacement vectors.

(3) 2

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

2.3. Temporal discretization: central difference explicit scheme

completely cut by the crack. The enrichment function is the generalized Heaviside function: { +1 if x ≥ 0, (29) Hgen (x) = −1 if x < 0.

The previous semi-discrete equations constitute a set of coupled ordinary differential equations that need to be integrated in time. We choose to use Newmark scheme for this purpose: Un+1 = Un + Δt U̇ n + Δt 2 (

1 − 𝛽)Ü n + Δt 2 𝛽 Ü n+1 , 2

̈ n + Δt 𝛾 Ü n+1 , U̇ n+1 = U̇ n + Δt (1 − 𝛾)U

As the approximation spaces change when the crack evolves as new enrichment functions are added with their associated DOFs, the stability of the time integrator has to be taken into account. This has been studied in details by Réthoré et al. [32] for discontinuous and singular tip enrichment. We use the same approach, that is all the kinematic fields are enriched in the same way. In this case, the enriched FE problem ca be written exactly in the same way as the standard FE one given in Equations (24)–(26). As shown by Réthoré et al. [32] the conservation of the energy during the simulation when the crack evolves, as well as the stability of the scheme, is preserved when the following method is adopted. When the crack evolves all the old enriched DOFs are kept and new ones are created if needed in accordance with the new crack position. Since we only use discontinuous enrichment, this means that new enriched DOFs are added if the corresponding nodes have their support entirely cut by the crack.

(22)

(23)

where 𝛾 and 𝛽 are the two parameters of the scheme. Stability and convergence order depend on the values of those two parameters. In this work, we use the central difference scheme which is equivalent to Newmark’s scheme for 𝛾 = 12 et 𝛽 = 0. This scheme is explicit, second order accurate in time and conditionally stable. The discrete equations are: Un+1 = Un + Δt U̇ n +

Δt 2 ̈ Un,

(24)

2

M Ü n+1 = Fn+1 − Fnint , +1 U̇ n+1 = U̇ n +

(25)

Δt ̈ (Un + Ü n+1 ).

(26)

2

3. Crack representation and propagation model

= KUn+1 . For a non-linear For a linear elastic material, we have Fnint +1 model, such as plasticity, internal forces must be calculated from the stress tensor. The time integrator being conditionally stable in time, the time step must respect the Courant-Freidrich-Lewy condition: Δt ≤ Δtc =

2

𝜔max

,

For 3D crack simulation using X-FEM, the representation of the crack surface and front, that can be non planar and curved, is much more complicated than in the 2D case. The use of an implicit representation of the crack with level sets has been adopted since the early work of Stolarska et al. [13], Moës et al. [14], Gravouil et al. [15], Sukumar et al. [16]. It consists in representing a signed distance field on the whole domain instead of representing the given interface Γ, which in our case is the crack surface, see Equation (32).

(27)

where 𝜔max is the maximal eigenvalue of the system. Looking at Equation (24), one can see that the update of the displacement field is explicit and therefore solving the balance of momentum does not require any iteration whatever the material model. Moreover if the mass matrix is lumped, solving Equation (25) is trivial and unexpensive which is why it is adopted in most commercial codes. The main drawback of this approach is that the time step has an upper bound called the critical time step Δtc - that can be very small. Since in this work we focus on the modeling of dynamic crack growth under transient loadings, the characteristic time of the phenomenon we are trying to capture is of the same order as the critical time step, this is not a big issue.

𝜙(x, t ) = sign((x − xΓ ).n).|x − xΓ |,

(30)

‖∇𝜙‖ = 1,

(31)

Γ(t ) = {x, 𝜙(x, t ) = 0}.

(32)

In order to represent the crack surface and front, two level sets functions are used: 𝜙t for the surface and 𝜙n for the front Γ that is represented by the intersection of the two iso-zeros, see Equation (33). This means that those two level set functions have to be orthogonal two each other.

2.4. Space discretization: the X-FEM

Γ(x) ≡ (𝜙t (x) ≤ 0) ∩ (𝜙n (x) = 0).

The X-FEM has been introduced for the first time by Black and Belytschko [1], Moës et al. [2]. It is a finite element method based on a local partition of unity (see Babuška and Melenk [4]) that introduces enrichment functions to improve the approximation properties of the discrete spaces. If one considers a set of nenr enrichment functions 𝜓 𝛼 (x), the displacement field is given by:

Gravouil et al. [15] used the same mesh as for the structure to represent the level sets. The two level sets evolve to represent the crack propagation using Hamilton-Jacobi equations solved with finite differences as introduced by Sethian [54], Osher and Fedkiw [55]. In order to simplify the resolution of these equations with finite differences, Prabel et al. [18] in 2D and Rannou et al. [19] in 3D have introduced an independent regular finite difference mesh to evolve the two level sets. The distance fields are projected onto the mesh of the structure when needed. The fracture model, which is the starting point of the evolution algorithm will be presented in the next section. In order to use this model, the crack front, which is a curve in the 3D case, is discretized with crack front points. This discretization is done by projecting the level set mesh nodes around the crack front onto the front itself following Rannou [56]. This method is very easy to implement and makes it straightforward to project back the crack velocity and angle coming front the fracture model onto the nodes of the level set mesh around the crack front. The relative size between the FE and auxiliary mesh suggested by Prabel et al. [18], Rannou et al. [19] is followed in the rest of the paper:

ui (x) =

∑ A ∈

uA N A (x) + i

∑

nenr

∑

B∈enr 𝛼=1

bB𝛼i N B (x)𝜓𝛼 (x),

(28)

where  is the set of all the FE nodes, NA (x) are the standard FE shape functions, uA is the ith component of the associated DOF A, enr is the i set of enriched nodes, and bB𝛼i is the ith component of the enriched DOF B for the 𝛼 th enrichment function. Based on previous work combining explicit time integrators with X-FEM (see Elguedj et al. [24], Gravouil et al. [25], Menouillard et al. [27,29]) in order to maximize the critical time step and limit the increase of DOFs as the crack evolves, we only use Discontinuous enrichment to represent the displacement jump of elements that are 3

(33)

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

∙ the auxiliary mesh size should be sufficiently small to take into

to compute accurately the crack local frame along the front as well as the enrichment functions. Geometric methods rely on level sets, but instead of solving partial differential equations to evolve the level sets, they introduce geometric equations that are approximate away from the crack surface or front but are accurate enough where needed, see for example Duflot [17], Ventura et al. [49], Agathos et al. [50], Colombo [52], Menouillard [58]. The main advantages of such methods is that the geometric equations are directly used to compute the new values of the level sets without requiring to solve any PDE. The method we introduce in the following paragraphs is very similar to the 𝜙𝜓 r𝜃 method introduced in 2D by Duflot [17] but somewhat simpler and is used in 3D coupled with a geometric velocity extension introduced by Rannou [56].

account curvature effect, we typically use a value of roughly one tenth fo the curvature radius, ∙ in order to keep accurate crack geometry information during the projection step, the two mesh size should be close to each other, we typically use an auxiliary mesh twice finer compared to the FE one.

3.1. The limitations of Hamilton-Jacobi based level set updating As briefly presented in the previous section, the propagation of an arbitrary 3D crack using level sets is generally done using HamiltonJacobi equations. Assuming that the appropriate crack propagation law is used, the input of the level set evolution algorithm is the velocity of the crack front given in vector form. The method is presented in details in Algorithm 1. Even though this method has been used quite successfully in the literature, it has several difficulties. Hamilton-Jacobi equations are difficult to solve and require algorithms (finite differences or stabilized finite elements) that can take a large number of iterations to obtain a satisfactory solution. Moreover, in the case of two level sets used to represent the crack, keeping those two functions orthogonal to each other and defined as signed distance functions at the same time is most of the time contradictory. This is particularly true for complex crack propagations in dynamics where the crack surface curvature can be important and the crack front have strong and frequent changes in direction. This results in a loss of robustness in the update algorithm and therefore a badly described crack. This is illustrated in the example given in Fig. 1. We consider an initial planar straight crack in a cuboid domain. We impose a non-planar crack propagation by prescribing an angle of +50◦ with respect to the initial crack plane on one side and of −50◦ on the other side. The iso-zero of 𝜙t is represented in red, and the iso-zero of 𝜙n in blue; other iso-values of both level sets are shown to better understand their evolution. Fig. 1(a) shows the initial crack and the imposed propagation directions, Fig. 1(b) shows the two level sets before solving the orthogonalization equation and Fig. 1(c) after solving this equation. We observe that 𝜙n changes a lot with respect to its initial value (this is illustrated by the black arrows) during the orthogonalization solving. We can see in Fig. 1(d) the two level sets after solving the reinitialization equation for 𝜙n . We can see that 𝜙n changes a lot again and more or less goes back to its initial position. Solving the orthogonalization equation imposes strong changes in the level sets values that are inconsistant with the reinitialization equations that will later cancel them. This means that a large number of iterations to converge those two equations are unnecessary. For such non planar cracks, satisfying at the same time the orthogonality condition ∇𝜙n .∇𝜙t = 0 and the signed distance function property (∥∇𝜙n ∥ = 1) cannot be reached at the same time. One of theses two equations must be favored, but this implies for long propagation that errors will be accumulated in the crack description and this will affect the final crack path. Similar observations have been made by Colombo and Massin [51].

3.2.1. Geometric velocity extension As for the algorithm based on Hamilton-Jacobi equations, the first step of the geometric update consists in extending the velocity field on the whole domain. For any given node M in the level set mesh, we look for the two closest interpolation points on the crack front. Node M is projected on the crack front segment formed by those two points, which gives us point H. If H is between the two front points, the velocity of point M is defined by the geometric projection on the crack front segment using Equation (42), see Fig. 2. If point H is not between the two front points, the velocity for point M is set to the velocity of the closest front point. ȧ M =

d1 d2 ȧ + ȧ . d1 + d2 2 d1 + d2 1

(42)

3.2.2. Level set geometric updating The geometric updating method that we present here was designed with the following constraints:

∙ ∙ ∙ ∙

robustness even for long crack propagations, numerically inexpensive, adapted to complex crack geometries, precise localization of the crack in the case of strong direction changes.

Similarly to other geometrical methods [17,49,50,52], in the present method, we geometrically decompose the space into different zones based on the new crack segments geometry. Different geometric formulas are used in each of these zones to update both level sets. As stated in the previous section, the orthogonality condition and the signed distance function cannot be enforced at the same time. In the present method, the orthogonality condition is favored in order to accurately represent the local crack frame in all situations. Once the crack velocity and direction has been computed by the fracture criterion and extended to the whole domain, we divide the domain into four zones that have different update equations, see Fig. 3. The scalar velocity and crack front direction are computed using the following equations: ) ( V𝜙t 𝜃 = arctan , (43) V𝜙n √ ȧ = V𝜙2 + V𝜙2 , (44)

3.2. Geometric updating of level sets In order to circumvent the problem illustrated in the previous section and to reduce the numerical cost of solving Hamilton-Jacobi equations, several methods have been used in the literature. They essentially fall into three classes: geometric, purely explicit or implicitexplicit methods. Implicit-explicit methods consist in representing the crack with a 2D mesh and level sets at the same time. The crack is evolved by adding new elements to the mesh and new level sets are created at each propagation step, see for example, Prabel et al. [20], Fries and Baydoun [21]. Purely explicit methods are similar to the previous ones, except that they only rely on the explicit representation of the crack with a 2D mesh, see for example Duarte et al. [47], Garzon et al. [48]. Nevertheless, they loose the advantages of using level sets

t

n

The definition of the four zones and the evolution equations are summarized in Table 1. The zones and their geometric evolution are defined such that the following holds:

∙ 𝜙t does not change in zone 1, ∙ 𝜙t is rotated by the angle 𝜃 in zones 2 and 3, i.e., parallel to the new crack segment,

∙ 𝜙t forms a circular arc in zone 4 to ensure continuity between zones 1 and 3, while preserving its initial sign.

∙ 𝜙n advances of ȧ Δt while rotating by the angle 𝜃 in all zones. 4

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Algorithm 1 Hamilton-Jacobi evolution equations to update the crack representation with two level sets. 1. Crack propagation model provides the velocity field ȧ and crack propagation angle for each node around the crack front. V𝜙∗ = ȧ sin(𝜃), V𝜙n = ȧ cos(𝜃), t

(34)

2. Velocity extension to the whole mesh.

𝜕 V𝜙∗ 𝜕𝜏

𝜕 V𝜙∗ 𝜕𝜏

t

+ sign(𝜙t )n𝜙t .∇V𝜙∗ = 0,

𝜕 V𝜙n + sign(𝜙n )n𝜙n .∇V𝜙n = 0 𝜕𝜏

(35)

+ sign(𝜙n )n𝜙 .∇V𝜙∗ = 0,

𝜕 V𝜙n + sign(𝜙t )n𝜙 .∇V𝜙n = 0. t 𝜕𝜏

(36)

t

t

n

t

3. Modification of the crack surface velocity, ∗

V𝜙 t =

𝜙n H (𝜙n ) V𝜙t ∇𝜙t . Δt V𝜙n

(37)

4. Update of the crack surface [17] and front [15,57],

𝜕𝜙n 𝜕𝜙t + V𝜙t · ∇𝜙t = 0, + ‖∇𝜙n ‖V𝜙n = 0. 𝜕t 𝜕t

(38)

5. Reinitialization of the crack surface,

𝜕𝜙t + sign(𝜙t )(‖∇𝜙t ‖ − 1) = 0. 𝜕𝜏

(39)

6. Orthogonalization of the crack front,

𝜕𝜙n + sign(𝜙t )n𝜙 ∇𝜙n = 0. t 𝜕𝜏

(40)

7. Reinitialization of the crack front,

𝜕𝜙n + sign(𝜙n )(‖∇𝜙n ‖ − 1) = 0. 𝜕𝜏

𝜙t =

min(‖(𝜙0t )‖, ‖(𝜙0t cos(𝜃) − 𝜙0n

sin(𝜃))‖)

𝜙t = 𝜙0t cos(𝜃) − 𝜙0n sin(𝜃), 𝜙t =

√

(𝜙0t )2 + (𝜙0n )2 sign(𝜙0t ),

𝜙n = (𝜙0n cos(𝜃) + 𝜙0t sin(𝜃)) − ȧ Δt .

sign(𝜙0t ),

(41)

(top side of the initial crack tip on Fig. 3). However, the two methods are different regarding 𝜙n . The 𝜙𝜓 r𝜃 method results in 𝜙n beeing orthogonal to 𝜙t everywhere, whereas in our method 𝜙n is only orthogonal to the new crack segment everywhere. These two differences are such that the new 𝜙n is parallel to the new crack segment in a bigger region with the present method and that the orthogonality condition of 𝜙t with the new crack segment is ensured even far away from the crack. In our experience this is particularly important when there are kinks or zigzag patterns in the propagation as will be seen in the next paragraphs. We can also compare the present method with the one of Colombo [52]. In Colombo’s method, a freezing condition is applied to 𝜙t for the region where 𝜙0n ≤ 0, which corresponds to zones 1 and 2. The present method also has this freezing condition but only in zone 1. In Colombo’s method, 𝜙t is parallel to the new crack segment outside the freezing zone which corresponds to zones 3 and 4. In the present method, we set 𝜙t to be parallel to the new crack segment in zones 2 and 3. In both methods, 𝜙n is updated everywhere so that it is orthogonal to the new crack segment. In Colombo’s methods, 𝜙t is fully discontinuous at the boundary of the two zones, whereas in the present method, it is only discontinuous where 𝜙t > 0. As pointed by Colombo [52], we have observed that this discontinuity of 𝜙t does not represent a problem for the fracture mechanics calculations.

(45) (46)

(47)

(48)

The proposed approach can be compared to the 𝜙𝜓 r𝜃 method by making the following observations. There are 5 areas in the 𝜙𝜓 r𝜃 method whereas we have only 4 zones. Area 1 is partly in zones 1 and 2, area 2 is entirely in zone 1, area 3 and zone 4 are identical, area 4 is entirely in zone 3, area 5 is partly in zones 2 and 3. In particular, it should be noted that the boundary between zones 1 and 2 is orthogonal to the new crack segment, whereas the boundary between area 1 and 5 is the bisector between the old and new crack segments. The update equations of 𝜙t produce similar results: orthogonal to the old crack segment backward (left side of the initial crack tip on Fig. 3), orthogonal the new one forward (right side of the initial crack tip on Fig. 3), circular connection on one side (bottom side of the initial crack tip on Fig. 3)) and discontinuous on the other side 5

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Fig. 1. Illustration of the lack of robustness for the update of 𝜙n using Algorithm 1.

Fig. 2. Projection of a level set node onto the crack front for velocity extension.

A graphical comparison of the final level sets for the example of Fig. 3 with the present method, the 𝜙𝜓 r𝜃 method of Duflot [17] and Colombo’s method [52] can be observed in Fig. 4.

3.3. Evolution examples: robustness of the approach In this section, we evaluate the performance and the robustness of the geometric update algorithm on several propagation cases. A cuboid domain with a regular 30 × 30 × 16 hexahedral mesh is used for all simulations. In all cases, several propagation steps are simulated in various configuration. For simplicity the crack extension Δa is equal to the mesh size at each step.

Fig. 3. Graphical definition of the four evolution zones. 6

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Table 1 Zone definition and update equations. Zones

Zone 1

Zone 2

Zone 3

Zone 4

Definition

− tan(𝜃)𝜙0t ≥ 𝜙0n 𝜙0n ≤ 0 𝜙0t

− tan(𝜃)𝜙0t < 𝜙0n 𝜙0n ≤ 0

− tan(𝜃)𝜙0t ≤ 𝜙0n 𝜙0n > 0

− tan(𝜃)𝜙0t > 𝜙0n 𝜙0n > 0

Equation (45)

Equation (46)

Equation (47)

Update of 𝜙t Update of 𝜙n

Equation (48)

Fig. 4. Graphical comparison of the geometric level set updating methods for the example of Fig. 3: (a) initial level sets, final level sets for (b) present method, (c) 𝜙𝜓 r𝜃 method of Duflot [17], (d) method of Colombo [52].

3.3.1. Helical propagation This example is identical to the one used in the previous section to demonstrate the issues with Hamilton-Jacobi based level set evolution. The results are shown in Fig. 5, with the level sets before propagation, after one step, and after ten steps. We can see that the propagation is well described and that even after a total propagation that doubles the crack length the iso-zeros of the two level sets are smooth and do not show any discrepancy.

gation with an angle of +50◦ followed by propagation with an angle of −50◦ and so on. The results are shown in Fig. 6 were we plot the two level sets before propagation, after the first step and after the tenth step. Again, we can see the the crack description remains smooth and does not show any discrepancy. We can see in Fig. 6(c) that 𝜙n is curved because of the geometric update equations away from 𝜙t but despite that all the good properties are preserved close to the crack front.

3.3.2. Zigzag propagation In this example we impose successive strong changes in the propagation direction of a planar crack. We successively impose a propa-

3.3.3. V-shaped front The next example consist in the straight planar propagation of a vshaped crack front. The results are shown in Fig. 7 were we plot the two

Fig. 5. Helical level sets evolution (a) before propagation, (b) after the first step, (c) after ten steps.5(c).

Fig. 6. Zigzag level set evolution (a) before propagation, (b) after the first step, (c) after ten step. 7

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

where 𝜎 ij are the components of the Cauchy stress tensor and 𝜔 is a weighting function than can take the following form:

𝜔(r ) = e

−𝛼

r2 R2

.

(50)

The integration domain used to compute the equivalent stress tensor is a half-disk in 2D and a half-sphere in 3D. As we want a 3D criterion that can take into account tensile and shear driven fracture, we use the approach introduced by Schollmann et al. [62] for the tensile plane stress case and generalized to full 3D and shear driven fracture by Haboussa [63]. The propagation is first triggered by the following equations: { 𝜎 ̃ < 𝜎Ic (𝜀) ̇ ȧ = 0 no propagation, (51) 𝜎 ̃ ≥ 𝜎Ic (𝜀) ̇ ȧ ≠ 0 propagation.

̃ is the maximum principal stress of the equivalent stress tenwhere 𝜎 sor and 𝜎 Ic is the analogous to the fracture toughness and has to be identified. Following Haboussa et al. [23], Schollmann et al. [62], we replace the stress intensity factors by the corresponding components of the equivalent stress tensor:

Fig. 7. V-shaped front evolution (a) after the first step, (b) after ten steps.

⎧KI ⎪ ⎨KII ⎪ ⎩KIII

→ → →

𝜎 ̃22 , 𝜎 ̃12 , 𝜎 ̃23 ,

(52)

We also introduce the normalized equivalent stress components:

Fig. 8. Bifurcation of a curved front with an initial angle of 70 (a) before propagation, (b) after ten steps.

level sets after the first step and after the tenth step. Again the results are very satisfactory (see Fig. 7).

|𝜎 ̃23 | |𝜎 ̃12 | ,𝜎 ̃n = . |𝜎 ̃22 | + |𝜎 ̃12 | + |𝜎 ̃23 | III |𝜎 ̃22 | + |𝜎 ̃12 | + |𝜎 ̃23 |

(53)

𝜎 ″ 22 = 𝜎zz sin2 (𝜓) + 2 cos(𝜓) sin(𝜓)𝜎𝜃z + 𝜎𝜃𝜃 cos2 (𝜓), 𝜎 ″ 33 = 𝜎zz cos2 (𝜓) − 2 cos(𝜓) sin(𝜓)𝜎𝜃z + 𝜎𝜃𝜃 sin2 (𝜓), 𝜎 ″ 12 = 𝜎rz sin(𝜓) + 𝜎r 𝜃 cos(𝜓), 𝜎 ″ 23 = (𝜎zz − 𝜎𝜃𝜃 ) sin(𝜓) cos(𝜓) + 𝜎𝜃z (2cos2 (𝜓) − 1), 𝜎 ″ 13 = 𝜎rz cos(𝜓) − 𝜎r 𝜃 sin(𝜓).

(54)

In the case of tensile driven fracture, the direction of propagation is the one that maximizes the tensile or hoop stress. This direction is obtained by solving the following equations:

4. Fracture model As we want to model the propagation of cracks under impact loadings mostly composed of metallic materials, we need a fracture criterion that has the following features:

⎧ ” ⎪ 𝜕𝜎22 ⎪ 𝜕𝜃 ⎨ ” ⎪ 𝜕𝜎22 ⎪ 𝜕𝜓 ⎩

∙ mode I, II and III influence the crack evolution law, ∙ propagation at low and high strain rates, that is tensile and shear driven fracture, To meet these requirements, we adopt a model based on the so called local approach to fracture proposed amongst others by Pineau [59]. We use a numerical adaption of this approach, similarly to many authors in the litterature, see for example Haboussa et al. [23], Menouillard et al. [27], Remmers et al. [60,61]. The model is based on the calculation of an equivalent stress tensor at the crack tip:

,

= 0,

(55)

= 0.

The following semi-analytical solution is obtained: ( ( )) √ 1 ̂ 𝜃ctensile = 2 sign(̃ 𝜎 12 ) arctan S − (̂ S )2 + 8 , 4 ) ( ̃𝜃z (𝜃ctensile ) 2𝜎 1 tensile 𝜓c = sign(̃ 𝜎 23 ) arctan , 2 𝜎 ̃𝜃𝜃 (𝜃ctensile ) − 𝜎 ̃zz (𝜃ctensile )

∙ generalized plasticity.

∫Ω 𝜔 dVM

𝜎 ̃nII =

𝜎 ″ 11 = 𝜎rr ,

3.3.4. Bifurcation of a curved front The last example consists in the bifurcation of planar crack with a circular front. We impose a propagation with an angle of 70◦ at the first step followed by several propagation steps in the same direction. The results are shown in Fig. 8 where we plot the two level sets before the first step and after the tenth step. Again the results show the robustness of the proposed method.

𝜎 ̃ij =

|𝜎 ̃22 | , |𝜎 ̃22 | + |𝜎 ̃12 | + |𝜎 ̃23 |

We consider an arbitrary 3D crack in an elastic medium and the asymptotic stress tensor around the front is given by Equation (54) in a frame parametrized by two angles 𝜃 and 𝜓 as shown in Fig. 9. The asymptotic stress tensor is given in the ″ frame that is obtained by successive rotation of the crack front frame of the angle 𝜃 around n then of the angle 𝜓 around t′1 , as shown in Fig. 9.

◦

∫Ω 𝜔𝜎ij (M ) dVM

𝜎 ̃nI =

(56) (57)

with

̂ S=

(49) 8

1+𝜎 ̃nI − (1 − 𝜎 ̃nIII )p(𝜈)

𝜎 ̃nII

and p(𝜈) =

1 √ ( 𝜋 − 5𝜈). 4

(58)

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Fig. 9. Various frames around the crack tip used to define the crack propagation law.

In the case of shear driven fracture, the direction of propagation is the one that maximizes the Von Mises stress. This direction is obtained by solving the following equations ⎧ ” ⎪ 𝜕𝜎VM ⎪ 𝜕𝜃 ⎨ 𝜕𝜎 ” ⎪ VM ⎪ 𝜕𝜓 ⎩

= 0,

The transition between shear and tensile fracture is based on an ̃eq as presented in equivalent strain measure around the crack tip 𝜀 Haboussa et al. [23], that can for example be taken as the equivalent plastic strain. Two parameters are considered: 𝜀 ̃shear and 𝜀̃tensile . ̃eq < 𝜀̃tensile the crack will follow the direction of maximum tenIf 𝜀 ̃eq > 𝜀̃shear , the crack will follow the direction of sile stress and if 𝜀 maximum shear stress. A linear transition between the two values is considered. The crack velocity is obtained using the well known relation introduced by Kanninen and Popelar [64] adapted to be used with the equivalent stress tensor at the crack tip:

(59)

= 0.

The following semi-analytical solution is obtained: 𝜎 ̃III ⎞⎞ n ⎛ ⎛𝜎 ) 3 ̃ (1 − 8 1+ = sign(̃ 𝜎 12 ) 𝜈 arctan ⎜ ⎜ I n 20 ⎟⎟ , ⎜5 ⎜ ⎟⎟ 4 145 𝜎 ̃II ⎝ ⎝ ⎠⎠ ) ( 𝜎 ̃zz (𝜃cshear ) − 𝜎 ̃𝜃𝜃 (𝜃cshear ) 1 𝜓cshear = sign(̃ 𝜎 23 ) arctan . 2 2𝜎 ̃𝜃z (𝜃cshear )

𝜃cshear

𝜋(

n

(60)

(61)

Fig. 10. Crack update for a propagation consistent with the mesh size.

Fig. 11. Experimental setup and geometry for Kalthoff’s experiment. 9

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

As we only use Heaviside enrichment in this work, the crack front can only be on element boundaries. It means that without any precaution, any small crack growth would entirely cut an element and lead to overestimated dissipation. To avoid this problem, we only extend the crack when the criterion produces a significant propagation consistent with the mesh size that we call lcar . As shown in Fig. 10 in 2D, many time steps can be needed to reach this value. To take into account the history of mechanical fields in the vicinity of the crack, the angle of propagation is not calculated with the equivalent stress tensor at the last time step, but with a ̃meanT given time weighted average of the equivalent stress tensor 𝜎 ij by: ( ) Δt −Δt + Δt 𝜎 1 𝜎 ̃meanT = − 𝜎 ̃meanT ̃ , (63) ij ij Tcumul Tcumul ij

Table 2 Material and fracture parameters of 18Ni900 for Kalthoff’s experiment. Young’s modulus Poisson ratio Density Rayleigh wave speed Yield stress Tangent modulus Critical stress Half sphere radius Spatial averaging radius Characteristic mesh length

( ) ̇ 𝜎 (𝜀) ȧ = cR 1 − Ic ,

𝜎 ̃

E

𝜈 𝜌 cR

𝜎y

Et

𝜎 Ic

R Rmoy lcar

190 GPa 0.3 7830 kg m−3 2800 m s−1 190 MPa 1600 MPa 35 MPa 7 mm 7 mm 1 mm

(62)

where cR is the Rayleigh wave speed.

Fig. 12. Kalthoff’s experiment. Von mises stress on the initial configuration and crack shape after 20 μs, 40 μs, 60 μs, 80 μs. 10

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Fig. 15. Kalthoff’s experiment. Normalized equivalent stress tensor component versus time in the center of the specimen.

Fig. 13. Kalthoff’s experiment. Crack length versus time in the centre and on a face of the 3D specimen superposed with 2D results from Gravouil et al. [25].

where Tcumul is a time counter that starts from the moment the criterion is reached and as long as it is reached, until lcar is achieved. This time counter is reset to zero after a propagation or when the criterion is not met anymore. This criterion is applied to each crack front node. Moreover, a spacial averaging of the predicted velocity field used to update the level sets is done in order to smooth the propagation along the crack front. For each crack front point, the velocity components V𝜙t and V𝜙n are averaged with the other front points that are closer to a given length Rmoy .

low strain rate, brittle failure, i.e., tensile stress driven, is observed with a simultaneous propagation of the two cracks with a global angle between 60◦ and 70◦ . If one increases the projectile’s velocity, a transition between brittle fracture and shear band propagation, i.e., shear stress driven, with a propagation angle of approximatively −10◦ occurs. We consider here only the brittle or tensile stress based propagation with an impact velocity of 16 m s−1 . As in Prabel et al. [18], Elguedj et al. [24], Menouillard et al. [29], Song et al. [30] we consider only the upper half of the plate with the appropriate symmetry boundary conditions. The material properties are those of a 18Ni1900 maraging type steel and are given in Table 2 along with the fracture parameters; the material is elasticplastic with isotropic linear hardening. We use a regular hexahedral mesh of 50 × 50 × 8 elements for the solid and a level set mesh twice finer.

5. Numerical examples of 3D dynamic crack growth 5.1. Kalthoff’s experiment This example deals with the numerical simulation of Kalthoff’s experiment of the failure mode transition under pure mode II loading [65]. The experimental configuration is presented in Fig. 11: a plate with two symmetrical edge cracks is impacted by a projectile at a given speed V0 . By modifying the projectile’s velocity, Kalthoff observed a transition in the type of failure. At low velocity, i.e., under

Fig. 14. Kalthoff’s experiment. Crack angle 𝜃 c as a function of time in the center and on a face of the specimen.

Fig. 16. Experimental setup and geometry for ZRR experiment. 11

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

center at the end of the simulation, which is expected since the stress state is close to plane strain whereas it is close to plane stress on the faces. In Fig. 14 we plot the evolution of the propagation angle versus time in the centre and on a face of the specimen. We can see that when the crack starts to propagate around t = 28 μs, the angle is roughly 65◦ followed by a straight propagation in the same direction. From t = 80 μs the angle starts to evolve again as the crack approaches the top free boundary of the specimen. This behavior is also observed in 2D in various results reported in the literature. In Fig. 15, we plot the normalized equivalent stress component versus time, which indicates the influence of each propagation mode. As expected, the mode III component is not significant during all the simulation since the problem is essentially 2D. We can also observe that, as the input wave propagate, the crack is mostly in a mode II state before propagation, which results in the initial angle observed in Figs. 12 and 14. After propagation, the mode I component is almost equal to 1, which is consistent with the crack propagation angle.

Table 3 Fracture parameters for ZRR experiment. Critical stress Half sphere radius Spatial averaging radius Characteristic mesh length

𝜎 Ic R Rmoy lcar

35 MPa 12 mm 4 mm 4 mm

Results are shown in Figs. 12–15. We show in Fig. 12 the crack surface and the von Mises stress at t = 20, 40, 60, 80 μs on the undeformed geometry. We can see that the crack propagates with an angle similar to the one observed in the experiment and in the literature. Since the problem is mostly 2D, we can see that the crack front is relatively straight. In Fig. 13, we plot the crack length versus time in the center of the specimen and on one of the faces, along with 2D plane strain results obtained from Gravouil et al. [25]. We can see that the three curves are very similar, and that in the 3D case, the crack propagates slightly faster in the

Fig. 17. ZRR experiment. Displacement norm on the initial configuration and crack shape at 0 μs, 20 μs, 40 μs, 60 μs, 70 μs. 12

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Fig. 18. Grégoire experiment. Experimental setup and geometry.

Fig. 19. Grégoire experiment. Close up view of the post-mortem specimen. Colors of the initial image from Grégoire [68] have been modified for better understanding.

5.2. Zhou Rhosakis Ravichandran experiment

5.3. 3D dynamic crack propagation: Grégoire’s experiment

Here, we are considering the experiment by Zhou et al. [66], which was modeled, amongst other authors, by Haboussa et al. [23], Li et al. [36], Armero and Linder [37], Zhou et al. [67]. The ZRR experiment is a variation of Kalthoff’s experiment in which the impacted specimen possesses only one notch. The geometry and boundary conditions are represented schematically in Fig. 16. In the ZRR experiment, a change of rupture mode can be observed in the course of the propagation. For a sufficiently high impact velocity (greater than 30 m.s−1 ), one can observe a rupture of the specimen by propagation of a shear band. At lower impact velocities, the crack first propagates in a shear band (with an angle of about −10◦ ), then bifurcates into a tensile mode with an angle of about 35◦ . For lower impact velocities, only tensile driven fracture is observed. As in Kalthoff’s experiment, we only consider the tensile driven propagation and therefore use an input velocity of V0 = 20 m s−1 . The material is the same as in the previous experiment, the fracture parameters are given in Table 3. We use a regular hexahedral mesh of 50 × 25 × 7 element with a characteristic length of about 4 mm, and a level set mesh 50% finer. The results are shown in Fig. 17, where the norm of the displacement field on the initial configuration as well as the crack faces during the simulation can be seen. The impact velocity being low, the propagation is driven by the tensile stress during the whole simulation. The obtained propagation angle is very similar to what was observed in the experiments by Zhou et al. [66] as well as other numerical simulations that can be found in the literature.

The last example is based on the experiment of Grégoire et al. introduced in [28, 68]. The experimental setup can be observed in Fig. 18. It consists of a Split Hopkinson Pressure Bar test rig. The test specimen is made out of PMMA and consists in a rectangular plate with a hole and a initial notch coming out of the hole. The hole provides a direct conversion of compressive waves coming from the impact of the input bar into tensile waves. In order to ensure mixed-mode I + II + III loading and crack orientation effects during the propagation, the initial notch is inclined in the width of the specimen with an angle of 22◦ as shown in the side view in Fig. 18. Strain gages are placed on the input and output bars. The measured signal can be interpreted in term of strength or velocity. An optical extensometer is used to measure the position of the crack tip during the experiment. Such measurement system can be used in these experiment since PMMA is transparent and with an appropriate lighting system, the contrast between the base material and the crack surface can be seen by the extensometer. It should be noted that the data acquired with this device corresponds to a mean value of the crack tip position across the width of the specimen.

Table 4 Material and fracture parameters of PMMA for Grégoire’s experiment. Young’s modulus Poisson ratio Density Rayleigh wave speed Critical stress for initiation Critical stress for propagation Half sphere radius Spatial averaging radius Characteristic mesh length Density of the bar 1D wave speed in the bar

E

𝜈 𝜌 cR

𝜎Ici 𝜎Icp R Rmoy lcar

𝜌bar

cbar

2.42 GPa 0.42 1180 kg.m−3 800 m.s−1 7 MPa 2.8 MPa 5 mm 6 mm 0.7 mm 1160 kg.m−3 1650 m.s−1

Fig. 20. Grégoire experiment. Crack length versus time: face and center of the specimen, experimental result from Grégoire [68]. 13

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Fig. 21. Grégoire experiment. Face and side view of the specimen, displacement norm on the initial configuration and crack shape at 0 μs, 200 μs, 225 μs, 250 μs, 275 μs and picture of the post-mortem specimen from Grégoire [68].

where 𝜌bar is the density of the output bar and cbar is the 1D wave speed in the output bar. The results are shown in Figs. 20–24. In Fig. 20, we plot the numerical crack length versus time in the center and on the face of the specimen as well as the experimental crack length measured with the optical extensometer. We can see that the crack initiation is well captured and that the initial crack velocity is very close to the experimental one. The numerical simulation predicts the crack arrest and restart. The duration of the crack arrest is correctly predicted, but it happens later in the simulation than in the experiment. This is probably due to the fact that the experimental velocity of this experiment is not given in Grégoire [68], therefore we have used the input data from another experiment. Although the velocity profile in the input bar is very similar, the impact velocity of the striker is slightly higher in this experiment than in the one we used for the input data (12.4 m.s−1 instead of 14 m.s−1 ). We can see that the crack is longer in the center of the specimen than on the faces of the specimen which is in accordance with the experimental observations and what was observed in the previous examples. Fig. 21 shows face and side view of the specimen with the crack surface and norm of the displacement field at t = 0, 200, 225, 250, 275 μs as well as a picture of the post-mortem specimen with a similar orientation. The twisting of the crack due to mode III can be seen, and as in the experiment, the crack tends to become planar and on the symmetry plane of the specimen. Fig. 22 shows a close up view around the crack with isocontours of the von Mises stress at the same times. In Figs. 23 and 24 we plot the normalized equivalent stress tensor components in the center and on the face of the specimen versus time. We can see that the mode III contribution is around 25% when the crack starts to propagate in both cases, and that the mode II contribution is non zero only on the faces of the specimen with a contribution of 18%, which is consistent with the the twisting of the crack. As can bee seen in Fig. 21, the crack surface is planar after t = 225 μs, which is consistent with the mode II stress component value close to zero on the

This experiment is composed of three phases: the crack initiates at around 220 μs and propagates at a constant velocity of approximatively ȧ ≃ 200 m s−1 ; then the crack stops for about 50 μs; finally a second propagation stage at the same velocity occurs. After a crack extension of 20 mm, the crack front is out of the scope of the optical extensometer and no further information can be obtained. Similar experiments have already been studied using X-FEM in the literature but with a planar crack and only mode I + II involved, see Prabel et al. [18], Gravouil et al. [25], Grégoire et al. [28]. The interest of this version of the experiment is that mode III is involved: the crack is expected to twist during the propagation in order to become planar and join the symmetry plane of the specimen. This can be observed on the images of the post-mortem specimen shown in Fig. 19 taken from Grégoire [68]. A mode I + III propagation of the crack is expected in the center of the specimen and a mode I + II + III propagation on the faces of the specimen. The material and fracture parameters used are given in Table 4. In the propagation zone, we use a regular hexahedral mesh with a mesh size of 1.4 mm. The level set mesh is regular with a size of 1 mm. Following Grégoire et al. [28], we use two values for the critical stress, one for the initiation, and a smaller one for the propagation. This is because the radius of the initial notch tip is larger than that of the propagating crack tip as it is machined with a saw. The boundary conditions used for the simulation are the following: the experimental velocity is imposed on the contact surface between the input bar and the specimen and an impedance condition is imposed on the contact surface between the specimen and the output bar. A similar approach was followed by Prabel et al. [18], Gravouil et al. [25], Grégoire et al. [28]. This allows us to only model the specimen and not the input and output bars. We impose a 1D impedance condition using the following equation: Fx = 𝜌bar cbar u̇ x ,

(64) 14

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

Fig. 22. Grégoire experiment. Close up view, von Mises stress on the initial configuration and crack shape at 0 μs, 200 μs, 225 μs, 250 μs, 275 μs.

Fig. 23. Grégoire’s experiment. Normalized equivalent stress tensor component in the center of the specimen.

Fig. 24. Grégoire’s experiment. Normalized equivalent stress tensor component on the face of the specimen.

6. Conclusion faces of the specimen that can be seen at that time in Fig. 24. After t = 250 μs the propagation is mostly governed by mode I in the whole specimen.

We have introduced in this work a numerical model based on the XFEM coupled with level sets for the simulation of dynamic propagation 15

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17

of arbitrary non-planar cracks under transient loadings. The method uses only discontinuous enrichments and is integrated in time with the the explicit central difference scheme. The propagation criterion uses a locally averaged stress tensor at the crack tip and can handle mode I + II + III propagation cases. One of the main contribution of this work is the demonstration of the lack of robustness of Hamilton-Jacobi based level set updating for crack propagation modeling. To circumvent these problems, we have introduced a simple geometric updating method that does not involve any equations solving. The efficiency of the method is demonstrated on two problems commonly used in the literature to validate planar dynamic fracture simulations and show very good results. Finally a complex experiment involving non-planar mode I + II + III propagation under impact loading is computed using the proposed method and shows very good results that compare well with the experimental one.

[21] T. Fries, M. Baydoun, Crack propagation with the extended finite element method and a hybrid explicit–implicit crack description, Int. J. Numer. Meth. Eng. 89 (2012) 1527–1558. [22] J.-P. Crété, P. Longère, J.-M. Cadou, Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM, Comput. Meth. Appl. Mech. Eng. 27 (2014) 204–233. [23] D. Haboussa, T. Elguedj, B. Leblé, A. Combescure, Simulation of the shear-tensile mode transition on dynamic crack propagations, Int. J. Fract. 178 (2012) 195–213. [24] T. Elguedj, A. Gravouil, H. Maigre, An explicit dynamics extended finite element method. Part 1: mass lumping for arbitrary enrichment functions, Comput. Meth. Appl. Mech. Eng. 198 (2009) 2297–2317. [25] A. Gravouil, T. Elguedj, H. Maigre, An explicit dynamics extended finite element method. Part 2: element-by-element Stable-Explicit/Explicit dynamic scheme, Comput. Meth. Appl. Mech. Eng. 198 (2009) 2318–2328. [26] A. Combescure, A. Gravouil, D. Grégoire, J. Réthoré, X-FEM a good candidate for energy conservation in simulation of brittle dynamic crack propagation, Comput. Meth. Appl. Mech. Eng. 197 (2008) 309–318. [27] T. Menouillard, J. Réthoré, N. Moës, A. Combescure, H. Bung, Mass lumping strategies for X-FEM explicit dynamics: application to crack propagation, Int. J. Numer. Meth. Eng. 74 (2008) 447–474. [28] D. Grégoire, H. Maigre, J. Réthoré, A. Combescure, Dynamic crack propagation under mixed-mode loading - comparison between experiments and X-FEM simulations, Int. J. Solid Struct. 44 (2007) 6517–6534. [29] T. Menouillard, J. Réthoré, A. Combescure, H. Bung, Efficient explicit time stepping for the extended finite element method, Int. J. Numer. Meth. Eng. 68 (2006) 911–938. [30] J. Song, P. Areias, T. Belytschko, A method for dynamic crack and shear band propagation with phantom nodes, Int. J. Numer. Meth. Eng. 67 (2006) 868–893. [31] J. Réthoré, A. Gravouil, A. Combescure, A combined space-time extended finite element method, Int. J. Numer. Meth. Eng. 64 (2005) 260–284. [32] J. Réthoré, A. Gravouil, A. Combescure, An energy-conserving scheme for dynamic crack growth using the extended finite element method, Int. J. Numer. Meth. Eng. 63 (2005) 631–659. [33] G. Zi, H. Chen, J. Xu, T. Belytschko, The extended finite element method for dynamic fractures, Shock Vib. 12 (2005) 9–23. [34] T. Belytschko, H. Chen, Singular enrichment finite element method for elastodynamic crack propagation, Int. J. Comput. Meth. 1 (2004) 1–15. [35] T. Belytschko, H. Chen, J. Xu, G. Zi, Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment, Int. J. Numer. Meth. Eng. 58 (2003) 1873–1905. [36] S. Li, W. Liu, A. Rosakis, T. Belytschko, W. Hao, Mesh-free galerkin simulations of dynamic shear band propagation and failure mode transition, Int. J. Solid Struct. 67 (2002) 868–893. [37] J. Armero, C. Linder, Numerical simulation of dynamic fracture using finite elements with embedded discontinuities, Int. J. Fract. 160 (2009) 119–214. [38] P. Krysl, T. Belytschko, The element free Galerkin method for dynamic propagation of arbitrary 3-D cracks, Int. J. Numer. Meth. Eng. 44 (1999) 767–800. [39] C. Duarte, O. Hamzeh, T. Liszka, W. Tworzydlo, A generalized finite element method for the simulation of three-dimensional dynamic crack propagation, Comput. Meth. Appl. Mech. Eng. 190 (2000) 2227–2262. [40] G. Ruiz, A. Pandolfi, M. Ortiz, Three-dimensional cohesive modeling of dynamic mixed-mode fracture, Int. J. Numer. Meth. Eng. 52 (2001) 97–120. [41] T. Rabczuk, S. Bordas, G. Zi, A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics, Comput. Mech. 40 (2007) 473–495. [42] T. Rabczuk, S. Bordas, G. Zi, On three-dimensional modelling of crack growth using partition of unity methods, Comput. Struct. 88 (2010) 1391–1411. [43] Q. Duan, J.-H. Song, T. Menouillard, T. Belytschko, Element-local level set method for three-dimensional dynamic crack growth, Int. J. Numer. Meth. Eng. 80 (2009) 1520–1543. [44] P. Rozycki, N. Moës, E. Béchet, C. Dubois, X-FEM explicit dynamics for contant strain elements to alleviate mesh constraints on internal or external boundaries, Comput. Meth. Appl. Mech. Eng. 197 (2008) 349–363. [45] N. Sukumar, D. Chopp, N. Moës, T. Belytschko, Modeling holes and inclusions by level sets in the extended finite element method, Comput. Meth. Appl. Mech. Eng. 190 (2000) 6183–6200. [46] N. Sukumar, N. Moës, B. Moran, T. Belytschko, Extended finite element method for three-dimensional crack modeling, Int. J. Numer. Meth. Eng. 48 (2000) 1549–1570. [47] C. Duarte, I. Babuška, J. Oden, Generalized finite element methods for three-dimensionnal structural mechanics problems, Comput. Struct. 77 (2000) 215–232. [48] J. Garzon, P. O’Hara, D.C.A.W. Buttlar, Improvements of explicit crack surface representation and update within the generalized finite element method with application to three-dimensional crack coalescence, Int. J. Numer. Meth. Eng. 97 (2014) 231–273. [49] G. Ventura, E. Budyn, T. Belytschko, Vector level sets for description of propagating cracks in nite elements, Int. J. Numer. Meth. Eng. 58 (2003) 1571–1592. [50] K. Agathos, G. Ventura, E. Chatsi, S. Bordas, Stable 3d X-FEM/vector level sets for non-planar 3D crack propagation and comparison of enrichment schemes, Int. J. Numer. Meth. Eng. 113 (2018) 252–276. [51] D. Colombo, P. Massin, Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling, Comput. Meth. Appl. Mech. Eng. 200 (2011) 2160–2180.

Acknowledgments R. Pelée de Saint Maurice was partially supported by a scholarship from the French Commissariat à l’Energie Atomique et aux énergies alternatives; R. Pelée de Saint Maurice, T. Elguedj and A. Combescure were partially supported by a research contract with CEA. This support is gratefully acknowledged.

References [1] T. Black, T. Belytschko, Elastic crack growth in finite elements with minimal remeshing, Int. J. Numer. Meth. Eng. 45 (1999) 601–620. [2] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Int. J. Numer. Meth. Eng. 46 (1999) 131–150. [3] T. Strouboulis, I. Babuška, K. Copps, The design and analysis of the generalized finite element method, Comput. Meth. Appl. Mech. Eng. 181 (2000) 43–69. [4] I. Babuška, J. Melenk, The partition of unity method, Int. J. Numer. Meth. Eng. 40 (1997) 727–758. [5] J. Dolbow, N. Moës, T. Belytschko, Discontinuous enrichment in finite elements with a partition of unity method, Finite Elem. Anal. Des. 36 (2000) 235–260. [6] J. Réthoré, R. de Borst, M.-A. Abellan, A two-scale approach for fluid flow in fractured porous media, Int. J. Numer. Meth. Eng. 71 (2007) 780–800. [7] J. Réthoré, R. de Borst, M.-A. Abellan, A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks, Comput. Mech. 42 (2008) 227–238. [8] S. Shao, L. Bouhala, A. Younes, P. Nunez, A. Makradi, S. Belouettar, An XFEM model for cracked porous media: effects of fluid flow and heat transfer, Int. J. Fract. 185 (2014) 155–169. [9] C. Daux, N. Moës, J. Dolbow, N. Sukumar, T. Belytschko, Arbitrary branched and intersecting cracks with the extended finite element method, Int. J. Numer. Meth. Eng. 48 (2000) 1741–1760. [10] J. Dolbow, N. Moës, T. Belytschko, An extended finite element method for modeling crack growth with frictional contact, Comput. Meth. Appl. Mech. Eng. 190 (2001) 6825–6846. [11] R. Ribeaucourt, M. Baietto-Dubourg, A. Gravouil, A new fatigue frictional contact crack propagation model with the coupled X-FEM/Latin method, Comput. Meth. Appl. Mech. Eng. 196 (2007) 3230–3247. [12] E. Giner, M. Tur, J. Tarancon, F. Fuenmayor, Crack face contact in X-FEM using a segment-to-segment approach, Int. J. Numer. Meth. Eng. 82 (2010) 1424–1449. [13] M. Stolarska, D. Chopp, N. Moës, T. Belytschko, Modelling crack growth by level sets and the extended finite element method, Int. J. Numer. Meth. Eng. 51 (2001) 943–960. [14] N. Moës, A. Gravouil, T. Belytschko, Non-planar 3D crack growth with the extended finite element and level sets - Part 1: mechanical model, Int. J. Numer. Meth. Eng. 53 (2002) 2549–2568. [15] A. Gravouil, N. Moës, T. Belytschko, Non-planar 3D crack growth with the extended finite element and level sets - Part 2: level set update, Int. J. Numer. Meth. Eng. 54 (2002) 2569–2586. [16] N. Sukumar, D. Chopp, B. Moran, Extended finite element method and fast marching method for three-dimensionnal fatigue crack propagation, Eng. Fract. Mech. 70 (2003) 29–48. [17] M. Duflot, A study of the representation of cracks with level-sets, Int. J. Numer. Meth. Eng. 70 (2006) 1261–1302. [18] B. Prabel, A. Combescure, A. Gravouil, S. Marie, Level set X-FEM non matching meshes: application to dynamic crack propagation in elastic-plastic media, Int. J. Numer. Meth. Eng. 69 (2007) 1553–1569. [19] J. Rannou, A. Gravouil, M.C. Baietto-Dubourg, A local multigrid X-FEM strategy for 3D crack propagation, Int. J. Numer. Meth. Eng. 77 (2008) 581–600. [20] B. Prabel, T. Charras, A. Simatos, T. Yuritzinn, 3d crack propagation in inelastic material with X-FEM in Cast3m, in: ECCOMAS Thematic Conference on Extended Finite Element Methods-XFEM 2011, 2011. 16

T. Elguedj et al.

Finite Elements in Analysis and Design 151 (2018) 1–17 [61] J. Remmers, R. de Borst, A. Needleman, A cohesive segments method for the simulation of crack growth, Comput. Mech. 31 (2003) 69–77. [62] M. Schollmann, H. Richard, G. Kullmer, M. Fulland, A new criterion for the prediction of crack development in multiaxially loaded structures, Int. J. Fract. 117 (2002) 129–141. [63] D. Haboussa, Modélisation de la transition traction-cisaillement des métaux sous chocs par la X-FEM, Ph.D. thesis, Institut National des Sciences Appliquées de Lyon, 2012. [64] M. Kanninen, C. Popelar, Advanced Fracture Mechanics, Oxford University Press, Oxford, UK, 1985. [65] J. Kalthoff, Modes of dynamic shear failure in solids, Int. J. Fract. 101 (2000) 1–31. [66] M. Zhou, A. Rosakis, G. Ravichandran, Dynamically propagating shearbands in impact-loaded prenoteched plates I. Experimental investigations of temperature signatures and propagation speed, J. Mech. Phys. Solid. 44 (1996) 981–1006. [67] M. Zhou, A. Rosakis, G. Ravichandran, Dynamically propagating shearbands in impact-loaded prenoteched plates II. Numerical simulations, J. Mech. Phys. Solid. 44 (1996) 1007–1032. [68] D. Grégoire, Initiation, propagation, arrêt et redémarrage de fissures sous impact, Ph.D. thesis, Institut National des Sciences Appliquées de Lyon, 2008.

[52] D. Colombo, An implicit geometrical approach to level sets update for 3D non planar X-FEM crack propagation, Int. J. Numer. Meth. Eng. 237–240 (2012) 39–50. [53] T. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, Mineola NY, 2000. [54] J. Sethian, Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, 1996. [55] S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surface, Springer, 2000. [56] J. Rannou, Prise en compte d’effets d’echelle en mécanique de la rupture tridimensionnelle par une approche X-FEM multigrille localisé non-linéaire, Ph.D. thesis, Institut National des Sciences Appliquées de Lyon, 2008. [57] D. Peng, B. Merriman, S. Osher, H. Zao, M. Kang, A pde-based fast local level set method, J. Comput. Phys. 155 (1999) 410–438. [58] T. Menouillard, Dynamique explicite pour la simulation numérique de propagation de fissure par la méthode des éléments finis étendus, Ph.D. thesis, Institut National des Sciences Appliquées de Lyon, 2007. [59] A. Pineau, Development of the local approach to fracture over the past 25 years: theory and application, Int. J. Fract. 138 (2006) 139–166. [60] J. Remmers, R. de Borst, A. Needleman, The simulation of dynamic crack propagation using the cohesive segments method, J. Mech. Phys. Solid. 56 (2008) 70–92.

17