Jessica Sari

Francesco Cremonesi

Mehdi Khalloufi

François Cauneau

Philippe Meliga

Youssef Mesri

Elie Hachem
email: elie.hachem@mines-paristech.fr

Anisotropic adaptive stabilized finite element solver for RANS models

Keywords: Anisotropic mesh adaptation, Turbulent flows, Stabilized finite element method, Spalart-Allmaras, Boundary layers, Reynolds Averaged Navier-Stokes

INTRODUCTION

In the aerospace industry, aerodynamic efficiency is one of the main components of the overall aircraft performance, and has thus a significant impact on the commercial success of aircraft.

Consequently, airplane and airship designers are often concerned with accurately predicting aerodynamic forces and moments as a function of design parameters. In particular, the capability to determine the drag and lift coefficients for a given configuration is instrumental in the process of aircraft design. For this purpose, it is essential to reach a high level of accuracy and reliability in the description of the physical phenomena involved. For instance, one drag count (i.e. 10 -4 in drag coefficient) is often cited as the desired accuracy for the drag prediction of subsonic transport aircraft [START_REF] Michael | Statistical Analysis of Computational Fluid Dynamics Solutions from the Drag Prediction Workshop[END_REF] or the shape optimization of boats [START_REF] Parolini | Mathematical models and numerical simulations for the americas cup[END_REF][START_REF] Formaggia | Fluidstructure interaction problems in free surface flows: Application to boat dynamics[END_REF]. However, the measurement of such quantities through experimental data can be often time and budget consuming as well as difficult to generalize, as such quantities depend strictly on the geometry at hand. For this reason, numerical simulations have become very important for the conceptual and design phases in the aerospace industry [START_REF] Jameson | A perspective on computational algorithms for aerodynamic analysis and design[END_REF].

The high accuracy requirement for the prediction of integral coefficients in aerodynamic simulations implies the necessity to correctly resolve local flow features such as boundary layers or flow separation. For this purpose, locally-refined meshes are usually used. They incur a computational cost that is several orders of magnitude lower than globally-refined meshes [START_REF] Rannacher | Finite element methods for the incompressible Navier-Stokes equations[END_REF]. However, the generation of locally-refined grids is an extremely time-and effort-consuming process that relies on the know-how and intuition of experts. It is not uncommon that engineers and technicians must re-create adapted meshes when simulations prove locally under-resolved. Mesh generation may amount to 80% of the overall analysis time in the automotive, aerospace, and ship building industries [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF].

One remedy to these challenges lies in algorithms for automatic adaptation of unstructured meshes. A large variety of methods were proposed to optimize the compromise between accuracy and computational cost in CFD (Computational Fluid Dynamics). One of the most successful approaches involves a metric field that locally modifies the evaluation of spatial lengths depending on an error estimator. With this technique, the mesh can be adapted in a classical manner, that is, the local mesh size can be varied while maintaining element shapes as regular as possible. Above all, it also makes it possible to perform anisotropic adaptation: the element shapes can be controlled to fulfill a local mesh size requirement that varies depending on the direction. This capability is necessary for obtaining optimal meshes in fluid dynamics problems, where anisotropic flow features such as boundary layers are ubiquitous. Nevertheless, it implies that the numerical scheme perform adequately with highly stretched elements or cells, which is not the case of all numerical methods.

Methods exploiting a priori error bounds through the evaluation of the Hessian or the gradient of the solution have been developed in both the isotropic and anisotropic setting in [START_REF] Frey | Anisotropic mesh adaptation for cfd computations[END_REF][START_REF] ¸ois Remacle | Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods[END_REF][START_REF] Castro | New progress in anisotropic grid adaptation for inviscid and viscous flows simulations[END_REF]. Another class of anisotropic error estimates have also been obtained in [START_REF] Formaggia | New anisotropic a priori error estimates[END_REF][START_REF] Formaggia | Anisotropic error estimates for elliptic problems[END_REF][START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF], leading to the standardization of the anisotropic metric construction procedure. However, dynamic anisotropic mesh adaptation remains a challenge in the simulation of high Reynolds number flows past complex geometries. In such cases, Reynolds-Averaged Navier-Stokes (RANS) models are most often used, due to the prohibitive cost of turbulent scale-resolving approaches. With such models, local flow features such as boundary layers and separation need to be well captured by mesh refinement, and therefore the use of a priori adapted meshes often proves inefficient.

In this paper, we propose to solve RANS models by a stabilized finite element method combined with an efficient dynamic anisotropic mesh adaptation procedure. In particular, we focus on the Spalart-Allmaras [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] turbulence model. The mesh adaptation, driven by an efficient error estimator, automatically captures thin boundary layers and flow separation points at high Reynolds number (Re ∼ 10 6 ). The procedure makes it possible to accurately simulate three-dimensional turbulent flows around complex geometries with a reasonable computational cost.

In this work, particular attention is paid to the automatic character of the mesh adaptation process, with the objective of providing a useful and practical tool for aerodynamic applications. The proposed algorithm enables adaptive meshing under the constraint of a fixed number of nodes, which spares the end user from dealing with a case-dependent parameter linked to the error analysis. Moreover, the choice of the appropriate input for the error estimator is usually difficult when unknowns of different natures are solved simultaneously. Thus, we describe in this paper a multicriteria functional that estimates the error on a combination of several variables including the velocity norm, all the velocity components as well as the turbulent viscosity.

In principle, the mesh adaptation procedure can be coupled to different spatial discretization schemes such as Finite Volumes or Finite Elements. However, the numerical method must be able to handle highly stretched simplicial elements. It is notoriously complicated to find schemes that exhibit the level of robustness required by automatic anisotropic mesh adaptation, especially among Finite Volume methods. We show here that stabilized Finite Element methods using appropriate stabilization parameters are adequate to deal with elements of aspect ratio up to the order of 1000. A Variational Multiscale formulation for the incompressible Navier-Stokes equations, as well as a Streamline Upwind Petrov-Galerking method used to solve the Spallart-Allmaras turbulence equation, are presented for this purpose [START_REF] Erb | Edgebased finite element implementation of the residual-based variational multiscale method[END_REF].

The proposed framework is first verified using the well known NACA0012 airfoil benchmark. The results obtained for turbulent flows around cylinders with three different geometries are also presented and compared to data from the literature. The solver is then validated by comparing numerical predictions for the flow past an Unmanned Aerial Vehicle (UAV) with experimental measurements.

This paper is organized as follows: we start in section 2 with the description of the anisotropic mesh adaptation technique. Then, we present in Section 3 the governing equations and their discretization. Section 4 describes the numerical test cases showing the efficiency and the accuracy of the proposed method. Finally, Section 5 is dedicated to the conclusion and the future work.

ANISOTROPIC MESH ADAPTATION

The mesh adaptation strategy presented in this work relies on the a posteriori definition of a metric field able, in some specified sense, to drive the re-meshing procedure so that the interpolation error on the finite element solution is minimized. It relies on a statistical representation of the distribution of edges sharing a vertex, a quantity we call length distribution tensor. In order to relate the length distribution tensor to the interpolation error, following [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF] we define an edge based error estimator based on a gradient recovery procedure. Once the optimal metric has been obtained, the mesh generation and adaptation procedure described in [START_REF] Mesri | On optimal simplicial 3d meshes for minimizing the hessian-based errors[END_REF], based on a topological representation, is applied to obtain the new mesh.

Definition of the length distribution tensor: a statistical representation

Let Ω ⊂ R d be a polygonal domain, we consider a triangulation Ω = K∈K K where K is a simplex such as a triangle or tetrahedron. Let Γ(i) be the "patch" associated to a mesh vertex x i , defined as the set of nodes x j which share one edge with x i , and let us denote by x ij the edge connecting x i to x j as in Figure 1. The problem of finding a unitary metric M i associated to the i th node can be formulated as the least squares problem [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF]

M i = argmin M∈R d×d sym   j∈Γ(i) Mx ij • x ij -|Γ(i)|   2 .
(

Provided that the vertices of |Γ(i)| form at least d non co-linear edges with vertex x i (which is the case in a valid mesh), then an approximate solution of (1) is given by

M i = 1 d X i -1 , (2)
where, denoting by ⊗ the tensor product between two vectors, we have introduced the length distribution tensor

X i = 1 |Γ(i)|   j∈Γ(i) x ij ⊗ x ij   (3)
whose purpose is to give an "average" representation of the distribution of edges in the patch.

Gradient recovery error estimator

Let v h be a P1 finite element approximation obtained by applying the Lagrange interpolation operator to a regular function v ∈ C 2 (Ω). At node x i , we seek the recovered gradient g i of v h defined by

g i = argmin g∈R d j∈Γ(i) (g -∇v h ) • x ij 2 = argmin g∈R d j∈Γ(i) g • x ij -v h (x j ) -v h (x i ) 2 . ( 4
)
The solution to (4) can be expressed in terms of the length distribution tensor introduced before as

g i = 1 |Γ(i)| X i -1 j∈Γ(i) v h (x j ) -v h (x i ) x ij . (5)
It can be shown [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF] that the quantity |g i • x ij | gives a second order accurate approximation of the second derivative of v along the edge x ij . Motivated by the fact that, for P 1 finite elements on anisotropic meshes, edge residuals dominate a posteriori errors, as proved in [START_REF] Kunert | Edge residuals dominate a posteriori error estimates for linear finite element methods on anisotropic triangular and tetrahedral meshes[END_REF], it is therefore suitable to define an error indicator function associated to the edge x ij as

e ij = g ij • x ij . (6)
Moreover, this quantity can be easily extended to account for several sources of error, instead of just the scalar field v h , by applying formula [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] to each component separately.

Metric construction

In order to relate the error indicator defined in (6) to a metric suitable for a mesh adaptation procedure, we introduce the concept of stretching the factor s ij defined as the ratio between the length of the edge x ij after the adaptation procedure and before the adaptation procedure. The new metric, denoted by M i , will then be given by

M i = 1 d X i -1 ; X i = 1 |Γ(i)|   j∈Γ(i) (s ij ) 2 x ij ⊗ x ij   (7)
To relate the metric to the interpolation error, following [START_REF] Jannoun | Anisotropic meshing with time-stepping control for unsteady convection-dominated problems[END_REF] we set a target total number of nodes N . Because of the quadratic behavior of the error as a function of the scaling factor, denoting by ẽij the quantity defined in (6) computed after the mesh adaptation process we have

s ij = ẽij e ij 1/2 . ( 8
)
Moreover, the number of nodes in the new mesh resulting from the application of the scaling factor s ij to the edge will be roughly equal to

N ij = 1 s ij , (9)
so that the total contribution of node i (in the old mesh) to the number of nodes in the new mesh will be given by

N i = det   X i -1 j∈Γ(i) N ij x ij ⊗ x ij   (10)
By combining [START_REF] Castro | New progress in anisotropic grid adaptation for inviscid and viscous flows simulations[END_REF] and [START_REF] Codina | Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods[END_REF], it is possible to consider N i in (10) is a function of ẽij . Assuming that the total error is equidistributed among all edges such that each edge contributes a constant error e to the total, one can see that

N ij = e ij e 1/2
, which results in the relation

N i (e) = e -d 2 N i (1). ( 11
)
This article is protected by copyright. All rights reserved.

By summing over all the nodes of the old mesh, an expression for the total error as a function of the number of nodes in the new mesh is obtained. Inverting this equality yields

e =   i N i (1) N   2 d , (12)
and injecting this into (8) provides an expression for the scaling factors

s ij = i N i (1) 1 d N 1 d (e ij ) 1 2
.

(13)

Mesh adaption criteria

In the preceding sections, we have explained how a metric filed driving the mesh adaptation can be derived from an edge-based error estimator under the constraint of a constant number of nodes. Nevertheless, this work is dedicated to the solution of Fluid Dynamics problems, which involve several unknowns, namely the velocity field v and the turbulent eddy viscosity ν (see Section 3). The common way to adapt a mesh to several variables is to compute the metrics corresponding to each of them and then to produce a unique metric by intersection. The resulting metric should yield the largest mesh size in any direction that still fulfills the size constraint given by all the metrics. The intersection operation between several metric fields is not uniquely defined. The most common technique, consisting of a sequence of simultaneous reductions of two metrics [START_REF] Alauzet | Adaptation de maillage anisotrope en trois dimension[END_REF], does not provide the optimal metric, is not commutative and incurs a relatively high computational cost. Robust techniques, based on the optimization problems, are even more costly.

In this work, we simplify this operation by using one field that accounts for different variables. Applying the error estimator (6) makes it possible to define a metric that to account for several sources of error. A vector of error sources is locally defined in two dimensions as

Y(x i ) =    v i x v i , v i y v i , v i max j v j , νi max j νj    . ( 14
)
In three dimensions, a component v i z / v i is added. The field used as input for the error estimator [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] is the norm

k Y 2 k 1/2 of Y, with Y k the components of Y.
This definition measures the error in the norm as well as in the direction of the velocity vector v, in addition to the error on the turbulent eddy viscosity ν. Because all fields are normalized (the velocity components v x , v y and v z by the local velocity norm, the velocity norm v and the turbulent eddy viscosity ν by their respective global maxima), a field that is much larger in absolute value does not dominate the error estimator, and the variations of all variables are fairly taken into account.

GOVERNING EQUATIONS AND NUMERICAL SCHEME

In this section, we present the governing equations, namely the unsteady incompressible Navier-Stokes equations coupled to the Spalart-Allmaras turbulence model, and their discretization. Time splitting is applied to the system, so that the Navier-Stokes and the Spalart-Allmaras transport equation can be solved separately using different numerical methods. At each time step, the Navier-Stokes equations are solved in a first stage using the value of the turbulent variable from the preceding time step, and the resulting velocity filed is used to solve the Spallart-Allmaras equation in a second stage. The stabilized Finite Element scheme for the Navier-Stokes equations is derived from a variational multiscale point of view. Both the velocity and the pressure spaces are enriched which cures the spurious oscillations in the convection-dominated regime and deals with the pressure instability. The Spalart-Allmaras transport equation is solved through a Streamline Upwind Petrov-Galerkin (SUPG) method. In both schemes, particular attention is paid to the determination of the stabilization parameters in the presence of mesh anisotropy, through the use of a directional element size.

The Navier-Stokes equations

To fix notations, let Ω ⊂ R d be the fluid domain, where d is the space dimension, and ∂Ω its boundary. The strong form of the incompressible Navier Stokes equations reads:

ρ (∂ t v + v • ∇v) -∇ • σ = f ∇ • v = 0 ( 15
)
where t ∈ [0, T ] is the time, v(x, t) the velocity, p(x, t) the pressure and ρ the density. The Cauchy stress tensor for a Newtonian fluid is given by:

σ = 2µ (v) -p I d , (16)
with I d the d-dimensional identity tensor and µ the dynamic viscosity. In order to close the problem, Equations ( 15) are subject to the boundary and initial conditions to be specified later.

The Reynolds-Averaged Navier-Stokes (RANS) equations are obtained by applying a timeaveraging procedure to the Navier-Stokes equations, which gives rise to an additional stress term (the Reynolds stress tensor) taking into account the effect of turbulent fluctuations on the averaged flow. The Boussinesq approximation can then be used to model this additional stress as a turbulent eddy viscosity µ t . The RANS equations, with Reynolds-averaged velocity and pressure as unknowns, have the same form as the Navier-Stokes equations, except that the molecular viscosity µ is replaced by the sum µ + µ t . The eddy viscosity µ t is obtained from a model involving additional partial differential equations. We choose for this work the one-equation model described in Section 3.2. Hereafter, we use the same notations for the Navier-Stokes and RANS equations and denote the total viscosity as µ for the sake of simplicity.

The weak form of problem [START_REF] Förster | Stabilized finite element formulation for incompressible flow on distorted meshes[END_REF] combined with ( 16) can be obtained by multiplication of a test function and integration by parts. Let H 1 (Ω) be the Sobolev space of square integrable functions whose distributional derivatives are square integrable, and let V ⊂ H 1 (Ω) d be a functional space properly chosen according to the boundary conditions. Finally, let Q = q ∈ L 2 (Ω) : Ω q = 0 . By denoting (•, •) the scalar product of the space L 2 (Ω), the weak form of problem [START_REF] Förster | Stabilized finite element formulation for incompressible flow on distorted meshes[END_REF] on ∂Ω reads, under the assumption of homogeneous Dirichlet boundary conditions:

       Find (v, p) ∈ V × Q such that: ρ [(∂ t v, w) + (v • ∇v, w)] + (2µ (v) : (w)) -(p, ∇ • w) = (f , w) , ∀w ∈ V (∇ • v, q) = 0, ∀q ∈ Q. ( 17
)
where ρ and µ are the density and the dynamic viscosity, respectively. Based on a mesh K h of Ω made of N el elements K, the functional spaces for the velocity V and for the pressure Q are approximated by the finite dimensional spaces V h and Q h respectively. It is well known that the stability of the semi-discrete formulation requires an appropriate choice of the finite element spaces V h and Q h , that must to fulfill a compatibility condition [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. Accordingly, the standard Galerkin method with the P1/P1 element (i.e. the same piecewise linear space for V h and Q h ) is not stable. Moreover, convection-dominant problems (i.e. problems where the convection term v • ∇v is much larger than the diffusion term ∇ • (2µ )) also lead to a loss of coercivity in formulation [START_REF] Hachem | Stabilized finite element method for incompressible flows with high reynolds number[END_REF]. This phenomenon manifests itself as oscillations that pollute the solution. In this work, we use a Variational MultiScale method [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origin of stabilized methods[END_REF] which circumvents both problems through a Petrov-Galerkin approach. The basic idea is to consider that the unknowns can be split into two components, a coarse one and a fine one, corresponding to different scales or levels of resolution. First, we solve the fine scales in an approximate manner and then we replace their effect into the large-scale equation. We present here only an outline of the method, and the reader is referred to [START_REF] Hachem | Stabilized finite element method for incompressible flows with high reynolds number[END_REF] for extensive details about the formulation. Let us split the velocity and the pressure fields into resolvable coarse-scale and unresolved finescale components: v = v h + v and p = p h + p . The same decomposition can be applied to the weighting functions: w = w h + w and q = q h + q . Subscript h is used hereafter to denote the finite element (coarse) component, whereas the prime is used for the so called subgrid scale (fine) component of the unknowns. The enrichment of the functional spaces is performed as follows:

V = V h ⊕ V , V 0 = V h,0 ⊕ V 0 and Q = Q h ⊕ Q .
Thus, the finite element approximation for the time-dependent Navier-Stokes problem reads:

                   Find(v, p) ∈ V × Q such that: ρ (∂ t (v h + v ), (w h + w )) + ρ ((v h + v ) • ∇(v h + v ), (w h + w )) + (2µε ε ε(v h + v ) : ε ε ε(w h + w )) -((p h + p ), ∇ • (w h + w )) = (f , (w h + w )) , ∀w ∈ V 0 (∇ • (v h + v ), (q h + q )) Ω = 0, ∀q ∈ Q. ( 18
)
To derive the stabilized formulation, we split Equations ( 18) into a large-scale and a fine-scale problem. The fine-scale problem is defined on element interiors. Under several assumptions about the time-dependency and the non-linearity of the momentum equation of the subscale system detailed in [START_REF] Hachem | Stabilized finite element method for incompressible flows with high reynolds number[END_REF], the fine-scale solutions v and p written in terms of the time-dependent largescale variables using residual-based terms that are derived consistently. Consequently, we can use static condensation, that consists in substituting directly v and p into the large-scale problem. This gives rise to additional terms in the Finite Element formulation, that are tuned by a local stabilizing parameter. These terms are responsible for the enhanced stability compared to the standard Galerkin formulation. The large-scale system finally reads:

             ρ (∂ t v h , w h ) Ω + (ρv h • ∇v h , w h ) Ω -K∈T h (τ 1 R M , ρv h ∇w h ) K + (2µε ε ε(v h ) : ε ε ε(w h )) Ω -(p h , ∇ • w h ) Ω -K∈T h (τ 2 R C , ∇ • w h ) K = (f , w h ) Ω , ∀w h ∈ V h,0 (∇ • v h , q h ) Ω -K∈T h (τ 1 R M , ∇q h ) K = 0, ∀q h ∈ Q h (19)
where (•, •) Ω represents the scalar product on the whole domain omega while (•, •) K is the scalar product on Element K. The quantities τ 1 and τ 2 are stabilization parameters defined hereafter. The momentum residual R M and the continuity residual R C are expressed as:

R M = f -ρ∂ t v h -ρv h • ∇v h -∇p h R C = -∇ • v h (20)
Compared to the standard Galerkin method, the proposed stable formulation involves additional integrals that are evaluated element-wise. These additional terms represent the stabilizing effect of the sub-grid scales and are introduced in a consistent way in the Galerkin formulation. They make it possible to avoid instabilities caused by both dominant convection terms and incompatible approximation spaces. All of these terms are controlled by the stabilization parameters τ 1 and τ 2 , for which we adopt the definition proposed in [START_REF] Codina | Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods[END_REF]:

τ 1 = 2ρ v h K h K 2 + 4µ h 2 K 2 -1 2 , ( 21
)
τ 2 = µ ρ 2 + c 2 v h K c 1 h K 2 1 2 (22)
where h K is the characteristic length of the element and c 1 and c 2 are algorithmic constants. We take them as c 1 = 4 and c 2 = 2 for linear elements [START_REF] Codina | Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods[END_REF]. Equations [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origin of stabilized methods[END_REF] are discretized in time by a semi-implicit scheme. The convective term, the viscous term and the pressure term in the momentum equation, as well the divergence term in the continuity equation, are integrated implicitly through a backward Euler scheme. All other contributions (i.e. the source term and the stabilization terms) are integrated explicitly by a forward Euler scheme. Due to the split between the Navier-Stokes and the Spalart-Allmaras equations, the dependence of the total viscosity µ on the turbulent viscosity µ t is also explicit. Among the implicit terms, the convection operator v h • ∇v h is non-linear, thus a Newton root finding algorithm is applied to the system at each time step. Using the superscript n to denote the value of the quantities at time t n and n + 1, i to denote the value of the quantities for time t n+1 at the i-th iteration of the Newton method, the system to be solved for (v n+1,i+1 h , p n+1,i+1 h ) at each iteration i reads:

                       ρ v n+1,i+1 h -v n h ∆t , w h Ω + ρ v n+1,i h • ∇v n+1,i+1 h + v n+1,i+1 h • ∇v n+1,i h -v n+1,i h • ∇v n+1,i h , w h Ω -K∈T h (τ n 1 R n M , ρv n h ∇w h ) K + 2µ n ε ε ε(v n+1,i+1 h ) : ε ε ε(w h ) Ω -p n+1,i+1 h , ∇ • w h Ω -K∈T h (τ n 2 R n C , ∇ • w h ) K = (f n , w h ) Ω , ∀w h ∈ V h,0 ∇ • v n+1,i+1 h , q h Ω -K∈T h (τ n 1 R n M , ∇q h ) K = 0, ∀q h ∈ Q h . (23)
We usually precondition the linear system with a block Jacobi method supplemented by a incomplete LU factorization, and solve it by means of a GMRES algorithm.

The Spalart-Allmaras turbulence model (SA)

The turbulence model chosen to compute the eddy viscosity is the one-equation Spalart-Allmaras (SA) turbulence model. The SA model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] represents the evolution of the kinematic eddy viscosity ν by a a non linear advection-diffusion-reaction equation:

∂ ν ∂t + v • ∇ν -c b1 (1 -f t2 ) S ν + c w1 f w - c b1 κ 2 f t2 ν d 2 - c b2 σ ∇ν • ∇ν - 1 σ ∇ • [(ν + ν)∇ν] = 0 (24)
The eddy viscosity can then be obtained from µ t = ρνf v1 , with:

f v1 = χ 3 χ 3 + c 3 v1 , χ = ν ν , f v2 = 1 - χ 1 + χf v1 f t2 = c t3 e -ct4χ 2 f w = g 1 + c 6 w3 g 6 + c 6 w3 1 6 , g = r + c w2 (r 6 -r), r = ν Sκ 2 d 2 , S = S + ν κ 2 d 2 f v2 , S = 2ε ε ε(v) : ε ε ε(v). ( 25
)
where d is the shortest distance to the wall, κ = 0.4 is the von Kármán constant, and the remaining model coefficients are:

c b1 = 0.1355, c b2 = 0.622, σ = 2/3, c v1 = 7.1, c v2 = 0.7, c v3 = 0.9, c w1 = c b1 κ + 1 + c b2 σ , c w2 = 0.3, c w3 = 2, c t3 = 1.2, c t = 0.5. (26)
In order to improve accuracy and convergence, some modifications may be performed, in particular to avoid the generation of negative viscosity values. There exist in the literature many variations for the SA model, most of which are collected in NASA's turbulence modeling resource webpage [START_REF] Rumsey | Nasa turbulence modeling resource page[END_REF]. In this work, the Negative Spalart Allmaras Model is used because of its capability to avoid the generation of negative turbulent viscosity without the use of clipping [START_REF] Allmaras | Modifications and clarifications for the implementation of the spalart-allmaras turbulence model[END_REF]. It consists in replacing [START_REF] Khurram | Multiscale finite element method applied to the spalart-allmaras turbulence model for 3d detached-eddy simulation[END_REF] when ν is negative by: [START_REF] Mesri | On optimal simplicial 3d meshes for minimizing the hessian-based errors[END_REF] with f n = (c n1 + χ 3 )/(c n1 -χ 3 ) and c n1 = 16. Moreover, the turbulent eddy viscosity µ t is set to zero when ν is negative.

∂ ν ∂t + v • ∇ν -c b1 (1 -c t3 )S ν -c w1 ν d 2 - c b2 σ ∇ν • ∇ν - 1 σ ∇ • [(ν + f n ν)∇ν] = 0,
A stabilized Finite Element discretization of the Spallart-Allmaras model is proposed in [START_REF] Khurram | Multiscale finite element method applied to the spalart-allmaras turbulence model for 3d detached-eddy simulation[END_REF]. Following a similar idea, we recast Equation ( 24) into a convection-diffusion-reaction form, and apply a backward Euler time discretization:

νn+1 -νn ∆t + v n+1 - c b2 σ ∇ν n+1 • ∇ν n+1 convection - 1 σ ∇ • (ν + νn+1 )∇ν n+1 diffusion -c b1 (1 -f n+1 t2 ) Sn+1 + c w1 f n+1 w - c b1 κ 2 f n+1 t2 νn+1 d 2 νn+1 reaction = 0, ( 28
)
where νn stands for the value of ν at discrete time t n . We recall that the value of the velocity field at time t n+1 is computed before the Spallart-Allmaras equation at the same time step, so the quantities v n+1 , Sn+1 , f n+1 t2 and f n+1 w are explicitly available when solving Equation [START_REF] Micheletti | Stabilized finite elements on anisotropic meshes: A priori error estimates for the advection-diffusion and the stokes problems[END_REF]. Equation ( 28) is non-linear. Following [START_REF] Khurram | Multiscale finite element method applied to the spalart-allmaras turbulence model for 3d detached-eddy simulation[END_REF], we apply the non-linear root finding algorithm at semi-discrete level, but we chose to work with a simpler Picard-like linearization instead of the Newton method employed in [START_REF] Khurram | Multiscale finite element method applied to the spalart-allmaras turbulence model for 3d detached-eddy simulation[END_REF]. Using the superscript νn+1,i for the value of νn+1 at the i-th iteration, an iteration of the non-linear root search reads:

νn+1,i+1 -νn ∆t + v n+1 - c b2 σ ∇ν n+1,i • ∇ν n+1,i+1 - 1 σ ∇ • (ν + νn+1,i )∇ν n+1,i+1 -c b1 (1 -f t2 ) Sn+1 + c w1 f w - c b1 κ 2 f t2 νi d 2 νn+1,i+1 = 0, (29)
Equation ( 29) is then discretized in space using a Streamline Upwind Petrov-Galerkin (SUPG) method. The Galerkin formulation is obtained by multiplying (29) by appropriate test functions ω h ∈ W h , applying the divergence theorem to the diffusion terms and integrating by parts. Following the lines in [START_REF] Badia | Analysis of a stabilized finite element approximation of the transient convection-diffusion equation using an ale framework[END_REF] on the use of stabilization methods for transient convection-diffusion-reaction equations, the stabilized weak form of equation ( 29) reads:

νn+1,i+1 -νn ∆t , ω h Ω + v n+1 - c b2 σ ∇ν n+1,i • ∇ν n+1,i+1 h , ω h Ω - 1 σ (ν + νi h )∇ν n+1,i+1 h , ∇ω h Ω - c b1 (1 -f n+1 t2 ) Sn+1 h + c w1 f n+1 w - c b1 κ 2 f n+1 t2 νn+1,i h d 2 νn+1,i+1 h , ω h Ω + K R(ν n+1,i ), τ n+1,i 3 v n+1 h - c b2 σ ∇ν n+1,i h • ∇ω h K = 0, ∀ω h ∈ W h . ( 30
)
This article is protected by copyright. All rights reserved.

where R(ν) is the finite element residual of [START_REF]On the performance of high aspect ratio elements for incompressible flows[END_REF]. The stabilization parameter τ 3 , is computed within each element as:

τ 3 = c 2 h α c K + c 1 h 2 K α d + α r -1 (31)
where α c , α d and α r are respectively the convection, diffusion and reaction coefficients in Equation ( 30), h K is the element size, α c K a characteristic norm of the convection term and c 1 = 4, c 2 = 2 for linear elements.

The linear system arising from Equation ( 30) is solved using the same numerical method as for the Navier-Stokes equations (see Section 3.1).

Element size measures in stabilization parameters

Recall that the coefficients (τ 1,2,3 ) weight the stabilization terms added to the weak formulations ( 19) and [START_REF] Parolini | Mathematical models and numerical simulations for the americas cup[END_REF]. They are defined for each element K of the triangulation and depend on the local mesh size h K . Many numerical experiments show that good results can be obtained when using the minimum edge length of K [START_REF]On the performance of high aspect ratio elements for incompressible flows[END_REF], while others always use the triangle diameter (see [START_REF] Micheletti | Stabilized finite elements on anisotropic meshes: A priori error estimates for the advection-diffusion and the stokes problems[END_REF] for details).

However, in the case of strongly anisotropic meshes with highly stretched elements, the definition of h K is still an open problem and plays a critical role in the design of the stabilizing coefficients [START_REF] Harari | What are c and h?: inequalities for the analysis and design of finite element methods[END_REF][START_REF] Codina | Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods[END_REF]. In [START_REF] Förster | Stabilized finite element formulation for incompressible flow on distorted meshes[END_REF], the authors examine deeply the effect of different element length definitions on distorted meshes. In [START_REF] Cangiani | The residual-free-bubble finite element method on anisotropic partitions[END_REF], anisotropic error estimates for the residual free bubble (RFB) method are developed to derive a new choice of the stabilizing parameters suitable for anisotropic partitions. In this work, we adopt the definition proposed in [START_REF] Tezduyar | Finite element stabilization parameters computed from element matrices and vectors[END_REF] to compute h K as the size of K in the direction of the velocity:

h K = 2|v h | Σ N K i=1 |v h • ∇ϕ i | ( 32
)
where N K is the number of vertices of K and ϕ 1 , ..., ϕ N K are the usual basis functions of P 1 (K) mapped onto K.

NUMERICAL EXPERIMENTS

This section describes the numerical tests conducted to assess the correctness and the efficiency of the method. We compare results obtained for steady and unsteady 2D classical benchmarks such as the flow around a prismatic cylinder and the NACA 0012 airfoil with data from the literature. For all test cases, we use different type of meshes and number of elements to make adequate comparisons. Finally, we present a more complex application, namely the flow past the 3D geometry of a UAV.

Evaluation of aerodynamic forces

The aerodynamic forces, which are the quantities of prime engineering interest in an aerodynamic simulation, are most often expressed in terms of the lift coefficient C L and the drag coefficient C D defined as

C D = 2F D ρ Ū 2 A C L = 2F L ρ Ū 2 A ( 33
)
where A is a reference area, F D is the drag force oriented in the direction e x parallel to the flow and F L is the lift force oriented in a direction e y normal to the flow. In 2D, these coefficients are usually expressed per unit of span length by substituting a reference length L for the reference area A.

The total aerodynamic force exerted by the flow on an object of boundary S ∈ ∂Ω is given by

F aero = S σn so F D = F aero • e x = S (σn) • e x F L = F aero • e y = S (σn) • e y (34)
We use the variational technique described in [START_REF] Volker | Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder[END_REF] to compute the components F D and F L of F aero . Starting with the drag force, let us consider the Sobolev space [H 1 ∂Ω\S (Ω)] D of functions supported in Ω that are equal to e x on S and that vanish on all other boundaries of Ω. A variational form of the momentum conservation Equation ( 15) can be obtained by multiplying by a test function and integrating by parts:

Ω ∂ t v • w + Ω σ : ∇w + ∂Ω (σn) • w + Ω (v • ∇) v • w = 0, ∀w ∈ [H 1 ∂Ω\S (Ω)] D ( 35
)
where n is the normal pointing inside the computational domain. Considering that w| ∂Ω\S = 0, the third term in [START_REF] Shimada | Application of a modified k -model to the prediction of aerodynamic characteristics of rectangular cross-section cylinders[END_REF] reduces to F D , and rearranging the terms yields

F D = - Ω ∂ t v • w - Ω σ : ∇w - Ω (v • ∇) v • w, ∀w ∈ [H 1 ∂Ω\S (Ω)] D . ( 36
)
In practice, we use for w the piecewise linear function that is equal to e x at the nodes lying on S and 0 at every other nodes, i.e. the function that vanishes everywhere except on the first layer of elements around S. The lift force can be obtained in the same manner, substituting e y for e x .

Flow past a 2D cylinder

As a first numerical test, we disregard the turbulence model (µ t = 0) and consider the classical benchmark proposed in [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF], that is, the laminar flow past a cylinder located slightly asymmetrically in a channel. The main point of this test is to demonstrate the ability of the mesh adaptation framework to capture local flow features. Also, it assesses the capacity of the method described The detailed geometry, parameters, flow conditions and boundary conditions are described in [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF].

The first set of simulations corresponds to the steady case with Reynolds number Re = 20. Two tests are performed, starting with a uniform solution and marching in time until the steady state is reached. In the first test, a fixed fine mesh is used. In the second one, dynamic mesh adaptation is triggered every a given frequency (i.e. here we fixed it at 10 iterations) under the constraint of a fixed number of nodes N = 40, 000, starting from an arbitrary mesh.

Recall that a multi-criteria adaptivity is used in this case. As explained in Section 2.4, each component of the velocity field and its norm were used to account for several sources of error.

The final drag and lift coefficients are listed in Table I, where they are compared with the reference data from [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF]. Very good agreement is found (about 0.2% in drag and 8% in lift).

To highlight the capability of the proposed adapted method, we repeated the same case with a Reynolds number equal to 100. The same dynamic mesh adaptation parameters as for the steady case are used. The velocity norm and the corresponding mesh at several points in time are shown in the upper part of Fig. 2. As expected, the boundary layers on the cylinder and on the channel walls, as well as the flow detachments, are captured automatically. The maximum drag and lift (i.e. 3.21 and 0.94 respectively) are again in good agreement with the reference.

Turbulent flow around prismatic cylinders

We turn now to a case with fully turbulent flow. Non-smooth cylindrical geometries are often employed to assess stabilized finite element methods, because the sharp and localized flow features that they develop challenge the stabilization techniques. For the same reason, they represent a good test for mesh adaptation procedures. In this Section, we focus on the test cases proposed in [START_REF] Bao | Numerical prediction of aerodynamic characteristics of prismatic cylinder by finite element method with spalart-allmaras turbulence model[END_REF]. They involve geometries of 2D primatic cylinders that are parametrized by lengths H 1 and H 2 as shown in Figure 3. The computational domain is 33H 1 long and 21H 1 wide. The cylinder is centered in width, and its front side is located 8H 1 downstream of the inlet boundary. Three configurations (a), (b) and (c), corresponding to aspect ratios H 2 /H 1 of 0, 0.6 and 1 respectively, are considered. The values of the inlet velocity V in and the dynamic viscosity µ are chosen to yield a Reynolds number based on H 1 equal to Re = 2.2 × 10 4 . The inflow boundary conditions are v = (V in , 0), ν = 3ν, as recommended in [START_REF] Philippe | Effective inflow conditions for turbulence models in aerodynamic calculations[END_REF]. Symmetry conditions are imposed on the outer boundaries parallel to the flow direction as v y = 0, ∂v x /∂y = 0 and ∂ ν/∂y = 0. For the outflow, p • n -(µ + µ t ) ∇v + ∇ T v • n = 0 and ∇ν • n = 0 are prescribed, whereas for the rigid body no slip conditions are imposed as v = 0 and ν = 0. The nondimensional time step is chosen equal to ∆t = 0.002. We are interested in assessing the evolution of the drag and lift coefficients over time, once a condition of periodic shedding of vortices has been reached. These quantities are computed by means of the method described in Section 4.1.

The plots in the upper part of Figure 4 depict the turbulent variable ν at a time in the periodic regime of the flow for the three different geometries with fixed meshes. It is clear that the turbulence model is only active in the restricted area of interest, and the effects of the averaging process and damping function seem to be correctly taken into account. The streamlines, shown in the bottom part of Figure 4, are similar to those given in [START_REF] Bao | Numerical prediction of aerodynamic characteristics of prismatic cylinder by finite element method with spalart-allmaras turbulence model[END_REF] in the cylinder wake. The results in Table II are in good agreement with the literature presented in [START_REF] Bao | Numerical prediction of aerodynamic characteristics of prismatic cylinder by finite element method with spalart-allmaras turbulence model[END_REF] and the references therein. Some differences may however be spotted (for instance in the lift coefficient C L with Configuration (c)). This may be due to the fineness of the spatial discretization, but also to the time discretization, as pointed out in [START_REF] Volker | Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder[END_REF]. In addition to the use of a first-order implicit time integration scheme, the time splitting error may contribute to this inaccuracy.

As in the previous case, the effects of mesh adaptation are investigated on Configuration (b). The adaptive process always starts with an arbitrary uniform mesh and is iterated every 5 time steps. It considers again a multi-criteria adaptivity taking into account the velocity components, its norm and the turbulent viscosity. The number of nodes N is set to 10 5 . As shown in Figure 5, the dynamic mesh adaptation capture well the evolution of the solution. Note the concentration of the resolution not only along the boundary layers but also in the wake regions. This reflects well the anisotropy of the solution caused by the nature of the flow past the cylinder. The snapshots in Figure 5 explains also, how for a controlled number of nodes, the mesh is refined around the newly appeared vortices and automatically is coarsened in the regions far from the cylinder.

The zoom-in in Figure 6 highlights how sharply the layers are captured. It shows that the shape but also the size and the orientation of the elements matches the directional features of the flow (boundary layers, flow detachments). The zoom-out in Figure 6 shows how the anisotropic adaptive procedure modifies the mesh in a way that the local mesh resolutions become adequate in all directions. The mean aspect ratio for the anisotropic elements is O(1:1000). Again, the developed stabilized finite element methods show to be very efficient and robust for solving flows at high Reynolds number using highly stretched elements.

The evolution of drag and lift coefficients for different aspect ratios H 2 /H 1 is shown in Figure 7. Both reach as expected a steady oscillating state, and compare well to the plots shown in [START_REF] Bao | Numerical prediction of aerodynamic characteristics of prismatic cylinder by finite element method with spalart-allmaras turbulence model[END_REF].

NACA0012 airfoil at

Re = 6 × 10 6
The objective of this test case is to demonstrate the capability of the method to simulate turbulent flows at Reynolds numbers typical of aeronautical applications, as well as to assess the multi-criteria adaptation feature presented in Section 2. We consider the flow around a NACA0012 airfoil [START_REF] Rumsey | Nasa turbulence modeling resource page[END_REF], that is a standard benchmark for aerodynamic simulations. The geometry of the NACA airfoil is placed in a domain that is 8 chord lengths wide and 20 chord lengths long, which is considered sufficiently large to avoid any influence of the boundaries on the flow in the vicinity of the airfoil (see Figure 8). The Reynolds number based on the unitary cord length is Re = 6 × 10 6 and the Mach number is M a = 0.2, thus ensuring the validity of the incompressibility hypothesis.

Flat velocity and turbulent viscosity (ν = 3ν) profiles are imposed at the inlet boundary. Slip conditions are prescribed on the boundaries parallel to the flow as v • n = 0. On the profile, noslip boundary conditions for the velocity as well as homogeneous Dirichlet conditions for the SA working variable ν are imposed. Time marching is performed until the steady-state is reached. The non-dimensional time step is set to ∆t = 10 -4 .

As for the preceding test cases, the effects of the mesh adaptation process described in Section 2 are investigated, in particular with regard to the variables driving the adaptation and the use of a fixed number of elements. Mesh adaptation is triggered every 5 time steps.

In a first phase, two cases are studied with a constant number of N = 100, 000 nodes in the mesh. In the first one, only the velocity (norm and directions) is included in the adaptation criterion, whereas in the second case, the combination of the velocity and the turbulent viscosity is considered. The meshes obtained are depicted in Figure 9. It is clear that adapting on both the velocity and the turbulent viscosity increases mesh resolution in the wake. In a second phase, the adaptation criterion including both velocity and turbulent viscosity is chosen, and three tests are carried out, with a constraint of N = 50, 000, N = 100, 000 and N = 150, 000 nodes respectively. The meshes obtained at the stationary state are shown in Figure 10. We notice that pressure gradients are too localized to drive the adaptation procedure in a way that would be beneficial to the resolution of the complete flow. As expected, a high degree of mesh anisotropy is reached in all the cases because of the steep local variations of the velocity and the turbulent viscosity. Note the high resolution in the wake: the element size increases anisotropically when moving away from the trailing edge.

The effect of the mesh resolution and adaptation on the drag coefficient is summarized and compared to results obtained by different codes [START_REF] Rumsey | Nasa turbulence modeling resource page[END_REF] in Table III. The influence of the mesh size is clearly visible. Under the constraint of a fixed number of elements, the anisotropic mesh adaptation algorithm, based on the principle of equi-distribution of the error, increases the mesh size at the interface to better capture and refine the mesh in the wake. As may be seen from Table III, all tests except the first one are in a 1% range of error, which confirms that it is appropriate to solve the RANS equations by combining stabilized Finite Element methods and mesh adaptation. The influence of the mesh adaptivity criteria is also well highlighted. Indeed, using only the velocity field as a criteria shows higher error for the drag coefficient. Finally, the solutions and meshes obtained in Test 2 are shown in Figure 11 and seem to be in agreement with similar tests in the literature. The mean aspect ratio for the anisotropic elements is O(1:1000).

Analysis of turbulent flow past a drone

In the last numerical test case, the aim is to assess the capability of our numerical method to simulate the turbulent flow past a complex 3D geometry, namely a UAV. We seek to reproduce experimental results obtained for several angles of attack (0, 5, 10, 15) with a set-up pictured in Figure 12. The The model is placed in a computational domain of length 3m and width 1m. Similarly to the preceding test cases, flat velocity and turbulent viscosity (ν = 3ν) profiles are prescribed at the inlet boundary, a pressure condition is imposed at the outlet, and slip conditions are used on the boundaries parallel to the flow. At the UAV wall, no-slip boundary conditions for the velocity as well as homogeneous Dirichlet conditions for the turbulent viscosity ν are specified. Again, time marching is performed until the steady-state is reached. The non-dimensional time step is set to ∆t = 10 -4 . Figure 13 gives a general view of the computational domain.

All the 3D simulations have been conducted in parallel on 64, 2.4Ghz Opteron cores and have required 12.5 hours of computational time. On average, 4 nonlinear iterations were required to reach convergence.

In the initial configuration, the UAV is simulated at zero angle of attack. The mesh is adapted every 10 time steps while marching in time. When steady state is reached, the configuration is instantaneously rotated to the following prescribed angle of attack. Figure 14 shows the steady-state mesh for zero angle of attack in the vicinity of the UAV. Results obtained for the drag are presented in Figure 15. The agreement between the numerical and the experimental results is satisfying despite the discrepancies between the experimental setup and the numerical simulations. Note that three experiments were made for each case and the averaged value is considered here. In particular, the surface roughness of the model obtained by additive manufacturing may be the cause of the higher drag, due to the increased turbulent skin friction. Figure 16 shows a visualization of the flow past the UAV by means of streamlines and pressure distribution.

Mesh Turbulent viscosity

CONCLUSION

An adaptive anisotropic mesh methodology for performing accurate numerical simulations of turbulent flows past complex geometries was presented. It couples a stabilized Variational Multiscale Navier-Stokes solver, a Spalart-Allmaras turbulent model with a dynamic anisotropic mesh adaptation algorithm. The benefit from using an edge-based error estimator with a multicriteria functional that estimates the error on a combination of several variables including the velocity norm, all the velocity components as well as the turbulent viscosity is shown on different challenging cases in 2 and 3 dimensions. All the main directional features, characterizing the turbulent flows past immersed solids, are well detected and captured. The comparison with the literature show that despite a low number of nodes, simulations do not suffer any instabilities and a good accuracy is obtained. Further investigations will take into account the performance of the presented adaptive method in terms of the computational cost, the speed-up and the frequency of remeshing in a massively parallel environment.

Figure 1 .

 1 Figure 1. Patch associated with node x i .

Figure 2 .

 2 Figure 2. Dynamic anisotropic mesh adaptation and fluid flow past a fixed cylinder.

Figure 3 .

 3 Figure 3. Geometry for cylinder case in section 4.3.

Figure 4 .

 4 Figure 4. Plots of ν (top) and streamlines (bottom) for Configurations (a) (left), (b) (center) and (c) (right).

Figure 5 .

 5 Figure 5. Subsequently adapted meshes and the corresponding velocity norm.

Figure 6 .

 6 Figure 6. Zoom-out with vortex shedding (left) anf zoom-in the adapted mesh (right).

Figure 7 . 2 Figure 8 .

 728 Figure 7. Drag (lines) and lift (dots) coefficients for aspect ratios (a) H 2 /H 1 =0.0 ,(b) H 2 /H 1 =0.6 and (c) H 2 /H 1 =1.0 .

Figure 9 .

 9 Figure 9. Adaptation driven by the velocity (left) and combination of velocity and turbulent viscosity (right).

Figure 10 .

 10 Figure 10. The obtained meshes using (a)50, 000, (b)100, 000 and (c) 150, 000 elements.

Figure 11 .

 11 Figure 11. Final plots for Test # 2.

Figure 12 .

 12 Figure 12. The experimental setup: the drone inside a wind tunnel

Figure 13 .Figure 16 .

 1316 Figure 13. The immersed drone inside the numerical wind channel

Table I .

 I Drag and lift coefficient values with literature for fixed mesh (Test 1) and adaptive mesh (Test 2). Section 4.1 to correctly evaluate the aerodynamic coefficients C D and C L on highly anisotropic meshes.

		C D	C L
	Schäfer and Turek [34] 5.58 ± 0.01 0.0107 ± 0.0003
	Test 1	5.57	0.0098
	Test 2	5.56	0.0099
	in

Table II .

 II Comparisons of drag and lift coefficients with the literature for Configurations (b) and (c).

		mean C D	r.m.s. C L
	Bearman and Obasaju (c) [5]	2.10	1.20
	Sohankar et al. (c) [36]	2.25	1.50
	Shimada and Ishihara (c) [35]	2.05	1.43
	Bao et al. (c) [4]	2.04	1.24
	Present work (c)	2.08	1.57
	Bao et al. (b) [4]	2.50 ± 0.125 1.7 ± 0.1
	Present work (b)	2.57	1.79

Table III .

 III Summary of test cases. span of 300 mm and a total length of 300 mm. The inlet velocity is set to 10 m/s, yielding a span-based Reynolds number of Re = 2 × 10 6 .

	Reference	Test #1 Test #2 Test #3 Test #4
	No. elements ×10 5	0.5	1.0	1.5	1.5
	Adaptation using
	35	7.90	8.12	8.16	7.7
	Relative error	3%	0.3%	0.1%	4.5%
	model has a

v & ν t v & ν t v & ν t v C D × 10 -3

8.155 ± 0.