Observable Divergence Theorem:
Evolution Equations for Inviscid Regularization of Shocks and Turbulence

For the last several years we have been promoting the idea of simultaneous inviscid regularization of shocks and turbulence in inviscid flows. While shocks and turbulence are often treated as separate problems ¹¹1Shocks are treated as discontinuity formation in hyperbolic equations in compressible flows while turbulence is considered as a closure problem which occur in both compressible and incompressible flows., they share some similar characteristics. The primary shared feature between them is a high wavenumber cascade processes created by the nonlinear convective term $\mathbf{u}\cdot\nabla\mathbf{u}$. Traditionally, this high wavenumber irregularity is regularized and controlled by the addition of a viscous term or a random walk in the Euler equations. While successful in describing many complex flow phenomena, this approach is still facing significant theoretical and computational challenges in high wavenumbers (or high Reynolds number flows)LionsPL:98ab . In this investigation, we take a different approach and explore a formal derivation of an inviscid regularization of this high wavenumber behavior in shocks and turbulence. The essence of the technique is in the re-derivation of the divergence theorem of Gauss for field quantities that are averaged. The divergence theorem of Gauss in calculus is revisited in order to prove an analogous theorem for a control volume where the fluxes at the faces of the fluid element are calculated based on averaged flow quantities over a certain length scale as opposed to a fully resolved one.

Our work here was motivated by our earlier results in inviscid regularization of the 1D and multidimensional Burgers equations Mohseni:08w ; Mohseni:09h where the nonlinear convective term is replaced by $\overline{\mathbf{u}}\cdot\nabla\mathbf{u}$; $\overline{\mathbf{u}}$ is an averaged velocity. We have shown that under certain conditions the solution to our regularized Burgers equation exists and is unique at all times and it converges to the entropy solution of the Burgers equation. Norgard and Mohseni Mohseni:10k has recently extended some of these results to 1D compressible Euler equations. Supporting numerical simulations were presented and it was shown that these equations share the same traveling wave solutions and entropy solutions with the Euler equations. This manuscript provides a formal derivation of our regularized Euler equations from basic principles and its extension to multi-dimensional flows. It should be noted, while a regular viscous term could be included in our analysis without any difficulty that is not necessary for our purposes here.

In the following sections we are required to calculate the flux of some conserved quantities which depends on an averaging process. A general class of averages for a function $f$ can be defined by the convolution operation as $\overline{f}=g\ast f,$ where $g$ is the kernel of the convolution. In such an averaging there are several guidelines that seem intuitively reasonable; see Mohseni:08w ; Mohseni:09h . As listed in Table 1, the averaging kernel should be non-negative, decreasing and radially symmetric. Such a filter $g$ can be associated with a characteristic wavelength $\alpha$. This parameter is introduced by scaling the filter as such, $g^{\alpha}=\frac{1}{\alpha}g\left(\frac{\mathbf{x}}{\alpha}\right)$, or in Fourier domain $\widehat{g^{\alpha}}(\mathbf{k})=\widehat{g}(\alpha\mathbf{k})$. Thus as $\alpha$ becomes smaller, the wavelength where the filter exerts influence also become smaller. With this scaling, the filter remains normalized, non-negative, decreasing and isotropic. Furthermore, as $\alpha\rightarrow 0,\,g^{\alpha}$ approaches the Dirac delta function. When convolving $g^{\alpha}$ with $f$, features in $f$ that have length scales less than $\alpha$ will be averaged out.

Our ability to observe and measure (physically or numerically) any fluid property is limited by the spatial resolution of our experimental or numerical technique. For instance, if one uses hotwire annemometry for measuring turbulence velocity components, the spatial resolution of the measurement is always limited by the size of the wire. In fact, what is actually measured is an averaged value of the velocity over a length scale of the order of the size of the hotwire. Similarly in numerical simulations, one is limited by scales of the order of the mesh size. As a result, what we often observe is an averaged quantity rather than the mathematical limit where the observation volume approaches zero. This observation has significant ramifications when one calculates the flux of a quantity at the face of a volume surrounding a point of interest. Before we investigate this, we will revisit the divergence of a vector field and the divergence theorem of Gauss. We will then extend these concepts to the case of a fluid volume with limited observability.

In writing the differential form of the conservation laws for a fluid, one considers the flux of the conserved quantities from the faces of a fluid volume. These fluxes are then related to the divergence of the conserved quantity inside the volume. Let $\Omega$ be a small and closed region centered at point $\mathbf{x_{0}}$ whose boundary is represented by the orientable surface $S$ and its volume is $\Delta V$. Let $\mathbf{F}$ be a vector function in $C^{1}(\Omega)$. We denote the flux of the vector field $\mathbf{F}$ through the surface $S$ by $\Phi(\Omega)$. The divergence of the vector field $\mathbf{F}$ around the point $\mathbf{x_{0}}\in\Omega$ is defined as the limit of the flux density of $\mathbf{F}$ over $\Omega$ around the point $\mathbf{x_{0}}$, that is

where $\mathbf{n}$ is the outer unit normal vector on $S$. Therefore, the divergence of a vector field is basically the volume averaged flux of the vector field out of $\Omega$; namely flux per unit volume. Note that by taking $\delta V\rightarrow 0$ we have already made the assumption that the vector field $\mathbf{F}$ is observable with infinite resolution. This assumption will be challenged later in this manuscript.

A well-known result from the vector calculus is the divergence theorem of Gauss reproduced below HildebrandFB:62a .

In order to extend the concept of the divergence of a vector field, discussed above, to the case of a finite volume, where the fluxes are calculated based on averaged quantities, we will adopt the definition of the divergence as defined in equation (1); that is a volume-averaged flux. We now follow the standard steps in the proof of the divergence theorem of Gauss with a more careful examination of the flux approximation at the faces of the volume. Before this, we first introduce the following definition.

As it will be clear shortly, this quantity represents the average flux of the vector field $\mathbf{F}$ out of a finite region with a representative length scale $\alpha$ defined by the averaging process $(\bar{.})$. For instance, for the Helmholtz filter $f=\overline{f}-\alpha^{2}\Delta\overline{f}$ where the representative scale of the finite averaging region is $\alpha$. The observation scale $\alpha$ could be imagined to be of similar size to the mesh size of the numerical simulation or the size of the measurement region (e.g., dimension of the wire in hotewire annemometry) in experiments. Note that, in a limit where the observable scale approaches zero this observable divergence reduces to the classical divergence defined in calculus.

We now prove the main result governing the divergence theorem for an observable conservation law.

We now start with the divergence definition in equation (1). The surface integral $\iint_{S}\mathbf{F\cdot n}\;dS$ over the rectangular parallelepiped can be regarded as a sum of six terms, one for each face of the parallelepiped. Now considering that $F_{1}=fu$ and using the average flux approximation at the surface $S_{1}$ and $S_{2}$, one can write

where we used the Taylor expansion and the definition of the Helmholtz filter to obtain an approximation for the flux on the face of the element²²2Since this is a rather important step here we reproduce the calculation of the values of a function at the edge of the parallelepiped based on the given values at the center of the volume. Using the Taylor expansion $$f(x+\frac{\Delta x}{2})\simeq f(x)+f^{\prime}(x)\dfrac{\Delta x}{2}+\mathcal{O}(\Delta x^{2})\vskip-2.84526pt$$ (4) Using the Helmholtz equation we obtain $$f(x+\frac{\Delta x}{2})\simeq\overline{f}(x)+f^{\prime}(x)\dfrac{\Delta x}{2}+\mathcal{O}(\Delta x^{2},\alpha^{2})\vskip-5.69054pt$$ (5) Notice that in the above relation we have assumed that the averaging process is performed on a scale representative of the volume $\Delta V$. Therefore, the averaging length $\alpha$ only decreases up to an observable scale. Beyond that we are not able to resolve smaller scales. At that point, while the observable scale remains constant, $\Delta x$ could independently approach a mathematically infinitesimal value as required in calculus.. Adding the contributions from the rest of the faces for the parallelepiped and taking the limit of the volume $\Delta V$ approaching zero, one can write

or

One should note that the limit of vanishing volume is performed independently and beyond what is observable with the averaging scale $\alpha$. That is the volume of the element continuously shrinks until it is not observable at the available numerical or experimental resultion. At such a limit there is a clear distinction between ${\rm odiv}\,\mathbf{F}$ and ${\rm div}\,\mathbf{F}$. At the mathematical limit where both the volume element $\Delta V$ and the length scale of the averaging kernel $\alpha$ approch zero independently, the difference between the ${\rm div}\,\mathbf{F}$ and ${\rm odiv}\,\mathbf{F}$ diminishes.

In this section we will develop the differential equations that must be satisfied by a fluid with an observability limit of $\alpha$. The equations are expected to resolve all relevant dynamical quantities up to this observable scale $\alpha$. As mentioned before, this scale is dictated by the resolution of a particular numerical simulation or the length scale resolution of a particular experimental equipment. For scales below $\alpha$ we assume that the medium still acts as a continuum and one can take the mathematical limit of $\Delta V\rightarrow 0.$ However, our observable scale remains at $\alpha$ even as this mathematical limit is taken to zero.

Before the differential form of the observable conservation laws for mass, momentum, and energy of a continuum are presented, we will formulate the general form of the observable differential conservation law for a quantity per unit volume, $f$. We assume $f$ satisfies a conservation law of the form

where $Q$ is the sum of all sources of $f$ inside the volume. For instance $Q$ is the total external force exerted on a control volume if one considers the conservation of momentum inside the volume. The rate of generation of $f$ is basically the outflow of $f$ minus the inflow at the boundaries of the volume plus the storage inside the volume. The rate of outflow minus the inflow at the boundaries was conveniently calculated in the previous section by the operator ${\rm odiv}\,f{\bf V}$ where ${\bf V}$ is the continuum velocity. Therefore, the observable conservation of $f$ reduces to the following differential equation

We now show that this equation satisfies the conservation of $f$ over a volume.

Now, the observable compressible Euler equations can be easily derived by using the conservation of mass ($\rho$), momentum ($\rho\mathbf{V}$), and energy ($\rho e_{T}=\rho(e+\frac{1}{2}V^{2})$) to obtain

where $\mathbf{\Sigma}$ is the sum of all external forces acting on the control volume or its surface and $S$ is the sum of all energy sources inside the control volume. Note that in the limit of $\alpha$ approaching zero, that is when the observable scale approaches a mathematical zero, one recovers the classical Euler equations.

In our previous publications Mohseni:10k , we presented some mathematical and numerical results for the existence and uniqueness of the solution in 1D Euler along with the convergence to the entropy solution when $\alpha\rightarrow 0$. We also showed that our 1D observable Euler equations support the same traveling shock solution as the Euler equations.

For the case of constant density, the observable continuity equation (11a) reduces to the classical case of a divergence free velocity field; that is ${\rm div}{\bf V}=0.$ In this case, the observable momentum equation (11b) can be simplified to

This is exactly the inviscid Leray equationLerayJ:34a , introduced without a formal derivation in 1934. A formal derivation of these equations from basic principles is offered here.

Several numerical experiments with increasing degrees of complexity were performed in order to numerically evaluate the performance of the observable Euler equations. These include benchmark test cases from Moin:10a using the Euler/Navier-Stokes equations; 1D Sod, 1D Shu-Osher, 2D shock-vorticity/entropy wave interaction, 3D Taylor-Green vortex, 3D isotropic compressible turbulence, and 3D shock-turbulence interaction. See Mohseni:10p ; Mohseni:10k for the details of our numerical experiments. In order to reduce potential numerical dissipation a pseudo-spectral discretization technique is used. The equations are advanced in time with a Runge-Kutta-Fehlberg predictor/corrector (RK45). Spatial derivatives and the inversion of the Helmholtz operator were computed in the Fourier domain. The terms were converted into the Fourier domain using a Fast Fourier Transform, multiplied by the appropriate term and then converted back into the physical domain.

Here we only present the results for the 2D interaction of a vorticity/entropy wave with a normal shock located at $x=1.5\pi$. Identical initial and boundary conditions and similar spatial resolutions as in Moin:10a is used here. Figure 2 shows the instantaneous vorticity contours and profiles for two different incident angles. The inflow vorticity/entropy wave interacts with the shock and changes its propagation direction and wavelength. In Figure 2(b) instantaneous vorticity profiles for different $\alpha$ are shown for the incident angle $\phi=75^{o}$. Compared to the solutions in Moin:10a , the vorticity profiles based on observable Euler equations do not show any spurious oscillations across the shock. Again note that this is done without any viscous terms. These results demonstrates the ability of the observable Euler equations to simulate 2D flows with shocks and resolve the shock without any viscous terms or side effects (numerical dissipation) of using one-sided schemes. Also note that the thickness of the shock is controlled by the parameter $\alpha$, the width of the averaging Krenel. For smaller $\alpha$, higher resolutions are required.

We expect to employ the observable divergence theorem in future to other field quantities in order to obtain regularized and observable evolution equations in other areas beyond the fluid dynamic problems considered here.