X-FEM Analysis of dynamic crack growth under transient loading in thick shells

Accepted Manuscript

X-FEM analysis of dynamic crack growth under transient loading in thick shells Thomas Elguedj, Yannick Jan, Alain Combescure, Bruno Leble, ´ Guillaume Barras PII: DOI: Reference:

S0734-743X(17)30860-6 https://doi.org/10.1016/j.ijimpeng.2018.08.013 IE 3158

To appear in:

International Journal of Impact Engineering

Received date: Revised date: Accepted date:

6 October 2017 7 August 2018 21 August 2018

Please cite this article as: Thomas Elguedj, Yannick Jan, Alain Combescure, Bruno Leble, ´ Guillaume Barras, X-FEM analysis of dynamic crack growth under transient loading in thick shells, International Journal of Impact Engineering (2018), doi: https://doi.org/10.1016/j.ijimpeng.2018.08.013

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

ACCEPTED MANUSCRIPT

Highlights • An existing Reissner Mindlin shell element is enriched using X-FEM to model dynamic crack propagation in thick shells. • A diagonal mass matrix is used for the enriched displacement and rotational DOFs

CR IP T

• An explicit representation of the crack lying on the shell?s mid-surface is done using a 1D mesh

AC

CE

PT

ED

M

AN US

• A crack propagation criterion able to account for tensile and shear driven fracture is used.

1

ACCEPTED MANUSCRIPT

X-FEM analysis of dynamic crack growth under transient loading in thick shells Thomas Elguedja,∗, Yannick Jana,b , Alain Combescurea , Bruno Lebl´eb , Guillaume Barrasc a Univ

Lyon, INSA-Lyon, CNRS UMR5259, LaMCoS, F-69621, France Group Research, D´ epartement Dynamique des Structures, F-44340, France c DGA Techniques Navales, F-83055 France

CR IP T

b Naval

Abstract

This paper is devoted to the modeling of crack propagation under impact loadings in thick shells composed of metallic materials using X-FEM. The proposed thick shell element is based on the existing Q4γ24 4-node

AN US

shell element for fast transient loadings that is enriched in the present work using X-FEM. The shell is considered to be always cut by a through the thickness crack. This element is hence enriched only with jump Heaviside functions on all degrees of freedom (displacements as well as rotations). The mass matrix corresponding to all added DOF is simply a copy of the usual continuous DOF diagonal mass matrix. The crack is discretized explicitly with a simple 1D mesh that lives on the shell’s mid-surface and independently of the finite element mesh. Plasticity as well as stress field used for crack propagation criterion is evaluated

M

using 5 Simpson points across the thickness. The crack propagation criterion is based on the measure of an equivalent stress at the crack tip and can predict both tensile and shear driven fracture as well as transitions

ED

between those two regimes. Comparisons with elastoplastic crack propagation experiments involving fracture under transient loadings show that the method is able to reproduce experimental fracture quite well.

CE

1. Introduction

PT

Keywords: X-FEM, shells, dynamic crack growth, explicit dynamics.

The eXtended Finite Element Method (X-FEM) developed in the early 2000s by Black and Belytschko

AC

[1], Mo¨es et al. [2] to perform crack propagation without remeshing has now achieved a certain maturity. Many contribution were since done for fatigue (see, e.g., [3–13]) and dynamic crack propagation in 2D or

5

3D volumetric problems (see, e.g., [14–24]). For the case of dynamic fracture, research work is often limited to brittle fracture simulations driven by the maximum hoop stress. In this context, fracture models require a relatively fine mesh in the vicinity of the crack tip in order to capture enough information and thus predict a correct crack propagation. For large scale thin-walled structures such as ship’s hull in which long ∗ Corresponding

author Email address: [email protected] (Thomas Elguedj)

Preprint submitted to International Journal of Impact Engineering

August 24, 2018

ACCEPTED MANUSCRIPT

propagation may occur, the use of a fine mesh with 3D elements would be too expensive. Therefore, it 10

seems relevant to consider an X-FEM shell element. Furthermore, experiments performed on shells made of ships’ structural steel indicate that the fracture model cannot be limited to brittle propagation. In this paper, we propose to create an enriched shell element for dynamic crack growth called X-Q4γ24 - based on the existing Q4γ24 shell element [25, 26] - using the X-FEM and a criterion which enables the simulation of crack propagation under mixed mode fracture driven by tensile or shear. The main idea behind the X-FEM is the incorporation of special enrichment functions into a standard

CR IP T

15

finite element approximation. This can be done using the partition of unity presented in Melenk and Babuska [27]. Thus, discontinuous basis functions are added to standard polynomial basis functions for nodes that belonged to elements that are intersected by a crack to provide a basis that included crack opening displacements. Other basis functions may also be added for example to locate the crack tip in an element. A key advantage of the X-FEM is that in such problems the finite element mesh does not need to

AN US

20

be updated to track the crack path.

Relatively little research has focused on developing crack growth in thin structures using the X-FEM. Mostly Kirchhoff-Love plates and shells have been considered, see for example, Areias and Belytschko [28], Lasry [29], Larsson et al. [30], Chau-Dinh et al. [31]. However, if we want to perform simulations 25

on structures such as submarines that have a significant thickness, the Reissner-Mindlin plate theory seems

M

mandatory. In comparison to the Kirchhoff-Love theory, the Reissner-Mindlin theory allows transverse shear strains through the thickness of the plate. Yet the complication of shear-locking has to be taken care of.

ED

The X-FEM was used with Reissner-Mindlin theory by Dolbow et al. [32] on the MITC4 element to consider elastic bidimensional plates. Therefore, the present work introduces, to the best of our knowledge, the first 30

X-FEM shell element to account for transverse shear.

PT

Concerning fracture modeling, we are dealing with a metallic material subjected to transient and intense loadings. This means we clearly are out of the scope of validity of linear elastic fracture mechanics. In

CE

this context, we need a fracture criterion in which the crack evolution law is influenced by all three modes, with tensile and shear driven fracture and also generalized plasticity. The local approach to fracture which 35

relies on the analysis of mechanical fields of highly solicited regions seems to be an appropriate choice. First

AC

contributions to this method can be found in Ritchie et al. [33]. The criterion was modified for dynamic purpose by evaluating mean quantities in the vicinity of the crack tip instead of local quantities (see [34–38]). In Haboussa et al. [39], the authors propose a new crack propagation criterion which is suitable for shear, and an algorithm which is capable of handling the transition from shear mode to tensile mode in the same

40

simulation. This paper is organized as follows. In the next section, we describe the standard shell element chosen and its modifications to be used in the X-FEM framework. Section 3 describes the fracture model for dynamic crack propagation in our shell element. Numerical examples are given in Section 4 and blast tests performed 2

ACCEPTED MANUSCRIPT

on metallic shells are compared to their simulations in Section 5. Finally some concluding remarks are 45

provided in Section 6.

2. X-FEM shell model 2.1. Shell model : the Q4γ24 element

CR IP T

A shell is a solid defined by a reference surface and a thickness small compared to the other dimensions of the structure. Let the local basis be (t1 , t2 , n) where n is the normal vector to the shell and (x, y, z) the 50

curvilinear coordinates of a point. For small displacements and small strains, the displacement field UM of a point M in the thickness of the shell is given by the displacements UG of a point G on the midsurface, and the rotation of this surface θ = θx t1 + θy t2 with respect to both local axis t1 et t2 .

AN US

The displacement field can thus be written :

UM = UG + θ ∧ z n = UG (x, y) + z θx (x, y)t2 (x, y) − z θy (x, y)t1 (x, y), 

UMx (x, y, z)

   UMy (x, y, z)  UMz (x, y, z)



UGx (x, y)

     =  UGy (x, y)   UGz (x, y)





−θy (x, y)

     + z  θx (x, y)   0



  . 

(2)

The standard shell element used in the present investigation is the exisiting Q4γ24 developed in Batoz and

ED

55



M

or also :

(1)

Dhatt [25], Zeng and Combescure [26]. Its approximation is based on the 4-node isoparametric quadrilateral with 6 degrees of freedom (DOF) per node : three local displacements (ux , uy , uz ) and three local rotations

PT

(θx , θy , θz ). It also has constant thickness.

This element is based on the Reissner-Mindlin theory that takes into account transverse shear strains. It is intended for thick plates in which the normal to the mid-surface remains straight but not necessarily

CE

60

perpendicular to the mid-surface. Thus, the local rotations are additional unknowns of the problem. When considering finite element approximations of relatively thin plates, it is necessary to address the

AC

phenomenon known as shear locking. As the plate becomes very thin, transverse shear strains vanish in the domain. When there are not enough functions in the finite element subspace to satisfy this constraint,

65

a poor approximation results for the plate displacements. To tackle this issue, we use the assumed strains method, see for example Dvorkin and Bathe [40], Belytschko and Bindeman [41], for our shell element. Let us assume that we have a triangulation of the domain Ω and that each element is a quadrilateral. On each element, both displacements and rotations are expressed in terms of the four nodal coefficients through

3

ACCEPTED MANUSCRIPT

the classical bilinear shape functions NA :

ui (x) =

X

NA (x)uiA ,

(3)

X

NA (x)θiA ,

(4)

A∈N

θi (x) =

A∈N

UM (x) =

X

(NA (x)UGA ) + z

A∈N

X

A∈N

CR IP T

with N the set of nodes of the element. With the discretization of eq. 1, we now have : (NA (x)θxA ) t2 − z

X

(NA (x)θyA ) t1 .

(5)

A∈N

Since the element can be curved, a local basis is defined at each node as shown on Figure 1. To simplify the problem, the surface normal at a random point is here obtained by interpolation of the four node basis

UM (x) =

X

(NA (x)UGA + zNA (x)θxA t2A − zNA (x)θyA t1A ) .

(6)

CE

PT

ED

M

A∈N

AN US

using the bilinear shape functions. As a result, the displacement field becomes :

Figure 1: Q4γ24 - 4-node curved element with 6 DOFs per node.

The numerical integration is performed with 4 Gauss points in the midsurface and 5 Simpson points

AC

70

through the thickness when plasticity is activated. 2.2. X-FEM enrichment : X-Q4γ24 element

Enriched approximation. For dynamic crack propagation, the loading conditions, the propagation parameters and the exact position of the crack are difficult to control. It has been shown, for example in Menouillard 75

et al. [19], Remmers et al. [21], Haboussa et al. [39], Areias and Belytschko [42], Song et al. [43], Menouillard and Belytschko [44], Menouillard et al. [45], that the use of near tip function (to completely locate the 4

ACCEPTED MANUSCRIPT

crack tip) is not mandatory as long as the mesh is enough refined near the crack tip. Moreover, using the classical near tip enrichment functions for dynamic crack growth requires proper handling of the enriched DOFs when the crack propagates in order to preserve energy stability, as demonstrated by R´ethor´e et al. 80

[17]. This requires to retain the enriched DOFs associated to the successive positions of the crack tip and has an important impact on the computational cost. Although this can be improved by retaining only one near tip function as shown by Elguedj et al. [22] and Menouillard et al. [45], the increased cost is still to be

CR IP T

considered. One can also note that, for the case of thick shells considered in the present work, the appropriate asymptotic field for the rotation as yet to be defined. Menouillard et al. [19, 45] have compared crack 85

path, crack length versus time and stress intensity factors versus time for the case of purely discontinuous enrichment. They show that for a sufficiently fine mesh the results are very similar. For most of the cases considered in the present work, a mesh size of approximately 1mm was used and proved to be sufficient to

AN US

obtain converged crack paths when an elastic-plastic material was used. For all these reasons, in the present work, we only use Heaviside enrichment to describe both displacement and rotation discontinuities. For all 90

the examples that will be shown at the end of the article, we made sure that the mesh was fine enough so as to obtained converged results. Convergence studies of the proposed method can be found in Haboussa et al. [39, 46]. The following X-FEM approximation is used for the displacement and rotation fields: X

NA (x)uiA +

X

NA (x)θiA +

A∈N

ED

θi (x) =

M

ui (x) =

PT

with H(x) =

NB (x)H(x)u∗iB ,

(7)

∗ NB (x)H(x)θiB ,

(8)

B∈N cut

X

B∈N cut

A∈N

 

X

1

if x is above the crack

 −1

if x is under the crack

,

(9)

and N the set of nodes of the element, N cut the set of nodes of cut elements, (u, θ) the standard DOF 95

CE

and (u∗ , θ∗ ) the enriched DOF. As we only use Heaviside enrichment, the crack tip can only be on element boundaries as shown on Figure 2. As can be seen from Equations (7) to (9) the classical Heaviside

AC

enrichment of the original X-FEM reference [2] was used . The use of this enrichment function may cause problems in blending element, i.e., partially enriched elements. The shifted enrichment function, introduced by Belytschko et al. [15, 47], can be used to overcome this issue when using the Heaviside step function. This approach was not used in the present work. In our experience, blending elements do not have a great

100

influence on accuracy when used for impact analysis with explicit time stepping.

The mass lumping technique introduced by Menouillard et al. [19] is used in the present study. The authors show in particular that the mass matrix for enriched DOF can be taken equal to the lumped mass 5

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 2: Enrichment strategy for a crack in the mesh.

matrix corresponding to standard DOF, thus conserving kinetic energy for rigid motions with or without a discontinuity.

M

For X-FEM shells, we not only use the same lumped mass matrix for enriched displacement, but also for enriched rotations. Finally, the [12 × 12] lumped mass matrix of the Q4γ 24 can be written as :

ED



 Mlump XF EM =

and :

CE

PT

with :

AC

105

Mlump standard

0

0

Mlump enriched



 Mlump standard = 

 Mlump enriched =

mfuem

0

0

mfθ em

em mxf u∗

0

0

em mxf θ∗





,

, 

,

(10)

(11)

(12)

em mxf = mfuem , u∗

(13)

em mxf = mfθ em . θ∗

(14)

By using a such diagonal mass matrix with enriched elements, Menouillard et al. [19, 48] showed that em the critical time step ∆txf to ensure stability is given by : c em ∆txf = c

6

∆tfc em , 2

(15)

ACCEPTED MANUSCRIPT

where tfc em is the critical time step to ensure stability with standard elements using an explicit scheme. Numerical integration. For elements cut by the crack and enriched with the jump function H(x), discon110

tinuous field appears in the element and the standard integration rule no longer works properly. To avoid a remeshing step, a subtriangle partitioning was first used in X-FEM (see Ref. [2]) which consists in partitioning the elements cut by the crack into subtriangles. With this technique, the projection of the stress

CR IP T

and internal variable fields is inevitable as the crack propagates, which can lead to problems such as the non-conservation of the energy. Indeed, if the crack propagates in a new direction, the subtriangles asso115

ciated with the previous crack position will also be partitioned into new subtriangles to be coherent with the new crack faces. Therefore, an efficient integration scheme was proposed by Elguedj et al. [49] : a simple subquadrangle partitioning technique is used and the elements cut by the crack are subdivided into 16 subquadrangles with 16 Gauss quadrature points each. With this technique the subelements edges are

120

AN US

not compatible with the crack faces, and integration errors may appear because of the enrichment functions that are discontinuous on the crack faces. However, this error remains small due to the great number of integration points and the technique as been successfully used in [20, 22]. In order to lower the integration error, Pel´ee de Saint Maurice [50] showed that the number of integration by subquadrangle matters more that the number of subquadrangles. In the present work, we use 16 subquadrangles with 4 Gauss quadrature

plasticity is activated.

AC

CE

PT

ED

125

M

points in each subquadrangle as shown Figure 3, keeping the 5 Simpson points through the thickness when

Figure 3: Numerical integration with 16 subquadrangles and 4 Gauss points in each.

It should be noted that a new approach using polynomial functions for X-FEM [51] and recently extended

to plasticity [52] seem to exhibit good results with no projection step and exact integration. This method was not employed in the present work but isn interesting way to reduce the computational cost of numerical integration while retaining the accuracy of the simulation. This is an area of future improvement as we

130

employ five integration points through the thickness when considering plasticity.

7

ACCEPTED MANUSCRIPT

Geometric non-linearity. As the X-Q4γ24 element introduced in the is work is simply an X-FEm enriched version of the existing Q4γ24 element, its extension to geometric non-linearity is straightforward based on the existing work on the Q4γ24 element by Batoz and Dhatt [25], Zeng and Combescure [26]. For brevity, this extension to geometric non-linearity os not given here, the interested reader should refer to the above 135

mentioned references for further insight. It is worth noting, however that an updated Lagrangian approach is used, as classically done in other explicit dynamics finite element codes. The reference configuration is

CR IP T

therefore updated at each time step resulting in arbitrarily curved surfaces, even if the initial geometry of the shell is planar, depending on the loading and boundary conditions. 2.3. Crack representation 140

The representation of the crack has drawn a lot of attention in the X-FEM litterature. An explicit representation using a 1D mesh was used in the 2D case [1, 2] but proved to be complex in the 3D case

AN US

where an implicit representation using level sets allows to represent arbitrary non planar 3D cracks, see e.g. Mo¨es et al. [53], Gravouil et al. [54]. However, the advantages of such an implict representation had to be balanced with the increased numerical cost to evolve the two level sets used to represent the crack surface 145

and front when the crack advances. In more recent works, hybrid methods, also called implicit-explicit crack representation, have been used for example by Fries and Baydoun [9]. In such approaches, a mesh is used

M

to evolve the crack and level sets are used to compute the local crack frame and the enrichment functions. In the case of shells [28, 30, 31, 55–57], most authors do not mention whether they use an implicit, explicit

150

ED

or hybrid representation of the crack. Song and Belytchko [38] claim to represent the discontinuity in the reference configuration by a level set implicit function in the shell midsurface without giving further details. Crack description. Since in the present contribution, we only consider through the thickness cracks, we

PT

choose to use an explicit representation of the crack with a simple 1D mesh. Although it requires some additional efforts to compute enrichment functions and the crack local frame, this extra cost is rather small

155

CE

in comparison to level sets evolution algorithms such as the ones used in [20, 53, 54]. In addition, as we only use Heaviside enrichment, the crack tip can only be on element boundaries, which implies that the crack only needs to be represented with one segment in each element. This implies that we only need a limited

AC

number of information to represent the crack within each element as shown in Figure 4: • element # iel, • input side #,

160

• relative position of the input point on the corresponding side αs , • exit side #, 8

CR IP T

ACCEPTED MANUSCRIPT

Figure 4: Crack representation using 1D elements and closeup view of the needed parameters to define the crack position within

AN US

a cut element.

• relative position of the exit point on the corresponding side βs .

A simple example of a crack crossing 9 elements of a given mesh is shown in Figure 5 as well as the data

AC

CE

PT

ED

M

structure and parameter values of the explicit representation in Table 1.

Figure 5: Crack crossing 9 elements of a larger mesh and described by the parameter values of Table 1.

165

Crack advance. When the crack grows in a planar thin structure, e.g. a plate, or when the elements in the mesh remain planar, advancing the crack by adding an extra segment is straightforward. It is more complex when the element is curved, as shown on Figure 6.This is of great important as geometric non-linearities are 9

ACCEPTED MANUSCRIPT

considered with an updated Lagrangian approach: even in the case of an initial planar shell, curvature as to

AN US

CR IP T

be taken into account as the shell deforms. The method we propose to follow is presented in this paragraph.

Figure 6: Representation of the case where a crack must grow on a curved non-planar shell mesh.

170

We consider a crack whose front is named A, with a local frame given as (~t1 , ~t2 , ~n) that has to grow on

M

a curved non-planar element, see Figure 6. The crack growth criterion gives us the crack propagation angle in the frame (A, ~t1 , ~t2 ) referenced with respect to the axis (A, ~t1 ) and the crack increment ∆a. As we can see on Figure 7, the new crack front M that is defined by these information does not lie on the shell surface,

175

ED

it is therefore necessary to project this new crack segment onto the shell midsurface. Since this element is non-planar, we have to deal with two issues: the normal vector to the shell midsurface varies inside this

PT

element and we need to pick one to define the projection plane ; the surface onto which we want to do the projection is non-planar.

The idea we suggest is to use the plane define by the entry side of the new element, that is the element

Figure 8, where the projection plane is named P = (N1 , ~e1 , ~e2 ) and the normal vector used for the projection ~n is the normal vector to the crack local frame.

AC

180

CE

side where the old crack tip is, and the next element side in the counter-clockwise direction. This is shown in

The next step consists in projecting a point P that belongs to the segment [AM ] onto plane P using

the normal vector ~n. This gives us a point P 0 thas is defined by its coordinates (α, β, 0) in the frame (N ; ~e1 , ~e2 , ~e3 ). It should be noted here that line (AP 0 ) cannot be used to define the new crack direction since it does not necessarily intersect one of the sides of the element. At this stage, we go back to the reference element using the Jacobian mapping to find the exit point of the crack inside the element. Since the reference

element is planar, finding the intersection is trivial. The image of point P 0 in the parent element is named

10

CR IP T

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

Figure 7: Crack growth direction and increment defining the new front M.

AC

Figure 8: Definition of the plane used to project the new crack increment.

0 Pref and is obtained thank to Equations (16) and (17).

Likewise, we can construct Aref

α αref −−−→ = −−→ . ||N1 N2 || ||I1 I2 ||

(16)

β βref (17) −−−→ = −−→ . ||N1 N4 || ||I1 I4 || which is the image of point A. Finally, we construct the half line

0 [Aref Pref ) and the new crack tip is the intersection between this half-line and one of the sides of the element.

11

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 9: Definition of the plane used to project the new crack increment.

Figure 10: Inverse mapping to the reference element and finding of the the intersection point between the new crack segment and the corresponding element side.

Then we use the Jacobian mapping to obtain the position of point B in the deformed configuration as the 185

image of point Bref in the parent configuration, see Figure 11. This projection algorithm slightly changes the crack increment, but as we use fine meshes, the curvature of the elements is small and consequently the 12

ACCEPTED MANUSCRIPT

ED

M

AN US

CR IP T

difference between the computed and the projected crack length is small.

PT

Figure 11: Mapping of the new crack segment from the reference to the current element.

3. Fracture model

metallic materials, we need a fracture criterion that has the following features : • mode I, II and III influence the crack evolution law,

AC

190

CE

As we want to model the propagation of cracks under impact loadings in shells mostly composed of

• propagation at low and high strain rates, that is tensile and shear driven fracture, • generalized plasticity.

To meet these requirements, we adopt a model based on the so called local approach to fracture proposed amongst others by Pineau [37]. We use a numerical adaption of this approach, similarly to many authors in the litterature, see for example Remmers et al. [21, 35], Haboussa et al. [39], Menouillard et al. [48]. The

13

ACCEPTED MANUSCRIPT

model is based on the calculation of an equivalent stress tensor at the crack tip : R ωσij (M ) dVM σ ˜ij = Ω R , ω dVM Ω

(18)

where σij are the components of the Cauchy stress tensor and ω is a weighting function than can take the r2

ω(r) = e−α R2 .

CR IP T

following form : (19)

The integration domain used to compute the equivalent stress tensor is a half-disk in 2D and a half-sphere 195

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. [58] for the tensile plane stress case and generalized to full 3D and shear driven fracture by Haboussa [59].

 σ ˜ ≥ σIc (ε) ˙

AN US

The propagation is first triggered by the following equations :   σ ˜ < σIc (ε) ˙ a˙ = 0 no propagation, a˙ 6= 0

(20)

propagation.

where σ ˜ is the maximum principal stress of the equivalent stress tensor and σIc is the analogous to the

M

fracture toughness and has to be identified.

Following Haboussa et al. [39], Schollmann et al. [58], we replace the stress intensity factors by the corresponding components of the equivalent stress tensor : → σ ˜22 , (21)

→ σ ˜12 , → σ ˜23 .

We also introduce the normalized equivalent stress components, following Haboussa [59] and Schollmann

CE

200

PT

ED

   KI    KII     K III

et al. [58] :

AC

σ ˜In =

|˜ σ22 | |˜ σ12 | |˜ σ23 | n n , rr˜ σII = , rr˜ σIII = . |˜ σ22 | + |˜ σ12 | + |˜ σ23 | |˜ σ22 | + |˜ σ12 | + |˜ σ23 | |˜ σ22 | + |˜ σ12 | + |˜ σ23 |

(22)

We consider an arbitrary 3D crack in an elastic medium and the asymptotic stress tensor around the

front is given by Equation (23) in a frame parametrized by two angles θ and ψ as shown in Figure 12. We 00

also consider plane stress conditions. The asymptotic stress tensor is given in the R frame that is obtain

205

0

by successive rotation of the crack front frame of the angle θ around n then of the angle ψ around t1 , as shown in Figure 12.

14

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 12: Various frames around the crack tip used to define the crack propagation law.

σ σ

00

11

=

σrr ,

22

=

2 cos(ψ) sin(ψ)σθz + σθθ cos2 (ψ),

33

=

−2 cos(ψ) sin(ψ)σθz + σθθ sin2 (ψ),

00

00

12

= σrz sin(ψ) + σrθ cos(ψ),

23

= −σθθ sin(ψ) cos(ψ) + σθz (2 cos2 (ψ) − 1),

13

= σrz cos(ψ) − σrθ sin(ψ).

PT

σ

00

M

σ

00

ED

σ

σ

00

(23)

CE

In the case of tensile driven fracture, the direction of propagation is the one that maximizes the tensile

AC

or hoop stress. This direction is obtained by solving the following equations :  ”   ∂σ22 = 0, ∂θ

”   ∂σ22

∂ψ

(24)

= 0.

The following semi-analytical solution is obtained : θctensile = 2sign(˜ σ12 ) arctan

ψctensile

   q 1 ˆ ˆ 2+8 S − (S) , 4

1 = sign(˜ σ23 ) arctan 2 15



2˜ σθz (θctensile ) σ ˜θθ (θctensile )



,

(25)

(26)

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 13: Transition between tensile and shear direction

with

n )p(ν) ˜III 1 √ 1+σ ˜In − (1 − σ and p(ν) = ( π − 5ν). Sˆ = n σ ˜II 4

(27)

In the case of shear driven fracture, the direction of propagation is the one that maximizes the Von Mises

ED

M

stress. This direction is obtained by solving the following equations  ”   ∂σV M = 0, ∂θ   ∂σV” M ∂ψ

(28)

= 0.

The following semi-analytical solution is obtained : π = sign(˜ σ12 ) 4

CE

PT

θcshear

ψcshear =

  8 1+ ν arctan 145

1 sign(˜ σ23 ) arctan 2



3 5

σ ˜In (1 − n σ ˜II

 −˜ σθθ (θcshear ) . 2˜ σθz (θcshear )

n σ ˜III 20

!!

,

(29)

(30)

The transition between shear and tensile fracture is based on an equivalent strain measure around the

AC

crack tip ε˜eq as presented in Haboussa et al. [39]. Two parameters are considered : ε˜shear and ε˜tensile . If

ε˜eq < ε˜tensile the crack will follow the direction of maximum tensile stress and if ε˜eq > ε˜shear , the crack will

210

follow the direction of maximum shear stress. A linear transition between the two values is considered, as shown in Figure 13. For the case of shells, we need two modifications of the method presented here. The first one is directly linked to the shell formulation. Indeed, the stress tensor at each Gauss point of the shell is given in the local frame that can be different in all points. An extra step is therefore needed in order to rotate the stress tensor

215

at each Gauss point in the same frame before computing its contribution to the equivalent stress tensor at 16

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 14: Equivalent stress area : mid-surface for the elastic case (top) and through the thickness integration for elastic-plastic case (bottom).

M

the crack tip. Secondly, when an elastic-plastic material is considered, our shell element has five integration points through the thickness. We adopt the following approach in this case : we select for the integration of

ED

the equivalent stress tensor all the Gauss point on the shell midsurface that belong to a half sphere centered at the crack tip, see Figure 14. We then propagate this information in all layers, and compute the equivalent stress tensor by taking into account the integration rule through the thickness :

PT

220

σ ˜ij =

R R h/2

ωσ (M ) dz dVM

ij Ω −h/2 R R h/2 ω Ω −h/2

dz dVM

,

(31)

CE

Finally, the crack velocity is obtained using the well known relation introduced by Kanninen and Popelar

AC

[60] adapted to be used with the equivalent stress tensor at the crack tip :   σIc (ε) ˙ a˙ = cR 1 − , σ e

(32)

where cR is the Rayleigh wave speed.

As we are using only Heaviside enrichment coupled with explicit time stepping, the relative time scale

between the computation in the bulk and for the crack have to be studied. The time step ∆t in the bulk is governed by the Courant condition and is in general very small as we use refined meshes. The crack 225

velocity is computed using Equation (32) is then multiplied by this time step to obtain the crack extension ∆a = a∆t ˙ . However, since the crack tip can only be on element boundaries because of the Heaviside 17

ACCEPTED MANUSCRIPT

function, the X-FEM approximation will produce a crack advance that can be significantly higher, which can be interpreted as artificial supersonic propagation, therefore unrealistic. To better understand this issue, we can perform a simple calculation considering the simulation results of the experiment by Kalthoff and 230

Winkler [61] reported for example in Belytschko et al. [15], Elguedj et al. [22], Song et al. [62]. The authors use a mesh of 80 × 80 elements with a characteristic length of ∆l = 1, 25mm. The time step used for such

simulation is about ∆t = 0, 1µs and the average crack velocity is a˙ = 1800m.s−1 which would result in an

CR IP T

average crack increment of ∆a = 0, 18mm. This value is roughly one tenth of the mesh size. Consequently a propagation that would be triggered in an element with the present method would result in a propagation 235

length of up to ten times the one obtained in the previous calculation, hence resulting in a crack velocity of up to 18000m.s−1 , largely above the Rayleigh wave speed of cR = 2120m.s−1 for this material.

We propose in this work to follow the same method as in Elguedj et al. [63] to ensure that the computed

AN US

crack velocity is physical and correct, and therefore impose crack advances that are compatible with the mesh size used in the numerical simulation. The idea in the proposed approach is to propagate the crack 240

only when the computed crack increment is close to the mesh size termed lcar in the remaining. Once the propagation criterion is met, the algorithm starts a step counter and computes the crack advance ∆a and compares it with the mesh size lcar . If the values are not close enough, the computation continues for the next time step. If the propagation criterion is met again at the next time step, a new tentative

∆a1 = 2∆ta. ˙ This value is again compared to lcar and the algorithm continues to retain the propagation

CE

PT

with ∆an = n∆ta˙ ≥ lcar .

ED

until the criterion is met as illustrated in Figure 15 where the propagation is performed after the n-th step

AC

245

M

crack extension is computed considering the crack velocity and twice the value of the mechanical time step

Figure 15: Crack update for a propagation consistent with the mesh size.

18

ACCEPTED MANUSCRIPT

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 time weighted average of meanT the equivalent stress tensor σ ˜ij given by : meanT σ ˜ij

=



1−

∆t Tcumul



meanT σ ˜ij +

∆t σ ˜ij , Tcumul

(33)

CR IP T

where Tcumul is a time counter that starts from the moment the criterion is met and as long as it stays valid, until lcar is achieved. This time counter is reset to zero after a propagation or when the criterion is 250

not met anymore. Finally this criterion has two geometrical parameters : R the size of the domain and lcar the distance above which we initiate propagation.

4. Numerical examples

AN US

The following examples have been chosen so as to illustrate several features of the proposed method. The first example is focused on the fracture model and demonstrates its ability to model tensile and shear 255

driven fracture as well as the transition from one mode to the other during the same simulation. The second example shows a simple example of mode I+III crack propagation to show the ability of the shell element to accurately describe such situations. The last two examples are devoted to comparisons of the

M

obtained simulations with complex experiments. In all the examples, but the second one, the material is considered elastic-plastic based on the level of stress and strain imposes by the external loading. As said 260

in the introduction, Generalized plasticity is considered in theses cases in order to be representative of the

ED

situation encountered in naval structures made out of high yield limit steel.

PT

4.1. Zhou-Rosakis-Ravichandran experiment

The Zhou-Rosakis-Ravichandran experiment [64, 65] is a variation of the well known experiment by Kalthoff and Winkler [61]. The specimen is not symmetrical and has only one crack. This modification allows to observe a transition between tensile and shear driven fracture in the same experiment,. For low

CE

265

impact velocities, the crack propagates in a tensile mode, and for high impact velocities, it propagates in a

AC

shear mode, as it was observed by Kalthoff and Winkler [61]. However, for intermediate velocities, the crack changes form tensile to shear driven fracture during its propagation. This example is particulary interesting to test the ability of the fracture model we use to represent the transition during the same simulation, as it

270

was shown in 2D under plane strain conditions by Haboussa et al. [39], Armero and Linder [66]. The geometry of the specimen as well as the mesh used for the simulation are shown in Figure 16(a). The specimen is a rectangular plate of width 2L and height 4L, the initial crack has a length of L = 50.8 mm, the plate’s width is e = 16 mm. The impact of the striker is modeled as an imposed velocity V0 . The mesh is composed of 203 × 100 X-Q4γ24 enriched shell elements. The material of the specimen is a maraging steel 19

ACCEPTED MANUSCRIPT

275

18Ni900 that is modeled by an elastic-plastic with isotropic hardening model. The material parameters are given in Table 2, the fracture parameters are given in Table 3 and are obtained using the same method as presented by Haboussa et al. [39]. The crack path obtained for three impact velocities V0 =20, 25, 30 m.s−1 are given on Figure 16(b) and Von Mises stress maps at various simulation times for V0 = 30 m.s−1 are shown in Figure 17.

280

We can see in Figure 16(b), that the proposed model can reproduced the different configurations observed

CR IP T

in the experiment. For a low impact velocity of 20m.s−1 , the propagation angle is roughly 70◦ as observed in the experiment and in other simulations [39, 66] for a tensile driven fracture. For high impact velocity of 30m.s−1 , a shear driven fracture is observed with an angle of around -10◦ , again as reported in the experiment [64, 65] and in the literature [39, 66]. For this case, Figure 17 shows Von Mises stress maps 285

at various simulation times. As can be observed in the images, the crack grows in the direction of the

AN US

maximum of the Von Mises stress which is the direction of the maximum shear stress. Finally, we can see in Figure 16(b) that for an intermediate impact velocity of 25m.s−1 , the crack starts to grow driven by the maximum shear stress and after a propagation of roughly half of the initial ligament, the crack bifurcates in the direction of the maximum tensile stress with a direction similar to the low velocity case. This behavior has been observed experimentally by Zhou et al. [64, 65]. It is important to note here that the transition is

CE

PT

ED

M

done automatically by the algorithm based on the threshold strain values.

AC

290

(a)

(b)

Figure 16: ZRR experiment : (a) geometry setup, boundary conditions and mesh used for the simulation, (b) obtained crack path for V0 =20, 25, 30 m.s−1 .

20

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

/ Pa

Figure 17: Simulation results showing Von Mises stress for the ZRR case with V0 =30 m.s−1 at times t=8, 16, 24, 32, 36, 40µs.

21

ACCEPTED MANUSCRIPT

4.2. Mode III case In order to asses the fracture model in mode III situations, we study an example similar to the one proposed by Li et al. [67]. A rectangular plate with an initial symmetrical notch along the largest side 295

is submitted to an out of plane constant velocity. The geometry and boundary conditions can be seen in Figure 18, the plate has a length of L = 5 m, height of H = 2 m and width e = 0.2 m, the notch has an

CR IP T

initial length l0 = 0.32 m. The imposed velocity is V0 = 5 m.s−1 and the mesh is composed of 126 × 51 XQ4γ24 elements. The material is considered linear elastic, the material and fracture parameters are given in Table 4. We do not consider an elastic-plastic material for this particular example, as its aim is to illustrate 300

the ability of the proposed method to predict mode III and mode I+III crack propagation.

We consider a relatively thick plate (L/e = 25 and H/e = 10) to ensure that transverse shear is significant and therefore a mode III loading of the crack. With the loading considered, we have observed that even for

AN US

higher slenderness values the plate has large out of plane displacements after a short time, which induces

CE

PT

ED

M

mostly membrane loading and consequently mode I at the crack tip.

Figure 18: Geometry, loading and boundary conditions for the mode III case.

The simulation results can be found in Figures 19 to 21. Figure 19 show the final crack path obtained

AC

305

at the end of the simulation. It can be seen that the crack propagates in the initial direction until roughly half of the specimen, and after a small perturbation symmetry is lost and the crack starts to zigzag. This is also observed in Figure 20 where we have plotted the normalized fracture mode contribution at the crack tip during the simulation. This corresponds to computing the appropriate components of the equivalent stress

310

tensor at the crack tip and calculating the respective contribution to the sum of those three components. The result is plotted as a total contribution in %, that is pure mode III corresponds to 100% for the mode III normalized component, etc. We can see at the very beginning of the simulation , a state of nearly pure 22

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 19: Crack path for the mode III case.

Figure 20: Normalized contribution of each fracture mode at the crack tip versus time for the mode III case.

23

ACCEPTED MANUSCRIPT

CR IP T

Uz Displacement / m

AN US

Time : 0.008s

Time : 0.024s

PT

ED

M

Time : 0.016s

AC

CE

Time : 0.032s

Time : 0.040s Figure 21: Deformed configuration and vertical displacement at various times for the mode III case. An amplification factor of 2 has been applied.

24

ACCEPTED MANUSCRIPT

crack element #

cut element #

input side #

αs

exit side #

βs

1

4

IV

0.75

II

0.25

2

5

IV

0.5

III

0.5

3

8

I

0.5

II

0.25

4

9

IV

0.75

III

0.75

CR IP T

Table 1: Needed parameter values to describe the crack shown in Figure 5.

E

200 GPa

Poisson ratio

ν

0.3

Density

ρ

7830 kg.m−3

Rayleigh wave speed

cR

2800 m.s−1

Yield stress

σy

200 MPa

Tangent modulus

Et

1600 MPa

AN US

Young’s modulus

Table 2: Material parameters of 18Ni900.

mode III is observed, producing a propagation in the direction of the initial crack. After 1ms, mode I and

315

M

mode III seem to oscillate while mode II is close to zero. According to Schollmann et al. [58], Haboussa [59] a state of mode I+III under traction produces a straight propagation, which is exactly what we obtain.

ED

After 10ms, a small perturbation in the crack path is observed, which results in a non zero contribution of mode II, then the mode III contribution is close to zero and the cracks continues to grow in a zigzag pattern. This behavior is also observed in Figure 21 where we have plotted at various times the vertical displacement

CE

4.3. ONERA experiment

In this section we consider the simulation of experimental results based on the device first presented in Potapov et al. [68] to test specimen made out of aluminium alloys. The test campaign has been conducted

AC

320

PT

component on the deformed configuration (amplified twice).

Critical stress

σIc

250 MPa

Tensile threshold strain

ε˜tensile

8.00×10−4

Shear threshold strain

ε˜shear

8.85 × 10−4

Characteristic mesh length

Lcar

0.001 m

Half disk radius

R

0.003 m

Table 3: Fracture model parameters of 18Ni900.

25

ACCEPTED MANUSCRIPT

E

200 GPa

Poisson ratio

ν

0.3

Density

ρ

8000 kg.m−3

Rayleigh wave speed

cR

2800 m.s−1

Critical stress

σIc

100 MPa

Characteristic mesh length

Lcar

0.04 m

Half disk radius

R

0.12 m

CR IP T

Young’s modulus

AC

CE

PT

ED

M

AN US

Table 4: Material and Fracture model parameters for the mode III case.

Figure 22: Drop tower experimental device : (a) schematic view, (b) picture of the various components if the setup.

on a drop tower of ONERA Lille, and consists in exhibiting a suddenly created crack growth in a piston-like vessel filled with fluid. The piston is impacted at its top by a dropped mass, which generates a sudden

325

rise of the internal pressure followed by deformation and rupture of the plate closing the piston’s bottom and leakage of the fluid from the piston. The main piece of the experimental setup is a thick-walled steel piston-type vessel filled with water at ambient temperature. A mass of several hundreds of kg is dropped 26

ACCEPTED MANUSCRIPT

from a pre-determined height impacting onto the piston. The impact velocity is in the range 1 − 10m.s−1 . The test rig is equipped with several measurement devices : high-speed cameras, HMI high-power cold330

light projectors, optical displacement transducers, strain gauges, piezoelectric pressure sensors and load cells, central computerized data acquisition system. Residual strains and deformed configuration after the experiment are obtained through a digitization device and software. The material we consider in this study is a marine construction steel with high yield limit, and its

335

CR IP T

behavior was characterized in Langrand et al. [69]. We consider an elastic-plastic with isotropic hardening behviour and neglect viscous effects. The parameters are given in Table 5 and the hardening curve in Figure 24. For confidentiality reasons, some of the values of the material parameters and simulation results have been normalized. Test samples were machined from a rectangular bloc so as to obtain a thiner circular or elliptic zone in the center of the plate. A hole and a crack was subsequently machined in the center of

ED

M

AN US

the shell as an initial defect.

AC

CE

PT

Figure 23: Geometry, boundary conditions and mesh used for the simulation of the EL30A6 experiment.

340

Young’s modulus

E

220 GPa

Poisson ratio

ν

0.3

Density

ρ

8000 kg.m−3

Rayleigh wave speed

cR

2800 m.s−1

Yield limit

σy

706 MPa

Table 5: Material parameters for the marine construction steel with high yield limit.

We are interested in this paper in studying the results of two experiments where an elliptical shell is used with a central hole and an inclined crack of two different initial length. The geometry and mesh used for the first specimen (EL30A6) can be seen in Figure 23. 27

CR IP T

ACCEPTED MANUSCRIPT

Figure 24: Hardening curve for the marine construction steel with high yield limit.

σIc

850 MPa

Characteristic mesh length

Lcar

0.0007 m

Half disk radius

R

0.0035 m

AN US

Critical stress

Table 6: Fracture model parameters for the EL30A6 experiment.

The fracture parameters used for the simulation are given in Table 6. As the fluid is not modeled in the

pressure obtained from the pressure sensor. The variation with respect to time of the experimental pressure

CE

PT

ED

as well as a simplified evolution used for the simulation are shown in Figure 25.

AC

345

M

simulation, the loading is applied through a normal pressure on the shell top surface using the experimental

Figure 25: Pressure loading for the EL30A6 experiment.

28

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 26: Experimental and numerical orthoradial strain at gage J1 versus time for the EL30A6 experiment.

Simulation results are shown in Figures 26 to 29. One can see that final crack paths are very similar between the experiment and the simulation : the cracks start to propagate along the initial direction and then

350

M

turn along the circumferential direction. The predicted final crack is slightly longer than the experimental one and propagates a little bit earlier based on what can be observed on the images acquired by the digital

ED

cameras and the sudden drop in strain measured by the gages. These differences are quite obvious when looking at the deformed configuration obtained in Figure 27 : when the crack propagates, the two parts of the specimen open and bend because of the maintained loading whereas in the experiment the same situation

355

PT

results in a local drop in the fluid pressure because of leakage. This shows that imposing the experimental pressure is note representative of the local loading when the crack grows and that a finer model would require to run a fluid-structure interaction simulation. A simple way to reduce the bending effect, as well

CE

as reduce the computational time, is to use mass scaling. This hypothesis can be accepted in such a case as the transient effects are not so fast in the considered case. Figure 28 presents the same results with an

360

AC

artificial increase of the density by a factor of 100 and a small reduction of σIc to 800 MPa. The final crack path is very similar to the previous case, and even closer to the experimental one as shown in Figure 29. we also compare strain evolution versus time at the position of gage J1 for the experiment and both simulations on Figure 26. One can see that without mass scaling the numerical result is closer to the experimental one before bending becomes important (arounds 1ms). For the case with mass scaling, the difference is more pronounced as expected, but the results are still of the same order of magnitude between the experiment 365

and the simulation.

29

ACCEPTED MANUSCRIPT

CR IP T

/ Pa

M

Time : 0.00084s

AN US

Time : 0.00076s

Time : 0.00100s

AC

CE

PT

ED

Time : 0.00092s

Time : 0.00108s Figure 27: Circumferential stress σθθ versus time and crack path on the undeformed (left) and deformed (right) configuration for the simulation of the EL30A6 experiment.

30

ACCEPTED MANUSCRIPT

CR IP T

/ Pa

M

Time : 0.00095s

AN US

Time : 0.00087s

Time : 0.00126s

AC

CE

PT

ED

Time : 0.00113s

Time : 0.00140s

Figure 28: Circumferential stress σθθ versus time and crack path on the undeformed (left) and deformed (right) configuration for the simulation of the EL30A6 experiment with mass scaling ρ = 100 × ρ0 .

31

CE

PT

ED

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(b)

AC

Figure 29: Comparison of numerical and experimental crack path for the EL30A6 experiment.

5. Blast tests

The last example consists in the simulation of blast experiments that were conducted jointly by DCNS

and DGA TN. The experimental setup consists in a circular plate of radius R =380 mm made out of the same material as in the previous case bolted in its circumference on a rigid frame; a spherical explosive 370

charge of a given mass m is suspended above the plate center at a given distance D and subsequently lit, see Figure 30. Because of the violence of the explosion, no experimental data can be acquired during the 32

ACCEPTED MANUSCRIPT

test and only post mortem crack path and deformation can be used for comparison. Several configurations with or without initial cracks of different length, with several explosive mass at various distances have been tested. We focus in this paper on the cases where a through the thickness initial straight crack is machined 375

on the plate with a length of L =40 mm, L =80 mm and L =160 mm, the same explosive mass m with a

AN US

CR IP T

distance D and 2D.

(a)

(b)

M

Figure 30: Blast test experimental setup : picture of the real setup (left) and simplified representation (right).

The pressure loading used in the simulation cannot be prescribed by analytical models such as CONWEP

ED

[70] as the explosive charge is very close to the sample and therefore the pressure wave cannot be considered as planar. Moreover, the CONWEP model does not take into account the effect of burnt gaz in the loading 380

which is important in the case considered. As the time scale of the shockwave and flow of burnt gaz is very

PT

small compared to the one of the deformation of the plate, we have chosen to use a decoupled approach. First, a simulation considering the plate as rigid is done in order to compute the equivalent pressure

CE

on the plate surface due to the explosion, and then this data is prescribed as the input pressure loading to simulate the plate deformation and the crack growth. The first simulation was performed using LS-Dyna 385

with a Multi-Material Arbitrary Lagrangian method, the burnt gases after detonation being modeled by the

AC

Jones Wilkins Lee state equation and considering that the air is a perfect gaz that undergoes an isentropic transformation. In this simulation, the plate is split in several concentric zones from which axisymmetric mean pressure values are extracted as a function of time, as presented in Figure 31 for an explosive charge of mass m at a distance D. We use 40 zones to apply the pressure loading that are numbered in a concentric

390

way starting form the center of the plate. The numerical simulation for the plate deformation and crack growth uses symmetry, and the geometry, mesh and pressure zones are given in Figure 32. The mesh is composed of X-Q4γ24 elements with a mean characteristic size of 5mm. This value is larger that the one used for the previous examples but is sufficient to obtain accurate and converged results as will be shown 33

ACCEPTED MANUSCRIPT

at the end of the paragraph . The material data is identical the previous example and is given in Table 5

AN US

CR IP T

and Figure 24, the fracture parameters are given in Table 7.

Figure 31: Mean normalized pressure versus normalized time for several zones for the blast experiment with an explosive mass

PT

38 0

crack

CE

Symmetric

m

m

ED

clamped

M

m and a distance D.

L/2

AC

395

Figure 32: Geometry, mesh and pressure zones for the blast test.

34

ACCEPTED MANUSCRIPT

Critical stress

σIc

910 MPa

Characteristic mesh length

Lcar

0.005 m

Half disk radius

R

0.01875 m

Table 7: Fracture model parameters for the blast test.

CR IP T

We have used the first test result with an initial crack of length L = 40 mm and a distance D between the plate and the explosive to adjust the fracture parameters used in the previous example so as to obtain the same final crack length of 180 mm between the numerical simulation and the experiment. A graphical comparison of the final deformed shaped is given in Figure 33. With these parameters, we have computed 400

the second case that has an initial crack of L = 80 mm, the results are given in Figure 34. In both cases, the crack nearly goes straight during the whole propagation, and we can observe a transverse opening of

AN US

the crack in the experimental case that is not seen in the numerical simulation. For the second case, the predicted crack extension is slightly longer than the experimental one : we obtain ∆L = 85 mm whereas ∆L = 55 mm can be measured form the post-mortem specimen. The third case considered has an initial 405

crack length of L = 160 mm with a distance between the plate and the explosive charge twice bigger than the two previous cases. The results can be seen in Figure 35. For this case the correlation between the final

M

crack length is very good between the experiment and the simulation as we obtain ∆L = 30 mm in both cases. The first test has also been simulated with a finer mesh with the same parameters. The results shown

410

ED

in Figure 36 exhibit similar crack paths and lengths although the elements are 2.5 times smaller.

6. Conclusion

PT

In this paper, we have presented a framework for the simulation of crack propagation under transient loadings in thick shells using the X-FEM and a crack propagation criterion that can represent both tensile

CE

and shear driven propagation. The underlying finite element is based on the enrichment of the existing Q4γ24 using Heaviside discontinuous functions. The propagation model is based on the local approach to 415

fracture and is the adaption to thin structures, i.e. plane stress, of the model introduced by Haboussa

AC

et al. [39], Schollmann et al. [58]. Thanks to the X-FEM, the crack does not have to conform to the shell’s midsurface finite element mesh and is represented explicitly with a separate 1D mesh. Through several examples, we have shown that the proposed model can represent in a single model tensile and shear driven fracture in shells as well as the transition form one regime to the other. Moreover, that propagation

420

criterion can accurately separate all fracture mode and is suitable for all mode combinations. Finally, we have simulated two complex experiments where the crack propagation is triggered by fast transient loadings and demonstrated that the model can reproduce accurately such complicated results. 35

CR IP T

ACCEPTED MANUSCRIPT

(a)

(b)

Figure 33: Blast test with an initial crack of L = 40 mm and distance D between the plate and the explosive charge, comparison

M

AN US

of experimental and numerical results.

(b)

ED

(a)

Figure 34: Blast test with an initial crack of L = 80 mm and distance D between the plate and the explosive charge, comparison

AC

CE

PT

of experimental and numerical results.

(a)

(b)

Figure 35: Blast test with an initial crack of L = 160 mm and distance 2D between the plate and the explosive charge, comparison of experimental and numerical results.

36

CR IP T

ACCEPTED MANUSCRIPT

Figure 36: Blast test with an initial crack of L = 40 mm and distance D between the plate and the explosive charge, comparison of numerical results with different mesh size.

7. Acknowledgments

425

AN US

Y. Jan was partially supported by a CIFRE fellowship of the French Association Nationale Recherche et Technologie and Naval-Group, formerly known as DCNS. Y. Jan, T. Elguedj and A. Combescure were also partially supported by a research contract with Naval-Group. This work was partially supported at Naval-Group by the French Direction G´en´erale de l’Armement trough the INFISCO program. The authors would also like to thank ONERA Lille and DGA TN for their support in obtaining the experimental results

430

M

shown in this paper. These supports are gratefully acknowledged.

References

ED

References

[1] T. Black, T. Belytschko, Elastic crack growth in finite elements with minimal remeshing, International Journal for

PT

Numerical Methods in Engineering 45 (1999) 601–620. [2] N. Mo¨ es, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal 435

for Numerical Methods in Engineering 46 (1999) 131–150.

CE

[3] D. Chopp, N. Sukumar, Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element/fast marching method, International Journal of Engineering Science 41 (2003) 845–869. [4] N. Sukumar, J.-H. Prevost, Modeling quasi-static crack growth with the extended finite element method. part i: Computer implementation, International Journal of Solids and Structures 40 (2003) 7513–7537. [5] N. Sukumar, D. Chopp, B. Moran, Extended finite element method and fast marching method for three-dimensionnal

AC

440

fatigue crack propagation, Engineering Fracture Mechanics 70 (2003) 29–48.

[6] E. Ferri´ e, J.-Y. Buffiere, W. Ludwig, A. Gravouil, L. Edwards, Fatigue crack propagation: In situ visualization using x-ray microtomography and 3d simulation using the extended finite element method, Acta Materiala 54 (2006) 1111–1122.

[7] J. Rannou, A. Gravouil, M. C. Baietto-Dubourg, A local multigrid x-fem strategy for 3d crack propagation, International 445

Journal for Numerical Methods in Engineering 77 (2008) 581–600. [8] J. Shi, D. Chopp, J. Lua, N. Sukumar, T. Belytschko, Abaqus implementation of extended finite element method using a level set representation for three-dimensional fatigue crack growth and life predictions, Engineering Fracture Mechanics 77 (2010) 2840–2863.

37

ACCEPTED MANUSCRIPT

[9] T. Fries, M. Baydoun, Crack propagation with the extended finite element method and a hybrid explicit–implicit crack 450

description, International Journal for Numerical Methods in Engineering 89 (2012) 1527–1558. [10] I. Singh, B. Mishra, S. Bhattacharya, R. Patil, The numerical simulation of fatigue crack growth using extended finite element method, International Journal of Fatigue 36 (2012) 109–119. [11] S. Kumar, A. Shedbale, I. Singh, B. Mishra, Elasto-plastic fatigue crack growth analysis of plane problems in the presence of flaws using xfem, Frontiers of Structural and Civil Engineering 9 (2015) 420–440. [12] S. Kumar, I. Singh, B. Mishra, A homogenized xfem approach to simulate fatigue crack growth problems, Computers and Structures 150 (2015) 1–22. [13] S. Kumar, I. Singh, B. Mishra, K. Sharma, I. Khan,

CR IP T

455

A homogenized multigrid xfem to predict the crack

growth behavior of ductile material in the presence of microstructural defects, https://doi.org/10.1016/j.engfracmech.2016.03.051 (2016). 460

Engineering Fracture Mechanics

[14] C. Duarte, O. Hamzeh, T. Liszka, W. Tworzydlo, A generalized finite element method for the simulation of threedimensional dynamic crack propagation, Computer Methods in Applied Mechanics and Engineering 190 (2000) 2227–2262. [15] T. Belytschko, H. Chen, J. Xu, G. Zi, Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous

AN US

enrichment, International Journal for Numerical Methods in Engineering 58 (2003) 1873–1905.

[16] T. Belytschko, H. Chen, Singular enrichment finite element method for elastodynamic crack propagation, International 465

Journal of Computational Methods 1 (2004) 1–15.

[17] J. R´ ethor´ e, A. Gravouil, A. Combescure, An energy-conserving scheme for dynamic crack growth using the extended finite element method, International Journal for Numerical Methods in Engineering 63 (2005) 631–659. [18] G. Zi, H. Chen, J. Xu, T. Belytschko, The extended finite element method for dynamic fractures, Shock and Vibration 12 (2005) 9–23.

[19] T. Menouillard, J. R´ ethor´ e, A. Combescure, H. Bung, Efficient explicit time stepping for the extended finite element

M

470

method, International Journal for Numerical Methods in Engineering 68 (2006) 911–938. [20] B. Prabel, A. Combescure, A. Gravouil, S. Marie, Level set X-FEM non matching meshes: Application to dynamic crack

ED

propagation in elastic-plastic media, International Journal for Numerical Methods in Engineering 69 (2007) 1553–1569. [21] J. Remmers, R. de Borst, A. Needleman, The simulation of dynamic crack propagation using the cohesive segments 475

method, Journal of the Mechanics and Physics of Solids 56 (2008) 70–92. [22] T. Elguedj, A. Gravouil, H. Maigre, An explicit dynamics extended finite element method. Part 1: mass lumping for

PT

arbitrary enrichment functions., Computer Methods in Applied Mechanics and Engineering 198 (2009) 2297–2317. [23] S. Kumar, I. Singh, B. Mishra, A. Singh, New enrichments in xfem to model dynamic crack response of 2-d elastic solids, International Journal of Impact Engineering 87 (2016) 198–211. [24] I. Asareh, Y.-C. Yoon, J.-H. Song, A numerical method for dynamic fracture using the extended finite element method

CE

480

with non-nodal enrichment parameters, International Journal of Impact Engineering 121 (2018) 63–76. [25] J.-L. Batoz, G. Dhatt, Mod´ elisation des structures par ´ el´ ements finis, Coques, volume 3, Herm` es, Paris, 1992.

AC

[26] Q. Zeng, A. Combescure,

A new one-point quadrature, general non-linear quadrilateral shell element with physical

stabilization, International Journal for Numerical Methods in Engineering 42 (1998) 1307–1338.

485

[27] J. M. Melenk, I. Babuska, The partition of unity finite element method : basic theory and applications, Computer methods in applied mechanics and engineering 139 (1996) 289–314.

[28] P. M. Areias, T. Belytschko, Non-linear analysis of shells with arbitrary evolvingcracks using xfem, International Journal for Numerical Methods in Engineering 62 (2005) 384–415. [29] J. Lasry, Calculs de plaques fissur´ ees en flexion avec la m´ ethode des ´ el´ ements finis ´ etendue (XFEM), Ph.D. thesis, INSA

490

de Toulouse, 2009. [30] R. Larsson, J. Mediavilla, M. Fagerstr¨ om, Dynamic fracture modeling in shell structures based on xfem, International

38

ACCEPTED MANUSCRIPT

Journal for Numerical Methods in Engineering 86 (2011) 499–527. [31] T. Chau-Dinh, G. Zi, P.-S. Lee, T. Rabczuk, J.-H. Song, Phantom-node method for shell models with arbitrary cracks, Computers and Structures 92-93 (2012) 242–256. 495

[32] J. Dolbow, N. Mo¨ es, T. Belytschko, Modeling fracture in mindlin-reissner plates with the extended finite element method, International Journal of Solids and Structures 37 (2000) 7161–7183. [33] R. O. Ritchie, J. F. Knott, J. Rice, On the relationship between critical tensile stress and fracture toughness in mild steel, Journal of the Mechanics and Physics of Solids 21 (1973) 395–410.

500

CR IP T

[34] M. Jirasek, Embedded crack models for concrete fracture, in: Computational Modelling of Concrete Structures, EURO C-98, volume 1, 1998, pp. 291–300.

[35] J. Remmers, R. de Borst, A. Needleman, A cohesive segments method for the simulation of crack growth, Computational Mechanics 31 (2003) 69–77.

[36] G. Wells, L. Sluys, A new method for modelling cohesive cracks using finite elements, International Journal for Numerical Methods in Engineering 50 (2001) 2667–2682. 505

[37] A. Pineau, Development of the local approach to fracture over the past 25 years: theory and application, International

AN US

Journal of Fracture 138 (2006) 139–166.

[38] J.-H. Song, T. Belytschko, Dynamic fracture of shells subjected to impulsive loads, Journal of Applied Mechanics 76 (2009) 051301.

[39] D. Haboussa, T. Elguedj, B. Lebl´ e, A. Combescure, Simulation of the shear-tensile mode transition on dynamic crack 510

propagations, International Journal of Fracture 178 (2012) 195–213.

[40] E. Dvorkin, K.-J. Bathe, A continuum mechanics four-node shell element for general non-linear analysis, Engineering Computations 1 (1984) 77–88.

M

[41] T. Belytschko, L. Bindeman, Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems, Computer Methods in Applied Mechanics and Engineering 88 (1991) 311–340. 515

[42] P. Areias, T. Belytschko, Analysis of fracture in thin shells by overlapping paired elements, Computer Methods in Applied

ED

Mechanics and Engineering 195 (2006) 5343–5360.

[43] J.-H. Song, P. Areias, T. Belytschko, A method for dynamic crack and shear band propagationwith phantom nodes, International Journal for Numerical Methods in Engineering 67 (2006) 868–893. [44] T. Menouillard, T. Belytschko,

International Journal for Numerical Methods in Engineering 84 (2010) 47–72.

PT

520

Smoothed nodal forces for improved dynamic crack propagation modeling in xfem,

[45] T. Menouillard, J.-H. Song, Q. Duan, T. Belytschko, Time dependent crack tip enrichment for dynamic crack propagation, International Journal of Fracture 162 (2010) 33–49.

CE

[46] D. Haboussa, D. Gr´ egoire, T. Elguedj, H. Maigre, A. Combescure, X-fem analysis of the effects of holes or other cracks on dynamic crack propagations, International Journal for Numerical Methods in Engineering 86 (2011) 618–636. 525

[47] T. Belytschko, N. Mo¨ es, S. Usui, C. Parimi, Arbitrary discontinuities in finite elements, International Journal for Numerical

AC

Methods in Engineering 50 (2001) 993–1013.

[48] T. Menouillard, J. R´ ethor´ e, N. Mo¨ es, A. Combescure, H. Bung, Mass lumping strategies for X-FEM explicit dynamics: Application to crack propagation, International Journal for Numerical Methods in Engineering 74 (2008) 447–474.

[49] T. Elguedj, A. Gravouil, A. Combescure, Appropriate extended functions for x-fem simulation of plastic fracture mechanics,

530

Computer Methods in Applied Mechanics and Engineering 195 (2006) 501–515.

[50] R. Pel´ ee de Saint Maurice, Extension de l’approche X-FEM en dynamique rapide pour la propagation tridimensionnelle de fissure dans des mat´ eriaux ductiles, Ph.D. thesis, Institut National des Sciences Appliqu´ ees de Lyon, 2014. [51] G. Ventura, On the elimination of quadrature subcells for discontinuous functions in the extended finite element method, International Journal for Numerical Methods in Engineering 66 (2006) 761–795.

39

ACCEPTED MANUSCRIPT

535

[52] A. Martin, J.-B. Esnault, P. Massin, About the use of standard integration schemes for x-fem in solid mechanics plasticity, Computer Methods in Applied Mechanics and Engineering 283 (2015) 551–572. [53] N. Mo¨ es, A. Gravouil, T. Belytschko, Non-planar 3D crack growth with the extended finite element and level sets - Part 1: Mechanical model, International Journal for Numerical Methods in Engineering 53 (2002) 2549–2568. [54] A. Gravouil, N. Mo¨ es, T. Belytschko, Non-planar 3D crack growth with the extended finite element and level sets - Part

540

2: Level set update, International Journal for Numerical Methods in Engineering 54 (2002) 2569–2586.

in Applied Mechanics and Engineering 195 (2006) 5343–5360.

CR IP T

[55] P. M. Areias, J. Song, T. Belytschko, Analysis of fracture in thin shells by overlapping paired elements, Computer Methods

[56] H. Bayesteh, S. Mohammadi, Xfem fracture analysis of shells: The effect of crack tip enrichments, Computational Materials Science 50 (2011) 2793–2813. 545

[57] N. Nguyen-Thanh, N. Valizadeh, M. Nguyen, H. Nguyen-Xuan, X. Zhuang, P. Areias, G. Zi, Y. Bazilevs, L. D. Lorenzis, T. Rabczuk, An extended isogeometric thin shell analysis based on kirchhoff–love theory, Computer Methods in Applied Mechanics and Engineering 284 (2015) 265–291.

[58] M. Schollmann, H. Richard, G. Kullmer, M. Fulland, A new criterion for the prediction of crack development in multiaxially

550

AN US

loaded structures, International Journal of Fracture 117 (2002) 129–141.

[59] D. Haboussa, Mod´ elisation de la transition traction-cisaillement des m´ etaux sous chocs par la X-FEM, Ph.D. thesis, Institut National des Sciences Appliqu´ ees de Lyon, 2012.

[60] M. Kanninen, C. Popelar, Advanced Fracture Mechanics, Oxford University Press, Oxford, UK, 1985. [61] J. Kalthoff, S. Winkler, Failure mode transition at high rates of shear loading, DGM Informationsgesellschaft mbH, Impact Loading and Dynamic Behavior of Materials 1 (1988) 185–195. 555

[62] J.-H. Song, P. Areias, T. Belytschko, A method for dynamic crack and shear band propagation with phantom nodes,

M

International Journal for Numerical Methods in Engineering 67 (2006) 868–893.

[63] T. Elguedj, R. Pel´ ee de Saint Maurice, A. Combescure, V. Faucher, B. Prabel, Extended finite element modeling of 3d dynamic crack growth under impact loading, Finite Element in Analysis and Design accepted. (2018).

560

ED

[64] M. Zhou, A. Rosakis, G. Ravichandran, Dynamically propagating shearbands in impact-loaded prenoteched plates I. Experimental investigations of temperature signatures and propagation speed, Journal of the Mechanics and Physics of Solids 44 (1996) 981–1006.

[65] M. Zhou, A. Rosakis, G. Ravichandran, Dynamically propagating shearbands in impact-loaded prenoteched plates II.

PT

Numerical simulations, Journal of the Mechanics and Physics of Solids 44 (1996) 1007–1032. [66] F. Armero, C. Linder, Numerical simulation of dynamic fracture using finite elements with embedded discontinuities, 565

International Journal of Fracture 160 (2009) 119–141.

CE

[67] T. Li, J. Marigo, D. Guilbaud, S. Potapov, Variational approach to dynamic brittle fracture via gradient damage models, Applied Mechanics and Materials 784 (2015) 323–341. [68] S. Potapov, B. Maurel, A. Combescure, J. Fabis, Modeling accidental-type fluid–structure interaction problems with the

570

AC

sph method, Computers and Structures 87 (2009) 721–734.

[69] B. Langrand, N. Lecomte, A. Menegazzi, T. Millot, Submarine hull integrity under blast loading, International Journal of Impact Engineering 36 (2009) 1070–1078.

[70] C. N. Kingery, G. Bulmash, Air blast parameters from TNT spherical air burst and hemispherical surface burst, Technical Report, U.S. Army Ballistic Research Laboratory, 1984.

40