# User:Oleg Schilling/Proposed/Turbulence: Monotone-Integrated and Implicit Large-Eddy Simulation

Figure 1: Your introductory color figure is here to "advertise" the article or summarize the main point.

The modeling of turbulent flows is immensely important for a large variety of fluid flows of scientific and engineering interest. The technique of large eddy simulation (LES) is a promising avenue to offer highly detailed simulations of turbulence in a time-accurate manner without excessive averaging inherent in less detailed approaches. In recent years, the technique of implicit LES has gotten attention and produced promising results on very challenging high Reynolds number flows. The use of high-resolution numerical methods to provide the effective turbulence model is not without significant controversy. Nonetheless, the impressive results speak for themselves, and analysis has emerged to provide a more systematic explanation for the success of the methods. Below, we discuss the history, use, and techniques used for ILES.

## Introduction

As we address issues of simulation of practical flows at moderate-to-high Re with very complex physics and geometry, Direct numerical simulation (DNS) – resolving all relevant space/time scales – is prohibitively expensive. On the other end of the CFD spectrum are the industrial methods such as the Reynolds-Averaged Navier-Stokes (RANS) approach, which simulate the mean flow and approximately model the effects of turbulent scales. In large eddy simulation (LES), the large energy containing structures are resolved whereas the smaller, presumably more isotropic, structures are filtered out and unresolved SGS effects are modeled [Sagaut]. By necessity, more than choice, LES becomes an effective intermediate approach between DNS and RANS. LES is capable of simulating key unsteady flow features that cannot be handled with RANS and provides higher accuracy. Adding to the physics based difficulties in developing and validating SGS models for LES, are truncation terms due to discretization that are comparable with SGS models in typical LES strategies [Ghosal]. LES resolution requirements thus become prohibitively expensive for practical flows and regimes. Implicit LES (ILES, MILES) [5] (hereafter referred to as IB, or IB-ChN when Chapter N is specifically cited), effectively addresses the seemingly insurmountable issues posed to LES by under-resolution, by relying on the use of SGS modeling and filtering provided implicitly by physics capturing numerics.

ILES adopts the same goals as LES, but follows a considerably different approach in pursuing their achievement. ILES is further distinguished by being applied to a wider and more challenging range of fluid dynamic applications than LES. We will examine how ILES is constructed to achieve the goal of turbulent simulation, and explore its application to a variety of challenging flows.

The chief difference in approach used in ILES compared with LES is the manner in which SGS modeling is pursued. In conventional LES an SGS model is added to the fluid dynamics calculation to account for the unresolved scales of turbulent motion. The SGS model is designed using the detailed knowledge of the structure and behavior of turbulent flows with an eye toward achieving these characteristics in the simulation. On the other hand, ILES relies upon the characteristics of a numerical method, which abides by a set of physical principles to achieve the same end as the explicit SGS modeling used in LES. At first, the approach taken in ILES may seem lazy and prone to error, but empirical experience has shown that ILES can achieve high quality simulations of turbulence. This article will introduce the reader to the basic concepts and the evidence supporting the above conclusion.

### Historical background

The history of LES and SGS models begins at the dawn of modern computational simulation for weather and climate, and connects back to the origin of computational methods for shock capturing, the concept of an artificial viscosity, which allows the wedding of a shock wave to a discrete grid.

Historically, the augmented dissipation was first implemented in Lagrangian simulations of shocked flows, where the equations were explicitly augmented by dissipative terms known as artificial viscosity (Richtmyer 1948; von Neumann and Richtmyer 1950). As we shall recount shortly, this idea was extended to turbulent flows, where the dissipative terms are SGS models. Development of explicit SGS models for turbulence has continued and grown ever more sophisticated without any acknowlegement of their origin. However, the evolution of artificial viscosity took a different turn in the early 1970s, when Jay Boris and Bram van Leer independently introduced the first NFV methods. In these methods, the entropy or physical realizability condition is satisfied implicitly as part of the numerical method. Since that time many new improvements to NFV methods have been developed with this class of methods becoming the norm in simulation codes. NFV methods for shock flows that become the basis of ILES exhibit many advantages over other approaches, including nonlinear stability, computational efficiency, ease of implementation and above all accurate and realistic results (physical realizability is built into these methods). For these reasons, such methods have become the preferred choice for many challenging problems in the field of computational fluid mechanics.

After World War II, John von Neumann worked to expand the role of numerical simulation in science extending efforts made in the development of the first atomic weapons. One of his efforts was in the area of numerical weather prediction, where he worked with Jules Charney at the institute for Advanced Study (IAS) during the early 1950's. In 1950 the first crude simulation of weather was completed at the IAS. By 1956 more advanced simulations were being attempted. Both von Neumann and Charney were present at a conference at IAS where Phillips presented his two-dimensional simulation of a month of weather of the Eastern North America (the Eastern seaboard). Also present at the conference was Joseph Smagorinsky a more junior member of the weather modeling team. It was observed that Phillip's calculation was polluted by oscillations (called “ringing”) late in the simulated month and Charney made the suggestion that von Neumann's viscosity could be used to eliminate the ring oscillations (Smagorinsky 1983). Smagorinsky was given the task of extending Phillip's results to three dimensions, including artificial viscosity to eliminate any unwanted oscillations.

Smagorinsky's implementation (Smagorinsky 1963) resulted in the first LES SGS model, and formed the basis for much future work. In fact, Smahorinsky’s name is given to the cnonical LES GSG model. After the fact, Lilly (1968) made a more rigorous connection of the Smagorinsky eddy viscosity to turbulence theory. Thereafter, SGS modeling has grown and evolved. However, the energy dissipation associated with the original Smagorinsky form persists in many more sophisticated models, e.g., the popular dynamic Smagorinsky models and mixed models. The important fact is that the SGS contributes not only to the modeling of unresolved physics, but also to the numerical stability of the integrations of equations.

Because von Neumann-Richtmyer artificial viscosity is the first shock capturing method, subsequent shock capturing methods owe their existence to its precedence. Nonetheless, the path to NFV methods was a little less direct, and begins with the early work of Peter Lax (Lax 1954; Lax 1973). In particular, the paper by Lax and Wendroff [LW62] first emphasized the importance of conservative methods (cf finite volume methods; flux methods). A later paper [LW64] describes a second-order accurate numerical approximation for the advective terms in which first-order diffusive errors are directly compensated. This Lax-Wendroff scheme produces oscillatory fields behind a shock wave, in contrast to the Lax-Friedrichs method, which produces monotone shock transitions, but is overly diffusive. Both methods differ in an essential way from the artificial viscosity methods in that they are "linear. " Specifically, a linear method uses the same stencil everywhere, whereas artificial viscosity is nonlinear in the sense that its magnitude depends on the flow variables themselves. However, it took nearly twenty more years for the importance of nonlinearity to be recognized in a substantial manner.

Sergei Godunov was a graduate student in the Soviet Union in the early 1950s. As part of his doctoral thesis, he was assigned to calculate shock wave propagation on a new computer at the Keldish Institute (he was studying under Keldish himself). At the time, existing algorithms in the Soviet Union were not sufficiently accurate or stable, and Godunov was not familiar with the work of Lax and Friedrichs (Lax 1954). Instead, he developed a new methodology, based on solving local Riemann problems, that not only satisfied his degree requirements; however, sewed the seeds of a computational revolution nearly twenty years later. Perhaps of greater importance, Godunov proved a fundamental theorem of numerical analysis, which states that no numerical method can be simultaneously linear, second-order accurate, and monotonicity preserving. This is a barrier theorem, which negatively tells one what cannot be done.

Godunov's method and theorem were published in 1959. However it was largely ignored in both the East and the West and laid largely dormant for over a decade. Then in 1971, two scientists, Jay Boris in the United States and Bram van Leer in the Netherlands, overcame the barrier of Godunov's theorem by recognizing the importance of giving up linearity. The nonlinearity simply means that the finite difference stencil used by the method is a function of the solution itself (i.e., a linear method uses the same stencil everywhere all the time). The flux-corrected transport (FCT) method (Boris and Book 1973) focused on eliminating unphysical oscillations. It is a hybrid scheme that mixes first-order and second-order accurate in a nonlinear manner. The MUSCL schemes of van Leer (van Leer 1974) are a more direct generalization of Godunov's method in which the initial states of the Riemann problem are modified by nonlinear limiting of the size of interpolating gradients used to define the approximation. The class of methods defined by Boris and van Leer formed the basis of the methods that would first demonstrate the ability to conduct ILES simulations.

(Aside: As fate would have it, a third scientist, Kolgan, produced an alternate generalization of Godunov's method that also overcame the barrier of Godunov's theorem. In today's parlance, his method would take the label of a second-order essentially nonoscillatory (ENO) method. Unfortunately, this work went unnoticed and Kolgan died before receiving any recognition.)

Boris describes the collision between the LES and ILES worlds in the following passages (from Ch. 1, ILES Book): “Professor John Lumley hosted a meeting at Cornell University in 1989, Whither Turbulence: Turbulence at the Crossroads" (Lumley 1989). The conference was focused on controversial questions in turbulence research that that John felt had not been satisfactorily resolved... As a responding presenter, my paper focused on giving the evidence and reasons for the surprising success of monotone no model LES methods. I named this approach to turbulence simulation Monotone Integrated Large Eddy Simulation (or MILES) as an up-front reminder that monotone methods come with an integrated LES turbulence capability built in… Looking back before the Cornell conference, my plunge into turbulence simulation actually dates to the early 1970’s when Flux-Corrected Transport (FCT) was invented. It soon became clear that FCT was good for much more than dynamic shocked flows but recognizing just how much was going on under the surface in the computed solutions took almost another two decades, until 1989.” By all accounts Boris was the first to begin attempting a detailed explanation for the success of these methods for turbulence, which he termed “a convenient conspiracy” (Boris and Oran 1993). In describing this conspiracy, he provides four key points, (1) turbulence is unsteady and should computed as such, (2) nonlinear non-oscillatory methods have show good agreement with experimental data without explicit SGS models, (3) these methods are well-suited to solving very difficult multi-physics problems and (4) more computational resolution is far more efficient in improving results that improvements is SGS models.

Moreover they observed that such methods impose important physical principles on the solution notably conservation, causality and locality. Conservation is built into the method by construction. Causality and locality result from the basic structure of the explicit numerical method. Beyond these properties, the non-oscillatory nature of the methods keeps solutions physical and provides enough dissipation to assure that solutions can be physically realized. Finally, they reasoned that the very nature of turbulence (the scale dependent transfer of energy) itself also lend themselves to being simulated by non-oscillatory methods.

Simple functional SGS models focus on dissipating energy at a physically correct rate. Newer SGS models focus more on scale-dependent effects in turbulence. The underlying idea of these new approaches has been to represent the effects of the unresolved dynamics by regularizing the larger scales of the flow. Such regularization may be based on physical reasoning resulting from an ab initio scale separation, or on the enforcement of numerical constraints that enforce monotonicity and entropy conditions. The latter is the underlying concept of ILES. As will be explained next, ILES may work simply because turbulence and shock physics scale in almost identical manners. As such the correct solution of shock problems may result in building the same scaling into methods used for turbulent flows.

MEA was first developed to assess the stability of numerical algorithms [hirt1968]. Basically, the analysis is Taylor series expansion applied to the discrete terms of the algorithm as if the equation (and its solution) were continuous. The main assumption is the use of Taylor analysis itself, which implies restrictions on the smoothness of the function. In particular, the Taylor analysis produces an infinite series of terms, but truncating this series, and keeping only the lowest order term, which is assumed to dominate the expansion, forms the modified equation. A discussion of the uses of MEA and caveats of its use can be found in [knoll].

Given the close dependence of ILES simulation upon NFV methods, it is not surprising that many of the pioneers of ILES were the pioneers who had also developed the underlying NFV methods. There is approximately a fifteen year interval before the introduction of the first NFV methods, FCT in 1973 [fct] MUSCL in 1973 [vanleer], and the first publications in ILES.

While it may not be possible to sort out the earliest efforts at simulating turbulent flows with these schemes, it is clear that credit for the first public documentation of the approach belongs to Jay Boris and colleagues [lumley]. Boris raised the idea that a minimal nonlinear subfilter and a suitable SGS reconstruction might be implicitly provided by discretization with a particular class of (monotonicity preserving) convection algorithms, when effectively building crucial SGS physics into the numerics to enforce conservation laws and positivity. He suggested that these built-in features appeared to ensure efficient transfer of the smallest grid-scale motions, generated by computationally resolved fluid dynamic mechanisms, smoothly off the resolved grid. These conjectures led to proposing the MILES approach [FDR], largely based on noting the {\it convenient conspiracy] of physical and computational circumstances that appeared to make such approach practical for high Reynolds number turbulent flows [conspiracy].

Later theoretical studies by Fureby and Grinstein used the modified equation (ME) framework to show that certain NFV (i.e., flux-limiting) algorithms with dissipative leading order terms can provide appropriate built-in (implicit) SGS models of tensorial (generalized) eddy-viscosity type [grin], [fg2002] [also 1999AIAAJ]. Key features in the comparisons with classical LES, leading to the identification of the said MILES implicit SGS model, were the use of the MEA framework, and the fact that the FV formulation readily allowed to recast the leading order truncation terms in divergence form. The MEA and FV frameworks are crucial basic elements in the analysis of ILES on which we elaborate further below. Extensive use of flux-limiting based MILES has been carried out for canonical and complex [REFS1], as well as reacting flows [REFS2].

The ILES work of Paul Woodward, one of the originators of PPM (piecewise parabolic method [ppm]) and colleagues, was applied to astrophysical problems [porter]. This is a regime of highly compressible flow with extremely high Reynolds' numbers. This method has been verified and validated on simple, classical flows, while being applied to very challenging astrophysical flows, which are characterized by extremely large Reynolds numbers. The performance of a PPM-based ILES on ideal turbulent flows provides distinct confidence in its ability to perform well in the applied context.

At about the same time, David Youngs' and colleagues [youngs] were applying van Leer methods to modeling the growth of turbulent regions and the mixing resulting from fluid instabilities, including Raleigh-Taylor and Kevin-Helmholtz. These calculations contain adjacent regions of very high and very low Reynolds' number, illustrating a very useful feature of ILES simulation -- that the same fluid solver can be used for smooth and for turbulent flows.

Margolin, Smolarkiewicz and colleagues published the first applications of ILES to geophysics. In [LES], they showed simulations of the atmospheric boundary layer using MPDATA [mpdata]. As in the astrophysics cases, these calculations involved very high Reynolds' number ($Re \sim 10^6$), but with stratified and nearly incompressible flow. This work compared two simulations using the same fluid solver, one with and one without an explicit SGS model. The correspondence of these two results demonstrates the adaptiveness of the NFVschemes, that the implicit dissipation of the numerical method does not simply add to the explicit dissipation of a turbulence model.

The papers above demonstrate the effectiveness of the ILES approach in a variety of applications, but they do not address the question of why the approach is successful. The first attempt to provide a rationale for implicit turbulence modeling is found in a 2002 paper by Margolin and Rider [rationale], which put forth three of the principal ideas of this paper; first, that a finite-scale equation can be derived by averaging and renormalizing the continuum PDE; second, that the MEA accurately describes the behavior of the NFV scheme; and third, that a comparison of these two equations would provide the justification for ILES. The results in [rationale] were restricted to Burgers' equation in one spatial dimension, and will be extended to the two dimensional Navier-Stokes equations in sections 4and 5 of this paper.

A second significant paper regarding rationale was published in 2003 [aiaa]. This paper had a more engineering perspective, comparing the implicit SGS models resulting from MEA of limiter based algorithms with the results applicable to TVD, FCT, Godunov-type, ENO and other similar algorithms drawing direct comparison to the explicit SGS models in use in the turbulence modeling community, We believe this is an important connection that may facilitate the acceptance of ILES in the broader community, and we will amplify these ideas in later in this report.

Although the early pioneers in ILES worked independently and in fact were mostly unaware of each other's accomplishments with computing turbulent flow, the past few years have seen the reversal of this isolation. At the 2001 Summer ECCOMAS conference, two entire sessions were devoted to VLES (very large eddy simulation), many of the papers addressing issues of implicit turbulence modeling. Special sessions on alternative methods for LES (including ILES) were held in AIAA conferences in 2002 (Reno) and 2003 (Orlando), as well as the February SIAM meeting on computational science and engineering in 2003, 2004 ILES workshops at LANL, and ILES minisymposia at the Summer 2006? ECCOMAS. Selected papers of these conferences appeared in special sections of archival journals (IJNMFD JFE 2002, and 2007 JFE). Recently, the first ILES book was published by Cambridge University Press (2007). This book represents a collaboration of many of the leading researchers in ILES presenting original research contributions covering theory, verification, validation and application of ILES techniques to a variety of practical applications.

High Reynolds' number flow has posed a continuing challenge to the computational fluid dynamics (CFD) community. The fact that the dissipative scales of motion cannot be represented on the computational mesh presents a significant problem, since dissipation is a necessary element of all continuum physical processes. The earliest solution to this problem was published in 1950 [VNR] for the simulation of fluid flows with shocks. In the framework of Lagrangian coordinates, von Neumann and Richtmyer proposed to supplement the physical viscosity by an artificial viscosity. This artificial viscosity was designed to eliminate unphysical oscillations behind the shock while preserving the shock speed and the jump conditions across the shock. An important feature of the artificial viscosity is that the viscous coefficient depends on the solution; this is commonly termed nonlinear viscosity. Artificial viscosity remains an enduring component of Lagrangian algorithms for high speed flows.

Lagrangian frameworks are not suitable for flows with significant shear or vorticity, and the CFD community has preferred Eulerian coordinates for such flows. In Eulerian coordinates (and indeed in all non-Lagrangian coordinates), one must take into account advection -- the exchange of mass and its associated properties between cells. The representation of the advective fluxes in flows with shocks was problematic for more than 20 years after artificial viscosity was introduced. First-order advective fluxes (e.g., donor cell) are too diffusive whereas higher-order advective fluxes (e.g., Lax-Wendroff) produced oscillatory solutions.

A solution to this impasse was forthcoming from two different sources in the early 1970s. Although put forward in somewhat different contexts, both the flux corrected transport methods [fct] and the MUSCL methods [vanleer] had two important similarities. Each was designed to eliminate the unphysical oscillations by construction, and each accomplished this by invoking nonlinear (i.e., solution dependent) approximations to the advective flux. Both methods enjoyed immediate success. In the succeeding years, many more NFV methods have been devised. Most of these methods are written in finite volume form, which guarantees exact conservation. NFV methods have become the mainstream of CFD.

Although our discussion so far has focused on shocks, the issues of adding artificial dissipation extends also to high Reynolds number turbulent flow. The methods developed for shock capturing are also essential in computing flows with under-resolved gradients and discontinuities. This similarity was recognized early by Smagorinsky, who extended the artificial viscosity of von Neumann and Richtmyer to multidimensions to suppress unphysical oscillations in early weather forecasts [smag]. This extension became the first explicit SGS model, now well known as the Smagorinsky model. An interesting account of the origin of this work can be found in [smag2].

An additional important similarity relates to the inherent anisotropy of high-Re turbulent flows in the high-wave-number end of the inertial subrange region, characterized by very thin filaments of intense vorticity and largely irrelevant internal structure, embedded in a background of weak vorticity.[jimenez]. Basing ILES on NFV numerics, can be viewed as requiring a sharp velocity-gradient capturing capability operating at the smallest resolved scales. By focusing on the inviscid inertial-range dynamics and on adaptive regularization of the under-resolved flow, ILES thus follows very naturally on the historical precedent of using this kind of schemes for shock capturing – in the sense that requiring emulation (near the cutoff) of the high wave-number-end features of the inertial subrange region of turbulent flows is analogous to spreading the shock width to the point that it can be resolved by the grid[fg2002].

### Modified equations for ILES

Formal analysis of LES can be based on a procedure called modified equation analysis (MEA) [6]. For the sake of the discussion, we focus here on a conceptually simple turbulent mixing case, that of incompressible flow with scalar mixing – an example of typical additional physics to be simulated with the flow. More detailed discussions can be found elsewhere ([7], IB-Ch5). The modified equations – satisfied by the numerically calculated solutions – are the following,

where the bar denotes spatial filtering, v is the solenoidal velocity field, θ is a conserved material scalar concentration, and denote momentum and material diffusivity, respectively, tv and tθ address effects of discretization and commutation between differentiation and filtering. To ensure closure of the equations in the filtered unknowns, models for and have to be provided; depending on based Reynolds number ((Re) and Schmidt number Sc= / , different SGS models for the scalar and the momenta may be needed in general.

As noted, a crucial practical LES computational aspect is the need to distinctly separate the effects of spatial filtering and SGS reconstruction models from their unavoidable implicit counterparts due to discretization. Indeed, it has been noted that in typical LES strategies, tv and tθ – truncation terms due to discretization and filtering, have contributions directly comparable with those of the explicit models [8]. Seeking to address the seemingly insurmountable issues posed to LES by under-resolution, the possibility of using the SGS modeling and filtering provided implicitly by the numerics has been considered as an option: this is generally denoted as numerical LES by Pope [9]. Arbitrary numerics will not work for LES: good or bad SGS physics can be built into the simulation model depending on the choice and particular implementation of the numerics.

The monotone integrated LES (MILES) approach first proposed by Boris [10], incorporates the effects of the SGS physics on the resolved scales through functional reconstruction of the convective fluxes using locally monotonic FV schemes. The more broadly defined ILES (IB) generally uses high-resolution non-oscillatory FV (NFV) algorithms – as defined by Harten [11] – to solve the unfiltered Euler or Navier-Stokes (NS) equations. Popular physics capturing methods have been used in ILES, such as flux-corrected transport, the piecewise parabolic method, Godunov, hybrid, and total variation diminishing algorithms. By focusing on inertially dominated flow dynamics and regularization of under-resolved flow, ILES follows on the precedent of using NFV methods for shock capturing – requiring weak solutions and satisfaction of an entropy condition.

MEA provides the framework to reverse-engineer desirable physics into the numerics design. MEA was used in the early formal comparisons [12] between MILES and traditional LES, to show that a class of NFV algorithms with dissipative leading order terms provides appropriate built-in (implicit) SGS models of a mixed tensorial (generalized) eddy-viscosity type. The analysis examines specific implementation aspects of the numerics, e.g., specific temporal integration schemes and synchronization issues, and can address how prescribed anisotropies introduced by non-uniform adaptive gridding – affecting any LES – contribute to the implicit SGS stress tensor in the modified equations [7,12]. Key features in [7,12] were the MEA framework and the use of the FV formulation. Volume integrals in the FV representation naturally link with the discrete spatial filtering operation in LES – the so-called top-hat filtering, and their use readily allow recasting leading-order truncation terms in divergence form. MEA of the implicit SGS models shows that FV discretizations are to be preferred over finite differences for ILES (IB–Ch5). Finally, FVs and MEA have been crucial ingredients in the analysis connecting ILES with the solution of finite scale equations for laboratory observables [13,14]. The latter connection provides a rationale for ILES: ILES works because it solves the equations that most accurately represent the dynamics of FVs of fluid, i.e., governing the behavior of measurable physical quantities on the computational cells [IB–Ch2].

## Numerical Methods for ILES

Overview and basics

At first blush, one might be tempted to simply use almost any modern numerical method to accomplish ILES. The working scientist should resist this temptation and observe a small number of guidelines in achieving a successful ILES. At the same time this necessitates that some methods be judged unsuitable for ILES. Below the basic guidelines for successful ILES will be stated and some of the better ILES methods will be introduced albeit briefly.

What make a method applicable for ILES? On the other hand, what make a method inapplicable for such use?

The principle guide for determining applicability for ILES is empirical results. More recently, detailed numerical analysis of methods has provided further understanding to compliment the empirical evidence. This numerical analysis is in the form of modified equation analysis, which can be used for nonlinear equations. The modified equation is usually used to determine the truncation error associated with a method, but can be used to elucidate the effective SGS model arising from a numerical method. Typically, this method provides a differential equation effectively solved by a numerical method as determined by the original differential equation plus higher order differential terms often termed as error, but in the context of ILES properly viewed as SGS modeling. The picture that emerges by considering the empirical evidence with the analysis shows that the successful ILES methods combine two key ideas: the methods use conservation form along side elements of high-order accuracy coupled with nonlinear stability commonly termed non-oscillatory methods. The conservation form produces the basic property of the energy cascade essential to turbulent phenomena, and the non-oscillatory methods produce stable, physically realizable solutions. Importantly, the energy cascade has the correct dimensional form provided by Kolmgorov’s theory of inertial range dissipation. It can be observed that conservation ideas applied in the inertial range is key to the proof of Kolmgorov’s theory of inertial range turbulence. The high-order aspects of successful methods assures that the linear numerical errors do not highly influence or corrupt the energy cascade arising from the conservation form.

Specific Methods

FCT

Flux-Corrected Transport (FCT) is one of the earliest non-oscillatory or high-resolution method [44]. It was developed by Jay Boris in the early 1970’s and further in collaboration with David Book through the mid-1970’s. The FCT method uses a combination of high-order accurate with low-order monotone method through a two-step process and a nonlinear correction, which provides oscillation-free (monotone) results. The standard form of the method uses a fourth-order finite difference method providing the high-resolution aspect of the scheme, and the nonlinear limited correction assures that oscillations are absent from the results. Another, quite popular form, of FCT was introduced by Zalesak in 1979, which is more genuinely multi-dimensional in character.

Not surprisingly, Boris found that FCT could be used to simulate turbulent flow to good effect. He had the practical problem of computing chemically reactive flows at very high Reynolds number not accessible to other methods. Nonetheless, the empirical observation by Boris was that the results were quite good. Boris accumulated more examples where turbulent flows were computed with greater success than could be justified by standard arguments. Boris eventually went public with this observation and dubbed the methodology, MILES for monotonically-integrated LES.

PPM

Paul Woodward was another pioneer in non-oscillatory high-resolution methods. In the late 1970’s Woodward traveled to Europe where he worked with Bram van Leer who developed a different class of method similar to FCT. Woodward extended Van Leer’s method in collaboration with Phil Colella and called the resulting method PPM for piecewise parabolic method [47]. In a sense, the method uses a parabola as the SGS model as the method produces a parabolic distribution in each zone to define the approximation. The parabola is constrained as not to cause new extrema in the solution through making certain that it is monotonically increasing or decreasing. Woodward is an astrophysicist and soon turned the method toward computing astrophysical phenomena. Astrophysical phenomena are almost intrinsically turbulent at very high Reynolds number not amenable to standard methods for simulating turbulence. Moreover, astrophysical flows also contain shocks and other discontinuous phenomena further complicating matters. Nonetheless, the PPM appears to be quite successful in both simulating astrophysics, but also standard classical turbulent test such as isotropic box turbulence. Today PPM remains the standard method for astrophysical fluid dynamics simulations.

MPDATA

Piotr Smolarkiewicz developed MPDATA for weather simulations in 1983 [48]. This method is much different than other ILES methods. MPDATA stands for multi-pass donor advective transport algorithm, and uses donor cell or upwind differencing in multiple passes to achieve a different non-oscillatory property, sign-preservation.

### Verification and validation

Figure 1. PDFs of the explicit and implicit SGS viscosities from the various LES models [12].

Challenges for LES (ILES) of high Reynolds-number ((Re) turbulent flow include emulating: • the dynamics of vortical coherent structures (CS) dominating the entrainment process, • the dynamics of transition and turbulence decay, • the proper physical interactions between resolved and SGS scales near the cutoff • the transient and long-term effects of initial conditions (ICs).

In what follows, we present an overview of research progress in addressing these challenges with ILES. The focus is on representative case studies ranging from canonical to moderately complex flows, including, isotropic turbulence, the Taylor-Green vortex, and inhomogeneous free jets.Extensive ILES applications in engineering, geophysics, meteorology, and astrophysics, are reported in IB, and selected representative examples are reviewed here; other important applications include wall-bounded flows [2,41,42] – where ILES is computationally competitive with the state-of-the-art LES strategies (e.g. [20,46] and IB–Ch16), and shock and acceleration driven turbulent mixing problems, where turbulence is generated by via Richtmyer-Meshkov (RM) and Rayleigh-Taylor instabilities (e.g., IB–Ch13). RM instabilities add the complexity of shock waves and other compressibility effects to the basic physics associated with turbulent mixing; the unique combination of shock and turbulence emulation capabilities makes ILES an effective simulation strategy for these applications – recently surveyed in [43]. In closing, a brief discussion of ILES application in the case of complex frontier flows is given.

We recognize two crucial aspects of the high-Re turbulent flows to be captured: (1) the dominant vortex interaction mechanisms underlying the cascade dynamics in the inertial subrange occur on convective time scales much shorter than the diffusive time scales, and are thus essentially inviscid; (2) when the smallest vortices of turbulent flows – so-called worm vortices (e.g., [15]) – are much thinner than the main flow scales, the details of their internal structure (and actual core diameters) may no longer be significant, and the strengths and positions of the centerlines of such characteristic regions of intense vorticity may be the most important features to be captured. These observations have suggested that nearly inviscid methods capable of handling vortices in a manner similar to that of shock capturing schemes might provide an efficient computational framework in such flow regimes. This has been generally the basis of ILES as well as the premise of the closely related contour dynamics [16] and vorticity confinement methods [17]. For regimes driven by large scale flow features – for which LES is designed, implicit models associated with locally adaptive dynamic NFV methods are capable by themselves (σ≡0) to effectively emulate the SGS physics effects on the statistics of turbulent velocity fluctuations (IB): • the small-scale anisotropies introduced by worm vortices and shocks, • the viscosity independent dissipation characteristic of the turbulent cascade, • the inherently discrete dynamics of finite scale laboratory observables.

Depending on the flow regimes, ICs, and grid resolution, additional modeling may be needed to further address SGS driven flow physics – e.g., near walls, mixing, chemical reaction, and when simulating backscatter. Robustness of the predictions is an important (frequently) unsettled issue. For example, if the IC information contained in the filtered-out smaller and SGS spatial scales can significantly alter the evolution of the larger scales of motion and practical integral measures, then the utility of any LES approach for their prediction as currently posed is dubious. A major research focus is on evaluating the extent to which particular SGS physical effects can be modeled as the turbulent velocity fluctuations, recognizing when additional explicit models and/or numerical treatments are needed and addressing how to ensure that mixed explicit and implicit SGS models effectively act in collaborative rather than interfering fashion.

Isotropic Turbulence The vorticity field in turbulent flows tends to be highly organized with intense vorticity concentrated in elongated filaments. These are the so-called worm vortices, characterizing the smallest coherent structures of turbulent flows,with typical cross-sectional diameters of the order of the Kolmogorovscale[15].The existence of worm vortices can be traced to an inherently anisotropic feature of the small-scale organization of turbulent flows:the fact that high-magnitude vorticity is preferentially aligned with the eigenvector corresponding to the intermediate(weakest)eigenvalue of the rate of strain tensor, while there is very little such preferential alignment for the lower-magnitude vorticity;this is regarded to be a kinematical property independent of the particular dynamical mechanism involved in the vorticity generation [18].Characteristic spectral content and probability density functions (PDFs) associated with the small-scale features of isotropic turbulence provide a useful basis for verification and validation of the various LES (ILES) strategies.

Early positive evaluations of ILES of forced and decaying isotropic turbulence by Porter et al. (IB–Ch7, [19]) and Fureby & Grinstein [12,20] were based on using locally monotonic NFV methods. Comparisons of instantaneous PDFs of explicit (LES) and implicit (ILES) SGS viscosities showed similar qualitative behaviors sensitive to the actual SGS models involved (Fig.1), turbulent spectra (Fig.2) and PDFs of the vorticity (Fig.3) and strain-rate magnitudes from MILES and LES were in good agreement with those of DNS of isotropic turbulence. More recently, further successful ILES assessments in this case were based on a wide range of NFV methods [21-23]; PDFs of velocity derivatives (Fig.4) were benchmarked with the JHU turbulence database ([21], IB–Ch5), and Domaradzki et al. [22] and Thornber et al. [23] reported well-behaved associated spectral eddy viscosities in agreement with theory. In contrast, poor performances in this fundamental context were reported by Garnier et al. [24] using finite difference discretizations of popular shock capturing schemes; the latter failings are attributed (IB–Ch5) to the inherent inadequacy of such implementations for ILES.

Figure 3. CDFs of vorticity magnitude for the various LES models [12] overlaid with DNS data [15]: a) ReT = 96 and b) ReT = 305.

Figure 4. The PDFs of the computed and experimental velocity increments are compared in this figure. The JHU experimentally measured PDFs are quite well reproduced by the ILES methods. This includes the magnitude of the velocity increments in the tails of the PDF (from [21]).

Transition and Turbulence Decay in the Taylor Green Vortex The Taylor Green Vortex (TGV) system is a well-defined flow that has been used as prototype for vortex stretching, instability and production of small-scale eddies to examine the dynamics of transition to turbulence based on DNS [25,26]. As such, it can also be effectively used as a convenient case to test the capabilities of explicit and implicit SGS modeling in simulating the basic empirical laws of turbulence (e.g., [27]), namely the existence of an inertial sub-range on the energy spectra for sufficiently high Re number and the finite (viscosity-independent) energy dissipation limit law. The TVG case has been previously used to assess presumably unwanted effects of numerical dissipation (e.g., Shu et al. [28]), or alternatively, to demonstrate how convective numerical diffusion effects of certain algorithms can be effectively used by themselves to emulate the dominant SGS physics in the high-Re applications [29]. The TGV configuration involves triple-periodic boundary conditions enforced on a cubical domain with box side length 2π cm. The flow is initialized with the solenoidal velocity components,

u = uosin(x)cos(y) cos(z), v = –uocos(x)sin(y)cos(z), w = 0,

the pressure is a solution of the Poisson equation for the above given velocity field, i.e.,

p=po+[ρ(uo)2/8][1+cos(2z)][cos(2x)+cos(2y)],

an ideal gas equation of state for air, and incompressible or low-subsonic compressible conditions have been typically considered. The mathematical flow simulation model used here is based on the conservation and balance equations of mass, momentum, and energy for the flow conditions in [29]. By design, ILES emulates the dynamics of convectively dominated flows characterized by high (but finite) Re ultimately determined by the non-vanishing residual dissipation of the numerical algorithms. ILES models tested on the TGV case used 643, 1283, or 2563 evenly spaced computational cells, and examined nominally inviscid flow (Euler based) or a linear viscous flow for which SGS effects are neglected (NS based). Among the various ILES and LES approaches tested in [29], we revisit here the ILES results based on using flux-corrected transport and Lagrange-Remap NFV algorithms (IB–Ch4a and IB–4c). Figure 5 shows the TGV flow dynamics based on instantaneous visualizations from the present TGV simulations. The figure shows the initial TGV at t*=0, depicts the later transition to increasingly smaller-scale (but organized) vortices (top row), and then to the fully developed (disorganized) decaying worm-vortex dominated flow regime (bottom row), as characteristic of developed turbulence. The particular results shown here were generated with ILES using a 4th order 3D FCT on the 1283 grid. The snapshots are based on (ray tracing) volume renderings of λ2 – the second-largest eigenvalue of the velocity gradient tensor [30]. The evolution in time of the kinetic energy dissipation , where and , denotes mean (volumetric average), has been demonstrated [29] using various ILES and LES strategies to be capable of robustly predicting the time and height of the dissipation peak for the large-Re limit suggested by DNS (Fig. 6). Studies of sensitivity of ILES to grid resolution carried out on the 643, 1283, and 2563 cases, are exemplified in Fig.7. Actual values of Re characterizing the flow at the smallest resolved scales (e.g., based on the Taylor microscale) are not a priori available in LES or ILES. The present predictions for the TGV dissipation rate of the kinetic energy further indicate that (explicit or implicit) SGS models act very similarly in predicting some sort of Re-independent regime asymptotically attained with increasing grid resolution.

Re dependent effects are clearly suggested in Fig. 7 associated with the lower predicted characteristic t* at dissipation peaks (as well as the wider peaks) with the coarser-grid simulations, a trend that is consistently exhibited by the DNS results as Re is lowered (Fig.7). Moreover, we also find noticeable correlation between profiles of mean kinetic energy dissipation rates and mean resolved enstrophy [29], from various different ILES, LES, and grids, with actual peak values increasing with grid resolution. This observed correlation is qualitatively consistent with the relation for an incompressible NS fluid with Reynolds number Re. Thus, the relevance of using a characteristic Re that effectively increases with grid resolution (for both ILES and conventional LES) is suggested by our analysis

Figure 5. MILES of the Taylor-Green vortex from [29]. Evolution of the flow visualizations (right) ranging from the initial at t*=0, to transition (t*~9) to increasingly organized smaller-scale vortices (top row), and then to fully developed (disorganized) worm-vortex flow characteristic of turbulence (t*>30).

Figure 6. Taylor-Green Vortex. Temporal evolution of the kinetic energy dissipation rate –dK/dt from the DNS by Brachet et al. [25,26].

Figure 7. Taylor-Green Vortex. Temporal evolution of the kinetic energy dissipation rate –dK*/dt* ; predictions from 643, 1283, and 2563 3D-FCT [29], and DNS data [25,26].

Inhomogeneous Free (Jet) Flows Large-scale CS dominate entrainment and mixing processes in free shear flows [33,34]. When properly designed, LES approaches are ideally suited to investigate these flows. Major aspects to be captured in the transition to turbulence from laminar ICs, are the dynamics of global instabilities and the complex three-dimensional vortex geometries (see IB–Ch8). Key computational requirements in this context include the ability to isolate the generation and propagation of acoustic disturbances associated with the CS dynamics, which demands adequately resolving the very low fluctuation levels typically involved. The discussion that follows exemplifies aspects of the development of the Kelvin-Helmholtz instability in free jet flows captured with ILES. We comment on the transitional vortex dynamics, and on the trends of the population of the small-scales vortices in the downstream portion of the simulated jets.

Popular passive jet entrainment control approaches are based on geometrical modifications of the jet nozzle that directly alter the flow development downstream; research on flow control with non-circular jets was reviewed by Gutmark & Grinstein [35]. Rectangular jet ILES studies [36,37] with aspect ratio AR = 1-4, were carried out in various non-reactive and reactive regimes. The jet simulations demonstrated the dynamical vorticity geometries characterizing the near field (Figure 12), involving: i) self-deforming vortex rings, ii) interacting ring and braid (rib) vortices, and, iii) a more disorganized flow regime in the far jet downstream, where the rotational fluid volume is occupied by a relatively weak vorticity background with strong, slender tube-like filament vortices filling a small fraction of the domain. Other distinct geometrical features are associated with the occurrence of vortex ring splittings due to reconnection, which are observed the larger values of AR. Underlying aspects of the vortex-ring bifurcation process – namely, reconnection, bridging, and threading– had been conjectured [38] but could not be captured in the laboratory studies – until more recently [39] – and were first demonstrated with ILES [40]. Various vortex interaction mechanisms contribute in a crucial way to breaking of large-scale CS into increasingly smaller ones in fully developed jets; this includes further complex reconnection and bridging mechanisms. Likely cascade scenarios were demonstrated [37,40], suggesting pathways for transition to turbulence based on successive vortex self-inductions, stretching, and reconnections. ILES has proven to be a promising capable simulation tool in this context.

Figure 13 from ILES-EU studies of square jets [36,37], is used to address issues of grid resolution and convergence. To facilitate the analysis of the results, the jet flow was organized using weak forcing with a (single) Strouhal frequency St=fDe/Uj=0.55, chosen within the range of observed laboratory preferred (large scale) jet frequencies, where De is the jet circular-equivalent diameter, and Uj is the jet exit velocity. The figure compares unfiltered data at selected representative times from simulations carried out with identical initial and boundary conditions in grids denoted F and C, with nominal smallest computational cell spacings Δ and 2Δ. The simulations were carried out with fixed Courant number 0.4, and estimated Re (based on Uj, De, and effective algorithm viscosity), were Re>220,00, and Re>78,000, respectively. The smallest grid size used was Δ=D/42, where D is the length of the minor side of the rectangular nozzle. These comparisons and further analysis have been used to examine dependence on spatial resolution (effective Re) of the simulations as well as issues of ILES convergence. The figure shows quite good visual agreement on the large-scale dynamics of ring and hairpin vortices near the jet exit, but their faster coherence breakdown downstream is apparent with increasing resolution (effective Re), as finer dynamical features can contribute and affect the larger scales. Quantitative statements can be made on the trends of the population of the small-scales vortices in the downstream portion of the simulated jets, where the flow is characterized by thinner filament vortices similar to those observed in fully developed turbulent flows. The goal is to quantify characteristic features of the “turbulence” to which the simulated flow transitions from laminar ICs to what is regarded as its established metrics.

Spectral analysis concentrating on the latter portion of the simulated developed jets can be used to investigate the small-scale jet behavior captured by the simulations. Analysis is based on the instantaneous velocity databases for the grid resolutions considered on (uniformly spaced) downstream subvolumes. Figure 14 shows time-averaged plots of the turbulent kinetic energy spectra in the case AR=1. In each case, spatial FFT analysis is based on datasets for 40 successive times separated by a time interval 0.1/f, on the downstream grid sub-volumes for the two resolutions. The largest wavenumber for which spectral amplitudes are plotted corresponds to a wavelength of four computational cells; the spectra depict short simulated inertial subranges consistent with the k–5/3 inviscid subrange of Kolmogorov’s K41 theory, and a longer inertial range for higher resolution case. The simulated inertial range is followed by the faster decay of the amplitudes due to the (fourth-order) FCT algorithm dissipation for wavelengths λ<10D. This limiting length-scale corresponds to approximately twice the smallest characteristic (full-width) cross-sectional length scales of the elongated vortex tubes in the transitional region of the jet [36].

The distribution of vorticity intensities was examined based on the cumulative distribution function (CDF, Fig.15), defined as the volume fraction occupied by vorticity magnitude values above a given threshold level [12,19]. Comparisons with isotropic turbulence DNS data by Jiménez et al. [15] show good agreement of the grid-F and DNS data even for fractional volumes of less than 1%, while differences between DNS and the grid-C data are significant much before, for fractional volumes between 10% and 1%. The structure of the resolved vorticity implies that high intensity regions tend to be organized mostly in elongated tubes. Because “worm” vortices typically involve fractional volumes of order 1% or less [15], of these two cases, only the grid-F simulation resolves these vortices, and this is supported by the volume visualizations of Figure 13. The possible dependence of resolved small-scale jet features on whether ILES was based on the Euler or NS equations were also assessed.

Figure 11. Comparative instantaneous visualizations based on volume visualizations of λ2, the second-largest eigenvalue of tensor S2 + Ω2, where S and Ω are the symmetric and antisymmetric components of the velocity gradient tensor, are shown on the left. Corresponding visualizations of the streamwise vorticity are shown on the right.

Figure 12. Comparative instantaneous volume visualization of the vorticity magnitude based on the database of square jets at the same time simulations on the finest grid left, and coarsest grid right [37].

Figure 13. Compressible (subsonic) jet flow: Time averaged plots of the turbulent kinetic energy spectra in the square jet case [37].

Figure 14. Cumulative distribution function of the vorticity magnitude [37].

FRONTIER FLOWS

For many extremely complex flows of interest in geophysics, astrophysics, and engineering, whole-domain scalable laboratory studies are impossible or very difficult. Moreover, deterministic simulation studies can be very expensive and are critically constrained by difficulties in modeling and validating all the relevant physical sub-processes, and acquiring all the necessary and relevant SPG information. Developing some sort of predictability is the challenge of unavoidably under-resolved simulations of these flows – and ILES is the most likely practical option in this context.

ILES studies of geophysics (IB-Ch14) have been used to address the difficulties of modeling the dynamics of the global ocean and atmosphere (i.e., geostrophic turbulence); the latter are compounded by the broad range of significant length scales (Re ~108) and by the relative smallness of the vertical height of these boundary layers in comparison to their horizontal extent, which accentuates the importance of the backscatter of energy to the larger scales of motion. ILES modeling based on high-order FV upwinding was demonstrated to be capable reproducing the complex features of global climate of the atmosphere and ocean. ILES studies of astrophysics based on using FV piecewise parabolic methods (IB–Ch15) examined homogeneous decaying and driven compressible turbulence, local area models of stellar convection, global models of red giant stars, and the Richtmyer-Meshkov mixing layer.

ILES studies of complex engineering turbulent flows (IB–Ch16), discussed the application of ILES to a variety of applications ranging from incompressible external flows around typical naval applications to external and internal supersonic flows in aerospace applications; cases studied included flows such as (i) the flow around a model scale submarine; (ii) the flow and wave pattern around a modern surface combatant with transom stern; and compressible flows such as (iii) weapon bay flows, and (vi) solid rocket motor flows.

Urban Flow and Dispersion

Effective use of ILES has been reported for large scale urban simulations of flow and dispersal for consequences management [51,52]. Hazardous chemical, biological, or radioactive (CBR) releases in urban environments may occur (intentionally or accidentally) during urban warfare or as part of terrorist attacks on military bases or other facilities. The associated contaminant dispersion is complex and semi-chaotic. Urban predictive simulation capabilities can have direct impact in many threat-reduction areas of interest, including, urban sensor placement and threat analysis, contaminant transport (CT) effects on surrounding civilian population (dosages, evacuation, shelter-in-place), assessing strategies, education and training of rescue teams and services.

Crucial technical issues include both, turbulent fluid transport and boundary condition modeling, and post-processing of the simulation results for consequences management. Relevant fluid dynamic processes to be simulated include, detailed energetic and contaminant sources, complex building vortex shedding, and flows in recirculation zones. Because of the short time spans and large air volumes involved, modeling a pollutant as well mixed globally is typically not appropriate. In typical urban scenarios, both particulate and gaseous contaminants behave similarly insofar as transport and dispersion are concerned, so that the contaminant spread can usually be simulated effectively based on appropriate pollutant tracers with suitable sources and sinks. In some cases, the full details of multi-group particle distributions are required. In such cases, additional physics to be modeled includes particulate fall-out, as well as deposition, re-suspension and evaporation of contaminants. Other crucial issues in the urban CT simulation process include, modeling building damage effects due to eventual blasts, addressing appropriate regional and atmospheric data reduction, and, feeding practical output of the complex combined simulation process into "urbanized" fast-response models.

The simulation model must also incorporate a consistent stratified urban boundary layer with realistic wind fluctuations. The upwind atmospheric boundary layer characterization directly affects the inflow condition prescription required in urban scenario simulations. Wind fluctuation specifics are major factors in determining urban contaminant transport. The flow data from actual field trials or wind-tunnel experiments are typically inadequate and/or insufficient to fully characterize the boundary layer conditions required in the urban flow simulation model. The sensitivity of turbulent flows to particular IC choices is now well recognized [50]. Far field portions of turbulent flows remember their particular near-field features, and the mechanism by which the transition from ICs to particular associated asymptotic flow occurs involves unsteady large-scale coherent-structure dynamics – which can be captured by LES but not by single-point closure turbulence modeling (e.g., RANS). As a very particular consequence, starting with the typical availability of single-point statistical data, there is no unique way to reconstruct a 3D unsteady velocity field with turbulent eddies to define realistic inflow BCs; such data is typically insufficient to parametrize turbulent inflow BCs for LES of inhomogeneous flows. In the recent simulations of flow and dispersal over an urban model using cube arrangements in wind tunnels [51], the available datasets from the laboratory experiments consisted of high quality, spatially dense (but not time-resolved) single-point statistical data. The ILES urban CT model used was the FAST3D-CT code, which involves a scalable, low dissipation, 4th order phase-accurate FCT convection algorithm, implementing direction-split convection, 2nd-order predictor-corrector temporal integration, and time-step splitting techniques (see [51], IB–Ch17, and references therein).

Relevant insights followed from [56], when comparing predicted and measured volume fractions of a tracer scalar in the first few urban model canyons using various different inflow condition models and grid resolutions (Figs. 9ab). The resolutions considered were on the fairly coarse side, e.g., if simulations of flow over a single (surface-mounted) such cube were to be performed. On the other hand, these are resolutions representative of what we can afford to resolve practically in urban simulations relative to typical building dimensions. A comparison of the average tracer concentration profiles from simulations and experiments is shown at selected stations located in the first three canyons. In all cases, the agreement is within a factor of two or better, with agreement somewhat worse in the first canyon – perhaps reflecting questions in resolving the precise details of the release there. Agreement gets better as we move downstream. This also reflects on the simulated and measured mean velocity and fluctuations agreeing better as we move downstream. Prescribing some reasonable inflow turbulence – as opposed to prescribing steady inflow, was found to be critical (Fig. 9cd). On the other hand, the fluid dynamics within the cube arrangement, i.e., beyond the first canyon, becomes somewhat insulated from flow events in the boundaries, i.e., it is less dependent on the precise details of the modeled inflow turbulence and largely driven by the urban geometry specifics within the urban arrangement. In the more recent work [52], a framework is being developed for detailed dispersal predictions in urban and regional settings based on effective linkage of strong motion (and other) codes capable of simulating detailed energetic and contaminant sources associated with the effects of a conventional or nuclear explosion, and ILES capable of emulating the CT due to wind and turbulence fields in the built-up areas.

Figure 22. MILES studies of flow and dispersion over an urban (cube arrangement) model, from [56]. Colors correspond to predicted and measured dispersal of tracer volume fractions in first few urban canyons indicated on the left.

### Subsection Outlook

Accurate predictions with quantifiable uncertainty are essential to many practical turbulent flow applications exhibiting extreme geometrical complexity and broad ranges of length and time scales. Under-resolved computer studies are unavoidable in such applications, and LES (and more typically ILES) becomes the effective simulation strategy mostly by necessity rather than by choice. SGS modeling issues have motivated intense research in the last four decades. A crucial practical computational issue is the need to distinctly separate the effects of explicit filtering and SGS reconstruction models from the unavoidable implicit ones due to discretization. ILES addresses the seemingly insurmountable issues posed to LES by under-resolution, by focusing on the use of SGS modeling and filtering provided implicitly by physics capturing numerics. ILES has been successfully applied to broad range of free and wall-bounded flows, ranging from canonical benchmark flows to extremely complex flows at the frontiers of current flow simulation capabilities. The performance of representative LES and ILES approaches is found to be equally good in emulating physical statistics of turbulent velocity fluctuations, and there is no discriminating characteristic favoring one or another. Looking towards practical complex flows and regimes, however, the ability of ILES to offer a simpler computational environment should be clearly emphasized.

Depending on grid resolution and flow regimes involved, additional explicit SGS modeling may be needed with ILES to address SGS driven flow physics – e.g., near walls, and when simulating backscatter, material mixing, or chemical reaction. A major research focus is on evaluating the extent to which particular SGS physical effects can be implicitly modeled as the turbulent velocity fluctuations, recognizing when additional explicit models and/or numerical treatments are needed, and when so, addressing how to ensure that the mixed explicit and implicit SGS models effectively act in collaborative rather than interfering fashion. An important challenge in this context is that of improving MEA - the mathematical and physical framework for its analysis and development, further understanding the connections between implicit SGS model and numerical scheme, and addressing how to reverse-engineer desirable SGS physics into the numerics.

SGS and SPG modeling issues become intertwined as modeling difficulties due to under-resolution become compounded with the inherent sensitivity of turbulent flows to choice of ICs and boundary conditions. For example, idealized fundamental problems, such as the decay of turbulence simulated in a box domain with mathematically well-defined periodic boundary conditions, may appear to avoid confronting questions of initial and BCs. However, periodic box simulations actually involve all the difficult SGS and SPG issues. On the one hand, it is now well recognized that turbulent flows remember their ICs; on the other, for sufficiently long simulation times, the integral scale of turbulence will eventually saturate since the larger simulated scales cannot have unaffected growth beyond the box size, and this will eventually distort the characteristic power law of the turbulence decay (e.g., [29,32]).

Since nature controls the flow physics independently of SGS and SPG constraints in the laboratory or numerical experimental process, a legitimate question to ask relates to whether instances exist for some sort of convergence of the larger scale observed flow features and dynamics, i.e., of scales much larger than characteristic numerical or instrumental resolution cutoff scales, but presumably small enough spatially and temporally, so that we can assume that they have not been heavily affected by SPG features. Convergence issues vs. resolution are typically problem dependent and very difficult to address in general.

Actual values of Re characterizing the flow at the smallest resolved scales – e.g., based on the Taylor microscale, are not a priori available in LES or ILES. Solutions associated with different resolutions are associated with correspondingly different values of some characteristic effective Re; this is clearly suggested by the comparisons between ILES and DNS results in Fig.7. This is actually inherent to any LES or ILES approach, where the smallest characteristic resolved scale is determined by the resolution cutoff wavelength prescribed by the explicit or implicit spatial-filtering process. These issues need to get directly projected into the process of establishing suitable procedures and metrics for LES verification and validation.

Sandia National Laboratories is a multi-program laboratory operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin company, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. Los Alamos National Laboratory is operated by the Los Alamos National Security, LLC for the U.S. Department of Energy NNSA under Contract No. DE-AC52-06NA25396.

## References

1. F.F. Grinstein 2009, On integrating large eddy simulation and laboratory turbulent flow experiments, Phil. Trans. R. Soc. A, 367, no. 1899 2931-2945.

2. Sagaut P. 2006, Large Eddy Simulation for Incompressible Flows, 3rd Ed., Springer, NY.

3. von Neumann, J. and Richtmyer, R.D., 1950, A method for the numerical calculation of hydrodynamic shocks. Journal of Applied Physics, 21, 232–237.

4. Smagorinsky, J. 1983, The beginnings of numerical weather prediction and general circulation modeling: early recollections. Advancements in Geophysics, 25, pp. 3–37.

5. Grinstein, F.F, Margolin, L.G., Rider, W.J., Editors, Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics, Cambridge University Press, 2nd printing, 2010.

6. Hirt C.W. 1968, Heuristic Stability Theory for Finite Difference Equations, J. Comput. Phys., 2, p 339.

7. Grinstein F.F. & . Fureby C. 2007, On Flux-Limiting-Based Implicit Large Eddy Simulation, J. Fluids Engineering, 129, 1483-1492.

8. Ghosal, S. 1996, An analysis of numerical errors in large-eddy simulations of turbulence, J. Comp Phys., 125, 187-206.

9. Pope, S.B. 2004, “Ten Questions Concerning the Large Eddy Simulation of Turbulent Flows”, New J. Phys., 6, 35.

10. Boris J.P. 1990, “On Large Eddy Simulation Using Subgrid Turbulence Models”, in Whither Turbulence? Turbulence At The Crossroads, J.L. Lumley Ed., Springer, New York, P 344.

11. Harten, A. 1983 High resolution schemes for hyperbolic conservation laws. J. Comput. Phys. 49, 357–393.

12. Fureby C. & Grinstein F.F., 1999, “Monotonically Integrated Large Eddy Simulation of Free Shear Flows”, AIAA Journal, 37, 544-56.

13. Margolin, L.G. & Rider, W. J., 2002: ”A rationale for implicit turbulence modeling,”Int. J. Num. Meth. Fluids, 39, 821–841.

14. L.G. Margolin 2009, Finite-scale equations for compressible fluid flow. Phil. Trans. R. Soc. A, 367, no. 1899 2861-2871.

15. Jimenez, J., Wray, A., Saffman, P., & Rogallo, R. 1993 The structure of intense vorticity in isotropic turbulence. J. Fluid Mech. 255, 65–90.

16. Oberman, E. A., & Zabusky, N. J. 1982 Evolution and merger of isolated vortex structures. Phys. Fluids 25, 1297–1305.

17. Steinhoff, J., & Underhill, D. 1994 Modification of the Euler equations for ‘Vorticity Confinement’: Application to the computation of interacting vortex rings. Phys. Fluids 6, 2738–2744.

18. Jimenez, J. 1992 Kinematic alignment effects in turbulent flows. Phys. Fluids, A4, 652–654.

19. Porter, D. H., Woodward, P. R., & Pouquet, A. 1998 Inertial range structures in decaying turbulent flows. Phys. of Fluids 10, 237–245.

20. Fureby C. & Grinstein F.F., 2002, Large Eddy Simulation of High Reynolds-Number Free and Wall Bounded Flows, J. Comput. Phys., 181, p 68.

21. Margolin, L. G., Rider W. J., & Grinstein, F. F. 2006 Modeling turbulent flow with implicit LES. Journal of Turbulence, 7 (15), 1–27.

22. Domaradzki, J. A., Z. Xiao, and P. Smolarkiewicz: 2003, ‘Effective eddy viscosities in implicit large eddy simulations of turbulent flows’. Phys. Fluids, 15, 3890–3893.

23. B. Thornber, A. Mosedale, and D. Drikakis, 2007, On the implicit large eddy simulations of homogeneous decaying turbulence, Journal of Computational Physics, 226, 1902-1929

24. Garnier, E., Mossi, M., Sagaut, P., Comte, P. & Deville, M., 1999 On the Use of Shock-Capturing Schemes for Large Eddy Simulation, J. Comput. Phys 153, 273-311.

25. Brachet, M.E., Meiron, D.I., Orszag, S.A., Nickel, B.G., Morg, R.H. & Frisch, U.J. 1983 Small-Scale Structure of the Taylor-Green Vortex, J. Fluid Mech. 130, 411.

26. Brachet, M.E. 1991 Direct Numerical Simulation of Three-Dimensional Turbulence in the Taylor-Green Vortex, Fluid Dynamics Research, 8, 1.

27. Frisch, U. 1995 Turbulence. Cambridge University Press Cambridge, p. 57.

28. Shu, C-W., Don, W-S., Gottlieb,D. Schilling, O. & Jameson,L. 2005, Numerical Convergence Sudy of Nearly Incompressible, Inviscid Taylor–Green Vortex Flow, Journal of Scientific Computing, 24, pp. 1-27.

29. Drikakis, D., Fureby, C., Grinstein, F.F. & Youngs, D. 2007, Simulation of Transition and Turbulence Decay in the Taylor-Green Vortex, Journal of Turbulence, 8, 020, 1-12.

30. Jeong, J., & Hussain, F. 1995 On the identification of a vortex. J. Fluid Mech. 285, 69–94.

31. Lesieur, M. and Ossia, S., 2000, 3D isotropic at very high Reynolds numbers: EDQNM study. Journal of Turbulence, 1, 007, 1–25.

32. Skrbek, L. and Stalp, S.R., 2000, On the decay of homogeneous isotropic turbulence. Physics of Fluids, 12, 370, 1997–2019.

33. Brown, G., & Roshko, A. 1974 On density effects and large structure in turbulent mixing layers. J. Fluid Mech. 64, 775–816.

34. Hussain, A. K. M. F. 1986 Coherent structures and turbulence. J. Fluid Mech., 173, 303–356.

35. Gutmark, E. J., & Grinstein, F. F. 1999 Flow control with noncircular jets. Annu. Rev. Fluid Mech. 31, 239–272.

36. Grinstein, F. F., & DeVore, C. R. 1996, Dynamics of coherent structures and transition to turbulence in free square jets. Phys. Fluids 8, 1237–1251.

37. Grinstein, F. F. 2001 Vortex dynamics and entrainment in regular free jets. J. Fluid Mech. 437, 69–101.

38. Hussain, F., & Husain, H.S. 1989 Elliptic jets. Part I. Characteristics of unexcited and excited Jets. J. Fluid Mech. 208, 257–320.

39. Leweke, T., & Williamson C. H. K. 1998 reconnection of a counterrotating vortex pair, in Advances in Turbulence VII, ed. U. Frisch. Dordrecht: Kluwer, 55–58.

40. Grinstein, F. F. 1995 Self-induced vortex ring dynamics in subsonic rectangular jets. Phys. Fluids 7, 2519–2521.

41. Fureby, C., Alin, N., Wikstr¨om, N., Menon, S., Persson, L., & Svanstedt, N. 2004 On large eddy simulations of high Re-number wall bounded flows. AIAA Journal 42, 457–468.

42. Fureby, C., 2007, ILES and LES of Complex Engineering Turbulent Flows, AIAA Journal 129, 1514

43. F.F. Grinstein, A.A. Gowardhan, J.R. Ristorcelli and A.J. Wachtor, “On coarse-grained simulations of turbulent material mixing”, Physica Scripta, 86 (2012) 058203.

44. Kuzmin, D.; Löhner, R.; Turek, S. (Eds.), Flux-Corrected Transport, Springer, 2012.

45. F.F. Grinstein, A.A. Gowardhan, and A.J. Wachtor, "Simulations of Richtmyer-Meshkov Instabilities in Planar Shocktube Experiments", Physics of Fluids, 23, 034106, 2011.

46. Grinstein, F. F., & Fureby, C. 2002 Recent progress on MILES for high Reynolds-number flows. J. Fluids Eng. 124, 848–886.

47. Porter D H, Woodward P R and Pouquet A 1998 Inertial range structures in decaying turbulent flows Phys. Fluids 10 237–45

48. P.K. Smolarkiewicz, and L.G. Margolin, MPDATA: A Finite-Difference Solver for Geophysical Flows, J. of Computational Physics, 140, 459–480 (1998).

49. P.E. Dimotakis, The mixing transition in turbulent flows, J. Fluid Mech. 409, 69, 2000.

50. George, W.K. and Davidson 2004, L., “Role of Initial Conditions in Establishing Asymptotic Flow Behavior”, AIAA Journal, 42, 438-446.

51. Patnaik, G., Boris, J.P., Young, T.R., Grinstein, F.F., F.F. 2007, Large Scale Urban Contaminant Transport Simulations with MILES, J. Fluids Eng., 129, pp. 1524-1532.

52. F.F. Grinstein, R. Bos, and T. Dey, "LES based Urban Dispersal Predictions for Consequence Management", by invitation, ERCOFTAC Bulletin, 78, March 2009, pp. 11-14.