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 molar production rate $\hat{\zeta}$ according to $$\begin{equation} \hat{c}^{\alpha}=\nu^{\alpha}\hat{\zeta}\,,\label{eq:cr-concentration-supply} \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{eq:cr-net-stoich-coeff} \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{eq:cr-reversible-reaction} \end{equation}$$ the relations of (4.11.1-2)-(4.11.1-3) still apply but the constitutive relation for $\hat{\zeta}$ would be different.

Example:

Consider the dissociation of CaCl$_{\mathrm{2}}$ into ions Ca$^{\mathrm{2+}}$ and Cl$^{\mathrm{-}}$ , $$\mbox{CaCl}_{2}\rightleftharpoons\mbox{Ca}^{2+}+2\mbox{Cl}^{-}.$$ The mixture contains three constituents. The stoichiometric coefficients of the reactants are $\nu_{R}^{\mbox{CaCl}_{2}}=1$ , $\nu_{R}^{\mbox{Ca}^{2+}}=0$ , $\nu_{R}^{\mbox{Cl}^{-}}=0$ , and those of the products are $\nu_{P}^{\mbox{CaCl}_{2}}=0$ , $\nu_{P}^{\mbox{Ca}^{2+}}=1$ , $\nu_{P}^{\mbox{Cl}^{-}}=2$ .

The reaction production rate $\hat{\zeta}$ enters into the governing equations of multiphasic mixtures via the mass balance relation for each solute, $$\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:cr-solute-mass-balance} \end{equation}$$ the mass balance for the mixture, $$\begin{equation} \divg\left(\mathbf{v}^{s}+\mathbf{w}\right)=\left(1-\varphi^{s}\right)\hat{\zeta}\overline{\mathcal{V}}\,,\label{eq:cr-mixture-mass-balance} \end{equation}$$ where $\overline{\mathcal{V}}=\sum\limits _{\alpha}\nu^{\alpha}\mathcal{V}^{\alpha}$ and $\mathcal{V}^{\alpha}=M^{\alpha}/\rho_{T}^{\alpha}$ is the molar volume of $\alpha$ , and the mass balance for solid-bound constituents, $$\begin{equation} D^{s}\rho_{r}^{\sigma}/Dt=\hat{\rho}_{r}^{\sigma}\,,\label{eq:cr-sbm-mass-balance} \end{equation}$$ where $\rho_{r}^{\sigma}$ is the referential apparent mass density (mass of $\sigma$ per mixture volume in the reference configuration), and $\hat{\rho}_{r}^{\sigma}$ is the referential apparent mass supply of solid constituent $\sigma$ , related to molar concentrations and supplies via $$\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{eq:concentration-density-relations} \end{equation}$$ Internally, the content of solid-bound species is stored in $\rho_{r}^{\sigma}$ and (4.11.1-8) is used to evaluate $c^{\sigma}$ when needed for the calculation of $\hat{\zeta}$ . If a solid-bound molecule is involved in a chemical reaction, equation (4.11.1-7) is integrated to produce an updated value of $\rho_{r}^{\sigma}$ , using $\hat{\rho}_{r}^{\sigma}=\left({J-\varphi_{r}^{s}}\right)M^{\sigma}\nu^{\sigma}\hat{\zeta}$ based on (4.11.1-2) and (4.11.1-8).

Evolving solid content due to chemical reactions implies that the referential solid volume fraction $\varphi_{r}^{s}$ may not remain constant. This value is updated at every time point using $$\begin{equation} \varphi_{r}^{s}=\varphi_{0}^{s}+\sum\limits _{\sigma}\frac{\rho_{r}^{\sigma}}{\rho_{T}^{\sigma}}\,,\label{eq:cr-referential-solid-vol-frac} \end{equation}$$ where $\varphi_{0}^{s}$ is the solid volume fraction specified by the multiphasic material parameter phi0 (Section 4.10.2↑). Thus, $\varphi_{0}^{s}$ may be used to account for the solid volume fraction not contributed explicitly by solid-bound molecules. Based on kinematics, the solid volume fraction in the current configuration is given by $\varphi^{s}=\varphi_{r}^{s}/J$ . Therefore, since $0\leqslant\varphi^{s}\leqslant1$ by definition, it follows that $0\leqslant\varphi_{r}^{s}\leqslant J$ , implying that the referential solid volume fraction may evolve to values greater than unity when growth leads to swelling of the multiphasic mixture.

Similarly, if solid-bound molecules are charged and their content evolves over time, the referential fixed charge density may also evolve with chemical reactions according to $$\begin{equation} c_{r}^{F}=c_{0}^{F}+\frac{1}{1-\varphi_{r}^{s}}\sum\limits _{\sigma}\frac{z^{\sigma}\rho_{r}^{\sigma}}{M^{\sigma}}\,,\label{eq:cr-referential-FCD} \end{equation}$$ where $c_{0}^{F}$ is the referential fixed charge density specified by the multiphasic material parameter fixed_charge_density (Section 4.10.2↑). Thus, $c_{0}^{F}$ may be used to account for the fixed charge density not contributed explicitly by solid-bound molecules.

A chemical reaction is properly balanced when $$\begin{equation} \sum\limits _{\alpha}\nu^{\alpha}M^{\alpha}=0\,,\label{eq:cr-balanced-reaction} \end{equation}$$ where $M^{\alpha}$ is the molar mass of $\alpha$ . This constraint implies that the net gain in mass of products must be the same as the net loss in mass of reactants. However, this constraint is not verified in the code, allowing users to model chemical reactions with implicit constituents (constituents that are neither explicitly modeled as solutes nor as solid-bound molecules, for which $\nu^{\alpha}$ and $M^{\alpha}$ are not given). For example, a chemical reaction where cells consume glucose to form a protein from amino-acids building blocks may have the form $$\mbox{cells}+\mbox{glucose}+\mbox{amino acids}\to\mbox{cells}+\mbox{protein}+\mbox{waste products}\,.$$ The user may opt to model only the glucose reactant and the protein product explicitly, while the presence of all other species in this reaction is implicit. In these types of analyses the user must beware of potential inconsistencies in the evolving mass of reactants and products since only some of those constituents are modeled explicitly. In particular, the evolution of $\varphi_{r}^{s}$ as given in (4.11.1-9) can only account for the explicitly modeled solid-bound molecules. Furthermore, when some reactants and products are implicit, the value of the reaction molar volume $\overline{\mathcal{V}}$ calculated in the code becomes inaccurate and may produce unexpected results in the evaluation of the mixture mass balance relation in (4.11.1-6). Therefore, the user is given the option to override the value of $\overline{\mathcal{V}}$ calculated in the code. In particular, if the precise molar volumes of all the species in a reaction are not known, assuming that $\overline{\mathcal{V}}\approx0$ is a reasonable choice equivalent to assuming that all the constituents have approximately the same density $\rho_{T}^{\alpha}$ , as may be deduced from (4.11.1-11).

Since the electroneutrality condition is enforced in multiphasic mixtures in FEBio, it follows that chemical reactions must not violate this condition. Enforcing electroneutrality in a chemical reaction is equivalent to satisfying $$\begin{equation} \sum\limits _{\alpha}z^{\alpha}\nu^{\alpha}=0\,.\label{eq:cr-electroneutrality} \end{equation}$$ This constraint is checked within the code and an error is generated when it is violated.

A constitutive relation must be provided for the molar production rate $\hat{\zeta}\left(\theta,\mathbf{F},c^{\alpha}\right)$ of each chemical reaction.