Youssef Mesri
email: youssef.mesri@ifpen.fr

Hervé Guillard
email: herve.guillard@inria.fr

Thierry Coupez
email: thierry.coupez@mines-paristech.fr

Automatic coarsening of three dimensional unstructured meshes for efficient reservoir simulations

Keywords: Mesh generation, Tetrahedral meshes, Coarsening, Anisotropy, Multigrid Algorithms, CFD, Boundary layers

This paper describes an algorithm designed for the automatic coarsening of three-dimensional unstructured simplicial meshes. This algorithm can handle very anisotropic meshes like the ones typically used to capture the boundary layers in CFD with Low Reynolds turbulence modeling that can have aspect ratio as high as 10 4 . It is based on the concept of mesh generation governed by metrics and on the use of a natural metric mapping the initial (fine) mesh into an equilateral one. The paper discusses and compares several ways to define node based metric from element based metric. Then the semi-coarsening algorithm is described. Several application examples are presented, including a full three-dimensional complex model of an aircraft with extremely high anisotropy.

INTRODUCTION

In recent years upscaling has become increasingly important for converting highly detailed geological models to computional grids. These geological models usually require fine-scale descriptions of reservoir porosity and permeability on grids of tens of millions of cells to honor the known and inferred statistics of these reservoir properties. The goelogical grids The geological grids of this order are far too fine to be used as simulation grids. Even with today's computing power, most of the full-field reservoir models are of the order of 100000 a factor 100 less than the geological grids. Upscaling has been developed to bridge the gap between these two scales. Given a fine reservoir description scale and a simulation grid, an upscaling algorithm is designed to obtain suitable values for the porosity, permeability, and other property data for use in the coarse grid simulation. Many upscaling methods have been devloped, such as pressure-solver [START_REF] Beggs | Zssigning effective values to simulator gridblock parameters for heterogeneous reservoirs[END_REF], renommalization [START_REF] King | The use of renormalization for calculating effective permeability[END_REF], effective medium [START_REF] King | The use of renormalization for calculating effective permeability[END_REF], power law averaging [START_REF] Whitaker | Flow in porous media I: A theoritical derivation of Darcy's law[END_REF] harmonic/arithmetic mean, local averaging, and homogenization [START_REF] Amaziane | Numerical simulation and homogenization of two-phase flow in heterogeneous porous media[END_REF]; see the reviews of upscaling and pseudoization techniques by Christie [START_REF] Christie | Upscaling for reservoir simulation[END_REF] and Barker and Thibeau [START_REF] Barker | A critical review of the use of pseudo relative permeabilities for upscaling[END_REF], for example. Here, we briefly mention a few of these techniques.

• For the single phase flow, the aim of upscaling is to preserve the gross features of flow on the simulation grid. An algorithm is needed to compute an effective permeability, which will result in the same total flow of the fluid through the coarse homogeneous grid as that obtained from the fine heterogeneous grid. In the pressure-solver method [START_REF] Beggs | Zssigning effective values to simulator gridblock parameters for heterogeneous reservoirs[END_REF] for example, we set up a single phase flow computation with specific boundary conditions and then ask what value of effective permeability generate the same flow rate as the fine-scale computation. The results obtained depend on the assumptions made, particularly with regard to the boundary conditions. If no-flow boundary conditions are used, a diagonal effective permeability tensor can be derived and entered directly into a reservoir simulator. Alternatively, if periodic boundary conditions are employed, a full effective permeability tensor can be obtained [START_REF] White | Computing absolute transmissibility in the presence of fine-scale heterogeneity[END_REF].

• For two-phase flow, it is generally believed that upscaling of the absolute permeability alone is not enough to capture the effects of het-erogeneous on two-phase fluid simulation [START_REF] Muggeridge | Generation of effective relative permeabilities from detailed simulation of flow in heterogeneous porous media[END_REF][START_REF] Durlofsky | A new method for scale up of displacement processes in heterogeneous reservoirs[END_REF], particularly when the correlation length of the heterogeneity not represented on the flow simulation grid is significant compared with the well spacing. A multiphase upscaling technique must be used. The most obvious technique is the use of pseudorelative permeabilities, i.e., pseudos [START_REF] Lake | Pseudofunctions-The key to practical use of reservoir description[END_REF]. The role of pseudorelative permeabilities is to determine the flow rate of each fluid phase out of a gridblock. They relate the flow rate to the pressure gradients between the gridblock and its neighbors, given the average saturation in each gridblock. Both the flow rate and the pressure gradient depend on the details of the saturation distribution within the gridblock. Hence, to obtain a pseudorelative permeability curve, it is necessary to determine the saturation distribution within the block for any given average saturation [START_REF] Barker | A critical review of the use of pseudo relative permeabilities for upscaling[END_REF]; see the review papers by Christie [START_REF] Christie | Upscaling for reservoir simulation[END_REF] and Barker and Thibeau [START_REF] Barker | A critical review of the use of pseudo relative permeabilities for upscaling[END_REF] for the generation of pseudorelative permeabilities.

A major limitation in upscaling is that it usually gives an answer whitout any indication of whether the assumptions made in th obtaining the answer hold. No rigorous theory exists behind the upscaling process. Furthermore, some factors give rise to a concern about whether the upscaled values are good approximations; these include large-aspect ratio gridblocks, significant transport at an angle to the grid lines, and upscaled gridblocks close in size to a correlation length of a heterogeneous reservoir. Compared with single phase upscaling, multiphase upscaling is far less developed and understood. Multigrid algorithms are among the most efficient methods to solve linear and non-linear algebraic problems in computational sciences. However, when non-structured simplicial meshes (triangles in 2D, tetrahedra in 3D) are used, these algorithms suffer from the necessity of constructing a sequence of meshes of different resolution of the same geometry. Two main strategies have been proposed to overcome this problem. The first family of methods uses an algebraic interpretation of the problems and relies only on the graph of the matricial problem corresponding to the PDE's that have to be solved. These methods that present several variations e.g. [START_REF] Ruge | Algebraic Multigrid Methods: in Multigrid Methods[END_REF][START_REF] Vanek | Algebraic Multigrid by Smoothed Aggregation for Second and Fourth Order Elliptic Problems[END_REF] are extremely efficient for linear problems. However for non-linear ones, their utilisation is not obvious. The second family of methods are the so-called geometric multigrid methods where a sequence of real meshes of different mesh sizes are generated. These methods can be used for linear as well as non-linear problems. Geometrical methods can be classified into several families:

• Embedded grids: The simplest manner to build a hierarchy of embedded grids consists in first choosing a coarse grid and then to generate the finer level by element subdivision. Unfortunately, this approach has serious limitations. Indeed, as finer and finer grids are built, they become to some extent less and less unstructured since artificial "macroelement" corresponding to the coarse mesh elements are present in the fine triangulation. These macro-elements may influence the computation and show themselves as artificial internal boundaries in the results [START_REF] Carre | Multigrid strategies for CFD problems on unstructured meshes in "Multigrid Methods VI[END_REF]. In three dimensions, another potential drawback of this method appears as the successive refinements of tetraedra produce a sequence of elements with deteriorating shape quality. Moreover, with this technique, the end user has a limited control on the final fine mesh.

• Volume agglomeration: Another idea closely related to algebraic methods is to use the volume agglomeration technique. In this approach, the discretisation on the fine grid is associated with control volume defined on a dual mesh. Coarser finite volume grids can then be built by agglomeration of the cells; different algorithms to realize this task are available(see e.g. [START_REF] Francescatto | A semi-coarsening strategy for unstructured multigrid method based on Agglomeration[END_REF][START_REF] Lallemand | Unstructured multigriding by Volume Agglomeration: Current Status[END_REF][START_REF] Koobus | Unstructured Volume Agglomeration MG : Solution of Poisson Equation[END_REF]). The volume agglomeration method is certainly one of the simpler and most efficient to deal with nonstructured meshes for first-order differential problems. A weakness of this approach lies in the difficulties to build consistent approximations for second-order differential operators on the coarse dual grids. Some solutions have been proposed [START_REF] Janka | Multigrid methods for compressible fluids[END_REF][START_REF] Guillard | P An Aggregation Multigrid Solver for Convection-Diffusion Problems on Unstructured Meshes[END_REF] but they still have to be fully validated for three dimensional real life problems.

• Non-nested coarse grid methods : An alternative to the agglomeration technique consists in building automatically a hierarchy of coarse grids whose elements are not embedded. In 2-D, the first algorithm proposed for this task was in [START_REF] Guillard | Node-Nested Multi-Grid with Delaunay Coarsening IN[END_REF]. It uses a Delaunay retriangulation associated with a MIS (Maximal Independent Subset) algorithm to delete nodes of the fine triangulation. Similar algorithms aimed at the same task have been proposed [START_REF] Chan | Domain decomposition and multigrid algorithms for elliptic problems unstructured meshes[END_REF], [START_REF] Mavriplis | Adaptive meshing techniques for viscous flow calculations on mixed-element unstructured meshes[END_REF][START_REF] Mavriplis | Multigrid strategies for viscous ow solvers on anisotropic unstructured meshes[END_REF], [START_REF] Adams | A parallel maximal independent set algorithm[END_REF]. The algorithm of [START_REF] Guillard | Node-Nested Multi-Grid with Delaunay Coarsening IN[END_REF] has been followed by 3-D extensions on non-structured isotropic tetrahedral meshes that proved to be very effective for inviscid steady computations [START_REF] Carre | Multigrid strategies for CFD problems on unstructured meshes in "Multigrid Methods VI[END_REF] and other algorithms using closely related idea have also been proposed [START_REF] Müller | Anisotropic adaptation and multigrid for hybrid grids[END_REF][START_REF] Ollivier-Gooch | Coarsening unstructured meshes by edge contraction[END_REF].

However, to the best of our knowledge, none of these geometrical algorithms are fully satisfactory when the fine meshes are highly anisotropic as is the case for instance for the meshes used for Navier-Stokes computations with Low Reynolds Turbulence modeling. The goal of the present paper is thus to extend coarsening algorithms to this type of meshes. For the meshes we will consider in this work, the maximum aspect ratio of the elements can be very large (typically, more than 10 4 for boundary layers problems) and the multigrid smoothers are quite inefficient for these cases. The remedy that is known since a long time is to use semi-coarsening in the direction perpendicular to the maximal stretching.

For non-structured grids, this is not an easy task as one has to devise an automatic way to define this direction and then to coarsen the mesh in this direction only. We will realize this task by using the concept of mesh generation controlled by a metric map. Many authors (for example [START_REF] Coupez | Génération de maillage et adaptation de maillage par optimisation locale Revue Européene des éléments finis[END_REF], [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF], [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications. Part 2 : Application examples[END_REF]) showed that the definition of a metric field simplifies the generation of adapted and anisotropic meshes : The metrics associated to the Riemannian space specify the size of the mesh and the direction of the stretching and then an adapted and anisotropic mesh in the Euclidean space can be represented as an isotropic and unitary mesh in the Riemannian space.

For semi-coarsening of non-structured meshes, the algorithm that we propose is therefore divided into two steps :

The first one is an analysis step where the current (fine) mesh is analysed to identify the directions and sizes of stretching. This is done by computing on each element, the metric tensor where this element is equilateral. Then, these tensors are averaged in such a way that on each node of the triangulation, a metric tensor representing the averaged sizes and stretching directions of the elements that contain this node is defined. Then, the eigenvalues of this metric tensor are modified in such a way that the edge lengths of the equilateral element for this tensor is twice the original edge length in the direction of minimal stretching.

The second step is a mesh generation step. We use the fine mesh and the metric field defined on its node as the background mesh and the size specification field as inputs of any existing mesh generation software using the concept of mesh generation governed by metric specifications. In this work, we have used the MTC software [START_REF] Coupez | Génération de maillage et adaptation de maillage par optimisation locale Revue Européene des éléments finis[END_REF][START_REF] Gruau | 3D tetrahedral unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF] but in principle any other mesh generation (Delaunay, frontal, etc) tool can be used for this purpose provided that it can handle anisotropic elements.

The remainder of this paper is organized as follows : Section 2 recall the main concept of a metric map and deal with the problem of metric interpolation. Then in section 3, we describe our directional semi-coarsening algorithm while section 4 gives a quick description of the mesh generator that we have used. Finally, in section 5, we demonstrate the effectiveness of our approach on several anisotropic meshes displaying huge aspect ratio (as large as 10 4 ).

MESHES AND METRICS

The concept of mesh generation governed by metric specifications is now of current use in this field [START_REF] Tam | Anisotropic mesh adaptation for 3D flows on structured and unstructured grids[END_REF][START_REF] Lohner | Adaptive h-refinement on 3D unstructured grids for transient problems[END_REF][START_REF] Mavriplis | Adaptive meshing techniques for viscous flow calculations on mixed-element unstructured meshes[END_REF][START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent cfd[END_REF][START_REF] Gruau | 3D tetrahedral unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF][START_REF] Li | 3D anisotropic mesh adaptation by mesh modification[END_REF][START_REF] Huang | Metric tensors for anisotropic mesh generation[END_REF]. In this section, we review the fundamental properties and definitions of a metric field that will be useful for our applications. First, we recall the definition of a metric ( [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF]).

Definition 1: A Metric (also Tensor) M in R d×d is a symmetric definite positive real matrix with d the dimension of space. It is then a d×d invertible matrix, its eigenvalues are real and positive, and its eigenvectors form an orthogonal basis of R d . such that

M = V (M) t Λ(M)V (M)
where V (M) denotes the orthonormal matrix corresponding to the eigenvectors of M while Λ(M) is the diagonal matrix of its eigenvalues

The orthogonalization property of the metric allows to define a scalar product of two vectors in R d with respect to a metric M :

(u, v) M = (u, Mv) = u t Mv ∈ R (1)
we now define, the associated norm of a vector in R d with respect to the metric M:

u M = ((u, u) M ) 1/2 (2)
A useful graphical representation of a metric M associates to it the ellipsoid that represents the unit ball for the scalar product (., .) M :

E M = {u/(u, u) 1/2 M = 1} (3)
The equation ( 3) represents a general ellipsoid form, where the main axis are given by the eigenvectors of M and the radius of each axis is given by the inverse of the square root of the associated eigenvalues. More details and definitions can be found in [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF].

Notations

We will use the following notations:

• V (M) an orthonormal matrix corresponding to eigenvectors of a Metric M

• λ 1 (M), λ 2 (M), .., λ d (M) the eigenvalues corresponding to a metric M.

• h 1 (M), h 2 (M), .., h d (M) the local size in each direction of the metric defined by: h i (M) = 1/λ i (M) for i = 1, ..., d.

• For any function f defined on R, we also define the matrix f (M) as the matrix

f (M) = V (M) t f (Λ(M))V (M)
where f (Λ(M)) is the diagonal matrix whose entries are f (λ i (M))

We will suppress the notation (M) when no ambiguity will possible

Metric associated with a simplicial element

Let T h be a tetrahedrization of a polygonal domain Ω defined by a set of nodes N and a set of simplicial elements ( triangles in 2-D, tetrahedra in 3-D) T whose nodes belong to N . It is quite obvious (and this fact is the basis of the concept of generation of a mesh governed by a metric field, see [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF] ) that the set of simplicial elements defined a P 0 field of metrics on Ω such that for this metric field each simplex is equilateral : for a given element T , define the set E(T ) of all edges of T ,

E(T ) = {(k, l) : k = l, k, l ∈ N , k, l ∈ T }.
Then the element metric M T on T should satisfy

x T kl M T x kl = 1 (4)
for all edges (k, l) ∈ E(T ) and x kl = x lx k where x k are the coordinates of the node k. This condition gives a system of 1, 3 or 6 linear algebraic equations (for 1D, 2D and 3D, resp.) for components (M T ) ij , i ≤ j of the metric matrix. The solution of this linear system is given in the following result :

Lemma 2.1. The element metric matrix M T can be computed as follows:

M T = C M •    d+1 i,j=1 i<j x ij x T ij    -1 , (5)
with C M = (d + 1)/2 and d the space dimension.

Metric associated to the nodes of a tetraedrization

A mesh generation algorithm governed by metric specifications uses a metric field to define at any point of the domain Ω the desired mesh size and stretching directions. However, this metric field has to be defined in some way. In some cases, it is possible to define this metric field by splitting Ω into several subdomains and defining analytically the metric field over the subdomains. However, in practice, it is more convenient to define this field on an initial mesh, sometimes called the background mesh. This background mesh is only useful for defining at any point of the domain the metric specification and do not have to be confused with the actual mesh that we want to generate. In our coarsening application, the background mesh will obviously be defined by the fine mesh. Using the result of section 2.2, it is clear that we can define a P 0 (constant by element) metric field on the background mesh. However, our coarsening algorithm will remove nodes of the fine triangulation and so thus will delete the elements of the background triangulation. It is thus much more convenient to define a P 1 metric field on the backgound mesh and thus to define the metric field at the nodes of the background (fine) triangulation. Moreover, to the best of our knowledge, the anisotropic mesh generation softwares currently in use or in development [START_REF] Borouchaki | Delaunay mesh generation governed by metric specifications[END_REF][START_REF] Coupez | Génération de maillage et adaptation de maillage par optimisation locale Revue Européene des éléments finis[END_REF] take as input a metric map defined only over the nodes of the background mesh because it is easier to transport P 1 fields during mesh generation process than P 0 fields. We thus have to consider the problem of defining from a P 0 tensor field a P 1 tensor field that, in some sense has to represent the same field : Given for each element, a metric that defines the length of edges and the stretching directions, we have to define a set of metrics defined on the nodes of the triangulation that will define for any point of the domain approximatively the same lengths and stretching directions. For this, we have considered three different strategies that are described and compared in the following sections.

Interpolation of element metrics to node

A first way is to consider a simple interpolation of the metrics defined on the elements that surround a common node. Let us denote T (i) the set of all elements T which contain a node i, i ∈ T . The nodal metric matrix of node i is simply obtained by:

M i =   1 card(T (i)) T ∈T (i) M -1/2 T   -2 (6)
where M -1/2 is a function of the matrix based on the diagonalization of M = V.D.V T calculated by taking M -1/2 = V.D -1/2 .V T with a similar definition for M -2 . The powers -1/2 and -2 have been chosen to give to the expression (6) the meaning of averaging the characteristic mesh size over the elements T ∈ T (i) adjacent to i.

The nodal metric defined by ( 6) has been constructed on the element metrics, and as such should take into account the local shape and direction of anisotropies of adjacent mesh elements. However, numerical experiments show that it does not give, in some cases, quite exact information on the average aspect ratio of the adjacent elements and that it oversizes the smallest characteristic mesh size (eg. the boundary layer thickness). This situation happens, for example, on a curved boundary layer with a high aspect order (order 10 4 or more). For such a case, we propose a correction of the eigenvalues of the nodal metric matrix M i computed in [START_REF] Carre | Multigrid strategies for CFD problems on unstructured meshes in "Multigrid Methods VI[END_REF]. For each metric matrix M, let us denote its diagonal matrix by D with eigenvalues ordered in an ascending order on the diagonal and an orthonormal matrix of corresponding eigenvectors by V ,

M = V.D.V T
With reference to expression (6) let us compute a modified diagonal matrix D i as follows

D i =   1 card(T (i)) T ∈T (i) D -1/2 T   -2 (7)
The final nodal metric is thus given by

M i = V i .D i .V T i (8)
where V i is the orthogonal matrix of eigenvectors of the original matrix M i used in formula [START_REF] Carre | Multigrid strategies for CFD problems on unstructured meshes in "Multigrid Methods VI[END_REF]. This procedure thus replaces the eigenvalues of M i by the corrected eigenvalues given by ( 7) that tries to have a better estimate of the aspect ratio of the elements at node i.

ANISOTROPIC COARSENING ALGORITHM

In this section we describe briefly the main steps of an automatic coarsening algorithm taking into account the mesh anisotropy and discuss some practical points regarding its implementation.

An anisotropic coarsening algorithm

The anisotropic coarsening algorithm follows three main steps :

1. Generate on each node of the finest mesh an initial nodal metric that reflects the size and stretching of elements belonging to this mesh.

2. Modify the initial metric to establish a corresponding coarsened mesh metric. This is done by modifying the eigenvalues λ i associated to the metric, which will modify the mesh size in the desired direction which is the corresponding eigenvector V i .

3. Provide the background (fine) mesh and the desired metric field to a mesh generation tool using metric specifications.

The modification of the eigenvalues necessary to specify the new metric is realized using the following algorithm :

Algorithm 1: Anisotropic semi-coarsening al- gorithm Input: λ 1 , • • • , λ d eigenvalues of a M k Output: Update λ 1 , • • • , λ d coarsened eigenvalues For each i ∈ N Set: h i k = (λ i k ) -1/2 , k = 1, . . . , d, h i 1 ≤ . . . ≤ h i d and h 0 = C cf .h 1 . h i k = max(h i k , min C cf • h i k , h i k-1 ) . Define new eigenvalues λ k ← h -2 k .
End for.

In this algorithm, C cf is the coarsening factor (in all our applications, taken equal to 2). It is seen that this algorithm doubles the mesh size only in the direction where the mesh size is minimal in the region of high anisotropy where h i 1 << h i 2 while it doubles it in the three space directions in the region of isotropy where

h i 1 ∼ h i 2 ∼ h i 3

Modified anisotropic coarsening algorithm

The previous anisotropic coarsening algorithm produces a semi-coarsening in the boundary layer region and a total coarsening in the isotropic region. However, as the different meshes are progressively coarsened, it may happen that this algorithm produces an extremely rapid change in the metric specifications. To understand this point, consider Figure ??. It represents on the left, a configuration of a mesh zone between boundary layer and isotropic mesh zone. On the right, we have shown a representation of the metric at the nodes i and j respectively. Observe that the x-radius of the two ellipses are equal. The previous algorithm will identify node i as being in the isotropic region and consequently, will ask for a doubling of the radius of the ellipsis in the x and y directions. However, point j will be identified as being in the boundary layer and therefore, the metric will be modified only in the ydirection. The metric specification for the new coarsened meshes is therefore represented in figure ??. It is clear that we have a conflict between the two opposite requirements given by these metrics. This situation has been studied in other context in [START_REF] Borouchaki | Mesh Gradation Control[END_REF] and mesh generation algorithms have difficulties to handle these cases. A solution is to allow a smooth evolution of the metric field. There are many different ways to obtain a regular metric map. For instance, one can use smoothing of the metric field by averaging the different metric on an enlarged stencil around a node using the techniques described in section 2.3. Here, we simply suggest an improved algorithm that allows to coarsen smoothly the intermediate region between anisotropic and isotropic mesh regions.

This algorithm takes into account the specification of the nodal metric on the neighbors of a given node and do not allow a fast variation of the metric specifications. This algorithm is the following :

Algorithm 2: Modified anisotropic coarsening algorithm Input:

λ 1 , • • • , λ d eigenvalues of a M k Output: Update λ 1 , • • • , λ d coarsened eigenvalues For each i ∈ N Set: h i k = (λ i k ) -1/2 , k = 1, . . . , d, h i 1 ≤ . . . ≤ h i d . Set h 0 = C CF • h 1 .
Store the initial size:

h i k,old = h i k .
// compare the current node size specifications with its neighbors

if max(h i d , min C cf • h i d , h i d-1 ) = C cf • h i d and ∃j 0 < i ∈ V (i) such as h j 0 1 < h i 1 . then h i k = h i k,old , k = 1, . . . d Define new eigenvalues λ k ← h -2 k .
End for

THE MTC MESH GENERATOR

The coarsening algorithm developped in this work can be used with any mesh generator using the concept of generation governed by a metric map (or alternatively by stretching directions). These mesh generators can use the Delaunay principle, advancing fronts or layers or any other kind of mesh generation method. However, in practice, we have found that for threedimensional geometries, the number of mesh generation softwares able to handle the extremely anisotropic meshes typical of CFD applications is extremely limited. In this section, we give a short description of the MTC mesh-generation tool that is indeed able to produce and modify meshes with extremely high aspect ratio. MTC is being developed at Ecole des Mines de Paris, Centre de Mise en Forme des Materiaux, Sophia Antipolis. It is based on the idea to improve iteratively, an initial unsatisfactory mesh by local improvements. The general algorithm can be expressed as follows, for further details see [START_REF] Coupez | Génération de maillage et adaptation de maillage par optimisation locale Revue Européene des éléments finis[END_REF][START_REF] Gruau | 3D tetrahedral unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF].

Meshing and re-meshing processes in MTC

MTC mesh generator re-mesh the initial mesh iteratively by a local mesh optimization technique. The mesh optimization technique consists in local re-meshing of cavities formed by small clusters of elements in order to increase the "quality" of the elements of the cluster.

In the re-meshing process, two principles are enforced :

• Minimal volume, that insures the conformity of the mesh, with no overlaps of elements: let T i (C) denote the i-th set of elements T filling the local cavity. Following the minimum volume principle we choose as an optimal (possibly not unique) re-triangulation of the cavity the one satisfying

T ∈T i (C) |(V olume)(T ))| → min, (9)
where the minimization is done over a small set i = {1, . . . , I} of possible triangulations T i (C) of elements (Fig. ?? right) connecting the nodes on the border of the cavity, or other nodes like the cavity barycenter, with all boundary faces.

• The geometrical quality Q(T ), which is evaluated for each element.

If the minimizer of ( 9) is not unique, this criterion picks among all admissible cavity re-triangulations the one improving the geometrical quality of the mesh by improving the quality of the worst element of the triangulation.

While the former criterion assures the conformity of the mesh, if the initial mesh was conform, the latter handles improvements of element shape, size, connectivity, etc., depending on the quality function Q(T ). Usually, the quality function Q(T ) is a function of the geometry of the element T and the prescribed background metric, which give together a measure for the element size and the element form (aspect ratio).

Definition of the quality function

Let C n = {T 1 , • • • , T n } ⊂ T h be a set of n elements. We define the quality of this set with respect to a given metric M ∈ R d×d as a real n-vector

Q n = {q(T 1 ), • • • , q(T n )}
The quality q(T ) of an element T of C n , measured in the metric M, is defined as the product of two factors :

q(T ) = Q F (T ) • Q S (T ) (10)
The first one Q F (T ) controls the shape of the element T . The best possible element being the equilateral element (in the metric M). It is defined by :

Q F (T ) = C 0 (V olume)(T ) M (T ) h d M (T ) (11)
where

• d is the space dimension,

• (V olume)(T ) M (T ) is the volume of T measured in the metric space, and is given by:

(V olume)(T ) M (T ) = (V olume)(T ) det(M(T )).
The matrix M(T ) is obtained by averaging the nodal metric matrices on nodes of the element T ,

M(T ) = 1 d + 1 i∈T M -1/2 i -2
, and correct its eigenvalues by the same process as in Section 2.3.1.

• h M (T ) is the average of lengths of edges of T measured by the metric M(T ),

h M (T ) =   2 d(d + 1) (i,j)∈T (M(T )( x j -x i ), ( x j -x i ))   1/2 (12)
• in [START_REF] Fletcher | Principal Geodesic Analysis on Symmetric Spaces: Statistics of Diffusion Tensors[END_REF] C 0 is such that 11 is equal to 1 when T is equilateral in the metric M(T ).

The second factor Q S (T ) controls the size of the element in the metric M. Its definition is given by :

Q S (T ) = min( 1 h M (T ) , h M (T ) ) d (13)
where x i ∈ R d are the coordinates of the node i.

With these definition of the quality of a set of elements, we just need now a way to compare two different sets in order to pick the "best" one. This is done by defining a lexicographic order " < ′′ between two different sets of elements.

Algorithm (MTC iteration)

With the concept previously defined, the mesh generation in MTC will proceed by successive improvements of given mesh. Let us denote T h , E h , N h , respectively the sets of all mesh elements, edges and nodes. Repeat for different cavities C 0 ∈ T h composed by a group of adjacent elements obtained as the nearest neighborhood of a node n ∈ N h or of an edge connecting 2 nodes

(n 1 , n 2 ) ∈ E h , n 1 , n 2 ∈ N h .
1. Denote N 0 ⊂ N h and E 0 ⊂ E h respectively, the set of all nodes and all edges of the cavity C 0 .

2. Denote ∂E ⊂ E 0 the set of all edges belonging to the border of the cavity C 0 .

3. Evaluate the quality Q 0 of the set C 0 of elements T ∈ C 0 by a given function, e.g. 10.

4. For each n ∈ N 0 do:

• Connect the node n to each edge e ∈ ∂E 0 , n ∈ e to get a set of elements C n = {T e , T e = (e, n), attempting to re triangulate the cavity C 0 }.

• Evaluate the quality Q n of the set C n by a given function, e.g. 10.

• If Q 0 < Q n (in the sense of Definition 4.2) then set C 0 ← C n , Q 0 ← Q n .
Until stagnation of the changes applied to the mesh.

Modified Quality function

In our application, we have found useful to re-define slightly some component of the MTC algorithm. This section details these modifications. First, we have re-defined the functions Q S (T ) and Q F (T ) in expression [START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent cfd[END_REF] in the following way. We take

Q S (T ) = min (i,j)∈T h M (T ) ij , 1 h M (T ) ij d , instead of Q orig S (T ) = min h M (T ) , 1 h M (T ) d ,
This new definition avoid taking averages, because the mesh lengths may be very different on an element due to the anisotropy of the mesh.

The measure of the shape quality remains the same as in the original MTC code,

Q F (T ) = C 0 (V olume) M (T ) (T ) h d M (T )
, Finally, we have found useful to add to the measure of the quality of an element [START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent cfd[END_REF] an additional function Q V (T ) to have a better control on the element characteristic volume. Thus our definition of the quality of an element is now

Q(T ) = Q S (T ) • Q F (T ) • Q V (T )
This increase the sensitivity of Q(T ) to the non-respect of prescribed aspect ratio. The definition of Q V (T ) that we use is

Q V (T ) = min 1 (V olume) M (T ) (T )
, (V olume) M (T ) (T ) .

APPLICATIONS

Coarsening synthetic mesh

The coarsening process and development detailed in the previous sections, is applied here on a model problem in order to show how it works. Then in the following section, we will apply the semi-coarsening tool to industrial meshes with a very large aspect ratio.

Our model problem is a regular hexaedral domain Ω divided into 8 hexaedrons whose volumes are V = δx×δy ×δz, we take δx = δy and δz << δx.

In the example given below, δz = 10 -3 δx but we have checked that the results does not depend on this ratio and that the results are unchanged even with δz = 10 -6 δx. Each hexaedron is divided into two prisms, each being subdivided into three tetrahedra. Figure 7 (a) shows this tetrahedrization that contains 27 nodes and is composed of three node layers in the z direction. On each of these 27 nodes, the eigenvectors of the corresponding metric are approximately the three x, y, z axis with corresponding eigenvalues respectively equal to (1/δx) 2 , (1/δy) 2 , (1/δz) 2 . The coarsening algorithm will change this metric for another one with the same eigenvectors but with eigenvalues corresponding to the length δx, δy, 2δz then the mesh generator will try to construct a mesh corresponding to this wished metric. The result is shown on figure [7 (b)]. As expected, this mesh contains 18 nodes and the coarsening process has eliminated the intermediate node layer thus performing a semi-coarsening by a factor 2 in the z direction.

Conclusion

The present work has studied the capability of a semi-coarsening algorithm to handle extremely stretched meshes used for Navier-Stokes computations with low Reynolds turbulence modeling. This algorithm is able to generate a sequence of meshes with decreasing number of nodes in the direction perpendicular to the geometry. This sequence of meshes is thus suitable for multigrid acceleration.

In particular in this work we have study the interpolation procedures and smoothing of the metric fields. Our work has relied on some recent results on metric interpolations [START_REF] Arsigny | Fast and Simple Computations on Tensors with Log-Euclidean Metrics[END_REF][START_REF] Fletcher | Principal Geodesic Analysis on Symmetric Spaces: Statistics of Diffusion Tensors[END_REF] as well as on some more simpler ideas based on averaging of the mesh edges length. In addition, we have proposed a way to handle the connection between isotropic regions in the far field and the anisotropic ones in the boundary layers. Then we have designed an anisotropic coarsening algorithm able to identify the directions of stretching and to remove nodes in the direction of minimal mesh size only. We have also performed several relevant tests for the semi-coarsening strategy. These experiments show that geometrical methods can now be used for multigrid acceleration even for highly anisotropic meshes employed for Navier-Stokes computations with turbulence modeling in an industrial context.

Acknowledgements

This work has been supported in part by a grant from TO BE COM-PLETED

Figure 1 : 2 -

 12 Figure 1: 2-D model boundary layer mesh

Figure 2 :

 2 Figure 2: Semi-coarsening of a model boundary layer mesh. Note that the y-axis is scaled in order to be able to visualize the result.