2.8 Biphasic-Solute Material

A biphasic-solute material is an extension of the biphasic material model that also includes transport and mechano-chemical effects of a neutral solute. Transport of a solute in a porous medium includes diffusion, resulting from gradients in the solute concentration, and convection of the solute by the solvent, as a result of fluid pressure gradients. Mechano-chemical effects describe phenomena such as osmotic pressurization and swelling.

Governing Equations

The governing equations adopted in this finite element implementation of neutral solute transport in deformable porous media are based on the framework of mixture theory ¹². A single solute is considered in this presentation for notational simplicity, though the extension of equations to multiple solutes is straightforward. Various forms of the governing equations have been presented in the prior literature ³⁴, though a presentation that incorporates all the desired features of this implementation has not been reported previously and is thus detailed here.

The fundamental modeling assumptions adopted in this treatment are quasi-static conditions for momentum balance (negligible effects of inertia), intrinsic incompressibility of all constituents (invariant true densities), isothermal conditions, negligible volume fraction of solute relative to the solid and solvent, and negligible effects of solute and solvent viscosities (friction within constituents) relative to frictional interactions between constituents. These assumptions are often made in studies of biological tissues and cells. External body forces and chemical reactions are not considered.

The three constituents of the mixture are the porous-permeable solid matrix (\(\alpha=s)\), the solvent (\(\alpha=w)\), and the solute (\(\alpha=u)\). The motion of the solid matrix is described by the displacement vector \(\mathbf{u}\), the pressure of the interstitial fluid (solvent \(+\) solute) is \(p\), and the solute concentration (on a solution-volume basis) is \(c\). The total (or mixture) stress may be described by the Cauchy stress tensor \(\boldsymbol{\sigma}=-p\mathbf{I}+\boldsymbol{\sigma}^{e}\), where \(\mathbf{I}\) is the identity tensor and \(\boldsymbol{\sigma}^{e}\) is the stress arising from the strain in the porous solid matrix. Because it is porous, the solid matrix is compressible since the volume of pores changes as interstitial fluid enters or leaves the matrix. Under the conditions outlined above, the balance of linear momentum for the mixture reduces to

\[ \begin{equation} \divg\boldsymbol{\sigma}=-\grad p+\divg\boldsymbol{\sigma}^{e}=\mathbf{0}\,.\label{eq108} \end{equation} \]

Similarly, the equations of balance of linear momentum for the solvent and solute are given by

\[ \begin{equation} \begin{aligned}\rho^{w}\grad\tilde{\mu}^{w}+\mathbf{f}^{ws}\cdot\left(\mathbf{v}^{s}-\mathbf{v}^{w}\right)+\mathbf{f}^{wu}\cdot\left(\mathbf{v}^{u}-\mathbf{v}^{w}\right) & =\mathbf{0}\,,\\ -\rho^{u}\grad\tilde{\mu}^{u}+\mathbf{f}^{us}\cdot\left(\mathbf{v}^{s}-\mathbf{v}^{u}\right)+\mathbf{f}^{uw}\cdot\left(\mathbf{v}^{w}-\mathbf{v}^{u}\right) & =\mathbf{0}\,, \end{aligned} \label{eq109} \end{equation} \]

where \(\rho^{\alpha}\) is the apparent density (mass of \(\alpha\) per volume of the mixture), \(\tilde{\mu}^{\alpha}\) is the mechano-chemical potential and \(\mathbf{v}^{\alpha}\) is the velocity of constituent \(\alpha\). \(\mathbf{f}^{\alpha\beta}\) is the diffusive drag tensor between constituents \(\alpha\) and \(\beta\) representing momentum exchange via frictional interactions, which satisfies \(\mathbf{f}^{\beta\alpha}=\mathbf{f}^{\alpha\beta}\). An important feature of these relations is the incorporation of momentum exchange term between the solute and solid matrix, \(\mathbf{f}^{us}\cdot\left(\mathbf{v}^{s}-\mathbf{v}^{u}\right)\), which is often neglected in other treatments but plays an important role for describing solid-solute interactions ³⁵⁶. These momentum equations show that the driving force for the transport of solvent or solute is the gradient in its mechano-chemical potential, which is resisted by frictional interactions with other constituents.

The mechano-chemical potential is the sum of the mechanical and chemical potentials. The chemical potential \(\mu^{\alpha}\) of \(\alpha\) represents the rate at which the mixture free energy changes with increasing mass of \(\alpha\). The mechanical potential represents the rate at which the mixture free energy density changes with increasing volumetric strain of \(\alpha\). In a mixture of intrinsically incompressible constituents, where the volumetric strain is idealized to be zero, this potential is given by \(\left(p-p_{0}\right)/\rho_{T}^{\alpha}\), where \(\rho_{T}^{\alpha}\) is the true density of \(\alpha\) (mass of \(\alpha\) per volume of \(\alpha)\), which is invariant for incompressible constituents, and \(p_{0}\) is some arbitrarily set reference pressure (e.g., ambient pressure).

From classical physical chemistry, the general form of a constitutive relation for the chemical potential is \(\mu^{\alpha}=\mu_{0}^{\alpha}\left(\theta\right)+\left(R\theta/M^{\alpha}\right)\ln a^{\alpha}\) ⁷, where \(R\) is the universal gas constant, \(\theta\) is the absolute temperature, \(M^{\alpha}\) is the molecular weight (invariant) and \(a^{\alpha}\) is the activity of constituent \(\alpha\) (a non-dimensional quantity); \(\mu_{0}^{\alpha}\left(\theta\right)\) is the chemical potential at some arbitrary reference state, at a given temperature. For solutes, physical chemistry treatments let \(a^{u}=\gamma c/c_{0}\), where \(c_{0}\) is the solute concentration in some standard reference state (an invariant, typically \(c_{0}=1\,\mbox{M})\), and \(\gamma\) is the non-dimensional activity coefficient, which generally depends on the current state (e.g., concentration) but reduces to unity under the assumption of ideal physico-chemical behavior ⁷. Since this representation is strictly valid for free solutions only, whereas solutes may be partially excluded from some of the interstitial space of a porous solid matrix, Mauck et al. ³ extended this representation of the solute activity to let \(a^{u}=\gamma c/\kappa c_{0}\), where the solubility \(\kappa\) represents the fraction of the pore space which is accessible to the solute (\(0<\kappa\leqslant1)\). In this extended form, it becomes clear that even under ideal behavior (\(\gamma=1)\), the solute activity may be affected by the solubility. Indeed, for neutral solutes, the solubility also represents the partition coefficient of the solute between the tissue and external bath ⁸⁹.

When accounting for the fact that the solute volume fraction is negligible compared to the solvent volume fraction ⁷¹⁰, the general expressions for \(\tilde{\mu}^{w}\) and \(\tilde{\mu}^{u}\) take the form

\[ \begin{equation} \begin{aligned}\tilde{\mu}^{w} & =\mu_{0}^{w}\left(\theta\right)+\frac{1}{\rho_{T}^{w}}\left(p-p_{0}-R\theta\,\Phi\,c\right)\,,\\ \tilde{\mu}^{u} & =\mu_{0}^{u}\left(\theta\right)+\frac{R\theta}{M}\ln\frac{\gamma c}{\kappa c_{0}}\,, \end{aligned} \label{eq110} \end{equation} \]

where \(\Phi\) is the osmotic coefficient (a non-dimensional function of the state), which deviates from unity under non-ideal physico-chemical behavior. Therefore, a complete description of the physico-chemical state of solvent and solute requires constitutive relations for \(\Phi\) and the effective solubility \(\tilde{\kappa}=\kappa/\gamma\), which should generally depend on the solid matrix strain and the solute concentration.

It is also necessary to satisfy the balance of mass for each of the constituents. In the absence of chemical reactions, the statement of balance of mass for constituent \(\alpha\) reduces to

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

The apparent density may be related to the true density via \(\rho^{\alpha}=\varphi^{\alpha}\rho_{T}^{\alpha}\), where \(\varphi^{\alpha}\) is the volume fraction of \(\alpha\) in the mixture. Due to mixture saturation (no voids), the volume fractions add up to unity. Since the volume fraction of solute is considered negligible (\(\varphi^{u}\ll\varphi^{s},\varphi^{w})\), it follows that \(\sum\nolimits_{\alpha}\varphi^{\alpha}\approx\varphi^{s}+\varphi^{w}=1\). Since \(\rho_{T}^{\alpha}\) of an incompressible constituent is invariant in space and time, these relations may be combined to produce the mixture balance of mass relation,

\[ \begin{equation} \divg\left(\mathbf{v}^{s}+\mathbf{w}\right)=0,\label{eq112} \end{equation} \]

where \(\mathbf{w}=\varphi^{w}\left(\mathbf{v}^{w}-\mathbf{v}^{s}\right)\) is the volumetric flux of solvent relative to the solid. The balance of mass for the solute may also be written as

\[ \begin{equation} \frac{\partial\left(\varphi^{w}c\right)}{\partial t}+\divg\left(\mathbf{j}+\varphi^{w}c\,\mathbf{v}^{s}\right)=0\,,\label{eq113} \end{equation} \]

where \(\mathbf{j}=\varphi^{w}c\left(\mathbf{v}^{u}-\mathbf{v}^{s}\right)\) is the molar flux of solute relative to the solid. This mass balance relation is obtained by recognizing that the solute apparent density (mass per mixture volume) is related to its concentration (moles per solution volume) via \(\rho^{u}=\left(1-\varphi^{s}\right)Mc\approx\varphi^{w}Mc\). Finally, it can be shown via standard arguments that the mass balance for the solid matrix reduces to

\[ \begin{equation} \varphi^{s}=\frac{\varphi_{r}^{s}}{J}\,,\label{eq114} \end{equation} \]

where \(\varphi_{r}^{s}\) is the solid volume fraction in the reference state, \(J=\det\mathbf{F}\) and \(\mathbf{F}=\mathbf{I}+\grad\mathbf{u}\) is the deformation gradient of the solid matrix.

Inverting the momentum balance equations in \eqref{eq109}, it is now possible to relate the solvent and solute fluxes to the driving forces according to

\[ \begin{equation} \begin{aligned}\mathbf{w} & =-\tilde{\mathbf{k}}\cdot\left(\rho_{T}^{w}\grad\tilde{\mu}^{w}+Mc\frac{\mathbf{d}}{d_{0}}\grad\tilde{\mu}^{u}\right)\,,\\ \mathbf{j} & =\mathbf{d}\cdot\left(-\frac{M}{R\theta}\varphi^{w}c\grad\tilde{\mu}^{u}+\frac{c}{d_{0}}\mathbf{w}\right)\,, \end{aligned} \label{eq115} \end{equation} \]

where \(\mathbf{d}\) is the solute diffusivity tensor in the mixture (solid \(+\) solution), \(d_{0}\) is its (isotropic) diffusivity in free solution; \(\mathbf{\tilde{k}}\) is the hydraulic permeability tensor of the solution (solvent \(+\) solute) through the porous solid matrix, which depends explicitly on concentration according to

\[ \begin{equation} \tilde{\mathbf{k}}=\left[\mathbf{k}^{-1}+\frac{R\theta c}{\varphi^{w}d_{0}}\left(\mathbf{I}-\frac{\mathbf{d}}{d_{0}}\right)\right]^{-1},\label{eq116} \end{equation} \]

where \(\mathbf{k}\) represents the hydraulic permeability tensor of the solvent through the solid matrix. The permeability and diffusivity tensors are related to the diffusive drag tensors appearing in \eqref{eq109} according to

\[ \begin{equation} \begin{aligned}\mathbf{k} & =\left(\varphi^{w}\right)^{2}\left(\mathbf{f}^{ws}\right)^{-1}\,,\\ \mathbf{d}_{0} & =R\theta\varphi^{w}c\left(\mathbf{f}^{uw}\right)^{-1}\equiv d_{0}\mathbf{I}\,,\\ \mathbf{d} & =R\theta\varphi^{w}c\left(\mathbf{f}^{us}+\mathbf{f}^{uw}\right)^{-1}\,, \end{aligned} \label{eq117} \end{equation} \]

though these explicit relationships are not needed here since \(\mathbf{k}\), \(\mathbf{d}\) and \(d_{0}\) may be directly specified in a particular analysis. Since the axiom of entropy inequality requires that the tensors \(\mathbf{f}^{\alpha\beta}\) be positive semi-definite (see appendix of ¹¹), it follows that \(d_{0}\) must be greater than or equal to the largest eigenvalue of \(\mathbf{d}\). Constitutive relations are needed for these transport properties, which relate them to the solid matrix strain and solute concentration. Note that the relations in \eqref{eq117} represent generalizations of Darcy's law for fluid permeation through porous media, and Fick's law for solute diffusion in porous media or free solution.

Continuous Variables

In principle, the objective of the finite element analysis is to solve for the three unknowns, \(\mathbf{u}\), \(p\) and \(c\), using the partial differential equations that enforce mixture momentum balance in \eqref{eq108}, mixture mass balance in \eqref{eq112}, and solute mass balance in \eqref{eq113}. The remaining solvent and solute momentum balances in \eqref{eq115}, and solid mass balance in \eqref{eq114}, have been reduced to relations that may be substituted into the three partial differential equations as needed. Solving these equations requires the application of suitable boundary conditions that are consistent with mass, momentum and energy balances across boundary surfaces or interfaces. When defining boundaries or interfaces on the solid matrix (the conventional approach in solid mechanics), whose outward unit normal is \(\mathbf{n}\), mass and momentum balance relations demonstrate that the mixture traction \(\mathbf{t}=\boldsymbol{\sigma}\cdot\mathbf{n}\) and normal flux components \(w_{n}=\mathbf{w}\cdot\mathbf{n}\) and \(j_{n}=\mathbf{j}\cdot\mathbf{n}\) must be continuous across the interface ¹⁰¹². Therefore, \(\mathbf{t}\), \(w_{n}\) and \(j_{n}\) may be prescribed as boundary conditions.

Combining momentum and energy balances across an interface also demonstrates that \(\tilde{\mu}^{w}\) and \(\tilde{\mu}^{u}\) must be continuous ¹⁰¹³, implying that these mechano-chemical potentials may be prescribed as boundary conditions. However, because of the arbitrariness of the reference states \(\mu_{0}^{w}\), \(\mu_{0}^{u}\), \(p_{0}\) and \(c_{0}\), and the ill-conditioning of the logarithm function in the limit of small solute concentration, the mechano-chemical potentials do not represent practical choices for primary variables in a finite element implementation. An examination of \eqref{eq110} also shows that continuity of these potentials across an interface does not imply continuity of the fluid pressure \(p\) or solute concentration \(c\). Therefore, pressure and concentration are also unsuitable as nodal variables in a finite element analysis and they must be replaced by alternative choices. Based on the similar reasoning presented by Sun et al. ¹⁴, an examination of the expressions in \eqref{eq110} shows that continuity may be enforced by using

\[ \begin{equation} \begin{aligned}\tilde{p} & =p-R\theta\,\Phi\,c\,,\\ \tilde{c} & =\frac{c}{\tilde{\kappa}}\,, \end{aligned} \label{eq118} \end{equation} \]

where \(\tilde{p}\) is the effective fluid pressure and \(\tilde{c}\) is the effective solute concentration in the mixture. Note that \(\tilde{p}\) represents that part of the fluid pressure which does not result from osmotic effects (since the term \(R\theta\,\Phi\,c\) may be viewed as the osmotic pressure contribution to \(p)\), and \(\tilde{c}\) is a straightforward measure of the solute activity, since \(a^{u}=\tilde{c}/c_{0}\). Therefore these alternative variables have clear physical meanings.

Since the unknowns are now given by \(\mathbf{u}\), \(\tilde{p}\) and \(\tilde{c}\), the governing partial differential equations may be rewritten in the form

\[ \begin{equation} \begin{aligned} & \grad\left(\tilde{p}+R\theta\,\Phi\,\tilde{\kappa}\tilde{c}\right)+\divg\boldsymbol{\sigma}^{e}=\mathbf{0}\,,\\ & \divg\left(\mathbf{v}^{s}+\mathbf{w}\right)=0\,,\\ & \frac{\partial\left(\varphi^{w}\tilde{\kappa}\tilde{c}\right)}{\partial t}+\divg\left(\mathbf{j}+\varphi^{w}\tilde{\kappa}\tilde{c}\,\mathbf{v}^{s}\right)=0\,, \end{aligned} \label{eq119} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{w} & =-\tilde{\mathbf{k}}\cdot\left(\grad\tilde{p}+R\theta\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}\cdot\grad\tilde{c}\right)\,,\\ \mathbf{j} & =\tilde{\kappa}\mathbf{d}\cdot\left(-\varphi^{w}\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}\mathbf{w}\right)\,,\\ \tilde{\mathbf{k}} & =\left[\mathbf{k}^{-1}+\frac{R\theta}{\varphi^{w}}\frac{\tilde{\kappa}\tilde{c}}{d_{0}}\left(\mathbf{I}-\frac{\mathbf{d}}{d_{0}}\right)\right]^{-1}\,. \end{aligned} \label{eq120} \end{equation} \]

Constitutive equations are needed to relate \(\boldsymbol{\sigma}^{e}\), \(\mathbf{k}\), \(\mathbf{d}\), \(d_{0}\), \(\tilde{\kappa}\) and \(\Phi\) to the solid matrix strain and effective solute concentration.

---

1. Truesdell, C.; Toupin, R.. "The classical field theories." Springer, vol. III/1 (1960). ↩

2. Bowen, R.M.. "Theory of mixtures." Continuum physics, vol. 3 (1976). ↩

3. 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). ↩↩↩

4. Ateshian, G. A.; Likhitpanichkul, M.; Hung, C. T.. "A mixture theory analysis for passive transport in osmotic loading of cells." J Biomech, vol. 39, pp. 464-75 (2006). ↩

5. Albro, M. B.; Chahine, N. O.; Li, R.; Yeager, K.; Hung, C. T.; Ateshian, G. A.. "Dynamic loading of deformable porous media can induce active solute transport." J Biomech, vol. 41, pp. 3152-7 (2008). ↩

6. Albro, M. B.; Li, R.; Banerjee, R. E.; Hung, C. T.; Ateshian, G. A.. "Validation of theoretical framework explaining active solute uptake in dynamically loaded porous media." J Biomech, vol. 43, pp. 2267-73 (2010). ↩

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

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

9. 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). ↩

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

11. Ateshian, G. A.; Weiss, J. A.. "Anisotropic hydraulic permeability under finite deformation." Journal of biomechanical engineering, vol. 132, pp. 111004 (2010). ↩

12. Eringen, A.C.; Ingram, J.D.. "Continuum theory of chemically reacting media -- 1." Int J Eng Sci, vol. 3, pp. 197 - 212 (1965). ↩

13. Katzir-Katchalsky, Aharon; Curran, Peter F.. "Nonequilibrium thermodynamics in biophysics." Harvard University Press, pp. x, 248 p. (1965). ↩

14. Sun, D. N.; Gu, W. Y.; Guo, X. E.; Lai, W. M.; Mow, V. C.. "A mixed finite element formulation of triphasic mechano-electrochemical theory for charged, hydrated biological soft tissues." Int J Numer Meth Eng, vol. 45, pp. 1375-402 (1999). ↩