Modeling of inter- and transgranular stress corrosion crack propagation in polycrystalline material by using phase field method

1 Introduction

Stress corrosion cracking (SCC) is a very common cause of failures for the engineering components and structures. At the microstructural scale, cracks often nucleate from corrosion pits and, depending on the kind of material and environment, it could be inter- or transgranular (IG-TG) in nature [1]. Especially, crystals orientation, grain boundary properties and mechanical loading strongly affect the switching between the inter- and transgranular stress corrosion cracking (IGSCC-TGSCC) [2], [3]. This competition provides a significant difference in the material resistance [4], [5], [6]. Therefore, understanding these aspects plays an important role for improvements in the development and application of advanced structural materials.

Many works in the literature have been proposed to deal with the simulation of microstructural effects on the stress corrosion assisted fracture growth in polycrystals. Cohesive Zone Models (CZM), based on the original idea in the work of Dugdale [7] and Barenblatt [8] to avoid the singularity of stress at the crack tip, is a popular technique used in this field. One of its application can be found in the study [9], where a grain-scale finite element (FE) model based on the CZM is proposed to study IGSCC in stainless steel wire. In the work [10], a 3D grain-level cohesive finite element model is proposed to study the IGSCC characteristics of AZ31 magnesium alloy due to the hydrogen embrittlement process. Another contribution based on cohesive zone domains can be found in [11] by using multiscale model to simulate the intergranular hydrogen-assisted cracking. The investigation of the microstructural effects (especially the grain boundaries) on the behavior of the IG stress corrosion assisted short cracks can be found in [12] for the alloy 600 by using the body force method, and [13] in the fully sensitized austenitic stainless steel by using the polycrystalline aggregate model, or in [14] by using a new model based on Markov Chain theory and Monte Carlo simulations. The polycrystalline model is also used in [15], [16] to simulate the IGSCC in stainless steel. Another technique, namely dislocation-based boundary element model is proposed in [17] to microstructurally study intergranular stress corrosion short cracks in an elastoplastic medium.

Development of the numerical models at the grain size scales has so far been mostly limited to simulate IGSCC. There is a very few studies dealing with both IGSCC and TGSCC such as the work of Musienko and Cailletaud [6]. The authors proposed a coupled approach (environment and plasticity) based on the CZM to model the IG-TG fracture transition. More recently, a numerical multiphysics peridynamic framework [18] based on the adsorption induced decohesion mechanism is proposed to this field. However, a detailed description of the electrochemical and mechanical effects on stress corrosion assisted damage is still missing, especially for the common situation of anodic dissolution induced SCC.

In this work, the phase field method augmented by a smeared description of interfaces is used to develop a fully coupled multiphysics model of stress corrosion assisted damage growth [19], [20], [21]. Starting from the ideas in [22], [23], we extend them to the case of polycrystalline materials, where the anisotropic effects for both elastic energy and fracture surface energy are taken into account. The behavior of grain boundaries is described by the smeared cohesive model following the work [24], [25]. Moreover, to account for the different behavior of grain boundary and grain interior, we propose an interpolation formulation to present the smoothed transition of the electrochemical-mechanical properties between the two regions. It is expected that the new model can model both IGSCC and TGSCC depending on the electrochemical-mechanical loading.

The overview of the paper is as follows. In Section 2 the fundamentals of the proposed model are presented. In Section 3, an extented formulation of phase field method able to model the stress corrosion induced crack propagation accounted the crystal effects is proposed. The variational principle for the proposed model is then detailed in Section 4. Finally, the potential of the proposed model is illustrated by several numerical examples in Section 5.

2 Fundamentals of the proposed method

Let Ω⊂ℝ^D be an open domain describing a corroding system at micro scale, which generally contains four domains: (i) metal or solid phase described by temporally and spatially constant metal atom concentration csalloy; (ii) electrolyte solution or liquid phase described by temporally and spatially constant electrolyte atom concentration clalloy; (iii) an interfacial region where the corrosion occurs with the atom concentration taken within the range (clalloy, csalloy); (iv) and the grain boundaries region described by metal atom concentration csalloy as well, but with another behavior.

The normalization procedure is then adopted following the work [22]. We define c_s as the normalized solid concentration and c_l as the normalized electrolyte concentration, reading:

The corrosion process is modeled by using the phase field method. A scalar parameter ϕ(x) is considered to describe the corroding system. ϕ(x)=1 is solid phase, ϕ(x)=0 is liquid phase and the interfacial region is within the range 0<ϕ(x)<1. Moreover, the grain boundaries are also approximated by a smeared representation defined by a temporally constant β(x), taking a unit value on grain boundaries and vanishing away from it.^[1] An illustration of a two-dimensional corroding system is depicted in Figure 1.

Figure 1:

2D representation of a typical corroding system due to the anodic dissolution of a polycrystalline body.

The dissolution of material (phase transformation: solid ϕ=1 to liquid ϕ=0) is provided by the diffusional transport of metal ions from solid domain to liquid domain at interfacial domain. Note that, the interfacial region (transition domain) is represented here as a mixture of both phases (solid–liquid) with different chemical compositions by using the assumption of (Kim Kim Suzuki) KKS model [26]. Then, by introducing the interpolation function, satisfying h(ϕ=0)=0, h(ϕ=1)=1 and ∂h(ϕ)∂ϕ|ϕ=0,ϕ=1=0; one choice: h(ϕ)=−2ϕ³+3ϕ², the ions concentration at any material point c can be described by the following:

As mentioned in the work [23], the phase transformation occurs when the ions concentration achieves the saturation value, denoted c_se and c_le for solid metal and electrolyte, respectively, in the normalized framework.

3 Phase field model

This work focuses on the study of polycrystalline effects to the anodic dissolution induced stress corrosion crack growth. To this aim, we extend the formulation proposed in [23] to the case of polycristalline materials. The new energy formulation is rewritten as follows:

where ψ_u(ε(u)) is the strain energy density, ψ_c(c, ϕ) is the electrochemical energy density and ψ_ϕ describes the fracture energy. The term ψ_v(〚u〛, κ) is added to account for the energy of grain boundaries, which is a function depending on the displacement jump across grain boundaries Γ_B and an history parameter κ.

In the phase field framework, the regularized description for strong discontinuities related to both fracture and grain boundaries is adopted (substituting 〚u〛 by v(x)). A sharp grain boundary is regularized by a smeared grain boundary (substituting Γ_B by Γ_β=∫_Ωγ_β(β, ∇β)dΩ). Then its decohesion will be described by the smeared cohesive law proposed in [25]. This approach is based on a smeared displacement jump approximation (see [24], [27]), in which the grain boundary decohesion is approximated as a smooth transition v(x). That provides a decomposition of infinitesimal strain tensor into a part related to the bulk and a part induced by smoothed jumps, denoted by ε^e and ε˜, respectively.

with

where nΓβ is the normal vector to Γ_β.

The strain energy density of damageable material ψ_u(ε^e(u), ϕ) is then constructed from elastic part ε^e by using unilateral contact assumptions introduced in [24]:

with the elastic tensor ℂ(ϕ) is defined by the following:

where ℂ⁰ denotes the initial elastic tensor of the material, possibly anisotropic; k₀ is the bulk modulus; the sign function sign⁻(x)=1 if x<0 and sign⁻(x)=0 if x≥0.

The energy formulation in the phase field framework is now defined by:

The third term in Eq. (8), ψ_ϕ(ϕ, ∇ϕ) describes the surface energy or the amount of energy released upon the creation of new fracture surfaces. It is chosen following the work [23]:

with ω_ϕ is the height of the imposed double-well energy barrier, α_ϕ is the second-order tensor describing gradient energy (included the anisotropic effects), which is an invariant with respect to rotations:

where α₀ is the classical gradient energy coefficient for the case of isotropic material; 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. This formulation can be extended for the case of many preferential cleavage planes. The detailed description can be found in the work [28].

Note that ω_ϕ and α₀ can be determined from the relation with the transition length ℓ and the interface energy σ_ϕ or the critical fracture energy g_c:

where α^*≈2.94 corresponds to the transition zone taken within the range 0.05<ϕ<0.95.

The electrochemical energy density ψ_c(c_s, c_l, ϕ) in Eq. 8 is chosen based on the KKS model [26], wherein the transition region is postulated as a mixture of the two phases with different chemical compositions, but equal chemical potentials:

where A is the free energy density curvature which is often identified from thermodynamical databases [29], [30]; c_se and c_le are the equilibrium concentrations of the coexisting phases; the concentration of solid c_s and liquid c_l phases can be expressed from the equilibrium condition at the transition zone as follows:

The total energy is then rewritten as E(u, v, ϕ, c)=∫_ΩψdΩ in which

is identified total energy density.

4 Variational principle

4.2 Variational principle for electrochemical problem

The Cahn-Hilliard and Allen-Cahn equation is used to described the electrochemical processes of phase transformation. The diffusion of ions from solid volume to electrolyte, or the spatial evolution of conserved field c is governed by a conservation law (Cahn-Hilliard). The material dissolution or the evolution of non conserved field ϕ is governed by Allen-Cahn equation:

(21){∂c∂t(x, t)=∇⋅M(∇∂ψ∂c),∂ϕ∂t(x, t)=−Lϕ(∂ψ∂ϕ−∇⋅∂ψ∂∇ϕ),

where L_ϕ is the interface kinetics parameter, that describes the transformation rate; M is the diffusion mobility, which can be evaluated based on an analogy between the Cahn-Hilliard equation and the Fick’s second law (see [22], [23], for more details), being written as:

(22)M=D2A,

with D is the diffusion coefficient. An important property of polycrystalline materials needs to be considered is the different SCC phenomena between grain boundary and grain interior. The grain boundary atoms are more easily and rapidly dissolved, or corroded, than the atoms within the grains. So that, the grain boundaries oxidize or corrode more rapidly. To model this phenomenon, the diffusion coefficient D is here defined as a function of the smeared description of grain boundary β by the following:

(23)D=DG(1−β)n(1−kD)+kD

where D_G is the diffusion coefficient for grain interior; n is an exponent; k_D is the diffusion mismatch ratio between grain interior and grain boundary. This definition ensures D=D_G in the grain interior, D=D_B=D_G/k_D in the center of grain boundary and D_G<D<D_B in the smeared region. An illustration of the variation of the diffusion coefficient crossing the grain boundary for the case of n=2 and k_D=0.1 is depicted in Figure 2. This approximation can be understood as an interpolation procedure for two different regions. By using the same principle, the fracture energy for whole domain g_c can be written:

Figure 2:

Explicit description of the different electrochemical behavior between grain boundary and grain interior.

(24)gc=[(1−β)n(1−kg)+kg]gcG

where gcG is the fracture energy in the grain interior, k_g is fracture energy mismatch between the grain boundary gcB and the grain interior gcG, defined by kg=gcB/gcG.

Using the variation in the concentration for Eq. (21)₁, and the variation in the phase field for Eq. (21)₂, the corresponding weak form was obtained^[2]:

(25){∫Ω∂c∂tδcdΩ+∫ΩM∇∂ψc∂c⋅∇δcdΩ=0∫Ω∂ϕ∂tδϕdΩ+∫ΩLϕδϕψdΩ=0

where δ_ϕψ reads

(26)δϕψ=(∂ψ∂ϕ−∇⋅∂ψ∂∇ϕ)δϕ=[∂ψc∂ϕ+[εe:∂ℂ(ϕ)∂ϕ:εe] +ωϕg′(ϕ)]δϕ+∇ϕ[α0(1+ξ(1−N⊗N))]∇δϕ.

In order to handle loading and unloading histories, we introduce the history function following the work of [31], [32]

(27)ℋ(x, t)=maxτ∈[0,t]{εe(x, τ): [ℂ0(x)−k0(x)1⊗1sign−(trεe(x, τ))]:εe(x, τ)}

From (26), (27), the final weak form for electrochemical processes is obtained by the following equation:

(28){∫Ω∂c∂tδcdΩ+∫ΩM∇∂ψc∂c⋅∇δcdΩ=0∫Ω∂ϕ∂tδϕdΩ+∫ΩLϕ[∂ψc∂ϕ+h′(ϕ)ℋ(x, t)+ωϕg′(ϕ)]δϕdΩ+∫ΩLϕ∇ϕ[α0(1+ξ(1−N⊗N))]∇δϕdΩ=0

5 Representative numerical examples

In this section, we present several numerical examples using the proposed method to study SCC of a zirconium alloy. Thanks to a tightly adherent and protective oxide film, zirconium is well known for its excellent corrosion resistance to a lot of acids (hydrochloric acid, hot organic acids, sulfuric acid,…) even when exceeding the normal boiling point temperature. Moreover, this material exhibits very good mechanical properties. However, zirconium can be corroded in the situation of hydrochloric acid containing oxidizing species such as cupric chloride, ferric chloride, or wet chlorine. Both Intergranular and Transgranular Stress Corrosion Cracking (IGSCC – TGSCC) were experimentally observed in the literature [6], [33].

In this study, anodic dissolution induced stress corrosion crack propagation of this material is investigated. The considered zirconium alloy is composed mostly of zirconium, tin, iron and chromium with a typical mass density of ρ_alloy=6.56 g/cm³. The complete chemical composition is given in Table 1. The metal concentration is finally obtained csalloy=72.23 mol/l (see the work [23] for more details).

The diffusion coefficient is assumed to be homogeneous for all chemical compositions, and taken as the value for pure Zr, with D_G=1.3×10⁻⁸ m²/s [34], [35]. The saturated concentration is taken as in the study [36]: c_sat=3 mol/l. Hence, the saturated concentration of solid and liquid phase are c_se=1, c_le=3/72.23. The free energy density curvature A is identified by fitting the chemical free energies obtained from thermodynamic databases. That gives A=8.02×10⁸ J/mol. The size of transition zone and the interface energy are chosen ℓ=0.5 μm and σ_ϕ=0.5 kN/m (giving the equivalent fracture energy gcG=10−3 kN/mm). The interface kinetics energy is L_ϕ=0.15 l/(J.s).

The grain boundary cohesive energy is taken following the work of Xu and Needleman [37]:

The normal traction separation law is obtained as:

with vn=gcB/(tuexp(1)), and t_u the fracture strength, which can be dependent on the angular mismatch between grains (see e.g. [38]) as follows:

where tuavg, Δt_u are the average and the maximal fracture strength deviation, respectively; the angular mismatch Δθ is defined as the difference between the orientations of two crystals θ1i, θ2i, along the considered grain boundary on plane i:

However, for the sake of simplicity and due to the lack of experimental data, the grain boundary fracture strength is considered here as independent for misorientation between two grains and taken as t_u=0.1 GPa. The fracture energy within the grain boundary is gcB=2.85×10−4 kN/mm. The elastic tensor of hcp symmetry is used. For 0° orientation, it is written as follows [39]:

The bulk modulus is defined as k0=(C110+2C120)/3=96.26 GPa. For easier numerical implementation, all parameters are normalized following the work by [29], [30].

We first consider a rectangular domain with dimension 40×30 μm² containing an initial crack and 10 grains. The geometry and mesh of this microstructure is generated using the open source software Neper [40] (developed by Romain Quey at CNRS and MINES Saint-Etienne, France). A refined mesh was constructed using triangular elements, with element size ~0.15 μm, to have 6–8 elements in the transition zone. The average grain size is about 40×3010≈11 μm. The details of the geometry and of boundary conditions are described in Figure 3. In each grain, both direction of anisotropy and preferential damage directions are generated randomly.

Figure 3:

Geometry and boundary conditions for the polycrystalline structure containing 10 grains: (A) grain boundary and loading description; (B) description of crystals with angle of the direction of anisotropy.

In the crystallographic orientation, the cleavage plane N=[1 0] is chosen. The structure will be numerically assigned to both electrochemical degradation and mechanical loading. The displacements are prescribed along the y–direction for upper edge (y=30 μm) and lower edge (y=0) while the displacement along x is here free. The rupture of passive film allowing electrochemical processes is manually activated from the beginning at the initial crack location. Plane strain condition is assumed.

The competition between IGSCC and TGSCC is analyzed by considering two situations:

1. SCC with low mechanical loading.

2. SCC with high mechanical loading.

The first case (C1) is set up to expect the IGSCC, because at low mechanical loading, SCC process, taking place preferentially at grain boundaries, prevails on mechanical process. Therefore IGSCC will be dominant in this situation. In the second case (C2), when the structure is subjected to higher mechanical loading, the mechanical process becomes dominant. So that, the cracks will propagate not only on grain boundaries but also within the grains depending on the direction which is suitable for fracture growth. The details of boundary conditions are given in Table 2.

The resulting crack propagation for the first case (C1) is depicted in Figure 4. The initial crack propagates towards the grain boundaries, according to the preferential direction (Figure 4 with U=0.21 μm). The strong influence of SCC is captured in this loading condition. The fracture width near the edge x=0 is much larger than the region nearly crack tip, because of material dissolution induced by SCC had more time to develop in this zone. After reaching the grain boundary, the crack has a tendency to propagate along the grain boundaries (Figure 4 with U=0.27, 0.30 and 0.32 μm). IGSCC is clearly observed.

Figure 4:

Crack propagation for the first situation (C1) corresponding to the low mechanical loading. The IGSCC is clearly captured. (A) U=0.21 μm (t=115 min); (B) U=0.27 μm (t=205 min); (C) U=0.30 μm (t=255 min); (D) U=0.32 μm (t=285 min).

A different fracture scenario is obtained in the second case (C2) (see Figure 5). With the higher loading rate, the mechanical process becomes predominant. The fracture morphology is shaper, but the effect of SCC is still captured at the edge x=0. The fracture growth in the grain boundary is always observed (see Figure 5 with U=0.28 μm), but then the crack in this situation mostly propagates within the grains (see Figure 5 with U=0.35 and 0.39 μm). Under the higher mechanical loading, fracture growth is essentially due to mechanical process, that supports for both IG and TG cracks (not like only IG fracture in loading condition C1). Thus cracks choose either IG or TG, depending on which one is more favorable energetically. This phenomenon has been proved in Figure 5 for U=0.41 μm, where the orientation of the last grain is not suitable for crack propagation, the fracture has switched from TG to IG. The obtained numerical simulation results show a very good qualitative agreement with experimental observations in the literature [6], [41], [42].

Figure 5:

Crack propagation for the second situation (C2) corresponding to the high mechanical loading. Both IGSCC and TGSCC are observed. (A) U=0.28 μm (t=95 min); (B) U=0.35 μm (t=150 min); (C) U=0.39 μm (t=220 min); (D) U=0.43 μm (t=285 min).

Note that the transition between IG and TG fracture is strongly affected by the interfacial properties at grain boundaries as well. By using the same mechanical loading condition as the case (C2) but taking a new grain boundary energy gcB=8.3×10−4 kN/mm (denoted here the case (C3)), TG crack propagation is obtained as in Figure 6C. The crack trajectory of this last case (C3) is quite similar to one obtained from the loading case (C2), but with lower grain boundary effects. Especially at the last loading step, the fracture switches from TG to IG for the case (C2), while it is still TG for the case (C3). This result confirms the influence of the fracture energy mismatch between grain boundaries and grain interior on the transition between IGSCC and TGSCC.

Figure 6:

Comparison of crack path for different cases: (A) case study (C1) with the low mechanical loading and gc/gcB=3.5; (B) case study (C2) with the high mechanical loading and gc/gcB=3.5; (C) case study (C3) with the high mechanical loading and gc/gcB=1.2.

The stress-displacement curves for all situations are plotted in Figure 7. The elastic period for the three cases exhibits the same behavior. Then, the IGSCC in the case (C1) shows the weakest resistance, while the fully TGSCC in the case (C3) presents the strongest resistance. These observations demonstrate the strong influence of the fracture features (IG or TG) on the macroscopic material properties, which plays an major role in the material design.

Figure 7:

Comparison of the stress-displacement curve for both different mechanical loading conditions.

In the next example, we perform the numerical simulation on a polycrystalline structure containing 50 grains (dimension: 160×120 μm²) is performed in order to demonstrate the efficiency of the proposed framework. The details of geometry and boundary conditions are depicted in Figure 8. The same values of the material parameters are considered as for the previous example. The boundary conditions are chosen here to capture both IGSCC and TGSCC as follows: the prescribed displacement rate V=0.1 μm/min for the period t=[0:5 min] and V=0.002 μm/min for t=(5:500 min]. Only the cleavage plane N=[1 0] is considered.

Figure 8:

Geometry and boundary conditions for the polycrystalline structure containing 50 grains: (A) grain boundary and loading description; (B) angle of the direction of anisotropy.

The results of fracture evolution corresponding to different loading times are plotted in Figure 9. Both crack in the grain interior and along the grain boundaries are observed. Note that the advantage of applying the phase field model is demontrated, the model is able to capture the nucleation of severals cracks when the structure is subjected to predominant mechanical loading (see Figure 9 with U=1.40 μm). These cracks often nucleate at the grain boundaries (the weakest regions) and then either propagate along grain boundaries or onset and penetrate the grain (see Figure 9 with U=1.42 μm), so that the process requires minimal energy. Then these new cracks may merge with the initial one.

Figure 9:

Crack propagation for polycrystalline structure containing 50 grains. The phase field is plotted corresponding to several loading steps. (A) U=1.10 μm (t=305 min); (B) U=1.31 μm (t=410 min); (C) U=1.40 μm (t=455 min); (E) U=1.42 μm (t=470 min).

A similar behavior has been observed experimentally in the SCC crack growth in stainless steels exposed to high temperature water [43]. The process is commonly cited as initiation dominant growth [44] and has been related to the nucleation of adjacent cracks in the vicinity of the primary crack tip (stress concentration zone), so the main defect become larger by coalescence with the recently nucleated cracks rather than by the direct growth of the primary crack.

6 Conclusion

In this contribution, a new multiphysics model is proposed to describe the transition between IGSCC and TGSCC at the microstructural level. The anodic dissolution induced SCC is modeled by using the anisotropic phase filed method based on KKS model [26] for the electrochemical process and the idea in the work [23] for coupling with mechanical processes. The anisotropic effects of oriented crystals are taken into account for both elasticity and fracture energies. Moreover, the behavior of grain boundaries are described by the smeared cohesive model in the phase field framework to avoid defining of new material phases for this region. More specifically, we propose an approximation to account the different SCC phenomena between grain boundary and grain interior by using the diffusion coefficient and fracture energy of phase concentration depending on a smeared grain boundary field. The new method is able to simulate the stress corrosion crack propagation accounting for polycrystaline effects in a very efficient manner.

The applications to Zirconium alloy demonstrate a very good qualitative agreement between numerical predictions and experimental observations for both situation: IGSCC and TGSCC when SCC system is subjected to different loading conditions. The transition from IG to TG fracture due to the distribution of crystallographic and fracture energy mismatch between grain boundaries and grain interior is also captured. The proposed method opens a new efficient way to model the interaction and competition between IGSCC and TGSCC.