9.8 Guidelines for Fluid Analyses¶
Currently, the fluid solvers of FEBio can be used to solve for isothermal fluid mechanics problems (“fluid” solver), isothermal fluid-structure interaction problems (“fluid-FSI” solver, Section Fluid-Structure Interactions) and isothermal fluid mechanics with solute (mass) transport (“fluid-solutes” solver). When solutes are included in an analysis, their treatment is very similar to the “multiphasic” solver described in Section Triphasic Multiphasic Materials. The fluid pressure \(p\) and solute concentrations \(c^{\iota}\) in a fluid-solutes domain are thus related to the effective fluid pressure \(\tilde{p}\) and effective solute concentrations \(\tilde{c}^{\iota}\) as per eq.\eqref{eq:mp-effective-p-c}, to account for the possibility that the solutes are electrically charged. The theoretical treatment of solutes in the “fluid-solutes” solver follows that of the “multiphasic” solver, with the following simplifications: (a) The effective solubility \(\hat{\kappa}^{\iota}\) appearing in the expression for the partition coefficient \(\tilde{\kappa}^{\iota}\) in eq.\eqref{eq:mp-effective-p-c} only embodies non-ideal physico-chemical behavior (whereas \(\hat{\kappa}^{\iota}\) in a multiphasic mixture with a porous solid matrix can also embody the steric volume exclusion of solutes from the pore space); therefore, one should always set \(\hat{\kappa}^{\iota}=1\) when modeling ideal behavior in a “fluid-solutes” analysis. (b) Since there is no porous solid matrix that can hinder solute transport, the diffusivity tensor \(\mathbf{d}^{\iota}\) of a solute in a “fluid-solutes” mixture is always equal to \(d^{\iota}_{0}\mathbf{I}\), where \(d^{\iota}_{0}\) is the solute diffusivity in free solution. Thus, even though FEBio requires users to specify a value for \(\mathbf{d}^{\iota}\), that value is essentially ignored in the “fluid-solutes” solver code, where it is replaced by \(d^{\iota}_{0}\mathbf{I}\). If users would like to neglect osmotic pressure in a “fluid-solutes” analysis, the osmotic coefficient \(\Phi\) should be set to zero, producing \(\tilde{p}=p\). Moreover, in the absence of solutes (\(c^{\iota}=0\) for all \(\iota\), as in the “fluid” and “fluid-FSI” solvers) this relation holds automatically. Therefore, when dealing with the “fluid” and “fluid-FSI” solvers, we do not need to add the “effective” modifier when mentioning the fluid pressure.
There are two broad categories of materials available for “fluid-FSI” analyses: (1) fluid-FSI” (Section General Fluid FSI Materials), and (2) biphasic-FSI (Section General Fluid BFSI Materials). A fluid-FSI material is used to describe a fluid domain whose mesh is deformable; though a solid material (typically neo-Hookean) must be assigned to the mesh of a fluid-FSI material, its mass density is ignored and it is the user's responsibility to ascribe a negligible elasticity to this solid domain, whose sole purpose is to regularize the mesh deformation. A biphasic-FSI material considers that the deformable domain is a biphasic material (see Section Biphasic Materials), with a porous-permeable solid domain that has non-negligible mass density, solid volume fraction, solid matrix stiffness, and hydraulic permeability. The interstitial fluid of such a domain may also be ascribed a viscosity (Newtonian or non-Newtonian).
Degrees of Freedom and Boundary Conditions¶
The nodal degrees of freedom for “fluid” analyses are the nodal fluid velocity \(\mathbf{v}^{f}\) and nodal fluid dilatation \(e^{f}\). For a “fluid-FSI” analysis the nodal fluid velocity \(\mathbf{v}^{f}\) is replaced with the relative fluid flux \(\mathbf{w}=\varphi^{f}\left(\mathbf{v}^{f}-\mathbf{v}^{s}\right)\) and supplemented with the nodal solid displacement, \(\mathbf{u}\), whose material time derivative equals the solid velocity \(\mathbf{v}^{s}\). Here, \(\varphi^{f}\) represents the fluid volume fraction (porosity) of the fluid-FSI domain, which is prescribed by the user in a biphasic-FSI material (via the complementary solid volume fraction) or automatically set to \(\varphi^{f}=1\) in a fluid-FSI material. Thus, the degrees of freedom of a fluid-FSI analysis reduce to those of a fluid (or fluid-solutes) analysis in the limit when \(\varphi^{f}\to1\) and \(\mathbf{v}^{s}\to\mathbf{0}\). For “fluid-solutes” analyses, the nodal degrees of freedom include the effective solute concentrations \(\tilde{c}^{\iota}\); as for “multiphasic” analyses, it is assumed that solutes occupy a negligible volume fraction in the mixture (thus \(\varphi^{f}=1\) for “fluid-solutes” analyses).
The fluid Cauchy stress is given by \(\boldsymbol{\sigma}^{f}=-p\mathbf{I}+\boldsymbol{\tau}\), where \(p\) is the actual (not the effective) fluid pressure, as per eq.\eqref{eq:mp-effective-p-c}. 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 Newtonian fluid.
In any of the fluid analysis types, the user may prescribe any of the components of \(\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 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 prescribe \(e^{f}\) (or \(\tilde{p}\left(e^{f}\right)\), see Section Fluid Pressure) as an essential boundary condition, 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.
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 (Section Fluid Velocity). 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 (Section Fluid Normal Velocity); 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 (Section Normal Fluid Velocity) to enforce this condition.
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}\) (Section Prescribed Displacement) and not specifying any value for \(e^{f}\) on that surface (which naturally enforces \(w_{n}=0\)).
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 (Section Fluid Backflow Stabilization) and fluid tangential stabilization (Section Fluid Tangential Stabilization) surface loads, which prescribe a velocity-dependent viscous traction \(\mathbf{t}^{\tau}\).
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.
Biased Meshes for Boundary Layers¶
In fluid analyses, boundary layers will 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 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 FEBioStudio. 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. Users should keep in mind that biased meshes don't necessarily require more nodes and elements, therefore there is no automatic computational cost for using biased meshes.
Computational Efficiency: Broyden's Method¶
Fluid 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 (Section Control Section) 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” (Section Control Section). 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.
Fluid-Structure Interactions¶
Fluid-structure interactions (FSI) may be modeled in FEBio using the fluid-FSI module (Section Module Section). In an FSI analysis the finite element mesh defining the fluid domain is deformable, so that the shape of the fluid domain does not have to remain fixed over time. The fluid flows through the deforming mesh, just as it does in standard computational fluid dynamics (CFD) analyses with fixed meshes. The deformation of the fluid domain may be controlled by prescribing suitable boundary conditions (so-called moving boundary problems) and by interactions with solid domains. For example, fluid may flow inside a flexible tube and the fluid pressure may cause the tube to expand or contract; the tube itself may be pinched at some location, reducing the cross-section through which fluid flows. The tube wall is a deformable solid domain which may be modeled in FEBio using any of the solid material models described in Sections Elastic Solids through Multigeneration Solids, or it may be modeled as a porous deformable domain containing a viscous fluid, in which case one should use the biphasic-FSI material for the tube wall (Section General Fluid BFSI Materials). The inner tube (the fluid domain) should be modeled as a fluid-FSI material (Section General Fluid FSI Materials).
In FEBio, a fluid-FSI material is modeled as a specialized fluid-solid mixture where the solid matrix density is zero and its stiffness should be set to a negligible value. In contrast, a biphasic-FSI material is a full-fledged fluid-solid mixture, where the solid may have mass and exhibit any desired nonlinear elastic or inelastic response. Both FSI materials require the user to define a fluid material from the same list as those available for CFD analyses (Section Fluids). It also requires the user to define a solid material to regularize the mesh deformation, which may be selected from Sections Elastic Solids through Multigeneration Solids for a biphasic-FSI material. In contrast, the solid material of a fluid-FSI material is massless (the density value entered by the user is ignored and reset internally to zero). Among all the choices available in Section Elastic Solids, it is recommended to use the neo-Hookean elastic solid (Section Neo Hookean), setting Poisson's ratio \(\nu\) to zero and Young's modulus \(E\) to a very low value (a few orders of magnitude smaller than the moduli of surrounding solid domains).
The degrees of freedom at each node of a FSI domain consist of the components of the solid displacement \(\mathbf{u}\) (which describe the mesh deformation), the fluid velocity \(\mathbf{w}\) relative to the solid, and the fluid dilatation \(e^{f}\). Essential boundary conditions may be prescribed on any of these degrees of freedom. For example, in a moving boundary problem, the motion of a boundary may be prescribed by setting the relevant components of \(\mathbf{u}\). A no-slip boundary between an FSI material and a solid domain may be defined by simply setting \(\mathbf{w}=\mathbf{0}\), regardless of the motion of that boundary. In general, all boundary conditions described for CFD analyses in Section Degrees of Freedom and Boundary Conditions may be applied identically to FSI analyses.
Currently in FEBio, the mesh of a fluid/biphasic-FSI material domain must be continuous with the mesh of surrounding solid domains. The surface between domains, which consists of the shared faces of adjoining finite elements from each domain, is called the FSI interface. Because the mesh is continuous, the solid displacement and the solid component of the traction are automatically continuous across FSI interfaces; when fluid is present on both sides of the interface, the relative fluid flux is also automatically continuous. In other words, when two FSI domains are adjoining (e.g., fluid-FSI with fluid-FSI, or fluid-FSI with biphasic-FSI), all boundary conditions are automatically satisfied across those interfaces. However, when an FSI domain interfaces with a solid domain, the fluid component of the traction in the fluid-FSI and biphasic-FSI domains must be explicitly transferred to the solid domain across that interface. This is done by prescribing a fluid-FSI traction or a biphasic-FSI traction on that interface (Section Fluid FSI Traction). This surface load does not require any user-defined parameters; its sole purpose is to identify those FSI interfaces where the fluid traction must be properly transferred to the surrounding solid domain. Omitting this boundary condition means that the fluid traction (i.e., the fluid pressure and the normal and tangential viscous components of the fluid stress) have no effect on this interface; the interface may still deform in response to the motion of the solid domain. For example, when a fluid interacts with a rigid solid domain whose motion is entirely prescribed, it is not necessary to apply a fluid-FSI traction at the interface between the fluid and rigid solid domains.
A fluid-FSI traction is also required on free fluid surfaces, such as the surface of a fluid in open channel flow, in order to properly capture the motion of that surface (such as waves). This condition is required since the net traction \(\mathbf{t}\) on a free surface is zero. As the fluid-FSI domain is modeled as a specialized solid-fluid mixture, the net traction consists of the sum of solid and fluid tractions, \(\mathbf{t}=\mathbf{t}^{s}+\mathbf{t}^{f}\). Setting \(\mathbf{t}\) to zero produces a traction \(\mathbf{t}^{s}=-\mathbf{t}^{f}\) on the solid mesh, allowing it to deform to the proper shape.
Here are a few cautionary notes to keep in mind when running FSI analyses: (1) Depending on the problem being analyzed, fluid-structure interactions may cause significant mesh deformations that lead to mesh distortion. Currently, adaptive meshing is not yet available to relieve this mesh distortion, therefore analyses and meshes need to be designed to minimize this effect or to terminate before the solution becomes too corrupted. (2) It is possible to pinch a solid domain that surrounds a fluid domain until the flow cross-section has reduced considerably; however it is not possible to completely shut off the flow, as this would require reducing the volume of some finite elements to zero. (3) While it is possible to model the deformation of free fluid surfaces, such as the formation of surface waves in open channel flow or the deformation of fluid blobs floating in a gravity-free environment, it is not possible to model the separation of such fluid domains into sub-domains such as spraying fluid droplets; this means that analyses where such phenomena are expected to happen will likely fail. (4) Whereas standard CFD analyses may be jump-started by prescribing a relatively large fluid velocity at the very first time step of an analysis, these types of boundary conditions may lead to dynamic oscillations and potential instabilities in an FSI analysis; therefore, users must be more considerate when prescribing boundary conditions for FSI analyses.
Fluid-Solutes Analyses¶
Fluid-solutes analyses introduce two types of molar fluxes into an analysis: (1) The diffusive molar flux \(\mathbf{j}^{\iota}_{d}=-\tilde{\kappa}^{\iota}d^{\iota}_{0}\grad\tilde{c}^{\iota}\) and (2) the sedimentation flux \(\mathbf{j}^{\iota}_{b}=s^{\iota}c^{\iota}\mathbf{b}^{\iota}\), where \(s^{\iota}\equiv d^{\iota}_{0}\frac{M^{\iota}}{R\theta}\) is the sedimentation coefficient and \(\mathbf{b}^{\iota}\) is the body force acting on solute \(\iota\). These molar fluxes represent the flux of the solute relative to the solvent. The natural boundary condition for solutes is \(\tilde{j}^{\iota}_{n}=0\), where
is the normal component of the effective molar flux \(\mathbf{j}^{\iota}+\sum_{\gamma}z^{\gamma}\mathbf{j}^{\gamma}\) on a boundary, with \(\mathbf{j}^{\iota}=\mathbf{j}^{\iota}_{d}+\mathbf{j}^{\iota}_{b}\) representing the sum of diffusive and sedimentation fluxes. The effective molar flux is constructed in this manner to enforce the electroneutrality condition,
as also done in “multiphasic” analyses (Section Prescribed Effective Solute Flux).
Zero natural flux \(\tilde{j}^{\iota}_{n}\) implies that the solute moves at the velocity of the solvent (the fluid). Therefore, if the fluid satisfies \(v^{f}_{n}=0\) on that boundary, the solute will not transport across it. But if the fluid is allowed to transport across the boundary (e.g., using a prescribed fluid pressure), the solute can also transport across the boundary via convection. Therefore, for open flows or flows through conduits, where the fluid pressure is prescribed (e.g., set to zero) at an outlet boundary, no specific boundary condition needs to be prescribed on the solute in order for the latter to convect freely with the solvent at that outlet boundary. However, if the solvent exhibits back flow on that boundary in a manner that requires the application of the fluid backflow stabilization (Section Fluid Backflow Stabilization) and fluid tangential stabilization (Section Fluid Tangential Stabilization) surfaces loads, then consider also using the solute backflow stabilization surface load (Section Solute Backflow Stabilization) on that boundary, to produce more physically-realistic solute flows on such truncated domains, and improve numerical convergence.
The fluid pressure \(p\) in a fluid-solutes analysis is given by eq.\eqref{eq:mp-effective-p-c}. To prescribe the fluid pressure \(p\) on the boundary of a fluid-solutes domain, use the fluid pressure boundary condition specified in Section Fluid Pressure BC.
By default, the fluid-solutes solver employs the staggered equation_scheme format for the stiffness matrix, and solves for the nodal fluid velocity \(\mathbf{v}^{f}\), nodal dilatation \(e^{f}\) and nodal effective solute concentrations \(\tilde{c}^{\iota}\) using a coupled solve_strategy. This default settings works well when the Peclet number (the ratio of solute convective velocity to solute diffusive velocity) is relatively low (e.g., less than 100). For these problems, using the default value of ctol=0.01 is generally appropriate. However, if users notice significant convergence difficulties on the solute concentration as the Peclet number approaches 100 or greater, consider setting ctol=0. When solving strongly convective solute transport problems with high Peclet numbers (e.g., greater than 1000), diffusion will be negligible compared to convection. In that case, users should switch to equation_scheme=block and solve_strategy=sequential and set ctol=0 to increase the chances of achieving a converged numerical solution. The rationale for these choices is explained in the FEBio journal article associated with this fluid-solutes solver 1.
As done in the multiphasic solver, the solute volume fraction is considered negligible in fluid mixtures modeled in the fluid-solutes solver. Therefore, the volume of a fluid mixture will not change with increasing solute concentrations. When a solute is convected by the solvent, it means that the solvent exerts a frictional force on the solute, whose magnitude is inversely proportional to the solute diffusivity \(d^{\iota}_{0}\). In the fluid-solutes solver, the acceleration of solutes is considered negligible, implying that the velocity of the solute in the solvent reaches its terminal velocity instantaneously (thus, inertial forces are negligible compared to the drag force between solute and solvent).
According to Newton's third law stating that every action force has an equal and opposite reaction force, the drag force exerted by the solvent on the solute is balanced by an equal and opposite reaction force exerted by the solute on the solvent. However, by default, this reaction force on the solvent is turned off in the fluid-solutes solver. In other words, the motion of the solvent affects that of the solute, but the motion of the solute does not affect the solvent. This default setting is adopted because it represent the common form of diffusion-reaction equations modeled in most of the mass transport literature. Therefore, when modeling solute diffusion problems in a stationary solvent (e.g., in the absence of pressure or velocity-driven solvent flows), the diffusive flux of a solute will not drag the solvent, keeping it in a stationary state. When the solvent motion is driven by solute concentration gradients, we call it osmosis. There is no common term describing solvent motion caused by solute sedimentation (motion of solute under the action of a body force), so we use the term osmosis to describe this phenomenon as well. Users can override the default setting of the fluid-solutes solver by checking the include osmosis box in the fluid-solutes material definition. When this box is checked, the solvent will have a force exerted on it by the solute (either due to diffusion or sedimentation, or both), causing it to move.
High-Peclet number flows are notoriously difficult to solve numerically. In FEBio we have reported valid solutions for a Peclet number as high as \(10^{11}\) 1, thanks to the specific finite element formulation employed in the fluid-solutes solver, and the use of the sequential solver when needed. Nevertheless, the fluid-solutes solver is not fool-proof. In some problems, spurious increases in solute concentration can occur, seemingly randomly in the mesh, even when using a highly regular mesh. In many of these cases, this type of numerical artifact may go away by simply refining the mesh. However, users should be aware that this artifact may emerge 'silently', meaning that it does not necessarily produce slower numerical convergence that could hint at this problem.
The explanation for this behavior is that the numerical solutions to the effective solute concentrations at large Peclet numbers are effectively solutions to the partial differential equation \(D^{f}\left(J^{f}c^{\iota}\right)/Dt=0\), as reviewed in the Theory Manual. This integration is performed based on the current solution for \(J^{f}=1+e^{f}\) and \(\mathbf{v}^{f}\) as well as the structure of the mesh. Any ill-conditioning in the mesh (such as a mesh which is too coarse) may introduce a numerical round-off error in the calculation of \(\tilde{c}^{\iota}\), and there is no other sufficiently dominant term in the solute mass balance equation that can dampen or correct for such numerical round-offs. Thus, they may grow over time. Mesh refinements are required when encountering these circumstances.
Other boundary conditions for the solute (such as the essential boundary condition \(\tilde{c}^{\iota}=0\)) may also improve predictions when using relatively coarse meshes, preventing the occurrence of spurious rises which often initiate on boundaries. Unfortunately, in some problems, this option forces users to adopt a non-physical boundary condition which may produce unacceptable results. This workaround may not be needed when using a sufficiently refined mesh, though finer meshes typically require longer solution times. Also note that highly irregular meshes near no-slip boundaries may cause spurious rises in solute concentration in high Peclet number flows. Therefore, it is recommended to use suitable meshing techniques for generating boundary layer meshes, such as stacking pentahedral elements atop tetrahedral elements when using unstructured tetrahedral meshes. For diffusion-dominated or reaction-dominated flows however, these numerical concerns may not arise. Nevertheless, a thorough understanding of the underlying physics of these problems is essential for properly prescribing boundary conditions, and the number and size of time increments in a finite element analysis.
The default settings for the fluid-solutes solver force it to reform the stiffness matrix at the start of each time step. This is a standard requirement in most of the FEBio solvers (including solid and structural mechanics, and standard biphasic and multiphasic domains), though it had previously been demonstrated that the fluid solver could avoid reforming the stiffness matrix even after achieving convergence at a given time step, producing a very efficient computational scheme. Unfortunately, this increased efficiency is not always an option with the fluid-solutes solver, because the stiffness matrix contains terms that depend significantly on the solute concentrations. As these concentrations evolve over time and space, it becomes essential to update the stiffness matrix to reflect those changes. Therefore, consider using the direct sparse solver from the Intel MKL library (linear_solver=Mkl dss) to improve computational efficiency during matrix factorization, in comparison to the Pardiso solver from the same library. Moreover, in applications where the initial solute concentration is non-zero, and non-negligible compared to the prescribed change in concentration on an inlet boundary, it may be possible to solve a fluid-solutes problem more efficiently without reforming the stiffness matrix at each time step.