# Simplified Eigenvalue Analysis for Turbomachinery Aerodynamics with Cyclic Symmetry

Eigenvalue analysis is widely used for linear instability analysis in both external and internal aerodynamics. It typically involves finding the steady state, linearizing around it to obtain the Jacobian, and then solving for its eigenvalues and eigenvectors. When the flow is modelled with Reynolds-averaged Navier–Stokes equations with a large boundary-layer-resolving mesh, the resulting eigenvalue problem can be of very high dimensions, and is thus computationally very challenging. To reduce the computational cost, a simplified approach is proposed to compute the eigenvalues and eigenvectors, by exploiting the cyclic symmetric nature of annular fluid domain for typical compressors. It is shown that via a rotational transformation, the Jacobian can be reduced to a block circulant matrix, whose eigenvalues and eigenvectors then can be computed using only one sector of the entire domain. This simplified approach significantly lowers the memory overhead and the CPU time of the eigenvalue analysis without compromising on the accuracy. The proposed method is applied to the eigenvalue analysis of an annular compressor cascade with 22 repeated sectors and it is shown the spectrum of the whole annulus can be obtained by using the information of 1 sector only, demonstrating the effectiveness of the proposed method.

## Authors

• 1 publication
08/07/2020

### Solving two-parameter eigenvalue problems using an alternating method

We present a new approach to compute selected eigenvalues and eigenvecto...
03/18/2020

### Computation of Tight Enclosures for Laplacian Eigenvalues

Recently, there has been interest in high-precision approximations of th...
03/01/2021

### Computation of transmission eigenvalues by the regularized Schur complement for the boundary integral operators

This paper is devoted to the computation of transmission eigenvalues in ...
10/22/2020

### One-shot Distributed Algorithm for Generalized Eigenvalue Problem

Nowadays, more and more datasets are stored in a distributed way for the...
02/02/2020

### Variational projector-augmented wave method: a full-potential approach for electronic structure calculations in solid-state physics

In solid-state physics, energies of crystals are usually computed with a...
04/04/2022

### Learning Linear Symmetries in Data Using Moment Matching

It is common in machine learning and statistics to use symmetries derive...
12/01/2020

### Enhancing Scalability of a Matrix-Free Eigensolver for Studying Many-Body Localization

In [Van Beeumen, et. al, HPC Asia 2020, https://www.doi.org/10.1145/3368...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Eigenvalue analysis based on the Reynolds-averaged Navier–Stokes (RANS) equations is a widely used approach to study the linear instability of external flow problems [theofilis2011global] such as the transonic shock buffet for airfoil [Crouch2007Predicting, Crouch2009Origin] and for wings [2018arXivTimme, paladini2019transonic, crouch2019global]. In those works, the RANS equations are linearized about the steady state and an eigenvalue problem is solved for the resulting large sparse linear system of equations. The stability of the system is then determined by the existence of eigenvalues with positive real parts. This approach is shown to be capable of determining the stability boundary that is consistent with unsteady RANS approach, but at orders-of-magnitude lower a cost.

Eigenvalue-based linear instability analysis has also been used for turbomachinery applications, mainly regarding the flow-instability inception towards the stall boundary of compressors. A general approach to predict the onset of such instability using the eigenvalue method is proposed in [Sun2013A] and has been successfully applied to both axial and centrifugal compressors [liu2014basic, sun2016flow]. Different from the application on airfoil and wing aerodynamics where the flow is modelled with three-dimensional RANS equations, these works on compressor stability study use a simplified approach by modelling the effect of the blades with body forces and ignoring the circumferential variation of the flow field. With such simplification, a small eigenvalue problem can be formulated and solved in the meridional plane to predict the linear instability of the system. Although the method is shown to predict the turbomachinery flow-instability onset with plausible accuracy, ignoring the three-dimensional details naturally brings error that is not known a priori. Besides, modern compressors are designed with increasingly higher loading and exhibit strong three-dimensional effect, which needs to be accounted for when studying the flow instability [pullan2015origins]. Furthermore, three-dimensional RANS calculation has become routine practice for turbomachinery industry and consequently the flow-instability inception prediction needs to be based on RANS equations as well in order to retain the same modelling capability.

One key factor limiting the application of the RANS-based eigenvalue analysis to turbomachinery is its high computational cost. Although the steady state performance of compressors can be computed using a single passage mesh, which typically has grid points on the order of one million, whole annulus domain needs to be employed in order to capture the instability onset. Consequently, a computational mesh one to three orders of magnitude larger than typical airfoil/wing applications is needed, and therefore computing eigenvalue (even a small subset of it) based on RANS equations for compressors is computationally much more demanding.

Turbomachinery components usually assume cyclic symmetric shape and so is the flow field inside, and the cyclic symmetry can be exploited in order to simplify the computation. It is in fact established practice for computing the natural frequencies and vibrational mode shapes of periodic structures with cyclic structure [thomas1979dynamics]. One first meshes one sector of the entire structure and specifies the nodal diameter and the natural frequencies and mode shapes for the specifies nodal diameters can be computed. Looping over all possible nodal diameters allows one to obtain all natural frequencies and mode shapes for the entire structure. Throughout the whole process, all computations are performed based on one sector of the domain, thus requires much less computational resource and much shorter CPU time, compared to the whole-annulus computation. Such simplification approach was revisited in [Schmid2017Stability] for flow problems with translational periodicity and was successfully demonstrated on linear cascade type of problems. However, the methodology proposed therein is not directly applicable to cases with rotational periodicity, which are representative of real compressors. Therefore, it is desirable to extend the methodology to flow problems with rotational periodicity in order to significantly simplify the eigenvalue computation in order to study the whole annulus compressor flow instability.

In this work, an extension of the method proposed in [Schmid2017Stability] is developed so that it can be directly applied to the eigenvalue calculation for cyclic symmetric domains using one sector of the geometry only. In the remaining part of the paper, the theoretical basis for the simplification of the eigenvalue computation is first reviewed in Sec. 2 with particular focus on the linearized RANS equations. The application of the proposed method is then applied to the eigenvalue calculation for an annular compressor cascade case to demonstrate its effectiveness in reducing the computational cost and results are discussed in Sec. 3. Conclusions are given in Sec. 4.

## 2 Theoretical framework

### 2.1 Eigenvalue analysis of block circulant matrices

A circulant matrix has the form as follows

 B=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣b0b1⋯bM−1bM−1b0⋯bM−2⋮⋮⋱⋮b1b2⋯b0⎤⎥ ⎥ ⎥ ⎥ ⎥⎦.

Due to its cyclic symmetric nature, it can be verified that matrix has the following eigenvalues and eigenvectors

 λm=M−1∑k=0bkρkm

and

 vm=[1,ρm,ρ2m,⋯,ρM−1m]T  (m=0,1,2,…,M−1),

where

 ρm=e j2πm/M  (m=0,1,2,…,M−1)

are the th complex roots of 1, and j is the imaginary unit. It can easily be verified that all eigenvectors are orthogonal to each other, i.e.,

 vm1⊥vm2  if  m1≠m2.

Matrix becomes a block circulant matrix when each is an square matrix. It has long been discovered that the computation of the eigenvalue and eigenvectors of a block circulant matrix can also be simplified by exploiting its cyclic symmetric nature[Friedman1961Eigenvalues]. It can be shown that all the eigenvector of the block circulant matrix can be divided into groups, with each group containing vectors. The th eigenvector in the th group, denoted as , should be of the following form

 wm,n=[1,ρm,ρ2m,⋯,ρM−1m]Tvm,n

where is a vector of length . In order for to be an eigenvector, the following equation has to be satisfied

 Bwm,n=λm,nwm,n

which can be expanded as follows

 (b0+ρmb1+⋯+ρM−1mbM−1)vm,n=λm,nvm,n(bM−1+ρmb0+⋯+ρM−1mbM−2)vm,n=λm,nρmvm,n⋮(b1+ρmb2+⋯+ρM−1mb0)vm,n=λm,nρM−1mvm,n.

It can be verified that any two equations in the linear system above are equivalent. Therefore and can be found by solving the first equation only. That is, they are the th eigenvalue and eigenvector of the matrix

 Bm=b0+ρmb1+ρ2mb2+⋯+ρM−1mbM−1.

Therefore, the eigen analysis of the large matrix of dimension can be obtained via the eigen analysis of smaller matrix of dimension , significantly simplifying the computation.

### 2.2 Jacobian matrix of the linearized RANS equations

For a typical finite volume solver based on the RANS equations, the vector that satisfies

 R(U)=0

is the steady state solution. The time-dependent behavior of the flow solution is governed by the unsteady form

 dUdt=R(U).

To study its linear stability, one could linearize the nonlinear residual vector about the steady state with respect to the flow solution vector to obtain the Jacobian matrix

 A:=∂R∂U

and compute its eigenvalues and eigenvectors.

To facilitate the discussion, it is assumed that the computational mesh is cyclic symmetric with -periodicity. Consequently, the full-annulus mesh can be divided into non-overlapping zones, each of which containing grid points, as illustrated in Fig. 1. In addition, the grid points in each zone is arranged in such a way that the th point in zone is the th point in zone rotated by an angle of .

Due to the cyclic symmetry of the mesh, the Jacobian matrix can be partitioned as a block matrix as follows

 A=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∂R0∂U0∂R0∂U1…∂R0∂UM−1∂R1∂U0∂R1∂U1…∂R1∂UM−1⋮⋮⋱⋮∂RM−1∂U0∂RM−1∂U1…∂RM−1∂UM−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

where and are the vectors of the flow variables and the nonlinear residual for the grid points in the th zone, each being a vector of length for Euler/laminar or for turbulent flow with a one-equation turbulence model, and each block is a or matrix.

Note that block matrices

are second-order tensors that are invariant to the coordinate system used, and can be expressed using tensor product as

 ∂Rm1∂Um2=[Mm1,m2]i,jeiej=[^Mm1,m2]i,j^ei^ej (1)

where are the basis vectors in the coordinate frame fixed to the computational mesh and is short for the tensor product . and are the matrix representation of the second-order tensor in the coordinate systems and respectively. Let be a coordinate system that is obtained by rotating by , as shown in Fig. 2

The basis vectors of the two coordinate systems can linked via the transformation matrix via

 ei[Tm1]i,j=^ej,

where the matrix denotes the rotational transform of . Substitute the in Eq. (1) yields

 Mm1,m2=Tm1^Mm1,m2T−m1.

In addition, it is straightforward to see that viewed in is identical to viewed in . This means the matrix representation of the former in is that of the latter in , i.e.,

 ^Mm1,m2=M0,m2−m1

and hence

 Mm1,m2=Tm1M0,m2−m1T−m1.

or equivalently,

 ∂Rm1∂Um2=Tm1∂R0∂Um2−m1T−m1.

The Jacobian matrix becomes

 A=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∂R0∂U0∂R0∂U1…∂R0∂UM−1T∂R0∂UM−1T−1T∂R0∂U0T−1…T∂R0∂UM−2T−1⋮⋮⋱⋮TM−1∂R0∂U1ˇ−(M−1)TM−1∂R0∂UM−2T−(M−1)…TM−1∂R0∂U0T−(M−1)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

Further more, change of variable is used for

 ~Um=T−mUm

so that

 ∂R0∂Um=∂R0∂~UmT−m

The Jacobian matrix can be further rewritten to be

 A =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∂R0∂~U0∂R0∂~U1T−1…∂R0∂~UM−1T−(M−1)T∂R0∂~UM−1T∂R0∂~U0T−1…T∂R0∂~UM−2T−(M−1)⋮⋮⋱⋮TM−1∂R0∂~U1TM−1∂R0∂~UM−2T−1…TM−1∂R0∂~U0T−(M−1)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=TBT−1

with matrices and defined as

Therefore, the Jacobian matrix is similar to a block-circulant matrix and thus they have the same eigenvalues. To compute the eigenvalues and eigenvectors of the matrix , we first compute and , th eigenvalue and eigenvector of the matrix

 Bm=∂R0∂~U0+ρm∂R0∂~U1+ρ2m∂R0∂~U2+⋯+ρM−1m∂R0∂~UM−1,

and the eigenvectors of , , can be computed as

 wm,n=[1,ρm,ρ2m,⋯,ρM−1m]Tvm,n.

Finally, the eigenvectors of the Jacobian matrix are , which can be calculated as follows

 Twm,n=[1,ρmT,ρ2mT2,⋯,ρM−1mTM−1]Tvm,n.

It can be seen that for any of the eigenvectors of , , the corresponding eigenvectors of , are modes with different nodal diameters. Same as for the circulant block matrix, to compute the eigenvalues and eigenvectors of the Jacobian matrix , one only needs to solve the eigenvalue problem of dimension , significantly reducing the computational cost, especially for large CFD analysis.

## 3 Results

The test case used to demonstrate the usefulness of the proposed method is an annulus compressor cascade case. The cascade is produced by taking the surface of revolution at the 50% height of the NASA Rotor 67 [strazisar1989laser]. The three-dimensional mesh has one cell in the radial direction and the subsequent analysis is thus a quasi-3D one. The compressor cascade rotates at 16,043 Revolutions per minute. The annular cascade has 22 passages in total. The mesh for the whole annulus is produced by meshing a single passage first and then copying the mesh to the whole annulus. The mesh for the entire domain has 126,808 grid points, with an average of 5,764 per zone.

The turbomachinery nonlinear flow solver NutsCFD is used to compute the performance of the cascade. The speedline, computed by gradually raising the back pressure until the steady state can no long converge, is shown in Fig. 3. The fully-converged steady state flow solution corresponding to the condition marked in red in Fig. 3 (to be called condition ‘A’ in the remainder of the paper) is visualized in Fig. 4, where the shock/boundary layer interaction can clearly be seen. All steady state calculations are performed using a single passage. The visualized flow fields are produced by copying the solution in the single passage to the whole annulus.

#### 3.0.2 Eigenvalue analysis

To perform the eigenvalue analysis for condition ‘A’, we first compute the Jacobian matrix for the single passage and then exploiting the periodic boundary condition in the steady solver to obtain and . Note that for RANS flow solver, are all zeros, so can be simplified as

 Bm=∂R0∂~U0+ρm∂R0∂~U1+ρM−1m∂R0∂~UM−1. (2)

It is well known both numerically and experimentally that the least stable modes when instability appears in rotating flows, the characteristic frequency of the instability is close to the rotating frequency, which is for the test case used. Consequently, the matrix is scaled by so that the imaginary part of the eigenvalue is equivalent to the frequency of the corresponding eigenmode expressed in terms of engine order (EO). Since the Jacobian matrix is of moderate size (of dimension , and with nearly million non-zero entries), the sparse matrix eigenvalue solver, “eigs”, in SciPy [jones2001scipy] is use on a single-processor computer. The “eigs” function in SciPy is a wrapper to the widely used eigenvalue solver library ARPACK [lehoucq1998arpack].

Two different methods (listed in Tab. 2) are used to compute the eigenvalues near the imaginary axis. Method 1 computes the eigenvalues of the Jacobian matrix for the whole annulus of 22 sectors, . We apply complex shifts of and and 30 interior eigenvalues are computed for each shift. A total of 90 eigenvalues can computed and they are shown in Fig. 5 as red crosses.

Method 2 computes the eigenvalues using 1 sector only. Nevertheless, it admits eigenvalues for the whole annular domain by setting in Eq. (2) to . Same as method 1, interior eigenvalues are computed using three different shifts for each . For each shift, 2 eigenvalues are computed, resulting in a total of 132 eigenvalues (some of them appear more than once for different shifts). All eigenvalues computed using method 2 are also plotted in Fig. 5 as blue dots. If can be seen that the least stable eigenvalues (closest to the imaginary axis) are located around and a zoomed view of the eigenvalues in the nearly region is shown on the right side in Fig. 5. It can be seen that both methods can capture a relevant subset of the eigenvalues and are both sufficient in order to study the linear stability of the underlying flow problem. However, method 2 is based only on one sector of the whole annulus domain, and thus has a significantly lower memory overhead and also costs much less CPU time than method 1.

Another advantage of method 1, besides the lower computational cost, is that each eigenvalue/eigenvector computed has a known nodal diameter, which is solely determined by the value it is computed with. The nodal diameter information reveals the circumferential spatial periodicity of each particular eigenmode.

## 4 Conclusions

A method to simplify the eigenvalue analysis for flow problems in circular cyclic symmetric domains are proposed. The method first transform the Jacobian matrix for the whole annulus into a block circulant matrices, and subsequently reduces the eigenvalue problem of the whole circular domain to one sector only. This reduces the size of the original eigenvalue problem by a factor of , being the number of passages/blades. The proposed method is applied to the eigenvalue analysis for an annular compressor cascade with sectors, and it is shown that by using sector only, spectrum of the whole domain can be obtained at a much lower memory as well as CPU time cost, demonstrating the advantage of the proposed method. Future work will explore the application of the method to realistic three-dimensional compressors to study the flow instability problems. In addition, the eigenvectors obtained can be used to construct reduced-order models for parametric study and optimization for rotating-flow instability.