# A physics-inspired alternative to spatial filtering for large-eddy
simulations of turbulent flows

Perry L. Johnson<EMAIL_ADDRESS>Department of Mechanical and Aerospace
Engineering, University of California, Irvine, CA, USA

###### Abstract

Large-eddy simulations (LES) are widely-used for computing high Reynolds
number turbulent flows. Spatial filtering theory for LES is not without its
shortcomings, including how to define filtering for wall-bounded flows,
commutation errors for non-uniform filters, and extensibility to flows with
additional complexity, such as multiphase flows. In this paper, the theory for
LES is reimagined using a coarsening procedure that imitates nature. This
physics-inspired coarsening (PIC) approach is equivalent to Gaussian filtering
for single-phase wall-free flows but opens up new insights for modeling even
in that simple case. Boundaries and nonuniform resolution can be treated
seemlessly in this framework without commutation errors, and the divergence-
free condition is retained for incompressible flows. An alternative to the
Germano identity is introduced and used to define a dynamic procedure without
the need for a test filter. Potential extensions to more complex physics are
briefly discussed.

##  1 Introduction

Because of the inherently multiscale nature of turbulence, the fine grids and
small time steps required for direct numerical simulation (DNS) of turbulent
flows are prohibitively expensive for many higher Reynolds number flow
applications. Large eddy simulation (LES) has become a popular technique for
computing and predicting turbulent flows. In LES, the computational grid may
be significantly coarser than the requirements of DNS, directly resolving only
the large-scale features of the flow. Sub-grid models are introduced to
approximate the net effect of unresolved fluctuations on the large-scale
dynamics.

The theory for LES is based on reducing the number of computational degrees of
freedom (DoFs) necessary for a digital representation of a turbulent flow,
i.e. increasing the required grid spacing, typically using a low-pass spatial
filter (Leonard, 1975; Germano, 1992; Sagaut, 2006; Moser et al., 2021). The
practical significance of spatial filtering theory for LES is not necessarily
obvious, because most approaches assume an implicit filter. That is, an
explicit calculation of a spatial filter is not necessary or even commonly
employed in LES practice. What the theory of spatial filtering does provide is
a definition of what the simulation aims to reproduce, and more importantly, a
partial differential equation for the evolution of the resolved flow features,
including a mathematical expression for the unclosed term(s) that need to be
modeled. The definitions provided by spatial filtering enable the direct
calculation of unclosed terms from DNS data, allowing for the a priori testing
of candidate models (Borue & Orszag, 1998) as well as the theoretical
derivation or justification of many proposed models (Clark et al., 1979;
Bardina et al., 1980; Verstappen, 2011; Rozema et al., 2015). Perhaps the most
notable fruit of spatial filtering theory for LES practice is the widely-used
dynamic procedure for determining model coefficients on the fly using an
explicit test filter (Moin et al., 1991; Germano et al., 1991; Lilly, 1992;
Vreman et al., 1994; Meneveau et al., 1996; Bou-Zeid et al., 2005). Dynamic
models are based on the Germano identity, which relates fields filtered at two
different filter widths (Germano, 1992).

However, spatial filtering LES theory is not without drawbacks, which have
motivated some proposed modifications and alternative formulations. The
optimal LES approach introduced by Langford & Moser (1999) addresses the
inherent loss of information by invoking conditional averages and seeking a
formulation for predicting single-time, multi-point statistics while
minimizing short-time errors in general. Ultimately, it still relies on
spatial filtering and assumes that no history information may be used for
modeling, e.g., by Lagrangian averaging (Meneveau et al., 1996; Bou-Zeid et
al., 2005) or evolving a sub-grid kinetic energy equation (Kim & Menon, 1996).

The self-conditioned fields approach introduced by Fox (2003) and Pope (2010)
redefines the objective of LES as the (unfiltered) velocity conditioned on the
filtered velocity field (or another lower dimensional representation). This
conception of LES has many positive characteristics but makes direct (a
priori) testing of models practically impossible. Self-conditioned fields LES
also introduces additional closure terms pertaining to gradients with respect
to conditioning variables (for which no models are proposed beyond ignoring
the terms for expediency), and it largely does not provide much insight into
how to construct models, mostly serving to justify existing models.

One significant difficulty with spatial filtering theory for LES is that
nonuniform spatial filters do not commute with spatial differentiation,
leading to commutation errors and loss of divergence-free velocity fields in
the case of incompressible flows (Ghosal & Moin, 1995; Langford & Moser, 2001;
Yalla et al., 2021). In practice, commutation errors are typically neglected
and a divergence-free condition is imposed without justification. Another
difficulty with spatial filtering arises close to a boundary, where the filter
operator requires flow information from outside the fluid domain (Drivas &
Nguyen, 2018), unless a nonuniform filter size is used (Bose & Moin, 2014)
leading to the commutation errors described above. The filtered field may also
be conceived of as the solution to an elliptic partial differential equation
(Germano, 1986a, b; Bull & Jameson, 2016), in which case the choice of
boundary conditions may alleviate this particular difficulty (Bae & Lozano-
Duran, 2017). However, common LES treatments near the wall simply revert to
RANS-based techniques without full justification in terms of spatial
filtering. This is true for hybrid RANS-LES models (Piomelli & Balaras, 2002;
Frohlich & von Terzi, 2008; Spalart, 2009; Mockett et al., 2012) as well as
wall-modeled LES (Larsson et al., 2016; Bose & Park, 2018). The self-
conditioned fields approach does manage to circumvent commutator errors
altogether and to provide a well-defined near-wall behavior matching the no-
slip boundary condition (Pope, 2010), though these positive characteristics
are counter-balanced by the challenges noted in the previous paragraph.

With interest growing in performing LES-like calculations with increasing
physical complexity, spatial filtering theory is encountering more challenges.
In multiphase flows, for instance, direct application of spatial filtering
blurs the interface (Labourasse et al., 2007), unnecessarily losing its
distinctiveness even though sharp interface methods exist for careful, robust
treatment such discontinuities. It is desirable to have an LES theory that
allows for retaining sharp interfaces while removing features such as small-
scale ripples that cannot be resolved on a coarse grid (Tryggvason & Lu,
2020). One alternative is the dual-scale approach of Herrmann (2013), but this
requires DNS-like resolution for the phase indicator field (or volume fraction
field) and a super-resolution enrichment for the velocity (Herrmann et al.,
2018).

In the simulation of turbulent flows laden with small particles, including
unresolved droplets or bubbles in multiphase flows, it is common to avoid the
need to resolve the flow around a particle (i.e., computationally resolve the
proper boundary or interface conditions at the particle surface) using a
Lagrangian tracking approach. A closure such as a drag law must be introduced.
When the particle is much smaller than the Kolmogorov scale, the Reynolds
number based on its relative velocity is small, and the particle is not near a
flow boundary, the Lagrangian tracking approach can be quite precise (Maxey &
Riley, 1983; Balachandar & Eaton, 2010). In other more complex scenarios,
Lagrangian tracking represents an approximate treatment motivated by
computational tractability not unlike spatial filtering for LES. In fact, if
the goal of reducing computational DoF is indeed central to LES, then the
mathematical theory of LES should be fundamentally compatible with an
Eulerian-Lagrangian approach to particle-laden flows in a way that spatial
filtering is not. This includes potential hybrid approaches to multiphase
flows that combine the direct resolution of large-scale interface features
with a Lagrangian tracking approach for small-scale features (Kim & Moin,
2020).

While LES may be practiced in isolation from specific concerns of a consistent
framework, a specific definition of that which an LES aspires to accurately
reproduce is required for advanced techniques such as data-driven closure
(Sarghini et al., 2003; Moreau et al., 2006; Vollant et al., 2017; Gamahara &
Hattori, 2017; Wang et al., 2018; Yang et al., 2019; Zhou et al., 2019; Cheng
et al., 2019; Beck et al., 2019; Sirignano et al., 2020; Xie et al., 2020b, a;
Yuan et al., 2020; Park & Choi, 2021; Bode et al., 2021; Freund & Ferrante,
2021; Stoffer et al., 2021; Portwood et al., 2021; Wang et al., 2021; Prakash
et al., 2021; Duraisamy, 2021) and super-resolution enrichment (Scotti &
Meneveau, 1999; Stolz & Adams, 1999; Domaradzki & Loh, 1999; Milano &
Koumoutsakos, 2002; Leonard, 2016; Ghate & Lele, 2017; Maulik & San, 2017;
Bassenne et al., 2019; Wang et al., 2019; Ghate & Lele, 2020; Liu et al.,
2020; Kim et al., 2021). For example, without a clear definition of what an
LES solution should represent, one cannot train a neural network to serve as a
sub-grid closure in a robust way. In this sense, the particular relationship
between an LES solution and the fully-resolved flow must be afforded more
direct attention in the age of data. This is especially true for more complex
scenarios with solid boundaries, interfaces, and multiphysics; where pure
spatial filtering is not likely the optimal approach.

The scope of this paper is as follows. This paper introduces the concept of
physics-inspired coarsening (PIC) as a proposed alternative to spatial
filtering as a basis for LES. In addition to laying the groundwork of PIC
theory, results for the simplest of turbulent flows are demonstrated with a
priori and a posteriori testing of representative PIC-based models. Beyond
that, the paper explores the extensibility of PIC theory for complex flows to
motivate future work in various directions. The basic theory of PIC is
established in §2. Following that, §3 explores in detail single-phase
unbounded flows with uniform resolution, for which PIC is equivalent to
Gaussian filtering. Even with this equivalence, PIC theory provides new
insights into energy cascade physics and model development. Preliminary models
based on these insights are demonstrated with a priori and a posteriori
testing. In §4, PIC theory is expanded to include a number of more complex
effects such as anisotropic and nonuniform resolution, heat and mass transfer,
and flows with solid boundaries and multiphase interfaces. Conclusions are
drawn in §5.

##  2 Physics-inspired coarsening (PIC)

The underlying philosophy of PIC is to view the removal of DoFs from a flow
more as an imitation of nature than as an image processing trick. As a primary
example, the viscosity, $\mu$, provides a natural mechanism that prevents the
creation of motions smaller than the Kolmogorov microscale,
$\eta=\nu^{3/4}\epsilon^{-1/4}$, where $\nu=\mu/\rho$ is the kinematic
viscosity, $\rho$ is the fluid mass density, and $\epsilon$ is the turbulent
dissipation rate. This natural mechanism may be imitated for LES theory by
conjuring an artificial viscous process to further smooth a frozen snapshot of
turbulence, removing small motions below a cutoff length scale, $\ell>\eta$.
The artificial physics-inspired process evolves in a pseudo-time that is
independent of physical time. The overall idea of physics-inspired coarsening
is illustrated in Figure 1 and described in detail in the remainder of this
section.

###  2.1 Navier-Stokes equation

For the present purposes, an incompressible flow of a Newtonian fluid is
considered. A velocity vector field, $\mathbf{u}=\mathbf{u}(\mathbf{x},t)$,
evolves as a function of space, $\mathbf{x}$, and time, $t$, according to the
Navier-Stokes equation,

$\frac{\partial u_{i}}{\partial t}+\frac{\partial}{\partial
x_{j}}\left[u_{i}u_{j}+p\delta_{ij}-\nu\left(\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\right]=0,\hskip
43.36464pt\frac{\partial u_{j}}{\partial x_{j}}=0.$ (1)

The pressure (divided by density), $p(\mathbf{x},t)$, enforces the divergence-
free condition. A Poisson equation for the pressure follows from the
divergence of the Navier-Stokes equation,

$\nabla^{2}p=-\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial
u_{j}}{\partial x_{i}}.$ (2)

Thus, a fully-resolved simulation (DNS) evolves Eq. (1) together with Eq. (2)
in physical time with pseudo-time fixed at zero, $\hat{t}=0$, indicated by the
green arrow in Figure 1.

Figure 1: For physics-inspired coarsening for LES, a generalized velocity
field is defined in as a function of physical time ($t$) and pseudo-time
($\hat{t}\sim\text{grid-size}^{2}$). The generalized Navier-Stokes equation
for the LES velocity in physical time requires a closure model.

###  2.2 Auxiliary Stokes equation

For the physics-inspired approach, a generalized velocity field,
$\mathbf{w}=\mathbf{w}(\mathbf{x},t;\hat{t})$, is defined as a function of
both time, $t$, and a pseudo-time, $\hat{t}$. The generalized velocity field
at $\hat{t}=0$ corresponds to the physical velocity field,

$\mathbf{w}(\mathbf{x},t;0)=\mathbf{u}(\mathbf{x},t).$ (3)

As $\hat{t}$ increases, small-scale motions are removed from the generalized
velocity field. Thus, the pseudo-time is an indicator of the resolution length
scale of the generalized velocity. To accomplish this, an auxiliary Stokes
equation may be introduced to govern the pseudo-time evolution,

$\frac{\partial w_{i}}{\partial\hat{t}}+\frac{\partial}{\partial
x_{j}}\left[\hat{p}\delta_{ij}-\hat{\nu}\left(\frac{\partial w_{i}}{\partial
x_{j}}+\frac{\partial w_{j}}{\partial x_{i}}\right)\right]=0,\hskip
43.36464pt\frac{\partial w_{j}}{\partial x_{j}}=0.$ (4)

The nonlinear term of the Navier-Stokes equation ($u_{i}u_{j}$) is removed for
Eq. (4), so there is no energy cascade in pseudo-time and hence no resupply of
energy to small-scales dissipated by pseudo-viscosity,
$\hat{\nu}(\mathbf{x})$. The length scale at which the flow is smoothed may be
heuristically estimated as $\ell\sim\sqrt{\hat{\nu}\hat{t}}$. A pseudo-
pressure, $\hat{p}(\mathbf{x},t;\hat{t})$, is also introduced to enforce a
divergence-free condition for the generalized velocity at all pseudo-times.
The divergence of the auxiliary Stokes equation is a Poisson equation for the
pseudo-presssure,

$\nabla^{2}\hat{p}=\frac{\partial^{2}}{\partial x_{i}\partial
x_{j}}\left[\hat{\nu}\left(\frac{\partial w_{i}}{\partial
x_{j}}+\frac{\partial w_{j}}{\partial x_{i}}\right)\right].$ (5)

Thus, a DNS may be artificially coarsened by, instead of applying a spatial
filter, advancing Eq. (4) in pseudo-time together with Eq. (5) having initial
conditions set by Eq. (3). This is represented in Figure 1 by the orange
arrows. In the case of a moving reference frame, the coordinate system only
advances in physical time. It is frozen in pseudo-time by definition, so that
Galilean invariance is satisfied. Note that the pseudo-viscosity is distinct
from the idea of an eddy viscosity (or turbulent viscosity) often used to
model turbulent flows. The pseudo-viscosity acts only in pseudo-time and
simply enables the coarsening procedure that is an alternative to spatial
filtering. An eddy viscosity acts in physical time and attempts to approximate
the evolution of a coarsened velocity field in physical time, which is
considered next.

###  2.3 Generalized Navier-Stokes equation

Explicit filtering is not typically used in practice for LES. Similarly, it is
not envisioned at present that PIC-based LES practice will (necessarily)
include the numerical solution of an auxiliary Stokes equation. The
significance of Eq. (4), rather, is its implication for the evolution of the
generalized velocity in physical time, keeping $\hat{t}\sim\ell^{2}$ fixed
based on the grid resolution. This is illustrated by the purple arrow in
Figure 1 and corresponds to establishing the effective equations for solving a
coarsened representation of the flow, which may in practice be done on a
coarser grid.

The outcome of the auxiliary Stokes equation ensures that the Navier-Stokes
equation must be altered to accurately describe the evolution of the coarsened
flow field in physical time. Thus, an effective evolution equation for LES may
be derived by finding the appropriate residual force, $\mathbf{F}$, such that,

$\frac{\partial w_{i}}{\partial t}+\frac{\partial}{\partial
x_{j}}\left[w_{i}w_{j}+\phi\delta_{ij}-\nu\left(\frac{\partial w_{i}}{\partial
x_{j}}+\frac{\partial w_{j}}{\partial x_{i}}\right)\right]=F_{i},\hskip
43.36464pt\frac{\partial w_{j}}{\partial x_{j}}=0.$ (6)

When $\hat{t}=0$, the generalized velocity does satisfy Eq. (1) because of the
initial condition, Eq. (3). This means that, $\mathbf{F}(\mathbf{x},t;0)=0$.

Because the auxiliary Stokes equation globally conserves momentum in the
absence of fluxes through domain boundaries, the residual force must likewise
conserve momentum, and it may be written as the divergence of a residual
stress tensor,

$F_{i}=\frac{\partial\sigma_{ij}}{\partial x_{j}}.$ (7)

Thus, the evolution of the generalized velocity in physical time is given by
the generalized Navier-Stokes equation,

$\frac{\partial w_{i}}{\partial t}+\frac{\partial}{\partial
x_{j}}\left[w_{i}w_{j}+\sigma_{ij}+\phi\delta_{ij}-\nu\left(\frac{\partial
w_{i}}{\partial x_{j}}+\frac{\partial w_{j}}{\partial
x_{i}}\right)\right]=0,\hskip 43.36464pt\frac{\partial w_{j}}{\partial
x_{j}}=0.$ (8)

Note that the divergence-free constraint on the auxiliary Stokes equation
implies a divergence-free condition for the physical time evolution as well.
This is enforced by a generalized pressure, $\phi(\mathbf{x},t;\hat{t})$, that
is distinct from the pseudo-pressure, $\hat{p}$, which acts only during
pseudo-time evolution. The generalized pressure may be found by solution of
the following Poisson equation,

$\nabla^{2}\phi=-\frac{\partial w_{i}}{\partial x_{j}}\frac{\partial
w_{j}}{\partial x_{i}}-\frac{\partial^{2}\sigma_{ij}}{\partial x_{i}\partial
x_{j}}.$ (9)

Equation (8), together with Eq. (9) is unclosed in general. The residual
stress tensor field must be approximated in terms of the generalized velocity
field in order to solve for its evolution. More specifically, only the
deviatoric component of $\boldsymbol{\sigma}$ is material to the closure
problem, because any isotropic component may be lumped in with the generalized
pressure and found via Eq. (9). For consistency, it is required that the
generalized velocity, $\mathbf{w}$, evolve according to the Navier-Stokes
equation for $\hat{t}=0$. That is, Eq. (8) must be identical to Eq. (1) when
$\hat{t}=0$. Thus,

$\phi(\mathbf{x},t;0)=p(\mathbf{x},t),\hskip
43.36464pt\sigma_{ij}(\mathbf{x},t;0)=0.$ (10)

###  2.4 Residual stress tensor equation

The significance of the above theory is its implications for the mathematical
definition of the residual stress tensor, which is found by requiring
consistency between the pseudo-time and physical time evolution of the
generalized velocity field sketched in Figure 1. Specifically,
$\boldsymbol{\sigma}$ may be defined by asserting that advancing the
generalized velocity field from $(t_{0},\hat{t}_{0})$ to $(t_{f},\hat{t}_{f})$
must be independent of the path in $t$-$\hat{t}$ space. Breaking the evolution
in $t$-$\hat{t}$ space into infinitesimal increments, this requirement may be
written as,

$\frac{\partial}{\partial\hat{t}}\left(\frac{\partial w_{i}}{\partial
t}\right)=\frac{\partial}{\partial t}\left(\frac{\partial
w_{i}}{\partial\hat{t}}\right),$ (11)

which is simply the symmetry of mixed partial derivatives. More detailed
treatment of the implications of this constraint for the form of the residual
stress tensor will be given later.

###  2.5 Kinetic energy equations

The scale-wise dynamics of kinetic energy is crucial to the representation of
turbulent flows. The generalized kinetic energy (per unit mass) in PIC is
given by $\tfrac{1}{2}w_{i}w_{i}$. First, consider how kinetic energy evolves
in pseudo-time. Multiplying the auxiliary Stokes equation, Eq. (4), by $w_{i}$
and using a product rule on the viscous term,

$\frac{\partial(\tfrac{1}{2}w_{i}w_{i})}{\partial\hat{t}}+\frac{\partial}{\partial
x_{j}}\left[\hat{p}w_{j}-2\hat{\nu}w_{i}S_{ij}\right]=-2\hat{\nu}S_{ij}S_{ij},$
(12)

where $S_{ij}=\tfrac{1}{2}(\partial w_{i}/\partial x_{j}+\partial
w_{j}/\partial x_{i})$ is the generalized strain-rate tensor. As pseudo-time
increases, kinetic energy is removed by the pseudo-viscosity at a rate
$2\hat{\nu}S_{ij}S_{ij}\geq 0$.

The dynamics of generalized kinetic energy in physical time (at fixed, finite
pseudo-time) can be written by multiplying the generalized Navier-Stokes
equation, Eq. (8) by $w_{i}$,

$\frac{\partial(\tfrac{1}{2}w_{i}w_{i})}{\partial t}+\frac{\partial}{\partial
x_{j}}\left[\frac{1}{2}w_{i}w_{i}w_{j}+w_{i}\sigma_{ij}+\phi w_{j}-2\nu
w_{i}S_{ij}\right]=\sigma_{ij}S_{ij}-2\nu S_{ij}S_{ij}.$ (13)

In addition to direct viscous dissipation, $2\nu S_{ij}S_{ij}\geq 0$, the
generalized kinetic energy is removed at a rate $\Pi=-\sigma_{ij}S_{ij}$. Note
that $\Pi$ may be positive or negative in general, but that the expectation
for turbulence dynamics is a preference toward $\Pi>0$ representing the
kinetic energy cascade. That is, $\Pi$ represents the rate at which kinetic
energy associated with motions larger than $\ell\sim\sqrt{\hat{\nu}\hat{t}}$
is passed to motions smaller than $\ell$. Thus, an important consideration for
any closure model for $\sigma_{ij}$ is its ability to remove energy at the
proper rate to mimic the energy cascade.

###  2.6 Juxtaposition with spatial filtering

The traditional approach to LES theory is based on developing equations for
the filtered velocity field (Leonard, 1975; Germano, 1992),

$\overline{u}_{i}^{\ell}(\mathbf{x},t)=\iiint
G^{\ell}(\mathbf{r})u_{i}(\mathbf{x}+\mathbf{r},t)d\mathbf{r},\hskip
21.68231pt\mathcal{F}\left\\{\overline{u}_{i}^{\ell}\right\\}=\left(2\pi\right)^{3}\mathcal{F}\left\\{G^{\ell}\right\\}\mathcal{F}\left\\{u_{i}\right\\},$
(14)

The relative width of the filter is $\ell$, with larger values of $\ell$
resulting in coarser fields. The notation $\mathcal{F}$ is used for a Fourier
transform in three dimensions. The definition used here for an arbitrary field
$a$ is,

$\mathcal{F}\left\\{a\right\\}(\mathbf{k})=\frac{1}{\left(2\pi\right)^{3}}\iiint
a(\mathbf{x})\exp\left(-i\mathbf{k}\cdot\mathbf{x}\right)d\mathbf{x}.$ (15)

The three-dimensional inverse Fourier transform for an arbitrary function $b$
is therefore,

$\mathcal{F}^{-1}\left\\{b\right\\}(\mathbf{x})=\iiint
b(\mathbf{k})\exp\left(i\mathbf{k}\cdot\mathbf{x}\right)d\mathbf{k}.$ (16)

This Fourier transform definition is commonly used for turbulence theory,
e.g., (Tennekes & Lumley, 1972; Monin & Yaglom, 1975; Pope, 2000).

LES practitioners rarely rely on a specific filter kernel shape. Common filter
shapes used for theoretical appraisals include top-hat, spectral cutoff, and
Gaussian kernels (Borue & Orszag, 1998). For example, the Gaussian filter
kernel is given by

$G^{\ell}(\mathbf{r})=\frac{1}{(2\pi\ell)^{3/2}}\exp\left(-\frac{|\mathbf{r}|^{2}}{2\ell^{2}}\right),~{}~{}~{}~{}~{}~{}~{}\left(2\pi\right)^{3}\mathcal{F}\\{G^{\ell}\\}(\mathbf{k})=\exp\left(-\frac{1}{2}|\mathbf{k}|^{2}\ell^{2}\right),$
(17)

The Gaussian filter provides a balanced trade-off between localization in
physical space and wavenumber space, with exponential-of-square dropoff in
both.

For a more general class of filter shapes, the Navier-Stokes equation, Eq.
(1), may be filtered to yield an equation for the filtered velocity field,

$\frac{\partial\overline{u}_{i}^{\ell}}{\partial t}+\frac{\partial}{\partial
x_{j}}\left[\overline{u}_{i}^{\ell}\overline{u}_{j}^{\ell}+\tau_{ij}^{\ell}+\overline{p}^{\ell}\delta_{ij}-\nu\left(\frac{\partial\overline{u}_{i}^{\ell}}{\partial
x_{j}}+\frac{\partial\overline{u}_{j}^{\ell}}{\partial
x_{i}}\right)\right]=0,\hskip
43.36464pt\frac{\partial\overline{u}_{j}^{\ell}}{\partial x_{j}}=0.$ (18)

Here, the subfilter stress tensor,
$\tau_{ij}^{\ell}=\overline{u_{i}u_{j}}^{\ell}-\overline{u}_{i}^{\ell}\overline{u}_{j}^{\ell}$,
plays a role analogous to that of the residual stress tensor in PIC,
$\sigma_{ij}$, though their mathematical definitions may differ. Also, the
filtered pressure, $\overline{p}^{\ell}$, is analogous to the generalized
pressure in PIC, $\phi$.

The kinetic energy of the filtered velocity field evolves as,

$\frac{\partial(\tfrac{1}{2}\overline{u}_{i}^{\ell}\overline{u}_{i}^{\ell})}{\partial
t}+\frac{\partial}{\partial
x_{j}}\left[\frac{1}{2}\overline{u}_{i}^{\ell}\overline{u}_{i}^{\ell}\overline{u}_{j}^{\ell}+\overline{u}_{i}\tau_{ij}^{\ell}+\overline{p}^{\ell}\overline{u}_{j}^{\ell}-2\nu\overline{u}_{i}^{\ell}\overline{s}_{ij}^{\ell}\right]=\tau_{ij}^{\ell}\overline{s}_{ij}^{\ell}-2\nu\overline{s}_{ij}^{\ell}\overline{s}_{ij}^{\ell}.$
(19)

The strain-rate tensor is $s_{ij}=\tfrac{1}{2}\left(\partial u_{i}/\partial
x_{j}+\partial u_{j}/\partial x_{i}\right)$ and the filtered strain-rate
tensor is
$\overline{s}_{ij}^{\ell}=\tfrac{1}{2}(\partial\overline{u}_{i}^{\ell}/\partial
x_{j}+\partial\overline{u}_{j}^{\ell}/\partial x_{i})$. The two sinks of the
filtered kinetic energy are analogous to those of the generalized kinetic
energy in Eq. (13).

Note that the preceding review of spatial filtering theory has implicitly
assumed a spatially uniform filtering operation, i.e.,
$\ell\neq\ell(\mathbf{x})$. Specifically, the derivation of Eq. (18) by
applying a filter operation to Eq. (1) relies on the commutative property of
filtering with spatial differentiation, which does not hold if the filter
width varies in space. In the case of spatially-varying filter width,
additional commutator errors arise which (theoretically) require additional
treatment (Ghosal & Moin, 1995; Langford & Moser, 2001; Yalla et al., 2021;
Moser et al., 2021). The above derivation also presumed an unbounded domain in
the filter definition. Section 3 explores the similarity between PIC and
filtering in the context of unbounded flows with uniform resolution. Even in
this case, PIC provides useful insight. Then, section 4 shows how PIC may be
extended to more complex flows, including those which represent a challenge to
the spatial filtering approach to LES theory.

##  3 Unbounded flows with uniform resolution

###  3.1 Equivalence of PIC with Gaussian filtering

For unbounded flows with uniform resolution (uniform $\hat{\nu}$), the pseudo-
pressure satisfies a Laplace equation in free space, so $\hat{p}=0$. That is,
the pseudo-pressure is not needed to enforce the divergence-free condition on
the generalized velocity. In this case, the auxiliary Stokes equation
simplifies to,

$\frac{\partial w_{i}}{\partial\hat{t}}=\hat{\nu}\nabla^{2}w_{i},\hskip
43.36464ptw_{i}(\mathbf{x},t;0)=u_{i}^{\ell}(\mathbf{x},t).$ (20)

The formal solution to Eq. (20) in an unbounded domain is readily obtained

$w_{i}(\mathbf{x},t;\hat{t})=\iiint\frac{1}{\left(4\pi\hat{\nu}\hat{t}\right)^{3/2}}\exp\left(-\frac{|\mathbf{r}|^{2}}{4\hat{\nu}\hat{t}}\right)u_{i}(\mathbf{x}+\mathbf{r},t)d\mathbf{r}\equiv\overline{u}_{i}^{\ell}(\mathbf{x},t),$
(21)

which is precisely a Gaussian filter with width
$\ell=\sqrt{2\hat{\nu}\hat{t}}$. Therefore, the physics-inspired coarsening
approach is precisely equivalent to spatial filtering with a Gaussian kernel
for unbounded flows with uniform resolution. All the strengths and theoretical
insights of spatial filtering naturally carry over, but a new perspective on
the residual stress tensor also emerges.

###  3.2 PIC expression for the residual stress tensor

With the equivalence to Gaussian filtering established for PIC of unbounded
flows with uniform resolution, the mathematical definition of the residual
stress tensor could be straightfowardly written as
$\sigma_{ij}=\tau_{ij}\equiv\overline{u_{i}u_{j}}^{\ell}-\overline{u}_{i}^{\ell}\overline{u}_{j}^{\ell}$.
However, insights beyond those arrived at via spatial filtering may be
obtained by following the PIC logic further in the form of Eq. (11).

In the following, flow subjected to an arbitrary (divergence-free) forcing
function, $\mathbf{f}(\mathbf{x},t)$ is considered, so Eq. (1) becomes,

$\frac{\partial u_{i}}{\partial t}+\frac{\partial}{\partial
x_{j}}\left[u_{i}u_{j}+p\delta_{ij}-\nu\left(\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\right]=f_{i},\hskip
43.36464pt\frac{\partial u_{j}}{\partial x_{j}}=0.$ (22)

The generalized Navier-Stokes equations are likewise supplemented with a
generalized force, $\mathbf{g}(\mathbf{x},t)$, that represents the impact of
the physical forcing on the coarsened flow representation,

$\frac{\partial w_{i}}{\partial t}+\frac{\partial}{\partial
x_{j}}\left[w_{i}w_{j}+\sigma_{ij}+\phi\delta_{ij}-\nu\left(\frac{\partial
w_{i}}{\partial x_{j}}+\frac{\partial w_{j}}{\partial
x_{i}}\right)\right]=g_{i},\hskip 43.36464pt\frac{\partial w_{j}}{\partial
x_{j}}=0.$ (23)

Substituting Eqs. (20) and (23) into Eq. (11),

$\frac{\partial}{\partial\hat{t}}\left(\nu\nabla^{2}w_{i}-\frac{\partial\phi}{\partial
x_{i}}-\frac{\partial w_{i}w_{j}}{\partial
x_{j}}-\frac{\partial\sigma_{ij}}{\partial
x_{j}}+g_{i}\right)-\frac{\partial}{\partial
t}\left(\hat{\nu}\nabla^{2}w_{i}\right)=0$ (24)

and rearranging in terms of time and pseudo-time derivatives,

$\nu\nabla^{2}\left(\frac{\partial
w_{i}}{\partial\hat{t}}\right)-\frac{\partial}{\partial
x_{i}}\left(\frac{\partial\phi}{\partial\hat{t}}\right)-\frac{\partial}{\partial
x_{j}}\left(w_{i}\frac{\partial w_{j}}{\partial\hat{t}}+w_{j}\frac{\partial
w_{i}}{\partial\hat{t}}\right)-\frac{\partial}{\partial
x_{j}}\left(\frac{\partial\sigma_{ij}}{\partial\hat{t}}\right)+\frac{\partial
g_{i}}{\partial\hat{t}}=\hat{\nu}\nabla^{2}\left(\frac{\partial
w_{i}}{\partial t}\right)$ (25)

and further substitution of Eqs. (20) and (23),

$\nu\nabla^{2}\left(\hat{\nu}\nabla^{2}w_{i}\right)-\frac{\partial}{\partial
x_{i}}\left(\frac{\partial\phi}{\partial\hat{t}}\right)-\frac{\partial}{\partial
x_{j}}\left(\hat{\nu}w_{i}\nabla^{2}w_{j}+\hat{\nu}w_{j}\nabla^{2}w_{i}\right)-\frac{\partial}{\partial
x_{j}}\left(\frac{\partial\sigma_{ij}}{\partial\hat{t}}\right)+\frac{\partial
g_{i}}{\partial\hat{t}}\\\
=\hat{\nu}\nabla^{2}\left(\nu\nabla^{2}w_{i}-\frac{\partial\phi}{\partial
x_{i}}-\frac{\partial w_{i}w_{j}}{\partial
x_{j}}-\frac{\partial\sigma_{ij}}{\partial x_{j}}+g_{i}\right).$ (26)

Now, further simplification leads to,

$-\frac{\partial}{\partial
x_{j}}\left(\frac{\partial\sigma_{ij}}{\partial\hat{t}}-\hat{\nu}\nabla^{2}\sigma_{ij}-2\hat{\nu}\frac{\partial
w_{i}}{\partial x_{k}}\frac{\partial w_{j}}{\partial
x_{k}}+\frac{\partial\phi}{\partial\hat{t}}\delta_{ij}-\hat{\nu}\nabla^{2}\phi\delta_{ij}\right)=\frac{\partial
g_{i}}{\partial\hat{t}}-\hat{\nu}\nabla^{2}g_{i}.$ (27)

Equation (27) simply expresses the condition necessary for Eq. (23) to
correctly describe the large-scale dynamics embodied in the generalized
velocity at finite pseudo-time. Because the residual stress should be
disentangled from the (arbitrary) forcing function, both sides for Eq. (27)
should be set to zero. Setting the right-hand side to zero, the generalized
forcing function may be recognized as the Gaussian-filtered force, c.f. Eq.
(20),

$\frac{\partial g_{i}}{\partial\hat{t}}=\hat{\nu}\nabla^{2}g_{i},\hskip
43.36464ptg_{i}(\mathbf{x},t;0)=f_{i}(\mathbf{x},t).$ (28)

The left-hand side of Eq. (27), once set to zero, may be simplified. Namely,
the terms involving the generalized pressure, $\phi$, may be removed on either
of two considerations. First, the generalized pressure may be identified with
the Gaussian-filtered pressure,

$\frac{\partial\phi}{\partial\hat{t}}=\hat{\nu}\nabla^{2}\phi,\hskip
43.36464pt\phi(\mathbf{x},t;0)=p(\mathbf{x},t).$ (29)

Alternatively, and more generally, the dynamics of the generalized velocity
depend only on the deviatoric part of the residual stress tensor, so any
isotropic contribution can be safely ignored.

A sufficient condition for satisfying Eq. (27) and thus Eq. (11) is

$\frac{\partial\sigma_{ij}}{\partial\hat{t}}=\hat{\nu}\nabla^{2}\sigma_{ij}+2\hat{\nu}\frac{\partial
w_{i}}{\partial x_{k}}\frac{\partial w_{j}}{\partial x_{k}}.$ (30)

This is a forced diffusion equation in pseudo-time for $\boldsymbol{\sigma}$,
and its formal solution in unbounded space is also readily obtained

$\sigma_{ij}(\mathbf{x},t;\hat{t})=\int_{0}^{\hat{t}}\left[\iiint\frac{1}{\left(4\pi\hat{\nu}(\hat{t}-t^{\prime})\right)^{3/2}}\exp\left(-\frac{|\mathbf{r}|^{2}}{4\hat{\nu}(\hat{t}-t^{\prime})}\right)\frac{\partial
w_{i}}{\partial x_{k}}\frac{\partial w_{j}}{\partial
x_{k}}(\mathbf{x}+\mathbf{r},t;t^{\prime})d\mathbf{r}\right]dt^{\prime}.$ (31)

The residual stress at pseudo-time $\hat{t}$ may thus be interpreted as an
integral over all earlier pseudo-times, $0\leq t^{\prime}\leq t$, of the
velocity gradient product at $t^{\prime}$ smoothed by a pseudo-viscosity from
its earlier pseudo-time up until the pseudo-time at which the residual stress
is evaluated, i.e., over the pseudo-time range $\hat{t}-t^{\prime}$.

###  3.3 Insights into the energy cascade

The equivalence of PIC with Gaussian spatial filtering may be invoked, in the
case of unbounded flows with uniform resolution, to rephrase the subfilter
stress tensor of filtering theory in terms of multiscale velocity gradients.
Dividing Eq. (30) by $2\hat{\nu}$, with the established relation
$\ell^{2}=2\hat{\nu}\hat{t}$,

$\frac{\partial\tau_{ij}}{\partial(\ell^{2})}=\frac{1}{2}\nabla^{2}\tau_{ij}+\frac{\partial\overline{u}_{i}^{\ell}}{\partial
x_{k}}\frac{\partial\overline{u}_{j}^{\ell}}{\partial x_{k}}$ (32)

therefore has a formal solution, Eq. (31), that can itself be written in terms
of an integral of Gaussian filtered velocity gradients at all scales
$0\leq\sqrt{\alpha}\leq\ell$,

$\tau_{ij}=\int_{0}^{\ell^{2}}\overline{\frac{\partial\overline{u}_{i}^{\sqrt{\alpha}}}{\partial
x_{k}}\frac{\partial\overline{u}_{j}^{\sqrt{\alpha}}}{\partial
x_{k}}}^{\beta}d\alpha$ (33)

where $\beta=\sqrt{\ell^{2}-\alpha}$ is the width of the complementary filter
which smooths from scale $\sqrt{\alpha}$ to scale $\ell$.

Invoking the definition of the generalized second moment from Germano (1992),

$\tau^{\ell}(a,b)=\overline{ab}^{\ell}-\overline{a}^{\ell}\overline{b}^{\ell},$
(34)

of which the subfilter stress tensor is one special case
$\tau_{ij}=\tau(u_{i},u_{j})$, the PIC-based phrasing of the Gaussian filter’s
stress tensor can be split into resolved scale and subfilter scale
contributions,

$\tau_{ij}=\ell^{2}\frac{\partial\overline{u}_{i}^{\ell}}{\partial
x_{k}}\frac{\partial\overline{u}_{j}^{\ell}}{\partial
x_{k}}+\int_{0}^{\ell^{2}}\tau^{\beta}\left(\frac{\partial\overline{u}_{i}^{\sqrt{\alpha}}}{\partial
x_{k}},\frac{\partial\overline{u}_{j}^{\sqrt{\alpha}}}{\partial
x_{k}}\right)d\alpha.$ (35)

The first of the two terms on the right hand side is the nonlinear gradient
model, which by itself performs well in a priori testing compared with eddy
viscosity models (Clark et al., 1979; Borue & Orszag, 1998). The second term
shows how such a model leaves out smaller-scale (unresolved) content. Further
decomposing the filtered velocity gradients in Eq. (35) into strain-rate and
vorticity, and forming the product
$\Pi^{\ell}=-\tau_{ij}^{\ell}\overline{s}_{ij}^{\ell}$, the local energy
cascade rate may be written in terms of vorticity stretching and strain self-
amplification,

$\Pi^{\ell}=\Pi_{s1}^{\ell}+\Pi_{\omega 1}^{\ell}+\Pi_{s2}^{\ell}+\Pi_{\omega
2}^{\ell}+\Pi_{c}^{\ell},$ (36)

where

$\Pi_{s1}^{\ell}=-\ell^{2}\overline{s}_{ij}^{\ell}\overline{s}_{jk}^{\ell}\overline{s}_{ki}^{\ell}=\text{strain-
rate self-amplification at scale~{}}\ell,$ (37) $\Pi_{\omega
1}^{\ell}=\frac{1}{4}\ell^{2}\overline{s}_{ij}^{\ell}\overline{\omega}_{i}^{\ell}\overline{\omega}_{j}^{\ell}=\text{vorticity
stretching at scale~{}}\ell,$ (38)
$\Pi_{s2}^{\ell}=-\overline{s}_{ij}^{\ell}\int_{0}^{\ell^{2}}\tau_{\beta}\left(\overline{s}_{jk}^{\sqrt{\alpha}},\overline{s}_{ki}^{\sqrt{\alpha}}\right)~{}d\alpha=\text{multiscale
strain amplification by strain at ~{}}\ell,$ (39) $\Pi_{\omega
2}^{\ell}=\overline{s}_{ij}^{\ell}\int_{0}^{\ell^{2}}\tau_{\beta}\left(\overline{\omega}_{i}^{\sqrt{\alpha}},\overline{\omega}_{j}^{\sqrt{\alpha}}\right)~{}d\alpha=\text{multiscale
vorticity stretching by strain at~{}}\ell,$ (40)
$\Pi_{c}^{\ell}=2\overline{s}_{ij}^{\ell}\int_{0}^{\ell^{2}}\tau_{\beta}\left(\overline{s}_{jk}^{\sqrt{\alpha}},\overline{\Omega}_{ki}^{\sqrt{\alpha}}\right)~{}d\alpha=\text{multiscale
vortex thinning by strain at~{}}\ell.$ (41)

In this way, the rephrasing of spatial filtering in terms of PIC for unbounded
flows with uniform resolution leads to a unique theoretical insight: the exact
relation of commonly invoked mechanisms to the energy cascade. For more
details, the reader is referred to Johnson (2020, 2021a). A treatment of the
topic for a broader audience is given in Johnson (2021b).

Direct numerical simulations of homogeneous isotropic turbulence demonstrate
that the multiscale gradient terms in $\boldsymbol{\tau}$ align closely with
the filtered strain-rate tensor (Johnson, 2021a). The eddy viscosity model is
thus a good approximation for the second term on the right hand side of Eq.
(35). Overall, this suggests the applicability of a mixed model, with
$\boldsymbol{\sigma}=\boldsymbol{\tau}$ being the sum of nonlinear gradient
and eddy viscosity terms (Vreman et al., 1994).

###  3.4 An alternative to the Germano identity

The dynamic procedure (Germano et al., 1991) based on the Germano identity
(Germano, 1992) is one of the most common direct uses of spatial filtering
theory in LES practice. The basic idea is that coefficients for a given model
may be estimated from the resolved scales using a test filter larger than the
grid size. A similar procedure may be developed with physics-inspired
coarsening via the dependence of $\boldsymbol{\sigma}$ on the cutoff scale
(i.e., pseudo-time).

This may be briefly demonstrated for a generic eddy viscosity model,

$\sigma_{ij}^{(d)}\approx-2\nu_{T}S_{ij}.$ (42)

First, Eq. (42) is substituted into Eq. (30) and the spatial variability of
$\nu_{T}$ is neglected (for convenience, as commonly done for the Smagorinsky
coefficient),

$-2S_{ij}\frac{\partial\nu_{T}}{\partial\hat{t}}=2\hat{\nu}\left(\frac{\partial
w_{i}}{\partial x_{k}}\frac{\partial w_{j}}{\partial
x_{k}}-\frac{1}{3}\frac{\partial w_{m}}{\partial x_{n}}\frac{\partial
w_{m}}{\partial x_{n}}\delta_{ij}\right).$ (43)

The above expression may be modified to the Smagorinski form and the spatial
variability of the filtered strain-rate may also be factored in, if desired.

Second, Kolmogorov (1941) scaling is assumed,
$\nu_{T}\sim\ell^{4/3}\sim\hat{t}^{2/3}$, so that the variation of the eddy
viscosity in pseudo-time may be evaluated,
$\frac{\partial\nu_{T}}{\partial\hat{t}}=\frac{2}{3}\frac{\nu_{T}}{\hat{t}}$.
Then Eq. (45) becomes,

$-\frac{4}{3}S_{ij}\nu_{T}=\ell^{2}\left(\frac{\partial w_{i}}{\partial
x_{k}}\frac{\partial w_{j}}{\partial x_{k}}-\frac{1}{3}\frac{\partial
w_{m}}{\partial x_{n}}\frac{\partial w_{m}}{\partial x_{n}}\delta_{ij}\right)$
(44)

where $\ell^{2}=2\hat{\nu}\hat{t}$. This equation is over-determined because,
in general, a scalar $\nu_{T}$ will not be found that can satisfy the full
tensor equation. Therefore, as is typical for dynamic procedures, a least-
squares procedure is employed (Lilly, 1992). This step leads to a projection
on $S_{ij}$ that carries the physical meaning of matching the energy cascade
rate. The resulting expression for the eddy viscosity is,

$\nu_{T}=\frac{3\ell^{2}}{4}\frac{A_{ik}A_{jk}A_{ij}}{S_{mn}S_{mn}}=\frac{3\ell^{2}}{4}\frac{\frac{1}{4}W_{i}S_{ij}W_{j}-S_{ij}S_{jk}S_{ki}}{S_{mn}S_{mn}},$
(45)

where $A_{ij}=\partial w_{i}/\partial x_{j}$ is the generalized velocity
gradient tensor, $S_{ij}=\tfrac{1}{2}(A_{ij}+A_{ji})$ is the generalized
strain-rate tensor, and $W_{i}=\epsilon_{ijk}A_{kj}$ is the generalized
vorticity. The denominator of this relation is positive-definite. The
numerator represents the sum of vortex stretching and strain-rate self-
amplification, which is positive on average for turbulent flows. Local
negative values of the numerator may necessitate averaging or clipping
strategies as done in the Germano-based dynamic procedure.

Johnson (2021a) demonstrated that the eddy viscosity approximation is more
physically accurate for the only part of the residual stress tensor,

$\epsilon_{ij}=\sigma_{ij}-2\hat{\nu}\hat{t}\frac{\partial w_{i}}{\partial
x_{k}}\frac{\partial w_{j}}{\partial x_{k}},$ (46)

which represents multiscale vortex stretching and strain-rate amplification in
Eq. (35). It may be shown from Eqs. (20) and (30) that,

$\frac{\partial\epsilon_{ij}}{\partial\hat{t}}=\hat{\nu}\nabla^{2}\epsilon_{ij}+4\hat{\nu}^{2}\hat{t}\frac{\partial^{2}w_{i}}{\partial
x_{m}\partial x_{n}}\frac{\partial^{2}w_{j}}{\partial x_{m}\partial x_{n}}$
(47)

Substituting the eddy viscosity approximation, Eq. (42), for
$\epsilon_{ij}^{(d)}$ rather than the full residual stress tensor,
$\sigma_{ij}^{(d)}$, and again assuming $\nu_{T}\sim\hat{t}^{2/3}$ with a
least squares approach leads to a ‘dynamic’ mixed model,

$\sigma_{ij}^{(d)}=\ell^{2}\left(A_{ik}A_{jk}\right)^{(d)}-2\nu_{T}S_{ij},\hskip
43.36464pt\nu_{T}=-\frac{3\ell^{4}}{4}\frac{B_{ikl}S_{ij}B_{jkl}}{S_{mn}S_{mn}},$
(48)

where $B_{ijk}=\partial^{2}w_{i}/\partial x_{j}\partial x_{k}$ and
$\ell^{2}=2\hat{\nu}\hat{t}$. The appearance of the second derivative of the
generalized velocity field in the mixed model is less practical when numerical
errors are inherent in assessing spatial derivatives. Nonetheless, Eq. (48)
provides a useful initial test for assessing the relative strengths of PIC
theory.

###  3.5 A priori testing

Direct numerical simulations of homogeneous isotropic turbulence are used for
a priori tests reported in this subsection. The incompressible Navier-Stokes
equation, Eq. (1), is solved using a pseudo-spectral method with $1024$
collocation points in each direction. The fully-resolved velocity field,
$\mathbf{u}$, is advanced in time with a second-order Adams-Bashforth scheme,
and the pressure $p$ simply enforces the divergence-free condition. The
$2\sqrt{2}/3$ rule for wave number truncation is used with phase-shift
dealiasing (Patterson & Orszag, 1971). The forcing, $\mathbf{f}$, is
specifically designed to maintain constant kinetic energy in the first two
wave number shells.

The simulation was initialized using a Gaussian random velocity field
satisfying a model turbulent energy spectrum (Pope, 2000). The simulation was
first run through a startup period to reach statistical stationarity. Then,
statistics are computed over $6$ large-eddy turnover times. The Taylor-scale
Reynolds number is approximately $Re_{\lambda}=400$ with grid resolution
$k_{\max}\eta=1.4$. The integral length scale is about $20\%$ of the periodic
box size of $2\pi$ and $L/\eta=460$. The skewness of the longitudinal velocity
gradient is $-0.58$, and the flatness of the longitudinal and transverse
velocity gradients are $8.0$ and $12.4$, respectively, in reasonable agreement
with previous simulations (Ishihara et al., 2007). The a priori assessment is
carried out for the ‘dynamic’ viscosity model, Eq. (45) and ‘dynamic’ mixed
model, Eq. (48), along with existing popular models including the standard
dynamic Smagorinsky model (Germano et al., 1991; Lilly, 1992) and nonlinear
gradient model (Clark et al., 1979; Borue & Orszag, 1998). For all three
dynamic models, spatial averaging is used for the numerator and denominator of
the Smagorinsky coefficient or eddy viscosity. The results are shown in Figure
2. For this case, PIC and Gaussian filtering are equivalent so the language of
each approach is used interchangeably.

(a) (b)
(c) (d)
(e) (f)

Figure 2: A priori comparison of new dynamic eddy viscosity (‘dyn visc’) and
dynamic mixed (‘dyn mixed’) models with DNS results and existing models:
dynamic Smagorinsky (‘dyn Smag’) and nonlinear gradient model (‘nonlin grad’):
(a) Mean cascade rate, (b) PDF of energy cascade rate at $\ell=46\eta$, (c-d)
correlation coefficient and R2 coefficient for the energy cascade rate,
respectively, (e-f) correlation coefficient and R2 coefficient for the full
residual stress tensor, respectively. Two vertical dashed gray lines indicate
rough bounds for the approximate inertial subrange of scales in the DNS,
$25\leq\ell/\eta\leq 75$.

Perhaps the most important practical consideration for LES models is providing
the proper rate of energy removal to unresolved scales. Figure 2a shows the
the average cascade rate, $\Pi$, as a function of filter width,
$\ell=\sqrt{2\hat{\nu}\hat{t}}$, normalized by the average dissipation rate in
the unfiltered DNS. The blue curve with star symbols shows that the cascade
rate is nearly equal to the dissipation rate for $25\lesssim\ell/\eta\lesssim
75$, an approximate inertial range for the simulation. The cascade rate decays
to zero as $\ell$ decreases because viscous dissipation becomes significant
near the Kolmogorov scale. The other four curves indicate the extent to which
the four models reproduce the correct cascade rate when applied directly to
coarsened DNS data. As is well-known, the nonlinear gradient model is under-
dissipative, especially in the inertial range, and the dynamic Smagorinsky
model provides a better approximation of the energy removal rate. The PIC-
based dynamic eddy viscosity model performs very similarly to the filter-based
dynamic Smagorinsky model. The PIC-based dynamic mixed model is confirmed to
remove energy at a higher rate in the inertial range than the nonlinear
gradient or eddy viscosity models on their own.

The probability density function (PDF) of the energy cascade rate, $\Pi$, for
each model at $\ell=46\eta$ is shown in Figure 2b along with the coarsened DNS
result. The DNS data indicate the existence of some backscatter, but the
distribution is strongly skewed toward positive cascade rates. The two eddy
viscosity models do not allow any backscatter by construction due to the
spatial averaging used as part of the dynamic procedure. The nonlinear
gradient and PIC-based dynamic mixed model provide relatively appropriate
probabilities of backscatter, with the mixed model having the best overall
agreement with DNS.

Figures 2c and 2d show the correlation coefficient and R-squared value based
on each models point-by-point prediction of the cascade rate, $\Pi$. An
R-squared value of one indicates a perfect model. All four models are highly
correlated with DNS in the inertial range, but the absolute error measured by
the R-squared value shows more variation in the models performances. Both
metrics show that the PIC-based dynamic mixed model is in closest agreement
with DNS. Figures 2e and 2f show the correlation coefficients and R-squared
values based on the full residual (subfilter) stress tensor. The stress tensor
is known to align poorly with the strain-rate tensor, so the eddy viscosity
models perform poorly for both metrics. The nonlinear gradient model performs
quite well and the dynamic mixed model shows only slight improvement. The
major benefit of the mixed model over the nonlinear gradient model has already
been demonstrated in Figure 2a.

###  3.6 A posteriori testing

Large eddy simulations were performed using three stress models: Smagorinsky,
PIC-based dynamic eddy viscosity (Eq. (45)), and PIC-based dynamic mixed (Eq.
(48)). The same pseudo-spectral code from the DNS is used. The Smagorinsky
coefficient is not determined dynamically for the LES, but was manually chosen
to produce accurate results. For the two PIC-based dynamic models, spatial
averaging is used for the numerator and denominator of the eddy viscosity, so
that the eddy viscosity does not vary in space for either. As such, the
simulation with the PIC-based dynamic eddy viscosity model is essentially a
DNS at a lower Reynolds number, with the PIC-based theory setting the
viscosity so as to match Gaussian-filtered DNS with a given filter width.

Figure 3: Energy spectra for a posteriori testing: (a) using PIC-based dynamic
eddy viscosity model with $\ell/\eta=24$ at three different resolutions, and
(b) using three different models with $k_{\text{max}}\ell=3$. A Kolmogorov
spectrum, $E(k)=1.6\epsilon^{2/3}k^{-5/3}$ is shown for reference in both
panels.

The filter width (pseudo-time) is chosen to be $\ell/\eta=24$, which is at the
small-scale end of the inertial range identified in the a priori tests. The
following grid resolutions were tested: $64^{3}$, $128^{3}$, and $256^{3}$.
These correspond to $k_{\text{max}}\ell=1.5$, $3.0$, and $6.0$. The DNS
results are filtered with a Gaussian kernel (equivalent to PIC) for direct
comparison with the LES results from each model. Figure 3a shows the energy
spectra from LES simulations at each of the three resultions with the PIC-
based dynamic eddy viscosity model. The spectra from each resolution overlap
until wavenumbers close to $k_{\text{max}}$, where a slight pileup occurs. The
results from $128^{3}$ and $256^{3}$ are very similar, indicating grid
convergence, but $64^{3}$ is not grid-converged. The spectra from the three
models are compared at $128^{3}$ resolution in Figure 3b. The Smagorinsky
model is known to produce a spectrum close to that of a Gaussian filter Pope
(2000), as affirmed by the results here. The PIC-based dynamic mixed model
produces a spectrum almost indistinguishable from the Gaussian filtered DNS
for the simulation shown here. All three models produce spectra in reasonable
agreement. The PIC-based dynamic procedure for determining the eddy viscosity
is thus shown to work well for both the pure eddy viscosity and mixed models.

In addition to testing the energy spectra produced by various LES models, it
is also useful to compare local flow topology statistics. Given the importance
of strain self-amplification and vorticity stretching to turbulence dynamics,
the LES models are compared to filtered DNS in terms of the following two
quantities,

$s^{*}=\frac{-\sqrt{6}S_{ij}S_{jk}S_{ki}}{\left(S_{mn}S_{mn}\right)^{3/2}}\hskip
43.36464pt\omega^{*}=\frac{\sqrt{6}W_{i}S_{ij}W_{j}}{2W_{k}W_{k}\left(S_{mn}S_{mn}\right)^{1/2}}.$
(49)

The first, $s^{*}$, was introduced by Lund & Rogers (1994) and quantifies the
efficiency of the strain self-amplification contribution to the energy
cascade, $-1\leq s^{*}=\Pi_{s1}/\Pi_{s1,\text{max}}\leq 1$. The second term
likewise indicates the efficiency of the vorticity stretching contribution to
the cascade, $-1\leq\omega^{*}=\Pi_{\omega 1}/\Pi_{\omega 1,\text{max}}\leq
1$. The concept of cascade efficiency was introduced by Ballouz & Ouellette
(2018) and extended to include the above definitions in Johnson (2021a). These
two quantities describe flow topology in a way most relevant to turbulent
cascade physics.

Figure 4: Probability density functions (PDFs) of $s^{*}$ and $\omega^{*}$,
Eq. (49), for LES with three different models compared to DNS and filtered
DNS.

Figure 4 compares the PDF of $s^{*}$ and $\omega^{*}$ from LES models and
filtered DNS. Results from unfiltered DNS are also shown to highlight the
similarity of the eddy viscosity models, particularly the PIC-based dynamic
viscosity model, with (unfiltered) DNS. The reason for this is that the eddy
viscosity models behave like DNS at lower Reynolds numbers, and thus produce
topology statistics similar to the Kolmogorov scales in DNS. However, the
results in Figure 4 demonstrate that the PIC-based mixed model captures subtle
physical differences between viscous-scale and inertial range flow topology
statistics.

Many other metrics can be used to judge the a posteriori accuracy of LES
models. It is not the goal at present to thoroughly explore all possibilities.
However, a number of highly relevant statistics are summarized in Table 1. It
may be noted that the PIC-based dynamic mixed model slightly over-predicts the
magnitude of velocity gradients, which may also be noticed in the spectra
shown in Figure 3b, where a slight over-prediction is noticeable for a range
of wavenumbers slightly smaller than the inverse of the filter scale. There
appears to be a slight (unphysical) build up of energy not entirely unlike the
bottleneck effect commonly observed in DNS and experimental data near the
Kolmogorov scale (Falkovich, 1994; Donzis & Sreenivasan, 2010).

All three models over-predict the fraction of the cascade rate due to single-
scale strain self-amplification and vorticity stretching. For the mixed model,
this shows that too much energy removal is done through the nonlinear gradient
term, and not enough through the eddy viscosity. The longitudinal and
transverse velocity gradient skewness and flatness values from the mixed model
LES are in fairly good agreement with the filtered DNS. Notice that the
dynamic eddy viscosity model leads to too much intermittency (higher skewness
and flatness), because it has a spatially uniform eddy viscosity, whereas the
local adjustments due to Smagorinsky’s eddy viscosity bring intermittency
metrics more in line with filtered DNS.

Overall, all three models perform reasonably well in a posteriori tests.
Furthermore, the relative advantages of the mixed model are clear, and the
success of the PIC-based dynamic scheme is demonstrated. It should be
appreciated that PIC theory for LES can facilitate model development well
beyond the specific models tested here. The success of these models, however,
does motivate future work to develop refined PIC-based models.

model | $\langle S_{ij}S_{ij}\rangle\tau_{\eta}^{2}$ | $\langle\Pi_{s1}\rangle/\langle\Pi\rangle$ | $\langle\Pi_{\omega 1}\rangle/\langle\Pi\rangle$ | $\langle A_{11}^{3}\rangle/\langle A_{11}^{2}\rangle^{3/2}$ | $\langle A_{11}^{4}\rangle/\langle A_{11}^{2}\rangle^{2}$ | $\langle A_{12}^{4}\rangle/\langle A_{12}^{2}\rangle^{2}$
---|---|---|---|---|---|---
fDNS |  $1.7$e-$2$ | $0.37$ | $0.12$ | $-0.41$ | $3.9$ | $4.6$
Smag. |  $1.8$e-$2$ | $0.42$ | $0.14$ | $-0.42$ | $3.6$ | $4.9$
dyn. visc. |  $1.9$e-$2$ | $0.50$ | $0.17$ | $-0.47$ | $4.7$ | $6.7$
dyn. mixed |  $1.9$e-$2$ | $0.45$ | $0.15$ | $-0.41$ | $4.0$ | $4.9$
Table 1: Statistical results from LES and filtered DNS at $\ell/\eta=24$ with
$128^{3}$ resolution ($k_{\text{max}}\ell=3$). These values were converged
when compared with $256^{3}$ resolution.

##  4 PIC for more complex flows

In section 3, it was shown that physics-inspired coarsening (PIC) is
mathematically equivalent to spatial filtering with a Gaussian kernel for
unbounded flows with uniform, isotropic resolution. Even so, the PIC-based
approach to LES provided theoretical and modeling insights, including an PIC-
based alternative to the Germano-based dynamic procedure that does not require
test filtering. Preliminary models showed promising a priori and a posteriori
results, but there is more room for future work developing and testing models
even for that simple case.

In this section, the extension of PIC to a number of more complex flow
scenarios is outlined. The focus is on how the basic theory of PIC provides
the flexibility and systematic approach to incorporate features needed for
practical application of LES to complex flows.

###  4.1 Anisotropic resolution

The first step toward LES of more complex flows is a consistent treatment of
grid anisotropy. Anisotropic grid resolutions are often unavoidable in
simulations with complex geometries. Various treatments of grid anisotropy
exist for spatial filtering theory (Bardina et al., 1980; Scotti et al., 1993;
Vreman, 2004; Rozema et al., 2015; Haering et al., 2019). In practice, simple
treatments of grid anisotropy are common, such as using a scalar measure for
effective resolution based on a suitable average of resolution in three
different directions on an orthogonal grid: $\Delta_{1}$, $\Delta_{2}$, and
$\Delta_{3}$. Deardorff (1970) suggested the use of
$\Delta=\sqrt[3]{\Delta_{1}\Delta_{2}\Delta_{3}}$ based on the cell volume and
Bardina et al. (1980) used
$\Delta=\sqrt{(\Delta_{1}^{2}+\Delta_{2}^{2}+\Delta_{3}^{2})/3}$.

For viscosity-based smoothing, anisotropic resolution is easily introduced
using an tensorial pseudo-viscosity aligned with the principle directions of
the grid. Thus, Eq. (20) for unbounded flows may be modified to incorporate
anisotropic resolution effects,

$\frac{\partial
w_{i}}{\partial\hat{t}}=\hat{\nu}_{jk}\frac{\partial^{2}w_{i}}{\partial
x_{j}\partial x_{k}},\hskip 43.36464pt\frac{\partial w_{j}}{\partial
x_{j}}=0.$ (50)

For uniform resolution, the pseudo-pressure will satisfy the Laplace equation
and hence be $\hat{p}=0$ in the absence of flow boundaries. In that case, the
result will be anisotropic Gaussian filtering aligned with the grid,

$w_{i}(\mathbf{x},t;\hat{t})=\iiint\frac{1}{\sqrt{(2\pi)^{3}\det\mathbf{C}}}\exp\left(-\frac{1}{2}r_{j}~{}C_{jk}^{-1}r_{k}\right)u_{i}(\mathbf{x}+\mathbf{r},t)d\mathbf{r}.$
(51)

Here, the covariance tensor, $C_{ij}=2\hat{\nu}_{ij}\hat{t}$, is equivalent to
the moment of inertial tensor of Bardina et al. (1980). The isotropic form,
Eq. (20), is recovered when $\hat{\nu}_{ij}=\hat{\nu}\delta_{ij}$. The
eigenvalues of $\mathbf{C}$, $\ell_{1,2,3}^{2}=2\hat{\nu}_{1,2,3}~{}\hat{t}$,
represent the (square of the) resolution in three directions given by their
respective eigenvectors (e.g., the local coordinate frame of an orthogonal
grid). The result of anisotropic PIC (i.e., anisotropic Gaussian filtering) is
illustrated in Figure 5. In this case, the coarsening procedure generates
anisotropy in the generalized velocity field from an initially isotropic
fully-resolved velocity field (Haering et al., 2019).

Figure 5: A coarsened velocity field using anisotropic pseudo-viscosity on DNS
data of homogeneous isotropic turbulence. The resolution length scale
($\sqrt{2\hat{\nu}\hat{t}}$) in the x direction (purple/magenta lines) is one-
fourth of that in the y direction (black/gray lines).

Substituting Eqs. (8) and (50) in Eq. (11), one obtains,

$\frac{\partial\sigma_{ij}}{\partial\hat{t}}=\hat{\nu}_{mn}\frac{\partial^{2}\sigma_{ij}}{\partial
x_{m}\partial x_{n}}+\hat{\nu}_{mn}\left(\frac{\partial w_{i}}{\partial
x_{m}}\frac{\partial w_{j}}{\partial x_{n}}+\frac{\partial w_{i}}{\partial
x_{n}}\frac{\partial w_{j}}{\partial x_{m}}\right).$ (52)

The details of the derivation leading to Eq. (52) follow the same steps as
those shown in section 3 leading to Eq. (30). Indeed, it is readily verified
that Eq. (30) is recovered in the isotropic limit,
$\hat{\nu}_{ij}=\hat{\nu}\delta_{ij}$. The PIC-based dynamic procedure for
anisotropic grids can be based on Eq. (52), which could lead to eddy viscosity
models that share some similarities with the anisotropic minimum dissipation
model (Rozema et al., 2015) and the M43 model (Haering et al., 2019), while a
PIC-based dynamic mixed model may share similarities with that of Vreman
(2004).

###  4.2 Nonuniform resolution

Uniform grid resolution is useful for some simple canonical turbulent flows,
but it is not practical for many naturally-occurring and engineered flows. As
discussed in the introduction, a non-uniform filter does not commute with
spatial differentiation, giving rise to additional terms in the governing
equations for filtered fields (Ghosal & Moin, 1995; Yalla et al., 2021). Most
notably, the filtered velocity field is no longer divergence-free (Langford &
Moser, 2001).

In contrast, the physics-inspired coarsening (PIC) approach to the LES
equations explicitly enforces a divergence-free condition on the generalized
velocity field, $\mathbf{w}$, using a pseudo-pressure, $\hat{p}$, in the
auxiliary evolution equation in pseudo-time, $\hat{t}$. The above theory for
uniform resolution was able to ignore the proposed pseudo-pressure due to the
form of its Poisson equation, Eq. (5), which may be alternatively written in
terms of pseudo-viscosity gradients,

$\nabla^{2}\hat{p}=2\frac{\partial\hat{\nu}}{\partial
x_{i}}\nabla^{2}w_{i}+\frac{\partial^{2}\hat{\nu}}{\partial x_{i}\partial
x_{j}}\left(\frac{\partial w_{i}}{\partial x_{j}}+\frac{\partial
w_{j}}{\partial x_{i}}\right).$ (53)

which is a Laplace equation for uniform pseudo-viscosity (i.e., uniform grid
resolution). Nonuniform resolution is represented as spatial variation of the
pseudo-viscosity, $\hat{\nu}=\hat{\nu}(\mathbf{x})$, which activates a non-
zero pseudo-pressure to enforce the divergence free condition. Note that, with
the nonzero pseudo-pressure, the PIC equations for pseudo-time evolution are
no longer purely parabolic, but include an elliptic nature as well. However,
the square-of-exponential falloff of the spectral content is still roughly
preserved, so that the removal of small-scales is still rather efficient
compared to other elliptic PDE based filters (Germano, 1986a, b; Bull &
Jameson, 2016).

Figure 6: A (divergence-free) coarsened velocity field using sinusoidally-
varying pseudo-viscosity applied to DNS data of homogeneous isotropic
turbulence. The resolution length scale, $\sqrt{2\hat{\nu}\hat{t}}$ varies by
a factor of $4$ between $y=0$ and $y=\pi$.

Figure 6 demonstrates the outcome of PIC with a spatially-varying pseudo-
viscosity for a velocity fied snapshot from DNS of homogeneous isotropic
turbulence with
$\hat{\nu}(y)=\hat{\nu}_{\max}-(\hat{\nu}_{\max}-\hat{\nu}_{\min})\cos(y)$,
where $0\leq y\leq 2\pi$ is the vertical axis as printed on the page. The
rightmost image is of a divergenge-free velocity field with variable
resolution. The effective resolution is isotropic, but varies in space, with
high resolution near the top and bottom ($\hat{\nu}\approx\hat{\nu}_{\min}$)
and low resolution near the center ($\hat{\nu}\approx\hat{\nu}_{\max}$).

The case of nonuniform (isotropic) resolution is useful for reflecting on the
physics-inspired coarsening (PIC) approach in general. The goal of PIC is to
create LES equations with well-defined closure terms that facilitate effective
modeling in ways that improve upon spatial filtering. It is not a central goal
of PIC to remain particularly faithful to physical processes during its
pseudo-time evolution. The physical processes inspiring the pseudo-time
dynamics are useful insomuch as they create helpful representations for LES.
This allows some freedom when designing the details of PIC, for example, one
may choose to simplify the auxiliary Stokes equation, Eq. (4), to,

$\frac{\partial w_{i}}{\partial\hat{t}}=-\frac{\partial\hat{p}}{\partial
x_{i}}+\frac{\partial}{\partial x_{j}}\left[\hat{\nu}\frac{\partial
w_{i}}{\partial x_{j}}\right],\hskip 21.68231pt\frac{\partial w_{j}}{\partial
x_{j}}=0,\hskip 21.68231pt\nabla^{2}\hat{p}=\frac{\partial}{\partial
x_{j}}\left(\frac{\partial\hat{\nu}}{\partial x_{i}}\frac{\partial
w_{i}}{\partial x_{j}}\right).$ (54)

The difference between Eq. (4) and Eq. (54) is only material for nonuniform
pseudo-viscosity. Either may be used for PIC, and it is not immediately
evident which will be most beneficial for LES. More research is needed to
elucidate the relative strengths of either.

Substituting Eqs. (54) and (8) into Eq. (11) and performing simplifications
like those in §3, one may arrive at the auxiliary equation for the residual
stress for general $\hat{\nu}=\hat{\nu}(\mathbf{x})$,

$\frac{\partial\sigma_{ij}}{\partial\hat{t}}=\frac{\partial}{\partial
x_{k}}\left(\hat{\nu}\frac{\partial\sigma_{ij}}{\partial
x_{k}}\right)+2\hat{\nu}\frac{\partial w_{i}}{\partial x_{k}}\frac{\partial
w_{j}}{\partial x_{k}}\\\ -\left(w_{i}\frac{\partial w_{k}}{\partial
x_{j}}+w_{j}\frac{\partial w_{k}}{\partial x_{i}}+\frac{\partial
w_{i}w_{k}}{\partial x_{j}}+\frac{\partial\sigma_{ik}}{\partial
x_{j}}-\frac{\partial\hat{p}}{\partial
x_{j}}\delta_{ik}\right)\frac{\partial\hat{\nu}}{\partial
x_{k}}+\left(w_{i}\delta_{jk}+w_{j}\delta_{ik}\right)\frac{\partial\hat{p}}{\partial
x_{k}}.$ (55)

Deviations from Eq. (30) are proportional to the pseudo-viscosity gradient, as
well as the pseudo-pressure. One way to view the extra terms in Eq. (55) is to
view them as consequences of the lack a commutivity between multiplication by
$\hat{\nu}$ and differentiation. The advantage of this PIC formulation is that
these commutator terms are built into the residual stress, and so provide a
path for extending the use of dynamic procedures and other approaches from
uniform resolution to more general cases. Beyond traditional (theory-based)
modeling approaches, data-driven techniques may stand to benefit even more
from PIC and Eq. (55). While isotropic nonuniform resolution has been
considered here, the same procedure may be used to combine anisotropic and
nonuniform grid effects into a single formulation.

###  4.3 Heat and mass transfer

In addition to complexities introduced by numerical grid effects, physics-
inspired coarsening provides a holistic approach for reducing the
computational DoFs needed to represent turbulence in regimes with additional
physical effects beyond unbounded single-phase incompressible flows. The next
few subsections touch briefly on a few examples, the simplest of which is
passive scalar transport.

For flows with heat or mass transfer described by an advection diffusion
equation,

$\frac{\partial T}{\partial t}+u_{j}\frac{\partial T}{\partial
x_{j}}=\alpha\nabla^{2}T$ (56)

the relevant scalar field $T$ (e.g., temperature or concentration) may be
coarsened in a physics-inspired manner using an pseudo-conductivity or pseudo-
diffusivity, $\hat{\alpha}$, to form a generalized temperature/concentration,
$\theta$, with pseudo-time evolution given by an auxiliary heat equation,

$\frac{\partial\theta}{\partial\hat{t}}=\hat{\alpha}\nabla^{2}\theta,\hskip
21.68231pt\frac{\partial\theta}{\partial
t}+w_{j}\frac{\partial\theta}{\partial
x_{j}}=\alpha\nabla^{2}\theta-\frac{\partial q_{j}}{\partial x_{j}}$ (57)

so that the physical time evolution of the generalized
temperature/concentration must include a residual heat/mass flux $\mathbf{q}$.
The equality of mixed partial derivatives for $\theta$,

$\frac{\partial}{\partial\hat{t}}\left(\frac{\partial\theta}{\partial
t}\right)=\frac{\partial}{\partial
t}\left(\frac{\partial\theta}{\partial\hat{t}}\right)$ (58)

may be satisfied by a forced diffusion equation for the residual heat/mass
flux,

$\frac{\partial
q_{i}}{\partial\hat{t}}=\hat{\alpha}\nabla^{2}q_{i}+2\hat{\alpha}\frac{\partial
w_{i}}{\partial x_{j}}\frac{\partial\theta}{\partial x_{j}}.$ (59)

Here, $\hat{\alpha}=\hat{\nu}$ is chosen (unity pseudo-Prandtl number) for
convenience. Other choices of pseudo-Prandtl number are possible, but unlikely
to be advantageous.

Similar extensions to nonuniform or anisotropic resolution may be found for
heat/mass transfer as shown for the momentum equation above. Equation (59) can
also provide insight into scale-wise scalar fluctuation dynamics. A turbulent
cascade of scalar variance may be linked to the multiscale squeezing of scalar
filaments using the formal solution in wall-free flows following the procedure
outline for the kinetic energy cascade in section 3. Furthermore, Eq. (59) can
serve as the basis for a PIC-based dynamic procedure for an eddy diffusivity
or mixed model for $\mathbf{q}$.

###  4.4 Wall-bounded turbulent flows

Turbulence near solid boundaries is very common in flows of interest for LES.
This is the source of some inconvenience for spatial filtering theory. A
uniform filter integral will require flow information outside the fluid domain
(Drivas & Nguyen, 2018). This may be avoided through the use of a filter width
which vanishes at the boundary is approached (Bose & Moin, 2014), giving rise
to commutation errors discussed in §4.2.

In contrast to spatial filtering, the physics-inspired approach naturally
extends to wall-bounded flows. To accomplish this, the auxiliary Stokes
equation, Eq. (4) must be given boundary conditions for $\mathbf{w}$. For
example, the use of user-defined boundary conditions for coarsening procedures
has been preliminarily investigated (Bae & Lozano-Duran, 2017, 2018). A number
of choices are possible depending on the level of near-wall resolution. If the
wall-normal grid spacing resolves the viscous sublayer, $\ell_{y}^{+}\sim 1$,
then the best choice of boundary conditions is likely Dirichlet,
$\mathbf{w}|_{\text{wall}}=0$. Depending on the wall parallel resolution, this
approach could facilitate wall-resolved LES ($\ell_{x}^{+}\sim\ell_{z}^{+}\sim
50$) or near-wall RANS for hybrid RANS-LES (coarser x-z resolution).

Alternatively, a coarser wall-normal resolution, $\ell_{y}\gg 1$, would
facilitate development of approximate closures for wall-modeled LES. In this
case, Dirichlet boundary conditions could still be used (Bae & Lozano-Duran,
2021), but other choices such as Robin (Bose & Moin, 2014) or Neumann boundary
conditions could be used. The choice of boundary conditions for the pseudo-
time dynamics would directly set the boundary conditions to be used in LES on
the generalized velocity field. Furthermore, complementary boundary conditions
on the residual stress tensor, $\boldsymbol{\sigma}$, would also need to be
specified in such a way as to preserve the integral of momentum transport
across the boundary during the pseudo-time evolution so as to facilitate
accurate calculations of forces in the LES. For example, in a wall-resolved
LES regime, Dirichlet boundary conditions would likely be the best choice for
the residual stress in conjunction with Dirichlet velocity boundary
conditions. Further research is needed to determine the best approach for
wall-modeled LES.

Data-driven techniques may be useful for testing the potential fruitfulness of
various choices for the boundary conditions. For example, the art of
constructing data-driven wall models may benefit substantially from the PIC-
based approach to the LES equations. Sufficiently far from the wall,
$y\gg\ell$, the physics-inspired coarsening naturally recovers Gaussian
filtering, but near the wall, the PIC theory outlined in this paper provides a
method for determining how the choice of boundary conditions influences the
residual stress tensor.

###  4.5 Particle-laden flows

Turbulent flows with small particles also represent a challenge to spatial
filtering theory. Lagrangian particle tracking methods, in combination with an
Eulerian representation of the carrier fluid, are increasingly common to
computing particle-laden flows. In the physics-inspired coarsening approach,
the flow disturbance in the vicinity of the particle would be smoothed out,
while boundary conditions may be freely chosen on the surface of the particle,
as in the case of wall-bounded flows above. For example, choosing Neumann
boundary conditions would, for particles smaller than the resolution length
scale, $\ell\sim\sqrt{2\hat{\nu}\hat{t}}$, cause the carrier flow velocity on
the surface of the particle to relax to that of the far-field coarsened
velocity. This would be similar to the idea of recovering the ‘undisturbed’
fluid velocity upon which drag law formulations are based (Horwitz & Mani,
2016, 2018; Balachandar et al., 2019) generalized to the case in which the
Kolmogorov scale is not captured by the grid. The considerations for particle-
resolved calculations would be the same as for wall-bounded turbulence the
previous subsection.

###  4.6 Multiphase flows with resolved interfaces

Multiphase turbulent flows include the additional complexity of moving
discontinuities at phase interfaces with the accompanying surface tension
force. Indeed, the interface is “a critical feature of such flows and it is
likely that coarsening of the flow must retain the interface, although
probably with a simplified structure” (Tryggvason & Lu, 2020). Difficulties
arise in spatial filtering across discontinuities (Sagaut & Germano, 2005;
Toutant et al., 2009), and spatially filtering across fluid-fluid interfaces
is not consistent with the typical (and desirable) sharp treatment of solid
boundaries (walls) in LES. Indeed, spatial filtering across interfaces blurs
unnecessarily blurs essential physics at the interface, given the
sophisticated numerical tools that have developed for treating sharp
interfaces.

Various alternatives to spatial filtering for LES of multiphase flows have
been proposed (see Cheng et al. (2019) for a more thorough, up-to-date
review). Herrmann (2013) proposed a dual-scale resolution procedure that
requires fine resolution of the interface embedded within coarser flow
resolution. McCaslin & Desjardins (2014) suggested, but did not pursue, the
idea of a surface filtering operation. Most recently, Chen et al. (2021) has
advocated for a procedure very similar to the physics-inspired approach
outlined here. Indeed, their approach, which specifically targets the use of
data-driven modeling, could be thought of as one possibility within the PIC
framework.

One possibility suggested by the physics-inspired philosophy is the use of an
artificial surface tension along with the auxiliary Stokes equation to remove
small-scale features. Physically speaking, surface tension limits the impact
of turbulent motions on multiphase features (drops, bubbles, etc) smaller than
the Hinze scale (Hinze, 1955). Thus, just as the pseudo-viscosity artificially
enhances the Kolmogorov scale, a pseudo surface tension acting during the
pseudo-time evolution would likewise smooth small scale interface curvature
while transforming more abrupt small-scale features into spherical particles
for Lagrangian particle tracking, see Figure 7. Note that this PIC process
would maintain a sharp interface treatment for features that can be captured
on the coarse LES grid while creating low DoF representations for unresolved
features. Such an approach provides a general framework to justify other
existing heuristic approaches, e.g., for physics-based conversion of subgrid
ligaments to point-particle representations (Kim & Moin, 2020). The use of an
artificial diffusion equation for Gaussian filtering has some precedence for
simulating particle-laden flows (Capecelatro & Desjardins, 2013). An
artificial surface tension may also be a useful approach for treating
multiscale (partially-resolved) wall roughness.

Figure 7: Sketch of how a pseudo surface tension could create low DoF
representations of a complex two-phase interface.

##  5 Conclusions

In turbulence, viscosity provides a natural mechanism preventing motions
smaller than the Kolmogorov scale. In this paper, it is proposed to view an
LES velocity field as the result of (artificial) viscous smoothing rather than
a spatial filtering operating. This approach is called physics-inspired
coarsening (PIC). In PIC, a pseudo-viscosity (and pseudo-pressure) acts in
pseudo-time on an initially fully-resolved snapshot of turbulence according to
an auxiliary Stokes equation. The generalized velocity field is a function of
both physical time and pseudo-time. The effective equations for the
generalized velocity field in physical time (at a fixed pseudo-time
corresponding to numerical resolution $\ell\sim\sqrt{2\hat{\nu}\hat{t}}$)
provide the governing equations for LES. The equality of mixed partial
derivatives in physical time and pseudo-time provides a consistency condition
that defines the residual stress tensor based on the chosen pseudo-time
evolution equation.

In the simple case of unbounded flow with uniform numerical resolution, PIC is
shown to be equivalent to Gaussian filtering. Thus, the advantages of spatial
filtering are retained while providing a basis for defining a more generally
applicable coarsening procedure for complex flows. Even in the simple case of
unbounded flows with uniform resolution, the PIC approach yields fresh
insight. First, PIC leads to a previously-hidden exact representation of the
energy cascade in terms of multiscale velocity gradient interactions. This
provides a more precise framework for representing residual stresses based on
vortex stretching and other similar processes. Furthermore, the resulting
parabolic PDE in pseudo-time may be used to construct a Germano-like dynamic
procedure that does not require a test filter. This is demonstrated in the
creation of a dynamic eddy viscosity and dynamic mixed model based on PIC
equations. The performance of both models in a priori and a posteriori testing
is demonstrated in homogeneous isotropic turbulence. The dynamic eddy
viscosity model performs similarly to the commonly-used Smagnorinsky model and
the relative advantages of the dynamic mixed model are illustrated.

For realistic flows where nonuniform grid resolution is desirable, the PIC
framework for LES provides a nature-inspired mechanism for maintaining a
divergence-free condition for incompressible flows while avoiding commutation
errors (extra terms in the LES equations). In the case of wall-bounded flows,
boundary conditions may be chosen, PIC provides a convenient definition of
what it means to coarsen a flow in the vicinity of a domain boundary, and can
facilitate various modes of near-wall resolution or modeling treatment. PIC
theory is also easliy extended for a consistent treatment of grid anisotropy
by using an anisotropic pseudo-viscosity.

Finally, it is suggested that the physics-inspired framework for LES extends
more naturally to flows with more complex physics. A simple example of such an
extension is heat/mass transfer, for which the PIC approach to momentum is
straightforwardly copied. More significant extensions are also possible with
PIC, such as particle-laden or (interface-resolved) multiphase flows. While
PIC provides theoretical developments that can aid model development for
complex LES, it also provides an important foundation for enabling robust
data-driven modeling approaches. Indeed, data-driven closure techniques may
benefit significantly from the advantages of the PIC approach to LES theory,
for example, the recovery of divergence-free coarsened velocity field for
incompressible flows.

## Declaration of Interests

The author reports no conflict of interest.

## References

  * Bae & Lozano-Duran (2017) Bae, H. J. & Lozano-Duran, A. 2017 Towards exact subgrid-scale models for explicitly filtered large-eddy simulation of wall-bounded flows. Center for Turbulence Research Annual Research Briefs .
  * Bae & Lozano-Duran (2018) Bae, H. J. & Lozano-Duran, A. 2018 Dns-aided explicitly filtered les of channel flow. Certer for Turbulence Research Annual Research Briefs .
  * Bae & Lozano-Duran (2021) Bae, H. J. & Lozano-Duran, A. 2021 Effect of wall boundary conditions on a wall-modeled large-eddy simulation in a finite-difference framework. Fluids 6 (3).
  * Balachandar & Eaton (2010) Balachandar, S. & Eaton, J. K. 2010 Turbulent dispersed multiphase flow. Annual Review of Fluid Mechanics 42 (1), 111–133, arXiv: https://doi.org/10.1146/annurev.fluid.010908.165243.
  * Balachandar et al. (2019) Balachandar, S., Liu, K. & Lakhote, M. 2019 Self-induced velocity correction for improved drag estimation in euler–lagrange point-particle simulations. Journal of Computational Physics 376, 160–185.
  * Ballouz & Ouellette (2018) Ballouz, J. G. & Ouellette, N. T. 2018 Tensor geometry in the turbulent cascade. Journal of Fluid Mechanics 835, 1048–1064.
  * Bardina et al. (1980) Bardina, J., Ferziger, J. & Reynolds, W. 1980 Improved subgrid-scale models for large-eddy simulation. AIAA 13th Fluid & Plasma Dynamics Conference AIAA-80-1357, arXiv: https://arc.aiaa.org/doi/pdf/10.2514/6.1980-1357.
  * Bassenne et al. (2019) Bassenne, M., Esmaily, M., Livescu, D., Moin, P. & Urzay, J. 2019 A dynamic spectrally enriched subgrid-scale model for preferential concentration in particle-laden turbulence. International Journal of Multiphase Flow 116, 270–280.
  * Beck et al. (2019) Beck, A., Flad, D. & Munz, C.-D. 2019 Deep neural networks for data-driven les closure models. Journal of Computational Physics 398, 108910.
  * Bode et al. (2021) Bode, M., Gauding, M., Lian, Z., Denker, D., Davidovic, M., Kleinheinz, K., Jitsev, J. & Pitsch, H. 2021 Using physics-informed enhanced super-resolution generative adversarial networks for subfilter modeling in turbulent reactive flows. Proceedings of the Combustion Institute 38 (2), 2617–2625.
  * Borue & Orszag (1998) Borue, V. & Orszag, S. A. 1998 Local energy flux and subgrid-scale statistics in three-dimensional turbulence. Journal of Fluid Mechanics 366, 1–31.
  * Bose & Moin (2014) Bose, S. T. & Moin, P. 2014 A dynamic slip boundary condition for wall-modeled large-eddy simulation. Physics of Fluids 26 (1).
  * Bose & Park (2018) Bose, S. T. & Park, G. I. 2018 Wall-modeled large-eddy simulation for complex turbulent flows. Annual Review of Fluid Mechanics 50 (1), 535–561.
  * Bou-Zeid et al. (2005) Bou-Zeid, E., Meneveau, C. & Parlange, M. 2005 A scale-dependent lagrangian dynamic model for large eddy simulation of complex turbulent flows. Physics of Fluids 17 (2), 025105, arXiv: https://doi.org/10.1063/1.1839152.
  * Bull & Jameson (2016) Bull, J. R. & Jameson, A. 2016 Explicit filtering and exact reconstruction of the sub-filter stresses in large eddy simulation. Journal of Computational Physics 306, 117–136.
  * Capecelatro & Desjardins (2013) Capecelatro, J. & Desjardins, O. 2013 An euler–lagrange strategy for simulating particle-laden flows. Journal of Computational Physics 238, 1–31.
  * Chen et al. (2021) Chen, X., Lu, J. & Tryggvason, G. 2021 Interface retaining coarsening of multiphase flows. Physics of Fluids 33 (7), 073316, arXiv: https://doi.org/10.1063/5.0058776.
  * Cheng et al. (2019) Cheng, Y., Giometto, M., Kauffmann, P., Lin, L., Cao, C., Zupnick, C., Li, H., Li, Q., Abernathey, R. & Gentine, P. 2019 Deep learning for subgrid-scale turbulence modeling in large-eddy simulations of the atmospheric boundary layer. arXiv:1910.12125 .
  * Clark et al. (1979) Clark, R. A., Ferziger, J. H. & Reynolds, W. C. 1979 Evaluation of subgrid-scale models using an accurately simulated turbulent flow. J. Fluid Mech. 91 (part 1), 1–16.
  * Deardorff (1970) Deardorff, J. W. 1970 A numerical study of three-dimensional turbulent channel flow at large reynolds numbers. Journal of Fluid Mechanics 41 (2), 453–480.
  * Domaradzki & Loh (1999) Domaradzki, J. A. & Loh, K.-C. 1999 The subgrid-scale estimation model in the physical space representation. Physics of Fluids 11 (8), 2330–2342.
  * Donzis & Sreenivasan (2010) Donzis, D. A. & Sreenivasan, K. R. 2010 The bottleneck effect and the kolmogorov constant in isotropic turbulence. Journal of Fluid Mechanics 657, 171–188.
  * Drivas & Nguyen (2018) Drivas, T. D. & Nguyen, H. Q. 2018 Onsager’s conjecture and anomalous dissipation on domains with boundary. SIAM Journal on Mathematical Analysis 50 (5), 4785–4811, arXiv: https://doi.org/10.1137/18M1178864.
  * Duraisamy (2021) Duraisamy, K. 2021 Perspectives on machine learning-augmented reynolds-averaged and large eddy simulation models of turbulence. Phys. Rev. Fluids 6, 050504.
  * Falkovich (1994) Falkovich, G. 1994 Bottleneck phenomenon in developed turbulence. Physics of Fluids 6 (4), 1411–1414, arXiv: https://doi.org/10.1063/1.868255.
  * Fox (2003) Fox, R. O. 2003 Computational models for turbulent reacting flows. Cambridge university press.
  * Freund & Ferrante (2021) Freund, A. & Ferrante, A. 2021 Large-eddy simulation of droplet-laden decaying isotropic turbulence using artificial neural networks. International Journal of Multiphase Flow 142, 103704\.
  * Frohlich & von Terzi (2008) Frohlich, J. & von Terzi, D. 2008 Hybrid les/rans methods for the simulation of turbulent flows. Progress in Aerospace Sciences 44 (5), 349–377.
  * Gamahara & Hattori (2017) Gamahara, M. & Hattori, Y. 2017 Searching for turbulence models by artificial neural network. Phys. Rev. Fluids 2, 054604.
  * Germano (1986a) Germano, M. 1986a Differential filters for the large eddy numerical simulation of turbulent flows. Physics of Fluids 29 (6), 1755–1757.
  * Germano (1986b) Germano, M. 1986b Differential filters of elliptic type. Physics of Fluids 29 (6), 1757–1758, arXiv: https://aip.scitation.org/doi/pdf/10.1063/1.865650.
  * Germano (1992) Germano, M. 1992 Turbulence: the filtering approach. Journal of Fluid Mechanics 238.
  * Germano et al. (1991) Germano, M., Piomelli, U., Moin, P. & Cabot, W. H. 1991 A dynamic subgrid‐scale eddy viscosity model. Physics of Fluids A: Fluid Dynamics 3 (7), 1760–1765, arXiv: https://doi.org/10.1063/1.857955.
  * Ghate & Lele (2017) Ghate, A. S. & Lele, S. K. 2017 Subfilter-scale enrichment of planetary boundary layer large eddy simulation using discrete fourier–gabor modes. Journal of Fluid Mechanics 819, 494–539.
  * Ghate & Lele (2020) Ghate, A. S. & Lele, S. K. 2020 Gabor mode enrichment in large eddy simulations of turbulent flow. Journal of Fluid Mechanics 903, A13.
  * Ghosal & Moin (1995) Ghosal, S. & Moin, P. 1995 The basic equations for the large eddy simulation of turbulent flows in complex geometry. Journal of Computational Physics 118 (1), 24–37.
  * Haering et al. (2019) Haering, S. W., Lee, M. & Moser, R. D. 2019 Resolution-induced anisotropy in large-eddy simulations. Phys. Rev. Fluids 4, 114605.
  * Herrmann (2013) Herrmann, M. 2013 A sub-grid surface dynamics model for sub-filter surface tension induced interface dynamics. Computers and Fluids 87, 92–101.
  * Herrmann et al. (2018) Herrmann, M., Kedelty, D. & Ziegenhein, T. 2018 A dual-scale subgrid closure for les of phase interfaces in turbulent flows. In Proceedings of the Summer Program, Center for Turbulence Research, Stanford University, pp. 45–54.
  * Hinze (1955) Hinze, J. O. 1955 Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. AIChE journal 1 (3), 289–295.
  * Horwitz & Mani (2016) Horwitz, J.A.K. & Mani, A. 2016 Accurate calculation of stokes drag for point–particle tracking in two-way coupled flows. Journal of Computational Physics 318, 85–109.
  * Horwitz & Mani (2018) Horwitz, J.A.K. & Mani, A. 2018 Correction scheme for point-particle models applied to a nonlinear drag law in simulations of particle-fluid interaction. International Journal of Multiphase Flow 101, 74–84.
  * Ishihara et al. (2007) Ishihara, T., Kaneda, Y., Yokokawa, M., Itakura, K. & Uno, A. 2007 Small-scale statistics in high-resolution direct numerical simulation of turbulence: Reynolds number dependence of one-point velocity gradient statistics. Journal of Fluid Mechanics 592, 335–366.
  * Johnson (2020) Johnson, P. L. 2020 Energy transfer from large to small scales in turbulence by multi-scale nonlinear strain and vorticity interactions. Physical Review Letters 124, 104501, arXiv: 1912.00293.
  * Johnson (2021a) Johnson, P. L. 2021a On the role of vorticity stretching and strain self-amplification in the turbulence energy cascade. Journal of Fluid Mechanics 922, A3.
  * Johnson (2021b) Johnson, P. L. 2021b The squeezes, stretches, and whirls of turbulence. Physics Today 74 (4), 46–51.
  * Kim & Moin (2020) Kim, D. & Moin, P. 2020 Subgrid-scale Capillary Breakup Model for Liquid Jet Atomization. Combustion Science and Technology 192 (7), 1334–1357.
  * Kim et al. (2021) Kim, H., Kim, J., Won, S. & Lee, C. 2021 Unsupervised deep learning for super-resolution reconstruction of turbulence. Journal of Fluid Mechanics 910, A29.
  * Kim & Menon (1996) Kim, W.-W. & Menon, S. 1996 A new dynamic one-equation subgrid-scale model for large eddy simulations, arXiv: https://arc.aiaa.org/doi/pdf/10.2514/6.1995-356.
  * Kolmogorov (1941) Kolmogorov, A. N. 1941 The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers. Dokl. Akad. Nauk SSSR 30, 299–303.
  * Labourasse et al. (2007) Labourasse, E., Lacanette, D., Toutant, A., Lubin, P., Vincent, S., Lebaigue, O., Caltagirone, J. P. & Sagaut, P. 2007 Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests. International Journal of Multiphase Flow 33 (1), 1–39.
  * Langford & Moser (1999) Langford, J. A. & Moser, R. D. 1999 Optimal LES formulations for isotropic turbulence. Journal of Fluid Mechanics 398, 321–346.
  * Langford & Moser (2001) Langford, Jacob A. & Moser, Robert D. 2001 Breakdown of continuity in large-eddy simulation. Physics of Fluids 13 (5), 1524–1527, arXiv: https://doi.org/10.1063/1.1358876.
  * Larsson et al. (2016) Larsson, J., Kawai, S., Bodart, J. & Bermejo-Moreno, I. 2016 Large eddy simulation with modeled wall-stress: recent progress and future directions. Mechanical Engineering Reviews 3 (1), 15–00418.
  * Leonard (1975) Leonard, A. 1975 Energy cascade in large-eddy simulations of turbulent fluid flows. In Turbulent Diffusion in Environmental Pollution (ed. F.N. Frenkiel & R.E. Munn), Advances in Geophysics, vol. 18, pp. 237–248. Elsevier.
  * Leonard (2016) Leonard, A. 2016 Large-Eddy Simulation of the Navier-Stokes Equations: Deconvolution, Particle Methods, and Super-Resolution, pp. 1–16. World Scientific.
  * Lilly (1992) Lilly, D. K. 1992 A proposed modification of the germano subgrid‐scale closure method. Physics of Fluids 4 (3), 633–635, arXiv: https://doi.org/10.1063/1.858280.
  * Liu et al. (2020) Liu, B., Tang, J., Huang, H. & Lu, X.-Y. 2020 Deep learning methods for super-resolution reconstruction of turbulent flows. Physics of Fluids 32 (2), 025105, arXiv: https://doi.org/10.1063/1.5140772.
  * Lund & Rogers (1994) Lund, T. S. & Rogers, M. M. 1994 An improved measure of strain state probability in turbulent flows An improved measure . of strain state probability in turbulent flows. Physics of Fluids 6, 1838–1847.
  * Maulik & San (2017) Maulik, R. & San, O. 2017 A neural network approach for the blind deconvolution of turbulent flows. Journal of Fluid Mechanics 831, 151–181.
  * Maxey & Riley (1983) Maxey, M. R. & Riley, J. J. 1983 Equation of motion for a small rigid sphere in a nonuniform flow. The Physics of Fluids 26 (4), 883–889, arXiv: https://aip.scitation.org/doi/pdf/10.1063/1.864230.
  * McCaslin & Desjardins (2014) McCaslin, J. O. & Desjardins, O. 2014 Theoretical and computational modeling of turbulence / interface interactions. Center for Turbulence Research Proceedings of the Summer Program 2014 pp. 79–88.
  * Meneveau et al. (1996) Meneveau, C., Lund, T. S. & Cabot, W. H. 1996 A lagrangian dynamic subgrid-scale model of turbulence. Journal of Fluid Mechanics 319, 353–385.
  * Milano & Koumoutsakos (2002) Milano, M. & Koumoutsakos, P. 2002 Neural network modeling for near wall turbulent flow. Journal of Computational Physics 182 (1), 1–26.
  * Mockett et al. (2012) Mockett, C., Fuchs, M. & Thiele, F. 2012 Progress in des for wall-modelled les of complex internal flows. Computers and Fluids 65, 44–55, sixth International Conference on Computational Fluid Dynamics (ICCFD6).
  * Moin et al. (1991) Moin, P., Squires, K., Cabot, W. & Lee, S. 1991 A dynamic subgrid‐scale model for compressible turbulence and scalar transport. Physics of Fluids 3 (11), 2746–2757, arXiv: https://doi.org/10.1063/1.858164.
  * Monin & Yaglom (1975) Monin, A. S. & Yaglom, A. M. 1975 Statistical Fluid Mechanics: Mechanics of Turbulence, , vol. II. Dover Publications.
  * Moreau et al. (2006) Moreau, A., Teytaud, O. & Bertoglio, J. P. 2006 Optimal estimation for large-eddy simulation of turbulence and application to the analysis of subgrid models. Physics of Fluids 18 (10), 105101, arXiv: https://doi.org/10.1063/1.2357974.
  * Moser et al. (2021) Moser, R. D., Haering, S. W. & Yalla, G. R. 2021 Statistical properties of subgrid-scale turbulence models. Annual Review of Fluid Mechanics 53 (1), 255–286.
  * Park & Choi (2021) Park, J. & Choi, H. 2021 Toward neural-network-based large eddy simulation: application to turbulent channel flow. Journal of Fluid Mechanics 914, A16.
  * Patterson & Orszag (1971) Patterson, G. S. & Orszag, S. A. 1971 Spectral calculations of isotropic turbulence: Efficient removal of aliasing interactions. The Physics of Fluids 14 (11), 2538–2541, arXiv: https://aip.scitation.org/doi/pdf/10.1063/1.1693365.
  * Piomelli & Balaras (2002) Piomelli, U. & Balaras, E. 2002 Wall-layer models for large-eddy simluations. Annual Review of Fluid Mechanics 34 (1), 349–374, arXiv: https://doi.org/10.1146/annurev.fluid.34.082901.144919.
  * Pope (2000) Pope, S. B. 2000 Turbulent flows.
  * Pope (2010) Pope, S. B. 2010 Self-conditioned fields for large-eddy simulations of turbulent flows. Journal of Fluid Mechanics 652, 139–169.
  * Portwood et al. (2021) Portwood, G. D., Nadiga, B. T., Saenz, J. A. & Livescu, D. 2021 Interpreting neural network models of residual scalar flux. Journal of Fluid Mechanics 907, A23.
  * Prakash et al. (2021) Prakash, A., Jansen, K. E. & Evans, J. A. 2021 Invariant data-driven subgrid stress modeling in the strain-rate eigenframe for large eddy simulation. arXiv:2106.13410 .
  * Rozema et al. (2015) Rozema, W., Bae, H. J., Moin, P. & Verstappen, R. 2015 Minimum-dissipation models for large-eddy simulation. Physics of Fluids 27 (8), 1–15.
  * Sagaut (2006) Sagaut, P. 2006 Large Eddy Simulation for Incompressible Flows: An Introduction. Springer Berlin Heidelberg.
  * Sagaut & Germano (2005) Sagaut, P. & Germano, M. 2005 On the filtering paradigm for LES of flows with discontinuities. Journal of Turbulence 6 (23), 1–9.
  * Sarghini et al. (2003) Sarghini, F., de Felice, G. & Santini, S. 2003 Neural networks based subgrid scale modeling in large eddy simulations. Computers and Fluids 32 (1), 97–108.
  * Scotti & Meneveau (1999) Scotti, A & Meneveau, C 1999 A fractal model for large eddy simulation of turbulent flow. Physica D 127 (3-4), 198–232.
  * Scotti et al. (1993) Scotti, A., Meneveau, C. & Lilly, D. K. 1993 Generalized Smagorinsky model for anisotropic grids. Physics of Fluids A 5, 2306–2308.
  * Sirignano et al. (2020) Sirignano, J., MacArt, J. F. & Freund, J. B. 2020 Dpm: A deep learning pde augmentation method with application to large-eddy simulation. Journal of Computational Physics 423, 109811\.
  * Spalart (2009) Spalart, P. R. 2009 Detached-eddy simulation. Annual Review of Fluid Mechanics 41 (1), 181–202, arXiv: https://doi.org/10.1146/annurev.fluid.010908.165130.
  * Stoffer et al. (2021) Stoffer, R., van Leeuwen, C. M., Podareanu, D., Codreanu, V., Veerman, M. A., Janssens, M., Hartogensis, O. K. & van Heerwaarden, C. C. 2021 Development of a large-eddy simulation subgrid model based on artificial neural networks: a case study of turbulent channel flow. Geoscientific Model Development 14 (6), 3769–3788.
  * Stolz & Adams (1999) Stolz, S. & Adams, N. A. 1999 An approximate deconvolution procedure for large-eddy simulation. Physics of Fluids 11 (7), 1699–1701, arXiv: https://doi.org/10.1063/1.869867.
  * Tennekes & Lumley (1972) Tennekes, H. & Lumley, J. L. 1972 A First Course in Turbulence. Cambridge, MA: MIT Press.
  * Toutant et al. (2009) Toutant, A., Chandesris, M., Jamet, D. & Lebaigue, O. 2009 Jump conditions for filtered quantities at an under-resolved discontinuous interface. Part 1: Theoretical development. International Journal of Multiphase Flow 35 (12), 1100–1118.
  * Tryggvason & Lu (2020) Tryggvason, G. & Lu, J. 2020 Direct numerical simulations of multiphase flows: Opportunities and challenges. AIP Conference Proceedings 2293 (1), 030002.
  * Verstappen (2011) Verstappen, R. 2011 When does eddy viscosity damp subfilter scales sufficiently? Journal of Scientific Computing 49 (1), 94–110.
  * Vollant et al. (2017) Vollant, A., Balarac, G. & Corre, C. 2017 Subgrid-scale scalar flux modelling based on optimal estimation theory and machine-learning procedures. Journal of Turbulence 18 (9), 854–878, arXiv: https://doi.org/10.1080/14685248.2017.1334907.
  * Vreman (2004) Vreman, A. W. 2004 An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications. Physics of Fluids 16 (10), 3670–3681.
  * Vreman et al. (1994) Vreman, B., Geurts, B. & Kuerten, H. 1994 On the formulation of the dynamic mixed subgrid‐scale model. Physics of Fluids 6 (12), 4057–4059, arXiv: https://doi.org/10.1063/1.868333.
  * Wang et al. (2019) Wang, Q., Zhao, X. & Ihme, M. 2019 A regularized deconvolution model for sub-grid dispersion in large eddy simulation of turbulent spray flames. Combustion and Flame 207, 89–100.
  * Wang et al. (2021) Wang, Y., Yuan, Z., Xie, C. & Wang, J. 2021 Artificial neural network-based spatial gradient models for large-eddy simulation of turbulence. AIP Advances 11 (5), 055216, arXiv: https://doi.org/10.1063/5.0053590.
  * Wang et al. (2018) Wang, Z., Luo, K., Li, D., Tan, J. & Fan, J. 2018 Investigations of data-driven closure for subgrid-scale stress in large-eddy simulation. Physics of Fluids 30 (12), 125101, arXiv: https://doi.org/10.1063/1.5054835.
  * Xie et al. (2020a) Xie, C., Wang, J. & E, W. 2020a Modeling subgrid-scale forces by spatial artificial neural networks in large eddy simulation of turbulence. Phys. Rev. Fluids 5, 054606\.
  * Xie et al. (2020b) Xie, C., Yuan, Z. & Wang, J. 2020b Artificial neural network-based nonlinear algebraic models for large eddy simulation of turbulence. Physics of Fluids 32 (11), 115101, arXiv: https://doi.org/10.1063/5.0025138.
  * Yalla et al. (2021) Yalla, G. R., Oliver, T. A., Haering, S. W., Engquist, B. & Moser, R. D. 2021 Effects of resolution inhomogeneity in large-eddy simulation. Phys. Rev. Fluids 6, 074604.
  * Yang et al. (2019) Yang, X. I. A., Zafar, S., Wang, J.-X. & Xiao, H. 2019 Predictive large-eddy-simulation wall modeling via physics-informed neural networks. Phys. Rev. Fluids 4, 034602\.
  * Yuan et al. (2020) Yuan, Z., Xie, C. & Wang, J. 2020 Deconvolutional artificial neural network models for large eddy simulation of turbulence. Physics of Fluids 32 (11), 115106, arXiv: https://doi.org/10.1063/5.0027146.
  * Zhou et al. (2019) Zhou, Z., He, G., Wang, S. & Jin, G. 2019 Subgrid-scale model for large-eddy simulation of isotropic turbulent flows using an artificial neural network. Computers and Fluids 195, 104319.