Chapter 5 Extending the model: simulating local pH evolution

This chapter is based on a manuscript prepared to be submitted: M. Barzegari, C. Wang, S.V. Lamaka, G. Zavodszky, and L. Geris, “Interface-coupled multiphysics computational modeling of local pH changes during the biodegradation of magnesium biomaterials.”

This chapter is partially based on a manuscript prepared to be submitted:
M. Barzegari, C. Wang, S.V. Lamaka, G. Zavodszky, and L. Geris, “Interface-coupled multiphysics computational modeling of local pH changes during the biodegradation of magnesium biomaterials.”

In this chapter, the basic biodegradation model developed in Chapter 3 is further developed to capture more advanced chemistry occurring at the corrosion interface in HBSS solutions, making it possible to simulate more complex environments and local pH changes during the biodegradation process.

5.1 Introduction

Mg is the most studied biodegradable metal [172, 299, 55, 295], on which many research groups have performed valuable biodegradation studies [77, 166, 19, 148]. The biodegradation behavior of Mg is investigated in in vitro corrosion tests, in which the selection of the corrosive media plays an important role since it affects the underlying chemical reactions [193]. By considering the main application of the biomaterial, which can be tissue engineering scaffolds, vascular stents, or orthopedic fixation implants, the corrosive media can be selected to be a representative of the service environment. The most basic form of the medium is a saline (NaCl) solution, in which the degradation rate is the highest possible [193]. More complex solutions can be used to mimic the behavior of the body environment by taking into account additional body fluid components, the most popular of which are Ringer’s solution, PBS (phosphate buffered saline), SBFs (simulated body fluids), HBSS (Hank’s balanced salt solution), and Earle’s balanced salt solution (EBSS) [193]. Adding more organic components to the solution will make it ready to simulate cell culture conditions. The common media for this purpose are MEM (Minimum Essential medium) and DMEM (Dulbecco’s modified Eagle’s medium) [193].

Various studies have already investigated the effect of different components in the aforementioned corrosive media on the degradation behavior of Mg materials [192, 292, 137, 154, 191]. A typical composition of ”simulated body fluid” solutions (such as SBF, HBSS, and EBSS) is chloride, carbonate, phosphates, sulfate, and calcium. The individual effect of these components on the rate of degradation of Mg has been extensively studied, in which it has been observed that carbonate and phosphate slow down the rate while the effect of sulfate is negligible [137, 191]. The concentration of $\mathrm{HCO}_{3}^{-}$ affects the pH buffering capacity and the degradation rate of Mg simultaneously [282].

The effect of calcium ions is more complex because it has been found that $\mathrm{Ca}^{2+}$ doesn’t contribute to the Mg corrosion directly [276], but a mixed effect of $\mathrm{Ca}^{2+}$ , $\mathrm{Mg}^{2+}$ , $\mathrm{HCO}_{3}^{-}$ , and $\mathrm{H}_{2}\mathrm{PO}_{4}^{-}/\mathrm{HPO}_{4}^{2-}$ forms a co-precipitation layer on the corroded surface of Mg, slowing down the corrosion rate of commercially pure Mg as well as of some Mg alloys [192, 154]. It has also been reported that although various Mg alloys show different intrinsic degradation behavior in NaCl solution, they possibly behave similarly in simulated body fluids [6, 191]. Since the humoral regulations inside the human body control the changes in pH of body fluids, it is common to use pH buffers to mimic a similar condition, but it should be taken into account that buffering solutions may affect the degradation rate of Mg [61, 142] and may also delay the formation of the precipitate layer [154]. An alternative solution to address this issue is to use natural pH buffers such as $\mathrm{HCO}_{3}^{-}/\mathrm{CO}_{2}$ , an option that is commonly used for immersion tests under cell culture conditions. In this situation, an equilibrium between $\mathrm{H}_{2}\mathrm{CO}_{3}\left(\mathrm{CO}_{2}\right)$ , $\mathrm{HCO}_{3}{}^{-}$ , and $\mathrm{CO}_{3}{}^{2-}$ keeps the pH constant. As a result, using simulated body fluids for corrosion tests without additional synthetic pH buffers is still acceptable [154, 191].

The major reactions occurring in simulated body fluids can be written as:

\mathrm{Mg}+2\mathrm{H}_{2}\mathrm{O}\rightarrow\mathrm{Mg}(\mathrm{OH})_{2}+% \mathrm{H}_{2}

(5.1)

2\mathrm{Mg}+2\mathrm{H}_{2}\mathrm{O}+\mathrm{O}_{2}\rightarrow 2\mathrm{Mg}(% \mathrm{OH})_{2}

(5.2)

5\mathrm{Ca}^{2+}+3\mathrm{PO}_{4}^{3-}+\mathrm{OH}^{-}\rightarrow\mathrm{Ca}_% {5}\left(\mathrm{PO}_{4}\right)_{3}\mathrm{OH}

(5.3)

\mathrm{Ca}^{2+}+2\mathrm{OH}^{-}\rightarrow\mathrm{Ca}(\mathrm{OH})_{2}

(5.4)

\mathrm{Mg}^{2+}+\mathrm{CO}_{3}^{2-}\rightarrow\mathrm{MgCO}_{3}

(5.5)

\mathrm{Ca}^{2+}+\mathrm{CO}_{3}^{2-}\rightarrow\mathrm{CaCO}_{3}

(5.6)

3\mathrm{Mg}^{2+}+2\mathrm{PO}_{4}^{3-}\rightarrow\mathrm{Mg}_{3}\left(\mathrm% {PO}_{4}\right)_{2}

(5.7)

\mathrm{CO}_{3}^{2-}(aq)+H^{+}(aq)\rightarrow\mathrm{HCO}_{3}^{-}(aq)

(5.8)

PO_{4}^{3-}(aq)+H^{+}(aq)\rightarrow\mathrm{HPO}_{4}^{2-}(aq)

(5.9)

The protection layer is formed on the corroded surface as a hydroxyapatite-like precipitation according to the following reaction [20, 243, 240, 136]:

	$\displaystyle m\mathrm{Mg}^{2+}+n\mathrm{Ca}^{2+}$	$\displaystyle+x\mathrm{H}_{2}\mathrm{PO}_{4}^{-}/\mathrm{HPO}_{4}^{2-}+y% \mathrm{HCO}_{3}^{-}+z\mathrm{OH}^{-}$		(5.10)
		$\displaystyle\rightarrow\mathrm{Mg}_{m}\mathrm{Ca}_{n}\left(\mathrm{PO}_{4}% \right)_{x}\left(\mathrm{CO}_{3}\right)_{y}\left(\mathrm{OH}^{-}\right)_{z}$		(5.10)

In fact, the similarity in corrosion behavior of various Mg alloys in SBF solutions is due to the similar composition of this quasi-protective layer, a mechanism that doesn’t occur in NaCl solution, leading to a more apparent difference in degradation rate between Mg alloys. The composition of the formed hydroxyapatite-like precipitation layer is close to the ones found in vivo [193]. Additionally, local pH measurements in HBSS and SBF show that, in contrast to saline solutions, the local pH value is not alkaline [154, 194]. This has been reported for hydrodynamics situations, under which the medium composition is kept constant by replacing the consumed ions by means of fluid flow.

Building mathematical and computational models of the biodegradation process in complex buffered solutions can help save the resources required to perform in vitro tests, but the details of the aforementioned chemistry are challenging to capture in a mechanistic model. Few attempts have been made to model the underlying chemical reactions in SBF solutions [119, 69, 290], but due to the complexity of the resulting mathematical models, it is not feasible to extend them to 3D and real-world cases, like for simulation of the degradation of biomedical implants. In the case of data-driven approaches [291], the applicability is limited to the studied conditions, making it difficult for developed models to achieve high extensibility and generalizability.

In this chapter, a detailed mathematical model is presented to extend the work discussed in Chapter 3 and [31], in which a mechanistic model of the biodegradation process is coupled with a thermodynamics-based code to predict local interfacial biodegradation of Mg in HBSS solutions. The local pH changes are the validation criterion to compare the simulation results with experimental ones. Besides other parameters affecting the Mg biodegradation mechanism, monitoring the pH changes at the degradation interface has proved to be significant due to its direct effect on the formation and stability of the degradation products layer [92].

5.2 Methods

5.2.1 Experimental setup

In this study, the corrosion tests were performed on ultra-high pure (UHP) and commercially pure (CP) Mg in hydrodynamics conditions. The elemental composition of the used materials is shown in Table 5.1, measured by Atomic Absorption Spectrometer (AAS). The samples were prepared as rod shapes with a diameter of $2\text{mm}$ mounted in epoxy resin with a disc shape. The electrolyte used for corrosion tests was commercial HBSS (ThermoFisher Scientific, no. 14025100), the composition of which is presented in Table 5.2.

Table 5.1: The elemental composition of ultra-high pure and commercial pure Mg used for performing corrosion experiments (percentages)

	Fe	Si	Mn	Al	Cu	Ni
CP-Mg	0.03420	0.0001	0.00237	0.00402	0.00037	¡0.0002
UHP-Mg	0.0012	¡0.0001	0.00037	0.00291	¡0.0001	¡0.0002

Table 5.2: Chemical composition of the HBSS electrolytes used to perform corrosion tests for local pH measurements.

Components	mM
$\mathrm{KCl}$	5.33
$\mathrm{KH}_{2}\mathrm{PO}_{4}$	0.44
$\mathrm{NaHCO}_{3}$	4.17
$\mathrm{NaCl}$	137.9
$\mathrm{Na}_{2}\mathrm{HPO}_{4}$	0.34
$\mathrm{CaCl}_{2}$	1.26
$\mathrm{MgCl}_{2}\cdot 6\mathrm{H}_{2}\mathrm{O}$	0.49
$\mathrm{MgSO}_{4}\cdot 7\mathrm{H}_{2}\mathrm{O}$	0.41
$\mathrm{D}-\mathrm{Glucose}$	5.56

Local pH measurements were performed using a glass-type pH microelectrode with a tip diameter of $10\mu\text{m}$ (Unisense, pH-10). The electrode was positioned above the sample at a distance of $50\pm 3\mu\text{m}$ . At this distance, a line scan mapping routine was performed to obtain the horizontal pH profiles, in which a specimen-centered area ( $3000\times 3000\mu\text{m}$ , which was used as the simulation domain as well) was mapped with the step length of $100\mu\text{m}$ . The sampling interval was 3 s, the result of which was a total duration of approximately 1 hour to scan the area. Similarly, the vertical pH profiles were scanned starting at the height of $50\mu\text{m}$ above the midpoint of the specimen up to $500\mu\text{m}$ in the bulk medium. Fig. 5.1 demonstrates the setup of the experiment schematically, in which the flow enters the chamber from the left with a rate of $1.0\text{mL.min}^{-1}$ (Medorex TL15E peristaltic pump) and leaves it to the right. The tests were performed independently at room temperature (RT, $22\pm 2$ in an air-conditioned lab) and $37^{\circ}\text{C}$ (Thermo Scientific SAHARA S13).

Figure 5.1: Experimental setup for validating the coupled biodegradation model

The in vitro cross-section morphology of the specimens was characterized using a dual beam FIB/SEM platform (LYRA3 TESCAN) equipped with an EDX system (Oxford Inca with a silicon drift detector), the result of which was used to compare qualitatively with the simulation predictions.

5.2.2 Computational model construction

The computational model in the current study comprises three coupled modules:

1.

An extended version of the mechanistic biodegradation model based on our previous work [31] (Chapter 3) to obtain Mg and hydroxide ions distributions and the initial formation of the protective layer based on the modeling hypotheses.
2.

A thermodynamics-based simulation to estimate the concentration of various components of the electrolyte in regions close to the surface of the sample based on the calculated pH of module 1.
3.

A module to link the former two, calculating the hydroxyapatite-like precipitation concentration, in which calculated pH values were transferred from module 1 to module 2 for each node of the computational mesh to calculate the concentration of ions depending on the computed local pH and transfer back the calculated precipitation concentration to module 1.

Fig. 5.2 shows a schematic representation of the coupled modules and the way that simulation data are being transferred between them. It should be noted that the results obtained from module 2 are calculated in an equilibrium state as it is a thermodynamics-based code, but module 1 works based on kinetics equations and calculations. As a result, it is assumed that the kinetics of the reactions is affected by the equilibrium state of involved chemical components in each time step, which can be considered as a small time window in which the chemical species are in equilibrium. This fact should be seen as an assumption made for simplifying the calculation of the stoichiometry of the protective layer (Eq. 5.10), in which the effect of the underlying kinetics on the change of pH is neglected.

Figure 5.2: Schematic representation of the coupled modules for modeling of local pH changes, showing how they are connected as well as the programming languages and tools used for the implementation.

The computational biodegradation model (module 1) was developed by deriving a set of reaction-diffusion-advection equations from the chemistry of the corrosion of Mg in hydrodynamics conditions. The following basic reactions are captured by this model, which is an extension of our previous contribution [31] (Chapter 3) as the overall process is described in more detailed equations:

\mathrm{Mg}+2\mathrm{H}_{2}\mathrm{O}\stackrel{{\scriptstyle k_{0}}}{{% \rightarrow}}\mathrm{Mg}^{2+}+\mathrm{H}_{2}+2\mathrm{OH}^{-}

(5.11)

\mathrm{Mg}^{2+}+2\mathrm{OH}^{-}\stackrel{{\scriptstyle k_{1}}}{{\rightarrow}% }\mathrm{Mg}(\mathrm{OH})_{2}+\mathrm{H}_{2}

(5.12)

\mathrm{Mg}(\mathrm{OH})_{2}+2\mathrm{Cl}^{-}\stackrel{{\scriptstyle k_{2}}}{{% \rightarrow}}\mathrm{Mg}^{2+}+2\mathrm{Cl}^{-}+2\mathrm{OH}^{-}

(5.13)

\mathrm{MgO}+\mathrm{Cl}^{-}+\mathrm{H}_{2}\mathrm{O}\stackrel{{\scriptstyle k% _{2}}}{{\rightarrow}}\mathrm{Mg}^{2+}+\mathrm{Cl}^{-}+2\mathrm{OH}^{-}

(5.14)

The following state variables hold the concentration of various basic ions involved in reactions described by Eqs. 5.11, 5.12, 5.13, and 5.14:

\begin{aligned} &C_{\mathrm{Mg}}=C_{\mathrm{Mg}}(\mathbf{x},t),\quad C_{% \mathrm{Mg}_{\mathrm{[s]}}}=C_{\mathrm{Mg}_{\mathrm{[s]}}}(\mathbf{x},t),\quad C% _{\mathrm{Mg}(\mathrm{OH})_{2}}=C_{\mathrm{Mg}(\mathrm{OH})_{2}}(\mathbf{x},t)% \\ &C_{\mathrm{Cl}}=C_{\mathrm{Cl}}(\mathbf{x},t),\quad C_{\mathrm{OH}}=C_{% \mathrm{OH}}(\mathbf{x},t)\quad\mathbf{x}\in\Omega\subset\mathbb{R}^{3}\end{% aligned},

(5.15)

Additionally, two more state variables are needed to couple the models, representing the calculated concentration of the hydroxyapatite-like precipitation as well as the cumulative layer concentration:

C_{\mathrm{Hydrox}}=C_{\mathrm{Hydrox}}(\mathbf{x},t),\quad C_{\mathrm{Film}}=% C_{\mathrm{Film}}(\mathbf{x},t)\quad\mathbf{x}\in\Omega\subset\mathbb{R}^{3},

(5.16)

where the total film concentration can be calculated as:

C_{\mathrm{Film}}=C_{\mathrm{Hydrox}}+C_{\mathrm{\mathrm{Mg}(\mathrm{OH})_{2}}}.

(5.17)

With the above state variables defined, the biodegradation model can be constructed by implementing the following PDEs:

\frac{\partial C_{\mathrm{Mg}_{\mathrm{[s]}}}}{\partial t}=-k_{0}C_{\mathrm{Mg% }_{\mathrm{[s]}}}

(5.18)

\frac{\partial C_{\mathrm{Mg}}}{\partial t}=\nabla\cdot\left(D_{\mathrm{Mg}}^{% e}\nabla C_{\mathrm{Mg}}\right)-\nabla\cdot\left({\mathbf{u}}C_{\mathrm{Mg}}% \right)+k_{0}C_{\mathrm{Mg}_{\mathrm{[s]}}}-k_{1}\alpha C_{\mathrm{Mg}}C_{% \mathrm{OH}}^{2}+k_{2}C_{\mathrm{Film}}{C_{\mathrm{Cl}}}^{2}

(5.19)

\frac{\partial C_{\mathrm{Mg}(\mathrm{OH})_{2}}}{\partial t}=k_{1}\alpha C_{% \mathrm{Mg}}C_{\mathrm{OH}}^{2}-k_{2}C_{\mathrm{Film}}{C_{\mathrm{Cl}}}^{2}

(5.20)

\frac{\partial C_{\mathrm{Cl}}}{\partial t}=\nabla\cdot\left(D_{\mathrm{Cl}}^{% e}\nabla C_{\mathrm{Cl}}\right)-\nabla\cdot\left({\mathbf{u}}C_{\mathrm{Cl}}\right)

(5.21)

\frac{\partial C_{\mathrm{OH}}}{\partial t}=\nabla\cdot\left(D_{\mathrm{OH}}^{% e}\nabla C_{\mathrm{OH}}\right)-\nabla\cdot\left({\mathbf{u}}C_{\mathrm{OH}}% \right)+k_{0}C_{\mathrm{Mg}_{\mathrm{[s]}}}-k_{1}\alpha C_{\mathrm{Mg}}C_{% \mathrm{OH}}^{2}+k_{2}C_{\mathrm{Film}}{C_{\mathrm{Cl}}}^{2}

(5.22)

showing the mathematical representation of reactions 5.11, 5.12, 5.13, and 5.14 in hydrodynamics condition in the form of a set of reaction-diffusion-advection equations. $k_{0}$ , $k_{1}$ , and $k_{2}$ are reaction rate constants corresponding to the Mg oxidation, film formation, and film elimination reactions, respectively. In Eqs. 5.19, 5.20, and 5.22, $\alpha$ is defined as:

\alpha=\left(1-\beta\frac{C_{\mathrm{Film}}}{[\mathrm{Film}]_{\max}}\right),

(5.23)

in which the protective film maximum concentration is calculated using its porosity ( $\epsilon$ ) [23]:

[\mathrm{Film}]_{\max}=\rho_{\mathrm{Film}}\times(1-\epsilon).

(5.24)

In Eqs. 5.19, 5.21, and 5.22, $\mathbf{u}$ is the velocity field from the surrounding fluid flow governed by the Stokes equation:

\left\{{\begin{array}[]{*{20}{l}}\displaystyle{-\nu\Delta\mathbf{u}+\nabla p+% \frac{\nu}{K}\mathbf{u}=0}\\ \displaystyle{\nabla\cdot\mathbf{u}+\varepsilon p=0.}\end{array}}\right.

(5.25)

in which $\mathbf{u}$ is the fluid velocity, $p$ is the pressure (which is actually pressure divided by the density), $\nu=\frac{\mu}{\rho}$ is the kinematic viscosity (with $\mu$ being the dynamic viscosity), and $K$ is the permeability function.

The local pH can be calculated using the simulated concentration of hydroxide (Eq. 5.22):

\mathrm{pH}=14+\log_{10}\left(C_{\mathrm{OH}}/\mathrm{MW}_{\mathrm{OH}}\times 1% 0^{6}\right),

(5.26)

with $\mathrm{MW}_{\mathrm{OH}}=17.01$ being the molecular weight (MW) of the hydroxide ions.

The concentration of the hydroxyapatite-like precipitation ( $C_{\mathrm{Hydrox}}$ in Eqs. 5.15 and 5.17) is calculated using the thermodynamics-based simulations (module 2) based on calculated local pH (Eq. 5.26) for each node of the desired domain. After solving the derived equations in each time step, the linking module passes the obtained local pH to the thermodynamics module to calculate the individual concentration of involved chemical components. Then, the individual concentrations are converted to the concentration of the hydroxyapatite-like precipitation by taking into account the stoichiometry of the formation reaction (Eq. 5.10), leading to the calculation of $C_{\mathrm{Hydrox}}$ for each node. After this, the total concentration of the film can be calculated according to Eq. 5.17 by passing back the calculated value to module 1.

The thermodynamics-based simulations were conducted using the Hydra-Medusa code [130, 271, 76], in which the input data of chemical equilibrium constants, solubility products, temperature, and involved chemical components are used to generate a set of equilibrium diagrams correlating pH to concentration or fraction of desired chemical components. The experimental conditions, including the initial composition of the electrolyte (Table 5.2) and evaluated temperatures ( $25^{\circ}C$ and $37^{\circ}C$ ), were given as input, and contributing components and solubility products were selected as represented in Table 5.3.

Table 5.3: Solubility products of related chemical reactions at

25^{\circ}C

(RT) and

37^{\circ}C

[266]

Chemical reaction	$pK_{sp}25^{\circ}C$	$pK_{sp}37^{\circ}C$
$\mathrm{Ca}_{5}\left(\mathrm{PO}_{4}\right)_{3}\mathrm{OH}\rightarrow 5\mathrm% {Ca}^{2+}+3\mathrm{PO}_{4}^{3-}+\mathrm{OH}^{-}$	54.46	58.77
$\mathrm{Mg}(\mathrm{OH})_{2}\rightarrow\mathrm{Mg}^{2}+2\mathrm{OH}^{-}$	11.25	11.25
$\mathrm{Ca}(\mathrm{OH})_{2}\rightarrow\mathrm{Ca}^{2+}+2\mathrm{OH}^{-}$	5.20	5.38
$\mathrm{MgCO}_{3}\rightarrow\mathrm{Mg}^{2+}+\mathrm{CO}_{3}^{2-}$	8.03	5.51
$\mathrm{CaCO}_{3}\rightarrow\mathrm{Ca}^{2+}+\mathrm{CO}_{3}^{2-}$	8.48	8.44
$\mathrm{Mg}_{3}\left(\mathrm{PO}_{4}\right)_{2}\rightarrow 3\mathrm{Mg}^{2+}+2% \mathrm{PO}_{4}^{3-}$	23.28	27.62
$\mathrm{H}_{2}\mathrm{O}(\mathrm{l})\rightarrow\mathrm{H}^{+}(aq)+\mathrm{OH}^% {-}(aq)$	14.00	13.61
$\mathrm{HCO}^{3-}(aq)\rightarrow\mathrm{CO}_{3}^{2-}(aq)+H^{+}(aq)$	10.33	10.24
$\mathrm{HPO}_{4}^{2-}(aq)\rightarrow\mathrm{PO}_{4}^{3-}(aq)+H^{+}(aq)$	12.35	12.32

The derived PDEs for the mechanistic model (Eqs. 5.19, 5.20, 5.21, and 5.22) were solved using a standard first-order finite element scheme. The open-source PDE solver FreeFEM [112] was used to implement the finite element model, resulting in a linear system of equations. The obtained equations were solved in parallel using efficient preconditioners and iterative solvers available in the open-source high-performance computing (HPC) toolkit PETSc [25]. The HYPRE BoomerAMG [78] and FieldSplit preconditioners were applied to the reaction-diffusion PDEs and the Stokes equations, respectively, and the GMRES solver [223] was used to solve the linear systems. Moreover, the computational mesh was partitioned and distributed among available computing resources using the HPDDM preconditioner [140]. Additionally, a Level-Set-based approach was employed to track the change of morphology of the degrading part, on the solution of which appropriate boundary conditions for the PDEs were applied via the penalization method. The details of this implementation are presented in our published works [31, 30] and Chapters 3 and 7. Since Eq. 5.22 is a non-linear PDE, a Picard-relaxation approach was followed to linearize this equation. The linking module (module 3) was implemented as a FreeFEM plugin in C++.

5.2.3 Simulation setup

Since the objective of the current study is to investigate local pH changes in regions close to the degrading metal, the computational domain was selected to include only a small portion of the chamber used in the experimental setup (Fig. 5.1). The domain was a $3\times 3\times 3\text{mm}$ cube on top of the degrading object, and the degrading block, represented as a disc with a diameter of $2\text{mm}$ and height of $0.3\text{mm}$ , was attached to the outside of the cube. The cube size was selected such that it represents the sample-centered area used in the experimental setup for line scan mapping. This setup is depicted in Fig. 5.3. The computational mesh was generated using a set of first-order tetrahedral elements, which was adaptively refined on the interface of the degrading sample to increase the numerical accuracy of the employed interface capturing technique. The Netgen mesh engine [233] in the SALOME platform [217] was used to generate the mesh. The generated mesh contained ${290997}$ elements with ${51757}$ degrees of freedom (DOF) for each equation.

Figure 5.3: Computational model setup for local pH simulations: a) the geometry of the desired domain, being a small portion of the full chamber on top of the disc-like sample, b) the dimensions of the domain and applied boundary conditions for the fluid flow simulation, c) the generated computational mesh, depicting the flow region in green and the biodegradable sample in yellow, d) a cross-section of the mesh showing the refined meshing on the interface of the degradable metal.

The boundary conditions include an inlet on the left of the box, an outlet on the right, a wall on the bottom except where the degrading part exists, and a free slip on top. The inlet velocity was selected to be a linear profile ranging from zero on the bottom to $0.01\text{mm.s}^{-1}$ on the top, a value coming from solving the Navier-Stokes equations in the full chamber model (described in Chapter 4).

The values of model parameters were set based on our previous work [31] (Chapter 3), but the following assumptions were also applied for selecting the proper parameters of the coupled computational model:

1.

Diffusion rate of hydroxide ions is 5-10 times more than the rate of Mg ions, so $D_{\mathrm{OH}}$ was set as $D_{\mathrm{OH}}=7.5D_{\mathrm{Mg}}$ [92].
2.

The hydroxyapatite-like precipitation does not form immediately at the beginning and emerges later [92, 266]. In the current model, it is set to start forming after the first hour of simulation time.
3.

The magnesium hydroxide layer is ticker at $37^{\circ}C$ [266], so the film elimination parameter ( $k_{2}$ ) was set to be 50 times lower in this temperature in comparison to the room temperate.

The density of the electrolyte was selected to be $1.085\times 10^{-3}\text{g.mm}^{-3}$ , and the dynamic viscosity was set to $1.28\times 10^{-3}\text{g.mm}^{-1}\text{s}^{-1}$ . The simulations were carried out on the Snellius supercomputer using 100 CPU cores on a thin node with 256 GB of total memory.

5.3 Results

5.3.1 Thermodynamics-based simulation

Inputting the experimental conditions in the Hydra-Medusa software, including the initial composition of the electrolyte, temperature and the contributing chemical components, results in a big and complex output containing all the possible occurring chemical reactions. From this output, relevant reactions for the current biodegradation systems were filtered out, and the filtered reactions were converted to desired concentration units. Fig. 5.4 depicts the output of this process separately for simulations performed at $25^{\circ}C$ (RT) and $37^{\circ}C$ , showing how the concentration of relevant components varies with changing the environment pH. The solubility products of these simulations were taken from Table 5.3, and equilibrium concentrations were set to be equal to the electrolyte composition, listed in Table 5.2. These results (module 2) were used by the linking module (module 3) to provide equilibrium information for the mechanistic model (module 1).

Figure 5.4: Selected relevant components from the Hydra-Medusa software output for given experimental conditions, showing how the concentration of various components vary with changing local pH

5.3.2 Biodegradation simulations

In the current study, the local pH profiles, i.e., the pattern of pH changes in the region close to the degradation surface, were used to validate the developed coupled computational model. The comparison between computational predictions and experimental results was made in 3 different ways: 1) qualitative comparison of pH distribution above the degradable metal in a sample-centered square with an edge size of $3\text{mm}$ , 2) comparing the horizontal pH profiles, where the pH was measured over a line parallel to the sample, and 3) comparing vertical pH profiles, in which the pH was measured over a distance from the surface of the sample to the bulk of the electrolyte.

Fig. 5.5 shows a visualization of the local pH profiles from the top and side views for simulations performed at $25^{\circ}C$ (RT) and $37^{\circ}C$ after 12 hours. These patterns are comparable with the local pH distribution measured in the experiments, as shown in Fig. 5.6. In these figures, the flow is from left to right, advecting the released ions in the flow direction.

Figure 5.5: Simulation results for local pH predictions, depicting the local pH in a) top view from a horizontal cross-section, b) side view from a vertical cross-section, and c) perspective view with both the top and side cross sections, for simulations performed at $25^{\circ}C$ (RT) and $37^{\circ}C$ .

Figure 5.6: Experimental results for the distribution of local pH above the surface of the sample, measured after 12 hours of degradation.

Quantitative profiles provide a more accurate comparison between the computational predictions and experimentally-obtained values. The horizontal profiles, also called line scans, are depicted in Fig. 5.8, in which the local pH changes are plotted over a horizontal line located $50\mu\mathrm{m}$ above the surface of the sample. The profiles are shown separately for the experiments performed with CP and UHP Mg and the computational model. The results were recorded after 3, 6, and 12 hours of immersion.

Figure 5.7: Comparing computational predictions and experimentally-obtained horizontal line scans of local pH changes, measured $50\mu\mathrm{m}$ above the surface of the sample after 3, 6, and 12 hours of degradation.

Similarly, Fig. 5.8 shows the vertical pH profile after 12 hours of degradation, which was the final time of the simulation. The local pH is measured over a distance of $0.5\mathrm{mm}$ vertically, starting from $50\mu\mathrm{m}$ above the sample.

Figure 5.8: Comparing computational and experimental results for the local pH profiles, measured over a vertical line above the center of the sample after 12 hours of immersion.

5.4 Discussion

In this study, an interface-coupled multiphysics biodegradation model was developed in order to predict the corrosion behavior in simulated body fluids in the presence of various chemical components. The quantity to compare with the experimental results was the predicted local pH, which reflects the capabilities of the model to capture the complex chemistry occurring near the biodegradation interface. For this end, our previous biodegradation model [31] (Chapter 3) was extended and coupled with a fluid flow model capturing the hydrodynamics condition. Then, the coupled model was linked to a thermodynamics-based code for computing the concentration of involved chemical components, the results of which were provided back to the biodegradation model via a linking module to compute the concentration of the precipitation layer.

It has been shown that the local pH changes can be a reflective characteristic of the biodegradation process in SBF-like solutions [92, 266]. Consequently, the vertical and horizontal local pH profiles can be reasonably used for the validation of the computational models of the biodegradation process. Capturing the complex interaction of various chemical components on the biodegradation interface in a mechanistic model can be a big challenge, especially for 3D cases with any arbitrary shape. That’s why several modules were coupled in the current study to deliver such a model.

In this work, while the experiments were performed using CP and UHP Mg, the computational model does not take into account the difference in the elemental composition of these materials. Instead, the biodegradation model was developed by ignoring the effect of alloying elements and impurities. There are noticeable quantitative and qualitative variations in the obtained experimental results, such as different behavior of CP and UHP for horizontal pH profiles in line scan mappings (Fig. 5.7) and distribution of local pH above the sample (Fig. 5.6). However, these differences are still roughly the same behavior and quantity, which support the performed assumption and simplification of ignoring the dissimilarity of elemental composition. As a consequence, the computational results, including the pH profiles (Fig. 5.8), line scans (Fig. 5.7), and local pH distributions (Fig. 5.5), lies between the values and profiles obtained for CP and UHP Mg.

Recent studies on measuring local pH changes on the biodegradation surface of Mg alloys show that the local values are different from the pH within the bulk of the electrolyte [92]. For example, in the tests performed by Mareci et al. [185], such observation was made when the vertical pH profiles were compared to the global pH values. Similar behavior was observed in the Wang et al. work [266]. Reproducing this behavior is problematic in mechanistic computational models due to the uniformity of the diffusion of hydroxide ions that change the pH. In other words, a uniform diffusion with a relatively high diffusion coefficient leads to the same local and global pH values. In the model presented here however, spatially dependent behavior was successfully reproduced in silico by letting a narrow film form at the beginning of the mechanistic biodegradation simulation and then computing the interaction of various ions in this narrow region by using the coupled thermodynamics-based code. Comparing the results of the predicted vertical pH profiles with the experimentally obtained curves (Fig. 5.8) shows good agreement, implying that the employed approach was able to mimic the complex chemistry on the biodegradation surface from a quantitative point of view. Yet, a different behavior can be observed further away from the surface where the local pH approaches the global value. In experimental results, this happens at a shorter distance from the implant compared to the computational predictions (Fig. 5.8) where yielding to the global pH seems to occur at longer distances from the biodegradation surface. This behavior can be confirmed by the visualization of computed pH (Fig. 5.5) and is due to the uniform diffusion mechanism of ions in the computational model, causing a gradient in the bulk part outside of the narrow formed region.

The results of horizontal line scan mapping on the surface of the sample and local pH distribution (Figs. 5.7 and 5.6) show that the pH starts to change in regions above the sample. In the computational predictions depicted in Fig. 5.7, the change of pH starts sharply when the line scan reaches the sample, a behavior rooted in the presence of the fluid flow model, which prevents ions from being diffused to the left, as shown in Fig. 5.5. Similarly, the advected ions in the direction of fluid flow (to the right in Fig. 5.5) prevent the pH value from decreasing dramatically where the sample ends in the line scan mapping. Moreover, the line scans show different behavior for experiments performed at room temperature compared to $37^{\circ}\text{C}$ , especially for CP Mg samples. At RT, the horizontal pH profile tends to descend slightly over time, meaning that the profile measured after 3 hours of degradation is placed above the one measured after 12 hours. This behavior occurs oppositely for the measurement done at $37^{\circ}\text{C}$ . The computational predictions of line scan profiles follow the same pattern, which could be due to how the coupled model computed the change of concentration of various chemical components and their contribution to the change of local pH.

Regarding the horizontal line scan results, a slight increase in the local pH values can be observed in the computational results in the region where the sample ends (Fig. 5.7). Although a similar behavior can be seen in the experimentally-obtained line scans (Fig. 5.6), it should be further investigated from the numerical implementation perspective to make sure that numerical artifacts do not contribute to it. One possible approach that may help remove numerical artifacts in this regard is to make the type of employed boundary conditions consistent in the employed fluid flow model. In the current model, a velocity condition was used for the inlet while the outlet condition was fixed pressure. Using the same type of boundary conditions for both inlet and outlet, meaning either assigning pressure conditions to the inlet or velocity to the outlet, can be considered in further extensions of the model to check their effect on the obtained local pH results.

As mentioned before, various chemical components of the electrolyte contribute differently to the change of pH, among which $\mathrm{Ca}^{2+}$ seems to have the most intricate effect [276]. It has been reported that in the absence of $\mathrm{Ca}^{2+}$ ions, the surface pH tends to be between 10 and 11, the typical range for pH value in biodegradation tests performed in saline solution [92]. This value is in line with the computational predictions of our previous study [31] (Chapter 3). However, in presence of $\mathrm{Ca}^{2+}$ cations, local pH between 7.8 and 8.5 and significantly lower degradation rates are reported [192, 90, 254, 153]. These findings are in line with the predictions made by the model presented in this chapter, obtained by coupling a mechanistic modeling approach and chemical equilibrium modeling, the latter of which considers the presence of $\mathrm{Ca}^{2+}$ .

The simplification assumptions made for developing the computational model, including the correlation between the diffusion rates of hydroxide and Mg ions, higher film elimination rate at room temperature, and the time at which the precipitation layer starts to form, can be seen as the limitations of the developed model. Regarding the latter mentioned assumption, the time was hard-coded into the model due to limitations of implementing the actual mechanisms of inducing the precipitation in the model. This can be improved in further developments of the model by incorporating a kinetics-based description of the precipitation layer formation, similar to the formation of magnesium hydroxide layer defined using a PDE (Eq. 5.20). This improvement will remove the necessity of the assumption made for linking the results obtained in equilibrium using the thermodynamics-based code (module 2) with the kinetics-based model (module 1), which was elaborated in section 5.2.2.

Knowing the mechanism of ionic activities at regions close to the biodegradation interface is crucial for comprehending the chemical process of Mg in complex solutions such as simulated physiological conditions. The developed computational model can be seen as an a facilitating tool for moving towards understanding these mechanisms by making it possible to virtually investigate the effective parameters such as the flow rate and environment composition. This knowledge will be helpful for biomedical applications where new coating and stabilization systems can be developed for different implantation environments.

5.5 Conclusion

In this chapter, the model developed in Chapter 3 was extended and combined with the fluid flow model developed in Chapter 4 to deliver a computational model of the biodegradation process in hydrodynamics conditions capable of predicting local pH changes close to the corrosion surface. In order to simulate local pH evolution in simulated body fluids, the precipitation of the hydroxyapatite-like protective model should be modeled. This was done by coupling the mentioned model with a thermodynamics-based code to compute the concentration of active chemical components contributing to the formation of the protective layer on the degradation interface. Results obtained from the coupled model show a good agreement with experimentally-obtained local pH measurements, demonstrating the effectiveness of the developed simulation workflow.