Eléonore Gauci
email: eleonore.gauci@inria.fr

Anca Belme
email: belme@dalembert.upmc.fr

Alexandre Carabias
email: alexandre.carabias@lemma-ing.com

Adrien Loseille
email: adrien.loseille@inria.fr

Frédéric Alauzet
email: frederic.alauzet@inria.fr

Alain Dervieux
email: alain.dervieux@inria.fr

A priori error-based mesh adaptation in CFD

Introduction

The main purpose of Computational Mechanics research is not to obtain more efficient computational methods on existing cases (or more accurate, which is somewhat equivalent). This purpose is to perform new computations which were not affordable before. Future researchers will hopefully be able to face the many challenges of the century (global change, health, energy). This remark holds (to a lesser extend) for mesh adaptation research. The purpose is to be able to compute new phenomena and also to master the numerical uncertainties which have been up to now unsufficiently controlled.

Many phenomena in continuum mechanics real world involve multiple scales or/and singularities and cannot computed with uniform meshes. Designing non uniform meshes will become so complex that automatic adaptation will be mandatory.

Second-order accurate anisotropic mesh adaptation has demonstrated that phenomena which were impossible to compute with non-adaptive meshes become computable even on medium size computers. The first section of this paper shows an example. But adaptative or not second-order methods are limited to routine computation for a low accuracy, typically with a 1% relative error. In many cases this is sufficient for applications and coherent with the accuracy of available models (e.g. basic statistical turbulence modeling). When a smaller error tolerance (0.1%, 0.01%) is necessary, second-order accurate methods become quite inefficient. Also, second-order accurate methods are definitively inefficient for many unsteady phenomena (e.g. wave propagation or with instability) The error accumulates while time is advanced. Due to multiple scales, explicit time-advancing can be prohibitive. Implicit time-advancing deteriotates the fastest waves. Then higher-order adaptive methods must be developed in order to efficiently adress the accuracy question.

To sum up this first discussion, there is a class of computations which are not possible without an increase of efficiency. Efficiency generally results from accuracy, algorithmics, and computers. We concentrate here on the increase of accuracy of approximation through mesh adaptation.

Since we need to improve the efficiency of mesh adaptation, we need to have tools to measure and predict it, replacing human intuition and heuristics by mathematical concepts. The existing mathematical roadmap for error control is the result of many scientific contributions. A mathematical definition of the mesh adaption purpose was proposed by the introduction of goal-oriented formulations [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF] dealing with the accuracy in the evaluation of a prescribed scalar output. But reducing the error of solely a scalar output is not always sufficient for the engineer who wants to increase its understanding by looking at the flow computed, or wants to use it for building reduced order models. Thanks to norm-oriented formulations [START_REF] Brèthe | Anisotropic Norm-Oriented Mesh Adaptation for a Poisson problem[END_REF][START_REF] Loseille | Anisotropic Norm-Oriented Mesh Adaptation for Compressible Flows[END_REF], error control applies to solution fields and Banach norms. At the other end of the process, namely the mathematical definition of the control parameter, the replacement of the mesh by a Riemannian metric allowed for transforming the mesh adaptation issue into a well-posed optimization problem. Finally we have the optimal mesh (according to our preferences) and, at least for the norm-oriented option described in the sequel, we have an idea of the numerical error. This paper presents a review of recents theoretical advances together with examples of applications. In particular the mesh adaptation roadmap which we define extends to higher-order accurate approximation and to unsteady phenomenas. This review is accompanied by new ensights. Among others, we revisit the application of a mathematical analysis predicting the asymptotic accuracy of the higher-order adaptation algorithm when the solution field contains singularities.

Features, goal, norms -oriented formulations

The three formulations are presented for the scond-order accurate particular case of P 1 approximations on triangles and tetrahedra.

Feature-based (FB) adaptation

The feature-based anisotropic mesh adaptation criterion is inspired by the classical finite-element a priori estimate. It assumes that a good mesh for the interpolation of the unknown (or of one or several sensors depending of it) is a good mesh for a Partial Differential Equation (PDE) approximation based on the same interpolation. This method is generally defined in a continuous framework by replacing the local interpolation error by the application of the Hessian of the exact solution times a local mesh size defined by a continuous metric, see [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF][START_REF] Loseille | Adaptive anisotropic simulations in aerodynamics[END_REF][START_REF] Castro-Díaz | Anisotropic unstructured mesh adaptation for flow simulations[END_REF][START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent CFD[END_REF][START_REF] Vasilevski | Error bounds for controllable adaptive algorithms based on a hessian recovery[END_REF][START_REF] Huang | Metric tensors for anisotropic mesh generation[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF]. It can be also defined on a discrete basis as in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF].

Given a numerical solution W h , a solution of higher regularity R h (W h ) is recovered, so that the following interpolation error estimate [START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in L p -norm[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF] holds:

R h (W h ) -Π h R h (W h ) L p ≤ N -2 3 Ω det |H R h (W h ) | p 2p+3 2p+3 3p
, where H R h (W h ) is the Hessian of the recovered solution and N an estimate of the desired number of nodes. Notation Π h holds for the usual elementwise affine interpolation relying on values at vertices. When working in the continuous mesh framework [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF], we also have the optimal metric field leading to the previous error estimate:

M L p = K L 1 ( H R h (W h ) ) = D L p (det H R h (W h ) ) -1 2p+3 H R h (W h ) .
If anisotropic mesh prescription is naturally deduced in this context, interpolationbased methods do not take into account the features of the PDE. However, in some simplified context and assumptions (elliptic PDE, specific recovery operator R h ), we have:

W -W h ≤ 1 1 -α R h (W h ) -Π h R h (W h ) with α > 1 ,
so that a good convergence to the exact solution may be observed [START_REF] Loseille | Achievement of global secondorder mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF]. Indeed, if R h (W h ) is a better approximate of W in the following meaning:

W -W h ≤ 1 1 -α R h (W h ) -W h where 0 ≤ α < 1,
and if the reconstruction operator R h has the property:

Π h R h (W h ) = W h ,
then we can bound the approximation error of the solution by the interpolation error of the reconstructed function R h (W h ):

W -W h ≤ 1 1 -α R h (W h ) -Π h R h (W h ) .
Note that from a practical point of view, R h (W h ) is never recovered, only its first and second derivatives are estimated. Standard recovery techniques include leastsquare, L 2 -projection, Green formula or the Zienkiewicz-Zhu recovery operator. For all the numerical examples, the L 2 norm is used with the L 2 -projection to recover derivatives. These estimates are very efficient to capture anisotropic physical phenomenon. A typical example is the prediction of the sonic boom signature of a supersonic aircraft. Indeed, in the case, every singularities of the geometry generates a shock waves that propagates down to the ground. All physical variables (density, pressure, . . . ) are discontinuous accross the shock. Special features dedicated to sonic boom mesh adaptation are discussed in [START_REF] Loseille | Adaptive anisotropic simulations in aerodynamics[END_REF]. To illustrate the use of FB mesh adaptation, we consider the C25D geometry of the 2nd Sonic-Boom Prediction Workshop [START_REF]2nd AIAA Sonic Boom Prediction Workshop[END_REF] represented in Figure 1. We use the Mach number as sensor. The sequence of generated meshes is reported in Table 1. Several cuts in the volume mesh and the imprint of the adapted mesh in the symmetry plane are depicted in Figure 2. The feature-based approach is particularly attractive due to its simplicity and to the interest in taking into account physical aspects through the choice of the sensors. However, for systems, the choice of sensors is extremely sensitive. Some additional mathematics for defining the best choice, through the goal-oriented approach, are required for most complex applications.

The goal-oriented formulation

The goal-oriented (GO) mesh adaptation, as stated before, focuses on deriving the optimal mesh for computing a prescribed scalar quantity of interest (QoI). This method is a priori different from the optimal mesh for a field/state. It is particulary suitable for a large class of problems in computational mechanics where one wants to better understand/compute such a quantity of interest (e.g. drag, lift coef-ficients, vorticity in a wake, a mass flow across a surface, etc,) while maintaining a reasonable cost. Such GO methods require however the computation of an adjoint state. Adjoint-based methods were known to be particulary useful in optimal control/design, where the computational cost for evaluating sensitivities can be importantly reduced using the adjoint-state. It has been also proved for a few years that using an adjoint is efficient for improving by a correction the accuracy of the evaluation of a quantity of interest ( [START_REF] Pierce | Adjoint Recovery of Superconvergent Functionals from PDE Approximations[END_REF]). As GO mesh adaptation is concerned, a remarkable quantity of papers deals with a posteriori goal-based error formulation to drive adaptivity, using adjoint formulations or gradients (see for example [START_REF] Ferro | Anisotropic mesh adaptation for crack propagation induced by a thermal shock in 2D. Accepted for the publication[END_REF][START_REF] Yano | An optimization framework for anisotropic simplex mesh adaptation: application to aerodynamics flows[END_REF]). Our group is investigating GO formulations for steady and unsteady problems in an a priori context. Loseille et al. [START_REF] Loseille | Anisotropic Norm-Oriented Mesh Adaptation for Compressible Flows[END_REF] derived the goal-based error estimate in a steady context for Euler flows, showing that the QoI error estimate is expressed as a weighted interpolation error on solution flow fields. This leads to an optimal metric computed as a sum of Hessians of Euler flux fields weighted by gradient components of the adjoint state.

A short description of our goal-based a priori error estimation for the compressible Navier-Stokes system is now presented.

For simplicity, the QoI j is assumed to be the scalar product of the state solution W by a given field g, j(w) = (g,W ) where W is the exact state solution to the continuous model Ψ (W ) = 0. We denote by W h the solution of a discretisation, Ψ h (W h ) = 0. Let W * be the adjoint state, satisfying (∂Ψ /∂W )W * = g. Under adhoc regularity assumptions concerning Ψ ,Ψ h ,W and W * , we have ( [START_REF] Loseille | Anisotropic Norm-Oriented Mesh Adaptation for Compressible Flows[END_REF]) the following compact form for our a priori error estimate:

|(Ψ (W ) -Ψ h (Π h W ), Π h W * )| Ω ∑ k G k (W, ∇W * , | ρ H (W * )|) S k (W ) -Π h S k (W ) dΩ , (1)
where A B holds for a majoration asymptotically valid, i.e. A ≤ B + o(A) when mesh size tends to zero. Expression | ρ H (ϕ) | holds for spectral radius of the Hessian of ϕ. The interest of such an estimate is that the error model is a sum of interpolation errors weighted by algebraic functions of state and adjoint. Let us consider the example of the Navier-Stokes system:

Ψ (W ) = 0 on Ω , with Ψ (W ) = ∇ • F E (W ) + ∇ • F V (W ) (2)
where W = t (ρ, ρu, ρE) is the conservative flow variables vector and vector F E represents the Euler fluxes, and F V the viscous fluxes defined hereafter:

F V (W ) = [0, σ , -(q -u • σ )] T ; σ = µ(∇u + ∇u T ) - 2 3 µ∇ • uI,
with µ representing the constant viscosity. The heat flux q is given by Fourier's law q = -λ ∇T , where λ is the heat conduction (assumed here to be constant). We have noted ρ the density, u = (u 1 , u 2 , u 3 ) the velocity vector, H = E + p/ρ is the total enthalpy, E = T + u 2 2 the total energy and p = (γ -1)ρT the pressure with γ = 1.4 the ratio of specific heat capacities, T the temperature, and (e 1 , e 2 , e 3 ) the canonical basis. Expliciting Formulation (1) for the Navier-Stokes equations leads to the following developement:

|(Ψ (W ) -Ψ h (Π h W ), Π h W * )| E
with:

E = Ω | ∇W * | • F E (W ) -Π h F E (W ) dΩ + d + 1 3 K d d ∑ i=1 Ω µ | ρ H (W * ρu i )| |u i -Π h u i | dΩ + 1 3 K d d ∑ i=1 d ∑ j=1 j =i Ω µ | ρ H (W * ρu i )| |u j -Π h u j | dΩ + K d Ω λ | ρ H (W * ρE )| |T -Π h T | dΩ + (d + 1 3 ) K d d ∑ i=1 Ω µ | ρ H (W * ρE )| |u i | |u i -Π h u i | dΩ + 1 3 K d d ∑ i=1 d ∑ j=1 j =i Ω µ | ρ H (W * ρE )| |u i | |u j -Π h u j | dΩ + 5 3 d ∑ i=1 d ∑ j=1 j =i Ω µ ∂ (Π h W * ρE ) ∂ x j (u j -Π h u j ) ∂ u i ∂ x i -(u i -Π h u i ) ∂ u j ∂ x i dΩ ,
in which K d = 3 in two dimensions, K d = 6 in three dimensions. Taking into account all the terms in E is a sufficient condition for an efficient GO adaptation. These terms are computed only once for each mesh considered in the adaptation loop.

As illustration, we consider the 3D computation of a turbulent shock-boundary layer interaction problem. A similar problem was addressed with FB mesh adaptation in [START_REF] Loseille | Boundary layer mesh generation and adaptivity[END_REF]. The shock is initiated by a diamond shape obstacle in side a rectangular channel. The zone of interest is the region where the shock hits the bottom of the channel. The incoming flow field is defined such that: M ∞ = 1.4, the angle of attack α = 0 and Re = 3.6 × 10 6 . The QoI is the drag on the bottom flate plate. A mesh adaptation loop has been applied for a total number of vertices of 1.5 million.

The complex shock waves pattern is accurately captured by the adaptive process in all the regions where it impacts the flat plate. We observe that at the end of the domain, the mesh is rather coarse. In contrast, the boundary layer is highly refined. A close-up view of the shock boundary-layer interaction region shows that we cap- ture accurately the separation bubble and the lift off of the boundary layer due to the impact of the oblique shock. We also capture the reflected shock in front of the bubble and the expansion fan at the foot of the oblique shock. All these waves are visible in the mesh.

Remark 1: During the presentation of this work at the conference, O. Pironneau raised some paradox in the proposed goal-oriented adaptation algorithm. The question was: how can you adequately adapt with a criterion expressed in terms of the state and adjoint state if you do not adapt for a better state and a better adjoint?

This question is related to a particular step in the design of the mesh adaptation algorithm. Indeed, the first step in our theory is to apply a continuous metric analysis, which produces a continuous optimality system involving continuous state, continuous adjoint and continuous stationarity with respect to the metric.

In a second step, we have to discretize the whole optimality system. For this, an option could be to use a sequence of uniformly refined meshes. State and adjoint would be well refined. This option would finally produce (with a prohibitive computational cost) the optimal metric.

Once we have a good approximation of the optimal metric, we can generate unit adapted meshes based on it. By construction, these meshes are adequate for the evaluation of the functional, through the evaluation of the state. In other words the state computed with an unit mesh of the metric is a good mesh for the functional evaluation. At the same time, with this adapted mesh, the state is not well computed since only its features influencing the functional are supposed to be well computed. Adapting for a better state would then be not optimal.

Adapting for a better adjoint is a somewhat second-order approach, since it does not influence directly the quality of the approximation of the functional, but solely how close is the metric from the perfectly optimal one : formally, an ε-large deviation of the (unconstrained) quasi-optimum provokes an ε 2 -large deviation of the (smooth) functional value. A similar argument applies to an adaptation for the stationarity equation.

Remark 2:

The GO method (as explained here) shows the following limitations: -It allows to minimize the error on a functional, but does not give an idea of the effective error on this functional, but the corrector of [START_REF] Pierce | Adjoint Recovery of Superconvergent Functionals from PDE Approximations[END_REF] can be used as an approximation of it.

-In many cases the engineer needs a good evaluation of the fields themselves. The norm-oriented technology described in next section addresses these two issues.

Numerical corrector

Given a discrete problem, a mesh (of metric) M h and the discrete solution W h computed with the mesh, we call "numerical corrector" a discrete field W h such that the sum W h + W h is a significantly more accurate approximation of the exact solution than W h . In other words, W h is a good approximation of the error. Clearly, W h is useful for estimating the error, for correcting it, and by the way for building the norm-oriented mesh adaptation algorithm defined in Section 2.4.

A trivial way to compute W h could be to first compute an extremely accurate and extremely costly W h/2 k (k large) solution computed on a mesh M h/2 k obtained by dividing k times the elements of M h , and finally to put W h = W h/k -W h . But the interesting feature of a numerical corrector should be that its computational cost is not much higher than the computational cost of W h . We describe now a corrector evaluation of low computational cost relying on the application of a Defect Correction principle and working on the initial mesh M h :

Ψ h (W h + W h,DC ) ≈ -R h/2→h Ψ h/2 (R h→h/2 W h ) ; W h,DC = 4 3 W h,DC -(π h W h -W h )
where the π h W h -W h is a recovery-based evaluation of the interpolation error (see [START_REF] Loseille | Anisotropic Norm-Oriented Mesh Adaptation for Compressible Flows[END_REF] for details). The notation R h/2→h holds for the transfer (extension by linear interpolation) operator from the twice finer mesh M h/2 to the initial mesh M h . The notation R h→h/2 holds for the transfer (restriction) operator from the the initial mesh M h to the twice finer mesh M h/2 . The finer-mesh residual Ψ h/2 (W h ) can be assembled by defining the sub-elements of M h/2 only locally around any vertex of M h . Applications of this method to the Navier-Stokes model can be found in [START_REF] Frazza | Nonlinear corrector for RANS equations[END_REF]. We present an application with the Euler model used for sonic boom prediction. We consider again the C25 geometry. The important input is the pressure signature at one-body length below the aircraft.

Figures 4,5 depict the pressure signal and the local error bar, from the non linear corrector, for a sequence of tailored meshes (meshes aligned with the Mach cone, Figure1) and for adapted meshes. The pressure perturbation on axis (z = 0) concentrates between abscissa x = 40 and x = 90. Mesh convergence does not clearly hold after x = 60 for both series of calculations. This is probably explained by the extreme complexity of the many shocks interacting in this region, which are captured in details with the adapted meshes (see Figure2). Let us concentrate on the interval 40 < x < 60. The left-side calculations seem converged as values are concerned but the corrector remains large. The right-side show a more coherent convergence, with a corrector also getting smaller.

The norm-oriented (NO) formulation

In order to obtain a good evaluation of a solution field, the user has to specify the field and also in which norm the approximated field should close to the exact one. Of course the choice of this norm is constrained by an admissibility statement: typically, a second-order accurate unstructured scheme will generally not provide a convergent prediction of the third-order derivative of the solution. Another important output of NO is an evaluation of the approximation error. We describe now a norm-oriented approach relying on the L 2 norm. It consists in the minimization of the following expression with respect to the mesh M :

j(M ) = ||W -W M || 2 L 2 (Ω ) . (3)
where the state W M is solution of the discretised PDE on a unit mesh for M :

Ψ M (W M ) = 0 (4)
Introducing g = W -W M , we get a formulation similar to a goal-oriented formulation:

j(M ) = (g,W -W M ). (5)
The central idea of NO is to replace g by a numerical corrector W h as defined in previous section. The rest of the NO process, follows the GO algorithm with g = W h . The whole NO adaptation algorithm finally writes:

Step 1: solve state equation

Step 2: solve corrector equation

Step 3: solve adjoint equation

Step 4: evaluate optimal metric

Step 5: generate unit mesh for M opt,norm and go to Step 1.

In order to give an idea of how this NO works, we consider as benchmark a test case from [START_REF] Formaggia | Anisotropic a priori error estimates for elliptic problems[END_REF] featuring a 2D boundary layer (Figure6) for comparing a featureoriented calculation with a norm-oriented one. The Laplace equation is solved with a RHS inducing the boundary layer depicted in the figure. When the feature-based method is applied, a tremendous improvement of the error is obtained with 128 vertices, then a uniform element division and a feature-oriented adaptation are applied in alternance. The L 2 error norm (ordinates) is measured from the analytic solution. While the element division is applied, the error is as expected divided by 4. In contrast, for 512, 2048,... vertices (abscissae), the effect of feature-oriented adaptation is to increase the error. This illustrates the fact that asymptotically the direction of a better interpolation deviates from the direction of a better approximation. As a consequence, with a feature-based adaptation, the second-order convergence can be lost. At the contrary, with this test case, each NO mesh-adaptation phase improves (even slightly) the error norm, producing an asymptotic numerical convergence of order two. When 30, 000 vertices are used, the error more than two orders of magnitude lower than with the uniform refinement.

Higher-order and unsteady

Due to their error size and characteristics (dispersion, dissipation), second-order accurate approximations are unable to compute many phenomena. For smooth contexts, high-order methods bring crucial improvements. But, as soon as singularities are involved in the model, they fail in producing a sufficient numerical convergence. Typical real life examples are sound propagation in media containing cracks, or flow with interfaces. Relative errors less than two-three orders of magnitude in continuous media with singularities is a challenge for computational mechanics. To address this challenge, high-order (HO) mesh adaptation methods extending to singular cases need be defined.

Higher-order estimates

Main existing HO schemes satisfy the so-called k-exactness property expressing the fact that if the exact solution is a polynomial of order k then the approximation scheme will give the exact solution as answer. Inside these schemes is involved a step of polynomial reconstruction (e.g. ENO schemes), or of polynomial interpolation (e.g. Continuous/discontinuous Galerkin). The main part of the error can then be expressed in terms of a (k + 1)-th term of Taylor series where the spatial increment is related with local mesh size.

We want to stress that this is the key of an easy extension of metric-based adaptation to HO schemes. We illustrate this with the computation of 2D Euler flows. Considering a triangulation of the computational domain and its dual cells C i built with triangle medians, the exact solution W of Euler equations verifies:

B(W,V 0 ) = 0, ∀ V 0 ∈ V 0 = {V 0 constant by cell},
with (omitting initial conditions):

B(W,V 0 ) = T 0 ∑ i E j (W,V 0 ) dt + T 0 ∂Γ F Γ (W ) • nV 0 dΓ dt E j (W,V 0 ) = C i v 0 ∂ π 0 W ∂t dΩ + 1 2 ∂C i V 0 (F (W )| ∂C i + F (W )| ∂C j ) • ndσ
where we denote by π 0 the operator replacing a function by its mean on each cell. Let us define a quadratic Central-ENO scheme [START_REF] Barth | Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction[END_REF][START_REF] Ivan | High-order solution-adaptive central essentially nonoscillatory (CENO) method for viscous flows[END_REF]. The computational cost of this scheme is rather large but acceptable for 2D calculations (its extension to 3D is even more computationally expansive). This scheme is based on a quadratic reconstruction on any integration cell C i using the means of the variable on cells around C i . Let us denote R 0 2 the global reconstruction operator mapping the constantby-cell discrete field into its quadratic-by-cell reconstruction. The CENO scheme writes in short:

Find W 0 constant by cell s.t. B(R 0 2 W 0 ,V 0 ) = 0, ∀ V 0 constant by cell.
A representative functional of goal-oriented error is:

δ j = (g, R 0 2 π 0 W -R 0 2 W 0 ).
Lemma [START_REF] Carabias | Anisotropic goal-oriented error analysis for a third-order accurate CENO Euler discretization[END_REF]: Introducing the adjoint state W * 0 ∈ V 0 , solution of:

∂ B ∂W (R 0 2 W 0 )(R 0 2 V 0 ,W * 0 ) = (g, R 0 2 V 0 ), ∀ V 0 ∈ V 0 . ( 6
)
we have the following equivalence:

(g, R 0 2 π 0 W -R 0 2 W 0 ) ≈ ∂ B ∂ u (W )(R 0 2 π 0 W -W,W * 0 ).
This estimate is typical of a k-exact variational scheme and permits to express the error as a Taylor term of rank k + 1, rank 3 in our case, with respect to directional mesh size δ x, which we replace by the power 3/2 of a quadratic term:

δ j sup δ x T(δ x) 3 ≈ sup δ x | H|(δ x) 2 3 2 ∀δ x ∈ R 2 .
In [START_REF] Carabias | Anisotropic goal-oriented error analysis for a third-order accurate CENO Euler discretization[END_REF] the equivalent pseudo-Hessian Hi is evaluated by least squares, so that the research of an optimal metric can be computed in a similar way to the second-order accurate case. An a priori better option, not tested here, is to apply the strategy of [START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D. 25th International Meshing Roundtable[END_REF].

High-order accurate unsteady mesh adaptation

We illustrate the use of the above estimate (with pseudo-Hessian) with an application to an unsteady flow. For some typical propagation phenomenon, the discretisation grid (space and time) necessary for a complete representation is very heavy. We consider here an acoustic wave propagation based on the Euler equations.

In order to aply an unsteady mesh adaptation, we adopt the so-called Global-Fixed-Point (GFP) algorithm, [START_REF] Belme | Time Accurate Anisotropic Goal-Oriented Mesh Adaptation for Unsteady Flows[END_REF]. In short, the time interval is divided in n adap sub-intervals. After a computation of state and adjoint on the whole time interval, an optimal space-time metric is evaluated as a set of spatial metrics for each of the time sub-intervals.

In Figure 7, the propagation of a noise from a road (bottom left) to a balcony (near top, right) around an anti-noise wall (middle of bottom) is computed. The functional is the pressure integral on an interval of the balcony. Since a few wavelengths are emitted at the noise source, the mesh adaptation process will concentrate on the part of the wave train which will hit the balcony. This dramatically reduces the region of the computational domain which needs to be refined. With 30980 vertices (mean of the 20 meshes used over the time interval) the resolution is about 10 points per half wave and would require 5 millions vertices if the mesh were an uniform mesh of same maximal fineness. Now, with the adaptive calculation, the error on the integral of pressure on the balcony is still about 1%. With this GFP algorithm, we have maximized the safety by always choosing meshes as fine as or finer than necessary. The GFP involves a delicate compromise between efficiency and accuracy. Indeed, if the number of meshes during the timeinterval is low, then each mesh is adapted to the evolution of solution over a large interval and is therefore heavier than with a smaller time-interval and of course even still heavier than a mesh adapted to a precise time step. This disadvantage amplifies if a thin detail (a discontinuity for example) in the solution is traveling, since the mesh will be extremely refined over the whole region swept by the detail. Conversely, if many adapted meshes are used in the time-interval, then there will be many transfers from a previous mesh to the next one. Now these transfers have been identified as deteriorating the solution, in an extent which justifies to avoid hundreds transfers in a computation. A way to partly avoid this compromise may be to move the mesh. Then mesh admissibility would be accounted by local repairing at some instant. It is anticipated that these repairings would bring less inaccuracy that global mesh-to-mesh transfers.

A discussion of singularity capturing

The sign of a complete success for a mesh adaptive algorithm is that for singular cases the numerical convergence order α is as high as the numerical convergence order α scheme on smooth cases, according to:

|u -u N | ≤ const. N α dim , α = α scheme
where N is the total number of nodes representing the approximate PDE solution and dim the computational dimension (space dimension for steady, space-time dimension for unsteady). Unfortunately, this convergence for singular cases faces many barriers related to the family of meshes considered in the adaptation step, rather independantly of the underlying approximation. Let us mention two important examples of such barriers. In the steady (2D/3D) case, mesh adaptation cannot converge at second-order on a discontinuity if it does not involve anisotropy in arbitrary direction. In the unsteady case (N being the total space-time nodes), convergence at second order is not obtained when time stepping uses a (spatially) uniform time step. For this second issue, some progress should be expected from multi-rate time advancing (see [START_REF] Seny | An efficient parallel implementation of explicit multirate Runge-Kutta schemes for discontinuous Galerkin computations[END_REF][START_REF] Itam | A Volume-Agglomeration Multirate Time Advancing Approach[END_REF]). The next two subsections concentrate on the high-order singularity capturing issue for steady calculations. The third one on unsteady case.

Low order interpolation error (1D)

When a second-order approximation is applied to a singular variable, the output generally show a more or less developed oscillatory zone close the the singularity. About one decade of research around the 80's produced TVD/LED schemes which are thought as fully second-order accurate far from the singularity and first-order near the singularity and non-oscillatory. This accuracy statement is true for local approximations like the interpolation of a given function. But the global error is at most of first order for example in L 1 norm (order zero in L ∞ ) for a series of uniform meshes. For the approximate solution of a PDE, this accuracy statement is also false: in many case it will be globally first-order accurate, in particular in L 1 . It is interesting to analyse how the above FB, GO, NO mesh adaptation technologies behave in presence of a singularity (typically a discontinuity through a smooth surface in 3D).

To discuss the feature-based option we consider the basic problem of meshadaptive interpolation of a the sum of a quadratic function with a Heaviside 1D function h(x) = 1/2(sign(x) + 1). The derivative of h at zero is a Dirac δ 0 and its Hessian is a δ 0 . Looking at the P 1 interpolation on a regular mesh, the absolute value of the Hessian is proportional to 1/∆ x 2 on the two points close to zero, that is on a support of O(∆ x width). After a short calculation of the optimal metric for L 1 interpolation, the optimal strategy divides by four the interval containing zero and by two the other ones, which produces a second-order accurate convergence. This illustrates that although the interpolation error analysis defining the Hessian term as a representative main term of the error is not rigorous, keeping this Hessian in the feature-based adaptation makes a reasonable job.

Let us mention that this analysis extends to the multi-dimensional context and to the goal/norm-oriented approach, and that second-order convergence also holds with curved singularity supports. Let us give an example. A typical problem in multiphase flow is the discontinuity of pressure gradient along an interface. Our model problem is a Poisson problem with a discontinuous coefficient. We apply a uniform refinement, a feature-based adaptation and a norm-oriented mesh adaptation, converging them by increasing the number of vertices. Second-order convergence is obtained while uniform mesh refinement show a desastrous first-order convergence (still 10% error with a 300 × 300 mesh).

Higher order interpolation error

In the case of an approximation order larger or equal 3, the discontinuity may affect the interpolation over an element or an approximation molecule and then can be more expansive. For order three, the main part of the error involves a third-order derivative, proportional to 1/∆ x 3 with a support width of a few ∆ x. Then with a quadratic central ENO scheme, it is interesting to use limiters which return to first order through the discontinuity in order to show a thin shock structure, easier to capture with a small number of nodes.

Passing to multidimensional, let us examine the configuration where the exact solution involves a discontinuity aligned with a rather regular surface. Geometrically, the discontinuity can be put between two lines(2D)/surfaces(3D) of nodes.. This results in a too large number of nodes is the mesh is isotropic, but a smaller number of nodes (allowing second-order convergence) if the mesh is with straight edges and anisotropic, and finally a reasonable number of nodes (allowing higher-order convergence) if the mesh is anisotropic and the edges are curved.

Due to (k + 1)-th order accuracy far from the singularity, a sequence of mesh which divide by two the regular part of the mesh should divide by at least 2 k+1 the elements crossed by the discontinuity where convergence accuracy is less or equal unity. Now, according to the type of high-order approximation, the numerical capturing will be spread over one or several high-order elements. Then a refinement region covering a time interval may cover several times 2 k+1 elements. It may become fairly computational cost consuming. Again an evident way to reduce this cost is to use an approximation with one cell per element, with a shock capturing mechanism which reduce locally the order to one and the numerical shock spread to be between only two lines(2D)/surfaces(3D) of cells.

This strategy is applied for obtaining a steady scram jet flow (Figure 9) with adaptive capturing of the shocks as thin (in terms of width and of cells) as possible. This test case does not involve curved discontinuities for which our adapted mesh with straight edges would face an accuracy barrier (accuracy limited to two). From this standpoint, the path for third-order mesh adaptation indeed involves the mastering of curved elements, see for example [START_REF] Dobrzynski | High order mesh untangling for complex curved geometries[END_REF].

Complexity analysis for the CENO-based Global Fixed Point

The space-wise third-order capture of the singularity can be addressed by spatial mesh adaptation: for each time level, we apply a metric-based anisotropic spatial mesh adaptation. Mesh size normal to the discontinuity will be reduced by a factor 8, solely in the vicinity of the discontinuity. In other directions near discontinuity and in any direction in other regions of the computational domain, the mesh size is divided by only 2. Then the increment in number of nodes near the discontinuity is asymptotically negligible and in the total, passing from N spatial degrees of freedom to 4N leads to a eight times reduction of the spatial error. Unfortunately, for decreasing accordingly time discretization error, the time step also needs to be reduced by a factor 8, leading to a total number of degrees of freedom multiplied by 32. Thus:

||u -u 64N || L 1 (Ω ×[0,T ]) ≤ 1 4 ||u -u N || L 1 (Ω ×[0,T ]) . (7)
By comparing 1 4 with 32 -α/n , this time with n = 3 (space-time dimension), we get a space time convergence order α = 8/5.

In 3D, with this strategy, passing from N spatial degrees of freedom to 8N leads to a eight times reduction of the spatial error. Unfortunately, the time step also needs to be reduced by a factor 8, leading to a total number of degrees of freedom multiplied by 64. Thus:

||u -u 64N || L 1 (Ω ×[0,T ]) ≤ 1 4 ||u -u N || L 1 (Ω ×[0,T ]) . (8)
By comparing 1 4 with 64 -α/n , this time with n = 4 (space-time dimension), we get α = 4/3 which is even smaller.

Lemma: (Barrier of convergence order for time advancing CENO3 discretization) for a time-advancing third-order flow solver coupled with an unsteady goal-oriented mesh adaptation at each time level applied to a traveling discontinuity, the spacetime convergence rate is at most α = 8/5 for 2D unsteady calculations and and α = 4/3 for 3D unsteady calculations. This analysis extends to our Global Fixed Point if the number of time steps is kept constant while the adapted mesh is frozen. In other words, in 3D for example, when we pass from N space-time nodes to 64N, the number of meshes over the time interval must be multiplied by 8, which may fastly amplify the mesh-to-mesh transfer error.

As already mentioned, this barrier can be overcome with a full space-time adaptive time advancing (with space-time elements) or by applying a multirate strategy.

Conclusions

This review of recent progress in mesh adaptation has validated the many promises made by this numerical technology. The promises concern the calculation of new complex (singular/multi-scale, coupled) phenomenas in continuum mechanics, which were not affordable before. This is conditioned by the design of algorithms which show a numerical convergence (in terms of number of nodes used) really better than shown by existing ones. For this purpose, we are increasing our knowledge of the conditions under which mesh adaptive steady high-order methods will enjoy high numerical convergence order for non-regular phenomenas.

Simultaneously, correctors are derived in order to improve our control on error level.

Clearly, for steady calculations, things are progressing rather fastly. The optimal order of convergence can be obtained and this results in the success of computations, like sonic boom, which were not possible without automatic mesh adaptation. But the highest expectation concerns unsteady phenomena. Some preliminary successes show that much can be expected in this direction. At the same time, many difficulties remain to be solved before unsteady mesh adaptation methods show the maximal numerical convergence order.

The roadmap for being able to compute new phenomenas includes of course the further developement of efficient unstructured higher-order approximations. This issue, not discussed here, is progressing (see an interesting discussion in [START_REF] Vincent | Facilitating the adoption of unstructured high-order methods amongst a wider community of fluid dynamicists[END_REF].

Once the number of nodes for the expected accuracy is made small, it remains to reduce the number of arithmetic operations to solve the discrete PDE system. Novel solution algorithms should be combined with the mesh-adaptation loop in order to produce the most performing global mesh-adaptation/system resolution algorithms. Algorithms combining O(N) complexity and parallel scalability are compulsory. Also, the iterative PDE solver should be a good corrective solver, that is a solver which leads efficiently from a rather good initial field to a better/just as necessary converged solution field. The steady case is again an easier case and we refer to [START_REF] Brethes | A mesh-adaptative metric-based full multigrid for the Poisson problem[END_REF] for a FMG-adaptative proposal.

Fig. 1

 1 Fig.1Lowboom C25 Geometry (left) and typical computational domain, around this geometry, aligned with Mach cone, used for the "tailored mesh"(right).

Fig. 2

 2 Fig. 2 Lowboom C25 example: Cut plane x = 30 (top left) and x = 50 (top right) and meshes on the symmetry plane (bottom left) and close view near the inlet (bottom right).

Fig. 3

 3 Fig. 3 Goal Oriented mesh adaptation of the viscous flow in a channel with diamond-shape obstacle. Global mesh (top), global pressure solution (middle) and zoom on the beginning and thickening of the bottom boundary layer (bottom).

Fig. 4

 4 Fig.4Lowboom C25 example: Pressure levels (red) and non-linear corrector (black) predictions of the pressure on tailored meshes.

Fig. 5

 5 Fig.5Lowboom C25 example: Pressure levels (red) and non-linear corrector (black) predictions of the pressure on adapted goal-oriented meshes .

Fig. 6

 6 Fig. 6 Elliptic test case of a 2D boundary layer. A comparison between uniform refinement ("FMG"), feature/Hessian-based, and norm-oriented mesh adaptation methods: error |uu h | L 2 in terms of number of vertices.

Fig. 7

 7 Fig. 7 Goal Oriented unsteady calculation of nonlinear acoustics propagation with third-order goal-oriented adaptation. Pressure at three different time levels and the corresponding meshes.

Fig. 8

 8 Fig.8Mesh-adaptive calculation of a pressure with discontinuous gradient due to an interface. The second order mesh convergence (in terms of number of vertices) of a feature-based adaptation and a norm-oriented adaptation are compared with the first-order uniform refinement.

Fig. 9

 9 Fig. 9 Mesh-adaptive third-order accurate calculation of a scramjet flow, Mach contours, global view, and zoom around the bottom leading edge (improved palette).

Table 1

 1 ). Lowboom C25 example: Statistics for the sequence of adapted meshes. Ratio holds for the aspect ratio defined on any element as max i h i /min i h i . Quotient is defined on any element as max i h 3 i /(h 1 h 2 h 3 ) and measures the gain with respect to an isotropic mesh.

	# vertices	# tetrahedra mean ratio mean quotient
	1 946 918	11 093 526	10	97
	5 529 916	32 726 736	18	297
	12 953 271	77 497 813	18	325
	27 644 454	166 148 147	19	400

Acknowledgements

Alain Dervieux's thesis was advised by Roland Glowinski. Happy birthday Roland and many thanks. This work has been supported by French National Research Agency (ANR) through project MAIDESC n o ANR-13-MONU-0010. This work was granted access to the HPC resources of CINES under the allocations 2017-A0022A05067 and 2017-A0022A06386 made by GENCI (Grand Equipement National de Calcul Intensif).


Mariam El Khatib

Philippe Ferrandis

Erwan Morvan

Cea Leti

Guillot Gerard

Mariam El-Khatib

Gérard Guillot

Georges Bremond

with a buried Mg doped layer characterized by deep level transient spectroscopy and photoluminescence

AlGaN/GaN metal-insulator-semiconductor capacitors

GaN material holds an advantageous position in the fabrication of power devices. This advantage is manifested by the possibility to perform GaN based devices working in high voltage, high current, high frequency and high temperature conditions. However, despite these theoretical forecasts, trapping mechanisms limit the performances of the GaN based devices revealed by the so-called "drain current collapse". Our study is based on a methodology to understand trapping mechanisms in GaN metal-insulator-semiconductor high-electron mobility transistors. This work was achieved by means of electrical and optical characterization techniques such as Fourier transform deep level transient spectroscopy and photoluminescence. The activation energy and the apparent capture cross section of eight traps were extracted in normally-off (Ids=0A when Vgs=0V) AlGaN/GaN heterostructure technology used for power conversion. Six of these traps, E1=0.16eV, E2=0.31eV, E3=0.46eV, E4=0.5eV, E5=0.64eV and E6=0.79eV are electron traps located in the channel. An identification has been proposed for each trap. Two holelike traps, H1=0.17eV and H2=0.74eV were assigned to the Mg and C doping of the GaN buffer layers, respectively. These traps might play a role in the current collapse which appears after the application of a large reverse voltage on the gate of the device. Furthermore, the results obtained using electrical and optical techniques allowed concluding that oxygen atoms and dislocations are incorporated in GaN layers during the growth.

I.Introduction

Efficiency and energy management, with the least losses, present today one of the most famous global concerns. Particularly, works are in progress to increase performances of energy conversion systems [1]. Unfortunately, the efficiency of these systems is limited by restricted properties of traditional switching silicon based devices [2]. Therefore, over the last years, due to their inherent material characteristics, gallium nitride has attracted a particular attention to perform power high-electron-mobility transistors (HEMTs) [3]. The large band gap of gallium nitride provides it promising properties such as a high breakdown voltage and a high temperature operation capability, useful for high-power applications [START_REF] Meneghini | Power GaN Devices[END_REF]. In high-power switching systems, normally-off characteristics are favored to realize the fail-safe requirement [START_REF] Ishida | [END_REF]. For this purpose, several architectures and processes have been proposed such as thin AlGaN barrier, Gate Injection Transistor [6], Fluorine ion treatment [START_REF] Hamady | International Symposium on Power Electronics, Electrical Drives, Automation and Motion (SPEEDAM)[END_REF], p-GaN gate structures [START_REF] Efthymiou | [END_REF] and metal-insulator-semiconductor (MIS) recessed gate structures [9,10]. Despite these technological efforts, the dynamic on-resistance degradation and the current collapse remain open issues to overcome because they seriously reduce the switching performance of devices [11][12]. These effects are reported to be related either to surface states [11,13] or to charge trapping at deep levels [14]. As a result, a great deal of interest has been generated in detecting and identifying these trapping centers and in determining in which part of the device structure they are located. In this work, we studied GaN MIS capacitors using a buried Mg doped GaN layer to isolate the channel from the semi-insulating C doped GaN buffer.

II. Experimental methods

Electrical study was carried out with MIS capacitors grown by metalorganic chemical vapor deposition on a 200nm diameter Si substrate, oriented (111). The epitaxy consists of an AlN nucleation layer, followed by AlGaN buffer layers, then a 1.8μm thick semi-insulating carbon doped GaN buffer layer, 100nm of Mg doped GaN (Na = 310 17 cm -3 ), a 200nm unintentionally doped GaN layer acting as the channel, a 1 nm unintentionally doped AlN layer and a 21.4nm unintentionally doped Al0.20Ga0.80N barrier layer. The Mg-activation anneal of the samples was done at 820°C. Si3N4 was used for the passivation of the stack (Fig. 1a). Interdigitated comb-like capacitors with six 151000μm 2 fingers were then formed using TiN/W as gate metal and Ti/Al annealed for 30s at 875°C as ohmic contact (Fig. 1b). Recessed gate electrodes were achieved by etching Si3N4 passivation and 35nm of the AlGaN/GaN bilayer using a dry etching process. 30nm of Al2O3 were used to form the gate oxide. Static electrical characteristics and deep level transient spectroscopy data were acquired with a Fourier transform deep level transient spectroscopy (FT-DLTS) system from PhysTech (FT 1030) using a Boonton 72B capacitance meter with a 100mV test signal at 1MHz. Capacitance transients were recorded between 80K and 450K in the warm up mode. Optical measurements were performed with four samples consisting of 1.5µm thick doped GaN layer grown on Si. One sample has a Mg doping at 110 17 cm -3 and three samples have a C doping at 510 15 cm -3 , 510 16 cm -3 and 510 17 cm -3 , respectively. The PL steady state measurements of the doped GaN layers grown on Si were performed by means of a solid state laser (CryLas -Berlin) using a continuous excitation source at 266nm wavelength and an excitation power density in the 1-10W.cm -2 range. PL emissions were dispersed by a singlegrating monochromator (Jobin Yvon HR 640). All spectra were corrected by the transfer function of the optical setup including cryostat window, lens and monochromator. Currentvoltage (I-V) measurements recorded at various temperature below 380K (Fig. 3) reveal a low reverse current of 8nA at a bias of 0V which emphasizes the good electrical quality of the MIS capacitors.

III. Experiments and results

A. Static measurements

B. DLTS analysis

In Fig. 4, we see DLTS spectra recorded with a reverse bias VR = 0V, a filling pulse height VP = 2V, a filling pulse width tP = 1ms and four period widths TW between 2ms and 5s. Taken into account the pinch-off voltage found in Fig. 2, we determine that the probe zone is located in the 2DEG channel of the device. Six electron traps labelled E1, E2, E3, E4, E5 and E6 emerge in Fig. 4. Arrhenius plots reported in Fig. 5 allowed extracting the features of each trap. Activation energies and apparent capture cross sections are summarized in Table 1. C related [31,32] Table 1: Activation energies and apparent capture cross sections of the deep levels found in this work and their proposed identification according to previous works.

Traps E5 and E6 are commonly observed in epitaxial GaN and are related to gallium vacancy (VGa) and interstitial nitrogen (NI), respectively. Traps E1 and E3 which are Nitrogen vacancy (VN) related, are expected to be formed by the etching process. The inductively coupled plasma (ICP) etched surface was found to be N-deficient, meaning that a number of VN were produced by ICP etching [33]. Trap E2 was associated to reactive ion etching (RIE) surface damage and has already been related to the effect of an ICP-RIE etched surface [34]. Despite the absence of carbon during the growth of the channel, residual carbon atoms could incorporate the epitaxial layer especially at the vicinity of the carbon doped GaN buffer layer. This trend could explain the appearance of trap E4. In Fig. 6, we see the evolution of the spectra as a function of the filling pulse duration between 100µs and 20ms. DLTS data were recorded with a reverse bias VR = 0V, a filling pulse height VP = 2V and a period width TW = 200ms. Figs. 7 and 8 report the evolution of the peak amplitude for traps E1, E2, E3, E4 and E6 versus the filling pulse duration. We observe a saturation of the DLTS signal above 5ms for traps E1, E2, E3 and E4, which is expected for point defects. However, trap E6 exhibits a logarithmic dependence of the peak amplitude on the range of tP used (Fig. 8b). Such a trend is consistent with a line defect behavior [35]. Previous works support this result by associating this trap to threading dislocations [36,37]. The nature of trap E5 was not determined due to the difficulty to find the peak maximum on the spectra of Fig. 6. Despite the detection of traps E1 to E6 in the channel, these deep traps are not expected to play a role in the current collapse [38]. Keeping a filling pulse height VP = 2V, we applied a reverse bias VR = -10V to increase the depletion region of the capacitor and probe the carbon doped GaN buffer layer. In Fig. 9, we see that using a filling pulse width tP = 1ms and a period width TW = 1s, two dominant holelike traps labelled H1 and H2 are observed. Arrhenius plots extracted from DLTS data in Fig. 9 are reported in Fig. 10 and the signature of the traps are summarized in Table 1. H1 and H2 are Mg and C related defects, respectively. Trap H1 is assigned to MgGa acceptors [39] and trap H2 to CN-ON complexes [31,32]. Minority carrier injection is necessary to observe these hole-like traps. It can come either from the Mg doped GaN which is p-type even if this layer is assumed to be largely depleted or from the semi-insulating C doped buffer layer which is expected to be slightly p-type [40]. The domination of traps H1 and H2 in Fig. 9 indicates that (i) the reverse bias VR = -10V allows extending the space charge region beyond the Mg doped GaN until the C doped GaN buffer layer, (ii) oxygen atoms are incorporated in the epitaxial GaN layers during the growth and (iii) these levels are more efficient to trap carriers than traps E1 -E6.

C. Photoluminescence study

Photoluminescence measurement of the Mg doped GaN sample carried out at 12K is reported in Fig. 11. ABE1 consists of the neutral A1 acceptor bound exciton line [41] and is related to MgGa acceptor [39]. DAP1 is a donor-acceptor pair luminescence involving MgGa acceptor.

The other lines at lower energies represent a set of replicas separated by the energy of LO phonon in GaN which is about 91meV. DBX3-LO represents the LO phonon interaction with a donor bound exciton line. DBX3 is an overlay of an ionized donor bound exciton (D + X) related to DBX2 and a deep donor bound exciton [42]. The neutral donor bound exciton DBX2 is related to the oxygen donor [42]. In Fig. 12, we can also see the Y7 line with the characteristic LO phonon replicas. The Y7 line is associated to an exciton bound to a structural defect [43]. Reshchikov and Morkoç [43] suggested that the Y7 line arises from recombination of an exciton bound to some point defect that is in turn trapped by the stress field of the threading-edge dislocation. The photoluminescence result reported in Fig. 11 is consistent with the DLTS analysis described in section B since both characterization techniques allow concluding to the presence of an acceptor level assigned to MgGa and the existence of oxygen atoms and dislocations in GaN layers. Figure 12: Photoluminescence spectra of C doped GaN samples recorded at room temperature and 0.6W/cm 2 . Oscillations are caused by interference effects in the Si/GaN/Air cavity.

Fig. 12 reveals photoluminescence spectra of C doped GaN layers carried out at room temperature. A yellow luminescence (YL) band centered at 2.2eV arises on all spectra and is C doping dependent. The YL band has been attributed to the result of electronic transitions implying the CN-ON complex [31,32]. This finding is consistent with the assignment of trap H2 observed by DLTS to CN-ON level (see section B) and confirms the presence of oxygen associated with substitutional carbon atoms in GaN layers.

IV. Conclusion

AlGaN/GaN MIS capacitors have been studied using FT-DLTS and photoluminescence. In the channel region, six electron traps E1 -E6 have been identified. Among them, traps E1, E2 and E3 are expected to be formed by the etching process. E4 could be related to residual carbon atoms at the vicinity of the C doped GaN buffer layer. Traps E5 and E6 are commonly observed in epitaxial GaN. Trap E6 exhibits a logarithmic dependence of the peak amplitude on the range of tP used. This result is consistent with previous works which associate this trap to threading dislocations. Extending the space charge region beyond the Mg doped GaN until the C doped GaN buffer layer, spectra are dominated by two hole-like traps H1 and H2. H1 is related to MgGa acceptors and H2 is ascribed to CN-ON complexes, highlighting the presence of oxygen in GaN layers. The photoluminescence study on Mg and C doped GaN samples supports the existence of traps involving MgGa acceptors, CN-ON complexes and dislocations in GaN epitaxial layers composing the MIS capacitors. The domination of traps H1 and H2 on DLTS spectra suggests that these levels might play a role in the current collapse which appears after the application of a large reverse voltage on the gate of the device.

Figure 1 :

 1 Figure 1: (a) Schematic of a normally-off AlGaN/GaN MIS capacitor and (b) planar image of MIS capacitors studied in this paper.

  Typical capacitancevoltage (C-V) characteristics of GaN MIS capacitors recorded at different temperatures are depicted in Fig.2. When the temperature increases, the pinch-off voltage position shifts towards lower voltages, suggesting charge trapping effects between GaN (200 nm) 80K and 400K. Room temperature C-V measurements (shown in the inset of Fig.2) point out that the diode exhibits a low voltage hysteresis (less than 40.7mV) which results of a very low amount of residual charges in the Al2O3 insulating layer.

Figure 2 :

 2 Figure 2: C-V characteristics carried out at different temperatures. The inset reports an hysteresis effect obtained by sweeping up and down the voltage at 299K.

Figure 3 :

 3 Figure 3: I-V characteristics carried out at different temperatures. The inset reports a low reverse current of 8nA with a bias of 0V at 297K.

Figure 4 :

 4 Figure 4: DLTS measurements of AlGaN/GaN MIS capacitors recorded with VR = 0V, VP = 2V, tP = 1ms and four TW varying from 2ms to 5s.

Figure 5 :

 5 Figure 5: Arrhenius plots extracted from DLTS data in Fig. 4.

Figure 6 :

 6 Figure 6: DLTS measurements of AlGaN/GaN MIS capacitors as a function of the filling pulse duration.

Figure 7 :

 7 Figure 7: Evolution of the peak amplitude of trap E1 versus the filling pulse duration extracted from data in Fig. 6.

Figure 8 :

 8 Figure 8: Evolution of the peak amplitude of traps (a) E2, E3, E4 and (b) E6 versus the filling pulse duration extracted from data in Fig. 6.

Figure 9 :

 9 Figure 9: DLTS measurements of AlGaN/GaN MIS capacitors recorded VR = -10V, VP = 2V, tP = 1ms and TW = 1s.

Figure 10 :

 10 Figure 10: Arrhenius plots extracted from DLTS data in Fig. 9.

Figure 11 :

 11 Figure 11: Photoluminescence spectrum of Mg doped GaN sample recorded at 12 K and 1W/cm 2 .

Electronic mail: philippe.ferrandis@univ-tln.fr

Acknowledgment

This work was done in collaboration with "CEA-LETI" and is supported by a doctoral fellowship (ADR to M. EL KHATIB) from the "Région Auvergne-Rhone-Alpes", so funding from the "Région Auvergne-Rhone-Alpes" via the "ARC4" is acknowledged.


Emmanuelle Itam
email: emmanuelle.itam@umontpellier.fr

Stephen Wornom
email: stephen.wornom@inria.fr

Bruno Koobus
email: bruno.koobus@umontpellier.fr

Alain Dervieux
email: alain.dervieux@inria.fr

COMBINING A DDES MODEL WITH A DYNAMIC VARIATIONAL MULTISCALE FORMULATION

Keywords: turbulence, RANS, DDES, hybrid RANS/LES, Variational Multi-scale, cylinder, tandem

This paper studies the improvement of DDES models by locally replacing the LESlike component with a Dynamic Variational Multiscale model. Improvements are measured for flow around cicular cylinders at Reynolds numbers ranging from 20000 to 1 million and for tandem cylinders (from a distance between the cylinders of 3.7 diameters up to 12 diameters).

Introduction

From the beginning of turbulence modelling, researchers have tried to find a universal model which they could propose to engineers. Today, hybrid RANS/LES models, among them DDES models, attempt to obtain higher universality by applying different theories (statistics, filters) in different regions of a complex flow. This quest for universality is illustrated in the litterature by applications of these models to a large set of turbulent flows. For the assembly of such composite models, the first difficulty is to a priori identify which part of the flow is best suited to each model. A second difficulty is to choose the models to combine, in particular, to choose which model in a given category (RANS, LES). The DDES method answers the first question and has been adopted by many teams. Concerning the RANS model to combine, the Spalart-Allmaras model, has often been replaced by a two-equation closure, like k -ω. Indeed, it is expected that the many improvements performed on the latter sophisticated RANS closure will have a positive influence on the predictivity of the final model. DDES involves a LES-like model, based on the turbulent viscosity, which is controlled indirectly by a measure of local mesh size. The turbulent viscosity then acts as a Smagorinsky viscosity. When the viscosity coefficient is too low, this model is equivalent to no modelling. We shall show in the sequel that this is not the best model. When the viscosity coefficient is too high, this introduces a too high damping and this results in a bad modelling. When the viscosity coefficient is of convenient strength, it keeps the usual disadvantage of damping large scales, the dissipative effects of the unresolved scales being taken into account on all resolved scales (small and large), which results in a too low backscatter in the simulation. Several LES models behave in most cases better than a baseline Smagorinsky model, while staying very simple (no extra equation). Among them, the dynamic subgrid-scale (SGS) model and the variational multiscale approach which is aimed at limiting the effects of the SGS closure model to the smallest resolved scales. In [START_REF] Moussaed | Impact of dynamic subgrid-scale modeling in variational multiscale large-eddy simulation of bluff body flows[END_REF] we define DVMS, a combination of both models which is not computer intensive and, at the same time, produces better prediction for the test cases we tried.

The focus of this work is two fold. First, to investigate the possible improvement of a DDES when we locally replace the LES-like component by DVMS. Second, to evaluate DDES as a quasi-universal tool, where quasi-universal in our case will more modestly restrict to a small collection of cylinder/tandem flows with a range of Reynolds numbers corresponding to a subcritical flow regime for the isolated cylinder or the first cylinder in the tandem problem, as well as a case of supercritical cylinder flow.

Modelling

The main ingredients adopted in our modelling platform are the following ones. As RANS component, we use the Goldberg-Peroomian-Chakravarthy [START_REF] Goldberg | A wall-distance-free k -ε model with enhanced near-wall treatment[END_REF] low Reynolds k-eps model which improves near-wall flow predictions, in particular the predictive capabilitiy of adverse pressure gradient flows, compared with more traditional k-eps models. This RANS model is integrated inside a DDES formulation in a similar way to what is done with several k -ω RANS, see for example [START_REF] Gritskevich | Development of DDES and IDDES formulations for the k -ω shear stress transport model[END_REF]. As concerns the LES component, previous investigations [START_REF] Moussaed | Impact of dynamic subgrid-scale modeling in variational multiscale large-eddy simulation of bluff body flows[END_REF] show that our DVMS model, which combines a Germano-type dynamic SGS model with a variational multiscale formulation, gives good predictions for a low computational cost (model of low complexity with a good predictivity on rather coarse meshes). Our focus is the action of DVMS in the wakes after the blunt bodies. In the present study, we restrict our work to cylindrical geometries and we locally apply DDES to the region defined as points located at a distance from the cylinder less than 20% of the cylinder radius, DVMS being otherwise applied. Although this combination is zonal and rather simple, this allows to evaluate the effects on the flow prediction of the application of DVMS in the wake region, instead of the LES-like mode of DDES. We call DDES/DVMS the resulting model. As for the numerics, a second-order accurate scheme stabilized by a tunable numerical diffusion proportional to sixth-order space derivatives is used. Cylinder flows have been studied and measured for at least half a dozen decades. Although many measurements are available in the litterature, it remains difficult to identify from them accurate estimates of the bulk coefficients. Computations also are difficult in that the predictions of bulk coefficients are relatively error-prone, even for Reynolds numbers relatively far from the drag crisis interval. However, these bulk coefficients are the main expected output in engineering studies.

Numerical experiments

C d -Cp b C rms L Θ sep L r S t Experiments Norberg [

Low Reynolds number calculations

The DES and DDES models have been introduced for addressing flows in which the boundary layer is turbulent. In principle, the flow around a cylinder at Reynolds number 20K, which is known as presenting laminar boundary layers, should not be computed with DDES. In all our calculations, RANS and DDES adimensioned inflow conditions (unit reference length and velocity) for k and ε are k = 5 × 10 -5 and ε = 5 × 10 -6 in such a way that the turbulent viscosity is µ t = 4.5 × 10 -5 . Of course the RANS model diminishes this viscosity level according to turbulence decay. The decreasing is accelerated by the DDES model at about 2.5 diameters before the stagnation point. The DDES viscosity becomes about 10 times smaller than for RANS in the boundary layer. Therefore the laminar boundary layer is not too badly addressed by DDES. Moreover, this DDES viscosity becomes really small in the wake, but of Smagorinsky type. Bulk coefficients are given in Table 1. It appears DDES works really better than RANS, but that locally substituting the LES-like component of DDES by DVMS brings a rather sensible improvement.

Flow near criticality

In Table 2 we examine the impact of the turbulence model on the prediction of a still subcritical flow around a cylinder at Reynolds number 140K. Assuming that the incoming flow is laminar, the boundary layer should also be laminar. A RANS model will represent a turbulent boundary layer and will also spuriously damp the turbulent wake. As a result, the bulk coefficients do not agree well with experiments. The behavior of the DDES model is now very similar to RANS in the boundary layer but very different in the wake. We observe a notable positive improvement of the prediction. This improvement is further amplified by DDES/DVMS, but not enough to get a model as good as DVMS.

C d -C p b L

Supercritical flow

The Reynolds number is now 1 million. Mean coefficients were long to obtain and only the DDES/DVMS result is presented. We expect to present also the DDES calculation in our talk.

We then compare our calculation to experiments and to a LES calculation by Kim and Mohan using 6.8M cells. For this test case two Strouhal numbers have been experimentally observed but with the present mesh, the DDES/DVMS computation predicts only the smaller one. [START_REF] Szechenyi | Supercritical Reynolds number simulation for two-dimensional flow over circular cylinders[END_REF] 0.25 0.32 0.35 Goelling [START_REF] Goelling | Experimental investigations of separating boundary-layer flow from circular cylinder at Reynolds numbers from 10 5 up to 10 7 ; three-dimensional vortex flow of a circular cylinder[END_REF] 0.35/0.10 Zdravkovich [START_REF] Zdravkovich | Flow around circular cylinders Vol 1: Fundamentals[END_REF] 0.2-0.4 .1-.15 .2-.34 0.50/0.18 Present simulation DDES/DVMS 0.24 0.12 0.21 0.09 Other simulation LES of Kim and Mohan [START_REF] Kim | Prediction of unsteady loading on a circular cylinder in high Reynolds number flows using large eddy simulation[END_REF] 0.27 0.12 0.28

C d C l -C pb St Experiments Szechenyi
Table 3: Bulk coefficients of the flow around a circular cylinder at Reynolds number 10 6 . The mesh involves 2.85M vertices.

Tandem cylinders with 3.7D distance

We now consider a tandem cylinder with a separation distance of 3.7 diameter. This test case is particularly challenging and useful for understanding particular flows in aeronautics. First, the flow over the first cylinder has to be correctly predicted, second the turbulent wake betwen the two cylinders has to be well captured, which of course conditions the quality of the prediction of the flow around the second cylinder. This implies to be able to have not too dissipative turbulence and numerical models, as well as an appropriate grid. In Table 4, we compare bulk coefficients obtained in our RANS and DDES/DVMS simulations with those of Aybay [START_REF] Lockard | Summary of the Tandem Cylinder Solutions from the Benchmark Problems for Airframe Noise Computations-I[END_REF] and Garbaruk [START_REF] Garbaruk | DDES and IDDES of Tandem Cylinders[END_REF], and measurements by Neuhart [START_REF] Neuhart | Measurements of the flowfield interaction between tandem cylinders[END_REF]. We observe that the RANS model gives rather poor predictions while the results obtained with the DDES/DVMS model are in overall good agreement with experimental data, although the prediction is less accurate for the second cylinder. It should also be noted that the DES and DDES results of Aybay and Garbaruk on finer meshes compare less well with experimental data.

Mesh C 1 d C 2 d C 1
As an example of distribution of a quantity of interest, that of the turbulent kinetic energy (TKE) predicted by the DDES/DVMS model along the wake centerline in the gap between the two cylinders is depicted in Figure 1, which shows a good agreement with BART experimental data.

Tandem cylinders with large distance

The accurate modelling of the wake is even a more important issue when the calculation is applied to the analysis of the interaction between two aligned cylinders with a rather large distance between them. This kind of flow are encountered when studying interaction between risers in water [START_REF] Levasseur | Basin Tests and Numerical Simulations of Wake-Induced Oscillations of Inline Risers in Subcritical/Critical Regime[END_REF]. In the case which we study, this distance is 12 diameters. The mesh (8M vertices) is designed in order to be a Cartesian one in most part of the interval between the two cylinders. This is done for enjoying the superconvergence (up to order 5) of the numerical scheme. In this region the DVMS applies and preserves as much as possible the medium scales of the flow. As a result, the influence of the first cylinder is accurately propagated and causes a 5-6 times smaller drag on the second one (Figure 2). A detailed study of the impact of the wake model on important quantities (in particular turbulence intensity which was observed as very sensitive to the LES model in other cases) will be presented at the conference.

Concluding remarks

In this work, a zonal combination of DDES and DVMS is investigated for the flow prediction around cylinders and tandem cylinders at various Reynolds numbers. This first study tends to show that the application of DVMS in the wake region instead of the LES-like mode of DDES has a positive impact on the considered flow predictions. Compared to more simple LES models, the DVMS model in which the dynamic SGS procedure and the variational multiscale formulation work in a complementary way, allows for a more accurate simulation of unsteady turbulent wakes, resulting in an improved prediction of the global flow. Nevertheless, further investigations need to be done in order to confirm this tendancy. Other flow problems will be considered (backward facing step, channel flow), and a non-zonal combination of DDES and DVMS is also under investigation.

Figure 1 :

 1 Figure 1: Distribution of the resolved TKE along the wake centerline in the gap between the two cylinders for the DDES/DVMS computations (continuous line).

Figure 2 :

 2 Figure 2: Flow around a tandem cylinder with an inter-cylinder distance of 12 diameters and a Reynolds of 140K. Drag history for both cylinders (in black, mean drag of a single cylinder).

Table 1 :

 1 Bulk flow parameters for Re = 20000 flow around a cylinder. The mesh involves 1.8M vertices. C d is mean drag, C p b is pressure coefficient at base, C rms L is lift rms, Θ sep is the separation angle, L r is recirculation length, S t is Strouhal number, Θ sep is the separation angle.

	5]	1.16	1.16	0.47	78	1.03	0.194
	Present simulations
	No model	1.27	1.35	0.61	82	0.96
	RANS k -ε	1.31	1.49	0.75	85	0.50	0.212
	DDES k -ε	1.25	1.21	0.58	86	0.90	0.194
	DVMS	1.18	1.20	0.46	81	0.96	0.196
	DDES/DVMS	1.17	1.13	0.46	82	1.05	0.200

Table 2 :

 2 Bulk quantities for Re = 140000 flow around a cylinder. The mesh involves 0.9M vertices.

	r	S t

Table 4 :

 4 Tandem cylinder with 3.7D distance at Re = 166000: mean drag coefficient and root mean square of the lift time fluctuations. Circled superscripts hold for cylinder 1 and cylinder 2.

	Lrms

Lrms

C 2