Finite Element Framework for Computational Fluid Dynamics in FEBio

Ateshian, Gerard A.; Shim, Jay J.; Maas, Steve A.; Weiss, Jeffrey A.

doi:10.1115/1.4038716

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

Introduction

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,¹ 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² 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³ is geared toward aerospace design, Fluidity⁴ is well suited for geophysical fluid dynamics, and REEF3D⁵ is focused on marine applications. Some CFD codes explore alternative solution methods: Palabos⁶ uses the lattice Boltzmann method, LIGGGHTS⁷ uses a discrete-element method particle simulation, the Gerris Flow Solver⁸ uses an octree finite volume discretization, and COOLFluiD⁹ 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.”¹⁰ This open-source code [1] is based on the original work of Taylor [2]; it uses the flow solver from the PHASTA project¹¹ [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.

Finite Element Implementation

Governing Equations.

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

ρ a = div σ + ρ b

(2.1)

where ρ is the density,

σ

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

a = \dot{v} = \frac{\partial v}{\partial t} + L \cdot v

(2.2)

where

L = grad v

is the spatial velocity gradient. The mass balance equation is

\dot{ρ} + ρ div v = 0

(2.3)

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

\dot{ρ} = \frac{\partial ρ}{\partial t} + grad ρ \cdot v

(2.4)

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

\dot{F} = L \cdot F

(2.5)

Let

J = \det F

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 - 1

⁠. Using the chain rule, J's material time derivative is

\dot{J} = J F^{- T} : \dot{F}

which, when combined with Eq. , produces a kinematic constraint between

\dot{J} and div v

\dot{J} = J div v

(2.6)

Substituting this relation into the mass balance, Eq. , produces

\dot{\bar{ρ J}} = 0

⁠, which may be integrated directly to yield

ρ = ρ_{r} / J

(2.7)

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

The Cauchy stress is given by

σ = - p I + τ

(2.8)

where

I

is the identity tensor,

τ

is the viscous stress, p is the pressure arising from the elastic response

p = - \frac{d Ψ_{r} (J)}{d J}

(2.9)

and $Ψ_{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 $Ψ_{r}$ cannot be a function of the rate of deformation $D = (L + L^{T}) / 2$ ⁠. In contrast, the viscous stress $τ$ is generally a function of $J and D$ ⁠.

Boundary conditions may be derived by satisfying mass and momentum balance across a moving interface Γ. Let Γ divide the material domain V into subdomains

V_{+} and V_{-}

and let the outward normal to

V_{+} on Γ

be denoted by n. The jump condition across Γ derived from the axiom of mass balance requires that

[[ρ u_{Γ}]] \cdot n = 0

(2.10)

where

u_{Γ} \equiv v - v_{Γ}

on

Γ and v_{Γ}

is the velocity of the interface Γ. Thus,

u_{Γ}

represents the velocity of the fluid relative to Γ. The double bracket notation denotes

[[f]] = f_{+} - f_{-}

⁠, where

f_{+} and f_{-}

represent the value of

f on Γ in V_{+} and V_{-}

⁠, respectively. This jump condition implies that the mass flux normal to Γ must be continuous. In particular, if

V_{+}

is a fluid domain and

V_{-}

is a solid domain, and Γ denotes the solid boundary (e.g., a wall), we use

ρ_{+} = ρ, v_{+} = v

for the fluid, and

v_{-} = v_{Γ}

for the solid, such that Eq. reduces to

ρ (v - v_{Γ}) \cdot n = 0

⁠. The jump condition derived from the axiom of linear momentum balance similarly requires that

[[σ - ρ u_{Γ} \otimes u_{Γ}]] \cdot n = 0

(2.11)

This condition implies that the jump in the traction

σ \cdot n

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

(I - n \otimes n) \cdot [[u_{Γ}]] = 0

(2.12)

which implies that the velocity component tangential to Γ is continuous across that interface.

In our finite element treatment, we use $v and J$ as nodal variables, implying that our formulation automatically enforces continuity of these variables across element boundaries, thus $[[v]] = [[u_{Γ}]] = 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, $[[ρ]] = 0 and [[p]] = 0$ ⁠. Thus, the mass jump in Eq. (2.10) is automatically satisfied, and the momentum jump in Eq. (2.11) reduces to $[[σ]] \cdot n = 0$ ⁠, requiring continuity of the traction, or more specifically according to Eq. (2.8), the continuity of the viscous traction $τ \cdot n$ ⁠, since p is automatically continuous.

Virtual Work and Weak Form.

The nodal unknowns in this formulation are

v and J

(or e), which may be solved using the momentum balance in Eq. and the kinematic constraint between

J and v

given in Eq. . The virtual work integral for a Galerkin finite element formulation [43] is given by

\begin{matrix} δ W = \int_{Ω} δ v \cdot (div σ + ρ (b - a)) d v \\ + \int_{Ω} δ J (\frac{\dot{J}}{J} - div v) d v \end{matrix}

(2.13)

where

δ v

is a virtual velocity and

δ J

is a virtual energy density, Ω is the fluid finite element domain, and

d v

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,

δ W = δ W_{ext} - δ W_{int}

⁠, where

\begin{matrix} δ W_{int} = \int_{Ω} τ : grad δ v d v + \int_{Ω} δ v \cdot (grad p + ρ a) d v \\ - \int_{Ω} (δ J \frac{\dot{J}}{J} + grad δ J \cdot v) d v \end{matrix}

(2.14)

and

δ W_{ext} = \int_{\partial Ω} δ v \cdot t^{τ} d a + \int_{Ω} δ v \cdot ρ b d v - \int_{\partial Ω} δ J v_{n} d a

(2.15)

Here, $\partial Ω$ is the boundary of $Ω and d a$ is a differential area on $\partial Ω, t^{τ} = τ \cdot n$ is the viscous component of the traction $t$ ⁠, and $v_{n} = v \cdot n$ is the velocity normal to the boundary $\partial Ω$ ⁠, with $n$ representing the outward normal on $\partial Ω$ ⁠. From these expressions, it becomes evident that essential (Dirichlet) boundary conditions may be prescribed on $v and J$ ⁠, while natural (Neumann) boundary conditions may be prescribed on $t^{τ}$ and $v_{n}$ ⁠. 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 = v_{n} n + v_{t}$ ⁠, where $v_{t} = (I - n \otimes n) \cdot v$ ⁠. In particular, for inviscid flow, the viscous stress $τ$ and its corresponding traction $t^{τ}$ are both zero, leaving $v_{n}$ as the sole natural boundary condition; similarly, J becomes the only essential boundary condition in such flows, since $v_{t}$ 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 $v_{n} = 0$ and $t^{τ} = 0$ ⁠, which represents a frictionless wall. Conversely, it is possible to prescribe $v_{n} and t^{τ}$ 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 $v_{t}$ 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, $t_{n}^{τ} = t^{τ} \cdot n$ ⁠, to completely determine the normal traction $t_{n} = t \cdot n$ (or else $t_{n}^{τ}$ naturally equals zero). Mixed boundary conditions represent common physical features: Prescribing $v_{n} and v_{t}$ completely determines the velocity $v$ on a boundary; prescribing $t^{τ} and J$ completely determines the traction $t = σ \cdot n$ on a boundary. Note that $v_{n} and J$ are mutually exclusive boundary conditions, and the same holds for $v_{t}$ and the tangential component of the viscous traction, $t_{t}^{τ} = (I - n \otimes n) \cdot t^{τ}$ ⁠.

Temporal Discretization and Linearization.

The time derivatives,

\partial v / \partial t

which appears in the expression for

a

in Eq. , and

\partial J / \partial t

which similarly appears in

\dot{J}

⁠, may be discretized upon the choice of a time integration scheme, such as the generalized-

α

method [44] (Appendix). The solution of the nonlinear equation

δ W = 0

is obtained by linearizing this relation as

δ W + D δ W [Δ v] + D δ W [Δ J] \approx 0

(2.16)

where the operator

D δ W [\cdot]

represents the directional derivative of

δ W at (v, J)

along an increment

Δ v of v

⁠, or

Δ J

of J [43]. Using the split form of

δ W

between external and internal work contributions, this relation may be expanded as

\begin{matrix} D δ W_{int} [Δ v] + D δ W_{int} [Δ J] - D δ W_{ext} [Δ v] \\ - D δ W_{ext} [Δ J] \approx δ W_{ext} - δ W_{int} \end{matrix}

(2.17)

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

δ W_{int}

along increments

Δ v

in the velocity and

Δ J

in the volume ratio, by simply bringing the directional derivative operator into the integrals of Eqs. and . The linearization of

\partial v / \partial t and \partial J / \partial t

is given by

D (\frac{\partial v}{d t}) [Δ v] = ξ \frac{Δ v}{Δ t}

(2.18)

D (\frac{\partial J}{\partial t}) [Δ v] = ξ \frac{Δ J}{Δ t}

(2.19)

Here, $Δ t$ is the current time increment and ξ is a parameter chosen as described in Ref. [44]; for example, $ξ = 1$ yields the backward Euler scheme, in which case $\partial v / \partial t \approx (v - v^{- Δ t}) / Δ t$ and $\partial J / \partial t \approx (J - J^{- Δ t}) / Δ t$ ⁠, where $v^{- Δ t} and J^{- Δ t}$ are the velocity and volume ratio, respectively, at the previous time $t - Δ t$ ⁠. More generally, ξ is evaluated from the spectral radius for an infinite time-step, $ρ_{\infty}$ (see Appendix).

The linearization of

δ W_{int}

along an increment

Δ v

is then

\begin{matrix} D (δ W_{int}) [Δ v] = \int_{Ω} grad δ v : C^{τ} : grad Δ v d v \\ + \int_{Ω} δ v \cdot ρ ((\frac{ξ I}{Δ t} + L) \cdot Δ v + grad Δ v \cdot v) d v \\ - \int_{Ω} (\frac{δ J}{J} grad J + grad δ J) \cdot Δ v d v \end{matrix}

(2.20)

where we have introduced the fourth-order tensor

C^{τ}

representing the tangent of the viscous stress with respect to the rate of deformation

C^{τ} = \frac{\partial τ}{\partial D}

(2.21)

Note that $C^{τ}$ exhibits minor symmetries because of the symmetries of $τ and D$ ⁠; in Cartesian components, we have $C_{i j k l}^{τ} = C_{j i k l}^{τ}$ and $C_{i j k l}^{τ} = C_{i j l k}^{τ}$ ⁠. In general, $C^{τ}$ does not exhibit major symmetry (⁠ $C_{i j k l}^{τ} \neq C_{k l i j}^{τ}$ ⁠), though the common constitutive relations adopted in fluid mechanics produce such symmetry as shown below.

The linearization of

δ W_{int}

along an increment

Δ J

is

\begin{matrix} D (δ W_{int}) [Δ J] = \int_{b} Δ J τ_{J}^{'} : grad δ v d v - \int_{Ω} δ v \cdot Δ J \frac{ρ}{J} a d v \\ + \int_{Ω} δ v \cdot (p^{'} grad Δ J + Δ J p^{″} grad J) d v \\ - \int_{Ω} \frac{δ J}{J} ((\frac{ξ}{Δ t} - \frac{\dot{J}}{J}) Δ J + grad Δ J \cdot v) d v \end{matrix}

(2.22)

where we have used

D J [Δ J] = Δ J

⁠;

p^{'}

and

p^{″}

⁠, respectively, represent the first and second derivatives of

p (J)

⁠. We have also defined

τ_{J}^{'}

as the tangent of the viscous stress

τ

with respect to J

τ_{J}^{'} = \frac{\partial τ}{\partial J}

(2.23)

For the external work, when

t^{τ}, b

⁠, and

v_{n}

are prescribed, these linearizations simplify to

D (δ W_{ext}) [Δ v] = 0

(2.24)

and

D (δ W_{ext}) [Δ J] = - \int_{b} δ v \cdot Δ J \frac{ρ}{J} b d v

(2.25)

We may define the fluid dilatation $e = J - 1$ 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 $grad J = grad e$ and $\partial J / \partial t = \partial e / \partial t$ ⁠. Therefore, the changes to the above-mentioned equations are minimal, simply requiring the substitution $J = 1 + e and Δ J = Δ e$ ⁠. Steady-state analyses may be obtained by setting the terms involving $Δ t^{- 1}$ to zero in Eqs. (2.18)–(2.20), and (2.22).

Constitutive Relations.

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

τ (J, D) = (κ - \frac{2}{3} μ) (tr D) I + 2 μ D

(2.26)

where

μ and κ

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

\dot{γ} = \sqrt{2 D : D}

⁠, since it is the only nonzero invariant in viscometric flows [45]. In this case, substituting Eq. into Eq. produces

\begin{matrix} C^{τ} = (κ - \frac{2}{3} μ) I \otimes I + 2 / \dot{γ} (\frac{\partial κ}{\partial \dot{γ}} - \frac{2}{3} \frac{\partial μ}{\partial \dot{γ}}) (tr D) I \otimes D \\ + 2 (2 / \dot{γ} \frac{\partial μ}{\partial \dot{γ}} D \otimes D + μ I \underline{\bar{\otimes}} I) \end{matrix}

(2.27)

The term containing $I \otimes D$ is the only one that does not exhibit major symmetry. In Newtonian fluids, $μ$ and $κ$ are independent of $D$ ⁠; in incompressible fluids, they are independent of J (since $J = 1$ remains constant and $tr D = 0$ ⁠). Thus, for both of these cases the term containing $I \otimes D$ drops out and $C^{τ}$ exhibits major symmetry.

Similarly, using Eq. , the tangent

τ_{J}^{'}

in Eq. reduces to

τ_{J}^{'} = (\frac{\partial κ}{\partial J} - \frac{2}{3} \frac{\partial μ}{\partial J}) (tr D) I + 2 \frac{\partial μ}{\partial J} D

(2.28)

Explicit forms for the dependence of $μ or κ on J$ are not illustrated here, since viscosity generally shows negligible dependence on pressure (thus J) over typical ranges of pressures in fluids, hence $τ_{J}^{'} \approx 0$ in most analyses.

Many fluid mechanics textbooks employ Stoke's condition (⁠

κ = 0

⁠) for the purpose of equating the elastic pressure p with the mean normal stress

- 1 / 3 tr σ

[46]; in FEBio,

κ

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

μ and κ

as functions of

(J, \dot{γ})

⁠, and others functions that return

C^{τ}

⁠, evaluated as shown in Eq. , and

τ_{J}^{'}

as shown in Eq. . For example, the Carreau model, where

τ = 2 μ (\dot{γ}) D

⁠, is a special case of Eq. , with

κ = 2 μ / 3

and

μ = μ_{\infty} + (μ_{0} - μ_{\infty}) {(1 + {(λ \dot{γ})}^{2})}^{(n - 1) / 2}

(2.29)

where $λ$ is a time constant, $n$ is a parameter governing the power-law response, $μ_{0}$ is the viscosity when $\dot{γ} = 0$ and $μ_{\infty}$ is the viscosity as $\dot{γ} \to \infty$ ⁠.

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

p (J) = K (1 - J)

(2.30)

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^{'} (J) = - K and p^{″} (J) = 0

in Eq. . 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

p_{a} = R θ ρ / M

⁠, where

R

is the universal gas constant,

θ

is the absolute temperature, and M is the gas molar mass. The ambient pressure

p_{r}

may be obtained from this formula when the gas is at some reference density

ρ_{r}

⁠. Thus, using Eq. , the gauge pressure

p = p_{a} - p_{r}

is given by

p (J) = \frac{R θ ρ_{r}}{M} (\frac{1}{J} - 1)

(2.31)

By definition, the effective bulk modulus of a fluid is $K_{e} = - J p^{'}$ ⁠, and K is the value of $K_{e}$ evaluated at $J = 1$ ⁠. Thus, for an ideal gas we have $K = R θ ρ_{r} / M$ ⁠, which is the absolute pressure in the reference state.

Spatial Discretization.

The velocity

v (x, t)

and Jacobian

J (x, t)

are spatially interpolated over the domain Ω using the same interpolation functions

N_{a} (x)

⁠, with

a = 1

to

n where n

is the number of nodes in an element)

v (x, t) = \sum_{a = 1}^{n} N_{a} (x) v_{a}, J (x, t) = \sum_{a = 1}^{n} N_{a} (x) J_{a}

(2.32)

Here, $v_{a} and J_{a}$ 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 using div v = 0$ instead of Eq. (2.6), equal-order interpolation is acceptable in this formulation since the governing equations for $v and J$ involve spatial derivatives of both variables (⁠ $grad v and grad J$ ⁠). The expressions of Eq. (2.32) may be used to evaluate $L, div v, a, grad J$ ⁠, $\dot{J}$ ⁠, etc. Similar interpolations are used for virtual increments $δ v and δ J$ ⁠, as well as real increments $Δ v$ and $Δ J$ ⁠.

When substituted into Eq. , we find that the discretized form of

δ W_{int}

may be written as

δ W_{int} = \sum_{a} δ v_{a} \cdot (f_{a}^{σ} + f_{a}^{ρ}) + f_{a}^{J} δ J_{a}

(2.33)

where

\begin{matrix} f_{a}^{σ} = \int_{Ω} (τ \cdot grad N_{a} + N_{a} grad p) d v \\ f_{a}^{ρ} = \int_{Ω} N_{a} ρ a d v \\ f_{a}^{J} = \int_{Ω} - (N_{a} \frac{\dot{J}}{J} + grad N_{a} \cdot v) d v \end{matrix}

(2.34)

Similarly, the discretized form of

D δ W_{int} [Δ v]

in Eq. becomes

\begin{matrix} D (δ W_{int}) [Δ v] = \sum_{a} δ v_{a} \cdot \sum_{b} (K_{a b}^{v v} + M_{a b}^{v v}) \cdot Δ v_{b} \\ + \sum_{a} δ J_{a} \sum_{b} k_{a b}^{J v} \cdot Δ v_{b} \end{matrix}

(2.35)

where

\begin{matrix} K_{a b}^{v v} = \int_{Ω} grad N_{a} \cdot C^{v} \cdot grad N_{b} d v \\ M_{a b}^{v v} = \int_{Ω} N_{a} ρ (N_{b} (\frac{ξ}{Δ t} I + grad v) + (grad N_{b} \cdot v) I) d v \\ k_{a b}^{J v} = \int_{Ω} - (\frac{N_{a}}{J} grad J + grad N_{a}) N_{b} d v \end{matrix}

(2.36)

whereas that of

D δ W_{int} [Δ J]

in Eq. becomes

\begin{matrix} D (δ W_{int}) [Δ J] = \sum_{a} δ v_{a} \cdot \sum_{b} (k_{a b}^{v J} + m_{a b}^{v J}) Δ J_{b} \\ + \sum_{a} δ J_{a} \sum_{b} k_{a b}^{J J} Δ J_{b} \end{matrix}

(2.37)

where

\begin{matrix} k_{a b}^{v J} = \int_{Ω} [N_{b} {τ^{'}}_{J} \cdot grad N_{a} + N_{a} (p^{'} grad N_{b} + N_{b} p^{″} grad J)] d v \\ m_{a b}^{v J} = \int_{Ω} - N_{a} N_{b} \frac{ρ}{J} a d v \\ k_{a b}^{J J} = \int_{Ω} - \frac{N_{a}}{J} ((\frac{ξ}{Δ t} - \frac{\dot{J}}{J}) N_{b} + grad N_{b} \cdot v) d v \end{matrix}

(2.38)

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

δ W_{ext} = \sum_{a} δ v_{a} \cdot (f_{a}^{t} + f_{a}^{b}) + δ J_{a} f_{a}^{v}

(2.39)

where

\begin{matrix} f_{a}^{t} & = \int_{\partial Ω} N_{a} t^{τ} d a \\ f_{a}^{b} & = \int_{Ω} N_{a} ρ b d v \\ f_{a}^{v} & = \int_{Ω} - N_{a} v_{n} d a \end{matrix}

(2.40)

The discretized form of

D (δ W_{ext}) [Δ J]

in Eq. is

D (δ W_{ext}) [Δ J] = \sum_{a} δ v_{a} \cdot \sum_{b} k_{a b}^{b} Δ J_{b}

(2.41)

where

k_{a b}^{b} = \int_{Ω} - N_{a} N_{b} \frac{ρ}{J} b d v

(2.42)

Condition Number.

The nonsymmetric tangent stiffness matrix

[K_{int}]

resulting from the linearization of the internal work may be constructed from Eqs. –

[K_{int}] = [\begin{matrix} K_{a b}^{v v} + M_{a b}^{v v} & k_{a b}^{v J} + m_{a b}^{v J} \\ k_{a b}^{J v} & k_{a b}^{J J} \end{matrix}]

(2.43)

The effective fluid bulk modulus $K_{e}$ only appears within $k_{a b}^{v J}$ ⁠, as a linear term, as may be construed from Eqs. (2.38) and (2.23). Furthermore, $k_{a b}^{J J}$ 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 $K_{e} \to \infty$ ⁠. Nevertheless, similar to the argument presented by Ryzhakov et al. [48], the stabilization provided by the nonzero $k_{a b}^{J J}$ submatrix is sufficient to produce a well-behaved solution even for very large values of $K_{e}$ ⁠.

In contrast, mixed formulations solve for the velocity $v$ and pressure p using the momentum balance, Eq. (2.1), supplemented by the mass balance $div v = 0$ for incompressible flow. The resulting tangent matrix has a zero entry in lieu of $k_{a b}^{J J}$ ⁠, 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 $\dot{p} = - K_{e} div v$ [48–50], which may be reproduced from our treatment by evaluating $\dot{p} = p^{'} \dot{J}$ (valid at constant temperature), substituting $\dot{J}$ from Eq. (2.6), and using the definition of $K_{e}$ above. However, it should be noted that these prior studies [48–50] used this relation with $\dot{p} = \partial p / \partial t$ ⁠, which is only appropriate in a material description, whereas the above-mentioned derivations show that $\dot{p} = \partial p / \partial t + grad p \cdot v$ should be evaluated in a spatial description for consistency.

Special Boundary Conditions

Backflow Stabilization.

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 = β ρ (v \otimes v) \cdot n

with a tensile normal component, that counters the backflow (only when

v_{n} < 0

⁠). Here, β is a nondimensional user-defined parameter; a value of

β = 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

t_{n}^{τ} = {\begin{matrix} β ρ_{r} v_{n}^{2} & v_{n} < 0 \\ 0 & v_{n} \geq 0 \end{matrix}

(2.44)

The choice of

ρ_{r}

in lieu of ρ is for convenience, to avoid the dependence of

ρ on J

(which is negligible for nearly incompressible flow). Then, the contribution of this traction to the virtual external work

δ W_{ext}

is

δ G = \int_{\partial Ω} δ v \cdot t_{n}^{τ} n d a

(2.45)

The linearization of δG along an increment

Δ v

in the velocity is given by

D δ G [Δ v] = \int_{\partial Ω} δ v \cdot K_{n} \cdot Δ v d a

(2.46)

where

K_{n} = {\begin{matrix} 2 β ρ_{r} v_{n} (n \otimes n) & v_{n} < 0 \\ 0 & v_{n} \geq 0 \end{matrix}

(2.47)

The discretized form of δG is

δ G = \sum_{a} δ v_{a} \cdot f_{a}^{n}, f_{a}^{n} = \int_{\partial Ω} N_{a} t_{n}^{τ} n d a

(2.48)

whereas the discretized form of

D δ G [Δ v]

is

D δ G [Δ v] = \sum_{a} δ v_{a} \cdot \sum_{b} K_{a b}^{n} \cdot Δ v_{b}, K_{a b}^{n} = \int_{\partial Ω} N_{a} N_{b} K_{n} d a

(2.49)

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 $v_{n}$ ⁠.

Tangential Flow Stabilization.

For certain outlet conditions, using the natural boundary condition $t_{t}^{τ} = 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

t_{t}^{τ} = - β ρ_{r} | v_{t} | v_{t}

(2.50)

This form shows that

t_{t}^{τ}

opposes tangential flow. The external virtual work for this traction is

δ G = \int_{\partial Ω} δ v \cdot t_{t}^{τ} d a

(2.51)

Its linearization along an increment

Δ v

is

D δ G [Δ v] = \int_{\partial Ω} δ v \cdot K_{t} \cdot Δ v d a

(2.52)

where it can be shown that

K_{t} = - β ρ_{r} | v_{t} | (I - n \otimes n + \frac{v_{t}}{| v_{t} |} \otimes \frac{v_{t}}{| v_{t} |})

(2.53)

The discretized form of δG is

δ G = \sum_{a} δ v_{a} \cdot f_{a}^{τ}, f_{a}^{τ} = \int_{\partial Ω} N_{a} t_{t}^{τ} d a

(2.54)

The discretized form of

D δ G [Δ v]

is

D δ G [Δ v] = \sum_{a} δ v_{a} \cdot \sum_{b} K_{a b}^{t} \cdot Δ v_{b}, K_{a b}^{t} = \int_{\partial Ω} N_{a} N_{b} K_{t} d a

(2.55)

Flow Resistance.

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

p = R Q, Q = \int_{\partial Ω} v_{n} d a

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

Solver.

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 $t_{n}$ when the vector $Δ U^{(k)}$ of nodal degree-of-freedom (DOF) increments (consisting of $Δ v$ or $Δ J$ values at all unconstrained nodes) at the kth iteration satisfies $‖ Δ U^{(k)} ‖ \leq ε_{r} ‖ Δ U^{(1)} ‖$ ⁠, where $ε_{r}$ represents the relative tolerance criterion (⁠ $ε_{r} = 10^{- 3}$ for both $Δ v and Δ 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].¹² Results reported later employ the Intel Math Kernel Library implementation of PARDISO.¹³

Verifications and Validations

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.¹⁴ All analyses use a Newtonian fluid with $κ = 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.

Lid-Driven Cavity Flow.

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 $\partial v_{1} / \partial x_{1} \to \infty$ at those points). The Reynolds number for this problem, $Re = ρ V L / μ$ ⁠, is based on the length $L$ of the square sides, and the velocity V of the lid. FEBio results are shown for $Re = 400 and Re = 5000$ (with $ρ = 1$ ⁠, $V = 1, L = 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 = 10^{9}$ ⁠. Transient analyses were performed using increasing time increments up to time 10³, 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.

Fig. 1

View large Download slide

Lid-driven cavity flow: Unstructured tetrahedral mesh (left column) and structured hexahedral mesh (right column). Vorticity contours are shown for $Re = 400$ (middle row) and $Re = 5000$ (bottom row).

Fig. 2

Velocity profiles in lid-driven cavity flow, for hexahedral and tetrahedral meshes, compared to Ghia et al. [56]: (a) v1 versus x2 along vertical centerline (x1=0.5); and (b) v2 versus x1 along horizontal centerline (x2=0.5)

View large Download slide

Velocity profiles in lid-driven cavity flow, for hexahedral and tetrahedral meshes, compared to Ghia et al. [56]: (a) $v_{1}$ versus x₂ along vertical centerline (⁠ $x_{1} = 0.5$ ⁠); and (b) $v_{2}$ versus x₁ along horizontal centerline (⁠ $x_{2} = 0.5$ ⁠)

Flow Past a Backward Facing Step.

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 $x_{1} - x_{2}$ 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 $x_{3}$ ⁠. The fluid properties were $K = 10^{9}, ρ_{r} = 1$ ⁠, and $μ = 1$ and the prescribed inlet velocity $v_{1}$ varied along $x_{2}$ according to experimental data, with a mean value V producing $Re = ρ_{r} V L / μ = 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.

Fig. 3

View large Download slide

Flow past backward facing step with $Re = 229$ ⁠. Top panel shows the geometry, tetrahedral mesh, and boundary conditions. Middle panel shows contours of horizontal velocity component $v_{1}$ ⁠; bottom panel shows pressure contours.

Fig. 4

View large Download slide

Velocity profiles $v_{1} (x_{1}, x_{2})$ in the flow past a backward facing step. Solid curves represent the finite element solution and symbols represent experimental data [57].

Inviscid Flow One-Dimensional Wave Propagation.

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 (⁠

- 0.5 \leq x \leq 0.5

⁠). The viscosity

μ

was set to zero. Assuming SI units, its bulk modulus was evaluated from the properties of air under ambient condition (⁠

ρ_{r} = 1.225

⁠,

K = 101 325

⁠). The normal velocity was naturally zero (⁠

v_{n} = 0

⁠) on all boundaries, except the leftmost face where the velocity was prescribed as

v_{n} = - 1 / 2 (1 - cos 2 π t / T)

for

0 \leq t \leq T

⁠, with

T = 3 \times 10^{- 4}

⁠. The analysis was run for 1200 uniform time-steps, to a final time

6 \times 10^{- 3}

⁠, using a relative convergence criterion

ε_{r} = 10^{- 5}

⁠. Euler time integration was used in one case; three additional analyses were performed using generalized α-integration with

ρ_{\infty} = 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

ρ_{\infty} = 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

E = \int_{V} ρ (ψ + \frac{1}{2} v \cdot v) d V

(3.1)

where $ρ ψ = J^{- 1} Ψ_{r}, Ψ_{r} = K / 2 {(J - 1)}^{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 $ρ_{\infty} = 1$ (the midpoint rule) produced zero dissipation (no damping, within six significant digits), as would be expected from theory. Moreover, $ρ_{\infty} = 0$ produced only a small amount of energy dissipation over the duration of this analysis (less than two percent), whereas $ρ_{\infty} = 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 for ρ_{\infty} = 0$ ⁠, $5964 for ρ_{\infty} = 1 / 2$ ⁠, and $5985 for ρ_{\infty} = 1$ ⁠.

Fig. 5

View large Download slide

1D wave propagation in an inviscid fluid: (a) pressure profile $p (x, t)$ as wave propagates rightward starting from $x = - 0.5$ ⁠, reflects at right wall (⁠ $x = 0.5$ ⁠), and reverses direction; finite element result using Euler integration is shown in orange, generalized α-method with $ρ_{\infty} = 1$ is shown in blue and (b) fluid internal + kinetic energy $E (t)$ ⁠, evaluated over entire domain $- 0.5 \leq x \leq 0.5$ ⁠, as it varies with time t, using Euler method and generalized α-method with various values of $ρ_{\infty}$ ⁠.

Flow Past Block in a Narrow Channel.

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 = 10^{9}$ ⁠, $ρ_{r} = 1, and μ = 0.0025$ ⁠. A uniform horizontal velocity V was prescribed at the inflow boundary, ramping up linearly from 0 to $1 from t = 0 to t = 1$ ⁠, and then maintained constant up to $t = 200$ ⁠, producing $Re = 400$ based on the block size. Time increments $Δ 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 \approx 30$ ⁠, after which the flow settled into a metastable state until $t \approx 140$ ⁠, then transitioned toward the steady-state response, achieved at $t \approx 190$ ⁠. Using generalized α-integration, the solution never achieved a steady-state; instead, it settled into a periodic response around $t \approx 30$ ⁠, 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 \approx 23$ with Euler integration, and $t \approx 11$ for generalized α-integration, soon after flow reversal began to develop at the outlet boundary.

Fig. 6

View large Download slide

Flow past block in a narrow channel. Top panel shows geometry, mesh, and boundary conditions. Remaining panels show vorticity contour plots and velocity vector plots (⁠ $ρ_{\infty} = 1$ ⁠) at four time points within a cycle of the periodic response achieved after $t \approx 30$ ⁠.

Carotid Bifurcation.

A bifurcated carotid artery model was obtained from the GrabCAD online community,¹⁵ converted to the length unit of meter, imported into SimVascular¹⁰ 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 (⁠ $ρ_{r} = 1060 kg / m^{3}, μ = 0.004 Pa \cdot s$ ⁠, $K = 2 \times 10^{9} Pa$ ⁠), and time increments (250 increments of $Δ t = 2 ms$ ⁠), 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 = v_{n} n$ was prescribed with a parabolic spatial profile, and an average value $v_{0} (t)$ whose time history is shown in Fig. 7(b), ranging from a minimum of $0.10 m / s$ to a maximum of $0.48 m / s$ (⁠ $Re = 165 - 800$ ⁠). A constant pressure $p_{0} = 13.3 kPa$ 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 $ρ_{\infty} = 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⁻⁴ for “continuity” and 10⁻⁵ 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 $ρ_{\infty} = 0$ or $ρ_{\infty} = 0.5$ ⁠, beyond the first 22 ms; similarly, no noticeable differences could be observed in the SimVascular responses with Euler integration, $ρ_{\infty} = 0 or ρ_{\infty} = 0.5$ beyond the first 14 ms.

Fig. 7

View large Download slide

(a) Model of bifurcated carotid artery using coarse mesh (top) and finer mesh (bottom) of four-node tetrahedral elements. Boundary conditions are shown on the coarse mesh; in addition, $v = 0$ on the arterial wall. The prescribed inlet velocity has a parabolic profile, with average value v₀ and (b) time history of the average inlet velocity v₀

Table 1

Bifurcated carotid artery mesh parameters. All elements are four-node tetrahedra. Solution times are based on single-processor and eight-processor computations for all three software programs. SimVascular times include the solution time + time to convert results to viewable plot files, whereas FEBio and Fluent automatically produce viewable plot/data files.

			Solution time, 1 thread (min)			Solution time, 8 threads (min)
Mesh	Elements	Nodes	FEBio	SimVascular	Fluent	FEBio	SimVascular	Fluent
Coarse	120,383	21,671	16	31 + 3 = 34	31	8	15 + 3 = 18	19
Medium	233,899	41,218	38	62 + 7 = 69	73	19	28 + 6 = 34	35
Fine	396,341	68,843	78	115 + 12 = 127	114	41	49 + 10 = 59	54
Finer	498,196	86,124	108	142 + 19 = 161	153	56	62 + 13 = 75	68

			Solution time, 1 thread (min)			Solution time, 8 threads (min)
Mesh	Elements	Nodes	FEBio	SimVascular	Fluent	FEBio	SimVascular	Fluent
Coarse	120,383	21,671	16	31 + 3 = 34	31	8	15 + 3 = 18	19
Medium	233,899	41,218	38	62 + 7 = 69	73	19	28 + 6 = 34	35
Fine	396,341	68,843	78	115 + 12 = 127	114	41	49 + 10 = 59	54
Finer	498,196	86,124	108	142 + 19 = 161	153	56	62 + 13 = 75	68

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 s$ (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₁ and P₂ in Fig. 7(a)) are reported here. An examination of the FEBio results for the fluid pressure p at point P₁ 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₁ also shows good agreement among the three software programs (Fig. 9(c)). At P₂, 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.

Fig. 8

View large Download slide

FEBio solution for bifurcated carotid artery model (coarse mesh) at time t = 0.2 s, showing the WSS distribution (top) and the velocity vector field in a sectioned transparent model (bottom); the largest velocity magnitude at this time point was 1.05 m/s.

Fig. 9

Bifurcated carotid artery results at point P1: (a) FEBio mesh convergence analysis for fluid pressure p. Comparison of SimVascular, Fluent and FEBio results for (b) p and (c) WSS, using finite element model with finer mesh

View large Download slide

Bifurcated carotid artery results at point P₁: (a) FEBio mesh convergence analysis for fluid pressure p. Comparison of SimVascular, Fluent and FEBio results for (b) p and (c) WSS, using finite element model with finer mesh

Fig. 10

Bifurcated carotid artery results at point P2: Mesh convergence analyses for WSS with (a) SimVascular, (b) Fluent, (c) FEBio, and (d) comparison of WSS results using finite element model with finer mesh

View large Download slide

Bifurcated carotid artery results at point P₂: Mesh convergence analyses for WSS with (a) SimVascular, (b) Fluent, (c) FEBio, and (d) comparison of WSS results using finite element model with finer mesh

Idealized Model of Arterial Stenosis.

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 $D_{0}$ ⁠, with a hole of diameter $D_{1}$ ⁠, inside a long tube of inner diameter $D_{0}$ ⁠. Using water and water–glycerol mixtures, they related the Euler number, $Eu = Δ p / ρ U_{0}^{2}$ ⁠, to the Reynolds number, $Re = ρ U_{0} D_{0} / μ$ ⁠, between two locations, 0.3 m upstream of the plug and 0.61 m downstream of it. In the expression for $Eu, Δ 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); $U_{0}$ 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 $D_{0} = 1.27 cm, D_{1} = D_{0} / 4$ ⁠, and $L / D_{0} = 2$ ⁠, which are reproduced here (Fig. 11).

Fig. 11

View large Download slide

Experimental and CFD results for simulated stenosis, plotting the Euler number against the Reynolds number: Experimental measurements using water and water + glycerol mixtures are from [58]. CFD results from this study include FEBio analyses at 18 increasing values of Re (from 25 to 1100), and Fluent analyses at a subset of 15 values of Re.

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

v_{n} = 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,¹⁶ 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

v_{n}

was also prescribed on that face, varying smoothly with time from

0 to - U_{0}

according to

v_{n} = {\begin{matrix} - \frac{U_{0}}{2} (1 - cos (π t / t_{0})) & 0 \leq t \leq t_{0} \\ - U_{0} & t > t_{0} \end{matrix}

For these analyses we used $U_{0} = 1, t_{0} = 0.05, ρ = 1$ ⁠, and varied $μ$ to achieve the desired values of Re. The bulk modulus was set to $K = 2.2 \times 10^{9}$ ⁠. A total of 500 time-steps were used, with uniform time increments $Δ t = 5 \times 10^{- 4}$ ⁠, 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 ≲ 400$ ⁠. 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.

Discussion

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 and D$ (i.e., excluding temperature); in particular, we used the integrated form of the mass balance to solve for the density $ρ = ρ_{r} / J$ in terms of the referential (invariant) density $ρ_{r}$ ⁠. To solve for $v$ and J we employed the linear momentum balance, Eq. (2.1), and the kinematic constraint between $J and v$ ⁠, 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 and J$ ⁠, 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 $10^{6} to 10^{12}$ 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 $ρ v = ρ_{r} J^{- 1} v$ is similarly continuous across element boundaries, consistent with mass balance requirements (recall that $ρ_{r}$ is constant). Had we chosen $w = J^{- 1} v$ 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 $\partial J / \partial t = J^{2} div w$ ⁠; since this expression does not contain $grad J$ ⁠, the numerical scheme would not produce the same stability as achieved in the above-mentioned formulation. Therefore, using $v and J$ 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.

Disclaimer

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).
National Institute of General Medical Sciences, U.S. National Institutes of Health (Award No. R01GM083925).

Appendix

Energy Balance

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

\begin{matrix} \frac{d}{d t} \int_{V} ρ (ε + \frac{1}{2} v \cdot v) d V = - \int_{S} ρ (ε + \frac{1}{2} v \cdot v) (v \cdot n) d S + \int_{S} t \cdot v d S + \int_{V} ρ b \cdot v d V \\ - \int_{S} q \cdot n d S + \int_{V} ρ r d V \end{matrix}

(A1)

where

S

is the control surface bounding

V, ε

is the specific internal energy,

q

is the heat flux across

S

⁠, and

r

is the heat supply per mass to the material in V resulting from other sources. Bringing the time derivative inside the integral on the left-hand side, and using the divergence theorem, this integral statement of the energy balance may be written as

\begin{matrix} \int_{V} [ρ (\dot{ε} + v \cdot a) + ρ (ε + \frac{1}{2} v \cdot v) (div v - \frac{\dot{J}}{J})] d V \\ = \int_{V} [σ : D - div q + ρ r + v \cdot (div σ + ρ b)] d V \end{matrix}

(A2)

This statement must be valid for arbitrary control volumes and arbitrary processes, from which we conventionally derive the differential form of the axioms of mass, momentum and energy balance.

For the specialized conditions of a viscous fluid at constant temperature assumed in our treatment, the only state variables for the functions of state $ε, σ, and q$ are $J and D$ (i.e., the temperature is not a state variable since it is assumed constant). Under these conditions, the entropy inequality shows that the specific entropy $η$ and the heat flux $q$ must be zero, and the Cauchy stress $σ$ must have the form of Eq. (2.8) where p is given by Eq. (2.9) as a function of J only, leaving the residual dissipation statement $τ : D \geq 0$ as a constraint that must be satisfied by constitutive relations for $τ$ ⁠. (For a Newtonian fluid, this constraint is satisfied when the viscosities $μ$ and $κ$ are positive.) From these thermodynamic restrictions, we conclude that $ε = ψ$ ⁠, where $ψ$ is the specific (Helmholtz) free energy, with $Ψ_{r} = ρ_{r} ψ$ ⁠.

For the conditions adopted here (isothermal viscous fluid), the axiom of energy balance reduces to

ρ \dot{ψ} = σ : D + ρ r

⁠; since

ψ

is only a function of J, this expression may be further simplified using Eqs. – to produce

τ : D + ρ r = 0

⁠. In other words, isothermal conditions may be maintained only if heat dissipated by the viscous stress is emitted in the form of a heat supply density

ρ r = - τ : D

(heat leaving the system). Now, the integral form of the energy balance in Eq. simplifies to

\int_{V} [v \cdot (div σ + ρ (b - a)) + ρ (ψ + \frac{1}{2} v \cdot v) (\frac{\dot{J}}{J} - div v)] d V = 0

(A3)

A comparison of this statement with the statement of virtual work in Eq. (2.13) establishes a clear correspondence between $δ v and v$ ⁠, and between $δ J$ and $ρ (ψ + 1 / 2 v \cdot v)$ ⁠, with the latter representing the sum of the internal (free) and kinetic energy densities. We now have a clearer understanding of the significance of virtual energy density $δ J$ in the context of this formulation.

Generalized α-Method

For the generalized α-method, we combine the degrees-of-freedom into

Y_{n} = {v, J}_{n}

⁠, where the subscript

n

denotes time

t_{n}

⁠; similarly, we let

{\dot{Y}}_{n} = {\partial v / \partial t, \partial J / \partial t}_{n}

⁠. According to this method [44], the virtual work is evaluated at

δ W ({\dot{Y}}_{n + α_{m}}, Y_{n + α_{f}})

⁠, where

t_{n + α} = t_{n} + α Δ t and Δ t = t_{n + 1} - t_{n}

⁠. The parameters

α_{f} and α_{m}

are evaluated from a single parameter

ρ_{\infty}

using

α_{f} = \frac{1}{1 + ρ_{\infty}}, α_{m} = \frac{1}{2} \frac{3 - ρ_{\infty}}{1 + ρ_{\infty}}

(A4)

where $0 \leq ρ_{\infty} \leq 1$ ⁠. This parameter is the spectral radius for an infinite time-step, which controls the amount of damping of high frequencies; a value of zero produces the greatest amount of damping, annihilating the highest frequency in one step, whereas a value of one preserves the highest frequency.

The linearization of

δ W

reported previously is effectively performed along

Δ Y_{n + α_{f}}

so that the solution to

δ W = 0

produces

Y_{n + α_{f}}

⁠. The parameter ξ used in Eqs. and is given by

ξ = \frac{α_{m}}{γ α_{f}}, γ = \frac{1}{2} + α_{m} - α_{f}

(A5)

In this scheme,

{\dot{Y}}_{n + α_{m}}

is evaluated from

{\dot{Y}}_{n + α_{m}} = (1 - \frac{α_{m}}{γ}) {\dot{Y}}_{n} + \frac{ξ}{Δ t} (Y_{n + α_{f}} - Y_{n})

(A6)

Given the solution

({\dot{Y}}_{n + α_{m}}, Y_{n + α_{f}})

⁠, the solution at

t_{n + 1}

is evaluated from

\begin{matrix} Y_{n + 1} & = Y_{n} + \frac{Y_{n + α_{f}} - Y_{n}}{α_{f}} \\ {\dot{Y}}_{n + 1} & = {\dot{Y}}_{n} + \frac{{\dot{Y}}_{n + α_{m}} - {\dot{Y}}_{n}}{α_{m}} \end{matrix}

(A7)

Four different options are presented in Ref. [44] for initializing $Y_{n + 1} and {\dot{Y}}_{n + 1}$ at the beginning of time-step $t_{n + 1}$ ⁠; the first three of these have been implemented in FEBio. For steady flows, these authors recommend disregarding $ρ_{\infty}$ and setting $α_{f} = α_{m} = γ = 1$ to recover the backward Euler scheme.

1

febio.org

2

openfoam.org

3

su2.stanford.edu

4

fluidityproject.github.io

5

reef3d.wordpress.com

6

freecode.com/projects/xflows-2

7

www.cfdem.com

8

gfs.sourceforge.net/wiki/

9

coolfluid.github.io

10

simvascular.org

11

secure.cci.rpi.edu/wiki/index.php/PHASTA

12

www.pardiso-project.org

13

software.intel.com/en-us/intel-mkl/

14

www.ansys.com

15

grabcad.com/library/carotid-bifurcation

16

cubit.sandia.gov

References

1.

Updegrove

,

A.

,

Wilson

,

N. M.

,

Merkow

,

J.

,

Lan

,

H.

,

Marsden

,

A. L.

, and

Shadden

,

S. C.

,

2017

, “

SimVascular: An Open Source Pipeline for Cardiovascular Simulation

,”

Ann. Biomed. Eng.

,

45

(

3

), pp.

525

–

541

.

Google Scholar

Crossref

PubMed

2.

Taylor

,

C. A.

,

Hughes

,

T. J.

, and

Zarins

,

C. K.

,

1998

, “

Finite Element Modeling of Blood Flow in Arteries

,”

Comput. Methods Appl. Mech. Eng.

,

158

(

1–2

), pp.

155

–

196

.

Google Scholar

Crossref

3.

Whiting

,

C. H.

, and

Jansen

,

K. E.

,

2001

, “

A Stabilized Finite Element Method for the Incompressible Navier-Stokes Equations Using a Hierarchical Basis

,”

Int. J. Numer. Methods Eng.

,

35

(

1

), pp.

93

–

116

.

Google Scholar

Crossref

4.

Maas

,

S. A.

,

Ellis

,

B. J.

,

Ateshian

,

G. A.

, and

Weiss

,

J. A.

,

2012

, “

FEBio: Finite Elements for Biomechanics

,”

ASME J. Biomech. Eng.

,

134

(

1

), p.

011005

.

Google Scholar

Crossref

5.

Ateshian

,

G. A.

,

Albro

,

M. B.

,

Maas

,

S.

, and

Weiss

,

J. A.

,

2011

, “

Finite Element Implementation of Mechanochemical Phenomena in Neutral Deformable Porous Media Under Finite Deformation

,”

ASME J. Biomech. Eng.

,

133

(

8

), p.

081005

.

Google Scholar

Crossref

6.

Ateshian

,

G. A.

,

Maas

,

S.

, and

Weiss

,

J. A.

,

2013

, “

Multiphasic Finite Element Framework for Modeling Hydrated Mixtures With Multiple Neutral and Charged Solutes

,”

ASME J. Biomech. Eng.

,

135

(

11

), p.

111001

.

Google Scholar

Crossref

7.

Ateshian

,

G. A.

,

Maas

,

S.

, and

Weiss

,

J. A.

,

2010

, “

Finite Element Algorithm for Frictionless Contact of Porous Permeable Media Under Finite Deformation and Sliding

,”

ASME J. Biomech. Eng.

,

132

(

6

), p.

061006

.

Google Scholar

Crossref

8.

Ateshian

,

G. A.

,

Maas

,

S.

, and

Weiss

,

J. A.

,

2012

, “

Solute Transport Across a Contact Interface in Deformable Porous Media

,”

J. Biomech.

,

45

(

6

), pp.

1023

–

1027

.

Google Scholar

Crossref

PubMed

9.

Ateshian

,

G. A.

,

Nims

,

R. J.

,

Maas

,

S.

, and

Weiss

,

J. A.

,

2014

, “

Computational Modeling of Chemical Reactions and Interstitial Growth and Remodeling Involving Charged Solutes and Solid-Bound Molecules

,”

Biomech. Model. Mechanobiol.

,

13

(

5

), pp.

1105

–

1120

.

Google Scholar

Crossref

PubMed

10.

Ateshian

,

G. A.

,

2015

, “

Viscoelasticity Using Reactive Constrained Solid Mixtures

,”

J. Biomech.

,

48

(

6

), pp.

941

–

947

.

Google Scholar

Crossref

PubMed

11.

Nims

,

R. J.

,

Durney

,

K. M.

,

Cigan

,

A. D.

,

Dusséaux

,

A.

,

Hung

,

C. T.

, and

Ateshian

,

G. A.

,

2016

, “

Continuum Theory of Fibrous Tissue Damage Mechanics Using Bond Kinetics: Application to Cartilage Tissue Engineering

,”

Interface Focus

,

6

(

1

), p.

20150063

.

Google Scholar

Crossref

PubMed

12.

Krittian

,

S.

,

Janoske

,

U.

,

Oertel

,

H.

, and

Böhlke

,

T.

,

2010

, “

Partitioned Fluid-Solid Coupling for Cardiovascular Blood Flow: Left-Ventricular Fluid Mechanics

,”

Ann. Biomed. Eng.

,

38

(

4

), pp.

1426

–

1441

.

Google Scholar

Crossref

PubMed

13.

Wolters

,

B. J. B. M.

,

Rutten

,

M. C. M.

,

Schurink

,

G. W. H.

,

Kose

,

U.

,

de Hart

,

J.

, and

van de Vosse

,

F. N.

,

2005

, “

A Patient-Specific Computational Model of Fluid-Structure Interaction in Abdominal Aortic Aneurysms

,”

Med. Eng. Phys.

,

27

(

10

), pp.

871

–

883

.

Google Scholar

Crossref

PubMed

14.

Zhang

,

W.

,

Herrera

,

C.

,

Atluri

,

S. N.

, and

Kassab

,

G. S.

,

2004

, “

Effect of Surrounding Tissue on Vessel Fluid and Solid Mechanics

,”

ASME J. Biomech. Eng.

,

126

(

6

), pp.

760

–

769

.

Google Scholar

Crossref

15.

Moghani

,

T.

,

Butler

,

J. P.

,

Lin

,

J. L.-W.

, and

Loring

,

S. H.

,

2007

, “

Finite Element Simulation of Elastohydrodynamic Lubrication of Soft Biological Tissues

,”

Comput. Struct.

,

85

(

11

), pp.

1114

–

1120

.

Google Scholar

Crossref

PubMed

16.

Masoumi

,

N.

,

Framanzad

,

F.

,

Zamanian

,

B.

,

Seddighi

,

A. S.

,

Moosavi

,

M. H.

,

Najarian

,

S.

, and

Bastani

,

D.

,

2013

, “

2D Computational Fluid Dynamic Modeling of Human Ventricle System Based on Fluid-Solid Interaction and Pulsatile Flow

,”

Basic Clin. Neurosci.

,

4

(

1

), pp.

64

–

75

.https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4202551/

Google Scholar

PubMed

17.

Sweetman

,

B.

, and

Linninger

,

A. A.

,

2011

, “

Cerebrospinal Fluid Flow Dynamics in the Central Nervous System

,”

Ann. Biomed. Eng.

,

39

(

1

), pp.

484

–

496

.

Google Scholar

Crossref

PubMed

18.

Sweetman

,

B.

,

Xenos

,

M.

,

Zitella

,

L.

, and

Linninger

,

A. A.

,

2011

, “

Three-Dimensional Computational Prediction of Cerebrospinal Fluid Flow in the Human Brain

,”

Comput. Biol. Med.

,

41

(

2

), pp.

67

–

75

.

Google Scholar

Crossref

PubMed

19.

Luo

,

H.

,

Mittal

,

R.

,

Zheng

,

X.

,

Bielamowicz

,

S. A.

,

Walsh

,

R. J.

, and

Hahn

,

J. K.

,

2008

, “

An Immersed-Boundary Method for Flow–Structure Interaction in Biological Systems With Application to Phonation

,”

J. Comput. Phys.

,

227

(

22

), pp.

9303

–

9332

.

Google Scholar

Crossref

PubMed

20.

Malvè

,

M.

,

Pérez del Palomar

,

A.

,

Chandra

,

S.

,

López-Villalobos

,

J. L.

,

Mena

,

A.

,

Finol

,

E. A.

,

Ginel

,

A.

, and

Doblaré

,

M.

,

2011

, “

FSI Analysis of a Healthy and a Stenotic Human Trachea Under Impedance-Based Boundary Conditions

,”

ASME J. Biomech. Eng.

,

133

(

2

), p.

021001

.

Google Scholar

Crossref

21.

Barakat

,

A. I.

,

2001

, “

A Model for Shear Stress-Induced Deformation of a Flow Sensor on the Surface of Vascular Endothelial Cells

,”

J. Theor. Biol.

,

210

(

2

), pp.

221

–

236

.

Google Scholar

Crossref

PubMed

22.

Fung

,

Y. C.

, and

Liu

,

S. Q.

,

1993

, “

Elementary Mechanics of the Endothelium of Blood Vessels

,”

ASME J. Biomech. Eng.

,

115

(

1

), pp.

1

–

12

.

Google Scholar

Crossref

23.

Fritton

,

S. P.

, and

Weinbaum

,

S.

,

2009

, “

Fluid and Solute Transport in Bone: Flow-Induced Mechanotransduction

,”

Annu. Rev. Fluid Mech.

,

41

, pp.

347

–

374

.

Google Scholar

Crossref

PubMed

24.

Kwon

,

R. Y.

, and

Frangos

,

J. A.

,

2010

, “

Quantification of Lacunar-Canalicular Interstitial Fluid Flow Through Computational Modeling of Fluorescence Recovery After Photobleaching

,”

Cell Mol. Bioeng.

,

3

(

3

), pp.

296

–

306

.

Google Scholar

Crossref

PubMed

25.

Scheiner

,

S.

,

Pivonka

,

P.

, and

Hellmich

,

C.

,

2016

, “

Poromicromechanics Reveals That Physiological Bone Strains Induce Osteocyte-Stimulating Lacunar Pressure

,”

Biomech. Model. Mechanobiol.

,

15

(

1

), pp.

9

–

28

.

Google Scholar

Crossref

PubMed

26.

Avanzini

,

A.

, and

Battini

,

D.

,

2014

, “

Structural Analysis of a Stented Pericardial Heart Valve With Leaflets Mounted Externally

,”

Proc. Inst. Mech. Eng. H

,

228

(

10

), pp.

985

–

995

.

Google Scholar

Crossref

PubMed

27.

Fiore

,

G. B.

,

Redaelli

,

A.

,

Guadagni

,

G.

,

Inzoli

,

F.

, and

Fumero

,

R.

,

2002

, “

Development of a New Disposable Pulsatile Pump for Cardiopulmonary Bypass: Computational Fluid-Dynamic Design and In Vitro Tests

,”

ASAIO J.

,

48

(

3

), pp.

260

–

267

.

Google Scholar

Crossref

PubMed

28.

Malvè

,

M.

,

Del Palomar

,

A. P.

,

Chandra

,

S.

,

López-Villalobos

,

J. L.

,

Finol

,

E. A.

,

Ginel

,

A.

, and

Doblaré

,

M.

,

2011

, “

FSI Analysis of a Human Trachea Before and After Prosthesis Implantation

,”

ASME J. Biomech. Eng.

,

133

(

7

), p.

071003

.

Google Scholar

Crossref

29.

Malkus

,

D. S.

, and

Hughes

,

T. J.

,

1978

, “

Mixed Finite Element Methods—Reduced and Selective Integration Techniques: A Unification of Concepts

,”

Comput. Methods Appl. Mech. Eng.

,

15

(

1

), pp.

63

–

81

.

Google Scholar

Crossref

30.

Hughes

,

T. J.

,

Liu

,

W. K.

, and

Brooks

,

A.

,

1979

, “

Finite Element Analysis of Incompressible Viscous Flows by the Penalty Function Formulation

,”

J. Comput. Phys.

,

30

(

1

), pp.

1

–

60

.

Google Scholar

Crossref

31.

Bercovier

,

M.

, and

Engelman

,

M.

,

1979

, “

A Finite Element for the Numerical Solution of Viscous Incompressible Flows

,”

J. Comput. Phys.

,

30

(

2

), pp.

181

–

201

.

Google Scholar

Crossref

32.

Reddy

,

J. N.

, and

Gartling

,

D. K.

,

2001

,

The Finite Element Method in Heat Transfer and Fluid Dynamics

, 2nd ed.,

CRC Press

,

Boca Raton, FL

.

33.

Brezzi

,

F.

, and

Bathe

,

K.-J.

,

1986

, “

The Inf-Sup Condition, Equivalent Forms and Applications

,”

International Conference Reliability of Methods for Engineering Analysis

, Swansea, UK, July 9–11, pp. 197–219http://web.mit.edu/kjb/www/Publications_Prior_to_1998/Studies_of_Finite_Element_Procedures_The_Inf-Sup_Condition_Equivalent_Forms_and_Applications.pdf.

34.

Brezzi

,

F.

, and

Bathe

,

K.-J.

,

1990

, “

A Discourse on the Stability Conditions for Mixed Finite Element Formulations

,”

Comput. Methods Appl. Mech. Eng.

,

82

(

1–3

), pp.

27

–

57

.

Google Scholar

Crossref

35.

Brooks

,

A. N.

, and

Hughes

,

T. J.

,

1982

, “

Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows With Particular Emphasis on the Incompressible Navier-Stokes Equations

,”

Comput. Methods Appl. Mech. Eng.

,

32

(

1

), pp.

199

–

259

.

Google Scholar

Crossref

36.

Hughes

,

T. J.

,

Franca

,

L. P.

, and

Balestra

,

M.

,

1986

, “

A New Finite Element Formulation for Computational Fluid Dynamics: V. Circumventing the Babuška-Brezzi Condition: A Stable Petrov-Galerkin Formulation of the Stokes Problem Accommodating Equal-Order Interpolations

,”

Comput. Methods Appl. Mech. Eng.

,

59

(

1

), pp.

85

–

99

.

Google Scholar

Crossref

37.

Hughes

,

T. J.

,

Franca

,

L. P.

, and

Hulbert

,

G. M.

,

1989

, “

A New Finite Element Formulation for Computational Fluid Dynamics—VIII: The Galerkin/Least-Squares Method for Advective-Diffusive Equations

,”

Comput. Methods Appl. Mech. Eng.

,

73

(

2

), pp.

173

–

189

.

Google Scholar

Crossref

38.

Tezduyar

,

T. E.

,

1991

, “

Stabilized Finite Element Formulations for Incompressible Flow Computations

,”

Adv. Appl. Mech.

,

28

, pp.

1

–

44

.

Google Scholar

Crossref

39.

Franca

,

L. P.

, and

Frey

,

S. L.

,

1992

, “

Stabilized Finite Element Methods—II: The Incompressible Navier-Stokes Equations

,”

Comput. Methods Appl. Mech. Eng.

,

99

(

2–3

), pp.

209

–

233

.

Google Scholar

Crossref

40.

Tezduyar

,

T. E.

,

Mittal

,

S.

,

Ray

,

S.

, and

Shih

,

R.

,

1992

, “

Incompressible Flow Computations With Stabilized Bilinear and Linear Equal-Order-Interpolation Velocity-Pressure Elements

,”

Comput. Methods Appl. Mech. Eng.

,

95

(

2

), pp.

221

–

242

.

Google Scholar

Crossref

41.

Tezduyar

,

T. E.

, and

Osawa

,

Y.

,

2000

, “

Finite Element Stabilization Parameters Computed From Element Matrices and Vectors

,”

Comput. Methods Appl. Mech. Eng.

,

190

(

3

), pp.

411

–

430

.

Google Scholar

Crossref

42.

Akin

,

J. E.

, and

Tezduyar

,

T. E.

,

2004

, “

Calculation of the Advective Limit of the SUPG Stabilization Parameter for Linear and Higher-Order Elements

,”

Comput. Methods Appl. Mech. Eng.

,

193

(

21

), pp.

1909

–

1922

.

Google Scholar

Crossref

43.

Bonet

,

J.

, and

Wood

,

R. D.

,

1997

,

Nonlinear Continuum Mechanics for Finite Element Analysis

,

Cambridge University Press

,

New York

.

44.

Jansen

,

K. E.

,

Whiting

,

C. H.

, and

Hulbert

,

G. M.

,

2000

, “

A Generalized-α Method for Integrating the Filtered Navier–Stokes Equations With a Stabilized Finite Element Method

,”

Comput. Methods Appl. Mech. Eng.

,

190

(

3

), pp.

305

–

319

.

Google Scholar

Crossref

45.

Reddy

,

J. N.

,

2008

,

An Introduction to Continuum Mechanics: With Applications

,

Cambridge University Press

,

New York

.

46.

Panton

,

R. L.

,

2006

,

Incompressible Flow

,

Wiley

, New York.

47.

Cho

,

Y. I.

, and

Kensey

,

K. R.

,

1991

, “

Effects of the Non-Newtonian Viscosity of Blood on Flows in a Diseased Arterial Vessel—Part 1: Steady Flows

,”

Biorheology

,

28

(

3–4

), pp.

241

–

262

.

Google Scholar

Crossref

PubMed

48.

Ryzhakov

,

P. B.

,

Rossi

,

R.

,

Idelsohn

,

S. R.

, and

Onate

,

E.

,

2010

, “

A Monolithic Lagrangian Approach for Fluid-Structure Interaction Problems

,”

Comput. Mech.

,

46

(

6

), pp.

883

–

899

.

Google Scholar

Crossref

49.

Zienkiewicz

,

O. C.

, and

Codina

,

R.

,

1995

, “

A General Algorithm for Compressible and Incompressible Flow—Part I: The Split, Characteristic-Based Scheme

,”

Int. J. Numer. Methods Eng.

,

20

(

8–9

), pp.

869

–

885

.

Google Scholar

Crossref

50.

Zienkiewicz

,

O. C.

,

Taylor

,

R. L.

, and

Nithiarasu

,

P.

,

2014

,

The Finite Element Method for Fluid Dynamics

, 7th ed.,

Butterworth-Heinemann

,

Oxford, UK.

.

51.

Bazilevs

,

Y.

,

Gohean

,

J.

,

Hughes

,

T.

,

Moser

,

R.

, and

Zhang

,

Y.

,

2009

, “

Patient-Specific Isogeometric Fluid–Structure Interaction Analysis of Thoracic Aortic Blood Flow Due to Implantation of the Jarvik 2000 Left Ventricular Assist Device

,”

Comput. Methods Appl. Mech. Eng.

,

198

(

45

), pp.

3534

–

3550

.

Google Scholar

Crossref

52.

Esmaily Moghadam

,

M.

,

Bazilevs

,

Y.

,

Hsia

,

T.-Y.

,

Vignon-Clementel

,

I. E.

, and

Marsden

,

A. L.

,

2011

, “

A Comparison of Outlet Boundary Treatments for Prevention of Backflow Divergence With Relevance to Blood Flow Simulations

,”

Comput. Mech.

,

48

(

3

), pp.

277

–

291

.

Google Scholar

Crossref

53.

Vignon-Clementel

,

I. E.

,

Figueroa

,

C. A.

,

Jansen

,

K. E.

, and

Taylor

,

C. A.

,

2006

, “

Outflow Boundary Conditions for Three-Dimensional Finite Element Modeling of Blood Flow and Pressure in Arteries

,”

Comput. Methods Appl. Mech. Eng.

,

195

(

29

), pp.

3776

–

3796

.

Google Scholar

Crossref

54.

Engelman

,

M.

,

Strang

,

G.

, and

Bathe

,

K.-J.

,

1981

, “

The Application of Quasi-Newton Methods in Fluid Mechanics

,”

Int. J. Numer. Methods Eng.

,

17

(

5

), pp.

707

–

718

.

Google Scholar

Crossref

55.

Schenk

,

O.

, and

Gärtner

,

K.

,

2004

, “

Solving Unsymmetric Sparse Systems of Linear Equations With PARDISO

,”

Future Gener. Comput. Syst.

,

20

(

3

), pp.

475

–

487

.

Google Scholar

Crossref

56.

Ghia

,

U.

,

Ghia

,

K. N.

, and

Shin

,

C.

,

1982

, “

High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method

,”

J. Comput. Phys.

,

48

(

3

), pp.

387

–

411

.

Google Scholar

Crossref

57.

Denham

,

M.

, and

Patrick

,

M.

,

1974

, “

Laminar Flow Over a Downstream-Facing Step in a Two-Dimensional Flow Channel

,”

Trans. Inst. Chem. Eng.

,

52

(

4

), pp.

361

–

367

.http://archive.icheme.org/cgi-bin/somsid.cgi?session=924344F&page=52ap0361-001&type=framedpdf

58.

Seeley

,

B. D.

, and

Young

,

D. F.

,

1976

, “

Effect of Geometry on Pressure Losses Across Models of Arterial Stenoses

,”

J. Biomech.

,

9

(

7

), pp.

439

–

448

.

Google Scholar

Crossref

PubMed

59.

Tezduyar

,

T. E.

,

1992

,

Stabilized Finite Element Formulations for Incompressible Flow Computations

,

Academic Press

, New York.

60.

Pelletier

,

D.

,

Fortin

,

A.

, and

Camarero

,

R.

,

1989

, “

Are Fem Solutions of Incompressible Flows Really Incompressible? (or How Simple Flows Can Cause Headaches!)

,”

Int. J. Numer. Methods Fluids

,

9

(

1

), pp.

99

–

112

.

Google Scholar

Crossref

Finite Element Framework for Computational Fluid Dynamics in FEBio

Introduction

Finite Element Implementation

Governing Equations.

Virtual Work and Weak Form.

Temporal Discretization and Linearization.

Constitutive Relations.

Spatial Discretization.

Condition Number.

Special Boundary Conditions

Backflow Stabilization.

Tangential Flow Stabilization.

Flow Resistance.

Solver.

Verifications and Validations

Lid-Driven Cavity Flow.

Flow Past a Backward Facing Step.

Inviscid Flow One-Dimensional Wave Propagation.

Flow Past Block in a Narrow Channel.

Carotid Bifurcation.

Idealized Model of Arterial Stenosis.

Discussion

Disclaimer

Funding Data

Appendix

Energy Balance

Generalized α-Method

References

Supplementary data

Get Email Alerts

Cited By

Journals

Conference Proceedings

eBooks

Resources

Opportunities

Finite Element Framework for Computational Fluid Dynamics in FEBio

Introduction

Finite Element Implementation

Governing Equations.

Virtual Work and Weak Form.

Temporal Discretization and Linearization.

Constitutive Relations.

Spatial Discretization.

Condition Number.

Special Boundary Conditions

Backflow Stabilization.

Tangential Flow Stabilization.

Flow Resistance.

Solver.

Verifications and Validations

Lid-Driven Cavity Flow.

Flow Past a Backward Facing Step.

Inviscid Flow One-Dimensional Wave Propagation.

Flow Past Block in a Narrow Channel.

Carotid Bifurcation.

Idealized Model of Arterial Stenosis.

Discussion

Disclaimer

Funding Data

Appendix

Energy Balance

Generalized α-Method

References

Supplementary data

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

Journals

Conference Proceedings

eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only