Numerical modeling of crack propagation with dynamic insertion of cohesive elements

Engineering Fracture Mechanics xxx (xxxx) xxxx

Contents lists available at ScienceDirect

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

Numerical modeling of crack propagation with dynamic insertion of cohesive elements ⁎

D. Uribe-Suáreza,b, , P.-O. Boucharda, M. Delbob, D. Pino-Muñoza a MINES ParisTech, PSL Research University, CEMEF-Centre de Mise en Forme des Matériaux, CNRS UMR 7635, CS 10207 rue Claude Daunesse, 06904 Sophia Antipolis Cedex, France b Université Côte d’Azur, Observatoire de la Cote d’Azur, CNRS, Laboratoire Lagrange, CS 34229, 06304 NICE Cedex 4, France

A R T IC LE I N F O

ABS TRA CT

Keywords: Crack propagation Cohesive zone models Mesh-independent cohesive elements Dynamic insertion Remeshing techniques

One of the most challenging issues in computational fracture mechanics is the propagation of a crack through a finite element mesh for arbitrary crack paths. In this work, this problem is approached by means of an advanced remeshing technique that propagates a crack using cohesive elements. The crack direction is computed using the maximal energy release rate criterion which is implemented using finite elements and the Gθ method. The remeshing procedure used here is composed of two stages. In the first step, a conforming mesh is obtained in the computed crack direction, ensuring that edges are placed over the sought direction. In the second stage, cohesive elements are dynamically inserted at the conforming edges previously remeshed. The combination of this remeshing technique with dynamic insertion of cohesive elements, leads to a mesh-independent crack propagation method. The effects of different numerical and physical parameters regarding the crack path and fracture energy is investigated.

1. Introduction The identification of the most likely mode of failure and the application of a suitable failure criterion are crucial in the design of structures and machine components. In addition, the modeling of crack growth plays a essential role in the assessment of engineering structures regarding more accurate prediction of structural damage and failure can prevent catastrophic failures. The fracture mechanics community recognises the role of inherent structural flaws that affect performance and life. These inevitable defects which combined with structural cracking result in high stress concentrations can lead to failure of engineering structures and threaten safety [1,2]. Fracture can be investigated using laboratory experiments, analytical studies and numerical simulations. The first option is usually expensive, but often required to some extent, whereas the second option is restricted to simple configurations [3]. Therefore, numerical simulations are often an effective strategy to deal with most complex cases. Different applications in civil engineering, aerospace or automotive industries require accurate prediction of crack propagation. A suited and proven approach for the study of fracture mechanics is the finite element method. However arbitrary crack propagation through a finite element mesh is a challenging and demanding task [4]. The fracture process can be decomposed into 2 steps: crack initiation and crack propagation. The first one is essential but not easy. If the goal is to describe this process for a body without an initial pre-crack, damage-based numerical models should be used. These models study the evolution of damage in a continuous way and, at a critical damage value, the crack is initiated as a result of a nucleation process [5–7]. In the fracture mechanics field, the simulation of crack initiation is a challenging task and for this reason, this process often is avoided. Efforts are oriented towards the

⁎

Corresponding author. E-mail address: [email protected] (D. Uribe-Suárez).

https://doi.org/10.1016/j.engfracmech.2020.106918 Received 30 July 2019; Received in revised form 11 December 2019; Accepted 28 January 2020 0013-7944/ © 2020 Elsevier Ltd. All rights reserved.

Please cite this article as: D. Uribe-Suárez, et al., Engineering Fracture Mechanics, https://doi.org/10.1016/j.engfracmech.2020.106918

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Nomenclature

δ δc δmax ν σc σij θ a E G Gc KI , KII p

t tmax Ui Vi vi CZMs DOF FE LEFM MERRC PZL SDA XFEM

effective opening displacement characteristic opening displacement maximum effective opening displacement Poisson’s ratio maximum cohesive normal traction stress tensor virtual direction of propagation crack length Young’s modulus energy release rate critical energy release rate stress intensity factors for mode I and mode II pressure

cohesive traction maximum cohesive traction displacement virtual displacement field velocity cohesive zone models degree of freedom finite element linear elastic fracture mechanics maximum energy release rate criterion process zone length Strong Discontinuity Approach Extended Finite Element Method

study of the evolution of a pre-existing crack [8]. When a pre-existing crack is present, its propagation introduces a displacement discontinuity. Dealing with this discontinuity in a finite element mesh is fundamental in fracture mechanics. For this reason, many approaches have been developed to handle it. Proposed approaches range from simple procedures like element erosion [9], where elements are removed once their load carrying capacity has been eroded, to more complex methodologies using enriched finite element methods and/or remeshing operations. Element erosion is perhaps the simplest way to take into account discontinuities in finite element simulations due to its low computational cost and its straightforward implementation. However, this method also presents several limitations, such as volume loss and mesh size and element shape sensitivity issues [10]. In the case of crack propagation without any remeshing, an interesting methodology developed by Belytschko [11], known as the element free Galerkin method, only needs nodal information and boundary descriptions. In this method, when evaluating the displacement in one point, any node that is not visible from this point is omitted (visibility criterion). A straightforward implementation of this criterion results in some interior discontinuities around the crack tip. One of the main drawbacks of this method is that the boundary conditions are difficult to impose [12]. Among methods without remeshing, there is another approach called Arbitrary Local Mesh Replacement method. Introduced by Rashid [13], this method uses two different meshes: one surrounding the domain close to the crack tip, and that moves with it, and an other one that describes the entire domain. In this method, additional degrees of freedom and a special treatment of partial elements are required. Another model is the Strong Discontinuity Approach (SDA) [14–16] in which displacement jumps due to the presence of the crack are embedded locally in each cracked finite element without affecting neighbouring elements. In this method, additional degrees of freedom must be added to the finite element model. Another famous method for crack propagation without remeshing is the well-known Extended Finite Element Method (XFEM) [17,18] in which the displacement-based approximation is enriched near a crack by incorporating both discontinuous and near tip asymptotic fields. In this method, discontinuities might be modeled independently from the mesh [10]. As a disadvantage, some burdensome modifications to the finite element code are required. Furthermore, this formulation might not be suitable in the case of large strain sollicitations. One of the most recently developed methods, among techniques without remeshing, is the phase field approach [19–21]. It is a numerical technique to deal with discontinuities based on energy minimization principles. The phase field method has the notable advantage of avoiding explicit front tracking by making material interfaces spatially diffuse. It does not require any prescription of shape geometry and allows crack nucleation and branching [22]. Among crack propagation with remeshing operations, the approach proposed by Bouchard et al. [8] is worth mentioning. Through the use of an advanced remeshing technique combined with nodal relaxation propagation of the crack is achieved. An other common alternative is the family of cohesive zone models (CZMs) [23,24]. This approach was developed in order to model the energy dissipation rate, an issue not tackled by the aforementioned methods. Through CZMs the fracture process is modelled as the transition from a sound material to a fully broken one. CZMs are independent from the mechanical behavior of the bulk material, the extent of the cracks and the size of the plastic zone [25]. They also offer other advantages, such as straightforward implementation inside conventional finite element codes. Even so, they can sometimes introduce artificial stiffness into the system [26]. Combinations of remeshing operations and cohesive zone models have been performed in the literature. Chiaruttini et al. [27] present an approach for the numerical simulation of crack propagation, based on cohesive models, in the case of structures subject to mixed mode loadings. In this work, the evolving crack path is remeshed as the crack evolves. This approach is able to tackle arbitrary crack paths. Another recent approach is presented by Choi and Park [28], in which, through a novel stress recovery technique, a domain integral and an element splits remeshing procedure, an accurate crack path under mixed-mode as well as a mesh bias reduction are obtained. In this method, a virtual mesh around the crack tip is generated in order to get an accurate stress evaluation to compute the crack propagation direction. Most CZMs approaches suffer from crack-path mesh dependency [10]. The work presented here will handle the propagation of a crack through the use of remeshing operations together with the dynamic insertion of cohesive elements in a mesh-independent way. The novelty of the present work relies on the dynamic insertion of cohesive elements in an arbitrary direction (mesh-independent). In addition, the cohesive elements, that are characterized by a phenomenological cohesive 2

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

law, ensure that the released energy due to the fracture process is controlled. Also, the impact of different numerical parameters on the dynamic insertion of cohesive elements process is studied. The rest of the paper is organized as follows. In Section 2, the remeshing operations and dynamic insertion of cohesive elements in a mesh-independent way are explained. This section also briefly reviews crack propagation criteria using an energetic approach. One benchmark example is presented in Section 3 to show the accuracy of the proposed method when calculating the crack propagation direction under mixed-mode loading. In Section 4, the impact of different numerical and physical parameters on the crack path and fracture energy is studied. Finally some concluding remarks are presented in Section 5. 2. Description of the approach In this section some of the fundamental concepts of cohesive zone models are presented. The method used to perform crack propagation, using remeshing operations and dynamic insertion of cohesive elements in a mesh-independent way is explained. An energetic criterion commonly used to calculate the crack propagation direction will also be described. 2.1. Mesh independent cohesive zone models A pioneering perspective stating that fracture is a phenomenon taking place accross an extended crack tip (cohesive zone) was proposed back in the 60’s [23,24]. The proposed concept generated cohesive zone models (CZMs). This simple approach asserts that, at the crack tip, it exists a region of finite size where there is a transition from a fully broken material to a sound material [29]. Fig. 1 shows a schematic representation of the cohesive region (process zone), where thr fracturing process is taking place in a brittle material. The cohesive zone corresponds to prospective fracture surfaces ahead of a crack that are permitted to separate under loading. Atomic or molecular forces are in charge of preventing the process of separation and creation of crack surfaces [30]. The force that exerts a resistance to the opening of new surfaces is known as cohesive force, and it is modeled through a phenomenological traction-separation law (cohesive law). Many cohesive laws might be used to model the fracture process. Fig. 2 depicts the irreversible exponential cohesive law used in this work [31] and described by Eq. (1). δ

t=

⎧ eσc δ e− δc δ c

⎨ ⎩

tmax δ δmax

(Loading ) (Unloading )

(1)

where t is the cohesive traction, δ is the effective opening displacement, σc is the maximum cohesive normal traction and δc is the characteristic opening displacement. The present work does not account for plastic deformation at the crack tip. Even though, when working in elasticity, the final crack path will be different depending on the loading path applied, this analysis is out of the scope of the present work, and will not

Fig. 1. Schematic representation of the cohesive zone. Red arrows represent the distribution of traction loading over the cohesive region. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 3

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Fig. 2. Cohesive law in terms of an effective opening displacement δ and traction t (loading-unloading envelop).

be addressed here. Readers interested in the effects caused by different loading paths can refer to Yi et al. [32]. In this work, a fracture analysis dedicated to steel pipelines, used for oil and gas transportation, with embedded cracks and subjected to biaxial loading conditions is performed. The critical energy release rate Gc , i.e. the fracture energy required to create the new free surfaces and break the material can be related to the cohesive law in the following way:

Gc =

∫0

∞

tdδ

(2)

where t is a scalar effective traction and δ an effective opening displacement. The total amount of energy dissipated by the cohesive elements per unit of length can be computed simply using (2). Thus, the degradation process can be fine tuned by picking the mathematical form of t. The control of the fracture energy in this method is one of the main reasons the use of the cohesive zone models (cohesive elements) was chosen for this work. When using cohesive elements, a common issue is that an artificial reduction of the stiffness of the material can be induced. This is due to the fact that, in most traction-separation laws, there is a initial region where the traction increases monotonically from zero up

Fig. 3. (a) Prescribed crack path (dashed blue line). (b) Insertion process of cohesive elements previously used in the literature. (c) Insertion process of cohesive elements used in this work. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 4

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

to a maximum value. This increase can be linear or not depending on the traction-separation law used. When the traction level increases as a function of the opening displacement, an artificial stiffness that modifies the macroscopic response of the material is introduced into the system [26]. This problem can be solved through the introduction of Lagrange multipliers in such a way that the opening of the element is only allowed once a critical traction is achieved [33]. This solution has the disadvantage that it requires the modification of the finite element formulation. Commonly, cohesive zone models have been used to solve problems involving interfaces, surfaces undergoing decohesion or problems where the crack path is known a priori [34,35]. In these cases, cohesive elements can be inserted in a limited zone (known crack path) in the finite element mesh, reducing, but not completely solving the problem regarding artificial reduction of the stiffness. In the case where the crack path is unknown, mesh dependency may appear as a consequence of the fact that cracks are only able to propagate across boundaries between bulk elements. Therefore, the crack path depends on the mesh. This issue is depicted in Fig. 3b, where the predicted crack path and the actual one are shown in blue dashed line and green line respectively. Once the direction of propagation (blue dashed line in Fig. 3a) is computed using a suitable criteria, Fig. 3b shows the procedure that has been widely used to insert cohesive elements. These are inserted through the closest edges and nodes to the computed direction. The insertion of cohesive elements in this way exhibits a mesh dependency behaviour, so if the mesh changes, the crack pattern changes slightly as well [10]. When the crack path is unknown, some authors choose to insert cohesive elements on each interface between bulk elements throughout the material [36,37]. This approach increases the number of degrees of freedom greatly and makes the predicted crack path mesh-dependent [27,38]. Fig. 3c shows the methodology used here. Before inserting the cohesive elements, a local remeshing procedure in the predicted direction is performed [39]. Contrary to the process shown in Fig. 3b, the one shown here is mesh-independent. Inserting cohesive elements on the fly, i.e., while the crack tip is propagating through the domain, avoids a well-known strong drawback in the implementation of cohesive zone models into a FE framework: artificial reduction of the stiffness of the material. Other implications regarding how cohesive elements can be used to model fracture are discussed in [10]. Once the remeshing has been performed, the insertion process of cohesive elements can take place. A graphical illustration of the insertion process implemented in this work is shown in Fig. 4. Pragmatically, two neighbor bulk elements are separated by inserting a cohesive element at their shared face. This insertion is carried out by duplicating the nodes that form the separating face and inserting new cohesive elements linking the original nodes to the new duplicated ones (Fig. 4c). In Fig. 4a the red dashed line showing the faces that will be split can be observed. Fig. 4b shows the blue dots corresponding to the nodes that have been duplicated (in this case there are two nodes at the same location). The red line corresponds to a initially flat cohesive element. Finally, in Fig. 4c, it is shown that, after loading, the new inserted cohesive elements open. Fig. 5a shows the crack tip position at time step n. Also, a red line depicts the computed crack propagation direction. Fig. 5b shows the crack tip position after propagation at time step n + 1. It also shows a dashed magenta line illustrating the performed remeshing operation and the previous crack tip position. Finally, Fig. 5c shows a zoomed-in view of the inserted cohesive elements (red) along the crack path.

Fig. 4. Schematic representation of the insertion of cohesive elements. 5

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Fig. 5. (a) Crack tip position at time step n. (b) Crack tip position at time step n + 1. (c) Inserted cohesive elements in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2.2. Crack propagation Once a crack has initiated, it is necessary to check whether or not the crack is going to propagate and in which direction. This work will use a criterion based on the distribution of the energy throughout the cracked part: the maximum energy release rate criterion [40].

2.2.1. Maximum energy release rate criterion (MERRC) According to this criterion, the crack propagation direction will be the one which maximises the energy release rate (G), i.e., the energy required to increase the crack length by one. It is evaluated using all virtual and kinematically admissible displacement fields at the crack tip neighborhood. Here, for the computation of the energy release rate (G), the numerical technique known as Gθ method [41] will be used. The Gθ method is based on a virtual displacement field V representing the virtual kinematics of the crack. According to Bouchard et al. [8], the Gθ method is very accurate and completely mesh independent. Additionally, its implementation is quite simple and multiple extensions are available. It should be noted that the energy release rate can also be defined as the decrease in the total potential energy during crack growth. Taking into account that the well-known path-independent parameter J-Integral is equivalent to the rate of change of the total potential energy with respect to the crack length, it is possible to state that, for elastic bodies, the J-integral is equivalent to the energy release rate (Gθ method) [42,43]. In this implementation of the Gθ method, we first defined two contours C1 and C2 around the crack tip; those contours divide this region into three domains Cint , Cring and Cext as it is shown in Fig. 6. The virtual displacement field, V(v1, v2) , is defined by Eq. (3).

Fig. 6. Contours and domains used to computed G using the Gθ method. 6

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

V=

( ) cos(θ) = (1 − ) sin(θ)

⎧ v1 = 1 −

AB AC

⎨v 2 ⎩

AB AC

(3)

where O is the crack tip, B is an integration point belonging to the ring, A and C are the intersections between OB and C1 and C2 respectively; and θ is the virtual direction of propagation. Hence the values of V in the three domains are:

• The norm of V in C • The norm of V in C • The norm of V in C

is constant and equal to 1. is equal to 0. ring varies continuously from 1 to 0. int

ext

When there is neither thermal strain nor load applied directly to the crack faces, the energy release rate may be expressed as:

G=

∫ring ⎛⎝σij Uj,k Vk,i − 12 σij Uj,i Vk,k ⎞⎠ dC

(4)

In Eq. (4), described in [8] Einstein’s notation is used, σij is the stress field, Ui is the displacement field and Vi is the virtual displacement field. Furthermore, as Vi varies only inside Cring , the integration will only be performed over Cring (Fig. 7a). The energy release rate is calculated for the ring of elements for a given angle θ belonging to [−70°, 70°] [44]. The Gθ method provides an efficient way of computing G for a virtual extension in the direction θ . Fig. 7a shows the virtual ring in which all the integration points are considered. G can be computed for each value of θ which makes the identification of the θ value that maximises G (θ) straightforward (Fig. 7b). The direction of propagation will be the value of θ that corresponds to the maximum value of G, i.e., the maximum of the curve G (θ) in Fig. 7b. 3. Numerical examples In this section, different benchmark problems are considered to assess the proposed methodology. The present work has been implemented in CimLib, a C++ in-house finite element library developed at the CEMEF [45]. The finite element framework in CimLib is an implicit formulation (mixed) with first-order elements using linear interpolation functions for both velocity and pressure. Both variables are P1. The P1/ P1 elements (linear in both → v and p) do not satisfy the inf − sup condition [46,47]. For this reason an extra degree of freedom (DOF) in the center of each element is added in → v . This extra DOF is condensed and thus written as a function of the remaining DOFs of the element. The use of this extra DOF is a stabilization technique often called “bubble” or P1 + / P1 which satisfies the inf − sup condition, and has been implemented in the CimLib [48]. In the work presented here, isotropic unstructured triangular meshes that are refined in the neighborhood of the crack tip are used. All the examples presented here consider an initial notch or a prescribed crack, nucleation is not treated. The simulations performed are quasi-static. Also, a single crack is considered. Although the crack propagation distance is indeed an important parameter involved in fracture process simulations, here, for the sake of simplicity, the crack growth distance is set to a fixed value equal to 8 times the average size of the elements attached to the crack tip. This length size remains small considering the process zone length G E (PZL = c 2 , see Section 4). The influence of this parameter will not be studied in this article. A first example regarding crack πσc

propagation under mixed-mode loading is examined to show the accuracy of the proposed method in terms of propagation direction. Another example that includes several holes, showing the influence of these on the crack path is subsequently addressed. 3.1. Crack propagation of an edge-crack under mixed-mode loading An edge crack geometry fixed at the bottom and subjected to a top unit shear load is considered here (Fig. 8). Material properties are also given in Fig. 8. This problem was also solved by Nguyen-Xuan et al. [49] and Liao et al. [3]. In order to prove that the computation of the crack propagation direction through Gθ method implemented here is accurate, the simulated crack propagation

Fig. 7. (a) Ring of elements around the crack tip. (b) G (θ ) curve for the maximum energy release rate criterion. 7

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Fig. 8. Geometry of the plate with an edge crack under shear (dimensions in cm).

direction was compared against the direction value computed using the theoretical crack growth methodology based on the local stress fields at the crack tip (LEFM). Therefore, knowing the values of KI , KII and using Eq. (A.4, See Appendix A), it is possible to compute the theoretical direction of propagation of a crack in a mixed-mode problem. KI and KII could be computed based on displacement or stress extrapolation techniques. However, the accuracy of such methods in a finite element framework requires high accuracy of the stress field of the crack tip. The use of quarter-point elements [50] is recommended for better accuracy. Energetic approaches such as the Gθ method are less sensitive to the accuracy of the stress singularity at the crack tip and can, therefore, be used with conventional elements. Using the exact values of the stress intensity factors (Fig. 8) and Eq. (A.4), it is possible to compute that the exact direction of propagation is θ = −14.74°. This simulation is carried out with the aforementioned crack propagation technique and using a θ step variation of 0.01° in the Gθ method. It is found that the angle of propagation is approximately θ = −14.92°. Even though this value differs slightly (0. 18° ) from the theoretical one, this is proof that the Gθ method is able to compute a quite accurate crack propagation direction. To calculate the aforementioned value using the approach presented here, it was only taken the direction calculated at the first time step once the crack has reached the threshold that allows its propagation. Fig. 9 shows the comparison in terms of crack path between 2 previous works and the approach presented here. It is shown that this work reproduces the crack path very well which means that, compared to the results presented in [49,3], a good agreement is obtained here. 3.2. Cracked beam with three holes This case consists in a cracked beam supported in two points and loaded in the center, as shown in Fig. 10. The material properties are: Young’s modulus E = 29. 106 psi , Poisson’s ratio ν = 0.3. The load P is equal to 1 lb . Here, three cases of the initial crack length a and its distance b from the left side of the beam are considered. The different values considered for these parameters are given in

Fig. 9. Comparison of crack path of the edge-cracked plate obtained here against different literature results. 8

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Fig. 10. Cracked beam supported in two points and loaded in the center. The beam has three holes to create complex crack paths. (units in inches). Table 1 Geometrical configurations.

Case I Case II Case III

a (inches)

b (inches)

1.0 1.5 1.5

4.0 5.0 6.0

Fig. 11. Comparison of the crack path between present work and Nguyen-Xuan et al. [49] for (a) Case I, (b) Case II and (c) Case III. It is shown only a section of the beam in the region. near the holes. 9

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Table 1. In practical applications holes are common an can affect crack propagation. The proposed benchmark example aims at dealing with mixed-mode crack paths that are affected by the presence of holes. This problem has also been addressed by [49,51,52], both numerically and experimentally. Fig. 11a shows the complete crack path for Case I. The crack reaches the second hole from the lower left side. Fig. 11b shows the complete crack path for Case II. Crack propagation takes place through the material region between the first and second hole. The crack reaches the second hole from the lower right side. Finally, Fig. 11c shows the complete crack path for Case III. In this situation the crack reaches the first hole from the lower right side. In all the three cases the comparison between this work and Nguyen-Xuan et al. [49] is shown. The crack paths predicted in this work show excellent agreement with cited work [53,51,52]. It is clear that, due to the presence of holes, the crack path experiences significant changes. The holes attract the crack path depending on the initial notch geometry. These two mixed mode examples validate this new crack propagation method that combines (i) crack propagation direction (Gθ method), (ii) automatic remeshing and (iii) dynamic insertion of cohesive elements on the crack path. In the next section we address a topic which is rarely studied in the literature: the influence of time step on the numerical solution. A sensitivity analysis to the cohesive law parameters is also presented and discussed.

4. Influence of numerical and physical parameters In this section the influence of different numerical and physical parameters on the solution is studied. This influence is assessed in terms of either the force-displacement curve or the crack path comparison. The benchmark example considered in Section 3.1 (Fig. 8) is used again. The sensitivity to the mesh size and the time step are considered. First, aiming to prove the convergence of the solution in terms of the force-displacement curve as a function of the element size of the mesh, three different meshes are considered. Table 2 shows the element sizes of the different meshes used here. Simulations using these meshes are performed. Even though, in the present work, the mesh is refined everywhere, it is worth mentioning that it is also possible to only remesh the neighborhood of the crack tip while it evolves. Fig. 12 shows the force-displacement curve for a case where the mesh size changes. The red “ ” indicates the start of the propagation process. Mesh convergence in the solution is obtained. The curve up to the red “ ”, i.e. when propagation process has not initiated, exhibits a totally linear relation between force and displacement. This behaviour is linear because cohesive elements have not been yet inserted into the mesh. Once they are inserted, the nonlinear behavior appears. Due to the presence of oscillations after the propagation process has started, the nonlinearity in the force-displacement curve after this point looks smoother. The above Figure exhibits oscillations after the crack starts its propagation. This is due to the fact that, in this work, once the energy released rate (G) exceeds its critical value (Gc , Eq. (2)), the propagation of the crack is allowed. After propagation, cohesive elements are inserted leading to a decrease of the stresses around the crack tip, leading to a drop in G. This drop of G eventually leads to a value lower than Gc and, thus, propagation is momentaneously arrested. Afterwards, an increase of the loading brings G to a value higher than Gc , the propagation then restarts. This is an essential point (rarely mentioned in the literature) which is directly related to the space and time discretizations associated with the FE method. The curve describing the evolution of G through time for the three simulated cases (Reference Mesh, Mesh 1 and Mesh 2) is presented in Fig. 13. The red “ ” indicate the start of the propagation process which are the same for all cases. This figure shows that G varies from values close to 0.25Gc up to values greater than 2Gc . This range of variation is thought to depend strongly on the time step. Therefore, the time step is the next numerical parameter studied. Taking into account that mesh convergence is achieved, from now on, and for the sake of simplicity, the “Reference Mesh” case will be used as a reference for the upcoming comparisons. It will be referred to as “Reference”. Fig. 14 shows the force-displacement curve obtained when the time step was decreased by half. In this case, the necessary force to start the propagation process ( ) is less than the one required in the ”Reference” case. In addition, when half of the time step is used, crack propagation starts sooner than in the ”Reference” case. Also, after crack propagation has started, for the same time, using half of the time step, the crack tip has propagated a greater distance than in the “Reference” case, so less force is needed to continue the propagation process. The evolution of G through time is presented in Fig. 15. As a result of using half of the time step value, the crack propagation process begins sooner and leads to faster propagation of the crack. It is also evident that decreasing the time step reduces the oscillation range of G. Thus, the crack propagation process has a strong dependence on the time step. This phenomenon is due to the fact that the same propagation distance is prescribed once G exceeds Gc , regardless of the time step used. In other words, this is due to Table 2 Elements size of the different meshes used (mm).

Reference Mesh Mesh 1 Mesh 2

Far away from the crack tip

Crack tip neighborhood

1.850 1.000 0.541

0.231 0.125 0.068

10

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Fig. 12. Force-displacement curve for different meshes.

Fig. 13. G-Time curve for different meshes.

Fig. 14. Force-displacement curve varying time step.

the fact that the crack velocity is unknown in these computations. If the goal is to avoid time step dependence, methodologies such as the πθ method (crack growth stability) [54], that is an extension of the Gθ method, and fatigue analysis (crack growth under cyclic conditions) [55], that allow computation of the crack velocity, should be used. In addition to the time step dependence, the current approach also depends on the crack growth distance. Even though the influence of this parameter was not studied here, the oscillations observed in the force-displacement curves (Figs. 12 and 14) come from the fact that crack growth distance was set arbitrarily. The computation of the crack velocity is required if one wants to define the crack increment distance for a given time step. Compared to LEFM, the current approach has the advantage of allowing fracture energy control. Compared to CZMs, the main advantage is that crack propagation proceeds through an arbitrary direction. However it must be noted that the dynamic insertion of such cohesive elements introduces a numerical dependence of the crack increment and the associated time step. According to the cohesive law implemented in this work [25], the relevant parameters involved in the model are the maximum cohesive stress (σc ) and the critical opening displacement of the cohesive elements (δc ). These values are directly related through the 11

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Fig. 15. G-Time curve varying time step.

following equation describing the fracture energy [25]: (5)

Gc = eσc δc Eq. (5) can also be related to Eq. (6), defining the process zone length (PZL), as follows [56,57]:

PZL =

Gc E πσc2

(6)

Aiming at finding the effect on the crack propagation process when the aforementioned cohesive parameters are varied, several simulations using different values for them are performed. Table 3 shows the different values that were tested for the cohesive parameters during the simulations. Not counting the reference case (A), five different simulations were performed. For the first three cases, B1, B2 and B3 the critical cohesive stress (σc ) is varied. The fracture energy (G) was kept constant, so the new critical opening displacement of the cohesive elements (δc ) was computed in each case using Eq. (5). Variations of the critical cohesive stress (σc ) were defined trying to keep a logical value for the process zone length (PZL). In the remaining two cases, C1 and C2 , the fracture energy was doubled. In the first case, the critical cohesive stress was doubled, while in the second one, the doubled parameter was the critical opening displacement. Fig. 16 shows the different cohesive laws used in this work. Fig. 17 shows the different force-displacement curves obtained when cohesive parameters were varied. It is evident that the variation of these parameters did not lead to significant changes in the results. When fracture energy is kept constant, the influence of varying the maximum cohesive stress (σc ) and the critical opening displacement (δc ) is minimal. In the other scenario, when fracture energy was doubled through the increasing of either the maximum cohesive stress or the critical opening displacement, the initiation of the propagation started later, and the required force for the fracture process was greater. The same remark can be asserted, if the fracture energy remains constant, variation of cohesive parameters do not lead to relevant changes in the solution. When the critical opening displacement (δc ) decreases while fracture energy is kept constant, the slope (stiffness) of the traction separation law is significantly increased, as can be seen in Fig. 16. Numerically, this high stiffness might heavily impact the conditioning of the linear system and, therefore impact the convergence of the numerical solver. Inversely, when δc increases, the artificial stiffness is reduced. The crack paths obtained using the studied cohesive parameters are shown in Fig. 18. In terms of crack path, there are no significant differences. These results prove that the insertion of cohesive elements in the crack propagation process plays an important role in terms of the energy involved in the propagation and more so than in the crack path itself, alowing to simulate a more realistic fracture process.

Table 3 Values of the different physical parameters used in the cohesive law. Legend

G [MPa mm]

σc [MPa]

δc [mm]

PZL [mm]

A (Reference)

8.378·10−7

0.200

1.541·10−6

2.000

B1 B2

8.378·10−7

0.141

2.179·10−6

4.000

8.378·10−7

0.282

1.089·10−6

1.000

B3

8.378·10−7

0.400

7.705·10−7

0.500

C1

1.675·10−6

0.400

1.541·10−6

1.000

C2

1.675·10−6

0.200

3.081·10−6

4.000

12

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

Fig. 16. Different cohesive laws used here.

Fig. 17. Force-displacement curves for different cohesive parameters.

Fig. 18. Crack paths comparison for different cohesive parameters.

5. Conclusion A new methodology for 2D crack growth simulations under mixed-mode loading in brittle materials has been presented. It combines an advanced remeshing technique, that enables remeshing exactly in a computed direction, with dynamic insertion of cohesive elements in the remeshed zone. Additionally, the method presented here allows for an accurate modeling of the energy dissipation rate through a traction-separation law, as well as crack propagation using remeshing operations and cohesive elements in a mesh-independent way. Over the various applications presented in this work, the robustness and the accuracy of the proposed approach in terms of crack path is shown. By varying the mesh size, it is found that the presented approach exhibits mesh convergence. Compared to classical crack propagation techniques based on LEFM, the current approach offers the possibility to control fracture energy. This approach improves on classical CZM approaches since it enables to model configurations with arbitrary crack paths. However, as shown in this paper, the 13

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

dynamic insertion of cohesive elements introduces oscillations in the post-peak regime. These oscillations are merely due to the crack increment length and the time step. The aforementioned oscillations are the result of a dependence of the current approach on the time step, as well as on the constant crack increment used once G exceeds Gc . In order to overcome this problem, methodologies that allow for the computation of the crack velocity should be used. The πθ method (extension of the Gθ method with a double ring of integration [54]) or any kind of fatigue crack propagation law could be a suitable solution to avoid this dependency. Variation of the maximum cohesive stress (σc ) and of the critical opening displacement (δc ) while the fracture energy was kept constant, showed that these parameters play a purely numerical role. Thus, the fracture energy is the most important parameter to take into account when tuning the model. Future work will address the generalization of the presented methodology to multiple cracks and to 3D configurations, as well as real application examples in the engineering fracture mechanics field. Also, the dependence of the present method on the crack growth distance, as well as on the time step, along with the development of an appropriate technique will be investigated in a future study. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements We acknowledge support from Academies of Excellence on Complex Systems and Space, Environment, Risk and Resilience of the Initiative d’EXcellence ”Joint, Excellent, and Dynamic Initiative” (IDEX JEDI) of the Université Côte d’Azur and the Center for Planetary Origin (C4PO) A Material approach (https://www-n.oca.eu/morby/C4PO/C4PO-home.html). We also acknowledge support from The DIGIMU ANR industrial chair (https://chaire-digimu.cemef.mines-paristech.fr/). Appendix A. Stress field at the crack tip under mixed-mode loading In Fig. A.19 it is sketched the stress state for a differential element in the neighborhood of the crack tip in polar coordinates. For a general mixed-mode problem, the asymptotic near-tip stress field for the above case is given by:

(A.1) where KI and KII are respectively the stress intensity factors for mode I and mode II, θ is the angle of the vector that goes from the crack tip to the interest point and r is the distance between the crack tip and the interest point. According to the maximum circumferential stress criterion [44], which states that a crack in an elastic material is going to propagate in the direction for which the circumferential stress (σθθ ) is maximum, we proceed to maximize σθθ in (A.1) doing

∂σθθ ∂θ

= 0 and

∂2σθθ ∂θ2

θ 3θ θ 3θ KI ⎡sin ⎛ ⎞ + sin ⎛ ⎞ ⎤ + KII ⎡cos ⎛ ⎞ + cos ⎛ ⎞ ⎤ = 0 ⎢ ⎢ ⎝ 2 ⎠⎥ ⎝ 2 ⎠⎥ ⎦ ⎦ ⎣ ⎝2⎠ ⎣ ⎝2⎠

< 0:

(A.2)

Eq. (A.2) can be simplified to: (A.3)

KI sinθ + KII (3cosθ − 1) = 0 which can be solved leading to [49]:

Fig. A.19. Stress components near the crack tip in cylindrical coordinates. 14

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

θ = 2arctan

2 2 ⎛ KI − KI + 8KII ⎞ ⎜ ⎟ 4KII ⎝ ⎠

(A.4)

where θ is the actual crack propagation direction, the one that maximizes the circumferential stress. Appendix B. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.engfracmech. 2020.106918.

References [1] Erdogan F. Fracture mechanics. Int J Solids Struct 2000;37(1):171–83. [2] Zerbst U, Klinger C, Clegg R. Fracture mechanics as a tool in failure analysis – prospects and limitations. Eng Fail Anal 2015;55:376–410. [3] Liao M, Deng X, Guo Z. Crack propagation modelling using the weak form quadrature element method with minimal remeshing. Theoret Appl Fract Mech 2018;93:293–301. [4] Bouchard P-O, Bay F, Chastel Y, Tovena I. Crack propagation modelling using an advanced remeshing technique. Comput Methods Appl Mech Eng 2000;189(3):723–42. [5] Alessi R, Marigo J-J, Vidoli S. Gradient damage models coupled with plasticity: variational formulation and main properties. Mech Mater 2015;80:351–67. materials and Interfaces. [6] Marigo J-J, Maurini C, Pham K. An overview of the modelling of fracture by gradient damage models. Meccanica 2016;51(12):3107–28. [7] Tanné E, Li T, Bourdin B, Marigo J-J, Maurini C. Crack nucleation in variational phase-field models of brittle fracture. J Mech Phys Solids 2018;110:80–99. [8] Bouchard P-O, Bay F, Chastel Y. Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria. Comput Meth Appl Mech Eng 2003;192(35):3887–908. [9] Wulf J, Steinkopff T, Fischmeister H. Fe-simulation of crack paths in the real microstructure of an Al(6061)/SiC composite. Acta Mater 1996;44(5):1765–79. [10] Shakoor M, Navas VT, Muñoz DP, Bernacki M, Bouchard P-O. Computational methods for ductile fracture modeling at the microscale. Arch Comput Meth Eng 2018:1–40. [11] Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Meth Eng 1994;37(2):229–56. [12] Ingraffea AR, de Borst R. Stein E, de Borst R, Hughes TJR, editors. Computational fracture mechanics, encyclopedia of computational mechanics. 2nd ed.John Wiley and Sons; 2017. [13] Rashid M. The arbitrary local mesh replacement method: an alternative to remeshing for crack propagation analysis. Comput Meth Appl Mech Eng 1998;154(1):133–50. [14] Simo J, Oliver J, Armero F. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Comput Mech 1993;12(5):277–96. [15] Garikipati K. On strong discontinuities in inelastic solids and their numerical simulations Ph.D. thesis Stanford University; 1996. [16] Oliver J, Huespe A. Continuum approach to material failure in strong discontinuity settings. Comput Meth Appl Mech Eng 2004;193(30):3195–220. computational Failure Mechanics. [17] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45(5):601–20. [18] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46(1):131–50. [19] Kuhn C, Müller R. A continuum phase field model for fracture. Eng Fract Mech 2010;77(18):3625–34. [20] Borden MJ, Verhoosel CV, Scott MA, Hughes TJ, Landis CM. A phase-field description of dynamic brittle fracture. Comput Meth Appl Mech Eng 2012;217–220:77–95. [21] Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoret Appl Fract Mech 2018;96:174–92. [22] Nguyen T, Yvonnet J, Zhu Q-Z, Bornert M, Chateau C. A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure. Eng Fract Mech 2015;139:18–39. [23] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Dryden HL, von Kármán T, Kuerti G, van den Dungen FH, Howarth L, editors. Advances in Applied Mechanics, vol. 7. Elsevier; 1962. p. 55–129. [24] Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech 1968;35(2):379–86. [25] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Numer Methods Eng; 1999. [26] Tomar V, Zhai J, Zhou M. Bounds for element size in a variable stiffness cohesive finite element model. Int J Numer Meth Eng 2004;61(11):1894–920. [27] Chiaruttini V, Geoffroy D, Riolo V, Bonnet M. An adaptive algorithm for cohesive zone model and arbitrary crack propagation. Revue Européenne de Mécanique Numérique/Eur J Comput Mech 2012;21:208–18. [28] Choi H, Park K. Removing mesh bias in mixed-mode cohesive fracture simulation with stress recovery and domain integral. Int J Numer Meth Eng 2019:1–24. [29] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8(2):100–4. [30] Rice JR. Mathematical analysis in the mechanics of fracture. In: Liebowitz H, editor. Fracture: an advanced treatise. N.Y.: Academic Press; 1968. p. 191–311. chap. 3. [31] Needleman A. Micromechanical modelling of interfacial decohesion. Ultramicroscopy 1992;40(3):203–14. [32] Yi D, Xiao ZM, Idapalapati S, Kumar SB. Fracture analysis of girth welded pipelines with 3D embedded cracks subjected to biaxial loading conditions. Eng Fract Mech 2012;96:570–87. [33] Lorentz E. A mixed interface finite element for cohesive zone models. Comput Meth Appl Mech Eng 2008;198(2):302–17. [34] Alfano M, Furgiuele F, Leonardi A, Maletta C, Paulino GH. Cohesive zone modeling of Mode I fracture in adhesive bonded joints. In: Advances in Fracture and Damage Mechanics VI, vol. 348 of Key Engineering Materials, Trans Tech Publications; 2007. p. 13–6. [35] Turon A, Camanho P, Costa J, Renart J. Accurate simulation of delamination growth under mixed-mode loading using cohesive elements: definition of interlaminar strengths and elastic stiffness. Compos Struct 2010;92(8):1857–64. [36] Xu X-P, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42(9):1397–434. [37] Tijssens MG, Sluys BL, van der Giessen E. Numerical simulation of quasi-brittle fracture using damaging cohesive surfaces. Eur J Mech A. Solids 2000;19(5):761–79. [38] Geißler G, Netzker C, Kaliske M. Discrete crack path prediction by an adaptive cohesive crack model. Eng Fract Mech 2010;77(18):3541–57 computational Mechanics in Fracture and Damage: A Special Issue in Honor of Prof. Gross. [39] Shakoor M, Bernacki M, Bouchard P-O. A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: analysis of void clusters and stress state effects on coalescence. Eng Fract Mech 2015;147:398–417. [40] Hussain M, Pu S, Underwood J. Strain energy release rate for a crack under combined Mode I and Mode II, fracture analysis, proc. 1973 national symposium on fracture mechanics, Part II, ASTM STP 560. Amer Soc Testing Mater 1973;1973:2–28.

15

Engineering Fracture Mechanics xxx (xxxx) xxxx

D. Uribe-Suárez, et al.

[41] Destuynder P, Djaoua M, Lescure S. Quelques remarques sur la mécanique de la rupture élastique. J Méca Théo Appl 1983;2(1):113–35. [42] Li F, Shih C, Needleman A. A comparison of methods for calculating energy release rates. Eng Fract Mech 1985;21(2):405–21. [43] Raju I, Shivakumar K. An equivalent domain integral method in the two-dimensional analysis of mixed mode crack problems. Eng Fract Mech 1990;37(4):707–25. [44] Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J Basic Eng 1963;85(4):519–25. [45] Digonnet H, Silva L, Coupez T. Cimlib: a fully parallel application for numerical simulations based on components assembly. In: Cesar de Sa JMA, Santos AD, editors, Materials processing and design; modeling, simulation and applications; NUMIFORM ’07, vol. 908 of American Institute of physics conference series; 2007. p. 269–74. [46] Arnold DN, Brezzi F, Fortin M. A stable finite element for the stokes equations. CALCOLO 1984;21(4):337–44. [47] Brezzi F, Fortin M. Mixed and hybrid finite element methods, Springer series in computational mathematics. Springer; 1991. [48] Cao TS. Modeling ductile damage for complex loading paths, Ph.D. thesis, Ecole Nationale Supérieure des Mines de Paris; 2013. [49] Nguyen-Xuan H, Liu G, Nourbakhshnia N, Chen L. A novel singular ES-FEM for crack growth simulation. Eng Fract Mech 2012;84:41–66. [50] Barsoum RS. On the use of isoparametric finite elements in linear fracture mechanics. Int J Numer Meth Eng 1976;10(1):25–37. [51] Azócar D, Elgueta M, Rivara MC. Automatic LEFM crack propagation method based on local Lepp-Delaunay mesh refinement. Adv Eng Softw 2010;41(2):111–9. [52] Bittencourt T, Wawrzynek P, Ingraffea A, Sousa J. Quasi-automatic simulation of crack propagation for 2D LEFM problems. Eng Fract Mech 1996;55(2):321–34. [53] Ingraffea A, Grigoriu M. Probabilistic fracture mechanics: a validation of predictive capability, Tech. Rep., Department of Structural Engineering, Cornell University; 1990. [54] Suo X-Z, Combescure A. Double virtual crack extension method for crack growth stability assessment. Int J Fract 1992;57(2):127–50. [55] Lemaitre J, Desmorat R. Engineering damage mechanics: ductile, creep, fatigue and brittle failures. Berlin, Heidelberg: Springer Berlin Heidelberg; 2005. [56] Ha K, Baek H, Park K. Convergence of fracture process zone size in cohesive zone modeling. Appl Math Model 2015;39(19):5828–36. [57] Hermes F. Process zone and cohesive element size in numerical simulations of delamination in bi-layers Master’s thesis Eindhoven University of Technology; 2010.

16