The mechanics of biological fluids is an important topic in biomechanics, often requiring the use of computational tools to analyze problems with realistic geometries and material properties. This study describes the formulation and implementation of a finite element framework for computational fluid dynamics (CFD) in FEBio, a free software designed to meet the computational needs of the biomechanics and biophysics communities. This formulation models nearly incompressible flow with a compressible isothermal formulation that uses a physically realistic value for the fluid bulk modulus. It employs fluid velocity and dilatation as essential variables: The virtual work integral enforces the balance of linear momentum and the kinematic constraint between fluid velocity and dilatation, while fluid density varies with dilatation as prescribed by the axiom of mass balance. Using this approach, equal-order interpolations may be used for both essential variables over each element, contrary to traditional mixed formulations that must explicitly satisfy the inf-sup condition. The formulation accommodates Newtonian and non-Newtonian viscous responses as well as inviscid fluids. The efficiency of numerical solutions is enhanced using Broyden's quasi-Newton method. The results of finite element simulations were verified using well-documented benchmark problems as well as comparisons with other free and commercial codes. These analyses demonstrated that the novel formulation introduced in FEBio could successfully reproduce the results of other codes. The analogy between this CFD formulation and standard finite element formulations for solid mechanics makes it suitable for future extension to fluid–structure interactions (FSIs).

# Finite Element Framework for Computational Fluid Dynamics in FEBio OPEN ACCESS

**Gerard A. Ateshian**,

**Jay J. Shim**

**Steve A. Maas**,

**Jeffrey A. Weiss**

Manuscript received May 4, 2017; final manuscript received December 5, 2017; published online January 12, 2018. Editor: Victor H. Barocas.

*J Biomech Eng*140(2), 021001 (Jan 12, 2018) Paper No: BIO-17-1193; doi: 10.1115/1.4038716 History: Received May 04, 2017; Revised December 05, 2017

The mechanics of biological fluids is an important topic in biomechanics, most notably for the study of blood flow through the cardiovascular system and related biomedical devices, cerebrospinal fluid flow, airflow through the respiratory system, biotribology by fluid film lubrication, and flow through microfluidic biomedical devices. Therefore, the application domain of multiphysics computational frameworks geared toward biomechanics and biophysics, such as the free finite element software FEBio,^{1} can be expanded by incorporating solvers for computational fluid dynamics (CFD). Biological fluids are generally modeled as incompressible materials, though compressible flow may be needed to analyze wave propagation, for example, in acoustics (airflow through vocal folds) and the analysis of ultrasound propagation for imaging blood flow.

Many open-source CFD codes are currently available in the public domain, though few are geared specifically for applications in biomechanics. OpenFOAM^{2} has many solvers applicable to a very broad range of fluid analyses, using the finite volume method. Other open-source CFD codes are typically more specialized: SU2^{3} is geared toward aerospace design, Fluidity^{4} is well suited for geophysical fluid dynamics, and REEF3D^{5} is focused on marine applications. Some CFD codes explore alternative solution methods: Palabos^{6} uses the lattice Boltzmann method, LIGGGHTS^{7} uses a discrete-element method particle simulation, the Gerris Flow Solver^{8} uses an octree finite volume discretization, and COOLFluiD^{9} uses either a spectral finite difference solver or a finite element solver. Currently, the CFD code most relevant to biomechanics is SimVascular, a finite element code specifically designed for cardiovascular fluid mechanics, “providing a complete pipeline from medical image data segmentation to patient-specific blood flow simulation and analysis.”^{10} This open-source code [1] is based on the original work of Taylor [2]; it uses the flow solver from the PHASTA project^{11} [3].

In contrast, FEBio was originally developed with a focus on the biomechanics of soft tissues, providing the ability to model nonlinear anisotropic tissue responses under finite deformation. It accommodates hyperelastic, viscoelastic, biphasic (poroelastic), and multiphasic material responses which can combine solid mechanics and mass (solvent and solute) transport, with a wide range of constitutive models relevant to biological tissues and cells [4–6]. It also provides robust algorithms to model tied or sliding contact under large deformations between elastic, biphasic, or multiphasic materials [7,8]. More recently, FEBio has incorporated reactive mechanisms, including chemical reactions in multiphasic media, growth mechanics, reactive viscoelasticity, and reactive damage mechanics [9–11]. Many of these features are specifically geared toward applications in biomechanics, including growth and remodeling and mechanobiology.

Our medium-term goal is to expand the capabilities of FEBio by including robust fluid–structure interactions (FSIs). Such interactions occur commonly in biomechanics and biophysics, most notably in cardiovascular mechanics where blood flows through the deforming heart and vasculature [12–14], diarthrodial joint lubrication where pressurized synovial fluid flows between deforming articular layers [15], cerebrospinal mechanics where fluid flow through ventricular cavities may cause significant deformation of surrounding soft tissues [16–18], vocal fold and upper airway mechanics [19,20], viscous flow over endothelial cells [21,22], canalicular and lacunar flow around osteocytes resulting from bone deformation [23–25], and many applications in biomedical device design [26–28]. Some FSI capabilities already exist in several of the open-source CFD codes referenced previously; these implementations were developed as addendums to sophisticated fluid analyses, often employing a limited range of solid material responses, such as small strain analyses of linear isotropic elastic solids. By formulating a CFD code in FEBio, we expect to capitalize on the much broader library of solid and multiphasic material models already available in that framework to considerably expand the range of FSI analyses in biomechanics. Our long-term goal is to expand these FSI features to incorporate chemical reactions within the fluid and at the fluid–solid interface, growth, and remodeling of the solid matrix of multiphasic domains interfaced with the fluid, active transport of solutes across cell membranes as triggered by mechanical, electrical, or chemical signals at fluid–membrane interfaces, cellular and muscle contraction mechanisms, and a host of other mechanisms relevant to biomechanics, biophysics, and mechanobiology, complementing many of the similar features already implemented in FEBio's elastic, biphasic, and multiphasic domain frameworks. The FEBio CFD formulation proposed in this study represents a significant milestone toward that long-term goal.

Since FEBio uses the finite element method, this study formulates a finite element CFD code. Finite element analyses of incompressible flow typically enforce the incompressibility condition using either mixed formulations [29], which solve for the fluid pressure by enforcing zero divergence of the velocity, or the penalty method [30,31], where the fluid pressure is proportional to the divergence of the fluid velocity via a penalty factor [32]. The mixed formulation results in a coefficient matrix that must pass the inf-sup condition to produce accurate results, using a finite element interpolation of the fluid pressure with a lower order than that of the fluid velocity [32–34].

The most broadly adopted CFD schemes today rely on stabilization methods, such as streamline-upwind/Petrov–Galerkin (SUPG), pressure-stabilizing/Petrov–Galerkin and Galerkin/least-squares [35–39]. These schemes have three principal benefits: They stabilize the flow and produce smooth results even on coarse meshes; they allow the use of equal interpolation for the velocity and pressure; and the resulting equations are amenable to iterative solving, with suitable preconditioning. They also have some drawbacks: The finite element formulation is more complex because it requires additional terms to supplement the standard Galerkin residual; these terms involve second-order spatial derivatives, which may be neglected in linear elements [40] but must be included in higher-order interpolations [39,41,42]; and they require the careful formulation of a weight factor, called the stabilization parameter *τ*, whose value depends on the local Reynolds number, some characteristic element length along the streamline direction, and some special considerations for the element shape and interpolation order [41,42]; this parameter critically controls the effectiveness of this stabilization scheme.

In this study, we overcame the constraint of the inf-sup condition by employing a novel approach to analyze isothermal compressible flow using fluid dilatation (the relative change in fluid volume) as an essential variable, instead of fluid pressure or density. This approach also allowed us to forgo the application of stabilization methods, reducing the complexity of the implementation and avoiding the need to formulate and use a stabilization parameter suitable for various element types and orders. Just like fluid velocity, dilatation is a kinematic quantity that may further serve as a state variable in the formulation of functions of state (such as viscous stress and elastic pressure). Another benefit of this approach is that the dependence of density on dilatation is obtained by analytically integrating the mass balance equation, as conventionally done in finite elasticity. In effect, we adapted an approach from finite deformation elasticity to describe the elastic compressibility of the fluid, using the dilatation as the only required component of the strain tensor, since fluids are isotropic and cannot sustain shear stresses under steady shear strains. In this approach, fluid velocity and dilatation may be obtained by satisfying the momentum balance equation as well as the fundamental kinematic relation between the material time derivative of the dilatation and the divergence of the velocity. The finite element implementation, using a standard Galerkin residual formulation, produces accurate numerical solutions even when using equal-order interpolation for the velocity and dilatation.

In a spatial (Eulerian) frame, the momentum balance equation for a continuum is

where *ρ* is the density, $\sigma $ is the Cauchy stress, **b** is the body force per mass, and $a$ is the acceleration, given by the material time derivative of the velocity $v$ in the spatial frame

where $L=grad\u2009v$ is the spatial velocity gradient. The mass balance equation is

where the material time derivative of the density in the spatial frame is

Let **F** denote the deformation gradient (the gradient of the motion with respect to the material coordinate). The material time derivative of **F** is related to $L$ via

Let $J=det\u2009F$ denote the Jacobian of the motion (the volume ratio, or ratio of current to referential volume, *J* > 0); then, the dilatation (relative change in volume between current and reference configurations) is given by $e=J\u22121$. Using the chain rule, *J*'s material time derivative is $J\u02d9=JF\u2212T:F\u02d9$ which, when combined with Eq. (2.5), produces a kinematic constraint between $J\u02d9\u2009and\u2009\u2009div\u2009v$

Substituting this relation into the mass balance, Eq. (2.3), produces $\rho J\xaf\u02d9=0$, which may be integrated directly to yield

where $\rho r$ is the density in the reference configuration (when $J=1$). Since $\rho r$ is obtained by integrating the above-mentioned material time derivative of $\rho J$, it is an intrinsic material property that must be invariant in time and space.

where $I$ is the identity tensor, $\tau $ is the viscous stress, *p* is the pressure arising from the elastic response

and $\Psi r$ is the free energy density of the fluid (free energy per volume of the continuum in the reference configuration). The axiom of entropy inequality dictates that $\Psi r$ cannot be a function of the rate of deformation $D=(L+LT)/2$. In contrast, the viscous stress $\tau $ is generally a function of $J\u2009and\u2009D$.

Boundary conditions may be derived by satisfying mass and momentum balance across a moving interface Γ. Let Γ divide the material domain *V* into subdomains $V+\u2009and\u2009V\u2212$ and let the outward normal to $V+\u2009on\u2009\Gamma $ be denoted by **n**. The jump condition across Γ derived from the axiom of mass balance requires that

where $u\Gamma \u2261v\u2212v\Gamma $ on $\Gamma \u2009and\u2009v\Gamma $ is the velocity of the interface Γ. Thus, $u\Gamma $ represents the velocity of the fluid relative to Γ. The double bracket notation denotes $[[f]]=f+\u2212f\u2212$, where $f+\u2009and\u2009f\u2212$ represent the value of $f\u2009on\u2009\Gamma \u2009in\u2009V+\u2009and\u2009V\u2212$, respectively. This jump condition implies that the mass flux normal to Γ must be continuous. In particular, if $V+$ is a fluid domain and $V\u2212$ is a solid domain, and Γ denotes the solid boundary (e.g., a wall), we use $\rho +=\rho ,\u2009v+=v$ for the fluid, and $v\u2212=v\Gamma $ for the solid, such that Eq. (2.10) reduces to $\rho (v\u2212v\Gamma )\xb7n=0$. The jump condition derived from the axiom of linear momentum balance similarly requires that

This condition implies that the jump in the traction $\sigma \xb7n$ across Γ must be balanced by the jump in momentum flux normal to Γ. In addition to jump conditions dictated by axioms of conservation, viscous fluids require the satisfaction of the no-slip condition

In our finite element treatment, we use $v\u2009and\u2009J$ as nodal variables, implying that our formulation automatically enforces continuity of these variables across element boundaries, thus $[[v]]=[[u\Gamma ]]=0$ and $[[J]]=0$. Based on Eqs. (2.7) and (2.9), it follows that the density and elastic pressure are continuous across element boundaries in this formulation, $[[\rho ]]=0\u2009and\u2009[[p]]=0$. Thus, the mass jump in Eq. (2.10) is automatically satisfied, and the momentum jump in Eq. (2.11) reduces to $[[\sigma ]]\xb7n=0$, requiring continuity of the traction, or more specifically according to Eq. (2.8), the continuity of the viscous traction $\tau \xb7n$, since *p* is automatically continuous.

The nodal unknowns in this formulation are $v\u2009and\u2009J$ (or *e*), which may be solved using the momentum balance in Eq. (2.1) and the kinematic constraint between $J\u2009and\u2009v$ given in Eq. (2.6). The virtual work integral for a Galerkin finite element formulation [43] is given by

where $\delta v$ is a virtual velocity and $\delta J$ is a virtual energy density, Ω is the fluid finite element domain, and $dv$ is a differential volume in Ω. This virtual work statement may be directly related to the axiom of energy balance, specialized to conditions of isothermal flow of viscous compressible fluids (see Appendix). Using the divergence theorem, we may rewrite the weak form of this integral as the difference between external and internal virtual work, $\delta W=\delta Wext\u2212\delta Wint$, where

and

Here, $\u2202\Omega $ is the boundary of $\Omega \u2009and\u2009da$ is a differential area on $\u2202\Omega ,\u2009t\tau =\tau \xb7n$ is the viscous component of the traction $t$, and $vn=v\xb7n$ is the velocity normal to the boundary $\u2202\Omega $, with $n$ representing the outward normal on $\u2202\Omega $. From these expressions, it becomes evident that essential (Dirichlet) boundary conditions may be prescribed on $v\u2009and\u2009J$, while natural (Neumann) boundary conditions may be prescribed on $t\tau $ and $vn$. The appearance of velocity in both essential and natural boundary conditions may seem surprising at first. To better understand the nature of these boundary conditions, it is convenient to separate the velocity into its normal and tangential components, $v=vnn+vt$, where $vt=(I\u2212n\u2297n)\xb7v$. In particular, for inviscid flow, the viscous stress $\tau $ and its corresponding traction $t\tau $ are both zero, leaving $vn$ as the sole natural boundary condition; similarly, *J* becomes the only essential boundary condition in such flows, since $vt$ is unknown a priori on a frictionless boundary and must be obtained from the solution of the analysis.

In general, prescribing *J* is equivalent to prescribing the elastic fluid pressure, since *p* is only a function of *J*. On a boundary where no conditions are prescribed explicitly, we conclude that $vn=0$ and $t\tau =0$, which represents a frictionless wall. Conversely, it is possible to prescribe $vn\u2009and\u2009t\tau $ on a boundary to produce a desired inflow or outflow while simultaneously stabilizing the flow conditions by prescribing a suitable viscous traction. Prescribing essential boundary conditions $vt$ and *J* determines the tangential velocity on a boundary as well as the elastic fluid pressure *p*, leaving the option to also prescribe the normal component of the viscous traction, $tn\tau =t\tau \xb7n$, to completely determine the normal traction $tn=t\xb7n$ (or else $tn\tau $ naturally equals zero). Mixed boundary conditions represent common physical features: Prescribing $vn\u2009and\u2009vt$ completely determines the velocity $v$ on a boundary; prescribing $t\tau \u2009and\u2009J$ completely determines the traction $t=\sigma \xb7n$ on a boundary. Note that $vn\u2009and\u2009J$ are mutually exclusive boundary conditions, and the same holds for $vt$ and the tangential component of the viscous traction, $tt\tau =(I\u2212n\u2297n)\xb7t\tau $.

The time derivatives, $\u2202v/\u2202t$ which appears in the expression for $a$ in Eq. (2.2), and $\u2202J/\u2202t$ which similarly appears in $J\u02d9$, may be discretized upon the choice of a time integration scheme, such as the generalized-$\alpha $ method [44] (Appendix). The solution of the nonlinear equation $\delta W=0$ is obtained by linearizing this relation as

where the operator $D\delta W[\xb7]$ represents the directional derivative of $\delta W\u2009at\u2009(v,J)$ along an increment $\Delta v\u2009of\u2009v$, or $\Delta J$ of *J* [43]. Using the split form of $\delta W$ between external and internal work contributions, this relation may be expanded as

In this framework, the finite element mesh is defined on the spatial domain Ω, which is fixed (time-invariant) in conventional CFD treatments. Thus, we can linearize $\delta Wint$ along increments $\Delta v$ in the velocity and $\Delta J$ in the volume ratio, by simply bringing the directional derivative operator into the integrals of Eqs. (2.14) and (2.15). The linearization of $\u2202v/\u2202t\u2009and\u2009\u2202J/\u2202t$ is given by

Here, $\Delta t$ is the current time increment and *ξ* is a parameter chosen as described in Ref. [44]; for example, $\xi =1$ yields the backward Euler scheme, in which case $\u2202v/\u2202t\u2248(v\u2212v\u2212\Delta t)/\Delta t$ and $\u2202J/\u2202t\u2248(J\u2212J\u2212\Delta t)/\Delta t$, where $v\u2212\Delta t\u2009and\u2009J\u2212\Delta t$ are the velocity and volume ratio, respectively, at the previous time $t\u2212\Delta t$. More generally, *ξ* is evaluated from the spectral radius for an infinite time-step, $\rho \u221e$ (see Appendix).

The linearization of $\delta Wint$ along an increment $\Delta v$ is then

where we have introduced the fourth-order tensor $C\tau $ representing the tangent of the viscous stress with respect to the rate of deformation

Note that $C\tau $ exhibits minor symmetries because of the symmetries of $\tau \u2009and\u2009D$; in Cartesian components, we have $Cijkl\tau =Cjikl\tau $ and $Cijkl\tau =Cijlk\tau $. In general, $C\tau $ does not exhibit major symmetry ($Cijkl\tau \u2260Cklij\tau $), though the common constitutive relations adopted in fluid mechanics produce such symmetry as shown below.

The linearization of $\delta Wint$ along an increment $\Delta J$ is

where we have used $DJ[\Delta J]=\Delta J$; $p\u2032$ and $p\u2033$, respectively, represent the first and second derivatives of $p(J)$. We have also defined $\tau J\u2032$ as the tangent of the viscous stress $\tau $ with respect to *J*

For the external work, when $t\tau ,\u2009b$, and $vn$ are prescribed, these linearizations simplify to

and

We may define the fluid dilatation $e=J\u22121$ as an alternative essential variable, since initial and boundary conditions $e=0$ are more convenient to handle in a numerical scheme than $J=1$. It follows that $\u2009grad\u2009J=grad\u2009e$ and $\u2202J/\u2202t=\u2202e/\u2202t$. Therefore, the changes to the above-mentioned equations are minimal, simply requiring the substitution $J=1+e\u2009and\u2009\Delta J=\Delta e$. Steady-state analyses may be obtained by setting the terms involving $\Delta t\u22121$ to zero in Eqs. (2.18)–(2.20), and (2.22).

The most common family of constitutive relations employed for viscous fluids, including Newtonian fluids, is given by

where $\mu \u2009and\u2009\kappa $ are, respectively, the dynamic shear and bulk viscosity coefficients (both positive), which may generally depend on *J* and, for non-Newtonian fluids, on invariants of $D$. In practice, most constitutive models for non-Newtonian viscous fluids only use a dependence on $\gamma \u02d9=2D:D$, since it is the only nonzero invariant in viscometric flows [45]. In this case, substituting Eq. (2.26) into Eq. (2.21) produces

The term containing $I\u2297D$ is the only one that does not exhibit major symmetry. In Newtonian fluids, $\mu $ and $\kappa $ are independent of $D$; in incompressible fluids, they are independent of *J* (since $J=1$ remains constant and $\u2009tr\u2009D=0$). Thus, for both of these cases the term containing $I\u2297D$ drops out and $C\tau $ exhibits major symmetry.

Similarly, using Eq. (2.26), the tangent $\tau J\u2032$ in Eq. (2.23) reduces to

Explicit forms for the dependence of $\mu \u2009or\u2009\kappa \u2009on\u2009J$ are not illustrated here, since viscosity generally shows negligible dependence on pressure (thus *J*) over typical ranges of pressures in fluids, hence $\tau J\u2032\u22480$ in most analyses.

Many fluid mechanics textbooks employ Stoke's condition ($\kappa =0$) for the purpose of equating the elastic pressure *p* with the mean normal stress $\u22121/3\u2009tr\u2009\sigma $ [46]; in FEBio, $\kappa $ is kept as a user-defined material property, which may be set to zero if desired. Several constitutive models for non-Newtonian viscous fluids have been implemented in FEBio to date, including the Carreau, Carreau–Yasuda, Powell-Eyring, and Cross models [47]. Their implementation is achieved using a C++ class that provides functions to return $\mu \u2009and\u2009\kappa $ as functions of $(J,\gamma \u02d9)$, and others functions that return $C\tau $, evaluated as shown in Eq. (2.27), and $\tau J\u2032$ as shown in Eq. (2.28). For example, the Carreau model, where $\tau =2\mu (\gamma \u02d9)D$, is a special case of Eq. (2.26), with $\kappa =2\mu /3$ and

where $\lambda $ is a time constant, $n$ is a parameter governing the power-law response, $\mu 0$ is the viscosity when $\gamma \u02d9=0$ and $\mu \u221e$ is the viscosity as $\gamma \u02d9\u2192\u221e$.

For nearly incompressible fluids, a simple constitutive relation may be adopted for the pressure

where *K* is the bulk modulus of the fluid in the limit when $J=1$; it is a physical property that may be measured or looked up in a handbook. It follows that $p\u2032(J)=\u2212K\u2009and\u2009p\u2033(J)=0$ in Eq. (2.22). This constitutive relation is adopted for nearly incompressible CFD analyses in FEBio, though alternative formulations may be easily implemented. For example, for an ideal gas, the equation of state for its absolute pressure is $pa=R\theta \rho /M$, where $R$ is the universal gas constant, $\theta $ is the absolute temperature, and *M* is the gas molar mass. The ambient pressure $pr$ may be obtained from this formula when the gas is at some reference density $\rho r$. Thus, using Eq. (2.7), the gauge pressure $p=pa\u2212pr$ is given by

By definition, the effective bulk modulus of a fluid is $Ke=\u2212J\u2009p\u2032$, and *K* is the value of $Ke$ evaluated at $J=1$. Thus, for an ideal gas we have $K=R\theta \rho r/M$, which is the absolute pressure in the reference state.

The velocity $v(x,t)$ and Jacobian $J(x,t)$ are spatially interpolated over the domain Ω using the same interpolation functions $Na(x)$, with $a=1$ to $n\u2009where\u2009n$ is the number of nodes in an element)

Here, $va\u2009and\u2009Ja$ are nodal values of $v$ and *J* that evolve with time. In contrast to classical mixed formulations for incompressible flow [32], which solve for the pressure $p\u2009using\u2009\u2009div\u2009v=0$ instead of Eq. (2.6), equal-order interpolation is acceptable in this formulation since the governing equations for $v\u2009and\u2009J$ involve spatial derivatives of both variables ($\u2009grad\u2009v\u2009and\u2009\u2009grad\u2009J$). The expressions of Eq. (2.32) may be used to evaluate $L,\u2009\u2009div\u2009v,\u2009a,\u2009\u2009grad\u2009J$, $J\u02d9$, etc. Similar interpolations are used for virtual increments $\delta v\u2009and\u2009\u2009\delta J$, as well as real increments $\Delta v$ and $\Delta J$.

When substituted into Eq. (2.14), we find that the discretized form of $\delta Wint$ may be written as

where

Similarly, the discretized form of $D\delta Wint[\Delta v]$ in Eq. (2.20) becomes

where

whereas that of $D\delta Wint[\Delta J]$ in Eq. (2.22) becomes

where

For the external work in Eq. (2.15), its discretized form is

where

The discretized form of $D(\delta Wext)[\Delta J]$ in Eq. (2.25) is

The nonsymmetric tangent stiffness matrix $[Kint]$ resulting from the linearization of the internal work may be constructed from Eqs. (2.35)–(2.41)

The effective fluid bulk modulus $Ke$ only appears within $kabvJ$, as a linear term, as may be construed from Eqs. (2.38) and (2.23). Furthermore, $kabJJ$ is not zero in general, as may be noted from Eq. (2.38). For this matrix structure, it can be shown that the condition number becomes proportional to the bulk modulus as $Ke\u2192\u221e$. Nevertheless, similar to the argument presented by Ryzhakov et al. [48], the stabilization provided by the nonzero $kabJJ$ submatrix is sufficient to produce a well-behaved solution even for very large values of $Ke$.

In contrast, mixed formulations solve for the velocity $v$ and pressure *p* using the momentum balance, Eq. (2.1), supplemented by the mass balance $\u2009div\u2009v=0$ for incompressible flow. The resulting tangent matrix has a zero entry in lieu of $kabJJ$, and this type of matrix must pass the inf-sup condition to prevent an overconstrained system of equations [32–34]. Stabilization methods such as the SUPG method, which have been introduced to minimize spurious velocity oscillations on coarse meshes [35], populate this matrix entry [2,36], allowing the use of equal-order interpolations for velocity and pressure. The characteristic-based split algorithm presented by Zienkiewicz and Codina similarly produces a nonzero entry for arbitrary velocity and pressure interpolations, when using a suitable dual time-stepping scheme [49,50]. Indeed, these and other investigators have solved for the pressure *p* using the relation $p\u02d9=\u2212Ke\u2009div\u2009v$ [48–50], which may be reproduced from our treatment by evaluating $p\u02d9=p\u2032J\u02d9$ (valid at constant temperature), substituting $J\u02d9$ from Eq. (2.6), and using the definition of $Ke$ above. However, it should be noted that these prior studies [48–50] used this relation with $p\u02d9=\u2202p/\u2202t$, which is only appropriate in a material description, whereas the above-mentioned derivations show that $p\u02d9=\u2202p/\u2202t+grad\u2009p\xb7v$ should be evaluated in a spatial description for consistency.

For arterial blood flow, backflow stabilization has been proposed previously to deal with truncated domains where the entire artery is not modeled explicitly [51,52]; for these types of problems, letting $t=0$ or prescribing a constant pressure at the outflow boundary may not prevent flow reversals that compromise convergence of an analysis. Instead, these authors proposed a velocity-dependent traction boundary condition, $t=\beta \rho (v\u2297v)\xb7n$ with a tensile normal component, that counters the backflow (only when $vn<0$). Here, *β* is a nondimensional user-defined parameter; a value of $\beta =0$ turns off this feature, while a value of *β* = 1 generally shows good numerical performance. We adapt this previously proposed formulation by letting the normal component of the viscous traction be given by

The choice of $\rho r$ in lieu of *ρ* is for convenience, to avoid the dependence of $\rho \u2009on\u2009\u2009J$ (which is negligible for nearly incompressible flow). Then, the contribution of this traction to the virtual external work $\delta Wext$ is

The linearization of *δG* along an increment $\Delta v$ in the velocity is given by

The discretized form of *δG* is

whereas the discretized form of $D\delta G[\Delta v]$ is

A (viscous) tangential traction is implemented as a separate flow stabilization method in Sec. 2.7.2, applicable to inlet or outlet surfaces, without a conditional requirement based on the sign of $vn$.

For certain outlet conditions, using the natural boundary condition $tt\tau =0$ may lead to flow instabilities. It is possible to minimize these effects by prescribing a tangential viscous traction onto the boundary surface, which opposes this tangential flow. Optionally, this condition may be combined with the backflow stabilization described previously.

Similar to Sec. 2.7.1, we introduce a nondimensional parameter *β*, with the tangential traction given by

This form shows that $tt\tau $ opposes tangential flow. The external virtual work for this traction is

Its linearization along an increment $\Delta v$ is

where it can be shown that

The discretized form of *δG* is

The discretized form of $D\delta G[\Delta v]$ is

Flow resistance is typically implemented when modeling arterial flow, where the finite element domain only describes a portion of an arterial network [53]. A flow resistance may be imposed on downstream boundaries to simulate the resistance produced by the vascular network with its branches and bifurcations. The resistance is equivalent to a mean pressure which is proportional to the volumetric flow rate $Q$

where $R$ is the resistance. Using the pressure–dilatation relation Eq. (2.30), equivalent to $p=\u2212K\xb7e$, this pressure may be prescribed as an essential boundary condition on the dilatation *e*.

The numerical solution of the nonlinear system of equations is performed using Newton's method for the first iteration of a time point, followed by Broyden quasi-Newton updates [54], until convergence is achieved at that time-step. The stiffness matrix is thus evaluated only 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 (e.g., 50 updates); this approach offers considerable numerical efficiency as illustrated in the results later. Relative convergence is achieved at each discrete time $tn$ when the vector $\Delta U(k)$ of nodal degree-of-freedom (DOF) increments (consisting of $\Delta v$ or $\Delta J$ values at all unconstrained nodes) at the kth iteration satisfies $\Vert \Delta U(k)\Vert \u2264\epsilon r\Vert \Delta U(1)\Vert $, where $\epsilon r$ represents the relative tolerance criterion ($\epsilon r=10\u22123$ for both $\Delta v\u2009and\u2009\u2009\Delta J$ in the problems analyzed later, unless specified otherwise). 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 which the most efficient is the PARDISO parallel direct sparse solver for the solution of a linear system of equations with nonsymmetric coefficient matrix [55].^{12} Results reported later employ the Intel Math Kernel Library implementation of PARDISO.^{13}

In the sections below, we report the results of standard benchmark problems for computational fluid dynamics, followed by some comparisons of FEBio simulations with Simvascular and the commercial code ANSYS Fluent.^{14} All analyses use a Newtonian fluid with $\kappa =0$; similarly, all analyses use the backward Euler time integration scheme unless specified. Additional analyses are described in the Supplementary Materials section, which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection, including two-dimensional (2D) flow past a cylinder at Re = 100 (Sec. S3 of the Supplementary Materials section, which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection.), exhibiting the characteristic Kármán vortex street and fluctuations in drag and lift coefficients, and flow in a model of an ascending and descending aorta, illustrating flow resistance, backflow stabilization, and tangential stabilization.

This is a benchmark 2D steady flow problem, defined over a square region with boundary conditions described in Fig. 1. In the current implementation, it was necessary to set the dilatation *e* to zero at the top corners, as they represent singularity points where *e* would otherwise blow up in a numerical analysis (since $\u2202v1/\u2202x1\u2192\u221e$ at those points). The Reynolds number for this problem, $Re=\rho VL/\mu $, is based on the length $L$ of the square sides, and the velocity *V* of the lid. FEBio results are shown for $Re=400\u2009and\u2009\u2009Re=5000$ (with $\rho =1$, $V=1,\u2009L=1$), using an unstructured mesh of linear (four-node) tetrahedral elements with mesh refinement near the boundaries (96,318 elements, 26,838 nodes), as well as a biased mesh of linear (eight-node) hexahedral elements (128 × 128 elements, 33,282 nodes). All analyses used a bulk modulus $K=109$. Transient analyses were performed using increasing time increments up to time 10^{3}, at which point a steady response was observed. Representative vorticity contours are shown in Fig. 1. A comparison of steady-state horizontal and vertical velocity profiles across the midsections of the domain is provided against the benchmark numerical solution of Ghia et al. [56] in Fig. 2. The results show equally good velocity agreement for structured hexahedral and unstructured tetrahedral meshes, though contour plots of pressure and vorticity are evidently smoother when using the structured mesh.

This benchmark 2D flow problem, described in Ref. [50], may be compared to experimental measurements of fluid velocity [57], thus serving as a code validation problem. The domain dimensions in the $x1\u2212x2$ plane and mesh are shown in Fig. 3, with a step height $L=1$; the unstructured mesh consisted of 27,008 four-node tetrahedral elements (28,872 nodes) with a thickness 0.1 along $x3$. The fluid properties were $K=109,\u2009\rho r=1$, and $\mu =1$ and the prescribed inlet velocity $v1$ varied along $x2$ according to experimental data, with a mean value *V* producing $Re=\rho rVL/\mu =229$. A steady-state analysis was performed in FEBio, producing velocity and pressure contours as shown in Fig. 3. Velocity profiles were compared to experimental results in Fig. 4, showing very good agreement.

A one-dimensional (1D) analysis of inviscid flow was analyzed to investigate the different time integration schemes. A rectangular domain of width 1, height 0.1 and depth 0.05 was meshed uniformly with 2000 eight-node hexahedral elements along its length ($\u22120.5\u2264x\u22640.5$). The viscosity $\mu $ was set to zero. Assuming SI units, its bulk modulus was evaluated from the properties of air under ambient condition ($\rho r=1.225$, $K=101\u2009325$). The normal velocity was naturally zero ($vn=0$) on all boundaries, except the leftmost face where the velocity was prescribed as $vn=\u22121/2(1\u2212cos\u20092\pi t/T)$ for $0\u2264t\u2264T$, with $T=3\xd710\u22124$. The analysis was run for 1200 uniform time-steps, to a final time $6\xd710\u22123$, using a relative convergence criterion $\epsilon r=10\u22125$. Euler time integration was used in one case; three additional analyses were performed using generalized *α*-integration with $\rho \u221e=0$, $1/2$, and 1. In all cases, a pressure wave was produced at the leftmost face, which propagated over time to the rightmost face, then reflected back (Fig. 5(a)). For Euler integration, the pressure wave markedly decreased in height and increased in width as the wave propagated to the right, then reflected back in the leftward direction; for $\rho \u221e=1$ the pressure wave profile showed much more subtle changes. To quantify numerical damping caused by the integration scheme, the total internal and kinetic energy $E(t)$ at time *t* was evaluated from the results using

where $\rho \psi =J\u22121\Psi r,\u2009\Psi r=K/2(J\u22121)2$, and *V* is the spatial domain volume. A plot of the total energy (Fig. 5(b)) shows that the Euler scheme caused considerable damping in this analysis, with $E(t)$ decreasing by 65% of its peak value, whereas $\rho \u221e=1$ (the midpoint rule) produced zero dissipation (no damping, within six significant digits), as would be expected from theory. Moreover, $\rho \u221e=0$ produced only a small amount of energy dissipation over the duration of this analysis (less than two percent), whereas $\rho \u221e=1/2$ produced less than 0.1% loss. In contrast, Euler's method required only 2978 iterations for the entire analysis (averaging 2.48 iterations for solving the nonlinear system of equations at each time-step); for the remaining analyses, the total number of iterations was $3645\u2009for\u2009\u2009\rho \u221e=0$, $5964\u2009for\u2009\u2009\rho \u221e=1/2$, and $5985\u2009for\u2009\u2009\rho \u221e=1$.

This 2D flow problem examined the effectiveness of backflow and tangential stabilization. A channel of length $L=10$ and height $H=2.5$ was used, with a 1 × 1 block placed at the channel bottom, a distance *D* = 3 from the inflow boundary; the dimensions and mesh (3292 eight-node hexahedral elements and eight 6-node pentahedral elements, with 6880 nodes) are shown in Fig. 6, along with boundary conditions. The fluid properties were $K=109$, $\rho r=1,\u2009and\u2009\u2009\mu =0.0025$. A uniform horizontal velocity *V* was prescribed at the inflow boundary, ramping up linearly from 0 to $1\u2009from\u2009\u2009t=0\u2009to\u2009t=1$, and then maintained constant up to $t=200$, producing $Re=400$ based on the block size. Time increments $\Delta t=0.1$ were used for this transient flow analysis. Both backflow stabilization (*β* = 1) and tangential stabilization (*β* = 1) were prescribed at the outlet boundary. Using Euler integration, significant vortex shedding was observed until $t\u224830$, after which the flow settled into a metastable state until $t\u2248140$, then transitioned toward the steady-state response, achieved at $t\u2248190$. Using generalized *α*-integration, the solution never achieved a steady-state; instead, it settled into a periodic response around $t\u224830$, with a Strouhal number of 0.153. This periodic response was characterized by vortices alternatively shedding from the block and top boundary, producing periodic flow reversals across the top and bottom halves of the outlet (Fig. 6). Without backflow and tangential stabilization, the solution failed to converge beyond $t\u224823$ with Euler integration, and $t\u224811$ for generalized *α*-integration, soon after flow reversal began to develop at the outlet boundary.

A bifurcated carotid artery model was obtained from the GrabCAD online community,^{15} converted to the length unit of meter, imported into SimVascular^{10} and meshed to include a boundary layer refinement, using four-node tetrahedral elements. Four meshes were created as listed in Table 1, with two of these meshes (“coarse” and “finer”) shown in Fig. 7(a). Finite element models with identical meshes, boundary conditions, material properties ($\rho r=1060\u2009kg/m3,\u2009\mu =0.004\u2009Pa\xb7s$, $K=2\xd7109\u2009Pa$), and time increments (250 increments of $\Delta t=2\u2009ms$), were created in FEBio, SimVascular (version 2.0.20624), and ANSYS Fluent (Version 16.2). The inlet has a diameter of 6.28 mm and the two outlets have respective diameters of 4.26 mm and 3.04 mm. An inlet velocity $v=vnn$ was prescribed with a parabolic spatial profile, and an average value $v0(t)$ whose time history is shown in Fig. 7(b), ranging from a minimum of $0.10\u2009\u2009m/s\u2009$ to a maximum of $0.48\u2009\u2009m/s$ ($Re=165\u2212800$). A constant pressure $p0=13.3\u2009kPa$ was prescribed at the outlet boundaries, as well as on the rim of the inlet boundary. SimVascular solutions were obtained using the default solver (svLS-NS) with $\rho \u221e=0.5$, residual criteria of 0.001 (as recommended by the SimVascular Simulation Guide), and step construction consisting of a minimum of three and a maximum of twelve nonlinear iteration sequences. The Fluent solutions used the pressure-based solver and SIMPLE scheme with least squares cell based gradient, second-order pressure, and second-order upwind momentum; the transient formulation was first-order implicit; absolute convergence criteria were used, with values of 10^{−4} for “continuity” and 10^{−5} for velocity components; iterations per time-step were set to a maximum of 20. All models were solved on the same desktop computer, using either a single processor or eight processors. For this problem, no noticeable differences were observed in the FEBio responses when comparing Euler integration with $\rho \u221e=0$ or $\rho \u221e=0.5$, beyond the first 22 ms; similarly, no noticeable differences could be observed in the SimVascular responses with Euler integration, $\rho \u221e=0\u2009or\u2009\rho \u221e=0.5$ beyond the first 14 ms.

Solution times for all models are reported in Table 1. A representative FEBio solution for the coarse mesh is presented in Fig. 8, displaying the wall shear stress (WSS) and flow velocity field at time $t=0.2\u2009\u2009s$ (when WSS has peaked). For each mesh, a comparison of the FEBio, SimVascular and Fluent responses (wall shear stress, pressure and velocity magnitude) was performed at six distinct locations. The pressure and WSS results for two of these points (*P*_{1} and *P*_{2} in Fig. 7(a)) are reported here. An examination of the FEBio results for the fluid pressure *p* at point *P*_{1} shows that all four meshes produce nearly identical responses (Fig. 9(a)). A comparison of FEBio with SimVascular and Fluent, using the finer mesh, shows near identical responses as well (Fig. 9(b)). The value of WSS at *P*_{1} also shows good agreement among the three software programs (Fig. 9(c)). At *P*_{2}, Fluent showed the least variation in the WSS response with increasing mesh density, followed by FEBio, then SimVascular (Fig. 10(a)–10(c)). However, a comparison of the three software programs, using the finer mesh, demonstrated better agreement between FEBio and SimVascular, with Fluent producing lower values of WSS (Fig. 10(d)). The Courant–Friedrichs–Lewy (CFL) number was evaluated for all elements in the coarse and finer meshes, at the instant when the input velocity peaked. The maximum CFL number among all elements was 7.0 for the coarsest mesh and 9.8 for the finer mesh. In both models, the median CFL number for all elements was 1.1.

An experimental study by Seeley and Young explored the pressure loss across an idealized model of arterial stenosis [58]. These authors simulated a stenosis by placing a cylindrical plug of length *L* and outer diameter $D0$, with a hole of diameter $D1$, inside a long tube of inner diameter $D0$. Using water and water–glycerol mixtures, they related the Euler number, $Eu=\Delta p/\rho U02$, to the Reynolds number, $Re=\rho U0D0/\mu $, between two locations, 0.3 m upstream of the plug and 0.61 m downstream of it. In the expression for $Eu,\u2009\Delta p$ is the pressure drop between the pressure ports in the presence of the plug, minus the pressure drop over the same distance in the absence of a plug (Poiseuille flow); $U0$ is the mean flow velocity upstream and downstream of the plug. An additional length of pipe, 1.5 m long, extended upstream of the first pressure port, to ensure fully developed flow at that location. Detailed experimental results for Eu versus Re were provided for the specific case of a concentric hole, where $D0=1.27\u2009cm,\u2009D1=D0/4$, and $L/D0=2$, which are reproduced here (Fig. 11).

This experimental configuration was modeled in FEBio by including the entire length of pipe for which dimensions were provided in that paper [58] (the actual experimental setup had an additional length of pipe and a valve downstream of the second pressure port, for which dimensions were not provided, pouring into a reservoir at ambient pressure). An axisymmetric analysis was simulated in 3D by only modeling a wedge of the circular domain, spanning a 6 deg angle. This type of axisymmetric analysis took advantage of the natural boundary condition $vn=0$ available in our formulation, which allowed frictionless flow along, and prevented flow normal to, the planar wedge faces. The geometry and mesh (linear hexahedral elements) were created using Cubit version 13.2,^{16} with strong mesh biases in the radial direction along the entire length of the pipe, and in the axial direction at the junctures of the plug hole and upstream and downstream pipes (Fig. S9 in Supplementary Materials which are available under “Supplemental Data” tab for this paper on the ASME Digital Collection). The no-slip condition was prescribed on all cylindrical walls and plug upstream and downstream faces. Zero dilatation was prescribed on the outlet (downstream) face, along with backflow and tangential stabilization (*β* = 1 for both). Zero dilatation was prescribed on the rim of the inlet (upstream) face; a uniform outward normal velocity $vn$ was also prescribed on that face, varying smoothly with time from $0\u2009to\u2212U0$ according to

For these analyses we used $U0=1,\u2009t0=0.05,\u2009\rho =1$, and varied $\mu $ to achieve the desired values of Re. The bulk modulus was set to $K=2.2\xd7109$. A total of 500 time-steps were used, with uniform time increments $\Delta t=5\xd710\u22124$, and convergence settings were identical to the carotid bifurcation analysis; Euler integration was employed, as it proved to be the most efficient and reliable scheme, and results are presented here for the steady-state response. (Additional results are available in Sec. S5 of the Supplementary Materials section which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection) Six mesh refinements were examined to verify mesh convergence for the parameter Eu, spanning from 3126 to 33,250 elements (25,440–267,248 degrees-of-freedom), analyzing the case Re = 1000. Results of the mesh convergence analysis are presented in Fig. S10 of Supplementary Materials, which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection, showing that the highest two mesh refinements produced nearly indistinguishable results. All subsequent analyses were performed with the penultimate mesh (17,990 elements, 144,892 degrees-of-freedom), each requiring 15 min of wall clock time or less on a vintage 2011 desktop computer (eight threads) or 7 min on a higher end, vintage 2017 machine (12 threads).

The results of the FEBio analyses (minus the pressure drop evaluated from Poiseuille flow) are plotted with the experimental data of Seeley and Young [58] in Fig. 11, showing good overall agreement with experimental trends and very good accuracy for $Re\u2272400$. With increasing Re, FEBio results tended to moderately overestimate the experimental data, though the lack of measurement uncertainty ranges in the experimental study precluded a definitive assessment of agreement in this higher range of Re. To further examine whether this small overestimation might result from the specific finite element formulation developed in this study, we used the same mesh to perform a 2D axisymmetric analysis in Fluent, with the solver settings described in the carotid bifurcation analysis. Though the prescribed tolerances could not be met exactly, even when allowing up to 4000 iterations, Fluent was nevertheless able to produce solutions similar to those of FEBio (Fig. 11). These findings suggest that the small deviation between CFD and experimental results in the higher range of Re are not peculiar to one type of CFD solver; they may be attributed to measurement or modeling uncertainties, perhaps associated with the length of pipe and valve downstream of the second pressure transducer, which were not modeled in these CFD analyses. More detailed results are presented in Sec. S5 of the Supplementary Materials, which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection.

The objective of this study was to implement a computational fluid dynamics solver in FEBio, with the long-term objective of modeling fluid–structure interactions, including flow alterations driven by active contractile elements in the solid domain, chemical reactions in fluid and solid or multiphasic domains and at their interfaces, as well as growth and remodeling of the solid matrix triggered by signals transduced by fluid flow, relevant to biomechanics, biophysics, and mechanobiology. A novel CFD finite element formulation was developed, using a thermodynamically exact formulation of isothermal compressible flow, which produced a stable numerical scheme without explicitly appealing to stabilization methods.

The code was successfully verified and validated against standard benchmark solutions, including the lid-driven cavity flow (Figs. 1 and 2), flow past a backward-facing step (Figs. 3 and 4), 2D flow past a cylinder (Figs. S3 and S4 in Supplementary Materials which are available under “Supplemental Data” tab for this paper on the ASME Digital Collection), along with comparisons of arterial flow against other open-source and commercial software (Table 1, Figs. 7–10, and Figs. S5 and S6 which are available under “Supplemental Data” tab for this paper on the ASME Digital Collection). We also validated our code against an experimental study of a simulated stenosis [58] (Fig. 11); to the best of our knowledge, this is the first reported CFD investigation of that experimental study. Features employed in vascular flow simulations, such as backflow stabilization [51,52] and flow resistance [53], were adapted to the current formulation. Backflow and tangential stabilization were shown to be effective for stabilizing outlet conditions when strong fluctuations were observed (Fig. 6).

As reviewed earlier, many solution methods have been proposed previously for CFD analyses of nearly incompressible fluids, which typically require special handling to produce stable numerical solutions [2,32,35,36,48–50,59]. In this study, we developed a novel Galerkin-based finite element formulation that uses velocity and dilatation degrees-of-freedom, instead of the velocity–pressure formulations employed traditionally. For this purpose, we specialized the equations of mass, momentum and energy balance to the case of isothermal conditions (Appendix), which allowed us to limit our state variables to $J\u2009and\u2009\u2009D$ (i.e., excluding temperature); in particular, we used the integrated form of the mass balance to solve for the density $\rho =\rho r/J$ in terms of the referential (invariant) density $\rho r$. To solve for $v$ and *J* we employed the linear momentum balance, Eq. (2.1), and the kinematic constraint between $J\u2009and\u2009\u2009v$, Eq. (2.6). This kinematic constraint is typically derived in solid mechanics using a material frame, starting from the definition of **F**. Few investigators appear to avail themselves of the validity of Eq. (2.6) in the spatial frame (i.e., using the spatial descriptions of *J* and $v$ and applying the material time derivative of *J* in the spatial frame). Thus, both governing equations involve spatial gradients of $v\u2009and\u2009\u2009J$, producing a well-conditioned, fully populated stiffness matrix, Eq. (2.43). We used a function of state, Eq. (2.9), to relate the fluid pressure *p* to the volume ratio *J*, which naturally introduced the bulk modulus *K* (a physical property of the fluid) into our formulation, Eq. (2.30). In standard penalty methods [30], the introduction of a similar parameter compromises the stiffness matrix conditioning, the accuracy of the inverted stiffness matrix and thus the convergence of the nonlinear solver [60]. In contrast, varying *K* in our formulation from $106\u2009to\u20091012$ had no effect on the nonlinear convergence behavior, as evidenced by analyses of the convergence rate presented in Sec. S2 of Supplementary Materials which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection.

From a theoretical perspective, it should be noted that our adoption of *J* as a nodal degree-of-freedom implies that we enforce continuity of fluid strain across element boundaries. Since $v$ is also a nodal degree-of-freedom, it follows that the mass flux $\rho v=\rho rJ\u22121v$ is similarly continuous across element boundaries, consistent with mass balance requirements (recall that $\rho r$ is constant). Had we chosen $w=J\u22121v$ as a nodal degree-of-freedom in lieu of $v$, mass balance across element boundaries would also be automatically satisfied. However, the kinematic constraint of Eq. (2.6) would have reduced to $\u2202J/\u2202t=J2\u2009div\u2009w$; since this expression does not contain $\u2009grad\u2009J$, the numerical scheme would not produce the same stability as achieved in the above-mentioned formulation. Therefore, using $v\u2009and\u2009\u2009J$ appears to be the better choice for numerical stability.

The introduction of dilatation as a nodal degree-of-freedom implies that essential or boundary conditions must be applied in relation to this variable. In some cases, the dilatation needed to be prescribed on a corner edge as an essential boundary condition. Had it not been prescribed, it would be naturally assumed that the normal velocity is zero on that edge. However, normal velocity is ambiguous on an edge (e.g., a corner) that separates two surfaces, since the normal is not defined uniquely. By not specifying the dilatation, the burden would be placed on the numerical scheme to resolve this ambiguity. Sometimes, this approach fails to produce good numerical convergence, whereas prescribing the dilatation succeeds, as shown in some of our illustrations. Therefore, users of this velocity–dilatation formulation should keep in mind the occasional necessity to employ this essential boundary condition, rather than relying on the natural boundary condition on a corner edge. This type of decision in the choice of boundary conditions is not unusual in fluid analyses where the real fluid domain is truncated for the purpose of creating a finite element domain with inlet and outlet boundaries. The true boundary conditions at the inlet and outlet are not known; they are approximated for the purpose of the finite element analysis, as also demonstrated in the scheme for stabilizing backflow conditions on a truncated downstream boundary (Sec. 2.7) [51,52].

The fact that the fluid is modeled as compressible implies that dilatational (pressure) waves travel at a finite speed in this material. In principle, for any transient problem being analyzed, the wave propagation may be examined explicitly by using time increments small enough to capture its temporal evolution, as illustrated in the 1D problem in Fig. 5. Conversely, for problems whose salient time scale is much larger than the time required for waves to propagate through the domain, time increments may be chosen to capture those salient phenomena, such as vortex shedding frequency (Fig. 6) or pulsatile velocity frequency (Fig. 9). For such analyses, wave propagation has a negligible effect and does not interfere with the investigation of those lower frequency phenomena. Thus, modeling the fluid as a real compressible material does not produce any complication in the analysis of problems where the effects of its compressibility are negligible.

As a consequence of this velocity–dilatation formulation, our finite element formulation for isothermal compressible flow employs a straightforward implementation of Galerkin's method, which does not require the explicit satisfaction of the inf-sup condition [32–34] or the use of stabilization methods such as SUPG, Galerkin/least-squares, and pressure-stabilizing/Petrov–Galerkin [35–39], or characteristic-based split [49,50]. Unlike methods that employ stabilization schemes, such as SUPG, no special consideration is required for linear or higher order element interpolations (see lid-driven cavity flow solved with quadratic elements in Sec. S1 of Supplementary Materials which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection), nor is it necessary to carefully develop ad hoc formulas for the stabilization parameter *τ* that serves as a weight factor between the Galerkin residual and least-squares terms [39,41,42]. We expect that the resulting simplicity of this formulation makes it more amenable for future extensions, such as the incorporation of chemical reactions as previously done in FEBio's multiphasic domain [9], and FSI.

Based on the variety of problems examined in this study, we observed that FEBio produced significantly more accurate results when mesh refinement was introduced near nonslip boundaries (boundary layer meshing), in comparison with global refinement of uniform meshes. It was also necessary to prescribe the dilatation *e* along edges where a large gradient in velocity might occur, as noted at corners of the lid-driven cavity flow (Fig. 1), or the edges of inlet boundaries where velocity was prescribed, as in the channel flow (Fig. 6) and the carotid bifurcation (Fig. 7). Additionally, FEBio analyses occasionally benefited from ramping up prescribed velocities or pressures from zero to the desired initial value using short time increments, as done in the lid-driven cavity flow, the steady-state analysis for flow past a backward-facing step, the flow past a block in a narrow channel, and the idealized model of arterial stenosis, rather than imposing a step jump in those boundary conditions at the very first time-step, since Broyden updates might fail under these conditions.

The generalized *α*-method [44] was implemented successfully in FEBio as part of this study, demonstrating considerably less energy dissipation than Euler's method as shown in the one-dimensional wave propagation analysis using inviscid flow (Fig. 5). In viscous flows, the natural energy dissipation resulting from the fluid viscosity may sometimes dominate over numerical damping, as suggested in the bifurcated artery analysis where little difference was noted between Euler integration and the generalized *α*-method; in other cases of viscous flows, numerical damping may still play a dominant role, as suggested in the flow past a block in a narrow channel (Fig. 6), where steady flow was eventually achieved when using Euler integration, whereas a periodic vortex shedding response was produced with the generalized *α*-method.

The comparison of FEBio to SimVascular and Fluent in the analysis of the carotid bifurcation also provided useful insights regarding the characteristics of our formulation. We noted that Fluent, which uses the finite volume method, exhibited less sensitivity to mesh refinements than FEBio and SimVascular. In particular, wall shear stresses did not change as much with increasing mesh refinement as with FEBio and SimVascular (Fig. 10). Nevertheless, overall, the agreement between the three programs was remarkably good for this carotid bifurcation, despite the fact that each of them uses a significantly different formulation. The solution times for SimVascular and Fluent were very similar (Table 1), whereas FEBio exhibited shorter solution times for this problem. The similarity between SimVascular and Fluent is consistent with the fact that they both use iterative linear solvers at each time-step. Their solution times could be reduced by lowering the maximum number of iterations, at the expense of not meeting the convergence tolerance over a larger number of time-steps. The significantly faster performance of FEBio in this series of carotid bifurcation analyses can be explained by the use of Broyden's method. Although quasi-Newton methods have been used in CFD codes [54], they have not been widely adopted. For fluid mechanics problems, we perform one direct linear solve at the beginning of the first time-step, followed by Broyden nonlinear iterations until convergence. For all subsequent steps, the approximation to the stiffness matrix inverse from the last (converged) Broyden iteration of the previous step is used to initialize the new Broyden iterations. It is rare for Broyden's method to require more than five nonlinear iterations. If the number of Broyden iterations exceeds a user-defined threshold, another direct linear solve is executed. Thus, only a handful of direct linear solves are needed for almost all fluid mechanics problems. Using this approach, we have solved problems with 3 million DOFs without difficulty. This is not a hard limit—it simply represents the largest problem that we have tried. Importantly, the comparison performed in this study was intended primarily to verify the FEBio solution against other standard fluid solvers available to the user community, not as a benchmark for speed of computation, which can vary significantly based on the choices of time increments and convergence tolerances, and differ based on the choice of direct versus iterative solvers. Furthermore, for users focused on vascular flows, alternative solvers such as SimVascular offer a far greater variety of options relevant to the vascular system than the FEBio implementation presented here, including reduced-order models of the circulation, such as lumped parameter models and impedance boundary conditions [1,52,53].

Nevertheless, as model size increases to the millions of DOFs, direct linear solvers, even when used sparingly, become a significant bottleneck due to memory and computation requirements. Furthermore, a closer examination of the trends in solution times (Table 1) shows that the relative performance boost of FEBio decreased with increasing mesh sizes, consistent with faster growth of the number of computations with the number of DOFs for direct solvers. This observation was further supported by the analysis of a larger model of an ascending and descending aorta, described in Sec. S4 of the Supplementary Materials, which is available under “Supplemental Data” tab for this paper on the ASME Digital Collection, where the SimVascular run completed faster than FEBio.

Although this observation provides a strong motivation to make use of iterative solvers and to develop customized preconditioners for use with our fluid mechanics formulation, iterative linear solvers are generally unreliable for problems in computational solid mechanics. Moderate levels of nearly incompressible behavior, contact, shell elements, etc., all provide stiffness matrices that are challenging to handle with even the most carefully constructed preconditioners. Because of the variety of ways that the global stiffness matrix may be rendered stiff by the physics of solid mechanics problems, it is impossible to develop a “one size fits all” approach to preconditioning. As we implement FSI capabilities in FEBio, we will work to develop an effective preconditioning strategy to allow the use of iterative solvers with our fluid mechanics formulation, but we will not rely on iterative solvers exclusively. We plan to develop a nonlinear solution strategy that combines the use of the Jacobian-free Newton–Krlyov method with the Broyden quasi-Newton method. For problems that are large enough to benefit from iterative linear methods, the nonlinear solution process will begin with Newton–Krylov methods and then automatically switch to quasi-Newton methods when the iterative linear solvers fail to converge. We are confident that this approach will serve the vast majority of our current and future users.

In summary, a novel CFD solver has been implemented in the open-source FEBio software suite which is suitable for the analysis of nearly incompressible flows of the type typically encountered in biomechanics and biophysics. This implementation has been verified and validated using a variety of benchmark and test problems and it appears to perform equally well as other, more standard, formulations. An important characteristic of this formulation is its simplicity, requiring a straightforward application of the Galerkin residual method. This simplicity makes it suitable for future extension to fluid–structure interactions and the incorporation of reactive mechanisms in the fluid and structural domains, as well as at their interfaces, for modeling mechanobiology, growth, and remodeling.

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.

Division of Graduate Education, U.S. National Science Foundation (Grant No. NSF GRFP DGE-16-44869).

National Institute of General Medical Sciences, U.S. National Institutes of Health (Award No. R01GM083925).

The energy balance for a continuum may be written in integral form over a control volume *V* as