A time-consistent stabilized finite element method for fluids with applications to hemodynamics

Several finite element methods for simulating incompressible flows rely on the streamline upwind Petrov–Galerkin stabilization (SUPG) term, which is weighted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\text{SUPG}}$$\end{document}τSUPG. The conventional formulation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\text{SUPG}}$$\end{document}τSUPG includes a constant that depends on the time step size, producing an overall method that becomes exceedingly less accurate as the time step size approaches zero. In practice, such method inconsistency introduces significant error in the solution, especially in cardiovascular simulations, where small time step sizes may be required to resolve multiple scales of the blood flow. To overcome this issue, we propose a consistent method that is based on a new definition of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\text{SUPG}}$$\end{document}τSUPG. This method, which can be easily implemented on top of an existing streamline upwind Petrov–Galerkin and pressure stabilizing Petrov–Galerkin method, involves the replacement of the time step size in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\text{SUPG}}$$\end{document}τSUPG with a physical time scale. This time scale is calculated in a simple operation once every time step for the entire computational domain from the ratio of the L2-norm of the acceleration and the velocity. The proposed method is compared against the conventional method using four cases: a steady pipe flow, a blood flow through vascular anatomy, an external flow over a square obstacle, and a fluid–structure interaction case involving an oscillatory flexible beam. These numerical experiments, which are performed using linear interpolation functions, show that the proposed formulation eliminates the inconsistency issue associated with the conventional formulation in all cases. While the proposed method is slightly more costly than the conventional method, it significantly reduces the error, particularly at small time step sizes. For the pipe flow where an exact solution is available, we show the conventional method can over-predict the pressure drop by a factor of three. This large error is almost completely eliminated by the proposed formulation, dropping to approximately 1% for all time step sizes and Reynolds numbers considered.

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 τ, to prevent nonphysical oscillations inherent to the Galerkin method in strongly convective regimes [21][22][23][24] .One of the most commonly used formulations of τ has been the one proposed in the streamline upwind Petrov-Galerkin formulations (SUPG) 21 , τ SUPG , which is also adopted in others such as the residual-based variational multiscale (RBVMS) formulation 25,26 .
Traditionally, the steady form of τ 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 numbers 21,22,27,28 .While this steady form of τ 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 (∆t).The traditional strategy to overcome this issue has been adding a ∆t-dependent term to the definition of τ SUPG 22, 28, 29 .The added ∆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 zero [30][31][32] (also see the results section).As a consequence, this time step size-dependent design of τ 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 modeling [33][34][35][36] .In fact, some of the most popular packages for cardiovascular simulation 37,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 τ SUPG 30, 39, 40 .T.E.Tezduyar and others introduced an element-vector-based τ SUPG that uses the relative elemental magnitude of terms in the weak form of the Navier-Stokes equations 32,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 time 31 .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 simulations 37,42 .
In this study, we propose a formulation for τ 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 τ 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 τ SUPG is replaced by the spectral mode number 43 .The same parameter, when expressed in a spatio-temporal setting, inspires the use of a flowdependent time scale in the τ SUPG definition, hence motivating the present design.
The idea of replacing the time step size in τ SUPG with an acceleration-to-velocity ratio has been proposed independently earlier by J.Evans and others 44 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 τ 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 τ 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 τ 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 twodimensional 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.

Formulation
The Navier-Stokes equations for incompressible flows are stated as where ρ is the density, u(x,t) is the velocity, f(x,t) is the external forcing, Ω × [0, T ] is the fluid computational spatio-temporal domain, and the stress tensor where p(x,t) is pressure and µ is the dynamic viscosity.The Dirichlet and Neumann boundary conditions are defined as respectively, where Γ g and Γ h are subsets of the boundary Γ where the Dirichlet and Neumann boundaries are prescribed, n is the outward unit normal vector, and g and 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 u h ∈ S h u and p h ∈ S h p such that for all w h ∈ V h u and q h ∈ V h p , In the above problem statement, S h u and S h p are the discrete solution spaces for the velocity and pressure, respectively, and V h u and V h p are the finite-dimensional test function spaces for the velocity and pressure, respectively.In Equation ( 4), the terms directly obtained from Equation (1) are supplemented with three elemental stabilization terms.The two terms multiplied by τ SUPG are the conventional SUPG and PSPG stabilizations to ensure stability in strongly convec- tive flows and allow for equal order interpolation functions for velocity and pressure, respectively 21,22 .The term involving ν 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 comparisons 28,45,48 .In Equation ( 5), ν = µ/ρ is the kinematic viscosity, ∆t is the time step size, ξ ξ ξ is the covariant tensor obtained from the mapping of the physicalparent elements, and C I is a shape-function-dependent constant, which is 3 in our study.
The inconsistency of the above-mentioned stabilized formulation is caused by ∆t in τ SUPG .In a steady flow, in which the solution should be independent of ∆t, τ SUPG and thus the overall added diffusion will change with ∆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 domain 43 .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 ω is the spectral mode appearing as a source term in the frequency formulation of the unsteady Stokes equations.This spectral formulation of τ closely resembles the conventional definition of τ SUPG in Equation (5) if 2/∆t is replaced by ω.
The u h • ξ ξ ξ u h term does not appear in Equation ( 7) as the convective acceleration term is not present in the unsteady Stokes equations.Adding this term into Equation (7) and incorporating the O(1) constant C I results in which is adopted in this study to replace the conventional definition of τ SUPG from Equation (5) when solving Equation (4).
The new definition of τ SUPG in Equation (8) becomes identical to the traditional formulation (Equation ( 5)) if ω = 2/∆t.This value is close to the largest frequency associated with the time discretization, namely π/∆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 ω than π/∆t.The distinction of these two frequencies inspires the proposed definition of τ SUPG .
It is straightforward to evaluate Equation (8) in a spectral formulation as ω is the computed frequency and readily available as an independent parameter.However, its adoption is not straightforward in a traditional time formulation, where ω does not appear as an independent parameter.Ideally, ω 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 τ 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 This formulation of ω is designed to go to zero as the flow reaches a steady state, where ∂ u h ∂t goes to zero.We will show in the results section that this formulation is consistent in both steady and unsteady flows as ∆t → 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 Equation ( 4) in an arbitrary Eulerian-Lagrangian framework, u h in the convective acceleration term is replaced by the fluid velocity relative to moving mesh u h − ûh , where ûh denotes mesh velocity.It is this velocity that is employed in the definition of τ SUPG and also ω in Equation ( 9), changing it to where the acceleration term is measured at the mesh node location x and integrals performed over the fluid domain Ω f .By subtracting the mesh velocity from the fluid velocity in Equation ( 10), the resulting definition of ω 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 ω 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 past 44 , fails to converge in our implicit formulation.This lack of convergence can be attributed to the velocity appearing in the denominator of Equation ( 9), thereby creating widely varying ω 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 Equation (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 ∂ u h ∂t L 2 term receives non-zero contribution only from regions near the obstacle whereas u h L 2 receives a contribution from the entire domain.In this setting, ω 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 ω 2 term in Equation ( 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 ∆t.Therefore, increasing the domain size will decrease a parameter in the definition of τ SUPG that is already small, thus hardly changing the solution.
Even though the ω 2 value in τ SUPG is very small, dropping it from its definition will create convergence issues, as it has been established in the past 27,28 .The reason that the inclusion of the ω 2 term in τ SUPG prevents such scenarios is that a widely varying solution in time will lead to a relatively large ω.As a result, the contribution of ω 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 Equation (9).To avoid such possibilities, we set ω = 2/∆t at the first time step, recovering the conventional definition of τ SUPG .
As detailed in a previous publication 49 , we use the implicit generalized-α time integration scheme 50 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, ∂ u h ∂t and u h in Equation ( 9) are explicitly computed from those variables in a single operation.As we will demonstrate in the results, ω 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-α method also provides access to variables at the intermediate time points (i.e., n + α m and n + α f for the acceleration and velocity, respectively), one may elect to use those intermediate variables to update ω 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 τ SUPG at the Gauss quadrature points based on the intermediate variables are left unchanged and are identical to the conventional formulation.

Simulations and Results
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 equations 49,[51][52][53] .The solver has been verified 54 and extensively employed for cardiovascular modeling in the past [55][56][57] .This solver is parallelized using a message passing interface (MPI).The workload is parallelized using spatial partitioning by employing ParMETIS library 58 .All computations are performed on a cluster of AMD Opteron 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) method 59 with a tolerance of 10 −2 .
Four cases are simulated using both the conventional formulation (Equation ( 5)) and the present formulation (Equation ( 8)) for τ SUPG : a pipe flow, flow in a modified Blalock-Taussig shunt geometry 55 , 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 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 τ SUPG definition, are kept the same when comparing the present formulation against its conventional counterpart.

Steady pipe flow
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 flows 60 .For all cases, the pressure drop must remain the same according to the Hagen-Poiseuille analytical solution 61 , 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 ∆t = 10 −1 , 10 −2 , 10 −3 , and 10 −4 seconds.This range of time step sizes, which is typical in cardiovascular simulations, results in Courant-Friedrichs-Lewy numbers (CFL = u∆t/∆x) ranging from 5.2 to 5.2 × 10 −3 , based on the mean element size of the mesh, ∆x, and the mean flow velocity, u.This range of CFL numbers, which captures under-resolved to overresolved 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 stability 35 .
All simulation cases were run in parallel with 32 cores.Equations are integrated for five seconds (approximately three flowthrough 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 Figure 1, which shows the predicted pressure drop, normalized by that of the Poiseuille solution 61 , as a function of the time step size for three different Reynolds numbers.The solutions calculated using the conventional τ SUPG formulation becomes exceedingly less accurate as the Reynolds number is increased and the time step size is reduced (Figure 1(a)).That produces a very large error at Re =1,000 and ∆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 τ SUPG is larger in comparison to the physical viscosity, thus creating a larger variation in the solution as τ SUPG is varied with ∆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 (Figure 1(b)).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 (Re = 1, 000 and ∆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 ω and 2/∆t term in τ SUPG definition.The large difference between the two is depicted in Figure 2 which shows the history of ω over the course of the simulation.Given that ω = 2/∆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 ω that differs by around fifteen orders of magnitude.Given these attractive results and the fact that ω → 0 in the present formulation, one may be tempted to entirely drop 2/∆t term from the definition of τ 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 ω = 0 (Figure 2).Therefore, incorporating ω into the definition of τ 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 (Figure 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.

Blood flow in vascular anatomy
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 procedure 55,56 (Figure 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 multidomain simulations 57,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).We performed simulations using time step sizes of ∆t = 2.5 × 10 −2 , 2.5 × 10 −3 , and 2.5 × 10 −4 seconds, which are within the range of actual time step sizes used for these types of studies [55][56][57] .The geometry is discretized using 400, 936 tetrahedral elements (Figure 4).All simulations are run in parallel with 288 processors.Simulations are continued for at least six cardiac cycles (3 seconds) 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 Figure 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 × 10 −2 to 2.5 × 10 −3 seconds.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 × 10 −3 to 2.5 × 10 −4 seconds.Such large variations in the results are rather alarming because 500 times steps per cardiac cycle or more (corresponding to ∆t ≤ 10 −3 ) is a very modest number and has been commonly used in the past cardiovascular CFD studies 18, 55-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%).ω values using the present formulation is shown in Figure 6 for three different time step sizes.At the beginning of the simulations, ω quickly drops from that of the conventional formulation values to physics-based periodic values.As expected, the converged periodic ω values are almost independent of the time step size, thus resulting in predictions that do not change as ∆t → 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.

9/21
Flow over a square 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 Figure 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 × 10 −3 kg/m 3 and a dynamic viscosity of 1.82 × 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, ∆t = 10 −3 , 4 × 10 −4 , and 10 −4 seconds.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 seconds 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 seconds of the simulation is shown in Figure 8.Consistent with what we observed in the previous cases, this case also shows significant improvement in the results when the present design of τ 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 ∆t → 0 (Figure 9a).In contrast, little change in these predictions is observed when the presented formulation is adopted (Figure 9b).The contrast between the two methods can also be observed when comparing the pressure contours in Figure 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.Similar to the steady pipe flow, the convergence of the solutions as the time step size decreases is tied to the convergence of ω.With the present formulation, the value of ω, although not steady, converges to periodic values around 3 regardless of ∆t.The value of ω 2 , which is approximately 9 s −2 , can be contrasted against (2/∆t) 2 that ranges from 4 × 10 6 to 4 × 10 8 s −2 for the simulated cases.That large change generates a significant variation in the solution as ∆t is reduced in the conventional formulation.
To better see the effect of 2/∆t term on τ SUPG , we have produced a snapshot of τ SUPG over the entire computational domain for the two formulations in Figure 11 at two different time step sizes.For the conventional formulation and in particular at ∆t = 10 −4 , τ SUPG is very small and approximately equal to ∆t/2 over the entire domain.In the vicinity of the obstacle,

12/21
where flow is the fastest, we observe some deviation from that baseline due to the contribution of the convective term in τ SUPG .Nevertheless, that variation is not as significant as that of the present formulation where we see large changes in τ 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 results [65][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 × 10 −6 kg/m-s.Simulations are repeated using both formulations at ∆t = 5 × 10 −3 seconds.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.Table 2. Results for the flow over a square object at Re = 22, 000 using the conventional and present formulations, and their comparison against the literature.
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 13/21 corresponding to the original mesh in Figure 7 remains unchanged so that the reported results are minimally affected by this change in the domain size.We repeated the original Re = 332 simulations for both formulations using ∆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 ω.Following what we argued earlier, extending the domain reduced ω to approximately half its original value.This change, nevertheless, has a negligible effect on the overall results given that the ω term is much smaller than the sum of the other terms appearing in τ SUPG .The relatively small value of ω, in comparison to the sum of the two other terms appearing in τ SUPG in Equation ( 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 ∆t (Figure 5), where the sum of the other two terms in τ SUPG is dominant.
In general, if we only consider u h • ξ ξ ξ u h relative to ω 2 , we can show the ratio of the two scales as the square of u/(ω∆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.

Fluid-structure interaction
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 codes 46,49 .10.

14/21
The solid domain is modeled as a St. Venant-Kirchhoff elastic solid 69 .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 d is the displacement vector.The density of the beam at t = 0 is ρ 0 s = 0.1kg/m 3 .The Young's modulus, E = 2.5 × 10 6 kg/(s 2 m), and Poisson's ratio, ν = 0.35 are used to calculate Lamé parameters λ and µ in Equation 12as The details of the FSI formulation are presented in a previous paper and not repeated here for brevity 49 .In summary, an arbitrary Lagrangian-Eulerian (ALE) approach and a quasi-direct FSI solution strategy are employed 9,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-meshing 9,72,73 .A back-flow stabilization scheme is employed at the Neumann boundaries to prevent simulation divergence caused by partial flow reversal at the outlet 74 .
The mesh used for all cases contains roughly 20 thousand triangular elements for the combined fluid and solid domains (Figure 12).The simulations were run for 10 seconds with two different time step sizes ∆t = 1 × 10 −3 and 5 × 10 −4 seconds using the conventional and present formulations of τ SUPG , resulting in four simulations in total.Figure 13 shows the vertical displacement of the beam tip (marked in Figure 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 ∆t = 5 × 10 −4 s to 0.307 s at ∆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 literature 46,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 versus 5.1 hours using 16 cores).Using a time step size smaller than ∆t = 5 × 10 −4 s causes linear solver convergence issues in the case of the present formulation.

Conclusion
In this study, we propose a new formulation for the stabilization parameter (τ SUPG ) that appears in the streamline upwind Petrov-Galerkin and pressure stabilizing Petrov-Galerkin method (SUPG/PSPG) to overcome the historical limitation of the conventional formulation that produces a large error at the small time step size.The proposed formulation uses a flow time scale instead of the time step size to account for the contribution of the acceleration term to τ SUPG .Using the new formula- tion of τ SUPG , we successfully produce an overall stable technique that is consistent with regard to the time step size.The new definition of τ SUPG (Equation (8)) is simple to implement in existing SUPG/PSPG formulated fluid solvers.Although the present formulation comes at the cost of increasing the number of linear solver iterations, it significantly improves the overall accuracy of the stabilized finite element methods for computational fluid dynamics, making it an attractive choice for cardiovascular simulations.

Figure 1 .
Figure 1.The predicted pressure drop normalized by the analytical solution (∆P/∆P ref ) for the steady pipe flow case as a function of the time step size (∆t) for (a) the conventional formulation and (b) the present formulation.The three Reynolds numbers are 10 (red circle), 100 (blue square), and 1,000 (black triangle).

Figure 3 .
Figure3.The average number of the linear solver (GMRES) iterations per time step ( Nlsitr ) as a function of the time step size (∆t) for the steady pipe flow case.The results correspond to the three simulated Reynolds numbers: 10 (red circle), 100 (blue square), and 1,000 (black triangle) using the conventional (solid line) and the present formulation (dot-dashed line) of τ SUPG .

Figure 4 .
Figure 4.The meshed geometry and the inlet condition (the ascending aorta flow rate, Q AoA ) for the modified Blalock-Taussig shunt simulation.

Figure 7 .
Figure 7.The mesh constructed for the flow over a square obstacle simulation.

Figure 10 . 15 Table 1 .
Figure 10.Pressure contour for the flow over a square obstacle (Figure 7) captured at the maximum lift.Re = 332.(a) and (c) are the results obtained from the conventional formulation while (b) and (d) are the results obtained from the present formulation.(a) and (b) are obtained using ∆t = 10 −3 and (c) and (d) using ∆t = 10 −4 .

Figure 11 .
Figure 11.The snapshot of τ SUPG obtained at peak lift for the conventional formulation (a,c) and the present formulation (b,d) with ∆t = 10 −3 (a,b) and 10 −4 (c,d) seconds.Re = 332.Note the range used for the color bars.

Figure 12 .
Figure 12.The schematic of the FSI case involving a flexible beam attached to a solid square in a cross-flow.

4 Figure 13 .
Figure 13.The beam tip vertical displacement, d 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 ∆t = 1 × 10 −3 (dashed line) and 5 × 10 −4 (solid line).

Table 3 .
Result comparison between the conventional formulation and present formulation using the original and the extended domains for the flow over a square case.∆t = 10 −3 .Re = 332.