Abstract

In biomechanics, solid–fluid mixtures have commonly been used to model the response of hydrated biological tissues. In cartilage mechanics, this type of mixture, where the fluid and solid constituents are both assumed to be intrinsically incompressible, is often called a biphasic material. Various physiological processes involve the interaction of a viscous fluid with a porous-hydrated tissue, as encountered in synovial joint lubrication, cardiovascular mechanics, and respiratory mechanics. The objective of this study was to implement a finite element solver in the open-source software febio that models dynamic interactions between a viscous fluid and a biphasic domain, accommodating finite deformations of both domains as well as fluid exchanges between them. For compatibility with our recent implementation of solvers for computational fluid dynamics (CFD) and fluid–structure interactions (FSI), where the fluid is slightly compressible, this study employs a novel hybrid biphasic formulation where the porous skeleton is intrinsically incompressible but the fluid is also slightly compressible. The resulting biphasic-FSI (BFSI) implementation is verified against published analytical and numerical benchmark problems, as well as novel analytical solutions derived for the purposes of this study. An illustration of this BFSI solver is presented for two-dimensional (2D) airflow through a simulated face mask under five cycles of breathing, showing that masks significantly reduce air dispersion compared to the no-mask control analysis. In addition, we model three-dimensional (3D) blood flow in a bifurcated carotid artery assuming porous arterial walls and verify that mass is conserved across all fluid–permeable boundaries. The successful formulation and implementation of this BFSI solver offers enhanced multiphysics modeling capabilities that are accessible via an open-source software platform.

1 Introduction

The fields of biomechanics and biophysics increasingly necessitate more sophisticated computational tools to describe complex physiological processes that span across the molecular, cellular, tissue, organ, and whole-body levels. Mixture theory is a general framework for modeling mixtures of solid, fluid, and solute constituents [14], leading to significant advances in modeling the mechanics of biological tissues and cells [59]. The mixture framework was first used to represent biological tissues modeled as a porous solid that was permeable to an interstitial fluid, where each constituent was intrinsically incompressible, though the mixture could lose or gain volume via fluid exchange with the pore domain [5,10,11]. This mixture formulation has been called the biphasic theory in the field of cartilage biomechanics. Bowen [11] demonstrated that this biphasic material reduced to Biot's incompressible poroelasticity theory under certain simplifying assumptions [12]. Biphasic theory was successfully used in many quasi-static applications including modeling the mechanics of articular cartilage [5,1316], brain [17], lung [18], liver [19], heart [20,21], cells and their extracellular matrix [2224], and the eye [2527], among many others. The increasing use of biphasic materials led to the development of related computational tools [2833].

The most commonly used biphasic material formulation developed by Mow and colleagues made a number of assumptions that simplified its usage when modeling biological tissues [5]. For standard biphasic materials, fluid and solid constituents were assumed to be intrinsically incompressible, intrinsic fluid viscosity was negligible, and only quasi-static analyses were typically examined. Neglecting inertia terms and intrinsic fluid viscosity also reduced the nodal degrees-of-freedom (DOFs) to the solid displacement vector and fluid pressure, eliminating the fluid velocity. Frictional drag between the fluid and the porous solid skeleton gave the mixture a hydraulic permeability that dictated the fluid flow solely based on the fluid pressure gradient, a relation known as Darcy's law in classical porous media mechanics. These models effectively assumed that viscous dissipation in the fluid was negligible compared to that produced by the frictional drag between the fluid and solid. Since both constituents were assumed to be intrinsically incompressible, the fluid pressure appearing in biphasic theory is a Lagrange multiplier that enforces this incompressibility constraint. Therefore, while the classical biphasic material has been used successfully to model many types of tissues, assumptions such as inviscid interstitial fluid, quasi-static analyses, and intrinsic incompressibility may not apply to all situations.

Various investigators have extended the classical biphasic theory to generalize its application to a variety of modeling situations. For example, intrinsic fluid viscosity has been included in the linear biphasic formulation by Hou et al. [34], to model quasi-static fluid transport between a biphasic domain and a surrounding viscous fluid body. Later, finite element codes were developed that implemented such analyses [3537]. These models have been used in applications such as arterial mechanics [38], cerebrospinal fluid hydrodynamics [39], lung respiratory mechanics [40], and articular cartilage squeeze-film lubrication [41].

In principle, existing frameworks for mixtures of intrinsically incompressible constituents already provide the necessary foundation for modeling a biphasic domain with dynamic effects [5,11] and fluid viscosity [42]. However, since the fluid is incompressible in those formulations, wave propagation in that medium has infinite speed, which places a limitation on the range of physical phenomena that can be modeled realistically. Alternatively, mixture formulations where all the constituents are compressible [43] require internal variables and evolution functions to distinguish the volumetric strains in each of the mixture constituents when subjected to hydrostatic pressure, which would add undesirable complexity to our formulation.

In this study, we extend these prior investigations by formulating and implementing a more general finite element framework where the biphasic material has a viscous Newtonian or non-Newtonian interstitial fluid, dynamic effects in the fluid and solid, and a hyperelastic or viscoelastic porous solid matrix that is undergoing finite deformation. In a recent study [44], we modified the standard biphasic theory to make the fluid compressible, while the porous–permeable solid skeleton remained intrinsically incompressible. This hybrid biphasic formulation allows us to model biological tissues that interface with a viscous fluid domain under dynamic conditions, accommodating fluid exchanges between these domains. This hybrid biphasic theory thus facilitates the analysis of a wider range of physiological processes that were previously difficult to model, such as dynamic interactions between synovial fluid and permeable cartilage, blood and permeable arterial or cardiac walls, cerebrospinal fluid and permeable brain, or air and permeable lungs. We may describe these interactions generically as fluid–structure interactions (FSI), where the structural domain is biphasic, which we call biphasic-FSI (BFSI). We implement this framework in the open-source finite element software febio1 to facilitate its dissemination and model sharing [33,45].

The finite element implementation of this framework represents a natural extension of our recently reported computational fluid dynamics (CFD) and FSI solvers [46,47]. In the CFD solver, the mesh is stationary and the viscous fluid flows through it. In the FSI solver, the mesh is deformable and the viscous fluid flows through it; the deformability of the mesh is regularized by assigning it a solid matrix with negligible (but nonzero) elastic stiffness, zero mass density, and zero solid–fluid frictional drag. Interactions of this deformable fluid domain with surrounding solid structures are achieved by implementing suitable boundary conditions at interfaces between those domains. In the BFSI formulation proposed here, the fluid similarly flows through a deformable mesh, which now represents the porous solid matrix of a hybrid biphasic material, whose material response, mass density, and frictional momentum exchange with the interstitial fluid are specified by the user. Hybrid biphasic domains may interact with solid and fluid domains that may be deformable. Throughout this presentation, we refer to stationary fluid, deformable fluid, and deformable hybrid biphasic domains as CFD, FSI, and BFSI domains, respectively.

2 Finite Element Implementation

The BFSI domain, denoted by Ωb, is a fluid–solid mixture that can be viewed as a generalization of the standard biphasic, CFD, and FSI domains. This hybrid biphasic domain may interface with other hybrid biphasic domains, fluid domains Ωf, and solid domains Ωs across interfaces Γbb,Γbf, and Γbs, respectively (Fig. 1). 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 equations presented below address the analysis of BFSI domains Ωb. A more detailed derivation of the hybrid biphasic theory describing the response of the BFSI domain can be found in a previous study [44]. The finite element analyses of solid, CFD, and FSI domains were described previously [33,46,47].

Fig. 1
Schematic description of hybrid biphasic domains (Ωb) interfacing with each other and with a fluid or FSI domain (Ωf) and a solid domain (Ωs) via interfaces Γbb, Γbf, and Γbs. For example, Ωf may represent the blood flowing within an artery, the inner Ωb may represent the tunica media, the outer Ωb may represent a portion of the tunica adventitia, and Ωs may represent some fatty tissue.
Fig. 1
Schematic description of hybrid biphasic domains (Ωb) interfacing with each other and with a fluid or FSI domain (Ωf) and a solid domain (Ωs) via interfaces Γbb, Γbf, and Γbs. For example, Ωf may represent the blood flowing within an artery, the inner Ωb may represent the tunica media, the outer Ωb may represent a portion of the tunica adventitia, and Ωs may represent some fatty tissue.
Close modal

2.1 Governing Equations.

We model the domain Ωb 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 α (α=s,f for solid and fluid, respectively) has its own apparent density ρα=dmα/dV, which represents the ratio of the elemental mass dmα of constituent α in the elemental mixture volume dV. This apparent density is related to the true density ρTα=dmα/dVα (mass of α per volume of α) via ρα=ϕαρTα, where ϕα=dVα/dV is the volume fraction of α in the mixture, satisfying ϕs+ϕf=1. Since the solid skeleton is intrinsically incompressible, ρTs is constant. The boundaries of a biphasic domain are defined on the porous solid matrix. The deformation gradient of the solid is denoted by Fs and its determinant, Js=detFs=dV/dVr, represents the ratio of the mixture elemental volumes in the current (dV) and reference (dVr) configurations. Thus, since the solid skeleton is intrinsically incompressible, Js 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 ρs=ρrs/Js, where ρrs=dms/dVr is the solid apparent density in the mixture reference configuration. Therefore, we may define the referential solid volume fraction as ϕrs=ρrs/ρTs. In the finite element analysis, ϕrs=dVs/dVr and ρTs are both specified by the user, then ρrs, ρs, ϕs, and ϕf are evaluated from the above relations, given a solution for Fs (and thus, Js). The solution for Fs=I+Gradu is obtained by solving for the nodal solid displacement vector u, where Grad(·) represents the gradient operator in the material frame.

For the compressible fluid phase of a hybrid biphasic material, the true density ρTf varies with the intrinsic fluid volumetric strain, or dilatation, ef=Jf1, where Jf is the volume ratio of the fluid in its current and reference configurations, dmf is the element mass of fluid, and dVf is the elemental volume of that fluid in the pores of the mixture, in the current configuration. It follows from this definition that:
Jf=dVfdVrf=dVfdmfdmfdVrf=ρTrfρTf
(2.1)
where dVrf is the elemental volume of this same fluid (having the same elemental mass dmf) in the fluid's reference configuration. Thus, ρTrf=dmf/dVrf is a constant representing the true fluid density in its reference configuration, which is specified by the user. Accordingly, Jf=1 in the limit when the fluid is idealized to be intrinsically incompressible. This definition of Jf is consistent with that for a pure fluid, as shown in our earlier formulation of computational fluid dynamics [46]. As shown above, ρf may be evaluated as ϕfρTf, which now expands into an expression involving solid and fluid volume ratios
ρf=(1ϕrsJs)ρTrfJf
(2.2)

The limiting case when the fluid within the solid matrix pores has been completely squeezed out corresponds to ρf=0 (or dmf=0). Based on the above relation, this is equivalent to having the entire mixture volume dV reduce to the solid volume dVs in the current configuration, or equivalently Js=ϕrs [48]. In this limiting case, the mixture becomes an intrinsically incompressible solid and the finite element formulation presented in this study no longer applies. Other investigators have proposed computational methods to account for this limiting case [49], though these are not contemplated in this study.

We assume that Fs=I (or equivalently, u=0) and ef=0 at the start of a finite element analysis. The fluid dilatation ef 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
DfJfDt=Jfϕfdiv(vs+w)
(2.3)
where vs is the solid velocity (the material time derivative of the solid displacement u), and
w=ϕf(vfvs)
(2.4)
is the volumetric flux of the fluid relative to the solid [12], with vf representing the fluid velocity; some authors also call it the filtration velocity [49]. Here, Df(·)/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 previous implementations of solute transport within deformable porous domains, to account for the motion of solutes relative to the porous solid matrix, as well as in our FSI formulation [47,50,51]. Thus, we substitute the following identity:
DfJfDt=Jft+gradJf·vf=DsJfDt+gradJf·(vfvs)
(2.5)
into the axiom of mixture mass balance in Eq. (2.3) to produce the final form
J˙f+1ϕfgradJf·w=Jfϕfdiv(vs+w)
(2.6)

In this expression, the dot operator in 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 vs=u˙. In the BFSI implementation, the fluid volumetric flux w relative to the solid is added to the list of nodal DOFs, giving us the complete set (u,w,ef). This choice implies that 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 [44], the momentum balance equations for the fluid and solid constituents reduce to
ρfaf=ϕfgradp+divτf+ρfbfϕfk1·w
(2.7)
and
ρsas=ϕsgradp+divσe+ρsbs+ϕfk1·w
(2.8)
where aα=Dαvα/Dt is the acceleration, bα is the body force per mass acting on constituent α, p is the fluid pressure, τf is the apparent fluid viscous stress, and σe is the apparent solid elastic stress. These stress tensors are called apparent because their associated traction vectors represent a force acting on constituent α per elemental mixture area. Here, k is the hydraulic permeability tensor which regulates frictional drag between the fluid and solid constituents; setting k1 to 0 implies that this frictional drag is nonexistent. 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 ef, and its simplest form is
p=Kef
(2.9)

where K is the fluid bulk modulus (a user-specified material property). In the limit when inertia and body forces are neglected (af=bf=0), the fluid momentum balance (2.7) produces the classical Darcy–Brinkman relation [52]. If the fluid viscous stress is also neglected (τf=0), we recover Darcy's law, w=k·gradp.

Since the finite element implementation requires all material time derivatives to follow the motion of the solid, we recognize that as=ü and we evaluate the fluid acceleration as af=v˙f+gradvf·(vfvs). However, since vf is not a nodal degree-of-freedom, we use Eq. (2.4) to substitute
vf=1ϕfw+vs
(2.10)
into this expression. It follows that the fluid velocity gradient Lf=gradvf may be evaluated as:
Lf=Ls+1ϕf(Lw1ϕfwgradϕf)
(2.11)
where Lw=gradw and Ls=gradvs. Now, the fluid acceleration takes the form
af=1ϕf(w˙+ϕfüϕsϕfJ˙sJsw+Lf·w)
(2.12)

2.2 Jump Conditions.

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 [44]. 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, Γbb, and the interfaces Γbf and Γbs. Thus, the velocity of the interface Γ is given by the velocity vs of the solid constituent of the biphasic material on Γ. We employ the notation [[f]]=f+f to denote the jump in the function f across the interface Γ, with f+ and f denoting the values of f on either side of Γ. The unit normal on Γ is n, which points away from the + side. A variable f which is continuous across Γ satisfies [[f]]=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 Γ, [[ρTfw]]·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), [[ψf+p/ρTf]]=0, where ψf is the fluid specific free energy. This jump condition applies only when there is fluid on both sides of the interface Γ. In an isothermal framework, the specific free enthalpy is a function of state that only depends on Jf; therefore, this energy jump condition implies that Jf must be continuous across Γ, thus
[[Jf]]=0
(2.13)
also implying that [[p]]=0. Given Eq. (2.1), it follows that [[ρTf]]=0 and the mass balance jump condition reduces to [[w]]·n=0, implying that the relative fluid flux component normal to Γ must be continuous. For the tangential component of w on Γ, we appeal to the analysis of Hou et al. [34], who showed that a valid pseudo-noslip condition requires this tangential component to be continuous. Combining these two jump conditions produces
[[w]]=0
(2.14)
The momentum jump condition requires that the mixture traction be continuous across Γ; thus, [[pI+σe+τf]]·n=0. Since [[p]]=0 based on the energy jump, this mixture momentum jump condition reduces to
[[σe+τf]]·n=0
(2.15)
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)
[[τfϕf]]·n=0
(2.16)

This jump condition (2.16), which also applies only if fluid is present on both sides of Γ, 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 τf=ϕfτ where τ would be the true fluid viscous stress. Thus, we can use febio's various constitutive relations for the viscous stress τ of Newtonian or non-Newtonian fluids and adapt those models to a biphasic mixture where τf is evaluated as ϕfτ. Accordingly, the contribution of τf would properly reduce to zero in the limit as fluid content reduces to zero (ϕf0), in which case the mixture momentum jump in Eq. (2.15) would reduce to [[pI+σe]]·n=[[σe]]·n=0, since [[p]]=0.

Letting w and Jf be nodal DOFs automatically enforces the jump conditions (2.13) and (2.14), acting as essential boundary conditions, along with the solid displacement u. Finally, subtracting Eq. (2.16) from Eq. (2.15), we obtain the momentum jump condition for the solid constituent
[[ϕsϕfτf+σe]]·n=0
(2.17)

A BFSI domain may be reduced to a FSI domain by letting ϕrs0 and k10; 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 u=0. Therefore, essential and natural jump conditions satisfied across Γbb are equally valid for interfaces Γbf.

However, when the biphasic domain interfaces with a nonporous solid domain across Γbs, the jump conditions (2.13) on Jf and (2.17) on τf do not apply. In that case, the jump condition (2.14) on the relative fluid volumetric flux should reduce to w=0 for the BFSI domain on Γ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 [[pI+σe+τf]]·n=0 implies that the mixture traction on the Ωb side of Γbs is equal to the solid traction on the Ωs side. This jump condition would also need to be enforced explicitly, as described below.

2.3 Virtual Work and Weak Form.

The virtual work statement is used to enforce the three governing equations needed to solve for the nodal DOFs u, w, and ef, namely, the mixture mass balance (2.6), the fluid momentum balance (2.7), and the solid momentum balance (2.8). We may rewrite the momentum balance equations to facilitate the enforcement of natural traction boundary conditions given in Eqs. (2.16) and (2.17). Using τf=ϕfτ, these become
ϕs(ρTfaf+ρTsas)=ϕs2ϕfτ·gradϕfϕs+div(ϕsτ+σe)+ϕs(ρTfbf+ρTsbs)+k1·wρTfaf=gradp+1ϕfτ·gradϕf+divτ+ρTfbfk1·w
(2.18)
The virtual work statement for a Galerkin finite element formulation [53] is δW=0, where
δW=Ωbδvs·[div(ϕsτ+σe)ϕs2ϕfτ·gradϕfϕs+k1·w]dv+Ωbδvs·ϕs(ρTf(bfaf)+ρTs(bsas))dv+Ωbδw·[divτgradp+1ϕfτ·gradϕfk1·w]dv+Ωbδw·ρTf(bfaf)dv+ΩbδJf[divw+J˙sJs1Jf(ϕfJ˙f+gradJf·w)]dv
(2.19)
These integrals are evaluated in the current configuration of Ωb. Here, δvs is the virtual solid velocity, δw is the virtual relative fluid volumetric flux, and δJf is the virtual fluid energy density. Integrating by parts and using the divergence theorem, the weak form of this statement may be written as δW=δWextδWint where the internal virtual work is
δWint=Ωb[(ϕsτ+σe):gradδvsδvs·k1·w]dv+Ωbδvs·ϕs(ϕsϕfτ·gradϕfϕsρTfaf+ρTsas)dv+Ωbτ:gradδwdv+Ωbδw·(gradp+k1·w)dv+Ωbδw·(ρTfaf1ϕfτ·gradϕf)dv+Ωbw·gradδJfdv+ΩbδJf[1Jf(ϕfJ˙f+gradJf·w)J˙sJs]dv
(2.20)
and the external part is
δWext=Ωbδvs·tσda+Ωbδvs·ϕs(ρTfbf+ρTsbs)dv+Ωbδw·tτda+Ωbδw·ρTfbfdv+ΩbδJfwnda
(2.21)
where
tσ=ϕstτ+te
(2.22)

Here, te=σe·n is the elastic traction, tτ=τ·n is the true fluid viscous traction, and wn=wn·n is the normal component of the relative fluid flux on the boundary Ωb, whose outward unit normal is n. The traction tσ emerges from the jump condition in Eq. (2.17). The integrands of the surface integrals represent the natural boundary conditions for this formulation. If boundary conditions are not set explicitly on Ωb, the natural boundary conditions are tσ=0, tτ=0, and wn = 0. These natural boundary conditions are consistent with the jump conditions presented above. Essential boundary conditions are prescribed on the solid displacement u, relative volumetric fluid flux w, and fluid dilatation ef, which are also consistent with the above jump conditions. In particular, an essential no-slip boundary condition may be prescribed on Γbs by setting w=0. A symmetry plane may be prescribed with the essential boundary condition unu·n=0 and the natural boundary conditions tτ=0 and wn = 0.

In this formulation, the mixture traction is defined as t=pn+te+ϕftτ, which may also be written as t=pn+tσ+tτ. Because of the way we chose to split the internal and external virtual work in Eqs. (2.20) and (2.21), t is not a natural boundary condition in this formulation. In this expression for t, tσ and tτ may be prescribed as natural boundary conditions, whereas p may be prescribed as an essential boundary condition on ef, using Eq. (2.9). However, there are two general scenarios where t needs to be prescribed on a region of the boundary Ωb with incomplete prior knowledge of p, tσ, or tτ: (1) When a BFSI boundary Γb represents a free surface (such as the fluid surface in channel flow), the mixture traction boundary condition requires that t=0, in which case it is necessary to explicitly enforce pn+tσ+tτ=0 as a constraint equation on that boundary, to impart the free surface Γb its natural shape. (2) At a biphasic–solid interface Γbs, t must balance the traction acting on the adjoining solid domain. Since u is continuous across Γbs due to shared nodes, the solid natural boundary condition tσ of t is already accounted for by the deformation, so that it is only necessary to prescribe the portion pn+tτ of t on the solid domain; thus, pn+tτ+ts=0 where ts is the (equal and opposite) traction acting on the solid domain. In both cases, the form of this traction boundary condition is the same, with te and ts representing the tractions acting on Γb and Γbs, respectively.

For both of these cases, the resulting virtual work on the free surface Γb or the interface Γbs takes the form
δF=Γδvs·(pn+tτ)da
(2.23)
where the elemental area da on Γ may be evaluated from the covariant basis vectors gα (α=1,2)
da=|g1×g2|dη1dη2
(2.24)
where
gα=x(η1,η2)ηα
(2.25)
and x(η1,η2) is the parametric representation of Γ, defined on the solid constituent. The outward normal n to Ωb on Γ is evaluated from
n=g1×g2|g1×g2|
(2.26)
As a result, the virtual work can be rewritten as
δF=Γδvs·(pI+τ)·(g1×g2)dη1dη2
(2.27)

In febio, this boundary condition is called BFSI traction, which the user must explicitly prescribe on free surfaces Γb and deformable interfaces Γbs. The code automatically determines which of these two types of interfaces is being considered.

2.4 Linearization and Spatial Discretization.

The solution to the nonlinear equation δW=0 is obtained by linearizing this relation along increments in the DOFs as
δW+DδW[Δu]+DδW[Δw]+DδW[Δef]0
(2.28)

where the operator DδW[·] represents the directional derivative of δW at (u,w,Jf) along increments Δu of u, Δw of w, and Δef of ef [53]. The details of this linearization, as well as the spatial discretizations of δW and its directional derivatives, are presented in the febio Theory Manual.1 The temporal discretizations of the solid velocity and solid and fluid accelerations use the generalized-α scheme that was adapted previously for CFD and FSI domains [46,47,54,55]. The linearization of the BFSI interface traction, along with specialized boundary conditions, which were defined for the CFD and fluid-FSI domain in our previous study [46], such as backflow and tangential stabilization and flow resistance, are also adapted to the BFSI formulation as shown in the Theory Manual.

2.5 Solver.

The finite element solver previously used for FSI domains is employed for BFSI domains as well [47]. The numerical solution of the overall nonlinear system of equations is obtained using Newton's method for the first iteration of a time point, followed by Broyden quasi-Newton updates [56], until convergence is achieved for that time-step. As a result, the stiffness matrix is only evaluated once for that time-step. Optionally, the stiffness matrix evaluation and Newton update may be postponed at the start of subsequent time-steps until Broyden updates exceed a user-defined value, such as 50 updates, to offer considerable numerical efficiency. Relative convergence is achieved at each discrete time tn when the vector ΔU(k) of nodal DOF increments, consisting of Δu,Δw, and Δef values at all unconstrained nodes, at the kth iteration satisfies ||ΔU(k)||εr||ΔU(1)||, where εr represents the relative tolerance criterion, which is typically set to εr=103 for all DOF increments. An absolute convergence criterion is also set based on the magnitude of the residual vector. For Newton updates, febio uses a variety of linear equation solvers, among the most efficient is the pardiso parallel direct sparse solver for the solution of a linear system of equations with nonsymmetric coefficient matrix [57].2 Results reported in this study employed the Intel Math Kernel Library implementation of pardiso.3febio uses openmp to parallelize finite element calculations to multiple processor cores and can efficiently be solved with high performance computing clusters.

3 Verification and Validation

We analyzed a variety of test problems to ensure that the BFSI solver produced expected results. We compared models using our previous CFD, FSI, and biphasic solvers with those obtained from limiting cases of the BFSI solver, as a means of verifying the code. In particular, we looked at one-dimensional (1D) confined compression of a biphasic material where the BFSI fluid viscosity was set to zero, two-dimensional (2D) unsteady flow past a cylinder using CFD where the BFSI solid displacements were fixed and its solid volume fraction and diffusive drag k1 were set to zero, and bubble inflation using FSI where the BFSI solid matrix stiffness was also prescribed to be negligible. We also examined fluid–biphasic interface benchmark problems described in the literature [34,35,37,58], which have steady-state analytical solutions for the case of a standard biphasic material subjected to infinitesimal strains with a Newtonian viscous interstitial fluid and isotropic constant permeability. These models assessed our BFSI solver's ability to capture continuity conditions along the tangential direction of interface Γbf under steady-state conditions. We also obtained a novel analytical solution for transient 1D ultrafiltration of a Newtonian viscous fluid through a hybrid biphasic material and verified the BFSI solver against that solution. This model assessed the BFSI solver's ability to capture continuity conditions along the normal direction of interface Γbf under quasi-static transient conditions. Then, we examined if the BFSI code could predict two dilatational wave propagation speeds in a 1D analysis, as predicated by the compressibility of the porous solid and interstitial fluid. This wave propagation analysis was used to verify the BFSI solver's ability to capture dynamic phenomena in the fluid and solid constituents of a hybrid biphasic material. In the next example, 2D airflow through a thin, porous BFSI barrier with different hydraulic permeabilities was tested to simulate breathing through a face mask or in the absence of one. This problem served as an illustration of configurations that cannot be easily modeled using previously existing computational methods. Finally, we investigated three-dimensional (3D) blood flow within a bifurcated carotid artery with porous BFSI walls, where we compared realistic and exaggeratedly porous arterial walls. This example showed that the BFSI solver can be used to model realistic 3D biological tissues with nonlinear and anisotropic constitutive models, and that mass is conserved in both cases, even with significant fluid leaving through the wall in the more porous case.

3.1 Biphasic One-Dimensional Confined Compression.

This 1D analysis verified that the BFSI domain reduces to the standard biphasic formulation [5]. A block of height h=1mm along the z-direction (0zh) was subjected to a prescribed linear ramp-and-hold displacement at z=1mm, producing a compression of 10% in 100s, and then allowed to relax to time t=1000s. The fluid pressure was set to zero at z = h to simulate loading with a free-draining rigid porous filter. The displacement was fixed at z =0 and on the lateral surfaces of the block. The fluid flux normal to those surfaces was naturally equal to zero. The material was set as either a standard biphasic or a BFSI material, with ϕrs=0.2. The porous solid matrix was a compressible neo-Hookean material [53] with E=1kPa, and ν = 0. The hydraulic permeability tensor was isotropic, k=kI, with constant k=0.001mm4/mNs. For the BFSI material, the fluid constituent had a bulk modulus K=1×1012kPa and zero shear and bulk viscosities, making the fluid inviscid and nearly incompressible to reproduce the behavior of the standard biphasic material. The mesh consisted of ten hexahedral eight-node elements in the z-direction and was biased near the porous filter to capture fluid pressure gradients. A time-step of Δt=2.5s was used for the analysis. The mesh had 44 nodes, with 80 DOFs for the standard biphasic material and 212 DOFs for the BFSI material. The comparison of axial displacements uz and fluid pressure p over the range of z at selected time points t demonstrated agreement between the two solvers (Fig. 2).

Fig. 2
1D biphasic confined compression analysis, comparing the standard biphasic and BFSI solvers. The graphs show the (a) normalized displacement and (b) fluid pressure at selected time points with respect to normalized z-position. The pressure at t=1000 s is not shown here as it is zero for all z/h.
Fig. 2
1D biphasic confined compression analysis, comparing the standard biphasic and BFSI solvers. The graphs show the (a) normalized displacement and (b) fluid pressure at selected time points with respect to normalized z-position. The pressure at t=1000 s is not shown here as it is zero for all z/h.
Close modal

3.2 Unsteady Flow Past Cylinder.

This problem analyzed 2D transient fluid flow past a cylinder at a Reynolds number of 100 and compared results between the CFD, FSI, and BFSI solvers in febio. The same geometry, material parameters, time increments, and boundary conditions were taken from the unsteady flow past cylinder example presented in the supplementary materials from our previous publication on febio's CFD solver [46], except for the flow disturbance at the inlet used to seed vortex shedding, which in this case was a prescribed wy of a single sawtooth profile of amplitude 0.01m/s and total duration of 2s. Here, we only looked at the case when the generalizedα spectral radius parameter was set to ρ=0. For the FSI and BFSI analyses, to replicate a CFD material, the domain was fixed, and the solid matrix was modeled as a compressible neo-Hookean material with Young's modulus E=109Pa and Poisson's ratio ν = 0. For the BFSI analysis, the fluid domain was additionally assigned ϕs=0 and k1=0m4/Ns. In all cases, the models had 27,100 nodes, 13,300 elements, and 80,394 DOFs, where the mesh was biased near the cylinder. The flow in this analysis evolved to a periodic state exhibiting the characteristic Karman vortex street where vortices were shed behind the cylinder. Results are shown for the drag and lift coefficients, CD and CL, respectively, in Fig. 3, where the results were the same between the FSI and BFSI solvers, which in turn were nearly identical to the results from the CFD solver. In all cases, the Strouhal number was 0.171.

Fig. 3
(a) Drag coefficient CD and (b) lift coefficient CL for unsteady flow past a cylinder at Re = 100, using the CFD, FSI, and BFSI solvers with ρ∞=0
Fig. 3
(a) Drag coefficient CD and (b) lift coefficient CL for unsteady flow past a cylinder at Re = 100, using the CFD, FSI, and BFSI solvers with ρ∞=0
Close modal

3.3 Bubble Inflation Problem.

This bubble inflation problem was described in Fig. 13 of Bathe and Ledezma [59] and in Sec. 3.3 of Shim et al. [47], which describes the corresponding febio FSI solution. Here, we compared the results from those prior analyses against the BFSI solver. In this problem, two fluid domains, Ωtopf and Ωbotf, were separated by a thin Mooney–Rivlin elastic solid domain, Ωs, discretized with solid (hexahedral) elements (Fig. 4). The left, right, and top surfaces of Ωtopf were exposed to ambient conditions (zero pressure), whereas the bottom surface of Ωbotf was subjected to a fluid pressure that increased linearly with time. No-slip conditions were imposed on the bottom wall of Ωtopf, the side walls of Ωbotf, and the fluid–solid interfaces Γtopfs and Γbotfs on both sides of Ωs. The model was run with time-step Δt=0.001s until t=0.25s. For the BFSI analysis, the fluid domains were modeled as BFSI materials with ϕs=0 and k1=0m4/Ns; the solid matrix was modeled as a compressible neo-Hookean material with Young's modulus E=109Pa and Poisson's ratio ν = 0. The febio model had 24,348 nodes, 11,882 elements, and 118,022 DOFs. The mesh deformation at the final time point is shown in Fig. 4 and may be compared to Fig. 13 of Bathe and Ledezma [59]. The time evolution of the maximum principal stretch and vertical displacement at the membrane center were compared between the FSI and BFSI solvers in febio and Fig. 14 of Bathe and Ledezma [59], showing agreement (Fig. 5).

Fig. 4
Bubble inflation problem, showing mesh in reference configuration (left) and final configuration at time t = 0.25 (right). The pressure was prescribed as a linearly increasing function of time p0(t) at the bottom face of Ωbotf.
Fig. 4
Bubble inflation problem, showing mesh in reference configuration (left) and final configuration at time t = 0.25 (right). The pressure was prescribed as a linearly increasing function of time p0(t) at the bottom face of Ωbotf.
Close modal
Fig. 5
Bubble inflation problem, comparing febio FSI [47] and BFSI results to those of Bathe and Ledezma [59] for (a) the vertical displacement and (b) the first principal stretch of the top of the bubble
Fig. 5
Bubble inflation problem, comparing febio FSI [47] and BFSI results to those of Bathe and Ledezma [59] for (a) the vertical displacement and (b) the first principal stretch of the top of the bubble
Close modal

3.4 Taylor Slip Problem.

This problem compared the analytical solution of the classical Taylor slip problem modeling Couette flow over a rigid porous–permeable layer [60] against the febio BFSI solver. This benchmark problem has been used previously in various studies to verify related finite element solvers [34,35,37]. Here, a fluid channel, Ωtopf, was located directly above a biphasic layer, Ωbotb, whose solid displacement was fixed in order to represent a rigid porous material. These domains shared an interface Γbf. A constant shearing velocity v0 was prescribed at the top of the fluid layer, while the fluid velocity was fixed at the bottom of the porous layer (Fig. 6). The steady-state fluid flux depended on the nondimensional parameters
δ=1h2μkϕf,ξ=h1h2,η=ϕfμfμ
(3.1)

where μf and μ are the true fluid viscosities in the fluid and biphasic layers, respectively, and h1 and h2 are the heights of the fluid and porous layers, respectively. Four cases were examined using different values for δ and η, with ξ=0.25 and ϕf=0.5. To replicate the analytical solution, inertia was neglected by setting all mass densities to zero, and the fluid was made nearly incompressible with K=1012Pa. The febio models had 2058 nodes, 960 elements, and 5664 DOFs, with the mesh biased at the bottom of the fluid layer and the top of the biphasic layer. In the results presented here, the relative fluid flux wx was normalized by v0. All febio solutions agreed with corresponding analytical solutions (Fig. 7).

Fig. 6
Couette flow of a viscous fluid Ωtopf over a rigid porous–permeable layer Ωbotb. All displacements were fixed, where wx=v0 at the top of the fluid layer and w=0 at the bottom of the biphasic layer.
Fig. 6
Couette flow of a viscous fluid Ωtopf over a rigid porous–permeable layer Ωbotb. All displacements were fixed, where wx=v0 at the top of the fluid layer and w=0 at the bottom of the biphasic layer.
Close modal
Fig. 7
Four different cases of the Taylor slip problem modeling Couette flow over a biphasic layer comparing febio results and analytical solutions for the normalized fluid flux versus normalized height
Fig. 7
Four different cases of the Taylor slip problem modeling Couette flow over a biphasic layer comparing febio results and analytical solutions for the normalized fluid flux versus normalized height
Close modal

3.5 Mixed Couette and Poiseuille Flow Over a Deformable Biphasic Layer.

This problem generalized the Taylor solution from Sec. 3.4 by also including a Poiseuille flow contribution and allowing the biphasic material to deform [58]. The boundary conditions were similar to Fig. 6, except that the displacement was only fixed at the bottom of the biphasic layer, and an additional pressure gradient dp/dx=10Pa/m was prescribed across the domains. For this solution, the nondimensional parameters δ and ξ retain the same form as in Eq. (3.1), whereas η is redefined as shown below; an additional nondimensional parameter α was needed to account for the solid matrix stiffness
η=μfϕfμα=ϕfμv0μsh2
(3.2)

where μs is the shear modulus of the porous solid in the biphasic layer. We examined cases where δ2=0.01,δ2=0.1, and δ2=1, while ξ=0.1, η = 1, and α=0.01. We also varied η=0.1, η = 1, and η = 10, while ξ=0.1, δ2=0.01, and α=0.01. We set ϕf=0.75 in all cases. The mesh was modified to also be biased across the interface Γbf. In all models, there were 6724 nodes, 3240 elements, and 20,008 DOFs. febio results for the normalized relative fluid flux, wx/v0, and the normalized solid displacement, ux/h2, were compared against the analytical solution, showing agreement (Fig. 8).

Fig. 8
Comparison of febio results and analytical solutions for mixed Couette and Poiseuille flow of a viscous fluid flowing tangentially over a deformable biphasic layer. The normalized fluid fluxes and displacements are compared while examining the influence of varying (a) δ2 and (b) η.
Fig. 8
Comparison of febio results and analytical solutions for mixed Couette and Poiseuille flow of a viscous fluid flowing tangentially over a deformable biphasic layer. The normalized fluid fluxes and displacements are compared while examining the influence of varying (a) δ2 and (b) η.
Close modal

3.6 Ultrafiltration of Viscous Fluid Through Biphasic Layer.

One-dimensional ultrafiltration of a viscous fluid through a biphasic layer was analyzed in febio, and the results were compared against novel analytical solutions derived in this study and presented in the Appendix. The configuration for this problem is shown in Fig. 9, where a rigid impermeable platen Ωtops meets a viscous fluid layer Ωtopf at the interface Γtopfs located at z=h1, where Ωtopf then interfaces a biphasic layer Ωbotb at Γbf located at z =0. The bottom of the biphasic layer is in contact with a stationary free-draining rigid porous platen Ωbots at the interface Γbotbs located at z=h2. The rigid platen Ωtops was prescribed either a step velocity, v0, or step pressure, p0, applied at t=0+.

Fig. 9
Boundary conditions for the 1D ultrafiltration of a viscous fluid through a biphasic layer problem. Ωtops is a rigid body domain, Ωtopf is a deformable fluid domain, Ωbotb is a BFSI domain, and Ωbots is a rigid porous platen.
Fig. 9
Boundary conditions for the 1D ultrafiltration of a viscous fluid through a biphasic layer problem. Ωtops is a rigid body domain, Ωtopf is a deformable fluid domain, Ωbotb is a BFSI domain, and Ωbots is a rigid porous platen.
Close modal

The problem was modeled in febio using K=1014Pa to nearly enforce fluid incompressibility, biphasic porosity ϕf=0.7, upstream true fluid viscosity μf=1Pa·s, biphasic true fluid viscosity μ=1Pa·s, biphasic solid Young's modulus E=9×105Pa, Poisson's ratio ν=0.2 (producing HA=106Pa), and hydraulic permeability k=1014m4/Ns. Either a step pressure p0=104Pa or a step velocity v0=108m/s was prescribed. The mesh had 1276 nodes, 570 elements, and 3675 DOFs, biased at all interfaces. The model ran up to t=1000s, where the time-step size started at Δt=0.001s and gradually increased to Δt=1s such that it increased by an order of magnitude approximately every 500 time-steps. The resulting normalized displacement, normalized fluid velocity, and pressure were compared with the analytical solutions at t=1s, t=10s,t=100s, and t=1000s, as reported in Fig. 10, showing good agreement. Note that in the step pressure case, the value of v0 used in calculating the normalized fluid velocity was obtained from Eq. (A15) as t.

Fig. 10
Comparison between analytical solutions and febio results for 1D ultrafiltration of viscous fluid through a biphasic layer, showing solid displacement uz, fluid velocity vzf, and fluid pressure p when prescribing (a) a step velocity and (b) a step pressure on the upstream fluid. The t=1000 s case is not shown for the step pressure graphs because the results were the same as the t=100 s cases.
Fig. 10
Comparison between analytical solutions and febio results for 1D ultrafiltration of viscous fluid through a biphasic layer, showing solid displacement uz, fluid velocity vzf, and fluid pressure p when prescribing (a) a step velocity and (b) a step pressure on the upstream fluid. The t=1000 s case is not shown for the step pressure graphs because the results were the same as the t=100 s cases.
Close modal

3.7 Dilatational Wave Propagation.

In our recent formulation of the hybrid biphasic theory [44], the governing equations for one-dimensional isothermal dilatational wave propagation were presented for the special case of a linear isotropic elastic solid, inviscid fluid, infinitesimal solid and fluid strains, zero body forces, and negligible frictional drag between the fluid and solid constituents, while inertia terms were included. These governing equations led to the identification of two wave propagation speeds, cf in the fluid and cs in the solid
cf=KρTf,cs=HAρs
(3.3)

where HA is defined in the Appendix. Here, we verified the existence of these two wave propagations, as well as the accuracy of the formulas for cf and cs in febio, using the BFSI solver.

We used a configuration similar to the dilatational wave propagation problem examined in our previous implementation of the CFD solver [46]. The rectangular hybrid biphasic domain has length 1m along the x-direction (0x1), height 0.1m along y, and depth 0.05m along z. The prescribed properties were ϕs=0.2, zero true fluid viscosity, ρTrf=1.225kg/m3, and K=101,325Pa. The porous solid was a compressible neo-Hookean hyperelastic material with ρTs=103kg/m3 and E=0.16GPa. Here, k1 was set to zero to eliminate frictional drag between the fluid and solid constituents. At x=1m, solid displacement uxwas fixed, while natural boundary conditions were enforced for the fluid (wx = 0). At x =0, we prescribed either a solid displacement ux=(1/2)u0(1cos2πt/T) or a normal fluid flux wx=(1/2)w0(1cos2πt/T), with u0=106m, w0=103m/s, and T=3×104s. In the case when wx was prescribed, ux was subsequently fixed at x =0 for t > T, to allow wave reflection on this boundary. The solid displacements along y and z were explicitly fixed, while the normal component of the fluid flux was naturally zero on the remaining boundaries.

Each analysis was run for 1200 uniform time-steps to a final time of 0.006s. The generalized-α parameter ρ was set to 1 in order to verify numerical energy conservation in this materially nondissipative problem. The mesh was composed of 2000 uniform hexahedral elements in the x-direction, resulting in 8004 nodes and 40,016 DOFs. febio results for both test cases are shown in Fig. 11 at four selected time points. When prescribing a solid dilatational disturbance, disturbances traveled in the fluid at two distinct speeds cs and cf, as manifested in the fluid pressure p(x,t). Similarly, when prescribing a fluid dilatational disturbance, disturbances traveled in the solid at the same distinct speeds, as manifested in the solid velocity u˙x(x,t). Disturbances reflected at both ends of the domain, with u˙x(x,t) inverting in sign for each reflection. In both analyses, the wave velocities calculated from the febio results were found to match the theoretical values of cf=287.6m/s and cs=894.4m/s evaluated from Eq. (3.3), throughout the duration of the analysis. No measurable numerical dissipation occurred in the febio solution, indicating that energy was properly conserved.

Fig. 11
One-dimensional dilatational wave propagation through hybrid biphasic material under prescribed dilatational disturbances in (a) the solid and (b) the fluid. Plots represent the fluid pressure p(x,t) in (a) and the solid velocity u˙x(x,t) in (b), at four selected time points t. Both analyses produced a superposition of two dilatational wave propagations, one traveling at a slower speed cf (line arrow) and the other at the faster speed cs (solid arrow).
Fig. 11
One-dimensional dilatational wave propagation through hybrid biphasic material under prescribed dilatational disturbances in (a) the solid and (b) the fluid. Plots represent the fluid pressure p(x,t) in (a) and the solid velocity u˙x(x,t) in (b), at four selected time points t. Both analyses produced a superposition of two dilatational wave propagations, one traveling at a slower speed cf (line arrow) and the other at the faster speed cs (solid arrow).
Close modal

3.8 Two-Dimensional Airflow Through a Simulated Face Mask.

In this example, we modeled 2D airflow through a simulated face mask, modeled as a hybrid biphasic material, at different values of its hydraulic permeability, to examine how much a face mask might mitigate dispersion of expiratory flow. Gravitational body force was not included in this analysis. Two fluid domains (air) were separated by a thin deformable BFSI domain (face mask) to simulate inspiratory and expiratory flow from the lungs into the ambient environment, through the nose and mouth. The boundary conditions and simplified geometry are described in Fig. 12. Briefly, the simplified lungs–mouth–nose domain, Ωlmnf, had a height h2=0.1m and length l=0.2m, connecting with the thin biphasic face mask domain, Ωmaskb, which had a thickness of t=0.002m. This mask domain interfaced with a larger ambient air domain, Ωambf, which had height H=1m and length L=2m. All domains had a uniform depth of 0.01m. A plug velocity profile vn was prescribed on the leftmost boundary of Ωlmnf to replicate cyclical breathing conditions. The exact temporal profile of vn was calculated to replicate the flow produced by a clinical mechanical ventilator (Getinge, Servo-u Ventilator). To produce 24 breaths per minute, the period of each cycle was set to T=2.5s. To produce a physiological inspiratory:expiratory ratio, the first 0.9s consisted of inspiration and the remainder of the cycle was expiration. A plot of the prescribed inlet velocity vn calculated for this geometry is presented for one cycle in Fig. 13. The febio model ran for a total of five cycles, with time-step Δt=0.00025s. Zero fluid pressure was prescribed on the boundaries of Ωambf as shown in Fig. 12. For the standard face mask, the hydraulic permeability k was set to 3×106m4/(Ns), calculated to produce a pressure drop of approximately 600Pa across the mask, which is a value typically found in a mechanical ventilator. This value of k was also compared to k=3×104m4/(Ns) and k1=0 to investigate the effects of a more porous mask and no mask, respectively. The air was modeled as Newtonian viscous with ρTrf=1.2047kg/m3,K=101,325Pa, and μ=1.8205×105Pa·s. For nonzero k1, the properties of the solid constituent of the mask were approximated as ρTs=80kg/m3, E=108Pa, and ν=0.4, with ϕrs=0.75. Tangential and backflow stabilization boundary conditions were prescribed on the zero fluid pressure boundaries of Ωambf. For all models, there were 109,066 nodes, 53,982 elements, and 540,284 DOFs, in a mesh that was biased across interfaces Γbf and normal to no-slip walls to capture any sharp gradients in the fluid flow or pressure.

Fig. 12
Boundary conditions for 2D airflow through a simulated face mask. Here, Ωlmnf is the fluid domain representing the lungs, mouth, and nose, Ωmaskb is the BFSI domain representing the face mask, and Ωambf represents the ambient air.
Fig. 12
Boundary conditions for 2D airflow through a simulated face mask. Here, Ωlmnf is the fluid domain representing the lungs, mouth, and nose, Ωmaskb is the BFSI domain representing the face mask, and Ωambf represents the ambient air.
Close modal
Fig. 13
Fluid velocity vn prescribed at leftmost boundary of Ωlmnf, for one respiratory cycle lasting 2.5 s. Inspiration (positive normal velocity vn, implying leftward velocity) occurs during the first 0.9 s, while expiration makes up the remainder of the cycle.
Fig. 13
Fluid velocity vn prescribed at leftmost boundary of Ωlmnf, for one respiratory cycle lasting 2.5 s. Inspiration (positive normal velocity vn, implying leftward velocity) occurs during the first 0.9 s, while expiration makes up the remainder of the cycle.
Close modal

The vorticity in the z-direction and the magnitude of w are displayed at two selected time points in Fig. 14. Time points t=8s and t=11.5s occur during inspiration of the fourth breath and expiration of the fifth breath, respectively. For the standard face mask, the fluid pressure ranged from 724Pa to 791Pa over the entire domain and duration of the analysis. For the face mask with higher k, the corresponding range was from 13.5Pa to 10.4Pa, and in the no-mask case it was 7.95Pa to 5.49Pa. These simulation results suggest that a face mask does indeed reduce dispersion of air during inspiration and expiration, for both values of k employed in this study, when compared to the no-mask case, though the effect is not particularly dramatic. Effectively, the face mask (which acts as a porous filter) produces more laminar flow during expiration. The mask also slows down the expiratory flow, which travels a slightly shorter distance at a given time point with decreasing mask permeability. However, as simulated here, the mask does not appear to reduce the range of air travel produced by respiratory expiration, at least over the 2 m distance examined here, in a 2D analysis.

Fig. 14
Comparison of 2D airflow through a simulated face mask, using different hydraulic permeabilities k, at two selected time points, during inspiration of the fourth breath (t=8 s) and expiration of the fifth breath (t=11.5 s). (a) Fluid z-vorticity and (b) magnitude of relative fluid velocity are shown.
Fig. 14
Comparison of 2D airflow through a simulated face mask, using different hydraulic permeabilities k, at two selected time points, during inspiration of the fourth breath (t=8 s) and expiration of the fifth breath (t=11.5 s). (a) Fluid z-vorticity and (b) magnitude of relative fluid velocity are shown.
Close modal

From a numerical perspective, it is interesting to note that small flow disturbances appear whenever shedded vortices cross the zero fluid pressure boundaries of Ωambb, as is apparent at the right boundary at t=8s, and the bottom and right boundaries at time t=11.5s, in Fig. 14. These disturbances are likely artifacts produced by the backflow and tangential stabilizations prescribed on those boundaries. However, removing those boundary conditions leads to premature termination of the computational analysis.

3.9 Carotid Bifurcation With Porous Walls.

This problem analyzed a 3D model of the bifurcated carotid artery, with thick and porous arterial walls modeled as BFSI materials. Similar to the bifurcated carotid artery example from our earlier study [47], the model was obtained from the GrabCAD community.4 It was converted into units of meters and meshed to include a boundary layer refinement using four-node tetrahedral elements, defining the blood domain, Ωbldf, as an FSI domain to account for its deformation. The inlet, Γinf, had a diameter of 6.28mm, and the two outlets, Γlgf and Γsmf, had respective diameters of 4.24mm and 3.04mm. The arterial wall was created by extruding the outer surface of the artery twice, with the tunica media given a thickness of 0.3mm and the tunica adventitia a thickness of 0.15mm, for a total arterial wall thickness of 0.45mm. These domains correspond to Ωmedb and Ωadvb. The intima layer was not modeled explicitly, as it was assumed to provide a negligible contribution to the mechanical behavior of the arterial wall. All elements in the arterial wall consisted of six-node pentahedra.

The solid constituent of the arterial wall was assigned the constitutive model by Holzapfel et al. [61], which describes the tissue as a fiber reinforced composite with two fiber families. The material properties used here were obtained from Sommer and Holzapfel [62], using the reported mean values of the human common carotid at the lower pressure range: c=122,300Pa, k1=24,700Pa, and k2=16.5 in the media, and c=59,600Pa, k1=180,900Pa, and k2=109.8 in the adventitia. The fiber families formed an angle of ±6.9deg in the media and ±30.1deg in the adventitia relative to the circumferential direction, which was determined at each element from the local direction of maximum principal curvature of the outer surface of each layer. In addition, we set ρTs=1000kg/m3 for the dynamic response of the solid constituents of the domains Ωmedb and Ωadvb, with solid bulk modulus 108Pa to enforce a nearly incompressible response. Both blood and the interstitial fluid within the arterial wall were modeled as Newtonian fluids where ρTrf=1000kg/m3, K=2.2×109Pa, and μ=0.004Pa·s.

For the boundary conditions, an inlet velocity w=wnn was prescribed on Γinf having a parabolic spatial profile, and an average value w0(t) that rises from a diastolic plateau of 0.10m/s to a systolic maximum of 0.48m/s in three consecutive cycles such that Re=165800, where the graph of the time history of w0(t) can be found in our previous publication [47]. A flow resistance boundary condition [46,63] was prescribed on both outlet faces Γlgf and Γsmf, with R=4×108kg/m4s and a pressure offset of p=104Pa. The inlet was initially ramped up smoothly from zero to the diastolic value in 0.05s, while the outlet pressure offset was ramped up smoothly over 0.25s, to minimize the effects of pulse waves propagating in the opposite directions from the inlet and outlets. Backflow and tangential stabilization were prescribed on both outlet faces Γlgf and Γsmf. A fluid dilatation corresponding to a pressure of p=104Pa was prescribed at the inlet rim between the media and blood domains to enhance the stability of the numerical solution and was ramped up smoothly in 0.25s. The outer surface of the adventitia, Γwallb, was prescribed zero fluid pressure. The artery was constrained by fixing the displacements of all nodes located on the inlet and outlet faces, including that of the blood and wall. Given these model configurations and boundary conditions, the model had 60,223 nodes, 225,467 elements, and 414,691 DOFs. The models were set to run 1550 time-steps with Δt=1ms.

In this problem, we compared the effects of applying different hydraulic permeability values for the arterial wall layers, where in one case we used a realistic value of k=5.26×1016m4/(Ns) [31], while in the second case we set k=5.26×1011m4/(Ns) to exaggerate the permeability of the arterial wall. In both cases, the referential solid volume fraction was set to ϕrs=0.5. The geometry of the artery as well as a vector plot of the fluid flux through the arterial wall in the exaggerated permeability case during systole (t=0.87s) is shown in Fig. 15. The fluid flux from the blood domain into the arterial wall was relatively uniform along the length of the artery and was higher during systole as expected, although for both cases the magnitudes were small compared to the overall blood velocity. In addition, the fluid pressure, shear stress, and fluid flux results were very similar between the different permeability values. We also investigated the conservation of mass in this problem, as reported in Fig. 16, by examining the mass flowrate, Γρfw·nda, across each of the boundaries of the fluid domain Ωbldf, including the arterial wall, Γwallb, along with the time rate of change of total fluid mass m=Ωρfdv in the overall arterial lumen Ωbldf. The sum of all mass flowrates, which should be zero due to the axiom of mass balance, produced an error less than 1.51% relative to the peak inlet mass flowrate across the three consecutive flow cycles, for both permeability values, verifying that mass was satisfactorily conserved. By design, the mass flowrate across the arterial wall had a significant contribution in the exaggerated permeability case.

Fig. 15
(a) Plane-cut geometry of bifurcated carotid artery defining each of the domains and surfaces. (b) Fluid flux vector plots are shown through the arterial wall during peak velocity at t=0.87 s.
Fig. 15
(a) Plane-cut geometry of bifurcated carotid artery defining each of the domains and surfaces. (b) Fluid flux vector plots are shown through the arterial wall during peak velocity at t=0.87 s.
Close modal
Fig. 16
Comparison of the mass flowrates across each of the boundaries of Ωbldf (inlet, outlet 1, outlet 2, and wall), the time rate of change m˙ of fluid mass in Ωbldf, and the sum of these measures (total m˙), in the case of (a) exaggerated and (b) realistic hydraulic permeability values
Fig. 16
Comparison of the mass flowrates across each of the boundaries of Ωbldf (inlet, outlet 1, outlet 2, and wall), the time rate of change m˙ of fluid mass in Ωbldf, and the sum of these measures (total m˙), in the case of (a) exaggerated and (b) realistic hydraulic permeability values
Close modal

4 Discussion

The objective of this study was to implement a novel BFSI solver in febio that enhances existing modeling capabilities for biological tissues and other fluid–solid mixtures. febio is an open-source finite element software specialized for biomechanics and biophysics applications that features a variety of solvers for deformable and rigid solid materials, as well as fluid, FSI, biphasic, and multiphasic materials. Our recently developed FSI solver [47] made it possible to model interactions between viscous fluids and surrounding deformable or rigid solids under dynamic conditions. Just like our earlier CFD solver [46], this FSI solver modeled the fluid as slightly compressible and used the fluid dilatation as a nodal DOF to simplify the finite element formulation, as explained in those prior publications. However, this FSI solver is not compatible with the biphasic and multiphasic solvers currently available in febio, because the standard febio biphasic material consists of intrinsically incompressible fluid and solid constituents, with no explicit modeling of fluid viscosity, undergoing quasi-static motion. Its nodal fluid pressure represents a Lagrange multiplier that enforces the intrinsic incompressibility constraint of the fluid and solid [5,11], producing an incompatibility of nodal DOFs when interfacing FSI and biphasic domains.

To overcome this incompatibility, we first formulated a novel theory for a hybrid biphasic material under isothermal conditions, whose interstitial fluid is viscous and compressible, without neglecting inertia and body forces [44]. This hybrid formulation was carefully crafted to reduce to the standard biphasic theory in the limit when the fluid is intrinsically incompressible (Jf=1), as described in Ref. [44]. This study describes the finite element modeling of this hybrid biphasic material and its integration with the existing CFD/FSI solvers in febio, to overcome modeling limitations and enhance users' ability to analyze a wider range of problems in biomechanics where dynamic interactions of a fluid with a porous deformable tissue play an important role.

The BFSI solver was verified using a variety of problems that examined special cases and compared their solutions with those of the CFD, FSI, and standard biphasic solvers. In particular, the BFSI solver was verified against a standard biphasic 1D confined compression problem (Fig. 2) in the limiting case of an inviscid interstitial fluid whose bulk modulus was set to a sufficiently large value, while also neglecting inertia. In this type of application, it should be noted that the BFSI solver runs slower than the standard biphasic solver, since the BFSI domain has more degrees-of-freedom (u, w, and ef, instead of just u and p). The unsteady flow past a cylinder problem verified that the BFSI solver can replicate viscous fluid simulations when compared to the CFD and FSI solvers in the limit when the solid displacement is fixed (Fig. 3). However, while the FSI and BFSI solvers produced identical results, the results between the BFSI and CFD solvers were very close but not exactly the same. This difference occurred because the generalized-α scheme was implemented for a second-order system in the FSI/BFSI solvers compared to the CFD solver, where it was implemented as a first-order system. Similarly, the bubble inflation problem originally described by Bathe and Ledezma [59] could be reproduced with the BFSI solver in the limit when its solid matrix was ascribed negligible stiffness, zero mass, and zero frictional interactions with the interstitial fluid, verifying the solver's ability to analyze very large deformations (Fig. 5). Additional BFSI verification problems which could not be simulated with existing solvers employed analytical solutions that were either available in the prior literature or derived for the purposes of this study. These include the Taylor slip problem (Fig. 7) and mixed Couette and Poiseuille flow over a biphasic layer (Fig. 8), which examined steady-state conditions, ultrafiltration of a viscous fluid through a biphasic layer which examined quasi-static transient conditions (Fig. 10), and dilatational wave propagation which examined dynamic effects under infinitesimal strains and nondissipative conditions (Fig. 11). Then, the BFSI solver was used to examine a 2D simulation of airflow through a face mask (Fig. 14), illustrating an application of this solver to a problem with no prior known results. Finally, the solver investigated 3D blood flow though a bifurcated carotid artery with a porous–permeable arterial wall (Fig. 15), whose solid constituents consisting of a heterogenous, anisotropic two-fiber family composite. This example showed that mass was satisfactorily conserved under physiologic and exaggerated hydraulic permeabilities (Fig. 16).

The finite element implementation of BFSI used the mixture theory framework, where material time derivatives were evaluated in the material frame for the solid, since the mesh is defined on that constituent, and in the spatial frame for the fluid, since it flows through the porous solid. The mixture mass balance and the solid and fluid momentum balance equations were used as governing equations and solved monolithically. The DOFs were the solid displacement u, relative fluid velocity w defined in Eq. (2.4), and fluid dilatation ef. The momentum equations were rewritten into the form of Eq. (2.18), such that the finite element formulation automatically satisfied all the jump conditions presented in Sec. 2.2 across interfaces Γbb and Γbf.

As a generalization of the CFD, FSI, and biphasic solvers in febio, the BFSI solver has many similarities when it comes to user experience and accessibility. Many of the same operating recommendations from these existing solvers apply here, such as the need for a biased mesh in the presence of large flow and pressure gradients. In particular, models with lower hydraulic permeability values require more refined meshes and take longer to complete, as observed in the 3D bifurcated carotid artery blood flow problem where the more permeable case finished in 10h while the realistic case took 24h. In addition, in the stiffness matrix, there are terms that reduce to zero when using finite elements with linear DOF interpolations, such as those that include gradJ. While the presence of such terms might suggest that quadratic elements may be more accurate or produce better convergence, to date we have not noticed a difference in solution quality when using linear versus quadratic elements. Furthermore, any putative advantages may be offset by the added computational costs of using quadratic meshes.

Since febio is an open-source software, users can take advantage of the solver simplicity to create their own plugins or extensions of this framework for their specific applications. All constitutive material models and boundary conditions are described in their own C++ classes, which are all children of the same parent classes for materials and boundary conditions. Virtual functions are defined that give each of these classes a basic template. For instance, all material model classes have functions for evaluating the stress and its tangent with respect to the relevant DOFs. This modular structure allows the user to easily create and switch between different boundary conditions or materials, such as between a Newtonian and non-Newtonian fluid. As a result, the BFSI implementation in febio is uniquely positioned to meet users' modeling needs for biphasic mixtures, with both functionality and customizability.

Validation of this BFSI solver against experimental data was difficult to achieve in this study, due to limited availability of benchmark problems. However, with the introduction of this BFSI solver, biological models that were previously impossible to analyze are now possible, most notably those involving fluid exchanges between dynamically deforming fluid and biphasic domains or wave propagation analyses in hydrated biological tissues. The BFSI solver can be utilized in a variety of problems such as cartilage lubrication or blood flow through and across porous arteries and cardiac wall, among many others. These types of studies may capitalize on the modular structure of the code and extensive library of solid and biphasic materials available in febio, which include anisotropic fiber distributions [64], viscoelasticity [65], active contraction, and damage mechanics [66], as well as the library of non-Newtonian viscous fluids [46,67].

The BFSI solver formulated in this study represents an important step in developing robust solvers that can be used for a wide range of biomechanics and biophysics applications. Various improvements of the current BFSI solver may be developed to improve its functionality. For example, tied contact interfaces may be formulated to accommodate dissimilar meshes at interfaces Γbb,Γbf, or Γbs. With the completion of this formulation, the natural next step is to implement solute transport and reactive processes among all constituents, expanding on similar features currently available in the febio multiphasic solver [51,68]. This type of multiphasic-FSI solver could accommodate charged solutes and solid matrix, osmotic effects, and reactions involving solutes and solid-bound molecules, such as those involved in growth and remodeling. This extended solver would potentially be capable of analyzing more elaborate processes, such as aerosol and droplet transport in airflow through a porous barrier such as face masks or ventilation filters, while also accounting for droplet evaporation, or binding of aerosols to filters and solid surfaces, as would be relevant to the study of SARS-Cov-2 transmission.

Acknowledgment

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the National Science Foundation.

Funding Data

  • Division of Graduate Education, U.S. National Science Foundation (Grant No. NSF GRFP DGE-16-44869; Funder ID: 10.13039/100000001).

  • National Institute of General Medical Sciences, U.S. National Institutes of Health (Award No. R01GM083925; Funder ID: 10.13039/100000002).

Footnotes

Appendix: Ultrafiltration of Viscous Fluid Through Biphasic Layer

In order to solve for this problem analytically, we assumed no body forces, linear, homogenous, and isotropic material models for both fluid and solid constituents, infinitesimal deformation, no inertial effects, and intrinsic incompressibility of solid and fluid constituents. Using governing equations from Sec. 2.1, applying these assumptions and simplifying to 1D, the normal stress components reduce to
σzzs=ϕsp+HAϵzzσzzf=ϕfp+43μaDzz
(A1)
where HA=2μs+λs is the aggregate modulus of the porous solid matrix, ϵzz=uz/z is the normal component of the infinitesimal strain tensor, Dzz=vzf/z is the normal component of the rate of deformation tensor, and μa is the fluid shear viscosity. The governing equations now reduce to
0=z(uzt+wz)
(A2)
0=ϕspz+HA2uzz2+ϕf1kwz
(A3)
0=ϕfpz+43μa2vzfz2ϕf1kwz
(A4)
where vzf=(1/ϕf)wz+(uz/t) is obtained from Eq. (2.4). Based on our assumptions, all volume fractions can be considered constant. These equations can also be obtained from Hou et al. [34,41]. From Eqs. (A2)(A4), in the fluid domain Ωtopf
0=vzfz0=pz
(A5)
Therefore, in Ωtopf, we find that vzf(t) and p(t) are uniform in space, satisfying vzf(t)=v0(t) and p(t)=p0(t). When a step pressure is applied, p(t)=p0H(t) where H(t) is the Heaviside unit step function, we need to solve for vzf(t)=v0(t), and vice versa. Using the jump conditions from Sec. 2.2 reduced to the case of an incompressible fluid (also given in Ref. [41]), the boundary conditions for the biphasic layer at z =0 are
v0(t)=ϕsuzt|z=0+ϕfvzf(z=0,t)p0(t)=p(z=0,t)0=43μavzfz|z=0+HAuzx|z=0
(A6)
and at z=h2 they are
0=uz(z=h2,t)0=p(z=h2,t)v0(t)=ϕfvzf(z=h2,t)
(A7)
The initial conditions are
0=uz(z,t=0)0=vzf(z,t=0)0=p(z,t=0)
(A8)
The analytical solution of these equations was obtained using the method of Laplace transforms for uz(z,t),vzf(z,t), and p(z,t) within the biphasic layer, for both loading cases. First, we report the analytical solution within the biphasic layer under steady-state conditions
uz(z)=12HAkv0(z2h22)p(z)=p0(11h2z)vzf=1ϕfv0
(A9)
Then, the transient solution in the case of a step velocity within the biphasic layer is given by
uz(z,t)=v0HAk[12(z2h22)+2h22n=1(1)ncos(γnh2z)γn3eγn2τ(1+43(ϕs)2δ2γn2)t]
(A10)
vzf(z,t)=1ϕfv0H(t)+2ϕsϕfv0n=111+43(ϕs)2δ2γn2(1)ncos(γnh2z)γneγn2τ(1+43(ϕs)2δ2γn2)t
(A11)
p(z,t)=83ϕsμv0h2n=1(1)nsin(γnh2z)1+43(ϕs)2δ2γn2eγn2τ(1+43(ϕs)2δ2γn2)t+v0k[z+2h2n=1(1)n+1sin(γnh2z)γn2eγn2τ(1+43(ϕs)2δ2γn2)t]+p0(t)p0(t)=83ϕsμv0h2n=111+43(ϕs)2δ2γn2eγn2τ(1+43(ϕs)2δ2γn2)t+h2v0k[12n=11γn2eγn2τ(1+43(ϕs)2δ2γn2)t]
(A12)
where γn=π(n12)
δ2=1h22μkϕf
as in the Taylor slip problem from Sec. 3.4, and
τ=h22HAk
which is the characteristic time constant for the response of a standard biphasic material. The transient solution in the case of a step pressure within the biphasic layer is given by
uz(z,t)=p0HA(z2h222h2+2h2n=11+(1)n+1cos(γnh2z)γn2(1+43ϕsδ2γn2)eγn2τ(1+43(ϕs)2δ2γn2)t)
(A13)
p(z,t)=43ϕsμp0k(2h22n=1γn1+43(ϕs)2δ2γn2(1)n+1sin(γnh2z)1+43ϕsδ2γn2eγn2τ(1+43(ϕs)2δ2γn2)t)p0(zh2+2n=1(1)nsin(γnh2z)γn(1+43ϕsδ2γn2)eγn2τ(1+43(ϕs)2δ2γn2)t)+p0H(t)
(A14)
vzf(z,t)=2ϕsϕfp0kh2n=111+43(ϕs)2δ2γn21+(1)ncos(γnh2z)1+43ϕsδ2γn2eγn2τ(1+43(ϕs)2δ2γn2)t+1ϕfv0(t)v0(t)=83(ϕs)2δ2p0kh2n=1γn21+43(ϕs)2δ2γn2(1)n+1cos(γnh2z)1+43ϕsδ2γn2eγn2τ(1+43(ϕs)2δ2γn2)t+p0kh2(1+2n=1(1)ncos(γnh2z)1+43ϕsδ2γn2eγn2τ(1+43(ϕs)2δ2γn2)t)2p0kh2n=111+43(ϕs)2δ2γn21+(1)ncos(γnh2z)1+43ϕsδ2γn2eγn2τ(1+43(ϕs)2δ2γn2)t
(A15)
where, in this case, γn=πn. It is important to note that fluid viscosity, which enters into δ2, only affects the transient response, with increasing values of δ2 producing a longer characteristic time constant for the displacement, pressure, and fluid velocity responses. The inviscid response may be recovered by setting δ2=0.

References

1.
Truesdell
,
C.
, and
Toupin
,
R.
,
1960
, “
The Classical Field Theories
,”
Encyclopedia of Physics
, Vol.
III/1
,
Springer-Verlag
,
Berlin, Germany
.
2.
Green
,
A. E.
, and
Naghdi
,
P. M.
,
1969
, “
On Basic Equations for Mixtures
,”
Q. J. Mech. Appl. Math.
,
22
(
4
), pp.
427
438
.10.1093/qjmam/22.4.427
3.
Bowen
,
R. M.
,
1976
, “
Theory of Mixtures
,”
Continuum Physics
,
Academic Press
, New York, pp.
1
127
.
4.
Bedford
,
A.
, and
Drumheller
,
D. S.
,
1983
, “
Theories of Immiscible and Structured Mixtures
,”
Int. J. Eng. Sci.
,
21
(
8
), pp.
863
960
.10.1016/0020-7225(83)90071-X
5.
Mow
,
V. C.
,
Kuei
,
S. C.
,
Lai
,
W. M.
, and
Armstrong
,
C. G.
,
1980
, “
Biphasic Creep and Stress Relaxation of Articular Cartilage in Compression: Theory and Experiments
,”
ASME J. Biomech. Eng.
,
102
(
1
), pp.
73
84
.10.1115/1.3138202
6.
Oomens
,
C. W.
,
van Campen
,
D. H.
, and
Grootenboer
,
H. J.
,
1987
, “
A Mixture Approach to the Mechanics of Skin
,”
J. Biomech.
,
20
(
9
), pp.
877
885
.10.1016/0021-9290(87)90147-3
7.
Lai
,
W. M.
,
Hou
,
J. S.
, and
Mow
,
V. C.
,
1991
, “
A Triphasic Theory for the Swelling and Deformation Behaviors of Articular Cartilage
,”
ASME J. Biomech. Eng.
,
113
(
3
), pp.
245
258
.10.1115/1.2894880
8.
Huyghe
,
J. M.
, and
Janssen
,
J.
,
1997
, “
Quadriphasic Mechanics of Swelling Incompressible Porous Media
,”
Int. J. Eng. Sci.
,
35
(
8
), pp.
793
802
.10.1016/S0020-7225(96)00119-X
9.
Humphrey
,
J. D.
, and
Rajagopal
,
K. R.
,
2003
, “
A Constrained Mixture Model for Arterial Adaptations to a Sustained Step Change in Blood Flow
,”
Biomech. Model. Mechanobiol.
,
2
(
2
), pp.
109
126
.10.1007/s10237-003-0033-4
10.
Kenyon
,
D. E.
,
1976
, “
The Theory of an Incompressible Solid-Fluid Mixture
,”
Arch. Ration. Mech. Anal.
,
62
(
2
), pp.
131
147
.10.1007/BF00248468
11.
Bowen
,
R. M.
,
1980
, “
Incompressible Porous Media Models by Use of the Theory of Mixtures
,”
Int. J. Eng. Sci.
,
18
(
9
), pp.
1129
1148
.10.1016/0020-7225(80)90114-7
12.
Biot
,
M. A.
, and
Temple
,
G.
,
1972
, “
Theory of Finite Deformations of Porous Solids
,”
Indiana Univ. Math. J.
,
21
(
7
), pp.
597
620
.10.1512/iumj.1972.21.21048
13.
Ateshian
,
G. A.
,
2009
, “
The Role of Interstitial Fluid Pressurization in Articular Cartilage Lubrication
,”
J. Biomech.
,
42
(
9
), pp.
1163
1176
.10.1016/j.jbiomech.2009.04.040
14.
Ateshian
,
G. A.
,
Warden
,
W. H.
,
Kim
,
J. J.
,
Grelsamer
,
R. P.
, and
Mow
,
V. C.
,
1997
, “
Finite Deformation Biphasic Material Properties of Bovine Articular Cartilage From Confined Compression Experiments
,”
J. Biomech.
,
30
(
11–12
), pp.
1157
1164
.10.1016/S0021-9290(97)85606-0
15.
Huang
,
C.-Y.
,
Soltz
,
M. A.
,
Kopacz
,
M.
,
Mow
,
V. C.
, and
Ateshian
,
G. A.
,
2003
, “
Experimental Verification of the Roles of Intrinsic Matrix Viscoelasticity and Tension-Compression Nonlinearity in the Biphasic Response of Cartilage
,”
ASME J. Biomech. Eng.
,
125
(
1
), pp.
84
93
.10.1115/1.1531656
16.
Park
,
S.
,
Krishnan
,
R.
,
Nicoll
,
S. B.
, and
Ateshian
,
G. A.
,
2003
, “
Cartilage Interstitial Fluid Load Support in Unconfined Compression
,”
J. Biomech.
,
36
(
12
), pp.
1785
1796
.10.1016/S0021-9290(03)00231-8
17.
Smith
,
J. H.
, and
García
,
J. J.
,
2009
, “
A Nonlinear Biphasic Model of Flow-Controlled Infusion in Brain: Fluid Transport and Tissue Deformation Analyses
,”
J. Biomech.
,
42
(
13
), pp.
2017
2025
.10.1016/j.jbiomech.2009.06.014
18.
Lande
,
B.
, and
Mitzner
,
W.
,
2006
, “
Analysis of Lung Parenchyma as a Parametric Porous Medium
,”
J. Appl. Physiol.
,
101
(
3
), pp.
926
933
.10.1152/japplphysiol.01548.2005
19.
Ricken
,
T.
,
Dahmen
,
U.
, and
Dirsch
,
O.
,
2010
, “
A Biphasic Model for Sinusoidal Liver Perfusion Remodeling After Outflow Obstruction
,”
Biomech. Model. Mechanobiol.
,
9
(
4
), pp.
435
450
.10.1007/s10237-009-0186-x
20.
Chapelle
,
D.
,
Gerbeau
,
J.-F.
,
Sainte-Marie
,
J.
, and
Vignon-Clementel
,
I. E.
,
2010
, “
A Poroelastic Model Valid in Large Strains With Applications to Perfusion in Cardiac Modeling
,”
Comput. Mech.
,
46
(
1
), pp.
91
101
.10.1007/s00466-009-0452-x
21.
Cimrman
,
R.
, and
Rohan
,
E.
,
2003
, “
Modelling Heart Tissue Using a Composite Muscle Model With Blood Perfusion
,”
Computational Fluid and Solid Mechanics
,
Elsevier
, Boston, MA, pp.
1642
1646
.
22.
Ateshian
,
G. A.
,
Costa
,
K. D.
, and
Hung
,
C. T.
,
2007
, “
A Theoretical Analysis of Water Transport Through Chondrocytes
,”
Biomech. Model. Mechanobiol.
,
6
(
1–2
), pp.
91
101
.10.1007/s10237-006-0039-9
23.
Barocas
,
V. H.
, and
Tranquillo
,
R. T.
,
1997
, “
An Anisotropic Biphasic Theory of Tissue-Equivalent Mechanics: The Interplay Among Cell Traction, Fibrillar Network Deformation, Fibril Alignment, and Cell Contact Guidance
,”
ASME J. Biomech. Eng.
,
119
(
2
), pp.
137
145
.10.1115/1.2796072
24.
Guilak
,
F.
, and
Mow
,
V. C.
,
2000
, “
The Mechanical Environment of the Chondrocyte: A Biphasic Finite Element Model of Cell–Matrix Interactions in Articular Cartilage
,”
J. Biomech.
,
33
(
12
), pp.
1663
1673
.10.1016/S0021-9290(00)00105-6
25.
Causin
,
P.
,
Guidoboni
,
G.
,
Harris
,
A.
,
Prada
,
D.
,
Sacco
,
R.
, and
Terragni
,
S.
,
2014
, “
A Poroelastic Model for the Perfusion of the Lamina Cribrosa in the Optic Nerve Head
,”
Math. Biosci.
,
257
, pp.
33
41
.10.1016/j.mbs.2014.08.002
26.
Chen
,
X.
,
Dunn
,
A. C.
,
Sawyer
,
W. G.
, and
Sarntinoranont
,
M.
,
2007
, “
A Biphasic Model for Micro-Indentation of a Hydrogel-Based Contact Lens
,”
ASME J. Biomech. Eng.
,
129
(
2
), pp.
156
163
.10.1115/1.2472373
27.
Tandon
,
P. N.
, and
Autar
,
R.
,
1991
, “
Biphasic Model of the Trabecular Meshwork in the Eye
,”
Med. Biol. Eng. Comput.
,
29
(
3
), pp.
281
290
.10.1007/BF02446710
28.
Spilker
,
R. L.
, and
Maxian
,
T. A.
,
1990
, “
A Mixed-Penalty Finite Element Formulation of the Linear Biphasic Theory for Soft Tissues
,”
Int. J. Numer. Methods Eng.
,
30
(
5
), pp.
1063
1082
.10.1002/nme.1620300508
29.
Suh
,
J.-K.
,
Spilker
,
R.
, and
Holmes
,
M.
,
1991
, “
A Penalty Finite Element Analysis for Nonlinear Mechanics of Biphasic Hydrated Soft Tissue Under Large Deformation
,”
Int. J. Numer. Methods Eng.
,
32
(
7
), pp.
1411
1439
.10.1002/nme.1620320704
30.
Almeida
,
E. S.
, and
Spilker
,
R. L.
,
1997
, “
Mixed and Penalty Finite Element Models for the Nonlinear Behavior of Biphasic Soft Tissues in Finite Deformation—Part I: Alternate Formulations
,”
Comput. Methods Biomech. Biomed. Eng.
,
1
(
1
), pp.
25
46
.10.1080/01495739708936693
31.
Simon
,
B. R.
,
Kaufmann
,
M. V.
,
McAfee
,
M. A.
, and
Baldwin
,
A. L.
,
1998
, “
Porohyperelastic Finite Element Analysis of Large Arteries Using ABAQUS
,”
ASME J. Biomech. Eng.
,
120
(
2
), pp.
296
298
.10.1115/1.2798315
32.
Ateshian
,
G. A.
, and
Ricken
,
T.
,
2010
, “
Multigenerational Interstitial Growth of Biological Tissues
,”
Biomech. Model. Mechanobiol.
,
9
(
6
), pp.
689
702
.10.1007/s10237-010-0205-y
33.
Maas
,
S. A.
,
Ellis
,
B. J.
,
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2012
, “
FEBio: Finite Elements for Biomechanics
,”
ASME J. Biomech. Eng.
,
134
(
1
), p. 011005.10.1115/1.4005694
34.
Hou
,
J. S.
,
Holmes
,
M. H.
,
Lai
,
W. M.
, and
Mow
,
V. C.
,
1989
, “
Boundary Conditions at the Cartilage-Synovial Fluid Interface for Joint Lubrication and Theoretical Verifications
,”
ASME J. Biomech. Eng.
,
111
(
1
), pp.
78
87
.10.1115/1.3168343
35.
Chan
,
B.
,
Donzelli
,
P. S.
, and
Spilker
,
R. L.
,
2000
, “
A Mixed-Penalty Biphasic Finite Element Formulation Incorporating Viscous Fluids and Material Interfaces
,”
Ann. Biomed. Eng.
,
28
(
6
), pp.
589
597
.10.1114/1.1305529
36.
Badia
,
S.
,
Quaini
,
A.
, and
Quarteroni
,
A.
,
2009
, “
Coupling Biot and Navier–Stokes Equations for Modelling Fluid–Poroelastic Media Interaction
,”
J. Comput. Phys.
,
228
(
21
), pp.
7986
8014
.10.1016/j.jcp.2009.07.019
37.
Unnikrishnan
,
G.
,
Unnikrishnan
,
V.
, and
Reddy
,
J.
,
2009
, “
Tissue–Fluid Interface Analysis Using Biphasic Finite Element Method
,”
Comput. Methods Biomech. Biomed. Eng.
,
12
(
2
), pp.
165
172
.10.1080/10255840802372045
38.
Bukac
,
M.
,
Yotov
,
I.
,
Zakerzadeh
,
R.
, and
Zunino
,
P.
,
2015
, “
Effects of Poroelasticity on Fluid-Structure Interaction in Arteries: A Computational Sensitivity Study
,”
MS&A
,
Springer International Publishing
, Cham, Switzerland, pp.
197
220
.
39.
Tully
,
B.
, and
Ventikos
,
Y.
,
2009
, “
Coupling Poroelasticity and CFD for Cerebrospinal Fluid Hydrodynamics
,”
IEEE Trans. Biomed. Eng.
,
56
(
6
), pp.
1644
1651
.10.1109/TBME.2009.2016427
40.
Berger
,
L.
,
Bordas
,
R.
,
Burrowes
,
K.
,
Grau
,
V.
,
Tavener
,
S.
, and
Kay
,
D.
,
2016
, “
A Poroelastic Model Coupled to a Fluid Network With Applications in Lung Modelling
,”
Int. J. Numer. Methods Biomed. Eng.
,
32
(
1
), pp.
1
17
.10.1002/cnm.2731
41.
Hou
,
J.
,
Mow
,
V.
,
Lai
,
W.
, and
Holmes
,
M.
,
1992
, “
An Analysis of the Squeeze-Film Lubrication Mechanism for Articular Cartilage
,”
J. Biomech.
,
25
(
3
), pp.
247
259
.10.1016/0021-9290(92)90024-U
42.
Mow
,
V. C.
, and
Lai
,
W. M.
,
1980
, “
Recent Developments in Synovial Joint Biomechanics
,”
SIAM Rev.
,
22
(
3
), pp.
275
317
.10.1137/1022056
43.
Bowen
,
R. M.
,
1982
, “
Compressible Porous Media Models by Use of the Theory of Mixtures
,”
Int. J. Eng. Sci.
,
20
(
6
), pp.
697
735
.10.1016/0020-7225(82)90082-9
44.
Shim
,
J. J.
, and
Ateshian
,
G. A.
,
2021
, “
A Hybrid Biphasic Mixture Formulation for Modeling Dynamics in Porous Deformable Biological Tissues
,”
Arch. Appl. Mech.
, epub, pp.
1
21
.10.1007/s00419-020-01851-8
45.
Maas
,
S. A.
,
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2017
, “
FEBio: History and Advances
,”
Annu. Rev. Biomed. Eng.
,
19
(
1
), pp.
279
299
.10.1146/annurev-bioeng-071516-044738
46.
Ateshian
,
G. A.
,
Shim
,
J. J.
,
Maas
,
S. A.
, and
Weiss
,
J. A.
,
2018
, “
Finite Element Framework for Computational Fluid Dynamics in FEBio
,”
ASME J. Biomech. Eng.
,
140
(
2
), p.
021001
.10.1115/1.4038716
47.
Shim
,
J. J.
,
Maas
,
S. A.
,
Weiss
,
J. A.
, and
Ateshian
,
G. A.
,
2019
, “
A Formulation for Fluid–Structure Interactions in FEBio Using Mixture Theory
,”
ASME J. Biomech. Eng.
,
141
(
5
), p.
051010
.10.1115/1.4043031
48.
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2010
, “
Anisotropic Hydraulic Permeability Under Finite Deformation
,”
ASME J. Biomech. Eng.
,
132
(
11
), p.
111004
.10.1115/1.4002588
49.
Pierce
,
D. M.
,
Ricken
,
T.
, and
Holzapfel
,
G. A.
,
2013
, “
A Hyperelastic Biphasic Fibre-Reinforced Model of Articular Cartilage Considering Distributed Collagen Fibre Orientations: Continuum Basis, Computational Aspects and Applications
,”
Comput. Methods Biomech. Biomed. Eng.
,
16
(
12
), pp.
1344
1361
.10.1080/10255842.2012.670854
50.
Ateshian
,
G. A.
,
Albro
,
M. B.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2011
, “
Finite Element Implementation of Mechanochemical Phenomena in Neutral Deformable Porous Media Under Finite Deformation
,”
ASME J. Biomech. Eng.
,
133
(
8
), p.
081005
.10.1115/1.4004810
51.
Ateshian
,
G. A.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2013
, “
Multiphasic Finite Element Framework for Modeling Hydrated Mixtures With Multiple Neutral and Charged Solutes
,”
ASME J. Biomech. Eng.
,
135
(
11
), p.
111001
.10.1115/1.4024823
52.
Brinkman
,
H.
,
1949
, “
A Calculation of the Viscous Force Exerted by a Flowing Fluid on a Dense Swarm of Particles
,”
Flow, Turbul. Combust.
,
1
(
1
), p.
27
.10.1007/BF02120313
53.
Bonet
,
D. J.
, and
Wood
,
D. R. D.
,
2008
,
Nonlinear Continuum Mechanics for Finite Element Analysis
,
Cambridge University Press
, Cambridge, UK.
54.
Chung
,
J.
, and
Hulbert
,
G. M.
,
1993
, “
A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized-Alpha Method
,”
ASME J. Appl. Mech.
,
60
(
2
), pp.
371
375
.10.1115/1.2900803
55.
Jansen
,
K. E.
,
Whiting
,
C. H.
, and
Hulbert
,
G. M.
,
2000
, “
A Generalized-Alpha Method for Integrating the Filtered Navier–Stokes Equations With a Stabilized Finite Element Method
,”
Comput. Methods Appl. Mech. Eng.
,
190
(
3–4
), pp.
305
319
.10.1016/S0045-7825(00)00203-6
56.
Engelman
,
M. S.
,
Strang
,
G.
, and
Bathe
,
K.-J.
,
1981
, “
The Application of Quasi-Newton Methods in Fluid Mechanics
,”
Int. J. Numer. Methods Eng.
,
17
(
5
), pp.
707
718
.10.1002/nme.1620170505
57.
Schenk
,
O.
, and
Gartner
,
K.
,
2004
, “
Solving Unsymmetric Sparse Systems of Linear Equations With PARDISO
,”
Future Gener. Comput. Syst.
,
20
(
3
), pp.
475
487
.10.1016/j.future.2003.07.011
58.
Ateshian
,
G. A.
,
1991
, “
Biomechanics of Diarthrodial Joints: Applications to the Thumb Carpometacarpal Joint
,” Ph.D. thesis,
Columbia University
, New York.
59.
Bathe
,
K.-J.
, and
Ledezma
,
G. A.
,
2007
, “
Benchmark Problems for Incompressible Fluid Flows With Structural Interactions
,”
Comput. Struct.
,
85
(
11–14
), pp.
628
644
.10.1016/j.compstruc.2007.01.025
60.
Taylor
,
G. I.
,
1971
, “
A Model for the Boundary Condition of a Porous Material—Part I
,”
J. Fluid Mech.
,
49
(
2
), pp.
319
326
.10.1017/S0022112071002088
61.
Holzapfel
,
G. A.
,
Gasser
,
T. C.
, and
Ogden
,
R. W.
,
2000
, “
A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models
,”
J. Elasticity
,
61
(
1/3
), pp.
1
48
.10.1023/A:1010835316564
62.
Sommer
,
G.
, and
Holzapfel
,
G. A.
,
2012
, “
3D Constitutive Modeling of the Biaxial Mechanical Response of Intact and Layer-Dissected Human Carotid Arteries
,”
J. Mech. Behav. Biomed. Mater.
,
5
(
1
), pp.
116
128
.10.1016/j.jmbbm.2011.08.013
63.
Vignon-Clementel
,
I. E.
,
Figueroa
,
C. A.
,
Jansen
,
K. E.
, and
Taylor
,
C. A.
,
2006
, “
Outflow Boundary Conditions for Three-Dimensional Finite Element Modeling of Blood Flow and Pressure in Arteries
,”
Comput. Methods Appl. Mech. Eng.
,
195
(
29–32
), pp.
3776
3796
.10.1016/j.cma.2005.04.014
64.
Hou
,
C.
, and
Ateshian
,
G. A.
,
2016
, “
A Gauss-Kronrod-Trapezoidal Integration Scheme for Modeling Biological Tissues With Continuous Fiber Distributions
,”
Comput. Methods Biomech. Biomed. Eng.
,
19
(
8
), pp.
883
893
.10.1080/10255842.2015.1075518
65.
Puso
,
M. A.
, and
Weiss
,
J. A.
,
1998
, “
Finite Element Implementation of Anisotropic Quasi-Linear Viscoelasticity Using a Discrete Spectrum Approximation
,”
ASME J. Biomech. Eng.
,
120
(
1
), pp.
62
70
.10.1115/1.2834308
66.
Nims
,
R. J.
,
Durney
,
K. M.
,
Cigan
,
A. D.
,
Dusséaux
,
A.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2016
, “
Continuum Theory of Fibrous Tissue Damage Mechanics Using Bond Kinetics: Application to Cartilage Tissue Engineering
,”
Interface Focus
,
6
(
1
), p.
20150063
.10.1098/rsfs.2015.0063
67.
Cho
,
Y. I.
, and
Kensey
,
K. R.
,
1991
, “
Effects of the Non-Newtonian Viscosity of Blood on Flows in a Diseased Arterial Vessel—Part I: Steady Flows
,”
Biorheology
,
28
(
3–4
), pp.
241
262
.10.3233/BIR-1991-283-415
68.
Ateshian
,
G. A.
,
Nims
,
R. J.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2014
, “
Computational Modeling of Chemical Reactions and Interstitial Growth and Remodeling Involving Charged Solutes and Solid-Bound Molecules
,”
Biomech. Model. Mechanobiol.
,
13
(
5
), pp.
1105
1120
.10.1007/s10237-014-0560-1