Bikash Sahoo
email: bikashsahoo@nitrkl.ac.in

Sébastien Poncet
email: sebastien.poncet@univ-amu.fr

Blasius flow and heat transfer of a fourth grade fluid with slip

Keywords: Blasius flow, Partial slip, Fourth grade fluid, Broyden's method

This investigation deals with the effects of slip, magnetic field and non-Newtonian flow parameters on flow and heat transfer of an incompressible, electrically conducting fourth grade fluid past an infinite porous plate subject to suction and blowing, known as the Blasius flow problem. The heat transfer analysis has been carried out for two heating processes, namely, when the plate is assumed to be at higher temperature than the fluid, and in the second case, the plate is assumed to be insulated. The order of the resulting differential equations exceeds the available adherence boundary conditions due to the presence of the higher order material derivatives of the strain tensor in the constitutive relation for the fourth grade fluid. The system of highly non-linear differential equations is solved using a fourth order Runge-Kutta method along with a shooting method for moderate values of the flow and heat transfer parameters. The effective Broyden's technique has been adopted in order to improve the initial guesses and to satisfy the boundary conditions at infinity. The combined effects of the non-Newtonian flow parameters and slip on the momentum and thermal boundary layers have been discussed in details. An exceptional crossovers is obtained in the velocity profile in presence of slip. The fourth grade fluid parameter is found to increase the momentum boundary layer thickness, whereas the slip parameter substantially decreases it. Similarly, the non-Newtonian fourth grade fluid parameters and the slip have opposite effects on the thermal boundary layer thickness.

Introduction

Despite their applications in many industrial processes, the boundary layer flow of higher order fluid of differential type, namely the fourth grade fluid here, has drawn the attention of only few researchers during the last few decades. The main reason is that the constitutive equation of the fourth grade fluid is very complex involving a wide number of parameters. Moreover, the order of the differential equation corresponding to the momentum equation is exceptionally large, and exceeds the number of available boundary conditions, even in the simplest boundary layer flow past a flat plate. Renardy [START_REF] Renardy | On the domain space for constitutive laws in linear viscoelasticity[END_REF] has shown that such fluids are the last to have a stable state of rest. They have been considered in few other configurations such as plane Poiseuille [START_REF] Hayat | On flow of a fourth-grade flud with heat transfer[END_REF][START_REF] Islam | The optimal solution for the flow of a fourth-grade fluid with partial slip[END_REF] or Couette [START_REF] Shah | Couette and Poiseuille flows for fourth grade fluids using optimal homotopy asymptotic method[END_REF] flows or recently, in an asymmetric channel [START_REF] Mehmood | Heat transfer on peristaltic flow of fourth grade fluid in inclined asymmetric channel with partial slip[END_REF].

The unsteady and steady flows of a fourth grade fluid, subject to a transverse uniform magnetic field past a porous plate has been investigated by Hayat et al. [START_REF] Hayat | Flow of a fourth grade fluid[END_REF]. The boundary layer approximation resulted in a highly nonlinear differential equation of order six corresponding to the momentum equation with only two adherence boundary conditions. The finite difference method has been adopted to solve the sixth order nonlinear differential equation by augmenting four extra asymptotic boundary conditions. Hayat et al. [START_REF] Hayat | The unsteady flow of a fourth-grade fluid past a porous plate[END_REF] have studied the unsteady flow of a fourth grade fluid past a porous plate. They have successively reduced the order and number of independent variables of the resulting partial differential equations using Lie point symmetries. Subsequently, Sajid et al. [START_REF] Sajid | On the analytic solution of the steady flow of a fourth grade fluid[END_REF] have studied the steady flow of a fourth grade fluid past a plate with suction and injection velocities at the surface of the plate. The homotopy analysis method (HAM) is used to find the solution of the fifth order highly nonlinear differential equation with two available boundary conditions. Hayat et al. [START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF] have used the HAM method to investigate the influence of the heat transfer on the flow of a fourth grade fluid past a porous plate. The unsteady MHD flow of a fourth grade fluid due to an oscillating plate with suction and blowing has been investigated numerically by Wang and Wu [START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF]. The corresponding sixth order nonlinear partial differential equation has been solved by finite difference technique after augmenting four asymptotic boundary conditions. Recently, Marinca et al. [START_REF] Marinca | An optimal homotopy asymptotic method applied to the steady flow of a fourth-grade fluid past a porous plate[END_REF] have used optimal homotopy asymptotic method to solve the resulting highly nonlinear momentum equation with no-slip boundary condition.

The foremost reason which motivated for the present study is the nonlinearity and the order of the differential equation corresponding to the momentum equation due to the flow of a fourth grade fluid. The order of the resulting differential equation exceeds drastically from the available adherence boundary conditions due to the presence of the higher order material derivatives of the strain tensor in the constitutive relation for the fourth grade fluid. The literature survey reveals that in the aforementioned studies, respective authors [START_REF] Hayat | Flow of a fourth grade fluid[END_REF][START_REF] Hayat | The unsteady flow of a fourth-grade fluid past a porous plate[END_REF][START_REF] Sajid | On the analytic solution of the steady flow of a fourth grade fluid[END_REF][START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF][START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF] have augmented additional boundary conditions at infinity. However, these are relatively weak and physically plausible. Few other interesting results from the literature survey, which arouse a curiosity to consider the MHD flow and heat transfer of a fourth grade fluid with partial slip boundary condition are as follows. Hayat et al. [START_REF] Hayat | Flow of a fourth grade fluid[END_REF] have reported that in presence of the magnetic field, the boundary layer thickness decreases with an increase in the fourth grade fluid parameter, and the magnetic field has a substantial effect on the boundary layer thickness. The conclusions of the study of the steady flow of a fourth grade fluid past a porous plate by Sajid et al. [START_REF] Sajid | On the analytic solution of the steady flow of a fourth grade fluid[END_REF] reveal that there does not exist a solution in the case of injection and the boundary layer thickness is rapidly decreased with an increase in the fourth grade fluid parameter for all values of the suction parameter exceeding one. Moreover, the recent study [START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF] of the flow and heat transfer of a fourth grade fluid past a porous plate depicts that the thermal boundary layer thickness increases with an increase in the fourth grade fluid parameter. We can cite also the recent contribution of Islam et al. [START_REF] Islam | The optimal solution for the flow of a fourth-grade fluid with partial slip[END_REF], who investigated the steady flow of fourth-grade fluid in a infinite plane Poiseuille channel with slippage using the optimal HAM. They highlighted that this technique provides more satisfactory results than the homotopy perturbation technique.

Hayat et al. [START_REF] Hayat | Partial slip effects on the flow and heat transfer characteristics in a third grade fluid[END_REF] considered partial slip effects on the flow and heat transfer characteristics of a third grade fluid past a porous plate using the homotopy analysis method. Nadeem et al. [START_REF] Nadeem | Effects of partial slip on a fourth-grade fluid with variable viscosity: An analytic solution[END_REF] studied the effects of partial slip on a fourth-grade fluid with variable viscosity contained in a stationary cylinder. Analytic solutions of velocity and temperature have been developed and discussed for the Reynolds and Vogels models. The present paper is the extension of the recent work of Sahoo [START_REF] Sahoo | Flow and heat transfer of an electrically conducting third grade fluid past an infinite plate with partial slip[END_REF], who investigated the flow and heat transfer of an electrically conducting third grade fluid past an infinite plate with partial slip. It deals in the present case with the numerical study of the slip flow and heat transfer of an incompressible, electrically conducting fourth grade fluid past an infinite porous plate subject to suction and blowing. The heat transfer analysis has been carried out for two heating processes, namely, when the plate is assumed to be at higher temperature than the fluid, and in the second case, the plate is assumed to be insulated. To the best of our knowledge, no attention has been given to the effects of partial slip on the MHD boundary layer flow and heat transfer of a fourth grade fluid over a flat plate. Even, the hydrodynamic boundary layer flow and heat transfer of a fourth grade fluid with slip boundary condition is not discussed so far. However, some non-Newtonian fluids like polymer melts exhibit macroscopic wall slip, which makes the no-slip condition fully inappropriate. It is noteworthy to point out that the quantitative behavior of higher grade fluids is strongly dependent on the choice of the material parameters. At the present time, there appears to be no experimental data determining these parameters, particularly the material parameters in the theory of fourth grade fluids. The present investigation is a step towards understanding the effects of the slip parameter on the flow and heat transfer of a fourth grade fluid. This paper is divided as follows: in Section 2, the problem is mathematically formulated to get a system of ordinary differential equations with appropriate boundary conditions and is solved using an effective second order numerical scheme in Section 3. The effects of the flow parameters (partial slip, heat transfer, non newtonian parameters, suction/injection velocity) on a Blasius flow of a fourth grade fluid are discussed in details in Section 4 before some concluding remarks in Section 5.

Formulation of the problem

We consider the steady two-dimensional flow of an incompressible fourth grade fluid past an infinite porous plate (infinite in the x-direction), coinciding with the plane y = 0. The flow far away from the plate is uniform, and is in the direction x parallel to the plate. A transverse uniform magnetic field (0, B 0 , 0) (Fig. 1) has been applied at the surface of the plate and the surface admits partial slip. Thus, we seek velocity and temperature fields under the form:

V = u(y)i + v(y)j (1) θ = θ(y) (2)
where u and v are the components of the velocity in the x and y directions respectively and θ denotes the temperature.

Flow analysis

Neglecting the body force, the equations of motion and continuity for an electrically conducting, incompressible fluid write:

ρ dV dt = ∇ • T + J × B (3)
∇ • V = 0 (4)
where d/dt denotes the temporal material derivative, ρ is the density of the fluid assumed to be constant, T is the Cauchy stress for an incompressible homogeneous fluid of grade four, and J is the electric current. The term due to the Lorentz force (J × B) is simplified by making some assumptions valid for small magnetic Reynolds numbers and no displacement current [START_REF] Shercliff | A text book of magnetohydrodynamics[END_REF]: the physical quantities are kept constant, the magnetic field is perpendicular to the velocity field, the induced magnetic field remains small compared to the applied magnetic one and the electrical field is assumed to be zero.

Under the consideration of the flow, it follows from ( 1) and ( 4) that:

v = -V 0 (5)
where V 0 > 0 is the suction velocity and V 0 < 0 is the injection velocity. V 0 is considered as being constant.

The Cauchy stress T for an incompressible fourth grade fluid has the constitutive equation (see in [START_REF] Mehmood | Heat transfer on peristaltic flow of fourth grade fluid in inclined asymmetric channel with partial slip[END_REF] for example):

T = -pI + µA 1 + α 1 A 2 + α 2 A 1 2 + β 1 A 3 + β 2 (A 2 A 1 + A 1 A 2 ) + β 3 (trA 1 2 )A 1 + γ 1 A 4 + γ 2 (A 3 A 1 + A 1 A 3 ) + γ 3 A 2 2 + γ 4 (A 2 A 1 2 + A 1 2 A 2 ) + γ 5 (trA 2 )A 2 + γ 6 (trA 2 )A 1 2 + [γ 7 (trA 3 ) + γ 8 tr(A 2 A 1 )]A 1 ( 6
)
where I is the identity tensor. The coefficient µ is the coefficient of viscosity (kg/m/s). α 1 and α 2 are the second order material constants (kg/m), β 1 , . . . , β 3 the third order material constants (kg.s/m) and γ 1 , . . . , γ 8 the fourth order material constants (kg.s 2 /m). The first to fourth order Rivlin-Ericksen tensors A i (i = 1, . . . , 4) are defined by the recursive relation:

A 0 = I (7) A n = dA n-1 dt + A n-1 ∇V + (∇V) T A n-1 , n ≥ 1 (8)
where (∇V) T is the transpose of ∇V.

Now substituting (1) and ( 5) in the constitutive equation ( 6) for the fourth grade fluid, we obtain the following physical components of the Cauchy stress tensor T:

T xx = -p+α 2 du dy 2 -2β 2 V 0 du dy d 2 u dy 2 +2γ 2 V 2 0 du dy d 3 u dy 3 +γ 3 V 2 0 d 2 u dy 2 2 +2γ 6 du dy 4 ( 9
)
T xy = µ du dy -α 1 V 0 d 2 u dy 2 + β 1 V 2 0 d 3 u dy 3 + 2(β 1 + β 3 ) du dy 3 -γ 1 V 3 0 d 4 u dy 4 -(6γ 2 + 2γ 3 + 2γ 4 + 2γ 5 + 6γ 7 + 2γ 8 )V 0 du dy 2 d 2 u dy 2 (10)
T yy = -p + (2α 1 + α 2 ) du dy 2 -2V 0 (3β 1 + β 2 ) du dy d 2 u dy 2 + 6γ 1 V 2 0 d dy du dy d 2 u dy 2 + 2(γ 1 + γ 2 )V 2 0 d 3 u dy 3 du dy + γ 3 V 2 0 d 2 u dy 2 2 + 2(2γ 3 + 2γ 4 + 2γ 5 + γ 6 ) du dy 4 (11)
T zz = -p (12) T xz = T zy = 0 (13) T xy = T yx , T xz = T zx , T yz = T zy (14)
Now, the three scalar components of the momentum equation (3), in view of Eqn. ( 5), become:

-
)
where σ is the electrical conductivity.

Inserting the above stress components in Eqns. ( 15)-( 17), we obtain:

µ d 2 u dy 2 + ρV 0 du dy -α 1 V 0 d 3 u dy 3 + β 1 V 2 0 d 4 u dy 4 + 6(β 2 + β 3 ) du dy 2 d 2 u dy 2 -γ 1 V 3 0 d 5 u dy 5 -2(3γ 2 + γ 3 + γ 4 + γ 5 + 3γ 7 + γ 8 )V 0 d dy du dy 2 d 2 u dy 2 -σB 2 0 u = ∂ p ∂x ( 18
)
σB 2 0 V 0 = ∂ p ∂y (19) 0 = ∂ p ∂z ( 20
)
where the modified pressure p is is given by:

p = -T yy =p -(2α 1 + α 2 ) du dy 2 + 2V 0 (3β 1 + β 2 ) du dy d 2 u dy 2 -6γ 1 V 2 0 d dy du dy d 2 u dy 2 -2(γ 1 + γ 2 )V 2 0 d 3 u dy 3 du dy -γ 3 V 2 0 d 2 u dy 2 2 -2(2γ 3 + 2γ 4 + 2γ 5 + γ 6 ) du dy 4 (21)
Eqn. (20) indicates that p is not a function of z. On eliminating the pressure gradient from ( 18) and ( 19) by cross differentiation, we finally obtain:

-ρV 0 d 2 u dy 2 =µ d 3 u dy 3 -σB 2 0 du dy -α 1 V 0 d 4 u dy 4 + β 1 V 2 0 d 5 u dy 5 + 6(β 2 + β 3 ) d dy du dy 2 d 2 u dy 2 -γ 1 V 3 0 d 6 u dy 6 -2(3γ 2 + γ 3 + γ 4 + γ 5 + 3γ 7 + γ 8 )V 0 d 2 dy 2 du dy 2 d 2 u dy 2 (22)
Equation ( 22) is the stationary form of the one used by Wang and Wu [START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF].

The appropriate partial slip boundary conditions are:

u(0) = λ 1 T xy y=0 u(y) → U ∞ as y → ∞ (23)
where λ 1 is the slip factor and U ∞ is the mainstream velocity.

By introducing the following non-dimensional quantities, with ν the kinematic viscosity of the fluid,

u = u U ∞ , y = U ∞ y ν , V 0 = V 0 U ∞ , M n = νσB 2 0 ρU 2 ∞ , α 1 = α 1 U 2 ∞ ρν 2 , β 1 = β 1 U 4 ∞ ρν 3 , γ 1 = γ 1 U 6 ∞ ρν 4 , β = 6(β 2 + β 3 )U 4 ∞ ρν 3 , γ = 2(3γ 2 + γ 3 + γ 4 + γ 5 + 3γ 7 + γ 8 ) U 6 ∞ ρν 4 (24)
the boundary value problem (22) takes the form:

M n du dy -V 0 d 2 u dy 2 - d 3 u dy 3 + α 1 V 0 d 4 u dy 4 -β 1 V 2 0 d 5 u dy 5 + γ 1 V 3 0 d 6 u dy 6 =
u = λ du dy -V 0 α 1 d 2 u dy 2 + V 2 0 β 1 d 3 u dy 3 + β 3 du dy 3 -V 3 0 γ 1 d 4 u dy 4 -V 0 γ du dy 2 d 2 u dy 2 at y = 0 u → 1 as y → ∞ ( 26
)
where λ is the non-dimensional slip factor indicating the relative importance of the partial slip, and for simplicity, the bars of the non-dimensional quantities have been omitted. M n is the magnetic parameter and α 1 , β 1 , γ 1 , β and γ are the non-Newtonian flow parameters.

Heat transfer analysis

We now turn our attention to the energy equation subject to two different sets of boundary conditions. In the first instance, we assume that the plate is at a temperature higher than the fluid. For the second case, we assume that the plate is insulated. Neglecting the radiant heating, the energy equation can be written as:

ρc p dθ dt = κ∇ 2 θ + T.∇V (27)
where c p is the mass specific heat (J/K/kg) and κ the thermal conductivity (kg.m/s 3 /K). Since the velocity and temperature depend only on y, using Eqns. ( 1) and (2) in the above energy equation ( 27), we get:

-ρc p V 0 dθ dy = κ d 2 θ dy 2 + T xy du dy (28)
The boundary conditions depend on the following two heating processes.

Constant surface temperature

In this case, the appropriate boundary conditions are:

θ(0) = θ w , θ(y) → θ ∞ as y → ∞ (29)
where θ w is the wall temperature and θ ∞ is the temperature of the ambient fluid.

Using the non-dimensional temperature θ = θ-θ∞ θw-θ∞ , the non-dimensional quantities (24) and the expression for T xy from Eqn. [START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF] in the energy equation (28), we obtain the following non-dimensional highly non-linear differential equation:

d 2 θ dy 2 + P r V 0 dθ dy + P r E c du dy 2 -α 1 V 0 d 2 u dy 2 du dy + β 1 V 2 0 d 3 u dy 3 du dy + β du dy 4 -γV 0 du dy 3 d 2 u dy 2 -γ 1 V 3 0 d 4 u dy 4 du dy = 0 ( 30
)
where P r = µc p κ is the Prandtl number and

E c = U 2 ∞ cp(Tw-T∞)
is the Eckert number and bars have been suppressed.

The boundary conditions for θ are:

θ(0) = 1, θ → 0 as y → ∞ (31)

Insulated surface

In this case, the relevant boundary conditions are:

dθ dy y=0 = 0, θ(y) → θ ∞ as y → ∞ (32)
Again, using the non-dimensional temperature g = θ-θ∞ θw-θ∞ , the non-dimensional quantities (24) and the expression for T xy from Eqn. [START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF] in the energy equation (28), we obtain the following highly non-linear differential equation:

d 2 g dy 2 + P r V 0 dg dy + P r E c du dy 2 -α 1 V 0 d 2 u dy 2 du dy + β 1 V 2 0 d 3 u dy 3 du dy + β du dy 4 -γV 0 du dy 3 d 2 u dy 2 -γ 1 V 3 0 d 4 u dy 4 du dy = 0 (33)
which is same as Eqn.(30) and all the other parameters are the same as before.

Bars have been omitted.

The boundary conditions for g are:

dg dy y=0 = 0, g → 0 as y → ∞ (34)
3 Numerical solution of the problem

The system of highly non-linear differential equations ( 25) and (30) subject to boundary conditions ( 26) and (31), and the system of equations ( 25) and (33) subject to boundary conditions ( 26) and (34) are solved using a shooting method along with the fourth order Runge-Kutta method for different moderate values of the flow and heat transfer parameters. The effective Broyden's technique has been adopted in order to improve the initial guesses and to satisfy the boundary conditions at infinity. Moreover, it eases the computations by avoiding to compute the Jacobian matrix in each integration step. In fact, the above system of differential equations are highly non-linear, whose order is much more higher than the available known boundary conditions, which debars the use of the finite difference method without augmenting asymptotic boundary conditions at infinity [START_REF] Hayat | Flow of a fourth grade fluid[END_REF][START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF]. It is noteworthy to mention that the shooting technique in calculating process does not require the entire starting boundary conditions, instead it only needs to match with the other end of the boundary condition to add amendment to its starting boundary conditions. Thus, it can process higher order differential equation problem that lacks of sufficient boundary conditions. Despite of few shortcomings of the shooting method, we found the shooting method quite suitable and aesthetically pleasing for the present investigation.

The semi-infinite integration domain y ∈ [0, ∞) is replaced by a finite domain y ∈ [0, y ∞ ). In fact, y ∞ usually depends upon the physical parameters of the problem and its value needs to be adjusted as the values of the parameters change. In practice, y ∞ should be chosen sufficiently large so that the numerical solution closely approximates the terminal boundary conditions at ∞. Here, y ∞ has been fixed to 10.

The above systems of differential equations are converted into an equivalent initial value problem in order to use the fourth order Runge-Kutta method.

For the sake of brevity, here we have explained the solution procedure for the system (25) and (30) subject to boundary conditions ( 26) and (31) concisely.

The solution procedure for the system (25) and (33) subject to boundary conditions ( 26) and ( 34) is just similar.

To reduce the equations to equivalent first order system, we set:

u = z 1 , du dy = z 2 , d 2 u dy 2 = z 3 , d 3 u dy 3 = z 4 , d 4 u dy 4 = z 5 , d 5 u dy 5 = z 6 , θ = z 7 , dθ dy = z 8 , ( 35
)
to finally get:

z 1 = z 2 , z 1 (0) = λ[s 1 -α 1 V 0 s 2 + β 1 V 2 0 s 3 + β 3 s 3 1 -γ 1 V 3 0 s 4 -γV 0 s 2 1 s 2 ] z 2 = z 3 , z 2 (0) = s 1 z 3 = z 4 , z 3 (0) = s 2 z 4 = z 5 , z 4 (0) = s 3 z 5 = z 6 , z 5 (0) = s 4 z 6 = -M n z + V 0 z 3 + z 4 -α 1 V 0 z 5 + β 1 V 2 0 z 6 + β[2z 2 z 2 3 + z 2 2 z 4 ] -γV 0 [2z 3 3 + 6z 2 z 3 z 4 + z 2 2 z 5 ], z 6 (0) = s 5 z 7 = z 8 , z 7 (0) = 1 z 8 = -P r V 0 z 8 -P r E c [z 2 2 -α 1 V 0 z 3 z 2 + β 1 V 2 0 z 4 z 2 + βz 4 2 -γV 0 z 3 2 z 3 -γ 1 V 3 0 z 2 z 5 ], z 8 (0) = s 6 . ( 36
)
Here, the prime indicates differentiation with respect to y. The essence of this method is to reduce the boundary value problem to an equivalent initial value problem and then use a shooting numerical technique to guess s i , i = 1, 2, . . . 6, until the terminal boundary conditions z 1 = 1, z 7 = 0 are satisfied at y = y ∞ . Because of the non-linearities, which appear in the differential equations, the solution by the above method has been iterated and with reasonably close initial guesses to start the iterations, the convergence to the actual values within an accuracy of O (10 -6 ) is attained after only 9 to 11 iterations.

Results and discussion

Keeping in view of a proliferation of the number of parameters, we have restricted the discussions of the obtained results to the effects of emerging important flow parameters α 1 , β, γ, V 0 , M n , P r , E c , and the newly introduced slip parameter λ on the velocity and temperature fields. It is found that β, β 1 , and similarly γ, γ 1 have similar effects on the velocity and temperature fields. Further, one can refer the works of Hayat et al. [START_REF] Hayat | Flow of a fourth grade fluid[END_REF][START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF] and Sajid et al. [START_REF] Sajid | On the analytic solution of the steady flow of a fourth grade fluid[END_REF] for a precise discussion with the conventional no-slip boundary conditions.

As a preliminary test [START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF] to check the accuracy of the solution, it may be noted that for the no-slip and non magnetic case, at y = 0, Eqn. (25) reduces to:

V 0 d 2 u dy 2 - d 3 u dy 3 + α 1 V 0 d 4 u dy 4 -β 1 V 2 0 d 5 u dy 5 + γ 1 V 3 0 d 6 u
For all solutions computed for different ranges of the flow parameters, the left hand side of Eqn. ( 37) is found to be less than 10 -4 .

The variations of the non-dimensional velocity u with different flow parameters have been plotted in Figs. 234567for β 1 = 1 and γ 1 = 1. Fig. 2 depicts that in presence of slip, as the normal stress coefficient α 1 increases from zero, the velocity increases near the surface and then starts decreasing away from the surface of the plate, resulting in a 'cross-over' in the profile. It is noteworthy to mention that such a cross-over is absent for the no-slip case [START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF]. Figs. 3 and4 show that velocity decreases, or in other words, the momentum boundary layer thickness increases with an increase in the non-Newtonian fluid parameters β and γ respectively. These findings are in agreement with those reported by Hayat et al. [START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF] for the flow past a porous flat plate or by Islam et al. [START_REF] Islam | The optimal solution for the flow of a fourth-grade fluid with partial slip[END_REF] for a plane Poiseuille flow. In general, the fourth grade non-Newtonian fluid exhibits a thicker boundary layer than the Newtonian fluid, and increasing the higher order material parameters of non-Newtonian fluid causes further thickening of the boundary layer. Fig. 5 shows the variation of u with the nondimensional suction velocity. It is again interesting to find that in presence of slip, as the suction increases, the velocity near the plate increases, and then decreases gradually. This results in a cross-over in the velocity profile. However, in absence of slip, such a cross-over in not found, as reported by Sajid et al. [START_REF] Sajid | On the analytic solution of the steady flow of a fourth grade fluid[END_REF] and Hayat et al. [START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF]. Fig. 6 reveals that an increase in the magnetic interaction parameter M n increases the velocity throughout the domain of integration. The effect of slip on the velocity profile, when other flow parameters are kept constant, can be seen from Fig. 7. As was expected, an increase in slip results in an increase in the velocity near the surface of the plate. The velocity approaches towards the free stream value 1 with an increase in the slip, and flow behaves as though it were inviscid.

Another interesting quantity is the skin-friction coefficient at the wall. It is related to the wall shear stress, given by the Eqn. [START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF], which in the dimensionless form becomes:

T xy | y=0 = du dy -V 0 α 1 d 2 u dy 2 +V 2 0 β 1 d 3 u dy 3 + β 3 du dy 3 -V 3 0 γ 1 d 4 u dy 4 -V 0 γ du dy 2 d 2 u dy 2 y=0 ( 38
)
The 12 respectively depict that even in presence of slip, an increase in the non-Newtonian parameters α 1 , β and γ increase the temperature θ(y). These findings are in agreement with those reported by Hayat et al. [START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF]. The effect of the fourth garde fluid parameter γ on θ is not prominent throughout the domain of integration, as shown in Fig. 12. Figs. [START_REF] Nadeem | Effects of partial slip on a fourth-grade fluid with variable viscosity: An analytic solution[END_REF][START_REF] Sahoo | Flow and heat transfer of an electrically conducting third grade fluid past an infinite plate with partial slip[END_REF][START_REF] Shercliff | A text book of magnetohydrodynamics[END_REF] show that the temperature, and hence the thermal boundary layer thickness, decreases with an increase in the suction, magnetic parameter and slip parameter respectively. The suction has a prominent effect on the temperature profile. Figs. [START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF]-21 depict the effects of various flow parameters on the non-dimensional temperature g(y) for the insulated wall case. It is clear that all the flow parameters have similar effects of θ(y) and g(y). Fig. 21 shows that the slip has a substantial effect on the thermal boundary layer for the insulated wall case.

The dimensionless surface temperature gradient θ (0) in the constant surface temperature case, and the dimensionless surface temperature g(0) in the insulated wall case against different flow parameters have been tabulated in Table 1. It is clear that the shear thickening parameter β decreases the magnitude of the rate of heat transfer θ (0), when the other parameters are kept constant. It can be of interest to find that the heat transfer rate increases with γ in presence of comparatively large value (λ = 2) of the slip parameter. The slip parameter λ increase the rate of heat transfer at the surface of the plate. Moreover, it is observed that for constant E c , the value of |θ (0)| increases with P r , since the thermal boundary layer thickness δ T decreases with P r due to the estimate δ T ∼ 1 √ P r . Apparently, the surface temperature g(0) in the insulated wall case increases with E c , when other flow parameters are kept constant.

Conclusions

Flow and heat transfer of an electrically conducting fourth grade fluid past an infinite plate with partial slip boundary condition, subject to uniform suction has been considered in this paper. Due to the presence of higher order material derivatives of the strain tensor in the constitutive equation of the fourth grade fluid, the momentum equation gives rise to a highly nonlinear differential equation of order six with only two adherence boundary conditions. The resulting system of differential equations with inadequate boundary conditions has been solved by the shooting method along with the fourth order Runge-Kutta method and the Broyden's method. The use of the Broyden's method as a zero-finding algorithm has enhanced the efficiency of the shooting technique by reducing the computational time. The combined effects of the emerging non-Newtonian flow parameters and slip on the momentum and thermal boundary layers have been discussed elaborately and shown graphically. It is interesting to find the exceptional cross-overs in the velocity profile in the presence of slip. The fourth grade fluid parameter γ increases the momentum boundary layer thickness, whereas the slip parameter substantially decreases it. Similarly, the non-Newtonian fourth grade fluid parameters and the slip have opposite effects on the thermal boundary layer thickness.

The present study will be extended in a close future to take into account additional effects such as the dependance of the slip velocity on the normal stress, Joule heating in the energy equation or the thermal slip condition. As the shooting method limited the present investigation to moderate values of the flow parameters, a more effective and suitable numerical technique would be also suitable to see the effects of higher values of the fourth grade and slip parameters on the Blasius flow problem.

Fig. 1 .

 1 Fig. 1. A sketch of the Blasius flow problem.

  Figs.[START_REF] Wang | Unsteady flow of a fourth grade fluid due to an oscillating plate[END_REF]-21 indicate the effects of different emerging flow parameters on the non-dimensional temperature profiles θ(y) and g(y) for the constant surface temperature and the insulated wall cases, when β 1 = γ 1 = 1, P r = 3, and E c = 0.5. Figs 10-12 respectively depict that even in presence of slip, an increase in the non-Newtonian parameters α 1 , β and γ increase the temperature θ(y). These findings are in agreement with those reported by Hayat et al.[START_REF] Hayat | Heat transfer analysis of the steady flow of a fourth grade fluid[END_REF]. The effect of the fourth garde fluid parameter γ on θ is not prominent throughout the domain of integration, as shown in Fig.12. Figs.[START_REF] Nadeem | Effects of partial slip on a fourth-grade fluid with variable viscosity: An analytic solution[END_REF][START_REF] Sahoo | Flow and heat transfer of an electrically conducting third grade fluid past an infinite plate with partial slip[END_REF][START_REF] Shercliff | A text book of magnetohydrodynamics[END_REF] show that the temperature, and hence the thermal boundary layer thickness, decreases with an increase in the suction, magnetic parameter and slip parameter respectively. The suction has a prominent effect on the temperature profile. Figs.[START_REF] Garg | Stagnation point flow of a non-Newtonian fluid[END_REF]-21 depict the effects of various flow parameters on the non-dimensional temperature g(y) for the insulated wall case. It is clear that all the flow parameters have similar effects of θ(y) and g(y). Fig.21shows that the slip has a substantial effect on the thermal boundary layer for the insulated wall case.

Fig. 2 .

 2 Fig. 2. Variation of u with α 1 at β = 2.0, γ = 1.0, V 0 = 0.5, M n = 1.0 & λ = 1.0.

Fig. 3 .

 3 Fig. 3. Variation of u with β at α 1 = 1.0, γ = 1.0, V 0 = 0.5, M n = 1.0 & λ = 0.1.

Fig. 4 .

 4 Fig. 4. Variation of u with γ at α 1 = 1.0, β = 2.0, V 0 = 0.5, M n = 1.0 & λ = 0.1.

Fig. 5 .

 5 Fig. 5. Variation of u with V 0 at α 1 = 1.0, β = 2.0, γ = 1.0, M n = 1.0 & λ = 1.0.

Fig. 6 .

 6 Fig. 6. Variation of u with M n at α 1 = 1, β = 2, γ = 1, V 0 = 0.5 & λ = 0.5.

Fig. 7 .

 7 Fig. 7. Variation of u with λ atα 1 = 1, β = 2, γ = 1, V 0 = 0.5 & M n = 1.

Fig. 8 .

 8 Fig. 8. Variation of T xy | y=0 with γ atα 1 = 1, β = 2, V 0 = 0.5, M n = 1 & λ = 1.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Variation of T xy | y=0 with λ atα 1 = 1, β = 2, γ = 1, V 0 = 0.5 & M n = 1.

Fig. 11 .

 11 Fig. 11. Variation of θ with β atα 1 = 1, γ = 1, V 0 = 0.5, M n = 1 & λ = 0.1.

Fig. 12 .

 12 Fig. 12. Variation of θ with γ at α 1 = 1, β = 2, V 0 = 0.5, M n = 1 & λ = 0.1.

Fig. 13 .

 13 Fig. 13. Variation of θ withV 0 at α 1 = 1, β = 2, γ = 1, M n = 1 & λ = 1.

Fig. 14 .

 14 Fig. 14. Variation of θ with M n at α 1 = 1, β = 2, γ = 1, V 0 = 0.5 & λ = 0.5.

Fig. 15 .

 15 Fig. 15. Variation of θ with λ atα 1 = 1, β = 2, γ = 1, V 0 = 0.5 & M n = 1.

Fig. 16 .

 16 Fig. 16. Variation of g with α 1 at β = 2, γ = 1, V 0 = 0.5, M n = 1 & λ = 0.1.

Fig. 17 .

 17 Fig. 17. Variation of g with β atα 1 = 1, γ = 1, V 0 = 0.5, M n = 1 & λ = 0.1.

 1Variations of θ (0) and g(0) with different flow parameters.


Sébastien Poncet
email: poncet@l3m.univ-mrs.fr

Stéphane Viazzo

Adrien Aubert

Riccardo Da Soghe

Cosimo Bianchini

Turbulent Couette-Taylor flows with endwall effects: a numerical benchmark

Keywords:

teaching and research institutions in France or abroad, or from public or private research centers.

The accurate prediction of fluid flow within rotating systems has a primary role for the reliability and performance of rotating machineries. The selection of a suitable model to account for the effects of turbulence on such complex flows remains an open issue in the literature. This paper reports a numerical benchmark of different approaches available within commercial CFD solvers together with results obtained by means of in-house developed or open-source available research codes exploiting an innovative Reynolds Stress Model (RSM) closure, Large Eddy Simulation (LES) and a direct numerical simulation (DNS). The predictions are compared to the experimental data of Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] in an original enclosed Couette-Taylor apparatus with endcap rings. The results are discussed in details for both the mean and turbulent fields. A particular attention has been turned to the scaling of the turbulent angular momentum G with the Reynolds number Re. By DNS, G is found to be proportional to Re α , the exponent α = 1.9 being constant in our case for the whole range of Reynolds numbers. Most of the approaches predict quite well the good trends apart from the k-ω SST model, which provides relatively poor agreement with the experiments even for the mean tangential velocity profile. Even though no approach appears to be fully satisfactory, the innovative RSM closure offers the best overall agreement.

Keywords: Couette-Taylor flow ; Direct Numerical Simulation ; Large Eddy

I. INTRODUCTION

The present investigation is concerned with the numerical modeling of fluid flow in an enclosed Couette-Taylor system. The turbulent flow is confined between two coaxial cylinders, with an inner rotating cylinder and an outer stationary one. This kind of Couette-Taylor flows is of great importance, since they found many applications in process engineering (dynamic membrane filtration, rheology, UV disinfection), in astrophysics (accretion disks) and most of all in turbomachineries for bearings, asynchronous motor with axial ventilation, rotating heat exchangers or gas turbine engines among others. In rotating machineries, a better knowledge of the fluid flow is required to better predict the heat transfer distribution and thus to optimize the performances of the device.

Up to now, only very few works have been dedicated to real industrial flows, which imply very high rotation rates and narrow gap cavity. The difficulty to perform accurate measurements in such clearances and especially in the thin boundary layers along the cylinders is obvious and has slowed down the development of specific turbulence models. For many years, the turbulent flow regime has then been treated, in computations but also in experiments, as though the flow patterns within the cavity were axisymmetric. Treating the turbulent flow as statistically stationary implicitly assumed that the coherent structures present in the transitional modes disappeared or that their effect when averaged out was accounted for in the model without explicit treatment. Owen [START_REF] Owen | Flow and heat transfer in rotating disc systems[END_REF] first suggested that the flow within such cavity might not be steady. In the case of rotating disk systems with throughflow, the RANS model he used failed to converge for a certain range of parameters and he speculated that some form of organized large-scale structure may have been present in the flow. These three-dimensional patterns might be responsible for some discrepancies between experiments and axisymmetric computations. Even if three-dimensional calculations may be required in some specific configurations, their highly expensive computational cost makes the selection of a suitable turbulence model for the study of rotating cavity flows still pertinent. We can cite, among others, the work of Bazilevs and Akkerman [START_REF] Bazilevs | Large Eddy Simulation of turbulent Taylor-Couette flow using isogeometric analysis and the residual-based variational multiscale method[END_REF], who presented an application of the residual-based variational multiscale turbulence modeling methodology to the computation of turbulent Couette-Taylor flow at high Reynolds number. Its formulation globally conserves angular momentum, a feature that is felt to be important for flows dominated by rotation, and that is not shared by standard stabilized formulations of fluid flow. We can cite also the DNS results of Dong [START_REF] Dong | Turbulent flow between counter-rotating concentric cylinders: a direct numerical simulation study[END_REF] for turbulent Couette-Taylor flows between counterrotating cylinders.

The aim of the present work is then to perform a numerical benchmarking of turbulence modellings and Large Eddy Simulations (LES) against three-dimensional direct numerical simulation (DNS) results and experimental data. The geometrical configuration corresponds to the Couette-Taylor apparatus developed by Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF], who offered reliable velocity measurements in a rather complex geometry. Enclosed Couette-Taylor systems have been widely considered since one of the pioneering work by Wendt [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylindern[END_REF] because they offer an academic and relatively simple geometry to investigate the influence of rotation on turbulence for fundamental purposes or applications in astrophysics. When the inner cylinder is rotating while the outer one is at rest, the base flow depends mainly on three control parameters:

• the rotational Reynolds number of the inner cylinder based on the gap width Re

= ΩR 1 (R 2 -R 1 )/ν,
• the aspect ratio of the cavity Γ = H/[2(R 2 -R 1 )] and,

• its radius ratio η = R 1 /R 2 ,
where R 1 and R 2 are the radii of the inner and outer cylinders respectively and H is the cavity height. Biage and Campos [START_REF] Biage | Visualization study and quantitative velocity measurements in turbulent Taylor-Couette flow tagging: a description of the transition to turbulence[END_REF] proposed a classification of the systems depending on the value of their radius ratio η as follows: narrow-gap cavities for η > 0.67, middle-gap cavities for 0.33 < η < 0.67 and wide-gap cavities for η < 0.33. Most of the experiments up to now including the seminal one of Taylor [START_REF] Taylor | Stability of a viscous liquid contained between rotating cylinders[END_REF] have considered cavities with an aspect ratio larger than 100. Experiments with small aspect ratios Γ < 10 are less common. Kageyama et al. [START_REF] Kageyama | Numerical and experimental investigation of circulation in short cylinders[END_REF] considered both experimentally and numerically the flow in a small aspect ratio Couette-Taylor apparatus for η = 0.256 and Γ = 0.45. For such parameters, the vertical boundaries have a dynamical role on the secondary flow, that is why these authors proposed a novel approach to increase control over the velocity profile by increasing the number of endcap boundary conditions. Thus, they designed a new system with multiple independent endcap rings (up to 5) rotating at different speeds. Compared to classical Couette-Taylor systems, these rings allow for a reduction of the Ekman circulation in order to create a Couette-like velocity profile: V θ ∝ 1/r 2 . Ji et al. [START_REF] Ji | Hydrodynamic turbulence cannot transport angular momentum effectively in astrophysical disks[END_REF] developed a similar Couette-Taylor apparatus characterized by a low aspect ratio Γ = 1.06 and a low radius ratio η = 0.35 with two independent end-rings between the inner and outer cylinder. Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] reported velocity fluctuations and the fluctuation-driven radial transport of angular momentum in the same apparatus. Fluctuation levels and the mean specific angular momentum were found to be roughly constant over radius, in accordance with previous studies featuring narrower gaps.

Couette-Taylor flows as well as Rayleigh-Bénard convection are used also as canonical flows to analyse the transport properties in closed turbulent systems with a strong analogy between them [START_REF] Dubrulle | Momentum transport and torque scaling in Taylor-Couette flow from an analogy with turbulent convection[END_REF]. Besides their primary concern with astrophysical flows, the main goal is to understand the interplay between the boundary layers and the bulk flow. By torque measurements, Wendt [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylindern[END_REF] proposed some correlations for the dimensionless torque per unit length G for η = 0.68, 0.85, 0.935:

G = 1.45 η 3/2 (1 -η) 7/4 Re 1.5 , 400 < Re < 10 4 (1)
G = 0.23 η 3/2 (1 -η) 7/4 Re 1.7 , 10 4 < Re < 10 5 (2)
Other experiments indicated the existence of two scaling regimes, one for small Re, where the exponent, denoted α in the following, is 1.5 and one for larger Re, where α remains in the range [1.7 -1.8] [START_REF] Taylor | Stability of a viscous liquid contained between rotating cylinders[END_REF][START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylindern[END_REF]. From their torque measurements, Tong et al. [START_REF] Tong | Anisotropy in turbulent drag reduction[END_REF] found indeed that G ∝ Re 1.8 for 4 × 10 4 < Re < 4 × 10 5 and η = 0.448. The experimental data of Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] for the turbulent angular momentum scale with the Reynolds number to the power α = 1.7, whatever the flow being laminar or turbulent. Dubrulle and Hersant [START_REF] Dubrulle | Momentum transport and torque scaling in Taylor-Couette flow from an analogy with turbulent convection[END_REF] identified three logarithmic regimes using an analogy with turbulent Rayleigh-Bénard convection:

• Regime 1 for 400 < Re < 10 4 :

G = 1.46 η 3/2 (1 -η) 7/4 Re 3/2
(3)

• Regime 2 is not always observed in Taylor-Couette systems:

G = 2.12 η 2/3 (1 -η) 5/3 Re 5/3 ln[η 2 (1 -η)Re 2 /20] 2/3 (4)
• Regime 3 for 10 3 < Re < 10 6 :

G = 0.5 η 2 (1 -η) 3/2 Re 2 ln[η 2 (1 -η)Re 2 /10 4 ] 3/2
(5)

On the other hand, the experiments of Lathrop et al. [START_REF] Lathrop | Transition to shear-driven turbulence in Taylor-Couette flow[END_REF] and Lewis and Swinney [START_REF] Lewis | Velocity structure functions, scaling, and transitions in high Reynolds number Couette-Taylor flow[END_REF] yielded no region of constant exponent. For 800 < Re < 1.23 × 10 6 , these authors found indeed that α can vary between 1.6 and 1.86 in the turbulent regime. Eckhardt et al. [START_REF] Eckhardt | Torque scaling in turbulent TaylorCouette flow between independently rotating cylinders[END_REF] argued that there is no a single scaling but a superposition of different scalings: from α = 1.5 for small

Re to α = 2 for larger ones. The exact value of α is then still debated since the pioneering work of Wendt [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylindern[END_REF].

The paper is organized as follows: the geometrical configuration and the numerical modelings are described in Sections II and III respectively. The results of the numerical benchmarking are presented and discussed in Section IV in terms of flow structures, mean and turbulent fields. Finally some conclusions and closing remarks are provided in Section V.

II. DETAILS OF THE GEOMETRICAL CONFIGURATION AND FLOW CONTROL PARAMETERS

The fluid is confined between two concentric cylinders of radii R 1 = 0.071 m and R 2 = 0.203 m and height H = 0.28 m (Fig. 1). The inner cylinder is rotating at a constant angular rate Ω, while the outer cylinder is stationary. The originality of the apparatus developed by Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] is that each endcap is divided into two independently driven rings of equal radial extension. The inner cylinder and the inner ring rotate at Ω, while the other walls are stationary. In the following, the radial and axial coordinates are normalized as:

r * = (r -R 1 )/(R 2 -R 1
) and z * = z/H. The junction between the rings is then at r * = 0.5.

Under isothermal conditions, the main flow is governed by three flow control parameters:

the rotational Reynolds number Re based on half the hydraulic diameter

D H /2 = R 2 -R 1 ,
the aspect ratio of the cavity Γ and the radius ratio η defined as:

Re = ΩR 1 (R 2 -R 1 ) ν = 5 × 10 4 , 10 5 , 2 × 10 5 , 4 × 10 5 Γ = H D H = 1.06 η = R 1 R 2 = 0.35
The radius ratio remains in the range 0.33 < η < 0.67, so that the cavity is considered as a middle gap cavity [START_REF] Biage | Visualization study and quantitative velocity measurements in turbulent Taylor-Couette flow tagging: a description of the transition to turbulence[END_REF]. For enclosed systems characterized by a low aspect ratio Γ, such as in the experiment of Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF], the choice of the boundary conditions on the endcap disks is primordial and can favour the development of large Ekman recirculations Red walls are rotating and grey ones are stationary.

in the meridian plane (r,z). Depending on the rotation rates of the rings, the flow can be either in solid body rotation (constant angular velocity V θ /r of the fluid) or present some characteristics of quasi-keplerian flows (V θ /r ∝ r -3/2 ) or of pure Couette flows (V θ /r ∝ r -2 ).

III. NUMERICAL MODELLING

The present flow configuration presents several complexities (high rotation rate, confinement effect and singularities in the boundary conditions among others), which are very challenging for numerical approaches. Different numerical approaches including DNS, RANS and LES models either developed at the M2P2 laboratory or at the University of Florence or available within commercial CFD codes have been compared against the velocity measurements of [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF]. A brief description of the various numerical methods employed is given in the following.

A. Direct numerical simulation (DNS)

The incompressible fluid motion is governed by the three-dimensional Navier Stokes equations written in primitive variables for cylindrical coordinates (r, θ, z). No-slip boundary conditions are applied at all walls so that all the near-wall regions are explicitly computed.

Conservation equations are solved using a Fourier approximation in the homogeneous tan-gential direction. In both non homogeneous radial and vertical directions the solutions are approximated using a fourth-order compact finite-difference scheme [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF] on a structured mesh grid. The time advancement is second-order accurate and is based on the explicit Adams-Bashforth time-stepping for the convective terms and an implicit backward-Euler scheme for the viscous terms. The velocity-pressure coupling is solved using a two-step fractional scheme (predictor-corrector) reducing the problem at each time step to a set of two-dimensional Helmholtz equations.

A multidomain algorithm has also been developed in the non-periodic radial and axial directions with the aim to refine the mesh grid in some specific flow regions. The calculation domain has been then divided into four subdomains, two in each direction. The decomposition in the radial direction is to take into account the singularity in the boundary conditions between the endcap rings around r * = 0.5. The decomposition in the axial direction enables to refine the mesh grid around midplane where an intense shear layer is expected to be created by the recirculations of the secondary flow. The continuity conditions required at the interfaces of each subdomain, are computed by using the continuity influence matrix technique. This technique developed by Abide and Viazzo [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF] in a Cartesian frame ensures a direct connection of C 1 -type for the solution between two neighboring domains. This technique adapted here to cylindrical coordinate system allows to link the solutions in the consecutive subdomains only by one matricial calculation.

The cartesian solver has been fully validated in [START_REF] Beaubert | Large Eddy Simulations of plane turbulent impinging jets at moderate Reynolds numbers[END_REF] and in [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF] for the multidomain version. The numerical validation of the cylindrical solver has been done by Viazzo et al. [START_REF] Viazzo | High-order LES of confined rotor-stator flows[END_REF] in an interdisk rotor-stator cavity in its monodomain version. The numerical details are provided in Table I. Computations have been runned on the M2P2 cluster composed of 2

Xeon quadcore 3 GHz.

B. The Reynolds Stress Model of Elena & Schiestel [9]

The present Reynolds Stress Model (RSM) is based on one-point statistical modeling using a low Reynolds number second-order full stress transport closure sensitized to rotation effects by Elena and Schiestel [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF]. This approach allows for a detailed description of near-wall turbulence and is free from any eddy viscosity hypothesis. Four additional terms account for the implicit effects of rotation on turbulence and act only when the flow is sub-jected to strong rotation. Φ (R) is a part of the pressure-strain correlation term sensitized to the dimensionality tensor. D R is an inhomogeneous diffusion term, which slows down the tendency to bidimensionalization for wall bounded flows. B is a homogeneous source term, which rectifies the pressure-strain correlation and which acts only in case of strong rotation.

It produces spectral phase scrambling (angular dispersion). The rotation also reduces the energy transfer from large to small turbulent scales. It is modeled through an inverse flux J considered as isotropic for high Reynolds numbers. The full model is given in [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF].

The RSM model has been implemented in a finite-volume code using staggered grids for mean velocity components with axisymmetry hypothesis in the mean and non staggered grids for the Reynolds stress tensor. The code is unsteady elliptic. The velocity-pressure coupling is solved using the SIMPLER algorithm. In order to overcome stability problems, several stabilizing techniques are introduced in the numerical procedure. Also, the stress component equations are solved using matrix block tridiagonal solution to enhance stability using nonstaggered grids.

This version of the RSM has been fully validated in various rotating flow arrangements:

rotor-stator cavity with throughflow [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF][START_REF] Poncet | Batchelor versus Stewartson flow structures in a rotor-stator cavity with throughflow[END_REF] and Taylor-Couette-Poiseuille flows with heat transfer [START_REF] Poncet | Numerical modeling of fluid flow and heat transfer in a narrow Taylor-Couette-Poiseuille system[END_REF] among others. Unsteady axisymmetric calculations using a structured mesh have been performed on the M2P2 cluster with 2 Xeon quadcore at 3 GHz. All the numerical details have been given in Table I. Note that the mesh grid has not been refined close to the junctions between the rings.

C. The scale-adaptive simulation method within CFX

The turbulence model considered here is the innovative k-ω SST SAS originally proposed by Menter and Egorov [START_REF] Menter | A Scale-Adaptive Simulation Model using Two-Equation Models[END_REF] and available within CFX 12.0. The Scale-Adaptive Simulation I for the computational details concerning the SAS approach.

D. Large Eddy Simulation within OpenFOAM

OpenFOAM is an open-source CFD code composed of a set of libraries implemented to resolve continuum mechanics problems in 3D on unstructured grids [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. the modeled turbulence kinetic energy [START_REF] Menon | Effect of subgrid models on the computed interscale energy transfer in isotropic turbulence[END_REF][START_REF] Yoshizawa | Bridging between eddy-viscosity-type and second-order models using a twoscale DIA[END_REF]:

∂k ∂t + ∂(U i k) ∂x i - ∂ ∂x i ((ν + ν SGS σ k ) ∂k ∂x i ) = ν SGS ∥S∥ 2 - c e k 3/2 ∆ , ν SGS = c k √ k∆ (6)
The LES filter width is a simple cube root of the cell volume. A multi-block hexahedral grid mesh has been used and calculations have been performed on 6 cores of an Intel Xeon 3430 at 2.4 GHz for the case Re = 10 5 (see the numerical details in Table I).

E. Large Eddy Simulation within STAR-CCM+

Classical RANS models (k-ε, k-ω SST, RSM, v 2 -f ) including both high and low-Reynolds number approaches, but also a Detached Eddy Simulation based on the Spalart-Allmaras model available within Star CCM+ 5.04 have been compared for the present test case in [START_REF] Guillerm | Numerical predictions of flow field in closed and opened Taylor-Couette cavities, 4th European Conference for Aerospace Sciences[END_REF].

The LES model considered here uses a standard Smagorinsky model as subgrid scale model with a model constant C s equal to 0.1 offering a good compromise between the value 0.165 for homogeneous isotropic decaying turbulence and 0.07 for channel flows. The synthetic eddy method has been used to provide an initial perturbated flow field. A polyhedral unstructured mesh grid has been generated. To take into account the thin boundary layers developed along the cylinders and the rings, 12 prismatic layers are imposed on each wall on just 2% of the gap width to guarantee acceptable values for the wall coordinates (see Table I). Computations have been conducted on a HP Z400 workstation with a quadcore Intel Xeon CPU W 3520 processor at 2.67 GHz.

F. Computational details

All the numerical approaches together with the corresponding computational details are summed up in Table I. Note that N , K and M represent the number of mesh points in the radial, tangential and axial directions respectively. Some computations have been performed considering a half cavity or only a 90 • sector, which is justified by the absence of large scale vortical structures as highlighted by the DNS of Fukushima et al. [START_REF] Fukushima | Dynamics of largeand small-scale vortical structures in turbulent Taylor-Couette flow[END_REF] for Re up to 12000.

IV. RESULTS

The numerical results are compared to the velocity measurements of Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] in terms of mean tangential velocity and angular momentum and of the cross-correlation coefficient and the turbulent angular momentum. The DNS results provide also useful data to investigate the three-dimensional structure of the flow field. All the computed data have been averaged both in time and in the tangential direction.

A. Flow structures

The critical Reynolds number Re for the transition to turbulence in such a flow depends upon a large number of factors, including the aspect ratio of the cavity, the curvature and the flow history. Nevertheless, one commonly admitted value deduced from the Reynolds number dependence of the mean torque at the inner cylinder is around 10 4 [START_REF] Lathrop | Transition to shear-driven turbulence in Taylor-Couette flow[END_REF]. In the present case, the transition to a state of featureless turbulence appears for lower values of Re. At very low Re values, the flow is laminar unstable with some circular rolls of different sizes propagating towards the gap (Re = 5 × 10 3 , Fig. 2a). They are somewhat different from classical Taylor vortices, which would have occupied the whole cavity with a corresponding axial wave number of 2 (close to 2Γ). At Re = 10 4 (Fig. 2b), some spiral patterns aligned with the tangential direction appear the inner rings and around the midplane (z * = 0.5).

They coexist with the circular rolls on the top and bottom stationary rings. Some evidence of a disorganized flow structure can be seen from this Q-criterion map. When one increases further the Reynolds number, the flow gets clearly turbulent at Re = 5 × 10 4 (Fig. 2c) with fine scale eddies elongated in the tangential direction. Some circular rolls (not shown here) persist along the inner cylinder, with a structure similar to the toroidal vortices observed by Biage and Campos [START_REF] Biage | Visualization study and quantitative velocity measurements in turbulent Taylor-Couette flow tagging: a description of the transition to turbulence[END_REF] for Re = 2.04 × 10 5 , η = 0.38 and Γ = 10. At Re = 10 5 , the iso-values of the Q-criterion obtained by DNS (Fig. 2d) show that the flow is clearly turbulent with even thinner turbulent structures. There is absolutely no evidence of large scale vortical structures embedded in the flow apart from large circular patterns, which appear along the inner cylinder. Turbulence is mainly concentrated at midplane (z * = 0.5) because of the shear produced by the intense recirculation bubbles (Fig. 4). Thus, the critical Reynolds number for the transition to turbulence is not significantly modified by the presence of the endcap rings compared to classical Couette-Taylor flows.

B. Mean hydrodynamic field

In such configuration, the main flow is directed tangentially due to the high rotation rate of the inner cylinder and rings. Due to centrifugal and confinement effects, a secondary flow develops in the meridional plane (r,z). The results for the mean field are first qualitatively analysed considering bidimensional streamlines. When one increases the Reynolds number up to Re = 10 5 , very different patterns arise from the different computations (Fig. 4). The DNS (Fig. 4a), LES (Fig. 4d&e) and RSM (Fig. 4b) predict a main big vortex for each of the two sides with some differences between these approaches. Looking at the streamline patterns obtained by DNS, the k-ω SST SAS (Fig. 4c) and the LES of StarCCM+, asymmetry can develop at this Reynolds number. The secondary flow structure in the meridional plane gets more complex with the presence of small vortices arising along the inner cylinder or in the corners between the rings and the cylinders. The k-ω SST SAS model predicts a mean velocity pattern characterized by four rotating structures at the meridional plane, which is confirmed by computations performed by other k-ω models, indifferently in their unsteady or steady-state formulations and whatever the flow solver: CFX, OpenFOAM or STAR CCM+ (some results can be found in [START_REF] Guillerm | Numerical predictions of flow field in closed and opened Taylor-Couette cavities, 4th European Conference for Aerospace Sciences[END_REF]). The fluid pumped along the rotating disks is deviated away from the wall almost as soon as it reattaches on the fixed rings creating a high tangential velocity zone around mid-radii (Fig. 5). Such movement towards the center of the cylinder creates two counter rotating bubbles in the inner and outer zone. While the outer bubble is almost symmetric in the axial direction, the inner one is centered quite close to the rotating disk showing a quite stretched shape. Smaller recirculations attached to the inner cylinder can also be observed.

These smaller vortical structures are formed close to mid-height of the cylinder and are convected towards the rings. the k-ω SST SAS model predict a radial profile of the mean tangential velocity divided into three areas. Two very thin boundary layers develop along the cylinders with a thickness of the order of 1% of the gap width. They are separated by a central region where the tangential velocity component scales like: V * θ (r) = 0.43 × r * 2 -0.78 × r * + 0.59 (from the DNS). Hollerbach and Fournier [START_REF] Hollerbach | End-effects in rapidly rotating cylindrical Taylor-Couette flow[END_REF] pointed out that, at this Reynolds number, the shear layers developed close to the discontinuity should penetrate into the cavity and create a small discontinuity in the velocity profile. Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] attributed the absence of this discontinuity to an inadequate resolution of their measurements. The present data do not indicate any discontinuity in the velocity profile confirming the results of [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF].

The k-ω SST SAS model as well as all the other k-ω SST models considered in [START_REF] Guillerm | Numerical predictions of flow field in closed and opened Taylor-Couette cavities, 4th European Conference for Aerospace Sciences[END_REF] fail to predict the expected profile with a large overestimation of the mean tangential velocity within the gap (Fig. 5b). While most of all the models predict V * θ (r * = 0.5) between 0.3 (LES OpenFOAM) and 0.36 (RSM), the SAS model predicts V * θ (0.5) ≃ 0.48. A local maximum is also observed around r * ≃ 0.36 from the SAS results, which is not the case from the experiments. The unsatisfying agreement between the k-ω SST and the experimental data in terms of tangential velocity is probably related to the presence of the four rotating structures on the meridian plane. The velocity peak close to the radius r * = 0.36 on Figure 5 can be indeed justified by the tangential momentum transport operated by these inner vortices. Its behavior is similar to what obtained both experimentally and numerically by Kageyama et al. [START_REF] Kageyama | Numerical and experimental investigation of circulation in short cylinders[END_REF] for Re ≃ 10 6 , η = 0.256 and Γ = 0.45. The tangential velocity profile is nearly constant over radius in their case and the authors attribute it to the reduction of the Ekman cells due to endcap ring effects. It is noteworthy that the DNS predicts a local maximum close to the outer cylinder but due to the lack of experimental data very close to the wall, no definitive conclusion can be drawn. 5a to 5d, a noticeable influence of the Reynolds number on the mean velocity profiles can be observed. For Re ≥ 10 5 , the mean tangential velocity slightly increases for increasing values of Re but this effect remains weak. It is besides not observed from the LES (OpenFOAM) results, which underpredicts the tangential fluid velocity in that cases. The SAS model still overestimates V * θ in the main part of the gap. For Re = 5 × 10 4 (Fig. 5a), the experimental profile is much lower, which is characteristic of much weaker turbulence levels. All the numerical approaches fail to reproduce the good profile. The DNS seems to offer the best agreement with the experiment. The overestimation of the tangential velocity of the fluid may be attributed to higher turbulence intensities and to an earlier transition to the fully turbulent regime in the numerics. Surprisingly, the LES of OpenFOAM strongly overestimates the mean velocity profile in that case, with a profile quite similar to the laminar circular Couette profile, indicating much lower turbulence levels than in the experiments or in the DNS for examples.

By comparing Figures

The same data can also be plotted in terms of the angular momentum L θ = rV θ of the mean flow. Thus, Figure 6 presents the radial profiles of L θ at the same axial location z * = 1/4. It has been made dimensionless using the value of L θ at the inner cylinder

L θ (R 1 ) = ΩR 2 1 .
For Re = 5 × 10 4 , the angular momentum is rather constant over the radius with a measured value around 0.5, whereas a slight increase may be observed from the numerical profiles with L * θ (r * = 0.5) ≃ 0.62 for the DNS. Surprisingly, the LES of OpenFOAM predicts an angular momentum, which can be slightly higher than one. For Re ≥ 10 5 , the influence of Re on L * θ remains weak. Apart very close to the walls, the numerical approaches predict an angular momentum almost constant over the radius, while a slight increase with r * may be noticed from the measurements. The value at mid-gap is around 0.65, which is much larger than the usual value 0.5 obtained in classical Couette-Taylor flows [START_REF] Smith | Turbulent Couette flow between concentric cylinders at large Taylor numbers[END_REF]. It can be explained by considering the additional torque from the inner rotating rings, which rotate at the same rate as the inner cylinder but at a larger radius. Even if most of the models predict quite well the mean tangential velocity, the variations of the angular momentum with the radius highlight larger discrepancies mainly because of the normalization. The DNS slightly underestimates the slope of the profile in the central region, where the RSM of Elena and Schiestel [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF] offers the best agreement with the experimental data of Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF]. The k-ω SST SAS model strongly overestimates the mean angular momentum in the whole gap.

C. Turbulent statistics

Cross-correlation coefficient

Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] provided also the variations with the Reynolds number of the crosscorrelation coefficient C rθ defined as:

C rθ = v ′ θ v ′ r √ v ′2 θ √ v ′2 r (7)
The variation of C rθ with the Reynolds number Re are provided on Figure 7. Burin et al.

[5] observed that C rθ slightly decreases with Re as the flow becomes increasingly random.

As it can be seen from Figure 7, the time averaged value obtained by the different models at the monitoring point is close to 20% and does not vary so much with Re in the range [500, 4 × 10 5 ], in quite good agreement with [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF]. For Re = 5 × 10 4 , Smith and Townsend [START_REF] Smith | Turbulent Couette flow between concentric cylinders at large Taylor numbers[END_REF] found C rθ = 30 -40%. This value has been found using the RSM of StarCCM+ (not shown here), whereas all the other approaches predict C rθ = 0.24 or less. The high value obtained by [START_REF] Smith | Turbulent Couette flow between concentric cylinders at large Taylor numbers[END_REF] is due, according to [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF], to the presence of large scale vortical structures embedded in the flow, which has not been observed here by none of the present models.

In fact, the LES results (OpenFOAM) show that the cross-correlation coefficient strongly varies between -0.4 and 0.4 depending on the computed position in a meridian plane, as it can be seen from the C rθ -map in Figure 8. Thus, those punctual values have to be considered carefully and might be not so representative of the time averaged cross correlation in that region.

Transport of turbulent angular momentum

The measured velocity correlation is discussed also in terms of a dimensionless torque per unit length G, in which the viscous term has been neglected as in [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF]:

G = 2πr 2 v ′ θ v ′ r ν 2 (8)
Its variation with the Reynolds number is shown in Figure 9 for all the approaches and compared to the measurements of [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF]. It is also compared to the scalings proposed by Dubrulle and Hersant [START_REF] Dubrulle | Momentum transport and torque scaling in Taylor-Couette flow from an analogy with turbulent convection[END_REF] (Eq.??). For this range of Reynolds numbers, the DNS, LES (OpenFOAM) and RSM results confirm that G depends on Re following an unique power law: G ∝ Re α with α = 1.9, which suggests that the transport of angular momentum may be seen as a diffusive process. This value remains slightly larger than the one α = 1.7 deduced from the measurements of Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] or Wendt [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylindern[END_REF]. Tong et al. [START_REF] Tong | Anisotropy in turbulent drag reduction[END_REF] obtained α = 1.8 for η = 0.448 and 4 × 10 4 < Re < 4 × 10 5 . It shows that only one regime is obtained here corresponding to the regime 3 of Dubrulle and Hersant [START_REF] Dubrulle | Momentum transport and torque scaling in Taylor-Couette flow from an analogy with turbulent convection[END_REF]. Equation 5 would fit indeed quite well our numerical results with a lower radius ratio η.

Most of the models underestimate the values of G compared to the predictions of Equation [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] or to the measurements of [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylindern[END_REF]. Two main reasons can be advocated in the present case.

Firstly, the aspect ratio of the Burin et al.'s [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] apparatus is low compared to classical Couette-Taylor systems, which increases the influence of the endcap disks on the flow. The fact that only the inner rings are rotating with the inner cylinder may partly explain the reduced values of the turbulent angular momentum. Secondly, looking at Figure 10, one may also argue whether comparing a punctual value is an efficient strategy to verify turbulence levels when, as this is the case here, the local variability of the resolved part of G is high.

Figure 10 shows also how the modeled contributions are important only very close to the outer walls and are of the same order of magnitude as the resolved part in that regions. A preliminary calculation on a 5 degree sector using the LES within OpenFOAM has been performed providing much lower values of both C rθ and G with at least one order of magnitude. This finding confirms the importance of large tangential vertical structures on the behaviour of the flow.

V. CONCLUSION

Numerical predictions have been performed using a DNS code, the RSM closure of Elena and Schiestel [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF] and some LES and RANS models available within OpenFOAM, StarCCM+ or CFX. The computed results have been compared to experimental data available in the literature for an original Couette-Taylor arrangement initially developed by Burin et al. [START_REF] Burin | Local measurements of turbulent angular momentum transport in circular Couette flow[END_REF] characterized by endcap rings enclosing the cavity.

In the present cavity, the difficulty arises from the presence of two endcap rings attached to both cylinders. It appeared then crucial to take into account accurately the brutal variations in the boundary conditions on the upper and lower walls. A multidomain DNS code has then been developed providing very valuable data for comparisons with the experimental data and other numerical results. DNS showed that there is no evidence of Taylor vortices embedded in the turbulent flow apart very close to the inner cylinder where circular vortices have been observed. It explains why the 2D unsteady RSM [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF], which has been sensitized

to the implicit effects of rotation on turbulence, provides the best overall agreement in that case, where the rotation rate is quite high. On the contrary, all the k-ω SST models including the innovative SAS model fail to predict the mean tangential velocity profiles with large discrepancies even in the core of the flow. It can be attributed to the existence of four large recirculations in a meridian plane, whereas all the other models predict only two bubbles. All the models predict the good tendency for the turbulent angular momentum transport G, which is found to vary as Re 1.9 in agreement with previous approaches [START_REF] Dubrulle | Momentum transport and torque scaling in Taylor-Couette flow from an analogy with turbulent convection[END_REF].

Nevertheless, most of the models including the RSM [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF], the LES (OpenFOAM) but also the DNS underpredict the shear stress tensor component. The cross correlation coefficient C rθ is well predicted by all the approaches with almost a constant value around 0.2 (at a given location) for a wide range of Reynolds numbers. The LES provides very useful data showing the high variability of these two parameters.

For industrial applications, the DNS and LES remain computationally expensive. The RSM model [START_REF] Elena | Turbulence modeling of rotating confined flows[END_REF] showed a good behaviour in this challenging test case providing fair results

for the mean and turbulent fields. It offers then a good compromise between accuracy and calculation cost. More comparisons at higher Reynolds numbers are now required to confirm these tendencies.

Acknowledgments: The French part of this work was financially supported by the

  Turbulent Couette-Taylor flows with endwall effects: a numerical benchmark Sébastien Poncet, * Stéphane Viazzo, and Adrien Aubert Laboratoire M2P2, UMR 7340 CNRS -Aix-Marseille Université IMT La Jetée, Technopôle Château-Gombert, 38 rue F. Joliot-Curie, 13451 Marseille cédex 20 -FRANCE -Fax. +0033(0)491118502 Riccardo Da Soghe and Cosimo Bianchini Energy Engineering Department "S. Stecco", University of Florence 50139, via di S. Marta 3, Florence -ITALY (Dated: February 23, 2012)

Figure 1 .

 1 Figure 1. Geometrical configuration corresponding to the system considered by Burin et al. [5].

(

  SAS) is an improved URANS formulation, which allows the resolution of the turbulent spectrum in unstable flow conditions. The SAS concept is based on the introduction of the von Kármán length-scale into the turbulence scale equation. The information provided by the von Karman length-scale allows SAS models to dynamically adjust to resolved structures in a URANS simulation, which results in a LES-like behavior in unsteady regions of the flow field. At the same time, the model provides standard RANS capabilities in stable flow regions. The governing equations of the k -ω SST SAS model differ from those of the SST RANS model by an additional source term in the transport equation for the turbulence eddy frequency. Further details regarding the k -ω SST SAS model are available in the literature [19, 20]. The calculations have been performed under unsteady isothermal conditions. The time step have been selected assuring, for each operating velocity, a RMS Courant number lower than unity. The evolving fluid is considered ideally incompressible with constant transport properties. The convective fluxes were solved using a second order upwinded scheme. Automatic wall treatment (i.e. automatic blending from low to high Reynolds treatment as a function of the local first node wall distance) have been considered. A multi-block hexahedral grid has been used. Computations have been performed on 4 cores of an Intel Xeon 3430 at 2.4 GHz. Please refer to Table

  A pressure based segregated solver implementing the PISO loop to solve the pressure-velocity coupling is used in transient formulation. Both the convective and diffusive flux schemes employed are based on a purely linear interpolation, time advancement is achieved by means of the implicit backward Euler scheme, guaranteeing second order accuracy both in space and in time. The time step was generally chosen to achieve a maximum Courant number inside the domain close to unity to guarantee stability and improve the convergence rate. The LES subgrid effects are evaluated by means of a one-equation model exploiting a transport equation for

Figure 2 .

 2 Figure 2. Iso-values of the Q-criterion obtained by DNS for (a) Re = 5 × 10 3 , (b) Re = 10 4 , (c) Re = 5 × 10 4 and (d) Re = 10 5 .

Figure 3

 3 Figure 3 presents the streamline patterns obtained by DNS and the RSM for Re = 5 × 10 4 . In the bulk of the fluid, pressure balances the centrifugal force. The pressure being approximately uniform axially, the difference in centrifugal force between the bulk flow and the endcap rings generates a radial secondary flow in the (r,z) plane. The fluid is pumped along the rotating disks and being not blocked by the fixed rings, continues to

Figure 3 .

 3 Figure 3. Mean streamline patterns obtained at Re = 5 × 10 4 by (a) DNS and (b) RSM. The rotation axis is located on the left of each subfigure.

Figure 4 .

 4 Figure 4. Mean streamline patterns obtained at Re = 10 5 by (a) DNS, (b) RSM, (c) k-ω SST SAS, (d) LES (OpenFOAM) & (e) LES (Star CCM+). The rotation axis is located on the left of each subfigure.

Figure 5 Figure 5 .

 55 Figure 5 presents the radial profiles of the mean tangential velocity component V * θ = V θ /(ΩR 1 ) at the axial position z * = 1/4 for four Reynolds numbers from Re = 5 × 10 4 to Re = 4 × 10 5 . At Re = 10 5 for example (Fig.5b), all the numerical approaches apart from

Figure 6 .

 6 Figure 6. Radial profiles of the normalized angular momentum at z * = 1/4 for (a) Re = 5 × 10 4 , (b) Re = 10 5 , (c) Re = 2 × 10 5 and (d) Re = 4 × 10 5 .

Figure 7 .

 7 Figure 7. Variations of the cross correlation coefficient C rθ with the Reynolds number Re at z * = 1/4.

Figure 8 .

 8 Figure 8. Map of the cross-correlation coefficient C rθ in a meridian plane obtained by the LES of OpenFOAM at Re = 10 5 .

Figure 9 .

 9 Figure 9. Variations of the turbulent angular momentum transport G with the Reynolds number Re at z * = 1/4.

Figure 10 .

 10 Figure 10. Map of the turbulent angular momentum G (resolved and modelled parts) in a meridian plane obtained by the LES of OpenFOAM at Re = 10 5 .

Institut Carnot Star through the TCPtherm project. The postdoctoral position of A. Aubert is supported by Liebherr Aerospace Toulouse. They are here gratefully acknowledged. The authors are also grateful to Dr. Raphaël Guillerm and Prof. Bruno Facchini for their useful suggestions.