text:
G Jeanmasson
email: guillaume.jeanmasson@onera.fr

I Mary

L Mieussens

On some explicit local time stepping finite volume schemes for CFD

Keywords: Explicit local time stepping schemes, Runge-Kutta methods, LES, turbulent flows

In this paper, two local time stepping schemes of order two and three in time are proposed. By construction, they are not mass conservative but a correction stage is added to make them conservative. These schemes are compared with some local time stepping schemes existing in the literature (schemes of Constantinescu and Sandu). The comparisons are carried out on various test cases. They prove that our schemes are efficient and our third order local time stepping has a higher time accuracy than the schemes based on the strategy of Constantinescu and Sandu. Our third order local time stepping scheme is used to perform an industrial-like test case: a 3D Large-Eddy Simulation over an airfoil.

Introduction

Numerical simulations are widely used to study fluid mechanics problems. Simulations are based on the discretized form of fluid mechanics equations (mainly Euler or Navier-Stokes equations). This discretized form is obtained by applying space and time integration methods. Finite volume approach is one of the most popular space integration method in CFD. For time integration, explicit Runge-Kutta methods have proven their efficiency for solving fluid mechanics problems. However, with these methods, the time step is strongly restricted by the numerical stability condition (CFL condition).This stability condition is defined in each cell. For each cell, it imposes that the time step is lower than the ratio between the grid size and the sum of the sound speed and the fluid velocity in the cell. The time step for the entire domain is chosen to be lower than the smallest 1 ratio in the mesh (generally imposed by the smallest cell in the mesh in our applications, since the fluid velocity is rather uniform). Consequently, the time step may be much smaller than necessary in the other cells, which is not efficient. It is particularly true in most CFD computations involving walls, which present some zones of mesh refinements where the stability condition is much more restrictive than in the other areas.

Implicit, unconditionally stable timestepping methods allow to use larger time steps, which reduces the computational time as compared to explicit methods. Note that the time step can not be too large in order to keep a good accuracy of the solution. For a Large Eddy Simulation of turbulent flows, the ratio between the time step and each cell characteristic length can not exceed 2-3 [START_REF] Choi | Effects of the computational time step on numerical solutions of turbulent flow[END_REF]. Moreover, implicit methods seem to have two drawbacks. First, an efficient parallel implementation of these methods may be difficult, especially if a strong scalability is required. In addition, implicit methods can be less accurate than explicit methods due to a residual numerical noise.

The use of an explicit local time stepping scheme seems to be a good approach to reduce the computational cost of explicit methods keeping their advantages with respect to implicit methods (solution accuracy and good compatibility with parallel programming). Indeed, an explicit local time stepping scheme allows a time step adaptation on the mesh to satisfy local stability conditions. Small time steps are used on cells which present the most severe CFL and larger time steps are applied on cells where CFL is less restrictive. This should result in a reduction of the computational cost with respect to explicit methods using a uniform time step.

Two approaches are possible to construct explicit local time stepping schemes. In the first approach, each cell of the computational domain is integrated with its own maximal allowable time step [START_REF] Unfer | An asynchronous framework for the simulation of the plasma/flow interaction[END_REF][START_REF] Toumi | Asynchronous numerical scheme for modeling hyperbolic systems[END_REF][START_REF] Semiletov | Cabaret scheme for computational aero acoustics: Extension to asynchronous time stepping and 3d flow modelling[END_REF]. In the second approach, each time step of the computational domain is defined as a fraction of the largest time step and each time step is associated to a group of contiguous cells [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF][START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF][START_REF] Tang | A class of high resolution schemes for hyperbolic conservation laws and convection-diffusion equations with varying time and space grids[END_REF][START_REF] Schlegel | Multirate runge-kutta schemes for advection equations[END_REF][START_REF] Liu | Nonuniform time-step runge-kutta discontinuous galerkin method for computational aeroacoustics[END_REF][START_REF] Liu | Nonuniform-time-step explicit runge-kutta scheme for high-order finite difference method[END_REF][START_REF] Seny | Multirate time stepping for accelerating explicit discontinuous galerkin computations with application to geophysical flows[END_REF][START_REF] Savcenco | A multirate time stepping strategy for stiff ordinary differential equations[END_REF][START_REF] Dawson | High resolution schemes for conservation laws with locally varying time steps[END_REF][START_REF] Kvaerno | Low order multirate runge-kutta methods in electric circuit simulation[END_REF][START_REF] Sexton | [END_REF][START_REF] Sandu | A class of multirate infinitesimal gark methods[END_REF]. In the litterature, these schemes are frequently called "multirate schemes". All the local time stepping schemes presented in this article are constructed by following this second approach. Indeed, this approach seems well suited to High Performance Computations (HPC) since the memory accesses are much more efficient due to the regular predictable pattern of the cells groups. However, local time stepping schemes seem to have an important drawback: it is difficult to guarantee both high order time accuracy and mass conservation. In the 1980's, Osher and Sanders [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF] proposed a locally varying time stepping scheme based on the forward Euler method. It is conservative but only first order accurate in time. Dawson and Kirby [START_REF] Dawson | High resolution schemes for conservation laws with locally varying time steps[END_REF] extended the work of Osher and Sanders by performing a conservative local time stepping scheme based on a second order Runge-Kutta scheme. Tang and Warnecke [START_REF] Tang | A class of high resolution schemes for hyperbolic conservation laws and convection-diffusion equations with varying time and space grids[END_REF] proved that the schemes in [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF][START_REF] Dawson | High resolution schemes for conservation laws with locally varying time steps[END_REF] have a lack of consistency and proposed a new second order time accurate local time stepping scheme. However, their scheme is not mass-conservative. The local time stepping schemes presented in [START_REF] Hundsdorfer | Analysis of explicit multirate and partitioned Runge-Kutta schemes for conservation laws[END_REF][START_REF] Liu | Nonuniform time-step runge-kutta discontinuous galerkin method for computational aeroacoustics[END_REF][START_REF] Liu | Nonuniform-time-step explicit runge-kutta scheme for high-order finite difference method[END_REF] are also second order time accurate but not mass-conservative. In [START_REF] Schlegel | Multirate runge-kutta schemes for advection equations[END_REF], a comparison is made between different multirate Runge-Kutta schemes and a multirate time stepping scheme of order three in time is proposed. This third order local time stepping scheme borrows some ideas of implicit-explicit (IMEX) schemes [START_REF] Knoth | Implicit-explicit runge-kutta methods for computing atmospheric reactive flows[END_REF]. However, it is not mass-conservative. Another third order non conservative explicit local time stepping scheme based on the theory of partionned Runge-Kutta methods (PRK) is proposed in [START_REF] Kvaerno | Low order multirate runge-kutta methods in electric circuit simulation[END_REF]. Recently, multirate Runge-Kutta schemes up to order four are developed by several authors [START_REF] Sexton | [END_REF][START_REF] Sandu | A class of multirate infinitesimal gark methods[END_REF]. They are derived from the Multirate Infinitesimal Step method, introduced in [START_REF] Wensch | Multirate infinitesimal step methods for atmospheric flow simulation[END_REF]. These methods are very interesting since they can deal with several time step rates ("telescopic approach"), unlike [START_REF] Wensch | Multirate infinitesimal step methods for atmospheric flow simulation[END_REF]. The authors also provide rigorous stability analyses and propose numerical applications to several ordinary differential equations that confirm their theoretical findings. However, the authors don't perform numerical test on partial differential equations (PDEs), and their method seems to be not mass-conservative when applied to PDEs.

1 Local time stepping scheme and Partitioned Runge-Kutta methods

Semi-discretization of a one-dimensional hyperbolic equation

We consider the one-dimensional hyperbolic equation :

∂y(t, x) ∂t + ∂f (y(t, x)) ∂x = 0, (1)
with y(0, x) = y 0 (x), x ∈] -∞, +∞[ and t > 0. We also consider a 1D computational domain divided into cells of variable length ∆x i . Cell i is bounded by points x i- . The finite volume method is applied to (1). After space integration over cell i, we obtain :

∂y i ∂t = - 1 ∆x i (F i+ 1 2 -F i-1 2 ), (2)
where

y i = x i+ 1 2 x i-1 2 y(t, x)dx. The numerical flux at point x i+ 1 2 is denoted by F i+ 1 2
, and depends on several values around y i . For instance, we have:

F i+ 1 2
(y i-1 , y i , y i+1 , y i+2 ) in the case of a flux function with a four cell stencil. For the sake of clarity for the reader, F i+ 1 2

(y i-1 , y i , y i+1 , y i+2 ) is denoted by F i+ 1 2
. For the sake of convenience, we introduce D i , the flux balance on cell i :

D i = 1 ∆x i (F i+ 1 2 -F i-1 2 ). (3)
The flux balance D i is obviously a function of y i , and several values around it (y i-2 , y i-1 , y i+1 ...). The arguments of D i are omitted for a better clarity. With this notation, (2) now reads:

∂y i ∂t = -D i . (4)
Now ( 2) and ( 4) can be viewed as ordinary differential systems that can be solved with an explicit Runge-Kutta (RK) method, as described in the following section. This is the so called method of lines (MOL).

Explicit Runge-Kutta methods

We consider a non-autonomous ordinary differential equation:

dy dt = f (t, y) y(t = 0) = y 0 (5)
The solution of (5) at time t n is y n . A s-step explicit Runge-Kutta method allows to compute y n+1 (solution of (5) at time t n+1 = t n + ∆t) by using y n and s -1 intermediate values. In [START_REF] Hairer | Solving ordinary differential equations I : Nonstiff Problems[END_REF], a s-step explicit Runge-Kutta method is defined by the formulas :

t n,i = t n + c i ∆t, (6)
y n,i = y n + ∆t i-1 j=1 a ij K n,j , (7)
K n,i = f (t n,i , y n,i ), (8)
where 1 ≤ i ≤ s, and K n,1 = f (t n,1 , y n,1 ) = f (t n , y n ). The i -th intermediate time, the i -th intermediate value, and the i -th intermediate evaluation of the f function between t n and t n + ∆t are denoted by t n,i , y n,i and K n,i . The value at time t n + ∆t is given by the following formula:

y n+1 = y n + ∆t s i=1 b i K n,i (9)
The method is defined by its coefficients A = {a ij }, b = {b i } and c = {c i } that can be represented in a Butcher tableau, Table 1 (see [START_REF] Hairer | Solving ordinary differential equations I : Nonstiff Problems[END_REF]) :

c A b c 1 = 0 0 c 2 a 21 c 3 a 31 a 32 . . . . . . c s a s1 a s2 . . . a s,s-1 b 1 b 2 . . . b s-1 b s
c i = i-1 j=1 a ij . ( 10
)

Partitioned Runge-Kutta methods

In [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF] and [START_REF] Hundsdorfer | Analysis of explicit multirate and partitioned Runge-Kutta schemes for conservation laws[END_REF], Partitioned Runge-Kutta (PRK) methods are applied to temporal multiscale problems.

Let us consider a 1D computational domain composed of n cells (figure 1). The finite volume method applied to the one-dimensional hyperbolic equation (1) on the whole domain leads to the autonomous system of ordinary differential equations:

dy dt = f (y), y ∈ R n . ( 11
)
As shown in figure 1, the computational domain is divided into two subdomains which have different characteristic times. The first subdomain has a slow characteristic time and contains cells from cell 1 to cell i 0 (large cells). In the following, this subdomain is referred to as the "slow subdomain". The second subdomain has a fast characteristic time and contains cells from cell i 0+1 to cell n (small cells). In the following, this subdomain is referred to as the "fast subdomain". We denote by Y S = [y 1 , y 2 , ..., y i 0 ] the vector of unknowns in the slow subdomain and by Y F = [y i 0+1 , y i 0+2 , ..., y n ] the vector of unknowns in the fast subdomain. System (11) can be rewritten using the variables Y S and Y F :

| x 1 2 y 1 cell 1 | x 3 2 | | | x i 0 -1 2 y i 0 cell i 0 x i 0 + 1 2 | x i 0 + 3 2 y i 0+1 cell i 0+1 | | | | | | | x n-1 2 x n+ 1 2 y n cell n ∆x i0 slow subdomain fast subdomain
     dY S dt = f S (Y S , Y F ) dY F dt = f F (Y S , Y F ) (12)
The dependency between vectors Y S and Y F in each equation of [START_REF] Savcenco | A multirate time stepping strategy for stiff ordinary differential equations[END_REF] occurs at the interface between the subdomains. The dissimilar time scale for Y S and Y F is a motivation for solving these equations with two Runge-Kutta methods that a use a different time step. In the following, we will show that these methods are, in fact, derived from the same base method. This is useful for the generalization of the scheme to more than two time steps. We denote by RKS the s-stage Runge-Kutta method for the slow subdomain. It is used to solve the first equation and its coefficients are : A S , b S , c S . Similarly, we denote by RKF the s-stage Runge-Kutta method for the fast subdomain. It is used to solve the second equation and its coefficients are : A F , b F , c F . According to [START_REF] Hairer | Solving ordinary differential equations I : Nonstiff Problems[END_REF], the PRK method using RKS and RKF is given by the following formulas:

Y n+1 S = Y n S + ∆t s i=1 b S i K n,i S Y n+1 F = Y n F + ∆t s i=1 b F i K n,i F Y n,i S = Y n S + ∆t i-1 j=1 a S ij K n,j S Y n,i F = Y n F + ∆t i-1 j=1 a F ij K n,j F K n,i S = f L (Y n,i S , Y n,i F ) K n,i F = f F (Y n,i S , Y n,i F ) (13)
where

1 ≤ i ≤ s with K n,1 S = f S (Y n S , Y n F ) and K n,1 F = f F (Y n S , Y n F ).
We remind that Y n S and Y n F denote the vectors of unknowns at time t n in the slow subdomain and in the fast subdomain, respectively. Y n+1 S and Y n+1 F are the vectors of unknowns at time t n + ∆t. The i -th intermediate slow solution and the i -th evaluation of the slow function f S between t n and t n + ∆t are denoted by Y n,i S and K n,i S .The i -th intermediate fast solution and the i -th evaluation of the fast function f F between t n and t n + ∆t are denoted by Y n,i F and K n,i F . In most explicit local time stepping schemes (for instance [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF], [START_REF] Tang | A class of high resolution schemes for hyperbolic conservation laws and convection-diffusion equations with varying time and space grids[END_REF], [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF] and [START_REF] Schlegel | Multirate runge-kutta schemes for advection equations[END_REF]), RKS and RKF methods are based on the same RK method. RKF method is composed of two successive applications of the base method with the time step ∆t 2 , and some additional values to make the transition with the RKS method. RKS method is composed of the base method with the time step ∆t, and some additional values to make the transtion with the RKF method. With this strategy, an explicit local time stepping scheme with time steps ∆t in the slow subdomain and ∆t 2 in the fast subdomain is obtained. To have a Partitioned Runge-Kutta method of order p and therefore, a local time stepping scheme of order p, both RKS and RKF methods must be of order p but this is not sufficient. The coefficients of the methods have to satisfy additional coupling conditions. Following the methodology of [START_REF] Hairer | Solving ordinary differential equations I : Nonstiff Problems[END_REF] we can derive the following conditions on the Runge-Kutta coefficients to obtain a PRK method of order two and three, respectively (tables 2 and 3).

RKS of order 2 RKF of order 2 Coupling conditions

s i=1 b S i = 1 s i=1 b F i = 1 s i=1 b S i c F i = 1 2 s i=1 b S i c S i = 1 2 s i=1 b F i c F i = 1 2 s i=1 b F i c S i = 1 2
Table 2: Conditions to obtain a second order accurate PRK method RKS of order 3 RKF of order 3 Coupling conditions

s i=1 b S i = 1 s i=1 b F i = 1 s i=1 b S i c F i = 1 2 s i=1 b F i c S i = 1 2 s i=1 b S i c S i = 1 2 s i=1 b F i c F i = 1 2 s i=1 i-1 j=1 b S i a F ij c S j = 1 6 s i=1 i-1 j=1 b F i a S ij c F j = s i=1 b S i (c S i ) 2 = 1 3 s i=1 b F i (c F i ) 2 = 1 3 s i=1 i-1 j=1 b S i a F ij c F j = 1 6 s i=1 i-1 j=1 b F i a S ij c S j = s i=1 i-1 j=1 b S i a S ij c S j = 1 6 s i=1 i-1 j=1 b F i a F ij c F j = 1 6 s i=1 i-1 j=1 b S i a S ij c F j = 1 6 s i=1 i-1 j=1 b F i a F ij c S j = s i=1 b S i c F i c S i = 1 3 s i=1 b F i c F i c S i = 1 3
Table 3: Conditions to obtain a third order accurate PRK method

The number of conditions quickly increases with the order of accuracy which makes the implementation of high order local time stepping schemes difficult.

Constantinescu and Sandu strategy

In [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF], Constantinescu and Sandu have developed a strategy, based on PRK methods, to obtain conservative and second order accurate local time stepping schemes. In this subsection, we briefly describe two schemes based on this strategy. The first one is the original scheme, based on a RK2 method [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. The second one follows the same strategy but is based on a RK3 method (see [START_REF] Schlegel | Multirate runge-kutta schemes for advection equations[END_REF][START_REF] Seny | Multirate time stepping for accelerating explicit discontinuous galerkin computations with application to geophysical flows[END_REF]). A comparison between these schemes and the schemes we propose (that are based on a different strategy, see part 3) is performed in part 4.

First scheme : original scheme of Constantinescu and Sandu

The original scheme developed by Constantinescu and Sandu is based on the second order Heun method:

0 0 1 1 0 1 2 1 2
Table 4: Base method for the original scheme of Constantinescu and Sandu

The RKF method (Table 6) is obtained by applying twice successively the base method with the time step ∆t 2 . The RKS method (Table 5) is obtained by applying twice successively the base method with the time step ∆t and dividing all b i coefficients by 2. Table 6: RKF method Both RKS and RKF methods are second order accurate and we have b S i = b F i for every i. These properties ensure that this local time stepping scheme is second order accurate [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. The equality of coefficients b i ensures that the scheme is mass conservative [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. In [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF], the authors demonstrate that the RKS method described in Table 5 is only applied in a buffer zone. It is an area of variable length (the length depends on the spatial stencil), adjacent to the interface. Outside this buffer zone, RKS reduces to the Heun method described in Table 4. Finally, our semi-discrete approximation ( 4) is solved by the algorithm shown in table 7. We introduce the following notations for this algorithm :

• for all n, the values at times t n and t n+1 (with t n+1 = t n + ∆t, ∆t corresponding to the largest time step) are denoted by y n and y n+1 ,

• for all n, the k -th intermediate value between times t n and t n+1 is denoted by y n,k .

According to (3) and the notations given above :

• for all n, D n i = 1 ∆x (F i+ 1 2 (y n i , y n i+1 ) -F i-1 2 (y n i-1 , y n i )),
• for all n, the k-th intermediate flux balance on cell i is given by :

D n,k i = 1 ∆x (F i+ 1 2 (y n,k i , y n,k i+1 )- F i-1 2 (y n,k i-1 , y n,k i )).

Slow subdomain

Fast subdomain We can notice that this scheme is internally inconsistent, since at the interface, some values estimated at different times are mixed in the flux balance D. For instance, at the interface, D n,2 i contains values estimated at time t n (y n,2 i from the buffer zone) and values estimated at time t n+ 1 2 (y n,2 i from the fast subdomain). Despite this inconsistency, the scheme is second order accurate and mass-conservative.

1 y n,1 i = y n i -∆tD n i y n,1 i = y n i -∆tD n i y n,1 i = y n i -∆t 2 D n i 2 y n,2 i = y n i y n,2 i = y n i -∆t 4 [D n i + D n,1 i ] 3 y n,3 i = y n i -∆tD n,2 i y n,3 i = y n,2 i -∆t 2 D n,2 i 4 y n+1 i = y n i -∆t 2 [D n i + D n,1 i ] y n+1 i = y n i -∆t 4 [D n i y n+1 i = y n,2 i -∆t 4 [D n,2 i + D n,3 i ] +D n,1 i + D n,2 i + D n,3 i ] Buffer zone
For each cell in the buffer zone, four flux balances have to be calculated, therefore the computational cost is the same as in the fast subdomain. In our applications, we use a flux function with a four cell stencil. In this case, it can be proven that the size of the buffer zone is equal to four cells.

Second scheme

In [START_REF] Schlegel | Multirate runge-kutta schemes for advection equations[END_REF] and [START_REF] Seny | Multirate time stepping for accelerating explicit discontinuous galerkin computations with application to geophysical flows[END_REF], a RK3 method is used as the base method and the strategy of Constantinescu and Sandu is applied. Here, we extend this approach using a RK3 low-storage scheme [START_REF] Williamson | Low-storage runge-kutta schemes[END_REF] as the base method. Low-storage Runge-Kutta schemes require less storage than "classical" Runge-Kutta schemes. The RK3 low-storage method applied to [START_REF] Seny | Multirate time stepping for accelerating explicit discontinuous galerkin computations with application to geophysical flows[END_REF] is shown in Table 8.

Step 1

K n = ∆tf (t n , y n ) y n,1 = y n + β 1 K n β 1 = 1 2 Step 2 K n,1 = α 2 K n + ∆tf (t n + c 2 ∆t, y n,1 ) α 2 = -0.6830127 y n,2 = y n,1 + β 2 K n,1 c 2 = 1 2 β 2 = 0.9106836 Step 3 K n,2 = α 3 K n,1 + ∆tf (t n + c 3 ∆t, y n,2 ) α 3 = -4 3 y n+1 = y n,2 + β 3 K n,2 β 3 = 0.3660254 c 3 = 0.7886751
Table 8: The RK3 low-storage scheme written in "low-storage form"

The values of the coefficients α and β have been proposed by Lowery and Reynolds [START_REF] Lowery | Numerical simulation of a spatially-developping, forced, plane mixing layer[END_REF]. Successive values of K and y overwrite the previous ones so that at any stage only 2N storage locations are required. It is to be noticed that this scheme is not total variation diminishing (TVD). Indeed, in [START_REF] Gottlieb | Total variation diminishing runge-kutta schemes[END_REF] the authors show that a Runge-Kutta method written as a combination of Euler steps is TVD if all coefficients are positive. It is possible to show that this scheme does not satisfy this condition (because of the negative value of α 2 , for instance). The same method can be written in the "classical Runge-Kutta form" and a Butcher tableau can be associated:

Step 1

K n = ∆tf (t n , y n ) y n,1 = y n + a 21 K n a 21 = β 1 Step 2 K n,1 = ∆tf (t n + c 2 ∆t, y n,1 ) a 31 = β 1 + β 2 α 2 y n,2 = y n + a 31 K n + a 32 K n,1 a 32 = β 2 c 2 = 1 2
Step 3

K n,3 = ∆tf (t n + c 3 ∆t, y n,2 ) b 1 = a 31 + β 3 α 2 α 3 y n+1 = y n + b 1 K n + b 2 K n,1 + b 3 K n,2 b 2 = β 2 + β 3 α 3 b 3 = β 3 c 3 = 0.7886751 0 0 c 2 a 21 c 3 a 31 a 32 b 1 b 2 b 3
Table 9: The RK3 low-storage scheme written in "classical form" and its associated Butcher tableau

In the following, we only use the "classical form" described in table 9 for our RK3 scheme. Indeed, its associated Butcher tableau is very useful for constructing and describing local time stepping schemes.

Following the strategy of Constantinescu and Sandu, the base method is duplicated to construct RKS (table 10) and RKF (table 11) methods.

0 0 c 2 a 21 c 3 a 31 a 32 0 0 0 0 c 2 0 0 0 a 21 c 3 0 0 0 a 31 a 32 b 1 2 b 2 2 b 3 2 b 1 2 b 2 2 b 3 2
Table 10: RKS method

0 0 c 2 2 a 21 2 c 3 2 a 31 2
a 32 2 1 2 b 1 2 b 2 2 b 3 2 1+c 2 2 b 1 2 b 2 2 b 3 2 a 21 2 1+c 3 2 b 1 2 b 2 2 b 3 2 a 31 2 a 32 2 b 1 2 b 2 2 b 3 2 b 1 2 b 2 2 b 3 2
Table 11: RKF method Even though both RKS and RKF methods are third order accurate, the local time stepping scheme is only second order accurate because some third order coupling conditions are not satisfied [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. However, the local time stepping is mass conservative [START_REF] Seny | Multirate time stepping for accelerating explicit discontinuous galerkin computations with application to geophysical flows[END_REF]. Like for the original scheme of Constantinescu and Sandu, there is a buffer zone at the interface, in the slow subdomain. Outside this buffer zone, the RKS method reduces to the base method described in table 9. Finally, our semi-discrete approximation (4) is solved by the algorithm shown in table 12, in which the notations introduced in subsection 1.4 are re-used.

Slow subdomain

Fast subdomain

1 y n,1 i = y n i -a 21 ∆tD n i y n,1 i = y n i -a 21 ∆tD n i y n,1 i = y n i -a 21 2 ∆tD n i 2 y n,2 i = y n i -a 31 ∆tD n i y n,2 i = y n i -a 31 ∆tD n i y n,2 i = y n i -a 31 2 ∆tD n i -a 32 ∆tD n,1 i -a 32 ∆tD n,1 i -a 32 2 ∆tD n,1 i 3 y n,3 i = y n i y n,3 i = y n i -b 1 2 ∆tD n i -b 2 2 ∆tD n,1 i -b 3 2 ∆tD n,2 i 4 y n,4 i = y n i -a 21 ∆tD n,3 i y n,4 i = y n,3 i -a 21 2 ∆tD n,3 i 5 y n,5 i = y n i -a 31 ∆tD n,3 i y n,5 i = y n,3 i -a 31 2 ∆tD n,3 i -a 32 ∆tD n,4 i -a 32 2 ∆tD n,4 i 6 y n+1 i = y n i -b 1 ∆tD n i y n+1 i = y n i -b 1 2 ∆t(D n i + D n,3 i ) y n+1 i = y n,3 i -b 1 2 ∆tD n,3 i -b 2 ∆tD n,1 i -b 3 ∆tD n,2 i -b2 2 ∆t(D n,1 i + D n,4 i ) -b 3 2 ∆t(D n,2 i + D n,5 i ) -b 2 2 ∆tD n,4 i -b 3 2 ∆tD n,5 i

Buffer zone

Table 12: Algorithm of the second local time stepping scheme

The algorithm is written in "classical form" but it can be easily implemented in the low storage form since these two formulations are equivalent. Note that with a flux function with a four cell stencil, the size of the buffer zone is equal to 6 cells.

New schemes

In this section, we present a different approach to construct local time stepping schemes. Like in the strategy of Constantinescu and Sandu, we use the theory of PRK methods but here we guarantee the internal consistency of the scheme. This allows to obtain higher order accuracy, since we are able to design a third order scheme. Note that our schemes require a correction step to make them mass conservative. In these new schemes, the ratio between two adjacent time steps is restricted to 2. Several authors [START_REF] Liu | Nonuniform time-step runge-kutta discontinuous galerkin method for computational aeroacoustics[END_REF][START_REF] Liu | Nonuniform-time-step explicit runge-kutta scheme for high-order finite difference method[END_REF][START_REF] Tang | A class of high resolution schemes for hyperbolic conservation laws and convection-diffusion equations with varying time and space grids[END_REF][START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF] develop local time stepping schemes that can deal with ratios between 1 and 2 (2, 3/2, 4/3...) and even superior to 2 (3, 4...). In most meshes dedicated to LES or CFD in general, the grid size variation is rather continuous to guarantee the accuracy of the solution. The use of time step ratios superior to 2 is not appropriate since such grid size variations are not present. The use of several time step ratios between 1 and 2 allows a better adaptability to all grid size variations than the use of ratios restricted to 2. However, it tends to increase the algorithmic complexity. Indeed, the algorithm must be able to handle all possible ratios, while it has only 3 possibilities when the time step ratio is limited to 2. In this case, the neighbouring time step can be identical, equal to twice or equal to half the current time step. The choice to limit the ratio between adjacent time steps to 2 is a compromise between efficiency and algorithmic complexity.

First scheme: a second order time accurate scheme

Again, we use the second order Heun method (table 4) as the base method. It is applied twice with the time step ∆t 2 to obtain the RKF method (table 14). Our RKS method is composed of the base method and some additional values to perform the transition with the RKF method and achieve the desired order of accuracy. Thus, both RKS and RKF methods are of order 2. Now, we impose internal consistency by the relation c S i = c F i , for every i. It can be easily seen that second order coupling conditions (see table 2) are satisfied with this condition. Note that the intermediate time coefficients of RKF are already fixed: we have c F = {0; 1 2 ; 1 2 ; 1}. For RKS, we have c S 1 = 0 and c S 4 = 1, but c S 2 and c S 3 are still undetermined. Then the internal consistency is obtained with c S 2 = c S 3 = 1 2 . This requires to define additional intermediate values at time t n + 1 2 ∆t at the interface, in the slow subdomain. The coefficients a ij needed for the calculation of these additional value are determined by the consistency relation [START_REF] Liu | Nonuniform-time-step explicit runge-kutta scheme for high-order finite difference method[END_REF]. We remind that these additional values are only calculated at the interface and are not used for the calculation of y n+1 in the slow subdomain. The Butcher tableau of RKS and RKF methods are shown in Table 13 and Table 14.

0 0 1 2 1 2 0 1 2 1 2 0 0 1 1 0 0 0 1 2 0 0 1 2 Table 13: RKS method 0 0 1 2 1 2 0 1 2 1 4 1 4 0 1 1 4 1 4 1 2 0 1 4 1 4 1 4 1 4
Table 14: RKF method However, this local time stepping scheme defined by RKS and RKF is not mass conservative. It can be shown by applying RKS and RKF methods to (2) on cells i 0 and i 0+1 , respectively. We remind that cell i 0 is the last cell of the slow subdomain while cell i 0+1 is the first cell of the fast subdomain (see figure 1). The application of RKS and RKF methods to (2) on cells i 0 and i 0+1 is shown in table [START_REF] Sexton | [END_REF]. Note that to clarify the presentation, we assume that the flux function has a 2 cell stencil:

F i+ 1 2 = F i+ 1 2
(y i , y i+1 ). The following notations are introduced for this algorithm :

• for all n, the values at times t n and t n+1 (with t n+1 = t n + ∆t, ∆t corresponding to the largest time step) are denoted by y n and y n+1 ,

• for all n, the k -th intermediate value between times t n and t n+1 is denoted by y n,k .

• for all n, F n

i+ 1 2 = F i+ 1 2 (y n i , y n i+1 ),
• for all n, the k-th intermediate flux at the interface x i+ 1 2 is given by : F n,k

i+ 1 2 = F i+ 1 2 (y n,k i , y n,k i+1 ). step Solution on cell i 0 Solution on cell i 0+1 1 y n,1 i0 = y n i0 -∆t 2∆x i0 (F n i0+ 1 2 -F n i0-1 2 ) y 1 i0+1 = y n,1 i0+1 -∆t 2∆x i0+1 (F n i0+ 3 2 -F n i0+ 1 2 ) 2 y n,2 i0 = y n,1 i0 y n,2 i0+1 = y n i0+1 -∆t 4∆x i0+1 (F n i0+ 3 2 + F n,1 i0+ 3 2 -F n i0+ 1 2 -F n,1 i0+ 1 2 ) 3 y n,3 i0 = y n i0 -∆t ∆x i0 (F n i0+ 1 2 -F n i0-1 2 ) y n,3 i0+1 = y n,2 i0+1 -∆t 2∆x i0+1 (F n,2 i0+ 3 2 -F n,2 i0+ 1 2 ) 4 y n+1 i0 = y n i0 -∆t 2∆x i0 (F n i0+ 1 2 + F n,1 i0+ 1 2 y n+1 i0+1 = y 2 i0+1 -∆t 4∆x i0+1 (F n,2 i0+ 3 2 + F n,3 i0+ 3 2 -F n i0-1 2 -F n,1 i0-1 2 ) -F n,2 i0+ 1 2 -F n,3 i0+ 1 2
) We analyze the expression of y n+1 i0 . We notice that, between t n and t n + ∆t, the total flux computed with the RKS method through the interface x i0+ 1 2 is :

tf RKS i0+ 1 2 = 1 2 [F n i0+ 1 2 + F n,1 i0+ 1 2 ] (14)
The same analysis can be made on the other side of the interface. The total flux computed with the RKF method through the interface x i0+ 1 2 is:

tf RKF i0+ 1 2 = 1 4 [F n i0+ 1 2 + F n,1 i0+ 1 2 + F n,2 i0+ 1 2 + F n,3 i0+ 1 2 ] (15)
It is clear that total fluxes tf RKS i0+ 1 2

and tf RKF i0+ 1 2

are different and consequently the local time stepping scheme is not mass-conservative at the interface between the subdomains.

We propose a correction step to make the scheme conservative. Following the notations given above, y n+1 i0 reads :

y n+1 i0 = y n i0 - ∆t ∆x i0 tf RKS i0+ 1 2 + ∆t 2∆x i0 (F n i0-1 2 + F n,1 i0-1 2 ) ( 16
)
This value is corrected by adding the total flux difference. This leads to the following expression for the corrected value y n+1,c i0 :

y n+1,c i0 = y n+1 i0 + ∆t ∆x i0 tf RKS i0+ 1 2 - ∆t ∆x i0 tf RKF i0+ 1 2 (17)
We substitute y n+1 i0 by its expression given by ( 16). The corrected value y n+1,c i0 now reads:

y n+1,c i0 = y n i0 - ∆t ∆x i0 tf RKF i0+ 1 2 + ∆t 2∆x i0 (F n i0-1 2 + F n,1 i0-1 2 ) (18)
In [START_REF] Knoth | Implicit-explicit runge-kutta methods for computing atmospheric reactive flows[END_REF], the total flux computed by the RKS method through the interface

x i0+ 1 2 , tf RKS i0+ 1 2
, has been replaced by the total flux computed by the RKF method through the interface

x i0+ 1 2 , tf RKF i0+ 1 2
. Consequently, the corrected scheme is mass conservative. Time order and spatial order studies made on the different test cases in section 3 show that this correction step does not alter the time order of accuracy nor the spatial order of accuracy of the solution.

The algorithm of this local time stepping applied to our semi-discrete approximation (4) is shown in table [START_REF] Sandu | A class of multirate infinitesimal gark methods[END_REF]. In this algorithm, the notations introduced in subsection 1.4 are re-used.

Slow subdomain

Fast subdomain

1 y n,3 i = y n i -∆tD n i y n,1 i = y n i -∆t 2 D n i y n,1 i = y n i -∆t 2 D n i 2 y n,2 i = y n,1 i y n,2 i = y n i -∆t 4 [D n i + D n,1 i ] 3 y n,3 i = y n i -∆tD n i y n,3 i = y n,2 i -∆t 2 D n,2 i 4 y n+1 i = y n i -∆t 2 [D n i + D n,3 i ] y n+1 i = y n i -∆t 2 [D n i + D n,3 i ] y n+1 i = y n,2 i -∆t 4 [D n,2 i + D n,3 i ] 5
Correction step (last cell)

Buffer zone

Table 16: Algorithm of the first proposed scheme (second order time accurate)

In our applications, with a four cell stencil flux function, the size of the buffer zone is equal to two cells, while it is equal to four cells for the original scheme of Constantinescu and Sandu. Moreover, only two flux evaluations (D n i and D n,3 i ) are needed in the buffer zone whereas the original scheme of Constantinescu and Sandu requires four flux evaluations in this zone. Although its buffer zone is larger, Constantinescu and Sandu scheme is mass-conservative and TVB (Total Variation Bounded) if both temporal and spatial scheme are TVD (Total Variation Diminishing) [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. Our scheme requires a correction step to be mass conservative and we didn't design our scheme to be TVB or TVD.

Second scheme: a third order time accurate scheme

The base method used to construct this local time stepping scheme is the RK3 low storage method described in table 9. The construction of RKF method is performed with two successive applications of the base method with the time step ∆t 2 . Consequently, RKF method has the following intermediate time coefficients: c F = {0; c 2 2 ; c 3 2 ; 1 2 ; 1+c 2 2 ; 1+c 3 2 }. For the construction of RKS method, the base method is applied once with a time step ∆t, which leads to the following intermediate time coefficients: c S = {0; c 2 = 1 2 ; c 3 }. As for the first proposed scheme, we impose internal consistency by the relations c S i = c F i , for every i. Indeed, this property implies that some of the third order coupling conditions of table 3 are necessarily satisfied (namely that of lines 1,4,5). To satisfy internal consistency, an additional value at time t + c 3 ∆t must be calculated in RKF method. In RKS method, additional values at times t + c 2 2 ∆t, t + c 3 2 ∆t, t + 1+c 2 2 ∆t and t + 1+c 3 2 ∆t must be computed. RKS and RKF methods are now two third order Runge-Kutta methods with 7 stages. Their intermediate time coefficients are the same and read: { 0; c 2 2 ; c 3 2 ; 1 2 ; 1+c 2 2 ;c 3 ; 1+c 3 2 }. The Runge-Kutta coefficients for the calculation of additional values, that are still to be determined, are denoted by α S in RKS method and by α F in RKF method. The Butcher tableaus of RKS and RKF methods are shown in tables 17 and 18, respectively. Table 18: RKF method

0 0 c 2 2 α S 1 c 3 2 α S 2 α S
The unknown coefficients α S i and α F i will be determined so that the coupling conditions of lines 2 and 3 of table 3 are satisfied. Some of them are equivalent, therefore the independent coupling conditions read:

s i=1 s j=1 b S i a F ij c F j = 1 6 , (19)
s i=1 s j=1 b F i a S ij c S j = 1 6 , (20)
The coefficients for the computation of additional values in RKS method have to satisfy [START_REF] Hairer | Solving ordinary differential equations I : Nonstiff Problems[END_REF] and four relations given by the consistency property [START_REF] Liu | Nonuniform-time-step explicit runge-kutta scheme for high-order finite difference method[END_REF]. This gives 5 relations that can be used to determine 5 unknowns only. The other eight coefficients can be set to any arbitrary value. A simple and natural choice is to set: α S 3 = α S 5 = α S 6 = α S 7 = α S 9 = α S 10 = α S 12 = α S 13 = 0, while the resolution of the system leads to the following values for the other coefficients:

                           α S 1 = c 2 2 α S 2 = c 3 2 α S 4 = 1+c 2 2 α S 11 = 2 3b 3 ≈ 1.821 α S 8 = 1+c 3 2 -2 3b 3 ≈ -0.9270 (21)
A similar analysis leads to the following values for the coefficients dedicated to the calculation of additional values in RKF method:

                           α F 2 = ( 2 3 -b 2 2 ) 1 b 3 ≈ 1.2440 α F 1 = c 3 -α F 2 ≈ -0.45534 α F 3 = 0 α F 4 = 0 α F 5 = 0 (22)
As for our previous second order scheme, the local time stepping scheme defined by RKS and RKF methods described in tables 17 and 18 is not mass conservative. The strategy based on the Total Flux Lost and the Total Flux Received at the interface between subdomains can be applied to make the scheme conservative. The algorithm of this local time stepping scheme applied to solve (4) is shown in table [START_REF] Wensch | Multirate infinitesimal step methods for atmospheric flow simulation[END_REF]. In this algorithm, the notations introduced in subsection 1.4 are re-used.

Slow subdomain

Fast subdomain This algorithm is written in the Runge-Kutta "classical form" but it can be implemented in the "low-storage form", since these forms are equivalent. In this scheme, there are two buffer zones because both RKS and RKF methods require one or several additional values. The size of each buffer zone is equal to two cells in our applications (flux function with a four cell stencil). We remind that for the scheme of Constantinescu and Sandu based on RK3, there is only one buffer zone of six cells and no correction step since this scheme is mass conservative. Moreover, if both temporal and spatial scheme are TVD, the scheme of Constantinescu and Sandu is TVB, which is not the case for our scheme.

1 y n,1 i = y n i -α S 1 ∆tD n i y n,1 i = y n i -a 21 ∆t 2 D n i y n,1 i = y n i -a 21 ∆t 2 D n i 2 y n,2 i = y n i -a 31 ∆t 2 D n i y n,2 i = y n i -a 31 ∆t 2 D n i y n,2 i = y n i -α S 2 ∆tD n i -a 32 ∆t 2 D n,1 i -a 32 ∆t 2 D n,1 i 3 y n,3 i = y n i -b 1 ∆t 2 D n i y n,3 i = y n i -b 1 ∆t 2 D n i y n,3 i = y n i -a 21 ∆tD n i y n,3 i = y n i -a 21 ∆tD n i -b 2 ∆t 2 D n,1 i -b 3 ∆t 2 D n,2 i -b 2 ∆t 2 D n,1 i -b 3 ∆t 2 D n,2
n i -α F 1 ∆tD n i -a 32 ∆tD n,3 i -a 32 ∆tD n,3 i -α F 2 ∆tD n,1 i 6 y n,6 i = y n i -α S 8 ∆tD n i y n,6 i = y n,3 i -a 31 ∆t 2 D n,3 i y n,6 i = y n,3 i -a 31 ∆t 2 D n,3 i -α S 11 ∆tD 4 i -a 32 ∆t 2 D 4 i -a 32 ∆t 2 D 4 i 7 y n+1 i = y n i -b 1 ∆tD n i y n+1 i = y n i -b 1 ∆tD n i y n+1 i = y n,3 i -b 1 ∆t 2 D n,3 i y n+1 i = y n,3 i -b 1 ∆t 2 D n,3 i -b 2 ∆tD n,3 i -b 3 ∆tD n,5 i -b 2 ∆tD n,3 i -b 3 ∆tD n,5 i -b 2 ∆t 2 D n,4 i -b 3 ∆t 2 D n,6 i -b 2 ∆t 2 D n,4 i -b 3 ∆t 2 D n

Numerical tests

In this section, some numerical test cases are performed to compare the new proposed schemes and the schemes proposed in the literature that are based on the strategy of Constantinescu and Sandu. The schemes based on RK2 and RK3 which follow the strategy of Constantinescu and Sandu are denoted by CSRK2 and CSRK3, respectively. The schemes that we propose, based on RK2 and RK3, are denoted by Prop. RK2 and Prop. RK3, respectively.

Two different schemes are used for the Euler flux discretization. The Roe scheme [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF] with a minmod limiter [START_REF] Roe | Characteristic-based schemes for the euler equations[END_REF] is retained for the Sod tube, whereas an hybrid centered/upwind version of the scheme AUSM+(P) proposed by Mary and Sagaut [START_REF] Mary | Large eddy simulation of flow around an airfoil near stall[END_REF] is employed in the others test cases. Viscous fluxes are discretized with a classical second order scheme.

Two-dimensional vortex advection

This test case is based on 2D compressible Euler equations for a perfect gas: The initial state is a 2D uniform flow (ρ = 1, u 1 = 1, u 2 = 0, p = 1) with the superposition of perturbations (δρ, δu 1 , δu 2 , δp) which simulate a vortex. The perturbations read:

           ∂ t ρ + ∂ x j (ρu j ) = 0 ∂ t ρu i + ∂ x j (ρu i u j + pδ ij ) = 0 ∂ t (ρe) + ∂ x j (ρeu j + pu j ) = 0
                   δρ = (1 + δT ) 1 γ-1 δu 1 = 2π e 0.5(1-r 2 ) (y c -y) δu 2 = 2π e 0.5(1-r 2 ) (x -x c ) δp = (1 + δT ) γ γ-1 ( 24
)
where δT = -(γ-1) 2 8γπ 2 e 1-r 2 , is the temperature perturbation around the value T = 1. Note that the corresponding initial value and perturbation for the entropy are S = 1 and δS = 0, respectively. The coordinates (x c , y c ) are the coordinates of the the vortex center, r = (x -x c ) 2 + (y -y c ) 2 is the vortex radius and = 5 is the vortex strength. The initial solution of this problem reads: scheme based on RK2 are identical and show second order accuracy. For the schemes based on RK3, our proposed approach is really third order accurate, while CS RK3 is only a second-order accurate scheme leading to a higher level of L 2 error. Figure 7 shows the results of a space order study, for each local time stepping scheme. For this study, the spatial steps ∆x and ∆y are uniform (no space refinement in the middle zone). There is only a temporal refinement in the middle zone. The finest solution is computed with ∆x = ∆y = 0.013. The time step is taken sufficiently small in order to compute only the spatial error. The time step ∆t = 9.10 -4 (CF L = 0.1) is used in the first and last zones, while the time step ∆t = 4, 5.10 -4 is used in the middle zone. Three other solutions are computed keeping the same time steps and with the spacial steps (2∆x,2∆y), (3∆x,3∆y), (4∆x,4∆y). The L 2 errors on the density field between the finest solution and the three other solutions are computed. The results show that the four local time stepping schemes have the same spatial order of accuracy. This order of accuracy is close to 2, the theoretical spatial order of accuracy for the AUSM+P scheme. Figures 6 and7 demonstrate that our flux correction stage used in schemes "Prop. RK2" and "Prop. RK3" does not alter the time order of accuracy nor the spatial order of accuracy.

Sod shock tube

The Sod tube problem is used to test if the proposed schemes are able to predict accurate and stable solution for a flow containing discontinuities. Note that this 1D problem is solved here in 2D. At the initial time the velocity field is equal to zero, while the density and pressure profiles are given by :

• ρ 0 (x, y) = 1, p 0 (x, y) = 1 for x ≤ 9, • ρ 0 (x, y) = 0.125, p 0 (x, y) = 0.1 for x > 9,
The rectangular computational domain (18 × 0.6) is divided into three zones of length 6. In the first and last zones, the grid spacings in the x direction (∆x) and in the y direction (∆y) are: ∆x = ∆y = 0.01. In the middle zone, ∆x = 0.005 and ∆y = 0.01. Thus, the time step ∆t 2 is used in the middle zone and the time step ∆t is applied in the first and last zones. The mesh partitioning and the initial density profile at y = 0.3 are plotted in figure [START_REF] Schlegel | Multirate runge-kutta schemes for advection equations[END_REF].

Figures 9 and10 show the density profile at y = 0.3 for each local time stepping scheme and for the reference solution at the final time, t= 4.8. The reference solution mentioned in these figures is obtained by using the classical RK4 explicit scheme with the time step ∆t = 5.10 -4 . The CF L number in the middle zone (the maximum CF L number over the mesh) is CF L = 0.2. The time step ∆t of the local time stepping schemes (the time step in the first and last zones) is their maximal allowable time step shown in table 21 below. Figures 9 and10 show a very good agreement between the reference solution and the solutions obtained with the local time stepping schemes. Table 21 also shows the speedup as compared to a RK2 global time stepping scheme, which uses a uniform time step satisfying the most restrictive CFL condition. Our local time stepping schemes have roughly the same speedup as their homologous schemes of Constantinescu and Sandu based on the same RK method.

The results of a time order study are plotted in figure 11. For each local time stepping scheme, a finest solution is computed with the time step ∆t = 2.5.10 -4 in the first and last zones. We remind that in the middle zone the grid spacing and the time step are refined by a factor 2, therefore the CF L number is uniform over the mesh (CF L = 0.05). Three other simulations, with time steps equal to 2∆t, 4∆t and 8∆t for the first and last zones are also run. The L 2 errors committed on the density field between the finest solution and the three other solutions are computed. The slope of each convergence curve is in agreement with the theoretical temporal order of the local time stepping schemes. Our local time stepping scheme "Prop. RK3" is third-order accurate and leads to the lowest error. Figure 12 shows the results of a space order study, for each local time stepping scheme. This study is performed by employing the same method as in subsection 3.1. The finest solution is computed with ∆x = ∆y = 0.0057. The time step ∆t = 5.10 -4 (CF L = 0.1) is used in the first and last zones, while the time step ∆t = 2, 5.10 -4 is used in the middle zone. The results show that the four local time stepping schemes have the same spatial order of accuracy. This order of accuracy is slightly superior to 1. This result is logical since our third order accurate Roe scheme uses a slope limitor which reduces the order of accuracy at discontinuities. Figures 11 and12 demonstrate that our flux correction stage used in schemes "Prop. RK2" and "Prop. RK3" does not alter the time order of accuracy nor the spatial order of accuracy.

Taylor-Green Vortex

This test case is used to validate the proposed schemes in the context of DNS or LES simulation of turbulent flows. At the initial time, large vortices are present in the computational box and during the simulation an energy transfer occurs from the large scales to smaller ones. When the energy cascade is sufficient to reach the Kolmogorov scale, the kinetic energy is dissipated. For that case, the compressible 3D Navier-Stokes equations are solved for a viscous fluid :

           ∂ t ρ + ∂ x j (ρu j ) = 0 ∂ t ρ + ∂ x j (ρu i u j + pδ ij -τ ij ) = 0 ∂ t (ρe) + ∂ x j (ρeu j + pu j -τ jk u k + q j ) = 0 (26)
with 1 ≤ i, j ≤ 3. The variable ρ is the fluid density, u i are the fluid velocity components, e is the total specific energy of the fluid, p is the fluid pressure (ideal gas law), q j is the heat flux, and τ ij is the tensor of viscous constraints. The computational domain is a cube of dimension 2πL 0 (L 0 = 1 is a reference length), with periodic boundary conditions. The initial conditions for variables u i , p and ρ read :

                         u 0 1 = U 0 sin(x/L 0 ) cos(y/L 0 ) cos(z/L 0 ) u 0 2 = -U 0 cos(x/L 0 ) sin(y/L 0 ) cos(z/L 0 ) u 0 3 = 0 p 0 = P 0 + ρ 0 U 2 0 /16(cos(x/L 0 ) + cos(y/L 0 ))(cos(2z/L 0 ) + 2) ρ 0 = p/(RT 0 ) (27)
where P 0 = 101183P a, ρ 0 = 1.2kg.m -3 , T 0 = 294K and U 0 = 34.38m.s -1 . The simulation is performed during a period of 14L 0 /U 0 . The computational domain is divided into three zones of equal length, as shown in figure 13. In the first and last zones, the grid spacings in the x direction (∆x), y direction (∆y) and z direction (∆z) are ∆x = ∆y = ∆z = 0.025, respectively. In the middle zone, the mesh is refined by a factor two in the x direction (∆x = 0.0125, ∆y = ∆z = 0.025). In this zone, the time step ∆t 2 is used while the time step ∆t is applied in the other zones. For each scheme, we compute at each time step the enstrophy , defined by

= L 0 U 2 0 ρ 0 Ω ρ w. wdΩ, (28)
where w denotes the vorticity, defined by : w = ∇ ∧ V , with V = (u 1 , u 2 , u 3 ). The enstrophy is known to be highly sensitive to the numerical method [START_REF] Debonis | Solutions of the taylor-green vortex problem using high-resolution explicit finite difference methods[END_REF]. Figures 14 and15 show the time evolution of the enstrophy for each scheme. The reference solution mentioned in these figures is obtained with a classical RK4 explicit scheme used with the time step ∆t = 0.00015. The CF L number in the middle zone (the maximum CF L number over the mesh) is CF L = 0.03. The time step ∆t in the first and last zones of each local time stepping scheme is the maximal allowable time step, shown in the first row of table 22. A very good agreement is observed between the reference solution and the solutions obtained with the local time stepping schemes.

Table 22 also shows the speedup obtained with the local time stepping schemes as compared to an RK2 explicit global time stepping scheme, and the conservation defects. Our local time stepping schemes provide the same speedup as their homologous schemes based on the strategy of Constantinescu and Sandu. The conservation defects of our schemes are very close to the ones obtained with the schemes of Constantinescu and Sandu. It confirms that our flux-correction strategy is efficient.

Finally, figure 16 shows the results of a time order study. For each local time stepping scheme, the "finest" solution is computed with the time step ∆t = 7.10 -5 in the first and last zones. We remind that in the middle zone the grid spacing and the time step are refined by a factor 2, therefore the CF L number is uniform over the mesh (CF L = 0.007). In this figure, each mark is associated to the use of the time steps : 2∆t, 4∆t and 8∆t, for the first and last zones. Once again the slope of each curve is in agreement with the theoretical temporal order of each local time stepping scheme. Our proposed scheme based on RK3 has the lowest L 2 error. Figure 17 shows the results of a space order study, for each local time stepping scheme. This study is performed by employing the same method as in subsection 3.1. The finest solution is computed with ∆x = ∆y = 0.022. The time step Prop. RK3 0.9 (∆t max = 0.0019) 59% 3, 8.10 -15

Turbulent flow over an airfoil

This test case is chosen to assess the interest of the local time stepping method for a flow configuration representative of current challenging configurations in the aerospace context. For that purpose a Large Eddy Simulation (LES) over a SD7003 profile is performed. Due to the angle of attack and the Reynolds number, a laminar-turbulent transition occurs in a laminar separation bubble (LSB).

It is followed by a turbulent reattachment and a classical turbulent boundary layer subjected to an adverse pressure gradient. As this flow configuration is very sensitive to physical perturbations due to the laminar-turbulent transition process [START_REF] Marxen | The effect of small-amplitude convective disturbances on the size and bursting of a laminar separation bubble[END_REF], it represents a relevent test case to highlight possible unwanted numerical perturbations introduced by our proposed temporal scheme. This test case is based on 3D Navier-Stokes equations [START_REF] Mary | Large eddy simulation of flow around an airfoil near stall[END_REF]. The Mach number, M , and the Reynolds number Re based on the airfoil chord c and the infinite velocity U ∞ are set to M = 0.1 and Re = 60000, respectively. The angle of attack α is set to α = 8 • . This Reynolds number is not so challenging for nowadays simulations, but it is sufficient to test the numerical method with a reasonable CPU cost for a comparison between different schemes.

The computational domain is a 3D structured mesh with approximately 90.10 6 cells. The spanwise extent is equal to 0.2c, while the outer boundaries are located at 100c from the airfoil. Characteristic non reflective boundary conditions are used for the farfield, an adiabatic wall is considered for the profile, and periodic boundary conditions are used in the spanwise z direction. To trigger the transition process in a deterministic way, the perturbation proposed in [START_REF] Marxen | The effect of small-amplitude convective disturbances on the size and bursting of a laminar separation bubble[END_REF] is imposed on the suction side at x/c = 0.03 with an amplitude A v = 0.001U ∞ . Figure 18 shows a mesh section in the ( x, y) map with a reduced number of points (one point over 10). An extruded 2D RANS solution is used to initialize the 3D flow solution. No reliable experimental data are available to validate the simulation since the size of the LSB is strongly influenced by the level of the freestream turbulence and also by sidewall effects of the windtunnel. For a computational cost reason, taking into account properly these phenomena is out of the scope of the present simulation. The present wall mesh resolution, ∆x + ≈ 25, ∆y + ≈ 1, ∆z + ≈ 8, is finer than that used in [START_REF] Mary | Large eddy simulation of flow around an airfoil near stall[END_REF] and equivalent to that used in [START_REF] Alferez | Study of stall development around an airfoil by means of high fidelity large eddy simulation[END_REF][START_REF] Laurent | Dns database of a transitional separation bubble on a flat plate and application to rans modeling validation[END_REF], in which grid converged results are presented for the same kind of flow: airfoil with laminar separation bubble and turbulent reattachment. As the numerical method employed in the present paper is the same as the one validated in [START_REF] Mary | Large eddy simulation of flow around an airfoil near stall[END_REF][START_REF] Larchevêque | Large-eddy simulation of a compressible flow past a deep cavity[END_REF][START_REF] Alferez | Study of stall development around an airfoil by means of high fidelity large eddy simulation[END_REF][START_REF] Laurent | Dns database of a transitional separation bubble on a flat plate and application to rans modeling validation[END_REF][START_REF] Daude | Self-adaptive newton-based iteration strategy for the les of turbulent multi-scale flows[END_REF][START_REF] Dandois | Large-eddy simulation of laminar transonic buffet[END_REF] thanks to windtunnel experiment comparisons, we are confident that our results obtained with the classical RK3 scheme and a global timestep are physically correct for the specific set of Reynolds number and boundary conditions. Therefore the solution obtained with the classical RK3 explicit scheme is considered as a reference to assess the accuracy and efficiency of the proposed temporal scheme.

Our most accurate local time stepping scheme "Prop. RK3" is used to compute this test case. It is compared to a classical explicit scheme (RK3) and one implicit scheme in terms of accuracy (see figures 23 -28) and CPU cost (see table 24). The implicit scheme is based on the second order accurate Gear scheme [START_REF] Gear | Numerical initial value problems in ordinary differential equations[END_REF] for the approximation of the time derivative. A self-adaptive Newton method is used to solve the non linear problem [START_REF] Daude | Self-adaptive newton-based iteration strategy for the les of turbulent multi-scale flows[END_REF]. This method takes advantage of Newton convergence rate heterogeneities by using a spatially varying number of iterations for the Newton process. A user defined convergence criterion drives the spatially varying number of iterations needed to solve the Newton process. Two implicit simulations are performed with different values for the convergence criterion of the Newton process: = 0.1, and = 0.06 (see [START_REF] Daude | Self-adaptive newton-based iteration strategy for the les of turbulent multi-scale flows[END_REF] for more information about this convergence criterion). Note that for this test case, all the schemes are implemented with parallel OpenMP directives and run on 64 threads on one Intel Xeon Phi 7230 node.

The initial 2D RANS solution is used for an automatic mesh partitioning based on the local value of the CFL number. Figure 19 shows a zoom on the airfoil to illustrate the spatial repartition of the different time steps for our local time stepping scheme. Only four time levels have been considered for this validation and the largest time step is such that ∆t × U ∞ /c = 1.10 -5 . This table 23 shows the relative number of cells for each time step. This table 23 illustrates the interest of using a local time stepping scheme for this case. Indeed, it seems inefficient to restrict the time step for the whole domain to ∆t 8 , like in a classical explicit scheme, while only 4% cells require to use this time step. Relative number of cells (%) 61 23 12 4 Table 23: Relative number of cells for each time step

The time step used for the classical RK3 scheme is such that ∆t × U ∞ /c = 1.25.10 -6 and the time step of the implicit simulations is such that ∆t × U ∞ /c = 6.10 -5 (largest CFL number of e for each scheme in the transition, reattachment and mid-chord zones.

Figure 23 shows a perfect agreement between all schemes for the mean streamwise velocity profiles. The same perfect agreement is observed for the Reynolds stress profiles between the local time stepping scheme, the implicit scheme with = 0.06 and the reference solution (global time stepping scheme). Small discrepancies between the reference solution and the implicit scheme run with = 0.1 are present for the Reynolds stresses profiles. It shows that this value of is too large to prevent an interaction between the flow physics and the numerical error due to the incomplete convergence of the Newton process. No significant influence of the temporal scheme is observed on the pressure spectrums of figures 26 -28. A slight offset between the implicit schemes spectrums and the explicit schemes spectrums is observed for high Strouhal numbers (St > 10 3 ) at points ( 2) and (3). The implicit schemes seem slightly less accurate than the explicit schemes to compute high frequencies, certainly because of the use of a larger time step and the dissipative nature of the implicit Gear scheme.

Table 24 shows the CPU time obtained for each scheme and the time reduction as compared to the explicit global time stepping scheme. It appears that our local time stepping scheme "Prop. RK3" provides a significant CPU time reduction (69%), whitout altering the accuracy of the solution. Same order of CPU time reduction is obtained with the implicit scheme using = 0.1 (73%), whereas the implicit simulation using = 0.06 only provides a 47% CPU time reduction. This test case demonstrates that our local time stepping scheme "Prop RK3" is able to compute a challenging, industrial-like test case with a very good accuracy and for a competitive CPU time. In order to reach the accuracy of our local time stepping scheme, the implicit scheme convergence criterion has to be reduced from = 0.1 to = 0.06 (see figures 24 -25). It leads to an important increase in the CPU time, which makes the implicit scheme less competitive than our explicit local time stepping scheme in this field.

Conclusion

The diversity of spatial scales in fluid mechanics problems leads to the use of several cell sizes in CFD meshes. It leads to a wide range of numerical stability conditions that must be satisfied by explicit time integration methods. Local time stepping schemes are well suited to these problems because the time step can vary across the mesh to satisfy several local stability conditions. It results in a computational speedup as compared to classical explicit schemes, where the time step for the entire domain is restricted by the most severe stability condition. However, the construction of local time stepping schemes requires a particular attention to ensure high order time accuracy and mass conservation.

Constantinescu and Sandu [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF] have developed an interesting strategy which allows to construct mass conservative and second order accurate local time stepping schemes. Their strategy is based on the theory of Partitioned Runge-Kutta (PRK) methods. We also use the theory of PRK methods as well as the computation of consistent additional values at the interfaces to develop two new local time stepping schemes of order two and three. Our schemes are made conservative with an additional flux-correction step. The academic tests show that very similar results are obtained with our local time stepping schemes and the local time stepping schemes of [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. Moreover, the tests confirm that our third order local time stepping scheme has a lower temporal error than our second order local time stepping scheme and than the local time stepping schemes of [START_REF] Constantinescu | Multirate timestepping methods for hyperbolic conservation laws[END_REF]. The final test case shows that our third order local time stepping scheme is able to compute complex, industrial-like test cases with a very good accuracy and for a CPU time close to the one obtained with our implicit scheme. Our strategy seems to be a good approach to construct efficient, accurate and conservative local time stepping schemes.

We believe that our method will be more efficient than implicit schemes on a larger number of CPU cores, thanks to MPI or hybrid MPI/OpenMP parallelization. This requires a specific work on load balancing for our local time stepping scheme, like in [START_REF] Seny | An efficient parallel implementation of explicit multirate runge-kutta schemes for discontinuous galerkin computations[END_REF] and [START_REF] Rietmann | Load-balanced local time stepping for large scale wave propagation[END_REF]. This is a work in progress that will be the subject of a future publication.

Figure 1 :

 1 Figure 1: Partitioned computational domain

1

 1 ≤ i, j ≤ 2. The variables ρ, u 1 , u 2 , p are the fluid density, the fluid velocity in directions x and y and the fluid pressure, respectively. Total energy is denoted by e and reads : e = 1

Figure 2 :

 2 Figure 2: Density field at initial time and view of the mesh partitioning for the local time stepping computation.Figure 3: Initial density profile at y = 10

Figure 3 :

 3 Figure 2: Density field at initial time and view of the mesh partitioning for the local time stepping computation.Figure 3: Initial density profile at y = 10

Figure 4 :

 4 Figure 4: Density profiles at y = 10 and final time for the schemes Prop.RK2 and CS RK2 and the reference solution

Figure 5 :

 5 Figure 5: Density profiles at y = 10 and final time for Prop.RK3, CS RK3 schemes and the reference solution

Figure 6 :Figure 7 : 20 CFLmax

 6720 Figure 6: L 2 error as a function of the time step ∆t for each local time stepping schemeFigure 7: L 2 error as a function of the spatial step ∆x for each local time stepping scheme

Figure 8 :

 8 Figure 8: Initial density profile and view of the mesh partitioning for local time stepping scheme.

Figure 9 :

 9 Figure 9: Density profile at the final time for the CS RK2 and Prop. RK2 schemes Figure 10: Density profile at the final time for the CS RK3 and Prop. RK3 schemes.

Figure 11 :Figure 12 :

 1112 Figure 11: L 2 error as a function of the time step ∆t for each local time stepping scheme

Figure 13 :

 13 Figure 13: Mesh used for the Taylor-Green Vortex test case

Figure 14 :Figure 15 :

 1415 Figure 14: Enstrophy obtained with the schemes "CS RK2" and "Prop. RK2" Figure 15: Enstrophy obtained with the schemes "CS RK3" and "Prop. RK3"

Figure 16 :Figure 17 :

 1617 Figure 16: L 2 error as a function of the time step ∆t for each local time stepping scheme Figure 17: L 2 error as a function of the space step ∆x for each local time stepping scheme

Figure 20 :

 20 Figure 20: Q-criterion iso-surface (Q = 1.5) coloured by the fluid velocity magnitude

Figure 21 :Figure 22 :Figure 23 :

 212223 Figure 21: Location of transition and reattachment zones with their pressure recording points Figure 22: Location of the mid-chord zone with its pressure recording point

Figure 24 :

 24 Figure 24: Mean Reynolds stress profiles < ρu u > /ρu 2 e for each scheme in the transition, reattachment and mid-chord zones.

Figure 25 :

 25 Figure 25: Mean Reynolds stress profiles < ρu v > ρu 2e for each scheme in the transition, reattachment and mid-chord zones.

Figure 26 :

 26 Figure 26: PSD spectrum as a function of the Strouhal number St (based on U ∞ ) for each scheme in the transition zone (point 1)

Figure 27 :

 27 Figure 27: PSD spectrum as a function of the Strouhal number St (based on U ∞ ) for each scheme in the reattachment zone (point 2)

Table 1 :

 1 Butcher tableau of an s-stage explicit Runge-Kutta method All Runge-Kutta methods in this article satisfy the property of consistency :

Table 7 :

 7 Algorithm of the original scheme of Constantinescu and Sandu

Table 15 :

 15 RKS and RKF methods applied on cells i 0 and i 0+1

Table 19 :

 19 Algorithm of the second proposed scheme (third order time accurate)

	,6
	i

Table 20 :

 20 Comparison of efficiency and conservation property indicators for the local time stepping schemes

Table 22 :

 22 Taylor-Green Vortex: comparison of some indicators for the different local time stepping schemes

Table 24 :

 24 Global time steppping Local time stepping Imp. scheme Imp. scheme Temporal scheme influence on the CPU cost of the SD7003 airfoil simulation on a single Intel Xeon Phi 7230 node

		(RK3)	(Prop. RK3)	= 0.1	= 0.06
	CPU time (hours)	653	197	175	347
	Time reduction (%)	0	69	73	47

Acknowledgment

This work has been supported by ONERA (Office Nationale d'Etudes et de Recherches Aérospatiales) and DGA (Direction Générale de l'Armement), grant No. CC1750701.

The solution of this problem is a simple advection of the vortex in the x direction. For more information about this test case, see [START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF]. The computational domain is rectangular (dimensions : 30 × 20) with periodic boundary conditions. It is divided into three zones. In the first and last zones, the grid spacings in directions x and y are ∆x = ∆y = 0.05. In the middle zone, the grid spacing in the x direction is refined by a factor 2 (∆x = 0.025 and ∆y = 0.05). Thus, the time step ∆t 2 is used in this zone and the time step ∆t is applied in the first and last zones. The mesh and the density field at initial time are shown in figure 2, while the density profile at initial time and y = 10 is plotted in Figure 3. The reference solution mentioned in figures 4 and 5 is obtained with a fourth order Runge-Kutta scheme (RK4) with ∆t = 0.0035. The CF L number in the middle zone (the maximum CF L number over the mesh) is CF L = 0.28. For each local time stepping scheme, the time step used in the first and last zones (∆t) is the maximum allowable time step shown in the row "CFLmax (∆t max )" of table [START_REF] Hairer | Solving ordinary differential equations I : Nonstiff Problems[END_REF].

Figures 4 and5 show the density profile along y = 10 for each local time stepping scheme and for the reference solution at the final time 43.5. At this time, the vortex has returned approximately to its initial position (we remind that we are using periodic boundary conditions). A very good agreement between each scheme and the reference solution is observed.

Some characteristics studied at the final time are listed in table 20. The "speedup" refers to the relative difference between the CPU time of an explicit global time stepping scheme based on RK2 (which uses a uniform time step satisfying the most restrictive CFL condition) and the CPU time of the local time stepping schemes. The row "Conservation defect" is the relative difference between the total fluid mass in the computational domain at final time m(t f ) and at initial time m(t 0 ) :

The local time stepping schemes based on RK3 have a better stability than the local time stepping schemes based on RK2, which allows to use larger time steps. It is the reason why the local time stepping schemes based on RK3 (CS RK3 and Prop. RK3) provide a better speedup than the local time stepping schemes based on RK2 (CS RK2 and Prop. RK2). The conservation defects obtained with our local time stepping schemes are similar to the ones obtained with the schemes of Constantinescu that are conservative. It shows that our flux-correction strategy to make our schemes conservative is efficient. Figure 6 shows the results of a time order study. For each local time stepping scheme, a finest solution is computed with the time step ∆t = 10 -4 in the first and last zones. We remind that in the middle zone the grid spacing and the time step are refined by a factor 2, therefore the CF L number is uniform over the mesh (CF L = 0.008). Three other simulations, with time steps equal to 2∆t, 4∆t and 8∆t for the first and last zones are also run. The L 2 errors on the density field between the finest solution and the three other solutions are computed. The slope of each convergence curve is in agreement with the theoretical temporal order of all the local time stepping schemes. We can notice that the L 2 errors of the original scheme of Constantinescu and Sandu and our proposed roughly 30). The simulations with the different schemes are performed until the final time t f such that t f × U ∞ /c = 2 is reached.

Figure 20 shows the iso-surface of the Q-criterion (Q = 1.5) coloured with the fluid velocity magnitude, obtained with our explicit local time stepping scheme at the final time t f . The Qcriterion is an indicator of turbulent structures, defined by the formula:

where U = (u 1 , u 2 , u 3 ) is the fluid velocity vector. This figure shows that the flow is laminar at the leading edge, then a laminar-turbulent transition zone is observed. Finally, the flow is fully turbulent. Figures 21 and22 show instantaneous isolines of the streamwise velocity component in the plane z/c=0. Figure 21 shows that the transition occurs in a laminar separation bubble, characterized by a fluid recirculation (negative streamwise velocity component). The turbulent structures created in the transition zone allow to reattach the fluid to the airfoil (reattachment zone).


Noureddine Oubeha
email: nourdine.oubeha@usmba.ac.ma

Karim El Khadiri

Hassan Qjidaa
email: qjidah@yahoo.fr

Li-Ion

Karim El Khadiri
email: karim.elkhadiri@usmba.ac.ma

Li-Ion battery charger for Systems RF wireless

Keywords: Li-Ion charger, Multi-mode, LDO

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

With the variety of electronic devices are generally existed in our daily life, many kinds of charger in accordance with those electronic devices emerged. Most chargers only provide fixed voltage. Hence users need to purchase a variety of rechargeable battery with specification which reduces the convenience greatly. Besides, the size of the portable products gets smaller.

The technology of global universal charger battery has been approved, most systems RF wireless have to use integrated charger. Implementation of unified chargers can be said a good way to solve multiple problems, it can not only eliminate the threat of health to consumers, but also benefits to environmental protection and the development of industries. When the battery charger is designed with switching-mode power supply (SMPS) converters, the passive elements are too huge to put them into chip. Also, the battery chargers based on the switched capacitor(SC) topology have the same problem, so the circuitsize is not suitable for portable products. There are some previous works designed with the architecture of low drop out(LDO) voltage regulator [START_REF] Rincon-Mora | Current Efficient, Low Voltage, Low Drop Out Regulators[END_REF]. Meanwhile, it also provides a reversible polarity charging mechanism. The proposed charger provides applications greatly and conveniently. In this paper, the system structure of the proposed circuit design is shown in Section II. Section III provides the simulation results to the performance and specification of the proposed charger. Finally, the conclusion is made in section IV.

CIRCUIT DESIGN

This paper can improve the [START_REF]Multi-Mode LDO-Based Li-Ion Battery Charger in 0.35μm CMOS Technology[END_REF], [START_REF] Hwang | New compact CMOS Li-Ion battery charger using charge-pump technique for portable applications[END_REF]. The structure of the proposed charger includes mode selector, current generating circuit, current sensing circuit, and automatic voltage detecting circuit. Function block diagram is showed in Fig. 2.

SIMULATION RESULTS

The simulation results of the proposed charger using Cadence Spectre is shown in Fig. 3.

CONCLUSION

A Li-Ion battery charger with multimode techniques is proposed in this paper. When charger charges battery, the proposed techniques solve battery polarity problems. Beside, utilizing variable current source to control three-stage charging mode is safe. The maximum efficiency of the proposed charger is about 80%. The charger can provide the output voltage from 2.2V to 4.2V precisely. The proposed charger provides applications greatly and conveniently.

Fig. 1 (

 1 a) is the architecture of traditional Li-Ion battery charger. The drawbacks of battery needs polarity of electric charging, this paper proposed a method to implement the problem. The architecture of the proposed Li-Ion battery charger shows in Fig.1(b). The charger solves the device characteristics problems.

Fig. 1 (

 1 Fig. 1(a) Architecture of traditional Li-Ion battery charger.

Fig. 1 (

 1 Fig.1(b) Architecture of the proposed Li-Ion battery charger.

Fig. 2

 2 Fig. 2 Function block diagram of the proposed charge

Fig. 3

 3 Fig. 3 Battery voltage and charging current.

Fig. 4

 4 Fig.4shows layout graph of the proposed charger, the layout area is 0.5mm*1mm (with PADs).

Fig. 4 .

 4 Fig. 4. Layout graph of the proposed charger.

BIBLIOGRAPHIE


Ikram Daanoune

Abdennaceur Baghdad
email: nasserbaghdad@yahoo.fr

Abdelhakim Balllouk
email: balloukhakim@gmail.com

Definition of an energy optimization protocol of a wireless sensor network

Keywords: WSN, Clustering, LEACH, residual energy, network lifetime, distance

Nowadays the WSN (Wireless Sensor Network) implemented in different domains such as environment, medical, military, etc. At the same time, the WSN has some limitations, among them and the major one is the energy resource and limitation in the battery's lifetime due to the fact of the small size of the sensor node. The clustering is the most efficient method for increasing the network's lifetime which most routing protocols based on it and LEACH (Low Energy Adaptive Clustering Hierarchy) is the first energy-efficient clustering hierarchical routing protocol. This clustering programs choose randomly the CH (Cluster Head) and just if the sensor node has a few energy. Thus in LEACH transmission, every CH sends data immediately to the BS (Base Station), so they consume more energy that decreases the network lifetime. In this paper, we propose a new approach based on LEACH protocol that purpose to extend the network's lifetime by balancing the energy of sensor nodes. Our model allows to chose CHs according to the current energy. Thus, data archives the BS by multi-hop method.

I. INTRODUCTION

Recently, WSN becomes an emerging topic in scientific research and it becomes increasingly omnipresent. It is an ad hoc network of hundreds of micro-sensors deployed in large area any-place in the environment. This micro-sensors include four units, detection data from the environment, processing unit, wireless communication, and power unit (usually a battery). The small size of the sensor node battery imposes a major problem in WSN that influences the network lifetime. To solve this constraint, routing protocols are one of the most solutions have appeared. LEACH is one of the famous energy-efficient hierarchical routing protocol, that divided network into numerous clusters which collect data from cluster members, aggregate them, then send to the BS [START_REF] Fu | An energy balanced algorithm of LEACH protocol in WSN[END_REF]. In this paper, we first present a synthesis of the LEACH protocol and its improved algorithms, then we propose a model of our new protocol based mainly on the residual energy, distance to the BS and multi-hop method. This document is orderly as follows: Part II presents a review of past researches. The proposed model was introduced in part III. Part IV discusses a discussion between the proposed protocol and the present one. Finally, Part V concludes the article.

II. REVIEW OF PAST RESEARCHES

A. LEACH (Low Energy Adaptive Clustering Hierarchy)

Heinzelman et al. [START_REF] Heinzelman | Energyefficient communication protocol for wireless microsensor networks[END_REF] proposed a Low Energy Adaptive Clustering Hierarchy (LEACH) for WSN. This protocol provides a conception of clusters. In every cluster, it elected a sensor node as a CH and the rest of cluster nodes become normal nodes which gather data from the environment then send it to the CH, the CH is responsible for collecting information from other cluster members, compressing it and sending it to the BS [START_REF] Fu | An energy balanced algorithm of LEACH protocol in WSN[END_REF]. Each round in LEACH is characterized by two states: cluster set-up phase and steady phase. The decision to be a CH in the set-up phase comes back to the node by choosing a random value between 0 and 1. If the value is less than a threshold value T (n) defined by equation ( 1) this node will be a CH in the next round, else it turns into a normal node.

T (n) = P 1-P * (r mod(1/P )) , n ∈ G 0 otherwise (1)
Where P introduces the percentage of the CHs in the network, the number of the current round is presented by r; G presents the collections of the nodes that have not yet been CH in the last 1/p rounds. Fig. 1 presents the architecture of LEACH protocol. Fig. 2 presents the energy consumption of the radio module to send an s bit to remote receiver d meter.

Where E elec is the energy consumption of the transmitter to transmit or to receive the signal and E amp is the energy consumption of the amplifier.

In clustering routing protocols, the CH consumes more energy than normal nodes due to the fact that it collects the information, merges it and sends to the BS which is too far away.

So it has to choose the CHs, also LEACH algorithm is based on the single hop process between CHs and BS no matter the distance and the current energy of these CHs, which influences the network lifetime.

B. V-LEACH [3]

V-LEACH presents an enhanced algorithm of the hierarchical LEACH protocol that draws a new organization where the normal node gather a physical quantity from the environment and transmits it to the CH, this last aggregated the data received from cluster members and forwards it to the BS and vice-CH is the component that is turn out to be a CH in case of death of the primary CH [START_REF] Ahlawat | An extended vice-cluster selection approach to improve v leach protocol in WSN[END_REF] [5] [START_REF]A Literature Survey of Energy-Efficient clustering in WSN using Different Variant of LEACH Routing Protocol[END_REF]. The creation of a vice-CH in WSN enhance the LEACH algorithm, but the adoption of the single hop and that vice-CH only replaces the CH after its death and does not replace the other nodes impact on the network lifetime.

C. E-LEACH [7]

E-LEACH protocol is used to improve the CH selection algorithm. The main objective of this algorithm is to avoid the nodes with less residual energy while keeping nodes with greater current energy. E-LEACH allows picking a root sensor node with the highest residual energy. [START_REF] Xu | Improvement of LEACH protocol for WSN[END_REF]. The E-LEACH protocol is based on the residual energy to determine the optimal way to the BS, although it does not take into consideration the distance to the BS.

D. Energy-LEACH and Multihop-LEACH

Xiangning et al. [START_REF] Xiangning | Improvement on LEACH protocol of wireless sensor network[END_REF] have proposed two improved protocols of LEACH as Energy-LEACH and Multihop-LEACH. Energy-LEACH enhances the method of selection of CH, it takes into account the residual energy after the first round. Each node calculates its residual energy, some nodes with more residual energy are transformed into CHs while the rest transforms into normal nodes [START_REF] Tyagi | A systematic review on clustering and routing techniques based upon LEACH protocol for wireless sensor networks[END_REF]. The multihop-LEACH protocol has been proposed to enhance the communication between the CH and the BS by selecting the nearest CH at the BS, that will aggregate the data of the rest of CHs [START_REF] Tyagi | A systematic review on clustering and routing techniques based upon LEACH protocol for wireless sensor networks[END_REF]. Energy-LEACH protocol improves the procedure for selecting CHs, it takes into account the residual energy of the node, but it does not take into account the distance between the CHs and the BS, each CH communicates with this last any distance and based on the single-hop which increases the energy consumption. Contrary to Multihop protocol which uses multi-hop process between CHs, but this protocol doesn't take the residual energy, the root CH may have low residual energy.

E. LEACH-TLCH(LEACH Protocol with Two Levels Cluster Head)

The LEACH-TLCH protocol serves to elect a node with maximum energy in each cluster as a secondary CH at the next round if the mainly CH has residual energy less than average energy or distance to the BS is more than average distance. In that case, the secondary CH collects data from cluster members, aggregates it then sends it to the CH. In contrary, the CH is responsible only to send the information received to the BS [START_REF] Fu | An energy balanced algorithm of LEACH protocol in WSN[END_REF]. On the other hand, while the current energy of the CH is more than the average energy and the distance to the BS is less than the average distance, so this cluster does not need to select a secondary CH [START_REF] Fu | An energy balanced algorithm of LEACH protocol in WSN[END_REF]. One advantage of the LEACH-TLCH protocol is that it deals with the problem of residual energy of the CH. Another advantage of LEACH-TLCH protocol is that taking the distance to the BS into consideration. However, randomly cluster formation and reliance on single-hop process impacts the network's lifetime.

III. THE PROPOSED MODEL

Some nodes in WSN deplete their energies previously than others due to dynamically and random picking of CH, because the CH gather numerous packets from cluster members and aggregates them then sends to the BS. So it consumes more energy than a normal node. Thus communication between CH and the BS by single-hop no matter the distance and the residual energy of CH impact on network lifetime. Therefore, it is important to the choice optimal path to reach the BS. In this paper, we propose an enhanced model based on LEACH and some improved protocols. It is based on three metrics: residual energy, distance and multi-hop to achieve the BS. Our model serves to select CHs according to current energy. We select a number of nodes with more residual energy to become CHs in the next round, to avoid low residual energy nodes turning into CHs and to give nodes with more energy residual to take part as CH.

Fig. 1 .

 1 Fig. 1. LEACH protocol architecture [2].

Fig. 2 .

 2 Fig. 2. Radio model.

A. Set-up phase CHs are selected according to the current energy. Cluster formation ends when all nodes in the network join their corresponding CHs, basing on the received signal from CHs.

B. Steady phase

TDMA (Time Division Multiple Access) schedules are created by each CH depending on the number of cluster members and the distance. Every node sends its data in its TDMA time slot to its CH [START_REF] Gowda | DUCA: An Approach to Elongate the Lifetime of Wireless Sensor Nodes[END_REF]. Afterwards, we pick the CH which more residual energy and less distance to the BS as a root CH. When clusters, the TDMA schedule, and the root CH are formed, the normal nodes starting to send the collected data to the CH. The farthest CHs achieve the root CH using multi-hop process. The root CH aggregates all data and by single-hop method sends it to the BS. The flowchart Fig. 3 shows the steps of the proposed model.

IV. DISCUSSION

The proposed algorithm is based on three routing metrics: residual energy, distance, and multi-hop. The proposed approach intended to improve network lifetime and balance power consumption based on LEACH and its improved protocols. Table.I presents a relative comparison of LEACH based protocols according to data transmission (single hop to BS or multi-hop), routing protocol, residual energy, and the distance to the BS.

V. CONCLUSION

The main goal of the current study was to present a comparative synthesis of the performances of the original hierarchical routing LEACH protocol and some of their improved algorithms such as E-LEACH, V-LEACH, Energy-LEACH, Multihop-LEACH, and LEACH-TLCH. The second aim of this study was to define a new improved LEACH protocol coming to extend the network lifetime, based mainly on residual energy, distance to the BS and multi-hop process. According to our hypothesis, this model is outperforming other protocols discussed because it is better suited to the constraint of energy consumption in the WSN.