Time-space collocation method

The semi-discretized integral form of Navier-Stokes equations can be written as

where V is the volume of control cell, Q is the vector of conservative variables, R(Q) is the residual vector which contains the spatial discretized convective flux, viscous flux, and the source terms.

For the temporal periodic flows in turbomachinery, the conservative variables Q can be approximated by a truncated temporal Fourier series

where Q¯ is the time-averaged quantity, Nt is the number of temporal harmonics to be retained. Ati and Bti are the temporal Fourier coefficients of the ith unsteady component with frequency ωti, vi is the wave speed of the unsteady perturbation, λi is the wave length of the unsteady perturbation, and vp is the rotating speed of reference coordinates. The solution of Equation 2 is obtained at 2Nt+1 time instants for a reference blade passage.

In order to consider the variation of time averaged flow field in different blade passages, the spatial harmonics are introduced. The time-averaged quantity Q¯ can be expanded by a spatial Fourier series

where Q^ is the time-space averaged quantity, Ns is the number of spatial harmonics being retained.Asj and Bsj are the spatial Fourier coefficients of the jth harmonic with wave number ωsj, λb is the pitch of a single passage, mλb denotes the circumferential position of the considered blade passage. Combing Equations 4 and 2, the conservative variables can be approximated by both temporal and spatial Fourier series and written as

In addition, Equation 6 can be written in a more compact form as

where N=Nt+Ns is the total number of temporal and spatial harmonics being retained. It can be seen that Equation 7 has the same form as Equations 2 and 7 can be reduced to Equation 2 when there is only temporal periodic unsteady perturbations. In order to retain the N harmonics, at least 2N+1 independent samples of different time and passage index should be used.

Time space sampling

In the time-space collocation method, one of the most important issue is the selection of time instants and passage index. According to the previous studies, the the condition number of the inverse Fourier transform matrix of the selected time instants and passage index should reach a small value, otherwise the calculation will be divergent. An optimization algorithms named modified Gram Schmidt process has been used (Guédeney et al., 2013). In this paper, the passage index ranges from 1 to the blade counts in the studied row. When determining the time instants, the longest period of considered disturbances in the studied blade row is expanded by 500 times, and then firstly divided into K=2N+1 evenly spaced time levels. The condition number of the inverse Fourier transform matrix was calculated, and the number of samples K are increased until the condition number of the inverse Fourier transform matrix is below 5.0. In this way, although the final number of time-passage samples would be greater than K, but the resulted condition number of the inverse Fourier transform matrix is small enough to ensure the stability of the calculation.

Boundary conditions

The phase-shift boundary condition is adopted in the time-space collocation method to use computational domain of a single blade passage. As shown in Figure 1, The flow variable at interior cells along the lower periodic boundary can be given as

Figure 1.

The periodic boundaries and dummy cells.

The flow variable at dummy cell along the upper boundary can be updated as

where δk is the phase lag angle of the kth harmonic. The flow variables at dummy cells along the lower boundary can be updated using the same method, but with a different phase lag angle. At the rotor-stator interface, the flow field on the donor side is reconstructed using the phase-shift periodicity for the whole annulus, and then interpolate to the receiver side for the K time-space samples.

The time-space collocation method (TSC) has been implemented in an in-house CFD code named MAP (Multi-purpose Advanced Prediction code for fluid dynamics), which has been extensively validated in simulating turbomachinery flows (Ning, 2014). MAP solves the Reynolds-Averaged Navier-Stokes equations on structured grids, based on cell-centered finite volume method. The advective flux is discretized by LDFSS scheme which is one of the AUSM-family schemes and performs quite well in the viscous and shock wave-containing flows (Edwards, 1997). The MUSCL interpolation with minmod limiter is used to achieve high orders of spatial accuracy and guarantee monotone solutions close to discontinuities. The diffusive fluxes are evaluated using the conventional central differencing scheme. And the one-equation Spalart-Allmaras turbulence model is applied for turbulent flows (Spalart and Allmaras, 1992).

Test case

The propagation of inlet total pressure distortion through a two-stage compressor is simulated to verify the proposed TSC method. As shown in Figure 2, the compressor contains four blade rows, including one tandem stator blade row. The rotating speed is 13,266 rpm. Figure 3 shows the computational mesh at the 70% span that is generated by an in-house code named Turbomesh, which is primarily designed to automatically generate structured mesh for turbomachinery flow simulation (Ning, 2014). The O-3H mesh topology is used for a single blade passage, and the grid size for each blade passage and blade counts are listed in Table 1. The overall mesh size of single passage is about 2.25 million.

Figure 2.

Computational model of the compressor.

Figure 3.

Computational mesh at the 70% span.

Figure 4 depicts two total pressure distortion patterns representing the typical distortion patterns in experiments, which are specified at inlet boundary. The distortion patterns compose circumferential and radial non-uniformity, which are generated by inserting a plate in a straight pipe. In order to study the effect of distortion intensity on the compressor's performance, the Pattern-H contains larger regions with low total pressure than the Pattern-L. Meanwhile, the total temperature remains uniform and the circumferential flow angle is set to be zero at the inlet.

Figure 4.

Inlet total pressure patterns in 3D configuration. (a) Pattern-L (b) pattern-H.

Quasi-3D flow at 70% span

Firstly, a quasi-3D section of 70% span of this compressor is simulated for fast comparison purpose. And the quasi-3D computational domain is depicted in Figure 5, including two monitoring stations named st1 and st2, respectively. Figure 6 depicts the total pressure distortion profile specified at the inlet which is extracted from the Pattern-L as displayed in Figure 4, and the total pressure profile reconstructed by different harmonics of Fourier series are also plotted. It implies that the total pressure distribution approximated by only considering the first three harmonics could matches the specified inlet total pressure distortion pattern very well.

Figure 5.

Computational domain at the 70% span.

Figure 6.

Total pressure circumferential distribution at the inlet.

One of the purposes of this work is to identify the influence of harmonic orders adopted in the TSC method, hence the TSC simulation with the first three, five, and seven harmonics of the inlet total pressure distortion considered and without considering blade row interactions were performed. Besides, the extra disturbances induced by adjacent blade row interactions was considered by using the first three harmonics of the blade passing frequency to further assess its impact on the distortion transfer process. As a comparison, an URANS simulation of the inlet total pressure distortion was also performed. The number of physical time steps per revolution was set to be 950, and 20 sub-iteration steps in one physical time step was adopted. The five quasi-3D simulation cases have been carried out on a personal computer with 24 cores (3.0 GHz). To investigate the computational time consuming of TSC as compared to URANS, a simple speed up factor can be defined as following

where CURANS and CTSC are the CPU time consuming by the URANS and TSC method respectively.

Firstly, the total pressure ratio and isentropic efficiency characteristic of the compressor under clean and distorted inlet flow are shown in Figure 7. There is no obvious differences between the overall performance calculated by the TSC method with different harmonics considered. At the same mass flow rate, the isentropic efficiency calculated by the TSC method deviate from the results of URANS simulation by only 0.3%. When the additional three time harmonics are employed to consider the rotor-stator unsteady interactions, the efficiency characteristic is very close to the results of URANS simulation.

Figure 7.

Comparison of performance characteristic with clean and distorted flow.

Figure 8 depicts the efficiency convergence history of these five simulations. Notably, the convergence history of URANS simulation only shows the results at the physical time steps. The URANS simulation reaches a well-posed periodical state after three revolutions, and the actual number of iteration steps is as high as 57,000. As expected, the TSC simulations need much fewer iteration steps to approach a quasi-steady state. As shown in Table 2, compared to the reference URANS simulation, the TSC simulations consume much less CPU time. And the TSC gains about two orders of speed-up when three spatial harmonics are considered.

Figure 8.

Convergence history of the isentropic efficiency.

Figure 9 depicts the circumferential distributions of the nondimensional axial velocity, static pressure, and swirl angle at two positions upstream of Rotor 1, and the location of these two axial planes are shown in Figure 5. The compressor faces an nonuniform axial velocity distribution around the annulus when subjected to inlet total pressure distortion. The TSC results are quite compatible with those obtained by URANS simulation. As can be seen, three spatial harmonics adopted in the TSC method are sufficient to capture the distortion flow patterns upstream of the compressor. And adopting five or seven spatial harmonics in the TSC method only slightly narrows the distortion region, which is to be expected when high order harmonics are contained.

Figure 9.

Normalized axial velocity, static pressure and swirl angle distribution at two axial position.

The distorted region contains lower mass flow, which pushes the local operating point to higher pressure ratio, inducing a static pressure field imbalance upstream of the compressor. The gradient of static pressure drives the flow moves from undistorted region to distorted region, and forces the mass flow distribution to be more uniform, that is clearly visible from the axial velocity and swirl angle distribution at st2. The same flow redistribution is revealed by TSC solution with different harmonic orders as compared to the reference URANS solution.

The instantaneous entropy contours are displayed in Figure 10, which can illustrate the transportation of distortion within the compressor blade rows. For the reference URANS simulation, the distortion region with high entropy has been shifted by a phase angle in the rotation direction at the outlet of compressor. As previously mentioned, the compressor must confront co-swirl and counter-swirl flow when it enters and exits the distortion region (as shown in Figure 9). As a result, the blade's local operating point is different around the annulus, which causes the blade loading to be varied across the blade passages. Thick wakes generated by the high loading blades in the distorted region interact with the distortion flows and expand the high entropy region.

Figure 10.

Instantaneous entropy contours with distorted flow.

The phase shift phenomena is revealed by the TSC solution with three and seven spatial harmonics, but it also tends to narrow the compressor's high loss region downstream of compressor. The wake convection can be resolved when adjacent rotor-stator interactions are taken into account, producing roughly the same range of high entropy region compared to the URANS solution. This also explains why the efficiency decreases when three harmonics of rotor-stator interactions are considered.

The rotor blade travels in and out of distorted region periodically, inducing the high static pressure perturbations on the blade surface. The instantaneous static pressure perturbation field is depicted in Figure 11. When the rotor enters into the distorted region, the co-swirl flow tends to reduce the incidence, resulting in lower static pressure rise in the blade passage. However, when the rotor leaves the distorted region, the counter-swirl flow tends to rise the incidence, which leads to higher static pressure rise. The above static pressure perturbation phenomenon around the annulus is captured by the TSC simulation with three and seven spatial harmonics. It is important to note that the high static pressure perturbations in Stator 1 has not been captured unless three extra harmonics of rotor-stator interactions were contained in the TSC simulations. This demonstrates that considering the adjacent rotor-stator interactions is important to accurately simulate the pressure perturbation, but not so important for the resolution of distortion transfer.

Figure 11.

Instantaneous static pressure perturbation distributions.

The instantaneous total pressure perturbation distributions are shown in Figure 12 and the instantaneous static pressure distributions on some blade surfaces in Rotor 1 are shown in Figure 13. As depicted in Figure 12, Blade-1 is in the undistorted region, Blade-6 is entering the distorted region, Blade-10 is in the middle, and Blade-15 is leaving out. It is visible that the incidence angle of the TSC simulations is larger compared to URANS, thus the static pressure at leading edge of suction side is a little lower. However, it can capture the shock position and static pressure in other regions precisely. And the static pressure distribution on Rotor 1 blades does not be improved when three harmonics of stator 1 potential interactions were considered. Figure 14 depicts the instantaneous static pressure distribution on blade surface in Rotor 2. As depicted in Figure 12, Blade-1 is in the undistorted region, Blade-11 is entering the distorted region, Blade-16 is in the middle, and Blade-23 is leaving out. Similar with the results in Rotor 1, the TSC simulation only with inlet total pressure disturbance tends to lower the static pressure at leading edge, and there is a deviation of shock position from the URANS solution. For the Rotor 2 blades, they feel the wakes and pressure disturbance from the upstream Rotor 1 and Stator 1 blades, together with the potential disturbance from the downstream Stator 2 blades. The resolution of the static pressure at leading edge and shock position can be improved when three temporal harmonics resulted from rotor-stator interactions are considered. However, minor deviations still exists compared to the URANS simulation due to the nonlinearity of these unsteady interactions.

Figure 12.

The blade number and instantaneous total pressure perturbation distributions.

Figure 13.

The instantaneous static pressure of blade surface in Rotor 1.

Figure 14.

The instantaneous static pressure of blade surface in Rotor 2.

The instantaneous distribution of total pressure and total temperature at three axial positions are displayed in Figures 15 and 16, respectively. The rotor wakes are clearly visible. The total pressure distortion interacts with the blades of Rotor 1 causing the blade loading varied across the annulus, thus producing obvious total pressure and total temperature distortions downstream of Rotor 1. Figure 17 shows the instantaneous stagnation enthalpy rise downstream of Rotor 2, in which the higher-order harmonics have been filtered out. It can be seen that the stagnation enthalpy rise is higher in the distorted region. And it is demonstrated that the blades of Rotor 2 in distorted region contribute more work than blades in the undistorted region, thus producing a more uniform total pressure distribution at the outlet of Rotor 2. However, the total temperature distortion is obvious at the outlet of Rotor 2. Although these TSC simulations capture the primary total pressure and total temperature distribution patterns, there are small deviations that need further investigation.

Figure 15.

Normalized total pressure distribution at three axial positions.

Figure 16.

Normalized total temperature distribution at three axial positions.

Figure 17.

Instantaneous stagnation enthalpy rise downstream of Rotor 2.

The harmonic spectrum of the total temperature and total pressure distribution downstream of Rotor 1 and 2 are obtained by performing Discrete Fourier Transform (DFT), and the amplitude is normalized by the local averaged values. The majority of perturbation energy is occupied by the first two harmonics, as shown in Figures 18 and 19. The amplitude of the first three harmonics obtained by the TSC with three spatial harmonics is close to the URANS solution. And the resolution accuracy is improved a little when additional spatial or temporal harmonics are retained.

Figure 18.

Normalized amplitude of harmonics of total pressure downstream of Rotor 1.

Figure 19.

Normalized amplitude of harmonics of total temperature downstream of Rotor 1.

Compared to Figure 18, the amplitude of the first two harmonics has been significantly reduced (as shown in Figure 20), indicating the total pressure profile downstream of Rotor 2 becomes more uniform. However, the generated total temperature distortion indicated by the first three harmonics remain essentially unchanged, as shown in Figure 21. Due to the interactions of wakes and potential disturbance, the amplitude of some perturbations with high wave number cannot be ignored. The TSC with three harmonics of adjacent rotor-stator interactions can capture some perturbations with high wave number, but the improvement is not significant.

Figure 20.

Normalized amplitude of harmonics of total pressure downstream of Rotor 2.

Figure 21.

Normalized amplitude of harmonics of total temperature downstream of Rotor 2.

Three dimensional configuration

The two stage compressor subjected to inlet total pressure distortion is studied in this section. Firstly, the steady simulation with clean inlet condition was carried out. The normalized steady characteristic shown in Figure 22 matches the experimental results very well. Then the TSC with only five spatial harmonics was used to simulate the compressor flow subjected to inlet total pressure distortion, as denoted as Pattern-L and Pattern-H (in Figure 4). For the purpose of comparison, the URANS simulations with these two distortion patterns were also performed. As similar with quasi-3D case, the number of physical time steps per revolution was set to be 950, and 20 sub-iteration steps in one physical time step was adopted. The five 3D simulation cases have been carried out on a service computer with 112 cores (3.8 GHz).

Figure 22.

Characteristic of two stage compressor.

As depicted in Figure 22, the choke mass flow decreases obviously when subjected to inlet total pressure distortion. At the same mass flow rate, the total pressure ratio also decreases but the isentropic efficiency rises a little. The total pressure distortion induces the compressor to stall at higher mass flow rate and a larger distortion intensity results in more degradation of performance. The performance results of TSC agree well with that obtained by URANS at the peak efficiency working point. There is a slight deviation in the mass flow rate, which is only 0.25%. And in the simulation case with Pattern-L, the URANS simulation obtains a well-posed periodic state after three revolutions, but TSC simulation only need 1,000 iteration steps to approach the quasi-steady state. As shown in Table 3, the TSC simulation consumes much less CPU time, and gains about an order of speed-up. Near the stall boundary, two URANS simulation cases with Pattern-L have been performed. The physical time steps per revolution and sub-iteration steps were consistent with these used at the peak efficiency working point, and it needed to simulate four revolutions to obtain a well-posed periodic state. The performance results of TSC are also very close to that of URANS results at the near stall point. The TSC tends to underestimate the total pressure ratio, but the percentage deviation is less than 1.2%. Therefore, the TSC method can obtain relatively accurate results in estimating the performance characteristics of this compressor with inlet distortion.

The distortion transfer phenomenon in the compressor is three dimensional, that is much more complicated than the quasi-3D flow analysis in the above section. The peak efficiency working point with Pattern-L is chosen for detail comparisons between the TSC and URANS solutions. The instantaneous mass flux normalized by inlet averaged value and static pressure normalized by the inlet total pressure at undistorted region are shown in Figure 23. It is obvious that the distorted region contains lower mass flux (region A), thus lower static pressure is induced (region B), especially around the hub (region C). The static pressure gradient drives the flow of undistorted region migrates into the distorted region, creating the co-swirl (region D) and counter-swirl (region E) flow regions, as shown in Figure 24. In this figure, the positive swirl angle represents the opposite direction of the rotor rotation. Meanwhile, the fluid in the hub (region H) convects radially to the higher span sector. From Figures 23 and 24, it can be seen that the solution of TSC with five spatial harmonics reveals the complicated flow redistribution upstream of Rotor 1 at st2 as accurate as the URANS solution.

Figure 23.

Comparisons of normalized mass flux and static pressure upstream of Rotor 1 at st2. (a) URANS (b) TSC (five spatial harmonics).

Figure 24.

Comparisons of swirl angle and radial angle upstream of Rotor 1 at st2. (a) URANS (b) TSC (five spatial harmonics).

Figure 25 depicts the instantaneous normalized total pressure and total temperature at the downstream side of the R1-S1 interface. Since the disturbance of rotor-stator interactions were not contained in the TSC simulations, the wakes of rotor blades are invisible at the downstream side of the R1-S1 interface in the TSC solution. Due to the interactions of distorted flows and rotor blades, the nonuniform total pressure and total temperature patterns are consequently generated downstream of Rotor 1. When the blades of Rotor 1 travels into the distorted region, it face the co-swirl flow and the incidence as well as work input were reduced, resulting in lower total temperature and total pressure (region A). As shown in Figure 26, the circumferential total pressure distortion is attenuated significantly downstream of Rotor 2, but the total temperature distortion still remains. It can be shown that due to the radial migration of the flow in the hub regions, the radial total pressure and total temperature distortions are obvious downstream of Rotor 2. It can be seen from Figures 25 and 26, the TSC solution with five spatial harmonics considered captures the similar complex flow patterns as the reference URANS solution, thus the performance degradation caused by inlet distortion are well estimated by the TSC method.

Figure 25.

Comparisons of normalized total pressure and total temperature downstream of Rotor 1. (a) URANS (b) TSC (five spatial harmonics).

Figure 26.

Comparisons of normalized total pressure and total temperature downstream of Rotor 2. (a) URANS (b) TSC (five spatial harmonics).

Due to only the spatial harmonics are considered in TSC simulation, the comparisons of instantaneous results indicate some differences in the flow details captured by TSC and URANS. The time-averaged results of URANS in the near stall working point are further compared to the results of TSC in Figure 27. It can be seen that the TSC captures almost the same time-averaged nonuniform total pressure and total temperature patterns downstream of Rotor 1 as the URANS. Figure 28 displays the total pressure and total temperature distributions downstream of Rotor 2. The distribution patterns are similar, but the regions with lower total pressure and total temperature are broader in the solution of TSC than that of URANS, which would be caused by the neglecting of the rotor-stator interactions in the TSC simulation.

Figure 27.

Comparisons of time-averaged results downstream of Rotor 1. (a) time-averaged results of URANS (b) TSC (five spatial harmonics).

Figure 28.

Comparisons of time-averaged results downstream of Rotor 2. (a) time-averaged results of URANS (b) TSC (five spatial harmonics).

Symbols

Acronyms

Superscripts and Subscripts