In the following analysis, the SST PANS model (Luo et al., 2014; Ranjan and Dewan, 2015; Pereira et al., 2015; Ranjan and Dewan, 2016; Pereira et al., 2018; Qian et al., 2020) is used. The transport equations of the SST PANS model are as follows: ∂ k u ∂ t + u j ∂ k u ∂ x j = P k u ∼ − β ∗ k u ω u + ∂ ∂ x j ( ( ν + σ k u ν u ) ∂ k u ∂ x j ) , (1) ∂ ω u ∂ t + u j ∂ ω u ∂ x j = γ ν u P k u ∼ − γ ν u β ∗ ( 1 - 1 f ω ) k u ω u − β ω u 2 f ω + 2 ( 1 − F 1 u ) σ ω 2 1 ω u ∂ k u ∂ x j ∂ ω u ∂ x j f ω f k + ∂ ∂ x j ( ( ν + σ ω u ν u ) ∂ ω u ∂ x j ) , (2) ν u = a 1 k u max ( a 1 ω u , F 2 u S ) , (3) where u is the partially averaged velocity; partial averaging corresponds to filtering a portion of the fluctuating scales (Girimaji, 2006), and throughout the study, the words filtering and averaging will be used synonymously. k _u is the unresolved turbulent kinetic energy, ω _u is the unresolved specific dissipation rate, f _ω is the modeled-to-total ratio of the specific turbulence dissipation and f ω = 1 / f k , υ is the kinematic viscosity, ν u is the unresolved eddy viscosity, S is the invariant measure of the strain rate, P ∼ k u = min ( P k u , 10 β ∗ k u ω u ) is the production term, P k u = ν u ∂ u i ∂ x j ( ∂ u i ∂ x j + ∂ u j ∂ x i ) , and a 1 , β, β ^* , σ_ω2, and γ are constant coefficients. Moreover, σ k u and σ ω u are Prandtl numbers, and their specific forms are as follows: σ k u = σ k f ω f k , σ ω u = σ ω f ω f k . (4)

In addition, F _1u and F _2u are two blending functions for the SST PANS model, which are defined as follows: F 1 u = tanh [ ( min [ max ( k u β ∗ ω u y , 500 ν y 2 ω u ) , 4 σ ω 2 u k u C D k ω y 2 ] ) 4 ] , (5) F 2 u = tanh [ ( max ( 2 k u β * ω u y ; 500 ν y 2 ω u ) ) 2 ] , (6) C D k ω = max ( 2 ρ σ ω 2 u ω u ∂ k u ∂ x i ∂ ω u ∂ x i , 10 − 10 ) , (7) where y denotes the distance to the next surface.

According to the Kolmogorov hypothesis (Pope, 2000), the energy spectrum without rotation is E ( κ ) = α ε 2 / 3 κ − 5 / 3 , where E ( κ ) is the energy spectrum, α is the Kolmogorov constant, ε is the dispassion of turbulent kinetic energy, and κ is the wavenumber. Several researchers derived different f _k equations (Song and Park, 2009; Foroutan and Yavuzkurt, 2014; Qian et al., 2020) based on the Kolmogorov hypothesis. The f _k expression deduced by Song and Park (2009) can be seen as follows: f k = k u k t = ∫ κ Δ κ η E ( κ ) d κ ∫ κ l κ η E ( κ ) d κ = κ η − 2 / 3 − κ Δ − 2 / 3 κ η − 2 / 3 − κ l − 2 / 3 ≈ ( Δ η ) 2 / 3 ( 1 − ( Δ η ) 2 / 3 ) ( Δ l t u r b ) 2 / 3 . (8)

Eq. 8 is named PANS-ES1. It is similar to that proposed by Song and Park (2009), and l _turb , η, and Δ denote the turbulent length scale for RANS, Kolmogorov length scale (length scale for DNS), and grid size, respectively. The corresponding coupled model with SST k-ω PANS is called SST k-ω PANS-ES1 hereafter. l t u r b = k T 3 / 2 ε , η = ( ν 3 ε ) 1 / 4 , Δ = ( Δ x Δ y Δ z ) 1 / 3 , (9) where k _T is the total turbulent kinetic energy, which is equal to k _r + k _u . In addition, k _r is the resolved turbulent kinetic energy, k r = 0.5 ( U i − U − ) 2 , where U i is the instantaneous velocity and U − is the time-averaged velocity.

Different from Song et al., Foroutan and Yavuzkurt (2014) adopted a von Kármán-like spectrum (Schiestel and Dejoan, 2005) and derived another f _k formulation, (Eq. 10), which here is referred to as PANS-ES2, and the corresponding coupled model with SST k-ω PANS is called SST k-ω PANS-ES2 hereafter. f k = 1 − [ ( l t u r b / Δ ) 2 / 3 0.23 + ( l t u r b / Δ ) 2 / 3 ] 4.5 (10)

For the rotation-corrected energy spectrum, the energy spectrum between κ _l and κ _η is divided into two parts by the Zeman number κ Ω (Zeman, 1994), which is defined as ( Ω 3 / ε ) 1 / 2 , where Ω denotes the rotation speed. More specifically, when κ c < κ Ω , E ( κ ) = β ( Ω ε ) 1 / 2 κ − 2 , while when κ c > κ Ω , E ( κ ) = α ε 2 / 3 κ − 5 / 3 (Figure 1).

From the definition of f _k , f k = k u / k t (different from k _T in Eq. 9), it is assumed that for PANS simulation, the grid spacing is located between the turbulent length scale (l _turb ) and the Kolmogorov length scale (η) (Song and Park, 2009). Therefore, the total turbulent kinetic energy k _t can be obtained as follows: k t = ∫ 0 ∞ E ( κ ) d κ = ∫ κ l κ Ω E ( κ ) d κ + ∫ κ Ω κ η E ( κ ) d κ = ∫ κ l κ Ω β ε 1 / 2 Ω 1 / 2 κ − 2 d κ + ∫ κ Ω κ η α ε 2 / 3 κ − 5 / 3 d κ = ( 3 2 α − β ) ε Ω − 1 + 1 2 π β ε 1 / 2 Ω 1 / 2 l t u r b − ( 1 2 π ) 2 3 3 2 α ε 2 / 3 η 2 3 . (11)

For the small turbulent length scale (large wavenumber, κ _Δ ∼κ _η ), the turbulent kinetic energy cannot be resolved directly; thus, if κ Δ < κ Ω , k u = ∫ κ Δ ∞ E ( κ ) d κ = ∫ κ Δ κ Ω E ( κ ) d κ + ∫ κ Ω κ η E ( κ ) d κ = ∫ κ Δ κ Ω β ε 1 / 2 Ω 1 / 2 κ − 2 d κ + ∫ κ l κ η α ε 2 / 3 κ − 5 / 3 d κ = ( 3 2 α - β ) ε Ω - 1 + 1 2 π β ε 1 / 2 Ω 1 / 2 Δ − ( 1 2 π ) 2 3 3 2 α ε 2 / 3 η 2 3 , (12)

and if κ Δ > κ Ω , k u = ∫ κ Δ ∞ E ( κ ) d κ = ∫ κ Δ κ η E ( κ ) d κ = ∫ κ Δ κ η α ε 2 / 3 κ − 5 / 3 d κ = ( 1 2 π ) 2 3 3 2 α ε 2 / 3 ( Δ 2 3 − η 2 3 ) . (13)

Based on the study proposed by Thangam et al. (1999), it was assumed that α = β = 1.8 . Consequently, the final f _k form for the rotation-corrected energy spectrum is as follows: f k = k u k t = { 0 . 9 ε Ω - 1 + 0.286 ( ε Ω ) 1 / 2 Δ − 0.793 ( ε η ) 2 / 3 0.9 ε Ω − 1 + 0.286 ( ε Ω ) 1 / 2 l t u r b − 0.793 ( ε η ) 2 / 3       ( κ c ≤ κ Ω ) 0.793 ε 2 / 3 ( Δ 2 / 3 − η 2 / 3 ) 0.9 ε Ω − 1 + 0.286 ( ε Ω ) 1 / 2 l t u r b − 0.793 ( ε η ) 2 / 3       ( κ c > κ Ω ) . (14)

It should be noted that when κ Δ = κ Ω , the above two equations are equivalent. Here, the f _k equation based on the rotation-corrected energy spectrum is referred to as PANS-RCES, and the corresponding coupled model with SST k-ω PANS is called SST k-ω PANS-RCES hereafter. In all PANS simulations, the dispassion of turbulent kinetic energy, ε , is calculated by ε = β k ω , where β is a constant coefficient.

In rotating channel flows, the channel rotates in the spanwise direction with a constant angular velocity and is significantly affected by rotating flow characteristics, such as typical asymmetric distribution, streaks near the pressure side, and Taylor–Götler (TG) vortices (Grundestam et al., 2008), which is suitable for verifying the new f _k expression based on the rotation-corrected energy spectrum. In the present study, rotating channel flows with a Reynolds number of 7,000 and rotation numbers of 0.3 and 0.6 are simulated. The Reynolds number is defined as Re = U m h / ν , where h is the half channel height and ν denotes the kinematic viscosity. In addition, rotation number Ro is defined as R o = 2 ω h / U m , where ω denotes the angular velocity of the rotation and U m = 1 / 2 h ∫ 0 2 h U ( y ) d y is the bulk velocity along the x coordinate. The simulation domain is exhibited in Figure 2, where 2 π h × 2 h × 2 π h are the dimensions along with the x (streamwise), y (normal), and z (spanwise) directions, respectively. The positive x coordinate is the inflow direction.

Subsequently, the SST k-ω PANS-ES1, SST k-ω PANS-ES2, and SST k-ω PANS-RCES models are used to simulate the rotating channel flow. The new turbulence models are compiled in OpenFOAM, and the simulation results are compared with the DNS data (Yang et al., 2012). For all PANS simulations, the PISO algorithm is applied to pressure–velocity coupling. A “Gauss linear” scheme with second-order accuracy is chosen for both the gradient term and the divergence term. A second-order implicit backward scheme is used for the time scheme.

The simulation results of the velocity and Reynolds stress analyzed in this study were time-, spanwise-, and streamwise-averaged: ϕ ¯ = 1 L x L z T ∫ 0 L x ∫ 0 L z ∫ 0 T ϕ d x d z d t , (15) where ϕ is a transient physical quantity, L _x and L _z are the streamwise and spanwise lengths, respectively, and T is the duration for time averaging. The half-height of the channel, h, is set as the reference scale for the length. Furthermore, the plane y/h = 0 is at the position of the pressure side, and y/h = 2 is at the position of the suction side.

The flow in a centrifugal pump is another typical rotating flow type. The large curvature and multiwall characteristics of the pump structure make the accurate prediction of the internal flow very difficult. In this study, a low specific speed centrifugal pump is used to verify the applicability of the new f _k expression based on the rotation-corrected energy spectrum. The pump has a rotating speed of 725 r/min, inlet diameter of D ₁ = 71 mm, and outlet radius of R ₂ = 95 mm, and the Reynolds number based on the inlet diameter D ₁ and the rated fiow rate Q ₀ was approximately 5.5×10⁴. The calculation domain is illustrated in Figure 10, and the grid points are approximately 2.78 million, which is determined based on previous simulations of this flow case (Byskov et al., 2003; Huang et al., 2015). In addition, the flow rate at the design condition was 3.06 L/s, and at the stall condition, it was 0.76 L/s. The corresponding experimental and simulation results are found in the studies by Pedersen et al. (2003) and Byskov et al. (2003).

Under the design conditions, the streamlines are smooth. Nevertheless, under part-load conditions, a stall happens in the impeller. More specifically, under stall conditions, a stall vortex is generated and blocks the channel. Therefore, the flow becomes more unstable, and the requirements for the turbulence model become higher (Pedersen et al., 2003; Zhou et al., 2014; Tao et al., 2014; Yao et al., 2016; Huang X.-b. et al., 2019). In this study, both the design-load condition, 1.0Q ₀, and the part-load condition, 0.25Q ₀, are used to verify the accuracy of the SST k-ω PANS-RCES (abbreviated as PANS hereafter) model. The mass flow rate is used as the inlet boundary condition, a pressure outlet boundary condition is used at the outlet, and a no-slip condition is set at the walls.

The Q = 1.0Q ₀ flow condition results of the velocity vector distributions in the impeller mid-height at the different radial positions of r/R ₂ = {0.65, 0.75, 0.90, 1.01} obtained by PANS and laser Doppler velocimetry (LDV) (Pedersen, 2000) are depicted in Figure 11. Under the design-load condition, the velocity vectors in the six flow channels are the same. The main feature is that the velocity at the suction surface (SS) of blades is large at 0.65 R ₂ and 0.75 R ₂, and the velocity at the pressure surface (PS) is small. The velocity distribution at 0.90 R ₂ is more uniform. At 1.01 R ₂, the velocity of the pressure surface is slightly greater than the velocity of the suction surface. This may be due to the secondary flow of the low specific speed centrifugal pump that causes the jet-wake flow at the exit to occur (Brun and Kurz, 2005; Zhang et al., 2019). The results of PANS and LDV show that PANS can capture the above phenomenon well. A more detailed comparison of the velocity distribution at 0.5 R ₂ and 0.75 R ₂ is shown in Figures 12, 13. From the figures, it is more obvious that the results of PANS are consistent with those of LES and PIV results, indicating that the velocity distribution can be well predicted in both near-wall and mainstream regions, and the newly developed f _k expression RCES can accurately predict the flow field in the design-load condition.

The velocity vector in the impeller mid-height under the part-load condition is shown in Figure 14. Under this working condition, the velocity vector of adjacent flow passages is different, which can be marked as non-stall Passage A and stall Passage B, respectively. It can be seen from Figure 14 that only one vortex appears on the suction surface in the non-stall Passage A. A smaller vortex appears at the inlet of the stall passage, and a larger one appears near the exit of the passage. The results of PANS are consistent with those of LES. Figure 15 shows the velocity vector results of r/R ₂ = {0.50, 0.65, 0.75, 0.90}. For stall passage, backflow occurs on the suction surface at the inlet, which is an important factor to generate stall vortices. The results of PANS and LES can accurately predict this result. For other flow direction positions, PANS also shows good prediction accuracy.