Riccardo Milani
email: riccardo.milani@edf.fr

Jérôme Bonelle
email: jerome.bonelle@edf.fr

Alexandre Ern
email: alexandre.ern@enpc.fr

Compatible Discrete Operator schemes for the steady incompressible Stokes and Navier-Stokes equations

Keywords: CDO, Stokes, Navier-Stokes, structure-preserving schemes MSC (2010): 65N08, 76D05, 76M12

Introduction

The Compatible Discrete Operator (CDO) schemes provide a low-order framework which is part of the so-called mimetic or structure-preserving methods. One of the main advantages of the CDO schemes is that they can handle polytopal, nonmatching (cf. left part of Fig. 1) and deformed meshes. Taking advantage of a dual mesh, invisible to the end user, discrete differential operators are carefully designed to satisfy conservation laws and properties typical of their continuous counterparts. This special treatment results in optimal order of convergence in space [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF] (first order for the reconstructed gradient and second for the original variable) and a cell-wise and fully parallelizable building stage ensures good performances of the overall method. Thanks to its flexibility, the CDO framework allows to define the main problem variables on different mesh entities, according to their physical nature. Hence, one can choose to use a cell-, vertex- [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF], edge- [START_REF] Cantin | Vertex-based compatible discrete operator schemes on polyhedral meshes for advection-diffusion equations[END_REF] or face-based [START_REF] Bonelle | Compatible discrete operator schemes on polyhedral meshes for elliptic and Stokes equations[END_REF] scheme. Here, the Stokes and Navier-Stokes equations (NSE) are discretized by means of face-based CDO (CDO-Fb) schemes. In this case, the velocity is defined at faces and cells, and the pressure is defined at cells only. CDO-Fb was introduced initially for the Poisson problem [START_REF] Bonelle | Compatible discrete operator schemes on polyhedral meshes for elliptic and Stokes equations[END_REF] and its key ingredient is a stabilized subcell gradient reconstruction, which can be bridged to the one used in the Hybrid Mixed Mimetic (HMM) framework [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF] and to a generalization of the Crouzeix-Raviart framework [START_REF] Di Pietro | An extension of the Crouzeix-Raviart space to general meshes with application to quasi-incompressible linear elasticity and stokes flow[END_REF] (GCR). A divergence operator is derived from this gradient and it is the tool on which the velocity-pressure coupling hinges. Finally, the discretization of the convection term is inspired by the lowest-order case of the Hybrid High-Order (HHO(k = 0)) method [START_REF] Di Pietro | A discontinuous-skeletal method for advectiondiffusion-reaction on general meshes[END_REF]. The Stokes problem in its curl formulation has been already treated by means of CDO with vertex-and cell-based schemes [START_REF] Bonelle | Analysis of Compatible Discrete Operator schemes for the Stokes equations on polyhedral meshes[END_REF] but here we retain the face-based one.

x c p f,c x f n fc
Let D ⊂ R d , d = 2, 3, be a bounded connected polyhedral domain and denote by ∂D its boundary. We consider the following model problem:

-ν∆ u + χ(u • ∇ )u + ∇ p = f , in D (1a) ∇ • u = 0 , in D (1b)
where ν > 0 is the viscosity, and χ = 0 for the Stokes equations or χ = 1 for the NSE. For the sake of simplicity, homogeneous Dirichlet boundary conditions (BC) are considered. The pressure is uniquely defined by requiring that D p = 0.

Space discretization

A mesh discretizing D is a finite collection C := {c} of nonempty, disjoint, open, polytopal elements of R d , d = 2, 3, usually referred to as cells c. The faces f are assumed to be planar and are gathered in the set F which may be subdivided in two disjoint sets: F b := {f | f ⊂ ∂D} collects the boundary faces and F i := F \ F b the interfaces. One associates with each face f a normal vector n f : if f ∈ F b , n f points outward D and, if f ∈ F i , the direction is chosen arbitrarily. For a mesh entity z = c, f, xz denotes its barycenter and |z| its measure. Consider now a generic cell c. Define the set of faces of the cell c as Fc := {f ∈ F |, f ⊂ ∂c}. For every f ∈ Fc, n fc := ±n f is the normal vector to the face f pointing outward c, the sign depends on the direction chosen for n f . The subpyramid obtained by joining the vertices of f to the barycenter xc of the cell (cf. the central part of Fig. 1) is denoted by p f,c .

Discrete functional spaces and differential operators

Given a generic mesh entity z = c, f, P 0 (z) ≡ R denotes the scalar-valued, zero-th order polynomials defined on z. We denote with πz : L 1 (z) → P 0 (z) ≡ R the L 2 -projection (average): for all s ∈ L 1 (z), πz(s) = z s/ |z|. For vector-valued functions s ∈ [L 1 (z)] d the projection is applied component-wise: πz(s) := (πz(si)) i=1,...,d .

In the CDO-Fb framework the velocity is hybrid, meaning that it has cell-and face-based degrees of freedom (DoFs). Hence, the global velocity space is

U h := × c∈C [P 0 (c)] d × × f∈F [P 0 (f)] d . (2)
An element of U h is denoted by u h := (uc) c∈C , (u f ) f∈F , where, for a generic z = c, f, uz is the z-based DoF. Notice that the value at the interfaces is uniquely defined. The velocity DoFs associated with a cell c are denoted uc

:= uc, (u f ) f∈Fc ∈ Uc := [P 0 (c)] d × × f∈Fc [P 0 (f)] d .
The pressure is defined at cells only: P h := × c∈C Pc p h := (pc) c∈C , where Pc := P 0 (c). In order to take into consideration the velocity BC and the constraint on the pressure average, one also needs

U h,0 := u h ∈ U h | u f = 0 ∀f ∈ F b , P h, * := p h ∈ P h | c∈C |c| pc = 0 .
The right part of Fig. 1 gives an example of local velocity and pressure DoFs for a cell.

Discrete velocity gradient and divergence

For each cell c ∈ C, the discrete local gradient Gc is piecewise constant on the pyramid partition {p f,c } f∈Fc (cf. central part of Fig. 1) and is defined as follows: Gc:

Uc → [P 0 ({p f,c } f∈Fc )] d×d such that for all f ∈ Fc Gc( uc) |p f,c := G 0 c ( uc) + β |f| |p f,c | (u f -uc) -G 0 c ( uc) (x f -xc) ⊗ n fc , (3)
where β > 0 is an arbitrary stability parameter and G 0 c ( u) is a P0-consistent gradient, constant inside each cell and defined as G 0 c ( uc) := 1/ |c| f∈Fc |f| (u f -uc) ⊗ n fc . The definition (3) is the vector-valued version of the gradient introduced in [START_REF] Bonelle | Compatible discrete operator schemes on polyhedral meshes for elliptic and Stokes equations[END_REF]. In the numerical tests, we will use β = 1, which recovers the GCR framework [START_REF] Di Pietro | An extension of the Crouzeix-Raviart space to general meshes with application to quasi-incompressible linear elasticity and stokes flow[END_REF]; the choice β = 1/ √ d gives the HMM one [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF]. For each cell c ∈ C, the discrete velocity divergence Dc : Uc → P 0 (c) is defined as follows

Dc( uc) := trace G 0 c ( uc) = 1 |c| f∈Fc |f| u f • n fc . (4)
Notice that only the face-based DoFs are used (since faces are planar). The discrete velocity divergence is the tool on which the velocity-pressure coupling hinges. This divergence operator can be found also in the HMM framework [START_REF] Droniou | Gradient Schemes for Stokes problem[END_REF].

Discrete advection scheme

The design of the advection scheme is inspired by HHO(k = 0) [START_REF] Di Pietro | A discontinuous-skeletal method for advectiondiffusion-reaction on general meshes[END_REF]. We aim at discretizing the classical advective trilinear form such that D ((w

• ∇)u) • v. Given u h , v h , w h ∈ U h,0 , we use a adv ( w h ; u h , v h ) := 1 2 c∈C f∈Fc |f| (w f • n fc )(u f -uc)(v f + vc) +θ upw f∈F i c∈C f |f| |w f • n f | (u f -uc)(v f -vc) , (5)
where θ upw := 1 in one wants a stabilization by upwinding, or θ upw := 0 for a centered scheme. Suppose, for now, that θ upw = 0. One has:

a adv ( w h ; u h , v h ) + a adv ( w h ; v h , u h ) = - c∈C |c| Dc( wc)uc • vc + f∈F b |f| (w f • n f )u f • v f (6)
obtained by using (4) and by discarding the internal face-defined DoFs since they sum to zero. The boundary DoFs are kept in order to better show that ( 6) is the discrete counterpart of a known integral-by-parts result. Plugging v h = u h into (6) one obtains

a adv ( w h ; u h , u h ) = - 1 2 c∈C |c| Dc( wc)u 2 c . (7)
Supposing there exists µ > 0 such that -1/2 Dc( wc) ≥ µ for all c ∈ C (this is a discrete counterpart of the well-known stability hypothesis of the continuous advection problem), then [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF] proves the coercivity of a adv ( w h ; •, •).

Discrete bilinear form

The discrete counterpart of problem (1) stemming from the CDO-Fb scheme writes: Find

( u h , p h ) ∈ U h,0 × P h, * such that, ∀ v h ∈ U h,0 and ∀q h ∈ P h, * c∈C c {νGc( uc) : Gc( vc) -pcDc( vc)} + χa adv ( u h ; u h , v h ) = c∈C c f • vc , ( 8a
) c∈C -Dc( uc)qc = 0 . (8b)
The Stokes problem (χ = 0 in (1a)) has been analyzed in [START_REF] Droniou | Gradient Schemes for Stokes problem[END_REF].

A static condensation procedure eliminating the cell-based velocity DoFs can be performed in order to reduce the size of the global system, which thus becomes d Card(F) + Card(C). The discarded DoFs are recovered after the solving stage, as a post-processing.

Numerical results

The proposed framework is validated on four test cases, two for the Stokes equations (in 2D and 3D), and two for the NSE (both in 2D). When considering the latter, the nonlinear equations are solved by Picard iterations, and the stopping criterion is evaluated using the cell-based, discrete L 2 -norm of the increment, namely

u k h -u k-1 h C / u k-1 h C < ε, where u h 2 C := c∈C |c| uc 2 2 .
When computing the errors, this velocity norm is considered, as well as the norm of the velocity gradient u h 2 G,C := c∈C |c| G( uc) 2 2 and the discrete L 2 -norm of the pressure p h 2 C := c∈C |c| p 2 c . The resulting error norms used in the analysis are:

erru := u h -π h (u) C π h u C , errgu := u h -π h (u) G,C π h (u) G,C , errp := p h -π h (p) C π h (p) C , (9)
where π h (u) := ((πc(u))c∈C, (π f (u)) f∈F ) and π h (p) := (πc(p)) c∈C . Let nuu (resp. npu) stand for the number of velocity (pressure) unknowns. They will be used to evaluate the orders of convergence in space. We will use the CDO implementation available via Code Saturne [START_REF] Archambeau | Code saturne: A finite volume code for turbulent flows -industrial applications[END_REF], an open-source multipurpose CFD solver developed at EDF R&D. The computations have been performed on a octa-core, Intel i7 laptop with 32GB RAM using PETSc and MUMPS libraries to solve the linear systems. 2D Bercovier-Engelman test case It is proposed in the test case 2.1 of the benchmark [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF]. The sequence of Cartesian meshes (denoted by Hn where n is the number of segments an edge of the domain is divided into) from [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF] and a 2D polygonal family (similarly denoted by Pn, cf. left part of Fig. 2) have been considered. The results are collected in Table 1. [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF]. The meshes used were the Cartesian (Hn) and prismatic with triangular bases (PrTn) sequences proposed in [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF], and one composed of tetrahedra (Tn, the refinement is achieved by dividing each tetrahedra into 8 subtetrahedra). The results are collected in Table 2.

Navier-Stokes equations

Two test cases are considered for the Navier-Stokes equations (χ = 1 in (1a)).

Burggraf flow It consists in a manufactured polynomial solution of the 2D NSE presented in [START_REF] Burggraf | Analytical and numerical studies of the structure of steady separated flows[END_REF]. The centered scheme was considered (θ upw = 0 in (5)). The viscosity is ν = 1/100. About 15 Picard iterations were needed to reach the prescribed tolerance ε = 10 -7 . Two sequences of meshes have been considered: the Cartesian one (Hn) from [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF], and one composed

Figure 1 :

 1 Figure 1: Example of cell compatible with CDO. Left cell with hanging nodes. Center cell with one of its subpyramids p f,c . Right cell with velocity (arrows) and pressure (circle) DoFs

Figure 2 :

 2 Figure 2: Examples of 2D meshes. Left polygonal. Right progressively refined Cartesian

Figure 3 :

 3 Figure3: Lid-driven cavity, vertical and horizontal velocity profiles at the axis of symmetry. Data: CDO H127 (dotted line), CDO H255 (dashed line), CDO H511 (solid line),[START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF] (circle),[START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF] (cross). Left ν = 1/400. Right ν = 1/1000

Table 1 :

 1 Errors for the 2D Bercovier-Engelman test case -Cartesian and polyg-10 -2 0.91 4.07 • 10 -3 2.05 3.36 • 10 -1 2.03 P40 10080 1681 3.13 • 10 -2 0.94 2.29 • 10 -3 2.05 1.91 • 10 -1 2.02 3D Taylor-Green vortex This test case corresponds to Sect. 2.2 of the benchmark

	onal meshes
	mesh	nuu	npu	errgu	order	erru	order	errp	order
	H32	4224	1024 9.15 • 10 -4	-	7.71 • 10 -4	-	1.06 • 10 -1	-
	H64 16640 4096 3.16 • 10 -4 1.55 1.93 • 10 -4 2.02 2.87 • 10 -2 1.98
	H128 66048 16384 1.35 • 10 -4 1.24 4.82 • 10 -5 2.01 7.36 • 10 -3 1.99
	H256 263168 65536 6.41 • 10 -5 1.07 1.21 • 10 -5 2.01 1.85 • 10 -3 1.99
	P10	720	121 9.99 • 10 -2	-	3.33 • 10 -2	-	2.63 • 10 0	-
	P20	2640	441 5.81 • 10 -2 0.83 9.07 • 10 -3 2.00 7.42 • 10 -1 1.96
	P30	5760	961 4.07 •

2D lid-driven cavity It is proposed in the test case 6 of the benchmark [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF]. Two values of the viscosity have been considered: ν = 1/400, 1/1000. Computations have been run on Cartesian meshes with edges divided into 127, 255, and 511 segments. The centered scheme was considered (θ upw = 0 in ( 5)). The prescribed tolerance for the Picard iterations is ε = 10 -7 , less than 25 iterations were needed for ν = 1/400 and less than 30 for 1/1000. In Fig. 3, one can find the plots of the computed vertical and horizontal velocity profiles on the symmetry axes for three Cartesian meshes as well as those from references [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF] and [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF]. Some computations have been run with an upwind scheme (θ upw = 1 in ( 5)) for the advection term, and the results on the velocity profiles were less accurate on the coarser meshes than those obtained with the centered one.