Interpolation of Transonic Flows Using a Proper Orthogonal Decomposition Method

A proper orthogonal decomposition (POD)method is used to interpolate the flow around an airfoil for variousMach numbers and angles of attack in the transonic regime. POD uses a few numerical simulations, called snapshots, to create eigenfunctions. These eigenfunctions are combined using weighting coefficients to create a new solution for different values of the input parameters. Since POD methods are linear, their interpolation capabilities are quite limited when dealing with flow presenting nonlinearities, such as shocks. In order to improve their performance for cases involving shocks, a new method is proposed using variable fidelity. The main idea is to use POD to interpolate the difference between the CFD solution obtained on two different grids, a coarse one and a fine one. Then, for any new input parameter value, a coarse grid solution is computed using CFD and the POD interpolated difference is added to predict the fine grid solution. This allows some nonlinearities associated with the flow to be introduced. Results for various Mach numbers and angles of attack are compared to full CFD results. The variable fidelity-based POD method shows good improvement over the classical approach.


Introduction
In spite of continuous improvement in computer resources, predicting the aerodynamic characteristics of aircraft using high fidelity computational fluid dynamics (CFD) remains a CPU intensive task.With the confirmed success of CFD in aerodynamic analysis, CFD technology has become a candidate for partially replacing wind tunnel testing in order to obtain aircraft stability and control characteristics.This calls for substantially more CPU capacity, however.For example, as discussed by Salas [1], a CFD database for use in 6-DOF flight simulation requires thousands of data points in the flight envelope.The same type of increase in CPU demand is observed in the aerodynamic and multidisciplinary design optimization communities, where thousands of CFD solutions are often required in typical design optimization studies.A viable alternative to meet this demand is surrogate models.Among the various methods available for constructing surrogate models from a reduced set of high fidelity solutions are reduced-order models which are based on POD.These models constitute a powerful tool for obtaining an approximate CFD solution over the entire domain.Such tools could be used in preliminary design phases because of their easiness of implementation and their low need in computational time.Engineers could easily use these tools to determine a few promising candidates to be tested using high fidelity CFD.
In the context of using POD methods for the interpolation of steady flows, LeGresley and Alonso [2] were the first to use geometric design variables as the varying parameters.More recently, Taeibi-Rahni et al. [21] have investigated the application of POD using Mach numbers and angles of attack as input parameters.However, since POD is a linear decomposition method, problems have been experienced when it is applied to nonlinear flows, such as transonic flows involving shocks.Attempts have been made to solve this problem, one of which was proposed by Lucia et al. [13] International Journal of Aerospace Engineering and is called domain decomposition.Here, the flow domain is divided into multiple zones to isolate the portion of the domain where it is difficult for POD to interpolate.In their case, this was the shock zone.Their remedy is to perform full CFD simulation in the shock zone and use POD interpolation for the rest of the flow domain, which gives very accurate results for the 2D blunt body problem.LeGresley [7] uses the same technique for airfoil analysis, which enables him to predict the position of the shock using only 5 snapshots.The disadvantage of this method is that it requires a very flexible CFD solver and is hardly ever applicable when a commercial CFD solver, such as Fluent, is used.Indeed, such approaches are using the solver to evaluate the weighting coefficients by the mean of residuals minimization.To do so, one must have full control over its solver which is not the case when working with commercial software.Another approach to dealing with transonic flows is to divide the snapshots into two groups, relative to the presence of a shock.The idea is inspired by the snapshot splitting procedure proposed by Cizmas et al. [19] for unsteady flows.Malouin [22] applied this concept to transonic flow, and showed some improvement in the subsonic solutions.In the same vein, Taeibi-Rahni et al. [21] proposed a filtering approach and reported improved results in the transonic regime.
In the present paper, an original combination of POD and the variable fidelity method is proposed to interpolate the transonic flow around an airfoil based on a set of snapshots generated with a commercial CFD solver.The input parameters are the Mach number and the angle of attack.The method is based on the creation of a POD decomposition for the difference between two CFD solutions obtained on two different grids of different resolutions.Then, in order to retrieve an approximation of the fine grid solution, a coarse grid solution is computed using CFD, and the POD method is used to approximate the difference between the grids and generate the fine grid solution.
The next section presents the classical POD method and the concept of our proposed approach based on variable fidelity.Results are presented in Section 3 for the interpolation of the flow field around an RAE2822 airfoil.The results show that the interpolation accuracy in the transonic regime of our new method is an improvement over that of the classical POD method.

POD for CFD Interpolation
In this section, we first briefly review the theory of POD.More details can be found in the work of LeGresley [2,7].Then, we describe our proposed approach, which combines POD and variable fidelity CFD.

Classical POD Theory.
The first task to perform when applying the POD method is to compute the eigenfunctions based on a finite set of CFD solutions.This task is described in the following algorithm.
(1) Obtain  CFD solutions, called snapshots, by varying the Mach number and the angle of attack.
(2) Construct the correlation matrix , as follows: where  is a matrix representing the collection of snapshots for a state variable  ∈ (, V, , , , ) of the flow solution.This matrix is of size  × , where , for a cell-centered scheme, is the number of cells in the mesh and  is the number of CFD solutions.Note that  is symmetric and must be obtained for each state variable.In our case, we repeat the procedure for all state variables, but it can be done at once as described in [23].
(3) Calculate eigenvalues  and eigenvectors  of the correlation matrix  by solving  = . ( (4) Using  and , sort the eigenvalues in descending order and construct  eigenfunctions  for each state variable , as follows: (5) Normalize the eigenfunctions, as follows: Once these steps have been completed, a new solution,  * , can be computed and expressed in terms of the eigenfunctions, as follows: where  is a weighting coefficient chosen to construct the new solution,  is any state variable, and   are the associated eigenfunctions.
In order to obtain the weighting coefficients associated with a new set of flow parameters (Mach number and angle of attack), two approaches have been proposed in the literature.The projection method seeks a least squares minimization of the residuals of a finite-volume formulation of the CFD problem, where the unknowns are the weighting coefficients [7].The interpolated POD method [9] first determines the weighting coefficients associated with the original snapshots and uses simple spline interpolation to estimate these coefficients for a different flow condition.In the present work, we use the interpolated POD approach for its simplicity and for its greater compatibility with commercial CFD software.It is briefly described below.
In order to compute the weighting coefficients associated with the original snapshots, we first rewrite (5) as follows: where    correspond to the weighting coefficient multiplying the eigenfunction ⃗    in order to reconstruct the th snapshot of the state variable .Multiplying ( 6) by    and remembering that the eigenfunctions are orthonormal, the weighting coefficients are obtained as follows: The coefficients    are computed for all the original snapshots, each corresponding to a given value of the parameters (Mach number and angle of attack).A two-dimensional parameter space is then defined to express the    coefficients as a function of the input flow parameters.Then, in order to retrieve the weighting coefficient for a new flow condition, a simple interpolation in the parameter space is performed.In the present work, the weighting coefficients are interpolated with a cubic spline in the two-dimensional parameter space (Mach number and angle of attack).The reader should notice that we use all the POD modes to derive a new solution since their weighting coefficients are simply obtained by a cubic spline interpolation.In other works [10,11,[24][25][26] these coefficients are obtained by optimization, Galerkin projection, and residuals minimization.In this case, the more POD modes you use, the more time consuming the optimization becomes.One has to choose the proper number of modes with respect to the desired level of error to minimize the computational time.This is not the case here, as all modes are used.

Combining POD with Variable Fidelity CFD.
As shown in the Results and Validation section, the classical POD approach described above is not appropriate for predicting the position of the shock in the transonic regime, as its performance is greatly limited by the nonlinearity of the transonic flow.Our approach is to try to reintroduce some information about the nonlinearities in the flow, which we obtain from a coarse grid CFD solution.Our rationale is that coarse grid solutions are much cheaper to obtain than fine grid solutions, but could still provide important information about the flow nonlinearities.This idea originates from [27], where oil flow inside a tank was simulated with two different codes of distinct levels of fidelity.
Furthermore, because we are using commercial CFD software, we do not have full control over it; thus, we cannot use it to evaluate the weightings coefficients.The aim of our paper is to propose a method that can be used with any CFD software without any prerequisite knowledge about the solver itself.In other words, the solver is treated as a black box and the method is thus solver independent.
Our methodology starts by computing two sets of CFD solutions: one set on a coarse grid and the other on a fine grid.The two sets contain the same number of CFD solutions, and these solutions are computed for the same values of the flow parameters.For each flow condition, the difference between the two CFD solutions (Δ) is computed, and the procedure described in Section 2.1 is applied with the Δ solutions as snapshots.This produces a set of eigenfunctions representing the difference between the coarse grid and fine grid solutions.Then, following the interpolated POD approach, weighting coefficients are computed and stored in the two-dimensional parameter space.The computation of a new solution ⃗  Fine Interpolated for different values of the flow parameters follows a two-step procedure.In order to better capture the nonlinearities of the transonic regime, a CFD solution, ⃗  Coarse CFD , is first computed on the coarse grid for the given flow condition.Then, an approximation of the solution on the fine grid is obtained by adding to the coarse grid solution the POD representation of the difference between the grids Δ ⃗  * POD , using weighting coefficients interpolated in the parameter space.This can be described by the following equation: This latter equation requires the evaluation of the coarse solution on the fine mesh.To achieve this task, a simple interpolation of "nearest neighbor" type was performed.
The fact that a CFD solution on a coarse grid is now computed for any new flow condition increases the CPU cost of the interpolation compared to that of the classical POD.However, as discussed above, the CPU cost for a coarse grid solution can be significantly lower than the cost for a fine grid solution.Now, the higher complexity of the algorithm will only be worthwhile if it is associated with an increase in the accuracy of the interpolation.In order to determine whether or not this is the case, we carefully compare the results of classical POD and our proposed approach in the next section.

Results and Validation
3.1.Description of the Problem.The transonic flow over an RAE2822 airfoil, as illustrated in Figure 1, is studied.The transonic regime is emphasized, since our aim is to improve the accuracy of POD interpolation for flow solutions involving shocks.The Mach number and the angle of attack are varied to generate a set of 9 snapshots, as illustrated in Figure 2. In this figure, the circles represent the position in the parameter space where the snapshots were generated, and the square dots represent the positions where the interpolation performance of the two methods will be compared.This Mach/Angle of attack plane was chosen in order to represent a typical transonic range of a commercial plane when flying at cruise conditions.
The coarse grid contains 13600 cells, and the fine grid contains 52650 cells.An overview of the coarse grid, produced with ICEM-CFD, is presented in Figure 3.The fine grid has the same topology and grid concentration.All the snapshots were computed using the commercial CFD flow solver, Fluent 12.1.4.For all the computations, the Reynolds number was 16.5 × 10 6 in average, and + was less than 5 on both grids.The Spalart-Allmaras turbulence model was used, with a 4% turbulence viscosity ratio at the inlet.The densitybased solver combined with the Roe implicit scheme was   selected, and second order upwind resolution associated with the Green-Gauss cell center algorithm was used.Boundary conditions were set to pressure far field.

Validation.
A verification of the implementation is performed first.In fact, as represented by ( 6), the eigenfunctions are theoretically able to perfectly reconstruct any snapshot used to derive them.Figure 4 shows a comparison of the original CFD solution and the reconstructed solution using (6) for a flow condition found in the snapshot set.We can see that the two solutions overlap perfectly, which leads us to conclude that our computer implementation of the POD algorithm is correct and accurate.

Results for the Classical POD Method.
To provide a basis for comparing our proposed method with the classical POD approach, and to better illustrate the difficulties encountered in interpolating flow in the transonic regime using classical POD, we apply this approach, as described in Section 2.1.Referring to Figure 2, the case selected corresponds to Mach number 0.7125 and angle of attack 3.125.The results of the interpolation of the pressure using the classical POD method are reproduced in Figure 5, where they can be compared to that of the full CFD solution for the same flow condition.We can see in Figure 5(b) that two shocks are present in the classical POD interpolated solution, while the CFD results in Figure 5(a) show a single shock.In practice, instead of interpolating the position of the shock, the classical POD method creates two smaller shocks.From a mathematical point of view, this makes sense, because that solution corresponds to a combination of the shocks from all the snapshots.However, from a physics point of view, it is unrealistic to do this.Similar behavior can be observed in ).This inability to interpolate the shock position is the main weakness of the classical POD approach when applied in the transonic regime.In the next section, we apply our new method to that specific situation.

Results
for the Proposed Method.Our proposed approach aims to predict the position of the shock, which classical POD is unable to do, by running a CFD case on a coarse grid and then using POD to interpolate the differences between the solutions obtained on the two different grids, a coarse one and a fine one.The results of our approach are presented in Figures 5(c) to 8(c) and 9, again using the pressure field as the selected variable for comparison.The performance of our new method for the case corresponding to Mach number 0.7125 and angle of attack 3.125 is given in Figures 5(c) and 9(a).In Figure 9(a), we can see that the interpolated solution contains one main shock, properly located when compared to the full CFD solution.In Figure 5(a), we can observe some perturbations on both the upstream and downstream sides of this shock, but they are weak, and the global behavior of the interpolated solution with the new method is a great improvement over that with the classical POD approach.Results for the three other test points defined in Figure 2 ((0.7125; 3.375 ∘ ), (0.7375; 3.125 ∘ ), and (0.7375; 3.375 ∘ )) are presented in Figures 6(c) to 8(c) and 9. Globally, the results of our new method are always better than the classical POD results.On Figures 7 and 8, one can see that, as the Mach number increases, the perturbations located upstream and downstream of the shock are stronger.However, the pressure distribution on the airfoil surface shown in Figure 9 clearly illustrates the superiority of our method over classical POD.The next section will provide a more quantitative measure of the interpolation error for the two methods.From the point of view of CPU requirements, a full CFD solution on the fine grid requires, on average, 1500 iterations and 600 seconds for a 1e-8 convergence level.On the coarse grid, the required CPU time goes down to 150 seconds for the same level of convergence.This represents an acceleration by a factor of 4, considering that the POD interpolation of the differences is almost negligible in terms of CPU effort.

Interpolation Error.
To provide a more quantitative comparison of the two methods, we compare the interpolation errors by computing the difference between the full CFD solution and the interpolated solution, using classical POD and our proposed variable fidelity modification.
Two different error norms are computed, each associated with a different region of the flow field.First, a global error is calculated over all the cells of the grid and is called the GE norm.Second, the error is computed only in the cells adjacent to the airfoil and is called the AE norm.The GE norm provides a good measure of the overall performance of the interpolation method.The AE norm is useful in that it helps us understand the impact of the error on the flow properties close to the airfoil surface, providing a better measure of the error in the flow-structure interaction.For example, the accuracy of aerodynamic coefficients, which are based on surface properties, is closely related to the AE norm.The error norms  are defined as follows: where  is the number of cells in the whole domain for the global error (GE) or around the airfoil for the airfoil error (AE).Table 1 compares the interpolation errors.When comparing the errors between the two methods, we can conclude that the proposed variable fidelity-based POD method always shows improvement over classical POD, since the errors are reduced for all cases.In the cases associated with the lower Mach number, the error reduction is significant, and the errors are reduced by a factor close to 5.However, as discussed in the previous section, the improvement is not as good for the case of the higher Mach number, and error reduction drops to a factor of less than 2.However, the error reduction is greater on the airfoil surface, which is an asset in many situations where surface values are needed.

Conclusion
In this paper, the performance of the classical proper orthogonal decomposition (POD) method is investigated for the interpolation of transonic flows around an airfoil for various Mach numbers and angles of attack.The main weakness of classical POD in the prediction of the shock position in the transonic regime is clearly illustrated, and the causes of this weakness are explained.A new variable fidelity-based POD interpolation method is proposed, which makes use of   CFD solutions on two grids, a coarse one and a fine one.
In this new method, classical POD is used to interpolate the difference between the CFD solutions obtained on these two grids.Then, for any new value of the input parameters, a coarse grid solution is computed using CFD, and the POD interpolated difference is added to predict the fine grid solution.The results show that the new variable fidelity-based POD is always more accurate than classical POD.However, International Journal of Aerospace Engineering the improvement tends to diminish as the strength of the shock increases.Furthermore, there is no limitation for using this method in 3D configurations.The correlation matrix size is only determined by the number of snapshots, and all operations are simple matrix multiplications.In no way the method would become unpracticable because of high number of cells associated with a 3D configuration.
Finally, it would be interesting to compare the performances of the present method with other approaches recently developed.For example, Alonso et al. [10] proposed a new method applied to viscous, transonic aerodynamics.

Figure 3 :
Figure 3: Partial view of the coarse grid used for the computations.

Figure 9 :
Figure 9: Pressure distribution over cord line.

Table 1 :
Comparison of the global error (GE) and airfoil error (AE) on pressure for classical POD (CL-POD) and variable fidelity POD (VF-POD) methods.