Computational framework for Launch, Ascent, and Vehicle Aerodynamics (LAVA)

Accepted Manuscript Computational framework for Launch, Ascent, and Vehicle Aerodynamics (LAVA)

Cetin C. Kiris, Jeffrey A. Housman, Michael F. Barad, Emre Sozer, Christoph Brehm, Shayan Moini-Yekta

PII: DOI: Reference:

S1270-9638(16)30178-X http://dx.doi.org/10.1016/j.ast.2016.05.008 AESCTE 3667

To appear in:

Aerospace Science and Technology

Received date: Revised date: Accepted date:

9 February 2015 7 April 2016 7 May 2016

Please cite this article in press as: C.C. Kiris et al., Computational framework for Launch, Ascent, and Vehicle Aerodynamics (LAVA), Aerosp. Sci. Technol. (2016), http://dx.doi.org/10.1016/j.ast.2016.05.008

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Computational Framework for Launch, Ascent, and Vehicle Aerodynamics (LAVA) Cetin C. Kiris$a , Jeffrey A. Housman,$$a , Michael F. Barad,$$a , Emre Sozer,$$b , Christoph Brehm,$$b , Shayan Moini-Yekta,$$b a NASA

b Science

Ames Research Center, Moffett Field, CA 9403 and Technology Corporation, Moffett Field, CA 94035

Abstract The Launch Ascent and Vehicle Aerodynamics (LAVA) framework, developed at NASA Ames Research Center, is introduced. The solver originated from addressing some of the key challenges that were present during the re-design of the launch infrastructure at Kennedy Space Center. The solver combines Computational Fluid Dynamics (CFD) capabilities with auxiliary modules, such as Conjugate Heat Transfer (CHT) and Computational Aero-Acoustics (CAA). LAVA is designed to be grid agnostic, i.e., it can handle block-structured Cartesian, generalized curvilinear overset and unstructured polyhedral grids either in stand-alone or by coupling different grid types through an overset interface. A description of the spatial discretizations utilized for each grid type, along with the available explicit and implicit time-stepping schemes, is provided. The overset grid coupling procedure for Cartesian and unstructured mesh types, as well as the CHT and CAA capabilities are discussed in some detail. Several NASA mission related applications are also presented to demonstrate the capabilities for large-scale applications, such as pressure, thermal and acoustic analyses of the geometrically complex launch environment, steady and unsteady ascent aerodynamics, plume-induced flow separation analyses of heavy lift launch vehicles and aeroacoustic applications. In addition, two validation cases related to NASA aeronautics applications are presented: the 1st AIAA Sonic Boom Prediction Workshop test cases and a computational study of slat noise. Keywords: Computational Fluid Dynamics, Unstructured Grids, Curvilinear, Immersed Boundary Methods, Overset Methodology, Hybrid Grids

$ Branch

Chief, Computational Aerosciences Branch, NAS Division, MS N258-2 Scientist, Computational Aerosciences Branch, NAS Division, MS N258-2

$$ Research

Preprint submitted to Elsevier

May 25, 2016

1. Introduction Computational Fluid Dynamics (CFD) has achieved global acceptance as a mature discipline that complements traditional wind tunnel and flight testing. In a modern, fast-paced design environment where decisions are tightly supported/driven by extensive simulation data, increasing pressure has been applied on CFD practitioners to improve geometric fidelity and reduce turnaround times. This trend represents a paradigm shift that favors efficient and versatile CFD frameworks in place of specialized legacy CFD solvers typically optimized for a limited set of problem types. The emphasis on CFD simulation turnaround time highlights several bottlenecks in the simulation process, the most significant of which is preparation of computational geometry representation and meshing. Several different meshing and numerical discretization strategies such as Cartesian[13, 14, 2, 10, 11, 107, 1, 70] and unstructured[62, 12] have emerged as alternatives to the classical structured curvilinear[100, 79, 55, 71, 106, 28, 39, 80, 29, 77] approach (see Figure 1 for examples). Yet, the latter method is still widely favored for the combined quality and numerical efficiency of the simulations once the relatively laborious mesh generation process has been completed. The inherent properties of structured meshes has led to the development of widely used CFD solvers, such as INS3D [79], CFL3D [55], OVERFLOW [71], and DPLR [106]. However, structured grid generation is a largely time-consuming and manual procedure that requires significant user expertise, particularly for complex geometries. Although semi-automation and scripting of common procedures through packages such as Chimera Grid Tools (CGT) [29] have helped ease the process, the main obstacle for using structured curvilinear grids is the ratio between grid generation time and solution time. In cases where a large number of simulations around a given geometry are required, the fast turnaround time of structured solvers may justify the initial investment in grid generation. However, for a limited number of computations, or a design environment that has rapid geometry changes, the cost of grid generation may become dominant. This motivates approaches that require little, if any, investment in grid generation. Unstructured grid technology has made significant strides in the past decade. Unstructured grid approaches are highly flexible and require low user input in terms of grid generation. With a variety of mainstream grid generation software packages available, complex geometries can be modeled in a fraction of the time required by traditional structured curvilinear techniques. However, critics of unstructured grid methods often cite three issues: (1) unstructured grids are computationally less

2

(a)

(b)

(c)

Figure 1: Computational grid strategies for a typical launch environment applications: (a) block-structured Cartesian AMR, (b) unstructured arbitrary polyhedral, and (c) generalized curvilinear meshes.

efficient to solve on, (2) retaining quality for the high-aspect ratio boundary layer cells is problematic, and (3) user control over the grid quality is reduced, especially for large-scale, complex geometries. Despite these shortcomings, successful implementations have been achieved by CFD solvers such as USM3D [36, 37], FUN3D [3], Loci/CHEM [59] and TAU [84]. Structured Cartesian methods are essentially free of manual volume mesh generation, are computationally efficient, and well-suited to implementation of higher-order methods and Adaptive Mesh Refinement (AMR)[13, 14]. However, resolving boundary layers in practical Reynolds number ranges is prohibitively expensive due to the isotropic nature of the mesh. Cartesian methods lend themselves particularly well for inviscid design optimization applications and have attained widespread use, notably Cart3D [1, 70]. It is clear that a general conclusion about which method is best does not exist, but the decision must be made on a problem-by-problem basis, or even at a sub-problem level, where a combination of these methods for different regions of the flow field may represent an optimal solution to efficiency and accuracy. With this in mind, the LAVA framework is designed to be grid-flexible, i.e., it can handle block-structured Cartesian/curvilinear or unstructured grids either in stand-alone mode or by coupling different grid types through an overset interface. A primary objective of the LAVA development is to increase the level of automation while maintaining high performance, with an emphasis on NASA mission critical applications. During the design of the 21st century launch infrastructure at Kennedy Space Center (KSC) different designs had to be evaluated in a very short period of time. Therefore, a solver was needed that provides a wide range capabilities that can be

3

used to study the pressure environment (ignition over-pressure and duct over-pressure) for highly complex geometries, conjugate heat transfer on the main flame deflector surface, and jet impingement noise. In this paper, the main capabilities of the LAVA solver will be demonstrated for relatively large scale applications. Basic verification and validation studies of LAVA were previously reported by Moini-Yekta et al.[67], and are continuously applied throughout the development and incorporation of additional features. In two companion papers that utilize the LAVA solver, an analysis and overview of higher-order discretization schemes for Cartesian grids are presented by Brehm et al.[18], and a study of gradient calculation methods for unstructured meshes is presented by Sozer et al.[93]. The organization of the paper is as follows. Section 2 presents the core algorithms used for each grid type in LAVA. Descriptions of the spatial and temporal discretization methods, overset techniques, and auxiliary modules are provided. Section 3 showcases several high-impact, NASA applications where the LAVA solver has been recently utilized, including results from the launch environment, ascent, and vehicle aerodynamics. Section 4 summarizes the paper and addresses the near and long term plans for the solver.

Figure 2: Schematic of LAVA infrastructure design illustrating current and future features.

4

2. Technical Approach The LAVA solver is highly flexible with respect to the computational mesh, and supports blockstructured Cartesian grids with Adaptive Mesh Refinement (AMR) and immersed-boundary capabilities, structured curvilinear overset, unstructured arbitrary polyhedral, and overset grid coupling. For all grid types, the unsteady compressible Navier-Stokes equations can be solved using a second-order, multi-species[45, 46], finite-volume formulation. Higher-order conservative finitedifference discretizations are also supported for structured Cartesian or Curvilinear grids[19]. The Spalart-Allmaras (SA) [97] and Shear Stress Transport (SST) [63] turbulence models, as well as their Detached Eddy Simulation (DES) variants[95, 94], are also available. Implicit time-integration methods for time-accurate unsteady simulations are implemented using a dual time-stepping approach. Standard and strong stability preserving higher-order Runge-Kutta schemes are available for explicit time-integration. Several flux difference and flux vector splitting schemes, such as the modified Roe scheme with preconditioning[45, 46] and AUSMPW+[54], are available. Standard boundary conditions such as slip and no-slip walls with adiabatic/isothermal/CHT conditions, periodic, characteristic based inflow/outflow, etc. are implemented. The overall infrastructure of the LAVA solver is illustrated in Figure 2. LAVA was developed in an object oriented framework and utilizes C++ and Fortran with MPI protocols and OpenMP for parallel processing. Computational performance is obtained using domain decomposition techniques and shared data between coupled modules. In Figure 2, current capabilities are shown in green, and solid lines indicate modules that have been fully integrated within LAVA. Future multi-physics capabilities are highlighted in yellow. The current LAVA infrastructure supports three grid types, i.e., block-structured Cartesian AMR, unstructured arbitrary polyhedral, and structured curvilinear. Overset coupling is available within the structured curvilinear solver module and the between blockstructured Cartesian and the unstructured arbitrary polyhedral solver modules. Moving boundary capabilities are available for 6-DOF body motion as well as for fully coupled fluid structure interaction within the Cartesian solver. The multi-physics capabilities include far-field acoustic solvers, conjugate heat transfer, as well as a structural finite element solver. The key focus of current developments is to include additional multi-physics capabilities, such as chemical reactions. All these capabilities are seamlessly available through a fully MPI-parallel C++ interface. More detailed information about the available modules are discussed in the following sections.

5

2.1. Governing Equations The conservation of mass, momentum, and energy of a viscous, N component mixture of arbitrary fluids with N − 1 individual species conservation equations can be written in vector form as:   ∂W  −F v = 0, +∇· F ∂t     v = F 1 , F 2 , F 3 .  and F F = F 1, F 2, F 3 v v v

(1)

where the conservative variable vector W = (ρ, ρu, ρv, ρw, ρH −p, ρY1 , · · · , ρYN −1 )T and the inviscid fluxes are:

⎡ ρu ⎢ ⎢ ⎢ ρu2 + p ⎢ ⎢ ⎢ ⎢ ρuv ⎢ ⎢ ⎢ ρuw ⎢ 1 F =⎢ ⎢ ρuH ⎢ ⎢ ⎢ ⎢ ρuY1 ⎢ ⎢ .. ⎢ . ⎢ ⎣ ρuYN −1

⎤

⎤

⎡

ρv ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ρvu ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ρv 2 + p ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 2 ⎢ ρvw ⎥,F = ⎢ ⎥ ⎢ ρvH ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ρvY1 ⎥ ⎢ ⎥ ⎢ .. ⎥ ⎢ . ⎥ ⎢ ⎦ ⎣ ρvYN −1

⎡

⎤

ρw ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ρwu ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ρwv ⎥ ⎢ ⎥ ⎢ 2 ⎥ ⎢ ⎥ 3 ⎢ ρw + p ⎥,F = ⎢ ⎥ ⎢ ρwH ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ρwY1 ⎥ ⎢ ⎥ ⎢ .. ⎥ ⎢ . ⎥ ⎢ ⎦ ⎣ ρwYN −1

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

where H = h + u · u/2. Here ρ is the density, (u, v, w) are the velocity components, p is the pressure, H is the total enthalpy and Yi are species mass fractions. The last species maintains conservation: YN = 1 −

N −1 i=1

6

Yi .

(2)

The viscous fluxes are (species mass diffusion is neglected here): ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 1 Fv = ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎤

⎡

0 ⎥ ⎢ 0 ⎥ ⎢ ⎢ τ τxx ⎥ ⎥ ⎢ yx ⎥ ⎢ ⎥ ⎢ ⎢ τyy τxy ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ τxz ⎥ 2 ⎢ ⎢ τyz ⎥ , Fv = ⎢ ⎢ f2 fv1 ⎥ ⎥ ⎢ v ⎥ ⎢ ⎥ ⎢ ⎢ 0 0 ⎥ ⎥ ⎢ ⎢ . .. ⎥ ⎥ ⎢ .. . ⎥ ⎢ ⎦ ⎣ 0 0

⎤

⎡

⎢ 0 ⎥ ⎢ ⎥ ⎢ τ ⎥ ⎢ zx ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ τzy ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ 3 ⎢ τzz ⎥ , Fv = ⎢ ⎢ f3 ⎥ ⎢ v ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 0 ⎥ ⎢ ⎥ ⎢ . ⎥ ⎢ .. ⎥ ⎢ ⎥ ⎣ ⎦ 0

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

where the viscous stress tensor is: 

τ = μ ∇u + (∇u)T + λI∇ · u

(3)

or fv1

=

fv2

=

fv3

=

∂T , ∂x ∂T , uτyx + vτyy + wτyz + κ ∂y ∂T , uτzx + vτzy + wτzz + κ ∂z uτxx + vτxy + wτxz + κ

where κ is the thermal conductivity coefficient. It is common, but not required, to use Stoke’s hypothesis; λ = −2μ/3. The solver supports arbitrary user-defined equations of state in the form:

ρ = f (p, T )

and

h = f (p, T ) .

(4)

2.2. Temporal Discretizations LAVA employs two main time-stepping schemes: explicit and implicit with a dual-time formulation. Explicit schemes benefit from the simplicity and performance of not requiring a left-hand-side (LHS) linearization, while they are restricted by CFL limits on time step size. Implicit schemes can significantly enhance convergence on meshes that contain high-aspect ratio cells (encountered with body-fitted meshes) or large variations in cell length scales. Implicit schemes may also be advan7

tageous when very large temporal scales are of interest, for example during the initial ascent of a rocket. Both approaches are described in this subsection, with a greater emphasis on the implicit scheme. The explicit time stepping scheme uses the generalized Shu-Osher formulation[87]: W (0) = W n , W (i) =

i−1 

 αi,k W (k) + Δtβi,k ∇ · F (W (k) ) ,

(5)

k=0

W n+1 = W (m) , which allows the user to specify any explicit RK scheme through the input file. This generalized formulation supports the standard Runge-Kutta (RK) schemes (1,2,4), as well as the Strong Stability Preserving RK (SSPRK) family[40]. Optimally smoothed multi-stage schemes [104] can be effectively used for steady-state convergence. The implicit time discretization described below follows Merkle and Choi [64] and Housman et al. [48]. Given the primitive variable vector, Q = (p, u, v, w, T, Y1 , · · · , YN −1 )T , the pseudo-time equations are as follows: Γp

  ∂Q ∂W  −F v = 0, + +∇· F ∂s ∂t

where W is the conservative variable vector and

∂Q ∂s

(6)

= 0 when Equation (1) is discretely satisfied.

Discretizing Equation (6) backwards in physical time yields:

Γp

∂Q W n+1 − W n + + ∂s

ψ n+1 2 (W

− 2W n + W n−1 )

Δt

n+1   −F v +∇· F = 0,

(7)

where ψ = 0 and ψ = 1 for first- and second-order implicit time differencing, respectively. Equation (7) can be compactly expressed as: Γp

∂Q = R(Q), ∂s

(8)

where the residual R(Q) contains the discrete deviation from Equation (1). The system of equations given in Equation (8) is marched in pseudo-time until a solution satisfying R(Q) ≈ 0 is found. Once R(Q) ≈ 0 is achieved, the solution satisfies Equation (1) at time nΔt. The inviscid contributions are computed using a first-order approximation. Only the face-adjacent contributions are included in the viscous Jacobian. It is crucial to study the effect of sub-iteration convergence and avoid error

8

propagation due to insufficient convergence in each time-step. The linearized system of equations, results in a block-septa-diagonal (in 3D), banded matrix (A) for structured grids and a blocked sparse matrix for unstructured grids. At pseudo-time m, the system of equations is: Am ΔQm = Rm .

(9)

Equation (9) is then solved using one of the available linear solvers (structured meshes: point and line Jacobi and Symmetric Gauss Seidel relaxation schemes as well as red-black variants; unstructured meshes: GMRES using PETSc [5, 6] or point Gauss-Seidel). For Cartesian meshes, grid cycling (successively adding levels) is utilized to accelerate convergence for steady-state simulations. For LES-type simulations, the time-step is chosen based on accuracy requirements for resolving the expected physical scales in the unsteady flow region, such as wakes and shear-layers. The time-step based on the explicit CFL restriction is usually orders of magnitude lower. The choice between explicit and implicit time-integration is up to the user and must be made on a case by case basis. 2.3. Spatial Discretization 2.3.1. Second-Order Discretization The finite-volume discretizations are formulated in divergence form using discrete fluxes at face centers as: ∇ · F ≈ ≈

1 V

 ∇ · F dV = V Nf aces

1 V

 F · ndA ∂V

1  Fi (QL , QR ) · ni Ai V i=1

(10)

where V is the cell volume, Ai is the face area and ni is the face normal vector. Following a secondorder MUSCL[102, 103] procedure, primitive variables at face centers are reconstructed from the left and right side of the face as:

QL = Qc,L + rL · ∇Qc,L

(11)

QR = Qc,R + rR · ∇Qc,R

(12)

where Qc and ∇Qc are the primitive values and their gradients at the cell center, and r is the position vector from cell center to face center. Gradient terms used in the MUSCL extrapolation can be 9

limited for non-smooth problems with an appropriate limiter function. Within LAVA, typically the minmod limiter is used. Convective numerical fluxes F at face centers are evaluated using either flux difference or flux vector splitting methods. The flux difference splitting evaluation of the numerical flux F is given by:  1  F (QR ) + F (QL ) − Γp |Γ−1 F (QL , QR ) = p A| (QR − QL ) , 2

(13)

where Γp = ∂W/∂Q is the change of variable matrix from conservative to primitive variables and A = ∂ F /∂Q is the convective flux Jacobian. The overbar represents evaluation of the matrices Γp −1 and |Γ−1 p A| at the density weighted face averaged state. The absolute value of the matrix |Γp A|

is given by R|Λ|R−1 where R and R−1 are the right and left eigenvectors of A, and |Λ| is the absolute value of the components of the eigenvalues of A. At low-speeds Γp becomes the low-Mach preconditioning matrix as described in [45, 46]. For flux vector splitting the numerical flux function is evaluated using AUSMPW+[54]. Diffusive fluxes are evaluated in a second-order accurate manner. 2.3.2. Higher-Order Discretization Several higher-order shock capturing schemes are available in the LAVA framework. The different spatial discretization schemes are discussed and analyzed in Brehm et al.[19] with respect to their accuracy and computational efficiency. Three classes of higher-order shock capturing schemes can be used: (1) central finite-difference schemes with explicit artificial dissipation, (2) compact centered finite-difference schemes with localized artificial diffusivity and (3) weighted essentially non-oscillatory (WENO) schemes in both explicit and compact finite difference forms. The WENO schemes are mostly utilized due to their increased robustness while being slightly more expensive than the two other types of shock capturing schemes schemes. For Cartesian and curvilinear the WENO schemes are implemented in conservative finite difference form. In the finite difference WENO scheme, the reconstruction is based on node values, while finite volume WENO schemes use cell averages for the interpolation. As shown in [87] the finite volume and the finite difference approaches can be linked by implicitly defining the primitive function h(x) of the flux f (x) as f (x) =

1 Δx



x+Δx/2

h(ξ)dξ, x−Δx/2

10

(14)

so that the spatial derivative ∂f /∂x is exactly defined by a conservative finite difference formula   ∂f 1 1 ˆ = (hi+1/2 − hi−1/2 ) = (fi+1/2 − fˆi−1/2 ) + O Δx2n−1 , ∂x Δx Δx

(15)

where n is the number of candidate stencils. Therefore, the numerical flux fˆi±1/2 should approximate hi±1/2 . In the LAVA solver, a distinction between the Cartesian and curvilinear implementation of the WENO schemes is made here. On Cartesian meshes, the classical flux reconstruction scheme is used following Equations 14 and 15. Generally, this approach only allows the use of flux vector splitting schemes. To obtain free-stream preservation on curvilinear meshes, variable interpolation is utilized. This approach follows the method outlined in Deng et al.[34]. To ensure higher-order accuracy with this approach the conservative finite difference formula in Equation 15 must be modified to higher-order accurate centered finite difference formula. The basic idea of WENO schemes is to combine lower order candidate stencils for the interpolation to the midpoint i + 1/2 to achieve a higher-order approximation in smooth flow regions and an essentially non-oscillatory interpolation at discontinuities[87, 58, 52, 7, 72]. In what follows, we will introduce the WENO scheme in a general form which applies to variable interface interpolation as well as interface flux reconstruction. Depending on the choice of the numerical scheme the actual stencil coefficients and weights will be different (see Nonomura et al.[75]). The variable φ is used to represent either, the variable or the flux. When referring specifically to a variable or flux, q and f ˆ at the interface can be written as are used instead of the placeholder φ. The final variable/flux, φ, φˆi+1/2 =

n−1

ωk φki+1/2 ,

(16)

k=0

where n is the number of candidate stencils, φki+1/2 is the interface variable/flux based on the k th candidate stencil, and ωk are convex weight coefficients. Choosing the optimal upwind biased weight coefficients ωk = ck for all k results in a (2n − 1)th -order accurate numerical approximation, where n is the number of candidate stencils. The WENO scheme only obtains these optimal coefficients in smooth flow regions. The idea is to use a smoothness indicator which automatically determines if the flow is smooth or contains large gradients/discontinuities. The smoothness indicator applied to

11

each candidate stencil essentially scales the optimal weights. An interim quantity αk is defined as

αk =

ck p, (βk + )

(17)

where βk is the smoothness indicator of the k th stencil. The power p, here only p = 1 and 2, in the denominator is used to achieve fast convergence to zero in non-smooth flow regions. To avoid division by zero (or values close to zero) a small number  is added to the denominator (typically  = 10−6 )[86]. To keep the weights αk in Equation (17) convex, they need to be normalized to obtain the final weight coefficients αk ωk = n−1 k=0

αk

.

(18)

For Cartesian and curvilinear meshes an alternative for calculation of the nonlinear weights in Equation (17) is generally employed. Borges et al.[16] presented a different way of computing the nonlinear weights that attempts to mimic the characteristic of the WENO-M scheme[43] at a reduced cost of roughly 25% by eliminating the mapping procedure. In the current paper we will refer to Borges’ scheme as WENO-Z, where the modified smoothness indicator, β Z , is defined as βkZ =

βk +  , β k + τ5 + 

(19)

with τ5 = |β2 − β0 |. The nonlinear weights are αZ ωkZ = 2 k i=0

, Z

αi

αrZ =

  ck τ5 . = c 1 + k βk +  βkZ

(20)

At this point all ingredients for the WENO schemes have been introduced. The flux reconstruction schemes being considered are limited to flux vector splitting schemes in the form f = f + + f − , i.e., Lax-Friedrichs, van Leer flux splitting, AUSM flux family[54], etc. In summary, the WENO flux reconstruction approach proceeds as follows: 1. Reconstruct fluxes at the midpoints to obtain f ± employing the following modified weights

12

and smoothness indicators for the variable interpolation as follows: 13 2 (φi−2 − 2φi−1 + φi ) + 12 13 2 β1L = (φi−1 − 2φi + φi+1 ) + 12 13 2 β2L = (φi − 2φi+1 + φi+2 ) + 12 β0L =

1 2 (φi−2 − 4φi−1 + 3φi ) , 4 1 2 (φi−1 − φi+1 ) , and 4 1 2 (3φi − 4φi+1 + φi+2 ) , 4

1 7 11 0 fˆi+1/2 = fi−2 − fi−1 + fi , 3 6 6 1 5 1 1 fˆi+1/2 = − fi−1 + fi + fi+1 , and 6 6 3 1 5 1 2 fˆi+1/2 = fi + fi+1 − fi+2 . 3 6 6

(21a) (21b) (21c)

(21d) (21e) (21f)

2. Apply conservative finite difference formula in Equation (15) to total flux, fˆ = fˆ+ + fˆ− . In order to provide superior freestream and vortex preservation properties on generalized curvilinear meshes a WENO variable interpolation/extrapolation approach can be summarized follows: 1. Interpolate/extrapolate primitive variables to the midpoints to obtain q ± employing the following modified weights and smoothness indicators for the variable interpolation as follows: 1 2 2 (qi−2 − 4qi−1 + 3qi ) + (qi−2 − 2qi−1 + qi ) , 4 1 2 2 β1 = (qi−1 − qi+1 ) + (qi−1 − 2qi+1 ) , 4 1 2 2 β2 = (−3qi + 4qi+1 − qi+2 ) + (qi − 2qi+1 + qi+2 ) , 4

β0 =

3 5 15 qi−2 − qi−1 + qi , 8 4 8 1 3 3 = − qi−1 + qi + qi+1 , and 8 4 8 3 3 1 = qi + qi+1 − qi+2 . 8 4 8

(22a) (22b) (22c)

0 qˆi+1/2 =

(22d)

1 qˆi+1/2

(22e)

2 qˆi+1/2

(22f)

2. Compute fluxes at the midpoints using a numerical flux in the form f + (q + )+f − (q − ) or f (q + , q − ). 13

3. Compute flux derivatives at the nodes by applying the generalized higher-order finite-difference formula[73],    n  ∂f  ∂f  1 ˆ ∂f  + b + c ≈ αk fi+k+1/2 − fˆi−k−1/2 +    ∂x i−1 ∂x i ∂x i+1 Δx

(23)

k=0

1 βk (fi+k − fi−k ) , Δx n

k=1

where in the current implementation we simply used an explicit finite-difference scheme assuming a = c = 0, b = 1, and n = 3 (or n = 5) as in Nonomura et al.[73, 74]. Note that for the solution of the Euler and Navier-Stokes equations the scalar quantities q and f can simply be replaced with vector quantities, Q and F . Both variable interpolation or flux reconstructions can be applied to the variable/flux itself or the variable/flux in characteristic space:

Φ = Lφ,

(24)

where Φ represents the variable/flux in characteristic space and L is the matrix containing the left eigenvectors of the Jacobian matrix A = ∂F/∂W. Another variation of these schemes comes from choosing between standard or compact interpolation/differences. While additional variants of the WENO family are available in LAVA the presented schemes are utilized for most applications. 2.4. Mesh Paradigms 2.4.1. Block-Structured Cartesian Mesh For Cartesian meshes involving geometry, a block-structured, Immersed-Boundary (IB) AMR method is utilized. This methodology is capable of automatically generating, refining, and coarsening nested Cartesian volumes given a closed surface triangulation, and hence offers the ability to dynamically track important flow features as they develop (see Figures 3a and 3b). In adaptive methods, one adjusts the computational effort locally to maintain a uniform level of accuracy throughout the problem domain. Cartesian AMR is a proven methodology for multi-scale problems, with an extensive existing mathematical and software knowledge base [14, 2, 10, 11, 107]. The LAVA code incorporates data structures and inter-level AMR operators from the high-performance Chombo AMR library [31].

14

(a)

(b)

Figure 3: (a) Block-structured Cartesian AMR feature tracking for model rocket launch simulations. (b) Adaptive mesh refinement for a parachute simulation. Yellow boxes mark the finest level on the mesh and are quadrant clipped to expose the parachute and passive particles for visualization.

AMR allows the simulation of a wide range of spatial and temporal scales through local refinement. For this work, the constraint that all AMR levels are advanced with the same time step is imposed; i.e., sub-cycling the levels as is done in the literature[2, 61] has not been implemented. In the AMR research community, methods that do not sub-cycle the levels are termed composite methods. See Martin et al. [61] for issues associated with sub-cycling. Spatially refined regions are organized into rectangular patches. Our approach is to express the AMR discretizations in terms of the corresponding uniform grid discretizations at each level, followed by the appropriate inter-level operators. As in Barad et al.[9], interpolation is used to provide ghost cell values for points in the stencil extending outside of the grids at that level and inter-level operators to coarsen and refine data when regridding. Refinement and coarsening is based on either entropy adjoint[35] or gradient/vorticity based tagging (|∇q| or | ω |) where the threshold is automatically tuned to provide a user-requested number of cells, either per level or for the entire AMR hierarchy (e.g., see Figure 3). Refinement can be constrained to occur within specified regions (e.g., points, boxes, or arbitrary surfaces), and can be specified based on surface components and several other fixed geometric features/metrics. Two sharp immersed boundary methods (IBM) are available in LAVA-Cartesian. The IB representation makes it possible to easily model complex geometry without having to manually generate a volume grid. The first IBM is based on a ghost cell boundary treatment[65, 68] and the second higher-order IBM utilizes locally stabilized finite difference stencils[23, 24]. The second IBM has been

15

(a)

(b)

Figure 4: Ghost cell filling using two different IB methods: (a) full image reflection [65], (b) fixed distance reflection into fluid (indicated by dotted line). Image-Points (IP) are interpolated. IP data is combined with boundary condition and the resulting ghost cell values are computed. With (a) IP can get arbitrarily close to the ghost cell, and the interpolation stencil for IP can be dominated by the ghost cell. With (b) the fixed stand-off distance forces the interpolation stencil to have significant contributions from fluid cells.

utilized in combination with up-to 11th -order accurate upwind-biased finite difference schemes[42]. In the following discussion both methods are briefly described. The ghost cell method is based on the idea of filling cells that lie inside the geometry with appropriate flow data. Figures 4a and 4b schematically display the procedure for filling the ghost cells. For a second-order accurate scheme only the first layer of grid points beyond the boundary of the geometry needs to be filled and are referred to as ghost cells while the rest of the grid points inside the geometry are referred to as solid points. IB ghost cells are covered cells that are within the stencil of a fluid cell, and are filled using Image-Points (IP) and boundary conditions. Two variants for defining IP locations are available: exact reflection (Figure 4a) and fixed distance, typically |Δx|/2 (Figure 4b). An advantage of the fixed distance IP is that the implicit nature of the linear interpolation stencil makes it well conditioned (i.e., when geometry nearly coincides with the ghost cell center). Interpolation to the IP uses either linear or weighted least-squares interpolation. The interpolation can be posed as either implicit or explicit. If implicit ghost cell filling is used, the interpolation can use more than one ghost cell (as in Figure 4a, top right stencil), but iterations are required for ghost cell boundary conditions. If explicit ghost cell filling is used, neighboring IB ghost cells are excluded from the least-squares stencil, thus the procedure is not iterative. Note that a higher-order extension of this method is theoretically straightforward by extending sizes of interpolation and difference stencils. In practice, however, the robustness of the numerical scheme suffers severely when applying this extension.

16

While the ghost-cell method is relatively easy to implement, the second IBM is higher-order accurate and does not require ghost cells. In this method, the irregular one-sided finite-difference stencils carry a free parameter that is adjusted locally to improve the stability characteristics of the finite difference scheme. The irregular finite difference stencils are treated isolated from the rest of the computational domain as shown in Figure 5c. The validity of isolating the local grid stencils from the rest of the computational domain has previously been demonstrated in Brehm et al.[23, 24]. Figure 5a shows the irregular finite difference stencil in the vicinity of the immersed boundary with a normalized distance to the grid line intersection point of ψ. The dependence on the stability of the update matrix1 on the free stencil coefficient α˜1 is illustrated in Figure 5a for the scalar advectiondiffusion equation considering the fourth-order accurate Runge-Kutta time integration scheme, two different boundary distances (ψ = 0.1 and ψ = 0.9), and a third-order accurate immersed boundary treatment. The boundary conditions are provided at grid line intersection points and they are obtained by weighted least-squares interpolation/extrapolation stencils. For slip wall boundary conditions the velocity field is also extrapolated towards the grid line intersection points while imposing the no-penetration condition (v · n = 0) at the solid wall. ρ(A) for ψ=0.9 ρ(A) for ψ=0.1 ρ(Aid) ∼ calculated α1 for ψ=0.1 ∼ for ψ=0.9 calculated α

1.1 1

ρ(A)

1

0.9 0.8 0.7 0

40

80

∼ α 1

(a)

120

160

200

(b)

(c)

Figure 5: (a) Dependence of spectral radius on free stencil coefficient α ˜ 1 . (b) Localization illustrated in 2D for circle immersed into Cartesian grid. (c) Computational setup in the vicinity of the immersed boundary including irregular grid points (), points in fluid (), boundary points () and points outside the computational domain are dropped.

Some of the key challenges for simulating extremely complex geometries, such as the launch site or a contra-rotating open rotor system[21], is the treatment of the thin geometry, accounting for the thin boundary layer near solid walls, and dealing with relative body motion. Several fast algorithms for discrete geometry manipulation are required for an efficient IB implementation. Specifically, geometric queries from the closed surface representation, including inside-outside, and signed distance 1 The

update matrix M updates the solution φ from time level tn to tn+1 , i.e. φn+1 = M φn .

17

irregular pnt. cloud 1 pnt. cloud 2 pnt. cloud 3 pnt.

(a)

pierce pnt. x−rays in x−dir. x−rays in y−dir.

(b)

(c)

Figure 6: (a) Cloud selection process for the boundary extrapolation stencil, (b) illustration of x-ray algorithm used for in-out testing and (c) freshly cleared cells after moving the immersed boundary in the time interval [tn ,tn+1 ].

(exact and approximate) functions. For simple analytic geometries, a Constructive Solid Geometry (CSG) methodology [57] is exact and efficient. For complex geometries where CSG is not available, geometries are represented with closed surface triangulations. In what follows, closed surface triangulations representing the surface (as CSG is relatively trivial) are assumed. The treatment of the thin geometries only affects the cloud selection process for the boundary condition extrapolation stencil since the current sharp immersed boundary method does not utilize ghost cells. The algorithm ensures that points are only selected for the point cloud if they reside on the same side of the boundary as the irregular grid point. Figure 6a schematically illustrates the cloud selection process at an irregular grid point. First, all irregular grid points compute their base clouds while the base clouds for the regular grid points are assumed to contain the direct neighbors, i.e., (i ± 1,j,k), (i,j ± 1,k), and (i,j,k ± 1). Then the point cloud at the irregular grid point is gradually extended by adding the base clouds from the valid neighbors. There are two additional issues that had to be addressed for complex geometry in particular for moving geometries: (1) speeding up the geometry queries, (2) regeneration of irregular stencils after moving the geometry, and (3) the treatment of Freshly Cleared Cells (FCC). An x-ray algorithm, as visualized in Figure 6b, is used for inside-outside (in-out) testing. Parallel geometry kernels were implemented to allow fast in-out testing and exact computation of the distance from the irregular grid points to the surface triangulation (including point to plane and point to edge cases). An optimized Boundary Volume Hierarchy (BVH) based algorithm was developed for this purpose. Figure 6c displays some FCCs after moving the immersed boundary in the time interval [tn ,tn+1 ]. Since the FCCs do not contain a valid time history their state vector is initialized considering valid grid points (not including other FCC) in the neighborhood of the FCC. Some more advanced strategies for updating the FCCs similar to what has been discussed in Brehm and Fasel[22] for the incompressible Navier-Stokes equations

18

are currently being explored. Distance function evaluation is undertaken with two main approaches: exact and approximate (both signed using an inside-outside test). The exact distance function is computed using a tree data-structure to locate the nearest candidate surface triangles, which are then queried for exact distances (including point-to-plane and point-to-edge cases). Exact distance calculation is expensive, therefore approximate distances are used when adequate (i.e., for turbulence models). Approximate distances on Cartesian grids are computed using a signed (i.e., in-out) version of the fast sweeping method[108]. One of the key challenges for immersed boundary methods is the treatment of viscous walls at high Reynolds numbers. For high-Reynolds number viscous flow simulations, current immersed boundary approaches are inefficient in resolving the viscous boundary layers since the Cartesian mesh layout generally does not allow the use of (wall-normal) high aspect ratio cells in the vicinity of the wall. This issue has been addressed in several research investigations mainly by employing wall models near the immersed boundary[82]. To incorporate the effect of the thin boundary layers a wall modeling approach was chosen. The wall modeling approach was first successfully used to simulate the flow around the partially dressed landing gear for the Third AIAA Workshop on Benchmark problems for Airframe Noise Computations (BANC-III). A vorticity visualization of the wake flow behind the landing gear is shown in Figure 7. Figure 7a shows the effect of the wall velocity boundary condition on the pressure distribution around the star board wheel of the landing gear used in the AIAA BANC-III workshop. The pressure distribution for no-slip wall and slip wall boundary conditions do not match the experimental data (BART and UFAFF) very well due to predicting too early separation for no-slip and failing to predict any separation at all for slip. An improved match between the numerical simulation and the experimental data can be obtained with partial slip wall boundary condition. Figure 7b shows a comparison between PIV measurements and CFD data in a cut-plane through the starboard wheel. A good comparison between the experimental and CFD results was obtained. 2.4.2. Unstructured Arbitrary Polyhedral Mesh Resolution of high Reynolds number boundary layers with the isotropic Cartesian mesh approach can be prohibitively expensive, or requires the use of wall modeling approaches that are still an ongoing topic of fundamental research in the context of immersed boundary methods. For capturing

19

(a)

(b)

(c) Figure 7: (a) Comparison of pressure distribution around starboard wheel of the landing gear used in the AIAA BANC-III workshop considering slip wall, no-slip wall, and wall modeling results. The data is compared to two sets of experimental data (BART and UFAFF). (b) Comparison of u-velocity, v-velocity, vorticity, and turbulent kinetic energy in PIV plane through wheels of landing gear. (c) Vorticity flow visualization in the wake of the partially dressed landing gear.

20

anisotropies inherent in fluid dynamics applications, body-fitted high-aspect ratio cells (often in the range of 105 ) are advantageous. For this purpose, unstructured and structured curvilinear methods are available within LAVA.






 

(a)

  





















































 

 







  




(b)

(c)

Figure 8: Different methods for gradient calculation via (a) Green-Gauss method, (b) weighted least squares stencil, and (c) curvilinear method.

LAVA’s unstructured solver offers the ability to handle arbitrary polyhedral cells (including standard cell types such as tetrahedral or hexahedral) as well as high-aspect ratio polygonal prism layers near wall boundaries. Parallel domain decomposition of unstructured meshes in LAVA is achieved via ParMETIS[53]. Different approaches for gradient calculation are available for unstructured grids. The difficulty in calculating gradients in an unstructured mesh stems from the lack of a consistent, inherent connectivity. The stencil for gradient calculation, as well as the corresponding coefficients vary cell-by-cell and are costly to compute. Hence, those are typically pre-computed and stored. Two of the most common methods for gradient calculation are the Green-Gauss and the least squares approaches (see Figure 8a and 8b). For the Green-Gauss gradient method, it is assumed that valid data on the face centroids are available to compute the cell gradients. Nf aces 1 ¯ ∇φ = ˆ f Af , φf n V

(25)

f =1

where Nf aces is the number of faces and φ¯f is the average of the scalar over the face f . Different choices for face averaging methods are: simple face averaging, inverse distance weighted face interpolation, weighted least squares face interpolation, and weighted tri-linear face interpolation. In the least squares gradient calculation, an overdetermined system of equations is solved in terms of the cell values φi included in the WLSQR stencil and the computed values at targeted cell centroid φc ,

21

and the derivatives φc,x and φc,y . s  f = Bφ AΦ ⎡



wi

 s = [φ1 , φ2 , φ3 , ...]T ,  s , with Φ  c = A−1 B φ  c = [φc , φc,x , φc,y ]T , φ Φ

→

⎢ ⎢ A=⎢ ⎢ wi Δxi ⎣ wi Δyi

  

wi Δxi wi Δx2i wi Δxi Δyi



⎤

⎡

wi Δyi

⎥ ⎥ wi Δxi Δyi ⎥ ⎥ , and ⎦  wi Δyi2 

B

⎢ w1 ⎢ =⎢ ⎢w1 Δxi ⎣ w1 Δyi

w2 w2 Δxi w2 Δyi

⎤ ...⎥ ⎥ ...⎥ ⎥. ⎦ ...

(26)

(27)

For the least squares weights wi = 1/|di |p is used and the distance vector is defined as di = ri − r0 . In addition to direct computation of gradients employing least squares stencils including neighboring grid cells, a curvilinear gradient method is also available (see Figure 8c). A selection procedure evaluates every possible combination of stencil point pairs to pick curvilinear directions ξ and η. Once the stencil is selected, central differencing yields ∂φ φi+1 − φi−1 = , ∂ξ 2Δξ

and

∂φ φj+1 − φj−1 = ∂η 2Δη

with Δξ = Δη = 1.

(28)

The gradient in Cartesian coordinates can then be arrived as ⎧ ⎫ ⎪ ⎪ ⎨ ∂φ ⎬ ∂x ⎪ ⎩ ∂φ ⎪ ⎭ ∂y

= J −1

⎧ ⎫ ⎪ ⎪ ⎨ ∂φ ⎬ ∂ξ ⎪ ⎩ ∂φ ⎪ ⎭ ∂η

⎤

⎡ , with

∂x ⎢ ∂ξ

J =⎣

∂x ∂η

∂y ∂ξ ⎥

⎡

1 ⎢ (xi+1 − xi−1 ) ⎦= ⎣ 2 ∂y (xj+1 − xj−1 ) ∂η

⎤ (yi+1 − yi−1 ) ⎥ ⎦ (yj+1 − yj−1 )

(29)

In the current implementation, the stencil pairs that produce the largest determinant of the Jacobian |J | are selected out of all the possible options. Considering that only a compact neighborhood consisting of face sharing cells are considered for the stencil selection. The largest Jacobian determinant criteria usually identifies the most orthogonal ξ and η directions. In Sozer et al.[93], the different methods for calculating gradients on unstructured mesh was tested for the Euler equations. No clear best method emerged in this study but strengths and shortcomings of the methods for different cell types were exposed. The current best practice is to utilize the curvilinear scheme in prism layers and the Green-Gauss method with LSQR face interpolation on arbitrary grid types. Note that the viscous terms are discretized by utilizing the WLSQR face interpolation procedure and subsequently applying the Green-Gauss formula to the viscous fluxes. For the implicit time integration scheme, i.e., BDF1 and BDF2, only first-order terms are considered in the linearization for the left hand side operator. The generalized minimal residual (GMRES) method is most com-

22

monly used for steady flow problems to solve the linear system, while other iterative schemes are also available. 2.4.3. Structured Curvilinear Mesh For capturing anisotropies inherent in fluid dynamics applications such as boundary layers, shocks, wakes and shear layers, body-fitted structured curvilinear grids are advantageous. Utilizing curvilinear (and Cartesian) grids has many computational benefits, including low memory requirements and efficient cache utilization. When structured curvilinear meshes are necessary, an overset grid methodology is adopted to reduce the grid generation constraints for complex geometries. In order to solve the governing equations on a curvilinear grid, the equations must be transformed before discretization. Below is the set of governing equations, including the pseudo-time derivative term, transformed from Cartesian to a non-orthogonal curvilinear coordinate system,   ˆi − Fˆ v 3 ∂ F ˆ ˆ i ∂Q ∂W Γp + + = 0, ∂s ∂t ∂ξ i i=1

with

ˆ = J −1 Q, W ˆ = J −1 W, Q

(30)

where s is pseudo-time, and the convective and viscous fluxes in curvilinear coordinates are 3 3 ∂ ξˆi ∂ ξˆi ∂ ξˆi + Fk , and Fˆiv = Fkv , Fˆi = ∂t ∂xk ∂xk k=1

with

k=1

∂ξi ∂ ξˆi ∂ξi ∂ ξˆi = J −1 , = J −1 . ∂t ∂t ∂xk ∂xk

(31)

The curvilinear transformation metric terms, ∂ξi /∂xk and J −1 , are evaluated using either the secondorder free-stream preserving method described in Vinokur [105] or higher-order centered difference formulas as discussed in Deng et al.[34]. The key for obtaining free-stream preservation for higherorder accurate schemes is that the same interpolation and difference operators have to be used for the metric terms as for the discretization of the field equations. The biggest advantange of the curvilinear solver is that it allows the use of efficient and well-established higher-order accurate schemes. For higher-order unstructured schemes the construction of the left hand side for implicit time integration schemes can be very memory and time expensive. This issue has been addressed by various researchers but has not reached a level of maturity where it is production ready yet. Hence, the curvilinear solver is currently the method of choice within LAVA for high-fidelity hybrid RANS/LES-type simulations with viscous wall treatment.

23

2.4.4. Overset Cartesian/Unstructured Coupling The LAVA solver features an overset scheme that inherits the benefits of the block-structured Cartesian mesh (for the off-body) and the body-fitted unstructured mesh (for the near-body). The AMR Cartesian background mesh provides an optimal data structure to capture flow physics and incorporate grid adaption techniques. The unstructured near-body approach reduces the grid generation time and the level of user expertise associated with the generation of viscous grids. The hybrid overset scheme provides the capability to use an IB Cartesian scheme to model the geometry, where boundary layer resolution is not required. To reduce user inputs and improve connectivity, the Cartesian grid is automatically refined to the representative length scale of the local near-body grid at the overlap region. These features are particularly advantageous for simulations that involve complex geometries. The overset coupling approach between the block-structured Cartesian mesh and the bodyfitted unstructured mesh is illustrated in Figures 9a-f. Overset connectivity methods are used for

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9: Overset coupling approach between the block-structured Cartesian mesh and the body-fitted unstructured mesh: (a) Airfoil geometry immersed into Cartesian mesh, (b) hole cutting of Cartesian mesh, (c) generation of unstructed body-fitted mesh, (d) hole cutting of unstructured mesh, (e) identification of fringe points on unstructured mesh and donor cells on Cartesian mesh, and (e) identification of fringe points on Cartesian mesh and donor cells on unstructured mesh.

24

communication between the Cartesian off-body and the unstructured near-body grids. The overset coupling approach is outlined for a simple 2D airfoil, while a more complicated application of this approach is provided in Section 3.2.3. The coupled simulation starts with immersing an arbitrary complex geometry into a Cartesian mesh (see Figure 9a). The Cartesian mesh is adjusted to a grid resolution that would be conducive for Euler-type solution. After the initial Cartesian grid is generated a minimum distance function is computed and all grid cells are blanked that fall inside the geometry plus an explicit hole cutting distance (see Figure 9b). Next, the body-fitted unstructured mesh is generated with prism layers generated at viscous wall boundaries (see Figure 9c). After a certain distance the prism layers transition into arbitrary polyhedral cells (this is a function of the mesh generator). Mixed mesh types can also be utilized to generate a valid mesh topology in cracks or other delicate geometric irregularities. After the initial body-fitted unstructured grid is generated an explicit hole cutting approach is utilized to remove dispensable unstructured cells that lie outside the hole cutting distance (see Figure 9d). An approximate distance function computation is used to reduce computational expense since the accuracy of the distance function is not very crucial for this process. Once step 1-4 have been completed the most difficult aspects of the coupling approach remain, namely the identification of fringe and donor cells on the body-fitted unstructured mesh and Cartesian mesh. In the connectivity scheme, the cells at the outer boundary of the unstructured near-body grid are marked as fringe (see Figure 9e). Trilinear interpolation is used to fill the unstructured fringe cells from the corresponding Cartesian donor dual cell. The Cartesian donor cells can be simply identified by determining the Cartesian cell encompassing the unstructured fringe point. This process is very efficient and straightforward due to the underlying Cartesian data layout. Cartesian fringe cells are determined by marking cells whose discretization stencil contains a hole point (see Figure 9f). An efficient kd-tree search algorithm is utilized to find the closest unstructured cell centroid and the associated cell neighbors to interpolate data to the Cartesian fringe cells using least squares. Note stencil walking strategies are employed to identify the best possible unstructured donor stencils. Parallelization has been accomplished using domain decomposition and MPI message passing protocols with custom communicators for improved performance. Although the LAVA solver currently supports the coupling of unstructured near-body and structured Cartesian off-body grids, additional configurations will be pursued. The flexibility of the hybrid overset scheme can be extended to incorporate different combinations of structured Cartesian, structured curvilinear and

25

unstructured grids. 2.5. Wall Modeled Large Eddy Simulation The LAVA solver supports various turbulence models, among the most common ones are the SST and SA models. For applications where the unsteadiness in the flow must be modeled with higher fidelity, because it is inherent to the problem in hand, LES-type capabilities are available. Within the realm of LES-type approaches, the implicit large eddy simulation approach may be considered as the most simplistic strategy among them. In contrast to classical LES, which employs an explicitly computed sub-grid scale (SGS) closure, ILES relies on the inherent regularization mechanism through the truncation error of the convective fluxes, which is also referred to as the implicit SGS model ([50]). The primary role of the SGS model is to remove the energy that would cascade down to the unresolved scales. Hence, the numerical scheme must be chosen carefully to provide a physically consistent implicit SGS model. Within LAVA, the modified sixth-order accurate WENO scheme ([51]) is commonly used due to its superior, physically motivated scale separation properties. Hu and Adams[50] have demonstrated that the modified sixth-order accurate WENO scheme reproduces the Kolmogorov range of the turbulent kinetic-energy spectrum at the limit of infinite Reynolds number independent of grid resolution, while maintaining the shock-capturing capabilities of the original WENO-CU6. It is important to point out that our experience, as well as the results from other research studies (see [27] and [89, 88]), confirms that for simulations of turbulent flows away from solid boundaries, the ILES approach is a well-suited simulation strategy. Moreover, the ILES approach is very efficient in utilizing the available grid resolution for capturing the shear layer transition process and modeling salient features of the turbulent flow. Preliminary studies and other researchers ([90, 89]) have even demonstrated that the transition of the shear layer can be delayed and the proper amount of turbulence generation was not generated as well when utilizing a sub-optimal explicit SGS model for these flow conditions. The ILES approach has been successfully applied within LAVA for various shear-layer dominated flow fields [47, 49, 8, 26]. In Brehm et al.[25], the computational simulation strategy with respect to nozzle inflow conditions and the sub-grid-scale modeling is closely aligned with the recommendations of [89, 88], who discuss and compared different strategies for turbulent jet flow simulations. They demonstrated that their simulation approach to predict the noise field is very accurate and computationally efficient. To address the need for a physically motivated SGS model LAVA provides the Detached Eddy

26

Simulation (DES) [95, 99] and Delayed DES [94, 91, 96] (DDES) turbulence model closures that are well-tested hybrid RANS/LES models for highly separated flows. In the original DES model, the transition between RANS and LES models was based strictly on local mesh size relative to the wall-distance. For geometries with a wide range of geometric length scales, such as a high-lift device with finite-thickness leading and trailing edges, the local mesh spacing may become small (but not small enough) to force transition from the RANS model to the LES model in places where the model is not appropriate. This brought about the modification of the model denoted DDES which attempts to maintain RANS mode in the attached boundary layer. Inspection of the switching function often shows a strange behavior of going from RANS near the wall, to LES, back to RANS just past the edge of the boundary layer, and subsequently back to LES. An alternative strategy appropriate for structured multi-block and overset grids is the Zonal DES (ZDES) approach [32] in which specific zones are designated to use the RANS, DDES, or LES models explicitly. In Section 3.3.2, this approach was employed on curvilinear meshes for modeling slat noise. 2.6. Auxiliary Modules 2.6.1. Conjugate Heat Transfer LAVA has the capability to perform conjugate heat transfer simulations for unsteady surface temperature and heat transfer rate predictions. The heat conduction into the solid is modeled through local 1D rays that extend from the surface face centroids into the solid. Temperature dependent solid material thermal properties can be used where data is available. Fluid-solid thermal coupling is done at the sub-iteration level and the unsteady heat conduction equation is solved along each 1D ray. A graphical representation of the conjugate heat transfer method is illustrated in Figure 10. The surface temperature can be capped at the melting point of the solid material ac-

Figure 10: Schematic of conjugate heat transfer

counting for phase changes. Both isothermal and adia-

method between body-fitted-unstructured and

batic boundary conditions at the back of the material are available.

27

1D solid grids.

2.6.2. Aero-Acoustics High-fidelity simulation of acoustics is still a challenging task using CFD techniques, as demonstrated in Housman et al. [44]. For the computation of acoustic wave propagation and scattering over large distances (x >> λac ) even higher-order methods become computationally expensive. The dissipative and dispersive nature of these methods make them impractical for large-scale applications. Acoustic propagation can often be accurately predicted using more efficient, simplified linear models such as solving the wave equation or employing an acoustic analogy in terms of integral equations. The LAVA framework provides a hybrid CFD/CAA coupling approach. First, a high-fidelity unsteady CFD simulation is performed with an adequate time-step (relative to the frequencies of interest or physical scales) and a highly refined grid in noise generating regions of the computational domain. Time accurate solution histories are recorded on a permeable acoustic surface embedded within the CFD volume grid and enclosing the noise source. Two linear acoustic solvers are available for acoustic wave propagation away from the source region: (1) the Ffowcs-Williams and Hawkins (FW-H) method and (2) the equivalent source method (ESM). Both methods are currently employed in the frequency domain. The implementation of the FW-H method follows the description in [60] and is implemented utilizing MPI for parallelization. The ESM is utilized when acoustic scattering effects are expected to play an important role. The ESM method proceeds as follows: Once a significant amount of time-integration has been completed, the recorded time history is transformed into the frequency domain using the Fast Fourier Transform, and a linear Helmholtz Boundary Value Problem (BVP) is solved for each frequency independently using a highly scalable implementation of the Equivalent Source Method, see Ochmann [76]. An outline of the Helmhotz BVP solver is described in the following sections. At each frequency the following radiating/scattering Helmholtz BVP is posed: ∂P  = −ˆıkρ∞ c∞ u ˆ · n, ∂n

∇2 P  + k 2 P  = 0, ∀x ∈ Ω;     ∂P  + ˆıkP R =0 lim ∂R R=| x|→∞

∀x ∈ ∂Ω;

(32)

Sommerfeld radiation condition.

where P  is the acoustic pressure, k = 2πf /c∞ is the wave number, ρ∞ is the free-stream density and c∞ is the free-stream sound speed. On perfectly reflecting scattering surfaces, u ˆ · n = 0 is enforced. For radiating surfaces, u ˆ · n is the Fourier transform of u · n (time averaged velocity has been subtracted out) obtained from the CFD simulation. An alternative Dirichlet boundary

28

condition for the acoustic pressure on radiating surfaces is also available. In order to solve the Helmholtz BVP, a set of M monopole source locations are generated inside the radiating/scattering surface {xm }m=1,M where M < N , and N is the number of surface elements used to discretize the radiating/scattering surface. Using the free-space Green’s function G(x, xm ) =

x− xm | 1 e−ˆık| 4π | x− xm | ,

an approximation to the acoustic pressure field, p , is constructed such

that: 

p (x, t) =

M

! am G(x, xm ) eˆıωt ,

(33)

m=1

where ω is the angular frequency and {am }m=1,M are complex-valued coefficients that remain to be determined. By constructing p as a linear combination of free-space Green’s functions, the wave equation and the Sommerfeld radiation condition are satisfied exactly in the entire domain. The boundary condition on the radiating/scattering surfaces is the remaining part of the Helmholtz BVP that must be satisfied. In order to approximately satisfy the boundary condition, the following linear overdetermined system of equations is constructed for the unknown coefficients {am }m=1,M . ! ∂G(xn , xm ) = −ˆıkρ∞ c∞ [ˆ am u · n] (xn ) for n = 1, · · · , N . ∂n m=1 M

(34)

Solution of the overdetermined system is obtained by applying a parallel implementation of the complex-valued Conjugate-Gradient algorithm to the normal equations A∗ Ax = A∗ b where A = ∂G(xn , xm )/∂n, x = (a1 , · · · , aM ) , b = −ˆıkρ∞ c∞ ([ˆ u · n] (x1 ), · · · , [ˆ u · n] (xN )) , and A∗ is the T

T

Hermitian transpose of A. M

Once the coefficients {am }m=1 are found, the solution can be evaluated at any observer location x ∈ Ω. For a large number of observer locations (e.g., evaluation on a plane), a parallel version of the solution evaluation has been implemented. Once the solution is evaluated at each observer location for all frequencies, it can be transformed back into the time domain, which has also been implemented. A detailed example of the CFD/CAA coupling procedure is provided in the application of launch acoustics (see Section 3.1.3). 2.6.3. Fluid Structure Interaction In order to model the structural deformation of solids a linear structural finite-element solver was implemented in the LAVA framework. Two different element types are currently available within LAVA. For 2D calculations, we use the Bernoulli-Euler beam element that contains three

29

w2

v2

v2,Θy

u2,Θz u2,Θx

v1

y x

u1,Θz

h z

w3

w1

a

v1,Θy u1,Θx

b

(a)

v3,Θy u3,Θx

(b)

Figure 11: (a) Bernoulli-Euler beam for 2D calculations. (b) Triangular shell elements for 3D calculations.

degrees of freedom at each of the two nodes allowing for end displacements (u1 ,v1 ,u2 ,v2 ) and end rotations (Θz,1 , Θz,2 ). For 3D calculations, triangular shell elements are available containing five degrees of freedom per node. At each node of the triangular shell elements in-plane (un ,vn ) and out-of-plane (wn ,Θxn ,Θyn ) deformations are supported. The triangular elements are formulated in area coordinates to ensure that the normal slope varies linearly along an edge. This formulation will guarantee continuity of the lateral displacement and both its first derivatives between elements. The global system is dynamically solved for in terms of the displacement vector u in the well-known matrix vector form M

d2 u du + Ku = f , +D dt2 dt

(35)

where M is the mass matrix, D is the damping matrix, K is the stiffness matrix, and f is the load vector. Load vector external forces and moments can be prescribed by the user in addition to the pressure and viscous forces/moments which are supplied by the CFD solver. For time-integration the full mass matrices are considered and the system of equations are solved in parallel employing Krylov subspace methods. Commonly, the generalized minimal residual method (GMRES) is employed to iteratively solve the linear system at each time step while other iterative solver are also available. A second-order accurate time-integration scheme is utilized. For more details, the reader is referred to Brehm et al.[20]. Figure 12a displays the computational setup of a fully coupled fluid-structure-interaction test case. This test case involves a shock wave impacting on a deforming panel as first presented in Giordano et al.[38]. The objective of this study is to numerically predict the deformation of the cantilever panel motion and the coupled fluid dynamics. The characteristics of the panel were chosen to obtain a strong interaction between the structure and the fluid. A quantitative comparison of

30

W=0.295m

u L=109.7 m/s

E=200GPa ρ =7870kg/m3

T L=332.04 K

h=0.015m

p =154,000 Pa L t=1mm

uR =0 m/s p =100,000 Pa R

H=0.065m

x1 =0.035 m

L=0.04 & 0.05m

x1 =0.03 m

TR =293.04 K

x1 =0.045 m

Figure 12: (a) Simulation setup for deforming panel problem. (b) Comparison of horizontal tip deflection for L = 40mm.

the computational results against experimental data is provided in Figure 12b. The horizontal tip deflection is plotted over time. The horizontal tip deflection agrees well with the two numerical simulations while some discrepancy with the experiment can be observed. Figure 13 shows a comparison of Schlieren contours from experiments and numerical Schlieren contours by plotting the magnitude of the density gradient for the numerical simulation. The exact time of the first snapshot was estimated to closely match the key features of the flow field at t = 0. The LAVA simulation results closely match the key flow features that appear in both the experiment and computation provided by Giordano et al.[38]. It must be pointed out that the strong starting vortex dissipates much more quickly in experiments due to the presence of 3D effects and viscosity.

3. Applications The motivation for developing this framework was to improve the predictive capability of space and aeronautical applications for NASA and the aerospace industry, specifically focusing on Launch, Ascent and Vehicle Aerodynamics (LAVA). The following sections highlight some of the available capabilities for pressure, thermal and acoustic analysis for the launch environment, steady and unsteady ascent aerodynamic analyses of launch vehicles, and vehicle aerodynamics applications. For all applications available experimental data is used for comparison to validate the presented simulation capabilities.

31

Experiment and simulation by Giordano et al.

LAVA-Cartesian simulation

t = 0s

t = 2.8 × 10−4 s

t = 4.2 × 10−4 s

t = 5.6 × 10−4 s

t = 7.0 × 10−4 s

t = 8.4 × 10−4 s Figure 13: Comparison of instantaneous Schlieren contours for L = 50mm and t = 1mm.

(a)

(b)

(c)

Figure 14: (a) Particle visualization of plume interaction with the geometrically complex launch environment colored by Mach number (high values are colored white). (b) IOP/DOP waves traveling from the flame trench up toward the SLS heavy lift vehicle. (c) Plume-tower interaction during SLS vehicle lift-off.

32

3.1. Launch Environment With the retirement of the Space Shuttle program, NASA is developing a new heavy-lift capability, the Space Launch System (SLS). In addition, commercial vehicles such as SpaceX’s Falcon, ATK’s Liberty and ULA’s Delta IV Heavy and Atlas V rockets, could be potential candidates to be launched from the Kennedy Space Center (KSC). Flow physics with widely varying scales are present during the launch of a space vehicle. The initial stage of the rocket launch contains highly complex fluid dynamics that need to be carefully analyzed to avoid damage to the launch site or even catastrophic failure of the launch vehicle. While the rocket engines are firing, the launch pad and the flame trench are subjected to extremely harsh flow conditions with highly unsteady pressure loads and high heat transfer rates (see Figure 14a). During the ignition of the Solid Rocket Booster (SRB), Ignition Over-Pressure (IOP) and Duct Overpressure (DOP) waves are generated traveling up the launch vehicle (see Figure 14b), which may cause damage to the vehicle, avionics or the payload. At the same time, the interaction of the rocket plume with the launch platform and tower generates intense pressure waves and high structural loads. All these scenarios share common features of unsteady and highly shocked flows. The LAVA solver was used extensively for the re-design efforts of the KSC launch pad. 3.1.1. Pressure Environment The Cartesian AMR solver was extensively utilized for characterizing the pressure environment during rocket launch. The focus of the effort was the assessment of unsteady pressure loads on critical components of the launch site and vehicle, such as side walls of the flame trench, main flame deflector, umbilical hardware of the launch tower, and payload compartments on the vehicle. The Cartesian immersed-boundary methodology with Detached Eddy Simulation (DES) and slip walls was applied to improve estimates of peak pressures (subsequently low-pass filtered for structural analysis). The generation of water-tight surface triangulations took several days to create, which is typical for a CFD approach. An advantage of the Cartesian method is the automated volume grid generation given only a water-tight surface triangulation. This eliminates the time-consuming grid generation procedure of structured curvilinear or unstructured grids. The grid has approximately 300 million cells and was run on 1000 computational cores for less than a week, where the duration of the simulation was sufficient to obtain statistics on tower impingement pressures. The Cartesian AMR solver also provide the capability of full moving body lift-off 6DOF simulations.

33

For these studies, time-dependent simulations were conducted with a time-step of Δt = 3.5 × 10−5 sec employing a second-order accurate dual-time-stepping approach. At least two orders of magnitude convergence in the right-hand-side residual, measured by the max-norm was ensured, before advancing to the next time-step. The second-order MUSCL scheme with a minmod flux limiter was used. The fully coupled multi-species formulation was employed considering three gas species: air, core-stage engine exhaust gas, and SRB exhaust gas. The chemical composition and the gas properties at the nozzle exit were determined by 1D chemical equilibrium analysis. It was assumed that once the gas exits the nozzle the composition of the plume does not change significantly. 80

200

Raw Measurement Filtered Measurement LAVA

70

150 P ressure (P SIA)

P ressure (P SIA)

60 50 40 30 20

100

50

10 0 0.0

0.5

1.0 T ime (s)

1.5

0 0.0

2.0

(a) Pressure at Bottom Sensor

0.5

1.0 T ime (s)

1.5

(b) Pressure at Top Sensor

Figure 15: LAVA surface pressure predictions vs. flight measurements at the two sensor locations.

Verification of this approach was performed using the MFD data collected during the Space Shuttle launch STS-135[66]. The simulation modeled the rocket plume development starting from ignition time. The time-dependent conditions at the rocket plenum were taken from STS-1 flight measurements (which differ from the STS-135). A non-reacting, single-species plume exhaust gas was assumed and both the exhaust gas and free air were modeled as ideal gas. The water sound suppression system was not modeled because its main purpose is to reduce acoustic loads on the vehicle. No significant effect on the peak surface heating at least at the primary and secondary impingement locations (the most critical regions) was expected because the water is not able to penetrate inside the jet. The water vapor environment mainly affects the phase speeds of pressure waves. Note that the refractory material coating on the MFD is highly irregular (containing roughness, valleys and peaks). The current simulation, in contrast, used a smooth MFD surface while further ignoring the surface reactions and recession during launch. This will be important below when comparing the 34

2.0

heat transfer predictions to the STS flight data. In this simulation, the vehicle was held stationary, ignoring its ascent, which is a valid assumption for the first second of launch where the vehicle moves only slightly. Figures 15a-b show a comparison of the predicted versus measured pressure data at the two sensor locations. Peak pressure (due to IOP) and the subsequent sustained pressure levels are closely predicted. However, a phase shift between the predictions and the flight data is present. Omission of the water sound suppression system (and hence the change in the speed of sound) is likely a major contributor to this discrepancy. It should also be emphasized, once again, that the unsteady rocket plenum data was borrowed from the STS-1 flight data, and the discrepancies with that of STS-135 data are not known to the authors at the time of this writing. 3.1.2. Thermal Environment During launch,

the MFD

must withstand the harsh environment of direct rocket plume impingement while deflecting the plume away from the launch vehicle.

Besides structural loads

due to impingement, large heat transfer rates on the MFD surface are also a major concern. The current Kennedy Space Center’s (KSC) MFD includes a coating with a refractory material that absorbs heat with significant erosion occuring over time.

Figure 16: Unstructured computational domain for STS-1 validation case.

While this approach has been successful, it requires frequent, high cost repairs after each launch. KSC is currently investigating refinements and/or alternative approaches. During the evolution of this new design, accurate CFD predictive capabilities are essential to perform design trade studies and to investigate MFD compatibility with several different launch vehicles (SLS and commercial). Traditional CFD wall boundary conditions (e.g., adiabatic and isothermal) fall short for the

35

purposes of MFD thermal environment characterization. The adiabatic approach may offer a conservative estimate of the surface temperature, while the isothermal approach (depending on the assumed wall temperature) may provide conservative heat transfer rates. However, both approaches prove to be inadequate, as they result in overly conservative and ambiguous computational predictions. The conjugate heat transfer capability of LAVA, as mentioned earlier (see Figure 10), is well suited for the task of predicting time-dependent surface material response to plume heating. A body-fitted unstructured mesh was generated to model the complex geometry including highresolution boundary layer prismatic cells. The computational mesh, along with two sensor locations on the MFD, is shown in Figure 16. The mesh consisted of approximately 21 million arbitrary polyhedral cells and features polygonal prisms at the boundary layers, with a resolution that meets the required wall spacing of y + < 1. Note that for heat transfer calculations the definition of the y + value may not be appropriate especially for high enthalpy flows. The sensors were embedded in stainless steel encasings and recorded presssure, temperature and heat transfer rates during the STS-135 launch. Figure 17 shows the instantaneous surface temperature and heat flux distributions on the MFD surface at two different time instances. The surface temperature is capped at the melting point of the material. At t = 0.63 seconds, most of the MFD surface refractory material has reached its melting temperature. Note that the refractory material has a significantly lower melting temperature than the stainless steel that was used to enclose the sensors. In addition, the higher heat capacity of the stainless steel corresponds to a slower temperature rise. As a result, the sensor location points are relatively colder (barely visible in this view). Figure 17 also reveals that the surface heat flux is highly localized to the primary and secondary impingement locations of the SRB plumes. For the computational results, the pressure peaks due to IOP (see Figure 15a-b) have a corresponding peak in the temperature and the heat transfer rates (see Figure 18a-d). The corresponding peaks are not as clearly visible in the measurements. Furthermore, the time-shift between the computed and measured results seems more severe than in the pressure comparison. These two issues led to the suspicion that the temperature and heat flux sensors’ delayed response as probable causes. In spite of these discrepancies, both the unsteady temperature and heat flux data predictions are favorable. Temperature at the bottom sensor location closely follows the measured time history. Considering the previously mentioned time shift, the computed temperature at the top sensor also reproduces the initial temperature rise rate as well as the subsequent degrading of the slope. However,

36

Figure 17: Instantaneous surface temperature (top row) and heat transfer rate (bottom row) contours at two different time instances.

the measured temperature at the top location exhibits a somewhat large spike at around t = 1.0 second, followed by a slight cooling and plateauing before recovering to the expected increase for t > 1 second. Heat flux measurements at both locations are also well within the experimentally measured range. Most importantly for the design of the new MFD, the maximum amplitudes are well captured. Note that a great number of simplifying assumptions, as mentioned earlier, had to be made in order to provide an engineering estimate of the MFD heating while preserving turnaround times that are conducive to integration with the actual design process. In perspective of these, the conjugate fluid flow/surface heat transfer capability of the LAVA solver proved to be an extremely useful design evaluation tool that KSC engineers relied on regularly. 3.1.3. Launch Pad Acoustics During liftoff, a severe acoustic environment is generated by the exhaust jet plumes and their impingement on the flame deflector and mobile launcher. An accurate, robust, and efficient procedure

37

Raw Measurement LAVA

2000

1500

T emperature (F )

T emperature (F )

2000

1000

500

0 0.0

1500

1000

500

0.5

1.0 T ime (s)

1.5

0 0.0

2.0

(a) Temperature at Bottom Sensor

1.0 T ime (s)

1.5

2.0

(b) Temperature at Top Sensor

700

2000

Raw Measurement LAVA Heat F lux (BT U/(f t2 · s))

600 Heat F lux (BT U/(f t2 · s))

0.5

500 400 300 200

1500

1000

500

100 0 0.0

0.5

1.0 T ime (s)

1.5

0 0.0

2.0

(c) Heat Rate at Bottom Sensor

0.5

1.0 T ime (s)

1.5

(d) Heat Rate at Top Sensor

Figure 18: LAVA temperature and heat transfer predictions vs. flight measurements at the two sensor locations.

has been developed using the LAVA framework to predict the jet acoustic environment. Validation of the procedure is performed for a perfectly expanded, Mach 1.8, cold jet impinging on a 45 degree flat plate. Experimental data [69] was made available through collaboration between NASA and the Japan Aerospace Exploration Agency (JAXA). A diagram of the problem setup along with the sensor locations is reproduced from the Nakanishi et al. [69] in Figure 19a, along with a slice of the overset grid system in Figure 19b. A brief description of the procedure is provided below: First, the nozzle interior flow was analyzed using steady Reynolds-averaged Navier-Stokes (RANS) simulations to obtain nozzle exit conditions assuming fully turbulent flow throughout the nozzle. An additional steady RANS computation is performed on the external jet impinging onto the plate. Results from this simulation are used in the next step, and also provide good information on necessary improvements to the computational 38

2.0

grid, such as wall and shear layer spacing. Second, acoustic source identification is performed following the procedure outlined in Tosh et al. [101], but modified to use steady RANS information, which is relatively inexpensive to compute.

(a)

(b)

(c)

(d)

Figure 19: (a) Diagram of experimental setup including far-field acoustic sensors at 40 nozzle diameters from jet impingement location, courtesy of Nakanishi [69]. (b) Slice of the curvilinear overset grid system cutting through the nozzle and flat plate including additional point sensors added to the computational simulation. (c) Contour plot of the scalar source indicator function illustrating where the major noise generating regions are located in the domain. (d) Acoustic surface that is embedded in the CFD domain and used to collect time-accurate data from the CFD simulation for acoustic propagation.

39

A Reynolds velocity decomposition is performed on the Lighthill stress tensor, ui = Ui + ui :

Tij

=

  ρui uj + γp − ρc20 δij /γ

=

ρUi Uj

mean velocity contribution

+

ρUi uj + ρUj ui

linear velocity fluctuation

+

ρui uj quadratic velocity fluctuation   γp − ρc20 δij /γ entropy violation

+

(36)

(37)

The mean velocity contribution is silent, while the linear velocity fluctuation are properly modeled by the linear acoustic solver. The quadratic velocity fluctuations must be captured using the full NavierStokes equations modeled using CFD. The Boussinesq approximation is now used to approximate the time-average of the quadratic velocity fluctuations:  T Rij = ρui uj ≈ −μT

∂ui ∂uj + ∂xj ∂xi

 (38)

This equation represents a tensor of quantities known from the steady RANS analysis performed in the first step. In order to identify the noise source generating region, a scalar source indicator is computed from the trace of Equation (38):

T Rdiag =

" 2 + T R2 + T R2 T R11 22 33

(39)

In Figure 19, the scalar source indicator shows that the major noise generating regions are confined to the jet shear layer, plate shock at impingement location, and the wall jet created from the deflected plume. Once the acoustic sources have been identified, acoustic surfaces are embedded into the CFD domain, which encloses the noise source regions as shown in Figure 19d. In the third step, time-accurate large eddy simulation (or detached eddy simulation) of the external jet flow is performed to predict the acoustic noise. The noise is generated by unsteady shear layers, oscillating plate shock, and the wall jet created by the deflected plume. A time history of the CFD solution is mapped to the embedded acoustic surface at each time-step for acoustic propagation. Finally, once a sufficient amount of solution time-history is gathered on the acoustic surface, a

40

Fast Fourier Transform (FFT) algorithm is applied at each surface location to transform the data into the frequency domain. Once in the frequency domain, a linear Helmholtz BVP is solved including both the radiating surface data provided by the CFD and scattering boundaries representing the exterior nozzle, deflector, and any other relevant geometry. The solution procedure of the Helmholtz BVP was described in Section 2.6.2. The solution of the acoustic wave propagation can then be evaluated at any location outside the radiating/scattering surfaces. The acoustic solution can also be transformed back into the time domain using the inverse FFT for visual post-processing.

(a)

(b)

Figure 20: Jet impingement at 45 degrees: (a) iso-contours of vorticity magnitude and dilatation contours (red-whiteblue) and (b) passive trace particle visualization colored by vorticity magnitude.

The four-step procedure, described above, was applied to the Mach 1.8 perfectly expanding jet impinging on a flat plate at 45 degrees. The steady RANS analysis was performed using body-fitted structured curvilinear overset grids in order to efficiently model the wall boundary layers. Two grid approaches have been applied in the third step, the Cartesian grid with immersed-boundary conditions and the body-fitted structured curvilinear overset grid approaches. For the present simulation, it turned out that the dominant noise sources could be well predicted with both approaches. A comparison of the velocity profiles in the wall jet showed close agreement between the immersed boundary and curvilinear simulations (see Brehm et al.[25]). Figure 20a displays an instantaneous snapshot of the flow field showing iso-contours of vorticity magnitude and dilation contours on a cutting plane through the center of the jet. Figure 20b displays the spectral energy of the u = u − U

41

disturbance velocity, Eu , at four different locations in the shear layer region of the jet and compares the spectra to the expected asymptotic behaviors of a turbulent flow in the inertial sub-range and dissipation range of the von K´arm´an energy spectrum (see Pope[78]). The point probes are located at four vertical positions, i.e., z/D = 1.67(A), 3.33(B), 3.89(C), and 4.44(D), in the shear layer region just below the lip of the nozzle exit. The spectrum of the disturbance velocity at point location A, achieves the asymptotic behavior of f −5/3 in the inertial sub-range of the von K´arm´an energy spectrum. The f −7 asymptotic behavior for the dissipation range is shown as a blue-dashed line. In the dissipation range, most of the turbulent kinetic energy is dissipated by the physical action of the molecular viscosity. For the current Reynolds number of 1.6 × 106 , it is expected that (even with more than 0.5 billion grid points) most of the dissipation is actually introduced by the numerical scheme and not by the physical action of the molecular viscosity. The energy spectra at point A suggests that the unsteadiness in the flow is relatively well resolved until approximately St = 3 to 150

Exp. CFD/CAA

PSD (dB/St)

130 120 110

60

80

100

120

Angle to inlet axis (deg.)

100 -2 10

140

(a)

10-1

St

100

(b)

150

150

Exp CFD/CAA

Exp CFD/CAA

PSD (dB/St)

140

PSD (dB/St)

140 130

130

120

120

110 100 -2 10

Exp CFD/CAA

140

OASPL (dB)

166 164 162 160 158 156 154 152 150 148

110 10-1

St

100 -2 10

100

(c)

10-1

St

100

(d)

Figure 21: (a) Comparison of the predicted OASPL and experiment at 40 nozzle diameters from impingement location, where dashed line represents ±3 dB increment on experimental data [69]. Comparison of the Power Spectral Density (PSD) spectrum at 75 (b), 90 (c), and 105 (d) degrees. See Figure 19a for sensor locations.

42

4. A complete description of the shear layer breakdown and jet impingement physics can be found in Brehm et al. [26, 25]. Figure 21a plots the predicted Overall Sound Pressure Level (OASPL) against experimental data at 40 nozzle diameters from the impingement location. The CFD/CAA prediction is within ±3 dB of the experiment data from 70 to 135 degrees, but overpredicts the OASPL for less than 70 degrees. This is likely caused by insufficient mesh resolution in the wall jet or lack of modeling for the end plate effects. The Power Spectral Density (PSD) spectra at 75, 90 and 105 degrees are shown in Figure 21b-d. Very good agreement is observed in the individual spectra indicating that the acoustic propagation code performs well in propagating the waves associated with each individual frequency to the far field. A more complete acoustic analysis including sensitivity to acoustic surface position, boundary conditions, propagation type, and CFD modeling is included in Housman et al. [49]. 3.2. Heavy Lift Vehicle Ascent Aerodynamics 3.2.1. Steady Analysis for Ascent Trajectory The design analysis cycle of NASA’s future heavy-lift vehicles, including trajectory optimization and structural analysis, relies heavily on predictions of the aerodynamic loads throughout the ascent trajectory. A combination of wind tunnel testing and CFD analysis are utilized during each design iteration. Aerodynamic databases are generated that encompass a wide range of Mach numbers from subFigure 22: Arbitrary polyhedral unstructured mesh used for LAVA simulations.

Lower-left:

Multi-

sonic to high supersonic, as well as pitch and yaw angles in the range of 10 degrees. CFD analysis is

purpose crew vehicle launch abort system close-up. Upper-right: Booster skirt close-up.

used to provide integrated and distributed loads for the flight scale vehicle, requiring resolution of geo-

metric protuberances with a wide range of length scales. Viscous analysis are commonly performed using NASA’s well-established OVERFLOW code employing structured overset grids. The grid generation process for unstructured meshes is significantly less time consuming that for structured overset grids. The structured overset grid approach, however, is a proven technology for providing reliable loads predictions.

43

To gain confidence in LAVA’s unstructured solver module, a thorough code-to-code comparison was performed for Design Analysis Cycle 2 (DAC2) of the SLS vehicle with an arbitrary polyhedral unstructured mesh. Results for a limited set of the ascent trajectory points are presented here. The unstructured mesh contained polygonal prism layers near the wall boundaries with a resolution to ensure that y + < 1, at the first cell center off the wall, was met. The mesh generation process was completed in two days by a single individual, in comparison to the structured overset grid approximately three weeks of effort by three people. The unstructured mesh contains approximately 40 million cells (see Figure 22), and the OVERFLOW grid system has about 145 million grid points. The Spalart-Allmaras turbulence model was used for both the LAVA and OVERFLOW analyses. Figure 23 compares the integrated force predictions from the two codes. LAVA consistently predicts the vehicle aerodynamic loads within 5% of the OVERFLOW values, with a few exceptions where the actual values are relatively low (hence amplifying the difference percentage). The trends of the force curves with respect to the Mach number are also qualitatively the same. Actual values for the axial and normal force coefficients have been removed due to public domain restrictions. Figure 24 offers a more detailed comparison via surface pressure coefficient distributions from OVERFLOW and LAVA simulations. The bulk flow features are similarly captured by both codes, with only slight differences observed. Results from this code-to-code verification study provide confidence in utilizing the unstructured mesh capabilities of LAVA in ascent aerodynamic database generation for the next design cycle.

LAVA-Uns OVERFLOW

Δ 0.75% Normal Force Coefficient

Axial Force Coefficient

Δ 2.22%

Δ 0.98%

Δ 0.92%

Δ 0.16% Δ 5.71%

Δ 2.18% Δ 1.74%

Δ 5.16%

Δ 4.50% 1

2

3 Mach

4

5

1

2

3 Mach

4

5

Figure 23: Integrated aerodynamic force coefficient in the axial (left) and normal (right) directions for a range of Mach numbers at an angle of attack of α = 6 degrees with percent difference between LAVA and OVERFLOW annotations.

44

(a) LAVA

(b) OVERFLOW Figure 24: Comparison of surface pressure contours for the SLS during supersonic ascent using (a) LAVA and (b) OVERFLOW.

(a)

(b)

Figure 25: (a) Block-structured Cartesian grid used for simulations of unsteady vortex shedding for a heavy lift SLS vehicle for NASA Langley TDT wind tunnel model. (b) Instantaneous snapshot of coefficient of pressure on surface.

3.2.2. Unsteady Aerodynamic Analysis Protuberances and attachment hardware on launch vehicles may cause significant aerodynamic unsteadiness, which may lead to cyclical loads and an undesirable acoustic environment. Timedependent viscous and inviscid predictive capabilities are crucial in order to understand the unsteady flow physics. Inviscid analysis using Cartesian and viscous analysis using unstructured topologies of LAVA have been performed to assist during the advanced design cycles. Rapid design analysis can

45

10−1

106 Wind Tunnel (Raw) Wind Tunnel (Filtered) Cartesian Inviscid (Filtered)

105

PSD ((Pa)2/Hz)

|Y(freq)|

104

10−2

103 102 101 100 10

10−3 1 10

102 Freq (Hz)

−1

10−2 10

103

(a)

Wind Tunnel (Raw) Wind Tunnel (Filtered) Cartesian Inviscid (Filtered) Unstructured Viscous DES (Filtered)

100

Freq (Hz)

1000

10000

(b)

Figure 26: (a) Spectral amplitudes (|Y (f req)|) of pressure data sampled just downstream of SRB attachment point, where red is LAVA inviscid solution using block-structured Cartesian AMR and black is wind tunnel data. Note, frequency binning was used to filter the PSD data (11Hz bin size) for numerical and experimental results [17]. (b) SLS buffet results for NASA Ames UPWT model. Spectral PSD of inviscid Cartesian, unstructured viscous DES and experimental pressure data [81] sampled downstream of attachment hardware. Note: Frequency binning was used to filter the PSD data (60Hz bin size) for numerical and experimental results.

be done with the Cartesian mesh, which eliminates time-consuming volume grid generation if the inviscid flow assumption is a valid approximation. The NASA Langley Transonic Dynamics Tunnel (TDT) test data [17] was utilized to validate the inviscid approach. All data are at wind tunnel scale (i.e., frequencies are much higher than flight scale). See Figure 25 for the computational grid and for the instantaneous surface pressure. Sensitivity studies of time-step, sub-iteration convergence, mesh resolution, and duration of the simulation (not shown here) were conducted in order to establish a best practice and confidence in the simulations. Comparison of the numerical simulations and experimental data focused on the oscillatory wake region behind the SRB forward attachment hardware. The spectral peak at around 400Hz in the PSD plot in Figure 26a is associated with the wake dynamics behind the SRB forward attachment. Comparison of simulated pressure time series on the test article and the TDT wind-tunnel data match well. Subsequent simulations and design iterations focused on the NASA Ames Unitary Plan Wind Tunnel (UPWT) test article. Initially, an unsteady inviscid adaptive Cartesian IB approach was applied for rapid design iterations to support the wind tunnel test. Subsequently, an unsteady viscous body-fitted unstructured mesh solution were used to explore the viscous boundary layer effects. A representative unstructured polyhedral grid for this analysis is shown in Figure 22. Spectral PSD plots comparing inviscid and viscous DES results to experimental data [81] are shown in Figure

46

26b. The numerical solutions closely predict the frequency location of the maximum PSD for several probe locations. 3.2.3. Plume-Induced Flow Separation Analysis Rockets at high-altitude are subject to a fluid dynamics phenomenon known as Plume-Induced Flow Separation (PIFS) [98]. One cause of a strong adverse pressure gradient during ascent is the expansion of the exhaust plume as the rocket gains altitude. In a low ambient pressure environment, the high pressure at the nozzle exit rapidly expands the exhaust jet in both downstream and radial directions. This produces an obstruction to the free-stream flow, which forms an adverse pressure gradient near the aft section of the rocket. Ultimately, the flow separates and recirculation from the base of the vehicle to the upstream separation point allows convective transport of hot exhaust gas along the surface of the vehicle. The distance between the end of the vehicle and the separation point of the surface is denoted as the PIFS distance. Figure 3.2.3 shows a frame from the Apollo 6 flight with PIFS visible as well as a schematic of the PIFS distance definition. Accurate prediction of the PIFS distance is critical to the design of the thermal protection system. The Saturn V rocket flight data is used as a validation case for the Cartesian-unstructured overset coupling capabilities within the LAVA solver. Saturn V is an Apollo-era launch vehicle, with a vertical-stack design that resembles design concepts for future heavy lift launch vehicles currently being developed at NASA. All flight trajectory information is derived from the flight evaluation report from the launch of Apollo 11 [83]. The free-stream conditions for each trajectory point are

(a) Photo from Apollo 6 flight

(b) Schematic

Figure 27: PIFS examples: (a) frame from chase plane footage of Apollo 6 flight with PIFS visible by extent of radiating exhaust gas, and (b) diagram of Saturn V station zero (located 2.84 m downstream from the core base) and measurement of PIFS distance.

47

M∞ 1.5 2.7 4.4 6.5

P∞ (Pa) 12111 2250 151 22

T∞ (o K) 217 221 264 247

ReD 6.1522x107 2.2623x107 1.6970x106 4.0600x105

Table 1: Free-stream conditions for the Saturn V PIFS simulations.

(a) Unstructured grid topology

(b) Hybrid grid topology

(c) Close-up of core

(d) Close-up of nozzle

Figure 28: Overview of (a) unstructured and (b) hybrid computational grid topology for PIFS application as well as close-ups of (c) core and (d) nozzle for the hybrid grid.

48

summarized in Table 3.2.3. PIFS distance measurements were made from video footage of the launch, and contain an uncertainty of approximately 10%. At Mach 4.4, the PIFS distance was measured at 15 meters. To reduce loads on the vehicle and crew, the center engine cut-off (CECO) event occurs at 136 seconds after launch, leading to a brief reduction in PIFS distance before it climbs back to approximately 33 meters at Mach 6.5. Gusman et al. [41] examined the numerical prediction of PIFS using OVERFLOW with structured curvilinear grids: a similar study was conducted by Deere et al. [33] using unstructured grids and USM3D. To establish best practices for PIFS analysis for both approaches, sensitivity to grid resolution (wall spacing and grid resolution), boundary conditions, fidelity of numerical discretization (single species vs. multi-species) and turbulence models were analyzed. The published results indicate that the PIFS distance was sensitive to both the turbulence model and the boundary conditions. The best practices for USM3D indicated that the k −  turbulence model and SA turbulence model had the best correlation to flight data. Using OVERFLOW, the computed PIFS distances were in good agreement with flight data using SST and underpredicted using SA. However, at the time of the study the standard SA turbulence model was not used. The goal of the current work is to examine the PIFS predictive capabilities using LAVA with the standard SA and SST turbulence models. The fully unstructured polyhedral grid capability of the LAVA solver was applied to the PIFS analysis. With the unstructured grid approach, high-fidelity body-fitted grids were generated with high-resolution surface and volume meshes in the predicted PIFS regions. High-aspect ratio prismatic cells, grown from the surface, are used to capture the boundary layer, and then transition to arbitrary polyhedra away from the body. To resolve the PIFS structure, volumetric controls are utilized to limit the cell size near the surface. The unstructured grid topology around the Saturn V rocket is shown in Figure 28a. The typical unstructured grid count for these simulations was 13 million mixed element cells. In order to resolve the plume of the jets with an unstructured grid approach, a large number of unstructured cells are required. To assess the sensitivity of the plume resolution to the predicted PIFS distance, the hybrid overset grid approach was applied, which couples the near-body unstructured grid with structured Cartesian off-body grids. The advantage of using a hybrid overset grid technique for this application is the ability to efficiently obtain viscous wall resolution and resolve flow structures. An unstructured grid with similar plume resolution can exceed 40 million mixed

49

(a) M∞ = 1.5

(b) M∞ = 2.7

(c) M∞ = 4.4

(d) M∞ = 6.5

Figure 29: Mach number distribution for four points in the trajectory using SST with the hybrid grid: (a) M∞ = 1.5, (b) M∞ = 2.7, (c) M∞ = 4.4 and (d) M∞ = 6.5.

element cells, which is computationally much more expensive than utilizing the Cartesian solver. Similarly, using a stand-alone Cartesian method to resolve the boundary layer would lead to an excessive number of cells. With the hybrid overset grid approach, the unstructured domain is limited to 2 meters from the Saturn V surface with efficient off-body Cartesian cells for plume resolution. An overview of the hybrid overset grid system and close-up of the Saturn V tip and nozzle are shown in Figure 28. Prismatic unstructured near-body grids resolve the viscous boundary layer development on the vehicle surface and transition to isotropic Cartesian off-body grids containing 98 million cells. To reduce grid generation efforts, a computational grid using the smallest wall spacing was utilized for the full range of Mach numbers. The wall spacing was previously determined for the Mach 6.5 trajectory point considering both SA and SST turbulence models. Grid convergence is achieved with a wall spacing value of 9.48 × 10−4 meters or smaller, with a predicted PIFS value of approximately 29.7 meters for SA and 34.0 meters with SST compared to the 33.0 meters observed in flight. Steady-state simulations were completed for each trajectory point using the AUSMPW+ flux vector splitting scheme to discretize the convective fluxes and the GMRES linear solver. To improve robustness, the boundary conditions were ramped over 500 iterations with a maximum local CFL of 20. The axial stress tensor τy on the surface of the rocket was recorded and monitored to assess

50

solution convergence. The Mach number distribution is shown in Figure 29 for selected trajectory points using the SST turbulence model. The oblique shock angles, which reduce as the free-stream flow velocity increases, are visible at each protuberance on the surface of the rocket. For Mach numbers of 1.5 and 2.7, no visible PIFS separation distance were observed on the core, however, flow recirculation up to 3.0 meters is observed in the base region with all computational approaches. At Mach 4.4, the predicted PIFS distance increases to 14.5 meters using SA and 18.2 meters using SST. The SA model slightly underpredicts the PIFS distance while the SST model overpredicts the PIFS distance in comparison to the observed flight data of 15.0 meters. For Mach 6.5, after CECO, the predicted PIFS distance increases to 29.6 meters and 34.0 meters for the SA and SST turbulence models, respectively. Again, the SA model underpredicts while the SST model overpredicts the PIFS distance in comparison to the observed flight data value of 33.0 meters. The hybrid overset grid approach using SA yielded 14.3 and 30.3 meters for Mach 4.4 and 6.5, respectively, which is similar to the fully unstructured approach predictions. Figure 30a shows a bar graph summarizing the PIFS prediction using the unstructured and hybrid overset grid approaches for the four trajectory points, and Figure 30b overlays them onto the recorded flight data (note that the SA and SST results overlap for Mach 1.5 and 2.7). Good agreement is visible for both SA and SST with comparison to the 10% spread of the flight data. 3.3. Aeronautics Applications 3.3.1. Supersonic Low-Boom Transport

(a)

(b)

Figure 30: PIFS distances for SA and SST turbulence models for selected trajectory points, with comparison to flight data.

51

NASA’s Commercial Supersonic Transport Project is developing new technologies to enable supersonic civilian x

aircraft with low sonic boom signatures for overland flight. To assist in accomplishing this goal, improved computational prediction tools are being investigated. In order to

h

ck

Off-tra

Φ

ack

On-tr

establish the state-of-the-art in computational prediction Figure 31: 69-Degree Delta Wing Body test

tools, NASA, in collaboration with industry partners and

case geometry for the 1st AIAA Sonic Boom

the American Institute of Aeronautics and Astronautics

Prediction Workshop.

Schematic illustrating

(AIAA) has organized the 1st AIAA Sonic Boom Predic- the signal extraction coordinate system where tion Workshop. As part of the workshop, participants are

h is the radial distance from the centerline of the model and Φ is the circumferential distance

asked to perform two simulations and extract several near

from the on-track line in degrees.

and mid-field pressure signatures for comparison. 0.10 0.01

0.05

0.00

dp/pinf

dp/pinf

0.00 -0.05

-0.10

-0.01

-0.15

-0.20

-0.25

Unstructured AUSM Structured Roe Structured Central 0.0

5.0

X

-0.02

10.0

25.0

(a) h = 0.5 inches Φ = 0 degrees

Unstructured AUSM Structured Roe Structured Central 30.0

X

35.0

40.0

(b) h = 21.2 inches Φ = 0 degrees

Figure 32: Near-field pressure signatures dp/pinf for the 69-Degree Delta Wing Body configuration at h = 0.5 inches (a) and h = 21.2 inches (b) with Φ = 0 degrees.

One of the test cases that was studied is the Mach 1.7 flow around a 69-Degree Delta Wing Body test article [56] at zero degrees angle of attack with a Reynolds number of 1.38 × 107 per meter. Figure 31 displays the geometry of the 69-Degree Delta Wing Body test article used in the simulations with a schematic illustrating the coordinate system used for signal extraction. Nearfield pressure signatures for the 69-Degree Delta Wing Body configuration using the unstructured AUSMPW+, structured Roe, and structured central discretization methods are plotted for h = 0.5 inches, Φ = 0 degrees in Figure 32a and h = 21.2 inches, Φ = 0 degrees in Figure 32b. At the closer extraction location, grid effects are observed in the aft part of the signal associated with the

52

0.020

0.015

0.015

0.010

0.010

0.005

0.005

0.000

dp/pinf

dp/pinf

0.020

-0.005 -0.010 -0.015 -0.020

0.000 -0.005 -0.010

Experiment Unstructured AUSM Structured Roe Structured Central

-0.025 30.00

35.00

X

40.00

-0.015 -0.020 -0.025 30.00

45.00

(a) Φ = 0 degrees

Experiment Unstructured AUSM Structured Roe Structured Central 35.00

X

40.00

45.00

(b) Φ = 30 degrees

Figure 33: Mid-field pressure signatures dp/pinf for the 69-Degree Delta Wing Body configuration at h = 24.8 inches and (a) Φ = 0 deg, and (b) phi = 30 deg.

trailing edge of the wing, the base, and the initial part of the sting. In this part of the signature both the Roe and central methods, which use the structured grid, match very closely, while the AUSMPW+ scheme, which is used on the unstructured grid, has slightly lower peaks and small oscillations on the pressure recovery portion of the sting. These mesh differences quickly vanish at the further extraction location of h = 21.2 inches. At h = 24.8 inches and Φ = 0, and Φ = 30 degrees experimental data is available for code validation, see Cliff et al. [85]. Figure 33 displays the comparison of the experimental data with the unstructured AUSMPW+, structured Roe, and structured central methods. Good comparison between the locations of the shocks and slopes of the expansions are observed between the experimental and numerical signatures. The numerical methods predict larger peaks than experiment throughout the signal. This is most likely caused by the averaging procedure applied to the experimental data outlined in Cliff et al. [85]. The numerical predictions also capture the trailing edge wing shock at each of the Φ extractions, while the experiment shows a very smoothed shock on the off-track locations and no shock at Φ = 0. This again is caused by the averaging procedure. Overall the comparison is in good agreement, indicating that LAVA’s unstructured and structured curvilinear overset methodologies can be used as tools for the design of quiet supersonic commercial transports. In fact, the low cost of the structured overset approach along with the viscous turbulent capabilities indicate a strong advantage of this approach for mature designs.

53

3.3.2. Slat Noise High-fidelity time-accurate simulations were performed to investigate slat noise generated by a modified version of the 30P30N high-lift model with the flap retracted. Computations of the model included open-jet installation effects to better match the conditions of the experiment [4]. Preliminary analysis using unsteady RANS on both the full-span model and a two inch centerline cross-section were performed to assess the effect of the side walls on the flow-field near the center of the model. It was found in Housman et al.[47] that the side-walls have little effect on the centerline fluid dynamics for the configuration studied in this work. Since the experiment utilizes an acoustic array which can focus on the centerline slat portion of the model, the two-inch wide cross-section was used for the aeroacoustic analysis. The reference conditions used for the analysis are Mach 0.17 based on the exit velocity of the open-jet and a Reynolds number of 1.7 million based on the stowed chord of 16.73 inches. Standard sea level atmospheric conditions were used for the reference pressure and temperature. A timestep of 1 microsecond was used for the time-integration of the aeroacoustic analysis and 3 orders of magnitude residual reduction of the discrete non-linear equations was enforced at each physical time-step. This equated to approximately 4 to 5 sub-iterations per physical time-step. In order to start the simulations, an unsteady RANS analysis with a larger time-step and coarser grid was run for 30 flow throughs, the solution was then mapped to the finer aeroacoustic mesh. The geometric model used for full-span threedimensional analysis consisted of a planar nozzle representing the open-jet, side walls which are attached to the nozzle and hold the high-lift model in-place, a modified (slat) 30P30N model with the flap retracted, and the collector plate (not shown) which deflects the flow turned by the high-lift model towards a pressure exhaust vent to remove it from the anechoic chamber, see Figure 34.

Figure 34: Half-body isometric CAD view of

Modifications to the original 30P30N slat were made to

the simplified nozzle colored red, the side-wall colored blue, the installed high-lift model col-

the cove region where a straight section was added near

ored grey, and the collector plate colored green

the leading-edge and a thickened trailing-edge for installation purposes. A structured overlapping grid system was generated for the installed high-lift configuration consisting of 39 zones and 900,000 grid points for a single plane of the constant span 54

(a)

(b)

(c)

(d)

Figure 35: Improvements to the resolution capacity of the LAVA structured overset grid solver are illustrated by iso-contours of Q-criteria colored by stream-wise velocity. (a) Original discretization using fifth-order ZWENO interpolation to half-points and standard DDES length scale-definition, (b) modified length-scale definition allowing two-dimensional Kelvin-Helmholtz instability to be captured, (c) use of optimal weights in the interpolation procedure which remains stable for this sub-sonic flow condition, and (d) increase to seventh-order interpolation in span-wise direction and use of the blended central/upwind biased interpolation based on local Mach number.

model. A span-wise length of 2 inches was used based on the work from the BANC-III workshop [30], and a total of three span-wise grid resolutions were assessed with 33, 65, and 129 grid points in the span-wise direction for a total of 30, 59, and 116 million grid points in the volume grid. Simulating aeroacoustic phenomenon requires high-resolution schemes which can accurately capture linear and nonlinear wave propagation with 5 − 7 points-per-wavelength to be computationally feasible. Figure 35 illustrates a series of improvements made to the LAVA solver during the course of this work. In each of the sub-figures (a)-(d) an iso-contour of the Q-criteria is plotted and colored by stream-wise velocity (note the limits on the contour axis are clipped) for the medium grid with 65 points in the span-wise direction. In Figure 35 (a) the sixth-order WCNS is combined with fifth-order upwind biased ZWENO [15] interpolation along with the original DDES length scale definition. The iso-contour shows very little three-dimensional content and almost no sign of the two-dimensional Kelvin-Helmholtz instability which occurs near the leading-edge of the slat. Figure 35 (b) retains the numerical discretization but utilizes the new DDES length scale [92]. Clearly

55

-5 QFF CFD slatAng 30 AOA 27 CFD slatAng 30 AOA 26

-4

Cp

-3 -2 -1 0 1 -0.2

0

0.2

0.4

x/c

0.6

0.8

1

(a)

(b)

Figure 36: (a) Plot of time-averaged Cp assessing an angle of attack correction to account for possible chamber circulation effects or installation measurement errors. (b) Diagram showing Kulite dynamic pressure transducer locations on the slat (transducers 1-5) and the main element (transducers 6-10). Note that transducers 9 and 10 are located at the same location on the airfoil profile, but transducer 9 is on the center-span of the model while transducer 10 is positioned one inch off the centerline.

the new length scale definition enables the model to capture the Kelvin-Helmholtz instability and reduces the delay in transitioning to resolved three-dimensional structures. Since the flow is subsonic everywhere, the computationally expensive and diffusive ZWENO limiters can be replaced by the optimal weight upwind-biased interpolation stencils. Elimination of the extra artificial dissipation allows more three-dimensional content to be resolved, as shown in Figure 35 (c) figure. Higherorder interpolation and derivative operators can easily be applied in the span-wise periodic direction without requiring additional fringe points since no overlap exists in this direction. In addition, a blending of the upwind-biased and central interpolation operators can be applied which enables even more resolvable scales on the same mesh. This is demonstrated in Figure 35 (d) which utilizes seventh-order interpolation and eighth-order central difference operators in the spanwise direction, along with a blended central/upwind biased interpolation in all three coordinate directions. Figure 36 plots the Cp comparison between the experiment and the computational model at both 26 and 27 degrees. The two different AoAs were used to assess an AoA correction to account for possible chamber circulation effects or installation measurement errors. The Cp profiles between the experiment and the 26 degree computational results fall almost on top of each other. The adjusted angle of attack of 26 degrees was used for the subsequent aeroacoustic analysis. Figure 36 (b) a plot of the layout for the Kulite dynamic pressure transducers on the slat and main element where PSD spectral analysis is performed. A comparison of the predicted PSD spectrum at dynamic pressure transducers 1, 4, and 9 using each of the three span-wise mesh resolutions is shown in Figures 37 (a)-(c). At dynamic pressure transducer location 1 the PSD spectrum between the 65 and 129 span-

56

Frequency (Hz)

102

Frequency (Hz)

103

104

102

103

104

Kulite 1

Kulite 4

100

PSD (dB/Hz)

PSD (dB/Hz)

100

80

60

40

10-1

nspan 33 (navg=13) nspan 65 (navg=11) nspan 129 (navg=5) 100

80

60

nspan 33 (navg=13) nspan 65 (navg=11) nspan 129 (navg=5)

40

101

St

102

10-1

100

St

(a)

101

102

(b) Frequency (Hz)

102

103

104 Kulite 9

PSD (dB/Hz)

100

80

60

40

10-1

nspan 33 (navg=13) nspan 65 (navg=11) nspan 129 (navg=5) 100

St

101

102

(c) Figure 37: Plot of Power Spectral Density (PSD) spectrum using different span-wise resolution grids for (a) dynamic pressure transducer 1, (b) dynamic pressure transducer 4, and (c) dynamic pressure transducer 9.

wise grid point meshes are nearly on top of each other, while the coarser mesh shows a faster decay at higher frequencies. Dynamic pressure transducer 4 shows a different trend in which the coarse and medium grid spectrum are very close up to a frequency of 10 kHz, then the medium and fine grids become almost identical. Dynamic pressure transducers 1 and 3 are located where the wake from the slat leadingedge impinges on the lower surface of the slat. This is a source of the broadband noise which is clearly visible in the spectrum, and all tonal content at these locations are completely covered by the magnitude of the broadband noise. The remaining dynamic pressure transducers on the slat show several tonal peaks starting at approximately 6 kHz and show large magnitudes of broadband noise at lower frequencies. The PSD spectrum for the dynamic pressure transducer 9 is located on the main-element. Tonal content is observed for this dynamic pressure transducers starting at

57

6 kHz as observed on the slat. The tones are much more pronounced in the 20 − 40 kHz range on the main-element when compared to the slat. These tones are likely caused by shedding from the finite-thickness trailing-edge of the slat which is much larger relative to the slat chord in the model than what would typically be found on a full-size aircraft wing.

4. Summary The defining methods of the Launch Ascent and Vehicle Aerodynamics (LAVA) code have been introduced. A grid-flexible approach has been taken for the CFD solver which is agnostic to the actual grid type such that block-structured (body-fit or immersed) or body-fit unstructured grids can be utilized. An immersed-boundary approach has been described in which volume grid generation is fully-automated and adaptive mesh refinement can be used for enhanced solution accuracy. Some of the current challenges with IBM approach such as speeding up geometry queries, recomputation of irregular grid stencils, and wall treatment for high Reynolds number applications have been addressed and are still part of ongoing research. When near wall accuracy is very crucial for the flow physics the higher-order structured curvilinear solver with overset connectivity can be utilized. Various LEStype approaches are also available to simulate high Reynolds number, unsteady flow applications, such as those needed for computational aeroacoustic problems. In addition, an overset connectivity interface was introduced which allows coupling of unstructured near-body and structured Cartesian grids in an accurate and efficient manner. Various other solver modules that are available within LAVA have also been introduced. A linear acoustic solver has been developed that is coupled to the CFD solver by embedding the acoustic surface in the CFD domain and mapping the solution to the surface, independent of grid type. A conjugate heat transfer module modeling the unsteady heat conduction within the solid was described. Moreover, a brief description of a shell element based finite element solver that is coupled to the Cartesian solver was provided. The LAVA code has been applied to several NASA applications that are integral parts of major NASA missions. The immersed-boundary capabilities of LAVA were utilized in redesigning the main flame trench deflector at KSC through pressure environment predictions for numerous deflector shapes and many different launch vehicles. Thermal analysis on the deflector surface was also performed with the unstructured CFD solver coupled with conjugate heat transfer to aid in downselecting the various deflector shapes. Far-field jet acoustic prediction capabilities are validated for a supersonic jet impinging on a flat plate in preparation for prediction of the launch acoustic 58

environment. The higher-order shock capturing schemes available within LAVA were utilized for simulations of the unsteady flow field in the jet acoustic source region. Steady ascent capabilities have been compared to existing aerodynamic database generation techniques with success. This work led to a major contribution to the design of the SLS to characterize unsteady loads. Advantages of coupling the fast and efficient Cartesian capabilities with the viscous body-fitted unstructured capabilities is explored in predicting the Plume-Induced Flow Separation (PIFS) environment of the Saturn V during ascent. Comparisons with flight data indicate that both the Spalart-Allmaras and SST turbulence models are reasonably good at predicting the extent of the PIFS location. Finally, the body-fitted structured curvilinear and unstructured approaches are applied to vehicle aerodynamic applications related to NASA’s next generation aircraft. Accurate and economical CFD prediction procedures are demonstrated for both the low speed and supersonic flows. The higher-order accurate curvilinear solver with a modified DDES method was employed to simulate slat noise.

5. Ongoing Research In order to keep this paper at a reasonable length, the fluid structure interaction capabilities that are available within LAVA were only briefly discussed in this paper and future publications will provide additional information and validation results. Furthermore, code optimization techniques are currently being implemented in LAVA, and increased performance will expand the utilization of the code for more NASA applications. Coupling of the unstructured and structured curvilinear grid topologies through the existing overset connectivity interface is also being pursued in order to increase overall efficiency by generating structured grids on portions of the geometry which can easily be gridded, while using unstructured meshes on more difficult parts of the geometry. In addition, multi-phase and combustion chemistry models are being developed within LAVA to improve the modeling capabilities for launch site and ascent simulations. Multi-phase models are necessary to incorporate the effect of the water sound suppression system and combustion in the launch environment. Both of these effects will increase the fidelity of the heat-transfer predictions and provide more accurate exhaust plume dynamics, improving the prediction of ignition overpressure and jet acoustic phenomenon. Development of these physical models will be reported in future publications.

59

6. Acknowledgements Some of the application efforts presented in this paper were supported by NASA Kennedy Space Center (KSC), the NASA’s Space Launch System (SLS), and the Advanced Air Vehicles Program (AAVP). Computer time has been provided by the NASA Advanced Supercomputing (NAS) facility at NASA Ames Research Center. [1] Aftosmis, M. J., Berger, M. J., Saltzman, J. S., June 1998. Robust and efficient Cartesian mesh generation for component-base geometry. AIAA Journal 36 (6), 952–960. [2] Almgren, A. S., Bell, J. B., Colella, P., Howell, L. H., Welcome, M. L., 1998. A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations. J. Comp. Phys. 142, 1–46. [3] Anderson, W. K., Bonhaus, D., 1994. An Implicit Upwind Algorithm for Computing Turbulent Flows on Unstructured Grids. Computers and Fluids 23 (1), 1–21. [4] Bahr, C., Hutcheson, F., Thomas, R., Housman, J., May 2016. A Comparison of the Noise Characteristics of a Conventional Slat and Krueger Flap. In: 22nd AIAA/CEAS Aeroacoustics Conference, Lyon, France. Abstract submitted. [5] Balay, S., Brown, J., Buschelman, K., Gropp, W. D., Kaushik, D., Knepley, M. G., McInnes, L. C., Smith, B. F., Zhang, H., 2013. PETSc web page. Http://www.mcs.anl.gov/petsc. [6] Balay, S., Gropp, W. D., McInnes, L. C., Smith, B. F., 1997. Efficient management of parallelism in object oriented numerical software libraries. In: Arge, E., Bruaset, A. M., Langtangen, H. P. (Eds.), Modern Software Tools in Scientific Computing. Birkh¨auser Press, pp. 163–202. [7] Balsara, D. S., Shu, C.-W., 2000. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. Journal of Computational Physics 160 (2), 405–452. [8] Barad, M., Brehm, C., Kirs, C., May 2015. Parallel Adaptive High-Order CFD Simulations Characterizing Sofia Cavity Acoustics. In: 27th International Conference on Parallel CFD. [9] Barad, M. F., 2006. An adaptive Cartesian grid projection method for environmental flows. Ph.D. thesis, University of California, Davis. 60

[10] Barad, M. F., Colella, P., October 2005. A Fourth-Order Accurate Local Refinement Method for Poisson’s Equation. J. Comp. Phys. 209 (1), 1–18. [11] Barad, M. F., Colella, P., Schladow, S. G., 2009. An Adaptive Cut-Cell Method for Environmental Fluid Mechanics. Int. J. Numer. Meth. Fluids 60 (5), 473–514. [12] Barth, T. J., 1992. Aspects of unstructured grids and finite-volume solvers for the euler and navier-stokes equations. In: In AGARD, Special Course on Unstructured Grid Methods for Advection Dominated Flows 61 p (SEE N92-27671 18-34). Vol. 1. [13] Berger, M., Oliger, J., Mar. 1984. Adaptive mesh refinement for hyperbolic partial differential equations. J. Comput. Phys. 53, 484–512. [14] Berger, M. J., Colella, P., May 1989. Local adaptive mesh refinement for shock hydrodynamics. J. Comput. Phys. 82 (1), 64–84. [15] Borges, R., Carmona, M., Costa, B., Don, W., 2008. An improved weighted essentially nonoscillatory scheme for hyperbolic conservation laws. Journal of Computational Physics 227 (6), 3191–3211. [16] Borges, R., Carmona, M., Costa, B., Don, W. S., 2008. An improved weighted essentially nonoscillatory scheme for hyperbolic conservation laws. Journal of Computational Physics 227 (6), 3191–3211. [17] Brauckmann, G. J., Rausch, R. D., Blevins, J. A., 2013. Private communication. [18] Brehm, C., Barad, M., Housman, J., Kiris, C., 2014. A Comparison of Higher-Order Shock Capturing Schemes Within the LAVA CFD Solver. In: AIAA Conference, Jan 13-17, National Harbor, Maryland. pp. 1–38, aIAA–2014–1278. [19] Brehm, C., Barad, M., Housman, J., Kiris, C., 2015. A comparison of higher-order finitedifference shock capturing schemes. Computers & Fluids 122, 184 – 208. [20] Brehm, C., Barad, M., Kiris, C., Jun 13-17 2016. An immersed boundary method for solving the compressible navier-stokes equations with fluid-structure interaction. In: 54th AIAA Aviation and Aeronautics Forum, Washington, DC.

61

[21] Brehm, C., Barad, M. F., Kiris, C. C., 2016. Open rotor computational aeroacoustic analysis with an immersed boundary method. In: 54th AIAA Aerospace Sciences Meeting. p. 0815. [22] Brehm, C., Fasel, H., 2011. Immersed interface method for solving the incompressible navier stokes equations with moving boundaries. In: 49th AIAA Aerospace Sciences Meeting. [23] Brehm, C., Fasel, H., 2013. A novel concept for the design of immersed interface methods. J. Comput. Phys. 242 (0), 234 – 267. [24] Brehm, C., Hader, C., Fasel, H., 2015. A Locally Stabilized Immersed Boundary Method for the Compressible Navier-Stokes Equations. J. Comput. Phys. 295, 475 – 504. [25] Brehm, C., Housman, J., Kiris, C., 2016. Noise Generation Mechanisms for a Supersonic Jet Impinging on an Inclined Plate. Journal of Fluid MechanicsAccepted for publication. [26] Brehm, C., Sozer, E., Moini-Yekta, S., Housman, J., Barad, M., Kiris, C., Bruce, V., Parlier, C., 2013. Computational Prediction of Pressure Environment in the Flame Trench. In: AIAA Conference, June 24-27, San Diego, CA. pp. 1–32, aIAA–2013–2538. [27] Brehm, C., Sozer, E., Moini-Yekta, S., Housman, J. A., Barad, M. F., Kiris, C. C., Vu, B. T., Parlier, C. R., 2013. Computational Prediction of Pressure Environment in the Flame Trench. In: AIAA Paper, June 24-27, San Diego, CA. pp. 1–32, aIAA–2013–2538. [28] Chan, W., Gomez, R., Rogers, S., Buning, P., Jun 24-26 2002. Best Practices in Overset Grid Generation. In: 32nd AIAA Fluid Dynamics Conference and Exhibit, St. Louis, Missouri. AIAA–2002–3191. [29] Chan, W. M., June 2011. Developments in Strategies and Software Tools for Overset Structured Grid Generation and Connectivity. In: 20th AIAA Computational Fluid Dynamics Conference, Honolulu, Hawaii. [30] Choudhari, M., Lockard, D. P., June 2015. Assessment of Slat Noise Predictions for 30P30N High-Lift Configuration from BANC-III workshop. In: 21st AIAA/CEAS Aeroacoustics Conference, Dallas, Texas. AIAA-2015-3139. [31] Colella, P., Graves, D. T., Ligocki, T. J., Martin, D. F., Modiano, D., Serafini, D. B., Straalen, B. V., 2000. Chombo Software Package for AMR Applications - Design Document, unpublished. 62

[32] Deck, S., 2012. Recent improvements in the zonal detached-eddy simulation (ZDES) formulation. Theoretical and Computational Fluid Dynamics 26, 523–550. [33] Deere, K., Elmiligui, A., Abdol-Hamid, K., Jan 4–7 2011. USM3D Simulations of Saturn V Plume Induced Flow Separation. In: 49th AIAA Aerospace Sciences Meeting, Orlando, Florida. AIAA–2011–1055. [34] Deng, X., Mao, M., Tu, G., Liu, H., Zhang, H., 2011. Geometric conservation law and applications to high-order finite difference schemes with stationary grids. Journal of Computational Physics 230 (4), 1100–1115. [35] Fidkowski, K. J., Roe, P. L., 2010. An entropy adjoint approach to mesh refinement. SIAM Journal on Scientific Computing 32 (3), 1261–1287. [36] Frink, N. T., 1992. Upwind Scheme for Solving the Euler Equations on Unstructured Tetrahedral Meshes. American Institute of Aeronautics and Astronautics 30 (1), 70–77. [37] Frink, N. T., 1998. Tetrahedral Unstructured Navier-Stokes Method for Turbulent Flows. American Institute of Aeronautics and Astronautics 36 (11), 1975–1982. [38] Giordano, J., Jourdan, G., Burtschell, Y., Medale, M., Zeitoun, D., Houas, L., 2005. Shock wave impacts on deforming panel, an application of fluid-structure interaction. Shock Waves 14 (1-2), 103–110. [39] Gomez, R., Ma, E., June 20-22 1994. Validation of a Large Scale Chimera Grid System for the Space Shuttle Launch Vehicle. In: 12th AIAA Applied Aerodynamics Conference, Colorado Springs, CO. AIAA–1994–1859. [40] Gottlieb, S., Ketcheson, D. I., Shu, C.-W., 2009. High order strong stability preserving time discretizations. Journal of Scientific Computing 38 (3), 251–289. [41] Gusman, M., Housman, J., Kiris, C., Jan 4–7 2011. Best Practices for CFD Simulations of Launch Vehicle Ascent with Plumes - OVERFLOW Perspective. In: 49th AIAA Aerospace Sciences Meeting, Orlando, Florida. AIAA–2011–1054. [42] Hader, C., Brehm, C., Fasel, H., 2013. Numerical Investigation of transition delay in a Mach 6 Boundary Layer using porous walls. In: AIAA 2013-2740, 24-27 June, San Diego, CA.

63

[43] Henrick, A. K., Aslam, T. D., Powers, J. M., 2005. Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points. Journal of Computational Physics 207 (2), 542–567. [44] Housman, J., Barad, M., Kiris, C., June 2011. Space-time accuracy assessment of CFD simulations for the launch environment. In: 29th AIAA Applied Aerodynamics Conference. pp. 1–18, aIAA-2011-3650. [45] Housman, J., Kiris, C., Hafez, M., February 2009. Time-Derivative Preconditioning Methods for Multicomponent Flows - Part I: Riemann Problems. Journal of Applied Mechanics 76 (2). [46] Housman, J., Kiris, C., Hafez, M., March 2009. Time-Derivative Preconditioning Methods for Multicomponent Flows - Part II: Two-Dimensional Applications. Journal of Applied Mechanics 76 (3). [47] Housman, J., Kirs, C., May 2016. Slat Noise Predictions using Higher-Order Finite-Difference Methods on Overset Grids. In: 22nd AIAA/CEAS Aeroacoustics Conference, Lyon, France. Abstract submitted. [48] Housman, J. A., 2007. Time-derivative Preconditioning Method for Multicomponent Flow. Ph.D. thesis, University of California, Davis. [49] Housman, J. A., Brehm, C., Kiris, C., 2013. Towards jet acoustic prediction within the launch ascent and vehicle aerodynamics framework. The Journal of the Acoustical Society of America 134 (5), 4057–4057. URL http://scitation.aip.org/content/asa/journal/jasa/134/5/10.1121/1.4830807 [50] Hu, X., Adams, N., 2011. Scale separation for implicit large eddy simulation. Journal of Computational Physics 230 (19), 7240–7249. [51] Hu, X., Wang, Q., Adams, N., 2010. An adaptive central-upwind weighted essentially nonoscillatory scheme. Journal of Computational Physics 229 (23), 8952–8965. [52] Jiang, G.-S., Shu, C.-W., 1996. Efficient implementation of weighted eno schemes. Journal of Computational Physics 126 (1), 202–228. [53] Karypis, G., Schloegel, K., Kumar, V., 1998. Parmetis: Parallel graph partitioning and sparse matrix ordering library. 64

[54] Kim, K. H., Kim, C., Rho, O.-H., 2001. Methods for the accurate computations of hypersonic flows: I. ausmpw+ scheme. Journal of Computational Physics 174, 38–80. [55] Krist, S. L., Biedron, R. T., Rumsey, C. L., June 1998. CFL3D User’s Manual (Version 4.0). NASA Technical Memorandum TM-208444, NASA Langley Research Center. [56] L. W. Hunton and R. M. Hicks and J. P. Mendoza, 1972. Some Effects of Wing Planform on Sonic Boom. Technical Report NASA-TN-D-7160, NASA. [57] Ligocki, T. J., Schwartz, P. O., Percelay, J., Colella, P., 2008. Embedded boundary grid generation using the divergence theorem, implicit functions, and constructive solid geometry. In: Journal of Physics: Conference Series. Vol. 125. Institute of Physics Publishing, p. 012080. [58] Liu, X.-D., Osher, S., Chan, T., 1994. Weighted essentially non-oscillatory schemes. Journal of computational physics 115 (1), 200–212. [59] Luke, E. A., Tong, X.-L., Wu, J., Cinnella, P., 2004. CHEM 2: A Finite-Rate Viscous Chemistry Solver - The User Guide. Tech. Rep. MSSU-COE-ERC-04-07. [60] Lyrintzis, A., 2003. Surface integral methods in computational aeroacousticsfrom the (cfd) near-field to the (acoustic) far-field. International Journal of Aeroacoustics 2 (2), 95–128. [61] Martin, D. F., Colella, P., Graves, D., 2008. A Cell-Centered Adaptive Projection Method for the Incompressible Navier-Stokes Equations in Three Dimensions. J. Comput. Phys. 227, 1863–1886. [62] Mavriplis, D., 1992. Three-dimensional unstructured multigrid for the euler equations. AIAA journal 30 (7), 1753–1761. [63] Menter, F., July 1993. Zonal Two Equation k-ω Turbulence Models For Aerodynamic Flows. In: 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference, Orlando, FL. AIAA–93– 2906. [64] Merkle, C., Choi, Y., 1985. Computation of Low Speed Compressible Flows with TimeMarching Methods. International Journal for Numerical Methods in Engineering 25, 292–311.

65

[65] Mittal, R., Dong, H., Bozkurttas, M., Najjar, F., Vargas, A., von Loebbecke, A., 2008. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. Journal of Computational Physics 227 (10), 4825 – 4852. [66] Moini-Yekta, S., Barad, M., Sozer, E., Brehm, C., Housman, J., Kiris, C., 2012. Towards Hybrid Grid Simulations of the Launch Environment. In: ICCFD7-3102, July 9-13, Big Island, Hawaii. pp. 1–22. [67] Moini-Yekta, S., Barad, M., Sozer, E., Brehm, C., Housman, J., Kiris, C., 2013. Verification and Validation Studies for the LAVA CFD Solver. In: AIAA Paper, June 24-27, San Diego, CA. pp. 1–30, aIAA–2013–2448. [68] Nakahashi, K., June 2011. Immersed boundary method for compressible euler equations in the building-cube method. In: 29th AIAA Applied Aerodynamics Conference. AIAA-2011-3386. [69] Nakanishi, Y., Okamoto, K., Teramoto, S., Okunuki, T., Tsutsumi, S., March 2012. Acoustic characteristics of correctly-expanded supersonic hot impinging on an inclined flat plate. In: Asian Joint Conferencee on Propulsion and Power, AJCPP2012-129. [70] Nemec, M., Aftosmis, M., Jan 7-10 2008. Adjoint-based adaptive mesh refinement for complex geometries. In: 46th AIAA Aerospace Sciences Meeting, Reno, NV. [71] Nichols, R., Buning, P., August 2008. User’s Manual for OVERFLOW 2.1, version 2.1t. [72] Nichols, R. H., Tramel, R. W., Buning, P. G., 2008. Evaluation of two high-order weighted essentially nonoscillatory schemes. AIAA journal 46 (12), 3090–3102. [73] Nonomura, T., Fujii, K., 2012. Robust explicit formulation of weighted compact nonlinear scheme. Computers & Fluids. [74] Nonomura, T., Iizuka, N., Fujii, K., 2010. Freestream and vortex preservation properties of high-order weno and wcns on curvilinear grids. Computers & Fluids 39 (2), 197–214. [75] Nonomura, T., Morizawa, S., Terashima, H., Obayashi, S., Fujii, K., 2012. Numerical (error) issues on compressible multicomponent flows using a high-order differencing scheme: Weighted compact nonlinear scheme. Journal of Computational Physics 231 (8), 3181 – 3210.

66

[76] Ochmann, M., 1995. The Source Simulation Technique for Acoustic Radiation Problems. Acustica 81, 512–527. [77] Pandya, S., Chan, W., Kless, J., Jun 22–25 2009. Automation of Structured Overset Mesh Generation for Rocket Geometries. In: 19th AIAA Computational Fluid Dynamics Conference, San Antonio, Texas. AIAA–2009–3993. [78] Pope, S. B., 2000. Turbulent Flows. Cambridge U. Press. [79] Rogers, S. E., Kwak, D., Kiris, C., Apr. 1991. Steady and unsteady solutions of the incompressible Navier-Stokes equations. AIAA Journal 29, 603–610. [80] Rogers, S. E., Suhs, N. E., Dietz, W. E., 2003. Pegasus 5: an automated preprocessor for overset-grid computational fluid dynamics. AIAA Journal 41 (6), 1037–1045. [81] Ross, J. C., Panda, J., Blevins, J. A., 2013. Private communication. [82] Ruffin, S. M., Lee, J., 2009. Adaptation of a k-epsilon model to a cartesian grid based methodology. International Journal of Mathematical Models and Methods in Applied Sciences 3, 238–245. [83] Saturn Flight Evaluation Working Group, 1969. Saturn V Launch Vehicle Flight Evaluation Report-AS-506 Apollo 11 Mission. NASA Technical Memorandum TMX-62558, NASA/MSFC. [84] Schwamborn, D., Gerhold, T., Heinrich, R., 2006. The DLR TAU-code: Recent Applications in Research and Industry. In: European conference on computational fluid dynamics, ECCOMAS CFD. [85] S.E. Cliff and D.A. Durston and A.A. Elmiligui and J.C. Jensen and W.M. Chan, January 13-17 2014. Computational and Experimental Assessment of Models for the First AIAA Sonic Boom Prediction Workshop. In: 52nd AIAA Aerospace Sciences Meeting. Accepted. [86] Shu, C.-W., 2009. High order weighted essentially nonoscillatory schemes for convection dominated problems. SIAM review 51 (1), 82–126. [87] Shu, C.-W., Osher, S., 1988. Efficient implementation of essentially non-oscillatory shockcapturing schemes. Journal of Computational Physics 77 (2), 439–471.

67

[88] Shur, M., Spalart, P., Strelets, M., 2005. Noise Prediction for Increasingly Complex Jets. Part 2: Applications. International Journal of Aeroacoustics 4 (4), 247–266. [89] Shur, M., Spalart, P., Strelets, M., 2005. Noise Prediction for Increasingly Complex Jets. Part I: Methods and Tests. International Journal of Aeroacoustics 4 (3), 213–246. [90] Shur, M., Spalart, P., Strelets, M., Travin, A., 2003. Towards the prediction of noise from jet engines. Int. J. Heat Fluid Flow 24 (4), 551 – 561, selected Papers from the Fifth International Conference on Engineering Turbulence Modelling and Measurements. URL http://www.sciencedirect.com/science/article/pii/S0142727X03000493 [91] Shur, M., Spalart, P. R., Strelets, M., Travin, A., 2008. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. International Journal of Heat and Fluid Flow 29, 1638–1649. [92] Shur, M., Spalart, P. R., Strelets, M., Travin, A., 2015. An enhanced version of DES with rapid transition from RANS to LES in separated flows. Flow, Turbulence and Combustion submitted. [93] Sozer, E., Brehm, C., Kiris, C., 2014. Gradient Calculation Methods on Arbitrary Polyhedral Unstructured Meshes for Cell-Centered CFD Solvers. In: AIAA Paper 2014-1440, Jan 13-17, National Harbor, Maryland. pp. 1–24, aIAA–2014–1440. [94] Spalart, P., Deck, S., Shur, M., Squires, K., Strelets, M., Travin, A., 2006. A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theo. and Comp. Fluid Dynamics 20(3), 181–195. [95] Spalart, P., Jou, W., Strelets, M., Allmaras, S., 1997. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. In: In Advances in DNS/LES, ed. C Liu, Z Liu, Columbus. OH: Greyden Press. [96] Spalart, P. R., 2009. Detached-Eddy Simulation. Annual Review Fluid Mechanics 41, 181–202. [97] Spalart, S., Allmaras, S., January 1992. A One-Equation Turbulence Model for Aerodynamic Flows. In: 30th Aerospace Sciences Meeting and Exhibit, Reno, NV. AIAA–92–0439.

68

[98] Springer, A., January 10-13 1994. Experimental investigation of plume-induced flow separation on the national launch system 1 1/2-stage launch vehicle. In: 32nd AIAA Aerospace Sciences Meeting & Exhibit, Reno, NV. AIAA-1994-0030. [99] Strelets, M., January 2001. Detached eddy simulation of massively separated flows. In: 39th Aerospace Sciences Meeting and Exhibit, Reno, Nevada. AIAA-2001-0879. [100] Thompson, J. F., Thames, F. C., Mastin, C. W., 1974. Automatic numerical generation of body-fitted curvilinear coordinate system for field containing any number of arbitrary twodimensional bodies. Journal of Computational Physics 15 (3), 299–319. [101] Tosh, A., Liever, P., Owens, F., Liu, Y., June 04-06 2012. A High-Fidelity CFD/BEM Methodology For Launch Pad Acoustic Environment Prediction. In: 18th AIAA/CEAS Aeroacoustics Conference (33rd AIAA Aeroacoustics Conference). AIAA-2012-2107. [102] Van Leer, B., 1976. Muscl, a new approach to numerical gas dynamics. In: Computing in Plasma Physics and Astrophysics. Vol. 1. [103] Van Leer, B., 1977. Towards the ultimate conservative difference scheme. iv. a new approach to numerical convection. Journal of computational physics 23 (3), 276–299. [104] Van Leer, B., Lee, W.-T., Roe, P. L., Powell, K. G., Tai, C.-H., 1992. Design of optimally smoothing multistage schemes for the euler equations. Communications in applied numerical methods 8 (10), 761–769. [105] Vinokur, M., 1974. Conservation Equations of Gasdynamics in Curvilinear Coordinate Systems. Journal of Computational Physics 14, 105–125. [106] Wright, M. W., White, T., Mangini, N., Oct. 2009. Data Parallel Line Relaxation (DPLR) Code User Manual Acadia Version 4.01.1. NASA Technical Memorandum TM-2009-215388. [107] Zhang, Q., Johansen, H., Colella, P., 2012. A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation. SIAM Journal on Scientific Computing 34 (2), 179–201. [108] Zhao, H., 2005. A fast sweeping method for eikonal equations. Mathematics of computation 74 (250), 603–627.

69