Phase-field modeling of pitting and mechanically-assisted corrosion of Mg alloys for biomedical applications

A phase-field model is developed to simulate the corrosion of Mg alloys in body fluids. The model incorporates both Mg dissolution and the transport of Mg ions in solution, naturally predicting the transition from activation-controlled to diffusion-controlled bio-corrosion. In addition to uniform corrosion, the presented framework captures pitting corrosion and accounts for the synergistic effect of aggressive environments and mechanical loading in accelerating corrosion kinetics. The model applies to arbitrary 2D and 3D geometries with no special treatment for the evolution of the corrosion front, which is described using a diffuse interface approach. Experiments are conducted to validate the model and a good agreement is attained against in vitro measurements on Mg wires. The potential of the model to capture mechano-chemical effects during corrosion is demonstrated in case studies considering Mg wires in tension and bioabsorbable coronary Mg stents subjected to mechanical loading. The proposed methodology can be used to assess the in vitro and in vivo service life of Mg-based biomedical devices and optimize the design taking into account the effect of mechanical deformation on the corrosion rate. The model has the potential to advocate further development of Mg alloys as a biodegradable implant material for biomedical applications.

Despite successful clinical studies, Mg-based implants have only been approved for a limited number of applications.Rapid corrosion in an aggressive chloride medium like human body fluids is the main reason limiting the widespread use of Mg alloys as a biodegradable material [7] .In vivo [8][9][10] and in vitro studies [11][12][13] have reported that biodegradable Mg alloys (as well as industrially used ones) generally corrode in a localized fashion (e.g., pitting), Fig. 1 .For instance, an in vivo study [14] has shown that Mg alloy plates and screws in the diaphyseal area of long bones in pigs were nearly completely degraded after 12 months of implantation.The screws showed a faster and nonhomogeneous degradation profile in the intramedullary cavity owing to enhanced exposure to interstitial fluids.Another study [15] has found that corrosion resistance is the major challenge for using Mg interference screws for anterior cruciate ligament reconstruction.Moreover, regions of the interference screw exposed to synovial fluids in the knee joint cavity suffer accelerated degradation.
Different strategies, mainly based on surface modification, have been developed to mitigate these risks and improve the corrosion resistance and biocompatibility of Mg alloys [16,17] .Yet, pitting corrosion could only be diminished to a certain extent [18,19] .Moreover, load-bearing biomedical devices are continuously subjected to various loading conditions during service.In such an environment, biodegradable Mg alloys show sensitivity to stress corrosion cracking (SCC) [19][20][21][22] .The synergistic effect of mechanical loading and a corrosive environment significantly reduces the corrosion resistance and mechanical integrity of Mg alloys.The concurrence of these two factors dramatically accelerates the corrosion rate and promotes crack propagation ( Fig. 1 ), leading to the sudden failure of implants [23] .A recent study [24] has reported that elastic strains decrease the corrosion potential, increase corrosion current and accelerate the degradation of WE43 Mg wires, while plastic strains enhance localized corrosion.Hence, SCC is a severe concern for thin implant applications like stents, membranes, and wires.
Before clinical trials, the corrosion performance of Mg alloys is usually determined in in vitro tests under various corrosive environments that resemble biological fluids.Ideally, this informa-tion can be used to calibrate numerical models to predict the in vivo performance of implants and to guide design taking into account their progressive degradation.A wide variety of computational tools have been developed to predict the corrosion behavior of biodegradable Mg alloys.Phenomenological models based on the continuum damage (CD) theory [25] for uniform and localized corrosion [26,27] are based on a scalar damage parameter that reduces the mechanical properties of the corroded regions of the material.The damage evolution law for stress corrosion depends on threshold stress above which corrosion progresses [26] while pitting corrosion is introduced via a dimensionless pitting parameter controlled by a distribution probability density function [27] .The CD approach is further improved by including several different features [28][29][30][31][32] .However, the diffusion process, as the underlying physical mechanism in the corrosion process, is not included.In more advanced phenomenological models [33,34] , the diffusion of Mg ions and the evolution of other species have been incorporated through physicochemical interactions and diffusion equations.
Although phenomenological models have been typically used to simulate Mg corrosion, there is a growing interest in the development of physically-based models that can resolve the physical processes governing corrosion and thus provide mechanistic predictions and insight [35][36][37][38][39][40][41] .While the underlying physics is relatively well understood, there are significant theoretical and computational challenges intrinsic to the coupled nature of the problem and the difficulty of tracking the evolution of complex corrosion interfaces in arbitrary domains.Regarding the former, a mechanistic model of Mg corrosion must account for Mg dissolution at the corrosion front (short-range interactions), diffusion of Mg ions in solution (long-range interactions) and the coupling with other physical phenomena such as capturing the role of mechanics in enhancing corrosion rates and the electrochemistry-corrosion interplay.Regarding the challenges associated with tracking an evolving corrosion front computationally, a number of numerical techniques have been recently proposed, including Arbitrary Lagrangian-Eulerian (ALE) approaches [35][36][37] , level set methods [38][39][40] , and peridynamics [41] .However, these are mainly used in the context of uniform corrosion as they are still limited in capturing localized corrosion, coupling with other physicomechanical phenomena, and handling geometric interactions in arbitrary di- mensions (2D/3D), such as the coalescence of corrosion pits.Phasefield formulations have emerged as a promising approach for modeling moving interfaces and handling topological changes at different length scales [42] .In phase-field models, the interface between two phases is smoothed over a thin diffuse region using a continuous auxiliary field variable (e.g., φ), see Fig. 2 .The phasefield variable φ has a distinct value in each phase (e.g., 0 and 1), while varying smoothly in between.The movement of the interface is implicitly tracked without presumptions or prescribing the interface velocity.Topological changes of arbitrary complexity (e.g., divisions or merging of interfaces) can be naturally captured in 2D and 3D without requiring any special treatments or ad hoc criteria.The phase-field method has been recently extended to several challenging interfacial phenomena relevant to corrosion in nonbiodegradable metallic materials [43][44][45][46][47][48][49][50] and to internal galvanic corrosion and damage localization induced by insoluble secondary phases in Mg alloys [51] .This work aims to extend this success to biomaterial degradation, presenting the first phase-field model for the surface-based localized (pitting) corrosion of biodegradable Mg alloys that incorporates material dissolution, Mg ionic transport, and mechano-chemical interactions.
The outline of the paper is as follows.The degradation mechanisms governing mechanically-assisted corrosion of biodegradable Mg alloys are presented in the following section and the phasefield model is subsequently formulated.The interplay between Mg dissolution, ionic transport in solution, and mechanical straining is captured by defining a generalized thermodynamic free energy functional that incorporates chemical, interfacial, and mechanical terms.The impact of mechanical fields in accelerating corrosion kinetics is integrated through a mechano-electrochemical mobility term that depends on local stress and strain distributions.The constructed model is calibrated and validated against in vitro corrosion data on WE43 Mg alloy wires immersed in simulated body fluid in Section 3 .Dedicated experiments are conducted to validate the model, both qualitatively (pitting patterns) and quantitatively (hydrogen gas released), demonstrating as well its ability to capture localized corrosion phenomena.After validation, the potential of the model to handle mechano-chemical effects during corrosion is demonstrated in Section 4 through two representative case studies: pitting corrosion associated with the local failure of a protective layer and the nonhomogeneous stress state of a bioabsorbable coronary stent.The potential of the model is discussed in Section 5 along with recommendations for future work.Conclusions of the investigation are summarized in Section 6 .

Degradation mechanisms
Magnesium dissolution in aqueous environments, such as biological fluids, is governed by an electrochemical reaction and the corrosion process can be summarized as follows [3] Mg (s ) → Mg 2+ (aq ) + 2 e − (anodic reaction) ( The last reaction in Eq. ( 1) is the precipitation reaction that leads to the formation of a passive layer of magnesium hydroxide (Mg(OH) 2 ) on the Mg surface, Fig.
undermining the integrity of the passive film.Fast corrosion rates and pitting corrosion are generally associated with aggressive chloride ions [52] .The presence of inorganic ions and organic compounds in body fluids further increases the complexity of the degradation process [53] .While the effect of certain organic compounds [54] and inorganic ions [55,56] on the corrosion rate have been identified, the degradation mechanisms of Mg alloys in body fluids are not fully understood [57,58] .Therefore, it is assumed that the primary degradation mechanism is driven by the bulk diffusion of Mg ions in the physiological environment.The complex composition of the porous protective layer, its negligible thickness compared to the size of the surrounding body fluid, and the high solubility of MgCl 2 in aqueous environments allow to neglect the product formation and layer dissolution reactions, an approach frequently followed in the literature [26,27,32,34,35,41] .The presence of mechanical stresses increases the corrosion susceptibility of Mg alloys.In vitro studies [24,59] have indicated that mechanical fields decrease the corrosion potential of Mg alloys, thereby increasing corrosion current densities and dissolution rates.Following Gutman's theory of mechano-electrochemical interactions [60] , the anodic dissolution kinetics is given as where i is the anodic dissolution current in the presence of mechanical stresses, i 0 the anodic dissolution current in the absence of mechanical stresses, ε p the effective plastic strain, ε y the initial yield strain, σ h the hydrostatic stress, V m the molar volume of the metal, R the universal gas constant, and T the absolute temperature.After rupture of the passive film and pit nucleation, local stress and plastic strain distributions intensify local material dissolution in the vicinity of the pit, promoting pit-to-crack transition and crack propagation, as schematically illustrated in Fig. 1 .
As shown below ( Section 2.4 ), the amplification factor in Eq. ( 3) is embedded into the model kinetics parameter that characterizes solid-liquid interface movement to incorporate the role of mechanical fields in accelerating corrosion.

Thermodynamics
The problem formulation is depicted in Fig. 2 and could be summarized as follows.The system consists of a biodegradable Mg alloy in contact with physiological environments that, by composition, mimic body fluids.The system domain includes both the Mg alloy and the corrosive environment.A continuous phase-field parameter φ is introduced to distinguish different phases: φ = 1 represents the solid phase (Mg alloy), φ = 0 corresponds to the liquid phase (physiological fluid), and 0 < φ < 1 indicates the thin interfacial region between the phases (solid-liquid interface).With vanishing normal fluxes ( n • J = 0) on the domain boundary ∂ , the independent kinematic variables necessary for model description are the non-conserved phase-field parameter describing the evolution of the corroding interface φ(x , t ) , the displacement vector to characterize deformation of the solid phase u (x , t ) , and the normalized concentration of Mg ions c Mg (x , t ) with respect to the concentration in the solid phase ( c Mg = c Mg /c s Mg ).More details regarding nondimensionalization are given in Section 2.5 .
The free energy functional for a heterogeneous system such as the one in Fig. 2 can be written as where f chem , f grad , and f mech are the chemical, gradient, and mechanical energy densities defined below.

Chemical free energy density
Following the phase-field model for phase transitions in binary alloys [61] , the chemical free energy density of a homogeneous system consisting of solid and liquid phases is decomposed into the chemical energy density associated with material composition and double-well potential energy where f chem l ( c l Mg ) and f chem s ( c s Mg ) are the chemical free energy densities within the liquid and solid phases as a function of normalized phase-concentrations c l Mg and c s Mg .In the above equation, g(φ) and h (φ) are the double-well potential energy and interpolation functions commonly expressed as ω in Eq. ( 5) is a constant that determines the energy barrier height at φ = 1 / 2 between the two minima at φ = 0 and φ = 1 .
The chemical free energy densities within each phase in Eq. ( 5) are approximated by simple parabolic functions with the same curvature parameter A as where c l,eq Mg = c l,eq Mg /c s Mg and c s,eq Mg = c s,eq Mg /c s Mg are the normalized equilibrium Mg concentrations in the liquid and solid phases (refer to Section 2.5 for dimensional analysis).Alternatively, the chemical free energy density can be approximated assuming a dilute solution [46][47][48] .Physically, the equilibrium concentration in the solid phase c s,eq Mg represents the average concentration of Mg ions within the material.Since the product formation and protective layer dissolution are neglected in the current work ( Section 2.1 ), c l,eq Mg is determined based on the mass density and molar mass of MgCl 2 formed on the exposed Mg surface.
The interfacial region is defined as a mixture of both phases with different concentrations but with the same diffusion chemical potential [61] c Using Eqs. ( 7) and ( 8) renders the following definition for the chemical free energy density of the system

Gradient energy density
The interfacial energy density is defined as where κ is the isotropic gradient energy coefficient.The phasefield parameters ω and κ are connected to the physical quantity (interfacial energy ) and computational parameter (interface thickness ).For the accepted double well potential g(φ) in Eq. ( 6) , the following relations are obtained [62]

Strain energy density
The mechanical behavior of the solid phase is assumed to follow the von Mises theory of plasticity [63] .Considering deformable elasto-plastic solids, the mechanical free energy density f mech in Eq. ( 4) is additively decomposed into elastic f mech e and plastic components (12) where h (φ) ensures the transition from the intact solid (uncorroded Mg alloy) to the completely corroded (liquid) phase.The elastic strain energy density f mech e is a quadratic form of the elastic strain where C is the rank-four elastic stiffness tensor and ε e is the elastic strain tensor obtained by subtracting the plastic strain tensor ε p from the total strain ε .For linearized kinematics, the total strain tensor is the symmetric part of the displacement gradient The elastic deformation of the solid is described by the isotropic linear elasticity theory so that the rank-four elastic stiffness tensor reads where λ and μ are the Lamé elastic constants.
The plastic strain energy density f mech p is incrementally computed from the plastic strain tensor ε p and the Cauchy stress tensor σ 0 for the intact configuration as

Governing equations
Using the balance of power and the principle of virtual power [64] , the following time-dependent governing equations for the independent kinematic fields φ(x , t ) , c Mg (x , t ) , and u (x , t ) are derived (see details in the Supplementary Materials) complemented with boundary conditions The resulting set of governing equations includes the Allen-Cahn equation [65] for the non-conserved phase-field parameter, the diffusion equation for the Mg concentration in the liquid and solid phases, and the linear momentum balance equation for quasistatic mechanical deformation.In the above equation, L is the kinetic coefficient that characterizes the interfacial mobility and D c Mg the effective diffusion coefficient interpolated with the phase-field parameter between the phases where D l c Mg and D s c Mg stand for the diffusion coefficients of Mg ions in the liquid (corrosive environment) and solid phases.The role of mechanical fields on the interface kinetics is incorporated by modifying the interface mobility parameter L , which includes a mechano-electrochemical contribution that amplifies the dissolution process, as shown in Section 2.4 .Thus, the mechani- 17) ).For an alternative way of incorporating the me-chanical contribution to interface kinetics, the interested reader is referred to Refs.[46][47][48] .

Mechano-electrochemical coupling
The role of mechanical fields in enhancing corrosion kinetics is incorporated by following Gutman's theory [60] .As shown in Eq. ( 3) , the anodic dissolution can be amplified by an amplification factor that depends on local stress and strain distributions.As the anodic dissolution kinetics dictates interface motion, the interfacial mobility coefficient L is analogously connected to mechanical fields.Using Eq. ( 3) and considering the linear relationship between L and i a (corrosion current density) [45] returns the following expression for the kinetic coefficient in Eq. ( 17) where L 0 is the interfacial mobility that physically corresponds to the anodic dissolution current i 0 in the absence of mechanical stresses and plastic strains, Eq. ( 3) .The interfacial mobility L 0 is determined in Section 3 considering stress-free corrosion experiments on Mg wires.

Dimensional analysis
To facilitate numerical simulations and improve convergence, the governing equations Eq. ( 17) are normalized using the interface thickness as the characteristic length, Mg concentration in the solid phase c s Mg , diffusion coefficients of Mg ions in the liquid phase D l c Mg , and the energy barrier height ω as the energy normalization factor.Thus, the nondimensional time t , nondimensional space coordinates x , and nondimensional gradient ∇ are given as The above nondimensional variables return the following governing equations along with the corresponding nondimensional boundary conditions.The details of the numerical implementation of Eq. ( 23) are given in Supplementary Materials.
The characteristic times for diffusion t d and interface reaction t φ are then given by and their ratio determines the rate-limiting process.For the case of τ 1 (i.e., t d t φ ), diffusion is slower than interface reactions and the process is driven by bulk diffusion.This situation is denominated diffusion-controlled corrosion.On the contrary, diffusion is faster if τ 1 (i.e., t d t φ ) so that there is no accumulation of Mg ions at the metal-fluid interface.Under that condition, the rate of material transport is interface reaction-controlled (commonly called activation-controlled corrosion).The criterion for the interfacial mobility coefficient for the two rate-limiting processes reads The effects of diffusion-and activation-controlled processes on the corrosion behavior are discussed in Section 5 .

Materials and experimental methods
An in vitro degradation study was carried out on WE43MEO Mg alloy wires of 0.3 mm in diameter to validate the proposed phasefield model.The Mg wires were manufactured by cold drawing at Meotec GmbH (Aachen, Germany) from WE43MEO Mg alloy with a nominal composition of 1.4-4.2%Y, 2.5-3.5% Nd, < 1% (Al, Fe, Cu, Ni, Mn, Zn, Zr) and balance Mg (in wt.%).The wires were annealed at 450 °C for 5 s after cold drawing to reduce the dislocation density induced during drawing and improve the ductility [66,67] .
Corrosion tests were carried out in wires of 120 mm in length immersed in Simulated Body Fluid (c-SBF) at 37 °C.The experimental setup is schematically depicted in Fig. 3 .The composition (per liter) of the c-SBF was 8.035 g NaCl, 0.355 g NaHCO 3 , 0.225 g KCl, 0.176 g K 2 HPO 4 , 0.145 g MgCl 2 , 0.292 g CaCl 2 , 0.072 g Na 2 SO 4 , and 50 mL of Tris buffer pH 7.5.The ratio of c-SBF volume to the wire surface area was > 0.5 mL/mm 2 , according to the ASTM G31-72 standard.
The degradation rate was assessed by measuring the amount of hydrogen gas released which is, according to Eq. ( 1) , equivalent to the mass loss of corroded Mg.To measure the evolved hydro-gen gas, the Mg wires were placed inside a glass burette, as illustrated in Fig. 3 .The burette was inserted into a sealed plastic bottle filled with SBF.The released hydrogen gas was captured in the burette and tracked by a eudiometer.Twelve samples were used in the immersion tests.After 24 h of immersion, four samples were taken out of the SBF to assess the extent of pitting corrosion.To this end, twenty random cross-sections of the corroded Mg wires were mounted in an epoxy resin, grounded, polished, and images were taken in an optical microscope.These images were analyzed using the PitScan Framework [68] to quantify the degree of pitting corrosion.

Statistical analysis
The experimental measurements and numerically obtained results in Section 3.2 are expressed as mean value ± standard deviation.Microsoft Excel software was used for the statistical calculations.

Model validation
Two types of simulations are performed to validate the phasefield model.First, experimental measurements of hydrogen release over the immersion time are used to calibrate the model kinematic parameter L 0 considering uniform corrosion.Second, the wire cross-sections measured after 24 h of immersion in SBF are compared with pitting corrosion predictions.In these simulations, the mechanical effect is not considered.The properties of the Mg alloy used in the simulations are listed in Table 1 .Due to the lack of experimental data on the diffusivity of Mg ions in biological fluids, the magnitude of D l c Mg is estimated as the average value utilized in various numerical studies [33][34][35][36][37][38][39][40] .Although the concentration of Mg ions in the solid phase ( c s Mg ) could be evaluated using the material data for pure Mg and the mass fraction of alloying elements and impurities, the value for pure Mg is used in this investigation for simplicity [69] .The physiological environment and the presence of chloride ions determine the equilibrium concentration of Mg ions in the liquid phase c l,eq Mg (saturated concentration).As the formation of the partly protective layer is neglected ( Section 2.1 ), the saturated concentration of Mg ions in the cor-

Table 1
Parameters common to all phase-field simulations.The phase-field parameters, energy gradient coefficient κ and energy barrier height ω, are connected to the interfacial energy and the interface thickness , Eq. (11) .While the interface thickness is a purely computational parameter whose choice is based on the scale of the problem, the interfacial energy is a physical quantity that depends on crystallographic orientation.The average value reported for pure Mg in Refs.[70,71] is used in this investigation.Lastly, the chemical free energy density curvature parameter A in Eq. ( 7) is assumed to have a similar value as in Refs.[43,45] for corrosion in metallic materials.

Phase-field simulations of uniform corrosion
The phase-field simulations of uniform corrosion are performed using an axisymmetric domain as illustrated in Fig. 3 .The nondimensional form of governing Eq. ( 23) is solved with accompanying initial and boundary conditions.A smooth equilibrium phase-field profile is prescribed as the initial solid-liquid interface.The interface thickness is selected to be significantly smaller than the diameter of the Mg wire ( = 4 μm).No flux boundary conditions for diffusion and phase-field are imposed at all the outer edges of the domain to simulate an unbounded environment.These boundary conditions preserve mass conservation and imply that no diffusion occurs across the domain boundary.The applied boundary conditions and the large domain size of SBF in the horizontal direction are selected to mimic the experimental setup and ensure that the solution does not saturate with Mg ions.The solid material and the surroundings are isotropic.Thus, the vertical dimension of the domain does not influence the results, and the analysis can be reduced to a one-dimensional axisymmetric problem.
The volume of hydrogen released per unit of exposed area is calculated in the simulations from the ideal gas law where P is the pressure (1 atm), A the exposed area, and n Mg the total amount of dissolved Mg in (mol) determined as The predicted hydrogen gas evolution per unit area of the Mg wire is plotted as a function of the immersion time in Fig. 4 , together with the experimental data obtained from the corrosion tests in c-SBF.The experimental results show that the corrosion rate was initially fast and approximately linear up to 24 h.The corrosion rate slowed down afterward to reach a plateau at 120 h.The phasefield simulations return the same trends and accurately reproduce The good agreement between experiments and simulations indicates that the first factor is dominant while the effect of the protective layer (that is not considered in the phase-field simulations) can be neglected in the presence of Cl − ions.

Phase-field simulations of pitting corrosion
A certain degree of randomness needs to be introduced in the system to simulate pitting corrosion.Even assuming uniform alloy composition and surface properties, pitting may occur due to nonuniform distributions of aggressive Cl − ions, as those ions undermine the protective layer ( Eq. ( 2) ).Following that analogy, pitting is introduced in the model through a spatially-dependent kinetic coefficient L 0 , correlating it to a random nonuniform distribution of Cl − ions.Areas with higher values of L 0 reflect the higher concentration of aggressive Cl − ions, thereby promoting pit-ting corrosion.Random distribution functions are used to define the nonuniform distribution of Cl − ions and to capture the stochastic nature of pitting.
Introducing randomness in terms of the nonuniform distribution of Cl − ions breaks the axial symmetry conditions and, consequently, requires 3D simulations.Performing such simulations for the geometry considered (very long and thin wires) is computationally expensive.Hence, without loss of generality, 3D simulations are replaced by multiple 2D simulations to provide statistical information about the pitting corrosion metrics that can be compared to the experimental data obtained from numerous crosssections examined in the optical microscope.Space-dependent 2D data is constructed using a sum of trigonometric functions combined with two uniform random distribution functions that introduce randomness in the system.The sum of trigonometric functions can be seen as the assembly of many spatial waves whose amplitudes and phase angles are defined through the two random distribution functions.Therefore, the spatially dependent interfacial mobility parameter L 0 can be written as where L 0 is the interfacial mobility parameter determined for uniform corrosion and f ( x , ȳ ) acts as a dimensionless pitting function.This function is represented by the double summation term (i.e., the sum of spacial waves in the x and y directions) and serves as a stochastic function to introduce randomness in the system.In the previous expression, x and ȳ are normalized spatial coordinates, m and n are spatial waves in the x and y directions.The number of spatial waves in both directions is N.The first uniform random distribution function γ (m, n ) is defined between zero and one and determines random amplitudes.High-frequency amplitudes are attenuated with the exponent β to generate smooth amplitude coefficients.Higher β values return a smooth (more uniform) pitting function.The second uniform random distribution function, which controls the phase angle of each wave ϕ(m, n ) , is defined between −απ / 2 and απ/ 2 and governs the spatial distribution of the data.Thus, different spatial distributions, periodicity, magnitudes, and smoothness of random data are controlled with the number of spatial frequencies N, exponent β, and the range of the distribution function ϕ(m, n ) varying the α value.For the purpose of pitting corrosion, the coefficients a 0 and a 1 are included to preserve the non-negativity of the interfacial mobility parameter ( L 0 > 0) and control the desired difference between the maximum and minimum amplitude to manage the pitting intensity.
Three pair sets of N and β are selected, such that the first pair is N = 2 and β = 0 . 1 , the second pair N = 2 .5 and β = 1 .5 , and the third pair N = 1 .25 and β = 0 .75 .For each N − β pair set, three different values of the α parameter are considered, i.e., α = 1 , α = 3 , and α = 5 .For each of these nine combinations, three different amplitudes are applied (by adjusting the coefficients a 0 and a 1 ) such that the spatially dependent interfacial mobility parameter lies in between 0 .very similar to the phase-field simulations but quantitative comparisons between experiments and simulations can be carried out through three different metrics parameters [68] .They are (i) the uniform corrosion radius in Fig. 6 (a) (the radius of the circular section that has the same area as the corroded cross-section), (ii) the average pit depth in Fig. 6 (b) (average distance from the degraded cross-section to the uniform corrosion circle), and (iii) the maximum pit depth in Fig. 6 (c) (maximum distance from the corroded cross-section to the uniform corrosion circle).The experimental uniform corrosion radius after 24 h of immersion in SBF was 108 ± 21 μm, which corresponds to approximately 48% mass loss and is in agreement with hydrogen gas evolution tests.The higher standard deviation indicates the variation of mass loss among different sections of wire.The experimental values of the average and maximum pit depth were 26 ± 12 μm and 56.62 ± 19.3 μm, which shows the severity of pitting corrosion at particular cross sections.
The average experimental values with standard deviations of these three parameters obtained from the analysis of ten different crosssections are shown in Fig. 6 (a)-(c), together with corresponding results obtained from the twenty-seven phase-field simulations.The agreement between experiments and simulations in terms of the uniform corrosion radius and the average pit depth is satisfactory, while the simulations slightly underestimate the maximum pit depth.Overall, the agreement between the experimental measurements and phase-field predictions for hydrogen gas evolution and pitting metrics indicates that the proposed model can be utilized to simulate uniform and pitting corrosion of biodegradable Mg alloys immersed in physiological environments.The model satisfactorily predicts hydrogen gas evolution and captures the experimental trend for pitting corrosion.In the following section, the model is used to ascertain the role of mechanical fields in accelerating the corrosion process.

Applications
The proposed framework to predict the degradation of Mg alloys in physiological environments is applied in this section to assess the evolution of corrosion in the presence of mechanical stresses in two different scenarios.The first deals with a wire loaded in tension in which the protective layer is damaged, leading to the formation of a pit.The second one analyzes the corrosion of a bioabsorbable Mg alloy coronary stent.

Stress-assisted corrosion of Mg wires
In this simulation, the Mg wire is simultaneously immersed in SBF and subjected to tensile deformation along the wire axis.It is further assumed that the wire surface is protected against corrosion by a thin surface layer locally damaged in a small area.The initial breakdown of the protective layer enables the ingress of aggressive Cl − ions leading to the nucleation of a pit that acts as a stress concentrator.The initial pit has a semi-circular shape with a radius of 10 μm around the whole diameter of the wire to maintain axisymmetric boundary conditions, Fig. 7 .
Due to symmetry, only half of the axisymmetric domain is considered in the simulation, as depicted in Fig. 7 .Similarly to the previous case, to represent an unbounded domain, no flux (Neumann) boundary conditions are enforced at all the outer boundaries of the computational domain for both the phase-field and the Mg concentration.The protective film is modeled as an impermeable layer with a thickness of 0.5 μm around the wire surface with the corresponding no flux boundary condition for both Mg concentration and phase-field.To investigate the influence of the mechanical fields on corrosion kinetics, additional constraints are enforced for the mechanical equilibrium equation.The normal component of the displacement vector along the vertical and horizontal symmetry axes is constrained ( n • ū = 0 ) while a non-zero remote tensile deformation ε ∞ is prescribed on the top surface, Fig. 7 .The remote deformation is prescribed at the beginning of the simulation and held fixed over a total simulation time.Its magnitude is varied to study the role of mechanical fields in enhancing corrosion kinetics and SCC behavior.
The material properties and the phase-field parameters used are the same as in the previous case study of uniform corrosion.
The interfacial mobility coefficient L incorporates the role of mechanical fields through Eq. ( 20) , whereas L 0 is previously determined in the comparison with load-free experiments.The mechanical properties of an AZ31 Mg alloy from the literature are used for the simulations [27] .The Mg alloy is assumed to behave as an isotropic, elasto-plastic solid.The Lamé elastic constants are λ = 38 GPa and μ = 16.3GPa.Plastic deformation is described using the J 2 flow theory with non-linear isotropic hardening, with a yield stress of 138 MPa and an ultimate tensile strength of 245 MPa at an engineering strain of 17%.
The obtained results in terms of phase-field contours, Mg concentration distribution, and mechanical fields for various remote deformations ε ∞ after 24 h of immersion in SBF are presented in Fig. 8 .In the absence of mechanical loading ( ε ∞ = 0 ), the pit grows uniformly, keeping the initial circular shape with a low and uni-   form concentration of Mg ions within the pit.The application of a relatively small axial deformation (e.g., ε ∞ = 0.096%) increases the magnitude of σ h in a small localized area and produces negligible plastic deformation.The hydrostatic stress distribution changes the pit morphology initiating a pit-to-crack transition.The Mg concentration increases near the tip of the pit, indicating that corrosion of Mg is localized in this region because of the stress concentration associated with the sharp tip.Further increase in the applied strain ( ε ∞ = 0.1%) raises the stresses high enough to trigger noticeable plastic deformations.The shape of the evolving defect is governed by both the hydrostatic stress and the plastic strain distribution.Longer and smoother cracks are observed compared to the previous case ( ε ∞ = 0.096%), as the hydrostatic stress and plastic strain distributions engage a more extensive area.The Mg concentration is significantly increased at the crack tip, but it is still well below the equilibrium value in the liquid phase ( c l,eq Mg = 0.57 mol/L), indicating an activation-controlled process.Model predictions in terms of pit depth and hydrogen gas evolution for the cases considered are given in Fig. 9 .Both pit kinetics and hydrogen gas production increase with an increase in applied strain.However, the pit kinetics is dramatically altered in the presence of mechanical loading, leading to rapid crack growth and fracture of the wires after a short time in SBF.

Bioabsorbable coronary Mg stent
The potential of the model is demonstrated in predicting the degradation of a bioabsorbable coronary Mg stent immersed in biological fluid.Mg alloy-based stents are attractive as temporary scaffolds to diseased blood vessels and exhibit good clinical performance [5] .However, premature failure due to fast corrosion rates limits their cardiovascular applications.A physically-based model for uniform corrosion [35] and phenomenological approaches [74][75][76] have been employed in simulating the degradation of Mg stents, considering various stent geometries and Mg alloys.The present phase-field corrosion model is applied in this section to simulate stent degradation taking into account both uniform and pitting corrosion.
An idealized stent geometry that resembles different stent designs frequently used in the literature [74][75][76] is considered in this work.The stent has an outer diameter of 2 mm and a total length of 7.5 mm.The whole geometry comprises six rings interconnected by a link with a length of 0.30 mm and a diameter of 0.125 mm.Each ring has six peak-to-valley struts with a total height of 1 mm and a diameter of 0.15 mm, Fig. 10 (a).
The computational domain consists of the stent and the surrounding corrosive environment.It is assumed that the stent is immersed in a physiologically representative blood vessel environment such that the initial Mg concentration corresponds to human blood (0.875 mmol/L).As in the previous examples, no-flux boundary conditions are imposed on all the external surfaces of the domain for phase-field and Mg concentration.The size of the corrosive environment is significantly larger than the stent geometry to avoid saturation effects.The material properties, the phase-field parameters, and the interfacial mobility parameter L 0 follow those used in the previous examples.
Two different case studies are considered.In the first study, the stent is mechanically loaded before immersion in the corrosive environment.It is radially expanded to an outer diameter of 2.25 mm, mimicking the balloon inflation stage during the deployment process.The balloon is modeled as a rigid cylindrical body.This step is followed by the stent recoil, which corresponds to the balloon deflation and extraction process.The final stent outer diameter following the recoil is 2.168 mm.The stent deployment process is summarized in Fig. 10 (b).The stress state in terms of von Mises stresses and equivalent plastic strains for the representative ring element after stent recoil is shown in Fig. 10 (c).The plastic strains are then incorporated into the subsequent corrosion simulation.In the second study, the as-manufactured stent is immersed in biological fluid in the absence of mechanical stresses.This case corresponds to uniform corrosion and is a reference study for comparison.
The results of the phase-field simulations for mechanically assisted and uniform corrosion are given in Fig. 11 .In the former case, plastic strains are localized at the union between rings and links (as shown in Fig. 10 (c)), providing hot spots for pitting nucleation.Mass loss ratio (computed using Eq. ( 28) as n Mg /n t=0 Mg ) in Fig. 11 (a) shows that pitting corrosion is initiated immediately after immersion in SBF due to the initial plastic strains, whereas uniform corrosion progresses more slowly.After 24 h of immersion in SBF, pitting corrosion returns a slightly higher mass loss ratio than uniform corrosion.Although Fig. 11 (a) indicates that the stent dissolves faster in the presence of mechanical fields, pitting corrosion notably deteriorates the structural integrity of the stent, as further elaborated.Phase-field isosurface plots after 24 h of immersion in solution for the first case study considered are presented in Fig. 11 (b).A pitting zone is observed in the vicinity of the union between rings and links.The dissolution rate within the pitting zone is much higher than in the remaining parts of the stent.This locally enhanced dissolution significantly reduces the thickness of the strut, as shown in two characteristic cross-sections close to the union point in Fig. 11 (c), indicating hot spots for early stent failure.The structural integrity of the stent at these locations is severely undermined.In the case of uniform corrosion, the contour plots in Fig. 11 (c) show that the stent gradually dissolves, covering the whole sample with a constant dissolution rate.This example demonstrates the importance of including mechanical fields in analyzing the degradation of coronary Mg stents.These structures inevitably experience complex stress states during deployment and service, and thus, uniform corrosion models would give overestimated service life predictions.

Discussion
The present diffuse interface model for assessing the in vitro corrosion of biodegradable Mg-based alloys captures different corrosion mechanisms.Uniform corrosion is included through a constant interfacial mobility parameter calibrated with in vitro corrosion data in terms of hydrogen gas evolution (or mass loss).A spatially-dependent interfacial mobility parameter ( Eq. ( 29) ) is The potential of the model for capturing the role of mechanical fields in enhancing the corrosion of Mg alloys is demonstrated by modeling the behavior of a circumferential sample containing a notch and undergoing tensile testing ( Section 4.1 ).The mechanical contribution is incorporated via a mechano-electrochemical effect that depends on local stress and strain distributions, Section 2.4 .Mechanical stresses have deleterious influences on the corrosion resistance of Mg alloys, as previously observed in several in vitro studies [19][20][21][22] , leading to the localization of damage and the formation of sharp cracks that accelerate failure.Changes in pit morphology, increased hydrogen gas production, pit-to-crack transition initiation, and faster crack propagation are noticed under external loads in Figs. 8 and 9 .More importantly, the model shows that once the pit-to-crack transition develops, it leads to rapid and uncontrollable crack growth.For practical purposes, the model can be utilized in designing and estimating the service life of load-bearing biomedical devices from Mg alloys.Moreover, it can serve as an effective tool to foresee the mechanical strength of body implants (i.e.scaffolds for bone tissue engineeringas well as fixation devices for bone fracture) after a certain period of degradation depending on their geometry and to preempt catastrophic implant failures.
The proposed framework is also used to assess the effect of complex stress conditions, which arise during stent deployment and service, on the corrosion of bioabsorbable Mg stents.Pitting corrosion, initiated due to local plastic strains developed during stent deployment ( Fig. 10 (c)), proves to have more detrimental effects on stent degradation than uniform corrosion, Fig. 11 .The model may serve as a cost-effective way of predicting the degradation of Mg stents and assessing their residual strength during the degradation process.The availability to foresee the locations of early break points in the sample and determine scaffolding capabilities during degradation is appealing for practical applica- tions in the design of biomedical devices such as bioabsorbable stents.Obtaining an optimized stent design is beyond the scope of the current paper.However, integrating the proposed model with optimization analysis would return more sophisticated stent designs with improved corrosion performance and help develop new bioabsorbable metallic stents.
The model potential thus being discussed, it is important to address the questions regarding the role of the equilibrium concentration in the liquid phase and the rate-limiting process on the corrosion behavior, showing additional model capabilities.The formation of the partly protective layer is neglected in the present formulation.The saturated concentration of Mg ions in the corrosive environment is determined based on the mass density and molar mass of MgCl 2 formed on the exposed Mg surface ( Section 2.1 ).This yields the equilibrium concentration of Mg ions in the liquid phase c l,eq Mg = 0.57 mol/L.Taking a higher value for the saturated concentration would lead to difficulties in forming the protective layer, thereby promoting corrosion.On the contrary, decreasing c l,eq Mg would physically represent easier precipitation of the protective layer on the metal surface, increasing corrosion resistance and consequently decelerating the degradation process.To show the effect of c l,eq Mg on the corrosion process, two additional case studies are considered with lower c = 0 .5 c l,eq Mg and higher c = 2 c l,eq Mg equilibrium concentrations in the liquid phase while keeping all the other parameters fixed as in Section 4.1 .Phase-field contours and Mg concentration distributions for the final pit shape are shown in Fig. 12 .As expected, the pit depth increases with the equilibrium concentration in the liquid phase.Hence, the proposed framework can be tweaked to capture the formation of the protective film or other phenomena related to Mg surface modifications by varying c l,eq Mg .
The rate-limiting process between diffusion-and activationcontrolled corrosion is defined in Eq. (25) .Considering the geometry and material properties as in Section 4.1 , two corrosion tests are conducted to illustrate the effect of the rate-limiting process on corrosion behavior.Phase-field contours for the final pit shape and Mg concentration distribution in the absence of mechanical load are shown in Fig. 13 for both rate-limiting processes.In agreement with expectations for the diffusion-controlled process ( τ 1), the pit growth is pronounced and Mg concentration around the interface is close to the equilibrium value in the liquid phase c l,eq Mg .On the contrary, pit growth is slower and Mg concentration stays significantly below c l,eq Mg for the activation-controlled process ( τ 1), Fig. 13 .
The present phase-field formulation overcomes the limitations in tracking the evolution of corrosion interfaces in arbitrary domains under complex physics and handling complex topological changes without requiring ad hoc criteria.The current paper focuses on biodegradable Mg alloys due to their high attractiveness as biomaterials.However, the framework developed is general and  easily extendable to other biodegradable metals, such as Fe and Znbased alloys, using the corresponding material properties and the composition of the corrosion layers.The advantages of the phasefield method can be further exploited by extending the present formulation and adding other physical phenomena, such as the electrochemistry-corrosion interplay.As emphasized in Section 2.1 , the reactions for product formation and layer dissolution ( Eqs. ( 1) and ( 2) ) are neglected in the current model.Thus, the degradation mechanism is based on the anodic reaction and diffusion of Mg ions in the physiological environment.The contribution of the electric field to material dissolution and species diffusion is also neglected in the present model, as the Mg ions are not considered as charged species.These limitations could be overcome by incorporating the reactions for product formation and layer dissolution along with the transport of charged ions, their interactions, and electric field distribution.This would make the model more advantageous and enhance its versatility.Such a model would contribute to understanding the underlying electrochemical process and disclose the effect of the composition of the environment on the corrosion process.The extension to incorporate the electrochemistrycorrosion interplay will be addressed in future works.Incorporating the above-mentioned ingredients could potentially solve the long-term open question of the mismatch of corrosion rates between in vivo and in vitro tests.In addition, future work should consider the effect of microstructural features, such as alloying elements, grain size/shape, grain boundaries, and interfacial energy dependence on grain orientation, on Mg corrosion.These features would deliver new scientific insight into other phenomena related to Mg corrosion, such as intra-and trans-granular corrosion.

Conclusions
A computational framework based on the phase-field method has been presented for assessing the corrosion of biodegradable Mg alloys in physiological environments that resemble biological media.Built upon thermodynamical principles, the model uses an Allen-Cahn equation to capture Mg dissolution, a diffusion equation to estimate the diffusion of Mg ions in solution, and a mechano-chemical enhancement of the phase-field mobility coefficient to capture the interplay between corrosion kinetics and mechanical fields.In addition to uniform corrosion, pitting corrosion is introduced assuming nonuniform distributions of chloride ions in solution.The proposed framework applies to arbitrary two-dimensional and three-dimensional geometries with no special treatment for the evolution of the corrosion front.
The model parameters for uniform corrosion are calibrated with in vitro corrosion data and predictions of pitting corrosion are compared with experiments conducted on Mg wires.A good agreement between experiments and simulations is retrieved.The importance of including the effect of pitting corrosion mechanism and mechanical loads in accelerating degradation is demonstrated in representative case studies: pitting corrosion associated with the local failure of a protective layer and the nonhomogeneous stress state of a bioabsorbable coronary stent.The following conclusions can be drawn: (i) Mechanical loading significantly alters corrosion kinetics and has deleterious influences on the corrosion resistance of Mg alloys.The application of tensile deformation changes the pit morphology and initiates a pit-to-crack transition.Further increase in mechanical loading may trigger rapid crack growth and premature fracture after a short time in SBF, as previously observed in in vitro studies.(ii) Local plastic strains developed during stent deployment act as initiators for pitting corrosion, indicating hot spots for early stent failure.The results show that pitting corrosion in the stent is initiated immediately after immersion in SBF due to the initial mechanical strains, whereas uniform corrosion progresses more slowly, covering the whole sample with a constant dissolution rate.In addition, pitting corrosion proves to have more detrimental effects on stent degradation than uniform corrosion and severely compromises the structural integrity of the stent.This study reveals that neglecting the mechanical effects on stent degradation and considering uniform corrosion would lead to unsafe design solutions and overestimated service life predictions of coronary Mg stents.The proposed framework can assist in designing and predicting the service life of biomedical devices after a certain period of immersion.The model may serve as a complementary tool for planning in vitro, tests and as a cost-effective way of assessing the residual strength and scaffold capabilities of temporary body implants such as bioabsorbable stents, bone implants or porous scaffolds for tissue regeneration.

Declaration of Competing Interest
The authors declare that they do not have competing financial interests or personal relationships that could have influenced the work reported in this paper.

Fig. 1 .
Fig. 1.Schematic illustration of the pitting corrosion (left) and stress corrosion cracking (right) mechanisms of Mg alloys during immersion in the physiological environment (simplified representation).Pitting is caused by breakages in the passive film exposing the Mg alloy to the corrosive environment.
is enforced to retard diffusion of Mg ions inside the solid phase.

( 21 )
Other dimensionless fields and parameters are c Mg = c Mg /c s Mg c l,eq Mg = c l,eq Mg /c s Mg c s,eq Mg = c s,eq Mg /c s Mg D c Mg = D c Mg /D l c Mg

Fig. 3 .
Fig. 3. Schematic disposition of the experimental setup (left) and the corresponding nondimensional computational domain (right) for WE43 Mg alloy wires immersed in SBF.The size of the nondimensional computational domain ( w s = 37.5, w l = 362.50,and h = 250) is normalized using the interface thickness = 4 μm as the characteristic length.
Value Unit Diffusion coefficient of Mg ions in the liquid phase D l c Mg 10 −10 m 2 /s Diffusion coefficient of Mg ions in the solid phase D calculated based on the mass density and molar mass of MgCl 2 formed on the exposed Mg surface.Assuming that the mass density of MgCl 2 is 54.20 g/L and its molar mass of 95.21 g/mol, the equilibrium concentration of Mg ions is c l,eq Mg = 0.57 mol/L.The role of the saturated concentration in the degradation process is further addressed in Section 5 .

Fig. 4 .
Fig. 4. Hydrogen gas evolution as a function of immersion time for Mg wires.Numerical results assuming uniform corrosion and experimental measurements.The light blue area stands for the standard deviation of the experiments.
Fig.5 (a)-(c).The corresponding 2D contour plots of the remaining cross-section of Mg after 24 h of immersion in SBF are given in Fig.5 (d)-(f).Pitting corrosion initiates and follows regions with high L 0 values (i.e., more Cl − ions).Three representative experimental cross-sections of the Mg wires after 24 h of immersion in SBF are plotted in Fig.5 (g)-(i) for the sake of comparison.They are

Fig. 5 .
Fig. 5. (a-c) Spatial distribution of the mobility parameter L 0 generated with N = 2 and β = 0 . 1 for three different α values.(d-f) Phase-field predictions of the crosssection of the Mg wire after 24 h of immersion in SBF using L 0 from (a-c).(g-i) Representative experimental cross-sections of the Mg wires after 24 h of immersion in SBF.The grey points and grey lines indicate the center and initial cross-section of the Mg wire before degradation in both experiments and simulations.The red circle stands for a uniform corrosion radius.The scale bar for all figures is 50 μm.

Fig. 6 .
Fig. 6.Comparison between experimental and simulated values for the three pitting metrics parameters.(a) Uniform corrosion radius (the radius of the circular section that has the same area as the corroded cross-section).(b) Average pit depth (average distance from the degraded cross-section to the uniform corrosion circle).(c) Maximum pit depth (maximum distance from the corroded cross-section to the uniform corrosion circle).

Fig. 7 .
Fig. 7. Simulation domain for the Mg alloy wire of 1 mm in length and 300 μm in diameter immersed in SBF and subjected to tensile deformation along the wire axis.The semi-circular pit created by the rupture of the protective layer has a radius of 10 μm.The size of the nondimensional computational domain ( w s = 37.5, w l = 362.50,and h = 250) is normalized using the interface thickness = 4 μm as the characteristic length.

Fig. 8 .
Fig. 8. Contour plots of the phase-field variable, Mg concentration distribution in SBF, hydrostatic stress σ h , and effective plastic strain ε p for various prescribed remote deformations ε ∞ after 24 h of immersion in SBF.The initial surrounding corrosive environment is not shown in the plots.

Fig. 9 .
Fig. 9. Model predictions of (a) pit depth and (b) hydrogen release as a function of immersion time for different applied axial strains.

Fig. 11 .
Fig. 11.Pitting and uniform corrosion of bioabsorbable coronary Mg stents.(a) Mass loss ratio as a function of immersion time.(b) Phase-field isosurface plots for pitting corrosion after 24 h of immersion.The pitting zone is observed in the areas of high plastic strains.(c) Phase-field contour plots of two characteristic cross-sections after 24 h of immersion.The red line indicates the initial cross-section of the Mg stent before degradation.

Fig. 12 .
Fig. 12. Contour plots of the phase-field variable and Mg concentration distribution in liquid for different equilibrium concentrations of Mg ions in the liquid phase ( c l,eq Mg ) after 24 h of immersion in SBF in the absence of mechanical loading.The simulation domain corresponds to Fig. 7 .The surrounding corrosive environment is not shown in the plots.

Fig. 13 .
Fig. 13.(a) Contour plots of the phase-field variable and Mg concentration distribution in liquid for diffusion-controlled ( τ = 10 6 ) and activation-controlled corrosion ( τ = 10 −6 ) after 10 h of immersion in SBF.The simulation domain corresponds to Fig. 7 .The surrounding corrosive environment is not shown in the plots.(b) Pit depth as a function of immersion time for the two rate-limiting processes.