Accepted: April 10, 2023 | Published Online: April 12, 2023

Citation: Mifsud M (2023) Modeling the Aerodynamic Load Using a Variable-Fidelity Model Based on Proper Orthogonal Decomposition of a Two-Fidelity Set of Data. Int J Astronaut Aeronautical Eng 8:063

Copyright: © 2023 Mifsud M. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

In this article, a variable-fidelity aerodynamic model is presented based on proper orthogonal decomposition (POD) of an ensemble of computational fluid dynamics (CFD) solutions at different combinations of a set of parameters. The ensemble of CFD solutions consists of two subsets of numerical solutions, or snapshots, that are computed on two different computational meshes. These are inviscid and viscous flow solutions, which are referred to as low-fidelity and high-fidelity data, respectively. In this model, the relatively inexpensive low-fidelity data and the more accurate but expensive high-fidelity data are considered altogether to devise an efficient prediction methodology involving as few high-fidelity analyses as possible while obtaining the desired level of detail and accuracy. The POD of this set of variable-fidelity data produces an optimal linear set of orthogonal basis vectors that best describes the ensemble of numerical solutions. These solutions are projected onto this set of basis vectors to provide a finite set of scalar coefficients that represent either the low-fidelity or high-fidelity solutions. Subsequently, a global response surface is constructed through this set of projection coefficients for each basis vector, which allows predictions to be made at values of parameters not in the original set of observations. This approach is used to predict the aerodynamically distributed load over an airfoil and a complete aircraft configuration. The numerical examples show that the proposed model is efficient and reasonably accurate. Apart from providing a way to integrate variable-fidelity computational data and a complete surface pressure distribution, this technique significantly reduces the inevitable up-front cost of surrogate modeling based on POD.

Keywords

Variable-fidelity modeling, Surrogate modeling, Proper orthogonal decomposition, Radial basis function, Cokriging, CFD

Introduction

Due to the large number of computations that are required to determine the loads imposed on an aircraft during the initial stages of a design process, an efficient computational tool for making rapid but reasonably accurate predictions of the aerodynamic loads acting on an aircraft throughout its flight envelope is sought after. One approach to achieving this is through variable-fidelity modeling (VFM), wherein a combination of several inexpensive to compute low-fidelity analyses and a small number of more accurate but expensive to compute high-fidelity analyses is considered. The basic idea behind VFM is that global trends are captured by the low-fidelity data, while the high-fidelity data is used to correct the global trends, resulting in a response that is close to the high-fidelity response. This combination is commonly used to create efficient prediction approaches that require as few high-fidelity analyses as possible while still achieving the desired level of accuracy.

The variable-fidelity modeling concept has been used by a number of researchers to solve design optimization problems. In their work, various strategies for the reduction of the number of high-fidelity analyses were suggested. For example, Haftka [1] introduced the concept of applying a scaling factor between the variable-fidelity data. In this strategy, a large number of data points are selected for the relatively inexpensive low-fidelity analyses, and from these data points, a subset is chosen where high-fidelity analyses are generated. The low-fidelity results are used to fit a response surface, while the high-fidelity analyses are used to fit a linear correction response surface. For a common design point between the two datasets, the ratio of the responses is evaluated, and this is used for establishing the variable-fidelity approximation. Subsequently, the notion of correction response surfaces was applied by others, such as Hutchinson, et al. [2]. Another strategy presented by Kaufmann, et al. [3] entails the use of low-fidelity models to reduce the region in the design space. Once this is reduced, a high-fidelity response surface is constructed over this reduced space. A slightly different strategy, particularly of interest in optimisation problems, uses a model of low-fidelity to conduct the optimization. Then, occasionally and systematically, one uses information from high-fidelity solutions to check and recalibrate the designs generated. In order to manage the approximations, various ways to decide when the fidelity is increased or decreased were suggested by Alexandrov, et al. [4-8]. An approach related to the one proposed by Haftka uses an additive correction method in which the difference between low- and high-fidelity data is evaluated. This correction is subsequently added to the low-fidelity response surface. A comparison between multiplicative and additive correction surfaces was conducted by Toropov and Markine [9], who suggested that the multiplicative correction leads to better approximations. Another method proposed by Knill, et al. [10] is to use low-fidelity data to identify unimportant response surface terms or to identify insignificant variables that reduce the problem dimensionality.

An autoregressive cokriging approach suggested by Kennedy and O’Hagan [11] was utilized by Huang, et al. [12] and Forrester, et al. [13]. The method has the advantage of being applicable to a dataset with more than two levels of fidelity. In that approach, an autoregressive model was used to calculate the covariances and the cross-covariances in the covariance matrix, and a Bayesian approach was used to predict the output from an expensive high-fidelity simulation code with the assistance of lower-fidelity simulation codes. This Bayesian approach is essentially like cokriging for computer experiments. In the same method, a multiplicative constant scale change and an additive random function were utilized as an adjustment between high- and low-fidelity simulation data. Leary, et al. [14] presented a knowledge-based variable-fidelity approach, where low-fidelity data are dealt with as a priori knowledge in the training process of artificial neural networks and kriging interpolation. This latter approach could achieve a significant reduction in the computational cost with sufficient accuracy. A space-mapping model between low-fidelity and high-fidelity data using POD was suggested by Robinson, et al. [15,16]. This model is based on the gappy POD method for the reconstruction of incomplete datasets developed by Everson and Sirovich [17]. A completely different approach, which utilises radial basis functions to fuse experimental and computational integrated data for a missile configuration, was proposed by Reisenthel, et al. [18] and to fuse surface pressure experimental and computational data for an aircraft wing by Rendall and Allen [19].

In [20] and [21], Han, et al. presented an efficient and accurate aerodynamic data modeling strategy to construct approximation models of aerodynamic force and moment coefficients in a multi-dimensional parameter space based on sampled data. In that approach, a set of CFD methods with varying degrees of fidelity and computational expense was exercised to reduce the number of expensive high-fidelity computations. A low-fidelity CFD method was used to automatically compute hundreds or thousands of solutions at points in the parameter space selected with design of experiments (DOE) techniques. The remainder of the parameter space was "filled in" using interpolation procedures. A few points in the parameter space were selected using DOE methods and computed using a high-fidelity CFD method. Subsequently, the VFM was conducted with low-cost data, indicating trends and a small number of high- cost simulations, correcting the absolute values. This was done using either the so-called bridge functions or cokriging interpolation [22]. The model was then adaptively refined by inserting additional samples based on different refinement criteria. In [23], Keane addresses the problem of robust design optimization using surrogate modelling based on cokriging to avoid sensitive designs that fail to meet the objectives when inevitable uncertainties of manufacture, operating conditions, and degradation are considered. Also, a weighted gradient enhanced kriging is proposed in [24] to address the problem associated with the large cost of training a gradient-enhanced kriging for high-dimensional design optimization problems. Also, Fernández-Godino, et al. discussed in [25] when it is useful to create and use multi-fidelity surrogates and if there are any indications that a substantial cost reduction is realized with multi-fidelity surrogates when compared to surrogate models using high-fidelity samples only and at the same accuracy.

Since industry is interested in very efficiently predicting the surface pressure distribution over an aerodynamic surface, say a complete aircraft configuration, a model that relates the surface pressure at a large number of grid points with respect to a number of independent variables such as the aerodynamic state and control variables is required. Moreover, it is necessary that the model be sufficiently accurate and of low computational cost, both in terms of online and offline costs. The offline computational cost refers to the up-front computational cost necessary to generate the ensemble of observations for VFM. Towards this end, the author embarked on the development of a model by which it is possible to model the pressure distribution Cp on the surface of an aircraft based on proper orthogonal decomposition (POD) of two sets of CFD methods of varying fidelity and hence computational cost, such as the Euler and Navier-Stokes equations. The POD methodology is considered because it provides the aerodynamic model with the desired additional efficiency since, instead of evaluating all of the surface data points, a very small number of coefficients are evaluated for each prediction, hence reducing the number of unknowns. The variable-fidelity ensemble of observations is considered because it further alleviates the up-front computational cost associated with modeling based on POD. In this article, two different variable-fidelity models are proposed. In one of these, the parametric space is augmented by an auxiliary parameter that is binary in nature and is used to denote whether the data is of low-fidelity or high-fidelity. Radial basis functions are used to model the variable-fidelity data while predictions are generated by interrogating the created model for parameter combinations not available in the dataset of observations. Essentially, this method is the same as that presented in [26] and [27], but in this article it is applied to a different flow regime. In the other approach, cokriging interpolation is used to model the cross-correlation between the low- and high-fidelity data and also used to predict the high-fidelity response again at parameter combinations not available in the dataset of observations. In both cases, the low-fidelity data was derived from solutions to the Euler equations, while the high-fidelity data was derived from solutions to the Reynolds Averaged Navier-Stokes equations.

The outline of this article is as follows: In Section 2, a description of the various numerical methods employed is presented. In particular, a description of POD and a detailed description of the two different variable-fidelity models based on POD are made. In Section 3, the results obtained from application of the POD based variable-fidelity aerodynamic model on an RAE 2822 airfoil and a complete aircraft configuration are presented. In Section 4, a conclusion is drawn.

Numerical Methods

In this section, a description of the various numerical methods employed in this work is presented. These include the CFD tool, the POD, which is in effect a dimensionality reduction technique, as well as the variable-fidelity model.

CFD method

In this work, the TAU flow solver [28,29] was used to generate the low- and high-fidelity dataset of observations.

Surrogate modelling based on proper orthogonal decomposition

Generally, surrogate modelling based on POD comprises three main elements; a technique for sampling the parameter space using a design-of-experiment (DoE) approach, the POD of the snapshots set generated, and the fitting of a response surface model through the set of projection coefficients derived from the POD and for each POD mode. Below is a detailed description of the latter two elements.

Proper orthogonal decomposition: In fluid mechanics, POD was first introduced by Lumley [30] in the context of stochastic turbulence. The same procedure is commonly referred to as Karhunen-Loéve expansion and principal component analysis. POD is also very closely related to the singular value decomposition. The POD provides a basis for the modal decomposition of a dataset, usually obtained from experiments or numerical simulations. The resulting basis vectors are called proper orthogonal modes and are the best possible uncorrelated and data-dependent linear set of basis vectors that describe the initial observations.

The POD method has been used extensively in various fields. In the field of fluid dynamics, it has been used in unsteady flow problems such as aeroelasticity and stochastic turbulence to capture the temporal variation of the flow. It was also used in steady flow problems to capture parametric variations [31-34]. The POD is usually described using the calculus of variations applied to a multi-dimensional spatio-temporal dataset. In such problems, the dataset is produced from solution vectors obtained at a particular instant of time. In this article, the POD is described for steady flow problems in terms of the singular value decomposition (SVD). Although these two approaches are equivalent, the SVD approach is preferred as it is more straightforward.

Let $A\in R^{mxn}$ denote the matrix whose rows are the snapshots with data centered about the origin, where m is the number of snapshots and n is the number of grid points. Each snapshot is constructed by placing in order the solution at each grid point for the complete grid. This order can be determined arbitrarily, but it must be consistent throughout the whole set of snapshots.

$A=\left(\begin{array}{ccccc}{\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_{\text{1}}{}^{\text{(1)}}& \text{.}& \text{.}& \text{.}& {\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_n{}^{\text{(1)}}\\ {\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_{\text{1}}{}^{\text{(2)}}& \text{.}& \text{.}& \text{.}& {\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_n{}^{\text{(2)}}\\ {\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_{\text{1}}{}^{\text{(3)}}& \text{.}& \text{.}& \text{.}& {\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_n{}^{\text{(3)}}\\ \text{.}& \text{.}& \text{.}& \text{.}& \text{.}\\ {\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_{\text{1}}{}^{\text{(}m\text{)}}& \text{.}& \text{.}& \text{.}& {\text{<mover accent="true"> a <mo>˜</mo> </mover>}}_n{}^{\text{(}m\text{)}}\end{array}\right)\text{ (1)}$

The SVD of A can be written as $U\sum V^T$ where $U\in R^{mxn}$ and $V\in R^{nxn}$ orthogonal matrices and these feature the left and right singular vectors as columns, respectively. $\sum \in R^{mxn}$ is a diagonal matrix whose diagonal elements consist of $q=\min \left(m,n\right)$ non-negative real numbers ${\sigma }_i$ arranged in decreasing order, that is, ${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3\ge \mathrm{.....}\ge {\sigma }_q$ where ${\sigma }_i$ are referred to as the singular values of A. Since Σ is a diagonal m x n matrix, then the above matrix decomposition can be written in reduced or thin form as follows if we assume that m < n,

$$A^{mxn}=U^{mxn}{\sum }^{mxn}{V}^{Tmxn}\text{ (2)}$$

where the matrices Σ and VT are reduced in size. The columns of V and hence the rows of VT are the proper orthogonal modes of the system. Writing the product of U and Σ as a matrix $\left[{\alpha }_{ij}\right]$, Eq. (2) can be written as ${\tilde{a}}_k^{(i)}={\displaystyle \sum _{j=1}^m{\alpha }_{ij}{v}_{jk}}$. The scalar coefficients αij are also referred to as projection coefficients because these are obtained by projecting the CFD solution onto the basis vectors. A complete reconstruction of the snapshots can be obtained from

$$a=〈a〉+\tilde{a}\text{ (3)}$$

Where $〈a〉$ is the mean vector by which the snapshots were centered. Now, a may represent a vector of scalar physical properties such as the primitive or conservative variables and therefore the method described can be applied to each variable in turn to form a distinct basis for each variable. However, an improvement in the ability of the basis to represent the system may be achieved by considering not only how the individual variables vary from one snapshot to another but also how variables change relative to one another. Hence a is developed from state variable vectors consisting of all the primitive or conservative variables simultaneously [32]. In this case, the POD modes are sensitive to the scaling of the flow variables as these are in different units and have significantly varying magnitudes. Consequently, appropriate scaling factors are necessary for each fluctuating flow variable which makes their magnitude of the same order [35].

When a problem is represented by a number of snapshots from which a suitably rich set of basis vectors is available, the singular values rapidly become small and a few basis vectors are adequate to reconstruct and approximate the snapshots. Assuming that p modes corresponding to the p largest singular values are dominant, then the relative energy E or variance in the data captured by the first p modes can be computed as

$$E\left(p\right)=\frac{{\displaystyle \sum _{i=1}^{i=p}{\sigma }_i^2}}{{\displaystyle \sum _{i=1}^{i=m}{\sigma }_i^2}}\text{ (4)}$$

If this relative energy is sufficiently high enough, say over 99.99% of the total relative energy, then p modes are adequate to capture the principal features and approximately reconstruct the dataset. Thus, a reduced subspace is formed which is only spanned by p modes.

Fitting a response surface model and snapshots prediction: To predict the projection coefficients for snapshots at intermediate parametric values not included in the original data ensemble, a response surface model can be fitted through the set of projection coefficients for each basis vector, as long as ${\alpha }_i$ varies as a smooth function with the change in parameters. In this work, radial basis functions namely thin-plate splines [36] and cokriging [37] that produce an interpolative fit through all of the sample points were considered.

The predicted solution vector $a\left(\beta \right)$ for any variable β within the parametric space is then given by,

$$a\left(\beta \right)=〈a〉+{\displaystyle \sum _{i=1}^{p^{\prime }}{\alpha }_i\left(\beta \right)}{v}_i\text{ (5)}$$

where p’ is normally greater than p and the weighting coefficients ${\alpha }_i\left(\beta \right)$ are found by interrogating the response surface.

The variable-fidelity surrogate model based on POD: Consider a two-fidelity dataset comprising of low-fidelity and high-fidelity CFD solutions. If these low-fidelity and high-fidelity solutions significantly correlate with one another, a variable-fidelity model based on POD can possibly be set up. In this article, two VFM methods based on POD are presented. In both methods, proper orthogonal decomposition is performed on the set of observations. In one method, the parametric space is augmented with an auxiliary parameter that is binary in nature and denotes whether the snapshot is either of low- or high-fidelity. From the resulting POD scalar coefficients, a multi-dimensional global response surface is generated, which is subsequently interrogated to generate solutions at parameter values not available in the dataset of observations. The second method is that of applying the cokriging technique to the low- and high-fidelity sets of POD coefficients from which high-fidelity coefficient predictions are generated. In contrast with the previous methodology, this technique inherently models the cross-correlation between the low- and high-fidelity POD coefficients rather than augments the parametric space using an additional variable. In both cases, the predicted coefficients are subsequently used together with the POD modes to obtain the required high-fidelity solution. In the next two sections, a detailed description of these two methods is presented.

Augmenting the parametric space with an auxiliary parameter: Considering an ordered ensemble of variable-fidelity data A, where $A\in R^{mxn}$ is obtained from solution vectors of the Euler and the Reynolds averaged Navier-Stokes equations at various parameter values with the total number of realizations or parameters combination $m=m_1+m_2$, where $m_1$ is the number of snapshots obtained from the high-fidelity CFD method and $m_2$ is the number of snapshots obtained from the low-fidelity CFD method. It is assumed that $m_1<m{}_2$ and n is the number of grid points over which the computational domain is sampled.

It is assumed that the rows of the matrix A are centered snapshots where the primed entries denote the high-fidelity solution vectors, from row 1 to row $m_1$. The non-primed entries represent the low-fidelity solution vectors, from row $m_1+1$ to m. In this variable-fidelity model, the parameter combinations between the low and high-fidelity CFD solutions can be different and it is not necessary to have common snapshots between the variable-fidelities.

$$A=\left(\begin{array}{l}{a^{\prime }}_1^{(1)}...{a^{\prime }}_n^{(1)}\\ .....\\ {a^{\prime }}_1^{(m_1)}...{a^{\prime }}_n^{(m_1)}\\ a_1^{(m_1+1)}...a_n^{(m_1+1)}\\ .....\\ a_1^{(m_1+m_2)}...a_n^{(m_1+m_2)}\end{array}\right)\text{ (6)}$$

By evaluating the SVD of A in reduced form and multiplying the left hand singular vectors with the singular values of A to obtain the projection coefficient matrix $\left[{\alpha }_{ij}\right]$ then,

$$A=\left(\begin{array}{l}{{\alpha }^{\prime }}_{11}.{{\alpha }^{\prime }}_{1m_1}{{\alpha }^{\prime }}_{1\left(m_1+1\right)}.{{\alpha }^{\prime }}_{1m}\\ ......\\ {{\alpha }^{\prime }}_{m_{1}1}.{{\alpha }^{\prime }}_{m_1{m}_1}{{\alpha }^{\prime }}_{m_1\left(m_1+1\right)}.{{\alpha }^{\prime }}_{m_{1}m}\\ {\alpha }_{\left(m_1+1\right)1}.{\alpha }_{\left(m_1+1\right)m_1}{\alpha }_{\left(m_1+1\right)\left(m_1+1\right)}.{\alpha }_{\left(m_1+1\right)m}\\ ......\\ {\alpha }_{m1}.{\alpha }_{m{m}_1}{\alpha }_{m_1\left(m_1+1\right)}.{\alpha }_{mm}\end{array}\right)\cdot \left(\begin{array}{l}{V}_{11}.V_{1n}\\ ...\\ V_{m_{1}1}.V_{m_{1}n}\\ V_{\left(m_1+1\right)1}.V_{\left(m_1+1\right)n}\\ ...\\ V_{m1}.V_{mn}\end{array}\right)\text{ (7)}$$

Note that the rows of the RHS matrix (equivalent to $V^T$ ) in Eq. (7) represent the proper orthogonal modes of the system. From this matrix, the set of scalar coefficients ${{\alpha }^{\prime }}_{i,j=const}$ and ${\alpha }_{i,j=const}$ are considered together to develop multi-dimensional response surfaces for variable-fidelity data representation. Note that the high-fidelity response surface is formed by the projection coefficients ${{\alpha }^{\prime }}_{i,j=const}$ for $1\le i\le m_1$ while the low-fidelity response surface is formed by the projection coefficients ${\alpha }_{i,j=const}$ for $m_1+1\le i\le m_1$. This is achieved by the introduction of an auxiliary variable $\epsilon \equiv x_{d+1}$ to the d-dimensional problem with variables $\left(x_1,x_2,\mathrm{......},x_d\right)$. This auxiliary variable simply denotes whether the data is of low-fidelity $\left(\epsilon =0\right)$ or high-fidelity $\left(\epsilon =1\right)$. From this, a global response surface is then computed in the d + 1 dimensions using thin-plate splines. By interrogating the newly developed response surface projected along $\epsilon =1$ at any arbitrary parameter value not in the original set of observations, predictions can be made. The resulting model representation respects the accuracy of the high-fidelity data while following the features of the low-fidelity data. In this model, the introduction of the low-fidelity sub-space by making use of an additional variable, allows an interpolation-based extrapolation to be performed. In other words, an extrapolation at high-fidelity based on interpolation at low-fidelity is conducted.

Following this, the predicted solution vector $\widehat{a}$ is determined by

$$\widehat{a}=〈a〉+{\displaystyle {\sum }_{k=1}^q{{\alpha }^{\prime }}_k{v}_k}\text{ (8)}$$

Where q < m and $〈a〉$ is the mean vector. Since m1 is normally a small number, all of the high-fidelity POD modes are considered. This model offers the advantage that it gives access to the full flow-field data in contrast with other variable-fidelity models reported in the literature which deal with integrated quantities only.

In this work, the values of the POD coefficients were normalised between 0 and 1 before conducting the interpolation process because these projected coefficients varied significantly between the low- and high-fidelity data. This normalization procedure provided an improved match between the predicted solutions and the high-fidelity solutions at the same parameter values, when a comparison was made. If ${\alpha }_i^p$ is the ith projection coefficient along POD mode p of the two-fidelity set of snapshots, where i varies between 1 and the number of used modes, then all of the coefficients along the POD mode p were normalized using

$${\alpha }_i^p=\left({\alpha }_i^p-{\alpha }_{\min }^p\right)/\left({\alpha }_{\max }^p-{\alpha }_{\min }^p\right)\text{ (9)}$$

where ${\alpha }_{\min }^p$ is the minimum projection coefficient value amongst the set ${\alpha }_i^p$ and ${\alpha }_{\max }^p$ is the maximum projection coefficient value amongst the set ${\alpha }_i^p$. Ultimately, the interpolated scaled coefficients are de-scaled before generating the predicted solution. Another technique by which cross-correlation can be modelled is through the use of cokriging, which is the theme of the next section.

Cokriging model: Both kriging and cokriging are statistical interpolation methods for determining the so-called best linear unbiased estimator for a given set of sample points with corresponding sample values. The kriging estimator depends on information about the spatial correlation between the given sample points. For cokriging interpolation, the correlations and cross correlations between both the sample points for a primary variable and those for one or more secondary variables are considered. The estimator is modeled as a random function, and the spatial correlations are determined by the co-variances between the function values at different sample points, which are then thought of as realizations of the random function at those specific points. In practice, the co-variances are unknown and therefore must be approximated. In most applications, the correlations are estimated via standard functions, such as the Gaussian or cubic spline correlation functions, to model spatial dependencies. In general, these functions are weighted distance functions, and optimization of the model is conducted to choose the distance weights such that the statistical model likelihood is maximized. Once the appropriate best linear unbiased estimator model is chosen, the optimization procedure is analogous to kriging and cokriging even when more sophisticated regression models are considered. Various cokriging models were considered, such as in [11,13,22]. The implementation follows the one in [37].

Let $y:R^d\to R,x\to y\left(x\right)$ be a function representing the high-fidelity POD projection coefficients associated with a particular POD mode. Also let $z:R^d\to R,x\to y\left(x\right)$ be a function representing the low-fidelity POD projection coefficients associated with the same POD mode. Hence, z is a low fidelity representation of y.

Suppose that and are sampled at certain locations in the parametric space $x^i\in R^d,i\in \left\{1,\mathrm{...}\mathrm{...}\mathrm{...},m_1\right\}$ and ${x^{\prime }}^j\in R^d,j\in \left\{1,\mathrm{...}\mathrm{...}\mathrm{...},m_2\right\}$ respectively and with corresponding responses $y_i=y\left(x^i\right)\in R,i\in \left\{1,\mathrm{...}\mathrm{...},m_1\right\}$, and $z_i=z\left({x^{\prime }}^j\right)\in R,j\in \left\{1,\mathrm{...}\mathrm{...},m_2\right\}$, obtained from the projection matrix of Eq. (7). Considering the POD coefficients associated with the jth POD mode in a vector as follows

$$\left(Y^T,Z^T\right)={\left({{\alpha }^{\prime }}_{1j},\mathrm{...},{{\alpha }^{\prime }}_{m1j},{\alpha }_{1j},\mathrm{...},{\alpha }_{m2j}\right)}^T\in R^{\left(m_1+m_2\right)}\text{ (10)}$$

Let $1_{m1}={\left(1,\mathrm{...},1\right)}^T\in R^{m_1}$ and $1_{m2}={\left(1,\mathrm{...},1\right)}^T\in R^{m_2}$ be the vectors with all entries equal to one of suitable dimensions and let

$$F=\left(\begin{array}{l}{1}_{m_1}0\\ 01_{m_2}\end{array}\right)\in R^{\left(\left(m_1+m_2\right)\times \left(m_1+m_2\right)\right)}$$

The cokriging best linear unbiased predictor becomes

$$\widehat{y}\left(x\right)=c^T\left(x\right)C^{-1}\left(\left(\begin{array}{l}Y\\ Z\end{array}\right)-F\left(\begin{array}{l}{\beta }^y\\ {\beta }^z\end{array}\right)\right)+{\beta }^y\text{ (11)}$$

Where c(x) is the auto- and cross-variance vector and C is the auto- and cross-covariance matrix and ${\beta }^y,{\beta }^z$ are the regression constants, all of which are explained below. Covariances are modeled by spatial correlation functions of the form

$$r\left(\theta ,x^i,x^j\right)={\prod }_{k=1}^{d}sc{f}_k\left(\theta ,x^i,x^j\right)\text{ (12)}$$

Where $\theta =\left({\theta }_1,\mathrm{...}\mathrm{...},{\theta }_d\right)\in R^d$ is a vector of hyper-parameters called distance weights. For example, the Gaussian correlation is

$$sc{f}_k\left(\theta ,x^i,x^j\right)=\exp \left(-{\theta }_k{\left|x_k^i-x_k^j\right|}^2\right),r\left(\theta ,x^i,x^k\right)=\exp \left(-{\displaystyle \sum _k{\theta }_k}{\left|x_k^i-x_k^j\right|}^2\right)$$

Following the simplified covariance fitting scheme proposed in [32], the auto- and cross-covariance matrix is modelled using

$$C\left(\theta ,\gamma \right)={\sigma }^2\left(\begin{array}{l}r{\left(\theta ,x^i,x^j\right)}_{ij}yr{\left(\theta ,x^i,{x^{\prime }}^k\right)}_{ik}\\ \gamma r{\left(\theta ,{x^{\prime }}^k,x^j\right)}_{kj}r{\left(\theta ,{x^{\prime }}^k,{x^{\prime }}^l\right)}_{kl}\end{array}\right)={\sigma }^2\rho \left(\theta ,\gamma \right)\text{ (13)}$$

for i, j = 1 ,...., m1, k , l = 1 ,......, m2. Hence, the covariance matrix C and the correlation matrix ρ are symmetric matrices of dimensions (m1 + m2) by (m1 + m2). The covariance vector at an untried location x is

$$c\left(x\right)={\sigma }^2\left(\begin{array}{l}r\left(\theta ,x^i,x\right)\\ \gamma r\left(\theta ,{x^{\prime }}^k,x\right)\end{array}\right)\in R^{m_1+m_2},i=1,\mathrm{...}\mathrm{...}\mathrm{..},m_1;k=1,\mathrm{...}\mathrm{...}\mathrm{...},m2\text{ (14)}$$

The value ${\sigma }^2$ denotes the simplified process variance and $0<\gamma <1$ is an additional hyper-parameter, which controls the correlation between the primary and secondary variables, for details, see [37]. The simplified process variance depends on the hyper-parameters as follows

$${\sigma }^2\left(\theta ,\gamma \right)=\frac{1}{m}{\left(\begin{array}{l}Y-1_{m_1}{\beta }^y(\theta ,\gamma )\\ Z-1_{m_2}{\beta }^z(\theta ,\gamma )\end{array}\right)}^T{\rho }^{-1}\left(\theta ,\gamma \right)\left(\begin{array}{l}Y-1_{m_1}{\beta }^y(\theta ,\gamma )\\ Z-1_{m_2}{\beta }^z(\theta ,\gamma )\end{array}\right)\text{ (15)}$$

where $m=m_1+m_2$. In Eqs. (11) and (15), the regression constants ${\beta }^y$ and ${\beta }^z$ are

$$\left(\begin{array}{l}{\beta }^y\left(\theta ,\gamma \right)\\ {\beta }^z\left(\theta ,\gamma \right)\end{array}\right)={\left(F^T{C}^{-1}\left(\theta ,\gamma \right)F\right)}^{-1}{F}^T{C}^{-1}\left(\theta ,\gamma \right)\left(\begin{array}{l}Y\\ Z\end{array}\right)\text{ (16)}$$

Note that ${\sigma }^2$ cancels in both equations (11) and (16). The cokriging model’s hyper-parameters $\theta =\left({\theta }_1,\mathrm{...},\mathrm{...},{\theta }_d\right)$, and $\gamma$ are determined by optimizing the associated maximum likelihood estimation

$${\max }_{\theta ,\gamma }MLE\left(\theta ,\gamma \right)=-\left(m\ln \left({\sigma }^2\left(\theta ,\gamma \right)\right)+\ln \det \left(\rho \left(\theta ,\gamma \right)\right)+m\right)\text{ (17)}$$

Maximizing the likelihood is a multidimensional nonlinear optimization problem. In [37,38], this problem is solved by using the modified direct search algorithm of Hooke and Jeeves [39]. The cokriging predictor of Eq. (11) interpolates the high-fidelity sample points $x^i,i=1,\mathrm{...}\mathrm{..}\mathrm{...},m_1$ and follows the trend indicated by the low-fidelity sample points ${x^{\prime }}^j,i=1,\mathrm{...}\mathrm{..}\mathrm{...},m_2$.

So, in a variable-fidelity modeling context, the low- and high-fidelity sets of POD coefficients are used in the aforementioned cokriging model, from which a new set of high-fidelity POD coefficients is calculated. These are then used together with the POD modes to predict the high-fidelity solutions. In contrast with the previous methodology, this technique inherently models the cross-correlation between the low- and high-fidelity POD coefficients rather than augments the parametric space using an additional variable.

Mapping the low-fidelity data onto the high-fidelity grid: Grids for viscous and inviscid computations are typically very different in terms of surface discretization, point count, and element type. Since the suggested POD-based variable fidelity model requires that the evaluation of the flow-field variables be made at common grid points for both low- and high-fidelity data so that the POD modes are evaluated at a common location, an interpolation method to map the inviscid grid results onto the viscous grid is necessary. Therefore, the mesh2mesh tool for interpolating between surface grids, which was developed at DLR as an add-on to the TAU flow solver, was used. In this code, various interpolation techniques are available, including linear interpolation, radial basis function interpolation, nearest point search, and inverse distance weighting of n nearest points, where n is specified by the user. It should be noted that if only surface pressure distributions are of interest, a shared surface grid between the inviscid and viscous flow computations would suffice.

Results

This article contains results obtained in the subsonic flow regime using either the parameter augmentation methodology or cokriging to model the low- and high-fidelity POD coefficients. In this work, viscous flow solutions were considered as the high-fidelity data, while inviscid flow solutions were considered as the low-fidelity data. Moreover, all variable-fidelity modeling was performed to predict the surface pressure distribution over an RAE 2822 airfoil and a DLR F12 aircraft configuration.

Two-dimensional case

In the first example, two-dimensional aerodynamic flow around an RAE 2822 airfoil is considered. Inviscid and viscous flow solutions were computed using the TAU flow solver. Figure 1 shows (a) The grid with 13124 points on which the inviscid flow computations were made and (b) The grid with 27874 points on which the viscous flow computations were made. The flow conditions considered in this case were a freestream Mach number of 0.3, a Reynolds number of 100,000, and a temperature of 288.15 K. All inviscid snapshots were computed at the corresponding viscous reference pressure and density. In this particular case, one or two aerodynamic parameters, namely, the angle of attack and Mach number, were considered. From Figure 1, it is evident that the surface discretization of the two grids is different. Since the POD-based VFM models as described in section 2.2.3 require that the field evaluations for both low- and high-fidelity be made at common grid points, an interpolation method to map the inviscid results onto the viscous grid is necessary. It is suggested that the inviscid solution be mapped onto the viscous flow grid rather than the other way round, since more sampling points are taken into account and hence better results are obtained. Consequently, the mesh2mesh interpolation tool described previously was utilized. In particular, the linear interpolation technique was used, which resulted in a very fast mapping process. Once the mapping was achieved, the interpolated inviscid flow solutions, or snapshots, were used together with the computed viscous flow solutions.

Modeling the surface pressure distribution over an RAE 2822 airfoil using the POD-based variable-fidelity model: Eleven inviscid and viscous flow snapshots were generated with the angle of attack varying from 0° to 10° in steps of 1°. Out of this set of CFD solutions, six inviscid flow snapshots from 0° to 5° in steps of 1° were considered together with one viscous flow snapshot at an angle of attack of 1°. Thus, seven snapshots were considered in total to build the model, and the variable-fidelity modeling was conducted by cokriging. By considering a variable-fidelity model, the problem was transformed from an extrapolatory problem into an interpolation-based extrapolation setting, and hence reasonably accurate estimates were possible. Figure 2 shows the resulting Cp distribution plots at different angles of attack. It must be emphasized that the Navier-Stokes (N-S) solution was included for verification purposes only. From this figure, it can be noticed that the predictions are accurate at angles of attack of 2°, 3°, and 4°. At an angle of attack of 0°, the surface pressure distribution at x/c less than 0.2 is reasonably accurate. However, it is not as accurate as the surface pressure distribution over the remaining part of the airfoil. Also, at an angle of attack of 5°, the peak suction pressure is slightly overestimated. It is expected that as one gets further away from the viscous flow solution at an angle of attack of 1°, the predictions will get less and less accurate. To this end, when considering five inviscid flow snapshots from 6° to 10° in steps of 1° and one viscous flow solution at 10°, the predicted results from the variable-fidelity model at an angle of attack of 8° and 9°, are not as good as in the previous example at lower angles of attack. See Figure 3. While on the pressure side, the predictions are reasonably good, on the suction side the prediction is unsatisfactory, because the two pressure distributions on the suction side do not correlate. In particular, the suction pressure immediately after stagnation of the inviscid flow solution is considerably different from the viscous flow solution. Also, the pressure distributions on the upper surface are quite different. Specifically, at the trailing edge some viscous flow separation is apparent.

When a variable-fidelity model was set up taking into consideration the set of snapshots as specified in Figure 4a and two parameters, namely, the angle of attack and Mach number, the results obtained are shown in Figure 4b, Figure 4c, Figure 4d and Figure 4e.

In this work, only variable-fidelity modeling of subsonic flow is presented because, in the transonic flow regime, there are quite large differences between inviscid and viscous flow solutions. Consequently, it is suggested that this methodology is not applied to transonic flow problems with inviscid and viscous flow solutions as the low- and high-fidelity data, respectively. However, it would be interesting to apply this methodology in the transonic flow regime when considering coarse and fine domain discretizations of viscous flow solutions, for example. Obviously, due to the highly nonlinear features in the transonic flow regime, such as shock waves moving over a surface, a large number of snapshots are necessary for effective POD modeling.

Three-dimensional case

In this case, three-dimensional aerodynamic flow around a DLR F12 aircraft is considered. Inviscid flow (low-fidelity) and viscous flow (high-fidelity) snapshots were computed using the TAU flow solver. Figure 5 shows (a) The grid on which the inviscid flow computations were conducted and (b) The grid on which the viscous flow computations were made. Note that both grids were adaptively refined, but this is not strictly necessary. The inviscid flow grid has 1.9 million grid points, while the viscous flow grid has 9.6 million grid points. The flow conditions considered in this case were a freestream Mach number of 0.21, a Reynolds number of 1.2 million, and a temperature of 288.15 K. As in the previous case, the inviscid snapshots were computed at the corresponding viscous reference pressure and density. In this case, one single aerodynamic parameter, the angle of attack, was considered. Since the inviscid flow and viscous flow grids are significantly different at the surface, and since the POD-based VFM models as described in section 2.2.3 require that the evaluations for both low- and high-fidelity solutions be made at common grid points, a linear interpolation method to map the inviscid results onto the viscous grid was adopted. Once the mapping was achieved, the interpolated inviscid flow solutions, were used together with a selection of the computed viscous flow solutions in a variable-fidelity context to model the surface pressure distribution over the aircraft.

Modeling the surface pressure distribution over a DLR F12 aircraft using the POD-based variable-fidelity model: Nine inviscid and viscous snapshots were generated with the angle of attack varying from 1° to 9° in steps of 1°. Out of this set of CFD solutions, six inviscid snapshots were considered at an angle of attack of 2°, 4°, 5°, 6°, 7° and 9°, together with two viscous snapshots at an angle of attack of 3° and 8°. Thus, 8 snapshots were considered in total, and the variable-fidelity modeling was conducted by either augmenting the parametric space with an auxiliary parameter or else by cokriging. In this case, it was observed that there were very small differences between using either radial basis functions or kriging interpolation methods to construct the response surface, and consequently only the results from one interpolation method are shown. Moreover, the Cp distributions at different section locations along the wing were compared. Figure 6 shows the section locations where the Cp distribution plots were taken into consideration.

Figure 7 and Figure 8 show a comparison of the surface Cp distribution at an angle of attack of 4° and 5°, respectively, between that obtained from the POD-based VFM, the Euler solution, and the Navier-Stokes solution. It is important to note that the Navier-Stokes solution is only included as a reference here and was not used to derive the POD-based VFM result. From a glance at the various Cp distributions, one would observe that the difference between the high- and low-fidelities is small. However, a closer look at the pictures reveals that there is significant resemblance between the VFM predicted solution and the Navier-Stokes one, most especially over the vertical and horizontal stabilizers.

On having a further detailed inspection of the predicted solution by looking at the Cp distribution at different section locations along the wing as shown in Figure 6, it can be observed that, in general, the POD-based VFM methodology corrects the difference between the Euler and Navier-Stokes solutions. See Figure 9 and Figure 10, which show the Cp distribution at different section locations and at an angle of attack of 5° and 6°, respectively. The difference between the low- and high-fidelity is at its maximum at the outboard part of the wing, i.e., eta is equal to 0.8.

A comparison between the POD-based VFM predictions using different response surface construction methods is also shown in Figure 11. In this figure, the surface pressure distribution over the complete aircraft at an angle of attack of 7° is considered. From this figure, it is evident that the two different response surface techniques produce similar results, though one can still come across regions where one is superior to the other one. For example, the cokriging technology produced better results over the fuselage and the leading half of the wing, while the parameter augmentation technique produced better results over the trailing edge of the wing and the tail plane. Some of these features are conspicuous in Figure 11. However, a more detailed and quantitative comparison is presented in Figure 12, where plots of the coefficient of pressure at different wing sections are shown. This latter figure shows that the accuracy of the approach depends to a certain extent on the interpolation method of the projection coefficients, but it was also observed that on further increasing the range between the high-fidelity snapshots, the quality of the results deteriorates.

Conclusions

In this article, a POD-based variable-fidelity model was applied to predict the surface pressure distribution over an RAE 2822 airfoil and a DLR F12 aircraft configuration. An aerodynamic flow at subsonic speed was considered. In this model, a POD of a set of snapshots made up of two different levels of fidelity was performed, from which the resulting low- and high-fidelity projection coefficients along each significant POD mode were mathematically modeled using a variable-fidelity model. Two different variable-fidelity models were considered, and a comparison was made between the two. In this work, the low-fidelity data was derived from inviscid flow computations while the high-fidelity data was derived from viscous flow computations. High-fidelity predictions were made for the projection coefficients at parameter values that were not available in the training dataset for each POD mode. The POD coefficients were modeled using a multi-dimensional global response surface, which is interrogated to generate solutions at the parameter values of interest. The model was applied to correct the low-fidelity surface pressure distribution over an airfoil and a complete aircraft, and good agreement was found between the model predictions and the viscous flow CFD solutions at the same parametric values, as long as there was correlation. Also, it was found that the model is sensitive to the ensemble of observations, and the accuracy of the approach depends to a certain extent on the continuous representation of the projection coefficients. Apart from offering a means to fuse computational data of variable-fidelity, this model also alleviates the inevitable up-front cost necessary for surrogate modeling based on POD.

Acknowledgements

The work presented in this article used and extended the Surrogate Modeling for Aero-data Toolbox (SMART), a code package developed by the Aero-Loads Prediction Group, Institute of Aerodynamics and Flow Technology of the German Aerospace Center (DLR), with contributions from the Institute ‘Computational Mathematics’ of Technische Universität Braunschweig. The work was co-funded by Airbus, DLR and the German Federal State of Niedersachsen as part of the C2A2S2E project, which in turn was sponsored by the European Regional Development Fund and Economic Development Fund of the Federal German State of Niedersachsen Contract/Grant no. W3-80026826. Also, the author would like to thank his ex-colleagues at DLR, Braunschweig, for assistance.

---