## Abstract

We review existence and non-uniqueness results for the Euler equation of fluid flow. These results are placed in the context of physical models and their solutions. Non-uniqueness is in direct conflict with the purpose of practical simulations, so that a mitigating strategy, outlined here, is important. We illustrate these issues in an examination of mesh converged turbulent statistics, with comparison to laboratory experiments.

## 1. Introduction

Practical simulations of fluid turbulence are inherently under resolved, owing to the very large range of length scale separating the finest vortices of turbulence and the scale typical of engineering interest. Simulations (of smaller, or more viscous flows) which do resolve the small vortex scales (the Kolmogorov scale) are called direct numerical simulations (DNS). Simulations which resolve some but not all of the turbulence range are called large eddy simulations (LES), while those that do not resolve the turbulence scales, and treat them via models are called Reynolds-averaged Navier–Stokes (RANS) simulations. RANS is in some sense the workhorse of practical engineering simulations, but the RANS models have undetermined parameters, which cannot be set universally. To fill in the gaps between experimental data, as is needed for the use of RANS simulations, LES are often used. LES, with all problem-specific details resolved, are more nearly universal, and have smaller data requirements, meaning that validation against experimental data for LES has a higher chance of being feasible.

LES thus have a fundamental role in engineering and scientific studies. To a good approximation, the physical transport parameters (viscosity, heat and species concentration) play only a small role in these simulations, and so in effect the solutions are properties of the Euler equation.

However, the Euler equation suffers from a disease of non-unique solutions, as we review here. Non-uniqueness of solutions is a fundamental obstacle to scientific predictions, preventing the success of such goals as predictive science, validation of solutions and quantification of the uncertainties on which engineers base their design conclusions. Some of the non-uniqueness of the Euler equation appears to be mathematical in nature, and in the optimistic tradition of practicing engineers and physicists, may be assumed to disappear after some future, deeper level of mathematical results have been obtained. However, some of the non-uniqueness is physical, or at least numerical in origin, owing to the under specification of the Euler equation on physical grounds. These deficiencies can be addressed only at the level of physics and engineering, with improved equations for LES, to fully specify the equations and thereby recover solutions which in principle might be unique.

## 2. Existence and non-uniqueness theories

Strong solutions are known to exist for short time and smooth data. We turn to weak solutions and beyond that (weaker yet), to measure valued (Young measure, or stochastic) solutions. The convergence in the LES framework is the convergence of a probability distribution function (PDF) to a probability distribution limit solution of the Euler or Navier–Stokes equations. Such a probability distribution limit, which is a function of space and time, is called a Young measure [1]. The existence of a Young measure limit for the compressible Euler equation was established in [2,3], although without identification of the equation satisfied by the limit. Starting from an assumed K41 upper bound on the velocity spectrum (the famous *k*^{−5/3} law), as an upper bound, a Sobolev norm bound with a fractional derivative can be derived. Assuming this bound, strong limits in an *L*_{p} space of the Navier–Stokes solutions (as the Navier–Stokes viscosity ) to a solution of the incompressible Euler equation were obtained [4].

Not all Euler solutions will obey a uniform *k*^{−5/3} spectral decay. For example, shock waves for the compressible Euler equation have a *k*^{−1} spectrum. Even for turbulent (and shock free) incompressible solutions, where the Young measure appears not to be needed, it is still a convenient method for the formulation of convergence of statistics, as we observe in §4.

The mathematical theory of non-uniqueness is uniformly negative. Starting with [5,6] and with recent developments reviewed in [7], non-unique solutions are known to abound, including solutions with non-zero values for a bounded time interval only. This non-uniqueness appears to be of a mathematical nature, but all efforts to remove it on this mathematical level (by imposing additional restrictions desired of a physical solution) have so far not eliminated the non-uniqueness. Still, regularity or admissibility restrictions appear to be a reasonable hope for an ultimate resolution. So the mathematical theories of non-uniqueness do not necessarily impede the use of the Euler equation, or LES of the Navier–Stokes equations in practical problems. Such pathological solutions have not been seen experimentally or in simulations.

More troubling, but also possible to address in a more direct fashion, is a type of non-uniqueness of a clear physical or numerical origin, as in this case a physical or numerical mitigation is possible.

In simulations, the solution value associated with a computational grid cell is not considered to be a point value at the centre of the grid cell. Rather it represents an average of the solution variables over the grid cell. For linear equations, this concept works well, but for nonlinear problems, averaging of the nonlinear flux term *F*(*U*) introduces an error between and . As the cell average is a primitive variable in the finite difference equation, we know the first of these quantities. But the averaged equation requires the second, for which no formula is provided. In practice, the cell block averages are replaced by filters, which are convolution averages of a few neighbouring cell blocks.

The strategy is to invent a model, usually proportional to ∇*U*, for the difference, with an unknown coefficient. By considering two adjacent grid levels, the current and a once coarsened level, and assuming that the solution is in an asymptotic (self-similar, scaling) regime, a space–time-dependent value for this missing coefficient can be determined [8–10]. The added term, with its coefficient, is called a subgrid scale (SGS) term. The physical/numerical role of the SGS term is to represent subgrid turbulent fluxes to the extent that they influence the resolved grid levels.

For conservative equations, the divergence of the SGS flux appears in the balance equations, and as the SGS term is modelled as a gradient, the result is an elliptic operator in the balance equations. For reasons of stability, the coefficient must be positive, and the net effect of the SGS turbulence model is to increase the diffusion in the solution.

There are two related issues which compound this picture. The first is numerical artefacts, which often also introduce diffusion into the solution. Even for higher order methods, these methods are invariably only first order near discontinuities or large gradients. The numerical truncation errors are of the same type and sign as the SGS terms. More serious is the unfortunate fact that they are typically large and can be larger than the SGS terms. The ultimate factor contributing to the solution morphology is not only these physical transport terms and the numerically based turbulent transport terms, but their dimensionless ratios, such as the Schmidt number *Sc*=*ν*/*D*, where *ν* is the viscosity and *D* is the species diffusion, with a similar definition for the dimensionless Prandtl number *Pr* governing thermal diffusion. The numerical problem arises from the fact that the dimensionless ratio of combined effects transport (physical, or molecular, turbulent or grid dependent; and numerical, or algorithmic dependent) is typically perturbed by the numerical truncation error contribution, and due to the sign restriction on the SGS turbulent transport term, this error cannot be compensated for by adjusting its coefficient.

Thus, we see that the SGS terms only address one half of the physics-based non-uniqueness. The other half is addressed by a high-quality numerical algorithm, which minimizes numerical error, especially in the equations where it is typically the largest: thermal and concentration diffusion. The front tracking addresses the numerical issue, which controls, or limits, diffusion in the presence of a sharp interface or steep thermal or concentration gradient. We call the algorithm FT/LES/SGS, as it combines and depends essentially upon front tracking, large eddy simulation and dynamically constructed SGS terms.

## 3. Numerical examples of non-uniqueness and its amelioration

Typical benchmark test problems for turbulent mixing are the Rayleigh–Taylor (RT), Richtmyer–Meshkov (RM) and Kelvin–Helmholtz (KH) instabilities. The first two involve flow directed normally to a density discontinuity, and driven by a force (e.g. gravity) (RT), or a shock wave (RM). The third is a shear layer or tangential velocity discontinuity across an interface. The FT/LES/SGS algorithm has achieved some of the most extensive validation results available, especially for RT (see [11,12] and references therein). In contrast to many codes [13], we find agreement between simulation and experiment for the dimensionless constant *α* that characterizes the growth rate of the outer edge of the mixing zone. Here we look at a second-order correlation
which characterizes a reaction rate between two species, with concentrations *f* and 1−*f*. The experiment [14] involves the flow of hot and cold water across a splitter plate. The heavier (cold) fluid enters the channel on the top and, due to gravity, it mixes (RT) as the flow passes through an observation window.

In figure 1, we compare values of *θ* from experiment and from two simulation codes. The comparison code, Miranda, is a high-quality turbulence code with 10th-order central differences, but no SGS terms. The grid resolution for the Miranda simulation is 1152×760×1280 and the grid size is Δ*x*=Δ*y*=0.025 cm and Δ*z*=0.01875 cm, which is approximately equal to the Batchelor scale *η*_{B}≈0.021 cm and finer than the Kolmogorov scale *η*_{K}≈0.055 cm. Thus, it is claimed in [14] to run at a nearly DNS level. DNS is properly defined as resolution sufficient to resolve all relevant physics. Thus, DNS depends on the code used as well as on what is relevant. Here, a normalized second moment is under study. There appears to be no theory or numerical experiment (other than the cited source) of mesh requirements. However, only one mesh is presented in [14], and as noted, it is not sufficient in the context of the Miranda code. So Miranda requirements for DNS are unknown and we surmise that the authors are comparing Δ*x* to the Batchelor scale in asserting nearly DNS convergence.

Compared to this is our FT/LES/SGS simulation, run at three grid levels, 8, 4 and 2 times coarser resolution than the Miranda simulation. In spite of this lower resolution, and with a lower order numerical algorithm, the FT/LES/SGS goes through all late time data points and shows a possibly converging trend. By contrast, Miranda misses all data points for the final two-thirds of the observational time. Our conclusion is that the extreme level of turbulence sophistication of Miranda (10th order, compact central differencing, nearly DNS mesh resolution) is less valuable than tracking in the context of this example. In our hot and cold water simulations, the interface between the two fluids is defined as an isotherm at the mid temperature between the ambient hot and cold temperatures.

## 4. Mesh convergence for flow statistics

In a series of papers [15–20], we have formulated a notion of convergence which appears to be useful in the study of turbulent simulations. The idea is based on the concept of Young measures. The idea is that the statistical fluctuations are rapid and stochastic in nature, while the probability distributions describing the statistics are slowly varying. Accordingly, this variation can be captured on a coarse grid, namely the grid of supercells. We describe the statistics in terms of probability density functions (PDFs) or their integrals, the cumulative distribution functions (CDFs). The convergence mesh data are organized into coarse-grained supercells. The coarse mesh supercell describes the spatial variation in the PDFs, while all the data in a single supercell represent a finite approximation to a Young measure concentrated there, and specifically to a PDF or CDF. A systematic test of these ideas is given in [17], where convergence norms are applied to the PDFs and CDFs. First, the differences (as the mesh is changed) in the PDF/CDF is measured by an *L*_{1} norm of the difference. These norm differences, one for each supercell, are compared (in *L*_{1} or ), to give an overall norm of the difference of PDF/CDF. We find convergence or converging trends under mesh refinement, with the order of convergence depending on the supercell size. From this fact, we conclude that the primary error is related to the statistical sample size as an approximation to the PDF/CDF, rather than being mesh related. These facts are illustrated with a table extracted from Kaman *et al*. [17] (table 1).

## 5. An application programming interface for front tracking

Having argued that FT/LES/SGS has desirable properties for the study of turbulent mixing, we turn to the question of the availability of this technology to others. Here, we address three issues: the application programming interface (API), a conservative version of the API and a graceful end to the tracking of the interface in overly complex frontal topologies.

Before going into the details, we describe the front or interface software, a product of Jiao and students [21–23], and replacing our older interface software. Conceptually the ideas are simple. A front is presented as an unstructured triangulation of a smooth surface (folds and self-intersections can be supported but are not discussed here). The triangles and adjacency information are stored. Associated with a front point (i.e. a triangle vertex) is a stencil of adjacent triangles extending a prescribed distance in terms of triangle adjacency. Relative to the tangent plane at the central vertex, the interface is represented by a height function *h*. *h* is approximated as *C*^{k} polynomials, leading to equations for the polynomial coefficients. The stencil is chosen too large, so that the equations are over determined, and they are then solved approximately by least squares. The result is a robust description of the curvilinear surface, which survives poor triangulation. From this *C*^{k} description, all differential geometry operations can be derived in a straightforward if tedious manner, including formulae for normal, curvature and quadrature, with the desired order of accuracy.

### (a) The application programming interface

The API [24] to support the current (well-validated) version of FT/LES/SGS is called front tracking interface (FTI). It has three essential operations. We first introduce the notion of a front state. Associated with a front point is a two sided set of states, namely solution state variables associated with the two sides of the interface, and connected or constrained by the laws of physics that the interface describes.

We are still experimenting with alternative versions of the algorithm, and here describe the simplest of the versions presently considered. The first essential step is to find the front points at the old and new time levels. The algorithm is second-order Runge–Kutta in time. Solution state values at front points are constructed from a one-sided interpolation (front aware interpolation). These are constructed at the old and (after advancing interior states) at the new time level. The normal velocity at the old and new levels is part of the Runge–Kutta point propagation.

The second essential step modifies the stencil states used in a one-dimensional sweep of the interior solver using a ghost cell extrapolation algorithm [25]. The modification, in case the stencil crosses a front, is to replace the stencil states on the wrong side of the front with front states from the right side. This step is removed in the conservative version of FTI (§5b).

The third essential step occurs when the cell centre crosses in time a moving front, so that the new and old time levels occur on different sides of the front. Here we use again a one-sided interpolation, i.e. interpolation that does not cross the front, to fill in the cell centre value. This step is removed in the conservative version of the API (see §5b). To illustrate the FTI operationally we list client functions needed by a physics code. Aside from utility I/O and initialization operations the main algorithmic additions required of a client code are a Riemann solver and a one-sided interpolation function, which uses states on one side of an interface only.

### (b) A conservative application programming interface

For the conservative API, we propagate the front points and the front itself to the new time level, as above. As the front moves over the time interval, it sweeps out in time a space–time three-manifold that cuts the space–time four-cells into the one-sided cut cells on each side of the tracked front. We refer to the new time level as the top and the old time level as the bottom, for a space time cell. Some of these cut cells may have small tops or no top at all. All such cells are merged with neighbours, until all cut cells have unit sized tops (space volume at the new time level). The conservative algorithm finds the average over the cut space cell at the top as the sum of a similar value from the bottom (old time level) plus the integral of *n*⋅*F* over the cut cell space–time swept interface surface. Here, *n* is the normal and *F* is the flux. This integration is performed to second-order accuracy using support functions from the interface package. To achieve second-order accuracy, it is sufficient to integrate *n*⋅*F* over the front at the half time level . Integration at this time level involves second-order accurate integration of a continuous function over a two-dimensional surface in three-space, a task which is supported by the current interface technology. The second-order accurate evaluation of the flux *F* on this surface needs a detailed analysis, while the second-order accurate evaluation of *n* is supported by the interface package. Likewise, the second-order evaluation of fluxes on regular cell surfaces near the interface, with stencils crossing the interface needs a detailed analysis. The algorithm is an improved version of the conservative tracking proposed in [26,27]. The primary improvement is in the simplified evaluation of *n*⋅*F*. Interface steps associated with repair of dynamically generated interface self-intersections are not addressed (yet) at a second order, conservative level by this algorithm.

### (c) A graceful termination to tracking

The step is simple: according to some criteria, we wish to declare a triangle to be untracked or passive. For example, the concentration or thermal gradient it is tracking may have become sufficiently diffuse that tracking is no longer advantageous. Or the interface may have become overly complex and resolving self-intersections is too difficult. It is not uncommon for the second condition to imply the first. In the interior state near front update, when we look for crossings of the front with a one-dimensional stencil line, we omit all passive triangles from the search. The interior stencil proceeds in its usual manner and ignores the passive portion of the front. Such passive triangles play no role in the balance of the computation (they are permanently passive), and optionally they can be removed from the front.

## 6. Non-uniqueness re-examined

We justified non-uniqueness of solutions as arising in numerical algorithms not aware of turbulent transport, or as due to numerical and algorithmic contributions to the turbulent transport (see §2).

A more fundamental argument is advanced in [28], namely that the *N*-species compressible Navier–Stokes equations have *N*+2 dimensionless parameters (Reynolds number *Re*, a ratio of the two viscosities (bulk and shear), a Prandtl number and *N*−1 independent Schmidt numbers). Transport properties of mixtures [29,30] depend on the mixture fraction and the individual species transport; the latter depends on the thermodynamic variables (pressure and density). To avoid complications not essential to present considerations (§3), we assume constant (and common) viscosity and thermal conductivity. Setting removes one dimensionless parameter and *N*+1 remain. So the solution can hardly be unique. This rather simple reasoning is connected to a renormalization group (RNG) derivation of the SGS terms, in which we recognize the SGS terms as RNG ‘essential variables’, whose value at each RNG step or mesh level is to be set from some external information. Uniquely for turbulence and not typical of RNG theories, this information is extracted from the dynamic simulation itself.

Both the SGS and the RNG reasoning are grid dependent and do not distinguish between non-uniqueness at the level of numerical approximations versus non-uniqueness of exact, fully resolved theories, except to the extent that the SGS terms satisfy scaling laws as observed in simulations (e.g. [20]).

Here, we take the point of view of the Kolmogorov scaling analysis, which does not assume any finite-scale numerical grid. To simplify matters, we consider incompressible flow coupled to a passive scalar, initialized with values 1 and 0 only, separated by some surface in three-space. Non-uniqueness is implied by a solution dependence of the order of limits as diffusion and viscosity . Agreement of the limits would suggest possible uniqueness of the resulting solutions of the governing equations. With the limit first, second, the mixing coefficient *θ*=0, as is easy to evaluate.

To gain insight into the limit in the opposite order, we appeal to the Kolmogorov scaling theory. Vortices with wavenumber *k* and radius *r*∼1/*k* have an RMS velocity *v*∼*k*^{−5/6} and a rotational period proportional to *r*/*v*∼*k*^{−1/6}. The spin up time for fully developed turbulence presumably involves wavenumber *k* vortices exciting the wavenumber *k*+1, etc., so that if this analysis is correct, the spin up time is infinite. Such a conjecture comes from Mailybaev and co-workers, who observed constant spin up time for one level of turbulence in starting from level *n* versus *n*+1 in large-scale DNS 31. For a single level of spin up the time difference between 1 and 2^{−1/6}=0.89 might be difficult to determine in a numerical simulation, so we regard his conjecture as consistent with our analysis. Now let us imagine that we have somehow initialized the fully developed Kolmogorov turbulent field, for example with some kind of fractal stirring rod. Even with fully developed turbulence, it would seem that the stirring would also have to progress successively, with the larger vortices stirring at their own length scales, and presenting smaller unmixed regions, which could be stirred by the smaller vortices afterwards. Large-scale stirring is needed to break up large regions into small ones. Small-scale stirring has to follow large-scale stirring as small-scale stirring cannot mix large regions in a bounded time. In other words, the stirring seems to involve the divergent integral , and with before , the mixing is not achieved and we still observe *θ*=0. In other words, within the context of an analysis based on an exact physical and mathematical problem formulation, the solution appears to be potentially unique, once numerical and mesh-dependent issues are removed from the analysis.

Of course, the mesh-dependent non-uniqueness survives and will be what governs practical considerations. Continuing on this theme, the above solution with only 1 and 0 as its values is no doubt a Young measure which does not have a PDF, and fails to satisfy a scaling law such as *k*^{−5/3}. However, in numerical practice, all solution fields (generally active rather than passive scalars) do satisfy scaling laws, the only contention being over the value for the scaling exponent. So we conclude that numerical solutions do not pick up the 0,1 solution in the double limit , . There is a validation counterpart to essentially the same issue. While the front tracking software will presumably give the 0,1 solution of this double limit problem, few if any other codes can be expected to do this. The issue has an experimental validation counterpart. Splitter plate experiments for fresh/salt water [32] model flow with a high Schmidt number *Sc*=600, with very narrow error bars for the measured value of *θ* as a function of time. Front tracking simulations [33] lie within these error bars, and no other researchers, using other codes, have even attempted the problem.

## 7. Conclusion

A goal of a theory of turbulent mixing is to provide input to various closure models which describe multiphysics problems in which mix has a role. Such models include RANS models of the mix itself and they include models for reactive flow and radiation in a mixture flow environment. Such goals will stretch existing knowledge of turbulent mixing and will require detailed properties (including mesh convergence) of the higher order single-point statistics and the length scales of two-point statistics. Results presented here can be viewed as a stepping stone in the direction of such a goal, but require further development before the goal is reached.

## Competing interests

We declare we have no competing interests.

## Funding

This research was supported in part by the US Department of Energy via Los Alamos National Laboratory contract no. 228022 and Stanford University contract no. 60541887107908 and the Army Research Organization grant nos. W911NF1310249 and W911NF1410482. This paper has been co-authored by Brookhaven Science Associates, LLC, under contract no. DE-AC02-98CH10886 with the US Department of Energy.

## Acknowledgements

The United States Government retains, and the publisher, by accepting this article for publication, acknowledges, a worldwide license to publish or reproduce the published form of this paper, or allow others to do so, for the United States Government purposes. This research used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the US Department of Energy under contract DE-AC02-06CH11357. Los Alamos National Laboratory Preprint LA-UR-14-29521. Stony Brook University Preprint SUNYSB-AMS-14-04.

## Footnotes

One contribution of 15 to a theme issue ‘Free boundary problems and related topics’.

- Accepted January 8, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.