Numerical Simulation of Taylor—Couette—Poiseuille Flow at Re = 10,000

: A fully developed turbulent ﬂow in a concentric annulus, Re = 10,000, r i / r o = 0.5, with an inner rotating cylinder in the velocity range N = U ω / U b = 0 ÷ 4, is studied via a large-eddy simulation. Also, for comparison, simulations by steady-state, unstatiounary RANS k - ω SST (URANS), and Elliptic Blending Model (EBM) were made. The main focus of this study is on the effect of high rotation on the mean ﬂow, turbulence statistics, and vortex structure. Distribution of the tangential velocity and the Reynolds stress tensor change their behaviour at N > 0.5 ∼ 1. With rotation increases, the production of tangential ﬂuctuation becomes dominant over axial ones and the position of turbulent kinetic energy maximum shifts towards the wall into the buffer zone. URANS and EBM approaches show good agreement with LES in mean ﬂow, turbulent statistics, and integral parameters. The difference in pressure loss prediction between LES and URANS does not exceed 20%, but the average difference is about 11%. The EBM approach underestimates pressure losses up to 9% and on average not more than 5%. Vortex structures are described well by URANS.


Introduction
The annular flow phenomenon is ubiquitous in engineering and can be observed in diverse applications such as heat exchangers, chemical mixing devices, sliding bearings, well drilling, and turbomachinery.As a result, numerous experimental and numerical studies have been conducted to investigate the flow characteristics associated with the internal and external rotation of the cylinder (wall).Despite the geometric simplicity of this configuration, annular flow exhibits a complex three-dimensional boundary layer, swirling, and turbulence enhancement when the inner wall rotates.Conversely, when the outer wall rotates, the flow is suppressed, leading to a tendency toward laminarization [1].
Most of the previous research has focused on the inner wall rotational flow, apparently because of the interest in the shear-driven three-dimensional boundary layer.Previous experiments conducted by Nouri et al. (1993Nouri et al. ( -1997) ) [2][3][4], Escudier and Gouldson (1995) [5], and Rothe and Pfitzer (1997) [6] have provided valuable insights into the influence of inner wall rotation on turbulence structures in the inner wall region.A DNS study by Chung et al. (2002) [7], while limited to a non-rotating annular flow at a Reynolds number of 8900, discussed the effect of convex and concave lateral curvature on turbulence dynamics in regions adjacent to the inner and outer wall, respectively.Another DNS investigation conducted by Jung and Sung (2006) [8] used DNS to explore coherent structures near the inner rotating wall and their modification due to the work of centrifugal forces.In contrast, investigations of annular flows with a rotating outer wall are relatively scarce.
In the literature, a common observation is that the rotation rate in annular flows is usually low to moderate, where the velocity of the moving wall is either slightly above or equal to the fluid bulk velocity.For instance, Chung and Sung (2005) [9] conducted LES of annular flow at Re = 8900 and investigated three rotation rates (N = 0.2145, N = 0.429, N = 0.858), which exhibited discernible but not significant variations for different values of N.
Annular flow in the narrow gap was investigated by LES in studies [10][11][12].The authors [10] investigated the influence of the rotation on the mean flow and on the turbulence statistics, and studied the nature of the coherent structures appearing in the boundary layer and their influences on the heat transfers.Authors [12] showed a large contribution of the turbulent transport term on the Nusselt number results.In article [11], the turbulent kinetic energy transport in the near wall region was investigated for heat recovery system application.
Hanjalić and Launder in book [13] note that RANS models using the eddy viscosity hypothesis cannot properly describe flows with a high degree of swirl.This was also shown in a list of publications [14][15][16][17] about the Taylor-Couette-Poiseuille annular flow.
Authors of [16] compare results of LES and URANS approaches for Newtonian and power-law fluid.They prove and explain the reason for the applicability of the URANS approach to simulate swirling flow in the annulus.
In study [14], authors show good agreement with the first and second moment of flow simulated by RANS RSM [18] with the experimental data of [5].Applied model [18] is inspired by studies [19][20][21].
Here in the introduction, we mentioned only key studies related to current work.The work of Lockett should also be noted [22] for the extensive review of experimental and numerical work on the subject done before the 1990s.An extensive and fresh literature review of heat transport in rotating annular ducts was done in [23].The review includes articles from the first work by Taylor to recent experimental and numerical studies.A wide range of numerical methods was observed.There are not only articles related to heat transfer but also about isothermal annulus flow.
As mentioned above in the literature review, the Taylor-Couette-Poiseuille system has been extensively investigated for N < 1.One of the purposes of the present study is to extend the investigated range of the Reynolds number to Re= 10, 000 and the rotation number up to N = 4 using LES.Another goal of this work is to compare LES results with results obtained by URANS k-ω SST, RANS k-ω SST, and RANS EBM [24] approaches.The first and second flow moments obtained by the EBM approach and URANS were compared with LES.Integral parameters such as pressure losses and torque were compared for all above-mentioned approaches.

Problem Statement and Numerical Algorithm
Consider the incompressible flow throw an annular axisymmetric channel.The cylinder axis is aligned with the z-axis of the coordinate system.The inner cylinder, with radius r i , rotates clockwise (viewed toward z-direction) with constant velocity U ω , while the outer cylinder, with radius r o , is at rest.The ratio of the inner to outer cylinder radius is r i /r o = 0.5.The Reynolds number Re= 10, 000 is based on fluid bulk velocity U b and hydraulic diameter D h = 2(r o − r i ).The dimensionless rotation velocity of the inner cylinder takes values N = U ω /U b = 0, 0.2, 0.5, 1, 2 and 4.
The simulation is carried out using an in-house CFD code, which has been verified by DNS of the flow in the pipe [25] and by LES in annulus [16].The numerical algorithm is based on the finite volume method for an unstructured mesh.The SIMPLE-C algorithm was applied for the pressure correction procedure and the collocated grid arrangement with the Rhie-Chow interpolation.The system of linear algebraic equations for the pressure correction equation is solved using an algebraic multigrid solver.Periodic boundary conditions are applied for the flow direction.
In the LES approach, here the filtered Navier-Stokes equations and the continuity equation were closed by the dynamic Germano-Lilly eddy viscosity sub-grid model [26,27].The convective and diffusive terms were approximated by a second-order central difference scheme.The time step discretization was carried out by the second-order Crank-Nicholson scheme.The cell size on the channel walls in wall units was chosen so that in the radial direction ∆r + < 1, in the axial direction ∆z + < 16, and in the tangential direction (∆θr) + < 20.The sizes of the wall cells were chosen on the basis of recommendations from the book [28] and previous works on the same subject [1,9,11].From the walls of the computational domain to the middle of the gap, the mesh sizes were increased with a growth factor of no more than 5% from mesh to mesh.Thus, the mesh is redundant for LES, since the modelled part of TKE did not exceed 9% in most cases (Figure A1), with recommendations up to 20% [29].The time step was chosen based on the condition CFL < 0.9.Discussion about choosing of the time step and meshing is placed in Appendix A. The channel length L z = 4.5D h was chosen based on the results of [8], where a two-point correlation was used to determine the sufficiency of the channel length.The parameters of the computational mesh for LES are shown in Table 1.Other non-stationary simulations were carried out using URANS with the k-ω SST turbulence model.The numerical discretization schemes used are the same as in the LES approach and have the same radial distribution of computational nodes.In the tangential and axial directions, the number of nodes is reduced but satisfied ∆z + < 30, (∆θr) + < 30.The time step is also chosen based on the condition CFL< 0.8 and channel length L z = 4.5D h .The parameters of the computational mesh are shown in Table 2.In the LES and URANS approaches, all statics were obtained by averaging over time and in uniform directions, along the channel, and in the angular direction.The averaging time in all transient simulations was about 200D h /U b .
For steady-state simulations, two RANS approaches were also used for comparison.The first one is widely used in engineering applications: the k-ω SST with eddy viscosity assumption.The second approach, based on the Reynolds Stress Model (RSM), which models the full Reynolds stress tensor, is Elliptic Blending Model (EBM) [24].For both of these approaches, mesh with two cells in an axial direction was used, and the same radial distribution of cells was used as in the LES approach.The convective terms were approximated by a QUICK scheme.

Mean Flow Properties and Turbulence Statistics
It should be noted that the flow features and redistribution mechanics of the turbulence characteristics have been sufficiently described in earlier works.In this section, we focus on flow features that have not been noted before and compare the LES results with the results of the EBM simulation.
An increase in rotation increases the axial velocity gradient on the walls, making the axial velocity distribution U z flatter in the central region of the gap.The EBM model slightly overestimates the axial velocity near walls (Figure 1a), thereby underestimating it in the central region.The azimuthal velocity component U θ at N = 0.2-0.5 rotations has small local maxima at a distance 0.2-0.4 of the gap from the inner wall(Figure 1b).The EBM does not capture this small effect but tends to the LES curves at larger rotations.Let us discuss the axial velocity profile of the inner and outer walls in wall unites.In the non-rotating case, the axial velocity profile fits well to the wall-law near the inner and outer walls of the channel.As in work [9], the inner wall axial velocity profile U z,i appreciably decreases with increasing rotation N (Figure 2a).This is explained by the increasing in the axial friction on the inner wall with N.However, near the outer wall, these differences are noticeably smaller (Figure 2b).It is also remarkable that with rotation of N = 0.2, the axial velocity is higher than in the case without rotation, but with increasing rotation, as near the inner wall, it decreases noticeably in the buffer and logarithmic zone.The EBM approach does not capture this small effect at 0.2 rotation, but the general trend in behaviour and numerical values are captured correctly.
Wide theoretical investigation of swirling turbulent flow in axially symmetric flow was made in [30].Here, we should shortly repeat some logic and conclusions to describe the results of the simulation.Let us briefly overview relations between different terms in axially symmetric systems with centrifugal force.The swirl effect on turbulent flow is determined by ∂(rU θ )/∂r and u r u θ .The classification made in the three cases according to the signs of these two terms are The first case realizes when the outer cylinder rotates; this suppresses turbulence and leads to flow laminarization, which is shown in the study [1].When only the inner cylinder rotates, the u r u θ value is positive along the entire annular gap, but the value of ∂(rU θ )/∂r can change sign, thus, here we should consider the next two cases.The mechanism of the rotation effect on momentum transport is not straightforward; therefore, we show it using a block diagram for Case 2 and Case 3 in Figure 3, which is inspired by [30].Expressions for the production of the various components of the stress tensor are given below: .Through redistribution, it promotes the radial velocity fluctuations u r u r , which also are increased due to their positive production term P rr = 2u r u θ U θ /r.
Stress u r u θ > 0 that is responsible for tangential shear stress and torque is produced by rotation.That stress u r u θ is part of the production of radial velocity fluctuation u r u r and stress u θ u z .Stress u θ u z is generated by the two production terms −u r u θ ∂U z /∂r and −(u z u r /r)∂(rU θ )/∂r.In the current simulations, term (P θz ) I = −(u z u r )/r∂(rU θ )/∂r changes sign along the radius, but at the same time is significantly smaller than the value of (P θz ) I = −u r u θ ∂U z /∂r.By these means, Case 2 and Case 3 in the present conditions give the same results, and term u θ u z is promoted.Thus, the generation term of u z u r ((P zz ) I I = 2u θ u z U z /r) is also enhanced.Stress u r u θ is also enhanced by the increase in 2u θ u z U θ /r.
The term responsible for the axial momentum transfer from the walls u z u r is generated by the two production terms −u r u r ∂U z ∂r and u θ u z U θ r .The first one depends on the gradient of the axial velocity and radial fluctuations promoted by the axial and rotational velocity.The second term is promoted by rotation.Thus, the term u z u r responsible for axial friction (pressure loss) increases with both axial speed and rotational velocity.
The distribution of root-mean-square velocity fluctuations and Reynolds shear stresses obtained by LES and EBM are shown in Figure 4a.Without rotation, the axial component of the velocity fluctuations dominates over the other components and has well-defined maximums near the channel walls, similar to the flow in the channel.Rotation leads to an increase in all velocity fluctuation components as centrifugal force leads to flow destabilization and greater turbulence; the angular component appears with maximums near the channel walls like the axial component, and at high values of rotation their maximums become comparable in magnitude.It should also be noted that the position of the fluctuation peak gets closer to the wall with rotation as the viscous sublayer becomes thinner.The EBM approach describes the behavior of velocity fluctuation curves well, except for the axial fluctuations u z,rms , underestimating it.Thus, the total kinetic energy of velocity fluctuations is underestimated by EBM.At the transition from N = 1 to N = 2, the stress shape u θ u z changes behavior along the radius.At the region of the inner wall and N < 2, stress u θ u z is strictly negative, but at N ≥ 2 it changes sign and has a local maximum (Figure 4b).The mechanism of production and suppression of Reynolds stress tensor is well-described in article [30].Production of u θ u z consist of two terms, P θz = − ∂U z ∂r u r u θ − ∂U θ r ∂r u z u r r .The second term changes the sign, but it is too small to change the sign of production.It would be logical to assume that the redistribution term is responsible for the behavioural change of the component u θ u z .EBM misses such a small effect; below, we will see that this is not reflected in the integral characteristics.
Production terms P zz and P θθ have same structure (Formula (1)) and they depend on stress responsible for the transfer of respective momentum in radial direction and respective velocity gradient.It is evident that the tangential gradient ∂U θ /∂r grows faster than the axial gradient ∂U z /∂r, with increase in rotation N at inner wall.At the same time, rotation leads to an increase in wall friction and friction velocity U τ .Thus, production P θθ of the rotational component increases, but production P zz of axial component decreases near the inner wall in wall units.Hence, the fluctuations of the rotational velocity component grow faster with rotation (Figure 4a).It is also worth noting that the maximum production shifts from the buffer zone closer to the viscous sublayer, probably due to the rotation of the inner cylinder and centrifugal force (Figure 5).The sum of productions P zz and P θθ remains approximately constant near the rotating wall in near-wall units (P rr significantly smaller than both of them).A comparison of TKE distribution obtained by LES, URANS, and EBM is shown in Figure 6.As mentioned above, the EBM underestimates the fluctuations of the axial velocity, hence the total kinetic energy of the fluctuations, but correctly captures the position of its maximums.
It was shown earlier that k-ω SST as an eddy viscosity model does not work properly in rotational flow, underestimating TKE.Thus, we do not show curves for it in Figure 6a.Let us try to describe why the URANS approach improves simulation results.Due to fluid ejection from the inner wall of a fluid layer with a low level of modelled turbulent energy, the thickness of the near-wall layer increases, shear stresses on the walls decrease, and the modelled turbulent energy dissipates.A decrease in the modelled turbulent viscosity leads to flow instability and formation of vortices in the near-wall region.The ejection process is responsible for suppression of the modelled part of the turbulence and an increase in the resolvable part of the turbulence caused by large-scale spiral vortices.Thus, the mechanism for increasing the resolvable fluctuation fraction is proportional to the rotation; this can be seen from Figure 6a.The URANS approach applied here reproduces values and a maximum of TKE well near the inner wall, but does not catch in at the outer wall of the annular channel.Thus, we can conclude that URANS probably resolves vortex structures at passive wall worth than at a wall that actively produces vortices by rotation.Also, we should note that URANS at rotation N = 0-0.2does not resolve any vortices, simulating turbulent steady-state flow, i.e., the resolvable TKE fraction is zero (Figure 6b).At rotation N = 0.5, vortex structures are resolved only near the inner cylinder by URANS.The resolved fraction of the TKE value tends to unity near the wall, the location of flow instability generation, and noticeably decreases towards the center of the gap.Interesting to note that the total resolved TKE increases at N = 0-1 and decreases at N > 1 (Figure 6b).In Chang's work [9], it was shown that as N grows, the structural parameter a 1 also grows, which indicates an increase in shear stresses relative to normal stresses, i.e., the efficiency of shear stress generation grows.In our simulations, this conclusion is confirmed for N < 1, but at larger rotations, the parameter a 1 decreases (Figure 7a).EBM cannot fully replicate the LES curves for parameter a 1 , but the general character of the dependence is captured, which is a very good result for turbulent flow with swirling flow (Figure 7b).Also, it should be noted that EBM is not sensitive to small rotations, because curves for N = 0 and N = 0.2 are the same for EBM.

Flow Patterns, Vortical and Turbulence Structures
Briefly, flow visualization by λ 2 criterion were shown in [9] and [8] for Re = 8900, r i /r o = 0.5, and N = 0, 0.429.A more detailed flow structures analysis was made by [11] for the radii ratio r i /r o = 0.809 at the maximum Reynolds number Re = 8776 with rotation N = 0.85.Here, we concentrate our attention on high rotation and comparison of flow structures obtained by LES and URANS.Classical Görtler instability was described when the boundary layer thickness is comparable to the radius of the curvature; under the action of centrifugal force, a pressure change in the boundary layer occurs [31].In our case, the rotation of the inner cylinder generates a centrifugal force sufficient to destabilise the flow.Centrifugal instability of the boundary layer leads to the subsequent formation of Görtler vortices.The instability leads to the ejection of fluid from the channel walls, forming a pair of vortices and a region of reduced friction on the wall between them.A region of increased friction is formed between two pairs of such vortices, where fluid injects back.Schematically, ejection/injection, zone of higher/lower friction, and vortex formation are shown in Figure 8a.Visualisation of such vortices are shown in Figure 8b.The velocity fluctuations is higher on the inner wall of the channel (Figure 6a), as it is an active generator of turbulence due to its rotation, while vortex formation on the outer wall is secondary, as it occurs due to rotation of the whole liquid near the resting wall.As a consequence, we should expect a higher density of vortices on the inner wall, as we can see from the LES and URANS visualizations (Figure 9a).The vortices are also distributed non-uniformly along the walls, forming helical clouds or clusters on both the inside and outside walls of the channel.
Vortices on different walls of the channel are inclined in opposite directions.On the inner cylinder, the spiral of the vortices is wrapped against the rotation of the flow and on the outer cylinder in the direction of the flow rotation.The eddies form in pairs and are carried downstream.As the inner wall is rotating, i.e., moving in a tangential direction relative to the liquid, the vortices formed will be carried away in the opposite direction to its movement.Thus, the direction of the vortices and the twisting of their spiral clusters is opposite to the rotation of the flow.In the case of the outer wall, the situation is exactly the opposite.Non-uniform distribution of vortices along the channel was observed in the Couette-Taylor flow in [32], where only the inner cylinder rotated without axial flow.If, in the case of purely rotational flow, the vortex clusters resemble a torus in case of a combination of rotational flow and axial flow, they appear as spirals.The angle of inclination of spirals in relation to flow direction increases with an increase in rotation speed, and at significant dominance of rotational flow spiral clots of vortices can possibly degenerate into toroidal ones.Such a phenomenon has already been observed, however, for laminar flows at Re ∼ 300 and rotation N > 3 when the phenomenon is observed for laminar regimes when the vortices of the Görtler-type degenerated into Taylor vortices with increasing rotation [33].
The distribution of axial friction at the inner and outer cylinder is shown in Figure 9b,c, respectively.Strips of inhomogeneity in the wall friction are formed by vortices.The increase in vortex angle and decrease in the vortices' size with increasing rotation N are reflected by inhomogeneities in the wall friction.At rotation N = 4, the vortices are deviated very far from the channel axis so that the axial friction can become negative; such areas are marked in blue in Figure 9b.b,c).
LES resolves significantly smaller vorticity structures than URANS, thus, URANS visualizations look like LES flow but filtered at a larger scale (Figure 9a).Groups of small vortices united into one larger and longer vortex.This is clearly visible in cases with N = 2 and N = 4, where clusters of small vortices at the outer wall (red-colored) obtained by LES transformed into larger and longer continuous spiral vortices.The axial shear stress distribution obtained with URANS saves the slope angles of the inhomogeneities formed by the vortex structures, but the distribution is greatly smoothed due to the approach itself (Figure 9b,c).Overall, it should be noted that URANS correctly shows the structure of the flow, albeit in an enlarged form and not in as much detail, but it saves computational resources considerably.

Integral Parameters
In this paragraph, we will consider the integral parameters, axial friction, responsible for pressure loss and tangential friction, responsible for the torque acting on the inner and outer cylinders.Values of axial friction on the inner and outer cylinder are shown in Table 3.The results of LES show that without rotation, the friction on the inner and outer cylinder is comparable.As the rotational velocity increases, the friction on the inner cylinder grows faster because vortices are more intensely formed on it, which increases energy dissipation.Also, at rotation N < 0.5, friction on both walls increases non-linearly with rotation, while at higher values of rotation the dependence is practically linear (Table 3).We have seen a rearrangement of the rotational velocity profile (Figure 1b) and the behaviour of structural parameter a 1 (Figure 7a) changes at N > 0.5.It changes the momentum transfer mechanism and friction dependence in rotation.Also, this is true for the total pressure loss (Figure 10a).Most authors disregard the analysis of the tangential friction or the torque acting on the cylinders, although this may be an important parameter both from a fundamental point of view for tuning turbulence models, and from a technical point of view for solving engineering problems.Here we show non-dimension torque T nondimensionalized using analytical solutions that are proportional to viscosity and rotation.With good accuracy, it can be assumed that the dimensionless torque is proportional to rotation at N > 0.5, then, its dimensional value is proportional to the square of the rotation frequency of the inner cylinder (Figure 10b).As has been shown many times, models with the eddy viscosity hypothesis are unable to adequately capture swirling flows [13].Thus, the k-ω SST model underestimates friction at the inner cylinder more where more vortices are formed (Table 3).At rotations, N < 0.5 gives an acceptable error in pressure loss (8-9%), but with increasing rotation, the underestimation grows to 37%.Average difference in pressure loss is 26%.The torque is underestimated by this approach from 20% up to 60%, with an average value of 46%.Thus, tangential friction is much worse predicted in a steady-state simulation by the k-ω SST model.
The disadvantages of the stationary k-ω SST approach are partially compensated for by the non-stationary URANS approach, which is capable of resolving non-stationary vortex structures [16,34,35].The behaviour of pressure losses by URANS at a low rotation (N ≤ 0.2) is similar to RANS, clearly seen in Figure 10a and Table 3.Thus, URANS at small rotations is not sensible to rotation.This is happening because high eddy viscosity suppresses the formation of instabilities and vortices.As rotation increases, the difference in URANS with LES in pressure loss increases, but not exceeding 20%, with an average value of 11%.The torque is predicted by URANS quite well, with maximum difference less than 12% and an average value around 8% (Figure 10b).
The EBM reproduces Reynolds stresses well (Figure 4) and thus proved to be suitable for describing swirling flows.It is interesting to note that the EBM approach overestimates the axial friction at the axial flow on the channel walls, as does the k-ω SST approach, but as the rotation increases, the difference first increases and then almost disappears.The maximum difference with LES for wall friction is 14% and for total pressure loss 9% (average 5%).The torque is predicted satisfactorily, with an error of up to 19% (average 16%).

Conclusions
LES simulations of a fully-developed turbulent annular flow at Re = 10, 000 and r o /r i = 0.5 with an axially rotating inner wall at a range of rotation rates N = U θ /U b = 0, 0.2, 0.5, 1.0, 2.0, and 4 was carried out to study the effects of rotation on the flow, turbulence properties, and eddy structures.The observation and analysis of the mean velocity, turbulent stress fields, the vortical structures visualized by the λ 2 , and integral parameters, as well as comparisons made between them and the results of URANS k-ω SST, RANS k-ω SST, and EBM, lead to the following conclusions: • Rotation leads to a thinning of the viscous sublayer and, as a consequence, a widening of the mean axial profile and an increase in gradients at the wall.The axial velocity component distribution in wall units decreases due to an increase in friction on the wall with increasing rotation N.

•
Rotation decreases axial fluctuation production and increases tangential fluctuation production in wall units, while the maximum value of total production is weakly dependent on rotation, but its position shifts towards the wall into the buffer zone.

•
With increase in rotation and N > 0.5, the following changes are observed: -The tangential velocity component changes its shape, eliminating the local maximum at about one-third of the channel width from the inner wall.

-
The profile of the Reynolds stress tensor component u θ u z changes its monotonicity and becomes positive in areas where it was not before.

-
The structural parameter a 1 begins to decrease with increasing rotation, i.e., the production of shear components becomes less efficient.

•
The applicability of URANS and EBM techniques for describing first and second statistical moments of velocity fluctuation as well as integral characteristics of the flow is shown.URANS describes the vortex structures of the flow well, however, in an enlarged form.
As already written in the introduction, flow in an annular channel is encountered in many applications such as chemical reactors, heat exchangers, turbomachines, and drilling.The fluids used for drilling have non-Newtonian rheology.The turbulent flow of such fluids is not widely studied.Thus, the turbulent flow of non-Newtonian fluids in an annular channel is both interesting from a fundamental point of view and has important applications.Research of this type should be a logical continuation of this study.where f is the numerical solution, f 0 is the extrapolated solution, and a, b, c, d are the coefficients of the series.The parameters f 0 , a, b, c, and d can be found by solving an optimization problem on the set of value f at different parameters, ∆h and ∆t.This extrapolation was applied, which gave extrapolated values of C f and T.These values are included in Tables A1 and A2 for ∆t = 0 and ∆h = 0.The error of f (∆h, ∆t) with respect to the relatively extrapolated value f 0 is determined as follows: and included in Tables A1 and A2.The error estimate of the LES method for the pressure loss coefficient C f is 2% and 8.6% for the torque coefficient T. For the URANS method, the corresponding error estimations are 8.3% and 17%, respectively.The Roache's grid convergence index (GCI) is also evaluated, assuming second order of convergence in space and time [36].The GCI for the fine grid numerical solution is defined as: where F s = 3 is a factor of safety, p = 2 order of accuracy, f 1 is value obtained using fine time or space resolution, and f 2 is value obtained using coarser time or space resolution.The calculation of the convergence index GCI requires a solution on two grids or with two different time steps and the order of discretisation accuracy.The estimates given in Tables A1 and A2 are made with respect to the smallest time step and space resolution.Since the GCI convergence index is given separately for space resolution and time resolution, we will take the larger one for estimation.The convergence index GCI of the LES method for the pressure loss is 3.6% and for the torque it is 6.2%.For the URANS method, the corresponding estimates are 6.8% and 11.5%, respectively.
LES for the selected parameters showed a small error in determining the pressure losses andan acceptable error in determining torque.With a sufficiently detailed mesh, applicable even for the unresolved LES, URANS shows a rather high error.Apparently, an excessively detailed mesh has a negative effect on the URANS solution.The problem of grid convergence of the URANS method in solving problems with centrifugal force requires a special study, which is beyond the scope of this paper.Thus, the choice of grid and time step can be considered reasonable for LES and satisfactory for URANS.

Figure 1 .
Figure 1.Axial (a) and tangential (b) velocity normalized by the bulk and inner cylinder velocity, correspondently.Solid curves show LES and the dashed ones shows EBM.

Figure 2 .
Figure 2. Axial velocity in wall units.Inner wall region (a) and outer wall region (b).Solid curves show LES and the dashed curved shows EBM.

Figure 3 .
Figure 3.The mechanism of swirl effects on momentum transport for Case 2 and Case 3. The axial fluctuations u z u z promoted by axial flow gradient ∂U z ∂r with u z u r because they change sign at the same time as the walls, making production positive (P zz = −2u z u r ∂U z ∂r ).Through redistribution, it promotes the radial velocity fluctuations u r u r , which also are increased due to their positive production term P rr = 2u r u θ U θ /r.Stress u r u θ > 0 that is responsible for tangential shear stress and torque is produced by rotation.That stress u r u θ is part of the production of radial velocity fluctuation u r u r and stress u θ u z .Stress u θ u z is generated by the two production terms −u r u θ ∂U z /∂r and −(u z u r /r)∂(rU θ )/∂r.In the current simulations, term (P θz ) I = −(u z u r )/r∂(rU θ )/∂r changes sign along the radius, but at the same time is significantly smaller than the value of (P θz ) I = −u r u θ ∂U z /∂r.By these means, Case 2 and Case 3 in the present conditions r i )/(r o − r i ) r i )/(r o − r i )

Figure 4 .
Figure 4. Distributions of root-mean-square velocity fluctuations (a) and Reynolds shear stress (b).Solid curves show LES, and the dashed one shows EBM.

Figure 5 .
Figure 5. Production of u z u z (a) and u θ u θ (b) in wall units near the inner wall.
r i )/(r o − r i ) r i )/(r o − r i )

Figure 8 .
Figure 8. Vorticity formation on inner wall, their structure, and axial wall stress.(a) Schematic.(b) Flow visualization by λ 2 at N = 2. Red-coloured vortices have a positive rotation, blue ones have a negative rotation.Greater axial stresses on the wall are shown in a darker tone.

Figure 9 .
Figure 9. Flow structure visualization by LES and URANS.(a) Vortices visualisation by λ 2 isosurfases.Coloured by radius, thus, blue corresponds to the inner cylinder and red corresponds to the outer one.(b,c) Axial shear stress at the inner and outer cylinders, respectively.Red is positive, blue is negative friction.The same scale is applied for Figures (b,c).

Figure 10 .
Figure 10.Dependence of the skin friction coefficient (a) and torque acting on cylinders (b) in rotation.

Table 2 .
Mesh resolution applied at URANS for different rotation rates N.

Table 3 .
Skin friction factor at the inner (C f ,i ) and outer (C f ,o ) cylinder.The difference with LES is written in brackets.

Table A1 .
Dependence of pressure loss coefficient C f and torque coefficient T on mesh and time step resolution for N = 1 obtained by LES.

Table A2 .
Dependence of pressure loss coefficient C f and torque coefficient T on mesh and time step resolution for N = 1 obtained by URANS.(∆h, ∆t) = f 0 + a∆h + b∆h 2 + c∆t + c∆t 2 , (A1) f