A note on the accuracy of the generalized-$\alpha$ scheme for the incompressible Navier-Stokes equations

We investigate the temporal accuracy of two generalized-$\alpha$ schemes for the incompressible Navier-Stokes equations. The conventional approach treats the pressure with the backward Euler method while discretizing the remainder of the Navier-Stokes equations with the generalized-$\alpha$ method. We developed a suite of numerical codes using inf-sup stable higher-order non-uniform rational B-spline (NURBS) elements for spatial discretization. In doing so, we are able to achieve very high spatial accuracy and, furthermore, to perform temporal refinement without consideration of the stabilization terms, which can degenerate for small time steps. Numerical experiments suggest that only first-order accuracy is achieved, at least for the pressure, in this aforesaid approach. Evaluating the pressure at the intermediate time step recovers second-order accuracy, and the numerical implementation, in fact, becomes simpler. Therefore, although the pressure can be viewed as a Lagrange multiplier enforcing the incompressibility constraint, its temporal discretization is not independent and should be subject to the generalized-$\alpha$ method in order to maintain second-order accuracy of the overall algorithm.


Introduction
The generalized-α method was initially proposed for structural dynamics [4] and possesses the desirable numerical properties of implicit schemes, as noted by Hilber and Hughes [12]. The method was later extended to address first-order systems, with an important example being the Navier-Stokes equations [16]. Interestingly, when used in conjunction with the first-order structural dynamics, the generalized-α method has been found not to suffer from the 'overshoot' phenomenon [17], making it an ideal time integration algorithm for low-frequency, inertial type problems. Over the years, this temporal scheme has gained popularity as an implicit scheme for fluid and fluidstructure interaction (FSI) problems. When it was originally introduced to computational fluid dynamics (CFD) problems, Jansen et al. considered the compressible Navier-Stokes equations written in the pressure primitive variables, which were evaluated at the intermediate time step [16]. Later, when dealing with the incompressible Navier-Stokes equations, the time-stepping scheme was used in a slightly different manner: the velocity alone was evaluated at the intermediate time step, while the pressure was treated by the backward Euler scheme. This choice of temporal discretization has in fact become quite popular in both CFD and FSI simulations (see e.g. [1,2,3,5,8,9,19,23]). Indeed, the pressure can be viewed as a Lagrange multiplier enforcing the divergence-free constraint and should therefore "not be subjected to a time integration scheme" [6]. While this may explain the widespread adoption of the backward Euler method for pressure, we note that there exist a few works evaluating the pressure at the intermediate time step [10,15,18,22]. In our view, treating the pressure by the backward Euler method is unnatural and inefficient for the following reasons. First, given that the incompressible Navier-Stokes equations constitute a limiting case of their compressible counterparts, treating the pressure by the backward Euler method is evidently inconsistent with the scheme proposed in [16]. The inconsistency is even more apparent when we recall that the pressure is the hydrostatic part of the Cauchy stress, suggesting that the hydrostatic and deviatoric parts of the Cauchy stress are evaluated differently in time in this conventional approach. Second, in the implementation of CFD algorithms, variables are typically stored in an 'interlaced' pattern, which is cache-friendly for high-performance computing [11]. Using two distinct solution updating strategies in the time integration algorithm necessitates separate extraction of the velocity and pressure degrees-of-freedom from the interlaced solution vector, thereby reducing the overall algorithmic efficiency. In addition to the inconsistency and inefficiency, we may also reasonably expect the backward Euler method to yield a loss of numerical accuracy. Second-order accuracy is one of the most important numerical attributes for a desirable implicit scheme [12], since the error accumulation from first-order schemes may be excessive for practical applications.
In this work, we compare the numerical accuracy of two generalized-α schemes. The first, which we denote as Scheme-1, is the popular choice in which the pressure is integrated by the backward Euler method; the second, denoted as Scheme-2, evaluates the pressure at the intermediate time step. We present the accuracy of quantities of interest in two benchmark cases. The Navier-Stokes equations are discretized in space by inf-sup stable NURBS element pairs. Our rationale behind this choice is twofold. We want to avoid invoking the stabilized formulation, which can degenerate for small time steps when used without special treatment [13] and is therefore not robust for temporal convergence studies. Furthermore, the robustness and exact geometric representation made possible by the NURBS-based isogeometric analysis technique allow us to eliminate the possibility of confounding errors from the spatial discretization and thereby focus on the temporal accuracy. wherein In the above, ρ is the fluid density, v is the velocity field, σ is the Cauchy stress, f is the body force, µ is the dynamic viscosity, ε is the rate-of-strain tensor, p is the pressure, and I is the second-order identity tensor. Both ρ and µ are assumed constant, and the kinematic viscosity is defined as ν := µ/ρ. Given a divergence-free velocity field v 0 , the initial condition is given by v(·, 0) = v 0 (·), inΩ.
The boundary Γ can be partitioned into two non-overlapping subdivisions, that is, Γ = Γ g ∪ Γ h and ∅ = Γ g ∩ Γ h . The subscripts g and h indicate the Dirichlet and Neumann partitions. The unit outward normal vector to Γ is denoted as n. Given the Dirichlet values g and the boundary traction h, the initial-boundary value problem is to solve for v and p that satisfy (1), (2), (4), and

Spline Spaces on the Parametric Domain
We start by reviewing the construction of B-splines and NURBS basis functions. Given the polynomial degree p and the dimensionality of the B-spline space n, the knot vector can be represented by Ξ := {ξ 1 , ξ 2 · · · , ξ n+p+1 }, wherein 0 = ξ 1 ≤ ξ 2 ≤ · · · ≤ ξ n+p+1 = 1. The B-spline basis functions of degree p, denoted as N p i (ξ), for i = 1, · · · , n, can then be defined recursively from the knot vector using the Cox-de Boor recursion formula [14]. Given a set of weights {w 1 , w 2 , · · · , w n }, the NURBS basis functions of degree p can be defined as Importantly, the knots can also be represented with two vectors, one of the unique knots {ζ 1 , ζ 2 , · · · , ζ m } and another of the corresponding knot multiplicities {r 1 , r 2 , · · · , r m }. As is standard in the literature of computer-aided design, we consider open knot vectors in this work, meaning that r 1 = r m = p + 1. We further assume that r i ≤ p for i = 2, · · · , m − 1. Across any given knot ζ i , the B-spline basis functions have α i := p − r i continuous derivatives. The vector α := {α 1 , α 2 , · · · , α m−1 , α m } = {−1, α 2 , · · · , α m−1 , −1} is referred to as the regularity vector. A value of −1 for α i indicates discontinuity of the basis functions at ζ i . We introduce the function space R p α := span{R p i } n i=1 , where the notation R p α is used to indicate that α i = α for i = 2, · · · , m − 1, suggesting continuity C α for the spline function spaces. The construction of multivariate B-spline and NURBS basis functions follows a tensor-product manner. For l = 1, 2, · · · , d, given p l , n l , the knot vectors Ξ l = {ξ 1,l , ξ 2,l , · · · , ξ n l +p l +1,l }, and the weight vectors {w 1,l , w 2,l · · · , w n l ,l }, the univariate NURBS basis functions R p l i l ,l are well-defined. Consequently, the multivariate NURBS basis functions can be defined by exploiting the tensor product structure, , for i l = 1, 2, · · · , n l and l = 1, 2, . . . , d. The tensor product NURBS space is denoted as

Mixed Formulation
In this work, we always consider three-dimensional problems (i.e. d = 3). Two discrete function spacesV h andQ h can be defined on the parametric domainΩ = (0, 1) 3 aŝ Assuming the physical domain can be parametrized by a geometrical mapping ψ :Ω → Ω X , the discrete functions on the physical domain can be defined through the pull-back operations, This pair of elements can be viewed as a smooth generalization of the Taylor-Hood element. We may then define the trial solution spaces for the velocity and pressure as and the corresponding test function spaces are defined as The semi-discrete formulation can be stated as follows. for Here v h0 and p h0 are the L 2 projections of the initial data onto the finite dimensional trial solution spaces.

The Generalized-α Schemes Under Comparison
In this section, we present the two generalized-α schemes to be compared as time integration schemes for the semidiscrete formulation in (6) and (7). Let the time interval [0, T ] be divided into a set of N ts subintervals delimited by a discrete time vector {t n } Nts n=0 . The time step size is defined as ∆t n := t n+1 − t n . The discrete approximations to the velocity, its time derivative, and pressure at the time step t n are denoted as v n ,v n , and p n . Let N A and M B represent basis functions for the velocity and pressure spaces, respectively, and let e i be the Cartesian basis vector with i = 1, 2, 3. We may define the residual vectors as follows, With these residual definitions, we may introduce the two generalized-α schemes as follows.
Scheme-1 At time step t n , givenv n , v n , p n and the time step size ∆t n , findv n+1 , v n+1 , and p n+1 such that, Although the quantityṗ n+1 does not appear explicitly in the above formulation, we may still obtain a consistent definition asṗ n+1 := (p n+1 − p n ) /∆t n .
Scheme-2 At time step t n , givenv n , v n , p n , and the time step size ∆t n , findv n+1 , v n+1 , and p n+1 such that, In this scheme, givenṗ n ,ṗ n+1 is given bẏ The definitions of the above two schemes can be completed by assigning values to the parameters α m , α f , and γ. It has been demonstrated that, for linear problems, second-order accuracy, unconditionally stability, and optimal high frequency dissipation can be achieved via the parametrization is the spectral radius of the amplification matrix at the highest mode [16]. In this study, ∞ is fixed to be 0.5.

Remark 1.
In conventional CFD calculations, it is unnecessary to calculateṗ. In hemodynamic simulations, however, bothv andṗ can be used as inputs to reduced-order models of the upstream and downstream vasculature, which are commonly used as boundary conditions coupled to the three-dimensional domain [20]. We therefore examine the accuracy ofv andṗ in addition to v and p.

Numerical Examples
The parameters and results reported below are presented in the centimeter-gram-second units.
Consistent with the original design of the benchmark, the domain is a cube Ω = (−1, 1) 3 . The constants are fixed as a = π/4 and d = π/2, and the fluid density ρ and viscosity µ are fixed to be 1.0 and 0.1. Traction boundary conditions are enforced on all faces, and initial conditions are obtained from the exact solution. Prior to examination of the temporal accuracy, a spatial convergence study was first performed, and the code was verified to yield optimal convergence rates. In the following, we report the errors for a fixed mesh with 25 elements along each dimension, polynomial degree p = 4, and continuity α = 3. The problem is simulated up to T = 1.0 with various uniform time step sizes. The relative errors of the velocity, pressure, and their time derivatives in the L 2 -and H 1 -norms at time t = 1.0 for Scheme-1 and Scheme-2 are presented in Tables 1 and 2, respectively. Based on the numerical results, we make the following salient observations. First, the pressure converges linearly with respect to the time step size in Scheme-1 and quadratically in Scheme-2. Second, both schemes exhibit first-order temporal accuracy for the time derivative of pressureṗ h , yet Scheme-2 yields significantly smaller errors forṗ h . Third, the velocity errors are identical in the two schemes, signifying that the degraded accuracy for pressure does not affect the accuracy for velocity. This phenomenon may, however, be attributed to the use of inf-sup stable elements. We do not rule out the possibility for velocity errors to be polluted by the use of Scheme-1 in stabilized formulations.  Table 1: Temporal convergence rates for the Ethier-Steinman benchmark using Scheme-1.

Pulsatile Flow in a Rigid Pipe
In the second example, we utilize the Womersley solution to evaluate the two schemes for flows under more physiological settings. The Womersley solution describes axisymmetric, fully developed flow subject to a periodically oscillating pressure gradient in a rigid cylindrical pipe with circular cross section [24,25]. Here, we consider the z-axis to be oriented along the length of the pipe and express the pressure in the following Fourier series, where p ref is the reference pressure at z = 0, k 0 is the steady component of the pressure gradient, k n is the amplitude of the oscillatory component of the pressure gradient at the n-th Fourier mode, ι is the solution of the equation ι 2 = −1, ω := 2π/T p is the fundamental frequency of the pressure wave, and T p is the period of the pressure. Per the assumptions of axisymmetric and fully developed flow, the velocity is identically zero in the radial and circumferential directions and takes the following analytical form in the axial direction, wherein r := x 2 + y 2 , R is the radius of the pipe, J 0 is the zeroth-order Bessel function of the first kind, and α n := R nω/ν is the Womersley number for the n-th Fourier mode. The complex form of the solutions indicates the existence of two sets of real independent solutions. We use the real parts of (8) and (9)     benchmark solution. In this study, the pipe radius R is 0.3, and the pipe length is 1.0. The fluid density ρ and viscosity µ are fixed to be 1.0 and 0.04, respectively, and the period T p is set as 1.1. The fundamental frequency ω and the corresponding Womersley number α 1 are approximately 5.71 and 3.59, respectively. The reference pressure p ref is 0. We choose to use a single mode in the Fourier series to represent the oscillatory part of the pressure (i.e. N = 1), and the values of k 0 and k 1 are −21.0469 and −33.0102 + 42.9332ι. The parameters above are chosen to represent physiological flows in arteries. The circular cross section is exactly represented by NURBS, and the parametrization is based on quadratic basis functions, which can be found in Example 6.2 of [21]. Using k-refinement, the spatial mesh for analysis consists of 10 × 10 elements in the x-y plane and 16 elements along the z axis, with polynomial degree p = 4 and continuity α = 3. We apply the no-slip boundary condition on the wall and traction boundary conditions on both ends of the pipe. The analytical forms of the tractions as well as the initial conditions are obtained from the exact solution in (8) and (9). Simulations are performed with uniform time step sizes ∆t n = ∆t = 1.00 × 10 −1 , 5.00 × 10 −2 , 2.50 × 10 −2 , 1.25 × 10 −2 , 5.00 × 10 −3 , 2.50 × 10 −3 , and 1.25 × 10 −3 . Similar to the previous example, we observed that the velocity errors in all norms are identical for the two schemes.
The velocity profiles at various time instances are illustrated in Figure 1. The relative errors of p h andṗ h in the L 2 -and H 1 -norms at time t = 0.8 are plotted in Figure 2. We observe that the pressure converges linearly in Scheme-1 and quadratically in Scheme-2 in both norms. Importantly, when ∆t = 1.25 × 10 −3 , the pressure errors in Scheme-2 are more than two orders of magnitude smaller than those in Scheme-1. Finally, while both schemes exhibit first-order temporal accuracy forṗ h , we again find thatṗ h calculated from Scheme-2 is more accurate.

Concluding Remarks
In this work, we implemented and compared two generalized-α schemes for the incompressible Navier-Stokes equations using inf-sup stable higher-order NURBS elements for spatial discretization. In the widely adopted scheme that is used in numerous existing CFD and FSI implementations, we have observed only first-order accuracy for pressure. While we have not observed degraded accuracy for velocity in these schemes, we cannot rule out this possibility for other spatial discretization methods, such as the stabilized method. By evaluating the pressure at the intermediate time step in the temporal scheme, second-order temporal accuracy for pressure is recovered. Moreover, the time derivative of the pressure is computed much more accurately. Considering pressure is one of the most important quantities for CFD and FSI, we recommend the use of Scheme-2 as the generalized-α scheme of choice for integrating the incompressible Navier-Stokes equations.