Chapter 3 Developing the core computational model

Due to their bio-friendly properties, biodegradable metallic biomaterials, including magnesium (Mg), iron (Fe), and zinc (Zn), are regaining attention in recent years [299]. These biomaterials find important applications in the design and manufacturing of supportive implants such as temporary devices in orthopedics and the cardiovascular field [55, 296]. In orthopedics, the biodegradable metallic biomaterials are used as fixation devices, providing adequate support in the early stages while being absorbed gradually during the bone healing process [209]. Implants fabricated using Mg and its alloys are being used for such a purpose [216] due to the similarity of the stiffness between natural bone and Mg, which helps to reduce the stress shielding induced by the implanted device. Additionally, Mg is reported to have a non-toxic contribution to the human body’s metabolism and the bone healing process, which makes the release and absorption of metallic ions safe and biocompatible [283].

Accumulation of mechanistic understanding of Mg degradation achieved by experimental approaches over the years gradually provided a mechanistic understanding of the biodegradation process. Combining these insights with in silico modeling approaches enables researchers to study the biodegradation properties and behavior of the implant in a virtual environment prior to conducting any in vitro or in vivo tests. When fully validated, computational modeling can (in part) replace certain stages of costly and time-consuming experiments verifying the expected degradation behavior of the designed implants. Additionally, the developed models can be efficiently combined with existing computational models to examine other related phenomena such as tissue growth or mechanical integrity.

Previous contributions to the computational modeling of the degradation process include a wide range of different approaches, from the basic phenomenological implementations to comprehensive mechanistic models that take into account various aspects of the degradation and resorption process.

Continuous damage (CD) modeling has always been a common approach for corrosion simulation, but from a physicochemical point of view, it focuses on the mechanical integrity of the degradation and neglects the diffusion process. As a result, its application in the degradation modeling of biomaterials, which includes various fundamental phenomena such as mass transfer through diffusion and reaction, is relatively limited. Despite this issue, a CD model proposed by Gastaldi et al. showed a good performance for simulation of bioresorbable Mg-based medical devices [86], in which geometrical discontinuities were interpreted as the reduction of material.

Alternatively, mathematical modeling using transport phenomena equations has shown great flexibility in capturing different mechanisms involved in the biodegradation process. As an example, in Ahmed et al., a set of mathematical equations in cylindrical and spherical coordinates was derived to model the chemical reactions of Mg degradation [7]. Despite the simplicity of their approach from the computational perspective, their model was able to demonstrate the contribution of various chemical components to the in vitro degradation of Mg. Similarly, Grogan et al. developed a mathematical model based on the Stefan problem formulation in 1D space to correlate the mass flux of metallic ions into the solution to the velocity of shrinkage of the material during degradation [98]. This was done by considering the mass diffusion and change of the concentration of $\mathrm{Mg}^{2+}$ ions, and then, employing an arbitrary Eulerian-Lagrangian (ALE) approach to extend the model to 3D on an adaptive mesh. A similar approach was taken by Shen et al. to develop a theoretical model of the degradation behavior of Mg-based orthopedic implants showing great consistency with in vitro test results [236].

An ultimate application of the computational modeling of the biodegradation process of biomaterials can be the prediction of how biodegradation affects the shape of the bulk material, medical device, or implant over time. One of the ways to achieve such a prediction is to capture the movement of the corrosion front mathematically using an appropriate method. The level set method (LSM) is a widely used example in this regard, which is an implicit mathematical way of representing the moving interfaces. This approach was used in Wilder et al. to study galvanic corrosion of metals [275]. They employed LSM on an adaptive mesh to track the moving corrosion interface, but their model lacked a thorough validation using experimental data. Gartzke et al. also worked on a simplified representation of the interface movement by developing a mechanochemical model of the biodegradation process, which helped them to study the effect of degradation on the mechanical properties [84]. They performed a basic qualitative validation on the predictions made by the model. Another similar study in this regard is the Sun et al. work [250], in which a detailed mathematical model was derived and validated to study the deposition of corrosion products on the surface of materials. This mathematical approach was also employed in the biomedical field by Bajger et al. to study the mass loss of Mg biomaterials during biodegradation [23]. They used LSM as well as a set of reaction-diffusion equations to track the change of geometry, which can be directly correlated to the loss of material. The derived equations were also able to capture the formation of the corrosion film that decreases the rate of degradation. Another comprehensive mathematical model was developed by Sanz-Herreraa et al. to investigate the role of multiple chemical components involved in the in vitro degradation of Mg implants [227]. One important drawback of this study was its 2D nature. Although the computational model was capable of studying the effect of multiple components, due to the high number of derived equations, it would be difficult to extend and use the same model for real 3D implants. Additionally, a 2D model cannot capture the full phenomenon of corrosion, and as a result, the validation of the model will be more qualitative. It was shown in the study conducted by Gao et al. [82], where they compared the results of a multi-dimension model of the degradation of cardiovascular stents with those of a single-dimension model, that the number of considered dimensions had an important effect on the model predictions. In the end, it is worth mentioning that no dedicated experiments were performed in the aforementioned studies to validate the constructed mathematical and computational models.

The current study focuses on developing a physicochemical model of the biodegradation process of commercially-pure (CP) Mg biomaterials by continuing the work of Bajger et al. [23]. In this model, a set of partial differential equations (PDE) was derived according to the underlying chemistry of biodegradation, described as reaction-diffusion processes taking place at the interface of the biomaterial and its surrounding environment. The formation of a protective layer, effects of the ions in the solution, and the change in the pH due to the corrosion phenomenon were taken into account in the mathematical model. The corresponding computational model was implemented in a fully parallelized manner. Model calibration and validation were executed using data obtained from the immersion tests performed in saline (NaCl) and simulated body fluid (SBF) solutions.

The biodegradation process can be considered as a reaction-diffusion system [270], in which the ions are released due to the chemical reactions on the surface, and the released ions diffuse through the surrounding solution and materials. These ions can interact with other ions and form new compounds [193]. As the reaction-diffusion systems have been studied in science and engineering for a couple of decades, the analogy with a reaction-diffusion system makes it convenient to construct a mathematical model of the biodegradation process based on the well-established transport phenomena equations [96]. From the mathematical perspective, a reaction-diffusion system is expressed by a set of parabolic PDEs that describe the conservation of contributing chemical species in the studied system.

Moving boundary problems, also called Stefan problems, are the general class of mathematical problems in which the boundary of the domain should also be calculated in addition to the solution of the other equations [59]. Coupling the reaction-diffusion system of biodegradation with a moving boundary problem constructs a mathematical model in which the change of the domain geometry due to the material loss can be correlated to the underlying reaction and diffusion processes of corrosion. As the geometry can be determined accurately, this approach provides a way to measure the mass loss directly by computing the change in the volume of the material. In such a system, the moving boundary is the material-solution interface (corrosion front).

For a 1D corrosion diffusion system, the position of the diffusion interface can be determined by [59]:

in which the $s(t)$ represents the position at any given time, and $s_{0}$ is the initial interface position. $\alpha$ coefficient can be calculated using:

where $[\mathrm{Mg}]_{\text{sol}}$ is the concentration in the solid bulk (i.e. materials density) and $[\mathrm{Mg}]_{\text{sat}}$ is the concentration at which the material is released to the medium. $[\mathrm{Mg}]_{0}$ represents the initial concentration of the metallic ions in the medium, which is usually zero for most corrosion cases.

Eqs. 3.1 and 3.2 can be used to simply track the movement of the corrosion front, which is the employed method in studies like the Gorgan et al. work [98], but apparently, the real-world corrosion problems are 3D and much more complex than the described system.

As will be described later, Eq. 3.1 is used strictly for the first time step of the simulations in low diffusion regimes for calculating the initial velocity of the interface. Generally speaking, a more sophisticated approach, such as the level set method, is required for tracking the interface of complex 3D geometries.

In the current study, the corrosion front is tracked using an implicit function such that the zero iso-contour of the function represents the metal-solution interface. As a common practice, this implicit function is expressed as a signed distance function that defines the distance of each point of space (the domain of interest) to the interface. Such a definition implies that the zero iso-contour of the function belongs to the interface. The level set method provides an equation to declare such an implicit function, $\phi=\phi(\mathbf{x},t),\mathbf{x}\in\Omega\subset\mathbb{R}^{3}$, which can be obtained by solving [219]:

in which $\overrightarrow{V^{\mathrm{E}}}$ is the external velocity field, and $\mathrm{V}^{\mathrm{N}}$ is the value of the normal interface velocity. The last term is related to the curvature-dependent interface movement and is omitted. As the effect of perfusion is neglected in the current study, the term containing the external velocity is also eliminated, resulting in the following simplified form of the level set equation:

By having the normal velocity of the interface ($\mathrm{V}^{\mathrm{N}}$) at each point and solving Eq. 3.4, the interface can be captured at the zero iso-contour of the $\phi$ function.

The chemistry of biodegradation of Mg depends considerably on the surrounding solution and the presence of certain ions [193]. In NaCl solutions, the anodic and cathodic reactions as well as the formation and elimination of side corrosion products can be considered as follows [299]:

Reaction 3.6 is not fully correct from the chemical point of view. In fact, Mg surface is always covered by MgO layer, and $\mathrm{Mg}(\mathrm{OH})_{2}$ forms on top of that either at atmospheric conditions or during the immersion. The integrity of this MgO layer is undermined by $\mathrm{Cl}^{-}$ ions, leading to an increase in degradation rate:

Although $\mathrm{Cl}^{-}$ formally does not participate in reaction 3.7, it reflects the dependence of Mg corrosion rate on $\mathrm{Cl}^{-}$ concentration. This effect on the rate of degradation has been widely expressed as the effect of $\mathrm{Cl}^{-}$ on the $\mathrm{Mg}(\mathrm{OH})_{2}$ in the literature [299, 296]. In the developed model, this effect is used interchangeably by omitting the MgO component, so the protective film formed on the corrosion interface is assumed to contain $\mathrm{Mg}(\mathrm{OH})_{2}$ only. Moreover, it has been shown recently that oxygen reduction reaction also takes place during corrosion of Mg [265, 248, 240]. However, this is a secondary reaction (complementing water reduction) contributing to 1-20% of the total cathodic current depending on the conditions. Hence, it is not taken into consideration in this model. Additionally, the involved chemical reactions are more complicated in SBF solutions due to the presence of further inorganic ions and the formation of a layered precipitate structure [193], but the effect of these ions is currently encapsulated in the reaction rates and the diffusion coefficients of the developed mathematical model. The summary of the considered chemistry to develop the mathematical model is depicted in Fig. 3.1.

To keep track of the concentration changes of various contributing chemical components, we define four state variables for the concentration of $\mathrm{Mg}^{2+}$ ions, protective film ($\mathrm{Mg}(\mathrm{OH})_{2}$), chloride ($\mathrm{Cl}^{-}$) ions, and the hydroxide ($\mathrm{OH}^{-}$) ions:

which are indeed 4 scalar functions of space and time. $\Omega$ denotes the whole region of interest, including both the Mg bulk and its surrounding medium. By doing this, the value of pH at each point of $\Omega$ can be calculated as:

where $C_{\mathrm{OH}}$ implies the activity of $\mathrm{OH}^{-}$. By having the definition of the state variables in Eq. 3.8, the biodegradation of Mg described by Eqs. 3.5 and 3.6 can be represented as a set of reaction-diffusion PDEs:

in which the maximum concentration of the protective film can be calculated according to its porosity ($\epsilon$) [23]:

$D^{e}$ is the effective diffusion coefficient for each component. Due to the formation of the protective film, the diffusion coefficient is not constant and varies from the actual diffusion coefficient of the ions to a certain fraction of it. This fraction can be defined as ${\epsilon}/{\tau}$ [94, 120], in which $\epsilon$ and $\tau$ are the porosity and tortuosity of the protective film, respectively. The effective diffusion coefficient can be then calculated by interpolating the two aforementioned values:

The $\beta$ coefficient is called momentum here and controls the effect of the saturation term $(1-\frac{C_{\mathrm{Film}}}{[\mathrm{Film}]_{\max}})$.

In order to take advantage of the level set method for tracking the corrosion front, the velocity of the interface at each point should be determined. Then, by solving Eq. 3.4, the interface is obtained at the points with a zero value of the $\phi$ function. The interface velocity in mass transfer problems can be calculated using the Rankine–Hugoniot equation [230], and by considering the transportation of $\mathrm{Mg}^{2+}$ ions, it can be written as:

where $\mathbf{J}$ is the mass flux at the interface. Rearranging Eq. 3.16 and inserting the value of the normal interface velocity into Eq. 3.4 yields:

which is the final form of the level set equation to be solved. In the case of simulations with a low diffusion rate, the interface moves slowly in the beginning, which results in a linear degradation, whereas based on the experimental results, the degradation rate is fast at the beginning and slows down eventually [192]. So, to mimic the same behavior in the low diffusion regimes, we took advantage of the theoretical Stefan formulation (Eqs. 3.1 and 3.2) to push the interface in the first time step. According to Eq. 3.1, the velocity of the interface can be calculated as $({2\alpha}/{\sqrt{t}})$, but as we are dealing with a 3D model and not a 1D one, we pick a fraction (denoted by $\gamma$) of this ideal value to be used as the driving force of the interface at the beginning of the simulation. So, the normal velocity of the interface can be written in the general form as:

in which the $\alpha$ value should be calculated from Eq. 3.2. By selecting $\gamma$ equal to zero, the Stefan formulation can be eliminated, and a value of 1 for $\gamma$ restores the ideal 1D velocity definition.

The implementation of boundary conditions is relatively challenging and complex for the developed model as they should be imposed inside the domain of interest on virtual interfaces defined by mathematical expressions (i.e. on the moving interface defined by the zero iso-contour of the level set equation). The penalty method was used to overcome this issue and define the desired boundary conditions on the moving corrosion front.

Fig. 3.2 demonstrates a schematic presentation of the boundary conditions and general considerations of each PDE of the biodegradation mathematical model. This figure is divided into 5 different parts, presenting the 5 PDEs of the model. The Mg block is depicted in the center, and the interface separates it from the surrounding medium. There is no specific boundary condition for the level set and film formation equations, but in comparison to the other 3 transport equations, it should be noted that diffusivity is not considered for $\mathrm{Mg}(\mathrm{OH})_{2}$, which is also reflected in Eq. 3.11. The level set function $\phi$ is defined in a way that is positive inside and negative outside the solid region. For the $\mathrm{Mg}^{2+}$ ions transport equation, a Dirichlet boundary condition is applied on the mathematical interface to make the concentration equal to the saturation concentration of $\mathrm{Mg}^{2+}$ ions, a value that was already used in Eq. 3.17. For the $\mathrm{Cl}^{-}$ and $\mathrm{OH}^{-}$ ions transport equations, a no-flux boundary condition is applied to the interface by making the diffusion coefficient equal to zero inside the Mg block, preventing ions to diffuse inside the solid material.

To simulate the developed mathematical model, which is comprised of Eqs. 3.10, 3.11, 3.12, 3.13, and 3.17, a combination of finite difference and finite element methods was used, leading to discrete forms of these equations, which were subsequently solved using appropriate linear solvers.

To discretize the temporal terms of the aforementioned parabolic PDEs, a first-order backward Euler finite difference scheme was used, whereas the spatial terms were converted to a weak form using a standard first-order finite element scheme. Then, the open-source PDE solver FreeFEM [112] was used to implement the weak form and obtain a linear system of equations for each PDE. The obtained linear systems were solved in parallel using the HYPRE preconditioner [78] and the GMRES solver [223] via the open-source high-performance computing (HPC) toolkit PETSc [25]. Additionally, to increase the efficiency of the computation and decrease the simulation execution time, the computational mesh was decomposed and distributed among available computing resources using the interface of HPDDM package in FreeFEM [140]. The details of this implementation are presented in Chapter 7. A simple iterative solver based on the Newton method was also developed to solve Eq. 3.2 to obtain the value of $\alpha$ parameter if it was required in the simulations.

The computational mesh was generated using a set of first-order tetrahedral elements and was adaptively refined on the metal-solution interface to increase the numerical accuracy of the simulation of the level set equation (Eq. 3.17). The Netgen mesh engine [233] in the SALOME platform [217] was used to generate the mesh.

Similar to the technique employed by Bajger et al. [23], the gradient of concentration of $\mathrm{Mg}^{2+}$ in Eq. 3.17 was calculated at a distance $h$ in the normal direction from the interface, with $h$ being the smallest element size of the mesh:

Considering the adaptively refined mesh, the $h$ value is very small, so the gradient is computed at the regions close enough to the interface. In addition to this technique, the mass lumping feature of FreeFEM was used to prevent the oscillation of concentration values on the diffusive metal-medium interface.

The degradation rate of CP Mg was evaluated based on the hydrogen evolution tests performed either in NaCl or SBF solutions with eudiometers. The composition of the electrolytes is shown in the following table (Table 3.1). 0.5 g metallic chips (with a surface area of $47.7\pm 5.0\mathrm{cm}^{2}/\mathrm{g}$ and chip thickness ca. 200 microns) of CP Mg were put in 500 ml electrolyte for 22-24 hours for monitoring the amount of evolved hydrogen. The method of measuring evolved hydrogen was chosen for monitoring the degradation rate because although such a measurement is prone to experimental errors such as relatively high solubility of hydrogen in water and volume change due to temperature and pressure variations, it provides a continuous assessment of the process, resulting in a continuous and smooth curve. Additionally, as small metallic chips were used for the tests, it was not possible to clean these pieces in chromic acid without losing them to measure the mass loss directly. The drawback of choosing the evolved hydrogen as the monitoring method is that it is not the only occurring reaction since oxygen reduction also takes place during the process [265, 248, 240]. As a result, measuring only hydrogen does not capture the totality of the degradation reactions. However, for CP Mg, the contribution of oxygen reduction is low (in contrast to high-purity Mg [265]) and can simply be ignored, meaning that the evolved hydrogen is an accurate equivalent for the mass loss. The bulk pH of electrolytes before and after corrosion was measured by a pH meter (Metrohm-691, Switzerland). Local pH was measured by positioning pH microprobes (Unisense, Denmark, pH-sensitive tip size 10x50 micron) 50 micron above the surface of Mg and monitoring the pH values either in one spot or by horizontal or vertical line-scans or mapping by following a horizontal grid. The electrolytes were not pH buffered additionally since SBF contains carbonates and phosphates that stabilize the pH at the approximate value of 8.5 instead of the 10.5 characteristic for pure NaCl solutions where pH is stabilized by precipitation of $\mathrm{Mg}(\mathrm{OH})_{2}$. Meanwhile, synthetic pH buffers, such as TRIS and HEPES were proven to affect the degradation mechanism rather significantly and should not be used for this purpose [193]. The measurements were performed at room temperature of $22\pm 2^{\circ}\mathrm{C}$ maintained by the laboratory climate control system. More detailed information about experimental set up and procedures can be found elsewhere [192, 191].

The constructed mathematical model contains some parameters that need to be calibrated prior to final validation of the model: diffusion coefficient of $\mathrm{Mg}^{2+}$ and $\mathrm{Cl}^{-}$ ions ($D_{\mathrm{Mg}}$ and $D_{\mathrm{Cl}}$ to be inserted into Eq. 3.15 to get effective diffusion coefficients), the reaction rates of Eqs. 3.5 and 3.6 ($k_{1}$ and $k_{2}$), the momentum parameter, $\beta$, for controlling the saturation term behavior (in Eqs. 3.10, 3.11, and 3.15), and the $\gamma$ parameter for the initial interface velocity (Eq. 3.18). An inverse problem setup was required to estimate the proper value of these parameters.

Performing a parameter estimation requires running the computational models several times. Considering the computationally-intensive model of the current study, a sensitivity analysis was performed prior to the parameter estimation to exclude non-essential parameters and reduce the time required to complete the inverse problem run. This sensitivity analysis was accomplished separately in the low diffusion (similar to the SBF solution) and high diffusion (similar to NaCl solution) regimes.

After determining the essential parameters to include, a Bayesian optimization approach [196] was used to construct the inverse problem and calibrate the parameters. The reason for choosing a Bayesian approach was to minimize the number of optimization iterations, in each of which the simulation should run once. The Bayesian optimization is a more efficient option for such computational intensive cases in comparison to gradient-based or fully-stochastic methods as it takes into account all the preceding iterations in a probability tree [190].

The objective function of the optimization problem was the difference between the predicted and experimentally obtained values of evolved hydrogen. In the computational model, the evolved hydrogen can be computed directly at any time through the mass loss as each mole of corroded Mg is correlated to one mole of released hydrogen (Eq. 3.5). The mass loss can be obtained using the following volume integral:

where $\Omega_{+}(t)=\{\mathbf{x}:\phi(\mathbf{x},t)\geq 0\}$, and then, the amount of produced hydrogen is calculated using the ideal gas law:

in which $R$, $P$, $T$, $\mathrm{Mg}_{\mathrm{mol}}$ are the universal gas constant, the pressure, the medium temperature, and the molar mass of Mg, respectively.

In order to simulate the developed mathematical model, the experimental setup was reconstructed in silico with some minor differences. As there is no perfusion in the solution chamber, the mixing effect was neglected, so, as can also be seen in the mathematical model, the advection terms were not considered. Furthermore, the experiments were conducted using small metallic chips, yet, as the biodegradation behavior heavily depends on the exposed surfaces, we represented these chips by a cuboid with the same surface-to-mass ratio. By considering the approximate surface-to-mass of $50\mathrm{cm}^{2}/\mathrm{g}$ and the total mass of $0.5\mathrm{g}$, the chips were replaced by a cuboid with the size of $60\mathrm{mm}\times 21\mathrm{mm}\times 0.2\mathrm{mm}$, which approximately has the same ratio, surface area, volume, and mass. Also, the solution chamber with a capacity of $500\mathrm{ml}$ was represented by a cubic container with an edge size of $80\mathrm{mm}$. Fig. 3.3 depicts the constructed geometry as well as the computational mesh generated to represent it. The mesh is refined on the interface and contains ${18049471}$ elements, resulting in ${3077227}$ degrees of freedom (DOF) for each PDE.

Simulations were carried out on the VSC (Flemish Supercomputer Center) supercomputer. Taking advantage of HPC techniques to parallelize the simulation is an inevitable aspect of such a computational-intensive model, so based on what described in the implementation section, the mesh was decomposed among 170 computing cores, i.e. ${24137}$ DOF per core (which includes the ghost nodes to satisfy the boundary condition in each sub-mesh). On the VSC supercomputer, we made use of 5 nodes, 36 cores each, each node holding CPUs with a clock speed of 2.6 GHz, with 960 GB of the total available memory.

The $\mathrm{OH}^{-}$ transport equation (Eq. 3.13) was not solved during the parameter calibration process. Afterwards, two full simulations (for the NaCl and SBF solutions) were conducted to calculate the pH changes based on the change of the concentration of $\mathrm{OH}^{-}$ ions in the medium. This acted as the validation of the numerical model because no calibration was performed on the output of this equation. The pH was calculated using Eq. 3.9, based on the solution of Eq. 3.13 and a reported value of $7.00e\times 10^{-5}\mathrm{cm}^{2}/\mathrm{s}$ ($25.2\mathrm{mm}^{2}/\mathrm{hour}$) for the diffusion coefficient of $\mathrm{OH}^{-}$ ions ($D_{\mathrm{OH}}$ to be used in Eq. 3.15) in aqueous solutions [159].

According to the experimental setup, the initial concentration of the $\mathrm{Mg}^{2+}$, $\mathrm{Cl}^{-}$, and $\mathrm{OH}^{-}$ ions were set to 0 (no $\mathrm{Mg}^{2+}$ ions at the beginning), $146\mathrm{mM}$ ($5.175\times 10^{-6}\mathrm{g}/\mathrm{mm}^{3}$), and $1\times 10^{-7}\mathrm{g}/\mathrm{mm}^{3}$, respectively. The porosity ($\epsilon$) and tortuosity ($\tau$) of the protective film were considered to be 0.55 and 1, respectively [250]. The saturation concentration $[\mathrm{Mg}]_{\mathrm{sat}}$ was set to the solubility of magnesium chloride in water, which is $134\times 10^{-6}\mathrm{g}/\mathrm{mm}^{3}$ at $25^{\circ}\mathrm{C}$ [93]. The density of Mg ($[\mathrm{Mg}]_{\mathrm{sol}}$) and $\mathrm{Mg(OH)}_{2}$ were set to $1735\times 10^{-6}\mathrm{g}/\mathrm{mm}^{3}$ and $2344\times 10^{-6}\mathrm{g}/\mathrm{mm}^{3}$, respectively [23]. A time step convergence study was performed to determine the implicit time step size. Based on the results, a time step with a size of $0.025$ hours was chosen. The overall simulated time is 22 hours in accordance with the experimental design of performed immersion tests.

To further investigate the predictions of the current model on more complex shapes, the biodegradation of a simple screw was studied in the SBF solution using the parameters obtained for the low diffusion regimes. Similar to the simulation of Mg cuboid, the mesh was refined on the metal-medium interface, and it consisted of ${1440439}$ elements with ${246580}$ DOFs for each PDE. All the simulation parameters and materials properties were identical to the simulation of biodegradation in the SBF solution, and the target was to simulate $42$ days ($1008$ hours) of the process. This was selected as a sufficiently long simulated time to observe the effects of biodegradation on larger time scales. The goal of this case study is to demonstrate the applicability of the developed model for any desired 3D shape with no geometrical limitations. As a result, although it was possible to consider a more complicated geometry for the screw (for example by considering threads around the cylindrical part of the screw or having a realistic geometry for the head), the screw geometry consists of basic 3D primitives, which are adequate for the mentioned purpose.

Based on the performed sensitivity analysis, two parameter sets were obtained for the high diffusion (in NaCl solution) and low diffusion (in SBF solution) simulations, respectively. These parameters are listed in Table 3.2. According to the results, the reaction rate of Eq. 3.5 ($k_{1}$), which demonstrates the rate of oxidation-reduction, has less contribution to the process in comparison to the rate of the weakening of the protective film ($k_{2}$). Because of this, the parameter $k_{1}$ was not selected for the parameter estimation. Also, the model was sensitive to the effect of parameter $k_{2}$ in different ranges of values and not on a specific point, and as a result, three constant values were chosen as the delegates of these ranges in the optimization process. The model was not sensitive to the diffusion rate of $\mathrm{Cl}^{-}$ ions, which was also expected because although $\mathrm{Cl}^{-}$ has an important role in the weakening of the partially protective MgO film, its transport equation (Eq. 3.12) is purely diffusive and does not include any reaction term.

The parameter optimization process was performed on the specified range of selected parameters, while the rest of the parameter values were obtained from the literature [23, 159]. Table 3.3 shows the output of this process, which was used to simulate the full model. For two estimation processes, 120 optimization iterations were taken cumulatively, which took 276 hours of simulation execution time using 170 computing nodes for each simulation.

Fig. 3.4 shows the model output for the predicted produced hydrogen, protective film formation, and the pH changes. The graph of the evolved hydrogen is used as input during the parameter optimization process, but the pH results are produced by the simulations using the optimized parameters to demonstrate the validation of the developed mathematical and computational models. The predicted pH result (Fig. 3.4-d) shows a difference of $5.35\%$ for the simulation in NaCl and $1.03\%$ for SBF simulation. Each simulation took about 3 hours to complete.

In Fig. 3.4, a post-processed view of the final shape of the Mg cuboid in the NaCl solution is presented, in which the degraded geometry is plotted on the $\mathrm{Mg}^{2+}$ ions (Fig. 3.4-b) and protective film concentration (Fig. 3.4-c) contours. A transparent contour of the pH values in the solution is depicted for both the NaCl (Fig. 3.4-e) and SBF (Fig. 3.4-f) solutions. The range of colors is kept equal for both contours to make it easy to compare the change of pH in both solutions.

The concentration values of the state variables of the derived transport PDEs ($\mathrm{Mg}^{2+}$, $\mathrm{Cl}^{-}$, $\mathrm{OH}^{-}$, and $\mathrm{Mg(OH)}_{2}$) are plotted along a diagonal line in the solution container in Fig. 3.5, showing how they change in the zones close to the corrosion surface and far from it.

The simulation of 42 days (${19200}$ time steps) of the degradation of the simple screw took 9 hours to run using 170 computing cores. Fig. 3.6 depicts the post-processed interface and $\mathrm{Mg}^{2+}$ ions release (similar to Fig. 3.4-b) as well as the mass loss during the degradation of the screw in the SBF solution. It is worth mentioning that the roughness observed on the surface of the screw geometry is related to the node-based visualization of the level-set function evolution and is not caused by either non-uniform corrosion or numerical error in the simulations.