# Detailed Unsteady Simulation of a Counterrotating Aspirated Compressor with a Focus on the Aspiration Slot and Plenum

Authors: Robert D. Knapke; Mark G. Turner
Journal: International Journal of Rotating Machinery
 (2013)

---

## Abstract

An unsteady analysis of the MIT counterrotating aspirated compressor (CRAC) has been conducted using the Numeca FINE/Turbo 3D viscous turbulent flow solver with the Nonlinear Harmonic (NLH) method. All three blade rows plus the aspiration slot and plenum were included in the computational domain. Both adiabatic and isothermal solid wall boundary conditions were applied and simulations with and without aspiration were completed. The aspirated isothermal boundary condition solutions provide the most accurate representation of the trends produced by the experiment, particularly at the endwalls. These simulations provide significant insight into the flow physics of the aspiration flow path. Time histories and spanwise distributions of flow properties in the aspiration slot and plenum present a flow field with significant temporal and spatial variations. In addition, the results provide an understanding of the aspiration flow path choking mechanism that was previously not well understood and is consistent with experimental results. The slot and plenum had been designed to aspirate 1% of the flow path mass flow, whereas the experiment and simulations show that it chokes at about 0.5% mass flow.


---

## Body

## 1. Introduction

The overall goal during the design process of a gas turbine engine is to reduce specific fuel consumption for a given level of thrust. Often great lengths are taken to obtain even small gains in engine efficiency. In the field of aviation, weight reduction is of equal importance as the weight of the engine affects fuel consumption. Two compressor design concepts for improving engine performance are counterrotation of neighboring stages and aspiration along the suction surface of the compressor blades.Counterrotation allows for a pressure ratio similar to two stages without the intermediate stator, thus decreasing weight. Early investigations of counterrotating designs were conducted by Curtiss-Wright in 1950 as discussed by Wennerstrom [1]. Although this early effort failed to reach design goals, improved aerodynamic design tools and the use of aspiration allowed the MIT Gas Turbine Laboratory (GTL) to develop more fruitful designs as discussed by Freedman [2], Kirchner [3], and Merchant et al. [4]. Two drawbacks of counterrotating designs are high relative mach numbers (Mrel) and high loading. The use of aspiration can mitigate these detrimental effects.The effectiveness of aspiration on highly loaded stators was initially investigated by Loughery Jr. et al. [7], which showed that aspiration improved performance. As discussed by Kerrebrock [8], aspiration removes the high entropy boundary layer fluid, which allows for higher loading at a given blade speed. Numerous aspirated compressor designs and experiments were conducted at the MIT GTL by Reijnen [9], Ziminsky [10], Merchant et al. [11–14], and Schuler et al. [15–17]. These studies showed potential for significant performance benefits using small amounts of aspiration through careful placement of the aspiration location. In addition, with transonic compressor designs the use of aspiration moves the shock downstream, thus reducing shock losses.The details of the CRAC design and testing are described by Kerrebrock et al. [5] and Parker [6]. The performance goals for the CRAC design were a pressure ratio of 3 and an adiabatic efficiency of 87%. The compressor includes an inlet guide vane (IGV) and two counterrotating rotors with blade counts of 35, 20, and 29, respectively. A parametric design was used to balance the Mrel and loading of the two rotors. A higher rotor one tip speed of 442.0 m/s was chosen to reduce the loading and a lower tip speed of 350.5 m/s was chosen for rotor two to reduce the Mrel. The rotor two Mrel was further reduced by designing the IGV to introduce a counter-swirl into rotor one. Aspiration of 1% of the inlet mass flow was included in the design of the second rotor to control both the shock position and the boundary layer growth. The aspiration slot is positioned at 60% chord and extends from the hub to approximately 80% span. The relevant parameters of the two rotors are listed in Table 1.Table 1
Aerodynamic parameters of rotor one and rotor two designs.
Rotor one
Rotor two
Tip speed (m/s)
442.0
350.5
Pressure ratio
1.92
1.60
Adiabatic efficiency (%)
91.2
86.6
AverageMrel
1.26
1.40The CRAC experiment was conducted in the MIT GTL blowdown test facility. The test section of that facility is depicted in Figure1(a). A mixture of CO2 and Argon was used as the working fluid to reduce the mechanical stresses on the blades while maintaining a ratio of specific heats (γ) similar to air. Before the test, this fluid was stored in a supply tank. The test section and the dump tank were evacuated and the rotors were accelerated to a set speed. Then, a fast-acting valve was opened and the working fluid flowed through the test section and into the dump tank. The aspirated fluid passes into a plenum and is removed through the root of the blade, as shown in Figure 1(b). Between 250 and 350 milliseconds, the axial and tip mach numbers were constant and the rotor corrected speeds were within acceptable limits to be representative of the given test point.Test section and aspiration flow paths.
(a)
Test section [5]
(b)
Rotor two aspiration flow path [5]Instrumentation of the CRAC rig consisted of an inlet plane just upstream of the IGV and a measurement plane downstream of rotor two. Both measurement locations consisted of 3 pairs of windows equally spaced around the annulus. Each of the pairs had total pressure, total temperature, and wall static pressure probes, and one window included a Pitot probe. Spanwise measurements were obtained using rakes of total pressure and total temperature at both the upstream and downstream locations. A variable area orifice was placed downstream of the aspiration plenum in an attempt to measure the aspirated mass flow, as depicted in Figure1(a). However, attempts to obtain the aspirated mass flow measurement failed. Further testing of the aspiration flow path showed that it choked at about 0.5% of the inlet mass flow [18]. Details related to the instrumentation are provided by Parker [6] and Onnee [19].The analysis in this paper provides insight into the flow physics of the CRAC compressor. The authors believe that this is the first detailed study of the unsteady characteristics of aspiration in a transonic, counterrotating compressor. The mechanism for the choking of the aspiration flow path is investigated, thus shedding light on details of the aspiration design that were not fully understood during the CRAC tests. Radial characteristics of the aspiration slot and unsteady interactions between the main passage and aspiration flow paths are discussed in detail. Through this analysis, a clear picture of the aspiration flow physics is realized.


## 2. Numerical Approach

The CRAC geometry was modeled with the Numeca International 3D viscous turbulent flow solver, FINE/Turbo [20, 21]. A finite volume 2 order accurate central scheme with artificial dissipation was used. Modeling turbulence was accomplished with the Spalart-Allmaras turbulence model. Unsteady blade row interactions were captured with the Nonlinear Harmonic (NLH) method, which is described by Vilmin et al. [22]. This and other comparisons between the NLH and time marching unsteady approaches have demonstrated the applicability of the NLH method to turbomachinery simulations [23, 24].The NLH method represents the solution variables as a summation of the time mean solution and a number of perturbations. The perturbations are periodic with a user-specified number of harmonics of the blade passing frequency. With this representation, the computational domain can be reduced to a single meshed passage per blade row because the solution can be reconstructed at any given time. In addition, the equations produced by the method are a series of steady state governing equations. Not only is this type of equation less computationally intensive to compute, but also acceleration methods, such as multigrid, can be used to speed up convergence.The computational efficiency of the NLH method is most evident when performing unsteady multistage turbomachinery simulations. With the traditional time marching approach, a periodic portion of the annulus must be included in the computational domain. For multistage simulations, this restriction often requires the inclusion of the entire annulus. Unsteady features of the CRAC geometry were captured with the NLH method while simulating a single passage in each blade row. Using a time accurate temporal scheme of a periodic domain would result in a computational domain nearly 30 times larger.


### 2.1. Computational Grid

The main passage grids were generated with the turbomachinery meshing tool AutoGrid. Constant axial inlet and outlet planes were positioned at the approximate locations of the upstream and downstream measurement planes. In addition, meridional control lines were placed along the aspiration slot to axially cluster the grids at that location. Figure2 depicts the meridional view of the computational domain. Previous experience with simulating the CRAC geometry by Knapke et al. [25] and grid dependency studies by Guidotti and Turner [26, 27] and Knapke [28] have shown that the results are sensitive to the spanwise grid resolution. For this reason, a spanwise cell count of 172 was chosen. An additional 12 spanwise cells were placed in the 10 mil rotor tip clearance gaps. A first cell wall spacing of 3×10-6 meters was applied to all blade rows which resulted in all y+ values below 5.Figure 2
AutoGrid meridional view (blades extruded beyond endwalls to ensure robust grid generation).An O4H block topology was used for each blade passage, which consists of an O-grid block around the blade and H-grid blocks at the inlet, outlet, upper, and lower sides of the passage. High stagger optimization and blunt trailing edge meshing features were used for the rotor passages. The trailing edge detail of both rotors is shown in Figure3. The blunt trailing edge was a design feature of the CRAC rotors. With the blunt trailing edge, the skin mesh is created with a C-grid. The tip clearance mesh for both rotors utilized an H-grid block surrounded by a C-grid block. A nonmatching boundary (NMB) is automatically applied to account for any mismatch in the number of cells between the suction and pressure sides. For NLH simulations, Numeca recommends the number of cells in the tangential direction be defined by the following equation:
(1)Ncellθ=30Nharmnadjn,
where Nharm is the number of harmonics and nadj and n are the adjacent and current blade counts, respectively. This guideline is established so the adjacent blade row has enough grid points (30) to define each harmonic mode. Following this guideline for Nharm=3, the tangential number of cells in the IGV, rotor one and rotor two blade rows were set to 72, 116 and 104, respectively. Small axial cell expansion ratios were maintained at the inlet and exit of each blade row to negate axial wave dissipation. The total cell counts are listed in Table 2.Table 2
Number of cells per grid.
Millions of cells
Total (aspirated)
17.53
Total (nonaspirated)
17.19
IGV
2.06
Rotor one
6.29
Rotor two
8.84
Aspiration slot and plenum
0.34Midspan rotor grids at the trailing edge.
(a)
Rotor one
(b)
Rotor twoIGG, a general structured multiblock 3D meshing tool, was used to generate the aspiration flow path grids. The aspiration flow path domains were held as close to the test geometry as possible, while making small modifications to simplify the grid generation. Special care was taken to maintain the areas of the slot and plenum. The area directly impacts the way in which the flow path chokes. The number of cells and cell spacings along the span of the slot were defined such that the spanwise grid densities were similar to the main passage grids. The width of the slot was resolved with 16 cells. Wall spacings of3×10-6 meters were applied throughout the aspiration flow path to maintain y+ values below 5. Decisions on the number of cells and clustering were driven by grid quality metrics, such as orthogonality, aspect ratio, and expansion ratio. Figure 4 shows the aspiration flow path grid and the clustering along the rotor two surface.Main passage clustering and aspiration flow path grid.
(a)
Rotor two blade surface grid
(b)
Aspiration slot and plenum grid


### 2.2. Boundary Conditions

The CRAC testing at MIT consisted of 14 test cases. Of these cases, Run 13 was chosen for two reasons. First, Run 13 was near 100% design speed for both rotors. Second, sufficient inlet and outlet data was available to provide accurate boundary conditions and ample data for comparison. From the range of acceptable data of the blowdown test, the 250 milliseconds data point was chosen. The details of Run 13 and the data used are presented by Parker [6].Total conditions with axial flow was applied at the inlet boundary. Inlet rake measurements for Run 13 were not available, however Run 10 data was used as an appropriate substitute given the similar mass flow and rotor speeds. A pressure screen was located upstream of the test section to improve the rig operability and provide uniform flow to the IGV (Figure1(a)). This screen mixes out the viscous boundary layer. Given the short distance between the pressure screen and the inlet rake and the low radial resolution of that rakes, the viscous boundary layer was not measured accurately. Therefore, TEXSTAN was used to calculate the viscous boundary layer that developed between the pressure screen and the computational inlet plane. TEXSTAN is a two-dimensional finite difference code for solving the boundary layer mass, momentum, and heat transfer equations [29]. Two annular duct TEXSTAN simulations were completed. One to simulate the total temperature from the supply tank to the test section and another to capture the viscous boundary layer between the screen and the test section. Both simulations assume constant isothermal walls at room temperature (293 K) and a uniform inlet flow boundary condition. A mixing-length model for the turbulent viscosity was used.The calculated total pressure profile is shown in Figure5(a). To check the accuracy of the TEXSTAN boundary layers, the total temperature profile is compared with the Run 10 test data in Figure 5(b). The TEXSTAN results agree well with the test data. Some difference occur near the hub, but that difference is likely the result of the heat produced by rotor one motor, which was not modeled in the TEXSTAN simulation. This is strong evidence that the blowdown rig is not adiabatic. The TEXSTAN total pressure and the test total temperature profiles were applied at the inlet boundary. An unintended consequence of the aforementioned pressure screen was that it produced inlet distortion, which resulted in a 4% variation of total pressure around the annulus, as described by Parker [6]. This distortion was not accounted for in the CRAC simulations because the distortion was not measured along the entire span.TEXSTAN inlet profiles and the Parker [6] inlet total temperature profile.
(a)
Total pressure profile
(b)
Total temperature profilesThe aspiration and main passage grids were combined in IGG. The connection between the slot and rotor two passage grids was made with a Full Nonmatching Boundary (FNMB). This FNMB was required because the slot grid boundary points do not correspond exactly to rotor two blade boundary points. Although a direct connection boundary would provide a more accurate transfer of flow details, the increase in complexity of a connecting topology greatly outweighed the benefits. Moreover, this topology would require modification of AutoGrid blocks in IGG. This type of modification is not recommended by Numeca as it removes the grid smoothing performed during the AutoGrid grid generation process.The FNMB uses interpolation to transfer information between the nonmatching boundaries. Although the grid cells are not coincident along the FNMB boundary, significant effort was made to have similar spacings on either side of the boundary. A poorly defined connection can result from cell spacings that vary greatly across the FNMB. Figure6 shows the results of a triangulation and projection procedure performed in IGG. This procedure visualizes the same connection process used by the solver.Main passage grid on the rotor two suction surface.
(a)
Full blade to slot FNMB
(b)
FNMB near hubAt the exit of the main passage, a simple radial equilibrium static pressure boundary condition was used. The design mass flow of 1% of the inlet mass flow was applied to aspiration plenum exit boundary and the aspiration flow path was allowed to choke, thus setting the mass flow. Attempts were made to switch to a static pressure boundary condition because it would more accurately represent the plenum. However, a low pressure was required to choke the flow, causing negative pressures and densities when switching from steady to NLH computations. Since the exit of the plenum is reasonably far from the location where the flow chokes, the mass flow exit boundary condition was deemed acceptable.The rotating speeds of the two rotors at the 250 millisecond data point were 12600 rpm and 10020 rpm. Simulations at these speeds resulted in an upstream shock position that is not consistent with the design intent. It was found that raising the rotating speeds to 103% of the experimental value moved the shock to the correct position, as depicted in Figure7. This increase of the rotor speed was deemed plausible for two reasons. The first reason is the potential uncertainties of the rotor twist during operation. Although the amount of twist was accounted for in the design, the actual rotor tip positions during the test were not measured. Second, the corrected speed, which sets the shock location, is dependent on temperature and R. During the experiment, the rotor could have been at a lower temperature than expected because of heat transfer. This difference would make the corrected speed higher. Additionally, variations in the gas mixture could affect R and further change the corrected speed.Static pressure contours at midspan of rotor two. Adiabatic three harmonic solutions.
(a)
100% rotor speed
(b)
103% rotor speedRecent research concerning the heat transfer boundary conditions at solid walls for compressor simulations, by Bruna and Turner [30] and Isomura et al. [31], has shown the adiabatic wall approximation to be inaccurate. The short duration of the CRAC blowdown test and the difference between the fluid and metal temperatures prior to the test results in more heat transfer than seen during more traditional test methods [6]. For these reasons both adiabatic and isothermal solid wall boundary conditions were simulated for the aspirated case. The surface temperatures were not measured during the experiment. Therefore, a conservative estimate of room temperature (293 K) was assumed for all surfaces for the isothermal simulation.At the periodic boundary, the time mean flow is treated as a traditional periodic boundary. For the harmonic perturbation equations, the periodic boundary condition comprises a phase-shift of the harmonic amplitudes related to the pitch of the given blade row and the perturbation periodicity. At the interface boundary, Fourier Transforms (FT) in both time and space are performed on both sides of the interface. The harmonics produced by the temporal decomposition are matched across the interface to the discrete number of perturbation frequencies. From this matching process and the spatial FT, the time mean and perturbation flows can be reconstructed. Continuity across the interface boundary cannot be exactly reached for a finite number of harmonics; however, increasing the number of harmonics reduces the discontinuity error.


### 2.3. Initialization

To improve convergence, steady mixing-plane solutions were used to initialize the NLH simulations. Steady computations were completed using full multigrid. The steady cases were initialized on a coarse grid from constant values and a lower back pressure. The back pressure and grid density were increased until the final grid and the Run 13 back pressure values were reached. This final steady solution was used to initialize an NLH simulation with one harmonic per blade row. For the adiabatic cases, the number of harmonics was incremented to three in successive simulations. The isothermal cases were converged with two harmonics.


### 2.4. Gas Model

The CO2 and Argon mixture used as the working fluid has a 44% higher molecular weight than air, which resulted in lower physical rotor speeds for a given corrected speed. The physical speed reduction in turn reduced the stresses in the rotating components of the compressor rig. For the Run 13 mixture, the molar percentages of CO2 and Argon were 48.91% and 51.09%. Over the temperature ranges of the experiment (250 K to 450 K), the mixture γ is similar to the γ of air. However, the higher molecular weight results in a lower gas constant (R), as shown in Table 3. Calculation of the mixture specific heat at constant pressure (Cp), viscosity (μ), and thermal conductivity (k) were also required to correctly define and simulate the fluid.Table 3
Gas constants of air, the fluid mixture, and its components.
Gas
MW(g/mol)
R  (J/(kg*K))
Air
28.96
287.1
C
O
2
44.01
188.9
Argon
39.95
208.1
Mixture
41.94
198.7The componentCp data over the experiment temperature range was obtained using the NASA Thermo Build Program [32]. The mixture Cp properties were calculated by mass averaging the component values. From Cp and R, the mixture γ properties were obtained with the well-known relation γ=Cp/(Cp-R). The viscosity of a gas can be defined by Sutherland’s formula. The constants for this equation for common gases, such as CO2, are widely available. These constants are not established for Argon. Instead, tabulated μ data were obtained for the pressure and temperature range [33]. The mixture viscosity was calculated using the Herning and Zipperer approximation [34], as defined in (2):
(2)μmix=∑‍(XiμiMWi)∑‍(XiMWi),
where Xi, μi, and MWi are the mole fraction, viscosity, and molecular weight of the ith component gas. The CO2 and Argon thermal conductivity data was obtained from the FINE/Turbo gas properties database. The Wassiljewa equation with the Mason and Saxena modification was used to calculate the mixture thermal conductivity [34]. Details of this procedure and tables of the component and mixture data are described by Knapke [28]. Comparisons between air, the mixture, and component Cp, γ, μ, and k are shown in Figure 8. The mixture γ and μ are similar to that of air, though the mixture γ does have a ±2.5% difference from air at the extremes of the experiment temperature range (250 K and 450 K). The Cp and k curves for the mixture and air are different across the entire temperature range.Gas properties for air, CO2, Argon and the Run 13 mixture.
(a)
C
p
(b)
γ
(c)
μ
(d)
kAlong with the common gases available in the FINE/Turbo fluid properties database, user defined fluids can be added. The CO2 and Argon mixture was created using the real gas option, which requires the definition of R, Cp, μ, and k as functions of temperature. Coefficients derived from a Sutherland’s formula curve fit of μ and quadratic curve fits of Cp and k were used to define the temperature-dependent properties. Although the data is input with the “real gas” option, it is actually an imperfect gas as the ideal gas law is used.


## 2.1. Computational Grid

The main passage grids were generated with the turbomachinery meshing tool AutoGrid. Constant axial inlet and outlet planes were positioned at the approximate locations of the upstream and downstream measurement planes. In addition, meridional control lines were placed along the aspiration slot to axially cluster the grids at that location. Figure2 depicts the meridional view of the computational domain. Previous experience with simulating the CRAC geometry by Knapke et al. [25] and grid dependency studies by Guidotti and Turner [26, 27] and Knapke [28] have shown that the results are sensitive to the spanwise grid resolution. For this reason, a spanwise cell count of 172 was chosen. An additional 12 spanwise cells were placed in the 10 mil rotor tip clearance gaps. A first cell wall spacing of 3×10-6 meters was applied to all blade rows which resulted in all y+ values below 5.Figure 2
AutoGrid meridional view (blades extruded beyond endwalls to ensure robust grid generation).An O4H block topology was used for each blade passage, which consists of an O-grid block around the blade and H-grid blocks at the inlet, outlet, upper, and lower sides of the passage. High stagger optimization and blunt trailing edge meshing features were used for the rotor passages. The trailing edge detail of both rotors is shown in Figure3. The blunt trailing edge was a design feature of the CRAC rotors. With the blunt trailing edge, the skin mesh is created with a C-grid. The tip clearance mesh for both rotors utilized an H-grid block surrounded by a C-grid block. A nonmatching boundary (NMB) is automatically applied to account for any mismatch in the number of cells between the suction and pressure sides. For NLH simulations, Numeca recommends the number of cells in the tangential direction be defined by the following equation:
(1)Ncellθ=30Nharmnadjn,
where Nharm is the number of harmonics and nadj and n are the adjacent and current blade counts, respectively. This guideline is established so the adjacent blade row has enough grid points (30) to define each harmonic mode. Following this guideline for Nharm=3, the tangential number of cells in the IGV, rotor one and rotor two blade rows were set to 72, 116 and 104, respectively. Small axial cell expansion ratios were maintained at the inlet and exit of each blade row to negate axial wave dissipation. The total cell counts are listed in Table 2.Table 2
Number of cells per grid.
Millions of cells
Total (aspirated)
17.53
Total (nonaspirated)
17.19
IGV
2.06
Rotor one
6.29
Rotor two
8.84
Aspiration slot and plenum
0.34Midspan rotor grids at the trailing edge.
(a)
Rotor one
(b)
Rotor twoIGG, a general structured multiblock 3D meshing tool, was used to generate the aspiration flow path grids. The aspiration flow path domains were held as close to the test geometry as possible, while making small modifications to simplify the grid generation. Special care was taken to maintain the areas of the slot and plenum. The area directly impacts the way in which the flow path chokes. The number of cells and cell spacings along the span of the slot were defined such that the spanwise grid densities were similar to the main passage grids. The width of the slot was resolved with 16 cells. Wall spacings of3×10-6 meters were applied throughout the aspiration flow path to maintain y+ values below 5. Decisions on the number of cells and clustering were driven by grid quality metrics, such as orthogonality, aspect ratio, and expansion ratio. Figure 4 shows the aspiration flow path grid and the clustering along the rotor two surface.Main passage clustering and aspiration flow path grid.
(a)
Rotor two blade surface grid
(b)
Aspiration slot and plenum grid


## 2.2. Boundary Conditions

The CRAC testing at MIT consisted of 14 test cases. Of these cases, Run 13 was chosen for two reasons. First, Run 13 was near 100% design speed for both rotors. Second, sufficient inlet and outlet data was available to provide accurate boundary conditions and ample data for comparison. From the range of acceptable data of the blowdown test, the 250 milliseconds data point was chosen. The details of Run 13 and the data used are presented by Parker [6].Total conditions with axial flow was applied at the inlet boundary. Inlet rake measurements for Run 13 were not available, however Run 10 data was used as an appropriate substitute given the similar mass flow and rotor speeds. A pressure screen was located upstream of the test section to improve the rig operability and provide uniform flow to the IGV (Figure1(a)). This screen mixes out the viscous boundary layer. Given the short distance between the pressure screen and the inlet rake and the low radial resolution of that rakes, the viscous boundary layer was not measured accurately. Therefore, TEXSTAN was used to calculate the viscous boundary layer that developed between the pressure screen and the computational inlet plane. TEXSTAN is a two-dimensional finite difference code for solving the boundary layer mass, momentum, and heat transfer equations [29]. Two annular duct TEXSTAN simulations were completed. One to simulate the total temperature from the supply tank to the test section and another to capture the viscous boundary layer between the screen and the test section. Both simulations assume constant isothermal walls at room temperature (293 K) and a uniform inlet flow boundary condition. A mixing-length model for the turbulent viscosity was used.The calculated total pressure profile is shown in Figure5(a). To check the accuracy of the TEXSTAN boundary layers, the total temperature profile is compared with the Run 10 test data in Figure 5(b). The TEXSTAN results agree well with the test data. Some difference occur near the hub, but that difference is likely the result of the heat produced by rotor one motor, which was not modeled in the TEXSTAN simulation. This is strong evidence that the blowdown rig is not adiabatic. The TEXSTAN total pressure and the test total temperature profiles were applied at the inlet boundary. An unintended consequence of the aforementioned pressure screen was that it produced inlet distortion, which resulted in a 4% variation of total pressure around the annulus, as described by Parker [6]. This distortion was not accounted for in the CRAC simulations because the distortion was not measured along the entire span.TEXSTAN inlet profiles and the Parker [6] inlet total temperature profile.
(a)
Total pressure profile
(b)
Total temperature profilesThe aspiration and main passage grids were combined in IGG. The connection between the slot and rotor two passage grids was made with a Full Nonmatching Boundary (FNMB). This FNMB was required because the slot grid boundary points do not correspond exactly to rotor two blade boundary points. Although a direct connection boundary would provide a more accurate transfer of flow details, the increase in complexity of a connecting topology greatly outweighed the benefits. Moreover, this topology would require modification of AutoGrid blocks in IGG. This type of modification is not recommended by Numeca as it removes the grid smoothing performed during the AutoGrid grid generation process.The FNMB uses interpolation to transfer information between the nonmatching boundaries. Although the grid cells are not coincident along the FNMB boundary, significant effort was made to have similar spacings on either side of the boundary. A poorly defined connection can result from cell spacings that vary greatly across the FNMB. Figure6 shows the results of a triangulation and projection procedure performed in IGG. This procedure visualizes the same connection process used by the solver.Main passage grid on the rotor two suction surface.
(a)
Full blade to slot FNMB
(b)
FNMB near hubAt the exit of the main passage, a simple radial equilibrium static pressure boundary condition was used. The design mass flow of 1% of the inlet mass flow was applied to aspiration plenum exit boundary and the aspiration flow path was allowed to choke, thus setting the mass flow. Attempts were made to switch to a static pressure boundary condition because it would more accurately represent the plenum. However, a low pressure was required to choke the flow, causing negative pressures and densities when switching from steady to NLH computations. Since the exit of the plenum is reasonably far from the location where the flow chokes, the mass flow exit boundary condition was deemed acceptable.The rotating speeds of the two rotors at the 250 millisecond data point were 12600 rpm and 10020 rpm. Simulations at these speeds resulted in an upstream shock position that is not consistent with the design intent. It was found that raising the rotating speeds to 103% of the experimental value moved the shock to the correct position, as depicted in Figure7. This increase of the rotor speed was deemed plausible for two reasons. The first reason is the potential uncertainties of the rotor twist during operation. Although the amount of twist was accounted for in the design, the actual rotor tip positions during the test were not measured. Second, the corrected speed, which sets the shock location, is dependent on temperature and R. During the experiment, the rotor could have been at a lower temperature than expected because of heat transfer. This difference would make the corrected speed higher. Additionally, variations in the gas mixture could affect R and further change the corrected speed.Static pressure contours at midspan of rotor two. Adiabatic three harmonic solutions.
(a)
100% rotor speed
(b)
103% rotor speedRecent research concerning the heat transfer boundary conditions at solid walls for compressor simulations, by Bruna and Turner [30] and Isomura et al. [31], has shown the adiabatic wall approximation to be inaccurate. The short duration of the CRAC blowdown test and the difference between the fluid and metal temperatures prior to the test results in more heat transfer than seen during more traditional test methods [6]. For these reasons both adiabatic and isothermal solid wall boundary conditions were simulated for the aspirated case. The surface temperatures were not measured during the experiment. Therefore, a conservative estimate of room temperature (293 K) was assumed for all surfaces for the isothermal simulation.At the periodic boundary, the time mean flow is treated as a traditional periodic boundary. For the harmonic perturbation equations, the periodic boundary condition comprises a phase-shift of the harmonic amplitudes related to the pitch of the given blade row and the perturbation periodicity. At the interface boundary, Fourier Transforms (FT) in both time and space are performed on both sides of the interface. The harmonics produced by the temporal decomposition are matched across the interface to the discrete number of perturbation frequencies. From this matching process and the spatial FT, the time mean and perturbation flows can be reconstructed. Continuity across the interface boundary cannot be exactly reached for a finite number of harmonics; however, increasing the number of harmonics reduces the discontinuity error.


## 2.3. Initialization

To improve convergence, steady mixing-plane solutions were used to initialize the NLH simulations. Steady computations were completed using full multigrid. The steady cases were initialized on a coarse grid from constant values and a lower back pressure. The back pressure and grid density were increased until the final grid and the Run 13 back pressure values were reached. This final steady solution was used to initialize an NLH simulation with one harmonic per blade row. For the adiabatic cases, the number of harmonics was incremented to three in successive simulations. The isothermal cases were converged with two harmonics.


## 2.4. Gas Model

The CO2 and Argon mixture used as the working fluid has a 44% higher molecular weight than air, which resulted in lower physical rotor speeds for a given corrected speed. The physical speed reduction in turn reduced the stresses in the rotating components of the compressor rig. For the Run 13 mixture, the molar percentages of CO2 and Argon were 48.91% and 51.09%. Over the temperature ranges of the experiment (250 K to 450 K), the mixture γ is similar to the γ of air. However, the higher molecular weight results in a lower gas constant (R), as shown in Table 3. Calculation of the mixture specific heat at constant pressure (Cp), viscosity (μ), and thermal conductivity (k) were also required to correctly define and simulate the fluid.Table 3
Gas constants of air, the fluid mixture, and its components.
Gas
MW(g/mol)
R  (J/(kg*K))
Air
28.96
287.1
C
O
2
44.01
188.9
Argon
39.95
208.1
Mixture
41.94
198.7The componentCp data over the experiment temperature range was obtained using the NASA Thermo Build Program [32]. The mixture Cp properties were calculated by mass averaging the component values. From Cp and R, the mixture γ properties were obtained with the well-known relation γ=Cp/(Cp-R). The viscosity of a gas can be defined by Sutherland’s formula. The constants for this equation for common gases, such as CO2, are widely available. These constants are not established for Argon. Instead, tabulated μ data were obtained for the pressure and temperature range [33]. The mixture viscosity was calculated using the Herning and Zipperer approximation [34], as defined in (2):
(2)μmix=∑‍(XiμiMWi)∑‍(XiMWi),
where Xi, μi, and MWi are the mole fraction, viscosity, and molecular weight of the ith component gas. The CO2 and Argon thermal conductivity data was obtained from the FINE/Turbo gas properties database. The Wassiljewa equation with the Mason and Saxena modification was used to calculate the mixture thermal conductivity [34]. Details of this procedure and tables of the component and mixture data are described by Knapke [28]. Comparisons between air, the mixture, and component Cp, γ, μ, and k are shown in Figure 8. The mixture γ and μ are similar to that of air, though the mixture γ does have a ±2.5% difference from air at the extremes of the experiment temperature range (250 K and 450 K). The Cp and k curves for the mixture and air are different across the entire temperature range.Gas properties for air, CO2, Argon and the Run 13 mixture.
(a)
C
p
(b)
γ
(c)
μ
(d)
kAlong with the common gases available in the FINE/Turbo fluid properties database, user defined fluids can be added. The CO2 and Argon mixture was created using the real gas option, which requires the definition of R, Cp, μ, and k as functions of temperature. Coefficients derived from a Sutherland’s formula curve fit of μ and quadratic curve fits of Cp and k were used to define the temperature-dependent properties. Although the data is input with the “real gas” option, it is actually an imperfect gas as the ideal gas law is used.


## 3. Results

### 3.1. Convergence History

The FINE/Turbo solver outputs iteration histories of residuals and inlet and exit mass flow. Conservation of mass from the inlet to exit was used to measure convergence. The residuals were used as a secondary measure of convergence as the residuals for each simulation are relative to the previous simulation. An example of the convergence history from initialization to NLH convergence is presented in Figure9. The Pb values are the main passage exit pressures, with 155 kPa being the Run 13 value. The fine grid level corresponds to the grid presented in Table 2. The oscillations of mass flow for the NLH simulations were less than 0.5% of the inlet mass flow, which was considered acceptable for convergence. The mass flow discontinuity across the NLH interface boundaries was within the mass flow variations during convergence.Figure 9
Example convergence history (adiabatic 100% rotor speed case).From Figure9 it can be seen how the steady mixing-plane solution is used to initialize the NLH solution. These steady solutions have not been shown in this paper because they are less accurate than the unsteady results for two reasons. The steady simulation relies on the mixing-plane interface boundary, which mixes out the unsteady features. In addition, the reflective nature of the steady interface boundary produced spurious shock reflections.


### 3.2. Comparison with Data

Stage performance for the simulations and the experiment are shown in Table4. The simulated efficiencies were calculated using an imperfect gas form of the efficiency calculation, as described by Knapke [28]. In addition, the stage performance quantities were calculated using the same area averaged approach (subscript AA) as used to calculate the experimental values, providing a direct comparison to the experiment.Table 4
Experiment and simulation pressure ratio, temperature ratio, adiabatic efficiency, and mass flow. MIT experiment mass flow comes from an average across the three windows at the inlet.
Aspiration
N
harm
Thermal BC
PR
TR
η (%)
P
R
AA
T
R
AA
η
AA (%)
m
˙
inlet kg/s
No
3
Adiabatic
2.889
1.438
80.0
2.884
1.437
80.2
23.86
Yes
2.887
1.417
84.2
2.895
1.417
84.6
24.63
Yes
2
Isothermal
2.887
1.390
90.2
2.899
1.395
89.6
24.88
MIT experiment
—
—
—
2.934
1.411
87.9
24.40The mass flow of the aspirated cases are within 2% of the experiment mass flow. The small increase in mass flow for the aspirated simulations, when compared to the experiment, corresponds to the low pressure ratio and 103% rotor speeds. The low simulated pressure ratios could be the result of not modeling the inlet distortion. The area averaged efficiency is lower than the experiment, whereas the isothermal case efficiency is higher. However, the isothermal case under predicts the total temperature ratio, suggesting that the simulated wall temperatures cause more heat transfer than occurred during the experiment. Some of the simulation to experiment discrepancy could be attributed to the simulation representing one point during the blowdown test, while the experiment values are an average over a small time range.Because the NLH perturbations are periodic, the time mean solution is also the time-averaged solution. The time mean variables were mass averaged in the circumferential direction and spanwise profiles were extracted just upstream of the exit plane. The exit profiles for the aspirated solutions are compared to the experiment and an APNASA simulation in Figure10. During the design process, the MIT APNASA simulation was completed as a prediction of the CRAC performance. The FINE/Turbo simulation boundary conditions are consistent with the experiment at 250 milliseconds.Spanwise exit profiles of efficiency, total temperature ratio, and total pressure ratio for comparison between the FINE/Turbo simulations, MIT simulations and the experiment.
(a)
(b)
(c)Both the adiabatic case and the predicted performance of the APNASA simulation fail to capture the increased efficiency at the tip, which the isothermal case captures. The adiabatic total temperature profile matches the experiment well between 30% and 80% span, but fails to capture the decrease near the endwalls. The isothermal case better matches the decreasing trends at the endwalls. However, between the endwalls the isothermal case under predicts the total temperature. This suggests that the isothermal case heat transfer is too high in those regions. Good agreement with the pressure ratio is shown for both the adiabatic and isothermal cases along most of the span.


### 3.3. Aspiration Mass Flow

Details of the time averaged mass flow through the aspiration plenum are listed in Table5. Fluid was not found to move from the plenum back into the blade passage in a time averaged or unsteady sense. Therefore, m˙plenum is equivalent to the mass flow passing through the slot. Because the mass flow was not measured during the experiment whose results provided the main passage boundary conditions, a direct comparison with data is not possible. Additional testing of the aspiration flow path suggested that the flow path choked at closer to 0.5% of the inlet mass flow. Both the adiabatic and isothermal cases have near 0.5% aspiration, which provides more evidence that the simulations accurately represent the experiment. The isothermal case mass flow is lower than the adiabatic case.Table 5
Plenum exit mass flow and aspiration percentage of inlet mass flow.
N
harm
Thermal BC
m
˙
inlet kg/s
m
˙
plenum kg/s
% ofm˙inlet
3
Adiabatic
24.63
0.1521
0.618
2
Isothermal
24.88
0.1376
0.553


### 3.4. Choking Mechanism

To illustrate the way in which the aspiration flow path chokes, cut-planes displayingMrel contours at different spanwise locations and an isosurface of Mrel are presented in Figure 11. The isosurface extends to 50% of the slot span, which indicates that the top half of the slot is not choked. The flow entering the plenum from the slot forms a jet that creates a fluid barrier between two regions, one upstream and one downstream. To illustrate the fluid barrier, streamlines were released near the tip of the plenum and are seen traveling downward (Figure 11(a)). The flow rotates on either side of the slot, but does not pass across the jet of fluid. The area of the downstream side constricts, causing the flow to choke at a quarter of the slot span. From that point to the base of the slot, the area opens and the flow accelerates to become supersonic. On the upstream side, the Mrel=1 isosurface begins at the base of the slot, indicating the place where the upstream region chokes (Figure 11(b)). The two regions join below the slot and the flow accelerates to the exit because of the increasing plenum area.Plenum and slot time meanMrel contours and isosurface for the adiabatic case.
(a)
M
rel contours with streamlines
(b)
M
rel
=
1 isosurfaceAlthough the slot was thought to choke and restrict the aspiration mass flow, the limiting area of the aspiration flow path is actually in the plenum. The jet entering from the slot creates two regions that choke at different points in the plenum. Therefore, the throat areas of the upstream and downstream regions limit the mass flow. The isothermal contours and isosurface are similar to those of the adiabatic case and are not shown for brevity.


### 3.5. Slot Profiles

To understand the interaction between the main passage and the aspiration flow path, it is important to investigate the flow entering the slot. Spanwise distributions of several time-averaged and mass-averaged quantities for the adiabatic and isothermal cases are shown in Figure12. Below 20% span the entropy increase for all cases has been attributed to the combination of two main passage flow features. The first feature is a rotor one secondary flow, as shown by the streamwise vorticity contours and streamlines in Figure 13, which passes through the rotor two passage. Second, the rotor two shock is upstream of the slot near the hub, as depicted in Figure 14. These main passage flow path features also contribute to the low relative total pressure below 20% span. The isothermal relative total temperatures are lower than the adiabatic case due to the heat transfer on the solid surfaces. This lower temperature causes the lower entropy distribution.Time-averaged spanwise distributions of several quantities at the slot inlet. Mass averaged across the slot width at each spanwise location.
(a)
m
˙
per
A
(b)
Mass averaged entropy
(c)
Relative total temperature
(d)
Relative total pressureFigure 13
Progression of a secondary flow along the rotor one suction surface. Stream ribbons and cut plans colored by streamwise vorticity. Time mean solution. Red is positive, blue is negative.Figure 14
Rotor two axisymmetric blade suction surface static pressure contours for the adiabatic aspirated case.Them˙perA distribution can be better understood by considering the equation for mass flow of an ideal compressible gas, as shown below:
(3)m˙=APtTtγRM(1+γ-12M2)-(γ+1)/2(γ-1).
In the case of the aspiration slot, the relative quantities are relevant. Along the slot span, the width is essentially constant and does not affect the mass flow. The low mass flow between 0% and 20% of the slot span is a result of low relative total pressure, high relative total temperature, and a sonic Mrel. The low Mrel at higher span, seen in Figure 11(a), causes the low mass flow in that region. The mass flow peak at 40% span is produced by an increased relative total pressure. Despite the overall decrease in temperature for the isothermal case, the mass flow is lower than the adiabatic case because of a decreased pressure.


### 3.6. Unsteady Results

Unsteady effects on the aspiration flow path were investigated by interrogating the NLH reconstruction of the flow field. A time history of the mass flow as a percentage of inlet flow is shown in Figure15 for the adiabatic and isothermal cases. As expected from the time-averaged mass flows in Table 5, the mass flows history of the adiabatic case is higher than the isothermal case. A peak to valley mass flow variation of 0.05% occurs at 9×10-5 seconds. This drop in mass flow is attributed to the rotor one wake. Figure 16 presents entropy contours along the blade and through the passage. The wake in the main passage, indicated by the high entropy region on the cut-plane, passes the slot at time 5.73×10-5 seconds (Figure 16(a)). However, the effect of the wake traveling through the boundary layer (the high entropy region on the blade near the slot in Figure 16(b)) does not reach the slot until time 9.54×10-5 seconds. The additional drops in mass flow are the result of the Fourier representation of the solution.Figure 15
Time history of the slot mass flow percentage relative to the inlet mass flow.The rotor one wake depicted by instantaneous entropy contours on the rotor two blade and a 40% span cut-plane (equivalent to 50% of the slot span). Red is high entropy, blue is low entropy.
(a)
t
=
5
.
73
×
1
0
-
5 seconds
(b)
t
=
9
.
54
×
1
0
-
5 secondsMass-averaged spanwise distributions were also calculated for the NLH reconstructed time steps. Contours of the unsteady spanwise distributions ofm˙perA and entropy for the adiabatic and isothermal case are depicted in Figure 17. The m˙perA is normalized by the overall time mean value. The adiabatic contours exhibit more unsteady variations than the isothermal case because the simulation used three harmonics to represent the unsteadiness, instead of two. The entropy field of the isothermal case is lower as a result of heat transfer to the solid walls.Unsteady spanwise distributions of normalizedm˙perA and entropy.
(a)
Adiabatic, normalizedm˙perA
(b)
Isothermal, normalizedm˙perA
(c)
Adiabatic, entropy(J/(kg*K))
(d)
Isothermal, entropy(J/(kg*K))The most striking feature of the mass flow contours is the sharp drop that varies from8×10-5 seconds near the hub to 1×10-4 seconds near the tip of the slot. As with the time history, this drop in mass flow corresponds to the passing rotor one wake. Throughout most of the contour, an increase in entropy corresponds to a decrease in mass flow. Unexpectedly, in the wake region below 40% span, the passing wake results in an area of decreased entropy. In a time mean sense, the forward shock position contributes to an increase in entropy below 40% span. However, the unsteady contours of static pressure near the hub, in Figure 18, show that the passing rotor one wake moves the shock downstream to a position just upstream of the slot. This movement downstream causes the decreased entropy, which is opposite the usual effect of a passing wake on the flow field.Instantaneous static pressure contours on a 10% rotor two span cut-plane.
(a)
t
=
0 seconds
(b)
t
=
7
.
64
×
1
0
-
5 secondsTime mean and unsteady plots of the slot inlet are shown in Figure19. Increased mass flow occurs along the downstream side of the slot over the entire span. The passage of the rotor one wake is visible as a reduction of m˙perA moving from the hub at time 0.6Tp to the tip at time 0.9Tp. The wake is similarly visible in the entropy plots. The entropy decreases at lower span between 0.4Tp and 0.6Tp as the wake moves the shock and the entropy increases at higher span between 0.6Tp and 0.8Tp. In general, the entropy is higher where m˙perA is lower.Time mean and unsteady contours of the slot inlet for the adiabatic case. The width to height aspect ratio was adjusted to better view the variation of properties along the width of the slot.Tp=1.273×10-4 seconds, which is the rotor one blade passing period. Left is upstream, right is downstream.
(a)
m
˙
per
A, (kg/s*m2)
(b)
Entropy,(J/kg*K)
(c)
Relative total pressure,(Pa)
(d)
Relative total temperature,(K)To understand the correlation between high entropy and lowm˙perA, a time mean contour of entropy at 50% slot span is displayed in Figure 20. At the upstream side of the slot the flow deepest in the boundary layer is ingested. Moving downstream, flow from progressively higher in the boundary layer is ingested. This higher flow has lower entropy and higher total pressure. As indicated by (3), a higher total pressure on the downstream side results in a higher mass flow.Figure 20
Time mean entropy contours with relative velocity vectors at 50% slot span.Along the choked region of the slot (below 50% slot span), the total pressure contours are similar to them˙perA. Because Mrel decreases at higher span, the m˙perA decreases even though the total pressure increases (again related to (3)). The total temperature across the width of the slot is nearly constant for the time mean flow. In an unsteady sense, the passing rotor one wake causes variations of total temperature from time 0.5Tp to 0.9Tp.


## 3.1. Convergence History

The FINE/Turbo solver outputs iteration histories of residuals and inlet and exit mass flow. Conservation of mass from the inlet to exit was used to measure convergence. The residuals were used as a secondary measure of convergence as the residuals for each simulation are relative to the previous simulation. An example of the convergence history from initialization to NLH convergence is presented in Figure9. The Pb values are the main passage exit pressures, with 155 kPa being the Run 13 value. The fine grid level corresponds to the grid presented in Table 2. The oscillations of mass flow for the NLH simulations were less than 0.5% of the inlet mass flow, which was considered acceptable for convergence. The mass flow discontinuity across the NLH interface boundaries was within the mass flow variations during convergence.Figure 9
Example convergence history (adiabatic 100% rotor speed case).From Figure9 it can be seen how the steady mixing-plane solution is used to initialize the NLH solution. These steady solutions have not been shown in this paper because they are less accurate than the unsteady results for two reasons. The steady simulation relies on the mixing-plane interface boundary, which mixes out the unsteady features. In addition, the reflective nature of the steady interface boundary produced spurious shock reflections.


## 3.2. Comparison with Data

Stage performance for the simulations and the experiment are shown in Table4. The simulated efficiencies were calculated using an imperfect gas form of the efficiency calculation, as described by Knapke [28]. In addition, the stage performance quantities were calculated using the same area averaged approach (subscript AA) as used to calculate the experimental values, providing a direct comparison to the experiment.Table 4
Experiment and simulation pressure ratio, temperature ratio, adiabatic efficiency, and mass flow. MIT experiment mass flow comes from an average across the three windows at the inlet.
Aspiration
N
harm
Thermal BC
PR
TR
η (%)
P
R
AA
T
R
AA
η
AA (%)
m
˙
inlet kg/s
No
3
Adiabatic
2.889
1.438
80.0
2.884
1.437
80.2
23.86
Yes
2.887
1.417
84.2
2.895
1.417
84.6
24.63
Yes
2
Isothermal
2.887
1.390
90.2
2.899
1.395
89.6
24.88
MIT experiment
—
—
—
2.934
1.411
87.9
24.40The mass flow of the aspirated cases are within 2% of the experiment mass flow. The small increase in mass flow for the aspirated simulations, when compared to the experiment, corresponds to the low pressure ratio and 103% rotor speeds. The low simulated pressure ratios could be the result of not modeling the inlet distortion. The area averaged efficiency is lower than the experiment, whereas the isothermal case efficiency is higher. However, the isothermal case under predicts the total temperature ratio, suggesting that the simulated wall temperatures cause more heat transfer than occurred during the experiment. Some of the simulation to experiment discrepancy could be attributed to the simulation representing one point during the blowdown test, while the experiment values are an average over a small time range.Because the NLH perturbations are periodic, the time mean solution is also the time-averaged solution. The time mean variables were mass averaged in the circumferential direction and spanwise profiles were extracted just upstream of the exit plane. The exit profiles for the aspirated solutions are compared to the experiment and an APNASA simulation in Figure10. During the design process, the MIT APNASA simulation was completed as a prediction of the CRAC performance. The FINE/Turbo simulation boundary conditions are consistent with the experiment at 250 milliseconds.Spanwise exit profiles of efficiency, total temperature ratio, and total pressure ratio for comparison between the FINE/Turbo simulations, MIT simulations and the experiment.
(a)
(b)
(c)Both the adiabatic case and the predicted performance of the APNASA simulation fail to capture the increased efficiency at the tip, which the isothermal case captures. The adiabatic total temperature profile matches the experiment well between 30% and 80% span, but fails to capture the decrease near the endwalls. The isothermal case better matches the decreasing trends at the endwalls. However, between the endwalls the isothermal case under predicts the total temperature. This suggests that the isothermal case heat transfer is too high in those regions. Good agreement with the pressure ratio is shown for both the adiabatic and isothermal cases along most of the span.


## 3.3. Aspiration Mass Flow

Details of the time averaged mass flow through the aspiration plenum are listed in Table5. Fluid was not found to move from the plenum back into the blade passage in a time averaged or unsteady sense. Therefore, m˙plenum is equivalent to the mass flow passing through the slot. Because the mass flow was not measured during the experiment whose results provided the main passage boundary conditions, a direct comparison with data is not possible. Additional testing of the aspiration flow path suggested that the flow path choked at closer to 0.5% of the inlet mass flow. Both the adiabatic and isothermal cases have near 0.5% aspiration, which provides more evidence that the simulations accurately represent the experiment. The isothermal case mass flow is lower than the adiabatic case.Table 5
Plenum exit mass flow and aspiration percentage of inlet mass flow.
N
harm
Thermal BC
m
˙
inlet kg/s
m
˙
plenum kg/s
% ofm˙inlet
3
Adiabatic
24.63
0.1521
0.618
2
Isothermal
24.88
0.1376
0.553


## 3.4. Choking Mechanism

To illustrate the way in which the aspiration flow path chokes, cut-planes displayingMrel contours at different spanwise locations and an isosurface of Mrel are presented in Figure 11. The isosurface extends to 50% of the slot span, which indicates that the top half of the slot is not choked. The flow entering the plenum from the slot forms a jet that creates a fluid barrier between two regions, one upstream and one downstream. To illustrate the fluid barrier, streamlines were released near the tip of the plenum and are seen traveling downward (Figure 11(a)). The flow rotates on either side of the slot, but does not pass across the jet of fluid. The area of the downstream side constricts, causing the flow to choke at a quarter of the slot span. From that point to the base of the slot, the area opens and the flow accelerates to become supersonic. On the upstream side, the Mrel=1 isosurface begins at the base of the slot, indicating the place where the upstream region chokes (Figure 11(b)). The two regions join below the slot and the flow accelerates to the exit because of the increasing plenum area.Plenum and slot time meanMrel contours and isosurface for the adiabatic case.
(a)
M
rel contours with streamlines
(b)
M
rel
=
1 isosurfaceAlthough the slot was thought to choke and restrict the aspiration mass flow, the limiting area of the aspiration flow path is actually in the plenum. The jet entering from the slot creates two regions that choke at different points in the plenum. Therefore, the throat areas of the upstream and downstream regions limit the mass flow. The isothermal contours and isosurface are similar to those of the adiabatic case and are not shown for brevity.


## 3.5. Slot Profiles

To understand the interaction between the main passage and the aspiration flow path, it is important to investigate the flow entering the slot. Spanwise distributions of several time-averaged and mass-averaged quantities for the adiabatic and isothermal cases are shown in Figure12. Below 20% span the entropy increase for all cases has been attributed to the combination of two main passage flow features. The first feature is a rotor one secondary flow, as shown by the streamwise vorticity contours and streamlines in Figure 13, which passes through the rotor two passage. Second, the rotor two shock is upstream of the slot near the hub, as depicted in Figure 14. These main passage flow path features also contribute to the low relative total pressure below 20% span. The isothermal relative total temperatures are lower than the adiabatic case due to the heat transfer on the solid surfaces. This lower temperature causes the lower entropy distribution.Time-averaged spanwise distributions of several quantities at the slot inlet. Mass averaged across the slot width at each spanwise location.
(a)
m
˙
per
A
(b)
Mass averaged entropy
(c)
Relative total temperature
(d)
Relative total pressureFigure 13
Progression of a secondary flow along the rotor one suction surface. Stream ribbons and cut plans colored by streamwise vorticity. Time mean solution. Red is positive, blue is negative.Figure 14
Rotor two axisymmetric blade suction surface static pressure contours for the adiabatic aspirated case.Them˙perA distribution can be better understood by considering the equation for mass flow of an ideal compressible gas, as shown below:
(3)m˙=APtTtγRM(1+γ-12M2)-(γ+1)/2(γ-1).
In the case of the aspiration slot, the relative quantities are relevant. Along the slot span, the width is essentially constant and does not affect the mass flow. The low mass flow between 0% and 20% of the slot span is a result of low relative total pressure, high relative total temperature, and a sonic Mrel. The low Mrel at higher span, seen in Figure 11(a), causes the low mass flow in that region. The mass flow peak at 40% span is produced by an increased relative total pressure. Despite the overall decrease in temperature for the isothermal case, the mass flow is lower than the adiabatic case because of a decreased pressure.


## 3.6. Unsteady Results

Unsteady effects on the aspiration flow path were investigated by interrogating the NLH reconstruction of the flow field. A time history of the mass flow as a percentage of inlet flow is shown in Figure15 for the adiabatic and isothermal cases. As expected from the time-averaged mass flows in Table 5, the mass flows history of the adiabatic case is higher than the isothermal case. A peak to valley mass flow variation of 0.05% occurs at 9×10-5 seconds. This drop in mass flow is attributed to the rotor one wake. Figure 16 presents entropy contours along the blade and through the passage. The wake in the main passage, indicated by the high entropy region on the cut-plane, passes the slot at time 5.73×10-5 seconds (Figure 16(a)). However, the effect of the wake traveling through the boundary layer (the high entropy region on the blade near the slot in Figure 16(b)) does not reach the slot until time 9.54×10-5 seconds. The additional drops in mass flow are the result of the Fourier representation of the solution.Figure 15
Time history of the slot mass flow percentage relative to the inlet mass flow.The rotor one wake depicted by instantaneous entropy contours on the rotor two blade and a 40% span cut-plane (equivalent to 50% of the slot span). Red is high entropy, blue is low entropy.
(a)
t
=
5
.
73
×
1
0
-
5 seconds
(b)
t
=
9
.
54
×
1
0
-
5 secondsMass-averaged spanwise distributions were also calculated for the NLH reconstructed time steps. Contours of the unsteady spanwise distributions ofm˙perA and entropy for the adiabatic and isothermal case are depicted in Figure 17. The m˙perA is normalized by the overall time mean value. The adiabatic contours exhibit more unsteady variations than the isothermal case because the simulation used three harmonics to represent the unsteadiness, instead of two. The entropy field of the isothermal case is lower as a result of heat transfer to the solid walls.Unsteady spanwise distributions of normalizedm˙perA and entropy.
(a)
Adiabatic, normalizedm˙perA
(b)
Isothermal, normalizedm˙perA
(c)
Adiabatic, entropy(J/(kg*K))
(d)
Isothermal, entropy(J/(kg*K))The most striking feature of the mass flow contours is the sharp drop that varies from8×10-5 seconds near the hub to 1×10-4 seconds near the tip of the slot. As with the time history, this drop in mass flow corresponds to the passing rotor one wake. Throughout most of the contour, an increase in entropy corresponds to a decrease in mass flow. Unexpectedly, in the wake region below 40% span, the passing wake results in an area of decreased entropy. In a time mean sense, the forward shock position contributes to an increase in entropy below 40% span. However, the unsteady contours of static pressure near the hub, in Figure 18, show that the passing rotor one wake moves the shock downstream to a position just upstream of the slot. This movement downstream causes the decreased entropy, which is opposite the usual effect of a passing wake on the flow field.Instantaneous static pressure contours on a 10% rotor two span cut-plane.
(a)
t
=
0 seconds
(b)
t
=
7
.
64
×
1
0
-
5 secondsTime mean and unsteady plots of the slot inlet are shown in Figure19. Increased mass flow occurs along the downstream side of the slot over the entire span. The passage of the rotor one wake is visible as a reduction of m˙perA moving from the hub at time 0.6Tp to the tip at time 0.9Tp. The wake is similarly visible in the entropy plots. The entropy decreases at lower span between 0.4Tp and 0.6Tp as the wake moves the shock and the entropy increases at higher span between 0.6Tp and 0.8Tp. In general, the entropy is higher where m˙perA is lower.Time mean and unsteady contours of the slot inlet for the adiabatic case. The width to height aspect ratio was adjusted to better view the variation of properties along the width of the slot.Tp=1.273×10-4 seconds, which is the rotor one blade passing period. Left is upstream, right is downstream.
(a)
m
˙
per
A, (kg/s*m2)
(b)
Entropy,(J/kg*K)
(c)
Relative total pressure,(Pa)
(d)
Relative total temperature,(K)To understand the correlation between high entropy and lowm˙perA, a time mean contour of entropy at 50% slot span is displayed in Figure 20. At the upstream side of the slot the flow deepest in the boundary layer is ingested. Moving downstream, flow from progressively higher in the boundary layer is ingested. This higher flow has lower entropy and higher total pressure. As indicated by (3), a higher total pressure on the downstream side results in a higher mass flow.Figure 20
Time mean entropy contours with relative velocity vectors at 50% slot span.Along the choked region of the slot (below 50% slot span), the total pressure contours are similar to them˙perA. Because Mrel decreases at higher span, the m˙perA decreases even though the total pressure increases (again related to (3)). The total temperature across the width of the slot is nearly constant for the time mean flow. In an unsteady sense, the passing rotor one wake causes variations of total temperature from time 0.5Tp to 0.9Tp.


## 4. Conclusion

Unsteady simulations of the CRAC compressor have been conducted using the FINE/Turbo 3D viscid turbulent flow solver with the Nonlinear Harmonic method. The aspiration slot and plenum were included in the computational domain. Adiabatic and isothermal thermal boundary conditions were simulated. The isothermal aspirated case best captures the trends of the experiment, especially near the endwalls.From these simulations, a detailed unsteady analysis of the CRAC aspiration slot and plenum has been completed allowing for a comprehensive understanding of the flow physics. The simulation agrees with experimental results that show the aspiration flow path choking at 0.5% of the inlet mass flow, as opposed to the design intent of 1%. Previously, the slot was believed to choke, thereby limiting the mass flow. However, this study shows that the limiting areas of the aspiration flow path are in the plenum. In addition, the unsteady interaction between the rotor one wake and the slot has been presented. Flow variations along the span and across the width of the slot were captured. The aspiration flow path has been shown to possess a complex 3D unsteady flow field, which is highly dependent on the characteristics of the main flow path.



---

*Source: 857616-2013-01-30.xml*