# A Data-Driven Dynamic Nonlocal Subgrid-Scale Model for Turbulent Flows

S. Hadi Seyedi  Department of Mechanical Engineering, Michigan State
University, 428 S Shaw Ln, East Lansing, MI 48824, USA  Mohsen Zayernouri
<EMAIL_ADDRESS>Department of Mechanical Engineering, Michigan State
University, 428 S Shaw Ln, East Lansing, MI 48824, USA  Department of
Statistics and Probability, Michigan State University, 619 Red Cedar Road,
Wells Hall, East Lansing, MI 48824, USA

###### Abstract

We developed a novel autonomously dynamic nonlocal turbulence model for the
large and very large eddy simulation (LES, VLES) of the homogeneous isotropic
turbulent flows (HIT). The model is based on a generalized (integer-to-
noninteger) order Laplacian of the filtered velocity field, and a novel
dynamic model has been formulated to avoid the need for tuning the model
constant. Three data-driven approaches were introduced for the determination
of the fractional-order to have a model which is totally free of any tuning
parameter. Our analysis includes both the a priori and the a posteriori tests.
In the former test, using a high-fidelity and well-resolved dataset from
direct numerical simulations (DNS), we computed the correlation coefficients
for the stress components of the subgrid-scale (SGS) stress tensor and the one
we get directly from the DNS results. Moreover, we compared the probability
density function of the ensemble-averaged SGS forces for different filter
sizes. In the latter, we employed our new model along with other conventional
models including static and dynamic Smagorinsky into our pseudo-spectral
solver and tested the final predicted quantities. The results of the newly
developed model exhibit an expressive agreement with the ground-truth DNS
results in all components of the SGS stress and forces. Also, the model
exhibits promising results in the VLES region as well as the LES region, which
could be remarkably important for the cost-efficient nonlocal turbulence
modelings e.g., in meteorological and environmental applications.

††preprint: AIP/123-QED

##  I Introduction

The prohibitively high computational cost of direct numerical simulations
(DNS) of realistic turbulent flows has motivated the community of research in
turbulence to develop coarse-grained techniques including Reynolds-averaged
Navier–Stokes (RANS) and large-eddy simulation (LES) methods to reduce the
intractably large degrees of freedom in DNS studiesPope (2001); Sagaut (2006);
Meneveau and Katz (2000); Germano  _et al._ (1991); Moin  _et al._ (1991);
Vasilyev, Lund, and Moin (1998). Using Reynolds averaging approach in the
Navier–Stokes (N-S) equations provide temporally averaged quantities, however,
in LES approaches one employs a subgrid-scale (SGS) model, which represents
the effects of the finer scalesPope (2004, 2001).

Turbulence experimental and DNS features have confirmed that the turbulence is
intrinsically nonlocal, and its statistics are non-Gaussian, which means
velocity increments have sharp peaks, heavy-skirts and also skewedWilczek,
Vlaykov, and Lalescu (2017); Akhavan-Safaei, Seyedi, and Zayernouri (2020).
However, most of the turbulence models have been built based on the
Boussinesq’s turbulent viscosity concept, in which one assumes turbulent
stress tensor is proportional to the local mean velocity gradient at any
point, and the proportionality coefficient is set to the turbulent viscosity.
Prandtl in 1942 aimed to disregard this local constraint by introducing the
extended mixing length concept for the first time. The new model was a
migration from locality to nonlocality, however, the model and its
implementation were not remarkably successful since the scale of nonlocality
was comparable with the differential length scale. Afterward, Prandtl
parametrized the model in a way that the mixing length was taken to be bigger
than the differential length. This strategy was the same as adding a weak
nonlocal concept to the model, hence called weak nonlocal. Bradshaw Bradshaw
(1973) in 1973 showed that Boussinesq’s hypothesis fails over curved surfaces
and noted that form of the stress-strain relations is responsible for this
failure. It should be mentioned that there were some important developments
mostly based on polynomial series compared to the Boussinesq type modeling
including the works done by Spencer and Rivlin Spencer and Rivlin (1958,
1959), Lumley Lumley (1970) and PopePope (1975); however, they all lacked the
accuracy that a “true” physical modeling should provide especially for the
second-order and higher tensor series development.

Using generalized-order derivatives is a relatively convenient approach to
bring in nonlocality concept from mathematical point of view. The generalized-
order operators represent the underlying heavy-tailed stochastic processes at
the continuum level, which can be properly utilized in incorporating the long-
range interactions in various mathematical models including but are not
limited to beam vibration analysis Suzuki  _et al._ (2021a), anomalous
rheology modeling Suzuki  _et al._ (2021b); Naghibolhosseini and Long (2018);
Naghibolhosseini (2015), damage modeling considering memory effects Suzuki
_et al._ (2021c) and visco-elasto-plastic models Suzuki  _et al._ (2016).
Moreover, harnessing the generalized-order models capabilities can be obtained
properly using the highly accurate numerical schemes for integer and
fractional-order PDEs D’Elia  _et al._ (2020); Kharazmi and Zayernouri (2019);
Samiee, Zayernouri, and Meerschaert (2019); Zhou  _et al._ (2020); Seyedi,
Nemati Saray, and Ramazani (2019); Lischke, Zayernouri, and Zhang (2019);
Seyedi (2018); Zayernouri and Karniadakis (2014); Kharazmi, Zayernouri, and
Karniadakis (2017); Zayernouri and Karniadakis (2013); Kharazmi and Zayernouri
(2018), which is also an active research topic.

In a pioneer work by Hinze et al. Hinze, Sonnenberg, and Builtjes (1974) in
1974, the authors described the memory effect in a turbulent boundary layer
flow. They used the experimental data, produced downstream of a hemispherical
cap attached to the lower wall of channel geometry and illustrated that, when
one computes eddy-viscosity using Boussinesq’s theory in the lateral gradient
of the mean flow, there is a significant non-uniform distribution that also
exists in the outer region of the boundary layer. Interestingly, a nonlocal
expression for the gradient of the transported field was proposed in a novel
approach by Kraichnan in the same year for the scalar quantity transport
Kraichnan (1964). Afterward, fractional-order models based on the RANS
approach were offered in Chen (2006); Epps and Cushman-Roisin (2018); Egolf
and Hutter (2017); Hamba (1995, 2005). One of the main contributions in the
development of nonlocal RANS models is by Egolf and Hutter Egolf and Hutter
(2017); Egolf and Kutter (2020). They started from Lévy flight statistics and
generalized the zero-equation local Reynolds shear stress expression to a
nonlocal and fractional type. The method is based on Kraichnanian convolution-
integral approach and utilizing different weighting functions. Using the
mentioned weighting functions, one can make a bridge between the first-order
gradient of the common eddy diffusivity models and the mean velocity
difference term. Their proposed model is called the Difference-Quotient
Turbulence Model (DQTM) and is based on the four distinct steps that can be
followed to change a local operator to a nonlocal one. In reality, the
proposed model is a more general version of Prandtl’s zero-equation mixing
length and shear-layer turbulence models.

There is also emerging attention for the nonlocal LES models. Samiee et al.
proposed a new model for the HIT flows based on the fractional Laplacian by
employing Lévy stable and tempered Lévy stable distributions in kinetic level
Samiee, Akhavan-Safaei, and Zayernouri (2020, 2021). They showed that the new
nonlocal models can recover the non-Gaussian statistics of subgrid-scale
stress motions. Laval et al. Laval, Dubrulle, and Nazarenko (2001) analyzed
the effects of the local and nonlocal interactions on the intermittency
corrections in the scaling properties. They observed that nonlocal
interactions are responsible for the creation of the intense vortices and on
the other hand, local interactions are trying to dissipate them. Akhavan-
Safaei et al. Akhavan-Safaei, Samiee, and Zayernouri (2021) proposed a
fractional LES approach for the subgrid-scale modelings of the scalar
turbulence. They utilized the two-point statistics for defining the optimal
fractional-order of the new nonlocal model, and by using a priori assessment
they showed that there is proper agreement between the probability
distribution function (PDF) of the SGS dissipation and the one that comes from
the filtered DNS data. Harmonious with this study, Akhavan-Safaei and
Zayernouri developed a corresponding nonlocal spectral transfer model and a
new scaling law for scalar turbulence in Akhavan-Safaei and Zayernouri (2021).
Their new analysis additionally reconciled the close similarities between this
work and their earlier development in Akhavan-Safaei, Samiee, and Zayernouri
(2021) when the filter scale approaches the dissipative scales of turbulent
transport. There are also other related studies that one can consult with,
including preliminary fractional modeling in wall-bounded turbulent flows
Keith, Khristenko, and Wohlmuth (2021), a priori survey of nonlocal eddy
viscosity-based model for the isotropic and anisotropic (channel flow
canonical test cases) turbulent flows Di Leoni  _et al._ (2021), hybrid
nonlocal model in the case of magnetically confined plasma Milovanov and
Rasmussen (2014), generalization of a deconvolution model with fractional
regularization for the rotational N-S equations Ali (2014), and fractional
Laplacian closure and its connection to Richardson pair dispersion Gunzburger,
Jiang, and Xu (2018). Going even beyond the scope of research in turbulence, a
new comprehensive survey on the nonlocal models for several crucial
applications, including anomalous subsurface dynamics, turbulence modeling,
and extraordinary materials, was recently performed by Suzuki et al. Suzuki
_et al._ (2021d).

Considering the nonlocal models in the literature, there are some important
imperfections including the sensitivity to the model constant and fractional-
order parameter, relatively low correlation coefficients
$\rho(\tau^{\Delta},\tau^{\Delta,Model})$, and no back-scatter prediction of
kinetic energy from small scales to large scales. However, the conventional
and frequently utilized local LES turbulent models are being improved over
time to be free from mentioned deficits. One of the methods in the improvement
process is using the dynamic procedure for the determination of the model
constant Germano  _et al._ (1991).

To fill the gap in the literature and provide an applicable and relatively
easy to implement nonlocal LES model, we have developed a new dynamic nonlocal
model that accounts for all the aforementioned downsides. In the new dynamic
fractional subgrid-scale model (D-FSGS), both nonlocality and dynamic features
have been leveraged together for the first time. This match between two
important features provides a unique and higher performance than the dynamic
local or static nonlocal models. Interestingly, the analysis showed that in
the new model, we have remarkably less sensitivity to the fractional-order,
which is needed to be specified in nonlocal models. This relative freedom is
obtained thanks to the novel coupling between the dynamic procedure and the
nonlocal nature of the base model. In the following, we derived and
implemented the D-FSGS model in both a priori and a posteriori stringent tests
and compared the results with the conventional local models including
Smagorisnky (SMG) and dynamic Smagorinsky (D-SMG) models.

This study is organized as follows: in section II we talk about the governing
equations and development of the new model. Section III starts with
introducing three main distinct approaches for the determination of the
fractional-order, and then we do a comprehensive a priori assessment along
with a comparative study with conventional models. In section IV, we study the
numerical stability of the proposed model and do different tests including
two-point diagnostics to have a complete overview of the model performances in
a a posteriori sense. Finally, in section V we recap the findings with a
conclusion.

##  II Model Development

Implementation of a low-pass filter on the N-S equation forms a closure term
on the right-hand side of the momentum equation, which needs to be modeled,
representing the (unknown) SGS dynamics. The filtered incompressible N-S
equations can be written as

$\displaystyle\frac{\partial\bar{u}}{\partial t}+\bar{u}\cdot\nabla\bar{u}$
$\displaystyle=-\frac{1}{\rho}\ \nabla\bar{p}+\nu\
\nabla^{2}\bar{u}-\nabla\cdot\tau,$ (1) $\displaystyle\nabla\cdot\bar{u}$
$\displaystyle=0.$

In this equation, $\bar{u}$ is the filtered velocity vector, $\rho$ denotes
the constant density, $\bar{p}$ represents the filtered pressure, and $\nu$
shows the kinematic viscosity. The effects of the small scales arise in the
so-called SGS stress tensor

$\displaystyle\tau_{ij}=\overline{u_{i}u_{j}}-\bar{u_{i}}\bar{u_{j}},$ (2)

forming the closure term to be modeled. Several models have been proposed
during the last decades to close the filtered N-S equation in both functional
and structural LES models Pope (2001); Meneveau and Katz (2000); Sagaut
(2006). In the classical Smagorinsky model (SMG)Smagorinsky (1963), the
deviatoric part of the stress tensor is written as

$\displaystyle\tau_{ij}^{SMG}=-2\ \nu_{t}\ \bar{S_{ij}},$ (3)

where $\nu_{t}$ denotes the eddy viscosity, and
$\bar{S_{ij}}=\frac{1}{2}(\frac{\partial\bar{u_{i}}}{\partial
x_{j}}+\frac{\partial\bar{u_{j}}}{\partial x_{i}})$ represents the filtered
(resolved) strain rate tensor. In this model, $\nu_{t}$ is constructed based
on the Prandtl’s Mixing Length hypothesis,

$\displaystyle\nu_{t}=C_{s}\ \mathcal{L}^{2}\ |\bar{S}|,$ (4)

where $\mathcal{L}$ shows the effective grid scale, and
$|\bar{S}|=\sqrt{2\bar{S_{ij}}\bar{S_{ij}}}$ exhibits the magnitude of the
resolved scale strain rate tensor Meneveau and Katz (2000). One of the main
drawback of this model is its flow-dependent feature, which means that the
model can not correctly predict final quantities with a single universal
constant in different scenarios such as shear flows, wall-bounded flows, or
transitional flows. To overcome this challenge in the Smagorisnky model,
Germano et al. Germano  _et al._ (1991) proposed a novel procedure for
evaluation of the model coefficient, which is called the dynamic Smagorinsky
model (D-SMG). The suggested procedure was a breakthrough in turbulence
modeling, and several researchers utilized the same concept afterward Lund and
Kaltenbach (1995); Najjar and Tafti (1996); Piomelli and Liu (1995); Ghosal
and Rogers (1997). The dynamic procedure was designed based on the classic
idea that one can extract useful information from the smallest resolved scales
for modeling the subgird-scales; however, the way it was applied was indeed
novel. It calculates the eddy-viscosity coefficient locally for each LES grid
point as the calculation progresses, and there is no need for any predefined
inputs in the model. This model was constructed based on the scale-invariance
hypothesis and calculates the model constant using the information from the
resolved section.

###  II.0.1 Fractional SGS Model Samiee, Akhavan-Safaei, and Zayernouri
(2020):

We have recently developed a nonlocal SGS model that we present here briefly
for making this work self-contained. Starting from the Boltzmann kinetic level
description of the flow, we employ the BGK model for the collision of the
particles and have

$\displaystyle\frac{\partial f}{\partial t}+u\cdot\nabla
f=-\frac{f-f^{eq}}{\mathcal{T}},$ (5)

where $f=f(x,u,t)$ is called the single-particle probability distribution
function and shows the particles density in the phase space $(x,u)$ at time
$t$. Moreover, $f^{eq}$ represents the local equilibrium distribution function
and is defined based on the Maxwell distribution.

$\displaystyle f^{eq}(\mathcal{K})=\frac{\rho}{U^{3}}F(\mathcal{K}),$ (6)

where $F(\mathcal{K})=e^{-\mathcal{K}/2}$,
$\mathcal{K}=\frac{|{u}-V|^{2}}{U^{2}}$ and ${U}$ corresponds for the thermal
agitation speed. The left-hand side of Eq. 5 correspond to the streaming of
the non-reacting particles and right-hand side shows the collision operator
with a relaxation time $\mathcal{T}$. One can solve Eq. 5 analytically by
method of characteristic to find distribution in terms of equilibrium state
Chen, Orszag, and Staroselsky (2010)

$\displaystyle f(t,x,u)=\int_{0}^{\infty}e^{-s}\,{f^{eq}(x-u\tau s,u,t-\tau
s)}\,ds.$ (7)

Moments of the $f$ would provide the macroscopic flow variables. Therefore,
one can write $\rho=\int f(t,x,u)du$, and $\rho V(t,x)=\int uf(t,x,u)du$ to
compute density and fluid velocity. Incorporating filtering procedure into the
Eq. 5 would provide filtered Boltzmann equation as

$\displaystyle\frac{\partial\bar{f}}{\partial
t}+u\cdot\nabla\bar{f}=-\frac{\bar{f}-\overline{f^{eq}(\Delta)}}{\mathcal{T}}.$
(8)

However, the collision term is highly-nonlinear and filtering kernel cannot
commute Girimaji (2007). Hence, a closure problem would be built by defining
$\overline{\mathcal{K}}=\frac{|{u}-\overline{V}|^{2}}{U^{2}}$ since

$\displaystyle\overline{f^{eq}(\mathcal{K})}\neq
f^{eq}(\overline{\mathcal{K}}).$ (9)

One of the approaches for handling this problem is using a power-law
distribution for modeling
$\overline{f^{eq}(\mathcal{K})}-f^{eq}(\overline{\mathcal{K}})$, i.e.,

$\displaystyle\overline{f^{eq}(\mathcal{K})}-f^{eq}(\overline{\mathcal{K}})\simeq
f^{Model}(\bar{\mathcal{K}})=\mathfrak{D}_{\beta}\,f^{\beta}(\bar{\mathcal{K}}),$
(10)

in which
$f^{\beta}(\bar{\mathcal{K}})=\frac{\rho}{U^{3}}F^{\beta}(\mathcal{K})$, where
$F^{\beta}(\mathcal{K})$ represents an isotropic Lévy $\beta$-stable
distribution. Also, $\mathfrak{D}_{\beta}$ is a constant which is going to be
addressed in the present dynamic model. In the fractional Laplacian model, the
SGS forces are defined as

$\displaystyle(\nabla.\tau)_{i}=\nu_{\alpha}(-\Delta^{\alpha})\bar{u}_{i},\quad\alpha\in(0,1],$
(11)

where the filtered velocity
$\bar{u}_{i}:\mathbb{R}^{d}\times(0,T]\rightarrow\mathbb{R}$, where $d=3$ is
the dimension of physical domain and $T$ represents the simulation time, in
addition, $(-\Delta^{\alpha})(\cdot)$ denotes the space-fractional Laplacian
of order $2\alpha\in(1,2]$, which can be defined as a singular integral
operator given by

$\displaystyle(-\Delta)^{\alpha}\bar{u}_{i}=c_{d,\alpha}\int\limits_{\mathbb{R}^{d}}{\frac{\bar{u}_{i}(\bm{x})-\bar{u}_{i}(\bm{y})}{|\bm{x}-\bm{y}|^{d+2\alpha}}\,d\bm{y}},$
(12)

where
$c_{d,\alpha}=\frac{4^{\alpha}\Gamma(d/2+\alpha)}{\pi^{d/2}|\Gamma(-\alpha)|}$.
Moreover employing the periodic boundary conditions, the corresponding Fourier
transform of the fractional Laplacian in (12) is given by (see e.g., Lischke
_et al._ (2020))

$\displaystyle\mathcal{F}\Big{[}(-\Delta^{\alpha})\bar{u}_{i}(\bm{x})\Big{]}=|\bm{k}|^{2\alpha}\mathcal{F}[\bar{u}_{i}](\bm{k}),$
(13)

in which $\mathcal{F}$ and $\bm{k}$ are the Fourier transform and Fourier
numbers, respectively. Evidently, the integer-order Laplacian operator is
recovered simply by putting $\alpha=1$. The Fourier transform 13 provides a
rather convenient way of handling the fractional operators Fourier space and
in our Fourier spectral method for simulating the problem. The corresponding
“eddy-viscosity"-like model coefficient in (11) is then obtained as

$\displaystyle\nu_{\alpha}=\mathcal{C}\,F(\alpha),$ (14)

where $\mathcal{C}$ represents an up-scaling model input being proportional to
$U^{{}^{2\alpha}}\mathcal{T}^{2\alpha-1}$ from the kinetics theory’s
perspective, yet to be dynamically computed in the subsequent (continuum-
level) simulations, moreover, $F$ denotes a deterministic univariate function
of the fractional-order $\alpha$, explicitly given by

$\displaystyle F(\alpha)=\frac{2^{2\alpha}\
\Gamma(\frac{2\alpha+3}{2})\,\,\Gamma(2\alpha+1)}{\pi^{\frac{3}{2}}|\Gamma(-\alpha)|},$
(15)

rendering (11) as
$(\nabla.\tau)_{i}=\mathcal{C}\,F(\alpha)\,(-\Delta^{\alpha})\bar{u}_{i},\,\alpha\in(0,1]$.
In what follows and in this generalized order context, we develop a new
dynamic procedure to automatically compute $\mathcal{C}$ from data on-the-fly.

###  II.0.2 Derivation of the Nonlocal Dynamic Model Procedure:

We write $(\nabla.\tau)_{i}=C_{1}F(\alpha)(-\Delta^{\alpha})\bar{u}_{i}$ for
the sake of simplicity. Implementing the second filtering process (test-level
filter), gives the divergence of the SGS stresses at the test-level filter
(subtest-scale stress) as

$\displaystyle(\nabla.T)_{i}=C_{2}F(\alpha)(-\Delta^{\alpha})\widehat{\bar{u}_{i}}.$
(16)

In the above equation, $\widehat{(\cdot)}$ indicates the test-level filtering,
which is commonly chosen as twice the grid-level filtering. Now, the Germano
identity Germano  _et al._ (1991), which relates the stresses at grid-level
$(\tau_{ij})$ and test-level $(T_{ij})$, is employed to make a bridge between
the resolved scales and the subgrid-scales,

$\displaystyle
G_{ij}=\widehat{\bar{u}_{i}\bar{u}_{j}}-\widehat{\bar{u}_{i}}\widehat{\bar{u}_{j}}=T_{ij}-\widehat{\tau_{ij}},$
(17)

which is a known quantity, and represents the resolved turbulent stress. In
the divergence form, we have

$\displaystyle\nabla.\big{(}\widehat{\bar{u}_{i}\bar{u}_{j}}-\widehat{\bar{u}_{i}}\widehat{\bar{u}_{j}}\big{)}=(\nabla.G)_{i}.$
(18)

Therefore, we construct the Germano identity in the divergence form based on
the previously introduced fractional Laplacian model. The Germano identity
extrapolates and parametrize the model constant for the subgrid-scale part
using the information in the smallest resolved scales,

$\displaystyle(\nabla.T)_{i}-(\nabla.\widehat{\tau})_{i}=C\bigg{(}F(\alpha)(-\Delta^{\alpha})\widehat{\bar{V}_{i}}-\ensurestackMath{\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("F(\alpha)(-\Delta^{\alpha})\bar{V}_{i}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{F(\alpha)(-\Delta^{\alpha})\bar{V}_{i}}{\tmpbox}}\bigg{)}$

Here, we assume the model constant $C$ being scale invariant around the grid
and test filter size, ($C=C_{1}=C_{2}$), and $C$ is in general a space-time
dependent tensor, yet, it has been assumed to be spatially uniform. One needs
to solve Eq. II.0.2 for $C$, which forms a tensorial equation. There are two
common methodologies to get a scalar model constant. First, one may contract
with $\bar{S_{ij}}$ as mentioned originally in Germano  _et al._ (1991).
Second approach is proposed by Lilly Lilly (1992), which is more preferable
and commonly used. This approach is based on contracting with the tensor,
which has been multiplied to the unknown coefficient in the right hand side of
Eq. II.0.2.

(a)

(b)

Figure 1: Difference between commuting of the low pass filter on the (a)
integer order derivatives, and (b) fractional-order derivative ($\alpha=0.6$).

One important point that should be mentioned in this section is that in the
conventional local LES modeling, the main assumption is that the operation of
filtering and integer-order differentiation commutes, then we get a set of
filtered equations. Commuting means commutation between the filtered quantity
and the spacial derivatives. However, the filtering procedure does not commute
with the fractional-order operators in general. One simple example regarding
this important matter is depicted in Figure 1 for $\alpha=0.6$, comparing the
derivatives of the filtered velocity at test-level and filtered derivative for
the integer and fractional-order cases.

We have used Lilly’s approach for contracting purposes based on the Least
Square Method that was described by Lilly in 1992 Lilly (1992). Therefore, one
can write Eq. II.0.2 as

$\displaystyle\big{(}(\nabla.G)_{i}-CN_{i})^{2}=e,$ (19)

in which $e$ is the squared error and the nonlocal term $N_{i}$ is defined as

$\displaystyle
N_{i}=\big{(}F(\alpha)(-\Delta^{\alpha})\widehat{\bar{V}_{i}}-\ensurestackMath{\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("F(\alpha)(-\Delta^{\alpha})\bar{V}_{i}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{F(\alpha)(-\Delta^{\alpha})\bar{V}_{i}}{\tmpbox}}\big{)}$

Finding the unknown in minimized form is conveniently possible by putting the
derivative equals to zero, considering and testing
$\frac{\partial^{2}e}{\partial C^{2}}>0$ for the error minimization scope.
Finally, the scalar model coefficient $C$ can be computed dynamically as

$\displaystyle C=\frac{\langle(\nabla.G)_{i}N_{i}\rangle}{\langle
N_{i}N_{i}\rangle}.$ (20)

Numerical instability may be occurred due to the negative eddy-viscosity in
prolonged periods of time. As a remedy, one can perform an averaging over the
directions of statistical homogeneity as suggested by Germano et al. Germano
_et al._ (1991). Figure 2 illustrates the variations of model constants in
D-SMG and D-FSGS models in the imaginary center-line of a periodic domain in
the first direction to have a comparison between the model constant variations
in the context of an example.

Figure 2: Comparing the model coefficients in dynamic Smagorinsky model
(D-SMG) and dynamic FSGS model (D-FSGS) in the middle imaginary line in a
periodic domain using ${\mathcal{L}_{\delta}}=4$.

##  III A Priori Analysis

We assess the performance of the D-FSGS and compare the results with the
results of conventional LES models including SMG and D-SMG. Also, we discuss
about the performance of the D-FSGS model using different filter sizes
$\mathcal{L}_{\delta}$ (covering both the LES and VLES regions). Before
starting the a priori tests, we need to study the effects of the fractional-
order parameter ($\alpha$) and its proper selection for different scenarios.

We utilize ten independent snapshots of forced HIT data distributed over
enough eddy turnover times on a triply periodic domain, $\Omega=[0,2\pi]^{3}$,
for each targeted Taylor Reynolds number. The covered Taylor Reynolds numbers
are $170$, $240$, and $355$, and the corresponding grid resolution for each
case is $320^{3}$, $520^{3}$, $1024^{3}$ , respectively. The turbulence
dataset comes from our open-access pseudo-spectral parallel code, which has
been elaborated in Akhavan-Safaei and Zayernouri (2020).

###  III.0.1 The Inference of Optimum Laplacian-Order $(\alpha^{opt})$

The D-FSGS model depends on the fractional-order, $\alpha$. We provide three
different data-driven approaches for the calculation of this parameter.

#### First Approach:

we perform a precise comparison between the obtained correlation coefficients,
considering the effects of all stress components in the SGS tensor and the
ground-truth stresses that come from the DNS results. Subsequently, we find
the $\alpha^{opt}$, where the maximum of $\rho_{ij}=\langle
Avg(\rho[\tau^{DNS}_{ij},\tau^{Model}_{ij}])\rangle$ is achieved, in which,
$\langle\cdot\rangle$ denotes the ensemble-averaged on different snapshots of
data, and for each snapshot of data, we consider the mean of all corresponding
components of the SGS tensor.

By sweeping the fractional-order values from zero to one and step size equals
to 0.01, we determined the $\alpha^{opt}$ for different characteristic filter
sizes in three different Reynolds numbers. In Figure 3, we have shown the
$\alpha^{opt}$ values for different filter sizes for three different Taylor
Reynolds numbers. Also $\mathcal{L}_{\delta}=\frac{\mathcal{L}}{2\delta_{X}}$,
where $\delta_{X}=\frac{2\pi}{N}$ shows the computational grid size associated
with $N=320,520,1024$ for three different Reynolds numbers. These results are
obtained based on the ensemble-averaged quantities of ten snapshots of data
for each Reynolds number to ensure meeting the necessary conditions regarding
performing LES studies Meneveau (1994). There is an interestingly noticeable
trend between the filter size and the obtained $\alpha^{opt}$ from the
ensemble-averaged relation for each Reynolds number (see Figure 3).

#### Second Approach:

we employ a logarithmic regression to come up with a correlation between the
$\alpha^{opt}$ and ${\mathcal{L}_{\delta}}$ for each Reynolds number. The
results are showing that $\alpha^{opt}$ values obey a power-law form, which
reads

$\displaystyle\ \alpha^{opt}=a\ {\mathcal{L}_{\delta}}^{b}.$ (21)

Table 1: Proper parameters for being used in second approach (Eq. 21) for determination of optimum Laplacian-order. $Re_{\lambda}$ |  | $a$ |  | $b$
---|---|---|---|---
170 |  | 1.53 |  | -0.62
240 |  | 1.39 |  | -0.49
355 |  | 1.08 |  | -0.30

Increasing the filter size would incorporate more nonlocality, and that is the
reason for the increase in the error bars in the VLES section; however, all
cases converge statistically to the realized power-law relations. Moreover,
there is a direct relationship between the Taylor Reynolds number and the
optimum fractional-order in each filter size. Using the above relation for
finding the $\alpha^{opt}$ would be the second method that can be applied
conveniently and confidently since the coefficient of determination, $R^{2}$,
for these equations are above $0.98$, and indicates that the relation
correlates very well with the measured data. The $a$, $b$ values in Eq. 21 are
being determined for different taylor Reynolds numbers. Table 1 shows the
proper parameters for being used in second approach of optimum Laplacian-order
determination. The values in this table obtained based on the ensemble-
averaged results of ten data snapshots for each $Re_{\lambda}$. Moreover,
interpolation methods can be conveniently utilized using this table to find
the proper value of $\alpha^{opt}$ for other Reynolds numbers.

(a)

(b)

(c)

(d)

Figure 3: Obtained optimum fractional-orders at different filter sizes based
on the first approach, (a) effect of $Re_{\lambda}$, amount of modeled
turbulent kinetic energy when (b) $Re_{\lambda}=170$, (c) $Re_{\lambda}=240$,
and (d) $Re_{\lambda}=355$.

#### Third Approach:

we propose a third approach for obtaining the proper fractional-order for a
quite wide range of (fine-to-very) LES region yet still in the a priori sense,
shown in Figure 3. The right vertical axes (secondary axis) in Figures 3, 3, 3
show the volumetric averaged values of the modeled kinetic energy in percent
($\phi$) at that specific filter size. To calculate this value, $\phi$, we
find the ratio between the resolved ($E_{r}$) and filtered ($E_{f}$) kinetic
energy as

$\displaystyle E_{f}=\frac{1}{2}\ \overline{u_{i}}\ \overline{u_{i}},$ (26)
$\displaystyle E_{r}=\frac{1}{2}\ (\overline{u_{i}u_{i}}-\overline{u_{i}}\
\overline{u_{i}}).$

Noting that almost all of the LES studies are designed to model about less
than ten percent of the total kinetic energy, one can conclude based on the
obtained plots for all $Re_{\lambda}$ numbers, taking $\alpha=0.5$ can be an
average corresponding factor for the LES region. We will show that even with
this rough estimation, the proposed model provides better correlations in all
of the stress components than the other conventional models. Therefore, using
the second or third approach, we will have a dynamic nonlocal model, which is
totally free of any tuning parameter including fractional-order.

###  III.0.2 Statistical Performance Assessment

We use ten three-dimensional snapshots of velocity fields, $u_{i}(x),i=1,2,3$,
corresponding to $Re_{\lambda}=240$ and perform a priori testing to show the
model performance. The DNS datasets with $Re_{\lambda}=240$ and $N=520$ will
be utilized for the rest of the paper.

(a)

(b)

Figure 4: A priori assessment in the context of the ensemble-averaged
correlation coefficient comparisons in different models for (a)
${\mathcal{L}_{\delta}}$ = 4, and (b) ${\mathcal{L}_{\delta}}$ = 8 .

We have depicted the ensemble-averaged correlation coefficients for all of the
stress components in Figure 4 for the filter size equals to $4$ and $8$. We
see a remarkable better performance in both filter sizes and all directions
for the new proposed model, D-FSGS. Interestingly, the coupling between the
new dynamic procedure and the nonlocal nature of the model helps significantly
to compensate for the directions in which we have low correlations. We see
that D-FSGS outperforms static nonlocal (FSGS) and dynamic local (D-SMG)
models in both filter sizes and all directions. Its well-known that, dynamic
local models do not necessarily increase the correlation coefficient
remarkably. In contrast to the D-SMG model, we see a significant raise in the
correlation coefficient of the D-FSGS versus FSGS. Therefore, the new dynamic
procedure not only adds the back-scatter prediction and being free from tuning
constant capabilities, it also improve the performance as well. In another
word, the dynamic procedure in the context of the nonlocal modeling seems to
be more advantageous and fruitful than the local modeling.

In Figure 5 we have plotted the $\prod=-\langle\tau_{ij}\bar{S_{ij}}\rangle$,
the SGS dissipation of kinetic energy for the same filter sizes mentioned
previously. The negative values of $\prod$ mean we have back-scatter, where
the kinetic energy is transferred from subgrid to resolved scales. As the SGS
dissipation of the filtered DNS results (ground-truth) show, a significant
amount of flow shows back-scatter phenomena. We see the new model, D-FSGS, and
D-SMG models are the only ones that truly predict the back-scattering in
energy flow. Moreover, comparing the PDFs of dissipation show that the D-SMG
model over-predicts the values for the tails of the distributions, which can
be a potential factor in facing unstable conditions in simulations with this
model necessitating the averaging operation in a posteriori testsMeneveau and
Katz (2000). However, the new model provides closer results to the DNS ones.
It should also be mentioned that a comparison of the percentage of grid points
with back-scatter in filtered-DNS, D-SMG, and D-FSGS revealed that there is a
better agreement for the D-FSGS model in the prediction of this quantity. As
has been reported by Piomelli et al Piomelli  _et al._ (1991), this percentage
is a function of the filter type that one is using. In this study, three-
dimensional box filtering for both the grid-level and test-level was utilized.
The percentage of grid points with back-scattering in DNS was $26\%$, in D-SMG
was $20\%$ and in D-FSGS was $24\%$ for ${\mathcal{L}_{\delta}}=4$. For bigger
filter size, ${\mathcal{L}_{\delta}}=8$, these numbers were $30\%$, $18\%$ and
$25\%$ for DNS, D-SMG and D-FSGS, respectively, which quantitatively indicates
the better prediction of back-scatter for the D-FSGS model.

(a)

(b)

Figure 5: PDF of ensemble-averaged SGS dissipation of kinetic energy using
different models for (a) ${\mathcal{L}_{\delta}}$ = 4, and (b)
${\mathcal{L}_{\delta}}$ = 8.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: PDF of ensemble-averaged SGS forces using different models at
${\mathcal{L}_{\delta}}$ =4 (first row), and ${\mathcal{L}_{\delta}}$=8
(second row) .

Figure 6 shows the PDF of the ensemble-averaged SGS forces in three directions
using ${\mathcal{L}_{\delta}}=4$ and ${\mathcal{L}_{\delta}}=8$. As it is
clear in all three directions, there is superior math between the results of
the D-FSGS model and DNS ones in comparison to the results of the Smagorinsky-
based, and FSGS model considering tail and peaks of the SGS forces.

Figure 7: Model performance comparison using ensemble-averaged correlation at
different filter sizes using three proposed approaches for the determination
of the fractional-order.

To employ the second and third approaches in the determination of the
fractional-order, we do a test using $\alpha=\alpha^{opt}$ based on the
precise determination method (first approach), then we apply the obtained
power-law relations (second approach) to get $\alpha^{opt}$ and finally
$\alpha=0.5$ (third approach) to show its effect on the performance of the
D-FSGS model. As has been depicted in Fig. 7, even using the rough-estimation
approaches for $\alpha$, the D-FSGS model provides a higher correlation for
all the filter sizes. Moreover, this plot again verifies the existence of
power-law relation between the $\alpha^{opt}$ and filter size. One of the main
disadvantages of nonlocal models in general and nonlocal turbulence models, in
particular, is their sensitivity to the fractional-order. Surprisingly, the
dynamic nonlocal model exhibits a less sensitivity to this parameter around
the value $0.5$, since its effect is compensated by the automatic tuning
parameter, thanks to the designed dynamic procedure. Moreover, the whole
closure term in the LES framework indeed originates from filtering the first-
order nonlinear convective term in the N-S equations. Here, the corresponding
$\alpha^{opt}=0.5$ in fact validates that the overall scaled up and nonlocal
behavior of the closure term in our dynamic model emerges of total
$2\alpha^{opt}=1^{st}$ order.

##  IV A Posteriori Analysis

To assess the practical ability of the proposed dynamic nonlocal model to
capture unsteady large-scale coherent structures and verifies numerical
stability, we carry out the corresponding a posteriori study here. In this
section, we test the performance of the new proposed model (D-FSGS) against
the static and dynamic Smagorinsky models (SMG, D-SMG) as well as the ground-
truth filtered DNS results. For this purpose, pseudo-spectral N-S solver,
which has been discussed in Akhavan-Safaei and Zayernouri (2020), was utilized
for a triply periodic domain as $\Omega=[0,2\pi]^{3}$. The resolution for this
dataset is $520^{3}$ and large-eddy turnover times is
$\tau_{\mathcal{L}}\simeq 2.7$. First, we set up a decaying HIT case in which
the initial condition is based on the statistically stationary data sets.
Thereafter, DNS and filtered-DNS results were gathered for the sake of
comparisons. Second, each LES model including D-FSGS, SMG, D-SMG, and also
unresolved numerical simulation (UNS), which is basically a DNS solver using
the LES grid, implemented using the same pseudo-spectral solver, initiated
with the filtered DNS initial condition for that specific filter size. In the
mean time, the time step of the problems is chosen to have the CFL less than
unity to ensure a stable time-integration.

(a)

(b)

Figure 8: Prediction of the kinetic energy spectra in different turbulence
models at (a) $\frac{t}{\tau_{\mathcal{L}}}\simeq 2$, and (b)
$\frac{t}{\tau_{\mathcal{L}}}\simeq 4$ .

We compared the kinetic energy spectrum obtained from each model with the
ground-truth ones. In this section we are just showing the results of using
${\mathcal{L}_{\delta}}=4$ since in bigger filter widths, there is not enough
resolution especially in the small scale sections. Fig. 8 illustrates energy
spectrum in two relatively high integration times
$\frac{t}{\tau_{\mathcal{L}}}\simeq 2,\ 4$ after initiation of decay. Higher
integration times were chosen to address the numerical stability of models,
which can be a concern especially in the dynamic models. It is well noted that
Smagorinsky model is one of the most dissipative models and this matter is
evident in both plots. Moreover, the UNS case as one may expect should have an
accumulation of noise due to the zero turbulent viscosity. Therefore, we can
expect that the D-FSGS model should be ideally placed between highly
dissipative SMG and UNS energy curves, which act as lower and higher bonds,
respectively.

We continued the implementation of the a posteriori analysis by comparing the
performance of the models in decaying homogeneous isotropic turbulence with
$Re_{\lambda}=240$ using LES filter width ${\mathcal{L}_{\delta}}=4,\ 20$ as
representative of LES and VLES regions, respectively. Since the DNS and LES
solvers have not been initiated with the same initial conditions, we are
comparing the results of filtered-DNS and LES models.

(a)

(b)

Figure 9: Decay of the resolved turbulent kinetic energy in different
turbulence models using (a) ${\mathcal{L}_{\delta}}=4$ as an LES case, and (b)
${\mathcal{L}_{\delta}}=20$ as a VLES case.

Fig. 9 shows the decay of resolved kinetic energy with time. As the plots
show, the D-FSGS model provides better agreement with the filtered-DNS results
in both LES and VLES cases. Interestingly, this difference is more remarkable
in using bigger filter sizes and can be due to the fact that enlargement of
the filter size incorporates more details to the solver and local models are
not successful enough to correctly handling them due to the nature of the
model. On the other hand, the D-FSGS model shows better performance by taking
into account the nonlocal features and fractional derivations.

(a)

(b)

Figure 10: Two-point diagnostics, velocity structure functions in different
models and comparison to the filtered-DNS using ${\mathcal{L}_{\delta}}=4$ at
$\frac{t}{\tau_{\mathcal{L}}}\simeq 2$ for (a) n=2, and (b) n=3 in Eq. 27

As the last and most stringent test in the context of the two-point
diagnostics, higher-order structure functions for the velocity increments are
compared with the ground-truth filtered-DNS results. The second and third
order structures functions are computed based on the following relation

$\displaystyle\langle(\delta_{r}\bar{u}_{L})^{n}\rangle=\langle\big{[}\bar{u}_{L}(x+re_{L})-\bar{u}_{L}(x)\big{]}^{n}\rangle,\quad
n=2,3.$ (27)

To get the structure functions, first, we compute the velocity fields at a
certain time for each model then we shift them based on ${u}_{L}(x+re_{L})$,
and finally, we do the filtering operation on each obtained field. The spatial
shift is related to the filter size as $r={\mathcal{L}_{\delta}}$. Comparing
the results for the second and third-order structure functions in Fig. 10 for
$\frac{t}{\tau_{\mathcal{L}}}\simeq 2$, reveals that all models are almost
preserving the second-order structures, however, in the third-order structures
we have more deviation for the Smagorinsky-based models. The new proposed
model is more successful in the prediction of the third-order structures and
its trend in the LES region.

##  V Concluding Remarks

The present study introduces a novel dynamic nonlocal turbulence model for the
isotropic turbulent flows, which can be applied for both LES and VLES
purposes. The model has been developed based on the fractional Laplacian
derivative of resolved velocity field, and a unique dynamic procedure defines
the model constant in the divergence form. The effect of the fractional-order,
$Re$ number, and characteristic filter size in LES and VLES cases are
scrutinized using multiple high-fidelity and well-resolved DNS datasets
belongs to forced homogeneous isotropic turbulence. Final decisions were made
based on the ensemble-averaged quantities gathered from ten separate three-
dimensional snapshots distributed over enough turnover times. The optimum
fractional-order for each scenario is chosen when the maximum correlation is
obtained for the ground-truth and predicted stresses of the model. The
obtained results indicate that the relation between the filter size and
fractional-order is closely obeying a power-law form for all the assessed $Re$
numbers. Also, there is a direct relationship between the $Re$ number and the
fractional-order. Considering all the attained statistical results, two other
data-driven straightforward methods for the determination of the fractional-
order suggested and tested successfully.

Our analysis included both a priori and a posteriori assessments. In the first
one, we showed that there is a higher correlation between the results of the
proposed method and the ground-truth DNS results in comparison to other
conventional LES models. Also, the capability of the model in the prediction
of the back-scatter is discussed in different filter sizes. In the a
posteriori assessments we tested the model performance in long-time
integration in the context of a real N-S solvers, and its ability in capturing
the large-scale coherent structures examined. Analysis were performed in a
decaying isotropic turbulence scenario, and the new D-FSGS, static Smagorinsky
and dynamic Smagorinsky models implemented separately. The results show that
the new model is more successful in the prediction of the resolved turbulent
kinetic energy in both LES and VLES studies. To test the models’ performance
in the long-time integration and its numerical stability, kinetic energy
spectra are compared with the filtered-DNS results. Finally, two-point
diagnostics were accomplished to compare different models’ performance in the
context of preserving second and third-order structures. The results
demonstrate that the new model behaves better in preserving the high-order
structures than the conventional ones.

This study shows the promising potential of using bigger filter sizes rather
than the conventional LES filter sizes. This important characteristic can be
utilized to compensate the fractional model’s high-cost demand. Moreover,
accelerated evaluation methods for the fractional operators including
learning-based approaches and fast solvers can be leveraged to pave the road
for the real and more practical applications of nonlocal models.

## Acknowledgement

This work was supported by the ARO YIP award (W911NF-19-1-0444), and partially
by the MURI/ARO award (W911NF-15-1-0562) and the NSF award (DMS-1923201). The
HPC resources and services were provided by the Institute for Cyber-Enabled
Research (ICER) at Michigan State University.

## Data Availability

The data that support the findings of this research are available from the
corresponding author upon reasonable request.

## References

  * Pope (2001) S. B. Pope, “Turbulent flows,” (2001).
  * Sagaut (2006) P. Sagaut,  _Large eddy simulation for incompressible flows: an introduction_ (Springer Science & Business Media, 2006).
  * Meneveau and Katz (2000) C. Meneveau and J. Katz, “Scale-invariance and turbulence models for large-eddy simulation,” Annual Review of Fluid Mechanics 32, 1–32 (2000).
  * Germano  _et al._ (1991) M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, “A dynamic subgrid-scale eddy viscosity model,” Physics of Fluids A: Fluid Dynamics 3, 1760–1765 (1991).
  * Moin  _et al._ (1991) P. Moin, K. Squires, W. Cabot, and S. Lee, “A dynamic subgrid-scale model for compressible turbulence and scalar transport,” Physics of Fluids A: Fluid Dynamics 3, 2746–2757 (1991).
  * Vasilyev, Lund, and Moin (1998) O. V. Vasilyev, T. S. Lund, and P. Moin, “A general class of commutative filters for les in complex geometries,” Journal of computational physics 146, 82–104 (1998).
  * Pope (2004) S. B. Pope, “Ten questions concerning the large-eddy simulation of turbulent flows,” New journal of Physics 6, 35 (2004).
  * Wilczek, Vlaykov, and Lalescu (2017) M. Wilczek, D. G. Vlaykov, and C. C. Lalescu, “Emergence of non-gaussianity in turbulence,” in  _Progress in Turbulence VII_ (Springer, 2017) pp. 3–9.
  * Akhavan-Safaei, Seyedi, and Zayernouri (2020) A. Akhavan-Safaei, S. H. Seyedi, and M. Zayernouri, “Anomalous features in internal cylinder flow instabilities subject to uncertain rotational effects,” Physics of Fluids 32, 094107 (2020).
  * Bradshaw (1973) P. Bradshaw, “Agardograph, no. 169,” Nato Science and Technology Organisation, USA (1973).
  * Spencer and Rivlin (1958) A. J. M. Spencer and R. S. Rivlin, “The theory of matrix polynomials and its application to the mechanics of isotropic continua,” Archive for rational mechanics and analysis 2, 309–336 (1958).
  * Spencer and Rivlin (1959) A. J. M. Spencer and R. S. Rivlin, “Further results in the theory of matrix polynomials,” Archive for rational mechanics and analysis 4, 214–230 (1959).
  * Lumley (1970) J. Lumley, “Toward a turbulent constitutive relation,” J. Fluid Mech 41, 413–434 (1970).
  * Pope (1975) S. Pope, “A more general effective-viscosity hypothesis,” Journal of Fluid Mechanics 72, 331–340 (1975).
  * Suzuki  _et al._ (2021a) J. L. Suzuki, E. Kharazmi, P. Varghaei, M. Naghibolhosseini, and M. Zayernouri, “Anomalous nonlinear dynamics behavior of fractional viscoelastic beams,” Journal of Computational and Nonlinear Dynamics 16, 111005 (2021a).
  * Suzuki  _et al._ (2021b) J. L. Suzuki, T. G. Tuttle, S. Roccabianca, and M. Zayernouri, “A data-driven memory-dependent modeling framework for anomalous rheology: Application to urinary bladder tissue,” Fractal and Fractional 5, 223 (2021b).
  * Naghibolhosseini and Long (2018) M. Naghibolhosseini and G. R. Long, “Fractional-order modelling and simulation of human ear,” International Journal of Computer Mathematics 95, 1257–1273 (2018).
  * Naghibolhosseini (2015) M. Naghibolhosseini,  _Estimation of outer-middle ear transmission using DPOAEs and fractional-order modeling of human middle ear_ (City University of New York, 2015).
  * Suzuki  _et al._ (2021c) J. Suzuki, Y. Zhou, M. D’Elia, and M. Zayernouri, “A thermodynamically consistent fractional visco-elasto-plastic model with memory-dependent damage for anomalous materials,” Computer Methods in Applied Mechanics and Engineering 373, 113494 (2021c).
  * Suzuki  _et al._ (2016) J. Suzuki, M. Zayernouri, M. Bittencourt, and G. Karniadakis, “Fractional-order uniaxial visco-elasto-plastic models for structural analysis,” Computer Methods in Applied Mechanics and Engineering 308, 443–467 (2016).
  * D’Elia  _et al._ (2020) M. D’Elia, Q. Du, C. Glusa, M. Gunzburger, X. Tian, and Z. Zhou, “Numerical methods for nonlocal and fractional models,” arXiv preprint arXiv:2002.01401 (2020).
  * Kharazmi and Zayernouri (2019) E. Kharazmi and M. Zayernouri, “Fractional sensitivity equation method: Application to fractional model construction,” Journal of Scientific Computing 80, 110–140 (2019).
  * Samiee, Zayernouri, and Meerschaert (2019) M. Samiee, M. Zayernouri, and M. M. Meerschaert, “A unified spectral method for FPDEs with two-sided derivatives; part i: a fast solver,” Journal of Computational Physics 385, 225–243 (2019).
  * Zhou  _et al._ (2020) Y. Zhou, J. L. Suzuki, C. Zhang, and M. Zayernouri, “Implicit-explicit time integration of nonlinear fractional differential equations,” Applied Numerical Mathematics 156, 555–583 (2020).
  * Seyedi, Nemati Saray, and Ramazani (2019) S. H. Seyedi, B. Nemati Saray, and A. Ramazani, “High-accuracy multiscale simulation of three-dimensional squeezing carbon nanotube-based flow inside a rotating stretching channel,” Mathematical Problems in Engineering 2019 (2019).
  * Lischke, Zayernouri, and Zhang (2019) A. Lischke, M. Zayernouri, and Z. Zhang, “Spectral and spectral element methods for fractional advection–diffusion–reaction equations,” Numerical Methods 157 (2019).
  * Seyedi (2018) S. Seyedi, “Multiresolution solution of burgers equation with b-spline wavelet basis,” arXiv preprint arXiv:1812.10117 (2018).
  * Zayernouri and Karniadakis (2014) M. Zayernouri and G. E. Karniadakis, “Fractional spectral collocation method,” SIAM Journal on Scientific Computing 36, A40–A62 (2014).
  * Kharazmi, Zayernouri, and Karniadakis (2017) E. Kharazmi, M. Zayernouri, and G. E. Karniadakis, “Petrov–Galerkin and spectral collocation methods for distributed order differential equations,” SIAM Journal on Scientific Computing 39, A1003–A1037 (2017).
  * Zayernouri and Karniadakis (2013) M. Zayernouri and G. E. Karniadakis, “Fractional Sturm–Liouville eigen-problems: theory and numerical approximation,” Journal of Computational Physics 252, 495–517 (2013).
  * Kharazmi and Zayernouri (2018) E. Kharazmi and M. Zayernouri, “Fractional pseudo-spectral methods for distributed-order fractional pdes,” International Journal of Computer Mathematics 95, 1340–1361 (2018).
  * Hinze, Sonnenberg, and Builtjes (1974) J. O. Hinze, R. Sonnenberg, and P. J. H. Builtjes, “Memory effect in a turbulent boundary-layer flow due to a relatively strong axial variation of the mean-velocity gradient,” Applied Scientific Research 29, 1–13 (1974).
  * Kraichnan (1964) R. H. Kraichnan, “Direct-interaction approximation for shear and thermally driven turbulence,” The Physics of Fluids 7, 1048–1062 (1964).
  * Chen (2006) W. Chen, “A speculative study of 2/ 3-order fractional Laplacian modeling of turbulence: Some thoughts and conjectures,” Chaos: An Interdisciplinary Journal of Nonlinear Science 16, 023126 (2006).
  * Epps and Cushman-Roisin (2018) B. P. Epps and B. Cushman-Roisin, “Turbulence modeling via the fractional Laplacian,” arXiv preprint arXiv:1803.05286 (2018).
  * Egolf and Hutter (2017) P. W. Egolf and K. Hutter, “Fractional turbulence models,” in  _Progress in Turbulence VII_ (Springer, 2017) pp. 123–131.
  * Hamba (1995) F. Hamba, “An analysis of nonlocal scalar transport in the convective boundary layer using the green’s function,” Journal of Atmospheric Sciences 52, 1084–1095 (1995).
  * Hamba (2005) F. Hamba, “Nonlocal analysis of the reynolds stress in turbulent shear flow,” Physics of Fluids 17, 115102 (2005).
  * Egolf and Kutter (2020) P. W. Egolf and K. Kutter, “Nonlinear, nonlocal and fractional turbulence,” Graduate Studies in Mathematics. Springer, (2020).
  * Samiee, Akhavan-Safaei, and Zayernouri (2020) M. Samiee, A. Akhavan-Safaei, and M. Zayernouri, “A fractional subgrid-scale model for turbulent flows: Theoretical formulation and a priori study,” Physics of Fluids 32, 055102 (2020).
  * Samiee, Akhavan-Safaei, and Zayernouri (2021) M. Samiee, A. Akhavan-Safaei, and M. Zayernouri, “Tempered fractional LES modeling,” Journal of Fluid Mechanics, in press, (arXiv preprint arXiv:2103.01481) (2021).
  * Laval, Dubrulle, and Nazarenko (2001) J. Laval, B. Dubrulle, and S. Nazarenko, “Nonlocality and intermittency in three-dimensional turbulence,” Physics of Fluids 13, 1995–2012 (2001).
  * Akhavan-Safaei, Samiee, and Zayernouri (2021) A. Akhavan-Safaei, M. Samiee, and M. Zayernouri, “Data-driven fractional subgrid-scale modeling for scalar turbulence: A nonlocal les approach,” Journal of Computational Physics 446, 110571 (2021).
  * Akhavan-Safaei and Zayernouri (2021) A. Akhavan-Safaei and M. Zayernouri, “A nonlocal spectral transfer model and new scaling law for scalar turbulence,” arXiv preprint arXiv:2111.06540 (2021).
  * Keith, Khristenko, and Wohlmuth (2021) B. Keith, U. Khristenko, and B. Wohlmuth, “A fractional pde model for turbulent velocity fields near solid walls,” Journal of Fluid Mechanics 916 (2021).
  * Di Leoni  _et al._ (2021) P. C. Di Leoni, T. A. Zaki, G. Karniadakis, and C. Meneveau, “Two-point stress–strain-rate correlation structure and non-local eddy viscosity in turbulent flows,” Journal of Fluid Mechanics 914 (2021).
  * Milovanov and Rasmussen (2014) A. V. Milovanov and J. J. Rasmussen, “A mixed soc-turbulence model for nonlocal transport and lévy-fractional fokker–planck equation,” Physics Letters A 378, 1492–1500 (2014).
  * Ali (2014) H. Ali, “Theory for the rotational deconvolution model of turbulence with fractional regularization,” Applicable Analysis 93, 339–355 (2014).
  * Gunzburger, Jiang, and Xu (2018) M. Gunzburger, N. Jiang, and F. Xu, “Analysis and approximation of a fractional Laplacian-based closure model for turbulent flows and its connection to richardson pair dispersion,” Computers & Mathematics with Applications 75, 1973–2001 (2018).
  * Suzuki  _et al._ (2021d) J. Suzuki, M. Gulian, M. Zayernouri, and M. D’Elia, “Fractional modeling in action: A survey of nonlocal models for subsurface transport, turbulent flows, and anomalous materials,” arXiv preprint arXiv:2110.11531 (2021d).
  * Smagorinsky (1963) J. Smagorinsky, “General circulation experiments with the primitive equations: I. the basic experiment,” Monthly weather review 91, 99–164 (1963).
  * Lund and Kaltenbach (1995) T. S. Lund and H.-J. Kaltenbach, “Experiments with explicit filtering for les using a finite-difference method,” (1995).
  * Najjar and Tafti (1996) F. M. Najjar and D. K. Tafti, “Study of discrete test filters and finite difference approximations for the dynamic subgrid-scale stress model,” Physics of Fluids 8, 1076–1088 (1996).
  * Piomelli and Liu (1995) U. Piomelli and J. Liu, “Large-eddy simulation of rotating channel flows using a localized dynamic model,” Physics of fluids 7, 839–848 (1995).
  * Ghosal and Rogers (1997) S. Ghosal and M. M. Rogers, “A numerical study of self-similarity in a turbulent plane wake using large-eddy simulation,” Physics of Fluids 9, 1729–1739 (1997).
  * Chen, Orszag, and Staroselsky (2010) H. Chen, S. A. Orszag, and I. Staroselsky, “Macroscopic description of arbitrary knudsen number flow using boltzmann–bgk kinetic theory. part 2,” Journal of fluid mechanics 658, 294–309 (2010).
  * Girimaji (2007) S. S. Girimaji, “Boltzmann kinetic equation for filtered fluid turbulence,” Physical review letters 99, 034501 (2007).
  * Lischke  _et al._ (2020) A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth, _et al._ , “What is the fractional Laplacian? a comparative review with new results,” Journal of Computational Physics 404, 109009 (2020).
  * Lilly (1992) D. K. Lilly, “A proposed modification of the germano subgrid-scale closure method,” Physics of Fluids A: Fluid Dynamics 4, 633–635 (1992).
  * Akhavan-Safaei and Zayernouri (2020) A. Akhavan-Safaei and M. Zayernouri, “A parallel integrated computational-statistical platform for turbulent transport phenomena,” arXiv preprint arXiv:2012.04838 (2020).
  * Meneveau (1994) C. Meneveau, “Statistics of turbulence subgrid-scale stresses: Necessary conditions and experimental tests,” Physics of Fluids 6, 815–833 (1994).
  * Piomelli  _et al._ (1991) U. Piomelli, W. H. Cabot, P. Moin, and S. Lee, “Subgrid-scale backscatter in turbulent and transitional flows,” Physics of Fluids A: Fluid Dynamics 3, 1766–1771 (1991).