domain: math phys info

text:
Victorien Menier

Adrien Loseille

Frédéric Alauzet

CFD Validation and Adaptivity for Viscous Flow Simulations

This article deals with the description and the validation of the unstructured viscous RANS flow solver Wolf and its compliance with anisotropic unstructured mesh adaptation. We first describe the numerical scheme of the flow solver and its validation on standard test cases. We then give a general description of the mesh generator and its main features for viscous mesh generation and adaptation. From these points, we define some simple strategies to perform adaptive computation for RANS solutions. Each technique is illustrated with a numerical example and we discuss its robustness, effectiveness and ability to handle complex geometries.

Introduction

Mesh adaptation provides a way to control the accuracy of the numerical solution by modifying the domain discretization according to size and directional constraints. When dealing with real life flow problems, Hessian-based unstructured mesh adaptation has already proved its efficiency to improve the ratio between the solution accuracy and the number of degrees of freedom (the problem complexity). [START_REF] Castro-Díaz | Anisotropic Unstructured Mesh Adaptation for Flow Simulations[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Gruau | 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF][START_REF] Li | 3D anisotropic mesh adaptation by mesh modification[END_REF][START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF][START_REF] Tam | Anisotropic mesh adaptation for 3D flows on structured and unstructured grids[END_REF] In addition, as a large number of physical phenomena are anisotropic by nature, anisotropic mesh adaptation improves even more this ratio.

In the context of high Reynolds compressible viscous flow simulations though, mesh adaptation remains challenging. It is notably due to steep gradients in the near-wall regions (boundary layers). Dramatic variations in the normal direction of variables such as the velocity require specific meshing techniques. The generation of quasi-structured meshes in the near wall regions has proved to be a reliable approach but their integration in the mesh adaptation loop is difficult. Indeed, to perform a mesh adaptation for a viscous flow simulation, one needs a reliable flow solver and suitable meshing technologies. The scope of this paper is to apply our flow solver and meshing technologies to adaptive simulations of viscous flows.

The paper is organized as follows. In Section I, we describe the Wolf flow solver. We focus on the CFD modeling, the turbulence model used, time step computation, numerical fluxes computation. Code verification and validation is based on comparison with other codes on standard test cases. In Section II, we describe the local remeshing strategies. In particular, we describe the algorithms involved for surface, volume, boundary layer and anisotropic cartesian mesh generation in an adaptive context. In Section III, we introduce four adaptive approaches for viscous simulations: (i) a fully unstructured, (ii) an hybrid with anisotropic element and a frozen boundary layer, (iii) an hybrid with anisotropic element and a regenerated boundary layer, and (iv) a fully adaptive cartesian approach.

I. Implementation and validation of the flow solver Wolf

This section details the in-house compressible flow solver Wolf to model viscous flows. We first describe a few key characteristics of its implementation: CFD modeling, turbulence model used, finite element/finite volume spatial discretization and the time implicit discretization. The intent of this section is not to describe the whole implementation in details but to give its main characteristics. Finally, we exhibit some results of the code verification and validation study that we led.

A. Modeling equations

The Reynolds Averaged Navier-Stokes (RANS) system relying to the Spalart-Allmaras model is composed of the compressible Navier-Stokes equations and the standard Spalart-Allmaras equation with no trip.

The compressible Navier-Stokes equations. The compressible Navier-Stokes equations for mass, momentum and energy conservation reads:

             ∂ρ ∂t + ∇ • (ρu) = 0 , ∂(ρu) ∂t + ∇ • (ρu ⊗ u) + ∇p = ∇ • (µT ) , ∂(ρe) ∂t + ∇ • ((ρe + p)u) = ∇ • (µT u) + ∇ • (λ∇T ) ,
where ρ denotes the density, u the velocity, e the total energy per mass, p the pressure, T the temperature, µ the laminar dynamic viscosity and λ the laminar conductivity. T is the laminar stress tensor:

T = (∇ ⊗ u + t ∇ ⊗ u) - 2 3 ∇. u I .
The variation of nondimensionalized laminar dynamic viscosity and conductivity coefficients µ and λ as a function of a dimensional temperature T is defined by the Sutherland law.

µ = µ ∞ T T ∞ 3 2 T ∞ + Su T + Su and λ = λ ∞ T T ∞ 3 2 T ∞ + Su T + Su ,
where Su = 110 is the Sutherland constant and the index ∞ denotes reference quantities. The relation linking µ and λ is expressed from the Prandtl laminar number: Pr = µc p λ with Pr = 0.72 for (dry) air .

This system can be rewritten under vectorial form:

W t + F 1 (W ) x + F 2 (W ) y + F 3 (W ) z = S 1 (W ) x + S 2 (W ) y + S 3 (W ) z
where W is the nondimensionalized conservative variables vector:

W = (ρ, ρu, ρv, ρw, ρE) T F(W ) = (F 1 (W ), F 2 (W ), F 3 (W )
) are the convective (Euler) flux functions:

F 1 (W ) = (ρu, ρu 2 + p, ρuv, ρuw, u(ρE + p)) T F 2 (W ) = (ρv, ρuv, ρv 2 + p, ρvw, v(ρE + p)) T F 3 (W ) = (ρw, ρuw, ρvw, ρw 2 + p, w(ρE + p)) T
and S(W ) = (S 1 (W ), S 2 (W ), S 3 (W )) the laminar viscous fluxes: where τ ij are the components of laminar stress tensor defined by:

S 1 (W ) = (0, τ xx ,
τ ij = µ ∂v i ∂x j + ∂v j ∂x i - 2 3 µ ∂v k ∂x k δ ij .
where (v i , v j , v k ) are the three components of the velocity and δ ij the Kroneker symbol.

Turbulence modeling. The chosen turbulence model is the one equation Spalart-Allmaras model:

7 ∂ ν ∂t + u • ∇ν = c b1 [1 -f t2 ] S ν -c w1 f w - c b1 κ 2 f t2 ν d 2 + 1 σ ∇ • ((ν + ν)∇ν) + c b2 ∇ν 2 + f t1 ∆u 2 .
where ν is the turbulent kinematic viscosity. In the standard model the trip term is being left out, i.e., f t1 = 0. Moreover, some implementations ignore also the f t2 term as it is argued that if the trip is not included, then f t2 is not necessary. [START_REF] Eca | Verification of RANS Solvers with Manufactured Solutions[END_REF] In Wolf, this simplified version has been considered and we prefered to write it under the following form which is more appropriate for its discretization with the finite element/finite volume method:

∂ρν ∂t + u • ∇ρν convection = c b1 Sρν production -c w1 f w ρ ν d 2 destruction + ρ σ ∇ • ((ν + ν)∇ν) dissipation + c b2 ρ σ ∇ν 2 dif f usion
Notice that this is not a conservative model. If conservative form of the Spalart-Allmaras is foreseen, we have to consider the variation of Catris and Aupoix. [START_REF] Catris | Density corrections for turbulence models[END_REF] The turbulent eddy viscosity is computed from:

µ t = ρνf v1
where

f v1 = χ 3 χ 3 + c 3 v1 and χ = ν ν with ν = µ ρ .
Additional definitions are given by the following equations:

f v2 = 1 - χ 1 + χf v1 and S = Ω + ν κ 2 d 2 f v2 where Ω = ∇ × u .
c w1 = c b1 κ + 1 + c b2 σ c w2 = 0.3 c w3 = 2 c v1 = 7.1 .
Finally, the function f w is computed as:

f w = g 1 + c 6 w3 g 6 + c 6 w3 1/6
with g = r + c w2 r 6 -r and r = min ν Sκ 2 d 2 , 10 .

B. Spatial discretization

The spatial discretization of the fluid equations is based on a vertex-centered finite element/finite volume formulation on unstructured meshes. It combines a HLLC upwind schemes for computing the convective fluxes and the Galerkin centered method for evaluating the viscous terms. Second order space accuracy is achieved through a piecewise linear interpolation based on the Monotonic Upwind Scheme for Conservation Law (MUSCL) procedure which uses a particular edge-based formulation with upwind elements. A specific slope limiter is employed to damp or eliminate spurious oscillations that may occur in the vicinity of discontinuities.

Finite Volume discretization

Let H a mesh of domain Ω, the vertex-centered finite volume formulation consists in associating with each vertex P i of the mesh a control volume or finite volume cell, denoted C i . Discretized domain Ω h can be written as the union of the elements or the union of the finite volume cells:

Ω h = N T i=1 K i = N S i=1 C i .
In order to discretize accurately the flow equations on highly anisotropic meshes, the dual finite volume cells are the containment sphere cells introduced in 2D by Barth 10 and generalized to 3D by Dervieux 11 instead of the classic median cells. It consists in subdividing each tetrahedron into four hexahedra cell around each vertex. The hexahedron cell vertices associated with vertex P i are (i) the middle of the three edges issued from P i , (ii) the containment circle center of the three faces containing P i , (iii) the containment sphere center of the tetrahedron and (iv) the considered vertex P i . The containment sphere cells of vertex P i is the union of all its hexahedra cells. The containment sphere center correspond the the sphere circumcenter if it falls inside the element.

The Reynolds Averaged Navier-Stokes equations integrated on the finite volume cell C i following a finite volume formulation reads (using the Green formula):

|C i | dW i dt + ∂Ci F(W i ) • n i dγ = ∂Ci S(W i ) • n i dγ , + Ci Q(W i ) dx + B.T. ( 1
)
where W i is the mean value of the solution W on the cell C i , n i is the outer normal to the finite volume cell surface ∂C i , F, S and Q are respectively the convective, viscous and source terms and B.T. represents boundary conditions terms.

Discretization of the convective terms

The integration of convective fluxes F of Equation ( 1) is done by decomposing the cell boundary in many facets

∂C ij : ∂Ci F(W i ) • n i dγ = Pj ∈V(Pi) F| ∂Cij • ∂Cij n i dγ ,
where V(P i ) is the set of all neighboring vertices linked by an edge to P i and F| ∂Cij represents the constant value of F(W ) at interface ∂C ij . The flow is calculated with a numerical flux function, denoted Φ ij :

Φ ij = Φ ij (W i , W j , n ij ) = F| ∂Cij • ∂Cij n i dγ ,
where n ij = ∂Cij n i dγ. The numerical flux function approximates the hyperbolic terms on the common boundary ∂C ij . We notice that the computation of the convective fluxes is performed mono-dimensionally in the direction normal to the boundary of the finite volume cell. Therefore, the numerical calculation of the flux function Φ ij at the interface ∂C ij is achieved by the resolution of a one-dimensional Riemann problem in the direction of the normal n ij by means of an approximate Riemann solver. In this work, the HLLC approximate Riemann solver is used for the mean flow -more details can be found [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF] -and linear upwind advection is used for the turbulent variable convection.

HLLC approximate Riemann solver. The idea of the HLLC flow solver is to consider locally a simplified Riemann problem with two intermediate states depending on the local left and right states. The simplified solution to the Riemann problem consists of a contact wave with a velocity S M and two acoustic waves, which may be either shocks or expansion fans. The acoustic waves have the smallest and the largest velocities (S L and S R , respectively) of all the waves present in the exact solution. If S L > 0 then the flow is supersonic from left to right and the upwind flux is simply defined from F(W l ) where W l is the state to the left of the discontinuity. Similarly, if S R < 0 then the flow is supersonic from right to left and the flux is defined from F(W r ) where W r is the state to the right of the discontinuity. In the more difficult subsonic case when S L < 0 < S R we have to calculate F(W * l ) or F(W * r ). Consequently, the HLLC flux is given by:

Φ hllc lr (W l , W r , n lr ) =          F(W l ) • n lr if S L > 0 F(W * l ) • n lr if S L ≤ 0 < S M F(W * r ) • n lr if S M ≤ 0 ≤ S R F(W r ) • n lr if S R < 0
Now, let us specify how W * l and W * r are evaluated. We denote by η = u • n. Assuming that η * = η * l = η * r = S M , the following evaluations are proposed [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF] (the subscript l or r are omitted for clarity):

W * = 1 S -S M    ρ (S -η) ρu (S -η) + (p * -p)n ρE (S -η) + p * S M -pη    where p * = ρ (S -η)(S M -η) + p .
A key feature of this solver is in the definition of the three waves velocity. For the contact wave we consider:

S M = ρ r η r (S R -η r ) -ρ l η l (S L -η l ) + p l -p r ρ r (S R -η r ) -ρ l (S L -η l ) ,
and the acoustic wave speeds based on Roe average:

S L = min(η l -c l , η -c) and S R = max(η r + c r , η + c) .
With such waves velocities, the HLLC Riemann solver has the following properties. It automatically (i) satisfies the entropy inequality, (ii) resolves isolated contacts exactly, (iii) resolves isolated shocks exactly, (iv) preserves positivity.

Linear convection. The turbulent variable ν is linearly convected:

Φ ρν ij (W i , W j , n ij ) = η ρν i if η > 0 η ρν j otherwise where η = 1 2 (u i • n ij + u j • n ij ) .
High-order accurate version. The MUSCL type reconstruction method has been designed to increase the order of accuracy of the scheme. [START_REF] Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF] The idea is to use extrapolated values W ij and W ji of W at the interface ∂C ij to evaluate the flux. The following approximation is performed:

Φ ij = Φ ij (W ij , W ji , n ij ) ,
with W ij and W ji which are linearly interpolated as:

W ij = W i + 1 2 (∇W ) ij • --→ P i P j and W ji = W j + 1 2 (∇W ) ji • --→ P j P i ,
where, in contrast to the original MUSCL approach, the approximate "slopes" (∇W ) ij and (∇W ) ji are defined for any edge and obtained using a combination of centered, upwind and nodal gradients. The centered gradient, which is related to edge P i P j , is defined as:

(∇W ) C ij • --→ P i P j = W j -W i .
Upwind and downwind gradients, which are also related to edge P i P j , are computed according to the definition of upstream and downstream tetrahedra of edge P i P j . These tetrahedra are respectively denoted

K ij and K ji . K ij (resp. K ji )
is the unique tetrahedron of the ball of P i (resp. P j ) the opposite face of which is crossed by the line defined by the edge P i P j . Upwind and downwind gradients are then defined for vertices P i and P j as:

(∇W ) U ij = (∇W )| Kij and (∇W ) D ij = (∇W )| Kji .
where (∇W )| K = P ∈K W P ∇φ P | K is the P 1 -Galerkin gradient on tetrahedron K. Parametrized nodal gradients are built by introducing the β-scheme:

(∇W ) ij = (1 -β)(∇W ) C ij + β (∇W ) U ij (∇W ) ji = (1 -β)(∇W ) C ij + β (∇W ) D ij ,
where β ∈ [0, 1] is a parameter controlling the amount of upwinding. For instance, the scheme is centered for β = 0 and fully upwind for β = 1.

Numerical dissipation of fourth-order: V4-scheme. The most accurate β-scheme is obtained for β = 1/3. Indeed, it can be demonstrated that this scheme is third-order for the two-dimensional linear advection on structured triangular meshes. On unstructured meshes, a second-order scheme with a fourthorder numerical dissipation is obtained. These high-order gradients are given by:

(∇W ) V 4 ij = 2 3 (∇W ) C ij + 1 3 (∇W ) U ij (∇W ) V 4 ji = 2 3 (∇W ) C ji + 1 3 (∇W ) D ij .
Numerical dissipation of sixth-order: V6-scheme. An even less dissipative scheme has been proposed. [START_REF] Debiez | Mixed-Element-Volume MUSCL methods with weak viscosity for steady and unsteady flow calculations[END_REF] It is a more complex linear combination of gradients using centered, upwind and nodal P 1 -Galerkin gradients. The nodal P 1 -Galerkin gradient of P i is related to cell C i and is computed by averaging the gradients of all the tetrahedra containing vertex P i :

(∇W ) Pi = 1 4 |C i | K∈Ci |K|(∇W )| K .
A sixth-order dissipation scheme is then obtained by considering the following high-order gradient:

(∇W ) V 6 ij = (∇W ) V 4 ij - 1 30 (∇W ) U ij -2 (∇W ) C ij + (∇W ) D ij - 2 15 (∇W ) Mi -2 (∇W ) Pi + (∇W ) Pj (∇W ) V 6 ji = (∇W ) V 4 ji - 1 30 (∇W ) D ij -2 (∇W ) C ij + (∇W ) U ij - 2 15 (∇W ) Mj -2 (∇W ) Pj + (∇W ) Pi ,
where (∇W ) Mi,j is the gradient at the points M i,j intersection of the line defined by P i P j and upwinddownwind tetrahedra. These gradients are computed by linear interpolation of the nodal gradients of faces containing M i,j .

Dervieux limiter. The previous MUSCL schemes are not monotone. Therefore, limiting functions must be coupled with the previous high-order gradient evaluations to guarantee the TVD property of the scheme. The gradient is substituted by a limited gradient denoted (∇W ) lim ij . Here, we consider the three-entries limiter introduced by Dervieux which is a generalization of the Superbee limiter:

15 if uv ≤ 0 then Lim(u, v, w) = 0 else Lim(u, v, w) = Sign(u) min(2 |u|, 2 |v|, |w|) ,
and we use:

Lim((∇W ) C ij , (∇W ) D ij , (∇W ) HO ij ) where (∇W ) HO ij is even (∇W ) V 4 ij or (∇W ) V 6 ij .

Discretization of the viscous terms

In Wolf, viscous terms are discretized with the finite element method (FEM). We evaluate viscous terms of the form:

∂Ci S(W i ) • n dγ = Pj ∈V(Pi) ∂Cij S(W i ) • n dγ + BT
where ∂C ij is the common interface between cells C i and C j , and BT are boundary terms. Let φ i the P 1 finite element basis function associated with vertex P i , we have: K ∇φ i dx = -∂Ci∩K n dγ and if we assume that S(W i ), which comes from a gradient, is constant on element K, then we get:

Pj ∈V(Pi) ∂Cij S(W i ) • n dγ = - K Pi K S(W i )| K • ∇φ i dx .
We notice that the finite element discretization is equivalent to the finite volume one. The effective computation of the previous integral then leads to the computation of integrals of the following form:

K ∇φ i ∇φ j dx = |K| ∇φ i | K ∇φ j | K .
In this expression, ∇φ i | K is the constant gradient of basis function φ i associated with vertex P i . This discretization is justified because the characteristic times associated with the diffusive terms are large as compared to the characteristic times associated with the hyperbolic (convective) terms. We now apply the FEM formulation to all convected variables that are averaged on the element and we easily verify that the components of the (Cauchy) stress tensor S(W ) are constant on each element K. For instance, the term u τ xy of the (Cauchy) stress tensor reads:

(u τ xy )| K = u| K µ| K pi∈K u i ∂φ i | K ∂y + v i ∂φ i | K ∂x .
The other terms are computed analogously.

The Spalart-Allmaras dissipation term is also discretized with the FEM:

Φ SA visc,K (W i , W j , W k , W l ) = |K| 1 σ ρ i (ν| K + ν| K ) ∇ν| K • ∇φ i | K .

Discretization of the source terms

Finally, the Spalart-Allmaras source terms (diffusion, production and destruction) are discretized by simple integration on each vertex cell:

Φ SA source (W i ) = |C i | ρ i c b1 Si νi -ρ i c w1 f w νi d i 2 + c b2 σ ρ i ∇ν i ,
where |C i | is the volume of the vertex cell and all variables are point-wise:

Si = Ω i + νi κ 2 d 2 i f v2 .

C. Boundary conditions

For presented flow simulations, three boundary conditions are required. Slip boundary conditions are imposed for bodies when the flow is considered inviscid or for symmetry. For viscous flow, no slip boundary conditions are considered for bodies. And finally, we used Steger-Warming flux to set up free-stream (external flow) conditions.

Slip condition For this boundary condition we impose weakly u • n = 0, which is done by imposing the following boundary flux:

Φ slip (W i , T ) = T Pi ∂Ci∩T F slip (W i ) • n T dγ with F slip (W i ) • n T = (0, p i n T , 0) t .
where T are boundary faces with normals n T . According to, [START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allamas turbulence model[END_REF] the slip boundary conditions for the turbulent equations is different if a wall or a symmetry plane is considered: νslipwall = 0 and ∂ νsymmetry ∂n = 0 .

No slip condition. Adiabatic conditions are considered, therefore only a null velocity is imposed strongly for this boundary condition: u = 0. The turbulent variable is also strongly imposed to zero: νnoslip = 0.

Free-stream condition. This condition imposes a free-stream uniform flow from the infinite. It applies when we have a boundary Γ ∞ for which the infinite constant state W ∞ is prescribed:

W ∞ = (ρ ∞ , (ρu) ∞ , (ρE) ∞ ) t and νfarfield ∈ [3ν ∞ , 5ν ∞ ] .
This state enables upwind fluxes at the infinite to be computed. The considered boundary fluxes are built from a decomposition following the characteristics values. We consider the Steger-Warming flux which is completely upwind on solution W i :

Φ ∞ (W i , T ) = A + (W i , n T )W i + A -(W i , n T )W ∞ where A + = |A| + A 2 and A -= |A| -A 2 .
D. Time integration

Implicit scheme

For an explicit time discretization, the semi-discretized RANS system reads:

|C i | δt n i W n+1 i -W n i = R i (W n ) where Ri(W n ) = - j∈V(i) Φ hllc (W n i , W n j , nij)+ K P i Φvisc,K (W n i , W n j , W n k , W n l )+Φsource(W n i )- T P i Φ bc,T (W n i , nT ).
For implicit time discretization, we have :

|C i | δt n i W n+1 i -W n i = R i (W n+1 )
which is linearized as:

|C i | δt n i I d - ∂R i ∂W (W n ) W n+1 i -W n i = R i (W n )
where ∂Ri ∂W (W n ) contributes the i th line of the matrix and the following linearization have been done (with the notation δW = W n+1 -W n ):

Φ hllc (W n+1 i , W n+1 j , nij) = Φ hllc (W n i , W n j , nij) + Φ hllc ∂Wi (W n i , W n j , nij)δWi + Φ hllc ∂Wj (W n i , W n j , nij)δWj Φvisc,K (W n+1 i , W n+1 j , W n+1 k , W n+1 l ) = Φvisc,K (W n i , W n j , W n k , W n l ) + Φvisc,K ∂Wi (W n i , W n j , W n k , W n l )δWi + Φvisc,K ∂Wj (W n i , W n j , W n k , W n l )δWj + Φvisc,K ∂W k (W n i , W n j , W n k , W n l )δW k + Φvisc,K ∂W l (W n i , W n j , W n k , W n l )δW l Φsource(W n+1 i ) = Φsource(W n i ) + Φsource ∂Wi (W n i )δWi Φ bc,T (W n+1 i , nT ) = Φ bc,T (W n i , nT ) + Φ bc,T ∂Wi (W n i , nT )δWi .
Terms ∂Φ ∂Wi and ∂Φ ∂Wi contribute to the matrix diagonal block and extra-diagonal block, respectively, and terms Φ contributes to the right-hand side.

We then rewrite the linearized system in compact form:

A n δW n = R n where A n = |C| δt n I - ∂R n ∂W and δW n = W n+1 -W n .

Newton's method

To solve the non-linear system, we follow the approach based on Lower-Upper Symmetric Gauss-Seidel (LU-SGS) implicit solver initially introduced by Jameson [START_REF] Jameson | Lower-Upper implicit schemes with multiple grids for the Euler equations[END_REF] and fully developed by Sharov et al. and Luo et al. [START_REF] Luo | A fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF][START_REF] Luo | An accurate, fast, matrix-free implicit method for computing unsteady flows on unstructured grids[END_REF][START_REF] Sharov | Reordering of hybrid unstructured grids for Lower-Upper Symmetric Gauss-Seidel computations[END_REF][START_REF] Sharov | Implementation of unstructured grid GMRES+LU-SGS method on shared-memory, cache-based parallel computers[END_REF] The Newton's method can be either the LU-SGS approximate factorization or the SGS relaxation or the GMRES method with LUSGS or SGS as preconditioner. The LU-SGS and SGS are very attractive because they use an edge-based data structure which can be efficiently parallelized with pthreads. [START_REF] Alauzet | On the use of space filling curves for parallel anisotropic mesh adaptation[END_REF][START_REF] Sharov | Implementation of unstructured grid GMRES+LU-SGS method on shared-memory, cache-based parallel computers[END_REF] From our experience, we have made the following -crucial -choices to solve the compressible Navier-Stokes equations.

Converging the Newton's method is important for the global convergence of the Navier-Stokes non-linear problem. Hence, an iterative method is required such as SGS or GMRES+LUSGS or GMRES+SGS a . Usually, the Newton's method iterates until the residual of the linear system is reduced by two orders of magnitude: 0.01.

a In comparison, the LU-SGS method works well for the compressible Euler equation

The choice of the renumbering also impacts strongly the convergence of the non-linear system. While Hilbert-type (space filling curve) renumbering is very efficient for cache misses and memory contention, [START_REF] Alauzet | On the use of space filling curves for parallel anisotropic mesh adaptation[END_REF][START_REF] Sharov | Implementation of unstructured grid GMRES+LU-SGS method on shared-memory, cache-based parallel computers[END_REF] Breadth-first search renumbering proves to be more effective for the convergence of the implicit method and the overall efficiency.

Luo et al. [START_REF] Luo | A fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF][START_REF] Luo | An accurate, fast, matrix-free implicit method for computing unsteady flows on unstructured grids[END_REF][START_REF] Sharov | Implementation of unstructured grid GMRES+LU-SGS method on shared-memory, cache-based parallel computers[END_REF] proposed to use a simplified flux function -a Rusanov approximate Riemann solver for the convective terms and the operator spectral radius for the viscous terms -to compute Jacobians while keeping the complex flux function for the right-hand side term. But, we observed that this modification slow down the convergence of the whole process. We found very advantageous to fully differentiate the HLLC approximate Riemann solver, [START_REF] Batten | Average-State Jacobians and Implicit Methods for Compressible Viscous and Turbulent Flows[END_REF] the FEM viscous terms and the Spalart-Allmaras source terms. [START_REF] Alauzet | Wolf user guide. An edge-based Navier-Stokes flow solver based on the MEV numerical scheme[END_REF] To achieve high efficiency, automation and robustness in the resolution of the non-linear system of algebraic equations to steady-state, it is mandatory to have a clever strategy to specify the time step. This is done by coupling local under-relaxation coefficient and local CFL, see Section D 4.

In this work, we have considered the symmetric Gauss-Seidel (SGS) relaxation. This linear system can be re-written:

(D + L)D -1 (D + U) δW n = R n + (LD -1 U) δW n
The following approximate system is used:

(D + L)D -1 (D + U) δW n = R n .
Matrix (D + L)D -1 (D + U) can inverted in two sweeps which correspond to the LU-SGS approximate factorization: Forward sweep:

(D + L) δW * = R Backward sweep: (D + U) δW = DδW * .
This sweeps can be written point wise:

δW * i = D -1 ii R i - j∈L(i) L ij δW * j δW i = δW * i -D -1 ii j∈U (i) U ij δW * j .
where L(i) (resp. U(i)) is the set of vertices with an index lower (resp. upper) that i. The lower and upper parts can be stored or not (i.e., matrix-free) as choice between efficiency or memory requirements.

In the SGS relaxation, we first zero the unknown: δW 0 = 0. Then, k max sub-iterations are made using forward and backward sweeps:

(D + L) δW k+1/2 = R -U W k (D + U) δW k+1 = R -L W k+1/2 .
or rewritten point wise:

δW k+1/2 i = D -1 ii R i - j∈L(i) L ij δW k+1/2 j - j∈U (i) U ij δW k j δW k+1 i = D -1 ii R i - j∈U (i) U ij δW k+1 j - j∈L(i) L ij δW k+1/2 j .
For one sub-iteration, the SGS method is equivalent to the LU-SGS method.

Time step computation

The local time step is computed at each vertex:

δt = CFL h 2 h (c + u ) + 2 γ ρ µ Pr + µ t Pr t
where Pr= 0.72 and Pr t = 0.9 are the Prandtl and the turbulent Prandtl constants and, µ and µ t are the (molecular dynamic) viscosity and the turbulent (molecular dynamic) viscosity.

CFL Law

Many CFL laws exist in the literature -linear, geometric, residual based, ... -but these laws generally require parameters that are difficult to establish optimally because they depend on the considered flow, the geometry and the size of the mesh. In other words, they dependent too much user's data. But, they are mandatory to achieve fast convergence in solving non-linear equations.

To avoid this issue, Luke et al. proposed a new approach [START_REF] Luke | A Step Towards Shape Shifting Algorithms: Reacting Flow Simulations Using Generalized Grids[END_REF] based on bounding the primitive variables, ρ, p and T , variations at each time step. More precisely, initially we allowed the maximal time step at each vertex, then this local time step is truncated such that the change in ρ, p and T are below a user given percentage η. But, the change in primitive variables during a given interval of time has to be estimated. A way to accomplish this is to solve an explicit time-integration step to describe a functional relationship between time and primitive variable. Notice that it is done before assembling the matrix and the truncated local time step is used to compute the mass matrix.

This method achieves a maximal efficiency as each vertex is progressing at its own optimal time step. But, that choice is made from an estimation before the linear system resolution, thus there is no guarantee on the convergence of the Newton's method.

Another attractive approach has been proposed by Burgess and Glasby 26 which couples under-relaxation coefficient and dynamic CFL. Here, the solution is analyzed at each step of the Newton's method (after solving the linear system) and before updating the solution. First, the change in primitive variables, ρ and p, is again controlled by a user given percentage η and defines a under-relaxation coefficient ω n at each step of the process. This global coefficient is then applied to the solution evolution: W n+1 = W n + ω n δW . Then, the CFL value is updated depending on that under-relaxation coefficient:

CF L n+1 =      0.1 CF L n if ω n < 0.1 CF L n if 0.1 ≤ ω n < 1 αCF L n + β if ω n = 1
where we choose α = 1 and β = 1 for a linear increase or α = 2 and β = 0 for a geometric increase. This adaptive CFL, thus time step, is attractive because it is based on the behavior of the Newton's method. To improve even more the robustness of the method, they propose to set the solution update to zero when the value of ω n is less than 0.1, i.e., ω n = 0. This approach is extremely robust because if the Newton's method diverges, the current step is cancelled and the time step, via the CFL, is automatically reduced. But the considered criterium is global and hence one bad vertex in the mesh can kill the overall efficiency by not allowing the CFL to grow.

In Wolf, we consider an hybrid method having the efficiency of the first method and the robustness of the second one. We proceed exactly like the second approach but the under-relaxation coefficient is set locally, i.e., vertex-wise, and each vertex is supplied with its own CFL coefficient which evolves with respect to its own under-relaxation coefficient. Thus, we have a local time step and a local CFL for each vertex.

E. Validation of the flow solver

In this Section, we compare Wolf, our in-house flow solver, to other codes for several Reynolds-Averaged Navier-Stokes (RANS) simulations. These comparisons ensure that the turbulence model has correctly been implemented in Wolf. This validation step is mandatory to confidently rely on the results of more complex simulations.

We compared Wolf to other codes for a large panel of verification and validation cases including flat plates, bumps, steps, airfoils... and we obtained relevant results. The material we used for comparison includes the test cases database released by NASA, [START_REF] Rumsey | Turbulence Modeling Resource[END_REF] which provides flow conditions, grids and results from CFL3D 28 for comparison purpose.

We exhibit one verification case and two validation cases: a turbulent bump-in channel, a NACA 0012 airfoil and a transsonic RAE 2822 airfoil.

Bump-in-channel

This test case is quite similar to a turbulent flat plate, except that the lower wall is curved. It was run at Mach M = 0.2 and at a Reynolds number of 3 Million based on length 1 of the mesh.

Speed isolines are depicted in Figure 1 and extracted velocity profiles in Figure 2. The results obtained with Wolf are similar to those from CFL3D and FUN3D. 3 and4. We compare pressure coefficients and skin friction coefficients for various angles-of-attack in Figure 5. RAE 2822 Here we study a transsonic flow (Mach 0.7) around the RAE 2822 airfoil. [START_REF] Slater | RAE 2822 Transonic Airfoil: Study 1[END_REF] The Reynolds number is 6.5 Million based on airfoil chord c = 1.5 and the angle-of-attack is α = 2.31 deg. Mach isolines are depicted in Figure 6. We compare [START_REF] Cook | Aerofoil RAE 2822: Pressure Distributions, and Boundary Layer and Wake Measurements[END_REF] the pressure coefficients to experimental data and to Wind-US in Figure 7.

II. Meshing technologies

In this section, we review the main features of the adaptive mesh generator feflo.a. It is a metricbased mesh generator meaning that a metric tensor is used to prescribe the anisotropy in the domain. All the operators are based on local mesh modifications in order to guarantee that an adapted mesh is always provided on output. In the sequel, we describe respectively how the volume, surface, boundary layer and anisotropic cartesian mesh are generated.

A. Volume and surface cavity-based operators

The generator is based on the concept of unit mesh, which extends the notion of uniform mesh to anisotropic meshes. The anisotropy is prescribed through a metric field, so we suppose that a couple mesh and metric (H, M) is always provided. In this context, the adaptive mesh procedure aims at generating a unit mesh with respect to M. A mesh is said to be unit when composed of almost unit-length edges. The length of an edge e = [A, B] in M is then evaluated with:

M (A, B) = 1 0 t AB M((1 -t)A + tB) AB.
A first step in the adaptive mesh generation process thus consists in performing point insertions and edge collapses recursively until all edges have a length comprised in [1/ √ 2, √ 2]. Then, an optimization procedure is added based on point smoothing, edges and faces swapping. During this step, we try to improve the overall quality of the mesh (in comparison to the unit mesh step where the convergence to a prescribed length distribution is the main concern). The quality is given in M by:

Q M (K) = 36 3 1 3 |K| 2 3 M 6 i=1 2 M (e i ) ∈ [0, 1],
and a perfect unit element has a quality of 1.

In feflo.a, all the aforementioned mesh modification operators are based on an unique cavity-based operator revisited in metric-based context. This operator is inherited from Delaunay incremental insertion and can be written in a compact form as:

H k+1 = H k -C P + B P ,
where C P is the cavity of P and B P the ball of P . In, [START_REF] Loseille | Serial and Parallel Mesh Modification Through A Unique Cavity-Based Primitive[END_REF] we derive several initializations and modifications of C p to perform insertions, swaps, collapses and point smoothing. This operator also performs surface-based operations. In that case, a local surface approximation is computed based on a provided background surface mesh.

B. Hybrid cavity-based operators

The generation of a quasi-structured boundary layer mesh is based on the insertion of layers of hybrid entities (prisms when a triangulated surface is provided). The creation of hybrid entities is also based on a modified cavity-based operator. The modification consists in providing a list of initial faces from which the boundary layer must be extruded. This additional information is preserved and propagated into the domain for the subsequent layers. More precisely, starting from the initial surface mesh S = (F i ) i , a set of normals are computed from (F i ) i . Multi-normals are added in the case of singularities (an open ridge for example). We then list all the faces associated to a normal: (n i , F k1 , . . . , F kn ). Then, the insertion of P is given by the following cavity operator:

H k+1 = H k -(C P -K) + B P S k+1 = S k -∪ i F ki + ∪ i Fki ,
where K is the set of hybrid entities composing the previous or current layers and ∪ i Fki are the new faces connected to P . This set of element defines a constraint on the cavity operator in order to avoid the removal of any elements defining the quasi-structured mesh.

In Figure 8, we represent the boundary layer mesh extrusion around a high lift geometry. In this example, 100 000 prisms/second are inserted on a Intel Core i7 at 2,7 Ghz. Note that no transition is applied, see Section B for detail smoothing and transition between the boundary layer and the remaining volume mesh.

C. Cartesian operators

Cartesian meshes are usually devised to generate quasi-structured meshes. In CFD, this meshes usually allows us to reduce transverse dissipation and to reduce the number of degrees of freedom by deactivating some edges/points when using appropriate finite volume cell definition (Barth cell for instance). However, the generation of unstructured cartesian meshes is usually combined with a frontal approach where a front is propagating into the domain while combining iteratively sequence of of (paving) tetrahedra in order to create regular structures. However, only uniform meshes are generally generated in 3D with this approach. In feflo.a, the cartesian operator relies on the previous cavity operator. The difference occurs in the insertion pattern. The new points are created along the eigenvector directions of the input metric. This choice allows us to favor orthogonally and alignment even with an anisotropic metric. When an isotropic metric is provided, alignment is along the natural

III. Application to mesh adaptation for viscous flow simulations

Contrary to inviscid flow simulations, viscous flow simulations require a specific meshing of the boundary layers in the near wall regions. This is due notably to the dramatic variation in the normal direction of variables such as the velocity. The generation of quasi-structured meshes in these regions has proven to be a reliable approach. [START_REF] Aubry | Generation of Viscous Grids With Ridges and Corners[END_REF][START_REF] Bottasso | A procedure for tetrahedral boundary layer mesh generation[END_REF][START_REF] Garimella | Boundary layer mesh generation for viscous flow simulations[END_REF][START_REF] Hassan | Unstructured Tetrahedral Mesh Generation for Three-Dimensional Viscous Flows[END_REF][START_REF] Ito | An approach to generate high quality unstructured hybrid meshes[END_REF][START_REF] Löhner | Generation of unstructured grids suitable for RANS calculations[END_REF][START_REF] Marcum | Adaptive Unstructured Grid Generation for Viscous Flow Applications[END_REF][START_REF] Pirzadeh | Viscous unstructured three dimensional grids by the advancing-layers method[END_REF][START_REF] Leicht | Error Estimation and Anisotropic Mesh Refinement for 3D Laminar Aerodynamic Flow Simulations[END_REF][START_REF] Park | Parallel Anisotropic Tetrahedra Adaption[END_REF][START_REF] Chalot | Industrial Aerodynamics[END_REF] Among the difficulties of boundary layer mesh generation include (i) having a smooth transition with the fully unstructured mesh used for the rest of the domain, (ii) its integration in the mesh adaptation loop described in the Introduction and (iii) the complex geometries often encountered in aeronautics.

In this section, we mainly focus on (ii) with a consideration for (i) and (iii). To this end, we combine the work introduced in Sections I and II to study how to couple anisotropic mesh adaptation with boundary layer mesh generation for complex geometries. Starting from a fully unstructured mesh adaptation scheme, improvements to generate and adapt quasi-structured boundary layer meshes are presented.

A. Fully unstructured mesh adaptation

In this first method, an initial semi-structured boundary layer mesh is mapped onto a metric tensor field. This metric field contains information on elements' orientation and size. Any standard mesh estimate is then used to generate a unit mesh with respect to the metric tensor field.

We apply this strategy to capturing shock/boundary layer interactions, see Figure 9. A supersonic flow (Mach 1.4) is applied across a double wedge wing. It generates a shock wave which interacts with a boundary layer at the bottom of the domain. We aim at capturing these interactions. This simulation leads to the following observations:

• Generating a semi-structured boundary layer mesh extruded from the surface mesh gives only the require accuracy for the smaller layers. Indeed, the distance of the bottom of the shock from the viscous plate is around 10 -3 whereas the initial height of the uniform boundary layer mesh was at 0.2. Consequently, it is not possible to keep a constant boundary layer mesh when an interaction occurs.

• This approach is completely generic and robust and can handle complex geometries. However, if the shock/boundary layer interaction is automatically handled, the impact of having a fully unstructured mesh is not yet analyzed in term of solution accuracy and solver stability in the viscous area. Consequently, it seems also interesting to derive a method to generate structured mesh for the smaller layers (at least) while preserving (upper) anisotropic refinements.

This method has the advantage of simplifying the mesh adaptation mesh: one does not need a specific strategy for the near-wall regions. However, if the shock/boundary layer interaction is automatically handled, the impact of having a fully unstructured mesh is not yet analyzed in term of solution accuracy and solver stability in the viscous area. Consequently, it seems also interesting to derive a method to generate structured mesh for the smaller layers (at least) while preserving (upper) anisotropic refinements: this work is presented in B and C.

B. Unstructured mesh adaptation with fixed boundary layer mesh

It is not clearly established wether meshing near-wall regions with fully-unstructured meshes is suited for capturing viscous phenomena. Therefore, one might want to use semi-structured boundary layer meshes which have proven to be reliable. In this Subsection, we introduce a method [START_REF] Park | Turbulent output-based anisotropic adaptation[END_REF] in which a boundary layer mesh is generated once at the beginning of the simulation. During the mesh adaptation loop, only the rest of the domain is adapted.

We describe the main steps of the adaptation process. First, the boundary layer mesh is generated according to a given wall spacing, a given number of layers and a given size gradation. Then, it is tagged and extracted: vertices on the interface between the two meshes are duplicated and new surface elements are created. We obtain two meshes: one boundary layer mesh and one we call domain mesh. The boundary layer mesh remains unchanged while the domain mesh is adapted. During the mesh adaptation of the domain, we make sure the vertices and triangles at the interface are not modified in order to be able to gather the two mesh parts.

We applied this method to a plume exhaust (see Figure 10) at Mach 2.2 with a Reynolds num-ber of 1.8 Million. The final mesh and mach isovalues are depicted in Figures 11 and12. See Figure 13 for cuts in the volume.

To adapt the domain mesh, the solution we based the adaptation on is interpolated on the domain mesh and an hessian-based metric is computed. If there are several solutions of interest we want to adapt, we compute a metric for each one of them and perform a metric intersection. The metric field obtained for the domain will be used for a mesh adaptation.

A smooth transition between the boundary layer mesh and the rest of the domain is mandatory because it ensures a good behavior of the flow solver. In other words, we want no discontinuity in the elements' sizes. To do so, the natural geometric field of the boundary layer is propagated in the domain and intersected with the hessian-based metric. This method has a major weakness: the skin of the plane is not adapted. Therefore it is not suited for capturing strong pressure variations on the wings for instance. Furthermore, this method requires too many mesh manipulations (splitting, gathering, metric gradation etc...). We want to devise a more simple approach with skin adaptation. Improvements are introduced in Subsection C.

C. Hybrid mesh with skin adaptation

This method is quite similar to the one introduced in Subsection B, except that the boundary layer mesh is re-generated at each iteration in the mesh adaptation loop. We proceed as follows. Given an initial fully unstructured mesh, a semi-structured boundary layer mesh is extruded. A flow solution and a metric tensor field are then computed. We then perform a fully-unstructured remeshing of the whole mesh according to that metric field.

The simulation to the following observations:

• This method allows to better capture anisotropic phenomena on the skin.

• It requires a strong robustness of the remesher. Re-generating a quasi-structured mesh from an anisotropic mesh presents challenging local configurations.

D. Cartesian approach

In the Cartesian approach, the boundary layer mesh is generated at the same time as the adapted volume mesh. To do so, a boundary layer metric is computed based on the distance to the viscous body. In the normal directions, the size is based on the desired size progression while in the tangential direction, the size is extrapolated from the surface mesh. We illustrate the cartesian operator with an isotropic metric field around a f117 geometry at high angle of attack (20 degrees). This mesh is depicted in Figure 14 and is composed of 7 532 632 vertices and 4 572 1814 tetrahedra. It was generated in 7 min on an Intel Core i7 at 2.7 Ghz with 16 Gb of RAM. It is used to study the vortical flows generated by the delta wings. We expect to have a high level of details when observing the wake generated by the aircraft. From a flow solver point of view, this meshes enhances the high-order features of the solver (6th order for linear advection on structured meshes) by removing some high-order truncation error terms. Gradients of the solution are also better approximated.

If the previous example is isotropic, this approach works as well in an anisotropic context. We consider a transonic flow computation around a generic Falcon geometry at mach 0.8 with an angle of attack of 3 degrees. We control the interpolation error of the Mach number in L 2 norms, the final mesh is obtained after 30 iteration and is composed of 1 110 735 vertices and 6 546 789, see Figure 15. The total cpu time is 40 minutes on an Intel Core i7 at 2.7Ghz.

IV. Conclusion and future work

The results of the code validation study led on the flow solver Wolf are satisfying. However, as a relatively young code several improvements remain to be achieved. We have recently implemented a full multigrid method and we are currently analyzing the results in order to increase the convergence gain. Performing mesh adaptations on every grid levels could be interesting. We intend to parallelize the solver in MPI in order to run simulations of cases of high complexity such as the high ligt prediction. [START_REF] Rumsey | High Lift Prediction Workshop[END_REF] We have experimented four different adaptive approaches for viscous simulations. The fully unstructured approach doesn't require any specific mesh nor metric manipulation (i.e. splitting, metric intersection, boundary layer extrusion etc.) for near-wall regions. That makes the adaptive process simple but it does not preserve the structure of the smaller layers. The structured mesh of the near-wall regions is preserved by the method presented in B and C, but at the cost of simplicity. We seek a completely generic method that preserves the structure. A first draw of this method is introduced inD. Based on the local cartesian operator presented in Section II, it

V. Acknowledgements

This work was partially funded by the Airbus Group Foundation.

  with d the distance to nearest wall which is computed for each vertex at the beginning of the simulation. The constants are σ = 2 3 c b1 = 0.1355 c b2 = 0.622 κ = 0.41

Figure 1 .

 1 Figure 1. Speed isolines around the bump Figure 2. Bump: Velocity profiles

Figure 3 .

 3 Figure 3. NACA0012: Pressure isolines (α = 0).Figure 4. NACA0012: Speed isolines (α = 0).

Figure 4 .

 4 Figure 3. NACA0012: Pressure isolines (α = 0).Figure 4. NACA0012: Speed isolines (α = 0).

Figure 5 .

 5 Figure 5. NACA0012: Pressure coefficient Cp and skin friction coefficient C f for α = 0, 10, 15 deg. Comparison with experimental data (Cp) and CFL3D (Cp and C f ).

Figure 6 .

 6 Figure 6. RAE 2822: Mach isolines.

Figure 7 .

 7 Figure 7. Comparison of the pressure coefficient around the RAE 2822 airfoil.

Figure 8 .

 8 Figure 8. Boundary layer mesh generation around a high lift geometry.

Figure 9 .

 9 Figure 9.Shock/boundary layer interaction test case. Top, from left to right: Mach number isolines and final adapted mesh. Bottom: Velocity isolines in the interaction region and corresponding adapted mesh.

Figure 10 .

 10 Figure 10. Initial domain of the plume exhaust.

Figure 11 .

 11 Figure 11. Final surface mesh of the domain.

Figure 12 .

 12 Figure 12. Surface mesh near the exhaust (left) and Mach iso-values (right).

Figure 13 .

 13 Figure 13. Cuts in the volume, red elements are part of the structured extruded layers

Figure 14 .

 14 Figure 14. Cartesian mesh around a f117 geometry: Top left, cut in the volume mesh, then snapshots of the Mach number at different time steps showing the vortices and the wake of the aircraft.

Figure 15 .

 15 Figure 15. Final adapted surface mesh on the falcon geometry (left) and density iso-values (right).

Figure 16 .

 16 Figure 16. Different cuts in the final adapted volume mesh showing the quasi-aligned anisotropic refinements: in the wake (top and bottom left), near the lambda shock (top right) and 100, behind the falcon (bottom right).

  τ xy , τ xz , uτ xx + vτ xy + wτ xz + λT x ) T S 2 (W ) = (0, τ xy , τ yy , τ yz , uτ xy + vτ yy + wτ yz + λT y ) T S 3 (W ) = (0, τ xz , τ yz , τ zz , uτ xz + vτ yz + wτ zz + λT z ) T