CONSERVATION LAWS IN DISCRETE GEOMETRY

. The small length scales of the dissipative processes of physical viscosity and heat conduction are typically not resolved in the numerical simu- lation of high Reynolds number ﬂows in the discrete geometry of computational grids. Historically, the simulations of ﬂows with shocks and/or turbulence have relied on solving the Euler equations with dissipative regularization. In this paper, we begin by reviewing the regularization strategies used in shock wave calculations in both a Lagrangian and an Eulerian framework. We exhibit the essential similarities with Large Eddy Simulation models of turbulence, namely that almost all of these depend on the square of the size of the computational cell. In our principal result, we justify that dependence by deriving the evo- lution equations for a ﬁnite–sized volume of ﬂuid. Those evolution equations, termed ﬁnite scale Navier–Stokes (FSNS), contain dissipative terms similar to the artiﬁcial viscosity ﬁrst proposed by von Neumann and Richtmyer. We describe the properties of FSNS, provide a physical interpretation of the dis- sipative terms and show the connection to recent concepts in ﬂuid dynamics, including inviscid dissipation and bi-velocity hydrodynamics.


Introduction.
"It is wrong to think that the task of physics is to find out how nature is. Physics concerns what we can say about nature." Niels Bohr [48] "It was a wondrous alternative reality governed solely by the eternal laws of pure mathematics, unsullied by the crass realities of the world around us." Amir Alexander [1] The Euler equations, expressing the conservation of mass, momentum and energy, are the predominant model for numerical simulations of high Reynolds number flows. There are mathematical issues concerning the existence and uniqueness of solutions of those partial differential equations (PDEs). It is our purpose in this paper to show that those issues can be bypassed when the equation are posed in their more fundamental form as integral equations over discrete (e.g., finite) volumes. Our development is strongly focused on numerical methods, reflecting the fact that most interesting fluid flow problems can only be solved by computer simulation and that much experience has been gained in ensuring solutions that are stable and convergent. However, our principal result is an augmented set of PDEs (as opposed to discrete equations) that describe the evolution of finite volumes of fluids. Thus in our title we emphasize discrete geometry as opposed to discrete solutions.
The study of nonlinear hyperbolic conservation laws is both sophisticated and difficult in the sense of rigorous mathematics. In particular, much attention has been focused on the compressible Euler equations which express the physical principles of conservation of mass, momentum and energy. There is a practical reason for that interest, namely that numerical simulations of compressible fluid flows almost always solve the Euler equations which must, however, be regularized to yield physical solutions. It is the form and the magnitude of the regularization that is of practical importance.
From a physical point of view, the conservation laws are formulated in integral form. Then it is the derivation of partial differential equations (PDEs) from the integral form that engenders most of the mathematical difficulties. It is ironic that modern numerical methods for solving compressible fluid flows, and in particular those high Reynolds number flows that feature shocks and/or turbulence, are based on finite volume (FV) methods [34] which also solve the conservation laws in integral form. In other words, the discrete variables in a FV method are the volumeaveraged fields as opposed to the field variables evaluated at a point within the computational cell.
It is not our purpose in this paper either to contribute or to criticize the mathematical approaches for regularizing the Euler equations. Rather we wish to show how this step may be effectively avoided for the purpose of developing model equations for numerical simulations of high Reynolds number fluid dynamics. Our strategy can be summarized as answering the following question: If every point of a volume obeys the Navier-Stokes equations, what equations do the volume averages of those solutions obey? The difficulty in executing this strategy is the nonlinearity of the Navier-Stokes equations, which appears in the advective terms as well as the work term. The issue is that integration and multiplication do not commute, i.e., the average of the product is not the product of the averages. However, those commutation relations have been derived [41,42]. Using the commutation relations allows us to transform the integral relations for point quantities into PDEs for the volume-averaged quantities. We term those PDEs finite scale Navier-Stokes (FSNS) equations, which are presented in section 4, eqs. (14)- (16) and are suitable for FV simulation.
The FSNS can be interpreted as a regularized form of Navier-Stokes in which both the form and magnitude of the regularization are predicted. The essence is that the regularization depends on the square of the mesh size. This is completely consistent with each of the regularizations discussed in section 3, i.e., in the artificial viscosity used for shock simulations and in the turbulence models used for large eddy simulation. A fuller discussion of those regularizations, describing their motivations, their differences and their commonalities can be found in [38].
A penetrating analysis of the FSNS offered in section 5 suggests a physical explanation for the source of the new dissipative terms that connects to several modern theories of inviscid dissipation and bi-velocity hydrodynamics. That analysis will lead to some reinterpretation of Lagrangian motion of finite sized fluid parcels.
The plan of the paper is as follows. In section 2, we will discuss the Euler equations, the Navier-Stokes equations and their theoretical relationship in terms of the kinetic theory of gases. In section 3, we will progress to computer simulation, the issues that arise when all the dynamical degrees of freedom cannot be resolved on practical computational grids and the numerical strategies that have been developed to deal with those issues. Section 4 contains our main results, the derivation of the finite scale Navier-Stokes (FSNS) equations. This section also contains a brief description of the closure theorem, which is an essential element of the derivation. In section 5, we wlll offer an analysis of the FSNS and suggest a physical basis for new dissipative terms that have appeared. In that section, we will also discuss several other theoretical results and computational models that exhibit inviscid dissipation and bi-velocity hydrodynamics. We will conclude the paper with a short summary of our results and some remarks on the opening quotes of Bohr and Alexander.
2. Fluids in motion. The Euler equations [15] and Navier-Stokes equations [53] were originally derived on the basis of classical mechanics and experiment. Here, we briefly describe those two sets of evolution equations and their relationship in the more modern context of the kinetic theory of gases [31].
2.1. Euler and Navier-Stokes equations. We start by writing the Navier-Stokes equations in conservation form and employing internal energy rather than temperature. In one spatial dimension, those equations are: Here, ρ,u, and I have their usual meanings of mass density, material velocity and specific internal energy. The total energy E ≡ I + 1 2 u 2 . The momentum flux is the sum of the thermodynamic pressure p plus the viscous stress P: where γ is the ratio of specific heats and µ is the longitudinal viscosity coefficient. The heat flux q is given by Fourier's law written in terms of the internal energy rather than the temperature, where κ is the coefficient of heat transport. Note that the only length scale in the Navier-Stokes equations is contained in the transport coefficients µ and κ; that length scale is the molecular mean free path . Both transport coefficients are dimensionally proportional to ρc where c is the speed of sound. is typically much smaller than the length scales of practical problems. In contrast, by direct numerical simulations (DNS) of turbulence or shocks, one implies computer simulations in which is well resolved by the computational mesh. We also define the specific thermodynamic (equilibrium) entropy where c v is the specific heat at constant volume. S 0 is a reference entropy as the entropy is only defined within an additive constant. The balance equation for entropy, in local form and written for open systems, is where the third term on the LHS is the entropy flux [44] and Γ is the dissipation function [54]. The second law of thermodynamics dictates that the dissipation function is nonnegative, Γ ≥ 0, which then leads to the form known as the Clausius-Duhem inequality ∂ρS ∂t Here we have assumed that there are no internal sources of energy production. The thermodynamic entropy of eq. (6) satisfies the Clausius-Duhem inequality in all Navier-Stokes simulations, i.e., in the evolution of nonequilibrium systems [54]. The Euler equations result 1 from setting the transport coefficients to zero, i.e., µ = κ = 0. Then, the case of the Euler equations, the heat flux q and the dissipation function Γ both vanish. Thus for the Euler equations, the Clausius-Duhem inequality becomes an equality: However, we note that this last equation also depends on the smoothness of the solutions of the Euler equations. In [16], Drivas and Eyink quantify sufficiently smooth in terms of Besov regularity. Their principal result states that solutions that are not sufficiently smooth will exhibit anomalous dissipation 2 and nonnegative entropy production. Those comments strike at the heart of this paper. From the mathematical point of view [13], the strong solutions of the Euler equations even from well-posed initial conditions will form discontinuities (shocks) and subsequently will exist only in the sense of weak solutions (i.e., generalized functions). However, those weak solutions are not unique and some "entropy" principle must be invoked to identify the physically relevant solution. It is the convergence of that limiting process (as the transport coefficients vanish) that leads to the mathematical difficulties alluded to in the introduction. It is our purpose to show how this limiting process can be avoided for the purpose of finding entropy stable numerical algorithms for finite volume schemes.

Kinetic theory.
In kinetic theory, the Euler and the Navier-Stokes equations represent different orders of approximation in the Chapman-Enskog hierarchy [31] of perturbation approximations of the Boltzmann equation. The Euler equations represent the lowest-order approximation in which the velocity probability distribution (PDF) is the equilibrium Maxwell-Boetzmann distribution. The next order of perturbation, formally in terms of the Knudsen number [31], yields the (compressible) Navier-Stokes equations which have proven to be an accurate model for studying most fluid flows since their inception [53]. Navier-Stokes has also proved a useful model for numerical simulation, but with the notable exception of high Reynolds numbers flows in which shocks typically appear. In this paper, we will be principally concerned with those high Reynolds number flows.
The Knudsen number is defined as a ratio of length scales, K n = /L, where is the molecular mean free path and L is a macroscopic length scale characterizing a problem. Thus K n is a microscopic analog of the Reynolds number that indicates when a continuum approximation is appropriate. However even in strong shocks, one finds K n ≤ 1/3. Experimentally measured shock widths are always more than three mean free paths wide, see e.g., Fig. 10 in [51]. On this basis, one might expect the Navier-Stokes equations to do a good job predicting shock structure and further expect that higher-order Chapman-Enskog approximations, e.g., the Burnett equations, will improve the prediction. In fact, neither expectation is correct; the Navier-Stokes equations fail to predict either the width or the shape of the shock profile as measured experimentally and use of the Burnett equations does not improve that lack of agreement. Specifically, the failure of Navier-Stokes equations to describe measured shock profiles was first predicted theoretically by Becker [3] and verified experimentally by Schmidt [51] and Alsmeyer [2] among many others. The critical issue is the additional assumption of local thermodynamic equilibrium upon which both the Euler and the Navier-Stokes equations are based. The further results of simulations based on the Burnett equations are discussed in [18]. A careful study of shock structure, and associated entropy issues can be found in [39].

Why Euler?
The Reynolds number estimates the ratio of two length scales, the large scales of advective processes to the smaller scales of dissipation. Representative Reynolds numbers in such diverse fields as weather prediction and astrophysics easily reach and exceed 10 6 ; thus it is not possible to fully resolve that range of scales on computers of today. Since achievable resolutions will not resolve the scales of dissipation, one is effectively solving the Euler equations even when the Navier-Stokes equations are programmed. Fortunately, in most cases it is only the larger scales of the flow that are of interest. Then the details of the dissipative processes can be ignored, although some form of dissipation must be added to the Euler equations to ensure a stable and convergent numerical solution. This is what is meant by regularization.
As described in the previous section, the strong solutions of the Euler equations even from well-posed initial conditions form discontinuities (shocks) and subsequently exist only in the sense of weak solutions (i.e., generalized functions). The issue is that those weak solutions are not unique, and a so-called entropy principle must be invoked to identify the unique physical solution. One such entropy principle is based on the method of vanishing viscosity [5] in which an infinitesimal amount of physical viscosity is added to the Euler equations and the physical solution is identified as the unique limit of a convergent series of solutions as the viscosity vanishes. An extension of the ideas of weak convergence is the subject of compensated compactness [14].
Unfortunately, as there are no infinitesimals in numerical methods, the ideas vanishing viscosity cannot be applied directly to numerical algorithms. In the next subsection, we will summarize two widely used approaches to regularizing the Euler equations for numerical simulation. We will carefully distinguish between the model equations and the discretized equations. With artificial viscosity [56], the model equations are augmented by an explicit dissipation. Artificial viscosity is still the main methodology for regularizing shocks in the Lagrangian framework. In nonoscillatory finite volume (NFV) methods, the Euler equations are the model equations and dissipation is introduced implicitly through the approximation of the advective terms. NFV methods are the consensus choice for high Reynolds number simulations in the Eulerian framework.

Discrete regularization.
"It is frequently desirable to solve the equations of fluid motion by stepwise numerical procedures, but the work is usually severely complicated by the presence of shocks." von Neumann & Richtmyer [56] "Just because we cannot prove that compressible flows with prescribed initial values exist doesn't mean that we cannot compute them." Peter Lax [32] The earliest applications of computer simulation reflected the interests of computational fluid dynamics (CFD) pioneer John von Neumann. The physical regimes of weather prediction and of nuclear weapons simulation may appear very far apart, but both are examples of high Reynolds number flows. The earliest shock calculations were Lagrangian, meaning that the cells of the computational grid moved with the local fluid velocity. It was quickly discovered that it is not sufficient to solve the discretized Navier-Stokes equations in those calculations. The velocity gradients estimated on grids with cells much larger than the width of the physical shock do not produce sufficient dissipation of the kinetic energy; then it is observed that unphysical oscillations of the velocity field build up behind a shock leading ultimately to numerical instability and possibly calculation failure.
In 1950, von Neumann and Richtmyer described a conceptually simple solution for dealing with unresolved shocks; in [56] they proposed to add an artificial dissipation, now termed artificial viscosity, to broaden the shock to the width of several computational cells. Although the artificial viscosity is often introduced as a numerical artifice, see e.g., [9], the artificial viscosity has some physical justification. It is noted in [56] that the shock speed and jump conditions do not depend on magnitude of the viscosity. It is only the detailed shape of the shock profile that varies.
The form of the von Neumann-Richtmyer artificial viscosity in 1D Lagrangian simulations is [57] where ρ is mass density, u is material velocity, c 2 is a dimensionless constant and ∆x is the size (length) of a computational cell. Q is added to pressure in the (discrete) momentum and energy equations. The dependence on the square of the velocity gradient and the square of the cell size is noteworthy and is not justified in [56]. Note that the FSNS momentum equation (15) derived in section 4 predicts both the form and the magnitude of the artificial viscosity. Artificial viscosity methods remain the principal Lagrangian methodology for flows with shocks. Lagrangian codes 3 are not conducive to simulations of turbulent flows due to large distortions of the Lagrangian grid; simulations of turbulent flows are mainly performed in the Eulerian framework where the computational grid remains fixed in space and fluid can move from cell to cell. Similar to the experience with shock calculations, the earliest weather simulations exhibited unphysical oscillations. Joseph Smagorinsky developed a three-dimensional version of the von Neumann-Richtmyer artificial viscosity to control those oscillations, giving rise to the first large eddy simulation (LES) subgrid scale model [50]. The development of subgrid scale models continues to be an area of active research, but the fundamental dependence on the square of the velocity gradient persists [49].
The earliest applications of Eulerian codes to shock flows also used artificial viscosity, but were considered inferior to the Lagrangian codes due to excessive smearing of the shock profiles. That situation changed dramatically with the innovation of nonoscillatory finite volume (NFV), beginning with FCT [6] in 1973 and rapidly followed by high-order Godunov, [55] , TVD [26] and MPDATA [52]. Today there is a virtual alphabet soup of NFV methods. The basic idea of NFV is to preserve monotonicity in the numerical solution. This is accomplished by mixing low-order and high-order approximations for the advective terms with the mixing coefficient depending on the local flow. NFV methods use the Euler equations as their model equations and gain dissipation implicitly through their approximation of the advective terms. This is in contrast to the artificial viscosity methods, which add dissipation explicitly in the model equations.
The success of NFV methods in representing shock flows led Jay Boris to speculate that his FCT methods might also prove effective as an implicit subgrid scale model for turbulent flows [47]. This idea, now termed implicit large eddy simulation (ILES), was initially strongly opposed by the traditional turbulence community. Although the quality of results produced by ILES simulations using many different NFV algorithms speaks for itself, the turbulence community required a theory to justify the technique. Such a rationale for ILES was provided in a paper by Margolin and Rider [41]. That paper is based in part on a modified equation analysis 4 (MEA) of MPDATA [52], one of the NFV scheme schemes that has successfully reproduced many experimental results and computational benchmarks of the LES community. In [41] it is shown that that MPDATA implicitly contains the artificial viscosity of eq. (10). The same authors extended that result to include many other NFV schemes in [23].
To summarize this section, • Numerical methods for solving high Reynolds number flows solve the regularized Euler equations. • The regularization takes principally two different forms, the explicit regularization of artificial viscosity methods and the implicit regularization of NFV methods. • In the sense of MEA, both the explicit and implicit regularizations have the form of eq (10), being quadratic in the velocity gradient and in a finite length scale.
However, the landmark paper of von Neumann and Richtmyer [56] does not derive this unique form. In the next section, we will show that this "magic" term arises naturally in the derivation of the finite scale equations that govern the motion of finite-sized volumes of fluid. In section 5, we will offer a physical interpretation of the source of the magic.

Finite scale Navier-Stokes.
Here we answer the question posed in the introduction, "what equations govern the evolution of volume-averaged variables". First we will define the volume-averaged variables. Then we will briefly describe the derivation of the closure theorem, which provides commutation relations for switching the order of multiplication with volume averaging. Application of the theorem then leads to the finite scale Navier-Stokes equations.
4.1. Volume averaging. The foundation of finite scale theory is the existence of a finite length scale that is not inherently determined by the physical processes being modeled. In the case of numerical methods, it is the computational cell size, ∆x; however, the theory is general and is not constrained to depend on any relation to numerical simulation. With a view toward use in developing models for numerical simulation, we define the volume-averaged variables by integrating over the computational cell; then, for example, the volume-averaged (macroscopic) density in one spatial dimension (1D) is defined: where volume-averaged quantities are indicated with a hat. Further, we will use a tilde to indicate density-averaged (Favre) variables, which will also naturally appear in the theory; e.g., the density-averaged velocity in 1D is: We note that the both the volume-averaged quantities and the Favre averaged quantities are smooth continuous functions.
To derive the finite scale equations, we will apply the averaging operator defined in eq. (11) to each of the Navier-Stokes equations (1)- (4). Note in the averaging process, nonlinear terms arise, e.g., in the advective terms, that must be resolved. For example, ρu 2 = ρ u 2 . That devolution is the subject of the following closure theorem.

4.2.
Closure theorem. The compressible FSNS that are derived below follow from the following general closure theorem, presented in 1D below, but readily extended to three spatial dimensions and to the time domain: Closure Theorem: for any (continuum) fields A and B that are sufficiently smooth at small scales, where A x ≡ ∂A ∂x , HOT are higher-order terms, e.g., O(∆x 4 ), etc., and averaging (overbars) are defined as in eq. (11).
In words, the derivation begins with the assumption that physical viscosity ensures that the flow field may be expanded in a convergent Taylor series at sufficiently small length scales. In that case, the integrals of the expansion terms may be explicitly evaluated; note that by symmetry, only the terms with even powers of ∆x survive. Equation (13) is easily verified for ∆x in this regime; see [41] for details.
The proof proceeds with a combination of induction and renormalization for larger values of ∆x. That is, we assume that eq. (13) holds for some large value of ∆x and prove that the same form then holds for 2∆x. It is important to recognize that the meaning of the 'hat' depends on the averaging length scale, and so one is effectively doing a change of variables in the renormalization. The invariance of the parametric dependence of the equations is termed form invariance. Details of the proof can be found in [41] where the time domain is considered and in [42,23] where multiple space dimensions are considered.

4.3.
Presenting the finite scale Navier-Stokes equations. The derivation of FSNS begins by applying the averaging operator of eq. (11) to the one-dimensional version of the Navier-Stokes equations (1)-(3) and then commuting the order of operations using the closure theorem (13) to resolve the nonlinear terms. We choose the averaged density ρ, the averaged momentum density M = ρu, the averaged total energy density E ≡ ρE and the averaged pressure p as the primary dependent variables. Then the following intermediate form of FSNS is derived: Note that the closure theorem indicates that terms of order ∆x 4 and higher will appear in general. The FSNS eqs. (14) -(16) represent a truncation to order ∆x 2 of an infinite series of terms in even powers of the cell size. As a practical matter, implementation of the FSNS in a FV computer program will also lead to higherorder terms in the sense of the modified equation [28]. It is for this reason that there are many useful formulations of NFV schemes for shock propagation and of subgrid scale models (both LES and ILES) for turbulent flow, as described in section 3.2.
We begin an analysis of these equations by considering the continuity equation (14). The appearance of u and its derivatives cause an immediate problem here. u cannot be a prognostic variable because we have chosen momentum as a primary variable.
To proceed, we use the closure theorem to write to order ∆x 2 : and then define a momentum velocity u by Note that u is just the Favre-averaged velocity, cf. eq. (12). Similarly, we define the Favre-averaged total energy and internal energy Using eq. (18), the continuity equation can be written in the more familiar form of a conservation law: The revised momentum equation can be rewritten in terms of u as follows. From eq. (18) we have M = ρ u (21) so that Also, since u and u differ by terms of order ∆x 2 , we can write the approximate inverse of eq. (18) Substituting for M and u using eqs. (22) and (23) in eq. (15) then leads to where The finite scale momentum flux P consists of the classical Newtonian viscosity and a new term of order O(∆x 2 ) that arises from the averaging of the momentum advection.
Similarly, the energy equation can be rewritten in terms of u, E and I: where the finite scale heat flux The finite scale heat flux q consists of the classical Fourier heat flux and two new terms. The first new term arises from averaging the nonlinear energy advection and is also O(∆x 2 ). The second new term arises from the averaging of the nonlinear work term. In [35], it is shown to have the same form as the first term in the case of a gamma law gas. More generally, this term will depend on the form of the equation of state, but will also be O(∆x 2 ). The averaged pressure is given by Here, the appearance of S ≡ 1 2 ∆x 2 2 ( u x ) 2 merits further discussion. In the Navier-Stokes equation (3), the total energy E = I + 1 2 u 2 so that after the finite scale averaging, one has However, the specific volume-averaged kinetic energy is defined

CONSERVATION LAWS IN DISCRETE GEOMETRY 197
The term S in eq. (28) represents unresolved kinetic energy. In other words, in the finite scale equations there are three partitions of the total energy, namely the averaged internal energy, the averaged kinetic energy and the unresolved kinetic energy A more complete derivation of these relations from kinetic theory can be found in [37]. We note the similarity of this equation to energy measured in the Sobolev norm which is a starting point for the LANS-alpha model [30].

A physical interpretation of artificial viscosity.
"Everything happens for a reason. It's usually physics." (George Takei) The first term in eq. (25) in fact is Q, compare with eq. (10). It has arisen straightforwardly in the mathematical derivation of the finite scale equations and the recognition that finite volume algorithms are solving those finite scale equations. However, it is now possible to provide a more physical interpretation.

Lagrangian parcels.
Let us distinguish between a Lagrangian parcel, meaning a finite volume of fluid whose boundaries move with velocity u, and a Lagrangian point meaning a mathematical point of fluid that moves with velocity u. Integration of eq. (20) over the fluid parcel leads to a more general notion of Lagrangian motion. If the boundary of the parcel moves with the velocity u, then eq. (20) implies that the total mass inside the volume remains constant. This does not preclude the flux of mass across the boundary, but rather means that the surface integral of the flux around the entire boundary vanishes. That is, the parcel exchanges mass with the surrounding fluid, but there is no net change of mass of the parcel. However, the flux of mass implies a concomitant flux of momentum and of energy associated with the exchange of mass. In general, those integrals do not vanish. In other words, the net zero exchange of mass with its neighbors will lead to a nonzero exchange of momentum and energy. It is the artificial viscosity that accounts for that nonzero exchange of momentum. Similarly, the artificial heat conduction of eq. (27) accounts for the nonzero exchange of heat accounts for the nonzero exchange of energy. Such a term has been proposed by Noh to deal with the troubling phenomenon of wall heating [46], although it is not presently used in Lagrangian simulations.
Another inference of the nonzero mass exchange is that all fluid particles initially within a particular parcel may eventually leave that parcel. That is, one cannot label the parcel with any of the particles initially within the parcel. Particle exchange is reminiscent of the particle relabeling symmetry [45].
One advantage of using FSNS rather than Navier-Stokes can be exposed in Eulerian simulations of turbulence by estimating the effective Reynolds number. In Navier-Stokes theory, the Reynolds number is defined [54] Re where U is a typical velocity of the flow and L is some (often nebulous) characteristic length scale of the flow. The Reynolds number is interpreted as a measure of the relative importance of advective and diffusive processes in a flow; large Re indicates diffusional instability and turbulence. However, if one identifies L with ∆x, it becomes clear that the Reynolds number is not only a property of the flow, but also of the discretization. Further, when considering the FSNS at scales where the observer is much larger than the molecular mean free path, ∆x >> , the physical viscosity can be neglected; then, an effective Reynolds number can be defined by replacing µ by 1 3 ∆x 2 Here, we have made the simplest order of magnitude estimates, e.g., U x ∼ U/∆x. Now since both the advective flux and the diffusive flux depend on the same length scale, ∆x, the effective Reynolds number is no longer a ratio of length scales, but is a constant. A critical Reynolds number for transition to turbulence, e.g., in pipe flow, is of the order of thousands [54]. Thus, the effective Reynolds number is small and the flow solutions of FSNS are not turbulent. This is the essence of why the ILES methodology produces stable results without requiring any explicit subgrid scale turbulence model. We note that the same scaling arguments apply for most conventional LES models when the filter length is chosen proportional to ∆x.

5.2.
Inviscid dissipation. The energy dissipation due to the momentum flux of eq. (25) aligns with several results of inviscid dissipation in the classical theory of high Reynolds number flows. From eqs. (26) and (25), the global rate of energy dissipation by the momentum fluxes in a one-dimensional domain D is; Integrating by parts and neglecting surface terms (work done by external forces) Here the brackets signify spatial integration over the domain. Then, in the limit of vanishing viscosity (µ → 0) and interpreting u x ∆x ≈ ∆ u, we have that the inviscid energy dissipation is proportional to (∆ u) 3 . This same proportionality was derived by Kolmogorov for (incompressible) turbulent flows and is known as the (4/5) th law, see e.g., [17]. The ability of an ILES calculation to reproduce the law is demonstrated in [42].
The dependence of energy dissipation on (∆ u) 3 has also been noted in shocks. Bethe [4] showed that the rate of entropy production across a shock is: where T is the temperature. Recall that the change in the internal energy I is given by the first law of thermodynamics Also, Frisch [17] notes a similar scaling for the dissipation of kinetic energy in the case of shocks in a Burgers' fluid. Note that the inviscid energy dissipation scales with ∆x 2 . As the finite volume grows, the inviscid dissipation increases as opposed to the viscous dissipation which grows smaller, cf. section 2.2. 5.3. Alpha models. The alpha models of Holm and many collaborators to which we add the older regularization of Leray [33] have many similarities to FSNS. In this collection of models designed to regularize turbulent flows, the advective terms are modified so that the point velocity u is advected by a "smoothed" velocity u (in our notation). The resulting equations depend on both those velocities so that an additional constitutive law is required to relate them. In simulations using the Leray regularization the averaging scale is taken to be order of the grid size ∆x as in FSNS. In LANS-alpha, the averaging is done on the length scale α which is the Lagrangian correlation length. The latter is a property of the flow whereas in FSNS, ∆x can be chosen independent of the flow. However, one would expect to choose the computational grid to resolve important macroscopic features of the flow so that the two choices may not be dissimilar in practice.
Both the Leray and LANS-alpha have rigorous theoretical derivation [33,30] which is essentially different from our derivation of FSNS in section 4. However, all three approaches depend on presuming a nonzero length scale over which to average, a necessary ingredient to ensure the existence of two distinct velocity fields. Both Leray and LANS-alpha have been successfully employed in numerical simulations, see [19,20,10,11] among many other examples. The relation between Leray and LANS-alpha is discussed in [24]. Also, we note a recent synthesis of Leray and LANS-alpha in [12]. 5.4. Volume transport models. The volume transport model of Howard Brenner [7] proceeds from yet a different direction, being based on classical mechanics and thermodynamics rather than volume averaging. Brenner questions the assertion that the velocity that appears in eq. (1) which he terms the mass velocity, v m , is the same velocity that appears in eq. (2), which is the material velocity v. In [8], he writes that the assumption v m = v is "a relation introduced into the foundations of rational fluid mechanics by Euler [15] in 1755, and seemingly beyond question ever since." Brenner then discards the assumption of equality and shows there is no logical contradiction in an alternative theory in which these velocities are not equal. Brenner does not perform numerical calculations, but relies on experimental validation that is centered in slow moving low Reynolds number flows where variations in density are primarily caused by gradients of temperature rather than of pressure.
If, as Brenner argues in [7] v m = v, then a constitutive relation between the two velocities is required, similar to the situation in the alpha models. Brenner proposes [22] a constitutive relation which in one spatial dimension is: Here α v is a constant termed the "volume diffusivity". Dimensionally, α v contains a length scale. Although Brenner works in the context of the Navier-Stokes equations, his theory has a nice interpretation within the FSNS. If we identify our velocity u ∼ v m and u ∼ v, then eq. (17) suggests an interpretation for the volume diffusivity Even though Brenner has been concerned mainly with low speed flows, his modifications would have large effect in flows with appreciable density gradients, e.g., shocks. As described in section 2.2, the Navier-Stokes theory does not capture the experimentally measured shock structure for shocks of strength Mach 2 and above. Greenshields and Reese [22] applied Brenner's model to simulate shock structure and write: "While it is important not to draw strong conclusions based on just one test case, our results are generally encouraging for the Brenner-Navier-Stokes equations." We add that those authors found a large sensitivity to the value of α v which they assumed was constant.
6. Discussion. In this paper we have been primarily concerned with deriving algorithms for the numerical simulation of high Reynolds number flows, e.g., flows with shocks and/or turbulence. To summarize the discussion in section 3, current algorithms approximate the Euler equations, regularized either by the explicit methodology of artificial viscosity, or on the implicit methodology of nonoscillatory approximation; in both cases, calculations can be made stable and convergent. However convergent is not the same as accurate, especially when the underlying model does not guarantee unique solutions. Such is the case for the Euler equations which are known to allow multiple (weak) solutions. The strategy then is to regularize the equations by adding dissipation. It is the form and magnitude of that regularization that is at issue. A recent overview of many regularization techniques as applied both to the continuous Euler equations and to discrete numerical algorithms can be found in [38].
In section 4, we have pursued an alternate strategy. Noting first that the conservation laws are fundamentally written in integral form and second, that most modern numerical schemes use finite volume approximation, it seems unnecessary (for practical purposes) to treat the equations in PDE form, i.e., not necessary to take the limit as the integral domains shrink to a mathematical point. Instead, we have derived the evolution equations for the averages of the field variables over finite domains. This was accomplished by integrating the Euler equations (or the Navier-Stokes equations) over a finite volume and then permuting the order of operations of spatial and temporal differentiation with the volume integration.
There are subtleties in dealing with nonlinear advective terms as integration and multiplication do not commute. Indeed, it is precisely that lack of commutativity that gives rise to inviscid dissipation -i.e., dissipation independent of the physical viscosity. The key to this derivation is knowing the commutation relations. We have used a closure theorem, briefly described in section 4.2 with several references to the detailed derivation in the literature. Application of the theorem to the Navier-Stokes equations leads to our principal result, the finite scale Navier-Stokes (FSNS) in section 4.3. Those commutation relations depend on the domain of the integrals ∆x, which is how that length scale enters into the FSNS.
In section 5, we provide some physical interpretation of the FSNS. The concepts of Lagrangian motion are generalized for discrete fluid parcels (finite volumes) which move with the Favre-averaged velocity. We show that the energy dissipated by the inviscid flux terms is congruent with previous theoretical results of inviscid dissipation in shocks and turbulence.
In section 5.3 we discuss similarities of FSNS to alpha models of turbulence which also result from adding a finite length scale through spatial averaging. One important practical difference is that the FSNS can be closed without reference to the particle velocity u. LANS-alpha models have been tested successfully in numerical simulations of turbulence, but have not, to our knowledge, been applied to flows with shocks.
In section 5.4, we discuss the volume transport model which has been developed in the context of classical continuum theory and thermodynamics. In the underlying theory, Howard Brenner questions an implicit assumption of Euler that the velocities that appear separately in the mass and momentum equations of Navier-Stokes are the same. Brenner does not explicitly introduce an additional length scale, but needs to define a volume-diffusivity parameter to develop a constitutive relations between his velocities. We show that his constitutive relation can be interpreted in the sense of the FSNS relation between the average velocity and the Favre average velocity, thus implicitly invoking the length scale ∆x.
We conclude this paper with some thoughts in the context of the quotes that open our introduction. Amir Alexander is a mathematician and a historian. In "Duel at Dawn" [1] he writes about the development of modern mathematics in Europe, describing the transition from the "gèométres" of the eighteenth century (e.g., Euler, Lagrange, D'Alambert) to the "algébristes" of the nineteenth century (e.g., Cauchy, Galois, Abel). Where the gèométres were motivated by the problems of the real world, the algébristes worked in the pursuit of abstraction and mathematical truth.
Physicists today debate whether nature at the finest length scales is continuous or discrete. However, Bohr reminds us that in either case, measurements of nature are invariably done with finite instruments and are limited by discrete scales of length and time. It would seem that discrete equations are the more effectual language for physics.
More specifically, the conservation laws of continuum mechanics are formulated in terms of integral balances. So long as the integral domains remain finite, the discrete scales appear in both the equations and in their solutions. It is only in the limiting process that leads to the PDEs where the mathematical difficulties arise. However, PDEs are easier to analyze. We believe that an important contribution of the finite scale theory is to demonstrate the equivalence of the PDEs that govern the evolution of the integral averages of the conserved field variables with the integral equations that govern the evolution of (experimentally unmeasureable) field variables at a point.

DEDICATION
For Darryl from Len: After more than 50 years of science and friendship, hiking and chocolate cake, we still grow old together.