2.16 Fluid-Solutes Analyses

Fluid-solutes analyses combine solute transport with fluid mechanics. The fluid-solutes domain, denoted by \(\Omega^{f}\), is a hybrid multiphasic domain without a deformable porous solid matrix ¹. Like the CFD domain, the CFD-solutes domain is fixed in space, and the fluid flows within it and possibly across its boundaries. The boundaries \(\partial\Omega^{f}\) of \(\Omega^{f}\), and interfaces within this domain (such as faces of adjoining finite elements), are denoted generically by \(\Gamma^{f}\). In the current implementation, finite element domains that meet at these interfaces share common nodes and element faces, though it is possible to consider future extensions where tied interfaces are used to connect dissimilar meshes. The formulation and finite element implementation of this type of domain are presented below, using only salient governing equations simplified from the more general hybrid multiphasic theory presented in our recent study ¹.

Governing Equations

We model the domain \(\Omega^{f}\) as a mixture of an isothermal, compressible, and viscous solvent which is electrically neutral, and isothermal, intrinsically incompressible, and dilute solute constituents which may be electrically charged. All solute accelerations are neglected, since they contribute insignificantly in comparison with diffusive forces, except in problems such as centrifugation, where they can be represented by an external inertial body force. Unlike the standard (Darcy flow through a porous-deformable) multiphasic domain ²³⁴, the fluid-solutes domain is fixed in space. We denote the mixture constituents using a superscripted \(\alpha\), where \(\alpha=f,\iota\) for solvent and solutes, respectively. The volume fraction \(\varphi^{\iota}\) of each solute species \(\iota\) is assumed to be negligible (\(\varphi^{\iota}\ll1\)), therefore the volume fraction of the solvent is \(\varphi^{f}=dV^{f}/dV\approx1\). The apparent density of each constituent \(\alpha\) is \(\rho^{\alpha}=dm^{\alpha}/dV\), which represents the ratio of the elemental mass \(dm^{\alpha}\) of constituent \(\alpha\) in the elemental fluid mixture volume. Under these assumptions, the true density of the solvent simplifies to \(\rho_{T}^{f}=\rho^{f}\) while the true density \(\rho_{T}^{\iota}\) of each solute \(\iota\) is constant due to its intrinsic incompressibility. Since the solvent is compressible, as was done in the fluid mechanics formulation ⁵ (Section Fluid Mechanics), the solvent volume ratio \(J^{f}\) was employed as a kinematic state variable to represent the fluid volumetric strain,

\[ \begin{equation} J^{f}=\frac{\rho_{r}^{f}}{\rho^{f}}\thinspace,\label{eq:Fluid-Dil-CFDSol} \end{equation} \]

where \(\rho_{r}^{f}\) is the referential solvent density (\(\rho^{f}\) in the reference configuration), which is a user-specified material constant. An alternative representation of fluid volumetric strain is the fluid dilatation, \(e^{f}=J^{f}-1\), which is used in the actual finite element code.

The axiom of mass balance for the mixture takes the form

\[ \begin{equation} \frac{D^{f}J^{f}}{Dt}=J^{f}\left(\divg\mathbf{v}^{f}-\sum_{\iota}\frac{\hat{\rho}^{\iota}}{\rho_{T}^{\iota}}\right)\thinspace,\label{eq:Mix-Mass-Bal-CFDSol} \end{equation} \]

where \(\mathbf{v}^{f}\) is the velocity of the fluid constituent and \(\hat{\rho}^{\iota}\) is the mass density supply to \(\iota\) due to reactions with all other solutes. It is assumed that reactions do not add or remove solvent mass (\(\hat{\rho}^{f}=0\)). Here, \(D^{f}\left(\cdot\right)/Dt=\partial\left(\cdot\right)/\partial t+\grad\left(\cdot\right)\cdot\mathbf{v}^{f}\) is the material time derivative operator in the spatial frame, following the fluid motion. The mixture mass balance reproduces the kinematic constraint between the fluid velocity and dilatation for the fluid mechanics formulation ⁵ (Section Fluid Mechanics), with the addition of reactive mass supply terms. The axiom of mass balance for each solute takes the form

\[ \begin{equation} \hat{c}^{\iota}=\frac{\partial c^{\iota}}{\partial t}+\divg\left(c^{\iota}\mathbf{v}^{\iota}\right)=\frac{1}{J^{f}}\frac{D^{f}\left(J^{f}c^{\iota}\right)}{Dt}+\divg\mathbf{j}^{\iota}\thinspace,\label{eq:Mass-Bal-Solute-CFDSol} \end{equation} \]

where \(\hat{c}^{\iota}=\hat{\rho}^{\iota}/M^{\iota}\) is the molar concentration supply and \(c^{\iota}=\rho^{\iota}/M^{\iota}\) is the molar concentration of solute \(\iota\), and \(\mathbf{j}^{\iota}=c^{\iota}\left(\mathbf{v}^{\iota}-\mathbf{v}^{f}\right)\) is the molar flux of solute \(\iota\) relative to the solvent. As shown previously ¹, the molar solute flux may be evaluated from the solute momentum balance when inertia terms have been neglected,

\[ \begin{equation} \mathbf{j}^{\iota}=\tilde{\kappa}^{\iota}d_{0}^{\iota}\left(-\grad\tilde{c}^{\iota}+\frac{M^{\iota}}{R\theta}\tilde{c}^{\iota}\mathbf{b}^{\iota}\right)\thinspace,\label{eq:Solute-Flux-CFDSol} \end{equation} \]

where \(\tilde{\kappa}^{\iota}\) is the partition coefficient of solute \(\iota\) relative to an ideal solution ⁶⁷. It is given by

\[ \begin{equation} \tilde{\kappa}^{\iota}=\frac{\kappa^{\iota}}{\gamma^{\iota}}\exp\left(-\frac{z^{\iota}F_{c}\psi}{R\theta}\right)\thinspace,\label{eq:Partition-Coeff-CFDSol} \end{equation} \]

where \(\kappa^{\iota}\) is the solubility of solute \(\iota\) in the mixture, or the fraction of fluid volume accessible to the solute ⁸ (\(\kappa^{\iota}=1\) in a fluid mixture), and \(\gamma^{\iota}\) is the activity coefficient of solute \(\iota\). The ratio \(\kappa^{\iota}/\gamma^{\iota}\) is a non-dimensional property that describes the deviation of the solute chemical potential from that of an ideal, dilute solution ⁹. We can represent this ratio as \(\hat{\kappa}^{\iota}\equiv\kappa^{\iota}/\gamma^{\iota}\), and call it the effective solubility of \(\iota\) in the solution ¹⁰. In \eqref{eq:Partition-Coeff-CFDSol}, \(z^{\iota}\) is the charge number of solute \(\iota\), \(F_{c}\) is Faraday's constant, \(\psi\) is the electrical potential, \(R\) is the universal gas constant, and \(\theta\) is the absolute temperature. In \eqref{eq:Solute-Flux-CFDSol}, \(d_{0}^{\iota}\) is the diffusivity of solute \(\iota\) in the fluid, \(\mathbf{b}^{\iota}\) is the body force per mass acting on solute \(\iota\), and \(\tilde{c}^{\iota}\) is the effective concentration of solute \(\iota\), such that

\[ \begin{equation} \tilde{c}^{\iota}=c^{\iota}/\tilde{\kappa}^{\iota}\,.\label{eq:effective-actual-concentration} \end{equation} \]

The effective solute concentration represents an alternative form of the solute electrochemical potential, expressed in units of molar concentration ²³. The solute flux in \eqref{eq:Solute-Flux-CFDSol} has contributions from diffusion (as represented by the term containing \(\grad\tilde{c}^{\iota}\)), sedimentation (\(\tilde{c}^{\iota}\mathbf{b}^{\iota}\) when \(\mathbf{b}^{\iota}\) is proportional to an inertial force), electrophoresis (\(\tilde{c}^{\iota}\mathbf{b}^{\iota}\) when \(\mathbf{b}^{\iota}\) is proportional to the Lorentz force acting on charged solute \(\iota\)), and convection (\(\tilde{c}^{\iota}\mathbf{v}^{f}\)).

The solute flux in eq.\eqref{eq:Solute-Flux-CFDSol} may also be expressed in terms of its components,

\[ \begin{equation} \mathbf{j}^{\iota}=\underbrace{-\tilde{\kappa}^{\iota}d_{0}^{\iota}\grad\tilde{c}^{\iota}}_{\mathbf{j}_{d}^{\iota}}+\underbrace{s^{\iota}c^{\iota}\mathbf{b}^{\iota}}_{\mathbf{j}_{b}^{\iota}}\label{eq:solute-flux-components} \end{equation} \]

where

\[ \begin{equation} s^{\iota}\equiv d_{0}^{\iota}\frac{M^{\iota}}{R\theta}\label{eq:sedimentation-coeff} \end{equation} \]

is the sedimentation coefficient. In the expression of eq.\eqref{eq:solute-flux-components} \(\mathbf{j}_{d}^{\iota}\) represents the diffusive flux, and \(\mathbf{j}_{b}^{\iota}\) is the sedimentation flux. Now the mass balance expression in eq.\eqref{eq:Mass-Bal-Solute-CFDSol} may be combined with eq.\eqref{eq:Mix-Mass-Bal-CFDSol} to produce

\[ \begin{equation} \hat{c}^{\iota}-c^{\iota}\sum_{\gamma}\mathcal{V}^{\gamma}\hat{c}^{\gamma}=\frac{1}{J^{f}}\frac{D^{f}\left(J^{f}c^{\iota}\right)}{Dt}+\divg\left(\mathbf{j}_{d}^{\iota}+\mathbf{j}_{b}^{\iota}\right)\label{eq:solute-mass-alt} \end{equation} \]

where \(\mathcal{V}^{\iota}\equiv\nicefrac{M^{\iota}}{\rho_{T}^{\iota}}\) is the molar volume of solute \(\iota\). Since we have neglected the volume fraction of solutes in the current mixture formulation, it follows that \(c^{\iota}\mathcal{V}^{\iota}\ll1\). Thus, we may neglect the second term on the left-hand side of the above equation, as well as the related term on the right-hand side of eq.\eqref{eq:Mix-Mass-Bal-CFDSol}.

The fluid momentum balance is given by

\[ \begin{equation} \rho^{f}\mathbf{a}^{f}=-\grad\tilde{p}+\divg\boldsymbol{\tau}+\rho^{f}\mathbf{b}^{f}+\sum_{\iota}\left(-R\theta\tilde{\kappa}^{\iota}\grad\tilde{c}^{\iota}+M^{\iota}c^{\iota}\mathbf{b}^{\iota}\right)\thinspace,\label{eq:Fluid-Mtm-Bal-CFDSol} \end{equation} \]

where \(\mathbf{a}^{f}=D^{f}\mathbf{v}^{f}/Dt\) is the fluid acceleration, \(\boldsymbol{\tau}\) is the fluid viscous stress, \(\mathbf{b}^{f}\) is the body force per mass acting on the fluid \(f\), and \(\tilde{p}\) is the effective fluid pressure, related to the fluid pressure \(p\) via

\[ \begin{equation} p=\tilde{p}+R\theta\Phi\sum_{\iota}c^{\iota}\thinspace.\label{eq:Total-fluid-p-CFDSol} \end{equation} \]

Here, \(\Phi\) is the osmotic coefficient, a non-dimensional property that describes the deviation of the osmotic pressure from the ideal behavior known as van't Hoff's law ¹¹, and \(R\theta\Phi\sum_{\iota}c^{\iota}\) is the osmotic or chemical contribution to the fluid pressure. Therefore, \(\tilde{p}\), which is an alternative form of the solvent mechano-chemical potential, may be interpreted as the mechanical or hydraulic contribution to \(p\).

Previously, we showed that the effective pressure must only depend on the fluid volume ratio, such that \(\tilde{p}=\tilde{p}\left(J^{f}\right)\) ¹. Based on mass, momentum and energy jump conditions for this type of mixture, it can also be shown that \(J^{f}\) (or alternatively \(e^{f}\)) is continuous across interfaces \(\Gamma^{f}\), and thus may be used as a nodal degree of freedom (DOF). We may choose a constitutive relation valid for arbitrary (but strictly positive) values of \(J^{f}\)

\[ \begin{equation} \tilde{p}=K\left(\frac{1}{J^{f}}-1\right)\,,\label{eq:Eff-Pressure-CFDSol} \end{equation} \]

where \(K\) is the solvent's bulk modulus (a user-specified material property).

The solute-dependent terms appearing on the right side of the solvent momentum balance \eqref{eq:Fluid-Mtm-Bal-CFDSol} are typically neglected in classical CFD-solutes formulations, but we included them here since they emerged naturally to describe the phenomena of osmosis (solvent flow due to solute diffusion or sedimentation), and electro-osmosis (solvent flow due to solute electrophoresis). To reproduce classical formulations, by default we turned off these terms (within the summation on the right-hand side of eq.\eqref{eq:Fluid-Mtm-Bal-CFDSol}), and we added a user-defined flag (“include osmosis”) to optionally turn them back on. In addition, users can turn off osmotic pressure by setting the osmotic coefficient \(\Phi\) to zero, in which case \(p=\tilde{p}\) according to eq.\eqref{eq:Total-fluid-p-CFDSol}. Both of these approaches were investigated in the test problems described below.

As done in the standard multiphasic formulation (Section Triphasic and Multiphasic Materials), we assumed there can be no electric charge accumulation in the mixture, so we enforce the electroneutrality condition

\[ \begin{equation} \sum_{\iota}z^{\iota}c^{\iota}=0\,.\label{eq:Electroneutrality-CFDSol} \end{equation} \]

Multiplying the mass balances from \eqref{eq:Mass-Bal-Solute-CFDSol} by \(z^{\iota}\) and taking the sum over all the constituents, making use of \eqref{eq:Electroneutrality-CFDSol} and adopting the assumption that reactive processes maintain electroneutrality (\(\sum_{\iota}z^{\iota}\nu^{\iota}=0\) where \(\nu^{\alpha}\) is the net stoichiometric coefficient of \(\alpha\) ⁴) produces the constraint

\[ \begin{equation} 0=\divg\sum_{\iota}z^{\iota}\mathbf{j}^{\iota}\,,\label{eq:Current-Density-Constraint} \end{equation} \]

or equivalently \(\divg\mathbf{I}_{e}=0\), where \(\mathbf{I}_{e}=F_{c}\sum_{\iota}z^{\iota}\mathbf{j}^{\iota}\) is the electric current density.

In the CFD-solutes implementation the fluid velocity \(\mathbf{v}^{f}\), fluid dilatation \(e^{f}\), and effective solute concentrations \(\tilde{c}^{\iota}\) are the nodal degrees of freedom \(\left(\mathbf{v}^{f},e^{f},\tilde{c}^{\iota}\right)\), as further addressed in Section Jump Conditions. These choices imply that \(\mathbf{\mathbf{v}}^{f}\), \(e^{f}\), and \(\tilde{c}^{\iota}\) are continuous across finite element boundaries, as justified below when we review jump conditions across interfaces. Note that in the absence of solutes, the governing equations \eqref{eq:Mix-Mass-Bal-CFDSol} and \eqref{eq:Fluid-Mtm-Bal-CFDSol} simplify to the original CFD equations, with \(p=\tilde{p}\). More information on calculating the supply terms \(\hat{c}^{\iota}\) and \(\hat{\rho}^{\iota}\) can be found in Section Chemical Reactions. The electric potential \(\psi\) can be determined as described in the previous standard multiphasic solver ³⁴.

Jump Conditions

Jump conditions on the axioms of mass, momentum and energy balance are needed to determine which variables may be selected as nodal degrees of freedom in the finite element implementation, and which tractions are naturally continuous across an interface. The jump conditions for the CFD-solutes material were determined from the jump conditions of the general hybrid multiphasic material ¹. Here, we summarize the relevant results, which apply to any interface \(\Gamma^{f}\). Here, since the domain is fixed in space, the interface \(\Gamma^{f}\) is stationary. We employ the notation \(\left[\left[f\right]\right]=f_{+}-f_{-}\) to denote the jump in an arbitrary function \(f\) across the interface \(\Gamma^{f}\), with \(f_{+}\) and \(f_{-}\) denoting the values of \(f\) on either side of \(\Gamma^{f}\). The unit normal on \(\Gamma^{f}\) is \(\mathbf{n}\), which points away from the \(+\) side. A variable \(f\) which is continuous across \(\Gamma^{f}\) satisfies \(\left[\left[f\right]\right]=0\).

Based on the jump condition on the axiom of mass balance for the fluid, the normal component of the mass flux of the fluid is continuous across \(\Gamma^{f}\), \(\left[\left[\rho^{f}\mathbf{v}^{f}\right]\right]\cdot\mathbf{n}=0\). Furthermore, a sufficient condition to satisfy the jump on the axiom of energy balance is to enforce continuity of the solvent specific free enthalpy (or Gibbs function), such that \(\left[\left[\psi^{f}+p/\rho_{T}^{f}\right]\right]=\left[\left[\mu^{f}\right]\right]=0\), where \(\psi^{f}\) is the fluid specific free energy. In the multiphasic mixture literature, the solvent specific free enthalpy \(\mu^{f}\) is also called the mechano-chemical potential ¹². This jump condition applies only when there is solvent on both sides of the interface, such as across \(\Gamma^{f}\). In an isothermal framework, the solvent mechano-chemical potential is a function of state that only depends on \(J^{f}\) such that \(\mu^{f}=\mu^{f}\left(J^{f}\right)\) ¹. Therefore this energy jump condition implies that \(J^{f}\) must be continuous across \(\Gamma^{f}\), thus

\[ \begin{equation} \left[\left[J^{f}\right]\right]=0\,,\label{eq:J-Jump-CFDSol} \end{equation} \]

also implying that \(\left[\left[\tilde{p}\right]\right]=0\). Given \eqref{eq:Fluid-Dil-CFDSol}, it follows that \(\left[\left[\rho^{f}\right]\right]=0\) and the fluid mass balance jump condition reduces to \(\left[\left[\mathbf{v}^{f}\right]\right]\cdot\mathbf{n}=0\), implying that the fluid velocity component normal to \(\Gamma^{f}\) must be continuous. For the tangential component of \(\mathbf{v}^{f}\) on \(\Gamma^{f}\) we use the no-slip condition for viscous fluids, which states that the velocity component tangential to \(\Gamma^{f}\) must be continuous across that interface, as is standard in CFD solvers. Combining these two jump conditions produces

\[ \begin{equation} \left[\left[\mathbf{v}^{f}\right]\right]=0\,.\label{eq:Mass-Jump-Final-CFDSol} \end{equation} \]

From the solute mass jump condition, we can show that the normal component of the molar solute flux is continuous across \(\Gamma^{f}\),

\[ \begin{equation} \left[\left[\mathbf{j}^{\iota}\right]\right]\cdot\mathbf{n}=0\,.\label{eq:Mass-Jump-Sol-CFDSol} \end{equation} \]

The momentum jump condition requires that the mixture traction be continuous across \(\Gamma^{f}\), thus \(\left[\left[-p\mathbf{I}+\boldsymbol{\tau}\right]\right]\cdot\mathbf{n}=\mathbf{0}\). Another relation which is thermodynamically sufficient to satisfy the jump condition on the energy balance is the continuity of the viscous fluid traction,

\[ \begin{equation} \left[\left[\boldsymbol{\tau}\right]\right]\cdot\mathbf{n}=0\,.\label{eq:Fluid-Mtm-Jump-CFDSol} \end{equation} \]

This jump condition \eqref{eq:Fluid-Mtm-Jump-CFDSol}, which also applies only if fluid is present on both sides of \(\Gamma\), implies that \(\left[\left[p\right]\right]=0\).

Finally, the last relation that is thermodynamically sufficient to satisfy the energy balance jump condition applies to the electrochemical potential for the solutes, where \(\left[\left[\tilde{\mu}^{\iota}\right]\right]=0\). For a solute, the general constitutive relation for the electrochemical potential is \(\tilde{\mu}^{\iota}=\frac{R\theta}{M^{\iota}}\ln\frac{c^{\iota}}{\tilde{\kappa}^{\alpha}c_{0}^{\iota}}\), and the jump can be simplified to

\[ \begin{equation} \left[\left[\tilde{c}^{\iota}\right]\right]=0\,.\label{eq:Conc-Jump-CFDSol} \end{equation} \]

Thus, letting \(\mathbf{v}^{f}\), \(e^{f}\), and \(\tilde{c}^{\iota}\) be nodal degrees of freedom automatically enforces the jump conditions \eqref{eq:J-Jump-CFDSol}, \eqref{eq:Mass-Jump-Final-CFDSol}, and \eqref{eq:Conc-Jump-CFDSol}, acting as essential boundary conditions. As stated earlier, the jump conditions are the same as the previous CFD formulation in the case of zero solute concentration (all \(\tilde{c}^{\iota}=0\)) ⁵.

---

1. Shim, Jay J; Ateshian, Gerard A. "A Hybrid Reactive Multiphasic Mixture With a Compressible Fluid Solvent." J Biomech Eng, vol. 144 (2022). ↩↩↩↩↩↩

2. Ateshian, G. A.. "On the theory of reactive mixtures for modeling biological growth." Biomech Model Mechanobiol, vol. 6, pp. 423-45 (2007). ↩↩

3. Gerard A. Ateshian; Steve Maas; Jeffrey A. Weiss. "Multiphasic Finite Element Framework for Modeling Hydrated Mixtures With Multiple Neutral and Charged Solutes." J. Biomech. Eng., vol. 135 (2013). ↩↩↩

4. Gerard A. Ateshian; Robert J. Nims; Steve Maas; Jeffrey A. Weiss. "Computational modeling of chemical reactions and interstitial growth and remodeling involving charged solutes and solid-bound molecules." Biomech. Model. Mechanobiol., vol. 13, pp. 1105--1120 (2014). ↩↩↩

5. Ateshian, Gerard A; Shim, Jay J; Maas, Steve A; Weiss, Jeffrey A. "Finite Element Framework for Computational Fluid Dynamics in FEBio." J Biomech Eng, vol. 140 (2018). ↩↩↩

6. Ogston, A. G.; Phelps, C. F.. "The partition of solutes between buffer solutions and solutions containing hyaluronic acid." Biochem J, vol. 78, pp. 827-33 (1961). ↩

7. Laurent, Torvard C.; Killander, Johan. "A Theory of Gel Filtration and its Experimental Verification." J Chromatogr, vol. 14, pp. 317-330 (1963). ↩

8. Mauck, R. L.; Hung, C. T.; Ateshian, G. A.. "Modeling of neutral solute transport in a dynamically loaded porous permeable gel: implications for articular cartilage biosynthesis and tissue engineering." J Biomech Eng, vol. 125, pp. 602-14 (2003). ↩

9. Tinoco Jr., I.; Sauer, K.; Wang, J. C.. "Physical chemistry : principles and applications in biological sciences." Prentice Hall (1995). ↩

10. Ateshian, G.A.; Albro, M. B.; Maas, S.A.; Weiss, J.A.. "Finite element implementation of mechanochemical phenomena in neutral deformable porous media under finite deformation.." Journal of Biomechanical Engineering, vol. 133, pp. 1005-1017 (2011). ↩

11. McNaught, Alan. "Compendium of chemical terminology : IUPAC recommendations." Blackwell Science (1997). ↩

12. Gerard A. Ateshian. "Mixture Theory for Modeling Biological Tissues: Illustrations from Articular Cartilage." Studies in Mechanobiology, Tissue Engineering and Biomaterials, pp. 1--51 (2016). ↩