2.15 Hybrid Biphasic Material

In FEBio, the standard biphasic material consists of a mixture of intrinsically incompressible solid and fluid constituents, as described in Section Biphasic Material. In this standard material, the solid and fluid dynamics are neglected, as well as the fluid viscosity. When a user wishes to consider these dynamics effects, as well as fluid viscosity, they may use a hybrid biphasic material, as described in this section. The complete theoretical framework for such materials can be found in ¹. A hybrid biphasic domain is a mixture of an isothermal compressible viscous fluid and a hyperelastic compressible porous solid whose solid skeleton is intrinsically incompressible. Unlike the fluid-FSI material, here the solid has non-negligible mass density and elasticity, and frictional interactions may occur between the fluid and solid constituents, modeled using a hydraulic permeability (Section Hydraulic Permeability). In FEBio, the hybrid biphasic material is called “biphasic-FSI”, since it allows dynamic fluid-structure interactions between this material and a fluid-FSI material, or a solid material. Here, we may abbreviate “biphasic-FSI” as BFSI. For BFSI domains, as with standard biphasic domains, the finite element mesh is defined on the porous solid material. The viscous fluid flows through this mesh, experiencing frictional drag caused by the porous solid. When two BFSI domains are interfaced, or when a BFSI domain is interfaced with a fluid-FSI domain, the pseudo-no slip condition ² is enforced automatically on those interfaces. When a BFSI domain interfaces with a solid domain, the no-slip boundary condition has to be prescribed explicitly on those interfaces, where applicable.

BFSI Governing Equations

We model the BFSI domain as an unconstrained mixture of an isothermal, compressible, and viscous fluid and an isothermal, deformable porous solid whose skeleton material is intrinsically incompressible. As usual in mixture theory, each constituent \(\alpha\) (\(\alpha=s,f\) for solid and fluid, respectively) has its own apparent density \(\rho^{\alpha}=dm^{\alpha}/dV\), which represents the ratio of the elemental mass \(dm^{\alpha}\) of constituent \(\alpha\) in the elemental mixture volume \(dV\). This apparent density is related to the true density \(\rho_{T}^{\alpha}=dm^{\alpha}/dV^{\alpha}\) (mass of \(\alpha\) per volume of \(\alpha\)) via \(\rho^{\alpha}=\varphi^{\alpha}\rho_{T}^{\alpha}\), where \(\varphi^{\alpha}=dV^{\alpha}/dV\) is the volume fraction of \(\alpha\) in the mixture, satisfying \(\varphi^{s}+\varphi^{f}=1\). Since the solid skeleton is intrinsically incompressible, \(\rho_{T}^{s}\) is constant. The boundaries of a biphasic domain are defined on the porous solid matrix. The deformation gradient of the solid is denoted by \(\mathbf{F}^{s}\) and its determinant, \(J^{s}=\det\mathbf{F}^{s}=dV/dV_{r}\), represents the ratio of the mixture elemental volumes in the current (\(dV\)) and reference (\(dV_{r}\)) configurations. Thus, since the solid skeleton is intrinsically incompressible, \(J^{s}\) purely represents the compressibility of the pore volume as fluid enters or leaves the pore space, or as the compressible fluid within the pores changes in volume. The axiom of mass balance for the solid may be integrated in closed form to produce \(\rho^{s}=\rho_{r}^{s}/J^{s}\), where \(\rho_{r}^{s}=dm^{s}/dV_{r}\) is the solid apparent density in the mixture reference configuration. Therefore, we may define the referential solid volume fraction as \(\varphi_{r}^{s}=\rho_{r}^{s}/\rho_{T}^{s}\). In the finite element analysis, \(\varphi_{r}^{s}=dV^{s}/dV_{r}\) and \(\rho_{T}^{s}\) are both specified by the user, then \(\rho_{r}^{s}\), \(\rho^{s}\), \(\varphi^{s}\) and \(\varphi^{f}\) are evaluated from the above relations, given a solution for \(\mathbf{F}^{s}\) (and thus, \(J^{s}\)). The solution for \(\mathbf{F}^{s}=\mathbf{I}+\Grad\mathbf{u}\) is obtained by solving for the nodal solid displacement vector \(\mathbf{u}\), where \(\Grad\left(\cdot\right)\) represents the gradient operator in the material frame.

For the compressible fluid phase of a hybrid biphasic material, the true density \(\rho_{T}^{f}\) varies with the intrinsic fluid volumetric strain, or dilatation, \(e^{f}=J^{f}-1\), where \(J^{f}\) is the volume ratio of the fluid in its current and reference configurations, \(dm^{f}\) is the element mass of fluid and \(dV^{f}\) is the elemental volume of that fluid in the pores of the mixture, in the current configuration. It follows from this definition that

\[ \begin{equation} J^{f}=\frac{dV^{f}}{dV_{r}^{f}}=\frac{dV^{f}}{dm^{f}}\frac{dm^{f}}{dV_{r}^{f}}=\frac{\rho_{Tr}^{f}}{\rho_{T}^{f}}\thinspace,\label{eq:bfsi-Jf} \end{equation} \]

where \(dV_{r}^{f}\) is the elemental volume of this same fluid (having the same elemental mass \(dm^{f}\)) in the fluid's reference configuration. Thus, \(\rho_{Tr}^{f}=dm^{f}/dV_{r}^{f}\) is a constant representing the true fluid density in its reference configuration, which is specified by the user. Accordingly, \(J^{f}=1\) in the limit when the fluid is idealized to be intrinsically incompressible. This definition of \(J^{f}\) is consistent with that for a pure fluid, as shown in our earlier formulation of computational fluid dynamics ³. As shown above, \(\rho^{f}\) may be evaluated as \(\varphi^{f}\rho_{T}^{f}\), which now expands into an expression involving solid and fluid volume ratios,

\[ \begin{equation} \rho^{f}=\left(1-\frac{\varphi_{r}^{s}}{J^{s}}\right)\frac{\rho_{Tr}^{f}}{J^{f}}\,.\label{eq:bfsi-rhof-Jf-Js} \end{equation} \]

The limiting case when the fluid within the solid matrix pores has been completely squeezed out corresponds to \(\rho^{f}=0\) (or \(dm^{f}=0\)). Based on the above relation, this is equivalent to having the entire mixture volume \(dV\) reduce to the solid volume \(dV^{s}\) in the current configuration, or equivalently \(J^{s}=\varphi_{r}^{s}\) ⁴. In this limiting case, the mixture becomes an intrinsically incompressible solid and the finite element formulation presented in this study no longer applies.

We assume that \(\mathbf{F}^{s}=\mathbf{I}\) (or equivalently, \(\mathbf{u}=\mathbf{0}\)) and \(e^{f}=0\) at the start of a finite element analysis. The fluid dilatation \(e^{f}\) is included as a nodal degree of freedom, implying that it is continuous across finite element boundaries. This assumption is verified below, when we review the jump conditions on axioms of mass, momentum and energy balance across interfaces.

Since the axiom of mass balance for the solid has already been solved in closed form, we only need to be concerned with the axiom of mass balance for the fluid, or alternatively, that of the mixture, which takes the form

\[ \begin{equation} \frac{D^{f}J^{f}}{Dt}=\frac{J^{f}}{\varphi^{f}}\divg\left(\mathbf{v}^{s}+\mathbf{w}\right)\thinspace,\label{eq:bfsi-mix-mass-balance-f} \end{equation} \]

where \(\mathbf{v}^{s}\) is the solid velocity (the material time derivative of the solid displacement \(\mathbf{u}\)), and

\[ \begin{equation} \mathbf{w}=\varphi^{f}\left(\mathbf{v}^{f}-\mathbf{v}^{s}\right)\thinspace\label{eq:bfsi-fluid-flux} \end{equation} \]

is the volumetric flux of the fluid relative to the solid, with \(\mathbf{v}^{f}\) representing the fluid velocity. Here, \(D^{f}\left(\cdot\right)/Dt\) is the material time derivative operator in the spatial frame, following the fluid motion. Since the finite element mesh is defined on the porous solid matrix of a biphasic mixture, and since the fluid flows relative to the solid, material time derivatives need to follow the solid motion in the finite element implementation. A similar scheme was used in our implementations of solute transport within deformable porous domains (Sections Triphasic and Multiphasic Materials), to account for the motion of solutes relative to the porous solid matrix, as well as in our FSI formulation (Section Fluid-Structure Interactions) ⁵⁶⁷. Thus, we substitute the following identity,

\[ \begin{equation} \begin{aligned}\frac{D^{f}J^{f}}{Dt} & =\frac{\partial J^{f}}{\partial t}+\grad J^{f}\cdot\mathbf{v}^{f}\\ & =\frac{D^{s}J^{f}}{Dt}+\grad J^{f}\cdot\left(\mathbf{v}^{f}-\mathbf{v}^{s}\right) \end{aligned} \,,\label{eq:bfsi-DfJf-DsJf} \end{equation} \]

into the axiom of mixture mass balance in \eqref{eq:bfsi-mix-mass-balance-f} to produce the final form

\[ \begin{equation} \dot{J}^{f}+\frac{1}{\varphi^{f}}\grad J^{f}\cdot\mathbf{w}=\frac{J^{f}}{\varphi^{f}}\divg\left(\mathbf{v}^{s}+\mathbf{w}\right)\,.\label{eq:bfsi-fluid-kin-constr} \end{equation} \]

In this expression, the dot operator in \(\dot{J}^{f}\) represents the material time derivative following the solid motion. In the solid material frame, this material time derivative reduces to the partial time derivative. Accordingly, we may now write \(\mathbf{v}^{s}=\dot{\mathbf{u}}\). In the BFSI implementation, the fluid volumetric flux \(\mathbf{w}\) relative to the solid is added to the list of nodal DOFs, giving us the complete set \(\left(\mathbf{u},\mathbf{w},e^{f}\right)\). This choice implies that \(\mathbf{w}\) is continuous across finite element boundaries, as justified further below when we review jump conditions across interfaces.

Based on the constitutive assumptions of our hybrid biphasic formulation ¹, the momentum balance equations for the fluid and solid constituents reduce to

\[ \begin{equation} \rho^{f}\mathbf{a}^{f}=-\varphi^{f}\grad p+\divg\boldsymbol{\tau}^{f}+\rho^{f}\mathbf{b}^{f}-\varphi^{f}\mathbf{k}^{-1}\cdot\mathbf{w}\,,\label{eq:bfsi-fluid-momentum} \end{equation} \]

and

\[ \begin{equation} \rho^{s}\mathbf{a}^{s}=-\varphi^{s}\grad p+\divg\boldsymbol{\sigma}^{e}+\rho^{s}\mathbf{b}^{s}+\varphi^{f}\mathbf{k}^{-1}\cdot\mathbf{w}\,,\label{eq:bfsi-solid-momentum} \end{equation} \]

where \(\mathbf{a}^{\alpha}=D^{\alpha}\mathbf{v}^{\alpha}/Dt\) is the acceleration, \(\mathbf{b}^{\alpha}\) is the body force per mass acting on constituent \(\alpha\), \(p\) is the fluid pressure, \(\boldsymbol{\tau}^{f}\) is the apparent fluid viscous stress, and \(\boldsymbol{\sigma}^{e}\) is the apparent solid elastic stress. These stress tensors are called apparent because their associated traction vectors represent a force acting on constituent \(\alpha\) per elemental mixture area. Here, \(\mathbf{k}\) is the hydraulic permeability tensor which regulates frictional drag between the fluid and solid constituents; setting \(\mathbf{k}^{-1}\) to \(\mathbf{0}\) implies that this frictional drag is non-existent. Since the fluid is compressible, its pressure must be given by a function of state. In the isothermal framework used here, this function only depends on \(e^{f}\), and the form adopted in FEBio is

\[ \begin{equation} p=-Ke^{f}\,,\label{eq:bfsi-pressure-state} \end{equation} \]

where \(K\) is the fluid bulk modulus (a user-specified material property). In the limit when inertia and body forces are neglected (\(\mathbf{a}^{f}=\mathbf{b}^{f}=\mathbf{0}\)), the fluid momentum balance \eqref{eq:bfsi-fluid-momentum} produces the classical Darcy-Brinkman relation ⁸. If the fluid viscous stress is also neglected (\(\boldsymbol{\tau}^{f}=\mathbf{0}\)), we recover Darcy's law, \(\mathbf{w}=-\mathbf{k}\cdot\grad p\).

Since the finite element implementation requires all material time derivatives to follow the motion of the solid, we recognize that \(\mathbf{a}^{s}=\ddot{\mathbf{u}}\) and we evaluate the fluid acceleration as \(\mathbf{a}^{f}=\dot{\mathbf{v}}^{f}+\grad\mathbf{v}^{f}\cdot\left(\mathbf{v}^{f}-\mathbf{v}^{s}\right)\). However, since \(\mathbf{v}^{f}\) is not a nodal degree of freedom, we use eq.\eqref{eq:bfsi-fluid-flux} to substitute

\[ \begin{equation} \mathbf{v}^{f}=\frac{1}{\varphi^{f}}\mathbf{w}+\mathbf{v}^{s}\label{eq:bfsi-fluid-velocity} \end{equation} \]

into this expression. It follows that the fluid velocity gradient \(\mathbf{L}^{f}=\grad\mathbf{v}^{f}\) may be evaluated as

\[ \begin{equation} \mathbf{L}^{f}=\mathbf{L}^{s}+\frac{1}{\varphi^{f}}\left(\mathbf{L}^{w}-\frac{1}{\varphi^{f}}\mathbf{w}\otimes\grad\varphi^{f}\right)\,,\label{eq:bfsi-Lf} \end{equation} \]

where \(\mathbf{L}^{w}=\grad\mathbf{w}\) and \(\mathbf{L}^{s}=\grad\mathbf{v}^{s}\). Now, the fluid acceleration takes the form

\[ \begin{equation} \mathbf{a}^{f}=\frac{1}{\varphi^{f}}\left(\dot{\mathbf{w}}+\phi^{f}\ddot{\mathbf{u}}-\frac{\varphi^{s}}{\varphi^{f}}\frac{\dot{J}^{s}}{J^{s}}\mathbf{w}+\mathbf{L}^{f}\cdot\mathbf{w}\right)\,.\label{eq:bfsi-af} \end{equation} \]

BFSI Continuous Variables

Jump conditions on the axioms of mass, momentum and energy balance are needed to determine which variables may be selected as nodal DOFs in the finite element implementation, and which tractions are naturally continuous across an interface. The full set of jump conditions for a hybrid biphasic material were derived in our recent study for the constitutive assumptions adopted in this formulation ¹. Here, we summarize the salient results, which apply to an interface defined on the porous solid matrix of the hybrid biphasic domain, which includes the shared faces of adjoining biphasic elements. Thus, the velocity of the interface is given by the velocity \(\mathbf{v}^{s}\) of the solid constituent of the hybrid biphasic material. We employ the notation \(\left[\left[f\right]\right]=f_{+}-f_{-}\) to denote the jump in the function \(f\) across the interface \(\Gamma\), with \(f_{+}\) and \(f_{-}\) denoting the values of \(f\) on either side of \(\Gamma\). The unit normal on \(\Gamma\) is \(\mathbf{n}\), which points away from the \(+\) side. A variable \(f\) which is continuous across \(\Gamma\) satisfies \(\left[\left[f\right]\right]=0\).

Based on the jump condition on the axiom of mass balance, the normal component of the mass flux of the fluid relative to the solid is continuous across \(\Gamma\), \(\left[\left[\rho_{T}^{f}\mathbf{w}\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 fluid specific free enthalpy (also known as the Gibbs function), \(\left[\left[\psi^{f}+p/\rho_{T}^{f}\right]\right]=0\), where \(\psi^{f}\) is the fluid specific free energy. This jump condition applies only when there is fluid on both sides of the interface \(\Gamma\). In an isothermal framework the specific free enthalpy is a function of state that only depends on \(J^{f}\), therefore this energy jump condition implies that \(J^{f}\) must be continuous across \(\Gamma\), thus

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

also implying that \(\left[\left[p\right]\right]=0\). Given eq.\eqref{eq:bfsi-Jf}, it follows that \(\left[\left[\rho_{T}^{f}\right]\right]=0\) and the mass balance jump condition reduces to \(\left[\left[\mathbf{w}\right]\right]\cdot\mathbf{n}=0\), implying that the relative fluid flux component normal to \(\Gamma\) must be continuous. For the tangential component of \(\mathbf{w}\) on \(\Gamma\) we appeal to the analysis of Hou et al. ², who showed that a valid pseudo-noslip condition requires this tangential component to be continuous. Combining these two jump conditions produces

\[ \begin{equation} \left[\left[\mathbf{w}\right]\right]=0\,.\label{eq:bfsi-fluid-flux-jump} \end{equation} \]

The momentum jump condition requires that the mixture traction be continuous across \(\Gamma\), thus \(\left[\left[-p\mathbf{I}+\boldsymbol{\sigma}^{e}+\boldsymbol{\tau}^{f}\right]\right]\cdot\mathbf{n}=\mathbf{0}\). Since \(\left[\left[p\right]\right]=0\) based on the energy jump, this mixture momentum jump condition reduces to

\[ \begin{equation} \left[\left[\boldsymbol{\sigma}^{e}+\boldsymbol{\tau}^{f}\right]\right]\cdot\mathbf{n}=0\,.\label{eq:bfsi-mix-mtm-jump} \end{equation} \]

Finally, another relation which is sufficient to satisfy the jump condition on the energy balance is the continuity of the true fluid traction (force acting on fluid per fluid area),

\[ \begin{equation} \left[\left[\frac{\boldsymbol{\tau}^{f}}{\varphi^{f}}\right]\right]\cdot\mathbf{n}=0\,.\label{eq:bfsi-fluid-mtm-jump} \end{equation} \]

This jump condition eq.\eqref{eq:bfsi-fluid-mtm-jump}, which also applies only if fluid is present on both sides of \(\Gamma\), is interesting because it implies that the viscous stress (and thus, the viscosity) of a fluid flowing in a porous solid matrix scales with the porosity of that medium, such that \(\boldsymbol{\tau}^{f}=\phi^{f}\boldsymbol{\tau}\) where \(\boldsymbol{\tau}\) would be the true fluid viscous stress. Thus, we can use FEBio's various constitutive relations for the viscous stress \(\boldsymbol{\tau}\) of Newtonian or non-Newtonian fluids and adapt those models to a biphasic mixture where \(\boldsymbol{\tau}^{f}\) is evaluated as \(\varphi^{f}\boldsymbol{\tau}\). Accordingly, the contribution of \(\boldsymbol{\tau}^{f}\) would properly reduce to zero in the limit as fluid content reduces to zero (\(\varphi^{f}\to0\)), in which case the mixture momentum jump in eq.\eqref{eq:bfsi-mix-mtm-jump} would reduce to \(\left[\left[-p\mathbf{I}+\boldsymbol{\sigma}^{e}\right]\right]\cdot\mathbf{n}=\left[\left[\boldsymbol{\sigma}^{e}\right]\right]\cdot\mathbf{n}=0\), since \(\left[\left[p\right]\right]=0\).

Letting \(\mathbf{w}\) and \(J^{f}\) be nodal DOFs automatically enforces the jump conditions eq.\eqref{eq:bfsi-Jf-jump} and eq.\eqref{eq:bfsi-fluid-flux-jump}, acting as essential boundary conditions, along with the solid displacement \(\mathbf{u}\). Finally, subtracting eq.\eqref{eq:bfsi-fluid-mtm-jump} from eq.\eqref{eq:bfsi-mix-mtm-jump}, we obtain the momentum jump condition for the solid constituent,

\[ \begin{equation} \left[\left[-\frac{\varphi^{s}}{\varphi^{f}}\boldsymbol{\tau}^{f}+\boldsymbol{\sigma}^{e}\right]\right]\cdot\mathbf{n}=0\,.\label{eq:bfsi-solid-mtm-jump} \end{equation} \]

A BFSI domain may be reduced to a FSI domain by letting \(\varphi_{r}^{s}\to0\) and \(\mathbf{k}^{-1}\to\mathbf{0}\); the solid matrix would still be ascribed a material response but its stiffness would need to be negligible. The number of nodal DOFs would remain the same. The CFD domain is a special case of the FSI domain where the mesh displacement is uniformly \(\mathbf{u}=\mathbf{0}\).

However, when the biphasic domain interfaces with a non-porous solid domain across \(\Gamma^{bs}\), the jump conditions eq.\eqref{eq:bfsi-Jf-jump} on \(J^{f}\) and eq.\eqref{eq:bfsi-solid-mtm-jump} on \(\boldsymbol{\tau}^{f}\) don't apply. In that case, the jump condition eq.\eqref{eq:bfsi-fluid-flux-jump} on the relative fluid volumetric flux should reduce to \(\mathbf{w}=\mathbf{0}\) for the BFSI domain on \(\Gamma^{bs}\), although the user would have to enforce this no-slip condition explicitly by prescribing it as an essential boundary condition. The mixture momentum jump \(\left[\left[-p\mathbf{I}+\boldsymbol{\sigma}^{e}+\boldsymbol{\tau}^{f}\right]\right]\cdot\mathbf{n}=\mathbf{0}\) implies that the mixture traction on the BFSI side of an interface with a solid is equal to the solid traction on the solid side. This jump condition would also need to be enforced explicitly.

---

1. Shim, Jay J; Ateshian, Gerard A. "A Hybrid Biphasic Mixture Formulation for Modeling Dynamics in Porous Deformable Biological Tissues." Arch Appl Mech, vol. 92, pp. 491-511 (2022). ↩↩↩

2. Hou, J S; Holmes, M H; Lai, W M; Mow, V C. "Boundary conditions at the cartilage-synovial fluid interface for joint lubrication and theoretical verifications." J Biomech Eng, vol. 111, pp. 78-87 (1989). ↩↩

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

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

5. Ateshian, Gerard A; Maas, Steve; Weiss, Jeffrey A. "Solute transport across a contact interface in deformable porous media." J Biomech, vol. 45, pp. 1023-7 (2012). ↩

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

7. Shim, Jay J; Maas, Steve A; Weiss, Jeffrey A; Ateshian, Gerard A. "A Formulation for Fluid Structure-Interactions in FEBio Using Mixture Theory." J Biomech Eng (2019). ↩

8. Brinkman, Hendrik C. "A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles." Flow, Turbulence and Combustion, vol. 1, pp. 27--34 (1949). ↩