3.5 Computational Fluid Dynamics

A more detailed description of the FEBio fluid solver can be found in ¹.

Weak Formulation

The nodal unknowns in this formulation are \(\mathbf{v}\) and \(J\) (or \(e\)), which may be solved using the momentum balance in eq.(2.13-1) and the kinematic constraint between \(J\) and \(\mathbf{v}\) given in eq.(2.13-6). The virtual work integral for a Galerkin finite element formulation ² is given by

\[ \begin{equation} \begin{aligned}\delta W & =\int_{\Omega}\delta\mathbf{v}\cdot\left(\divg\boldsymbol{\sigma}+\rho\left(\mathbf{b}-\mathbf{a}\right)\right)dv\\ & +\int_{\Omega}\delta J\left(\frac{\dot{J}}{J}-\divg\mathbf{v}\right)dv\,, \end{aligned} \label{eq:virtual work} \end{equation} \]

where \(\delta\mathbf{v}\) is a virtual velocity and \(\delta J\) is a virtual energy density; \(\Omega\) is the fluid finite element domain and \(dv\) is a differential volume in \(\Omega\). 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 Section Energy Balance). Using the divergence theorem, we may rewrite the weak form of this integral as the difference between external and internal virtual work, \(\delta W=\delta W_{ext}-\delta W_{int}\), where

\[ \begin{equation} \begin{aligned}\delta W_{int} & =\int_{\Omega}\boldsymbol{\tau}:\grad\delta\mathbf{v}\,dv+\int_{\Omega}\delta\mathbf{v}\cdot\left(\grad p+\rho\mathbf{a}\right)\,dv\\ & -\int_{\Omega}\left(\delta J\frac{\dot{J}}{J}+\grad\delta J\cdot\mathbf{v}\right)\,dv\,, \end{aligned} \label{eq:virtual-work-internal-1} \end{equation} \]

and

\[ \begin{equation} \delta W_{ext}=\int_{\partial\Omega}\delta\mathbf{v}\cdot\mathbf{t}^{\tau}\,da+\int_{\Omega}\delta\mathbf{v}\cdot\rho\mathbf{b}\,dv-\int_{\partial\Omega}\delta J\,v_{n}\,da\,.\label{eq:virtual-work-external-1} \end{equation} \]

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

Temporal Discretization and Linearization

The time derivatives, \(\partial\mathbf{v}/\partial t\) which appears in the expression for \(\mathbf{a}\) in eq.(2.13-2), 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-\(\alpha\) method ³ (Section Generalized \(\alpha-\)Method). In this scheme, \(\delta W\) is evaluated at an intermediate time step \(t_{n+\alpha}=\alpha t_{n+1}+\left(1-\alpha\right)t_{n}\) between the current time step \(t_{n+1}\) and previous time step \(t_{n}\), though different values of \(\alpha\) are used for the primary variables and their time derivatives. The velocity and volume ratio are evaluated as \(\mathbf{v}_{n+\alpha_{f}}\) and \(J_{n+\alpha_{f}}\) at the intermediate time step \(t_{n+\alpha_{f}}\), whereas their time derivatives are evaluated as \(\left(\partial\mathbf{v}/\partial t\right)_{n+\alpha_{m}}\) and \(\left(\partial J/\partial t\right)_{n+\alpha_{m}}\) at the intermediate time step \(t_{n+\alpha_{m}}\). The parameters \(\alpha_{f}\) and \(\alpha_{m}\) are evaluated from the spectral radius for an infinite time step, \(\rho_{\infty}\), as described in Section Generalized \(\alpha-\)Method. The solution of the nonlinear equation \(\delta W=0\) is obtained by linearizing this relation as

\[ \begin{equation} \delta W+D\delta W\left[\Delta\mathbf{v}\right]+D\delta W\left[\Delta J\right]\approx0\,,\label{eq:linearized-virtual-work} \end{equation} \]

where the operator \(D\delta W\left[\cdot\right]\) represents the directional derivative of \(\delta W\) at \(\left(\mathbf{v},J\right)\) along an increment \(\Delta\mathbf{v}\) of \(\mathbf{v}\), or \(\Delta J\) of \(J\) ². The aim of this analysis is to solve for the velocity \(\mathbf{v}_{n+1}\) and volume ratio \(J_{n+1}\) at the current time step \(t_{n+1}\). Using the split form of \(\delta W\) between external and internal work contributions, this relation may be expanded as

\[ \begin{equation} \begin{aligned} & D\delta W_{int}\left[\Delta\mathbf{v}\right]+D\delta W_{int}\left[\Delta J\right]-D\delta W_{ext}\left[\Delta\mathbf{v}\right]\\ & -D\delta W_{ext}\left[\Delta J\right]\approx\delta W_{ext}-\delta W_{int}\,. \end{aligned} \label{eq:linearized-work-split} \end{equation} \]

In this framework the finite element mesh is defined on the spatial domain \(\Omega\), which is fixed (time-invariant) in conventional CFD treatments. Thus, we can linearize \(\delta W_{int}\) along increments \(\Delta\mathbf{v}\) in the velocity \(\mathbf{v}_{n+1}\) and \(\Delta J\) in the volume ratio \(J_{n+1}\), by simply bringing the directional derivative operator into the integrals of eqs.\eqref{eq:virtual-work-internal-1}-\eqref{eq:virtual-work-external-1}. The linearization of \(\mathbf{v}_{n+\alpha_{f}}\) and \(J_{n+\alpha_{f}}\) is given by

\[ \begin{aligned}D\mathbf{v}_{n+\alpha_{f}}\left[\Delta\mathbf{v}\right] & =\alpha_{f}\Delta\mathbf{v}\\ DJ_{n+\alpha_{f}}\left[\Delta J\right] & =\alpha_{f}\Delta J \end{aligned} \]

whereas that of \(\left(\partial\mathbf{v}/\partial t\right)_{n+\alpha_{m}}\) and \(\left(\partial J/\partial t\right)_{n+\alpha_{m}}\) is given by

\[ \begin{equation} D\left(\frac{\partial\mathbf{v}}{dt}\right)_{n+\alpha_{m}}\left[\Delta\mathbf{v}\right]=\frac{\alpha_{m}}{\gamma}\frac{\Delta\mathbf{v}}{\Delta t}\,,\label{eq:discretized-a} \end{equation} \]

\[ \begin{equation} D\left(\frac{\partial J}{\partial t}\right)_{n+\alpha_{m}}\left[\Delta\mathbf{v}\right]=\frac{\alpha_{m}}{\gamma}\frac{\Delta J}{\Delta t}\,.\label{eq:discretized-J} \end{equation} \]

Here, \(\Delta t\) is the current time increment and \(\gamma\) is the Newmark integration parameter ³.

The linearization of \(\delta W_{int}\) along an increment \(\Delta\mathbf{v}\) is then

\[ \begin{equation} \begin{aligned}D\left(\delta W_{int}\right)\left[\Delta\mathbf{v}\right] & =\alpha_{f}\int_{\Omega}\grad\delta\mathbf{v}:\boldsymbol{\mathcal{C}}^{\tau}:\grad\Delta\mathbf{v}\,dv\\ & +\alpha_{f}\int_{\Omega}\delta\mathbf{v}\cdot\rho\left(\left(\frac{\alpha_{m}}{\alpha_{f}\gamma}\frac{\mathbf{I}}{\Delta t}+\mathbf{L}\right)\cdot\Delta\mathbf{v}+\grad\Delta\mathbf{v}\cdot\mathbf{v}\right)\,dv\\ & -\alpha_{f}\int_{\Omega}\left(\frac{\delta J}{J}\grad J+\grad\delta J\right)\cdot\Delta\mathbf{v}\,dv \end{aligned} \label{eq:Wint-linearization-v} \end{equation} \]

where we have introduced the fourth-order tensor \(\boldsymbol{\mathcal{C}}^{\tau}\) representing the tangent of the viscous stress with respect to the rate of deformation,

\[ \begin{equation} \boldsymbol{\mathcal{C}}^{\tau}=\frac{\partial\boldsymbol{\tau}}{\partial\mathbf{D}}\,.\label{eq:stress-tangent-D} \end{equation} \]

Note that \(\boldsymbol{\mathcal{C}}^{\tau}\) exhibits minor symmetries because of the symmetries of \(\boldsymbol{\tau}\) and \(\mathbf{D}\); in Cartesian components, we have \(\mathcal{C}_{ijkl}^{\tau}=\mathcal{C}_{jikl}^{\tau}\) and \(\mathcal{C}_{ijkl}^{\tau}=\mathcal{C}_{ijlk}^{\tau}\). In general, \(\boldsymbol{\mathcal{C}}^{\tau}\) does not exhibit major symmetry (\(\mathcal{C}_{ijkl}^{\tau}\ne\mathcal{C}_{klij}^{\tau}\)), though the common constitutive relations adopted in fluid mechanics produce such symmetry as shown below.

The linearization of \(\delta W_{int}\) along an increment \(\Delta J\) is

\[ \begin{equation} \begin{aligned}D\left(\delta W_{int}\right)\left[\Delta J\right] & =\alpha_{f}\int_{b}\Delta J\,\boldsymbol{\tau}_{J}^{\prime}:\grad\delta\mathbf{v}\,dv-\alpha_{f}\int_{\Omega}\delta\mathbf{v}\cdot\Delta J\frac{\rho}{J}\mathbf{a}\,dv\\ & +\alpha_{f}\int_{\Omega}\delta\mathbf{v}\cdot\left(p^{\prime}\grad\Delta J+\Delta J\,p^{\prime\prime}\grad J\right)\,dv\\ & -\alpha_{f}\int_{\Omega}\frac{\delta J}{J}\left(\left(\frac{\alpha_{m}}{\alpha_{f}\gamma}\frac{1}{\Delta t}-\frac{\dot{J}}{J}\right)\Delta J+\grad\Delta J\cdot\mathbf{v}\right)\,dv \end{aligned} \label{eq:Wint-linearization-J} \end{equation} \]

where we have used \(DJ\left[\Delta J\right]=\Delta J\); \(p^{\prime}\) and \(p^{\prime\prime}\) respectively represent the first and second derivatives of \(p\left(J\right)\). We have also defined \(\boldsymbol{\tau}_{J}^{\prime}\) as the tangent of the viscous stress \(\boldsymbol{\tau}\) with respect to \(J\),

\[ \begin{equation} \boldsymbol{\tau}_{J}^{\prime}=\frac{\partial\boldsymbol{\tau}}{\partial J}\,.\label{eq:stress-tangent-J} \end{equation} \]

For the external work, when \(\mathbf{t}^{\tau}\), \(\mathbf{b}\) and \(v_{n}\) are prescribed, these linearizations simplify to

\[ \begin{equation} D\left(\delta W_{ext}\right)\left[\Delta\mathbf{v}\right]=0\,,\label{eq:Wext-linearization-v} \end{equation} \]

and

\[ \begin{equation} D\left(\delta W_{ext}\right)\left[\Delta J\right]=-\alpha_{f}\int_{b}\delta\mathbf{v}\cdot\Delta J\,\frac{\rho}{J}\mathbf{b}\,dv\,.\label{eq:Wext-linearization-J} \end{equation} \]

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 equations are minimal, simply requiring the substitution \(J=1+e\) and \(\Delta J=\Delta e\). Steady-state analyses may be obtained by setting the terms involving \(\Delta t^{-1}\) to zero in eqs.\eqref{eq:discretized-a}-\eqref{eq:discretized-J}, \eqref{eq:Wint-linearization-v} and \eqref{eq:Wint-linearization-J}.

Spatial Discretization

The velocity \(\mathbf{v}\left(\mathbf{x},t\right)\) and Jacobian \(J\left(\mathbf{x},t\right)\) are spatially interpolated over the domain \(\Omega\) using the same interpolation functions \(N_{a}\left(\mathbf{x}\right)\), with \(a=1\) to \(n\) where \(n\) is the number of nodes in an element),

\[ \begin{equation} \mathbf{v}\left(\mathbf{x},t\right)=\sum_{a=1}^{n}N_{a}\left(\mathbf{x}\right)\mathbf{v}_{a}\,,\quad J\left(\mathbf{x},t\right)=\sum_{a=1}^{n}N_{a}\left(\mathbf{x}\right)J_{a}\,.\label{eq:v-J-interpolation} \end{equation} \]

Here, \(\mathbf{v}_{a}\) and \(J_{a}\) are nodal values of \(\mathbf{v}\) and \(J\) that evolve with time. In contrast to classical mixed formulations for incompressible flow ⁴, which solve for the pressure \(p\) using \(\divg\mathbf{v}=0\) instead of eq.(2.13-6), equal order interpolation is acceptable in this formulation since the governing equations for \(\mathbf{v}\) and \(J\) involve spatial derivatives of both variables (\(\grad\mathbf{v}\) and \(\grad J\)). The expressions of eq.\eqref{eq:v-J-interpolation} may be used to evaluate \(\mathbf{L}\), \(\divg\mathbf{v}\), \(\mathbf{a}\), \(\grad J\), \(\dot{J}\), etc. Similar interpolations are used for virtual increments \(\delta\mathbf{v}\) and \(\delta J\), as well as real increments \(\Delta\mathbf{v}\) and \(\Delta J\).

When substituted into eq.(4.2-32), we find that the discretized form of \(\delta W_{int}\) may be written as

\[ \begin{equation} \delta W_{int}=\sum_{a}\delta\mathbf{v}_{a}\cdot\left(\mathbf{f}_{a}^{\sigma}+\mathbf{f}_{a}^{\rho}\right)+f_{a}^{J}\delta J_{a}\,,\label{eq:discretized-Wint} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{f}_{a}^{\sigma} & =\int_{\Omega}\left(\boldsymbol{\tau}\cdot\grad N_{a}+N_{a}\grad p\right)\,dv\,,\\ \mathbf{f}_{a}^{\rho} & =\int_{\Omega}N_{a}\rho\mathbf{a}\,dv\,,\\ f_{a}^{J} & =\int_{\Omega}-\left(N_{a}\frac{\dot{J}}{J}+\grad N_{a}\cdot\mathbf{v}\right)\,dv\,. \end{aligned} \label{eq:discretized-internal-forces} \end{equation} \]

Similarly, the discretized form of \(D\delta W_{int}\left[\Delta\mathbf{v}\right]\) in eq.\eqref{eq:Wint-linearization-v} becomes

\[ \begin{equation} \begin{aligned}D\left(\delta W_{int}\right)\left[\Delta\mathbf{v}\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\left(\mathbf{K}_{ab}^{vv}+\mathbf{M}_{ab}^{vv}\right)\cdot\Delta\mathbf{v}_{b}\\ & +\sum_{a}\delta J_{a}\sum_{b}\mathbf{k}_{ab}^{Jv}\cdot\Delta\mathbf{v}_{b}\,, \end{aligned} \label{eq:discretized-DWint-v} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab}^{vv} & =\alpha_{f}\int_{\Omega}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{v}\cdot\grad N_{b}\,dv\,,\\ \mathbf{M}_{ab}^{vv} & =\alpha_{f}\int_{\Omega}N_{a}\rho\left(N_{b}\left(\frac{\alpha_{m}}{\alpha_{f}\gamma}\frac{\mathbf{I}}{\Delta t}+\grad\mathbf{v}\right)+\left(\grad N_{b}\cdot\mathbf{v}\right)\mathbf{I}\right)\,dv\,,\\ \mathbf{k}_{ab}^{Jv} & =-\alpha_{f}\int_{\Omega}\left(\frac{N_{a}}{J}\grad J+\grad N_{a}\right)N_{b}\,dv\,, \end{aligned} \label{eq:discretized-stiffness-v} \end{equation} \]

whereas that of \(D\delta W_{int}\left[\Delta J\right]\) in eq.\eqref{eq:Wint-linearization-J} becomes

\[ \begin{equation} \begin{aligned}D\left(\delta W_{int}\right)\left[\Delta J\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\left(\mathbf{k}_{ab}^{vJ}+\mathbf{m}_{ab}^{vJ}\right)\,\Delta J_{b}\\ & +\sum_{a}\delta J_{a}\sum_{b}k_{ab}^{JJ}\,\Delta J_{b} \end{aligned} \label{eq:discretized-DWint-J} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{k}_{ab}^{vJ} & =\alpha_{f}\int_{\Omega}\left[N_{b}\boldsymbol{\tau}_{J}^{\prime}\cdot\grad N_{a}+N_{a}\left(p^{\prime}\grad N_{b}+N_{b}p^{\prime\prime}\grad J\right)\right]\,dv\,,\\ \mathbf{m}_{ab}^{vJ} & =-\alpha_{f}\int_{\Omega}N_{a}N_{b}\frac{\rho}{J}\mathbf{a}\,dv\,,\\ k_{ab}^{JJ} & =-\alpha_{f}\int_{\Omega}\frac{N_{a}}{J}\left(\left(\frac{\alpha_{m}}{\alpha_{f}\gamma}\frac{1}{\Delta t}-\frac{\dot{J}}{J}\right)N_{b}+\grad N_{b}\cdot\mathbf{v}\right)\,dv\,. \end{aligned} \label{eq:discretized-stiffness-J} \end{equation} \]

For the external work in eq.(4.2-34), its discretized form is

\[ \begin{equation} \delta W_{ext}=\sum_{a}\delta\mathbf{v}_{a}\cdot\left(\mathbf{f}_{a}^{\text{t}}+\mathbf{f}_{a}^{\text{b}}\right)+\delta J_{a}f_{a}^{v}\,,\label{eq:discretized-Wext} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{f}_{a}^{\text{t}} & =\int_{\partial\Omega}N_{a}\mathbf{t}^{\tau}\,da\,,\\ \mathbf{f}_{a}^{\text{b}} & =\int_{\Omega}N_{a}\rho\mathbf{b}\,dv\,,\\ f_{a}^{v} & =-\int_{\Omega}N_{a}v_{n}\,da\,. \end{aligned} \label{eq:discretized-external-forces} \end{equation} \]

The discretized form of \(D\left(\delta W_{ext}\right)\left[\Delta J\right]\) in eq.\eqref{eq:Wext-linearization-J} is

\[ \begin{equation} D\left(\delta W_{ext}\right)\left[\Delta J\right]=\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{k}_{ab}^{\text{b}}\,\Delta J_{b}\,,\label{eq:discretized-DWext-J} \end{equation} \]

where

\[ \begin{equation} \mathbf{k}_{ab}^{\text{b}}=-\alpha_{f}\int_{\Omega}N_{a}N_{b}\frac{\rho}{J}\mathbf{b}\,dv\,.\label{eq:discretized-stiffness-ext} \end{equation} \]

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 ⁵⁶; for these types of problems, letting \(\mathbf{t}=\mathbf{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, \(\mathbf{t}=\beta\rho\left(\mathbf{v}\otimes\mathbf{v}\right)\cdot\mathbf{n}\) with a tensile normal component, that counters the backflow (only when \(v_{n}<0\)). Here, \(\beta\) is a non-dimensional user-defined parameter; a value of \(\beta=0\) turns off this feature, while a value of \(\beta=1\) generally shows good numerical performance. We adapt this previously proposed formulation by letting the normal component of the viscous traction be given by

\[ \begin{equation} t_{n}^{\tau}=\begin{cases} \beta\rho_{r}v_{n}^{2} & v_{n}<0\\ 0 & v_{n}\ge0 \end{cases}\,.\label{eq:BFS-normal-traction} \end{equation} \]

The choice of \(\rho_{r}\) in lieu of \(\rho\) is for convenience, to avoid the dependence of \(\rho\) on \(J\) (which is negligible for nearly incompressible flow). Then, the contribution of this traction to the virtual external work \(\delta W_{ext}\) is

\[ \begin{equation} \delta G=\int_{\partial\Omega}\delta\mathbf{v}\cdot t_{n}^{\tau}\mathbf{n}\,da\,.\label{eq:BFSckflow-virtual-work} \end{equation} \]

The linearization of \(\delta G\) along an increment \(\Delta\mathbf{v}\) in the velocity is given by

\[ \begin{equation} D\delta G\left[\Delta\mathbf{v}\right]=\int_{\partial\Omega}\delta\mathbf{v}\cdot\mathbf{K}_{n}\cdot\Delta\mathbf{v}\,da\,,\label{eq:BFS-linearization} \end{equation} \]

where

\[ \begin{equation} \mathbf{K}_{n}=\begin{cases} 2\beta\rho_{r}v_{n}\left(\mathbf{n}\otimes\mathbf{n}\right) & v_{n}<0\\ \mathbf{0} & v_{n}\ge0 \end{cases}\,.\label{eq:BFS-stiffness} \end{equation} \]

The discretized form of \(\delta G\) is

\[ \begin{equation} \delta G=\sum_{a}\delta\mathbf{v}_{a}\cdot\mathbf{f}_{a}^{n}\,,\quad\mathbf{f}_{a}^{n}=\int_{\partial\Omega}N_{a}t_{n}^{\tau}\mathbf{n}\,da\,,\label{eq:BFS-discretized-work} \end{equation} \]

whereas the discretized form of \(D\delta G\left[\Delta\mathbf{v}\right]\) is

\[ \begin{equation} D\delta G\left[\Delta\mathbf{v}\right]=\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{n}\cdot\Delta\mathbf{v}_{b}\,,\quad\mathbf{K}_{ab}^{n}=\int_{\partial\Omega}N_{a}N_{b}\mathbf{K}_{n}\,da\,.\label{eq:BFS-discretized-tangent} \end{equation} \]

A (viscous) tangential traction is implemented as a separate flow stabilization method in the next section, 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 \(\mathbf{t}_{t}^{\tau}=\mathbf{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 above.

Similar to the previous section, we introduce a non-dimensional parameter \(\beta\), with the tangential traction given by

\[ \begin{equation} \mathbf{t}_{t}^{\tau}=-\beta\rho_{r}\left|\mathbf{v}_{t}\right|\mathbf{v}_{t}\,.\label{eq:tangential-stabilization} \end{equation} \]

This form shows that \(\mathbf{t}_{t}^{\tau}\) opposes tangential flow. The external virtual work for this traction is

\[ \begin{equation} \delta G=\int_{\partial\Omega}\delta\mathbf{v}\cdot\mathbf{t}_{t}^{\tau}\,da\,.\label{eq:TS-external-work} \end{equation} \]

Its linearization along an increment \(\Delta\mathbf{v}\) is

\[ \begin{equation} D\delta G\left[\Delta\mathbf{v}\right]=\int_{\partial\Omega}\delta\mathbf{v}\cdot\mathbf{K}_{t}\cdot\Delta\mathbf{v}\,da\,,\label{eq:TS-linearization} \end{equation} \]

where it can be shown that

\[ \begin{equation} \mathbf{K}_{t}=-\beta\rho_{r}\left|\mathbf{v}_{t}\right|\left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}+\frac{\mathbf{v}_{t}}{\left|\mathbf{v}_{t}\right|}\otimes\frac{\mathbf{v}_{t}}{\left|\mathbf{v}_{t}\right|}\right)\,.\label{eq:TS-stiffness} \end{equation} \]

The discretized form of \(\delta G\) is

\[ \begin{equation} \delta G=\sum_{a}\delta\mathbf{v}_{a}\cdot\mathbf{f}_{a}^{\tau}\,,\quad\mathbf{f}_{a}^{\tau}=\int_{\partial\Omega}N_{a}\mathbf{t}_{t}^{\tau}\,da\,.\label{eq:TS-discretized-work} \end{equation} \]

The discretized form of \(D\delta G\left[\Delta\mathbf{v}\right]\) is

\[ \begin{equation} D\delta G\left[\Delta\mathbf{v}\right]=\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{t}\cdot\Delta\mathbf{v}_{b}\,,\quad\mathbf{K}_{ab}^{t}=\int_{\partial\Omega}N_{a}N_{b}\mathbf{K}_{t}\,da\,.\label{eq:TS-discretized-tangent} \end{equation} \]

Flow Resistance

Flow resistance is typically implemented when modeling arterial flow, where the finite element domain only describes a portion of an arterial network ⁷. 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=RQ\,,\quad Q=\int_{\partial\Omega}v_{n}\,da\,, \]

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

---

1. Ateshian, Gerard A; Shim, Jay J; Maas, Steve A; Weiss, Jeffrey A. "Finite Element Framework for Computational Fluid Dynamics in FEBio." J Biomech Eng, vol. 140 (2018). ↩

2. Bonet, Javier; Wood, Richard D.. "Nonlinear continuum mechanics for finite element analysis." Cambridge University Press (1997). ↩↩

3. Jansen, Kenneth E; Whiting, Christian H; Hulbert, Gregory M. "A generalized-\(\alpha\) method for integrating the filtered {Navier}--{Stokes} equations with a stabilized finite element method." Comput. Methods Appl. Mech. Engrg., vol. 190, pp. 305--319 (2000). ↩↩

4. Reddy, J. N.; Gartling, David K.. "The finite element method in heat transfer and fluid dynamics." CRC Press, pp. 469 p. (2001). ↩

5. Bazilevs, Y; Gohean, JR; Hughes, TJR; Moser, RD; Zhang, Y. "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. Engrg., vol. 198, pp. 3534--3550 (2009). ↩

6. Esmaily Moghadam, Mahdi; Bazilevs, Yuri; Hsia, Tain-Yen; Vignon-Clementel, Irene E; Marsden, Alison L. "A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations." Comput. Mech., vol. 48, pp. 277--291 (2011). ↩

7. Vignon-Clementel, Irene E; Figueroa, C Alberto; Jansen, Kenneth E; Taylor, Charles A. "Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries." Comput. Methods Appl. Mech. Engrg., vol. 195, pp. 3776--3796 (2006). ↩