Introduction
Reynolds-averaged Navier-Stokes (RANS) simulation is the industrial workhorse for predicting turbomachinery flows, but users can occasionally find inconsistencies in RANS solvers such as:
Using the “same” numerical settings from different RANS solvers can predict different results.
Each RANS solver has a “best-practice” list of numerical settings based on previous user experience, but this list may vary from one solver to another.
When using the “best-practice” settings of a solver, the quantity of interest may be over-predicted in one case and under-predicted in another, despite the flow phenomena of the cases being similar.
The source of RANS solver inconsistencies can be rooted back to the uncertainties and errors in experiments and simulations. In terms of experiments, uncertainties and errors in geometries and test instruments vary for different test cases and test rigs; in some cases, the experimental uncertainties and errors may be pronounced enough to overwhelm the difference in RANS results. In terms of simulations, the prediction is subjected to many user-specified numerical settings; comparisons among different RANS solvers at the same numerical settings are seldom reported due to user preference (i.e., following the “best-practice” guideline of a specific solver), coding preference (e.g., an in-house variant of a model is implemented instead of the original one) and potential coding errors.
To tackle the issue of RANS solver inconsistency, a series of validation and verifications (V&V) studies of RANS solvers are crucial. Following the definitions from the American Institute of Aeronautics and Astronautics (AIAA) (Computational Fluid Dynamics Committee, 1998):
Validation is the process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended uses of the model.
Verification is the process of determining that a model implementation accurately represents the developer's conceptual description of the model and the solution to the model.
V&V studies have drawn attention from the general aerospace sector for the past ten years (Rumsey et al., 2010). In response, the NASA Turbulence Modeling Resource (TMR) website (Rumsey, 2021) was established, which serves as an open repository of RANS turbulence model documentation and V&V test cases. However, V&V studies have not yet received as much attention in the turbomachinery community as in the external flow community. The last major V&V campaign for compressor flows dates back to the 1994 International Gas Turbine Institute (IGTI) blind test on the NASA Rotor 37 (Denton, 1997), along with many individual validation efforts on other NASA compressor geometries (Suder, 2021) over the years. From today's perspective, the measurement accuracy of these data may no longer present state-of-the-art. Online and version-controlled documentation of these geometries and measured data was also not achieved due to the limitations of the times—users today need to find the most updated geometries and data from pieces of literature by themselves (e.g., characteristics and tip gap size of NASA Rotor 67 (Grosvenor, 2007)), and some details of the geometry are still not known (e.g., running tip gap size of rotors, hub cavity geometry of stators). In addition, it is challenging to conclude from the early V&V campaign due to a lack of solver verification (Denton, 1997). Therefore, a new V&V campaign of turbomachinery RANS solvers based on recent common research models of compressors is of high interest.
In the 2020 GPPS Chania conference, the Institute of Gas Turbines and Aerospace Propulsion (GLR) of Technische Universität Darmstadt released the TUDa-GLR-OpenStage transonic axial compressor including the geometry in digitized form and an abundant set of measurement data (Klausmann et al., 2021). This is an ideal compressor test case for RANS V&V purposes. In this work, established commercial turbomachinery RANS solvers are performed on this case to provide the benchmark solutions. These benchmark solutions along with the associated grids and boundary conditions will be released to the public for future verification of other RANS solvers. Based on the commercial solvers, the effects of grid density, advection scheme, turbulence model and rotor-stator interface model on predicting the TUDa-GLR-OpenStage case are investigated. These results will be presented after introducing the experimental and numerical setup of the case.
Methodology
Case description
The investigated TUDa-GLR-OpenStage is a single-stage high-speed axial compressor, which represents a typical front stage of a high-pressure compressor in a commercial turbofan engine. The compressor stage consists of a blisk rotor with radially stacked CDA-airfoils, a stator with optimized 3D shapes, and an outlet guide vane (OGV) that straightens the flow. The rotor is originally designed by MTU and first tested in 1994. It has been investigated extensively in a series of previous experimental and numerical works (Hoeger et al., 1999; Bergner et al., 2006; Muller et al., 2011). The rotor blades are highly cambered near the hub and thin near the tip. The running tip gap size is approximately 0.75 mm at the design condition, and the hub fillet radius is 5 mm. The stator is designed conjointly between the GLR and the German Aerospace Center (DLR). It is optimized via an automated multi-objective optimization process to suppress the separation size (Bakhtiari et al., 2015). The stator blades have a forward sweep feature near both endwalls and a bow feature towards the pressure surface near the casing. The shapes of the fillets at both endwalls are prescribed by a digitized geometry file. More details about the investigated stage are summarized in Table 1.
Table 1.
Test rig
The TUDa-GLR-OpenStage was tested on the transonic compressor rig of GLR, which can measure the steady-state performance, the aerodynamic instabilities and the blade vibration levels of the compressors. The schematic of the test facility is shown in Figure 1. The inlet flow passes through an inlet throttle, a settling chamber and a mass flow measurement section before reaching the compressor core. The compressor is driven by a DC motor with a gearbox. Shaft input torque and rotor speed is measured using a torque meter.
The measurement sections of the compressor are also marked in Figure 1. At the stage inlet section (ME15), radial profiles of total pressure and total temperature are measured. At the rotor exit section (ME21), radial profiles of total pressure and absolute Mach number components are available. At the stage exit section (ME30), 2D distributions of total pressure and total temperature are acquired. The total pressure ratio and total temperature ratio characteristics are calculated by the area-averaged probe data at sections ME15 and ME30; the isentropic efficiency is calculated using the probe-based total pressure ratio and the shaft power. For a fair comparison between the experiment and CFD, the isentropic efficiency from CFD is calculated based on the area-averaged pressure ratio and the mass-averaged temperature ratio. Calculations of the measurement uncertainty along with the measured data are detailed in the concurrent experimental work (Klausmann et al., 2021).
The measurements were performed at 100%, 80% and 65% design speed. This work mainly focuses on the 100% design speed at the peak-efficiency (PE) condition
CFD solvers
In this work, two widely used commercial RANS flow solvers are adopted: Ansys CFX 20.1 and Numeca FineTurbo 14.1.
Ansys CFX uses an element-based finite volume method to discretize the RANS equations based on unstructured meshes. The diffusion term is discretized by the second-order accurate central scheme, and the advection term is discretized by the second-order accurate high resolution scheme. The temporal discretization is achieved by the second-order accurate backward Euler scheme. The linearized equations are solved by the Incomplete Lower Upper (ILU) factorization technique, which is accelerated by an algebraic multigrid method. A pseudo time-stepping algorithm with an automatic time scale technique is used for steady-state calculations. For best-practice of turbomachinery simulation, the updated Menter
Numeca FineTurbo discretizes the RANS equations using a cell-centered control volume approach based on structured meshes. The diffusion term is discretized by the second-order accurate central scheme. The advection term is discretized by either the central scheme with a Jameson-Schmidt-Turkel (JST) type dissipation term (Jameson et al., 1981) or the upwind scheme using flux difference splitting from Roe (1981), both of which are second-order accurate. The temporal discretization is achieved by the explicit fourth-order accurate Runge-Kutta integration scheme. For steady-state cases, the pseudo time-stepping algorithm is used. Local time stepping, implicit residual smoothing, and multigrid techniques are applied to accelerate convergence. The modified SA turbulence model (Ashford, 1996) (i.e., “SA-fv3-noft2”1 [1] (Rumsey, 2021)) is set as the default for turbomachinery simulations. Other turbulence models including the modified SARC model (Shur et al., 2000) (i.e., “SA-RC-fv3-noft2” (Rumsey, 2021)), the original SST model (Menter, 1994) (i.e., “SST” (Rumsey, 2021)) and the simplified baseline EARSM model (Menter et al., 2012) are also available. More details of the solver can be found in its user guide (NUMECA International, 2021).
Flow domain, boundary conditions and grids
The simulation is conducted by using a single-passage flow domain with the circumferentially periodic boundary condition, as illustrated in Figure 2a. The inlet plane is set at the stage inlet section ME15. The measured total pressure and total temperature profiles are used to prescribe the inlet boundary condition. In addition, the inlet flow direction is axial; the turbulence intensity and turbulence length scale are estimated as 4% and 0.09 mm respectively, which are equivalent to an eddy viscosity ratio of 35. At the rotor-stator and the stator-OGV interface, a mixing plane boundary condition (specifically, “Stage (Mixing-Plane)” with “Constant Total Pressure” as downstream velocity constraint for Ansys CFX; “Full Non-Matching Mixing Plane” for Numeca FineTurbo) is used. The outlet boundary is set at the location about 1.5 times the compressor core axial length downstream of the OGV. A radial equivalent backpressure is assigned at the outlet boundary, but a constant backpressure value can also be used since the outlet flow is turned to the axial direction by the OGV. When approaching the stall limit, a step increase of backpressure of 200 Pa is used. The last converged point is taken as the numerical “stall” point.
The flow domain is discretized by hexahedron grids with an O4H topology, which are generated by using the Numeca AutoGrid5. In total, five sets of rotor and stator grids with different grid densities are generated. The mesh refinement is achieved by uniformly increasing the grid number to about 1.5 times in all dimensions, while the first layer grid displacement is kept the same at an area-averaged
Table 2.
Results
Effect of grid density
In this section, five sets of uniformly refined grids (detailed in Table 2) are adopted to investigate the effect of grid density. Unless stated elsewhere, the rest of the numerical settings are kept the same, including the Ansys CFX solver, the SST turbulence model and the mixing plane rotor-stator interface model.
The stage performance characteristics at the design speed are presented in Figure 3, where the simulation results are represented by curves with open symbols, and the experimental data are denoted by solid squares with error bars. As the grid density increases, the characteristic curves gradually translate towards the top-right corner (i.e., higher pressure ratio, efficiency and mass flow rate; similar slopes of pressure ratio and efficiency curve). Excluding the UltraCoarse grid result, the stable flow range
To examine the effect of grid density on the flow fields, the circumferentially mass-averaged radial profiles downstream of the stator are presented in Figure 4. In these plots, curves represent the numerical results obtained from the five grids, and solid squares denote the experimental results. For simplicity, only the PE results are presented. Results show that grid refinement leads to an increased total pressure ratio and total temperature ratio, especially at the upper 20% span. The Fine and UltraFine results are almost overlapping with each other, indicating grid independence has been achieved. However, they over-predict the total pressure ratio in the upper and lower 20% spans near the endwalls. Such a deficiency can be caused by the deficiency of the turbulence model, the deficiency of the mixing plane model, or plausible geometric uncertainties and errors. These factors will be examined later.
Following the concept of Roache (1998), the grid discretization error scales with
Comparing the investigated five grids, the choice of grid density yields a 7% relative (to experiment) variation in total pressure ratio and a 7% absolute variation in isentropic efficiency. The discretization errors of the Fine grid, i.e.,
Comparing the investigated four turbulence models implemented in the two solvers, the choice of turbulence model yields a 5% relative variation in total pressure ratio and a 5% absolute variation in isentropic efficiency for the discretization-error-free values
Effect of Advection Scheme
In this section, the JST central scheme (Jameson et al., 1981) and the Roe upwind scheme (Roe, 1981) are compared with each other, both of which are second-order accurate. The rest of the numerical settings are kept the same, including the Numeca solver, the Fine grid, the SST turbulence model and the mixing plane rotor-stator interface model.
The discretization error related to a second-order accurate advection scheme scales with the square of the grid spacing. Thus, the difference between the schemes is expected to be evident only when a coarse grid is used. In Figure 6a, the grid convergence of different advection schemes on predicting the compressor total pressure ratio at the PE condition is presented. It is confirmed that the choice of advection scheme has a negligible effect when refining the grid to the Fine grid level.
For a more detailed comparison, the grid convergence of different advection schemes on predicting the rotor suction surface shock front location is presented in Figure 6b. Here, the shock front is defined at the grid point where
For best practice, it is the choice of grid density rather than the choice of advection scheme that should be paid more attention to. Although any advection scheme can achieve an accurate prediction once the grid independence requirement has been met, a high-order scheme is still favorable in terms of computational cost.
Effect of turbulence model
In this section, four turbulence models namely SST, SA, SARC and EARSM are compared with each other. The rest of the numerical settings are kept the same, including the Numeca solver with the central advection scheme, the Fine grid and the mixing plane rotor-stator interface model.
The stage performance characteristics at the design speed are compared with experimental data in Figure 7, where the measured data are represented by solid squares with error bars, and the simulation data are denoted by lines with scatters. In general, the SA model and the EARSM model predict the best agreement with the measured data. The SA model reproduces the stall margin2 [2] (SFR = 13.0% for the experiment; SFR = 14.3% for the SA prediction), the NS pressure ratio and the PE efficiency; however, the PE pressure ratio is under-predicted and the NS efficiency is over-predicted. The EARSM model accurately predicts the PE pressure ratio, but it under-predicts the stall margin (SFR = 9.6%) and slightly over-predicts the NS pressure ratio and the overall efficiency. The SST model predicts the highest pressure ratio and efficiency, but it is close to the EARSM results. The SARC model predicts the lowest pressure ratio and efficiency, which shows the least agreement with the measured data.
For further comparisons between the four turbulence models, the radial profiles at the rotor exit and the stage exit are examined in Figure 8. In these plots, simulation results are represented by lines, and the measured data are denoted by solid squares.
At the PE condition, both the SA model and the SARC model under-predict the pressure ratio and the flow angle at the rotor exit below 80% span, as shown in Figure 8a. Such deficiencies generally propagate to the stage exit shown in Figure 8b, which accounts for the under-prediction of the overall stage pressure ratio at the PE condition. The EARSM model and the SST model slightly over-predict the rotor exit pressure ratio above the 80% span and the rotor exit flow angle at full span, as presented in Figure 8a. At the stage exit shown in Figure 8b, the near-tip over-prediction of pressure ratio in the rotor passage is propagated to the stage exit; in addition, the near-hub over-prediction of pressure ratio arises, which is induced in the stator passage. These deficiencies of the endwall pressure ratio prediction, however, are too small to be observed in the overall stage pressure ratio shown in Figure 7a.
At the NS condition, the SA model predicts the most accurate pressure ratio profile at the stage exit, as shown in Figure 8d; the SARC model still under-predicts the pressure ratio at most spans, whereas the SST model and the EARSM model over-predict it especially near the endwalls. The shapes of all the pressure ratio profiles at the stage exit are similar to that at the rotor exit shown in Figure 8c, except that the SST result shows a total pressure loss near 80% span at the stage exit plane. This drop of total pressure corresponds to an excessive prediction of the corner separation size in the stator, as will be shown later. In terms of the stage exit total temperature ratio shown in Figure 8d, the EARSM model and the SST model show the best agreement with the measured data, whereas the SA model and the SARC model slightly under-predict it.
To complement the validation on the radial profiles, the total pressure ratio contours at the stage exit are presented in Figure 9. In these plots, the top row and the bottom row represent the PE condition and the NS condition, respectively; the columns from left to right represent the experiment, SST, SA, SARC and EARSM results.
For the experimental results, a hub-suction surface corner separation is observed at both the PE and NS conditions, which is featured by low values of total pressure. As the compressor operates from PE to NS, the separation size increases. This corner separation, however, is not reproduced by any of the simulation results, leading to the near-hub over-prediction of the total pressure ratio observed previously in Figure 8b and d.
Among the simulation results, the SARC model and the SST model predict the lowest and highest overall total pressure ratio, respectively, which deviates from the measured data the most. The SST model also over-predicts the near-tip trailing edge separation size at the NS condition. Despite the EARSM model and the SA model predicting the best agreement with the measured data, evident differences still persist. Further investigations regarding the effect of the rotor-stator interface model and the effect of geometric uncertainties and errors (e.g., stator hub cavity leakage flow) are needed to improve the prediction accuracy.
To compare the four turbulence models in predicting the rotor flow field, the static pressure contours superimposed with the surface limiting streamlines on the rotor suction surface are presented in Figure 10. In general, the rotor flow field is featured by a near-hub corner separation and a near-tip shock-induced separation. When the compressor is throttled towards stall, the shock front moves further upstream, and the sizes of both separation cells increase. Comparing the effect of turbulence models, the SST model predicts the largest separation size, followed by the EARSM, SA and SARC model. In particular, the near-tip shock-induced separation is not captured by the SA model at the PE condition, as well as the SARC model at both the PE and NS conditions. Such results imply that the SA and SARC models predict higher eddy viscosity in the rotor passage than the other models, which will be examined latter.
To compare the four turbulence models in predicting the stator flow field, the static pressure contours superimposed with the surface limiting streamlines on the stator suction surface are presented in Figure 11. In general, the stator flow field is featured by a near-hub corner separation and a near-tip trailing edge separation. As the compressor operates toward stall, the sizes of both separation cells increase. Similar to the observations in the rotor flow fields, the SST model predicts the largest separation size, followed by the EARSM, SA and SARC model. In particular, at the NS condition, the SST model predicts a near-tip corner separation that connects with the trailing edge separation. This excessive prediction of separation size is consistent with previous observations in Figures 8d and 9.
To help understand the loss mechanism predicted by different turbulence models, the loss of efficiency
The radial profiles of
The investigated four turbulence models are all based on the concept of eddy viscosity. To help understand the behavior of the turbulence models further, the turbulent-to-laminar viscosity ratio contours are presented in Figure 13, where each row corresponds to an operating condition, and each column corresponds to a turbulence model. For the compactness of the figure, the inlet duct and the OGV are not shown.
For the two-equation turbulence models SST and EARSM, they are insensitive to the freestream eddy viscosity boundary condition (i.e.,
For one-equation turbulence models SA and SARC, they are more sensitive to the freestream eddy viscosity as
Effect of rotor-stator interface model
In this section, the mixing plane rotor-stator interface model typically used in steady simulations is compared with the sliding plane model used in unsteady simulations. The investigation is performed by using the Ansys CFX solver, the Fine grid topology and the SST turbulence model. To reduce the required computational resources, the number of stator blades has been changed from 29 to 32 so that a periodic sector of one rotor blade and two stator blades (1R2S) can be used for the sliding plane model. The time step used in the unsteady simulation corresponds to
The stage exit (ME30) radial profiles obtained from the mixing plane method and the time-averaged sliding plane method are compared in Figure 14a. In the first place, the steady-state performances before and after changing the number of stator blades are compared with each other. It shows that both cases have almost the same profiles, indicating the scaling of the stator geometry barely changes the flow physics. Then, the sliding plane result is compared with the mixing plane result. No evident difference is found except that the sliding plane result shows slightly better agreement with the measured total temperature ratio. Such results indicate the mixing plane model is reasonably accurate at the investigated PE and NS conditions. Finally, comparing all CFD results with the measured data, the over-prediction of the stage pressure ratio in the upper and lower 20% span still persists. Such deficiencies are potentially caused by the geometric errors (as will be discussed in the later section) and the non-ability of RANS to sufficiently capture the turbulence effect of endwall flows.
To unveil the effect of the rotor-stator interface model on predicting the eddy viscosity field, the turbulent-to-laminar viscosity ratio contours are presented in Figure 14b. Results show that the mixing plane results have an abrupt increase of
Discussions on geometric uncertainty and error
The above analyses are based on the assumption that the geometry used in the simulation is exactly the same as the one used in the experiment. In this section, the effects of some real geometry features are briefly discussed, including the rotor running tip gap shape, the rotor casing pinch and the stator hub cavity.
Effect of rotor running tip gap shape
In the experiment, the rotor running tip gap shape was measured by a capacitive blade tip timing system, and the tip gap size was found to increase by 0.2%
Effect of rotor casing pinch
In the experiment, a small pinch exists on the rotor casing endwall, but it is simplified as a smooth straight line in previous simulations. The shape of the pinch is illustrated in Figure 15a. In this section, the compressor with the realistic pinched casing is simulated by using the CFX solver, the Fine grid, the SST turbulence model and the mixing plane rotor-stator interface model. The radial profiles of the total pressure ratio are presented in Figure 15b. In general, the rotor casing pinch reduces the total pressure ratio at all spans, which improves the agreement with the measured data. This is because the pinched casing reduces the rotor throat area, increases the axial velocity (when comparing at the same mass flow), reduces the work input (by changing the rotor exit velocity triangle) and eventually reduces the total pressure ratio. Such a reduction effect is more evident in the PE condition than in the NS condition, because the slope of the total pressure ratio characteristic is steeper near the PE condition.
Effect of stator hub cavity
In the test rig, a secondary flow passage connects the stator inlet hub to the stator outlet hub, and a leakage flow may exist. To explore its effect on the prediction, a cavity geometry shown in Figure 16a is investigated numerically in this section. The simulation is conducted using the CFX solver, the Fine grid, the SST turbulence model and the mixing plane rotor-stator interface model. The total pressure ratio contours at the stage exit (ME30) are presented in Figure 16b. It shows that the stator hub cavity increases the near-hub total pressure loss at both the PE and the NS conditions, leading to slightly better agreement with the measured data. However, because the leakage mass flow rate of the investigated cavity geometry is only about 0.2% of the main flow rate, its effect is not strong enough to eliminate the near-hub disagreement with the measured data. According to previous research (Seshadri et al., 2015), a leakage mass flow rate of about 0.5% can potentially help match the measured data better. A parametric study of the cavity leakage flow rate will be pursued for future research.
Conclusions
A comprehensive validation and verification study of turbomachinery RANS flow solvers is conducted for the transonic axial compressor TUDa-GLR-OpenStage. Two commercial solvers namely Ansys CFX and Numeca FineTurbo are adopted to provide the benchmark solutions. Based on these solvers, five sets of grids, two advection schemes (i.e., central difference and second-order upwind), four turbulence models (i.e., SA, SA-RC, SST and EARSM) and two rotor-stator interface models (i.e., mixing plane and sliding plane) are investigated to quantify their effects on predicting the performance and the flow field of the compressor stage. Several conclusions are drawn as follows.
Among the investigated numerical setups, the choice of grid density and turbulence model has the most significant effects on the prediction. Specifically, the choice of grid density between
Regarding the choice of grid density, the method from Roache (1998) has been demonstrated useful in estimating the discretization error of the current compressor, and it is recommended for future research. The current grid independence study shows that using about 1 million and 3 million grid points per blade passage can achieve a small discretization error in the global performance characteristics and the local flow fields, respectively. These numbers can be used as a rule of thumb for future RANS simulations. It is also demonstrated that the grid independence obtained from one turbulence model is transferable to another.
Regarding the choice of turbulence model, the EARSM model achieves overall the best accuracy and thus is recommended among the four models. Despite the SST model showing similar global performance characteristics to the EARSM model, it predicts the largest separation size and over-predicts the stator tip corner separation when operating near stall. The SA model over-predicts the eddy viscosity in the inlet duct, leading to the excessive loss of efficiency in the blade passages and thus the under-predicted total pressure ratio at the peak efficiency condition. The SARC model is the least accurate model due to the excessive prediction of eddy viscosity in the inlet duct, which leads to the smallest separation size, the largest friction loss and the lowest total pressure ratio and isentropic efficiency.
When using the grid density and the advection scheme with a small discretization error, none of the investigated turbulence models and the rotor-stator interface models can perfectly capture the measured compressor flow fields. The deficiency mainly comes from the near-tip flow of the rotor and the near-hub flow of the stator, which are likely related to the turbulence modeling deficiency in wakes/tip leakage flows/corner separations and the unknown geometric error near the stator hub. Future V&V works should focus on the verification and evaluation of RANS turbulence models with advanced features (e.g., nonlinear constitutive relation, transition, roughness, etc.), scale-resolving simulations and the real-geometry effects (e.g., stator hub leakage flow).