A Physical Approach to Simulate the Corrosion of Ceramic-Coated Magnesium Implants

: Magnesium-based biodegradable materials are currently of great interest in various biomedical applications, especially those related to the treatment of bone trauma and the manufacturing of bone implants. Due to the complexity of the degradation process of magnesium, several numerical models were developed to help predict the change of the implant’s integrity in the body using in vitro tests. In this study, experimental in vitro tests and ﬁnite element methods are combined to calibrate a diffusion-based model of the uniform galvanic corrosion of high purity magnesium (HP-Mg). In addition, and for the ﬁrst time, the impact of a porous coating layer generated by the Micro Arc oxidation (MAO) method is investigated and incorporated into the model. The calibrated model parameters are validated using the same immersion test conditions on a near-standard of treatment screws geometry made of HP-Mg.


Introduction
Bone fractures have a high incidence and are treated by osteosynthesis, which includes implanting fixation devices to hold the adjacent bone fragments still until the completion of the healing process [1]. The market size for these devices in the United States reached $6.9 billion in 2017 [2]. Bio-inert metals with high corrosion resistance such as stainless steel and titanium alloys have been the conventional materials used to make these fixation devices due to their strength and biocompatibility. The downfall is that another surgery is usually needed to remove the implant after healing as they become redundant and may cause severe complications. In other cases, where bone defects are difficult to treat by osteosynthesis alone, bone grafting techniques are used [3]. In auto-grafting for instance, bone tissue is extracted from a different healthy bone to fill the defected site and provide support to the healing bone tissue [3]. However, it becomes challenging as the defect size increases. Therefore, bioabsorbable alternatives grabbed the attention in the past two decades as base materials for osteosynthesis devices as well as bone scaffolds for bone grafting. Magnesium alloys are given more research interest due to their superior biocompatibility and closest range of mechanical properties compared to natural bone tissue and other used metals for orthopedic implants (see Figure 1) [4].
However, magnesium is the least noble metal in the galvanic series with an electrode voltage of −2.372 V against standard hydrogen electrode [5]. Hence, alloying and surface treatments were investigated in the literature to improve the corrosion resistance and strengths of pure magnesium [6]. Moreover, ceramic coatings made by the promising Micro Arc oxidation coating method have been shown to improve the corrosion resistance significantly for magnesium and its alloys [7][8][9]. As a result, the corrosion behavior of magnesium-based implants is challenging to predict in an aqueous environment, and it requires extensive in vitro and in vivo corrosion experiments to best predict this behavior. In order to provide an inexpensive method to predict this corrosion behavior, numerical simulations have been introduced over the past decade to describe the different corrosion mechanisms found in magnesium alloys. For instance, a phenomenological approach based on the continuum damage theory was used to describe the uniform corrosion as well as the heterogeneous pitting corrosion combined with the stress corrosion cracking [10][11][12][13]. Physical modeling of the corrosion was also used to model a reactioncontrolled corrosion behavior as well as transport-controlled corrosion using diffusion law [14][15][16][17][18]. In the model by Grogan et al. [18], a transport-controlled corrosion was assumed to carry out the model for a uniform corrosion behavior in pure magnesium. The transport was assumed mainly diffusive and follows Fick's law of diffusion in Equation (1).
where c Mg is the magnesium ions concentration, and D Mg is the diffusion coefficient of the magnesium ions in the environment, which was assumed isotropic. In this model, diffusion is assumed following the saturation of the ions near the corrosion surface. Subsequently, the corrosion surface is moved inward to model the loss in the volume due to corrosion. The surface movement was modeled in an arbitrary Lagrangian-Eulerian adaptive mesh and finite element method. Bajger et al. extended the diffusion-controlled model by introducing a coupled system of partial differential equations to capture the interaction between the corrosion products and chloride species with magnesium ions [14]. The concept of the effective diffusivity of the magnesium ions through the corrosion products was also introduced. The level set method was used to model the corrosion surface movement in a 2D finite element model. Besides the diffusion coefficient in Equation (1), the controlling parameters in the model include the saturation limit of the magnesium ions in the environment and the initial ions concentration in the solid magnesium and the environment.
The initial concentration of the magnesium ions in the environment is well established and is used as the concentration level in the blood plasma. The concentration in the solid is assumed equivalent to the partial density of the magnesium in an alloy or the total density for pure magnesium. The saturation limit is used following the assumption that the transport of the ions follows the saturation of the environment on the corrosion boundary. This value was approximated in the previous models by the saturation limit of the anhydrous magnesium chloride in water at room temperature [14,18]. The diffusion coefficient was not available prior to the work by Bajger et al. [14] where they calibrated the model against mass loss data collected using hydrogen evolution testing in the study by Abidin et al. [19]. Although the hydrogen evolution test can provide more data points along with the test duration, it was found to highly underestimate the corrosion rate as compared to weight loss data [19].
In this study, a physical and diffusion-controlled corrosion model was adopted to describe the uniform galvanic corrosion observed in pure magnesium. The saturation limit of the magnesium ions was investigated experimentally for the first time in a simulated body environment to be used in the subsequent calibration of the diffusion coefficient. The calibration of the diffusion coefficient was based on the weight loss data which was found to best represent the actual corrosion in magnesium [20]. Furthermore, the effective diffusivity concept was utilized to study the effect of adding a porous coating layer of a given thickness.

Saturation Tests
The modified simulated body fluid (m-SBF) prepared according to the procedure described in [21] is one of the solutions that best mimics the level of inorganic compounds in the blood plasma. The solution is prepared using ultrapure water and reagent grade salts including NaCl (5.403 g/L), NaHCO 3 (0.504 g/L), Na 2 CO 3 (0.426 g/L), KCl (0.225 g/L), K 2 HPO 4 .3H 2 O (0.23 g/L), MgCl 2 .6H 2 O (0.311 g/L), HEPES (buffering agent, 17.892 in 100 mL of 0.2M NaOH), CaCl2 (0.293 g/L), Na 2 So 4 (0.072 g/L), and 1M NaOH (15 mL/L). To prepare a liter of the solution, all glassware were cleaned using neutral detergent followed by 1M HCL and eventually ultrapure water. The salts were dissolved in the same order in 500 mL of ultrapure water under heating and stirring using a magnetic stirrer and heating plate. The pH was adjusted to 7.4 at 37 • C, and then ultrapure water was added to complete the 1 L volume.
Ten milliliters of this solution was tested for its saturation limit of MgCl 2 under a body temperature of 37 • C. Reagent grade Anhydrous MgCl 2 was used as an approximation to magnesium ions, the species under concern, which was also used to provide an estimation for the solubility limit in water in the literature. Ten milliliters of the solution was put in an open glass beaker, and a heating and magnetic stirrer plate was used to hold the temperature at 37 • C and dissolve the added increments of the salt. An amount of the salt was placed in an air-tight HDPE 250 mL bottle to be used for adding the salt until saturation. The bottle and the salt weight were measured before and after running the test using a precision scale with an accuracy of 1 mg. MgCl 2 salt has an exothermic dissolution nature in water, which increases the solution temperature above the target temperature at each increment. Hence, the test solution was left to cool down before adding the next increment. The process was repeated until visible crystals of the salt were no longer dissolving under stirring. The test was repeated for 10 runs. The separation of the precipitates was not possible, and their weight was not deducted due to the precipitation of other compounds than the added salt.
Another approach was also used based on thermal feedback (see Figure 2), rather than visual feedback that relies on precipitation of non-dissolving crystals. The exothermic dissolution nature of the MgCl 2 was utilized by monitoring the temperature change above 37 • C after each salt increment. When the change was insensible by the used thermometer, saturation was assumed. This approach was repeated for another 10 runs.

Micro Arc Oxidation Coating
Round coupons of 10 mm in diameter and 3 mm thickness were cut from a rolled HP-Mg rod (Goodfellow, Coraopolis, PA, USA). A hole of~1.5 mm in diameter and~3 mm depth was machined radially to connect the samples with a copper wire that delivers the coating current. The connection point was covered and fixed by a non-conductive epoxy. An in-house setup shown in Figure 3 was used to for coating the samples. The coating process parameters developed in a previous study for an Mg-1.2Zn-0.5Ca-0.5Mn alloy were used except for a higher current density of 100 mA/cm 2 [22]. A pulsed DC current was applied at 5000 Hz and a constant amplitude, which lead to an increase of the applied voltage up to 450 ± 30 V. The electrolyte was 3 g/L (NaPO 3 ) 6 and 8 g/L KF-2H 2 O prepared using high purity reagent grade chemicals and distilled water.

In Vitro Immersion Test
The immersion test in m-SBF was used to collect mass loss (ML) data of high purity magnesium (HP-Mg) due to corrosion for both the bare and coated sample groups. Round coupons of 10 mm in diameter and 3 mm thickness were cut, dry sanded using SiC sandpapers of 400-2000 grit size, cleaned in an ethanol bath, and then dried using hot air for 2 min. The weight measurements for both groups were recorded using a 1 mg accuracy precision scale. The samples were then immersed in 50 mL/cm 2 of m-SBF to eliminate the solution volume effect [20]. The buffering reagent is HEPES in the m-SBF to hold the pH around the physiological value of 7.4. An automatic pH controller (Bluelab, New Zealand) was also used to adjust the rising pH by dosing 1M HCl. The test was conducted inside a thermal incubator (Thermo Scientific Heratherm) to hold the temperature at 37 • C. The test was conducted for 7, 14, 21, and 28 days (n = 3) for both groups. By the end of each period, the samples were collected and cleaned in chromic acid (150 g/L chromium trioxide, 10 g/L silver chromate, and ultrapure water) according to the ASTM G1-03 standard to remove the corrosion products [23]. Samples were then rinsed in distilled water followed by ethanol cleaning and drying before recording the new weight. In order to guarantee symmetric exposure of the round samples in the test solution, the samples were mounted floating by two concentric nylon screws along the sample axis. The mounting setup is shown in Figure 4 for two samples. The jigs and screws were made of inert nylon to eliminate alteration of the solution chemicals and galvanic corrosion with the test samples.

Numerical Model Development
The corrosion process in magnesium alloys can be summarized as two stages; the first short-term stage is the reaction-controlled period governed by the total chemical reaction in Equation (2). At this stage, corrosion products of magnesium hydroxide form on the surface and provide partial protection against the corrosive environment. Following this stage, the rate of the transport of the ions through the corrosion surface becomes slower than the rate of the corrosion product formation. Hence, the corrosion becomes transport controlled afterwards [18], see Figure 5. Following the model by Grogan et al. [12], Fick's law of diffusion was also used to describe the uniform corrosion in HP-Mg in this model. The diffusivity coefficient was considered constant and isotropic, leading to Equation (3). The thickness of the developed corrosion products was assumed constant, and the kinetics of the electrochemical reaction were neglected. The volume of the solid magnesium was assumed constant and the concentration inside the volume was allowed to decrease following the diffusion. Initial values in the solid and fluid domain were equivalent to the pure magnesium density of 1732 Kg/m 3 and the magnesium concentration in the blood plasma of 0.03 Kg/m 3 [21]. Dirichlet-type boundary conditions were imposed on the solid-fluid interface equivalent to the saturation limit of the magnesium ions in the environment, calibrated in this study, and on the fluid outer boundary equivalent to the magnesium concentration in the blood plasma of 0.03 Kg/m 3 . A schematic of the model domains and the boundary conditions is shown in Figure 6. The equation was solved using the finite element method included in Ansys APDL, with the magnesium ions concentration as the scalar field variable. The numerical mass loss data was calculated according to Equation (4).
where V 0 is the sample volume, ρ 0 is the initial density of the pure magnesium, ρ t is the temporal density which was calculated by the element area-weighted average of the magnesium concentration in the 2D model or by the volume-weighted average of a 3D model, and A 0 is the exposed surface area of the implant. From the uniform corrosion nature and the controlled symmetric exposure in the immersion test, the model could be reduced to a 2D axisymmetric model, Figure 7a, to reduce the element count and improve the mesh resolution at the corrosion surface. The mesh was built using rectangular elements with 8 nodes. The solution was found to converge at the minimum element size of 10 µm at the high concentration gradient region of the sample boundary. The m-SBF domain dimensions were chosen to replicate the solution volume to surface ratio in the immersion test. Following the calibration of the bare HP-Mg model for C sat and D Mg , the impact of the coating layer was modeled using a fine layer of 10 µm thickness surrounding the sample region having a reduced diffusivity coefficient, Figure 7b. The concept was previously used in the model by Bajger et al. [14] for modeling the porous corrosion product's impact on the diffusivity coefficient in a transport-controlled corrosion model. The coating layer was assumed completely porous through its thickness, and the transport of the ions follows the effective diffusivity coefficient (Equation (5)) [24].
where is the porosity and τ is the tortuosity of the porous ceramic coating layer. Due to the lack of information on these parameters in the literature, they were grouped in a single parameter β to be calibrated by experiment.

Model Validation
Further immersion test data of a near-standard ankle and foot orthopedic screw (3.5 mm diameter), shown in Figure 8, was used to validate the calibrated values from the round coupons' geometry. The screws were machined out of the same HP-Mg as rolled rod using a CNC lathe. A 2D axisymmetric finite element model was used again with 6 triangular node elements to resolve the model geometry with acceptable elements quality.

Model Parameters Calibration
The saturation limit measurements of the MgCl 2 measured by the two approaches for 10 runs are compared using a normal distribution curve as shown in Figure 9 to show the difference in the variation of the measurements. The thermal feedback approach yielded a lower standard deviation from the mean value of 342.27 Kg/m 3 compared to the mean of 377.47 Kg/m 3 from the visual feedback approach. The difference can be attributed to the excess added salt after the saturation of the solution that is not visible as precipitated crystals and was suspended in the solution. The measured value is~255% of the used value in the literature, that is, 134 Kg/m 3 in water at room temperature [14,18]. The differences can be attributed to the difference in the temperature and the composition between water and the m-SBF solution of ultrapure water and the dissolved inorganic salts. The measured saturation limit by the thermal feedback approach was used in the next step of the model calibration for the magnesium ions diffusion coefficient in the m-SBF. Experimental mass loss per exposed surface area, Equation (6), of bare HP-Mg from the immersion test over 4 weeks of immersion was compared to the numerical mass loss calculated from the 2D axisymmetric model.
where m i is the initial mass before immersion, m f is the final mass after immersion, and A is the initial exposed surface area and assumed constant. Given the boundary conditions and the initial values of the concentration in the solid and the fluid domains, the last parameter is the diffusion coefficient to be calibrated to fit the experimental data. Since it is a one-parameter fitting process, a good match is only feasible to a portion of the data by shifting the numerical data up and down following the change in the diffusion coefficient. Contour plots of the concentration are shown in Figure 10. Figure 11 shows the experimental data of the mass loss compared to the numerical data at three levels of the diffusion coefficient; upper level (most conservative) D Mg = 7 × 10 −10 m 2 /day, lower level (least conservative) D Mg = 1 × 10 −10 m 2 /day, and an average level of D Mg = 4 × 10 −10 m 2 /day. It was shown that the upper level of the diffusion coefficient provides a good match to the experiment by the end of week 3 and along week 4. This result confirms the fact that the corrosion in magnesium is reaction controlled at the early stage and shifts to transport-controlled corrosion in the long term [18].  Given the saturation value Csat, and the upper limit diffusion coefficient from the bare HP-Mg model calibration, the impact of the coating layer (CL) is investigated. The mass  Given the saturation value C sat , and the upper limit diffusion coefficient from the bare HP-Mg model calibration, the impact of the coating layer (CL) is investigated. The mass loss data of coated samples are compared against the numerical model, and the reduction factor β is varied to fit the data. Figure 12 shows the result of fitting the experimental data at a reduction factor of 0.0714. The data were found to best match again towards the end of the test duration from the start of week 3 to the end of week 4. The used coating layer thickness was assumed equivalent to the average value of the previous study using the same process parameters [22]. It was found that varying the process parameters of the Micro Arc oxidation coating can alter the developed ceramic layer thickness [25]. From the calibrated values, the effect of increasing the coating layer thickness was investigated by calculating the mass loss per exposed surface area at the end of 4 weeks at different coating layer thicknesses (10-100 µm). Figure 13 shows that the mass loss data for the simple geometry of the round coupon follows the power rule. This trend was expected to vary according to the geometry being modeled with a general reduction behavior. In the case of the coated round coupons, it was expected that the mass loss after 28 days would be reduced by 45% by increasing the coating layer thickness up to 50 µm. A further increase in the thickness was found to have an insignificant effect on the mass loss value. Figure 13. The mass loss per exposed surface area after 28 days of immersion plotted vs. the thickness of the coating layer.
The collected immersion mass loss data for the coupons are made available in Table 1. The calculated corrosion rate P W using Equation (7) [19] from the mass loss rate mg/cm 2 .day after 28 days was found to be 5.2 mm/year and 3.4 mm/year for bare and coated samples, respectively. Assuming the same rates apply in vivo and that uniform corrosion takes place, it is expected that for a standard orthopedic ankle and foot screw to be absorbed completely by 4 months and 5 months for the bare and MAO-coated screws, respectively.  Figure 14 shows the corroded screws of bare and coated HP-Mg screws under the same conditions of the immersion test in a simulated body environment. Visual inspections show the superior performance of the MAO-coated screws in maintaining the head and thread geometry at the same duration. The FEA model results of the screws are shown in Figure 15. The bare screw model shows wider spread in the contour plot of the concentration near the screw surface than the coated screw model. This result reflects that more magnesium ions had diffused into the environment more than for the coated screw. The numerical mass loss data for bare screw were compared with the collected experimental mass loss data in Figure 16. The transport numerical model was found to best match the experimental data towards the end of the test period of 4 weeks. This result validates the ability of the model to predict the mass loss of the ceramic-coated HP-Mg implants to good accuracy on the long term.

Conclusions
The diffusion-based physical model of the uniform galvanic corrosion in HP-Mg is calibrated in this study. Saturation tests were conducted to find the saturation limit of the magnesium ions in simulated body fluid under physiological conditions. The magnesium ions diffusion coefficient was calibrated for the first time against mass loss data under controlled pH and physiological conditions. Furthermore, the impact of the added porous coating layer on the effective diffusion coefficient in the coating region was measured. The calibrated model was validated by using the repeated immersion test of near-standard orthopedic screws. The calibrated diffusion-based model predictions are limited to longterm corrosion which is approximately after the first two weeks of corrosion. The model by itself cannot predict the geometry change directly; yet, it offers an estimation for the massloss rate for the geometry under question. This rate can be easily converted to the corrosion front speed under the assumption of uniform corrosion to monitor the geometry change at the sections of concern. Further investigation utilizing ALE or the level set method and the developed parameters in this work is needed to provide both accurate qualitative and quantitative simulations for the corrosion in magnesium-based implants. Geometric change due to corrosion is crucial to the mechanical analysis of the implant structure and fatigue studies. Parallel experimental assessment of the mechanical properties of the corroded samples is of great importance to the development of empirical equations for the mechanical properties as functions in time.