Riccardo Milani

Jérôme Bonelle

Alexandre Ern

Artificial compressibility methods for the incompressible Navier-Stokes equations using lowest-order face-based schemes on polytopal meshes

Keywords: incompressible Navier-Stokes, artificial compressibility, polytopal meshes, lowest-order hybrid schemes MSC 65M12, 65M22, 76D05, 76M10

We investigate artificial compressibility (AC) techniques for the time discretization of the incompressible Navier-Stokes equations. The space discretization is based on a lowest-order face-based scheme supporting polytopal meshes, namely discrete velocities are attached to the mesh faces and cells, whereas discrete pressures are attached to the mesh cells. This face-based scheme can be embedded into the framework of hybrid mixed mimetic schemes and gradient schemes, and has close links to the lowest-order version of hybrid high-order methods devised for the steady incompressible Navier-Stokes equations. The AC timestepping uncouples at each time step the velocity update from the pressure update. The performances of this approach are compared against those of the more traditional monolithic approach which maintains the velocity-pressure coupling at each time step. We consider both first-order and second-order time schemes and either an implicit or an explicit treatment of the nonlinear convection term. We investigate numerically the CFL stability restriction resulting from an explicit treatment, both on Cartesian and polytopal meshes. Finally, numerical tests on large 3D polytopal meshes highlight the efficiency of the AC approach and the benefits of using second-order schemes whenever accurate discrete solutions are to be attained.

Introduction

The goal of this work is to study the accuracy and efficiency of artificial compressibility techniques for the time discretization of the incompressible Navier-Stokes equations. These equations are encountered in a wide range of industrial applications, ranging from aeronautics to the simulations of flows in micro-fractures, to cite two salient examples. In this work, the space discretization is realized by means of lowest-order face-based schemes supporting classical (simplicial or Cartesian) as well as general (polytopal) meshes. The support of polytopal meshes is motivated by the applications, especially in the industrial context. Indeed, such meshes often alleviate substantially the burden of mesh generation resulting from the complexity of the geometry (as, e.g., the shape of industrial equipments) and/or the heterogeneity of the physical properties (requiring local mesh refinements producing hanging nodes).

The unsteady incompressible Navier-Stokes equations read as follows:

     ∂u ∂t -ν∆ u + (u • ∇ )u + ∇ p = f , ∇ • u = 0 , in Ω × (0, T ) , (1)
together with suitable boundary and initial conditions on the velocity. Here, Ω ⊂ R d , d = 2, 3, is an open, polytopal, bounded, Lipschitz domain and T > 0 denotes the observation time. The unknowns are the velocity u and the pressure p (by convention, vector-valued quantities are underlined). The problem data are the viscosity ν > 0, the mass density is normalized to unity, and the body force is denoted by f . For simplicity, we consider smooth solutions to [START_REF] Aghili | Hybridization of mixed high-order methods on general meshes and application to the Stokes equations[END_REF], and we enforce Dirichlet boundary conditions over the whole boundary ∂Ω at all times. Recall that the first equation in [START_REF] Aghili | Hybridization of mixed high-order methods on general meshes and application to the Stokes equations[END_REF] expresses the momentum balance in the flow, and the second equation the mass conservation. Let us briefly describe the time-marching schemes considered in this work without introducing any space discretization yet. These schemes are either of monolithic form, thereby requiring to solve a saddle point problem at each time step, or they uncouple the velocity and the pressure at each time step by means of an artificial compressibility technique. Let ∆t > 0 be the time step (taken to be constant for simplicity) and let (u n , p n ) denote the approximate solution at the discrete time node t n := n∆t for all n ≥ 1. The monolithic approach is the traditional time-marching scheme for the incompressible Navier-Stokes equations and it reads as follows: For all n ≥ 1, given u n-1 from the initial condition or the previous time node, find (u n , p n ) by solving the saddle point problem in Ω,

     u n -u n-1 ∆t -ν∆ u n + (u n • ∇ )u n + ∇ p n = f n := f (t n ) , ∇ • u n = 0 .
(

) 2
The monolithic approach is well-known and is considered here as the reference time-stepping scheme. The only approximation introduced in (2) with respect to [START_REF] Aghili | Hybridization of mixed high-order methods on general meshes and application to the Stokes equations[END_REF] is the discretization of the velocity time-derivative (here, by means of the implicit Euler scheme to fix the ideas). In contrast, the artificial compressibility (AC) approach introduces an additional approximation at each time step in that the velocity is first updated and the pressure is then corrected. In its first-order form, the scheme reads as follows: For all n ≥ 1, given (u n-1 , p n-1 ) from the initial condition or the previous time node, find (u n , p n ) by solving a parabolic-like problem on the velocity and then updating the pressure:

u n -u n-1 ∆t -ν(∆ u n + η∇ ∇ • u n ) + (u n • ∇ )u n = f n -∇ p n-1 , ( 3a
)
p n = p n-1 -νη∇ • u n , ( 3b
)
where η > 0 is a user-defined adimensional parameter. The AC method was introduced in the late sixties [START_REF] Chorin | Numerical Solution of the Navier-Stokes Equations[END_REF][START_REF] Temam | Sur l'Approximation de la Solution des Equations de Navier-Stokes par la Méthode des Pas Fractionnaires (I)[END_REF] in the Western literature, although the seminal ideas can be traced independently in the Russian literature as well [START_REF] Vladimirova | Numerical calculation of the symmetrical flow of viscous incompressible liquid around a plate[END_REF][START_REF] Ladyzhenskaya | The mathematical theory of viscous incompressible flow[END_REF][START_REF] Yanenko | The method of fractional steps[END_REF]. For a recent analysis of the method, the reader is referred to [START_REF] Guermond | High-order time stepping for the incompressible Navier-Stokes equations[END_REF]. Notice the appearance of the additional grad-div term in (3a), whereas (3b) shows that the time discrete velocity field u n is no longer divergence-free. Notice also that an approximation to the initial pressure is needed. Furthermore, we observe that higher-order versions of the AC scheme are available; see [START_REF] Guermond | High-order time stepping for the incompressible Navier-Stokes equations[END_REF][START_REF] Guermond | High-order adaptive time stepping for the incompressible Navier-Stokes equations[END_REF]. The computational study performed in this work considers both the monolithic and the AC approaches, either in the above form which is first-order accurate in time, or in second-order form (as described below). Moreover, we examine either the implicit treatment of the convection term (as above) leading to a nonlinear problem to be solved at each time step, or a (semi-)explicit treatment leading to a linear problem to be solved at each time step. The natural choice when using the AC scheme is to resort to an explicit treatment since the goal is to reduce as much as possible the computational cost per time step. However, for completeness, we also consider the implicit treatment of the convection term for the first-order AC schemes, but not for the second-order ones. Indeed, in this latter case, there are two nonlinear problems to be solved at each time step, which is deemed to be too expensive.

Space discretization by the CDO-Fb scheme

In this section, we recall the CDO-Fb scheme for the space discretization of the steady incompressible Navier-Stokes equations introduced in [START_REF] Bonelle | Compatible Discrete Operator schemes for the steady incompressible Stokes and Navier-Stokes equations[END_REF]. In weak form, and considering homogeneous Dirichlet conditions for simplicity (and fixing the pressure mean value to zero), one seeks

(u, p) ∈ H 1 0 (Ω) × L 2 0 (Ω) (standard notation for Lebesgue and Sobolev spaces is employed) such that νa(u, v) + t(u; u, v) + b(v, p) = l(v) , ∀v ∈ H 1 0 (Ω) , b(u, q) = 0 , ∀q ∈ L 2 0 (Ω) , (4)
with the bilinear and trilinear forms

a(u, v) := Ω ∇ u : ∇ v , b(v, p) := - Ω (∇ •v)p , t(w; u, v) := Ω (w•∇ )u •v , ( 5
)
and the linear form l(v)

:= Ω f •v. x c p f,c
x f n fc = vector (velocity) = scalar (pressure)

Figure 1: Mesh, notation and DoFs. Left: Example of a polyhedral cell with a hanging node. Center: The barycenter x c of the cell c and a face f ∈ F c with its barycenter, x f , its normal n fc and its sub-pyramid p f,c . Right: Full set of velocity and pressure DoFs for the considered cell (dashed arrows are used for velocity DoFs on hidden faces).

Mesh entities and degrees of freedom

The mesh discretizing Ω is a finite collection C := {c} of nonempty, disjoint, open, polytopal subsets of R d , usually referred to as cells. The mesh faces are assumed to be planar and are gathered in the set F := {f} which can be subdivided in the two disjoint sets F b := {f | f ⊂ ∂Ω} collecting the mesh boundary faces and F i := F\F b collecting the mesh interfaces. One associates with each mesh face f a normal vector n f such that if f ∈ F b , n f points outward Ω and, if f ∈ F i , the direction is chosen arbitrarily and fixed once and for all. For a generic mesh entity z = c or z = f, x z denotes its barycenter, |z| its measure and h z := max x 1 ,x 2 ∈z |x 1 -x 2 | its diameter. Moreover, h := max c∈C h c is called the size of the mesh. We consider the shape-regularity setting of [START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF] for polytopal mesh families, and we additionally require that every mesh entity z = c or z = f is star-shaped with respect to its barycenter x z . The faces composing the boundary of a cell c ∈ C are collected in the set F c := {f ∈ F | f ⊂ ∂c}. For each face f ∈ F c , a unit normal vector pointing outward c is denoted by n fc = ±n f , depending on the orientation chosen for n f . We are going to make use of a subdivision of the cell c as P c := {p f,c } f∈Fc , where the subsets p f,c are the nonempty, disjoint subpyramids (or subtriangles if d = 2) obtained by considering the cell barycenter x c as apex, and a face f ∈ F c as basis. For a generic mesh entity z = c or z = f, P q (z) (resp., P q (z) and P q (z)) is composed of the restriction to z of the scalar-valued (resp., R d -valued and R d×d -valued) polynomials of degree at most q. Moreover, for a collection Z of mesh entities, P q (Z), P q (Z), and P q (Z) refer to piecewisepolynomial functions; for instance, P 1 (F c ) := × f∈Fc P 1 (f). Notice also that for z = c or z = f, P 0 (z) ≡ R.

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

U h := × f∈F P 0 (f) × × c∈C P 0 (c) . ( 6
)
A generic element of U h is denoted by v h := (v f ) f∈F , (v c ) c∈C . Notice that the velocity DoFs at the interfaces are uniquely defined. Moreover, the velocity DoFs associated with a generic cell c ∈ C are denoted by

v c := (v f ) f∈Fc , v c ∈ U c := × f∈Fc P 0 (f) × P 0 (c) . ( 7
)
The pressure has only cell-based DoFs, so that the global pressure space is

P h := × c∈C P c p h := (p c ) c∈C , P c := P 0 (c) . ( 8
)
In order to account for the velocity boundary conditions and the constraint on the pressure average, we consider the subspaces

U h,0 := { v h ∈ U h | v f = 0 ∀f ∈ F b }, P h, * := {p h ∈ P h | c∈C |c| p c = 0}. (9)
Finally, for z = c or z = f, π z (resp., π z ) denotes the L 2 -orthogonal projection onto P 0 (z) (resp., P 0 (z)). The L 2 -orthogonal projection onto P h is defined cellwise so that π h (q) := (π c (q| c )) c∈C for all q ∈ L 1 (Ω), whereas the projection onto the hybrid discrete space

U h is defined such that π h (v) := (π f (v| f )) f∈F , (π c (v)) c∈C for all v ∈ H s (Ω), s > 1 2 .
Similarly, for the local hybrid space U c , we set π c (v) := (π f (v| f )) f∈Fc , π c (v) for all v ∈ H s (c).

Discrete diffusion-like bilinear form

Consider a cell c ∈ C and the velocity DoFs U c associated with c. We define locally in the mesh cell c a velocity gradient reconstruction operator, G c : U c → P 0 (P c ), i.e., for all v c ∈ U c , G c ( v c ) is a piecewise constant tensor-valued field on the subpyramids collected in P c . In the context of CDO schemes, this type of operator was introduced in [START_REF] Bonelle | Compatible Discrete Operator schemes on polyhedral meshes for elliptic and Stokes equations[END_REF]; see also [START_REF] Bonelle | Low-order reconstruction operators on polyhedral meshes: application to compatible discrete operator schemes[END_REF]. It is composed of a consistent part, G 0 c , plus a stabilization part, G s c , so that

G c ( v c ) := G 0 c ( v c ) + G s c ( v c ) , (10)
where G 0 c ( v c ) ∈ P 0 (c) ⊂ P 0 (P c ) and G s c ( v c ) ∈ P 0 (P c ) are defined as follows:

G 0 c ( v c ) := 1 |c| f∈Fc |f| (v f -v c ) ⊗ n fc , G s c ( v c )| p f,c := θ |f| |p f,c | (v f -v c ) -G 0 c ( v c )(x f -x c ) ⊗ n fc , ∀f ∈ F c . ( 11
)
The positive stabilization parameter θ is user-defined (with the only requirement of being positive): choosing θ := 1 recovers the generalized Crouzeix-Raviart scheme from [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], whereas θ := 1 √ d leads to the HFV scheme from [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF].

The global version of the gradient reconstruction operator, G h : U h → × c∈C P 0 (P c ), is defined cellwise so that, for any v h ∈ U h and any c ∈ C, we have

G h ( v h )| c := G c ( v c ).
With the above operators in hand, the discrete diffusion-like bilinear form a h :

U h × U h → R is defined as follows: a h ( u h , v h ) := Ω G h ( u h ) : G h ( v h ) = c∈C c G c ( u c ) : G c ( v c ) . ( 12
)
The above-defined gradient reconstruction operators enjoy several important properties. First, G c is consistent for affine functions, meaning that

G c π c (v) = ∇v , ∀v ∈ P 1 (c) . ( 13
)
(More precisely, we have G 0 c π c (v) = ∇v and G s c π c (v) = 0.) Moreover, the consistent and stabilization parts of the reconstructed gradient are

L 2 (c)-orthogonal, c G 0 c ( v c ) : G s c ( v c ) = 0 , ∀ v c ∈ U c . ( 14
)
Finally, and most importantly, the following stability and boundedness properties hold true: There is δ > 0 such that for all c ∈ C and all v c ∈ U c ,

δ v c 2 1,c ≤ G c ( v c ) 2 L 2 (c) ≤ δ -1 v c 2 1,c , (15)
with the discrete H 1 -like seminorm defined on

U c such that v c 2 1,c := f∈Fc h -1 c |f| |v f -v c | 2 .
The lower bound [START_REF]Analysis of Compatible Discrete Operator schemes for the Stokes equations on polyhedral meshes[END_REF] is the main reason for introducing the stabilization part of the reconstructed gradient. Notice that v h 2 1,h := c∈C v c 2 1,c defines a norm on U h,0 .

Discrete velocity-pressure coupling

The discrete velocity-pressure coupling hinges on a discrete divergence operator D h :

U h → P h that is defined cellwise. For all c ∈ C, the local discrete divergence operator D c : U c → P c is such that for all v c ∈ U c , D c ( v c ) := tr G 0 c ( v c ) = 1 |c| f∈Fc |f| (v f -v c ) • n fc . ( 16
)
The global operator is then defined by setting

D h ( v h )| c = D c ( v c ) for all v h ∈ U h and all c ∈ C. Notice that D h ( v h ) ∈ P h, * for all v h ∈ U h,0 .
With the above divergence operator in hand, the discrete bilinear form b h : U h × P h → R handling the velocity-pressure coupling is defined as follows:

b h ( v h , q h ) := - Ω D h ( v h )q h = - c∈C |c| D c ( v c )q c . ( 17
)
The same discrete bilinear form is found in HMM schemes [START_REF] Droniou | Gradient Schemes for Stokes problem[END_REF][START_REF]Benchmark: Two hybrid mimetic mixed schemes for the lid-driven cavity[END_REF] and the lowest-order HHO scheme [START_REF] Aghili | Hybridization of mixed high-order methods on general meshes and application to the Stokes equations[END_REF][START_REF] Di Pietro | A discontinuous skeletal method for the viscosity-dependent Stokes problem[END_REF]. A crucial property is the following discrete inf-sup condition: There is β > 0, only depending on the mesh shape-regularity, such that inf

q h ∈P h, * sup v h ∈ U h,0 |b h ( v h , q h )| q h h v h 1,h ≥ β , ( 18
)
where v h 1,h is defined above whereas q h

2 h := q h 2 L 2 (Ω) = c∈C |c| |q c | 2
for all q h ∈ P h . Another important property of the discrete divergence operator is its commuting property with the L 2 -orthogonal projection in the sense that D

c π c (v) = π c (∇ • v) for all c ∈ C and all v ∈ H 1 (c).

Discrete convection operator

Let us finally devise a discrete CDO-Fb counterpart of the trilinear form t defined in [START_REF] Arioli | An iterative generalized Golub-Kahan algorithm for problems in structural mechanics[END_REF]. To this purpose, we define the discrete trilinear form

t h : U h × U h × U h → R such that t h ( w h ; u h , v h ) := c∈C t c ( w c ; u c , v c ) + f∈F b t f (w f ; u f , v f ) , ( 19
) with t c ( w c ; u c , v c ) := 1 2 f∈Fc |f| (w f • n fc )(u f -u c ) • (v f + v c ) , t f (w f ; u f , v f ) := |f| (w f • n fc ) -u f • v f , (20)
where (x) -:= 1 2 (|x| -x) denotes the negative part of any real number x ∈ R. The trilinear form t f is related to the weak enforcement of inflow boundary conditions and vanishes whenever at least one of the arguments is in U h,0 . The discrete trilinear form t h is inspired from the lowest-order HHO method for scalar advection-diffusion equations introduced in [START_REF] Di Pietro | A discontinuous-skeletal method for advection-diffusion-reaction on general meshes[END_REF]; see also [START_REF] Botti | A Hybrid High-Order method for the incompressible Navier-Stokes equations based on Temam's device[END_REF]Remark 9] for the Navier-Stokes equations.

The first important property of the discrete trilinear form t h is its positivity and skewsymmetry. Indeed, assume that w h ∈ U h is discretely divergence-free, i.e., D c ( w c ) = 0 for all c ∈ C. Then, one readily verifies that

t h ( w h ; u h , u h ) ≥ 0 , ∀ u h ∈ U h . ( 21
)
Moreover, if additionally the normal component of w h vanishes at the boundary, namely

w f •n f = 0 for all f ∈ F b , then t h ( w h ; •, •) is skew-symmetric, i.e. t h ( w h ; u h , u h ) = 0 , ∀ u h ∈ U h . ( 22
)
Property ( 21) is crucial in the context of the Navier-Stokes equations since it is instrumental to establish the dissipativity of the discrete problem. Notice also that ( 22) is reminiscent of the so-called Temam's trick on the discrete trilinear form. A further relevant property of the discrete trilinear form t h is its limit-conformity. This notion is described, e.g., in [START_REF] Eymard | Family of convergent numerical schemes for the incompressible Navier-Stokes equations[END_REF], and is verified for the above CDO-Fb setting in [START_REF] Milani | Compatible Discrete Operator schemes for the unsteady incompressible Navier-Stokes equations[END_REF]Remark 2.50]; see also [START_REF] Di Pietro | A Hybrid High-Order method for the steady incompressible Navier-Stokes problem[END_REF]Prop. 6].

Remark 1 (Upwinding). An upwinding stabilization can be considered by adding to [START_REF]Finite Volumes for Complex Applications VIII -Methods and Theoretical Aspects[END_REF] the discrete trilinear form

t u h ( w h ; u h , v h ) := 1 2 f∈F i c∈C f |f| |w f • n fc | (u f -u c )•(v f -v c
). This option is not considered further in the present work.

Time discretization by monolithic and artificial compressibility schemes

We present in this section the fully discrete problems obtained by using either a monolithic or an AC scheme for the time discretization and by considering either an implicit or an explicit treatment of the convection operator. Moreover, we consider both first-order and second-order time schemes. For a hybrid velocity field v h := (v f ) f∈F , (v c ) c∈C ∈ U h , we denote its cell-based components by v C := (v c ) c∈C . Then, to discretize the time-derivative of the velocity, we only use the cellbased components and consider the mass bilinear form such that

m(u C , v C ) := c∈C |c| u c •v c . ( 23
)
We also use the cell-based components to discretize the source term and we set

l n (v C ) := c∈C c f n •v c for all n ≥ 1.
For any hybrid velocity field v h ∈ U h,0 , we define its discrete kinetic energy as

E kin,h ( v h ) := 1 2 v C 2 L 2 (Ω) = 1 2 c∈C |c| |v c | 2 2 . ( 24
)
3.1 First-order schemes

Monolithic scheme

Let us first consider an implicit treatment of the convection term which is dealt with by means of a Picard algorithm. Then, the first-order monolithic time-stepping scheme reads as follows.

For all n = 1, . . . , N , iterate on k ≥ 1 until convergence:

Find ( u n,k h , p n,k h ) ∈ U h,0 × P h, * such that      1 ∆t m(u n,k C -u n-1 C , v C ) + νa h ( u n,k h , v h ) + t h ( u n,k-1 h ; u n,k h , v h ) + b h ( v h , p n,k h ) = l n (v C ) , b h ( u n,k h , q h ) = 0 . ( 25
)
for all v h ∈ U h,0 and all q h ∈ P h, * . Here, u n-1 h denotes the solution given by the Picard algorithm at the time step n -1. The time-stepping is initialized with the initial condition u 0 h := π h (u 0 ). Moreover, at each time step, the Picard iterations have to be initialized: a suitable choice is to take the solution at the previous time step, i.e., u n,0 25) and using simple arithmetic manipulations together with t h ( u n,k-1 h ; u n,k h , u n,k h ) = 0 owing to [START_REF] Cantin | An edge-based scheme on polyhedral meshes for vector advectionreaction equations[END_REF], we obtain the following discrete kinetic energy balance:

h := u n-1 h . Taking v h = u n,k h in (
E kin,h ( u n,k h ) -E kin,h ( u n-1 h ) + E kin,h ( u n,k h -u n-1 h ) + ν∆t G h ( u n,k h ) 2 L 2 (Ω) = ∆t l n (u n,k C ) . ( 26
)
Notice that the third and fourth terms on the left-hand side are non-negative, i.e., they only contribute to energy dissipation. Notice also that ( 26) remains valid as k → ∞.

With an explicit treatment of the convection term, the first-order monolithic time-stepping scheme reads as follows: For all n = 1, . . . , N , find (

u n h , p n h ) ∈ U h,0 × P h, * such that    1 ∆t m(u n C -u n-1 C , v C ) + νa h ( u n h , v h ) + b h ( v h , p n h ) = l n (v C ) -t h ( u n-1 h ; u n-1 h , v h ) , b h ( u n h , q h ) = 0 , ( 27
)
for all v h ∈ U h,0 and all q h ∈ P h, * . Unfortunately, a dissipative discrete kinetic energy balance cannot be derived for [START_REF] Di Pietro | A discontinuous-skeletal method for advection-diffusion-reaction on general meshes[END_REF] because the term t h ( u n-1 h ; u n-1 h , u n h ) cannot be given an a priori sign. Thus, the scheme ( 27) is subject to a first-order CFL restriction on the time step for its stability. This condition will be investigated numerically in Section 4.

AC scheme

Let us first consider an implicit treatment of the convection term. As in the monolithic approach, we consider Picard iterations to deal with the resulting nonlinearity. Then, the first-order AC time-stepping scheme reads as follows. For all n = 1, . . . , N , iterate on k ≥ 1 until convergence:

Find u n,k h ∈ U h,0 such that for all v h ∈ U h,0 , 1 ∆t m(u n,k C -u n-1 C , v C ) + ν a h ( u n,k h , v h ) + ηd h ( u n,k h , v h ) + t h ( u n,k-1 h ; u n,k h , v h ) = l n (v C ) -b h ( v h , p n,k-1 h
) , [START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF] and then set

p n,k h = p n,k-1 h -νηD h ( u n,k h ) , (29)
recalling that η > 0 is a user-defined parameter (choices are discussed in Section 4). In (28), we introduced the discrete div-div bilinear form such that

d h ( u h , v h ) := Ω D h ( u h )D h ( v h ) = c∈C |c| D c ( u c )D c ( v c ) . ( 30
)
Notice that in ( 28)-( 29), we took advantage of the Picard iterations to update progressively the pressure. Letting k → ∞, we observe that the velocity and pressure updates stemming from the monolithic and the AC schemes are the same. Differently from the monolithic coupling, the initialization of the Picard iteration requires to specify the velocity and the pressure, which are here taken from the previous time step or the initial condition.

With an explicit treatment of the convection term, the first-order AC scheme reads as follows: For all n = 1, . . . , N , find (u n h , p n h ) ∈ U h,0 × P h, * such that for all v h ∈ U h,0 ,

1 ∆t m(u n C -u n-1 C , v C ) + ν a h ( u n h , v h ) + ηd h ( u n h , v h ) = l n (v C ) -b h ( v h , p n-1 h ) -t h ( u n-1 h ; u n-1 h , v h ) , ( 31
)
and then set

p n h = p n-1 h -νηD h ( u n h ) . ( 32
)
As for the monolithic scheme, the explicit AC scheme is subject to a first-order CFL restriction on the time step for its stability. This condition will be investigated numerically in Section 4. Remark 2 (Initial pressure). As is classical with the AC technique, an initial pressure needs to be specified. In all the numerical experiments reported in this work, the initial pressure is indeed known. If this were not the case, a possibility proposed in [START_REF] Guermond | High-order time stepping for the incompressible Navier-Stokes equations[END_REF] (see also [START_REF] Ern | Finite Elements III[END_REF]Remarks 74.4 and 75.7]) is to recover the initial pressure by solving the following steady problem:

     ∆p 0 = ∇ • f 0 , ∂p 0 ∂n | ∂Ω = f 0 -(-ν∆ u 0 + (u 0 • ∇ )u 0 ) | ∂Ω •n ∂Ω . ( 33
)
Obviously, the initial pressure is zero is the fluid is initially at rest and the source term is divergence-free and has a zero normal component on the boundary.

Second-order schemes

Monolithic scheme

One standard technique is to use a second-order backward differentiation formula (BDF2) to devise second-order time-stepping schemes within the monolithic approach. In general, BDF2 is employed for every time step n ≥ 2, and an implicit Euler step can still be considered for n = 1. With an implicit treatment of the convection term by means of a Picard algorithm, the second-order monolithic scheme reads as follows. For all n = 2, . . . , N , iterate on k ≥ 1 until convergence:

Find ( u n,k h , p n,k h ) ∈ U h,0 × P h, * such that            1 2∆t m(3u n,k C -4u n-1 C + u n-2 C , v C ) + νa h ( u n,k h , v h ) + t h ( u n,k-1 h ; u n,k h , v h ) + b h ( v h , p n,k h ) = l n (v C ) , b h ( u n,k h , q h ) = 0 , ( 34
)
for all v h ∈ U h,0 and all q h ∈ P h, * . Recall that u n-1 h (resp. u n-2 h ) denotes the solution given by the Picard algorithm at the time step n -1 (resp. n -2). The initialization of the iterative procedure is done with an approximation of the solution, for instance u n-1 h (as for the first-order scheme), or by using the second-order extrapolation formula (2 u n-1 h -u n-2 h ). Finally, for the BDF2 time discretization, it is possible to derive by means of classical algebraic manipulations (detailed, e.g., in [START_REF] Ern | Finite Elements III[END_REF]Lemma 68.1] for the heat equation) a time-discrete kinetic energy balance with dissipation, in the same spirit as in [START_REF] Delcourte | Développement de méthodes de volumes finis pour la mécanique des fluides[END_REF].

With an explicit treatment of the convection term, the second-order monolithic scheme reads as follows: For all n = 2, . . . , N , find (

u n h , p n h ) ∈ U h,0 × P h, * solving          1 2∆t m(3u n C -4u n-1 C + u n-2 C , v C ) + νa h ( u n h , v h ) + b h ( v h , p n h ) = l n (v C ) -2t h ( u n-1 h ; u n-1 h , v h ) -t h ( u n-2 h ; u n-2 h , v h ) , b h ( u n h , q h ) = 0 , ( 35
)
for all v h ∈ U h,0 and all q h ∈ P h, * . Notice that a second-order extrapolation formula is used for the convective term on the right-hand side; here, we applied the extrapolation formula to the operator, but it is also possible to consider applying the extrapolation formula to the discrete velocity and then form the convection operator. Finally, as for the first-order time discretization, a dissipative kinetic energy balance is not available with an explicit treatment of the convection term, and a CFL restriction on the time step is required for stability.

AC scheme

It has been shown in [START_REF] Guermond | High-order time stepping for the incompressible Navier-Stokes equations[END_REF] that arbitrary order in time can be achieved by combining the AC method with a bootstrap technique. In order to obtain the k-th order, k linear systems similar to (3) have to be solved. An alternative technique to reach arbitrary order in time, which is also discussed in [START_REF] Guermond | High-order time stepping for the incompressible Navier-Stokes equations[END_REF], is based on a defect correction procedure, leading to a similar computational cost. To avoid the proliferation of variants, we focus only on the bootstrap technique.

As motivated in the introduction, we only consider an explicit treatment of the convection term. The second-order AC time-stepping scheme reads as follows: For all n ≥ 1, find

( u n 1,h , p n 1,h ) ∈ U h,0 × P h, * such that for all v h ∈ U h,0 ,            1 ∆t m(u n 1,C -u n-1 1,C , v C ) + ν a h ( u n 1,h , v h ) + ηd h ( u n 1,h , v h ) = l n (v C ) -b h ( v h , p n-1 1,h ) -t h ( u n-1 1,h ; u n-1 1,h , v h ) , p n 1,h = p n-1 1,h -νηD h ( u n 1,h ), δp n 1,h := p n 1,h -p n-1 1,h . (36a) (36b)
Then, for all n ≥ 2, find (

u n 2,h , p n 2,h ) ∈ U h,0 × P h, * such that for all v h ∈ U h,0 ,            1 2∆t m(3u n 2,C -4u n-1 2,C + u n-2 2,C , v C ) + ν a h ( u n 2,h , v h ) + ηd h ( u n 2,h , v h ) = l n (v C ) -b h ( v h , p n-1 2,h + δp n 1,h ) -2t h ( u n-1 2,h ; u n-1 2,h , v h ) -t h ( u n-2 2,h ; u n-2 2,h , v h ) , p n 2,h = p n-1 2,h + δp n 1,h -νηD h ( u n 2,h ) , (36c) (36d)
with the following initialization choice: u 1 2,h := u 1 1,h , u 0 2,h := u 0 1,h := u 0 h for the velocity and p 1 2,h := p 1 1,h for the pressure.

Numerical results

In this section, we present numerical results to assess the performances of the AC method together with the CDO-Fb discretization. We first consider the 2D Taylor-Green vortex test case with the goal to verify the convergence rates in time for the first-and second-order timestepping schemes and both treatments of the convection term. We also use this test case to study the CFL restriction on the time step in the case of an explicit treatment of the convection term. Then, we consider a modified 3D Taylor-Green vortex to compare all of the above strategies in terms of efficiency, that is, we compare the reached error thresholds on the velocity and the pressure with the needed CPU time. All the results are computed with the CDO module available in the open-source, single-phase CFD solver Code_Saturne [START_REF] Archambeau | Code Saturne: A Finite Volume Code for Turbulent flows -Industrial Applications[END_REF]. For the 2D test cases, the linear systems are solved with a sparse direct solver available in the MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF] library. For the 3D test cases, iterative solvers are considered and further discussed in Section 4.2.

2D Taylor-Green vortex

The Taylor-Green vortex [START_REF] Taylor | Mechanism of the production of small eddies from large ones[END_REF] is a well-known 2D test case usually considered to evaluate the performances of an unsteady Navier-Stokes solver. In fact, it constitutes an analytic solution of the 2D Navier-Stokes given by

        
u TGV (x, y) := exp(-2νt) sin(x) cos(y) -cos(x) sin(y) , p TGV (x, y) := 1 4 exp(-4νt)(cos(2x) + cos(2y)) .

(

) 37
The computational domain is Ω := [0, 2π] 2 , whereas the observation time T depends on the numerical experiment. The viscosity ν is uniform but different values are considered in order to modify the Reynolds number Re. Since the reference length and velocity for this test case can be set to L := 1 and U := 1, we obtain Re = 1 ν . Cartesian and polygonal meshes are considered.

Convergence results

The convergence in time is first verified for all the schemes presented in Section 3. Let v(t, x) (resp. q(t, x)) be the analytical velocity (resp. pressure), and { v n h } n=1,...,N ∈ [ U h,0 ] N (resp. {q n h } n=1,...,N ∈ [P h, * ] N ) its space-time discretization. The following space-time error norms are used:

π h (v) -{ v n h } n 2 2 (L 2 ) := N n=1 ∆t π h (v(t n , •)) -v n h 2 C := N n=1 ∆t c∈C |c| |π c (v(t n , •)) -v n c | 2 , π h (v) -{ v n h } n 2 2 (H 1 ) := N n=1 ∆t c∈C f∈Fc |p f,c | G c π c (v(t n , •)) -v n c | p f,c 2 , π h (q) -{q n h } n 2 2 (L 2 ) := N n=1 ∆t π h (q(t n , •)) -q n h 2 C := N n=1 ∆t c∈C |c| (π c (q(t n , •)) -q n c ) 2 . ( 38
)
These errors are then normalized by the value obtained for the corresponding norm of the exact solution.

Figure 2 shows the results for first-(top row) and second-order time-stepping schemes (bottom row) on a Cartesian mesh composed of 512 2 cells. The Reynolds number is first set to Re = 1 (i.e., we take ν = 1), and the AC parameter is set to η = 10. We consider both the implicit and the explicit treatments of the convection term, and observe that in the present setting, no instability is observed for the considered values of the time step. This favorable situation is attributed to the moderate value of the Reynolds number, and we will see below (see Section 4.1.2) that a CFL stability restriction needs to be enforced for larger Reynolds numbers. The results presented in Figure 2 show that the monolithic and AC approaches lead to very similar errors, both for the velocity and the pressure. The convection treatment has an effect only on the pressure error (right column): the implicit treatment (squares) is more accurate than the explicit one (circles). Indeed, due to the peculiar construction of the Taylor-Green vortex for which the pressure gradient and the convection term balance out, a good approximation of the convection term leads also to an accurate approximation of the pressure.

Additional computations with higher Reynolds numbers are performed to study the dependence of the results delivered by the AC method on the parameter η. The results are reported in Table 1. As expected, the predictions of the AC method get closer to those of the monolithic approach when larger values of η are used. However, using larger values also leads to higher computational costs because it affects the conditioning of the linear system through the presence of the grad-div term. The results altogether indicate that an appropriate choice is to set η = 10Re. Indeed, this choice leads to a satisfactory trade-off between the accuracy and the conditioning of the linear systems.

CFL conditions with explicit convection

We investigate numerically the stability of the monolithic and AC approaches when an explicit treatment of the convection term is used. For these tests, we consider (i) three Reynolds numbers, Re ∈ {200, 500, 700}, (ii) an observation time T such that T Re = 10 4 , (iii) η = 10 Re whenever the AC method is considered. We flag a computation as having diverged if, for some n ≥ 1, we have

E kin,h ( u n h ) > 1.1 E kin,h ( u 0 h ) = 1.1 E kin,h ( π h (u 0 )) , ( 39
)
10 -3 Table 1: 2D Taylor-Green vortex. Space-time velocity and pressure errors. Explicit convection. AC method with η ∈ {Re, 10Re, 100Re}. The errors obtained with the monolithic approach are given in parenthesis as reference. First group of columns: Re ≈ 33, Cartesian mesh composed of 128 2 cells, T = 40, ∆t = T 64 = 0.625. Second group of columns: Re = 100, Cartesian mesh composed of 512 2 cells, T = 120, ∆t = T 64 = 1.875. where the definition of the discrete kinetic energy E kin,h (•) is given in [START_REF] Cockburn | Local discontinuous Galerkin methods for the Stokes system[END_REF]. Notice that the solution (37) goes exponentially towards 0 with respect to time. Thus, a failure to satisfy ( 39) is a symptom of stability issues. We denote by ∆t s the critical time-step, that is, the largest ∆t ensuring that the computation does not diverge. We seek a resolution of 1%, meaning that the gap between ∆t s and the smallest ∆t leading to divergence is less than 1% of ∆t s . Whenever ∆t > ∆t s , we define the divergence time, T d , as the smallest t n for which (39) is satisfied. Two meshes have been considered, a Cartesian mesh and a polygonal mesh (illustrated in Fig. 3), which have, respectively, 16,384 and 15,129 cells. In Fig. 4, we plot the divergence time, T d , for different time-step values. We observe that the divergence time sharply increases as the time step is decreased towards the critical value. This confirms that using an observation time such that T Re = 10 4 is reasonable to flag stable/unstable computations. Figure 5 shows the critical time-step values for the various time-stepping schemes and the two meshes. For each case, a dependence of the critical time-step on the inverse of the Reynolds number is observed. For convenience, the results are also reported in Tables 2 and3 for the Cartesian and the polygonal meshes, respectively, including a quantitative comparison between the firstand the second-order schemes. Generally speaking, the AC method turns out to be as stable as the monolithic scheme, and no significant differences on the critical time-step are observed. Moreover, the first-order schemes are more stable than the second-order ones, and allow one to choose time-step values more than two times larger than those possible with second-order schemes. Finally, the use of a polygonal mesh only marginally affects the value of the critical time-step.

Re ≈ 33 Re = 100 η Re 10Re 100Re (MONO) Re 10Re 100Re (MONO) 2 (L 2 ) velocity 1.05e-2 1.41e-3 1.24e-3 (1.24e-3) 2.91e-2 3.92e-3 1.38e-3 (1.14e-3) 2 (H 1 ) velocity 1.26e-2 1.04e-2 1.03e-2 (1.03e-2) 2.01e-2 8.25e-3 7.83e-3 (7.80e-3) 2 (L 2 ) pressure 2.81e-2 1.19e-2 1.03e-2 (1.00e-2) 6.02e-2 1.69e-2 1.25e-2 (1.19e-2)

3D modified Taylor-Green vortex

We consider the 3D modified Taylor-Green vortex solution of [19, Benchmark case 2.2] and adapt it so that it is a solution to the Navier-Stokes equations. In particular, we make the solution time-dependent by considering a sinusoidal amplitude. The exact solution reads as follows:

                             u(x,
and the resulting source term is such that

                   f (x, y, z) := α(t)f (x, y, z) + 8 π cos(8 πt)u (x, y, z)+ - α 2 (t) 4     -2 sin(4πx)(cos(4πy) + cos(4πz) -2)
sin(4πy)(cos(4πx) -2 cos(4πz) + 1) sin(4πz)(cos(4πx) -2 cos(4πy) + 1)

    , f (x, y, z) := [-36π 2 cos(2πx) sin(2πy) sin(2πz), 0, 0] T . ( 41
)
The domain is the unit cube Ω := (0, 1) 3 and the viscosity is ν := 1. Hence, since it is again reasonable to set L := 1 and U := 1 for the reference length and velocity, respectively, we obtain a Reynolds number of Re = 1. The observation time is set to T := 2 and the considered time-step values are ∆t ∈ { T 32 , T 64 , T 128 , T 256 }. The convection term is always treated explicitly. A quite refined mesh composed of prisms with polygonal bases is considered. Namely, the mesh obtained from the PrG sequence proposed in the benchmark [START_REF]Finite Volumes for Complex Applications VI -Problems & Perspectives[END_REF] by gluing six of the finest meshes side by side and then rescaling the resulting mesh to the desired dimension (an example is shown in Fig. 6). The resulting mesh has more than 14M cells and leads to a final coupled system of more than 170M DoFs (after static condensation).

Only iterative solvers are considered in this 3D test case. In particular, the linear systems obtained with the AC method are solved by means of a Jacobi-preconditioned Conjugate Gradient (CG) solver. In order to deal with the saddle point systems obtained with the monolithic approach, an iterative Golub-Kahan Bidiagonalization (GKB) [START_REF] Arioli | Generalized Golub-Kahan bidiagonalizationand stopping criteria[END_REF][START_REF] Arioli | An iterative generalized Golub-Kahan algorithm for problems in structural mechanics[END_REF] procedure without augmentation is chosen; the internal solver is a Conjugate Gradient with a K-cycle Algebraic Multigrdid preconditioner [START_REF] Notay | Recursive Krylov-based multigrid cycles[END_REF]. All these solvers are natively available in Code_Saturne. The computations have been run on the GAIA cluster of EDF1 on up to 700 cores. Two values of the AC parameter have been considered: η = 10Re and η = 50Re. Fig. 7 reports the convergence results for the normalized space-time velocity and pressure errors as defined in [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF]. These results corroborate the optimal convergence in time for both first-and second-order schemes. A slight decrease of the convergence rate for the velocity 2 (H 1 )error is observed for the finest ∆t and the second-order schemes, which indicates that the space discretization error tends to become dominant. Moreover, we remark that the higher the value of η, the more accurate the results obtained with the AC method. Actually, for η = 50Re, the errors are essentially superimposed to those obtained with the monolithic approach.

Our aim is now to compare the computational efficiency of the monolithic and AC methods by studying the achieved velocity and pressure errors as a function of the computational cost (defined as the product between the elapsed time and the number of computing cores). As a preliminary step, computations are run in order to find reasonably optimized values for the tolerances of the iterative solvers. The results for the first-order AC and monolithic approaches are, respectively, gathered in Tables 4 and5. A tolerance of 10 -4 appears to be a reasonable choice for the CG iterative solver in the AC context, for both considered values of η. In order to anticipate the greater accuracy of second-order schemes, the tolerance is decreased to 10 -5 for the second-order AC time-stepping scheme. Similarly, for the monolithic approach, the tolerances for the GKB and CG iterative solvers are both set to 10 -4 and 10 -5 , respectively, for the first-and second-order time-stepping schemes.

We can now proceed to the assessment of the computational efficiency of the various schemes with the tolerances prescribed as above. The performances of the three strategies (monolithic, AC(η = 10), and AC(η = 50)), for both first-and second-order schemes, are compared in Fig. 8. The AC method turns out to be more efficient: having fixed the order in time for the three compared methods, the AC(η = 10) and AC(η = 50) methods can achieve a given error threshold (y-axis) in less computation time (x-axis) than the monolithic method. The performance of AC depends on the parameter η: the higher η, the more numerical effort (hence, Velocity -2 (H 1 ) 1.27e-3 5.17e-4 2.4 7.73e+6 1.71e+7 0.5 computational time) is required, and altogether the choice η = 10 turns out to be more effective than the choice η = 50. When considering the second-order schemes, the difference between the AC and monolithic approaches is somewhat mitigated. One explanation is that the bootstrap procedure employed in the second-order AC requires two linear system resolutions per time step, whereas only one is needed with the BDF2 and monolithic coupling. Finally, the results show an advantage of the second-order schemes over the first-order ones, especially when accurate calculations are considered. A comparison at the fixed time-step ∆t = T 128 (corresponding to the finest time step before the space discretization errors tend to become influential, see the discussion of Fig. 7) is presented in Table 6. Focusing first on the bottom right part of the table related to the computational times, one can compare the performances of monolithic and AC methods, with the latter being consistently two times faster than the former, while still providing satisfactory errors. Moreover, the errors obtained with second-order schemes are more than two times smaller than those obtained with first-order schemes, while the computation time are less than two times larger.

Conclusions

We have extended the CDO-Fb schemes to the unsteady incompressible Navier-Stokes problem and investigated the accuracy and efficiency of the artificial compressibility (AC) technique for the time discretization using either a first-order time-stepping scheme or a second-order one. The assessment of the AC technique is made by means of systematic comparisons with the standard monolithic approach, using either a first-order or a second-order version. First, optimal convergence rates in time have been recovered for all the time-stepping schemes in various norms for the velocity and the pressure. Moreover, the stability of the time-stepping schemes when the nonlinear convection term is treated explicitly has been investigated on Cartesian and polygonal meshes, showing in both cases a linear dependency on the reciprocal of the Reynolds number, with only a slightly tighter restriction when using polygonal meshes. Finally, an assessment of the performance of all the schemes has been carried out using large 3D polytopal meshes. The AC method proved to be an efficient alternative to the classical monolithic approach: it ensures accuracy levels close to those obtained with the monolithic approach while using only half of the computational time. Moreover, a comparison between first-and second-order time-stepping schemes indicates an advantage for the latter: second-order schemes provide errors as much as four times smaller than those obtained with first-order schemes, while taking only up to twice the computation time.

As an outlook, the improvement of the preconditioning for large linear systems arising from the AC technique is paramount to achieve scalable and higher efficiency since the presence of the ∇ ∇• term prevents the usual preconditioners such as algebraic multigrid techniques to perform well. The construction of a specific preconditioner adapted from the seminal work in [START_REF] Arnold | Preconditioning in H(div) and applications[END_REF] along with the more recent ones in [START_REF] Benzi | Modified augmented Lagrangian preconditioners for the incompressible Navier-Stokes equations[END_REF] and [START_REF] Farrell | An Augmented Lagrangian preconditioner for the 3D stationary incompressible Navier-Stokes equations at high Reynolds number[END_REF] will be the subject of further investigations.

Figure 2 :

 2 Figure 2: 2D Taylor-Green vortex. Convergence results for the space-time velocity and pressure errors. Re = 1, T = 1.2. Top: first-order schemes; bottom: second-order schemes. Cartesian mesh composed of 512 2 cells.

Figure 3 :

 3 Figure 3: 2D polygonal mesh used for CFL study.

Figure 4 :Figure 5 :

 45 Figure 4: 2D Taylor-Green vortex. Monolithic approach and explicit convection. Divergence time T d (scaled by Re), for different choices of ∆t (scaled by Re). Left: Cartesian mesh composed of 16,384 cells. Right: polygonal mesh composed of 15,129 cells.

Figure 6 :

 6 Figure 6: Coarsest mesh of the PrG sequence of [41] composed of prisms with polygonal basis.

- 2

 2 y, z) := α(t)u (x, y, z) , p(x, y, z) := α(t)p (x, y, z) , α(t) := sin(8πt) , u (x, y, z) cos(2πx) sin(2πy) sin(2πz) sin(2πx) cos(2πy) sin(2πz) sin(2πx) sin(2πy) cos(2πz) , y, z) := -6π sin(2πx) sin(2πy) sin(2πz) ,

Figure 7 :

 7 Figure 7: 3D modified Taylor-Green vortex. Convergence results for the space-time velocity and pressure errors. Re = 1.

Pressure - 2 (L 2 )Figure 8 :

 228 Figure 8: 3D modified Taylor-Green vortex. Comparison of computational efficiency for the various first-and second-order schemes. Legend: see Fig. 7

Table 2 :

 2 2D Taylor-Green vortex. Critical time-step for stability, ∆t s , as a function of the Reynolds number, Re. Comparison of first-and second-order schemes. Cartesian mesh composed of 16,384 cells.

			Monolithic			AC(η = 10Re)
	Re	1 st	2 nd	1 st 2 nd	1 st	2 nd	1 st 2 nd
	200 2.98e-2 1.15e-2 2.60 2.98e-2 1.14e-2 2.61
	500 1.03e-2 4.16e-3 2.48 1.04e-2 4.16e-3 2.51
	700 7.27e-3 2.97e-3 2.45 7.25e-3 2.95e-3 2.46

Table 3 :

 3 2D Taylor-Green vortex. Critical time-step for stability, ∆t s , as a function of the Reynolds number, Re. Comparison of first-and second-order schemes. Polygonal mesh composed of 15,129 cells.

			Monolithic			AC(η = 10Re)
	Re	1 st	2 nd	1 st 2 nd	1 st	2 nd	1 st 2 nd
	200 2.71e-2 1.23e-2 2.21 2.71e-2 1.23e-2 2.21
	500 1.06e-2 4.55e-3 2.33 1.06e-2 4.55e-3 2.33
	700 7.07e-3 3.00e-3 2.36 7.10e-3 3.00e-3 2.37

Table 4 :

 4 3D modified Taylor-Green vortex. Space-time errors for different tolerances of the CG linear solver. First-order AC time-stepping scheme, ∆t = T 128 . The adopted tolerances for further calculations are circled.

	(CG)	2 (L 2 ) u	2 (H 1 ) u	2 (L 2 ) p
			η = 10
	1e-2	6.58e-2 6.76e-4 3.17e-3
	1e-4	5.20e-2 4.91e-4 3.00e-3
	1e-6	5.20e-2 4.91e-4 3.00e-3
			η = 50
	1e-2	3.83e-2 4.22e-4 1.49e-3
	1e-4	3.00e-2 2.92e-4 1.27e-3
	1e-6	3.00e-2 2.92e-4 1.27e-3

Table 5 :

 5 3D modified Taylor-Green vortex. Space-time errors for different tolerances of the GKB and CG linear solvers. First-order monolithic time-stepping scheme, ∆t = T 128 . The adopted tolerances for further calculations are circled.

	(GKB)	(CG)	2 (L 2 ) u	2 (H 1 ) u	2 (L 2 ) p
		1e-2	4.27e-2 8.31e-4 2.80e-3
	1e-2	1e-4	2.84e-2 3.21e-4 1.28e-3
		1e-6	2.83e-2 3.20e-4 1.28e-3
	1e-4	1e-4 1e-6	2.86e-2 2.80e-4 1.13e-3 2.87e-2 2.80e-4 1.13e-3
	1e-6	1e-4 1e-6	2.86e-2 2.80e-4 1.13e-3 2.87e-2 2.80e-4 1.13e-3

Table 6 :

 6 3D modified Taylor-Green vortex. Comparison of accuracy and performance for the various first-and second-order schemes. Time step set to ∆t = T

	128

243 rd of the TOP500 list of November

Acknowledgments

The PhD fellowship of R. Milani was partially supported by EDF R&D and ANRT.


Yara Yaacoub

F Nouvel

S Haese

J.-Y Baudais

VLC Modelization for VLC-PLC System: Evaluation of Optical, Electrical, and Frequency Behavior

Keywords: channel, LED frequency response, optical power, photo-current

The integration of powerline communication (PLC) with visible light communication (VLC) has attracted the attention of many researchers. Cascading these two techniques offers the possibility of performing a wired transmission over existing powerline infrastructure and a last-mile energy-efficient wireless transmission through pre-installed LEDs for indoor illumination. Therefore, VLC should be deeply analyzed to optimize its integration with other systems. So far, most research has only modeled the channel and modulation bandwidth of the VLC system. However, the electrical and optical aspects are rarely taken into account and there is no explicit expression of the received signal as a function of that emitted. So in this paper, An explicit expression is proposed allowing a direct calculation of the photo-current as a function of the transmitted current taking into account the optical, electrical, and frequency behavior of the VLC system. This expression is validated using measurement results.

I. INTRODUCTION

VLC is a relatively new technology based on LEDs. It can transmit on a very large regulation-free bandwidth using intensity modulation and direct detection. Indeed, it is an energy efficient communication system because of the LED capability to support simultaneously communication and lighting [START_REF] Bao | Protocol design and capacity analysis in hybrid network of visible light communication and OFDMA systems[END_REF]. Moreover, VLC is preferred to be used in communicationsensitive indoors like hospitals, as it satisfies electromagnetic compatibility requirements [START_REF] Ding | A hybrid power line and visible light communication system for indoor hospital applications[END_REF]. However, combining VLC with another technology is necessary to increase the network coverage. One of the most suitable technologies is the powerline communication (PLC) because it serves as a natural backbone for the high power LEDs used for illumination [START_REF] Nlom | A simplistic channel model for cascaded PLC-VLC systems[END_REF].

The first hybrid PLC and VLC system prototype was proposed in [START_REF] Komine | Integrated system of white LED visible-light communication and powerline communication[END_REF]. It was based on using a single carrier binary phase-shift keying modulation to provide a low rate transmission. Further works studied OFDM to compensate for the effect of multipath fading [START_REF] Komine | Performance evaluation of narrowband OFDM on integrated system of power line communication and visible light wireless communication[END_REF]. Later, many other scenarios were built to optimize the hybrid system in terms of complexity, cost, rate, etc. [START_REF] Wang | Demonstration of 575-Mb/s downlink and 225-Mb/s uplink bi-directional SCM-WDM visible light communication using RGB LED and phosphor-based LED[END_REF]- [START_REF] Mapfumo | Performance comparison of FFT and DHT hybrid PLC-VLC systems based on DCO-OFDM in impulsive noise[END_REF].

In hybrid communication, there are two types of relaying systems: i) the decode and forward (DF), ii) the amplify and forward (AF). The DF technique consists on decoding the PLC signal and re-encoding it before being transmitted by the VLC channel. Whereas, in the AF the PLC signal is amplified and re-transmitted as it is by the optical system. This technique is less complicated but, at the same time, it can create a big challenge [START_REF] Mossaad | Amplifyand-forward integration of power line and visible light communications[END_REF]. More precisely, when VLC is cascaded to a highly attenuated system like PLC, the available modulation bandwidth of the whole system is reduced. This is owing to the limited frequency response of the LED (few MHz) which can lead to a serious communication failure [START_REF] Ding | A hybrid power line and visible light communication system for indoor hospital applications[END_REF]. Therefore when adopting AF technique, the parameters of the VLC system should be well adjusted and its components should be carefully chosen to ensure successful PLC-VLC communication. However, the LED modulation bandwidth varies according to the technology adopted for generating a white light. Actually, in lighting industry, there are two commonly used techniques: the first one uses a blue LED with a yellow phosphor coat. This coat has a very slow response that reduces the modulation bandwidth up to 1 MHz. The second technique is based on combining three colors: red, green, and blue (RGB LEDs). The RGB LEDs provide a larger bandwidth than that of the phosphor LED. However, balancing their intensities to get the desired white color can be complicated in manufacturing [START_REF] Ding | A hybrid power line and visible light communication system for indoor hospital applications[END_REF].

Unlike VLC, PLC technology has already been deployed and many standards have been elaborated like HomePlug 1.0.1 (2005), HomePlugAV (2005)(2006)(2007), HD-PLC (2007), IEEE 1901 (2010), G3-PLC (2010) or HomePlugAV2 (2010) [START_REF] Laere | Powerline communication transmission performance study in field conditions: Cases of G3-PLC for smart metering and Homeplug-AV2 for railway signaling[END_REF].

Our objective is to be able to transmit HomePlugAV-like signals through a trichromatic lighting system (commercial RGB LEDs) using the AF method. But as explained before, when adopting AF, the already attenuated PLC signal is remarkably affected by the VLC system. So, a preliminary study of the VLC system should be made. Unfortunately, there is no explicit expression making it possible to evaluate the received VLC signal as a function of the transmitted signal taking into account all of the optical, electrical and frequency behavior of the system. In fact, most research has only modeled the VLC channel [START_REF] Ding | A hybrid power line and visible light communication system for indoor hospital applications[END_REF], [START_REF] Nlom | Cascaded plc-vlc channel: An indoor measurements campaign[END_REF]. In [START_REF] Mossaad | Visible light communications using OFDM and multiple LEDs[END_REF] the optical and electrical nonlinearity of LEDs is investigated. However, the spectral nonlinearity of the photo-current generated at the receiver is not taken into account. Further, in [START_REF] Chun | Led based wavelength division multiplexed 10 gb/s visible light communications[END_REF] the optical side of the VLC system is covered without considering the effect of the LED optical nonlinearity. In [START_REF] Zhou | Visible light communication system evaluations with integrated hardware and optical parameters[END_REF], a theoretical analysis of the relationship between the current energy of the hardware circuit and the optical power is performed taking into account all nonlinearties (optical, electrical and frequency) of the VLC system. Due to the complexity of this model, an explicit expression to simplify the calculation of the received signal is not obtained. However, This expression is very important for a quick evaluation of the hardware system. Thus, in our paper, an attempt is made to compensate for these shortcomings. Consequently, an explicit expression permitting the calculation of the photo-current as a function of the transmitted signal is proposed. This expression takes into account the electrical, optical and frequency behavior of the VLC system in order to built and optimize PLC-VLC system.

The remainder of this paper is organized as follow: In Section II, VLC channel is discussed and the LED frequency response expression is also figured out based on experimental results. In Section III-A, an analytical expression is provided to describe the relationship between the maximum transmitted and generated photo-current. Measurements are carried out in Section III-B to validate the accuracy of this expression. In Section III-C, a comparison was made between the measurement results and the values obtained using our model with and without taking into account the effect of the optical nonlinearity and the LED spectral width. Finally, conclusions are given in Section IV.

II. VLC CHANNEL MODEL

The VLC channel is actually composed from a direct line of sight (LOS) and a non direct line of sight (NLOS). A NLOS has a very low power compared to the LOS. According to [START_REF] Komine | Fundamental analysis for visible-light communication system using LED lights[END_REF], the NLOS can be neglect and the channel model considers only the LOS. Hence, the channel is represented by the following expression:

H C = (m+1)A cos m (Φ) cos(ψ) 2πD 2 G, ψ < ψ FOV 0 ψ > ψ FOV ( 1
)
where m is the Lambertian index, A is the active area of the photodiode, Φ and ψ are the angles of irradiance and incidence respectively, D is the distance between the LED and the photodiode, ψ FOV is the irradiance angle of the field of vision of the receiver and G is the gain of the lens that concentrates the light to the photodiode. Moreover, the LED limited bandwidth is an important factor that can affect the VLC and the hybrid channel as well. Hence, the overall system attenuation can be expressed as follows:

H V = H L H C (2)
where H L is the LED frequency response (LFR).

In order to model the attenuation of the VLC system, the LFR should be measured and then modeled. Hence, a wideband VLC transmitter was manufactured. The latter consists of monochromatic LED (red, green or blue) from the RGB module, LUXEON MultiColor LED Module 2.5 W [START_REF]Luxeon multicolor module 2.5W, DS244[END_REF], and a wide-band bias tee circuit up to 100 MHz at -3 dB. Figure 1a illustrates the realized bias tee which is a three-port network used to set the DC bias point of certain electronic components without disturbing other components. The equivalent circuit model is represented in Figure 1b where the low-frequency port is used to define the polarization; the high-frequency port passes the data signals but blocks bias levels; the combo port is connected to the device, which can provide both polarization and data signals [START_REF] Stratil | Broadband over visible light: High power wideband bias-t solution[END_REF] (see Figure 1b). Back to LFR measurements, in our study, each LED are biased at 100 mA. As shown in Figure 2, the red LED has the widest modulation band of 20 MHz at -3 dB. However, the green and blue LEDs have approximately the same response with a -3 dB cutoff frequency of 5 MHz. Hence, the LFR can be modeled by a first order low-pass filter:

H L (f ) = 1 1 + j f fc (3)
where f c is the -3 dB cutoff frequency of the LED LFR. Therefore, the red LED is chosen to perform the first PLC-VLC transmission tests, and thus it is used to do the expression validation measurement.

III. ANALYTIC EXPRESSION OF THE TRANSMITTED CURRENT

A. Proposed model

When dimensioning and optimizing a PLC-VLC system, it is important to have an estimate of the maximum current amplitude that the photodiode would generate when transmitting signals through the LED, without the need to repeat the measurements each time the amplitude of the transmitted signal changes. Thus, an analytical relationship between the transmitted and the generated currents must be found. In this section, our proposal is described and its accuracy is validated by measurement results. In LED datasheet, the typical transmitted power density at the LED's dominant wavelength P is not given explicitly. Although, all the curves presented in the LED data sheet are normalized to this value. However, it can be calculated using the luminous flux expression [START_REF] Jabban | Theoretical and experimental optimization of DMT-based visible light communication under lighting constraints[END_REF] as following:

P = φ V max 780 380 v(λ)p(λ) dλ (4)
where φ is LED the luminous flux, v(λ) is normalized to its maximum eye efficacy V max (V max = 683 lm/W at λ = 555 nm) and p(λ) is the optical power density normalized to P . φ and p(λ) are given in the LED datasheet. To be able to calculate the integral, p(λ) and v(λ) are modeled by the sum of two Gaussian functions [START_REF] Zhou | Visible light communication system evaluations with integrated hardware and optical parameters[END_REF]:

p(λ) = ae -( λ-b c ) 2 + de -( λ-e f ) 2 (5)
and

v(λ) = a 1 e -( λ-b 1 c 1 ) 2 + d 1 e -( λ-e 1 f 1 ) 2 (6)
Moving to the relationship between the light output and the current crossing the LED L(i), it is given in the datasheet normalized to its typical and it could be expressed by a polynomial function as following:

L(i) = k 1 i 3 + k 2 i 2 + k 3 i + k 4 (7)
So, the LED transmitted power density in function of wavelength at a determined current can be described by:

P tr (λ, i) = P p(λ)L(i) (8)
However, this transmitted power is affected by the channel attenuation and the LFR. Hence, the received power density is written as follows:

p rec (λ, i) = P tr H V (9)
The photodiode receives the light and generates the photocurrent depending on its sensitivity. The following equation shows the current spectral density as a function of the received optical power density:

i P (λ, i) = Ss(λ)p rec (λ, i) ( 10
)
where S is the maximum sensitivity of the photodiode, s(λ) is given in the photodiode's datasheet. In order to reduce the complexity of calculation in the following steps, we modeled only the part of the spectrum where the RGB LEDs radiate, as follows:

s(λ) = s 1 λ 2 + s 2 λ + s 3 (11)
Hence, the overall generated current is the integral of the current density over the all spectrum:

I(i) = 780 380 i P (λ, i)dλ ( 12
)
where the integration boundries are the extremums of the visible spectrum. However, since the boundaries of p(λ) are too attenuated for wavelengths greater than 780 nm and less than 380 nm. the integral of I(i) can be considered from 0 to +∞ to simplify the final expression The final expression of the generated current is:

I(i) = P H B H V SL(i) +∞ 0 s(λ)p(λ)dλ = P H B H V SL(i) √ π ca s 1 ( c 2 2 + b 2 ) + s 2 b + s 3 +df s 1 ( f 2 2 + e 2 ) + s 2 e + s 3 ) (13)
where H B is the bias tee frequency attenuation. P , H V and L(i) are respectively the typical power, the VLC channel, and the light output described in (4), [START_REF] Ding | A hybrid power line and visible light communication system for indoor hospital applications[END_REF], and [START_REF] Mossaad | Amplifyand-forward integration of power line and visible light communications[END_REF]. The expression results has been compared to measurements performed on our testbed described in Section II.

B. Current measurement setup

In this part we validate (13) using the red LED of the previously described transmitter in Section II. The LED is biased at 110 mA, the receiver used involves the BPW34 photodiode [START_REF]Silicon PIN photodiode[END_REF] and a trans-impedance amplifier (TIA) of 50 kΩ. The transmitter and the receiver are aligned and separated by a distance of 13 cm. The 50 kHz transmitted current has a maximum amplitude of 10 mA. The bias tee attenuation at this frequency is 0.6. The measurement result shows that the maximum voltage amplitude at the output of TIA amplifier is 10 mV which means that the maximum amplitude of the generated photo-current is 200 nA.

C. Results analysis

As Figs. 3a to 3c show, the proposed model of p(λ) (5), L(i) [START_REF] Mossaad | Amplifyand-forward integration of power line and visible light communications[END_REF] and s(λ) [START_REF] Mossaad | Visible light communications using OFDM and multiple LEDs[END_REF] fits properly the curves of the LED and photodiode datasheets. In order to show the improvement brought by this model when it takes into account the optical nonlinearity L(i) and the LED spectral power width p(λ) with the spectral sensitivity of the photodiode s(λ). The maximum generated photocurrent I and the TIA's output voltage V T IA were calculated based on the test-bed parameters described in Section III-B at three different cases: i) when all nonlinearities are considered (our proposed model), ii) when the light output L(i) is considered linear, iii) when the photodiode sensitivity s(λ) and p(λ) are considered only at the dominant wavelength (λ dominant = 625 nm in the case of the LED red) . The results in Table I show the precision enhancement in calculation that our model provides in comparison to the results obtained in the last two cases. The calculated V T IA is compared to the measurement (10 mV), resulting in the relative error ERR. It has decreased from 41 % in case ii and 27 % in case iii to 22 % in case i. The 22 % error obtained in our model is acceptable regarding the test-bed characteristics. In fact, during measurements, the junction temperature cannot be controlled despite the great influence of the junction temperature on LED characteristics. Especially on the red LED which seems the most affected by the increase in junction temperature [START_REF]Luxeon multicolor module 2.5W, DS244[END_REF]. In fact, high power LEDs used for illumination run at high drive current leading to high junction temperature which cause signal power degradation over the time [START_REF] Thomas | Experimental study on the effect of junction temperature on power led's[END_REF]. However, the thermal effect is not treated in this paper. Thus, this simple expression provides a good estimate of the maximum current amplitude that a photodiode could generate, which can help when dimensioning the PLC-VLC system to maintain the linearity of the receiver and not end up being saturated.

IV. CONCLUSION

In this paper, the optical, electrical, and frequency aspects of the VLC system are studied. The optical channel model is detailed by taking into consideration the LFR. The relationship between the transmitted and the generated current is described by an explicit expression. The efficacy of this expression is tested by a VLC system manufactured in our lab using the red LED of the RGB module. The precision improvement of the model is highlighted by comparing calculation results with and without considering the optical nonlinearity and the LED spectral width. Therefore, a VLC system model is proposed replicating as much as possible its actual optical, electrical, and frequency behavior. A simple expression is established allowing a fast calculation of the photo-current. This expression facilitates the dimensioning and the optimization of the VLC system, especially when it is going to be integrated with other systems like PLC.

Figure 1 :

 1 Figure 1: Bias Tee. (a) Our manufactured circuit, (b) explanatory scheme.

Figure 2 :

 2 Figure 2: Red, green, and blue LED frequency responses.

Figure 3 :

 3 Figure 3: LED and photodiode datasheet curves models. (a) The spectral power density of red, green, and blue LEDs compared to their modeled curves, (b) the three LEDs' electrical/optical characteristics compared to our polynomial model, (d) spectral sensitivity response of the photodiode compared to our polynomial model.

Table I :

 I comparison between V T IA calculated in 3 different cases: considering all the non-linearities, considering L(i) linear, and considering s(λ) and p(λ) at the dominant wavelength of the LED .

	Assumption	I(nA)	V T IA (mV )	ERR(%)
	i) proposed model	244	12.2	22
	ii) L(i) linear	281	14.1	41
	iii) p(λ) and s(λ) at λ dominant	254	12.7	27

ACKNOWLEDGMENT

The authors would like to thank the Direction générale de l'armement maîtrise de l'information (DGA-MI) for its general supervision and financial support.