# Probabilistic Analysis of Aircraft Using Multi-Fidelity Aerodynamics
Databases

Jayant Mukhopadhaya

PROBABILISTIC ANALYSIS OF AIRCRAFT USING MULTI-FIDELITY AERODYNAMICS DATABASES

Jayant Mukhopadhaya

August 2024

###### Contents

  1. 1 Introduction
    1. 1.1 Motivation
    2. 1.2 Uncertainty Quantification
      1. 1.2.1 Turbulence Modeling Uncertainty
      2. 1.2.2 Numerical Discretization Error
    3. 1.3 Multi-Fidelity Analysis
    4. 1.4 Aerodynamic Databases
    5. 1.5 Certification by Analysis
      1. 1.5.1 Uncertainty Propagation
      2. 1.5.2 Aircraft Maneuver Simulation
    6. 1.6 Aircraft Configurations
      1. 1.6.1 NASA Common Research Model
      2. 1.6.2 Generic T-tail Transport
    7. 1.7 Contributions
  2. 2 Quantifying Uncertainties in RANS CFD Simulations
    1. 2.1 Turbulent Flows
    2. 2.2 Uncertainty and Error in RANS Simulations
    3. 2.3 Eigenspace Perturbation Methodology
    4. 2.4 Implementation in SU2
    5. 2.5 Validation
      1. 2.5.1 Turbulent flow over a backward facing step
      2. 2.5.2 Turbulent flow through an asymmetric diffuser
      3. 2.5.3 Jet efflux of the NASA Acoustic Response Nozzle
      4. 2.5.4 NACA 0012 airfoil at a range of angles of attack
      5. 2.5.5 Turbulent flow over a three-element High-lift airfoil
    6. 2.6 Numerical Discretization Error vs. Turbulence Modeling Uncertainty
      1. 2.6.1 NACA0012 Airfoil
      2. 2.6.2 ONERA M6 Wing
    7. 2.7 Application to NASA CRM
  3. 3 Multi-Fidelity Gaussian Process Regression
    1. 3.1 Gaussian Process Regression
    2. 3.2 Multi-Fidelity Gaussian Process Regression
    3. 3.3 Application to NASA CRM
  4. 4 Probabilistic Aerodynamics and Controls Databases
    1. 4.1 Aerodynamics and Controls Databases
    2. 4.2 Multi-fidelity Data Generation
      1. 4.2.1 Lowest Fidelity: Vortex Lattice Simulations in AVL
      2. 4.2.2 Medium Fidelity: RANS CFD Simulations in SU2
      3. 4.2.3 Highest Fidelity: Wind/Water Tunnel Experiments
    3. 4.3 Databases for the Generic T-tail Transport Aircraft
      1. 4.3.1 Single-fidelity Databases
      2. 4.3.2 Multi-fidelity Databases
  5. 5 Certification by Analysis
    1. 5.1 Maneuver of Interest
    2. 5.2 Simulation Procedure
      1. 5.2.1 Maneuver Simulation
      2. 5.2.2 Engine Failure Simulations
      3. 5.2.3 Evaluating Success or Failure
      4. 5.2.4 Monte Carlo Analysis
      5. 5.2.5 Modifications to the Simulation
    3. 5.3 Results
      1. 5.3.1 Simulations using Single- vs. Multi-Fidelity Databases
      2. 5.3.2 Simulations using Reduced Data
    4. 5.4 Certification-based Design Decisions
  6. 6 Conclusions
    1. 6.1 Summary
      1. 6.1.1 Uncertainty Quantification
      2. 6.1.2 Multi-Fidelity Modeling
      3. 6.1.3 Probabilistic Aerodynamics and Controls Databases
      4. 6.1.4 Certification by Analysis
    2. 6.2 Contributions
    3. 6.3 Future Work
      1. 6.3.1 Uncertainty Quantification
      2. 6.3.2 Multi-Fidelity Modeling
      3. 6.3.3 Certification by Analysis

##  Chapter 1 Introduction

Engineering design is a predictive discipline involving the estimation of an
object’s future real-world behavior. The object, at the time, may only exist
as a thought in the brain, as a sketch on paper, or more recently, as a file
on a computer. The design process starts with a definition of requirements
that the object in question must fulfill. It ends with the manufacturing of
prototypes that, hopefully, confirm those requirements’ fulfillment. Depending
on the object’s complexity, the interim can be as short as a few hours or as
long as multiple years involving thousands of hours of engineering work..

Those hours are spent predicting the real-world behavior of an object that
doesn’t physically exist yet. Numerous analysis techniques progressing from
basic back-of-the-envelope calculations, through computational numerical
simulations, to prototyping and experimental testing of subsystems are
employed in this endeavor. Due to the intricacies of real-world physics,
almost none of these techniques are perfect. Each method has some uncertainty
associated with its predictions that must be taken into account by the
engineers employing them. Quantifying these uncertainties can be highly
specific to the method in question.

The word fidelity is used to refer to how closely a method can mimic real-
world behavior. High-fidelity methods are better at predicting real-world
behavior, whereas low-fidelity methods employ simplifications that introduce
uncertainties into the analyses. As an example, consider estimating the weight
of an object. A low-fidelity method would be to pick up the object and
estimate the weight based on how heavy it feels. Factors such as personal
bias, left vs. right-hand usage, and muscle soreness, would contribute to the
estimate’s uncertainty. A high-fidelity method would be to use a weighing
scale that is accurate up to one milligram to measure the weight. Direct
measurement introduces significantly fewer sources of uncertainty, and the
weight would be accurate up to $0.5$ milligrams.

Fidelity comes at a cost. This is to be expected. If high-fidelity analyses
were less expensive than low-fidelity ones, there would be no reason to use
low-fidelity analyses. Continuing with the weight estimation example, the low-
fidelity method’s only cost is the time taken to pick up and estimate the
object’s weight. The high-fidelity method incurs the additional cost of the
weighing scale. Cost minimization is often a priority, and a mix of low- and
high-fidelity analyses are needed. Greater emphasis is placed on lower-
fidelity methods in earlier stages of the design process, where rough
estimations are sufficient to make design decisions. This emphasis transfers
to higher-fidelity methods as designs progress, and more certainty in
performance metrics is required before the significant investment of creating
a prototype is made.

This dissertation focuses on combining results from analyses of varying
fidelity to create a superior prediction of an engineered object’s real-world
behavior. These predictions incorporate the analyses’ uncertainties to create
a probabilistic representation rather than a single, deterministic result. A
method to quantify the uncertainties due to modeling simplifications,
particularly computational fluid dynamics (CFD) simulations, is implemented
and validated. Finally, statistical analysis that allows for the explicit
calculation of the likelihood of meeting/failing a particular design
requirement is presented. The engineered object of choice to showcase the
work’s real-world impact is an aircraft. It represents one of the most
complex, engineered objects, and it requires years of development to design.
Standard industrial analysis techniques are employed. The design requirements
are real-world flight certification tests created by the Federal Aviation
Administration (FAA) that are in use today.

###  1.1 Motivation

Aircraft design is a complex, non-linear, multi-disciplinary problem that
requires many years and thousands of engineering hours to solve. An aircraft
comprises numerous subsystems that work together to make flight possible. The
sheer size, complexity, and number of parts make aircraft manufacturing a very
long and expensive process. It is imperative that all aspects of the
aircraft’s real-world performance are thoroughly investigated before the
manufacturing step is taken in the design process. While historical experience
in designing prior aircraft is useful, more precise analyses are required to
push the boundaries of performance.

The rapid improvement of computational capabilities in the recent past has
increased computer simulations’ use to predict various aspects of the
aircraft’s real-world behavior in a virtual setting. Coupled with advances in
understanding the underlying physics of these phenomena, this has led to the
development of computer simulations of varying sophistication and
computational cost that can describe the relevant quantities of interest (QoI)
at different levels of fidelity. These simulations allow for the critical
assessment of engineering designs significantly earlier in the design process
than previously possible with purely experimental design campaigns. These
analyses have been used for aerodynamic shape optimizations [53, 4, 18],
structural optimizations [9, 59, 130], and more recently, combined aero-
structural optimizations [48, 15].

While these simulations can use simplifying assumptions that introduce
uncertainties in their analyses, they have significantly improved the ability
to predict the satisfaction of performance-based design requirements, such as
range, passenger capacity, and weight, early in the design process. However,
performance-based metrics are not the only design requirements on an aircraft.
Governing bodies, such as the Federal Aviation Administration (FAA), set
stringent flight certification requirements that test for an aircraft’s air-
worthiness [111]. These are fundamentally different from performance-based
requirements as the outcome is binary; either the aircraft passes or fails the
test. Consequently, the ramifications of not meeting the certification
requirements are worse than not meeting performance-based requirements. Flight
certification suggests that these tests can only be performed with a full-
scale prototype. Nevertheless, learning from the trend of increased reliance
on computational analyses, virtual representations of the aircraft design can
be put through simulated air-worthiness testing to estimate the likelihood of
passing or failing a requirement.

The current methodology involves building a virtual representation of the
aircraft using aerodynamic databases that contain the force and moment
coefficients experienced by the aircraft across its expected flight envelope.
These databases are created using a single information source and at specific
milestones during the design process. They do not contain any information
about the uncertainties introduced due to the particular information source
used. Simulating a flight certification maneuver using these databases yields
a single deterministic outcome. Such a result belies the uncertainty present
in the database and wrongly assumes a $100\%$ certainty in the force and
moment coefficient data that populates the database. With rigorous handling
and propagation of these uncertainties, the deterministic result can be
converted to a more realistic likelihood of success or failure.

This work aims to perform statistical analyses on simulations of flight
certification maneuvers for commercial aircraft. This is achieved by tackling
the problem on multiple fronts. A method to quantify the uncertainties arising
from modeling assumptions in the widely used Reynolds-Averaged Navier-Stokes
(RANS) computational fluid dynamics (CFD) simulations is presented,
implemented, and validated. These simulation results are combined with data
from other information sources into aerodynamic databases that utilize multi-
fidelity models to create a stochastic representation of the aircraft’s
potential performance. These models are sampled to create hundreds of
individual representations of the databases that have small variations due to
the uncertainties in the underlying data. These samples are used to propagate
the uncertainties through the flight simulation of a current air-worthiness
test maneuver. Statistical analysis of these maneuver simulation results
yields the likelihood of an aircraft passing or failing the given test
maneuver.

This enables the consideration of flight certification metrics earlier in the
aircraft design process. It also allows for assessing the adequacy of the
planned control systems. Quantifying the risks involved with failing a
certification maneuver equips the engineers with the necessary information
required to mitigate them.

###  1.2 Uncertainty Quantification

Computer simulations and real-world experiments alike carry some inaccuracies
in their predictions. The field of uncertainty quantification (UQ) does
exactly as the name suggests; it provides methodologies to efficiently
characterize, propagate, and, in some cases, minimize uncertainties that
plague analyses. UQ has been adopted for a wide range of applications,
including modeling climate change [55, 104], understanding uncertainties in
numerical simulations [89, 39, 115], and even predicting the economic effects
of COVID-19 [6]. UQ benefits are realized when the process is taken one step
further, and the quantified uncertainties are used to make better decisions
[61]. Information on the underlying uncertainties of analysis techniques can
be used to make better medical decisions [8], create reliable and robust
designs [50, 103, 37], and create safer autonomous driving algorithms [36, 13,
127]. This work aims to propagate uncertainties in design analyses through
flight simulations to calculate the likelihood that an aircraft will succeed
or fail in meeting flight certification requirements.

Uncertainties are often divided into two categories: aleatoric and epistemic.
Aleatoric uncertainties arise due to natural variation in the parameters that
describe a situation. Epistemic uncertainties arise due to incomplete
knowledge of the situation. To exemplify the differences between the two
categories, consider a ball being thrown multiple times using the same amount
of force. Slight changes in the wind direction and speed will result in slight
differences in the distance that the ball travels. The uncertainty in the
distance traveled due to this natural variation in the wind is aleatoric. It
can be quantified by throwing the ball thousands of times, recording each
distance, and calculating statistical information, such as the mean and
standard deviation, about the result. These uncertainties are easier to
quantify through Monte Carlo simulations [40, 77] and polynomial chaos
expansions [94, 90]. However, they are irreducible. Now consider that to
remove the effects of natural variation, the same ball throw is simulated
using simple calculations employing Newton’s laws of motion. The calculated
distance is the same every time, but simplifications in the equations, such as
ignoring wind resistance, introduce uncertainty in the result. These
simplifications reflect the lack of knowledge that defines epistemic
uncertainties. Contrary to aleatoric ones, epistemic uncertainties are
reducible (through better modeling/measurement) but are difficult to calculate
and propagate as they are extremely problem-dependent [38].

Analyses of varying fidelity are used in this work. There are uncertainties
associated with each. For some techniques, the uncertainties are provided by
subject matter experts (SME) who rely on historical data and their experience
in using the analyses. Explicit quantification of epistemic uncertainties is
performed for one commonly used analysis technique: Reynolds-averaged Navier-
Stokes (RANS) computational fluid dynamics (CFD) simulations. Specifically,
the uncertainties introduced due to turbulence models and numerical error due
to insufficient discretization are quantified.

####  1.2.1 Turbulence Modeling Uncertainty

The Navier-Stokes equations are a set of non-linear partial differential
equations that describe the motion and behavior of fluids. CFD is the field of
study that involves using computers and numerical analysis techniques to solve
the non-linear Navier-Stokes equations to simulate the flow of a fluid over a
domain of interest. This is made more difficult due to turbulence, which is a
state of fluid flow characterized by chaotic, small-scale fluctuations in the
density and velocity of the fluid. The range of length and time scales that
need to be resolved through spatial and temporal discretization makes it
computationally intractable to solve exactly, i.e., without simplifying
models. Most flows of engineering interest are plagued by turbulence. The
difficulty in solving these flows has paved the way for developing a hierarchy
of solution techniques that trade computational cost for prediction accuracy.

The most widely used industry method is the Reynolds-averaged Navier-Stokes
(RANS) simulation. Steady RANS simulations are very computationally efficient
and can be used for expensive undertakings such as iterative aerodynamic shape
optimization [72, 57, 18] and aircraft database generation [125]. The
computational efficiency comes at the cost of modeling inaccuracies. The
simulations assume that the flow is steady (no time-dependent variation in the
flow) and require simplifying turbulence models that aggregate the effects of
the turbulent eddies that would be present in the flow.

The inadequacy of turbulence models in predicting certain flow features has
been well documented [118]. These models have various parameters and constants
calibrated using canonical flow conditions, such as the flow over a flat
plate. Cheung et al. in [19] treat these parameters as random variables and
use high fidelity data from direct numerical simulations (DNS) to learn
posterior distributions for these parameters. Similarly, Dow et al. in [27]
use DNS data to solve an inverse RANS problem to determine the eddy viscosity
field that would result in a flow field closest to the DNS data. Both methods
lean on computationally expensive DNS results, which are limited to simple
geometries such as flat plates [51], and channels [67, 73].

This work focuses on the eigenspace perturbation methodology [33, 52], a
physics-based UQ method that does not require any high-fidelity information
and can provide uncertainty estimates by running 6 RANS simulations. This has
been applied to wind-engineering problems [45] and to perform design
optimizations under uncertainty [84]. Chapter 2 details this methodology,
provides various validation cases, and explores the relationship between
uncertainty introduced by turbulence models and the numerical errors
introduced due to insufficient discretization.

####  1.2.2 Numerical Discretization Error

Turbulence models are not the only source of uncertainty in RANS CFD
simulations. Solving continuous equations on a discrete domain creates
numerical discretization error. In RANS CFD simulations, the continuous RANS
equations that define fluid flow are solved on a discrete domain known as the
mesh or grid. Numerical discretization error is classified as an epistemic
uncertainty as it can be reduced by increasing the number of discrete points
in the domain. As the discretization increases and the grid spacing tends
towards $0$, the numerical error approaches $0$. This is the basis for the
"Grid Convergence Study" method to quantify numerical discretization error
[3]. This is an important tool used in the verification and validation (V&V)
of CFD codes and turbulence models [113].

Section 2.6 delves into the details of quantifying numerical discretization
error. It goes on to explore its relationship with the turbulence modeling
uncertainty estimate by applying both uncertainty quantification techniques to
two benchmark simulations: a 2D NACA0012 airfoil and a 3D ONERA M6 wing.

###  1.3 Multi-Fidelity Analysis

During the course of the typical aerospace design process, different kinds of
performance analysis tools are used at different stages. Low-fidelity computer
simulations accept lower accuracy for faster computations. They are useful at
the very early stages of the design process when the aircraft’s geometry is
not well defined and is subject to significant change. For example, AVL [28]
solves the incompressible potential flow equations, taking mere seconds. It is
used to rapidly explore a large multi-dimensional design space, changing
variables like the number of engines or wing placement [12, 11]. These are
often replaced by higher-fidelity techniques, such as RANS CFD simulations, as
the design progresses, and more details of the design are finalized.
Experimental data, normally obtained through a costly wind-tunnel test,
typically provides the most accurate representation of the phenomena analyzed.
It is obtained late in the design process as expensive prototypes of
subsystems need to be manufactured.

Instead of discarding the low-fidelity simulation data when higher-fidelity
data is available, there are methods to combine data from multiple fidelity
levels to better represent the quantities of interest (QOI). Multi-fidelity
extensions to popular surrogate modeling techniques have been developed.
Polynomial chaos expansion (PCE) [94, 10], which is a popular technique for
sensitivity analysis [120, 23], can be used to combine multi-fidelity
information by learning an additive correction on the low-fidelity data [91,
97]. Gaussian processes (GP) [62, 75, 105], popular for their direct
estimation of the error in its modeling, combines multiple information sources
by learning additive and multiplicative corrections on low-fidelity data [56].
This correction of the lower-fidelity data $f_{i-1}$ to a higher-fidelity
$f_{i}$ is represented as

$f_{i}(\mathbf{x})=\rho_{i}(\mathbf{x})f_{i-1}(\mathbf{x})+\delta_{i}(\mathbf{x}),$
(1.1)

where $\rho_{i}(\mathbf{x})$ is the multiplicative term, and
$\delta_{i}(\mathbf{x})$ is the additive term. These are sometimes referred to
as bridge functions [49]. One disadvantage of using GP vs. PCE is that the GP
assumes a Gaussian distribution across the uncertainty intervals. PCE can
handle different probability distributions, but the superior multi-fidelity
handling and the direct error estimation of the GP regression make multi-
fidelity GP the tool of choice.

Multi-fidelity GP regression has been used extensively in aerospace
applications [65, 77], modeling oil reservoir production and pressures [56],
hydrodynamic simulations [47], and biological tissue growth [69]. Numerous
extensions to the framework have been proposed. Ghoreishi et al. introduced
the ability to have non-hierarchical information sources by relating each
fidelity level to the highest fidelity, as opposed to each other [42]. This is
unnecessary for this work as there is a well-defined hierarchy based on the
physics that is modeled. Perdikaris et al. [99] significantly improved
prediction capability by using deep Gaussian processes [26]. The high
computational cost and the loss of a Gaussian posterior precludes its use in
this effort. Han et al. [49] proposed using gradient information to improve
predictive capability. Unfortunately, gradient information is only available
for RANS CFD simulations through adjoint simulations [53], but not for the
other information sources.

This work does employ important improvements to the multi-fidelity GP process
that were proposed by Gratiet [46]. The new equations significantly reduced
the GP’s computational cost by splitting the data into individual information
sources. This results in the inversion of smaller matrices which is more
computationally efficient. He also extends the multiplicative term to be a
function of $\mathbf{x}$ instead of a constant. Section 3.2 presents these
equations in detail and extends them for use with noisy data for design sets
that are not nested.

###  1.4 Aerodynamic Databases

Aerodynamic databases are a representation of the aircraft’s behavior in-
flight. It contains all the forces and moments that are experienced by the
aircraft as a function of the aircraft’s configuration (control surface
deflections), orientation (angle of attack, angle of sideslip), and operating
conditions (dynamic pressure, Mach number, altitude). Calculating these forces
and moments at various points in its operating envelope is paramount to
predicting real-world performance. Most aerodynamic analyses, be it
computational or experimental, are geared towards creating a database that
catalogs these values as a function of the aircraft’s orientation and
operating conditions.

The industry standard is to have a lookup-table populated by data from
aerodynamic analyses that are performed during the design process. They get
updated as the design progresses, and the results from the higher-fidelity
analysis techniques replace the lower-fidelity data. The forces and moments
are described as multi-dimensional functions depending on up to 5 input
variables: angle of attack, sideslip angle, Mach number, dynamic pressure, and
altitude. Often only a subset of the 5 input variables is used. Discrete
analyses in this multi-dimensional domain provide data points that are used to
interpolate values between analysis locations. These databases are
deterministic and have no characterization of the uncertainties present in the
analysis techniques.

Engelund et al. [35] created databases for the Hyper-X scramjet using wind
tunnel data to analyze and predict the expected behavior before flight testing
was performed. Keshmiri et al. [58] used a mixture of CFD analyses and wind
tunnel experiments to create aerodynamic databases for the generic hypersonic
vehicle. Instead of having a lookup table, the coefficients were described
using analytic polynomial functions of arbitrary order. Databases for the Mars
Pathfinder aerodynamics created using CFD data by Gnoffo et al. [43] were
validated using flight measurements and were found to be accurate within
reason.

Each of these previous works has created deterministic expressions of the
databases. There is no quantification of the uncertainties affecting the
analyses used. Previous work by Wendorff et al. [125] introduces the concept
of probabilistic aerodynamic databases that uses multi-fidelity data and its
associated uncertainties in a Gaussian Process regression framework to create
a non-deterministic representation of the database. Using a combination of
sensitivity and uncertainty analysis, an adaptive sampling technique was
developed to find the best location to perform the next analysis to minimize
the objective function’s uncertainty at minimum analysis cost. This was
extended to include physics-based uncertainty quantification for the RANS CFD
simulations [88].

The current work follows from this and extends the databases from exclusively
describing longitudinal dynamics to include lateral-directional information.
The databases are concerned with the force coefficients of lift, drag, and
side-force, $C_{L},C_{D},$ and $C_{SF}$, respectively, and the moment
coefficients of pitch, roll, and yaw, $C_{m},C_{l},$ and $C_{n}$,
respectively. These coefficients are treated as functions of the angle of
attack ($\alpha$) and angle of sideslip ($\beta$). The databases are divided
into two parts, the aerodynamics databases that describe the bare airframe
aerodynamics, and the controls databases that describe the effect of controls
surface deflections. Simple linear combinations of the two are used to
determine the aircraft’s final forces and moments. The generic coefficient
$C_{i}$ is calculated as a function of $\alpha,\beta,$ and the control surface
deflections $\delta_{j}$ as

$C_{i}\left(\alpha,\beta,\delta_{1},...,\delta_{N}\right)=C_{i0}\left(\alpha,\beta\right)+\sum_{j}^{N}C_{i_{\delta_{j}}}\left(\alpha,\beta,\delta_{j}\right),$
(1.2)

where $i0$ refers to the bare airframe coefficient and $C_{i_{\delta_{j}}}$
refers to the incremental effect the control surface $\delta_{j}$ on
coefficient $C_{i}$:

$C_{i_{\delta_{j}}}\left(\alpha,\beta,\delta_{j}\right)=C_{i}\left(\alpha,\beta,\delta_{j}\right)-C_{i0}\left(\alpha,\beta\right).$
(1.3)

Deflections of the ailerons, rudder, elevator, flaps, and spoilers are used.

Chapter 4 delves into the details of creating multi-fidelity aerodynamics and
controls databases using multi-fidelity GP regression. This creates
probabilistic databases that can be sampled to create hundreds of individual
databases with slight variations due to the uncertainties present in the
underlying data. AVL simulations, RANS CFD simulations, and wind and water
tunnel experiments provide three different information sources of increasing
fidelity. The required data is generated, and databases are created for the
generic T-tail transport (GTT) aircraft [25].

###  1.5 Certification by Analysis

To fly a new aircraft design, it needs to go through rigorous air-worthiness
testing to ensure it is safe and can perform predictably in a variety of
different situations. These tests are defined and carried out by the Federal
Aviation Administration (FAA) [111] in the US. They occur at the very end of
the design process, once a functional full-scale aircraft is built. Failing a
certification test at this stage would require a redesign that would be
incredibly expensive. To prevent this from occurring, safety factors are
employed to account for potential uncertainty and error in the design
analyses. The quantification of uncertainties in analyses can replace
arbitrary safety factors with the explicitly calculated design margins
required to account for the errors. This is the cornerstone of reliability-
based design processes [50].

Certification by analysis (CbA) purports that the passing of certification
requirements can be done using purely simulation-based analyses. It is the
logical conclusion to the current trend of increased reliance on computer
simulations for design analysis. CbA has garnered interest from the aerospace
community with the maturing of simulation procedures. To achieve this,
simulations would have to be as accurate, if not more accurate, than what is
possible with flight testing. Many required improvements to simulation
capabilities [118] will take years to develop. The effects of uncertainties on
flight performance predictions and predicted performance in flight testing
could be quantified in the interim. This provides aircraft designers with a
method to estimate the likelihood that a design will pass or fail a
certification test.

####  1.5.1 Uncertainty Propagation

A common method for uncertainty propagation is to perform a Monte Carlo
analysis. This is a brute-force method of characterizing the effect of input
uncertainties on an output quantity of interest (QoI) [54, 122]. It involves
running multiple deterministic calculations where input variables are randomly
sampled from their respective probability distributions. The results are
analyzed to determine the effect of the variation in the input variables on
the posterior probability distribution of an output quantity of interest
(QoI). In this work’s context, the deterministic calculations are the flight
simulations, the input values are the aerodynamics and controls databases, and
the QoI is the success/failure of the certification test.

A significant advantage of using GP regression to represent the probabilistic
aerodynamic databases is the ability to take deterministic samples of the
database. Each sample is a potential candidate aircraft that could explain the
data used to create the databases. They represent aircraft that have a
slightly different performance from each other. These differences respect the
uncertainty associated with the underlying data. Since multi-fidelity GP
regression is used to model the databases, the input probability distributions
are Gaussian. The flight simulation does not represent a linear transfer
function. Therefore the posterior probability distribution of the QoI is
arbitrary. As such, cumulative density functions are used to describe the
posterior probabilities.

####  1.5.2 Aircraft Maneuver Simulation

Aerodynamic databases contain all the information needed to perform high-
fidelity flight simulations. With multi-fidelity aerodynamics and controls
databases, a virtual representation of the GTT aircraft can fly through real-
world flight certification maneuvers. With the help of industry experts at The
Boeing Company, a maneuver was picked from the FAA’s Flight Test Guide for
Certification of Transport Category Airplanes [111]. Within Chapter 5.3
Directional and Lateral Control of the guide, the Lateral Control: Roll
Capability §25.147(d) maneuver was chosen. The testing procedure is
paraphrased as:

  1. 1.

The airplane starts in a trimmed state for steady straight flight at maximum
takeoff speed.

  2. 2.

Establish a steady $30^{\circ}$ banked turn.

  3. 3.

Roll the airplane to a $30^{\circ}$ bank angle in the other direction.

  4. 4.

Aircraft must have sufficient roll authority to perform the $60^{\circ}$
change in bank angle in no more than $11$ seconds.

  5. 5.

The aircraft must do this with one engine inoperative, specifically the one
that makes this maneuver more difficult.

This part of the work leans heavily on The Boeing Company’s expertise in
flight simulation and control law mixing. Due to proprietary and patent
restrictions, the exact implementation of the flight simulation code is
unavailable, but enough information is provided to outline the simulations’
overarching methodology and workflow. This maneuver is simulated using a 5
degree of freedom (DOF) flight simulator that The Boeing Company uses in their
design processes. Details of the maneuver simulation process are discussed in
Section 5.2.

###  1.6 Aircraft Configurations

The methodologies presented in this thesis were developed and demonstrated
using two aircraft configurations: the NASA Common Research Model (CRM) and
the Generic T-tail Transport (GTT).

####  1.6.1 NASA Common Research Model

The NASA CRM is a well-investigated full-configuration aircraft [109, 108]
that was developed to create a baseline geometry upon which numerous
experimental and computational studies could be performed and compared. It was
created for the 4th Drag Prediction Workshop [86] and was used for the
subsequent 5th and 6th editions of the workshop as well [70, 87, 112, 123].
The wealth of experimental and computational data lends itself well to
showcase the performance of the uncertainty quantification and multi-fidelity
data fusion techniques developed in this thesis.

The model, designed by The Boeing Company, is based on the Boeing 777 aircraft
with a modified wing. It is a conventional tube-wing configuration designed
for a cruise Mach number of 0.85. The NASA CRM was built to be modular such
that additional components could be attached to the baseline geometry.
Consequently, different configurations of the aircraft were tested in the wind
tunnel:

  1. 1.

Baseline wing + fuselage model,

  2. 2.

Wing + fuselage + pylon and nacelle,

  3. 3.

Wing + fuselage + horizontal tail mounted at either $-2^{\circ},0^{\circ},$ or
$+2^{\circ}$.

Figure 1.1(a) shows one such configuration in a wind tunnel. The configuration
with the wing + fuselage + horizontal tail mounted at $-2^{\circ}$ is used for
this work.

(a) NASA Common Research Model

(b) Generic T-tail Transport

Figure 1.1: Wind tunnel models for the two aircraft configurations that are
investigated in this work.

####  1.6.2 Generic T-tail Transport

The GTT aircraft was developed in another collaborative effort between NASA
and The Boeing Company. This effort aimed to provide experimental and
computational data to inform a flight simulator that provides stall recovery
training to pilots [25]. The aircraft configuration derives from a remotely
piloted airplane configuration designed by NASA and built by Area-I [63]. It
is a $16\%$ scaled version of the Bombardier CRJ700, which is currently in-
service. The configuration resembles a generic short to medium range single-
aisle regional jet with a T-tail and twin aft engines. Tables 1.2 and 1.2
lists the aerodynamic reference dimensions and the mass properties of the GTT
aircraft, respectively.

Mean Aerodynamic Chord | $3.374~{}\mathrm{m}$
---|---
Wing Span | $23.159~{}\mathrm{m}$
Wing Area | $70.079~{}\mathrm{m^{2}}$
Table 1.1: Aerodynamic reference dimensions for the GTT aircraft.

Weight | $25,332~{}\mathrm{kg}$
---|---
$I_{xx}$ | $238,419~{}\mathrm{kg~{}m^{2}}$
$I_{yy}$ | $1,510,624~{}\mathrm{kg~{}m^{2}}$
$I_{zz}$ | $1,717,539~{}\mathrm{kg~{}m^{2}}$
Table 1.2: Mass properties of the GTT aircraft.

The primary motivation for using the GTT configuration for this work is the
wealth of experimental wind tunnel data that is available [24]. Three
experimental campaigns in different facilities explored various aspects of the
aircraft performance characteristics:

  1. 1.

NASA Langley Research Center 12-Foot Low-Speed Tunnel (12-Foot LST): Focus on
stability derivatives and $\alpha$ sweeps,

  2. 2.

Boeing North American Aviation Research Tunnel (NAART): Focus on control
derivatives and $\alpha$ and $\beta$ sweeps,

  3. 3.

Boeing Flow Visualization Water Tunnel (FVWT): Focus on stability derivatives
at various values of $\alpha$.

Figure 1.1(b) shows the GTT model mounted in the NAART facility. This work
only uses data from the NAART and the FVWT experimental campaigns. The data
from the 12-Foot LST is a subset of the data from the other two sources. Table
1.3 lists the scales of the models used for the various wind tunnel campaigns
and the CFD campaign performed using SU2 for this work.

Data Source | Scaling factor
---|---
NASA 12-Foot LST | $0.057$
Boeing NAART | $0.020$
Boeing FVWT | $0.019$
SU2 | $0.158$
Table 1.3: Scale factors for the different sources of data.

The NAART and 12-Foot LST models tested individual control surface parts to
define the control characteristics of the aircraft. Table 1.4 lists the
control surfaces and their associated deflection ranges. In addition to these,
the entire horizontal stabilizer could be rotated to change the its incidence
angle. However, since most of the wind tunnel tests mounted the horizontal
stabilizer at an angle of $-6^{\circ}$, only that particular incidence angle
is used.

Control Surface | Deflection Range
---|---
Ailerons | $-25^{\circ}\leq\delta_{a}\leq 25^{\circ}$
Elevator | $-30^{\circ}\leq\delta_{e}\leq 30^{\circ}$
Rudder | $0^{\circ}\leq\delta_{r}\leq 35^{\circ}$
Flaps | $0^{\circ}\leq\delta_{f}\leq 60^{\circ}$
Spoilers | $0^{\circ}\leq\delta_{s}\leq 60^{\circ}$
Table 1.4: Control surface deflection ranges for the GTT aircraft.

###  1.7 Contributions

This thesis establishes a framework to perform virtual flight testing of an
aircraft early in the design process. While the focus is on aircraft design,
the principles of multi-fidelity modeling and uncertainty quantification, and
certification testing can be applied to any design problem that requires
significant analysis.

Starting with UQ, Chapter 2 reviews the eigenspace perturbation methodology to
quantify epistemic uncertainties introduced by turbulence models into RANS CFD
simulations. The methodology is implemented in an open-source CFD solver, SU2,
to enable its widespread use in the research community. It is validated on a
bevy of test cases that range from commonly used benchmark flow conditions to
those of specific aerospace interest. The relationship between turbulence
modeling uncertainty and numerical discretization error, another common source
of uncertainty in CFD simulations, is investigated. The UQ methodology is
applied to two aircraft, the NASA Common Research Model and the Generic T-tail
Transport, to create aerodynamic databases with physics-informed
uncertainties.

Since CFD is not the only analysis technique used in the aircraft design
process, Chapter 3 presents multi-fidelity Gaussian processes (GP) that are
used to combine data from different information sources. Existing equations
for multi-fidelity GP regression are extended to use noisy data when design
sets are not nested. Multi-fidelity aerodynamic databases are modeled using
these equations. AVL simulations are the lowest fidelity level, RANS CFD
simulations are the middle-fidelity, and wind tunnel experimental data are the
highest fidelity. The benefits of multi-fidelity data fusion are elucidated
using one-dimensional aerodynamic databases created for the NASA CRM aircraft.
This is extended in Chapter 4, where the first comprehensive, multi-fidelity,
multi-dimensional aerodynamics and controls databases representing a full-
configuration aircraft’s lateral and longitudinal dynamics are presented.
These are created for the Generic T-tail Transport (GTT) aircraft.

With all aspects of an aircraft’s performance defined, Chapter 5 delves into
the simulation of an FAA flight certification maneuver used to test commercial
jets’ air-worthiness. The Monte Carlo method is used to propagate the
uncertainties in the design analyses through the flight simulation. The effect
of the input uncertainty on the aircraft maneuver is analyzed, and a
probability of succeeding/failing the certification test is explicitly
quantified. The virtual flight testing is also performed using databases based
on lower-fidelity data, representing what would be available at earlier stages
in the design process. By quantifying the likelihood of success/failure of a
flight certification maneuver, potential problems in the aircraft design can
be identified, and the risk of failure can be mitigated.

These contributions are discussed in further detail in Chapter 6. Potential
avenues for future research are presented as well.

##  Chapter 2 Quantifying Uncertainties in RANS CFD Simulations

The Navier-Stokes equations are a set of non-linear partial differential
equations that describe the motion and behavior of fluids. Broadly, there are
two classes of forces acting on a fluid: inertial and viscous. Inertial forces
refer to those exerted by the fluid’s momentum or its lack thereof. Consider
the force exerted by water coming out of a tap versus water coming out of a
fire hose. The water flowing out of the fire hose has significantly larger
momentum by virtue of its flow velocity. This results in higher inertial
forces being exerted by the water flowing out of the fire hose. Viscous forces
are a result of the fluid’s resistance to deformation. Viscosity is a measure
of this resistance. It is a property specific to the fluid in question. Honey
is more viscous than water. Pouring both fluids out of a container
demonstrates the effect of these viscous forces. Honey is more resistant to
this deformation and flows out more slowly than water.

The ratio of the inertial and viscous forces on a fluid determines the fluid’s
behavior in motion. Reynolds number $\left(Re=\frac{\rho uL}{\mu}\right)$ is a
measure of this ratio. When the inertial forces acting upon a fluid are small
compared to the viscous forces, i.e., the Reynolds number is low, the fluid
flow is laminar. This flow is smooth, without significant cross-currents or
time dependence. Laminar flow usually occurs at lower fluid velocity, is
relatively easy to predict computationally, and does not require any
simplifying models.

Fluid flow turns turbulent when the ratio of inertial to viscous forces is
high, i.e., at high Reynolds numbers. Turbulence refers to the chaotic, time-
dependent, varied-scale fluctuations in flow variables that characterize
turbulent flows. Its unsteady nature, and the range of length and time scales
of the flow eddies, make it computationally intractable to predict precisely.
Simplifying assumptions are required to make the computations feasible and
productive in an engineering situation.

Note the ambiguous use of “low” and “high” to describe the Reynolds numbers
for laminar and turbulent flow. This ambiguity is intentional. The point of
transition between the two regimes is highly dependent on the fluid and the
geometry of the flow itself. It is an active area of research [5]. This thesis
focuses on turbulent flows.

###  2.1 Turbulent Flows

In the realm of engineering, there is often a need to predict fluid flow
behavior in and around engineered products. These situations can range from
simulating the flow around an airplane to simulating the flow through laptop
cooling systems. Computational Fluid Dynamics (CFD) is the field of study that
involves using computers and complex numerical analysis techniques to solve
the non-linear Navier-Stokes equations on the flow domain of interest. As
computational capabilities have increased over the past few decades, so has
the reliance on the predictive capability of CFD simulations.

CFD simulations are made more difficult due to turbulence. The range of length
and time scales that need resolving through spatial and temporal
discretization makes it computationally intractable to solve precisely, i.e.,
without simplifying models. Turbulence plagues most flows of engineering
interest. The difficulty in solving these flows has paved the way for
developing a hierarchy of solution techniques that trade computational cost
for prediction accuracy.

Figure 2.1: Hierarchy of CFD solution techniques in order of increasing
computational cost and decreasing model influence.

Figure 2.1 presents this hierarchy in the form of a pyramid. Sitting at the
top of this pyramid are Direct Numerical Simulations or DNS. These
calculations do not employ any mathematical models and resolve all spatial and
temporal time scales to give an exact reproduction of the fluid behavior.
These time-varying (unsteady) simulations are computationally expensive and
memory intensive. At the time of writing, DNS calculations are only possible
at low Reynolds numbers and for simple geometries such as flat plates [51],
and channels [67, 73]. These limitations disqualify the use of DNS, in its
current state, for practical engineering applications.

Below DNS, in Figure 2.1, lie Large Eddy Simulations (LES) and Hybrid RANS-LES
simulations. These unsteady simulations resolve some, but not all, of the
scales of turbulence. Unresolved time and length scales are modeled using
simplifying assumptions [101]. These solution methodologies allow for the
fine-tuning of the desired influence of simplifying models through the chosen
level of spatial and temporal discretization. They are in the process of being
adopted by industry for specific use cases such as performance predictions for
high-lift configuration aircraft [114], or engine simulations.

The pyramid’s base represents the most widely used method in the industry,
Reynolds-averaged Navier-Stokes (RANS) simulations. These simulations suffer
the most from modeling inaccuracies. They assume that the flow is steady (no
time-dependent variation in the flow features) and require simplifying
turbulence models that aggregate the effects of the turbulent eddies that
would be present in the flow. These simplifications significantly reduce the
computational cost but inhibit the flow features that the simulations can
capture. While its unsteady counterpart can resolve some of the time
variations of the flow, turbulence modeling is still required, affecting
prediction accuracy. Steady RANS simulations are very computationally
efficient and can be used for expensive undertakings such as iterative
aerodynamic shape optimization [72], and aircraft database generation [125].
This chapter will focus on quantifying the uncertainties injected by
turbulence models employed in RANS simulations.

###  2.2 Uncertainty and Error in RANS Simulations

In order to use steady RANS simulations to design engineering products, it is
essential to understand its strengths and shortcomings. The extensive use of
mathematical models to accelerate these simulations makes it feasible to run
multiple function evaluations at different test conditions. Concurrently, the
inadequacy of the models in predicting real-world behavior of fluids limits
its use to flow conditions where the models are valid. Over the years, these
restrictions have been found through thorough the validation of CFD
simulations and turbulence models [92] against experimental data.

Precise language is required to have a productive discussion of uncertainties
and errors in these simulations. For this purpose, definitions and guidelines
established by the American Institute of Aeronautics and Astronautics (AIAA)
in [93] are followed. In a non-academic setting, error and uncertainty are
often used interchangeably. Ambiguity is avoided in the current discussion by
defining uncertainty as: A potential deficiency in any phase or activity of
the modeling process that is due to lack of knowledge [93].

In contrast, errors are defined as: A recognizable deficiency in any phase or
activity of modeling and simulation that is not due to lack of knowledge [93].
Two key differences surface from these definitions:

  1. 1.

Potential vs. recognizable deficiency - the effects of the uncertainty in some
model may or may not cause a knowable deficiency in the resulting prediction.
On the other hand, deficiencies introduced due to modeling errors are
identifiable upon examination.

  2. 2.

Lack of knowledge - uncertainties arise due to inadequate modeling of the
real-world physics in the simulation process. Errors exist even with complete
knowledge and often have established practices and methods to reduce them.

To solidify these definitions with examples, consider the uncertainty
introduced by turbulence models in contrast to discretization errors caused by
insufficient grid resolution. Turbulence models are required to solve the
closure problem with the RANS equations. There are numerous models and
numerous variations built upon those models. Their hyper-parameters are often
calibrated to work well in certain flow conditions. This means that the
deficiencies in the employed model may or may not present themselves,
depending on the flow regime. There is a fundamental lack of knowledge that
necessitates the development and use of so many different turbulence models.
For these reasons, the deficiencies introduced by turbulence models are
considered uncertainties.

Solving non-linear PDEs, such as the RANS equations, requires discretizing the
continuous spatial and temporal domains into discrete ones. The discrete
spatial domain is known as the grid or mesh. An insufficient number of grid
points introduces discretization errors. As the number of grid points
increases, i.e., the size of the grid spacing tends to zero, the discrete
representation approaches the original continuous domain. Correspondingly, the
discretization error approaches zero. Thus, discretization error results from
insufficient grid quality and limited computational resources rather than a
lack of knowledge. It can be quantified using grid convergence studies [3].
For these reasons, deficiencies due to insufficient discretization are
considered errors. A similar argument can be made for temporal discretization.

Numerous sources of errors and uncertainties can plague a CFD simulation.
Table 2.1 lists a sampling of these sources. Errors due to iterative
convergence, programming, and usage are easily mitigated. Computer round-off
error is inevitable yet inconsequential when compared to discretization error.
This chapter focuses on turbulence modeling uncertainties and their
relationship with discretization errors. Section 2.3 introduces the eigenspace
perturbation methodology, which provides a framework to estimate the
uncertainties introduced by turbulence models. This methodology is implemented
in the open-source CFD software, SU2 [30]. Section 2.5 presents the validation
of this methodology on a wide variety of test cases. Section 2.7 goes on to
present the methodology’s performance on the NASA CRM aircraft.

Errors | Uncertainties
---|---
Discretization | Turbulence modeling
Iterative convergence | Simulation conditions
Computer round-off | Complexity of physics
Programming |
Usage |
Table 2.1: Sources of errors and uncertainties in CFD simulations. This is not
an exhaustive list.

###  2.3 Eigenspace Perturbation Methodology

The balance of computational cost and fidelity provided by Reynolds Averaged
Navier-Stokes (RANS) simulations makes it the tool of choice in most
engineering CFD applications. Reynolds averaging starts with the decomposition
of the chaotic velocity field ($u_{i}$) into its mean ($\overline{u_{i}}$) and
fluctuating ($u_{i}^{\prime}$) velocity components such that
$u_{i}=\overline{u_{i}}+u_{i}^{\prime}$. Here $i$ denotes a coordinate
direction, $i=1,2,3$. This decomposition allows for the ensemble-averaging of
the Navier-Stokes equations [101], leading to a closure problem due to the
resulting non-linear $\overline{u_{i}^{\prime}u_{j}^{\prime}}$ term. This term
is also known as the Reynolds stress tensor, $R_{ij}$ and requires modeling.

Turbulence models are concerned with predicting the behavior of the Reynolds
stress tensor throughout a flow field of interest, without any a priori high-
fidelity data, in a computationally tractable manner. To this end, these
models often make simplifying assumptions about the tensor that can inject
significant uncertainties into their predictions. For example, a popular
assumption employed in numerous turbulence models is the Boussinesq
approximation (also known as the linear eddy viscosity hypothesis). It assumes
that $R_{ij}$ can be defined by a combination of the mean rate of strain
($S_{ij}$), an eddy viscosity ($\nu_{t}$), and the turbulent kinetic energy
($k$):

$R_{ij}=\nu_{t}S_{ij}-\frac{2}{3}k\delta_{ij},$ (2.1)

where $S_{ij}=\left(\frac{\partial\overline{u_{i}}}{\partial
x_{j}}+\frac{\partial\overline{u_{j}}}{\partial x_{i}}\right)$, $x_{i}$ is a
coordinate direction, and $\delta_{ij}$ is the Kronecker delta. This linear
eddy viscosity model purports a simplified proportional relationship between
the Reynolds stress tensor and the mean rate of strain tensor. It assumes that
the turbulent fluid is an isotropic medium where Reynolds stresses are always
aligned with the mean strain rate. Although reasonable for simple flows
without adverse pressure gradients, this assumption can severely limit the
flow features predicted by the turbulence model, thus introducing uncertainty
into the QoIs predicted by RANS simulations. The shortcomings of the
Boussinesq assumption are particularly evident in areas of the operating
envelope of an aircraft where separated flows and shock-boundary layer
interactions exist.

To better understand the effects of such an assumption, the eigenspace
perturbation methodology was developed by Emory et al. [33] and Mishra et al.
[52]. It aims to quantify model-form uncertainties arising from the use of
turbulence models in RANS simulations by introducing perturbations in the
eigenvalues and eigenvectors of the Reynolds stress anisotropy tensors
$\left(b_{ij}\right)$ that the models predict. This methodology does not rely
on any higher-fidelity data. The stress tensor is decomposed into its
anisotropic and deviatoric components as

$R_{ij}=2k(b_{ij}+\frac{\delta_{ij}}{3}).$ (2.2)

Here, $k~{}(=\frac{R_{ii}}{2})$ is the turbulent kinetic energy and $b_{ij}$
is the Reynolds stress anisotropy tensor. The anisotropy tensor can be further
decomposed into its eigenvalues and eigenvectors and represented as

$b_{ij}=Q\Lambda Q^{T},$ (2.3)

where $\Lambda$ is a diagonal matrix that contains the eigenvalues
$\lambda_{i}\in\mathbb{R}$ [41], and $Q$ is a matrix where the $i$-th column
represents the eigenvector corresponding to $\lambda_{i}$. The matrices $Q$
and $\Lambda$ are ordered such that
$\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$. To understand the eigenspace
perturbations, it helps to visualize the Reynolds stress tensor as an
ellipsoid where the eigenvectors define the ellipsoid axes, the eigenvalues
define the relative lengths along each of the axes, and the turbulent kinetic
energy defines the size of the ellipsoid. Figure 2.2 (d) shows an example of
such an ellipsoid, represented in a coordinate system defined by the
eigenvectors of $S_{ij}$.

The perturbations exercise the limits of the physical realizability
constraints placed on the Reynolds stress tensor [116, 119, 82]. One way to
graphically represent the constraints placed on the anisotropic part of the
stress tensor is to use the eigenvalues ($\Lambda$) to project the stress
tensor onto an anisotropy-invariant map. The map often takes the form of a
triangle. The triangle’s vertices represent the one-, two- and three-component
limiting states of turbulence (referred to as the $1C$, $2C$, and $3C$
states). All of the physically realizable states of the stress tensor lie
within the triangle. One such representation of an anisotropy-invariant map is
the barycentric map [7]. This mapping allows the writing of the projection as
a convex combination of the limiting states of turbulence:

$\textbf{x}=\textbf{x}_{1C}(\lambda_{1}-\lambda_{2})+\textbf{x}_{2C}(2\lambda_{2}-2\lambda_{3})+\textbf{x}_{3C}(3\lambda_{3}+1),$
(2.4)

where $\textbf{x}_{1C}$, $\textbf{x}_{2C}$, and $\textbf{x}_{3C}$ are the
coordinates of the vertices that represent the one-, two-, and three-component
limiting states of turbulence. For example, the stress ellipsoid shown in
Figure 2.2 (d) would be mapped on to the barycentric map as shown in Figure
2.2 (a).

The eigenvalue perturbations are performed such that the limiting states of
turbulence are simulated. This involves perturbing the state of the tensor,
sequentially, to each of the vertices of the triangle. The eigenvalues
resulting from these perturbations are:

$\Lambda^{*}_{1C}=\begin{bmatrix}2/3&0&0\\\ 0&-1/3&0\\\
0&0&-1/3\end{bmatrix},~{}\Lambda^{*}_{2C}=\begin{bmatrix}1/6&0&0\\\ 0&1/6&0\\\
0&0&-1/3\end{bmatrix},~{}\Lambda^{*}_{3C}=\begin{bmatrix}0&0&0\\\ 0&0&0\\\
0&0&0\end{bmatrix}.$ (2.5)

where $*$ represents the perturbed state. The eigenvector perturbation
involves changing the alignment such that the perturbed state is

$Q^{*}=Qv$ (2.6)

where $v$ is either:

$v_{max}=\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{bmatrix},~{}v_{min}=\begin{bmatrix}0&0&1\\\ 0&1&0\\\
1&0&0\end{bmatrix}.$ (2.7)

Combining the eigendecomposition from Equation (2.3), the eigenvalue
perturbations from Equation (2.5), and the eigenvector perturbations from
Equations (2.6, 2.7), the perturbed anisotropy tensor can be built as

$b^{*}_{ij}=Q^{*}\Lambda^{*}Q^{*T}.$ (2.8)

In the case shown in Figure 2.2 (b), the eigenvalues are perturbed to the 1C
limiting state. This state corresponds to changing the shape of the ellipsoid
to be infinitely long in one direction, as shown in Figure 2.2 (e).
Consequently, the eigenvector alignment is changed to maximize turbulent
kinetic energy production ($v_{max}$). This is represented in Figure 2.2 (c)
and (f). The change in alignment of the eigenvector can not be represented in
the barycentric map, which is why Figures 2.2 (c) and (b) are identical.

Figure 2.2: Schematic outline of Eigenspace perturbations from an arbitrary
state of the Reynolds stress.

The combinations of 3 eigenvalue perturbations (towards the $1C,2C$, and $3C$
states) and 2 eigenvector perturbations ($v_{min}$ and $v_{max}$) result in 5
unique sets of perturbations. Figure 2.3 shows the resulting Reynolds stress
ellipsoid shapes. This number of perturbations is five and not six because the
$3C$ eigenvalue perturbation results in a rotationally symmetric stress
ellipsoid, which means the eigenvector perturbation does not change anything.
Table 2.2 lists the eigenvalue and eigenvector combinations that make up each
perturbation. These five perturbed states correspondingly require five
different RANS simulations, in addition to the baseline simulation with the
unmodified version of the turbulence model, to provide information about the
uncertainty introduced by the turbulence model. For each simulation, the same
perturbation is applied uniformly across the entire computational domain, at
each pseudo-time step, until the simulation converges. There are now five new
realizations of the flow field and the flow field predicted by the unperturbed
turbulence model. The maximum and minimum values of any QoI predicted by these
6 (5 perturbed + 1 baseline) simulations create the interval bound for that
QoI. These interval bounds are used in the multi-fidelity framework as
uncertainty estimates for CFD simulations. Note that the additional five
simulations for the perturbed turbulence model can be run simultaneously. If
sufficient computational resources are available, the entire process of RANS
UQ can run in the same wall-clock time as a typical deterministic RANS
simulation.

$\mathbf{\Lambda^{*}}$ | $\mathbf{Q^{*}}$
---|---
$\Lambda_{1C}$ | $Qv_{max}$
$\Lambda_{2C}$ | $Qv_{max}$
$\Lambda_{3C}$ | $Qv_{max}$
$\Lambda_{1C}$ | $Qv_{min}$
$\Lambda_{2C}$ | $Qv_{min}$
Table 2.2: Combinations of eigenvalue and eigenvector perturbations that are
performed.

An inherent assumption is that the maximum and minimum value of any QoI occurs
at the corners of the barycentric triangle and not somewhere in the interior.
In [31] the bounds were monotonically increasing as the perturbation magnitude
increases. This monotonicity means that the interval bound is at its largest
at the corner when perturbing along that direction. It is still possible that
there could be a different direction of perturbation that would result in
larger interval bounds. For example, the largest interval could result from a
perturbation towards one of the sides of the barycentric triangle. For this
work, the perturbations occur towards the triangle’s vertices to exercise the
limits of the componentiality of the turbulence field.

The theoretical underpinnings of the eigenspace perturbation methodology have
been discussed in detail by Mishra, and Iaccarino [80]. They prove that the
eigenspace perturbations extend the isotropic eddy viscosity assumption to an
anisotropic relation between the mean velocity gradients and the Reynolds
stresses. This anisotropy represents the most general relationship between the
mean flow gradients, and the Reynolds stresses,
$R_{ij}=\nu_{ijkl}S_{kl}+\mu_{ijkl}W_{kl}$, where $S_{kl}$ and $W_{kl}$
represent the mean rates of strain and rotation respectively. This extension
enables the perturbed models to account for the effects of complex flow
features such as flow separation, secondary flows, and highly anisotropic
flows.

The eigenspace perturbation methodology has exhibited substantial success in
various engineering applications. This methodology has been successfully
applied to the analysis of uncertainty in flows through scramjets [34],
contoured aircraft nozzles [85, 83, 81, 2], and turbomachinery designs [32].
The methodology has been used for the design optimization under uncertainty of
turbine cascades [106] and in aerospace designs [22]. In civil engineering
applications, this methodology has been applied in the design of urban
canopies [39, 107]. Owing to the structural similarity between the RANS and
LES paradigms, this methodology has been extended and successfully applied for
the uncertainty estimation of Large Eddy Simulations and scalar flux models
[44] as well.

(a) $2C$ state with $v_{max}$ eigenvector alignment.

(b) $2C$ state with $v_{min}$ eigenvector alignment.

(c) $3C$ isotropic turbulence state.

(d) $1C$ state with $v_{max}$ eigenvector alignment.

(e) $1C$ state with $v_{min}$ eigenvector alignment.

Figure 2.3: Graphical visualization of each eigenspace perturbation as
Reynolds stress ellipsoids.

It is important to note that this methodology provides no probability
distribution information for the QoIs within these interval bounds and
assuming any particular distribution would be inconsistent with the
methodology. The multi-fidelity GP requires data to be jointly normally
distributed. Consequently, to use these interval bounds in the multi-fidelity
GP framework, a Gaussian distribution of the QoIs within the bound is assumed.
The QoI is given by $y\sim\mathcal{N}(\mu,\sigma^{2})$ where $\mu$ is the
center of the bound, and $\sigma^{2}$ is calculated such that $95\%$ of the
interval bound predicted by the RANS UQ methodology lies at $2\sigma$ from the
mean, $\mu$.

###  2.4 Implementation in SU2

SU2 is an open-source software suite used to solve non-linear Partial
Differential Equations (PDE) along with PDE-constrained optimization problems
while utilizing unstructured meshes [95]. The framework is general and meant
to be extensible to any governing equations for multi-physics problems. This
work focuses on its Reynolds Averaged Navier-Stokes (RANS) solver capable of
simulating compressible, turbulent flows found in aerospace engineering and
design problems. Extensive validation studies of the SU2 platform conducted
across a diverse assortment of turbulent flows [96] ensures the solver’s
reliability.

SU2 is actively developed around the world. It has been released under an
open-source license and is freely available to the community so that
developers may contribute to the source code and further improve the accuracy
and capabilities of the suite. To accomplish ease of use, the SU2 includes
industry-standard solver technology for turbulent flows while also developing
numerical solution algorithms that result in robust convergence rates.
Finally, to aid with gradient-based optimization, SU2 includes a discrete
adjoint solver implementation for efficiently computing the gradient of any
quantity of interest with respect to any design variable.

This section outlines details of implementing the eigenspace perturbation
methodology into the CFD solver. The first part of the discussion outlines the
algorithms and the steps required to perform the eigenspace perturbation
within the CFD solver. The second part focuses on its implementation via
computer code and its integration with the SU2 solver.

Figure 2.4: Flow chart showing the implementation of EQUiPS within SU2

Fig. 2.4 shows where the eigenspace perturbation fits within the context of
the CFD solver. SU2 carries out the perturbation for each cell of the grid at
each pseudo-time step. This process repeats until the flow simulation
converges and outputs a perturbed solution.

Figure 2.5: Schematic outlining interaction between CFD solver and
perturbation methodology, at each cell and iteration.

Fig. 2.5 focuses on the steps required to perform the eigenspace perturbation.
The first step is to use the flow field to calculate the Reynolds stress
tensor, using Equation 2.1, at each cell in the computational domain. The
anisotropic part of the tensor, $b_{ij}$ in Equation 2.2, is decomposed into
its eigenvalues and eigenvectors, as shown in Equation 2.3. The decomposition
is done efficiently by leveraging the anisotropic tensor’s symmetry. First,
the symmetric tensor is reduced to a symmetric tridiagonal form using the
Householder’s transform and accumulating orthogonal similarity
transformations. This step uses the TRED2 procedure outlined by [74]. [102]
presents the details of the routine. Then, the eigenvalues and eigenvectors of
this symmetric tridiagonal matrix are computed and sorted using the TQL2
procedure outlined in [126]. These two steps provide a complete
eigendecomposition of the Reynolds stress tensor at a cell in the domain.

The eigenvalues are used to project the stress state onto the barycentric
anisotropy invariant map using Equation 2.4. According to the type of
perturbation, the stress state, denoted by the coordinates $\mathbf{x}$, is
perturbed towards one of the vertices of the barycentric map. An under-
relaxation factor ensures solution stability. For example, if the stress state
is perturbed towards the 1-component state, the perturbed stress state is
represented by the coordinates,

$\mathbf{x^{*}}=\mathbf{x}+r\left(\mathbf{x}_{1C}-\mathbf{x}\right)$ (2.9)

where $r$ is the under-relaxation factor. $r$ is a user-tunable parameter. All
the results shown in this thesis use the default value of $r=0.1$.

Using the barycentric coordinates of the perturbed stress tensor,
$\mathbf{x^{*}}$, Equation 2.4 is refactored to get the perturbed eigenvalues.
This results in 2 equations, one for each coordinate, with three unknown
eigenvalues. The fact that the anisotropy tensor has a zero trace i.e.

$\lambda_{1}+\lambda_{2}+\lambda_{3}=0$ (2.10)

provides the third equation required to solve the system of equations for the
perturbed eigenvalues.

Similarly, for the eigenvector perturbation, the eigenvectors of the Reynolds
stress tensor at a cell are permuted to $v_{min}$ or $v_{max}$, depending on
the specific perturbation. The Reynolds stress tensor at each cell is
reconstructed using these perturbed eigenvalues and eigenvectors.
Additionally, dependent quantities like the turbulence production tensor are
re-constituted using this perturbed Reynolds stress tensor. These are input
back into the CFD solver for the next iteration. This process continues till
numerical convergence as determined by the convergence criteria.

SU2’s modular architecture enables integrating this framework without
significant alterations to the main code. The module is split into two parts:
a high-level Python script that sequentially specifies the type of
perturbation and C++ code that performs the perturbations during the execution
of simulations. For smooth operation, it is best to have performed a baseline
simulation with SU2 and have achieved sufficient convergence. It ensures that
the mesh file and the input configuration file are well-posed and can provide
converged perturbed solutions if run through the Python script.

The Python script takes an input configuration file and a mesh file identical
to the ones used to run the baseline CFD simulation in SU2. The script sets
the necessary configuration options to run the EQUiPS framework. These include
the direction of the eigenvalue perturbation (aligned towards one of the
$1C,2C$ or $3C$ state on the barycentric triangle for a specific simulation),
the magnitude of the eigenvalue perturbation ($\Delta_{B}\in[0,1]$) and the
perturbed alignment of the Reynolds stress eigenvector ($v_{min}$ or
$v_{max}$, as detailed earlier in this section). It sequentially runs through
the simulations for each perturbed state, creating a new directory for each
new simulation, and outputting the results in the respective directories. In
conjunction with the converged baseline solution, these results can be post-
processed to provide the necessary model-form uncertainty information.

The implementation within the existing code base is limited to the two areas
where the perturbations are injected into the simulation: the viscous and
turbulent flux calculations. The perturbed stress tensors replace the original
ones in the flux calculations. The new fluxes advance each node to the next
pseudo-timestep. In the SST turbulence model, [78], the effects of the
perturbations manifest themselves in the turbulence production term.

As the solution converges, the Reynolds Stress also converges to its perturbed
state. Once a perturbed solution is converged, the Python script moves on to
the next eigenspace perturbation. It creates a new directory and configuration
file to specify the new perturbation options. This process continues until all
the specified perturbed simulations are complete.

In the spirit of versatility, the implementation was designed to minimize the
number of user-defined inputs. The implementation allows the use of the
framework without in-depth knowledge of its mechanics. The details of the
perturbations (componentiality, eigenvector permutation) are abstracted away
to provide a clean user interface that does not deviate from the workflow of
running regular CFD simulations. All the options used for the module can also
be specified in the configuration file for the potential expert user without
needing the Python script.

###  2.5 Validation

Before the eigenspace perturbation methodology can be applied to flow
simulations involving complex geometries, the methodology must be tested on a
range of flow configurations often used to validate CFD codes. Validation
refers to The process of determining the degree to which a model is an
accurate representation of the real world from the perspective of the intended
uses of the model [93]. In the context of the eigenspace perturbation
methodology, validation involves ensuring the predicted uncertainty estimates
are a reasonable representation of the uncertainty injected by the turbulence
model in use. Flow configurations with either experimental or high-fidelity
simulation data are required to validate the methodology. It should meet the
following expectations:

  1. 1.

If there is a discrepancy between the RANS simulation predictions and the
experimental data, the uncertainty estimates should be large.

  2. 2.

If the RANS simulation predictions agree with the experimental data, the
uncertainty estimates should be small.

In an ideal case, if the turbulence model was the only source of uncertainty
and error, all the experimental data should lie within the uncertainty
estimates from the methodology. In practice, this is rarely the case. Often
errors stemming from insufficient numerical discretization or a discrepancy in
experimental and simulation conditions affect the quality of the simulation
prediction.

This section presents the uncertainty estimation module results for a range of
test cases detailed in Table 2.3. The first two flows are benchmark cases
widely used to test turbulence models. They test the ability of the
methodology to predict the behavior of highly separated flows. The last four
cases are flows of aerospace interest. The jet flow cases test the
methodology’s performance for shear-layer flows, whereas the airfoil flows
test its performance over a range of flow regimes. All the simulations use the
SST turbulence model[78].

Table 2.3:  Details of testing and validation cases Test Case | Rationale | Notes
---|---|---
Turbulent flow over a backward facing step | Benchmark flow | 2D, incompressible simulation.
Turbulent flow through an asymmetric diffuser | Benchmark flow | 2D, compressible simulation.
Jet effluxes from the NASA Acoustic Response Nozzle | Engineering case | 3D, sub-sonic simulation.
NACA 0012 airfoil at a range of angles of attack | Engineering case | 2D, sub-sonic simulation.
Turbulent flow over a three-element High-lift airfoil | Engineering case | 2D, sub-sonic simulation.

####  2.5.1 Turbulent flow over a backward facing step

The inlet and outlet boundary conditions are chosen to yield the correct Mach
number $\left(M_{ref}=0.128\right)$ upstream from the step. The figure is
scaled by the step height $\left(H\right)$. Note that the $x$ and $y$ axes
have different scales. The different scales help visualize the geometry since
the geometry is $\approx 20\times$ longer in the $x$ direction when compared
to the channel height. The channel height after the step is $9H$.

The flow is allowed to develop before flowing over the backward-facing step.
The sudden expansion in the channel height results in flow separation and
subsequent reattachment downstream from the step. The size of the resulting
recirculation region or the reattachment length is the critical quantity that
must be predicted accurately by a turbulence model. Classical two-equation
turbulence models underpredicted the re-attachment length by 10-25% [121]. The
underprediction occurs due to inaccurate predictions for the Reynolds stresses
arising from the use of an isotropic eddy viscosity [121]. This configuration
has been investigated experimentally by [29] and their data is used in our
investigation.

The uncertainty estimates can account for the discrepancy between model
predictions and experimental data at most locations. While there is some
discrepancy between the unperturbed RANS simulation and the experimental data,
the uncertainty intervals can account for most of this discrepancy, with
almost all experimental data points lying within the uncertainty estimates.

####  2.5.2 Turbulent flow through an asymmetric diffuser

Turbulent flow in an asymmetric diffuser has some interesting features, such
as separation over a smooth wall, subsequent reattachment, and redevelopment
of the downstream boundary layer, which offer considerable challenges to eddy-
viscosity-based models. The flow is allowed to develop in the channel leading
up to the diffuser. The channel height, $H$, is used to normalize the
geometry’s length scale. As before, the $x$ and $y$ axes have different scales
to make it easier to represent the geometry. The reference pressure,
$P_{ref}$, was adjusted to ensure the channel centerline velocity just
upstream of the diffuser is $22.6m/s$. The simulation results are validated
using experimental data of [17]. The data include mean velocities at various
stations in the diffuser and skin friction data on the upper and lower walls.

####  2.5.3 Jet efflux of the NASA Acoustic Response Nozzle

Reliable predictions of turbulent jets exhausting from aircraft nozzles are
essential for various aerospace design applications. However, these exhaust
jets involve many complications, such as the shear layer between the jet
efflux and the ambient air, complicated nozzle geometries, compressibility
effects, and under or over-expanded flows. These pose significant challenges
to eddy-viscosity-based models. Focusing on the mixing between the jet and the
ambient fluid in the vicinity of the nozzle exit, RANS models predict a
significantly lower rate of jet mixing compared with high-fidelity data.
Farther downstream of the jet’s potential core, RANS models predict the far-
field mixing rate to become substantially higher than is observed in
experiments. Similarly, the fidelity of RANS predictions is highly
inconsistent, having higher fidelity for cold jets than heated, for
axisymmetric than non-axisymmetric geometries, and even varying significantly
over different quantities of interest.

In this test case, we investigate sub-sonic jets from the NASA Acoustic
Response Nozzle. These jet efflux cases have been studied experimentally by
[14], and extensive Particle Image Velocimetry (PIV) datasets are available.
The NASA Glenn Research Center generated the PIV dataset using the Small Hot
Jet Acoustic Rig. The tests were repeated on numerous instances during
2001-2007 with varying PIV configurations. In addition to a robust mean, this
corpus provides insight into the data uncertainty. The raw data set consists
of over 23,000 velocity fields available for statistical analysis at each
location. The researchers also shared a final consensus dataset, with the
averaged mean over the raw data, and the standard deviations over the
measurements at each point, termed as error bars in the data set by [14]. This
data uncertainty attempts to account for any statistical bias and random
errors in the measurements arising from PIV calibration, image analysis
procedure, rig flow instrumentation, minor errors in the instantaneous
measurements, and the size and number of ensembles used.

Outlined here are the results for different subsonic jet conditions. Table 2.4
details the inlet ratios and the resulting Mach number at the jet throat.

Table 2.4:  Reference conditions for subsonic jet flow cases considered Test Case | Classification | $Ma$ | ${T_{T}}/{T_{ref}}$ | ${P_{T}}/{P_{ref}}$
---|---|---|---|---
Case I | heated, perfectly expanded | 0.376 | 1.764 | 1.102
Case II | cooled, perfectly expanded | 0.513 | 0.950 | 1.197

The solid black dots represent the mean experimental data, and the associated
error bars represent the reported standard deviation of the PIV data. The
solid black line represents the baseline result using the unperturbed SST
turbulence model, and the gray area is the uncertainty estimated by the
eigenspace perturbation methodology.

The differences in the solid black line and black dots show that the baseline
turbulence model has great difficulty in accurately capturing the flow physics
in the jet wake. The uncertainty estimates, represented by the grey areas,
account for a significant proportion of the discrepancy between the RANS
simulations and the high-fidelity experimental data. If the experimental data
points themselves do not lie within the grey area, there is a non-trivial
intersection between the error bars and the grey area.

####  2.5.4 NACA 0012 airfoil at a range of angles of attack

In this test case, we consider the flow over a NACA 0012 airfoil at a range of
angles of attack, varying from 0∘ to 20∘. This design was chosen specifically
due to its widespread adoption in the industry. These include in conventional
aircraft, for instance, the wingtips of the Cessna 140A, 207; helicopter
designs such as the inboard and outboard blades of the Aerospatiale AS365,
Boeing 600N, Lockheed 475; in addition to numerous horizontal and vertical
axis wind turbines. The simulation conditions are shown in Table 2.5. The
experimental data [64] report the coefficient of lift ($C_{L}$) and the
surface pressure coefficient ($C_{P}$) for different angles of attack.

Table 2.5: Simulation conditions for the NACA0012 test case. Mach Number | $0.15$
---|---
Reynolds Number | $6\times 10^{6}$
Reference chord length |  $1$ m
Freestream Temperature | $300~{}\text{K}$
$\alpha$ | $0^{\circ}\leq\alpha\leq 20^{\circ}$
Figure 2.6: Variation in the coefficient of lift $C_{L}$ with the angle of
attack $\alpha$

In Fig 2.6, we outline the variation in the coefficient of lift with the angle
of attack $\alpha$. There is almost no discernible difference between the RANS
predictions and the experimental data at low angles of attack. Accordingly,
here the uncertainty bounds are negligible. There is a substantial discrepancy
between the RANS predictions and the high fidelity data at higher angles of
attack. For these values of $\alpha$, the uncertainty estimates are
significant. At all values of $\alpha$, the uncertainty estimates envelop the
experimental data.

Away from the stall, at $\alpha=10^{\circ}$, the RANS predictions are in good
agreement with the experimental data. Accordingly, the uncertainty bounds are
negligibly small. At $\alpha=15^{\circ}$, there is a significant discrepancy
between the RANS predictions and the experimental data. Here, the uncertainty
estimates are sizable, and they envelop the experimental data. However, the
uncertainty estimates account for this discrepancy, and all the experimental
data points intersect the shaded interval estimates.

####  2.5.5 Turbulent flow over a three-element High-lift airfoil

In this test case, we investigate the turbulent flow over a McDonnell Douglas
Aerospace (MDA) single-flap, three-element airfoil. The flap rigging used
corresponds to the 30P/30N designated by MDA. The results correspond to the
case with Mach number, $Ma=0.2$; Reynolds number, $Re=5\times 10^{6}$ with an
angle of attack of $\alpha=8^{\circ}$. We use the experimental data of [20].

We outline this case as a test against the false positive. In cases where
there is a significant discrepancy in the RANS predictions, the uncertainty
bounds should exhibit the same. However, in cases where the RANS predictions
are accurate, having spurious uncertainty bounds that are significant in their
extent would be misleading and would correspond to a false positive. [60] have
tested this flow case for a range of RANS models and have found the RANS
predictions to be accurate. In such a scenario, ideally, we would expect the
uncertainty bounds to be negligible at most locations along the airfoil
sections. Figure 2.7 outlines the distribution of the coefficient of pressure
($C_{p}$) on the surfaces of the different sections. The discrepancy between
the RANS predictions and the experimental data is minimal over the main
element and the airfoil flap. Accordingly, the uncertainty bounds are
negligibly small over these zones. The upper surface of the slat exhibits an
appreciable amount of discrepancy between RANS predictions and the high-
fidelity data. The uncertainty estimates over this surface are substantial,
and they envelop the experimental data.

Figure 2.7: Pressure coefficient for the 30P30N configuration

###  2.6 Numerical Discretization Error vs. Turbulence Modeling Uncertainty

Solving continuous equations on a discrete domain creates numerical
discretization error. RANS CFD simulations involve solving the continuous RANS
equations that define fluid flow on a discrete domain known as the mesh or
grid. Increasing the number of discrete points in the domain reduces the
numerical discretization error. As the discretization increases and the grid
spacing tends towards $0$, the numerical error approaches $0$ as well. This
property is the basis for the "Grid Convergence Study" method for
quantification of numerical discretization error [3]. Details of the
methodology are reproduced here for clarity.

The methodology requires using a family of sequentially coarser grids using
the same grid generation parameters. The easiest method to create the grid
family involves generating a very dense structured mesh which would result in
minimal discretization error. Then each successive coarser grid level removes
every other grid line in each direction. In 2D and 3D computations, this
process results in the number of grid points reducing by a factor of 4 and 8,
respectively, at each grid level. It results in uniformly refined grids that
isolate the effect of the discretization on the simulations.

Three grids of successive refinement are required to calculate all the
numerical error metrics. First, a representative grid size for the $i$-th grid
is defined as

$h_{i}=\left(\frac{1}{N_{i}}\right)^{A}$ (2.11)

where $N_{i}$ refers to the number of degrees of freedom in the $i$th grid and
$A=1/2$ for 2-D simulations and $A=1/3$ for 3-D simulations. SU2 uses node-
centered numerics, so $N_{i}$ refers to the number of nodes in the grid. The
grids are ordered such that $i=1$ refers to the most refined grid (most number
of nodes), and $i=3$ refers to the coarsest grid. Correspondingly, a smaller
value of $h$ indicates a finer mesh with more nodes. Next, grid ratios are
defined as

$r_{21}=\frac{h_{2}}{h_{1}},~{}r_{32}=\frac{h_{3}}{h_{2}},$ (2.12)

and solution differences are defined as

$\epsilon_{21}=\phi_{2}-\phi_{1},~{}\epsilon_{32}=\phi_{3}-\phi_{2},$ (2.13)

where $\phi_{i}$ is the solution on the $i$-th grid. Often for aerodynamics-
related problems $\phi$ represents a force or moment coefficient such as
$C_{L}$ or $C_{m}$.

These quantities are used to compute the observed order of convergence, $p$,
by solving the following equations using a fixed point iteration:

$\displaystyle p$
$\displaystyle=\frac{1}{\ln{\left(r_{21}\right)}}\left(\ln{\left|\frac{\epsilon_{32}}{\epsilon_{21}}\right|}+q(p)\right)$
(2.14) $\displaystyle q(p)$
$\displaystyle=\ln{\left(\frac{r_{21}^{p}-s}{r_{32}^{p}-s}\right)}$ (2.15)
$\displaystyle s$ $\displaystyle=1\times
sign\left(\frac{\epsilon_{32}}{\epsilon_{21}}\right).$ (2.16)

The order of convergence should be close to the order of the numerical method
used to solve the simulations. All of the RANS calculations made using SU2 use
second order numerical methods, so we expect the observed order to be close to
$2$.

Using the grid ratios and the apparent order of convergence, the infinite grid
solution is extrapolated as,

$\phi_{ext}^{21}=\frac{r_{21}^{p}\phi_{1}-\phi_{2}}{r_{21}^{p}-1}.$ (2.17)

The approximate relative fine-grid error:

$e_{a}^{21}=\left|\frac{\phi_{1}-\phi_{2}}{\phi_{1}}\right|,$ (2.18)

can be used to calculate the fine-grid convergence index:

$GCI_{fine}^{21}=\frac{1.25e_{a}^{21}}{r^{p}_{21}-1}.$ (2.19)

Here $1.25$ is an empirically recommended Factor of Safety (FoS) based on
hundreds of CFD case studies [110]. Furthermore, the GCI can be used to
express the $95\%$ confidence interval on the fine grid solution. The solution
can be expressed as

$\phi\approx\phi_{1}\pm\left(GCI_{fine}^{21}\right)\left|\phi_{1}\right|$
(2.20)

Armed with these metrics to quantify the numerical discretization error in the
CFD simulations, comparisons to the uncertainty introduced by turbulence
models can be made. The eigenspace perturbation methodology discussed in
Section 2.3 only estimates uncertainties introduced by turbulence models. The
RANS CFD simulations required to quantify the uncertainties are run on grids
that introduce some degree of discretization error. The following sections
explore the relationship and the relative magnitudes of the two quantities.

####  2.6.1 NACA0012 Airfoil

The same NACA0012 case presented in Section 2.5.4 is used here. The 5th and
6th Drag Prediction Workshops use this case as a verification study [70, 112].
The grids used for those verification studies are used here. Table 2.6 details
the mesh metrics.

Mesh Level | Nodes | Surface Nodes | Wall spacing | Approx. $y^{+}$
---|---|---|---|---
L1 | $14,687,744$ | $4,097$ | $1.0\times 10^{-7}~{}m$ | 0.025
L2 | $3,673,856$ | $2049$ | $2.0\times 10^{-}7~{}m$ | 0.05
L3 | $919,424$ | $1,025$ | $4.0\times 10^{-7}~{}m$ | 0.1
L4 | $230,336$ | $513$ | $8.0\times 10^{-7}~{}m$ | 0.2
L5 | $57,824$ | $257$ | $1.6\times 10^{-6}~{}m$ | 0.4
L6 | $14,576$ | $129$ | $3.2\times 10^{-6}~{}m$ | 0.8
L7 | $3,584$ | $65$ | $6.4\times 10^{-6}~{}m$ | 1.6
Table 2.6: Details of the grid family used to perform numerical discretization
error quantification for the NACA0012 case.

To be consistent with the results form the eigenspace perturbation
methodology, the SST turbulence model is used for these simulations. The
$x$-axis for the figures is the representative mesh size $h$ (Equation 2.11).
Going from right to left on the $x$-axis corresponds to increasing grid
refinement, from the L7 mesh to the L1 mesh. As the grid refinement increases
the coefficients converge towards the infinite grid solution.

Ideally, to compare the numerical discretization error and the turbulence
modeling uncertainty, the finest three meshes should be in the asymptotic
region of convergence for the quantities of interest. For this case the
$GCI_{fine}^{43}=5.0\times 10^{-5}$ which represents a $95\%$ confidence
interval of $1$ drag count ($10^{-}4$). This value is a suitably low numerical
discretization error that would not adversely affect the RANS CFD simulations.
The finer meshes, L2 and L1, are computationally expensive to perform
simulations. So to compare the numerical discretization error and the
uncertainties predicted by the eigenspace perturbation methodology, we use the
L3, L4, and L5 meshes. The L6 mesh, which has significant numerical error as
evidenced by its under-prediction of $C_{L}$ and over-prediction of $C_{D}$,
is also used to see how the RANS UQ methodology handles insufficient
discretization that might not capture relevant flow features.

The RANS UQ methodology is applied to each of the meshes in question, L3-L6.
The same $\alpha$ sweep that was used in Section 2.5.4,
$0^{\circ}\leq\alpha\leq 20^{\circ}$, is used for these simulations. As the
mesh refinement increases, going from L6 to L3, the simulation results tend to
agree more with the experimental data. As expected, increasing the number of
grid points in the mesh reduces the numerical discretization error in the
simulations. The color of the uncertainty area corresponds to the mesh level
represented. It is hard to distinguish between the uncertainty areas from the
L3-L5 meshes because all of these uncertainty areas are almost coincident. The
notable exception is the uncertainty estimate for the L6 mesh in blue. These
two observations indicate that while a minimum grid refinement that captures
the relevant flow features is required, further grid refinement does not
significantly impact the uncertainty estimate predicted by the eigenspace
perturbation methodology. This conclusion is significant as it means that
uncertainty quantification can be performed using coarser meshes that reduce
the computational cost of performing the six simulations required to engender
the uncertainty estimate.

The L3-L5 meshes are used to compute the numerical error metrics. The error
bars on the simulation data points represent the $95\%$ confidence interval on
the numerical discretization error as defined by Equation 2.20. The blue
shaded area is the turbulence uncertainty estimate. For most of the points,
the discretization error bars are too small to be visible at this scale. The
turbulence uncertainty is greater than the discretization error bars for most
simulations, even those at lower angles of attack. However, talking about the
simulations at higher angles of attack requires nuance. The mesh family is
unlikely to be in the asymptotic range of convergence in the stall region of
the $\alpha$ sweep, and so the numerical discretization error metrics are not
necessarily valid for those simulations. The solutions for $\alpha=12^{\circ}$
suffered from oscillatory convergence, and the ones for $\alpha=20^{\circ}$
yielded a negative apparent order. The error bars for these two cases are
omitted as they are not valid.

From this exploration of a 2D case, it seems that, given sufficient grid
refinement to capture the relevant flow physics, the uncertainty arising from
turbulence modeling remains independent of further grid refinement. The
following section explores a 3D case to confirm these observations.

####  2.6.2 ONERA M6 Wing

A classic CFD validation case, the ONERA M6 wing, is chosen as the 3D case to
investigate the relationship between numerical discretization errors and
turbulence modeling uncertainties. It is a swept, untwisted, semi-span wing
that uses a symmetric airfoil, the ONERA D, as its cross-sectional shape. The
simulations conditions are presented in Table 2.7. A family of structured
meshes is available from the NASA Turbulence Modeling Resource website. Table
2.8 details the subset of the structured meshes used. The finest mesh level
from this family, the L1 mesh, has $\geq 64$ million points. It would be very
computationally expensive to perform the eigenspace perturbation simulations,
and so it is omitted from this study.

Table 2.7: Simulation conditions for the ONERAM6 wing. Mach Number | $0.84$
---|---
Reynolds Number | $14.6\times 10^{6}$
Reference chord length |  $0.805$ m
Freestream Temperature | $300~{}\text{K}$
$\alpha$ | $3.06^{\circ},6.06^{\circ}$
Mesh Level | Nodes | Approx. $y^{+}$
---|---|---
L2 | $8,677,681$ | $1.0$
L3 | $1,087,665$ | $2.0$
L4 | $136,793$ | $4.0$
Table 2.8: Details of the grid family used to perform numerical discretization
error quantification for the ONERA M6 case.

The transonic flow over the wing creates a strong $\lambda$-shock on the
wing’s upper surface. Accurate predictions for the relative locations of the
two shocks and their meeting point requires sufficient mesh refinement. It is
also a problematic feature for RANS CFD solvers to predict correctly as the
turbulence model can introduce significant uncertainty in the results. The
black squares represent the baseline RANS CFD predictions, and the gray shaded
region represents the estimated uncertainty that the turbulence model injects.
The finest mesh data point, shown on the left of these plots, has the
numerical discretization error bars associated with it. Using just three
meshes results in error metrics for the finest mesh only. The uncertainty
estimated by the EQUiPS module is significantly larger than the errors
introduced due to insufficient numerical discretization. Additionally, the
mesh refinement does not significantly impact the size of the uncertainty
estimates.

When the angle of attack is increased to $\alpha=6.06^{\circ}$, a different
story unfolds. The error metrics for the discretization error are no longer
valid as the grids are not sufficiently large to capture all the relevant flow
physics. Consequently, increasing the grid refinement changes the uncertainty
estimates drastically.

In summary, the 2D and 3D investigations into the relationship between
numerical discretization errors and the turbulence modeling uncertainty
estimates yield similar results and provide guidelines for the grids that
should be used with the EQUiPS module. A minimum level of grid refinement that
can capture all the relevant flow physics is required. More specifically, the
grids should be in the asymptotic region of grid convergence for the
particular flow configuration in question. This requirement ensures that the
numerical discretization error is small and that the turbulence modeling
uncertainty estimate is the dominant source of error in the simulations. Once
a grid of sufficient refinement is used, additional refinement does not
significantly impact the turbulence modeling uncertainty estimate. These
conclusions are kept in mind when creating grids for use with the EQUiPS
module.

###  2.7 Application to NASA CRM

Section 2.5 demonstrated the eigenspace perturbation methodology has been
demonstrated on various test cases. Here, the methodology is applied to a
full-configuration aircraft. The results in this section were first published
in [88] and are reproduced here. The aircraft configuration is identical
(without accounting for aeroelastic deflections of the model) to the one
tested in the wind tunnel. It closely resembles a Boeing 777 aircraft with a
fully redesigned wing. Table 2.9 describes the transonic simulation
conditions. Note that this range of angles of attack at a free-stream Mach
number of 0.85 lead to non-linear physical phenomena and flow separation that
can significantly increase the uncertainties in the RANS predictions.
Separated flow exists for $\alpha>4^{\circ}$ at this Mach number. All of the
necessary RANS CFD simulations were conducted using the SU2 [30] solver. The
SST turbulence model [78, 79] was used and the previously defined
perturbations, Equations (2.5) - (2.7), were applied to it.

Table 2.9: Simulation conditions for the NASA CRM. Mach Number | $0.85$
---|---
Reynolds Number | $5\times 10^{6}$
Reference chord length |  $7.00532$ m
Freestream Temperature | $310.928~{}\text{K}$
$\alpha$ | $-2^{\circ}\leq\alpha\leq 12^{\circ}$
$N$ (mesh elements) | $11.8\times 10^{6}$

Figure 2.8 shows details of the unstructured mesh used for the CFD
simulations. The computational domain consists of $11.8\times 10^{6}$ mixed
elements ($4.6\times 10^{6}$ nodes), which corresponds to a coarse mesh based
on the grid convergence studies performed for multiple solvers and grid
topologies [124]. The difference between the grid-converged drag at a fixed
$C_{L}$ of 0.5 and the present study is $\approx 4\%$; this difference can be
used as a rough metric for the discretization error. This error could be
reduced through additional grid resolution studies, but it is not a primary
factor for the present demonstration.

(a) Surface mesh of the NASA CRM.

(b) Close up of the nose cone showing boundary layer cells on the symmetry
plane.

(c) Details of the wing surface mesh.

Figure 2.8: Images of the NASA CRM mesh that was used for the CFD simulations.

Figure 2.9 shows the convergence history for the simulations required to
characterize the RANS uncertainty at a particular operating condition
($\alpha=2.35^{\circ}$). Most cases achieve six orders of magnitude reduction
in the density residual. For the $1C,v_{max}$ case, the residuals stall, but
the solution is said to be converged once the force coefficients stabilize to
six orders of magnitude (less than $10^{-6}$ change in force coefficients over
$100$ non-linear iterations).

Figure 2.9: Convergence history of the density residual for RANS UQ CFD
simulations at $\alpha=2.35^{\circ}$, Baseline SST turbulence model + $5$
perturbed simulations.

In this figure, the solid black line represents the predictions made by the
baseline SST turbulence model, the grey area represents the uncertainty
estimates predicted by the eigenspace methodology, and the black crosses
represent the wind tunnel data. These wind tunnel data points have error bars
associated with them, but these are barely discernible on the scale of the
plot.

The performance of the UQ module is illustrated by comparing the variation of
the coefficients of lift ($C_{L}$), drag ($C_{D}$), and longitudinal pitching
moment ($C_{m}$) due to the angle of attack ($\alpha$), as predicted by the
CFD simulations to those experimentally determined. At low angles of attack,
the flow remains well attached to the aircraft body; therefore, the turbulence
model does not introduce significant uncertainty in its predictions.
Accordingly, the interval bounds predicted by the UQ module are relatively
small. At higher angles of attack when there is flow separation over portions
of the aircraft, turbulence models struggle to make accurate flow predictions
due to the unsteady nature of the flow features and the structural limitations
of the isotropic eddy viscosity enforced by the Boussinesq assumption. The
growing uncertainty estimates predicted by the module reflect this
shortcoming.

This plot illustrates the paradigm change that RANS UQ methodologies such as
the one described in this paper can bring about. A drag polar is often used in
aerospace engineering to understand an aircraft’s behavior across its
operating range at a given Mach number. Design engineers use the deterministic
values represented by the solid black line in the drag polar plot to determine
the optimal operating condition for an aircraft and the aircraft performance
characteristics. Traditionally, only the solid black line is available to the
aircraft designer, and they use conservative factors of safety to build
operating margins into the design. The grey areas representing the possible
variability of the drag polar explicitly quantify the uncertainties injected
by the turbulence model. Instead of relying on a single deterministic value to
design the aircraft around, the uncertainty in the performance prediction can
inform the most robust optimal operating condition, and the reliability of the
design choices can be quantified.

Since the model-form uncertainty introduced by the turbulence model is
negligible at low angles of attack, the CFD predictions of the force and
moment coefficients should adhere more closely to the experimental data. Error
in RANS CFD simulations can originate from multiple sources, including, but
not limited to: discretization error due to insufficient mesh resolution,
modeling error due to the turbulence models, or any errors due to geometric
discrepancies in simulated and real-world objects. If there were only
turbulence-model-related uncertainties in the simulations, we would see all
the experimental data points lie within the grey uncertainty bounds predicted
by the RANS UQ methodology. However, we observe a significant deviation of the
simulation data from the wind tunnel experiments. This deviation manifests as
a bias and a difference in the slope of the lift coefficient curve. As was
discovered at the 4th Drag Prediction Workshop (DPW) [70], the wind tunnel
model of the NASA CRM underwent aeroelastic deformation that affected the
values recorded for the force and moment coefficients. The deformations were
more extensive than expected, and as a result, the shape of the model analyzed
using numerical simulations was different from that of the wind tunnel model.
New computational geometries that accurately reflected the deformed wind
tunnel shape were used for the subsequent workshop [87]. For this work, we use
the older computational geometry (without the wind-tunnel-induced
deformations) to exercise the ability of the multi-fidelity method to learn
biases existing in the data.

On a swept wing like the one on the NASA CRM model, the added aeroelastic
twist results in a decrease in the angle of attack of the wing-tip region
relative to the numerically analyzed rigid shape. This discrepancy leads to a
lower $C_{L}$ than the one calculated using CFD. Moreover, the increase in the
overall twist of the wing unloads the wing-tip sections, effectively
displacing the center of lift of the wing upstream and leading to the observed
increased values of $C_{M}$ when compared to those calculated using CFD. These
differences introduced new uncertainties (of an aeroelastic nature) in the
numerical predictions that escape the UQ analysis.

This discussion serves as an important reminder that regardless of the level
of model (in)adequacy of a simulation method, there may be unforeseen
uncertainties and errors introduced in the predictions that can cause results
to deviate from real-world experiments. In this particular case, the
unexpected aeroelastic deformation of the CRM model in the wind tunnel
resulted in slightly different geometries being numerically and experimentally
analyzed. The biases and errors introduced due to such unknowns are not
explicitly derived here, but they are present in the data. The auto-regressive
formulation of the multi-fidelity GP will learn these biases and errors from
the data. If the high-fidelity data is truly the most accurate representation
of the modeled physical system, it is good that the multi-fidelity GP can
learn the bias between the lower-fidelity simulations and the high-fidelity
data and compensate for it. However, if there is an error in the high-fidelity
data, the multi-fidelity GP will learn based on the erroneous data and still
assume it to be the most accurate source. This observation highlights the
importance of the hierarchy of data sources in this formulation. This work
employs a clearly defined hierarchy based on the physics captured by each
information source. Wind tunnel experiments form the highest fidelity level,
followed by RANS CFD simulations, while vortex-lattice methods form the lowest
fidelity level. In general, the hierarchy could be defined based on the
uncertainty of the predictions or expert-informed trust in the information
source.

Another critical point is that the interval predictions from the RANS UQ
methodology only estimate the uncertainty in simulations due to the turbulence
model used. The methodology cannot estimate other sources of error as it does
not rely on any high-fidelity data. For example, discretization error due to
insufficient mesh quality can bias the performance predictions from CFD
simulations. The RANS UQ methodology does not predict this bias. Instead, we
rely on the multi-fidelity GP framework to learn this and any other biases
from the differences between the low- and high-fidelity data.

As mentioned at the end of Section 2.3, these bounds contain no probability
distribution information. Nonetheless, for multi-fidelity modeling, it is
assumed that the distribution of the QoIs within the bounds is Gaussian and
symmetric about the middle of the interval. Additionally, the standard
deviation ($\sigma$) of the Gaussian distribution is defined such that the
extents of the interval estimates are $2\sigma$ away from the middle of the
interval. This definition means that for any CFD data point with RANS UQ
uncertainty estimate, the middle of the predicted interval is regarded as the
mean of the Gaussian distribution of the prediction, and the extent of the
interval bounds are $\pm 2\sigma$ away from the mean.

The RANS UQ methodology does not only provide interval estimates on integrated
quantities like $C_{L}$, $C_{D}$, and $C_{m}$. Since the eigenspace
perturbations result in different realizations of the entire flow field, the
data can be post-processed to provide valuable insight into the mechanics of
the turbulence model and the regions of the flow field that contribute to the
resulting uncertainties. Figure 2.10 depicts iso-surfaces of areas where the
local Mach number varies by greater than $0.2$ across all the perturbed
simulations. This Mach variability $(M_{v})$ is defined at every point in the
computational domain as $M_{v}=max(M_{i})-min(M_{i})$ where $i$ refers to each
realization of the flow field ($5$ perturbed + $1$ baseline flow fields) and
$M_{i}$ represents the Mach number at each point in that flow field.

At low angles of attack, Figure 2.10(a), the Mach variability is low and
limited to small regions in the flow field. The small regions indicate that
the eigenspace perturbations do not cause significant changes in the flow,
resulting in smaller uncertainty bounds. As the angle of attack increases, as
shown in Figures 2.10(b) and 2.10(c), larger areas of variability appear where
the shock is expected: at the upper surface of the wing and away from the
leading edge. These regions denote an uncertainty in the shock location. This
area grows rapidly until it reaches the leading edge in Figure 2.10(d),
signaling large uncertainty estimates and reduced confidence in the CFD
predictions. Such visualizations allow us to analyze the relationship between
the dominant flow features and their role in introducing uncertainty in the
turbulence models.

(a) $\alpha=1^{\circ}$.

(b) $\alpha=2.35^{\circ}$.

(c) $\alpha=3^{\circ}$.

(d) $\alpha=4^{\circ}$.

Figure 2.10: Isosurfaces representing areas where local Mach variability
$M_{v}=0.2$ at various angles of attack.

Similarly, the variability in any other flow quantity can be analyzed. Knowing
which areas contribute to the uncertainty in performance predictions can aid
design decisions. These results can inform sensor placement when moving to
experimental campaigns. For example, a higher density of pressure sensors can
be used in areas with large variability in the pressure coefficient ($C_{P}$).
Flow visualization techniques can be focused on areas with large velocity
variability. From the perspective of turbulence modeling, these variability
visualizations can shed light on the types of flow features that are hard to
predict using RANS simulations. This additional data processing provides more
qualitative applications for the RANS UQ methodology.

##  Chapter 3 Multi-Fidelity Gaussian Process Regression

Most computational or experimental analysis techniques provide realizations of
a quantity of interest (QoI) at discrete points in the domain of interest.
Each analysis, referred to as a function evaluation, has a monetary and
computational cost. In practice, higher-fidelity function evaluations are more
costly than lower-fidelity ones. If the converse were true, there would be no
reason to use a lower-fidelity analysis in place of a higher-fidelity one.
Fidelity here refers to how closely an analysis technique mirrors real-life
physics.

In theory, the discretization of the domain can be fine enough that it results
in a nearly continuous representation of the QoI. Additionally, using only the
highest-fidelity analysis techniques would minimize the uncertainty or error
in the QoI predictions. In reality, this is not tractable. Often cost
minimization is a priority. This constraint equates to representing QoIs as
sparsely as possible and lower-fidelity analysis techniques replacing higher-
fidelity ones wherever valid.

Numerous statistical methods that use these discrete realizations to create
continuous representations of the QoI, have been developed. These
representations, called surrogate models, assume that the underlying function
of interest varies smoothly between the discrete data points. It allows us to
squeeze the most out of the available, limited function evaluations. Popular
surrogate modeling techniques include radial basis functions [98], Gaussian
processes (GP) [62, 75, 105], stochastic collocation [71], and polynomial
chaos expansions (PCE) [94, 10].

Of these, the use of multi-fidelity data has been developed for GP [56, 46]
and, more recently, for PCE [91, 97]. Both methods handle the multi-fidelity
data by using correction terms trained on the difference between the low-
fidelity and high-fidelity data. GP have a slight advantage in this regard as
they use an additive and a multiplicative term in the correction, whereas PCE
only uses an additive term. PCE have a distinct advantage in sensitivity
analyses as Sobol indices can be directly post-processed [120, 23]. On the
other hand, GP can handle uncertain inputs and directly estimate the error in
its modeling. This feature can be used for adaptive sampling techniques that
suggest additional function evaluations to reduce the uncertainty in the model
[128]. For this work, GP are the surrogate model of choice due to its advanced
handling of multi-fidelity data that have associated uncertainties and the
direct estimation of uncertainty in its predictions.

This chapter introduces the fundamental equations of GP regression used to
handle single-fidelity data that can have uncertainties associated with it. It
builds upon this by introducing multi-fidelity GP regression that combines
data from different sources to build a single, superior surrogate model for
the QoI. These multi-fidelity GP equations are then used to create
probabilistic aerodynamic databases representing the performance
characteristics of a model aircraft, the NASA CRM. The benefits of using
multi-fidelity data vs. single-fidelity data are emphasized as well.

###  3.1 Gaussian Process Regression

The basic building block of the multi-fidelity framework is Gaussian Process
(GP) regression [105]. It is a supervised learning technique used to build a
surrogate model for an unknown function $y=f(\mathbf{x})$ given $n$ observed
input-output pairs $\mathcal{D}=(\mathbf{x}_{i},y_{i})$ for
$i\in\\{1,...,n\\}$. The function can be non-deterministic and have Gaussian
noise, $\sigma$, associated with its observations. The unknown function can
have multi-dimensional inputs but is assumed to have a scalar output. These
input-output pairs can be arranged in matrices $X$, $\mathbf{y}$ and $\sigma$.
If the function has an $m$ dimensional input then $X$ is an $\left(n\times
m\right)$ matrix of inputs, and $\mathbf{y}$ and $\sigma$ are $\left(n\times
1\right)$ vectors of outputs and associated uncertainties, respectively.

In the context of this work, while computer simulations are deterministic,
they can have modeling errors and uncertainties. These are treated as Gaussian
noise in the function of interest with a $\sigma$ proportional to the errors
and uncertainties. Experimental data is not deterministic as factors such as
natural variation in environmental conditions, or sensor limitations, can
introduce uncertainties in the data. These are also estimated as Gaussian
noise.

Since these observations can be imperfect, each observation is assumed to
carry some Gaussian noise associated with it such that
$y_{i}\sim\mathcal{N}(E(f(\mathbf{x}_{i})),\sigma^{2}(\mathbf{x}_{i}))$.
Assuming that all the observations in $\mathcal{D}$ have a joint Gaussian
distribution, a GP can be used as a surrogate model for the data. A GP is
completely defined by its mean function, $\mu(\mathbf{x})$, and a kernel
function $k(\mathbf{x,x^{\prime};\theta})$ that is parameterized by some
hyperparameters $\theta$. For this study, the squared exponential function is
used as the kernel function:

$k\left(\mathbf{x,x^{\prime}}\right)=\sigma_{f}^{2}\exp\left(-\sum_{d=1}^{m}\frac{\left(x_{d}-x^{\prime}_{d}\right)^{2}}{2l_{d}}\right),$
(3.1)

where $m$ is the dimension of the input. The hyperparameters for this kernel
function are the signal variance $\sigma_{f}^{2}$, and the length scales
$l_{d}$. The kernel function is used to create a kernel matrix
$K\in\mathbb{R}^{n\times n}$ where
$K_{ij}=k\mathbf{\left(x_{i},x_{j}\right)}$.

To enable the GP to estimate functions with a non-zero mean, the mean of
$f(\mathbf{x})$ is represented using $p$ fixed basis functions,
$\mathbf{h(x)}$, and learned regression coefficients $\beta$. At a minimum,
these basis functions include a constant term but can have multiple polynomial
terms. With these in mind, the surrogate model evaluated at a location of
interest $\mathbf{x}_{*}$ can be represented as a mean value,
$\mathbf{h(\mathbf{x}_{*})}^{T}\beta$, plus a zero-mean GP.

The $n_{*}$ sample locations and the basis functions can also be arranged in
matrices $X_{*}\in\mathbb{R}^{n_{*}\times m}$ and $H\in\mathbb{R}^{p\times
n_{*}}$ such that each row of $X_{*}$ is a $m$-dimensional sample location and
each column of $H_{*}$ is a $p$-dimensional result of the basis functions at
the locations in $X_{*}$.

Combining the GP regression equations for noisy observations with those
incorporating explicit basis functions and writing in the matrix notation, the
surrogate model is defined its mean and error estimate:

$\mu(X_{*})=H_{*}^{T}\hat{\beta}+K(X_{*},X)[K(X,X)+\text{diag}(\sigma_{i})]^{-1}(\mathbf{y}-H^{T}\hat{\beta}),$
(3.2)
$\sigma^{2}(X_{*},X_{*})=K(X_{*},X_{*})-K(X_{*},X)[K(X,X)+\text{diag}(\sigma_{i})]^{-1}K(X,X_{*}),$
(3.3)

where $\hat{\beta}=(H^{T}V^{-1}H)^{-1}H^{T}V^{-1}y$ is the best linear
estimator for the regression coefficients and
$V=K(X,X)+\text{diag}(\sigma_{i})$ represents the kernel matrix at the
observed points $\left(K(X,X)\right)$ and includes the Gaussian noise that is
associated with each observation $\left(\sigma_{i}\right)$. The prediction
from the surrogate model is defined by the mean $\mu(X_{*})$ and the
uncertainty associated with these predictions is represented by the diagonal
of the $\sigma^{2}(X_{*},X_{*})$ function.

To fully define the GP, the hyperparameters of the kernel function need to be
learned from the data. The marginal likelihood of the model as a function of
the hyperparameters is expressed as:

$p(\mathbf{y}|\theta,X)=\int p(\mathbf{y}|f,X)p(f|\theta,X)df.$ (3.4)

Often the log of the marginal likelihood is used. In the case of a Gaussian
process with a Gaussian prior the log marginal likelihood function is

$\log~{}p(\mathbf{y}|\theta,X)=-\frac{1}{2}\log|V|-\frac{1}{2}\alpha^{T}V^{-1}\alpha-\frac{n}{2}\log
2\pi,$ (3.5)

where $\alpha=\left(\mathbf{y}-H^{T}\hat{\beta}\right)$ [105]. The
hyperparameters are chosen to maximize this log marginal likelihood.

For consistency across sections, the following low- and high-fidelity analytic
functions will be used to show the functioning of the GP regression process:

$\displaystyle f_{LF}(x)$
$\displaystyle=0.5\left(6x-2\right)^{2}\sin{\left(12x-4\right)}+10\left(x-0.5\right)-5.$
(3.6) $\displaystyle f_{HF}(x)$
$\displaystyle=2f_{LF}(x)-20x+20+\sin{\left(10\cos{\left(5x\right)}\right)}.$
(3.7)

The high-fidelity function differs from the low-fidelity function by a few
polynomial terms and has a high-frequency variation absent in the low fidelity
approximation.

Using only the high-fidelity function, we show a basic example of how single-
fidelity GP regression uses discrete function evaluations to create a
continuous representation of the QoI. The solid red line represents this
function. The data points, shown as black circles, are uniformly distributed
between $0$ and $1$. The mean prediction of the GP is the solid black line,
and the gray area represents the $2\sigma$ error estimate. The error estimate
from the GP regression goes to zero near these data points, where the
uncertainty in the value of the underlying function is zero. The error
estimate increases between the data points, where the uncertainty in the
modeling parameters introduces uncertainty in the surrogate model prediction.
The larger, near-constant size error estimate shows that GP regression
respects these uncertain data inputs, and the GP incorporates it into the
prediction of the underlying function.

An essential feature of GP regression used extensively in this work is the
ability to create samples of the GP that are potential representations of the
underlying function and respect the error estimates from the GP. A sample
mean, $\mu_{S}(X_{*})$ at $X_{*}$ locations is generated as

$\mu_{S}(X_{*})=\mu(X_{*})+\sigma^{2}(X_{*},X_{*})U,$ (3.8)

where $U$ is a $n_{*}\times 1$ vector of random variables drawn from a
standard normal distribution $\mathcal{N}(0,1)$. Each colored line is a
separate sample and represents a potential candidate for the underlying
function being estimated, based on the limited provided data.

###  3.2 Multi-Fidelity Gaussian Process Regression

It is often the case that simulations or experiments of sufficiently high
fidelity are too expensive to perform over the entire domain of interest for a
modeled problem. Often, there are lower-fidelity approximations available that
can be evaluated quickly to perform parameter studies. The multi-fidelity GP
aims to use data from different fidelity levels to create a surrogate model
that can best approximate the highest-fidelity function and its uncertainty
while reducing the required number of high-fidelity function evaluations.

Assume there are $s$ information sources $f_{t}(\mathbf{x})$, where
$t\in\\{1,2,...,s\\}$, and the function at the highest fidelity level,
$f_{s}(\mathbf{x})$, is being approximated using a Gaussian Processes,
$Z_{s}(\mathbf{x})\sim\mathcal{N}(\mu_{s}(\mathbf{x}),\sigma_{s}^{2}(\mathbf{x}))$.
An auto-regressive formulation of the multi-fidelity framework is used. This
was first put forward in [56] and was improved upon by [46] to reduce
computational cost and improve predictions. The GP approximation at the $t$
fidelity level is modeled as

$Z_{t}(\mathbf{x})=\rho_{t-1}(\mathbf{x})Z_{t-1}(\mathbf{x})+\delta_{t}(\mathbf{x}),$
(3.9)
$\rho_{t-1}(\mathbf{x})=\mathbf{g}_{t-1}^{T}(\mathbf{x})\beta_{\rho_{t-1}},$
(3.10)

where $\mathbf{g}_{t-1}(\mathbf{x})$ is a set of $q$ basis functions, similar
to $\mathbf{h}(\mathbf{x})$ in the previous section, $\beta_{\rho_{t-1}}$ is
the learned regression coefficients, and $\delta_{t}(\mathbf{x})$ is modeled
using a GP. A way to interpret these terms is to consider
$\delta_{t}(\mathbf{x})$ the additive bias and $\rho_{t-1}(\mathbf{x})$ the
multiplicative bias between fidelity levels $t$ and $t-1$. To account for the
different fidelity levels and their corresponding data, the subscript $t$ is
added to the notation introduced in Section 3.1. For example, $X_{t}$ refers
to all the input data at level $t$. Additionally, the term
$\Sigma_{t}=\text{diag}(\sigma^{2}_{i,t})$ is introduced, which refers to the
noise in the outputs $\mathbf{y_{t}}$.

In Appendix B of [46], Gratiet presents the predictive equations for the case
when the design sets are not nested ($\mathcal{D}_{t}\notin\mathcal{D}_{t-1}$)
and the data has no process noise, such that $\Sigma_{t}$ is a null matrix.
This work extends those equations to include process noise
$\Sigma_{t}\neq\emptyset$, which produces the following representations for
the mean and covariance equations for fidelity level $t\neq 1$ as

$\begin{split}\mu_{t}(X_{*})=&~{}\rho_{t-1}\left(X_{*}\right)\mu_{t-1}\left(X_{*}\right)+H_{*}^{T}\beta_{t}\\\
&+\left[\left(\rho_{t-1}\left(X_{*}\right)\rho_{t-1}\left(X_{t}\right)^{T}\right)\odot\sigma^{2}_{t-1}\left(X_{*},X_{t}\right)+K_{t}\left(X_{*},X_{t}\right)\right]\\\
&\times\left[\left(\rho_{t-1}\left(X_{t}\right)\rho_{t-1}\left(X_{t}\right)^{T}\right)\odot\sigma^{2}_{t-1}\left(X_{t},X_{t}\right)+V_{t}\right]^{-1}\\\
&\times\left(\mathbf{y}_{t}-\rho_{t-1}\left(X_{t}\right)\odot\mu_{t-1}\left(X_{t}\right)-F_{t}^{T}\beta_{t}\right),\end{split}$
(3.11)

and

$\begin{split}\sigma^{2}_{t}(X,\tilde{X})=&~{}\left(\rho_{t-1}\left(X\right)\rho_{t-1}(\tilde{X})^{T}\right)\odot\sigma^{2}_{t-1}(X,\tilde{X})+K_{t}(X,\tilde{X})-\\\
&\left[\left(\rho_{t-1}\left(X\right)\rho_{t-1}\left(X_{t}\right)^{T}\right)\odot\sigma^{2}_{t-1}\left(X,X_{t}\right)+K_{t}\left(X,X_{t}\right)\right]\\\
&\left[\left(\rho_{t-1}\left(X_{t}\right)\rho_{t-1}\left(X_{t}\right)^{T}\right)\odot\sigma^{2}_{t-1}\left(X_{t},X_{t}\right)+V_{t}\right]^{-1}\\\
&\left[\left(\rho_{t-1}\left(X_{t}\right)\rho_{t-1}(\tilde{X})^{T}\right)\odot\sigma^{2}_{t-1}(X_{t},\tilde{X})+K_{t}(X_{t},\tilde{X})\right],\end{split}$
(3.12)

where $(X,\tilde{X})$ are generic input arguments,
$V_{t}=K_{t}\left(X_{t},X_{t}\right)+\Sigma_{t}$, and
$\rho_{t-1}(X)=G_{t-1}(X)^{T}\beta_{\rho_{t-1}}$. $G_{t-1}(X)$ is a $q\times
n$ matrix where each column is a $q$-dimensional result of the basis functions
for the corresponding $m$-dimensional row of input samples in $X$, and
$\beta_{\rho_{t-1}}$ are learned regression coefficients. For the lowest
fidelity level, $t=1$, the regular GP regression equations, Equations 3.2 and
3.3, are used. For a set of sample locations $X_{*}$, the mean predictions at
fidelity level $t$ is given by $\mu_{t}(X_{*})$ and the variance in the
predictions is given by the diagonal of $\sigma^{2}_{t}(X_{*},X_{*})$.

To fully define the GP of each fidelity level, the regression coefficients
($\beta_{\rho_{t-1}}$ and $\beta_{t}$) and the hyperparameters of the kernel
functions of each fidelity level need to be learned from the data. The
parameter estimation equations from [46] are extended for noisy observations:

$\begin{bmatrix}\beta_{t}&\beta_{\rho_{t-1}}\end{bmatrix}=\left[J_{t}^{T}\left(K_{t}(X_{t},X_{t})+\Sigma_{t}\right)^{-1}J_{t}\right]^{-1}\left[J_{t}^{T}\left(K_{t}(X_{t},X_{t})+\Sigma_{t}\right)^{-1}\mathbf{y_{t}}\right],$
(3.13)

with $J_{1}=H_{1}$ and for $t>1$,
$J_{t}=\begin{bmatrix}G_{t-1}\odot\left(\mu_{t-1}\left(X_{t}\right)\mathbf{1}_{q_{t-1}}\right)&F_{t}\end{bmatrix}$.
$\mathbf{1}_{q_{t-1}}\in\mathbb{R}^{q_{t-1}\times n_{t}}$ is a matrix of ones.
The hyperparameters of the kernel functions are learned by minimizing the log
marginal likelihood of each fidelity level:

$\log~{}p(\mathbf{y}_{t}|\theta_{t},X_{t})=-\left(\frac{1}{2}\log|V_{t}|+\frac{1}{2}\alpha_{t}^{T}V_{t}^{-1}\alpha_{t}+\frac{n_{t}}{2}\log
2\pi\right),$ (3.14)

where
$\alpha_{t}=\left(\mathbf{y_{t}}-\rho_{t-1}\beta_{\rho_{t-1}}-F_{t}\beta_{t}\right)$.

The potential advantage of using multi-fidelity data to estimate a function of
interest is shown by using the analytic functions defined by Equations 3.6 and
3.7 as the low-fidelity approximation and the high-fidelity function of
interest, respectively. For this case, the number of high-fidelity data points
is restricted to $4$. There are not enough data points for the GP regression
to learn all the nuanced trends in the high-fidelity function of interest. The
low-fidelity data does not have some of the higher-frequency information
present in the high-fidelity data, but it approximates the general trends
reasonably well. In this case, the low-fidelity data bolsters the scant $4$
high-fidelity data points, and the multi-fidelity GP that combines both sets
of data provides a more accurate representation of the underlying function of
interest. It is also essential to notice the large error estimates in the GP
prediction for both cases. More high-fidelity data would be required to reduce
the uncertainty in the GP modeling.

As mentioned earlier in this section, the recursive formulation put forth by
Gratiet [47] improves on the work initially done by Kennedy and O’Hagan [56]
by reducing the computational complexity of the training and sampling steps of
the multi-fidelity GP. The Gratiet implementation achieves this improvement by
splitting the dataset into individual fidelity levels instead of agglomerating
the data from all levels into one set of equations. This step results in
inverting smaller matrices, which dramatically improves the computational cost
of the process. The same analytic functions used thus far, Equations 3.6 and
3.7, are used to train the GP regressions. In this case, the number of high-
fidelity data points ($n_{HF}$) and low-fidelity data points ($n_{LF}$) had a
constant ratio: $\frac{n_{HF}}{n_{LF}}=0.2$. These savings in computational
time increase with more fidelity levels and higher-dimensional functions.

###  3.3 Application to NASA CRM

To provide more quantitative applications of the RANS UQ methodology
introduced in Chapter 2, we demonstrate its inclusion in the multi-fidelity GP
framework. In the multi-fidelity setting, the CFD data is augmented with low-
fidelity simulations using the Athena Vortex Lattice (AVL) code [28], and
high-fidelity experimental data from the wind tunnel campaigns [109, 108] used
in the preceding section. For the vortex-lattice simulations, the uncertainty
information is provided by subject matter experts (industry users and
academics), while for the experimental data, the uncertainty intervals used
are those described in the wind-tunnel campaign reports.