Wall-Modeled Large-Eddy Simulation of Autoignition-Dominated Supersonic Combustion

Wall-Modeled Large-Eddy Simulation of Autoignition-Dominated

Supersonic Combustion

Graham V. Candler

∗

University of Minnesota, Minneapolis, Minnesota 55455

and

Niccolo Cymbalist

†

and Paul E. Dimotakis

‡

California Institute of Technology, Pasadena, California 91125

DOI: 10.2514/1.J055550

Simulations of combustion in high-speed and supersonic flows need to account for autoignition phenomena,

compressibility, and the effectsof intense turbulence.In the presentwork,the evolution-variable manifoldframework

of Cymbalist and Dimotakis (

“

On Autoignition-Dominated Supersonic Combustion,

”

AIAA Paper 2015-2315,

June 2015) is implemented in a computational fluid dynamics method, and Reynolds-averaged Navier

–

Stokes and

wall-modeled large-eddy simulations are performed for a hydrogen

–

air combustion test case. As implemented here,

the evolution-variable manifold approach solves a scalar conservation equation for a reaction-evolution variable that

represents both the induction and subsequent oxidation phases of combustion. The detailed thermochemical state of

the reacting fluid is tabulated as a low-dimensional manifold as a functionof density, energy, mixture fraction, and the

evolution variable. A numerical flux function consistent with local thermodynamic processes is developed, and the

approach for coupling the computational fluid dynamics to the evolution-variable manifold table is discussed.

Wall-modeled large-eddy simulations incorporating the evolution-variable manifold framework are found to be in

good agreement with full chemical kinetics model simulations and the jet in supersonic crossflow hydrogen

–

air

experiments of Gambaand Mungal(

“

Ignition,FlameStructureandNear-WallBurningin Transverse HydrogenJets

in Supersonic Crossflow,

”

Journal of Fluid Mechanics

, Vol. 780, Oct. 2015, pp. 226

–

273). In particular, the evolution-

variable manifold approach captures both thin reaction fronts and distributed reaction-zone combustion that

dominate high-speed turbulent combustion flows.

Nomenclature

= speed of sound,

∕

;

= progress variable, nonfuel mass fraction, and fuel

mass fraction

v;s

= mixture and

-species specific heats,

∕

⋅

= diffusion coefficient,

∕

= total energy per unit volume,

∕

e; e

= mixture and

-species specific energy,

∕

= convective flux vector

h; h

= enthalpy and total enthalpy,

∕

= grid index

= jet momentum ratio

= diffusive enthalpy flux

= kinetic energy,

∕

L; R

= values obtained from left and right data

-species molar mass,

∕

kmol

= number of species in detailed kinetics model

= element face unit normal vector and components

= pressure,

= heat flux vector

R; R

−

= eigenvector matrices

= universal (molar) gas constant,

∕

kmol

⋅

= mixture (specific) gas constant,

∕

⋅

= species index

= temperature,

= time,

= vectors of conserved and primitive variables

;

u; v; w

= velocity vector and components,

∕

= face-normal velocity component,

∕

= position vector and its components,

-species mass fraction

= dissipative flux factor

= reaction-zone thickness,

= dissipation rate

= evolution-variable source term,

∕

= Kolmogorov length scale,

;

= diagonal matrix of eigenvalues and eigenvalue

;

= turbulence field variable and kinematic viscosity

;

= density and

-species density,

∕

= viscous and Reynolds stress

= evolution variable

= stoichiometric fuel

–

air ratio

= subgrid-scale relaxation-rate parameter

-species chemical source term,

∕

= flux directions

Subscript

= equilibrium state

I. Introduction

ELIABLE simulation of turbulent combustion in high-speed

airbreathing engines is important for their design and

optimization. The most promising approach for representing these

complex flows is large-eddy simulation (LES), in which large-scale

turbulent motion is resolved and small-scale effects are modeled.

There are several widely used approaches for the LES of combustion,

but most were originally developed for combustion in low-speed

flows and in combustion regimes that can be approximated as a

collection of thin flamelets [1,2]. In high-speed flows, several effects

Presented as Paper 2015-3340 at the 45th AIAA Fluid Dynamics Conference,

Dallas, TX, 22

–

25 June 2015; received 28 July 2016; revision received 15

December 2016; accepted for publication 3 January 2017; published online 13

American Institute of Aeronautics and Astronautics, Inc., with permission. All

requests for copying and permission to reprint should be submitted to CCC

at www.copyright.com; employ the ISSN 0001-1452 (print) or 1533-385X

(online) to initiate your request. See also AIAA Rights and Permissions

www.aiaa.org/randp.

*Russell J. Penrose and McKnight Presidential Professor, Aerospace

Engineering and Mechanics. Fellow AIAA.

†

Associate, Thermal Sciences, Exponent. 5401 McConnell Ave., Los

Angeles CA 90066.

‡

John K. Northrop Professor of Aeronautics and Professor of Applied

Physics. Fellow AIAA.

2410

AIAA J

OURNAL

Vol. 55, No. 7, July 2017

Downloaded by CALIFORNIA INST OF TECHNOLOGY on July 20, 2017 | http://arc.aiaa.org | DOI: 10.2514/1.J055550

become important: compressibility that is coupled to heat release,

ignition delay, and combustion in a distributed reaction-zone (DRZ)

regime.

To illustrate these effects, consider a convected fluid element

undergoing combustion. First, the fuel gradually decomposes into

radicals; then, ignition occurs, followed by rapid oxidation and heat

release. With high convection speeds, the initial decomposition step

and ignition delay may span a large physical distance in the

combustor, and dominate overall engine performance. For example,a

1 ms ignition delay in a flow traveling at

1000 m

∕

results in a 1 m

induction length. The notional convected fluid element also

undergoes mixing with the surrounding fluid, resulting in a partially

premixed state. The intense turbulence generated by the high-speed

flow yields Kolmogorov length scales of a few micrometers. Thus, the

Kolmogorov-scale eddies are approximately an order of magnitude

smaller than the typical reaction-zone thickness of a flame, and they are

expected to penetrate reaction zones and break down local flamelet

structures into locally distributed reactions (e.g., Law [3]). This is

supported by direct numerical simulation results by Aspden et al. [4]

and subsequently experimentally confirmed by Zhou et al. [5]. This

combustion regime occurs in the near-wall burning of a fuel jet in a

supersonic crossflow: for example,asin [6] andaswillbedemonstrated

with a posteriori analysis of the current simulations. Thus, the goal of

the current work is to develop a combustion modeling approach for the

LES of high-speed combustion flows that can capture both ignition-

delay and DRZ combustion effects.

There have been many reviews of turbulent combustion modeling

approaches for the LES of high-speed flows. Fureby [7], Foster and

Miller [8], Oefelein [9], and Gonzalez-Juez et al. [10] provide recent

comprehensive surveys of the most widely used approaches. These

include the family of flamelet and flamelet/progress variable (FPV)

approaches, the linear eddy model (LEM), and transported

probability density function (TPDF), and filtered mass density

function (FMDF) methods. Each approach has its own set of physical

assumptions and limitations, resulting invarying levels offidelity and

simulation costs.

The FPV approach [11] and many other related methods [12]

assume that the combustion process occurs on a low-dimensional

manifold that can be derived from one-dimensional flames. The

manifold is precomputed and provides the thermochemical state

during a simulation, assuming a probability density function for

subgrid-scale (SGS) interactions between turbulence and com-

bustion.

The LEM [13] approach is based on a one-dimensional stochastic

representation of turbulent mixing within each grid element. It is

often used in conjunction with approaches to reduce the cost of

evaluating the chemical source term, such as in situ adaptive

tabulation (ISAT) [14].

TPDF [15] and FMDF [16] methods are related, in that they solve

transport equations for the probability density functions. Typically,

these methods represent the statistics of the thermochemical state

with an ensemble of Monte Carlo particles for which the transport is

governed by a stochastic difference equation. An important

advantage of the FMDF approach is that the chemical source term has

a closed form.

These methods have been used for the LES of supersonic

combustion flows. Recent examples include compressible FPV

simulations by Chan and Ihme [17], Saghafian et al. [18], Wang et al.

[19], and Cao et al. [20]; LEM simulations by Génin and Menon [21];

and TPDF and FMDF simulations by Delarue and Pope [22], Cocks

et al. [23], and Irannejad et al. [24].

Recently, Cymbalist and Dimotakis [25,26] proposed an approach

to model high-speed combustion flows, including the effects of

ignition delay and combustion in the DRZ regime. Their evolution-

variable manifold (EVM) approach explicitly represents ignition

delay and subsequent heat release as separate phases of combustion.

Combustion is represented by a well-stirred reactor that is given an

initial condition and is integrated in time until it reaches equilibrium.

The time domain is then converted to the evolution variable

that

tracks the variation of the gas state from the initial condition (

)

as it evolves to equilibrium (

). In the flow simulation,

computed as a field variable, and its value represents the Lagrangian

evolution of a fluid element as it flows through the computational

domain. The EVM approach precomputes and tabulates the gas state

as a low-dimensional manifold in terms of energy, density, mixture

fraction, and evolution variable. A well-stirred reactor model is

expected to represent the small-scale structure of DRZ combustion

and is used to generate the manifold. With this approach, the

complete thermodynamic state of thegas is accessible to the flowfield

computation, making it possible to develop a numerical method that

conserves energy and is consistent with the evolving thermochemical

state. Because the manifold that describes the combustion process is

precomputed, the EVM approach significantly reduces the

simulation cost relative to full kinetics models, and even to ISAT-

based methods. Thus, the EVM approach has the potential to greatly

reduce the cost and improve the accuracy of large-eddy simulations

of high-speed combustion flows.

Reliable simulation of high-speed turbulent combustion flows

requires four modeling elements: compressible flow, near-wall

turbulent motion, free or detached turbulent motion, and combustion

including ignition-delay effects. In the present work, we use standard

approaches for the first three elements, and we focus on the accurate

and efficient modeling of turbulent combustion at conditions relevant

to high-speed flight. To this end, we modify an existing

computational fluid dynamics code [27] for the wall-modeled LES

of compressible reacting flows to represent the detailed

thermochemical state of the reacting fuel

–

air mixture. The focus is

on the development of a numerical method that is rigorously

consistent with the thermochemical model so that the numerical

method correctly reflects the variations in the thermodynamic state.

We apply the proposed numerical flux function to the supersonic

reacting hydrogen injection experiments of Gamba and Mungal [6]

using the hydrogen

–

air chemical kinetics model of Burke et al. [28].

II. Evolution-Variable Manifold Approach

In this section, the evolution-variable manifold approach is briefly

discussed. Additional details may be found in the work of Cymbalist

and Dimotakis [25,26]. In the EVM framework, the state of a reacting

gas is represented in terms of four flowfield variables. These are the

thermodynamic state

, the fuel mixture fraction [2],

, and a

reaction-evolution variable

. Although the thermodynamic

quantities and mixture fraction are conserved scalar quantities, the

evolution variable has a source term

for which the value can be

determined based on ignition-delay behavior before ignition and

detailed chemical kinetics after ignition, as discussed in the

following. The detailed state of the fluid is precomputed and stored in

a four-dimensional table (manifold), spanning the thermodynamic

state space of the reacting flow in question. At the LES runtime, the

table is used by the computational fluid dynamics (CFD) code to

provide the detailed state of the reacting flow.

For high-speed combustion flows, we assume that Lagrangian

fluid elements in the DRZ regime can be modeled as well-stirred

reactors (WSRs) convected by the flow. These convected fluid

elements can be characterized by their thermodynamic state,

stoichiometry (or mixture fraction), and extent of reaction evolution.

We determine the detailed state of the fluid as a function of these

variables using WSRs with constant mass, energy, and stoichiometry:

(1a)

where

is the reactor time, and

(1b)

Species within the WSR react with mass-fraction evolution

equations given by

CANDLER, CYMBALIST, AND DIMOTAKIS

2411

Downloaded by CALIFORNIA INST OF TECHNOLOGY on July 20, 2017 | http://arc.aiaa.org | DOI: 10.2514/1.J055550

;:::;Y

;:::

;

−

;

;:::;N

(1c)

where the

are the

-species mass fractions,

is the species

’

chemical-production rate, and

is the total number of reacting

species in the WSR. The role of

is discussed in the following. These

WSRs are initialized across the range of thermodynamic state and

stoichiometry anticipated in the flow and run to equilibrium. The

results are converted from the time domain to the evolution-variable

domain as follows.

We can define a product mass fraction and its equilibrium value in

the WSR as

prod

;

prod

(2)

where the

“

”

subscript indicates that the mass fraction is evaluated

at the local thermodynamic equilibrium state. The evolution variable

is defined as the ratio of the product mass fraction to its value at

equilibrium, and its rate of change is denoted by

; these have the

form

;

prod

(3)

Thus,

tracks how far product formation has progressed relative to

its final equilibrium condition. A similar definition of the progress

variable in the FPV context was used by Cao et al. [20].

The relaxation parameter

provides a simple model for SGS

processes that promote the evolution of the gas mixture toward

equilibrium, especially where the chemical reactions are slow. In the

simulations presented in this paper, it is set to a constant value:

∕

. Appendix A examines the effect of

on the evolution rate

of the WSR; it is shown that, under robust combustion conditions,

has little effect; whereas at low-density conditions, it has the intended

effect of increasing the rate of product formation to represent the

subgrid-scale processes that drive the system toward equilibrium.

These equations provide the rate of change of the

thermodynamic and chemical state of a fluid element undergoing

combustion. A question then follows: How can constant-energy,

-density, and -stoichiometry WSRs capture the behavior of

convected reacting fluid elements, for which the energy, density,

and stoichiometry may be changing as they progress from

unburned to burned? The answer to that question lies in the

assumption of path independence, i.e

., that the state of the reacting

fluid at some location in

;Z;

space does not depend on the

path (in state space) the fluid ha

s taken to that location. For the

hydrogen

–

air chemical-kinetic system

in the present simulations,

the solution to the convected WSR equations is found to be

approximately independent of the reaction trajectory taken

through the state space [

;

]; see Appendix B for details. This

trajectory independence can be illustrated by considering two

WSRs that are initialized with different values of

;

and

integrated in time (corresponding to

). If the two reactors happen

to pass through the same point in

;

space and their

thermochemical state at that point is found to be close to one

another, then their state is traj

ectory independent. Appendix B

shows that, for the hydrogen

–

air system studied here, the

approximation of trajectory independence is very good for lean

and near-stoichiometric WSRs, and it is valid within 10% for rich

conditions. This approximation greatly simplifies the population

of the EVM manifold because the WSRs may be integrated on

constant-

;

trajectories.

Before ignition, the evolution variable

is defined by the extent of

the WSR

’

s progress toward ignition. For the hydrogen

–

air system,

this phase of the reaction is reliably predicted by detailed chemical

kinetics and requires no special handling. However, when the

induction period is poorly predicted by detailed chemical kinetics,

the EVM methodology offers an alternative data-driven approach

incorporating a characteristic ignition-delay time, and it has been

applied to ethylene [25,26]. A similar approach was used by Li et al.

[29] to represent the induction process in the simulation of

detonations. Postignition,

is defined as the ratio of the total

combustion product mass fraction to its local equilibrium mass

fraction, and its source term

is given by the chemical source term of

combustion products [Eq. (3)]. This approach has the potential to be

extended to recent lumped-reduced reaction models for longer-chain

hydrocarbons, in which fuel cracking is treated as an initial step

before oxidation [30].

Solution trajectories in the state space of the WSR equations are

used to populate a table of the chemical state of the gas mixture and

as a function of all values of

;

and

encountered in the flow

regions of interest. Once the table has been constructed, the partial

derivatives of

with respect to the table coordinates are also

computed and tabulated for use in the implicit numerical method.

Thus, during a flow simulation at any given grid element, the local

values of

;

are used to determine the corresponding values

, and the

derivatives by interpolation on the table. Then, all

additional thermodynamic variables, such as pressure, temperature,

speed of sound, gas constant, and specific heats, are constructed from

the local density and energy and the tabulated mass fractions. Here,

we use the McBride and Gordon thermodynamic curve fits [31] for

the thermodynamics computations.

Figure 1 plots

and the water (

) mass fraction as a function of

;

at fixed

for the EVM table constructed to represent the

hydrogen

–

air experiments of Gamba and Mungal [6]. The detailed

chemical kinetics model of Burke et al. [28] was used to construct

the table.

Fig. 1 Portion of the EVM table for the Gamba and Mungal experiment [6] at a constant value of energy (

30 MJ

∕

(left) and mass fraction of

(right).

2412

CANDLER, CYMBALIST, AND DIMOTAKIS

Downloaded by CALIFORNIA INST OF TECHNOLOGY on July 20, 2017 | http://arc.aiaa.org | DOI: 10.2514/1.J055550

III. EVM

–

LES Implementation

In this section, we discuss the integration of the EVM state

tabulation with the wall-modeled LES. We focus on how the

thermodynamic state and transport properties of the gas mixture are

computed, which variables are tabulated, and how the diffusive

enthalpy transport term is evaluated.

First, consider the calculation of the mixture state. Rather than

tabulating all possible thermodynamic variables and transport

properties (e.g., Saghafian et al. [18] in the context of a compressible

FPV model), we follow the approach suggested by Oevermann [32]

and used by Wang et al. [19]. That is, we tabulate the mass fractions

obtained from the detailed chemistry WSR model and compute the

local thermodynamic state from the density and energy computed by

the CFD and the mass fractions interpolated from the table. This

requires the use of thermodynamic data for the chemical species, such

as the McBride

–

Gordon enthalpy curve fits [31]. Then, we can

compute the mixture transport properties from the thermodynamic

state and tabulated mass fractions.

One additional term requires attention. In the energy conservation

equation [Eq. (6e)], the diffusive enthalpy flux

−

∂

(4)

involves the computation of the species enthalpy and gradients of the

detailed model mass fractions (assuming Fickian diffusion). Quinlan

et al. [33] proposed a simplification to this term of the following form:

≅−

∂

(5)

when this approximation is adequate, only the term

∂

must be tabulated. Their a priori analysis showed relatively minor

differences using this approach compared to the full expression [33].

In our simulations, however, we found that this approximation is not

sufficiently accurate, resulting in spurious temperature variations and

poor numerical stability. Thus, in the work presented here, we use the

full expression [Eq. (4)]. This requires storage of the values of

element centroids and the computation of

gradients.

Also note that we solve for

, which is the product mass fraction,

but tabulate the product state with

∕

. However,

is a

function of

;

only, and it can be obtained from the table by

extracting the product mass fractions at

;

. The local

value of

is then computed from Eq. (3) and used to obtain the local

mixture state.

The EVM table is four-dimensional, and a general search in this

domain would be expensive. However, it is stored as an ordered array,

and the relationship between the values of the variables

;

and their indices are known a priori. Thus, integer math can be used to

determine the bounding indices of the location in the table, and then a

tetralinear interpolation is performed to obtain thevalues of

, and

the

derivatives with respect to the tabulation variables. The table is

stored as single-precision (32 bit) real values to reduce its size; we

have verified that there are no detectable differences in the

simulations with 32 and 64 bit precision.

To summarize, the LES provides the value of

;

at each

time step to the EVM table; the table then returns the detailed fluid

composition vector

; :::;Y

; :::

and product evolution

rate,

, and its partial derivatives from the precomputed Lagrangian

well-stirred reactor and ignition-delay model. Figure 2 illustrates the

exchange of information between the LES and EVM modules during

runtime.

IV. Governing Equations

In this section, we discuss the governing equations used for the

EVM approach. As discussed previously, the detailed thermochemi-

cal state of the reacting fluid is determined at runtime from

;

]. Note that

is not a conserved variable, and we solve for

the progress variable mass fraction

[Eq. (2)]. Thus, the

compressible CFD code must solve conservation equations for the

total density

, mixture fraction density

, reaction product density

, momentum

, total energy

∕

⋅

, and a

turbulence field variable

ρν

. These equations are as follows:

∂

(6a)

∂

(6b)

∂

(6c)

∂

−

(6d)

∂

−

(6e)

∂

ρν

∂

−

∂

−

∂

ρν

−

∂

(6f)

Here,

and

are the molecular and turbulent transport mass-

diffusionvelocities. The source term in the

equation represents the

progress variable production rate

, where

is the source

term of the evolution variable, as described previously. The stress

includes the viscous stress and the Boussinesq approximation for the

turbulent Reynolds stress. Likewise, the heat flux

includes

molecular and turbulent transport. The equation for the turbulence

field variable

has the form taken from Catris and Aupoix [34] with

turbulence modeling constants from the original Spalart and

Allmaras paper [35]. See also the NASA Langley Turbulence

Modeling Resource [36] for full details and for the form and values of

the turbulence model parameters:

;

S; c

. For the

purposes of the wall-modeled LES discussed here, the improved

delayed detached-eddy simulation (IDDES) approach of Shur et al.

[37] is used to modify the source term of the turbulence transport

equation to reduce dissipation away from walls; this results in a

Smagorinsky-like SGS model [38] in regions of free or detached

turbulent motion. Thus, the near-wall Reynolds-averaged Navier

–

Stokes (RANS) turbulence model provides surface shear and heat

flux boundary conditions for the outer large-eddy simulation. This

type of wall-modeled LES is a well-developed approach and has been

used for many high-speed turbulent flow simulations.

In the preceding,

is the density of species comprising the fuel

and its reaction products. We can define a complementary partial

Fig. 2 Exchange of information between the LES and EVM modules

during runtime.

CANDLER, CYMBALIST, AND DIMOTAKIS

2413

Downloaded by CALIFORNIA INST OF TECHNOLOGY on July 20, 2017 | http://arc.aiaa.org | DOI: 10.2514/1.J055550

density of the gas mixture as

−

, which is the density of

the gas mixture that is not composed of fuel and its reaction products.

The conservation equation for

∂

(7)

with

−

for mass conservation. In the following, we use the

variable because it simplifies the derivation of the flux Jacobian

and the formulation of the numerical flux function. Also, the resulting

flux Jacobian has a form similar to that obtained for a mixture of

reacting gases [39].

Thus, the vector of conserved quantities and the corresponding

convective-flux vector for a three-dimensional flow are

;

;E;

ρν

(8)

;

ρν

(9)

where

⋅

is thevelocity normal to a cell facewith unit normal

vector

V. Numerical Method

In this section, we develop a numerical flux function for the

coupled LES

–

EVM framework that uses the tabulated thermody-

namic data. First, we develop the convective flux Jacobian required

for an upwind-biased flux; then, we discuss an approach for obtaining

second-order accurate fluxes, and finally provide a brief description

of a low-dissipation higher-order centered flux function suitable for

the LES of compressible flows.

A. Upwind Numerical Flux Formulation

Let us first consider an upwind flux formulation such as modified

Steger

–

Warming flux vector splitting [40] or Roe flux-difference

splitting [41]. Such a numerical flux is appropriate for RANS

simulations of high-speed combustion flows. For LES applications,

we use the dissipative portion of the upwind flux to stabilize a low-

dissipation centered flux.

For the upwind flux, we must diagonalize the flux Jacobian:

∂

∕∂

. This is straightforward, except for the derivatives of the

pressure with respect to the conserved variables. Following the

approach of Candler et al. [39], it can be shown that

∂

−

∕

−

∕

−

(10)

Where the diagonal matrix of eigenvalues is given as

diag

;

−

;

(11)

where

−

is the speed of sound,

and

are defined as

−

and

−

(12)

and

∕

⋅

is the total enthalpy of the gas

mixture.

The pressure derivatives

∂

(13)

must be computed holding the conserved variables

fixed. For a

multispecies mixture of thermally perfect gases, the pressure may be

written as

(14)

To calculate the pressure derivatives, we can use the EVM table

and the chain rule of differentiation to obtain the full expression for

these quantities. For example,

∂

(15)

where these derivatives are taken holding

;

;E;

ρν

fixed,

and the

derivatives are evaluated with finite differences on the

EVM table. A simpler approach is to assume that

is frozen during

the differentiation, resulting in the closed-form expressions

∂

;

;E;

ρν

∂

;

;E;

ρν

−

⋅

(16a)

∂

;

ρν

(16b)

where

, and

are the mixture internal energy per unit mass, gas

constant, and specific heat:

;

(17)

We implemented both approaches and found negligible

differences, and therefore the simpler closed-form expressions are

used. (This finding makes sense because the correct value of pressure

is obtained when the frozen

approximation is made.)

The preceding equations provide a complete expression for the

convective flux Jacobian, which can be used to formulate an upwind

numerical flux function. For example, the Roe flux is [41]

∕

−

(18)

where the tilde variables are Roe averaged using left and right data.

For modified Steger

–

Warming flux vector splitting [40], the

convective fluxes may be written as follows:

∕

−

∕

−

∕

(19)

where

are the positive and negative eigenvalue matrices.

2414

CANDLER, CYMBALIST, AND DIMOTAKIS

Downloaded by CALIFORNIA INST OF TECHNOLOGY on July 20, 2017 | http://arc.aiaa.org | DOI: 10.2514/1.J055550