Thanh-Tung Nguyen

José Bolivar

Julien Réthoré
email: julien.rethore@ec-nantes.fr

Marie-Christine Baietto

Marion Fregonese

A Phase field method for modeling stress corrosion crack propagation in a nickel base alloy

Keywords: Stress Corrosion Cracking, Crack propagation, Phase field, Digital Image Correlation, Slip dissolution, Crack coalescence, Complex crack networks

Stress Corrosion Cracking (SCC) is a very common failure mechanism characterized by a slow, environmentally influenced crack propagation in structural components. The mechanisms proposed to explain, at the microscopic scale, the cracking propagation processes are not able to elucidate all aspects of this phenomenon in different metal/environment systems. This work is concerned with the development of a new multiphysics model for understanding the phenomena of crack propagation under the effect of SCC. This new model is based upon: (i) a Phase field method, based on a variational formulation of brittle fracture with regularized approximation of discontinuities; (ii) a robust algorithm capable to prescribe the displacements (over the boundary of a small sub-volume) and crack onset obtained by image processing based on digital image correlation in the sample during the numerical simulations; (iii) a coupling with a diffusion model informed with first-principles

Introduction

Predicting the strength and durability of engineering components and structures using numerical simulations of fracture phenomena induced by Stress Corrosion Cracking (SCC) is a very challenging problem. SCC can often be regarded as the result of localized oxidation enhanced by stress/strain: a combined action of a sustained mechanical loading and a chemically aggressive environment. The occurrence of SCC in a structure can lead to catastrophic failures, it has thus long been recognized as potentially dangerous. One of the mechanisms responsible for this process noted in the literature is the embrittlement of the material due to local anodic dissolution of fresh metal created by slip band emergence at the crack tip. In this area recent developments and corresponding issues have been reviewed in [START_REF] Newman | Stability, validity, and sensitivity to input parameters of the slip-dissolution model for stress-corrosion cracking[END_REF], followed by some original contributions related to different aspects of the problem: loading rate effects, harmonization of different formulations of the model, and sensitivity of the model to various parameters. Developing numerical models for the quantitative evaluation of stress corrosion cracking growth rates in terms of key engineering parameters is highly required.

Various experimental methods have been proposed in the literature to study cracks induced by SCC. The measurements of corrosion contribution is usually based on electrochemical methods such as Electrochemical Noise (EN) [START_REF] Hladky | The measurement of localized corrosion using electrochemical noise[END_REF]; [START_REF] Kearns | Electrochemical noise measurement for corrosion applications[END_REF]; [START_REF] Kovac | Correlations of electrochemical noise, acoustic emission and complementary monitoring techniques during intergranular stress-corrosion cracking of austenitic stainless steel[END_REF]. Thanks to its capacity of distinguishing between different corrosion types, the (EN) method is mostly applied to detect and evaluate SCC activity. Acoustic Emission (AE) [START_REF] Bellenger | Use of acoustic emission technique for the early detection of aluminum alloys exfoliation corrosion[END_REF]; [START_REF] Jomdecha | Study on source location using an acoustic emission system for various corrosion types[END_REF]; [START_REF] Shaikh | Evaluation of stress corrosion cracking phenomenon in an aisi type 316ln stainless steel using acoustic emission technique[END_REF] is also widely used because of its ability to detect material deformation, cracking and fracture. Many microscopy techniques have been applied

to characterize SCC at different scales. For example, magnetic force microscopy (MFM) has been used to detect the chromium depleted regions of type 304 stainless steel [START_REF] Takaya | Estimation of stress corrosion cracking sensitivity of type 304 stainless steel by magnetic force microscope[END_REF], the early stages of pitting and localized corrosion have been studied by combining in situ observations and electromechanical atomic force microscopy (EAFM) [START_REF] Williford | In situ observations of the early stages of localized corrosion in type 304 ss using the electrochemical atomic force microscope[END_REF], and scanning electron microscopy (SEM) is used to study the pitting corrosion of Al [START_REF] Richardson | A study of the pitting corrosion of al by scanning electron microscopy[END_REF]. This technique is also used to observe the fracture surface of SCC in high strength steel fasteners [START_REF] Abhay | Stress corrosion cracking of high strength steel fasteners for space applications[END_REF]. The Diffraction Contrast Tomography (DCT) is another technique used to characterize SCC, for example the work of [START_REF] King | Observations of intergranular stress corrosion cracking in a grain-mapped polycrystal[END_REF] study the interaction between intergranular stress corrosion cracking and microstructure in a Grain-Mapped Polycrystal. In the work of [START_REF] Babout | X-ray microtomographic observation of intergranular stress corrosion cracking in sensitised austenitic stainless steel[END_REF], the phenomenon of intergranular stress corrosion cracking in a sensitised type 302 stainless steel wire has been observed in situ using high resolution X-ray microtomography. An important step after image acquisition by OM/SEM is the processing of microstructural images to obtain a qualitative/quantitative analysis of the sample. Many techniques have been developed in the literature with this aim. Among these techniques, digital image correlation (DIC), or digital volume correlation (DVC) in 3D, is one of the most efficient and versatile tool. This procedure is based on the comparison of images acquired at different stages of a mechanical test and provides quantitative descriptions of local deformation. DIC could be used instead of, or in addition to the prevailing methods to detect defects, by observing the singularity due to the change in the strain distribution, or a step jump in the displacement field across the cracks.

In the literature, most of the models proposed to predict SCC rates are based on linear elastic fracture mechanics concepts without deeply considering corrosion kinetics. A numerical simulation scheme of SCC is developed in the work of [START_REF] Choi | Stress corrosion cracking in plastic pipes: Observation and modeling[END_REF] for a thermoplastic material based on crack layer theory. SCC is here considered as a superposition of creep, it induces aging and chemical degradation. Then, they use the phenomenological power laws to describe the loss of material toughness due to these two mechanisms.

In addition, this work proposed a generic kinetic parameter to represent the rate of chemical corrosion. Both SCC and mechanically driven crack growth were simulated. In the work of [START_REF] Saito | Mechanochemical model to predict stress corrosion crack growth of stainless steel in high temperature water[END_REF], a predictive methodology for SCC crack growth using a mechanochemical model based on a slip formation/dissolution mechanism is presented. The mechanochemical model consists of the combined kinetics of the plastic deformation process as a mechanical factor and the slip dissolution-repassivation process as an environmental factor at the crack tip. The role of the passive film has also been studied in many works, its physical degradation is usually expressed by a unique parameter namely the "crack tip strain rate". Many works in literature use this concept that is applied to study SCC in many cases, such as those by [START_REF] Vermilyea | A theory for the propagation of stress corrosion cracks in metals[END_REF]; [START_REF] Vermilya | Concerning strain-enhanced corrosion mechanisms of scc[END_REF], [START_REF] Parkins | Predictive approaches to stress corrosion cracking failure[END_REF][START_REF] Parkins | Current topics in corrosion: Factors influencing stress corrosion crack growth kinetics[END_REF], [START_REF] Ford | Development and use of a predictive model of crack propagation in 304/316 l, a 533 b/a 508 and inconel 600/182 alloys in 288 deg c water[END_REF],... However, the results of these models are rarely compared with experimental ones, especially concerning crack morphology. With recent advances in numerical simulation methods new studies are now possible, allowing the development of fracture models at microscale, in particular for complex cracks morphology including phenomena of crack initiation, coalescence and propagation. The objective of this proposed paper is to construct a numerical model to simulate the failure of a nickel base alloy under SCC in an acidified tetrathionate environment. The reactiondissolution process is simulated with a kinetics model, and integrated into a mechanical analysis. The anodic dissolution, that is involved in the decrease of fracture property is considered to be dependent on the extent of corrosion, and such dependence is implemented within a fracture model. For this purpose, we use the phase field method based on the variational formulation of the crack evolution problem due to [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. The regularized setting of their framework has been considered in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]. Another contribution for the phase field method based on LandauG-inzburg type phase-field evolution equations can be found in [START_REF] Karma | Phase-field model of mode iii dynamic fracture[END_REF]. The phase field method may be approached both theoretically and practically, providing an excellent framework for taking into account environmental effects on the failure phenomenon. This framework is also extremely suitable for establishing parametric trends besides enabling the assessment of safety margins. This work is concerned with the development of a new multiphysics model for simulating the phenomena of crack propagation under the effect of SCC based on the phase field method and the numerical implementation proposed by [START_REF] Miehe | A phasefield model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. Direct comparisons between the simulationed crack paths and experimental data are performed. In addition, a strategy is proposed to identify the parameters of the constitutive model by using DIC results.

The overview of the paper is as follows. A brief introduction of experimental methods is described in section 2. Then, we present a short review about the slip dissolution model for SCC in section 3. The new model proposed for SCC propagation based on the phase field method is presented in section 4. The proposed algorithm is described in section 5. Finally, the method is evaluated and illustrated by damage benchmarks and practical examples involving crack onset and propagation provided by SCC, and then we present the direct comparisons between the model and the experiments in section 6.

Experimental method for observation of crack propagation

Experimental procedures

This study focuses on Inconel 600 alloy material. Its chemical composition is given in Table 1. The presented experiments were performed on specimens with a gauge dimension of 117 mm length, 8 mm width and 2 mm thickness as depicted in Fig. 1(a). (see [START_REF] Bolivar | Evaluation of multiple stress corrosion cracks interactions by in-situ digital image correlation[END_REF] for more details of sample preparation)

Elements Ni

Fe Cr Mn Cu Co Ti-Al C S Wt. Pct. bal. 9.25 15.52 0.12 0.1 0.1 0.25 0.03 0.002 To use DIC, a specific procedure has been proposed to improve the grey level contrast on the sample surface. Finally, the samples exhibit many random heterogeneities that give local contrast of grey levels in the images, what is promising for DIC. The details of experimental setup has been described

in [START_REF] Bolivar | Evaluation of multiple stress corrosion cracks interactions by in-situ digital image correlation[END_REF] Constant load tests were performed for this study. First, the samples are immersed in the test solution (10 mM K 2 S 4 O 6 , pH = 3) for about 45 minutes.

Then, the specimens are subjected to a constant applied load corresponding to 80 % of the yield stress (263 MPa) by using an electromechanical tensile machine (LLOYD LR30K Plus).

The sample surface was observed in-situ through a 2/3" CCD Stingray ton et al. (1983) and this technique is currently used for many applications.

We refer to references [START_REF] Sutton | Book-image correlation for shape, motion and deformation measurements[END_REF] for a detailed description of this technique. The specific 2D-DIC procedure used in this work is similar to the one used in reference [START_REF] Besnard | finite-element' displacement fields analysis from digital images: Application to portevin-le châtelier bands[END_REF].

The DIC analysis is performed using the in-house software. The procedure is based on four main steps: • DIC analysis: the correlation analysis is performed. We use here a mesh size of 16 pixels (55 µm). After this step, we obtain the displacement fields for series of deformed images (see Fig. 2(a)).

(a) (b) (c)
• Postprocessing displacement fields: crack detection by using only displacement field is not straightforward see Fig. 2(a). A method has been proposed to reduce the noise and increase the contrast due to displacement jumps. This method is based on the removal of Rigid Body Movement (RBM) and the displacement due to the overall homogeneous strain. A median filter is subsequently applied to the resulting displacement to remove the fluctuation that doesn't have a discontinuity feature. Then, the element-wise displacement variations in the loading direction are calculated. The obtained result is depicted in • Crack extraction: a segmentation procedure is performed to extract the cracked area from the post-processed displacement field. In this example, we perform the extraction of cracks by using the classical threshold criterion for the post-processed displacement field (threshold value 0.45 µm). The obtained result is depicted in Fig. 2(c).

• Parameters identifications: Finally, cracks are labeled and followed during the image sequence. We also identify several fracture parameters such as: crack length, crack growth rate, crack initiation...(see [START_REF] Bolivar | Evaluation of multiple stress corrosion cracks interactions by in-situ digital image correlation[END_REF] for more details)

In this work, DIC crack extraction results will be compared with those of the numerical simulation. For this purpose, an algorithm is proposed to construct input modeling from these data, this procedure will be introduced in Section 6.2.1.

Slip dissolution model for stress corrosion cracking

One of the models that successfully rationalizes the phenomenology and the kinetics of stress-corrosion cracking (SCC) in this situation, is the slip dissolution model (also called the film rupture or slip-oxidation model). In this model, the propagation of cracks is mainly due to the rupture of a passive film at the crack tip. In the initial state, the passive film forms over the crack faces and at the crack tip. SCC results from the dynamic synergy between passive film formation, crack tip strain rate, depassivation, repassivation...The slip dissolution model postulates that the rupture of passive film at the crack tip is due to the local dynamic plastic strain; this process allows the dissolution of metal in this region and provides the advance of cracks.

The main ideas of this model are illustrated in Figs. 3,4. Crack velocity in this case is directly related to the rate of metal dissolution. However, this assumption leads to consider another aspect that is repassivation of the crack faces. The works of [START_REF] Scully | The theory of stress corrosion cracking in alloys[END_REF][START_REF] Scully | The current situation in stress-corrosion studies[END_REF][START_REF] Scully | The interaction of strain-rate and repassivation rate in stress corrosion crack propagation[END_REF] proposed that crack propagation can only be observed when re-passivation rate and crack tip strain rate are close.

The slip dissolution model can be divided into three steps: (S1) the crack tip is covered by the passive film and the loading is applied; (S2) the rupture of the passive film occurs when the strain rate reaches a threshold value, then metal dissolution at the crack tip takes place and crack propagation occurs;

(S3) the passive film is reformed. The detail of these steps are illustrated in Four main "postulates" can be enhanced and re-inserted in the model that will be developed in this work. (i) the strain rate can be considered as the rate-limiting step of SCC in this system due to the rupture of the passive film at the crack tip (P1). At higher strain rate, the time available for corrosion process to develop is too short and the ductile fracture remains predominant. At lower strain rate, the passive film at crack tip doesn't fail and dissolution of the metal at the crack tip is inhibited. So, ductile failure due to insufficient corrosion is observed.

At a smaller scale, SCC is related to the micro-structure of the material as dissolution occurs preferentially within specific features, for example along the grain boundaries (see Fig. 3,4). Because they are areas of high energy, many impurities segregate to them, and they are a preferential site for precipitation of carbides. However, in this work, the macroscopic scale (the scale of the homogeneous material) is used. Hence, (ii) only the SCC at the crack tip will be considered (P2).

On the other hand, the experimental observations show that, SCC cracks present a ratio between crack length and crack opening that can be found within the range 100 to 1000, some time above 1000. This result means that, when the crack propagates, the crack faces are expected to show neither corrosion nor plastic deformation. Moreover, we observe no stress appearing in these regions (crack faces), whereas a high-stress concentration is found at the crack tip. Thus, (iii) the environment is in contact with a material over a short distance close to the crack tip (P3).

In addition, one has to make a difference between material points along the crack front that are located on the specimen surface (in direct contact with the solution) or in the bulk material. It provides (iv) the different in-fluences of SCC on the crack tip with respect to the depth (P4). This aspect have to be accounted for in the model.

Modeling

Description of film rupture phenomena by using diffusion of degradation

In this subsection, we propose a model to describe the influence of film rupture on the fracture problem. As discussed above, the major impact of slip dissolution is the change of the material behavior. In more details, material embrittlement occurs and reduces the resistance of material. To introduce these phenomena in the modeling, we use here a diffusion model to describe the evolution of the softening of material around the crack tip. In the activated crack, i.e when the crack tip strain rate is higher than a critical value, the passive film fails. Then, with the presence of chemical products and SCC conditions in this region, the dissolution of the material occurs. However, modeling these phenomena by using numerical simulations is not straightforward, so we use here another way: the dissolution of material is replaced by the diffusing degradation c. In addition, the criterion of strain rate ε (see P1) is replaced by the crack propagation rate ḋ because the definition of strain at the crack tip is a difficult question.

The main idea of the proposed model is depicted in Fig. 5. We describe here the equivalent meaning of the two models that yields a pronounced softening behavior of material at the crack tip (due to the slip dissolution or the appearing of "micro-cracks").

The distribution of the diffusing degradation is described by using the Fick's laws of diffusion. The first law states that the diffusing degradation flux q is proportional to the diffusing concentration gradient. q represents the source of SCC (slip dissolution in the crack tip) that induces the propagation of "micro-cracks" or material degradation. In mathematical terms this can be expressed as:

q = -D(x)∇c (1)
where D(x) is the diffusion coefficient and ∇ is the gradient operator.

Fick's second law relates the diffusing flux q at a given point to its degradation rate by means of the equation

∂c ∂t = -∇q (2)
We apply this principle only for activated crack tips (see P2,3), i.e. we impose c = 1 at the crack tip by using Dirichlet boundary conditions and then the diffusion problem is solved for the whole sample. Note that the crack tips are updated for each time steps corresponding to crack propagation.

Following (P4), the influence of film rupture phenomena in the surface and within the volume are different, that provides the tendency of crack propagation preferentially on the surface. To model this effect, we use here an in-homogeneous diffusion coefficient, i.e depending on the depth. We present an example of a structure containing a semi elliptical crack (see for the case of structure containing a semi elliptical crack this case). From these properties, we can construct a simple solution for this example:

           D(x) = D 0 ∀x ∈ ∂Ω x=0 D(x) = D 0 B -x B m ∀x ∈ Ω/∂Ω x=0 (3)
Where the coefficient m 2

For an in-homogeneous isotropic medium, (D(x) depends only on the position x) the diffusion equation is rewritten as

∂c ∂t = D(x) c + 3 i=1 ∂D(x) ∂x i ∂c ∂x i (4)
For more complex sample geometries, a levelset function (the distance to the surface in contact with the aggressive medium) is used to compute the value of D at any point inside the sample. To model macroscopic material embrittlement, the fracture toughness g c directly depends on the diffusing degradation. The increasing of c induces the decreasing of g c . However, when full diffusing degradation occurs (c = 1), the material is still resistant, thus g c should not vanish for c = 1. We define the relation between g c and c

by the following expression:

g c = (1 -k c ) (1 -c) n + k c .g 0 c (5)
An example of this typical function is depicted in Fig. 7, for several values of n = 2, 3, 6. We use here k c = 10%, i.e for c = 1, 10% of the material resistance remains.

Phase field modeling of stress corrosion cracking

In the following, the basic concepts of the phase field method are briefly summarized. For more details and practical implementation aspects, the interested reader may refer to [START_REF] Miehe | A phasefield model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]; [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF]. The phase field method is based on a regularized formulation of a sharp crack description. A regularized variational principle describing both the evolution of the mechanical problem and of an additional field d describing the damage (called phase field), is discretized by a finite element procedure and a staggered algorithm. The method alleviates the shortcomings of re-meshing crack geometry by using a fixed mesh and a regularized description of the discontinuities. In addition, crack initiation can be modeled in a straightforward manner. In contrast to volume damage models, usually implemented in non-linear codes, such regularized approach is directly connected to the brittle failure theory of crack propagation. In the present work, the phase field method has been implemented in a house-made code both in 2D and 3D.

In the phase field method, assuming small strains, the regularized form of the energy describing the cracked structure is expressed by:

E(u, d) = Ω W u (ε(u), d)dΩ + Ω g c (x, τ )γ(d, ∇d)dΩ, (6)
where ε is the linearized strain tensor, W u is the elastic strain energy density, depending on the displacements u(x) and on the phase field d(x) describing the damage of the solid. g c is a function of the fracture toughness at initial state g 0 c and the diffusing degradation c. Note that, the irreversibility condition is applied for the fracture toughness, by the following expression:

g c (x, t) = min τ ∈[0,t] g c g 0 c , c(x, τ ) . ( 7
)
For the sake of simplicity, we write g c = g c (x, τ ). In ( 6

γ(d, ∇d) = 1 2l d 2 + l 2 ∇d • ∇d, ( 8
)
where l is the regularization parameter. l can be considered as a pure numerical parameter of the regularized model of brittle fracture or seen as a real material parameter for a gradient damage model. In the first case, it is recommended to take l as small as possible to better approximate brittle fracture, with regards to the size of the mesh. In the second case, l should be identified from experimental data. Such analysis with experimental validations has been recently published in [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF]. It is demonstrated that this kind of regularized model converge to the brittle failure theory when l vanishes. In the case considered herein when g c is heterogeneous, this convergence is maintained if l is smaller than the shortest wave length of the g c field. It will checked latter on that this condition is satisfied and thus that the ability of the phase field model to approximate brittle failure is effective.

The total energy is then rewritten as E = Ω W dΩ with

W = W u (ε(u), d) + g c γ(d, ∇d) (9)
The unilateral contact formulation following the work of Miehe et al.

(2010) is chosen to ensure damage induced by traction only, (assuming isotropic elastic behavior of the body) through:

W u (u, d) = Ψ + (ε(u)) {g(d) + k} + Ψ -(ε(u)) . ( 10
) With Ψ + (ε) = λ 2 T r(ε) + 2 + µT r ε + 2 , ( 11
) Ψ -(ε) = λ 2 T r(ε) - 2 + µT r ε -2 , ( 12
)
where ε + and ε -are the extensive and compressive modes of strain field ε = ε + + ε -. In ( 11), ( 12) x + = (x + |x|) /2 and x -= (x -|x|) /2. The degradation function g(d) in Eq. ( 10) is assumed to have the simple form:

g(d) = (1 -d) 2 + k. ( 13
)
The function g(d) has been chosen such that g (1) = 0 to guarantee that the strain energy density function takes a finite value as the domain is locally cracked (see e.g. [START_REF] Braides | Approximation of Free Discontinuity Problems[END_REF]) and g(0) = 1 to guarantee that the material is initially undamaged. g(1) = 0 is the limit for fully damaged material. The small parameter k << 1 is introduced to maintain the well-posedness of the system for partially broken parts of the domain [START_REF] Miehe | A phasefield model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF].

In the following, a crack phase field evolution law is derived that can guarantee the irreversibility of the process. Assuming isothermal processes, without external micro forces, the reduced form of the Clausius-Duhem inequality is written as (see [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF] for more details):

A ḋ ≥ 0, (14)
where A = -δW δd = -∂W ∂d + ∇. ∂W ∂∇d is the variational derivative of W with respect to the phase field d. At this stage, a threshold function F (A) such that

F (A) ≤ 0 ( 15
)
is introduced. Assuming the principle of maximum dissipation then requires the dissipation A ḋ to be maximum under the constraint (15). Using the method of Lagrange multipliers and the following Lagrangian:

L = -A ḋ + λF (A), (16)
yields the Kuhn-Tucker equations:

∂L ∂A = 0, λ ≥ 0, F ≤ 0, λF = 0. ( 17
)
The first equality in (17) gives:

ḋ = λ ∂F (A) ∂A . (18)
Without loss of generality, the threshold function F (A) is assumed in the form F (A) = A. From (18) and using the second inequality in (17), we obtain:

ḋ = λ ≥ 0 (19)
For ḋ > 0, and from ( 14), ( 19) and the third equality in (17), which give F = 0, implying:

F = - δW δd = - ∂W u (u, d) ∂d -g c δγ(d, ∇d) = 0. ( 20
)
With (see e.g. [START_REF] Miehe | A phasefield model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]):

δγ(d, ∇d) = d l -l∆d (21)
From ( 10), ( 13) and ( 20) we obtain

2(1 -d)Ψ + -g c δγ(d, ∇d) = 0 (22)
that is the evolution law for the phase field d. As 2(1 -d)Ψ + ≥ 0, then

δγ(d, ∇d) ≥ 0. ( 23
)
We can thus check that due to (23), the variation of crack length:

Γl = Ω δγ(d, ∇d) ḋdΩ ≥ 0, ( 24
)
satisfying irreversible evolution of cracks.

The strain energy density function H(t) is introduced to describe a dependence on history [START_REF] Miehe | A phasefield model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and to make possible loading-unloading.

In this work, this function can be defined by the following:

H(x, t) = max τ ∈[0,t] l g c Ψ + (x, τ ) , (25)
which is substituted to Ψ + in (22) and using (21) it yields the following phase field problem to be solved to evaluate the field d(x, t) at time t:

         2(1 -d)H -{d -l 2 ∆d} = 0 in Ω, d(x) = 1 on Γ, ∇d(x) • n = 0 on ∂Ω, (26)

Numerical formulation and algorithm

Diffusion problem

The strong form defined in (2) can be expressed as:

∂c ∂t + ∂q i ∂x i = 0 (27)
where q i represents the components of the flux vector expressed in an orthonormal basis with coordinates x i . The corresponding boundary conditions are given by

         c(x) = c with x ∈ ∂Ω 1 q i (x) = q i with x ∈ ∂Ω 2 ∂Ω 1 ∪ ∂Ω 2 = ∂Ω and ∂Ω 1 ∩ ∂Ω 2 = ∅ (28)
Multiplying the strong form by a test function δc and integrating over Ω, we obtain:

Ω δc ∂c ∂t + δc ∂q i ∂x i dΩ = 0. ( 29
)
where δc satisfy:

Ω ∂δc ∂x i dΩ < ∞ and δc = 0 on ∂Ω 1 (30)
Using the integral by parts and the divergence theorem, we finally obtain:

Ω δc ∂c ∂t - ∂δc ∂x i q i dΩ + ∂Ω δcq i n i dS = 0. (31)

Phase field problem

Starting from (26) 1 , multiplying by a test function δd and integrating over Ω, we obtain:

Ω 2(1 -d)Hδd -d -l 2 ∆d δd dΩ = 0. ( 32
)
Using the property:

(∆d) δd = ∇ • (∇dδd) -∇d • ∇(δd) (33)
and the divergence theorem, Eq. ( 32) is rewritten as:

Ω {2(1 -d)H -d} δd - Ω l 2 ∇d • ∇(δd) dΩ + ∂Ω l 2 ∇d • nδd dΓ = 0 (34)
Using ( 263 ), we finally obtain:

Ω (2H + 1) dδd + l 2 ∇d • ∇(δd) dΩ = Ω 2HδddΩ. ( 35
)
In the present work, the computations are performed in quasi-static conditions. Then, the time steps introduced in the following actually refer to load increments. Introducing a time stepping, the problem to be solved at time t n+1 is expressed by seeking d(x) ∈ S d , such that:

Ω (2H n + 1) d n+1 δd + l 2 ∇d n+1 • ∇(δd) dΩ = Ω 2H n δddΩ. ∀δd(x) ∈ H 1 0 (Ω). ( 36
)
where H n = H(u n ) and the fracture toughness g c n are computed from the previous time step (load increment) by:

   g c n (x) = g c n (g 0 c , c n (x)) if g c n (x) ≤ g c n-1 (x), g c n (x) = g c n-1 (g 0 c , c n (x)) if g c n (x) > g c n-1 (x). (37)
and

           H n (x) = l g c n Ψ + n (x) if Ψ + n (x) > Ψ + n-1 (x), H n (x) = l g c n Ψ + n-1 (x) if Ψ + n (x) ≤ Ψ + n-1 (x). (38)
Note that Eqs. ( 38),( 37) is the algorithmic counterpart of Eqs. ( 25), (7).

Displacement problem

The weak form associated with the displacement problem is found by solving the variational problem:

u(x) = Arg inf u∈Su E(u, d) -W ext (39)
where

S u = {u|u(x) = ū on ∂Ω u , u ∈ H 1 (Ω)} and W ext = Ω f • udΩ + ∂Ω F F
• udΓ with f and F body forces and prescribed traction over the boundary ∂Ω F . We obtain the classical weak form for u(x) ∈ S u :

Ω σ : ε(δu)dΩ = Ω f • δudΩ + ∂Ω F F • δudΓ ∀δu ∈ H 1 0 (Ω), (40)
where the second-order Cauchy stress tensor σ = ∂Wu ∂ε is given using ( 10) and ( 13), by:

σ = (1 -d) 2 + k λ T rε + 1 + 2µε + + λ T rε -1 + 2µε -.
(41)

Crack onset by using experimental data

Using the experimental procedures introduced in section 2, we present here a typical result in Fig. 8, six steps are depicted. We recognize the crack onset at different points in the sample surface. Faults of micro-structure lead to stress concentrations in the surrounding region, in which the passive film fails and provides new crack nucleation sites. Then, crack merging is observed, and finally a continuous crack network is created, the coalescence process becomingxs an essential characteristics.

In this work, we study the fracture phenomena at meso scale. The sources (for example: the fault of micro-structure) of these phenomena come from The main idea of the method is the following: at a loading step n, experimental results are used to initiate cracks. Then, cracks are free to propagate by themselves. We present here an example in Fig. 9. At loading step 1, the structure is subjected to a prescribed displacement U 1 . This time step shows a crack (black color) and two crack onsets, which are depicted by two blue squares in Fig. 9 The initial crack can be modeled as either a discrete crack in the geometry or as an induced crack in the phase-field. However, the crack initiation is promoted by a micro-structurally induced stress concentration. When crack initiation is detected experimentally from DIC data, stress concentrations have already played their role of crack initiation precursor. By using the discontinuous geometry, the discrete crack always requires a period of time (several loading steps) to induce stress concentration, and then propagate. The propagation of cracks obtained in the simulations is thus delayed if discrete cracks are introduced at the time when they are detected experimentally.

We propose here a solution to treat this issue, the crack is initiated at earlier state (for example: nucleation) by using the strain history function defined in Eq. ( 25). This method has been introduced in the work of [START_REF] Borden | A phasefield description of dynamic brittle fracture[END_REF] to initiate cracks. We use this method to nucleate cracks (d = 0.25, see [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]) in the numerical model before they can be detected experimentally.

For the induced crack nucleation, an initial strain history field is specified such that an initial crack in the phase-field is defined. To define the initial strain-history field we let L be a line that represents the discrete crack we wish to include and DisT (x, L) is the distance from x to the line L. The strain history field is then defined as:

H 0 (x) =            d 0 2(1 -d 0 ) 1 - DisT (x, L) l for DisT (x, L) ≤ l 0 for DisT (x, L) > l (42)
where d 0 is the initial prescribed value of the phase field. We can use this principle to initiate cracks anywhere in the domain without reference to the mesh. This method proves highly advantageous in specifying complex surface cracks in 3D bodies. We also set the intial state of crack by defining a corresponding value for d 0 which is here set to the crack nucleation threshold, i.e d 0 = 0.25. Crack nucleation introduces stress concentration, providing a source to crack onset later. This mechanism is similar to stress concentrations induced by the fault in micro-structure.

However, this method requires a new parameter, that is the delay between the time to nucleate cracks and the time when they are detected from DIC data. We define ∆t as the delay between crack nucleation and crack onset (observed by DIC). This parameter will be determined by using inverse analyzes (introduced in the next section). Moreover, in order to account for experimental observations (rupture of passive films at a region in the surface provides the stress concentration and then crack nucleation), cracks are nucluated within an finite element in the surface, and we immediately activate SCC by means of the diffusing degradation i.e. prescribing c = 1 at the surface nodes of this element.

Algorithm

The overall algorithm, involving the coupled SCC -damage mechanics model, is described as follows:

Initialization: Initialize the displacement field u 0 (x), the phase field d 0 (x), and the degradation field c 0 (x).

FOR all loading increments (pseudo time t n+1 ):

Given u n (x), d n (x) and c n (x):

1. Compute phase field d n+1 (x)

1.1 Update fracture toughness g c (x, τ ) from ( 5) and ( 7

END

In this algorithm, the overall coupled problem is solved in a staggered manner. This is made possible by the use of the history function in the phase field problem. This strategy avoids iterating between the displacement and the phase field problems but the load increment has to be small enought to ensure converge of in the space-time domain.

Results and discussion

6.1. Benchmark test of a structure containing a semi elliptical crack

The main purpose of this first example is to consider the influence of SCC on fracture mechanics. For this aim, we consider the benchmark problem of the propagation of a semi elliptical crack. A domain whose edge lengths are

L × B × H = 1.0 × 0.4 × 1.0 mm 3 contains an initial crack (semi elliptical
form), as depicted in Fig. 10. The lower surface (z = 0) of the domain is constrained along z-directions, while the displacement along x, y-directions are free. On the upper surface (z = H), the uniform z-displacement U increases with time (the displacement along x, y-directions are also free). Due to this tensile loading, the semi elliptical crack propagates. A slice of the mesh in the plane of the semi elliptical crack is depicted in Fig. 10(c). The typical size of an element in the crack propagation zone is

h min ≈ 5 µm.
In the rest of the domain the typical mesh size is h max = 0.05 mm.

The material is supposed to be homogeneous, elastic and isotropic with properties E = 210 GPa and ν = 0.3 GPa. The fracture toughness is g 0 c = 100 N/m, and the length scale parameter is chosen as l = 0.01 mm to satisfy the condition h min ≤ l/2. Monotonic tensile displacement increments of ∆U = 8.10 -5 mm have been prescribed as long as d < 0.9 in all elements. To ensure, the convergence of the staggered scheme, we use ∆U = 1.10 -5 mm as soon as d > 0.9 at one integration point.. The evolution of the crack during the simulation is shown in Fig. 11(a In the following, SCC will be considered through the proposed model.

The parameter for diffusion is chosen as a monotonically increasing function of depth:

D(x) = D 0 B -x B 6 (43)
The influence of diffusing degradation on the the fracture toughness is defined by the following expression: The crack evolution predicted by the two models are strongly different (see Fig. 12). In the case of classical phase field model (purely mechanical), the propagation tendency agrees with the theory of semi-elliptic surface crack growth [START_REF] Newman | An empirical stress-intensity factor equation for the surface crack[END_REF]; [START_REF] Smith | A finite element model for the shape development of irregular planar cracks[END_REF]; [START_REF] Kikuchi | Surface crack growth simulation under mixed mode cyclic loading condition[END_REF]. Depending on the ratio of crack depth to crack length, the crack first propagates to obtain a reasonable ratio, then propagates in both directions (depth and surface) in order to keep this ratio. In the case when SCC is considered, we observe that the crack propagates in the surface more than in depth. We recognize that, the propagation tendency of crack in the SCC model is strongly dependent on the diffusing direction of degradation. The evolution of the diffusing degradation on a plane of observation is plotted in Fig. 13. Due to the expression adopted for D (see Eq. ( 43)), the degradation diffuses more along y-direction than along x-direction. This result reproduces very well the phenomena mentioned in experiments in the literature, that is promising for direct comparisons of cracks morphology between experiments and numerical simulations. As mentioned above, in this frame work, the SCC is accounted by using a spatially varying fracture energy g c , depending on the diffusing degradation.

g c = (1 -k c ) (1 -c) 3 + k c .g 0 c (44)
To see better this phenomenon, in Fig. 14 we plot the variation of g c along a line of investigation located at the intersection of two planes z = 0.5 and y = 0.5. The fracture energy reaches the minimum value g min c = k c g 0 c = 10 -5 kN/mm at crack tip and the maximum value g min c = g 0 c = 10 -4 kN/mm away from crack tip. The variation take place in a region about of 0.15 mm that is much higher than the regularization parameter l = 0.01 mm. This ensures that the two diffusion problems, governing the evolution of c (or g c ) and d have well separate characteristic length. The Fig. 15 plots the surface energy along the same line of investigation, for one loading step when damage has not yet occurred and three cracked steps (depicted in Fig. 13). This figure allows to check that there is no interaction between the variation of g c which is very smooth as depicted on Fig. 14 and the phase field evolution that is related to the distribution of dissipated energy.

The stress-displacement curves obtained for both cases are depicted in Fig. 18. SCC crack propagates at a stress level below the one observed for classical case (without SCC).

In the following, we illustrate the performance of the new model in the case of constant load tests as usually performed experimentally to study SCC.

Neumann boundary conditions are used in this test. The same structure as in the previous example is loaded until 80% of the critical stress, and then the load is kept constant. Then, SCC is activated. Because of the diffusing degradation, that reduces the resistance of the material, the crack propagates. The result is depicted in Fig. 17

Experimental validation

In this section, the crack network obtained from the experiments and the one predicted by the simulations will be compared. The proposed procedure to identify some parameters of the constitutive relations of the model is presented. This procedure is based on the inverse analysis combining experiments and simulations.

Identification of model parameters by inverse analysis

The elastic parameters have been taken from literature (e.g. [START_REF] Mills | On the relationship between stretch zone formation and the j integral for high strain-hardening materials[END_REF]; [START_REF] Mills | Fracture toughness of alloy 600 and an en82h weld in air and water[END_REF]; Al-Rubaie et al. ( 2007)). We have E = 210 GPa and ν = 0.3, the fracture toughness at initial state (without SCC) is g 0 c ≈ 300 N/m, the regularization length is considered as a pure numerical parameter, it is chosen corresponding to mesh size l = 0.11 mm. To perform the numerical simulations, some parameters are still to be defined. The first one is the diffusing degradation coefficient D 0 in Eq. ( 3), (we assume m = 6). The second parameter is the coefficient n in Eq. ( 5) (we use here k c = 1%). Another parameter is ∆t, the delay between crack nucleation and crack onset, as discussed above in Subsection 5.4. For the sake of simplicity, the identification is restricted to a comparison between experimental and simulation results within a small region. The main idea is described in Fig. 19. At time step s, the experimental data allows us to know, the displacement prescribed on the boundary U s and the initiation of several cracks, as depicted in Fig. 19(b). We perform a numerical simulation with the boundary conditions defined by DIC data until time step s. The regular mesh used for DIC is used to built a 3D mesh of 8-node elements. Note that, the crack is nucleated at time step s -∆t. In addition, the displacement field of DIC clearly shows fluctuations due to noise (see The cost function to minimize for identifying the model parameters is defined as: are the phase field values obtained for given D i 0 , n i , ∆t i at the iterations i. We perform this minimization over several regions, and then take the average of the converged parameter values, we obtain: D 0 ≈ 3.32 × 10 -6 mm 2 /s, n ≈ 3 and ∆t ≈ 4800s.

f (D i 0 , n i , ∆t i ) = N k d exp k -d num k (D i 0 , n i , ∆t i ) 2 (

Surface comparison of crack morphology between numerical simulations and experiments

In the following, three sub-volumes are considered for a detailed comparison between the experimentally observed cracks from DIC processing and the numerical predictions using the proposed model. The sub-volume dimensions we first consider are 2 × 2.2 × 1.925 mm 3 , the material and model parameters are taken from the previous section 6.2.1.

The regular mesh used for DIC is used to built a 3D mesh of 8-node elements (see Fig. 21). By this way, we can apply the same segmentation procedure 6.2.3. 3D observations of cracks morphology by using numerical simulations and comparison with experiments

In this section, we will observe the crack propagation not only on the surface, but also within the volume. Moreover, in order to demonstrate the performance of the new procedures, an finer unstructured mesh is used. The have been identified by an inverse approach. Finally, the simulation and experiments are compared for the same load.

Figure 1 :

 1 Figure 1: Description of experimental procedures: (a) geometry of specimen and counter electrode (b) schematic illustration of experimental set-up

Figure 2 :

 2 Figure 2: Crack detection by using image processing based on DIC: (a) displacement field; (b) Postprocessed displacement fields; (c) crack paths obtained by segmentation

Fig

  Fig. 2(b)

Figure 3 :Figure 4 :

 34 Fig.3

  Simulation of slip dissolution induced SCC is an open question in the literature. The mechanical properties of the passive films forming on the specimen surface are difficult to obtain. Therefore the simulation of film rupture is not possible. In this work, we propose a new way to model our four principal points (P1,P2,P3 and P4) of SCC processes, in which corrosion occurs when the dynamic at crack tip (called activated crack) satisfies a condition related to strain rate (introduced in the next section). SCC provides here the removal of material, so the control volume just ahead of crack tip is softening. Consequently, crack propagation occurs at a stress intensity factor K ISCC below the fracture toughness K IC or at a stress σ SCC below the ultimate tensile stress σ C .

Figure 5 :

 5 Figure 5: Description of the influence of slip dissolution on fracture mechanics by using an equivalent model based on diffusion of degradation

Fig. 6 )Figure 6 :

 66 Fig. 6) subjected to SCC. The degradation coefficient is taken D(x) = D 0 in the surface, and D(x) is a monotonically increasing function of depth (x in

Figure 7 :

 7 Figure 7: Influence of diffusing degradation on the fracture toughness

  ) , γ(d, ∇d) is the crack density function per unit volume (see e.g. Miehe et al. (2010); Nguyen et al. (2015)), defined by:

Figure 8 :

 8 Figure 8: Crack extraction by using DIC: crack propagation for several loading steps

Figure 9 :

 9 Figure 9: Crack onset by using result of DIC for three different loading steps: (a) first loading step; (b) second loading step; (c) third loading step

  Compute diffusing degradation field c n+1 (x) 2.1 Update the crack tips 2.2 Check crack propagation rate if ḋ ≥ d c at the crack tips Solve the diffusion problem (31) end 3. Compute displacement field u n+1 (x) 3.1 Update the elastic parameters: E(x) = g(d(x), k)E 0 and ν(x) = g(d(x), k)ν 0 3.2 Solve the displacement problem (40) 4. (.) n ←-(.) n+1 and go to (1).

Figure 10 :

 10 Figure 10: Tensile test of a structure containing a semi elliptical crack: (a) geometry and boundary conditions; (b) surface view of FEM mesh; (c) mesh view of process zone

  )(b)(c) for 3D view and Fig. 12(a)(b)(c) for a plane of investigation.

Figure 11 :

 11 Figure 11: Tensile test of an homogeneous sample with a semi elliptical crack: the iso surface of phase field d(x) = 0.95 is plotted. Figures (a), (b) and (c) depict the cracks onset and propagation by using the classical phase field method and correspond to U = 0.0031 mm, U = 0.0032 mm, and U = 0.0033 mm, respectively. Figures. (d), (e) and (f) depict crack propagation for the model including both phases and SCC correspond to U = 0.002 mm, U = 0.0022 mm, and U = 0.0023 mm, respectively.

Figure 12 :

 12 Figure 12: Crack propagation in a chosen plane (z = 0.5H). The phase field d(x) is plotted. Figures (a), (b) and (c) depict the cracks onset and propagation by using the classical phase field method and correspond to U = 0.0031 mm, U = 0.0032 mm, and U = 0.0033 mm, respectively. Figures. (d), (e) and (f) depict crack propagation for the model including both phases and SCC correspond to U = 0.002 mm, U = 0.0022 mm, and U = 0.0023 mm, respectively.

Figure 13 :

 13 Figure 13: Observation of SCC in a chosen plane (z = 0.5). The diffusing degradation c(x) is plotted. Figures (a), (b) and (c) depict the degradation correspond to U = 0.00205 mm, U = 0.0021 mm, and U = 0.00215 mm, respectively

Figure 14 :

 14 Figure 14: The influence of SCC on fracture energy: distribution of g c along a line of investigation corresponding to U = 0.0021 mm. We observe here the variation in zone about of 0.15 mm above crack tip

Figure 15 :

 15 Figure 15: Evolution of surface energy along a line of investigation, corresponding to one undamaged step and three damaged steps.

Figure 16 :

 16 Figure 16: Comparison of stress-displacement curve for two models.

Figure 17 :Figure 18 :

 1718 Figure 17: The iso surface of phase field d(x) = 0.95 is plotted. Figures (a), (b) and (c) depict the cracks onset and propagation by using the classical phase field method and correspond to ε zz = 0.0022, ε zz = 0.0025, and ε zz = 0.0027, respectively

Figure 19 :

 19 Figure 19: Main principle to define cost function of inverse analysis algorithm

Fig. 20 )

 20 Fig.20). These spurious local displacement fluctuations could induce non physical damage. Hence, in order to avoid this issue, we propose to filter the DIC measurements along the sub-volume boundary at each experimental load step by a Savitzky-Golay filtering (see[START_REF] Orfanidis | Introduction to signal processing[END_REF]), as illustrated in Fig.20(black points).

Figure 20 :

 20 Figure 20: Displacement of boundary for an experimental load step obtained by DIC and its Savitzky-Golay filtering

  45) Where k = 1..N are the number of crack nodes obtained from experiments; d exp k = 1, and d num k

Figure 21 :

 21 Figure 21: Experimental validation, first sub-volume (2 × 2.2 × 1.925 mm 3 ): (a) geometry; (b) FEM mesh.

  considered sub-volume dimensions are 2.2 × 2.2 × 2 mm 3 (seeFigs. 26), the material and model parameters are taken from the previous section 6.2.1.The mesh of the sub-volume consists in 668972 tetrahedal elements with linear interpolation (4-node elements). The mesh is refined in the vicinity of the surface x = 0 in contact with the aggressive solution in the experiments, as shown in Fig.26(b). The typical size of an element in the crack propagation zone is about h min ≈ 0.0225 mm and about h max = 0.06 mm elsewhere in the domain. Note that, the length scale l is kept to the same value as for the previous tests (in the case of regular meshes). The evolving boundary conditions obtained from DIC data are applied on the sub-volume boundary during 100 steps. Before applying these boundary conditions both spatial and temporal interpolation and filtering as described in the previous section are used. Numerical and experimental crack propagations are compared for several loading steps. For the sake of clarity, in the 3D view ,we plot only the iso surface of phase field (d = 0.95). The results are presented in Fig.29for several loading steps. We observe the initiation of several cracks in the step presented in Fig.29(a). In the next step, these cracks are merging and their propagation is more pronounced at the surface. A complex crack network

Figure 26 :Figure 27 :

 2627 Figure 26: Experimental validation, sub-volume (2 × 2.2 × 2.2 mm 3 ): (a) geometry; (b) FEM mesh.

Figure 28 :

 28 Figure 28: Comparison between experimental crack obtained from DIC procedures and from phase field method of the Inconel 600 alloy sample in sub-volume size 2 × 2.2 × 2.2 mm 3 , for several loading steps (bottom: experiment; top: numerical simulation.)

Table 1 :

 1 Composition of the Inconel 600 alloy

on the simulation results to construct a binary map (cracked or undamaged) of elements on the sample surface. The mesh of the sub-volume consists in 50400 elements with linear interpolation. The evolving boundary conditions obtained from DIC data are applied on the boundary of the sub-volume during 95 steps. Before applying these boundary conditions both spatial and temporal interpolation and filtering as described in the previous section are used. Note that, the DIC displacement has been interpolated in time into 10 sub-increment to ensure the convergence of the staggered scheme. We compare the crack network obtained from the simulation with that one from the experiment at the same loading step. The result is depicted in Fig. 22, We perform the same analysis on other sub-volumes of different dimensions. The results are depicted in Figs. 23, 24 for two different sub-volume sizes. We obtain here a good agreement between the experiment and the

This error is bounded between 0 (if f exp and f num are identical) and 1

(when fluctuations of f exp and f num are opposite). The obtained result is plotted in Fig. 25. The sub-volume size decrease induces a reduction in the problem due to lack of micro-structure details. Hence, we obtain the better agreement, i.e the crack morphology error will be smaller.

Aknowledgements

The authors would like to acknowledge the National French Research Agency (ANR) for its financial support under contract MATETPRO ANR-12-RMNP-0020 (ECCOFIC project). The authors would like also to thank our partners AREVA, ANDRA, Institut de la Corrosion, ACXCOR and MISTRAS by their participation into the fruitful discussion during this work