A time-consistent stabilized finite element method for fluids with ...

Studies of the human cardiovascular system have greatly benefited from the advances in computational fluid dynamics (CFD) since the end of the twentieth century1,2,3. Among the numerical methods for solving the Navier–Stokes equations4,5,6, the finite element method has found popularity for cardiovascular CFD simulations because it is a convenient framework for dealing with complex geometries and modeling fluid–structure interaction7,8,9,10. Together with advancements in clinical imaging techniques11,12,13,14, CFD simulations using the finite element method have taken a significant part in in-vitro studies, clinical diagnosis, and surgical planning for cardiovascular diseases15,16,17,18,19,20.

The finite element method for solving the unsteady Navier–Stokes equations relies on an upwind term that adds an artificial diffusion along the stream-wise directions, weighted by a stabilization parameter \(\tau\), to prevent nonphysical oscillations inherent to the Galerkin method in strongly convective regimes21,22,23,24. One of the most commonly used formulations of \(\tau\) has been the one proposed in the streamline upwind Petrov–Galerkin formulations (SUPG)21, \(\tau _{\text{SUPG}}\), which is also adopted in others such as the residual-based variational multiscale (RBVMS) formulation25,26.

Traditionally, the steady form of \(\tau _{\text{SUPG}}\) is derived from a 1D steady advection-diffusion model problem, such that the added diffusion to the Galerkin’s formulation is just enough to recover the exact solution, thereby eliminating numerical oscillations at higher element Peclet numbers21,22,27,28. While this steady form of \(\tau _{\text{SUPG}}\) works well for strongly advective steady-state flows, it is not readily applicable to time-varying flows, as it exhibits poor convergence behavior, particularly at smaller time step sizes (\(\Delta t\)). The traditional strategy to overcome this issue has been adding a \(\Delta t\)-dependent term to the definition of \(\tau _{\text{SUPG}}\)22,28,29. The added \(\Delta t\)-dependent term, which is based on the discrete approximation of the inverse of the strong differential operator, dominates the contributions associated with the advection and diffusion terms at small time step sizes. This design, which has found widespread use for its excellent convergence characteristics, produces a strong solution dependency on the time step size, such that the solution becomes less accurate as the time step size is reduced toward zero30,31,32 (also see the “Results” section). As a consequence, this time step size-dependent design of \(\tau _{\text{SUPG}}\) produces an inconsistent method with regard to the time step size. This inconsistency issue is particularly exacerbated in a subset of cardiovascular simulations that demand a small time step size for their multiscale behaviors, such as those involving lumped parameter network modeling33,34,35,36. In fact, some of the most popular packages for cardiovascular simulation37,38 are based on a stabilized finite element formulation that also suffers from the described issue, thereby motivating the present study.

There has been some effort in the past to overcome the inconsistency issue associated with this design of \(\tau _{\text{SUPG}}\)30,39,40. Tezduyar and others introduced an element-vector-based \(\tau _{\text{SUPG}}\) that uses the relative elemental magnitude of terms in the weak form of the Navier–Stokes equations32,41. Although time step size dependency is not entirely eliminated with the proposed formulation, a notable reduction in the solution variation with the time step size was reported in turbulent channel flow simulations. Furthermore, this method relies on sets of integrals on each element that must be performed before the evaluation of the discrete form at the Gauss quadrature. Thus, it may not be simple to implement this technique in an existing finite element program that is not already designed based on the element-vector-based method. In another study, R. Codina and others proposed a subscale-tracking approach that solves a time-dependent ordinary differential equation at each Gauss integration point to evolve its stabilization parameter in time31. This method eliminates the time step size dependency for steady-state solutions. The additional computational cost to solve an ordinary differential equation at each Gauss point, however, is non-negligible for this method. Despite the relatively better time step size consistency shown by these methods in the numerical experiments, they have not found widespread implementation in cardiovascular simulations37,42.

In this study, we propose a formulation for \(\tau _{\text{SUPG}}\) that eliminates the solution’s dependency on the time step size. This method is simple to implement in existing CFD solvers that are based on the streamline upwind Petrov–Galerkin and pressure stabilizing Petrov–Galerkin (SUPG/PSPG) method. More specifically, we replace the inverse of the time step size in the definition of the \(\tau _{\text{SUPG}}\) with a measure of the flow frequency. The motivation behind such formulation lies in the spectral formulation of the unsteady Stokes equations where the time step size dependent parameter in the \(\tau _{\text{SUPG}}\) is replaced by the spectral mode number43. The same parameter, when expressed in a spatio-temporal setting, inspires the use of a flow-dependent time scale in the \(\tau _{\text{SUPG}}\) definition, hence motivating the present design.

The idea of replacing the time step size in \(\tau _{\text{SUPG}}\) with an acceleration-to-velocity ratio has been proposed independently earlier by Evans and others44 to simulate turbulent flows. Our present study nonetheless distinguishes from the previous study both in the formulation and the application. First, the previous study adopts an elemental measure of velocity and acceleration, while we use the global-averaged values in defining \(\tau _{\text{SUPG}}\). As we discuss later in the formulation section, this choice was made for stability reasons. While a local measure produces stable results when used in conjunction with an explicit time integration, it creates instabilities in implicit solvers, which is the concern of the present study. Second, the present study investigates the behavior of the new \(\tau _{\text{SUPG}}\) using a range of canonical and realistic anatomical cases, thereby evaluating its potential for cardiovascular simulations.

The article is organized as follows: We first present the formulation of the stabilized finite element method for the Navier–Stokes equations, the motivation behind the proposed \(\tau _{\text{SUPG}}\), and its formulation. We present four cases to compare the present formulation against the conventional one: a pipe flow with steady boundary conditions, a time-periodic blood flow in a complex cardiovascular geometry (modified Blalock–Taussig shunt), a two-dimensional external flow over a square, and a two-dimensional flow over a square with an attached flexing beam. We will also discuss the convergence and computational cost of the present formulation. Lastly, we conclude our study and discuss the future outlook in the broader fields of cardiovascular simulations and stabilized finite element methods.

The Navier–Stokes equations for incompressible flows are stated as

where \(\rho\) is the density, \(\textbf{u}(\textbf{x},t)\) is the velocity, \(\textbf{f}(\textbf{x},t)\) is the external forcing, \(\Omega \times [0,T]\) is the fluid computational spatio-temporal domain, and the stress tensor

where \(p(\textbf{x},t)\) is pressure and \(\mu\) is the dynamic viscosity. The Dirichlet and Neumann boundary conditions are defined as

respectively, where \(\mathrm {\Gamma _g}\) and \(\mathrm {\Gamma _h}\) are subsets of the boundary \(\Gamma\) where the Dirichlet and Neumann boundaries are prescribed, \(\textbf{n}\) is the outward unit normal vector, and \(\textbf{g}\) and \(\textbf{h}\) are the given Dirichlet and Neumann boundary conditions, respectively.

The discrete form of the Navier–Stokes equations we used in this study is stated as finding \(\textbf{u}^h \in S_\textbf{u}^h\) and \(p^h \in S_p^h\) such that for all \(\textbf{w}^h \in V_\textbf{u}^h\) and \(q^h \in V_p^h\),

$$\begin{aligned}{} & {} \int _\Omega \textbf{w}^h\cdot \rho \left( \frac{\partial \textbf{u}^h}{\partial t}+\textbf{u}^h\cdot \varvec{\nabla }\textbf{u}^h-\textbf{f}\right) d\Omega + \int _\Omega \varvec{\varepsilon }\left( \textbf{w}^h\right) :\varvec{\sigma }\left( p^h,\textbf{u}^h\right) d\Omega -\int _{\Gamma ^h}\textbf{w}^h\cdot \textbf{h}^h d\Gamma +\int _\Omega q^h\varvec{\nabla }\cdot \textbf{u}^h d\Omega \nonumber \\{} & {} \quad +\sum _{e=1}^{n_{el}}\int _{\Omega ^e}\left[ \tau _{\text{SUPG}}\left( \textbf{u}^h\cdot \varvec{\nabla }\textbf{w}^h+\frac{1}{\rho }\varvec{\nabla } q^h\right) \cdot \mathbf {R_{\text{M}}}^h + \rho \nu _{\text{C}} \nabla \cdot \textbf{w}^h R_{\text{C}}^h\right] d\Omega = 0. \end{aligned}$$

(4)

In the above problem statement, \(S_\textbf{u}^h\) and \(S_p^h\) are the discrete solution spaces for the velocity and pressure, respectively, and \(V_\textbf{u}^h\) and \(V_p^h\) are the finite-dimensional test function spaces for the velocity and pressure, respectively.

In Eq. (4), the terms directly obtained from Eq. (1) are supplemented with three elemental stabilization terms. The two terms multiplied by \(\tau _{\text{SUPG}}\) are the conventional SUPG and PSPG stabilizations to ensure stability in strongly convective flows and allow for equal order interpolation functions for velocity and pressure, respectively21,22. The term involving \(\nu _{\text{C}}\) comes from the residual-based variational multiscale methods (VMS)45,46,47. While there are some variations in defining these terms, we consider

and

as the common conventional definition of these stabilization parameters for later comparisons28,45,48. In Eq. (5), \(\nu =\mu /\rho\) is the kinematic viscosity, \(\Delta {t}\) is the time step size, \(\varvec{\xi }\) is the covariant tensor obtained from the mapping of the physical-parent elements, and \(C_{\text{I}}\) is a shape-function-dependent constant, which is 3 in our study.

The inconsistency of the above-mentioned stabilized formulation is caused by \(\Delta t\) in \(\tau _{\text{SUPG}}\). In a steady flow, in which the solution should be independent of \(\Delta t\), \(\tau _{\text{SUPG}}\) and thus the overall added diffusion will change with \(\Delta t\), creating a time step size dependent solution. In the later sections, we will demonstrate that this inconsistency issue is not unique to steady-state flows and also occurs for unsteady flows.

In an earlier study, we introduced a pressure-stabilized technique for solving the unsteady Stokes equations expressed in the frequency domain rather than the time domain43. The resulting complex-valued stabilization parameter was derived systematically by taking the divergence of the momentum equation and estimating the Laplacian in the diffusion term using a characteristic element size. The modulus of that PSPG-type stabilization parameter is

where \(\omega\) is the spectral mode appearing as a source term in the frequency formulation of the unsteady Stokes equations. This spectral formulation of \(\tau\) closely resembles the conventional definition of \(\tau _{\text{SUPG}}\) in Eq. (5) if \(2/\Delta {t}\) is replaced by \(\omega\). The \(\textbf{u}^h\cdot \varvec{\xi }\textbf{u}^h\) term does not appear in Eq. (7) as the convective acceleration term is not present in the unsteady Stokes equations. Adding this term into Eq. (7) and incorporating the O(1) constant \(C_{\text{I}}\) results in

which is adopted in this study to replace the conventional definition of \(\tau _{\text{SUPG}}\) from Eq. (5) when solving Eq. (4).

The new definition of \(\tau _{\text{SUPG}}\) in Eq. (8) becomes identical to the traditional formulation (Eq. 5) if \(\omega = 2/\Delta {t}\). This value is close to the largest frequency associated with the time discretization, namely \(\pi /\Delta {t}\) that occurs when the solution oscillates between consecutive time steps. In practice, especially in cardiovascular simulations, the solution is a much smoother function of time and has a frequency content that peaks at a much smaller \(\omega\) than \(\pi /\Delta {t}\). The distinction of these two frequencies inspires the proposed definition of \(\tau _{\text{SUPG}}\).

It is straightforward to evaluate Eq. (8) in a spectral formulation as \(\omega\) is the computed frequency and readily available as an independent parameter. However, its adoption is not straightforward in a traditional time formulation, where \(\omega\) does not appear as an independent parameter. Ideally, \(\omega\) must satisfy several properties. First, it must produce a scheme that remains stable under a variety of conditions. Second, it must be extracted from physical variables, such as the velocity and acceleration, rather than the time step size, so that \(\tau _{\text{SUPG}}\) converges to a unique quantity as the time step size goes to zero. Third, it should be simple to implement and cost-efficient to calculate. Given these criteria, we propose

where \(\Vert \textbf{f}\Vert ^2_{L^2} = \int _\Omega \Vert \textbf{f}\Vert ^2 d\Omega\). This formulation of \(\omega\) is designed to go to zero as the flow reaches a steady state, where \(\frac{\partial {\textbf{u}^h}}{\partial t}\) goes to zero. We will show in the results section that this formulation is consistent in both steady and unsteady flows as \(\Delta t \rightarrow 0\). We will also show that the present formulation is relatively robust even though it increases the computational cost compared to the conventional method.

For moving domain simulations that express Eq. (4) in an arbitrary Eulerian–Lagrangian framework, \(\textbf{u}^h\) in the convective acceleration term is replaced by the fluid velocity relative to moving mesh \(\textbf{u}^h - \hat{\textbf{u}}^h\), where \(\hat{\textbf{u}}^h\) denotes mesh velocity. It is this velocity that is employed in the definition of \(\tau _{\text{SUPG}}\) and also \(\omega\) in Eq. (9), changing it to

where the acceleration term is measured at the mesh node location \(\hat{\textbf{x}}\) and integrals performed over the fluid domain \(\Omega ^f\). By subtracting the mesh velocity from the fluid velocity in Eq. (10), the resulting definition of \(\omega\) will be Galilean invariant. This results in a scheme that produces a unique solution if all velocities were to be measured from a moving inertial reference frame. Later in the result section, we demonstrate the consistency of this method in a moving domain configuration using a fluid structure interaction simulation test case.

Ideally, one would compute velocity and acceleration locally at the Gauss quadrature point when evaluating \(\omega\) so that the solution becomes a function of the local dynamics of the problem. Unfortunately, this choice, which has been successfully employed with explicit time integration in the past44, fails to converge in our implicit formulation. This lack of convergence can be attributed to the velocity appearing in the denominator of Eq. (9), thereby creating widely varying \(\omega\) in regions where flow is temporarily stagnant. This convergence issue is avoided for all cases tested here by using a global measure of velocity and acceleration through integrating their norms over the entire domain as in Eq. (9).

In theory, using a global measure of velocity and acceleration could produce a solution that depends on the domain size. Consider an external flow over an obstacle in which the \(\Vert \frac{\partial {\textbf{u}^h}}{\partial t}\Vert _{L^2}\) term receives non-zero contribution only from regions near the obstacle whereas \(\Vert {\textbf{u}^h}\Vert _{L^2}\) receives a contribution from the entire domain. In this setting, \(\omega\) goes to zero as the size of the computational domain grows, resulting in a domain size-dependent value.

As we will show later in the results section, this domain-size dependency issue does not translate to inconsistency of the formulation in practice. That is because the \(\omega ^2\) term in Eq. (8) is much smaller relative to the sum of the other two terms. In fact, we show that the solution obtained from the proposed formulation is very similar to that of the conventional formulation with a very large \(\Delta t\). Therefore, increasing the domain size will decrease a parameter in the definition of \(\tau _{\text{SUPG}}\) that is already small, thus hardly changing the solution.

Even though the \(\omega ^2\) value in \(\tau _{\text{SUPG}}\) is very small, dropping it from its definition will create convergence issues, as it has been established in the past27,28. The reason that the inclusion of the \(\omega ^2\) term in \(\tau _{\text{SUPG}}\) prevents such scenarios is that a widely varying solution in time will lead to a relatively large \(\omega\). As a result, the contribution of \(\omega\) grows as the solution becomes more unstable, thereby creating a recovery effect that stabilizes the simulation.

In CFD applications, the velocity field may be initialized from zero, thus creating a divide-by-zero operation in the code when evaluating Eq. (9). To avoid such possibilities, we set \(\omega = 2/\Delta {t}\) at the first time step, recovering the conventional definition of \(\tau _{\text {SUPG}}\).

As detailed in a previous publication49, we use the implicit generalized-\(\alpha\) time integration scheme50 in our solver. The use of this time integration scheme significantly simplifies the implementation of the proposed formulation since velocity and acceleration are readily available as discrete state variables. Thus, in our implementation, \(\frac{\partial {\textbf{u}^h}}{\partial t}\) and \(\textbf{u}^h\) in Eq. (9) are explicitly computed from those variables in a single operation. As we will demonstrate in the results, \(\omega\) is a slowly varying parameter, thus we perform this operation only once in each time step using the solution at the previous time step. Given that the generalized-\(\alpha\) method also provides access to variables at the intermediate time points (i.e., \(n+\alpha _m\) and \(n+\alpha _f\) for the acceleration and velocity, respectively), one may elect to use those intermediate variables to update \(\omega\) within each Newton–Raphson iterations. This choice, however, is not adopted in this study as it entails extra computations and has little effect on the overall stability of the solver. The rest of the implementation, including the computation of \(\tau _{\text{SUPG}}\) at the Gauss quadrature points based on the intermediate variables are left unchanged and are identical to the conventional formulation.

The above formulation is implemented in our in-house finite-element solver, multi-physics finite-element solver (MUPFES)35,49,51. A specialized iterative algorithm, preconditioner, and parallelization strategy are employed for an efficient and scalable solution of the linear system of equations49,51,52,53. The solver has been verified54 and extensively employed for cardiovascular modeling in the past55,56,57. This solver is parallelized using a message passing interface (MPI). The workload is parallelized using spatial partitioning by employing ParMETIS library58. All computations are performed on a cluster of AMD Opteron\(^\text{TM}\) 6378 processors that are interconnected via a QDR Infiniband.

At each time step, several Newton–Raphson iterations are performed to ensure the residual falls by over three orders of magnitude. At each Newton–Raphson iteration, a linear system is solved using the generalized minimal residual (GMRES) method59 with a tolerance of \(10^{-2}\).

Four cases are simulated using both the conventional formulation (Eq. 5) and the present formulation (Eq. 8) for \(\tau _{\text {SUPG}}\): a pipe flow, flow in a modified Blalock–Taussig shunt geometry55, a two-dimensional external flow over a square, and a flow over a square with an attached flexible beam. These numerical experiments are designed to stress-test various aspects of the two formulations in canonical and physiologic settings. More specifically, these cases represent three classes of flow simulations where (1) the boundary conditions and the solution are both steady, (2) the boundary condition and the solution are both unsteady, (3) the boundary conditions are steady but the solution is unsteady (in this case due to vortex shedding), and (4) the fluid domain is not fix with unsteady solution. All of the simulations are initialized using \(\mathbf {u_0} = 0\) and continued in time to reach cycle-to-cycle convergence or steady-state solutions. For apple-to-apple comparison, all parameters, except for the \(\tau _{\text {SUPG}}\) definition, are kept the same when comparing the present formulation against its conventional counterpart.

We first consider the case of flow in a straight pipe with steady boundary conditions, which can be considered the most simple and fundamental flow in the cardiovascular system with an existing analytical solution for comparison. In our simulations, a steady flow rate is imposed at the inlet and the pressure drop across the pipe is predicted using the present and conventional methods. The pipe has a length of 15 cm and a radius of 1 cm. A parabolic velocity profile is imposed at the inlet with an amplitude that results in a 10 mL/s flow rate. A zero Neumann boundary condition is imposed at the outlet. The dynamic viscosity is fixed at 1 g/cm-s. Three different densities, 1.571, 15.71, and 157.1 g/mL are selected to produce three Reynolds numbers (Re), 10, 100, and 1000, respectively. This way, we capture a wide range of Reynolds numbers that occur in cardiovascular flows60. For all cases, the pressure drop must remain the same according to the Hagen-Poiseuille analytical solution61, thereby allowing us to measure the accuracy of each method.

The mesh generated for this case contains 207,063 tetrahedral elements, which are used for the velocity and pressure interpolation as well as their test functions. Each Reynolds number is simulated using four different time step sizes of \(\Delta {t} = 10^{-1}\), \(10^{-2}\), \(10^{-3}\), and \(10^{-4}\) s. This range of time step sizes, which is typical in cardiovascular simulations, results in Courant–Friedrichs–Lewy numbers (\(\text{CFL} = {\overline{u}}\Delta t/\overline{\Delta x}\)) ranging from 5.2 to \(5.2\times 10^{-3}\), based on the mean element size of the mesh, \(\overline{\Delta x}\), and the mean flow velocity, \({\overline{u}}\). This range of CFL numbers, which captures under-resolved to over-resolved time discretizations, can be encountered in a typical simulation due to the differences in velocity and mesh resolution in different cardiovascular branches. Furthermore, multi-domain simulations can require the use of a smaller time step size (hence CFL) to ensure stability35.

All simulation cases were run in parallel with 32 cores. Equations are integrated for 5 s (approximately three flow-through times) to ensure steady conditions are reached. The \(l^2\)-norm of the residual is dropped by 3.5 orders of magnitude at each time step using Newton–Raphson iterations. Considering three Reynolds numbers, four time step sizes, and two formulations, we ran a total of 24 simulations for this case.

The results for this case are condensed in Fig. 1, which shows the predicted pressure drop, normalized by that of the Poiseuille solution61, as a function of the time step size for three different Reynolds numbers. The solutions calculated using the conventional \(\tau _{\text{SUPG}}\) formulation becomes exceedingly less accurate as the Reynolds number is increased and the time step size is reduced (Fig. 1a). That produces a very large error at \(\text{Re} =\)1000 and \(\Delta {t} = 10^{-4}\) where the predicted pressure drop is three times that of the analytical prediction. This large deviation from the reference solution confirms the inconsistency of the conventional formulation that we discussed earlier. This effect is amplified at higher Reynolds numbers when the artificial viscosity introduced through \(\tau _{\text{SUPG}}\) is larger in comparison to the physical viscosity, thus creating a larger variation in the solution as \(\tau _{\text{SUPG}}\) is varied with \(\Delta t\).

The present formulation, on the other hand, produces predictions that are independent of the time step size, confirming that it is a consistent formulation for steady-state flows (Fig. 1b). The overall error, which is primarily associated with spatial discretization, is negligible in comparison to the conventional formulation. The error for the case discussed above (\(\mathrm Re=1000\) and \(\Delta t=10^{-4}\)) drops from 300% to 0.7% when one uses the present rather than the conventional formulation.

Such a large difference in the solutions can be attributed to the large difference between \(\omega\) and \(2/\Delta t\) term in \(\tau _{\text {SUPG}}\) definition. The large difference between the two is depicted in Fig. 2 which shows the history of \(\omega\) over the course of the simulation. Given that \(\omega = 2/\Delta t\) at \(t=0\), the two methods are equivalent at the beginning of the simulation. However, as time progresses, the conventional formulation will significantly deviate from the present formulation by producing an \(\omega\) that differs by around fifteen orders of magnitude.

Given these attractive results and the fact that \(\omega \rightarrow 0\) in the present formulation, one may be tempted to entirely drop \(2/\Delta t\) term from the definition of \(\tau _{\text{SUPG}}\). As we argued earlier, such a modification results in a method that fails to converge particularly at a relatively small time step size and in cases where the flow is highly unsteady. Such regime corresponds to the initial stage of the pipe flow simulation, where the flow is rapidly evolving and \(\omega \ne 0\) (Fig. 2). Therefore, incorporating \(\omega\) into the definition of \(\tau _{\text{SUPG}}\) plays a crucial role in improving the stability of the overall scheme.

Lastly, we must note that although the present method converges for all cases considered, it produces a linear system that is stiffer than that of the conventional formulation. As a consequence, the solution of the linear system through an iterative solver will require more iterations, which can be over an order magnitude per time step when compared against the conventional method (Fig. 3). This larger number of iterations translates to a higher overall cost of these calculations, which is on average twice higher than that of the conventional formulation.

In this case, we compare the performance of the conventional formulation and our present formulation in a realistic cardiovascular geometry with a wide range of Reynolds and CFL numbers. The adopted anatomy represents that of an infant who has undergone the modified Blalock–Taussig shunt procedure55,56 (Fig. 4). The geometry contains multiple branches, among which an unsteady flow with a parabolic profile is imposed at the ascending aorta, which is interpolated from earlier multi-domain simulations57,62. Zero Neumann boundary conditions are imposed on all other branches, which are non-physiological, but nevertheless selected to highlight the difference between the two formulations. Since the flow rate through the pulmonary arteries is critical in understanding the performance of the shunt for this procedure, it is used here for simulation results comparison. The blood is assumed Newtonian with a density of 1.06 g/mL and a dynamic viscosity of 0.04 g/cm-s. The Reynolds number ranges from 50 to 1300 depending on the branch and the time within the cardiac cycle. Although the primary source of flow unsteadiness can be traced to the unsteady boundary condition, the complex geometry may also induce more local variations in the solution as a function of time, in this case (namely, generating frequency content in the solution that is not present in the boundary condition).

Figure 4

The meshed geometry and the inlet condition (the ascending aorta flow rate, \(Q_{\text {AoA}}\)) for the modified Blalock-Taussig shunt simulation.

Full size image

We performed simulations using time step sizes of \(\Delta {t} = 2.5 \times 10^{-2}\), \(2.5 \times 10^{-3}\), and \(2.5 \times 10^{-4}\) s, which are within the range of actual time step sizes used for these types of studies55,56,57. The geometry is discretized using 400,936 tetrahedral elements (Fig. 4). All simulations are run in parallel with 288 processors. Simulations are continued for at least six cardiac cycles (3 s) to ensure cycle-to-cycle convergence.

The predicted flow rate through the pulmonary arteries is extracted as a parameter of interest and plotted over one cardiac cycle in Fig. 5. For the conventional formulation, we can see a similar trend of solution deterioration as the time step size gets smaller. There is a 10% change in flow rate as the time step size is reduced from \(2.5 \times 10^{-2}\) to \(2.5 \times 10^{-3}\) s. There is another 15% deviation in the prediction of the conventional formulation when the time step size is further reduced by another order of magnitude, from \(2.5 \times 10^{-3}\) to \(2.5 \times 10^{-4}\) s. Such large variations in the results are rather alarming because 500 times steps per cardiac cycle or more (corresponding to \(\Delta t\le 10^{-3}\)) is a very modest number and has been commonly used in the past cardiovascular CFD studies18,55,56,57,63,64. According to our numerical experiment with the conventional method, such a commonly used time step size is already too small that it produces substantial error in the results.

In contrast to the conventional formulation, the method proposed is consistent with results that are almost independent of the time step size (the observed changes were less than 0.1%). \(\omega\) values using the present formulation is shown in Fig. 6 for three different time step sizes. At the beginning of the simulations, \(\omega\) quickly drops from that of the conventional formulation values to physics-based periodic values. As expected, the converged periodic \(\omega\) values are almost independent of the time step size, thus resulting in predictions that do not change as \(\Delta t \rightarrow 0\).

The total CPU time of the computations performed with the present formulation is around 1.5 times that of the conventional formulation, which is acceptable but warrants future improvements. Simulation results generated with the present formulation will enable researchers and clinicians to have high resolution in terms of time discretization without concerns about solution accuracy.

In this case, we will present a more textbook unsteady CFD numerical experiment to demonstrate that the consistency issue is not unique to cardiovascular applications. We will demonstrate that the present formulation solves the inconsistency issue for this case but also demonstrate where the limitation of our formulation currently lies. We will also use this case to discuss the domain size dependency of (or lack thereof) the present formulation for external flows.

We consider a two-dimensional unsteady flow over a square object with a steady inflow boundary condition. The geometry and mesh used for this case are shown in Fig. 7. The square has a side length of 1 m in a 12 m by 29.2 m fluid domain. The square obstacle is centered vertically and placed at a distance of 5 m from the inlet on the left side of the domain. Uniform horizontal flow with a velocity magnitude of 51.3 m/s is prescribed at the inlet and the outlet is a zero Neumann boundary. The top and bottom of the domain are both no-penetration boundaries (\(u_y\) = 0) with zero traction in the horizontal direction (\(h_x = 0\)). The fluid has a density of \(1.18 \times 10^{-3}\) kg/m3 and a dynamic viscosity of \(1.82 \times 10^{-4}\) kg/m-s. The Reynolds number of this case is 332, which results in vertex shedding downstream of the obstacle.

The mesh generated for this case contains 28, 502 triangular elements. Three different time step sizes are considered, \(\Delta {t} = 10^{-3}\), \(4 \times 10^{-4}\), and \(10^{-4}\) s. At each time step, the time integration residual is dropped by more than three orders of magnitude through Newton–Raphson iterations. The simulations are continued for 5 s to ensure statistically stationary conditions are established. We used 16 cores to perform these calculations.

The lift exerted on the square obstacle during the last 0.5 s of the simulation is shown in Fig. 8. Consistent with what we observed in the previous cases, this case also shows significant improvement in the results when the present design of \(\tau _{\text {SUPG}}\) is adopted. The conventional formulation prediction of the oscillation amplitude and period strongly depends on the time step size, with both values decreasing as \(\Delta t \rightarrow 0\) (Fig. 9a). In contrast, little change in these predictions is observed when the presented formulation is adopted (Fig. 9b). The contrast between the two methods can also be observed when comparing the pressure contours in Fig. 10. The snapshots shown in this figure are taken when the obstacle experiences a maximum lift. The dependency and lack of dependency of the conventional and present formulation on the time step size are evident in this figure.

The results from the two formulations are compared in terms of mean drag coefficient \(\overline{C_\text{d}}\), root mean square of drag coefficient fluctuation \(C_\text{d}'\), mean lift coefficient \(\overline{C_\text{l}}\), root mean square of lift coefficient fluctuation \(C_\text{l}'\), and Strouhal number St in Table 1. This comparison confirms our earlier observation that the present formulation nearly eliminates the dependence of the bulk flow parameters on the time step size.

Similar to the steady pipe flow, the convergence of the solutions as the time step size decreases is tied to the convergence of \(\omega\). With the present formulation, the value of \(\omega\), although not steady, converges to periodic values around 3 regardless of \(\Delta t\). The value of \(\omega ^2\), which is approximately 9 s\(^{-2}\), can be contrasted against \((2/\Delta t)^2\) that ranges from \(4\times 10^6\) to \(4\times 10^8\) s\(^{-2}\) for the simulated cases. That large change generates a significant variation in the solution as \(\Delta t\) is reduced in the conventional formulation.

To better see the effect of \(2/\Delta t\) term on \(\tau _{\text {SUPG}}\), we have produced a snapshot of \(\tau _{\text {SUPG}}\) over the entire computational domain for the two formulations in Fig. 11 at two different time step sizes. For the conventional formulation and in particular at \(\Delta t= 10^{-4}\), \(\tau _{\text {SUPG}}\) is very small and approximately equal to \(\Delta t/2\) over the entire domain. In the vicinity of the obstacle, where flow is the fastest, we observe some deviation from that baseline due to the contribution of the convective term in \(\tau _{\text {SUPG}}\). Nevertheless, that variation is not as significant as that of the present formulation where we see large changes in \(\tau _{\text {SUPG}}\) depending on the flow velocity, mesh resolution, and flow orientation relative to the element edges.

We also simulated the present flow over a square object case at a higher Reynolds number of 22, 000 to benchmark the two formulations against previously published results65,66,67,68. All simulation parameters are kept the same as in the case discussed above except for the dynamic viscosity which is reduced to \(2.75 \times 10^{-6}\) kg/m-s. Simulations are repeated using both formulations at \(\Delta t = 5 \times 10^{-3}\) s. The obtained results as well as those from the literature are summarized in Table 2. Note that the results from the literature were obtained from three-dimensional simulations whereas our computations are performed in two dimensions. Despite such differences, we observe a relatively good agreement with the published results. That is particularly the case for the present formulation that produces predictions closer to the reported range in comparison to the conventional formulation.

In order to investigate the effect of domain size for the present formulation, we extended the domain so that the overall domain size is four times that of the original domain. In doing so, we made sure that the mesh for the subset of the domain corresponding to the original mesh in Fig. 7 remains unchanged so that the reported results are minimally affected by this change in the domain size. We repeated the original \(\text{Re} = 332\) simulations for both formulations using \(\Delta t = 10^{-3}\) with the extended domain. The results are summarized in Table 3. The variation in results is almost identical for the two formulations indicating that the change is more likely a result of moving the locations of the boundaries rather than the change in the value of \(\omega\). Following what we argued earlier, extending the domain reduced \(\omega\) to approximately half its original value. This change, nevertheless, has a negligible effect on the overall results given that the \(\omega\) term is much smaller than the sum of the other terms appearing in \(\tau _{\text{SUPG}}\).

The relatively small value of \(\omega\), in comparison to the sum of the two other terms appearing in \(\tau _{\text{SUPG}}\) in Eq. (8), is a general feature of the present formulation rather than being unique to the case above. Evidently, it was also relatively small for the vascular model discussed earlier as the results obtained from the proposed formulation closely resembled that of the conventional formulation at large \(\Delta t\) (Fig. 5), where the sum of the other two terms in \(\tau _{\text{SUPG}}\) is dominant.

In general, if we only consider \(\textbf{u}^h \cdot \varvec{\xi } \textbf{u}^h\) relative to \(\omega ^2\), we can show the ratio of the two scales as the square of \(u/(\omega \Delta x)\). For the present method to be domain-size-dependent, that ratio must be smaller than one, implying that the flow at a node must do a full oscillation before it has the time to advect the fluid across a single element. That is an extremely fast oscillating flow, which even if occurs in reality, requires a much smaller time step to resolve, thereby indicating that the present method will do much better than the conventional method.

We should note that the linear solver convergence, which was not an issue for the previous two cases, posed a problem in this case. In general, the linear solver takes longer to converge as the time step sizes are reduced. In extreme cases, the linear solver may not converge at all. One such scenario occurs when the time step size is very small and there is a significant change in the solution between time steps, e.g., the present simulation starting from a zero velocity field. To overcome such issues, one may start a simulation with a larger time step size and then reduce it to the target value after a few time steps. Nevertheless, this convergence issue is a shortcoming of the present formulation that warrants future research.

In this section, we will demonstrate the generalization of the present formulation to moving domain configurations using a fluid–structure interaction (FSI) test case. The chosen case is a 2D flow over a fixed square with an attached flexible beam, which has been used in the past for verification of the FSI codes46,49.

The specifications of the test case are shown schematically in Fig. 12. The computation domain is \(12\times 19.5\) m. The square is \(1\times 1\) m with the center of the square placed 5 m from the inlet. The beam is \(4\times 0.06\) m centered behind the square. The fluid domain setup is identical to the previous section’s flow over a square case, where the fluid density is \(1.18 \times 10^{-3}\) kg/m3 and the dynamic viscosity is \(1.82 \times 10^{-4}\) kg/m-s with uniform inlet horizontal velocity at 51.3 m/s, resulting in Re\(=332\). The \(\omega\) in the stabilization constant, \(\tau _{\text {SUPG}}\), is calculated in an arbitrary Eulerian–Lagrangian framework, as specified in Eq. 10.

The solid domain is modeled as a St. Venant-Kirchhoff elastic solid69. The discretized weak form of the solid problem is given as

where

is the second Piola-Kirchhoff stress, in which

is the Green-Lagrange strain tensor, with

as the Cauchy-Green deformation tensor, where

is the deformation matrix and \(\textbf{d}\) is the displacement vector. The density of the beam at \(t=0\) is \(\rho _s^0 = 0.1 \mathrm {kg/m^3}\). The Young’s modulus, \(E = 2.5 \times 10^6\) kg/(s2m), and Poisson’s ratio, \(\nu = 0.35\) are used to calculate Lamé parameters \(\lambda\) and \(\mu\) in Eq. 12 as

The details of the FSI formulation are presented in a previous paper and not repeated here for brevity49. In summary, an arbitrary Lagrangian–Eulerian (ALE) approach and a quasi-direct FSI solution strategy are employed9,46,70,71. The increments of the fluid and structure solution are monolithically computed. Jacobian-based stiffening is used for the elastic mesh motion without re-meshing9,72,73. A back-flow stabilization scheme is employed at the Neumann boundaries to prevent simulation divergence caused by partial flow reversal at the outlet74.

The mesh used for all cases contains roughly 20 thousand triangular elements for the combined fluid and solid domains (Fig. 12). The simulations were run for 10 s with two different time step sizes \(\Delta t = 1\times 10^{-3}\) and \(5\times 10^{-4}\) s using the conventional and present formulations of \(\tau _{\text {SUPG}}\), resulting in four simulations in total.

Figure 13

The beam tip vertical displacement, \(d_\text{y}\), for the case shown in Figure 12 computed using the conventional (a; left) and present formulation (b; right). The two curves shown as a function of time for each case correspond to the simulations run with \(\Delta t = 1\times 10^{-3}\) (dashed line) and \(5\times 10^{-4}\) (solid line).

Full size image

Figure 13 shows the vertical displacement of the beam tip (marked in Fig. 12). The conventional formulation produces an oscillation period that slightly varies as the time step size is changes by a factor of two. The change in period is approximately 0.7%, increasing from 0.305 s at \(\Delta t = 5\times 10^{-4}\) s to 0.307 s at \(\Delta t = 10^{-3}\) s. On the contrary, the present formulation produced significantly closer results for both time step sizes. The predicted period, in this case, is approximately 0.299 s, which agrees well with the literature46,49,70.

In terms of stability and computational cost, we observe a behavior similar to the flow over the square case from the previous section. For cases reported above, the total simulation cost of the present formulation is around 4 times that of the conventional formulation (1.3 vs. 5.1 h using 16 cores). Using a time step size smaller than \(\Delta t = 5\times 10^{-4}\) s causes linear solver convergence issues in the case of the present formulation.

Taylor, C. A., Hughes, T. J. & Zarins, C. K. Finite element modeling of blood flow in arteries. Comput. Methods Appl. Mech. Eng. 158, 155–196. https://doi.org/10.1016/S0045-7825(98)80008-X (1998).

Article  ADS  MathSciNet  MATH  Google Scholar 

Taylor, C. A., Hughes, T. J. & Zarins, C. K. Finite element modeling of three-dimensional pulsatile flow in the abdominal aorta: Relevance to atherosclerosis. Ann. Biomed. Eng. 26, 975–987. https://doi.org/10.1114/1.140 (1998).

Article  CAS  PubMed  Google Scholar 

Taylor, C. A. et al. Predictive medicine: Computational techniques in therapeutic decision-making. Comput. Aid. Surg. 4, 231–247 (1999).

Article  CAS  Google Scholar 

Eymard, R., Gallouët, T. & Herbin, R. Finite volume methods. Handb. Numer. Anal. 7, 713–1018. https://doi.org/10.1016/S1570-8659(00)07005-8 (2000).

Article  MathSciNet  MATH  Google Scholar 

LeVeque, R. J. Finite difference methods for ordinary and partial differential equations. Soc. Ind. Appl. Math.https://doi.org/10.1137/1.9780898717839 (2007).

Article  MATH  Google Scholar 

Hughes, T. J. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis (Prentice-Hall Inc, 1987).

MATH  Google Scholar 

Van de Vosse, F. et al. Finite-element-based computational methods for cardiovascular fluid–structure interaction. J. Eng. Math. 47, 335–368. https://doi.org/10.1023/B:ENGI.0000007985.17625.43 (2003).

Article  MathSciNet  MATH  Google Scholar 

Bazilevs, Y., Calo, V. M., Zhang, Y. & Hughes, T. Isogeometric fluid–structure interaction analysis with applications to arterial blood flow. Comput. Mech. 38, 310–322. https://doi.org/10.1007/s00466-006-0084-3 (2006).

Article  MathSciNet  MATH  Google Scholar 

Bazilevs, Y., Hsu, M.-C., Benson, D. J., Sankaran, S. & Marsden, A. L. Computational fluid–structure interaction: methods and application to a total cavopulmonary connection. Comput. Mech. 45, 77–89. https://doi.org/10.1007/s00466-009-0419-y (2009).

Article  MathSciNet  MATH  Google Scholar 

Quarteroni, A., Veneziani, A. & Vergara, C. Geometric multiscale modeling of the cardiovascular system, between theory and practice. Comput. Meth. Appl. Mech. Eng. 302, 193–252. https://doi.org/10.1016/j.cma.2016.01.007 (2016).

Article  ADS  MathSciNet  MATH  Google Scholar 

Markl, M. et al. Time-resolved three-dimensional phase-contrast MRI. J. Magn. Reson. Imag. 17, 499–506. https://doi.org/10.1002/jmri.10272 (2003).

Article  Google Scholar 

Dyverfeldt, P. et al. 4d flow cardiovascular magnetic resonance consensus statement. J. Cardiovasc. Magn. Reson. 17, 1–19. https://doi.org/10.1186/s12968-015-0174-5 (2015).

Article  Google Scholar 

Schulz-Menger, J. et al. Standardized image interpretation and post-processing in cardiovascular magnetic resonance-2020 update. J. Cardiovasc. Magn. Reson. 22, 1–22. https://doi.org/10.1186/s12968-020-00610-6 (2020).

Article  Google Scholar 

Markl, M., Frydrychowicz, A., Kozerke, S., Hope, M. & Wieben, O. 4d flow MRI. J. Magn. Reson. Imag. 36, 1015–1036. https://doi.org/10.1002/jmri.23632 (2012).

Article  Google Scholar 

Antiga, L. et al. An image-based modeling framework for patient-specific computational hemodynamics. Med. Biol. Eng. Comput. 46, 1097–1112. https://doi.org/10.1007/s11517-008-0420-1 (2008).

Article  PubMed  Google Scholar 

Taylor, C. & Figueroa, C. Patient-specific modeling of cardiovascular mechanics. Ann. Rev. Biomed. Eng. 11, 109–134. https://doi.org/10.1146/annurev.bioeng.10.061807.160521 (2009).

Article  CAS