Fully-turbulent adjoint method for the unsteady shape optimization of multi-row turbomachinery

The possibility of taking into account unsteady ﬂow effects if performing turbomachinery shape optimization is attractive to accurately address inherently time dependent design problems. The harmonic balance method is an eﬃcient solution for computational ﬂuid dynamics problems of turbomachinery characterized by quasi-periodic ﬂows. If applied in combination with adjoint methods, it enables the possibility to deal with unsteady ﬂuid-dynamic design in a cost effective manner, opening the way towards multi-disciplinary applications. This paper presents the development of a novel fully-turbulent discrete adjoint based on the time domain harmonic balance method and its application to the constrained ﬂuid dynamic optimization of an axial turbine stage. As opposed to previous works, the proposed method does not require any assumption on the turbulent eddy viscosity and on the set of input frequencies. The results show that the method provides accurate gradients, if compared with second order ﬁnite differences, and signiﬁcant deviation with respect to the sensitivity computed with the constant eddy viscosity approximation. The application of the method to the ﬂuid-dynamic shape optimization of the exemplary stage leads to improve the total-to-static eﬃciency of 0 . 8%. The eﬃciency increase is found to be higher than that obtained by means of a steady state optimization method.


Introduction
Adjoint-based shape optimization methods are increasingly becoming essential for automated design. Due to their efficiency in obtaining design sensitivities irrespectively of the number of design variables, these methods have allowed the possibility to effectively tackle multi-disciplinary optimization problems characterized by a high number of design variables and discretized on large domains [1].
Although originally formulated for aircraft design [2,3], adjointbased methods have been successfully extended to turbomachinery design problems. However, the majority of the methods currently adopted is based on steady state flow computations, essentially because this enables the reduction of the computational cost for design optimization [4][5][6][7][8][9].
Given the inherently unsteady nature of turbomachinery flows, the use of unsteady design methods is expected to provide steps forward in fluid dynamic performance as compared to steady state methods. Furthermore, if transient flow effects are accounted for, intrinsically unsteady multi-disciplinary optimization problems of * Corresponding author. turbomachinery can be effectively addressed. Examples include the minimization of tonal noise in transonic fans [10], the minimization of structural excitations caused by dynamic fluid-structure interaction phenomena [11,12], and the aero-thermal performance improvement of cascades subject to unsteady heat transfer mechanisms [13][14][15]. Performing adjoint-based unsteady design for multi-row turbomachinery problems is however a formidable challenge because of i) the presence of multiple interacting blade rows, which results in very costly CFD computations ii) the inherent difficulty of attaining sufficiently converged flow and adjoint solutions, which ultimately affects the accuracy of the calculated gradients.
Due to the high computational cost and memory storage requirements associated with unsteady adjoints [16], several methods have been proposed to improve the efficiency of obtaining time-accurate design sensitivities. These methods mainly target the reduction of memory storage requirements. The algorithm proposed have resulted in less accurate gradient computations by time and space coarsening techniques [17] or in higher I/O overhead if checkpointing algorithms are adopted [18]. Recently, a discreteadjoint method has been applied to time-accurate turbomachinery optimization in combination with time coarsening [19].
Reduced order models have been investigated as a possible effective alternative to time-accurate simulations, in order to de-crease the computational cost and storage requirements of the primal solver. The harmonic balance (HB) method, based on spectral discretization in time of the unsteady flow equations, is a cost effective option for non-linear dynamic problems dominated by a known set of frequencies.
Past work has been conducted to obtain HB-based adjoint design gradients for turbomachinery applications. Nevertheless, these studies are limited to the design of a single blade row, constant eddy viscosity and the inability to solve for spectral gaps, i.e. the inability to deal with frequencies that are not harmonically related [20,21]. A design algorithm based on a fully-turbulent HB adjoint has been recently developed and applied to the optimization of problems characterized by quasi-periodic flows [22]. This algorithm is restricted to a single computational domain, thus only suited for the automated design of isolated turbomachinery cascades.
This paper documents the extension of the novel HB-based design method proposed in Ref. [22] to fully-turbulent multi-row simulations, enabling the solution of quasi-periodic unsteady optimization turbomachinery problems, without any restriction on the turbulent eddy viscosity and on the set of input frequencies to be resolved. The method is based on a duality-preserving approach [23] and it is implemented in the open source code SU2 [24,25].
The design gradients obtained from the HB adjoint equations are verified using second-order central finite differences and applied to the constrained shape optimization of a gas turbine stage. Two expansion ratios are considered for the selected stage, corresponding to subsonic and transonic flow conditions. Furthermore, the baseline stage shape is optimized by means of both the proposed HB-based unsteady method and of a steady state method based on the mixing plane (MP) row interface. The objective of this comparison is to assess whether the HB-based automated design provides a gain in computed fluid dynamic performance over the MP-based one, and if these are dependent on the operating conditions. Finally, the computational performance of the method is analyzed in detail in terms of computational cost, memory, and storage requirements.

Flow solver
Let ρ be the density, E the total specific energy, t time and v the velocity vector in a Cartesian frame of reference, the semidiscrete form of the Navier-Stokes equations can be written as The convective fluxes are and the viscous fluxes are Here, p and T are the static pressure and temperature, κ the thermal conductivity, μ the dynamic viscosity and τ the viscous stress tensor. More in general, for RANS equations, the vector of the conservative variables U can be redefined as and U t is the vector of the conservative variables associated to the selected turbulence model. For example, in case of the Menter Shear Stress Transport (SST) model [27], U t = (ρκ, ρω) with κ the turbulent kinetic energy and ω the specific dissipation.
Using an implicit Euler scheme for time-discretization of (1) leads to where q is the physical time step index, and D t is the timederivative operator. After time-integration and linearizing the residual operator one can obtain the following expression applying the harmonic balance method [22] with dual-time stepping [28] of with U = U q+1 − U q . N is the total number of resolved time instances, linked to the number of input frequencies K by N = 2K + 1. The operator R n is defined as is the harmonic balance operator, calculated as E and E −1 are the direct and inverse Fourier matrix, and D is the diagonal matrix containing the K input frequencies, i.e., D = diag(0, iω 1 , ..., iω K , −iω −K , ..., −iω 1 ). A more detailed description is given in Ref. [29,22]. For a steady-state calculation (7) reduces to

Fully-turbulent discrete adjoint method
The equations for the time-domain HB method formulated in Ref. [22] are here extended to account for multi-row fullyturbulent optimization of turbomachinery. The general formulation given here can be applied to any design problem involving multiple time-zones and geometrical zones and it does not require any restrictive assumption on the eddy viscosity.
The expression of (7) written as a fixed-point iteration is in which U z,n and G z,n are the vector of conservative variables and the iteration operator of the pseudo time-stepping relative to the physical zone z and for time instance n. Each physical zone z corresponds to a blade row.
Using the definition of the operator G z,n given in (12), the optimization problem for the generic objective function J can be written as For each physical zone z and time instance n, α z is the set of design variable, X z,n the physical grid, and M z is a differentiable function representing the mesh deformation algorithm. The objective function J is obtained as the spectral average over the resolved time instances and the physical zones with J * z,n = E * −1 (J z,n E) , (15) in which E * −1 is the extended inverse discrete Fourier transform matrix of size N × N * calculated for N * time instances, whereas E is the N × N discrete Fourier transform matrix computed for the N input time instances (N < N * ). Equation (15) allows one, by means of Fourier interpolation on uniformly spaced samples, to reconstruct the trend of the objective function in time.
The Lagrangian of the constrained optimization problem is given by The differential of the Lagrangian is from which the adjoint equations can be obtained as Equation (19) can be solved directly once the solution of (18) is known. Similar to the flow solver, (18) can be seen as a fixed-point where U * n is the numerical solution for the flow equation (12) and N is the shifted Lagrangian defined as Finally, the gradient of the objective function J with respect to the vector of the design variables α z can be computed from the converged flow and adjoint solutions using The right hand side of (20) is obtained using Algorithmic Differentiation applied to the underlying source code, including the boundary conditions and the stator-rotor sliding-mesh interface. The AD tool adopted [30] makes use of the Jacobi taping method in combination with the Expression Templates feature of C++, leading only to a runtime overhead in the order of 10 − 20% as compared to the flow solver.

Application
The two-dimensional axial stage depicted in Fig. 2 was chosen in this work in order to test the proposed method. The blade geometries correspond to the mid-span profiles adapted from the 1.5 stage experimental setup of the Institute of Jet Propulsion and Turbomachinery at RWTH Aachen [31]. Compared to the original geometry, the stator-rotor blade count ratio has been modified from 36 :41 to 41 :41. In order to resemble the flow characteristics of a typical gas turbine stage, the test case is simulated under the operating conditions given in Table 1, which correspond to subsonic (C1) and transonic (C2) flow conditions.  The fluid dynamic simulations are carried out using the opensource code SU2 [24,25], extended in this work to allow for multirow HB based flow solutions and unsteady constrained optimization using the method discussed in Sec. 2. For both C1 and C2 the Roe scheme [32] is used to discretize the convective fluxes and second order accuracy is obtained by means of the MUSCL [33] approach with gradient limitation based on the Venkatakrishnan limiter. The SST turbulence model [27] is employed for both test cases with a hybrid quad-triangular mesh of approximately 80, 000 elements, in order to ensure a value of y + lower than 1 all along the blade surface. Non-reflective boundary conditions are imposed at the inlet and outlet of the turbine cascade according to the formulation described in Refs. [9,34]. Total conditions and flow direction are imposed at inlet, while the static back-pressure is imposed at outlet. Non-reflectivity treatment is also applied to the stator-rotor interface.
The mesh convergence study for the C2 case using the harmonic balance method with two harmonics is displayed in Fig. 1. A mesh of about 80k elements is deemed suited for optimization purposes and it is then used for both steady and unsteady computations.
The selected objective function is the dimensionless stage entropy generation calculated as in which s s,in and s s,out are the entropy values calculated as mixed-out average [35] over the stator inlet and outlet. The same procedure is used to retrieve the average entropy at the rotor inlet, i.e. s r,in , and at the rotor outlet, i.e. s r,out . T 0s,in is the total temperature at the stator inlet and v 0 the spouting velocity de- where h 0,in is the total enthalpy at the inlet of the stage and h is,out is the isentropic stage outlet static enthalpy. The stage entropy generation is calculated as the summation of the stator and rotor generation separately in order to prevent spurious entropy drops across the interface resulting from numerical accuracy and truncation errors.

Flow field analysis
The results from the harmonic balance (HB) simulations are first verified by comparison with the results obtained using a second-order accurate in time (TA) simulation and those obtained using a steady state mixing plane (MP) model at the stator-rotor interface. The time-accurate simulations are based on the dual time stepping method [28], using 50 time steps per period (corresponding to the blade pitch) and 80 inner iterations for a total of 10 periods.    4. Dimensionless static pressure distribution over the stator and rotor blade surfaces for C1. The total inlet pressure P 0 is used as reference, for both stator and rotor. Furthermore, Table 2 reports the root mean square error of the HBbased and the TA-based stage performance as a function of time (RMSE TA−HB ). Note the RMSE is in the order of 10 −4 , meaning that the two models provide results well in agreement. As a result, two harmonics are used for computing the HB-based adjoint sensitivities described in the following. For both C1 and C2, the total-to-static stage efficiency given by the steady state simulations is characterized by a low deviation compared to the HB timeaveraged results. The relative deviation is approximately 0.05% for C1 and 0.02% for C2. However, the total-to-total efficiency exhibits a larger relative difference between the HB-based and the MP-based simulation results, with the transonic configuration having a relative difference of 0.49% and the subsonic configuration of 0.17%. Finally, Fig. 4 and Fig. 5 show the dimensionless static pressure distribution along the blade profiles retrieved from the MP and HB simulation results. For both C1 and C2, the steady state blade loading differs from the time-averaged harmonic balance blade loading. Furthermore, the shock wave intensity and the location of the associated flow discontinuity computed by the MP-based simulation deviate from the time-averaged HB results, for the transonic configuration C2. The main reasons for this difference are: i) a steady-state model with the MP interface cannot simulate the unsteady potential effects generated by the stator-rotor interaction; ii) the stator wake is not transported to the rotor when using the mixing plane interface; iii) for transonic calculations, the HB method is able to model the unsteady non-linear effects deriving from the shock waves appearing due to the imposed flow conditions and to the time-dependent mutual position of the blade rows.

Adjoint-based design sensitivities
The design gradients of the objective function with respect to the design variables, as defined by (22), are calculated for both the stator and the rotor of the selected test case. To this purpose, two free form deformation (FFD) boxes [36] containing the stator and rotor blade profiles are employed. The design variables (DVs) correspond to the control points of the FFD box as shown in Fig. 6 for a set of twelve DVs per row. Therefore, there are in total 48 degrees of freedom in the optimization problem.
The gradients obtained using the HB-based adjoint equations are first verified using second order finite differences (FD). The results of this verification for the sensitivities in the y-direction are reported in Fig. 7a together with the results obtained using a steady-state adjoint solver [9] with the mixing-plane MP interface.   7. Validation of the C1 normalized adjoint-based design gradients using second-order finite differences, for both MP and HB (Fig. 7a). Rotor design gradients ( y-direction) obtained with the MP and the HB method for a different number of input frequencies (Fig. 7b). The number of the design variables corresponds to the rotor FFD box given in Fig. 6b.
For both HB and MP results there is a very good agreement between adjoint-based and finite differences gradients, with a RMSE of approximately 4e−3. The same results are found for the sensitivities in the axial direction. Fig. 7b depicts the normalized values of the design gradients for an increasing number of resolved frequencies as well as those computed with the steady state MP approach. The design variable numbering corresponds to the DV labels given in Fig. 6b. There are two main observations that can be drawn by analyzing Fig. 7b: i) the value of the gradients computed with the HB-based adjoint is comparatively the same for more than two frequencies; ii) the largest discrepancy between HB and MP simulation results occurs in the proximity of the rotor leading edge. This can be attributed to the effects of wake-rotor interaction, that are not captured by the steady-state model. All gradients converge to the same values towards the rotor outlet.

Constant eddy viscosity (CEV) assumption
As discussed in Sec. 2, the proposed HB-based adjoint method allows one to avoid the use of any restrictive condition on the turbulent eddy viscosity. Past work focused in adopting a constant eddy viscosity (CEV) approximation to ease the development process of the adjoint solver and to improve its computational efficiency but at the cost of a lower gradient accuracy [37,21]. Because of this consideration, the impact of the CEV assumption on the design gradients is assessed. The aim of the analysis described here is to quantify the importance of adopting fully-turbulent adjoint methods for unsteady turbomachinery design. Fig. 8 reports the design gradients of the entropy generation obtained by using a fully-turbulent SST model and the CEV assumption. For both the C1 and C2 operating conditions, the largest differences between the two computed gradients are those relative to the rotor design variables. In addition, the deviations between CEV and SST-based gradients are more marked for the C1 configuration, in case the flow is subsonic. This can be arguably attributed  to the larger share of viscous losses in the subsonic case, and therefore to a higher sensitivity of the objective function to the turbulent flow quantities. Such dependence is somewhat not properly modeled when using a frozen turbulence approximation in the solution of the adjoint equations.
In order to quantify these differences, the relative deviation between the sensitivities computed by the CEV and SST model is presented in Fig. 9. Relative differences in excess of 20% can be observed for the C1 configuration whereas they are up to 12% for the C2 configuration. Differently from what has been reported in Ref. [37], the results of this analysis show that the CEV approximation has a relevant effect on the final design gradients if compared to the fully-turbulent adjoint solution. This outcome confirms the results of a previous similar study, but limited to a comparison between the solutions obtained from a steady-state adjoint method [38].
These findings show that the proposed adjoint method can provide gradients that are significantly more accurate than those obtained by using the CEV approximation with minimal additional runtime cost. This way, it is possible to obtain fast convergence to the actual local optimum with any gradient-based optimization algorithm.

Constrained optimization
The fully-turbulent design sensitivities, computed with the steady state MP and the unsteady HB method, are employed in a gradient-based design procedure to perform a constrained shape optimization of the selected turbomachinery test case. The modified version of the nonlinear least-squares method (SLSQP) [39] was selected as gradient-based optimization algorithm.
The constrained optimization problem of the turbine stage is formulated as minimize α s gen (α) , α = {α 1 , α 2 } subject to:P * = P * 0 , δ t,z = δ t 0 ,z , z = 1, 2 n = 1, 2, ..., N U z,n = G z,n , X z,n = M z,n , (24) in which the objective function is given by the entropy generation of the stage, s gen , averaged over all the N resolved time instances.
with w the Euler work, ṁ the 2D mass flow rate based on the blade pitch y p , ρ 0,in the total density at the stage inlet, and u b the blade speed. Five time instances are selected to perform the optimization, from the analysis of the spectrum of the objective function and from the design sensitivities given by an increasing number of input frequencies (see e.g. Fig. 7b). Fig. 10 shows the evolution of the optimization for both C1 and C2. The entropy generation and the stage power output are scaled in order to better visualize the deviations between the steady and the unsteady results. In the case of the C1 configuration (Fig. 10a)  Fig. 11. Shape optimization for the C1 configuration.  the computed stage entropy generation is reduced by about 7% with the steady state optimization method and by approximately 20% with the HB-based optimization. The constraint on the power output is satisfied within 0.6% in both cases. Fig. 10b depicts the optimization convergence for the C2 configuration. In this case, the objective function is reduced by approximately 11% for the steady optimization and 14% for the unsteady one. The equality constraint on the non-dimensional stage power is maintained within 0.8%.
Figs. 11 and 12 report the baseline and the optimized blade profiles. For both C1 and C2 the largest deformations are located in the area of the rotor leading edge, as consequence of a reduction of the rotor incidence angle. In other words, in both cases, the optimization process succeeded in better aligning the rotor nose to the incoming flow direction. The HB-based optimized shape is furthermore characterized, for C1 only, by a significant shape deformation on the stator suction side and on the rotor rear area. For the stator, this is due by the higher impact of potential stator-rotor interaction effects on the suction side pressure distribution, as can be observed in Fig. 4. From the analysis of the final design, the largest differences between the steady and the unsteady optimization results are associated with the subsonic operative conditions, i.e. C1.
This occurs despite the fact that the discrepancy in performance between the steady and the HB simulation results, computed on the baseline profile, is lower for the C1 configuration (Table 2).
Additionally, a time-accurate simulation is performed using the optimized shapes obtained with both the steady and the unsteady optimization method. Figs. 13 and 14 display the total-to-total stage efficiency and the non-dimensional stage power as a function of time computed with URANS simulations as well as the timeaveraged values. Table 3 reports a summary of the final optimization results based on the total-to-total efficiency. The total-to-total efficiency of the C1 configuration is increased by approximately  one percentage point using the HB-based optimization. The increase in efficiency achieved by means of the steady MP optimization is of 0.2 percentage points. Thus, the unsteady optimization method results in a higher performance improvement. Furthermore, the unsteady-based optimization better satisfies the power constraint, as depicted in Fig. 14a. The final timeaveraged non-dimensional power from the HB-optimized stage differs by 0.3% from the baseline one whereas the MP-optimized solution differs by 2.0%.
The time-accurate simulation results for the C2 configuration indicate that the HB-based optimization leads to an efficiency increase of about 0.4% compared to the 0.2% obtained by the MPbased one. For both final design solutions the power constraint is satisfied within 1%.
Finally, for both operating conditions the optimization leads to a reduced amplitude of the time-dependent efficiency and the dimensionless power. The decrease in amplitude of the total-to-total efficiency is 24.4% for C1 and 14.1% for C2. This demonstrates that the application of the proposed unsteady method intrinsically affects the variation in time of the objective function. As consequence, the variability can be optimized, if needed, by simply reformulating the objective function, i.e. including the amplitude of the quantity of interest in the optimization problem.

Performance assessment
The performance of the proposed HB-based design method is assessed in terms of computational cost, memory and storage requirements. This analysis is conducted for both the primal and the adjoint solver on a 2D and on a 3D geometry. with N the number of input frequencies. This can be explained by recalling that, with the proposed HB method, in order to solve N frequencies 2N + 1 time instances are required as expressed, e.g., in (8)). However, because of the semi-implicit HB formulation adopted in (7), for a number of frequencies higher than 4 a deterioration of the convergence rate is observed. This explains why the computational cost increases at a higher rate if the number of resolved frequencies is greater than four (Fig. 15a). In the case of two input frequencies the HB simulation is approximately 6.5 faster than the time-accurate (TA) simulation. The TA and the HB simulations feature nearly the same CPU time for 10 frequencies.

2D stage
The computational time associated to the steady state simulation is approximately 3 times lower than that of the HB simulation obtained for one input frequency. Fig. 15b depicts the memory and storage requirements for an increasing number of resolved frequencies. The results are normalized using the values obtained from the TA simulations. Both storage and memory increase linearly at a rate of 2N + 1. For 2 input frequencies the memory requirement of the TA simulations is about 4 times lower than the HB simulations, whereas the necessary storage is 41 times higher. These results outline the performance advantage of using HB-based over TA-based adjoint methods, for unsteady turbomachinery design.
The CPU time and memory requirements of the adjoint solver are 1.2 and 4.5 times higher if compared to the primal flow solver.

3D stage
The performance is analyzed for a 3D turbine stage operating at the same working conditions of C1. The goal of this analysis is to assess the capability to obtain adjoint-based sensitivities for a three-dimensional geometry and to evaluate its scalability in terms of computational cost, memory requirements and storage. Given the objective of the analysis, the calculations are conducted by assuming shrouded blades with and free-slip boundary conditions are applied to the hub and shroud. Based on these model assumptions for the problem at hand, a structured mesh of about 300k elements was selected after a mesh independence study. For this test case the numerical schemes employed are those used for C1 and C2.
A HB simulation based on 3 time instances is considered in order to compute the design sensitivities. Fig. 16a shows the geometry of the stage as well as the mid-span contours of the entropy generation normalized with the inlet conditions. The results are relative to the time instance t = T /3. Fig. 16b depicts the HB adjoint-based sensitivity corresponding to t = T /3. Furthermore, the primal and the adjoint solver are tested considering a varying number of input frequencies to investigate the   computational performance. Fig. 17 reports the performance results of the primal flow solver obtained with a number of resolved harmonics ranging from 1 to 5. Similarly to the 2D test case the computational cost, the memory, and storage requirements increase linearly at a rate of 2N +1. When 2 frequencies are resolved, the computational cost and the storage required by the TA simulation are approximately 3.5 and 42 times higher than the HB-based simulation. The memory required by the HB solver for 2 harmonics is about 4.7 times larger than that of the TA solver.
The computational cost of the adjoint solver is approximately 1.2 times higher than the primal solver, whereas the required memory of adjoint solver is about 4.9 times that of the flow solver. The memory required by the adjoint solver for the 250 000 elements mesh was 32.2Gb. The total simulation time for the adjoint solver was of approximately 450 minutes using a 10 cores Intel Xeon E5-2687W v3 CPU with hyper-threading.

Conclusions
This work documents the development of a fully-turbulent harmonic balance (HB) discrete adjoint method for multi-row tur-bomachinery design. The method was applied to the constrained shape optimization of an exemplary axial turbine stage.
The key findings of this study can be summarized as follows 1. The design sensitivities can be accurately calculated with the proposed HB discrete adjoint method. These sensitivities were verified using second order finite differences, without any assumption on the turbulent eddy viscosity. 2. For fluid dynamic design purposes the number of relevant input frequencies to be resolved can be lower than those necessary to accurately model the flow behavior. 3. The assumption of constant eddy viscosity (CEV) was found to significantly affect the accuracy of the design sensitivities.
Relative differences in excess of 20% between the CEV-based and the fully-turbulent results were calculated. 4. Computational cost, memory and storage requirements increase linearly at a rate proportional to the number of time instances. For a number of input frequencies higher than 4 the computational cost featured a slower convergence due to the semi-implicit formulation adopted for the HB flow solver. 5. The HB-based simulations exhibited higher memory requirements but lower storage if compared to time-accurate (TA) RANS simulations. For the analyzed test case, if 2 input frequencies are considered, the memory requirements were approximately 4 times larger than that of the TA simulations, whereas the storage required was about 41 times smaller. 6. The HB adjoint solver featured a computational cost approximately 1.2 higher when compared to the primal flow solver. The ratio between the memory required by the adjoint and flow solver was approximately 4.5. 7. The optimization results achieved by the proposed HB adjoint show remarkable differences when compared with steady state optimization results. Differences in the optimized totalto-total stage efficiency up to 0.8 percentage points were obtained for the exemplary 2D test case. 8. For the analyzed test case, the use of the unsteady optimization method always led to better fluid dynamic performance if compared to the steady state optimization results.
The focus of the present analysis was on the unsteady adjointbased fluid dynamic optimization of a turbine stage characterized by an equal number of blade count per row. Future developments are devoted to extend the periodic boundary condition of the flow solver in order to simulate a single blade passage having unequal azimuthal blade pitch. This would enable the resolution of multistage unsteady problems characterized by a set of frequencies that are not integer multiple of one fundamental harmonic.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.