T Balarac

R Berthelon

M Mercier

P Bernard

G Bénard

O Lartigue

Métais

A Grenouilloux

J Leparoux

V Moureau

G Balarac

T Berthelon

R Mercier

M Bernard

P Bénard

G Lartigue

O Métais

Toward the use of LES for industrial complex geometries. Part I: Automatic mesh definition

Keywords: Large-Eddy Simulation, Mesh Adaptation, Automatic mesh generation

Toward the use of LES for industrial complex geometries. Part I: automatic mesh definition

Introduction

Owing to the large range of scales in high Reynolds number turbulent flows, the Direct Numerical Simulation (DNS) of realistic applications is not yet available because of the significant computational cost. To overcome this limitation, the Large Eddy Simulation (LES) technique proposes to explicitly solve the large scales of the flow and to model the effect of the smallest ones thanks to a sub-grid scale model (SGS). For the past years, the strong increase in computational power has allowed the numerical study of industrial configurations with this methodology [START_REF] Moureau | From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flamepdf modeling[END_REF][START_REF] Grenson | Large-eddy simulation of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF][START_REF] Papadogiannis | Large eddy simulation of a high pressure turbine stage: Effects of sub-grid scale modeling and mesh resolution[END_REF][START_REF] Benard | Large-eddy simulation of the lean-premixed preccinsta burner with wall heat loss. Proceedings of the Combustion Institute[END_REF][START_REF] Arroyo | Towards the large-eddy simulation of a full engine: Integration of a 360 azimuthal degrees fan, compressor and combustion chamber. part i: Methodology and initialisation[END_REF]. To gain a deeper insight on state of the art applications of LES as well as its long term vision, the reader may refer to the latest reviews of Moreau et al [START_REF] Moreau | The third golden age of aeroacoustics[END_REF] or NASA's CFD vision for 2030 [START_REF] Slotnick | Cfd vision 2030 study: A path to revolutionary computational aerosciences[END_REF]. LES relies on an explicit filtering of the Navier-Stokes equations. The filter width ∆ defines the cut-off frequency of the scale separation. In practice, as it is the case of most CFD codes, the filtering operation is rather implicitly carried out by the mesh and is related to the local cell-size. Therefore, the quality of the computed solution is directly linked to the quality of the mesh. However, the definition of an adequate mesh is still an issue when LES is used to predict the flow in an industrial complex geometry.

Two main types of meshes are used in LES. On the one hand, structured meshes are considered in a wide range of applications [START_REF] Blazek | Chapter 11 -principles of grid generation[END_REF][START_REF] Stich | Jet noise prediction for chevron nozzle flows with wall-modeled large-eddy simulation[END_REF][START_REF] Renaud | Computation of goahead configuration with chimera assembly[END_REF] and use a set of quadrilateral (2D) or hexahedral (3D) elements. While ensuring a lighter mesh, lines must be extended to the limits of the domain. This means that a particular discretization of a given region forces the presence of cells in other parts of the domain when not required [START_REF] Grenson | Large-eddy simulation of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF][START_REF] Luo | Numerical simulation of turbulent flow over a backward facing step using partially averaged navier-stokes method[END_REF]. In most practical cases, a conjunction of several mesh-blocks of different characteristics is used, but a proper transition is not always guaranteed and often leads to hanging nodes or lines. On the other hand unstructured meshes use a wide variety of shapes, ranging from polygonal structures of Voronoï meshes to combination of various elements of different nature, as prism layer and simplicial elements such as tetrahedra. They have proven to be particularly suited for the meshing of complex topologies as in industrial cases.

Authors have established a series of guidelines for the correct mesh generation [START_REF] Piomelli | Large-eddy simulations: Theory and applications[END_REF], or often perform the same numerical simulation with different mesh to quantify the sensitivity of the solution [START_REF] Hadžiabdić | Vortical structures and heat transfer in a round impinging jet[END_REF]. Nevertheless, in the LES framework an improvement of the mesh resolution will decrease both the overall numerical discretization error and the effects of the SGS model, as the cut-off frequency is shifted, making it particularly challenging to separate the effect of these two contributions. For this reason, the concept of mesh independency can not be achieved in LES, as the solution evolves until DNS is eventually reached. In an effort to estimate the influence of the mesh quality on the computed solution, authors have proposed different approaches. The so-called systematic grid-variation methods require the use of several geometrically similar meshes but of different resolutions. Klein et al. [START_REF] Klein | An attempt to assess the quality of large eddy simulations in the context of implicit filtering[END_REF] proposed to verify the quality of the solution thanks to a Richardson's extrapolation of the error. By comparing the solution on a coarse, intermediate and a refined mesh, Klein et al. were able to define a threshold for the acceptance of a solution. Richardson's extrapolation has also been proposed by Celik et al. [START_REF] Celik | Index of Resolution Quality for Large Eddy Simulations[END_REF] to formulate an Index of Quality, often referred as LES-IQ. Both methodologies require the use of two to three meshes, thus increasing the overall cost of the numerical run. Single grid estimators methods aim at quantifying errors on a unique mesh and therefore appear as more affordable. Geurts and Frohlich [START_REF] Geurts | A framework for predicting accuracy limitations in large-eddy simulation[END_REF] presented an indicator of the subgrid-scale activity 0 ≤ s < 1 defined as the ratio of the volume averaged turbulent dissipation over the volume averaged overall dissipation (ie turbulent and molecular). Values of s = 0 correspond to a simulation that can be qualified as DNS and s = 1 to LES. Yet, Celik et al. [START_REF] Celik | Index of Resolution Quality for Large Eddy Simulations[END_REF] reported that s is very little sensitive to the mesh resolution and might not be the most suited quantity to use. Pope [START_REF] Pope | Turbulent flows[END_REF] proposed to consider the ratio between the resolved and the total turbulent kinetic energy, and stated that well-resolved LES should resolve at least 80% of the total turbulent kinetic energy. These two approaches, namely systematic grid-variation and single grid estimators, present themselves as a posteriori indicators for the quality of the mesh used.

Mesh adaptation appears as one of the most promising techniques for optimized mesh generation [START_REF] Park | Unstructured grid adaptation: Status, potential impacts, and recommended investments toward cfd vision[END_REF]. While extensive work has been done for Reynolds-Averaged Navier-Stokes simulations (RANS) [START_REF] Michal | Comparing anisotropic error estimates for the onera m6 wing rans simulations[END_REF][START_REF] Alauzet | 3d rans anisotropic mesh adaptation on the high-lift version of nasa's common research model (hl-crm)[END_REF][START_REF] Alauzet | Feature-based and goal-oriented anisotropic mesh adaptation for rans applications in aeronautics and aerospace[END_REF][START_REF] Alauzet | Periodic adjoints and anisotropic mesh adaptation in rotating frame for high-fidelity rans turbomachinery applications[END_REF][START_REF] Frazza | Unstructured anisotropic mesh adaptation for 3d rans turbomachinery applications[END_REF], the unsteady nature of the LES equations challenges the generation of the optimal mesh. Two main feature-based approaches have been developed, namely Dynamic or Static Mesh Refinement (hereinafter called DMR and SMR). DMR is applied to unsteady cases and allows to carefully refine regions with strong moving fronts, such as in two-phase flows [START_REF] Pertant | A finite-volume method for simulating contact lines on unstructured meshes in a conservative level-set framework[END_REF][START_REF] Janodet | A massively parallel accurate conservative level set algorithm for simulating turbulent atomization on adaptive unstructured grids[END_REF] or reactive cases [START_REF] Moureau | Parallel and dynamic mesh adaptation of tetrahedral-based meshes for propagating fronts and interfaces: application to premixed combustion[END_REF]. SMR, also known as the h-refinement method, is often used to improve the mesh in regions where time-averaged Quantities of Interest (QoI) show strong gradients. Bénard et al. [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF] proposed to improve the discretization of a meso-combustor chamber based on a combination of two criteria: an interpolation error based on the mean velocity and an error analogous to Pope's criterion [START_REF] Pope | Turbulent flows[END_REF] to ensure the resolution of 80% of the turbulent kinetic energy spectrum. Daviller et al. [START_REF] Daviller | A mesh adaptation strategy to predict pressure losses in les of swirled flows[END_REF] proposed to target the dissipation of the kinetic energy to match the correct pressure loss in an injector. This work has recently been extended to reactive cases by Agostinelli et al. [START_REF] Agostinelli | Static mesh adaptation for reliable large eddy simulation of turbulent reacting flows[END_REF], where two methodologies for the optimisation of an exiting mesh or the generation of the optimal mesh from an arbitrary initial mesh are presented. Yet, this iterative mesh adaptation process relies on an a posteriori comparison with experimental data to define a stop criterion. Therefore, such a stopping criterion is hard to establish when the user has a priori no knowledge of the flow characteristics.

Another important key ingredient for automatic mesh determination is the mesh generation itself. In many industrial cases, this latter is still a tedious task and a bottleneck. The difficulties arise often from the CAD quality required to perform the mesh generation. Most meshing algorithms are not robust to holes in the geometry, small surface mismatches or intersecting surfaces. For this reason, automatic mesh generation often starts with CAD cleaning and repairing [START_REF] Ito | Unstructured mesh generation for viscous flow computations[END_REF], which hinders the whole process. Implicit meshing based on level set functions [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF] has gained momentum over the last decade [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF]. In this method, a body-fitted mesh is created from an existing mesh and an implicit surface representation. The boundary condition is created at the location of an iso-contour defined on the existing mesh, which is subsequently adapted to keep a good element quality. This type of method is very well suited for an automatic mesh determination process as it is very robust. However, determining the level set function from an existing CAD is difficult, especially in a parallel computing context.

A methodology for the automatic mesh determination in complex industrial geometries is then proposed in this work. This methodology does not require any a priori flow knowledge. Such a process relies on two main steps. First, an Eulerian unstructured mesh is generated from a set of STL files. Such a format is commonly accepted by most CAD software. The proposed method is based on level set functions and implicit meshing. Then, an iterative mesh convergence process relies on an extension of the work initiated by Bénard et al. [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF] and aims at generating an optimal mesh guaranteeing a mesh independent mean kinetic energy (MKE) balance. Such a methodology, presented in this paper, referred as Part I, is an important step to extend the use of LES as a decision support tool. Another important challenge is to be able to reduce the LES "time-to-solution". This is the objective of the second paper, referred as Part II, which extends the automatic mesh determination with an automatic time-step determination by using a linearized implicit time advancement.

Numerical tool: LES solver based on unstuctured body-fitted mesh

In this work, LES were performed using the YALES2 flow solver [START_REF] Moureau | Design of a massively parallel cfd code for complex geometries[END_REF]. This code solves the incompressible and low-Mach number Navier-Stokes equations for turbulent flows on unstructured meshes using a projection method for pressure-velocity coupling [START_REF] Chorin | Numerical solution of the navier-stokes equations[END_REF]. It relies on central finite-volume schemes and on highly efficient linear solvers [START_REF] Malandain | Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines[END_REF]. The time integration is explicit for convective terms using a fourth-order modified Runge-Kutta scheme [START_REF] Kraushaar | Application of the compressible and low-mach number approaches to largeeddy simulation of turbulent flows in aero-engines[END_REF], with a semi-implicit integration for the diffusive terms. The YALES2 solver is able to deal with unstructured meshes composed of tetrahedral elements allowing to perform LES or DNS of complex geometries in the context of massively parallel computations. To manage mesh adaptation, the re-meshing step is based on MMG, a sequential isotropic and anisotropic mesh adaptation library for tetrahedral (3D) and triangle elements (2D) [START_REF] Dobrzynski | Anisotropic delaunay mesh adaptation for unsteady simulations[END_REF] by using a parallel mesh adaptation strategy [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF]. The solver has been validated for various applications such as combustion [START_REF] Moureau | From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flamepdf modeling[END_REF][START_REF] Guedot | Design of implicit high-order filters on unstructured grids for the identification of large-scale features in large-eddy simulation and application to a swirl burner[END_REF], bio-mechanics [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF], hydro-electricity [START_REF] Wilhelm | Analysis of head losses in a turbine draft tube by means of 3d unsteady simulations[END_REF][START_REF] Véras | Reconstruction of proper numerical inlet boundary conditions for draft tube flow simulations using machine learning[END_REF], or wind energy [START_REF] Benard | Large-eddy simulation of wind turbines wakes including geometrical effects[END_REF]. For all the presented cases, a constant Courant-Friedrichs-Lewy (CFL) number is set to 0.9 and the dynamic Smagorinsky subgrid-scale model is used [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF].

Mesh adaptation strategy for LES based on mean kinetic energy equilibrium

Definition of a mesh-converged LES

As recalled in the introduction, in LES, the separation between resolved large scales and modeled small ones is performed by a filtering operation. This filtering operation is characterized by the filter width. In most CFD codes, the filtering is performed implicitly by the mesh and it is therefore stated that this width is directly linked to the local cell size. It is of paramount importance to define a mesh to guarantee both the correct discretization of the flow and that a sufficient part of turbulent scales are explicitly solved for reliable LES.

Considering the velocity field, ⃗ u, based on Reynolds decomposition ⃗ u = ⟨⃗ u⟩ + ⃗ u ′ , the mean kinetic energy, K = 1/2⟨u 2 i ⟩ can be decomposed as the kinetic energy of the mean field (MKE), K M = 1/2⟨u i ⟩ 2 , and the turbulent kinetic energy (TKE),

K T = 1/2⟨u ′2 i ⟩.
Here the bracket symbol stands for the statistical average over a suitable ensemble. In the LES framework, ergodicity of the flow is often assumed and statistics are computed thanks to a time-averaging procedure. This decomposition leads then to K = K M +K T . In LES, considering now that a part of the fluctuating field is unresolved, the TKE can be decomposed as the resolved part, K > T = 1/2⟨ū ′2 i ⟩ and the unresolved subgrid-scale (SGS) part, K < T = K T -K > T , with the resolved velocity fluctuation ⃗ ū′ = ⃗ ū -⟨ ⃗ ū⟩ and ⃗ ū the resolved velocity field. Here the bar sign stands for the classical LES filtering operator. The part of resolved and SGS (unresolved) TKE are directly linked to the mesh filter size. To be relevant, the LES technique assumes that the filtering operator acts only at turbulent scales, and hence that the mean field is not affected, i.e. ⟨ ⃗ ū⟩ = ⟨⃗ u⟩. The kinetic energy decomposition in LES then reads

K = K M + K > T + K < T .
The proposed mesh adaptation strategy relies on this natural decomposition of the mean kinetic energy field, proper to the LES framework. Indeed, it is based on the premise that accurate LES should, at least, lead to mesh-independent mean fields. In particular, the kinetic energy of the mean field (MKE) should be mesh-independent. This will be obtained if the MKE balance is mesh-independent [START_REF] Grenouilloux | On the use of kinetic-energy balance for the feature-based mesh adaptation applied to large-eddy simulation in complex geometries[END_REF]. With LES formalism and assuming a classic SGS modeling strategy based on an eddy viscosity, ν T , the MKE balance writes

∂K M ∂t = ⟨ū ′ i ū′ j ⟩⟨ Sij ⟩ -2⟨ν T Sij ⟩⟨ Sij ⟩ Transfer to TKE -2ν⟨ Sij ⟩⟨ Sij ⟩ Dissipation + ∂ ∂x j -⟨ū j ⟩K M + ⟨ū j p * ⟩ + ⟨ū i ⟩ 2ν⟨ Sij ⟩ + 2⟨ν T Sij ⟩ + ⟨ū ′ i ū′ j ⟩ Diffusion , (1)
with p * the modified resolved pressure due to the SGS model. The energetic balance is then controlled by the transfer to TKE, i.e. the energy transfer between the mean and turbulent fields, the molecular dissipation, i.e. the pure loss of kinetic energy due to molecular viscosity effect, and diffusion, i.e. the spatial energy redistribution. Note that this last term does not participate to the global energy balance. Indeed, it cancels out when integrated over the whole domain.

The mesh adaptation strategy will be then defined to guarantee the global energy balance in a given computational domain, V, i.e. to guarantee that the global transfer to TKE,

P V TKE = V ⟨ū ′ i ū′ j ⟩⟨ Sij ⟩ -2⟨ν T Sij ⟩⟨ Sij ⟩dV, (2)
and the global molecular dissipation,

D V ν = -2ν V ⟨ Sij ⟩⟨ Sij ⟩dV, (3)
are mesh independent. As D V ν is only dependent on the mean field, this quantity is then a good candidate to assess the accuracy of the mean field discretization. P V TKE will allow to assess the relevance of the scales separation. Note that the two terms of P V TKE have to be considered together. Indeed, mesh refinement leads to a decrease of the filter size, and hence a decrease of the SGS part (the second term) and an increase of the resolved part (the first term). It is then not expected that the terms taken independently are mesh independent.

Metric definition

The objective of the proposed strategy is to start from an initial and arbitrary mesh, and to adapt it until P V TKE and D V ν are mesh independent. Knowing time-averaged quantities on the initial mesh, the objective is then to be able to define the new cell size of the adapted mesh, ∆ a (⃗ x), from the initial mesh size, ∆ i (⃗ x). Two different criteria have already been proposed by Bénard et al. [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF] and they are recalled here.

The first criterion QC 1,i (⃗ x) ensures a correct discretization of the mean field. It is defined from an error estimation of a given computed field f [START_REF] Habashi | Anisotropic mesh adaptation: towards user-independent, mesh-independent and solver-independent CFD. Part I: general principles[END_REF]. It can be shown that for a mesh composed by isotropic tetrahedra with typical size ∆ i (⃗ x), this error is limited by the quantity,

QC 1,i (⃗ x) = ∆ 2 i (⃗ x)H f (⃗ x), (4)
where H f is the highest absolute value of the eigenvalue of the Hessian of f . The adapted local mesh size is then defined based on a constant target value QC 1,T such as,

∆ a (⃗ x) = QC 1,T H f (⃗ x) , (5)
to control and make uniform the error in the overall computational domain. In this work, this criteria is applied to the kinetic energy of the mean field, K M , and the Hessian is computed using the high-order framework proposed by Bernard et al. [START_REF] Bernard | A framework to perform high-order deconvolution for finite-volume method on simplicial meshes[END_REF]. Note that the adequate value of QC 1,T is dependent on the case. For a given case, the strategy will be then to determine the adequate QC 1,T to guarantee a mesh independent D V ν field, that only depends of the mean field discretization.

The second criterion has to guarantee that enough turbulent scales are solved to perform reliable LES. For fully turbulent flows, this will be the case if at least 80% of the TKE is explicitly solved [START_REF] Pope | Turbulent flows[END_REF]. From an initial mesh, the part of unresolved TKE can be defined as

QC 2,i (⃗ x) = K < T (⃗ x) K T (⃗ x) = K < T (⃗ x) K > T (⃗ x) + K < T (⃗ x) , (6)
where K < T (⃗ x) is evaluated from the eddy viscosity [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF]. Assuming that the eddy viscosity varies as the filter size with a 4/3 power law [START_REF] Lesieur | Turbulence in fluids[END_REF], it can be shown that the adapted local mesh size should be

∆ a (⃗ x) = ∆ i (⃗ x) QC 2,T QC 2,i 3/2 , (7)
where QC 2,T is a constant target value. For fully turbulent flows, QC 2,T = 0.2 [START_REF] Pope | Turbulent flows[END_REF]. Yet, this value might be case-dependent for more complex configurations. The methodology will then determine the adequate QC 2,T to guarantee a mesh independent P V TKE field, that includes the effect of the modeled scales on the MKE budget. Note that the application of both criteria will lead to two different metrics for the adaptation process, but only the smallest one will be considered.

Automatic mesh convergence procedure

The mesh convergence procedure proposes to sequentially determine the target values for QC 1,T and QC 2,T , respectively. This procedure does not require any a priori flow knowledge or any comparison to reference results, which is mandatory for the use of LES for complex flows as a predictive tool. Note also that the procedure has been fully automated and it does not require any user intervention. In particular, as the mesh is adapted there is a rescaling of the number of CPU cores based on the mesh size to ensure the code performances even when massively distributing the mesh on a super-computer. Then, the number of CPU cores is set to roughly obtain 50,000 elements per core [START_REF] Gava | Optimisation of the parallel performances of a cfd solver for emerging computational platforms[END_REF].

A sketch of the overall algorithm is shown in Fig. 1. Starting from an initial and arbitrary mesh, flow statistics are computed (State = 0) once the initial, and sometimes unphysical, transient is evacuated. Keeping the theoretical value QC 2,T = QC 0 2,T = 0.2, the optimal value QC 1,T is found by an iterative process. The initial QC 1,T = QC 0 1,T value that allows to homogenize the MKE field discretization error by conserving the number of cells of the initial mesh [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF] is chosen. The initial mesh is adapted given this initial value QC 0 1,T . Flow statistics are computed on the new mesh, and a new adaptation is performed keeping QC 0 2,T and QC 0 1,T . Indeed, because QC 1,i (⃗ x) and QC 2,i (⃗ x) are defined from flow quantities, an iterative procedure is needed to converge the mesh at fixed values of QC 1,T and QC 2,T . The procedure is then repeated until the number of cells between two succesive adaptations varies by less than 10%. This allows a first estimation for the global molecular dissipation, D V,0 ν , given the first targeted value QC 0 1,T . Then a decreased value QC 1,T = QC 1 1,T = QC 0 1,T /1.6 is considered. Note that in 3D, if the adaptation is purely based on the criterion QC 1 , dividing QC 1,T by 1.6 = 2 2/3 allows to compute a refinement ratio that multiplies by 2 the number of cells in all the computational domain [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF]. The mesh is converged for this new value, i.e. it is adapted until the overall cell counts varies less than 10%. Then an estimation of the global molecular dissipation, D V,1 ν given the target value QC 1 1,T is obtained. The procedure is repeated k times, until ∆D

V ν = |D V,k ν -D V,k-1 ν |/D V,k ν < 0.05. The final value for QC 1,T is fixed equal to QC k 1,
T . This part is identified as State = 1 on Fig. 1. At the end of this State, the mesh has been adapted in such a way that viscous dissipation is now mesh-independent.

Keeping this final value for QC 1,T , a reduced value QC 2,T = QC 1 2,T = QC 0 2,T /1.17 is tested. As for the first criterion, the division by 1.17 = 2 2/9 allows to compute a refinement ratio that multiplies by 2 the number of cells in all the computational domain if the adaptation is purely controlled by the second criterion [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF]. Mesh is converged for these given values, and the global transfer to TKE, P V,1 TKE is obtained given QC 1 2,T . The procedure is repeated l times, until ∆P V TKE = |P V,l TKE -P V,l-1 TKE |/P V,l TKE < 0.05 and ∆D V ν < 0.05. Then, the final value for QC 2,T is fixed equal to QC l 2,T . This part is identified as State = 2 in Fig. 1. The final mesh with both a mesh independent transfer to TKE and molecular dissipation of the MKE is obtained. At last, main statistics of interest are computed (State = 3).

Note that before each adaptation step, statistics need to be computed. In this work, the time window for statistical accumulation is fixed, in a compromise between the "time-to-solution" of the simulation and the quality of statistics convergence. Future works will be devoted to define an automatic procedure to define the optimal time windows for statistical accumulation.

Canonical test cases

To assess the performance of the described automatic mesh convergence procedure two canonical test cases are first considered. The first flow configuration consists of a round jet in transition to turbulence, whereas the second configuration is a turbulent pipe flow.

Turbulent round jet

The round jet flow configuration is defined by its inlet velocity profile. At the inlet, the mean velocity field is non zero only for the streamwise component, which is given by a hyperbolic tangent profile [START_REF] Michalke | On the inviscid instability of a circular jet with external flow[END_REF],

u ref (x) = U 1 + U 2 2 - U 1 -U 2 2 tanh R 4Γ 0 r R - R r ,
where U 1 is the centerline velocity, U 2 is a small co-flow, Γ 0 is the momentum thickness of the initial shear layer, r the radial coordinates (taking the origin at the center of the jet), and R the initial jet radius. The Reynolds number is fixed at a moderate value, Re = (U 1 -U 2 )R/ν = 4, 000. To accelerate the transition two forcing terms are added. First, an isotropic turbulence is injected at the inlet in the shear layer regions with an amplitude set to 3% of U 1 -U 2 [START_REF] Guedot | Design of implicit high-order filters on unstructured grids for the identification of large-scale features in large-eddy simulation and application to a swirl burner[END_REF][START_REF] Santoso | Hybrid particle-grid methods for the study of differential diffusion in turbulent flows[END_REF]. For the streamwise component, the forcing is then complemented by a deterministic part, which consists in a varicose (axisymetric) excitation [START_REF] Da Silva | Vortex control of bifurcating jets: A numerical study[END_REF],

u df (x, t) = ϵu ref (x) sin 2πSt R U 1 + U 2 2Γ 0 t ,
with a forcing amplitude ϵ = 2.5% , and St R fixed at 0.033 to trigger the frequency predicted by the linear stability theory [START_REF] Michalke | On the inviscid instability of a circular jet with external flow[END_REF]. The computational domain is a parallelepiped with size 30R × 28R × 28R in the streamwise and the two transverse directions, respectively. The size of the transverse directions is large enough to allow the flow rate added by the co-flow to be greater than the flow rate entrainment of the jet. This initial mesh has been coarsen on purpose. It is generated as a cartesian uniform structured mesh where each hexahedral cell is cut into 6 tetrahedral cells. The edges of the hexahedral cells are equal to R, the jet's radius, meaning that only two cells discretized the initial diameter of the jet. The initial tetrahedral mesh is composed of 307, 806 elements and is shown in Fig. 2 (left-top). The jet's dynamics predicted with this initial mesh (Fig. 2, left-bottom) do not reproduce the well-known scenario of transition in free shear layers. Indeed, even with the forcing, the Kelvin-Helmholtz instability does not properly arise, preventing the formation of Kelvin-Helmholtz vortices. Indeed, this instability depends on the thickness of the shear layer, and the initial coarse mesh does not allow to describe it properly. The proposed automatic mesh convergence procedure is then performed starting with this very coarse mesh, as described by the sketch in Fig. 1. During the procedure, the statistical time window has been fixed to T w = 170R/(U 1 -U 2 ), which appears as a good compromise between statistics quality and the "time-to-solution" of the simulation. At the end of the procedure the final statistics have been computed on a time window equal to 600R/(U 1 -U 2 ). ν (top), the global transfer to TKE, P V TKE (middle), and the number of elements, N elem (bottom) during the different states of the automatic mesh convergence procedure in the turbulent round jet configuration. Note that the resolved, P V res , and the modeled P V mod , parts of the transfer to TKE are also shown.

The evolution of the main quantities characterizing the procedure are shown by Fig. 3. Based on the statistics on the initial mesh (State 0), the sequential determination of the target values for QC 1,T (State 1) and QC 2,T (State 2) starts. The number of elements has a strong increase at the first steps, leading to a final mesh composed of 8, 125, 037 elements. The global molecular dissipation, D V ν , and the global transfer to TKE, P V TKE also have a strong variation during the first steps. It can also be observed that the resolved part of the transfer to TKE increases and the modeled part of the transfer to TKE decreases at the end of the mesh convergence procedure as expected with the refinement of the mesh. Note that the contribution of the transfer to TKE dominates the contribution of the molecular dissipation in the global MKE balance, and the modeled part of the transfer to TKE is weak for this moderate Reynolds number case. The final mesh and the flow dynamics predicted with it are shown in Fig. 2 (right). As expected, at the on-start of the jet, the mesh elements are concentrated in the shear layer. This allows to properly predict the Kelvin-Helmholtz vortices generation down to the end of the potential core. After, in agreement with the classical scenario of transition in free shear layers, there is an abrupt increase in the level of small-scales turbulence. The jet's width increases due to the turbulent diffusion of the momentum. This leads to a radial expansion of the refined regions in the second part of the computational domain.

Finally, the mesh influence on the statistics prediction is shown in Fig. 4. The mean and root mean square (rms) axial velocities obtained with the initial and the final meshes are compared with a reference case. This latter is based on an homogeneous mesh refinement [START_REF] Rivara | Mesh refinement processes based on the generalized bisection of simplices[END_REF] of the final mesh obtained with the proposed methodology, thus leading to a mesh of 65, 020, 936 elements. The strong numerical diffusion due to the very coarse mesh leads to poor prediction even at the on-start of the jet. Conversely, the final mesh is in good agreement for both mean and rms velocity profiles at the different stages of the round jet transition. This confirms that the proposed automatic mesh convergence procedure built to guarantee a mesh-independent MKE balance is able to converge towards a mesh that correctly predicts shear-layer transition even starting from a very coarse mesh. To illustrate the ability of the proposed mesh convergence procedure to lead to accurate LES, this procedure has been applied by starting from another initial mesh (Appendix A), and also by using the LIKE criterion [START_REF] Daviller | A mesh adaptation strategy to predict pressure losses in les of swirled flows[END_REF] to define the adapted mesh instead of the QC 1 and QC 2 criteria (Appendix B). In both cases, good statistics of the turbulent jet are obtained at the end of the mesh convergence procedure.

Turbulent pipe flow

The proposed automatic mesh convergence procedure is now applied to a turbulent wall-bounded flow. A turbulent pipe flow with a Reynolds number based on the bulk velocity, U b , and the pipe diameter, D, equal to Re = U b D/ν = 5, 300 is considered. The flow configuration is similar to the DNS performed by Wu and Moin [START_REF] Wu | A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow[END_REF] on a structured mesh composed of 67.7 millions cells. This DNS results will be used as reference. A periodic condition is used between the inlet and the outlet, with a forcing to ensure a constant bulk velocity. The length of the simulated pipe is 15R as suggested by Wu and Moin [START_REF] Wu | A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow[END_REF]. The initial mesh is composed by tetrahedral cells with edge sizes increasing linearly with a growth factor equal to 1.1 from 0.08R at the wall, to 0.15R at the center of the pipe, where R is the pipe's radius. This initial tetrahedral mesh is composed of 498, 798 elements. This mesh and the dynamics of the flow predicted with it are shown in Fig. 5

(left).

Figure 6 shows the evolution of quantities during the mesh convergence procedure. Here, the time window for statistical computation is fixed to T w = 75D/U b and the final statistics are computed on a time window equal to 600D/U b . As for the turbulent round jet case, strong mesh variations occur during the first mesh adaptation steps. The initial mesh is so coarse that the part of modeled transfer to TKE exceeds the part of resolved transfer to TKE. This implies that the transfer to TKE dominates the molecular dissipation. The modeled part of the transfer to TKE strongly decreases during the mesh convergence procedure to be a less than 10% of the overall transfer to TKE. This leads to a global reduction of the transfer to TKE. Finally the molecular dissipation caused by the friction at the wall dominates the MKE energy balance. The final mesh is composed of 58, 556, 190 elements. The final mesh and the dynamics of the flow predicted with it are shown by Fig. 5 (right). As expected, the elements are mainly gathered close to the wall, leading to the prediction of a thinner boundary layer with higher turbulent activities. Note that the number of elements of the final mesh is close to the number of cells used in the reference case [START_REF] Wu | A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow[END_REF]. This is because a structured mesh with anisotropic elements close to the wall is used in the DNS simulation. The remeshing tool used in this work manages only isotropic tetrahedral elements. Future work will be devoted to the extension to anisotropic elements to avoid too strong constraint close to walls for example.

Statistical quantities obtained with the initial and final meshes are now compared with the DNS results in Fig. 7. As expected, the near wall behavior is not correctly predicted with the initial mesh. Due to a too coarse mesh close to the wall, the viscous sub-layer is over-predicted leading to a shorter log-law region. This is also visible on rms quantities, where the peak of turbulent intensity is predicted further from the wall in comparison with the DNS. The final mesh corrects this behavior. A good agreement is observed for the mean velocity profile, and the peak of turbulent intensity is correctly predicted. Note that when rms quantities are considered, an exact agreement between DNS and LES results is not expected because a part of the turbulent field is not resolved in LES. This second test case confirms that the proposed automatic mesh convergence procedure is able to define a correct mesh also for wall-bounded flows.

Complete mesh determination for the LES of the PRECCINSTA burner

The LES of the PRECCINSTA swirl burner [START_REF] Meier | Detailed characterization of the dynamics of thermoacoustic pulsations in a lean premixed swirl flame[END_REF] is now considered as a representative configuration of LES in a complex geometry. This configuration is representative of an aeronautical injection system developed to stabilize aerodynamically a swirl flame in a helicopter combustion chamber. The computational domain comprises a plenum, an air swirler, a square combustion chamber (86 × 86 × 110mm 3 ) and an exhaust pipe. This configuration is well-known for being a semi-industrial test case with advanced experimental characterization on several non-reacting and reacting operating conditions. It has been extensively used for combustion model validation during the last 18 years [START_REF] Moureau | From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flamepdf modeling[END_REF][START_REF] Lartigue | Experimental and numerical investigation of self-excited combustion oscillations in a scaled gas turbine combustor[END_REF][START_REF] Mercier | Les of turbulent combustion: On the consistency between flame and flow filter scales[END_REF]. The non-reacting operating condition is considered in this study with an air flow rate set to 1.114 × 10 -2 m 3 /s. The complete mesh determination strategy is composed by first an automatic mesh generation from STL files followed by the automatic mesh adaptation strategy presented above.

Mesh generation from STL files

Constraints

The aim of this work is to provide high-quality adapted meshes of complex geometries in order to perform high-fidelity LES. One tedious task in the simulation workflow is the cleaning of the geometry before meshing. Here, the geometry is assumed to be a triangulation stored in STL format for example and obtained from parametric surfaces or any other geometry representation. The STL format has gained a lot of momentum with additive manufacturing and is the chosen format here. The proposed meshing methodology aims at alleviating the cleaning step, which means that the methodology needs to be robust to low-quality triangulations with holes, intersecting triangles, infinitely thin walls, high-skewness or near-zero aera triangles, etc... Thus, the proposed methodology makes no assumption on the STL topology, i.e. if the triangles are connected or not, and on the shape or size of the triangles. The only quantity of interest that will be computed from the STL is the distance of points near the STL surface to the triangles, i.e. the distance of points to their projection onto the STL. As a consequence, the only requirement is to have an STL which correctly describes the geometry especially in regions of high curvature or corners and feature edges. The methodology is illustrated by generating a body-fitted mesh for the PRECCIN-STA burner. The starting STL is plotted in Fig. 8. It consists of approximately 100,000 triangles and 50,000 nodes. The triangle size and quality, which can be measured by the skewness of the triangles, are highly variable especially in the swirler, i.e. in the channels of the injector which make the fluid rotate. A cut of the final generated mesh is given in Fig. 9, in which the plenum on the left-hand side, the swirler in the center and the chamber on the right-hand side are clearly visible. Fig. 9 illustrates also the strong cell size gradient at the wall, which is imposed to have a good description of the geometry.

Methodology

Step 1: Loading of the uniform volume mesh and of the STL The first step consists in loading the uniform Eulerian mesh and of the STL file in parallel. Each processor loads a fraction of the STL files so that the number of triangles is evenly distributed on the processors. The triangles are then relocalized on the mesh: each triangle is moved and associated to the mesh element which contains the center of the triangle. This operation may require to send the triangle to another processor. Step 2: Creation of particles from the STL The loaded triangles are then subdivided into four smaller triangles recursively until each edge of the triangle is smaller than the local cell size and relocalized. The aim is to obtain a set of triangles, which can be transformed into point particles located at the barycenter of each triangle as shown in Fig. 10. The generated particles are clearly visible and their arrangement illustrates the subdivision process. The motivation of using particles instead of keeping the subdivided triangles is that the particles can be easily localized in the elements of the mesh and data structures for Lagrangian particles available in many CFD codes can be reused.

Step 3: Volume refinement of the mesh The created particles carry data such as the coordinates of the nodes of the associated subdivided triangle but also the desired cell size on the STL. This information enables to refine the Eulerian mesh as depicted in Fig 11. This refinement process relies on several kernels to ensure robustness, speed and quality of the mesh. First, starting from the mesh in Fig. 11a, the cells which contain a particle with a desired cell size smaller than the current local cell size are flagged and their edges are cut recursively in half as shown in Fig. 11b. Second, a prescribed cell size gradient is imposed by cutting recursively the edges of the elements around the STL as presented in Fig. 11c. Third, the MMG library is used to remesh and ensure a good quality of the elements as shown in Fig. 11d. The quality is quantified here by computing the skewness of the tetrahedral elements. If required, the particles are split into several sub-particles based on the subdivision of the carried triangle.

Step 4: Creation of a level set Once the Eulerian mesh is sufficiently refined, the distance d(x) to the particles or more specifically to the triangles carried by the particles can be computed. This distance function is not suitable for the identification of the inner and outer parts of the mesh and to localize the wall. To this aim, a new distance function d offset (x) is defined by substracting the local cell size: d offset (x) = d(x) -α∆x, where α takes values around 1.5. The iso-contour d offset (x) = 0 defines two surfaces on each side of the walls in most cases. If an interior point is prescribed, a single iso-contour can be identified and kept: d offset,inside (x) = 0. With this signed distance function, it is simple to move back the iso-contour close to the wall: d final (x) = d offset,inside (x) + β∆x. To avoid any crossing of the level set with the wall position, β is chosen smaller than α and some further corrections of the boundary position will be required. The final distance function and the iso-contour are illustrated in Fig. 12. Step 5: Cutting along the level set and removal of exterior cells From the level set, the edges of the Eulerian mesh, which cross the iso-contour, are cut and the outer cells are removed. The cutting kernel is depicted in Fig. 13. The slight offset of the level set from the wall is clearly visible.

Step 6: Projection onto the STL In order to correct the offset of the generated boundaries, a projection of the nodes of the new boundaries onto the STL triangles is performed. Nodes are moved if and only if the connected tetrahedral elements keep a positive volume. This step is performed several times until all the nodes have been projected onto the STL surface. This kernel is illustrated in Fig. 14 where the improvement of the surface mesh is noticeable. The projection step therefore guarantees that the surface nodes of the final body-fitted mesh are located onto triangles of the STL. However, the methodology does not guarantee that the sharp angles or feature edges of the STL mesh are represented in the final body-fitted mesh. The error in this representation depends on the prescribed cell size of the final body-fitted mesh. Features of the STL mesh smaller that the prescribed cell size cannot be represented. Step 7: Volume and surface adaptation before coarsening After projection, a volume and surface adaptation step is carried out with the MMG library to guarantee the quality of the mesh before further processing such as coarsening. This whole process is performed with very fine meshes at the walls to limit the error on the boundaries curvature and sharp edges. The resulting body-fitted mesh is presented in Figs. 15 and 16. In Fig. 15, the sliced mesh is colored by the edge length to show that the cell size and its gradient are well imposed. In Figs. [START_REF] Geurts | A framework for predicting accuracy limitations in large-eddy simulation[END_REF], the surface mesh is plotted along with the featured edges of the STL. A last adaptation is performed with the MMG library to coarsen the mesh far from these regions. The mesh has finally a minimum cell size of 0.8 mm to correctly describe the boundaries curvature and sharp angles, and the cell size growths linearly with a maximum cell size close to 30 mm in the combustion chamber. This allows to start with a mesh composed of only 446, 661 elements, but with a correct enough geometrical description of the boundaries (Figs. 9 and 17-left).

Mesh adaptation

The automatic mesh convergence procedure is then launched starting from the initial coarse mesh given by the procedure described in the previous section. Fig. 17 (left) shows the capture of the flow dynamics with this first mesh. This procedure is performed with a time window for statistical computation fixed to T w = 0.02s. The final statistics are computed on a time windows equal to 0.2s. The evolution of the main quantities characterizing the procedure are shown in Fig. 18. The evolutions of the quantities are similar to the turbulent jet and the pipe flow test cases. The stronger mesh variations occur during the first mesh adaptation steps. During the adaptation, the molecular dissipation increases due to a better capture of the strong velocity gradients (State 1), and the part of the modeled TKE transfert decreases due to a larger part of turbulent scales explicitly resolved (State 2). Due to the fully turbulent state of the flow, the global TKE transfert dominates the MKE balance in comparison with molecular dissipation contribution. The final mesh is composed of 178, 586, 310 elements (Fig. 17, right). To limit the computational cost, the plenum and the outlet regions are "masked", i.e. the initial mesh metric is conserved during the procedure, assuming that accurate simulation of these regions is not needed to obtain accurate results in the combustion chamber. The center planes colored by the instantaneous velocity norm are also shown for both the initial and the final meshes in Fig. 17 mesh. It is known that this swirling flow exhibits two distinct flow structures [START_REF] Bénard | Large-eddy simulation of the lean-premixed preccinsta burner with wall heat loss[END_REF] with a central and an outer recirculation zone, respectively. As expected, the mesh elements are mainly concentrated in these regions. Statistical quantities obtained with the initial and final meshes are now shown in Fig. 19. Axial mean and rms velocity profiles are compared at five different stations in the combustion chamber, from 1.5 mm to 35 mm downstream the swirler. These profiles are compared with experimental results [START_REF] Meier | Detailed characterization of the dynamics of thermoacoustic pulsations in a lean premixed swirl flame[END_REF][START_REF] Lartigue | Experimental and numerical investigation of self-excited combustion oscillations in a scaled gas turbine combustor[END_REF] and a reference simulation composed of 330 millions elements by using uniform mesh in the combustion chamber [START_REF] Moureau | From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flamepdf modeling[END_REF]. Regarding the mean axial velocity profiles, the simulation on the initial mesh completely fails to correctly predict the outer recirculation zone close the walls. The final mesh of the adaptation procedure leads to results close to the reference mesh composed by two times more elements. About the rms axial velocity profiles, the final adapted mesh allows to strongly improve the poor prediction obtain with the initial mesh. Note that rms values are slightly under-evaluated with the final adapted mesh in the center of the flow in comparison with the reference case. This is probably because a larger part of the turbulent field is resolved on the reference case due to smaller elements in this region. However, it is not necessary expected to recover all the rms fields in LES, since the LES strategy is to accept not to explicitly resolve all the TKE to reduce the computational cost.

Conclusion

In this paper, a user-independent strategy is proposed to define a mesh allowing to guarantee LES accuracy in complex geometries. This approach does not need any a priori knowledge of the flow dynamics. The strategy is composed of the unstructured mesh generation step and the automatic mesh convergence step. A robust and parallel mesh generation process has been designed based on i) a level set function computed from STL files, which describe the geometry, ii) implicit meshing, and iii) a projection of the boundaries onto the STL geometry. The automatic mesh convergence is based on two criteria to guarantee both the correct discretization of the mean field, and the determination of the cut-off of the resolved scales. These criteria are sequentially constrained until a mesh independent mean kinetic energy field is obtained since it is expected that only turbulent scales are unresolved in LES. This strategy has been successfully applied to three different configurations, namely a turbulent round jet, a turbulent pipe flow and finally the non-reactive PRECCINSTA burner. Designed to be user-independent, this approach opens the door to the use of LES approach as a decision support tool for various applications. Note that, considering the overall simulation process, the computational cost is distributed as 50% to 75% for the automatic mesh adaptation procedure, and 25% to 50% for the final statistics computed on the final mesh. To avoid too prohibitive "time-to-solution" and computational cost in LES approach, various extensions should be proposed with the optimization of the time windows of the statistical accumulation, or by using anisotropic cells to reduce the number of mesh elements. This can be also obtained using a larger time step value. This is the objective of a second paper, referred as Part II, which extend the automatic mesh determination with an automatic time-step determination by using a linearized implicit time advancement.

Appendix A. Turbulent round jet case starting from another initial mesh

To confirm that the procedure is little dependent of the flow computed on the initial mesh, the mesh convergence procedure has been applied to another initial mesh for the turbulent round jet case (section 3.4.1). This second mesh is composed of 9, 807, 676 elements, with a mesh refinement surrounding the jet (see Fig. A1-left). The automatic mesh convergence procedure is then started by using the same value to initialize QC 1,T . The evolution of main quantities during the overall procedure is shown by Fig. A2. The beginning of the State 1 leads first to a decrease of the number of elements due to the coarsening of the artificially refined regions. Then, the procedure is quasi similar to the procedure shown by Fig. 3 with a growth of the number of elements until a final mesh composed of 8, 387, 262 elements is obtained. Note that this number of cell counts is close to the final mesh presented in section 3.4.1. The similarities between the final meshes and the instantaneous velocity fields can be seen by comparing the Figs. A1 and 2 (right part). As expected, the statistics obtained on these two final meshes are in agreement, as shown by Fig. A3. This confirms that the proposed procedure is not very dependent of the initial mesh used. For a practical use, the strategy is to start the procedure with a mesh as coarse as possible, just refined to accurately represent the geometry details during the automatic mesh generation procedure (section 4.1). It is then expected that the automatic mesh convergence will allow to design the needed mesh in the core of the computational domain. criteria. One reason is that the LIKE criterion is not able to coarsen the mesh during the adaptation steps.

The mesh influence on the statistics prediction is shown in Fig. B2. The mean and root mean square (rms) axial velocities obtained with the final meshes are compared with the reference case. The final meshes obtained with both procedures are in good agreement with the reference case. For the procedure based on LIKE criterion, results on the second to last mesh are also considered. This mesh is composed by 13, 902, 402 elements. For this mesh, first and second order statistics are in good agreement except for the rms axial velocity at the beginning of the jet which are significantly underpredicted. This under-prediction is corrected by the last adaptation step. This confirms that the proposed mesh convergence procedure built to guarantee a mesh-independent MKE balance is able to lead to accurate LES for various mesh adaptation criteria.

3 Figure 1 .

 31 Figure1. Overview of the automatic mesh convergence procedure. ∆N CELL is the relative variation of the number of mesh cells before and after an adaptation.

Figure 2 .

 2 Figure 2. Turbulent round jet configuration. Top: Initial (left) and final (right) meshes of the center plane of the jet flow. Bottom: Center plane colored by the non-dimensional instantaneous velocity norm, |⃗ u|/(U 1 -U 2 ), computed on the initial (left) and the final (right) meshes.

NelemFigure 3 .

 3 Figure 3. Evolution of the global molecular dissipation, D Vν (top), the global transfer to TKE, P V TKE (middle), and the number of elements, N elem (bottom) during the different states of the automatic mesh convergence procedure in the turbulent round jet configuration. Note that the resolved, P V res , and the modeled P V mod , parts of the transfer to TKE are also shown.

Figure 4 .

 4 Figure 4. Comparison of turbulent round jet statistics: reference case (black), initial mesh (red) and final mesh (blue). Mean axial velocity profile (top) and rms axial velocity profile (bottom) at three sections: x/D = 1 (left), 5 (middle) and 10 (right).

Figure 5 .

 5 Figure 5. Turbulent pipe flow configuration. Mesh of a cross section and cross section colored by the nondimensional instantaneous velocity norm, |⃗ u|/U b , for the initial (left) and final (right) meshes.

NelemFigure 6 .

 6 Figure 6. Evolution of the global molecular dissipation, D V ν (top), the global transfer to TKE, P V TKE (middle), and the number of elements, N elem (bottom) during the different states of the automatic mesh convergence procedure in the turbulent pipe flow configuration. Note that the resolved, P V res , and the modeled P V mod , parts of the transfer to TKE are also shown.

13

 13

Figure 7 .

 7 Figure 7. Comparaison of turbulent pipe flow statistics: DNS (black), initial mesh (red) and final mesh (blue). Mean velocity profile (left) and rms velocity profiles (right) for the three directions: axial (solid line), wall-normal (dashed-dotted line) and azimuthal (dashed line).

Figure 8 .

 8 Figure 8. STL of the PRECCINSTA burner.

Figure 9 .

 9 Figure 9. Final mesh of the PRECCINSTA burner.

Figure 10 .

 10 Figure 10. Starting Eulerian mesh with the STL triangles and the created particles.

Figure 11 .

 11 Figure 11. Volume refinement kernel.

Figure 12 .Figure 13 .

 1213 Figure 12. Iso-contour of the level set function and mesh cut colored by the distance to the iso-contour.

Figure 14 .

 14 Figure 14. Projection of boundary nodes to the STL.

Figure 15 .

 15 Figure 15. Mesh cut colored by the edge length before coarsening.

Figure 16 .

 16 Figure 16. Mesh with the featured edges of the STL before coarsening.

Figure 17 .

 17 Figure 17. PRECCINSTA configuration. Top: Initial (left) and final (right) meshes of the center plane. Bottom: Center plane colored by the non-dimensional instantaneous velocity norm, |⃗ u|/Umax, computed on the initial (left) and the final (right) meshes.

NelemFigure 18 .Figure 19 .

 1819 Figure 18. Evolution of the global molecular dissipation, D V ν (top), the global transfer to TKE, P V TKE (middle), and the number of elements, N elem (bottom) during the different states of the automatic mesh convergence procedure in the PRECCINSTA burner configuration. Note that the resolved, P V res , and the modeled P V mod , parts of the transfer to TKE are also shown.

Figure A1 .

 A1 Figure A1. Turbulent round jet configuration starting from another initial mesh. Top: Initial (left) and final (right) meshes of the center plane of the jet flow. Bottom: Center plane colored by the non-dimensional instantaneous velocity norm, |⃗ u|/(U 1 -U 2 ), computed on the initial (left) and the final (right) meshes.

NelemFigure A2 .

 A2 Figure A2. Evolution of the global molecular dissipation, D V ν (top), the global transfer to TKE, P V TKE (middle), and the number of elements, N elem (bottom) during the different states of the automatic mesh convergence procedure in the turbulent round jet configuration starting from another initial mesh.

Figure B2 .

 B2 Figure B2. Comparison of turbulent round jet statistics: reference case (black), final mesh for the procedure based on QC 1 and QC 2 criteria (blue), final mesh for the procedure based on LIKE criterion (red, solid line) and the mesh obtained just before (red, dashed line). Mean axial velocity profile (top) and rms axial velocity profile (bottom) at three sections: x/D = 1 (left), 5 (middle) and 10 (right).

Acknowledgements

The authors gratefully acknowledge support from NETHUNS project under grant ANR-21-CHIN-0001-01. This work was granted access to the HPC resources of CINES/TGCC/IDRIS under the projects 2B06880 and 2A00611 made by GENCI. Part of this work has been initiated during the Extreme CFD Workshop & Hackathon (https://ecfd.coria-cfd.fr).

Appendix B. Turbulent round jet case using LIKE criterion

The proposed methodology is a mesh convergence procedure to guarantee an accurate mesh for LES. It is based on the premise that accurate LES should, at least, lead to mesh-independent mean fields, and then, a mesh-independent MKE balance. This methodology is user-independent and it does not require any a priori knowledge of the flow dynamics (as experimental data) conversely to previously proposed strategies [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF][START_REF] Daviller | A mesh adaptation strategy to predict pressure losses in les of swirled flows[END_REF][START_REF] Agostinelli | Static mesh adaptation for reliable large eddy simulation of turbulent reacting flows[END_REF]. In this paper, this methodology is coupled with the two criteria proposed by Benard et al. [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF] to determine how the new mesh has to be generated. However, the proposed methodology is expected to work with other mesh definition criteria. The same mesh convergence procedure is then applied using the LIKE criterion [START_REF] Daviller | A mesh adaptation strategy to predict pressure losses in les of swirled flows[END_REF] for the turbulent round jet case (section 3.4.1) to confirm that.

The LIKE criterion proposes to base the mesh-refinement on

which is the molecular and turbulent dissipation of the average of the overall resolved kinetic energy in LES, K = 1/2⟨ū 2 i ⟩. From this quantity, the adapted local mesh size, ∆ a (⃗ x), is then defined from the initial local mesh size, ∆ i (⃗ x), as

where Φ min and Φ max correspond to the minimum and the maximum of Φ in the whole computational domain, and where α and ϵ are two user-parameters of the metric definition corresponding to a smoothing parameter and to the maximum refinement ratio allowed. This criterion applies this maximum refinement ratio where Φ is maximum, and keeps the initial metric where Φ is minimum. The recommended range of values to use for the criterion parameters are 0.3 ≤ ϵ ≤ 0.7 and 30 ≤ α ≤ 100. In the present case, the LIKE criterion is applied with ϵ = 0.5 and α = 65. The proposed methodology adapted to the use of LIKE criterion is finally simpler than the methodology proposed for the two criteria (section 3.3). Starting from the same coarse mesh composed 307, 806 elements, the mesh is adapted based on the metric definition of the LIKE criterion, Eq.(B2), after each statistical time window, T w = 170R/(U 2 -U 1 ). This is repeated until the deviations of the global molecular dissipation, ∆D V ν , and the global turbulent production, ∆P V TKE , are both smaller than 0.05. Figure B1 shows the evolution of the global molecular dissipation and the global turbulent production as a function of the number of elements during the mesh convergence procedure. The procedure based on LIKE criterion is compared to the procedure based on QC 1 and QC 2 criteria presented in section 3.4.1. At the end of the mesh convergence procedures, the global molecular dissipation and the global turbulent production are very similar, showing that the same global MKE balance is recovered. However, the procedure based on LIKE criterion leads to a final mesh composed by 102, 354, 009 elements, which is important in comparison with the final mesh composed by 8, 125, 037 elements obtained with the procedure based on QC 1 and QC 2