Resolution and implementation of the nonstationary vorticity velocity pressure formulation of the Navier–Stokes equations

This paper deals with the iterative algorithm and the implementation of the spectral discretization of time-dependent Navier–Stokes equations in dimensions two and three. We present a variational formulation, which includes three independent unknowns: the vorticity, velocity, and pressure. In dimension two, we establish an optimal error estimate for the three unknowns. The discretization is deduced from the implicit Euler scheme in time and spectral methods in space. We present a matrix linear system and some numerical tests, which are in perfect concordance with the analysis.


Introduction
The nonlinear Navier-Stokes equations model the flow of a viscous and incompressible fluid such as water, air, and oil in stationary or nonstationary states. Those equations were and are the subject of a large number of research papers. The modification of any of the parameters associated with these equations (the domain on which the equations are posed, boundary conditions, nature of the data, variational formulation, time dependence, choice of the approximation method, etc.) leads to new research problems. In the initiating paper [1] the authors handle the Stokes and Navier-Stokes equations with nonstandard boundary conditions on the velocity and the pressure for a convex or regular domain. Our interest concerns the nonstationary Navier-Stokes equations with boundary conditions on the normal component of the velocity and the tangential components of the potential vector vorticity. Such a problem allows us to model, for instance, two fluids separated by a membrane or the flow in a network of pipes. The equivalent variational formulation of the Navier-Stokes equations provided with these boundary conditions admits three unknowns: the vorticity, velocity, and pressure [2][3][4][5]. This formulation has been studied in several works for the finite element discretization of the Stokes and Navier-Stokes problems in the stationary case (see [3,6]). We cite in the same context the works of Bernardi et al. [7,8], which present a posteriori error analysis of time-dependent Stokes and Navier-Stokes problems. The extension to spectral discretization has been handled in [9,10] for the stationary Stokes and Navier-Stokes problems and in [11,12] for the nonstationary case.
In this paper, we propose a discretization of such a formulation by the implicit Euler scheme in time and the spectral method in space in the square ] -1, 1[ 2 for dimension two and in the cube ] -1, 1[ 3 for dimension three. The spectral method can be easily extended to more complex geometries thanks to the arguments in [13,14]. In dimension two, we prove an optimal error estimate for the vorticity and the velocity and a quasioptimal error for the pressure, using the theorem of Brezzi, Rappaz, and Raviart [15]. However, the extension to dimension three remains a difficult problem.
We describe a numerical algorithm used to solve the discrete nonlinear problem. We also present clearly the matrices and the linear system derived from the discrete problem. This linear system is solved using the GMRES method since the global matrix is not symmetric [16].
Finally, we present some numerical experiments, which confirm a good convergence of our algorithm and the benefit of this formulation. These numerical results are coherent with the theoretical ones. This paper is organized as follows: • In Sect. 2, we present a continuous problem and some regularity results.
• Sect. 3 is devoted to a description of time and full discrete problems.
• An error estimate is derived in Sect. 4.
• In Sect. 5, we describe an iterative algorithm used to solve a nonlinear discrete problem and a linear matrix system. We conclude by presenting some numerical experiments.

A continuous problem and some regularity results
We consider an open bounded simply connected domain of R d (d = 2, d = 3) with Lipschitz and connected boundary . Let x = (x, y) for d = 2 or x = (x, y, z) for d = 3 be the Cartesian coordinates. In this paper, we mainly focus on the following nonstationary Navier-Stokes system: where v and P are the unknowns velocity and pressure, f represents the density of the body forces, ν > 0 is the viscosity, and n is the unit outward normal vector on the boundary . We define the boundary operator ς such that ς(curl v) is the boundary of curl v in dimension d = 2 or the boundary of the tangential components of the curl v in dimension We introduce the unknown vorticity θ = curl v (see [2,3]), and since v.∇v is equal to θ × v + 1 2 grad |v| 2 , system (1) is equivalent to the system The dynamical pressure p is defined as We assume that the following condition is satisfied by the initial velocity v 0 and the initial vorticity θ 0 = θ (x, 0): We define the space and its subspace H 0 (div, ) = u ∈ H(div, ); u.n = 0 on .
Let B a separable Banach space. We need to define the following space to handle the nonstationary Navier-Stokes system: which is a Banach space with the norm We also introduce the Banach space L(B) of the continuous linear functions from B to R with the norm If the data f belongs to the space L 2 (0, T; H 0 (div, ) ), where (H 0 (div, ) is the dual space of H 0 (div, ) (see [17] for more detail ), then problem (2) is equivalent to the following variational formulation: where ≺ ·, · is the duality product between H 0 (div, ) and H 0 (div, ). The bilinear forms l(·, · ; ·), b(·, ·) and t(·, · ; ·) are defined as follows: In another way, we define the nonlinear term Z(·, · ; ·) by For proving the existence of solution for problem (2), we need to define the following two kernels: the kernel of the bilinear form b(·, ·), which coincides with the space of divergence-free functions in H 0 (div, ), and W = (ϑ, ϕ) ∈ H 0 (curl, ) × K; ∀ψ ∈ H 0 (curl, ), t(ϑ, ϕ; ψ) = 0 the kernel of the bilinear form t(·, · ; ·). From the continuity of the bilinear forms b(·, ·) and t(·, · ; ·) we deduce that K and W are Hilbert spaces. If (θ , v, p) is a solution of problem (4), then (θ , v) is solution of the following reduced problem: Find (θ , v) ∈ L 2 (0, T; W) such that In dimension two, it is simple to show that problem (5) has a solution. However, in dimension three, giving a sense to the nonlinear term Z(·, · ; ·) is related to the following Assumption 1. In that case, the spaces H 0 (div, ) ∩ H(curl, ) and H(div, ) ∩ H 0 (curl, ) are compactly embedded in H 1 ( ) see ( [18], Thm 2.17).

Assumption 1
In dimension three, we suppose that the boundary is C 1,1 or convex.
We recall the uniform inf-sup condition on bilinear form b(·, ·): There exists a constant γ > 0 such that see [19] or [20, Chap. I, Cor. 2.4] for its proof. When Assumption 1 and the inf-sup condition (6) are satisfied, then problems (5)  Finally, we establish some regularity properties of the solution of problem (4). These regularity results can be easily derived from [18,Chap. 2], [23], and [24] by using a bootstrap argument.
if Assumption 1 holds and there exists a constant such that where c is a positive constant independent of i.
Hereinafter, for the spectral discretization, we assume that is a square or cube. Using the same idea of Nédélec's finite elements (see [27,Sect. 2]), we introduce our discrete spaces.
Let N ≥ 2 be an integer. The velocity discrete space V N is defined as The vorticity discrete space T N is defined as Finally, the pressure discrete spaces M N are defined as We have the following inequality [28]: For continuous functions ϕ and ψ on¯ , we define the discrete scalar product Hereinafter, we suppose that f is continuous on × [0, T]. The full discrete problem is constructed from problem (7)-(8) by using the Galerkin method combined with numerical integration. If The bilinear formsl N (·, · ; ·), b N (·, ·), and t N (·, ·; ·) are defined as follows: From (11) combined with the Cauchy-Schwarz inequality it follows that the bilinear formsl N (·, ·; ·), b N (·, ·). and t N (·, ·; ·) are continuous respectively on ( N is linear and continuous on V N . As a consequence of the exactness property (10), the bilinear forms b(·, ·) and b N (·, ·) coincide on V N × M N . The discrete nonlinear term Z N (·, ·; ·) is defined as We introduce the kernel of the discrete bilinear form b N (·, ·) which is equal to the space of divergence-free polynomials in D N . We also define the discrete kernel of the bilinear form t N (·, · ; ·) We remark that the discrete kernel W N is not included in the continuous kernel W; see [10,Cor 3.2],. So the full discrete problem (12) is reduced as follows: We consider the inf-sup condition proved in [10,Lemma 3.9]. There exists a positive constant β independent of N such that the discrete bilinear form b N (·, ·) satisfies The arguments used to prove the existence of a solution of problems (13) and (12) are exactly the same as those for the continuous problems (9) and (8) where c is a positive constant independent of N and i..
Remark 1 Note that the previous existence result still holds when Z N (·, · ; ·) is replaced by Z(·, · ; ·) in problem (12). In practice, this means that a more precise quadrature formula, exact on P 3N-1 ( ), is used to evaluate the integrals that appear in the treatment of the nonlinear term. The corresponding discrete problem reads: In the same way the discrete reduced problem (13) is written as:

Error estimates
This section is devoted to the proof of the error estimate between the solution of problems (7)- (8) and (15) in dimension two since the proof is difficult in dimension three. This proof is based on the Brezzi-Rappaz-Raviart theorem [15].

2.4])
, where SL is the solution (θ i , v i ) of the following reduced problem: Let We also define the mapping G from the space X = H 0 (curl, ) × K into the dual space of H 0 (div, ) by So we conclude that problem (9) is equivalent to the problem We proceed in the same way for the discrete case. Let X N = W N × K N . Consider the discrete Stokes operator S N such that S N L is the solution (θ i N , v i N ) of the following problem: We also remind from [10] the following properties of the discrete operator S N : • The stability property • If (θ i , v i ) belongs to H s+1 ( ) × H s ( ) 2 for any 1 ≤ i ≤ I and s ≥ 1, then we have the following error: We also define the discrete mapping G N from X N into the dual space of V N by Then we conclude that the problem (13) is equivalent to the problem Let D be a differential operator. We make the following assumption.
We start by proving the following continuity property using the discrete implicit function defined in the theorem of Brezzi, Rappaz, and Raviart [15].

Lemma 1 For any
where c is a positive constant.
Proof Using the Hölder inequality for all r > 2 and s > 2 such that 1 r + 1 s = 1 2 , we obtain Then by the inverse inequality (see [29]) and the fact that the embedding of H(curl, ) into L s ( ) is continuous with a norm bounded by s 1 2 (see [30]), we get inequality (21) when s = ln(N).
Let L the space of linear operators from X into X. Proof Writing we prove the desired result in three steps.
1) The last term in the right-hand side of equality (22) 2) The second term in the right-hand side of equality (22) tends to zero as N tends to infinity. For any (ϑ, ω), we have Then from the property [9, (2.37)] we obtain We conclude that the term SDG(θ , v) · (ϑ, ω) is in H 2 ( ) × H 1 ( ) 2

and satisfies
Finally, from (20) we have (3) Using Assumption 2, take η = (Id + SDG(θ , v)) -1 L for N large enough. Then the quantity in (23) is bounded by 1 2η , which gives the desired result with the norm of the inverse bounded by 2η. Now we prove the following Lipschitz property of the operator S N .

Lemma 3
There exists a constant k > 0 such that for any (θ ,v) ∈ X, Proof For any w N ∈ K N , we have This leads to the desired Lipschitz property using the same idea as in the proof of Lemma 1 and (19).
Proof From equality (17) we derive that Based on (20), we bound the first term in the right-hand side of inequality (25). To bound the second term, we write, thanks to the definition of G and G N , The approximation properties of the operators N-1 and I N (see [28]) give the bound for this last term, which concludes the proof of (24).
Consequently, in the following theorem, we state an error estimate.

Theorem 2
Assume that the data function f belongs to the space L 2 (0, T; H μ ( ) 2 ), μ > 3 2 , and that the solution 2 × H s ( ), s > 1, and satisfies Assumption 2. Then, there exist an integer N * and a positive real number h * such that for all N ≥ N * and |h| ≤ h * , problem (15) has a unique solution. Moreover, this solution satisfies for all 1 ≤ i ≤ I the following error estimate: Proof Using Lemmas 2-4 together with the Brezzi-Rappaz-Raviart theorem (see [15]), for N large enough, we obtain that for all 1 ≤ i ≤ I, problem (16) has a unique solution (θ i N , v i N ), which satisfies Besides, using the discrete inf-sup condition (14), for all 1 ≤ i ≤ I, there exists a unique pressure p i N in M N such that Furthermore, for any q N ∈ M N , having we deduce the estimate for p ip i N L 2 ( ) from (14) and the triangle inequality. (26) is fully optimal for the vorticity and velocity, whereas it is quasioptimal for the pressure.

Resolution algorithm
Considering the result of the error estimate, we will establish numerical tests just in the two-dimensional case on the square =] -1, 1[ 2 . We elaborate the following iterative algorithm to solve problem (12). For simplicity, we omit the indices N .
Step 1. We start by solving the linear Stokes problem: Step 2. We suppose that the (k -1)th iteration (θ i k-1 , v i k-1 , p i k-1 ) is known. Then we solve the problem: We do the iterations until the following condition is satisfied: In the following, we start by presenting the linear system deduced from the discrete problem (27). We build a basis of the discrete spaces T N , V N , and M N .
We consider the Lagrange polynomial ψ p in P N (-1, 1) associated with the nodes p , 0 ≤ p ≤ N , such that where δ pq is the Kronecker symbol. We fix the integer p * between 0 and N equal to N 2 or to (N+1) 2 . We denote the set P * = {0, . . . , N} \ {p * } and consider the polynomial ψ * p ∈ P N-1 ([-1, 1]) such that Then the discrete unknowns θ i N , v i N and a pseudopressurep i N are written as The functionp i does not belong to L 2 0 ( ). However, the real pressure p i N is obtained from the formula The components of the unknowns θ i N , v i Nx , v i Ny , and p i N allow us to form the unknowns vectors denoted by i , V i x , V i y , and P i . Their dimensions are equal to (N -1) 2 , N(N -1), N(N -1), and N 2 -1, respectively. We consider V 0 = (V 0 1 , V 0 2 ), where the components of the vectors V 0 1 and V 0 2 are respectively v 0 Nx ( p , q ) and v 0 Ny ( p , q ) such that v 0 = (v 0 Nx , v 0 Ny ). Consequently, we formulate the discrete problem (27) as the following equivalent linear system: where B T is the transpose of a matrix B. The matrices , and C ω are the same as for the Stokes problem (see [10,Sect. 6]).
The matrices D 1 = (D 11 , D 12 ), D 2 = (D 21 , D 22 ), and N 1 = (N 1 , N 2 ) are made respectively from the terms Z N (θ i Since the global matrix of the linear system (28) is not symmetric, the GMRES method [16] is used for the resolution.

Numerical results
In this section, we start by studying the time convergence. We consider a given solution obtained from the formulas v = curl ϕ and θ = curl v, where ϕ is the stream function. We handle the following two cases.
Case (1): Assume that the steam function ϕ and the pressure p are C ∞ are related to the time and space so that (θ , v; p) is a solution of problem (4): ϕ(t, x, y) = e t sin(πx) sin(πy), p(t, x, y) = e -t (x + y).
Case (2): Assume that the steam function ϕ and the pressure p are less regularly related to the time and space so that (θ , v; p) is a solution of problem (4):    Figures 1(a) and 1(b) represent the convergence in time for the continuous solutions defined in (29) and (30), respectively. We remark that in the two cases (regular solution or less regular solution) the time convergence order is almost equal to 1, which confirms the result of Theorem 2.
In Fig. 2(a), for the solution issued from (29), we present the spectral convergence curves on the vorticity θ in norm H(curl, ), the velocity v in norm H(div, ), and the pressure in norm L 2 ( ). These error curves are provided in semilogarithmic scales, as functions of log(N), for N varying from 5 to 30. As can be expected from Theorem 2, the convergence is exponential for the solution, and the slope for the error curve on the pressure is the same as that for the vorticity and velocity.   H(div, ), and the pressure in norm L 2 ( ) in semilogarithmic scales, as functions of log(N), for the solution issued from (30). We note that the error is much larger for the singular solution (30) than that for the regular solution (29), which confirms the results of Theorem 2. Figure 3 corresponds from top to bottom and left to right to the discrete vorticity, the two components of the discrete velocity, and the discrete pressure for the data f = (f x , f y ) = txy 2 , 0 , v 0 = (0, 0), homogeneous boundary conditions v · n = g = 0 on , and N = 35. Now we handle the influence of the viscosity ν on the number of iterations. We take ξ = 10 -12 and the regular solution issued from (29). Figure 4 presents the number of iterations processed by the algorithm as a function of N as it varies from 5 to 25. This number of iterations is greater when the viscosity ν decreases.

Conclusion
This paper deals with the resolution and implementation of the implicit Euler scheme in time and spectral discretization in space of the nonstationary vorticity velocity pressure formulation of the Navier-Stokes problem with nonstandard boundary conditions. We prove using the Brezzi-Rappaz-Raviart theorem that the new discrete formulation has a unique local solution. In dimension two, we show an optimal error estimate for the vorticity and velocity and a nearly optimal for the pressure,.