## 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.

## 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

## 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, $\lambda c=2\pi /kc$, disturbance superimposed with a long-wavelength, $\lambda L=2\pi /kL$, and perturbation having the latter, a much longer wavelength than the former ($\lambda L\u226b\lambda 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 $\lambda L\u223c\lambda 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.

The question that is addressed in this work is whether this problem can be tackled solving just a few small blocks of size $\lambda 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 $\lambda c$ discretized with $NP$ mesh points each. The block size, which is equal to $2\pi /NB$, is set to contain the short-wavelength $\lambda 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\pi /(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.

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

The mean value or the $0$th harmonic.

$N=(NP\u22121)/2$ block-embedded harmonics, $k=nkc$, where $n$ is an integer such that $n\u2264N$. 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$.

$M$ external block modulating harmonics with wavenumber $k=m$, where $m$ is an integer such that $m\u2264M$. 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.

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

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 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 $j$th 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.

## 6 Spectral View of Block Reconstruction

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^\u2212k,m=u^k,m*$ but $u^k,\u2212m\u2260u^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$.

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\xb1m$, since the method is designed in such a way that $kc\u226b1$. 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\Delta t/\Delta x=0.1$ and $\mu =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.

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 $j$th 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.

### 7.2 Outer Harmonic Superimposed With Embedded Harmonics Case.

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\u2264(NP\u22121)/2$.

FULL | $PSM(M=1)$ | $|FULL\u2212BSM|$ | |
---|---|---|---|

$kL$ | $0.199212$ | $0.199213$ | $5.02\xd710\u22126$ |

$kc$ | $0.433871$ | $0.433869$ | $4.61\xd710\u22126$ |

$2kc$ | $0.215419$ | $0.215421$ | $9.28\xd710\u22126$ |

FULL | $PSM(M=1)$ | $|FULL\u2212BSM|$ | |
---|---|---|---|

$kL$ | $0.199212$ | $0.199213$ | $5.02\xd710\u22126$ |

$kc$ | $0.433871$ | $0.433869$ | $4.61\xd710\u22126$ |

$2kc$ | $0.215419$ | $0.215421$ | $9.28\xd710\u22126$ |

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.

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.

FULL | $PSM(M=1)$ | $|FULL\u2212PSM|$ | |
---|---|---|---|

$kL$ | $0.198187$ | $0.198186$ | $5.05\xd710\u22126$ |

$kc$ | $0.434122$ | $0.434119$ | $6.91\xd710\u22126$ |

$2kc$ | $0.214259$ | $0.214257$ | $9.33\xd710\u22126$ |

FULL | $PSM(M=1)$ | $|FULL\u2212PSM|$ | |
---|---|---|---|

$kL$ | $0.198187$ | $0.198186$ | $5.05\xd710\u22126$ |

$kc$ | $0.434122$ | $0.434119$ | $6.91\xd710\u22126$ |

$2kc$ | $0.214259$ | $0.214257$ | $9.33\xd710\u22126$ |

### 7.4 “Scattered” Mode Case.

A case to verify the accuracy of the scattered or sideband harmonics $nkc\xb1m$ 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 $\alpha =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=kc\u2212m$, 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.

### 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, $\alpha =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).

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.

Figure 11 shows the solution of Eq. (23) for $kL=1$, $kc=6$, and $\beta =0.5$. The solution can be expressed as a function of Bessel functions of the first kind $Jn(\beta )$ (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.

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 $\u03f5=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\pi /kc$ and the harmonic $kL$ injected in the initial condition.

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 $\beta $. In the cases studied in this work, the modulation index is very small, $\beta \u226a1$, since we are targetting the limit $kL\u226akc$, and then the amplitude of the carrier harmonic $kL$ is $J0(\beta )$, 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\xb1m$ whose amplitude in the limit $\beta \u226a1$ is $J1(\beta )\u2243\beta /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\xb1mkL$, are very small because their amplitudes decay as $O(\beta n)$ if $\beta \u21920$. 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 $\beta \u223cO(1$), though such cases have not been illustrated here.

Table 3 compares the amplitude of the fundamental, $kL$, and the lower scattered harmonics, $kc\xb1kL$, of the solutions obtained using the full-annulus domain and the BSM with $M=1$ for the case $\u03f5=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.

FULL | $PSM(M=1)$ | $|FULL\u2212BSM|$ | |
---|---|---|---|

$kL$ | $0.990302$ | $0.990304$ | $2.02\xd710\u22126$ |

$kc\u2212kL$ | $0.0661679$ | $0.0661682$ | $4.53\xd710\u22126$ |

$kc+kL$ | $0.0661696$ | $0.0661698$ | $3.02\xd710\u22126$ |

FULL | $PSM(M=1)$ | $|FULL\u2212BSM|$ | |
---|---|---|---|

$kL$ | $0.990302$ | $0.990304$ | $2.02\xd710\u22126$ |

$kc\u2212kL$ | $0.0661679$ | $0.0661682$ | $4.53\xd710\u22126$ |

$kc+kL$ | $0.0661696$ | $0.0661698$ | $3.02\xd710\u22126$ |

## 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 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.

The solution of this equation is sped-up using steady-state acceleration techniques, marching in $\tau $ instead of $t$. The preconditioning matrix $Pi\u22121$ 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.

*N*passages, hence the name PSM.

### 9.2 Three-Dimensional Block Spectral Method.

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.

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.

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\xaf0+\Delta P0sinkL\theta $, being $kL=1$ the wavenumber of the harmonic perturbation in the azimuthal direction. The amplitude of the inlet stagnation pressure variation is $\Delta P0=5kPa$ which corresponds to about 5% of the mean inlet total pressure, $P\xaf0$. The inlet flow nonuniformity gives rise to a strong forcing on the fan blades. Inlet distortion can also severely compromise fan stability [18–20]. 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\xd7106$ nodes per passage. The size of the full-annulus mesh is approximately $50\xd7106$ nodes. The Reynolds number of the simulations is very high ($Re\u223c106$), and the standard $k$–$\omega $ 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)).

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\xb1m=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.

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 $\beta =\u03f5kL/kc$ is rather small. It is important to recall that the amplitude of the first sideband in the varying speed linear equation is $\beta /2$ and hence, even if $\u03f5=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\xb1m$, 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

*n*th Bessel function of the first kind. It can be seen that the carrier spatial harmonic term, $kL$, has an infinite number of sidebands $kL\xb1nkc$ where $kc$ is the modulating spatial harmonic. The amplitude of the sidebands is given by the coefficients $Jn(\beta $) and its variation is not a monotonic function with the Fourier index $n$ unless the modulation index is small. If $\beta \u226a1$, the amplitude of the sidebands decreases as $\beta n$ and the trend of the sideband amplitudes can be easily predicted. Moreover, in the limit of $\beta $ small $J0(\beta )\u22431$ and $J1(\beta )\u2243\beta /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 $\beta <\pi /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.