Numerical simulation of turbulent flow between shrouded contra-rotating disks

Physical model

The schematic of coordinate and physical system adopted in the present simulation is shown in Figure 1, and the computational geometry matches with that of the experimental study carried out by Kilic et al. ¹⁶ and Gan and MacGregor. ¹⁷ The rotating disk cavity system is formed by two disks (disk 1 and disk 2 separated by a spacing of s), with cylindrical shrouds attached to each disk at the outer disk radii. A small axial clearance, s_c, separates the shrouds at the axial midplane allowing each disk/shroud set to rotate freely from the other. Disk 1 rotates clockwise at speed Ω₁, and disk 2 rotates counterclockwise at speed Ω₂. The rotation ratio of disk 2 to disk 1 is denoted by Γ(Γ = Ω₂/Ω₁), which is varied in a range of −1 ≤ Γ ≤ 0. In order to present the rotational effects on the turbulent flow between the shrouded contra-rotating disks, the rotational Reynolds number is defined as Re_ϕ = Ω₁b²/v, where v is the kinematic viscosity of cooling air.

OPEN IN VIEWER

Figure 1.

Schematic of coordinate and physical system.

Governing equations and turbulence models

General formulation

The flow in the shrouded contra-rotating disk cavity is assumed to be steady rotating axisymmetric flow, and it is also assumed that the cooling air is an incompressible fluid with constant properties. The experimental results in the work of Gan and MacGregor ¹⁷ show that the flow in the core is turbulent for Re_ϕ as low as 2.16 × 10⁴. For the parameter range of 6 × 10⁴ < Re_ϕ < 4 × 10⁶ used in the calculation, the flow in the core region inside contra-rotating disk cavity is turbulent. The Reynolds stresses are related to the mean strain rate via an isotropic scalar turbulent viscosity, v_t

V ′ z V ′ r ¯ = − v t ( ∂ V z ∂ r + ∂ V r ∂ z )

(1)

The steady-state Reynolds averaged mass and momentum equations with respect to a reference frame at a rotating velocity Ω can be depicted as

∇ · V = 0

(2)

ρ ( V ∇ ) V = − ∇ P + ( μ + μ t ) ∇ 2 V − ρ Ω × ( Ω × r ) − 2 ρ Ω × V

(3)

where V is the velocity vector (V_r, V_ϕ, V_z), P is pressure, ρ is density, µ is the dynamic viscosity, − ρ Ω × ( Ω × r ) is the centrifugal force, and − 2 ρ Ω × V is the Coriolis force.

Two k-ε turbulence models were employed in the present simulation, which are the TL turbulence model that includes the standard k-ε turbulence model for the flow in the core region, a single-equation model for the flow near walls, and the modified low-Reynolds number k-ε model of Launder–Sharma ¹⁵ that is able to model the transition from the turbulent core to viscous sublayer close to the wall without the use of special near-wall models.

The TL turbulence model

For the core region flow, the standard k-ε turbulence model is used to calculate the turbulent viscosity µ_t as follows

∇ · ( ρ V k ) = ∇ · [ ( μ + μ t σ k ) ∇ k ] + ρ G k − ρ ε

(4)

∇ · ( ρ V ε ) = ∇ · [ ( μ + μ t σ ε ) ∇ ε ] + ε k [ c ε 1 G k − ρ c ε 2 ε ]

(5)

μ t = c μ ρ k 2 ε

(6)

where G_k represents the generation of turbulence kinetic energy due to the mean velocity gradients calculated as

G k = ν t [ 2 ( ∂ V z ∂ z ) 2 + 2 ( ∂ V r ∂ r ) 2 + 2 ( V r r ) 2 + ( ∂ V z ∂ r + ∂ V r ∂ z ) 2 + ( ∂ V ϕ ∂ z ) 2 + ( ∂ V ϕ ∂ r − V ϕ r ) 2 ]

(7)

All the constants c_ε1, c_ε2, c_µ, σ_k, and σ_ε in governing equations have the generally agreed values ¹⁸

c ε 1 = 1 . 44 , c ε 2 = 1 . 92 , c μ = 0 . 09 , σ k = 1 . 0 , σ ε = 1 . 3

For the near-wall regions flow, Wolfshtein’s single-equation model ¹⁹ is used to account for the wall proximity effects. In Wolfshtein’s approach, Equation 4 is used to compute k near the wall, and ε and µ_t are obtained as

ε = k 3 / 2 l ε , μ t = ρ c μ k l μ ,

(8)

by estimating the length scales l_ε and l_µ from

l ε = 2 . 55 y [ 1 − exp ( − 0 . 263 y * ) ]

(9)

l μ = 2 . 55 y [ 1 − exp ( − 0 . 016 y * ) ]

(10)

y * = y k 1 / 2 v

(11)

respectively, where y denotes the distance of the closest grid point to the wall. At the wall, k = 0 is set, and at the shroud clearance, ∂ k / ∂ z = ∂ ε / ∂ z = 0 is set.

Modified Low-Reynolds turbulence model of Launder–Sharma

With respect to the performance of the Launder–Sharma turbulence model, there is a tendency for it to predict excessively high levels of near-wall turbulence. ²⁰ The incorporation of the so-called Yap correction for the turbulence length scale and a rotating-related modification for the dissipation rate has led to improved predictions of the velocity distribution for a rotor-stator system, ²¹ but as far as the authors are aware, these modifications have not been applied to contra-rotating disks.

In the Launder–Sharma turbulence model with the incorporation of the Yap correction (LSY), the turbulent viscosity, µ_t, is obtained from

μ t = c μ f μ ρ k 2 ε ~

(12)

∇ · ( ρ V k ) = ∇ · [ ( μ + μ t σ k ) ∇ k ] + ρ G k − ρ ( ε ~ + D )

(13)

∇ · ( ρ V ε ~ ) = ∇ · [ ( μ + μ t σ ε ) ∇ ε ~ ] + ε ~ k [ c ε 1 G k − ρ c ε 2 f 2 ε ~ ] + ρ E + Y a p

(14)

where

f 2 = 1 − 0 . 3 exp ( − R ~ e t 2 )

(15)

f μ = exp [ − 3 . 4 ( 1 + 0 . 02 R ~ e t ) 2 ]

(16)

R ~ e t = ρ k 2 ε ~ μ

(17)

D = 2 v [ ( ∂ k 1 / 2 ∂ x i ) 2 ]

(18)

E = 2 v v t [ ( ∂ 2 u j ∂ x i ∂ x k ) 2 ]

(19)

Yap = max [ 0 . 83 ρ ε ~ 2 k ( k 1 . 5 2 . 55 y ε ~ − 1 ) ( k 1 . 5 2 . 55 y ε ~ ) 2 , 0 ]

(20)

At the wall, ∂ k / ∂ z = ε ~ = 0 is set, and at the shroud clearance, ∂ k / ∂ z = ∂ ε ~ / ∂ z = 0 is set.

Boundary conditions

The no-slip boundary conditions on velocity are applied on all the solid walls. The axial and radial velocities are taken to be zero, and the azimuthal velocity is assumed to vary linearly at the shroud clearance

z = 0 , 0 ≤ r ≤ b : V r = 0 , V z = 0 , V ϕ = Ω 1 r

(21)

z = s , 0 ≤ r ≤ b : V r = 0 , V z = 0 , V ϕ = Ω 2 r

(22)

0 ≤ z ≤ ( s − s c ) / 2 , r = b : V r = 0 , V z = 0 , V ϕ = Ω 1 b

(23)

( s + s c ) / 2 ≤ z ≤ s , r = b : V r = 0 , V z = 0 , V ϕ = Ω 2 b

(24)

( s − s c ) / 2 < z < ( s + s c ) / 2 , r = b : V r = 0 , V z = 0 , V ϕ = 1 2 Ω 1 b ( ( 1 + Γ ) + ( 1 − Γ ) s − 2 z s c )

(25)

Flow structures at various rotation ratio of the two disks

Figure 2(a) and (b), respectively, presents the contours of stream function Ψ computed with LSY model and TL model for Re_ϕ = 1.25 × 10⁶ and Γ = −0.6. In the Ekman layer attached on each of the two rotating disks, the fluid is pumped radially outward by the combination of the centrifugal force and the radial component of Coriolis force. The outward radial flow along the faster rotating disk 1 is deflected by the shroud, then rolls radially inward along the slower rotating disk 2, pushes against the outward radial flow along disk 2 because of the unbalance radial forces, and then the fluid is forced downward into the core region as well as axially toward the two disks separately. Consequently, a larger clockwise circulation and a weaker counterclockwise circulation separated by the stagnation streamline originating from the stagnation point (at a nondimensional radius x_st = r_st/b near the slower rotating disk) are formed in the r-z plane. Both turbulence models can describe the flow structure except the location of the stagnation point (to be discussed in more detail in the section “The critical radius of flow structure transition for −1 < Γ < 0”).

OPEN IN VIEWER

Figure 2.

Contours of stream function modeled using the (a) LSY and (b) TL models for Re_ϕ = 1.25 × 10⁶ and Γ = −0.6.

Figure 3 shows the axial profiles of V_r/Ωr and V_ϕ/Ωr corresponding to the conditions in Figure 2. From Figure 3(a) and (b), it can be seen that radial inward flow occurs across the entire interior core, and the tangential velocity varies almost linearly with the axial position (z/s) in the core region at r/b = 0.6. The flow structure is of Stewartson type. ²³ At r/b = 0.85, The region near the axial midplane (0.3 < z/s < 0.7) is characterized by a quasi zero radial velocity and a constant tangential velocity. The angular speed of the rotating core increases slightly with increasing radial distance (V_ϕ/Ωr = 0.2 for r/b = 0.85 and V_ϕ/Ωr = 0.25 for r/b = 0.9). This is consistent with the Batchelor-type flow characteristic. ²⁴ In fact, the Stewartson-type flow holds near the center region, while the Batchelor-type flow is limited to the periphery region of the disk cavity for −1 < Γ < 0. The flow structure transition is associated with the two-cell structure as described above, that is, Stewartson-type flow occurs radially inward of the stagnation point and Batchelor-type flow radially outward. Both turbulence models yield qualitative agreement with Kilic’s experiment predictions ¹⁶ for this transition.

OPEN IN VIEWER

Figure 3.

Axial profiles of V_r/Ωr and V_ϕ/Ωr for Re_ϕ = 1.25 × 10⁶, Γ = −0.6: (a) V_r/Ωr and (b) V_ϕ/Ωr.

Figure 4(a) and (b), respectively, present the contours of stream function Ψ computed with LSY model and TL model for Re_ϕ = 3.93 × 10⁵ and Γ = −1. Owing to the centrifugal force and the radial component of Coriolis force, the air near contra-rotating disk walls flows radially outward to the shrouds and turns inward at the shrouds clearance, then meets to form an inflowing free shear layer in the midplane (z/s = 1/2), and ultimately flows axially to the adjacent disk from the free shear layer. The mean secondary flow in the r-z plane thus consists of two large counter recirculation cells, with the centers of circulating cells occurring near the shrouds (x = r/b = 0.95 for LSY model and x = r/b = 0.935 for TL model).

OPEN IN VIEWER

Figure 4.

Contours of stream function modeled using the (a) LSY and (b) TL models for Re_ϕ = 3.93 × 10⁵ and Γ = −1.

Figure 5 shows the axial profiles of V_r/Ωr and V_ϕ/Ωr corresponding to the conditions in Figure 4. From Figure 5(a), it is seen that there is no radial inward “jet” at the midplane same as Batchelor-type flow. Instead, radial inward flow occurs over the entire region between the two boundary layers. From Figure 5(b), it is seen that the tangential velocity in the core region varies almost linearly with the axial position, z/s, and there is no rotating core between the disk boundary layer and the free shear layer as Batchelor ²⁵ predicted. Obviously, the flow between the contra-rotating disks at the same speeds is of a Stewartson type. ²³ Both the predictions of the two turbulence models agree with the Stewartson-type flow structure as obtained in Gan’s experiment. ¹⁷

OPEN IN VIEWER

Figure 5.

Axial profiles of V_r/Ωr and V_ϕ/Ωr for Re_ϕ = 3.93 × 10⁵, Γ = −1: (a) V_r/Ωr and (b) V_ϕ/Ωr.

The critical radius of flow structure transition for −1 < Γ < 0

As discussed above, Stewartson-type flow occurs radially inward of the stagnation point and Batchelor-type flow radially outward for −1 < Γ < 0. The location of the stagnation point can be obtained by inspecting the sign change of the radial velocity component near the slower rotating disk 2. Figure 6 shows the axial distribution of V_r/Ωr at x = r/b = 0.7 and x = r/b = 0.75 calculated with the LSY model for Re_ϕ = 1.25 × 10⁶ and Γ = −0.6. The V_r/Ωr at x = 0.7 is positive near disk 2, indicating that the flow there is radially outward. The radial velocity profile at x = 0.75 is negative near disk 2, indicating that the flow is radially inward at that location. Thus, the stagnation point occurs between 0.7 < x_st < 0.75 for the LSY turbulence model. Figure 7 shows the axial distribution of V_r/Ωr at x = r/b = 0.8 and x = r/b = 0.85 calculated with TL model corresponding to the conditions in Figure 6. The V_r/Ωr at x = 0.8 is positive near disk 2, and the radial velocity profile at x = 0.85 is negative near disk 2. Thus, the stagnation point occurs between 0.8 < x_st < 0.85 for the TL turbulence model. Calculations for the locations of the stagnation points were made at other values of Γ in addition to those shown here for Re_ϕ = 1.25 × 10⁶, and the results are shown in Table 1, accompanied by the experimental data from Kilic et al. ¹⁶ Both the computations and measurements predict the fact that the nondimensional critical radius of flow structure transition x_st increases as the rotation ratio Γ decreases. However, the numerical results from LSY turbulence model match the experimental data much better than that for TL turbulence model, although the location of the stagnation point for Γ = −0.4 could not be confirmed by Kilic because no measurements were made for x < 0.6 in their experiments. Obviously, TL turbulence model predicts larger values of x_st, that is, TL turbulence model predicts larger domains of Stewartson-type flow.

OPEN IN VIEWER

Figure 6.

Axial profiles of V_r/Ωr at x = 0.7 and x = 0.75 modeled using the LSY model, Re_ϕ = 1.25 × 10⁶, Γ = −0.6.

OPEN IN VIEWER

Figure 7.

Axial profiles of V_r/Ωr at x = 0.8 and x = 0.85 modeled using the TL model, Re_ϕ = 1.25 × 10⁶, Γ = −0.6.

OPEN IN VIEWER

Table 1.

The locations of the stagnation points.

	Γ = −0.4	Γ = −0.5	Γ = −0.6	Γ = −0.8
Experimental data 16	xst < 0.6	xst = 0.6	xst = 0.75	xst = 0.85
Results of LSY	0.45 < xst < 0.5	0.55 < xst < 0.6	0.7 < xst < 0.75	0.85 < xst < 0.9
Results of TL turbulence model	0.55 < xst < 0.6	0.65 < xst < 0.7	0.8 < xst < 0.85	0.9 < xst < 0.95

LSY: Launder–Sharma turbulence model with the incorporation of the Yap correction; TL: two-layer.

The quantitative prediction of V_r/Ωr and V _ϕ /Ωr

From Figures 3 and 5 it can be seen that at all values of x = r/b considered in this study, there is good agreement between the computed and measured velocities for both the radial and tangential components of velocity for LSY turbulence model. However, for TL turbulence model, the peak value of V_r/Ωr in the boundary layers is higher than that measured, and the tangential velocity in the rotating core is lower than that measured even though it is well predicted in the Stewartson-type flow region. The computed axial distributions of V_r/Ωr and V_ϕ/Ωr for Γ = −0.2, −0.4, and −0.8 shown in Figures 8–10, respectively, support the result. When the flow becomes more turbulent, the peak value of V_r/Ωr is increased as a result of stronger momentum exchange between the boundary layer and the core flow, and V_ϕ/Ωr in the rotating core is decreased as a result of energy dissipation. This indicates that the viscous effects are slightly underestimated by the TL turbulence model because the underestimated viscous effects cause the flow to become fully turbulent earlier than the actual experimental results.

OPEN IN VIEWER

Figure 8.

Computed axial profiles of V_r/Ωr and V_ϕ/Ωr for Re_ϕ = 1.25 × 10⁶, Γ = −0.2: (a) V_r/Ωr and (b) V_ϕ/Ωr.

OPEN IN VIEWER

Figure 9.

Computed axial profiles of V_r/Ωr and V_ϕ/Ωr for Re_ϕ = 1.25 × 10⁶, Γ = −0.4: (a) V_r/Ωr and (b) V_ϕ/Ωr.

OPEN IN VIEWER

Figure 10.

Computed axial profiles of V_r/Ωr and V_ϕ/Ωr for Re_ϕ = 1.25 × 10⁶, Γ = −0.8. (a) V_r/Ωr and (b) V_ϕ/Ωr.

The quantitative prediction of turbulent quantities

The turbulent viscosity fields predicted by the two turbulence models are shown in Figures 11 and 12. The maximum turbulent viscosity is located at the center at about x = r/b = 0.7 for Γ = −0.6 and moves radially outward for Γ = −1. The locations of the minimum turbulent viscosity are at the wall for all conditions. The distributions of the turbulent viscosity predicted by the LSY and TL models are similar; however, the maximum value of the turbulent viscosity predicted by TL model is larger than that predicted by the LSY model by 72% and 60% for Γ = −0.6 and Γ = −1, respectively. The axial distributions of turbulent viscosity µ_t at x = r/b = 0.6 for different conditions are compared in Figure 13 for the two turbulence models. It is also clear that the TL model predicts much larger values of µ_t both in the boundary layer and in the core region than the LSY model. Therefore, it is verified that TL turbulence model predicts a more turbulent flow inside contra-rotating disk cavity compared to the LSY model. In conjunction with the forgoing result that TL model predicts a larger Stewartson-type flow region compared to the LSY model and experimental measurements, it can be concluded that a more turbulent flow in the core has destroyed the Batchelor-type flow and creates a larger Stewartson-type flow region.

OPEN IN VIEWER

Figure 11.

Contours of µ_t computed with the (a) LSY and (b) TL models for Re_ϕ = 1.25 × 10⁶ and Γ = −0.6.

OPEN IN VIEWER

Figure 12.

Contours of µ_t modeled using the (a) LSY and (b) TL models for Re_ϕ = 3.93 × 10⁵ and Γ = −1.

OPEN IN VIEWER

Figure 13.

Axial distributions of µ_t at x = r/b = 0.6 for (a) Re_ϕ = 1.25 × 10⁶, Γ = −0.6 and (b) Re_ϕ = 3.93 × 10⁵, Γ = −1.

The mechanisms of flow structure transition

In equation (3), with V  = V_re_r + V_ϕe_ϕ + V_ze_z, Ω = Ωe_z, and e_r, e_ϕ, and e_z denoting unit vectors in radial, tangential, and axial directions, respectively, the centrifugal force and the Coriolis force can be expressed, respectively, as follows

− ρ Ω × ( Ω × r ) = ρ Ω 2 r e r

(26)

− 2 ρ Ω × V = 2 ρ Ω V ϕ e r − 2 ρ Ω V r e ϕ

(27)

from which one can see the centrifugal force lies in radially outward direction and the Coriolis force has one radial component 2ρΩV_ϕe_r and one tangential component −2ρΩV_re_ϕ. With the same sign of Ω and V_ϕ, the radial component 2ρΩV_ϕe_r is positive and radially outward.

For the case of Γ = −1, the radial velocities are positive in the Ekman layer on the counter-rotating disks due to the radially outward centrifugal force and Coriolis force, and are negative in the core region between the two boundary layers due to the deflecting effects of the shrouds. So the tangential viscous forces near the two disks are in opposite direction with Ω₁ = −Ω₂. The tangential component of Coriolis force is also of opposite directions both in the boundary layers on the two disks and in the core region with Ω₁ = −Ω₂ and V_r1 = V_r2. As a result, the inward radial flow penetrates the core region and the opposite rotating flows on either hand side of the axial midplane are neutralized, which is the forming mechanism of the Stewartson-type flow for the case of Γ = −1.

As illustrated schematically in Figure 14, the tangential component of Coriolis force on the fluid on left hand side of the axial midplane can be calculated as −2ρΩ₁V_r1, which is positive with positive Ω₁ and negative V_r1, while the tangential Coriolis force on the fluid on right hand side can be calculated as −2ρΩ₂V_r2, which is negative with negative Ω₂ and V_r2. So the tangential Coriolis forces on both sides of the axial midplane are of the same magnitude but in opposite direction. The coupling of the opposite tangential Coriolis forces generates tangential velocity gradient in axial direction and tears fluid circumferentially. As a result, the flows between the counter-rotating disks appear more chaotic and tumbling, which provided the interpretation for Gan and MacGregor’s ¹⁷ experiment result that the flow in counter-rotating disk cavity transits earlier to turbulent compared with rotor-stator disk cavity.

OPEN IN VIEWER

Figure 14.

Sketches of the coupling effects of the opposite tangential Coriolis forces for Γ = −1: (a) left hand of the axial midplane and (b) right hand of the axial midplane.

For the case of −1 < Γ < 0, the fluid near the faster rotating disk 1 moves toward disk 2 as well as radially inward. As a consequence, the radial velocity is quasi-zero in the core region and the fluid is brought on rotation with disk 1 like a sandwich via the combined effects of tangential viscous force near disk 1 and tangential components of the Coriolis force near disk 2 for x ≥ x_st. This is known as the Batchelor flow structure. With the sandwich-like rotating core acting in the same way as the rotating shroud, the forming mechanism of the Stewartson-type flow structure follows the same philosophy as that for Γ = −1 for x ≤ x_st.