Abstract

The block-wise spectral reconstruction of turbomachinery flows has been successfully used in the past to reduce the computational cost of full-annulus turbomachinery simulations. The rationale of the method is based on an explicit representation of the long-wavelengths of the problem using a Fourier expansion of the blocks, which in turn contain the short-wavelength content of the spectrum. This work analyzes the accuracy of this approach and shows that this approximation retains not only the explicit Fourier representation of the long wavelengths and finite volume approximation of the short wavelengths but their interaction as well. This interaction due to the scattering associated with the nonuniformity of the flow is referred to as Tyler–Sofrin modes by turbomachinery practitioners and is highly relevant for many practical applications. The solution’s harmonic content is illustrated by resorting to the linear wave equation with varying propagation speed and the NASA rotor 67 fan.

Graphical Abstract Figure
Graphical Abstract Figure
Close modal

1 Introduction

The rapid growth of computer power, and advancements in numerical techniques have made computational fluid dynamics (CFD) methods powerful and essential tools for modeling turbomachinery flows. This development of computer power has enabled the simulation of the flow in full-annulus domains, although the computing requirements are significantly larger than those of single-passage simulations. Unsteady Reynolds-averaged Navier–Stokes (URANS) simulations can improve gas turbine designs by increasing the physical understanding of flow phenomena in complex systems. Modeling the unsteady aerodynamics of turbomachinery using full-annulus unsteady simulations with multiple blade rows, though feasible [1,2], is computationally intensive when using standard time-marching methods. Physical computational time often involves several engine rotations and time-marching methods require many time-steps to achieve accurate solutions. Therefore, developing efficient and accurate numerical methods for turbomachinery unsteady flows, which avoid the computation of the whole annulus domain, is an active research field.

Discrete Fourier series are an efficient means to represent periodic signals. Many turbomachinery flows exhibit temporal and/or spatial periodicity, and are well suited for their Fourier representation. Flow unsteadiness can arise from various sources be periodic or non-periodic. A truncated Fourier series is, in essence, a reduced-order model when the number of harmonics is limited to bound the computational effort and only some of the scales are resolved. The intrinsic filtering capability of a Fourier series can be regarded as a strength rather than a limitation when the most relevant scales are a priori known and are retained explicitly.

Extensive reviews of Fourier-based methods for turbomachinery applications can be found in Refs. [3,4], where a daunting array of modeling strategies are presented. The use of spatial Fourier methods for the computation of nonlinear flows in turbomachinery was first proposed by He [5]. In his original work, He used a block decomposition to study the steady 2D interaction between a fan outlet guide vane (OGV) and a Pylon. Moreover, simplified unsteady cases with two-dimensional [5] or axisymmetric [6] configurations undergoing longwavelength unsteady perturbations were expanded in spatial Fourier series. It was shown that the agreement between the spectral method and the full-annulus computation was good, with significant computational savings. He applied the same idea to film cooling configurations [7], dimples [8], and lined acoustic intakes [9]. Time-domain spatial Fourier methods have been used for general non-axisymmetric flows across multiple blade rows by Stapelfeldt and Di Mare [10], who extended the technique by exploiting the time/space periodicity of the flow and retaining circumferential perturbations due to interactions within the same frame of reference. More recently, He [11] also extended the methodology to unsteady flows of short temporal and spatial scales (e.g., self-excited trailing-edge vortices and turbulent disturbances), where a source term-based approach was adopted to facilitate a two-way coupling in terms of time-averaged flow solutions.

The authors have tested the computational efficiency of the block spectral method (BSM) in the past [12,13]. In all the cases in which the approach was employed to analyze the NASA Rotor 67 subject to inlet perturbations, the transient (not the final periodic solution) was reproduced correctly, provided that the number of spatial harmonics retained in the simulation was high enough. The convergence toward the periodic state is not spoiled using the block spectral methods. The computational efficiency is proportional to the reduction in the simulation’s points or blade passages. However, the block spectral method cannot compete in computational efficiency with harmonic methods that reduce the size of the computational domain and the time required to converge to the periodic state. Nevertheless, the block spectral approach is robust and can easily handle strongly separated flows.

The present work is aimed at understanding the modeling limitations of the block spectral method. The one-dimensional linear wave equation (LWE) has been chosen as the main vehicle for the analysis because of its simplicity. The approach is meant to be used in non-time periodic problems displaying a significant disparity of spatial scales. In turbomachinery problems, the short scale is dictated by the airfoil pitch, which gives rise to a fundamental spatial harmonic with a wavenumber associated with the blade count. This fundamental and its higher-order harmonics are resolved within the blade passage. This fundamental block is, on many occasions, modulated with a slowly varying spatial variation in the circumferential direction due to the air-foil interaction with upstream or downstream blade rows [5], or sometimes with long-wavelength spatial nonuniformities such as the engine intake [12].

Ideally, all the homologous points within the different blocks could be Fourier transformed to represent this spatial variation. This approach is not efficient since it would require performing as many Fourier transforms as the number of points per passage per time-step. Instead, if the baseline solver is fully explicit only the information of a few neighboring points is needed to advance a certain point in time. The same is valid on the circumferential boundaries of the block, which in this particular case represents the blade passage, and therefore, only the information of the neighboring nodes to this boundary needs to be Fourier-reconstructed to complete the numerical scheme.

It will be shown that the method can retain all the harmonics self-contained within the block and the slow varying harmonics reconstructed in the vicinity of the neighboring blocks using a Fourier series. Additionally, a third family of scattered or sideband spatial harmonics (also known as Tyler–Sofrin models) is captured by this approach. These modes are natively retained by the method. The ability of the method to retain the Tyler–Sofrin modes is not apparent.

In many non-time periodic practical problems, a few low-order spatial harmonics suffice to predict the solution accurately. This means that only a few blocks out of the tens or hundreds forming a full wheel of airfoils need to be solved with a remarkable reduction in the computational resources with respect to full-annulus simulations, giving the best compromise between accuracy and computational effort. The authors have verified the method using full-annulus simulation for rotating stall and force response cases in an isolated fan-stage [12]. It has also been modified recently to efficiently deal with a single spatial harmonic and then applied to assess the effect of inlet distortion on nonlinear fan stability [13].

The numerical formulation of the method is introduced first using the one-dimensional linear wave equation. This model problem is used to describe the harmonic content retained by the method since the spatial representation is a mixture of a finite volume method with a spectral reconstruction. Finally, the implementation of the so-called block spectral method (named by the authors sometimes as the passage-spectral method (PSM) when applied to the solution of turbomachinery blade rows) in an unstructured edge-based finite volume Reynolds-averaged Navier–Stokes solver is briefly described and applied to the inlet distortion of an isolated fan rotor.

2 One-Dimensional Linear Wave Equation Method

The one-dimensional Euler equations in conservative form can be written as
t(ρρvρE)+x(ρvρv2+pρvH)=0
(1)
where ρ is the density, v is the velocity, p is the static pressure, E=e+(1/2)v2 is the total energy, and H=h+(1/2)v2 is the total enthalpy. This system of equations has to be completed by an equation of state defining the thermodynamic properties of the fluid. Considering a perfect gas, p=ρRgT, e=cvT, and h=cpT. Equation (1) can be written in compact form as tU+xF(U)=0, where U={ρ,ρv,ρE}T is the vector of the conservative variables and F is the flux.
Let us consider the linearization of the one-dimensional Euler equations about a uniform steady flow U0. In this case, the flow can be expressed as U=U0+u(x,t), where u is a small perturbation (uU0). The linearized Euler equations, ut+(F/U)0ux=0, can be diagonalized and written in characteristic form as
t(p+ρ0a0vpρ0a0vpρa02)+[v0+a0000v0a0000v0]x(p+ρ0a0vpρ0a0vpρa02)=0
(2)
where a0 and v0 are the sound speed and the convection velocity of the mean flow, respectively. The resulting system of equations represents three independent linear wave equations for the characteristic variables, each of them with its own propagation velocity, c0(j). The jth characteristic equation can be written as ut(j)+c0(j)ux(j)=0, and therefore the problem reduces to show how to solve the LWE
ut+c0ux=0,x(0,2π),t>0
(3)
with a constant propagation speed, c0, subject to the initial condition
u(x,t=0)=u0(x)
(4)
and periodic boundary conditions at both ends of the computational domain
u(0,t)=u(2π,t)
(5)
The choice of these boundary conditions is based not only on its simplicity but on the fact that they are the natural boundary conditions in the circumferential direction for a full wheel of blades.

3 Motivation and Qualitative Description of the Method

The one-dimensional wave equation is an excellent vehicle to illustrate the whole purpose of the method due to its simplicity. The main idea of the method is sketched in Fig. 1 where the solution of the LWE has been initialized with a short wavelength, λc=2π/kc, disturbance superimposed with a long-wavelength, λL=2π/kL, and perturbation having the latter, a much longer wavelength than the former (λLλc). Nevertheless, such a condition does not restrict the method, since this hypothesis is not used in its derivation and the technique also works for λLλc. However, in such a case, the computational savings of the method would be much smaller. The wavenumber of the characteristic short-wavelength perturbation will be referred to as kc. The spatial harmonics of both perturbations propagate at the same speed and do not interact with each other since the problem is linear and non-dispersive. However, the longest spatial scale imposes a constraint on the minimum size of the computational domain. This problem frequently appears in turbomachinery.

Fig. 1
Comparison of the BSM and the analytical solution (kL=1 and kc=NB=9)
Fig. 1
Comparison of the BSM and the analytical solution (kL=1 and kc=NB=9)
Close modal

The question that is addressed in this work is whether this problem can be tackled solving just a few small blocks of size λc and reconstructing the solution between them using a Fourier series. Figure 1 shows, in white, blocks that are actually solved, and in gray the interpolated or reconstructed regions that are used just for displaying purposes. The black squares in such figure represent the grid which is actually solved.

The computational domain is subdivided into NB identical blocks of size λc discretized with NP mesh points each. The block size, which is equal to 2π/NB, is set to contain the short-wavelength λc and the higher spatial harmonics included within each block k=nkc, where n is an integer. On the other hand, the number of uniformly distributed blocks, along with the whole domain, with spacing s=2π/(2M+1), determines the number of spatial harmonics, M, of the block modulation reconstruction and enables the resolution of the long-wavelength harmonics with wavenumber kL.

Every single block is computed independently as a separate instance of the solver, and the blocks play the role of the sampling points in a Fourier transform. Different numerical schemes (e.g., finite-difference, finite volume, spectral methods, etc.) can be used to compute the solution inside each block but boundary conditions have to be provided at the interfaces with the neighboring blocks to account for the long-wavelength variations. In this particular case, the first spatial harmonic (see Fig. 1), with kL=1, and a short-wavelength perturbation of the same size as the block (kc=NB=9) are computed using only three blocks (j=1, 2, and 3 white regions in Fig. 1). The existence of a disparity of length scales in the problem is explicitly exploited in the method to improve the computational efficiency.

The very same idea appears in several turbomachinery problems though. The short wavelength is always associated with airfoil perturbations, and therefore, has the size of the blade pitch. The generation of the long characteristic length can be due to an explicit perturbation like intake distortion or can be more subtle as it is the case when the long scale is associated with flow instabilities such as rotating stall or rotating disk instabilities. The size of these perturbations is typically of the order of the rotor circumference, and the ratio between the long and short characteristic lengths of the order of the number of blades of the full wheel oscillates between 20 and 150.

4 Harmonic Content of the Solution

This section aims to foresee the harmonic content of the solution retained by the method. It is apparent that since only some of the nodes or blocks of the full computational domain are solved, the solution is necessarily filtered out. If each block is discretized with NP points and NB=2M+1 blocks, the total number of harmonics is hence NPNB/2.

Figure 2 shows the harmonic content resolution of the block Fourier transform in a schematic way, where the bars represent the harmonics that the method can capture and are within the resolution range of the method.

Fig. 2
Schematic view of BSM harmonic content
Fig. 2
Schematic view of BSM harmonic content
Close modal

Four families of spatial harmonics are then included in the reconstruction of the full-domain solution using the block spectral method:

  1. The mean value or the 0th harmonic.

  2. N=(NP1)/2 block-embedded harmonics, k=nkc, where n is an integer such that nN. These inner modes are self-contained in the block domain and are not affected by any Fourier transform. They are solved inside each block with the numerical method of choice. The fundamental short-wavelength perturbation is represented and named in Fig. 2 as kc, while the higher harmonics of the solution are represented by nkc.

  3. M external block modulating harmonics with wavenumber k=m, where m is an integer such that mM. These harmonics are reconstructed using 2M+1 uniformly spaced blocks to retain the fundamental long-wavelength harmonics kL, shown as bars with circles in Fig. 2.

  4. 2MN scattered harmonics, k=nkc±m, which are implicitly calculated by the method. These modes are shown in Fig. 2 as side bands surrounding the nkc harmonics.

The method was conceived originally to reproduce just the long- and short-wavelength harmonics, represented by the wavenumbers kL and nkc, respectively. However, the block Fourier reconstruction approach followed by the block spectral method has shown to be better than initially expected, and the scattered modes are natively retained in the method as well. This fact turns the technique into a powerful tool for turbomachinery applications where the scattered modes play a central role.

Turbomachinery configurations exhibit the so-called Tyler–Sofrin modes [14] or spatial harmonics k=mNb+nNv, where Nb represents the number of rotor blades, and Nv refers to the total number of vanes. These are visible in forced response or acoustic simulations increasing the size of computational models consisting of multiple airfoil rows. The fundamental harmonic of the wake rotor with wavenumber kL=Nb interacts with the vane perturbation, with kc=Nv, to generate spinning modes (or Tyler–Sofrin modes). These are nothing else than scattered modes. Depending on the actual blade count, these scattered modes can also exhibit very long-wavelengths. What is relevant here is that the structure of the response of the physical and the numerical models are fully compatible.

5 Numerical Method Outline

The method is illustrated by discretising spatially the LWE, ut+c0ux=0, using a first-order finite-difference method because of its simplicity. However, the idea behind the method can be implemented in any solver with a small numerical stencil regardless of the spatial accuracy of the scheme. The temporal discretization does not interfere with the method at all as long the scheme is explicit or, if implicit, is based at least on marching in a pseudo-time, as is the case of the dual time-step method. The semi-discrete first-order form of the wave equation at the node i and block b can be written in the inner points of the mesh as
ui(b)t=c0[ui+1(b)ui1(b)2ΔxμΔx(ui+1(b)2ui(b)+ui1(b)Δx2)]
(6)
where central differences have been used to discretise ux, and numerical diffusion of the form μc0Δxuxx has been added to give a certain upwinding character to the discretization. The baseline code is second-order accurate in space [15] but follows a similar structure. For uniformly distributed grid points, the mesh spacing is Δx=xixi1, and μ is a dimensionless constant. Every block is computed in the solver separately, and the index i is local to the block.
Figure 3 introduces the nomenclature used in Eq. (6) where ui(b) represents the unknown in the ith node of the block b. The index i is local to the block, and nodes with identical i index are physically located at the same relative position within the block. These nodes are referred to as homologous nodes in the blocks. However, the block index b is global and spans from b=1,,NB. As it has been mentioned before, the idea of the method is not to discretize all the blocks but just those corresponding to the indexes:
b=1+(j1)NB/(2M+1)
(7)
where j spans from j=1,,2M+1. Therefore, the block index j refers to the local and consecutive list of the blocks which are actually solved. The spatial discretization described in Eq. (6) can be directly applied to all the inner points of the block. If the block is discretized using NP points, the block length is L=NPΔx. The node index i spans from i=1,,NP and hence the scheme in the block boundaries is
i=1:u1(b)t=c0Δx[12(u2(b)u0(b))μ(u2(b)2u1(b)+u0(b))]i=NP:uNP(b)t=c0Δx[12(uNP+1(b)uNP1(b))μ(uNP+1(b)2uNP(b)+uNP1(b))]
(8)
where the nodes u0(b) and uNP+1(b) do not form part of the block domain. If the block b were solved with periodic boundary conditions, the problem would be closed by setting u0(b)=uNP1(b) and uNP+1(b)=u2(b). In this case, the solution of the block is self-contained and the numerical method does not need to resort to out-of-block information.
Fig. 3
Sketch of the node and block index notation used in the one-dimensional numerical scheme
Fig. 3
Sketch of the node and block index notation used in the one-dimensional numerical scheme
Close modal
Since the actual problem does not involve periodic blocks, boundary conditions have to be provided. If the information at all the nodes of the full-domain were available, these boundary conditions would be
u0(b)=uNP1(b1),uNP+1(b)=u2(b+1)
(9)
However, the nodes uNP1(b1) and u2(b+1) belong to the neighboring blocks (gray region in Fig. 3) and their value is unknown since these blocks are not contained in the discretization. In the present method, the solution at the nodes uNP1(b1) and u2(b+1) is reconstructed using the existing information in the 2M+1 existing blocks. Among the many different forms of reconstructing the signal in the full-domain, the most natural, and efficient in this case, is the use of discrete Fourier transforms. The Fourier coefficients of the inner homologous points at both periodic boundaries are computed as
u^NP1m(t)=12M+1j=12M+1uNP1(j)(t)eI2πmj/(2M+1)u^2m(t)=12M+1j=12M+1u2(j)(t)eI2πmj/(2M+1)
(10)
and used as interpolant functions to reconstruct the solution in the block boundaries. Figure 4 sketches the idea of the block-wise spectral reconstruction strategy. The Fourier coefficients in Eq. (10) of the solution at the inner nodes of both periodic boundaries are calculated using the uNP1(j) and u2(j) homologous boundary nodes in all the existing 2M+1 blocks (represented as filled red and blue squares in Fig. 4). The solution at the nodes u0(b) and uNP+1(b), located at the neighboring blocks b1 and b+1, respectively, are computed using the Fourier coefficients as
u0(b)=m=1Mu^NP1m(t)eI2πm(b1)/NBuNP+1(b)=m=1Mu^2m(t)eI2πm(b+1)/NB
(11)
Taking this into account, the solution at the periodic boundaries of each block b can be expressed formally as a function of the inner solution of all the homologous points of the sampling blocks
u0(b)=f(uNP11,,uNP11+(j1)K,,uNP11+(2M)K)uNP+1(b)=f(u21,,u21+(j1)K,,u21+(2M)K)
(12)
with spacing K=NB/(2M+1), where M is the number of harmonics in the Fourier transformation. The discrete Fourier transform used to reconstruct the solution in the block boundaries in Eq. (12) is a linear combination of homologous points. Once the fluxes are computed, and the variables in the blocks are updated, the solutions uNP+1(b) and u0(b) can be recomputed, and the whole process is repeated.
Fig. 4
Schematic of the block-wise spectral reconstruction of the first harmonic (kL=1) in a domain consisting of six blocks, NB=6. Solid diamonds: data used to reconstruct the solution on the boundaries. Empty diamonds: interpolated data outside of the computational domain. Circles: computational node.
Fig. 4
Schematic of the block-wise spectral reconstruction of the first harmonic (kL=1) in a domain consisting of six blocks, NB=6. Solid diamonds: data used to reconstruct the solution on the boundaries. Empty diamonds: interpolated data outside of the computational domain. Circles: computational node.
Close modal

The local point-wise scheme within the uniformly spaced blocks to accommodate the in-block high-harmonic content with wavenumbers k=nkc enables the reconstruction of the solution without any loss of information. In other words, the harmonics contained within the blocks are always well represented within the accuracy of the method employed to solve the blocks.

If M spatial harmonics are retained in the block Fourier transform, then just 2M+1 blocks out of the total number of NB blocks need to be computed. Since only a neighboring node at each periodic boundary of the jth block is needed to construct the numerical scheme due to the compact stencil of the spatial discretization, only the Fourier coefficients of the adjacent inner nodes to the original periodic boundaries need to be evaluated. Since the number of points involved in the Fourier transform is proportional to the number of spatial harmonics, the method can eventually become inefficient if the number of retained harmonics is high. However, the facts that just a few harmonics are usually required in turbomachinery problems, and that the Fourier transforms need to be performed only in a small subset of the mesh points, make the penalty associated with the block-wise spectral reconstruction at the boundaries negligible.

The computational efficiency of the method is then proportional to the ratio between the total number of blocks and the number of the retained spatial harmonics in the simulations. The ideal speed-up factor, R, with respect to the full-domain solution is then
R=NB2M+1
(13)
The method is especially efficient when the number of blocks or blade passages is high. This number oscillates between 20 and 150 in realistic cases.

6 Spectral View of Block Reconstruction

The method can be more easily understood if it is regarded as a sequence of two aligned Fourier expansions in the same direction. Let us consider that the governing equations within each block were solved using a spectral method with periodic boundary conditions. In this case, the solution in the block b could be expanded in Fourier series as
u(b)=k=NP/2NP/2u^k(b)eIkx
(14)
where x is the block local coordinate in the span x[0,2π/NB]. If the solution is block periodic and all the blocks are identical, in the global coordinate system, the solution can be expressed as
u=k=NP/2NP/2u^k(b)eI(kNB)x
(15)
If the solution is not block periodic and varies from block to block, the complex Fourier coefficients, u^k(b), can be block-wise Fourier transformed to obtain
u=k=NP/2NP/2m=MMu^k,meI(kNB+m)x
(16)
which is a double Fourier expansion in the same direction, x. This procedure is analogous to the one that would be conducted in a two-dimensional spectral solver with periodic boundary conditions where the variables are first Fourier transformed in the x-direction and then transformed again in the y-direction. Here, the solution is transformed twice in the same direction.

In this context, the wavenumbers kNB represent the block inner harmonics while the wavenumbers m account for the out-of-block harmonics, i.e., harmonics with long wavelengths that are not self-contained inside the block and modulate the block solution. There are NP/2 inner harmonics and M outer harmonics. The simulation contains the spatial harmonics kNB+m, which are a particular subset of all the potentially existent modes. Among all the existing harmonics, MNP modes are obtained as the combination of inner and outer wavenumbers. These are referred to here as scattered modes since they are obtained physically from the interaction between the inner and outer scales.

It is important to highlight that since the first Fourier transform is real and the second is complex, u^k,m=u^k,m* but u^k,mu^k,m*. This means that both positive and negative m indexes contribute to the simulation and are independent.

In the present method, the “inner” Fourier transform does not exist since the method does not have a spectral discretization. In practice, the inner block is solved using a second-order finite-difference or finite volume method (illustrated in this section using a first-order discretization for the sake of simplicity). However, the outer Fourier reconstruction is exact and can exactly replicate the inexact solution that would be obtained if the whole computational domain had been simulated.

Figures 5(a)5(c) illustrate the spectrum content as a function of the number of reconstructed harmonics, M, used by the Fourier transformation for M=1, 2, and 3, respectively. Using 2M+1 uniformly spaced blocks, M harmonics can be reconstructed. The figure represents these harmonics as kL, close to the origin, and their number increases with M.

Fig. 5
Sketch of the spectrum content of the method for different values of the number of spatial harmonics, M, retained in the Fourier reconstruction; (a) M = 1, (b) M = 2, and (c) M = 3
Fig. 5
Sketch of the spectrum content of the method for different values of the number of spatial harmonics, M, retained in the Fourier reconstruction; (a) M = 1, (b) M = 2, and (c) M = 3
Close modal

Embedded harmonics, which are self-contained inside each block, are noted as nkc. Two of them, kc and 2kc, are represented in Fig. 5. Also, the scattered harmonics appear in the spectrum as sidebands of the inner harmonics, nkc±m, since the method is designed in such a way that kc1. By increasing the number of outer harmonics, the sidebands about the harmonics nkc widen, as shown in Figs. 5(b) and 5(c).

7 Numerical Results for the Linear Wave Equation

This section presents the application of the block spectral method to the resolution of the one-dimensional wave equation with periodic boundary conditions. The objective is to validate and illustrate the conclusions outlined in the previous section. The ability of the scheme to solve different initial value problems with varying numbers of initial spatial harmonics, reconstructed spatial harmonics, M, and blocks NB is displayed through examples.

The method solution is compared to the full-domain solution in various cases to verify the results. Snapshots of the solution are taken after several “flow-through” times. The computational mesh consists of NP=100 points per block. Numerical parameters for the resolution of the LWE using the first-order scheme are c0Δt/Δx=0.1 and μ=0.5.

The numerical scheme used to describe the method is simple and rudimentary but very high resolutions have been employed here to avoid contaminating the solution with issues related to the low quality of the underlying scheme, which has nothing to do with the main idea of the method presented in this work.

7.1 Single Outer Harmonic Case.

First, a validation case is carried out to verify the correct reconstruction of a spatial harmonic kL spanning outside of the block boundaries. Let us consider the initial condition
u(x,t=0)=sin(kLx)
(17)
with kL=2. The BSM method is compared to the full-domain solution (full-annulus solution in the turbomachinery jargon) in the whole domain, being NB=20. The problem is solved using the BSM with M=2, which means that five uniformly spaced computational blocks are used to calculate the second harmonic of the initial value problem in Eq. (17).

Figure 6(a) represents the full-domain and the BSM solutions. The gray region displays the solution in the full computational domain, while the white areas plot the solution of the BSM at each jth block using the block-wise spectral reconstruction previously described. As shown in Fig. 6(a), both the full-domain and the BSM solutions are identical. The BSM reproduces the full-domain solution with the corresponding reduction in grid points and computational time since the method only solves five out of the NB=20 blocks to reconstruct the whole domain solution, which corresponds to a speed-up factor of 4.

Fig. 6
Comparison of the BSM and full-domain solutions. (a) Initial solution with kL=2, M=2 (five computational blocks in white), and NB=20. (b) Initial solution containing the harmonics kL=2 and kS=15, M=2, and total number of blocks NB=15.
Fig. 6
Comparison of the BSM and full-domain solutions. (a) Initial solution with kL=2, M=2 (five computational blocks in white), and NB=20. (b) Initial solution containing the harmonics kL=2 and kS=15, M=2, and total number of blocks NB=15.
Close modal

7.2 Outer Harmonic Superimposed With Embedded Harmonics Case.

Once it has been verified that the method can reconstruct solutions containing wavelengths longer than the block size, a second perturbation with a wavelength contained within the block size is superimposed to the first. The initial conditions are
u(x,t=0)=sin(kLx)+αsin(kSx)
(18)
being α=0.4, kL=2, and kS=15. The second term of the initial condition has been enforced to have a wavelength equal to the block size since the number of blocks in which the domain is subdivided is NB=15. The full-domain and the BSM solutions with two harmonics (M=2) are compared in Fig. 6(b). Both solutions show the same sinusoidal variation as the imposed initial value with spatial variations inside the blocks as a consequence of the presence of the embedded harmonic kc. The BSM solution can reproduce the two different spatial scales. The fundamental large-scale harmonic, kL, is reproduced with the two spatial harmonics retained in the Fourier transform. The second, related to the short-wavelength pattern, is solved by the inner block discretization since the harmonic kc is spatially periodic in each block. This case corroborates the idea that the method can reproduce simultaneously wavelengths contained in the block and wavelengths longer than the block size.
The aim now is to verify that the proposed BSM can reproduce any number of higher periodic harmonics inside each block domain as long as the internal discretization of the block supports them. The LWE is solved with the following initial condition:
u(x,t=0)=u0(x)=αsin(kLx)+b=1NBea[x(2b1)(π/NB)]2
(19)
which represents a sinusoidal wave modulating a compact solution within the blocks (second term in Eq. (19)) if the parameters are properly chosen. The amplitude of the first spatial harmonic (kL=1) has been set to α=0.2. The second term in Eq. (19) represents a summation of Gaussian functions. The value of a has been set to a=10 to ensure that the tails of the Gaussian do not extend beyond the limits of the blocks.

It could be argued that this case is irrelevant since if the problem is linear, the results could be presented in a harmonic-by-harmonic basis. However, it was thought that displaying results with a more “physical” initial condition could help to grasp the engineering applications of the method.

Figure 7 shows the comparison between the BSM (M=1) and full-domain solutions in a case where NB=9. The white regions in Fig. 7(a) represent the computational domain for the BSM. In these zones, the solution of the BSM is compared to that of the full-domain solution. Both solutions are shown to be identical. Figure 7(b) shows the harmonic content of both BSM and full-domain solutions. The block-to-block deviation introduced by the first spatial harmonic of the initial solution is seen in Fig. 7(a). Figure 7(b) displays the harmonic content of the initial solution and the solution after several wave passings. The initial condition (Eq. (19)) contains high-order harmonics embedded inside each block. It is interesting to verify if the method can capture all these block-embedded harmonics nkc. These correspond, in Fig. 7(b), with the Fourier indexes 9, 18, 27, etc. Table 1 compares the harmonic content of the solutions obtained using the full-annulus domain and the BSM with M=1 for this case. The results indicate the accuracy of this approach. It is concluded that the block-wise spectral reconstruction inherently reproduces all the block-embedded spatial harmonics nkc for any integer n(NP1)/2.

Fig. 7
Comparison of BSM (M=1, three computational blocks) and full-domain solutions (a) and harmonic Fourier content (b) with an initial solution described in Eq. (19), where kL=1, NB=9, α=0.2, and a=10
Fig. 7
Comparison of BSM (M=1, three computational blocks) and full-domain solutions (a) and harmonic Fourier content (b) with an initial solution described in Eq. (19), where kL=1, NB=9, α=0.2, and a=10
Close modal
Table 1

Harmonic content of the solutions obtained with the full-annulus domain and the BSM (M=1) for the case kL=1, NB=9, α=0.2, and a=10 displayed in Fig. 7 

FULLPSM(M=1)|FULLBSM|
kL0.1992120.1992135.02×106
kc0.4338710.4338694.61×106
2kc0.2154190.2154219.28×106
FULLPSM(M=1)|FULLBSM|
kL0.1992120.1992135.02×106
kc0.4338710.4338694.61×106
2kc0.2154190.2154219.28×106

It is apparent that the numerical accuracy of the computed high-order harmonics will be low, especially if a first-order method is used. However, the BSM is transparent in this regard, and the accuracy of the original method implemented in the blocks is preserved.

7.3 Nonmatching Block Case.

For an arbitrary number and size of physical blocks, the use of uniformly spaced discrete Fourier transforms may lead to the use of computational blocks which do not fit the physical blocks. This fact may contribute to increasing the error with respect to the cases presented up to now where the number of blocks was chosen in such a way that the blocks used to reconstruct the signal corresponded to physical blocks. This degree-of-freedom does not exist in practice where the number of blocks is given by the number of blades.

To shed some light about the case when the physical position of some of the 2M+1 block samples in the BSM does not coincide with the real position of the NB blocks in the full-domain, a case with three block samples (M=1) to retain the first harmonic and a total number of blocks NB=10 is studied. In this case, the collocation of harmonic kL=1 with a single-mode (M=1) requires using “unphysical” non-fitting blocks with the original mesh. Figure 8(a) shows this case’s BSM and full-domain solutions. It can be seen that in the BSM, the physical positions of the second (j=2) and third (j=3) blocks do not correspond with any of the blocks of the full computational domain. The position of these virtual blocks is imposed by the need of using an equispaced Fourier transform with a minimum number of blocks (three in this case). The solution of the BSM appears displaced with respect to one of the full-domain solution since the blocks are in a different position. Nevertheless, the local amplitude of the first harmonic is retained. However, the matching is excellent when the solution is appropriately reconstructed and compared to the full-domain solution.

Fig. 8
Comparison of BSM (M=1, three computational blocks) and full-domain solutions (a), BSM-reconstructed solution (b) and harmonic Fourier content (c) with the initial solution described in Eq. (19), where kL=1, NB=10, α=0.2, and a=10
Fig. 8
Comparison of BSM (M=1, three computational blocks) and full-domain solutions (a), BSM-reconstructed solution (b) and harmonic Fourier content (c) with the initial solution described in Eq. (19), where kL=1, NB=10, α=0.2, and a=10
Close modal

Once the solution in the uniformly spaced blocks is calculated, the full-domain solution is reconstructed by Fourier transforming all the homologous nodes of the blocks. In the example presented here, only the first harmonic can be obtained. The discrete Fourier series of each mesh node is then used as an interpolant function to derive the solution in the physical nodes. This Fourier series can be seen as an expansion of the first virtual block (j=1), always chosen to match the first physical block.

Figure 8(b) compares the reconstructed BSM solution to the full-domain solution using the three uniformly spaced blocks to reconstruct the complete solution. Even though the blocks are not placed in the same physical position, both solutions are identical, as shown in Fig. 8(b) in the physical plane and Fig. 8(c) in the Fourier domain. Table 2 compares the harmonic content of the solutions obtained using the full-annulus domain and the BSM with M=1. The same level of accuracy as in the case where the collocation blocks coincide with the physical ones has been obtained. It is also concluded that the BSM can calculate the exact solution of the domain-embedded harmonics even if the virtual blocks do not match the physical blocks.

Table 2

Harmonic content of the solutions obtained with the full-annulus domain and the BSM (M=1) for the case kL=1, NB=10, α=0.2, and a=10 displayed in Fig. 8 

FULLPSM(M=1)|FULLPSM|
kL0.1981870.1981865.05×106
kc0.4341220.4341196.91×106
2kc0.2142590.2142579.33×106
FULLPSM(M=1)|FULLPSM|
kL0.1981870.1981865.05×106
kc0.4341220.4341196.91×106
2kc0.2142590.2142579.33×106

7.4 “Scattered” Mode Case.

A case to verify the accuracy of the scattered or sideband harmonics nkc±m is presented next. The aim is to verify the accuracy of the block-wise spectral reconstruction for mixed modes obtained by a combination of Fourier-reconstructed and block-resolved parts of the solution. These modes correspond to the sidebands sketched in Fig. 2.

Since the LWE is non-dispersive and cannot scatter modes, these are explicitly included in the initial condition. In this particular case, the initial state consists of two sines with different wavelengths, as in the previous case (see Eq. (18)), being α=0.4 and kL=2, where the domain is subdivided into NB=15 blocks. However, in this case, the second term of the initial value in Eq. (18) corresponds to the spatial harmonic kS=13. This mode is neither directly included in the Fourier reconstructions nor fitting inside a single block. Therefore, the ability of the method to deal with this type of problem deserves to be pursued.

The full-domain and BSM solutions with two harmonics (M=2) are compared in Fig. 9. It can be appreciated that both solutions are identical. This is consistent with the basic idea that the method can deal not only with the block-embedded modes, in this case, kc=NB=15 and all its multiples, but also their sidebands. In this case, since M=2, the mode kS=kcm, with m=2, is retained in the solution. The modes 13, 14, 15, 16, and 17 would be reproduced accurately if they were present in the initial solution. This case corroborates the idea that the method can reproduce as well scatter or sideband harmonics.

Fig. 9
Comparison of the BSM (M=2, five computational blocks) and the full-domain solutions with the initial condition described in Eq. (18) where kL=2, kS=13, NB=15, and α=0.4
Fig. 9
Comparison of the BSM (M=2, five computational blocks) and the full-domain solutions with the initial condition described in Eq. (18) where kL=2, kS=13, NB=15, and α=0.4
Close modal

7.5 Out-of-Spectrum Modes.

The limitations of the proposed method are explored in the following lines. Let us simulate now an initial value problem consisting of two different sinusoidal perturbations with different wavelengths again. The form of the perturbation is exactly the same as the one described in Eq. (18) being, in this case, α=0.4 and kL=2. The domain is subdivided into NB=15 blocks. The wavenumber of the second term in Eq. (18) is set to kS=5. Figure 10 shows that the BSM solution with two harmonics (M=2) is quantitatively different than that of the full-domain solution. The BSM cannot reproduce the mode k=5 since it is not included in the list of potential scattered modes and the solution of the BSM has little to do with that of the full-domain in this case. Thus, the mode kS=5 is aliased in other low- and high-order harmonics, as shown in the harmonic content of the BSM solution in Fig. 10(b).

Fig. 10
Comparison of the BSM (M=2) and the full-domain solutions with the initial solution described in Eq. (18) where kL=2, kS=5, NB=kc=15, and α=0.4. (a) Snapshot of the solution. (b) Harmonic content as a function of the wavenumber.
Fig. 10
Comparison of the BSM (M=2) and the full-domain solutions with the initial solution described in Eq. (18) where kL=2, kS=5, NB=kc=15, and α=0.4. (a) Snapshot of the solution. (b) Harmonic content as a function of the wavenumber.
Close modal

The cases displayed in Figs. 9 and 10 represent the same numerical setting but include different initial conditions: the modes kS=13 and kS=5, respectively. The mode kS=13 is higher than the kS=5, but it is represented much better because the BSM’s harmonic content is not uniform (see Fig. 2).

8 Solution of the Wave Equation With Varying Propagation Speeds

The homogeneous linear wave equation, apart from the domain size, lacks any built-in characteristic length. The method has been artificially illustrated, setting different initial conditions providing different characteristic lengths. Block sizes have been arbitrarily set to demonstrate the technique since there is no natural way to choose them. This is not the case in turbomachinery problems where the blade pitch determines the block size, and this is no longer a degree-of-freedom. In the absence of incoming perturbations, rotor blades give rise to a nonuniform steady flowfield whose fundamental wavelength is the pitch of the blades. When spatially periodic unsteady perturbations are injected into the domain, they interact with the nonuniformity, creating new spatial harmonics that were not included in the original perturbation.

The wave equation with variable propagation speed c(x)
ut+c(x)ux=0
(20)
is an interesting vehicle for testing the method since it enables the possibility of including an intrinsic characteristic length in the spatial variation of the propagation speed, precisely as in the case of interest. Moreover, the variable propagation speed, c(x), gives rise to scattering, i.e., pure sinusoidal waves are no longer a solution to the problem. Incoming or existing waves in the domain are dispersed by c(x), creating waves with wavelengths different from the original one. The variable speed c(x) plays the role of a nonuniform field created supposedly by a row of blades.
Among the many potential choices for c(x) in this work, the following expression has been chosen:
c(x)=c01+ϵcos(kcx)
(21)
where ϵ controls the level of the nonuniformity and λc=2π/kc is the wavelength of the variation of the propagation velocity. The reason for selecting this form of c(x) is threefold. First, it eases the analytical solution of the problem, second, the variation is spatially periodic, and third, enables the analytical determination of the Fourier content of the solution (see  Appendix).
Finally, the problem is solved with the following initial condition:
u(x,t=0)=up(x,t=0)
(22)
where up(x,t) is a solution to the problem. The idea behind the use of a solution of the problem as an initial condition is to avoid a long transient leading to such solution. The solution of Eq. (20) with the propagation speed given by Eq. (21) is readily obtained by performing the following change of variables dη=(1+ϵcoskcx)dx to reduce the problem to that of solving the LWE with constant velocity in the transformed space, i.e., ut+c0uη=0.
The solution can be obtained after substituting η=x+(ϵ/kc)sin(kcx) in the solution of the transformed wave equation. An exact solution of the wave equation with the propagation velocity defined in Eq. (21) is
up(x,t)=sin(kL(xc0t)+βsin(kcx))
(23)
being β=ϵkL/kc. This analytical solution is used to compare the results obtained using the BSM and the full-domain computation. The harmonic content of the solution is derived in Appendix, where it can be seen how the modulation of the harmonic kL in Eq. (23) gives rise to sidebands around the harmonic kc, as shown in Eq. (A4) and Fig. 17 in Appendix for different values of the modulation index β.

Figure 11 shows the solution of Eq. (23) for kL=1, kc=6, and β=0.5. The solution can be expressed as a function of Bessel functions of the first kind Jn(β) (see Eq. (A4) in Appendix). The modulation of the mean flow velocity c(x) in Eq. (21) gives rise to the formation of scattered harmonics. It is essential to highlight that the energy of the system is constant and equal to that of an initial condition containing just the fundamental harmonic, i.e., u0=sin(kLx). This can be demonstrated by making use of the properties of the Bessel functions of first kind described in Appendix.

Fig. 11
Snapshot of the analytical solution for the varying speed linear wave equation in Eq. (23)(a) and harmonic content as a function of the Fourier index (b) with kL=1, kc=6, and β=0.5
Fig. 11
Snapshot of the analytical solution for the varying speed linear wave equation in Eq. (23)(a) and harmonic content as a function of the Fourier index (b) with kL=1, kc=6, and β=0.5
Close modal

The main objective of this section is to elucidate if the method can accurately capture these scattered harmonics as it has been previously claimed. The block spectral solution with one spatial harmonic (three blocks) is compared against the full-domain and the analytical solutions. In this case, kc=NB to force the blocks to contain the periodic spatial variation of the propagation speed. Figure 12 compares the solutions obtained by the block spectral method and the full-domain simulation to the analytical solution for M=1 and two different wavenumbers of the spatial perturbation, kc=4and6, with a fixed value of ϵ=0.8 (see Eq. (21)). The excellent agreement between all of them allows to conclude that the block spectral method can capture the interaction between the nonuniform propagation speed with wavelength 2π/kc and the harmonic kL injected in the initial condition.

Fig. 12
Comparison of the BSM (M=1, three computational blocks), the full-domain, and the analytical solutions (sub-figures (a) and (b)) and their harmonic contents (sub-figures (c) and (d)) for the case ϵ=0.8 and kL=1, for two different wavenumbers of the propagation speed, kc=4 and kc=6. The wavelengths of the propagation speed and the block size are identical, NB=kc.
Fig. 12
Comparison of the BSM (M=1, three computational blocks), the full-domain, and the analytical solutions (sub-figures (a) and (b)) and their harmonic contents (sub-figures (c) and (d)) for the case ϵ=0.8 and kL=1, for two different wavenumbers of the propagation speed, kc=4 and kc=6. The wavelengths of the propagation speed and the block size are identical, NB=kc.
Close modal

Figures 12(c) and 12(d) show the harmonic content of the three solutions displayed in Figs. 12(a) and 12(b) for kc=4 and 6, respectively. Both cases exhibit sidebands of the modulating harmonic kc, and the spectra contains a carrier harmonic of amplitude close to 1 for the wavenumber k=1, due to its presence in the initial solution. These results can be directly compared to the ones obtained in Appendix for the corresponding values of the parameter β. In the cases studied in this work, the modulation index is very small, β1, since we are targetting the limit kLkc, and then the amplitude of the carrier harmonic kL is J0(β), which is approximately unity. On the other hand, since we have a modulated wave, two more harmonics appear around the modulated harmonic kc, in the form of sidebands, as shown in Figs. 12(c) and 12(d) for kc=4 and 6, respectively. These sidebands in the spectra correspond to the scattered harmonics ks=kc±m whose amplitude in the limit β1 is J1(β)β/2. In principle, an infinite number of harmonics is generated by the interaction of the initial condition with the nonuniform mean flow perturbation c(x). However, since the modulation index is small, the higher-order sidebands, ks=nkc±mkL, are very small because their amplitudes decay as O(βn) if β0. This is seen, for instance, in Fig. 12(c) for kc=4, where the harmonics 7 and 9 are hardly visible in the spectrum. Although the amplitude of the sidebands decays very quickly, the block spectral method can retain all the scattered harmonics. This applies as well to the cases where βO(1), though such cases have not been illustrated here.

Table 3 compares the amplitude of the fundamental, kL, and the lower scattered harmonics, kc±kL, of the solutions obtained using the full-annulus domain and the BSM with M=1 for the case ϵ=0.8, kL=1, and kc=6 whose solution is displayed in Figs. 12(b) and 12(d). It is concluded that the accuracy of the solutions is very high, even for the scattered harmonics and that the BSM is able to obtain a good agreement against full-annulus and analytical solutions.

Table 3

Error between the solutions of the full-annulus simulation and the BSM (M=1) in the first scattered harmonic for the case ϵ=0.8, kL=1, and kc=6

FULLPSM(M=1)|FULLBSM|
kL0.9903020.9903042.02×106
kckL0.06616790.06616824.53×106
kc+kL0.06616960.06616983.02×106
FULLPSM(M=1)|FULLBSM|
kL0.9903020.9903042.02×106
kckL0.06616790.06616824.53×106
kc+kL0.06616960.06616983.02×106

9 Generalization to Turbomachinery Applications

The block spectral method has been implemented in an unstructured edge-based second-order finite volume compressible URANS solver for turbomachinery applications [15,16]. The details of the implementation in the context of the block spectral method are discussed in Ref. [12] and applied to investigate the nonlinear stability of a fan [13]. However, for the sake of completeness, the method is briefly described next.

9.1 Baseline Flow Solver.

The Navier–Stokes equations in conservative form for an arbitrary control volume can be written in compact form as
ddtΩUdΩ+ΣF(U).dA=0
(24)
where U is the vector of conservative variables, F is the sum of the inviscid and viscous fluxes, Ω is the volume of the control volume, Σ is its boundary, and dA is the differential area pointing outward to the boundary.

The baseline solver, Mu2s2T, discretizes the spatial domain using hybrid unstructured grids and can contain cells with arbitrary number of faces. The governing equations are discretized following a finite volume approach based on the dual mesh where the solution vector stores the variables at the vertices of the mesh. The control volume associated with a node is formed by connecting the median dual of the cells surrounding it, using an edge-based data structure. The spatial discretization is second-order accurate, and the numerical dissipation terms’ matrix form is used to damp high-frequency modes.

Once the spatial domain is discretized, the semi-discrete version of nonlinear equation (24) for the internal node i can be written in compact form as
dUidt=R(U)
(25)
and the resulting equations are integrated in time using a second-order backward implicit scheme
3Uin+14Uin+Uin12Δt=R(Un+1)
(26)
where n is the time level. The resulting system of nonlinear equations is solved using a dual-time-step, τ, adding the pseudo-derivative Uτ [17], and preconditioning the system to speed up the convergence to the pseudo-steady-state Un+1
Pi1dUidτ=R*(U)=R(U)3Ui4Uin+Uin12Δt
(27)

The solution of this equation is sped-up using steady-state acceleration techniques, marching in τ instead of t. The preconditioning matrix Pi1 can include block-Jacobi, multi-grid, and any other mechanism to accelerate the convergence to the pseudo-steady-state Un+1. The acceleration mechanisms used to remove the time-step limitations of the explicit scheme can be regarded as an inner iteration of the method and are transparent to the block spectral method.

The vectors of unknowns, U, and residuals, R, in a full-annulus simulation can be split in a block-wise manner as
U=(Ui1Ui2...UiNB),R=(Ri1Ri2...RiNB)
(28)
where the column vector U(j) contains just the unknowns of the jth block and NB is the number of blocks. The node i represents homologous nodes in different blocks or passages. In turbomachinery problems, the computational domain of one blade row is typically decomposed into N passages, hence the name PSM.

9.2 Three-Dimensional Block Spectral Method.

The homologous nodes of the different passages, Ui(j), can be Fourier transformed in the circumferential direction to obtain their Fourier coefficients
U^im(x,t)=12M+1j=12M+1Ui(j)(x,r,θj,t)eImθj
(29)
If all the passages of the wheel are used, then (NB1)/2 different complex Fourier coefficients, U^im(x,t), are obtained at each node i and the exact Fourier representation is obtained. Instead, a truncated Fourier series with MNb/2 harmonics is computed using just a few blocks to reduce the computational time.

As it has been mentioned before, the rationale is that the wavelength of the circumferential nonuniformities are long, and therefore, a few harmonics suffice to represent the flow field accurately. The idea of the method is sketched in Fig. 13, where it can be seen that the first harmonic of the circumferential variations of the flowfield is reconstructed using only three passages (in dark gray). The solution in the rest of the passages (in light gray) is interpolated using a Fourier reconstruction containing solely the information of a few blocks. The method is aimed at reproducing long-wavelength patterns of the size of the whole circumference and short-wavelength disturbances of the order of the blade pitch, which are absorbed within the block-passage mesh, as shown in Fig. 13 by dark gray passages. The long-wavelength variations are explicitly retained using the Fourier coefficients obtained by sampling uniformly spaced blocks along the circumference. It should be recalled that the technique does not require the wavelength of the explicitly retained modes to be long compared to the airfoil pitch.

Fig. 13
Scheme of whole annulus sinusoidal variation (left) and its approximation (right) with three reduced passages samples (dark gray) uniformly spaced and nonmatching with the physical position of the airfoils
Fig. 13
Scheme of whole annulus sinusoidal variation (left) and its approximation (right) with three reduced passages samples (dark gray) uniformly spaced and nonmatching with the physical position of the airfoils
Close modal

The computation of the residual Ri typically involves two layers of neighboring nodes due to the fourth-order artificial viscosity terms which require at least the neighbors of the neighbors in its evaluation. If the sensors for the limiters are taken into account, the stencil of the spatial discretization is even larger, and more neighboring nodes are required. Nevertheless, the numerical scheme is local and the inner nodes that are located at a certain topological distance from the neighboring block can be computed disregarding the block spectral reconstruction.

A phantom mesh at the periodic patches is used to reconstruct the flow variables with the corresponding phase-shift in the nonexistent neighboring block. This halo of cells is obtained duplicating the points and the connectivity in the vicinity of the periodic boundaries and rotating them to locate new elements of the grid on the opposite boundary (see Fig. 14). Once the inner part of the halo cells is computed (red and blue regions in Fig. 14), the Fourier coefficients for these points are obtained. A control volume located in the passage periodic boundary contains inner nodes (red region in the control volume on the left of Fig. 14) and nodes located in neighboring blocks which are reconstructed using the previously computed Fourier coefficients. Once the conservative variables in these nodes are updated, the corresponding correction to the gradients and fluxes at each time-step is performed.

Fig. 14
Dual mesh in the passage block boundary (left) and the halo cells (dashed lines) needed to reconstruct the numerical scheme in the boundary using Fourier series
Fig. 14
Dual mesh in the passage block boundary (left) and the halo cells (dashed lines) needed to reconstruct the numerical scheme in the boundary using Fourier series
Close modal

The baseline numerical scheme works in a point-wise manner and just a small fraction of the total number of nodes needs to be Fourier transformed and, therefore, the computational penalty associated with the specifics of the block spectral method is very small.

9.3 Method Verification in Turbomachinery Applications.

This section aims to verify if the block spectral method can capture the Tyler–Sofrin modes in a realistic configuration. The basic validation of the method for intake-distortion-induced compressor instability has been presented in Refs. [12,13]. The NASA rotor 67, under inlet distortion, has been chosen as the test vehicle. The fan operating at 100% of the design speed and close to the peak efficiency operating point (point A (□) in Ref. [12]) has been subject to a steady sinusoidal variation of total pressure at the inlet of the form P0=P¯0+ΔP0sinkLθ, being kL=1 the wavenumber of the harmonic perturbation in the azimuthal direction. The amplitude of the inlet stagnation pressure variation is ΔP0=5kPa which corresponds to about 5% of the mean inlet total pressure, P¯0. The inlet flow nonuniformity gives rise to a strong forcing on the fan blades. Inlet distortion can also severely compromise fan stability [1820]. However, this operating point under this forcing is stable, and the flow unsteadiness reduces to a periodic forcing in the rotating frame of reference [13].

The semi-unstructured mesh [21] used to discretize the domain is shown in Fig. 15(b). First, a hybrid grid in the blade-to-blade plane with 118 points around the airfoil and 24 quad layers in the boundary layer has been constructed. Next, 151 layers in the span-wise direction are extruded. The semi-unstructured grid of the single-passage domain consists of about 2.3×106 nodes per passage. The size of the full-annulus mesh is approximately 50×106 nodes. The Reynolds number of the simulations is very high (Re106), and the standard kω turbulence model including wall functions was used to alleviate mesh requirements in the boundary layers. More details can be found in Refs. [12,13]. The block spectral method is implemented using the blade passage as the fundamental unit. In this particular case, just the first circumferential mode is retained (M=1), and therefore, just three reduced passages are employed to Fourier transform the flow signal (see Fig. 15(b)).

Fig. 15
NASA rotor 67 full computational domain with three reduced passages (M=1) (marked regions) (a) and semi-unstructured grid view in the meridional plane (b)
Fig. 15
NASA rotor 67 full computational domain with three reduced passages (M=1) (marked regions) (a) and semi-unstructured grid view in the meridional plane (b)
Close modal

Unsteady simulations have been conducted using 400 time-steps per rotor revolution (the NASA rotor 67 has 22 blades). The temporal resolution was studied in a separate work, and the discussion can be found in Ref. [12]. The full-annulus simulations are compared to the single harmonic PSM solution.

The main objective of this work is to verify if the PSM can capture the whole harmonic content of the intake distortion problem; the block-embedded harmonics created by the 22 rotor blades (kc=Nb=22), the fundamental long-scale harmonic, kL=1, imposed by the distorted inflow and the inlet, and the scattered harmonics, k=nkc±m=21 and 23, produced by the interaction between the inlet long wavelength perturbation and the blades.

Figure 16(a) compares the total pressure signal in a numerical probe located near the outer casing, in the relative frame of reference, 0.5 axial chords downstream of the fan trailing-edge (see Fig. 15(b)) obtained by a full-annulus simulation, and the PSM with a single harmonic. The initial solution corresponds to the steady flow in the relative frame of reference obtained for a clean inlet. It can be seen that after approximately four rotor revolutions, the solution becomes periodic in time, exhibiting a synchronous excitation. The PSM with M=1 accurately represents the forcing due to the inlet distortion downstream of the fan for the time evolution of the absolute total pressure in the frame of reference of the rotor. As could be foreseen, the most significant disturbance is associated with the first harmonic, which is the one imposed at the inlet. It can be noticed that the M=1 and the full-annulus case display slight differences, which nearly vanish when more spatial harmonics are included (M=3). The authors believe that this slight variation of the solution with the number of harmonics raised by one of the reviewers is due to the generation of higher-order circumferential harmonics at the inlet since the total pressure is a nonlinear combination of the conservative variables.

Fig. 16
Unsteady signal (a) and harmonic content (b) of an absolute total pressure numerical probe (non-dimensionalized with the inlet total pressure variation, ΔP0), in the rotor frame of reference located 0.5 chords downstream of the fan (see Fig. 15) obtained by the PSM (M=1) (solid line) and a full-annulus simulation (dashed line) for the case kL=1 and NB=kc=22
Fig. 16
Unsteady signal (a) and harmonic content (b) of an absolute total pressure numerical probe (non-dimensionalized with the inlet total pressure variation, ΔP0), in the rotor frame of reference located 0.5 chords downstream of the fan (see Fig. 15) obtained by the PSM (M=1) (solid line) and a full-annulus simulation (dashed line) for the case kL=1 and NB=kc=22
Close modal

Figure 16(b) shows the harmonic content of the two converged solutions in a periodic sense. The harmonic content of both simulations exhibits sidebands, with Fourier indices ks= 21 and 23, due to the scattering of the modulating (kc=22) and the carrier (kL=1) harmonics. The carrier harmonic has an amplitude close to 1 due to the distorted inflow. The amplitude of these scattered harmonics is small since the modulation index β=ϵkL/kc is rather small. It is important to recall that the amplitude of the first sideband in the varying speed linear equation is β/2 and hence, even if ϵ=1, the amplitude of the side bands is estimated in about a fortieth of the fundamental harmonic. It is worth noting that unlike in the LWE with varying speed in this realistic case, the amplitude of the two sidebands is not the same. The higher-order sidebands, ks=nkc±m, are very small for their rapid amplitude decay.

The good agreement between both simulations allows us to conclude that the block spectral method captures the Tyler–Sofrin modes created by the interaction between the rotor and the inlet distortion though these modes are neither explicitly modeled nor fit exactly within a single blade passage. In the authors’ opinion, this is a remarkable result that stresses the powerfulness of the block spectral method. It is concluded as well that in this case a single spatial harmonic suffices to reproduce the full-annulus solution.

10 Conclusions

An analysis of the block spectral method for turbomachinery applications has been conducted. The details of the method and its analysis have been first described by resorting to the one-dimensional wave equation. It is concluded that the reconstructed solution of the full-domain using the block spectral method has the same level of accuracy as the original full-annulus simulation if enough spatial harmonics are retained.

It has been shown that the passage-spectral method can accurately capture three different families of spatial harmonics. The first two, namely, the fundamental harmonics, explicitly retained in the Fourier representation of the long wavelengths, and the embedded harmonics self-contained in the blocks were expected and come as no surprise. However, the Tyler–Sofrin modes (i.e., the scattered harmonics) are natively retained by the method as well. This fact had not been noticed before by other authors. This means that the method is fully compatible with the spatial harmonic spectrum created in rotor/stator interaction problems (including the Tyler–Sofrin modes), further extending the applicability of the block spectral method. These scattered modes are longer than the airfoil pitch (and therefore not contained in the passage) but do not coincide with the Fourier representation of the longest wavelengths. They are a hybridization of a spectral and a finite volume representation. Adding higher harmonics of the fundamental mode of the spectral representation does not lead directly to including the Tyler–Sofrin modes.

The block spectral method has been finally implemented in an existing implicit finite volume Navier–Stokes solver and its generalization for turbomachinery problems described in some detail. The method does not require any hypothesis about the temporal periodicity of the flow, as shown in Refs. [12,22], where the method was applied to successfully predict the instability onset and the transient of a fan subject to periodic perturbations. Moreover, the method can be easily implemented in time-marching CFD solvers with minimum modifications on the flux calculation routine since the method needs to reconstruct the information just in the boundaries with the neighboring blocks.

Finally, it has been verified, using 3D simulations, that the Tyler–Sofrin modes are natively retained in the method, and that the matching with full-annulus simulations is very good. This observation is new and surprising at first glance.

Acknowledgment

This work has been funded by Ministerio de Ciencia e Innovación under the grant RTC2019-007194-4 (TASTE) from RETOS-Colaboración program.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix: Harmonic Content of a Particular Solution of the Varying Speed Wave Equation

The particular solution derived for the varying speed linear wave equation
u(x,t)=u0sin(kL(xc0t)+βsin(kcx))
(A1)
has a large degree of similarity with the signals created in frequency modulation (FM) where the signals are described as y(t)=sin[Θ(t)], where Θ is the instantaneous phase. The instantaneous frequency in a FM signal consisting of a constant carrier frequency, ωc, plus an usually sinusoidal modulation with frequency, ωm, is
ω(t)=ωc+ϵωccos(ωmt+ϕ)
(A2)
where ϵ is the dimensionless amplitude of the modulation. The phase is given by
Θ=0tω(t)dt=ωct+βsin(ωmt+ϕ)
where β=ϵωc/ωm is the modulation index. Therefore, an FM signal has the form
y(t)=sin[ωct+βsin(ωmt+ϕ)]
(A3)
which has a very similar form to the particular solution of the linear wave equation with varying speed (Eq. (A1)). The Fourier transform of the function y(t) is given in Ref. [23]. Replacing ωct by kL(xc0t), ωmt by kcx, β by β, setting ϕ=0, and using Eq. (111) in [23], u(x,t) can be written as
u(x,t)u0=J0(β)sin[kL(xc0t)]+J1(β){sin[(kLkc)xkLc0t]sin[(kL+kc)xkLc0t]}+J2(β){sin[(kL2kc)xkLc0t]+sin[(kL+2kc)xkLc0t]}+Jn(β){sin[(kLnkc)xkLc0t]+(1)nsin[(kL+nkc)xkLc0t]}
(A4)
which is nothing else than its Fourier transform where β=ϵkL/kc and Jn is the nth Bessel function of the first kind. It can be seen that the carrier spatial harmonic term, kL, has an infinite number of sidebands kL±nkc where kc is the modulating spatial harmonic. The amplitude of the sidebands is given by the coefficients Jn(β) and its variation is not a monotonic function with the Fourier index n unless the modulation index is small. If β1, the amplitude of the sidebands decreases as βn and the trend of the sideband amplitudes can be easily predicted. Moreover, in the limit of β small J0(β)1 and J1(β)β/2, whereas the rest of the sidebands are much smaller and can be disregarded. This regime is known as well as the narrowband FM limit.

Figure 17 compares the harmonic content of the solution obtained by a numerical Fourier transform of Eq. (A1) with the coefficients of Eq. (A4). The first observation is that the matching between both is exact as could be expected. On the second hand, it can be seen that if β<π/2, the amplitude of the sidebands decreases monotonically while for higher values the trend is non-monotonic, and the amplitude of the sidebands can be higher than that of the carrier or fundamental harmonic.

Fig. 17
Spatial harmonic content of Eq. (A1) for kL=1 and kc=6 and different values of the modulation index. Black: numerical evaluation of the Fourier transform. Gray: coefficients of Eq. (A4).
Fig. 17
Spatial harmonic content of Eq. (A1) for kL=1 and kc=6 and different values of the modulation index. Black: numerical evaluation of the Fourier transform. Gray: coefficients of Eq. (A4).
Close modal
An interesting property of the solution of the varying speed wave equation (Eq. (A4)) is given by Neumann’s addition theorem [24], which states that
1=J02(0)=J02(β)+2n=1Jn2(β)
(A5)
This effectively means that the energy of the scattered particular solution is the same as the one of an initial condition containing just the fundamental harmonic. This is a direct consequence of the fact that the system conserves energy.

References

1.
Choi
,
M.
,
Smith
,
N. H.
, and
Vahdati
,
M.
,
2013
, “
Validation of Numerical Simulation for Rotating Stall in a Transonic Fan
,”
ASME J. Turbomach.
,
135
(
2
), p.
021004
.
2.
Lee
,
K.-B.
,
Wilson
,
M.
, and
Vahdati
,
M.
,
2017
, “
Numerical Study on Aeroelastic Instability for a Low-Speed Fan
,”
ASME J. Turbomach.
,
139
(
7
), p.
071004
.
3.
Dowell
,
E. H.
, and
Hall
,
K. C.
,
2001
, “
Modeling of Fluid–Structure Interaction
,”
Annu. Rev. Fluid Mech.
,
33
(
1
), pp.
445
490
.
4.
He
,
L.
,
2010
, “
Fourier Methods for Turbomachinery Applications
,”
Prog. Aerosp. Sci.
,
46
(
8
), pp.
329
341
.
5.
He
,
L.
,
2006
, “
Fourier Modeling of Steady and Unsteady Nonaxisymmetrical Flows
,”
J. Propul. Power
,
22
(
1
), pp.
197
201
.
6.
He
,
L.
,
2011
, “
Efficient Computational Model for Nonaxisymmetric Flow and Heat Transfer in Rotating Cavity
,”
ASME J. Turbomach.
,
133
(
2
), p.
021018
.
7.
He
,
L.
,
2011
, “
Block-Spectral Approach to Film-Cooling Modeling
,”
ASME J. Turbomach.
,
134
(
2
), p.
021018
.
8.
He
,
L.
,
2013
, “
Block-Spectral Mapping for Multi-Scale Solution
,”
J. Comput. Phys.
,
250
, pp.
13
26
.
9.
He
,
L.
,
2013
, “
Fourier Spectral Modelling for Multi-Scale Aero-Thermal Analysis
,”
Int. J. Comput. Fluid Dyn.
,
27
(
2
), pp.
118
129
.
10.
Stapelfeldt
,
S. C.
, and
Di Mare
,
L.
,
2015
, “
Reduced Passage Method for Multirow Forced Response Computations
,”
AIAA J.
,
53
(
10
), pp.
3049
3062
.
11.
He
,
L.
,
2018
, “
Multiscale Block Spectral Solution for Unsteady Flows
,”
Int. J. Numer. Methods Fluids
,
86
(
10
), pp.
655
678
.
12.
Romera
,
D.
, and
Corral
,
R.
,
2020
, “
Efficient Passage-Spectral Method for Unsteady Flows Under Stall Conditions
,”
ASME J. Turbomach.
,
142
(
12
), p.
121007
.
13.
Romera
,
D.
, and
Corral
,
R.
,
2021
, “
Nonlinear Stability Analysis of a Generic Fan With Distorted Inflow Using Passage-Spectral Method
,”
ASME J. Turbomach.
,
143
(
4
), p.
061001
.
14.
Tyler
,
J. M.
, and
Sofrin
,
T. G.
,
1962
, “Axial Flow Compressor Noise Studies,” Technical Report, SAE Technical Paper.
15.
Burgos
,
M. A.
,
Contreras
,
J.
, and
Corral
,
R.
,
2011
, “
Efficient Edge-Based Rotor/Stator Interaction Method
,”
AIAA J.
,
49
(
1
), pp.
19
31
.
16.
Corral
,
R.
,
Gisbert
,
F.
, and
Pueblas
,
J.
,
2017
, “
Efficient Execution of a Parallel Edge-Based Navier–Stokes Solver on Graphics Processing Units
,”
Int. J. Comput. Fluid Dyn.
,
31
(
2
), pp.
1
16
.
17.
Jameson
,
A.
, “
Time Dependent Calculations Using Multigrid, With Applications to Unsteady Flows Past Airfoils and Wings
,”
AIAA 10th Computational Fluid Dynamics Conference
,
Honolulu, HI
,
June 24–26
, p.
1596
.
18.
Zhang
,
W.
, and
Vahdati
,
M.
,
2018
, “
A Parametric Study of the Effects of Inlet Distortion on Fan Aerodynamic Stability
,”
ASME J. Turbomach.
,
141
(
1
), p.
011011
.
19.
Lee
,
K.-B.
,
Wilson
,
M.
, and
Vahdati
,
M.
,
2019
, “
Effects of Inlet Disturbances on Fan Stability
,”
ASME J. Eng. Gas Turbines Power
,
141
(
5
), p.
051014
.
20.
Wenqiang
,
Z.
,
Stapelfeldt
,
S.
, and
Vahdati
,
M.
,
2020
, “
Influence of the Inlet Distortion on Fan Stall Margin at Different Rotational Speeds
,”
Aerosp. Sci. Technol.
,
98
(
4
), p.
105668
.
21.
Burgos
,
M. A.
,
Chia
,
J. M.
,
Corral
,
R.
, and
López
,
C.
,
2010
, “
Rapid Meshing of Turbomachinery Rows Using Semi-unstructured Multi-block Conformal Grids
,”
Eng. Comput.
,
26
(
4
), pp.
351
362
.
22.
Romera
,
D.
, and
Corral
,
R.
,
2021
, “
Nonlinear Stability Analysis of a Generic Fan With Distorted Inflow Using Passage-Spectral Method
,”
ASME J. Turbomach.
,
143
(
6
), p.
061001
.
23.
Hartmann
,
W. M.
,
1995
,
Handbook of Perception and Cognition
, 2nd ed.,
Hearing, Academic Press Inc.
,
San Diego, CA
, p.
32
.
24.
Abramowitz
,
M.
, and
Stegun
,
I. A.
,
1965
,
Handbook of Mathematical functions: With Formulas, Graphs, and Mathematical Tables
,
Dover Publications
,
Washington, DC
.