T-T Nguyen

J Bolivar

Y Shi

J Réthoré

A King

M Fregonese

J Adrien

J-Y Buffiere

M-C Baietto

A Phase field method for modelling anodic dissolution induced stress corrosion crack propagation

Keywords: Stress corrosion cracking, anodic dissolution, crack propagation, phase field

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 modelling 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.

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. [START_REF] Burleigh | The postulated mechanisms for stress corrosion cracking of aluminum alloys: A review of the literature 1980-1989[END_REF] 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 [START_REF] Magnin | Mechanism of brittle fracture in a ductile 316 alloy during stress corrosion[END_REF] 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 [START_REF] Sarkar | A numerical framework for the modeling of corrosive dissolution[END_REF]. 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. [START_REF] Duddu | Numerical of corrosion pit propagation using the combined extended finite element and level set method[END_REF] for a detailed review. As a consequence, the Arbitrary Lagrangian-Eulerian (ALE) method [START_REF] Sarkar | A numerical framework for the modeling of corrosive dissolution[END_REF], 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 [START_REF] Kota | Microstructure-based numerical modeling of pitting corrosion in 316 stainless steel[END_REF], or [START_REF] Sun | An arbitrary lagrangian-eulerian model for modelling the time-dependent evolution of crevice corrosion[END_REF]. Another alternative technique based on Finite Volume Method (FVM) is proposed in the work of [START_REF] Onishi | Numerical method for time-dependent localized corrosion analysis with moving boundaries by combining the finite volume method and voxel method[END_REF] 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 [START_REF] Gavrilov | Finite element calculation of crack propagation in type 304 stainless steel in diluted sulphuric acid solutions[END_REF] 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 [START_REF] Wojtan | Animating corrosion and erosion[END_REF], or cellular automata [START_REF] Pidaparti | Cellular automation approach to model aircraft corrosion pit damage growth[END_REF], 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 [START_REF] Kim | Phase-field model for binary alloys[END_REF][START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF][START_REF] Boettinger | Phase-field simulation of solidification[END_REF][START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Miehe | A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. 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 [START_REF] Falkenberg | Simulation of environmentally-assisted material degradation by a thermodynamically consistent phase-field model[END_REF] to model the hydrogeninduced stress-corrosion cracking, or in [START_REF] Wu | A phase-field approach to fracture coupled with diffusion[END_REF] to account the effects of diffusion of chloride ions in hardened cement paste. The application of this model can also be found in [START_REF] Wen | Phase-field modeling of corrosion kinetics under dual-oxidants[END_REF] 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 [START_REF] Liang | Phase-field modeling of diffusion-induced crack propagations in electrochemical systems[END_REF] The application of the phase field method to simulate material dissolution can be found in the work of [START_REF] Chen | Quantitative phase field modeling of diffusioncontrolled precipitate growth and dissolution in Ti-Al-V[END_REF]. 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 [START_REF] Kovačević | Solution of a phase-field model for dissolution of primary particles in binary aluminum alloys by an r-adaptive mesh-free method[END_REF], 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 [START_REF] Xu | Phase-field modeling of solute precipitation and dissolution[END_REF] by using the phase field method. In the work by [START_REF] Nguyen | A phase field method for modeling stress corrosion crack propagation in a nickel base alloy[END_REF], SCC induced crack propagation is studied, however, the electrochemical contribution is simplified and not investigated in details. More recently in the work by [START_REF] Mai | A phase field model for simulating the pitting corrosion[END_REF][START_REF] Mai | A phase field model for simulating the stress corrosion cracking initiated from pits[END_REF], 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 [START_REF] Mai | A phase field model for simulating the pitting corrosion[END_REF][START_REF] Abubakar | Phase field modeling of V2O5 hot corrosion kinetics in thermal barrier coatings[END_REF] where the KKS (Kim-Kim-Suzuki) model [START_REF] Kim | Phase-field model for binary alloys[END_REF] 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 [START_REF] Mai | A phase field model for simulating the stress corrosion cracking initiated from pits[END_REF], 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 [START_REF] Miehe | A phasefield model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF]. This ensures that the dissolution can be accelerated under tensile mechanical states, and that effect of material microstructure can be considered as mentioned in [START_REF] Nguyen | A phase field method for modeling stress corrosion crack propagation in a nickel base alloy[END_REF].

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:

f M,i = f m,i /M i j f m,i /M i , (1)
where f m,i [%] is mass fraction of species i. Following the work by [START_REF] Scheiner | Stable pitting corrosion of stainless steel as diffusioncontrolled dissolution process with a sharp moving electrode boundary[END_REF], the average molar mass is computed as:

M = i f M,i M i , (2)
and the concentration c alloy solid of metal atoms in the solid phase now reads:

c alloy solid = ρ alloy M . (3)
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:

Me → Me n+ + ne -. (4)
That provides a diffusional transport of the metal ions from solid domain to liquid domain across the corrosion boundary. As noticed in [START_REF] Tester | Diffusional effects in simulated localized corrosion[END_REF][START_REF] Scheiner | Stable pitting corrosion of stainless steel as diffusioncontrolled dissolution process with a sharp moving electrode boundary[END_REF], this corrosion boundary is represented with temporally and spatially constant saturation concentration of dissolved metal ions c sat . This value describes the maximum value of the overall ionic concentration that can be reached c ions = i c ions i . Actually, due to the diffusional transport of the metal ions, the overall ionic concentration in solution increases until a saturation exists with respect to the chloride or the hydroxychloride complex of the cations in the case of a saline solution, i.e. until chemical equilibrium between the ionic solution and the solid metal is achieved. This leads to the formation of a salt film with saturation concentration c sat .

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 [START_REF] Scheiner | Stable pitting corrosion of stainless steel as diffusioncontrolled dissolution process with a sharp moving electrode boundary[END_REF]. 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 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 [START_REF] Mai | A phase field model for simulating the pitting corrosion[END_REF] who introduce the normalized concentration:

c = c ions c alloy solid (5)
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: The assumption of KKS model [START_REF] Kim | Phase-field model for binary alloys[END_REF] 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:

c s = c alloy
c = h(φ)c s + 1 -h(φ) c l , (7)
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.

Metal

Electrolyte

Transition zone

Solid phase Electrolyte

Line of investigation Interface (transition zone) Film rupture 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 [START_REF] Chen | The influence of passive film damage on pitting corrosion[END_REF]. In most of situations, the volume of the electrolyte is much larger than the one occupied by the solid metal and where the initiation of SCC occurs. Therefore the ions concentration and phase field in

Phase field model

Under the framework of variational approaches, we introduce the following energy functional for a cracked body:

E(u, φ, c) = E u (u, φ) + E c (c, φ) + E φ (φ, ∇φ). ( 10
)
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 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:

E u (u, φ) = Ω ψ u ε(u), φ dΩ, (11)
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 [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF][START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF], 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.

ψ u ε(u), φ =            1 2 h(φ) ε : C 0 : ε if tr ε ≥ 0 1 2 ε sph : C 0 : ε sph + 1 2 h(φ) ε dev : C 0 : ε dev if tr ε < 0, ( 12
)
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:

C(φ) = h(φ)C 0 + k 0 1 ⊗ 1 1 -h(φ) sign -(tr ε), (13)
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:

ψ u ε(u), φ = 1 2 ε : C(φ) : ε (14)
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 [START_REF] Duda | A phase-field/gradient damage model for brittle fracture in elastic-plastic solids[END_REF], 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-

ical energy E c (c, φ). E c (c, φ) = Ω ψ c (c, φ) dΩ, (15)
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 [START_REF] Wheeler | Phase-field model for isothermal phase transitions in binary alloys[END_REF][START_REF] Wheeler | Phase-field model of solute trapping during solidification[END_REF],

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 [START_REF] Tiaden | The multiphase-field model with an integrated concept for modelling solute diffusion[END_REF]; or some others works, such as Khachaturyan's model [START_REF] Wang | Shape instability during precipitate growth in coherent solids[END_REF] or Folch-Plapp's model [START_REF] Folch | Quantitative phase-field modeling of two-phase growth[END_REF].

In this contribution, the KKS model [START_REF] Kim | Phase-field model for binary alloys[END_REF] 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:

ψ c (c, φ) = h(φ)ψ cs (c s ) + 1 -h(φ) ψ cl (c l ), (16)
where ψ cs , ψ cl are local chemical energy densities of solid metal and solution, respectively.

Being evaluated in the work of [START_REF] Hu | Thermodynamic description and growth kinetics of stoichiometric precipitates in the phase-field approach[END_REF], the expressions of ψ cs , and ψ cl do not impact the evolution of coexisting phases, hence the simplest form is used as an approximation:

ψ cs (c s ) = A(c s -c se ) 2 ψ cl (c l ) = A(c l -c le ) 2 . (17)
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 [START_REF] Hu | Thermodynamic description and growth kinetics of stoichiometric precipitates in the phase-field approach[END_REF][START_REF] Abubakar | Phase field modeling of V2O5 hot corrosion kinetics in thermal barrier coatings[END_REF].

According to the KKS model [START_REF] Kim | Phase-field model for binary alloys[END_REF], the condition for equilibrium at the transition zone is:

∂ψ cs (c s ) ∂c s = ∂ψ cl (c l ) ∂c l . (18)
From ( 7) and ( 18), the expressions of concentrations in solid and solution are written as the following:

c s = c + h(φ) -1 c le -c se c l = c + h(φ) c le -c se . (19)
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:

E φ (φ, ∇φ) = ω ψ φ (φ, ∇φ) dΩ. ( 20
)
With ψ φ (φ, ∇φ) is surface energy density function, reads:

ψ φ (φ, ∇φ) = ω φ g(φ) + α φ (∇φ) 2 , ( 21
)
where

g(φ) = φ 2 (1 -φ 2
) 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:

σ φ = √ α φ ω φ 3 √ 2 and ℓ = α * 2α φ ω φ , (22)
where α * ≈ 2.94 corresponds to the transition zone taken within the range 0.05 < φ < 0.95, see e.g. [START_REF] Ammar | Modelling and simulation of phase transformation-mechanics coupling using a phase field method[END_REF] for more details.

Note that, these two parameters ( ω φ and α φ ) are also directly related to the parameters of variational principle proposed to classical fracture mechanics [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF]:

ω φ = g c 2ℓ α φ = g c ℓ 2
, and then

ψ φ (φ, ∇φ) = g c g(φ) 2ℓ + ℓ 2 (∇φ) 2 , (23)
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 [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Nguyen | Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF].

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 [START_REF] Kim | Phase-field model for binary alloys[END_REF], 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

g(φ) 2ℓ + ℓ 2 (∇φ) 2 .
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:

ψ = ψ u + ψ c + ψ φ , (24)
is identified as total density energy.

Variational principle for mechanical problem

The displacement solution is found by solving the variational problem:

u(x) = Arg inf u∈Su E(u, φ, c) -W ext (25)
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).

Ω σ : ε(δu) dΩ = Ω f • δu dΩ + Ω F F • δu dΓ. ( 26
)
With the above described strain energy, the Cauchy stress defined as σ = ∂ψ u ∂ε is given by using ( 12) and ( 13):

σ = h(φ)C 0 + k 0 1 ⊗ 1 1 -h(φ) sign -(tr ε(u)) : ε(u) (27
∂c ∂t (x, t) = ∇ • M ∇ ∂ψ ∂c , ( 28
)
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:

M = D 2A , ( 29
)
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:

∂c ∂t (x, t) = M ∆ ∂ψ c ∂c (30)

Strong form of phase problem

The phase field evolution is described by the Allen-Cahn equation:

∂φ ∂t (x, t) = -L φ ∂ψ ∂φ -∇ • ∂ψ ∂∇φ = -L φ ∂ψ c ∂φ + 1 2 h ′ (φ) ε : C h (φ) : ε + ω φ g ′ (φ) -α φ △φ , (31)
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 [START_REF] Miehe | A phasefield model for rate-independent 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] H

(x, t) = max τ ∈[0,t] ε (x, τ ) : C 0 (x) -k 0 (x) 1 ⊗ 1 sign -(tr ε (x, τ )) : ε (x, τ ) , (32)
which will be instituted 1 2 ε : C h (φ) : ε in (31), yields the phase field problem to be solved to evaluate the phase field φ(x, t) at time t:

∂φ ∂t (x, t) = -L φ ∂E c ∂φ + h ′ (φ)H(x, t) + ω φ g ′ (φ) -α φ △φ (33)

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:

                   Ω ∂c ∂t δc dΩ + Ω M ∇ ∂ψ c ∂c • ∇δc dΩ = 0 Ω ∂φ ∂t δφ dΩ + Ω L φ ∂ψ c ∂φ + h ′ (φ)H(x, t) + ω φ g ′ (φ) δφ dΩ - Ω L φ α φ ∇φ∇δφ dΩ = 0 (34)

Numerical implementation

Finite Element Method is used for the space discretization whereas the θ-method is adopted for time integration. In order to reduce the computational cost, staggered scheme is used and two separate problems are solved: (1) electrochemical and phase field problems solved using a classical Newton-Raphson method; (2) linear equations of the displacement problem. The summary of overall algorithm is described in the following:

At time t 0 : initialize the displacement field u 0 (x), the concentration field c 0 (x), the phase field φ 0 (x), and the strain-history functional H 0 = 0.

At time t n+1 : given u n (x), c n (x) and φ n (x) (a) Compute the strain history functional H n , see Eq. ( 32) (b) Compute c n+1 (x) and φ n+1 (x) by solving [START_REF] Duda | A phase-field/gradient damage model for brittle fracture in elastic-plastic solids[END_REF] (c) Compute u n+1 (x) by solving [START_REF] Mai | A phase field model for simulating the stress corrosion cracking initiated from pits[END_REF] 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 [START_REF] Mai | A phase field model for simulating the stress corrosion cracking initiated from pits[END_REF]. However, the determination of strain rate is still an open question for both numerical model and experimental technique [START_REF] Shoji | Formulating stress corrosion cracking growth rates by combination of crack tip mechanics and crack tip oxidation kinetics[END_REF], 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.

Results and discussion

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 [START_REF] Little | An electrochemical framework to explain the intergranular stress corrosion behavior in two Al-Cu-Mg-Ag alloys as a function of aging[END_REF].

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 computed following section 2. Using Eq. 2, we obtain M = 27.747 g/mol. The metal atom concentration is then determined following Eq. 3, giving c alloy solid = 100.2 mol/l. The saturated concentration is taken as in the work by [START_REF] Guseva | Modelling the early stage time dependence of localised corrosion in aluminium alloys[END_REF]: c sat = 3 mol/l. Hence, the saturated concentration of solid and liquid phase are computed as in Eqs. ( 8), [START_REF] Wojtan | Animating corrosion and erosion[END_REF]: 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 [START_REF] Buhler | A thermodynamic assessment of the Al-Cu-Mg ternary system[END_REF][START_REF] Hu | Thermodynamic description and growth kinetics of stoichiometric precipitates in the phase-field approach[END_REF], 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).

c se = 1, c le = 3/100.2

. The diffusion coefficient is assumed to be homogeneous for all

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 [START_REF] Hu | Thermodynamic description and growth kinetics of stoichiometric precipitates in the phase-field approach[END_REF][START_REF] Abubakar | Phase field modeling of V2O5 hot corrosion kinetics in thermal barrier coatings[END_REF].

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.

A rectangular domain with size L × B = 200 × 140 µm 2 is analyzed, where the film rupture region is initiated as a semi-circular area with a diameter a = 10 µm on the upper edge y = B, as depicted in Fig. 4. In order to ensure the stability of numerical solution, the mesh size must be chosen fine enough, to have about 5-6 elements within the transition zone. Here, we take the element size h e = ℓ/10 = 0.5 µm. The material parameters are defined from the previous section 6. Figs. 5 depicts the evolution of the concentration (phase composition) and the phase field at several time steps for the first case of strong passive film. The result shows that, 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 [START_REF] Ernst | Pit growth studies in stainless steel foils. i. introduction and pit growth kinetics[END_REF]. The concentration is recomputed in the real context (unnormalized), and then plotted along a line of investigation in (Fig. 6), The obtained results demonstrated that the ions's concentration in liquid phase is always bellow the saturation value c sat = 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 [START_REF] Ernst | Pit growth studies in stainless steel foils. i. introduction and pit growth kinetics[END_REF].

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 Fig. 10 presents a polarization curve of L-T plane of the sample in the testing solution (3.5% wt. NaCl). For AA 2024 alloy, SCC is associated to the existence of two pitting potentials corresponding to pitting localized at grain boundaries and within the grains respectively, the latter being higher than the former [START_REF] Urushino | Stress-corrosion cracking of aged al-cu-mg alloys in nacl solution[END_REF]. 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 [START_REF] Urushino | Stress-corrosion cracking of aged al-cu-mg alloys in nacl solution[END_REF]) is more difficult, but it can be assumed from literature results [START_REF] Urushino | Stress-corrosion cracking of aged al-cu-mg alloys in nacl solution[END_REF] 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 [START_REF] King | Observations of intergranular stress corrosion cracking in a grain-mapped polycrystal[END_REF], 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 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. The structure will be numerically assigned to both electrochemical degradation and mechanical loading according to the previous experiments. The mechanical loading (ydisplacements) is applied along the upper edge (y = B). On the other hand, along the lower edge (y = 0), the y-displacements are fixed. Along these edges, x-displacements are free. The investigated structure is assumed to be covered by a strong passive film.

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 sample (S2))

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 Indeed, the crack begins to propagate only from the time when the chemical process is activated and it can be noticed that the liquid domain is growing around the initial notch what is not obtained for a purely mechanical loading.

The corresponding von Mises stress is computed and plotted in Fig. 18 at a chosen loading step. In the first situation (C1), with more shallow crack, the maximal von Mises stress at crack tip is about 190 MPa. This value is much less than yield stress, i.e the damage is here mainly due to electrochemical process. For the second case (C2), with shaper damage morphology, the maximal von Mises stress is observed about 320 MPa at crack tip, which is comparable with yield stress. This result confirms the shape of the 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 [START_REF] Asaro | Interface morphology development during stress corrosion cracking: Part i. via surface diffusion[END_REF][START_REF] Ståhle | Phase field modelling of stress corrosion[END_REF]. 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 [START_REF] Asaro | Interface morphology development during stress corrosion cracking: Part i. via surface diffusion[END_REF][START_REF] Grinfel'd | Instability of the equilibrium of a nonhydrostatically stressed body and a melt[END_REF] 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 Fig. 21 for the second situation (C2) when mechanical loading prevails on electrochemical process. Can be noted that the numerical simulations reproduce qualitatively the damage morphology obtained in the experiments for both situations. In our knowledge, that is ). Yet, once the suitable plastic deformation is reached (132%YS in this case), the load is maintained constant and cracks propagate with no increase of plasticity. On the other hand, the present model gives very good results for crack propagation simulation. It can therefore be concluded that if plasticity is involved in crack initiation (which is not accounted for by the model here), it is not a first order parameter 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.

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: C 0 = P T C 0 g P, (A. [START_REF] Sarkar | A numerical framework for the modeling of corrosive dissolution[END_REF] where P is the transformation tensor in Voigt's notation.

The numerical example presented in section 6.

  concentration values of solid metal and solution, respectively. The details of this description are depicted in Fig.1.

noted c alloy sat and c solution sat ,

 sat

Figure 1 :

 1 Figure 1: Schematic of a typical corrosion: (top) 2D representation; (bottom) illustration in a line of investigation

  the electrolyte region are assumed to be zero, i.e the film-rupture process is simplified in the numerical model by imposing the Dirichlet boundary conditions: c = 0, φ = 0 over this region. Otherwise, in the regions where the metal is still covered by passive film, the zero-flux ∂c ∂n = 0, and ∂φ ∂n = 0 are imposed.

Figure 2 :

 2 Figure 2: Schematic of a typical corrosion for different kinds of the passive film using normalized variables

) 4 . 2 .

 42 Variational principle for electrochemical and phase field problem 4.2.1. Strong form of electrochemical problem Cahn-Hilliard equation is used to describe the diffusion of ions from solid volume to electrolyte:

Figure 3 :

 3 Figure 3: Variation of composition with chemical driving force fitting with thermodynamic databases.

1 .

 1 Three kinds of passive film are investigated: (a) strong passive film, where the passive film is not dissolved together with metal; (b) the passive film is dissolved together with the metal dissolution, denoted weak passive film; (c) only half of the passive film fails and form the lacy cover. For time t > 0, the breakdown of passive film is activated by imposing the Dirichlet boundary conditions as follows: the film rupture region is defined by c = 0 and φ = 0. Neumann conditions ∂φ ∂n = 0, and ∂c ∂n = 0 apply along the rest

Figure 4 :

 4 Figure 4: Geometry and boundary conditions of Benchmark test, different schemes used to model the break down of passive film is illustrated

  (a) t = 60 min (b) t = 100 min (c ) t = 140 min

Figure 5 :

 5 Figure 5: Semi-circular pit growth simulation, evolution of concentration and phase field for the first case of strong passive film (size of x, y-axes in µm)

Figure 6 :

 6 Figure 6: Evolution of real concentration (phase composition) along the line of investigation

  (a) weak film t = 60 min (b) Lacy film t = 60 min (c) weak film t = 120 min (d) Lacy film t = 120 min (e) weak film t = 160 min (f) Lacy film t = 160 min

Figure 7 :

 7 Figure 7: Corrosion morphology for the case of weak passive film (left) and Lacy cover form of passive film (right) for several time steps. The semi-elliptical damage form is observed.

Figure 8 :

 8 Figure 8: AA 2024 tested samples: (a) general view of sample geometry; (b) loading description

Figure 9 :

 9 Figure 9: Experimental device: (a) general view; (b) disposition of the electrodes

Figure 10 :

 10 Figure 10: Polarization curve of the AA 2024 in a solution of 3.5 % wt. NaCl (plane L-T)

(

  E=25 keV) on a Hamamatsu Detector (2048 × 2048 pixels). The scan duration is about 5 minutes (1020 images, average of 10 radiographs for each, exposure time of 25 ms) with a 2.39 µm voxel size. Initiations of corrosion defects are observed by tomography from the beginning of the test.

.

Figure 11 :

 11 Figure 11: Loading and polarization procedures for the sample S1 (dissolution contribution enhanced)

Figure 12 :

 12 Figure 12: Loading and polarization procedures for the sample S2 (mechanical contribution enhanced)

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Morphology of the corrosion damage in the T-L and T-S directions for specimen S1 (scan performed at 130% of Yield Stress after 7h of test)

  . A classical crack shape (mechanical

Figure 16 :

 16 Figure 16: Case (C1) electrochemical processing prevails on mechanical loading, the damage evolution for time: t a = 10 min, t b = 50 min and t c = 115 min

Figure 17 :

 17 Figure 17: Case (C2) mechanical loading prevails on electrochemical processing, fracture evolution for time: t a = 20s min, t b = 60 min and t c = 120 min

  (a) damage morphology -case (C1) (b) corresponding von Mises stress filed (C1) (c) damage morphology -case (C2) (d) corresponding von Mises stress filed (C2)

Figure 18 :

 18 Figure 18: Comparison of damage morphologies and corresponding stress states between two situations (C1) electrochemical contribution enhanced at time t = 50 min and (C2) mechanical contribution enhanced at time t = 80 min

Figure 19 :

 19 Figure 19: Effects of mechanical and electrochemical contributions on damage phenomena. Detailed comparison of damage morphology between two situations.

Figure 20 :

 20 Figure 20: Case (C1) electrochemical processing prevails on mechanical loading, comparison of crack morphology of: (a) experimental observation and (b) predicted simulation

Figure 21 :

 21 Figure 21: Case (C2) mechanical loading prevails on electrochemical processing, comparison of crack morphology of: (a) experimental observation and (b) predicted simulation

3 . 2

 32 Figs.A.[START_REF] Kovačević | Solution of a phase-field model for dissolution of primary particles in binary aluminum alloys by an r-adaptive mesh-free method[END_REF] shows the crack propagation for both cases. In the first situation, we observe the material dissolution following the preferential plane. In the second case, the orientation dependent failure energy leads to a zigzag crack profile, which is already known in the literature, for example the work[START_REF] Wu | Zig-zag crack propagation in MgAl2O4 crystals[END_REF] on MgAl 2 O 4 crystals and the work[START_REF] Kitagawa | Crack-morphological aspects in fracture mechanics[END_REF] on the mild, stainless steel and copper. This model is therefore promising for modelling inter and transgranular SCC at smaller scale, when the crystallographic orientation and grain boundaries have to be accounted for.

Figure A. 22 :

 22 Figure A.22: Fracture morphology with anisotropic effect, the phase field is plotted corresponding to: (a) electrochemical processing prevails on mechanical loading; (b) mechanical loading prevails on electrochemical processing

Table 1 :

 1 Chemical composition of the AA 2024-T3 alloy

	Alloy	Element Mass fraction, Molar mass, Mole fraction
			f m,i [%]	M[g/mol]	f M,i [%]
		Al	93.21	26.58	95.860
		Cu	4.36	63.55	1.904
		Mg	1.49	24.31	1.701
		Mn	0.46	54.94	0.232
	AA 2024-T3	Fe	0.25	55.85	0.124
		Si	0.14	28.09	0.138
		Zn	0.07	65.38	0.030
		Ti	0.01	47.87	0.006
		Cr	<0.01	<52.00	<0.005

chemical compositions, and taken as the value for pure Al, with D = 5.41 × 10 -10 m 2 /s

[START_REF] Guseva | Modelling the early stage time dependence of localised corrosion in aluminium alloys[END_REF]

.

Table 2 :

 2 Loading details of two situations: (C1) and (C2), where V is the applied displacement rate for mechanical loading, SCC-D is diffusion coefficient controlling the electrochemical process

		Period [0 : 40 min]	Period [40 : 60 min]	Period [60 : 120 min]
		V (µm/s) SCC -D(m 2 /s) V (µm/s) SCC -D(m 2 /s) V (µm/s) SCC -D(m 2 /s)
	(C1)	0.02	5.41.10 -10	0.03	5.41.10 -10	0.03	5.41.10 -10
	(C2)	0.05	5.10 -11	0.03	Deactivated	0.03	5.41.10 -10

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 [START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF], 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 [START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF][START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF] for more details.

The elastic moduli mentioned in [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF] 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


Yewan Wang
email: yewan1.wang@orange.com

David Nörtershäuser
email: david.nortershauser@orange.com

Stéphane Le Masson
email: stephane.lemasson@orange.com

Jean-Marc Menaud
email: jean-marc.menaud@imt-atlantique.fr

Experimental Characterization of Variation in Power Consumption for Processors of Different generations

Keywords: Manufacturing variability, CPU Temperature, Power Estimation, Thermal Interface Material, Leakage current

Data centers are energy-hungry facilities. Building energy consumption predictive models for servers is one of the solutions to use efficiently the resources. However, physical experiments have shown that even under the same conditions, identical processors consume different amount of energy to complete the same task. While this manufacturing variability has been observed and studied before, there is lack of evidence supporting the hypotheses due to limited sampling data, especially from the thermal characteristics. In this article, we compare the power consumption among identical processors for two Intel processors series with the same TDP (Thermal Design Power) but from different generations. The observed power variation of the processors in newer generation is much greater than the older one. Then, we propose our hypotheses for the underlying causes and validate them under precisely controlled environmental conditions. The experimental results show that, with the increase of transistor densities, difference of thermal characteristics becomes larger among processors, which has non-negligible contribution to the variation of power consumption for modern processors. This observation reminds us of re-calibrating the precision of the current energy predictive models. The manufacturing variability has to be considered when building energy predictive models for homogeneous clusters.

I. INTRODUCTION

Today, data centers are estimated to consume 200 terawatt hours (TWh) each year, which is more than the need of some countries, such as Iran [START_REF] Jones | How to stop data centres from gobbling up the world's electricity[END_REF]. Building energy efficient data centers and reducing environmental impact become a great concern. Therefore, it is essential to have a preview of the energy consumption before planning the construction, for both economic and environmental benefits. Energy predictive power model is one of the approaches [START_REF] Vor Dem Berge | Modeling data center building blocks for energy-efficiency and thermal simulations[END_REF]. We are currently working on an energy predictive model for servers in using corresponding indicators, such as system usage, hardware configurations, power management technologies applied, as well as environmental information, in order to optimize the energy usage for servers in data centers and to avoid wasting energy on un-used part of the servers.

Much work has been done on building accurate models. Some of them propose high-accuracy software-level solutions [START_REF] Costa | Methodology of measurement for energy consumption of applications[END_REF] [4] [START_REF] Canuto | A methodology for full-system power modeling in heterogeneous data centers[END_REF]. The question is, if we build the power model for one server and validate the model precision based on the model outputs and the real measurements, can we trust the model and rely on the same precision to other servers in the homogeneous cluster or even in the whole data centers? In the experimentation of [START_REF] Wang | Potential effects on server power metering and modeling[END_REF], a same test suite composed of several CPU intensive workloads is executed one by one to 12 identical servers in a rack, 7.8% power variation is observed among them. Thermal effect is turned out to be one of the causes: the rise of the ambient temperature increases the power consumption of servers in two ways: the consumption of cooling system (integrated fans) and leakage current. Further experimental study showed that, the contribution of leakage current is impressive. Taking a Gigabyte server as an example, even get rid of the integrated fan's power, server power can have a 16% rise by only raising CPU temperature. Fabrication process discrepancy can be another cause for the power variation observed. However, authors admit that they can hardly guarantee the same ambient temperature for each server in a rack, which constrains them from studying the variation introduced by fabrication process discrepancy. But recent studies altered that, the tiny fabrication discrepancy between the printed transistors can result in visible difference in terms of both performance and power consumption among high-performance microprocessors. The variation is becoming worse in modern processors [START_REF] Acun | Variation among processors under turbo boost in hpc systems[END_REF] [START_REF] Marathe | An empirical survey of performance and energy efficiency variation on intel processors[END_REF]. John C. McCullough et al. [START_REF] Mccullough | Evaluating the effectiveness of model-based power characterization[END_REF] found that when applying a power model trained on Intel Core i5-540M-1 to an identical processor 540M-2, mean prediction errors could be increased from 10% to 23%. They suggest using only power instrumentation for accurate power characterization. Among all of the components, processors are responsible for most of the power consumption and the variations [START_REF] Diouri | Your cluster is not power homogeneous: Take care when designing green schedulers![END_REF] [START_REF] Kistowski | Variations in cpu power consumption[END_REF]. Recent experimental studies have identified several sources of variation among which are: frequency variation introduced by advanced performance enhancement technologies such as Turbo Boost and Multi-Threading [START_REF] Marathe | An empirical survey of performance and energy efficiency variation on intel processors[END_REF] [START_REF] Acun | Variation among processors under turbo boost in hpc systems[END_REF], within die parameter [START_REF] Sarangi | Varius: A model of process variation and resulting timing errors for microarchitects[END_REF] and aging [START_REF] Diouri | Your cluster is not power homogeneous: Take care when designing green schedulers![END_REF].

However, the variation brought by different thermal features of processors has not been accurately characterized according to our knowledge. Addressing this lack of study, we compare the power variation for two processors with the same TDP but issued of different generations: Intel Xeon E5345 and Intel Xeon E5-2603v2. We test 30 identical samples for each type by switching processors in a same platform in order to eliminate the influences introduced by platform design. Thermal parameters such as ambient and CPU temperature are controlled and varied with the help of a climatic cabin and an external powered fan. The paper is organized as follows: Firstly, we introduce briefly the related work in Section II. In Section III, we explained how physical measurements are performed to compare the power consumption variation among identical processors. Then, in IV we propose and evaluate two hypotheses from the perspective of thermal characteristics: switch of TIM (Thermal Interface Material) and leakage current parameters. Conclusions and perspectives are given in V. Results of our experiments can help identifying and characterizing the influence of manufacturing variability on the power variation of processors. These observations remind us of reflecting correcting the precision of existing power predictive models. The major contributions of this paper are:

• We compare the power consumption variation for two generations of Intel processors. 30 identical processors are evaluated for each generation. • We explore two potential causes based on thermal characteristics hypotheses that may contribute to the power variation of identical processors: the switch of TIM and the variation of leakage current. Each hypothesis is validated in physical experiment. • Results of our experiments can help identifying and characterizing the influence of manufacturing variability on the power variation of processors. These observations remind us of reflecting correcting the precision of existing power predictive models.

II. RELATED WORK

In our previous work, we studied an energy predictive model for servers in using corresponding indicators, such as system usage, hardware configurations, power management technologies applied, as well as environmental information. Many studies relating to real word empirical measurement found that beside IT load applied to the components (CPU, memory, network and storage), the power of servers can be affected by external factors, such as original fabrication process [START_REF] Coles | Comparing server energy use and efficiency using small sample sizes[END_REF] [11], ambient temperature [START_REF] Orgerie | Demystifying energy consumption in grids and clouds[END_REF] [15] [START_REF] Sampath | Thermal analysis of high end servers based on development of detail modeland experiments[END_REF]. John C. McCullough et al. [START_REF] Mccullough | Evaluating the effectiveness of model-based power characterization[END_REF] found that when applying a power model trained on Intel Core i5-540M-1 to an identical processor 540M-2, mean prediction errors could be increased from 10% to 23%. They suggest using power instrumentation only for accurate power characterization. Marathe et al [START_REF] Marathe | An empirical survey of performance and energy efficiency variation on intel processors[END_REF] performed several tests to compare both performance and energy efficiency variation among identical nodes on Sandy Bridge, Ivy Bridge and Broadwell clusters. The variations are compared separately with and without hardware-enforced power limit. They found that processor performance and energy efficiency variation is becoming worse with the evolution of computation capacity on modern Intel processors. Balaji et al. [START_REF] Balaji | Accurate characterization of the variability in power consumption in modern mobile processors[END_REF] compare the power consumption variation for modern mobile processors. Their data shows power consumption variation among processors ranging from 5% to 17% when processors operate at the lowest and highest frequency respectively. Different power management settings such as Turbo Boost and C-state can also affect the value of variation. El Mehdi Diouri et al. [START_REF] Diouri | Your cluster is not power homogeneous: Take care when designing green schedulers![END_REF] find that different nodes from a homogeneous cluster have different power consumption at idle state. The power consumptions of two nodes stay unchanged even after exchanging the positions. They blame the causes to the age of the processors, as the server equipped with older processors shows more variation. Acun et al. [START_REF] Acun | Variation among processors under turbo boost in hpc systems[END_REF] investigate the processors under Turbo Boost in HPC systems. They point out that dynamic overclocking feature of processor is responsible for substantial frequency difference among the processors, which explains the up to 16% of core-to-core performance variation. The faster processors usually consume more than the slower ones. Jóakim et al. [START_REF] Kistowski | Variations in cpu power consumption[END_REF] characterize the variation on CPU power consumption. Experiments are performed on three different platforms and different processors are picked for each platform. Identical processor samples are exchanged after each run to guarantee the identical conditions. The power consumption can differ as much as 29.6% in idle and 19.5% at full load for identical samples. Their observations also show that CPU power directly influence system power, and most of the variation is caused by CPU power variations, other than by the other components [START_REF] Diouri | Your cluster is not power homogeneous: Take care when designing green schedulers![END_REF]. The authors use worklets in SERT (Server Efficiency Rating Tool) [START_REF] Lange | The design and development of the server efficiency rating tool (sert)[END_REF] as workloads to stress the SUT (System Under Test). To our knowledge, during the phase of calibration, SERT will identify the maximum rate at which transaction can be executed for each worklet. This value is highly reproducible for one server in run-to-run test, but may vary from server to server calibration [START_REF] Wang | Etude de l'influence des aspects thermiques sur la consommation et l'efficacité énergétique des serveurs[END_REF]. However, as the performance variation among processors has been observed in previous studies, and the authors has not mentioned the coreto-core performance variation reported by SERT in this paper, we cannot tell if the workloads (worklets after the phase of calibration) used to stress each sample of processor are exactly the same. Apart from frequency variability, S. R. Sarangi et al [START_REF] Sarangi | Varius: A model of process variation and resulting timing errors for microarchitects[END_REF] emphasize that, within-die parameter variation can result in process variation including both random and systematic effects, can also negatively impact a processor's frequency and leakage power. Most of the studies focus on the difference at the server level, since processors are placed in different servers or different sockets, the results cannot eliminate the influences caused by system noise, such as the influence of the nearby processors, the platform bugs [START_REF] Marathe | An empirical survey of performance and energy efficiency variation on intel processors[END_REF] [20] [START_REF] León | System noise revisited: Enabling application scalability and reproducibility with smt[END_REF]. Moreover, thermal control strategy is rarely mentioned. As processors are temperature sensitive components, different operational temperature can affect the results.

III. DO IDENTICAL PROCESSORS CONSUME THE SAME?

Newer technologies bring diverse features in improving performance for processors, and in the same time, introduce more variation regarding performance and energy usage to processors with the same design. We compare the power variation for Intel Xeon E5345 and Intel Xeon E5-2603 v2. Details about the processors can be found in Table I. The two processors both have 4 cores and the same TDP (Thermal Design Power) of 80W. Otherwise, they dispose different designs on architecture and operate on different base frequencies.

A. Experiments setup

The whole platform (excluding the power supply) is placed in a climatic cabin, where the ambient temperature is configured at 35 • C. The ambient set point is little higher than ASHRAE typical requirement limits (15-32 • C, for products require a stable and more restrictive environment) [START_REF] Committee | Data center power equipment thermal guidelines and best practices whitepaper[END_REF] in order to exposure as much as possible the thermal characteristics difference between samples. Fan is placed on top of the processor's heat sink and powered by an external DC (Direct Current) source so as to avoid power variation brought by motherboard's fan control. Homemade CPU intensive applications "bzip2" and "pi_calculator" are chosen to maximize the usage and heat dissipation of CPU. The application assign a specific task to the server, and active all the threads available in the system at full load to complete the task. This allows analyzing both the power and performance (execution time) variations among samples. The power consumption variation is analyzed by comparing the average power of the server in the last minute, where the thermal condition is tending towards stable. Performance variation. Samples of processors are tested and exchanged one by one in the same platforms. When exchanging the samples after each run, we try our best to uniform the heat paste applied between processor and fan. The platform test diagram shown in Figure 1 includes three parts: SUT, control and measurement system.

• SUT: Server Under Test, installed with Ubuntu Server OS. • The control system:

Controller: normal PC (personal computer) with Linux system installed, in the same local network with SUT and controls SUT remotely by SSH. Controller also gathers and synchronizes measurement data recorded by power and temperature analyzers. Climatic chamber: Servathin RC01, controls ambient temperature.

Fan: controls CPU surface temperature, powered by external DC source. • The measurement system: Power: measured by Wattsup-Pro [START_REF]Wattsup pro power meter[END_REF]. Measures contain the power consumption of the whole server except the fan's. Ambient temperature: measured by a thermometer with a thermocouple connected (type K, with 0.1mm diameter). CPU Temperature: collected by command line tool lm-sensors. Sampling frequency for all the measures are configured at 1Hz.

B. Power variation for Intel processors of different generations

After analyzing all the measurement data, we choose 18 representative processors for each platform: 9 samples that consume the most and another 9 samples that consume the least among all the 30 samples. Rank of consumption is shown in the following figures. Figure 2 and Figure 4 show the power by running application "pi_calculator" on the 18 processor samples of two types. It can be observed that, power variation among samples of Xeon E5-2603v2 (newer generation) between the best and worst can be more than 30% (16.1W / 51.9W). That is much larger than the variation observed from samples of Xeon E5345: 2.8% (2.8W / 98.4W). In addition, we verify the power variation of Xeon E5-2603v2 by running another application "bzip2", and the rank of power consumption for samples stays the same. Therefore, the power variation observed in Xeon E5-2603v2 is not supposed to be affected by workload. In addition, we can notice that the server power of all the samples of Xeon E5-2603v2 are below the TDP (80W, refer to Table I). Actually, TDP of CPU gives a guidance to design cooling system. It represents for the maximum amount of heat generated by CPU. The cooling system is designed to be capable to dissipate heat under any workload without exceeding thermal envelope. Thus, TDP is not reflecting the actual power of CPU, it can be lower or higher than actual power usage [START_REF] Andrew | The technical details behind intel's 7 watt ivy bridge cpus[END_REF]. Figure 3 and Figure 5 present the CPU temperature distributions along the test (around about 1 hour), values are calculated as the mean temperature of all 4 cores. Comparing to Xeon E5345, CPU temperatures of Xeon E5-2603v2 spread over a wider range. However, no obvious correlation can be seen between CPU temperature and power consumption.

We also compare the performance variation (execution time) among all the samples, little performance variations are observed for both Xeon E5-2603v2 and Xeon E5345, 0.3% and 0.6% respectively. We find this normal, according to previous studies, operating frequency is responsible for most of the performance variation. While none of them support Turbo Boost or Multi-Threading technologies, operating frequency remains stable at maximum rate when applying the same workload. Therefore, the power variation presented here has no relationship with the performance variation.

IV. HOW THERMAL EFFECTS IMPACT ON THE CONSUMPTION OF SERVERS?

Results of the experiments in Section III show that, processors of Xeon E5-2603v2 comes from newer generation have much more power variation than Xeon E5345. CPU temperature distributions along the test cannot perfectly explain the power variation. Even more, we confirm that the power variation is not introduced by the performance variation. Therefore, we can benefit from the situations to focus on the difference of thermal features between samples. The consumption of processor depends on IT load and cooling system. As the power increases with the load, processor works harder and dissipates more heat, if the heat is not evacuated in time by the cooling system, the temperature surrounding processor becomes higher and leads to the rise of leakage current, which will in reverse increase the power of CPU [START_REF] Kursun | Temperature variation characterization and thermal management of multicore architectures[END_REF] [START_REF] Moss | Energy impact of increased server inlet temperature[END_REF]. In the study of [START_REF] Wang | Potential effects on server power metering and modeling[END_REF], the authors observed that leakage current of processor can contribute to over 30% power variation for server. However, there is doubt whether the samples have the same ability to evacuate the heat. We are wondering if different thermal features can also result in the variation among processors. In this section, we propose two hypotheses regarding different thermal features between different samples of Xeon E5-2603v2. In subsection IV-A, we study the influence of TIM (Thermal Interface Material), in order to verify if the new PTIM (Polymer Thermal Interface Material) can result in the rise of the power consumption variation. In subsection IV-B, we investigate the parameters of leakage current among processors. We start from the representation equation of static power and analyses the relationship between the leakage current and the temperature. Fans, DC power generator and climatic cabins are provided to help controlling precisely the temperature.

A. Hypothesis 1 : The switch of Thermal Interface Material (TIM)

TIM is applied to fill the air gap between processor Integrated Heat Spreader (IHS) and silicon die, so as to dissipate heat produced by processor more efficiently [27] [28]. A thermal model proposed by Huang W et al. [START_REF] Huang | The need for a full-chip and package thermal model for thermally optimized ic designs[END_REF] shows that, the thickness variation of TIM can affect a lot silicon die temperature distributions across processors. Figure 6 demonstrates a simple structure of the thermal package design. Since Ivy Bridge processor generation, Intel decides to adopt Polymer TIM (PTIM) as a replacement of Solder TIM (STIM), for reasons of economics, environment and better cooling performance [START_REF] Mukul | Material technology for environmentally micro-electronic packagingl perspective, status, and future directions[END_REF]. Xeon E2603v2 applied with PTIM is turned out to have more power variation than Xeon E5345. We try to clear this suspect by removing PTIM from the original structure. However, TIM is hidden in IHS from exposing to outside, that arises a lot the challenge of manipulation: we have to remove TIM along with the entire IHS and this manipulation is irreversible. Three samples consuming less: P14, P23, P9 and two samples consuming more P30, P20 are involved in this study. After the manipulation, the processor's silicon die was cooled directly with the heat spreader of fan (as shown in Figure 7). We applied a little heat paste between them to facilitate the heat exchange. Then we repeat the test as described in Section III and re-analyze the power variation. Table II shows the results of power consumption and average CPU temperature in two circumstances: with and without PTIM. Unfortunately, removing the PTIM does not help correcting the power variation between samples. After removing PTIM, most of the samples (except P9) just have little power variation as well as with PTIM. Therefore, PTIM is not the source of power variation affected by fabrication process.

B. Hypothesis 2 : The leakage current variation

The overall power dissipation in today's microprocessors is composed principally by two sources: dynamic power and static power. Other power loss such as short circuit occurs Fig. 7. TIM remove schematic diagram at whenever gate switch is relative small, it can be absorbed by dynamic power [START_REF] Kim | Leakage current: Moore's law meets static power[END_REF]. In this sub-section, we firstly review the CMOS technology, as theoretical guideline for our experiments. Then we present the details and results of the test.

1) CMOS Technology review: Dynamic power results from charging and discharging the processor's capacitive loads. It can be described by the equation ( 1) [START_REF] Su | Low power architecture design and compilation techniques for high-performance processors[END_REF] [33]:

P dynamic = α C i f V 2 dd ( 1
)
where f is the switching activity (operating frequency), C i is the sum of gate and interconnection capacitances, and V dd is the supply voltage and α is the activity factor of the overall circuit. Static power is the product of voltage supply and leakage current. There are different kinds of leakage modes in MOS transistor, and the most dominant leakage mechanism is sub-threshold leakage I Dsub [START_REF] Keshavarzi | Intrinsic leakage in low power deep submicron cmos ics[END_REF]. I Dsub is the current flow between source and drain at off-state. Off-state current becomes now a limitation factor for down-scaling the threshold voltage, since it determines the power consumption of a chip in its idle state. Therefore, the static power dissipation representation can be simplified and described by equation ( 2). I Dsub can be described by equation ( 3), according to previous studies [START_REF] Horowitz | Low-power digital design[END_REF] [START_REF] Butts | A static power model for architects[END_REF].

P static = V dd I Dsub ( 2
)
I Dsub = e -qV th akT (3)
where q, k, a and k are physical related constants, T is the absolute temperature and V th is the threshold voltage of the transistor, which sits between ground and the supply voltage. In earlier years, traditional low-power microprocessor design focus mainly on reducing dynamic power consumption. At that time, static power consumption is not a limitation and is negligible compared to dynamic power [32] [35]. In pursuing higher performance and lower consumption, CMOS technologies scaled the chips from generation to generation by following Moore's law "The number of transistors and resistors on a chip doubles every 24 months" [START_REF] Moore | Cramming more components onto integrated circuits[END_REF]. The chips then have more transistors, denser CMOS circuitry and smaller dimension. In 2007, in applying 45nm process technology, there were 3.3 million transistors per square millimeter (MTr/ mm2) for chips of Intel. Ten years after, in 2017, Intel announced its latest chip generation: 10nm technology, with density as high as 100.8 MTr/ mm2, which is 30 times denser than in 2007 [START_REF] Rachel | Intel now packs 100 million transistors in each square millimeter -ieee spectrum[END_REF]. Note that, the designation like "45nm", "10nm" refers to commercial name for certain lithography process technology.

The number doesn't indicate the size of any particular feature of the chip and can vary significantly between manufactures [START_REF]Technology node[END_REF]. As a general rule, the smaller the number is, the denser the circuitry becomes.

As the chip dimension scaled, for the purpose of reliability required by constant field scaling, supply voltage V dd has to be decreased by the same factor for chips to keep the electric fields the same across different generations [START_REF] Dennard | Design of ion-implanted mosfet's with very small physical dimensions[END_REF]. That brings additional benefit of dynamic power saving as suggested by equation [START_REF] Jones | How to stop data centres from gobbling up the world's electricity[END_REF]. The threshold voltage (V th ), has also to be scaled down along with the V dd in order to avoid performance degradation as shown by equation ( 4) [START_REF] Sakurai | Alpha-power law mosfet model and its applications to cmos inverter delay and other formulas[END_REF].

Delay = 1 f = C gate V dd I Dsat ∝ V dd (V dd -V th ) 1.3 (4)
Yet, as a result, the I Dsub increases exponentially with the V th decreases [START_REF] Kim | Leakage current: Moore's law meets static power[END_REF] [41] according to equation [START_REF] Costa | Methodology of measurement for energy consumption of applications[END_REF]. Therefore, in modern processors, static power takes greater and greater part in consumption and becomes increasingly dominant. Moreover, advanced lithography process with thinner xnm may have more fabrication deviations to impact the parameter values of leakage current among identical processor samples, then leads to eventually power consumption variation. Comparing the two processors, the newer generation E5-2603v2 with 22nm technology has 10 times more power variation than E5345 with 65nm technology. This provides a great opportunity to empirically evaluate the impact of fabrication brought by CMOS technology evolution on the processors variability.

However, there is no way to measure directly the leak current of the processor by physical device, without "opening" some key components hidden and protected by Intel packing technology. Such manipulation is too risky and may cause permanent damage to the motherboard. The only measure of consumption accessible is the whole consumption of server composed of consumption of CPU (P cpu ) and of the other components including motherboard (P chip ). The processor samples are switched one by one in the same motherboard, P chip remains the same, the only item can vary the consumption of server is P cpu . As discussed before, P cpu is composed by dynamic and static power. It can be seen from the equation (1) (2) and (3) that, dynamic power is frequency dependent value but independent from temperature variation. On the contrary to dynamic power, static power is temperature sensitive value, but cannot be affected by frequency scaling. The consumption of the whole server can be then simply represented by the equation ( 5) and ( 6):

P server (V dd , T cpu , f ) = P chip + P cpu (V dd , T cpu , f ) (5)
P cpu (V dd , T cpu , f ) = P dynamic (V dd , f ) + P static (V dd , T cpu ) (6) Our test platform has no support in BIOS to regulate manually V dd . f is governed by frequency governor "performance" and adjusted by frequency driver intel_pstate. f varies according to system current load. In this case, the data obtained is too limited to separate and quantify the two sources of consumption by statistical analyses. Deriving the models of static and dynamic power of processor is interesting, but it is beyond the scope of the study discussed here. If interested, Goel et al [START_REF] Goel | A methodology for modeling dynamic and static power consumption for multicore processors[END_REF] present a systematic methodology for modeling static and dynamic power consumption of individual cores and uncore components in their work. In our cases, for a given processor sample, we fix f and V dd , and suppose that the value of P dynamic is constant and independent from temperature variation. The parameters of leakage current as expressed in equation ( 2) and (3) can be simply identified by varying the T cpu whiling keeping the f and V dd unchanged.

2) Experiment setup: Same platform as described in section III-A is adopted. CPU intensive benchmark cpuburn [START_REF] Robert | Ubuntu manpage: cpuburn, burnbx, burnk6, burnk7, burnmmx, burnp5, burnp6 -a collection[END_REF] is selected to stress the SUT in this test. Comparing other CPU intensive benchmarks, cpuburn is better at maximizing the heat dissipation of processor. During the execution of cpuburn, each thread is stressed at full load and operating frequency (f ) is then maintained constantly at (1.8GHz). Besides ambient and CPU temperature, we monitor as well the CPU voltage (in using command-line tool lm-sensors [START_REF] Archwiki | lm sensors -archwiki[END_REF]) and CPU frequency (in using python library psutil [START_REF] Psutil | [END_REF]). Climatic cabin is used to vary the surface temperature of CPU. Some of the processors samples have been used in previous test to validate the hypotheses 1 in section IV-A. Their TIM and HIS are removed. This manipulation is irreversible and may introduce unexpected influence to the consumption, therefore, these samples are excluded from this test. Four samples consume less: P12, P21, P6, P4 and four samples consume more P3, P10, P5, P8 are chosen to validate the second hypotheses.

For each sample, we repeat the test procedures as follows:

• Keep the server on idle state for 3 minutes at 22 8 shows the evolution of CPU temperatures with the ambient temperature increases from 22 • C to 50 • C. It can be observed that, cooling abilities of CPU samples differ from each other. CPU temperature increases gradually with the rise of ambient temperature but ends at different values where the ambient temperature is supposed to be at 50 • C. Figure 9 shows how server power varies with the rise of CPU temperature. Processor samples are marked with different colors. The points represent for the real measure data. The server power can be represented by CPU temperature in using an exponential function of form exp(-k/T ) as mentioned in equation [START_REF] Costa | Methodology of measurement for energy consumption of applications[END_REF]. We can find that the change rates of power with the increase of CPU temperature are different among samples. Samples consume more (P3, P10, P5 and P8) have also higher change rates than samples consume less (P12, P21, P6, P4). In another word, processor samples P3, P10 P5 and P8 are harder to be cooled down than the others. Change rates of samples can be identified clearly in presenting the exponential function without the intercept, as shown in Figure 10 b) On idle state: On idle state, only OS is running in the server. Power consumptions on idle state for different samples are compared at 22 • C and 50 • C ambient temperature respectively. Results are shown in Figure 11. We can observe that, there is more power variation on idle state at 50 • C ambient temperature. Moreover, with the same rise of ambient temperature, idle power goes up to a higher value for samples that consume more.

V. CONCLUSIONS AND PERSPECTIVES

With the decrease of lithography size for modern processors, technology today can hardly control precisely the variability between processor samples in fabrication process. Apart from the performance variability brought by the operating frequency variation, our experimental results show surprisingly high variability of power consumption within modern processor. 30% power consumption variation is observed among 30 identical processor samples, which has no correlation with the performance (frequency variation). In the first time, we confirm and characterize by means of physical experiments, the power consumption variation introduced by within die parameter variation. Fabrication variability has random effects on processor samples. Samples can have different heat dissipation rates under the supposed same conditions (ambient temperature and load level), which results in different leakage current distributions then finally affect the static power consumption among samples. Our demonstration highlights the challenges of modeling techniques posed by rising processors variability. Even for servers in homogeneous clusters cannot consider being exactly the same. As a result, the precision of existing power models can be questioned when applying to other identical SUT samples. Taking these concerns into account, power modeling techniques based on pure software solutions such as linear regression or performance counters, are not reliable enough for modern servers in predicting power from one server to several servers in a rack or data centers. In addition, the findings present in this paper can be also applied to optimize the energy management strategies in data center, such as VM migration, shut down technology. Power consumption variation occurs at server level could be worse than at processor level. Processor samples with worse thermal features will not only have more static power, but also request more energy in associated cooling system. We will further investigate to filter out the servers with both better and worse thermal features in a data center, in order to better orient energy management strategies. Processors are becoming smaller, more powerful and less consuming from generation to generation. On the other side, processors become more complex than ever. Leakage current variations arise from imperfections in the fabrication process among modern processors, such as lithographic length aberration. According to our observations, leakage current variations have impact on cooling ability, samples dispose of different increase rate with the rise of temperature.

Fig. 1 .

 1 Fig. 1. Platform test diagram

Fig. 2 .

 2 Fig. 2. Server power comparison for samples of Xeon E5-2603v2

Fig. 4 .Fig. 5 .

 45 Fig. 4. Server power comparison for samples of Xeon E5345

Fig. 6 .

 6 Fig. 6. Processor thermal package structure

Fig. 8 .

 8 Fig. 8. CPU temperature as function over time

Fig. 10 .

 10 Fig. 10. Illustrations of server powers functions for different processor samples

TABLE II POWER

 II AND CPU TEMPERATURE FOR XEON E5-2609V2: WITH AND WITHOUT PTIM (WORKLOAD: PI CALCULATOR)

	Processors With PTIM	Without PTIM	∆ (Watt & • C)
		P(W) T( • C) P(W) T( • C) ∆ P(W) ∆ T( • C)
	P14	50.0	38.2	50.3	40.0	+0.3	+1.8
	P23	50.6	40.3	51.0	36.0	+0.4	-4.3
	P9	50.6	39.5	56.4	36.7	+5.8	-2.8
	P30	63.9	42.9	63.3	38.9	-0.6	-4
	P20	66.8	42.7	68.3	47.7	+1.5	+4.7

  Execute cpuburn for 30 minutes. Along the test, configure the climatic cabin to increase the ambient temperature from 22 • C to 50 • C. • After the execution of cpuburn, keep the server on state idle for 3 minutes at 50 • C ambient temperature. • Shutdown the server, switch for the next processor sample and wait for the whole platform cool down at 22 • C ambient temperature. 3) Experiment results and analyzes: a) On high load state: High load state occurs at the second step, where server is stressed by cpuburn. Figure

• C ambient temperature. •