A phase field model to simulate crack initiation from pitting site in isotropic and anisotropic elastoplastic material

A multiphysics phase field framework for coupled electrochemical and elastoplastic behaviors is presented, where the evolution of complex solid-electrolyte is described by the variation of the phase field variable with time. The solid-electrolyte interface kinetics nonlinearly depends on the thermodynamic driving force and can be accelerated by mechanical straining according to the film rupture-dissolution mechanism. A number of examples in two- and three- dimensions are demonstrated based on the finite element-based MOOSE framework. The model successfully captures the pit-to-crack transition under simultaneous electrochemical and mechanical effects. The crack initiation and growth has been demonstrated to depend on a variety of materials properties. The coupled corrosion and crystal plasticity framework also predict the crack initiation away from the perpendicular to the loading direction.


Introduction
Structural components are often simultaneously subject to both mechanical loading and harsh environments in aerospace, biomedical, and marine applications. Environmental influence plays a significant role in the failure and fracture mechanism of the components in these applications [1,2]. The combination of mechanical loading and corrosive environments facilitate the nucleation and growth of cracks [3][4][5]. A range and complex mechanisms may involve in this process, including local metallic dissolution and rupture of the protective passive film [6,7], and uptake of pernicious species such as hydrogen generated by cathodic reactions during the corrosion process [8]. The need for capturing multiple chemical and mechanical phenomena, and their interplay, is one of the most major challenges in predicting environmentassisted failure [9,10]. Another important challenge is the interfacial nature of the corrosion and fracture problem, that requires tracking the evolution of corrosion front and/or crack propagation at the solid/environment interface [11][12][13].
Phase field-based method has recently received significant attention in addressing the evolution of complex geometries and have been applied to various problems including solidification, solid phase transformation, crack propagation, etc [14,15]. In this framework, the free energy of the system consists of contributions from chemical, interfacial, mechanical terms, etc. Thus, the phase field methodology developed in a thermodynamically consistent manner can be conveniently extended to coupled multiphysics phenomena. The dynamic process is driven by the reduced total free energy. Recently, phase field-based method has been applied to corrosionrelated phenomena [16][17][18][19][20][21][22][23][24]. Mai developed a Kim-Kim-Suzuki (KKS)-based model to simulate the activation-controlled, diffusion-controlled and mixed-mode pit growth kinetics [17] and further extended the framework to dissolution driven crack initiation from pitting [18]. Ansari developed a multiphysics phase-field model by incorporating the diffusion of chemical species and charge conservation to simulate pitting corrosion with [19] and without [20] insoluble corrosion product, and the intergranular corrosion [21]. A phase field model has also been developed by Martínez-Pañeda to understand the hydrogen assisted fracture and the unstable crack growth in the presence of hydrogen penetration [22]. The microstructure (polycrystalline structure) on the heterogeneous pitting corrosion has been discussed by Brewick [25,26] and applied to additively manufactured 316L stainless steels [27]. There have been relatively few studies to understand the transition from corrosion pits to initial crack under combined corrosion and mechanical effects using phase field framework. Mai discussed the dissolution driven crack initiation from initial elliptical pitting [18]. Cui investigated the pitting and stress corrosion cracking based on the film rupture-dissolution-repassivation mechanism where mechanical straining is proposed to accelerate corrosion kinetics [23], and later developed a framework capable to consider the hybrid effects of dissolution-driven and hydrogen embrittlement [28]. However, the elastic and plastic anisotropy of the materials are mostly not considered in these studies.
In this work, we aim to develop a multiphysics phase field model to simulate the crack initiation from a semi-circular initial pitting site under mechanical loading, that is capable to describe the coupled electrochemical and mechanical behaviors. A nonlinear dependence of the interface kinetics on the electrochemical driving force will be assumed in the electrochemical process. The mechanical behavior takes into account of both the isotropic and anisotropic materials, where the isotropic and anisotropic mechanical behaviors are described by J 2 plasticity and crystal plasticity, respectively. The pit-to-crack transition is simulated by coupling the mechanical and electrochemical process through film-rupture-dissolution mechanism. The rest of this manuscript is structured as follows. The theoretical formulation of the phase field based multiphysics model that governs the evolution of the metal-environment interface will be discussed in section 2. The theoretical model will be implemented using the open-source finite element based MOOSE (Multiphysics Object-Oriented Simulation Environment) [29] framework. Section 4 demonstrates the capacities of the model through two representative benchmark problems followed by a number of two dimensional (2D) and three dimensional (3D) numerical examples, where the alloy will be first modeled as an isotropic elastoplastic material, followed by an anisotropic description through crystal plasticity. Finally, the conclusion and future work will be presented in section 5.

Phase field governing equation for corrosion process
The corrosion system to be considered is composed of a metallic solid phase and corrosive environment. The breakdown of the passive film can lead to the dissolution of the metal solid. A simple electrode reaction is considered as follows, The metal M atom produces the ions M z+ in the solution and electrons. This electrochemical reaction is driven by the overpotential, i.e. ∆G = −zFη, where F is the Faraday constant and the overpotential η = E − E eq [30]. Here E is the potential of the system at which the reaction occurs, and E eq is the equilibrium potential. It is apparent that when E = E eq the forward and backward reactions have equal rates and thus no reaction occurs overall. At η > 0, ∆G < 0, the metal oxidation/dissolution (forward reaction) is favored. At η < 0, ∆G > 0, the ion reduction (backward reaction) is favored.
In the phase field framework, a phase field variable ξ is introduced to distinguish the metal (ξ = 1) and environment (ξ = 0), and varies continuously from 1 to 0 in the interfacial region. The total free energy consists of the chemical bulk and gradient terms: where the gradient term is given by f grad = κ 2 (∇ξ ) 2 , and κ is the gradient energy coefficient. The bulk term represents the local free energy density of the mixture of metal and solution, given by f bulk = h (ξ ) ∆G + g (ξ ) that. Here the interpolation function h (ξ ) = ξ 3 6ξ 2 − 15ξ + 10 interpolates the free energy at the interface and the double-well free energy function g (ξ ) = Wξ 2 (ξ − 1) 2 represents the energy barrier at the interface which are commonly used in the phase field-based framework [14,15].
Following Liang's framework [31] and the usual kinetic rate theory for chemical reaction, a nonlinear dependence of the interface kinetics on the electrochemical driving force (∆G) is assumed. The dependence on the interfacial energy is still considered as the Allen-Cahn form for non-conserved variable. Thus, the governing equation for ξ is as follows: where L int is the interface mobility, F int,total is the total interfacial energy of an inhomogeneous system given by F int,total =´V g (ξ ) + κ 2 (∇ξ ) 2 dV, and α, β are the charge transfer coefficient that satisfy α + β = 1. This equation is considered as the diffuse-type of Butler-Volmer equation.
The analytical solution of equation (3) under the sharp-interface limit of the phase field model can be obtained by multiplying the equation by ∂ξ /∂u and integrate over −∞ < u < ∞, where u is a coordinate normal to the interface. This will lead to the expression below, which is a similar form as the Butler-Volmer equation for the electrode reaction kinetic, following a linear dependence at a low positive overpotential and nonlinear dependence at a high positive overpotential.

Mechanical constitutive laws
Two different model materials in terms of elastoplastic responses will be considered: isotropic and anisotropic. The isotropic material is simulated as ideal elastoplastic with linear work hardening. The Young's modulus E, Poisson's ratio ν, yield strength σ Y and hardening modulus H are the model parameters of the mechanical behaviors. The mechanical strain ϵ mech is assumed to be the sum of elastic ϵ e and plastic strain ϵ pl as, A Radial Return approach is employed [32], where a trial stress is calculated by assuming the new strain increment is elastic strain ∆ϵ assumed,e via equation (6) and check if the trial stress σ (t+1) trial is outside of the von Mises yield surface. If the stress is outside, the amount of effective plastic strain ∆dϵ p,eff required to return the stress state to the yield surface will be computed via equation (7).
where C iso is the stiffness tensor for isotropic material, σ trial, eff is the scalar von Mises trial stress, and G is the shear modulus. For more details of radial return stress update of the J 2 plasticity model, interested readers are referred to [33,34]. The anisotropic material considers a linear elastic stress strain relation through the stiffness tensor C aniso . The plastic behavior is simulated using crystal plasticity material, where deformation gradient F is multiplicatively decomposed in its elastic and plastic parts as F = F e F p , and the second Piola-Kirchhoff tensor S is related to its conjugate Lagrangian strain tensor E as S = C aniso E e , with E e as the elastic part of Lagrangian strain tensor given by E e = 1 2 F eT F e − I . S can be computed from the Cauchy stress σ via equation (8), and the resolved shear stress τ α on each slip system α with associated slip direction s α and slip plane normal unit vectors m α can be calculated as given by equation (9), The power law slip rateγ α and the resistance to slip (i.e. slip system strength) g α are respectively given as,γ whereγ 0 is a reference slip rate, m is the strain rate sensitivity exponent, q αβ is a hardening coefficient matrix that accounts for the different self and latent hardening, h 0 is an initial hardening term, g 0 and g sat are the initial and saturated slip system strength, and a is the hardening exponent. The sum of the slip increments on all of the slip systems is given by For more details of the crystal plasticity model, interested readers are referred to [35].

Coupling between mechanical field and corrosion phase field
The effect of corrosion caused materials degradation on the mechanical field is considered by reformulating the moment balance equation ∇ · σ = 0 as to take into account the degraded solid, where a small positive number ζ = 1.0 · 10 −6 is introduced to ensure the numerical stability in the degraded area. The effect of mechanical field on the corrosion degradation is simulated using the film rupture-dissolution model adopted in this study [36][37][38], where stress/strain is to rupture the protective passive film at the crack tip and expose the susceptible bare metal underneath. The interface velocity will be increased due to the dissolution of metal upon the film rupture. As the metal can corrode in the absence of stress effect, the L int0 in equation (3) will be modified from a constant to incorporate mechanical straining in the form of equations (14) and (15) for isotropic and anisotropic model, respectively.
where ϵ p,eff is the effective plastic strain and ϵ f is the effective plastic strain required for film rupture, γ acc and γ acc,f are the accumulative slip and that required to rupture the film. equations (14) and (15) considers an increased interfacial mobility due to plastic straining from L int,passive that corresponds to the state with the full protection under passive film. On the other hand, the interfacial mobility cannot exceed that of the active dissolution, L dis , when the passive film does not exist.

Numerical implementation
The model parameter W and κ are related to the surface energy σ and interface thickness l, given by σ = √ κW 4 and l = α * 2 κ W , where α * = 2.94 corresponding to the definition of interface region 0.05 < ξ < 0.95 [39]. Introducing the characteristic length scale l c = 1.6l and time scale t c , the governing equation (3) can be rearranged into the form as given below, where the reduced parameters with respect to l c are given as∇ The mechanical field governing equation can be similarly reformatted employing the reduced length scale, The phase field and mechanical field governing equations are solved using the open-source MOOSE (Multiphysics Object-Oriented Simulation Environment) framework [29]. Adaptive mesh scheme is adopted so that the evolving interface will be captured with refined mesh whereas the domain without significant field variable variation maintains coarser mesh. An adaptive time scheme is also adopted to improve the computational efficacy.

Validation and verification studies
In this section, three examples are demonstrated to verify and demonstrate the capacities of the multiphysics model: (1) pitting in the pencil electrode, and crack initiation from an initial pit in (2) isotropic elastoplastic material and (3) anisotropic elastoplastic material.

Pencil electrode test
The pencil electrode test that represents the pure corrosion condition with no applied mechanical load is discussed first. The metal wire of diameter of 32 µm and of length of 128 µm is mounted in an epoxy coating with only one end exposed to the solution. The initial metal and solution domains are differentiated by the phase field variable ξ = 1 and ξ = 0, respectively (figure 1). A diffusive interface where 0 < ξ < 1 locates between the pure metal and solution domains. The boundary conditions are also illustrated in figure 1, where Dirichlet conditions ξ = 0 and ξ = 1 are applied at the solution boundary and at the left end of the metal wire, respectively. The Neumann condition ∂ξ /∂t = 0 is applied along the length of the metal wire insulated by the epoxy coating. The model parameters for the phase field evolution governing equation (15) are summarized in the table below following previous studies on a pencil electrode made of 304 SS in 1 M NaCl solution at an overpotential of 100 mV SCE [40]. The characteristic time scale is chosen as t c = 1 s. At the applied positive overpotential 100 mV, the uniform corrosion occurs to the metal wire. The evolution of the phase field ξ at different time instants is demonstrated in figure 2, where the refined mesh at the interface is shown in figure 2(e). The interface remains flat as  it advances into the metal, thus the variation of ξ extracted along the center line is presented in figure 2(f), showing that the diffusive nature of the interface captured with approximately ten elements at the interface. The advancing corrosion front with time is computed by finding the location where ξ is just below 0.5. Then the front velocity v n is estimated through linear fitting of the time-depth data. The current density can thus be calculated using the Faraday's second law: i n = zFv n c solid . A good consistency of both the corrosion depth and current density is demonstrated for the study of activation-controlled pit growth (figure 3(a)); thus, our model can exactly replicate the simulation results from Scheiner and Hellmich's study using a finite volume implementation [40]. The values of the phase field parameters presented in table 1 will be used for all the examples discussed in the following sections unless stated otherwise.  To further examine the activation-controlled interface mobility, a series of overpotential η varying from 10 to 100 mV is applied. A linear dependence of the interface movement on time indicating a constant velocity is identified for all conditions. The extracted interface velocity through linear fitting is shown in figure 3(b), where a transition from linear to exponential overpotential dependence is identified. It should be noted that the pencil electrode test is also discussed by Mai [18] and Cui [23], using the KKS model where the corrosion process is activation controlled at a low overpotential and diffusion controlled at a high overpotential. The model in our study is focused on the activation-controlled corrosion process so that the corrosion depth at a high overpotential still linearly vary with time.
In the following, the model will be demonstrated through two different examples including (1) initiation of crack from pit in isotropic materials in the two-dimensional (2D) and threedimensional (3D) conditions, and (2) initiation of crack from pit in anisotropic materials in the 2D and 3D conditions.

Crack initiation from pit in isotropic material
This example aims to understand the crack initiation from a pre-existing semicircle pit under applied mechanical loading. A model material of Young's modulus E = 200 GPa, Poisson's ratio v = 0.28, initial yield strength σ Y = 400 MPa and hardening modulus H = 1 GPa is employed. The constitutive model is first examined using a single-element 3D model and the von Mises stress-effective plastic strain relation is shown in figure 4. Then, a 2D plane strain case is considered for the combined electrochemical-mechanical phenomena, with the out-ofthe plane direction as z. The block dimension in the x − y plane is of width 64 µm and height 32 µm, respectively (figure 5). The initial pit is a semicircle of radius 1.6 µm, representing the pit size an early stage [41][42][43] (figure 5). The Neumann condition ∂ξ /∂t = 0 is applied at the top, left and right boundaries and Dirichlet condition ξ = 1 is employed at the bottom boundary. The left boundary is constrained with its x-direction displacement u x as zero and the bottom boundary is constrained with its y-direction displacement u y as zero. The top surface is kept stress-free and the displacement is applied to the right boundary at a strain ratė ε x = 2.5 × 10 −8 s −1 . The passive interface mobility L int,passive is taken as 1·10 4 times smaller than that at the dissolution stage, corresponding to a current density of 2.8·10 −5 A cm −2 . The characteristic time scale is chosen as t c = 1·10 4 s for this example and the other ones discussed in the following.
To clarify the effect of mechanical loading in enhancing localized corrosion, pit growth without applied displacement is first studied and the evolution of phase field variable is demonstrated in figure 6(a). The pit maintains a semicircle shape and has its radius increase linearly with time. With the applied displacement, enhanced corrosion is observed at the bottom of the pit where the plastic deformation is the highest and a transition from pit to crack initiation is observed (figures 6(b)-(d)). The effective plastic strain required for film rupture strain ϵ f is set as 0.05%, 0.1% and 0.2%, respectively. The different plastic strain values correspond to the     that close to the pit/crack tip. The pit/crack length advances with time have been plotted in figure 9, showing that the mechanical loading greatly enhances the localized corrosion rate and promotes the pit-to-crack transition. A smaller film rupture strain ϵ f indicates a lower threshold for crack to initiate from the ruptured passive film. Thus, the crack has been observed to initiate earlier and grows faster at a smaller ϵ f as expected ( figure 9). The effect of applied mechanical loading has recently been investigated in the SCC of Q345R in hydrofluoric acid by Dai et al [44]. In their study, it has been found that crack propagation rate increases with the applied stress increases in the elastic regime. In addition, the increase in stress promotes the transition of pits to cracks [44].
We proceed to understand the pit-to-crack growth in a plate of finite thickness in the 3D condition. The initial setup is shown in figure 10 with the semispherical shaped pit in the center of the same radius as the 2D counterpart. The block dimension is of width 64 µm, height 32 µm and depth 32 µm, respectively. The left (normal x), bottom (normal y), and back (normal z) surfaces are constrained in their displacement by setting u x , u y , and u z to be zero, respectively. The Neumann condition ∂ξ /∂t = 0 is applied at the all boundaries except for the bottom surface where the Dirichlet condition ξ = 1 is employed. The displacement is applied to the right boundary (normal x) at a strain rateε x = 2.5 × 10 −8 s −1 . Similar to the 2D case, the crack starts to initiate after the loading is significant enough to cause plastic deformation exceeding the film rupture strain. Noticeable crack initiation can be observed on the y-z plane perpendicular to the loading direction after 24 500 s (figure 11). Compared with the 2D plane strain condition, the crack initiation is slower at finite thickness due to lower plastic strain resulted from the reduced displacement constraint in the z direction ( figure 12).  Demonstration of the pit-to-crack from an initial semicircle pit at in the isotropic elastoplastic material; von Mises stress field is overlapped with the contour of the phase field variable at ξ = 0.5.

Crack initiation from pit in anisotropic material
In this last example, the material with anisotropic mechanical property is investigated. The initial and boundary conditions for phase field and mechanical field are the same as described above for the isotropic material model in the 3D condition. The 316L stainless steel is chosen as the model material and its crystal plasticity parameters [45,46]   An observation of the pit-to-crack growth in the crystal shown in figure 14 demonstrates similar crack initiation as the case for isotropic materials. To more clearly examine the dependence of the initial crack feature on crystal orientation, the stress and pit/crack contour are projected onto different coordinate planes. In the case of {100} y-plane model, the crack front starts to bifurcate into two branches as crack grows and the crack front also moves slower due to the smaller stress concentration ( figure 14(a)), whereas in the {110} y-plane model, pit/crack depth evolution is comparably faster and continue to propagate on the y-z plane without bifurcation ( figure 14(b)). It is possible that the initial orientation of the {100}-y plane crystal favors the simultaneous operation of multiple slip systems at the bottom of pit. An estimate based on the Schmid's law gives that eight out of the 12 slip systems will be in favor, which doubles the case for other two crystal orientations. This estimate is based on the assumption that local stress state remains as a uniaxial condition along the x-direction as the global condition. Due to complex interface morphology and thus the local stress state, the estimate above may not be accurate, yet still indicate that more slip systems may act simultaneously. This may also be attributed to the fact that the mechanical properties are symmetric respect to the plane perpendicular to the loading direction in the {100}-y plane crystal. The locations that surround the pit bottom thus have higher effective plastic strain instead of the pit bottom, which leads to bifurcation. In the {111} y-plane model, the crack tip growth is also relatively fast but no longer on the y-z plane that is perpendicular to the loading direction as the crack grows ( figure 14(c)). This could be ascribed to the anisotropic elastoplastic properties, the stress and accumulative slip distribution vary along the different crystal direction, and could deviate from the perpendicular to applied stress direction. Note here that the different corrosion propensities of the crystallographic planes are not considered here.
The simulations are then conducted in the 2D case under plane strain condition and using the same model setup as descried above for the isotropic case (i.e. figure 5), the key features found in the 3D condition (bifurcation, propagation rate difference, crack front deviation from perpendicular plane) are still captured in the 2D case (figures 15 and 16). It has been also observed experimentally that crack could travel at an angle from the perpendicular or even parallel to the applied load direction in SCC, and crack growth direction is reported to closely related to the crystallographic texture existent in the material [47,48].
Finally, the effects of corrosion related degradation on the mechanical responses of the material are analyzed in terms of the stress-strain relation along the loading direction in the 2D situation ( figure 17). The 'nodegrade' condition is conducted without corrosion effect considered. It can be observed that the deduction of stress due to corrosion starts to be prominent even though most of the domain is still in the elastic domain. The initiation and degree of stress reduction also strongly depend on the crystal orientation: a 40% deduction in stress is found at ε x = 0.12% in the {100} y-plane model, whereas a tremendous reduction (∼60%) is already observed at ε x = 0.06% in {110} and {111} y-plane model.

Conclusions and future work
In this paper, we present a new multiphysics phase field model to simulate the crack initiation from corrosion pits in isotropic and anisotropic materials. A modification to the Allen Cahn equation is made to incorporate the nonlinear dependence on the thermodynamic driving force of the interface kinetics, so that the Butler-Volmer kinetics of the electrode-electrolyte interface migration can be captured. The mechanical effect on the corrosion kinetics is considered through the film dissolution mechanism, by enhancing the interface kinetics when the mechanical strain is sufficient to rupture the surface film. The anisotropic mechanical behavior is described by the crystal plasticity model. The applicability of this computational framework to understand the effects of materials parameters and model configuration is also demonstrated in 2D and 3D examples. In all conditions, a pit-to-crack transition has been observed. It has also been shown that the crack initiation could strongly depend on crystal orientation and deviate from the perpendicular to the applied loading direction in materials with anisotropic mechanical behaviors.
There are a number of model parameters that can be varied including the electrochemical conditions, mechanical properties, initial pit size and shape, and loading condition. As an introduction to the model, we tested limited model parameters in this manuscript to understand the pit-to-crack transition. In future, we will also investigate the model through more benchmark problems such as those involving more complex structure geometry and loading conditions. Furthermore, structures with inhomogeneous or gradient materials properties are also of great interest. In addition, we will focus on further improving the model capacities. The current model does not include the diffusion-controlled interface kinetics, and this will be improved by introducing the field variables related to chemical species and modifying the governing equations of corrosion processes. Last but not least, this multiphysics model capable to simulating coupled electrochemical and mechanical phenomena can be modified and extended to similar processes such as that occurred in solid fuel cell battery [49][50][51][52][53].

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://clmmd.wordpress.com/ [54]. Data will be available from 28 April 2023.