Inertial range statistics of the Entropic Lattice Boltzmann in 3D turbulence

Turbulence is common in nature and its unpredictable behavior has fundamental consequences on the understanding and control of various systems, from smaller engineering devices dumitrescu2004rotational ; Davidson2015 ; Pope2001 , up to the larger scales geophysical and astrophysical flows gill2016atm ; alexakis2018 ; barnes2001 . Turbulent flows are described by the Navier-Stokes equations (NSE),

$$\begin{cases}\partial_{t}\bm{u}+\bm{u}\cdot\bm{\nabla}\bm{u}=-\bm{\nabla}p+\nu\nabla^{2}\bm{u}+\bm{f}\\ \bm{\nabla}\cdot\bm{u}=0\end{cases}$$

(1)

which give the evolution of the incompressible velocity field $\bm{u}(\bm{x},t)$, with kinematic viscosity $\nu$, subject to a pressure field $p$ and to an external forcing $\bm{f}$. However, even though the equations of motion are known since almost two hundred years a direct analytical approach remains elusive Frisch1995 . To overcome mathematical difficulties, scientists, helped by the exponential growth of the computational power, have tried to search for approximate solutions using numerical algorithms fox2003computational ; Pope2001 ; ishihara2009study . Unfortunately, also in this direction not every effort were successful. Indeed, the NSE have a very rich non-linear dynamics, where a large range of scales from the domain size up to the small scales fluctuations, are coupled together. This results in a very high dimensional problem, with the dimensionality proportional to the range of active scales Frisch1995 , and with highly intermittent statistics dominated by the presence of extreme and rare fluctuations yeung2015extreme ; benzi2010inertial ; iyer2017reynolds . As a consequence, no matter how powerful new-supercomputers are, numerical algorithms cannot handle all the degrees of freedom involved in the dynamics Davidson2015 . The way out is to introduce a scale separation and compute only the dynamics of degrees of freedom belonging to a subset of scales while neglecting the other ones Meneveau2000 ; filippova2001multiscale ; Dong2008a ; Dong2008b . However due to the non-linearity of NSE there is never a real scale separation in the equations of motion bohr2005dynamical ; Frisch1995 , and the small-scale effects on the scales of interest need to be compensated by the introduction of a model. In other words, the benefit of multi-scale modeling is to achieve a scale separation, and the main challenge is to find a “closure”, which guaranties the numerical stability being at the same time the most accurate as possible in reproducing the coupling of the missing scales on the resolved ones. This is the principle behind the celebrated Large-Eddy simulations (LES), which actually solve the flow only on a subset of “large” scales by filtering each term of the NSE and replacing with a closure the non-linear coupling term between the resolved and the sub-grid scales (SGS) Pope2001 ; Lesieur2005 . One of the most important differences between the real coupling term coming from the filtered NSE and the common LES closures used in the literature is for the latter to be purely dissipative to ensure the simulation stability, see Meneveau2000 . As a consequence, it is impossible for the closures to reproduce the backscatter events of energy going from the SGS to the large scales, with important consequences on the statistics of the resolved velocity field. Another possible numerical approach, who has gained particular popularity, consists in solving the flow’s macroscopic hydrodynamical properties as an approximation of its mesoscopic behaviour de2013non . The Lattice Boltzmann Method (LBM) falls into that category benzi1992lattice ; Succi2001 . In LBM, the flow is simulated by evolving the Boltzmann equation for the single phase density function, $f(\bm{x},t)$. The idea is to evolve the streaming and collision of particles distribution functions, where the possible velocities are restricted on a subset of discrete lattice directions frisch1986lattice ; Wolf-Gladrow2000 . It is crucial to choose the collision operator so that macroscopically, (in the limit of small Knudsen number), the dynamics described by the NSE is recovered qian1995recent ; chen1998lattice . The most common collision operator is the Bhatnagar-Gross-Krook model (BGK), see Bhatnagar1954 , corresponding to the relaxation towards an equilibrium distribution, $f^{eq}_{i}(\bm{x},t)$, taken to be a discrete Maxwellian, with a fixed relaxation time, $\tau_{0}$,

$\displaystyle\begin{split}f_{i}(\bm{x}+\bm{c}_{i}\Delta_{t},t+\Delta_{t})-f_{i}(\bm{x},t)&=\\ =-\frac{1}{\tau_{\rm 0}}\left[f_{i}(\bm{x},t)-f^{eq}_{i}(\bm{x},t)\right]&+F_{i}.\end{split}$

(2)

Here, $F_{i}$, is a term introduced to model a macroscopic external forcing Succi2001 , and $i=0,..,q-1$, indexes the different velocity directions on the lattice. Eq. (2) is obtained discretizing the Boltzmann equation and selecting the lattice spacing, $\Delta\bm{x}$, such as divided by the time step, $\Delta_{t}$, they are equal to the lattice velocity $\bm{c}_{i}=\Delta\bm{x}/\Delta_{t}$. From those mesoscopic quantities, it is then possible to recover the macroscopic velocity and density fields by following the perturbative Chapman-Enskog expansion Wolf-Gladrow2000 . It can be shown that evolving Eq. 2 is equivalent, up to approximations, to evolving the weakly compressible NSE for a flow with a density $\rho(\bm{x},t)=\sum_{i=0}^{q-1}f_{i}(\bm{x},t)$, a velocity $\bm{u}(\bm{x},t)=\sum_{i=0}^{q-1}f_{i}(\bm{x},t)\bm{c}_{i}/\rho(\bm{x},t)$ and a viscosity directly related to the relaxation time Wolf-Gladrow2000 ,

$$\nu_{0}=c_{s}^{2}\Delta_{t}\left(\tau_{0}-\frac{1}{2}\right)\text{,}$$

(3)

where, $c_{s}$, is the speed of sound, and $\tau_{0}$ is the adimensional relaxation time. A numerical validation of the hydrodynamic recovery of BGK-LBM in the context of 2D homogeneous isotropic turbulence (HIT) was performed in Tauzin2018 showing good agreement with DNS either in decaying that in forced regimes. Although this method is adapted to describing various physics of multi-phase and flows with complex boundaries in a highly scalable fashion, the BGK-LBM model suffers of numerical instabilities when, $\tau_{\rm 0}\rightarrow\frac{1}{2}$, i.e. $\nu_{\rm 0}\rightarrow 0$, which has made the study of turbulent flows highly prohibitive for this method Tauzin2018 . To push LBM towards more turbulent regimes a number of collision operators have been proposed, see eggels1996direct ; yu2006turbulent ; sagaut2010toward ; malaspinas2012consistent . Here we focus on the Entropic LBM (ELBM) Karlin1999 ; ansumali2002single , which tackles the stability issue by equipping LBM with an H-theorem. To achieve this results the ELBM differs from BGK-LBM by two major aspects. First, the equilibrium distribution $\bm{f}^{eq}(\bm{x},t)$ is not anymore a discretization of the Maxwell-Boltzmann distribution, but it is calculated as the extremum of a discretized H-function defined as;

$$H[\bm{f}]=\sum_{i=0}^{q-1}f_{i}\log\left(\frac{f_{i}}{w_{i}}\right),\,\,\bm{f}=\left\{f_{i}\right\}_{i=0}^{q-1},$$

(4)

where $w_{i}$ are the weights associated to each lattice direction, under the constraints of mass and momentum conservation, see Ansumali2003 . The second difference in ELBM, is that the relaxation time is not a constant anymore but is modified at every time step in order to enforce the non-monotonicity of H after the collision. This results in an apparent unconditional stability as $\nu_{0}\rightarrow 0$ Karlin2015 . It follows that ELBM evolution equations are,

$\displaystyle\begin{split}f_{i}(x+c_{i}\Delta_{t},t+\Delta_{t})-f_{i}({\bm{x}},t)&=\\ =-\alpha({\bm{x}},t)\beta[f_{i}({\bm{x}},t)&-f_{i}^{eq}({\bm{x}},t)],\end{split}$

(5)

where $\beta=1/(2\tau_{0})$ is constant, while the new relaxation time, $\tau_{\text{eff}}(\bm{x},t)=1/(\alpha(\bm{x},t)\beta)$, fluctuates in time and space through the definition of an entropic parameter $\alpha(\bm{x},t)$. More recently the ELBM method has been extended to a family of multi-relaxation time (MRT) lattice Boltzmann models dorschner2016entropic ; dorschner2017transitional ; dorschner2018fluid . Note that $\alpha$ can be computed as the solution of the entropic equation $H[\bm{f}]=H[\bm{f}-\alpha\bm{f}^{neq}]$ which represents the maximum H-function variation due to a collision, with $\bm{f}^{neq}=\bm{f}-\bm{f}^{eq}$. Following this approach the computation of $\alpha(\bm{x},t)$ can be performed via an expensive Newton-Raphson algorithm for every grid and at every time step of ELBM. To alleviate this problem, after the original ELBM formulation Karlin1999 , a new version has been proposed where the computation of the entropic parameter is based on an analytical formulation derived as an first order expansion of the original model karlin2014gibbs ; bosch2015entropic . However, to our knowledge, a study of high-order structure functions in the context of forced 3D HIT, has never been attempted before using the ELBM original model. In this regard, and also aiming to measure high-order, extremely sensitive statistical observables, we implemented the ELBM original formulation, relying on the least number of approximations, even though computationally more expensive. More details about ELBM will be given in section II. Let us notice that BGK-LBM is recovered from Eq. (5) with $\alpha=2$ and the specific Maxwellian expression of $f_{i}^{eq}$. It is important to stress that the bridge relation described in Eq. (3) connecting viscosity and relaxation time still holds for fluctuating quantities, hence we can write,

$\displaystyle\begin{split}\nu_{\text{eff}}\left(\alpha\right)&=c_{s}^{2}\Delta_{t}\left(\frac{1}{\alpha\beta}-\frac{1}{2}\right)=\\ &=c_{s}^{2}\Delta_{t}\left(\frac{1}{2\beta}-\frac{1}{2}\right)+c_{s}^{2}\Delta_{t}\frac{2-\alpha}{2\alpha\beta}=\\ &=\nu_{0}+\delta\nu_{\alpha},\end{split}$

(6)

where $\nu_{0}$ represents the constant kinematic viscosity and $\delta\nu_{\alpha}$ is the fluctuating term. Following this idea, it has been shown that ELBM is implicitly enforcing a SGS model of an eddy-viscosity type Malaspinas2008 ; Karlin2015 . In particular, as initially done in Malaspinas2008 , and then rederived in chapter 4 of tauzin2019implicit , performing a third-order Chapman-Enskog perturbative expansion in the limit of small Knudsen number (Kn), it is possible to obtain a macroscopic approximation of $\delta\nu_{\alpha}$, which can be written as,

$$\delta\nu_{\alpha}^{M}=-c_{s}^{2}\Delta t^{2}\frac{S_{\ell j}S_{ij}S_{i\ell}}{S_{ij}S_{ij}},$$

(7)

where $S_{ij}=1/2(\partial_{j}u_{i}+\partial_{i}u_{j})$ is the strain rate tensor. The entropic eddy viscosity in Eq. (7) is particularly interesting because it is not positive definite and can reproduce events of energy backscatter, which is generally not the case among the other LES closures, i.e. see the Smagorinsky eddy viscosity Smagorinsky1963 . After the introduction of the new LES model, see section II, we compare it with standard Smagorinsky closure and with fully resolved data obtained from DNS. At the same time, we also present for the first time a quantitative investigation of the inertial range statistics provided by the mesoscopic ELBM approach in the context of 3d turbulence. Results provide evidence that ELBM is a good approximation of the 3d flows up to turbulent regimes never reachable with the standard BGK-LBM. We found that ELBM guaranties the simulation stability producing a considerable extension of the inertial range scaling, i.e. the extension of the energy spectrum power law. Measuring statistical properties of higher order observables such as the structure functions, we found that ELBM is also able to reproduce qualitatively the intermittent features of real 3d turbulent flows with an accuracy comparable to the Smagorinsky model. To get a more accurate estimation of the anomalous exponents more refined models are required biferale2019self .

The paper is organized as follows. In section II we introduce the details of the ELBM and the LES models considered in this work. In section  III we present the set of simulations performed. In sect. IV we evaluate the quality of the two LES closures by comparing them with the real SGS energy transfer measured from fully resolved DNS. In sect. V we focus on the intermittent properties by analysing high-order inertial range statistics Arneodo2008 ; benzi2010inertial . In sect. VI we discuss our conclusions.

In this section we give a description of the SGS modelling approaches considered in this work. We start discussing the mesoscopic ELBM approach highlighting the differences with respect to the standard BGK-LBM. Following we discuss the new hydrodynamic LES closure inspired by the ELBM macroscopic approximation first derived in Malaspinas2008 . In the end of this section we briefly recall the well known Smagorinsky model.

All simulations performed in this work are intended to model HIT on a three dimensional domain with periodic boundary conditions. In the following we provide some details about the sets of simulations performed with the different modelling techniques. Concerning the lattice Boltzmann simulation with entropic formulation of the relaxation time, ELBM, we have conducted a set of 3d simulations with a number of $512$ collocation points along each spatial direction. To reach stationarity the flow is forced at large scales, $1\leq|\bm{k}|\leq 2$ with a constant and isotropic forcing. More precisely we have used in all simulations the same forcing, defined in Fourier space with constant phases and amplitudes, added isotropically to all wave-vectors at large-scales. To ensure incompressibility the forcing is projected on its solenoidal component. The ELBM simulation uses a lattice with 27 discrete velocities (see Fig. 1), the D3Q27 Succi2001 ; Wolf-Gladrow2000 ; Kruger2017 . The spectral forcing is implemented using the exact-difference method forcing scheme Kuperstokh2004 for a relaxation time $\tau_{0}=0.5003$ corresponding to $\beta\approx 0.9994$.

Considering the LES we have performed two sets of pseudo-spectral fully dealiased simulations on a domain $\Omega=[0,2\pi]^{3}$ with periodic boundary conditions both at the resolution of $512^{3}$ grid points. The first LES is equipped with the Smagorinsky model (S-LES), see Eq. (18), and a second with the macroscopic formulation of the entropic eddy viscosity, see Eq. (19), (E-LES). As discussed above, all simulations are forced with the same isotropic constant forcing mechanisms acting only on the larger system scales ($1\leq|\bm{k}|\leq 2$). In the expression of Smagorinsky eddy viscosity Eq. (18), we use the standard value of $C_{S}=0.16$, while for the entropic eddy viscosity Eq. (19), we use $C_{E}=0.45$, found to be optimal values for the best compromise between the maximization of inertial range extension and the minimization of spurious effects produced by the model buzzicotti2020synch . For both we have $\Delta=\pi/k_{max}\approx 0.0184$ with $k_{max}=171$, which comes from the $2/3$ rule for the dealiasing projection Patterson71 . Additionally, as a reference, we have run a pseudo-spectral fully resolved DNS of the NSE with the same forcing scheme on the same 3d domain $\Omega=[0,2\pi]^{3}$, using a number of $512^{3}$ (DNSx1) and $1024^{3}$ (DNSx2) collocation points. The resolution in both DNS is kept such as $\eta_{\alpha}/dx\simeq 0.7$, where $dx$ is the grid spacing and $\eta_{\alpha}=(\nu^{3}/\epsilon)^{1/4}$ is the Kolmogorov microscale Borue95 with $\varepsilon$ denoting the mean energy dissipation rate. In order to create ensembles of statistically independent data all the ELBM, LES and DNS are sampled on time intervals of one large-scale eddy turnover time after reaching a statistically stationary state. It is worth mentioning that both the ELBM and E-LES simulations remain stable even though their modeling terms are not purely dissipative being their eddy viscosities not positive definite. In Fig. 2 we show the time-averaged energy spectra, $E(k)$, for all simulations. It is visible that all modelled simulations have an extended inertial range with respect to DNS at the same resolution (DNSx1) with an inertial range slope close to the Kolmogorov prediction of $k^{-5/3}$ Frisch1995 . Let us stress that the ELBM spectrum reaches a maximum wavenumber higher than the pseudo-spectral data, this is because in the ELBM simulations there is not any dealiasing operation. Anyway also in the ELBM case the spectrum loses the Kolmogorov inertial range scaling at wavenumber larger that the pseudo-spectral dealiasing cutoff.

In this section we provide a statistical comparison of the SGS energy transfers measured in the modelled simulations together with the original SGS transfer measured a priori from higher resolution DNS, see Table 1. Before staring the comparison between macroscopic and mesoscopic simulations, we analyse the quality of the approximation made in the Chapman-Enskog expansion to obtain the macroscopic formulation of the ELBM eddy viscosity Malaspinas2008 . In this direction, we have computed the SGS energy transfer defined in Eq. (21) using the two different definitions of the fluctuating viscosity. Namely, either the correct definition, $\delta\nu_{\alpha}$, depending on the entropic parameter or its third order expansion in the limit of small Kn, $\delta\nu_{\alpha}^{M}$ see Eq. (7). Their statistical comparison is shown in Fig. 3, where on the left panel we show the probability density functions (PDF) measured from the ELBM SGS energy transfer using the two different formulations. Here we can see that the PDFs once re-scaled by their standard deviations have almost an identical shape. From the center and right panels of the same figure, we can qualitatively see two visualizations of the SGS energy transfers measured by selecting the same plane of the velocity field. From these visualizations we can appreciate that there is a very high spatial correlation between them. These results suggest that the approximation of neglecting the extra third order terms coming from Chapman-Enskog expansion is a good approximation of the ELBM eddy viscosity.

Let us now analyse the statistics of the SGS energy transfers, comparing them also with the statistics of the real SGS energy transfer, see Eq. (17), measured a priori from fully resolved DNSx2. To obtain the a priori $\overline{{\Pi}}^{\Delta}$ we filter the velocity field with a sharp projector in Fourier space with a cutoff at the maximum wavenumber allowed in the modelled simulations, which corresponds to the dealiasing cutoff ($k_{max}=171$). As known, the presence of a forward energy cascade, as in 3d turbulence, reflects in a skewed PDF of the a priori $\overline{{\Pi}}^{\Delta}$, see gotoh2005statistics ; buzzicotti2018energy ; buzzicotti2018effect , see Fig. 4. Instead, the negative tail describes the presence of intense backscattering events with fluctuations up to two orders of magnitude larger than the standard deviation. The main remarkable difference between the different models considered here is that, as expected, the ELBM and E-LES produce backscatter events, while the Smagorinsky model is positive definite in its energy formulation and it produces a zero tail in the negative region. Let us notice that ELBM mesoscopic model shows qualitatively a better overlap with respect to the DNS data.

In this last section we analyse the inertial range statistics by measuring the longitudinal velocity increments defined as $\delta_{r}u=(\bm{u}(\bm{x}+\bm{r})-\bm{u}(\bm{x}))\cdot\bm{r}/r$. In this way we can quantify the effects produced by the different models at different scales, $r=|\bm{r}|$, hence, we can get an accurate estimation of the quality of the models in capturing the correct intermittent properties of the NSE. We study the scaling properties of the longitudinal structure functions (SF) defined as,

$$S_{p}(r)\equiv\langle[\delta_{r}u]^{p}\rangle$$

(22)

where the angular brackets indicate the ensemble average, that assuming spatio-temporal ergodicity can be evaluated averaging over space and time, $\langle(...)\rangle=\frac{1}{V}\frac{1}{T}\int_{V}\int_{t_{0}}^{t_{0}+T}(...)\,d\bm{x}dt$. In the limit of large Reynolds number, where $r$ can be taken arbitrarily small the structures function follows a powerlaw scaling behavior, $S_{p}(r)\sim r^{\xi_{p}}$ Frisch1995 , where a p-th order scaling exponent that, according to the phenomenological theory of Kolmogorov (K41) Kolmogorov1941 , is $\xi_{p}=p/3$. Nevertheless, both experimental and numerical studies have highlighted as bi-product of intermittency the presence of anomalous exponents in turbulent data, with important deviations from the K41 predicted values Gotoh2002 ; sinhuber2017dissipative ; benzi2010inertial ; biferale2019self . On the other hand, to get an accurate measurement of these exponents is extremely difficult. The reason is twofold, first it is required to have large scaling range (very well resolved simulations) and second it is simultaneously required to have large statistical ensemble. The first question we ask here is connected to the first of the aforementioned problems, namely whether those models are able or not to extend the length of inertial range of scales in our simulations. To answer this in the left column of Fig. 5 we show the 2nd and 4th order structure functions measured from all modelled and fully resolved simulations. In the right column on the same figure we show the local scaling exponents,

$$\xi_{p}(r)=\frac{d\,\log S_{n}(r)}{d\,\log(r)}\,,$$

measured form the structure functions shown on the left side. From these plots it is evident that all modelled simulations present an extension of the inertial range with respect to the same resolution, similar to the inertial range observed in the simulations with a number of grid points two times larger along each spatial direction, DNSx2.

Which means that the models allows to save a factor $8$ in terms of the number of degrees of freedom required in the simulations. Considering that also the time-step needs to be changed accordingly in the higher resolution DNS to resolve the smaller time-scales, it means that the modelled simulations are more than an order of magnitude cheaper than a DNS with the same inertial range. On the other hand, by measuring the local scaling exponents as reported on the right column of Fig. 5, we can see that the inertial range extension produced by the model is not as accurate as the fully resolved DNS. This is a problem if we want the model to correctly describe intermittency. Let us stress that the correction to the K41 prediction at the level of the second and fourth order exponents is very small, hence a model needs to be extremely accurate to correctly capture the intermittent scaling biferale2019self . In both right panels of Fig. 5 we report with solid line the $\xi_{p}$ value of the K41 prediction, and with dotted lines the values measured from DNS as reported in Gotoh2002 ; sinhuber2017dissipative . To highlight intermittency we can look at the ratio between SF at different orders. In particular, any systematic non-linear dependency of $\xi_{p}$ vs $p$, will introduce a scale-dependency in the Kurtosis, defined by the dimensionless ratio among fourth and second order SF,

$$K(r)=\frac{S_{4}(r)}{\left[S_{2}(r)\right]^{2}}.$$

(23)

In Fig. 6 we see that in all simulations, at large scales the increments are Gaussian ($K\sim 3$), while the Kurtosis quickly increases, decreasing the scale. This observation shows non-self similarity in the statistics of all data. It is interesting to observe that at this level the inertial range of scale observed in the DNSx2 simulation are well captured by all closures up to the dissipative scales $r\approx 0.1$ where deviations from the DNSs and models arise.

Going further, we measure the most refined quantity we can observe to quantify intermittency, namely the local scaling exponent in Extended Self-Similarity Benzi1993 ,

$$\zeta(r)=\frac{\xi_{4}}{\xi_{2}}.$$

(24)

A linear K41 behavior would recover in the inertial range a plateau value of $\zeta$ equal to $2$. The correction, accounting for intermittency, measured in both experimental and DNS data gives the plateau for $\zeta$ at the value of $1.86$  Gotoh2002 ; sinhuber2017dissipative ; biferale2019self . As we can see in Fig. 7, all models show deviations from the K41 self-similar prediction meaning that they all capture the non-self-similarity of the turbulent inertial range. However none of them is really accurate enough to extend the length of the plateau displayed, hence to improve the prediction obtainable from the DNS at the same resolution. Indeed, if we compare the modelled data with the fully resolved simulations we can see that former are not showing any flat plateau in the inertial range and we cannot estimate precisely the correction to K41 of the structure function scaling exponents.

It is interesting to point out that the models show a very similar accuracy up to this last analysis. This suggest the backscatter events of energy introduced by the entropic closures are not accurate enough to improve quantitatively the statistics of the Smagorinsky model. This results supports the observation that intermittency in turbulent flows comes as a result of highly non-trivial correlations among all degrees of freedom at different scales buzzicotti2016phase ; buzzicotti2016lagrangian ; lanotte2015turbulence . The observation that for all models we have a very similar inertial range statistics goes in agreement with the common property of these models to have the same scaling, proportional to the strain rate tensor.

In this paper, we performed a quantitative assessment of the ELBM capabilities in the modelling of 3d homogeneous isotropic turbulent flows by comparing the inertial range statistics of ELBM data with the one of high resolution DNS of the NSE. We also compared the quality of ELBM with respect to the hydrodynamical Smagorinsky model, popular in the realm of LES. Furthermore, in this work we have proposed and investigated for the first time, a new hydrodynamical closure for LES simulation inspired from the macroscopic approximation of the ELBM model introduce by Malaspinas2008 . We found that ELBM extends the length of the inertial range with respect to the fully resolved DNS, allowing to reduce the computational cost by an order of magnitude and at the same time preserving the simulation stability. Results showed that, in both the macroscopic and mesoscopic formulations, ELBM is able to reproduce an inertial range with a non self-similar dynamics. ELBM captures the correct deviations from the large scales Gaussian statistics as observed in the fully resolved DNS, with an accuracy comparable to the one produced by the Smagorinsky model. From the measure of the structure functions scaling exponents in ESS, we have highlighted the limitations of these models to get with high accuracy the turbulence corrections to the Kolmogorov scaling. In this context we found that the modelled data are not producing the same inertial range plateau as observed in the fully resolved DNS and experiments. To conclude, we found that ELBM suffers in the modeling of extreme and rare intermittent fluctuations, while on the other hand it is very efficient in modeling the mean properties of 3d turbulence. Which makes ELBM a good candidate for the modeling of 3d turbulent flows in complex geometries.

The authors thank Prof. Luca Biferale for inspiration and many useful discussions. This work was supported also by the European Unions Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie grant [grant number 642069] and by the European Research Council under the ERC grant [grant number 339032]. The authors would like to thanks Prof. Dirk Pleiter as well as the Juelich Supercomputing Center for providing access to the JURON cluster.