Adjoint-based fluid dynamic design optimization in quasi-periodic unsteady flow problems using a harmonic balance method

a duality-preserving approach is proposed. The method is implemented by leveraging algorithmic differentiation and is applied to two test cases: the constrained shape optimization of both a pitching airfoil and a turbine cascade. A key advantage of the current approach is the accurate computation of gradients as compared to second order ﬁnite differences without any approximation in the linearization of the turbulent viscosity. The shape optimization results show signiﬁcant improvements for the selected time-dependent objective functions, demonstrating that design problems involving almost-periodic unsteady ﬂows can be tackled with manageable computational effort.

Shape optimization in unsteady flow problems enables the consideration of dynamic effects on design. The ability to treat unsteady effects is attractive, as it can provide performance gains when compared to steady-state design methods for a variety of applications in which time-varying flows are of paramount importance. This is the case, for example, in turbomachinery or rotorcraft design. Given the high computational cost involved in time-accurate design problems, adjoint-based shape optimization is a promising option. However, efficient sensitivity analysis should also be accompanied by a significant decrease in computational cost for the primal flow solution, as well. Reduced-order models, like those based on the harmonic balance concept, in combination with the calculation of gradients via adjoint methods, are proposed for the efficient solution of a certain class of aerodynamics optimization problems. The harmonic balance method is applicable if the flow is characterized by discrete finite dominant flow frequencies that do not need to be integer multiples of a fundamental harmonic. A fully-turbulent harmonic balance discrete adjoint formulation based on a duality-preserving approach is proposed. The method is implemented by leveraging algorithmic differentiation and is applied to two test cases: the constrained shape optimization of both a pitching airfoil and a turbine cascade. A key advantage of the current approach is the accurate computation of gradients as compared to second order finite differences without any approximation in the linearization of the turbulent viscosity. The shape optimization results show significant improvements for the selected time-dependent objective functions, demonstrating that design problems involving almost-periodic unsteady flows can be tackled with manageable computational effort.

Introduction
The advancement of computational resources has enabled the application of CFD-based design methods to complex shapes, discretized on large domains, often in combination with high fidelity models [1,2].
Optimization methods for design purposes have significantly improved, offering the possibility to deal with complex problems at a reduced computational cost [3][4][5]. In particular, adjoint-based optimization methods [6,7] provide a very efficient approach for computing design sensitivities irrespective of the number of design variables. In applications where the number of design variables is considerably greater than the number of objective functions, adjoint methods allow the computation of optimal solutions in the most cost-effective way, making them well suited for complex industrial applications [8][9][10].
To date, most of the work on adjoint methods has been based on the assumption of steady flow. Obtaining the adjoint solution of an unsteady flow problem can pose a challenge because of the associated large computational and memory requirements [11]. However, accounting for time-dependent flow phenomena in the optimization process is often essential in applications characterized by intrinsic unsteady effects, such as the aerodynamic design of rotorcraft, turbomachinery, open rotors, and wind turbines, to name a few. Furthermore, unsteady adjoint-based optimization methods can pave the way to the solution of multidisciplinary problems characterized by time-dependent phenomena such as those encountered in aeroelasticity or noise reduction [12], for example.
Given the computational cost of accurately obtaining an unsteady solution and thus its adjoint, sufficiently accurate reduced order methods [13] must be used, if the objective is to solve design problems routinely. Among others, the Harmonic Balance method (HB) is a promising option for applications involving quasi-periodic flows characterized by a finite amount of dominant frequencies that need not be harmonics of each other [14][15][16].
Considerable development has been dedicated to HB methods for turbomachinery applications, whereby the flow spectrum is dominated by the blade passing frequencies [15]. One technique that is a subset of the harmonic balance method, sometimes referred to as the Time Spectral method [14,17,15], has been originally formulated for periodic flow problems, and associated adjoints have been derived [18][19][20][21].
HB adjoint-based optimization opens up the possibility to deal with problems featuring flow frequencies that are not harmonically related, thus without the restriction of resolving harmonics of a single fundamental flow frequency. Currently available methods in the literature are limited to inviscid flow problems or frozen turbulence assumptions during design [22,23], which can have a strong impact on the final optimization result [24]. In this work, a fully-turbulent HB adjoint-based shape optimization method is proposed and implemented within the SU2 open-source software environment [25,26]. The algorithm is based on the duality-preserving approach [27][28][29][30], which allows the adjoint solver to inherit the same convergence properties of the primal flow solver. The chosen algorithm ensures robust convergence of the numerical solution of the turbulent adjoint equations.
The method is described in detail first, followed by the illustration of two application cases. An airfoil pitching at a rate characterized by two frequencies that are not harmonically related is first considered, followed by a turbine cascade subject to unsteady inlet conditions.

Time discretization
The semi-discrete form of the Navier-Stokes equations, for a generic cell volume , is (1) is the vector of conservative variables, with ρ the density, v the velocity vector, and E the total specific energy. and its boundary ∂ are assumed to vary their position in time, with velocity u , without deforming. R is the residual operator for the spatial integration of the convective and viscous fluxes, i.e., F c and F v . Using an Arbitrary Lagrangian-Eulerian (ALE) formulation [31], where p is the static pressure. The viscous fluxes are given by where T is the static temperature, κ the thermal conductivity, μ the dynamic viscosity and τ the viscous stress tensor The turbulence modeling is considered, according to the Boussinesq hypothesis, by defining μ = μ l + μ t and κ = κ l + κ t . μ l and μ t are the laminar and turbulent dynamic viscosities, whereas κ l and κ t are the laminar and turbulent thermal conductivities.
The application of time-discretization to (1), using an implicit Euler scheme, leads to where q is the physical time step index, and D t is the time-derivative operator. For a dual time stepping approach [32,33] with pseudo-time τ , one would obtain the following discretization where the additional term is used to relax the solution at each physical time step.

Harmonic balance operator
As given in Ref. [15], the Fourier coefficients resulting from the application of the discrete Fourier transform (DFT) arê , is the vector of the conservative variables evaluated at N time instances Hence the corresponding Fourier coefficients arê with N = 2K + 1, ω k = 2πkf k and K being the number of frequencies. An odd number of time instances is used in this work in order to prevent numerical instabilities [34]. The ensemble of the K resolved frequencies is denoted by ω = [0, ω 1 , ..., ω K , ω −K , ..., ω −1 ]. By defining the DFT matrix as and its inverse discrete Fourier transform (IDFT) one can calculate the Fourier coefficients aŝ with the corresponding vector of conservative variables If the frequencies f k (and hence ω k ) are not multiples of f 1 , once the DFT matrix from (11) or the IDFT matrix from (12) are constructed, it is not possible to obtain an analytical expression for the corresponding inverse matrix. The time operator of (6) can be approximated using spectral interpolation. Applying the spectral operator to the vector of conservative variables Ũ , evaluated at N time instances, yields Given that û is independent from time, using (14) and (13), one can write From (12) where D k,n = iω k δ k,n . (18) D is the diagonal matrix given by One can combine (16) with (17) and define the spectral operator matrix H as Here, E −1 is given analytically by (12), and E is computed by inverting E −1 using Gaussian elimination. It eventually follows that

Time-domain harmonic balance
Considering Ũ , the set of the conservative variables evaluated at N time instances, one can write (7) for a single time Linearization of the residual yields From (21), D t is a linear operator, so the following manipulations are possible: Substitution of the linearized expressions transforms (22) into In this work, a semi-implicit approach is used to solve (26) as whereR Equation (27) is solved for each time instance in a segregated manner. Therefore, an unsteady flow problem characterized by K frequencies requires that the solution of 2K + 1 nonlinear systems of equations need to be computed. The current harmonic balance approach is implemented in the open-source SU2 code [25,26].

Governing equations of the adjoint solver
The adjoint equations are now derived for the proposed Harmonic Balance formulation described in Sec. 2.3 for both laminar and turbulent flows.
Equation (27) is reformulated in terms of a fixed-point iteration for each U n as U q+1 n = G n (U q ) , (29) where G n is the iteration operator of the pseudo time-stepping at time instance n. If G n is contractive, i.e., || ∂G n ∂ U n || < 1, according to the Banach fixed-point theorem [35], (29) admits a unique fixed-point solution U * n such that The aerodynamic design problem, including a possible explicit dependence of the objective function J on the vector of the design variables α, can be expressed as X n are the physical grids constructed for each time instance, and M n (α) is a differentiable function representing the mesh deformation algorithm. The objective function J is obtained as the spectral average (using (42)) over the resolved time The Lagrangian of the constrained optimization problem is given by with λ and μ being the adjoint variables. Given the constraint equations in (31) and omitting in the notation the explicit dependence from the independent variables, one can express the differential of the Lagrangian as from which the adjoint equations can be obtained as and Equation (37) can be solved directly once the solution of (36) is known. Similar to the flow solver, (36) can be seen as a fixed-point iteration in λ n , namely where Û * n is the numerical solution for the flow equation (30) and N is the shifted Lagrangian defined as Therefore, according to the fixed-point theorem, (38) will converge at the same rate as the primal flow solver.
The right hand side of equation (38) is obtained using Algorithmic Differentiation applied to the underlying source code of the program that computes G n . The AD tool adopted [36] makes use of the Jacobi taping method in combination with the Expression Templates feature of C++, leading only to a small runtime overhead.
Finally, the gradient of the objective function J with respect to the vector of the design variables α can be computed from the converged flow and adjoint solutions using

Results
The proposed HB-based adjoint method has been applied to two test cases: the fluid dynamic shape optimization of both a pitching airfoil and a turbine cascade. For both cases, a harmonic balance flow solution has been first obtained and verified against a fully-unsteady simulation using a second-order dual time stepping method [33]. For the spatial discretization, second-order accuracy has been obtained for the convective fluxes using a centered scheme or the MUSCL approach [37] and for the viscous fluxes using a corrected average-of-gradients method. The reader is referred to Ref. [26] for details on the numerical methods and tools available in SU2.
The results obtained with the HB solver at N time instances for a generic quantity of interest are interpolated to a larger time vector t * of length N * using * = E * −1 (E ) , (42) where E * −1 is the larger interpolated IDFT matrix of size N * × N given by E * −1 n,k = e iω k t * n , and E −1 is the N × N IDFT matrix.
The constrained shape optimization problem is solved once the gradient is obtained from the adjoint solution, using a modified version of the nonlinear least-squares method (SLSQP) [38]. Accurate gradient values are obtained by leveraging Algorithmic Differentiation (AD). The reverse mode of the open-source AD tool CoDiPack [29,39] is used in this work to linearize the primal solver.

NACA64A010 pitching airfoil
A two-dimensional pitching NACA 64A010 airfoil in inviscid flow is considered first. The rigid body motion is imposed by assigning a time-varying angle of attack computed as α = 1.01[sin(ωt) + sin(2.6 ωt)] . (43) The flow is transonic with shocks appearing along the upper and lower surfaces of the airfoil while pitching, as can be observed in Figs. 5a, 5b and 5c. The convective fluxes are computed with the JST scheme [32]. The mesh is composed entirely of triangular elements. It contains approximately 11 000 grid points with 200 points along the airfoil and 100 points on the far-field boundary. The main simulation parameters are summarized in Table 1. A reference time-accurate simulation was performed with 160 time steps per smallest period (corresponding to 2.6 ω 1 ) in order to get a well-resolved solution in time. Fig. 1a and 1b show the mesh close to the airfoil and the frequency spectrum of the drag coefficient. Fig. 1b shows that the overall dominant frequencies are the pitching frequencies ω, 2.6ω and their linear combinations. This analysis is performed to select the relevant frequencies for the HB computation.
The test case is then simulated using the harmonic balance method with several choices for the number of time instances with the two input frequencies as pitching frequencies. Fig. 2 shows the convergence of the harmonic balance solution to the fully-unsteady solution with an increase in the number of input frequencies. The converged solution for the selected time instances of both the lift and the drag coefficient is interpolated using (42). According to the spectrum of the drag coefficient (Fig. 1b), the selected input frequencies are reported in Table 2.

Table 2
Optimal time period ratio T opt /T 0 and input frequencies for the pitching airfoil test case, with ω 2 = 2.6ω 1 . T 0 is the period corresponding to the lowest frequency value.
In order to ensure convergence for any set of frequencies, an optimal time pseudo-period is selected using the algorithm proposed in Ref. [40]. A uniform time sampling within the pseudo-period T is adopted, and the values of the optimal time period are reported in Table 2.
Overall, the agreement of the lift coefficient is excellent, resulting in a root-mean-square-error (RMSE) of 0.00275 with just two input frequencies, due to the fact that the movement of the airfoil is primarily in the direction of the lift (Fig. 2a). Conversely, in the case of drag, a larger number of input frequencies is required to capture the fully-unsteady result due to the dominant non-linearities in the flow (shocks). Hence, more than five time instances are necessary for the accurate determination of the drag coefficient, as can be seen in Fig. 2b. Since the selected objective function is the drag coefficient, a suitable number of input frequencies should be identified. The time-averaged drag coefficient is 0.0027 for 5 time instances, 0.0028 for 9, 11 and 13 time instances compared to 0.0028 for the time-accurate unsteady solution. It can be inferred that 9 time instances are sufficient for shape optimization. Fig. 2 reports also the steady-state values of Cl and Cd obtained with After the analysis of the flow solution, the adjoint-based shape optimization problem is considered. The time-averaged drag coefficient c d is minimized over the pseudo-period. Equality constraints on the lift coefficient (c l ) and the maximum thickness (θ max ) are imposed. The resulting optimization problem is set as follows   The gradients of the discrete adjoint are first compared with the gradients obtained with a second-order accurate finite difference (FD) method. This verification is performed over a set of 8 Hicks-Henne bump functions with the first 4 located on the suction side of the airfoil and the others on the pressure side. Fig. 3a depicts the comparison between the adjoint-based gradients and the FD gradients. Fig. 3b shows the convergence of both the flow and the adjoint solver as a function of the number of iterations. As discussed in Sec. 2.4, the convergence rate of the adjoint solver is inherited from the primal flow solver. The computational time ratio between the adjoint solution and the flow solution is about 1.3. The average computational time for one iteration of the primal solver was approximately 0.28 seconds, using the 4 cores Intel Xeon E5-1620 Processor with hyper-threading. Fig. 4a reports the reduction in drag coefficient with the number of optimizer iterations as well as the lift coefficient constraint. The shape optimization process (Fig. 4) achieves a decrease in drag coefficient of about 50%, while maintaining the lift coefficient and the maximum thickness within 1% of the constraint value.
The Mach contours of the baseline airfoil compared with the contours of the optimized airfoil at different time instances are reported in Fig. 5. The optimized airfoil shape leads to an attenuation of the strong shocks characterizing the baseline configuration, resulting in drag reduction.
Finally, the unsteady optimization results are compared with those obtained from a similar constrained steady state optimization of the airfoil at α = 0. The final shape given by the unsteady design method differs from the steady one and  it is shown in Fig. 4b. Furthermore, when the steady-optimized airfoil undergoes the pitching motion prescribed by the test case, the reduction of the time-averaged drag coefficient is 44%, whereas it was reduced by more, i.e. 50%, for the unsteady optimization. This shows that the inclusion of unsteadiness is worthwhile in this case.

T106D-EIZ turbine cascade
The aim of this test case is to assess the capability of the HB-based design method for fully-turbulent flows. In the experimental setup [42], the unsteadiness in the wake of an upstream blade row is approximated by a moving bar, as depicted in Fig. 6a. The moving bars are located at x b /l = 0.7 upstream of the cascade inlet plane, having a velocity v b = 21.4 m/s parallel to the inlet. A schematic representation is shown in Fig. 6a, and the main operating conditions for the test case are reported in Table 3.
The two-dimensional flow domain is discretized with approximately 40 000 elements (Fig. 6b), and the Roe scheme is selected for the convective flux discretization. A suitable spacing of quadrilateral elements is used to cluster the near wall cells such that y + is less than 1. This test case is a benchmark for the study of laminar-to-turbulent transition. However, Table 3 Simulation parameters of the T106D-EIZ test case, as described in [42]. Re 2 is the Reynolds number based on the exit velocity and density; Ma 2 the exit Mach number; T u 1 the inlet turbulence intensity. since the present work aims to assess the methodology for design optimization only, the computations are performed assuming fully-turbulent conditions, employing the SST turbulence model [43]. In order to calculate the cascade performance, the total pressure loss coefficient is evaluated as where P tot,1 ∂ 1 and P tot,2 ∂ 2 are the inlet and outlet total pressure averaged over their corresponding boundary, respectively. P 2 ∂ 2 is the average static pressure at the cascade outlet. The averages at the boundaries are calculated using a mixed-out averaging procedure [44].
First, a validation is performed using the experimental data [42] of the turbine cascade operating at steady state (Fig. 7). The simulation results show a very good agreement with the experimental data. As expected, the main deviation between CFD and experiments occurs at about x/l = 0.75, possibly due to transition not being accurately modeled.
For this test case, two configurations are considered for design: i) a spatially non-uniform, time-dependent inlet boundary condition; ii) a spatially uniform, time-dependent inlet boundary condition. The terminology OptC1 and OptC2 is used to refer to the first and the second shape optimization problem, respectively.

OptC1 configuration
In this configuration, an inlet boundary condition is imposed in order to reproduce the wakes generated by the moving bars. The imposed values of the total pressure, temperature, and flow direction at the boundary are interpolated from the results of a steady state simulation of the flow past the bars. With this boundary condition, only multiples of the blade passing frequency are expected. The fundamental blade passing frequency, given a ratio between the blade pitch and the bar pitch y p /y b = 3, is defined as To verify the HB solution, a second-order time-accurate URANS simulation using the dual time stepping method is performed with a time-step 150x smaller than the lowest period (1/ f 1 ). The total pressure loss coefficient from this simulation is compared in Fig. 8a with the HB solution obtained with 3, 5, and 7 time instances. The selected time instances correspond to the solution for the frequency vectors ω N3 = [0, ± ω 0 ], ω N5 = [0, ± ω 0 , ± 2ω 0 ] and ω N7 = [0, ± ω 0 , ± 2ω 0 , ± 3ω 0 ].
The resolved frequencies are, therefore, multiples of the fundamental blade passing frequency only. The total pressure loss coefficient, defined in (45), and shown in Fig. 8a as function of time, is obtained by spectral interpolation of the harmonic balance results using (42). The RMSE of the total pressure loss coefficient for the solution obtained with 5 time instances is equal to 0.010. The harmonic balance solution obtained with 5 time instances is about 9x faster than the time-accurate  solution calculated over a total simulation time of five periods, which includes the initial transient before reaching convergence to a periodic flow field solution. 5 time instances are used for shape optimization, as a trade-off between accuracy and computational cost.
In Figs. 11a, 11b, and 11c, the Mach number contours from the HB simulation are reported for 3 different time instances with the simulation period given by T = 1/ f 1 . The results show the bar wakes entering the cascade and a separation area occurring at about x/l = 0.7.
Next, the shape optimization problem of the cascade configuration is considered. It can be expressed as minimize α ζ P (U n , X n , α) where the time-averaged total pressure loss coefficient ζ P , obtained from (42), is selected as objective function. Inequality constraints on the absolute exit flow angle (α out ) and trailing edge thickness (δ t ) are imposed. The optimization is performed using an ensemble of 16 geometrical design parameters α based on a free-form deformation (FFD) approach [45]. The gradients of the objective function are again obtained with the proposed adjoint technique and compared with the same gradients obtained by second-order central finite differences (FD). The results of this comparison are reported in Fig. 9a, showing excellent agreement between AD and FD gradients (RMSE = 2 · 10 −5 ). The ratio between the computational time of the adjoint solution and the primal flow solution is approximately 1.7. The average primal solver CPU time per iteration was about 1.41 seconds on a 4 cores Intel Xeon E5-1620 Processor with hyper-threading. Fig. 10 shows that the convergence of the optimization to the minimum objective is nearly reached after only 7 evaluations, although satisfying the constraint requires more evaluations. Fig. 10b highlights that the performance of the optimized  blade is significantly improved, as the total pressure loss coefficient is approximately 38% lower, while the constraint on the absolute outlet flow angle is satisfied. The separation area, as seen in Fig. 11, is considerably smaller with the optimized blade shape. The unsteady optimization leads to a decrease in the peak of the total pressure loss coefficient of 44% and a reduction of 54% of the signal amplitude (Fig. 12a). Furthermore, the objective function spectrum obtained from a URANS simulation of the optimized blade (Fig. 12a) does not contain additional frequencies when compared with the baseline configuration.

OptC2 configuration
In order to investigate the capabilities of the HB-based design method to deal with problems characterized by frequencies that are not multiples of one fundamental harmonic, a second configuration of the T106D-EIZ cascade is considered. In analogy with previous work [46], a time fluctuating inlet total pressure is prescribed as
As for the OptC1 configuration, a time-accurate simulation is performed to verify the HB solution and select the relevant input frequencies by analyzing the spectrum of the total pressure loss coefficient. Fig. 13a depicts the evolution in time of the total pressure loss coefficient resulting from both the URANS and HB computations. Table 4 reports the input frequency vectors for the HB simulation with the associated RMSE between the HB and URANS values of the total pressure loss coefficient.
Figs. 15a, 15b, and 15c show Mach number contour plots, obtained for 3 of the 7 resolved time instances. In this case, as opposed to the OptC1 configuration, the inlet flow field at the upstream boundary is uniform in space with no incoming wakes, given that time-varying, but uniform-in-space, values of total pressure and total temperature are imposed.
After analysis of the flow solution, the shape optimization is performed based on 7 time instances, and the corresponding averaged total pressure loss coefficient is chosen as the objective function. The adjoint-based gradients needed for the optimization are first compared with second-order finite differences. The results of this comparison are reported in Fig. 9b, and are characterized by a RMSE lower than 2 · 10 −5 . Furthermore, as seen in Fig. 13b, the adjoint solver again inherits the convergence rate from the primal solver and its computational cost is about 1.7 higher than that associated with the primal solver. The average computational cost of the primal solver, for a single iteration, was about 1.98 seconds on a 4 cores Intel Xeon E5-1620 Processor with hyper-threading.
The design problem considered in this test case is analogous to that for the OptC1 configuration, and it is formally expressed by (47). Fig. 14a shows that the final design is reached after approximately 14 optimizer evaluations, satisfying the constraint on the outlet flow angle. The Mach contour plots for the optimized shape are shown in Fig. 15. Again, the flow separation on the rear suction side is significantly mitigated for all of the resolved time instances, due to the lower camber angle of the optimized blade profile (Fig. 14b). In this case, a reduction in the average pressure loss coefficient of approximately 14% is achieved with both the amplitude and maximum peak reduced by 44% and 29%, respectively. The time-dependent simulation of the optimized blade confirms these results (Fig. 12b) and reveal that no additional frequencies in the objective function spectrum are present when compared with those appearing in the solution for the baseline shape.

Conclusions
A fully-turbulent harmonic balance discrete adjoint formulation has been developed and applied to the shape optimization of two test cases in unsteady flows: i) a pitching airfoil with a moving grid, and ii) an axial turbine cascade subject to unsteady flow conditions at the inlet. The proposed method is based on a duality-preserving algorithm, which enables the adjoint solver to inherit the convergence properties from the primal flow solver. Due to its efficiency, the framework enables shape optimization for quasi-periodically forced nonlinear fluid problems characterized by a set of frequencies that are not necessarily integer multiple of one fundamental harmonic.
The results of the two test cases have clearly demonstrated that the method is capable of providing accurate gradients in the unsteady setting, as compared to sensitivities computed by second-order finite differences. The gradient-based unsteady optimization has led to improvements of practical significance. The mean and the amplitude of the time-varying aerodynamic losses have been minimized with respect to the baseline configuration. This has been accomplished by considering the minimum number of input frequencies, according to a spectral analysis of the flow field, which results in significant computational cost savings.
The development of a fully-turbulent adjoint optimization framework based on HB paves the way to the solution of additional unsteady design problems that are encountered in numerous advanced applications, such as the multi-disciplinary optimization of turbomachinery, including novel concepts of propulsion systems based on boundary layer ingestion. Future efforts will be devoted to the comparison between fully-turbulent, time-accurate and HB adjoint-based shape optimization methods in order to assess advantages and drawbacks.