## Abstract

Over the past years, nonlinear frequency-domain methods have become a state-of-the-art technique for the numerical simulation of unsteady flow fields within multistage turbomachinery as they are capable of fully exploiting the given spatial and temporal periodicities, as well as modeling flow nonlinearities in a computationally efficient manner. Despite this success, it still remains a significant challenge to capture nonlinear interaction effects within the context of configurations with multiple fundamental frequencies. If all frequencies are integer multiples of a common fundamental frequency, the interval spanned by the sampling points typically resolves the period of the common base frequency. For configurations in which the common frequency is very low in relation to the frequencies of primary interest, many sampling points are required to resolve the highest harmonic of the common fundamental frequency and the method becomes inefficient. In addition, when a problem can no longer be described by harmonic perturbations that are integer multiples of one fundamental frequency, as it may occur in two-shaft configurations or when simulating the nonlinear interaction in the context of forced response or flutter, then the standard discrete Fourier transform is no longer suitable, and the basic harmonic balance (HB) method requires extension. One possible approach is to use almost periodic Fourier transforms with equidistant or nonequidistant time sampling. However, the definition of suitable sampling points that lead to well-conditioned Fourier transform matrices and small aliasing errors is an intricate issue and far from straightforward. To overcome the issues regarding multifrequency problems described earlier, a new harmonic balance approach based on multidimensional Fourier transforms in time is presented. The basic idea of the approach is that instead of defining common sampling points in a common time period, separate time domains, one for each base frequency, are spanned, and the sampling points are computed equidistantly within each base frequency's period. Since the sampling domain is now extended to a multidimensional time-domain, all time instant combinations covering the whole multidimensional domain are computed as the Cartesian product of the sampling points on the axes. In a similar fashion, the frequency-domain is extended to a multidimensional frequency-domain by the Cartesian product of the harmonics of each base frequency, so that every point defined by the Cartesian product is an integer linear combination of the occurring frequencies. In this way, the proposed method is capable of fully integrating the nonlinear coupling effects between higher harmonics of different fundamental frequencies by using multidimensional discrete Fourier transforms within the harmonic balance solution procedure. The aim of this paper is to introduce the multidimensional harmonic balance method in detail and demonstrate the capability of the approach to simultaneously capture unsteady disturbances with arbitrary excitation frequencies. Therefore, the well-established aeroelasticity testcase standard configuration 10 in the presence of an artificial inflow disturbance, which mimics an upstream blade wake, is investigated. The crucial aspect of the proposed testcase is that a small ratio of the frequency of the inflow disturbance and the blades vibration frequency is chosen. To demonstrate the advantages of the newly proposed multidimensional harmonic balance approach, the results are compared to unsteady simulations in the time-domain and to state-of-the-art frequency-domain methods based on one-dimensional discrete Fourier transforms.

## Introduction

In today's turbomachinery design process, highly developed frequency-domain methods have become a state-of-the-art technique for the numerical simulation of unsteady flow fields within multistage turbomachinery. Its widespread success is related to the fact that the aerodynamic flow fields within multistage turbomachinery are mostly dominated by deterministic, periodic phenomena due to the relative motion between the adjacent blade rows. In the numerical simulation of such periodic flow phenomena, frequency-domain methods are an attractive alternative to the computationally expensive nonlinear time-domain simulations as they are capable of fully exploiting the given spatial and temporal periodicity, as well as modeling flow nonlinearities.

Over the past years, many contributions have been made to the development of frequency-domain methods that are capable of simulating time-periodic flow fields, i.e., time-linearized methods, where the governing equations for the flow perturbations decouple from those governing the mean flow field. To include the nonlinear interaction with the time-mean solution through the nonlinear stress terms, the nonlinear harmonic method that solves a system of a time-mean solution and a single complex harmonic perturbation has been proposed [1]. Furthermore, this method has been successfully extended to integrate the coupling between the higher harmonics [2]. Nevertheless, there is a need to explicitly model the nonlinear stress terms as they are essential for the simulation quality of the nonlinear harmonic method. To include the coupling relationship between the higher harmonics without modeling the nonlinear stress terms explicitly, a full harmonic balance (HB) formulation based on the use of the discrete Fourier transform (DFT) has been introduced [3]. The evaluation of the nonlinear residual of the unsteady Reynolds-averaged-Navier–Stokes equations is carried out in the time-domain, so that conventional routines from the steady solver can be used and the development is greatly simplified. The resulting equation system can be solved directly in the time-domain [3], as well as in the frequency-domain [4].

Despite this success, it still remains a significant challenge to capture nonlinear interaction effects within the context of configurations with multiple fundamental frequencies. If all frequencies are integer multiples of a common fundamental frequency, the interval spanned by the sampling points typically resolves the period of the common base frequency. For configurations in which the common frequency is very low in relation to the frequencies of primary interest, many sampling points are required to resolve the highest harmonic of the common fundamental frequency and the method may become inefficient. In addition, when a problem can no longer be described by harmonic perturbations that are integer multiples of one fundamental frequency, as it may occur in two-shaft configurations or when simulating the nonlinear interaction in the context of forced response or flutter, then the standard discrete Fourier transform is no longer suitable and the basic harmonic balance method requires extension. One possible approach is to use almost periodic Fourier transforms with equidistant or nonequidistant time sampling [5,6]. However, the definition of suitable sampling points that lead to well-conditioned Fourier transform matrices and small aliasing errors is an intricate issue and far from straightforward. To overcome the issues regarding multifrequency problems described earlier, a new harmonic balance approach based on multidimensional Fourier transforms in time is presented. The basic idea of the approach is that instead of defining common sampling points in a common time period, separate time domains, one for each base frequency, are spanned and the sampling points are computed equidistantly within each base frequency's period. Since the sampling domain is now extended to a multidimensional time-domain, all time instant combinations covering the whole multidimensional domain are computed as the Cartesian product of the sampling points on the axes. In a similar fashion, the frequency-domain is extended to a multidimensional frequency-domain by computing the Cartesian product of the harmonics of each base frequency, so that every point defined by the Cartesian product is an integer linear combination of the occurring frequencies. In this way, the proposed method is capable of fully integrating the nonlinear coupling effects between higher harmonics of different fundamental frequencies by using multidimensional discrete Fourier transforms within the harmonic balance solution procedure.

The aim of this paper is to introduce the multidimensional harmonic balance method in detail and demonstrate the capability of the approach to simultaneously capture unsteady disturbances with arbitrary excitation frequencies. To generate a multifrequency testcase, the well-established aeroelasticity testcase standard configuration 10 in the presence of an artificial inflow disturbance, which mimics an upstream blade wake, is investigated. The wake is defined as a sinusoidal time-periodic vorticity wave at the inflow of the computational domain, and the interaction of the time-dependent inflow condition with the unsteady flow field generated by the vibrating blade is analyzed. To demonstrate the advantages of the newly proposed multidimensional harmonic balance approach, the results are compared to unsteady simulations in the time-domain and to frequency-domain methods based on one-dimensional discrete Fourier transforms.

## Numerical Method

### Harmonic Balance Technique.

*t*denotes the physical time. If the spectrum of the unsteady flow field consists of arbitrary integer combinations of a finite number of fundamental frequencies, then the flow field can be termed quasi-periodic [7]. The frequency basis vector is denoted by

*k*is bounded by some upper and lower sidebands

^{m}*K*and $\u2212Km$. The

^{m}*m*th set of harmonic indices is given by

*ω*and the interblade phase angle

*σ*. Based on this, the unsteady Eq. (1) is reformulated in the frequency domain

*ω*with interblade phase angle

*σ*. Observe that since $R(q(t))$ is nonlinear, each term $R(q)\u0302(\omega ,\sigma )$ generally depends on all harmonics in $S$. The modeling of the nonlinearities in the function $R(q)\u0302(\omega ,\sigma )$ distinguishes the various HB methods proposed in the literature over the past decades [1,3,4,8]. The harmonic balance method proposed within this work is an alternating frequency-time domain approach, Fig. 1, using the discrete Fourier transform and its inverse to incorporate the nonlinearities implicitly [9]

where $F$ denotes the discrete Fourier transform and $F\u22121$ its inverse.

### Multidimensional Harmonic Balance Approach.

To avoid the intricate issue of defining appropriate sampling points to resolve higher harmonics of incommensurable fundamental frequencies and, at the same time, to include the nonlinear cross-coupling effects between higher harmonics of different fundamental frequencies, a new approach based on the multidimensional Fourier transform is proposed.

The fundamental idea of the multidimensional time harmonic balance approach is to exploit the fact that the one-dimensional Fourier coefficients $R(q)\u0302(\omega ,\sigma )$ coincide with the multidimensional Fourier coefficients of $R(q\u0303)$. To illustrate this fact, Fig. 2 shows a diagram with the one- and multidimensional flow and residual fields as well as their Fourier coefficients. Rather than proceeding along the lower row, the equality of the Fourier coefficients allows one to proceed along the upper one. Observe that $q\u0303$ and, thus, $R(q\u0303)$ are periodic in each time variable. Therefore, one can apply the MDFT to approximate the residual harmonic in the alternating frequency-time domain approach.

*M*fundamental frequencies. For every base frequency $\omega bm$, a separate frequency dimension is introduced. Each harmonic is defined by an

*M*-dimensional integer vector $k$ in the

*M*-ary Cartesian product over

*M*sets

*t*is spanned in a hypertime-domain $t$. The

^{m}*M*-dimensional time is represented by a hypertime vector $\tau $ of length

*M*

by introducing auxiliary $2\pi $-periodic phase variables *τ ^{i}*, one for each base frequency $\omega bm$.

*N*sampling points

^{m}*M*-dimensional time-domain is given by the

*M*-ary Cartesian product of the sampling points for each time axis

**using the reconstructed vector of conservative variables $q(x,r,t)$ is then given by**

*t*and $N=\u220fM\u22121Nm$ for every $m\u2208(0,\u2026,M\u22121)$.

For a two-dimensional time-domain, where the involved signals are periodic within their time dimension, the bi-dimensional time domain can be represented as a torus, Fig. 3. The unfolded view is depicted in Fig. 4 for two different frequency ratios. Figure 4(a) represents a frequency ratio of $f1/f2=3/2$. For each base frequency, a separate time axis is spanned; the period of each base frequency is marked on its corresponding axis. According to the Nyquist–Shannon sampling theorem, three sampling points are positioned within each base frequencies period on the axes. To cover the whole two-dimensional time-domain marked by the gray background, the Cartesian product of the sampling points on each axis is computed. Then, a two-dimensional Fourier transform can be applied. To highlight the difference of the multidimensional sampling, the one-dimensional sampling according to the Nyquist–Shannon sampling theorem based on the greatest common divisor of the two frequencies is as well depicted. In the case of a one-dimensional time-domain, the two time axes are equal with *t*^{1} = *t*^{2}, which corresponds to a line starting at $t1=t2=0$ with slope 1. Both frequencies are sampled at common time instants which by definition lay on the depicted *t*^{1} = *t*^{2} line. One can see that for the frequency ratio $f1/f2=3/2$, the one-dimensional DFT leads to less sampling points than the multidimensional DFT, although one has to sample a time period corresponding to the common period $Tbase=3T1=2T2$. This is as well visible in the torus view, Fig. 3. The trajectory indicates the progress of time in two dimensions. Since the frequency ratio is rational, after certain revolutions, the end of the trajectory reaches the origin. Figure 4(b) shows the same approach for a frequency ratio of $f1/f2=7/2$. Whilst the overall number of sampling points for the multidimensional DFT is still the same, the number of sampling points for the one-dimensional DFT is significantly increased since the ratio of the highest frequency *f*_{1}, to the greatest common divisor of the two frequencies is higher than in the case $f1/f2=3/2$ depicted in Fig. 4(a). Since one has to distribute many sampling points along the *t*^{1} = *t*^{2} line for multiple periods until all harmonics are sampled within a common period, the sampling distribution of the one-dimensional DFT is covering more and more of the whole multidimensional domain depicted in Fig. 4(b). This highlights that the MDFT becomes computationally attractive for cases where the greatest common divisor is small in comparison to the frequencies of interest. But the fundamental advantage of the MDFT is the case when the frequency ratio is irrational since the approach is still capable of resolving both frequencies with the same number of sampling points as in the rational case. In contrast to this, the one-dimensional DFT requires an infinite number of sampling points because the frequencies are not periodic within a common interval and therefore the time vector will never reach its origin on the torus again.

A further crucial advantage of the proposed MDFT method is that the cross-coupling terms between the higher harmonics of each base frequency are directly incorporated within this approach through the computation of the Cartesian product of the frequencies on each axis. In the case of a rational frequency ratio, it is always possible to include the coupling frequencies in a one-dimensional DFT framework as well, possibly at high computational costs if the greatest common divisor is small in comparison to the frequencies of interest. If the frequencies are incommensurable, the definition of appropriate sampling points for a one-dimensional DFT to incorporate the cross-coupling terms is more intricate. To ensure minimal error amplification during the resolution of the system of equations, the condition number of the Fourier transform matrix, which is highly dependent from the chosen instants of time, has to be small. Furthermore, the harmonics of the fundamental frequencies as well as the arising cross-coupling frequencies have to be resolved with an appropriate number of sampling points to ensure small aliasing errors. When using a one-dimensional DFT, the simultaneous consideration of these aspects is far from straightforward and can be avoided by using the MDFT, which directly incorporates the cross-coupling terms.

Furthermore, the MDFT approach is capable of simulating simultaneously two harmonics, which only differ by the interblade phase angle on a single-passage domain. In contrast, the one-dimensional harmonic balance method requires one fundamental interblade phase angle. Therefore, for the one-dimensional approach, it can be necessary to resort to multipassage configurations.

### Underlying Flow Solver.

The multidimensional harmonic balance approach introduced in this work is implemented in the framework of the computational fluid dynamics code TRACE. TRACE is a parallel Navier Stokes flow solver developed at the Institute of Propulsion Technology of the German Aerospace Center (DLR) to model and investigate turbomachinery flows [12]. TRACE has several solver modes for the steady, unsteady, linearized, and adjoint RANS equations. In recent years, it has been extended to solve the harmonic balance equations in an alternating frequency-time-domain approach [9]. The code uses the finite volume discretization of the compressible unsteady Favre- and Reynolds-averaged Navier–Stokes equations in the relative frame of reference in a highly parallelized framework by using a multiblock approach and modern communication structures. Inviscid fluxes are computed using Roe's flux-difference-splitting method. The upwind states are evaluated using the monotonic upstream-centered scheme for conservation laws family of schemes with a modified van Albada limiter to smooth the solution in the vicinity of shocks. Viscous terms are discretized using second-order accurate central differences.

### Application.

To demonstrate the capability of the multidimensional harmonic balance approach to simultaneously capture unsteady disturbances with arbitrary excitation frequencies, the well-established aeroelasticity testcase standard configuration 10 in the presence of an artificial inflow disturbance, which mimics an upstream blade wake, is investigated.

### Standard Configuration 10.

Periodicity | Translational | ||

Pitch | 1 | m | |

Chord | c | 1 | m |

Height | dz | 0.03 | m |

Stagger angle | γ | 45 | deg |

Pitching axis | (0.5, 0.05) | ||

Torsion angle | 0.1 | deg |

Periodicity | Translational | ||

Pitch | 1 | m | |

Chord | c | 1 | m |

Height | dz | 0.03 | m |

Stagger angle | γ | 45 | deg |

Pitching axis | (0.5, 0.05) | ||

Torsion angle | 0.1 | deg |

Inflow Mach number | M_{1} | 0.7 | |

Inflow angle | β_{1} | 55 | deg |

Inflow total pressure | $pt,1$ | 101325 | Pa |

Inflow total temperature | $Tt,1$ | 288.15 | K |

Inflow velocity | v_{1} | 227.387 | m/s |

Inflow pressure | p_{1} | 73047.5 | Pa |

Outflow pressure | p_{2} | 88250 | Pa |

Inflow Mach number | M_{1} | 0.7 | |

Inflow angle | β_{1} | 55 | deg |

Inflow total pressure | $pt,1$ | 101325 | Pa |

Inflow total temperature | $Tt,1$ | 288.15 | K |

Inflow velocity | v_{1} | 227.387 | m/s |

Inflow pressure | p_{1} | 73047.5 | Pa |

Outflow pressure | p_{2} | 88250 | Pa |

*k*=

*0.5 and an interblade phase angle of $\sigma =90deg$ has been simulated. Defining the reduced frequency by*

*k*=

*0.5 is*

where $c\u0302p,\omega $ is computed with the compressible dynamic inlet pressure $pdyn=pt,1\u2212p1$. The standard approach to simulate turbomachinery blade flutter is to predict the unsteady flow in the presence of a prescribed blade vibration. Therefore, it is crucial to accommodate time-dependent boundaries and time-dependent meshes. Detailed information about the implementation and validation of time-dependent meshes in the framework of the harmonic balance solver used in this work can be found in Ref. [14].

### Testcase Extension With Artificial Wake.

*π*. To apply the one-dimensional DFT to the proposed multifrequency testcase, the interval spanned by the sampling points should typically resolve the common base frequency of the combined signals to avoid spectral leakage. Furthermore, the highest harmonic of the common base frequency should be sampled by at least three instants according to Nyquist–Shannon [10] to avoid signal aliasing, but more instants are recommended by Orszag [11]. Since the common frequency of the sinusoidal wake and the frequency of the eigenmode is very low

a factor of $\u223c2.7$ more sampling points than for the MDFT is required to resolve the highest harmonic of the common fundamental frequency. This shows that even in the case where a greatest common divisor exists, simulations with a small greatest common divisor can become inefficient with one-dimensional DFTs and highlight the advantage of the MDFT method. As an alternative to one- or multidimensional DFTs, in industrial applications, the influence of the vibrating blade is often linearized about a previously computed mean flow. In a first step, a computation only with the artificial inflow wake is carried out. Afterward, the time-mean flow field in the presence of the inflow wake is extracted, and the flutter simulation is carried out as a linearized simulation about this new mean flow. This leads to the situation that both base frequencies in the system are completely decoupled in the simulation and no unsteady nonlinear interactions between the harmonics can be resolved. As a further alternative, simulations based on the harmonic set (HS) concept, introduced by the authors, are presented. Detailed information about the implementation and validation of the harmonic set concept can be found in Refs. [9] and [14–17]. To assess the validity and accuracy of the frequency-domain results, all HB simulations are compared to unsteady simulations in the time-domain using the BDF2 scheme with 3018 time steps per segment passing period and 40 subiterations. A summary of all simulations carried out in this work can be found in Table 3. To ensure comparability between all simulation methods, consistent 2D nonreflecting boundary conditions formulated in the frequency-domain have been used for the HB simulations [9,18,19] and for the time-domain simulations [20,21]. The 2D boundary conditions presented in Ref. [18] carry over to the MDFT approach. For the definition of the mean boundary value, the flux averages are defined on the basis of the zeroth Fourier coefficient of the reconstructed fluxes along the boundary.

# | Approach | Base frequency | Harmonics + (extended harmonics) | #Sampling points | #Passages |
---|---|---|---|---|---|

1 | HB MDFT | $feigenmode$ | 0 1 + (2 3) | 65 (273) | 1 |

$fwake$ | 0 1 2 3 + (4 5) | ||||

2 | HB HS | $feigenmode$ | 0 1 + (2 3) | 5 (13) | 1 |

$fwake$ | 0 1 2 3 + (4 5) | 13 (21) | |||

3 | HB DFT | $gcd(feigenmode,fwake)$ | 0 1 7 8 14 15 21 22 29 36 43 | 173 | 4 |

+ (6 9 13 …30 31 35 …66 72 73 …101) | (405) | ||||

4a | HB on mean | $feigenmode$ | 1 + (2 3) | 5 (13) | 1 |

4b | $fwake$ | 0 1 2 3 + (4 5) | 13 (21) | ||

5 | Time-domain | $gcd(feigenmode,fwake)$ | — | — | 4 |

# | Approach | Base frequency | Harmonics + (extended harmonics) | #Sampling points | #Passages |
---|---|---|---|---|---|

1 | HB MDFT | $feigenmode$ | 0 1 + (2 3) | 65 (273) | 1 |

$fwake$ | 0 1 2 3 + (4 5) | ||||

2 | HB HS | $feigenmode$ | 0 1 + (2 3) | 5 (13) | 1 |

$fwake$ | 0 1 2 3 + (4 5) | 13 (21) | |||

3 | HB DFT | $gcd(feigenmode,fwake)$ | 0 1 7 8 14 15 21 22 29 36 43 | 173 | 4 |

+ (6 9 13 …30 31 35 …66 72 73 …101) | (405) | ||||

4a | HB on mean | $feigenmode$ | 1 + (2 3) | 5 (13) | 1 |

4b | $fwake$ | 0 1 2 3 + (4 5) | 13 (21) | ||

5 | Time-domain | $gcd(feigenmode,fwake)$ | — | — | 4 |

For the simulations based on one-dimensional DFTs, the number of sampling points is chosen according to the extended sampling rule proposed by Orszag [11]. In a previous publication by the authors [22], it has been shown that this number of sampling points is recommended to suppress spurious alias terms in the residual functions that are cubic polynomials in the flow fields and its derivatives. For the harmonic set concept, Table 3, case 2, the number of sampling points is computed for each harmonic set separately. This leads to a significant reduction of the required sampling points in comparison to the one-dimensional DFT based on the greatest common divisor of the involved frequencies, Table 3, case 3. For the MDFT method, the number of sampling points on each axis is well-chosen according to Orszag [11]. The overall amount of sampling points is then given by the Cartesian product of the points on each axis, Table 3, case 1. To assess the influence of an increased number of harmonics on the simulation quality, all frequency-domain simulations have been carried out in a basic setup and with an extended number of harmonics.

The artificial wake, prescribed at the entry of the computational domain, has a circumferential wavenumber of $m=\pi $, which means an interblade phase angle of $\sigma =180deg$. Since the eigenmode has an interblade phase angle of $\sigma =90\u2009deg$, the computational domain of the HB simulation based on the greatest common divisor, case 3, and of the time-domain simulation, case 5, have been duplicated to four passages to ensure spatial periodicity. In contrast to this, the MDFT approach is as mentioned previously capable of simulating arbitrary interblade phase angle ratios on a single-passage domain. The harmonic set concept is well able to handle arbitrary interblade phase angle on a single-passage domain, as shown in Ref. [16], but neglects the cross-coupling between the higher harmonics of different harmonic sets, which is included within the MDFT method.

## Results

Figure 6 shows the unsteady imaginary pressure coefficient of the first harmonic of the eigenmodes frequency on the blades surface for harmonic balance simulations based on one-dimensional and multidimensional DFTs resolving the eigenmodes frequency; no artificial wake is prescribed at the entry. The aerodynamic damping for both simulation approaches and different interblade phase angles is depicted in Fig. 7. All presented simulation results show very good agreement to well-known linear reference data from the literatures [23,24].

Figure 8 shows the unsteady imaginary pressure coefficient of the first harmonic of the eigenmodes frequency for different simulation methods in the presence of the artificial wake. At relative chord length $s/c=0.2$, the MDFT simulation resolves a noticeable pressure coefficient drop with a strong gradient on the blade's suction side. The reference harmonic balance computation based on a standard one-dimensional DFT with a small base frequency computed as the greatest common divisor of the harmonics of the wake and the eigenmode confirms the existence of the pressure drop. Both simulations show exact agreement in terms of the pressure coefficient over the relative chord length of the blade, but the simulation based on the standard one-dimensional DFT requires substantially more sampling points than the MDFT, Table 3, which highlights the computational efficiency of the MDFT for configurations with frequencies that have a small greatest common divisor. To verify the general existence of the pressure drop, the results of the unsteady reference simulation carried out in the time-domain are as well depicted in Fig. 8. One can see that the time-domain simulation resolves the pressure drop as well, but with a smaller amplitude than the corresponding frequency-domain simulations. Furthermore, the position of the peak amplitude is shifted downstream toward a higher relative chord length. Despite the differences in the course of the pressure coefficient over the blade, the time-domain simulation generally confirms the existence of the pressure drop. In contrast to this, the simulation based on two harmonic sets and the harmonic balance simulation linearized about the previously computed mean flow in the presence of the wake are not able to resolve the pressure drop. Although a strong nonlinearity is prescribed at the inlet region of the computational domain, the course of the pressure coefficient of the first harmonic over the relative chord shows nearly identical behavior to the linear reference computation without artificial wake [23].

Figure 9(a) shows a snapshot of the instantaneous Mach number distribution over the domain computed with the MDFT method. For clarity, the single-passage frequency-domain simulation has been reconstructed into the time-domain and plotted on multiple passages in a postprocessing step. Furthermore, the amplitude $\alpha =0.1$ of the pitching motion has been amplified by a factor of 100 to emphasize the blades motion. The computed interblade phase angle is $\sigma =90deg$; therefore, the pitching motion of a blade is repeated every fourth blade. The prescribed wake enters the computational domain at the left side of the snapshot. After convecting through the inlet region, the wakes impinge on the suction side of every second blade. At the leading edge of the blades, the interaction of the incoming unsteady disturbance with the surface of the blade generates a local region where the flow state becomes transonic. The additional black contour line indicates the instantaneous $Ma=1$ region, which varies in the translational direction because of the different phase positions of the blades. Although the flow response is expected to be almost linear for the chosen amplitude of the pitching motion, the incoming wakes introduce a time-periodic inflow condition with a strong nonlinear interaction with the moving blade. To verify that the MDFT method correctly simulates the interaction of the incoming wake with the vibrating blade, Fig. 9(b) shows the results of the unsteady reference simulation in the time-domain. The comparison shows good agreement between the MDFT results and the unsteady time-domain simulation. The convection of the incoming wake is captured nearly identically, and the existence of the local transonic flow regimen generated by the interaction of wake and blade is confirmed by the time-domain simulation. Figure 9(c) shows the instantaneous Mach number distribution of the simulation based on the harmonic set concept. The harmonics of the wake and the harmonics of the eigenmode are computed in separate harmonic sets, cf. Table 3, case 2, and a standard DFT is applied to each harmonic set separately. In consequence, the number of sampling points is reduced compared to the MDFT method, case 1, and to the one-dimensional DFT based on the greatest common divisor, case 3, but the cross-coupling effect between the higher harmonics of the incoming disturbance and the flutter harmonic is neglected. One can see that the harmonic set approach is well capable of accurately predicting the convection and interaction of the incoming time-periodic disturbance with the blades. The $Ma=1$ region is computed nearly identical to the MDFT method and to the unsteady simulation in the time-domain. The HB simulation on mean, cases 4a and 4b, computes the interaction of the wake with the blades surface and the interaction of the unsteady flow field in the presence of the eigenmode in two separated simulations. Although the wake can be captured in the simulation 4a properly, the flow state is not transonic in the extracted time-mean solution. Therefore, in the flutter simulation 4b, no interaction between the eigenmode and the $Ma=1$ region can be resolved.

To emphasize how the flow field in the leading region of the blade transitions between a subsonic and a transonic flow state, the maximal Mach number on the surface of one blade for several snapshots in time is depicted in Fig. 10. In most snapshots, the Mach number is in a subsonic region; only in those instants of time when the wake impinge on the blades surface, the Mach number becomes transonic. The periodicity of this transonic region is determined by the time period of the incoming sinusoidal wake. The MDFT simulation, as well as the simulation based on two harmonic sets, captures the periodic behavior very well. Furthermore, both simulations show satisfactory agreement with the unsteady simulation in the time-domain. As mentioned earlier, the HB on mean approach is not able to capture the transonic region in terms of the mean solution in the primary wake simulation, and therefore, the region is not resolved in the subsequent flutter simulation. This is represented by the nearly constant maximal Mach number over the depicted snapshots.

With this knowledge in mind, it becomes clear that the simulation based on the harmonic set concept as well as the HB on mean simulation is not able to accurately predict the pressure coefficient of the first flutter harmonic in the presence of the transonic region generated by the incoming wake, cf. Fig. 8. Within these methods, the first flutter harmonic is not numerically coupled with the harmonics of the incoming wake through the previously mentioned cross-coupling terms. The interactions between the perturbations only occur through their influence on the common harmonics, in this case only through the zeroth harmonic. In consequence, these simulation approaches are not able to resolve the pressure coefficient drop.

To validate the existence of the cross-coupling terms, Fig. 11 shows the amplitude of the complex pressure for the cross-coupling frequency $f=feigenmode+fwake$ simulated with the MDFT method, the HB DFT method based on the greatest common divisor, and the time-domain simulation. Both frequency-domain methods are capable of resolving the coupling frequency, in the case of the HB method based on the greatest common divisor with high computational costs, since the common frequency is low and many sampling points are required to resolve all higher harmonics of interest. Furthermore, it has to be mentioned that in the case where no common divisor exist, the one-dimensional DFT is not able to resolve the sum and differences of the involved frequencies as multiples of a common base frequency. One possible approach to capture the coupling frequencies is to use almost periodic Fourier transforms with equidistant or nonequidistant time sampling. However, as discussed previously, the definition of suitable sampling points that lead to well-conditioned Fourier transform matrices and small aliasing errors is an intricate issue and far from straightforward. This highlights the advantage of the MDFT method to directly incorporate the cross-coupling terms through the computation of the Cartesian product of the harmonics of each base frequency, cf. Eq. (13).

The pressure amplitude depicted in Fig. 11 shows a noticeable influence of the cross-coupling frequency $f=feigenmode+fwake$ to the unsteady flow field. Furthermore, both frequency-domain simulations show very good agreement to the unsteady simulation in the time-domain; the contour lines are computed nearly identical in all depicted simulation methods. This emphasizes the similarity of the computed Fourier coefficients and validates the overall existence and influence of the cross-coupling frequencies. Additionally, one can see, for all simulation methods, that the contribution of the depicted cross-coupling frequency to the unsteady flow field is mainly positioned in the leading region of the blade, the region where the drop of the pressure coefficient in Fig. 8 occurs. This confirms the previous assumption that the pressure coefficient drop is only visible in the simulation methods that are capable of resolving the cross-coupling terms.

To further validate that the existence or nonexistence of the pressure coefficient drop is not an effect of underresolved frequency-domain simulations, Fig. 12 shows the unsteady imaginary pressure coefficient of the first harmonic of the eigenmodes frequency for the extended simulation setups with an increased number of harmonics according to Table 3. One can see that although the number of harmonics is increased, the simulation based on the harmonic set concept is still not capable of resolving the pressure coefficient drop. This emphasizes the restrictions of this approach and highlights that even if harmonic convergence is ensured, important flow physics might be unresolved if the cross-coupling terms are neglected. In contrast to this, the MDFT simulation as well as the HB simulation based on the greatest common divisor still exhibits the pressure coefficient drop. Both simulations show perfect agreement to each other, in case of the DFT based on the greatest common divisor again combined with higher computational costs since more sampling points are required in comparison to the MDFT approach. Furthermore, one can see that increasing the number of harmonics slightly shifts the position of the pressure coefficient drop to a higher relative chord length and in consequence closer to the unsteady simulation in the time-domain. Nevertheless, the peak amplitude of the pressure coefficient in the respective frequency-domain simulations is still smaller than in the time-domain simulation.

To compare the computational costs of the different simulation methods carried out in this work, Fig. 13 shows the CPU hours of the frequency-domain simulations in logarithmic scale for the basic setup as well as the extended setup with an increased number of harmonics. For reference, the CPU hours of the time-domain simulation are as well depicted. In alternating frequency-time-domain approaches, there are two main influencing factors for the computational costs of a simulation: first, the number of sampling points for which the discrete Fourier transforms and its inverse are performed; second, the number of cells in the domain for which the discretized RANS residual $R$ is evaluated, cf. Eq. (10). In time-domain simulations, the computational costs are mainly dependent on the number of time steps per period that are necessary to resolve the time history of the unsteady flow physics and as well on the overall number of cells in the computational domain. The HB simulation based on the harmonic set concept as well as the HB on mean simulation exhibits the smallest number of sampling points, cf. Table 3. Additionally, these approaches are performed on a single-passage domain. Therefore, the computational costs for these approaches are low, Fig. 13, but as shown before, both simulation methods neglect the cross-coupling between the higher harmonics of different fundamental frequencies and in consequence relevant flow physics may be unresolved. Using the extended setup with an increased number of harmonics only leads to a slight increase of the simulation duration for these two approaches since the number of sampling points is not significantly higher than in the basic setup, cf. Table 3.

The simulation based on the greatest common divisor requires the highest number of sampling points. Furthermore, for this specific testcase consisting of an eigenmode with an interblade phase angle of $\sigma =90deg$ and an incoming wake with $\sigma =180deg$, the computational domain has to be duplicated to four passages to ensure spatial periodicity and one fundamental interblade phase angle, which is mandatory for this approach. It has to be mentioned that this is testcase specific and there may be interblade phase angles ratios for which the spatial periodicity can be ensured by less passages or on a single-passage domain. Because of the increased domain in combination with the highest number of sampling points for this testcase, the simulation based on the greatest common divisor exhibits the highest number of CPU hours of all frequency-domain approaches, which becomes even worse for the extended setup. Nevertheless, the computational costs are still an order of magnitude lower than for the time-domain simulation, which highlights that for unsteady periodic flow phenomena, frequency-domain methods are an attractive alternative to the computationally expensive time-domain simulations.

The MDFT method requires a number of sampling points in between the other frequency-domain approaches, cf. Table 3. Furthermore, the approach is always capable of simulating arbitrary interblade phase angle ratios on a single-passage domain. In consequence, the computational costs are an order of magnitude lower than for the HB method based on the greatest common divisor, but still an order of magnitude higher than for the HB method based on the harmonic set concept and more than an order of magnitude higher than for the HB on mean simulation. But in contrast to the latter two approaches, the MDFT method directly incorporates all cross-coupling frequencies, which, otherwise only the expensive simulation based on the greatest common divisor, is able to resolve. Although the true advantage of the MDFT approach is more prominent in multifrequency setups with incommensurable base frequencies, since no other frequency-domain approach considered in this work is able to resolve harmonics of incommensurable frequencies in a straightforward manner, the presented results highlight that the MDFT is well computationally attractive for cases where the greatest common divisor is small in comparison to the frequencies of primary interest.

## Conclusion

A new harmonic balance approach using multidimensional time has been introduced and validated within this work. After a detailed introduction into the theory of the proposed method and its implementation, the method has been successfully applied to a multifrequency testcase.

It has been shown that the interaction of the higher harmonics of the incoming wake with the harmonics of the blades eigenmode leads to a severe drop in the imaginary part of the pressure coefficient, which is proportional to the local damping coefficient. The new harmonic balance approach using multidimensional discrete Fourier transforms is able to capture the pressure drop correctly. Since, in this work, a rational frequency ratio has been chosen, it is possible to compare the results to simulations carried out with one-dimensional discrete Fourier transforms based on the greatest common divisor of the involved frequencies. The results show that this approach is a well capable of accurately predicting the pressure coefficient drop, but for this specific testcase with high computational cost since many sampling points are necessary and the computational domain has to be duplicated to four passages to ensure spatial periodicity. In contrast to this, the MDFT method is capable of simulating arbitrary frequency and interblade phase angle ratios on a single-passage domain and if the greatest common divisor is small in comparison to the frequencies of interest with less sampling points than a standard one-dimensional DFT.

Furthermore, it has been shown that resolving the pressure coefficient drop is strongly related to the ability of the frequency-domain method to incorporate the cross-coupling between the harmonics of the incoming wake and the eigenmodes frequency. It has been validated that the MDFT as well as the much more expensive one-dimensional DFT based on the greatest common divisor is capable of accurately predicting the cross-coupling terms. The comparison to an unsteady simulation in the time-domain leads to good agreement and confirms the overall existence and influence of the cross-coupling frequencies. Additionally, harmonic balance simulations based on the harmonic set concept introduced by the authors and flutter simulations linearized about a previously computed mean flow in the presence of the wake have been carried out. In both methods, the cross-coupling terms are not considered, and therefore the first flutter harmonic is numerically not coupled with the higher harmonics of the incoming wake. In consequence, the drop of the pressure coefficient as a nonlinear interaction effect is not resolved, and the results exhibit an almost linear behavior.

Although the true advantage of the MDFT is more prominent in multifrequency setups with incommensurable base frequencies, it has been shown that the MDFT is well computationally attractive and accurate for cases where the greatest common divisor is small in comparison to the frequencies of the primary interest. In summary, the results highlight that the simulation of multifrequency phenomena such as the interaction of inlet disturbances with vibrating blades requires highly advanced frequency-domain techniques such as the MDFT approach to ensure an accurate prediction of the underlying interaction effects.

## Nomenclature

*c*=chord

- $c\u0302p,\omega $ =
Fourier coefficient of the pressure coefficient

*f*=frequency

- $gcd$ =
greatest common divisor

*i*=complex unit

- $I\u0302$ =
unsteady residual in the frequency-domain

- $Im$ =
imaginary part of a complex quantity

*k*=harmonic index

*K*=highest harmonic

*m*=circumferential wavenumber

*M*=number of time dimensions

- Ma =
Mach number

- $q$ =
vector of conservative flow variables

- $q\u0302$ =
Fourier coefficient of $q$

- $q\u0303$ =
multidimensional reconstruction of $q$

- $R$ =
flow residual

- $Re$ =
real part of a complex quantity

*t*=physical time

- $x,r,\u03d1$ =
cylindrical coordinates

### Greek Symbols

- α =
pitching amplitude

- $F$ =
Fourier transform

- $F\u22121$ =
inverse Fourier transform

- $S$ =
frequency set

- Σ =
interblade phase angle

- $\tau $ =
hypertime vector

- Ω =
angular frequency

### Superscripts and Subscripts

*b*=base frequency

*M*=base frequency index, time dimension index

### Abbreviations

- $DFT$ =
discrete Fourier transform

- $HB$ =
harmonic balance

- $HS$ =
harmonic set

- $MDFT$ =
multidimensional discrete Fourier transform