Assessment of advanced RANS turbulence models for prediction of complex flows in compressors

Wei SUN

China Academy of Aerospace Science and Innovation, Beijing 100176, China

KEYWORDS

Abstract Reynolds-Averaged Navier-Stokes (RANS) Computational Fluid Dynamics (CFD) has been widely used in compressor design and analysis.However,reasonable prediction of compressor flow and its impact on compressor performance remains challenging.In this study,Menter’s Shear Stress Transport(SST)model and its variants,as well as the ω-based Reynolds stress model(Stress-BSL) are assessed.For a single rotor(Rotor 67),under the peak efficiency operating condition, all studied turbulence models predict its performance with reasonable accuracy; under the off-design conditions, SST with Helicity correction (SST-Helicity) shows superiority in predicting the effect of flow on the spanwise distribution of aerodynamic parameters.For Darmstadt’s 1.5-stage transonic axial compressor, SST-Helicity outperforms SST, SST with the Quadratic Constitutive Relation(SST-QCR)and Stress-BSL in predicting the performance as well as the spanwise distribution of aerodynamic parameters.At the design rotating speed, the stall margin given by SST-Helicity(20.90%) is the closest to the experimental measurement (24.81%), which is more than twice that by SST (8.71%) and 1.5 times that by SST-QCR (14.14%).This paper demonstrates that SSTHelicity model,together with a good quality and sufficiently refined grid,can capture the compressor flow features with reasonable accuracy,which results in a credible prediction of compressor performance and stage matching.

1.Introduction

In the past few decades, with the continuous progress of CFD as a predictive tool for complex flows in fluid machinery,CFD has been widely used in almost all fields of fluid machinery scientific research, technology development, engineering design,manufacturing, production, operation and maintenance.The progress of CFD as a design and analysis tool for fluid machinery is, to a large extent, attributable to the continued efforts for refining the RANS turbulence models,and the continued advent of faster and more capable computer as well as solver algorithms.The widespread application of CFD has fundamentally changed the related industries.One of the manifestations is that RANS/Unsteady RANS(URANS)-based CFD is now routinely used by turbomachinery designers to aid the compressor design.

Generally, under the design operating condition, RANS usually predicts the compressor performance and related flow field with reasonable accuracy.However,as the operating condition departs from the close-to-optimum design condition,the flow becomes increasingly complex, which results in a nonlinear growth in the uncertainty of RANS simulations.This is because under such operating conditions,the turbulent flows in compressors are characterized by much more complex physics,1such as stronger tip clearance flows, stronger shock-induced separations, stronger endwall secondary flows and large-scale Three-Dimensional (3D) separations.However, the commonly used RANS turbulence models are all based on oversimplified assumptions relating to the turbulence energy equilibrium, and the isotropic distribution of turbulence intensity.This indicates that for RANS CFD,reasonable prediction of compressor flows is still challenging, due to the gap between the complex nature of these flows and the relatively simple nature of the assumptions made to most of RANS models.

An exemplary flow that highlights the turbulence modelling deficiencies is the 3D corner separation formed at the junction of the blade suction surface and endwall of an axial compressor.Corner separation is a common flow phenomenon in compressors that can significantly affect compressor performance.When a compressor blade row operates away from its intended design condition, corner separation appears to be large and localized.This causes increased blockage in blade passages,limiting the pressure rise and efficiency achievable in the compressor.2–4Further, in the multi-blade row environment, the predictive discrepancy of the blockage due to corner separation in one blade row causes discrepancies of the endwall flow deviation at the exit of this row.This in turn changes the inlet flow conditions of the next blade row and can ultimately lead to the stage-loading mismatches across the whole compressor.The mismatches of stage loading in multi-stage compressors,as commonly observed in RANS simulations, are in many cases due to CFD overprediction of the corner separations throughout the whole compressor.5Therefore, it is necessary for turbulence models to be able to predict corner separation and its consequences with good accuracy, in order to make CFD a reliable design and analysis tool.

Along with the corner separation flow,the tip leakage flow plays a major role in setting the axial compressor stall onset,and the efficiency potential.6The passage blockage, caused by the co-existing corner separation and tip leakage vortex(generated by the free shear layer roll-up in the tip leakage jet), impairs the pressure-rise capacity, and if the blockage is large enough, it causes the pressure rise characteristic to peak and turn over.7Especially as the flow coefficient is reduced,for example, under higher off-design conditions, the normal component of the tip leakage jet velocity increases relative to its component along the blade, and thus the leakage jet turns upwind against the freestream direction.8In these scenarios,the blockage due to the tip leakage vortex becomes larger and may dominate in the blade passage.However, as far as can be evaluated from the available numerical studies on tip leakage flows, reasonable prediction of the tip-leakage vortex core size and growth rate, center position, vorticity strength and turbulent flow field near the blade tip remains a major challenge for RANS turbulence models.9–11

In view of the impact of corner separation flow and tip leakage flow on compressor performance,improving the modelling accuracy of the two flows and their interaction is of crucial importance, not only for compressor design, but also for further understanding of the flow mechanisms in compressors.To help understand what is missing from RANS turbulence models, a series of scale resolved numerical simulations, e.g.,Large Eddy Simulation (LES) and Hybrid RANS-LES, have been conducted on the compressor stator cascade,12–15rotor cascade,16–19compressor stage,20,21and multi-stage compressor22–25for investigating the compressor flow physics.The scale-resolving simulation results show that high turbulence anisotropy exists both in the tip leakage vortex and in the corner separation.Most RANS turbulence models cannot reasonably capture turbulence anisotropy because the Boussinesq hypothesis (linear Reynolds stress–strain constitutive relation)on which the model is based does not account for turbulence anisotropy.Furthermore, the Turbulent Kinetic Energy(TKE) budgets show that high turbulence non-equilibrium(unbalance of TKE production and dissipation) exists both in the hub/casing endwall boundary layer and in the corner separation.Again,this is out of reach of most of RANS turbulence models as they were formulated based on the turbulent characteristics of equilibrium turbulent boundary layer.

To improve RANS modelling accuracy for compressor flows,many studies have focused on improving the turbulence equation source terms and the Reynolds stress–strain relation.Liu et al.26demonstrated that in 3D vortical flows,turbulence energy backscatter (i.e., energy transferring from the smaller scales of turbulence to large scales) is significant, and there exists positive correlation between velocity helicity and turbulence energy backscattering.They first introduced helicity correction to improve turbulence models to capture the turbulence energy backscattering physics.26The improved Spalart-Allmaras (S-A) model and SST model (hereafter referred to as SA-Helicity and SST-Helicity)have been demonstrated to improve the prediction of corner separation.27Given the effectiveness of helicity correction,this correction has been incorporated as an expert model option into the AESim simulation software developed by Aero Engine Corporation of China (AECC).

To address the deficiency of Boussinesq hypothesis,Spalart28and Menter et al.29separately proposed the Reynolds-stress anisotropy model (QCR) and coupled it to SA model(SA-QCR)and SST model (SST-QCR).The principle of QCR is that the added quadratic strain-vorticity term induces counter-rotating vortex pairs at the corner.The counter-rotating vortices extract the high-momentum fluid from the main flow into the corner, energizing corner flow to be resistant to separation.Li and Liu3,4validated SA-QCR model against the Cambridge Prescribed Velocity Distribution(PVD) compressor cascade.They found that with the QCR corrections for considering turbulence anisotropy,the blending SA models perform better than the original SA model but most of them still overpredict the separation region.Given the effectiveness of QCR within SA model, SST-QCR model is expected to improve the prediction of compressor corner separation in this study.

In this paper,Menter’s SST model was chosen as the benchmark turbulence model to validate the effectiveness of helicity correction and QCR in predicting compressor flow.Detailed flow field analysis was conducted to determine the physical reasons for the improved modelling accuracy while revealing modelling defects that still exist.

2.Numerical methods

2.1.Flow solver

Numerical simulations were conducted by using ANSYS Fluent 19.2, where model corrections were implemented via userdefined macros.The default pressure–velocity coupling algorithm was chosen for solving the flow governing equations and turbulence equations in a strongly coupled manner.The second-order upwind scheme is used for spatial discretization of the advection terms of all transport equations.The default second-order upwind scheme settings in ANSYS Fluent are the ‘‘Least Squares Cell Based” method to evaluate cellcentered gradients and‘‘TVD-Minmod function”gradient limiter to clip the reconstructed solution oscillations on cell faces.While choosing the spatial discretization scheme, the pseudo transient option is activated to accelerate solution convergence without compromising numerical stability.

2.2.Turbulence models

In this section, Menter’s SST model and its variants are presented.The seven-equation full Reynolds stress model(Stress-BSL) and its potential advantages in predicting complex flows are also outlined.Due to limitation of the article length, only key variables and terms are explained.

2.2.1.Menter’s SST 2003 model (SST)

In Eq.(1)and Eq.(2),the first three source terms represent the production, dissipation, and diffusion of TKE k and specific dissipation rate ω, respectively.The fourth term in Eq.(2) is the cross-diffusion term generated by transforming the k-ε model to its k-ω formulation, where F1is a blending function that blends k-ε model and k-ω model in the boundary-layer wake region.F1equals one in the log layer and below, where the original k-ω model recovers; from the log layer up to the boundary layer edge, F1gradually reduces to zero, indicating the region being covered by both k-ε and k-ω model; in the freestream, F1equals zero, and the model becomes the k-ω form of k-ε model.

where the TKE production termis defined as30

The kinematic eddy viscosity is given as30

where the second term in bracket is the eddy-viscosity limiter used to limit the growth of eddy viscosity when encountering the boundary layer region of strong adverse pressure gradients.This limiter is particularly useful when predicting the shockinduced boundary layer separation as the original definition(k/ω) overpredicts νtin the region immediately after shock and leads to premature boundary layer reattachment and thus smaller separation bubble.F2is a blending function, which remains as one in most regions of boundary layer and rapidly decreases to zero close to the boundary layer edge.This means that F2only serves to limit the growth of eddy viscosity inside the boundary layer.

2.2.2.SST with velocity helicity correction (SST-Helicity)

2.2.3.SST with QCR (SST-QCR)

In SST-QCR model, Reynolds stress anisotropy is partially accounted by adding higher-order strain/vorticity terms to the linear Reynolds stress–strain constitutive relation.The modified constitutive relation induces counter-rotating streamwise vortices at the corner in ducts,31,32which entrains fluids and inhibits unphysically large corner separation.In this paper, the QCR model is regarded as the corner flow correction, and it is expressed as

As shown in Eq.(7),the first part(the first two terms on the right-hand side) constitutes the Boussinesq hypothesis while the second part (the third term on the right-hand side) represents the anisotropic part of Reynolds stresses.The cellcentered Reynolds stresses on the left & right cells sharing a common face are mass-weighted averaged on the cell face.The face-centered Reynolds stresses are then substituted into the diffusion terms of momentum equations to close the flow governing equations.Ccorneris called anisotropy coefficient that represents the strength of turbulence anisotropy, with the default value being 1.0.

In addition, as seen in Eq.(3), the change in the Reynolds stresses changes the production of TKE, which indicates that the turbulence non-equilibrium may be induced by Eq.(7).

2.2.4.ω-based Reynolds stress model (Stress-BSL)

The Reynolds Stress Model (RSM) is known to date the most elaborate RANS turbulence model.33Abandoning the Boussinesq hypothesis,RSM solves the six Reynolds stresses directly to close the flow governing equations.It contains six Reynolds stress transport equations, together with an equation for the dissipation rate.Thus,for 3D flows,there are seven turbulence transport equations to be solved to obtain Reynolds stresses at the end of each iteration or time step.The exact Reynolds stress transport equation is shown below for relevant analysis.

where the seven terms on the right-hand side of Eq.(8)denote turbulent diffusion,laminar diffusion,stress production,buoyancy production, pressure strain, dissipation, and production by system rotation, respectively.Of the various terms in Eq.(8), DL,ij, Pijand Fijdo not require modelling while DT,ij,Gij,φijand εijneed to be modelled to close Eq.(8).The differences between Stress-BSL model and other types of RSM models exist in the modelling of the pressure-strain term φijand the dissipation term εij.In Stress-BSL,the dissipation term is solved by the Baseline(BSL)ω-scale equation.This gives the Stress-BSL model advantages over other RSM models by removal of sensitivity to freestream turbulence variable ω and avoidance of the use of wall functions.

Since the effects of streamline curvature, rotation, swirl,and sudden changes in the mean strain rate are mathematically taken into account,i.e.,effects contained in the seven terms on the right-hand side of Eq.(8), RSM models have greater potential to give accurate predictions for complex flows compared with the simpler one or two-equation models.However,the fidelity of RSM is still limited by the assumptions made to model the pressure-strain term and the dissipation term.This means that compared with simpler models, RSMs may not always yield improved results in all classes of flows that are worthy of the additional computational cost.However, for flow features which are the result of anisotropy in Reynolds stresses, e.g., stress-induced vortices at duct corner, the use of RSM is often necessary.

In this paper, the Stress-BSL model is employed to predict Rotor 67 flow.The prediction results are used for model comparison and flow field analysis.

3.Cases description

3.1.NASA Rotor 67

The first test case is NASA fan Rotor 67, which is the firststage rotor of a transonic two-stage fan.34Detailed test rig setup and flow field measurements were conducted by Strazisar et al.34,and the relevant parameters and measurement data can be found in their paper.34In this paper, model validation works include the global performance at full rotating speed,the radial distribution of aerodynamic parameters (i.e., absolute total pressure,total temperature,and circumferential flow angle) at the outlet measurement plane; flow field analysis includes the blade surface limiting streamlines, the near-tip vortex flow field and eddy viscosity field.

The flow domain along with the measurement planes is shown in Fig.1.In numerical practices, one single passage was set as the computational domain for steady RANS simulation,with the blade-to-blade interfaces set as rotational periodic boundaries.The performance characteristics were obtained based on mass-weighted averaged total pressure and total temperature at inlet measurement plane and outlet measurement plane.In the experiment, the total pressure was work-weighted averaged while the total temperature was mass-weighted averaged.Since the difference between the mass-weighted average and the work-weighted average is quite small (the relative error of total pressure ratio being less than 0.01% in this case), the mass-weighted average is used throughout this case.The total pressure, total temperature and flow angle were extracted from five-hole probes located at both measurement planes, and then circumferentially mass-weighted averaged to obtain radial profiles.For reasonable comparative analysis, the CFD data were extracted and post-processed in the same way as the experiment.

As for the inlet mean-flow conditions,a uniform total pressure and total temperature are set to be equal to the values of standard sea level atmosphere (101325 Pa and 288.15 K).The inlet flow direction is axial.

As for the inlet turbulence conditions,the freestream turbulence intensity is set as 1.5%, consistent with the measured value.At the outlet, the static pressure is given.Starting from a low pressure(60000 Pa under mass-flow choking condition),a gradually increased pressure is set to obtain rotor performance under each operating condition.When approaching stall, a step increase of 0.2 kPa is adopted, and the last converged point is taken as the numerical stall point.

An O4H grid topology was used, and hexahedra structural cells were generated using NUMECA Autogrid.In the spanwise direction, 81 mesh layers were arranged.The distance from the center of cells adjacent to the endwall to the endwall is set as 10–5m,ensuring that most of y+are less than 2.0.This ensures turbulence model’s resolution of boundary layer down to the endwall.Hence, no wall function is needed and used in this case.On the S1 surface,there are 209 nodes in the streamwise direction in which 109 nodes are inside the blade passage,65 nodes across one blade pitch and 21 nodes along the bladenormal direction within the skin O block.The O block mesh is clustered onto blade to ensure that most of y+on the cell centers adjacent to the blade are less than 2.0.

Fig.1 Flow domain for NASA Rotor 67.

As for rotor tip clearance and hub fillet,the dimensions are set as 0.61 mm (tip gap size) and 1.78 mm (hub fillet radius),respectively.For tip gap meshing, as shown in Fig.2(a), the HO grid topology was used, with 109 nodes in the streamwise direction, 13 nodes in clearance O mesh and 21 nodes in the spanwise direction.There are overall about 26000 cells within tip gap, which is considered fine enough to resolve tip leakage flow in RANS framework.For hub fillet meshing,as shown in Fig.2(b), the simple H grid topology was adopted, and 17 nodes were placed along the spanwise direction to ensure sufficient mesh resolution for hub fillet.

Overall, the computational domain consists of around 1.52 million cells.The computational mesh has proven to be fine enough to satisfy grid independence.35

3.2.Darmstadt’s 1.5-stage axial compressor

Fig.2 Grid topology of tip clearance/hub secondary geometry(Rotor 67).

The second test case is Darmstadt’s 1.5-stage transonic axial compressor.At 2020 Global Power and Propulsion Chania Technical Conference (GPPS Chania20), the Institute of Gas Turbines and Aerospace Propulsion (GLR) of TU Darmstadt University disclosed the geometry and part of experimental data of the compressor and named it as TUDa-open-stage.36TUDa-open-stage consists of a blisk rotor with 16 blades, a stator with 29 optimized blades, and 5 Outlet Guided Vanes(OGVs) used to guide flow to the axial direction.Since 1994,a series of experimental works have been conducted by GLR to investigate its steady-state performance, unsteady aerodynamics (stall inception mechanism, pre-stall disturbances)and aeroelastic phenomena (non-synchronous vibration,forced response).The quantities used for comparison in this paper include the performance characteristics at full rotating speed,the radial profiles of the absolute total pressure and flow angle at rotor exit (ME21), the radial profiles of the absolute total pressure and total temperature at stator exit (ME30),the rotor/stator surface limiting streamline and the near-tip flow field.

The flow domain along with the measurement planes(ME21&ME30) and illustration of test rig are shown in Fig.3 and Fig.437, respectively.In numerical practices, one rotor passage, one stator passage and one OGV passage were set as the computational domain for steady RANS simulation,with the blade-to-blade interfaces set as rotational periodic.In the experiment, the mass flow was corrected to the standard sea-level condition, and the total pressure ratio and total temperature ratio were calculated based on the area-weighted averaged data.36In the CFD calculations, the same areaaveraged method was used.One thing to note is that in the experiment, the isentropic efficiency was calculated based on the area-averaged total pressure ratio and measured shaft power, rather than the area-averaged total temperature ratio.This is due to the relatively lower accuracy of area-averaged total temperature ratio in this case.36Thus in the CFD validation, a combination of area-averaged total pressure ratio and mass-averaged total temperature ratio was recommended for use to calculate isentropic efficiency.36At the rotor exit(ME21), the total pressure and flow angle were measured by five-hole probes.Both were circumferentially area-averaged at the corresponding axial & radial location to obtain radial One-Dimensional (1D) profiles.At the stage exit (ME30), the total pressure and total temperature were obtained by the rakes distributed along the circumferential direction.The obtained total pressure and total temperature were circumferentially area-averaged to obtain the radial profiles.

For the inlet mean-flow conditions (at ME15), the experimental total pressure and total temperature over different operating points were averaged,with averaged profiles of both specified as the numerical inlet conditions.The averaged profiles of total pressure and total temperature at the inlet are shown in Fig.5.It can be seen that the boundary layer at the casing is captured by the rakes placed at ME15,and a uniform total temperature is given by thermocouple at the inlet.The inlet flow direction is axial.

Fig.3 Flow domain for TUDa-open-stage.

Fig.4 Illustration of TUDa-open-stage test rig.37

For the inlet turbulence conditions,a freestream turbulence intensity of 4% and turbulent length scale of 0.09 mm are set at ME15,consistent with the measured results.37At stage exit,a uniform static pressure is given.Starting from low back pressure(97295 Pa under choked condition),a gradually increased pressure is set to obtain the stage performance at each operating point.When approaching stall,a step increase of 0.2 kPa is employed,and the last converged point is taken as the numerical stall point.

Fig.5 Experimental and numerical inlet boundary conditions.

An O4H grid topology was used, and hexahedra structural cells were generated using NUMECA Autogrid.In the spanwise direction, 93 mesh layers were arranged.The distance from the center of the first cell layer to the endwall is set as 3×10-6m, ensuring that most of y+are less than 2.0.This ensures that turbulence model can resolve the entire endwall boundary layer.Hence, no wall function is needed and used in this case.On the S1 surface, for rotor, there are 129 nodes in the streamwise direction in which 97 nodes are inside the blade passage, 65 nodes across one blade pitch and 21 nodes along the blade-normal direction within the skin O block.The O block mesh is clustered onto blade to ensure that most of y+on the cell centers adjacent to the blade are less than 2.0.For stator,there are 137 nodes on the blade suction surface,81 nodes on the pressure surface,73 nodes across one blade pitch and 21 nodes along the blade-normal direction within O block.Again, O block mesh is clustered onto blade to ensure that most of y+on the cell centers adjacent to the blade are less than 2.0.For OGV, it is meshed with 0.8 million hexahedra grid cells.Since OGV has a negligible effect on the upstream flow field, the grid density is considered enough.

For rotor, the tip clearance is set to 0.75 mm to meet the average value measured at full speed.37The rotor hub fillet radius is set as 5 mm, a constant value based on the open experimental database.For tip gap meshing, as shown in Fig.6(a), the HO grid topology was used, with 97 nodes in the streamwise direction, 17 nodes in clearance O mesh and 21 nodes in the spanwise direction.There are about 30000 grid cells distributed within tip gap,which is considered fine enough to resolve the tip leakage flow in RANS framework.For hub fillet meshing, as shown in Fig.6(b), the H grid topology was adopted,and 16 cell layers were placed to ensure sufficient mesh resolution for hub fillet.

For stator, the radius of the hub and casing fillets is set according to the shape of the stator geometry.There are 16 spanwise mesh layers covering hub and casing fillets to ensure sufficient mesh resolution for these geometric details.

Overall, the one-passage domain consists of around 2.71 million grid elements, of which 1.56 million are for rotor,1.15 million for stator,and 0.8 million for OGV.In the present study, the grid independence test has been done for this case.The test results show that the above grid satisfies the grid independence requirements.

4.Results and discussion

4.1.NASA Rotor 67

The calculated choking mass flow of Rotor 67 at the full design speed(16043 r/min)is shown in Table 1.The calculated values by the four turbulence models are all around 34.67 kg/s,while the measured value is 34.96 kg/s.The discrepancy between them is less than 1%.This demonstrates that the selected turbulence models give high predicted accuracy for the choking flow condition of Rotor 67.

Fig.6 Grid topology of tip clearance/hub secondary geometry(TUDa-open-stage rotor).

Table 1 Mass flow rate (choking condition).

Fig.7 shows the predicted performance characteristics of Rotor 67.The mass flow rates (m) are nondimensionalized by their respective choking mass flow rate(mc).Both the calculated total pressure ratio characteristics and the isentropic efficiency characteristics are in good agreement with the measurements.The stall margin predicted by all four models is slightly smaller than that of experiment.Under near-stall conditions, the total pressure ratio is slightly underpredicted.

To further validate the models’ predicted performance,Fig.8 shows the radial distributions of the circumferentially mass-averaged total pressure, total temperature, and flow angle at the outlet measurement plane.Under Peak Efficiency(PE)condition(=0?989±0?004),the predicted distributions are in good agreement with the measured ones.However, it can also be seen that at midspan, the calculated yaw angle overshoots the measured yaw angle.Arnone38and He et al.39also gave similar results and attributed these to the lower mass flow rate predicted(see Table 2)while the predicted total pressure and total temperature remain almost the same as the measured ones (see around midspan in Fig.8(a)).From the perspective of flow field,the larger absolute flow angle at midspan is caused by the joint action of two factors.One is lower mass flow rate predicted under PE condition, compared with the measured result; the other is overpredicted suction surface boundary layer blockage by CFD.According to the velocity triangle rule,a combination of these two factors leads to larger absolute yaw angle at midspan.

Fig.7 Characteristics of Rotor 67 at design speed.

Furthermore, under PE condition, the total pressure and total temperature distributions predicted by Stress-BSL model near the casing (normalized span between 0.8 and 1.0) are in good agreement with measurements, while the distribution trend given by Stress-BSL model ranging from 0 to 0.2 span is opposite to that measured in experiment.In contrast, the total pressure and total temperature in the upper span are underpredicted by SST-type models.This indicates the total pressure loss predicted by SST-type models being higher, and the work done on the airflow near the blade tip being relatively smaller.In the lower span ranging from 0 to 0.3, although the distribution trend of both total pressure and total temperature are reasonably predicted, the SST-type models all predict lower total pressure and total temperature values.This indicates the larger corner separation being predicted by SSTtype models near the hub.Overall,it can be shown that under the PE condition,Stress-BSL model performs better than SSTtype models in terms of the total pressure and total temperature predictions in the upper blade span while performing the worst near the hub where intense secondary flow occurs.

Under Near Stall (NS) condition=0?924±0?004), as shown in Fig.8(b), the SST-Helicity model performs the best in predicting the radial distributions of total pressure and total temperature, and the predicted results are in good agreement with the experimental measurements.This demonstrates SST-Helicity model’s superiority in predicting complex flows in rotor under the NS conditions.The physical reasons for the advantages of helicity correction under off-design operating conditions are revealed by analyzing the flow field near both the blade tip and the blade surface (see Fig.9 and Fig.10).

Fig.8 Spanwise distributions of total pressure, total temperature, and yaw angle (α) at outlet measurement plane.

Table 2 Mass flow rate (PE condition).

As shown in Fig.9, the effect of the helicity correction on the growth of the tip leakage flow is revealed by comparing the velocity helicity field and the turbulent eddy viscosity field near the blade tip.Eight slices are arranged axially from the blade leading edge to trailing edge, which are numbered from 1 to 8 (see Fig.9(c)).

From Fig.9(a), it can be shown that in the vicinity of casing, high relative velocity helicity exists, indicating strong three-dimensionality of the casing boundary layer.In the vortex core region,due to vortex stretching,the angle between the velocity vector and the vorticity vector tends towards zero.The absolute value of normalized velocity helicity is close to 1.0.This means that the normalized velocity helicity can be used to identify local concentrated vortex cores.It can also be used to represent vortex stretching and its effect.That is, the larger the normalized velocity helicity value, the more intense the vortex stretching, and the stronger the entrainment of the tip leakage vortex on the surrounding mainstream fluid.In the process of fluid entrainment, intense mixing takes place in between the tip leakage vortex fluid and the mainstream fluid,accompanied by energy transfer from the mainstream to the tip leakage flow.This increases the TKE of tip leakage flow,as evidenced by the near-tip vortex eddy viscosity field shown in Fig.9(d).

By comparison between Fig.9(b) and Fig.9(c), it can be found that the tip leakage vortex predicted by SST and SSTQCR migrates faster to the pressure surface of the adjacent blade.This results in a larger tip-leakage blockage, which in turn reduces the total pressure in the near-tip region downstream of the trailing edge (see Fig.8(b)).In contrast, SSTHelicity and Stress-BSL more reasonably capture the migration path of tip leakage vortex.This can be demonstrated by comparing the location of the vortex core on Slice 5 predicted by the four models with the location of the vortex core on Slice 5 in Fig.9(a).The more reasonable predictions of tip leakage blockage by SST-Helicity model and Stress-BSL model make their predictions in terms of total pressure and exit flow angle more consistent with the experiment.

Fig.9 Near-tip vortex flow field.

Fig.10 Blade suction surface limiting streamlines of Rotor 67 (NS condition).

Moving from the near-tip region downwards, the effect of velocity helicity correction on the flow field is shown in Fig.10.Compared with SST and SST-QCR, the shockinduced separated shear layer predicted by SST-Helicity and Stress-BSL reattaches earlier,the separation bubble is smaller,and the shock is pushed further downwards.This constitutes one reason for SST-Helicity and Stress-BSL models’more reasonable prediction of total pressure distributions from 0.6 span upwards (see Fig.8(b)).In addition, compared with SSTHelicity model,SST,SST-QCR and Stress-BSL models all predict significantly larger corner separation, and thus more lowenergy fluid within corner separation moves radially upwards.This results in a more pronounced total pressure deficit area located at higher span,as shown from the total pressure distributions by SST, SST-QCR and Stress-BSL models around mid-span in Fig.8(b).

Fig.11 Streamwise vorticity contours around hub corner of Rotor 67 (NS condition).

Furthermore, by comparison between SST and SST-QCR results,it is found that QCR does not improve the corner separation prediction in this flow.Since QCR has proven to be able to improve corner flows in the simple rectangular ducts,31it is interesting to investigate the reason for QCR’s ineffectiveness for the flow in this study.

Fig.11 shows the distribution of streamwise vortices around the corner immediately upstream of the corner separation region.The location of the corner separation onset can be identified by the surface limiting streamlines shown in Fig.11.The streamwise vortex distribution predicted by SST-QCR is almost the same as that by SST.The Reynolds stress anisotropy due to QCR correction fails to induce the counterrotating streamwise vortex pairs at the corner.This is in direct contrast to the performance of QCR on the simple rectangular duct flow.In the rectangular duct flow, QCR successfully induces the counter-rotating streamwise vortex pairs at the duct corner.31The induced streamwise vortices entrain highmomentum fluid from the mainstream into the corner,energizing corner flow to be resistant to separation.Conversely, the failure of QCR to induce counter-rotating corner vortices accounts for the QCR’s ineffectiveness to inhibit the unphysically large corner separation.This study reveals the limitation of QCR in capturing the complex flow physics around the corner.One possible way to improve QCR is to use high-fidelity LES database to construct the Reynolds stress–strain constitutive relation, so as to accurately capture the anisotropic characteristics of corner flow.

In summary, compared with SST and SST-QCR, SSTHelicity gives predictions of the rotor performance and related flow field with higher accuracy.The success of SST-Helicity is attributed to the helicity correction reasonably capturing the non-equilibrium turbulence behavior of tip leakage and corner separation flows,i.e.,the energy transfer between them and the main flow.Through energy transfer with the main flow, both flows are energy equipped to be more resistant to the bulk adverse pressure gradients, thus improving the prediction accuracy under near-stall condition.In contrast, Stress-BSL does not always yield higher accuracy than the simpler twoequation eddy viscosity models.It performs even worse than SST-Helicity under the near-stall condition.Considering the much higher computational cost (solving seven turbulence equations in a coupled manner), numerical robustness, and modelling accuracy,Stress-BSL model is not recommended for routine compressor performance prediction and analysis.

Next, the SST model, SST-Helicity model and SST-QCR model are employed to validate the aerodynamic performance of a high-speed axial compressor.The purpose is to validate the capability of SST-Helicity to predict the aerodynamic performance of the compressor stage and related flow field.

4.2.Darmstadt’s 1.5-stage transonic compressor

Fig.12 shows the prediction results of the TUDa-open-stage performance characteristic lines.It can be shown that the choking mass flow rate predicted by the three models is very close (around 16.345 kg/s) and is 0.13 kg/s higher than the measurement value(16.213 kg/s).The relative error of choking mass flow rate is 0.81%.This shows that the original SST model and improved ones possess high accuracy in predicting the choking flow condition.Helicity correction and QCR correction have almost no effect on the results under this condition.In addition, all turbulence models reasonably predict the trend of the total pressure ratio.Among them, SST model prematurely predicts the stall point,the stall margin being too small, while the improved SST models significantly improve the prediction accuracy of the position of stall point.Among them, the stall margin predicted by SST-Helicity is more than double that by SST and is the closest to the actual stall margin.

Fig.12 Characteristics of TUDa-open-stage at design speed.

To further validate models’ performance, Table 3 and Table 4 show the prediction results of the global aerodynamic parameters under the PE condition(mcorr=16?00±0?02 kg/s)and NS condition (mcorr=14?78±0?02 kg/s), respectively.Under the PE condition, the predicted total pressure ratio is consistent with the experimental results.The predicted isentropic efficiency is about 1% higher than the experimental value.Here the higher isentropic efficiency indicates the predicted stage loss being smaller than the actual loss.Under the NS condition, the SST-Helicity model shows significant advantages in predicting the stall-point performance and the stable operating range.The stall margin given by SSTHelicity model more than doubles that by SST model and is the closest to the experimental stall margin.As for the isentropic efficiency, the experimental measurement deviates more from the predicted result as the NS condition approaches,indicating that the actual loss is much larger compared to the predicted loss.There shall be sources of loss existing in rig tests that are beyond the reach of current CFD tests.Detailed investigation of the loss sources will be presented in the later part.

Fig.13 shows the spanwise distribution of total pressure and circumferential flow angle at rotor outlet (ME21) under PE condition.It can be seen that near the casing (above 70% of the rotor span), the SST model and its variants overpredict the total pressure ratio, with the prediction results almost overlapping.This indicates the deficiency of SST model in predicting tip leakage flow or near-tip geometric uncertainty, while velocity helicity correction and QCR correction still fail to improve the predictions.Between 20% and 80%span,the absolute flow angle predicted by the models is somewhat larger than the experimental values.This suggests that the predicted boundary layer blockage is smaller than the actual blockage(‘‘thinner”blade suction side boundary layer).

Fig.14 shows the spanwise distribution of total pressure and total temperature at stage outlet(ME30)under PE condition.It can be seen that above 80%of stator span,a somewhat higher total pressure ratio is given by SST model and its variants.The near-tip total pressure difference originates from the rotor and is propagated from the near-tip region in the rotor passage.Near the hub(below 20%of stator span),SST model and its variants overpredict the total pressure.The prediction deviation is believed to be caused by two factors.One is the modelling error of the stator corner separation by turbulence models; the other is the stator hub leakage flow through the secondary passage below stator hub (see Fig.4), which mixes with the main flow in the stator passage and causes additional endwall losses.Since both the helicity correction and QCR have been demonstrated to be effective in predicting the corner separations,26,32the author is more inclined to the second factor,that is,the neglection of the stator hub leakage flow in the current CFD simulation.

In order to more intuitively determine the difference near both the stator endwalls, the predicted and measured total pressure distributions are shown in Fig.15.Near the casing,the predicted total pressure in the passage vortex core area is smaller than the measured.The predicted low-pressure area is more compact and concentrated.This indicates that RANS underestimates the mixing between the passage vortex and the mainstream as the passage vortex develops downstream.Near the hub, the three models give similar corner separation area,with extent of which being smaller than the measured.This indicates that the stator hub leakage flow is one probable source that leads to the hub corner difference.The secondary flow goes through the stator hub cavity and leaks from the hub between rotor and stator rows.As it develops downstream, the mixing between the leakage flow and the endwall secondary flow in the stator passage leads to an increase in the total pressure loss near the hub endwall.However, thenumerical simulations did not account for the secondary flow path beneath the stator hub, the stator hub leakage flow and its effect.The author wishes that in the future GLR would disclose the geometry details of the secondary path and the measurement data associated with the stator hub leakage flow.

Table 3 Aerodynamic parameter prediction results (PE condition).

Table 4 Aerodynamic parameter prediction results (NS condition).

Fig.13 Spanwise distribution of circumferentially area-averaged total pressure and yaw angle at ME21 (PE condition).

Fig.14 Spanwise distribution of circumferentially area-averaged total pressure and total temperature at ME30 (PE condition).

Next, the predictive performance of SST-Helicity model and SST-QCR model under the NS condition is demonstrated.The SST results are in absence here because the calculation using SST model undergoes divergence at mass flow rate higher than 14.78 kg/s (see Fig.12).Since GLR does not disclose the 1D radial distribution of total pressure and circumferential flow angle at ME21 under NS condition,Fig.16 only shows the 1D spanwise distribution of total pressure and total temperature at ME30 under NS condition.

Fig.15 Total pressure contours at ME30 (PE condition).

From Fig.16, it can be seen that near the casing (above 70% of stator span), the variation trend of the total pressure predicted by two models is similar to the experimental results.The SST-Helicity model overpredicts total pressure while SSTQCR model underpredicts it.The predicted deviation of the total pressure near casing is attributed to the predicted deviation of the blockage due to separation near the stator casing.Fig.17 shows the particle traces under the NS condition.It can be shown that the flow separates locally and the SSTQCR model gives much larger separated flow area.The vortices given by the SST-QCR model induce large reverse flow and thus blockage in the stator passage.In contrast, no big vortex is given by SST-Helicity model,indicating much smaller amount of reverse flow and blockage to the main flow.

Below the mid-span,the deviation of the total pressure predictions from the experiment becomes more significant compared to the PE condition (see Fig.14).The higher predicted total pressure combined with the lower predicted total temperature results in significantly higher isentropic efficiency (see Fig.12(b)).This is because under the NS condition, due to the higher back pressure imposed, more low-energy endwall fluid enters the hub gap downstream of the stator trailing edge,creating a larger amount of stator hub leakage flow.The stator hub leakage flow re-enters into the main flow path from the inter-blade row hub gap and mixes with the main flow, generating more endwall losses in the stator passage.

Fig.18 shows the effect of the predicted blockage on the downstream total pressure distribution.Compared with PE condition, the casing corner separation under NS condition becomes more intensive (see Fig.15 for comparison).The low total pressure area (dark blue region) is caused by the accumulation of low-energy fluid in the corner separation region.Compared with SST-QCR model, the total pressure distribution near the casing given by SST-Helicity model is more consistent with the experiment, despite the predicted low-pressure core being compact and concentrated.Near the hub,a large low-pressure area can be seen from the experimental results,which cannot be found in the CFD results.The reason is as what has been observed, that is, the stator hub leakage flow and its influence.

In summary, the SST-Helicity model has been demonstrated to be able to improve the prediction accuracy of compressor flow and thus performance.This is mainly reflected in the extension of the stall margin to be close to the real stall margin.Further analysis of the flow field reveals that the improved accuracy is largely due to its more reasonable prediction of the corner separation near the stator casing and the resulting blockage.The validation results of this case show that the commonly used two-equation turbulence models,when coupling with the velocity helicity correction, can significantly improve the prediction accuracy of compressor flow.

Fig.16 Spanwise distribution of circumferentially area-averaged total pressure and total temperature at ME30 (NS condition).

Fig.17 Particle traces of TUDa-open-stage stator (NS condition).

Fig.18 Total pressure contours at ME30 (NS condition).

5.Conclusions

In this study,SST model and its variants(SST-Helicity and SSTQCR)have been validated for the compressor rotor flow and the compressor stage flow.Detailed flow field analysis is also carried out to demonstrate the effectiveness of model corrections in capturing the complex flow features in compressors.The conclusions obtained from this study are drawn as follows:

(1) In the vortex core region, the velocity vector and the vorticity vector tend to parallel, and the absolute value of normalized velocity helicity tends to be 1.0.The modification of turbulence production source terms by using normalized velocity helicity better models the turbulent energy transfer from the main flow to 3D vortical flows(e.g., rotor tip leakage flow and corner separation flow in compressors).The 3D vortical flows are energy equipped to be resistant to diffusion and passage shock.The SST-Helicity model is regarded as one of the nonequilibrium turbulence models.

(2) For Rotor 67 validation case, all the turbulence models investigated show quite good predictive performance.Particularly under the near-stall operating condition,compared with SST model, SST-QCR model and Stress-BSL model, SST-Helicity model shows its significant advantage in predicting 3D vortical flow and its impact on the rotor performance.One of the manifestations is that the non-physical large corner separation is inhibited, which in turn weakens the radial migration of low-energy fluid in the separation zone.The radial distribution of the outlet total pressure is improved and is closer to reality.Another manifestation is the earlier reattachment of the shock-induced separated shear layer, resulting in smaller separation bubble over the entire blade span.Hence, the unphysically large total pressure loss induced by shock is avoided.

(3) For Darmstadt’s 1.5-stage axial compressor case, the advantages of the SST-Helicity model are more obvious,in terms of predicting stall margin and the near-stall performance.At the design speed, the stall margin given by the SST-Helicity model (20.90%) is the closest to the experimental measurement result(24.81%),which is more than double that by the original SST model(8.71%)and 1.5 times that by the SST-QCR model(14.14%).This success is attributed to the SST-Helicity model’s reasonable prediction of the growth of corner separation near the stator casing,where the large corner separation as given by SST model and SST-QCR model is avoided.The improved prediction of the blockage due to stator corner separation contributes to numerical convergence,thereby extending the predicted stable operating range.

(4) For the 1.5-stage axial compressor,under both its peakefficiency and near-stall conditions, one common difference between CFD and experiment exists in the lowpressure region around the stator hub corner.This is due to the existence of the leakage flow beneath the stator hub in the test rig.There can be a substantial effect on the blade row loss coefficient associated with even a small amount of leakage flow while its existence and interaction with mainstream flow are not considered in current CFD simulations.In the future, one research direction is to collaborate with GLR to obtain more test rig details and the measured data associated with the secondary flow, so as to finely validate the capability of CFD for simulating real compressor geometries.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The author would like to thank Dr.Jiahuan CUI and Prof.Dingxi WANG for their technical advice during paper drafting.