{'fraction_non_alphanumeric': 0.0750786038682631, 'fraction_numerical': 0.034744497729221585, 'mean_word_length': 3.4625850340136055, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A discontinuous Galerkin pressure correction scheme for the incompressible Navier-Stokes equations: stability and convergence

A discontinuous Galerkin pressure correction numerical method for solving the incompressible Navier-Stokes equations is formulated and analyzed. We prove unconditional stability of the propose scheme. Convergence of the discrete velocity is established by deriving a priori error estimates. Numerical results verify the convergence rates.


Introduction
The numerical simulation of the incompressible Navier-Stokes equations presents a challenging computational task primarily because of two reasons: (a) the coupling of the velocity and pressure by the incompressibility constraint and (b) the nonlinearity of the convection term [14,18]. The development of splitting schemes aims to overcome these difficulties by decoupling the nonlinearity in the convection term from the pressure term. For an overview of such methods, we refer to the works of Glowinski [15] and of Guermond, Minev, and Shen [18]. In this paper, we will focus on pressure correction schemes. The basic idea of a non-incremental pressure correction scheme in time was first proposed by Chorin and Temam [5,28]. This scheme was subsequently modified by several mathematicians leading to two major variations: (1) the incremental scheme where a previous value of the pressure gradient is added [16,30] and (2) the rotational scheme where the non-physical boundary condition for the pressure is corrected by using the rotational form of the Laplacian [29].
The main contribution of our work is the theoretical analysis of a discontinuous Galerkin (dG) discretization of the pressure correction approach. We derive stability and a priori error bounds on a family of regular meshes. The discrete velocities are approximated by discontinuous piecewise polynomials of degree k 1 and the discrete potential and pressure by polynomials of degree k 2 . Stability of the solutions is obtained under the constraint k 1 −1 ≤ k 2 ≤ k 1 +1 whereas the convergence of the scheme is obtained for the case k 2 = k 1 − 1 because of approximation properties. The proofs are technical and rely on several tools including special lift operators.
The semi-discrete error analysis of pressure correction schemes has been extensively studied, see for example the work by Shen and Guermond [21,27]. The use of finite element approximations in conjunction with pressure correction schemes is also well studied and a priori error estimates are established. Without being exhaustive, we refer to the work by Guermond and Quartapelle [19,20], and to the work by Nochetto and Pyo [24] for the analysis of finite element methods.
More recently, the combination of dG spatial approximations with pressure correction formulation in time to solve the incompressible Navier-Stokes equations has been the subject of several computational papers. These dG methods have multiple important features including local mass conservation and high order convergence.
In what follows, we mention a non-exhaustive list of such publications. Botti and Di Pietro employ a dG approximation to the velocity and a continuous Galerkin approximation to the pressure [2]. Liu et al. formulate an interior penalty dG method with the pressure correction approach for the computational estimation of rock permeability [23]. Piatkowski et al. use a modified upwind scheme based on Vijayasundaram numerical flux, postprocess the projected velocity so that it is discretely divergence free, and perform numerical experiments of large scale 3D problems [25]. We also mention the work by Fehn et al. where the stability of pressure correction and velocity correction dG methods is numerically investigated for small time step sizes, and the robustness of the methods is demonstrated for laminar high Reynolds number flow problems [9].
Error analysis of dG methods for the steady state Navier-Stokes equations has been established, see the work by Girault, Riviere and Wheeler [13]. The authors also prove a priori error estimates for the time dependent Navier-Stokes equations where an operator splitting scheme is combined with dG methods in [14]. To the best of our knowledge, the a priori error analysis for a pressure correction dG approximation of the incompressible Navier-Stokes equations is missing from the literature.
The outline of this paper is as follows. In Section 2, the model problem and the splitting scheme in time are presented. Section 3 introduces the dG forms and summarizes their properties. We present the fully discrete scheme and show existence and uniqueness of solutions in Section 4. Unconditional stability is established in Section 5, and convergence of the discrete velocity is shown in Section 6. Conclusions follow.

Model problem and time discretization
Consider the incompressible time dependent Navier-Stokes equations with homogenous Dirichlet boundary condition for the velocity and a zero average constraint for the pressure. viscosity. Let τ denote the time step size and consider a uniform partition of the time interval (0, T ] into N T subintervals. Throughout the paper, we use the notation g n = g(t n ) and q n = q(t n ) for given functions g and q evaluated at t n = nτ . We now formulate a variation of the consistent splitting scheme that splits the operators and introduces additional velocity v n and potential φ n [18,23]. Initially, set p 0 = 0. For n = 1, . . . , N T , given u n−1 and p n−1 , compute an intermediate velocity v n such that v n − τ µ∆v n + τ u n−1 · ∇v n + τ ∇p n−1 = u n−1 + τ f n , in Ω, (2.6) v n = 0, on ∂Ω.
In (2.11), δ is a positive parameter that can be chosen in the interval ]0, 1/(4d)], as stated in the stability and convergence results proved in this work. The vector n is the unit outward normal vector to ∂Ω.

Notation and discontinuous Galerkin forms
For a given non-negative integer m and real number r ≥ 1, define the Sobolev space on a domain O ⊂ Ê d : where α is a multi-index and D α w is the corresponding weak partial derivative of w. The usual Sobolev norms and semi-norms are denoted by · W m,r (O) and |·| W m,r (O) respectively. If r = 2, we denote . The L 2 inner-product over Ω is denoted by (·, ·) and the resulting L 2 norm by · . Let E h = {E k } denote a family of regular partitions of the domain Ω [6]. That is, there exists a constant ρ independent of h such that where h E = diam(E) and ρ E is the maximum diameter of a ball inscribed in E.
Let h denote the maximum diameter of the elements in E h . We define the following broken Sobolev spaces.
Denote by Γ h the set of all interior faces of the subdivision E h . For an interior face e ∈ Γ h , we associate a normal n e and we denote by E 1 e and E 2 e the two elements that share e, such that n e points from E 1 e to E 2 e . Define the average and jump for a function θ ∈ X as such, For a boundary face, e ∈ ∂Ω, the vector n e is chosen as the unit outward vector to ∂Ω. The definition of the average and jump in this case are extended as such, Similar definitions are used for scalar valued functions, q ∈ M . We now introduce the dG formulations for the spatial discretization of (2.6)-(2.12). For the convection term in (2.6), we use the same discretization form, a C , as in [13]. To define this form, we use the following notation: the vector n E denotes the outward normal to E, the trace of a function v on the boundary of E coming from the interior (resp. exterior) of E is denoted by v int (resp. v ext ). By convention, v ext | e = 0 if e is a boundary face (e ⊂ ∂Ω). We also introduce the notation for the inflow boundary of E with respect to a function z ∈ X: {z(x)} · n E < 0}. With this notation, we define for z, w, v, θ ∈ X, We recall the positivity property satisfied by a C (see (1.18) in [13]): In the analysis, it will be useful to separate the upwind terms from the form a C (w; w, v, θ). To this end, for z, w, v, θ ∈ X, we define: We also define the upwind terms It follows that: The discretization for the elliptic operator, −∆v, is given as follows [26]. For v, θ ∈ X, In the above form, h e = |e| 1/(d−1) , ǫ ∈ {−1, 0, 1}, σ > 0 is a user specified penalty parameter. The discretization for the term −∇p is given as follows. For θ ∈ X and q ∈ M , define To approximate u and p, we introduce discrete function spaces X h ⊂ X and In the above, for k ∈ AE, È k (E) denotes the space of polynomials of degree at most k. We will assume that (for reasons that will be evident in the analysis) To discretize the elliptic operator −∆φ, we define for φ h , q h ∈ M h , Here,σ > 0 is a penalty parameter. For θ ∈ X, define the energy norm as follows: For q ∈ M , the energy semi-norm is defined as such: Clearly, | · | DG is a norm for the space M h0 . We recall the following coercivity properties.
It is shown that (3.21) holds with κ = 1 if ǫ = 1 and with κ = 1/2 if ǫ ∈ {−1, 0} and σ is large enough [26]. Property (3.22) holds forσ large enough [26]. In what follows, we will assume that (3.22) and (3.21) are satisfied. Further, we recall that a ǫ is continuous on X h × X h with respect to · DG and a ellip is continuous on M h × M h with respect to | · | DG [26]. There exists a constant C independent of h, θ h , z h , q h , and ξ h such that We now state an equivalent expression for the form b.
Lemma 3.1. We have the following equivalent form for b(θ, q), Proof. The result is obtained by applying Green's theorem to the first term in (3.14) and by using the identity We will also make use of lift operators, which are useful tools in the theoretical analysis of dG methods. Given e ∈ Γ h ∪ ∂Ω, we introduce a new lift operator r e : (L 2 (e)) d → M h as follows We next recall the lift operator g e introduced in [7]: given an interior face e ∈ Γ h , the operator g e : L 2 (e) → X h satisfies With these definitions, we construct two operators, R h : It is easy to check that R h and G h are linear operators. We next show boundedness of the lift operators.
Lemma 3.2. There exist constants M k2 ,M k1 > 0 independent of h but depending on the polynomial degrees k 2 and k 1 respectively, such that the following bounds hold: Proof. We will show (3.31). We use the definitions of R h (3.29) and of r e (3.27). We have Cauchy-Schwarz's inequality yields With a local trace inequality and the fact that h E ≥ h e for any element E having a face e, we obtain the result with M k2 depending on a trace constant. The proof of (3.32) follows a similar argument and is omitted for brevity.
Further, with the definitions of the lift operators (3.29) and (3.30), we have the following equivalent forms to (3.14) and (3.25) where ∇ h and ∇ h · denote the broken gradient and divergence operators respectively.

Numerical scheme
We start with setting p 0 h = φ 0 h = 0. We let u 0 h be the local L 2 projection of u 0 onto X h .
Next, compute φ n h ∈ M h0 such that for all q h ∈ M h0 , . Finally, compute p n h ∈ M h and u n h ∈ X h such that for all q h ∈ M h and for all θ h ∈ X h , Using (3.33) and (3.34), steps (4.4) and (4.5) can be written as follows. For all q h ∈ M h and θ h ∈ X h , Lemma 4.1. Let n ≥ 0 and p n h ∈ M h be defined by (4.6), then p n h ∈ M h0 . Proof. We present a proof by induction on n. For n = 0, the statement trivially holds. Assume p n−1 h ∈ M h0 . Let q h = 1 in (4.4) and use (3.34). We obtain: We conclude the result by using the induction hypothesis and the fact that φ n h ∈ M h0 . , χ n h , χ n h ) + τ µa ǫ (χ n h , χ n h ) = 0. With the positivity property of a C (3.9) and the coercivity of a ǫ (3.21), we conclude that (4.9) χ n h 2 + κτ µ χ n h 2 DG ≤ 0. This implies that χ n h = 0. Thus, the solution to (4.2) is unique. Existence of u n h and p n h : The existence of φ n h ∈ M h0 follows by similar arguments, the coercivity property of a ellip (·, ·), and the fact that | · | DG is a norm for the space M h0 . The existence of p n h ∈ M h0 follows directly from (4.6) and Lemma 4.1. The existence of u n h ∈ X h follows from (4.7).

Stability
In this section, we present a stability analysis of the proposed scheme. We proceed by noting the following identity.
. For all q h ∈ M h and n ≥ 1, the following holds. ).
In addition, we have Proof. Let q h ∈ M h be arbitrary but fixed. Take θ h = ∇ h q h in (4.7). We have, By applying (4.3), we obtain (5.2).
We recall the following L q bound for the broken Sobolev spaces. Since X ⊂ H 1 (E h ) d , the proof of this lemma can be found in Lemma 6.2 in [13] in 2D and in Corollary 2.2 in 3D in [22]. Lemma 5.2. There exists a constant C P independent of h and τ but depending on q such that To state and prove our stability result, for n ≥ 0, we introduce auxiliary functions Note that S n h ∈ M h0 since with the equivalent definition (3.33), we have (5.10) This implies that ξ n h ∈ M h0 since p n h ∈ M h0 , see Lemma 4.1. We are now ready to prove the main result of this section.
, and u 0 ∈ L 2 (Ω) d . Then, the pressure correction scheme (4.2)-(4.7) is unconditionally stable. For all τ > 0 and 1 ≤ m ≤ N T , . We use the positivity property of a C (3.9) and the coercivity of a ǫ (3.21). We obtain the following.
Similarly, using (5.2) in Lemma 5.1, we obtain Substituting the above expressions in (5.13), we obtain (5.14) The last term is estimated using (3.32): Using the above estimate in (5.14) yields With the assumption on the penalty parameter,σ ≥M 2 k1 , we obtain We add (5.17) to (5.12), To handle the first term on the right-hand side of the above inequality, we use the auxiliary function S n h defined in (5.8) and write With the equivalent expression (3.33), the second term in (5.19) is rewritten: Recall the definition of ξ n h in (5.9) and note that (4.6) implies Since a ellip (·, ·) is symmetric, we have With the above expression, (5.18) becomes To handle the second term on the right-hand side of (5.25), recall that from the definition of S n h , see (5.8), Using (3.31), the assumption that δ ≤ κ/(2d) and Substituting the above bound in (5.25) yields, Using Cauchy-Schwarz's inequality, Poincare's inequality (5.7), and Young's inequality, we have We substitute the above abound in (5.27), multiply by 2 and sum the resulting inequality from n = 1 to n = m. We obtain Recall that ξ 0 h = S 0 h = 0. With the coercivity property of a ellip (3.22) and the stability of the L 2 projection, we attain the result.
Remark 5.4. The particular choice for the upper bound on δ in the statement of Theorem 5.3 is chosen to simplify the writeup. However, it can be altered and a similar result would hold. For example, the condition on δ can be chosen as δ < κ/d. This will lead to a larger constant C resulting from Young's inequality in (5.28).

Error estimates for the velocity
In this section, we denote k = k 1 and consider the case k 2 = k − 1. To simplify notation, for any function v ∈ (W 1,3 (Ω) ∩ L ∞ (Ω)) d , we define the norm: We assume that u 0 ∈ H 1 0 (Ω) d and that ∇ · u 0 = 0. Throughout the rest of this paper, we denote by C a constant independent of µ, h, and τ and by C µ a constant independent of h and τ but depending on e 1/µ . These constants may take different values when used in different places and may depend on certain norms of the true solution (u, p). These norms are associated with the spaces specified by the regularity assumptions of Theorem 6.7. 6.1. Approximation. We will make use of the operator Π h : H 1 0 (Ω) d → X h that preserves the discrete divergence [4,26]: This operator satisfies the following approximation properties [4].
In the subsequent sections, the following stability bound for the operator Π h will be used.
Then, we have the following bound.
Minkowski's inequality, an inverse estimate, stability and approximation properties of the Lagrange interpolant and (6.3), we have ) and the above bound is uniform in E, we obtain the result.
Define the local L 2 projection π h : L 2 (Ω) → M h as follows. For 0 ≤ t ≤ T , a given function p(t) ∈ L 2 (Ω), and any E ∈ E h , For 0 ≤ t ≤ T and p(t) ∈ H s (Ω), the following estimate holds. For all E ∈ E h , We also recall that the local L 2 projection is stable in the dG norm [10].

Error equations.
We recall that the dG discretization is consistent, in the sense that the true solution (u, p) of (2.1)-(2.5) satisfies for n ≥ 1 For readability, we define the following discretization errors,ẽ n h ∈ X h , e n h ∈ X h : In the above, we set v 0 h = u 0 h . Multiplying (6.11) by τ and subtracting it from (4.2), we have for all n ≥ 1 In the above, we set With (6.2) and (3.33), we have for all n ≥ 1 Using (6.16) and (4.3), we obtain for all n ≥ 1 Inserting Π h u n in (4.7) yields for all n ≥ 1 In addition, for n ≥ 1, where δ n,1 is the Kronecker delta and A n 1 and A n 2 for n ≥ 1 are given by Proof. Using (6.16) in (5.1), we obtain (6.19) for n ≥ 1. Substituting (6.17) in (6.19), we obtain (6.20) for n ≥ 1. To show (6.21), we use (6.18): For the second term, we use the equivalent definition (3.34).
We proceed by presenting some bounds for the forms a C , U, and C. The proof of these bounds is inspired from the proof of proposition 4.1 in [14]. However, the proof differs since in our scheme, we do not have that b(e n h , q h ) = 0 for all q h ∈ M h .
Then, there exists a constant C, independent of h, τ , w h , z, v and θ h , such that the following estimates hold.
(i) If v ∈ (W 1,3 (Ω) ∩ H 1 0 (Ω) ∩ L ∞ (Ω)) d , then for any z ∈ X and θ h ∈ X h : Proof. (i) Since v belongs to H 1 (Ω) d and vanishes on the boundary, we have: With Hölder's inequality and (5.7), we have To bound the second term in (6.30), we use a similar argument as in [14]. Define c 1 , c 2 as piecewise constant vectors where on each element E ∈ E h , (6.32) We then write, The first term is bounded in Proposition 4.1 in [14] in 2D domains. We follow a similar technique to obtain a bound in 2D/3D domains. We have, For the first term, note that by Hölder's, inverse and triangle inequalities, we have: ). To proceed, we recall Poincare's inequality, see Theorems 1 and 2 in Section 5.8 in [8]. For r ∈ [1, ∞], E ∈ E h and θ in W 1,r (E), Applying (6.35), we find To handle the second term in the right-hand side of (6.34), consider a face e ∈ Γ h and let E 1 e and E 2 e denote the elements sharing e. Terms involving faces e ∈ ∂Ω are handled similarly. With Cauchy-Schwarz's inequality, we have Let us consider the case i = 1 and j = 2; the other terms are handled similarly: . With a trace inequality and Hölder's inequality, we have . (6.37) Therefore, with discrete Hölder's inequality and (5.7), we obtain We bound the other terms similarly, and combining the bounds, we obtain It remains to handle the second term in the right-hand side of (6.33). To this end, we use (6.27).
Using Cauchy-Schwarz's inequality and (3.32), we have (6.40) To bound T 1 , we apply a trace inequality and (6.35), as it was done in (6.37). Consider the elements E 1 e and E 2 e that share e. We have, for i = 1, 2: . Hence, we obtain the following bound on T 1 .
To handle T 2 , we split it as such.
Thus, we have To bound T 1 2 , we will use a trace inequality, Hölder's inequalities, (6.35), and (5.7). We skip the details as the argument is similar to the one used in (6.37).
(ii) First, we remark that by a Sobolev embedding, v ∈ (W 1,3 (Ω) ∩ L ∞ (Ω)) d since k ≥ 1. We have, The first term is bounded exactly like in the proof of (i) (see the derivation of bound (6.31)) with v replaced by v − Π h v. By the stability of the interpolant, we have The second term is split in the same way as in the proof of (i).
where c 2 is defined in (6.32) andc 1 is a piecewise constant vector defined on each element E ∈ E h asc Similar to (6.38) and using the stability of the interpolant (6.7), we have the bound With similar arguments as equations (6.40), (6.41) and (6.43), we have HereT 1 ,T 1 2 andT 2 2 are defined in a similar way to T 1 , T 1 2 and T 2 2 with c 1 replaced byc 1 and v replaced by v − Π h v. The termsT 1 andT 1 2 are bounded exactly like T 1 (6.42) and T 1 2 (6.44) with v replaced by v − Π h v. With the stability of the interpolant, we have The only term that differs isT 2 2 . We have, Consider an edge e and an adjacent element E 1 e . Using trace and Hölder's inequalities, we have , With the approximation property (6.3), we have: Hence, with discrete Hölder's inequality, (5.7), (6.4), and (6.7),T 2 2 is bounded as follows.
In addition, we have the following bounds on a C and C presented in the next lemma.
Lemma 6.5. There exists a constant C, independent of h, τ, w, z, v, w h , and θ h such that the following bounds hold.
) d , then for any z ∈ X and θ h ∈ X h : Proof. (i) Note that since v has zero jumps and vanishes on the boundary, the upwind term vanishes. We have The first term is bounded similarly as before (with (5.7) and (6.5): Recall the definition of c 1 , c 2 in (6.32). Since b(Π h w − w, q h ) = 0 for any q h ∈ M h , we have Using Hölder's inequality and (6.35), we have the following bounds.
Hence, with the help of the discrete Hölder's inequality and (5.7), the volume term in the right-hand side of (6.62) is bounded as follows.
The face terms in the right-hand side of (6.62) are handled in a slightly different way as in (6.37), with w h replaced by Π h w − w. For instance, one of the terms is bounded as follows: ). With (6.3), (6.4), we have ). Hence, with discrete Hölder inequality and (5.7), the face terms are bounded as follows.
(ii) To show (6.58), we write With Hölder's inequality, stability of the interpolant and (6.5), the first term is bounded by: To handle the second term, we will use (3.25). We have: We bound B 1 with Hölder's inequality, stability of the interpolant (6.7), and the approximation properties (6.5)-(6.6). We have For B 2 , let e be an interior face in Γ h shared by E 1 e and E 2 e . We have, With Hölder's inequality and a trace estimate, we have for i = 1, 2: . We apply the approximation properties, sum over the faces, use (5.7), and obtain Bounds (6.67), (6.69), and (6.70) yield (6.58). To handle the upwind term, we consider the faces' contributions to the upwind terms. For details on such contributions, see Propositon 4.10 in [12]. We bound it by: Consider the case e ∈ Γ h . With a trace inequality and approximation properties (6.3)-(6.4), we have . With similar bounds for the case e ∈ ∂Ω, discrete Hölder's inequality, and (6.7), we conclude that (6.59) holds.
(iii) This bound can be found in Proposition 4.1 in [14] for k = 1 and d = 2. Here, we extend it to any polynomial degree k and for d = 2, 3. We will use local trace and inverse estimates. For i = 1, 2, we have: We also consider the faces' contributions to handle the upwind term. With Hölder's inequality, we bound the upwind term, U(z; w h , v − Π h v, θ h ) by: Consider the case when e ∈ Γ h and j = k. The other cases are handled similarly. We use (6.71), (6.72), and (6.3)-(6.4): , we obtain This bound along with similar bounds to (6.73), discrete Hölder's inequality, and (5.7) show estimate (6.60).
There exist positive constants γ and C independent of h, τ and µ such that if τ ≤ γµ, the following bound holds. For n ≥ 1, Proof. We follow a similar technique to [14]. Since u n belongs to H 1 0 (Ω) d , we can write a C (u n ; u n , u n ,ẽ n h ) = a C (u n−1 h ; u n , u n ,ẽ n h ). Therefore, we have ; Π h u n−1 − u n , u n ,ẽ n h ) = W n 1 + W n 2 + W n 3 . With the help of (3.12), we write W n 1 as follows. ; Π h u n−1 , Π h u n − u n ,ẽ n h ). For n = 1, with (6.48) and (6.49), we have: Since (6.20) holds, we can apply (6.29) in Lemma 6.4 to bound the first term in W n 1 for n ≥ 2. We also use estimates (6.58) for the second term, (6.60) for the third term, and (6.59) along with (5.7) for the fourth term. For n ≥ 1, we obtain Further, for k ≥ 1, with Morrey's inequality and a Sobolev embedding (see Theorem 1.4.6 in [3] and Theorem 2 in section 5.6.1 [8]), we have for n ≥ 0: With the help of (6.75), the assumption that u ∈ L ∞ (0, T ; H k+1 (Ω) d ), and Young's inequality, we obtain for constants C and C 1 : With (6.30) and (6.31), we have Since (6.20) holds, we can apply (6.28) to bound W n 2 for n ≥ 2. Thus, we have With (6.75) and Young's inequality, we obtain for constants C and C 2 To bound W n 3 , we split it into W n 3 = a C (u n−1 h ; Π h u n−1 − u n−1 , u n ,ẽ n h ) + a C (u n−1 h ; u n−1 − u n , u n ,ẽ n h ). Since b(u n−1 , q h ) = 0, ∀q h ∈ M h , ∀n ≥ 1, we apply (6.57) and the assumption that u ∈ L ∞ (0, T ; H k+1 (Ω) d ) to bound the first term.
Since u n−1 and u n belong to H 1 0 (Ω) d and their divergence is zero, we have By Hölder's inequality, (5.7), and (6.75), we obtain The result follows by combining all the above bounds, and by assuming (C 1 +C 2 )τ ≤ κµ/96.
We are now ready to present and prove the main convergence result.
Theorem 6.7. Assume that σ ≥ M 2 k−1 /d,σ ≥ 4M 2 k , and δ ≤ κ/(2d). There exists a positive constant γ 1 , independent of h, τ and µ, such that if τ ≤ γ 1 µ, the following estimate holds. For 0 ≤ m ≤ N T , The above estimate holds under the following regularity assumptions: Proof. Let θ h =ẽ n h in (6.13). Using the positivity property of a C (3.9) and the coercivity property of a ǫ (3.21), we have . Substituting (6.21) in the above expression, we obtain We begin by handling the fifth and sixth terms. With (3.32) and the assumption thatσ ≥ 2M 2 k , we have: Using Cauchy-Schwarz's inequality, Young's inequality, (3.32), and the assumption thatσ ≥ 4M 2 k , we obtain Next, let θ h = e n h in (6.18) and use (3.34 Similar to the stability proof, we let θ h = e n h −ẽ n h in (6.18) and use (6.19): Hence, (6.81) reads We add (6.83) to (6.77). With bounds (6.78) and (6.79), we obtain (6.84) . We now focus on the second term on the right hand side of (6.84). We write Since the operator π h preserves cell averages and p n satisfies (2.5), we know that π h p n ∈ M h0 . By (6.17), continuity of a ellip (·, ·) on M h × M h see (3.24), Young's inequality, and (6.10), we obtain Since ∇ h ·ẽ n h ∈ M h , by the definition of π h p n we have: By using trace inequalities and (6.9), we obtain Note that by (6.2), b(ẽ n h , p n−1 ). Hence, from (5.24), we have With this expression, (6.87) becomes (6.89) We now consider R t (ẽ n h ) (see (6.15)). Using Cauchy-Schwarz's inequality, Young's inequality and (5.7), we obtain With a Taylor expansion and approximation property (6.5), we have With the definition of S n h (5.8) and (6.16), we have With the assumptions that δ ≤ κ/(2d) and σ ≥ M To handle a ǫ (Π h u n − u n ,ẽ n h ), we write: The terms Q 1 , Q 3 and Q 4 are handled via standard arguments. By using trace inequalities, we have The term Q 2 is more delicate. We follow an approach similar to the one in Chapter 6 of [26] and [13]. By Cauchy-Schwarz's inequality, we have Consider one element E 1 e adjacent to e. By the trace theorem, we have ). Letũ n ∈ X h be an interpolant of u n , see Theorem 2.6 in [26], satisfying: By the triangle inequality, an inverse estimate, the above bound and (6.6), we obtain: Hence, (6.95) becomes The term involving the neighboring element E 2 e is handled similarly. Thus, we can conclude that (6.97) We employ the above bounds on Q 1 -Q 4 , the bound on the nonlinear terms |R C (ẽ n h )| (6.74), (6.92), and the bound on |R t (ẽ n h )| (6.90) in (6.89). With these bounds and the assumptions that u ∈ L ∞ (0, T ; H k+1 (Ω) d ) and p ∈ L ∞ (0, T ; H k (Ω)), (6.89) becomes: We sum (6.98) from n = 1 to n = m, multiply by 2, and use the regularity assumptions. We obtain We recall that φ 0 h = 0 and the definitions of A n 1 and A n 2 (see (6.22) -(6.23)). We use (3.32) with the assumption thatσ ≥M 2 k to obtain m n=1 With (3.34), (3.32), and approximation properties we have With Hölder's inequality, trace inequality, stability of the interpolant, and (5.7), we obtain The above bounds, approximation properties, the coercivity property of a ellip (3.22), discrete Gronwall's inequality, and the triangle inequality yield bound (6.76).
Two immediate results of the above arguments are as follows. Under the same regularity assumptions of Theorem 6.7 and for τ ≤ γ 1 µ and for 1 ≤ m ≤ N T , ≤ C µ 1 + µ + 1 µ (τ + h 2k ).
Remark 6.8. Assume that the mesh is quasi-uniform. If we use (4.7) and the following inverse inequality.

Numerical experiments
In this section, we use the manufactured solution method to compute the convergence rates. Let Ω = (0, 1) 3 denote the computational domain and let T = 1 be the end time for all numerical experiments. The domain Ω is partitioned into cubic elements. We select the following Beltrami flow with parameter µ = 1 as the prescribed solution. u x (x, y, z, t) = − exp (−t + x) sin (y + z) − exp (−t + z) cos (x + y), u y (x, y, z, t) = − exp (−t + y) sin (x + z) − exp (−t + x) cos (y + z), u z (x, y, z, t) = − exp (−t + z) sin (x + y) − exp (−t + y) cos (x + z), p(x, y, z, t) = − exp (−2t) exp(x + z) sin (y + z) cos (x + y) + exp(x + y) sin (x + z) cos (y + z) + exp(y + z) sin (x + y) cos (x + z) Herep = 7.63958172715414 is used to force the average pressure over Ω equal to zero (up to machine precision). The Beltrami flow has the property that the velocity and vorticity vectors are parallel to each other and that the nonlinear convection term is balanced by the pressure gradient, namely, ∂ t u − ∆u = 0 and u · ∇u + ∇p = 0, which implies f = 0. The initial condition and Dirichlet boundary condition on ∂Ω for velocity are imposed by the exact solution. First, we derive temporal rates of convergence by computing the solutions with time step size τ ∈ {1/2 3 , 1/2 4 , ..., 1/2 7 }. The È2 − È1 dG scheme with mesh resolution h e = 1/160 is employed to guarantee the error in space does not dominate the error in time. We set ǫ = −1, σ = 64 on Γ h and σ = 128 on ∂Ω. The other penalty parameterσ is set equal to 2. If err τ denotes the error when time step size is equal to τ , the convergence rate is defined by ln(err τ /err τ /2 )/ ln (2). Table 1 displays the errors and convergence rates for the velocities and pressure. We observe second order convergence rate in time, which is higher than the expected first order rate. Next, we obtain spatial rates of convergence by computing the solutions on a sequence of uniformly refined meshes (see the second column of Table 2 for h e ). We fix τ = 1/2 10 for the È1 − È0 dG scheme, we fix τ = 1/2 13 for the È2 − È1 dG scheme, and we fix τ = 1/2 15 for the È3 − È2 dG scheme to guarantee the spatial error dominates. For È1−È0, we setσ = 1, ǫ = −1, σ = 8 on Γ h and σ = 16 on ∂Ω. For È2−È1, we setσ = 2, ǫ = −1, σ = 64 on Γ h and σ = 128 on ∂Ω. For È3−È2, we setσ = 8, ǫ = −1, σ = 128 on Γ h and σ = 256 on ∂Ω. If err he is the error on a mesh with resolution h e , then the rate is defined by ln(err he /err he/2 )/ ln (2). We show the errors and rates in  Table 1. Errors and temporal convergence rates of velocity and pressure.

Conclusion
The main contribution of this work is a theoretical error analysis of a dG pressure correction scheme for solving the incompressible Navier-Stokes equations. We establish stability of the method and derive a priori error estimates for the discrete velocity that are optimal in the broken gradient norm. The convergence analysis is technical and relies on lift operators and appropriate local approximations. Future work will consist of deriving error bounds for the pressure and optimal error estimates for the discrete velocity in the L 2 norm.  Table 2. Errors and spatial convergence rates of velocity and pressure.