Simulation of Indexing and Clocking with Harmonic Balance †

The aim of this paper is to demonstrate how the harmonic balance method can be used to predict rotor–rotor and stator–stator interactions in turbomachinery. These interactions occur in the form of clocking and indexing. Whereas clocking refers to the dependency of the performance on the relative circumferential positioning of the rotors or stators, the term indexing is used when different blade (or vane) counts lead to an aperiodic time-averaged flow. The approach developed here is closely related to the one presented by He, Chen, Wells, Li, and Ning, who generalised the Nonlinear Harmonic method to zero-frequency disturbances. In particular, configurations with only one passage per blade row are used for the simulations. We validate the methods by means of the simulation of a fan stage configuration with rotationally asymmetric inlet conditions. It is demonstrated that the harmonic balance solver is able to accurately predict the inhomogeneity of the time-averaged flow field in the stator row. Moreover, the results show that the approach offers a considerable gain in computational efficiency.


Introduction
Since the development of the nonlinear harmonic (NLH) method by He and Ning [1], various frequency domain methods have been presented in the literature that allow for an efficient computational fluid dynamics (CFD) prediction of time-periodic unsteady flows in turbomachinery, in particular the so-called harmonic balance (HB) method [2,3].The fact that the flow residuals in an HB solver are evaluated at discrete time instants means that existing routines from a steady solver can be employed and the development is greatly simplified, even if robust solution methods and accurate blade-row interfaces for HB are still challenging.In the NLH approach, the development of the cross coupling terms between the harmonics can be quite intricate [4].Therefore, only the nonlinear coupling through the zeroth harmonic by deterministic stresses is usually taken into account.The NLH approach, however, has no limitations as to the choice of frequencies.The extension of harmonic balance solvers for multiple base frequencies is the subject of ongoing research [5][6][7].
Similar to the issue of multiple base frequencies is the prediction of rotational aperiodicities in the time-averaged flows, which, in the terminology of frequency domain solvers, can be viewed as additional harmonics with zero frequency and non-zero interblade phase angle.These occur, e.g., when rotor-rotor or stator-stator interactions are to be simulated.He et al. [8] have shown that the NLH approach can be carried over in a straightforward manner to zero-frequency harmonics.As for the harmonic balance method, Subramanian et al. [5] use a small, fictitious shift of the rotation speed and obtain a configuration with several base frequencies.The aim of this paper is to demonstrate that the harmonic set approach presented by the authors [7] is capable of accurately predicting rotationally asymmetric time-averaged flows, without introducing fictitious rotational speeds.
The paper is organised as follows.The harmonic set concept used in the present HB solver is summarised and its generalisation to zero-frequency harmonics is explained.It is shown that the discrete Fourier transform with respect to the passage index leads naturally to a harmonic balance formulation of the possibly aperiodic time-averaged flow.One thereby introduces a generalised notion of sampling points for zero-frequency modes that correspond to circumferential positions rather than time instants.It is shown that, if the sampling points correspond to the true positions of the blades (or vanes), then, in contrast to the method of He et al., the zero-frequency harmonic balance system is derived without using simplifications such as linearisation hypotheses.An academic test case is used to demonstrate the accuracy of the method.Finally, the method is applied to predict the impact of inlet distortions in a fan stage.For this configuration, the results are validated against unsteady simulations in the time domain and HB simulations on full annulus configurations.

The Harmonic Balance Method
The methods described below are integrated into the DLR flow solver TRACE [9,10].TRACE has several solver modes for the steady, unsteady, linearised, and adjoint RANS equations.In recent years, it has been extended to solve the HB equations in an alternating frequency time domain (AFT) framework [7].The underlying spatial discretisations employed here are based on the finite volume approach.We use Roe's upwind scheme [11] in combination with a MUSCL extrapolation to achieve second order accuracy [12].In this work, a van Albada type limiter is employed to smoothen the solutions in the vicinity of shocks [13].
The HB solver is based on the concept of so-called harmonic sets.A harmonic set is defined by a base frequency f , a base interblade phase angle σ and a set of harmonic indices k 1 , . . ., k n .Each harmonic index k j thus corresponds to a harmonic with angular frequency k j ω, ω = 2π f , and interblade phase angle k j σ.There may be several harmonic sets, and different harmonic sets may have some harmonics in common.Keep in mind that a harmonic is determined by both its frequency and its interblade phase angle.Thus, when dealing with two harmonic sets with base frequencies ω 1 , ω 2 and interblade phase angles σ 1 , σ 2 , there may exist integers k 1 , k 2 such that both k 1 ω 1 = k 2 ω 2 and k 1 σ 1 = k 2 σ 2 .Note that, in the solver, these harmonics are then treated as one harmonic flow field.In the following, denote by the i-th harmonic set and the union of all harmonics sets, respectively.The reconstruction in time reads where q (ω,σ) denotes the Fourier coefficient associated with the angular frequency ω and interblade phase angle σ.The harmonics can be defined in two steps.First, the full annulus time harmonic is given by where α 0 = 1 for ω = 0 and α ω = 1 2 for ω > 0.Then, the discrete (complex) Fourier transform with respect to the passage index yields where L denotes the number of passages.The HB method in TRACE is formulated in the frequency domain, i.e., it solves iω q (ω,σ) + R(q for all harmonics (ω, σ).Here, the second term denotes the discrete Fourier transform of the flow residual evaluated at an appropriate number of sampling points.When generalising a single base frequency harmonic balance code to multiple base frequencies, the crucial difficulty is to define the discrete forward Fourier transform.The HB formulation in TRACE uses, for each base frequency f i , equidistant sampling points The number of sampling points N i is set to Here, the number of sampling points per period for the highest harmonic, n hh , must be at least three in order to ensure that the highest harmonic is smaller or equal to the Nyquist frequency.n hh = 4 is the value defined by Orszag [14] to guarantee that taking products of two harmonics in S i does not yield a mode that is indistinguishable from some original harmonic in S i (aliasing).The harmonic balance residual for a mode (ω, σ) with respect to a certain harmonic set is now defined as Here, q S i denotes the unsteady flow field using all harmonics of S i , i.e., the sum in Equation ( 1) is taken over the modes in S i .The multiple harmonic set approach is to approximate the (ω, σ)-harmonic of the residual by a weighted sum where the HB residual for each harmonic set is computed with an individual, appropriate set of equidistant sampling points.The weights w i are assumed to satisfy where χ S i is the characteristic function of the set S i .This means that, in Equation ( 3), a harmonic may be accounted for more than once, but its weights will add up to 1.As explained in [7], the resulting total HB residual is accurate up to second order terms, i.e., the error can be estimated by where R(q) (ω,σ) denotes the exact Fourier coefficient of R(q(t)), viewed as an almost periodic function.
The coupling terms within one harmonic set are even third order accurate, cf.[7].

Example 1: Up-and Downstream Disturbances within a Component
To illustrate the concept of harmonic sets, consider the disturbances caused by up-and downstream stators with N S1 and N S2 blades on a rotor with N R blades.
The disturbances in the rotor can be represented by the harmonics sets which are depicted in Figure 1 by red and violet dots.Here, the circumferential mode order is plotted instead of the interblade phase angle.Keep in mind that, since two angles that differ by a multiple of 2π are identified, mode orders m, m for which m − m is a multiple of N R , are represented by the same mode in S .

Rotor Stator Stator
Figure 1.Red and violet harmonic sets represent the unsteady disturbances caused by two neighbouring blade rows with equal rotational speeds, e.g., wakes and potential effects in a rotor system due to neighbouring stators whose vane count ratio is 2:3.
Note that, since the relative rotational speeds are equal, the disturbances lie on the same line.S 1 and S 2 have a mode in common if there are integers k 1 , k 2 such that the k i -th harmonic is in S i and k 1 N S1 = k 2 N S2 , i.e., k 2 /k 1 equals the vane count ratio.These considerations carry over almost verbatim to a rotor-stator-rotor configuration.

Example 2: Up-and Downstream Disturbances with Different Relative Rotational Speeds
Consider a stator between two rotors with different wheel speeds, e.g., a mid turbine frame consisting of N S struts that are arranged between high-pressure and low-pressure turbines.The disturbances coming from the neighbouring rotors can be represented using harmonic sets with angular base frequencies ω i = |Ω Ri |N Ri and interblade phase angles Observe that the shaft speeds, Ω R1 , Ω R2 , are no longer equal in this case.Typical harmonic sets are depicted in Figure 2.This example is a configuration with a turning strut between high and low pressure turbine blade rows.Here, the high and low pressure components are contra-rotating.
Note that, in Figure 2, a point in the half planes corresponds to a mode rather than to a harmonic.Two modes with different circumferential mode orders m, m are contained in one harmonic (ω, σ) if m∆ϑ ≡ σ ≡ m ∆ϑ mod 2π.
Hence, two modes (ω, m), (ω , m ) in these plots correspond to the same harmonic if ω = ω and m − m is a multiple of N S .Therefore, even for unequal rotor speeds, the direct disturbance harmonic sets may intersect.Furthermore, note that the method of harmonic sets does not resolve possible modes that are generated by the interaction of different harmonic sets.Due to the nonlinear terms in the flow equations, two modes in the plots of Figure 2 will generate modes that correspond to integer linear combinations of the vectors in the half plane.

Strut
Rotor Rotor

Zero-Frequency Harmonics
The aim of this section is to demonstrate that the harmonic set approach can be carried over to the case where a base frequency vanishes.This allows one to simulate rotor-rotor or stator-stator interactions on computational domains that consist of a single passage per row.As an example, consider a blade-row configuration with blade counts N S1 , N R , N S2 and rotor speed Ω R .Then, the transmission and scattering of the wakes of the first stator in the rotor will cause a spectrum of modes (ω k,l , m l ) in the second stator with where k ∈ Z, l ∈ N.These modes can be resolved by harmonics with angular frequency ω l and interblade phase angle In particular, if N S1 is not an integer mutliple of N S2 , there will be modes with zero frequency but non-zero m and thus non-zero σ k,0 .Their interblade phase angles are 2πkN S1 /N S2 , where k is the integer corresponding to the k-th higher harmonic of the original wake.Since, for zero-frequency harmonics, the disturbances with interblade phase angles σ and −σ coincide (their harmonic flow fields are simply complex conjugates of each other), it suffices to consider positive interblade phase angles.Therefore, it follows that one can represent the resulting aperiodicity in Stator 2 by zero-frequency harmonics with interblade phase angle The resulting harmonic set, together with a harmonic set for the downstream disturbance of the rotor, is depicted schematically in Figure 3.For a harmonic set of the form S = (0, 0), (0, σ), (0, 2σ), . .., the reconstructed flow field solves the steady flow equations in all passages if and only if for all passage indices l = 0, . . ., N S2 − 1.For this, it suffices to require Equation ( 5) for the passage indices l = 0, . . ., N S2 / gcd(N S1 , N S2 ) − 1.The crucial observation is that this is equivalent to the standard harmonic balance equation (Equation ( 2)) for ω = 0.Moreover, the notion of a sampling point t j is replaced with the corresponding sampling phase ϕ l .Instead of ϕ j = ωt j for ω > 0, one sets This implies that the flow, reconstructed at the l-th sampling point, Re ∑ k q(0,kσ) e ikϕ l , equals the mean flow in the l-th passage.This approach is therefore exact in that, for the time-averaged flow, it is equivalent to the full annulus nonlinear steady equations.Elementary algebra shows that the ϕ l corresponds to an equidistant distribution of N S2 / gcd(N S1 , N S2 ) points in the interval [0, 2π].To reduce the computational costs, however, one can also approximate the aperiodicity by a harmonic set with K < K max higher harmonics and use N = 2K + 1 equidistant sampling phases ϕ j = 2πj/N.It should be pointed out that the interactions with other harmonic sets are not incorporated in this approach.In particular, the influence of indexing or clocking on the propagation of unsteady disturbances is not captured.

Mode Coupling and Solution Method
The blade row coupling method outlined in [7,15] is based on circumferential Fourier coefficients.In this algorithm, modes with equal frequency but different interblade phase angles can be distinguished.The generalisation to zero frequency modes is therefore straightforward.Similarly, the implicit time-marching procedure to solve the resulting harmonic balance systems is formulated on the basis of decoupled complex linear systems for the harmonics and thus generalises to the case ω = 0.

Numerical Example
The harmonic set approach for zero-frequency harmonics is illustrated by means of a numerical test case consisting of three two-dimensional duct segments through which an artificial jet-like axial velocity disturbance, with amplitude u 0 = 10 m/s is propagated.The duct segments consist of empty cylinder segments and mimic a stator-rotor-stator configuration.The Fourier coefficients of δu can be computed from the Fourier transform of the Gaussian distribution and are given by where m is an integer multiple of 2π/∆ϑ.The unsteady velocity and its reconstruction with up to four harmonics are depicted in Figure 4a.The harmonic balance solver is tested with two harmonic set configurations, summarised in Table 1.The results of the harmonic balance computations are plotted in Figure 4b,c, which show the velocity disturbance at the entry (Row 1), the "rotor" exit (Row 2) and the second "stator" exit (Row 3).The plots for the limited harmonic set (Figure 4b) show that, as expected, the resulting disturbance at the rotor exit corresponds to the reconstruction with two harmonics, whereas, at the second stator exit, it merely consists of the first harmonic.It can be seen that the velocity disturbance is propagated very accurately if it is fully resolved (Figure 4c).It should be pointed out that the harmonic set in the rotor filters the modal content of the disturbance.Hence, it is necessary to resolve the unsteady disturbance accurately in the rotor.(0, 0) (0, 0), (ω, 0), (2ω, 0) (0, 0), (0, σ) Harmonic Set (Conf.2) (0, 0) (0, 0), (ω, 0), . . ., (7ω, 0) (0, 0), (0, σ), . . ., (0, 4σ)

Fan Simulation with Inlet Distortion
In the following, the above methodology is applied to predict the impact of inhomogeneous inflow conditions on a fan stage.The motivation for such a CFD simulation is due to the fact that boundary layer ingesting propulsion systems offer considerable advantages regarding system efficiency and mission fuel burn [16].Whereas in these configurations the rotor is subject to a periodic variation of the inflow condition, the stator has a fixed relative position with respect to the inflow distortion and therefore exhibits indexing.
Here, a modern low pressure fan for a civil aero-engine, derived from a previous DLR design study [17], is used.The fan features a meridional inflow Mach number of 0.65 and a design total pressure ratio of Π t = 1.30.Both rotor and stator are comprised of 15 blades/vanes.The inflow distortion (depicted in Figure 5b) models an ingested boundary layer of an airplane fuselage and is similar to the ones presented by [18,19].It has the form of an ideal turbulent boundary layer that covers 50% of the inlet duct height in the lower part of the annulus.The highest total pressure distortion in this boundary layer is 12% relative to the far-field inflow total pressure.For this numerical study, a setup consisting of three domains is used: inlet duct, fan rotor and fan stator (see Figure 5a).To obtain reference results for the HB solver, a time domain full annulus simulation is performed using the BDF2 time integration scheme with 64 time steps per rotor pitch rotation.Two harmonic balance simulations are run, both using a full annulus inlet duct and a single rotor passage.The setups differ in the stator row where the first one consists of a full annulus while the second one uses a single passage mesh in combination with an additional zero-frequency harmonic set.Harmonic convergence studies have shown that a reasonable agreement between time-integration and HB results can be expected as long as eight harmonics are used for the downstream disturbances.For the upstream disturbances, four harmonics have proven to be sufficient.The meshes used are suitable for wall function boundary layer models and are, up to the number of passages, identical for all setups.The numbers of cells and the computational efforts for the different setups can be found in Table 2. Compared to both time-domain and HB full annulus simulations, the harmonic balance setup with a single passage stator shows a significant reduction in computational time.It should be pointed out, however, that the authors expect further speed-ups of the HB simulations from the use of fast Fourier transformations.In contrast to single passage HB setups, the numerical computation of the circumferential Fourier coefficients at the blade row interfaces can turn into a serious bottleneck when many passages are simulated.
Figure 6a shows the design characteristic with undistorted inflow, together with one operating point under distorted inflow conditions, as predicted by the different simulation methods.
For comparison, results obtained with a steady mixing plane approach with and without inlet distortion are added.Figure 6b,c show a well-known influence mechanism of ingested boundary layers on fan stages: due to reduced inflow velocity in the distorted area, the overall mass flow rate is reduced while the global total pressure ratio rises due to an increased inflow and hence incidence angle [20].For the same reason, the distorted inflow reduces the isentropic efficiency of the fan stage [18].The results show that the decrease in mass flow and the increase in total pressure ratio are predicted correctly by the HB simulations as well as by the mixing plane results with inlet distortion.Moreover, the drop in efficiency predicted by the HB simulations is roughly the same as the one obtained from the unsteady reference results.It is slightly overpredicted by the mixing plane simulations.In order to assess the capability of the HB method to predict local flow phenomena, the specific entropy distributions are compared in a blade-to-blade cut at 90% relative radial height.Figure 7 depicts the part of the slice in which the ingested boundary layer interacts with the fan stage.The distorted inflow appears in the form of increased entropy levels and is deflected in the direction of rotation of the rotor.Furthermore, the increased profile loading in both rotor and stator results in an increased entropy production in the profile boundary layers during the interaction with the distorted inflow.This results in thicker and more intense entropy wakes.These phenomena are correctly reproduced by the HB simulations.The fact that the single stator passage HB simulation does not take into account the impact of the indexing on the unsteady part of the flow in the stator can be seen by comparing Figure 7b,c.There are small differences in the rotor wakes in the middle of the inflow distortion.While in the full annulus simulations (unsteady and HB) the local reduction in axial velocity results in shorter distances between the wakes, this effect is missing in the harmonic balance simulation with a single passage stator.

Conclusions
The harmonic balance approach can be generalised to zero-frequency harmonics in order to resolve rotational aperiodicities of the time-averaged flow in turbomachinery.The harmonic set approach, developed by the authors in recent years, is capable of simulating simultaneously the zero-frequency and unsteady disturbances.While this approach neglects the nonlinear interactions between different harmonic sets, it captures correctly the nonlinear coupling terms within one harmonic set.The approach allows harmonic balance simulations of configurations including two stators or two rotors with different blade counts to be performed using only one passage per blade row.Thereby, computational costs can be reduced considerably.In this paper, it is demonstrated that the approach can also be used to accurately predict the indexing effect of a fan inlet distortion on the mean flow in a fan stage.Comparisons with time-domain and harmonic balance simulations on full annulus configurations show that the time-averaged flow in the stator row as well as the decrease in isentropic efficiency and increase in total pressure ratio are accurately predicted.

Figure 2 .
Figure 2. Green and blue harmonic sets represent the unsteady disturbances caused by two neighbouring blade rows with different rotational speeds whose blade count ratio is 2:3.

Figure 3 .
Figure 3. Two harmonic sets representing the disturbances caused by a rotor and a stator.The blade count ratio of the first two rows is 2:3.

Figure 4 .
Figure 4. Density distribution over circumferential coordinate at inlet, outlet, and third interface.analytic solution (a); numerical solutions with HB configuration 1 (b); and HB configuration 2 (c).

Figure 5 .
Figure 5. Full annulus setup of the fan stage (a) and distorted relative pressure distribution on fan stage inlet plane (a).

Figure 6 .
Figure 6.Compressor map for design speed (a) and impact of the inlet distortion on the operating point.Zoom of total pressure ratio (b) and isentropic efficiency (c) characteristics.

Figure 7 .
Figure 7. Entropy distribution at 90% relative radial height for time domain simulation (a); harmonic balance simulation with full annulus stator (b); and harmonic balance simulation with single stator and zero-frequency harmonic set (c).

Table 1 .
Harmonic set configurations for the numerical experiment of jet propagation through duct segments.

Table 2 .
Numerical setups and computational efforts for time-domain and HB simulations.