A time-spectral Stokes solver for simulation of time-periodic flows in complex geometries
aa r X i v : . [ c s . C E ] S e p A time-spectral Stokes solver for simulation of time-periodic flows in complexgeometries
Chenwei Meng a , Mahdi Esmaily a a Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY 14850, USA
Abstract
Simulation of unsteady creeping flows in complex geometries has traditionally required the use of a time-steppingprocedure, which is typically costly and unscalable. To reduce the cost and allow for computations at much largerscales, we propose an alternative approach that is formulated based on the unsteady Stokes equation expressed inthe time-spectral domain. This transformation results in a boundary value problem with an imaginary source termproportional to the computed mode that is discretized and solved in a complex-valued finite element solver usingBubnov-Galerkin formulation. This transformed spatio-spectral formulation presents several advantages over the tra-ditional spatio-temporal techniques. Firstly, for cases with boundary conditions varying smoothly in time, it providesa significant saving in computational cost as it can resolve time-variation of the solution using a few modes rather thanthousands of time steps. Secondly, in contrast to the traditional time integration scheme with a finite order of accuracy,this method exhibits a super convergence behavior versus the number of computed modes. Thirdly, in contrast to thestabilized finite element methods for fluid, no stabilization term is employed in our formulation, producing a solutionthat is consistent and more accurate. Fourthly, the proposed approach is embarrassingly parallelizable owing to theindependence of the solution modes, thus enabling scalable calculations at a much larger number of processors. Thecomparison of the proposed technique against a standard stabilized finite element solver is performed using two- andthree-dimensional canonical and complex geometries. The results show that the proposed method can produce moreaccurate results at 1% to 11% of the cost of the standard technique for the studied cases.
1. Introduction
Fast and cost-e ff ective simulation of flow in complex geometries is highly desirable, at least for two sets ofapplications. The first are the time-critical problems where the results must be obtained in a given time window, thusimposing an upper bound on the time-to-solution. The second are those with a bound on the computational cost, wherethe number of CPU-hours spent on the simulation must be limited. While these two restrictions on time and cost areidentical for sequential simulations, they may well diverge for parallel implementations depending on the scalabilityof the application.A prominent area where all the above conditions and restrictions apply is cardiovascular modeling. The flowis unsteady and time-periodic, the geometry is often complex, the cost must remain low, and the available time forsimulation is often limited. A specific example within this realm is patients with congenital heart disease who mustundergo an operation right after birth [1]. Utilizing cardiovascular modeling as a part of patient-specific surgicalplanning will require the simulation to be completed between diagnosis (birth) and the time of operation (a few daysafter birth at most) [2]. Moreover, performing shape optimization to improve surgical anatomy relies on havinga highly e ffi cient computational framework since a formal optimization study typically requires tens to hundredsof simulations [3]. A similar constraint exists on the e ffi ciency of the underlying computational method if suchcomputations were to be performed outside of the academic setting and at an industrial scale, where access to high-performance computing resources for a large number of patients remains limited due to the cost considerations [4].While there exists a broad range of numerical methods for simulating flow in complex geometries, here we con-sider stabilized finite element techniques to establish the relative accuracy and cost-e ffi ciency of our method sincethese methods are widely adopted for cardiovascular blood flow simulation [5] and closest to the method proposedin this study. Within the realm of stabilized finite element methods (FEM), we particularly focus on residual-based Preprint submitted to Elsevier September 16, 2020 ariational multiscale method (RBVMS), which is constructed on top of streamwise upwind Petrov-Galerkin (SUPG)technique [6, 7, 8].As is the case with the majority of the time-dependent partial di ff erential equations solvers, the application ofRBVMS technique to a time-periodic problem (e.g., blood flow in the circulatory system) requires discretization ofthe underlying governing equations in time. Thus, the cost of these methods is inherently dependent on the timestep size, increasing as time step size decreases. While the maximum allowable time step size is limited by thestability considerations for explicit time integration schemes, it is often more relaxed for implicit schemes, where itis determined solely based on the desired accuracy [9, 10]. The larger time step size for implicit schemes does notnecessarily imply their lower cost, as the fewer time steps can be o ff set by the larger number of Newton-Raphsoniteration per time steps. Furthermore, the high cost associated with the time integration is exacerbated by the transientnature of the solution that requires simulation of many cycles to achieve cycle-to-cycle convergence. In addition tothese setbacks, one well-known issue associated with the RBVMS formulation is its lack of consistency with regard totime step size that, for instance, lead to the time step size dependency of the final solution obtained from a steady-stateproblem.A remedy for the aforementioned issue is to use a time-spectral technique [11, 12, 13]. The underlying idea ofa time-spectral method is to transfer the temporal problem to a frequency domain using Fourier transformation andsubsequently solve the transformed boundary value problem mode-by-mode. Since the spectral methods bypass thetime integration, their accuracy and stability are not a ff ected by the time integration scheme. We should note that thenumber of modes in the case of a time-spectral solver, at least in theory, must be equal to the number of time stepsfor all time scales in the problem to be captured. However, for cardiovascular flows, the majority of those modes canbe neglected with minimal error, as we will show later in this article. In practice, a few modes (less than ten) su ffi cesfor su ffi ciently smooth and time-periodic boundary conditions, thus o ff ering a significant reduction in overall cost incomparison to a traditional solver that requires thousands of time steps for simulation of a single cardiac cycle [14].On top of this saving, the spectral technique can save on the computation required for cycle-to-cycle convergence asit is based on a boundary value problem that does not include a transient solution. Additionally, for the linear Stokesequation that is the focus of this study, all the solution modes are orthogonal and independent. Thus, the spectraltechnique can be implemented as an embarrassingly parallelizable scheme with each group of processors computingthe solution associated with only one of the many modes. Therefore, by formulating the original temporal problemas a set of orthogonal boundary value problems, one can significantly reduce the total computational cost and, at thesame time, scale computations to a much larger number of processors.In this paper, we propose a time-spectral technique named Complex-valued Stokes Solver (SCVS) for cardiovas-cular blood flow simulations. In the past, the time-spectral technique has been employed for simulation of flow overblu ff bodies, such as blades in a compressor rotor and a pitching airfoil, using a finite volume approach [15, 16, 17, 18].In the present study, we introduce this technique for the simulation of blood flow in the cardiovascular system. Whatdistinguishes cardiovascular systems from prior applications is the geometrical complexity, as such models typicallyinvolve multiple inlets and outlets that may be coupled to lower-order models as boundary conditions [19]. Secondly,this technique is implemented using finite element methods, which is a prominent technique for cardiovascular mod-eling [20, 21]. Additionally, in this study, we discuss in detail the convergence, accuracy, and computational e ffi ciencyof the SCVS as a function of element size, number of modeled modes, and the underlying linear solver tolerance.The paper is organized as follows: we introduce the governing equations, namely Stokes equation with a complex-valued source term, in the spectral domain. Second, we present the discrete formulation of those equations based onfinite element formulation. Third, we evaluate the introduced technique using four cases, where we compare thesolutions obtained from the SCVS and a standard RBVMS solver. The convergence and accuracy of the SCVS areestablished for a canonical case where an analytical solution is available. The computational e ffi ciency of the SCVSis demonstrated and contrasted against that of the RBVMS before drawing conclusions at the end.
2. Complex-valued Stokes solver formulation
In what follows, we first briefly describe the complex-valued Stokes equation as the governing equation for lowReynolds number flows, discuss its discrete formulation using FEM and then present the SCVS algorithm in detail.2 .1. Governing equations
Consider fluid flow in a domain Ω with boundary Γ = ∂ Ω that is governed by the incompressible Navier-Stokesequations as ρ ∂ u ∂ t + ρ u · ∇ u = −∇ p + ∇ · ( µ ∇ u ) in Ω , ∇ · u = Ω , u = g on Γ g , ( − p I + µ ∇ u ) · n = h on Γ h , (1)where u ( x , t ) refers to the velocity of the fluid at location x and time t , p ( x , t ) is the pressure, ρ is the density, n is theboundary Γ outward normal vector, and µ is the viscosity. g and h refer to the imposed velocity and traction on theDirichlet Γ g and Neumann Γ h boundaries where Γ = Γ g ∪ Γ h . Additionally, suppose that these boundary conditions aretime-periodic with a period of T . At su ffi ciently low Reynolds number, the nonlinear term ρ u · ∇ u can be neglected inEq. (1), producing ρ ∂ u ∂ t = −∇ p + ∇ · ( µ ∇ u ) in Ω . (2)To represent these equations in the spectral domain, we make use of u ( x , t ) = X ω i ˜ u ω i ( x ) e ˆ j ω i t , p ( x , t ) = X ω i ˜ p ω i ( x ) e ˆ j ω i t , (3)where ˆ j = √−
1. Here the frequency ω i is defined in terms of period T as ω i = π i / T for i = , , . . . , N m with N m denoting the largest computed mode. For time-periodic flows, such as those encountered in cardiovascular modeling, T is the cardiac cycle, which is the time it takes for boundary conditions to repeat themselves. For non-periodicflows, one may take T → ∞ and select it such that the results become insensitive to the choice of T . Similar tothe state variables, Dirichlet and Neumann boundary conditions can be expressed in the spectral domain as g ( x , t ) = P ω i ˜ g ω i ( x ) e ˆ jt ω i t and h ( x , t ) = P ω i ˜ h ω i ( x ) e ˆ jt ω i t , respectively. With these definitions, Eq. (2) and continuity equation canbe written as ˆ j ω i ρ ˜ u ω i = −∇ ˜ p ω i + ∇ · ( µ ∇ ˜ u ω i ) in Ω , ∇ · ˜ u ω i = Ω , ˜ u ω i = ˜ g ω i on Γ g , ( − ˜ p ω i I + µ ∇ ˜ u ω i ) · n = ˜ h ω i on Γ h . (4)These equations resemble those of the steady Stokes equation with a nonzero source term. The only di ff erence isthat the first (source) term in Eq. (4) is complex-valued unless ω i =
0, for which the steady Stokes equations in realdomain are recovered.
Since the solution to Eq. (4) at ω i and ω j for i , j are independent, we only need to provide a solution procedureat a single frequency ω i . Thus, for the sake notation brevity, we drop subscript i and denote our formulation in termsof variable ω below.The weak form of Eq. (4) is stated as follows. Given the frequency ω , find ˜ u ω ∈ S and ˜ p ω ∈ P such that for any w ∈ W and q ∈ Q B G ( w , q ; ˜ u ω , ˜ p ω ) = F G ( w , q ) , B G = Z Ω h ˆ j ωρ w · ˜ u ω + ∇ w : ( − ˜ p ω I + µ ∇ ˜ u ω ) + q ∇ · ˜ u ω i d Ω , F G = Z Γ h w · ˜ h w d Γ , (5)3olds. In Eq. (5), w and q are test functions for velocity and pressure, respectively, and S = n ˜ u ω | ˜ u ω ( x ) ∈ ( H ) n sd , ˜ u ω = ˜ g ω on Γ g o , W = n w | w ( x ) ∈ ( H ) n sd , w = on Γ g o , P = n ˜ p ω | ˜ p ω ( x ) ∈ L o , Q = n q | q ( x ) ∈ L o . (6) Denoting the finite-dimensional subspace of S , W , P , and Q by S h , W h , P h , and Q h , respectively, we attempt tosolve the discrete form of Eq. (5) using Galerkin’s approximation. Namely, we seek ˜ u h ω ∈ S h and ˜ p h ω ∈ P h such thatfor any ˜ w h ω ∈ W h and ˜ q h ω ∈ Q h B G (cid:16) w h , q h ; ˜ u h ω , ˜ p h ω (cid:17) = F G (cid:16) w h , q h (cid:17) , (7)holds. In writing Eq. (7), we assumed Ω h = Ω , i.e., the computational domain after discretization remains unchanged.If the two di ff er, then one must perform the integrals over Ω h rather than Ω when computing B G and F G in Eq. (7).The Galerkin’s formulation of incompressible flow has a saddle-point nature, which produces a singular systemif one were to adopt similar shape functions for velocity and pressure. Various techniques have been proposed toovercome this issue, ranging from mixed-element [22] and penalty techniques [23] to stabilized finite element methods[24]. In the present study, we adopt the mixed finite element method, which allows us to reduce the number of variablesthat influence the accuracy of the SCVS formulation, thus simplifying the measurement of its numerical properties.More specifically, we employ linear and quadratic shape functions for velocity and pressure, respectively, to satisfyinf-sup condition (also known as LBB condition) [25, 26, 27]. In 2D and 3D, we use linear triangular and tetrahedralelements with 3 and 4 nodal points, respectively, for pressure. For velocity, we use quadratic shape functions thatproduce 6 and 10 nodes in each of the aforementioned elements in 2D and 3D, respectively. Denoting these linearand quadratic shape functions at global node A by M A ( x ) and N A ( x ), respectively, the test functions and unknowns areinterpolated in space using w h ( x ) = X A ∈ η \ η g M A ( x ) W A , q h ( x ) = X A ∈ ˆ η N A ( x ) Q A , ˜ u h ω ( x ) = X A ∈ η \ η g M A ( x ) U A + X A ∈ η g M A ( x ) G A , ˜ p h ω ( x ) = X A ∈ ˆ η N A ( x ) P A , (8)where η , η g , and ˆ η refers to the velocity nodes, velocity nodes on the Dirichlet boundaries, and pressure nodes,respectively. U A , P A , W A , Q A in Eq. (8) are the velocity and pressure unknowns and their respective test functions. G A is the prescribed velocity defined on the Dirichlet boundaries after discretization such that˜ g h ω ( x ) = Π h ˜ g ω ( x ) = X A ∈ η g M A ( x ) G A , (9)where Π h ˜ g ω is an operator that projects ˜ g ω to the finite-dimensional discrete space.Substituting for the variables appearing in Eq. (7) using Eq. (8) while ensuring the results holds for any W A and Q A produces the following system of linear equations AX = R , (10)where A = " K DD T , X = " UP , R = " B , (11)4nd K and D matrices and B vector are computed as K AB = Z Ω (cid:16) ˆ j ωρ M A M B + µ ∇ M A · ∇ M B (cid:17) I d Ω , D AB = − Z Ω ∇ M A N B d Ω , B A = Z Γ h M A ˜ h ω d Γ − K AB G B . (12)The solution to this linear system is obtained using a brute-force approach. Namely, we treat the entire linear systemas a single sparse matrix and solve it iteratively using the Generalized minimal residual method (GMRES) method[28]. The GMRES algorithm used for this purpose is slightly modified from its real counterpart owing to K beingcomplex-valued. Provided that the left-hand side matrix is symmetric and has a block structure, one may adopt alinear solution strategy based on Schur complement to reduce the overall cost of solving this linear system in thefuture [29, 30, 31, 32]. In this section, we discuss the implementation of the SCVS algorithm. Overall, the procedure involves representingthe boundary conditions in the spectral domain, then solving for flow at each mode, and reconstructing the solution inthe real domain.1. Take the Fourier transformation of the boundary conditions g ( x , t ) and h ( x , t ) as˜ g ω i ( x ) = T Z T g ( x , t ) e − ˆ j ω i t d t , ˜ h ω i ( x ) = T Z T h ( x , t ) e − ˆ j ω i t d t . (13)2. Select the number of modes to truncate the above series based on the smoothness of the boundary conditions.More specifically, one can select a tolerance ǫ m to obtain a priori estimate of the number of required modes N m + e = k h − ˆ h k L ( Γ h × T ) k h k L ( Γ h × T ) + k g − ˆ g k L ( Γ g × T ) k g k L ( Γ g × T ) < ǫ , (14)in which ˆ h ( x , t ) = P N m ω i = ˜ h ω i ( x ) e ˆ j ω i t and ˆ g ( x , t ) = P N m ω i = ˜ g ω i ( x ) e ˆ j ω i t are the truncated Neumann and Dirichletboundary conditions using N m + k h k L ( Γ h × T ) = R Tt = P Γ i ∈ Γ g k h k L ( Γ i ) d t and k g k L ( Γ g × T ) = R Tt = P Γ i ∈ Γ g k g k L ( Γ i ) d t .3. Construct and solve the linear system in Eq. (10) for ω , ω , . . . , and ω N m given the boundary conditionscomputed from the previous step, i.e., ˜ g ω i ( x ) and ˜ h ω i ( x ).4. Given the solution U A and P A , reconstruct the solution in the spectral domain using Eq. (8) and then in time as u h ( x , t ) = N m X i = ˜ u h ω i ( x ) e ˆ j ω i t , p h ( x , t ) = N m X i = ˜ p h ω i ( x ) e ˆ j ω i t , (15)The number of computed modes N m +
1, which was selected based on a priori error estimate, can be refinedfollowing the computation of each solution mode. This way, one can compute the di ff erence between the solutionscomputed from Eq. (15) using N m − N m modes and decide whether to continue computing the next mode ifthat di ff erence is larger than a prespecified tolerance. Later in Section 3.5, we will show the rate at which a priori(Eq. (14)) and a posteriori (Eq. (16)) errors drop as N m increases are the same, implying that one can simply relyon Eq. (14) to select N m . In our experience, N m ≤ ffi cient for cardiovascular applications where theboundary conditions are smooth and periodic in time. 5 .5. Error estimation To ensure consistency and accuracy of the proposed formulation, we will consider a set of canonical cases wherean analytical solution to the governing equations is available. The overall error of a conventional scheme stemsfrom spatial discretization error, the linear solver error, time integration error, and in the case of nonlinear equations,Newton-Raphson iteration error. If we consider the solution at a single mode, i.e., boundary conditions oscillating ata single frequency, then SCVS is only prone to the first two types of error enumerated above (discretization and linearsolver). Thus, in what follows, we establish the overall accuracy of the SCVS algorithm as the mesh size or linearsolver tolerance is varied. Also, for scenarios in which the boundary conditions are generic functions of time, we willstudy the convergence of the SCVS solution as a function of N m .The overall error e ( t ) is defined based on the relative di ff erence between the computed velocity u h ( x , t ) and thereference velocity u ( x , t ) over the domain Ω as e ( t ) = k u ( x , t ) − u h ( x , t ) k L ( Ω ) k u ( x , t ) k L ( Ω ) , (16)where the L ( Ω ) norm is defined as k u k L ( Ω ) = R Ω u · u d Ω . As discussed earlier, this total error can be decomposed tothe error due to spatial discretization e H and the error due to the linear solver e L . Denoting the best approximation ofthe solution in the finite dimensional solution space S h by Π h u , L -norm of the interpolation error for the quadraticshape functions employed here can be estimated as [33] e H ( t ) = k u − Π h u k L ( Ω ) k u k L ( Ω ) C h k u k H ( Ω ) k u k L ( Ω ) , (17)where C is a constant and h is the diameter of the smallest element-bounding circle in 2D and sphere in 3D.The error related to the linear solver e L is produced by the approximate nature of the iterative solution procedureemployed in solving Eq. (10). This error can be related to the linear solver tolerance ǫ L as e L ( t ) = C k R − A b X k k R k ≤ C ǫ L , (18)in which b X is the approximate solution obtained from the iterative linear solver. In Eq. (18), the constant C dependson A − and thus is independent of the boundary conditions. How the overall error e changes with e H and e L for a pipeflow will be investigated in Section 3.4.
3. Results
The solution procedure and four case studies selected for evaluating the SCVS algorithm are discussed in thissection. To establish the relative accuracy and cost of the SCVS, we will compare it against a standard RBVMSsolver. The details pertaining to the RBVMS algorithm are included in Appendix A. The four test cases consideredin this study are (1) a 2D channel, (2) a 2D diverging nozzle, (3) a 3D cylinder, (4) a patient-specific model of Glennprocedure with an anastomosis [34]. Since the analytical solution is available for the first and third cases, these twocases will be employed to validate our implementation and evaluate its numerical characteristics. While the boundaryconditions are periodic and contain a single frequency for those two cases, they are generic functions of time for thecase two and four, thus allowing us to study the convergence of solution with regard to N m and its relationship to e M .At the end of this section, the performance of the SCVS and the RBVMS are compared in terms of total CPU timeand wall-clock time to demonstrate the performance of the SCVS. The 2D and 3D models are initially discretized using triangular and tetrahedral elements. For this purpose, we usea combination of Tetgen [35] and Simvascular [36] software. An in-house script was developed to generate quadraticelements from the linear triangular and tetrahedral elements through a node insertion process. For boundary nodes,special care was taken to ensure that inserted nodes are properly projected to be located on the curved boundaries.6 igure 1: The schematic of the simulated 2D channel flow. A cosinusoidal Neumann boundary condition is imposed at the inlet. The shown contouris normalized velocity magnitude at t = T / W = π case obtained from the SCVS. Since the SCVS simulations are performed using quadratic and the RBVMS using linear elements, we reduce thenumber of elements for the quadratic meshes such that the total number of degrees of freedom is roughly the same asthe linear meshes, thus allowing for a one-to-one comparison between the two. A zero initial condition is used for allthe RBVMS simulations. The time step size for the RBVMS is selected to su ffi ciently resolve the time variation ofboundary conditions. More specifically, we employ 2,000 time steps for cases with time-varying boundary conditions.For the steady simulations, the time step size is selected such that the di ff usive CourantFriedrichsLewy (CFL) number ∆ t ν/ h =
1. The time integration is continued for steady cases until the relative residual falls below 10 − . A total offive cycles are simulated to achieve cycle-to-cycle convergence for the RBVMS simulations. All the results shownbelow correspond to the last simulated cycle.For all cases, the Reynolds number is Re < − to ensure nonlinear e ff ects are negligible in the case of RBVMSsolver. The Reynolds number is defined based on the maximum flow rate through boundaries, the area of the boundarywith the maximum flow rate, and the kinematic viscosity ν = µ/ρ . We only perform a single Newton-Raphson iterationin the RBVMS simulations, given that the governing equations are linear.The GMRES method is adopted for solving the resulting linear system for both the SCVS and RBVMS using thesame tolerance ǫ L = − . We used our in-house linear solver for this purpose [32, 31, 37]. As discussed earlier,the linear system in the case of the SCVS is complex-valued, requiring us to make a necessary adjustment to thereal-valued version of our GMRES solver to accommodate for a complex-valued linear system. Both the RBVMSand the SCVS solvers are also implemented in our in-house finite element solver, which is written in object-orientedFortran and parallelized using the MPI library. A 2D channel with a length-to-height aspect ratio of 5 is considered as our first case study (Fig. 1). The geometryis discretized using 3,762 linear triangle elements for the RBVMS solver and 882 quadratic elements for the SCVSsolver, producing 2,000 and 1,881 nodes and 6,000 and 4,262 degrees of freedom for two cases, respectively. Thenon-slip boundary condition is imposed for the top and bottom walls. Either a nonzero constant or a cosinusoidalNeumann boundary condition is imposed on the inlet, and zero traction is imposed at the outlet. The amplitude of thewave is selected such that the Reynolds number is su ffi ciently small. The oscillation frequency ω is varied to simulatea wide range of conditions. More specifically, the Womersley number W = ω H /ν with H denoting the half-channelheight varies from 0 to 20 π for the results shown below.An analytical solution is available for the oscillatory flow in the 2D channel, expressed in the spectral domain asa function of local height y and ω according to [38]˜ u x ( y , ω ) = ˜ h x µ L ( H + y )( H − y ) , ω = , ˆ j ˜ h x ρ L ω (cid:20) − cosh ( Λ ) − cosh (cid:18) Λ yH (cid:19)(cid:21) , ω , , (19)where Λ = ˆ jW , ˜ u x is the streamwise velocity in the spectral domain, L is the channel length, and ˜ h x is the traction7 able 1: Comparison of error in the solution obtained from the SCVS and the RBVMSsolver as a function of the Womersley number W computed using Eq. (16) for the 2Dchannel flow shown in Fig. (1). The errors and relative figures are in percent. W = ω H /ν SCVS (%) RBVMS (%) Relative (%) e ( T / e ( T / e ( T / e ( T / T / T /
20 1 . × − . . × − π .
01 0 .
031 0 .
27 0 .
40 3 . . π .
12 0 .
46 0 .
84 8 . . π .
29 1 . . ω imposed at the inlet acting in the streamwise direction. This solution can also be expressed in time as u x ( y , t ) = real n ˜ u x ( y , ω ) e ˆ j ω t o . (20)The SCVS and the RBVMS simulation results are compared against the analytical solution for W = π , W = π ,and W = π at two time points t = T / t = T / W ). This larger error can be attributed to the velocity profile developing sharper gradientsnear the walls. At a similar W and time point, the SCVS provides more accurate predictions in comparison to theRBVMS. This higher accuracy is despite the fact that a larger number of degrees of freedom were employed in theRBVMS simulations. The primary reason for this improved accuracy is the use of quadratic shape function for theSCVS. The time integration scheme and the stabilization terms also contribute to the larger errors in the RBVMSresults.A more condensed version of these results is provided in Table 1, where the relative error at time point T / / e = . × − is solely due to the linear solver. Repeating this computation with a smaller ǫ L confirms that e for the SCVS canbe arbitrarily lowered without refining the grid. The steady solution for the RBVMS, on the other hand, shows arelatively larger error of 1.93% that is grid-dependent given that it utilizes linear shape functions. The same improvedaccuracy is observed at higher modes, where for instance, the RBVMS produces a solution with a 13.4% error at T / W = π , whereas the SCVS error remains as low as 1.83%. In general, in comparison to the RBVMS method,the SCVS is an order of magnitude more accurate for a mesh with a similar number of nodes. In the second example, we consider a 2D diverging nozzle with an expansion ratio of 2 (Fig. 3). This geometryis discretized using 2,400 linear triangle elements for the RBVMS solver and 600 quadratic elements for the SCVSsolver, producing 1,281 nodes in both cases with 5,124 and 4,184 degrees of freedom, respectively. A time-periodicNeumann boundary condition is imposed on the inlet, and zero traction is imposed on the outlet (Fig. 3). The timevariation of the inlet Neumann boundary is selected to resemble a physiologic pressure waveform, with its time-average being set to zero to emphasize the role of unsteady modes in the solution. The top and bottom walls areconsidered as non-slip boundaries.The inlet Neumann boundary contains a wide range of frequencies. However, to construct the SCVS solution,we consider up to N m = ω , . . . , ω . The Womersley number,i.e., W = ω H /ν with H being the half-inlet height, for these simulated modes ranges from 0 to 10 π . To show theconvergence of SCVS solution with N m , we reconstructed solution at t = T using two ( N m = N m = N m =
5) modes. The results of these computations normalized by u ∗ = | ˜ h x ( ω = | H /µ for N m =
1, 3, and 5 areshown in Fig. 4 along with the results obtained from the RBVMS. This figure shows the qualitative convergence ofthe SCVS at N m = Q / Q ∗ , with Q ∗ = Hu ∗ beingthe characteristic flow rate, is computed and shown in Fig. 5 for the RBVMS and the SVC with N m =
1, 3, and 5.Taking the flow computed from the SCVS solution with N m =
10 as the reference, the relative error in the prediction8 igure 2: Normalized streamwise velocity as a function of channel height y / H obtained from the SCVS (solid black), the RBVMS (dashed red),and the analytical solution (circles) for the 2D channel flow shown in Fig. 1. The velocity is normalized using analytical u max . The results on theleft and right columns are extracted at t = T / t = T /
2, respectively, and those on the first, second, and third row correspond to W = π , 10 π ,and 20 π , respectively. igure 3: The schematic of the 2D diverging nozzle geometry. A physiologic time-periodic Neumann boundary condition is imposed on the inlet.Figure 4: Velocity in the x -direction for 2D diverging nozzle shown in Fig. 3 at t = T computed using the RBVMS (a) and the SCVS (b-d) with N m =
1, 3, and 5, respectively. Figure 5: Normalized flow rate for the case shown in Fig. 3 obtained from the RBVMS (dashed red) and the SCVS with N m = N ele , N nds , and N dof denote the numbers of elements, nodes, and degrees of freedom, respectively. Mesh SCVS RBVMSQM1 QM2 QM3 QM4 LM1 LM2 LM3 LM4 N ele N nds N dof N m =
1, 3, and 5 is 25%, 5.0%, and 1.5%,respectively. The fast rate of convergence of the SCVS with respect to N m confirms our earlier argument that onlya few modes are needed to obtain a converged solution using the SCVS algorithm when boundary conditions varysmoothly in time. For the third case study, we consider an oscillatory laminar pipe flow. A cylinder with an aspect ratio of L / R = W = ω R /ν = , π, . . . , π with R denoting the radiusof the pipe. As we discussed earlier, the domain is discretized in space using linear and mixed quadratic-lineartetrahedral elements for the RBVMS and the SCVS simulations, respectively. A wide range of meshes is constructedto investigate the convergence of both solvers (Table 2). The number of elements in corresponding meshes is selectedsuch that the number of nodes is roughly the same for both solvers. In overall, the element size normalized by thepipe radius varies from 0.0340 to 0.228 among simulated cases.Similar to the first case considered above, an analytical solution is available for this case that permits us to establish11 igure 6: Schematic of the simulated oscillatory laminar flow in a pipe. The contour of normalized radial velocity magnitude is shown for W = π case at t = T / the accuracy of each solver. For an oscillatory flow in a pipe, the solution in the spectral domain is expressed as [39]˜ u x ( r , ω ) = ˜ h x µ L ( R − r ) , ω = , ˆ j ˜ h x ρ L ω (cid:20) − J ( Λ ) − J ( Λ rR ) (cid:21) , ω , , (21)where Λ = − ˆ jW and J is the zero order Bessel function of the first kind. This solution can be expressed in timeusing u x ( r , t ) = real n ˜ u x ( r , ω ) e ˆ j ω t o . (22)For a qualitative evaluation of the SCVS and the RBVMS solutions, simulated u x ( r , T /
4) and u x ( r , T /
2) along withthe analytical prediction of Eq. (22) at W = π , 40 π , and 80 π are shown in Fig. 7. These results are obtained using thefinest grids, namely QM4 and LM4 in Table 2. Similar to what we observed earlier in Section 3.2 for the 2D channelflow, the SCVS provides a better approximation, particularly at larger Womersley numbers or at t = T /
2, in whichthe solution exhibits sharper gradients. The error in the steady state solution ( W =
0) obtained from the RBVMSand SCVS solver using the finest grids is e = . × − and e = . × − , respectively. The superior accuracy ofthe SCVS, which is achieved despite using fewer degrees of freedom, is primarily a result of utilizing a higher-ordershape function. The numerical integration and stabilization terms are the secondary contributors to the larger error inthe case of the RBVMS.To study mesh convergence, the error e ( T /
2) defined in Eq. (16) is computed for the SCVS and the RBVMS as afunction of h at W = π (Fig. 8). Third-order accuracy is observed for the SCVS, indicating that the measured erroris dominated by the discretization error e H . This third order accuracy is in agreement with the estimate in Eq. (17),which also predicted e H ∝ h . A similar relationship can be obtained for the RBVMS, which utilizes linear shapefunctions, as e H C h k u k H ( Ω ) k u k L ( Ω ) . (23)The second order accuracy of RBVMS is also confirmed by the results shown in Fig. 8.To further analyze each method’s accuracy for various flow conditions, we next study the overall error as a functionof the Womersley numbers W . As we observed earlier, the velocity profile develops sharper gradients at higherfrequencies, thereby producing larger error at higher W . The same observation also holds in this case, where e ( T / W (Fig. 9).Since the linear solver tolerance ǫ L is selected to be su ffi ciently small in this case and e L is negligible in comparisonto e H , we can explain e ∝ W / for the SCVS and e ∝ W for the RBVMS by further analyzing the behavior of e H foreach solver. To utilize Eqs. (17) and (23) for this purpose, we first need to establish how k u k L , k u k H , and k u k H vary12 igure 7: Normalized velocity profiles as a function of radius for an oscillatory laminar pipe flow (shown in Fig. 6) predicted using the SCVS (solidblack), the RBVMS (dashed red), and analytical solution (circles). The results on the left and right columns are extracted at t = T / t = T / W = π , 40 π , and 80 π , respectively. Computations are performed on QM4and LM4 grids for the SCVS and the RBVMS, respectively. All the results are normalized using the maximum velocity from the analytical solution u max . .03 0.06 0.09 0.12 0.1510 -3 -2 -1 RBVMSSCVS
Figure 8: The relative error e (defined in Eq. (16)) at t = T / W = π .The shown data correspond to the SCVS (solid black circles), the RBVMS (hollow red circles), the analytical estimate of the SCVS error fromEq. (17) (solid black line), and the analytical estimate of the RBVMS error from Eq. (23) (dashed red line).
50 150 25010 -3 -2 -1 RBVMSSCVS
Figure 9: The relative error e (defined in Eq. (16)) at t = T / W for the oscillatory pipe flow shown in Fig. 6.The shown results correspond to the SCVS using QM4 grid (solid black circles), the RBVMS using LM4 grid (hollow red circles), the analyticalestimate of the SCVS error from Eqs. (17) and (26) (solid black line), and analytical estimate of the RBVMS error from Eqs. (23) and (26) (dashedred line). W . Considering only cases in which W , u x = u x ( r , t ), from Eq. (21) we can write k u k L = Z Ω u x d Ω ! / = ˜ h x ρω Z Rr = π rz ¯ z d r ! / , k u k H = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ u x ∂ r (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L = ˜ h x ρ w Λ R ! Z Rr = π rz ¯ z d r ! / , k u k H = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ u x ∂ r (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L = ˜ h x ρ w Λ R ! Z Rr = π rz ¯ z d r ! / , (24)where z ( r ; Λ ) = − J ( Λ ) − J ( Λ rR ), z ( r ; Λ ) = h J ( Λ rR ) − J ( Λ i rR ) i / (2 J ( Λ )), z ( r ; Λ ) = h J ( Λ rR ) − J ( Λ rR ) i / (4 J ( Λ )), and ¯ z denotes complex conjugate of z . Neglecting the dependence ofthese norms on the integrals of z , z , and z , it is straightforward to show that k u k L ∝ ˜ h x R µ W , k u k H ∝ ˜ h x R µ , k u k H ∝ ˜ h x R µ W , (25)indicating k u k L , k u k H , and k u k H should vary proportional to W − , W , and W / . Comparing these approximatepredictions against reference quantities (obtained from 1D numerical integration of Eq. (24)) show that all exponentsare over-predicted by 0.2 (Fig. 10). The ratio of these norms, however, is correctly predicted using Eq. (25) as k u k H k u k L ∝ R − W , and k u k H k u k L ∝ R − W / , (26)indicating that the RBVMS and the SCVS error should grow proportional to W and W / at su ffi ciently small ǫ L owingto Eqs. (17) and (23), respectively. These analytical predictions are in agreement with our numerical results that wereshown earlier in Fig. (9).Next, we investigate the e ff ect of linear solver tolerance ǫ L on the overall accuracy of the SCVS, utilizing QM4mesh and the solution corresponding to W = π evaluated at t = T /
2. The left and right tail of the results shown inFig. 11 indicate that the overall error e is independent of e L for su ffi ciently small e L and proportional to ǫ L when e L is su ffi ciently large. This observation confirms that (as expected) the overall error is dominated by the discretizationerror e H and that of the linear solver e L at small and large ǫ L , respectively. The transition between dominance of e H and e L occurs at approximately ǫ L = − , which is a value specific to this case study and a function h among others.Note that this tolerance is the optimal tolerance with regard to the computational cost, as further decrease in ǫ L is notmet with an improved overall solution accuracy. Lastly, note that the slope of 1 observed on the right tail of Fig. 11 isin agreement with Eq. (18), which predicts e L ∝ ǫ L .If we take the main objective of designing a computational method to be making a viable prediction at a specifiedaccuracy with the lowest cost possible, then it is desirable to compare the SCVS and the RBVMS methods on anerror-versus-cost plot. Thus, we have extracted the cost (in terms of CPU time) and error (in terms of e ( T / W . The results are shown in Fig. 12. Comparing thesymbols with the same color, i.e., for a flow at a given W , the SCVS always provides more accurate predictions at alower cost. Focusing on W = π cases (blue symbols), the SCVS yields a similar accuracy at the coarsest grid (solidtriangle) to that of the RBVMS at the finest grid (hollow circle) while reducing the cost by three orders of magnitude.If one holds the computational cost relatively similar (solid blue circle versus hollow blue triangle, for instance), theSCVS o ff ers over an order of magnitude improvement in accuracy relative to the RBVMS. For the last test case, we will consider a complex patient-specific geometry acquired from a patient undergoingGlenn operation [34]. As shown in Fig. 13, the superior vena cava (SVC) is anastomosed to the left pulmonary artery15 -4 -3 -2 -1 Figure 10: L (circles), H (squares), and H (triangles) norms of the analytical solution u for an oscillatory flow in cylinder as a functionWomersley number W . The symbols are obtained by 1D numerical integration of Eq. (24) whereas the lines with given slopes are curve fits. -8 -6 -4 -3 -2 -1 Figure 11: The overall relative error in the SCVS solution e ( T /
2) as a function of the linear solver tolerance ǫ L . The case considered herecorresponds to an oscillatory pipe flow shown in Fig. 6 at W = π and using QM4 mesh. A line with a slope 1 is shown for the reference. -4 -3 -2 -1 Figure 12: Relative error e ( T /
2) as a function of the computational cost t C for the pipe flow example shown in Fig. 6. Symbols with the same colorshould be compared. The solid symbols correspond to the SCVS and hollow to the RBVMS. The circle, square, diamond, and triangle symbolscorrespond to QM4 / LM4, QM3 / LM3, QM2 / LM2, and QM1 / LM1 meshes listed in Table 2. The black, red, and blue colors correspond to W = π ,40 π , and 80 π , respectively. These results show that the SCVS, compared to the RBVMS, always provides a higher accuracy at a lower cost. (LPA) and right pulmonary artery (RPA) in this operation. The geometry is discretized using 988,747 linear tetrahedralelements for the RBVMS solver and 132,066 quadratic tetrahedral elements for the SCVS. This discretization resultsin 163,791 and 183,708 nodes and 655,164 and 574,401 degrees of freedom for the linear and quadratic meshes,respectively. The boundary conditions are selected based on physiologic data and a Windkessel model [3, 40] for theRPA and LPA faces (Fig. 13). The wall is considered as a non-slip boundary.In this example, twelve modes ( N m =
11) are simulated in total for the SCVS. The case with N m =
11 is usedas the reference to measure error for the remaining cases, given that the solution accuracy of the SCVS was superiorto that of the RBVMS based on the cases discussed earlier. The Womersley number W = ω D /ν is based on thehydraulic diameter of the SVC ( D h ), which is defined such that π D / W = π , 293 π , . . . , 1611 π for simulated ω , ω , ω , . . . , ω . The results of these computations are normalized by u ∗ = | ˜ h SVC ( ω = | D h /µ and Q ∗ = π D h u ∗ / W increases (Fig. 14). While the velocity peaks at the center of the vessels for the steady solution (Fig. 14-(a)), itdevelops two peaks closer to the walls at a larger W . This behavior resembles the trend that we observed earlier in thecanonical geometries (e.g., cylinder in Section 3.5). This behavior is a result of a change in the relative importance ofvarious terms that appear in the Stokes equation. While only the viscous and pressure terms are active at W =
0, threeterms (acceleration, pressure, and viscous) balance each other for W >
0. As W increases, the relative importance ofthe acceleration term increases, leading to a solution that exhibits peaks near the walls rather than at the center of thegeometry.The SCVS solution in the time domain is reconstructed using Eq. (15) and qualitatively compared against thatof the RBVMS (Fig. 15). This qualitative comparison shows that a few modes (in this case 6, including the steadysolution) are su ffi cient for resolving the solution in this complex geometry.For a more quantitative comparison between the SCVS and the RBVMS solutions, the predicted flow throughthree branches is shown in Fig. 16. This figure confirms that the di ff erence between the two solvers’ solution reducesas more modes are employed. Using the case with N m =
11 as a reference, the relative errors of flux computed fromthe RBVMS are 0 . . .
34% on the SVC, LPA, and RPA, respectively. The relative errors in the SCVS17 igure 13: The geometry and boundary conditions employed for the patient-specific case study. h max denotes the maximum value of tractionimposed on the SVC.Figure 14: The velocity magnitude normalized by u ∗ obtained from the SCVS solver at t = T / W =
0, (b) W = π , (c) W = π , (d) W = π , (e) W = π , (f) W = π , (g) W = π , (h) W = π , (i) W = π , (j) W = π , (k) W = π , (l) W = π . igure 15: Normalized velocity magnitude at t = / T obtained from the SCVS with N m = Q / Q ∗ through the SVC (a), LPA (b), and RPA (c) predicted by the RBVMS (dashed red) and the SVC with N m = N m = N m = solution with N m = . . .
12% for the SVC, LPA, and RPA, respectively. Note that, due to theerror in the RBVMS solution (as well as discretization error in the SCVS), the two solutions do not converge even atvery large N m . Nevertheless, the SCVS does converge to a solution as N m → ∞ . The rate of this convergence dependson the smoothness of the boundary conditions.The relationship between the smoothness of the boundary condition time variation and the convergence rate ofthe SCVS with respect to N m is demonstrated in Fig. 17. e M from Eq. (14), which measures the truncation errorin the imposed boundary conditions associated with the finite N m , follows the same trend as the overall error e . Asmore modes are included in the solution, the boundary conditions approach their temporal reference profile, therebyimproving the overall accuracy of the SCVS solution. Note that the two trends slightly diverge on the right tail of thisplot. This divergence is an artifact of taking the case with N m =
11 as the reference when measuring k e k , whereas e M is measured against the reference solution. If the exact solution were available and employed as the reference solution,then k e k for very large N m would have converged to the larger of discretization e H and linear solver e L errors ratherthan going to zero. The primary reason for the development of the SCVS algorithm was to obtain a more computationally e ffi cientsolution procedure. We discussed some of the performance metrics of the SCVS algorithm in the discussion pertaining19 -4 -3 -2 -1 Figure 17: The overall error k e ( t ) k L ([0 , T ]) (solid circles) as a function of N m for the case shown in Fig. 13, using the case with N m =
11 as thereference. The truncation error for the imposed Neumann boundaries e M (red squares) is also shown, which follows the same trend as k e k . to Fig. 12. In this section, we provide a more comprehensive account of the performance of these two solvers in termsof computational costs as well as wall-clock-time for all cases studied above.For cases with multiple grids, we considered the finest grid in computing performance values. Also, we consider N m = ffi cient for su ffi ciently smooth boundary conditions. Thus, the computationalcost t C for the SCVS is computed as the cumulative cost of simulating the first 6 modes. Its wall-clock-time t W , on theother hand, is taken as the maximum of the wall-clock-time of those 6 simulated modes since they are embarrassinglyparallelizable. The simulation results presented in Sections 3.2 and 3.4 for the canonical 2D channel and 3D cylindercases were for individual modes. Nevertheless, we consider the sum of the cost of the first six modes for those casesas well. This choice was made to represent a situation in which the flow in these models is simulated with an arbitraryboundary condition rather than one with a unimodal oscillation. Furthermore, note that both t C and t W for the RBVMSscale linearly with the total number of time steps . The results reported here are what is considered typical, namelysimulating five cycles with each including 2,000 time steps for a total of 10,000 time steps. Finally, the performancemetrics are in general implementation- and machine-dependent. However, given that both the RBVMS and the SCVSsolvers are implemented by our group, written and compiled using the same language and compiler, linked againstthe same libraries, and ran on the same machine, the figures in Tables 3 provide an apple-to-apple comparison of thetwo formulations.As was hypothesized earlier, the SCVS reduces the overall computational cost significantly when it is comparedagainst the standard RBVMS formulation. This reduction in cost that ranges from less than 1% for 2D cases to 11%for the 3D Glenn geometry is primarily achieved by reducing the number of linear solves. While this number is N m + = = , / ffi ce for solving the linear system obtained for the RBVMS formulation. The reason for the large numberof iterations for the SCVS is mostly attributed to the high condition number of A matrix that has a zero diagonalsub-block (Eq. (11)). As discussed in Section 4, reducing the cost of solving this linear solver presents an opportunityfor significantly reducing the cost of the SCVS algorithm in the future. This statement applies only to the linear equations considered here, since a large ∆ t increases the number of Newton-Raphson iterations athigher Reynolds numbers, thus changing the cost of advancing the solution for a time step. able 3: Comparison of the computational performance of the SCVS and RBVMS solvers. N p , t C , t W denote thenumber of processors, computational costs in CPU-hours, and the simulation wall-clock-time in hours, respectively.The last two columns are t C and t W for the SCVS relative to those of the RBVMS. All the SCVS results correspondto N m =
5. For cases with multiple grids, we adopted the finest grid in computing these performance figures. Therelative figures are in percent.
Case SCVS RBVMS Relative (%) N p t C (hr) t W (hr) N p t C (hr) t W (hr) t C t W
2D channel 6 1 . × − . × −
16 1.9 0.12 0.79 3.42D nozzle 6 7 . × − . × −
16 1.0 0.064 0.74 2.93D cylinder 384 2 . × . ×
72 5.1 1.33D Glenn 768 3 . × . ×
23 11 2.2The performance gap between the SCVS and RBVMS widens when wall-clock time is concerned as t W for theSCVS is less than 4% of that of the RBVMS for all cases considered here (Table 3). As we discussed earlier, themodes in the SCVS algorithm are independent and can be solved concurrently. Therefore, the number of processorsutilized in the SCVS computation can far exceed that of the RBVMS without loss of parallel e ffi ciency. Di ff erentlyput, the wall-clock time for the SCVS is roughly independent of the number of computed modes as long as there issu ffi cient parallel computational resources available. The cumulative e ff ect of this added dimension for parallelizationand lower overall computational cost is a formulation that allows for a much faster turn-around time for a givenproblem.
4. Future work
The biggest hurdle to be overcome in the future is the extension of the SCVS formulation to high Reynolds numberflows. Such an extension will not be trivial, as one can not exploit the linearity of the governing equations to solve forvelocity and pressure at di ff erent modes independently.The e ffi ciency of the SCVS algorithm can be significantly improved in the future by making minor adjustmentsin its formulation or its underlying linear solver. Relaxing the incompressibility constraint by using a penalty methodor including stabilization terms in its formulation is expected to reduce the condition of the underlying linear systemsignificantly. Experimenting with more e ffi cient linear solver algorithms (e.g., bi-partitioned method [32]) as well aspreconditioners can lead to additional improvements in the overall performance of the SCVS algorithm in the future.
5. Conclusions
In this paper, we proposed the SCVS as an alternative approach for fast simulation of time-periodic flow at lowReynolds numbers in complex geometries. Starting from the unsteady Stokes equation, we showed that it could beexpressed as a steady Stokes equation with an imaginary source term in the time-spectral domain. The resultingequation can then be discretized and solved as a boundary value problem at a few selected modes, avoiding the use ofcostly and unscalable time integration schemes. As a proof of concept, we showed how this boundary value problemcould be solved using Galerkin’s formulation with mixed-elements. We later employed this formulation for simulatingflow in a variety of 2D and 3D geometries. To provide a point of reference, all these simulations were also performedusing the standard RBVMS formulation.For cases with an analytical solution available, the SCVS showed about an order of magnitude improvementin accuracy relative to the RBVMS at a similar number of degrees of freedom. This di ff erence was attributed tothe use of quadratic shape functions for the SCVS and to a lesser extent to the lack of stabilization terms and timeintegrator in our formulation. The variation of the overall error in the SCVS solution as a function of grid size, linearsolver tolerance, and the mode number was measured through our numerical test cases and were shown to followanalytically-derived power-law formulas. For cases with boundary conditions changing arbitrarily in time, we showedthe overall error follows the same trend as the truncation error associated with simulating a finite number of modes.21or smooth time-varying boundary conditions, such as those encountered in cardiovascular flows, we showed thatusing as few as 6 modes is su ffi cient for reducing the overall error to O (10 − ).The SCVS led to significant improvement in performance in comparison to the RBVMS. The total computationalcost of the SCVS for the cases considered in this study varied between 0.74% and 11% of that of the RBVMS, whereasits wall-clock time was consistently lower than 4% of that of the RBVMS. Further improvement in performancethrough the use of stabilized or penalty formulation as well as the extension to higher Reynolds number remain to beexplored in the future. Appendix A. RBVMS formulation
A brief account of the RBVMS formulation employed in this study is provided below. A more detailed descriptionof this algorithm is provided in reference [6, 7, 32].The weak formulation of the RBVMS is stated as follows. Find u ∈ S and p ∈ P , such that for all w ∈ W and q ∈ Q B G ( w , q ; u , p ) + B S ( w , q ; u , p ) = F ( w , q ) , (A.1)where B G = Z Ω (cid:2) ρ w · ( ˙ u + u · ∇ u ) + ∇ w : ( − p I + µ ∇ s u ) + q ∇ · u (cid:3) d Ω B S = X e ∈ I e Z Ω e h ρ ∇ w : (cid:16) ¯ τ u p ⊗ ( u p · ∇ u ) − u ⊗ u p + τ C ∇ · uI (cid:17) + ρ w · (cid:16) u p · ∇ u (cid:17) − u p · ∇ q i d Ω , F = Z Γ h w · h d Γ . (A.2)In this equation, B G contains Galerkin’s term whereas B S are the stabilization added to allow for equal-order veloc-ity and pressure functions and prevent convective instability associated with Galerkin’s method. Other parametersappearing in (A.2) are defined as u p = − τ M ˙ u + u · ∇ u + ρ ∇ p − µρ ∇ u − f ! ,τ M = c ∆ t ! + u · ξ u + c µρ ! ξ : ξ − , ¯ τ = (cid:16) u p · ξ u p (cid:17) − ,τ C = (cid:2) tr ( ξ ) τ M (cid:3) − , (A.3)in which c = c = ξ ∈ R n sd × R n sd is covariant tensor obtained from a mappingbetween the physical and parent domains, and ∆ t is the time step size. The resulting equations are discretized usingtriangular and tetrahedral elements in 2D and 3D. To relate ˙ u to u and integrate the resulting equations in time, thesecond order generalized- α method is adopted [41]. This algorithm, which is implicit in time, requires linearizationof the nonlinear terms in Eq. (A.2). Thus, we employ the Newton-Raphson method for linearization of the resultingequations and a predictor-corrector procedure to converge on the solution at the next time step. The spectral radiusof the infinite time step, which appears in the generalized- α time integration scheme, is set to 0.2 in this study. Weverified that this parameter has no discernible influence on the results reported in this study. The resulting linearsystem is solved iteratively using the GMRES technique. Owing to the last term in B G , the sub-block matrix thatrelates pressure to continuity is not zero in this formulation (in contrast to Eq. (11)), thereby allowing for fastersolution of the resulting linear system. Overall, the RBVMS algorithm requires three nested iteration loops: the outerloop for time-stepping, the first inner loop for Newton-Raphson iterations, and the innermost iterative loop for thelinear solver. For the problems investigated in this study, however, the nonlinear terms in Eq. (A.2) can be neglected,resulting in two nested loops for the linear solver and time stepping.22 eferences [1] A. Blalock, H. B. Taussig, The surgical treatment of malformations of the heart: in which there is pulmonary stenosis or pulmonary atresia,Journal of the American Medical Association 128 (3) (1945) 189–202.[2] A. L. Marsden, M. Esmaily-Moghadam, Multiscale modeling of cardiovascular flows for clinical decision support, Applied MechanicsReviews 67 (3).[3] M. Esmaily Moghadam, F. Migliavacca, I. E. Vignon-Clementel, T.-Y. Hsia, A. L. Marsden, Optimization of shunt placement for the norwoodsurgery using multi-domain modeling, Journal of biomechanical engineering 134 (5).[4] K. Pekkan, B. Whited, K. Kanter, S. Sharma, D. De Zelicourt, K. Sundareswaran, D. Frakes, J. Rossignac, A. P. Yoganathan, Patient-specificsurgical planning and hemodynamic computational fluid dynamics optimization through free-form haptic anatomy editing tool (SURGEM),Medical & biological engineering & computing 46 (11) (2008) 1139–1152.[5] D. A. Steinman, Y. Hoi, P. Fahy, L. Morris, M. T. Walsh, N. Aristokleous, A. S. Anayiotos, Y. Papaharilaou, A. Arzani, S. C. Shadden, et al.,Variability of computational fluid dynamics solutions for pressure and flow in a giant aneurysm: the ASME 2012 Summer BioengineeringConference CFD Challenge, Journal of biomechanical engineering 135 (2).[6] A. N. Brooks, T. J. Hughes, Streamline upwind / Petrov-Galerkin formulations for convection dominated flows with particular emphasis on theincompressible Navier-Stokes equations, Computer methods in applied mechanics and engineering 32 (1-3) (1982) 199–259.[7] Y. Bazilevs, V. Calo, J. Cottrell, T. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for large eddysimulation of incompressible flows, Computer methods in applied mechanics and engineering 197 (1-4) (2007) 173–201.[8] M.-C. Hsu, Y. Bazilevs, V. M. Calo, T. E. Tezduyar, T. J. Hughes, Improving stability of stabilized and multiscale formulations in flowsimulations at small time steps, Computer Methods in Applied Mechanics and Engineering 199 (13-16) (2010) 828–840.[9] J. Hsu, A. Jameson, An implicit-explicit hybrid scheme for calculating complex unsteady flows, in: 40th AIAA Aerospace Sciences Meeting& Exhibit, 2002, p. 714.[10] Z. Zhang, F. Liu, D. M. Schuster, Calculations of unsteady flow and flutter by an Euler and integral boundary-layer method on Cartesiangrids, in: 22nd Applied Aerodynamics Conference and Exhibit, 2004, p. 5203.[11] H. M. Blackburn, S. J. Sherwin, Formulation of a Galerkin spectral element–Fourier method for three-dimensional incompressible flows incylindrical geometries, Journal of Computational Physics 197 (2) (2004) 759–778.[12] C. H. Amon, Spectra element-Fourier method for transitional flows in complex geometries, AIAA journal 31 (1) (1993) 42–48.[13] G. Karniadakis, Spectral element-Fourier methods for incompressible turbulent flows, Computer Methods in Applied Mechanics and Engi-neering 80 (1-3) (1990) 367–380.[14] M. Rosenfeld, Utilization of Fourier decomposition for analyzing time-periodic flows, Computers & fluids 24 (4) (1995) 349–368.[15] K. C. Hall, J. P. Thomas, W. S. Clark, Computation of unsteady nonlinear flows in cascades using a harmonic balance technique, AIAAjournal 40 (5) (2002) 879–886.[16] M. McMullen, A. Jameson, J. Alonso, Acceleration of convergence to a periodic steady state in turbomachinery flows, in: 39th aerospacesciences meeting and exhibit, 2001, p. 152.[17] M. S. McMullen, The application of non-linear frequency domain methods to the Euler and Navier-Stokes equations, Citeseer, 2003.[18] A. Gopinath, A. Jameson, Time spectral method for periodic unsteady computations over two-and three-dimensional bodies, in: 43rd AIAAaerospace sciences meeting and exhibit, 2005, p. 1220.[19] M. Esmaily-Moghadam, I. E. Vignon-Clementel, R. Figliola, A. Marsden, A modular numerical method for implicit 0D /
3D coupling incardiovascular finite element simulations, Journal of Computational Physics 244 (2013) 63–79.[20] C. A. Taylor, C. Figueroa, Patient-specific modeling of cardiovascular mechanics, Annual review of biomedical engineering 11 (2009) 109–134.[21] Y. Bazilevs, K. Takizawa, T. E. Tezduyar, Computational fluid-structure interaction: methods and applications, John Wiley & Sons, 2013.[22] F. Brezzi, K.-J. Bathe, A discourse on the stability conditions for mixed finite element formulations, Computer methods in applied mechanicsand engineering 82 (1-3) (1990) 27–57.[23] T. J. Hughes, W. K. Liu, A. Brooks, Finite element analysis of incompressible viscous flows by the penalty function formulation, Journal ofcomputational physics 30 (1) (1979) 1–60.[24] T. E. Tezduyar, Stabilized finite element formulations for incompressible flow computations, in: Advances in applied mechanics, Vol. 28,Elsevier, 1991, pp. 1–44.[25] O. A. Ladyzhenskaya, The mathematical theory of viscous incompressible flow, Vol. 2, Gordon and Breach New York, 1969.[26] I. Babuˇska, Error-bounds for finite element method, Numerische Mathematik 16 (4) (1971) 322–333.[27] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers, Publicationsmath´ematiques et informatique de Rennes (S4) (1974) 1–26.[28] Y. Saad, M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM Journal onscientific and statistical computing 7 (3) (1986) 856–869.[29] P. M. Gresho, R. L. Sani, Incompressible flow and the finite element method. volume 1: Advection-di ff usion and isothermal laminar flow.[30] M. A. Olshanskii, Y. V. Vassilevski, Pressure Schur complement preconditioners for the discrete Oseen problem, SIAM Journal on ScientificComputing 29 (6) (2007) 2686–2704.[31] M. Esmaily-Moghadam, Y. Bazilevs, A. L. Marsden, A new preconditioning technique for implicitly coupled multidomain simulations withapplications to hemodynamics, Computational Mechanics 52 (5) (2013) 1141–1152.[32] M. Esmaily-Moghadam, Y. Bazilevs, A. L. Marsden, A bi-partitioned iterative algorithm for solving linear systems arising from incompress-ible flow problems, Computer Methods in Applied Mechanics and Engineering 286 (2015) 40–62.[33] S. Brenner, R. Scott, The mathematical theory of finite element methods, Vol. 15, Springer Science & Business Media, 2007.[34] G. Arbia, C. Corsini, M. E. Moghadam, A. L. Marsden, F. Migliavacca, G. Pennati, T.-Y. Hsia, I. E. Vignon-Clementel, M. O. C. H. A. M.Investigators, et al., Numerical blood flow simulation in surgical corrections: what do we need for an accurate analysis?, journal of surgicalresearch 186 (1) (2014) 44–55.
35] H. Si, TetGen, a Delaunay-based quality tetrahedral mesh generator, ACM Transactions on Mathematical Software (TOMS) 41 (2) (2015)1–36.[36] A. Updegrove, N. M. Wilson, J. Merkow, H. Lan, A. L. Marsden, S. C. Shadden, SimVascular: an open source pipeline for cardiovascularsimulation, Annals of biomedical engineering 45 (3) (2017) 525–541.[37] M. Esmaily-Moghadam, Y. Bazilevs, A. Marsden, Impact of data distribution on the parallel performance of iterative linear solvers withemphasis on CFD of incompressible flows, Computational Mechanics 55 (1) (2015) 93–103.[38] C. Loudon, A. Tordesillas, The use of the dimensionless Womersley number to characterize the unsteady nature of internal flow, Journal oftheoretical biology 191 (1) (1998) 63–78.[39] J. R. Womersley, Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known, TheJournal of physiology 127 (3) (1955) 553.[40] I. E. Vignon-Clementel, C. Figueroa, K. Jansen, C. Taylor, Outflow boundary conditions for 3D simulations of non-periodic blood flow andpressure fields in deformable arteries, Computer methods in biomechanics and biomedical engineering 13 (5) (2010) 625–640.[41] K. E. Jansen, C. H. Whiting, G. M. Hulbert, A generalized- α method for integrating the filtered Navier–Stokes equations with a stabilizedfinite element method, Computer methods in applied mechanics and engineering 190 (3-4) (2000) 305–319.method for integrating the filtered Navier–Stokes equations with a stabilizedfinite element method, Computer methods in applied mechanics and engineering 190 (3-4) (2000) 305–319.