D De Santis
email: dante.de_santis@inria.fr

High-order linear and non-linear residual distribution schemes for turbulent compressible flows

Keywords: High order schemes, residual distribution, gradient reconstruction, RANS equations, implicit method

A high-order Residual Distribution scheme for the solution of the compressible RANS equations is developed. The one-equation Spalart-Allmaras turbulence model is solved in a fully coupled approach; the mean flow equations as well as the turbulence equation are solved with high-order of accuracy. A continuous approximation of the solution is adopted and standard Lagrangian basis functions are used to construct the discrete space. Since viscous terms involve the gradient of the numerical solution which has a discontinuous normal component across the faces of the elements, a continuous approximation of the gradient of the numerical solution is recovered at each degree of freedom of the grid and then interpolated with the same basis functions used for the solution. The non-linear system of equations resulting from the numerical discretization is solved with the non-linear LU-SGS method. Several numerical experiments are performed to verify the accuracy of the numerical method, these are based on both subsonic and transonic flows in two and three spatial dimensions.

Introduction

The majority of the CFD tools used for aerodynamic simulations in the aeronautical industry is based on theoretically second-order Finite Volume (FV) methods [START_REF] Ekaterinaris | High-order accurate, low numerical diffusion methods for aerodynamics[END_REF], whose order of accuracy in practice ranges between one and two due to irregular and highly stretched meshes used. In order to obtain mesh independent solutions very fine meshes with a very large number of degrees of freedom are required, which lead to enormous computational times. High-order methods have the potential to achieve a faster reduction of the discretization errors, maintaining limited the increment of the computational effort.

Probably, one of the most popular scheme in the high-order community is the Discontinuous Galerkin (DG) method [START_REF] Cockburn | Discontinuous Galerkin methods: theory, computation and application[END_REF], which is very attractive for its compactness and flexibility. On the other hand DG methods suffer from the serious drawback of a very fast growth of the number of degrees of freedom (DOFs) with the cell polynomial degree, making these schemes very expensive. Furthermore, the treatment of discontinuities in the DG framework is still not clear and the non-oscillatory properties of DG schemes are still not well understood. An interesting alternative to DG schemes is represented by Residual Distribution (RD) methods. In the RD framework the solution is continuous and the number of DOFs is much lower than DG methods for the same order of approximation of the solution. Furthermore, in the RD framework is very clear how to construct non-linear schemes, which overcome the Godunov's barrier and can be at the same time high-order accurate in smooth regions and monotone across the discontinuities.

The basic idea of the RD approach is to define for each element an integral quantity (the total residual) which is the integral over the element of the governing equations. Such quantity is then distributed to the DOFs of the elements. The way the distribution process is performed determines the behavior of the numerical scheme. For example, in the past, truly multidimensional schemes have been developed. In these schemes, for each element, only upwind nodes can receive a contribution in the distribution process, while downwind nodes receive no contribution [START_REF] Deconinck | Residual distribution schemes: Foundations and analysis[END_REF]. Due to this multidimensional character, these schemes are much more accurate that standard second-order FV methods.

The drawback of the multidimensional upwind schemes is that the extension to orders of accuracy more than two and to non simplex elements is not straightforward, due to the fact that the construction of the numerical schemes relies on geometrical considerations. More flexible approaches make use of centered RD schemes in which, for each element, all the DOFs receive a contribution during the distribution process, regardless if they are in the upwind or downwind region. For example, finite element approaches like SUPG [START_REF] Hughes | Recent progress in the development and understanding of supg methods with special reference to the compressible euler and navier-stokes equations[END_REF] or Lax-Wendroff [START_REF] Selmin | Finite element methods for nonlinear advection[END_REF] can be interpreted as centered RD schemes. Recently, Abgrall [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] proposed a non-linear RD method in which a simple low-order centered scheme is used as base for the construction of a higher-order scheme that is still monotone near shocks. An additional advantage of this approach is that the numerical method can be easily extended to any order of approximation and to any type of element.

Although the RD discretization of advective problems is very clear, the extension of RD methods to advectiondiffusion problems is not straightforward. The simple blending of a RD scheme for advection terms with a Galerkin discretization of the diffusive terms usually produces a sub-optimal scheme in regions where advection and diffusion have the same importance. A different strategy consists in using a single RD scheme for the whole advection-diffusion equation. However, a new difficulty arises due to the fact that, for polynomial piecewise approximation of the solution, the normal component of the gradient of the numerical solution is discontinuous on the face of two adjacent elements. Instead of using a numerical flux for the viscous terms, leading to the introduction of additional parameters that might be difficult to understand, one could think of reconstructing an unique set of values of the gradient of the numerical solution at the DOFs. This set of values is used to interpolate the gradient with the same continuous functions used to interpolate the solution, making straightforward the calculation of the total residual of the whole advection-diffusion equations.

In [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] this approach was studied in detail considering the scalar advection-diffusion problem. In particular different gradient reconstruction strategies were considered and the effect of the gradient reconstruction on the accuracy of the numerical solution was studied. It was found that, in order to preserve the overall accuracy of numerical scheme, the gradient of the numerical solution must be reconstructed with a procedure such that the gradient has the same order of accuracy of the solution. Several gradient reconstruction techniques were considered and it was found that the so-called Super-convergent Patch Recovery method of Zienkiewicz and Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF] (SPR-ZZ) was able to reconstruct the gradient with the optimal accuracy without enlarging the stencil of the numerical scheme.

In [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF] the scheme has been extended to the Navier-Stokes equations. Numerical experiments confirmed that the optimal accuracy of the numerical scheme is lost if the gradient of the numerical solution is not properly reconstructed. The SPR-ZZ method has been used in the contest of the Navier-Stokes equations to construct linear and non-linear RD schemes for the discretization of steady laminar flows at the second and third order of accuracy.

In this work the RD numerical solver is extended to RANS equations, showing how high-order RD schemes might be very effective to get highly accurate solutions with a reduced number of DOfs than a lower-order discretization. In particular, non-linear RD schemes might be very attractive for simulations of transonic turbulent flows, because no artificial viscosity is required to get monotone solutions near discontinuities.

The structure of the paper is as follows. In Section 2, the multidimensional RANS equations with the Spalart-Allmaras turbulence model fully coupled with the mean flow equations are described. In Section 3, the proposed RD method for advection-diffusion problems is described. The detail of the construction of linear and non-linear schemes are given. In Section 4, the strategy used to reconstruct the gradient is described in detail. In Section 5, is described the construction a non-linear LU-SGS method for a robust and effective implicit scheme. In Section 6, numerical experiments are reported; both subsonic and transonic flows in two and three spatial dimensions are considered. Finally, a few conclusions are drawn in Section 7.

Governing equations

The compressible Reynolds Averaged Navier-Stokes (RANS) equations coupled with the one equation Spalart-Allmaras turbulence model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] can be written in the following conservative form

∂ρ ∂t + ∇• m = 0 ∂m ∂t + ∇• m ⊗ m ρ + PI = ∇• S ∂E t ∂t + ∇• E t + P m ρ = ∇• S • m ρ -q ∂µ ⋆ t ∂t + ∇• µ ⋆ t m ρ = ∇• η σ S A ∇ µ ⋆ t ρ + S S A , (1)
where the unknowns are the density ρ, the momentum vector m, the total energy per unit volume E t and the turbulent working variable µ ⋆ t , with m = ρv and

E t = ρe + m 2 2ρ ,
where v is the velocity vector and e is the specific (i.e., for unit of mass) internal energy. Furthermore, in the system (1) P is the static pressure, I ∈ R N dim is the identity matrix, with N dim number of spatial dimensions, S is the viscous stress tensor and q is the heat flux vector. The viscous stress tensor can be written as

S = - 2 3 (µ + µ t )(∇• v)I + (µ + µ t )(∇ T v + ∇v),
where µ is the fluid dynamic viscosity and µ t is the turbulent eddy viscosity, which is given by

µ t =        µ ⋆ t f v1 , µ ⋆ t > 0, 0, µ ⋆ t ≤ 0, (2)
where

f v1 = χ 3 χ 3 + c 3 v1 , with χ = µ ⋆ t µ t , c v1 = 7.1.
The heat flux, using the Fourier law, can be written as: q = -(κ + κ t )∇T , where T is the temperature and the thermal conductivity is evaluated using the constant Prandtl number hypothesis

κ = c p µ Pr , and
κ t = c p µ t Pr t ,
with c p the specific heat at constant pressure, Pr the laminar Prandtl number and Pr t the turbulent Prandtl number.

Here, Pr = 0.72 and Pr t = 0.9.

In the Spalart-Allmaras equation η = µ + µ t and σ S A = 2/3, while the source term is defined as follows

S S A = c b1 ωµ ⋆ t -ρc w1 f w µ ⋆ t ρ d min 2 + 1 σ S A ρc b2 ∇ µ ⋆ t ρ 2 , (3)
the three terms in the source term represent respectively the production, the destruction and the diffusion. In the production term, the modified vorticity function is defined as

ω = ω + ω
where ω is the magnitude of the vorticity vector, and

ω = µ ⋆ t f v2 ρk 2 d 2 min , with f v2 = 1 - χ 1 + χ f v1 .
The remaining closure functions are

f w = g       1 + c 6 w3 g 6 + c 6 w3       1 6 , g = r + c w2 r 6 -r , r = µ ⋆ t ρ ωk 2 d 2 min ,
with d min the minimum distance from the wall, c b1 = 0.1355, c b2 = 0.622, k = 0.41,

c w1 = c b1 k 2 + (1 + c b2 ) σ S A , c w2 = 0.3
and c w3 = 2.

The system (1) is closed by providing a thermodynamic model for the fluid. In this work only the polytropic ideal gas model is considered, for which the pressure and the temperature can be written as P(e, ρ) = (γ -1)ρe, and

T (e) = (γ -1) R e,
with γ = c p /c v the specific heats ratio and R the gas constant.

If the vector of the conservative variables u ∈ R N eq is introduced, with N eq the number of equations, the system of equations ( 1) can be written in a more compact form as follows

∂u ∂t + ∇• f a (u) -∇• f v (u, ∇u) = S(u, ∇u), (4)
with

u =                              ρ m E t µ ⋆ t                              , f a =                               m m ⊗ m ρ + PI E t + P m ρ µ ⋆ t m ρ                               , f v =                              0 S S • m ρ -q η σ S A ∇ µ ⋆ t ρ                              and S =                              0 0 0 S S A                              .
where f a and f v are the advective and the viscous flux functions, respectively, and S is the source term. For later convenience, it is pointed out here that, since the viscous flux function is homogeneous of degree one with respect to the gradient of the conservative variables, it can be also written as follows

f v (u, ∇u) = K(u)∇u,
where the viscous flux Jacobian K(u) ∈ R N eq ×N dim ×N eq ×N dim is defined as

K(u) = ∂f v (u, ∇u) ∂∇u , with K i j (u) = ∂f v i (u, ∇u) ∂ ∂u ∂x j .
The steady state case, is obtained by removing the time derivative in Eq. ( 4). The equations are made dimensionless by taking as reference quantities the density, the speed of sound and the temperature at the free-stream and a characterizing length of the problem.

Modification of the Spalart-Allmaras model

Although the use of the Spalart-Allmaras turbulence model is quite common in numerical simulations of aerodynamic flows, getting a robust solver is still a challenging task. A very important issue of the Spalart-Allmaras model is the behavior of the turbulent governing equation in zones with a high gradient of the working variable. High values of the gradient are generally found near solid walls, where the turbulent variable rises rapidly from zero to a high value, and in the outer zone of the boundary layer where the turbulent variable rapidly decrease to the free-stream value.

While the space discretization is generally such fine in the near-wall zones that the turbulent variable is computed monotonously, the mesh resolution in the outer part of the boundary layer may be insufficient to resolve the turbulent field. In this case, the turbulent field shows numerical oscillations with negative values of the turbulent variable, which can make the algorithm unstable.

Negative values of the turbulent working variable are disregarded in the definition of the eddy viscosity (see Eq. ( 2)), but they can compromise the stability of the numerical solver. Several modifications have been proposed in literature to remedy to the problem caused by negative values of µ ⋆ t . Some of them are based on simple clippings in the internal functions of the Spalart-Allmaras equations. These approaches, however, might have a negative effect on the implicit iterative solver. In addition, while these approaches might be effective for FV methods, they are usually not robust enough for high-order schemes.

When µ ⋆ t is negative, the turbulence model equation should be modified in order to remove the destabilizing effects. Some used modifications are based on the idea suggested by Allmaras in [START_REF] Allmaras | Multigrid for the 2-D Compressible Navier-Stokes Equations[END_REF]. In this approach, one writes the equation for the energy of the negative values of the turbulent working variable and the turbulence model is modified in such a way that the rate of change of the unphysical turbulent energy is negative, see [START_REF] Oliver | High-Order, Adaptive, Discontinuous Galerkin Finite Element Method for the Reynolds-Averaged Navier-Stokes Equations[END_REF] for more details. This approach makes the Spalart-Allmaras equation stable in the zone with negative µ ⋆ t and makes the equation continuous and differentiable at µ ⋆ t = 0. However, the modifications remove the negative contribution from the source term by making the production very strong in the negative µ ⋆ t region. Although the modification is mathematically stable, it has been observed in numerical simulations that the strong effect induced by the modifications of the production term can induce strong variations of the residual during the iterative convergence, with a possible solver failure. In order to have a more robust solver, the elimination of negative values of µ ⋆ t must be addressed directly. In the actual implementation of the Spalart-Allmaras equation, two modifications have been adopted which have shown to be very effective. The first modification consists in replacing the variable χ with the following function [START_REF] Moro | Navier-Stokes solution using hybridizable discontinuous Galerkin methods[END_REF]

ψ µ ⋆ t =        0.05 log 1 + e 20χ , χ ≤ 10, χ, χ > 10,
which has the role to deactivate the production and the destruction terms of the turbulence model equation when µ ⋆ t becomes negative. The function ψ µ ⋆ t is continuous and differentiable for all values of χ as can be observed in Fig. 1. The form of the turbulence model would be identical to the original one with the exception that the variable χ is replaced by ψ µ ⋆ t . The second modification consists in getting rid of the diffusion contribution in the source term when µ ⋆ t ≤ 0.

In conclusion, in the improved version of the Spalart-Allmaras turbulence model equation, the eddy viscosity is given by

µ t = µψ µ ⋆ t f v1 , (5)
and the turbulent model equations is written as

∂µ ⋆ t ∂t + ∇• µ ⋆ t m ρ - η σ S A ∇ µ ⋆ t ρ = S S A , with η = µ(1 + ψ µ ⋆ t )
and the source term given by

S S A =                  c b1 ωµψ µ ⋆ t -ρc w1 f w µψ µ ⋆ t ρ d min 2 + 1 σ S A ρc b2 ∇ µ ⋆ t ρ 2 , χ > 0, c b1 ωµψ µ ⋆ t -ρc w1 f w µψ µ ⋆ t ρ d min 2 χ ≤ 0, (6)
the following modified form of ω is used [START_REF] Allmaras | Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model[END_REF] ω that is identical to the original expression for ω > 0.3 ω , but remains positive for all nonzero ω and is C 1 continuous. The other functions involved in the turbulence model are re-defined as follows

=              ω + ω, ω > -c v2 ω , ω + ω c 2 v2 ω + c v3 ω (c v3 -2c v2 ) ω - ω , ω ≤ -c v2 ω ,
ω = µψ µ ⋆ t f v2 ρ k 2 d 2 min , f v1 = ψ 3 µ ⋆ t ψ 3 µ ⋆ t + c 3 v1 , f v2 = 1 - ψ µ ⋆ t 1 + ψ µ ⋆ t f v1 , f w = g       1 + c 6 w3 g 6 + c 6 w3       1 6 , g = r + c w2 r 6 -r , r = µψ µ ⋆ t ρ ω k 2 d 2 min .
Note that with the definition (6), the source term tends to zero when µ ⋆ t becomes negative. In [START_REF] Crivellini | A Spalart-Allmaras turbulence model implementation in a discontinuous Galerkin solver for incompressible flows[END_REF], similar effects on the source term are obtained by directly setting to zero the source term for negative values of µ ⋆ t .

RD space discretization

In this section the RD space discretization of the steady Navier-Stokes equations is described. The numerical method used to discretize the system of equations follows from the scheme already proposed in [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] to discretize the scalar advection-diffusion problem and has subsequently been extended to discretization of the Navier-Stokes equations in [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF].

The computational domain Ω is discretized with N e non-overlapping elements with characteristic length h, the set of all the elements is denoted by E h , the list of the DOFs is denoted by N h , and the total number of DOFs is N dof . The solution is approximated on each element by k-th order polynomials which are assumed to be continuous within the elements and on the faces of the elements. In this work Lagrangian basis functions are used to approximate the solution as follows

u h (x) = i∈N h ψ i (x) u i , x ∈ Ω,
with u i the numerical solution at the node, i.e., degree of freedom (DOF), i and ψ i the Lagrangian basis functions. Due to the approximation of the solution, the integral of the Eq. ( 4) calculated on each element e will be not null, but it will give rise to a residual on each element, namely

Φ e u h , ∇u h = Ω e ∇• f a u h -f v (u h , ∇u h ) dΩ + Ω e S(u h , ∇u h ) dΩ, (7)
where

∇u h | e = i∈N e h ∇ψ i u i ,
is the internal gradient of the element e, computed with the gradient of the basis functions as in the Finite Element approach, and N e h is the set of DOFs of the element e. The integral quantity Φ e u h , ∇u h is called total residual of the element e.

By using the divergence theorem, the first integral in Eq. ( 7) can be expressed as a surface integral and the divergence of the flux function is no longer required. This procedure is legitimate for the divergence of the advective flux, since the numerical solution by definition is continuous across the face of two neighbor elements. However, in handling the divergence of the viscous flux one has to cope with the fact that the normal component of internal gradient of the numerical solution is discontinuous across the face of two adjacent elements. The continuity of the normal flux across the element boundary is also necessary to show the consistency of the numerical scheme with respect to the original governing equations [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF].

Suppose that ∇u i i=1,..., N dof is the set of unique values of the gradient of the numerical solution, reconstructed at each DOF. The gradient can be now interpolated with the same basis functions used for the solution, namely

∇u h (x) = i∈N h ψ i (x) ∇u i , x ∈ Ω.
The quantity ∇u h is now continuous on the faces of the elements and the total residual can be written as follows

Φ e = ∂Ω e f a u h -f v (u h , ∇u h ) • n d∂Ω + Ω e S(u h , ∇u h ) dΩ, = ∂Ω e f a u h -K u h ∇u h • n d∂Ω + Ω e S(u h , ∇u h ) dΩ. ( 8
)
Note that the reconstructed gradient has been used also in the definition of the source term. This term could be computed also with the internal gradient, but the reconstructed gradient is used for sake of accuracy and robustness. The procedure to recovery the gradient at each DOF, i.e., how to compute ∇u i i=1,..., N dof , is detailed in Section 4.

Note that the idea to use a recovery procedure to remove the discontinuities at the element interfaces has been also used in the contest of DG methods for solving viscous flow problems [START_REF] Luo | A reconstructed discontinuous Galerkin method for the compressible NavierStokes equations on arbitrary grids[END_REF][START_REF] Van Leer | Discontinuous galerkin for diffusion[END_REF]. In these works, starting from the piece-wise discontinuous solution, a continuous approximation of the solution at the cell interface is reconstructed on the union of adjacent cells to the interface, and such continuous solution is used to evaluate the viscous flux function at the element interface.

In order to handle only nodal values, the total residual ( 8) is first distributed, in some way, to the DOFs of the element as follows

Φ e i = β e i u h Φ e u h , ∇u h , ∀i ∈ N e h , (9)
where β e i are the distribution coefficients, which can be in general function of u h . It is easy to see that the nodal residuals, Φ e i , must satisfy the following conservation constrain [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF]

i∈N e h Φ e i = Φ e , ∀e ∈ E h . (10)
Following [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF], formal accuracy is achieved when, for a smooth exact solution u exact which is interpolated by u h , the residuals Φ e i evaluated for that particular function satisfy

Φ e i (u h ) = O(h k+N dim ),
where N dim is the geometrical dimension (N dim = 2, 3) and k is the polynomial degree used to interpolate the solution.

The form (9) satisfies this constraint, provided that the coefficients β e i are uniformly bounded with respect to the numerical solution and the mesh. If a set of residuals satisfies (9), then the scheme is said to be linearity preserving.

To obtain an equation for each nodal value of the numerical solution, the following relations are written for each DOF

e∈E i h Φ e i = 0, ∀i ∈ N h , (11)
where E i h is the set of the elements which share the DOF i. The previous relations define a set of non-linear equations that must be solved for the nodal values of the solution [u i ] i=1,..., N dof . In practice, the solution with a RD method is obtained by the means of an iterative method, which in the simplest form reads

u n+1 i -u n i ∆t n i = - e∈E i h Φ e i , ∀i ∈ N h ,
with ∆t n i a scaled pseudo-time step. The change of the nodal values of the solution during the iterative process is driven by the non-zero total residuals on the elements; for n → ∞ the total residual on each element vanishes and the steady state solution is obtained.

Construction of distribution schemes for advection-diffusion problems

The distribution process in Eq. ( 9) has been expressed through the use of generic distribution coefficients, the focus now is on how to actually construct a distribution scheme for the RANS equations. In this work, the attention is focused on the construction of centered schemes. Linear and non-linear schemes are considered.

Linear scheme

Consider a numerical scheme composed by an isotropic distribution plus a streamline term. This scheme has been successfully used in the context of the RD schemes to discretize the scalar advection-diffusion equation [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] and the Navier-Stokes equations [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF]. The nodal residual, for a generic DOF i of the element e, can be written as follows

Φ e i = Φ e N e dof + Ω e A•∇ψ i Ξ A • ∇u h -∇• K ∇u h + S dΩ, (12)
where N e dof is the number of DOFs of the element e, A = ∂f a /∂ u is the inviscid flux Jacobian, and Ξ is a scaling matrix which is taken as

Ξ = 1 N dim |Ω e |          N e dof i=1 R n i (ū) Λ + n i (ū) L n i (ū) + N dim j=1 K j j (ū)          -1 , ( 13
)
where ū is the arithmetic average of the vector of the conservative variables on the element, R n , L n are respectively the matrices of the right and left eigenvectors along the direction of the vector n, and Λ n = diag(λ n ) is the corresponding diagonal matrix of the eigenvalues. The operator (•) + selects only the positive values and sets the negative ones to zero. The vector n i is taken as follows

n i = 1 N dim Ω e ∇ψ i dΩ,
such that it has the dimensions of a length (surface) in two (three) spatial dimensions. The scheme is consistent since it achieves the conservation constraint [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF]. Formal accuracy is satisfied too: when the exact solution is injected in Eq. ( 13), the residual is zero because the total residual vanishes, and the integral term vanishes as well due to the fact that the terms in the integral at the right of the scaling matrix represent the governing equations.

The scheme is linearity preserving but not positive, which means that it is not guaranteed that it is monotone for discontinuous solutions. Due to the integral formulation, the scheme is valid for any type of element and for any order of approximation of the solution.

Non-linear scheme

Non-linear schemes are needed to combine the non-oscillatory behavior with the higher-order discretization. The basic idea to construct a non-linear scheme is to start with a first-order, monotone scheme, and to map its nodal residuals onto a set of positive and non-linear residuals [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF]. To see in practice how to construct a non-linear scheme, consider the first order accurate and positive Rusanov's scheme (also know as the local Lax-Friedrichs scheme) defined as

Φe i = Φ e N e dof + 1 N e dof α e j∈N e h j i u i -u j ), ∀i ∈ N e h ,
with α e large enough for stability, in practice it is taken as

α e = max j∈N e h       |λ n j | + µ + µ t ρ + κ + κ t c v j       .
Since the Rusanov's scheme is first-order accurate, its distribution coefficients, β e i , are unbounded. The construction of the non-linear scheme consists in mapping the distribution coefficients of the low order scheme onto non-linear bounded distribution coefficients As described in [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF], the mapping for a system of equations is constructed in the characteristic space, consequently the nodal and the total residuals are first rewritten as

Φ e i ⋆ = L n Φe i and Φ e⋆ = i∈N e h Φ e i ⋆ ,
where the mean fluid velocity vector on the element is used as direction vector to compute the eigenvectors. The nodal high-order residuals are obtained as follows Φ e⋆ i = β e⋆ i Φ e⋆ , with the bounded distribution coefficients, β e⋆ i , obtained by applying a non-linear mapping to the original unbounded distribution coefficients, β e⋆ i = Φ e i ⋆ /Φ e⋆ . The map is constructed as follows

β e⋆ i = Φ e i ⋆ Φ e ⋆ + j∈N e h Φ e i ⋆ Φ e ⋆ + ,
which correspond to the so-called PSI limiter. Finally, the high-order nodal residuals are projected back on the physical space: Φe i = R n Φ e⋆ i . The use of a centered scheme, like the Rusanov's scheme, in combination with the limiting technique, produces undamped spurious modes and a poor iterative convergence to the steady state solution [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] due to the fact that no upwind mechanism is included. The cure to this problem consists in adding a filtering term by means of a streamline dissipation term

Φ e i = Φe i + ε e h (u h ) Ω e A•∇ψ i -∇• K∇ψ i Ξ A • ∇u h -∇• K ∇u h + S dΩ. ( 14
)
where the scaling matrix Ξ is defined as in Eq. ( 13). The term ε e h (u h ) is a smoothness sensor which assures that the filtering term is added only in the smooth regions of the solution, namely ε e h (u h ) ∼ 1 in smooth regions and ε e h (u h ) ∼ 0 near discontinuities. A pressure-based shock sensor is used for ε e h (u h ), see [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF],

ε e h (u h ) = 1 -max i∈N e h       max e∈E i h |P i -Pe | |P i | + | Pe | + ǫ       ,
where P i is the pressure at the DOF i, Pe is the arithmetic average of the pressure on the element e, and ǫ = 10 -12 .

Improved discretization of the diffusion terms

The schemes ( 12) and ( 14) might have a high numerical error. Following the previous works [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF][START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF], in order to improve in the numerical schemes suppose to write the governing system of equations in the following form

       ∇• f a (u) -∇• Kw = 0 w -∇u = 0 ( 15
)
where for sake of simplicity the source term has been omitted, since it is irrelevant to the analysis. In the previous system, the additional variable w has been introduced for the gradients, in this way the original second order governing equations can be recasted as a first order system. The first order formulation for viscous problems has become a very standard approach in the DG framework, where numerical fluxes for the solution, u, and the gradient variable, w, are introduced to handle the discontinuities at interface of the elements [START_REF] Cockburn | The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems[END_REF]. Here, the discontinuity of the variable w at the faces of the elements is taken into account via the gradient reconstruction procedure. Consider a numerical scheme for the system (15) composed by a centered contribution plus a streamline stabilization term, namely

Ω e ψ i ∇• f a (u h ) -∇• (Kw) w -∇u h dΩ + Ω e A•∇ψ i T ∇• f a (u h ) -∇• (Kw) w -∇u h dΩ = 0, (16)
where

A•∇ψ i = A • ∇ψ i -K • ∇ψ i -∇ψ i 0 ,
and the scaling term is taken as follows

T = Ξ 0 0 Ξ d .
Supposing, now, that the gradient of the numerical solution has been reconstructed at each DOF, one can replace the equation on the second line of the system ( 16) with the approximation w ≃ ∇u h and consider only the first line, which now reads

Ω e ψ i ∇• f a (u h ) -∇• (K ∇u h ) dΩ + Ω e A•∇ψ i Ξ A•∇u h -∇• K ∇u h dΩ + Ω e K ∇ψ i • Ξ d ∇u h -∇u h dΩ = 0. ( 17
)
The first two integrals in the previous equation represent the discretization for the advection-diffusion problem by the means of a centered scheme plus a streamline stabilization term. The last integral represents an additional stabilization term, for the diffusive terms only, which vanishes in the advective limit, and the parameter matrix Ξ d is dimensionless.

It is interesting to note that the role of this additional term is to penalize the difference between the discontinuous and the reconstructed gradients, on each element. As already observed in the case of the Navier-Stokes equations [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF], taking Ξ d as the identity matrix is a good choice. Further studies on the effect of Ξ d on the numerical scheme are left for future work. With the Eq. ( 17) in mind, the schemes ( 12) and ( 14) are modified in order to include the extra stabilization term for the diffusive part of the equation. Namely, the linear scheme reads

Φ e i = Φ e N e dof + Ω e A •∇ψ i Ξ A • ∇u h -∇• K∇u h + S dΩ + Ω e K ∇ψ i • ∇u h -∇u h dΩ, (18)
while the non-linear scheme becomes

Φ e i = Φe i + ε e h (u h ) Ω e A•∇ψ i -∇• K∇ψ i Ξ A • ∇u h -∇• K ∇u h + S dΩ + Ω e K ∇ψ i • ∇u h -∇u h dΩ. ( 19
)

Gradient recovery strategies

As explained in Section 3, one has to assume that a continuous value of the gradient of the numerical solution is available on the faces of the elements. The strategy adopted here to achieve this goal consists in recovering the gradients at the DOFs, then the nodal values of the gradient are interpolated with Lagrangian basis functions. It is evident that the gradient reconstruction represents a key point to get higher-order accuracy. On the other hand, the computational effort in the gradient reconstruction should be as limited as possible.

In [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] different reconstruction strategies were considered, and it was found that the so-called Super-convergent Patch Recovery method proposed by Zienkiewicz and Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF] (SPR-ZZ) offers the possibility to reconstruct the gradients with the same order of accuracy of the solution, without spoiling the compactness of the numerical scheme. It was shown that the SPR-ZZ is much more accurate than other common approaches like Green-Gauss, L 2 -projection or least-square, which fail to reconstruct the gradient with more than second-order of accuracy.

Furthermore, in [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF] it is shown that if the gradient of the numerical solution is not reconstructed with same order of accuracy of the numerical solution, the overall accuracy of the numerical scheme is spoiled in regions where diffusive effects prevail over the advective ones or where advection and diffusion are equally important.

For sake of clarity Super-convergent Patch Recovery method is briefly recalled here.

Super-convergent Patch Recovery

The gradient of the numerical solution sampled at certain points on the element may have the same order of accuracy of the discrete solution [START_REF] Zienkiewicz | Finite Element Method (5th Edition) Volume 1 -The Basis[END_REF]. It can be shown that in the case of a segment element such particular points correspond to the Gauss-Legendre points [START_REF] Herrmann | Interpretation of finite element procedure as stress error minimization procedure[END_REF]. By tensor product such points can be defined also for quadrangles and hexahedrons. For triangles or tetrahedrons such property cannot be rigorously shown, but numerical experiments confirm it.

Using this fact, it is possible to compute a high-order accurate gradient within the whole element. Indeed, if at the super-convergent points the value of the gradient is accurate to order k + 1, by using a polynomial of degree k (the same of the solution) it is possible to obtain a high-order approximation everywhere within the element if this polynomial fits the super-convergent gradient values in a least square manner within a patch of the elements. Such a technique is called super-convergent patch recovery, introduced by Zienkiewicz and Zhu (SPR-ZZ) [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF].

Consider, for simplicity, a two dimensional scalar problem and assume that u h is the piecewise continuous polynomial interpolation of degree k of the solution u. The aim here is to reconstruct the gradient of the numerical solution at the DOFs, with (k + 1)-th order of accuracy. For each vertex of the grid a patch is defined as the union of all the elements that share that vertex, see Fig. 2. On the patch, the components of the reconstructed gradient are written in a polynomial form as follows ∂u h ∂x = p T a x , and

∂u h ∂y = p T a y , (20)
with p T (x) = (1, x, y, x 2 , . . . , x k , x k-1 y, . . . , y k ), a x = (a x 1 , a x 2 , . . . , a x m ) and a y = (a y 1 , a y 2 , . . . , a y m ).

On each element of the patch the gradient of the numerical solution has the highest accuracy at the super-convergent points. Thus, the high-order accurate expression of the reconstructed gradient on the patch is obtained by a least square fitting of the polynomial expressions [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF] to the values of the gradient at the super-convergent points of the patch. Assuming that N s i super-convergent points, x i,ℓ , ℓ = 1 . . . N s i , are available over the patch built around the vertex i, the values of the coefficients a x i and a y i are obtained by introducing the following functions

F x i = N s i k=1 ∂u h ∂x (x i,k ) -p T i,k a x i 2
and

F y i = N s i k=1 ∂u h ∂y (x i,k ) -p T i,k a y i 2 ,
with p i,k = p(x i,k ), and by minimizing them with respect to the polynomial coefficients, that is

∂F x i ∂a x i = 0 and ∂F y i ∂a y i = 0.
The minimization problem for the patch around the vertex i requires the solution of the following linear systems

A i a x i = b h x i and A i a y i = b h y i , (21)
in a least-square sense, where

b h x i =                                  ∂u h ∂x (x i,1 ) ∂u h ∂x (x i,2 )
. . .

∂u h ∂x (x i,N s i )                                  , b h y i =                                   ∂u h ∂y (x i,1 ) ∂u h ∂y (x i,2 )
. . .

∂u h ∂y (x i,N s i )                                   and A i =                       1 x i,1 y i,1 . . . y k i,1 1 x i,2 y i,2 . . . y k i,2
. . . . . . . . . . . . . . .

1 x i,N s i y i,N s i . . . y k i,N s i                       .
Once the vectors a x i and a y i are known, the reconstructed gradient ∇u h can be evaluated at any point within the patch around the vertex i.

The dimensions of the matrix A i are determined by the number of super-convergent points N s i and by the degree of the polynomials used to express the reconstructed gradient, that is A i ∈ R N s i ×m , where m is the number of the coefficients in the vectors a x i and a y i . The problem admits an unique solution if Rank A i = m, which is always satisfied in the case in which N s i ≥ m. Note that the SPR-ZZ method requires only the solution of small linear systems instead of a single large system as that obtained by a global L 2 projection. Note that, since the matrices A i depend only on the geometry, they are computed and factorized only once.

Generally, the number of elements contained in the patch is such that the condition N s i ≥ m is always satisfied, this means that the gradient reconstruction is compact because it involves only the elements contained within the support of a grid vertex. It is important to remark that the SPR-ZZ method solves for the vector of the polynomial coefficients, a x and a y over each patch, instead of solving for the nodal values of the reconstructed gradient as in the classical least-square gradient reconstruction. Once the polynomial coefficients are computed on the patch associated to a grid vertex, they can be used to evaluate the gradient at any point on that patch. As matter of fact, for the extra nodes introduced by the higher order approximation of the solution, the gradient is reconstructed by simply evaluating on the patch, at the coordinates of the nodes, the polynomial function constructed for the nearest grid vertex. Note that these nodes may belong to more than one patch, each of them is equally valid to evaluate the gradient. An arithmetic average is used to uniquely define the reconstructed gradient at these nodes.

The same idea is used to get the value of the reconstructed gradient at the boundary of the domain. In fact, for the nodes belonging to the boundary, the condition N s i ≥ m might not be satisfied without enlarging the stencil. To avoid the use of larger stencils for boundary nodes, it is possible to obtain the value of the reconstructed gradient at the boundary by evaluating, at the boundary node, the polynomial expansion already computed for the nearest domain vertex, i.e., the interior patch is used to evaluate the gradient at the boundary.

Another advantage of the SPR-ZZ approach is that being based on an element-wise evaluation (the super-convergent points inside the element of the patch where the gradient is evaluated), the stencil of the SPR-ZZ method is smaller than that required by the direct used of the least-square gradient reconstruction method, which is based on a nodal evaluation of the solution. In fact, it was observed in [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] that if the stencil used for the gradient reconstruction is made only by the nearest neighbor elements, the least square gradient reconstruction is at the most second order accurate.

For a quadrangle, the super-convergent points are defined uniquely. Considering a reference segment defined as x = [-1, 1], the super-convergent point is the point x = 0 in the case of a linear element, while in the case of a quadratic element the super-convergent points have coordinates ±1/ √ 3. The super-convergent points on the reference quadrangle are simply obtained by a tensorial product of the points defined on the reference segment.

For a linear triangle, the super-convergent point is the point with barycentric coordinates λ = (1/3, 1/3, 1/3), while in the case of a quadratic triangle the super-convergent points are not unique. In [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF], it was shown that the following choice of the super-convergent points coordinates

λ 1 = 2 3 , 1 3 , 1 3 , λ 2 = 1 3 , 2 3 , 1 3
, and

λ 3 = 1 3 , 1 3 , 2 3 .
is optimal, where λ 1,2,3 are the coordinates of the super-convergent points in the barycentric coordinate system.

Non-Linear LU-SGS Implicit Solver

The RANS equations are generally stiffer than the Navier-Stokes equations, in addition, since in turbulent flows the variation of the solution in the direction normal to wall is much higher than that in the tangential direction, very stretched grids are usually used in numerical simulations of turbulent flows. These factors, makes the solution of turbulent problems much more difficult than laminar ones and the robustness of the numerical solver is drastically reduced.

In [START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF], a matrix-free approach with the LU-SGS preconditioner was developed for the Navier-Stokes equations, showing very good performance in the simulation of laminar flows. Later, in this work it has been observed that the robustness and the effectiveness of the matrix-free solver are drastically reduced in the simulation of turbulent flows. In fact, this approach usually fails with a sudden breakdown of the solver or it stagnates at high values of the residual.

Since high-order methods are generally less robust than standard low-order methods the construction of a reliable implicit solver becomes crucial. For these reasons, instead of using a matrix-free based Newton-Krylov method a different strategy has been adopted for the solution of the non-linear system obtained from the discretization of the RANS equations.

Consider the solution of the non-linear system Eq. ( 11) with the Newton's scheme, which is written as

I ∆t n + ∂R ∂u ∆u n h = -R(u n h ), (22)
where R(u n h ) is the spatial residual of the numerical scheme at the time level n. The previous equation can be recasted in the following form

I ∆t n + ∂R i ∂u i ∆u n i = -R i (u n ) - j∈N i h j i ∂R i ∂u j ∆u n j , ∀ i ∈ N h , (23)
where N i h is the set of the DOFs contained in the stencil of DOF i. When the Gauss-Seidel method is applied to solve [START_REF] Herrmann | Interpretation of finite element procedure as stress error minimization procedure[END_REF], the forward and backward sweep steps read

I ∆t n + ∂R i ∂u i ∆u * i = -R i (u n ) - j∈N i h j<i ∂R i ∂u j ∆u * j , i = 1, 2, . . . , N dof , I ∆t n + ∂R i ∂u i ∆u n i -∆u * i = - j∈N i h j>i ∂R i ∂u j ∆u n j , i = N dof , . . . , 2, 1.
A better scheme can be obtained by applying a symmetric variation of the Gauss-Seidel method, in which both forward and backward steps are applied alternatively. The forward and backward sweep steps can be written in concise form as

I ∆t n + ∂R i ∂u i ∆u (k+1) i = -R i (u n ) - j∈N i h j i ∂R i ∂u j ∆u ( * ) j , ∀ i ∈ N h , ( 24
)
where k is an inner iteration index, the superscript ( * ) indicates the most recently updated solution, and ∆u (0) = 0. At the end, u n+1 = u n + ∆u (k+1) . The right-hand side of the previous equation can be further manipulated as follows. Denoting the latest available solution as u ( * ) = u n + ∆u ( * ) , it is possible to linearize R i u ( * ) as

R i u ( * ) ≈ R i u n + j∈N i h ∂R i ∂u j ∆u ( * ) j = R i u n + ∂R i ∂u i ∆u ( * ) i + j∈N i h j i ∂R i ∂u j ∆u ( * ) j ,
and the following relation is obtained

R i u n + j∈N i h j i ∂R i ∂u j ∆u ( * ) j = R i u ( * ) - ∂R i ∂u i ∆u ( * ) i .
Substituting the previous relation in Eq. ( 24), one obtains the following scheme

I ∆t n + ∂R i ∂u i ∆u (k+1) i - ∂R i ∂u i ∆u ( * ) i = -R i u ( * ) ,
which at last can be recasted as

I ∆t n + ∂R i ∂u i ∆u (k+1) i -∆u ( * ) i = -R i u ( * ) + ∆u ( * ) i ∆t n . ( 25
)
The Eq. ( 25) is solved with the forward and backward sweeps. At the beginning of each step the small diagonal blocks of the Jacobian matrix in the left-hand side can be inverted using the LU decomposition. The last term in the right hand side of the ( 25) is usually dropped for steady simulations, in order to accelerate the convergence rate. Note that the right-hand side of the previous equation is nothing but the residual evaluated with the latest available solution. When the symmetric Gauss-Seidel converges, one is actually solving the original set of equations instead of the linearized version [START_REF] Zienkiewicz | Finite Element Method (5th Edition) Volume 1 -The Basis[END_REF], for this reason the algorithm is called non-linear LU-SGS [START_REF] Sun | Efficient implicit non-linear lu-sgs approach for compressible flow computation using high-order spectral difference method[END_REF]. The expression [START_REF] Parsani | An implicit high-order spectral difference approach for large eddy simulation[END_REF] only requires the diagonal blocks of the Jacobian matrix. In some works [START_REF] Sun | Efficient implicit non-linear lu-sgs approach for compressible flow computation using high-order spectral difference method[END_REF][START_REF] Parsani | An implicit high-order spectral difference approach for large eddy simulation[END_REF][START_REF] Haga | An Implicit LU-SGS Scheme for the Spectral Volume Method on Unstructured Tetrahedral Grids[END_REF] the diagonal blocks of the Jacobian are approximated with a finite difference method, however this strategy has been found unsatisfactory for RANS simulations, thus the approximated diagonal blocks are computed directly.

The evolution of time step ∆t n used to dampen the solution, is controlled by the CFL number that is chosen according to the following law [START_REF] Tidiriri | Preconditioning techniques for the Newton-Krylov solution of compressible flows[END_REF] CFL

n = CFL n-1 R(u n-2 ) ∞ R(u n-1 ) ∞ ,
starting from a low CFL number, typically CFL 0 = 1 and limiting the maximum value of the CFL to CFL max = 1 000. The iterative process is stopped when the residual of the equations becomes small enough with respect to the initial value.

Numerical Results

The linear and non-linear RD schemes previously described are now used to perform numerical simulations of subsonic and transonic turbulent flows.

In all the simulations, on the wall hybrid weak-strong boundary conditions are applied: the velocity (v = 0) and the turbulent working variable (µ ⋆ t = 0) are strongly imposed on the wall, while the adiabatic condition is weakly enforced. The in/out-flow, far-field boundary and symmetric conditions are imposed in a weak way. The procedure used to impose the boundary conditions in the RD framework is detailed in previous works [START_REF] Abgrall | High order residual distribution scheme for rans equations[END_REF][START_REF] Abgrall | Linear and non-linear high order accurate residual distribution schemes for the discretization of the steady compressible navier-stokes equations[END_REF].

High Reynolds Number Flow over a Flat Plate

The two dimensional turbulent flow over a flat plate is here considered. The free-stream Mach number is M = 0.2 and the Reynolds number, computed with the free-stream conditions and an unity length of the plate, is Re = 5 × 10 6 , and the length of the plate is L = 2. The range of the computational domain in the x-direction is [-0.33, 2] with the leading edge of the flat plate at x = 0. The size of the computational domain in the y-direction is 1, which is about 66 times the boundary layer thickness at x = 2. At the inlet, the inflow boundary condition is imposed, while at the top and the exit, the static pressure is imposed. Along the plane y = 0, the symmetry boundary condition is imposed for -0.33 ≤ x ≤ 0 and the adiabatic no-slip wall boundary condition is imposed for 0 ≤ x ≤ 2. The free-stream value of the turbulent working variable is such that µ ⋆ t /µ ∞ = 3. A sequence of uniformly refined grids is considered. The coarsest grid is obtained by the triangulation of a 35 × 25 structured grid; the first number refers to the number of nodes on the horizontal boundaries, with 29 nodes along the plate, the second number refers to the number of nodes on the vertical boundaries. The grid is stretched in the wall-normal direction, and is also clustered near the plate leading edge; the height of the elements along the wall is 8.32 × 10 -6 . Finer grids are obtained by applying successive uniform refinements to the coarsest grid.

Except for the singularity of the solution at the leading edge of the plate, the flow is smooth and for this reason the linear scheme is used. For each grid, simulations with linear and quadratic approximation of the solution are performed, and the L 2 norms of the residuals of the mean flow and turbulence equations are reduced by ten orders of magnitude.

The skin friction distributions along the plate with linear and quadratic approximation of the solution, for grids with different resolutions, are reported in Fig. 3 as function of the quantity Re x = x v ∞ /µ ∞ . For reference, the experimental measures of Wieghardt [START_REF] Wieghardt | On the turbulent friction layer for rising pressure[END_REF] are also reported.

Note that the linear approximation of the solution gives a quite poor resolution of the friction coefficient along the plate and the solution is still far from a grid independent value, even for the finest grid. The quadratic approximation of the solution, on the other hand, gives very regular profiles of the friction coefficient. The values agree well with the experimental data and it is also clear that a reasonable accurate solution is already obtained on the coarsest mesh.

A quantitative comparison between the linear and quadratic approximations of the solution is reported in Fig. 4-a, where the values of the drag coefficient are reported as function of the number of DOFs. A reference value of the drag coefficient was obtained with a Richardson extrapolation [START_REF] Roache | Verification and Validation in Computational Science and Engineering[END_REF] of the drag coefficients computed with quadratic elements.

It is evident how, for the same number of degrees of freedom, the solution computed with the quadratic elements converges faster to the reference value, than the solution computed with the linear approximation. Thus, a second order scheme requires much more DOFs than the third order scheme to get a grid independent and accurate solution. This can be observed also in figure Fig. 4-b, where the error for the drag coefficient is reported as function of the number of the degrees of freedom. The error is computed as the absolute value of the difference between the numerical value of the drag coefficient and the reference one, i.e., err

CD = |c D h -c D ref |.
In Fig. 5 are also reported the profiles of the axial velocity in the boundary layer at the station x = 0.97, in terms of dimensionless quantities. The dimensionless distance from the plate is defined as y + = yρu τ /µ, with the friction velocity given by u τ = µ ∂v x ∂y wall , the dimensionless velocity is defined as u + = v x /u τ . The figures report also the law-of-wall profiles computed according to the following equations where κ = 0.41 is the von Karman constant and B = 5. Note that, since the friction velocity is based on the computed skin friction values at the wall, the dimensionless quantities depend on the accuracy of the solution. For linear approximation of the solutions, a good agreement with the theoretical profiles in the viscous sub-layer and in the log-layer regions is obtained. Also in the outer part of the boundary layer the velocity profiles dependent only a little on the grid resolution. Differences in the velocity profiles for quadratic elements can hardly be observed from the figure.

[31]            u + = y + , in the viscous sublayer u + = 1 κ log y + + B, in the loglayer (26)
In Fig. 6 are reported the profiles of the quantity χ = µ ⋆ t /µ at the station x = 0.97, to show the turbulent working variable along the boundary layer. Note that when the spatial resolution is poor (generally at the outer part of the boundary layer) negative values of the turbulent working variable might be obtained, however this effect is canceled out with the mesh refinement and/or the order increment of the numerical discretization.

Subsonic Flow over a RAE2822 Airfoil

The second test case considered is the subsonic flow over a RAE2822 airfoil at angle of incidence α = 2.79 • , the Mach number based on the free-stream conditions is M = 0.4 and the Reynolds number based on the free-stream conditions and the airfoil chord is Re = 6.5 × 10 6 . No-slip adiabatic boundary conditions are applied on the airfoil and far-field boundary conditions are applied on the outer boundary.

The problem is solved, with the linear scheme on a sequence of unstructured grids of triangles obtained with a sequence of uniform refinements of the coarsest grid shown in Fig. 7, which is made of 4 048 elements. As standard practice in this work, an iso-parametric representation is used for the solution and the geometry. Due to the strong stretching of the boundary layer elements, some of the quadratic elements have curved boundaries non only on the airfoil wall but also at the interior. Solutions are considered to be at convergence when the L 2 norms of the residuals of all the variables are reduced by ten orders of magnitude.

In Fig. 8 are reported the contours of the computed Mach number and of the turbulent variable, χ = µ ⋆ t /µ, for the third order simulation on a fine grid, and Fig. 9 shows the pressure and friction coefficients along the airfoil computed with linear and quadratic elements on two different grids, such that the number of DOFs of the two simulations is the same (N dof = 32 784); the respective iterative convergence histories are reported in Fig. 10.

For completeness, in Fig. 11 are reported the values of the lift and drag coefficients computed on different grids with linear and quadratic elements. Note how the use of quadratic elements allows a faster convergence with the mesh refinement of the lift and drag coefficients to their respective asymptotic values. In addition, Fig. 12 for the force coefficients; the error for the lift coefficient is computed as

err CL = |c L h -c L ex |,
where c L h is the value of the lift coefficient computed numerically and c L ex is a reference value obtained with a Richardson extrapolation [START_REF] Roache | Verification and Validation in Computational Science and Engineering[END_REF] of the lift coefficients computed with quadratic elements. The same procedure is also used to compute the error on the drag coefficient. For the lift coefficient, it is clear the advantage of using a high-order approximation over a second order approach. For the drag coefficient, the difference between second and third order schemes reduces, but nevertheless the quadratic approximation of the solutions gives a smaller error for the same number of DOFs.

Transonic Flow over a RAE2822 Airfoil

In this test case a transonic, turbulent flow over a RAE2822 airfoil is considered. The free-stream Mach number is M = 0.734, the angle of incidence and the Reynolds number are the same of the previous test case, namely α = 2.79 • and Re = 6.5 × 10 6 . In these conditions, the upper surface of the airfoil is characterized by the presence of a shock wave which impinges the boundary layer, generating a small recirculation bubble behind the shock. No-slip adiabatic boundary conditions are applied on the airfoil and far-field boundary conditions are applied on the outer boundary.

Due to the presence of the shock, the problem is discretized with the non-linear scheme. The L 2 norms of the residuals of the mean flow and of the turbulence equation are reduced at least by nine orders of magnitudes.

The same type of grids used in the previous test case is adopted here and three levels of refinement are considered. In Fig. 13 are reported the Mach number contours for the solution on the fine grid (64 768 elements) with quadratic elements, together with the iterative convergence history. Note that the shock is sharply captured with no oscillations.

In Fig. 14 are shown the pressure and the friction coefficients along the airfoil with linear elements on a finer grid and quadratic elements on a coarse grid, such that the number of DOFs in the two cases is the same: 32 784. The values of the computed pressure coefficient agree very well with the experimental data, especially with the quadratic discretization of the solution. The leading edge suction peak of the pressure is slightly under-predicted by the numerical simulations, due to the fact that the simulations have been performed in a fully turbulent regime. Also the friction coefficient agrees well with the experimental data. On the upper part of the airfoil, note the negative values of the friction coefficient just behind the shock obtained with the quadratic approximation of the solution, that identify the separation bubble induced by the shock.

Finally, in Fig. 15 are reported the values of lift and drag coefficients for linear and quadratic approximation of the solution with three levels of grid refinement.

High-lift Multi-element Airfoil Configuration L1T2

The turbulent flow over the AGARD L1T2 high-lift multi-element airfoil is now considered. The geometry consists of a three-element airfoil configuration, the free-stream Mach number is M = 0.197, the angle of incidence is α = 20.18 • , and the Reynolds number based on the free-stream conditions and the retracted airfoil chord is Re = 3 520 000. No-slip adiabatic boundary conditions are applied on the airfoil and far-field boundary conditions are applied on the outer boundary.

The problem is solved with the linear scheme on a mesh of 33 338 triangles, shown in Fig. 16. This test case is particularly challenging due to the very high angle of attach and the flow complexity. The solution is considered to be at the steady state when the L 2 norms of the residuals of all the equations are reduced at least by nine orders of magnitude. In Fig. 17 is reported the converged history of the simulations with linear and quadratic elements. Note that order sequencing has been used to accelerate the convergence of the high-order simulation, i.e., the solution computed with linear elements is used as initial solution for the simulation with quadratic elements. Despite the complexity of the flow-field, the non-linear LU-SGS scheme successfully converges to the steady state with both linear and quadratic elements.

In Fig. 18 are reported the Mach number contours computed with linear and quadratic elements. Clearly the high-order discretization drastically improves the accuracy of the solution compared to the second-order simulation. In Fig. 19 are reported the computed pressure and skin friction coefficients over the elements of the airfoil. The computed pressure coefficient is also compared against the experimental data [32]: a very good agreement between and numerical experimental data con be observed. For sake of completeness, in Table . 1 are listed the computed lift and drag coefficients, together with some values reported in literature.

ONERA M6 wing

The simulation of the transonic turbulent flow over the ONERA M6 wing is performed. The free-stream Mach number is M = 0.8395, the angle of incidence is α = 3.06 • and the Reynolds number based of the free-stream values and the mean aerodynamic chord is Re = 11.72 × 10 6 . The grid is composed by 123 444 tetrahedrons, and the mean height of elements on the wing in dimensionless units is approximately y + 1 = 5. Far-field boundary conditions are applied on the outer part of the domain, symmetry boundary conditions are applied on the vertical plane that intersects the wing root and the wing surface is modeled as a non-slip adiabatic wall. In Fig. 20 are reported the grid used for the simulation and the pressure coefficient contours for the third-order simulation. [START_REF] Burgess | High-order Discontinuous Galerkin Methods for Turbulent High-lift Flows[END_REF] 265 311 4.055891 0.067782 DG P2, SA [START_REF] Burgess | High-order Discontinuous Galerkin Methods for Turbulent High-lift Flows[END_REF] 553 707 4.049804 0.066554 DG P3, SA [START_REF] Burgess | High-order Discontinuous Galerkin Methods for Turbulent High-lift Flows[END_REF] 945 930 4.048049 0.066471 DG P4, SA [START_REF] Burgess | High-order Discontinuous Galerkin Methods for Turbulent High-lift Flows[END_REF] 1 441 980 4.047551 0.066402 FV, k-ω EARSM [START_REF] Rumsey | Prediction of high-lift flows using turbulent closure models[END_REF] 114 908 4.08 0.068 Experiment [START_REF] Research | A Selection of Experimental Test Cases for the Validation of CFD Codes[END_REF] -4.11 0.068 The non-linear solver has been used to discretize the governing equations. Convergence is considered achieved when the L 2 norms of the residuals of all the equations are reduced by eight orders of magnitude. Simulations have been performed in parallel on 64 processors.

In figure Fig. 21, the values of the pressure coefficient for the third-order simulation at different stations along the span-wise direction of the wing are compared against the experimental data. Despite the quite coarse mesh, a good agreement with the experimental data is obtained and even the lambda shock is well represented. However, a finer grid would be required to better capture the shock structure at wing tip and also to better resolve the complex shock-boundary layer interaction occurring here.

NASA Delta wing

The simulation of the subsonic turbulent flow over the NASA 65 • sweep delta wing with a medium radius leading edge is now considered. The free-stream Mach number is M = 0.4, the angle of attach is α = 13.3 • and the Reynolds number based on the free-stream conditions and the mean aerodynamic chord is Re = 3 × 10 6 . This test case has been computed with the linear scheme using linear and quadratic approximation of the solution, on a grid of 1 145 797 tetrahedrons (Fig. 22-a). For this test case only a linear approximation of the geometry is available. Simulations has been performed in parallel using 96 cores. With linear elements the L 2 norms of the residuals have been reduced by 8 orders of magnitude, while with quadratic elements they could non be reduced by more than 5 orders of magnitude. This is issue is probably due to to the poor representation of the wing geometry.

In Fig. 22-b are shown, for the solution with P2 elements, the contours of the pressure coefficients and slices of the turbulent working variable at different stations along and behind the wing. Starting from the apex of the wing, the flow around the rounded leading edge is initially attached. Further downstream the flow separates and the inner vortex is generated, eventually a separation of the shear layer occurs further downstream at the leading edge and the outer vortex is generated. The resulting vortex structure is then transported behind the wing.

The comparison with the experimental results [START_REF] Chu | Experimental surface pressure data obtained on 65 • delta wing across Reynolds number and Mach number ranges[END_REF] is shown in Fig. 23. The agreement between the numerical results and the experimental data is mildly good, considering the quite coarse mesh and the simple turbulent model used. Note also that capturing the correct position of the vortices separation and reattachment points is quite challenging for this test case, as observed in other numerical simulations [START_REF] Schutte | Numerical investigations on the VFE-2 65-degree rounded leading edge delta wing using the unstructured DLR TAU-Code[END_REF].

Conclusions

This paper presents a high-order Residual Distribution discretization method for the solution of the multidimensional RANS equations. A modified form of the Spalart-Allmaras turbulence model has been introduced, which is strongly coupled with the mean flow equations. The numerical scheme relies on a continuous approximation of the solution and its gradient. A key aspect that emerged from previous works is that the reconstructed gradient of the numerical solution must have the same accuracy of the discrete solution. The so-called Super-convergent Patch Recovery method has been used to obtain accurate values of the gradient of the numerical solution, keeping the computational stencil compact.

The possibility to construct non-linear schemes within the Residual Distribution framework makes these schemes very attractive, because shocks and smooth solution can be handled within the same numerical scheme without introducing any tuning parameter or additional stabilization terms. On the other hand, linear schemes can be still useful when smooth solutions are expected.

The non-linear LU-SGS solver has been used to construct a robust implicit scheme which has been successfully used to make the scheme converge to the steady state for both the mean flow equation and the turbulent model. The accuracy of the numerical schemes has been tested with standard test cases: the turbulent flow over the flat plate, the subsonic and transonic flow over the RAE2822 airfoil. The proposed numerical approach has been subsequently used to perform more challenging simulations like the two-dimensional flow over the high-lift multielement airfoil and three-dimensional cases like the the transonic flow over the M6 wing and the subsonic flow over a delta wing.

Figure 1 :

 1 Figure 1: Comparison between ψ µ ⋆ t (solid line) and χ (dashed line) as function of µ ⋆ t /µ.

Figure 2 :

 2 Figure 2: Interior super-convergent patches for quadrilateral and triangular elements: (a-b) linear elements, (c-d) quadratic elements. The symbols (•) indicate the patch assembly points, the symbols (•) indicate the points where the gradient is reconstructed and the symbols (△) indicate the super-convergent points.

Figure 3 :

 3 Figure 3: Friction coefficient along the plate for the turbulent flow over a flat plate, with different grid resolutions and with linear (a) and quadratic (b) approximation of the solution.

Figure 4 :

 4 Figure 4: Effects of the grid and polynomial refinement on the solution accuracy for the turbulent flow over a flat plate: (a) drag coefficients (b) error on the computation of the drag coefficient.

Figure 5 :

 5 Figure 5: Dimensionless axial velocity profiles in the boundary layer at the point x = 0.97 for the turbulent flow over a flat plate, with different grid resolutions and with linear (a) and quadratic (b) approximation of the solution.

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Dimensionless turbulent working variable profiles in the boundary layer at the point x = 0.97 for the turbulent flow over a flat plate, with different grid resolutions and with linear (a) and quadratic (b) approximation of the solution.

Figure 9 :Figure 10 :

 910 Figure 9: Pressure (a) and friction (b) coefficient along the RAE2822 airfoil with linear and quadratic elements on two different grids, such that N dof = 32 784.

Figure 11 :

 11 Figure 11: Values of the lift (a) and drag (b) coefficients as function of the number of the degrees of freedom, for the subsonic turbulent flow over the RAE2822 airfoil, with linear and quadratic elements.

Figure 12 :

 12 Figure 12: Errors with linear and quadratic elements on the lift (a) and drag (b) coefficients for the subsonic turbulent flow over the RAE222 airfoil.

Figure 13 :

 13 Figure 13: Mach number contours (a) and iterative convergence history (b) for the transonic turbulent flow over a RAE2822 airfoil computed with quadratic elements.

Figure 14 :Figure 15 :

 1415 Figure 14: Pressure (a) and friction (b) coefficients over the RAE2822 airfoil, for the transonic turbulent flow computed with linear and quadratic elements on two grids, such that N dof = 32 784.

Figure 16 :

 16 Figure 16: Global view (a) of the computational mesh around the L1T2 multi-element airfoil, and zoom (b) near the slat and the main element.

Figure 18 :

 18 Figure 18: Mach number contours for the flow over the L1T2 multi-element airfoil with linear (a) and quadratic (b) elements.

Figure 19 :Figure 20 :

 1920 Figure 19: Pressure (a) and skin friction (b) coefficients for the flow over the L1T2 multi-element airfoil.

Figure 21 :Figure 22 :

 2122 Figure 21: Pressure coefficient distribution at different span-wise locations over the M6 wing, for the third-order simulation.

Figure 23 :

 23 Figure 23: 3D turbulent subsonic flow over a delta wing. Pressure coefficients distribution at different wing sections. Comparison between the quadratic solution and the experimental data.

Table 1 :

 1 Computed force coefficients and reference values reported in literature for the simulation of the flow over the L1T2 multi-element airfoil.


Jean-Charles Bricola

Michel Bilodeau

Serge Beucher

DEPTH MAP ESTIMATION: A REGION-BASED APPROACH

This paper presents an approach to the estimation of depth from stereo images which exploits correspondences of image segments. A key motivation is to bypass pixel correspondences which are ambiguous in texture-less regions. One of the main difficulty though is to come up to an equivalent partitioning of both stereo images that facilitates region matching. The reference image which the depth is estimated for is first segmented using the watershed algorithm. Regions are then transferred to the other image of the stereo pair according to a correlation analysis of region contours. This transfer yields regional disparities which are significant as long as the contours for a given region are not due to an occlusion. Based on that observation, we propose an algorithm that detects occlusion contours along a region and that rectifies disparity accordingly.

Introduction

This article focuses on stereoscopic configurations that involve two images of a scene taken at different viewpoints. A depth map is computed for one image chosen from the stereo pair, namely the reference image. We refer to the other image as the second image. For a pixel in the reference image, there exists an epipolar line in the second image along which a correspondence may be found, provided the transformation between the coordinate systems of the stereo cameras is composed of a non null translation and no occlusion prevents the scene point from being projected in both stereo images [START_REF] Hartley | Multiple view geometry in computer vision[END_REF][START_REF] Horaud | Vision par ordinateur: outils fondamentaux, chapter Vision stéréoscopique[END_REF]. In most scenarios, rectified images will be considered, for which a pixel and its correspondence always share the same ordinate. The disparity is the measured offset between the pixels abscissa and is inversely proportional to the depth we aim at estimating. The problem of depth map computation from stereo dates back to the 70s and has until now relied mostly on pixel correspondences between stereo images. Early approaches considered feature points and contour matching in that respect [START_REF] Jo | Computer Matching of Areas in Stereo Images[END_REF][START_REF] Ohta | Stereo by intra-and inter-scanline search using dynamic programming[END_REF]. Whilst the established correspondences were reliable, the measurements suffered from sparsity. Current methods on the contrary attempt to find a match for every pixel of the reference image. The correspondence is typically selected according to the smallest dissimilarity cost measured between image patches centred at the considered pixels. An exhaustive list of such measures is given in [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. Because of matching uncertainties in homogeneous regions, the resulting disparity maps usually require a post-processing step that consists of minimizing an energy function. The latter usually leads to smooth transitions of disparities among neighbouring pixels except on edges where discontinuities are permitted [START_REF] Fua | A parallel stereo algorithm that produces dense depth maps and preserves image features[END_REF] whilst ensuring that the disparities still agree with the reliable initial measurements.

The use of regions within depth estimation algorithms becomes increasingly popular and differentiates into super-pixels and layer-based categories. For instance, in [START_REF] Zitnick | Stereo for image-based rendering using image oversegmentation[END_REF], the reference image is over-segmented and disparities are estimated at the segment level whilst smoothness is enforced across neighbouring segments having close colour distribution. In [START_REF] Bleyer | A layered stereo matching algorithm using image segmentation and global visibility constraints[END_REF] however, disparities are computed across image layers for which planar equations are estimated. Layers are obtained by means of colour segmentation and regions belonging to the same disparity plane according to an initial disparity map are merged. Our work falls into the layer-based category. We propose an algorithm that relies on the disparity information along region contours and inside subregions to estimate depth maps. The gradient and markers generation required for the watershed segmentation of the reference image is first detailed in section 2. Section 3 introduces the concept of regional disparities. We show that the latter can be used to produce markers suitable for segmenting the second image of the stereo pair and producing instantaneous region matchings. We also discuss their relevance with respect to stereopsis. Section 4 explains how to overcome the shortcomings of regional disparities computed for initial segments and presents an algorithm that correct most disparity errors due to occluding contours. Results are finally exposed and evaluated in section 5.

Reference Image Segmentation

The watershed algorithm [START_REF] Beucher | Segmentation d'Images et Morphologie Mathématique[END_REF] partitions an image into a series of connected components satisfying each an homogeneity criterion and can be implemented effectively using hierarchical priority queues [START_REF] Beucher | Hierarchical queues: general description and implementation in mamba image library[END_REF]. It has moreover the ability to be controlled by markers. This remarkable property enables the transfer of markers from the reference to the second image of the stereo pair, hence producing equivalent stereo partitions together with their matchings as explained in section 3. The watershed-based segmentation constitutes, therefore, a very suitable choice with respect to our application. In short, watersheds are obtained by flooding a topological surface, typically the colour gradient of the image, from a set of markers. Each marker gives birth to a lake which a unique label is assigned to. Each lake in turn will lead to an image region. As the flooding goes on and uniformly increases in altitude, watersheds are constructed as soon as two lakes with different labels are on the verge of merging. More precisely the geodesic skeleton by influence zone of the lakes just before merging is computed inside the new lake obtained just after merging. The segmentation result depends on the choice of appropriate markers. As far as the reference image is concerned, markers are chosen from the h-minima of the gradient in order to capture the regions of interest.

Colour and thick gradients

Colour is taken into account within a watershed segmentation by computing a colour gradient. Our current implementation just considers the supremum of gradients for the red, green and blue channels. However, there exists more meaningful representations which are worth mentioning [START_REF] Angulo | Morphologie mathématique et indexation d'images couleur: application à la microscopie en biomédecine[END_REF]: in particular, the linear combination of the hue and luminance gradients according to the saturation could be used as a colour gradient instead as well as perceptual differences in CIE-La*b* space. Transitions between two objects which span several pixels are captured by computing a thick gradient. The gradient values have indeed a significant impact when extracting markers from the h-minima of the gradient and when a hierarchical transform is applied on the watersheds. The thick gradient of size λ is computed Symbol Description Morphological operators B and H are structuring elements H is an isotropic structuring element of elementary size

δ B (f ) Dilation of f : δ B (f ) [x] = sup x∈B f [x] ε B (f ) Erosion of f : ε B (f ) [x] = inf x∈B f [x] D 1 g (f ) Geodesic dilation of f under mask g: D 1 g (f ) = inf (δ H (f ) , g) R g (f ) Geodesic reconstruction of g from marker f : R g (f ) = D +∞ g (f ) = D 1 g . . . D 1 g (f ) . . . R * g (f ) Dual geodesic reconstruction of g from marker f : R * g (f ) = -R -g (-f )
Table 1: Notation according to equation 1:

g λ (f ) = (δ λH (f )) -(ε λH (f )) (1)
Hence, if a transition between two objects requires λ pixels to complete, the thick gradient g λ at the middle of the transition equals the difference between the smallest and the highest values of the transition. In comparison with the morphological gradient of size λ = 1, pixels on the edges are assigned a gradient magnitude similar to the one expected for a perfect transition between two regions whilst the magnitude remains small for irrelevant variations of intensities. Figure 1(b) shows an example of thick gradient of size λ = 4 computed for an image composed of 1920×1080 pixels.

Markers from h-minima

Markers control the behaviour of the watershed segmentation: ideally they should lie at the bottom of the topological surface being flooded and discriminate every region to be segmented as well. To this end, we use the h-minima of the gradient [START_REF] Beucher | Maxima and minima: a review[END_REF]. Let g be the gradient of the reference image and h a positive constant. The h-minima are computed according to equation 2.

min(g, h)[x] = 1 if R * g (g + h) -g [x] > 0 0 otherwise (2)
In other words, x is a h-minima of g if and only if there exists on g a strictly descending path from x to a minima of g, say m, such that g(x)g(m) > h. By definition, h-minima contain the local minima of the function as well. h is chosen to be sufficiently large so as to ignore small catchment basins that are due to noise in the gradient. This choice implies that h-minima occasionally overlap contours at low gradient values.

In order to solve that problem, we compute a distance function p by construction of successive erosions applied over the binary mask corresponding to the h-minima. Thin structures appearing on the mask will appear as valleys in function p. New markers are given by the connected components resulting from the indicator function p -R p (αp) > 0 for 0 α < 1. Small values of α encourage markers splitting at very narrow valleys only and reduce slightly the overall size of the markers, getting them away from the borders as shown in figure 1(c). Note that contrary to an erosion or an opening, that approach to process h-minima preserves small markers.

gions occluded in the second image as well as a correlation measure based on the intersection volume between matched regions.

Disparity refinement

The refinement algorithm proposed below aims at detecting occlusion contours and rectifying mistaken regional disparities accordingly. It first assigns to every region two disparities: one that is representative of the disparity near the left contour of the region, the other of the right contour. To do so, regions obtained for the reference image are vertically split into left and right subregions and regional disparities for these subregions are estimated. Second, neighbourhood relationships between regions are defined: N l (R i ) includes all regions R j ∈ R that are adjacent to the left subregion of R i and the same goes for N r (R i ) with respect to the right subregion of R i . Third, algorithm 1 is executed for every region R i :

Algorithm 1 Regional disparity refinement 1: function REFINEDISPARITIES(R, i) 2: d l ← getLeftContourDisparity(R i ) 3: d r ← getRightContourDisparity(R i ) 4: if d l > d r + τ then 5:
for all R j ∈ N l (R i ) do

r ← getRightContourDisparity(R j )

7:

if |d l -d (j) r | < τ then 8: d l ← d r 9:
end if for all R j ∈ N r (R i ) do

l ← getLeftContourDisparity(R j )

14:

if |d r -d (j) l | < τ
return (d l , d r ) 20: end function
The intuition behind this algorithm is the following: when a region is assigned to very different left and right disparity values (i.e. when algorithm condition 4 or 11 is fulfilled for a given threshold τ ), either the highest disparity value is due to an occluding contour or the region is simply not fronto-parallel. To verify whether we fall into the first category, we check whether condition 7 or 14 is true. In the affirmative, we assign to the entire region the disparity that is the smallest between the left and right disparities. Sample results are shown in figures 2(c) and 2(g). The second image markers (e) are obtained by translating each region from the reference image according to its regional disparity and hence produce an equivalent segmentation (f) for the second image.

Results and Discussion

Figures 1(a-c) show the main stages of the watershed segmentation procedure for a reference image (namely colour gradient computation and markers generation). Figure 1(d) represents the resulting partitioning of the reference image. Thanks to the regional disparities computed for that image, markers are generated in order to capture the regions of the second image, i.e. they fall into the catchment basins of the second image gradient (figure 1(e)). Since each marker is assigned to the label of the reference image region from which it originates, we obtain an equivalent labelled partitioning (figures 1(d) and 1(f)) from which region matchings are obtained. Figure 2 presents brute regional and refined disparity measurements mapped onto the regions of the reference image. It clearly appears that refined disparities give more insight on the stereopsis than regional disparities from initial image segments alone. The mistaken assignments are mainly due to occluding contours: on figure 2(f), one can see that the disparity value assigned to the pig has been transferred to the surrounding region. The same goes for the cone as shown in figure 2(b). The disparity refinement algorithm splits those ambiguous regions in two and discovers brutal changes of disparity values for the left and right resulting subregions. Since one of the subregion seems to inherit from the disparity value of a neighbour region, the disparity of that subregion is hence corrected using the disparity found for the other subregion. Figures 2(c) and 2(g) result from the refinement algorithm. Considering direct neighbourhood only is however a strong limitation of the refinement algorithm as noticed for the region between the jar and the brush shown in figure 2(c): since that region is homogeneous and completely surrounded by occlusion contours, both left and right subregions have indeed invalid disparities. Future work on our framework will aim at considering more adequate models to overcome this kind of issue.

Conclusion

We have presented an approach that exploits regions and the disparities along their contours to estimate depth maps from stereo. The proposed algorithm relies on the concept of regional disparities and yields visually acceptable results with respect to ground truth disparities. Still it only partially solves two problems arising out of this approach: the need to accurately detect contours which are due to occlusions and a way to compute the evolution of disparities across non fronto-parallel surfaces. Our future work is bound to concentrate on the resolution of these shortcomings.

Figure 1 :

 1 Figure 1: Watershed-based segmentation of stereo images. From the reference image (a) a colour thick gradient is computed (b). The h-minima of the gradient are rescaled using a distance function and serve as markers (c) used to drive the watershed segmentation of the reference image (d).The second image markers (e) are obtained by translating each region from the reference image according to its regional disparity and hence produce an equivalent segmentation (f) for the second image.

Figure 2 :

 2 Figure 2: Disparity maps estimated over Middlebury stereo datasets art and dolls [13]. (a) and (e): Reference images. (b) and (f): Regional disparity maps. (c) and (g): Refined disparities mapped onto image subregions. (d) and (h): Ground truth disparities.

  then

	15:	d r ← d l
	16:	end if
	17:	end for
	18:	end if
	19:

Acknowledgments

This work has been performed in the project PANORAMA, co-funded by grants from Belgium, Italy, France, the Netherlands, and the United Kingdom, and the ENIAC Joint Undertaking.

Region transfer and disparities

In this section, we explain how regional disparities are estimated for every region discovered in the reference image. We also discuss their relevance and usefulness with respect to depth estimation and equivalent stereo segmentation.

Regional disparities

Let R = {R 1 , • • • , R n } be the set of partitions obtained for the reference image. R i (x, y) is an indicator value that equals 1 for every pixel (x, y) belonging to region R i , 0 otherwise. Also let g R and g S be the colour gradients for the reference and second images respectively. For each region R i , the translation t (i) that provides the best overlap between the gradients of the reference image and the shifted second image within R i is estimated. As stereo images are assumed to be rectified, the translation only occurs along the x-axis, so t (i) = (d (i) , 0). Optimal translations are determined by equation 3.

From equation 3, we finally obtain for each region R i its regional disparity |d (i) |.

The sum of absolute differences (SAD) is considered here as a cross-similarity measure and the gradients as inputs for that measure. It is of course possible to take the luminance or a hue channel in place of the gradient. Doing so is preferable when source images have no illumination changes and do not suffer from noise as this will increase the support, hence the accuracy of the correlation. Whilst regional disparities are relevant to fronto-parallel regions, they only provide an average estimate of the disparities among other regions. Furthermore regional disparities are influenced by contour disparities themselves. This is a shortcoming with respect to regions for which contours are due to occlusions as shown in figures 2(b) and 2(f). These observations justify the disparity refinement step exposed in section 4.

Equivalent stereo segmentation

Estimating regional disparities and generating markers for segmenting the second image of the stereo pair are interestingly two equivalent problems. A simple method for obtaining an equivalent segmentation of the second image alongside a matching between regions is described below. Let M i be an indicator function representing the pixels of the second image which are spanned by the marker corresponding to region R i in the first image. Then:

Each marker is resized using a distance function as described in section 2 with the exception of markers labelling which remains unchanged, even if splitting occurs. Finally, the watershed segmentation driven by these markers and the the second image colour gradient is performed. Regions inherit from the marker labels which yields region matchings between the two images. These matchings potentially provide a mean to detect re-