Highly turbulent solutions of LANS-alpha and their LES potential
Highly turbulent solutions of LANS−α and their LES potential
Jonathan Pietarila Graham,1, 2 Darryl D. Holm,3, 4 Pablo D. Mininni,1, 5 and Annick Pouquet1
1National Center for Atmospheric Research,
∗ P.O. Box 3000, Boulder, Colorado 80307, USA
2currently at Max-Planck-Institut für Sonnensystemforschung,
37191 Katlenburg-Lindau, Germany
3Department of Mathematics, Imperial College London, London SW7 2AZ, UK
4Computer and Computational Science Division,
Los Alamos National Laboratory, Los Alamos, NM 87545, USA
5Departamento de Fı́sica, Facultad de Ciencias Exactas y Naturales,
Universidad de Buenos Aires, Ciudad Universitaria, 1428 Buenos Aires, Argentina
(Dated: March 2, 2022)
∗ The National Center for Atmospheric Research is sponsored by the National Science Foundation
http://arxiv.org/abs/0704.1928v3
Abstract
We compute solutions of the Lagrangian-Averaged Navier-Stokes α−model (LANS−α) for significantly
higher Reynolds numbers (up to Re ≈ 8300) than have previously been accomplished. This allows suf-
ficient separation of scales to observe a Navier-Stokes inertial range followed by a second inertial range
specific to LANS−α. Both fully helical and non-helical flows are examined, up to Reynolds numbers of
∼ 1300. The analysis of the third-order structure function scaling supports the predicted l3 scaling; it cor-
responds to a k−1 scaling of the energy spectrum for scales smaller than α. The energy spectrum itself
shows a different scaling which goes as k1. This latter spectrum is consistent with the absence of stretching
in the sub-filter scales due to the Taylor frozen-in hypothesis employed as a closure in the derivation of
LANS−α. These two scalings are conjectured to coexist in different spatial portions of the flow. The l3
(E(k) ∼ k−1) scaling is subdominant to k1 in the energy spectrum, but the l3 scaling is responsible for
the direct energy cascade, as no cascade can result from motions with no internal degrees of freedom. We
demonstrate verification of the prediction for the size of the LANS−α attractor resulting from this scaling.
From this, we give a methodology either for arriving at grid-independent solutions for LANS−α, or for
obtaining a formulation of a Large Eddy Simulation (LES) optimal in the context of the alpha models. The
fully-converged grid-independent LANS−α may not be the best approximation to a direct numerical sim-
ulation of the Navier-Stokes equations since the minimum error is a balance between truncation errors and
the approximation error due to using LANS−α instead of the primitive equations. Furthermore, the small-
scale behavior of LANS−α contributes to a reduction of flux at constant energy, leading to a shallower
energy spectrum for large α. These small-scale features, however, do not preclude LANS−α to reproduce
correctly the intermittency properties of the high Reynolds number flow.
PACS numbers: 47.27.ep; 47.27.E-; 47.27.Jv; 47.50.-d
I. INTRODUCTION
Since the degrees of freedom for high Reynolds number (Re) turbulence, as can be encountered
in geophysical and astrophysical flows, can be very large, the implementation of their numerical
modeling can easily exceed technological limits for computations. Furthermore, since truncation
of the omitted scales removes important physics, e.g., of multi-scale interactions, the only ap-
proach to a numerical study of such flows is to employ subgrid modeling of those scales. This
is frequently accomplished with Large Eddy Simulations (LES–see [1, 2, 3] for recent reviews).
This is of importance for geophysical, astrophysical and engineering applications and can have
consequences for meteorological [4] and climate prediction simulations [5], for instance. While
realistic Reynolds numbers will remain out of reach for the foreseeable future, subgrid modeling
can be an extremely useful tool in the computation of simulations for such applications.
The incompressible Lagrangian-averaged Navier-Stokes equations (LANS−α, α−model, or
also the viscous Camassa-Holm equation) [6, 7, 8, 9, 10, 11] is one possible subgrid model. It can
be derived, for instance, by temporal averaging applied to Hamilton’s principle (where Taylor’s
frozen-in turbulence hypothesis is applied as the closure, and also as the only approximation of the
derivation) [12, 13, 14]. For this reason, the momentum-conservation structure of the equations
are retained. For scales smaller than the filter width, LANS−α reduces the steepness of steep
gradients of the Lagrangian mean velocity and limits how thin vortex tubes become as they are
transported (the effect on larger length scales is negligible) [9]. The α−model may also be derived
from smoothing the transport velocity of a material loop in Kelvin’s circulation theorem [11]. Con-
sequently, there is no attenuation of resolved circulation, which is important for many engineering
and geophysical flows where accurate prediction of circulation is highly desirable. LANS−α
has previously been compared to direct numerical simulations (DNS) of the Navier-Stokes equa-
tions at modest Taylor Reynolds numbers (Rλ ≈ 72 [15], Rλ ≈ 130 [9], and Rλ ≈ 300 [16]).
LANS−α was compared to a dynamic eddy-viscosity LES in 3D isotropic turbulence under two
different forcing functions (Rλ ≈ 80 and 115) and for decaying turbulence with initial conditions
peaked at a low wavenumber (Rλ ≈ 70) and at a moderate wavenumber (Rλ ≈ 220) [17]. In
these comparisons, LANS−α was preferable in that it demonstrated correct alignment between
eigenvectors of the subgrid stress tensor and the eigenvectors of the resolved stress tensor and
vorticity vector. LANS−α and a related regularization, the Leray model, were contrasted with a
dynamic mixed (similarity plus eddy-viscosity) model in a turbulent mixing shear layer (Re ≈ 50)
[18, 19]. LANS−α, with relatively high subfilter resolutions, was the most accurate of these three
LES tested at this moderate Re, but it was found that the effects of numerical contamination can
be strong enough to lose most of this potential. This could pose some limitations on its practical
use. Quantifying those limitations is one of the goals of this present work. We will also find in this
study that, even with sufficient subfilter resolution, LANS−α fails to represent all the neglected
physics in a more turbulent regime (higher Re).
The α−model also describes an incompressible second-grade non-Newtonian fluid (under a
modified dissipation) [11]. In this interpretation, α is a material parameter which measures the
elastic response of the fluid. Either from this standpoint, from its status as a regularization of
the Navier-Stokes equations, or, independently of any physically motivation, as a set of partial
differential equations with proven unique regular solutions, we may analyze LANS−α without
any LES considerations. Analyzing inertial-range scaling for LANS−α for moderate and large α,
as well as identifying different scalings at scales larger and smaller than α is another of the goals
of this work. In this context we also study the numerical resolution requirements to obtain well-
resolved solutions of LANS−α (i.e., grid-independent solutions) which leads to a verification of
the predictions of the size of the attractor in LANS−α [11, 20]. Section II presents the LANS−α
model, our numerical experiments and technique. In Section III we analyze inertial-range scaling
for LANS−α. In Section IV we determine the numerical resolution requirements to obtain well-
resolved solutions of LANS−α. In Section V we address the LES potential of LANS−α by
comparing α−model simulations to a 2563 DNS (Re ≈ 500, Rλ ≈ 300), a 5123 DNS (Re ≈ 670,
Rλ ≈ 350), a 5123 DNS (Re ≈ 1300, Rλ ≈ 490), a 10243 DNS (Re ≈ 3300, Rλ ≈ 790), and a
20483 DNS (Re ≈ 8300, Rλ ≈ 1300). (The Re ≈ 3300 simulation has been previously described
in a study of the imprint of large-scale flows on local energy transfer [21, 22].) In Section VI,
we compare and contrast in more detail LANS−α solutions with DNS at Re ≈ 3300. Finally, in
Section VII we summarize our results, present our conclusion, and propose future directions of
investigation.
II. TECHNIQUE
We consider the incompressible Navier-Stokes equations for a fluid with constant density,
∂tvi + vj∂jvi = −∂ip+ ν∂jjvi + Fi
∂ivi = 0, (1)
where vi denotes the component of the velocity field in the xi direction, p the pressure divided
by the density, ν the kinematic viscosity, and Fi an external force that drives the turbulence (in
all results, the time, t, is expressed in units of the eddy-turnover time). The LANS−α equations
[6, 7, 8, 9, 10, 11] are given by
∂tvi + uj∂jvi + vj∂iuj = −∂iπ + ν∂jjvi + Fi
∂ivi = ∂iui = 0, (2)
where ui denotes the filtered component of the velocity field and π the modified pressure. Filtering
is accomplished by the application of a normalized convolution filter L : f 7→ f̄ where f is any
scalar or vector field. By convention, we define ui ≡ v̄i. We choose as our filter the inverse of a
Helmholtz operator, L = H−1 = (1 − α2∂kk)−1. Therefore, u = gα ⊗ v where gα is the Green’s
function for the Helmholtz operator, gα(r) = exp(−r/α)/(4πα2r) (i.e., the well-known Yukawa
potential), or in Fourier space, û(k) = v̂(k)/(1 + α2k2).
We solve Eqs. (1) and (2) using a parallel pseudospectral code [23, 24] in a three-dimensional
(3D) cube with periodic boundary conditions. In most of the runs, we employ a Taylor-Green
forcing [25],
sin k0x cos k0y cos k0z
− cos k0x sin k0y cos k0z
(generally, with k0 = 2), and employ dynamic control [26] to maintain a nearly constant energy
with time. This expression Eq. (3) is not a solution of the Euler’s equations, and as a result small
scales are generated fast when the fluid is stirred with this forcing. The resulting flow models
the fluid between counter-rotating cylinders [27] and has been widely used to study turbulence,
including studies in the context of the generation of magnetic fields through dynamo instability
[28]. We also consider some runs with random and ABC [22] forcing. We define the Taylor mi-
croscale as λ = 2π
〈v2〉/〈ω2〉, and the mean velocity fluctuation as vrms =
E(k)dk
The Taylor microscale Reynolds number is defined by Rλ = vrmsλ/ν and the Reynolds number
based on a unit length is Re = vrms × 1/ν.
III. INERTIAL RANGE SCALING OF LANS−α
A. l3 scaling of third-order structure function derived from the Kármán-Howarth theorem for
LANS−α
For LANS−α, the H1α(u) norm is the quadratic invariant to be identified with the energy,
= −2νΩα, (4)
where
(u− α2∇2u) · ud3x = 1
v · ud3x, (5)
ω · ω̄d3x. (6)
As usual, we define the (omni-directional) spectral energy density, Eα(k), from the relation
Eα(k)dσdk =
Eα(k)dk (7)
where
dσ represents integration over the surface of a sphere. The α−model possesses a theorem
corresponding to the Kármán-Howarth theorem [29] for the Navier-Stokes equations and, as in the
Navier-Stokes case, scaling of the inertial range energy spectra may be derived from it [30]. We
summarize here the dimensional analysis argument for the LANS−α inertial range scaling that
follows from this theorem, beginning from Equation (3.8) in Ref. [30]. We use the short notation
vi ≡ vi(x), u
i ≡ u
, t) and r ≡ x′ − x. In the statistically isotropic and homogeneous case,
without external forces and with ν = 0, taking the dot product of Eq. (2) with u
j we can obtain
the equation
∂tQij =
T mij − α2Smij
. (8)
The trace of this equation is the Fourier transform of the detailed energy balance for LANS−α.
Qij =
j + vju
is the second-order correlation tensor while
T mij =
j + vju
i + v
iuj + v
jui)u
, (10)
Smij =
(∂mul∂iul)u
j + (∂mul∂jul)u
i + (gα ⊗ τ ′
j )vi + (gα ⊗ τ ′
i )vj
, (11)
are the third-order correlation tensors for LANS−α and τ ji is the sub-filter scale stress tensor.
For α = 0 this reduces to the well-known relation derived by Kármán and Howarth. The energy
dissipation rate for LANS−α , εα, satisfies εα ∝ ∂tQij . By dimensional analysis in Eq. (8) we
arrive at
(vu2 +
u3). (12)
For large scales such that l ≫ α, the second right hand term is ignored, u ≈ v, εα ≈ ε, and we
arrive at the scaling of the four-fifths law, < (δv‖(l))
3 >∼ εl [31]. Here, δv‖(l) ≡ [v(x+ l) −
v(x)] · l/l is the longitudinal increment of v. The four-fifths law expresses that the third-order
longitudinal structure function of v, Sv3 ≡ 〈(δv‖)3〉, is given in the inertial range in terms of the
mean energy dissipation per unit mass ε by
Sv3 = −
εl, (13)
or, equivalently, that the flux of energy across scales in the inertial range is constant. We also obtain
the Kolmogorov 1941 [32, 33, 34] (hereafter, K41) energy spectrum, E(k)k ∼ v2 ∼ ε2/3l2/3, or,
equivalently,
E(k) ∼ ε2/3k−5/3. (14)
For small scales such that l ≪ α, however, v ∼ α2l−2u and both right hand terms are equivalent
in Eq. (12), and our scaling law becomes
Su3 ≡< (δu‖(l))3 >∼ εαα−2l3. (15)
Note that this scaling differs in a substantial way from the Kolmogorov scaling (∼ l). For our
small scale energy spectrum we then have
Eα(k)k ∼ uv ∼ ε2/3α α2/3, (16)
where we used u ∼ α−2l2v. The energy spectrum for scales smaller than α is then
Eα(k) ∼ ε2/3α α2/3k−1. (17)
This spectrum can also be derived from phenomenological arguments originally introduced by
Kraichnan [35], and it differs from the Navier-Stokes spectrum due to the fact that the fluid is
advected by the smoothed velocity u which does not directly correspond to the conserved energy
Eα [11].
FIG. 1: Third-order longitudinal structure function of the smoothed velocity field u, Su3 , versus l for large
α LANS−α (α = 2π/3 indicated by the vertical dotted line). The scales identified with an inertial range
are marked by vertical dashed lines and the scaling predicted by Eq. (15), l3, is indicated by a solid line.
The fitted scaling exponent ζu3 (S
3 (l) ∼ lζ
3 ) is found to be ζu3 = 2.39 ± .04. This is more consistent with
the scaling given by Eq. (15) than K41 scaling, l1 Eq. (13), or other proposed LANS−α scalings (indicated
by dotted lines, see text).
We test this prediction for LANS−α scaling at a resolution of 2563 (ν = 1.2×10−4) by moving
both the forcing (k0 = 1) and α (kα ≡ 2π/α = 3) to large scales in order to increase the number
of resolved scales for which kα > 1. In so doing, we are assuming that the scaling for large α is
the same as for small α and large k (for evidence to this effect, see [36]). Confirmation as given
by Eq. (15) is presented in Fig. 1 where we plot Su3 as a function of l (by convention, we plot
Su3 =< |δu‖(l)|3 > to reduce cancellation in the statistics). The scales identified with an inertial
range k ∈ [6, 10] are marked by vertical dashed lines and the predicted scaling, l3, is indicated by
a solid line. We fit a scaling exponent (Su3 (l) ∼ lζ
3 ) and find ζu3 = 2.39± .04. This is significantly
steeper than the classical Kolmogorov scaling given by Eq. (13); it can thus be viewed as more
consistent with the scaling given by Eq. (15). It is also more consistent with l3 than with other
possible LANS−α scalings: under the assumption that the turnover time scale of eddies of size ∼ l
is determined by the unsmoothed velocity v, we find Su3 (l) ∼ l5, and if it is determined by
v · u,
we find Su3 (l) ∼ l4 (see, e.g., Refs. [16, 36, 37, 38]). The observed scaling corresponds to none
of these cases, and is actually closer to an evaluation of the turnover time tl at the scale l given by
tl ∼ l/ul (with Su3 (l) ∼ l3). Note that for 2D LANS−α, however, it is the case that the scaling
is determined by the unsmoothed velocity v [36]. We note that this is one of many differences
between the 2D and 3D cases (e.g., ideal invariants and cascades). Another difference, which we
shall show in Section VI, is that in 2D vorticity structures decrease in scale as α increases while
in 3D there is a change in aspect ratio with structures getting both shorter and fatter. This may, in
fact, be related to the shallower LANS−α energy spectrum for kα > 1 which we show in Section
VI. While differences are observed between the scaling shown in Fig. 1 and Eq. (15), the error
bars deny a K41 scaling (as well as the l4 and l5 scalings) at scales smaller than α. We believe the
discrepancy between the observed and predicted scaling can be due to lack of resolution to resolve
properly the inertial range at sub-filter scales. We have less than a decade of inertial range and
only 2563 points for the statistics. As more computational resources become available, this scaling
should be re-examined.
B. Subdominance of the k−1 energy spectrum and rigid-body motions
As a consequence of LANS−α’s Taylor’s frozen-in hypothesis closure, scales smaller than α
can phase-lock into coherent structures and be swept along by the larger scales (see, e.g., [30]).
If we assume, formally, that this “frozen-in turbulence” takes the form of “rigid bodies” in the
smoothed velocity field (no stretching), we arrive at a much different spectrum than k−1, Eq.
(17). All scales smaller than α are subject to the frozen-in hypothesis and we expect to find such
rigid bodies at these scales. We note that collections of “rigid” portions of the flow (rotating or
non-rotating) reduce the total degrees of freedom (dof) and make physical sense with LANS−α’s
relation to second-grade fluids: these rigid bodies can be envisioned as polymerized portions of
FIG. 2: Spectral energy density, E(k), versus wavenumber, k, for large−α LANS−α solution. Here forcing
(k0 = 1) and α (kα ≡ 2π/α = 3, vertical dotted line) are set at the largest scales to increase the number of
scales for which kα > 1. Spectra are plotted for three norms: H1α(u) norm (solid line), L
2(u) norm (dotted
line), and the L2(v) norm (dashed line). As these last two norms are not quadratic invariants of LANS−α,
we employ the H1α norm for all following results. All three spectra correspond to that derived from the
assumption of rigid bodies in the smoothed velocity u, Eq. (19). The vertical dashed lines are at the same
scales as those in Fig. 1.
the fluid. As a matter of fact, in such structures all internal dof are frozen. These “rigid bodies”
follow as well from the consideration of LANS−α as an initial value problem in Fourier space,
for which we have û(k) = v̂(k)/(1+α2k2). In the limit as α approaches infinity, all wavenumber
(and spatial) dependence for v̄ is eliminated and the entire flow is advected by a uniform velocity
field (advection without internal degrees of freedom).
For a rigid body there can be no stretching and, therefore, all the longitudinal velocity incre-
ments, δu‖, must be identically zero (δu(l) = Ω × l from basic mechanics with Ω the rotation
vector and, hence, δu‖(l) = δu(l) · l/l = 0). Note that in LANS−α Eq. (2) the vj∂iuj term
contributes only a rotation and not a stretching of u. Such polymerization would have two con-
sequences. Firstly, since there is no stretching, these rigid bodies would not contribute to the
turbulent energy cascade,
< (δu‖(l))
3 >= 0. (18)
Secondly, the energy spectrum from dimensional analysis (u2 ∼ const, for large α/l: u = (1 +
α2/l2)−1v ∼ l2v, and Eα(k)k ∼ uv ∼ k2) is
Eα(k) ∼ k. (19)
This is, in fact, the observed LANS−α spectrum for kα ≫ 1 as is shown in Fig. 2. We verified
that the spectrum is not the result of under-resolved runs, as is the case, e.g., in the k2 spectrum
observed in truncated Euler systems [39] or in extremely under-resolved spectral simulations of
the Navier-Stokes equations. Indeed, equipartition of the energy among all modes in a truncated
Euler−α system should also lead to a k2 spectrum. Along with several experiments with different
viscosities and also with statistically homogeneous and isotropic forcing (not shown here), these
are assurances that the observed spectrum is not a result of inadequate numerical resolution. It
should be noted that this is the same computation for which the third-order structure function is
shown in Fig. 1. The third-order structure function is consistent with a l3 scaling (corresponding
to a k−1 energy spectrum) while the spectrum itself is k1. (Also shown in Fig. 2 are the L2(u) ≡
〈u2〉/2 and the L2(v) ≡ 〈v2〉/2 norms which (through u ∼ α2v/k2 for kα ≫ 1) correspond
to k−1 and k3 spectra, respectively. Since the analytical properties of the LANS−α solution are
based on the energy balance, dEα/dt = −2νΩα, in the H1α(u) norm, we employ this norm for all
following results.) These two different scalings, l3 and k1, are consistent with a picture where a
fluid has both rigid-body portions at scales smaller than α (wherein there is no turbulent cascade)
and spatial regions between these where the cascade does take place. For the structure functions,
a non-cascading rigid body does not contribute to the scaling and consequently the cascading
contribution, Eq. (15), dominates. The energy spectrum, however, for the limit of k very large, is
dominated by the k+1 term, and hence the k−1 component is subdominant.
We further explore the validity of this picture by examining the spatial variation of the cubed
longitudinal increment, (δv‖(l))
3 in DNS, and (δu‖(l))
3 in LANS−α for α/l ≫ 1, which in each
case is proportional to the energy flux across a fixed scale l. (The presence of the hypothesized
“rigid bodies” should be evident as significant portions of the flow where there is no energy flux.)
In Fig. 3 we show visualizations of these quantities corresponding to l = 2π/10 (k = 10) for
both the large-α LANS−α simulation and a highly turbulent DNS (k0 = 2, ν = 3 × 10−4). The
scale (k = 10) is chosen as it is in the inertial ranges of both flows. We note that for LANS−α,
a significant portion of the flow is not contributing to the flux of energy to smaller scales (the
filling factor for (δu‖(2π/10))
3 < 10−2 is 0.67 as compared to 0.26 for the Navier-Stokes case).
These regions can be identified as “polymerized” or “rigid bodies” in u and their locations are
found to be robust when the l used for (δu‖(l))
3 is varied over a factor of 2. Moreover, this is
highlighted in the probability distribution functions (pdfs), see Fig. 4, where we see the LANS−α
pdf is more strongly concentrated around zero than the DNS. This is consistent with the idea that
FIG. 3: Two-dimensional slice of the cubed longitudinal increment (δu‖(2π/10))
3for LANS−α and
(δv‖(2π/10))
3 for DNS. For all black pixels, the cubed longitudinal increment is less than 10−2 (approxi-
mately consistent with rigid bodies). On the top is the large-α simulation (k0 = 1, kα = 3, ν = 1.2×10−4)
where the filling factor (computed over the entire 3D domain) is 0.67. On the bottom is a DNS of Navier-
Stokes (k0 = 2, ν = 3 × 10−4) where the filling factor is 0.26. Thus, a much greater portion of the flow is
consistent with collections of rigid bodies for the large−α simulation.
FIG. 4: Pdfs of (δv‖(2π/10))
3 for DNS (N = 1024, solid line), and of (δu‖(2π/10))
3 for LANS−α
(N = 256, dashed line), and of the DNS downgraded to lower resolution (N = 256, dotted line). See Fig.
3 for simulation parameters. Note that both pdfs have a slight positive asymmetry consistent with a positive
dissipation rate ε(α). The LANS−α pdf is more strongly concentrated around zero consistent with the idea
that portions of the flow (at scales smaller than α) are acting as rigid bodies.
the internal dof of large portions of the flow (at scales smaller than α) are frozen. We point out that
this comparison is not a LES validation, but, rather, a comparison between the dynamics of two
different fluids at similar Reynolds numbers. One flow is a well-resolved numerical solution of the
Navier-Stokes equations, and the other is a well-resolved solution of the LANS−α equations with
large α. For this reason a reduced resolution (N = 256) representation for the DNS (for which
N = 1024) is not depicted in Fig. 3. We have performed such a down-sampling, however, and
find the filling factor is reduced even more, to 0.14, and the tails of the pdf increase over the full-
resolution analysis (dotted line in Fig. 4). No inverse Helmholtz filtering, H−1 is applied to the
DNS data. Note that this would amount to computing (δu‖(l))
3 in the DNS, which has no meaning
in the dynamics of the Navier-Stokes equations (the energy flux is proportional to (δv‖(l))
We end this section with further evidence of coexistent energy spectra, k−1 and k1, in separate
spatial portions of the flow. We mask out all portions of the flow that we identify with rigid
bodies ((δu‖(2π/10))
3 < 10−2, a 2D slice of which is shown in Fig. 3). The energy spectrum
of the remaining portion of the flow is shown in Fig. 5 as a dashed line to be compared with the
spectrum of the entire flow shown as a solid line. The operation of spatially filtering the flow
before computing the spectrum serves to “smear out” the energy spectrum by convolving it with
the spectrum of the filter. Deconvolution in 3D with N = 256 is intractable and we are, therefore,
unable to remove this “smearing” of the energy spectrum of the cascading portions of the flow.
FIG. 5: Spectral energy density, E(k), versus wavenumber, k, for large−α LANS−α solution. The solid
line indicates the spectrum as given in Fig. 2 but for a single snapshot (the same as selected for Fig. 3).
The dashed line indicates the spectrum wherein all portions of the flow associated with “rigid bodies” (a
2D slice of which is shown in Fig. 3) are removed. This provides further evidence that the flow spatially in
between the “rigid bodies” possesses a negative power law energy spectrum (the predicted k−1 power law
is shown as a solid line).
Nonetheless, after conducting what tests we could with the filtering process (not shown here),
we conclude that the power law of the energy spectrum of these portions is negative and, thus,
distinctly different from that of the rigid bodies.
IV. RESOLUTION REQUIREMENTS FOR GRID-INDEPENDENT LANS−α SOLUTIONS:
SIZE OF ATTRACTOR
It is useful to make a distinction between the quality of a subgrid model and effects arising from
nonlinear interactions with discretization errors at marginal spatial resolutions (which are more
characteristic of the discretization employed than of the subgrid model) [19, 40, 41]. Before doing
this, we require an estimate for the total degrees of freedom for the LANS−α attractor which as
we show, unlike for the 2D case (see [36]), for the 3D case is reduced compared to Navier-Stokes.
The subdominant l3 scaling is associated with the flux of energy to small scales and thus must
be used to estimate the degrees of freedom of the LANS−α attractor, dofα. For dissipation the
large wavenumbers dominate and, therefore, combining the LANS−α energy balance, Eq. (4),
with its sub-filter scale energy spectrum, Eq. (17), allows us to implicitly specify its dissipation
wavenumber, kαη , by
∫ kαη
k2Eα(k)dk ∼
∫ kαη
2/3α2/3k−1dk ∼ εα2/3α2/3(kαη )2. (20)
Then we have,
kαη ∼
ν1/2α1/3
. (21)
Using that the linear numerical resolution, N , must be proportional to the dissipation wavenumber
(N ≥ 3kαη ) and that Re ∼ ν−1, we arrive at
N = C0kα
1/3Re1/2, (22)
or, equivalently,
dofα =
Re3/2, (23)
where C0 is an unknown constant (for further details see [11]). We verify this prediction and
determine the constant C0 through the use of a database stemming from studies in which both the
free parameter, α (or, equivalently, kα) and the linear resolution, N , for a set of DNS flows with
Re ≈ 500, 670, 1300, and 3300 are varied. In so doing, we establish the necessary numerical
resolution for convergence to a grid-independent solution.
Convergence to the grid-independent solution is determined by comparison of the energy spec-
trum, Eα(k), between runs with a constant filter and varying resolution. In Fig. 6(a), we make
such a comparison for Re ≈ 500 (N = 256 for DNS) and kα = 14 (N = 84, 96, 108, 128, and
192 for LANS−α). We plot energy spectra compensated by k5/3 so that a K41 k−5/3 spectrum
would be flat. We see, based on comparing the energy spectra at wavenumbers smaller than kα
to the 1923 LANS−α spectrum, that simulations at resolutions of 963 and less are not converged
while the one at 1283 is. That is, except for the very small scales at the end of the dissipative range,
there is very little difference between the spectra at 1283 and at 1923 (i.e., the solution is “grid-
independent”). Meanwhile, for resolutions of 963 and less the spectra vary greatly with resolution
(i.e., they are “unresolved”). In Fig. 6(b), we collect all the results of similar studies (Re ≈ 500) in
a plot of resolution, N , versus inverse filter width, kα. (We change N for a given α, then change α
and iterate.) Pluses correspond to grid-independent solutions, X’s to under-resolved solutions, and
squares to “undecided” runs (i.e., that are neither clearly resolved nor clearly under-resolved). The
FIG. 6: (Color online.) Plots for Re ≈ 500 simulations demonstrating convergence to the grid-independent
LANS−α solution. (a) Average energy spectra (t ∈ [20, 33], t is time in units of eddy turn-over time)
compensated by K41 for LANS−α simulations, kα = 14: 1923 (black solid), 843 (red dotted), 963 (green
dashed), 1083 (blue dash-dotted), and 1283 (pink dash-triple-dot). The vertical dashed line denotes kα.
Inset is a blow-up near kα where convergence can be clearly seen. LANS−α at a linear resolution of
1283 is approximately converged to the grid-independent solution while resolutions of 963 and less are
clearly not. (b) The linear resolution of α−model simulations, N , is plotted versus kα. Simulations with
inadequate resolution are plotted as X’s, those with approximately grid-independent solutions as +’s, and
experiments that are neither clearly resolved nor clearly unresolved as boxes. The dashed lines represent
N = Ck
α indicating that a constant in the range 43.2 < C < 50.2 agrees with our data. This partially
confirms the prediction of Eq. (22) and provides a reliable method to determine the needed resolution for a
grid-independent LANS−α solution at a fixed Re.
dashed lines represent Eq. (22) with the minimal and maximal choice of C (where C0 = CRe
1/2),
that agrees with our results (i.e., 43.2 < C < 50.2). In Fig. 7 we conduct similar studies for
Re ≈ 670. We find 49.5 < C < 51.4 and again validate the predictive power of Eq. (22) for the
necessary numerical resolution for grid-independent solutions.
FIG. 7: As Fig. 6(b) but for Re ≈ 670 simulations. The dashed lines represent N = Ck1/3α indicating that
a constant in the range 49.5 < C < 51.4 agrees with our data. Note also that any power law, N ∝ kβα, with
0.30 < β < 0.46 also agrees with the data.
FIG. 8: Acceptable choices of C = C0Re
1/2, versus Reynolds number, Re, for grid-independent LANS−α.
Error bars are not confidence levels, but depict the range of values consistent with our database (N = Ck
at the four Reynolds numbers we tested. The dashed line depicts the least-squares fit with slope 0.54±0.14.
This completes the validation of Eq. (22) which predicts 0.5.
The greatest utility of the prediction, however, is due to the single constant C0 which is inde-
pendent of Reynolds number. A determination of this constant can cheaply be achieved repeating
this process for several runs for low and moderate Re, and determines the resolution requirement
for the highest Re attainable. The ranges of acceptable constants, C = C0Re
1/2, for the four
Reynolds number flows studied are plotted versus Re in Fig. 8. A power law C = C0Re
γ fits our
data with γ = 0.54± 0.14 demonstrating the final validation of the prediction, γ = 0.5, Eq. (22).
The value of the constant is found to be C0 = 2.0 ± 0.2. We made one study for the maximally-
helical ABC forcing at Re ≈ 1600 and α = 2π/25. It is consistent with a value of C0 = 1.8±0.1.
We therefore conclude that the constant C0 is not a strong function of the forcing employed or of
the scale at which the system is forced. As a result, and unlike in 2D LANS−α [36], we verify
that the size of the attractor in 3D LANS−α is smaller than that in Navier-Stokes, which is a
promising result if the LANS−α equation is going to be used as an LES. However, before doing
this, an assessment of the truncation errors introduced in discretized systems (as used to solve the
equations numerically) and a study of the optimal choice for α to capture the properties of a DNS
is needed. We consider these problems in the following section.
V. CAN LANS−α BE CONSIDERED AS A LARGE EDDY SIMULATION?
In this section, we consider the LANS−α equations as a means to an end, and consider the
solutions to their discretized equations as approximations to the Navier-Stokes solutions. We seek
numerical approximations of LANS−α that minimize the difference to a fully resolved or direct
numerical solution (DNS) of Navier-Stokes (i.e., we analyze the behavior of LANS−α solutions
in the LES framework, and call here the model a “LANS−α LES”, or in short “α-LES”). In the
LES framework, LANS−α’s turbulent stress tensor, τ̄αij , is given by (see, e.g., [42])
τ̄αij = H−1α2(∂kui∂kuj + ∂kui∂juk − ∂iuk∂juk). (24)
Previous studies have not made the distinction between grid-independent LANS−α and LANS−α
LES, though one did study convergence to grid-independent solutions at moderate Re [19]. We
find, however, a definite difference between the two approaches. We show in this section that, in
fact, LANS−α combined with truncation error yields a better fit to DNS than grid-independent
LANS−α. The resolution that yields an optimal α-LES (a terminology to be defined below) is
also found to follow Eq. (22). In the Section V A, we then address the quality and usability of the
predictions of the LANS−α model viewed as an LES.
A remark about nomenclature may be in order at this point. Traditionally, and for good reasons,
LES attempt at capturing the large-scale properties of a flow with a huge Reynolds number, as
found, e.g., in the atmosphere. In that case, the wavenumber at which the DNS is truncated
is, at best, in the inertial range and it might even be in the energy-containing range, as for the
atmospheric boundary layer with a Taylor Reynolds number Rλ ∼ 104. Of a different nature are
the modeling methods sometimes called quasi-DNS. Here, the idea is to model a flow at a given,
moderate Reynolds number but with an expense in computing resources lesser than if performing a
DNS. Under-resolved DNS fall in that category; in that case, the large-scales are presumably well
reproduced but the small scales are noisy. It is in that spirit that we now examine the properties of
the LANS−α model. We thus qualify a model as optimal in the sense of being optimal for the class
of LANS−α models examined herein; in order to avoid repetition, we also use the terminology of
alpha-optimal.
FIG. 9: (Color online.) Plot of Re ≈ 670 simulations. Average compensated energy spectra: DNS (solid
black line) and LANS−α simulations, kα = 41: N = 162 (red dotted), N = 192 (green dashed), and
N = 216 (blue dash-dotted). LANS−α at a linear resolution of 192 is approximately converged to the grid-
independent solution while a resolution of 162 is not. N = 162 does correspond, however, more closely to
the DNS spectrum. We observe, in general, that a combination of LANS−α and truncation error yields the
optimal α-LES.
In Fig. 9 with kα = 41, we plot the Re ≈ 670 DNS spectrum (solid black line) and LANS−α
spectra at three different resolutions. We observe that, while the N = 162 solution (dotted line,
red online) is not converged, it is a better approximation to the DNS than the grid-independent
LANS−α solution. For all simulations we studied, the grid-independent LANS−α solution is not
the best approximation to the DNS. Another example is given in Fig. 10 where we plot the mean
square spectral error normalized to make fair comparisons between large and small kα results,
Esq =
(Eα(k)−E(k))2
E2(k)
, (25)
where kF is the wavenumber for the forcing scale, E(k) is the DNS spectrum (in the L
2(v) norm),
Eα(k) is the LANS−α spectrum (in the H1α(u) norm), and n is the number of terms in the sum.
These errors are calculated for spectra averaged over turbulent steady-state solutions: t ∈ [16, 19]
for Re ≈ 670. We see that for a given filter or a given simulation resolution, there is a local
minimum in the error. This minimum is a balance between truncation errors and the approximation
error due to using LANS−α instead of the full Navier-Stokes equations. Due to these errors being,
in some sense, in opposition, the optimal α-LES solution is found at a lower resolution than the
grid-independent solution. Indeed, we see by examining Fig. 10 (a) that for a given filter the
combination of truncation error and the LANS−α solution is a better approximation to the DNS.
For fixed resolution, Fig. 10 (b), the optimal value for α is not zero but has some finite value.
This local minimum error shown in the figure keeps α from going to zero (kα → ∞) in dynamical
models [15]. We note, also, that the error is low for a finite range of N and kα near the minimum,
indicating that an α-LES solution may perform well for a range of parameters near the optimal
ones. We find the resolution for an optimal α-LES is also predicted by Eq. (22) (with C ≈ 47 for
Re ≈ 670, or C0 ≈ 1.8). That is, optimal α-LES resolution is just below that for grid-independent
LANS−α solutions. Having demonstrated the predictability of the resolution for grid-independent
LANS−α and of LANS−α LES given a Reynolds number and a filter, in the following section we
seek to determine sufficient conditions on the free parameter α for LANS−α to be a successful
A. Free parameter α and quality of the α-LES
In this section, we make an analysis of the LES potential of LANS−α by considering only
the grid-independent LANS−α solutions identified using Eq. (22). Note that from the results
discussed in the previous section, we expect LANS−α optimal grid-dependent α−LES approx-
imations to have better performance. In the limit of α going to zero, LANS−α Eq. (2) re-
covers the Navier-Stokes equations, Eqs. (1), but the question we address now is how small
must α be for LANS−α solutions to be good approximations to Navier-Stokes solutions. There
are several length scales that α could be related to: the forcing scale lF , the integral scale
L = 2π
E(k)k−1dk/
E(k)dk, the Taylor microscale λ, or the Kolmogorov dissipation
scale ηK . Plots of the mean square spectral errors to DNS (see Eq. (25)) versus these scales
are shown in Fig. 11. While the general trend of errors decreasing with α is apparent in all
cases, in Fig. 11(a) we see a large difference between errors at varying Reynolds numbers and
FIG. 10: Plots for Re ≈ 670 simulations. (a) Error (see Eq. (25)) versus simulation resolution for kα = 20.
The optimal (grid-dependent) LES is for a resolution of N ≈ 128 and has a much smaller error compared
to the DNS than the grid-independent LANS−α solution at higher resolution. (b) Error versus kα for
N = 128. At a given resolution the optimal value for α is not zero but occurs at a local minimal error. Any
kα ∈ [15, 25] has an error near the minimum indicating that an LES solution may perform well for a range
of parameters near the optimal ones. A constant of C = C0Re
1/2 ≈ 47 in Eq. (22) is found to correspond
with optimal α-LES approximations.
similar ratios of α to the forcing scale, lF . For a linear least-squares fit, the goodness-of-fit,
(Eactualsq − Efitsq )2, was found to be χ2 = 6.2× 10−2. The errors for Re ≈ 3300 are much
larger than for the same ratio lF/α as results at both Re ≈ 500 and Re ≈ 670. This is also the case
for the integral scale. However, the quality of the α-LES appears to be closely tied to the ratio of
α to the Kolmogorov dissipation scale. In Fig. 11(b) the errors are plotted versus the ratio of the
dissipation scale, ηK , to α. We see a very strong dependence (χ
2 = 2.5 × 10−2) between errors
for several runs with four different Reynolds numbers indicating that the quality of the LANS−α
LES approximation is a function of the ratio of α to the dissipative scale. Finally, in Fig. 11(c)
FIG. 11: Plot of errors, Eq. (25), of grid-independent solutions compared to DNS. Asterisks are for Re ≈
8300, squares for Re ≈ 3300, triangles for Re ≈ 670, and diamonds for Re ≈ 500. The single right-
most triangle in all plots corresponds to a value of α in the dissipative range (kα = 60). The norm we
employ to measure the error, Eq. (25), is no longer a good norm when dissipative scales are considered. (a)
Errors versus lF /α. No clear correlation between LES quality and the ratio of the forcing scale to α holds
independently of Reynolds numbers. (b) Errors versus ratio of dissipative scale, ηK , to α. The quality of
the LES appears to be closely tied to this ratio. (c) Errors versus ratio of Taylor wavenumber, λ, to α. The
Re ≈ 8300 experiment (asterisk) indicates that the quality of the α-LES is not tied to the Taylor scale.
the errors are plotted versus the ratio of the Taylor Scale, λ, to α. We find χ2 = 3.1 × 10−2 for a
linear least-squares fit. We note that a single experiment conducted at Re ≈ 8300 (the asterisks)
confirms that the maximal value of α is tied to the dissipation scale and not the Taylor scale. This
is more clearly demonstrated in Fig. 12 where we plot compensated energy spectra for a nearly
constant ratio λ/α at three Reynolds numbers. We see that the maximum deviation from the DNS
spectrum increases with Re. As λ/α is the same in all cases, the optimal α is not dependent on
the Taylor scale.
These findings were not accessible at lower Reynolds numbers due to inadequate separation
of scales. For example, we give in Fig. 13(a) spectral flux for DNS at Re ≈ 500, 670, and
3300 respectively. We define the kinetic energy transfer function, T (k), in Fourier space as
T (k) = −
v̂k · ̂(ω × v)dV , where (̂·) represents the Fourier transform. For LANS−α we have
Tα(k) = −
ũk · ̂(ω × u)dV where ω = ∇ × v. The flux is defined as usual from the transfer
function as
Π(α)(k) =
T(α)(k
. (26)
Only Re ≈ 3300 (and Re ≈ 8300 not pictured here) demonstrates a range of nearly constant flux
(a well-defined inertial range) before the dissipation scales. Following the scaling arguments in
Ref. [11], one effect of the α−model is to increase the time scale for the cascade of energy to small
scales. This reduces the flux as α increases (kα decreases) as do the hypothesized “rigid bodies;”
this can be seen in Fig. 13(b). (Note that in DNS at high resolution, 80% of the flux is from local
interactions which is strongly suppressed at scales smaller than α [21].) As dissipation dominates
the flux for low and moderate Reynolds number, the reduced flux of the α−model has little conse-
quence for these simulations. With a substantial inertial range, however, this reduced flux results
in a pile-up of energy for scales larger than the dissipative scale and the spectrum approaches the
k1 spectrum discussed in Section III. As a consequence of the integral conservation of energy
(Eα =
u · v) there is a corresponding decrease of energy at large scales. Consequently, as the
inertial range increases, α must be moved to smaller and smaller scales in order for LANS−α not
to alter scales larger than α. In summary, the α−model’s reduced flux of energy to small scales is
more crucial when the dissipation scale is farther away from α.
FIG. 12: Compensated averaged grid-independent energy spectra for DNS (solid) and LANS−α (dotted)
holding the ratio of Taylor scale λ to α nearly constant. Vertical dotted lines indicate kα. (a) Re ≈ 670
and kα = 35 (λ/α = 18). (b) Re ≈ 3300 and kα = 70 (λ/α = 17). (c) Re ≈ 8300 and kα = 110
(λ/α = 17). We see that the maximum deviation from the DNS increases with Re. This is due to the
greater distance between α and the dissipative scale ηK . (Note that scales larger than k = 3 are affected by
numerical truncation issues.)
FIG. 13: (Color online.) (a) Energy flux, Eq. (26), for three DNS with Re ≈ 3300 (black, solid), Re ≈ 670
(red, dotted), and Re ≈ 500 (green, dashed). No inertial range is discernible on the flux functions except
for the highest Reynolds number case. The initial plateau followed by a bump and another plateau (for the
case at the highest Reynolds number) is a result of the forcing employed. (b) Energy flux at Re ≈ 3300 for
both DNS and α−model runs; DNS is the black, solid line. See inset for LANS−α parameters. LANS−α
gives a reduced flux which is linked to the significant pile-up of energy at high wavenumber as visible in the
energy spectrum (see Fig 14). Plots of εα versus t (not shown) also show that flux decreases (on average, at
long times) with increasing α.
B. Numerical savings from employing LANS−α
If α must be directly proportional to the Kolmogorov dissipation scale, we can estimate the
LES computational savings of the LANS−α model. For the Navier-Stokes equations we have
dofNS ∝ Re9/4 and, as we verified in Section IV, for LANS−α we have dofα = C30kαRe3/2/27.
If kα is directly proportional to the Navier-Stokes dissipation wavenumber, kη, we arrive at
kη ∝ Re3/4, (27)
and, consequently,
dofLESα ∝ Re9/4. (28)
Note that for free α, dofα (dof of LANS−α) is much smaller than dofNS . But, to obtain an optimal
LES, α is tied to kη; then the resolution requirements (dof
α ) are different and the decrease in
necessary computational resolution from employing LANS−α is fixed. In fact, for the forcing and
boundary conditions employed, we find
dofLESα ≈
dofNS. (29)
We note that Eq. (28) is consistent with theoretical predictions given in Ref. [20]. Other LES
such as the similarity model [43] and the nonlinear (or gradient) model [44, 45] have also exhib-
ited the characteristic that they achieve only moderate reductions in resolution and are, therefore,
frequently used in mixed models with a Smagorinsky term (see, e.g., [3]). That such additional
terms will be required for LANS−α to reproduce the energy spectrum of high Re flows, may not
be a significant factor in its usability. Note that the usual addition of extra dissipative subgrid-
stress terms (as in the Smagorinsky model) also introduces a stronger dependence of the system
of equations with the spatial resolution, since the filter width in such models is often associated
to the maximum wavenumber in the box, kmax. In that case, it can make more sense to use grid-
dependent solutions of LANS−α (discussed at the beginning of Section V) which give an optimal
LANS−α LES, and can as a result give an extra gain in the computational costs.
We also conclude that, with the scale α being tied to the dissipation scale ηK , the model
LANS−α behaves more like a quasi-DNS by opposition to a traditional LES. Note however that
a factor of ≈ 2.3 in resolution gain translates into a factor 27 in CPU and a factor 12 in memory
savings, still a substantial gain.
FIG. 14: (Color online.) Compensated energy spectra averaged over t ∈ [8, 9], Re ≈ 3300. DNS is the solid
black line and grid-independent LANS−α solutions are shown as (red online) dotted (kα = 70), (green)
dashed (kα = 40), and (blue) dash-dotted (kα = 13) lines, respectively. A single LANS−α LES is shown
as a (pink) dash-triple-dotted line (kα = 40, N = 384). The LES is seen to better approximate the DNS
spectrum than the grid-independent solution for the same value of α (2π/40). As α is increased the energy
spectrum approaches the k1 spectrum discussed in Section III B.
VI. LANS−α AT VERY HIGH REYNOLDS NUMBER
In this section, we compare and contrast LANS−α and Navier-Stokes solutions at high
Reynolds number. Using results of previous sections for optimal resolution and the necessary
value of α to approximate DNS, we now evaluate both grid-independent LANS−α solutions and
a single LANS−α LES for a highly turbulent flow (Re ≈ 3300, Rλ ≈ 790). We calculate grid-
independent solutions for kα = 70 (N = 512), for kα = 40 (N = 512), and for kα = 13
(N = 384). A LANS−α LES solution is computed for kα = 40 (N = 384). Averaged com-
pensated energy spectra are shown in Fig. 14. We see that the optimal LANS−α LES is a better
approximation of the DNS spectra than the grid-independent LANS−α for the same value of α
(2π/40). We also see that if α is increased further, the energy spectrum approaches the k1 spectrum
discussed in Section III B.
Fig 15 is a perspective volume rendering of the enstrophy density ω2 (ω · ω̄ for LANS−α) for
the DNS, kα = 70 LANS−α, and kα = 13 LANS−α. Due to the late time depicted here (t = 9,
longer than a Lyapunov time) there can be no point-by-point comparison between the simulations.
However, we note that the helical structure of the vortex tubes is preserved by the α−model but
that the tubes themselves are shorter and somewhat thicker for large values of α. As was noted
for moderate Reynolds numbers, this is due to LANS−α suppressing vortex stretching dynamics
FIG. 15: (Color online.) Rendering of enstrophy density ω2 (ω · ω̄ for LANS−α). Due to the late time
depicted here (t = 9, longer than a Lyapunov time) there can be no point-by-point comparison between
the simulations. Instead, regions with approximately the same dimensions are selected around vortex tubes.
Velocity v field lines are also shown illustrating the helical nature of the tubes which is seen to be captured
by LANS−α. (a) DNS. The thick bars represent, from top to bottom, the Taylor scale λ and the dissipative
scale ηK , respectively. For LANS−α results the scale α is depicted between these two. (b) kα = 70,
N = 512. (c) kα = 13, N = 384. We see that, for large values of α, the vortex tubes become shorter and
somewhat thicker.
without changing its qualitative features [9]. This is in contrast to 2D LANS−α where the vorticity
structures are seen to get thinner as α increases [36]. This could also be related to the scaling
differences between 2D and 3D LANS−α. It has been claimed that the development of helical
structures in turbulent flows can lead to the depletion of nonlinearity and the quenching of local
interactions [46, 47]. The depletion of energy transfer due to local interactions at some cutoff in
wavenumber is also believed to bring about the bottleneck effect [22, 48, 49, 50]. Consistent with
these results, in 2D LANS−α (where the vorticity structures are more fine than Navier-Stokes) the
spectrum is steeper and in 3D LANS−α (where the vorticity structures are shorter but fatter than
Navier-Stokes) the spectrum is shallower.
FIG. 16: (Color online.) Compensated 3rd-order structure function versus length l (a horizontal line scales
with l). Structure functions corresponding to the Kármán-Howarth theorem are depicted (S3 for DNS,
3 ≡ 〈(δu)2δv〉 for LANS−α). Labels are as in Fig. 14. The dotted vertical lines indicate the various α’s.
A small inertial range for the DNS near l = 1 is reproduced by LANS−α. The largest α (2π/13) exhibits
a second inertial range at scales just smaller than α (〈(δu)2δv〉 ∼ l is consistent with Eq. (15)).
Figure 16 shows the third-order (mixed) structure functions corresponding to the Kármán-
Howarth theorems versus length l. For the DNS, we show S3 ≡ 〈δv3〉 and Sα3 ≡ 〈(δu)2δv〉
for LANS−α. The dotted vertical lines indicate the various α’s. A small inertial range for the
DNS near l = 1 is reproduced by all LANS−α results. The largest α (2π/13) exhibits a second
inertial range at scales just smaller than α (〈(δu)2δv〉 ∼ l is consistent with Eq. (15)). We note
this is the first demonstration of third-order structure functions in LANS−α consistent with a K41
inertial range followed by an α inertial range and finally a dissipative range. Next, we observe the
scaling of the longitudinal structure functions,
Sp(l) ≡ 〈|δv‖|p〉, (30)
where we again replace the H1α norm for the L
2 norm in the case of LANS−α,
Sαp (l) ≡ 〈|δu‖δv‖|p/2〉. (31)
We utilize the extended self-similarity (ESS) hypothesis [51, 52, 53] which proposes the scaling
Sp(l) ∝ S3(l)ξp (32)
or, for LANS−α,
Sαp (l) ∝ 〈(δu)2δv〉ξp. (33)
We display our results in Fig. 17. We note that for LANS−α, the third-order exponent is not equal
to unity, contrary to the Navier-Stokes case. The Kármán-Howarth theorem implies 〈(δu)2δv〉 ∼ l,
not Sα3 (l) ∼ l. We measured the deviation from linearity for each experiment (not depicted here)
and found that LANS−α becomes more intermittent as α increases (kα = 13 is slightly more
intermittent than the DNS). As artificially dropping local small-scale interactions gives enhanced
intermittency [54, 55], this increased intermittency is the expected result of LANS−α reducing in-
teractions at scales smaller than α. We note, however, that even with such a large filter, LANS−α
is a good approximation to the intermittency properties of the DNS. This is surprising given its en-
ergy spectrum and reduced flux in the inertial range. It is probably linked to the fact that LANS−α
preserves global properties (in an H1 sense) of the Navier-Stokes equations and that these proper-
ties are important to the dynamics of small scales as measured by high-order structure functions.
VII. CONCLUSIONS
We computed solutions of the Lagrangian-Averaged Navier-Stokes α−model (LANS−α) in
three dimensions for significantly higher Reynolds numbers (up to Re ≈ 8300) than have previ-
ously been accomplished and performed numerous forced turbulence simulations of LANS−α to
study their equilibrium states. The results were compared to DNS for Re ≈ 500, 670, 3300, and
8300 , the last performed on a grid of 20483 points. We note that there are two ways to view the
LANS−α simulations: as converged or “grid-independent” solutions of the LANS−α equations
or as large-eddy simulations (α−LES) which include grid effects. We found a definite difference
between the two approaches in that the fully-converged grid-independent LANS−α is not the best
FIG. 17: (Color online.) Structure function scaling exponent ξp versus order p. Black X’s are shown
for the DNS. Grid-independent LANS−α are shown as (red online) boxes (kα = 70), as (green) triangles
(kα = 40), as (blue) diamonds (kα = 13). LANS−α LES (kα = 40, N = 384) is shown as (pink) asterisks.
The dashed line indicates K41 scaling and the solid line the She-Lévêque (SL) formula [56].
approximation to a DNS of Navier-Stokes. Instead, the minimum error is a balance between trun-
cation errors and the approximation error due to using LANS−α instead of the full Navier-Stokes
equations. Due to these errors being, in some sense, in opposition, the optimal α-LES solution
was found at a lower resolution than the grid-independent solution (the error was low for a finite
range of N and α near the minimum, indicating that a LANS−α viewed as an LES solution may
perform well for a range of parameters). Unlike the 2D case [36], 3D LANS−α has been shown
to be a subgrid model (i.e., it reduces the resolution requirements of a given computation). This
difference between 2D and 3D LANS−α indicates that other α−models (as the LAMHD−α Eqs.
[57, 58] or the BV−α Eqs. [42]) may behave differently and studies of these systems at high
resolution may be required.
We confirm the presence of the theoretically predicted l3 scaling of the third-order structure
function (corresponding to a k−1 scaling of the energy spectrum) [11, 16, 37] through its bound on
the number of degrees of freedom for LANS−α [11], in the structure functions of the smoothed
velocity in simulations with large α, and in the spectrum of specific spatial portions of the flow.
In so doing, we have validated the predictive power of the bound dofα < Cα−1Re3/2, for the nu-
merical resolution for grid-independent LANS−α solutions and for optimal LANS−α LES (with
a separate constant of proportionality). The great utility of the prediction is that the single con-
stant can cheaply be determined at low and moderate Reynolds number and predicts the resolution
requirement for the highest Reynolds numbers attainable. We further found no great change in
this single constant when employing the non-helical Taylor-Green or the maximally-helical ABC
forcings.
However, the small scale (kα ≫ 1) LANS−α spectrum was observed to be k+1. We attribute
this to the frozen-in-turbulence closure employed in deriving the α−model. For scales smaller
than α, portions of the smoothed flow u are locked into “rigid bodies.” By “rigid bodies,” we
mean the internal degrees of freedom are frozen and these portions give no contribution to the
energy cascade. This is consistent both with the observed k+1 spectrum and with field increments
δu‖ being observed to be approximately zero over a large portion (compared to Navier-Stokes)
of the flow. The turbulent energy cascade occurs in the space between these “rigid” portions.
While the k−1 portions are subdominant to the k+1 portions in the energy spectrum, they prevail
in the cascade and hence both the structure functions and the degrees of freedom of the LANS−α
attractor.
We find that both of these scalings (k+1 and k−1) contribute to a reduction of flux at constant
energy (i.e., the dissipation is reduced as has previously been observed in 2D calculations [59]).
This leads to a shallower (or even growing) energy spectrum as α increases. Thus, for LANS−α
viewed as an LES to reproduce the Navier-Stokes energy spectrum it is necessary that α be not
much larger than the dissipation scale (α / 4ηK independent of Reynolds number); in that sense,
it can be considered as a quasi-DNS as opposed to a traditional LES, substantially larger Reynolds
numbers being modeled in the latter case, leading to substantially larger gain in resolution. As
a consequence, the computational savings of LANS−α is fixed and not a function of Reynolds
number. (However, and unlike the 2D case, the 3D α−model does give a computational saving
when used as a LES.) This result was not accessible at lower Reynolds numbers due to inadequate
separation of scales. However, in one previous study for decaying turbulence with energy initially
mostly at low wavenumbers (k = 3), it was evident that as time evolved and energy moved to
smaller scales, the resolution requirements of LANS−α increased [17]. Other LES such as the
similarity model [43] and the nonlinear (or gradient) model [44, 45] have also exhibited the char-
acteristic that resolution may be decreased only modestly and are, therefore, frequently used in
mixed models with a Smagorinsky term (see e.g., [3]). That such additional terms will be required
for LANS−α to reproduce the energy spectrum of high Re flows, may not be a significant factor
in its usability.
We compared and contrasted LANS−α to a DNS at Re ≈ 3300 considering both structures
and high-order statistics such as the longitudinal structure functions which are related with inter-
mittency. With an appropriate choice of α we were able to observe a Navier-Stokes inertial range
followed by LANS−α inertial range at scales smaller than α. For this second inertial range we
again observed a k+1 energy spectrum. As α increased, we noted a change in the aspect ratio of
vortex tubes (they became shorter and fatter). This can be related to quenching of local small-
scale interactions at scales smaller than α and, thus, to the shallower spectrum for 3D LANS−α
[22, 46, 47, 48, 49, 50]. Therefore, in 2D LANS−α (where the vorticity structures are more fine
than Navier-Stokes) the spectrum is steeper [36] and in 3D LANS−α (where the vorticity struc-
tures are shorter but fatter than Navier-Stokes) the spectrum is shallower. Finally, an examination
of the longitudinal structure functions indicate that intermittency is increased as the parameter α
is increased consistent with the suppression of local small-scale interactions at scales smaller than
α [54, 55].
The elimination of the faster and faster interactions among smaller and smaller scales through
the modified nonlinearity in LANS−α (together with the discrepancy between its solutions and
Navier-Stokes solutions) highlights the importance of these interactions down to scales only
slightly larger than the dissipative scale. That is, by removing these interactions anywhere in
the inertial range (e.g., α ' 4ηK), the resulting energy spectrum was found to differ from the DNS
at scales larger than α. The intermittency properties of the DNS, however, were well reproduced
even with large filters. Noting this, if LANS−α’s k1 energy spectrum is not important for a given
application, much greater reductions in resolution can be achieved. Future work should address
whether this may be remedied in a LANS−α LES by the inclusion of another (dissipative) model
for these interactions, or (in the case of magneto-hydrodynamics [57, 58] whether this problem is
less significant because of the presence of greater spectral nonlocality [60, 61, 62]. The effect of
LANS−α on the detailed scale-by-scale energy transfer should also be investigated as our results
indicate that a model for local small-scale interactions would improve the α−model. Another
direction of future research is to explore other reduced LANS−α models, Clark−α and Leray-α,
which break the frozen-in-turbulence closure and, also, the conservation of circulation. Finally,
note that because of its greater mathematical tractability, LANS−α possibly allows for a better un-
derstanding of multi-scale interactions in turbulent flows thus modeled; therefore, detailed studies
such as the one presented here may, in fine, allow for a better understanding of turbulence itself.
Acknowledgments
Computer time was provided by NCAR and by the National Science Foundation Terascale
Computing System at the Pittsburgh Supercomputing Center. The NSF Grant No. CMG-0327888
at NCAR supported this work in part and is gratefully acknowledged. Three-dimensional visual-
izations of the flows were done using VAPOR. The authors would like to express their gratitude
for valuable discussions with Bob Kerr.
[1] P. J. Mason, Quarterly Journal of the Royal Meteorological Society 120, 1 (1994).
[2] M. Lesieur and O. Metais, Annual Review of Fluid Mechanics 28, 45 (1996).
[3] C. Meneveau and J. Katz, Annual Review of Fluid Mechanics 32, 1 (2000).
[4] J. R. Kulkarni, L. K. Sadani, and B. S. Murthy, Boundary-Layer Meteorology 90, 217 (1999).
[5] G. Heinemann, Theoretical and Applied Climatology 83, 35 (2006).
[6] D. D. Holm, J. E. Marsden, and T. S. Ratiu, Adv. in Math. 137, 1 (1998).
[7] S. Chen, C. Foias, D. D. Holm, E. Olson, E. S. Titi, and S. Wynne, Physical Review Letters 81, 5338
(1998).
[8] S. Chen, C. Foias, D. D. Holm, E. Olson, E. S. Titi, and S. Wynne, Physica D Nonlinear Phenomena
133, 49 (1999).
[9] S. Chen, D. D. Holm, L. G. Margolin, and R. Zhang, Physica D Nonlinear Phenomena 133, 66 (1999).
[10] S. Chen, C. Foias, D. D. Holm, E. Olson, E. S. Titi, and S. Wynne, Physics of Fluids 11, 2343 (1999).
[11] C. Foias, D. D. Holm, and E. S. Titi, Physica D Nonlinear Phenomena 152-153, 505 (2001).
[12] D. D. Holm, J. E. Marsden, and T. S. Ratiu, Physical Review Letters 80, 4173 (1998).
[13] D. D. Holm, Chaos 12, 518 (2002).
[14] D. D. Holm, Physica D Nonlinear Phenomena 170, 253 (2002).
[15] H. Zhao and K. Mohseni (2004), arXiv:physics/0408113.
[16] A. Cheskidov, D. D. Holm, E. Olson, and E. S. Titi, Proceedings of the Royal Society of London
A461, 629 (2005).
[17] K. Mohseni, B. Kosović, S. Shkoller, and J. E. Marsden, Physics of Fluids 15, 524 (2003).
[18] B. J. Geurts and D. D. Holm, in Turbulent Flow Computation, edited by D. Drikakis and B. J. Geurts
(Kluwer Academic Publishers, London, 2002), pp. 237+.
[19] B. J. Geurts and D. D. Holm, Journal of Turbulence 7, 1 (2006).
[20] J. D. Gibbon and D. D. Holm, Physica D Nonlinear Phenomena 220, 69 (2006), nlin/0603059.
[21] A. Alexakis, P. D. Mininni, and A. Pouquet, Physical Review Letters 95, 264503 (2005),
physics/0507144.
[22] P. D. Mininni, A. Alexakis, and A. Pouquet, Phys. Rev. E 74, 016303 (2006), physics/0602148.
[23] D. O. Gómez, P. D. Mininni, and P. Dmitruk, Physica Scripta p. 123 (2005).
[24] D. O. Gómez, P. D. Mininni, and P. Dmitruk, Advances in Space Research 35, 899 (2005).
[25] G. I. Taylor and A. E. Green, Proceedings of the Royal Society of London A158 (1937).
[26] P. D. Mininni, Y. Ponty, D. C. Montgomery, J.-F. Pinton, H. Politano, and A. Pouquet, Astrophys. J.
626, 853 (2005).
[27] M. Brachet, Academie des Sciences Paris Comptes Rendus Serie Sciences Mathematiques 311, 775
(1990).
[28] Y. Ponty, P. D. Mininni, D. C. Montgomery, J.-F. Pinton, H. Politano, and A. Pouquet, Physical Review
Letters 94, 164502 (2005).
[29] T. de Kármán and L. Howarth, Proceedings of the Royal Society of London A164, 192 (1938).
[30] D. D. Holm, Journal of Fluid Mechanics 467, 205 (2002).
[31] U. Frisch, Turbulence, The Legacy of A. N. Kolmogorov (Cambridge University Press, Cambridge,
UK, 1995).
[32] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR 30, 299 (1941), reprinted in Proc. R. Soc. Lond. A (1991)
434, 9-13.
[33] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR 31, 538 (1941).
[34] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR 32, 16 (1941), reprinted in Proc. R. Soc. Lond. A (1991)
434, 15-17.
[35] R. H. Kraichnan, Physics of Fluids 10, 1417 (1967).
[36] E. Lunasin, S. Kurien, M. Taylor, and E. Titi, ArXiv Physics e-prints (2007), physics/0702196.
[37] C. Cao, D. D. Holm, and E. S. Titi, Journal of Turbulence 6, 19 (2005).
[38] A. A. Ilyin, E. M. Lunasin, and E. S. Titi, Nonlinearity 19, 879 (2006).
[39] C. Cichowlas, P. Bonaı̈ti, F. Debbasch, and M. Brachet, Physical Review Letters 95, 264502 (2005),
nlin/0410064.
[40] B. J. Geurts and J. Fröhlich, Physics of Fluids 14, L41 (2002).
[41] J. Meyers, B. J. Geurts, and M. Baelmans, Physics of Fluids 15, 2740 (2003).
[42] D. D. Holm and B. Nadiga, J. Phys. Oceanogr. 33, 2355 (2003).
[43] J. Bardina, J. H. Ferziger, and W. C. Reynolds, eds., Improved subgrid-scale models for large-eddy
simulation (1980).
[44] A. Leonard, in Turbulent diffusion in environmental pollution; Proceedings of the Second Symposium,
Charlottesville, Va., April 8-14, 1973. Volume A. (A75-30951 13-47) New York, Academic Press, Inc.,
1974, p. 237-248. NASA-supported research. (1974), pp. 237–248.
[45] R. A. Clark, J. H. Ferziger, and W. C. Reynolds, Journal of Fluid Mechanics 91, 1 (1979).
[46] H. K. Moffatt and A. Tsinober, Annual Review of Fluid Mechanics 24, 281 (1992).
[47] A. Tsinober, An Informal Introduction to Turbulence (Kluwer Academic Publishers, Dordrecht, 2001).
[48] J. R. Herring, D. Schertzer, M. Lesieur, G. R. Newman, J. P. Chollet, and M. Larcheveque, Journal of
Fluid Mechanics 124, 411 (1982).
[49] D. Lohse and A. Müller-Groeling, Physical Review Letters 74, 1747 (1995), 1994chao.dyn..5002L.
[50] D. O. Martı́nez, S. Chen, G. D. Doolen, R. H. Kraichnan, L.-P. Wang, and Y. Zhou, Journal of Plasma
Physics 57, 195 (1997).
[51] R. Benzi, S. Ciliberto, C. Baudet, G. Ruiz Chavarria, and R. Tripiccione, Europhysics Letters 24, 275
(1993).
[52] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Massaioli, and S. Succi, Phys. Rev. E 48, R29
(1993).
[53] R. Benzi, L. Biferale, S. Ciliberto, M. V. Struglia, and R. Tripiccione, Phys. Rev. E 53, R3025 (1996),
1995chao.dyn..9013B.
[54] J.-P. Laval, B. Dubrulle, and S. Nazarenko, Physics of Fluids 13, 1995 (2001), physics/0101036.
[55] B. Dubrulle, J.-P. Laval, S. Nazarenko, and O. Zaboronski, Journal of Fluid Mechanics 520, 1 (2004),
physics/0304035.
[56] Z. She and E. Lévêque, Physical Review Letters 72, 336 (1994).
[57] J. Pietarila Graham, P. D. Mininni, and A. Pouquet, Phys. Rev. E 72, 045301(R) (2005).
[58] J. Pietarila Graham, D. D. Holm, P. Mininni, and A. Pouquet, Physics of Fluids 18, 045106 (2006).
[59] D. Biskamp and E. Schwarz, Physics of Plasmas 8, 3282 (2001).
[60] A. Alexakis, P. D. Mininni, and A. Pouquet, Phys. Rev. E 72, 046301 (2005).
[61] P. Mininni, A. Alexakis, and A. Pouquet, Phys. Rev. E 72, 046302 (2005).
[62] A. Alexakis, P. D. Mininni, and A. Pouquet, Astrophys. J. 640, 335 (2006), physics/0509069.
	Introduction
	Technique
	Inertial range scaling of LANS-
	l3 scaling of third-order structure function derived from the Kármán-Howarth theorem for LANS-
	Subdominance of the k-1 energy spectrum and rigid-body motions
	Resolution requirements for grid-independent LANS- solutions: Size of attractor
	Can LANS- be considered as a Large Eddy Simulation?
	Free parameter  and quality of the -LES
	Numerical savings from employing LANS-
	LANS- at very high Reynolds number
	Conclusions
	Acknowledgments
	References
