R Abgrall

D De Santis

High order preserving residual distribution schemes for the laminar and turbulent Navier Stokes on arbitrary grids

This paper deals with the construction of a class of high order accurate Residual Distribution schemes for the Navier Stokes equations using conformal meshes. The approximation of the solution is obtained using standard Lagrangian finite elements, and the total residual of the problem is constructed taking into account both the advective and the diffusive terms in order to discretize within the same scheme both parts of the governing equation. To cope with the fact that the normal component of the gradients of the numerical solution is discontinuous across the faces of the elements, the gradient of the numerical solution is recovered at each degree of freedom of the grid and then interpolated with the same shape functions used for the solution. The procedure is fully described for the scalar case, and formaly extended to the system case. Linear and non-linear schemes are constructed and their accuracy is first tested with the help of manufactured solutions, and then applied to several (2D and 3D) test cases.

I. Introduction

In the last years different high order schemes have been developed to obtain a higher order (more than two) discretization of the Navier-Stokes equations. One of the most attractive method seems to be the discontinuous Galerkin (DG) scheme. [START_REF] Cockburn | Discontinuous Galerkin methods: theory, computation and application[END_REF] Residual Distribution (RD) schemes [START_REF] Struijs | Fluctuation splitting schemes for the 2d Euler equations[END_REF][START_REF] Abgrall | Toward the ultimate conservative scheme: Following the quest[END_REF][START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF] represent a very interesting alternative to DG schemes. While computationally compact and probably more flexible, DG schemes suffer from the serious drawback of a very fast growth of the number of degrees of freedom (DOF) with the cell polynomial degree. In RD schemes the formulation remains local, as in DG, but the number of DOFs growths less quickly because the numerical solution is assumed to be continuous. Furthermore, in the case of discontinuous solutions, the non oscillatory properties of the numerical scheme are probably better understood for RD scheme than for DG methods.

RD schemes have been developed mainly for advection problems due to possibility to construct multidimensional upwind schemes which guarantees a small discretization error compared to the standard Finite Volume schemes, but the discretization of advection-diffusion problems with the RD schemes is still an open problem. One of the main issue concerns the possibility to take into account within the same scheme advective terms, by the means of upwind mechanism, and diffusive phenomena, which on the other hand have an isotropic behavior. To address this problem mixed upwind/central schemes have been developed, in which RD methods for the advection terms are combined with central schemes, usually based on the Galerkin 1 discretization of the diffusion terms. For such type of schemes a proper blending between the RD and the Galerkin schemes must be constructed otherwise the accuracy of the resulting schemes is spoiled. [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF] The approach used in this work is based instead on the construction of a RD method in which the advection and the diffusion are handled within the same scheme. Unfortunately this introduces a new complication because, for piecewise polynomial 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 considering a numerical flux for the viscous term along the faces of the elements, as happens in the DG or Finite Volume schemes, the approach adopted in this work consists in recovering a unique set of values for the gradient of the numerical solutions at each degree of freedom (DOF) of the grid. Then, these values are interpolated with the same continuous functions used to interpolate the solution. It is evident that gradients have to be recovered with higher order of accuracy to construct a high order scheme: the crucial point is the strategy used to recover the gradients at the DOFs. The problem of the gradient recovery has been addressed in [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] together with the construction of accurate and robust, linear and non-linear RD schemes.

Here, the procedure of the gradient reconstruction is briefly recalled, and results about the accuracy of the reconstruction are reported for a scalar problem. The method is formally extended to the Navier-Stokes and RANS equations.

The structure of the paper is the following. In Section II, starting from the scalar advection problem, the basic ideas of the RD methods are recalled. The issues related to the discretization of advection-diffusion problems in the RD framework are discussed in Section III, while in Section IV is described in detail the construction of a class of linear and non-linear RD schemes. In Section V is evaluated the gradient recovery strategy we propose. In Section VII the accuracy of the proposed numerical schemes is tested through the use of the manufactured solution, and then the numerical method is used to perform simulation of two and three dimensional laminar and turbulent steady flows.

II. Basics of the residual distribution method

In this section, the main idea of the RD method for the scalar advection problems is briefly recalled, and the fundamental properties of conservation, consistency and accuracy of the numerical scheme are also reported. They extend immediately to systems with straightforward adaptations. Furthermore, the notation used through the paper is introduced.

Consider the steady conservation law for the quantity u

∇•f (u) = 0, (1a)
where f (u) ∈ R d is a given flux function of the unknown u(x) ∈ R, x ∈ Ω ⊂ R d , with d the number of the spatial dimensions (here d = 2 or d = 3). The Eq. (1a) must be supplemented with the proper boundary conditions on some part Γ of the boundary ∂Ω

B(s, u) = g(s), s ∈ Γ, (1b)
where the function g is known abd B is some functional. The domain Ω is discretized with N e non-overlapping elements with characteristic length h, the set of all the elements a is denoted by E h , the list of the DOFs is denoted by Σ h , the set of all the boundary faces is denoted by F h b , 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 accross the faces of the elements. If the standard Lagrangian shape functions are used, the approximated solution u h can be written as u h (x) = i∈Σ h ψ i (x) u i , for x ∈ Ω, with u i the numerical solution at the generic DOF i, and ψ i the Lagrange basis function at the DOF i.

a assumed to be polyhedral b a boundary face f is the intersection of an element e and ∂Ω. The mesh is assumed such that a boundary face f = e ∩ ∂Ω is also a face of e.

The approximated solution will give rise to a residual on each element, namely

Φ e u h = Ωe ∇•f u h dΩ = ∂Ωe f u h • n d∂Ω.
The integral quantity Φ e u h is called the total residual of the element e. Similarly, for any element f on the boundary ∂Ω of Ω, one can define a total residual

Φ f u h = f f u h • n -F u h , g d∂Ω,
where F is a numerical flux consistant with (1b).

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

Φ e i = β e i u h Φ e u h , ∀i ∈ Σ e h ,
where Σ e h is the list of the DOFs of the element e and β e i are the distribution coefficients, which can be in general function of u h . A similar relation is written on the boundary residuals

Φ f i = β f i u h Φ f u h ∀i ∈ Σ f h ,
where Σ f h is the list of DOFs of the boundary face f . It is easy to see 8 that the following conservation constraints must be satisfied for any element e and any boundary face

f i∈Σ e h Φ e i = Φ e ∀e ∈ E h , and
i∈Σ f h Φ f i = Φ f ∀f ∈ F h ,
in addition to the standard assumptions of the Lax-Wendroff theorem in order to garanty that the limit solution, if it exists, is a weak solution of (1a). To obtain an equation for each nodal value,the following relations are written for each DOF

e∈E h,i Φ e i u h + f ∈F h,i Φ f i u h = 0, ∀i ∈ Σ h , (2)
where E h,i (resp. F h,i ) is the set of the elements (resp. faces) 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 an 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 h,i Φ e i u h = 0, ∀i ∈ Σ h , (3)
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.

A. Consistency and accuracy

The fundamental properties of consistency and accuracy for RD schemes have been analyzed in [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF] and are briefly reported here for sake of completeness.

Assuming that a sequence u h is bounded in L ∞ when h → 0 and if exist w, such that u h → w in L 2 when h → 0, then w is a weak solution of (1). In the proof, the continuity of the interpolant across the faces is assumed, although this constrain may be alleviated and RD schemes with discontinuous elements can be constructed. [START_REF] Abgrall | Development of residual distribution schemes for the discontinuous Galerkin method: the scalar case with linear elements[END_REF][START_REF] Abgrall | A residual distribution method using discontinuous elements for the computation of possibly non smooth flows[END_REF][START_REF] Hubbard | A framework for discontinuous fluctuation distribution[END_REF] To analyze the accuracy of the RD the following truncation error is introduced, for any C 1 compactly supported function ϕ,

E u h , ϕ = i∈Σ h ϕ(x i ) e∈E h,i Φ e i + f ∈F h,i Φ f i . (4)
If the solution u is smooth enough and the residuals, computed with the numerical solution u h , are such that

Φ e i = O h k+d , ∀e ∈ E h and i ∈ Σ h ; Φ f i = O h k+d-1 , ∀f ∈ F h and i ∈ Σ h (5)
and if the approximation f u h is accurate with the order k + 1, then the truncation error satisfies the following relation

|E u h , ϕ | ≤ C(ϕ, f , u)h k+1 ,
with C a constant which depends only on ϕ, f , and u. It can be shown, under the previous hypothesis, that Φ e = O h k+d and if there exists a constant β E i , such that Φ e i = β e i Φ e , then the condition ( 5) is satisfied provided that β e i is uniformly bounded. Such a condition is historically called linearity preserving. To determine the conditions that must be satisfied by the numerical scheme in order to have nonoscillatory solutions, the distributed residual on a generic element e is re-written in the following form

Φ e i = j∈Σ e h j =i c e ij (u i -u j ),
with the coefficients c e ij that in general depend on the solution. The same would hold for the face residuals. It is well known that if the coefficients c e ij are positive numbers, then a maximum principle holds under a CFL-like condition. The systematic construction of formaly high order RD schemes is recalled later in the text.

III. Extension to the diffusion terms

When in the governing equation (1a) diffusive phenomena are considered together with the advective terms, the following advection-diffusion equation is obtained

∇•f (u) = ∇• K∇u on Ω ⊂ R d , d = 2, 3 (6)
where the tensor K, generally function of u and ∇u, is such that the system is dissipative. As well know, the relative importance of the advection and the diffusion is described by a non-dimensional parameter, the Peclet number, Pe = a h/ρ(K), with ρ(K) the spectral radius of K. In the advection and diffusion limits Pe → ∞, Pe → 0, respectively, while Pe ∼ 1 when advection and diffusion are equally important.

To extend RD methods to advection-diffusion problems, different strategies have been considered to compute and to distribute the residual associated with the diffusion terms. On a first attempt, based on the physical intuition that the diffusion has an isotropic behavior in the space, RD schemes for the advection terms were coupled with the Galerkin discretization of the diffusion terms, [START_REF] Van Der Weide | A parallel, implicit, multi-dimensional upwind, residual distribution method for the Navier-Stokes equations on unstructured grids[END_REF][START_REF] Paillere | Multidimensional upwind residual distribution schemes for the Euler and Navier-Stokes equations on unstructured grids[END_REF] but a truncation error analysis revealed that this simple approach results in a first order accurate scheme when advection and diffusion have the same order of magnitude. [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF] A different approach, which was developed for two-dimensional schemes on triangular grids, considered a hybridization of the RD method with a Petrov-Galerkin scheme by the means of a scaling parameter, function of the Peclet number. [START_REF] Ricchiuto | On uniformly high-order accurate residual distribution schemes for advection-diffusion[END_REF] A key aspect is that a RD scheme with an uniform order of accuracy in all the range of the Peclet numbers should not consider two different distribution schemes for the advection and diffusion terms, but only one distribution process has to be performed for the residual of the whole equation, namely

Φ e = Ωe ∇•f u h -∇• K∇u h dΩ.
To put the previous expression in term of a boundary integral, one has to cope with the fact that the normal component of the gradient of the numerical solution, ∇u h •n, is in general discontinuous on the faces of the elements. Suppose, now, that an unique value of the gradient is available at each DOF, the gradients can be interpolated with the same shape functions used for the solution and the total residual on the element can be written as follows

Φ e = ∂Ωe f u h -K ∇u h • n d∂Ω, (7)
where ∇u h is the interpolated gradient of the numerical solution, which is now continuous on the faces of the elements. Note that K is defined as K(u h , ∇u h ).

Once the total residual is evaluated, it can be distributed to the DOFs of the elements by the distribution coefficients, β e i . This strategy has been adopted in [START_REF] Paillère | Multidimensional upwind residual distribution schemes for the convection-diffusion equation[END_REF] to construct a second order RD scheme for advectiondiffusion problems on triangular grids an has been extended to the third order in. [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF] In, [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF] this has been extended to third order also with in addition a strict control of numerical oscillations. In these works the distribution process is done with purely advective distribution coefficients, which is not appropriate in the diffusion limit. A more general scheme consists in using distribution coefficients which are function of the local Peclet number in order to recover an isotropic scheme in the diffusion limit and an upwind scheme in the advection limit. [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF][START_REF] Chou | High order residual distribution conservative finite difference WENO schemes for convection-diffusion steady state problems on non-smooth meshes[END_REF] Another attempt in that direction is given by, 5 the scheme give satisfactory results except in the region P e ≈ 1, which is typical of a boundary layer. Hence the present contribution can be viewed as an improvement over the previous references.

The key idea of the Eq. ( 7) is the reconstruction of the gradient of the numerical solution at each DOF of the grid and is one of the issue analyzed in this work. Indeed, numerical experiments show that in order to obtain a high order accurate solutions, the gradients must be recovered with the same order of accuracy of the solution.

An alternative approach has been proposed by Nishikawa for diffusion problems [START_REF] Nishikawa | A first-order system approach for diffusion equation. I: Second-order residual-distribution scheme[END_REF] and advection-diffusion problems, [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF] it consists in reinterpreting the advection-diffusion scalar equation as an equivalent hyperbolic first order system, in this way the gradient recovery is no longer necessary, but the price to pay is the increment of the unknowns of the problem due to the fact that a system of equations must be solved instead of a single scalar equation.

IV. Residual distribution discretization of advection-diffusion problems

In the previous sections the distribution process of the total residual is expressed through the use of generic distribution coefficients, in this section is described how actually to perform this step, for a scalar case.

In the past years different RD schemes were developed with the objective to construct upwind schemes for linear, quadratic and cubic triangles/quadrangles and tetrahedra/hexahedra, see e.g. [START_REF] Abgrall | Toward the ultimate conservative scheme: Following the quest[END_REF][START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF][START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF] In this work, the attention is focused on the construction of central schemes which can be formulated on every type of element and which can be easily extended to high order approximations. Linear and non-linear schemes are considered.

A. Central linear and non-linear RD schemes

Linear scheme

The linear scheme proposed in this work is the extension to the integral formulation of the classical Ni's Lax-Wendroff scheme, [START_REF] Ni | A multiple grid scheme for solving the Euler equations[END_REF] namely

Φ e i = Φ e N e dof + Ωe a•∇ψ i τ a • ∇u h -∇• K∇u h dΩ, (8)
where N e dof is the number of DOFs on the element e, the scaling parameter τ is defined as follows

τ = 1 2 |Ω e | j∈Σ e h max(k j , 0) , with k j = 1 2 ā•n j and n j = Ωe ∇ψ j dΩ,
where ā represents the arithmetic average of the advection velocity on the element. The scheme ( 8) is linearity preserving but not positive, and due to the integral formulation, it is valid for any type of element and for any order of approximation. The scheme is conservative since i∈Σ e h Φ e i = Φ e , due to the fact that i∈Σ e h ∇ψ i = 0, and it is consistent; when the exact solution is injected in the Eq. ( 8) the residual is zero because the total residual vanishes, by definition, and the integral term vanishes as well due to the fact that the term between the brackets is exactly the governing equation.

Non-linear scheme

Non-linear schemes are needed to combine the non-oscillatory behavior of the numerical solution with the high order discretization. The basic idea to construct a non-linear scheme is to start with a first order, positive scheme, and to map its distributed residuals onto a set of positive and non-linear residuals. To see in practice how to construct of a non-linear scheme, consider the first order accurate and positive Rusanov's scheme (also know as Lax-Friedrichs scheme) c defined as

Φ e i = Φ e N e dof + 1 N e dof α j∈Σ e h j =i u i -u j ), ∀i ∈ Σ e h , (9)
with α ≥ max j∈Σ e h |k j | + ν > 0. Since the Rusanov's scheme is first order accurate, its distribution coefficients, β e i = Φ e i /Φ e , 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 βe i . A common choice for the map is the following

26 βe i = max β e i , 0 j∈Σ e h max β e j , 0 . ( 10
)
The use of a central 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] The cure to this problem consists in adding a filtering term by means of a streamline dissipation term

Φe i = βe i Φ e + θ e h (u h ) Ωe a•∇ψ i -∇• ν∇ψ i a • ∇u h -∇• ν∇u h dΩ. ( 11
)
The role of the parameter θ e h (u h ) is double. It provides the correct scaling of the streamline filtering and it makes sure that the filtering term is added only in the smooth regions of the solution. The following definition is used here

θ e h (u h ) = ε(u h )      j∈Σ e h |ā•n j | + ν 2|Ω e |      -1
, with ε(u h ) a smoothness sensor such that ε(u h ) ∼ 1 in smooth regions and ε(u h ) ∼ 0 around discontinuities.

B. Improved discretization of the diffusion terms

Numerical experiments reveal that the schemes ( 8) and ( 11) applied to the discretization of the advectiondiffusion problem are unsatisfactory from the point of view of the accuracy and the robustness. In order to obtain a better discretization of the diffusive terms, the advection-diffusion equation ( 6) is written in the form of a first order system as follows

   ∇•f (u) -∇•(νq) = 0 q -∇u = 0 (12)
Consider now a numerical scheme for the previous system obtained by writing the weak form of the system plus a streamline stabilization term:

Ωe ψ i ∇•f (u h ) -∇•(νq) q -∇u h dΩ + Ωe A•∇ψ i τ ∇•f (u h ) -∇•(νq) q -∇u h dΩ = 0, (13)
where

A•∇ψ i =        a•∇ψ i -ν ∂ψ i ∂x -ν ∂ψ i ∂y - ∂ψ i ∂x 0 0 - ∂ψ i ∂y 0 0        ,
and the term τ is assumed to be of the following form

τ =    τ a 0 0 0 τ d 0 0 0 τ d    ,
with τ a and τ d are strictly positive coefficients. Supposing, now, that the gradient of the numerical solution has been recovered at each DOF, one can replace the second equation of the system (12) with the approximation ∇u u ≃ ∇u h and consider only the first equation, which now reads

Ωe ψ i ∇•f (u h ) -∇•(ν ∇u h ) dΩ + Ωe a•∇ψ i τ c a•∇u h -∇• ν ∇u h dΩ + Ωe ν ∇ψ i • τ d ∇u h -∇u h dΩ = 0. ( 14
)
The first two integrals of the previous equation represent a discretization for the scalar advection-diffusion equation by the means of a central scheme plus a streamline stabilization term, in the same way as shown in equation Eq. (8). The last integral represents an additional stabilization term, for the diffusive part only, which vanishes in the advective limit and the parameter τ d is dimensionless. It is interesting to note that the additional term penalizes the difference between the discontinuous and the interpolated gradients, on each element. Note that with a slightly different procedure, a similar stabilization term for the diffusive part has been obtained by Nishikawa [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme I: Basic principle[END_REF] for the RD discretization of the diffusion problem.

With the Eq. ( 14) in mind, it is proposed here a modification of the schemes ( 8) and (11) in order to include the extra stabilization term for the diffusive part of the equation. In practice, the linear scheme reads

Φ e i = Φ e N e dof + Υ(Pe e ) Ωe a•∇ψ i τ a • ∇u h -∇• ν∇u h dΩ + 1 -Υ(Pe e ) Ωe ν ∇ψ i • ∇u h -∇u h dΩ, (15)
while the non-linear scheme becomes

Φe i = βe i Φ e + Υ(Pe e ) θ e h (u h ) Ωe a•∇ψ i -∇• ν∇ψ i a • ∇u h -∇• ν∇u h dΩ + 1 -Υ(Pe e ) Ωe ν ∇ψ i • ∇u h -∇u h dΩ, (16)
With the local Peclet number, defined as Pe e = a h e /ν, where h e is the characteristic length size of the element e. The function Υ(Pe e ) is defined such that Υ(Pe e ) → 0 in the diffusive limit and Υ(Pe e ) → 1 in the advective limit. In the numerical simulations the following definition is used

Υ(Pe e ) = max 0, 1 -1/ Pe e .
Note that in the schemes ( 15) or ( 16), the use of the blending function Υ(Pe e ), makes possible to recover, in the case of the pure advection, the same scheme used for the discretization of pure advective problems, while in the case of pure diffusion problems only the stabilization term for the diffusive terms is taken into account.

V. Gradient recovery strategy

As explained in Section III, 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 every DOF of the grid, then the nodal values of gradients are interpolated with Lagrangian functions on each element.

The key point is the recovery of the gradients at the DOF. Attention is focused on the possibility to obtain a high order gradient recovery, i.e., the gradient is recovered with the same order of accuracy of the solution. For simplicity, the description is always limited to the two-dimensional case, but the extension to the three spatial dimensions is straightforward.

In, [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] different way of reconstructing the gradient have been studied and compared in great details, namely the area-weighted method, the L 2 -projection method, the least-square approximation and the Super Patch Recovery technique proposed 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] In, [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] it is shown that only the last technique is able to guaranty high order approximation of gradient without enlarging the stencil, for this reason only this approach is described here.

In this paper, only on the approximation properties are considered, a more detailed analysis, including in particular a convergence study of the convection diffusion problem for high Peclet number, with linear and non linear problems, is described in. [START_REF] Abgrall | High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF] Surprisingly, the same method also works for pure ellitic problems.

A. Super-convergent patch recovery

In the field of the Finite Element method applied to the mechanical structures, it is know that the stresses (gradients of the displacements) sampled at certain points in a element possess a super-convergent property. This means that the stresses have the same order of accuracy of the displacements. [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] obviously by tensor product such points can be defined also in the case of quadrangles and hexahedra. For triangles or tetrahedra such property cannot be rigorously shown, at least up to our knowledge, but numerical experiments confirm that the stresses sampled at certain points have high order of accuracy.

Accepting the fact that gradients are sampled with high order accuracy in certain points of the element, it is possible to compute gradients which are high order accurate within all the element. Indeed, if at sampling points the value of gradients is accurate to order k + 1, by using a polynomial of degree k (the same order used to interpolate the solution), it is possible to obtain an approximation which has high order accuracy everywhere within the elements, if this polynomial is made to fit the values of the sampled gradients in a least square manner. 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] Assume that the numerical solution u h of the problem is known at each DOF of the grid to the (k+1)-th order of accuracy. The aim is to obtain the values of the solution gradient, ∇u h at all the DOFs with same order of accuracy of the solution. The components of the recovered gradient, at the generic DOF i, are written in a polynomial form as follows ∂u h ∂x i = p T a x and ∂u h ∂y i = p T a y , with p T (x) = (1, x, y, x 2 , . . . , x k+1 , x k y, . . . , y k+1 ), a x = (a x1 , a x2 , . . . , a xm ) and a y = (a y1 , a y2 , . . . , a ym ). Assuming that N s sampling points, ξ ℓ , ℓ = 1 . . . N i s , are available for each DOF i, the objective is to minimize the following functions

F x = N i s k=1 ∂u h ∂x (ξ k ) -p T k a x 2 and F y = N i s k=1 ∂u h ∂y (ξ k ) -p T k a y 2 ,
with p k = p(ξ k ). What the sampling points are, will be pointed out later.

To compute the coefficients a x and a y , a small least square problem, with matrix A i , must be solved for each DOF i of the grid. The dimension of the matrix A i are determined by the number of sampling points N i s and by the degree of the polynomials used to express the recovered gradient, that is

A i ∈ R N i s ×m
, where m is the number of the coefficients in the vector a x or a y . The problem admits a unique solution if Rank A i = m, which is always satisfied in the case in which N i s ≥ m. It is worth also noticing that since the matrix A i depends only on the geometry, for a given grid, the matrix coefficients needs to be computed only once at beginning of the simulation.

Generally, the number of elements which share the same node within the domain is such that the condition N s ≥ m is always satisfied, this means that the gradient recovery is compact because it involves only the elements contained within the support of a grid node. For the nodes belonging to the boundary of the grid, the condition N s ≥ m might not be satisfied without enlarging the stencil, otherwise the problem is ill conditioned. In this case, to avoid the use of larger stencil for a boundary node, it is possible to obtain the value of the recovered gradient with the same polynomial expansion used for nearest domain node.

The structure of the sampling points for gradient reconstruction technique, is now discussed. For each degree of freedom i, the set S i of sampling points is defined by

S i = e∈E h,i S e
i , where the sampling point for the element e depend on the structure of the element. Here, details about the sampling points are reported only for quadrilateral and triangular elements (the extension to element in three spatial dimensions is rather obvious). In Fig. 1 are shown examples of patches used to recover the gradient for a domain node in the case of quadrangular and triangular elements. For a quadrangle the sampling points are defined uniquely. Considering a reference segment defined as x = [-1, 1], the sampling point is the point x = 0 in the case of a linear element, while in the case of a quadratic element the sampling points have coordinates ±1/ √ 3. The sampling 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 sampling point is the point with barycentric coordinated λ = (1/3, 1/3, 1/3), while in the case of a quadratic triangle the sampling points are not uniquely defined. In Fig. 2 are shown three examples of sampling points used. In the first option, see Fig. 2-(a), are used three points with barycentric coordinates

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

λ 3 = 1 3 , 1 3 , 2 3 ,
In the second option (Figure 2-(b)) four points are used with barycentric coordinates

λ 1 = 1 3 , 1 3 , 1 3
, λ 2 = (0.6, 0.2, 0.2) , λ 3 = (0.2, 0.6, 0.2) , λ 4 = (0.2, 0.2, 0.6) . Another option (Fig. 2-(c)) consists in taking as sampling points the three points with barycentric coordinates

λ 1 = 1 2 , 1 2 , 0 , λ 2 = 0, 1 2 , 1 2 , λ 3 = 1 2 , 1 2 , 0 . (a) (b) (c)

B. Results and discussion

To study the accuracy of the presented gradient recovery strategies the following function is used

u = -cos(2πη) exp ξ 1 - √ 1 + 16π 2 ν 2 2ν , (17)
with η = a y x -a x y and ξ = a x x + a y y. Here a x = 0.5, a y = √ 3/2 and ν = 0.01. The solution is infinitely differentiable with continuous gradients. The computations of the recovered gradients are performed on four different kind of grids, shown in Fig. 3, namely unstructured grids of triangles, quadrangles and hybrid elements and highly deformed unstructured meshes of triangles, obtained randomly perturbing a regular grid.

The error of the recovery procedure is computed as the L 2 norm of the difference between the computed gradient, ∇u, and the exact gradient, ∇u ex normalised by ||∇u ex || L 2 . These errors and the orders of convergence of different gradient recovery methods on triangular grids are shown in Table . 1. In the case of linear elements, the differences between the recovery methods is small and all the schemes reach almost the second order accuracy. The L 2 -Projection and the SPR-ZZ methods have the smallest level of error, but the former scheme is much more expensive. In the case of quadratic elements, the methods have an order of accuracy no more than two except for the SPR-ZZ method which show almost third order accuracy. It is also worth noticing that the errors obtained with this method are one order of magnitude smaller than those obtained with other methods. In Table . 2 are reported the errors obtained with the SPR-ZZ procedure on quadratic triangular elements for the three different sampling strategies shown in Fig. 2. The first strategy clearly guaranties the smallest level of error, while the four-points strategy is unsatisfactory. Note that the results reported in the previous paragraph are obtained with the first strategy.

The errors of the recovery methods on unstructured grids of quadrangles and of hybrid elements are reported in Table . 3 andTable. 4, respectively. The behavior of the recovery methods is the same as observed in the case of triangular grids. Table . 5 shows the errors computed on a sequence of highly distorted triangular grids, the performance of the recovery methods is not optimal anymore due to very poor quality of the meshes, nevertheless the errors obtained with the SPR-ZZ methods are always much smaller than those obtained with other procedures.

Linear elements

Quadratic elements

N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2

Linear elements

Quadratic elements Table 3. Accuracy study of the SPR-ZZ method gradient recovery methods with linear and quadratic elements on unstructured grids of quadrangles.

N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2
Linear elements Quadratic elements Linear elements Quadratic elements Table 5. Accuracy study of the SPR-ZZ method with linear and quadratic elements on grids of randomly distorted triangles.

N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2
N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 100 1.2067E-01 --

VI. Extension to systems of equations

Consider now the RANS equations with the Spalart-Allmaras turbulent model, the governing equations read

∇•f a (u) -∇•f v (u, ∇u) = S(u, ∇u),
where u is vector of the conservative variables, f a and f v are the advective and diffusive flux functions, respectively, and S is the source term

u =          ρ m E t µ ⋆ t          , f a =          m m ⊗ m ρ + P I E t + P m ρ m ρ µ ⋆ t          , f v =           0 τ t τ t • m ρ -q t η∇ µ ⋆ t ρ           , S =          0 0 0 S SA         
, where ρ is the density, m is the momentum vector, E t is the total energy for unit volume, µ ⋆ t is the turbulent working variable, P is the pressure, I is the identity matrix, S SA is the source term of the Spalart-Allmaras equation

S SA = c b1 ωµ ⋆ t -ρc w1 f w µ ⋆ t ρ d min 2 + 1 σ SA ρc b2 ∇ µ ⋆ t ρ 2 .
In the production term the modified vorticity function is expressed as

ω =      |ω| + ω, ω > -c v2 |ω| |ω| + |ω| c 2 v2 |ω| + c v3 ω (c v3 -2c v2 )|ω| - ω , ω ≤ -c v2 |ω|,
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 v1 = Ψ 3 Ψ 3 + c 3 v1 , with Ψ =    0.05 log 1 + e 20χ , χ ≤ 10, χ, χ > 10, and χ
= µ ⋆ t µ . 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 and η = µ 1 + Ψ σ SA
where d min is the distance to the nearest wall, c b1 = 0.1355,

σ SA = 2/3, c b2 = 0.622, k = 0.41, c w1 = c b1 k 2 + (1 + c b2 ) σ SA , c w2 = 0.3, c w3 = 2, c v1 = 7.
1, c v2 = 0.7, and c v3 = 0.9.

The modified viscous stress tensor and the heat flux vector read, respectively

τ t = - 2 3 (µ + µ t )(∇•v)I + (µ + µ t )(∇ T v + ∇v), q t = c p µ P r + µ t P r t ∇T,
where v is velocity vector, c p is the specific heat at constant pressure, T is the temperature, P r and P r t are the Prandtl number and the turbulent Prandtl number, respectively, µ is the fluid viscosity and µ t is the turbulent viscosity µ t = µ ⋆ t f v1 .

The equation for laminar flows are obtained by setting µ t = 0 and by disregarding the last equation in the governing equations.

Note that the sub set of equation without the turbulent viscosity equation can be rewritten as The system of the Navier-Stokes equations can be easily rephrased as follows

∇•f a (u) -∇•(K∇u) = 0 with K∇u = K i,j ∂u ∂x j .

A. RD formulation for systems of equations

The numerical schemes introduced for the scalar advection-diffusion problem can be easily extended to case of a RANS or Navier-Stokes equation.

The gradient of the conservative variables is first reconstructed at each DOF, as shown in section V so the total residual on the generic element e can be computed as follows

Φ e = ∂Ωe f a (u h ) -f v (u h , ∇u h ) d∂Ω - Ωe S(u h ) dΩ.
To distribute the total residual a simple generalization of the Lax-Wendroff scheme can be used

Φ e i = Φ e N e dof + Υ(Re e ) Ωe A•∇ψ i τ A • ∇u h -∇• K ∇u h -S(u h ) dΩ + 1 -Υ(Re e ) Ωe K∇ψ i • ∇u h -∇u h dΩ with τ = 1 2 |Ω e |   j∈Σ e h K + j   -1 , with K j = 1 2 Ā•n j and n j = Ωe ∇ψ j dΩ,
and Re e the local Reynolds number and where Ā = A ūh is the Jacobian of the advective flux evaluated at the arithmetic mean state ūh , and K = K(u h ) is the Jacobian of the viscous flux. On the other hand a non-linear scheme can be constructed, based on the procedure illustrated in section IV,

Φ e i = Φ e i + Υ(Re e ) Ωe A•∇ψ i -∇• K∇ψ i τ A • ∇u h -∇• K ∇u h -S(u h ) dΩ + 1 -Υ(Re e ) Ωe K∇ψ i • ∇u h -∇u h dΩ.
To obtain the high order distributed residual Φ e i , the low order residuals computed by the means of the Rusanov's scheme are first projected on the space of the left eigenvectors, then the same procedure shown for the scalar equation is directly applied and finally the limited residual are project back to the physical space, see, 2 for more details.

B. Implementation details

From a practical point of view the non linear system of equations ( 2) is solved via a backward Euler method with a pseudo-transient continuation technique

u n+1 -u n ∆t n = -R(u n+1 ), (18)
where the right-hand side of the previous equation takes into account both the domain and the boundary residuals, and n = 0, 1, . . . , ∞. As standard practice the Eq. ( 18) is linearized at each time step n and the following linear system is obtained

I ∆t n + ∂R(u n ) ∂u ∆u n = -R(u n ), (19)
with ∆u n = u n+1u n . Different technique has been proposed to solve the previous linear system, one the most popular is the GMRES algorithm. An interesting alternative is the matrix-free approach, [START_REF] Knoll | Jacobian-free newton method: a survey of approaches and application[END_REF] which has been found to be very effective for laminar flows simulations. For turbulent simulations, however, this approach is unsatisfactory; usually the cannot reach the steady state solutions or diverge in the worst cases. An effective and robust technique for the solution of the problem (19), is the non-linear LU-SGS method, 27 which is briefly described in the following. Instead of using a Krylov-based solver, the Eq. ( 19) is solved with symmetric Gauss-Seidel method as follows

I ∆t n + ∂R i ∂u i ∆u (k+1) i = -R i (u n ) - j∈Σ h,i j =i ∂R i ∂u j ∆u ( * ) j , ∀i ∈ Σ h , (20)
where Σ h,i is the set of the DOFs which belong to the stencil of the DOF i, the superscript * indicates the most recently updated solution, k is an inner iteration index, 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 ( * ) , is possible to linearize R i u ( * ) as

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

R i u n + j∈Ωi j =i ∂R i ∂u j ∆u ( * ) j = R i u ( * ) - ∂R i ∂u i ∆u ( * ) i .
Substituting the previous relation in (20) 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 . (21)
The ( 21) is solved with the forward and backward sweeps and at the beginning of each step the small diagonal block of the Jacobian matrices in the left hand side can be inverted using LU decomposition. Note that the right hand side of the previous equation is nothing but the residual evaluated with the latest available solutions. The last term in the right hand side of the ( 21) is usually dropped for steady simulations, in order to accelerate the convergence rate. When the symmetric Gauss-Seidel converges one is actually solving the original equation ( 18) instead of the linearized version (19), for this reason the algorithm is called non-linear LU-SGS. It is worth noticing that only the diagonal part of Jacobian matrix is actually required for this implicit method.

For the simulation of laminar flows, the matrix-free approach has been used, while the non linear LU-SGS method has been used for turbulent problems. In all the simulations, the solution is considered to be at the steady state when the L 2 norm of the density residual is reduced by ten orders of magnitude respect to the initial value. The laminar flow over the NACA-0012 airfoil is computed with the proposed RD scheme. The Reynolds number based on the airfoil cord is 5 000, the Mach number is 0.5 and the angle of incidence is zero. Unstructured grids of triangles are consider, second and third order computation are performed. In Fig. 6 are reported the solutions computed with the linear scheme with P1 and P2 elements. The solution with the P1 elements has been performed on a grid obtained from that with P2 elements and splitting each P2 triangle with four P1 triangles, in such a way the number of DOF for the second and third order simulation is exactly the same. Note, in Fig. 6, that although there is not much difference in the Mach number contours between the second and the third order simulations, the streamlines near the trailing edge are very different and only the third order scheme is able to reproduce the symmetric recirculation bubble which is the main feature of this test case.

In Fig. 7 are reported the values of the lift coefficient for the second and third order linear scheme, with two different types of gradient reconstruction. Since the airfoil and flow conditions are symmetric, the absolute value of the lift coefficient can be taken as error estimation. Note, how it was expected, that non only the third order scheme gives a smaller error respect to the second order scheme for a given number of DOF, but the accuracy of the solution is also strictly related to the accuracy of the gradient reconstruction method.

C. Laminar flow over a delta wing

Consider now the laminar flow over a delta wing. The Reynolds number based on the root cord is 4 000, the Mach number is 0.3 and the angle of incidence is 12.5 • . The flow filed is caractherized by a big detached vortex over the wing which is convected far away behind the wing, and a small vortex which originates near the top edge of the wing. In Fig. 8 are depicted a coarse mesh used for the simulation and an example of third order solution computed with the linear scheme.

In Fig. 9 are reported the drag and the lift coefficients, for the third and second order linear scheme, on a sequence of unstructured grids of tetrahedra. The SPR-ZZ method has been used to reconstruct the gradients. The results are compared against the reference value obtained with a DG solver on an adapted grid. [START_REF] Leicht | Error estimation and anisotropic mesh refinement for 3d laminar aerodynamic flow simulations[END_REF] The results converge to the reference value refining the grid, the advantage of the third order scheme is evident for the lift coefficient. While for the drag coefficient the gain of using third order scheme respect to the second order one is reduced, probably due to the poor quality of the mesh.

D. Turbulent flow over a RAE airfoil

A turbulent flow over a RAE airfoil is considered. The Reynolds number, based on the airfoil cord, is 6.5 × 10 6 , the Mach number is 0.734 and the angle of attack is 2.79 • . The Spalart-Allmaras equation has been used to model the turbulent effects. An example of mesh used is depicted in Fig. 10 together with a third order solution obtained with the non-linear scheme and the SPR-ZZ method. In Fig. 12 the pressure and the friction coefficient along the airfoil of the second and third order solution are compared against the experimental values. The agreement is very good, the shock is sharply captured and in the third order solution the small boundary layer separation behind the shock is also visible.

In Fig. 11 are reported the convergence histories using the matrix-free and the non-linear LU-SGS methods for a third order simulations in the medium grid, it is evident how the LU-SGS is very effective to make the scheme converge to the steady state solution, while the matrix-free approach, after an initial reduction on the non-linear residual, is not able to converge. This behavior is probably due to the fact the matrixfree approach still require the computation of an approximated Jacobian for the preconditioning, here the LU-SGS has been used as preconditioner. The discrepancy in the quality of actual Jacobian (computed numerically) and the approximated one used for the preconditioning prevents the solver to converge to the steady state. In the non-linear LU-SGS method only the diagonal block of the approximated Jacobian is required, making the solver less sensitive to the accuracy with which the Jacobian is computed.

In Fig. 13 a convergence study on three grids is performed for the second and third order non-linear scheme. The third order solution on the finest grid is takes as reference solution. Note how the reference 20 value is reached very soon by the third order scheme, while the second order scheme require much more degrees of freedom to get the same value.

E. Turbulent ONERA M6 wing

As last case, a turbulent flow over the ONERA M6 wing is considered. The inflow Mach number is M = 0.8395, the Angle of Attack is α = 3.06 • . The Reynolds number is set to Re = 11.72 × 10 6 , based on the mean aerodynamic chord. The grid, provided by Dassault and shown in Fig. 14, consists of 77 061 points and 443 458 tetrahedra, although quite coarse, it is enough for the third order solver. The non-linear scheme has been used in combination with the SPR-ZZ method for the gradient reconstruction and the non-linear LU-SGS method for the solution on the non-linear system. The contours of the pressure coefficient are shown in Fig. 15 and in Fig. 16 are reported the value of the pressure coefficients at different station along the wingspan compared with the experimental data.

VIII. Conclusions

In this work some advances in the RD method have been presented and numerical experiments showed that such kind of schemes can be used also for the simulation of complex turbulent flows. Several features have been introduced respect to the classical RD method traditionally used; the new treatment of the viscous term handle advection and diffusion contributions within the same method. Furthermore, the SPR-ZZ technique for the reconstruction of the gradient of the numerical solution proved to be much more accurate respect to other standard approaches, without increasing the computational cost. Finally, the use of the non-linear LU-SGS method, for the solution of the non linear system of equations, has proved to be more effective and robust respect to Krylov based methods for the RD schemes.

Figure 1 .

 1 Figure 1. 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 recovered and the symbols (△) indicate the super-convergent sampling points.

Figure 2 .

 2 Figure 2. Three different examples of sampling points for quadratic triangles.

Figure 3 .

 3 Figure 3. Example of different kinds of grid used to test the accuracy of the gradient recovery procedures.

Strategy Fig. 2 -

 2

Figure 5 .

 5 Figure 5. Observed order of accuracy for the manufactured solution problem with different Reynolds number. Top: linear scheme bottom: non-linear scheme. Left: order of accuracy of the solution, right: order of accuracy of the gradient

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. Mach Number contours (top) and streamlines neat the trailing edge (bottom) for the third (left) and second (right) order linear scheme.

Figure 9 .Figure 10 .Figure 11 .

 91011 Figure 9. Drag (left) and lift (right) coefficients, for the delta wing test case, obtained with the second and third order schemes.

Figure 12 .

 12 Figure 12. Pressure (left) and friction (right) coefficient along the airfoil for the second and third order solution. Experimental results are also reported

Figure 13 .Figure 14 .

 1314 Figure 13. Drag (left) and lift (right) coefficients for the second and third order solution for successively refined grids.

Figure 15 .Figure 16 .

 1516 Figure 15. Contours of the pressure coefficient of the third order solution.

Table 1 .

 1 Accuracy study of different gradient recovery methods with linear and quadratic elements on an unstructured grids of triangles.

	∂u ∂y	O L 2

Table 2 .

 2 Accuracy study of SPR-ZZ recovery methods with quadratic triangular elements for different sampling strategies.

				3.04	7.4 10 -4	2.96	7.6 10 -4	3.0	1.5 10 -3	2.32	1.5 10 -3	2.4
	6673 9.55 10 -5	2.9	1.0 10 -4	2.82	1.1 10 -4	2.72	1.2 10 -4	2.7	3.7 10 -4	2.09	3.7 10 -4	2.1
	26401 1.36 10 -5	2.8	1.6 10 -5	2.62	2.1 10 -5	2.47	2.2 10 -5	2.4	9.2 10 -5	2.03	9.2 10 -5	2.0
	105025 2.23 10 -6	2.6	3.0 10 -6	2.44	4.4 10 -6	2.26	4.7 10 -6	2.2	2.3 10 -5	2.01	2.3 10 -5	2.0

Table 4 .

 4 Accuracy study of the SPR-ZZ method with linear and quadratic elements on unstructured grid of hybrid elements.

	∂u ∂y	O L 2

  Downloaded by Biblioteca Dipartimento Ingegneria on January 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/6.2013-2562

	M	M
	0.55	0.55
	0.5	0.5
	0.45	0.45
	0.4	0.4
	0.35	0.35
	0.3	0.3
	0.25	0.25
	0.2	0.2
	0.15	0.15
	0.1	0.1
	0.05	0.05

Downloaded by Biblioteca Dipartimento Ingegneria on January 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/6.2013-2562

Downloaded by Biblioteca Dipartimento Ingegneria on January 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/6.2013-2562

c Other low order, non oscillatory schemes can considered, like for example a Finite Volume scheme written as a RD scheme, see e.g. 16 Downloaded by Biblioteca Dipartimento Ingegneria on January 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/6.2013-2562

Downloaded by Biblioteca Dipartimento Ingegneria on January 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/6.2013-2562

Downloaded by Biblioteca Dipartimento Ingegneria on January 23, 2014 | http://arc.aiaa.org | DOI: 10.2514/6.2013-2562

Acknowledgements.

R.A. is partialy funded by the ERC Advanced Grant "ADDECCO" (contract #226316). D. de Santis has been financed by the FP7 strep "IDIHOM" (contract #265780).

VII. Numerical experiments

A. Manufactured Solution for the laminar NS equations

To assess the order of accuracy of the RD scheme, the method of the manufactured solutions is used. The numerical domain is a square of length L and the form of the primitive variables is chosen as a function of sines and cosines The conservative variables are first computed from the previous definition of the primitive variables and then are made dimensionless using as reference quantities the values at the point of coordinates (0, 0). In figure 4 are reported the contours of the conservative variables used as manufactured solution.

To assess the order of accuracy of the numerical scheme, a sequence of four unstructured grids of triangles is considered. On each grid the L 2 norm of the error of the solution and of the gradient is computed. Numerical experiments are performed with different Reynolds numbers, in order to verify that the formal order of accuracy is preserved in the whole advection-diffusion spectrum. The observed order of accuracy is computed as the mean slope of the curve which represents the error as function of the parameter h = 1/ √ N DOF . In Fig. 5 are reported the observed orders of accuracy of the solution and of the gradients for the second and third order formally accurate schemes. Two type of gradient reconstructions are used, namely the Green-Gauss and the SPR-ZZ methods. Results for both the linear and non linear scheme are reported.

Consider first the linear scheme, at the second order there is not much difference between the two type of gradient reconstructions in the observed accuracy order of the solution: both methods guarantee formal second order accuracy for the solution, however only the SPR-ZZ method gives also second order accuracy for the gradient. At the third order, situation changes drastically, only the SPR-ZZ method gives third order accuracy for all the range of the Reynolds number, for both solution and gradient. With the use of the Green-Gauss method, the formal accuracy of the solution is lost in the diffusion dominate regime. It is worth noting also that the Green-Gauss method gives at most second order accuracy for the gradient, to be compared with the third order accuracy guaranteed by the SPR-ZZ method.

The same considerations done above for the linear scheme are also valid for the non-linear scheme, except for the fact that the loose of accuracy, for the non-linear scheme, is more severe if gradients are not reconstructed properly.