## 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 [1–4], leading to significant advances in modeling the mechanics of biological tissues and cells [5–9]. 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,13–16], brain [17], lung [18], liver [19], heart [20,21], cells and their extracellular matrix [22–24], and the eye [25–27], among many others. The increasing use of biphasic materials led to the development of related computational tools [28–33].

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 [35–37]. 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 febio^{1} 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 $\Omega 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 $\Omega f$, and solid domains $\Omega s$ across interfaces $\Gamma bb,\u2009\Gamma bf$, and $\Gamma 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 $\Omega 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].

### 2.1 Governing Equations.

We model the domain $\Omega 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 *α* ($\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 *α* in the elemental mixture volume *dV*. This apparent density is related to the true density $\rho T\alpha =dm\alpha /dV\alpha $ (mass of *α* per volume of *α*) via $\rho \alpha =\varphi \alpha \rho T\alpha $, where $\varphi \alpha =dV\alpha /dV$ is the volume fraction of *α* in the mixture, satisfying $\varphi s+\varphi f=1$. Since the solid skeleton is intrinsically incompressible, $\rho 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 (*dV _{r}*) configurations. Thus, since the solid skeleton is intrinsically incompressible,

*J*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 rs/Js$, where $\rho rs=dms/dVr$ is the solid apparent density in the mixture reference configuration. Therefore, we may define the referential solid volume fraction as $\varphi rs=\rho rs/\rho Ts$. In the finite element analysis, $\varphi rs=dVs/dVr$ and $\rho Ts$ are both specified by the user, then $\rho rs$,

^{s}*ρ*, $\varphi s$, and $\varphi f$ are evaluated from the above relations, given a solution for $Fs$ (and thus,

^{s}*J*). The solution for $Fs=I+Gradu$ is obtained by solving for the nodal solid displacement vector

^{s}**u**, where $Grad(\xb7)$ represents the gradient operator in the material frame.

*J*is the volume ratio of the fluid in its current and reference configurations,

^{f}*dm*is the element mass of fluid, and

^{f}*dV*is the elemental volume of that fluid in the pores of the mixture, in the current configuration. It follows from this definition that:

^{f}*dm*) in the fluid's reference configuration. Thus, $\rho 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

^{f}*J*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 $\varphi f\rho Tf$, which now expands into an expression involving solid and fluid volume ratios

^{f}The limiting case when the fluid within the solid matrix pores has been completely squeezed out corresponds to $\rho 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 *dV ^{s}* in the current configuration, or equivalently $Js=\varphi 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 *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.

**u**), and

In this expression, the dot operator in $J\u02d9f$ 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\u02d9$. 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.

*α*,

*p*is the fluid pressure, $\tau f$ is the apparent fluid viscous stress, and $\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

*α*per elemental

*mixture*area. Here,

**k**is the hydraulic permeability tensor which regulates frictional drag between the fluid and solid constituents; setting $k\u22121$ 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

*e*, and its simplest form is

^{f}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 ($\tau f=0$), we recover Darcy's law, $w=\u2212k\xb7gradp$.

### 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, $\Gamma bb$, and the interfaces $\Gamma bf$ and $\Gamma 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+\u2212f\u2212$ to denote the jump in the function *f* across the interface Γ, with $f+$ and $f\u2212$ 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$.

*ψ*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

^{f}*J*; therefore, this energy jump condition implies that

^{f}*J*must be continuous across Γ, thus

^{f}**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

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 $\tau f=\varphi f\tau $ where $\tau $ would be the true fluid viscous stress. Thus, we can use febio's various constitutive relations for the viscous stress $\tau $ of Newtonian or non-Newtonian fluids and adapt those models to a biphasic mixture where $\tau f$ is evaluated as $\varphi f\tau $. Accordingly, the contribution of $\tau f$ would properly reduce to zero in the limit as fluid content reduces to zero ($\varphi f\u21920$), in which case the mixture momentum jump in Eq. (2.15) would reduce to $[[\u2212pI+\sigma e]]\xb7n=[[\sigma e]]\xb7n=0$, since $[[p]]=0$.

**w**and

*J*be nodal DOFs automatically enforces the jump conditions (2.13) and (2.14), acting as essential boundary conditions, along with the solid displacement

^{f}**u**. Finally, subtracting Eq. (2.16) from Eq. (2.15), we obtain the momentum jump condition for the solid constituent

A BFSI domain may be reduced to a FSI domain by letting $\varphi rs\u21920$ and $k\u22121\u21920$; 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 $\Gamma bb$ are equally valid for interfaces $\Gamma bf$.

However, when the biphasic domain interfaces with a nonporous solid domain across $\Gamma bs$, the jump conditions (2.13) on *J ^{f}* and (2.17) on $\tau 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 $\Gamma bs$, although the user would have to enforce this no-slip condition explicitly by prescribing it as an essential boundary condition. The mixture momentum jump $[[\u2212pI+\sigma e+\tau f]]\xb7n=0$ implies that the mixture traction on the $\Omega b$ side of $\Gamma bs$ is equal to the solid traction on the $\Omega s$ side. This jump condition would also need to be enforced explicitly, as described below.

### 2.3 Virtual Work and Weak Form.

**u**,

**w**, and

*e*, 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 $\tau f=\varphi f\tau $, these become

^{f}Here, $te=\sigma e\xb7n$ is the elastic traction, $t\tau =\tau \xb7n$ is the true fluid viscous traction, and $wn=wn\xb7n$ is the normal component of the relative fluid flux on the boundary $\u2202\Omega b$, whose outward unit normal is **n**. The traction $t\sigma $ 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 $\u2202\Omega b$, the natural boundary conditions are $t\sigma =0$, $t\tau =0$, and *w _{n}* = 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

*e*, which are also consistent with the above jump conditions. In particular, an essential no-slip boundary condition may be prescribed on $\Gamma bs$ by setting $w=0$. A symmetry plane may be prescribed with the essential boundary condition $un\u2261u\xb7n=0$ and the natural boundary conditions $t\tau =0$ and

^{f}*w*= 0.

_{n}In this formulation, the mixture traction is defined as $t=\u2212pn+te+\varphi ft\tau $, which may also be written as $t=\u2212pn+t\sigma +t\tau $. 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\sigma $ and $t\tau $ may be prescribed as natural boundary conditions, whereas *p* may be prescribed as an essential boundary condition on *e ^{f}*, using Eq. (2.9). However, there are two general scenarios where

**t**needs to be prescribed on a region of the boundary $\u2202\Omega b$ with incomplete prior knowledge of

*p*, $t\sigma $, or $t\tau $: (1) When a BFSI boundary $\Gamma 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 $\u2212pn+t\sigma +t\tau =0$ as a constraint equation on that boundary, to impart the free surface $\Gamma b$ its natural shape. (2) At a biphasic–solid interface $\Gamma bs$,

**t**must balance the traction acting on the adjoining solid domain. Since

**u**is continuous across $\Gamma bs$ due to shared nodes, the solid natural boundary condition $t\sigma $ of

**t**is already accounted for by the deformation, so that it is only necessary to prescribe the portion $\u2212pn+t\tau $ of

**t**on the solid domain; thus, $\u2212pn+t\tau +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 $\Gamma b$ and $\Gamma bs$, respectively.

*da*on Γ may be evaluated from the covariant basis vectors $g\alpha $ ($\alpha =1,2$)

**n**to $\Omega b$ on Γ is evaluated from

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

### 2.4 Linearization and Spatial Discretization.

where the operator $D\delta W[\xb7]$ represents the directional derivative of *δW* at $(u,w,Jf)$ along increments $\Delta u$ of **u**, $\Delta w$ of **w**, and $\Delta ef$ of *e ^{f}* [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 *t _{n}* when the vector $\Delta U(k)$ of nodal DOF increments, consisting of $\Delta u,\u2009\Delta w$, and $\Delta ef$ values at all unconstrained nodes, at the

*k*th iteration satisfies $||\Delta U(k)||\u2264\epsilon r||\Delta U(1)||$, where $\epsilon r$ represents the relative tolerance criterion, which is typically set to $\epsilon r=10\u22123$ 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.

^{3}febio 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 $k\u22121$ 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 $\Gamma 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 $\Gamma 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=1\u2009mm$ along the $z$-direction ($0\u2264z\u2264h$) was subjected to a prescribed linear ramp-and-hold displacement at $z=1\u2009\u2009mm$, producing a compression of $10\u2009%$ in $100\u2009s$, and then allowed to relax to time $t=1000\u2009s$. 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 $\varphi rs=0.2$. The porous solid matrix was a compressible neo-Hookean material [53] with $E=1\u2009kPa$, and *ν* = 0. The hydraulic permeability tensor was isotropic, $k=kI$, with constant $k=0.001\u2009mm4/mN\u2009s$. For the BFSI material, the fluid constituent had a bulk modulus $K=1\xd71012\u2009kPa$ 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 $\Delta t=2.5\u2009s$ 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 *u _{z}* and fluid pressure

*p*over the range of

*z*at selected time points

*t*demonstrated agreement between the two solvers (Fig. 2).

### 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 *w _{y}* of a single sawtooth profile of amplitude $0.01\u2009m/s$ and total duration of $2\u2009s$. Here, we only looked at the case when the generalized$\u2212\alpha $ spectral radius parameter was set to $\rho \u221e=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=10\u22129\u2009Pa$ and Poisson's ratio

*ν*= 0. For the BFSI analysis, the fluid domain was additionally assigned $\varphi s=0$ and $k\u22121=0\u2009m4/N\u2009s$. In all cases, the models had $27\u2009,100$ nodes, $13\u2009,300$ elements, and $80,\u2009394$ 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,

*C*and

_{D}*C*, 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.

_{L}### 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, $\Omega topf$ and $\Omega botf$, were separated by a thin Mooney–Rivlin elastic solid domain, $\Omega s$, discretized with solid (hexahedral) elements (Fig. 4). The left, right, and top surfaces of $\Omega topf$ were exposed to ambient conditions (zero pressure), whereas the bottom surface of $\Omega botf$ was subjected to a fluid pressure that increased linearly with time. No-slip conditions were imposed on the bottom wall of $\Omega topf$, the side walls of $\Omega botf$, and the fluid–solid interfaces $\Gamma topfs$ and $\Gamma botfs$ on both sides of $\Omega s$. The model was run with time-step $\Delta t=0.001\u2009\u2009s$ until $t=0.25\u2009\u2009s$. For the BFSI analysis, the fluid domains were modeled as BFSI materials with $\varphi s=0$ and $k\u22121=0\u2009m4/N\u2009s$; the solid matrix was modeled as a compressible neo-Hookean material with Young's modulus $E=10\u22129\u2009Pa$ and Poisson's ratio *ν* = 0. The febio model had $24\u2009,348$ nodes, $11,\u2009882$ elements, and $118\u2009,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).

### 3.4 Taylor Slip Problem.

*v*

_{0}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

where *μ _{f}* and

*μ*are the true fluid viscosities in the fluid and biphasic layers, respectively, and

*h*

_{1}and

*h*

_{2}are the heights of the fluid and porous layers, respectively. Four cases were examined using different values for

*δ*and

*η*, with $\xi =0.25$ and $\varphi 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=1012\u2009Pa$. 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

*w*was normalized by

_{x}*v*

_{0}. All febio solutions agreed with corresponding analytical solutions (Fig. 7).

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

*δ*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

where *μ _{s}* is the shear modulus of the porous solid in the biphasic layer. We examined cases where $\delta 2=0.01,\u2009\delta 2=0.1$, and $\delta 2=1$, while $\xi =0.1$,

*η*= 1, and $\alpha =0.01$. We also varied $\eta =0.1$,

*η*= 1, and

*η*= 10, while $\xi =0.1$, $\delta 2=0.01$, and $\alpha =0.01$. We set $\varphi f=0.75$ in all cases. The mesh was modified to also be biased across the interface $\Gamma 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).

### 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 $\Omega tops$ meets a viscous fluid layer $\Omega topf$ at the interface $\Gamma topfs$ located at $z=\u2212h1$, where $\Omega topf$ then interfaces a biphasic layer $\Omega botb$ at $\Gamma bf$ located at *z *=* *0. The bottom of the biphasic layer is in contact with a stationary free-draining rigid porous platen $\Omega bots$ at the interface $\Gamma botbs$ located at $z=h2$. The rigid platen $\Omega tops$ was prescribed either a step velocity, *v*_{0}, or step pressure, *p*_{0}, applied at $t=0+$.

The problem was modeled in febio using $K=1014\u2009Pa$ to nearly enforce fluid incompressibility, biphasic porosity $\varphi f=0.7$, upstream true fluid viscosity $\mu f=1\u2009Pa\xb7s$, biphasic true fluid viscosity $\mu =1\u2009Pa\xb7s$, biphasic solid Young's modulus $E=9\xd7105\u2009Pa$, Poisson's ratio $\nu =0.2$ (producing $HA=106\u2009Pa$), and hydraulic permeability $k=10\u221214\u2009m4/N\u2009s$. Either a step pressure $p0=104\u2009Pa$ or a step velocity $v0=10\u22128m/s$ was prescribed. The mesh had $1276$ nodes, 570 elements, and $3675$ DOFs, biased at all interfaces. The model ran up to $t=1000\u2009s$, where the time-step size started at $\Delta t=0.001\u2009s$ and gradually increased to $\Delta t=1\u2009s$ 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=1\u2009s$, $t=10\u2009s,\u2009t=100\u2009s$, and $t=1000\u2009s$, as reported in Fig. 10, showing good agreement. Note that in the step pressure case, the value of *v*_{0} used in calculating the normalized fluid velocity was obtained from Eq. (A15) as $t\u2192\u221e$.

### 3.7 Dilatational Wave Propagation.

*c*in the fluid and

^{f}*c*in the solid

^{s}where *H _{A}* is defined in the Appendix. Here, we verified the existence of these two wave propagations, as well as the accuracy of the formulas for

*c*and

^{f}*c*in febio, using the BFSI solver.

^{s}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 $1\u2009m$ along the $x$-direction ($0\u2264x\u22641$), height $0.1\u2009m$ along *y*, and depth $0.05\u2009m$ along *z*. The prescribed properties were $\varphi s=0.2$, zero true fluid viscosity, $\rho Trf=1.225\u2009kg/m3$, and $K=101,325\u2009Pa$. The porous solid was a compressible neo-Hookean hyperelastic material with $\rho Ts=103\u2009kg/m3$ and $E=0.16\u2009GPa$. Here, $k\u22121$ was set to zero to eliminate frictional drag between the fluid and solid constituents. At $x=1\u2009\u2009m$, solid displacement *u _{x}*was fixed, while natural boundary conditions were enforced for the fluid (

*w*= 0). At

_{x}*x*=

*0, we prescribed either a solid displacement $ux=\u2212(1/2)u0(1\u2212cos\u20092\pi t/T)$ or a normal fluid flux $wx=\u2212(1/2)w0(1\u2212cos\u20092\pi t/T)$, with $u0=10\u22126\u2009m$, $w0=\u2212103\u2009m/s$, and $T=3\xd710\u22124\u2009s$. In the case when*

*w*was prescribed,

_{x}*u*was subsequently fixed at

_{x}*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.006\u2009s$. The generalized-*α* parameter $\rho \u221e$ 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 *c ^{s}* and

*c*, 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\u02d9x(x,t)$. Disturbances reflected at both ends of the domain, with $u\u02d9x(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.6\u2009m/s$ and $cs=894.4\u2009m/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.

^{f}### 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, $\Omega lmnf$, had a height $h2=0.1\u2009m$ and length $l=0.2\u2009m$, connecting with the thin biphasic face mask domain, $\Omega maskb$, which had a thickness of $t=0.002\u2009m$. This mask domain interfaced with a larger ambient air domain, $\Omega ambf$, which had height $H=1\u2009m$ and length $L=2\u2009m$. All domains had a uniform depth of $0.01\u2009m$. A plug velocity profile *v _{n}* was prescribed on the leftmost boundary of $\Omega lmnf$ to replicate cyclical breathing conditions. The exact temporal profile of

*v*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.5\u2009s$. To produce a physiological inspiratory:expiratory ratio, the first $0.9\u2009s$ consisted of inspiration and the remainder of the cycle was expiration. A plot of the prescribed inlet velocity

_{n}*v*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 $\Delta t=0.00025\u2009s$. Zero fluid pressure was prescribed on the boundaries of $\Omega ambf$ as shown in Fig. 12. For the standard face mask, the hydraulic permeability

_{n}*k*was set to $3\xd710\u22126\u2009m4/(N\u2009s)$, calculated to produce a pressure drop of approximately $600\u2009Pa$ across the mask, which is a value typically found in a mechanical ventilator. This value of

*k*was also compared to $k=3\xd710\u22124\u2009m4/(N\u2009s)$ and $k\u22121=0$ to investigate the effects of a more porous mask and no mask, respectively. The air was modeled as Newtonian viscous with $\rho Trf=1.2047\u2009kg/m3,\u2009K=101,325\u2009\u2009Pa$, and $\mu =1.8205\xd710\u22125\u2009Pa\xb7s$. For nonzero $k\u22121$, the properties of the solid constituent of the mask were approximated as $\rho Ts=80\u2009kg/m3$, $E=108\u2009Pa$, and $\nu =0.4$, with $\varphi rs=0.75$. Tangential and backflow stabilization boundary conditions were prescribed on the zero fluid pressure boundaries of $\Omega ambf$. For all models, there were $109,066$ nodes, $53,982$ elements, and $540,284$ DOFs, in a mesh that was biased across interfaces $\Gamma bf$ and normal to no-slip walls to capture any sharp gradients in the fluid flow or pressure.

The vorticity in the $z$-direction and the magnitude of **w** are displayed at two selected time points in Fig. 14. Time points $t=8\u2009s$ and $t=11.5\u2009s$ occur during inspiration of the fourth breath and expiration of the fifth breath, respectively. For the standard face mask, the fluid pressure ranged from $\u2212724\u2009Pa$ to $791\u2009Pa$ over the entire domain and duration of the analysis. For the face mask with higher *k*, the corresponding range was from $\u221213.5\u2009Pa$ to $10.4\u2009Pa$, and in the no-mask case it was $\u22127.95\u2009Pa$ to $5.49\u2009Pa$. 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.

From a numerical perspective, it is interesting to note that small flow disturbances appear whenever shedded vortices cross the zero fluid pressure boundaries of $\Omega ambb$, as is apparent at the right boundary at $t=8\u2009\u2009s$, and the bottom and right boundaries at time $t=11.5\u2009\u2009s$, 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, $\Omega bldf$, as an FSI domain to account for its deformation. The inlet, $\Gamma inf$, had a diameter of $6.28\u2009mm$, and the two outlets, $\Gamma lgf$ and $\Gamma smf$, had respective diameters of $4.24\u2009mm$ and $3.04\u2009mm$. The arterial wall was created by extruding the outer surface of the artery twice, with the tunica media given a thickness of $0.3\u2009mm$ and the tunica adventitia a thickness of $0.15\u2009mm$, for a total arterial wall thickness of $0.45\u2009mm$. These domains correspond to $\Omega medb$ and $\Omega 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,\u2009300\u2009Pa$, $k1=24,700\u2009Pa$, and $k2=16.5$ in the media, and $c=59,600\u2009Pa$, $k1=180,900\u2009Pa$, and $k2=109.8$ in the adventitia. The fiber families formed an angle of $\xb16.9\u2009deg$ in the media and $\xb130.1\u2009deg$ 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 $\rho Ts=1000\u2009kg/m3$ for the dynamic response of the solid constituents of the domains $\Omega medb$ and $\Omega advb$, with solid bulk modulus $108\u2009Pa$ to enforce a nearly incompressible response. Both blood and the interstitial fluid within the arterial wall were modeled as Newtonian fluids where $\rho Trf=1000\u2009kg/m3$, $K=2.2\xd7109\u2009Pa$, and $\mu =0.004\u2009Pa\xb7s$.

For the boundary conditions, an inlet velocity $w=wnn$ was prescribed on $\Gamma inf$ having a parabolic spatial profile, and an average value $w0(t)$ that rises from a diastolic plateau of $0.10\u2009m/s$ to a systolic maximum of $0.48\u2009m/s$ in three consecutive cycles such that $Re=165\u2212800$, 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 $\Gamma lgf$ and $\Gamma smf$, with $R=4\xd7108\u2009kg/m4\u2009s$ and a pressure offset of $p=104\u2009Pa$. The inlet was initially ramped up smoothly from zero to the diastolic value in $0.05\u2009s$, while the outlet pressure offset was ramped up smoothly over $0.25\u2009s$, 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 $\Gamma lgf$ and $\Gamma smf$. A fluid dilatation corresponding to a pressure of $p=104\u2009Pa$ 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.25\u2009s$. The outer surface of the adventitia, $\Gamma 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 $\Delta t=1\u2009ms$.

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\xd710\u221216\u2009m4/(N\u2009s)$ [31], while in the second case we set $k=5.26\xd710\u221211\u2009m4/(N\u2009s)$ to exaggerate the permeability of the arterial wall. In both cases, the referential solid volume fraction was set to $\varphi 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.87\u2009s)$ 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, $\u222b\Gamma \rho fw\xb7nda$, across each of the boundaries of the fluid domain $\Omega bldf$, including the arterial wall, $\Gamma wallb$, along with the time rate of change of total fluid mass $m=\u222b\Omega \rho fdv$ in the overall arterial lumen $\Omega 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\u2009%$ 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.

## 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 *e ^{f}*, 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 *e ^{f}*. 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 $\Gamma bb$ and $\Gamma 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 $\u223c10\u2009h$ while the realistic case took $\u223c24\u2009h$. 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 $\Gamma bb,\u2009\Gamma bf$, or $\Gamma 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

*μ*is the fluid shear viscosity. The governing equations now reduce to

_{a}*z*=

*0 are*