2.12 Chemical Reactions

Chemical reactions may be incorportated into a multiphasic mixture by adding a mass supply term to the equation of mass balance,

\[ \begin{equation} \frac{\partial\rho^{\alpha}}{\partial t}+\divg\left(\rho^{\alpha}\mathbf{v}^{\alpha}\right)=\hat{\rho}^{\alpha}\,,\label{eq150} \end{equation} \]

Where \(\hat{\rho}^{\alpha}\) is the volume density of mass supply to \(\alpha\) resulting from chemical reactions with all other mixture constitutents. Since mass must be conserved over all constituents, mass supply terms are constrained by

\[ \begin{equation} \sum\limits_{\alpha}\hat{\rho}^{\alpha}=0\,.\label{eq151} \end{equation} \]

In a mixture containing a solid constituent (denoted by \(\alpha=s\) ), it is conveniemt to define the mixture domain (and thus the finite element mesh) on the solid and evaluate mass fluxes of constituents relative to the solid,

\[ \begin{equation} \mathbf{m}^{\alpha}=\rho^{\alpha}\left(\mathbf{v}^{\alpha}-\mathbf{v}^{s}\right)\,.\label{eq152} \end{equation} \]

Substituting eq.\eqref{eq152} into eq.\eqref{eq150}, the differential form of the mass balance may be rewritten as

\[ \begin{equation} \frac{D^{s}\rho_{r}^{\alpha}}{Dt}+J\divg\mathbf{m}^{\alpha}=\hat{\rho}_{r}^{\alpha}\,,\label{eq153} \end{equation} \]

Where \(D^{s}\left(\cdot\right)/Dt\) represents the material time derivative in the spatial frame, following the solid, \(J=\det\mathbf{F}\), where \(\mathbf{F}\) is the deformation gradient of the solid matrix; \(\rho_{r}^{\alpha}\) is the apparent density and \(\hat{\rho}_{r}^{\alpha}\) is the volume density of mass supply to \(\alpha\) normalized to the mixture volume in the reference configuration,

\[ \begin{equation} \rho_{r}^{\alpha}=J\rho^{\alpha},\quad\hat{\rho}_{r}^{\alpha}=J\hat{\rho}^{\alpha}.\label{eq154} \end{equation} \]

Since \(\rho_{r}^{\alpha}\) is the mass of \(\alpha\) in the current configuration per volume of the mixture in the reference configuration (an invariant quantity), this parameter represents a direct measure of the mass content of \(\alpha\) in the mixture, which may thus be used as a state variable in a framework that accounts for chemical reactions. A distinction is now made between solid and solute species in the mixture, since they are often treated differential in an analysis.

Solid Matrix and Solid-Bound Molecular Constituents

For constituents constrained to move with the solid (denoted generically by \(\alpha=\sigma\) and satisfying \(\mathbf{v}^{s}=\mathbf{v}^{\sigma}\), \(\forall\sigma)\), the statement of mass balance in eq.\eqref{eq153} reduces to the special form

\[ \begin{equation} D^{s}\rho_{r}^{\sigma}/Dt=\hat{\rho}_{r}^{\sigma}\,.\label{eq155} \end{equation} \]

This representation makes it easy to see that alterations in \(\rho_{r}^{\sigma}\) can occur only as a result of chemical reactions (such as synthesis, degradation, or binding). In contrast, as seen in eq.\eqref{eq153}, alterations in \(\rho_{r}^{\alpha}\) for solutes or solvent (\(\alpha\ne\sigma)\) may also occur as a result of mass transport into or out of the pore space of the solid matrix. Therefore, \(\rho_{r}^{\sigma}\) is the natural choice of state variable for describing the content of solid constituents in a reactive mixture.

When multiple solid species are present, the net solid mass content may be given by \(\rho_{r}^{s}=\sum\limits_{\sigma}\rho_{r}^{\sigma}\) whereas the net mass supply of solid is \(\hat{\rho}_{r}^{s}=\sum\limits_{\sigma}\hat{\rho}_{r}^{\sigma}\) such that \(D^{s}\rho_{r}^{s}/Dt=\hat{\rho}_{r}^{s}\). The referential solid volume fraction, \(\varphi_{r}^{s}\), may be evaluated from

\[ \begin{equation} \varphi_{r}^{s}=\varphi_{0}^{s}+\sum\limits_{\sigma}\rho_{r}^{\sigma}/\rho_{T}^{\sigma},\label{eq:referential-solid-volume-fraction} \end{equation} \]

where \(\rho_{T}^{\sigma}\) is the true density of solid constituent \(\sigma\) (mass of \(\sigma\) per volume of \(\sigma\)) and \(\varphi_{0}^{s}\) is the referential solid volume fraction of solid constituents not explicitly modeled by solid-bound molecules (a user-defined parameter). According to eq.\eqref{eq154}, it follows that the solid volume fraction in the current configuration is given by \(\varphi^{s}=\varphi_{r}^{s}/J\). Note that \(0\leqslant\varphi^{s}\leqslant1\) under all circumstances, while \(0\leqslant\varphi_{r}^{s}\leqslant J\), implying that \(\varphi_{r}^{s}\) may exceed unity when solid growth occurs. In this study, it is assumed that all mixture constituents are intrinsically incompressible, implying that their true density is invariant.

At the start of an analysis, the formula of eq.\eqref{eq:referential-solid-volume-fraction} takes the form

\[ \begin{equation} \varphi_{r}^{s}\left(0\right)=\varphi_{0}^{s}+\sum\limits_{\sigma}\rho_{r}^{\sigma}\left(0\right)/\rho_{T}^{\sigma}\label{eq:initial-referential-solid-fraction} \end{equation} \]

where \(\rho_{r}^{\sigma}\left(0\right)\) represents the initial apparent mass densities of solid constituents. This initial value of \(\varphi_{r}^{s}\) must satisfy \(0\le\varphi_{r}^{s}\left(0\right)\le1\), since \(J=1\) at the start of an analysis.

The various constituents of the porous solid matrix of a multiphasic mixture may be electrically charged. The charge density in the current configuration, normalized by the mixture volume in the current configuration, is given by

\[ \begin{equation} \check{c}^{F}=\frac{z^{F}dn^{F}+\sum_{\sigma}z^{\sigma}dn^{\sigma}}{dV}\,,\label{eq:FCD-mixture-volume} \end{equation} \]

where \(z^{F}\) is the charge number (equivalent charge per mole) and \(dn^{F}\) is the elemental number of moles associated with fixed charges present in the initial volume fraction \(\varphi_{0}^{s}\) of the solid matrix. Similarly, \(z^{\sigma}\) is the charge number of evolving solid constituent \(\sigma\). By convention however, the fixed charge density \(c^{F}\) of a multiphasic material is normalized by the fluid volume of the mixture, thus it can be shown that

\[ \begin{equation} c^{F}=\frac{\check{c}^{F}}{1-\varphi^{s}}=\frac{1}{J-\varphi_{r}^{s}}\left(\frac{z^{F}dn^{F}}{dV_{r}}+\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\right)\,.\label{eq:FCD-fluid-volume} \end{equation} \]

In particular, in the reference configuration this formula produces

\[ \begin{equation} c^{F}\left(0\right)=\frac{1}{1-\varphi_{r}^{s}\left(0\right)}\left(\frac{z^{F}dn^{F}\left(0\right)}{dV_{r}}+\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}\left(0\right)}{M^{\sigma}}\right)\,.\label{eq:FCD-initial-time} \end{equation} \]

We define the first term on the right-hand-side of this equation as the initial value of the user-specified fixed charge density \(c_{0}^{F}\) (which is not associated with solid species \(\sigma\)),

\[ \begin{equation} c_{0}^{F}\equiv\frac{1}{1-\varphi_{r}^{s}\left(0\right)}\frac{z^{F}dn^{F}\left(0\right)}{dV_{r}}\,.\label{eq:user-specified-cF0} \end{equation} \]

Therefore, in the current configuration, we may rewrite

\[ \begin{equation} c^{F}=\frac{1-\varphi_{r}^{s}\left(0\right)}{J-\varphi_{r}^{s}}c_{0}^{F}+\frac{1}{J-\varphi_{r}^{s}}\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\,.\label{eq158} \end{equation} \]

where the user-specified \(c_{0}^{F}\) may optionally be associated with a load curve representing the scale factor \(dn^{F}/dn^{F}\left(0\right)\), in case the user would like to allow \(c_{0}^{F}\) to evolve. By analogy, we may now define the referential fixed-charge density \(c_{r}^{F}\) of the mixture as

\[ \begin{equation} c_{r}^{F}\equiv\frac{J-\varphi_{r}^{s}}{1-\varphi_{r}^{s}\left(0\right)}c^{F}=c_{0}^{F}+\frac{1}{1-\varphi_{r}^{s}\left(0\right)}\sum_{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\label{eq157} \end{equation} \]

Here, \(c_{r}^{F}\) may evolve with time, but it represents the number of fixed equivalent charges in the current configuration, per fluid volume in the reference configuration. Hence we may also write

\[ c^{F}=\frac{1-\varphi_{r}^{s}\left(0\right)}{J-\varphi_{r}^{s}}c_{r}^{F}\,. \]

The molar concentration of a solid-bound molecular constituent, which may be needed in a reactive process involving solutes, is given by

\[ \begin{equation} c^{\sigma}=\frac{1}{J-\varphi_{r}^{s}}\frac{\rho_{r}^{\sigma}}{M^{\sigma}}\,.\label{eq:sbm-molar-concentration} \end{equation} \]

An alternative to using solid-bound molecules is to define solutes \(\sigma\) whose diffusivity \(\mathbf{d}^{\sigma}\) in the mixture is set to \(\mathbf{0}\) (however, the free diffusivity \(d_{0}^{\sigma}\) should not be set to zero, to prevent a division by zero; its exact value is not important, thus let \(d_{0}^{\sigma}=1\)). Then, FEBio treats these solutes as equivalent to solid-bound molecules: (1) Their concentration does not contribute to the osmolarity of the interstitial fluid; (2) their concentration contributes to the fixed charge density \(c^{F}\) in the current configuration; (3) these solutes do not contribute to the effective hydraulic permeability \(\tilde{\mathbf{k}}\) of the porous multiphasic mixture; (4) these solutes can be involved in chemical reactions; (5) while their initial effective concentration may be prescribed, it must not contribute to the prescribed initial effective fluid pressure, and the user should not prescribe any boundary conditions on the effective concentration of these solutes. Using these 'solid-bound' solutes increases the number of degrees of freedom in an FEBio analysis; however, unlike solid-bound molecules, the solution for the effective concentration of these solutes remains as accurate as all other degrees of freedom in an analysis.

To evaluate the contribution of these solutes to the referential solid volume fraction \(\varphi_{r}^{s}\) and to chemical reactions, we use eq.\eqref{eq:sbm-molar-concentration} above and substitute it into eq.\eqref{eq:referential-solid-volume-fraction} to produce

\[ \varphi_{r}^{s}=\varphi_{0}^{s}+\left(J-\varphi_{r}^{s}\right)\sum\limits_{\sigma}\frac{M^{\sigma}c^{\sigma}}{\rho_{T}^{\sigma}} \]

which can be solved for \(\varphi_{r}^{s}\) using the concentrations \(c^{\sigma}\),

\[ \begin{equation} \varphi_{r}^{s}=\frac{\varphi_{0}^{s}+J\sum\limits_{\sigma}\frac{M^{\sigma}c^{\sigma}}{\rho_{T}^{\sigma}}}{1+\sum\limits_{\sigma}\frac{M^{\sigma}c^{\sigma}}{\rho_{T}^{\sigma}}}\,.\label{eq:phirs-solutes-as-SBMs} \end{equation} \]

Solutes

Solutes are denoted generically by \(\alpha=\iota\). In chemistry solute content is often represented in units of molar concentration (moles per fluid volume). It follows that solute molar concentration \(c^{\iota}\) and molar supply \(\hat{c}^{\iota}\) are related to \(\rho^{\iota}\) and \(\hat{\rho}^{\iota}\) via

\[ \begin{equation} c^{\iota}=\frac{\rho^{\iota}}{\left(1-\varphi^{s}\right)M^{\iota}},\quad\hat{c}^{\iota}=\frac{\hat{\rho}^{\iota}}{\left(1-\varphi^{s}\right)M^{\iota}}\,.\label{eq159} \end{equation} \]

The molar flux of constituent \(\iota\) relative to the solid is given by

\[ \begin{equation} \mathbf{j}^{\iota}=\left(1-\varphi^{s}\right)c^{\iota}\left(\mathbf{v}^{\iota}-\mathbf{v}^{s}\right)\,,\label{eq160} \end{equation} \]

where it may be noted that \(\mathbf{m}^{\iota}=M^{\iota}\mathbf{j}^{\iota}\). Combining these relations with Eqs.\eqref{eq153}-\eqref{eq154} produces the desired form of the mass balance for the solutes,

\[ \begin{equation} \frac{1}{J}\frac{D^{s}\left[J\left(1-\varphi^{s}\right)c^{\iota}\right]}{Dt}+\mbox{div}\mathbf{j}^{\iota}=\left(1-\varphi^{s}\right)\hat{c}^{\iota}\,.\label{eq161} \end{equation} \]

This form is suitable for implementation in a finite element analysis where the mesh is defined on the solid matrix.

Mixture with Negligible Solute Volume Fraction

The volume fraction of each constituent is given by \(\varphi^{\alpha}=\rho^{\alpha}/\rho_{T}^{\alpha}\). In a saturated mixture these volume fractions satisfy \(\sum\limits_{\alpha}\varphi^{\alpha}=1\). Substituting \(\rho^{\alpha}=\varphi^{\alpha}\rho_{T}^{\alpha}\) into eq.\eqref{eq150}, dividing across by \(\rho_{T}^{\alpha}\) (invariant for intrinsically incompressible constituents), and taking the sum of the resulting expression over all constituents produces

\[ \begin{equation} \divg\left(\sum\limits_{\alpha}\varphi^{\alpha}\mathbf{v}^{\alpha}\right)=\sum\limits_{\alpha}\frac{\hat{\rho}^{\alpha}}{\rho_{T}^{\alpha}}\,.\label{eq162} \end{equation} \]

This mass balance relation for the mixture expresses the fact that the mixture volume will change as a result of chemical reactions where the true density of products is different from that of reactants. Indeed, assuming that \(\rho_{T}^{\alpha}\) is the same for all \(\alpha\) would nullify the right-hand-side of eq.\eqref{eq162} based on eq.\eqref{eq151}. We now adopt the assumption that solutes occupy a negligible volume fraction of the mixture (\(\varphi^{\iota}\ll1)\), from which it follows that \(\varphi^{s}+\varphi^{w}\approx1\) and \(\sum\limits_{\alpha}\varphi^{\alpha}\mathbf{v}^{\alpha}\approx\mathbf{v}^{s}+\mathbf{w}\), where \(\mathbf{w}=\varphi^{w}\left(\mathbf{v}^{w}-\mathbf{v}^{s}\right)\) is the volumetrix flux of solvent relative to the solid. Thus, the mixture mass balance may be reduced to

\[ \begin{equation} \divg\left(\mathbf{v}^{s}+\mathbf{w}\right)=\sum\limits_{\alpha}\frac{\hat{\rho}^{\alpha}}{\rho_{T}^{\alpha}}\,.\label{eq163} \end{equation} \]

In the special case of the solvent (\(\alpha=w\)), FEBio uses a solvent supply, \(\hat{\varphi}^{w}=\hat{\rho}^{w}/\hat{\rho}_{T}^{w}\), which may be incorporated in eq.\eqref{eq163} as

\[ \begin{equation} \divg\left(\mathbf{v}^{s}+\mathbf{w}\right)=\hat{\varphi}^{w}+\sum\limits_{\alpha\ne w}\frac{\hat{\rho}^{\alpha}}{\rho_{T}^{\alpha}}\,.\label{eq163b} \end{equation} \]

From the mass balance equation for the solvent we have

\[ \begin{aligned}\hat{\varphi}^{w} & =\frac{D^{s}\varphi^{w}}{Dt}+\grad\varphi^{w}\cdot\left(\mathbf{v}^{w}-\mathbf{v}^{s}\right)+\varphi^{w}\divg\mathbf{v}^{w}\end{aligned} \]

so that the mixture mass balance takes the form

\[ \begin{aligned}\divg\left(\mathbf{v}^{s}+\mathbf{w}\right) & =\frac{D^{s}\varphi^{w}}{Dt}+\grad\varphi^{w}\cdot\left(\mathbf{v}^{w}-\mathbf{v}^{s}\right)+\varphi^{w}\divg\mathbf{v}^{w}+\sum\limits_{\alpha\ne w}\frac{\hat{\rho}^{\alpha}}{\rho_{T}^{\alpha}}\,\end{aligned} \]

Using \(\varphi^{w}+\varphi^{s}=1\) this relation further simplifies to

\[ \varphi^{s}\divg\mathbf{v}^{s}+\frac{D^{s}\varphi^{s}}{Dt}=\sum\limits_{\alpha\ne w}\frac{\hat{\rho}^{\alpha}}{\rho_{T}^{\alpha}}\,. \]

Now we use the kinematic relation \(\divg\mathbf{v}^{s}=\frac{1}{J}\frac{D^{s}J}{Dt}\) to reduce this equation further to

\[ \frac{D^{s}\left(J\varphi^{s}\right)}{Dt}=J\sum\limits_{\alpha\ne w}\frac{\hat{\rho}^{\alpha}}{\rho_{T}^{\alpha}}\,. \]

But recall that

\[ J\varphi^{s}=\varphi_{r}^{s} \]

thus the mixture mass balance now simplifies to

\[ \frac{1}{J}\frac{D^{s}\varphi_{r}^{s}}{Dt}=\sum\limits_{\alpha\ne w}\frac{\hat{\rho}^{\alpha}}{\rho_{T}^{\alpha}} \]

Since

\[ \hat{\rho}^{\alpha}=\left(1-\varphi^{s}\right)M^{\alpha}\hat{c}^{\alpha} \]

it follows that

\[ \frac{D^{s}\varphi_{r}^{s}}{Dt}=J\left(1-\varphi^{s}\right)\sum\limits_{\alpha\ne w}\frac{M^{\alpha}\hat{c}^{\alpha}}{\rho_{T}^{\alpha}} \]

Another way of presenting these results is to restate the mixture mass balance as

\[ \begin{equation} \divg\left(\mathbf{v}^{s}+\mathbf{w}\right)=\hat{\varphi}_{0}^{w}+\frac{1}{J}\frac{D^{s}\varphi_{r}^{s}}{Dt}\label{eq:mixture-mass-reactive} \end{equation} \]

where \(\hat{\varphi}_{0}^{w}\) is the solvent supply from user-specified sources (not from the chemical reactions).

Chemical Kinetics

Productions rates are described by constitutive relations which are functions of the state variables. In a biological mixture under isothermal conditions, the minimum set of state variables needed to describe reactive mixtures that include a solid matrix are: the (uniform) temperature \(\theta\), the solid matrix deformation gradient \(\mathbf{F}\) (or related strain measures), and the molar content \(c^{\alpha}\) of the various constituents. This set differs from the classical treatment of chemical kinetics in fluid mixtures by the inclusion of \(\mathbf{F}\) and the subset of constituents bound to the solid matrix. To maintain a consistent notation in this section, solid-bound molecular species are described by their molar concentrations and molar supplies which may be related to their referential mass density and referential mass supply according to

\[ \begin{equation} c^{\sigma}=\frac{\rho_{r}^{\sigma}}{\left(J-\varphi_{r}^{s}\right)M^{\sigma}},\quad\hat{c}^{\sigma}=\frac{\hat{\rho}_{r}^{\sigma}}{\left(J-\varphi_{r}^{s}\right)M^{\sigma}}\,.\label{eq164} \end{equation} \]

Consider a general chemical reaction,

\[ \begin{equation} \sum\limits_{\alpha}\nu_{R}^{\alpha}\mathcal{E}^{\alpha}\to\sum\limits_{\alpha}\nu_{P}^{\alpha}\mathcal{E}^{\alpha}\,,\label{eq165} \end{equation} \]

where \(\mathcal{E}^{\alpha}\) is the chemical species representing constituent \(\alpha\); \(\nu_{R}^{\alpha}\) and \(\nu_{P}^{\alpha}\) represent stoichiometric coefficients of the reactants and products, respectively. Since the molar supply of reactants and products is constrained by stoichiometry, it follows that all molar supplies \(\hat{c}^{\alpha}\) in a specific chemical reaction may be related to a production rate \(\hat{\zeta}\) according to

\[ \begin{equation} \hat{c}^{\alpha}=\nu^{\alpha}\hat{\zeta}\,,\label{eq166} \end{equation} \]

where \(\nu^{\alpha}\) represents the net stoichiometric coefficient for \(\mathcal{E}^{\alpha}\),

\[ \begin{equation} \nu^{\alpha}=\nu_{P}^{\alpha}-\nu_{R}^{\alpha}\,.\label{eq167} \end{equation} \]

Thus, formulating constitutive relations for \(\hat{c}^{\alpha}\) is equivalent to providing a single relation for \(\hat{\zeta}\left(\theta,\mathbf{F},c^{\alpha}\right)\). When the chemical reaction is reversible,

\[ \begin{equation} \sum\limits_{\alpha}\nu_{R}^{\alpha}\mathcal{E}^{\alpha}\rightleftharpoons\sum\limits_{\alpha}\nu_{P}^{\alpha}\mathcal{E}^{\alpha}\,,\label{eq168} \end{equation} \]

the relations of Eqs.\eqref{eq166}-\eqref{eq167} still apply but the form of \(\hat{\zeta}\) would be different.

Using the relations of Eqs.\eqref{eq159}, \eqref{eq164} and \eqref{eq166}, it follows in general that \(\hat{\rho}^{\alpha}=\left(1-\varphi^{s}\right)M^{\alpha}\nu^{\alpha}\hat{\zeta}\), so that the constraint of eq.\eqref{eq151} is equivalent to enforcing stoichiometry, namely,

\[ \begin{equation} \sum\limits_{\alpha}\nu^{\alpha}M^{\alpha}=0\,.\label{eq169} \end{equation} \]

Thus, properly balancing a chemical reaction satisfies this constraint.

The mixture mass balance in eq.\eqref{eq163b} may now be rewritten as

\[ \begin{equation} \divg\left(\mathbf{v}^{s}+\mathbf{w}^{w}\right)=\hat{\varphi}_{0}^{w}+\left(1-\varphi^{s}\right)\hat{\zeta}\overline{\mathcal{V}}\,,\label{eq170} \end{equation} \]

where \(\overline{\mathcal{V}}=\sum\limits_{\alpha\ne w}\nu^{\alpha}\mathcal{V}^{\alpha}\) and \(\mathcal{V}^{\alpha}=M^{\alpha}/\rho_{T}^{\alpha}\) is the molar volume of \(\alpha\). (Currently in FEBio, \(\hat{\varphi}_{0}^{w}\) is specified independently of \(\hat{\zeta}\), because users may choose to neglect the contribution from \(\overline{\mathcal{V}}\) in eq.\eqref{eq170}; therefore, if one desires to model chemical reactions, eq.\eqref{eq165} or eq.\eqref{eq168}, that involve the solvent, it is necessary to explicitly provide a solvent supply function compatible with the above relations, namely \(\hat{\varphi}_{0}^{w}=\left(1-\varphi^{s}\right)\hat{\zeta}\nu^{w}\mathcal{V}^{w}\).) Similarly, the solute mass balance in eq.\eqref{eq161} becomes

\[ \begin{equation} \frac{1}{J}\frac{D^{s}\left[J\left(1-\varphi^{s}\right)c^{\iota}\right]}{Dt}+\divg\mathbf{j}^{\iota}=\left(1-\varphi^{s}\right)\nu^{\iota}\hat{\zeta}\,.\label{eq:solute-mass-balance} \end{equation} \]

These mass balance equations reduce to those of non-reactive mixtures when \(\hat{\zeta}=0\).