domain: spi phys

text:
Laurent Muscat

Guillaume Puigt

Marc Montagnac

Pierre Brenner

A coupled implicit-explicit time integration method for compressible unsteady flows

This paper addresses how two time integration schemes, the Heun's scheme for explicit time integration and the second-order Crank-Nicolson scheme for implicit time integration, can be coupled spatially. This coupling is the prerequisite to perform a coupled Large Eddy Simulation / Reynolds Averaged Navier-Stokes computation in an industrial context, using the implicit time procedure for the boundary layer (RANS) and the explicit time integration procedure in the LES region. The coupling procedure is designed in order to switch from explicit to implicit time integrations as fast as possible, while maintaining stability. After introducing the different schemes, the paper presents the initial coupling procedure adapted from a published reference and shows that it can amplify some numerical waves. An alternative procedure, studied in a coupled time/space framework, is shown to be stable and with spectral properties in agreement with the requirements of industrial applications. The coupling technique is validated with standard test cases, ranging from one-dimensional to three-dimensional flows.

. However, explicit schemes suffer from small time steps that are limited by the Courant-Friedrichs-Lewy (CFL) number. Hence, computing boundary layers with a fine mesh at the wall with LES becomes very demanding from the CPU time point of view.

. Therefore, for many stiff problems, the use of explicit methods requires very small time steps for the calculation to be feasible or stable.

.

 .

. The remaining question is how to couple this second-order explicit scheme with a second-order implicit scheme, making a hybrid version of the underlying schemes that still remains second-order accurate.

 and the analysis performed by Ascher et al. [31]). However, in our case, the coupling of RANS and LES leads to spatially couple explicit and implicit schemes as implicit schemes are suitable for RANS equations

Introduction

Turbulence occurring in an unsteady compressible fluid flow is characterized by the cascade of vortices distributed over a large range of temporal and spatial scales [START_REF] Chapman | Computational Aerodynamics Development and Outlook[END_REF]. When the Reynolds number Re is large (typically for Re 10 6 -10 9 ), the resolution of the entire turbulence spectrum using Direct Numerical Simulation is still out of reach today for industrial configurations, even on the most powerful supercomputers. Currently, two major classes of approaches are mainly used to model turbulence effects in these situations.

The first technique consists in modeling all effects of the turbulence on the flow and leads to Reynolds-Averaged Navier-Stokes (RANS) equations. All turbulence effects are represented by their averaged effects on the mean flow, and turbulence fluctuations are accounted for through a turbulence model. The RANS approach is today the preferred technique for industrial applications, even if it is unable to represent unsteady effects of the turbulence. In particular, its accuracy strongly depends on the turbulence model used to represent the turbulence effects on the mean flow. The main interest of this approach remains its relatively low CPU cost for reasonable accuracy, especially for computing boundary layers. Implicit time integrators are generally considered in that case since they allow either fast convergence to steady-state solutions or large time steps for unsteady flows. Indeed, implicit schemes are generally chosen so to be unconditionally stable. These schemes generally need the linearization of the implicit residual, which makes the Jacobian of the residual appear. Finally, this leads to look for the solution increment as the solution of a linear system. Examples of implicit time integrators are the family of implicit Runge-Kutta methods (IRK) proposed by Kunstmann [START_REF] Kuntzmann | Neure Entwicklungen der Methoden von Runge und Kutta[END_REF] and Butcher [START_REF] Butcher | Implicit Runge-Kutta processes[END_REF] or the Gear's scheme [START_REF] Gear | Numerical Initial Value Problems in Ordinary Differential Equations[END_REF]. Implicit methods have the drawback to require extra computational costs compared with explicit techniques, because a linear system must be solved at all time steps.

Another way to account for turbulence effects is Large Eddy Simulation (LES). Indeed, the principle of LES is to compute the largest turbulent scales and to model the smallest scales that cannot be captured by the mesh. Of course, the LES approach introduces additional terms to account for the intractable turbulence spectrum but the universal nature of the smallest scales makes the models more general than for the RANS approach. In practice, LES requires a very fine mesh to capture turbulence and explicit time integration for LES. To our knowledge, much less attention was paid on techniques to couple explicit and implicit schemes spatially. At this point, it must be highlighted that for practical applications, the monitoring of the convergence of the implicit formulation has a strong influence on the quality of the results , especially for LES. This is a reason why a full implicit formulation for both RANS and LES areas is not retained.

In 1986, Fryxell et al. [START_REF] Fryxell | An implicit-explicit hybrid method for lagrangian hydrodynamics[END_REF] built an hybrid implicit-explicit method to perform one-dimensional Lagrangian hydrodynamics simulation based on the extension of a Godunov-type scheme to implicit regime. Dai and Woodward [START_REF] Dai | A second-order iterative implicit-explicit hybrid scheme for hyperbolic systems of conservation laws[END_REF] proposed an iterative approach of this framework for Euler equations. Their scheme is second-order accurate in time and space but quite complicated to implement. This is especially a consequence of the requirement of solving the solution by an iterative process. In addition, the authors mention that the extension to multi-dimensional system is still an open question. Collins et al. [START_REF] Collins | An Implicit-Explicit Eulerian Godunov Scheme for Compressible Flow[END_REF] proposed the hybrid explicit-implicit scheme for Eulerian hydrodynamics, denoted CCG, that ensures the total variation diminishing (TVD) property at all CFL numbers for a linear equation. This hybrid scheme allows a second order of accuracy in space and time on pure explicit mesh cells. Unfortunately, the CCG scheme leads to incorrect numerical solutions and fails to maintain the TVD property in case of nonlinear equations. Men'shov and Nakamura [START_REF] Men'shov | Hybrid Explicit-Implicit, Unconditionally Stable Scheme for Unsteady Compressible Flows[END_REF] proposed an extension of this formulation to non-linear hyperbolic equation with Maximal Norm Diminishing (MND) property. There still remains a strong limitation: the switching procedure of CCG and Men'shov and Nakamura schemes (from explicit to implicit and vice versa) cannot maintain a second order of accuracy in space and time over the whole mesh. Timofeev and Norouzi [START_REF] Timofeev | Hybrid, explicit-implicit, finite-volume schemes on unstructured grids for unsteady compressible flows[END_REF] overcame this difficulty by blending second-order accurate schemes (in space and time) for explicit, implicit and hybrid cells using a smooth parameter to couple the considered time integrators. Their method is one of the most advanced ones since it maintains TVD property. In 2006, second order accurate explicit and implicit time integrators were hybridized by Tóth et al. [START_REF] Tóth | A parallel explicit/implicit time stepping scheme on block-adaptive grids[END_REF] and coupled with an adaptive mesh refinement (AMR) strategy for three-dimensional applications in magnetohydrodynamics using structured grids. The procedure is quite complex to implement and needs to blend the first and second order fluxes on the interfaces.

In 2017, May and Berger [START_REF] May | An Explicit Implicit Scheme for Cut Cells in Embedded Boundary Meshes[END_REF] proposed a technique to switch between explicit and implicit approaches called flux bounding. The coupling of explicit and implicit schemes is performed by the coupling of explicit and implicit flux contributions on the interface in the context of finite volume discretization procedure. Unfortunately, it led to a result that is not second-order accurate in time and space for "transient" cells for which residuals were computed using explicit and implicit flux contributions. Finally, it is important to notice that Persson [START_REF] Persson | High-Order LES Simulations using Implicit-Explicit Runge-Kutta Schemes[END_REF] also provided a Runge-Kutta based IMEX scheme that allows to couple spatially explicit and implicit time integration schemes. However, this technique does not respect our requirement to focus on Heun's scheme i.e. the explicit time integrator in FLUSEPA 1 . In this paper, a new way to couple spatially second-order accurate explicit and implicit time integrators for compressible unsteady flows is introduced. The goal is to propose a coupling technique that ends to a second-order accurate scheme over the whole computational domain. The remainder of this paper unfolds as follows. In Sec. 2, the standard form of the Navier-Stokes is recalled and the explicit and implicit schemes chosen for the present study are introduced. Our explicit and implicit basic schemes, Heun and implicit Runge Kutta schemes, are defined in Sec. 3. In Sec. 4, the method presented in [START_REF] Timofeev | Hybrid, explicit-implicit, finite-volume schemes on unstructured grids for unsteady compressible flows[END_REF] is recalled as it will serve as the basis of the new scheme. In Sec. 5, a first version of coupling Heun and implicit Runge Kutta (IRK2) time integrators is introduced with its dedicated space-time stability analysis. The first hybridized scheme is not acceptable and a second version of the coupling scheme is presented in Sec. 6. The space-time stability analysis is performed again. Before concluding, Sec. 7 is dedicated to the validation of the new coupling scheme for 1D, 2D and 3D test cases.

Discretization of the Navier-Stokes Equations

The Navier-Stokes conservation laws can be written in the following compact form,

∂W ∂t + ∇ • F (W ) + ∇ • G(W, ∇W ) = 0, ( 1
)
where W is the vector of conservative variables and ∇W its gradient, F and G are respectively the flux density for convection and viscous effects. The computational domain Ω is split into N non-overlapping rigid immobile cells Ω j and Eq. ( 1) are integrated over every mesh cell. The Gauss relation ties the volume integrals of the divergence terms to the interface fluxes, leading to d dt

Ω j W dΩ = - A j F (W ) • n dS - A j G(W, ∇W ) • n dS, (2)
where Ω j is the j-th control volume with its border A j and n is the outgoing unit normal vector. From now on, it is possible to define the averaged quantity of the conservative variables,

W j = 1 |Ω j | Ω j W dΩ. (3)
This formula allows to rewrite the conservation laws (2) discretized with a finite-volume formulation in the following compact form,

dW j dt = R(W j ), (4)
with R(W j ) is the residual computed using the averaged quantities in cell j. Actually, the residual is defined as the sum of the flux over the whole boundary of a cell. Here, convection and diffusion fluxes are discretized by means of a k-exact formulation coupled with successive corrections [START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible RANS-LES simulations on fully unstructured grids[END_REF][START_REF] Haider | Parallel Implementation of k-Exact Finite Volume Reconstruction on Unstructured Grids[END_REF]. Indeed, any high-order scheme needs a polynomial representation of the unknowns locally around a cell of interest. Theoretically, the polynomial approximation can also be defined as the Taylor expansion of the local solution. Standard k-exact schemes need the definition of a large stencil in order to link Taylor expansion coefficients and averaged quantities (see [START_REF] Ollivier-Gooch | A High-Order-Accurate Unstructured Mesh Finite-Volume Scheme for the Advection-Diffusion Equation[END_REF] as an example of standard k-exact scheme). This large stencil is a drawback to cope with parallel computations since the number of data to exchange between domains depends on the partitioned (initial) mesh. The successive corrections were designed in order to avoid geometrical reconstructions for parallel computations. The number of ghost cells are limited to the minimum of one ghost cell per volume owing an interface on the boundary. Starting from the solution, the polynomial reconstruction is built using the truncation error of the Taylor expansion by increasing the order of the approximation. The starting point is the availability of the averaged quantities in the mesh cells and the definition of a discrete gradient scheme. The discrete gradient scheme applied to the averaged quantities leads to a given order of accuracy. A local analysis enables to estimate the truncation error performed on the gradient and by removing the last non-zero term in the truncation error, the gradient scheme accuracy is made more accurate by one order. This gradient is then usable to define a more accurate representation of the unknown through a local polynomial reconstruction. The process can then be repeated for the Hessian in order to increase accuracy on the computation of the successive derivatives [START_REF] Haider | Parallel Implementation of k-Exact Finite Volume Reconstruction on Unstructured Grids[END_REF][START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible RANS-LES simulations on fully unstructured grids[END_REF]. Of course, since the truncation error must be computed, it is necessary to locally exchange the set of local derivatives (first, the gradients, then the second-order derivatives, then the third-order derivatives and so on). The TVD property is ensured with a MUSCL reconstruction. For a second-order scheme, the flux on an interface is estimated using an affine extrapolation of the cell-centered unknowns onto the interface. If f denotes a mesh interface with a unit normal vector n directed from the left cell j to the right cell i, and if C i , C j and C f respectively denote the cell centers and the interface center, the standard reconstruction on the vector of primitive variables V , is simply,

V L =V n j + (∇V ) n j • ---→ C j C f , V R =V n i + (∇V ) n i • ---→ C i C f , (5)
where the gradient ∇V needs to be computed with the k-exact reconstruction. These gradients are also used to compute the viscous flux. Finally, the finite-volume formulation of Eq. ( 4) can be expressed as the following Cauchy problem,

   dW j (t) dt = R(W j ), ∀t ∈ R + , W j (0) = W 0 j , ∀j ∈ {1, ..., N }, (6)
where W 0 j 1≤j≤N denotes the initial solution. For the sake of clarity, the averaging symbol will be dropped, and W will represent the averaged quantities over the control volumes.

Heun and Crank-Nicolson/IRK2 Time Integration Methods

Heun and Crank-Nicolson/IRK2 are considered as the two basic time integrators of the present study from which a new method will be developed. The Heun's explicit scheme [START_REF] Heun | Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen veränderlichen[END_REF] is a second-order accurate predictor-corrector method. The state W n+1 is first predicted with a forward Euler scheme, and then corrected by a standard trapezoidal rule,   

Predictor step: W = W n + ∆tR(W n ),
Corrector step:

W n+1 = W n + ∆t 2 R(W n ) + R( W ) . (7)
The Crank-Nicolson/IRK2 implicit scheme [START_REF] Crank | A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type[END_REF] is second-order accurate, and it can be expressed as,

W n+1 = W n + ∆t 2 R(W n ) + R(W n+1 ) . (8)
The residual being non linear in W , Eq. ( 8) is solved with an iterative Newton method. Starting from the iteration n in Eq. ( 8), the next solution is defined as the root X of f (X) = 0 with

f (X) = X -W n - ∆t 2 R(X) + R(W n ) . (9)
At the p-th step of the Newton algorithm, the next iterative X p+1 is such that

f (X) | X p (X p+1 -X p ) = -f (X p ). (10)
The final algorithm is simply

         X 0 = W n (I - ∆t 2 J) δX p = -f (X p ) X p+1 = δX p + X p (11)
with J = ∂R ∂W the Jacobian matrix and I the identity matrix. The linear system in the iterative loop is solved with a Krylov iterative method. The proposed Newton procedure is motivated by the strong non-linearities encountered during reentry and take-off. For small non-linear effects, the residual R(W n+1 ) can be linearized around

W n .
The open question concerns how to spatially couple both second-order time integrator schemes, while maintaining the second order of accuracy and ensuring the overall stability, and taking into account the fact that both time integrators need two residuals computed at two different times. To this aim, the hybrid time integrator proposed by Timofeev and Norouzi [START_REF] Timofeev | Hybrid, explicit-implicit, finite-volume schemes on unstructured grids for unsteady compressible flows[END_REF] is first considered. In the following, all numerical schemes will be expressed in a one-dimensional domain without loss of generality, and for convective terms. Diffusive terms will be discussed if required.

Principle of the Hybrid Time Integration

The spatially-coupled explicit and implicit scheme proposed by Men'shov and Nakamura [START_REF] Men'shov | Hybrid Explicit-Implicit, Unconditionally Stable Scheme for Unsteady Compressible Flows[END_REF] and extended by Timofeev and Norouzi [START_REF] Timofeev | Hybrid, explicit-implicit, finite-volume schemes on unstructured grids for unsteady compressible flows[END_REF] is recalled hereafter to make an introduction to the hybrid time integration. The hybrid scheme [START_REF] Timofeev | Hybrid, explicit-implicit, finite-volume schemes on unstructured grids for unsteady compressible flows[END_REF] applied to the Cauchy problem of Eq. ( 4) has a predictorcorrector formulation:

     Predictor step: W j = W n j + ω j ∆tR(W n j )
Corrector step:

W n+1 j = W n j + ∆t |Ω j | (F Hybrid j+ 1 2 -F Hybrid j-1 2 ) for 0 ≤ ω j ≤ 1. ( 12
)
The main components of the scheme are presented in the following using a face f , with the face center C f , separating two cells, i and j, with respectively cell centers C i and C j .

During the predictor step, the residual computation does not use any Riemann solver, and the predicted state W j is computed using a fully upwind scheme. So, the flux F (V f j ), computed from the vector of the primitive variables V and involved in the computation of the residual R(W n j ) of the cell j uses the reconstruction on the primitive state [START_REF] Norouzi | A hybrid, explicit-implicit, second order in space and time TVDscheme for one-dimensional scalar conservation laws[END_REF][START_REF] Timofeev | Application of a new hybrid explicit-implicit flow solver to 1d unsteady flows with shock waves[END_REF][START_REF] Norouzi | A hybrid,explicit-implicit,second order in space and time TVD scheme for two-dimensional compressible flows[END_REF]. This prediction step is not conservative since the flux F (V f j ) = F (V f i ). Moreover, the gradient ∇V is computed from the standard gradient ∇V and a slope limiter.

V f j = V n j + ( ∇V ) n j • ---→ C j C f
For any cell j, the hybrid flux F Hybrid is computed using a Riemann solver with the following reconstructed left and right primitive states on a face:

V L = ω j V n j + V j 2 + ( ∇V ) n j • ---→ C j C f + (1 -ω j ) V n+1 j + ( ∇V ) n+1 j • ---→ C j C f - 1 -ω j 2 ( ∆ t V ) n+1 j , V R = ω i V n i + V i 2 + ( ∇V ) n i • ---→ C i C f + (1 -ω i ) V n+1 i + ( ∇V ) n+1 i • ---→ C i C f - 1 -ω i 2 ( ∆ t V ) n+1 i . (13)
The parameter ω j , subsequently named cell status, is defined for each cell j such that 0 ≤ ω j ≤ 1. It is used to switch smoothly from a time-explicit (ω j = 1) to a time-implicit (ω j = 0) scheme. The cell status ω j is locally adapted to the flow and the stability constraint. First, the global time step is chosen to be the maximum allowable time step on a part D of the whole computational domain without source of stiff terms: ∆t = min

j∈D ν h j v j + c j , (14)
where ν < 1, h j is a reference length scale, v j the velocity vector, and c j the speed of sound in cell j.

Then, the parameter ω j is defined as:

ω j = min 1, 1 ν j , with ν j = ( v j + c j ) ∆t h j . (15)
In Eq. ( 13), the term ( ∆ t V )

n+1 j
comes from a generalization of the TVD property to the coupled space/time reconstruction. It can be seen as a time slope limiter. This term is given by a Newton's algorithm which at any step s reads:

( ∆ t V ) n+1,s j =                    max 0, min min j [β(V n+1,s j -V n j ) + ( ∆ t V ) n+1,s-1 j ], V n+1,s j -V n j if V n+1,s j -V n j ≥ 0, min 0, max max j [β(V n+1,s j -V n j ) + ( ∆ t V ) n+1,s-1 j ], V n+1,s j -V n j if V n+1,s j -V n j < 0. ( 16
)
The parameter β is computed from the local CFL number ν j and ω j at each step s:

β = 2 1 -ω s j ν s j (2 -ω s j ν s j ) ν 2 j (1 -ω s i ) 2 (17)
This method is second-order accurate in both space and time in fully explicit, fully implicit and explicit-implicit domains. Moreover, the scheme was proved to be TVD for 1D linear/non linear scalar conservation laws. It was applied successfully to problems including shocks [START_REF] Timofeev | Hybrid, explicit-implicit, finite-volume schemes on unstructured grids for unsteady compressible flows[END_REF]. This scheme will be referred as TN in the following.

It is clear that the underlying explicit (where ω j = 1) and implicit (where ω j = 0) schemes are not Heun's and Crank-Nicolson/IRK2 time integrators. In section 5, a first straightforward attempt to couple these latter schemes is explained.

First Hybrid Coupled Scheme

This section presents the coupling of three schemes: Heun for the explicit part, IRK2 for the implicit part, and the hybrid formulation TN [START_REF] Timofeev | Hybrid, explicit-implicit, finite-volume schemes on unstructured grids for unsteady compressible flows[END_REF] for a smooth transition between the latter two schemes. The resulting numerical time integration scheme is defined to ensure a unique flux on each interface, leading to a conservative formulation.

An interface separates two cells that may have a different status given by the parameter ω j . The cell status is named Heun for ω j = 1, IRK2 for ω j = 0, and Hybrid otherwise, as shown by the 1D example on Fig. 1. The flux formula used on an interface separating two cells depends on the status of these cells as shown in Tab. 1.

Table 1: Flux formula depending on neighbour cells. The table is "symmetrical": the flux between two cells is independent of the direction of information propagation.

Left cell status

Right cell status Heun Hybrid IRK2

Heun

F Heun F Heun × Hybrid F Heun F Hybrid F IRK2 IRK2 × F IRK2 F IRK2
According to Tab. 1, for the 1D example on Fig. 1, the first version of the coupling scheme, referred as HCS1, is designed as,

                                     Predictor step: W j = W n j + ω j ∆tR(W n j )
Corrector step:

                               W n+1 j-1 = W n j-1 + ∆t 2|Ω j-1 | (F n j-1 2 + F j-1 2 -F n j-3 2 -F j-3 2 ) W n+1 j = W n j + ∆t |Ω j | (F Hybrid j+ 1 2 - F n j-1 2 + F j-1 2
2 )

W n+1 j+1 = W n j+1 + ∆t |Ω j+1 | ( F n j+ 3 2 + F n+1 j+ 3 2 2 -F Hybrid j+ 1 2 ) W n+1 j+2 = W n j+2 + ∆t 2|Ω j+2 | (F n j+ 5 2 + F n+1 j+ 5 2 -F n j+ 3 2 -F n+1 j+ 3 2
).

(

) 18
The hybrid flux F Hybrid is computed with the state reconstruction method introduced in Eq. ( 13). The scheme Eq. ( 18) ensures the unique definition of the interface flux, and then flux conservation.

The differences with the TN scheme are:

• In the predictor step, the residual is computed from the flux balance over all cell faces. The flux is computed using a linear extrapolation of the unknowns between cell centers and cell interfaces and accounting for a slope limiter. Since two different extrapolations are obtained on a face shared by two cells, the flux is computed by means of a (approximated) Riemann solver. As a consequence, W j has a physical meaning, something which was lost with TN scheme.

• The scheme recovers the Heun's scheme in the explicit regime (ω j = 1).

• The scheme recovers the IRK2 scheme in the implicit regime (ω j = 0).

Remark 1: It is very important to make the transition between explicit and implicit time integration as fast as possible. This is the reason why in practice, the IRK2 formulation is applied when ω j ≤ 0.5.

Remark 2: For the Navier-Stokes equations, the standard diffusion scheme available in the code was used for the implicit and explicit regions. The interface gradient is defined from the limited cellcentered gradients by means of an average taking into account local metrics. In the hybrid regions, the definition of the cell-centered gradient is modified as

(∇V ) hyb j = ω j ( ∇V ) n j + ( ∇V ) j 2 + (1 -ω j )( ∇V ) n+1 j for 0.5 < ω j < 1. ( 19
)
Every gradient in Eq. ( 19) is computed from the standard gradient and corrected using a slope limiter.

Space-Time von Neumann Analysis

The von Neumann analysis is a powerful tool to analyse the spectral properties of any scheme. This analysis is generally performed separately either for spatial or temporal schemes but it was shown in many recent papers (see [START_REF] Sengupta | Spurious waves in discrete computation of wave phenomena and flow problems[END_REF][START_REF] Sengupta | Space-time discretizing optimal DRP schemes for flow and wave propagation problems[END_REF][START_REF] Sengupta | Analysis of central and upwind compact schemes[END_REF][START_REF] Vanharen | Revisiting the spectral analysis for high-order spectral discontinuous methods[END_REF] for instance) that the true scheme behaviour is only obtained by a coupled analysis involving both space and time discretizations.

Starting from the one-dimensional linear advection equation coupled with a harmonic initial condition and with periodic boundary conditions,

         ∂y ∂t + c ∂y ∂x = 0 y(x, 0) = exp(ikx) x ∈ [0, L] y(0, t) = y(L, t) t ∈ R + (20)
with i 2 = -1 and c > 0, the comparison between the exact theoretical solution and the numerical approximation gives information on both dissipation and dispersion of the fully discrete scheme. The fully discrete relation obtained from Eq. ( 20) using a second-order finite volume scheme reads:

y n+1 j = G j y n j , (21)
where y n j represents the harmonic solution at discrete time n and in the center of cell j. The complex coefficient G j represents the transfer function between two consecutive time solutions. The coefficient G j can be expressed as G j = |G j | exp(i arg(G j )), which highlights the dissipation of the scheme related to |G j |, and its dispersion related to exp(i arg(G j )). In the following, let µ j = |G j | be the dissipation coefficient, and φ j = arg(G j )/CFL the dispersion coefficient.

For standard finite element, difference or volume formulations, and focusing on a single scheme for space and time, the transfer coefficient does not depend on the space position. However, in our case, several time discretization schemes are blended, and the spectral properties of the global scheme will depend on the cell status ω j . Accounting for the cell status variations can be done through the change in the mesh spacing. But taking into account for a local change in space discretization in the spectral analysis is cumbersome (see for instance [START_REF] Vichnevetsky | Energy and group velocity in semin discretizations of hyperbolic equations[END_REF]). Once the mesh size ∆x is fixed, a manufactured choice of ω j keeps the analysis quite simple to perform. The number of 1D segments is fixed to N = 300 and the function ω j is defined as

           ω j = α ω j-1 for N 2 -50 ≤ j ≤ N 2 ω j = 1 α ω j-1 for N 2 + 1 ≤ j ≤ N 2 + 50 ω j = 1 elsewhere. ( 22
)
with α = 0.90. The maximum CFL value for the Heun explicit scheme without amplification at any wavenumber is 0.6. Then, the spectral analysis will be performed for CFL=0.1 and CFL=0.6. Figures 2 and3 give an overview of the dissipation and the dispersion of the coupled time/space HCS1 scheme depending on the cell index 1 ≤ j ≤ N and the non dimensionalized wave number 0 ≤ k∆x ≤ π for CFL=0.1. They clearly show that some phenomena occur at the transitions between the different temporal schemes. As a consequence, the standard curves (µ j and φ j as a function of k∆x) are introduced in Figs. 456for the specific cells at which there is a time integration transition. As indicated in Tab. 2, six specific discretization cells are analysed. In addition, Tab. 2 also gives the meaning of the legends for Figs. 456. Figure 4 shows the dissipation coefficient. The focus for k∆x ∈ [0, 0.6] in Fig. 5 shows that the scheme amplifies some wavenumbers for Heun/Hybrid cells. The level of amplification is very low and one could assume that since the treatment is local (on one cell in 1D), the amplification could not lead to global instability and failure of the computation. Fig. 6 shows that cells with Hybrid/Heun and Hybrid/IRK2 fluxes produce the highest levels of dispersion, while dispersion is close to the reference Heun's scheme for other choices.

The same analysis is performed at CFL=0.6, i.e. at the limit without amplification at any wavenumber for the Heun's scheme. The dissipation and the dispersion curves are shown in Figs. 78. The dispersion curve for the Heun/Hybrid cell looks like the one for Heun's scheme and presents a change of phase sign for k∆x 1.8. The dissipation curve presents again amplification but the level of amplification is much higher. Since the value of ω j may be defined to follow some physical property, numerical amplification of a physical phenomena could follow the phenomena itself, leading to an unacceptable flow or simulation break. Notice that the discrete jumps on dispersion curves in Fig. 8 at k∆x 1.2, 1.6, 1.8, 2.2, 2.4, and 2.6. These values of k∆x are maxima of dispersion error, and correspond to values of wave number k∆x when Re(G j ) = 0.

The first coupled scheme HCS1 based on Heun and IRK2 schemes presents amplification. In order to remove it, a new way to couple the desired explicit and implicit schemes is proposed in Sec. 6.

Figure 7: Dissipation µ j of the scheme HCS1 for all kinds of cells at CFL=0. [START_REF] Sengupta | Space-time discretizing optimal DRP schemes for flow and wave propagation problems[END_REF] 6 Another tIme integratiON (AION) scheme

In order to get better spectral characteristics than the HCS1 scheme, a new hybrid scheme still based on both Heun and IRK2 schemes is designed. This new hybrid time integrator is called AION (Another tIme integratiON). It is applied to the Cauchy problem Eq. ( 4).

This hybrid scheme has a predictor-corrector formulation as TN's and Heun's schemes. According Figure 8: Dispersion φ j of the scheme HCS1 for all kinds of cells at CFL=0.6

to Tab. 1, for the 1D example on Fig. 1, it is designed as:

                                    
Predictor step:

W j = W n j + ∆tR(W n j )
Corrector step:

                               W n+1 j-1 = W n j-1 + ∆t 2|Ω j-1 | (F n j-1 2 + F j-1 2 -F n j-3 2 -F j-3 2 ) W n+1 j = W n j + ∆t |Ω j | (F Hybrid j+ 1 2 - F n j-1 2 + F j-1 2
2 )

W n+1 j+1 = W n j+1 + ∆t |Ω j+1 | ( F n j+ 3 2 + F n+1 j+ 3 2 2 -F Hybrid j+ 1 2 ) W n+1 j+2 = W n j+2 + ∆t 2|Ω j+2 | (F n j+ 5 2 + F n+1 j+ 5 2 -F n j+ 3 2 -F n+1 j+ 3 2 ). (23)
Remark: A situation not adressed in Tab. 1 nor in Fig. 1 must be described. For any hybrid cell (according to the value ω j ) with hybrid flux for all the faces, the predictor state W j is time-integrated by accounting for ω j : W j = W n j + ω j ∆tR(W n j ). For the mixed explicit / implicit regime, the hybrid flux differs due to the different reconstruction methods which is the key point to maintain stability without amplification. The demonstration of the stability without amplification is the purpose of the next section.

The reconstructed states are defined as:

V L = ω j V n j + V j 2 + 1 2 ( ∇V ) n j • ---→ C j C f + ω j - 1 2 ( ∇V ) j • ---→ C j C f + (1 -ω j ) V n+1 j + ( ∇V ) n+1 j • ---→ C j C f - 1 -ω j 2 ( ∆ t V ) n+1 j , V R = ω i V n i + V i 2 + 1 2 ( ∇V ) n i • ---→ C i C f + ω i - 1 2 ( ∇V ) i • ---→ C i C f + (1 -ω i ) V n+1 i + ( ∇V ) n+1 i • ---→ C i C f - 1 -ω i 2 ( ∆ t V ) n+1 i . ( 24
)
For ω j = 1 and ω i = 1, the following reconstruction is obtained:

V L = V n j + V j 2 + 1 2 ( ∇V ) n j • ---→ C j C f + 1 2 ( ∇V ) j • ---→ C j C f , V R = V n i + V i 2 + 1 2 ( ∇V ) n i • ---→ C i C f + 1 2 ( ∇V ) i • ---→ C i C f , (25)
and the hybrid part leads to the Heun's scheme for the one-dimension linear advection equation Eq. ( 20).

Space-Time Stability Analysis

The analysis introduced in Sec. 5.1 is now applied to the AION scheme. The spectral behaviour is shown in Figs. 9 and 10 for each cell index 1 ≤ j ≤ N and wave number k∆x ∈ [0, π]. Using the isocontours of iso-µ j and iso-φ j , it is clear that the spectral behaviour of the scheme depends on the cell index j via the cell status ω j . As before, the standard curves are now presented for the 6 values of j introduced in Tab. 2 at CFL=0.1 and CFL=0.6. At CFL=0.1, the overall view of the dissipation coefficient in Fig. 11 shows that the AION scheme seems to have the same behaviour for any cell index. In addition, a zoom for k∆x ∈ [0, 0.6] does not show any amplification of waves and stability is kept for any kind of cell (Fig. 12). Finally, all schemes present the same kind of dispersion in Fig. 13. This shows that the limitation of the initial HCS1 scheme is now removed by our new coupling technique at CFL=0.1. In addition, as shown in Figs. 141516, the same conclusion is obtained at CFL=0.6, which is the limit without amplification at any wavenumber for the pure explicit scheme.

For this example of linear advection, the proposed scheme is shown to be stable. The value of ω j is set to vary between 0 and 1. For ω j = 1, it was shown that the reconstructed states at the interface enable to recover the standard explicit Heun's scheme for the linear equation. So, the idea is to keep a smooth transition from explicit to hybrid regions. However, it would be convenient to switch to implicit regions as fast as possible. In the next section, attention is paid on the transition between implicit and hybrid cells.

Analysis of the Hyb/IRK2 and IRK2/Hyb Transitions

This analysis is motivated by the desire to switch as fast as possible to the implicit formulation, while maintaining stability. The computational domain with 300 cells is split into two parts. Half of the domain is dedicated to the hybrid part of the AION scheme while the rest of the domain is dedicated to the IRK2 scheme. The spectral analysis is performed at the transition between the chosen schemes. Two configurations are analysed, one for the wave going from the hybrid domain to the implicit domain, and one for the wave going from the implicit domain to the hybrid domain.

Figs. 17 and 19 show the dissipation coefficient for CFL=0.5 and for CFL=0.6. For the transition between implicit and hybrid schemes. Amplification (positive dissipation) is represented by a grey area. Amplification occurs for the transition Hyb/IRK2. In the reverse mode, amplification does not occur, as highlighted by Figs. 18 and20. Moreover, a stable formulation is obtained for ω j ≤ 0.6, and this motivates our final choice to switch at ω j = 0.6 between the hybrid flux computation and the implicit Runge-Kutta scheme. As a consequence and for the following applications, ω j = 1 will lead to the Heun's scheme, ω j < 0.6 will lead to the IRK2 scheme. The new hybrid reconstruction method based on the reconstructed states Eq. ( 24) is applied for 0.6 ≤ ω j ≤ 1.

Figure 17: Isocontours of the dissipation µ j for Hyb/IRK2 transition at CFL=0.5 as a function of ω j and k∆x. The grey area represents the area of amplification.

Analysis of q-waves

Focusing on the one-dimensional linear advection equation with a constant velocity c, several authors mentioned and explained the existence of numerical waves travelling in the direction opposite to c, introducing the notion of q-waves. Indeed, a p-wave is a wave travelling in the same direction as c, eventually not at the same speed. q-waves are travelling in the opposite direction. Existence of q-waves was shown to be independent of the stability analysis of numerical scheme and was highlighted for instance in the seminal work of Vichnevetsky, summarised in [START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF], by Poinsot and Veynante for combustion [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF] or by Trefethen [START_REF] Trefethen | Group velocity in finite difference schemes[END_REF]. In this context, Sengupta et al. [START_REF] Sengupta | Spurious waves in discrete computation of wave phenomena and flow problems[END_REF] linked the numerical group velocity with the existence of q-waves. Following the definition proposed in [START_REF] Sengupta | Spurious waves in discrete computation of wave phenomena and flow problems[END_REF], new quantities related to the one-dimensional linear advection equation Eq. ( 20) are introduced. The numerical phase speed The coupled approach is always stable. is defined as

c N j = arg(G j ) k ∆t , ( 26
)
from which is deduced the numerical group velocity,

V gN j = d arg(G j ) dk c CF L ∆x . (27)
Indeed, for a positive advection velocity c > 0, q-waves are propagating upstream and their group velocity V gN is negative [START_REF] Sengupta | Spurious waves in discrete computation of wave phenomena and flow problems[END_REF]. It is important to notice that the negative group velocity is only a necessary condition for apparition of q-waves. Indeed, observation of q-waves also depends upon the real and imaginary part of the amplification factor of the space-time discretization and with an excessive dissipation, filtering or damping, the q-waves can be removed very quickly. Even if it was mentioned in [START_REF] Sengupta | Spurious waves in discrete computation of wave phenomena and flow problems[END_REF] that initial condition, grid resolution and multi-dimensional case have effects on q-waves, the analysis presented in the following is restricted to the one-dimensional case. A mesh composed of 300 cells is considered with the manufactured choice of ω j introduced in Sec. 5.1 and used again in Sec. 6.1. The numerical phase speed c N j and the numerical group velocity V gN j are computed from G j (Eq. 21) for any cell j (and therefore for a given ω j ). The analysis of q-waves is performed for the AION scheme and for CFL ∈ [0, 1]. Attention is paid on the computation of the group velocity V gN j for the different intersection of time integrators in AION scheme. In Figs. [START_REF] Limare | An adaptive remeshing strategy for unsteady aerodynamics applications[END_REF][START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible RANS-LES simulations on fully unstructured grids[END_REF][START_REF] Charrier | Hybrid RANS/LES simulation of a supersonic coaxial He/air jet experiment at various turbulent lewis numbers[END_REF][START_REF] Carpaye | Design and analysis of a task-based parallelization over a runtime system of an explicit finite-volume cfd code with adaptive time stepping[END_REF], dashed curves correspond to |V gN |/|c| >> 1, which represents discontinuity of dispersion φ j , and grey zones correspond to negative group velocity V gN j < 0. Remark: The choice of the CFL domain is dictated by the desire to compare several schemes. Of course, amplification always occurs for the Heun's scheme at CFL> 0.6 but the analysis can still be performed for CFL> 0.6.

Concerning the areas of negative-V gN j group velocity, it can be noticed that the zone is slightly more expanded for Heun/Hyb cells than for the pure Heun scheme, according to Figs. 21 and22. Indeed, negative group velocities can be found up to CFL=0.57 at Heun/Hyb cells, but the full Heun's scheme negative-group-velocity zone is present up to CFL=0.47. In addition, this negative group-velocity zone is larger for the other transitions according to Figs. 23 and24. So, the presence of q-waves is possible in a larger domain of CFL for AION scheme than for the full Heun's scheme. The last question concerns the capability of the AION scheme to damp q-waves. To answer it, it is mandatory to couple the observations on negative group velocity with the dissipation property of the AION scheme, as shown in Figs. 2627282930. According to Fig. 26, the CFL limit of the negative-V gN zone (CFL=0.47) due to Heun time integration corresponds to µ j < 0.2. In the case of the time integrator transitions of AION scheme, the CFL limit of their negative-V gN zone corresponds to µ j < 0.3 for cell Heun/Hyb (Fig. 27), µ j < 0.5 for hyb/IRK2 and IRK2/Hyb cells (Fig. 28 29), and finally µ j < 0.2 for cell Hyb/Heun cells (Fig. 30). Fortunately, it seems that the effect of extension of the q-waves zone due to our coupling scheme can be attenuated by dissipation. All this information allows us to conclude that the AION scheme seems to be as stable as the full Heun scheme on uniform mesh. Moreover, the possible additional q-waves that could propagate would be quite dissipated. The good numerical behaviour obtained through the proposed theoretical analysis must now be confirmed by numerical simulations, which is the goal of the next sections.

Validation

This section is devoted to the validation of the AION scheme using several test cases of increasing complexity, starting from 1D advection solution and Euler solutions to 3D Navier-Stokes computations.

Advection of a Gaussian hump

The first test case is the linear advection of a Gaussian hump (Eq. 1 with F (W ) = W and G(W, ∇W ) = 0 in a periodic domain) and serves to assess both stability and accuracy of the AION, Heun and IRK2 schemes. Introducing the periodic domain [0, 1] and the parameters of the initial hump x 0 = 0.5m, σ = 0.1, the analytic theoretical initial solution is simply:

y = 1 √ 2πσ . exp - (x -x 0 ) 2 2σ 2 . (28)
The initial solution must be recovered at any time t = k s with k ∈ N and here, the analysis is performed at t = 10s. The mesh is composed of N = 333 cells and is irregular but manufactured using the following rule:

∆x j = α ∆x j-1 for N 2 -50 ≤ j ≤ N 2 , ∆x j = 1 α .∆x j-1 for N 2 + 1 ≤ j ≤ N 2 + 50.
∆x j = ∆x max elsewhere [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] with α = 0.9. Then the parameter ω j is controlled as :

ω j = ∆x j ∆x max (30)
Such a variation in mesh size enables a strong control of the maximum CFL value authorized for a computation using an explicit time integration. In addition, the coupling scheme is fully explicit in the regular part of the grid (ie. when ∆x j = ∆x max ). The spatial scheme used is a second order spatial scheme for regular and irregular grids. The time step is fixed in order to obtain a CFL constant in regular grid (= CFL min ) and leads to higher values in the smaller cells of the irregular grid. The first computations are performed using IRK2, Heun and AION scheme at CFL min = 0.005 in the explicit part. The largest CFL value in the implicit part is 1. Fig. 31 illustrates the fact that the AION scheme conserves the time accuracy of Heun's scheme but the IRK2 scheme involves more dissipation. At this level, it must be mentionned that the Newton algorithm for the convergence of IRK2 and AION scheme is guaranteed: the final error (at any time step) is very low. This is to avoid any confidence in the results due to a bad convergence of the Newton algorithm. Now, the same test case is performed with a higher CFL min = 0.1. In that case, the maximum CFL number in the implicit part reaches 20, out of the stability region of Heun's scheme. In this configuration, both AION and IRK2 schemes are stable. Here again, the IRK2 scheme is more dissipative than the AION scheme. Moreover, for the present case, the convergence of the implicit Figure 31: Advection of a Gaussian hump at CF L min = 0.005, effect of the time integration scheme Figure 32: Gaussien advection at CF L min = 0.1 schemes (AION and IRK2) is performed until a low error on the solution. The number of Newton steps is much larger for the IRK2 scheme than for the AION scheme.

The last question of importance concerns the interest of using the AION scheme at moderate CFL numbers versus the IRK2 scheme at larger CFL value. The computation is now performed with the IRK2 scheme at an associated CFL min = 1 in the explicit part. At this CFL number, the AION computation cannot be performed since Heun's scheme is now unstable. As expected (Fig. 33), the IRK2 scheme is usable in this case but the number of iterations to converge Newton's algorithm has a strong influence on the solution accuracy. In our example, the reference solution needs 90 steps for Newton's algorithm and the solution is still dissipated and dispersed. In order to try to represent an industrial computation, the same simulation is now performed but only 60 steps of Newton's algorithm are allowed. The solution is more dispersive. For the sake of consiceness, the AION solution at CFL min = 0.1 is recalled. For the implicit part of the scheme, the convergence of Newton's algorithm does not require more than 10 steps.

Indeed several articles mentionned that the time integration of LES and DNS simulation by implicit time integrators involve a damping of the high-frequency solution modes for large time steps [START_REF] Cinnella | High-order implicit residual smoothing time scheme for direct and large eddy simulations of compressible flows[END_REF][START_REF] Martn | A parallel implicit method for the direct numerical simulation of wall-bounded compressible turbulence[END_REF]. Furthermore the number of Newton subiteration is not easy to control in order to obtain convergence of the error. It remains more interesting to perform implicit time integration in stiff and/or stationnary

Sod's tube

The Sod's tube is an unsteady inviscid one-dimensional problem used to characterize the properties of both time integration and convection schemes. The computational domain of length L x = 1 m is split into two parts separated by a membrane located initially at x = 0.5L x . The initial flow is defined by

     ρ L p L U L      =      1.0 1.0 0.0      ,      ρ R p R U R      =      0.125 0.1 0.0      , (31)
where L refers to the left side and R to the right side of the membrane. At t = 0, the membrane blows up and waves are travelling inside the computational domain. At the final time t = 0.2s, the solution is composed of a rarefaction wave, a contact discontinuity and a shock. The Euler equations are solved using a 1-exact (second order) upwind scheme using the Roe approximate Riemann solver. In order to avoid spurious oscillation and to compute TVD solution, our limiting strategy is provided with the minmod slope limiter [START_REF] Roe | Characteristic-based schemes for the euler equations[END_REF]. The computational domain is composed of N = 300 cells located in regular parts with a uniform mesh size and in irregular parts with non-uniform mesh size such that

for N 2 -50 ≤ j ≤ N 2 : ∆x j = α.∆x j-1 for N 2 + 1 ≤ j ≤ N 2 + 50 : ∆x j = 1 α .∆x j-1 (32)
with α = 0.98. Then the parameter ω j is controlled as:

ω j = ∆x j ∆x max . ( 33
)
This definition of ω j leads to a fully explicit time integration on the regular part of the grid, where ∆x j = ∆x max . In the following, ν j is defined by ν j = ∆t ∆x j . For ν min = 0.1 in the explicit part, the maximum value of ν in the implicit part reaches 0.3. Figures 34 and35 show the density and velocity profiles, with zooms in the region of the rarefaction wave and near the contact discontinuity. The AION scheme appears to be as accurate as the second order time accurate standard schemes and in agreement with the solution obtained using the Euler explicit scheme. A second set of computations is performed for ν min = 0.4 in the explicit part, which leads to a maximum value of ν in the implicit part of 1.01. The Heun' scheme is unstable, but Figs 36 and 37 show that the AION and IRK2 schemes are still stable. Both solutions are very close. Paying attention to the velocity and density near the shock, it seems that the AION scheme slightly dissipates the overshoots given by the IRK2 implicit scheme. From this case, it can be concluded that the reconstruction procedure to define the AION scheme is in agreement with the requirements to handle shock, contact discontinuity and rarefaction wave. Moreover, the AION scheme improves the explicit scheme with an enhanced stability and is more computationally efficient than the fully-implicit scheme. The next test case is dedicated to the analysis of the scheme accuracy.

Two-dimensional linear advection of an isentropic vortex

The advection of an isentropic vortex is a famous test case of the High Order Workshop. This problem is designed in order to test the solver capability to preserve vorticity in an unsteady simulation and this verification case is suitable for verifying total order of accuracy. The total error E total may be expressed as:

E total = Ah p + B∆t q + O(h p+1 , ∆t q+1 ) (34)
with (A, B) ∈ R 2 . The total order of accuracy is estimated as min(p, q). For self-consistence, the test case definition is summarized below. The computational domain is the square domain [-L 2 , L 2 ] 2 (L = 0.1) with periodic boundary conditions. On a uniform flow of pressure P ∞ , temperature T ∞ and Mach number M ∞ is superimposed an isentropic vortex defined by its characteristic radius R and strength β. The vortex is initialized in the center of the computational domain (x c , y c ) = (0, 0). The initial state is defined by:

δu = -U ∞ β (y -y c ) R e -r 2 2 , δv = U ∞ β (x -x c ) R e -r 2 2 , δT = (U ∞ β) 2 2 .e -r 2 , u 0 = U ∞ + δu, v 0 = δv, (35)
with:

r = (x -x c ) 2 + (y -y c ) 2 R , U ∞ = M ∞ γ R gas T ∞ . (36)
with R gas = 287.15 J/kg/K the gas constant and a constant ratio of specific heats γ equal to 1.4.

The isentropic relation leads to the complete set of unknowns

T 0 = T ∞ -δT, ρ 0 = P ∞ R gas T ∞ T 0 T ∞ 1 γ-1 , P 0 = ρ 0 R gas T 0 (37)
Here, the "fast vortex" test case is considered and defined by

P ∞ = 10 5 N/m 2 , T ∞ = 300 K, M ∞ = 0.5, β = 1 5 , R = 0.005. ( 38
)
The solution is time-marched during three rotations inside the periodic box. The computation is performed with a Successive-Correction 2-exact formulation for the spatial scheme (order three) and time integrated by Heun, IRK2 and AION schemes on the baseline Cartesian grids of 64 2 , 128 2 and 256 2 degrees of freedom (DOF). The time integration of theses computations was performed at CFL=0.1. For the AION scheme, it is mandatory to define the hybrid parameter ω j ; it is configured according to following equation:

ω j = 1 -e -r 2 2 with r 2 = (x j -x c ) 2 + (y j -y c ) 2 R 2 . (39)
Fig. 38 represents the distribution of ω j on the 256 2 mesh. According to the pressure and velocity field obtained with Heun, AION and IRK2 time integrator on grid 256 2 , it appears that time integrators have slightly same properties of dissipation and dispersion. The AION scheme tends to less dissipate than other second-order time integrator (according to Figs. 40 and42).

This feature is quite different regardless the ones found on one-dimensional test case with second-order space reconstruction, where the AION was more dissipative than other time integrators. Here the third-order accurate space reconstruction tends to improve the space-time spectral characteristics of the computation, time integrated with AION scheme. According to the computation of L 2 -norm of the error using the two velocity-vector components (see Fig. 40), it appears, as expected, that the AION time integrator conserves second-order space-time accuracy (see Tab. 4). Nevertheless, according to the slope of log(E total ), it seems that the time integrator has a significant impact on the error computation (see Tabs. [START_REF] Sengupta | Analysis of central and upwind compact schemes[END_REF][START_REF] Sengupta | Space-time discretizing optimal DRP schemes for flow and wave propagation problems[END_REF]. Indeed the value of the slope is quite different between log(h) = [-2.4, -2.1] and log(h) = [-2.1, -1.8], particularly in case of the Heun scheme. The computational cost of the AION and the IRK2 time integrator at same CFL is also compared in Tab. 3. The ratio between cells with ω j < 1 and cells with ω j = 1 (I/E cells) is equal to 18% for all grids. It appears that the AION scheme reaches an overall cost that is about nearly 50% lower than the pure implicit IRK2 scheme for grid with 256 2 DOFs at same CFL.

(h) = [-2.1, -1.8]
The convection of a vortex on an irregular grid of 260 2 DOFs is also performed. The time integration is performed according to CFL condition of the biggest cells of the domain such as

∆t = CF L max j (h) v j , (40)
In this irregular domain, the size ratio between the largest and the smallest cell is equal to 11. Hence, the constant time step of computation will be determined according to the CFL number of the largest cell, which leads to the CFL number 11 time higher for the smallest cells than for the larger ones.

The computation is performed at several CFL values, during one time period, in order to compare stability properties of the several time integrators. The parameter ω j is configured such as

ω j = |Ω j | max j (|Ω j |) , (41)
The parameter ω j is defined according to the size of cells in the domain (see Fig. 44 with grey zone corresponding to ω j = 1). In this case the proportion of cells with ω j < 1 corresponding to hybrid and implicit cells is equal to 73%. The final computational times for the several time integrator and CFL numbers is provided in Tab. 7. The values are normalized with respect to computational cost of the Heun time integrator at CFL=0.1. According to Tab. 7, as expected, the explicit time integrator is unstable for high CFL number due to small cells. Then, the AION scheme, as expected, tends to have better stability properties than full Heun time integrator, due to implicit time treatment in the smallest cells. The computational gain compared to full implicit computation is equal to 7%, due to high percentage of hybrid and implicit cells. In the following a focus is made on testing the hybrid treatment of the gradient for the diffusion scheme.

Three-dimensional Taylor Green Vortex

A three-dimensional laminar flow computation is now performed to assess the performance of the AION integration scheme for a viscous test case. The proposed test case is the Taylor Green vortex at Reynolds number (Re = ρ∞U∞L µ ) equal to 1600 and Mach number equal to 0.1. The initial solution is defined inside the periodic cubic domain [-πL, πL] 3 (L = 1):

u = U ∞ sin x L cos y L sin z L , v = -U ∞ cos x L sin y L cos z L , w = 0, p = P ∞ + ρ ∞ U 2 ∞ 16 cos 2x L + cos 2y L cos 2z L + 2 (42)
This test case is defined in order to observe a kind of vortex cascade, as encountered in turbulent flow computations. For the computation, the fluid is assumed to be a compressible perfect gas with γ = 1.4 and the Prandtl number (P r = µ.cp κ ) is equal to 0.71, with c p and κ, respectively, the heat capacities at constant pressure and the heat conductivity. The initial temperature is kept constant T ∞ = P∞ Rgasρ∞ with R gas the perfect gas constant. From the local pressure and temperature is deduced the density ρ = p RgasT∞ . The simulation is performed until the physical time of t = 20t c with t c = L U∞ the characteristic convective time is attained. The computation is performed with a 2-exact formulation for the spatial discretisation (three order accurate spatially) and time-integrated by Heun and AION schemes on a 256 3 grid at CF L = 0.9. The following outputs are performed:

• The temporal evolution of the kinetic energy integrated over the whole domain:

E k = 1 ρ ∞ Ω Ω ρ U • U 2 dΩ. (43)
• The temporal evolution of the kinetic energy dissipation rate:

ε = - dE k dt . (44)
• The temporal evolution of the kinetic energy dissipation rate on the whole domain:

E = 1 ρ ∞ Ω Ω ρ w • w 2 dΩ, (45)
with w the vorticity.

These outputs are compared with the DNS results obtain on a 512 3 grid with pseudo-spectral method, reference for the First International Workshop on High-Order CFD Methods held at the 50th AIAA Aerospace Meeting (https://www.grc.nasa.gov/hiocfd/). It appears that the AION scheme seems to be quite performing for this kind of viscous case. Indeed the AION scheme fits better the DNS results than Heun scheme for time range t = [0t c , 8t c ] and t = [12t c , 20t c ] according to the E k and ε (see Fig. 45 -46) when Heun scheme seems to be more dissipative in these time ranges. Nevertheless the Heun scheme seems to better reach the maximum of ε than the AION scheme on the time range t = [8t c , 12t c ]. On the contrary, the AION scheme gives better result for the enstrophy E, on the whole time range than Heun scheme (see Fig. 47). This results confirm the particularity that the Heun time integrator has a significant importance on accuracy of the numerical results. Indeed this conclusion was already observed on the slope of log(E total ) obtained in the previous case of the convected vortex. This comparison allows us to validate our hybridation of the viscous gradients presented previously.

Figure 45: The temporal evolution of the kinetic energy coupling procedure. In addition, in order to minimize the CPU cost, transition to explicit and implicit schemes was performed by reducing the transition area, playing with the values of ω. After the spectral analysis, attention was paid on simulations of increasing complexity, from the 1D Gaussian hump advection to 3D Taylor Green Vortex. Three aspects are of paramount importnace. First, the results with the new scheme were shown to have same or better quality than the standard basic schemes. Moreover, starting from the reference full-implicit time integration, the hybrid formulation enables to reduce the CPU cost. Finally, the AION formulation seems less dependent on the Newton's algorithm convergence than the standard IRK2 scheme. Moreover, the convergence of the IRK2 in a fully unsteady simulation needs many steps, something which was not necessary with the proposed AION scheme. We are working today on two new improvements. The first one deals with the use of an adaptive time-step procedure. In this case, explicit cells are separated into several classes depending of their local maximum stable time step. All the cells inside the same class are time-integrated using the same time step and for two adjacent classes, time step for the smallest cells is half the one for the largest cells. Attention must be paid on the computation of the flux for faces having the cell of one side in a given class and the cell on the other side in another one. This work is ongoing and will be summarized in a new paper. The last improvement concerns the application of the proposed technique to a coupled RANS/LES computation for demonstrating the capability in an industrial environment. We can assume that the corner stone will be the scaling of the hybrid parameter ω with the function that switches between RANS and LES models.

Figure 1 :

 1 Figure 1: 1D example to explain cell status and the flux conservation property of the coupling scheme HCS1

Figure 2 :

 2 Figure 2: Isocontours of the dissipation coefficient µ j for the HCS1 scheme at CFL=0.1

Figure 3 : 1 Figure 4 :

 314 Figure 3: Isocontours of the dispersion coefficient φ j for the HCS1 scheme at CFL=0.1

Figure 5 :

 5 Figure 5: Dissipation µ j of the scheme HCS1 for all kinds of cells at CFL=0.1, zoom from Fig. 4

Figure 6 :

 6 Figure 6: Dispersion φ j of the scheme HCS1 for all kinds of cells at CFL=0.1

Figure 9 :Figure 10 :

 910 Figure 9: Isocontours of the dissipation µ j for the AION scheme and 1 ≤ j ≤ N at CFL=0.1

Figure 11 : 1 Figure 12 :

 11112 Figure 11: Dissipation µ j of the scheme AION depending on the kind of cell at CFL=0.1

Figure 13 :

 13 Figure 13: Dispersion φ j of the scheme AION depending on the kind of cell at CFL=0.1

Figure 14 : 6 Figure 15 :

 14615 Figure 14: Dissipation µ j of the scheme AION depending on the kind of cell at CFL=0.6

Figure 16 :

 16 Figure 16: Dispersion φ j of the scheme AION depending on the kind of cell at CFL=0.6

Figure 18 :Figure 19 :

 1819 Figure 18: Isocontours of the dissipation µ j for Hyb/IRK2 transition at CFL=0.5 as a function of ω j . The coupled approach is always stable

Figure 20 :

 20 Figure 20: Isocontours of the dissipation µ j for IRK2/Hyb transition at CFL=0.5 as a function of ω j .The coupled approach is always stable.

Figure 21 :Figure 22 :

 2122 Figure 21: Isocontours of V gN for full Heun scheme

Figure 23 :

 23 Figure 23: Isocontours of V gN for intersection Hyb/IRK2 scheme

Figure 25 :Figure 26 :

 2526 Figure 25: Isocontours of V gN for intersection Hyb/Heun scheme

Figure 27 :

 27 Figure 27: Isocontours of the dissipation µ j for a cell with Heun/Hyb reconstruction and grey zone for negative group velocity

Figure 29 :

 29 Figure 29: Isocontours of the dissipation µ j for a cell with IRK2/Hyb reconstruction and grey zone for negative group velocity

Figure 33 :

 33 Figure 33: Comparaison with IRK2 scheme at CF L min = 1

Figure 34 :

 34 Figure 34: Sod's shock tube at ν min = 0.1. Global view of the density profile at t = 2s and close-up views near the rarefaction wave and near the shock.

Figure 35 :

 35 Figure 35: Sod's shock tube at ν min = 0.1. Global view of the velocity profile at t = 2s and close-up views near the rarefaction wave and near the shock.

Figure 36 :

 36 Figure 36: Sod's shock tube at ν min = 0.4. Global view of the density profile at t = 2s and close-up views near the rarefaction wave and near the shock.

Figure 37 :

 37 Figure 37: Sod's shock tube at ν min = 0.4. Global view of the velocity profile at t = 2s and close-up views near the rarefaction wave and near the shock.

Figure 38 :

 38 Figure 38: Behaviour of ω j for the manufactured solution ; the grey zone corresponds to ω j = 1

Figure 39 :Figure 40 :

 3940 Figure 39: Velocity field (CFL=0.1)

Figure 41 :Figure 42 :

 4142 Figure 41: Pressure field (CFL=0.1)

Figure 44 :

 44 Figure 44: Value of ω j

Table 2 :

 2 Cell position, legend, and associated flux types. Note that left and right sides concern a 1D domain with a positive advection speed.

	Legend	Flux type on left interface Flux type on right interface Cell position
	Heun	explicit	explicit	j = 31
	Heun/Hyb	explicit	hybrid	j = 100
	Hyb/IRK2	hybrid	implicit	j = 104
	IRK2	implicit	implicit	j = 150
	IRK2/Hyb	implicit	hybrid	j = 196
	Hyb/Heun	hybrid	explicit	j = 200

Table 3 :

 3 Elapsed CPU time to reach 3 time periods with AION and IRK2 time integrators .

	DOFs	64 2 128 2 256 2
	I/E Cells	18% 18% 18%
	AION/IRK2 0.75 0.66 0.52

Table 4 :

 4 Total error slopes

	Time integrator log(E total ) slope
	Heun	1.87
	AION	2.22
	IRK2	2.22

Table 5 :

 5

	Heun	2.54
	AION	2.59
	IRK2	2.54

Total error slopes for log(h) = [-2.4, -2.1]

Figure 43: Order of accuracy Time integrator log(E total ) slope

Table 6 :

 6 Total error slopes for log

Table 7 :

 7 Elapsed CPU time to reach 1 time period with AION and IRK2 time integrators on a 2602

	irregular mesh .
		CFL=0.1 CFL=0.4 CFL=0.9
	Heun	1	×	×
	AION	3.405	0.897	0.412
	IRK2	3.64	0.948	0.434

FR4009261,

27/07/2013

RANS and LES equations have almost the same shape and their own pros and cons. In an industrial environment, it could be useful to perform LES far from the wall and RANS near the wall in order to account easily for the turbulence effects, while keeping an acceptable mesh size. Standard time integration schemes for LES and RANS equations do not follow the same constraints. For LES, an explicit time integration is chosen in order to control simply and efficiently spectral properties (dissipation and dispersion). For RANS equations, it is of paramount importance to reach the steady state as fast as possible, using an implicit time integration procedure.

Indeed, the current paper specifically addresses how explicit and implicit time integrators can be coupled spatially. Here, two standard time integration schemes (Heun and second order implicit Runge-Kutta schemes -IRK2-) are hybridized / blended using a transition function ω, while keeping the standard expected properties (spectral behaviour). In the first part of the paper, a way to couple the proposed schemes, adapted from the literature, was first introduced but lead to instability for some wavenumbers and CFL values. A new alternative approach, named AION scheme, was proposed and designed in order to correct the unexpected behaviour of the first coupling procedure. The spectral analysis performed on the coupled space / AION schemes enabled us to check the stability of the