Mina R. Mankbadi, Nicholas J. Georgiadis, and James R. DeBonis

Glenn Research Center, Cleveland, Ohio

Comparison of High-Order and Low-Order

Methods for Large-Eddy Simulation of a

Compressible Shear Layer

NASA/TM—2015-218741

June 2015
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated 

to the advancement of aeronautics and space science. 

The NASA Scientifi c and Technical Information (STI) 

Program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI Program provides access 

to the NASA Technical Report Server—Registered 

(NTRS Reg) and NASA Technical Report Server—

Public (NTRS) thus providing one of the largest 

collections of aeronautical and space science STI in 

the world. Results are published in both non-NASA 

channels and by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 

professional papers, but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or of 

specialized interest, e.g., “quick-release” reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 757-864-6500

• Telephone the NASA STI Information Desk at

 757-864-9658

• Write to:

NASA STI Program

 Mail Stop 148

 NASA Langley Research Center

 Hampton, VA 23681-2199
Mina R. Mankbadi, Nicholas J. Georgiadis, and James R. DeBonis

Glenn Research Center, Cleveland, Ohio

Comparison of High-Order and Low-Order

Methods for Large-Eddy Simulation of a

Compressible Shear Layer

NASA/TM—2015-218741

June 2015

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

15th Aviation Technology, Integration, and Operations Conference

sponsored by the American Institute of Aeronautics and Astronautics

Dallas, Texas, June 22–26, 2015
Acknowledgments

This work was supported by the Rotary Wing project and the Aeronautical Sciences project under the Fundamental Aeronautics 

Program. Special thanks to D.A. Yoder and A.A. Ameri for reviewing this report.

Available from

Level of Review: This material has been technically reviewed by technical management. 

This report is a preprint of a paper intended for presentation at a conference. 

Because changes may be made before formal publication, this preprint is made available 

with the understanding that it will not be cited or reproduced without 

the permission of the author.

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

National Technical Information Service

5285 Port Royal Road

Springfi eld, VA 22161

703-605-6000

This report is available in electronic form at http://www.sti.nasa.gov/ and http://ntrs.nasa.gov/

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.
NASA/TM—2015-218741 1 

Comparison of High-Order and Low-Order Methods for 

Large-Eddy Simulation of a Compressible Shear Layer 

Mina R. Mankbadi, Nicholas J. Georgiadis, and James R. DeBonis 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract 

The objective of this work is to compare a high-order solver with a low-order solver for performing 

Large-Eddy Simulations (LES) of a compressible mixing layer. The high-order method is the 

Wave-Resolving LES (WRLES) solver employing a Dispersion Relation Preserving (DRP) scheme. The 

low-order solver is the Wind-US code, which employs the second-order Roe Physical scheme. Both 

solvers are used to perform LES of the turbulent mixing between two supersonic streams at a convective 

Mach number of 0.46. The high-order and low-order methods are evaluated at two different levels of grid 

resolution. For a fine grid resolution, the low-order method produces a very similar solution to the highorder method. At this fine resolution the effects of numerical scheme, subgrid scale modeling, and 

filtering were found to be negligible. Both methods predict turbulent stresses that are in reasonable 

agreement with experimental data. However, when the grid resolution is coarsened, the difference 

between the two solvers becomes apparent. The low-order method deviates from experimental results 

when the resolution is no longer adequate. The high-order DRP solution shows minimal grid dependence. 

The effects of subgrid scale modeling and spatial filtering were found to be negligible at both resolutions. 

For the high-order solver on the fine mesh, a parametric study of the spanwise width was conducted to 

determine its effect on solution accuracy. An insufficient spanwise width was found to impose an 

artificial spanwise mode and limit the resolved spanwise modes. We estimate that the spanwise depth 

needs to be 2.5 times larger than the largest coherent structures to capture the largest spanwise mode and 

accurately predict turbulent mixing.

Nomenclature 

ܽ speed of sound 

ܾ mixing layer thickness 

ܧ total energy per unit mass 

݂ frequency 

ܪ tunnel height 

݇ turbulent kinetic energy 

ܮ domain length 

ܯ ௖convective Mach number 

ܰ number of grid points 

ܲ pressure 

ݍ heat flux 

ܴ݁ Reynolds number 

ܶ temperature 

 time ݐ

ܷ mean velocity 

Δܷ ൌ ܷଶെܷଵ reference velocity 
NASA/TM—2015-218741 2 

ܷ∗ dimensionless mean velocity 

ݑ ௜෤Favre filtered velocity 

ܷഥ mean streamwise velocity 

ݑ′ݑ′ തതതതത Reynolds normal stress 

ݑ′ݒ′ തതതതത Reynolds shear stress 

ݑ௠௥௦ ൌ ඥݑᇱݑ തതതതതതᇱ rms streamwise velocity 

ݒ௠௥௦ ൌ ඥݒᇱݒ തതതതതതᇱ

 rms transverse velocity

ݓ௠௥௦ ൌ ඥത

തതᇱݓതݓതതതᇱ rms spanwise velocity

ݔ ௜Cartesian coordinates 

directions Cartesian ݖ ,ݕ ,ݔ

߬̅௜௝ spatially filtered shear stress 

߬௜௝

ௌீௌ subgrid scale shear stress 

ߩ̅ spatially filtered density 

ߜ ∗displacement thickness 

ߤ ௟laminar dynamic viscosity 

 reference velocity 

ߤ௧ turbulent dynamic

Superscripts: 

~ Favre-filtered 

െ Spatially-filtered 

ܵܩܵ subgrid scale 

Subscripts: 

ݐ total condition 

1 top stream 

2 bottom stream 

I. Introduction 

Use of the Reynolds Averaged Navier-Stokes (RANS) equations has become standard practice for 

aerodynamic analysis. However the Reynolds averaging procedure introduces additional correlation terms 

that must be modeled. These models work reasonably well for the equilibrium wall-bounded and free 

shear layer flows for which they were calibrated, but often require retuning of model coefficients for other 

flows. These models also have severe limitations in predicting unsteady, separated, or transitioning flow. 

Increasing computational throughput has enabled the use of more sophisticated methods such as Direct 

Numerical Simulations (DNS) and Large-Eddy Simulations (LES). These tools are more deterministic 

and, unlike RANS, are less reliant on tunable constants. Consequently, they offer promising physical 

insights that RANS methods cannot. Despite this, there remains a great many issues to ensure the 

accuracy of LES. 

The primary objective of this work is to compare the benefits and shortcomings of a high-order 

method and a low-order method for performing LES of a compressible mixing layer. The high-order 

solver is the Wave-Resolving LES (WRLES) code (Refs. 1 to 3) that employs a Dispersion Relation 

Preserving (DRP) scheme. The low-order solver is the Wind-US code (Ref. 4) which employs the 

second-order Roe scheme. It is well known that second-order schemes are more convenient to use for 

complicated geometries due to their small stencil size. Additionally, second-order schemes can be easily 
NASA/TM—2015-218741 3 

adapted to unstructured grids whereas high-order schemes frequently are limited to structured grids and 

simple geometries. The advantage of high-order methods is that they can achieve the same level of 

accuracy using significantly fewer grid points compared to the low-order schemes. The work herein sets 

to examine whether low-order solvers, given sufficiently fine enough resolution, can produce the same 

level of accuracy as high-order methods. This effort also attempts to address the current state-of-the-art 

practice in performing LES of a compressible mixing layer. 

Several technologies stand to benefit from an enhanced understanding of compressible mixing in free 

shear flows. Examples include supersonic combustion, exhaust nozzles, and internal flows present in jet 

engines and scramjets, to mention only a few. High speed turbulent flows are often investigated using 

Reynolds Averaged Navier-Stokes (RANS) models, and at times with dilitation-based compressibility 

corrections. It is widely known, however, that these compressibility corrections are empirical at best and 

do not replicate the actual physics of compressible turbulent mixing layers (Refs. 5 and 6). Therefore, the 

longer term objectives of this work are to better understand the fundamentals of compressible turbulent 

free shear flows using LES and to apply that insight towards improving existing RANS turbulence models 

for engineering applications. 

Brown & Roshko’s (Ref. 7) classical work on low-speed mixing layers made a key discovery in its 

noteworthy evidence of coherent large-scale spanwise oriented vortical structures. The formation of these 

large-scale “rollers” is triggered by the Kelvin-Helmholtz instability. As the spanwise rollers convect 

downstream they combine through pairing and tearing processes. This work led to a series of 

experimental studies investigating compressibility effects such as Chinzei et al. (Ref. 8), Papamoschou & 

Roshko (Ref. 9), Goebel and Dutton (Ref. 10), Samimy & Elliot (Ref. 11), Hall et al. (Ref. 12), and 

Clemens & Mungal (Ref. 13). The experimental studies observed that the most significant effect of 

compressibility on a turbulent mixing layer is suppression of the growth rate relative to an incompressible 

mixing layer at the same velocity and density ratios. The experiments demonstrate that compressibility 

effects are greatest between convective Mach numbers of 0.5 and 1. The particular case examined here is 

the turbulent mixing between two supersonic streams at a moderate convective Mach number of 0.46. 

Similarly, numerical simulations have been performed by several researchers including Foysi and Sarkar 

(Ref. 14), Freund et al. (Ref. 15), and Pantano and Sarkar (Ref. 16). The numerical simulations showed 

that dilatation terms are relatively small and not affected by increasing convective Mach number, whereas 

the pressure-strain terms in the Reynolds stress transport decrease monotonically with increasing 

convective Mach number. The redistribution of the Reynolds stresses in turn lowers turbulent production, 

which results in reduced mixing layer growth rate. 

In this work, LES is employed to numerically investigate a reference turbulent mixing layer at 

moderate convective Mach number, 0.46. Two different flow solvers are used. The low-order solver 

employs a second-order upwind scheme while the high-order method implements Bogey & Bailly’s 

(Refs. 17 and 18) 11-point DRP scheme. The effects of numerical scheme, subgrid scale (SGS) modeling, 

and filtering are also investigated. Furthermore, a previous hybrid RANS-LES work concluded that an 

insufficient spanwise domain limited the development of the full three-dimensional structures (Ref. 19). 

Therefore, the work herein, in addition to its primary purpose of comparing high-order and low-order 

methods, also seeks to analyze the effects of increasing the spanwise extent while maintaining the 

spanwise resolution. This type of parametric study was recently used with the same low-order technique 

to investigate sensitivity of separated bluff body flows to spanwise width (Ref. 20). 

Ideally, LES would be performed not only in the mixing region but also along the upstream channels 

to resolve their turbulent boundary layers. Recently, Inoue & Pullin (Ref. 21) investigated the zero 

pressure gradient turbulent boundary layer through LES over momentum thickness Reynolds numbers of 

ܴ݁ఏ ൌ ܱሺ10ଷሻ െ ܱሺ10ଵଶሻ. This required a large number of grid points to resolve the boundary layer. In 

the present work focusing on the shear layer, the Reynolds numbers were ܴ݁ఏ ൌ 5.1	x	10଺ and ܴ݁ఏ ൌ

4.7	x	10଺ for the top and bottom channel flows, respectively. Fully resolving these turbulent boundary 

layers while also performing LES in the mixing region is computationally prohibitive. Therefore, the LES 

simulations are designed to focus on resolving the shear layer itself. This will also enable a side-by-side 

comparison between the high-order method WRLES with the low-order method Wind-US. 
NASA/TM—2015-218741 4 

The paper is organized as follows: A description of the flow is provided in Section II. In Section III 

the formulation is discussed, while the details of the computational approach are given in section IV. A 

brief summary of the low-order solver Wind-US and the high-order solver WRLES is provided. Details of 

the grid, boundary conditions, and experimental conditions are also presented. Comparisons of the highorder and low-order method are presented in Section V. The section begins with a general description of 

the flow and then considers specifics of numerical method, SGS modeling, and filtering. The importance 

of spanwise width is highlighted. The results section finishes with a discussion comparing the low-order 

and high-order method again but at coarser grid resolutions. Conclusions are presented in Section VI. 

II. Description of the Flow 

The current simulations recreate case 2 of the experimental mixing layer tested by Goebel and Dutton 

(Ref. 10). The schematic in Figure 1 shows two supersonic streams separated by a thin splitter plate 

exiting into a long mixing section. The tunnel has a height of 48 mm and a width of 96 mm. The splitter 

tip separating the two flows has a thickness of 0.5 mm (an important length scale affecting vortex 

shedding) and the entrances of both streams into the mixing section have equal heights of 24 mm. The 

high-speed flow is located on the top with a Mach number of 1.91, while the bottom flow is at Mach 1.36. 

The convective Mach number calculated from the flow conditions is ܯ ௖ൌ 0.46, where ܯ ௖ൌ ሺܷଵ െ

ܷଶሻ/ሺܽଵ ൅ ܽଶሻ. The two flows are pressure matched at 49kPa. The total temperature of the top flow is 

578 K while the bottom flow is at 295 K. The velocities are ܷଵ ൌ 700	݉/ݏ and ܷଶ ൌ 399	݉/ݏ. 

The experiments used Laser Doppler Velocimetry (LDV) to measure velocities in the mixing section. 

The total viewing length of the optical instrumentation was 500 mm. Transverse velocity profiles at several 

streamwise locations were measured. The Reynolds normal stresses in the streamwise and transverse 

directions were measured; however, the normal stress in the spanwise direction was not measured. The LDV 

system was used to measure properties of the incoming boundary layers at the splitter tip. The boundary 

layer thickness, displacement thickness, and momentum thickness in each of the two streams feeding the 

mixing layer were measured just upstream of the splitter plate trailing edge. The spectra of the turbulent 

kinetic energy and the fundamental frequency of vortex shedding were not measured. 

Figure 1.—Schematic of splitter tip separating the top flow from the bottom flow. 

TABLE 1.—EXPERIMENTAL FLOW CONDITIONS 

Parameter Top flow Bottom flow 

 1.36 1.91 ܯ

ܷ	ሺm/s) 700 399 

ܶ௧	ሺK) 578 295 

ܲ	ሺkPa) 49 49 

ߩ	ሺkg/mଷሻ 0.51 0.79 
NASA/TM—2015-218741 5 

III. Formulation 

The Favre-filtered Navier-Stokes equations are solved numerically with both the high- and low-order 

methods. The viscosity and thermal conductivity are allowed to vary with temperature. Favre filtering is a 

density-weighted spatial filtering procedure that is applied to the time dependent Navier-Stokes equations 

to remove small-scale fluctuations that are too small to be resolved. The continuity, momentum, and 

energy equations are: 

̅ߩ߲

൅ ݐ߲

߲

௜ݔ߲

ሺߩ̅ݑ௜෤ሻ ൌ 0, (1)

ሻ෤௜ݑ̅ߩሺ߲

൅ ݐ߲

߲

௝ݔ߲

ത߲ܲ ൅ ൯෤௝ݑ௜෤ݑ̅ߩ൫

௜ݔ߲

ൌ ߲

௝ݔ߲

ሺ߬̅௜௝ ൅ ߬௜௝

ௌீௌሻ, (2)

෨ܧ߲

൅ ݐ߲

߲

௝ݔ߲

߲ ൌ ത൯ܲ෤௝ݑ ൅ ෨ܧ௝෤ݑ൫

௝ݔ߲

௝௜߬௜෤ݑ ൅ ௝௜߬̅௜෤ݑ൫

ௌீௌ൯ ൅ ߲

௝ݔ߲

௝ݍ ൅ ௝തݍ൫

ௌீௌ൯	. (3)

The density-weighted spatially-filtered velocity is denoted as ݑ௜෤, while the spatially-filtered density and 

pressure are ߩ̅ and	ܲത, respectively. The stress tensor is ߬̅௜௝. For the scales that are too small for the scheme 

and grid to resolve, an optional subgrid scale model specifies ߬௜௝

ௌீௌ and ݍ௝ௌீௌ. 

IV. Computational Approach 

Two different numerical approaches are employed herein with the purpose of investigating whether a 

low-order scheme is capable of reproducing the results of a high-order scheme. High-order methods 

require less grid to resolve turbulent structures, whereas low-order methods require additional grid points 

to achieve the same level of accuracy. However, the large stencil sizes and complex numerics typically 

limit high-order methods to simple geometries. If the low-order method is capable of producing the same 

answer as the high-order method, then the low-order method may be more attractive for other geometries 

that are more complex than the relatively simple shear layer that is the focus of this paper. 

A. Flow Solvers 

1. WRLES

The high-order solver is the Wave-Resolving LES (WRLES) code (Refs. 1 to 3). It is a finite 

difference code which employs the Dispersion Relation Preserving (DRP) scheme developed by Bogey & 

Bailly (Refs. 17 and 18). The stencil size is 11 points and is formally fourth-order accurate. While a 

traditional fourth-order scheme uses fewer points in the stencil, the additional points in the DRP scheme 

are used to minimize the dispersion error. It is an explicit central scheme similar to the classic tenth-order 

central scheme but with weights carefully chosen to minimize the dispersion error. Since it is a central 

differencing scheme, it is not inherently dissipative. A spatial filter, developed specifically to match the 

resolution of the differencing scheme, is used to preserve numerical stability. A coefficient multiplying 

the effect of the filter has been added, so that the magnitude of the dissipation can be controlled. The 

filter’s purpose and effects are strictly to help improve numerical stability. The DRP filter coefficient used 

was 0.5. Temporal discretization is performed using a low-dissipation and low-dispersion forth-order 

Runge-Kutta algorithm summarized in Berland et al. (Ref. 18). Three subgrid scale models have been 

implemented: Smagorinksy’s, Vreman’s, and a dynamic Smagorinsky. In particular, Vreman’s (Ref. 22) 

algebraic SGS model is used herein since it is less dissipative than the Smagorinsky subgrid model and 

not as computational intensive as the dynamic model. 
NASA/TM—2015-218741 6 

WRLES uses the Message Passing Interface (MPI) standard to divide the domain across nodes and 

the Open Multi-Processing (OpenMP) directive to parallelize work on each node. In order to maintain 

order accuracy across interfaces, WRLES uses overlapping, point-matched zone coupling. 

2. Wind-US

The low-order solver is Wind-US, the production flow solver of the NPARC Alliance. While it has 

historically been used primarily in RANS mode, the code may be used in LES or hybrid RANS-LES 

mode. It is capable of working with unstructured meshes as well as curvilinear structured meshes. It has 

several choices for temporal and spatial discretization as well as different options for flux splitting that are 

well documented (Ref. 4). For this particular problem, the numerical scheme used in the Wind-US CFD 

code is the second-order Roe Physical upwinding finite-volume scheme for spatial discretization and an 

implicit first-order time-stepping scheme. The Roe Physical scheme is based upon the original Roe 

scheme but modified for stretched grids. Note that while implicit temporal schemes require a matrix 

inversion, they are much more stable than explicit time schemes. When run in LES mode, Wind-US 

employs the Implicit (or Numerical) LES (ILES) philosophy. As discussed by Pope (Ref. 23), this 

approach assumes that the numerical scheme and grid act like a subgrid scale model, and so no explicit 

SGS model is used. 

For parallelization, Wind-US was run with the PVM (Parallel Virtual Machine) directive to create a 

parallel computing system. Wind-US uses abutting, point-matched zone coupling that is second-order 

accurate across interfaces. Within Wind-US, the MPI directive is also available for parallel computing. 

B. Grid 

The baseline grid used is structured as shown in Figure 2. The baseline mesh uses 14.4 million points 

to resolve the LES mixing section and has a thin spanwise width of 6 mm. Forty one points resolve the 

0.5 mm splitter tip yielding a Δݕ ൌ Δݔ ൌ 0.0125	݉݉ in this region. Grid details such as the number of 

points, domain size, and minimum and maximum grid spacing are provided in Table 2 for the mixing 

region. The total length of the mixing section is ܮ ൌ 500	݉݉ as may be observed in Figure 2. 

Downstream of the splitter plate towards the outflow the grid is stretched in the x-direction. The extent of 

grid stretching is deduced by comparing the minimum and maximum grid spacing. 

Figure 2.—Boundary conditions and dimensions of the LES mixing section (only every 8th point is plotted). 

TABLE 2.—GRID PROPERTIES OF THE MIXING SECTION 

Parameter Direction 

ݖ ݕ ݔ

Number of points 1025 425 33

Length of domain ሺ݉݉ሻ 500 48 6

Minimum spacing ሺ݉݉ሻ 0.0125 0.0125 0.1875

Maximum spacing ሺ݉݉ሻ 5 0.2 0.1875

 
NASA/TM—2015-218741 7 

C. Boundary Conditions 

1. General Boundary Treatment

The flow conditions are taken from the experiments of Goebel & Dutton (Ref. 10). The flow in the 

computational domain is supersonic throughout, thereby simplifying inflow and outflow boundary 

treatments (Figure 2). Inflows are fixed and the outflow is extrapolated. The upper and lower tunnel walls 

are assumed to be sufficiently far away from the shear layer such that they do not affect the mixing and 

are set as inviscid walls. Periodic boundaries are imposed on the spanwise sides of the domain. Their 

purpose is to simulate an infinite span and avoid the cost of modeling the sidewalls. The experiments had 

a finite span of 96 mm. As mentioned previously and discussed in greater detail later in this paper, the 

width of the spanwise dimension was investigated as a significant part of this work. 

2. Inflow Boundary Treatment 

In the experiment, turbulent boundary layers were present on the splitter plate upstream of the mixing 

section. Dutton & Goebel (Ref. 10) provided the values of the boundary layer thickness, displacement 

thickness, and momentum thickness. However, the actual velocity profiles and turbulent intensities were 

not reported. The best method for simulating these turbulent boundary layers would be a complete LES of 

the upstream region, producing unsteady turbulent boundary layers that match the experiment's integral 

thickness parameters. However, this type of analysis is cost prohibitive given the extremely small 

turbulent scales that would exist in a boundary layer at these Reynolds numbers. 

In this study, the inflow conditions were obtained from a separate two-dimensional RANS calculation 

of the mixing layer and upstream channels using the Menter Shear Stress Transport (SST) turbulence 

model. To do so, the lengths of the upstream ducts were varied until they produced RANS profiles that 

matched the experimental values of the displacement thicknesses. Table 3 compares the RANS-obtained 

profiles with the experiment by listing the computed values of displacement thickness ߜ∗, mean 

streamwise velocity ܷ, and the area-averaged pressure P. Notice that the RANS solution indicates that a 

pressure-matched flow of 49 kPa is being delivered as in the experiments of Goebel-Dutton (Ref. 10). The 

incoming velocity profiles are depicted in Figure 3 and have displacement thicknesses matching the 

TABLE 3.—PROPERTIES OF THE INCOMING FLOW DELIVERED TO THE LES MIXING SECTION 

Solution Top flow Bottom flow 

݇ܲܽ ܲ, ݏ݉/ ܷ, ݉݉ ∗ߜ ݇ܲܽ ܲ, ݏ݉/ ܷ, ݉݉	∗ߜ

RANS upstream 0.9318 700.3 48.898 0.4681 398.8 48.940 

Experiment 0.90 700 49 0.44 399 49 

Figure 3.—Velocity profile of the flow entering the shear layer. 
NASA/TM—2015-218741 8 

experiments to within 3.5 percent for the top flow and within 6 percent for the bottom flow. The resulting 

RANS mean-flow profiles were used as fixed conditions at the inflow plane to the mixing section. A 

limitation of this boundary treatment is that the RANS method only provides a steady state profile to an 

unsteady LES mixing region (Ref. 24). 

D. Time Averaging and Sampling Procedure 

To ensure numerical stability, the maximum time-step for the implicit low-order solver was found to 

be	݀ݐ ൌ 0.05ݏߤ, while the explicit high-order solver was more stringent and the maximum time step was 

݀ݐ ൌ 0.01ݏߤ. At the beginning of the simulations, there is a transient wave that travels through the 

domain. This originates from a normal shock produced at solution start-up. The time for all initial 

transient effects to be eliminated and before statistics can be gathered is no less than 5.5 Flow Thru Times 

(FTTs). The Flow Thru Time is defined as the time it takes a particle inside the shear layer to travel from 

the splitter tip to exit: 1	ܨܶܶ ൌ ܮ/ܷ௏஺ீ ൌ 0.51 ms. ܷ஺௏ீ is an estimate of the speed inside the shear 

layer taken to be the average velocity of the top and bottom streams. Extending the start-up wait time 

until 9.5 FTTs was found to be unnecessary, since statistics gathered after employing these two different 

start-up times were identical. 

The effect of sampling duration was examined by gathering second-order statistics over periods of 

9.5 and 33 FTTs that were then time averaged. Note that both of those sampling durations are 

considerably long. The resulting rms plots were compared and it was found that extending the sampling 

duration from 9.5 to 33 FTTs did not affect the rms values. Turbulent statistics were gathered at three 

frequencies corresponding to every 5, 1, and 0.4 ݏߤ. Those frequencies were all extremely high and the 

results were found to be identical. A similar study varying sampling duration and frequency was 

performed by Mankbadi & Georgiadis (Ref. 20). After the turbulent statistics have been gathered, the 

power spectra is obtained by sampling the data every 0.05	ݏߤ after which a fast Fourier transform is 

performed to identify the fundamental vortex shedding frequency. 

V. Results and Discussion 

All the computations were performed on the NASA Pleiades supercomputing platform which is ranked 

amongst the top ten U.S. supercomputers. Pleiades is made up of tens of thousands of nodes capable of 

4.49 petaflop/s. Each node has two processors and each processor has either 6, 8, 10, or 12 cores. A typical 

simulation adhered to the following procedure: (1) initialize the flow field and wait for the start-up transient 

wave to exit the domain, (2) begin collecting data required for computing time-averaged turbulent 

intensities, and (3) at a much higher sampling rate (about one hundred times the fundamental) gather data to 

compute the power spectra of the flow. This procedure was used for both flow solvers. Wind-US 

simulations typically employed 39 nodes with a total of 468 cores and the turnaround time for the entire 

procedure was about a week. WRLES utilized 129 nodes with a total of 1536 cores and also had a 

turnaround time of about a week. 

Figure 4 shows dimensionless density contours of the mixing layer obtained using the two flow solvers 

for the baseline grid with 6 mm spanwise width. The high-order method, WRLES, is compared against the 

low-order method, Wind-US. Note that the size of the structures grows as the vortices pair and propagate 

downstream. There are more organized structures that can be seen in the WRLES plot whereas in the 

Wind-US plot the structures appear to be smeared (Figure 4(A) and (B)). The shock structure also slightly 

differs between the two. There are more reflections present in Wind-US than with WRLES. Wind-US’s 

small stencil and Riemann solver are adept at shock capturing whereas WRLES reduces the order to avoid 

smearing the shock with its larger stencil. The von Karman shedding is evident in Figure 4(C) and (D) 

where vortices are being produced just downstream of the splitter tip. The close-up views of the splitter tip 

demonstrate that the high-order method better resolves vortex shedding than the low-order method. Further 

downstream at ݔ ൌ 100	݉݉, the Wind-US solution appears to maintain the regular von Karman pattern 

while the WRLES solution appears to have broken down into a more irregular pattern. 
NASA/TM—2015-218741 9 

Figure 4.—Dimensionless density contours of the mixing section for (A) the high-order solver WRLES and 

(B) low-order solver Wind-US. Close up views of the splitter tip showing von Karman shedding for (C) WRLES 

and (D) Wind-US. 

 

(A) 

WRLES 

(B) 

Wind-US 

Wind-US 

(D) 

WRLES 

(C) 
NASA/TM—2015-218741 10 

A. Flow Solver Comparison 

1. Effects of Flow Solvers

Streamwise mean and rms velocities from the high-order method WRLES and the low-order method 

Wind-US are compared in Figure 5 (the transverse and spanwise rms velocities are presented in 

subsection B). Recall that the Wind-US approach employed a second-order finite volume solver and 

WRLES is a formally fourth-order accurate finite difference method employing an eleven-point DRP 

stencil. Both methods employ the ILES approach where no subgrid scale model is used. Despite the 

significant differences in numerical methodology, notice the relatively close agreement between the two 

numerical methods shown in Figure 5. The rms profiles at three different axial stations are very similar 

(Figure 5D to F), and the mean velocity profiles match fairly well (Figure 5(B) and (C)). The largest 

difference in mean velocity is in the wake at ݔ ൌ 50݉݉ in Figure 5(A). 

Figure 6 compares plots of the transverse velocity spectra of both numerical methods. The frequency 

is scaled by the fundamental frequency of vortex shedding ݂௩௦. The vortex shedding frequency, 

݂௩௦,	predicted is 124.51 and 134.28 kHz by Wind-US and WRLES, respectively. The data from the WindUS solution shows a sub-harmonic exists at ݂/݂௩௦ ൌ 0.5. The vortices produced at the tip recombine 

further downstream but at a lower frequency than that at which they were produced (Ref. 25). The 

vortices coalesce at the frequency of the sub-harmonic, ݂/݂௩௦ ൌ 0.5. The sub-harmonic peak is not only 

present very close to the tip at ݔ ൌ 10݉݉ but persists until ݔ ൌ 100݉݉. The fundamental frequency is 

the strongest at ݔ ൌ 10݉݉ close to the tip and diminishes downstream, until it has disappeared by 

 100݉݉. ൌ ݔ

The sub-harmonic cannot be seen in the spectra coming from WRLES. In addition, the spectra 

coming from WRLES shows that all the low frequency peaks have been damped out by ݔ ൌ 50݉݉ and 

there is more energy in the smaller-scale, higher-frequency structures. The high-order numerics in 

WRLES can resolve smaller scales of turbulence, and these scales may be responsible for damping out 

the large-scale structures that persist in the Wind-US solution. Because of the differences in the 

fundamental frequency predicted and overall spectra, it is a little counterintuitive that two numerical 

methodologies would agree to such an extent in terms of the velocity and rms profiles shown in Figure 5. 

2. Effects of SGS Model

In the previous section, the effects of the numerical scheme used by the high-order and low-order 

methods were examined. This section examines the effects of subgrid scale modeling on the 11 point 

DRP scheme solution. An additional case was computed with WRLES using Vreman’s subgrid model. 

The WRLES with SGS and WRLES without SGS cases are compared in Figure 7 and Figure 8. As 

evident from the data, the mean profile and rms profile are unaffected. The shedding frequency predicted 

by WRLES is 134.28 kHz for both cases and only minor differences are seen in the spectra. Therefore, 

the effect of SGS modeling is negligible for the 14.4 million point grid. In Figure 9, a plot of the SGS 

viscosity normalized by the laminar viscosity is provided. It is clear that the subgrid scale model is only 

active within the shear layer. Closer inspection of the scale reveals that ߤ௧/ߤ does not exceed unity 

anywhere in the mixing section. This explains why SGS effects are negligible for this problem. The grid 

is fine enough such that the vast majority of turbulence is computed directly and the SGS model does not 

add a significant amount of dissipation. If the grid was coarsened, then perhaps the SGS model would 

come into play. The fact that ߤ௧/ߤ does not exceed unity also helps explain why the second-order method 

performed very well. The grid was fine enough in the mixing region, making any benefit of a higher order 

scheme irrelevant. The coarser grid and SGS effects will be explored in the next subsection. 
NASA/TM—2015-218741 11 

Figure 5.—Comparison between the low-order solver Wind-US and the high-order solver WRLES on the 

14.4 million grid. For three different axial stations, the shear layer’s mean streamwise velocity (A) at 

ݔ ൌ 50݉݉, (B) at ݔ ൌ 100݉݉, and (C) at ݔ ൌ 150݉݉. The rms streamwise velocity (D) at ݔ ൌ 50݉݉, 

(E) at ݔ ൌ 100݉݉, and (F) at ݔ ൌ 150݉݉. 

(A) (D) 

(B) (E) 

(C) (F) 
NASA/TM—2015-218741 12 

Figure 6.—Spectra of the transverse velocity for both 

numerical methods on the 14.4 million point grid. 
NASA/TM—2015-218741 13 

Figure 7.—Effect of subgrid scale modeling on the shear layer’s streamwise mean and rms velocity at three 

different locations. The streamwise velocity (A) at ݔ ൌ 50݉݉, (B) at ݔ ൌ 100݉݉, and (C) at ݔ ൌ 150݉݉. 

The rms streamwise velocity (D) at ݔ ൌ 50݉݉, (E) at ݔ ൌ 100݉݉, and (F) at ݔ ൌ 150݉݉. 

(A) (D) 

(B) (E) 

(C) (F) 
NASA/TM—2015-218741 14 

Figure 8.—Effects of SGS modeling on the spectra 

of the transverse velocity. 

Figure 9.—The turbulent subgrid scale viscosity normalized by the laminar viscosity for the 14.4 million point grid. 

 
NASA/TM—2015-218741 15 

3. Effect of Grid Density

Recall in Section A-1, we demonstrated that the low-order solver could produce similar mean flow 

and turbulent statistics to the high-order solver on a fine grid. This section compares the two codes on a 

coarser grid. The coarsened grid used every other point in the streamwise and transverse directions while 

retaining the same spanwise width and points. 

The results of the grid study using the high-order solver WRLES are presented in Figure 10. Overall, 

the high-order solver performs nearly as well on the coarser grid. Close to the splitter tip, even the peak 

streamwise rms velocity values differ very little between the two resolutions. However, the mean flow 

predictions in the wake of the splitter tip are sensitive to resolution. Further away from the splitter tip, the 

mean flow and turbulent intensities predicted by the 3.6 million point grid are almost identical to the fine 

14.4 million point grid. Figure 10 also demonstrates that effects of SGS modeling are still very minor, 

despite the coarser grid employed. Comparing the SGS viscosities of Figure 9 for the finer grid to those of 

the coarsened grid in Figure 11, demonstrates that as the resolution was coarsened the SGS viscosity 

increases, but is still on the same order as the laminar viscosity. 

The results of the grid study using the low-order solver Wind-US are presented in Figure 12. These 

results demonstrate that the low-order method changes significantly at the coarse resolution. The loworder method cannot resolve the complex flow physics of vortex shedding, mixing, and shear layer 

growth. The inherent dissipation of the low-order scheme dampens out the instabilities and inhibits proper 

turbulent growth. There are substantial differences in both the mean flow velocities and the rms 

quantities. Not only is the peak rms value diminished but the extent of the profile across the shear layer is 

reduced significantly from that observed with the higher resolution. In other words, very little turbulent 

mixing is occurring. 

Figure 13 shows the transverse velocity spectra comparing the high-order and low-order methods at 

the coarse resolution. Note that the high-order method has well-defined sharp peaks at the fundamental 

frequency and its multiples. However, the low-order method has a blunt peak at the fundamental 

frequency. Recall the fundamental frequency predicted by WRLES on the fine grid was 134.28 kHz, here 

for the coarse grid it is higher at 141.91. Likewise recall Wind-US predicted a fundamental frequency of 

124.51 kHz, for the coarser grid it is also higher 128.17 kHz. 

These findings demonstrate the expected result that high-order algorithms need less resolution than 

low-order methods to achieve the same level of accuracy. Unlike high-order methods, low-order methods 

are expected to fail at coarse resolutions and indeed did in this study. The increase in accuracy that highorder methods offer means that high Reynolds number problems become tractable at more reasonable grid 

sizes. With regards to computational wait time, high-order methods typically incorporate multi-stage time 

stepping and use large stencils that drives up the computational cost. Nevertheless, in many cases highorder methods can deliver that same level of accuracy at a cheaper computational cost since less grid is 

required. However, for the current problem, the two methods’ computational costs were about the same. 

This is primarily due to two reasons: (1) the explicit time-stepping of the high-order code requires a 

smaller time step than the implicit time stepping of the low-order method, and (2) the low-order method 

needs to regularly issue a checkpoint where the entire data is written to disk then an external process 

strips down the file to keep only the desired data. If not for this limitation in I/O performance, the loworder method would likely have been much faster than the higher-order method for the same grid. 
NASA/TM—2015-218741 16 

Figure 10.—Grid study using WRLES with a spanwise width of 6mm. The mean streamwise velocity (A) at 

ݔ ൌ 50݉݉, (B) at ݔ ൌ 100݉݉, and (C) at ݔ ൌ 150݉݉. The rms streamwise velocity (D) at ݔ ൌ 50݉݉, 

(E) at ݔ ൌ 100݉݉, and (F) at ݔ ൌ 150݉݉. 

 

(A) (D) 

(B) (E) 

(C) (F) 

WRLES WRLES
NASA/TM—2015-218741 17 

Figure 11.—The turbulent subgrid scale viscosity normalized by the laminar viscosity for the 3.6 million point grid. 
NASA/TM—2015-218741 18 

Figure 12.—Grid study using Wind-US with a spanwise width of 6mm. The mean streamwise velocity 

(A) at ݔ ൌ 50݉݉, (B) at ݔ ൌ 100݉݉, and (C) at ݔ ൌ 150݉݉. The rms streamwise velocity (D) at 

ݔ ൌ 50݉݉, (E) at ݔ ൌ 100݉݉, and (F) at ݔ ൌ 150݉݉. 

 

(A) (D) 

(B) (E) 

(C) (F) 

Wind-US Wind-US 
NASA/TM—2015-218741 19 

Figure 13.—Spectra of the transverse velocity for both 

numerical methods on the 3.6 million point grid. 

B. Effect of Spanwise Width 

Periodic boundary conditions are imposed in the spanwise direction to simulate an infinite span and 

avoid modeling the sidewalls. Selection of too small of a spanwise width cannot guarantee accurate 

capturing of the three-dimensional nature of vortex dynamics. To avoid imposing a nonphysical spanwise 

mode, the spanwise extent and resolution must be large enough to capture the wavelengths of the largest 

turbulent structures in the domain. 

An insufficient spanwise depth has been shown to artificially dampen the spanwise formation of 

vortices and resemble a two-dimensional instability (Ref. 20). In some cases, a two-dimensional 

instability predicts higher growth rates than a three-dimensional instability (Ref. 26). It is necessary to 

have a sufficient spanwise width to allow for the formation of three-dimensional coherent structures. Four 

simulations are performed here using the high-order method, with various spanwise widths in order to 

determine if: (1) artificial suppression is present when the spanwise width is smaller than the length scale 

of the largest coherent structure, (2) this artificial suppression is eliminated as the spanwise width is 

extended to be several times larger than the dominant length scale, and (3) as the mixing layer grows 

downstream, the size of the coherent structure grows thereby requiring a larger spanwise depth to be 

resolved. 

WRLES with the SGS model was used to study the effect of the spanwise width. Table 4 shows the 

grid resolution as the spanwise domain is doubled three times consecutively. Georgiadis et al. (Ref. 19) 

simulated this same flow and noted significant differences in rms statistics when examining much smaller 

spanwise depths ranging from 1݉݉ to 6݉݉. The baseline grid was chosen to have the largest spanwise 

depth considered by Georgiadis et al. (Ref. 19) of 6݉݉ and had a spanwise grid spacing of Δܼ ൌ

6݉݉/32. Keeping the spanwise grid spacing fixed, the spanwise width was extended up to 48݉݉, 

utilizing a total of 112 million grid points. The peak values of ݓ௠௥௦ were also obtained. The streamwise 

and transverse resolution was unchanged. 
NASA/TM—2015-218741 20 

TABLE 4.—SPANWISE WIDTH STUDY USING WRLES WITH SGS 

Spanwise width Resolution Grid points (million) ܯܽݔ ܹݎ݉ݏ/Δܷ

150݉݉ 100݉݉ 50݉݉ ൌ ݔ

6mm 1025  425  33 14.4 0.0947 0.09343 0.07341 

12mm 1025  425  65 28.3 0.1005 0.1248 0.1139 

24mm 1025  425  129 56.2 0.1056 0.1362 0.1315 

48mm 1025  425  257 112.0 0.1106 0.1356 0.1253 

Figure 14.—Effect of spanwise width on the dimensionless streamwise vorticity. 

Figure 14 demonstrates the effects of increasing the spanwise width by plotting the streamwise 

vorticity at the three axial stations for spanwise widths of 6, 12, 24, and 48݉݉. Evidently the size of the 

structures is small at ݔ ൌ 50݉݉, but as they propagate downstream the structures grow in size and by 

ݔ ൌ 150݉݉ the structures have grown considerably. The effect of increasing the spanwise width is most 

easily observed at ݔ ൌ 150݉݉. For a spanwise width of 6mm, the development of the structures is 

suppressed, but as the spanwise width is increased to 48mm, the structures are allowed to develop and 

thus appear to be bigger. The size of the turbulent structures relative to the mixing layer thickness can be 

seen by comparison with the value of b that is shown to the right of Figure 14. 

Although the mean flow is not affected by spanwise widths, the turbulence intensities are sensitive. 

The turbulence intensities in all three directions and turbulent kinetic energy are compared for all four 

spanwise widths at ݔ ൌ 50݉݉ (Figure 15), ݔ ൌ 100݉݉ (Figure 16), and ݔ ൌ 150݉݉ (Figure 17). At 

axial stations of 50 and 100mm, the peak values of ݑ௠௥௦ and ݒ௠௥௦ are almost unchanged as the spanwise 

width was gradually increased from 6mm to 48mm. At ݔ ൌ 150݉݉ (Figure 17(A) and (B)), the peak 

value in ݒ௠௥௦ is also unchanged but the peak ݑ௠௥௦ value slightly increases. From examining ݑ௠௥௦ and 

ݒ௠௥௦, the need to extend the spanwise domain is not apparent since ݑ௠௥௦ and ݒ௠௥௦ are not changing 

appreciably. However, close inspection of ݓ௠௥௦ and ݇ shows that ݓ௠௥௦ is artificially being suppressed for 

the 6mm case (Figure 15(C) and (D), Figure 16(C) and (D), Figure 17(C) and (D)). Expanding the width 

of the domain increases the peak ݓ௠௥௦ at the ݔ ൌ 50݉݉ and ݔ ൌ 100݉݉ locations and ݇

correspondingly increases. This differs from a similar study of spanwise width for a square cylinder 

computation (Ref. 20). For that problem ݓ௠௥௦ also increased with an increase in spanwise width but the ݇

remained constant, forcing a rebalance of the energy between ݑ௠௥௦, ݒ௠௥௦, ܽ݊݀	ݓ௠௥௦. 

Note that at ݔ ൌ 50݉݉, where the vortex size is small, the effect of the width on ݓ௠௥௦ is relatively 

small. However, as the turbulent structures increase in size moving further downstream, a larger spanwise 

width is required. Hence, at ݔ ൌ 100݉݉, note that the peak rms value is identical between the 24 and 

48mm case. Further downstream at ݔ ൌ 150݉݉, the peak value of ݓ௠௥௦ of the 48mm case lies between 

the 12 and 24mm cases. This demonstrates that extending the spanwise width eliminates artificial 

suppression of the spanwise rms value. 

Mixing

Layer

Thickness

4.4݉݉

ܾ ൌ

ܾ ൌ

X = 50mm

X = 100mm

X = 150mm

Spanwise Width48mm

Spanwise Width48mm

Spanwise Width48mm

6mm 12mm 24mm

6mm 12mm 24mm

6mm 12mm 24mm

Streamwise Vorticity

2.5݉݉

3.5݉݉

ܾ ൌ
NASA/TM—2015-218741 21 

Figure 15.—At ݔ ൌ 50݉݉, the effect of spanwise width on the shear layer’s (A) streamwise rms velocity, 

(B) transverse rms velocity, (C) spanwise rms velocity, and (D) the turbulent kinetic energy. 

 

(A) (B) 

(C) (D) 
NASA/TM—2015-218741 22 

Figure 16.—At ݔ ൌ 100݉݉, the effect of spanwise width on the shear layer’s (A) streamwise rms velocity, 

(B) transverse rms velocity, (C) spanwise rms velocity, and (D) the turbulent kinetic energy. 

In the aforecited work of Mankbadi & Georgiadis (Ref. 20) on a square cylinder, a spanwise depth ten 

times greater than the cylinder’s diameter demonstrated the ݓ௠௥௦ convergence. For a cylinder at low 

Reynolds number, the shed vortices are about the size of the cylinder’s diameter. However, for the shear 

layer examined here, the relevant length scale is the shear layer width, which grows downstream. The 

Kelvin-Helmholtz instability sets off the vortex shedding process that then becomes turbulent mixing 

downstream. As a consequence, the size of the largest structure increases axially. Hence, at ݔ ൌ 50݉݉

where the length scale is small, a spanwise depth of 6mm is sufficient and spanwise width convergence is 

observed in Figure 15(C). Further downstream at ݔ ൌ 150݉݉, the size of the length scale has grown and 

hence a variation in ݓ௠௥௦ is observed in Figure 17(C). 

The relevant length scale is the shear layer thickness which correlates with the size of the largest 

coherent structure. A closer examination of the shear layer thickness is necessary to assess the ratio 

between the thickness and spanwise width at which spanwise width independence is reached. The shear 

layer thickness was defined by Goebel & Dutton (Ref. 10) as ܾൌݕ଴.ଽ଴ െ ݕ଴.ଵ଴, where ݕ଴.ଽ଴ is the 

location corresponding to ܷ∗ ൌ 0.90. Recall ܷ∗ is normalized such that at the top flow it is one and at the 

bottom flow it is zero. From Table 5, note that at ݔ ൌ 50݉݉	the shear layer thickness is ܾ ൌ 2.5݉݉, 

while it has grown to the shear layer thickness of	4.4݉݉ at ݔ ൌ 150݉݉ (see Table 5). 

(A) (B) 

(C) (D) 
NASA/TM—2015-218741 23 

Figure 17.—At ݔ ൌ 150݉݉, the effect of spanwise width on the shear layer’s (A) streamwise rms velocity, 

(B) transverse rms velocity, (C) spanwise rms velocity, and (D) the turbulent kinetic energy. 

TABLE 5.—SHEAR LAYER THICKNESS 

Location 

ݔ	ሺ݉݉ሻ

LES ሺݖ ൌ 48݉݉ሻ

ܾሺ݉݉ሻ

.݌ݔܧ

ܾሺ݉݉ሻ

50 2.549 3.241 

100 3.480 4.508 

150 4.426 5.284 

Note that the ݑ௠௥௦ for the 6mm case in Figure 15(A) is identical to the 48mm case, thus spanwise 

width independence has been reached at 6mm. When the 6mm width is normalized by the 2.5mm 

thickness at ݔ ൌ 50݉݉	it yields a ratio of 2.3. Likewise in Figure 17(C) the peak ݓ௠௥௦ value for the 

48mm case lies between the 12mm and the 24mm cases, suggesting that they are relatively close at a 

spanwise width of 12mm. When the 12mm width is normalized by the 4.4mm thickness at ݔ ൌ 150݉݉	it 

yields a ratio of 2.7. These results suggest that the solution becomes domain independent when the 

spanwise width is 2.5 times larger than the mixing layer thickness. 

For the shear layer case investigated herein, it was possible to extend the domain out until 48mm. 

This should ensure uninhibited growth of the turbulent structures, even at the furthest downstream 

measurement location of 450mm. For stations much beyond that, the spanwise width would likely need to 

be extended even more. However, at some point the spanwise domain cannot be extended without taking 

sidewall effects into consideration. The proposed factor of 2.5 cannot necessarily be generalized to other 

more complex flows. For other similar problems where resolution is scarce, a spanwise width that is only 

2.5 times larger than the relevant length scale may be sufficient to prevent artificial suppression of ݓ௠௥௦. 

(A) (B) 

(C) (D) 
NASA/TM—2015-218741 24 

VI. Conclusions 

This work compared a high-order method with a low-order method for performing large-eddy 

simulations (LES). The high-order method (WRLES) employed an eleven-point DRP scheme that is 

formally fourth-order accurate in space. For time stepping, it used the four-stage Runge-Kutta method, 

which is fifth-order accurate locally. For subgrid scale modeling, it uses Vreman’s model. The low-order 

method (Wind-US) is second-order accurate spatially and uses the Roe Physical scheme. For time 

stepping, it used an implicit backward Euler step that is first-order accurate locally, and does not utilize 

any subgrid scale modeling. LES were performed using both methods to capture the turbulent mixing 

between two supersonic streams at a convective Mach number of 0.46. The high-order and low-order 

methods were evaluated at two different levels of grid resolution. In order to compare with the 

experiment, the inflow profiles for LES of the shear layer were obtained using Reynolds Averaged 

Navier-Stokes (RANS) modeling over the upstream splitter plate such that the displacement thickness 

matched the experimental value reported. The results presented in this paper are summarized as follows: 

 For a fine grid resolution, the low-order method produced a similar solution to the high-order method. 

Both methods predict turbulent stresses that are in good agreement with experimental data. 

 Effects of subgrid scale modeling were found to be negligible at both resolutions using the higher 

order method and the Vreman SGS model. At even coarser resolutions, the SGS model might have a 

larger effect. 

 WRLES and WIND-US density contours were compared. For both methods, vortex shedding was 

observed to originate at the splitter tip and quickly transition to turbulence. More organized structures 

were seen with the high-order method whereas those structures appeared to be smeared with the loworder method. 

 The time for all initial transient effects to be eliminated and before statistics can be gathered is no less 

than 5.5 FTTs (Flow Through Times). The recommended sampling duration is 9.5 FTTs or longer. 

Any time-averaging window that meets these two conditions will not produce any substantial 

uncertainty. 

 When the DRP filter was nearly turned off, it was affirmed that its purpose is primarily for numerical 

stability and does not act as a SGS model. Thereby confirming that the reason the low-order and highorder method agreed was because the grid was sufficiently fine enough. 

 For a coarse grid of 3.6 million points, the high-order flow solver arrived at the same solution it 

reached with the 14.4 million point fine grid. However, the low-order solver did not demonstrate grid 

independence between the coarse and fine grid. Thus, high-order algorithms need less resolution than 

low-order methods to achieve the same level of accuracy. 

 We conducted a parametric study for the effect of the spanwise width on the accuracy of the solution 

by varying the width from a 6mm depth to a 48mm depth. As the shed vortices travel downstream and 

grow in size, the spanwise width may not be large enough to resolve all the spanwise modes and may 

introduce artificial ones. The peak root-mean-square of the spanwise fluctuation, ݓ௠௥௦ , was found to 

increase by fifty percent, when the spanwise width was increased from 6mm to 48mm. A spanwise 

width that is at least 2.5 times larger than mixing layer thickness was found to be necessary to capture 

three-dimensionality of vortex shedding. 

 Despite the differences in parallelization and time-stepping scheme, the two method’s computational 

costs were about the same. The larger stencil size and explicit time-stepping of the high-order code 

resulted in a larger computational cost for a given grid size compared to the low-order, implicit timestepping code. However, the input/output performance of the low-order code was slowed by the need to 

regularly checkpoint the solution and spawn an external tool to extract the unsteady flow information. 

 
NASA/TM—2015-218741 25 

The high-order method herein proved to be more accurate on a coarse mesh when compared against 

the low-order method. The small stencils of low-order methods are advantageous for their ability to 

handle complicated geometries whereas high-order methods are constrained by a large stencil size. The 

work herein establishes that low-order methods can be used for LES but require a larger number of grid 

points than their counterpart. Historically, most high-order methods do not retain their accuracy with 

unstructured grids. Recently, unstructured grids are favored for the inherent ability to redistribute mesh 

points where needed and remove unnecessary mesh points. Ideally the benefits of low-order and highorder methods would be combined by creating stable high-order unstructured algorithms. Such a method 

may push the state of the art by enabling high-order unstructured LES of practical applications. 

References 

1. DeBonis J.R., “A High-Resolution Capability for Large-Eddy Simulation of Jet Flows,” 

AIAA Paper 2010-5023, 2010. 

2. DeBonis J.R., “An Examination of the Spatial Resolution Requirements for LES of a Compressible 

Jet,” Quality and Reliability of Large-Eddy Simulations II, Salvetti M.V., Geurts B., Meyers H., and 

Sagaut P., editors, Springer, pp. 329–338, 2010. 

3. DeBonis J.R., “Solutions of the Taylor-Green Vortex Problem Using High-Resolution Explicit 

Finite Difference Methods,” AIAA Paper 2013-0382, 2013. 

4. Nelson C., “An Overview of the NPARC Alliance’s Wind-US Solver,” AIAA Paper 2010-0027, 

2010. 

5. Yoder D.A., “Initial Evaluation of an Algebraic Reynolds Stress Model for Compressible Turbulent 

Shear Flows,” AIAA Paper 2003-0548, 2003. 

6. Yoder D.A., DeBonis J.R., Georgiadis N.J., “Modeling of Turbulent Free Shear Flows,” 

AIAA Paper 2013-2721, NASA/TM—2013-218072, 2013. 

7. Brown G.L., Roshko A., “On Density Effects and Large Structure in Turbulent Mixing Layers,” 

Journal of Fluid Mechanics, Vol. 64, pp. 775–816, 1974. 

8. Chinzei N., Masuya G., Komuro T., Murakami A., Kudou K., “Spreading of Two-Stream 

Supersonic Turbulent Mixing Layers,” Physics of Fluids, Vol. 29, pp. 1345–1347, 1986. 

9. Papamoschou D., Roshko A., “The Compressible Turbulent Shear Layer: An Experimental Study,” 

Journal of Fluid Mechanics, Vol. 197, pp. 453–477, 1988. 

10. Goebel S.G. , Dutton J.C. , “Experimental Study of Turbulent Compressible Mixing Layers,” 

AIAA Journal, Vol. 29, pp. 538–546, 1991. 

11. Samimy M., Reeder M.F., Elliott G.S., “Compressibility Effects on Large Structures in Free Shear 

Flows,” Physics of Fluids, Vol. 4, pp. 1251–1258, 1992. 

12. Hall J.L., Dimotakis P.E., Rosemann H., “Experiments in Nonreacting Compressible Shear Layers,” 

AIAA Journal, Vol. 31, pp. 2247–2254, 1993. 

13. Clemens N.T., Mungal M.G., “Large-Scale Structure and Entrainment in the Supersonic Mixing 

Layer,” Journal of Fluid Mechanics, Vol. 284, pp. 171–216, 1995. 

14. Foysi H., Sarkar S., Friedrich R., “Compressibility Effects and Turbulence Scalings in Supersonic 

Channel Flow,” Journal of Fluid Mechanics, Vol. 509, pp. 207–216, 2004. 

15. Freund J.B., Lele S.K., Moin P., “Compressibility Effects in a Turbulent Annular Mixing Layer, 

Part 1., Turbulence and Growth Rate,” Journal of Fluid Mechanics, Vol. 421, pp. 229–267, 2000. 

16. Pantano C., Sarkar S., “A Study of Compressibility Effects in the High-Speed, Turbulent Shear 

Layer using Direct Simulation,” Journal of Fluid Mechanics, Vol. 451, pp. 329–371, 2002. 

17. Bogey C., Bailly C., “A Family of Low Dispersive and Low Dissipative Explicit Schemes for Flow 

and Noise Computations,” Journal of Computational Physics, Vol. 194, pp. 194–214, 2004. 

18. Berland J., Bogey C., Bailly C., “Low-Dissipation and Low-Dispersion Fourth-Order Runge–Kutta 

Algorithm,” Computers and Fluids, Vol. 35, pp. 1459–1463, 2006. 

19. Georgiadis N.J., Alexander J.I.D., Reshotko E., “Hybrid Reynolds-Averaged Navier–Stokes/LargeEddy Simulations of Supersonic Turbulent Mixing,” AIAA Journal, Vol. 41, pp. 218–229, 2003. 
NASA/TM—2015-218741 26 

20. Mankbadi M.R., Georgiadis N.J., “An Examination of Parameters Affecting Large-Eddy 

Simulations of Flow Past a Square Cylinder,” AIAA Paper 2014-2089, 2014. 

21. Inoue M., Pullin D.I., “Large-Eddy Simulation of the Zero-Pressure-Gradient Turbulent Boundary 

Layer up to Reθ = O (1012),” Journal of Fluid Mechanics, Vol. 686, pp. 507–533, 2001. 

22. Vreman A. W., “An Eddy-Viscosity Subgrid-Scale Model for Turbulent Shear Flow: Algebraic 

Theory and Applications,” Physics of Fluids, Vol. 10, pp. 3670–3681, 2004. 

23. Pope S. B., Turbulent Flows, Cambridge University Press, 2000. 

24. Jarrin N., Prosser R., Uribe J.-C., Benhamadouche S., Laurence D., “Reconstruction of Turbulent 

Fluctuations for Hybrid RANS/LES Simulations using a Synthetic-Eddy Method,” International 

Journal of Heat and Fluid Flow, Vol. 30, pp. 435–442, 2009. 

25. Ongoren A., Rockwell D., “Flow Structure from an Oscillating Cylinder Part 1, Mechanisms of 

Phase Shift and Recovery in the Near Wake,” Journal of Fluid Mechanics, Vol. 191, pp. 197–223, 

1988. 

26. Mankbadi M.R., Balachandar S., “Compressible Inviscid Instability of Rapidly Expanding Spherical 

Material Interfaces,” Physics of Fluids, Vol. 24, 2012.