The spatial statistics of turbulent dissipation rates
James Glimm Vinay Mahadeo





   Department of Applied Mathematics and Statistics,

   Stony Brook University,

   Stony Brook, NY 11794
   


   Department of Applied Mathematics and Statistics,

   Stony Brook University,

   Stony Brook, NY 11794
   


We study the spatial statistics of velocity gradient volatility 
(i,e., the energy dissipation rate) in turbulent flow.
We extend 
the Kolmogorov-Obukhov theory but also narrow its scope. The models are
log normal, with verification from finely resolved large eddy and
direct numerical simulations. They are parameterized by
a mean and a covariance operator. Addressing applications to large eddy simulations,
the mean and the covariance depend on the resolved scale solution. 
Removing this resolved scale dependence by a locally defined rescaling of
the turbulent statistics yields a universal theory for
the subgrid statistics, specifically for the mean and variance of 
turbulent fluctuations, i.e., the log of the energy dissipation 
rate ϵ. 
The variance of the
velocity gradient statistics is found to be log normal, in accordance with
Kolmogorov 1962. Rescaling by the resolved scale mean and variance
removes the influence of the resolved scales on the subgrid
scales, justifying conceptually the universality we observe.
The universality is the basis of new power
scaling laws.
The new power laws, in turn, allow a simple parameterization of
the mean and covariance
of the subgrid rescaled dissipation rate, i.e., velocity gradient
volatility.
We treat the coarse grid
resolved space-time location as a random variable.
A restriction of the theory to a small range of unresolved but
stochastically modelled scales is proposed, with renormalization
group ideas offered to overcome this restriction.
 

The spatial statistics of turbulent dissipation rates
    Vinay Mahadeo
    December 30, 2023
=====================================================





§ INTRODUCTION

The scaling law of Kolmogorov <cit.>, although simple, is among the
deepest contributions to the understanding of turbulent flow. Turbulent flow is
generally considered to be the major outstanding problem of classical physics.
We propose new parameterized spatial statistics and
scaling laws for turbulent fluctuations which extend
ideas of the Kolmogorov 1962 (K62) theory but also narrow its scope. Verification is
from comparison to Direct Numerical Simulation (DNS).
Because we do not assume 
homogeneity or isotropy
of the flow,
the theory is applicable to the unresolved scales of
Large Eddy Simulations (LES), and is
a stochastic subgrid scale (SGS) model. A major contribution of the theory
presented here is to capture SGS fluctuations, in addition to the means
computed from traditional SGS theories.

We assume single fluid, constant density, incompressible flow, governed by the
Navier-Stokes equation.
The scalar dissipation rate is defined as

    ϵ = ν/2∑_i,j( ∂𝐮_i/∂𝐱_j +∂𝐮_j/∂𝐱_i)^2 
    = ν/2S+S^*_2^2  ,

with ν the dynamic viscosity, 𝐮 the fluid velocity, and

    S = ∂𝐮_i/∂𝐱_j

the strain rate.

Properties of ϵ are fundamental to a number of deep
theories of turbulence, including
notions of turbulent intensity, turbulent intermittancy <cit.>,
corrections to the Kolmogorov 5/3 exponent and scaling laws for
velocity moments. The analysis of ϵ has been applied to
turbulent diffusion, short distance asymptotics
for the velocity two point correlation function, and prediction of a fractal or
multifractal spatial distribution of regions of high turbulent intensity,
with the turbulent regions concentrated in a fractal set of dimension
D < 3 <cit.>. Scaling laws for wall bounded turbulence (a topic
not addressed here) are discussed in <cit.>. There is an extensive
literature on the above topics, which we do not attempt to document
in this article.
See the survey article <cit.> for elaborations and further references.
Kolmogorov <cit.> and Obukhov <cit.>
conjectured that ϵ obeys a lognormal distribution at sufficiently 
high Reynolds numbers.

Our main results are new scaling law for the 
mean and covariance of χ = ln(ϵ), 
which we designate as turbulent intensity.







The scaling laws originate
in a universality principle, of the mean μ and the
variance of χ for inertial
degrees of freedom. The connection between universality, which would
seem to require isotropic homogeneous turbulence, and subgrid scales
of LES, which can be very far from homogeneous or isotropic is one of
our central results. Depending on the local properties of the resolved
scales of turbulence, we introduce a rescaling of the turbulent statistics,
and it is only after this rescaling that the subgrid scales of turbulence
become universal. In fact, universality
does have limits, in terms of the
number of subgrid scales it can accommodate. In our numerical tests, we
consider 3 levels of mesh refinement, that is a change of length scales
by a factor of 8. For a large number of length scales, we have in mind
a renormalization group methodology, with successive integration over
smaller length scales, and a resetting of parameters using the
(statistical) knowledge of the already integrated scales, a point of view to
be elaborated in future work.


In Sec. <ref>, we specify the simulation study used
for verification of the theories proposed here. We define the
statistics and the ensemble used. We show that ϵ
has log normal statistics (and χ has normal statistics).

The spatial statistics for
ϵ and χ are characterized by the mean 
μ and covariance Σ^2
of χ. We compare a coarsely gridded 
LES to a finely gridded LES or DNS. The fine grid
extension of the coarse grid solution is of course unique and not stochastic,
but we add statistics in the dependence of this extension on the resolved
coarse grid cell from which it is derived. The statistics of the fine grid
extension for χ depends on the coarsely gridded resolved solution.

A detailed analysis of the χ statistics is
carried out in Fourier space,
in Sec. <ref>.
The analysis uses  the coarse grid resolved scales as
a random varible, to define the statistical properties of the subgrid
scales. 
Our premise is that the spatial aspects of this
dependence are largely captured by a resolved
scale variance Σ_0^2 and mean μ_0 for χ.
When the Σ_0^2 and μ_0
dependency is removed from the
refined solution extension by rescaling, the resulting statistics for χ,
i.e., Σ^2/Σ^2_0, is 
modeled as universal, with Σ^2/Σ^2_0 diagonal and
subject to simple modeling assumptions, including new power laws,
with only a few O(1) dimensionless parameters to govern 
its spatial statistics, all physically understandable.
The power laws result from a universality principal for the
rescaled statistics.
These approximations are confirmed by the fine grid comparison.

Sec. <ref>
extends this log normal
theory to individual components of the strain matrix S
and to the tensor dissipation rate.
A discussion of results and an outlook for future work is
given in Sec. <ref>. 



§ THE SPATIAL STATISTICS





 §.§ The simulation study


All verification tests are based on a series of coarse and finely gridded
simulations.
The finely gridded  LES/DNS study has over two decades of resolution,
with the coarsely gridded LES fraction of the
solution varying with the simulation
parameters. The coarsely gridded LES includes portions of the inertial
range and has up to a decade of unresolved inertial range scales
to which our modeling applies.

Our fine and coarse grid simulations are
conducted with the code <cit.>, using periodic boundary conditions
and forced turbulence in a domain size L^3, L = 0.05m, with the
viscosity (of air), 1.9e^-5m^2/s. The mesh parameters for three
simulations are given in Table <ref>. 
The resolved levels of turbulence in this numerical study
are not extreme.

Let Δ_c denote the coarse grid mesh size, and Δ_f
the fine grid mesh size, so that Δ_c = 8Δ_f.
As we are concerned with velocity gradients, which require two mesh cells
as a minimum stencil size for their evaluation, we also introduce
what we call the resolved coarse and fine grid mesh scales, with
Δ_r = 2Δ_c as the resolved coarse grid mesh and

    Δ_r = 2Δ_c = 8  (Resolved Δ_f) = 16Δ_f  .

Here the resolved fine grid level is 2Δ_f, needed to evaluate
the fine grid velocity gradients occurring in ϵ and χ.






 §.§ The statistical ensemble



As with turbulent flow in general, ϵ and χ 
can be regarded as random fields. In the LES context, 
ϵ and χ are known (deterministic) at the resolved scales,
while their subgrid scale fluctuations can be modeled using spatial statistics.
We write the fine grid resolved spatial mesh coordinate
i = i_r,i_s, where i_r is the coarse grid
resolved mesh portion of i. As the
resolved coarse mesh spacing is 2Δ_c, i_r spans all the coarse grid
LES cells other than the last or finest of them. i_s
is the remainder, occurring in the fine but not the coarse
LES grid resolution and also within the
finest level of the LES simulation. There are then n_r = 16^3 i_r values and
8^3 i_s values. Note that 2^3 Δ_f mesh indices are missing,
as we discuss resolved, not simulation mesh indices.

Our plan is to treat χ, regarded as a fiction of i_s, as random, 
depending on the random variable i_r.
The natural 
domain for its mean and covariance μ and Σ^2 is the
Hilbert space H of subgrid functions,
i.e., functions indexed by i_s but restricted to a single resolved 
coarse grid cell
and having mean zero there. μ is a vector in H and
Σ^2 is an i_s × i_s positive definite matrix
acting on H.

The statistics are collected as an average over all the possible n_r
values of i_r.



 §.§ The log normal property



We present numerical evidence to support the hypothesis that
the fine grid extension of the
resolved χ defines a Gaussian random field, 
with the mean and covariance defined from
averages of distinct resolved coarse grid cells.


The time dependent spatial average of χ has been tested for normality,
see <cit.> and references cited there, based on experimental and 
DNS evidence.



Our verification method to establish the Gaussian property
is also based on DNS evidence.
A multivariate random variable is Gaussian if
its inner product with any fixed vector is univariate Gaussian. Here we 
consider, as representative choices, the
fixed vectors defined by single DNS cell indices i_s. 

We use QQ plots to assess univariate normality. In these, the trial distribution
is transformed first by a change of the χ independent
variables to have zero mean and unit variance.
Then we apply an inverse Gaussian change of trial distribution
dependent variables, with the 
result compared to a unit slope straight line.

We create a single QQ plot by holding the fine grid mesh resolved cell indices
i_s fixed (16^3 choices), and using the multiple
i_r values from the 16^3 resolved coarse grid cells to form the ensemble. 
In Fig. <ref> 
we superimpose all QQ plots for all choices of i_s
as typical tests of normality,
demonstrating agreement to a Gaussian law up to ± 2 σ, 
a range of agreement consistent with the 4% accuracy expected
from the ensemble size.





 §.§ The mean and covariance



The coarse grid resolved level variance of χ is  basic to our rescaling
strategy. 
At the level of each coarse grid resolved grid cell, the variance
Σ_0^2 of χ(i_r) is defined as


    Σ_0^2 = (χ(i_r) - μ_0)^2

with the definitions 
μ_0 = χ(i_r) = n_r^-1∑_i_rχ(i_r).

Consider the fine grid values χ(i_r,i_s) for the resolved coarse
grid χ(i_r).
As a normal random variable, the statistics of χ is characterized
by its mean and its covariance. To obtain universal statistics,
we rescale χ by subtracting the resolved coarse grid mean
and dividing by the standard deviation
defining χ(i_r,i_s) = (χ(i_r,i_s) - μ(i_s))/Σ_0(i_r).
The mean μ of χ
is 


    μ(i_s) = n_r^-1∑_i_rχ(i_r,i_s)  .


The universal, rescaled, statistics are the statistics of
χ. 
We follow a similar definition for the covariance, Σ^2 of χ,


    Σ^2(i_s,j_s)
    = χ(i_r,i_s)χ(i_r,j_s) .


where the overbar denotes averaging over all resolved coarse grid cells, namely


    ( ·) = n_r^-1∑_i_r( ·)  .


Furthermore, we will normalize Σ^2 by Σ^2(i_r,i_r) = Σ^2(i_r) = Σ_0^2, the resolved scale variance, reducing the Reynolds number dependence of Σ.




§ SPATIAL STATISTICS IN FOURIER SPACE



We propose a Fourier space analysis for μ and Σ^2.
As these variables are independent of the resolved variables,
they are, in this approximation,  periodic within each resolved
cell. Thus their Fourier transforms depend on the
Fourier modes within a single resolved cell.

In preparing figures in this section, we present functions of 𝐤,
unless otherwise noted,
after binning all 𝐤 values into bins 
labeled by the scalar wave number
k with 200 bins whose size is linear in ln k.




 §.§ μ in Fourier space



In Fig. <ref>, we plot μ(k) vs. k in
ln ln scaled variables.
Observe the approximate scaling law μ(k) ∼ k^-d, with corrections coming from
the dissipation range of scales. Included in the plot is a model
curve of the form

    μ(k) ∼ k^-d

The coefficients in this model fit are tabuleted in
Table <ref>.












 §.§ Σ^2 in Fourier space




In Fig. <ref>, we plot the diagonal and off diagonal elements of
the covariance matrix 
Σ^2(𝐤,𝐥) for Reynolds number 2549.
The smooth curve in the upper part of the figure
represent the plot of the diagonal elements as a 
function of the frequency k, while the noisy data near
the bottom of the plot are the collection of all off diagonal elements,
plotted vs. (k+l)/2. The other Reynolds numbers give similar plots.
We conclude that the covariance is, within a good
approximation, diagonal.

Our modeling for Σ^2 (and for related operators describing the
strain matrix S and its individual matrix entries in 
Sec. <ref>) is based on a simple approximation, which appears
to be consistent with the DNS data as examined here. We note that the
two random variables contributing to the covariance Σ^2 are
orthogonal for 𝐤_𝐬≠𝐥_s and so the 
𝐤_s ≠𝐥_s terms vanish. In other words, the Fourier
coefficient of Σ^2(i_s,j_s) is Σ^2(𝐤_𝐬,𝐤_s),
which we write as Σ^2(𝐤_s).  In our DNS case,
Σ^2 is a 16^3 × 16^3 matrix before binning.

The Fourier coefficients for χ for distinct angular variables
𝐤/𝐤 are statistically independent. For distinct
values of the scalar
Fourier amplitude k = 𝐤,
we also assume statistical independence. These assumptions lead to a simple model for 
Σ^2, 
tractable for practical use.

We regard the wave number k as an observable on 
the random field χ, with the variance 
Σ^2(k) for
this observation. We derive k dependent scaling laws.
Approximately, the subspace of the Hilbert space in which
Σ^2 (k) is defined has dimension proportional to 4 π k^2,
namely the area of
a sphere in momentum space of radius k. Assuming equipartion of 
the variance of χ within this subspace, 
the inverse covariance is a multiple of the
identity there. We extend the variance equipartition hypothesis,
meaning a uniform level of variance,
to distinct k values. 
The subspace variance is then proportional to k^-2.

From a more probabilistic point of view, we regard the observation as the
result of asking O(k^2) questions, postulated to be independent:
what is the variance of
χ associated with each 𝐤 lying on the sphere of radius k?
Combining the variance of O(k^2) 
independent quantities gives the same answer. 

Volume scaling laws apply to the volume dependent coefficients in
the scaling law for Σ^2, and use a similar reasoning.
Thus, again assuming
equipartition, the mesh cell averages over a single coarse grid resoled
cell scale with Δ_r^3.
In summary, we propose scaling laws 
for Σ^2(k)dk in terms of
the scalar wave number k and spatial 
volume, with the dispersion relation

    Σ^2(k)/Σ_0^2 = 
    ( Δ_r^3 4 π C_2^2 k^2 exp(C_4^2 η^2 k^2) )^-1

with dimensionless constants C_2, C_4 to be determined. 
The factors Δ_r^3 and 4π k^2dk are volumes for physical space
and Fourier space. The C_4 term,
not previously discussed, is postulated to capture the leading order
viscous corrections to the otherwise purely inertial range theory.





In Fig. <ref>, we plot
Σ^2(k) vs. k
on ln ln scales with the
model curves (<ref>) superimposed, based on optimal choices
of the coefficients C_i. We first fit C_2 in the inertial range
(small k values), and then fit C_4 over a 1/2 decade at the end of the
interial range. We observe that most of the Reynolds number
dependence of the C_i 
has been removed by the scaling introduced here, see 
Table <ref>. 



§ A LOG NORMAL MODEL FOR VELOCITY GRADIENTS





We repeat the analysis of Secs. <ref> 
– <ref> for the strain matrix restricted to subspaces. Let E_1
and E_2 be self adjoint projection operators on R^3, and consider the
reduced strain matrix E_1 SE_2 = S_1,2.
In more detail, we consider 

    ϵ_1,2 =
    ( ν / 2 )  E_1∂𝐮_i/∂ x_j E_2 _2^2

and the associated 
χ_1,2 = ln (ϵ_1,2 ).
Thus S_1,2 is the strain associated with fluid velocities in directions
spaned by E_1 and fluctuations occurring with spatial variation in 
directions spanned by E_2.



With μ_1,2 the rescaled resolved grid level mean for 
the rescaled χ_1,2, 
Σ^2_1,2 is the rescaled covariance 

    Σ_1,2^2(i_s,j_s) =
        n_r^-1∑_i_r (χ_1,2(i_r,i_s) -μ_1,2(i_r)) 
    (χ_1,2(i_r,j_s) - μ_0,1,2(i_r))  ,

With these changes, we repeat the analysis of 
Figs. <ref>–<ref>
with the modified dispersion relation

    Σ^2_1,2(k)= 
    
        (dim E_1 dim E_2Δ_r^3 4 π C_2,1,2^2 k^2 exp( C_4,1,2^2 η^2 k^2 ) )^-1 ,

and modified coefficients C_2,1,2 and C_4,1,2. 



We assume that E_1 and E_2 are both one dimensional.



We also show that the tensor dissipation operators

    ϵ' = (ν/2) S^*S  ,     ϵ” = (ν/2)SS^*

satisfy log normal statistics. Here the exponential and the log
in the definition of log normal 
(χ' = log(ϵ') and
χ” = log(ϵ”)) 
is considered in the sense of matrix operations.
See Figs. <ref>,
<ref> and <ref>.

We note that the matrix QQ plots, and to some extent, all the QQ plots, raise
the possibility of nonGaussian statistics as a correction to the normal and
log normal statistics assumed here.
Especially in Fig. <ref>,
systematic corrections to the Gaussian property show smaller positive
excursions and stronger negative excursions for χ
than Gaussian for the diagonal
matrix elements. Both the negative and positive differential
excursions reflect a smaller ϵ and
a possible bias toward a laminar rather than turbulent flow. In contrast, the
off diagonal QQ plots suggest stronger χ excursions than Gaussian
for both positive and negative excursions,
a phenomena known as heavy tails in the statistical data analysis literature.
These are the matrix elements which occur in vortex stretching and folding,
commonly believed to be an important element of fully developed turbulence.







§ DISCUSSION



Verification by computational studies will typically encounter only
a limited range of LES unresolved scales, compared to a finely gridded
simulation. But the previous analysis, when applied over many
scales, has a conceptual flaw. If there is no resetting of the
Σ_0^2 parameters within the statistics model of a single
resolved scale, then we could apply this theory in the absence of
the coarsely gridded LES to the entire computation. Doing this would
remove the rescaling already seen to be important. For this reason, we
restrict the model to a small range of unresolved scales. In 
a future publication 
we reformulate the model in a renormalization group
manner in the case of multiple scales.

In an LES context, we analyze the unresolved scales statistically, with 
an emphasis on the strain rate ϵ and its logarithm 
lnϵ = χ. Resolved scale adjustments to the mean
and variance for χ lead to a theory for the resulting covariance Σ^2.
We find that Σ^2 is approximately diagonal in Fourier space.
An equipartition hypothesis for the variance of χ 
leads to a new power scaling law and 
a simple parameterization for Σ^2, verified by comparison to DNS
over about a decade of LES unresolved inertial scales. For reasons discussed,
we restrict this theory to a single or limited range within the inertial
range of scales.
This pure dissipation rate model is extended to the individual velocity
gradient components
ϵ_1,2 = (ν / 2) (∂u_i / ∂ x_j)^2
and to the tensor dissipation rates S^*S and SS^*.

Open questions, to be addressed in future work, concern temporal statistics of
turbulent intensity, statistical models for the velocity gradients and 
strongly nonhomogeneous flows, such as boundary layers.

plain


