Jonathan Passerat-Palmbach

Jonathan Caux

Yannick Le Pennec

Romain Reuillon

Ivan Junier

François Kepes

David R C Hill

§ Murex

Parallel stepwise stochastic simulation: Harnessing GPUs to Explore Possible Futures States of a Chromosome Folding Model Thanks to the Possible Futures Algorithm (PFA)

Keywords: GPU, parallel simulation, chromosome folding, Possible Futures Algorithm, PFA, stochastic simulation, stepwise simulation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Two major design choices are available for code parallelization. One possibility is to adapt the existing application to fit a parallel template, which mainly consists in breaking loops to spread their workload across several parallel computing elements, be they threads or processors. This approach tends to be limited by the so called Amdahl's law [Amdahl, 1967, Hill and[START_REF] Hill | [END_REF] that computes the maximum speed-up of an application, taking into account the part of it that can effectively be parallelized. Another possibility consists in redesigning the application from scratch, to use the most possibility of the target parallel platform. In this case, the application is more likely to benefit from performance gain laid down by the Gustafson's law [Gustafson, 1988]. Both situations present advantages: adapting an existing application will make it easier to write; redesigning an application from scratch often opens new parallelization perspectives.

For instance, when dealing with particular hardware accelerators such as GPUs, which are constrained by their underlying architecture, it is usually more efficient to provide a dedicated parallel implementation. Such an approach allows this implementation to better fit the platform's requirements. In the case of GPUs, applications are expected to favour computing intensive algorithms rather than those containing a high number of branches instructions.

Parallelization of an existing code, even from scratch, can be severely limited by the sequential nature of the underlying algorithms, which have been originally thought for Turing-like architectures. In this context, it may be useful to redesign the algorithmic principles themselves, just as new algorithmic schemes have emerged from the possibility of quantum computing. In this article, we present the redesign of a traditional model of polymer folding, which has been extensively used in the field of chromosome modelling [START_REF] Langowski | Computational modeling of the chromatin fiber[END_REF]. These models are known to suffer from poor efficiency when considered in situations that are similar to in vivo because of the strong confinement of chromosomes. So far, "simple" code parallelization has mainly been set up using tools such as OpenMP [START_REF] Dagum | Openmp: an industry standard api for shared-memory programming[END_REF] or MPI [Message Passing Interface Forum, 1993]. More recently, parallelization of dynamical models have been realized on GPU architecture, demonstrating the usefulness of such technology for the modeling of long polymers [START_REF] Reith | Gpu based molecular dynamics simulations of polymer rings in concentrated solution: Structure and scaling[END_REF]. In the same spirit, our algorithm is mainly meant to exploit massively parallel processors, more particularly NVIDIA CUDA-enabled GPUs.

Our algorithm is based on a two-step procedure: i) generation of several potential future states, and ii) selection of the next state as a random choice amongst the set of valid futures. We hence call this approach the "Possible Futures Algorithm" (PFA). PFA can be compared to the branch prediction feature available in modern CPUs, which enables CPUs to load the instructions following the branches of a conditional statement, prior to knowing the result of this statement. This avoids offloading the pipeline mechanism, and consequently waste precious instruction cycles [Lee andSmith, 1984, Smith, 1998]. As a result, we show that PFA highly increases the acceptance rate of the new model. While the initial algorithm only accepts ∼ 5% of the proposed states in certain configurations, our solution maintains a satisfying acceptance rate over 60% in the same conditions, regardless of the input parameters of the model.

The paper is organized as follows. First, we describe the initial chromosome folding model on which our parallel declination is based. Next, we point out the limitations that led us to design a new model. We then introduce our new parallel chromosome folding model and the PFA approach. Finally, we discuss some implementation choices and study the performances of our algorithm.

A classical Monte-Carlo simulation of polymer models

Large-scale properties of polymers have been efficiently captured thanks to simplified models that ignore the details of both the polymer composition and the exact nature of the surrounding solvent [ de Gennes, 1988]. In biology, various models have been used to model the behavior of DNA molecules and of chromosomes [START_REF] Strick | Stretching of macromolecules and proteins[END_REF]. Chromosomes have been modeled as simple flexible chains that cannot overlap with themselves (self-avoidance effect) -see Figure 1. These so-called worm-like chains provide a coarse-grained description of protein-coated DNA that is simple enough so that, using Monte-Carlo simulations, they can be used to investigate the folding properties of rather long biomolecules [START_REF] Langowski | Computational modeling of the chromatin fiber[END_REF]. Thus, they have been used to address both the problem of the structuration of chromosomes in vivo [START_REF] Fudenberg | Higherorder chromatin structure: bridging physics and biology[END_REF] and, more recently, the interplay between structuring and function [START_REF] Junier | CTCF-mediated transcriptional regulation through cell type-specific chromosome organization in the -globin locus[END_REF].

The worm-like chain is a continuous polymer model that needs to be discretized to be simulated on a computer. To this end, the polymer chain is divided into contiguous cylinders. Two consecutive cylinders are jointed together so that they can rotate freely around their joint (Figure 1), which is used to make the polymer conformations evolve dynamically. To mimic the persistence of the directional order of the polymer chain (see Figure 1), the junctions between each cylinder (hereafter called joints) carry an energy called "bending energy". This energy is exchanged with the solvent using a Monte-Carlo algorithm (see below) and depends on the angle between the two cylinders connecting the joint. We further consider chromosomes that are confined into cells that are more or less narrow, so that our approach can address the computational properties of chromosomes in confined geometries [ de Gennes, 1988]. In the following, we call "cell" the embedding volume of the chromosomes during the simulations, i.e. in silico. In the most general case, this may not correspond to the nucleus in which chromosomes are embedded in vivo.

In summary, a chromosome in silico can be viewed as a long string of contiguous "cylinders". The chromosome can be circular (as in bacteria) or linear (as in human). Its statistical properties are the result of its jiggling motion due to the interaction with the surrounding solvent. In these models, the solvent Figure 1: Worm-like chain model of chromosomes. On the top is indicated the original continuous model, which is characterized by the length beyond which the chromosome looses the memory of its directional order (persistence length l p ). The actual situation to be modelled is a situation where the chromosome is confined to a small volume (on the right). To be simulated on a computer, this model is discretized into cylinders (bottom). In this discretized version, two consecutive cylinders can rotate around their common joint but have to pay some energy to do so because as they are semi-flexible (captured by the persistence length). In the simulations, we use five cylinders per persistence length (figures adapted from [START_REF] Junier | Spatial and topological organization of dna chains induced by gene co-localization[END_REF]).

is not taken explicitly. Instead, we used a standard stochastic algorithm, also called a Monte-Carlo procedure, to make the chain move step by step. Each of these steps is conceptually simple, and altogether, they form Algorithm 1.

We used the standard Metropolis rate for accepting the rotation depending on its energy, i.e. if the new bending energy is lower, the transition is accepted; otherwise, the transition is accepted with probability e -∆E/k B T where ∆E stands for the increase of energy and k B T for the thermal energy characterizing the solvent (k B is the Boltzmann constant and T is the temperature). This dynamics is known to satisfy the required properties to simulated the thermodynamic behaviour of polymer models. In practice, it favours movements whose energy variation is dictated by the value of the temperature T . As a result, the chromosome is prompted to a state of equilibrium between the energies of the conformations and the number of equivalent conformations at a given energy, also called entropy.

This stochastic algorithm slowly makes the chromosomes evolve following a combination of random and physical constraints, over hundreds of millions of runs. The limiting factor for this evolution happens when most rotations are Algorithm 1 Monte Carlo procedure involved in the simulation process Select a section of consecutive cylinders of random size (called "block" in the following); Perform a rotation of the block by a random angle around the axis that intercepts the first and last joint of the block; if The rotated block does not overlap with the rest of the polymer then Compute the new bending energy of the joints at the edges of the block; Accept the rotation with some probability depending on the variation of the bending energy; end if Increase time +1; rejected due to the overlapping conditions, which occurs very often when the embedding cell is small, as in biological situations. In particular, simulations of large chromosomes become extremely time-consuming, which make this type of algorithms poorly efficient to tackle in detail the folding problems of multiple chromosomes in vivo.

3 Limitations of the sequential model

The collisions bottleneck

The original algorithm picks up a "block" at random on the chromosome, and randomly rotated it around an axis. And so on. By definition, it is therefore sequential. As a consequence, for string confinement constraints, i.e. for small cell sizes, it is difficult for the algorithm to find valid configurations towards which to evolve. Thus, because valid configurations are hard to find through random attempts, a large number of (wasted) steps is required to provide a new valid configuration. In situation of string confinement, 95% of the generated rotations can produce collisions (see below). It is consequently a major bottleneck that needed to be tackled in order for the model to scale up with chromosomes sizes. A first possibility would be to reduce the amplitude of the rotations. However, in this case, motions of the polymer are very small and both a large number of blocks and a large number of tries are necessary to thermalize the system (i.e. to make it visit all most likely conformations). Thus, in any case, the use of parallel motions can facilitate the dynamical evolution of the polymer.

The Possible Futures Algorithm (PFA)

We propose a parallel simulation model which considers two parallel optimizations. The first one consists in handling several blocks in parallel, instead of only one block as in the sequential algorithm. Processing blocks in parallel is possible in our case because it does not break any physical law at the heart of the model. In particular, it does not break the so-called ergodicity, i.e. the unbiased sampling of the configuration space. This new parallel rotation step is presented in Figure 2. The second optimization consists in considering for each single block a set of several possible rotations to maximize the chances of producing a correct rotation. We call each of these rotations a "possible future" of a block. In Figure 3, we can see a chain of cylinders called "Base". It represents the chromosome as it is at the beginning of an iteration. Yellow cylinders belong to different blocks labelled Block 0 through Block 2. Starting from blocks defined on the base chromosome, the algorithm will generate several possible futures for each block, applying each time a different rotation angle to the block. In Figure 3, three futures are generated for each of the three blocks defined.

Cylinders in grey do not belong to any block, and therefore remain static during this iteration. It is important to note that at each step, the current state of a given block is also viewed as a possible future. Consequently, the fact of not rotating a block is valid, and is always a possible solution for the dynamical evolution of the polymer (which is here a no-motion).

Parallel model

The Possible Futures Algorithm induces extra computation time, since futures need to be created before each step of the simulation, and merged as a complete new state for the whole chromosome at the end of each simulation step. However, such an organization reminds of a frequently used pattern in parallel computing: map/reduce [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF]. This pattern is particularly adapted to GPU computing and some frameworks can even simplify its implementation Figure 3: Example of generation of three possible futures for three different blocks on GPU [START_REF] He | Mars: a mapreduce framework on graphics processors[END_REF]. In our case, the map stage consists in generating and processing new futures; this stage is processed in parallel, while the reduce stage is performed sequentially to build the resulting state of the chromosome at the end of the current simulation step. Next section will detail the stages of our parallel simulation model.

We describe now the consecutive steps allowing the parallel algorithm to generate and transform possible futures, before combining them as a whole new state for the chromosome. In view of the massive number of repeated parallel tasks that result from the possible future approach, we have chosen to exploit GPU accelerators to run this algorithm.

Generate possible futures

First, the map stage must be initiated sequentially by generating random blocks in the chromosome. Blocks are equally sized arrays of consecutive joints from the chromosome. Each block owns a different, non-overlapping chain of joints. Blocks also differ from the random angle describing the rotation which will be applied to them. Blocks cannot overlap to prevent the same cylinder from being involved in two different rotations.

To generate futures, a "master" thread runs alone on the GPU to define random blocks and to generate random rotation angles. Then, a burst of threads is launched to clone the original blocks and change their rotation angles, thus making them into possible futures of the original block. Upon completion, possible future blocks are ready to be processed in parallel by thousands of GPU threads.

Calling F the number of possible futures per block and N the number of blocks, one can see that this stage operates in parallel the equivalent of what the sequential algorithm used to do over F × N time step. The potential performance gain is hence on the order of F × N since while the sequential simulation checks after each step that the block can be moved, the parallel implementation fully transforms each possible future without wasting time on intermediate verifications.

Select valid futures

Once all the futures have been transformed and represent a potential future state of their base block, it is necessary to filter them: i.e. dismiss those that are not valid. The algorithm calculates the bending energy carried by the joints at the edges of a future block. To be stated as valid the resulting bending energy of a future must be lower than the one of its original block. Should the new bending energy be greater than the original one, then the considered future has gained energy. In this case, there is a probability that the future may be valid, determined by a random draw according to the physical properties of chromosomes. In any other case, the physical properties of the future make it unsatisfying, and it is consequently withdrawn. The remaining futures must still meet two other requirements, referred to as static and dynamic conditions, to be compatible with their environment and to be potentially included in the global solution.

Static conditions are met when all cylinders of the future are within the virtual cell in which the chromosome must remain confined. Moreover, none of the cylinders of the future can collide with parts of the chromosome that did not evolve during this simulation step. By doing so, we eliminate all futures that are not physically possible.

For dynamic conditions to be valid, a second evaluation determines which futures are suitable with respect to the other futures. A future is suitable as long as it does not collide with any of the other valid future generated at this step of the simulation.

Futures acceptance tables are displayed in green in Figure 3, while dismissed futures appear in red.

To sum up, we have introduced extra operations when delaying the selection of valid futures, but this stage offers many potential solutions to increase the acceptance rate of rotations that make the chromosome evolve. This stage justifies on its own the use of a massively parallel architecture such as GPUs to speed-up the computation of the algorithm. Let's now study the operations which compose the reduction stage of a simulation step.

Determine compatible futures

We consider a two-round reduction in our algorithm. The first round models the problem as a graph to determine the compatible futures. The vertices of the graph represent possible futures, and the presence of an edge between two vertices represents the compatibility between two possible futures.

Compose the global result

The second reduction step now aims at composing a global result from the compatible futures identified by the graph. From this graph hundreds of random cliques are built from the whole set of possible futures at hand. A constant number of vertices is simply picked at random in the graph to build equallysized cliques, with the particular care to always draw a vertex from each group of possible futures.

We now compose a global result by finding one maximum clique from the set of cliques issued previously. A maximum clique owns as many possible block futures as there are blocks in the original chromosome. Thus, any maximum clique forms a solution to our problem and its nodes are the new state of each block defined at the beginning of the simulation step. Please note: should several maximum cliques be found during the first round of the reduction, one would be chosen at random, since they all represent equivalent satisfying solutions to make the chromosome evolve.

The elements from Figure 3 are displayed as a graph in Figure 4. It shows that the possible future 0 from Block 1 is suitable in terms of physics criteria, but is not compatible with any possible future from Block 2 (Although this future is compatible with the current state of the chromosome, for the sake of simplicity, it is not shown in Figure 4).

In Figure 4, we see a clique compounded of 3 vertices: {Block0.Future2, Block1.Future1, Block2.Future1}. This means these three are inter-compatible futures, and together form a global consistent future for the chromosome. To validate the whole operation, each possible future block is stored in the base chromosome, replacing its original counterpart. GPU applications often suffer from the excessive latency induced by memory accesses. It is strategic to optimize them as much as possible in order to obtain the best performance (memory transfers from the host to the device, or memory accesses from threads on the device) [START_REF] Ryoo | Optimization principles and application performance evaluation of a multithreaded gpu using cuda[END_REF].

In our case, we have decided to cut any communication between the host and the device, except synchronizations. Indeed, GPUs cannot have all their threads synchronized at one point, at the time of writing, due to hardware constraints. On the other hand, some stages of the simulation need the data to be in a particular state before they can achieve their part of the computation. The only way to synchronize all the threads on a GPU is to wait for the completion of the kernel that has launched them. Each simulation step is consequently divided into several parallel stages, launched consecutively on the GPU. Such a design enables us to finely tune the parallelism grain of the simulation: at each stage the parallelism grain is the most adapted to the kernel that is to be launched.

As a result, communications are shrunk to their utter minimum from the host to the device. Simulations data are transferred back to the host from the GPU only when a checkpoint of the state of the simulation is reached.

Approximate cylinders to avoid branch divergence

In the original model, the elements that compose the chromosomes are represented by cylinders with rounded edges, as shown in Figure 5. In this diagram, the point mid is the middle of the cylinder, and e1 the direction vector that indicates the orientation of the cylinder. rin and rfin are two points which denote the edges of the cylinder. This representation has the advantage of being accurate according to the observed physics, but on the other side, it highly complicates the calculation of collisions. Indeed, computing an intersection between two cylinders is rather slow: the corresponding code is about 200 lines and contains many conditions, which makes it less suitable for execution on GPUs.

Instead, cylinders can be approximated by groups of spheres, as illustrated in Figure 6. Points sphere0, sphere1 and, sphere2 are the centres of the three spheres composing the cylinder. Points joint0 and joint1 are the edges of the cylinder. Vector3 v e c t o r = t h i s ->P o s i t i o no t h e r S p h e r e . p o s i t i o n ; // A l l s p h e r e s have a r a d i u s RADIUS . // We compare t h e s q u a r e d l e n g t h s // t o a v o i d an e x p e n s i v e s q u a r e r o o t . r e t u r n v e c t o r . squaredLength ( ) -(RADIUS ✯ RADIUS ) ; } Three spheres (rather than two or four) were chosen because it allows a very good approximation of the original cylinders and collisions observed. Moreover, using more than three spheres does not increase the accuracy of the approximation in a significant way.

Results

Thanks to the log files generated by the application, which contain the 3D coordinates of each joint, we can obtain a 3D representation of chromosomes using third-party software such as gnuplot. Figure 7 is the output drawing generated by gnuplot for a circular chromosome after 100 million iterations. All the results presented in this section have been obtained through simulations run on homogeneous machines equipped with the following configuration: Intel Core i7-2600k 3.40GHz CPUs and NVIDIA GeForce GTX590 GPUs.

Efficiency of the parallelization approach

The first metric that is interesting to study from the new parallel model is the efficiency of the possible futures approach regarding the acceptance rate. Figure 8, shows the efficiency of the PFA approach in increasing the acceptance rate. As we can see, the acceptance rate increases with the number of possible futures involved in each simulation step.

For clarity's sake, the acceptance rate depicted in Figure 8 is obtained by averaging the acceptance rates of several configurations benefiting successively from 4, 8, 64 and 128 possible futures per block. The configuration test set is as follows:

❼ number of blocks → {1; 3; 7; 15} ❼ maximal block size → {32; 64; 256; 510} Eventually, it appears that the acceptance rate becomes trickier to improve when the number of blocks increases. Indeed, the number of cliques grows exponentially with the number of blocks. Thus, it is more and more difficult for the algorithm to find a maximal clique amongst the set of randomly generated cliques.

Performance on the cutting-edge Kepler architecture K20 GPU

Parallel architectures and especially GPUs evolve very quickly. In the particular case of NVIDIA GPUs, a new generation is shipped every two years. Thus, an application running on this architecture must take advantage of these evolutions and display improved performance when run on the forefront GPUs.

Having run our benchmark on Fermi architecture GTX 590 GPUs, we have tested some identical configurations on a recent Kepler architecture K20 GPU [NVIDIA, 2012]. Figure 9 depicts the speed-up obtained at no cost, only by switching from the GTX to the K20 GPU.

Conclusion

In this paper, we have described a parallel model of chromosome folding and its implementation on GPU, using the NVIDIA CUDA technology.

The purpose of this model is to make chromosomes evolve in the confined space of a cell, by applying successive rotations on random chunks of the chromosome. This new parallel model is based on an initial sequential one, which encounters difficulties when the cell space becomes more confined. The acceptance rate of the simulated rotations then quickly falls, making the sequential model evolve slower.

In order to increase the acceptance rate of simulated rotations, especially when the space becomes more and more confined, we designed a Possible Futures Algorithm (PFA). This approach proposes different new configurations for each block at each step, which have to be checked for inter-compatibility. All the futures can be generated and processed in parallel, thus inducing a massively parallel workload of repetitive tasks.

This approach has shown satisfying results in terms of acceptance rate. While the original model was likely to dismiss about 95% of the rotations, the possible future approach displays a satisfying acceptance rate of more than 60%. Such an acceptance rate enables the chromosome to evolve to a new configuration after almost each simulation step.

We now plan to challenge this parallel model with an implementation of the sequential one. Such a benchmark will have to be conceived very carefully, as long as the two models operate in two different ways and the same metrics cannot be easily compared between the two models. For instance, a simulation step does not perform the same actions in the two models. That is why we focused on the acceptance rate in this study. Still, the parallel implementation should display far better performance than the sequential one thanks to its possible futures approach. The more the chromosome becomes confined within the cell space, the more the parallel algorithm will be efficient. In fact, the efficiency of the parallel implementation is due to its ability to maintain a high acceptance rate throughout the execution, whereas the sequential model will see its acceptance rate decrease with the various stages and the confinement.

The interested reader will find a free version of both sequential and parallel implementations in the git repository of the French Institute of Complex Systems1 .

  Figure 2: A parallel rotation step

Figure 4 :

 4 Figure 4: Compatibility graph of the possible futures

Figure 5 :

 5 Figure 5: Cylindrical element that is part of the chromosome

Figure 6 :

 6 Figure 6: Pseudo-cylindrical approximation of a segment Thanks to this representation, calculating the presence of a collision between two cylinders consists in calculating the collisions between the three spheres of the first cylinder and the three spheres of the second cylinder. Furthermore calculating collisions between two spheres is a trivial operation, as shown in Listing 1: Listing 1: Calculation of collisions between spheres

Figure 7 :

 7 Figure 7: Gnuplot 3D representation of a chromosome

Figure 8 :

 8 Figure 8: Increase in the acceptance rate with the number of possible futures

Figure 9 :

 9 Figure 9: Model improved performance when run on a Kepler architecture K20 GPU

https://forge.iscpif.fr/projects/dna/repository

Acknowledgments

Ivan Junier is supported by a Novartis grant (CRG). This work is also partly funded by the Auvergne Regional Council and the FEDER. The authors would like to thank Alice Chalet for her careful reading and useful suggestions on the draft.

partially funded by the Auvergne Regional Council and the FEDER


R Oguic

S Viazzo

S Poncet

Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow

Introduction

Turbulent flows in an opened Taylor-Couette system with an axial throughflow is studied here by the means of large eddy simulations. The ultimate industrial application is the effective cooling of the rotor-stator gap of an electrical motor. The recent review of Fénot et al. [START_REF] Fénot | A review of heat transfer between concentric rotating cylinders with or without axial flow[END_REF] pointed out the lack of reliable data in such configurations for the hydrodynamic fields and most of all for the thermal one. Among others, 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] and Escudier and Gouldson [START_REF] Escudier | Concentric annular flow with centerbody rotation of a Newtonian and a shear-thinning liquid[END_REF] provided detailed velocity measurements in very elongated middle-gap cavities. In a recent paper, Poncet et al. [START_REF] Poncet | Numerical modeling of fluid flow and heat transfer in a narrow Taylor-Couette-Poiseuille system[END_REF] used a second-order turbulence model assuming the flow as being steady and axisymmetric, whose results compared quite well with the ones of [START_REF] Escudier | Concentric annular flow with centerbody rotation of a Newtonian and a shear-thinning liquid[END_REF]. Little is known about the near-wall turbulent structures when an axial Poiseuille is superimposed. We can cite the numerical work 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], who highlighted the destabilization of the near-wall turbulent structures due to the rotation of the inner wall giving rise to an increase of sweep and ejection events. The main objective of the present work is to validate the present LES approach against reliable data available in the literature [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF][START_REF] Nouri | Flow of newtonian and non-newtonian fluids in a concentric annulus with rotation of the inner cylinder[END_REF], which can be seen as a first step towards the simulation of the real flow conditions.

Geometrical configuration and flow control parameters

The fluid is confined between two concentric cylinders of height h (see Fig. 1). The inner cylinder of radius R 1 is rotating at the rate Ω , while the outer one of radius R 2 is stationary. The cavity may be characterized by two geometrical parameters: S. Viazzo • S. Poncet • R. Oguic Laboratoire M2P2 UMR 7340 CNRS -Aix-Marseille Université, Marseilles -France, e-mail: {viazzo,poncet,oguic}@l3m.univ-mrs.fr 1 its aspect ratio Γ = h/(R 2 -R 1 ) and its radius ratio η = R 1 /R 2 , fixed here to 9 and 0.5 respectively. An axial throughflow is imposed within the gap at a constant inlet axial velocity W m . The hydrodynamic fields are governed by the rotational and bulk Reynolds numbers defined as

Re Ω = Ω R 1 (R 2 -R 1 )/ν and Re Q = 2W m (R 2 - R 1
)/ν = 8900 (ν the fluid kinematic viscosity). The rotational Reynolds number may vary such that the rotation parameter N = Ω R 1 /W m will take the values 0.2145, 0.429 and 0.858.

Numerical method

The numerical method is based on the work of Abide and Viazzo [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF], who developed a 2D compact fourth-order projection decomposition method in cartesian coordinates. It has been recently validated in the case of turbulent rotor-stator flows in cylindrical coordinates by Viazzo et al. [START_REF] Viazzo | High-order Large Eddy Simulations of confined rotor-stator flows[END_REF]. The time advancement is second order accurate and based on the explicit Adams-Bashforth scheme for the convective terms and an implicit backward Euler scheme for the viscous ones. The derivatives are approximated using fourth order compact formula in the radial and axial directions. The time splitting scheme is an improved projection method ensuring the incompressibility at each time step. The projection decomposition method is based on a direct non-overlapping multidomain Helmholtz/Poisson solver, which provides the solution of each Helmholtz/Poisson problem. The multidomain solver ensures All the numerical details used in the present work are summed in Table 1 and compared to the LES calculations 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]. It is noticeable that the domain is decomposed into 2 subdomains in the axial direction. In order to save computational ressources, the solution is calculated by assuming a π/2-periodicity, though the grid resolutions are given for a π-periodicity to enable direct comparisons with [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF]. After reaching the statistical convergence (about 2000 s of computed physical time), the mean and rms quantities are sampled the last 80 s in each case and averaged both in the tangential and axial directions. All the calculations have been performed on the M2P2 cluster composed of 2 xeon quadcore 3 GHz.

Results

The present LES results are compared to the velocity measurements 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] and to the LES 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] for three values of the rotation parameters N = 0.2145, 0.429 and 0.858.

Flow structures

To highlight the presence of near-wall structures, Figures 2b-d present isovalues of the Q-criterion colored as a function of the radius. Whatever the value of the rotation parameter N, coherent structures are well aligned with the axial flow along the stator. Rotation of the inner cylinder tilts the coherent structures appearing as spiral rolls and forming a negative angle (as they roll up in the opposite sense of the rotor) with the axial direction. These spirals are not observed in 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] for N = 0.429.

Mean flow field

Figure 2 shows the radial distribution of the mean axial and tangential velocity components normalized respectively by W m and Ω R 1 . The axial velocity is weakly influenced by the rotation parameter N and the profile resembles the turbulent Poiseuille profile usually observed in pipe flows. The tangential velocity profile is similar to the Batchelor profile observed in interdisk rotor-stator flows with two thin boundary layers developed on each disk separated by a central region, where V θ slightly decreases with the radius. The present results are in good agreement with 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] for N = 0.429 and improve the previous ones of [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF]. The variation of V θ with 1/r is recovered confirming that the angular momentum is almost constant along a radius. The average friction coefficients C f and the friction Reynolds numbers Re τ obtained along both cylinders compare very favorably to 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] as listed in Table 1.

Turbulent flow field

The radial variations of the R θ θ and R rz components of the Reynolds stress tensor normalized by W 2 m are displayed in Figure 2. Velocity fluctuations are particularly higher for the tangential component with very intense peak values close to the walls. For the two other normal components (not shown here), the same behavior is observed with almost constant turbulence intensities within the gap, which vanish towards the cylinders. As expected, increasing the rotation parameter N (ie the rotation rate Ω ) induces higher turbulence levels especially for the radial and tangential components. Furthermore, the distributions of the normal components are asymmetric, which may be attributed to the destabilizing effect of the centrifugal force. The R rz shear component of the Reynolds stress tensor is shown in Figure 2. When N increases, the shear stresses increase in strength too. A good quantitative agreement has been obtained between our present LES results and the 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]. It slightly improves the LES results 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 the high order schemes.

Figure 3 presents the anisotropy invariant map for the Reynolds stress tensor. It is noticeable that the present LES results respect the realizability diagram of Lumley. As expected, turbulence is mainly at two-component close to the cylinders and tends to the isotropic state within the center of the gap. There is only a weak effect of the rotation parameter, the flow getting more isotropic in the center region for high N values. Our results and the ones of [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF] do not support their conclusions about a "disk-like" state at mid-radius. Even if the third invariant gets weakly negative, all the normal components of the Reynolds stress tensor are indeed all of the same order.

Conclusion

The LES code developed by Viazzo et al. [START_REF] Viazzo | High-order Large Eddy Simulations of confined rotor-stator flows[END_REF] has been extended here to its multidomain version using the methodology described in [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF]. This new approach has been fully validated against the velocity measurements 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] and the numerical results of Chung et al. [START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF] in an opened Taylor-Couette system with an axial throughflow. Near-wall structures appearing as elongated spiral rolls are observed along the inner cylinder. A very good agreement has been obtained for both the mean and turbulent fields and for the three values of the rotation parameter considered highlighting the importance of the order of the spatial schemes. Some numerical developments are still required to simulate the whole problem of the effective cooling in the rotor-stator gap of an electrical motor.

Fig. 1

 1 Fig. 1 (a) Schematic representation of the Taylor-Couette-Poiseuille system with relevant notations; Isosurfaces of the Q-criterion (Q = 0.5) for (b) N = 0.2145, (c) N = 0.429 and (d) N = 0.858

Fig. 2

 2 Fig. 2 Radial distributions of the mean tangential and axial velocity components and of the two components of the Reynolds stress tensor for the three values of N. Present LES (straight lines), LES of Chung and Sung [2] (dashed lines) and measurements of Nouri and Whitelaw [5] (symbols)

Fig. 3

 3 Fig. 3 Anisotropy invariant map for the Reynolds stress tensor: N = 0.2145 (circles), N = 0.429 (triangles) and N = 0.858 (squares)

Table 1

 1 Numerical details and mean flow parameters. The number of mesh points N r in the radial direction is fixed to 65 in all calculations. Comparisons with the LES of[START_REF] Chung | Large-eddy simulation of turbulent flow in a concentric annulus with rotation of an inner cylinder[END_REF]

	Parameters [2]	present	[2]	present	[2]	present
	N	0.2145	0.2145	0.429	0.429	0.858	0.858
	∆ z + i ∆ z + o ∆ r + i ∆ r + o (R 1 ∆ θ ) + 8.01 22.95 -0.27 0.23 (R 2 ∆ θ ) + 13.86	21.36 19.4 0.47 0.43 7.34 13.33	16.8 -0.3 0.24 8.8 14.34	15.87 13.41 0.33 0.45 8.22 13.9	9.93 -0.35 0.26 10.4 15.8	9.31 7.05 0.63 0.47 9.7 14.7
	(N θ , N z ) (128, 128) (128, 130) (128, 192) (128, 194) (128, 384) (128, 386) δt(s) -9 × 10 -3 -7 × 10 -3 -4 × 10 -3
	CFL	-	0.2	-	0.22	-	0.25
	CPU time -	5.2 (s/it)	-	9 (s/it)	-	25.5 (s/it)
	C f ,i
	C f ,o	-	7.44 × 10 -3 -	8.05 × 10 -3 -	9.05 × 10 -3
	Re τ,i	163.21	149.53	179.24	167.52	211.85	197.51
	Re τ,o	141.19	135.8	146.1	141.6	160.91	149.71

8.91 × 10 -3 8.91 × 10 -3 9.86 × 10 -3 1.08 × 10 -2 1.2 × 10 -2 1.58 × 10 -2


Sébastien Poncet

Stéphane Viazzo

Thermo-hydrodynamic instabilities in a high aspect ratio Couette-Taylor system using Direct Numerical Simulation

Thermo-hydrodynamic instabilities in a high aspect ratio

The numerical method is based on the 2D compact fourth-order projection decomposition method of Abide and Viazzo 2 , extended to cylindrical coordinates on non-staggered grids. The time advancement is second-order accurate. The derivatives are approximated using fourth-order compact formula in the radial and axial directions and Fourier series in the azimuthal direction. Due to the large aspect ratio, the domain is axially decomposed into 8 subdomains using the influence matrix technique. The thermal effects are considered using the Boussinesq approximation.

The aim of this work is to characterize all the flow patterns occuring in the system for various Rayleigh and Taylor numbers. At low Ra, the first instability appears in the form of regular spirals (Fig. 1a). For higher values of Ra, dislocations may occur (Fig. 1b) and when one increases T a, spirals coexist with a wavy vortex flow (Fig. 1c) and finally only the wavy vortex flow subsists at high T a (Fig. 1d). The spatiotemporal properties of these patterns have been fully determined and compared to the work of Guillerm

  Couette-Taylor system using Direct Numerical Simulation S. Poncet, S. Viazzo * The influence of a radial temperature gradient on the stability of Couette-Taylor flow has been investigated by direct numerical simulation. The cavity is characterized by a narrow gap η = R i /R o = 0.8 and a high aspect ratio Γ = H/(R o -R i ) = 80, R i and R o being the radii of the inner and outer cylinders and H their height. The flow depends on two physical parameters: the Taylor number T a = ΩR i d/ν(d/R i ) 1/2 and the Rayleigh number Ra = α∆T d 3 /(νκ), with d = R o -R i and Ω the rotation rate of the inner cylinder. 15 cases have been computed for T a = [11 -150] and Ra = [1166 -13228] to enable comparisons with the experiments of Guillerm 1 .

Figure 1 :

 1 Figure 1: Temperature maps at (r-R i )/d = 0.8 in the (θ, z) plane for: (a) Ra = 2063, T a = 50, (b) Ra = 7150, T a = 40, (c) Ra = 7150, T a = 75, (d) Ra = 7150, T a = 120. 61 × 64 × 488 mesh points are used in the (r, θ, z) directions with δt = 0.002 s.


N Peres

S Poncet

E Serre

A 3D pseudospectral algorithm to simulate rotating flows in cylindrical cavities

When simulating flows in cylindrical rotating cavities a difficulty arises from the singularities appearing on the axis. Its singularities are due to the presence of terms 1/r n (n = 1, 2) in the Navier-Stokes equations, where r is the radial distance. To avoid evaluating differential equation coefficients which are infinite at that point, the grid must exclude the origin or specific pole conditions must be imposed. An efficient and accurate pseudo-spectral method has been here developed using collocation Chebyshev polynomials in the radial and axial directions and Fourier approximation in the azimuthal direction. To avoid the difficulty on the axis without prescribing any pole conditions, a new approach based on the work of Heinrichs [W. Heinrichs J. Comp. Phys. 199 (2004) 66-86] has been developed. The calculation domain is defined as (r, θ , z) ∈ [-1, 1] × [0, 2π] × [-1, 1] using an even number N of collocation points in the radial direction. Thus, r = 0 is not a collocation point. The clustering of collocation points around the rotation axis is also avoided due to the utilization of a Gauss-Lobatto distribution. The flow being indeed laminar close to the axis in most of the rotating flows. In the azimuthal direction, the overlap in the discretization is avoided by introducing a shift equal to π/2K (K the number of mesh points in θ -direction) for θ > π in the Fourier transform. The accuracy of the method was checked on the exact steady and unsteady analytical solutions and the capability of the method to simulate complex flows is illustrated considering the well documented case of the vortex breakdown phenomenon.

Introduction

The use of a cylindrical coordinate system introduces a numerical difficulty due to the so-called cylindrical coordinate singularity. Some studies have considered the special behaviour of flow equations at the origin of cylindrical coordinate systems and have proposed different techniques. Apparently, these approaches depend greatly on whether the spectral, finite-volume or finite-difference method is used but in fact there are some similarities between them.

The singularity is due to the presence of terms 1/r n (n = 1, 2) in the Navier-Stokes equations governing the flow (r the radial distance), although the flow field itself does not exhibit any singular behaviour on the axis. The treatment of the cylindrical coordinate singularity can be done adequately by means of several techniques taking into account the spatial discretization method and the flow characteristics.

In the context of spectral methods, Boyd and Yu [START_REF] Boyd | Comparing seven spectral methods for interpolation and for solving the Poisson equation in a disk: Zernike polynomials, Logan-Shepp ridge polynomials, Chebyshev-Fourier series, cylindrical Robert functions, Bessel-Fourier expansions, square-to-disk conformal mapping and radial basis functions[END_REF] listed several options of radial basis sets and radial grids for a finite cylindrical domain. These methods developed with the objective of capturing the solution near the coordinate singularity include methods that expand the solution in specially designed basis functions, such as spherical harmonics, parity-modified Fourier series, modified Robert functions and eigenfunctions of singular Sturm-Liouville problems among others.

Heinrichs [START_REF] Heinrichs | Spectral collocation schemes on the unit disc[END_REF] proposed spectral collocation schemes to solve the Poisson problem on the unit disc by discretizing the equations over the whole diameter and choosing an even number of collocation points in that direction. This method has been here extended to solve the full Navier-Stokes equations. An angular shift is introduced in the Fourier transform that avoid pole and parity conditions usually required. It presents some advantages, among others: it keeps the spectral convergence that reduces the number of grid points with respect to lower-order numerical methods. The grid-points distribution densifies the mesh only near the boundaries that makes the algorithm well-suited to simulate rotating cavity flows where thin layers develop along the walls. The reader can refer to Bernardi et al. [START_REF] Bernardi | Spectral Methods for Axisymmetric Domains[END_REF] for a complete overview of the different methods treating this singularity using spectral methods.

Mathematical modelling and numerical approximation

Cylindrical cavities of radius R and height H will be considered in the following. The motion of an incompressible fluid within such a cavity is governed by the 3D Navier-Stokes equations written in the velocity-pressure formulation and the continuity equation. The components of the velocity vector V are denoted u, v and w in the cylindrical coordinate system (r, θ , z) and P is the pressure. The momentum and continuity equations, in dimensionless form, can be written in the domain D of border ∂ D as:

∂ V ∂t + (V.∇)V = -∇P + 1 Re ∆ V + F in D, ∇.V = 0 in D ( 1
)
where t is the time and F represents a given body force. Appropriate boundary conditions of Dirichlet type for the velocity vector can be written as V = W on ∂ D, with ∆ the vectorial Laplacian operator. The resolution of Equations ( 1) together with appropriate boundary conditions leads to a Stokes problem coupling velocity and pressure. The collocation method requires boundary conditions in the physical space that can lead to an indetermination at r = 0 when the domain is discretized over the radius for r ∈ [0, R] as in [START_REF] Serre | A three-dimensional pseudospectral method for rotating flows in a cylinder[END_REF]. Here, the computational domain is discretized in the radial direction over the whole diameter between -R and R with an even number N of Gauss-Lobatto collocation points, r i = cos(iπ/N), 0 ≤ i ≤ N. Thus, no collocation point is located at r = 0 and no pole conditions are required.

For a given even number of points K, a shift equal to π/ (2K) for θ > π is introduced in the angular discretization over [0, 2π], as formerly proposed by Heinrichs [START_REF] Heinrichs | Spectral collocation schemes on the unit disc[END_REF] to solve two-dimensional Poisson equations on an unit disc. According to the angular shift, the following discrete angles θ j are defined by:

θ j = ( j -1) 2π K if j = 1, ..., K 2 , π + π K + j -K 2 + 1 2π K if j = K 2 + 1, ..., K. (2)
Figure 1(a) represents an example for K = 8. Such discretization clearly leads to a better resolution in the azimuthal direction with twice the amount of points although the number of Fourier modes remains equal to K/2. An example of the 3D grid is shown on Figure 1(b). The distribution of Gauss-Lobatto collocation points over the diameter tightens the points within the boundary layers along the walls and avoids unnecessary clustering of points around the axis.

According to this angular discretization, the following real Fourier basis can be introduced:

Φ k (θ ) =      sin (kθ ) if k = 1, ..., K 2 -1, cos (kθ ) + sin (kθ ) if k = K 2 , cos k -K 2 + 1 θ if k = K 2 + 1, ..., K. (3)
that corresponds to the K × K transformation matrix in the spectral space P θ = (Φ k (θ j )), with j, k = 1, ..., K. According to this basis, the first derivation of any variable is given by the matrix product

D θ P -1 θ Ψ , where D θ is the K × K full matrix D θ = Φ k (θ j ) 1≤ j,k≤K
. For 1 ≤ j ≤ K. Φ k (θ j ) can be written as:

Φ k (θ j ) =      k cos kθ j , k = 1, ..., K 2 -1, -k sin kθ j + k cos kθ j , k = K 2 , -k -K 2 + 1 sin k -K 2 + 1 θ j , k = K 2 + 1, ..., K.
As referred by Peyret [START_REF] Peyret | Spectral methods for incompressible viscous flow[END_REF], numerical calculations have shown that the presence of the term at k = K/2 may lead to instabilities in time-dependent problems mainly because the first derivative cannot be represented in this basis. Therefore the mode K/2 must be filtered. The second-order derivative in Fourier space can be written as:

D2 θ = -diag k 2 if k = 1, ..., K 2 , -diag k -K 2 + 1 2 if k = K 2 + 1, ..., K. (4)
Fast Fourier Transform algorithm (FFT) based on Equation [START_REF] Bernardi | Spectral Methods for Axisymmetric Domains[END_REF] is not yet available. Consequently, the Fourier transform in the azimuthal direction is performed here by matrix products using the subroutine DGEMM of the Blas (Basic Linear Algebra Subprograms) library.

The velocity-pressure coupling is solved using the 3D projection algorithm proposed by Raspo et al. [START_REF] Raspo | A spectral projection method for the simulation of complex three-dimensional rotating flows[END_REF]. It is an improved version of the Goda's algorithm, which introduces a preliminary step for the pressure in order to allow a temporal evolution of the normal pressure gradient at the boundaries. It is based on a semi-implicit second-order scheme that combines an implicit backward Euler scheme for the diffusive terms and an explicit Adams-Bashforth extrapolation for the convective nonlinear terms.

All variables Ψ are spatially discretized using a collocation-Chebyshev method for both non-homogeneous directions (r, z) and a Fourier-Garlekin approximation in the azimuthal direction θ . Numerical solution is thus sought under the form of the following truncated series:

Ψ NKM (r, θ , z) = K ∑ k=1 N ∑ n=0 M ∑ m=0 Ψnkm T n (r) T m (z) Φ k (θ ) (5)
where T n and T m are Chebyshev polynomials of highest degrees N and M respectively. (r, z) denote the dimensionless space variables (r, z) mapped into the square

[-1, 1] × [-1, 1]
, a requisite for the use of Chebyshev polynomials. N and M correspond to the number of Gauss-Lobatto collocation points defined as the extrema of the Chebyshev polynomials of highest degrees N and M. Φ k (θ ) are the trigonometric basis functions defined in Equation ( 3) and K defines the cut-off frequency of the Fourier series that corresponds to K/2 modes. The explicit non-linear terms are calculated using the pseudospectral techniques described in [START_REF] Peyret | Spectral methods for incompressible viscous flow[END_REF], where the derivatives in each space direction are calculated in the spectral space and the products are performed in the physical one. The derivatives of implicit diffusion terms are all performed in the physical space using matrix products.

The coupling between the first two components of the vectorial cylindrical Laplacian usually requires a splitting of the operator in an implicit and an explicit part that leads to time-step limitations. Therefore, the variable transformation u + = u+iv and u -= uiv (i 2 = -1) is introduced [START_REF] Orszag | Secondary instability of wall-bounded shear flows[END_REF], that makes the Laplacian diagonal:

(∆V ) + = ∇ 2 - 1 r 2 + 2i r 2 ∂ ∂ θ u + (6)
(∆V ) -= ∇ 2 - 1 r 2 - 2i r 2 ∂ ∂ θ u - (7)
(∆V ) z = ∇ 2 w (8)
Such transformation requires the reorganization in the complex plane of the Fourier coefficients matrix Ψ as follows:

Ψ =    Ψ j, K 2 +k , k = 1, Ψ j, K 2 +k + i Ψ j,k-1 , k = 2, ..., K/2 (9)
for 1 ≤ j ≤ N × M. The Fourier mode K/2 has been set to zero [START_REF] Peyret | Spectral methods for incompressible viscous flow[END_REF] as previously mentioned. In accordance with the variable transformation previously referred, the spectral matrices u + and u -can be written as: u + = u + i v and u -= ui v. Thus,

u + =    u j, K 2 +k + i v j, K 2 +k , k = 1, u j, K 2 +k -v j,k-1 + i u j,k-1 + v j, K 2 +k , k = 2, ..., K/2 (10)
u -=    u j, K 2 +k + i v j, K 2 +k , k = 1, u j, K 2 +k + v j,k-1 + i u j,k-1 -v j, K 2 +k , k = 2, ..., K/2 (11)
Finally, the coefficients u j, K 2 +k , u j,k-1 , v j, K 2 +k and v j,k-1 write for 1 ≤ k ≤ K/2:

   u j, K 2 +k = Re u + + u - 2 , u j,k-1 = Im u + + u - 2 v j, K 2 +k = Im u + -u - 2 , v j,k-1 = Re u --u + 2 (12)
The time-discretization and the Fourier-Galerkin approximation described above lead to a set of elliptic equations for each Fourier wave to be solved in a 2D domain depending on the two non-periodic directions (r, z). Successive 2D Helmholtz and Poisson equations can be written in the matrix form as:

AΨ +Ψ B T = S (13) with A = D (2) r -1 r D (1)
rσ (Ψ ) I and B = D

(2) z

. It is noted that only B does not depend on the harmonic k. The Dirichlet boundary conditions are taken into account in S, while the Neumann boundary conditions appear in A or B according to the direction. These 2D algebraic linear systems are solved directly using a full diagonalization technique for each Fourier mode [START_REF] Raspo | A spectral projection method for the simulation of complex three-dimensional rotating flows[END_REF].

Accuracy and validation of the method

The method is validated for flows in shrouded cylindrical rotor-stator cavities of aspect ratio G = 2R/H (R the disc radius and H the interdisc spacing). The fluid of kinematic viscosity ν is driven by the rotation of the top disc at constant speed Ω , the other walls being at rest. The characteristics scales for space, time and velocity are H, Ω -1 and Ω R respectively. The main flow is governed also by the rotational Reynolds number defined as Re = Ω R 2 /ν. The accuracy of the method is checked by using the steady and unsteady divergence-free velocity field V e = (u e , v e , w e ) for a cavity of aspect ratio G = 10 and at a Reynolds number Re = 250. The velocity field V e satisfies homogeneous Dirichlet boundary conditions. The space and time accuracy is evaluated by computing the discrete error at the inner collocation points, fully described in Raspo et al. [START_REF] Raspo | A spectral projection method for the simulation of complex three-dimensional rotating flows[END_REF].

For the exact steady solution, the computations are initialized using the disturbed velocity field V 0 = V e + 1/(200π 2 ) * V e . The time step is fixed at δt = 5 × 10 -3 . The steady solution is assumed to be obtained when the residuals

Res(Ψ ) = max i,k, j |Ψ n NKM (r i , θ k , z j ) -Ψ n-1 NKM (r i , θ k , z j )|/δt are below 10 -10
. The spectral convergence of the method is obtained for all variables. The error reported in Figure 2(a), decays exponentially to the zero machine for N = M = K = 28 then slightly increases due to the round-off error as usually observed [START_REF] Raspo | A spectral projection method for the simulation of complex three-dimensional rotating flows[END_REF].

The exact time-dependent solution is time periodic with a period T = 0.125 and the computations are initialized with V 0 = V e (t = 0). The mesh is fixed at

Results

The capability of the method to simulate complex flows is illustrated considering the well documented case of the vortex breakdown phenomenon, which may appear in small aspect ratio cavities depending on the Reynolds number. Such flow produces physical phenomena of weak intensity with respect to the base flow in the vicinity of the axis, that still today challenge all numerical solvers. Velocity boundary conditions are no-slip on all walls and v = r at the rotating disc. The singularity at the junction between the stationary cylinder and the rotor is appropriately regularized by using an exponential smoothing of the v-profile as proposed in [START_REF] Serre | A three-dimensional pseudospectral method for rotating flows in a cylinder[END_REF].

The confined vortex flow can undergo breakdown for certain combinations of (G, Re). The vortex breakdown phenomenon occurs as an abrupt change in the structure of the vortex core around the cylinder axis and typically develops downstream into a stagnation zone in the flow. Escudier [START_REF] Escudier | Observations of the flow produced in a cylindrical container by a rotating end wall[END_REF] showed experimentally that several recirculation bubbles may exist on the main vortex axis. The accuracy of the algorithm to simulate secondary vortices around r = 0 is shown on simulations covering the 3D vortex breakdown. The mesh size is fixed at 126 × 48 × 65 in the radial, azimuthal and axial directions respectively and the time step is equal to 10 -2 .

At Re = 6500, the flow is fully three-dimensional with some asymmetries in the breakdown structure (Fig. 3). The vortex breakdown is related here to the 'S-shaped' regime first noted by Escudier [START_REF] Escudier | Observations of the flow produced in a cylindrical container by a rotating end wall[END_REF]. The solution is similar to the one obtained by DNS [START_REF] Serre | Three-dimensional swirling flow with a precessing vortex breakdown in a rotor-stator cylinder[END_REF] and characterized by a precessing structure around the axis in the direction of the main flow. At this Reynolds number, the temporal regime is chaotic due to the interaction of boundary layer instability patterns with the precessing bubble. The precession motion of the structure is still not fully understood and it is possibly related to a symmetry breaking of the structure with respect to the m = 1 mode as suggested by Marques and Lopez [START_REF] Marques | Precessing vortex breakdown mode in an enclosed cylinder flow[END_REF]. Some simulations have been here performed at Re = 5500 together with a nonlinear normal form analysis based on linear global modes by P. Meliga (private communication). Figure 4 presents the kinetic energy of the mode m = 1 for several Reynolds number Re, showing the consistence of the results with the theoretical analysis of P. Meliga (private communication). Snapshots of w = 0 isosurface and isocontours of w show on Figure 5 the full and filtered solutions. In agreement with the theoretical analysis of Meliga, the results show that 3D modes of centrifugal instability within the vertical boundary layers are related to azimuthal wave numbers strictly larger than one since by keeping the modes m = 0 and m = 1 only, the main flow becomes almost axisymmetric (Fig. 5(a)-5(c)). In contrast, the vortex breakdown topology does not change and the precession seems related to a m = 1 mode of instability as confirmed by the new filtered solution on Figures 5(g)-5(i). Moreover, transition of the solution to the m = 1 mode as well as its amplitude match well with the theoretical predictions.

Conclusions

A spectrally convergent method has been proposed to solve the three-dimensional Navier-Stokes equations in cylindrical coordinates without any requirement of polar or parity condition, as in most of the methods proposed in the literature. To perform that, the equations have been discretized over the whole cavity diameter using a Gauss-Lobatto collocation point distribution that avoids any clustering of points in the stable region around the axis. In addition, an angular shift has been introduced in the azimuthal discretization that allows a more efficient discretization of the physical space in this direction. Numerical solutions for the vortex breakdown phenomenon match previous results available in the literature.

Fig. 1

 1 Fig. 1 Grid structure: (a) Discrete angle θ j distribution defined by Equation (2) in a (θ , z) plane showing the angular shift. Black and dashed diameters correspond to discretization indices j = 1, ..., K 2 and j = K 2 + 1, ..., K, respectively. Grey circles show the radial mesh distribution; (b) Example of 3D mesh discretization in a thin interdisc cavity.

Fig. 2

 2 Fig. 2 Space and time accuracy of the method for exact solutions. Evolution of the error at the inner collocation points Er for the radial velocity component u and for the pressure P: (a) versus the polynomial degrees N = M = K, (b) versus the time step δt.

Fig. 3

 3 Fig. 3 Unsteady and three-dimensional 'S-shaped' vortex breakdown for (G = 0.5, Re = 6500). Snapshots of w = 0 isosurface and isolines of the axial velocity component w at different times.

Fig. 4

 4 Fig. 4 Kinetic energy of the mode m = 1, E (m=1) , versus the Reynolds number Re for G = 0.5.

1 Fig. 5

 15 Fig. 5 Role of the m = 1 mode in the precession. Snapshots of w = 0 isosurface and isocontours of w in a (r, z) plane for (G = 0.5, Re = 5500) at different times t = 22Ω -1 , t = 30Ω -1 and t = 50Ω -1 . (a,b,c) full solution, (d,e,f) filtered solution for all modes m > 1 and (g,h,i) solution where only the mode m = 1 is filtered.


S Poncet

POLYGONS ON A ROTATING THIN-WATER LAYER: A COMBINED EX-PERIMENTAL AND NUMERICAL APPROACH

Keywords: rotating flow, free surface, shear-layer instability, DNS, flow visualizations

The flow of a thin fluid layer over a rotating disk is studied by direct numerical simulation and flow visualizations. The primary bifurcation appears as spectacular sharp-cornered polygonal patterns located along the outer casing. The stability diagram is established experimentally and confirmed numerically with a strong hysteresis between the spin-up and spin-down phases. These structures are induced by the appearance of a secondary flow. They rotate at a frequency, which slightly increases with the number of vortices. For larger Reynolds numbers, the polygons may coexist with a crossflow instability appearing as positive spirals.

Introduction

The present work is concerned with the stability of a water layer over a rotating disk, leading to steady modes of a shear-layer instability. Shear layers in rapidly rotating systems are of primary importance in geophysics, as they can be observed in oceans, hurricanes or planetary atmospheres but also in computer hard drives. The stability of a fluid layer over a rotating disk is governed by many parameters whose the aspect ratio of the cavity G = e/R o (e the water depth at rest and R o the disk radius). For G > 0.15, the polygons appear at high rotation rates such that a dry central zone may occur close to the rotation axis and the Froude number F r = Ω 2 R o /g (Ω the disk rotation rate) pilots the instability mechanisms. The surface flow may be divided into two regions: an inner center in solid body rotation and an outer annulus, where a weak secondary flow occurs [START_REF] Vatistas | Experimental confirmation of Kelvin's equilibria[END_REF][START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF]. For G < 0.15, the transition mechanisms are quite different [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF][START_REF] Kahouadji | Analyse de stabilité linéaire découlements tournants en présence de surface libre[END_REF]. The experimental data of Poncet and Chauve [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] are here revisited and underpinned by new simulations in this case.

Experimental and numerical methods

The cavity consists of a smooth disk (radius R o = 140) mm rotating at Ω enclosed by a cylindrical shroud (radius R o + j, j = 0.85 mm). A central hub (radius R i ) may be attached to the rotor. The cavity is filled up by water (kinematic viscosity ν = 10 -6 m 2 /s), seeded with reflective particles of kalliroscope. The water depth at rest is denoted e. The flow is mainly controlled by three parameters: the rotational Reynolds number Re = ΩR 2 o /ν, the aspect ratio G = e/R o of the cavity and its radius ratio η = R i /R o . Images (768 × 576 pixels) are taken at a video frequency of 25 Hz. Between two observations, Ω is increased (spin-up) or decreased (spin-down) by step of 1 rpm and one respects the Ekman time scale τ E = √ Re/Ω [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]. The numerical approach is based on a pseudospectral technique using a collocation-Chebyshev method in the radial and axial directions and a Galerkin-Fourier method in the azimuthal direction. This approximation is applied at the Gauss-Lobatto collocation points in the radial and axial directions. In the azimuthal direction, a uniform distribution is considered. The time scheme is semi-implicit and second-order accurate. The solution method is based on an efficient projection scheme to solve the velocity-pressure coupling. Finally, for each Fourier mode, a full diagonalization technique is used and yields simple matrix products for the solution of successive two-dimensional uncoupled Helmholtz and Poisson equations at each time step. The free surface is assumed to be flat and undeformable as the Froude number F r remains small (here max(Fr)= 0.027 << 1).

Results

Modes m = 8 to 3 have been observed for different combinations of Re and G (Fig. 1). The influence of the radius ratio η is discussed in detail in [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]. The flow is divided into two regions: the inner region in solid body rotation and the outer region where the vortices occur, separated by a polygonal boundary. Figure 2a presents the marginal stability diagram of the first polygonal mode in the plane (Re, G) for η = 0 during spin-up. The present experimental and numerical results compare quite well with the stability analysis of Kahouadji [START_REF] Kahouadji | Analyse de stabilité linéaire découlements tournants en présence de surface libre[END_REF] for G ≥ 0.07. When G decreases, Re and m increase. Kahouadji [4] developed also a curvilinear formulation to take into account the deformation of the free surface and tried to explain the discrepancies between the flow visualizations and the stability analysis for G ≤ 0.06. But it did not improve the present results for a flat free surface. The secondary bifurcations are shown in Figure 2b for η = 0. Strong hysteresis cycles are observed as illustrated in Figure 2c. The same kind of sequence is observed for all values of G: N polygonal patterns are observed during spin-up and N -1 during spin-down. From DNS, the m = 5 instability threshold is Re = 18000 during spin-up, while this mode remains stable up to Re = 15000 during spin-down, highlighting the subcritical nature of this instability. Figure 3 present the flow structure of the pentagon. It confirms that the flow is in solid body rotation in an inner region and that a secondary flow developed at the periphery. The fluid moves radially outwards along the rotor and inwards along the free surface. The map of the axial velocity clearly shows this recirculation zone at the origin of the mode m = 5. The polygonal structures do not rotate at the same rotation speed as the rotor. One can define the frequency ratio F = Ω f /Ω between the polygon and disk rotation rates. Figure 4 presents its variation with the Reynolds number for all aspect ratios and all modes. As observed also by Vatistas et al. [START_REF] Vatistas | Experimental confirmation of Kelvin's equilibria[END_REF] for large aspect ratios, F is rather constant with Re but slightly increases with m. Spiral patterns coexisting with a pentagonal structure have been observed experimentally and numerically (Fig. 5). Figures 5b-d show that the spirals are characterized by a relatively high turbulence kinetic energy, and that they flow radially inward from the rotor to the free surface.The spirals move slower than the vortices at around 16.6% of the pentagon frequency. These are positive spiral patterns as they are rolled up in the same sense as the rotor. These spirals are very similar to the SRJ2 spirals of Poncet and Chauve [START_REF] Poncet | Crossflow instability in a rotor-stator cavity with axial inward throughflow[END_REF] in a rotor-stator cavity with throughflow as they correspond also to a crossflow instability due to the inflexion point in the axial profile of the mean radial velocity. Their main characteristics are summed up in Table 1. Their inclination angle remains small compared to the values obtained by [START_REF] Poncet | Crossflow instability in a rotor-stator cavity with axial inward throughflow[END_REF] (ε → 70 • ). The number of spiral arms 28 ≤ n ≤ 40 is also much lower: n → 90 in [START_REF] Poncet | Crossflow instability in a rotor-stator cavity with axial inward throughflow[END_REF]. A weak hysteresis is observed on the thresholds of the spiral instability. n = 29 spiral arms with a positive angle ε ≃ 25 • have been observed numerically (Fig. 5b-d

Conclusion

One has reported a shear-layer instability developed on a rotating disk with a free surface as a sharp-cornered polygonal pattern with m vortices. The stability diagram has been established highlighting large hysteresis cycles. Modes up to m = 8 have been observed with a particular emphasis for the pentagonal structure. This pattern produced by the appearance of a secondary flow and the destabilization of the Stewartson layer along the casing rotates at 73.5% of the rotor speed. It may coexist with a crossflow instability appearing a positive spirals located along the shroud. Further DNS calculations are now required to understand the discrepancies between the experiments and the linear stability analysis of Kahouadji [START_REF] Kahouadji | Analyse de stabilité linéaire découlements tournants en présence de surface libre[END_REF] for low aspect ratios G ≤ 0.06.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Flow visualizations of the modes m = 8 to 3 for η = 0. The disk rotates counterclockwise.

Figure 4 :

 4 Figure 4: Dimensionless azimuthal velocity F of the polygon patterns deduced from the flow visualizations (black symbols) and the DNS (blue symbols) for 0.0143 ≤ G ≤ 0.143.

Figure 5 :

 5 Figure 5: Crossflow instability: (a) experimental results for η = 0, Re = 24630 (spin-up) and G = 0.0714; DNS results (73×192×65 mesh points in the (r, θ, z) directions with δt = 5×10 -5 Ω -1 ) for η = 0.5, Re = 27000 and G = 0.0714 at z/e = 0.97: (b) turbulence kinetic energy k * , (c) instantaneous radial V * r and (d) axial V * z velocities.

Table 1 :

 1 Characteristics of the spirals obtained experimentally.


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

Sébastien Poncet

Stéphane Viazzo

Ébastien Poncet

St Éphane Viazzo

TAYLOR-COUETTE-POISEUILLE FLOWS : FROM RANS TO LES

published or not. The documents may come

INTRODUCTION

Since Taylor's first theoretical results obtained one century ago, researchers have dedicated obstinate efforts in modelling the rich dynamical behaviour exhibited by the flow moving between two concentric cylinders in differential rotation. In spite of the substantial literature on flows through rotating concentric annuli, little is known about the near-wall turbulent structures in such systems. This kind of flow is often encountered in turbomachinery. In addition, it provides also insights into the general problem of threedimensional turbulent boundary layers. The combined axial and rotational flow in an annulus was studied experimentally by [START_REF] Escudier | Concentric annular flow with centerbody rotation of a newtonian and a shearthinning liquid[END_REF]. Most numerical studies until now used turbulence modelling, which provided rather limited informations. It was shown that increasing the rotation rate of the inner tube amplified the turbulent kinetic energy, resulting in the enhancement of heat transfer [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]. The development of computational methods (including Direct and Large Eddy Simulations) has led to an increase in numerical studies of rotating flows, but little concern turbulent Taylor-Couette-Poiseuille flows. As example, [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. Compared to this latter, we focus our attention here on the narrow gap case, which models quite faithfully flows in the rotor-stator gap of electrical motors. The aim of the present work is twofold : (i) providing some reference LES data for this case, using an in-house code and (ii) questioning the capabilities of a hybrid RANS / LES method, as well as the underlying RANS model, in predicting this kind of flow.

FLOW CONFIGURATION AND PARAMETERS

The fluid is confined between two concentric cylinders of radii R 1 and R 2 and of height h (see Figure 1). The inner cylinder represents a rotor, thus rotating at a constant rate Ω, while the outer one is stationary. The cavity may be 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 )/ν,

NUMERICAL METHODS

Concerning LES reference calculations, the numerical method of the in-house code is based on the work of [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF], which devised a 2D compact fourth-order projection decomposition method in cartesian coordinates. It has been recently validated in the case of turbulent interdisk rotor-stator flows in cylindrical coordinates by [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 middle gap cavity [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF]. The time advancement is second order 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 Figure 1. Sketch of the narrow gap cavity 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 three-dimensional extension uses Fourier series. The set of problems is thus reduced to a 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 LES grid contains 65 × 144 × 130 points respectively in the radial, azimuthal and axial directions. The CFD code used for RANS and hybrid RANS/LES simulations is the open-source Code Saturne, developped 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 pressurevelocity 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 modelling 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 [START_REF] Manceau | Of the interpretation of DES as a hybrid RANS / temporal LES method[END_REF][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 developped on a rather phenomenological basis. First, the PITM was generalized by Fadai-Ghotbi et al. (2010b) 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. (2010b) 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 [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, 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 SFS energy or subfilter stresses τ i j SFS , damping the mod- elled energy 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 (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
)
HBD β 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 [START_REF] Manceau | Of the interpretation of DES as a hybrid RANS / temporal LES method[END_REF], 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 nearwall regions, to avoid the very fine resolution required by LES. In that aim, the RANS model proposed by [START_REF] Manceau | Elliptic blending model: A new near-wall Reynolds-stress turbulence closure[END_REF] and [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 , (6)
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 Fadai-Ghotbi et al. (2010a) for details) :

φ * i j -ε i j = α 3 (φ * i j -ε i j ) wall +(1-α 3 )(φ * i j -ε i j ) homogeneous .
(7) Under its RANS form, EB-RSM was successfully applied to several flows, such as 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], impinging jets [START_REF] Thielen | Predictions of flow and heat transfer in multiple impinging jets with an elliptic-blending second-moment closure[END_REF] or mixed and natural convection flows [START_REF] Dehoux | Modelling turbulent heat fluxes using the elliptic blending approach for natural convection[END_REF]. It was recently adapted by Fadai-Ghotbi et al. (2010a) 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 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 w is a constant, ∆ the grid step, S i j the filtered strainrate tensor, and S d i j the deviatoric part of the squared filtered velocity gradient tensor.

RESULTS AND DISCUSSION

Hereafter, some sample results are shown. Both rotation rates are represented on same graphs, in order to show the influence of N. Comparison between the reference LES, hybrid RANS/LES and RANS is proposed. Velocity profiles are first presented on Figure 2. Good agreement is observed, especially inside boundary layers, as expected, since near-wall treatment without empirical laws is a feature of EB-RSM. On the contrary, the classical turbulent Poiseuille profile observed for the axial velocity is not captured by the hybrid approach, which is not surprising (see e.g. Fadai-Ghotbi et al. (2010a)). At N = 2.24, the slope of azimuthal velocity V θ is also particularly badly predicted by RANS and hybrid RANS/LES in the free flow. Figures 3 and4 show the normal and shear components of the Reynolds stress tensor. As expected, R θ θ increases strongly with N, as well as R θ r and R θ z shear stresses intensities. Other components increase in a more marginal way, only due to redistribution. Another interesting feature, is that all components are far more symmetric than in the middle-gap (η = 0.5) studied by [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF], confirming that the flow field near the stator side is significantly influenced by the motion of the rotor, especially at those quite high values of N. Predictions by both RANS and hybrid RANS/LES approaches do not show clear hierarchy between them, but nevertheless RANS seems more efficient than hybrid RANS/LES for N = 2.24, and vice-versa for N = 4.47. Peak locations and especially intensities are not very well reproduced, but in a quite satisfactory way. This is the cause for the weak prediction of velocity profiles in the free flow, evoked above. resolved k res parts of the turbulent kinetic energy. For the same mesh, it seems that r is globally lower at N = 4.47 than at N = 2.24, which may explain why hybrid RANS/LES predictions are slightly better than RANS at N = 4.47, since the level of explicitly resolved energy is then bigger. Note that RANS mode r = 1 is enforced near walls, in order to allow EB-RSM to fully play its role. Nevertheless, the observed ratio does not obey it exactly. Indeed, the outer flow fluctuations influence the behaviour even inside the viscous sublayer.

Figures 6 and7 show isocontours of arbitrary Qcriterion values, seen from the rotor side. Obviously, the LES simulations show far more structures than hybrid RANS/LES ones. Nevertheless, an interesting behaviour is exhibited by both approaches : the streaks form a kind of net, because of the significant bidimensionality of the flow, and this net's angle changes with N, tending to place the main structures in the (r, θ ) plane.

Eventually, Table 2 summarizes the wall quantities and numerical resolutions in wall units, for all approaches. While good global agreement is found for Re τ at N = 2.24, there is more dispersion at N = 4.47, especially at the inner wall. Another interesting observation, is that Re τ is always equal at the inner and the outer walls, because the velocity gradient is the same at both walls.

SUMMARY

The present work is a first step in numerical simulation of Taylor-Couette-Poiseuille flows in narrow rotorstator cavities. Assuming that Large Eddy Simulation results can be considered as reference as this approach has been already validated in a middle gap cavity [START_REF] Oguic | Numerical simulations of a middle gap turbulent Taylor-Couette-Poiseuille flow[END_REF]

Figure 2 .

 2 Figure 2. Azimuthal and axial velocity profiles, for both rotation rates. Profiles for N=4.47 are shifted up, for clarity.

Figure 5 Figure 3 .

 53 Figure 5 compares observed and targetted values of the energy ratio r. "Targetted" means that r is evaluated by Eq.(3), while "observed" refers to post-processed r, in the fashion of Eq. (2), requiring to compute the modelled k m and the

Figure 4 .Figure 5 .

 45 Figure 4. Shear components of the Reynolds stress tensor for both rotation rates. Profiles for N=4.47 are shifted up, for clarity.

Figure 6 .

 6 Figure 6. Isocontours of Q-criterion (arbitrary value), for N=2.24. Up : hybrid RANS/LES, down : reference LES

Table 1 .

 1 ν being the fluid's kinematic viscosity. The present study considers two values of the ratio N = Re Ω /Re Q (see Table1), while other parameters remain unchanged. Geometrical and flow parameters for the two cases

	Case	A	B
	Γ	10	10
	η	8/9	8/9
	Re Q 3745 3745
	N	2.24 4.47

Table 2 .

 2 , it seems Wall resolution and mean flow parameters.

	HBD

ACKNOWLEDGEMENTS

This work was financially supported by Liebherr Aerospace (project "Entrefer moteur") and by the LABEX Figure 7. Isocontours of Q-criterion (arbitrary value), for N=4.47. Up : hybrid RANS/LES, down : reference LES MEC (project "Hydrex"), which are gratefully acknowledged.