Martin David
email: martin.david@inria.fr

Mahitosh Mehta

Remi Manceau
email: remi.manceau@inria.fr

Rémi Manceau

Rémi Manceau

R Manceau

A self-adaptive strategy for hybrid RANS/LES based on physical criteria

Keywords:

A self-adaptive strategy for hybrid RANS/LES is proposed.

• The mesh and the model are adapted based on physical criteria, with no intervention from the user.

• Resolved structures are rapidly generated in the LES region utilizing an active forcing.

• The approach is assessed on a flow over a backward facing-step.

• The method converges and the results are improved when compared to RANS k-SST.

Introduction

The Reynolds-Averaged Navier-Stokes (RANS) method makes it possible to investigate industrial cases due to its low computational cost, but often lacks reliability and only produces statistical results. On the other hand, LES, by only modeling the small turbulence scales, generally provides more accurate results and gives access to unsteady information. However, its computational cost remains too high in wall-bounded high Reynolds number industrial configurations. Even if the computing power is expected to continue increasing, Large Eddy Simulation (LES) should remain too expensive in many domains in a foreseeable future [START_REF] Pope | Ten questions concerning the large-eddy simulation of turbulent flows[END_REF].

In most cases, it is desirable to resolve a larger range of the turbulence spectrum to get more information on the flow and handle strongly swirling flows, large-scale separation, acoustics sources, etc. Hybrid RANS/LES models are thus very promising since they can drastically reduce the computational cost of LES, give access to time-dependent information, and provide more reliable results than RANS. In these approaches, modeled motions, associated with the RANS and the LES models, and resolved motions, associated with the unsteady modes captured in the LES region, coexist. Hybrid RANS/LES methods may be divided into two categories: zonal, or discontinuous, and non-zonal, or continuous, approaches. In the first case, two different models are used in two subdomains. The main challenge with these approaches is to deal with the boundary conditions at the interface between the two modes. In non-zonal, or continuous approaches, only one set of equations is used and the transition from RANS to LES relies on model parameters updated during the computation, which makes them more flexible and is the reason for growing interest. Some of the most usual non-zonal hybrid RANS/LES approaches are the family of Detached Eddy Simulation (DES) models [START_REF] Spalart | Strategies for turbulence modelling and simulations[END_REF], the Partially-Integrated Transport Model (PITM) [START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF], the Partially-Averaged Navier-Stokes (PANS) method [START_REF] Girimaji | Partially-Averaged Navier-Stokes Model for Turbulence: A Reynolds-Averaged Navier-Stokes to Direct Numerical Simulation Bridging Method[END_REF], and the Scale-Adaptive Simulation (SAS) [START_REF] Menter | The Scale-Adaptive Simulation Method for Unsteady Turbulent Flow Predictions. Part 1: Theory and Model Description[END_REF]. The reader is referred to Chaouat [START_REF] Chaouat | The State of the Art of Hybrid RANS/LES Modeling for the Simulation of Turbulent Flows[END_REF] or [START_REF] Heinz | Hybrid RANS-LES Methods with Continuous Mode Variation[END_REF] for recent reviews of hybrid RANS/LES methods.

In most hybrid RANS/LES approaches, the model is explicitly mesh size-dependent, making the position of the RANS and LES zones highly dependent on the mesh design chosen by the user. Based on his intuition and/or experience of the flow under consideration, and his knowledge of the limitations of RANS models, the user will choose to refine the mesh where he feels it is necessary to switch to LES mode. This is illustrated by Fig. 1 for the example of a backstep flow: two users may generate very different meshes, yielding different LES regions and, consequently, different results. Such user influence is highly dangerous in an industrial context, and in order to ensure the reliability of CFD for the design, dimensioning, and certification of industrial systems [START_REF] Park | Unstructured Grid Adaptation: Status, Potential Impacts, and Recommended Investments Towards CFD[END_REF][START_REF] Spalart | On the role and challenges of CFD in the aerospace industry[END_REF], mitigating the user influence is of major importance.

Some approaches, often referred to as self-adaptive, break the explicit link between model and local grid step. The switch from RANS to LES is then driven by a physical criterion evaluated during the computation, such as the ratio between the integral and von Kármán length scales [START_REF] Menter | The Scale-Adaptive Simulation Method for Unsteady Turbulent Flow Predictions. Part 1: Theory and Model Description[END_REF][START_REF] Krumbein | A Scale-Resolving Elliptic-Relaxation-Based Eddy-Viscosity Model: Development and Validation[END_REF] or the ratio of modeled to total turbulent energy [START_REF] Girimaji | Partially-Averaged Navier-Stokes Model for Turbulence: A Reynolds-Averaged Navier-Stokes to Direct Numerical Simulation Bridging Method[END_REF][START_REF] Perot | A self-adapting turbulence model for flow simulation at any mesh resolution[END_REF][START_REF] Perot | A stress transport equation model for simulating turbulence at any mesh resolution[END_REF][START_REF] Baglietto | STRUCT: A Second-Generation URANS Approach for Effective Design of Advanced Systems[END_REF]. This type of models, called second generation URANS models by Fröhlich and von Terzi [START_REF] Fröhlich | Hybrid LES/RANS methods for the simulation of turbulent flows[END_REF], introduces a partial decoupling between the model and the mesh, enabling the model, for example, to decide to remain in RANS mode even if the mesh is fine. The physical criteria used enable the model to detect whether turbulent structures appear locally in the resolved field, and to reduce the turbulent viscosity or the Reynolds stresses accordingly. However, it is clear that the decoupling is not complete, since the size of the resolved structures is still constrained by the grid, and the model cannot switch to LES mode if the mesh is too coarse, such that the user's influence remains very significant.

A way of avoiding the latter problem is that the mesh size is not set by the user, but by an adaptive meshing process. A few studies have explored hybrid RANS/LES methods focused on investigating mesh self-adaptivity within a predetermined and fixed LES region. Limare et al. [START_REF] Limare | Adaptive Mesh Refinement with an Automatic Hybrid RANS/LES Strategy and Overset Grids[END_REF], using DES, proposed remeshing based on physical criteria, to adapt the mesh in the LES zone using a simulation strategy combining an octree-AMR and overset grids. The grid size is partly controlled by one of the three following criteria, depending on the case study: the laminar-to-turbulent viscosity ratio, the von Kármán length scale-to-grid step ratio, or the Taylor length scale-to-grid step ratio. Assessing local LES-resolution sensors for hybrid RANS/LES simulations, Reuß et al. [START_REF] Reuß | Assessment of Local LES-Resolution Sensors for Hybrid RANS/LES Simulations[END_REF] performed Improved Delayed Detached-Eddy Simulation (IDDES) of a flow over a backward-facing step. They studied different formulations of a mesh sensor that detect under-resolved regions in LES region and can be used as input for automatic and local mesh refinement. Specifically, they proposed to refine the mesh in the LES region if the resolved turbulent kinetic energy is less than an arbitrary target of 80% of the total turbulent kinetic energy outside of the RANS region.

One step further is the work of Woodruff [START_REF] Woodruff | Adaptive Embedded LES of the NASA Hump[END_REF], in which the extent of the LES region is also self-adaptive. The location and size of the LES region are imposed by the user at the beginning of the workflow, then the region is allowed to grow until the solution becomes insensitive to this growth, and the grid size is adapted based on the local value of the strain rate. In this approach, the position of the LES zone is still determined by the user, based on his a priori knowledge of the flow, but the size of this zone and the mesh refinement are determined adaptively.

The aim of the present study is to develop a fully selfadaptive approach in which the user does not influence the definition of the RANS and LES regions and the resolution in the LES region. Starting from a RANS computation with a mesh obtained after a grid convergence study, the proposed method autonomously determines both the LES and RANS/LES transition regions, along with the size of the cells in the LES region based on physical criteria, as illustrated in Fig. 1. As RANS results are grid converged, they are considered here to be independent of the user. However, it is possible to rely on a RANS calculation with automatic mesh refinement, so as to completely avoid any user influence on the initial state of our algorithm. Successive hybrid RANS/LES computations are conducted and the characteristics of the LES and transition regions are updated after each simulation until convergence is achieved. As a demonstrator, the strategy is applied to the backward-facing step case [START_REF] Lamballais | Direct numerical simulation of a turbulent flow in a rotating channel with a sudden expansion[END_REF] using the HTLES approach [START_REF] Duffal | Development and Validation of a new formulation of Hybrid Temporal Large Eddy Simulation[END_REF], in its active version [START_REF] Mehta | An active hybrid Reynolds-averaged Navier-Stokes/large eddy simulation approach for gray area mitigation[END_REF].

The paper is organized as follows: our self-adaptive strategy is described in section 2, focusing on the possible physical criteria to identify the region where LES is necessary (section 2.1), the way to determine the size a RANS/LES transition region (section 2.2), and criteria used to fix the targeted grid step in the LES region (section 2.3). The hybrid RANS/LES model used for testing this strategy, the active HTLES of Mehta et al. [START_REF] Mehta | An active hybrid Reynolds-averaged Navier-Stokes/large eddy simulation approach for gray area mitigation[END_REF] is described in section 3. Finally, in section 4, the self-adaptive strategy is evaluated in the case of a backward-facing step, by investigating the evolution of the RANS and LES zones and the associated mesh during successive iterations of the algorithm, as well as comparing of the results with Direct Numerical Simulation (DNS) data.

Self-adaptive strategy for hybrid RANS/LES

The proposed self-adaptive strategy is summarized in Fig. 2. This method aims at autonomously determining the LES and transition zones as well as the cell size in the LES zone. It consists in successive simulations to iteratively update the zones and the cell size until the convergence of the algorithm is reached. More precisely, the algorithm is initialized by a grid-converged RANS calculation. The RANS results are used to evaluate, based on a physical criterion, in which region to switch to LES mode. Two nested boxes are defined, delimiting the LES region and the RANS/LES transition region. The mesh in these two boxes is then refined to satisfy a resolution criterion. A new calculation is performed in hybrid RANS/LES using this automatically refined mesh. The physical criterion for determining the LES zone is updated on the basis of the average solution obtained in this hybrid RANS/LES simulation. If the LES zone is larger than previously, then a new iteration of the algorithm is performed. Otherwise, the algorithm stops and the solution obtained at this iteration is the result of the auto-adaptive process.

The three steps highlighted in yellow correspond to pivotal aspects of the present study and are detailed in subsections 2.1, 2.2, and 2.3.

Where to switch to LES?

The physical criteria used to define the refinement region in the framework of self-adaptive hybrid RANS/LES must be (i) based on Reynolds-averaged quantities, (ii) independent of the grid step, (iii) non-dimensional, (iv) objective (independent of the reference frame).

Based on these constraints, several physical criteria have been investigated from RANS results obtained in the case of a backward-facing step [START_REF] Lamballais | Direct numerical simulation of a turbulent flow in a rotating channel with a sudden expansion[END_REF], with a bulk Reynolds number = ∕ = 5000, where is the half-width of the channel and the bulk velocity downstream of the expansion. The expansion ratio ℎ∕ is 2∕3. Only the lower half of the domain is considered with a symmetry applied at ∕ = 1. The extent of the computational domain is 4 upstream and 16 downstream. In the spanwise direction, the domain size is and the flow is periodic. A periodic RANS channel flow computation is used as the inlet boundary condition. The RANS model is the --SST model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF].

Note that in the present work, all the computations are carried out using code_saturne, a general CFD solver developed by EDF [START_REF] Archambeau | Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows -Industrial Applications[END_REF], based on a finite-volume method and a fully collocated arrangement for all the variables. A standard predictor-corrector scheme (SIMPLEC) is used to solve the velocity-pressure system. The Crank-Nicholson time scheme and a second-order (upwind in the RANS zone, centered in the LES zone) convection scheme are used. For hybrid RANS/LES, the statistical quantities are evaluated using both temporal averaging and spatial averaging in the periodic direction ( ).

Numerous physical criteria corresponding to the constraints listed above have been considered and evaluated on the basis of RANS results. Those that seem most relevant, at least for the case of the backward-facing step, are presented in Fig. 3. They involve the turbulent kinetic energy, , the production of turbulent kinetic energy, , the dissipation, , the strain rate ̃ = √ 2 , the Taylor length scale, = √ 10 ∕ , the Von Kármán length scale, = ̃ ∕| |, where = 0.41 and | | is the Froebenius norm of the Hessian matrix of the velocity. Overall, the criteria studied show a sharp increase downstream of the step, then the region colored in red extends downstream and moves away from the bottom wall. This is explained by the flow separation and the development of a reattached boundary layer. In Figs. 3a and3b, the criteria exhibit high values in the detached shear layer, the near wall region, and the beginning of the recirculation zone. As the near-wall region is intended to stay in RANS mode, they are not adopted for this self-adaptive strategy. The third criterion presented (Fig. 3c) shows three thin red zones: the first follows the wall and is mainly due to the term ∕ 3∕2 . The second and the third zones are induced by the first-order derivative-tosecond-order derivative ratio involved in the Von Kármán length scale. This criterion, which is used in SAS to identify regions where turbulent structures appear in the resolved field, is therefore unsuitable for identifying detached shear layers from a Reynolds-averaged solution. The last two criteria quite satisfactorily identify the detached shear layer. For these reasons, any of the two first criteria could be used in the case of the backward-facing step. Combining different criteria to make the adaptive approach ready for a wide range of applications is beyond the scope of the present paper.

Therefore, in order to demonstrate the applicability of the adaptive strategy, the criterion chosen herein is simply the production-to-dissipation ratio. This criterion makes it possible to separate the zones close to equilibrium, where RANS models are generally sufficiently accurate, and the out-of-equilibrium zones where it is preferable to switch to LES mode to improve the quality of the results. Thus, regions where to switch to LES are identified via the threshold

(a) 2 ̃ 2 (b) 3∕2 (c) 3∕2 (d) ̃ (e)
> , ( 1
)
where and are the total production and total dissipation, respectively, and the threshold is chosen as = 1.4. Figure 4 shows the region identified by this criterion. The total production is obtained from

= -′ ′ + m . ( 2
)
′ ′ is the Reynolds-averaged (in practice, averaged in time and in span), which is zero in a RANS computation. m is the Reynolds-averaged sub-filter stress tensor, equal to the Reynolds-stress tensor in RANS mode, and is the mean strain rate tensor. The total dissipation is computed as

= m ∕ , ( 3
)
where m is the Reynolds-averaged sub-filter dissipation. is the ratio of the modeled dissipation m to the total dissipation , which is equal to one in RANS mode and also generally considered very close to one in LES mode, since the small scales are not resolved [START_REF] Duffal | Development and Validation of a new formulation of Hybrid Temporal Large Eddy Simulation[END_REF]. However, in low Reynolds number regions, such as the recirculation region downstream the backstep, a small part of the dissipation can be resolved. This must be accounted for in order not to bias the evaluation of the criterion [START_REF] Archambeau | Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows -Industrial Applications[END_REF]. The determination of in the framework of the HTLES model used herein will be presented in section 3.

How to ensure the RANS/LES transition?

Bounding boxes are defined around the region satisfying the physical criterion given by Eq. ( 1) in order to adapt the mesh, as shown in Fig. 4. For practical reasons, the RANS, RANS/LES transition, and LES zones are denoted A, B, and C, respectively. They are obtained using the following rules:

(i) The LES zone (zone C) is the smallest rectangle containing the ∕ > region. (ii) The transition zone (zone B) is used to smoothly link the meshes used in LES and RANS zones. It is characterized by a rectangular shape bounding the LES zone with a thickness = 0.5 ⟨ ⟩ ∕ > where ⟨ ⟩ ∕ > is the mean of the integral length scale in the zone satisfying Eq. ( 1). If at the end of a simulation the region corresponding to ∕ > extends until the bounds of the transition zone, the thickness of this zone is multiplied by 10 in this direction for the next iteration of the algorithm, to reduce the total number of iterations to reach the convergence of the adaptive method. Upstream, the beginning of the transition zone is fixed with a length of

= , ( 4
)
where = 5 is a constant, which is determined by a preliminary study presented in Sec. 4.1. This is the minimal length required for the volume forcing presented in Sec. 3.2 for generating turbulent structures. Note that, since a shielding function is applied in the hybrid RANS/LES model until + = 20 to enforce the RANS mode in the near wall region, this shielded region is excluded from zone B. (iii) The RANS zone (zone A) is the rest of the domain.

The RANS mode of the hybrid RANS/LES model is imposed. In this zone, the mesh will not be refined.

What is the targeted cell size in each region?

In order to provide a user-independent method, the mesh size in the LES region must now be set by an objective criterion. The criterion chosen here is the ratio of modeled energy to total turbulent energy . Assuming a Kolmogorov spectrum, we have

= m = 1 ∫ ∞ ( )d = 3 2 3∕2 -2∕3 , (5)
with = ∕Δ the cutoff wavenumber and = + . The cell size to reach the targeted energy ratio is thus computed as

Δ = 2 3 3∕2 3∕2 . ( 6
)
Since it is linked to the integral length scale 3∕2 ∕ , Δ is variable in space and is given as a local target to the meshing tool.

The rules followed to build the new mesh are:

(i) In the RANS zone (zone A): the mesh is kept unchanged. (ii) In the LES zone (zone C): the new unstructured mesh is built, composed of cells whose local target size is given by Eq (6). (iii) In the transition zone (zone B): the targeted cell size in the transition zone is equal to the maximum cell size obtained in the LES zone and the smooth transition between the structured RANS and unstructured LES meshes is ensured by the meshing tool.

In this study, the GMSH software [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] is used to generate the grids. The unstructured zones are meshed with the Packing of Parallelograms algorithm.

Turbulence modeling

Although the self-adaptive methodology can be associated to any continuous hybrid RANS/LES approach, HTLES---SST [START_REF] Duffal | Development and Validation of a new formulation of Hybrid Temporal Large Eddy Simulation[END_REF] is used in the present study, in its active version recently proposed by Mehta et al. [START_REF] Mehta | An active hybrid Reynolds-averaged Navier-Stokes/large eddy simulation approach for gray area mitigation[END_REF]. For more details, the reader is referred to these two papers.

Hybrid Temporal LES (HTLES)

Derived from the PITM model [START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF][START_REF] Fadai-Ghotbi | Temporal filtering: A consistent formalism for seamless hybrid RANS-LES modeling in inhomogeneous turbulence[END_REF], the HTLES is a continuous hybrid approach in which the hybridization term in the energy equation is based on a time scale driven by the energy ratio . The equations resolved are the following:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ sfs =
= ( ) m + m * m , ( 12
)
where is the coefficient imposing the internal consistency constraint: in the RANS zone, resolved energy due to the penetration of turbulent structures from the LES region is not counted in the total turbulent energy in order to strictly tend to the RANS model in this zone,

= 0 if = 1, if < 1. ( 13
)
( ) is the hybridization function defined by

( ) = + ( - ) . ( 14
)
This expression goes to unity when = 1, such that the RANS model is recovered. The energy ratio

= (1 -) + × min 1, . ( 15
)
involves the shielding function introduced by Duffal et al. [START_REF] Duffal | Development and Validation of a new formulation of Hybrid Temporal Large Eddy Simulation[END_REF] to avoid grid-induced separation and the log-layer mismatch:

= 1 -tanh max 1 , 2 , (16)
with = 1 3 ∕ 1∕4 ∕ , and = 2 Δ ∕ , where is the distance to the wall and Δ the local cell elongation (longest edge of the cell). The energy ratio , active far from the wall where the shielding function goes to one, is obtained analytically from an equilibrium Eulerian temporal spectrum [START_REF] Tennekes | Eulerian and Lagrangian time microscales in isotropic turbulence[END_REF],

= 1 0 √ 2∕3 -2∕3 (17)
where

= m ( ) * m , (18)
= min , Δ is the cutoff frequency, Δ = Ω 1∕3 with Ω the cell volume, and = + √ is the sweeping velocity. Note that in the present study, the time-step is smaller than Δ∕ everywhere in the domain, such that = ∕Δ and is given by Eq. 5. The constants used are summarized in Tab. 1 and more details can be found in [START_REF] Duffal | Development and Validation of a new formulation of Hybrid Temporal Large Eddy Simulation[END_REF].

As mentioned in Sec. 2.1, the criterion ∕ > is biased if the resolved part of the dissipation is not taken into account, particularly in the low Reynolds number recirculation region. In order to avoid the cumbersome computation of the resolved part of dissipation during the simulation, Eq. ( 3) is used, where the modeled-to-total dissipation ratio is evaluated, assuming a Kolmogorov spectrum, from the definition of the energy and dissipation ratios, and , respectively, and the integration of the energy and dissipation spectra between the integral and the Kolmogorov length scales. Denoting by a hat variables made non-dimensional based on , and , we have

= ̂ m ̂ = 1 ̂ ∫ ̂ ̂ 2 ̂ 1∕3 d ̂ , ( 19
)
where = ∕ is the turbulent Reynolds number, yielding

= 1 - 2 Δ4∕3 . ( 20
)
Similarly, from

= ̂ m ̂ = 1 ̂ ∫ ̂ ̂ ̂ -5∕3 d ̂ , (21)
one can obtain

= Δ2∕3 2∕3 1∕2 . ( 22
)
Combining Eqs. ( 20) and ( 22) gives

= 1 -1 + 1∕2 -2 . ( 23
)
Since dissipative scales are smaller than energetic scales, ≥ , such that the final expression is

= max ⎡ ⎢ ⎢ ⎢ ⎣ 1 - ⎛ ⎜ ⎜ ⎝ 1 + 1∕2 ⎞ ⎟ ⎟ ⎠ -2 , ⎤ ⎥ ⎥ ⎥ ⎦ . ( 24
)
With this expression, = 1 in the RANS mode since = 1 and goes to zero at the DNS limit where = 0. = 0.42 is a coefficient calibrated by comparing dissipation fields obtained in hybrid RANS/LES and RANS.

Active approach

In continuous hybrid RANS/LES, during the transition from RANS mode to LES mode, the modeled turbulent energy is reduced by the model. However, the growth of the resolved energy is slow, leading to a strong local underestimation of the total turbulent energy and the turbulent stresses. To mitigate this gray area issue, Mehta et al. [START_REF] Mehta | An active hybrid Reynolds-averaged Navier-Stokes/large eddy simulation approach for gray area mitigation[END_REF] have developed the active approach, which injects the proper amount of energy into the resolved motion via a body force in the momentum equation written as

= ̃ + , ( 25
)
where ̃ is the resolved velocity. and have to satisfy two constraints: (i) energy injection in the resolved scales is imposed by the rate of energy removal in the modeled scales; (ii) the forcing should not affect the mean flow. It can be shown that these constraint are satisfied if and are determined by solving at each point of the domain the system

⎧ ⎪ ⎨ ⎪ ⎩ ′ ′ + ′ ′ = -m 1 d d , = -̃ , ( 26
)
where ′ is the fluctuating part of the resolved velocity and ⋅ stands for the Reynolds average. The force is thus a function of the rate of change of the energy ratio in the flow direction d ∕d and goes to zero where the mesh is uniform. This active version of the HTLES is used for all the simulations in the present paper.

Results and discussion

Influence of the length of transition

The intensity of the fluctuating body force added to the momentum equation is determined by Eq. ( 26) and relies on the material derivative of the energy ratio . Thus, increasing the length of the RANS-to-LES transition upstream the step corner (Fig. 4) has different and contradictory consequences. On the one hand, since the refinement of the mesh in the flow direction is less rapid, the energy ratio = ∕ decreases more slowly and the intensity of the forcing is reduced. On the other hand, the resolved structures have more space for developing, since the force is active over a longer region. The influence of is in theory moderate, since these two effects compensate for each other: (1∕ ) × d ∕d ≃ (1∕ ) × ∕ is of the order of ∕ , where is the streamwise mean velocity, and since the distance travelled is , the work done by the force is globally independent of . However, Mehta et al. [START_REF] Mehta | An active hybrid Reynolds-averaged Navier-Stokes/large eddy simulation approach for gray area mitigation[END_REF] showed, in the case of a channel flow, that the realistic turbulent strucures cannot develop if the transition region is too short. Therefore, a preliminary study is performed here to determine the value of to be used in Eq. ( 4). In this preliminary study, a structured mesh of 5.7M cells is used, which is fine enough to be in LES mode everywhere. However, the RANS mode is artificially imposed in the inlet channel, as well as a transition in the region ∈ [ trans ; 0] upstream of the step. To achieve this, the energy ratio related to the grid step given by Eqs. ( 15) and ( 17) is modified by imposing mod = 1 in the inlet channel up to the position = trans and then a gradual decrease of the form

mod = + (1 -) with ( ) = trans , ( 27
)
such that the usual value of is recovered at the step corner ( = 0). Two different cases are presented here, for trans ∕ = -1 and trans∕ = -2, as shown in Fig. 5. For both cases, the resolved structures are rapidly developed, as soon as the energy ratio is lower than unity and the structures are advected downstream. The streamwise velocity profiles in Fig. 6 show that, for the two cases, the profiles are nearly superimposed, showing that trans ∕ = -1 is sufficient. Reducing further the transition zone does not allow to generate realistic structures (not shown here). Since a slightly better agreement with DNS is found for the case with the shorter transition length, and with the purpose of reducing the global cost of the computations, this order of magnitude of ∕ ≈ 1 is retained for building the zone B (see Fig. 4). Based on the RANS results, this approximately correspond to calibrating the constant in Eq. (4) as = 5.

Evolution of the zones and the mesh

The region where switching to LES mode is desirable obtained from Eq. ( 1), using = 1.4. The cell size to reach the targeted energy ratio is given by Eq. ( 6). The targeted value = 0.1 is used, except for the first iteration, called HTLES1 below, for which the less strict value = 0.15 is applied to save computational time, since it is only the first mesh determined by the probably inaccurate RANS results . The mesh and zone obtained for the successive computations are presented in Figs. 7 and8.

A RANS computation, performed on a fully structured grid composed of 95 800 cells obtained after a grid convergence study, is used to initialize the algorithm. Note that the mesh consists of only one cell in the direction due to spanwise periodicity. Following the rules presented in Sec. 2.2, the region corresponding to ∕ > is extracted and used to determine the boundaries of the LES and transition zones of the first HTLES computation (purple boxes). The cell size in the LES region is computed following the rules described in Sec. 2.3. The HTLES1 mesh is generated by GMSH according to these zones and cell size requirements. In the HTLES computations, the grids in the transition and LES zones are also refined in the -direction to obtain cells which are as isotropic as possible. The HTLES1 computation is then performed and the region corresponding to ∕ > is updated based on the time-averaged results. Here, it extends until the end of the transition zone thus, following the rules established in Sec. 2.2, the thickness of the transition zone utilized for HTLES2 (largest red box at the bottom of Fig. 7) is increased in this direction to mitigate the total number of algorithm iterations. The same procedure is repeated for HTLES2 and HTLES3 (Fig. 8). The results of HTLES3 show that the area of the LES zone is reduced when compared to the one of HTLES2. This stops the algorithm, since the LES region is sufficiently large, and trying to optimize its size would lead to extra iterations of the algorithme that would actually increase the total computational cost. Here, we still performed a fourth HTLES run, using the same zones but cell sizes computed from HTLES3 results, to confirm the convergence of the results.

The cell number, the averaged value of Δ, the cell size to reach the targeted energy ratio in the region where the criterion (1) is satisfied and the area of zone C are shown in Fig. 9. Except for HTLES2 whose transition zone size has been increased, as explained above, to reduce the total number of algorithm iterations, the grid size exhibits a monotonic increase from RANS to HTLES4. The mesh obtained after the computation of HTLES4, which would have been used for HTLES5 if the algorithm were not stopped, is reduced when compared to the one obtained after HTLES3, highlighting the convergence of the method. After a monotonic increase, the value ⟨Δ⟩ ∕ > obtained with HTLES3 is reduced when compared to HTLES2, and remains virtually constant between HTLES3 and HTLES4. A similar behavior is observed in the graph showing the evolution of the area of Zone C, with the stabilization of the LES zones. From these graphs, it is clear that the proposed method reaches convergence.

Comparison with DNS

In this section the results obtained after the different iterations of the algorithm are compared with DNS. The results of HTLES4 computation are not shown since they are virtually superimposed with those of HTLES3.

In Fig. 10, the target energy ratio (i.e., used in the model, given by Eq. ( 15)), the effective energy ratio (i.e., computed during the simulation), the streamwise velocity, and the turbulent kinetic energy profiles are plotted. Since the DNS and RANS energy ratios are zero and one, respectively, throughout the entire domain, they are not shown. The HTLES energy ratio profiles depict the interface between RANS and LES regions. The area of the LES region significantly increases from HTLES1 to HTLES2, then stabilizes, reflecting the convergence of the self-adaptive method. The effective energy ratio profiles (Fig. 10b) globally follow those of the energy ratio even if some discrepancies appear, reflecting the fact that the partition of energy approximately reaches the expected level. The results indicate that downstream of the LES region, some resolved scales remain. Mehta et al. [START_REF] Mehta | An active hybrid Reynolds-averaged Navier-Stokes/large eddy simulation approach for gray area mitigation[END_REF] made the same observation for the case of a periodic hill, and concluded that it was not a problem, since their results are actually better than pure RANS results. Note that close to the wall, the effective energy ratio is not relevant since both and go to zero. The spikes in the profiles are due to the sharp variations of the total turbulent kinetic energy (Fig. 10d).

Regarding the streamwise velocity (Fig. 10c), upstream the transition region, the RANS mode is imposed and periodic RANS solution imposed at the inlet is preserved. The mesh refinement from = 0.15 (HTLES1) to = 0.10 (others) in the LES region and the definition of the zones has little influence on the streamwise velocity and turbulent kinetic energy (Fig. 10d) profiles, which highlights the robustness of the method. The expansion of the LES zone in the recovery region (HTLES2 and HTLES3) improves the turbulent kinetic energy profiles at ∕ = 3 and, for HTLES2, which has the largest LES zone, until ∕ = 4. The spikes of turbulent kinetic energy, principally observed at ∕ = 0, are due to the rapid RANS/LES transition in the vertical direction and to the forcing term that depends on the variations of the energy ratio. The overshoot of turbulent kinetic energy observed near the step in all HTLES rapidly fades out. A similar behaviour was observed by Mehta et al. [START_REF] Mehta | An active hybrid Reynolds-averaged Navier-Stokes/large eddy simulation approach for gray area mitigation[END_REF] for periodic hills and channel flows. Comparing RANS with HTLES3, it appears that RANS significantly underestimates the amount of turbulent kinetic energy from ∕ = 0 to ∕ = 2, while it is well estimated by the HTLES3. A close agreement with DNS is found between ∕ = 1 and ∕ = 3, which makes possible an accurate prediction of the streamwise velocity profiles between ∕ = 2 and ∕ = 4. Downstream of the location ∕ = 4, the HTLES streamwise velocity profiles are slightly improved when compared to the RANS profiles.

The total Reynolds stresses, ′ = ( -) + m , and the corresponding anisotropy components, defined as

= ′ 2 - 1 3 , ( 28
)
are shown in Fig. 11 for HTLES3, at the streamwise locations ∕ = 0 (top) and ∕ = 1.5 (bottom). Figs. 11a and 11b show that the fluctuations generated by the forcing do not exhibit the correct anisotropy upstream of the step. This is explained because the volumic forcing acts as an amplifier of the preexisting fluctuations based on Eq. ( 26), in which m is the anisotropy tensor given by the linear eddy viscosity models. Correct anisotropy could only be obtained if the RANS model itself was capable of representing it, as with a Reynolds-stress model or an explicit algebraic stress model. However, at ∕ = 1.5, the Reynolds stresses is quite well predicted, specifically for the ′ 11 and ′ 12 components. In Fig. 11d, despite the poor prediction of the anisotropy in ∕ = 0, the results observed in the LES region are remarkably close to those of the DNS.

The skin friction distribution on the lower wall, shown in Fig. 12, highlights the benefits associated with the selfadaptive method. Indeed, the HTLES3 profile leads to a relatively correct prediction of the negative peak of and gives a very accurate estimate of the reattachment location. Note that the results of HTLES1, which mesh is significantly coarser than that of HTLES3 in the transition and LES regions (see Figs. 7 and8), are quite accurate, underlying the robustness of the HTLES to mesh coarsening. HTLES gives remarkably better results than RANS in the recirculation zone and downstream of it. Moreover, the peak of observed in DNS at ∕ ≈ 0.1, associated with the secondary counterclockwise recirculation zone, is well reproduced by HTLES while it is not predicted by RANS. Far in the recovery region, HTLES tends toward RANS profiles since the RANS mode is reactivated.

Mean streamlines computed using the line integral convolution (LIC) are displayed in Fig. 13. The two classically observed recirculation zones [START_REF] Chen | A review of Backward-Facing Step (BFS) flow mechanisms, heat transfer and control[END_REF] are visible. The main recirculation zone (clockwise) extends up to ∕ = 2.1 and the secondary recirculation zone (counterclockwise) to ∕ = 0.4, to be compared to the DNS values, ∕ = 2.13 and ∕ = 0.63, respectively.

The PSD are computed from instantaneous velocity saved every time step, d = 0.0012 second, during 165 seconds of simulations, resulting in 137 500 instants. For each streamwise location, monitoring points are saved at nine evenly spaced locations in the periodic direction. The data are concatenated leading to a sample size of = 1 237 500 velocities for each PSD. The results are smoothed to improve the readability of the figures. The power spectral densities are defined by

∫ +∞ 0 ( )d = .
(29)

In Fig. 14, the PSD reflect the influence of the forcing on the flow, aiming to rapidly develop turbulence structures to compensate for the decrease of the modeled turbulent energy. The peaks observed for > 3 Hz highlight the effect of the forcing on the relatively high-frequency structures while the range corresponding to < 3 Hz is not affected. Note that a Strouhal number of 0.1 is found which is in agreement with the study of Celenligil and Mellor [START_REF] Celenligil | Numerical Solution of Two-Dimensional Turbulent Separated Flows Using a Reynolds Stress Closure Model[END_REF]. The inertial range is scaled with the classical -5∕3 slope. Because of its lower energy ratio, the probe located at ∕ = 0 shows a higher amplitude of the PSD and the sharp decay of the spectrum occurs for higher frequencies than that of ∕ = -0.5.

The PSD computed downstream of the step are plotted in Fig. 15. In the region of the downstream sensors, the forcing is virtually null since the variations of the energy ratio are low. Hence, contrary the probes located in the forcing region, the PSD decrease monotonously from the inertial to the dissipation ranges. The closer to the wall, the lower the amplitude of the spectra, which is in agreement with the total turbulent kinetic energy profiles. In the inertial range, the -5∕3 slope is well reproduced for the three + locations. In the high-frequency region, the dissipation scale plays a major role, inducing a fast decay of the turbulent spectrum. The sharp decreases of the slope for very high frequency should not be analyzed since it is beyond the cutoff frequency (colored in light gray). Overall, the PSD express the physical consistency of the turbulent structures generated in the LES region. At ∕ = 4, the RANS mode is activated and the resolved energy is reduced. The turbulent structures are thus partially filtered out which explains the lower PSD amplitude and the sharper slope in the decay region when compared to the other locations.

The Q-isosurfaces presented in Fig. 16 for HTLES3 show that the active HTLES rapidly develops resolved scales, thanks to the volume forcing. Indeed, turbulent structures are remarkably observed right from the beginning of the LES region. However, downstream of the LES region, the structures are advected into the RANS zone and are only gradually dissipated by the RANS eddy viscosity. As can be seen in Fig. 12, the solution transitions in this region from a LES to a RANS behavior.

Conclusions and future work

In this paper, a self-adaptive strategy is proposed to mitigate the influence of the user on the results of hybrid RANS/LES computations. To avoid the definition of the RANS and LES zones by the user choices made at the time of the construction of the mesh, the zones where to switch to LES are determined by physical criteria, and the mesh is also refined in this zone to target a resolution obtained based on a physical criterion. The method is assessed on a backward-facing step with the HTLES---SST model.

Several possible physical criteria within the context of self-adaptive hybrid RANS/LES are investigated and the ∕ criterion, which assesses the turbulence equilibrium, is selected. Starting from a RANS computation for initialization, successive HTLES are carried out using the results of the previous computation to generate the mesh of the next simulation.

The outcomes of the study show that the self-adaptive strategy converges in four iterations. The statistically-averaged profiles and the friction coefficient are significantly improved when compared to the RANS computation and are in good agreement with the reference DNS results. The assessment of the components of the Reynolds stress anisotropy tensor emphasizes the reliability of the active HTLES and its benefits when associated with the proposed adaptive approach. The combination of the quite accurate prediction of the turbulent kinetic energy, the reliable estimate of the turbulence anisotropy, and the conformity of the spectra expresses the physical consistency of the turbulent structures generated. Downstream of the LES region, the turbulent structures are only gradually dissipated by the RANS eddy viscosity as observed on the Q-isosurfaces and power spectral densities.

Although many parameters can be adjusted and/or generalized, this study demonstrates that a self-adaptive strategy can be used effectively. As the algorithm is initialized from a RANS calculation, the user's influence on the results is reduced to the same level as for a RANS calculation. The user's subjective choices, based on his intuition and a priori knowledge of the flow physics, may lead to very different zones being defined depending on the user, as illustrated in Fig. 1. Obviously, the criteria used in the self-adaptive algorithm have a definite influence on the results, but if these criteria are fixed once and for all, the results are reproducible and independent of the user, which is of great importance for industrial applications.

Since the method has been developed validated on a specific flow, future work must be devoted to the application of the self-adaptive strategy to other test cases, possibly with the application of additional physical criteria, intending to develop a robust methodology useful for a wide range of applications.

Figure 1 :

 1 Figure 1: Schematic view of the problem (here non-conformal grid are represented for the sake of clarity). Top: Usual user-dependent hybrid RANS/LES. Bottom: Proposed self-adaptive strategy.

Figure 2 :

 2 Figure 2: Algorithm used for the self-adaptive approach.

Figure 3 :

 3 Figure 3: Isocontours of some candidate physical criteria for defining the interface between RANS and LES.

Figure 4 :

 4 Figure 4: Visualization of the region where ∕ > based the RANS computation performed with the -SST model. The two bounding boxes delimit zones A, B, and C corresponding to the RANS, RANS/LES transition, and LES zones, respectively.

  1, 1 and 2 the blending functions, and 1 , 2 , 1 , and 2 are the usual coefficients of the k-SST model. m is the turbulent time scale that drives the transition from the RANS mode to the LES mode m

Figure 5 :

 5 Figure 5: Q-isosurfaces colored by the velocity magnitude ( = 0.5 2 ∕ ) and targeted turbulent energy ratio.

Figure 6 :

 6 Figure6: velocity profiles in the direction for 12 streamwise locations. DNS results are from Lamballais[START_REF] Lamballais | Direct numerical simulation of a turbulent flow in a rotating channel with a sudden expansion[END_REF].

Figure 7 :

 7 Figure 7: Visualization of the mesh, the regions where ∕ > and the bounding boxes, for the RANS and HTLES1 computations. Zones A, B, and C correspond to the RANS, transition, and LES zones, respectively. The red line on the mesh view corresponds to the plane ∕ = 0.2. On the right side, the mesh in this plane is displayed. The values of ∕ 2 indicated are the non-dimensional area of the region identified by the criterion.

Figure 8 :

 8 Figure 8: Same as Fig. 7 for HTLES2 and HTLES3.

Figure 9 :

 9 Figure 9: Top: number of cells and averaged value of Δ in the region where ∕ > . Bottom: area of Zone C. Note that HTLES5 has not been run, but the mesh, built with the results of HTLES4, is displayed to highlight the convergence of the algorithm.

  (a) Energy ratio (b) Effective energy ratio (c) Streamwise velocity (d) Total turbulent kinetic energy

Figure 10 :

 10 Figure 10: Profiles in the direction for 12 streamwise locations. DNS results are from Lamballais [14].

Figure 11 :

 11 Figure 11: Profiles of the Reynolds stresses (left) and anisotropy components (right) in the direction for ∕ = 0 (top) and ∕ = 1.5 (bottom). The white, light green and medium green zones show the RANS, transition, and LES zones, respectively. DNS results are from Lamballais [14].

Figure 12 :

 12 Figure12: Distribution of skin friction on the wall downstream of the step. DNS results are from Lamballais[START_REF] Lamballais | Direct numerical simulation of a turbulent flow in a rotating channel with a sudden expansion[END_REF].

Figure 13 :

 13 Figure 13: Mean streamlines of the HTLES3 colored with velocity magnitude.

Figure 14 :

 14 Figure 14: Power spectral densities of the streamwise, wallnormal, and spanwise velocities obtained with HTLES3, upstream of the step, at ∕ = 0.557 (corresponding to + ≈ 84).The results are normalized with the velocity, ℎ = 3 ∕2, and the frequency step, Δ . The cutoff frequency, = 50 Hz, is higher than the maximum plotted frequency.

Figure 15 :

 15 Figure 15: Power spectral densities obtained with HTLES3 for ∕ = 1∕3, ∕ = 0.2, and ∕ = 0.04 (corresponding to + ≈ 84, + ≈ 50, and + ≈ 10, respectively), from left to right. The results are normalized with the velocity, , and the frequency step, Δ . The ranges in light gray color show frequencies higher than the cutoff frequency.

Figure 16 :

 16 Figure 16: Q-isosurfaces colored by the velocity magnitude ( = 0.5 2 ∕ 2 ) and targeted turbulent energy ratio for the HTLES3.

Table 1

 1 Coefficients of the HTLES---SST model.

	0	1	2	1	2
	0.48 2∕3 13.5 1.2 8	6

Acknowledgments

This study was funded by E2S UPPA (ANR-16-IDEX-0002) and INRIA in the framework of the ASTURIES project (grant E2S-20-ScientificChallenges-01). The computing resources were provided by MCIA (Mésocentre de Calcul Intensif Aquitain) and GENCI-IDRIS (Grants 2022-A0122A10980 and 2023-A0142A10980).

CRediT authorship contribution statement


Flora Clement
email: flora.clement@cea.fr

Monika Hospodiuk-Karwowski

Anastasia Papoz

Amandine Pitaval

Emily Tubbs

Vincent Haguet

Nathalie Bertacchi

Stéphanie Combe

Xavier Gidrol

Characterizing biomarkers for early diagnosis of pancreatic cancer

The pancreas, essential for both endocrine and exocrine functions, is associated with serious conditions such as diabetes mellitus (DM) and pancreatic cancer (PC). PC is one of the most lethal tumors, largely due to late-stage diagnoses [1]. Individuals with DM, particularly those with new-onset diabetes, face an elevated risk of PC, yet the molecular mechanisms remain unclear [2].

Our hypothesis posits that a diabetic microenvironment facilitates the initiation of PC [3]. The project aims to develop an amphicrine pancreas-on-chip for investigating the links between DM and PC. This chip comprises a chamber with a duct lined by exocrine cells, surrounded by islets of Langerhans in hydrogel (A). 3D cultured islets and various hydrogels enable perfusion for extended studies (B-D). Monitoring involves multi-OMICs approaches (secretome analysis, proteomics, transcriptomics), and 3D imaging.

We have already established two distinct devices (D and E), and our preliminary results demonstrate effective modeling of endocrine function using various islets of Langerhans models (F-J).

Figure 1 Amphicrine pancreas-on-chip. A: Schematic view of the device: a channel is generated through the use of a sacrificial 3D printed ink. It is then surrounded by hydrogel containing islets of Langerhans. After sacrificial ink removal, the channel is lined by exocrine cells. The tumoral initiation will be visualize thanks to 3D imaging including Light Sheet Fluorescence Microscopy available in the lab (3D reconstructed samples have already been performed). B: 4 devices paralleled in the incubator. C: one device with a zoom on a channel lined by exocrine cells. D: perfusion flow from the entrance needle to the exit hole. E: new version of 4 chambers included in one device, with peristaltic pump and a collection box. This new version follows the ISO compatibility norms of medical devices. Different types of islet of Langerhans models available in the lab: F: native islet from deceased donor, G: Langerhanoids made from the reaggregation of native islet dissociation supplemented in stromal and endothelial cells[4], H: spheroid of EndoC-βH5 cells, I: pseudo-islet made from the aggregation of EndoC-βH5 cells supplemented in stromal and endothelial cells, J: commercially available InSphero® pseudo-islet.

The objectives encompass refining the device and modeling the diabetic microenvironment by using (i) a decellularized matrix obtained from adipose tissue of both healthy and diabetic patients, (ii) diabetic islets, (iii) diabetic adipose stem cells, and exploring neoplastic transformation for early diagnostic biomarkers. Anticipated outcomes involve gaining novel molecular insights into the links between DM and PC, as well as identifying early PC biomarkers which are desperately needed. The microphysiological device developed will make a substantial contribution to comprehending tumor initiation PC. Furthermore, it has the potential to extend beyond pancreas, serving as a valuable tool for studying disorders in other glands.


C Ghouila-Houri

T Arnoult

A Mazzamurro

S Kern

R Viard

D Teillet

E Palfi

E Garnier

A Talbi

P Pernod

Flight-testing of a MEMS wall shear stress thermal sensor on a microlight aircraft

Keywords: MEMS thermal sensor, Wall shear stress sensor, Flight-testing I

This paper reports the realization of thermal MEMS (Micro Electro-Mechanical Systems) sensors designed for wall shear stress measurement during real flight on a microlight aircraft. Based on constant temperature anemometry and calorimetric principle, the micro-sensors measure the viscous shear stress occurring at the wall of a structure interacting with a fluid flow and detect the flow direction. A proper packaging was realized for ensuring a flush but robust integration of the devices on a pod fixed on a microlight aircraft. Real flight-testing demonstrated the ability of the micro-sensors to realize aerodynamic measurements up to 250 km/h on the microlight aircraft.

INTRODUCTION

Aircraft instrumentation is mainly composed of average measurement sensors indicating the mean velocity, the mean altitude the mean exterior temperature and so on, of the aircraft. However, the localized values of velocity and wall shear stress on the plane is of great interest, particularly for flow control applications. Flow control aims at locally manipulating the turbulent flow field to induce favorable changes in the main performances such as decreasing the drag or the noise radiation. Controlling turbulence then requires the measurement of its localized properties and the today instrumentation of aircraft is not sufficient. A network of small sensors, placed strategically on the plane, could give important information on the interaction of the main flow field with the structure.

The measurement of aerodynamic flow led to the development of micromachined sensors based on thermal transduction. The basic principle consists in electrically heating a metallic or semiconductor wire and measuring the heat losses due to forced convection with the fluid flow. Two main kinds of sensors were developed: hot-wires and hotfilms. In hot-wires, the wire is free in the flow field and only fixed on both extremities ([1], [START_REF] Bailey | Turbulence measurements using a nanoscale thermal anemometry probe[END_REF]). Suitable for aerodynamic measurement in wind tunnel environment even in highly turbulent flows [START_REF] Vallikivi | Turbulence measurements at high Reynolds numbers using a new inclined Nano-Scale Thermal Anemometry Probe[END_REF], these sensors are too fragile to be implemented on an aircraft. In hot-films, the wire is deposited on a film or a membrane. More robust, these sensors can be flush mounted on the surface ([4], [START_REF] Miau | Mems thermal film sensors for unsteady flow measurement[END_REF]) but the measured heat loss is a combination of conduction in the wall and forced convection with the flow. Since they are bot robust and non-invasive (flush-mounted), hot-film probes are widely used for measuring the wall shear stress on models ([6]- [START_REF] Tabrizian | An experimental study on boundary layer transition detection over a pitching supercritical airfoil using hot-film sensors[END_REF]).

The integration and use of thermal micro-sensors on real vehicles is another step of technology maturity. The literature gives few examples of thermal sensor embedded on real vehicles. Demonstrations were provided on a unmanned aerial vehicle [START_REF] Xu | Flexible shear-stress sensor skin and its application to unmanned aerial vehicles[END_REF] and on an airplane [START_REF] Hausmann | Coated Hot-Film Sensors for Transition Detection in Cruise Flight[END_REF] both using coated hot-film probes.

Developed as a compromise between hot-films and hotwires and with the objective of real flight application, the high temperature gradient micro-sensors designed at IEMN is composed of wall suspended wires, mechanically supported by periodic perpendicular micro-bridges ( [START_REF] Ghouila-Houri | High temperature gradient microsensor for wall shear stress and flow direction measurements[END_REF], [START_REF] Ghouila-Houri | High temperature gradient calorimetric wall shear stress micro-sensor for flow separation detection[END_REF]). The sensors were successfully tested in various configurations: in low [START_REF] Ghouila-Houri | Unsteady flows measurements using a calorimetric wall shear stress micro-sensor[END_REF] and high-speed wind tunnels [START_REF] Ghouila-Houri | MEMS high temperature gradient sensor for skin friction measurements in highly turbulent flows[END_REF] and on models [START_REF] Ghouila-Houri | High temperature gradient microsensors array for flow separation detection and control[END_REF]. In this paper, we present the utilization of such micro-sensors for real flight experiments realized on a microlight aircraft. The first part of the paper will describe the design of the sensor and the packaging developed for ensuring the flush and robust integration of the chip on a pod fixed on the microlight aircraft. The second part will present the setup on the microlight aircraft and discuss the results obtained during the flight campaign.

II. DESIGN AND PACKAGING OF THE SENSOR

A. Design of the thermal micro-sensor

The thermal micro-sensor developed for wall shear stress measurement is designed as a compromise between fragile but highly sensitive hot-wire probes and robust but slower hot-film probes. Indeed, the micro-sensors is composed of 3 free-standing wires mechanically supported by perpendicular SiO2 based micro-bridges. The wires are separated from the substrate by a 20 µm deep cavity ensuring a high thermal insulation for a hot-wire like behavior. The micro-bridges provide a mechanical robustness of the structure without enhancing thermal loss by conduction The micro-sensor is designed to work simultaneously as a hot-wire probe using the middle wire and as a calorimetric sensor using the difference of temperature between the two lateral wires. Calorimetric sensors present the advantage of a differential measurement that allows measuring both the amplitude and the sign of the wall shear stress.

B. Packaging for flight experiments

For the previous experiments done on low and high speed wind tunnels and on models, a flexible PCB packaging was realized to have a mounting robust and as flush as possible, as shown on Fig. 2 (a). However, the sensor chip has front side electrical contacts meaning that the sensor chip was not perfectly flush mounted, but 100 µm below the surface of the PCB. Nevertheless, this integration was the best suited for the model experiments on which the area of integration was limited and a flexible solution was require.

Since the flight experiments allowed more surface for the sensor integration, we improved the packaging by using a rigid PCB, as shown on Fig. 2

III. FLIGHT EXPERIMENTS

A. Set-up on the microlight aircraft

Real flights tests were performed using a Flight Design CTSW microlight aircraft. It is equipped with a 3 blades propeller and a Rotax 912s engine, characterized by a gear box ratio of 1/2.43. PODs are set up on the airplane landing gear. One is instrumented with a calorimetric micro sensor, as illustrated in Fig. 4. The POD to ground communication is realized using a radio module embedded inside the POD along with the conditioning electronics of the calorimetric micro-sensor.

B. Results and discussion

The experiments were conducted at 1500 ft of altitude. Four speed steps were investigated: 100 km/h, 150 km/h, 150km/h and 250km/h. For the last step of speed, the pilot encountered difficulties maintaining the 250 km/h and consequently this step has a velocity comprised between 230 km/h and 250 km/h. The series of velocity steps was realized twice. Fig. 5 presents the results obtained during the flight experiments: it displays the output voltage corresponding to the difference of voltage from the two lateral wires, as a function of time. Several phases are distinguishable on the measurements. First, up to about 600s, the microlight aircraft is at rest with no velocity. From 600s to about 800s corresponds the taxiing phase with several turns done by the airplane. Then the sudden increase corresponds to take off. After that, the velocity steps study began. As shown on Fig. 5, the calorimetric micro-sensor detected the steps of velocity with a good repeatability since the experiments were conducted twice. In each case, the four steps are easily identifiable.

The time resolution of the calorimetric micro-sensor enabled also the detection of the propeller signature. Indeed, Fig. 6 displays the power spectral density of the signal at rest and on each airspeed step. The measured fundamental frequencies are 76 Hz, 87.5 Hz, 92 Hz and 106 Hz respectively for 100 km/h, 150 km/h, 200 km/h and 250 km/h. These frequencies are in agreement with the engine constructor documentation. Harmonics of these frequencies can also be observed in the spectra.

Fig. 1 .

 1 Fig. 1. Scanning electron microcopy image of the sensitive part of calorimetric micro-sensor unlike hot-film probes. Fig. 1 shows a Scanning Electron Microscopy (SEM) image of the sensor sensitive part. The central wire is the heater. Made of Ti/Au it is heated by an Joules effect. The two lateral wires are Ni/Pt based coldwires measuring the temperature on both sides of the heater.The micro-sensor is designed to work simultaneously as a hot-wire probe using the middle wire and as a calorimetric sensor using the difference of temperature between the two lateral wires. Calorimetric sensors present the advantage of a

  (b). Front side contacts were realized using conductive epoxy. Wind tunnel facilities in ONERA Lille allowed testing the operation of the micro-sensors before real flight experiments. The micro-sensors in both packaging were integrated in a turbulent boundary layer wind tunnel with a fluid velocity going up to 40 m/s. The novel packaging demonstrated a high sensitivity to wall shear stress, for mainstream velocities going up to 40 m/s, as shown in Fig.3. Compared to what was obtained with the flexible PCB flip chip integration, the sensitivity is about 10 times improved.

Fig. 2 .Fig. 3 .

 23 Fig. 2. (a) Flexible flip-chip packaging used for integrationg the calroimetric micro-sensors in wind tunnel facilities and on models (b) Rigid PCB based packaging developed for the microlight aircraft flight experiments

Fig. 4 .

 4 Fig. 4. Set up o fthe microlight aircraft flight experiments: phtograph of the Flight Design CTSW microlight aircraft with a zoom on the PODs and on the calorimetric micro-sensor chip

Fig. 5 .Fig. 6 .

 56 Fig. 5. Voltage difference of the two lateral wire as a function of time during the real flight experiments with identification of the different moments of the flight.

Authorized licensed use limited to: UNIVERSITE DE LILLE. Downloaded on January 16,2024 at 09:44:24 UTC from IEEE Xplore. Restrictions apply.

ACKNOWLEDGMENT

This work was funded by the French National Research Agency (ANR) in the framework of ANR ASTRID MATURATION "CAMELOTT-MATVAL" project. It is also supported by the CPER ELSAT 2020 project and the Innovation Defense Agency (AID). The authors also thank RENATECH, the French national nanofabrication network, and FEDER.


MD Christophe Milesi
email: c-milesi@chu-montpellier.fr

Erika Nogue

Julien Baleine

Lionel Moulis

Robin Pouyau

Arthur Gavotto

David Brosssier

Guillaume Mortamet

Gilles Cambonie

Erika Nogue2

| Julien

Baleine Md1

| Robin

Pouyau Phd3

| Arthur

Gavotto Phd1

| David

Brosssier Phd4

| Guillaume

Mortamet Phd5

MD, PhD1 | GFRUP Respiratory Study Group | Gilles Cambonie

ROX (Respiratory rate-OXygenation) index to predict early response to high-flow nasal cannula therapy in infants with viral bronchiolitis

Keywords: bronchiolitis, failure, high-flow nasal cannula, infant, risk prediction AUC, area under the curve, AVB, acute viral bronchiolitis, CPAP, continuous positive airway pressure, EDIN, [Échelle de douleur et d

ROX (Respiratory rate-OXygenation) index to predict early response to high-flow nasal cannula therapy in infants with viral bronchiolitis

Result:

HFNC failure occurred in 111 of 286 (39%) infants, and for 56 (50% of the failure) of them within the first 6 h. The area under the curve of ROX indexes at H0 and H1 were, respectively, 0.56 (95% confidence interval [CI] 0.48-0.63, p = 0.14), 0.56 (95% CI 0.49-0.64, p = 0.09). ROX-HR performances were better but remained poorly discriminant. HFNC failure was associated with higher mWCAS score at H1 (p < 0.01) and lower decrease in EDIN scale during the first hour of HFNC delivery (p = 0.02). In the multivariate analyses, age and mWCAS score were found to be independent factors associated with HFNC failure at H0. At H1, weight and mWCAS were associated factors.

Conclusion:

In this study, neither ROX index, nor physiological variables usually collected in infants with AVB had early discriminatory capacity to predict HFNC failure.

Acute viral bronchiolitis (AVB) is a leading cause of lower respiratory tract infection and hospitalization in young infants, notably those aged less than 1 year.1 Management of moderate-to-severe AVB is based on respiratory support, with as the first step administration of a heated and humidified mixture of air and oxygen with high-flow nasal cannula (HFNC).2 While several studies have confirmed HFNC effectiveness to reduce the work of breathing in AVB,3-5 failure occurs in 30%-40% of these patients.6-8 In these infants, evolving respiratory failure requires escalation in therapeutic measures, including transition to continuous positive airway pressure (CPAP), then if necessary to noninvasive ventilation or intubation.9

Early identification of patients who are most likely to fail with HFNC is critical for care organization in the pediatric emergency department.10 Indeed, these infants will need to be referred to a pediatric intensive care unit, necessary for the monitoring of any noninvasive ventilation technique including CPAP and bi-level positive airway pressure, while those who will improve with HFNC can potentially be transferred to a general pediatric ward.11 Currently, some patient characteristics have been individualized as predictors of respiratory deterioration, notably younger age or initial severity.12,13 However, neither isolated physiologic parameters, such as respiratory rate (RR), fraction of inspired oxygen (FiO2), or venous/capillary partial pressure of carbon dioxide (pvCO2), nor clinical scales that incorporate different vital signs have demonstrated a consistent association with the risk of HFNC failure and are discriminating enough to be used as triage tools. [START_REF] Kelly | High-flow nasal cannula use in children with respiratory distress in the emergency department: predicting the need for subsequent intubation[END_REF][START_REF] Betters | High flow nasal cannula use outside of the ICU; factors associated with failure[END_REF][START_REF] Er | Early predictors of unresponsiveness to high-flow nasal cannula therapy in a pediatric emergency department[END_REF][START_REF] Santis | Uses of highflow nasal cannula on the community paediatric ward and risk factors for deterioration[END_REF][START_REF] Suessman | Clinical factors associated with intubation in the high flow nasal cannula era[END_REF] Recently, Roca et al. developed a tool to assess the risk for HFNC failure in adult patients with hypoxemic acute respiratory failure.19 ROX (Respiratory rate-OXygenation) index corresponds to the ratio of patient oxygenation, which has been associated to HFNC success, over RR, which has been associated to HFNC failure. Subsequent studies confirmed it was a good predictor of HFNC failure in lower respiratory tract infections, including those caused by virus.20,21 In pediatric patients with acute respiratory failure, ROX index application 24 and 48 h after hospital admission also appeared a good marker for predicting the risk of HFNC failure. [START_REF] Yildizdas | Predicting nasal high-flow therapy failure by pediatric respiratory rate-oxygenation index and pediatric respiratory rate-oxygenation index variation in children[END_REF] In the specific context of AVB in <2 years infants managed in a pediatric emergency department, patients with ROX index in the lowest quartile at HFNC initiation were three times more likely to require CPAP compared to those in the highest quartile. [START_REF] Kannikeswaran | Association between respiratory rate oxygenation index and need for positive pressure ventilation in children on high flow nasal cannula for bronchiolitis[END_REF] The single center nature of this study, with lack of standardization for failure criteria, incited to test the relevance of ROX index in an homogeneous population of <6 months infants requiring HFNC for severe AVB, recruited in the framework of a multicenter study with predefined HFNC failure criteria.24 ROX-HR index is defined as the ratio of ROX index over HR (beats/min), multiplied by a factor of 100. This index could be a useful tool for early prediction of HFNC outcomes.25

The primary objective of this study was to assess the relationship between ROX index collected early (i.e., before HFNC initiation and 1 h after), and HFNC failure occurring in the following 48 h in patients admitted for severe AVB. This is an ancillary study to the multicenter randomized controlled trial TRAMONTANE 2, which compared two settings (2 vs. 3 L/kg/ min), in young infants with severe AVB (defined by a modified Wood's Clinical Asthma Score [mWCAS] >3) supported with HFNC. [START_REF] Milési | A multicenter randomized controlled trial of a 3-L/kg/min versus 2-L/kg/min high-flow nasal cannula flow rate in young infants with severe viral bronchiolitis (TRAMONTANE 2)[END_REF] The study protocol was approved by the South Mediterranean IV Ethics Committee (2016-A00900-51), and recorded on the National Library of Medicine registry (NCT02824744). Written authorization was obtained from the two parents.

The study, performed in 16 French university hospital centers between November 2016 and March 2017, enrolled 142 infants in the 2 L/min/kg group, 144 infants in 3 L/min/kg group, and found the same failure rate of 39% in each group, allowing to include the 286 TRAMONTANE 2 patients in the current analysis.

| Data collection

Heart rate (HR), RR, FiO2, venous/capillary pH and pvCO2, mWCAS, and the neonatal pain and discomfort scale (EDIN) were collected at baseline (H0), in the 15 min following admission, while the infant received a blend of air/oxygen using nasal cannula at a maximal flow rate of 1 L/min. The same physiological variables, except for capillary blood gas analysis, were collected 1 h after HFNC initiation at the flow rate allocated by randomization (H1).

At H0 and H1, ROX and ROX-HR were calculated as and respectively,

Predefined HFNC failure criteria included the occurrence, within 48 h of randomization, of one or more of the following: increase, compared to H0, in respiratory distress score (i.e., 1 point in mWCAS) or RR (i.e., >10 breath per min with RR >60 bpm), increase in discomfort (1 point in EDIN score), and severe apnea episodes requiring bag and mask ventilation.

| Outcome

The primary study outcome was discriminatory capacity of ROX and ROX-HR index values, at H0 and H1, and their change during the first hour of HFNC delivery (H1 -H0), to predict HFNC failure. We also tested, as secondary outcomes, the physiological variables collected at H0 and H1 and their changes between the two time points. We also analyzed predictive factors of HFNC failure among sociodemographic and clinical data.

| Statistical analyses

Patients' characteristics were presented using median and interquartile range for quantitative variables and frequencies with proportions for categorical ones.

Groups, defined as HFNC failure or success, were compared using Student or Wilcoxon Mann-Whitney test for continuous variables and Chi-square or Fisher test for categorical ones.

The accuracies of ROX index, ROX-HR index, and physiological variables to predict HFNC failure were assessed using receiver operating curve analyses. The area under the curve (AUC) was calculated by the Hanley method and compared to the 0.

| RESULTS

HFNC failure occurred in 111 of 286 (39%) infants, after a median time (Q25-Q75) of 6.25 (1.67-17.00) hours (Figure 1). In 96 of 111 (86%), failure occurred within 24 h, including 56 (50%) within the first 6 h.

| Primary outcome

ROX and ROX-HR indexes at H0 and H1, and their changes during the first hour were not different according to HFNC failure or success (Table 1). Their respective AUC were presented in in Table 2. These results did not support that early ROX indexes predicted HFNC failure. Only the ROX-HR at H1 appeared significantly associated with HFNC failure but the performance associated were poorly (Figure 2).

| Secondary outcome

In univariate regression, HFNC failure was associated with higher mWCAS score at H1 (p < 0.01) and lower decrease in EDIN scale during the first hour of HFNC delivery (p = 0.02) (Table 1). However, none of the physiological variables collected at H0 and H1 and their changes during the first hour of HFNC delivery were predictive of HFNC failure, according to their respective AUC (Table 2). In the multivariate analysis at H0, age was identified as a protective factors against HFNC failure, whereas mWCAS score was found to be a risk factor (Table 3). At H1, weight was found to be protective against HFNC failure, while mWCAS at H1 and EDIN change were retained in the model as potential risk factors (with EDIN change being nonsignificant).

In this multicenter population of young infants with severe AVB, the failure rate of HFNC was 39%. Failure of this first step of respiratory support occurred early, within 6 h after HFNC initiation for half of failing patients and within 24 h for nearly 90% of them. Neither ROX, ROX-HR index, nor physiological variables usually collected in infants with acute respiratory failure had early discriminatory capacity to predict failure of management with HFNC.

In a previous study, Kannikeswaran et al. observed The ROX index takes into account only two characteristics of a respiratory distress, namely oxygenation and tachypnea, plus HR for the ROX-HR. While these parameters are critical in patients with acute hypoxemic respiratory failure, they do not take into account all of the determinants of HFNC failure in AVB. Indeed, AVB present different phenotypes: sometimes as a restrictive parenchymal disease, but a majority of these infants demonstrate a severe obstructive lung disease, with markedly increased work of breathing and frequent apneas.4,5,26 These two elements, as well as comfort, are not integrated in the ROX index, whereas they often intervene in the clinician's decision to upgrade respiratory support in a patient with AVB.

No variable, observed or calculated, was able to predict HFNC failure in this work. However, half of failure occurred within 6 h, which suggests that such delay may be a relevant criterion for triage. Insofar as the volume of patients in the emergency room allows it, it may be consistent to maintain the infant with HFNC and adapted monitoring during this timeframe before deciding admission to pediatric intensive care unit.

This study has several strengths, including a large, multicentric and homogeneous population, in terms of severity and age, of infants with AVB. In addition, predefined HFNC failure criteria had been validated by a panel of experts providing from the 16 participating centers in the trial.

We acknowledge some limitations. Half of the children received a flow rate of 2 L/min/kg, and the other half 3 L/min/kg, but we are uncertain whether this affects the interpretation of our results because the failure rate was comparable between groups: 38.7% (2 L/min/kg) vs. 38.9% (3 L/kg/min; p = 0.98). Some data points for the ROX index at baseline and 1 h after HFNC initiation are missing.

The predictive value of the ROX index beyond the first hour of HFNC delivery could have been evaluated, but it seemed more obvious that the relevance of an urgent triage tool can hardly exceed 1-2h.

  5 value using Wilcoxon's W statistic. For the analyses at H1 and for changes during the first hour of HFNC delivery, all HFNC failure patients at H1 or before were excluded. To analyze risk factors of HFNC failure, we used a multivariate logistic regression. A backward model selection of Akaike's information criterion was implemented on variable with a p < 0.20 in univariate regression. A model was performed at H0 and another at H1. All statistical tests were two-tailed with an alpha set at 0.05. Statistical analyses were conducted with SAS (version 9.04; SAS Institute) and R software (R 4.3.1).

  an association between ROX index and HFNC failure, in other words the need for positive pressure ventilation, in <2 years infant with bronchiolitis.23 This result was particularly relevant from the perspective of directing these patients to the most suitable units downstream of the emergency department. In this perspective, our study aimed to identify a threshold for this index, associated with actual predictive capacities. Our AUC results suggested a weak and nonsignificant relationship between ROX or ROX-HR and HFNC failure, which currently does not confirm the interest of this tool in a clinical decision rule. The main difference with studies suggesting that ROX index may be a good marker to predict the risk of HFNC failure probably comes from different patient's characteristics.22,23 Our trial involved much younger patients and, as our results indicate, age is a key risk factor for respiratory failure in this population.12,13 In addition, the TRAMONTANE 2 study included patients probably affected by more severe forms of the disease, as mWCAS ≥ 3 was required to be eligible, which signals unambiguous respiratory distress. In Kannikeswaran et al. study, the regression model to estimate the odds ratio of PPV requirement was based on the highest ROX quartile, suggesting marked heterogeneity in the severity of bronchiolitis.[START_REF] Kannikeswaran | Association between respiratory rate oxygenation index and need for positive pressure ventilation in children on high flow nasal cannula for bronchiolitis[END_REF]

| CONCLUSION

In this population of patient with severe AVB supported with HFNC, half of failures occurred within the first 6 h. Use of ROX or ROX-HR index before or 1 h after delivery of this support was not able to predict this event.

ACKNOWLEDGMENTS

We want to thank the members of the GFRUP ventilation section: Anne FLorence Pierre; Anna Deho; Jean Michel Liet; Camille Guillot; Anne-Sophie Guilbert; Jerome Rambaud; Astrid Millet; Mickael Afaneti; Julie Guichoux; Mathieu Genuini; Thierry Mansir; Jean Bergounioux; Fabrice Michel: Marie-Odile Marcoux; Sabine Durand; Philippe Durand; Stéphane Dauger; Etienne Javouhey; Stéphane Leteurtre; Olivier Brissaud; Sylvain Renolleau, Aurélie Portefaix; Aymeric Douillard.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

CONFLICT OF INTEREST STATEMENT

The authors delcare no conflict of interest.

Funding information

None