Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 3, pp. 685-704, Warsaw 2007 INSTABILITY OF THE FLOW IN ROTATING CAVITY Ewa Tuliszka-Sznitko Artur Zieliński Technical University of Poznań , Institute of Thermal Engineering, Poland e-mail: sznitko@sol.put.poznan.pl The transitional flow in rotating cavity is investigated numerically by Direct Numerical Simulation (DNS), Large Eddy Simulation (LES) and theoretical (LSA) methods. LSA results coupled with accurate numeri- cal DNS and LES computations based on an efficient pseudo-spectral Chebyshev-Fouriermethod, brings new insight into the spatio-temporal characteristics of the isothermal and not-isothermal flows in the rotating cavities. DNS and LES computations have been performed for a wide range of Reynolds numbers to analyze different stages of the transition process. Computations have been performed for different geometrical configurations, including co- and counter-rotating cavitieswith through- flow. Key words: laminar-turbulent transition, rotating cavity, direct method 1. Introduction The flow in rotating disks system is not only a subject of fundamental in- terest but it is also a topic of practical importance. Typical configurations are cavities between the rotating compressors and turbines’ disks. Numerous works have been recently devoted to investigation of the instabilities, associa- ted with a single disk flow and with a differently rotating disks flow, cf. Serre et al. (2001a,b, 2004), Lingwood (1997), Tuliszka-Sznitko and Soong (2000), Tuliszka-Sznitko et al. (2002), Itoh (1991). For the high rotation rate, the flow in a rotor/stator cavity consists of two boundary layers, of the Ekman type on the rotating disk and of the Bödewadt type on the stationary disk, sepa- rated by an inviscid rotating core. The transition process in both layers is related to the type I and type II generic linear instabilities. Type I instability is the inviscid instability. Themechanismof type II instability is related to the combined effects of the Coriolis and viscous forces. The instability structures in the rotor/stator cavity were investigated numerically and experimentally 686 E. Tuliszka-Sznitko, A. Zieliński among others by Dijkstra and van Heijst (1983), Serre et al. (2001a, 2004), Gauthier et al. (2002), Schouveiler et al. (1999, 2001). The stability of co- and counter-rotating disks’ cavity was studied byGauthier et al. (2002) andMoisy et al. (2004). Non-isothermal flow conditions were also considered (Mochizuki andYang Wen-Jei, 1986; Tuliszka-Sznitko and Soong, 2000; Tuliszka-Sznitko and Zie- liński, 2006a), showing that the thermal effects and the rotation-induced bu- oyancy become influential on the stability characteristics and on the critical conditions. The most interesting from point of view of applications are cases with a superimposed radial outflow of cooling air, which are often used to remove heat from co-rotating turbine discs. Many experimental and numeri- cal investigations have been conducted in an attempt to understand this kind of flows and to use expensive cooling air more effectively (Bohn et al., 1994; Owen and Rogers, 1995; Chew et al., 1984). DNS method is mostly used to investigate transitional flow between ro- tating disks by solving 3D Navier Stokes, continuity and, in the case of not- isothermal flow, energy equations. With DNS all relevant length-scales and time-scales of motion are resolved. However, the required computational re- courses are large because CPU time increases as Re3R. Due to the CPU time requirements, computations are limited to flows of low or moderate Reynolds number. It means that only transitional or weakly turbulent flow is computed by DNS. For high Reynolds number the RANS method is commonly used in which only themean flow is computed by solvingReynolds – averagedNavier- Stokes equations and the effect of turbulence on the mean flow is modeled. The computational cost is small in comparison with the cost of DNS, however capability of themethod is limited. The existingmodels of turbulence,mostly using the isotropy assumption, are not sufficient for numerical investigations of the rotating flows. Large Eddy Simulation (LES) is an intermediate ap- proach in which the effect of large scales is directly computed and only the small subgrid scales are modeled. Since small scales are more isotropic than the large ones, it should be possible to parameterize them using simpler and more universal models than the RANSmodels. In the paper, we present our theoretical and numerical investigations on the transitional isothermal and non-isothermal flows in rotating cavity. Our early LSAcomputationswere used to enlighten theDNS andLES resultswith respect to type I and II instabilities.Moreover, the absolute instability regions which are supposed to be strongly connected with the turbulent breakdown process were identified by the LSA method. The three-dimensional DNS and LES computations have been performed for an annular roror/rotor cavity of aspect ratio L = 9 and curvature parameters Rm = 1.5, 3.0. DNS results obtained for the cylindrical cavity of aspect ratio L = 5 (published in Serre et al., 2004) are also discussed. Instability of the flow in rotating cavity 687 For the not-isothermal flow the Boussinesq approximation was used. Non- isothermal computations allowed us to analyze the influence of the thermal Rossby number on the instability structures and on critical parameters. The distributions of the localNusselt numbers along the radius of the disks are pre- sented. The isothermal and non-isothermal flows with a superimposed radial outflow of cooling air are analyzed. In the present paper the geometrical andmathematical models are presen- ted in Section 2. The LSA, DNS and LES numerical solution techniques are described in Section 3. In Sections 4 and 5 the obtained results are analyzed. In Section 6 concluding remarks are given. 2. Geometrical and mathematical model The geometricalmodel is a cavity between co- and counter-rotating diskswith and without throughflow (Fig.1). The cavity, bounded by an outer cylinder of radius R∗1 and height 2h is called cylindrical and the cavity bounded by inner and outer cylinders of radiuses R∗0 and R ∗ 1 is called annular. The faster rotating disk rotates at uniform angular velocity Ω∗1 = Ω ∗ 1ez, ez being the unit vector. The slower rotating disk rotates at angular velocity Ω∗2 = sΩ ∗ 1. Positive smeans that both disks rotate in the same direction and negative s means that disks rotate in opposite directions. The flow is controlled by the following physical parameters: the Reynolds number based on the external radius of the disks and on the angular velocity of the faster rotating disk Ω∗1, ReR =R ∗2 1 Ω ∗ 1/ν ∗, the local Reynolds number Reδ = r ∗/δ = √ r ∗2Ω∗1/ν ∗, the aspect ratio L = (R∗1 −R ∗ 0)/2h and the curvature parameter Rm = (R ∗ 1 + R∗0)/(R ∗ 1 −R ∗ 0). In the case of similarity solution considered in Section 3, radiuses of the disks are infinite (L→∞,Rm =1). The temperatures of the upper and lower disk are denoted by T∗1 and T ∗ 2 , respectively. Fig. 1. Schematic picture of a rotating cavity 688 E. Tuliszka-Sznitko, A. Zieliński The flow is described by the continuity, Navier-Stokes and energy equ- ations. The governing equations are written in a cylindrical polar coordinate system (r∗,z∗,ϕ), with respect to the rotating frame of reference (it is appro- priate for explicit presentation of the Coriolis and centrifugal forces) ∇·V ∗ =0 ρ∗ ∂V ∗ ∂t∗ +ρ∗(V ∗ ·∇)V ∗+ρ∗Ω∗1× (Ω ∗ 1×r ∗)+2ρ∗Ω∗1×V ∗ = =−∇p∗+µ∗∆V ∗ (2.1) ∂T∗ ∂t∗ +(V ∗ ·∇)T∗ = a∗∆T∗ where: t∗ is time, T∗ is temperature, p∗ is pressure, V ∗ is the velocity vector, a∗ is the thermal diffusivity. The time, space and velocity are normalized as follows: (Ω∗1) −1,hand Ω∗1R ∗ 1.Thedimensionless axial co-ordinate is z= z ∗/h, z ∈ [−1,1]. The radial coordinate is normalized additionally to obtain the domain [−1,1] requested by the spectral method based on the Chebyshev polynomials r = (r∗/Lh−Rm). The dimensionless component of velocity in radial, azimuthal and axial directions are u, v and w, respectively. The dimensionless temperature is defined as Θ = (T∗ −T∗1 )/(T ∗ 2 −T ∗ 1). To take into account the buoyancy effects induced by the involved body forces the Boussinesq approximation is used. No slip boundary condition is applied to all rigid walls, so u = w = 0. For the azimuthal velocity component, the boundary conditions are v = (Rm + r)/(Rm +1) on the top of the faster rotating disk, and v= s(Rm+r)/(Rm+1) on the slower rotating disk. 3. Numerical methods 3.1. Linear stability theory In LSA method the flow parameters are decomposed into the stationary basic state andnon-stationary disturbancefield.A similaritymodel of thermal flow with assumption of the Boussinesq fluid was formulated for generating basic solutions of axially symmetric flows (Tuliszka-Sznitko and Soong, 2000). InLSAwe assume that the perturbation quantities have the following normal- mode form [u′,v′,w′,p′,τ ′]⊤ = [û, v̂, ŵ, p̂, τ̂]⊤exp(α∗r∗+mϕ−ω∗t∗)+ cc (3.1) where û, ŵ, v̂, p̂, τ̂ are the dimensional amplitudes of the three components of velocity (in r∗, ϕ, z∗ directions), pressure and temperature, respective- Instability of the flow in rotating cavity 689 ly, α∗ and β∗ = m/r∗ are the components of wave number k∗ in the ra- dial and azimuthal directions, respectively, m is the number of spiral vorti- ces, ω∗ is the frequency and t∗ is time. Asterisk denotes dimensional values. The co-ordinate system is located on the disk under consideration. The li- near stability analysis equations plus the homogeneous boundary conditions (û(z∗) = v̂(z∗) = ŵ(z∗) = τ̂(z∗) = 0 at z∗ = 0 and z∗ = 2h) constitute an eigenvalue problem which is solved in a global manner (Tuliszka-Sznitko and Soong, 2000). 3.2. Direct numerical simulation In the DNS method, the numerical solution is based on a pseudo-spectral collocationChebyshev-Fourier-Galerkin approximation. The approximation of the flow variables Ψ = (u,w,v,p,Θ) is given by a development in truncated series (Serre et al., 2001a) ΨNMK(r,z,ϕ,t) = K/2−1∑ p=−K/2 N∑ n=0 M∑ m=0 Ψ̂nmp(t)Tn(r)Tm(z)e ipϕ −1¬ r z¬ 1 0¬ϕ¬ 2π (3.2) where N,M and K are the numbers of collocation points in the radial, axial and azimuthal directions, respectively. Tn(r) and Tm(z) are the Chebyshev polynomials. The time scheme is semi-implicit and second-order accurate. It corresponds to a combination of the second-order backward differentiation formula for the viscous diffusion terms and the Adams-Bashforth scheme for the non-linear terms. The method uses a projection scheme to maintain the incompressibility constraint. In the case of cylindrical cavity, the singularity introducedby the axis at r∗ =0, requires a dependentvariable transformation Ṽ ∗ = r∗V ∗, p̃∗ = r∗p∗, as proposed by Serre et al. (2001a). 3.3. Large eddy simulation We used the dynamic model for large eddy simulation to optimalize the subgrid-scalemode coefficient. Thismethod is used in conjunctionwith a new version of the dynamic Smagorinski eddy viscosity model proposed by Mene- veau et al. (1996). In this version, called aLagrangian subgrid-scalemodel, the Germano et al. (1991) identity was averaged for some time along fluid pathli- nes. In our incompressible computations we proceeded as follows: the filtering operationwas performed only in azimuthal direction using theGaussian func- tion. Coefficient field is computed from the following equation (Meneveau et al., 1996) C2S = LLM LMM (3.3) 690 E. Tuliszka-Sznitko, A. Zieliński where Ln+1LM (x)=H { ε[LijMij] n+1(x)+(1−ε)LnLM(x−u n∆t) } Ln+1MM(x)= ε[MijMij] n+1(x)+(1−ε)LnMM(x−u n∆t) (3.4) ε= ∆t/Tn 1+∆t/Tn and where x is the position along the pathline, ∆t is a period of time over whichaveraging is performed, n and n+1mean two consecutive time sections, u is an averaged velocity along the fluid-particle trajectory, H(x) is the ramp function (H(x) = x if x > 0, H(x) = 0 if x < 0). T is the relaxation time scale, which controls thememory of the Lagrangean averaging. Some possible choices for T were given by Meneveau et al. (1996). In our computations we define T as follows Tn = 3 2 ∆(MnijM n ij) − 1 4 (3.5) The finally obtained field of CS is averaged in azimuthal direction. 4. Linear stability theory results LSA computations were performed for the co- and counter-rotating disks and for the thermalRosby number from −0.1 to 0.1 (Tuliszka-Sznitko and Soong, 2000). For the rotor/stator case the flow consists of two disjoined boundary layers on both disks, separated by an inviscid core which rotates as a solid body. As in Itoh (1991), the solid-body angular velocity is constant and equ- al to v∗/Ω∗1r ∗ = −0.687 in the rotating frame of reference. In order to take into account the buoyancy effects, the Boussinesq approximation is invoked. In the case of non-isothermal flow, the influence of the thermal effects and the rotational-induced buoyancy on stability characteristics and the critical conditions becomes important.We have found two types of instability in both boundary layers: type II and type I. In the Bödewadt layer on the stationary disk, the onset of the type II instability has been found at ReδcII = 34.7. The type II instability exists only in a narrow range of Reδ, disappearing at Reδ =68. The type I instability occurs at a slightly larger Reynolds number, ReδcI = 47.5. The exemplary iso-lines of the temporal amplification rate ωi obtained at different local Reynolds numbers Reδ =65, 80 and 130 are shown in the plane of the wave-angle and wave-number (ε,k) in Fig.2. Our results show that the Ekman layer on the rotating disk is much more stable than the Bödewadt layer with respect to type I and type II instabilities. The on- set of type II instability has been found at ReδcII = 90.23 and type I at Instability of the flow in rotating cavity 691 ReδcI = 278.6. Our isothermal linear stability results of the Bödewadt and Ekman layers characterizing the type I and type II instabilities are in very good accordance with the results of Itoh (1991). Fig. 2. Iso-lines of ωi = const in the stationary disk boundary layer and at different local Reynolds numbers, in the plane of the wave-angle and wave-number (ε,k) (Tuliszka-Sznitko and Soong, 2000; Tuliszka-Sznitko et al., 2002) In order to analyze the effect of thermal conditions we have extended our investigations (Tuliszka-Sznitko and Soong, 2000) to the non-isothermal class of flow. Calculations have been performed for different thermal Rossby num- bers B=β(T∗2 −T ∗ 1); however, for validity of the Boussinesq approximation, the values of B have been limited to the range |B| ¬ 0.1. From the definition of the thermal Rossby number, the positive and negative values of B stand for T∗2 >T ∗ 1 (hotter stator) and T ∗ 2 0), the buoyancy driven secondary flow enforces the basic rotation driven flow. For B > 0, the fluid near the rotating disk is cooler than the fluid near the stationary disk. We have found that cooling of the stationary disk B< 0, stabilizes the flow with respect to both the type I and type II instabilities (Fig.3a) by increasing the critical Reynolds numbers. Weuse theBriggs (1964) criterionwith afixedwave number in the spanwi- se direction β to determine the region of absolute instability (Tuliszka-Sznitko et al., 2002).We have found that almost the entire layer on the stationary disk is absolutely unstable.The criticalReynoldsnumberof theabsolutely unstable flow has been found at Reδca = 48.5. On the rotating disk, the critical Rey- nolds number of the absolutely unstable flowwas determined at Reδca =562. In the next step, we extend our absolute/convective calculations to the non- isothermal class of flow (Tuliszka-Sznitko and Soong, 2000; Tuliszka-Sznitko et al., 2002). Fig.3b shows the neutral curves of absolutely unstable flow obta- ined for the stationary disk boundary layer and for different thermal Rossby 692 E. Tuliszka-Sznitko, A. Zieliński numbers (solid line). Fig.3b presents also the second families of branch points (dashed line). Fig. 3. Comparison of the neutral curves ωr = f(Reδ) in the stationary disk boundary layer obtained at different thermal Rossby number B; (a) convectively unstable areas, (b) absolutely unstable areas (dashed lines indicate second family) (Tuliszka-Sznitko and Soong, 2000; Tuliszka-Sznitko et al., 2002) 5. DNS and LES results 5.1. Cylindrical cavity In this Section we will analyze shortly the results obtained using the DNS method for cylindrical rotor/stator cavity of aspect ratio L = 5, with outer cylinder attached to the stator (Serre et al., 2004). The basic state consists of two disjoint boundary layers on each disk and of a central inviscid core flow. The fluid is pumped radially outwards, along the rotating disk (the upper one) and radially inwards, along the stator. In order to obtain an insight into the laminar-turbulent transition process in the rotor/stator flow, numerical simulations have been performed for theReynolds number close to the critical Instability of the flow in rotating cavity 693 Reynolds number of transition to unsteadiness. We gradually increased the Reynolds number and over a certain ReR, and we have observed 2D cylindri- cal vortices propagating towards the axis of cavity in the stator’s boundary layer. Cylindrical vortices are interpreted as the type II instability. Above a second critical Reynolds number, 3D spiral structures appeared (interpreted as type I instability) in the area near the outer end-wall. Spiral vortices pro- pagate towards the outer cylinder i.e. in opposite direction to the direction of the basic state. This behavior could suggest (in accordance with the LSA results) that the area of dominance of type I in the stator boundary layer may be absolutely unstable. Exemplary structures obtained at ReR = 13200 at the transient times t = 60 and t = 320 are presented in Fig.4. This pro- cedure was repeated for the flow additionally disturbed by superimposing on the initial condition at every consecutive ReR, the 3D perturbation function of the general form ηsin(pϕ), where p is an arbitrary number corresponding to an azimuthal wavelength and η is the amplitude growth rate. Calculations have been performed for η=0-3.5 and for p=2π/4.We have found that for η=0 the critical Reynolds number of transition to unsteadiness equals 12300. The amplitudes of disturbances in the stationary disk layer are much larger than the corresponding amplitudes in the rotating disk’s boundary layer. This result is in good agreement with the similarity solution presented in the pre- vious section. From the time history of the axial component of velocity in the stationary disk boundary layer (Fig.5), we can see that disturbances are first damped and the flow reaches a steady state. Then we observe the beginning of the exponential growth of disturbances. Finally, these oscillations reach an asymptotic finite-amplitude periodic state with a constant angular frequency σ=ω∗Reδ ≈ 1.1. Fig. 4. Iso-surface of fluctuations of the axial component of velocity at ReR =13200, η=0; (a) coexistence of annular and spiral structures related to Bödewadt layer instability during the transient time t=60, (b) final state showing only 12 spiral arms (Serre et al., 2004) 694 E. Tuliszka-Sznitko, A. Zieliński Fig. 5. Time histories of the axial component of velocity in the stationary disk boundary layer; ReR =13200, η=0 (Serre et al., 2004) 5.2. Annular cavity In our investigations we performed computations for the annular cavities with co- and counter-rotating disks of different aspect ratios L and different curvature parameters Rm (Tuliszka-Sznitko and Zieliński, 2006a,b). However in this sectionwe have restricted our analysis to the results obtained for L=9 and Rm =1.5 and 3.0. 5.2.1. Rotor/stator cavity Wefirst focused on the steady axisymmetric basic state of the rotor/stator case. Fig.6 shows the velocity field in the meridional section (r∗/h,z∗/h,0) obtained for the rotor/stator case (s=0,ReR =70000,B=0, the outer end- wall attached to the rotor and inner end-wall attached to the stator). From Fig.6 we can see that the flow consists of two disjoint boundary layers on each disk and of a central core flow. For arbitrary chosen positive value of s and very small negative value, themeridional structure is similar to that obtained for the rotor/stator case. Fig. 6. Flow in the meridional section, rotor/stator cavity, B=0,L=9,Rm =1.5, ReR =70000 (Tuliszka-Sznitko and Zieliński, 2006a,b) In the annular cavity we observed the same instability structures as in the cylindrical cavity. In stationary disk’s boundary layer we observed two types of vortices: 2D cylindrical vortices interpreted as type II instability and 3D spiral vortices interpreted as type I instability. Fig.7a shows the iso-lines of the azimuthal velocity component disturbances, in azimuthal section (r∗/h, z∗/h=−0.95, ϕ) obtained for L= 9, Rm =1.5, ReR =36000 and for outer Instability of the flow in rotating cavity 695 end-wall attached to the stator and the inner one attached to the rotor. For this configuration in ourDNS computations we observed 3 cylindrical vortices and 23 3D spiral vortices. The iso-lines of the axial velocity component di- sturbances in themeridional section (r∗/h, z∗/h,ϕ=0) are shown in Fig.7c. The critical Reynolds number of transition to unsteadiness was estimated at ReR =34000. Our results were compared with the experimental results obta- ined by Schouveiler et al. (1999, 2001) for the cavity of aspect ratio L=8.75, Rm =1,ReR =20900 (the same outer end-wall condition). Schouveiler (1999, 2001) observed also 3 cylindrical vortices and 18 spiral vortices (Fig.7b). Sim- lar computations, have been performed for the outer cylinder attached to the rotor and the inner one attached to the stator. For this configuration the cri- tical Reynolds number of transition to unsteadiness equaled ∼ 70000 and we observed 36 spiral vortices in the stationary disk boundary layer (Fig.7d and Fig.7f). These results were compared to the results obtained by Gauthier et al. (2002) for L = 20.9, Rm = 1 and for the outer cylinder attached to the rotor (Fig.7e). We can see that for the considered cavity (L=9, Rm =1.5), the end-wall conditions have significant influence on the critical parameters. However, this influence is expected to be negligible in the limit of a large aspect ratio L. Due to the CPU time requirements, DNS computations are limited to flows of low ormoderateReynolds numbers.For higherReynolds numbers, the computations were performed using LESmethod. Fig.8 shows the iso-lines of axial velocity componentdisturbances in themeridional section (r∗/h, z∗/h, 0) obtained for the rotor/stator case and for ReR = 36000, 50000 and 75000 (s=0,L=9,Rm =1.5,B=0, the outer end-wall attached to the stator and the inner end-wall attached to the rotor). For ReR = 50000 and 75000 the results were obtained using LES method. The iso-lines of azimuthal velocity component disturbances in the azimuthal section obtained for ReR = 75000 are presented in Fig.9. We can see that for ReR = 75000, the flow is at the last stage of the transition process. 5.2.2. Rotor/rotor cavity In this section we consider the instability patterns of the flow between two counter-rotating disks and two cylinders. In the counter-rotating case, the centrifugal flow induced by the faster disk (upper one) recirculates along the slower rotating disk (lower one) towards the inner end-wall. This inward flow meets the outward radial flow induced by the slower rotating disk, leading to a stagnation area where the radial component of the velocity vanishes. In our computations we increased slowly the rotation of the upper disk ReRupper from 6000 to 60000, with a fixed value of the bottom one ReR lower = 3000. The exemplary results are presented in Fig.10. We can see that with incre- 696 E. Tuliszka-Sznitko, A. Zieliński Fig. 7. (a) The iso-lines of azimuthal velocity component disturbances in the azimuthal section, L=9,Rm =1.5, ReR =36000, the outer end-wall attached to the stator, (b) Schouveiler et al. (2001), experimental result, Rm =1,L=8.75, (c) the iso-lines of azimuthal velocity component disturbances in the azimuthal section, z=−0.95,L=9,Rm =1.5, ReR =75000, the outer end-wall attached to the rotor, (d) Gauthier et al. (2002) experimental result, Rm =1,L=20.9, (e) the iso-lines of axial velocity component disturbances in the meridional section, L=9, Rm =1.5, ReR =36000, (f) the iso-lines of axial velocity component disturbances in the meridional section, L=9,Rm =1.5, ReR =75000 (Tuliszka-Sznitko and Zieliński, 2006a,b) asing ReRupper (with decreasing |s|), the stagnation circle moves from the outer cylinder toward the inner one. In Figs.10a,b for ReRupper = 6000 and 20000 we can see that the boundary layers on the disks are separated by the shear layer. This free shear layer, which separates two regions of opposite an- gular rotations, breaks the azimuthal symmetry of the flow. The instability structures of free shear layer are observed for much lower Reynolds numbers Instability of the flow in rotating cavity 697 Fig. 8. The iso-lines of azimuthal velocity component disturbances in the meridional section, L=9,Rm =1.5, ReR =36000, 50000, 75000 Fig. 9. The iso-lines of azimuthal velocity component disturbances in the azimuthal section z=−0.95, obtained for ReR =75000,L=9,Rm =1.5, s=0,B=0 Fig. 10. The flow and the corresponding iso-lines of axial velocity component disturbances in meridional section obtained for ReR lower =3000, and for ReRupper =6000, 20000, 600000, L=9,Rm =1.5 (Tuliszka-Sznitko and Zieliński, 2006a) 698 E. Tuliszka-Sznitko, A. Zieliński than the critical Reynolds numbers of type II and I instability in the rota- ting disks boundary layers. For some combination of parameters (ReR lower , ReRupper), interaction between the Ekman layer and the still existing she- ar layer results in negative spirals. The negative spiral vortices obtained for ReRupper =60000 arepresented inFig.11a (azimuthal section). InFig.11b the experimental results obtained by Gauthier for L = 20.9 are shown. We can see very good agreement between our numerical results and the experimental results obtained by Gauthier et al. (2002), Fig.11. Fig. 11. (a) The iso-lines of azimutal velicity component disturbances in azimuthal section obtained for: ReR lower =3000 and for ReRupper =60000, z=0.844 (Tuliszka-Sznitko and Zieliński, 2006a), (b) experimental results of Gauthier et al. (2002) 5.3. Rotor/rotor configuration with throughflow In this section we analyze the flow in a rotor/rotor cavity of s=1 (both disks rotate in the same direction with the same rotational speed) with thro- ughflow. We consider a radial outflow of fluid from a source at the inner cy- linder to a sink at the outer cylinder. Computations with throughflow require changes in the boundary conditions at the cylinders (in comparisonwith cases of impermeable cylinders considered in previous chapters). For the inner and outer cylinders we have the following boundary conditions: — for r=−1.0 u= CwL(Rm+1) 2 4πReR(Rm+r) v=w=0 — for r=1.0 u= CwL(Rm+1) 2 4πReR(Rm+r) v=w=0 where Cw = V̇ ∗/R∗1ν ∗ is adimensionless volumeflowrate. Inourcomputations we proceeded as follows: with the fixed Reynolds number we slowly increased Instability of the flow in rotating cavity 699 Cw and the result obtained for smaller Cw was used as the initial condition for higher Cw. The basic state is steady and axisymmetric. The exemplary meridional velocity field obtained for ReR = 200000 and Cw = 700 (L = 9, Rm =3) is displayed inFig.12a.We can see that for isothermal radial outflow from a uniform source at R∗0 to a uniform sink at R ∗ 1, themeridional flow can bedivided into three regions: an inner source, separatedEkman layers on each disk and an outer sink layer. In our computations, in order to accelerate the transition to unsteadiness, the axisymmetric basic flowwas perturbednear the inner cylinder by superimposing a disturbance of the following form: ηsin(pϕ) (as in Section 5.1). Results presented in Fig.12b were obtained for p = 1/4 and η = 3.4. We can see that the most unstable regions are sink and source areas. Four pairs of 2D counter-rotating rolls in both the Ekman layers are visible. Superimposed disturbance near the inner end-wall produced the rise of weak 8 spiral vortices in this area. Fig. 12. The mean flow (a) and the iso-lines of axial velocity component disturbances in meridional section obtained for ReR =200000,Cw =700,B=0, L=9,Rm =3 (Tuliszka-Sznitko and Zieliński, 2006a) 5.4. Non-isothermal flow in annular cavity The heat transfer in flow between two rotating disks is a problem of great importance to the gas-turbine air cooling designer. In this chapter we present our exemplary results obtained for the non-isothermal flow to estimate the influence of the thermal Rossby number on the basic state, distribution of the local Nusselt numbers and on the instability structures. 5.4.1. Rotor/stator configuration The influence of the thermal Rossby number on the basic state is small. However, we have found a significant influence of B on the critical Reynolds number of transition to unsteadiness. As we expected, the critical Reynolds number of transition to unsteadiness in the stationary disk’s boundary layer increases with decreasing B. The instability structures obtained for the non- isothermal flow are roughly the same as those for the incompressible cases. In 700 E. Tuliszka-Sznitko, A. Zieliński the rotor/stator cases we observed 2D cylindrical vortices for the lower Rey- nolds number and positive 3D spiral vortices for the higher Reynolds number, however we have found that the number of spiral vortices decreases with in- creasing B. To verify our non-isothermal results we compared distribution of the local Nusselt numbers on the rotating disk (L= 9, Rm = 1.5, outer end wall attached to the stator and inner one to the rotor) with the experimental data ofMochizuki andYangWen-Jei (1986) andwith the results obtained for a single heated rotating disk using the similarity analysis. Agreement of the results is good. Exemplary distributions of Nu along the radius r∗/h obta- ined fromDNS computations and from the similarity equations for the single heated rotating disk for B =−0.1, ReR =25000, 35000, 47000 are presented in Fig.13. The local Nusselt number is calculated from the following equation Nu= α∗r∗ λ∗ =− (∂Θ ∂z ) w r∗ h 1 Θw−Θc (5.1) where ΘC is the dimensionless temperature of inviscid core, λ ∗ is the coeffi- cient of thermal conductivity, α∗ is the coefficient of heat transfer and index w a indicates parameter at the wall. In Fig.13 the areas near the end-walls we- re cut off because of difficulties with definition of temperature ΘC near the cylinders. From Fig.13 we can see that the DNS and the similarity solutions (straight broken lines in Fig.13) coincide well along a large part of the rotor. Fig. 13. The exemplary distributions of the local Nusselt number on the rotor obtained for B=−0.1, ReR =25000, 35000, 47000 (Tuliszka-Sznitko and Zieliński, 2006a) 5.4.2. Rotor/rotor configuration with throughflow Themost interesting from the point of view of turbomachinery air cooling devices are examples with the superimposed radial outflow. In Fig.14 the iso- lines of temperature anddisturbances of theaxial velocity componentobtained for L=9, Rm =3, ReR =200000, B =0.1, t=16 and Cw from 300 to 700 Instability of the flow in rotating cavity 701 are displayed. Both disks rotate in the same direction with the same angular speed and on both disks Θ = 1. The boundary conditions at the inner and outer cylinders are as follows: — for r=−1.0 u= CwL(Rm+1) 2 4πReR(Rm+r) v=w=0 Θ=0 — for r=1.0 u= CwL(Rm+1) 2 4πReR(Rm+r) v=w=0 Θ=1 Figures 14 show the effectiveness of radial cooling; the areas dominated by coming cooling air are laminar for all considered mass flow rates. Fig. 14. The iso-lines of temperature and disturbances of the axial velocity component obtained for L=9,Rm =3, ReR =200000,B=0.1, t=16 and Cw =300 and 700 6. Conclusions In the present paper results of numerical simulation of the transitional flow with heat transfer in a rotating cavity were reviewed.The isothermal andnon- isothermal 3D fluid flows between co- and counter-rotating disks enclosed by two rotating cylinders, were investigated. Computations have been performed using theDNS, LES andLSAmethods. Special attention has been paid to the basic laminar flow in order to obtainmore insight into the onset of thedifferent instability patterns and their regions of existence. Three different patterns 702 E. Tuliszka-Sznitko, A. Zieliński have been found: axisymmetrically propagating vortices interpreted as type II instability, positive spiral vortices interpreted as type I instability andnegative spirals. The first two vortices, i.e. cylindrical vortices and positive spirals were present in all the considered configurations; cylindrical and annular cavities of the different aspect ratios and curvature parameters. Negative spirals, which exist in the counter-rotating configurations, differ significantly from the 2D and positive 3D vortices. First of all, the negative spirals rotate around the central axis and the mechanism is not of a cross-flow type; it is probably the result of interaction between free shear layer and the Ekman layer. We have found that the negative spirals can co-exist with the positive vortices. These results were discussed using the experimental results ofGauthier et al. (2002). Comparison shows a good agreement. Our LSA computations also turned out to be very useful in interpretation of 2D and positive 3D vortices in terms of the type II and type I instability. Additionally, the LSA allowed us to identify the areas of absolutely unstable flow. The configurations which are most interesting from the point of view of possible applications are those with the superimposed flow (flow with the source and sink). For the incompressible configuration we have found similar instability structures as those in paper by Serre et al. (2001b). We extended these computations to the non-isothermal cases. Our non-isothermal results showed that influence of the thermal Rossby number on instability structure is not large, however we have found significant influence B on the critical Reynolds number of transition to unsteadiness. Distributions of the local Nusselt numbers along the radius of the disk show a good agreement with the experimental data of Mochizuki and Yang Wen-Jei (1986) and with theoretical solutions. Computations showed effectiveness of the radial cooling. Acknowledgements The computations were carried out in the Poznań Supercomputing and Networ- king Center and in the frame of grants 4T07A04528 and COST/78/2006. References 1. BohnD., EmundsR., GorzelitzV., KrligerU., 1994,Experimental and theoretical investigations of heat transfer in closed gas filled annuli II, ASME Int. Gas. Turbine and Aeroengine Cong., The Hague, Paper 94-GT-175 2. Briggs R., 1964,Electron-Stream Interaction with Plasmas, MIT Press 3. Chew J., Owen J., Pincombe J., 1984, Numerical prediction for laminar source-sink flow in a rotating cylindrical cavity, J. Fluid Mech., 143, 451-466 Instability of the flow in rotating cavity 703 4. DijkstraD., and vanHeijstG., 1983,The flowbetween twofinite rotating discs enclosed by a cylinder, J. Fluid Mech., 128, 123, 123-154 5. Gauthier G., Gondret P., Moisy F., Rabaut M., 2002, Instabilities in the flow between co- and counter-rotating disks, J. Fluid Mech., 473, 1-21 6. GermanoM.,PiomelliU.,MoinP.,CabotW., 1991,Adynamic sub-grid eddy viscositymodel,Phys. of Fluids A, 3, 1760-1765 7. Itoh M., 1991, On the instability of the flow between coaxial rotating disks, In: Boundary Layer Stability and Transition to Turbulence, ASME FED 114, 83-89 8. Lingwood R., 1997, Absolute instabilities of the Ekman layer and related rotating flows, J. Fluid. Mech., 331, 405-428 9. Meneveau C., Lund T., Cabot W., 1996, A Lagrangian dynamic subgrid- scale model of turbulence, J. Fluid Mech., 319, 353-385 10. Mochizuki S., Yang Wen-Jei, 1986, Flow friction on co-rotating parallel- disk assemblies with forced through-flow,Experiments in Fluids, 4, 56-60 11. Moisy F., Doare O., Pasutto T., Daube O., Rabaud M., 2004, Experi- mental and numerical study of the shear layer instability between two counter- rotating disks, J. Fluid Mech., 507, 175-203 12. Owen J., Rogers R., 1995,Heat Transfer in Rotating Disc Systems, Vol. 2: Rotating Cavity, W.D.Morris (Edit.), Wiley 13. Schouveiler L., Le Gal P., Chauve P., Takeda Y., 1999, Spiral and circularwaves in the flowbetweena rotating anda stationarydisc,Experiments in Fluids, 26, 179 14. Schouveiler L., Le Gal P., Chauve P., 2001, Instabilities of the flow between a rotating and stationary disk, J. Fluid Mech., 443, 329-350 15. Serre E., Crespo del Arco E., Bontoux J., 2001a, Annular and spiral patters in flowbetween rotating and stationary disks, J. Fluid. Mech., 434, 65, 65-100 16. Serre E., Hugues S., Crespo del Arco E., Randriamampianina A., Bontoux P., 2001b, Axisymmetric and three-dimensional instability in an Ekman boundary layer flow, Int. J. Heat Fluid Flow, 22, 82-93 17. Serre E., Tuliszka-Sznitko E., Bontoux P., 2004, Coupled numerical and theoretical study of the flow transition between a rotating and stationary disk,Phys. of Fluids, 16, 3 18. Tuliszka-Sznitko E., Soong C.Y., 2000, Instability of non-isothermal flow between coaxial rotating disks, Proceedings of European Congress on Compu- tational Methods in Applied Sciences and Engineering, Barcelona 19. Tuliszka-Sznitko E., Serre E., Bontoux P., 2002, On the nature of the boundary layers instabilities in a flowbetween a rotating and a stationary disc, C.R. Mecanique, 30, 90-99 704 E. Tuliszka-Sznitko, A. Zieliński 20. Tuliszka-SznitkoE., Zieliński A., 2006a, Instability of the isothermal and not-isothermal flow in rotating cavity of aspect ratio 9,Proceedings of ISTP17, Japan 21. Tuliszka-Sznitko E., Zieliński A., 2006b, Numerical investigation of the transitional flow in co- and counter-rotating annular cavity, Journal of Theore- tical and Applied Mechanics, 44, 2, 405-421 Niestabilność przepływów w wirujących przestrzeniach Streszczenie Wpracy dokonano analizywynikówdotychczasowychpracwłasnychnad statecz- nością przepływów w wirujących przestrzeniach. Badania te prowadzono metodą li- niową,metodąDNS i LES.Do trójwymiarowychobliczeń numerycznych zastosowano metodę spektralnej kolokacji bazującąna szeregachCzebyszewa i szereguFouriera.Do metody LES zastosowanometodę Lagrangea zaproponowanąprzezMeneveau. Głów- nym celem było zbadanie struktur niestabilnościowych występujących w wirujących przestrzeniach.Wyniki te porównano zwynikami prac eksperymentalnych innych au- torów. Badanowpływ termicznej liczbyRossbiego na rozkłady liczbNusselta wzdłuż promienia tarcz oraz efektywność chłodzenia promieniowego tarcz za pomocą napły- wającego chłodnego czynnika. Manuscript received January 3, 2007; accepted for print March 19, 2007