Computational Fluid Dynamics (CFD) in FEBio
Introduction
Modeling the behavior of fluids is an important topic in biomechanics, particularly when studying blood flow through the cardiovascular system and medical devices, cerebrospinal fluid flow, airflow through the respiratory system, biotribology by fluid film lubrication, and flow through microfluidic biomedical devices. An example is shown below in Figure 1 where blood flow is modeled through a carotid artery bifurcation.
The CFD finite element solver in FEBio takes a different approach compared to other CFD codes. This tutorial gives an overview on the fundamentals of the CFD solver in FEBio and provides the user practical information on how to use it correctly. Please see the tutorial “Creating Your First CFD Model in FEBioStudio” to create a simple CFD model.
Degrees of Freedom and Boundary Conditions
FEBio’s fluid solver differs from most computational fluid dynamics programs by its use of fluid dilatation $e^f$, instead of fluid pressure $p$, as a nodal degree of freedom. In the fluid module, these two variables are related by a simple constitutive relation, $p=-K e^f$, where $K$ is the fluid’s bulk modulus. The reason for selecting $e^f$ instead of $p$ is that the dilatation is a kinematic measure, just like the fluid velocity $\mathbf{w}$ (the other nodal degree of freedom), whereas the pressure is a function of state (it needs to be formulated as a function of the deformation, just like the viscous fluid stress $\boldsymbol{\tau}$).
The fluid Cauchy stress is given by $\boldsymbol{\sigma}^{f}=-p\mathbf{I}+\boldsymbol{\tau}$. The viscous stress $\boldsymbol{\tau}$ depends on the fluid rate of deformation tensor $\mathbf{D}^{f}$ according to a user-selected constitutive relation, such as a Newtonian fluid. FEBio supports several options for viscous stress constitutive relations, including Newtonian fluid and several non-Newtonian fluids (e.g. Carreau, Carreau-Yasuda, Powell-Eyring, Cross). If a desired constitutive model is not available, users can add their own using a plugin.
In a fluid analysis, the user may prescribe any of the components of the fluid velocity $\mathbf{w}$ as an essential boundary condition, or the corresponding component of the viscous traction $\mathbf{t}^{\tau}=\boldsymbol{\tau}\cdot\mathbf{n}$ as a natural boundary condition ($\mathbf{n}$ being the unit outward normal to the boundary surface). When neither a component of $\mathbf{w}$ nor the corresponding component of $\mathbf{t}^{\tau}$ is prescribed explicitly, the latter is naturally assumed to be zero. Similarly, the user may either prescribe the dilatation $e^{f}$ as an essential boundary condition (which is equivalent to prescribing $p=-K e^{f}$), or $w_{n}=\mathbf{w}\cdot\mathbf{n}$ as a natural boundary condition. When neither is prescribed explicitly, it is naturally assumed that $w_{n}=0$. Natural boundary conditions are useful for creating symmetry planes in fluid analyses, where the normal fluid velocity and the viscous traction are zero: On those symmetry planes, no boundary conditions should be specified. An example where such a boundary condition is useful is shown in Figure 3, which is an axisymmetric wedge model representing a circular cylinder; when the flow is in the cylinder’s axial $Z-$direction, there is no need to prescribe a zero-flow condition normal to the side surfaces of the wedge.
It is noteworthy that the fluid formulation in FEBio allows the prescription of the nodal value of $\mathbf{w}$ as an essential boundary condition, and the surface value of w_n as a natural boundary condition. On boundaries where the fluid velocity is completely prescribed, i.e., when normal and tangential components are known (and not zero), prescribing $\mathbf{w}$ and $w_{n}$ instead of just $\mathbf{w}$ produces the best computational outcome. In those cases, the user should use the fluid velocity surface load, which combines both boundary conditions. If the user only wants to prescribe a known normal velocity $w_{n}$ and leave the tangential velocity unspecified, use the fluid normal velocity surface load; this surface load also provides the option of setting the tangential fluid velocity to zero. Finally, if the user wants the velocity to remain normal to the selected surface (e.g., an inlet or outlet surface on which $e^{f}$ is prescribed or fixed), but does not know the velocity magnitude a priori, use the normal fluid velocity constraint to enforce this condition (unless that surface is aligned with one of the coordinate planes, in which case the appropriate components of $\mathbf{w}$ may be set to zero).
A viscous fluid satisfies the no-slip condition on a physical boundary surface. This means that $\mathbf{w}=\mathbf{0}$ on no-slip boundaries. This boundary condition can be enforced simply by fixing the three components of $\mathbf{w}$ and not specifying any value for $e^{f}$ on that surface (which naturally enforces $w_{n}=0$).
Specialized Surface Loads
Computational fluid dynamics analyses are often performed over a mesh that truncates the real fluid domain. Thus, fluid may enter the finite element domain across some upstream inlet boundary and leave the domain across some downstream outlet boundary. The exact flow conditions on these boundaries may not be known, thus they have to be approximated using best guesses. Arguably, the viscous traction $\mathbf{t}^{\tau}$ is the least intuitive boundary condition to guess. On an outlet boundary, it is necessary to prescribe or fix the dilatation $e^{f}$ to overcome the natural boundary condition $w_{n}=0$. However, in most cases the tangential components of $\mathbf{w}$ on that boundary are not known, therefore they cannot be fixed or prescribed, leading to the natural condition $\mathbf{t}^{\tau}=\mathbf{0}$. This natural condition may work fine in some cases, but it will often fail numerically when vortices being shed inside the finite element domain try to cross the outlet boundary. In those cases, consider using the fluid backflow stabilization and fluid tangential stabilization surface loads, which prescribe a velocity-dependent viscous traction $\mathbf{t}^{\tau}$. Figure 4 shows the flow past a block example, where fluid backflow stabilization and fluid tangential stabilization prescribed at the (right-most) outlet boundary are necessary for the problem to converge.
Similarly, on an inlet boundary that shares an edge with a no-slip surface, prescribing the fluid velocity $\mathbf{w}$ on the inlet boundary may not accurately reproduce the velocity profile in the boundary layer that would be produced on the adjoining no-slip surface. This potential mismatch could result in numerical instabilities; in such cases, prescribing a dilatation $e^f$ on the shared edge may mitigate this issue. Figure 5 shows where the dilatation was prescribed to avoid this issue at the inlet rim for the bifurcated carotid artery.
Biased Meshes for Boundary Layers
In fluid analyses, boundary layers form in the vicinity of no-slip boundaries, where the velocity magnitude varies rapidly with the distance from a no-slip surface. The thickness of boundary layers decreases with increasing Reynolds number. To capture these boundary layers, and resulting wall shear stresses, accurately in a fluid finite element analysis, it is necessary to refine the mesh to produce thinner elements closer to the no-slip boundaries. This is done most conveniently by using mesh biasing tools, available in various meshing software programs, or post-hoc boundary-layer meshing tool that modify an existing mesh. Both of these options are available in FEBio Studio. Boundary layer meshing should always be used in fluid analyses in FEBio. Whereas other finite element codes employ numerical stabilization techniques that partially alleviate the need for biased meshes, FEBio does not employ such techniques. Therefore, using uniform finite element meshes may fail to produce numerical convergence, even in some seemingly simple problems.
One particular tool in FEBio Studio that is useful for creating biased meshes, is the Boundary Layer tool. An example of using this tool is shown in Figure 6, where the no-slip surface has to be pre-selected, before the tool is applied. While the tool can create a biased mesh for existing (linear) tetrahedral, pentahedral, and hexahedral meshes, doing so may greatly increase the number of nodes and elements, particularly for tetrahedral meshes.
Computational Efficiency: Broyden’s Method
Fluid finite element analyses produce a non-symmetric stiffness matrix and the resulting system of equations may be efficiently solved using Broyden’s quasi-Newton method, where an approximation to the matrix inverse is produced at each iteration after the first full-Newton iteration of a time step. In fluid analyses, it is often possible to achieve convergence at each time step without having to perform additional full-Newton updates. Therefore, it is recommended to set max_updates to a large number (e.g., 50), along with setting diverge_reform to 0 (false), which leads to a more efficient solution scheme. A further advantage arises in fluid analyses since the mesh remains invariant over time: it is often possible to continue using Broyden updates even across time steps, without performing a full-Newton iteration at the start of that step, by setting reform_each_time_step to 0 (false). In that case, the solver will continue using Broyden updates up to the value of max_updates, before performing a full Newton update.
Dynamic versus Steady-State Analyses
FEBio runs fluid analyses in dynamic mode by default, because many fluid flow problems do not have a steady-state solution, since the governing equations are nonlinear. For example, many flows may exhibit vortex shedding as fluid flows across an obstacle, such as a flow constriction or a solid body. Even if some form of periodicity emerges under specific flow conditions, as in the von Karman vortex street, the corresponding solution is not in steady state. However, for some low Reynolds number laminar flows, a steady-state response may exist and FEBio will attempt to find that solution when using analysis type “steady-state.” The user should be aware that, under steady-state analyses, the solution may not necessarily converge as efficiently to the final solution as would a dynamic analysis that allows the solution to reach steady state after a sufficiently long time. The best choice of analysis type for steady-state problems may need to be determined by trial and error.
Isothermal Compressible Flow versus Acoustics
The governing equations for fluid flow in FEBio represent isothermal conditions, thus eliminating temperature as a variable in the numerical analysis. Since the fluid solver accommodates some measure of fluid compressibility (via the fluid dilatation variable $e^f$), it is possible to solve for pressure wave propagation in fluids using this formulation. However, from a theoretical perspective, it should be recognized that the speed of sound (i.e., the wave speed) in isothermal (constant temperature) flow is different from the speed of sound in isentropic (constant entropy) flow. The field of acoustics typically examines pressure wave propagations under isentropic conditions, which are found to be more consistent with experimental measurements of the speed of sound in air. Therefore, the FEBio fluid solver may not produce realistic wave speeds when simulating acoustics. A fluid solver for isentropic flow may be developed for FEBio in the future.