Christophe Friess
email: friess@l3m.univ-mrs.fr

Ébastien Poncet

Éphane St

Viazzo

AN EXTENSIVE STUDY ON LES, RANS AND HYBRID RANS/LES SIMULATION OF A NARROW-GAP OPEN TAYLOR-COUETTE FLOW

Keywords: Simulation, EB-RSM Elliptic Blending Reynolds Stress Model, LES Large-Eddy Simulation, RANS Reynolds-Averaged

The present paper concerns a numerical benchmark of various turbulence modelings, from RANS to LES, applied to Taylor-Couette-Poiseuille flows in a narrow gap cavity for six different combinations of rotational and axial Reynolds numbers. Two sets of refined Large-Eddy Simulation results, using the WALE and the Dynamic Smagorinsky subgrid scale models available within an in-house code based on high-order compact schemes, hold for reference data. The efficiency of a RANS model, the Elliptic Blending Reynolds Stress Model (EB-RSM) [1], and a hybrid RANS/LES method, the so-called "Equivalent DES" [2], both run with Code Saturne, is then questioned. Thin coherent structures appearing as negative (resp. positive) spiral rolls are obtained by the LES but also the hybrid RANS/LES along the rotor (resp. stator) sides. More quantitatively, the hybrid RANS/LES does not improve the predictions of the EB-RSM for both the mean and turbulent fields, stressing the need for further theoretical development. NOMENCLATURE C, C m model parameters. C ε1 ,C ε2 model parameters. δt time step. h length of the cylinders. k total (resolved+modeled) turbulence kinetic energy. k m modeled part of the turbulence kinetic energy.

INTRODUCTION

Since Taylor's first theoretical results obtained one century ago, huge research efforts have been done to better describe and understand the rich dynamical behavior exhibited by the flow induced by the differential rotation of two concentric cylinders. This kind of flow has found many applications in chemical engineering or in the turbomachinery industry for examples. In the present case, the main motivation arises from the cooling of electrical motors, which can be modeled quite faithfully by considering a narrow-gap Taylor-Couette system with an axial Poiseuille flow. A better knowledge of the hydrodynamic field is absolutely necessary before considering the heat transfer aspect as pointed out recently by Fénot et al. [START_REF] Fénot | A review of heat transfer between concentric rotating cylinders with or without axial flow[END_REF]. These last authors provided a very useful and exhaustive literature survey on Taylor-Couette flows with or without an axial Poiseuille flow and including or not heat transfer processes. They explained the difficulty to establish universal correlations for the heat transfer coefficients by the large numbers of parameters involved in the problem. Moreover, most of the previous studies were mainly concerned with temperature measurements without any idea of what the hydrodynamic flow was and especially what the inlet flow conditions were. In the present work, one proposes a numerical benchmark for various combinations of the flow parameters in the isothermal case. It is a step forward the turbulence modeling of the full problem with heat transfer.

Nouri and Whitelaw [START_REF] Nouri | Flow of Newtonian and non-Newtonian fluids in a concentric annulus with rotation of the inner cylinder[END_REF] then Escudier and Gouldson [START_REF] Escudier | Concentric annular flow with centerbody rotation of a newtonian and a shearthinning liquid[END_REF] provided very useful experimental databases for middle-gap cavities (η ≃ 0.5) of large aspect ratios (Γ = 98 and 244 respectively). Until now, most numerical studies used turbulence modeling, which provided rather limited informations. It was shown by [START_REF] Poncet | Numerical modeling of fluid flow and heat transfer in a narrow Taylor-Couette-Poiseuille system[END_REF][START_REF] Rothe | The influence of rotation on turbulent flow and heat transfer in an annulus between independently rotating tubes[END_REF] that increasing the rotation rate of the inner cylinder amplified the turbulence kinetic energy, resulting in the enhancement of heat transfer along the rotor. More recently, the development of computational methods (including Direct and Large Eddy Simulations) has led to an increase in numerical studies of rotating flows, but few works were concerned by turbulent Taylor-Couette-Poiseuille flows and especially the near-wall turbulent structures in such systems. As example, in the middlegap configuration, Chung and Sung [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF] established the destabilization of the near-wall turbulent structures due to rotation of the inner wall giving rise to an increase of sweep and ejection events. In the same way, Hadziabdic et al. [START_REF] Hadziabdic | LES of turbulent flow in a concentric annulus with rotating outer wall[END_REF] studied by LES a fully-developed turbulent flow in a concentric annulus of radius ratio η = 0.5, with the outer wall rotating at a range of rotation rates N = Re Ω /Re Q from 0.5 up to 4. They focused their attention on the effect of the rotation parameter N on the turbulence statistics and coherent structures in the near-wall regions. To our knowledge, there is no reference numerical or experimental data for the narrow-gap case (η > 0.8). Thus, the present in-house code based on fourth-order compact schemes using the dynamic Smagorinsky subgrid scale model has been first validated by Oguic et al. [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF] against the experimental data of Nouri and Whitelaw [START_REF] Nouri | Flow of Newtonian and non-Newtonian fluids in a concentric annulus with rotation of the inner cylinder[END_REF]. Its slightly improved the LES results of Chung and Sung [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF] based on second-order numerical schemes highlighting the importance of the order of the spatial schemes. It will be considered here as providing the reference data to discuss the capability of RANS and hybrid RANS/LES of predicting the mean and turbulent flow fields in the narrow-gap case.

The present work is an extension of Friess et al. [START_REF] Friess | Taylor-Couette-Poiseuille flows : from RANS to LES[END_REF] to more combinations of the rotational and axial Reynolds numbers and to more turbulence modelings. Its goal remains twofold: (i) providing some reference LES data using two subgrid scale modes available within an in-house high-order solver and (ii) questioning the capabilities of a hybrid RANS / LES method, as well as the underlying RANS model, in predicting this kind of flow. The paper is organized as follows: the flow configuration is first described. The numerical approaches are then presented. Afterwards, the results about the hydrodynamic fields are discussed in details in terms of the coherent structures and the mean and turbulent flow fields, before some final conclusions.

FLOW CONFIGURATION AND PARAMETERS

The fluid is confined between two concentric cylinders of radii R 1 and R 2 and height h (see Figure 1). The inner cylinder rotates at a constant rate Ω, while the outer cylinder is stationary. The cavity is characterized by two geometrical parameters : its aspect ratio Γ = h/(R 2 -R 1 ) and its radius ratio η = R 1 /R 2 . An axial throughflow is imposed within the gap at a constant bulk velocity W m . The main flow parameters are the rotational Reynolds number Re Ω = ΩR 1 (R 2 -R 1 )/ν and the bulk Reynolds number Re Q = W m (R 2 -R 1 )/ν, ν being the fluid kinematic viscosity. The present study considers two values of the axial flow rate and three rotation rates of the inner cylinder, resulting in five distinct values of N = Re Ω /Re Q (see Table 1

NUMERICAL METHODS

Concerning LES reference calculations, the numerical method of the in-house code is based on the former work of Abide and Viazzo [START_REF] Abide | A 2D compact fourthorder projection decomposition method[END_REF], who devised a 2D compact fourth-order projection decomposition method in cartesian coordinates. It has been recently validated in the case of turbulent interdisk rotorstator flows in cylindrical coordinates by Viazzo et al. [START_REF] Viazzo | High-order Large Eddy Simulations of confined rotor-stator flows[END_REF] but also for turbulent Taylor-Couette-Poiseuille flows in a middlegap cavity by Oguic et al. [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF]. The time advancement is secondorder accurate and is based on the explicit Adams-Bashforth scheme for the convective terms and an implicit backward Euler scheme for the viscous terms. The derivatives are approximated using fourth-order compact formula in the radial and axial directions. The time splitting scheme is an improved projection method, which ensures the incompressibility at each time step. The projection decomposition method is based on a direct nonoverlapping multidomain Helmholtz/Poisson solver, which provides the solution of each Helmholtz/Poisson problem resulting of the velocity-pressure coupling. The multidomain solver ensures the continuity of the solution and its first normal derivative across the conforming interface. These continuity conditions are cleared up using an influence matrix technique. By assuming the periodicity of the solution in the azimuthal direction, the threedimensional extension uses Fourier series. The set of problems is thus reduced to series of two-dimensional problems associated with each wave number. Periodic boundary conditions are applied in the axial and circumferential (angle of 3π/4) directions and no-slip boundary conditions are imposed on the walls. The cavity is divided into 2 subdomains in the axial direction. All the numerical details including the mesh grid, the time step and the wall coordinates are provided in Table 2.

The CFD code used for RANS and hybrid RANS/LES simulations is the open-source Code Saturne, developed by EDF [START_REF] Archambeau | Code Saturne: A finite volume code for the computation of turbulent incompressible flows -Industrial applications[END_REF]. It is a finite volume solver, written in cartesian coordinates. A SIMPLEC algorithm, with the Rhie and Chow interpolation, is used for pressure-velocity coupling. Concerning hybrid RANS/LES calculations, convective fluxes are approximated by a second-order centered scheme, for momentum, and a first-order upwind scheme for subfilter quantities. Time marching uses a Crank-Nicholson, second-order scheme. The hybrid RANS/LES mesh contains 60 3 cells, with an angular periodicity of π/2 and an axial periodicity. No-slip boundary conditions are imposed on the walls. The RANS mesh is a 1D grid of 60 cells. In RANS mode, convective fluxes are approximated by a first-order upwind scheme, for all computed quantities, and the time scheme is first-order.

TURBULENCE MODELLING

Various levels of modeling are presented. Emphasis is first put on the hybrid RANS/LES method. Afterwards, the subfilter closure and RANS model will be briefly presented.

Hybrid Method

The hybrid method used for the present work is the "Equivalent DES" of Friess and Manceau [START_REF] Manceau | Of the interpretation of DES as a hybrid RANS / temporal LES method[END_REF] and Manceau et al. [START_REF] Friess | Investigation of the equivalence of two hybrid temporal-LES methods based on elliptic blending in separated flows[END_REF]. This approach was first derived for the purpose of bridging the PITM (Partially Integrated Transport Model) [START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF] and the DES (Detached Eddy Simulation) method. Indeed, the first is fully justified from a theoretical point of view, while the second one was developed on a rather phenomenological basis. The PITM was first generalized by Fadai-Ghotbi et al. [START_REF] Fadai-Ghotbi | Temporal filtering : A consistent formalism for seamless hybrid RANS/LES in inhomogeneous turbulence[END_REF] to inhomogeneous flows, considering temporal filtering, rather than spatial filtering. From a pragmatical point of view, considering RANS / TLES (temporal LES) hybridization does not cause any difficulty in implementing models or special terms, since Fadai-Ghotbi et al. showed that applying Temporal PITM (T-PITM) to an inhomogeneous stationary flow is equivalent to applying PITM to homogeneous, statistically unsteady flow. Actually, one must just keep in mind that "hybrid RANS/LES" is to be understood in a general way, including temporal LES.

The advantage of Detached Eddy Simulation lies in its simplicity and robustness. The idea of Manceau et al. [START_REF] Manceau | Of the interpretation of DES as a hybrid RANS / temporal LES method[END_REF] was thus to derive an approach bridging it with T-PITM. This latter is, in spite of its theoretical justification, not very easy to implement in any code. An equivalence criterion was then determined between T-PITM and DES, providing some theoretical justification to the latter, and allowing to interpret it as a hybrid RANS / TLES method. This equivalence criterion was derived analytically for equilibrium flows, but was successfully tested on a flow over a periodic hill, involving massive separation [START_REF] Friess | Investigation of the equivalence of two hybrid temporal-LES methods based on elliptic blending in separated flows[END_REF].

The principle of DES is to magnify the dissipation term of the transport equation for either subfilter turbulence kinetic energy k SFS or subfilter stresses τ i j SFS , damping the modeled en- ergy and allowing large scale eddies to be resolved:

ε DES i j = k 3/2 SFS ε SFS L ε i j SFS (1)
While classical DES uses the local grid step to determine the length scale L in Eq.( 1), equivalent DES uses the ratio r k defined as:

r k = k m /k,
where k is the total turbulence kinetic energy and k m its counterpart contained in the modeled (subfilter) scales. At the RANS limit, r tends to 1, and at the DNS limit, it tends to 0. The advantage of r k is that it can be estimated using a spatial or a temporal energy spectrum. As shown by Friess and Manceau [START_REF] Friess | Investigation of the equivalence of two hybrid temporal-LES methods based on elliptic blending in separated flows[END_REF], r k can be evaluated as :

r k = 1 k ∞ ω c E T (ω) dω = min 1, 1 β 0 U s √ k 2/3 ω c k ε -2/3 (2)
β 0 being a constant derived from the Kolmogorov constant (β 0 = 0.3 in the present case), U S a sweeping velocity and ω c the cutoff frequency of the considered filter. It can be defined as :

ω c = min π δt ; U s π ∆ (3)
where δt and ∆ are the time and grid steps, respectively. Finally, following [START_REF] Manceau | Of the interpretation of DES as a hybrid RANS / temporal LES method[END_REF], the length scale L entering Eq.( 1) is :

L = r 3/2 k 1 + C ε2 -C ε1 C ε1 (1 -r C ε1 /C ε2 ) k 3/2 ε (4)

RANS model and subfilter closure

In the temporally-filtered Navier-Stokes equations, the subfilter stresses must be modeled such a way that, making the temporal filter width go to infinity, the equations tend to the RANS equations. In particular, one of the main objectives of hybrid methods is to use RANS closures in the near-wall regions, to avoid the very fine resolution required by LES. In that aim, the RANS model proposed by Manceau and Hanjalic [START_REF] Manceau | Elliptic blending model: A new near-wall Reynolds-stress turbulence closure[END_REF] and Manceau [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF], the so-called elliptic blending Reynolds-Stress Model (EB-RSM), is adapted to the hybrid temporal LES context. In this model, an elliptic relaxation equation is solved for a scalar α :

α -L 2 SFS ∇ 2 α = 1 (5)
which is a sensor of the distance to the wall (α = 0 at the wall, and 1 far away), and is used to blend near-wall and homogeneous formulations for the redistribution and dissipation terms of the transport equations for subfilter stresses (see [START_REF] Fadai-Ghotbi | A seamless hybrid RANS/LES model based on transport equations for the subgrid stresses and elliptic blending[END_REF] for details) :

φ * i j -ε i j = α 3 (φ * i j -ε i j ) wall + (1 -α 3 )(φ * i j -ε i j ) homogeneous . (6)
Under its RANS form, EB-RSM was successfully applied to several flows, e.g. non-rotating and rotating channels [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF], and mixed and natural convection flows [START_REF] Dehoux | Modelling turbulent heat fluxes using the elliptic blending approach for natural convection[END_REF] among other things. It was recently adapted by Fadai-Ghotbi et al. [START_REF] Fadai-Ghotbi | A seamless hybrid RANS/LES model based on transport equations for the subgrid stresses and elliptic blending[END_REF] to serve as a subfilter-stress model in the framework of T-PITM. In the present work, it is also applied as a model for the equivalent DES, by simply substituting the dissipation in the fashion of Eq. ( 1) in the Reynolds-stress transport equations.

Subgrid-Scale Models

Reference data are here provided by Large Eddy Simulations, using two different subgrid-scale models implemented in the in-house code based on fourth-order compact schemes previously described. The Wall-Adapting Local Eddy Viscosity (WALE) model of Nicoud and Ducros [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF] has already been considered by Friess et al. [START_REF] Friess | Taylor-Couette-Poiseuille flows : from RANS to LES[END_REF] in the same configuration for particular values of the rotation parameter N. The study is extended to other values of N and to the use of the Dynamic Smagorinsky model developed by Germano et al. [START_REF] Germano | A dynamic subgrid-scale eddy-viscosity model[END_REF].

The Dynamic Smagorinsky model In the Smagorinsky model, the eddy viscosity is assumed to be proportional to the subgrid scale characteristic length ∆ = (r∆r∆θ ∆z) 1/3 and to a characteristic turbulent velocity taken as the local strain rate |S| = 2S i j S i j :

ν t = (C S ∆) 2 |S| (7)
with C S = 0.2. The classical Smagorinsky model appeared to be too dissipative in a previous work on turbulent rotor-stator flows [START_REF] Viazzo | High-order Large Eddy Simulations of confined rotor-stator flows[END_REF]. That is the reason why one uses it in its dynamical version, where C s is replaced by C d , which depends on both space and time. C d is evaluated with a least-squares approach as a part of the solution at each time step using a test filter equal to 2∆. C d is then averaged along the tangential direction (see in [START_REF] Germano | A dynamic subgrid-scale eddy-viscosity model[END_REF][START_REF] Lilly | A Proposed Modification of the Germano Subgrid Scale Closure Method[END_REF]).

The WALE model This choice is motivated by a better numerical stability than the dynamic Smagorinsky model, and its quality in near-wall treatment without damping functions. The subgrid viscosity is given by :

ν t = (C m ∆) 2 (S d i j S d i j ) 3/2 (S i j S i j ) 5/2 + (S d i j S d i j ) 5/4 , (8)
where C 2 m = 10.6C 2 s = 0.424 is constant, ∆ the grid step, S i j the filtered strain-rate tensor, and S d i j the deviatoric part of the squared filtered velocity gradient tensor. For the same mesh grid, the WALE model enables to save around 12% of the computational cost compared to the dynamic Smagorinsky model as shown in Table 2 mainly because it does not use any test filtering procedure.

RESULTS AND DISCUSSION

A new experimental set-up has been recently designed at IR-PHE for the same sets of parameters (η = 0.89, Re Ω and Re Q values) but for a relatively small aspect ratio Γ = 50, such that the hydrodynamic flow is never established whatever the values of the axial flow and rotation rates (see extensive comparisons in [START_REF] Poncet | Instabilities, turbulence and heat transfer in confined rotating flows -habilitation thesis[END_REF]). As the in-house code has already been validated for a middle-gap cavity [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF], it will hold for reference data in the following.

From Figures 2 to 5, a strong agreement between the WALE and the Dynamic Smagorinsky models, is obtained whatever the value of the rotation parameter. The profiles are indistinguishable most of the time with slightly higher intensity levels in the core of the flow provided by the WALE model. For the same mesh grid and so the same wall coordinates, the Dynamic Smagorinsky model provides slightly higher values of the friction Reynolds numbers along both walls. Both LES will be used in the following as reference data for comparisons between the RANS and the hybrid RANS/LES approaches. One can recall that the WALE enables to save about 12% of computational cost compared to the Dynamic Smagorinsky model (see Table 2).

The choice of the rotation parameter N = Re Ω /Re Q , as the main parameter governing the hydrodynamics of Taylor-Couette-Poiseuille flows, is first discussed. Indeed, cases B and F have the same value of N = 4.47, though different rotation and axial flow rates. The same normalization as used by Chung and Sung [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF] is introduced: ΩR 1 for the azimuthal mean and root mean square velocities, and W m (or W 2 m ) for other quantities. Figure 2 shows the first and second moments, exhibiting no significant differences between those two cases, whether modeling approach (LES, hybrid RANS/LES or RANS) is considered. It confirms that the rotation parameter pilots the hydrodynamics in such open Taylor-Couette system for this range of N values [1.49 -6.71]. One can expect that rotation effects dominate compared to those of the axial flow as N remains larger than unity. Figs. 3 through 5 compare all four other cases, put two by two on the same graphs, such as the left side is for N = 1.49 (case D) and N = 2.24 (case A), while the right side shows profiles for N = 2.98 (case E) and N = 6.71 (case C). The LES provide rather the same mean velocity distributions whatever the value of N. The LES predict a mean flow divided into three flow regions: two thin boundary layers on each wall separated by a central core at mid-gap in almost solid body rotation. For N = 1.49, the mean tangential velocity slightly decreases with the radius at mid-gap, highlighting weak turbulence intensities. For larger values of N, the slope of ∂ (V θ /ΩR 1 )/∂ r in the free stream, asymptotically tends to 0, such that V θ gets constant, which is typical of turbulent regimes in rotating flow configurations. It may be explained by considering an effective Reynolds number based on the effective velocity V e f f = W 2 m + β (ΩR 1 ) 2 (β a constant generally put to 0.5). At a given axial flow rate, increasing the rotation parameter leads to an increase in the effective velocity or Reynolds number and as a consequence to higher turbulence levels. Furthermore, the axial dimensionless velocity V z is also weakly influenced by N, at least concerning LES and RANS calculations. This latter behavior is different from that of middle-gap (η ≃ 0.5) Taylor-Couette-Poiseuille flows, for comparable values of N. For example, Hadziabdic et al. [START_REF] Hadziabdic | LES of turbulent flow in a concentric annulus with rotating outer wall[END_REF] have shown that increasing N puts the axial velocity towards a laminar Poiseuille-like shape. Here, the radial distribution of V z is typical of what is encountered in turbulent pipe flows regarding the LES and RANS profiles. On the other hand, the hybrid RANS/LES approach exhibits larger radial variations with a profile similar to a laminar Poiseuille-like profile indicating lower turbulence intensities at mid-gap. Regarding the radial distributions of all six Reynolds stress tensor components obtained by LES (Fig. 4 and5), turbulence is mainly concentrated in the thin boundary layers developed along the walls with remarkably high peak values. All profiles in absolute values are slightly asymmetric. The rotor boundary layer appears besides to be more unstable than those of the stator high-lighting the destabilizing effect of rotation for this range of N values. The normal components R θ θ and R zz are twice lower in the core region. The shear stress components R rz and R θ z show strong stresses along both walls. Increasing the value of N leads to an increase of the turbulence intensities in all flow regions. 0 0,2 0,4 0,6 0,8 1 0 The prediction capabilities of both RANS and hybrid RANS/LES approaches used here, are now discussed. Both approaches provide roughly the good profiles with larger discrepancies for some components of the Reynolds stress tensor. Clearly, the hybrid RANS/LES approach does not improve the predictions of the RANS model, which behaves quite well. The EB-RSM model has previously shown some satisfactory results on channel flows with or without global rotation (see e.g. [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF]). But in its hybridization, under the Equivalent DES method, it shows heterogeneous prediction quality, with strong dependence on N. Generally, they both overestimate the peak intensities for the R rr and cross components. Thus, they also overpredict the values of the friction Reynolds number as shown in Table 2.

0,05 0,1 0,15 0,2 (r -R1)/(R2 -R1) R 1/2 rr /Wm 0 0,2 0,4 0,6 0,8 1 0 0,1 0,2 0,3 0,4 (r -R1)/(R2 -R1) R 1/2 rr /Wm 0 0,2 0,4 0,6 0,8 1 0 0,05 0,1 0,15 0,2 (r -R1)/(R2 -R1) R 1/2 θ θ /ΩR1 0 0, 2 0,4 0,6 0,8 1 0 0
,05 0,1 0,15 0,2 (r -R1)/(R2 -R1) R 1/2 θ θ /ΩR1 0 0,2 0,4 0,6 0,8 1 0 0,1 0,2 0,3 0,4 0,5 (r -R1)/(R2 -R1) R 1/2 zz /Wm 0 0,2 0,4 0,6 0,8 1 0 0,2 0,4 0,6 0,8 (r -R1)/(R2 -R1) R 1/2 zz /Wm
For instance, for case D (N = 1.49), the profile of the energy ratio r k = k m /k, shown by Figure 6, betrays an unexpected RANS state, which can even be assimilated to a relaminarization, according to the azimuthal velocity profile (Figure 3). This weakness can be compared to the difficulty of any DES-like ap-

(r -R1)/(R2 -R1) Rrz/W 2 m 0 0, 2 0,4 0,6 0,8 1 -0
(r -R1)/(R2 -R1) R θ z /W 2 m 0 0, 2 0,4 0,6 0,8 1 0 0,1 0,2 (r
-R1)/(R2 -R1) R θ z /W 2 m FIGURE 5
. Shear stress components of the Reynolds stress tensor. See legend of Figure 3.

proach, to sustain a fluctuating state in quite stable configurations. Even case A (N = 2.24) shows a similar tendency. For higher values of N, statistics of the flow seem to be better predicted, however with flat profiles for the highest value (case C, N = 6.71). Indeed, while Figure 6 shows that the energy ratio r k is well driven in the free flow, second moment profiles (see Figs. 5 and4) seem to exhibit which could be considered as unphysical resolved structures. This would be the result of an unsufficient dissipation rate, thus suggesting that the value of β 0 (see Eq.( 2)) should be lowered. On the fringe of the previous observations, one can have a look upon the isovalues of the Q-criterion shown on Figures 7 and8. The isovalues of the Q-criterion obtained by the LES-WALE are shown on the rotor and stator sides, while for the hybrid RANS/LES, the isovalues are shown only in front of the rotor, since there are far less coherent structures on the stationary cylinder. The LES highlights the presence of 3D unsteady coherent structures within the two boundary layers. They appear as thin negative spiral patterns along the rotor as they roll up in the opposite sense of the inner cylinder rotation. The same spiral network is obtained along the stator with a positive angle. When the axial throughflow is more important (weak values of N), these structures are more aligned with the axial direction. They get progressively inclined with the tangential direction when the rotation rate increases. As examples, along the rotor, the angle formed by the spirals with the tangential direction is equal to -45 • for N = 1.49 and to -16 • for N = 6.71. Despite the hybrid approach captures less structures, these spiral rolls form rather the same angle with the tangential direction as in the LES, except for case C, where the streaks captured by the hybrid RANS/LES are almost horizontal. These 3D unsteady structures appearing as thin spiral rolls along both walls at a relatively small scale in the LES are responsible for the large peak values observed in the Reynolds stress tensor profiles. They may play a key role in the wall heat transfer process. They could explain why most of the RANS models, which assume the base flow as being stationary and axisymmetric, fail to predict the right distributions of the heat transfer coefficient.

Finally, Table 2 contained in Appendix A, shows some numerical parameters, as well as wall quantities, for all six cases, and for all four approaches (both LES, hybrid RANS/LES and RANS). Globally, this table confirms that the RANS approach provides satisfactory results for these quantities. Another observation that can be made, is that the friction Reynolds numbers are poorly predicted by the hybrid approach on both walls. This may be linked to the fact that the RANS mode (r k = 1) is not well enforced on either wall (Figure 6).

CONCLUDING REMARKS

The present work is a further step for Friess et al. [START_REF] Friess | Taylor-Couette-Poiseuille flows : from RANS to LES[END_REF], on six different values of the rotation parameter N providing a wider numerical data set of Taylor-Couette-Poiseuille flows in a narrow-gap rotor-stator cavity (η = 0.89). Since the LES code was validated in a middle-gap (η = 0.5) configuration [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF], and since there is no available experimental data set for a fully established narrow-gap flow (see e.g. [START_REF] Poncet | Instabilities, turbulence and heat transfer in confined rotating flows -habilitation thesis[END_REF]), its results may be con- sidered as reference. One used here two subgrid scale models, the Dynamic Smagorinsky and WALE models, providing very similar results for both the mean and turbulent fields, while the WALE enables to save about 12% of computational effort. The N values are in the range [1.49; 6.71], such that a parametric study can be performed. Three observations may be put forth: (i) as for middle-gap Taylor-Couette-Poiseuille flows, N is the main parameter governing the flow, (ii) the dimensionless azimuthal velocity slightly varies with N, (iii) unlike middle-gap Taylor-Couette-Poiseuille flows, the axial velocity profile does not significantly vary with N, even for high N values.

The conclusions about the performance of the so-called "Equivalent DES" hybrid RANS/LES approach, are as follows: (i) no significant improvement in predicting flow statistics, by comparison with the underlying RANS (EB-RSM) model, but (ii) Q-criterion isocontours exhibit net-shaped streaks, as well as the changing of the angle of the main directions of these structures, in relatively good accordance with the reference LES.

Eventually, one can say that this work is a new illustration of the fact that seamless hybrid RANS/LES methods are better suited for flows far from equilibrium, for instance involving separation. However, further theoretical work on "grey zones" treatment, could shake things up.

ACKNOWLEDGMENT

This work was financially supported by Liebherr Aerospace (project "Entrefer moteur") and by the LABEX MEC (ANR-11-LABX-0092) through the HYDREX program, which are gratefully acknowledged. Furthermore, this work was granted access to the HPC resources of IDRIS under the allocation 2013-090242 made by GENCI (Grand Equipement National de Calcul Intensif).

Appendix A: Wall resolution and mean flow parameters

2 mFIGURE 2 .

 22 FIGURE 2. Azimuthal and axial velocity components and all six Reynolds stress tensor components for N = 4.47. Cases F (shifted up) and B. : Dynamic Smagorinsky model, • : WALE model, ---: hybrid RANS/LES, --: RANS.

FIGURE 3 .

 3 FIGURE 3. Azimuthal (top) and axial (bottom) mean velocities. Left : cases A (shifted up) and D. Right : cases C (shifted up) and E. Same legend as Figure 2.

FIGURE 4 .

 4 FIGURE 4. Normal components of the Reynolds stress tensor. See legend of Figure 3.

5 (FIGURE 6 .

 56 FIGURE 6. Energy ratio r = k m /k tot . Left : cases F (shifted up) and B. Right : cases A (shifted up) and D. Bottom : cases C (shifted up) and E. ---: targetted, -observed.

FIGURE 7 .

 7 FIGURE 7. Isocontours of the Q-criterion in the full computational domain, colored by the radial coordinate r (blue on the rotor and red on the stator). Left : Case B, right : case F. From top to bottom : LES-WALE, respectively rotor and stator side, and hybrid RANS/LES, rotor side. The inner cylinder is rotating counter-clockwise and the axial throughflow is from the bottom to the top.

FIGURE 8 .

 8 FIGURE 8. Isocontours of the Q-criterion, colored by the radial coordinate r. Left : Case A, right : case C. Same legend as Figure 7.

  ).

	Case	A	B	C	D	E	F
	Γ			10
	η			8/9
	Re Q		3745			5617
	Re Ω 8378 16755 25133 8378 16755 25133
	N	2.24	4.47	6.71	1.49	2.98	4.47

TABLE 1 .

 1 Geometrical and flow parameters

FIGURE 1. Sketch of the narrow-gap cavity

TABLE 2 .

 2 Computational details and mean flow parameters.

Copyright c 2014 by ASME


G Fontaine

S Poncet

E Serre

Multidomain extension of a pseudospectral algorithm for the direct simulation of wall-confined rotating flows

In this work, we improve an existing pseudospectral algorithm, in order to extend its properties to a multidomain patching of a rotating cavity. Viscous rotating flows have been widely studied over the last decades, either on industrial or academic approaches. Nevertheless, the range of Reynolds numbers reached in industrial devices implies very high resolutions of the spatial problem, which are clearly unreachable using a monodomain approach. Hence, we worked on the multidomain extension of the existing divergence-free Navier-Stokes solver with a Schur approach. The particularity of such an approach is that it does not require any subdomain superposition: the value of a variable on the boundary between two adjacent subdomains is treated as a boundary condition of a local Helmholtz solver. This value is computed on a direct way via a so-called continuity influence matrix and the derivative jump of an homogeneous solution computed independently on each subdomain. Such a method is known to have both good scalability and accuracy. It has been validated on two well documented three-dimensional rotating flows.

Numerical modelling

Let's introduce the numerical fundamentals of the present method. A pseudospectral method is used to solve the Navier-Stokes PDE system in an annular cavity, where incompressibility is assured through a projection method.

Pseudospectral methods

Let Ω be the inner points of an annular cavity and Γ the domain boundary. The spatial approximation is of Chebyshev type in the axial z and radial r directions, of Fourier-Galerkin type in the azimuthal θ direction. Let Ψ be a variable (u, v, w, p, φ), which can be written for any point (r, θ , z) ∈ Ω ∪ Γ as:

Ψ N r N θ N z (r, θ , z) = N θ /2 ∑ k=-N θ /2 N r -1 ∑ n=0 N z -1 ∑ m=0 Ψ nkm T n (r) T m (z) e ikθ ( 1
)
where T n and T m are the Chebyshev polynomials of degrees n and m respectively. N r , N z and N θ are the approximation degrees in the radial, axial and azimuthal directions, respectively. Ψ nkm is given by :

Ψ nkm = 1 K 1 c k c ′ m K-1 ∑ q=0 N ∑ i=0 M ∑ j=0 1 c i c ′ j Ψ (r i , θ q , z j ) cos ( inπ N
) cos

( imπ M ) e -ikθ q
(2)

c 0 = c ′ 0 = c N = c ′ M = 2 and c k = c
′ m = 1 and for n = 1, ..., N -1 and m = 1, ..., M -1. This approximation is done using a Gauss-Lobatto point distribution in the radial and axial directions with Fourier-Galerkin points in the azimuthal direction. Thus derivatives can be estimated in the spectral space with a good precision, allowing the use of efficient FFT algorithms. These methods are called pseudospectral because of the non-linear diffusive terms, which must be computed in the physical space.

Geometry

The velocities are made dimensionless according to the Reynolds number Re = ωR 2 1 /ν, where ω is the angular velocity, R 1 the outer radius of the domain, ν the kinematic viscosity. Hence the dimensionless geometry is defined by two parameters: the aspect ratio of the cavity L = R 1 -R 0 2h and the curvature parameter

Rm = R 1 +R 0 R 1 -R 0 ,
where R 1 and R 0 are the outer and inner radii of the cavity, respectively, and 2h its height [START_REF] Serre | Instabilité de couche-limite dans les écoulements confinés en rotation. Simulation numérique directe par une méthode spectrale de comportements complexes[END_REF].

Projection method

Let's consider the Navier-Stokes equations in primitive variables :

( ∂ V ∂t + (V.∇)V ) = - 1 ρ ∇p + ν∆ V + F in Ω (3)
∇.V = 0 in Ω ( 4
)
where V is the velocity vector with (u, v, w) its components in the cylindrical basis, p the pressure, ρ the density and F a force term.

Time discretization

A semi-implicit second order scheme has been chosen for the time discretization with an implicit retarded Euler scheme of 2nd order for the diffusive terms and an explicit Adams-Bashforth evaluation for the non-linar convective terms.

Projection scheme

For incompressible viscous flows, an equation is missing to describe the pressure evolution. We deal with this particularity by using an improved projection method [START_REF] Raspo | A spectral projection method for the simulation of complex three-dimensional rotating flows[END_REF] based on the Goda's projection method. It requires 3 steps:

• First step: the computation of a preliminary pressure pn+1 through this Poisson equation:

           ∇ 2 pn+1 = ∇.[-2V.∇V n + V.∇V n-1 + F n+1 ] in Ω ∂ pn+1 ∂ n = n.[ -3W n+1 + 4V n -V n-1 2δt -2V.∇V n + V.∇V n-1 + ν(2∆ V n -∆ V n-1 ) + F n+1 ] on Γ (5)
where W represents the boundary conditions. • Second step: prediction step. Using the gradient of pn+1 , a preliminary velocity field V * is computed through 3 Helmholtz solvers:

       3V * -4V n + V n-1 2δt + 2V.∇V n -V.∇V n-1 = -∇ pn+1 + ν∆ V * + F n+1 in Ω V * = W n+1 on Γ (6)
This velocity field does not a priori satisfy the incompressibility constraint in Ω . The principle of the projection method is namely to project this field on a divergence-free field.

• Third step: pseudo-pressure calculation and correction. An intermediate variable called pseudopressure φ is computed through a Poisson solver:

φ = 2δt 3 (p n+1 -pn+1 )
As ∇.V n+1 = 0, the Poisson problem to solve for φ is :

     ∇ 2 φ = ∇.V * ∂ φ ∂ n = 0 (7)
Corrected pressure and velocity fields may be then evaluated at time step n + 1:

     p n+1 = p n+1 + 3 2δt φ V n+1 = V * -∇φ (8)

Resolution algorithm : complete matricial diagonalisation technique

For both Poisson and Helmholtz solvers, complete matricial diagonalisation technique is used. For the Poisson solver, this technique exhibits a null-eigenvalue problem, which is treated by a "source term reset" technique ( [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF]), for indexes in the operators corresponding to the null-eigenvalue index. According to the variable spatial evaluation, each Helmholtz/Poisson solving operation in the physical tridimensional domain is equivalent to N θ bidimensional Helmholtz/Poisson solving operations in the (r, z) plane, each k ∈ [[1; N θ ]] being the azimuthal mode in the flow spectrum [START_REF] Serre | Instabilité de couche-limite dans les écoulements confinés en rotation. Simulation numérique directe par une méthode spectrale de comportements complexes[END_REF].

For each mode k and for each variable Ψ = (u, v, w, p, φ), the following bidimensional system has then to be solved :

{ ∆ Ψ -σ k Ψ = S k in [-1, 1] × [-1, 1] AΨ = b on Γ (9)
Ψ being the Fourier transform of Ψ in the azimuthal direction. The azimuthal properties of the differentiation matrixes are treated in σ so as to allow the ∆ operator to be independant of the azimuthal mode. The azimuthal dependance of what will follow is only linked to σ 's one. Hence we will limit the multidomain approach to bidimensional problems, because the extension to tridimensional flows is quite immediate, thanks to spectral properties in the Fourier-Galerkin direction.

Multidomain approach: influence matrix technique

This technique is a direct Schur multidomain technique used by Raspo [START_REF] Raspo | A direct spectral domain decomposition method for the computation of rotating flows in a T-shape geometry[END_REF] for rotating flows using the vorticity-stream function formulation. It requires the patching of the cavity into subdomains without subdomain covering. The values of each variable on the frontier between two subdomains are treated as boundary conditions in the Helmholtz local solvers, the particularity of the influence matrix technique being the determination of this condition through a direct matrix computation.

Multidomain geometry

We will limit the discussion to a radial multidomain decomposition (Fig. 1), because the curvature terms vary only along this direction. The generalization of this technique to an axial decomposition is quite immediate. We introduce the local aspect ratio of the subdomain m denoted L (m) and its local curvature parameter Rm (m) for m ∈ [[1; M]] (M the number of subdomains), which satisfy : Local derivation matrixes are deduced from these, in order to have a good approximation for the curvature terms from one subdomain to another and to adapt the local mapping to the one which would be used in a monodomain approach.

M ∑ m=1 L (m) = L ( 10
)
Rm (m) -1 = Rm (m+1) + 1 (11)

Multidomain decomposition of the solutions

Let Ψ be either (u, v, w, p, φ), Ω (m) and Ω (n) two adjacent subdomains, ξ the border between these two subdomains and λ the value of Ψ on ξ . For both subdomains, the local problems to be solved may be written as :

{ ∆ (m) Ψ (m) -σ (m) Ψ (m) = S (m) in Ω (m) A (m) Ψ (m) = b (m) on Γ (m) (12) { ∆ (n) Ψ (n) -σ (n) Ψ (n) = S (n) in Ω (n) A (n) Ψ (n) = b (n) on Γ (n) (13)
The resulting problem is that the boundary conditions to be imposed on the parts of Γ (m) and Γ (n) corresponding to ξ are unknown. To find λ , we choose to ensure both C 0 and C 1 continuities through ξ . Ψ is written as the combination of an homogeneous solution Ψ and a boundary solution Ψ : Ψ = Ψ +Ψ .

On Ω (m) :

     ∆ (m) Ψ (m) -σ (m) Ψ (m) = S (m) in Ω (m) A (m) Ψ (m) = b (m) on Γ (m) Ψ (m) = 0 on ξ (m) (14)      ∆ (m) Ψ (m) -σ (m) Ψ (m) = 0 in Ω (m) A (m) Ψ (m) = 0 on Γ (m) Ψ (m) = λ (m) on ξ (m) (15)
On

Ω (n) :      ∆ (m+1) Ψ (n) -σ (n) Ψ (n) = S (n) in Ω (n) A (m) Ψ (n) = b (n) on Γ (m) Ψ (n) = 0 on ξ (16)      ∆ (n) Ψ (n) -σ (n) Ψ (n) = 0 in Ω (m) A (n) Ψ (n) = 0 on Γ (n) Ψ (n) = λ on ξ ( 17
)
If λ is assumed to be known, one can verify easily that (eq. 12)=(eq. 14)+(eq. 15) and (eq. 13)=(eq. 16)+(eq. 17).

The influence matrix technique

Let's consider the boundary solution Ψ . It can be written as the linear combination of Green's elementary solutions G (m) kξ , defined for each subdomain Ω (m) by :

         ∆ (m) G (m) kξ -σ (m) G (m) kξ = 0in Ω (m) A (m) G (m) kξ = 0on Γ (m) G (m) kξ (η l ∈ ξ (m) ) = δ kl ∀l ∈ [[1; N ξ ]] (18)
Assuming that ξ is the only boundary (i.e. there are only 2 subdomains Ω (m) and Ω (n) ), the boundary solution should be written as:

           Ψ (m) = Ψ (m) + N ξ ∑ k=1 λ k G (m) kE in Ω (m) Ψ (n) = Ψ (n) + N ξ ∑ k=1 λ k G (n) kS in Ω (n) (19) Ψ (m) ∪ Ψ (n) is obviously continuous through ξ , but not ∂ Ψ (m) ∂ r ∪ ∂ Ψ (n) ∂ r .
As the boundary solution is continuous too, Ψ (m) ∪Ψ (n) should be continuous for any value of λ . The influence matrix technique aims to find λ in order to make it C 0 and C 1 through ξ . We denote ∂ r = ∂ ∂ r . This C 1 -continuity problem on ξ writes:

∂ r Ψ (1) (ξ ) -∂ r Ψ (2) (ξ ) = N ξ ∑ l=1 λ l [∂ r G (2) lS (ξ ) -∂ r G (1) lE (ξ )] (20)
This can be written in a matrix form:

D = M λ (21)
where D k is the derivative (time-dependent) jump vector and M the continuity influence matrix of the problem. Note that it depends only on the time-independent Green solutions, so it just has to be computed in pre-processing. This matrix is diagonal-dominant. If there are more than 2 frontiers in the domain, the influence matrix is built by blocks. The block dimension is then N f ront , the number of frontiers. Each diagonal block is a Green derivative jump vector along each frontier. Some non-diagonal blocks appear, resulting locally of the influence of 2 frontiers on one another through a single subdomain, as shown on Figure [START_REF] Raspo | A spectral projection method for the simulation of complex three-dimensional rotating flows[END_REF].

If M is inversible, we can find λ as:

λ = M -1 D (22)
This is achieved using LAPACK subroutines. This vector is then used as a boundary condition on ξ in the local Helmholtz solvers to get a C 0 and C 1 Ψ (m) ∪Ψ (n) solution.

Singularity of the Poisson-problem

The Neuman-Poisson problem has an infinity of solutions defined up to an additive constant. As Dirichlet boundary conditions are implemented on the frontiers, this problem no longer exists locally. Nevertheless, it is transposed to the influence matrix of the Poisson-problem of the k = 0 Fourier mode, which has a null-eigenvalue.

It is treated by a diagonalisation technique of this modal matrix. The derivative jump is expressed in the diagonalisation basis and its i 0 -th component is set to zero, if i 0 is the null eigenvalue index, as proposed by Abide and Viazzo [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF].

Figures 3 and4 show the decrease of the quadratic truncature error L 2 when one increases N r and M respectively. In the first case, the spectral convergence is obtained. As we impose only a C 1 continuity through the frontier, the spatial accuracy increases indeed faster with N r than with M, because in its last case, it multiplies the number of interfaces and so the number of C d discontinuities (d ≥ 2). We consider two multidomain rotating flow configurations: a high aspect ratio Taylor-Couette system and an interdisk rotor-stator flow.

Axial decomposition : Taylor-Couette flow

The numerical parameters are fixed to M = 7, (N r , N θ , N z ) = (35, 24, 85) (in each subdomain) with a time step set to δt = 10 -3 . We consider a Taylor-Couette cavity with rotating terminal disks characterized by L = 0.025 and Rm = 12.33. We introduce the parameter ε = Ta/Ta c , where Ta c is the critical Taylor number at which "Taylor rolls" appear. For ε = 1.36, the flow becomes tridimensional (Wavy Vortex Flow, Fig. 6) with a dominant Fourier mode k = 3, which is in perfect agreement with the previous DNS results of Serre et al. [START_REF] Serre | Stability of Taylor-Couette flow with radial throughflow[END_REF].

Radial decomposition: interdisk rotor-stator flow

We consider the rotor-stator interdisk cavity considered by Poncet et al. [START_REF] Poncet | Revisiting the two first instabilities of the flow in an annular rotor-stator cavity[END_REF] for which L = 6.26 and Rm = 1.8. The domain is decomposed into 4 subdomains with (N r , N θ , N z ) = (45, 48, 45) in each subdomain and a time step equal to δt = 10 -4 . The transition to tridimensional flow appears at Re = 25000 with the appearance of 17 spiral arms, whose characteristics fully agree with the ones obtained by [START_REF] Poncet | Revisiting the two first instabilities of the flow in an annular rotor-stator cavity[END_REF].

Fig. 1

 1 Fig. 1 Example of a multidomain decomposition in the radial direction with 3 subdomains

Fig. 2

 2 Fig. 2 Block definition of the influence matrix for an 1-D multidomain decomposition

Fig. 3 Fig. 4

 34 Fig. 3 Spatial accuracy evolution with the total number of points (N r N z N θ M) when N r varies : spectral convergence

Fig. 5

 5 Fig. 5 Iso-value v = 0.5 for a Wavy Vortex Flow at ε = 1.36

Fig. 6

 6 Fig. 6 Iso-value w = 0.005 for a 3D rotor-stator flow at Re = 25000: 17 spiral arm structures


C Friess

S Poncet

S Viazzo

BENCHMARK OF TURBULENCE MODELINGS IN A TAYLOR-COUETTE-POISEUILLE SYSTEM

Introduction

The present investigation is concerned with the numerical modeling of the turbulent flow in a narrow gap Taylor-Couette-Poiseuille (TCP) system. The flow is confined between two coaxial cylinders of length h, with an inner rotating cylinder (radius R 1 ) and a outer stationary one (radius R 2 ). An axial thoughflow of fluid is supplied within the gap. Such turbulent flows have many applications in process engineering (dynamic membrane filtration) and also in the turbomachinery industry (rotating heat exchangers).

The effect of an axial throughflow in a Taylor-Couette system has been considered experimentally by [START_REF] Kaye | Modes of Adiabatic and Diabatic Fluid Flow in an Annulus With an Inner Rotating Cylinder[END_REF], who showed the existence of four flow regimes depending on the incoming flow velocity W m and on the rotation rate Ω: laminar and turbulent flows, with or without Taylor vortices. Most of the experimental works so far have been performed in systems characterized by a high aspect ratio Γ = h/∆R > 50 (∆R = R 2 -R 1 ) and a moderate radius ratio 0.33 < η = R 1 /R 2 < 0.66. Thus, [START_REF] Nouri | Flow of Newtonian and Non-Newtonian Fluids in a Concentric Annulus With Rotation of the Inner Cylinder[END_REF] then [START_REF] Escudier | Concentric annular flow with centerbody rotation of a Newtonian and a shear-thinning liquid[END_REF] provided very useful velocity measurements for (Γ = 98, η = 0.496) and (Γ = 244, η = 0.506) respectively. [START_REF] Naser | Prediction of Newtonian and Non-Newtonian flow through concentric annulus with centerbody rotation[END_REF] compared the predictions of a kε model with the experimental data of [START_REF] Escudier | Concentric annular flow with centerbody rotation of a Newtonian and a shear-thinning liquid[END_REF]. The model showed large discrepancies for the mean velocity components, with a strong dependence on the axial position not observed in the experiments. [START_REF] Kuosa | Numerical and experimental modelling of gas flow and heat transfer in the air gap of an electric machine[END_REF] performed a numerical benchmark in a Taylor-Couette system with a radial inlet. None of the two-equation models offer satisfactory results even for the mean tangential fluid velocity. [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF] performed Large Eddy Simulation (LES) using second-order spatial schemes. They compared their numerical data to the experimental ones of [START_REF] Nouri | Flow of Newtonian and Non-Newtonian Fluids in a Concentric Annulus With Rotation of the Inner Cylinder[END_REF] for η = 0.5 and three values of the rotation parameter N = ΩR 1 /W m = [0.2145; 0.429; 0.858]. Recently, [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF] validated their code based on fourth-order compact schemes in the inhomogeneous directions on the same reference data. Their results slightly improve those of [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF] highlighting the importance of high-order schemes.

One will focus our attention here on the narrow gap case, which models quite faithfully flows in the rotorstator gap of electrical motors. The aim of the present work is twofold : (i) providing some reference LES data and (ii) questioning the capabilities of a hybrid RANS / LES method, as well as the underlying RANS model, in predicting such flow.

Flow configuration and parameters

The cavity shown in Figure 1 may be characterized by two geometrical parameters : its aspect ratio Γ = h/∆R = 10 and its radius ratio η = 8/9. The inner cylinder is rotating at a constant rate Ω, while the outer one is stationary. An axial throughflow is imposed within the gap at a constant bulk velocity W m . The main flow parameters are the rotational Reynolds number Re Ω = ΩR 1 (R 2 -R 1 )/ν (in the range [8.4 × 10 3 -2.5 × 10 4 ]) and the bulk Reynolds number Re Q = W m (R 2 -R 1 )/ν (= 3745 or 5617), ν being the fluid kinematic viscosity. The results will be discussed in terms of the rotation parameter N = Re Ω /Re Q , which will take four different values: N = 1.49, 2.98, 4.47 and 6.71.

Numerical methods

The in-house code is based on semi-implicit and second-order accurate temporal schemes. The derivatives are approximated using fourth-order compact formula in the radial and axial directions and a Fourier approximation in the tangential direction. The time splitting scheme is an improved projection method ensuring the incompressibility at each time step. The multidomain solver ensures the continuity of the solution and its first normal derivative across the conforming interface using an influence matrix technique [START_REF] Abide | A 2D compact fourthorder projection decomposition method[END_REF]. Periodic boundary conditions are applied in the axial and circumferential (3π/4 periodicity) directions and no-slip boundary conditions are imposed on the walls. The cavity has been decomposed into two subdomains in the axial direction. Two subgrid scale models will be tested in the following: a dynamic Smagorinsky model (denoted LES-SD) and the WALE model (LES-WALE). Two mesh grids containing (N r ×N θ ×N z )=(65×144×130) and (65 × 144 × 130) points have been used for N = 6.71 and N ≤ 4.47 respectively. The time step δt remains in the range [3-8]×10 -5 s. The present code has been fully validated by [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF] for η = 0.5.

The code used for RANS and hybrid RANS/LES simulations is the open-source Code Saturne, developed by EDF [START_REF] Archambeau | Code Saturne: a finite volume code for the computation of turbulent incompressible flows -Industrial applications[END_REF]. It is a finite volume solver written in cartesian coordinates. A SIMPLEC algorithm, with the Rhie and Chow interpolation, is used for pressure-velocity coupling. Concerning hybrid RANS/LES calculations, convective fluxes are approximated by a second-order centered scheme for momentum and a first-order upwind scheme for subfilter quantities. Time marching uses a Crank-Nicholson second-order scheme. In RANS mode, convective fluxes are approximated by a firstorder upwind scheme and a first-order time scheme. For both approaches, the elliptic blending Reynolds Stress Model of Manceau and Hanjalic (2005) is used. The hybrid method is the equivalent DES of Manceau et al. (2010). The hybrid RANS/LES mesh contains 60 3 cells assuming a π/2-periodicity in the tangential direction and a periodic flow in the axial direction. The time step (normalized by t P = ∆R/W m ) is fixed to 1.4 × 10 -2 . The RANS mesh is a 1D grid of 60 cells. No-slip boundary conditions are imposed on the walls.

Turbulence models

Various levels of modelling are presented. Emphasis is first put on the hybrid RANS/LES method. Afterwards, the subfilter closure and RANS model will be briefly presented, as well as the subgrid-scale models used for LES.

Hybrid method

The hybrid method used for the present work is the "Equivalent DES" of Manceau et al.(2010). This approach was first derived for the purpose of bridging the PITM (Partially Integrated Transport Model) of [START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF], and the DES (Detached Eddy Simulation) method. Indeed, the first is fully justified from a theoretical point of view, while the second one was developped on a rather phenomenological ba-sis. First, the PITM was generalized by Fadai-Ghotbi et al.to inhomogeneous flows, considering temporal filtering, rather than spatial filtering. Indeed, inhomogeneous flows are more frequent. Since, in this kind of flow, the RANS operator corresponds, by ergodicity, to temporal averaging. Since any seamless (continuously transitioning) RANS/LES method must tend to RANS, at large filter width, inhomogeneous flow studies need then to consider hybrid RANS / temporal LES. Nevertheless, from a pragmatical point of view, considering RANS / TLES hybridization does not cause any difficulty in implementing models or special terms, since Fadai-Ghotbi et al.showed that applying Temporal PITM (T-PITM) to a inhomogeneous stationary flow is just equivalent to applying PITM to homogeneous, statistically unsteady flow. Actually, one must just keep in mind that "hybrid RANS/LES" is to be understood in a general way, including temporal LES. The advantage of Detached Eddy Simulation lies in its simplicity and robustness. The idea of Manceau et al.(2010), was thus to derive an approach bridging it with T-PITM. This latter is, in spite of its theoretical justification, not very easy to implement in any code, showing some numerical issues. An equivalence criterion was then determined between T-PITM and DES, providing some theoretical justification to the latter, and allowing to interprete it as a hybrid RANS / temporal LES method. This equivalence criterion was derived analytically for equilibrium flows, but was successfully tested on a flow over a periodic hill, involving massive separation [START_REF] Friess | Investigation of the equivalence of two hybrid temporal-LES methods based on elliptic blending in separated flows[END_REF]. The principle of DES is to magnify the dissipation term of the transport equation for either subfilter turbulent kinetic k SF S energy or subfilter stresses τ ij SF S , damping the modelled energy and allowing large scale eddies to be resolved:

ε DES ij = k 3/2 SF S ε SF S L ε ij SF S . (1)
While classical DES uses the local grid step to determine the length scale L in (1), equivalent DES uses the ratio r defined as :

r = k m k , ( 2
)
where k is the total turbulent kinetic energy, and k m its counterpart contained in the modelled (subfilter) scales. At the RANS limit, r tends to 1, and at the DNS limit, it tends to 0. The advantage of r is that it can be estimated using a spatial or a temporal energy spectrum. As shown by [START_REF] Friess | Investigation of the equivalence of two hybrid temporal-LES methods based on elliptic blending in separated flows[END_REF], r can be evaluated as :

r = 1 k ∞ ωc E T (ω) dω = min 1, 1 β 0 U s √ k 2/3 ω c k ε -2/3 . (3)
β 0 being a constant derived from the Kolmogorov constant (we used 0.3 here), U S a sweeping velocity and ω c the cutoff frequency of the considered filter. It can be defined as :

ω c = min π dt ; U s π ∆ , (4)
where dt and ∆ are the time and grid steps, respectively. Finally, following Manceau et al.(2010), the length scale L entering (1) is :

L = r 3/2 1 + C ε2 -C ε1 C ε1 (1 -r Cε1/Cε2 ) k 3/2 ε , (5)

RANS model and subfilter closure

In the temporally-filtered Navier-Stokes equation, the subfilter stresses must be modeled such a way that, making the temporal filter width go to infinity, the equations tend to the RANS equations. In particular, one of the main objectives of hybrid methods is to use RANS closures in the near-wall regions, to avoid the very fine resolution required by LES. In that aim, the RANS model proposed by Manceau & Hanjalić (2002), the so-called elliptic blending Reynolds-Stress Model (EB-RSM), is adapted to the hybrid temporal LES context. In this model, an elliptic relaxation equation is solved for a scalar α :

α -L 2 SF S ∇ 2 α = 1 , (6)
which is a sensor of the distance to the wall (α = 0 at the wall, and 1 far away), and is used to blend nearwall and homogeneous formulations for the redistribution and dissipation terms of the transport equations for subfilter stresses (see Fadai-Ghotbi et al., 2010, for details) :

φ * ij -ε ij = α 3 (φ * ij -ε ij ) wall +(1-α 3 )(φ * ij -ε ij ) homogeneous .
(7) Under its RANS form, EB-RSM was successfully applied to several flows, such as non-rotating and rotating channels, impinging jets or mixed and natural convection flows. It was recently adapted by Fadai-Ghotbi et al.. (2010) to serve as a subfilter-stress model in the framework of T-PITM. In the present work, it is also applied as a model for the equivalent DES, by simply substituting the dissipation in the fashion of Eq. ( 1) in the Reynolds-stress transport equations.

Subgrid-scale model

Eventually, LES calculations were performed using the Wall-Adapting Local Eddy Viscosity (WALE) model of [START_REF] Nicoud | Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor[END_REF]. This choice was motivated by a better numerical stability than the dynamic Smagorinsky model, and its quality in near-wall treatment without damping functions. The subgrid viscosity is given by :

ν t = (C w ∆) 2 (S d ij S d ij ) 3/2 (S ij S ij ) 5/2 + (S d ij S d ij ) 5/4 , (8)
where C w is a constant, ∆ the grid step, S ij the filtered strain-rate tensor, and S d ij the deviatoric part of the squared filtered velocity gradient tensor. Figure 2 presents the iso-values of the Q criterion obtained by LES-WALE for the two extrema values of the rotation parameter N . It highlights the presence of 3D unsteady coherent structures within the two boundary layers. They appear as thin negative spiral patterns along the rotor as they roll up in the opposite sense of the inner cylinder rotation. They are very similar to those obtained in a middle gap cavity by [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF]. The same spiral network is obtained along the stator with a positive angle. For N = 1.49, these structures are more aligned with the axial direction as the effect of the axial flow is comparable to the one of rotation (N ≃ 1). They get progressively inclined with the tangential direction when the rotation rate (or N ) increases. As examples, along the rotor, the angle formed by the spirals with the tangential direction is equal to 45 • for N = 1.49 and to 16 • for N = 6.71.

Results

a numerical point of view, there is no evidence of structures at the interface between the two subdomains produced by a numerical artifact, which valids the multidomain approach. These structures may play a key role in the wall heat transfer process. They could explain why most of the RANS models, which assume the base flow as being stationary and axisymmetric, fail to predict the right distributions of the heat transfer coefficient.

Figure 3 presents the iso-values of the Q criterion, obtained by the hybrid RANS/LES approach, for the two highest values of the rotation parameter N . Due to the use of a coarser mesh, only the rotor side is N = 4.47 N = 6.71 shown, contrary to LES results. Despite the structure density is lower than for LES, net-shaped streaks are captured, with angles very similar to the ones given by LES. It is noticeable that some streaks are also obtained within the boundary layers, but this approach is unable to capture so well the thin coherent structures, unlike LES. The reason why N = 1.49 is not shown for hybrid RANS/LES, is that the calculation ended up to a RANS level, thus without capturing a ny coherent structures.

Figures 4 through 9 show a comparison between LES, RANS and hybrid RANS/LES statistics of the flow, for each value of the rotation parameter N . Hereafter, the LES results will be commented first, then the RANS results. Eventually, the hybrid RANS/LES will be paid more attention to.

The two LES provide the same mean velocity profiles (Fig. 4). The azimuthal velocity profiles, normalized by the rotor's speed, V * θ = V θ /(ΩR 1 ) vary very weakly with the rotation parameter N such that V * θ = 0.5 at mid-gap. One can just notice that, at low N values, V * θ varies with 1/r in the gap and gets constant at high N . In that case, the profile resembles the Batchelor profile obtained in unmerged boundary layer interdisk flows. The axial velocity profiles is close to the turbulent Poiseuille profiles encountered in pipe flows and do not vary with N . Regarding the distributions of the two Reynolds stress tensor components, turbulence is mainly concentrated along the walls where remarkably high turbulence intensity peaks are observed. The LES-WALE may be considered as the reference data as it has been fully validated by [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF] for η = 0.5. The LES-SD provides here very similar results. It slightly overestimates the R θθ component within the gap at N = 6.71. The LES-WALE is moreover less time consuming as the time per iteration is 16.09 s for N = 6.71 to be compared to 18.53 s for the LES-SD, with similar wall coordinates.

The RANS model predicts reasonably the statistics of the flow, though the normalized tangential velocity profiles V * θ seem quite insensitive to the value of N . However, the near-wall peaks predicted by the reference LES, for most of Reynolds stresses, are rather well captured by the EB-RSM model, especially their radial location. The asymmetry of the aforementionned peaks, is not well reproduced by EB-RSM, but this may be due to numerics, rather than the model itself. A further study is needed, on that point. Eventuially, those quite satisfactory results may be linked to the fact that the EB-RSM model is supposed to be well suited for wall-bounded flows involving attached boundary layers and exhibiting three-dimensional features.

The hybrid RANS/LES approach's success is mixed, with a significant dependance, not only on N , but also on the effective Reynolds number Re eff . This latter is based on the effective velocity V eff = W 2 m + 0.5 (ΩR 1 ) 2 , which is a rough estimate of the resulting bulk velocity, in the plane tangential to the walls.

Conclusions

The present work is a step forward in the modeling of turbulent Taylor-Couette-Poiseuille flows in a narrow gap cavity. Four values of the rotation parameter N have been considered in the framework of a benchmark including two LES implemented in an in-house code and RANS and hybrid RANS/LES approach implemented in Code Saturne. Assuming that Large Eddy Simulation results can be considered as reference data [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF], the hybrid RANS/LES approach does not really improve the predictions of the RANS model for both the mean and turbulent fields, especially for the cases with lower values of the Reynolds number Re eff based on the effective velocity. This tendency has already been observed by Fadai-Ghotbi et al. (2010) for turbulent channel flows at moderate Reynolds numbers. The LES models highlight the presence of 3D coherent structures along both walls, appearing as two spiral networks whose angles strongly vary with N . In spite of the contrasted results of the hybrid RANS/LES approach on the statistics of the flow, instantaneous fields show similar features as LES, especially net-shaped streaks with acceptable angles.

Figure 1 :

 1 Figure 1: Sketch of the TCP system.

Figure 2 :

 2 Figure 2: Iso-values of the Q criterion along the rotor (a,c) and along the stator (b,d) obtained by LES-WALE for: N = 1.49 (a,b); N = 6.71 (c,d).

Figure 3 :

 3 Figure 3: Iso-values of the Q criterion, colored by the radial coordinate r, obtained by hybrid RANS/LES, for N = 4.47 and 6.71.

Figure 4 :

 4 Figure 4: Radial distribution of the mean tangential V* θ = V θ /(ΩR 1 )velocity component for all four cases ; Comparisons between the LES-WALE (squares), the LES-SD (circles), RANS (dashed lines) and hybrid RANS/LES (solid lines).

Figure 5 :Figure 6 :

 56 Figure 5: Radial distribution of the axial V * z = V z /W m velocity component. Same legend as Figure 4.

Figure 7 :Figure 8 :

 78 Figure 7: Radial distribution of the root-mean square Reynolds stress tensor component R 1/2 * θθ =

Figure 9 :

 9 Figure 9: Radial distribution of the Reynolds stress tensor component R * rz = v ′ r v ′ z /W 2 m . Same legend as Figure 4.

Table 1

 1 are better described by the hybrid approach, than by RANS, for cases N = 2.98 and N = 4.47. Similarly, there is a real gain in the prediction of V seems rather sensitive to the effective Reynolds number Re eff . It is worth mentioning that the overstimation of Reynolds stress components is certainly due to the value of the constant β 0 entering Eq (3). Eventually, one can guess there are thresholds for those Reynolds numbers. Below those thresholds, the hybrid RANS/LES approach does not improve the prediction of the flow statistics, in comparison with RANS. Above those thresholds, significant bettering is observed. Concerning Re Q , the value of 5617 seems above, and 3745 unsufficient, whereas, for Re Ω , 16755 seems above the threshold, and 8378 below it.

	summarizes the values of the three Reynolds
	numbers : axial, rotational and effective, Reynolds
	numbers. It is important to note that, despite 6.71 is
	the highest value for N , this case has the lowest ax-
	ial Reynolds number, which is 3745. This value is
	rather low. Recalling that hybrid RANS/LES meth-
	ods are rather designed for predicting high-Reynolds
	flows, one can explain why the case N = 1.49, hav-
	ing the lowest effective Reynolds number, the hybrid
	RANS/LES approach exhibits a bad overall prediction
	quality, for the flow statistics. Concerning other cases,
	the prediction quality of the hybrid RANS/LES ap-
	proach is more contrasted.
	Indeed, a low value of the axial Reynolds number Re Q
	(case N = 6.71) leads to a quite bad description
	of statistics involving the axial component z. Simi-
	larly, a low value of the rotational Reynolds number
	Re Ω results in a particularly poor prediction of statis-
	tics involving the azimuthal component θ. Recipro-
	cally, more satisfactory results are obtained for higher
	zz Reynolds values : e.g. the peaks of R *	1/2

* θ and R * θθ 1/2 , for cases N = 4.47 and N = 6.71. The prediction of R * rr 1/2

Table 1 :

 1 Values of axial, rotational and effective Reynolds numbers, for all four cases.

Acknowledgments

This work was financially supported by Liebherr Aerospace (project "Entrefer moteur") and by the LABEX MEC (ANR-11-LABX-0092).


Brigitte Bidegaray-Fesquet

Marianne Clausel
email: marianne.clausel@imag.fr

Data driven sampling of oscillating signals

Keywords: Level crossing sampling, oscillations, monoHölderian functions. 2000 AMS Mathematics Subject Classification -26A16, 60G15, 60G18, 94A12

HAL is

Introduction

Autonomy, size and weight are very important issues in the design of mobile systems. One possibility to reduce the power consumption of signal processing systems is the reduction of the number of samples. Non uniform sampling is a way to have few samples for a large class of signals, especially sporadic signals, while still describing correctly the active parts of the signal. This leads to a smaller number of samples compared to Nyquist sampling [START_REF] Guan | Opportunistic sampling by level-crossing[END_REF][START_REF] Mark | A nonuniform sampling approach to data compression[END_REF][START_REF] Marvasti | Nonuniform Sampling. Theory and Practice[END_REF][START_REF] Sayiner | A level-crossing sampling scheme for A/D conversion[END_REF]. Specific system architectures, such as event-driven architectures, allow the implementation of this specific sampling. These architectures take samples each time some specific event occurs, e.g. specific voltage levels are crossed. Simple, low power, analog circuits can be designed to acquire information, possibly at high speed.

In this paper we want to relate the signal regularity to the number of non uniform samples obtained via a level crossing technique. Indeed, intuitively, the more the signal is oscillating the more often the signal is sampled. This is of course a local property: the number of samples at the neighborhood of some point may then be related to the local smoothness of the signal, or more precisely to its Hölder regularity. This relationship will be tested on signals whose smoothness properties are perfectly known at each point. This can be useful to predict the processing complexity of e.g. biological signals such as EEG signals or fMRI data which are well-known to be both highly irregular and non stationary.

We introduce here an algorithm which is slightly different from the usual level crossing technique. In [START_REF] Allier | Asynchronous level crossing analog to digital converters[END_REF] the amplitudes are selected thanks a M -bit asynchronous analog-to-digital converter (AADC) that corresponds to 2 M predefined levels in the voltage range. Level crossing sampling consists in taking a sample each time the predefined levels are crossed. Each amplitude has to be associated to a time. More precisely we store a delay elapsed since the last sample was taken, the local clock that enables this is then reset to zero and ready to measure the next delay. In Figure 1 we display the case when the captured time is that of the next clock tick. The sample is displayed with disks, and the value of the signal at the capture times with circles (Figure 1, left). This leads to a few number of samples, especially for sporadic signals. This procedure is refined decimating the samples by keeping only the last one when a level has been crossed many times successively (Figure 1, right).

Our goal here is not to study the approximation of the signal but the number of non uniform samples, given the regularity of the signal, the clock precision and the level quantum. We introduce another sampling algorithm which is slightly different from the AADC, but easier to analyze mathematically, and which yields essentially the same number of samples. In Section 2, we describe this sampling algorithm and rephrase it in mathematical terms. In Section 3 we define the functions that we will use in the numerical experiments of Section 4. These functions are chosen because we are able to control exactly their Hölder regularity.

Figure 1: Non uniform sampling with an AADC [START_REF] Allier | Asynchronous level crossing analog to digital converters[END_REF]. Left: the disks correspond to the samples and the circle to the values at the clock ticks. Right: the disks are the only samples that are kept after decimation.

Algorithm and mathematical interpretation 2.1 Step 1: Generation of an oversampled signal

Even in event-driven systems, where the signal is not sampled at each clock tick, there are clocks that measure time, and specifically the time elapsed since the last event. These clocks have a certain precision, and all measured times are multiples of some basis time t b . Up to some re-scaling of time we suppose that t b = 2 -j , for some j ∈ N. We never have a complete knowledge of the original signal f (t), but only its samples f j,k = f (k2 -j ), for all k ∈ Z (Figure 2). Let V j be the space of continuous functions, which are linear on intervals

I j,k = [k2 -j , (k + 1)2 -j [ for all k ∈ Z.
The Faber-Schauder hierarchical basis, defined in [START_REF] Cohen | Wavelet Methods in Numerical analysis[END_REF], yields a natural basis of V j . Let ϕ(x) = max{0, 1 -|x|}. The functions ϕ j,k = ϕ(2 j • -k), for all k ∈ Z, form the Faber-Schauder basis. We can uniquely define the linear interpolation f j ∈ V j of f at scale 2 -j by imposing f j (k2 -j ) = f j,k , for all k ∈ Z, and

f j = k∈Z f j,k ϕ j,k . (1)
In the sequel we suppose that f is compactly supported in [0, 1] and therefore k = 0, . . . , 2 j -1 in Equation (1).

Step 2: Level crossing

We consider that levels are uniformly spaced by some quantum 2 -M . In applications where the range of the signal is [0, 1], the sample can then be stored with a M -bit register. The second step consists in approximating the samples f j,k by the nearest level below (Figure 3). We denote ⌊x⌋ the integer part of x, namely ⌊x⌋ = sup{n ∈ N, n ≤ x}. The function fj ∈ V j which coincides with the new samples is uniquely defined by fj =

2 j -1 k=0 2 -M 2 M f (k2 -j ) ϕ j,k .

Step 3: Decimation

Next, we decimate the samples so as to keep only one sample when consecutive samples have the same amplitude. We choose to keep the last sample to be compatible with the causality principle (Figure 4). We are interested in the number of samples after the three steps. Comparing Figures 1 (right) and 4, we notice that the number of final samples (the disks in both figures) are comparable. This is a generic situation. In fact the differences are mainly due to the extreme upper and lower levels.

From the mathematical point of view, decimation consists in keeping a subsequence of k = 0, . . . , 2 j -1, defined by induction: k 0 = 0 and

k i+1 = min{k ≥ 1 + k i / 2 M f (k2 -j ) = 2 M f (k i 2 -j ) }.
We only store the couples (δt i , a i ) where δt i = (k i -k i-1 )2 -j , i ≥ 1 is the delay since the last sample, and

a i = 2 -M 2 M f (k i 2 -j ) )
is the amplitude of the sample. Step 3 leads to a reduction of the number of samples, but does not introduce any approximation.

Application to monoHölderian functions

Our goal is to relate the number of non uniform samples to the regularity of the signal. In particular, we address the Hölderian regularity.

MonoHölderian functions

Before introducing the Hölderian regularity, we first recall a few definitions. They allow a weaker definition of pointwise smoothness. The final goal is to define strongly monoHölderian functions, a notion that formalizes the idea of a function which has the as uniformly as possible regularity.

Let x 0 ∈ R and r > 0, for all 0 ≤ h ≤ r we define

B h (x 0 , r) = {x : [x, x + h] ⊂ B(x 0 , r)},
where B(x 0 , r) is the ball of center x 0 and radius r, and denote, as usual, the oscillation of a function f : R → R at x 0 on the ball B(x 0 , r) as

Osc(f )(x 0 , r) = sup |h|≤r f (x + h) -f (x) L ∞ (B h (x 0 ,r)) .
Definition 1 Let f : R → R be a locally bounded function, let x 0 ∈ R and α ∈ (0, 1). The function f is Hölderian of exponent α at x 0 (f ∈ C α (x 0 )) if there exist C and R > 0 such that

Osc(f )(x 0 , r) ≤ Cr α , ∀r ≤ R. (2)
A function f is uniformly Hölderian of exponent α (f ∈ C α (R)) if C and R in Equation (2) are uniform in x 0 ∈ R.
The irregularity of a function can be studied through the notion of anti-Hölderianity.

Definition 2 Let f : R → R be a locally bounded function, let x 0 ∈ R and α ∈ (0, 1). The function f is anti-Hölderian of exponent α at x 0 (f ∈ I α (x 0 )) if there exist C and R > 0 such that

Osc(f )(x 0 , r) ≥ Cr α , ∀r ≤ R. (3)
Let us notice that the statement (3) is stronger than just negating the Hölderian regularity. Indeed such a negation only yields the existence, for any C > 0, of a subsequence (r n ) n (depending on C) for which

Osc(f )(x 0 , r n ) ≥ Cr α n .
Strongly monoHölderian functions naturally arise in the study of the regularity of mappings such as Weierstrass-type or random processes (see e.g. [START_REF] Geman | Occupation densities[END_REF][START_REF] Heurteaux | Weierstrass function with random phases[END_REF]). Indeed, many results only hold for such mappings.

Definition 3 Let α ∈ (0, 1). A function f : R → R is strongly monoHölderian of exponent α (f ∈ SM α (R)) if f ∈ C α (R) ∩ I α (R), i.e.
if there exists C and R > 0 such that, for any x 0 ∈ R,

r α /C ≤ Osc(f )(x 0 , r) ≤ Cr α ∀r ≤ R. (4)

Approximation properties

As already mentioned, only Steps 1 and 2 lead to approximations. To state our approximation results, we need some preliminary definitions. According to our application, we now restrict to functions defined on [0, 1]. For any continuous function f on [0, 1], we define its uniform regularity modulus by

ω f (r) = sup |h|≤r, x/[x,x+h]⊂(0,1) |f (x + h) -f (x)|.
The function ω f is a modulus of continuity in the sense that ω f (0) = 0 and that there exists some C > 0 such that ω f (2r) ≤ Cω f (r) (see [START_REF] Jaffard | Wavelets Methods for Pointwise Regularity and Local Oscillations of Functions[END_REF]).

In what follows we need the notion of strong modulus of continuity introduced in [START_REF] Clausel | Quelques notions d'irrégularité uniforme et ponctuelle[END_REF][START_REF] Jaffard | Wavelets Methods for Pointwise Regularity and Local Oscillations of Functions[END_REF]. The modulus of continuity θ is said to be strong if there exists C > 0 such that for any positive integer J one has

J j=0 2 j θ(2 -j ) ≤ C2 J θ(2 -J ) and ∞ j=J θ(2 -j ) ≤ Cθ(2 -J ).
It is well-known [START_REF] Jaffard | Wavelets Methods for Pointwise Regularity and Local Oscillations of Functions[END_REF] that if there exists some strong modulus of continuity θ such that

ω f (2 -j ) ≤ Cθ(2 -j ), then f -f j L ∞ ≤ Cθ(2 -j ),
where f j is defined by Equation (1). In particular if f ∈ C α (0, 1) there exists some C > 0 such that for any j ≥ 0

ω f (2 -j ) ≤ Cθ(2 -j ),
with θ(2 -j ) = 2 -jα , and then there exists a constant C (which depends on f but not on the scale j) such that

f -f j L ∞ ≤ C2 -jα .
Assume now that in addition

ω f (2 -j ) ≥ θ(2 -j )/C,
then, following [START_REF] Clausel | Quelques notions d'irrégularité uniforme et ponctuelle[END_REF][START_REF] Clausel | A wavelet characterization for the upper global Hölder index[END_REF], there exists a constant C and β > 1 (which depend on f but not on the scale j) such that

f -f j L ∞ ≥ j -β θ(2 -j )/C.
In particular, applying these results with the strong modulus of continuity θ(r) = r α , we deduce that if the function f is assumed to be uniformly monoHölderian with exponent α there exist some C > 0 and

β > 1 such that 2 -jα j -β /C ≤ f -f j L ∞ ≤ C2 -jα .
This yields estimates on the error due to Step 1. The approximation made at Step 2 clearly does not depend on the regularity of function f , and we have

f j -fj L ∞ ≤ 2 -M .

Theoretical number of samples in the case of a monotonous function

If f is a monoHölderian function with exponent α, by definition there exists C 1 , C 2 > 0 and for any scale j ≥ 0 and 0 ≤ k ≤ 2 j -1

C 1 2 -jα ≤ sup (x,h)/[x,x+h]⊂ k 2 j , k+1 2 j |f (x + h) -f (x)| ≤ C 2 2 -jα .
If the function is additionally supposed to be monotonous, we further have exactly sup

(x,h)/[x,x+h]⊂ k 2 j , k+1 2 j |f (x + h) -f (x)| = f k + 1 2 j -f k 2 j .
Hence

C 1 2 j(1-α) ≤ |f (1) -f (0)| = 2 j -1 k=0 f k + 1 2 j -f k 2 j ≤ C 2 2 j(1-α) .
A monoHölderian signal crosses equi-spaced levels with quantum 2 -M at most C 2 2 M +(1-α)j times. The worst case is that of monotonous signals.

Besides, initial sampling (Step 1) takes exactly 2 j samples. This is hence the first natural upper bound for the number of samples. Together with monoHölderianity we know that the number of samples is less than the minimum of these two bounds. For large values of M (or small values of α), we indeed keep almost all of the 2 j samples. Otherwise we can expect some reduction of the number of samples. For C = 1, the threshold is M ≃ αj.

Observe that the proof is based on the fact, that in the monotonous case, we can estimate in a very simple way the oscillations sup

(x,h)/[x,x+h]⊂ k 2 j , k+1 2 j |f (x + h) -f (x)|
of the function. Of course in the general case, the situation can be much more complicated. Nevertheless, generic results in the sense of prevalence as stated in [START_REF] Clausel | Some prevalent results about strongly monoHölder functions[END_REF] are expected to hold. In what follows, we illustrate through numerical simulations what happens in two cases.

4 Numerical simulations

Fractional Brownian motion and the Weierstrass function

We test level crossing on two toy models: sample paths of fractional Brownian motion B H and the Weierstrass function W H , which are indexed by the Hurst index H ∈ (0, 1). The choice of these two cases is guided by the fact that their smoothness properties are related to the Hurst index. The fractional Brownian motion (fBm) B H is the unique Gaussian H-self-similar process with stationary increments. It is defined from its covariance function

E B H (x)B H (y) = 1 2 |x| 2H + |y| 2H -|x -y| 2H
for all (x, y) ∈ [0, 1] 2 . The classical Brownian motion corresponds to H = 1/2. The sample paths of fBm are well-known to be almost surely continuous. Further, its Hurst index H is directly related to the roughness of its sample paths. More precisely the classical law of the iterated logarithm ensures that

B H ∈ C H-ε ([0, 1]) ∩ I H+ε ([0, 1]) almost surely.
Roughly speaking, a.s. for all (x, y) The Weierstrass function W H is a classical example of monoHölderian function with exponent H as proved in [START_REF] Hardy | Weierstrass's non differentiable function[END_REF]. It is defined as

∈ [0, 1] 2 , sup (u,v)∈[x,y] 2 |B H (u) -B H (v)| ∼ |x -y| H .
W H (x) = ∞ j=0 2 -jH cos(2 j x),

Tests

The tests are performed within the SPASS Matlab toolbox [START_REF] Bidégaray-Fesquet | Signal Processing for ASynchronous Systems toolbox[END_REF] (Signal Processing for ASynchronous Systems toolbox). It has been originally designed to treat non uniform signals produced by asynchronous systems, but can be used for a large variety of signals. To generate fractional Brownian motions, we make use of the genFBMJFC.m function [START_REF] Coeurjolly | Simulation and identification of the fractional Brownian motion: a bibliographical and comparative study[END_REF].

We use two values of j (10 and 13) and two values of M (4 and 5). These small values of M are sufficient for most mobile applications. Our output is the number of samples after decimation (Step 3). For the fractional Brownian motion, we perform 1000 realizations and average the number of samples obtained for each realization to obtain an average number n. We perform the same tests on the Weierstrass function (deterministic function, only one realization).

We perform this for values of the Hurst number H in the (0, 1) range and obtain the plots (in semi-log scale, with log-basis 2) in Figures 7 and8 for j = 10 and 13 respectively. We also plot the number of samples computed in the worst case (monotonous function i.e. maximum total variation) for C = 1: N worst = min(2 j , 2 M +(1-H)j ) = 2 j min(1, 2 M -Hj ). We distinguish two regimes: below some value of the Hurst number H ∼ M/j the algorithm more or less keeps all the original samples, above this value the decimation is efficient and yields a significant reduction of the number of samples. For the different curves these "critical" values of H are given in Table 1. The plots associated to fBm are much more regular than those associated to the Weierstrass function because they are obtained by an averaging procedure. Besides, for the Weierstrass function, the constant involved in Equation ( 4) is a priori not equal to 1 and indeed depends on H. It can be explained using the range of the fractional derivative of order H of Weierstrass function W H which is not reduced to a constant and depends on H ( [START_REF] Kui | On the fractional calculus of fractional functions[END_REF][START_REF] Samko | Fractional integrals and derivatives. Theory and applications[END_REF][START_REF] Zähle | Fractional Derivatives of Weierstrass-Type Functions[END_REF] for more details).

Conclusion

We have predicted for monoHölderian functions and shown numerically that there is a strong relationship between the smoothness properties of a signal and the number of samples that can be obtained by the crossing level algorithm presented in this paper. This is rigorously proved in the case of monotonous monoHölderian functions. The next step, which will be the purpose of a forthcoming paper, will then be to consider signals whose regularity may change from point to point such as multifractional or multifractal signals.

Figure 2 :

 2 Figure 2: Regular sampling of the input signal (Step 1).

Figure 3 :

 3 Figure 3: Reduction to predefined levels (Step 2). Samples from Step 1 are disks and the new samples are the circles.

Figure 4 :

 4 Figure 4: Non uniform samples after decimation (Step 3).

Figure 5

 5 Figure 5 presents three realizations of sample paths of fractional Brownian motions for H = 0.5, H = 0.7, H = 0.9, and j = 10 (1024 samples).

Figure 5 :

 5 Figure 5: Three realizations of fractional Brownian motions for H = 0.5, H = 0.7, and H = 0.9 (from left to right).

  and, for all (x, y)∈ [0, 1] 2 , sup (u,v)∈[x,y] 2 |W H (u) -W H (v)| ∼ |x -y| H .

Figure 6

 6 Figure6presents the graphs of the Weierstrass functions for H = 0.5, H = 0.7, H = 0.9, and j = 10 (1024 samples).

Figure 6 :

 6 Figure 6: The Weierstrass function for H = 0.5, H = 0.7, and H = 0.9 (from left to right).

Figure 7 :

 7 Figure 7: Number n of samples in terms of the Hurst number in the log scale for j = 10, and M = 4 (left) and M = 5 (right). Solid lines correspond to the averaged number for the fractional Brownian motion, the dashed lines to the Weierstrass function, and the dotted lines to the worst case j + max(1, M -Hj).

Figure 8 :

 8 Figure 8: Number n of samples in terms of the Hurst number in the log scale for j = 13, and M = 4 (left) and M = 5 (right). Same plotting conventions as Figure 7.

Table 1 :

 1 "Critical" values of the Hurst number.

ACKNOWLEDGEMENT

LJK is partner of the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir and this work has been partially supported by the MathSTIC Project OASIS of the Grenoble University. The authors wish to thank Jean-François Coeurjolly for fruitful discussions.