Optimal subgrid scheme for shell models of turbulence

TL;DR

This paper develops a theoretical framework for optimal sub-grid closures in shell models of turbulence, demonstrating that low-order approximations can accurately reproduce large-scale dynamics and properties, with some scale-dependent discrepancies.

Contribution

It introduces a systematic approach to define and approximate optimal sub-grid closures based on scale correlations, extending ideas from Kolmogorov's hypotheses.

Findings

Low-order closures reproduce known turbulence properties.

Discrepancies occur near the sub-grid threshold scale.

Connections to Large Eddy Simulations are discussed.

Abstract

We discuss a theoretical framework to define an optimal sub-grid closure for shell models of turbulence. The closure is based on the ansatz that consecutive shell multipliers are short-range correlated, following the third hypothesis of Kolmogorov formulated for similar quantities for the original three-dimensional Navier-Stokes turbulence. We also propose a series of systematic approximations to the optimal model by assuming different degrees of correlations across scales among amplitudes and phases of consecutive multipliers. We show numerically that such low-order closures work well, reproducing all known properties of the large-scale dynamics including anomalous scaling. We found small but systematic discrepancies only for a range of scales close to the sub-grid threshold, which do not tend to disappear by increasing the order of the approximation. We speculate that the lack of convergence might be due to a structural instability, at least for the evolution of very fast degrees of freedom at small scales. Connections with similar problems for Large Eddy Simulations of the three-dimensional Navier-Stokes equations are also discussed. Postprint version of the article published on Phys. Rev. E 95, 043108 (2017) DOI: 10.1103/PhysRevE.00.003100

Equations

$${\dot{u}}_n=k_n{Q}_n-\nu k_n^2{u}_n,n=1,\dots ,N,$$

$$Q_n=Q(u_{n-1},u_n,u_{n+1})=\sum _{i,j\in \{-1,0,1\}}{a}_{ij}{u}_{n+i}{u}_{n+j}.$$

$$u_0=u_0(t).$$

$$n_K\approx -\frac{3}{4}{\log }_{\lambda }\nu$$

$$\frac{\partial P}{\partial t}+\sum _{n=1}^N\frac{\partial }{\partial u_n}\left({\dot{u}}_{n}P\right)=0.$$

$${\dot{u}}_n=k_n{Q}_n,n=1,\dots ,s.$$

$$P_s(u_1,\dots ,u_s;t)=\int P(u_1,\dots ,u_N;t)\prod _{m=s+1}^{N}d{u}_m.$$

$$\frac{\partial P_s}{\partial t}+\sum _{n=1}^N\int \frac{\partial }{\partial u_n}\left({\dot{u}}_{n}P\right)\prod _{m=s+1}^{N}d{u}_m=0.$$

$$\frac{\partial P_s}{\partial t}+\sum _{n=1}^s\frac{\partial }{\partial u_n}\left(k_n{R}_n{P}_s\right)=0,$$

$$R_n(u_1,\dots ,u_s;t)=\int Q_n\frac{P}{P_s}\prod _{m=s+1}^{N}d{u}_m.$$

$${\dot{u}}_n=k_n{R}_n,n=1,\dots ,s.$$

$$R_{n}=\left\{\begin{array}[]{ll}Q_{n},&n=1,\ldots,s-1;\\[7.0pt] \displaystyle\int Q_{s}\frac{P_{s+1}}{P_{s}}\,du_{s+1},&n=s.\end{array}\right.$$

$$\frac{P_{s+1}}{P_s}d{u}_{s+1}=P_{\text{cond}}(u_{s+1}|u_s,\dots ,u_1;t)d{u}_{s+1},$$

$$P_{\text{cond}}(u_{s+1}|u_s,\dots ,u_1;t)d{u}_{s+1}={\tilde{P}}_{\text{cond}}(z_{s+1}|z_s,\dots ,z_1;t)d{z}_{s+1}.$$

$${\tilde{P}}_{\text{cond}}(z_{s+1}|z_s,\dots ,z_1;t)={\tilde{P}}_{\text{uni}}(z_{s+1}|z_s,z_{s-1},\dots )$$

$$R_s(u_s,u_{s-1},\dots )=\int Q(u_{s-1},u_s,u_{s+1}){\tilde{P}}_{\text{uni}}(z_{s+1}|z_s,z_{s-1},\dots )d{z}_{s+1}.$$

$$Q(u_{s-1},u_s,u_{s+1})=u_s^{2}Q(z_s^{-1},1,z_{s+1}).$$

$$R_s(u_s,u_{s-1},\dots )=u_s^2{\tilde{R}}_s(z_s,z_{s-1},\dots ),$$

$${\tilde{R}}_s(z_s,z_{s-1},\dots )=\int Q(z_s^{-1},1,z_{s+1}){\tilde{P}}_{\text{uni}}(z_{s+1}|z_s,z_{s-1},\dots )d{z}_{s+1}.$$

$$Q(u_{n-1},u_n,u_{n+1})=u_{n-1}^2-\lambda u_n{u}_{n+1}.$$

$$u_n\approx a{k}_n^{-1/3}$$

$$z_n=\frac{u_n}{u_{n-1}}={\lambda }^{-1/3}$$

$${\tilde{P}}_{\mathrm{uni}}(z_{s+1}|z_s,\dots )=\delta (z_{s+1}-{\lambda }^{-1/3}).$$

$${\tilde{R}}_s=\int \left(\frac{1}{z_s^2}-\lambda z_{s+1}\right)\delta (z_{s+1}-{\lambda }^{-1/3})d{z}_{s+1}=\frac{1}{z_s^2}-{\lambda }^{2/3}.$$

$$R_s(u_s,u_{s-1})=u_{s-1}^2-{\lambda }^{2/3}{u}_s^2.$$

$$\dot{u}_{n}=\left\{\begin{array}[]{ll}k_{n}\left(u_{n-1}^{2}-\lambda u_{n}u_{n+1}\right),&n=1,\ldots,s-1;\\[5.0pt] k_{s}\left(u_{s-1}^{2}-\lambda^{2/3}u_{s}^{2}\right),&n=s.\end{array}\right.$$

$$Q_n=Q(u_{n-2},\dots ,u_{n+2})=i\left(\lambda u_{n+2}{u}_{n+1}^{\ast }+b{u}_{n+1}{u}_{n-1}^{\ast }-c{\lambda }^{-1}{u}_{n-1}{u}_{n-2}\right),$$

$$\frac{\partial P}{\partial t}+\sum _{n=1}^N\left[\frac{\partial }{\partial a_n}\left({\dot{a}}_{n}P\right)+\frac{\partial }{\partial b_n}\left({\dot{b}}_{n}P\right)\right]=0,$$

$$\frac{\partial P_s}{\partial t}+\sum _{n=1}^N\left[\frac{\partial }{\partial a_n}\left(k_n{A}_n{P}_s\right)+\frac{\partial }{\partial b_n}\left(k_n{B}_n{P}_s\right)\right]=0,$$

$${\dot{u}}_n=k_n{R}_n,n=1,\dots ,s.$$

Full text

Postprint version of the article published on Phys. Rev. E 95, 043108 (2017) DOI: 10.1103/PhysRevE.00.003100

I Introduction

Three-dimensional turbulence is a multiscale phenomenon triggered when the nonlinear transport terms in the Navier–Stokes (NS) equations are much more intense than the viscous linear damping Frisch (1995). The control parameter is given by the Reynolds number, $Re=u_{0}l_{0}/\nu$, made out of the typical root mean square velocity, $u_{0}$, the typical length scale, $l_{0}$ and the kinematic viscosity, $\nu$. It is an empirical fact that in the turbulent regime the flow develops a dissipative anomaly: a $Re$-independent energy transfer, from the scale where the external forcing is acting till the viscous range. The energy transfer mechanism is characterized by anomalous scaling laws and by a highly non-Gaussian and intermittent statistics Frisch (1995). It is fair to say that we do not yet possess neither the analytical nor the numerical tools to fully quantify turbulence for three-dimensional flows.

Shell models provide a natural playground for fundamental studies of developed turbulence Frisch (1995); Bohr et al. (1998); Biferale (2003); Ditlevsen (2010). These models allow accurate numerical simulations and possess nontrivial properties of the Kolmogorov–Obukhov theory for turbulence at high Reynolds numbers: a forward energy transfer, a dissipative anomaly and intermittency with anomalous scaling similar to what observed for the original three-dimensional NS equations. The idea is to build simple models sharing the key statistical properties of the turbulent energy cascade. In this paper, we focus on the Sabra shell model L’vov et al. (1998) (a modified version of the Gledzer–Ohkitani–Yamada model Gledzer (1973); Ohkitani and Yamada (1989); Bohr et al. (1998)), which is obtained by reducing dynamics to a discrete sequence of shells $|\mathbf{k}|=k_{n}$ in the Fourier space for the geometric progression of wavenumbers $k_{n}=k_{0}\lambda^{n}$, $n=1,2,3,\ldots$ (we use $k_{0}=1$ and $\lambda=2$). The turbulent “flow” is described by complex velocity variables $u_{n}(t)$, which mimic the velocity increments at the corresponding shells, $u_{n}\sim\delta_{\ell}v=v(\ell)-v(0)$. Thus, the shell variable $u_{n}$ characterizes the velocity fluctuation at scale $\ell\sim 1/k_{n}$.

One of the main theoretical and applied challenges in the theory of turbulence consists in closing the equations of motion on a coarser grid, i.e. to derive a model for the small-scale degrees of freedom to be used to evolve the variables at large scales. The problem is key for Large Eddy Simulations (LES), a set of applied numerical tools meant to reduce the computational costs to simulate high Reynolds number turbulence Meneveau and Katz (2000); Pope (2001); Sagaut (2006); Lesieur (1987). The problem is also key from a theoretical point of view, because, if successful, would imply a complete control on the energy-transfer mechanism at all scales. The main difficulties to accomplish the goal for the three-dimensional NS case are connected to the extremely complicated functional and statistical dependency of the unresolved sub-grid variables from the resolved ones, the legacy of the strong non-linear character of the dynamical evolution together with the strong non-local coupling in both real and Fourier space of the original equations. In fact, despite many advancements, the problem of finding an optimal sub-grid model to be applied in LES is considered still open.

Our aim here is to show that this task can be accomplished for the Sabra shell model in a way that accurately describes the statistics of subgrid scales. The good news is that the simplified structure of the non-linear terms allows for a precise theoretical and numerical analysis of the statistical coupling among resolved and unresolved shells. As a result, it is possible to define what would be the optimal closure, in theory. The bad news is that the problem is not of easy implementation even in this case and that it is difficult to figure out a systematic protocol of more and more complex sub-grid models which converge toward the “optimal” one. The main idea is to close the sub-grid dynamics in terms of multi-scale correlations among multipliers, i.e. ratios among consecutive shell variables Benzi et al. (1993); Eyink et al. (2003); Benzi et al. (2004); Friedrich and Grauer (2016); Renner et al. (2001); Ragwitz and Kantz (2001); Schmitt (2003). The approach goes back to the third hypothesis of Kolmogorov Frisch (1995), made to disentangle universal small-scale fluctuations from non-universal coupling with the large-scale motion. Differently from the original case of NS equations, multipliers in shell models follow a simple non-linear dynamical evolution. It is therefore possible to manipulate them and to make predictions Benzi et al. (1993); Eyink et al. (2003). It turns out that it is crucial to distinguish the correlations among their amplitude and their phases. In this paper we first show how to define a formal optimal sub-grid model. The model is still too complicated to be implemented in practice, being defined in terms of the conditional probability of a few sub-grid variables with all resolved degrees of freedom, a task out of reach even for simplified dynamics as for the case of shell models. Then, we show how to develop a series of simple approximations for the sub-grid closure that work well, i.e. they are able to quantitatively reproduce the large-scale dynamics except for a short range of shells close to the cutoff. We also show that the observed deviations are Reynolds independent, i.e. the discrepancies remain localized to a limited number of scales close to the cut-off independently of the intensity of turbulence. Unfortunately, numerics demonstrates that the proposed systematic protocol of more and more refined closures denies a controllable convergence to the optimal model at small scales. We speculate that this might be due to non-trivial strong sensitivity of the structure of the attractor on the small-scale closure, a sort of breaking of ergodicity at fast small-scale degrees-of-freedom. A comment on the potential connections with the equivalent problem to find an optimal sub-grid closure for LES of turbulence is also proposed.

The paper is organized as follows. In Sec. II we discuss the set-up on how to define the optimal sub-grid model for a general shell model. In Sec. III we show how to implement the third hypothesis of Kolmogorov to define a systematic universal closure for the sub-grid model. In Sec. IV we show how this procedure works in simple shell models where the dynamical evolution is not intermittent. In Sec. V we formulate it for the case of the Sabra model, one of the most popular and studied shell models for turbulence. In the same section, we propose and apply a set of approximations to the optimal closure for the Sabra model and discuss their pluses and minuses. Conclusions follow in Sec. VI.

II Reduced system for a probability density

Shell models are dynamical systems which mimic the fluid dynamics by considering a geometric progression of wavenumbers, $k_{n}=k_{0}\lambda^{n}$, for some fixed $\lambda>1$ and $n=1,2,\ldots,N$. Each wavenumber defines a shell $|k|=k_{n}$ in Fourier space represented by one or several shell variables, which describe intensity of the flow at a corresponding scale. Characteristic scale in physical space can be defined as $\ell\sim 1/k_{n}$. Thus, $n\sim 1$ corresponds to large scales $\ell\sim 1/k_{0}$, while $n\sim N$ yields the smallest scales of the system.

For simplicity, we start by assuming real shell variables $u_{n}$ and considering a model with only the nearest-shell interaction. These assumptions are made in order to present the derivations in a simple and clear form, and then we extend the results to general shell models in Sec. 5. Equations of our simple shell model read

[TABLE]

with the quadratic nonlinear term coupling only the nearest neighbors:

[TABLE]

A boundary condition must be supplied for the initial shell

[TABLE]

The total number of shells $N$ is assumed to be large enough leading to the strong decay due to viscosity at small scales, i.e., $u_{N}\approx 0$. Note that we use no explicit forcing term in Eq. (1), with the excitation performed by the boundary condition (3) as it is typical for realistic flows. The nonlinear term in (2) must be chosen such that the system possesses an inviscid invariant $E=\frac{1}{2}\sum u_{n}^{2}$ called the energy.

The number of shells involved in the dynamics depends on viscosity $\nu$. Considering the integral scales of the system $L\sim 1/k_{0}\sim 1$ and $T\sim 1$, the Reynolds number is defined simply as $\mathrm{Re}=1/\nu$. In statistically stationary regime with large Reynolds numbers, one can distinguish three ranges of scales with qualitatively different behavior Frisch (1995). The range of large scales, $n\sim 1$, is called the forcing range, as it is influenced by the boundary conditions producing the energy input into the system. The energy dissipates at small scales $n\gtrsim n_{K}$ of the viscous range. The estimate

[TABLE]

can be obtained by comparing $\ell\sim 1/k_{n}$ with the Kolmogorov scale $\eta=(\nu^{3}/\varepsilon)^{1/4}$, where $\varepsilon\sim 1$ is the rate of energy dissipation Frisch (1995). For large Reynolds numbers (small viscosity) the forcing range, where energy is injected, is separated from the viscous range, where it dissipates. The intermediate range with $L\gg 1/k_{n}\gg\eta$ is called the inertial interval. In the inertial interval, both forcing and viscosity can be neglected leading to a positive mean energy flux $\varepsilon$ from larger to smaller scales, called the energy cascade.

We will consider the evolution of a statistical ensemble, corresponding to some probability distribution as initial condition. We denote by $P(u_{1},\ldots,u_{N};t)$ a probability density of the shell variables at time $t$. Time dependence of this distribution is governed by the continuity equation

[TABLE]

Our goal is to derive a reduced model for a given sequence of shells variables, $u_{1},\ldots,u_{s}$, where $s$ is any shell number from the inertial interval. The latter means that the viscous term in Eq. (1) can be neglected for the corresponding shells with

[TABLE]

The reduced probability distribution is defined as the result of integration over all shells with $n>s$:

[TABLE]

Similar integration applied to Eq. (5) yields

[TABLE]

The terms with the derivatives $\partial/\partial u_{n}$ for $n=s+1,\ldots,N$ vanish after the integration with respect to $u_{n}$. For the other terms, we write $P=(P/P_{s})P_{s}$ and substitute in (6). The resulting equation becomes

[TABLE]

where

[TABLE]

Here we specified that the functions $R_{n}$ may depend on all shell variables $u_{1},\ldots,u_{s}$ and time, due to the corresponding dependence of $P$ and $P_{s}$. The key point is to realize that equation (9) describes the evolution of the probability density for a reduced dynamical system:

[TABLE]

Eq. (11) will be our coarse grained system, when we obtain closed expressions for the right-hand sides in (10) as functions of the variables $u_{1},\ldots,u_{s}$.

The very same approach can be followed for the full Navier–Stokes equations, see (Pope, 2001, Chap. 13.5.6). The main advantage given by shell models is that they have only local or quasi-local interactions among consecutive shells. Indeed, the factor $Q_{n}=Q(u_{n-1},u_{n},u_{n+1})$ does not depend on the integration variables in (10) for those shells with $n<s$, while $Q_{s}=Q(u_{s-1},u_{s},u_{s+1})$ depends on $u_{s+1}$. Hence, the integration in (10) can be carried out using (7) and leading to the explicit expressions

[TABLE]

Here $P_{s+1}(u_{1},\ldots,u_{s+1};t)$ is defined by the expression analogous to (7). We see that the original system (6) and the reduced system (11), (12) differ only by the last equation. This is natural because the nonlinear term $Q_{s}$ is the only one that depends on the unknown shell variable $u_{s+1}$. Thus, the only missing component of the reduced system is the unknown integral expression in Eq. (12). In the jargon of LES the sub-grid model is influencing the explicit dynamical evolution of only one resolved variable (but still depends on the correlations with all of them).

In general, one needs to know the whole distribution $P_{s+1}(u_{1},\ldots,u_{s+1},t)$ to compute $R_{s}$ in (12). The main idea of this paper is that the form of the function $R_{s}$ is in fact universal in the developed turbulent dynamics, as suggested by numerical simulations and some theoretical considerations described below. This observation is central for our work and provides the subgrid model (12) in closed form.

III Kolmogorov’s third hypothesis and universality of the reduced equations

In 1962, Kolmogorov Kolmogorov (1962) conjectured that the statistics of velocity increment ratios (multipliers) $\delta_{\ell}v/\delta_{\ell^{\prime}}v$ is universal and depends only on the scale ratio $\ell/\ell^{\prime}$ in the inertial interval of homogeneous isotropic hydrodynamic turbulence. This conjecture, called the third Kolmogorov hypothesis, was confirmed both numerically and experimentally Frisch (1995); Chen et al. (2003); L’vov and Procaccia (1996); Benzi et al. (1998); Chhabra and Sreenivasan (1992); Pedrizzetti et al. (1996); Jouault et al. (1999); Nelkin and Stolovitzky (1996). For shell models, this conjecture implies that the probability distribution of multipliers $z_{n}=u_{n}/u_{n-1}$ is universal and does not depend on $n$ in the inertial interval, which agrees very well with numerical simulations for the Sabra shell model Eyink et al. (2003). Furthermore, Kolmogorov assumed that the multipliers for widely separated shells are statistically independent. Indeed, the distribution of multipliers appears to be short-range, i.e., correlations between $z_{n}$ and $z_{n+j}$ decay rapidly with increasing $j$.

The factor in the integral expression (12),

[TABLE]

is by definition the conditional probability of $u_{s+1}$ for given $u_{s},\ldots,u_{1}$ at time $t$. Note that there is a one-to-one correspondence between the shell variables $u_{1},\ldots,u_{s}$ (with given boundary condition for $u_{0}$) and the multipliers $z_{1},\ldots,z_{s}$ in the case when all of them are nonzero. The singular subset, when one of the variables vanishes, has zero measure and it is not important for our probabilistic analysis. Similarly, there is one-to-one correspondence between the shell variable $u_{s+1}$ and the multiplier $z_{s+1}$ for given $u_{1},\ldots,u_{n}$ (or $z_{1},\ldots,z_{n}$). Hence, the change of variables from $u_{n}$ to $z_{n}=u_{n}/u_{n-1}$ yields the conditional probability for shell variables in terms of the conditional probability for multipliers as

[TABLE]

The third Kolmogorov hypothesis for the developed turbulent regime implies that the function

[TABLE]

is universal and time-independent, such that it is uniquely determined for a given shell model. Also, $\widetilde{P}_{\textrm{uni}}$ must have short-range dependence on its arguments, i.e., it depends essentially only on a few neighboring shells $z_{s},z_{s-1},\ldots$ with very weak dependence on $z_{n}$ for smaller $n$. From now on, we will use the arguments written as $z_{s},z_{s-1},\ldots$ to indicate such a short-range dependence.

Note that Eqs. (14) and (15) do not necessarily imply the universality of $P_{\textrm{cond}}(u_{s+1}|u_{s},\ldots,u_{1};t)$. This is because the correlations between the shell variables $u_{n}$ extend to the whole range of scales and, hence, $P_{\textrm{cond}}$ may depend on the first (large scale) shells and on the boundary conditions.

Using Eqs. (13)–(15) and (2) in (12), yields

[TABLE]

Since $Q$ in (2) is a quadratic function of its arguments, we can write

[TABLE]

Using this expression in (16) provides the final expression

[TABLE]

[TABLE]

The universality of $\widetilde{P}_{\textrm{uni}}$ automatically implies the universality of the reduced system function $R_{s}$. Furthermore, the expressions show that $R_{s}$ is a homogeneous function of its arguments of degree $2$, just like the original nonlinearity $Q_{s}$. Short-range dependence of $\widetilde{P}_{\textrm{uni}}(z_{s+1}|z_{s},z_{s-1},\ldots)$ on its arguments leads to the similar property for $R_{s}(u_{s},u_{s-1},\ldots)$: this function depends essentially on a few variables $u_{s},u_{s-1},\ldots$, while the dependence on $u_{n}$ becomes very weak with decreasing $n$.

We arrived to the important and rigorous conclusion that the third Kolmogorov hypothesis yields the universal law (18), (19) describing the dynamics of the last shell $u_{s}$ in the system (11), (12). This deterministic dynamical system governs the evolution of the reduced probability density $P_{s}$ in Eq. (9) in the developed turbulent regime.

Note that the standard Large Eddy Simulations (LES) formulation of the Navier–Stokes equations involves modeling of the (effective) turbulent eddy viscosity Meneveau and Katz (2000); Pope (2001); Sagaut (2006); Lesieur et al. (2005). Such viscosity can be introduced explicitly in terms of velocity field or defined with renormalization techniques, see e.g. Zhou (2010); Verma and Kumar (2004). For example, it is defined in terms of the rate-of-strain tensor in the Smagorinsky model Smagorinsky (1963). As a result, the effective viscous term is a homogeneous function of degree 2 in the velocity field. Such an observation puts our closure for shell models in direct relation with the LES approach: the homogeneous function $R_{s}$ can be seen as a generalized term that includes the turbulent eddy viscosity. The idea of our work is to go beyond the concept of effective viscosity focused on the process of energy dissipation, by modeling a closure that describes the actual statistics at subgrid shells. As such, the proposed analysis of subgrid closures in shell models becomes a useful theoretical tool for testing optimal strategies, which may be potentially extended to the LES schemes for the NS equations.

IV Application to the Desnyansky–Novikov shell model

Typical models with simple first-neighbor coupling as (2) develop a non-chaotic (non-turbulent) behavior. However, we still can use such models for demonstrating basic principles of the reduction, before considering more sophisticated models in the next section. Let us consider the Desnyansky–Novikov model Desnyansky and Novikov (1974); Bell and Nelkin (1977) defined by the nonlinear term of the form

[TABLE]

The model possesses the energy $E=\frac{1}{2}\sum u_{n}^{2}$ as an inviscid invariant and has a non-chaotic time evolution. Solutions of equations (1) develop a power-law tail

[TABLE]

in the inertial interval, Fig. 1. Here $a>0$ is an arbitrary factor generally depending on time, which is related to the energy flux from large to small scales. Such a tail can be interpreted as a shock wave for a continuous representation of the model Mailybaev (2015), in close analogy with the Burgers equation.

Though the shell model dynamics is regular, the Kolmogorov hypothesis holds in the inertial interval. As follows from Eq. (21), all the multipliers

[TABLE]

are constant. The corresponding universal probability density becomes the Dirac delta-function as

[TABLE]

Relation (19) with $Q$ from (20) yields

[TABLE]

Finally, we find the reduced system function from (18) as

[TABLE]

The reduced system (11), (12) becomes

[TABLE]

For a numerical test, we consider the Desnyansky–Novikov model with $n=20$ shells, viscosity $\nu=10^{-7}$ and wavenumbers $k_{n}=\lambda^{n}$ with $\lambda=2$. As a boundary condition, we take $u_{0}=0$ and consider a decaying solution from the initial data $u_{n}(0)=k_{n}^{-1/3}$. The solution is shown by black lines with crosses in Fig. 1. It has viscous range around the shell $17$ and the power-law dependence (21) in the inertial interval $4\lesssim n\lesssim 15$. The reduced system (26) was integrated with $s=8$ shells, and the numerical results are presented in Fig. 1 by red circles demonstrating an excellent match with the full model solution.

V Application to the Sabra shell model

It is straightforward to extend the results of Sections II and III to more general shell models, where shell variables are complex numbers as in the GOY or Sabra shell models Gledzer (1973); Ohkitani and Yamada (1989); L’vov et al. (1998) or defined in term of sets of variables (vectors) Eggers and Grossmann (1991); Benzi et al. (1996); De Pietro et al. (2015); Lessinnes et al. (2009). Equations of motion for such models have the same structure (1), but the nonlinear term $Q_{n}$ may depend on several shells from each side.

In this section, we formulate the reduction for the Sabra shell model L’vov et al. (1998), which is characterized by complex shell variables $u_{n}\in\mathbb{C}$. The nonlinear term $Q_{n}$ in (1) describes the interaction with two neighbors given by

[TABLE]

where the shells $u_{0}$ and $u_{-1}$ must be specified by boundary conditions. The choice $c=-1/\lambda$ and $b=-1-c$ corresponds to the so-called 3D regime and leads to the two inviscid invariants: the energy $E=\frac{1}{2}\sum|u_{n}|^{2}$ and the helicity $H=\frac{1}{2}\sum(-1)^{n}k_{n}|u_{n}|^{2}$.

The probability density of this system evolves under the continuity equation written similarly to (5) as

[TABLE]

where we denoted $u_{n}=a_{n}+ib_{n}$. The analogous derivation as in Section II yields the description for the reduced probability distribution $P_{s}(u_{1},\ldots,u_{s})$ in the form

[TABLE]

where $s$ is a fixed shell number from the inertial interval, and the functions $A_{n}+iB_{n}=R_{n}(u_{1},\ldots,u_{s})$ are defined below. For complex variables $u_{n}$, this equation corresponds to the evolution of probability density for the reduced dynamical system

[TABLE]

Similarly to Eqs. (12) and (13) one derives

[TABLE]

where the first $s-2$ equations of the shell model remain unchanged. The last two equations are given by the conditional averages

[TABLE]

Again, our construction shows that the subgrid scheme is given by a deterministic system of equations. The time-independence and universality of this conditional average follows from the Kolmogorov hypothesis formulated for multipliers as we demonstrate below.

Let us introduce the complex multipliers $z_{n}\in\mathbb{C}$ as Benzi et al. (1993); Eyink et al. (2003)

[TABLE]

Here the phases $\Delta_{n}$ are chosen to be invariant under the phase symmetry

[TABLE]

which is an analog in the Sabra model of the physical space homogeneity L’vov et al. (1998). The combination of phases given by $\Delta_{n}$ is important because it is strictly connected to the existence of a mean forward energy cascade (see expression (48) and discussion thereof). It is easy to see that there is one-to-one correspondence between the multipliers $z_{1},\ldots,z_{s}$ and the shell variables $u_{1},\ldots,u_{s}$ with $u_{0}$ and $u_{-1}$ given by boundary conditions (except for a zero-measure subset when some $u_{n}=0$). Thus, an argument similar to the one used in Eqs. (14) and (15) can be applied to the conditional probability $P_{\textrm{cond}}$ in (32). Namely, one can use the third Kolmogorov hypothesis for expressing this function as

[TABLE]

where we denoted $z_{n}=x_{n}+iy_{n}$ and $\widetilde{P}_{\textrm{uni}}$ is a universal function describing the conditional distribution for multipliers in the developed turbulent regime. The function $\widetilde{P}_{\textrm{uni}}$ is expected to have a short-range dependence on its arguments, i.e., the dependence on $z_{n}$ for $n\leq s$ gets very weak with decreasing $n$. The universality property was thoroughly studied numerically in Benzi et al. (1993); Eyink et al. (2003).

According to (31) the reduced system (30) contains the two unknown functions on the right-hand sides, $R_{s-1}(u_{s},u_{s-1},\ldots)$ and $R_{s}(u_{s},u_{s-1},\ldots)$. Using the explicit form (27) of the nonlinear term in the conditional average (32) for $n=s-1$ yields

[TABLE]

where we used (33) and the fact that the conditional average for given $u_{s},u_{s-1},\ldots$ is equal to the conditional average for given $z_{s},z_{s-1},\ldots$. Similar computation for $R_{s}$ yields

[TABLE]

Combining (36) and (37), we write the final expressions for the reduced system as

[TABLE]

where

[TABLE]

By construction, $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ are functions of $z_{s},z_{s-1},\ldots$ or, equivalently, homogeneous functions of zero degree with respect to shell variables $u_{s},u_{s-1},\ldots$. As a result, the functions $R_{s-1}$ and $R_{s}$ in (38) are homogeneous of degree 2 with respect to shell variables $u_{n}$. Note that the multipliers in our approach are used for justification of the universality, but the final expressions can be written in terms of the original shell variables $u_{n}$.

V.1 Simple subgrid models

For computing the statistics of multipliers, we performed a single long-time simulation with the inter-shell ratio $\lambda=2$ and the viscosity $\nu=10^{-12}$ (the viscous range starts around $n_{K}\approx 30$) and constant boundary conditions with $u_{0}=u_{-1}/2=1+i$. All the results in this section are obtained using $10^{6}$ samples $(u_{1},\ldots,u_{N})$ equally spaced in time. Fig. 2 (a,b) shows the PDFs for the absolute values $w_{n}=|z_{n}|$ and phases $\Delta_{n}=\arg z_{n}$ for different $n$, confirming the scale-invariance of the distribution for multipliers in agreement with earlier results Benzi et al. (1993); Eyink et al. (2003). Fig. 2 (c,d) demonstrates the corresponding correlation functions confirming the short-range property, i.e., a rapid decay of correlations with the shell separation. The correlation is larger for phases than for absolute values. The main consequence is that, with such rapidly decaying correlations, one may hope to obtain a reasonably accurate approximate subgrid model by keeping very few multipliers in the conditional averages (39).

In this section, we introduce three low-order approximations for the functions (39) that uniquely determine the corresponding subgrid models given by (30), (31) and (38). The first model will be denoted by $\mathrm{SM_{K}}$ and called the Kolmogorov closure. It is given by the probability function

[TABLE]

Here the product of Dirac delta-functions simply means that both multipliers have deterministic absolute values $w_{s+2}=w_{s+1}=\lambda^{-1/3}$ according to the Kolmogorov scaling law, and their phases $\Delta_{s+2}=\Delta_{s+1}=\pi/2$ are fixed at the most probable values, see Fig. 2(b). In this case, expressions (39) yield

[TABLE]

The next model, denoted by $\mathrm{SM_{0}}$, is considered to be a zero-order approximation based on numerical information obtained from the unclosed original equations. Namely, we take the values

[TABLE]

by supposing to consider the minimal degree of correlation in (39). In particular, the value for $\mathcal{A}_{1}$ has been estimated from an unconditional average of the full viscous unclosed model, while for $\mathcal{A}_{2}$ we imposed the minimal constraints by fixing the first resolved multiplier to its Kolmogorov value. The conditioning in the latter expression is necessary, because the unconditional average $\langle|z_{s+1}|^{2}z_{s+2}\rangle$ would diverge. The divergence is related to shell variables $u_{s}$ passing close to the origin, which leads to large $|z_{s+1}|=|u_{s+1}/u_{s}|$ simultaneously with small $|z_{s}|=|u_{s}/u_{s-1}|$, see Eyink et al. (2003). Thus, this defect can be avoided by excluding the events with small $z_{s}$. Note that the numerical values in (42) are quite different from the Kolmogorov prediction in (41), caused by cancellations due to a large spread of the phases.

Finally, we consider the model denoted by $\mathrm{SM_{1}}$, which is obtained as a first-order approximation by considering the averages conditioned to the single multiplier $z_{s}=|u_{s}/u_{s-1}|e^{i\Delta_{s}}$ in (39), i.e.,

[TABLE]

The averages are determined numerically from the results of a direct numerical simulation of the full unclosed equations as mentioned above. Due to the “defect” of the averages at the origin, it is convenient to use the ansatz

[TABLE]

Here, the variable $0\leq\varphi_{s}<\pi/2$ is used for representing the infinite semi-interval $w_{s}\geq 0$. The functions $f(\varphi_{s},\Delta_{s})$ and $g(\varphi_{s},\Delta_{s})$ are universal in the inertial interval. We compute them by fitting the simulation data with the Fourier expansion in $\Delta_{s}$, where coefficients are polynomial functions of $\varphi_{s}$, see Fig. 3. The averaged data are obtained by taking mean values of the simulation results in every grid cell $(\delta\Delta_{s},\delta\varphi_{s})$ in Fig. 3. This yields

[TABLE]

with

[TABLE]

In these approximations, we kept the Fourier modes for $\Delta_{s}$, which had the dominant contribution, while the coefficients depending on $\varphi_{s}$ were approximated with low-order polynomials neglecting the coefficients smaller than $0.02$. All functions in (46) appear to be real due to the inherent symmetries of the Sabra model. Fig. 3 demonstrated the comparison of expressions (45) with the same functions found numerically by averaging expressions (43) in the inertial interval.

V.2 Numerical tests

The numerical tests are carried out for the three subgrid models from the previous section with $s=15$ and $20$ shells. We do this in the time interval $0\leq t\leq 10^{3}$ with the constant boundary conditions

[TABLE]

For initial conditions we take the Kolmogorov state $u_{n}=k_{n}^{-1/3}e^{i\phi_{n}}$ with random phases $\phi_{n}$. The comparison is also made with the simulation of the full viscous model (1), (27) with total $N=40$ shells and viscosity $\nu=10^{-12}$ (the viscous range starts around shell $n_{K}\sim 30$).

Simulation results for some representative interval of time are compared in Fig. 4. One can see that the dynamics at large scales (first row) is qualitatively similar for all models, but small scales (second row) demonstrate some qualitative differences. In particular, we notice a tendency to lock among the last three variables much more pronounced than in the unclosed case, an indication that phase correlation is probably not fully correct. Note that obviously we do not expect a detailed correspondence of the solutions, since subgrid models are designed to describe a probability distribution rather than a particular solution.

The differences among the models can be seen more clearly and systematically in Fig. 5 presenting the time-averaged energy spectra:

[TABLE]

One can see that, indeed, large deviations are observed in the region of large shell numbers $n\sim s$ (small scales near cutoff), while good agreement is attained at smaller $n$ (larger scales). Surprisingly, more elaborated subgrid models do not show improved results, with a better match given by the simplest SMK model. We postpone a discussion on why this happens to the next section, and concentrate on describing different aspects of solutions for different models now. From Fig. 5 one can see that the deviations demonstrated by every model do not depend on the total number of the shells $s$: they repeat the same pattern for $n\sim s$ that converges to the full model results for smaller $n$. Similar results are observed for other velocity moments as well. In particular, Fig. 6 presents the results for the flatness:

[TABLE]

demonstrating similar type of discrepancies at final shell numbers.

Though showing rather large deviations for average values, the subgrid models describe probability of non-Gaussian rare events (intermittency) reasonably well. Fig. 7(a,b) shows the normalized PDFs in the logarithmic vertical scale, where fat tails are well reproduced. Still, also in this case, the models SM1 or SM0 are less accurate than a simple model SMK based on the Kolmogorov closure. Fig. 7(c) provides the comparison of PDFs for the phase variable $\Delta_{s}$ at the last subgrid shell number $s=20$, demonstrating a considerable variation among different models. Note that the phases $\Delta_{n}$ determine the direction of energy flux, see Eq. (48) below.

The next test is related to the energy flux. Recall that the energy $E=\frac{1}{2}\sum|u_{n}|^{2}$ is an inviscid invariant for the Sabra model. A general expression for the energy flux across shell $n$ is given by $\Pi_{n}^{E}=\mathrm{Im}\,(k_{n+1}u_{n+2}u_{n+1}^{*}u_{n}^{*}-ck_{n}u_{n+1}u_{n}^{*}u_{n-1}^{*})$, see e.g. Biferale (2003); Ditlevsen (2010); Kumar and Verma (2015). With the phases $\Delta_{n}$ from (33), this expression is written as

[TABLE]

Expression (48) determines the non-linear contributions to the total change of the energy in the shells up to $n$. Thus, for inviscid dynamics, the energy balance takes the form

[TABLE]

where $\Pi_{0}^{E}$ is the work per unit time done by the boundary due to nonzero values of the shells $u_{0}$ and $u_{-1}$. For subgrid models, expression for the energy flux is different for the modified shells $n=s-1$ and $s$. Direct computations for these fluxes using (49) and (30), (31), (38) yields

[TABLE]

[TABLE]

In particular, for the SMK model given by (41) with $\lambda=2$ and $c=-1/2$, one obtains the strictly positive flux at last shell as

[TABLE]

Similarly, the inviscid invariant called helicity is introduced as $H=\frac{1}{2}\sum(-1)^{n}k_{n}|u_{n}|^{2}$. In this case the helicity flux of the inviscid model is given by the expression

[TABLE]

The time-averaged values of the energy and helicity fluxes computed at shell corresponding to the small-scale range $n\ll s$ are in very good agreement with the unclosed model. For example, in Fig. 8 we compare the PDFs of energy fluxes for different models, computed at the large-scale shell $n=3$ and for the last two shells $s-1=19$ and $s=20$. One can see that a very good convergence is attained for large scales, but again the some discrepancies are observed for the small scales. Here the models SMK and SM0 show no energy backscattering events with a strictly positive flux as predicted by Eq. (52). On the other hand, the full model and the model SM1 show energy backscattering events (negative energy flux).

For interpreting this result, it is necessary to recall that the subgrid models were designed for the statistical distribution averaged over the shell numbers $n>s$. Thus, the reference quantity for checking the validity of subgrid models must be the energy flux averaged in the same manner. We performed such an averaging of the energy flux numerically for the full viscous model: using the long-time simulation results, the flux function $\Pi_{s}^{E}$ was averaged over the points with nearby values of shell speeds $u_{s}$ and $u_{s-1}$ but arbitrary $u_{s+1}$ and $u_{s+2}$. The PDF for the resulting (conditionally averaged) flux is shown in Fig. 8(c) by bold grey dotted line. This results shows a rather surprising result that, after averaging over shells $n>s$, the energy flux $\Pi_{s}^{E}$ becomes strictly positive. This also means that the models SMK and SM0 have the correct behavior, while the backscattering events demonstrated by the model SM1 (green line in Fig. 8(c)) should not be interpreted in favor of this subgrid model.

It is also instructive to compute the flatness in terms of the energy flux as

[TABLE]

The corresponding numerical results are given in Fig. 6 (c), demonstrating a slightly more regular behavior. Deviations are still present at cutoff shells, but this time the model SM0 shows a better match with the full model.

V.3 Why improved subgrid models do not work better?

In subgrid models discussed in Sections V.1 and V.2, we used direct averages for equations of the last shells, or the average conditioned on one multiplier $z_{s}$. Though the closures constructed in this way are rather simplistic, the very fast decay of correlations of multipliers with a shell separation suggests that even such simple models should be reasonably accurate, see Fig. 2 (c,d). On the contrary, we saw that more accurate models do not demonstrate any improvement for the statistics at small scales. In order to verify this observation with higher-order approximations, we also constructed the subgrid model that depends on two multipliers, $z_{s}$ and $z_{s-1}$. This was done by using expansions in multi-dimensional spherical harmonics, which is possible due to a homogeneity property of $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ as functions of shell speeds $(u_{s},u_{s-1},u_{s-2})$. Details of this study are rather lengthy and will not be presented here. The conclusion is, however, the same: no considerable improvement was observed for the statistics at small scales, in comparison with the simplest model SMK. In this section we propose a possible explanation to this unfortunate lack of convergence.

The distribution limited to a finite number of shells in the inertial region, $P(u_{1},\ldots,u_{s})$, can be expected to be a smooth positive function, as one can infer from Kolmogorov’s third hypothesis (see Section III). Our subgrid models are constructed with the purpose of recovering this regular distribution in a statistical sense, i.e., as an attractor. However, the regularity of the distribution imposes a very strong requirement on the subgrid model, which is reminiscent to the mixing property for measure-preserving dynamical systems: an infinitely long trajectory of the subgrid model must be dense everywhere in the configuration space. As we know from the dynamical system theory, such a property is structurally unstable for a non-conservative dynamical system, like our subgrid model. This implies that an arbitrarily small change of the “ideal” subgrid system may drastically change its long-time statistical behavior; see also Berkooz (1994), where this issue was investigated for the Lorenz system. In a general case, one can expect that the subgrid model possesses a chaotic (fractal) attractor, therefore, occupying only a zero measure subset in configuration space. This structural instability may be the main cause of the persistent divergence from the full model statistics in our subgrid models. It is important to notice that this high sensitivity to the closure is nevertheless limited to a fixed (Reynolds independent) number of shells, indicating that the large scale dynamics is robust and universal with respect to the small-scale closure, i.e. we do not have a strong sensitivity of the global attractor on the fine details of high frequency fluctuations.

It is possible to demonstrate some quantitative evidence in favor of our hypothesis. The inviscid Sabra model has unstable time-independent solutions of Kolmogorov type, which up to phase-symmetry factors have the form

[TABLE]

with arbitrary period-3 real coefficients $a_{n}$, see e.g. Biferale et al. (1995); Mailybaev (2016). In this case, all multipliers (33) are fixed numbers, $z_{n}=-i\lambda^{-1/3}a_{n}/a_{n-1}$, similarly to the Desnyansky–Novikov shell model in Section IV. For our subgrid models, the factors $a_{n}$ are not arbitrary and can be found by substituting (54) into the equations $R_{s-1}=0$ and $R_{s}=0$, see (30) and (31). The elementary computations with expressions (38) yield

[TABLE]

For the two simplest models, where $\mathcal{A}_{1,2}$ are constants (see Section V.1), we have

[TABLE]

Figure 9 shows the PDFs for the absolute values of multipliers, $w_{s}=|z_{s}|$ and $w_{s-1}=|z_{s-1}|$, in these two models obtained from numerical simulations. Blue arrows indicate the positions given by the time-independent solutions (54)–(57) as $w_{s}=\lambda^{-1/3}|a_{s}/a_{s-1}|$ and $w_{s-1}=\lambda^{-1/3}|a_{s-1}/a_{s-2}|$. One can clearly see the correlation between the PDFs and time-independent solutions, which can be explained by the intermittent dynamics alternating between the chaotic and regular behavior, see, e.g., Ott (2002). One can also see some footprints of this temporary “locking” to a time-independent solution for the model SMK in Fig. 4 (f), where a small plateau is developed with slowly changing amplitudes.

We conclude that the chaotic attractor may be influenced by the time-independent solution at small scales. As we explained, this defect is generic and, hence, hardly can be removed by using more accurate subgrid models, unless some special extra conditions are imposed. For example, noise can be used as a mechanisms to improve the dynamics at subgrid scales, by counteracting the attraction to specific solutions. In order to see how large is the effect, we applied random phase perturbations $\Delta_{s+1}=\pi/2+x_{1}$ and $\Delta_{s+2}=\pi/2+x_{2}$ in the models SMK, see Section V.1. Here $x_{1}(t)$ and $x_{2}(t)$ are obtained as solutions of the Langevin equation $\dot{x}=-\tau^{-1}x+\sigma\tau^{-1/2}\xi$ with the Kolmogorov time scale $\tau=k_{s}^{-2/3}$ and white noise $\xi$ with $\langle\xi(t)\rangle=0$ and $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$. To have a moderate noise level we choose $\sigma=0.4$. Figure 10 shows the PDFs of the multipliers $w_{s}$ and phases $\Delta_{s}$ at the last shell of the model with $s=20$. One can see that the noise has some effect improving the distributions, as compared with the bold red PDFs of the full model, but this effect is small. The results suggest that adding an uncorrelated noise to the evolution of phase is not an effective mechanism for our subgrid model.

VI Conclusions

We have discussed a theoretical framework to define the optimal sub-grid closure for shell models of turbulence. The theoretical framework would predict a very complicated sub-grid models which depends on the conditional probability of sub-grid variables on all resolved scales, a task unrealistic even for the simple structure of shell models. We have proposed a series of approximate closures based on the ansatz that consecutive shell multipliers are short-range correlated, following the third hypothesis of Kolmogorov formulated for similar quantities for the three-dimensional Navier–Stokes turbulence. Different approximations assume different degrees of correlations across scales among amplitudes and phases of consecutive multipliers. We show numerically that such low order closures work well, reproducing all known properties of the large-scale dynamics including anomalous scaling. We found small but systematic discrepancies only for a range of scales close to the sub-grid model, which do not tend to disappear by increasing the order of the approximation. We speculate that the lack of convergence might be due to a breaking of ergodicity at least for the evolution of very fast degrees of freedom at small scales. Effects of the sub-grid closure on the resolved range of scales must be quantified also for real LES of three-dimensional Navier–Stokes equations. Correlations between the sub-grid stress tensor and velocity increments at the resolved scales can be estimated on the basis of fusion-rules L’vov and Procaccia (1996); Benzi et al. (1998). They are supposed to be sub-leading with respect to the scaling of single-scale velocity increments Meneveau (1994), i.e. reproducing the same kind of sensitivity to the particular closure only for a range of separations close to the sub-grid cutoff. A quantitative assessment of the importance of such a feedback is nevertheless missing.

Acknowledgments. The authors are grateful to R. Benzi, M. Cencini and G.L. Eyink for useful comments. The work was supported by the European Research Council under the European Union’s Seventh Framework Programme, ERC Grant Agreement No 339032. A.A.M. was also supported by the CNPq (Grant No. 302351/2015-9) and FAPERJ (grant E-26/210.874/2014).