A phase field method for modeling anodic dissolution induced stress corrosion crack propagation

Abstract The phase field method is a powerful tool for studying microstructural evolution in various domains of material sciences, including phase change, initiation and propagation of fracture. In this work, a new formulation is developed based on the phase field method for modeling stress corrosion cracking (SCC) induced by anodic dissolution. This method was applied for modeling SCC of an aluminum alloy (2xxx series) in a saline medium (NaCl), which allows considering the effects of both electrochemical and mechanical processes. The classical phase transition model for material dissolution is coupled with the mechanical problem in a robust manner, providing an efficient tool for studying the competition between electrochemical and mechanical contributions to fracture. A numerical implementation based on finite elements is elaborated. The numerical results are compared to experimental data obtained by in situ microtomography.


Introduction
Stress-corrosion cracking (SCC) is a common failure mechanism characterized by a combination of a mechanical load, a corrosive environment and a susceptible material.
The occurrence of SCC induces premature failure of the material, which is recognized as potentially dangerous. In that context, a detailed investigation regarding SCC role is necessary to accurately predict the strength and failure of materials and structures.
Email address: julien.rethore@ec-nantes.fr (J. Réthoré ) 1 Now at GEM, CNRS UMR 6183 CNRS / Ecole Centrale de Nantes / Université de Nantes Various mechanisms have been proposed in the literature to explain the synergistic stress-corrosion, such as adsorption model, film rupture model, pre-existing active path model, and embrittlement model,...(see e.g. [1] for a detailed review). Among them, the mechanisms based on the local anodic dissolution of fresh metal created by slip bands emergence at the crack tip [2] is often regarded as one of the most model. It has thus long been recognized as an indispensable step in the production of pure metal salts, energy production in batteries, etc. As a consequence, developing numerical model to study its effect in terms of key engineering parameters is highly required.
A first group of methods based on the moving of sharp interface have been proposed to this aim. In these methods, the interface is represented explicitly to apply the boundary conditions. Finite Element Method (FEM) by using re-meshing technique is widely used within this context [3]. However, the solution is mesh dependent. The remeshing step at each loading step reveals to be computationally expensive in practice and becomes extremely difficult for complex material and/or 3D problems, see e.g. [4] for a detailed review. As a consequence, the Arbitrary Lagrangian-Eulerian (ALE) method [3], has been proposed as an advanced FEM technique. ALE technique shows some improvements in comparison with classical FEM, such as its ability to simulate rapid dissolution processes as in the case of pitting corrosion. This method gives a stable solution when the pit growth is modelled in the work of [5], or [6]. Another alternative technique based on Finite Volume Method (FVM) is proposed in the work of [7] to study localized corrosion.
In this work, a voxel method is used to represent the anodic front. The ionic transport and complete electrostatic problems are solved in the FVM context. However, lack of corrosion progresses due to the changes in the electronic nature of boundaries (from anode to cathode and vice versa) is the drawback of this method. Some other contributions in the context of sharp interface methods are: the work of [8] to study SCC, in which the average current density (computed from the transport equation) is used as a criterion to model the crack growth. However the shape of the crack front must be predefined and it is difficult to simulate complex dissolution morphologies.
The level-set technique [9], or cellular automata [10], are also used in order to handle the complex shape of corrosion surfaces. Moreover, in these works, the electrochemical problem is often simplified.
Recently, the phase field method, where the sharp interface is described by a smeared interface, has been proposed [11,12,13,14,15,16]. This method becomes an efficient tool for the simulation of microstructural evolution in material science. With this method, the explicit tracking of interfaces is not required, i.e. all numerical difficulties encountered in the sharp-interface models are avoided. Therefore, the phase field method is capable to handle arbitrarily complex geometries, even in 3D. The phase field model has been used to investigate the effects of SCC in many works, such as in [17] to model the hydrogeninduced stress-corrosion cracking, or in [18] to account the effects of diffusion of chloride ions in hardened cement paste. The application of this model can also be found in [19] to simulate the corrosion kinetics under a dual-oxidant atmosphere, or to investigate the crack propagation in brittle materials induced by the variations of concentration (electrolyte-LiFePO 4 ) in the work [20] The application of the phase field method to simulate material dissolution can be found in the work of [21]. The authors have proposed a quantitative scheme, where the phase transformation is controlled by multicomponent diffusion. This model has been applied to study the dissolution in Ti-Al-V alloy. In the work of [22], the authors use a phase field model to investigate the dissolution of primary particles in aluminum phase under isothermal diffusion-controlled processes. The dynamics of liquid-solid interfaces resulting from precipitation and/or dissolution processing is studied by [23] by using the phase field method. In the work by [24], SCC induced crack propagation is studied, however, the electrochemical contribution is simplified and not investigated in details. More recently in the work by [25,26], this method is applied to simulate the pitting corrosion. However, a method to describe all aspects of electrochemical and mechanical processes inducing material dissolution or fracture growth is still missing. This work started from ideas presented in the work of [25,27] where the KKS (Kim-Kim-Suzuki) model [11] is adopted with new normalization procedures to describe the dissolution processes (phase transformation from solid to liquid). We propose herein a coupling with mechanical equilibrium to model dissolution induced SCC. It is worth mentioning that, differing from [26], the mechanical process is here directly accounted in the free energy functional, providing a performance scheme to investigate the competition between electrochemical and mechanical contribution. We also adopt a robust numerical framework by using the unilateral contact formulation and history functional following [16,28]. This ensures that the dissolution can be accelerated under tensile mechanical states, and that effect of material microstructure can be considered as mentioned in [24].
It is expected that the proposed model can simulate both the anisotropy in the elastic behaviour and in fracture evolution in the case of SCC.
The outline of the paper is as follows: first some fundamental basics of the model are introduced in Section 2 where several assumptions will be proposed. In Section 3, a new phase field framework able to model effects of both electrochemical and mechanical processes on crack growth is proposed. Then, the details of the numerical implementation are presented in Sections 4 and 5. Finally, the performances of the new model are illustrated through several numerical examples, which are compared to recent experimental observations in Section 6.

Fundamentals of the proposed method
Let Ω ⊂ R D be an open domain describing a corrosion system, which generally contains three domains: (i) metal or solid phase; (ii) electrolyte solution or liquid phase; (iii) an interfacial region where the corrosion occurs.
The characteristic of first metal solid domain is described by temporally and spatially constant metal atom concentration c alloy solid . That is identified based on its chemical composition. For the example of an alloy, it is composed from chemical element Me i with molar masses M i , regarding to be identical with molar masses of corresponding ions Me n+ i . Denoted f M,i [%], the mole fraction of species i, sometimes also referred to as "molar fraction" is computed by the following: where f m,i [%] is mass fraction of species i. Following the work by [29], the average molar mass is computed as: 4 and the concentration c alloy solid of metal atoms in the solid phase now reads: The corrosion progresses due to the dissolution of Me into aqueous solution as Me n+ ions at the anode boundary. Thus, the anodic reaction involves oxidation of Me as follows: That provides a diffusional transport of the metal ions from solid domain to liquid domain across the corrosion boundary. As noticed in [30,29], this corrosion boundary is After that, the evolution of corrosion is then controlled by the diffusion of metal ions away from the corrosion boundary (diffusion-controlled corrosion). In case of instantaneous dissolution of a metal alloy, dissolution takes place by orders of magnitude faster than the diffusion of dissolved ions, so that the metal ion concentration remains constantly at the saturation limit c sat , i.e. chemical equilibrium is always achieved at the electrode boundary.
According to the Rankine-Hugoniot condition, the rate of the metal dissolution depends on the ions concentration and its gradient along the interfacial region. This is usually described by the equilibrium between the dissolved metal atoms flux and the velocity of the moving interface [29]. In this work, we use the phase field method to model the dissolution of a metal which is corroding. The corrosion system is now described by a phase field variable φ. φ = 1 is solid phase, φ = 0 is liquid phase, and the corrosion boundary (sharp interface) is replaced by a thin interface within the range 0 < φ < 1, (it is often chosen that φ = 0.05 − 0.95 defines transition zone). The corrosion process 5 is then simulated by the evolution of auxiliary field φ over the entire system. The effect of the ions diffusion induced metal dissolution is modeled following the work by [25] who introduce the normalized concentration: The ions (or atoms) concentration now reads c s for metal solid and c l for electrolyte, which are defined in the normalized context by the following: where c alloy The assumption of KKS model [11] is used, in which the transition region is represented as a mixture of the two phases with different chemical compositions. Therefore, the ions concentration at any material point can be described by the following: where h(φ) can be understood as the interpolation function, satisfying h(φ = 0) = 0 to ensure that c = c l in liquid phase, h(φ = 1) = 1 and ∂h(φ) ∂φ | φ=0, φ=1 = 0, to ensure that c = c s in solid phase; one choice: h(φ) = −2φ 3 + 3φ 2 . This function also plays an important role to describe the loss of material strength discussed in the next section.
In the normalized framework, the metal solid and electrolyte ions concentration, de- where c se , c le are the saturated concentration of solid and liquid phase after being normalized. The details of the corrosion system in a normalized context is depicted in Fig. 2. In this work, the dissolution of material is not only due to electrochemical processes, but also to mechanical loading. Therefore, the phase field evolution φ = 1 → φ = 0 represents both interpretations: material dissolution and crack creation. In addition, different schemes based on different diffusion paths are used to describe the breakdown of the passive film, as depicted in Fig. 2, following the work in [31]. In most of situations, the volume of the electrolyte is much larger than the one occupied by the solid metal and

Phase field model
Under the framework of variational approaches, we introduce the following energy functional for a cracked body: The term E u (u, φ) represents the elastic energy stored in the cracked body, E c (c, φ) describes the electrochemical energy, and E φ (φ, ∇φ) is the energy required for material dissolution and/or crack creation. Then, the state variables are the displacement field u, the concentration field c and the phase field φ.
The proposed model is based on a small strain framework incorporating three fields: displacement, concentration and phase field. Hence the energy functional contains three 8 main contribution as described in Eq. 10. In this multiphysics process, there are three main interactions: • Diffusion with mechanical and material dissolution/fracture effects • Mechanics with diffusion and material dissolution/fracture effects • Dissolution/fracture with diffusion and mechanical effects However, in this work we try to capture the effect of material dissolution induced SCC by the simplest model. Hence, only the following most important factors are considered: (1) coupling between diffusion and dissolution/fracture (2) coupling between dissolution/fracture and mechanics. Note that, as the "bridging interactions" 2 , we can always capture the interactions among diffusion, mechanic and dissolution/fracture

Elastic energy
The elastic energy E u (u, φ) in the phase field framework is defined as the following: where ψ u ε(u), φ is the strain energy density of a damageable material. Small strain assumption is used, hence ε = 1 2 ∇u + ∇u T . To ensure that dissolution of the material is only accelerated by the tensile part of strain energy, and to prevent cracking in compression modes, we use here a unilateral contact formulation. Following the work of [32,33,28], the strain is decomposed into deviatoric ε dev and spheric ε sph parts. Then, it is assumed that phase evolution (dissolution or crack propagation) is sensitive to volumetric expansion and shear.
2 An example to illustrate this property, the diffusion process will affect the dissolution/fracture and the mechanical process is directly influenced by dissolution/fracture, i.e the diffusion will affect mechanical behavior via "bridging interactions" where C 0 denotes the initial elastic tensor of the material, possibly anisotropic. If k 0 denotes the corresponding bulk modulus for the undamaged material (relating the spherical part of the strain to the spherical part of the stress), hence the effective elastic tensor is: with the sign function sign − (x) = 1 if x < 0 and sign − (x) = 0 if x ≥ 0. The elastic density energy is now rewritten as: The interpolation function h(φ) can be understood as the degradation function. h(φ = 0) = 0 for fully material dissolution/damage describes the loss of material strength in this region and h ′ (φ = 0, φ = 1) = 0 ensures that the elastic energy converges to a finite value where the domain is locally intact/fully damaged.
Remark The effects of plastic strain can be easily accounted in this model by using e.g.
the idea in the work [34], where the authors have proposed a new formulation of a small-strain continuum theory for brittle fracture in elasticplastic solids. It could be done by using an additive decomposition of the strain ε into its elastic ε e and plastic ε p parts. In this case, contained plastic deformation can take place close to notch roots, and the new model will be a system of coupled equations governing the evolution of displacement, phase-field, concentration and accumulated plastic strain. However, for the sake of simplicity the effects of plastic deformation are not considered here. We keep it in the future work.

Electrochemical energy
Crack creation affected by the electrochemical process is described by the electrochem- where ψ c (c, φ) is the electrochemical energy density. Many models have been proposed to describe this effect, e.g. the WBM (Wheeler-Boettinger-McFadden) model [35,36], based on the assumption that any point within the interfacial region is assumed to be a mixture of solid and liquid both with the same composition; or under hypothesis that the transition zone is a mixture of solid and liquid with different compositions, but constant in their ratio [37]; or some others works, such as Khachaturyan's model [38] or Folch-Plapp's model [39].
In this contribution, the KKS model [11] is adopted. It postulates that the transition region is represented as a mixture of the two phases with different chemical compositions, but equal in chemical potentials: where ψ cs , ψ cl are local chemical energy densities of solid metal and solution, respectively.
Being evaluated in the work of [40], the expressions of ψ cs , and ψ cl do not impact the evolution of coexisting phases, hence the simplest form is used as an approximation: In these equations, c se and c le are the equilibrium concentrations of the coexisting phases defined in Eqs. 8,9. A is the free energy density curvature which is often identified from thermodynamical databases [40,27].
According to the KKS model [11], the condition for equilibrium at the transition zone is: From (7) and (18), the expressions of concentrations in solid and solution are written as the following: Finally, by using the Eqs. (16), (17) and (19), the electrochemical energy is properly described by two state variables, namely the phase field φ and the ions concentration c.

Surface energy
The surface energy is the amount of energy released upon the dissolution of metal surface or the creation of new fracture surface. It is described by the following expression: With ψ φ (φ, ∇φ) is surface energy density function, reads: where is the double well potential, ω φ is the height of the imposed double-well energy barrier and α φ is the gradient energy coefficient. In the present phase field corrosion model, α φ and ω φ are correlated to the interface energy σ φ and its thickness ℓ of the anodic dissolution model as: where α * ≈ 2.94 corresponds to the transition zone taken within the range 0.05 < φ < 0.95, see e.g. [41] for more details.
Note that, these two parameters ( ω φ and α φ ) are also directly related to the parameters of variational principle proposed to classical fracture mechanics [14,15]: where g c is fracture energy (amount of energy dissipated upon the creation of a unit fracture surface). Using this expression of the surface energy ensures that the Griffth criterion is fulfilled. ℓ denotes the size of transition zone or interface between solid and liquid phases. It 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 ℓ as small as possible to better approximate brittle fracture, with regards to the size of the mesh. In the second case, ℓ should be identified from experimental data. Such analysis is presently under study with experimental validations in [42,43].
Remark In the present model, the surface energy implies two different physical interpretations: (SE1) energy of the material dissolution provided by electrochemical processing or mass diffusion; (SE2) energy of the fracture creation provided by mechanical loading. In the first case (SE1), the surface energy parameters are understood as the classical parameters of phase transformation model [11], including (ω φ and α φ ). In the second case (SE2), the surface energy is now based on the mechanical fracture energy, defined as the product of Griffith energy g c with the crack surface density functional The definition of the surface energy allows to couple both phenomena but it could lead to some difficulties in its experimental identification. A better way is to divide into two kinds of parameter, one related to mechanical and the other one related to electrochemical processing. However, this scheme also need to define multi-phase-field (each phase-field corresponding to one process). This will change the physical meaning of damage (material dissolution/crack propagation). More investigations in future works need to be conducted on this point to make the model stronger.

Variational principle
The total energy is now rewritten as E = Ω ψ dΩ in which: is identified as total density energy.

Variational principle for mechanical problem
The displacement solution is found by solving the variational problem: where S u = {u|u(x) =ū on ∂Ω u , u ∈ H 1 (Ω)} and W ext = Ω f · u dΩ + ∂Ω F F · u dΓ with f and F body forces and prescribed traction over the boundary ∂Ω F . We obtain the classical weak form for u(x).
With the above described strain energy, the Cauchy stress defined as σ = ∂ψ u ∂ε is given by using (12) and (13): where M is the diffusion mobility, which can be evaluated based on an analogy between the Cahn−Hilliard equation (28) and the Fick's second law, being written as: where D is the diffusion coefficient. In this work, the diffusion coefficient D is assumed to be constant for both phases (i.e M is constant), so the diffusion problem to be solved to evaluate the concentration c(x, t) at time t is:

Strong form of phase problem
The phase field evolution is described by the Allen−Cahn equation: where C h (φ) = C 0 (x) − k 0 (x) 1 ⊗ 1 sign − (tr ε (x, τ )) and L φ is the interface kinetics parameter, that describes the transformation rate.
In addition, to handle loading and unloading histories, we introduce the history function following the work of [16,44] H( which will be instituted (31), yields the phase field problem to be solved to evaluate the phase field φ(x, t) at time t:

Weak form
Using the variation in the concentration for Eq. (30), and the variation in the phase field for Eq. (33), the corresponding weak form is obtained:

Numerical implementation
Finite Element Method is used for the space discretization whereas the θ-method is  (26) In the initial state, the passive film forms over metal solid, crack faces and also 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 when the plastic strain reaches a critical value; this process allows the dissolution of metal in this region and promotes the advance of cracks. The breakdown of the passive film can be simulated by using the interface kinetics parameter L φ as a function of crack tip strain rate (or in some cases creep rate) e.g. the approach proposed in the work [26]. However, the determination of strain rate is still an open question for both numerical model and experimental technique [45], which requires further efforts to make it applicable. Therefore, for sake of simplicity the breakdown of passive film is manually activated in this study and we do not consider the repassivation process.

Material parameters
The aluminum-copper alloys are known to be prone to both Intergranular Corrosion (IGC) and Intergranular Stress Corrosion Cracking (IGSCC). IGSCC initiates from localized corrosion sites then propagation is associated to the existence of preferential anodic path along the grain boundaries, thought to be caused by either solute depleted zones or anodic precipitates [46].
Due to its X-ray absorption properties allowing easy 3D synchrotron radiation imaging and its sensitivity towards dissolution induced SCC, an aluminum alloy AA 2024-T3 was chosen for confronting experimental and numerical results. It is composed mostly of aluminum, copper and magnesium, with a typical mass density of ρ alloy = 2.78 g/cm 3 .
The complete chemical composition is given in Table. 1, where the mole fraction is com-  chemical compositions, and taken as the value for pure Al, with D = 5.41 × 10 −10 m 2 /s [47].
As mentioned in Section 3.2, the free energy density curvature A could be identified by fitting the chemical free energies obtained from Eq. 17 with ones obtained from thermodynamic databases. Following [48,40], the chemical driving force for the phase transformation is ∆G = 6.9 × 10 8 J/m 3 . That gives A = 7.33 × 10 8 J/mol. (See Fig. 3).
Note that the interface energy can be either taken as the dissolution parameter related to electrochemical processing, or fracture parameter related to the mechanical loading.
In this work, we chose σ φ = 120 N/m. Cause of lack of experimental data, the size of transition zone is here considered as a pure numerical parameter of the regularized model and taken ℓ = 2.5 µm. The interface kinetics energy is chosen L φ = 0.15 l/(J.s).
For easier numerical implementation, all parameters will be normalized following the work by [40,27].

Case of pure local dissolution (semi-circular pitting growth)
In the first example, the phase field model is used to simulate the localized dissolution of metal in an homogeneous structure when the passive film is broken. This mode of corrosion is described in the literature as a semi-circular pit growth. An initial semicircular shaped pit is subjected to the corrosion process. Note that only electrochemical process is considered in this first example, i.e the mechanical loading is taken null and the evolution of phase field is assumed to be isotropic for the whole simulation. for the case of strong passive film, wherein it does not dissolve with metal, the pit keeps its semi-circular shape.This prediction is qualitatively in good agreement with experimental observation [49]. The concentration is recomputed in the real context (unnormalized), and then plotted along a line of investigation in (Fig. 6)  In the situation, where the passive is dissolved (totally for the case of weak passive film and partially for the case of Lacy cover film) along with metal, the obtained results are plotted in Fig. 7. The corrosion process is more promoted in the surface than in depth, a semi-elliptical damage form is observed. More interestingly, the Lacy cover form seems to protect the material better than weak passive film. Hence, at the same time step the case of Lacy cover passive film, the material is dissolved less than for the case of the weak passive film. These results reproduce very well th experimental observations in [49].

Cases of stress corrosion cracking
SCC tests were performed on AA 2021-T3 aluminum alloy in order to confront numerical simulation of damage growth and real propagation of stress corrosion damage, for two configurations for both numerical and experimental approaches: dissolution enhanced damage or mechanics enhanced damage.

Experimental procedure and results
The samples were machined by Electrical Discharge Machining from a 2 mm plate.
The loading of the sample was perpendicular to the rolling direction of the plate. The sample geometry is depicted in Fig. 8. In order to obtain high resolution 3D geometrical information on cracks, in situ SCC experiments are performed on Psiche beamline at Synchrotron SOLEIL. The experimental device consists in a miniaturized tensile machine on which samples can be surrounded by a saline media (NaCl 35 g/L at room temperature) contained in a Perspex cylinder which also allows electrochemical measurements and control (Fig. 9).
The tensile specimen is used as the Working Electrode (WE), a saturated calomel electrode (SCE) is used as a reference and a Platinum wire as the counter electrode.
Reference and counter electrodes are placed in the upper part of the cell so that they dont interact with the X-ray beam (Fig. 9). Metallic parts (grips) of the cell are protected by a varnish in order to avoid any galvanic coupling (Fig. 9). SCC initiation and propagation  pitting potentials corresponding to pitting localized at grain boundaries and within the grains respectively, the latter being higher than the former [50]. The pitting potential related to intragranular pitting can be identified as EpI in Fig. 10. On this figure, the identification of the pitting potential related to intergranular pitting (i.e. dissolution of the solute denuded zones along the grain boundaries [50]) is more difficult, but it can be assumed from literature results [50] that polarizing the specimen at potentials ranged between EpII and EpI will enhance intergranular pitting and SCC, whereas applying a potential equal to EpIII will prevent corrosion.
In order to enhance local dissolution, a first sample (S1) is polarized at -570 mV/SCE (EpI) from the beginning of the test. The tensile stress was applied progressively by steps of different durations from 110% to 132% of yield stress (σ yield = 290 ± 5 MPa), with a strain rate of 1.25 × 10 −4 s −1 between each step (Fig. 11). Scans were regularly performed during deformation from 0% to 132% of yield stress. As already performed on stainless steel [51], the displacement is stopped during each scan and the potential is decreased towards less anodic domain where neither the matrix, nor the copper-depleted zone can dissolve (-800 mV/SCE (EpIII)) ( Fig. 11). Images are obtained in pink beam mode   Figure 11: Loading and polarization procedures for the sample S1 (dissolution contribution enhanced)

25
For the second sample, in order to promote cracking rather dissolution (pitting and/or galvanic corrosion), the potential is fixed at -640 mV/SCE (EpII) and the stress is fixed at 125% of yield stress in the first part of the test. During this first stage, the dissolution processes should therefore be limited. No corrosion damage is actually observed. Then the potential is progressively increased up to -570 mV/SCE (EpI) and the tensile stress up to 132% of yield stress (Fig. 12). In this second stage, the first cracks can be observed as soon as the final conditions of potential and stress are reached. The procedure for scanning was the same than for sample S1. Results are exploited at 130% of yield stress and summarized in Fig. 13 for specimen S1, for which the dissolution is enhanced. The resulting "cracks" are shallow and branched.
For the second case (specimen S2), the results are depicted in Fig. 14. It is observed that the cracks are longer and tighter than for sample S1. 28 The rupture of passive film allowing electrochemical processes is manually activated from the beginning at the initial crack location. Two cases are investigated: • (C1) Electrochemical process prevails on mechanical loading (according to the sample (S1)) • (C2) Mechanical loading prevails on electrochemical process (according to the sam- In the first case (C1), the structure is subjected to mechanical loading (until critical state) faster than in the second case (C2). Moreover, the SCC is deactivated in the first period of time, according to the loading condition of sample (S2). The numerical simulation is performed in 2D model, so that it shows a shorter of SCC processes time when compared to the experiment (3D in reality). The details of boundary conditions are given in Table 2.  The material is supposed to be homogeneous, elastic and isotropic with E = 73.1 GPa and ν = 0.33. The result of damage evolution for the first case (C1) is depicted in Fig. 16. The damage evolution in the first period is quite similar as in the example of pure local dissolution, i.e. the semi-circular shape remains. During the second period, when the loading is increased, the fracture morphology is sharper as expected for mechanically driven failure.
In the second case (C2), when mechanical loading prevails on electrochemical processes, the damage evolution is depicted in Fig. 17. A classical crack shape (mechanical The corresponding von Mises stress is computed and plotted in Fig. 18 Figure 19 represents a detailed comparison of the damage geometry between two cases. The damage surfaces (crack edge) are extracted, showing the famous wavy-like shape as reported in the literature [52,53]. More interestingly, when the mechanical contribution is increased, the amplitude of this wave decreases, providing a shaper damage morphology.
These results are in agreement with the theoretical predictions in [52,54] The comparison with experimental observations in similar loading conditions is presented in Fig. 20 for the first case (C1) of the electrochemical contribution enhanced and  (a) damage morphology (C2) (b) damage morphology (C1)    for crack propagation, which is mainly governed by dissolution.

Conclusion
A phase field framework is introduced in this work to simulate the fracture growth induced by stress corrosion. The mechanical problem is coupled with a dissolution model in a robust manner. Moreover, the impact of material anisotropy has also been incorporated leading to directional effects, possibly along several preferential directions as described for a simple situation in Appendix A. The analyses of fracture morphology provided by various SCC conditions, qualitatively agree with experimental observations. That demonstrated the performance of the proposed framework to model all aspects of SCC.
Among the prospects of this work, one can mention the further extension to the heterogeneous materials or polycrystalline systems by coupling with a cohesive model for interfaces or grain boundaries. Also, the discussion concerning the parameters of the surface energy is difficult because they could be related either to the electrochemical parameters of the dissolution model, or to the mechanical failure process. This point is still an open question and needs to be investigated for the further development of the model. Quantitative validation using experimental results is also an essential step for further developments.

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 are also grateful for the beam time awarded by synchrotron SOLEIL on the Psiche beamline (20140951 and 20150856 accepted proposals) and acknowledge warmly the technical and scientific support they benefited. The authors would like also to thank their partners Institut de la Corrosion, AREVA NP, Andra, MISTRAS Group and Pierre Combrade for their participation in the fruitful discussions during this work.

Appendix A. Extending to the anisotropic behavior
The proposed model can be easily extended to anisotropic materials. In that case, an anisotropic fracture surface density function is written by the following expression: where ω is a second-order structural tensor being invariant with respect to rotations (characterizing the material anisotropy).
As considered in [33], to make the energy release rate orientation-dependent, the tensor ω can be defined by: where N denotes the unit vector normal to the preferential cleavage plane (with respect to the material coordinates), and β ≫ 0 is used to penalize the damage on planes not normal to N. Hence, in the case of isotropic material β = 0. Note that, for systems with many preferential directions, one phase field variable is defined for each direction. We suggest to follow the work in [28,55] for more details.
The elastic moduli mentioned in (12) possibly depend on the material orientation. An example of the case of polycrystalline material, the rotation of the grains is taken into account to determine the elastic stiffness tensor for each grain. This can be achieved by using the reference transformation tensor and grains coordinate system. Suppose that the components of the elastic stiffness tensor matrix C 0 g (in Voigt notation) are given in grains coordinates, then the position-dependent elastic stiffness tensor with respect to the reference coordinate system is given by: where P is the transformation tensor in Voigt's notation.
The numerical example presented in section 6.3.2 is considered to study the influence of material orientation. The anisotropic fracture evolution described in Eqs.