Warning: MathJax requires JavaScript to correctly process the mathematics on this page. Please enable JavaScript on your browser.

Prev Section 2.13: Hybrid Biphasic Material Up Section 2.13: Hybrid Biphasic Material Subsection 2.13.2: BFSI Continuous Variables Next

### 2.13.1 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}},\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 [49]. 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} [8]. 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),\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)\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 2.7↑), to account for the motion of solutes relative to the porous solid matrix, as well as in our FSI formulation (Section 2.12↑) [11, 52, 53]. 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 (2.13.1-3) 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 [90], 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 \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 (2.13.1-7) produces the classical Darcy-Brinkman relation [55]. 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 .

*mixture*area. Here, 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.(2.13.1-4) 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}