CFD methodology

A detached eddy simulation (DES) code was developed for investigating the internal flow. The Computational Fluid Dynamics (CFD) code was constructed to solve the three-dimensional compressible Navier-Stokes equations. The convective flux was evaluated using the flux difference splitting (FDS) method and extended to the third-order MUSCL interpolation. The viscous flux was evaluated as a second-order central difference. The Matrix Free Gauss-Seidel (MFGS) implicit algorithm was employed for time integration. The shear stress transport (SST) k-ω turbulence model was used in this code and applied to DES. This model depends on the local turbulent length scale and grid spacing size. The dissipation term in the k-equation of the SST k-ω turbulence model is modified as follows:

where l is the turbulent length scale, Δ is the local grid spacing, and the model’s constant coefficient C_DES was set to the recommended value by the report (Strelets, 2001). The several studies reported that DES analysis is more suitable than URANS for detailed flow field investigations near stall in both axial and centrifugal compressor (Im et al., 2012; Yamada et al., 2017). Our In-house code was successfully validated for unsteady flow fields of various turbomachinery flows (Fujisawa et al., 2019b; Zhang et al., 2020).

The computational domains applied in the numerical simulations are illustrated in Figure 2. The grid system included 14 impeller passages, vaneless diffuser, and a volute. The impeller (including clearance region) and vaneless diffuser systems had 32.0 and 11.0 million cells, respectively. The volute region had 22.0 million cells. In total, the computational grid had 65.0 million cells. The wall cell width was 0.1 µm, corresponding to a y+ parameter, equal to approximately one along all solid surfaces. At the inflow boundary, the total pressure and temperature were fixed, while mass flow rates were fixed at the outflow boundary. Across a sliding boundary separating the moving impeller and stationary diffuser frames, the most recent data on one side were interpolated to obtain opposite side data by using an unsteady simulation sliding mesh. Nonslip and adiabatic conditions were adopted for the wall conditions. The main CFD setup was summarized in Table 2.

Figure 2.

Numerical Mesh.

Compressor performance and stall characteristics

The experimental and numerical compressor performance results are shown in Figure 3. The numerical total pressure-rise characteristics were obtained from the time-averaged results of unsteady DES analysis and steady RANS analyses. The flow and total pressure-rise coefficients are defined as follows:

Figure 3.

Compressor Performance.

The unsteady simulations were conducted at the operational design point (ϕ = 0.24), and two off-design points (ϕ = 0.10 and 0.08). The steady simulations were conducted at three operational points from ϕ = 0.24 to ϕ = 0.14. The steady RANS analysis results were in good agreement with measured results. The unsteady CFD results were overestimated particularly at the off-design points as compared with the test results. The mixing loss especially within scroll passages was assumed to be underpredicted at near stall point, but the main factor of overprediction of pressure rise coefficient still was not confirmed. However, the relative error between experimental and unsteady CFD results at off design points was within 5%, as shown in Figure 3. Therefore, it was no problem to investigate the characteristics of a diffuser stall because the relative error was small.

To investigate the unsteady rotating stall characteristics, the radial velocity fluctuation of impeller discharge flow was measured with a split-film anemometer at ϕ = 0.08. The radial velocity was measured at D.I.0, 90, 180 and 240, which were indicated in Figure 1. The fast Fourier transformation (FFT) spatial distribution results at each measurement point in the spanwise direction are shown in Figure 4. The vertical axis represents the radial velocity fluctuation magnitude. From Figure 4, two types of velocity fluctuations were observed. One occurred at around 58 Hz and the other occurred at around 25–30 Hz. Based on previous report (Fujisawa et al., 2016), the unsteady phenomena at around 58 Hz were caused by a rotating stall within the impeller passages. On the other hand, the lower frequency fluctuation at around 25–30 Hz, was equivalent to 25–30% of the impeller rotating speed, and was caused by a rotating stall within the vaneless diffuser. Particularly, focusing on the diffuser stall, its characteristics varied depending on the circumferential position. From Figure 4(i), the small velocity fluctuation magnitudes caused by diffuser stall were observed near the shroud wall at D.I.0. Whereas, from Figure 4(ii), the large velocity fluctuation magnitude around 28 Hz was observed near the hub wall at D.I.90. Furthermore, Figure 4(iii) and (iv) show that the largest velocity fluctuation magnitude caused by the diffuser stall was observed near both shroud and hub walls at D.I.180 and 240, located near the cutoff. The mass flow fluctuations at the diffuser exit obtained from CFD analysis are shown in Figure 5. The red colored line indicated low mass flow regions and this rotational speed was about 25% of the impeller rotation speed. Therefore, the red colored low mass flow regions was assumed to be caused by the rotating diffuser stall, which rotated at 25% of the impeller rotation speed. This low mass flow region appeared at 90° and disappeared after passing the cutoff. Both experiment and CFD analysis indicated that the diffuser rotating stall initiated at 90° and developed further as the diffuser stall cell approached the cutoff. After passing the cutoff, the diffuser rotating stall was attenuated. The detailed unsteady diffuser stall behavior is discussed later.

Figure 4.

Distributions of Velocity Fluctuation in Spanwise Direction at Diffuser Inlet (EXP:ϕ = 0.08).

Figure 5.

Mass Flow Fluctuation at Diffuser Exit (CFD:ϕ = 0.08).

Generation mechanism of diffuser stall

To investigate the detailed unsteady diffuser stall behavior, an unsteady DES analysis was conducted. First, to understand the flow field, Figure 6 shows the distributions of radial velocity fluctuation within the diffuser passage obtained by CFD analysis. The radial velocity fluctuation was defined as follows:

Figure 6.

Distribution of Radial Velocity Fluctuations (CFD:ϕ = 0.08).

V_{r_ave} indicates the time-averaged radial velocity. The red color indicates the large magnitudes of radial velocity fluctuations. The region indicated by the white dashed circle is the diffuser stall cell that induced the largest mass flow fluctuations in Figure 5. Furthermore, the time-averaged mass flow rate of impeller discharge flow obtained from CFD analysis at ϕ = 0.08 is shown in Figure 7. The impeller exit plane mass flow rate was calculated on sections obtained by dividing the diffuser inlet into 32 circumferential sections. At t* = 2.2, the region with the large velocity fluctuations, as shown in Figure 6, was first generated at around 45°, where a strong circumferential adverse pressure gradient is observed, as depicted in Figure 8. Then, the large fluctuations region was formed at the diffuser exit area by the accumulation of several small backflow regions at t* = 2.9. After that, the diffuser stall cell expanded both radially and circumferentially as it approached the cutoff at t* = 3.7 and 4.2. Finally, the stall cell size reduced after passing the cutoff at t* = 5.4 due to increased impeller discharge mass flow rate at near the cutoff, as shown in Figure 7.

Figure 7.

Distribution of Mass Flow at Impeller Exit (CFD:ϕ = 0.08).

Figure 8.

Distribution of Static Pressure in Circumferential Direction (CFD:ϕ = 0.24 and 0.08).

The circumferential static pressure distributions at ϕ = 0.24 and ϕ = 0.08 obtained from the CFD analysis are shown in Figure 8. The static pressure was measured at shroud wall D.M. From this figure, at around 300°, the static pressure at off-design point ϕ = 0.08 was the same level as that at design point ϕ = 0.24. However, from 0° to the cutoff, the static pressure at ϕ = 0.08 was higher than that of ϕ = 0.24. The magnitude of circumferential adverse pressure gradient was strong near 0° for ϕ = 0.08. Therefore, the diffuser stall was initiated at around 45°. Furthermore, the time-averaged meridional distribution of radial velocity and streamlines at 0°, 45°, and 240° are shown in Figure 9. The boundary layer separation occurred at hub and shroud wall by turns at 240deg, with numerous papers reporting the same phenomena (Van den Braembussche, 2019). From this paper, the spanwise pressure gradient induced the flow concentration on the side with the highest pressure. In addition, the spanwise pressure gradient inversion was generated once separation occurred, and this pressure gradient change continued until flow reversal took place on the other side. Therefore, the separations generated at both hub and shroud by turns. First, the flow concentrated on the shroud side where the total pressure was highest and boundary layer separation occurred on the hub side at 240deg (point (a)), as shown in Figure 9. At point (b), the pressure gradient inversion occurred and the flow concentrated on the hub side. Therefore, the next boundary layer separation occurred at the opposite wall. Particularly, at 45deg, the separation region was observed near the diffuser inlet shroud wall. Accordingly, the hub side boundary layer separations at around 45°, triggering the diffuser stall cell generation, were caused by the shroud side separation region at the diffuser inlet and a large circumferential adverse pressure gradient. Figure 10 shows the time-averaged limiting streamlines on hub wall obtained by measurement and CFD analysis. In the experiment, the limiting streamlines were visualized by using an oil film method. The limiting hub wall streamlines were bent at around 45° where a small backflow region generated.

Figure 9.

Meridional Distribution of Radial Velocity and Streamlines (CFD:ϕ = 0.08).

Figure 10.

Limiting streamlines on hub wall (EXP and CFD:ϕ = 0.08).

For improved understanding of the diffuser stall cell generation mechanism, the unsteady flow field was visualized. The instantaneous flow field at two different times are shown in Figure 11. The distributions of radial velocity fluctuation within the diffuser passage at each time are shown. The developing behavior of diffuser stall cell was investigated from Figures 12 and 13. Figures 12(i) and 13(i) show the meridional radial velocity distribution around the region A and B at t* = 2.13 and 3.00. Figures 12(ii) and 13(ii) show the iso-surface of V_r = −10 m/s colored by the distance from the hub wall. From Figure 12, the hub wall boundary layer separation initiated at around 45°. Then, this hub wall separation induced the next shroud wall separation because of a spanwise pressure gradient, as shown in Figure 13. The boundary layer separations on the shroud wall developed, as it can be seen by the low velocity region at diffuser exit. The hub side Backflow Region I was the impeller wake. From numerical results, the impeller wake mainly expanded on the hub side from 90° to the cutoff. The shroud side Backflow Region II consisted of the boundary layer separations. The region indicated by the black circle was the low velocity region. Backflow Region I rotated faster than Backflow Region II, as shown in Figure 13a and b. Therefore, Region I caught up with Region II and induced the boundary layer separation enlargement on the shroud wall (Region II). Furthermore, the shroud wall boundary layer separations (Region II) induced the next hub wall boundary layer separations (Region III) because of a spanwise pressure gradient. The low velocity region within vaneless diffuser developed due to the boundary layer separations occurring on the shroud and hub wall by turns. Finally, the low velocity region formed the entire diffuser passage span blockage.

Figure 11.

Instantaneous Radial Velocity Fluctuation at t* = 2.13 and 3.00 (CFD:ϕ = 0.08).

Figure 12.

Generation Mechanism of Diffuser Stall Cell at t* = 2.13 (CFD:ϕ = 0.08).

Figure 13.

Development of Diffuser Stall Cell at t* = 2.88 and 3.00 (CFD:ϕ = 0.08).

In order to validate the numerical results, the radial velocity fluctuation of impeller discharge flow was measured at the meridional direction (D.I.150, D.M.150 and D.E.150), as shown in Figure 14. The diffuser stall fluctuations were observed near hub side at D.I.150. On the other hand, the region with large magnitudes of the diffuser stall was shifted to the shroud side at D.M.150 and D.E.150. It is because the diffuser stall was consisted of the boundary layer separations occurring on the shroud and hub wall by turns. Therefore, the boundary layer separations occurring on the shroud and hub wall were considered as the key factors contributing to the diffuser stall generation by both experimental and CFD results. In our future work, the general generation mechanism of the boundary layer separations, which occurred at the hub and then the shroud end wall, will be revealed by several CFD investigations.

Figure 14.

Distributions of Velocity Fluctuation in Meridional Direction (EXP:ϕ = 0.08).