Flow governing equations

The unsteady Euler equations for an ideal gas in a two-dimensional Cartesian coordinate system have been used in the investigation. They are written in a conservative and differential form as

Q is the conservative variable vector and F,G are convective flux vectors in x and y directions, respectively. p,ρ,u,v,et,ht denote pressure, density, velocity components, total energy, and total enthalpy. For a perfect gas

To simplify the following presentation, Equation (1) is rewritten as the following semi-discrete form

Harmonic balance method

The harmonic balance method is a computationally efficient reduced-order technique for time-periodic flows within turbomachinery (Hall et al., 2002). Its analysis time cost can be 1-2 orders of magnitude lower compared with the dual time-stepping method using a whole annulus computational domain for turbomachinery applications. With the harmonic balance method, the time-dependent flow variables are approximated using a truncated temporal Fourier series

where Q¯ is the vector of time-averaged conservative variables, QA,n, and QB,n denote the corresponding vectors of Fourier coefficients for the nth harmonic, N is the number of frequency components to be retained in an analysis, and ωn=nω is the angular frequency of the nth frequency component. For implementation, Equation (3) can be written in a discrete matrix form as

where Q is a column vector containing the conservative variable vectors at 2N+1 different time instants, Q~ is the corresponding Fourier coefficient vector. Their definitions are given as follows

and D is given by D=D⊗I, where D is the inverse discrete Fourier transform matrix, ⊗ is the Kronecker product, and I is an identity matrix with the size of nvar×nvar. For two-dimensional Euler equations, the number of conservative variables nvar=4. The inverse discrete Fourier transform matrix D is given as

Accordingly the time derivative term of the conservative variable vectors at 2N+1 time instants can be approximated by

where Dt=Dt⊗I and Dt is given by

The matrix E is defined as the time spectral operator, which turns the physic-time derivative term into the time spectral source term. Substituting Equation (6) into Equation (2), the harmonic balance equation system with a pseudo-time derivative term is given by

The pseudo-time derivative term is introduced to facilitate the use of time-marching methods to solve the harmonic balance equations.

Runge-Kutta scheme and LU-SGS method

The time integration is accomplished through a classic four-stage Runge-Kutta scheme together with the LU-SGS method. The LU-SGS method was first proposed by Yoon and Jameson (1988) as a time integration method. It was later used as a residual smoother at each stage of a Runge-Kutta scheme to improve stability and convergence rate (Rossow, 2006). The solution update at the kth stage of a Runge-Kutta scheme is defined as

where Q(0)=Qn−1, Qn=Q(K), and ΔQ is a column vector containing the solution increment at 2N+1 time instants. For the explicit scheme, the solution increment can be obtained directly as follows

While the explicit scheme is conditionally stable, the stability and convergence rate of the solution can be further improved by employing the implicit method, the LU-SGS method, as a preconditioner of the Runge-Kutta method. Subsequently, the solution increment is obtained by solving the following equations

where ∂R/∂Q=I2N+1⊗∂R/∂Q is a block diagonal matrix and I=I2N+1⊗I is an identity matrix. I2N+1 is the identity matrix of size (2N+1)×(2N+1) and ∂R/∂Q is the Jacobian matrix of the residual vector R(Q), which only contains the contribution of the convective flux for simplification. The pseudo-time step Δτ is calculated by the ratio of a cell volume to the sum of spectral radii of flux Jacobian matrices. In Equation (11), ε denotes the relaxation factor, αk denotes the kth Runge-Kutta coefficient and σ is the Courant number. The combination of ε, αk, and σ plays a significant role in determining the stability and damping behavior of a given scheme. Ideally, these coefficients should be optimized by using stability analysis. To solve Equation (11), an approximate LU decomposition is used to replace the left-hand side of the equation, given as follows

here B=I2N+1⊗B, L=I2N+1⊗L, and U=I2N+1⊗U, where B is a scalar, L is the lower operator, and U is the upper operator, given by

where matrices Am+ and Am− with m being i or j are defined as follows

where Am is the Jacobian matrix of convective flux through a cell face, given as

in which Sx and Sy are the bounding surface projections in the x and y coordinate directions, respectively. rm is the spectral radius of the Jacobian matrix Am. The shift operators Ei±l and Ej±l are defined as

More details about the derivation of the operators BLU can be found in Wang and Huang (2017).

So far the Jacobian matrix of the time spectral source term is not included on the left-hand side of Equation (11), implying that the time spectral source term is integrated explicitly. As the maximum grid-reduced frequency becomes high, the maximum allowable Courant number will be limited or even solution instability can occur in the harmonic balance equation system.

Block Jacobi method

To mitigate the impaired solution stability caused by the time spectral source term, it needs to be integrated implicitly together with other terms as follows

However, the exact inversion of the system Jacobian matrix of [I+εσΔτ/Ω(∂R/∂Q+ΩE)] is quite complicated, as E is a dense matrix. To avoid this complexity, the block Jacobi(BJ) method was proposed by Sicot et al. (2008) to approximately invert the Jacobian matrix of the time spectral source term through iterations. Wang and Huang (2017) combined the block Jacobi method with the LU-SGS method, leading to the LU-SGS/BJ method:

where M is the number of block Jacobi iterations and the initial solution increment ΔQ0 is set to 0. At least two Jacobi iterations are required to have the time spectral source term integrated implicitly. The advantage of the block Jacobi method is its easy implementation and minimal alterations to existing codes with a baseline LU-SGS method.

When using the LU-SGS/BJ method to solve the harmonic balance equation system, it is crucial to consider the influence of the time spectral source term in scalar B. Otherwise, the full potential of the LU-SGS/BJ method cannot be exploited. The scalar B can be modified as follows

where NΩω is the spectral radius of the time spectral source term.

Boundary conditions

In general, there are four kinds of boundary conditions used in this numerical analysis: subsonic inlet, subsonic outlet, slip wall, and periodic boundary conditions. At a subsonic inlet, the total pressure, total temperature, and flow angle are specified, while the static pressure is extrapolated from the interior domain to the boundary. Additionally, one can also extrapolate the one-dimensional Riemann invariant to the inlet boundary. The Riemann invariant R− is defined as follows

where V is the velocity magnitude and c is the speed of sound.

The well-posedness of a boundary condition is important for a numerical analysis as it can affect solution convergence characteristics and even solution accuracy. Extrapolating the Riemann invariation at an inlet boundary has been found to offer better convergence properties for nozzle flows (Jameson and Caughey, 2001). Nevertheless, there is currently a lack of stability analysis supporting this conclusion. This paper will examine the better convergence characteristics of extrapolating the Riemann invariant through stability analysis.

At a subsonic outlet boundary, the static pressure is specified, while the density and momentum in the x and y directions are extrapolated from the interior domain to the boundary. To obtain a meaningful unsteady solution for a Laval nozzle and the NASA rotor 37, a sinusoidal static pressure at the outlet is prescribed

where p¯ is the time-averaged static pressure, Δp is the amplitude of static pressure perturbation with the value of Δp=0.1p¯, and ω is the angular frequency of unsteadiness.

von Neumann method

Though Swanson et al. (2007) performed a stability analysis of the RK/LU-SGS methods for two-dimensional steady Euler equations, there has been no such effort of the RK/LU-SGS methods for solving the harmonic balance equation system in the open literature. The stability and damping behaviours of the RK/LU-SGS methods will change dramatically due to the time spectral source term in the harmonic balance equation system. To improve the solution stability of a harmonic balance equation system, the block Jacobi method is often used. Therefore, there is also a need to perform a stability analysis of the RK/LU-SGS/BJ method for solving the harmonic balance equation system.

In the von Neumann analysis, we assume the Euler equations are solved in a finite domain with periodic boundary conditions. Let us define a Cartesian grid consisting of im×jm cells, where the length and width of each cell are denoted by Δx and Δy, respectively. Q′i,j denotes the numerical error at the cell (i,j) and can be approximated using a truncated spatial Fourier series as follows

where J=−1 is the imaginary unit and Q^n,m is the amplitude of a spatial harmonic. φx and φy are the phase angles in the x and y directions of the spatial harmonic, covering the range of [−π,π] respectively.

The numerical error Q′i,j satisfies the linear form of Equation (14), given as

where L=I2N+1⊗L is the block diagonal matrix and L is the linearized discrete residual operator, consisting of the convective part Lc and dissipative part Ld as follows

As the Cartesian grid is uniform, Lc and Ld can be given by

where AF=(∂F/∂Q) and AG=(∂G/∂Q) are the Jacobian matrices of the convective fluxes in the x and y direction. ε(4) is the linear dissipative coefficient. Λi and Λj are the modified spectral radii

in which ΛF and ΛG are the spectral radii of the Jacobian matrices AF and AG.

Substitute a single spatial harmonic Q^e[J(iφx+jφy)] from Equation (17) into Equation (18), we can obtain

where L^ is the transformed linear residual operator, given by

where

and L^ and U^ are the transformed lower and upper operators, given by

Through M Jacobian iterations, the transformed increment of the solution can be simplified as

where P^ is the transformed preconditioner. For the LU-SGS/BJ method,

which can be obtained as follows

in which the initial value P^0 is set to 0. P^LU is the transformed preconditioner of the LU-SGS method, given by

When M=1, the preconditioner P^=P^LU, meaning that only the LU-SGS method is applied. When M≥2 the block Jacobi method is applied to deal with the time spectral source term implicitly. As the K-stage Runge-Kutta method is used for pseudo-time integration, the updated numerical error in the next pseudo-time step can be written as

where the amplification matrix GRK is given by

in which Z is the Fourier symbol of the spatial operator, defined by

Typically, the Fourier symbol Z is complex. When the phase angles φx and φy vary from −π to π, the Fourier symbol forms a locus on the complex plane. The locus will scale up with an increase in the Courant number or a decrease in the relaxation factor. The magnitude contour of the amplification factor can also be plotted on the complex plane according to Equation (24) with the contour level ranging from 0 to 1. The contour lines, often referred to as the stability region, depend on Runge-Kutta coefficients. The amplification factor measures the growth of the numerical error over iteration. For a stable scheme, the numerical error must not grow with iteration. So the maximum modulus of the amplification factor must be less than or equal to one for all phase angles φx and φy varying from [−π,π]. In other words, the locus of the Fourier symbol Z in the complex plane should be inside the stability region of the Runge-Kutta scheme. The maximum allowable Courant number and the minimum allowable relaxation factor thus depend on the interval between the locus of the Fourier symbol and stability domain along the imaginary axis as well as the negative real axis.

Matrix method

The matrix method, also known as the eigenvalue analysis, is a more general approach to stability analysis. It allows for considering various boundary conditions, grid nonuniformities and flow field nonuniformities. In the case of a non-uniform mesh, the total residual function R(Q) for each control volume consists of two parts: the convective part Rc(Q) and the artificial dissipative part Rd(Q). They are expressed as follows

where the convective fluxes through a bounding surface with the index i+(1/2), omitting the index j, are defined by

Fluxes at a bounding surface are calculated by simply averaging their values on two sides of the bounding surface. For example, Fi+(1/2)=1/2(Fi+Fi+1). The artificial dissipative fluxes are given by

where Λi+(1/2)=(1/2)(Λi+Λi+1), and ε(2), ε(4) are the coefficients of the nonlinear and linear dissipative terms, respectively. The nonlinear dissipative term is designed to provide dissipation at discontinuities such as a shock wave. The linear one is used to suppress spurious solutions, leading to a fast convergence to a steady state. The coefficients ε(2) and ε(4) depend on the distribution of static pressure from a solution and the constant parameters k(2) and k(4) (Jameson et al., 1981).

In the stability analysis based on the matrix method, the solution iteration needs to be written in a matrix form. For an explicit scheme, the solution increment is given by

where the index i=1,2,…,im and j=1,2,…,jm. The matrix form of Equation (29) can be written as follows

where (Q∗)(k−1) is the solution vector, ΔQ∗ represents the solution increment vector, and (R∗)(k−1) is the residual vector containing the convective fluxes and dissipative fluxes at each cell under various boundary conditions. (Δτ/Ω)∗ is a diagonal matrix that contains the ratios of pseudo-time step to volume and (ΩE)∗ is a block diagonal matrix containing the time spectral source terms at each cell. Let us define the column vector Q′∗ containing the numerical errors at each cell. The numerical error Q′∗ satisfies the linear equation system, given by

where L∗=(∂R∗/∂Q∗) is the system Jacobian matrix of the total residual R∗, which has been derived manually. In the matrix method, the impact of various boundary conditions on the stability of the solution can be considered as follows.

Firstly, the periodic boundary conditions in the x and y directions are applied in the matrix method, as illustrated in Figure 1a. The use of periodic boundary conditions in the y direction introduces non-zero off-diagonal elements to the system Jacobian matrix L∗. For the x direction, the periodic boundary condition is implemented in the submatrix along the main diagonal of L∗. Similarly, this introduces non-zero off-diagonal elements to this submatrix. The inlet, outlet, and slip wall boundary conditions are taken into account in the stability analysis by modifying relevant off-diagonal elements of L∗ and the submatrix along the diagonal as shown in Figure 1b.

Figure 1.

Visualization of the Jacobian matrix of the spatial discrete operator for the 2-D Euler equation system. Nonzero elements are displayed as filled blue dots. (a) Periodic B.C. in x and y directions. (b) Inlet,outlet B.C. in x direction and slip wall B.C. in y direction.

The Jacobian matrix L∗ is usually frozen with an initial or fully converged steady solution. After the Jacobian matrix is given, the updated numerical error at the next pseudo-time step can be written as

where GRK∗ is the amplification matrix of the equation system, given by

where I∗ is an identity matrix of size [nvar⋅im⋅jm⋅(2N+1)]×[nvar⋅im⋅jm⋅(2N+1)] and the transform matrix Z∗ is defined as

In the matrix method, if the spectral radius of GRK∗ is not larger than one, then the combination of the pseudo-time integration scheme and the spatial discretization scheme will be numerically stable.

The stability of an explicit Runge-Kutta scheme with different boundary conditions

First, the stability of an explicit Runge-Kutta scheme for solving steady Euler equations with periodic boundary conditions is examined using both the von Neumann method and the matrix method. The computational domain has a length of 0.2 meters and a height of 0.08 meters as shown in Figure 2. It is discretized using a uniform Cartesian grid with 33 points in the x direction and 9 points in the y direction. The resulting grid aspect ratio is AR=0.625. A uniform flow field with a Mach number of Ma=0.5 and a flow angle of β=0∘ is used as a baseline for linearization. Periodic boundary conditions are assumed in each direction for the von Neumann analysis and matrix analysis. Figure 3 presents the stability analysis results obtained from these two methods.

Figure 2.

The computational domain and grid for stability analysis.

Figure 3.

The locus of the Fourier symbol and the eigenvalue spectrum of the transform matrix with periodic boundary conditions, uniform flow fields and grids (AR=0.625,Ma=0.5,β=0∘,k(4)=0.21,σ=3) and the magnitude contours of the amplification factor of the classic four-stage Runge-Kutta scheme.

The locus of the Fourier symbol Z obtained by the von Neumann analysis aligns precisely with the distribution of eigenvalues of the transform matrix Z∗ obtained through the matrix method. The good agreement between the two methods verifies the correct implementation of the matrix method. Figure 3 also presents magnitude contours of the amplification factor of the classic four-stage Runge-Kutta scheme with coefficients α=[1/4,1/3,1/2,1]. The contour levels are |GRK|=0.002,0.01,0.05,0.20,0.35,0.50,0.65,0.85,1.00. The outermost contour corresponds to amplification factors with a magnitude of 1. For a given spatial discretization scheme, the locus of the Fourier symbol or the distribution of the eigenvalues should lie within this outermost contour, referred to as the stability domain. It ensures that the maximum magnitude of the amplification factor is not greater than 1 so that the combination of the pseudo-time integration scheme and the spatial discretization scheme is numerically stable.

To assess the impact of various boundary conditions by the matrix method, these boundary conditions are categorized into four combinations. The first combination serves as a reference, which includes periodic boundary conditions in the x and y directions. The second combination involves inlet and outlet boundary conditions in the x direction and periodic boundary conditions in the y direction. The third combination includes wall boundary conditions in the y direction and periodic boundary conditions in the x directions; The fourth combination consists of inlet and outlet boundary conditions in the x direction and wall boundary conditions in the y direction. The eigenvalue spectra of these four combinations of boundary conditions are presented in Figure 4.

Figure 4.

The eigenvalue spectra of the transform matrix for various boundary conditions. (k(4)=0.021,σ=3). (a) First combination. (b) Second combination. (c) Third combination. (d) Fourth combination.

Comparing Figure 4a and b, it can be seen that the inlet and outlet boundary conditions change the distribution of the eigenvalues considerably in the complex plane. Figure 4c presents the eigenvalue spectrum of the third combination, which is quite similar to that in Figure 4a. It implies that the impact of the wall boundary conditions is smaller than that of the inlet and outlet boundary conditions. Figure 4d presents the eigenvalue distribution of the fourth combination, which resembles that of the second one. In conclusion, different combinations of boundary conditions result in varying distributions of eigenvalues in the complex plane, leading to different stability and damping properties. Among these boundary conditions, inlet and outlet boundary conditions have the most significant impact on the eigenvalue spectrum, leading to enhanced stability. Note the impact of boundary conditions diminishes with the increase of grid points.

Next, the impact of different baseline flow fields and inlet boundary conditions on solution stability is investigated by the matrix method. The stability analysis is performed on a Laval nozzle as shown in Figure 5a. A mesh with 65 grid points in the x direction and 17 grid points in the y direction is used to discretize the domain. At the inlet, the total pressure Pt of 101,325Pa, the total temperature Tt of 288.15K, and the flow angle of 0∘ are given. At the outlet, the static pressure Pb of 80,000Pa is specified. The resulting steady solution is used as the baseline for linearization, as shown in Figure 5b. The flow field is transonic with a shock wave near the throat. The Mach number in front of the shock is 1.5. This indicates a strong non-linearity in the solution.

Figure 5.

The computational mesh and Mach number contours within a Laval nozzle obtained from solving the 2-D steady Euler equations. (a) mesh. (b) Mach number contour.

The stability analysis results are presented in Figure 6. Figure 6a and b present the eigenvalue spectra of the transform matrix obtained from the uniform flow field and the nozzle flow field with the fourth combination of boundary conditions, respectively. At the inlet boundary, the static pressure is extrapolated, denoted as inlet1. When using the steady solution of the nozzle as the baseline, the eigenvalues move along the real axis towards the negative direction, approaching the stability boundary. The resulting stability is reduced. To mitigate the instability caused by the nozzle flow field, the Riemann invariant extrapolation is applied at the inlet, denoted as inlet2. Figure 6c shows the eigenvalue spectrum obtained from inlet2 boundary condition. It can be seen that some leftmost eigenvalues in Figure 6b disappear in Figure 6c, indicating that the Riemann invariant extrapolation improves the solution stability and damping property.

Figure 6.

The eigenvalue spectra of the transform matrix for different baseline flow fields and inlet boundary conditions. (k(4)=0.021,σ=3). (a) uniform flow field, inlet1. (b) nozzle flow field, inlet1. (c) nozzle flow field, inlet2.

For an unsteady flow analysis using the harmonic balance method, the solution stability is compromised by the time spectral source term. In the von Neumann analysis or the matrix analysis, the time spectral source term introduces a pure imaginary value into the Fourier symbol or the eigenvalues, resulting in a shift of them along the imaginary axis in the complex plane, as illustrated in Figure 7. The extent of this shift is associated with a dimensionless parameter, the maximum grid-reduced frequency, defined as

Figure 7.

The locus of the Fourier symbol and the eigenvalue spectrum of the transform matrix (ωx=1,σ=1.4) and the contours of the magnitude of the amplification in the case of the classic four-stage Runge-Kutta scheme. (a) von Neumann method, periodic. (b) matrix method, inlet1. (c) matrix method, inlet2.

where ωmax is the maximum angular frequency of the unsteadiness, Ω is the cell volume, ri and rj are the spectral radii of the Jacobian matrix of convective fluxes in x and y directions respectively. The maximum grid-reduced frequency is a critical parameter with a substantial impact on the stability of the harmonic balance equation system. If the maximum grid-reduced frequency exceeds a certain value, the resulting locus of the Fourier symbol or the distribution of the eigenvalues will shift outside the stability region. To ensure the stability of the solution, a smaller Courant number has to be used to scale down the eigenvalues so that they stay within the stability domain of the time integration scheme.

To verify the results obtained from the stability analysis, the flow field within the nozzle configuration was analyzed by our in-house solver. A sinusoidal static pressure is prescribed at the outlet. The periodic unsteady flow induced by the pressure perturbation is solved through the harmonic balance method with a single harmonic. The maximum grid-reduced frequency varies from 0 to 8 though altering the frequency of the unsteadiness. Two subsonic inlet boundary conditions, named inlet1 and inlet2, are implemented in the solver. inlet1 represents the extrapolation of static pressure and inlet2 represents the extrapolation of the Riemann invariant. When the maximum grid-reduced frequency is set to ωx=0, the steady Euler equations are solved. The steady solution reveals that inlet2 has a higher maximum allowable Courant number than that of inlet1, as illustrated in Figure 8. However, the stabilizing effect of inlet2 does not hold for a harmonic balance solution when the maximum grid-reduced frequency is no less than 1. The results of the numerical experiments agree well with the stability analysis results.

Figure 8.

The maximum Courant number varying with the maximum grid-reduced frequency.

As the maximum grid-reduced frequency increases, the maximum allowable Courant number decreases quickly. Moreover, The maximum allowable Courant numbers with the two different inlet boundary conditions, obtained by numerical experiments, are in good agreement with the matrix analysis results. When it comes to the von Neumann analysis with the assumption of periodic boundary conditions, the predicted maximum allowable Courant numbers are slightly smaller than those obtained by the matrix method at some maximum grid-reduced frequencies, although they are still in good agreement with the results from the numerical experiments.

The stability of the LU-SGS method with the block Jacobi method

The explicit scheme encounters stability issues when it is used to solve the harmonic balance equation, due to the stiffness of the time spectral source term, especially with a large grid-reduced frequency. To stabilize a harmonic balance solution, two implicit methods are presented. The LU-SGS method was first proposed as an implicit preconditioning of the Runge-Kutta method to improve the stability of a steady solution. It can be applied directly to solve the harmonic balance equations without incorporating the Jacobian matrix of the time spectral source term on the left-hand side. It can also be combined with the Block Jacobi method to integrate the time spectral source term implicitly (Wang and Huang, 2017). To assess the stability of these two implicit methods, a von Neumann analysis is conducted.

In the stability analysis, a uniform flow with a Mach number of 0.5 and a flow angle of 0∘ is used as the baseline flow for linearization. At first, the stability of the LU-SGS method is investigated for solving the steady Euler equation system (the corresponding grid-reduced frequency ωx=0), the Courant number σ=1000 and the relaxation factor ε=0.6. In Figure 9a, the Fourier symbol locus of the LU-SGS method takes a kidney shape in the complex plane, which is notably different from that of the explicit method. In addition, this locus lies within the stability domain of the four-stage Runge-Kutta scheme, even though the Courant number is increased to the order of 103. Therefore, compared to the explicit scheme, the LU-SGS method effectively enhances the stability of a steady solution.

Figure 9.

The Fourier symbol loci for different maximum grid-reduced frequencies (σ=1000,ε=0.6) and the magnitude contours of the amplification factor of the classic four-stage Runge-Kutta scheme. (a) ωx=0, LU-SGS. (b) ωx=0.2, LU-SGS. (c) ωx=0.5, LU-SGS. (d) ωx=0.5, LU-SGS/BJ.

For the harmonic balance equation system, the locus of the Fourier symbol of the LU-SGS method is influenced by the time spectral source term. When the maximum grid-reduced frequency is increased from 0 to 0.2, the locus near the imaginary axis moves to the right side, as shown in Figure 9b. When the maximum grid-reduced frequency reaches 0.5, the locus extends beyond the stability boundary of the classic four-stage Runge-Kutta scheme, shown in Figure 9c. Consequently, the magnitude of the amplification factor is larger than unity and the scheme becomes unstable. To stabilize the solution, the block Jacobi method is used with the LU-SGS method to deal with the time spectral source term implicitly. In Figure 9d, the Fourier symbol locus is restricted within the stability domain when the maximum grid-reduced frequency is set to ωx=0.5. Therefore, the LU-SGS/BJ method can effectively stabilize the solution of the harmonic balance equations.

To further examine the stabilization effect of the LU-SGS/BJ method, we consider a sufficiently large maximum grid-reduced frequency in the stability analysis. In Figure 10a, the Fourier symbol locus of the LU-SGS method lies outside the stability domain when the maximum grid-reduced frequency ωx is set to 2 and the Courant number σ is set to 1,000. However, the locus of the LU-SGS/BJ method is restricted inside the stability domain as shown in Figure 10b. Thus, the LU-SGS/BJ method can effectively stabilize a harmonic balance solution even when the maximum grid-reduced frequency is large, for which the LU-SGS method is unstable.

Figure 10.

The Fourier symbol loci for ωx=2,σ=1000,ε=0.6. (a) LU-SGS method. (b) LU-SGS/BJ method.

To verify the stability results obtained through the von Neumann method, the NASA rotor 37 is employed as a more realistic test case. The computational domain consists of one blade passage. The domain is discretized using an H mesh with 57 grid points in the radial direction, 49 grid points in the circumferential direction, and 137 grid points in the streamwise direction, among which 72 grid points are on the blade surface (one side only). At the inlet, the total pressure Pt of 101,811.36Pa, the total temperature Tt of 288.03K, and the two flow angles of 0∘ are given. At the outlet, the back pressure Pb of 115,000Pa is specified. The 3-D RANS equations, coupled with the Spalart-Allmaras (S-A) turbulence model, are solved by our in-house solver (Wang and Huang, 2017; Huang et al., 2018). The calculated mass flow rate is 20.84kg/s, the total pressure rate is 1.74 and the isentropic efficiency is 82.72%. The Mach number contours from a steady solution are presented in Figure 11.

Figure 11.

The Mach number contours of the NASA rotor 37 at 50% blade span from the solution of the 3-D steady RANS equation.

Then an unsteady analysis is conducted by varying the static pressure periodically in time at the outlet. The harmonic balance method is used to analyze the unsteady flow with a single harmonic. The maximum grid-reduced frequency varies from 0.2 to 0.5 by altering the frequency of the perturbation. The convergence history of the density equation residual for different implicit methods and maximum grid-reduced frequencies are shown in Figure 12a. Both the LU-SGS and the LU-SGS/BJ method can achieve good convergences when the maximum grid-reduced frequency is not larger than 0.2. However, when the maximum grid-reduced frequency exceeds 0.2, the LU-SGS method fails to achieve convergence with the Courant number of σ=1000 and relaxation factor of ε=0.6. In contrast, the solutions can converge effectively by using the LU-SGS/BJ method. Figure 12b illustrates the time cost comparison between the LU-SGS method and the LU-SGS/BJ method. As expected, the LU-SGS/BJ method requires more CPU time per iteration than that of the LU-SGS method alone. Based on the analysis, we conclude that the block Jacobi method is unnecessary when the maximum grid-reduced frequency is below 0.2. However, if the maximum grid-reduced frequency exceeds this value, the time spectral source term could compromise the stability of the harmonic balance equation system. In such cases, it is recommended to use the block Jacobi method to stabilize the solution so that a larger Courant number can be used to achieve fast convergence.

Figure 12.

The convergence histories of Res(ρ) solved using the LU-SGS method and the LU-SGS/BJ method (ε=0.6). (a) iteration. (b) CPU time.

Furthermore, the numerical experiments of NASA rotor 37 with the maximum grid-reduced frequency of 2 are conducted by increasing the unsteady frequency. The convergence historyies of the density equation residual for different implicit methods are shown in Figure 13a. It can be observed that the solution solved with the LU-SGS method can only converge with a very small Courant number, approximately equal to 1. However, with the LU-SGS/BJ method, the maximum Courant number can reach up to1000, or even higher, resulting in a fast convergence rate. A finer mesh of NASA rotor 37 is employed to confirm the convergence performance. The total number of grid points of the finer mesh is twice that of the coarse mesh. The convergence performance of the LU-SGS/BJ method remains effective, as shown in Figure 13b. Notably, the maximum grid-reduced frequency is reduced due to the reduction in cell volumes of the finer mesh.

Figure 13.

The convergence histories of Res(ρ) for the LU-SGS method and the LU-SGS/BJ method (ε=0.6). (a) Coarse mesh. (b) Fine mesh.

The LU-SGS/BJ method demonstrates a significant reduction in total CPU time, requiring approximately 23.5% less time compared to the LU-SGS method to achieve the density equation residual level of 10−7 for a coarse mesh. This acceleration is primarily attributed to the better stability of the LU-SGS/BJ method to allow for a much larger Courant number.