Nicolas Guy

Darius M Seyedi

François Hild
email: hild@lmt.ens-cachan.fr

A probabilistic nonlocal model for crack initiation and propagation in heterogeneous brittle materials

Keywords: Brittle materials, Crack networks, Damage mechanics, Nonlocal model, Probabilistic approach, Weibull model

A probabilistic damage model is developed to study crack initiation and growth in quasi-brittle materials. Two different thresholds are considered to describe these mechanisms. A Weibull model is used to account for the randomness of crack initiation(s) and then a fracture mechanics based threshold is considered to model crack propagation. The model is integrated in a finite element code via a nonlocal damage approach. A regularization operator based on a stress regularization is introduced. Both damage thresholds are checked using the 'regularized' stress field to avoid mesh dependence. The interaction between propagating cracks and potential initiation sites is accounted for.

INTRODUCTION

Accurate modeling of damage and cracking of quasi-brittle materials such as concrete and rock is of great interest in different engineering applications. These materials exhibit important heterogeneities on their mechanical parameters because of their natural components. Their failure has usually a random character. The failure stress is scattered and the average level decreases with the volume of the sample. The effect of this scatter on the damage probability needs to be taken into account.

The description of strain localization is another crucial condition for simulating a damaged region with continuum models. In numerical simulations with standard softening models the strain field localizes in a band with a thickness of one element [START_REF] Pietruszczak | Finite Element Analysis of Deformation of Strain-softening Materials[END_REF]. Every mesh refinement modifies the global response of the structure. The results exhibit a pathological dependence on the mesh fineness, the orientation of elements and in a more general form on the spatial discretization. As an extreme case, it can be shown that for a very fine discretization the predicted fracture energy approaches zero [START_REF] Pietruszczak | Finite Element Analysis of Deformation of Strain-softening Materials[END_REF][START_REF] Bažant | Wave propagation in strain-softening bar: exact solution[END_REF], which is unrealistic from a physical point of view. The corresponding mathematical problem has been identified as a loss of ellipticity of the mechanical problem [START_REF] Triantafyllidis | A gradient approach to localization of deformation: I. Hyperelastic materials[END_REF][START_REF] Lasry | Localization limiters in transient problems[END_REF][START_REF] Borst | Fundamental issues in finite element analysis of localization of deformation[END_REF]. As early as 1976, Hillerborg et al. proposed a Fictitious Crack Model [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF], which provides a mesh-independent energy release upon crack propagation. This concept was adapted by Bažant and Oh [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF] to a smeared crack formulation in which the fracture energy is smeared out over the area where the damaged zone localizes.

Different regularization schemes have been used to avoid the ill-posedness of the mathematical description at a certain level of accumulated damage. Nonlocal damage models [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bažant | Nonlocal continuum damage, localization instability and convergence[END_REF] provide a suitable framework to have mesh objective results when dealing with the post-localized behavior. Nonlocal principals have also been used to provide enhanced plasticity models [START_REF] Borst | Gradient dependent plasticity: formulation and algorithmic aspects[END_REF][START_REF] Manzari | Gradient plasticity modeling of geomaterials in a meshfree environment. Part I: Theory and variational formulation[END_REF]. The higher-order deformation gradients are used to include a nonlocal effect regularizing strain localization, and thus removing the sensitivity of the model to spatial discretizations [START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF][START_REF] Lorentz | Gradient constitutive relations: numerical aspects and application to gradient damage[END_REF].

The main idea of nonlocal models is to consider that the behavior of each material point depends on the state of surrounding points. Let us assume that damage is governed by strains. In conventional local models, damage is usually related to a scalar measure of the local strain state. As an example, in a nonlocal model, damage is related to a weighted volume average of this equivalent strain through a regularization operator. The choice of such a regularization operator does not rely on physical observations. The implicit gradient enhanced scheme proposed by Peerlings et al. [START_REF] Peerlings | Computational modelling of gradient-enhanced damage for fracture and fatigue problems[END_REF][START_REF] Peerlings | Enhanced damage modelling for fracture and fatigue[END_REF] provides a regularization of the strain by resorting to a minimum square method and reads for a domain Ω

R (ϵ ϵ ϵ) = min ϵ ϵ ϵ ∫ Ω ( 1 2 (ϵ ϵ ϵ -ϵ ϵ ϵ) 2 + 1 2 (ℓ c ∇ ∇ ∇ϵ ϵ ϵ) 2 ) dΩ (1)
where ϵ ϵ ϵ and ϵ ϵ ϵ are strain and regularized strain tensors respectively, ℓ c a characteristic length, and ∇ ∇ ∇ the gradient operator. The gradient term introduces the interaction between each material point and its neighborhood. Minimizing the integral of Equation ( 1) is equivalent to solving the following differential equation

ϵ ϵ ϵ -ℓ 2 c ∆ϵ ϵ ϵ = ϵ ϵ ϵ (2)
where ∆ denotes the Laplacian operator, with the following boundary conditions

(∇ ∇ ∇ ϵ ϵ ϵ) • n = 0 (3)
with n a vector normal to the boundary. It can be noted that such regularization scheme will increase the strain level in a point in the vicinity of the damaged zone and limits the strain gradient. This procedure will enlarge the band of strain localization. Some modified versions include a projector canceling out the contribution of the completely damaged elements to the regularization procedure to prevent a physically non-realistic enlargement of the localization band [START_REF] Godard | Modélisation de l'endommagement anisotrope du béton avec prise en compte de l'effet unilatéral : Application à la simulation des enceintes de confinement nucléaires[END_REF].

The effect of a statistical description of the continuum on the behavior of a nonlocal model at and after the onset of localization has already been studied [START_REF] Carmeliet | Stochastic approaches for damage evolution in standard and non-standard continua[END_REF][START_REF] Carmeliet | Probabilistic nonlocal damage model for continua with random field properties[END_REF]. The considered studies were performed using a nonlocal model based on a regularized strain. The material macro-scale heterogeneity is described based on the framework of the random field theory [START_REF] Vanmarcke | Random Fields: Analysis and Synthesis[END_REF]. It has been shown that a stochastic description of the damage threshold does not solve the difficulties associated with strain localization and the use of a nonlocal scheme is necessary. Two internal lengths associated respectively with the internal length scale of the nonlocal continuum and the correlation length of the random field must be considered. The authors pointed out the proper modeling of the relationship between these internal lengths and their experimental identification as an important issue [START_REF] Carmeliet | Stochastic approaches for damage evolution in standard and non-standard continua[END_REF].

The statistical theory of strength proposed by Weibull [START_REF] Weibull | A Statistical theory of the strength of materials[END_REF] is largely used to describe the random nature of the failure of quasi-brittle materials. However it is shown that the classical Weibull theory cannot be used for situations where the failure of the structure is caused by a stable propagation of a macrocrack [START_REF] Bažant | Statistical size effect in quasi-brittle structures: I. Is Weibull theory applicable[END_REF]. More particularly, the stress field used to calculate the failure probability in each point must take into account the stress redistribution due to the presence of macrocracks and reflect the localization of strains occurring prior to failure. A nonlocal Weibull theory was proposed by Bažant and Xi for improving this shortcoming [START_REF] Bažant | Statistical size effect in quasi-brittle structures: II. Non-local theory[END_REF]. In this setting, an average stress value is used in the fracture-process zone around the crack tip. The local (i.e., Cauchy's) stress is used elsewhere. An extension of the proposed analytical solution to a finite element nonlocal damage model is also proposed. Considering the elastic analysis of fracture, the nonlocal strains are used to calculate the average stresses in the process zone as follows [START_REF] Bažant | Statistical size effect in quasi-brittle structures: II. Non-local theory[END_REF][START_REF] Bažant | Probabilistic nonlocal theory for quasibrittle fracture initiation and size effects. I: Theory[END_REF]

P f = exp [ - ∫ Ω n ∑ i=1 ( ⟨Eε i (x) -σ u ⟩ σ 0 ) m dΩ(x) Ω r ] (4)
where P f denotes the failure probability, n number of dimension (n = 1, 2 or 3), E Young's modulus, ε(x) nonlocal strain, σ u , σ 0 and m Weibull parameters, and Ω r a reference volume.

As it can be seen, the product of Eε(x) is considered as a nonlocal measure of the stress field in this setting. Two other alternatives of nonlocal averaging are also proposed by Bažant and Novák [START_REF] Bažant | Probabilistic nonlocal theory for quasibrittle fracture initiation and size effects. I: Theory[END_REF]. The first one consists in considering only some measures of damage or inelastic strains [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bažant | Nonlocal continuum damage, localization instability and convergence[END_REF] for the nonlocal constitutive model. The inelastic stress defined as 'the elastic stress corresponding to the given strain minus the actual stress for the same strain according to the given stressstrain diagram' [START_REF] Bažant | Probabilistic nonlocal theory for quasibrittle fracture initiation and size effects. I: Theory[END_REF] is subjected to spatial averaging. The second one consists in spatially averaging the inelastic strain, representing the difference of the current strain and the elastic strain corresponding to the same stress [START_REF] Bažant | Probabilistic nonlocal theory for quasibrittle fracture initiation and size effects. I: Theory[END_REF]. A uniform weighting function that vanishes out of a certain representative volume is used for averaging the nonlocal variables. More recently, this formalism is used by Bažant [START_REF] Bažant | Probability distribution of energetic-statistical size effect in quasibrittle fracture[END_REF] to derive the probability distribution of the size effect on the nominal strength of heterogeneous quasi-brittle materials.

It is worth noting that the aforementioned works are based on the failure probability defined as the crack initiation probability. Even though an example is presented [START_REF] Bažant | Statistical size effect in quasi-brittle structures: II. Non-local theory[END_REF] in which the asymptotic stress field around a crack tip is used to fit the averaging functions, no link is made to the stress intensity factor (SIF) calculation or crack propagation. Fracture mechanics is a powerful tools for describing the propagation of initiated cracks. Different numerical schemes are used to model crack propagation, which is based upon fracture mechanics principles. Among them, the double noding technique [START_REF] Kanninen | The numerical simulation of crack growth in weldinduced residual stress fields[END_REF][START_REF] Liaw | Double nodding technique for mixed mode crack propagation studies[END_REF], local or general remeshing [START_REF] Rashid | The arbitrary local mesh replacement method: An alternative to remeshing for crack propagation analysis[END_REF][START_REF] Seyedi | Numerical modeling of crack propagation and shielding effects in a striping network[END_REF], boundary elements [START_REF] Aliabadi | A new generation of boundary element methods in fracture mechanics[END_REF], cohesive zone models [START_REF] Charlotte | Initiation of cracks with cohesive force models: a variational approach[END_REF] and the extended finite element method [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF] can be mentioned. The models based on the traditional finite element method require mesh refinement in order to model the singularity of the stress field around the crack tip and to calculate correctly the value of stress intensity factors. However, the numerical description of propagating cracks remains a scientific and technical challenge.

In the present paper a new nonlocal damage model is proposed. The model aims at simulating crack initiation and propagation within a heterogeneous material. Crack initiation follows a Poisson-Weibull process. A local brittle damage model is introduced. It provides results following the weakest link hypothesis. The propagation of the initiated cracks is then modeled through fracture mechanics principles. The stress intensity factor is considered as the governing parameter for crack growth. It will be shown that the use of a simple local damage model with a nonlocal gradient enhanced regularization enables for the simulation of brittle failure of weakest link type, and provides mesh-objective results for crack initiation. In this view, a new regularization scheme based on a stress regularization operator is proposed. The stress gradient regularization scheme is very close to the deformation regularization proposed by Peerlings et al. [START_REF] Peerlings | Computational modelling of gradient-enhanced damage for fracture and fatigue problems[END_REF][START_REF] Peerlings | Enhanced damage modelling for fracture and fatigue[END_REF]. However, the proposed stress gradient regularization is a more realistic modeling of damage localization since no regularization is performed on the strain field. The regularization yields results that are independent of the spatial discretization with respect to the Weibull distribution properties. The developed model calculates the initiation probability of each element following a Poisson-Weibull model while accounting for the stress redistribution due to the presence of the propagating cracks. Furthermore, it will be shown that with the stress gradient enhancement of a simple local damage model the stress intensity factor and crack opening displacement are evaluated correctly without mesh refinement around tips of crack networks. A unique variable, i.e., the regularized stress will be used for evaluating crack initiation and propagation.

DESCRIPTION OF THE NONLOCAL MODEL

The main idea of the present approach is to use the same parameter for evaluating two different thresholds for crack initiation and propagation. The inception of propagation on an existing defect is considered as crack initiation. Once a crack has initiated, its growth is governed by a regularized stress intensity factor. A gradient enhanced stress is used to check initiation and propagation. The considered regularized stress is obtained by solving the following equation

σ σ σ -ℓ 2 c ∆σ σ σ = σ σ σ (5)
where σ σ σ is the regularized stress tensor and σ σ σ Cauchy's stress tensor. Furthermore, natural boundary conditions are considered

(∇ ∇ ∇ σ σ σ) • n = 0. (6)
Equation ( 5) is a Helmholtz's equation. Considering the introduced boundary conditions (6), the regularized field has been shown to correspond to that obtained by using a nonlocal model with Green's function as a weight function [START_REF] Peerlings | A critical comparison of non-local and gradient-enhanced softening continua[END_REF]. In the present model, the regularized stress is used for both the crack initiation and growth. The main objective of this section is to introduce both thresholds and to compare them.

CRACK INITIATION THRESHOLD

Crack initiation is described by an inverse Weibull law [START_REF] Weibull | A Statistical theory of the strength of materials[END_REF]. ) m [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF] where m stands for Weibull modulus, and σ m 0 /λ 0 is a scale parameter. A small value of m corresponds to a very heterogeneous material. Conversely, high Weibull moduli lead to more deterministic values. The probability of finding N d critical defects (i.e., the defects that may initiate a crack) in a considered domain Ω of size |Ω| is described by a Poisson distribution

P N d (Ω, σ I ) = [|Ω|λ t (σ I )] N d N d ! exp [-|Ω|λ t (σ I )] (8)
By making the weakest link hypothesis for crack initiation, the initiation probability of a domain is written as a function of the probability of finding at least one critical defect in this domain [START_REF] Freudenthal | Statistical approach to brittle fracture[END_REF]

P i (|Ω|, σ I ) = P N d ≥1 (Ω, σ I ) = 1 -P N d =0 (Ω, σ I ) = 1 -exp [-|Ω|λ t (σ I )] (9)
By substituting the definition of the total defect density [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF] in the expression of the initiation probability ( 9), a threshold value of the regularized stress for crack initiation is obtained. In the present model, a crack will initiate in a considered domain (i.e., in a finite element e) of size |Ω e | if the regularized maximum principal stress reaches the initiation threshold

S i = σ 0 (λ 0 |Ω e |) 1 m [-ln (1 -P i )]
1 m [START_REF] Borst | Gradient dependent plasticity: formulation and algorithmic aspects[END_REF] where

P i ∈ [0; 1]
. This threshold can be compared with the regularized maximum principal stress only if the initial defect size is less than characteristic length ℓ c . Furthermore, crack initiation is considered to occur at the center of the element, and the crack propagation direction is normal to that of the maximum eigen regularized stress. It is worth noting that the introduced crack initiation threshold leads to the same behavior for a domain independently to its spatial discretization. Let us consider a domain Ω uniformly stressed, and a set of n independent sub-domains Ω i such that ∪ n i=1 Ω i = Ω. For a given loading level, the probability for the first crack to initiate in the union of the n independent sub-domains Ω i and in the domain Ω are equal. This equality is obtained by

P i (| ∪ n i=1 Ω i |, σ I ) = 1 - n ∏ i=1 P N d =0 (|Ω i |, σ I ) = 1 -exp [ - n ∑ i=1 |Ω i |λ t (σ I ) ] = P i (|Ω|, σ I ) . (11)

CRACK GROWTH THRESHOLD

Let us consider the 2-dimensional case of a crack in an elastic medium submitted to a far field loading. In that case, Westergaard's solution [START_REF] Kanninen | Advanced Fracture Mechanics[END_REF] (i.e., the asymptotic solution in the vicinity of the crack tip) gives a good approximation of the stress field. The stress field given by Westergaard's solution is expressed in a basis located at the crack tip (figure 1)

σ σ σ (r, θ) = K I √ 2πr f I (θ) + K II √ 2πr f II (θ) (12)
where K I and K II denote the stress intensity factors corresponding to modes I and II, and f I and f II associated tensors. Tensors f I and f II read

f I (θ) =   cos θ 2 ( 1 -sin θ 2 sin 3θ 2 ) sin θ 2 cos θ 2 cos 3θ 2 0 sin θ 2 cos θ 2 cos 3θ 2 cos θ 2 ( 1 -sin θ 2 sin 3θ 2 ) 0 0 0 f I33   (13)
for mode I with f I33 = 0 under plane stress hypothesis and f I33 = 2 ν cos θ/2 under plane strain hypothesis, and ν is Poisson's ratio,

f II (θ) =   -sin θ 2 ( 2 + cos θ 2 cos 3θ 2 ) cos θ 2 ( 1 -sin θ 2 sin 3θ 2 ) 0 cos θ 2 ( 1 -sin θ 2 sin 3θ 2 ) sin θ 2 cos θ 2 cos 3θ 2 0 0 0 f II33   (14)
for mode II with f II33 = 0 under plane stress hypothesis and f II33 = -2 ν sin θ/2 under plane strain hypothesis. Let us consider Helmholtz' Equation ( 5) and turn it into its integral form that has an analytical solution for infinite plane problems [START_REF] Polyanin | Handbook of linear partial differential equations for engineers and scientists[END_REF]. For an infinite 2D medium, the regularized stress for a given point P is written in a cylindrical coordinate system as

σ σ σ(P ) = 1 2π ∫ ∞ 0 ∫ π -π σ σ σ (r, θ) ℓ 2 c K bζ ( r ℓ c ) rdrdθ ( 15
)
where K bζ is modified second order Bessel function

K b0 (a r ) = ∫ ∞ 0 cos a r t √ 1 + t 2 dt ( 16
)
for ζ = 0 [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs and mathematical tables[END_REF]. Modified second order Bessel function plays the role of a weight function in this formulation. For a point P at the tip of a crack with a sufficiently small characteristic length, Equation ( 12) becomes

σ σ σ(P ) (K I , K II ) = II ∑ i=I [ K i (2π) 3 2 ℓ 2 c ∫ ∞ 0 K b0 ( r ℓ c ) r 1 2 dr ∫ π -π f f f i (θ) dθ ] . ( 17
)
Copyright c ⃝ 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1-6 Prepared using nmeauth.cls so that

σ σ σ(P ) (K I , K II ) = Γ 2 ( 3 4 ) 5π 1 √ πℓ c   4K I 4K II 0 4K II 6K I 0 0 0 0   + σ 33 e 3 ⊗ e 3 (18)
where σ 33 = -ν (σ 11 + σ 22 ) under plane strain condition, σ 33 = 0 under plane stress condition, and Γ denotes Euler (Gamma) function of the second kind. Equation [START_REF] Carmeliet | Probabilistic nonlocal damage model for continua with random field properties[END_REF] shows that the proposed operator provides a regularized stress at the crack tip corresponding to Westergaard's asymptotic solution. This solution is relevant when the characteristic length is small compared to the crack length. It is worth noting that the proposed formulation allows us to calculate directly the stress intensity factor without any mesh refinement. Therefore a natural transition is made from a damage mechanics based model for crack initiation to a fracture mechanics based description for crack propagation. Furthermore, both models use the same variable, namely, the regularized stress.

In order to describe crack growth in the opening mode, the following crack growth threshold is compared to the maximum principal stress

S g = 6Γ 2 ( 3 4 ) 5π K c √ πℓ c ( 19
)
where K c is the fracture toughness of the studied material. Furthermore, crack growth is assumed to occur at an existing crack tip and its direction is perpendicular to the regularized maximum eigen stress direction. For a pure mode I loading case, the given threshold is equivalent to

K c ≤ K I ( 20
)
when the characteristic length ℓ c is small compared to the crack length. It is to be noted that the introduced threshold is a mixed mode criterion because for a mixed mode loading the regularized maximal principal stress is

σ I = 6Γ 2 ( 3 4 ) 5π 1 √ πℓ c [ 5K I + √ K 2 I + 16K 2 II ] 6 . ( 21
)
The introduced threshold is therefore equivalent to the following criterion

K c ≤ [ 5K I + √ K 2 I + 16K 2 II ] 6 ≈ K I + 2 3 K II . ( 22
)

LOCAL CONSTITUTIVE LAW

A simple damage law describing a perfectly brittle behavior of an isotropic material is used as a local damage law. It is derived by using Helmholtz' state potential [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]

ρψ e = 1 2 (1 -d) ϵ ϵ ϵ : C : ϵ ϵ ϵ ( 23
)
where C is the elasticity tensor of the virgin material, ρ the mass density, and d the damage variable. From the state potential, the elasticity law is obtained

σ σ σ = ρ ∂ψ e ∂ϵ ϵ ϵ = (1 -d) C : ϵ ϵ ϵ. ( 24
)
The thermodynamic force associated with the damage variable d is defined as

Y = -ρ ∂ψ e ∂d = 1 2 ϵ ϵ ϵ : C : ϵ ϵ ϵ. ( 25
)
Damage growth is driven by

d = H (⟨σ I -S i ⟩ + ⟨σ I -S g ⟩) when ḋ ≥ 0 ( 26
)
where H denotes Heaviside step function, σ I the maximal principal stress, S i and S g the initiation and crack growth thresholds.

THRESHOLD COMPARISON

The initiation threshold was introduced by considering that a crack initiation stress is associated with each initial defects. In order to compare the initiation threshold with the crack growth threshold, let us consider that the most critical defect in a considered domain has a threshold stress corresponding to a crack of half length a, which is orientated perpendicular to the eigen stress direction. In that case, the initiation threshold is expressed as

S i = K c √ πa ( 27
)
Assuming that the characteristic length is greater than a leads to

S i = K c √ πa ≥ K c √ πℓ c > 6Γ ( 3 4
) 5π

K c √ πℓ c = S g . (28)
A characteristic length greater than the size of the largest initial crack (or defect) provides a crack initiation threshold greater than or equal to the crack growth threshold. It is worth noting that this condition is not initially satisfied because of the conservative nature of a two parameter Weibull law. Using a two-parameter Weibull law (i.e., without threshold) is equivalent to considering the possibility of having initial cracks (i.e., defects) of any size, even larger than the characteristic length. In the proposed setting, the characteristic length stands for a micro-structural parameter that is the size of the largest initial crack modeling a defect.

NUMERICAL IMPLEMENTATION

The finite element method is used to solve the global problem using the discretization described in sub-section 3.1. A Newton-Raphson scheme is chosen. The constitutive law is solved implicitly [START_REF] Abbas | Algorithme non-linéaire quasi-statique[END_REF] and the damage field is considered uniform in each element. The numerical model is introduced in Code Aster [39].

SPATIAL DISCRETIZATION OF THE FIELD EQUATIONS

A specific finite element, initially developed for a strain gradient enhancement model [START_REF] Galenne | Modélisation non locale à gradient de déformation[END_REF], is adapted for the spatial discretization. As shown in Figure 2 the regularized stress field is added as nodal degrees of freedom (DOFs). The introduced finite element is thus quadratic for the displacement field and linear for the regularized stress. P 2 and P 1 shape functions u u, σ σ are used for displacement and regularized stress fields, respectively. Two types of unknowns are defined on the corner nodes and only displacements are calculated at middle nodes. A triangular element has thus 21 DOFs for a plane problem. Two sets of shape functions and derivatives are associated with displacement (N and B) and regularized stress ( N and B). A detailed description of the shape functions and integration points can be found in Ref. [START_REF] Desroches | Fonctions de forme et points d'intégration des éléments finis[END_REF]. Two equations govern the global problem, namely, classical equilibrium and those corresponding to the calculation of regularized stress. The integral formulation of the problem reads ∀u * ∈ V ad ,

∫ Ω σ σ σ : ϵ ϵ ϵ (u * ) dΩ = ∫ G u * • T dG + ∫ Ω u * • f V dΩ (29)
where V ad is the space of admissible displacements, T the tractions, f V the body forces and

∀s s s * ∈ [ H 1 (Ω) ] 6 , ∫ Ω ( s s s * σ σ σ + ∇ ∇ ∇s s s * • ℓ 2 c σ σ σ ) dΩ = ∫ Ω s s s * σ σ σ dΩ. (30)
In the following, σ σ σ, ϵ ϵ ϵ et σ σ σ stand for vector forms of stress, strain and regularized stress tensors, u denotes the displacement vector. The global resolution consists in minimizing a residual force associated with displacements that can be calculated for each element as

F u = F int + D T λ λ λ -F ext ( 31
)
where F int = ∫ Ω B T σ σ σ dΩ is the resultant of the internal forces, F ext = ∫ Γ N T T dΓ is the resultant of the external loads with T the nodal forces. Operator D is defined so that Du = u d , with u d the prescribed displacements, and λ λ λ Lagrange multiplier corresponding to Dirichlet boundary conditions. A regularized stress residual is also minimized

F σ = K σσ σ σ σ -F σ ( 32
)
where

F σ = ∫ Ω NT σ σ σ dΩ and K σσ = ∫ Ω ( NT N + ℓ 2 c BT B) dΩ.
The following tangent matrix is used for minimizing the above mentioned residuals

K = [ ∂F u ∂u ∂F u ∂σ σ σ ∂F σ ∂u ∂F σ ∂σ σ σ ] (33)
with

∂F u ∂u = ∫ Ω B T ∂σ σ σ ∂ϵ ϵ ϵ B dΩ ( 34
)
∂F u ∂σ σ σ = ∫ Ω B T ∂σ σ σ ∂σ σ σ N dΩ = 0 ( 35
)
∂F σ ∂σ σ σ = ∫ Ω ( NT N + ℓ 2 c BT B) dΩ (37)
After global inversion, displacements and regularized stresses are known on the associated nodes. In the proposed setting, the growth of the state of each element is evaluated by taking into account the regularized stresses, the geometry of the element and its neighborhood. A constant damage field per element is considered. When a damage variable is calculated for each element (and constant for all Gauss points of the element), the terms ∂σ σ σ/∂ϵ ϵ ϵ are calculated for each Gauss point. Finally, the global matrix is assembled. Since the damage variable has a binary nature, the material is described either by an elastic behavior with initial rigidity or by a completely damaged material with a residual rigidity, which is negligible compared to the initial rigidity (d = 1 -10 -6 ). The residual rigidity is introduced since it is impossible to invert the global matrix using a null rigidity with conventional methods. It can be noted that if damage does not change during a loading step, the global problem in displacement or regularized stress is linear. The tangent matrix is thus updated at each Newton iteration and not for the elastic prediction at the beginning of a time step.

NUMERICAL CRACK DESCRIPTION

Each undamaged element can be damaged in two different situations during the cracking process. An element may be located in the direct neighborhood of a crack tip. In that case, the element is marked and the crack propagation threshold will be considered for that element. Conversely, i.e., there is no propagating crack in the direct neighborhood of the element, the initiation threshold is considered. Consequently, the geometrical position of crack tips and their corresponding elements are stored during the calculation. Figure 3 shows the developed procedure for defining the neighboring elements during initiation (Figure 3(a)) and propagation phases (Figure 3(b)). Let us consider that the applied load is sufficient for initiating a macrocrack in element A el . A crack will initiate at the center of the element, and the damage variable for the considered element will be d = 1. Let us define n IApt as a vector corresponding to the direction of the maximum principal 'regularized' stress σ IApt at point A pt when the failure threshold is reached. It is assumed that the initiated crack propagates perpendicular to n IApt . Two crack tips are defined, namely, B pt and C pt that are the intersections between a line ∆ A el with the element edges. Line ∆ A el is defined as the segment that passes by A pt and is perpendicular to vector n IApt . Each crack tip is located on a neighboring element, i.e., B el and C el . These neighboring elements are thus subjected to the propagation threshold and not the initiation one because their failure is associated with crack propagation.

As shown in Figure 3(b), if the propagation threshold is reached on one of the crack tips D pt located on element D el , the damage variable in this element will be d = 1. The crack tip will move along a line ∆ D el . Line ∆ D el passes by D pt and is perpendicular to n IDpt , which is the direction associated with the maximum principal regularized stress σ IDpt at point D pt when the failure threshold is reached. The crack tip moves to E pt in element E el . The developed method for the treatment of neighboring elements enables us to handle several cracks at the same time. Each crack direction may also change during the propagation due to structural effects or to the interaction with other cracks. The calculations are performed based on the nodal values of the regularized stress for each element. The strain is estimated at Gauss points from the nodal values of displacements each time that the damage variable is updated. The tangent behavior is then evaluated at the Gauss points.

In this method, a finite element is marked for propagation if and only if it is marked during previous Newton iterations. Each crack tip may thus induce the failure of one element during each iteration. The material behavior is implemented in an implicit manner to ensure the stability of converged states. During each time step, a crack tip may cause the failure of as much elements that is necessary for reaching equilibrium. The unstable propagation of the macrocracks thus can be simulated by the present model. A time step subdividing criterion based on the damage growth on the whole structure is also introduced. Applying this criterion does not allow for the initiation of more than one crack during each Newton iteration. This function enables us, for stable crack propagation conditions, to obtain results independent of the time step selected by the user.

COMPARISON OF NONLOCAL AND USUAL THRESHOLDS

In Section 2.2 a crack growth threshold that is equivalent to a critical stress intensity factor when the characteristic length ℓ c becomes small with respect to the crack length is proposed. In this part, a simulation aims at comparing the nonlocal threshold associated with crack growth with the Weibullian initiation threshold and a stress intensity factor based threshold. The objectives are to estimate the accuracy of the nonlocal threshold associated with crack growth for characteristic lengths that are small compared to the crack length, and to evaluate the influence of this threshold for longer characteristic lengths. The analysis is based on the case of a beam of width 2b and of infinite (i.e., very large) length that contains a central crack of length 2a and is submitted to a horizontal stress of σ ∞ . The treated case is described in Figure 4. For the considered case, an approximate analytical solution exists. The stress intensity factor is estimated as [START_REF] Murakami | Stress Intensity Factor Handbook[END_REF]

K Ia ≈ σ ∞ √ πa cos ( πa 2b ) . ( 38
)
For the numerical simulations, half of the domain is meshed. The vertical displacements are set to zero on the lower edge (i.e., on the horizontal plane of symmetry). The horizontal displacements are set to zero on the left edge and the right edge is submitted to a horizontal stress. The half crack is modeled by a line of broken elements (i.e., elements where d = 1). For a given characteristic length, a stress intensity factor is deduced from the nonlocal threshold

K I (ℓ c ) ≈ 5π 6Γ 2 ( 3 4
) σ I √ πℓ c ( 39
)
with σ I the maximum principal regularized stress at the crack tip (see Section 3.2). In order to compare the 'nonlocal based' stress intensity factor to the usual one, the ratio

R I (ℓ c ) = K Ia /K I (ℓ c ) is introduced.
The results of simulations performed for three characteristic lengths ℓ c of 20 mm, 40 mm and 80 mm are shown in Figure 5. As shown in Figures 5(a analytical solution. Furthermore, for each characteristic length the nonlocal stress intensity factor tends to be equal to the usual stress intensity factor when the crack length increases. It also appears that the value of the maximum principal regularized stress is larger for lower characteristic lengths. For each case, an approximation of R I (ℓ c ) is plotted in Figure 5(b). This approximation reads

R I (ℓ c ) = 1 -exp [ - ( a ℓ c ) 2 3 ] , ( 40
)
and gives a good estimation of the characteristic length effect when a ≥ ℓ c . Further, the nonlocal threshold associated with crack growth is such that lim a/ℓc→0

K I (ℓ c ) = 5π 6Γ 2 ( 3 4 ) σ ∞ √ πℓ c . ( 41
)
Considering the approximations ( 40) and ( 41) an estimation of the nonlocal threshold associated with crack growth is plotted in Figure 6. In order to compare different thresholds, the usual stress intensity factor threshold and some examples of Weibull-based threshold describing crack initiation are plotted in Figure 6. Micro and macro scales are defined in Figure 6. The numerical strategy proposed in Section 3 is such that the microcracks (i.e., cracks with a ≤ ℓ c ) are not described by the damage field but associated with a Weibull threshold. The threshold of microcrack initiation is always greater than the minimum initiation threshold plotted in Figure 6 because the microcracks are assumed to be small (i.e., a ≤ ℓ c ). Furthermore, in the developed numerical strategy, the macrocracks (i.e., cracks with a ≥ ℓ c ) are represented by the damage field. When a crack initiates, the crack description changes and crack growth is then simulated by considering the nonlocal propagation threshold that tends to be equal to the usual stress intensity factor threshold as the crack grows. The numerical treatment of the threshold change is described in Section 3.2.

The crack opening displacement is another parameter to be checked. If the crack of Figure 4 is located in an infinite plate, its opening displacement ∆u x is calculated by resorting to the analytical solution proposed by Muskhelishvili [START_REF] Kanninen | Advanced Fracture Mechanics[END_REF]

∆u x = 4σ ∞ a E √ 1 - ( a y ) 2 (42)
where y denotes the coordinate of each point on the crack mouth from the crack center. The crack opening at its center can be also calculated as a function of the applied stress and crack geometry (i.e., a/b ratio) [START_REF] Murakami | Stress Intensity Factor Handbook[END_REF]. The crack opening calculated by the model is compared to the two solutions for a large plate (a/b = 0.14) on figure 7. A very good agreement is observed between the model prediction and the analytical results. Such accurate results are important in view of coupling the developed model with an external flow simulator for coupled hydromechanical applications.

SIMULATION RESULTS

In this section, two parts presenting simulations using the introduced model are described. The first one aims at validating the independence of the developed model to spatial discretizations. This property is studied with a particular attention since it has been the main reason for using a nonlocal formulation. The second part aims at demonstrating the ability of the numerical model to describe the formation and propagation of a crack network by considering simultaneously both initiation and growth of several cracks.

Ref. [START_REF] Murakami | Stress Intensity Factor Handbook[END_REF] Numerical simulation Ref. [START_REF] Kanninen | Advanced Fracture Mechanics[END_REF] Figure 7. Comparison between the simulated and analytical crack opening displacement (a/b = 0.14, ℓc = 20 mm).

VALIDATION OF SPATIAL DISCRETIZATION INDEPENDENCE

In Section 2, it has been assumed that thresholds based on the regularized stress describe crack initiation and growth without a spatial discretization dependence. The aim of this section is to validate this characteristic of the model.

CRACK INITIATION

The crack initiation threshold is studied considering the failure stress obtained for a domain of volume V submitted to a uniform tensile stress for various discretizations. For each considered volume, the initiation threshold is defined with respect to an inverse Weibull law (Equation ( 10)). For the different studied cases, the considered domain is divided into 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, and 1024 elements. For each spatial discretization, several calculations are performed by considering different numbers of random selections for the inverse Weibull law (i.e., 10, 100, 1,000 and 10,000 realizations). The numerical results in terms of normalized mean failure stress and standard deviation of the initiation stress for a set of simulations are plotted in Figure 8. The mean initiation stress is normalized with respect to the value σ i associated with a Weibull model

σ i = σ 0 (V λ 0 ) 1/m Γ ( 1 + 1 m + 1 ) (43)
and the standard deviation is normalized by

σ i = σ 0 (V λ 0 ) 1/m √ Γ ( 1 + 2 m + 1 ) -Γ 2 ( 1 + 1 m + 1 ) . ( 44
)
The convergence is faster for the mean initiation stress than the corresponding standard deviation. These results are compared with those given by the Weibull model. The results do not depend on spatial discretization thanks to the volume effect accounted for in Equation [START_REF] Borst | Gradient dependent plasticity: formulation and algorithmic aspects[END_REF].

It is also a direct consequence of the weakest link concept described by Equation [START_REF] Manzari | Gradient plasticity modeling of geomaterials in a meshfree environment. Part I: Theory and variational formulation[END_REF], an underlying hypothesis used in a Weibull model. Furthermore, it can be seen that the numerical results converge to the analytical results when the number of considered random selections increases. These results show that the main influence is given by the number of realizations, and to get a reliable estimate of statistical parameters (e.g., number of initiated cracks, size of propagating cracks), a convergence study is desirable.

CRACK GROWTH

In this section, the propagation of a single crack is simulated using different spatial discretizations. The chosen meshes allow to show that the results depend neither on the size nor on the shape and orientation of the the spatial discretization. However, the damage field depends on the spatial discretization because it is constant in each finite element. The effect of the spatial discretization is evaluated by considering both global and local criteria. The local criterion relies on the study of the crack tip position and the global one on the study of the elastic energy stored in the structure. The studied case is shown in Figure 9. The considered domain is a square that has an edge length of l = 1 m and contains a vertical crack with an initial length of a i = l/3. The initial crack is modeled by a set of broken elements (i.e., elements with d ≈ 1). The vertical displacements are set to zero on the lower boundary and the lateral boundaries are submitted to a horizontal displacement that varies linearly with the height. On the left boundary, the horizontal displacement reads u x (x = 0, y) = -cu d (1 -y/l) with c characterizing the loading level that is 0 at the initial state and 1 for the maximum loading. The horizontal displacement applied to the right boundary reads u x (x = l, y) = -u x (x = 0, y). We consider u d = 0.2 mm. The elastic parameters are E = 20 GPa for Young's modulus, and ν = 0, 25 for Poisson's ratio. The characteristic length is ℓ c = 0.1 m and the material toughness K c = 1 MPa √ m. The studied crack is submitted to a mode I loading because of the chosen boundary conditions. Crack growth is globally stable (i.e., the numerical sample failure is quasi-brittle). Figure 10 shows some of the simulation input and results, namely, two different meshes (M1 and M4, the corresponding horizontal displacement and damage fields for a loading level such that c = 1. The first three considered meshes (M1, M2 and M3) are of same shape but the element sizes are different. The length of an element edge is l/123 for the coarsest mesh, l/165 for the middle one and l/201 for the finest. For the fourth mesh, the elements are tilted at 45 • in the area where the crack is assumed to grow; the length of the element edge in that area is l/165. The horizontal displacement and damage fields shown in Figure 10 are very close for two of the meshes. For the first simulation, the length of the damaged area is equal to that of the crack, and the width of this area is equal to the element size. For the mesh containing tilted elements, the damaged area has a zigzagged shape that corresponds to the mesh shape. However the horizontal displacement field is very similar to the other ones.

To compare the results, the change of the crack length with the applied load is shown in Figure 11. The crack length a is deduced from the coordinate of the crack tip position given by the numerical model. The crack length appears to be independent of the spatial discretization. The observed differences are directly related to the time step of the calculation. For the first three simulations and for every loading step, the abscissa of the crack tip is equal to 0.5 m with an uncertainty of 1 µm. For the fourth case, the abscissa of the position of the crack tip ranges between 0.5 ± 8 × 10 -6 m for all steps. It is to be noted that the use of a tilted mesh does not significantly affect the direction of crack growth. For the simulation with the tilted mesh, the maximum distance the crack tip and the plane of symmetry is 8 µm, which represents 1.3 % of the length of an element.

In order to characterize the global response of the structure, the elastic energy of the structure at equilibrium is calculated. This quantity is evaluated instantaneously for each loading step without taking into account the structure history. The elastic energy of the structure Ω reads

E elas = 1 2 ∫ Ω ϵ ϵ ϵ : C : ϵ ϵ ϵ dΩ. ( 45
)
The change of the elastic energy with the load level is shown in Figure 12. Once again the .

results are independent of the spatial discretization. The change of the elastic energy is stairshaped because crack growth is incremental.

SIMULATION OF A CRACK NETWORK FORMATION

In this section the ability of the damage model to describe the formation of a crack network and subsequent growth is demonstrated. A rectangle of height h = 1 m and length l = 4 × h is meshed. Plane strain conditions are considered. As shown in Figure 13, the model is submitted to a vertical compressive stress of σ v = 10 MPa. A horizontal displacement, set to u x (x = 0, y) = -cu d1 (1 -y/l) + u d2 (y/l) is applied to the left side of the model. The parameter c characterizes the displacement level. It is equal to 0 for the initial state and to 1 at the end of loading. The horizontal displacement of the right side is expressed as the applied load leads to the initiation and growth of a crack network made of vertical cracks. The cracks initiate close to the lower edge and grow almost vertically. The location of crack initiation is affected by the material heterogeneity and by the regularized stress field shape. The initiation location is close to the lower edge where the principal regularized stress has the highest value, at least before the initiation of the first crack. The crack network quickly reaches a saturation state (i.e., a state in which no new cracks initiate). As shown in Figure 14(a) the first initiated crack is short. Figure 14(b) shows that the cracks that initiate just after the first one reach a length almost equal to that of the previously initiated cracks. Once saturation is reached, the degradation is only due to crack growth (Figure 14(c)). The influence of the cracks on the horizontal displacement field is observed in Figure 14(d). The cracks are not exactly straight, since the growth of each crack is influenced by the other ones.

u x (x = l, y) = -u x (x =
The change of the crack number N b with the applied load for the nine studied cases and the average crack length a are plotted in Figure 15. The results are very close for the nine considered cases (Figure 15(a)). The crack network saturation is always reached when 11, 12, or 13 cracks are formed (i.e., 11 cracks in 2 cases, for 12 in 6 cases and for 13 cracks in 1 case). For the nine considered random selections, the first initiation occurs when c ≈ 0.16, and the saturation is reached for c ≈ 0.30. It can be seen in Figure 15(b) that the cracks continue to propagate even when no new crack initiates.

SUMMARY AND PERSPECTIVES

A new nonlocal model that can be used to describe crack network formation and growth in the opening mode within a heterogeneous material is developed. Two different thresholds, both of them based on the value of regularized stress, are introduced for crack initiation and propagation. A regularization operator based on the stress gradient enhancement, which considers a characteristic length, is used to obtain a nonlocal framework.

A probabilistic damage mechanics threshold is used for crack initiation and a fracture mechanics based one for crack propagation. A natural transition between the two thresholds is performed through a 'consistency' criterion. Furthermore, the maximum size of the existing defects is considered as a 'microstructural' length that controls the characteristic length of the model. The threshold associated with crack initiation is based on a Weibull probabilistic model and takes into account scale effects and material heterogeneity. In this setting the onset of propagation of pre-existing cracks is considered as the crack initiation. The characteristic length must be greater than the size of heterogeneities to model. The crack propagation threshold is equivalent to a usual stress intensity factor criterion for small characteristic lengths compared to the crack length. The propagation threshold is expressed in terms of regularized stress by considering the stress field at the crack tip given by Westergaard's solution.

The developed method enables us to calculate the stress intensity factor in a crack network without any mesh refinement. An implicit time discretization scheme that ensures the stability of the converged states is used. The fact that the introduced regularization operator ensures mesh objectivity of the obtained results has been validated on a simple case by considering both local and global criteria. The ability of the model to describe the formation of a crack network has finally been demonstrated.

The developed model can also be extended to mode II crack propagation. Further, since a simple damage variable is introduced and because crack opening displacements are found to be in good agreement with fracture mechanics solution, the model can easily be coupled with a flow simulator to provide a hydro-mechanical damage model.

Figure 1 .

 1 Figure 1. Definition of a basis at the crack tip.

Figure 2 .

 2 Figure 2. Nodal unknowns for two types of elements.

Figure 3 .

 3 Figure 3. Numerical discretization of crack initiation (a) and growth (b).

Figure 4 .

 4 Figure 4. Schematic of the considered theoretical case.

Figure 5 .

 5 Figure 5. Stress intensity factors (a) and difference between analytical and nonlocal estimations (b) as functions of normalized crack lengths for different characteristic lengths. The proposed approximation (40) is also drawn (b).

Figure 6 .

 6 Figure 6. Change of the considered thresholds with normalized crack (or defect) size.

Figure 8 .

 8 Figure 8. Normalized mean (a) and standard deviation (b) of the initiation stress as functions of the number of elements for the numerical model and the number of random selections. The straight solid line corresponds to initiations described by a Weibull model.

Figure 9 .

 9 Figure 9. Schematic of the studied case.

Figure 10 .

 10 Figure 10. Meshes, horizontal displacement and damage fields for c = 1 and for two of the simulations corresponding to the case shown in Figure 9

Figure 12 .

 12 Figure 12. Change of elastic energy stored in the structure with a load increase from c = 0 to c = 1 for different meshes.

Figure 13 .

 13 Figure 13. Schematic representation of the studied case.

Figure 14 .Figure 15 .

 1415 Figure 14. Finite element state for different loading levels c = 0.16 (a), c = 0.25 (b), c = 1 (c) and horizontal displacement field for c = 1 (d).

  Figure 11. Crack growth with load c for the different meshes.

	1
	0.8
	0.6
	0.4
	0.2
	0
	0	0.2	0.4	0.6	0.8	1
	Copyright c ⃝ 2000 John Wiley & Sons, Ltd.			Int. J. Numer. Meth. Engng 2000; 00:1-6
	Prepared using nmeauth.cls

Copyright c ⃝ 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 00:1-6 Prepared using nmeauth.cls

ACKNOWLEDGEMENTS

This work was supported by BRGM through an 'Institut Carnot' research fund. The authors wish also to thank the Code Aster development team, in particular Mr. Clement Chavant and Dr. Sylvie Granet from EDF R&D for their invaluable help in implementing the model.