# Derivation of a new LES model approximated from exact two-point equations and evaluation in a Taylor-Green flow

Paul Beaumard, Jean Philippe Laval, J. Vassilicos

► To cite this version:

Paul Beaumard, Jean Philippe Laval, J. Vassilicos. Derivation of a new LES model approximated from exact two-point equations and evaluation in a Taylor-Green flow. Physical Review Fluids, 2026, 11 (1), pp.014607. 10.1103/pkyq-vtd6. hal-05517564

HAL Id: hal-05517564

https://hal.science/hal-05517564v1

Submitted on 18 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Derivation of a new LES model approximated from exact two-point equations and evaluation in a Taylor-Green flow

P. Beaumard

Univ. Lille, CNRS, ONERA, Arts et Metiers Institute of Technology, Centrale Lille,
UMR 9014 - LMFL - Laboratoire de Mécanique des Fluides de Lille - Kampé de Fériet, F-59000 Lille, France. and
DAAA, ONERA, Institut Polytechnique de Paris, 92190, Meudon, France.

J. P. Laval

Univ. Lille, CNRS, ONERA, Arts et Metiers Institute of Technology, Centrale Lille,
UMR 9014 - LMFL - Laboratoire de Mécanique des Fluides de Lille - Kampé de Fériet, F-59000 Lille, France.

J. C. Vassilicos

Univ. Lille, CNRS, ONERA, Arts et Metiers Institute of Technology, Centrale Lille,
UMR 9014 - LMFL - Laboratoire de Mécanique des Fluides de Lille - Kampé de Fériet, F-59000 Lille, France.

(Dated: December 10, 2025)

The Kármán-Howarth-Monin-Hill equation has been used by various authors to quantify and understand the small scale turbulent energy cascade without assumptions and therefore even in the case of non-homogeneous turbulence. A complementary equation was introduced by Germano [1] to describe the large scales dynamics of the flow. In this publication, the links between the Germano equation and Large Eddy Simulation (LES) are examined and an approximation of an exact equation is used to design a new LES model which captures the physical energy transfer between filtered scales and residual subfilter scales. This new model, based on simple simplifications of exact equations, is tested both a priori and a posteriori in a Taylor-Green flow. The model is found to better represent the local energy transfers compared to a similar model introduced in [2] especially for both extreme forward and backscatter events. The Bardina et al. scale similarity model [3] over-performs our model in the a priori analysis, meaning that refining the mathematical approximations used for the derivation of our model could be useful, but our model's performance is comparable if not even better than the scale similarity model in a mixed configuration including a turbulent viscosity model. Therefore, our model, derived with simple simplifications of the exact two-point equations, already performs similarly to LES models classically used for simulations. This publication is a proof of concept that two-point equations can be used to develop new LES models and it may be the right direction to develop more efficient models derived from the Navier-Stokes equations.

# I. INTRODUCTION

Numerical simulation is now routinely used for the design and analysis of flows in practical applications. However, in most of these applications, the simulation of all turbulent scales is not possible with accessible computational resources. Unfortunately, turbulent flows are defined by their wide range of scales co-existing and interacting together. In many applications, the small scales of the turbulence (even significantly larger than the Kolmogorov microscale [4]) can be very small compared to the scales of interest of the simulation. Following the description of the turbulent energy cascade from Richardson [5], formulated mathematically for homogeneous, stationary turbulence by Kolmogorov [4], the small eddies of the flow play a significant role in turbulent flow behavior as they dissipate the kinetic energy into heat. As a consequence, the energy cascades from large scales to small scales to feed the dissipated energy at small scales. This robust description of the physics underlines the impossibility of neglecting the dynamics of the small scales of the flow and not taking into account their effect on the large scales.

In Large Eddy Simulations (LES), only the large scales, or more precisely the low-pass filtered velocity field components, are resolved. This reduces significantly the computational cost compared to Direct Numerical Simulations (DNS) where all the scales are taken into account in the simulation via the well-resolved Navier-Stokes equations. Large Eddy Simulations rely on a subfilter scale model whose role is to model the non-resolved scales which are expected to, at least, dissipate the energy of the flow.

In the LES concept the modelling and the numerical methods should be independent however they are closely related in practice [6-8]. The turbulent models are usually designed for a specific filter used to filter the Navier-Stokes equation.

$$ \frac{\partial \bar{\mathbf{u}}}{\partial t} + (\bar{\mathbf{u}} \cdot \nabla) \bar{\mathbf{u}} = - \frac{\partial}{\partial x_j} \tau_{ij} - \frac{1}{\rho} \nabla \bar{p} + \nu \nabla^2 \bar{\mathbf{u}}, \quad (1) $$
where the residual stress $\tau_{ij} = \overline{u_i u_j} - \overline{u_i} \overline{u_j}$ must be modeled.

The modelling errors are associated with the difference between $\tau_{ij}$ and its model. On the other hand, the numerical error comes from the numerical resolution of equation (1) where $\tau_{ij}$ is modelled and can be evaluated through mesh convergence for a given filter. In this study we focus on the modelling error only.

An energy equation for the filtered velocity can be derived from equation (1):

$$ \frac{\partial}{\partial t} \left( \frac{\overline{u_i u_j}}{2} \right) + \overline{u_j} \frac{\partial}{\partial x_j} \left( \frac{\overline{u_i u_j}}{2} \right) = - \overline{u_i} \frac{\partial \tau_{ij}}{\partial x_j} - \frac{1}{\rho} \frac{\partial (\overline{u_i p})}{\partial x_i} + \nu \overline{u_i} \frac{\partial^2 \overline{u_i}}{\partial x_j \partial x_j} \quad (2) $$

where $T^e = -\overline{u_i} \frac{\partial \tau_{ij}}{\partial x_j}$ is the rate of energy exchange between the filtered scales and the residual sub-filter scales. It represents a subfilter stress energy contribution to equation (2) and is usually re-written as:

$$ T^e = -\overline{u_i} \frac{\partial \tau_{ij}}{\partial x_j} = -\frac{\partial}{\partial x_j} (\overline{u_i} \tau_{ij}) + \overline{S_{ij}} \tau_{ij}, \quad (3) $$

where $\overline{S_{ij}} = \frac{1}{2} \left( \frac{\partial \overline{u_i}}{\partial x_j} + \frac{\partial \overline{u_j}}{\partial x_i} \right)$ is the deformation tensor computed with filtered velocity. The term $\frac{\partial}{\partial x_j} (\overline{u_i} \tau_{ij})$ is usually interpreted as a spatial flux and $\overline{S_{ij}} \tau_{ij}$ as the transfer of energy across scales between filtered scales and residual scales: therefore closely related to the energy cascade [9]. A negative value of $\overline{S_{ij}} \tau_{ij}$ is usually interpreted as a direct energy transfer from large scales to small scales and the opposite is called backscatter [10]. In order to evaluate the performance of the LES models, we will compare directly the rate of energy exchange $T^e$ in the energy equation (2) with its modelled counterpart.

Following the concept of turbulent eddy viscosity (Boussinesq hypothesis), the non-resolved scales may be modeled as an additional viscosity acting on the flow. The most commonly used model based on the concept of turbulent viscosity was defined by Smagorinsky [11]. In this model the subgrid scale tensor is modelled as:

$$ \tau_{ij} - \frac{1}{3} \tau_{kk} \delta_{ij} = -2 \nu_T \overline{S_{ij}} \quad (4) $$

where:

$$ \nu_T = (C_S \Delta)^2 |\overline{S}|. \quad (5) $$

is the turbulent viscosity, $|\overline{S}| \equiv \sqrt{2\overline{S_{ij}}\overline{S_{ij}}}$, $\Delta$ is the filter width and $C_S$ is a constant to adapt.

This rather simple model is able to compensate part of the unresolved dissipation from the subfilter scales despite important limitations having been identified (see e.g. [8, 12]). In these publications, the Smagorinsky model is found to be unable to ensure regularisation of the solution, meaning that the solution does not converge to the right solution when the mesh used to discretise the equation is refined keeping the filter unchanged. Moreover, the modelled subgrid scale tensor $\tau_{ij}$ is known to not be representative of the real $\tau_{ij}$ evaluated by a-priori tests from DNS results [13]. Finally, this model is known to produce no backscatter given that the local energy flux $\tau_{ij} \overline{S_{ij}}$ in (3) is modelled as $-2\nu_T \overline{S_{ij}} \overline{S_{ij}}$ which is negative as $\nu_T > 0$. Backscatter events, already described in [10], appear locally in turbulent flows [2] and can even prevail on average over dissipation in strongly non-homogeneous turbulent flows such as in the buffer layer close to the wall of a channel flow [14].

A dynamic procedure [15–17] can be used to modulate the effect of any model with a constant depending on local flow properties by replacing the global coefficient, for example $C_S$, by a dynamic coefficient computed locally. This method was designed to extend the use of the Smagorinsky model to non-homogeneous flows. It improves significantly the results in specific situations. For example, it reduces dissipation near walls where the original Smagorinsky model produces an excess of dissipation [15]. However, the dynamic procedure is just a correction but does not solve all the problems associated with the essential hypothesis of the model. For instance, a Smagorinsky model with such a dynamic coefficient still cannot produce backscatter.

Another strategy is refered to as implicit LES in which a regularisation is introduced as a numerical dissipation through the discretization of the governing equations or extra discrete operator that damps selectively the smallest scales. For instance, in the implicit model developed by Dairay et al. [8], the numerical dissipation is introduced by the discretization of the viscous term specifically designed to introduce the desired level of dissipation. This method is equivalent to the use of spectral vanishing viscosity with an implementation in real space. This method was shown to be more efficient in reducing the number of degrees of freedom compared to the Smagorinsky model which is sensitive to numerical errors and which requires refinement of the mesh with respect to the filter width.
Other models were introduced dedicated to specific flows. This is for example the case of wall-adapted LES which has been widely developed these past years to capture non dissipative effects close to walls. However, we focus in this study on universal models such as the scale similarity model introduced by [18] where $\tau_{ij}$ is modelled as:

$$ \tau_{ij} = C_3(\overline{\tilde{u}_i \tilde{u}_j} - \tilde{u}_i \tilde{u}_j) \qquad (6) $$

where $\tilde{\cdot}(\cdot)$ is a filter applied on the resolved quantities. This model is a generalization of the Bardina's model [3] with two consecutive filters of different shape and filter width.

According to Meneveau and Katz [13], the scale similarity model (6) correlates much better than the Smagorinsky model with the exact subfilter stress once evaluated on real turbulence data. However, the mean dissipation is slightly under-predicted which makes the model practically unusable in this formulation. This is why this model is often used as a complement to a dissipative model such as Smagorinsky in a mixed configuration so that $\tau_{ij}$ is modelled as:

$$ \tau_{ij} = C_3(\overline{\tilde{u}_i \tilde{u}_j} - \tilde{u}_i \tilde{u}_j) - 2\nu_T \overline{S_{ij}} \qquad (7) $$

where, $\nu_T$ is defined in (5). The principle of the scale similarity model is to evaluate a subfilter stress with resolved quantities by over-filtering the resolved terms and assume a similarity between the subfilter stress to be modelled and the equivalent one built on the resolved scales. This method is expected to give good results when the filter width is in the inertial range. The mixed model configuration allows to tune the two parts of the model separately: one mainly taking into account most of the dissipation and the other where the local energy transfer includes local dissipation and backscatter. The scale similarity model (6) is an efficient and practical model. However, the physical justification of this model remains weak and does not improve our understanding of the physics captured by the model.

Another class of models can be obtained by truncated Taylor expansion of the subgrid stress tensor (which may be valid only for very small $\Delta$). Applying the procedure with a Gaussian filter leads to the gradient model [19] analysed in [20]. This model can also be combined in a mixed configuration with the Smagorinsky model so that $\tau_{ij}$ is modelled as:

$$ \tau_{ij} = \frac{\Delta^2}{12} \frac{\partial \overline{u_i}}{\partial x_k} \frac{\partial \overline{u_j}}{\partial x_k} - 2\nu_T \overline{S_{ij}}. \qquad (8) $$

This model is known to allow backscatter through the first gradient term. The recent work from Cimarelli [21] continues the development of this model and relates it to the structure functions following Germano's framework [22] presented in section III of the present paper. In this study, a new decomposition between small and large scales of the subgrid stress tensor is proposed and a modelling particularly adapted for non-homogeneous grid meshes is tested with promising results in a channel flow.

The need to model backscatter events is questioned by Vela-Martín [9] who suggests that backscatter does not need to be modelled in homogeneous turbulent flows. The argument is to say that, from an energetic perspective, addition of any divergence-free tensor to obtain an alternative representation of the subgrid stress tensor is physically valid and produces the same average SGS flux but in general different spatial distributions of the local SGS energy fluxes. The backscatter is defined, as usual, as negative events of $\tau_{ij} \overline{S_{ij}}$ which does not take into account the full energy transfer between filtered scales and residual scales. The non-homogeneous term $-\frac{\partial}{\partial x_j}(\overline{u_i} \tau_{ij})$ in equation (3) may affect the sign of the full contribution, $T^e$, locally or globally in non-homogeneous flows. Indeed, positive and negative values are produced by the non-homogeneous term $-\frac{\partial(\overline{u_i} \tau_{ij})}{\partial x_j}$ so that the full subfilter stress energy contribution $T^e$ has positive and negative contributions to equation (2) even when $\tau_{ij} \overline{S_{ij}}$ does not produce backscatter. The lack of connection between the physics of the local energy transfer and subgrid scale modelling is crucially missing (see [23–25]) and leads to very fundamental questions such as: do we need to model local backscatter effects? The Germano equation [1] might give an avenue towards a potential answer to this question. This exact equation describes the larger scale behavior and is analyzed experimentally in [26]. To our knowledge, the physics associated with this equation has not been analysed except by Germano. However, various studies of its complementary small scale equation exist [24, 25, 27–30]. We refer to this equation as the Kármán-Howarth-Monin-Hill (KHMH) equation and we introduce it in the following section. The KHMH equation is very useful to understand the energy transfer in non-homogeneous flows without any assumption. In [31], the relations between the interscale energy transfer rate of the KHMH equation (integrated in scale) and the rate of energy transferred between filtered scales and residual scales are analyzed. An analysis based on a DNS of isotropic turbulence at high Reynolds number ($Re_\lambda = 1250$) shows that both quantities share many properties but are not identical. The use of the larger scale Germano's equation [1] instead of the KHMH
<figure><img src="image_2.png" /><figcaption>FIG. 1: Schematic of fluid velocities at points $\zeta^{-} = X - r$ and $\zeta^{+} = X + r$ used for the derivation of the two-point equations (9) and (11).</figcaption></figure>

equation may allow to explain these common and different properties. These equations were also used in [32] to derive the different terms of the subgrid stress gradient which should be modelled, including some non-homogeneous terms. These authors show exactly why the links are not exact between KHMH interscale energy transfer rate and energy exchanges between filtered scales and residual scales. The Germano equation is also likely to be a powerful tool for understanding some physical constraints in subfilter scale models demonstrating the potential of the two-point equation for the simulation of practical applications for engineering.

## II. TWO-POINT STATISTICS THEORY

Interscale turbulence transfers for incompressible turbulence can be studied in the presence of all other co-existing turbulence transfer/transport mechanisms in terms of two-point equations exactly derived from the incompressible Navier-Stokes equations (see [1, 28, 33]) without any hypotheses or assumptions, in particular no assumptions of homogeneity or periodicity. The incompressible Navier-Stokes equation is written at two points $\zeta^{-} = X - r$ and $\zeta^{+} = X + r$ in physical space (see figure 1) where $X$ is the centroid and $r$ is the two-point separation vector. One defines the two-point velocity half difference $\delta u(X, r, t) \equiv \frac{u^{+} - u^{-}}{2}$ where $u^{+} \equiv u(\zeta^{+})$ and $u^{-} \equiv u(\zeta^{-})$ are the fluid velocities at each one of the two points and the two-point pressure half difference $\delta p(X, r, t) \equiv \frac{p^{+} - p^{-}}{2}$ where $p^{+} \equiv p(\zeta^{+})$ and $p^{-} \equiv p(\zeta^{-})$ are the pressure over density ratios at each one of the two points. Incompressibility immediately imposes $\nabla_X \delta u = \nabla_r \delta u = 0$ and the momentum equation leads to (see [28, 33])

$$ \frac{\partial \delta u}{\partial t} + (u_X \cdot \nabla_X) \delta u + (\delta u \cdot \nabla_r) \delta u = - \nabla_X \delta p + \frac{\nu}{2} \nabla_X^2 \delta u + \frac{\nu}{2} \nabla_r^2 \delta u \quad (9) $$

where $u_X(X, r, t) \equiv \frac{u^{+} + u^{-}}{2}$; $\nabla_X$ and $\nabla_X^2$ are the gradient and Laplacian in $X$ space; $\nabla_r$ and $\nabla_r^2$ are the gradient and Laplacian in $r$ space; and $\nu$ is the kinematic viscosity.

An energy equation is readily obtained by multiplying equation (9) by $2\delta u$:

$$ \frac{\partial |\delta u|^2}{\partial t} + \nabla_X \cdot (u_X |\delta u|^2) + \nabla_r \cdot (\delta u |\delta u|^2) = -2 \nabla_X \cdot (\delta u \delta p) + \frac{\nu}{2} \nabla_X^2 |\delta u|^2 + \frac{\nu}{2} \nabla_r^2 |\delta u|^2 - \frac{1}{2} \epsilon^+ - \frac{1}{2} \epsilon^- \quad (10) $$

where $\epsilon^{+} = \nu \frac{\partial u_{+}^{+}}{\partial \zeta_{k}^{+}} \frac{\partial u_{+}^{+}}{\partial \zeta_{k}^{+}}$ and $\epsilon^{-} = \nu \frac{\partial u_{-}^{-}}{\partial \zeta_{k}^{-}} \frac{\partial u_{-}^{-}}{\partial \zeta_{k}^{-}}$. In [30], equation (10), is integrated in scale space $(\iiint_{0}^{R} (.) dr^3)$ to clearly demonstrate that it represents the energy equation of the scales smaller than $R$. Therefore, $|\delta u|^2$ is related to a small scale quantity.

A complementary equation can be derived for the evolution for the two-point velocity half sum $u_X(X, r, t) \equiv \frac{u^{+} + u^{-}}{2}$. This equation was first obtained in [1]:

$$ \frac{\partial u_X}{\partial t} + (u_X \cdot \nabla_X) u_X + (\delta u \cdot \nabla_r) u_X = - \nabla_X p_X + \frac{\nu}{2} \nabla_X^2 u_X + \frac{\nu}{2} \nabla_r^2 u_X \quad (11) $$
where $p_X \equiv \frac{p^+ + p^-}{2}$, and the incompressibility means $\nabla_X \cdot u_X = \nabla_r \cdot u_X = 0$. An energy equation, also derived for the first time in [1], is readily obtained by multiplying equation (11) with $2u_X$:

$$ \frac{\partial |u_X|^2}{\partial t} + \nabla_X \cdot (u_X |u_X|^2) + \nabla_r \cdot (\delta u |u_X|^2) = -2 \nabla_X \cdot (u_X P_X) \\ + \frac{\nu}{2} \nabla_X^2 |u_X|^2 + \frac{\nu}{2} \nabla_r^2 |u_X|^2 - \frac{1}{2} \epsilon^+ - \frac{1}{2} \epsilon^- \quad (12) $$

In order to interpret the energy term $|u_X|^2$, the following equality can be derived:

$$ \iiint_0^R |\delta u|^2 dr^3 + \iiint_0^R |u_X|^2 dr^3 = \iiint_{|r|<R} \frac{1}{2} u^2 (x+r) dr^3. \quad (13) $$

The sum of the energy contribution of the scales smaller than $R$ ($\iiint_0^R |\delta u|^2 dr^3$) with $\iiint_0^R |u_X|^2 dr^3$ balances the total kinetic energy in a sphere of radius $R$. Therefore, $\iiint_0^R |u_X|^2 dr^3$ can be interpreted as the energy of the scales larger than $R$ and equation (12) as a equation for large scales. The connection between large eddy simulation framework and the two point statistics framework is described in more details in section III.

The two-point energy equations above involve important interscale and interspace transport terms. Germano interpreted the equations (11) and (12) in the context of large eddy simulations [1]. He showed that the term $(\delta u \cdot \nabla_r) u_X$ in (11) can be interpreted as the gradient of a subgrid stress. This term gives rise to the term $\nabla_r \cdot (\delta u |u_X|^2)$ in (12) which is therefore an energy transfer rate across the scales of $|u_X|^2$. Germano also derived the kinematic equation [1]

$$ \nabla_r \cdot (\delta u |u_X|^2) + \nabla_r \cdot (\delta u |\delta u|^2) = 2 \nabla_X \cdot (\delta u (\delta u \cdot u_X)) \quad (14) $$

which relates $\nabla_r \cdot (\delta u |u_X|^2)$ to $\nabla_r \cdot (\delta u |\delta u|^2)$ in (10) and where $\nabla_r \cdot (\delta u |\delta u|^2)$ accounts for non-linear interscale energy transfer and the turbulence cascade, e.g. see [30].

In homogeneous/periodic turbulence, the large scale interscale transfer rate $\nabla_r \cdot (\delta u |u_X|^2)$, where $\langle \cdot \rangle$ represents space averaging over all space or a triply periodic domain, is exactly equal to the small-scale interscale transfer rate $\nabla_r \cdot (\delta u |\delta u|^2)$ and balanced exactly with the turbulence dissipation rate:

$$ \nabla_r \cdot (\delta u |u_X|^2) = - \nabla_r \cdot (\delta u |\delta u|^2) = \langle \epsilon \rangle \quad (15) $$

These results cannot be derived for non-homogeneous turbulence but the small-scale and large-scale interscale transfer rate estimates measured in a mixing tank and averaged in time and locally in space have been shown to collapse with dissipation rate even though the flow is non-homogeneous [26]. These results are supported theoretically in the same publication (based on the theory of Chen et al [30]) for the small-scale interscale transfer but only observed experimentally for the large-scale interscale transfer. This means that the mean values of the time-averaged interscale transfer rates seem to be of the same order of magnitude as the dissipation rate even for some non-homogeneous flows. However, Yasuda and Vassilicos [34] quantified the small scale interscale transfer spatio-temporal fluctuations in a Direct Numerical Simulation of homogeneous turbulence. These fluctuations are very large with rare extreme events much larger than the mean dissipation value. These results are reminiscent of the backscatter already described in [10, 35–39]. Moreover, in [34], these rare event contributions cannot be neglected as they contribute significantly to the mean interscale transfer rate value (50% of the mean value is contributed by events whose probabilities are less than few percent of the probabilities of the most likely events for $|r|$ of the order of the Taylor scale $\lambda$).

The probability distribution function of the truncated estimate of the small-scale interscale transfer rate normalized by dissipation is computed from the experimental data presented in [26] and plotted in figure 2a. This dataset comes from 2D2C PIV measurements in the bulk of a mixing tank where baffles are used to break the flow rotation. These results are computed with 150 000 velocity fields of 2D2C PIV measurements with interrogation window size between 3.4$\eta$ and 5.1$\eta$. Very large tails of the PDF are measured, with positive and negative events much larger than dissipation (larger than 1000$\langle \epsilon' \rangle$ for $r_z > \lambda$). The local behavior is therefore very different from its mean value which is of the same order of magnitude as the mean dissipation (see [26]). These results confirm that, locally, the turbulence behavior is far from the Kolmogorov picture with both direct and backward energy transfers across scales much larger than dissipation. The spatio-temporal probability distribution function of the large-scale interscale transfer rate is also measured and plotted in figure 2b. The fluctuations of this term are similar and even larger than for the small-scale interscale transfer rate with very large tails and values larger than 1500$\langle \epsilon' \rangle$ for the various separation lengths plotted. In this publication we attempt to incorporate this property of large-scale interscale transfer rates in Large Eddy Simulation modeling. The very large fluctuations of the large-scale interscale transfer rates, both direct and backward, needs to be reproduced by a model which captures the main local properties of turbulence.
<figure><img src="image_3.png" /></figure>

$$ (a) \frac{\frac{\partial}{\partial r_x} (\delta u_x' (\delta u_x'^2 + \delta u_z'^2)) + \frac{\partial}{\partial r_z} (\delta u_z' (\delta u_z'^2 + \delta u_x'^2))}{\langle \epsilon' \rangle} $$

<figure><img src="image_4.png" /></figure>

$$ (b) \frac{\frac{\partial}{\partial r_x} (\delta u_x' (u_{Xx}^{\prime 2} + u_{Xz}^{\prime 2})) + \frac{\partial}{\partial r_z} (\delta u_z' (u_{Xx}^{\prime 2} + u_{Xz}^{\prime 2}))}{\langle \epsilon' \rangle} $$

FIG. 2: Spatio-temporal probability distribution function of the normalized small scale (a) and large scale (b) truncate estimates of the interscale transfer at different separation length $r_z$. These results are measured during experiments in a mixing tank with four rotating blades and with baffles (vertical bars that break the flow rotation). More details about the experiments are available in [26]. The statistics are computed with 150 000 velocity fields of 2D2C PIV measurements with interrogation window size between 3.4η and 5.17η. The red dashed line corresponds to the theoretical mean value in homogeneous turbulence (−1 in figure (a) and +1 in figure (b)).

### III. LINKS BETWEEN TWO-POINT EQUATIONS AND LES MODELLING

In the following sub-section, the mathematical links between the two-point equations and the Large Eddy Simulation framework are derived following [40].

#### A. Derivation of an exact form of the subfilter stress

The notations of $u_X$ and $\delta u$, defined in section II, are slightly different to the ones used in [40]. Therefore, Germano's results are converted into our set of notation. Using relations from [1] with the new notations, the large scale two-point equation (11) is rewritten as:

$$ \frac{\partial u_{Xi}}{\partial t} + u_{Xj} \frac{\partial u_{Xi}}{\partial X_j} = - \frac{\partial (\delta u_i \delta u_j)}{\partial X_j} - \frac{1}{\rho} \frac{\partial P_X}{\partial X_i} + \nu \frac{\partial^2 u_{Xi}}{\partial X_k \partial X_k} \quad (16) $$

Following [22], the term $-\frac{\partial(\delta u_i \delta u_j)}{\partial X_j}$ represents the subfilter stress gradient associated to the two-point filter:

$$ u_X = \frac{u(X+r) + u(X-r)}{2} = \int G_0(X-\xi)u(\xi)d\xi \quad (17) $$

where

$$ G_0(X - \xi) = \frac{1}{2} (\delta(X + r - \xi) + \delta(X - r - \xi)), \quad (18) $$

and $\delta$ is the Dirac delta function $\delta(x - \xi) = \prod_{k=1:3} \delta(x_k - \xi_k)$.

An exact form of the subfilter stress (which can, however, not be used as a subfilter stress model) derived from the two-point large scale equation (16) for the basic filter (18) is

$$ \tau_{ij}^{G_0} = \delta u_i \delta u_j \quad (19) $$
which is the basis of second-order structure functions. Therefore, equation (16) connects exactly large eddy simulation framework with the two point statistics framework introduced in section II.

A more practically useful equation for the subfilter stress $\tau_{ij}$ is derived in [40] and [22] directly from the Reynolds stress formulation. The filtered velocity $\bar{u}$ is obtained as

$$ \bar{u}_i = \int_{-\infty}^{\infty} G(\boldsymbol{x} - \boldsymbol{\xi}) u_i(\boldsymbol{\xi}) d\boldsymbol{\xi}. \qquad (20) $$

where the discrete weight function

$$ G(\boldsymbol{x} - \boldsymbol{\xi}) = \sum_{\alpha} g_{\alpha} \delta(\boldsymbol{x} + \boldsymbol{r}_{\alpha} - \boldsymbol{\xi}) \qquad (21) $$

is used, with the weights $g_{\alpha}$ such that $\sum_{\alpha} g_{\alpha} = 1$ and $\int G(\boldsymbol{r})d\boldsymbol{r} = 1$. The choice of $\alpha$ and $r_{\alpha}$ enable the implementation of different discrete filters but the following derivation from [22] does not depend on these parameters. Following [22] and [41] (equations (33) to (38)) the exact subfilter stress formula is derived as:

$$ \begin{aligned} \tau_{ij}(\boldsymbol{x}) &= \bar{u}_i \bar{u}_j - \bar{u}_i \bar{u}_j \\ &= \int_{-\infty}^{\infty} G(\boldsymbol{x} - \boldsymbol{\xi}) u_i(\boldsymbol{\xi}) u_j(\boldsymbol{\xi}) d\boldsymbol{\xi} - \int_{-\infty}^{\infty} G(\boldsymbol{x} - \boldsymbol{\xi}) u_i(\boldsymbol{\xi}) d\boldsymbol{\xi} \int_{-\infty}^{\infty} G(\boldsymbol{x} - \boldsymbol{\xi}') u_j(\boldsymbol{\xi}') d\boldsymbol{\xi}' \\ &= \frac{1}{2} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(\boldsymbol{x} - \boldsymbol{\xi}) G(\boldsymbol{x} - \boldsymbol{\xi}') (u_i(\boldsymbol{\xi}) - u_i(\boldsymbol{\xi}')) (u_j(\boldsymbol{\xi}) - u_j(\boldsymbol{\xi}')). \end{aligned} \qquad (22) $$

Using the change of variables $\boldsymbol{r} = \frac{\boldsymbol{\xi} - \boldsymbol{\xi}'}{2}$ and $\boldsymbol{s} = \frac{\boldsymbol{\xi} + \boldsymbol{\xi}'}{2}$

$$ \tau_{ij}(\boldsymbol{x}) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(\boldsymbol{x} - \boldsymbol{s} - \boldsymbol{r}) G(\boldsymbol{x} - \boldsymbol{s} + \boldsymbol{r}) \delta u_i(\boldsymbol{r}, \boldsymbol{s}) \delta u_j(\boldsymbol{r}, \boldsymbol{s}) dr ds \qquad (23) $$

where $\delta u_i(\boldsymbol{r}, \boldsymbol{s}) \delta u_j(\boldsymbol{r}, \boldsymbol{s})$ has the same definition as in (19), introduced in section II. The associated subfilter force is also derived (see [1], where the mathematical relation: $\frac{\partial \delta u_i}{\partial x_j} = \frac{\partial u_{X_i}}{\partial r_j}$ is used and $G$ is assumed to be homogeneous in space) and rewritten as:

$$ \frac{\partial \tau_{ij}}{\partial x_j} = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(\boldsymbol{x} - \boldsymbol{s} - \boldsymbol{r}) G(\boldsymbol{x} - \boldsymbol{s} + \boldsymbol{r}) \nabla_r \cdot (u_{X_i} \delta \boldsymbol{u}) dr ds \qquad (24) $$

so that the term $\nabla_r \cdot (u_{X_i} \delta \boldsymbol{u})$ which is at the origin of the large-scale interscale transfer rate $\nabla_r \cdot (\delta \boldsymbol{u} | u_{X_i} |^2)$ measured in [26] with PIV experiments, appears explicitly. Indeed, the large scale interscale transfer rate appears in equation (12) through the relation:

$$ \nabla_r \cdot (\delta \boldsymbol{u} | u_{X_i} |^2) = 2 u_{X_i} \nabla_r \cdot (u_{X_i} \delta \boldsymbol{u}) . \qquad (25) $$

The exact results (23) and (24) are very meaningful for the understanding of the subfilter stress properties in LES equations. Firstly, the subfilter stress has a direct link with the structure functions which are analyzed in many flows and for which Kolmogorov's 1941 theory was developed, thereby providing a direct connection with small-scale turbulence physics. Secondly, the term $\nabla_r \cdot (u_{X_i} \delta \boldsymbol{u})$ present in the exact subfilter stress gradient formulation (24) is closely related to the large-scale interscale transfer rate as already mentioned. Again, the physical understanding developed with the small-scale $(\delta \boldsymbol{u})^2$ two-point equation and large-scale $(u_{X_i}^2)$ two-point equation can be directly used for subfilter stress modelling in the LES framework.

## B. Existing LES models close to Germano exact subfilter stress equation

Several models close to Germano's formulation (23) have been proposed and can be used for comparison. For instance, the following increment model was introduced in [2]:

$$ \frac{\partial \tau_{ij}}{\partial x_j} = C_b \frac{\partial (\delta u_i \delta u_j)}{\partial x_j} \qquad (26) $$
where

$$ \tilde{\delta}u_i = \sum_{k=1}^{3} (\bar{u}_i(\boldsymbol{x} + m\Delta_k e_k) - \bar{u}_i(\boldsymbol{x} - m\Delta_k e_k)), \qquad (27) $$

and $C_b$ is a dimensionless coefficient. In the case of a homogeneous grid (as in the present paper's simulations), $\Delta_k = \Delta$ for all $k, m$ is an integer number and $e_k$ is the unit vector in direction $k$. This large eddy simulation model is not explicitly inspired by Germano [22] as it was published earlier but it has a striking similarity with the exact Germano equation (23). A dynamic procedure was defined by the authors to compute the model coefficient $C_b$ in (26). The model was applied to a channel flow and tested both a priori and a posteriori. Very good prediction of the probability distribution function of $\tau_{ij}\bar{S}_{ij}$, including backscatter, was observed in a priori analysis in the buffer layer of a channel flow at $Re_{\tau} = 180$. However, the filtering was weak (top hat filter of $4\Delta_{DNS}$, where $\Delta_{DNS}$ is the mesh resolution size of the DNS used for the a-priori analysis) compared to practical LES resolution at high Reynolds numbers. In a posteriori results, the model led to good predictions of mean quantities such as profiles of mean velocity and mean Reynolds stresses.

A similar model was defined in [42] with a slightly different formulation:

$$ \frac{\partial \tau_{ij}}{\partial x_j} = C_f \frac{\partial}{\partial x_j} (\delta u_i \delta u_j) \qquad (28) $$

where:

$$ \tilde{\delta}u_i(\boldsymbol{x}, \Delta)\delta u_j(\boldsymbol{x}, \Delta) = \frac{1}{2} [ (\bar{u}_i(\boldsymbol{x} + \Delta e_i) - \bar{u}_i(\boldsymbol{x})) (\bar{u}_j(\boldsymbol{x} + \Delta e_j) - \bar{u}_j(\boldsymbol{x})) \\ + (\bar{u}_i(\boldsymbol{x}) - \bar{u}_i(\boldsymbol{x} - \Delta e_i)) (\bar{u}_j(\boldsymbol{x}) - \bar{u}_j(\boldsymbol{x} - \Delta e_j))] \qquad (29) $$

The model coefficient $C_f$ was derived for the model (28) under several assumptions but a dynamic procedure was also tested. The model was also tested a posteriori in a channel flow and the mean velocity, mean turbulent kinetic energy and mean Reynolds stresses were in agreement with DNS predictions. Moreover, a priori analysis of the non-dimensionalized subfilter dissipation tensor components evaluated with fluctuating velocity, shows that the model is able to predict mean backscatter in some specific directions and regions of the flow where it is supposed to exist according to the DNS results from [43]. This is a promising result as classical dissipative models are not able to capture this mean backscatter.

More recently, Cimarelli et al introduced the mixed model [14]:

$$ \tau_{ij} = C_{\Delta} \tilde{\delta} u_i \delta u_j - 2\nu_T \overline{S_{ij}} \qquad (30) $$

This model was inspired by the increment model in [2] and the small-scale KHMH energy equation. The difference with the increment model is the combination with the Smagorinsky model in (30). The exact subfilter stress formulation (23) which involves $\delta u_i \delta u_j$ was also used to justify the first term on the right hand side of the model (30) but only as an inspiration as the model is not derived explicitly from (23). The increment model (26) is expected to capture backscatter events and the Smagorinsky model is introduced to model most of the dissipative effects of subfilter scales. A dynamic procedure was set to compute $C_{\Delta}$, and Germano's dynamic procedure was used for the Smagorinsky constant. Very good results were obtained for the averaged subfilter dissipation in a channel flow compared to filtered DNS results. Again, the regions of the flow containing mean backscatter are found to be well predicted by the model. The Smagorinsky model added to the modified model of Brun et al [2] seems to be important to capture the right energy transfer. The proposed model gives better results than classical mixed models such as the scale similarity model used in combination with the Smagorinsky model [18]. This is especially true for backscatter predictions. The model (30) proposed by Cimarelli et al is also close to but not exactly the same as the gradient model [19] introduced in section I.

Subfilter stress models based on $\delta u_i \delta u_j$ seem to be promising for large eddy simulation because they can capture backscatter events both statistically (probability distribution function) and on average. The models in [2], [42] and [14] are also very interesting as they are close to the exact equation (24). However, the physical and mathematical links with two-point equations can be developed further to support physically the model and bring all the physical understanding of the two-point equations into LES modelling. Such analysis requires to use the large scale two-point equations introduced in [1] instead of the small-scale KHMH equation. The models (26), (28) and (30) have only been tested extensively in channel flows and without local analysis of the subfilter stress. Therefore, these models may have some unknown limitations that will appear when testing them in other flows or by evaluating them locally.
Moreover, if the models are not adapted to the small-scale turbulence physics they will have to be adjusted to each flow limiting their universality. The main goal of this study is to try to develop a LES model from the fundamental equations and well-defined approximations in order to obtain a more universal model based on physical findings as opposed to models that need to be adapted for each flow and each specific simulation.

### C. Simplification of the exact subfilter stress formulation

In the present section, the Germano exact subfilter stress gradient formulation (24) is used to derive a new subfilter stress model. Usually, the subfilter stress is modelled rather than its gradient. This is the case for turbulent viscosity models which are based on the intuitive hypothesis that the effect of small scales is analogous to a molecular mechanism represented by a diffusion term. However, the new association of the subfilter contribution with the interscale transfer rate suggests to model directly the gradient. The change of variable $s \leftarrow x - s$ is used in equation (24) leading to

$$ \frac{\partial \tau_{ij}}{\partial x_j} = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(s-r)G(s+r)\nabla_r \cdot (u_{\mathbf{X}_i}(\mathbf{x}-s, \mathbf{r})\delta\mathbf{u}(\mathbf{x}-s, \mathbf{r}))drds. \quad (31) $$

The three-dimensional top hat filter of filter width $\Delta$ is used to simplify computations:

$$ G: \mathbb{R}^3 \rightarrow \mathbb{R} \\ \mathbf{x} \mapsto \begin{cases} 1: & \text{if } |\mathbf{x}|_{\infty} \le \Delta \text{ with } |\mathbf{x}|_{\infty} = \max\{|x_1|, |x_2|, |x_3|\} \\ 0: & \text{otherwise} \end{cases} \quad (32) $$

Our approach is restricted for now to this filter which is known to have a weak scale-selectivity. The derivation may be extended in the future to other filters such as the Gaussian filter. With this filter (32), (31) can be re-written as:

$$ \frac{\partial \tau_{ij}}{\partial x_j} = \frac{1}{V_{ref}} \int_{V_{ref}} \nabla_r \cdot (u_{\mathbf{X}_i}(\mathbf{x}-s, \mathbf{r})\delta\mathbf{u}(\mathbf{x}-s, \mathbf{r})) dr ds \quad (33) $$

where the integration is performed over a domain $V_{ref}$ of volume $V_{ref}$ defined as:

$$ V_{ref}: \begin{cases} -\Delta \le s_i - r_i \le \Delta \\ -\Delta \le s_i + r_i \le \Delta \end{cases} \quad (34) $$

One can define 6-dimensional volume regions $V_1$ and $V_2$

$$ V_1: \begin{cases} -\frac{\Delta}{2} \le s_i \le \frac{\Delta}{2} \\ -\frac{\Delta}{2} \le r_i \le \frac{\Delta}{2} \end{cases} \quad (35) $$

$$ V_2: \begin{cases} -\Delta \le s_i \le \Delta \\ -\Delta \le r_i \le \Delta \end{cases} \quad (36) $$

so that $V_1 \in V_{ref} \in V_2$ (see figure 3). In the following, $V_{ref}$ is approximated by $V_2$ for simplification. This approximation leads to an over-filtering compared to the initial top hat filter. Using this approximation, (33) becomes:

$$ \frac{\partial \tau_{ij}}{\partial x_j} \approx \frac{1}{V_2} \int_{V(\Delta)} \int_{V(\Delta)} \nabla_r \cdot (u_{\mathbf{X}_i}(\mathbf{x}-s, \mathbf{r})\delta\mathbf{u}(\mathbf{x}-s, \mathbf{r})) d^3r d^3s \equiv \Gamma_i \quad (37) $$

where $V_2 = \int_{V(\Delta)} \int_{V(\Delta)} d^3r d^3s$ and $\int_{V(\Delta)}(\cdot)d^3x = \int_{-\Delta}^{\Delta} \int_{-\Delta}^{\Delta} \int_{-\Delta}^{\Delta}(\cdot)d^3x$. The purpose of this approximation is the decoupling of the space filtering (s space) and the scale filtering (r space):

$$ \Gamma_i = \frac{1}{(2\Delta)^3} \int_{V(\Delta)} \left( \frac{1}{(2\Delta)^3} \int_{V(\Delta)} \nabla_r \cdot (u_{\mathbf{X}_i}(\mathbf{x}-s, \mathbf{r})\delta\mathbf{u}(\mathbf{x}-s, \mathbf{r})) d^3r \right) d^3s. \quad (38) $$
<figure><img src="image_5.png" /><figcaption>FIG. 3: Volume representation in one dimension of the volume $V_1$, $V_2$ and $V_{ref}$ defined respectively in equations (34),(35) and (36).</figcaption></figure>

The effect of this first approximation (37) is analysed and quantified in section IV B). The decoupling between the $x$ space and the $r$ space is useful because the divergence theorem can now be used in $r$ space:

$$ \Gamma_i = \frac{1}{(2\Delta)^3} \int_{V(\Delta)} \frac{1}{(2\Delta)^3} \left[ \oint_{|r|_{\infty}=\Delta} (u_{\mathbf{X}i}(\mathbf{x}-\mathbf{s}, \mathbf{r})\delta\mathbf{u}(\mathbf{x}-\mathbf{s}, \mathbf{r})) \cdot \mathbf{n} \, d^3\mathbf{r} \right] d^3\mathbf{s} \quad (39) $$

where $\mathbf{n}$ is the outward pointing unit normal vector at each point on the closed surface of integration. In this formulation, the closed surface integration is evaluated with the component $r_j$ in the direction $n_j$ such that $r_j = \Delta$. This is interesting as only quantities with a separation distance of $\Delta$ or larger (in multiples of $\Delta$) can be evaluated on the LES resolution mesh. However, the values of the other components of $\mathbf{r}$, $r_i \le \Delta$ where $i \ne j$ ($n_j$ is the direction of projection) are still needed during the surface integration in (39) but are neglected later in the derivation. The integration in physical space ($d^3\mathbf{s}$) needs to be replaced because it requires values not available during LES computations. The filtered two-point velocity half sum and two-point velocity increment using the three-dimensional top hat filter are:

$$ \begin{aligned} \overline{u}_{\mathbf{X}i}(\mathbf{x}, \mathbf{r}) &= \frac{1}{(2\Delta)^3} \int_{V(\Delta)} u_{\mathbf{X}i}(\mathbf{x}-\mathbf{s}, \mathbf{r}) d^3\mathbf{s} = \overline{u}_{\mathbf{X}i}(\mathbf{x}, \mathbf{r}) \\ \overline{\delta u}(\mathbf{x}, \mathbf{r}) &= \frac{1}{(2\Delta)^3} \int_{V(\Delta)} \delta u(\mathbf{x}-\mathbf{s}, \mathbf{r}) d^3\mathbf{s} = \overline{\delta u}(\mathbf{x}, \mathbf{r}) \end{aligned} \quad (40) $$

which suggests to rewrite equation (39) as follows:

$$ \Gamma_i = \frac{1}{(2\Delta)^3} \oint_{|r|_{\infty}=\Delta} (\overline{u}_{\mathbf{X}i}(\mathbf{x}, \mathbf{r})\delta\overline{u}(\mathbf{x}, \mathbf{r})).\mathbf{n}d^2\mathbf{r} \\ + \frac{1}{(2\Delta)^3} \left[ \oint_{|r|_{\infty}=\Delta} \left( \frac{1}{(2\Delta)^3} \int_{V(\Delta)} u_{\mathbf{X}i}(\mathbf{x}-\mathbf{s}, \mathbf{r})\delta\mathbf{u}(\mathbf{x}-\mathbf{s}, \mathbf{r})d^3\mathbf{s} - \overline{u}_{\mathbf{X}i}(\mathbf{x}, \mathbf{r})\delta\overline{u}(\mathbf{x}, \mathbf{r}) \right).\mathbf{n}d^2\mathbf{r} \right]. \quad (41) $$

Equation (41) can also be written as:

$$ \Gamma_i = \frac{1}{(2\Delta)^3} \oint_{|r|_{\infty}=\Delta} (\overline{u}_{\mathbf{X}i}(\mathbf{x}, \mathbf{r})\delta\overline{u}(\mathbf{x}, \mathbf{r})).\mathbf{n}d^2\mathbf{r} \\ + \frac{1}{(2\Delta)^3} \left[ \oint_{|r|_{\infty}=\Delta} \left( \frac{1}{(2\Delta)^3} \int_{V(\Delta)} u'_{\mathbf{X}i}(\mathbf{x}-\mathbf{s}, \mathbf{r})\delta\mathbf{u}'(\mathbf{x}-\mathbf{s}, \mathbf{r})d^3\mathbf{s} \right).\mathbf{n}d^2\mathbf{r} \right] \quad (42) $$
where $u'_{X_i}(x-s,r) = u_{X_i}(x-s,r) - \bar{u}_{X_i}(x,r)$ and $\delta u'(x-s,r) = \delta u(x-s,r) - \bar{\delta u}(x,r)$. The first part is the resolved part because it is computed with the resolved velocity field and the second part is the unresolved part. In section IV B, this second unresolved part is found to be much smaller than the first resolved part. However it should not be neglected as it contains important turbulent behavior at small scales and is expected to be related to dissipation. Indeed, the resolved part depends only on filtered quantities similarly to the Leonard term in [44] (though being mathematically different to the Leonard term). Therefore, this term is not expected to contain the majority of the dissipation so that it has to be present in the unresolved term. Moreover, the unresolved term depends only on residual velocity (without dependency on filtered velocity field) which again supports that this term is related to small scale statistics and should therefore contain most of the dissipation. Hereafter, we refer to this unresolved term as the residual term.

### D. Formulation of a new subfilter stress model

The resolved part of the model is expected to capture the local behavior of the energy transfer with large fluctuations mainly driven by the resolved scales but it probably does not capture the average dissipation associated with small scales. As already mentioned, the dissipation is likely to be related to the residual term of equation (42). This interpretation may only be valid for filtering width $\Delta \gg \eta$ so that the dissipative scales are completely unresolved. Under such hypothesis the residual term can be modelled with existing dissipative models such as the Smagorinsky model or with a spectral vanishing viscosity model which can be implemented in physical space as in [8]. In this study, the Smagorinsky model is used to model dissipation for simplicity. Some limitations of this model are well known but the dissipative model is not the main object of the present study. The conclusions on the resolved part are not expected to be significantly dependent on the choice of dissipation model. Using the Smagorinsky model for the unresolved part of (42) in a mixed configuration leads to the following model:

$$ \Gamma_i = \frac{C_1}{(2\Delta)^3} \oint_{|r|_{\infty}=\Delta} (\bar{u}_{X_i}(x,r)\delta\bar{u}(x,r)).nd^2r + \frac{\partial(-2\nu_T\bar{S}_{ij})}{\partial x_j} \quad (43) $$

where $C_1$ is a dimensionless coefficient which may depend on $\Delta$ and is introduced to rescale the resolved part of the model if needed. The value of both $C_1$ and the Smagorinsky constant $C_S$ in $\nu_T$ is discussed later in this paper. The first part refered as the resolved term is Galilean invariant (see appendix A) and the second part introduced to compensate the residual term in (42) is refered as the dissipative term. The resolved term in (43) is written with a closed surface integral over a mesh cell. However, in practice, this integral can be computed only with the available grid points. The simplest option is therefore to approximate the surface integral of the resolved term as

$$ \frac{C_1}{(2\Delta)^3} \oint_{|r|_{\infty}=\Delta} (\bar{u}_{X_i}(x,r)\delta\bar{u}(x,r)).nd^2r \\ \approx \frac{C_1}{(2\Delta)^3} \frac{6(2\Delta)^2}{6 \times 9} \left[ \begin{aligned} & \sum_{r_y \in I_\Delta} \sum_{r_z \in I_\Delta} \bar{u}_{X_i}(x, r_x = \Delta, r_y, r_z) \delta\bar{u}_x(x, r_x = \Delta, r_y, r_z) \\ & - \sum_{r_y \in I_\Delta} \sum_{r_z \in I_\Delta} \bar{u}_{X_i}(x, r_x = -\Delta, r_y, r_z) \delta\bar{u}_x(x, r_x = -\Delta, r_y, r_z) \\ & + \sum_{r_x \in I_\Delta} \sum_{r_z \in I_\Delta} \bar{u}_{X_i}(x, r_x, r_y = \Delta, r_z) \delta\bar{u}_y(x, r_x, r_y = \Delta, r_z) \\ & - \sum_{r_x \in I_\Delta} \sum_{r_z \in I_\Delta} \bar{u}_{X_i}(x, r_x, r_y = -\Delta, r_z) \delta\bar{u}_y(x, r_x, r_y = -\Delta, r_z) \\ & + \sum_{r_x \in I_\Delta} \sum_{r_y \in I_\Delta} \bar{u}_{X_i}(x, r_x, r_y, r_z = \Delta) \delta\bar{u}_z(x, r_x, r_y, r_z = \Delta) \\ & - \sum_{r_x \in I_\Delta} \sum_{r_y \in I_\Delta} \bar{u}_{X_i}(x, r_x, r_y, r_z = -\Delta) \delta\bar{u}_z(x, r_x, r_y, r_z = -\Delta) \end{aligned} \right] \quad (44) $$

where $I_\Delta = \{-\Delta, 0, \Delta\}$. The available points used for the approximation of the surface integral are plotted in figure 4 where the black arrows represent the projections of $\bar{u}_{X_i}(x,r)\delta\bar{u}(x,r)$ used during surface integration. As opposed to the previous increment models (27) and (29), this derivation leads to a natural definition of $\delta u$ with a prescribed separation vector $r$ related to $\Delta$ in all directions.
<figure><img src="image_6.png" /><figcaption>FIG. 4: Representation of the points available in r space for the surface integration in (44). Each arrow describes a point used during surface integration (with a coefficient 1) with the projection direction.</figcaption></figure>

<caption>TABLE I: Parameter for the reference DNS of forced Taylor Green flows</caption>

$$ \frac{\text{Number of points} \quad \text{Size of the domain} \quad Re_{\lambda} \quad k_{\text{trunck}} \eta}{512 \times 512 \times 512 \quad 7.1L \times 7.1L \times 7.1L \quad 91 \quad 3} $$

The error introduced by replacing equation (39) with equation (44) is quantified in section IV B). An interpolation for a better estimation of the integral in (44) might improve the results and could be used in future work. Also, for simplicity, the mesh is considered here to be isotropic but the results can in principle be generalized if necessary.

# IV. EVALUATION OF THE MODEL

## A. Reference DNS

An in-house pseudo-spectral code is used to solve the Navier-Stokes equations with periodic boundary conditions in all three dimensions. The model is introduced in physical space with the non-linear term of Navier Stokes. A de-aliasing is applied using the 2/3 truncation method (truncation for wavenumbers $k$ such as $|k| > N/3$, where N is the number of grid point in each dimension) and a second order Runge-Kutta scheme is used for the time integration. More information about the original code can be found in [45]. A Taylor-Green flow (where the forcing is described in [46] and [8]) is simulated with a 512³ grid until statistical convergence. The simulation parameters are described in table I where $L = \frac{3\pi}{4} \frac{E(k)}{dk}$ is the integral scale. The simulation is highly resolved (the mesh resolution, $\Delta x = \frac{2\pi}{N} \approx 0.7\eta$) and the Reynolds number is $Re_{\lambda} \approx 91$.

## B. A priori validation of the results

Based on a priori analysis, an accurate subfilter scale model should predict the occurrence of dissipation and backscatter at the right location in space and time and should generate the right level of energy transfer (dissipation or backscatter). These two properties are evaluated separately. The ability of the model to predict the correct sign and the locations of the energy transfers between residual and filtered scales is evaluated first. The quantitative analysis about the energy transfer which is linked to the definition of the model coefficients is carried out in a second stage.
<figure><img src="image_7.png" /><figcaption>FIG. 5: PDF of the subfilter stress contribution evaluated *a priori* from the DNS with a 12³ top hat filter: exact (blue), Smagorinsky model (red), increment model (26) (yellow), scale similarity model (6) (green), and for the resolved part of the present model (43) (purple).</figcaption></figure>

1. *Qualitative comparison*

The resolved part (44) of the new mixed model (43) is first evaluated *a priori* on the DNS velocity field and compared to the Smagorinsky model (4) with $C_S = 0.18$ and the model (26) of Brun et al. A 12³ points cubic top-hat filter corresponding to a volume of $\Delta^3 = (8\eta)^3$ is used to separate the large filtered scales from the subfilter scales. The full energy transfer between the residual scales and the filtered scale $T^e$, introduced in (2), is evaluated as function of space and time. The negative values of $T^e$ reduce the energy of $\bar{u}_i\bar{u}_i$ so they can be seen as dissipative while positive values represent backscatter. The PDF of $T^e$ is presented in figure 5. The choice of $\Delta = 8\eta$ is quite close to the Taylor length $\lambda = 10.7\eta$ which demarcates the inertial from the dissipative ranges (e.g. see [34] and [47]) so that most (though not all) dissipative scales are unresolved. The subfilter scales are mainly in the dissipative range and the energy transfer is to a significant extent between the inertial and the dissipative ranges. In fact, in [34] and [47], the non-linear interscale energy transfer defined in section II, averaged over all directions in scale space, is found to behave at scales $r$ smaller than but close to $\lambda$ similarly to scales in the inertial range with large positive and negative events.

For this first test, the model constants $C_1$ for our model (43), $C_2$ for the model (26) and $C_3$ for scale similarity model (6) are optimized to obtain the best match of the PDF of energy transfer $T^e$ with the reference PDF obtained *a-priori* by the DNS. The mean $T^e$ which can be computed from these PDFs corresponds to the mean contribution of the subfilter stress (expected to be negative) and the tails correspond to the rare, very intense events.

The PDF of energy transfer for the Smagorinsky model (with $Cs = 0.176$ obtained after numerical optimization similarly to the other models) has a significantly different shape compared to the reference PDF from DNS with a skewness that is too negative. This result is consistent with the well known over-dissipative behavior of the model. Backscatter events are scarce and their magnitude are relatively small. The PDF for the model (26) of Brun et al has a more realistic shape with both dissipative and backscatter events but it doesn't fit very well the reference results from DNS even after optimization of the model coefficient $C_2$. Both the subfilter stress contribution of the scale similarity model and the resolved part of our model are closer to the reference once the associated coefficients $C_1$ and $C_3$ are properly adjusted. These *a-priori* results seem to indicate a similarity in behavior between our model directly approximated from exact equations and related to the interscale turbulent energy transfer and the scale similarity model, which is commonly used in a mixed configuration for LES. This similarity is encouraging as it suggests a possible practical use of our new model.

The results for the resolved part with the coefficient $C_1$ used in figure 5 combined with the Smagorinsky model in a mixed configuration (43) and the scale similarity model also used in a mixed configuration as in (6) are presented in figure 6. Both models are over-dissipative with skewnesses that are too negative compared to exact results. However, it has been shown that the resolved part of our model does capture the large range of variation of the energy transfer rate and contains backscatter. Therefore, this model can in principle be used in Large Eddy Simulation but there is a need for future work to design a new dissipation model which would better complement the resolved part.
<figure><img src="image_8.png" /><figcaption>FIG. 6: PDF of the subfilter stress contribution evaluated *a priori* from the DNS with a 12³ top hat filter: exact (blue), scale similarity model in a mixed configuration (7) (green) and the present model (43) including its Smagorinsky model part (purple).</figcaption></figure>

The energy transfers $T^e$ are now analyzed in space for the same models. The values of $T^e$ in a randomly selected plane are compared in Fig. 7 with the reference for a first visual inspection of the spatial correlation between exact and modeled values. The same adjusted constants as for the PDFs in figure 5 and 6 are used. In figure 7a, clear structures are identified in the exact energy transfer function with well identified large events and filaments in the background. Large dissipative events are often observed to be close to large backscatter events which is consistent with the theoretical analysis of interscale flux energy transfer in [30]. As can be seen in figure 7b, the backscatter events are not captured by the Smagorinsky model. Some correlation with the reference seems to exist for the large events and the filaments compared to exact ones but the sign is often wrong and large backscatter events are even sometimes replaced by large dissipative events.

The results of the increment model (26) in figure 7c contain both dissipative and backscatter events as wanted. The background filaments are comparable to the reference but the large events do not correlate well with the reference except at specific locations. A better correlation with the reference is observed with the resolved part (44) of our model (43) where some of the large events predicted are representative of real events (see figure 7e). The background filaments are also well predicted even if the magnitude of some energy transfer events or their exact locations may be incorrect. The best results are obtained with the scale similarity model (6) in figure 7d where the large events modeled are very close to the real events. These results suggest that our approximations to derive the resolved part of our model can be improved further but the results are already close to the well known scale similarity model while having an exact connection to the exact equations (9) and (11).

As the previous analysis is restricted to a single snapshot in space, the correlation over the full domain is now analyzed for a more reliable conclusion. The joint PDFs of the modelled subfilter stress contribution for the same models are presented in figures 8(a-d). This test is more sensitive than a simple visualizations of a spatial map because a small position uncertainty of a significative energy transfer event is clearly identified as a strong decorrelation.

The Smagorinsky model has nearly no local correlation with the exact values which confirms the impressions given in figures 7a. The increment model (26) in figure 8b has a small correlation but the values around the axes suggest also strong anti-correlation. These results are probably due to the small position shift of the strong energy transfer structures as a correlation seems to be observed when looking at the map in figure 7c. The resolved part of our model (44), with $C_1$ optimized, has a better correlation compared to the other models (see figure 8d). The results are not fully satisfactory and can probably be improved in future work. The scale similarity model has the best results as observed in figure 8c. These results are consistent with the analysis of the local results in figure 7 and suggests the possibility of further refinement of the approximations used for the derivation of our model. Our approach being based on a progressive simplification of Germano's exact subfilter stress equation, it is possible to investigate the origin of the errors by checking the main steps where approximations are introduced.

For the present *a priori* analysis, a filter width $\Delta \approx 8\eta$ is used in all the steps of the validation process. The evaluation of the error introduced when approximating (24) requires a lot of memory and computational resources. Indeed, the size of $\Delta$ is the limiting parameter of this analysis because it is involved in the complex Germano filtering used in (24). This is the reason why the full analysis presented in the following section is conducted with a single
<figure><img src="image_9.png" /><figcaption>(a) Exact from DNS</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Smagorinsky model</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Increment model [2]</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Scale similarity model</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(e) Present model (resolved part)</figcaption></figure>

<caption>FIG. 7: 2D slice of $T^e = -\overline{u_i \frac{\partial \tau_{ij}}{\partial x_j}}$ evaluated *a-priori* from the DNS dataset: (a) exact, (b) Smagorinsky model, (c) increment model (26), (d) scale similarity model (6) (e) present model (with the resolved part only).</caption>
<figure><img src="image_14.png" /><figcaption>FIG. 8: Joint PDFs of $T^e = -\overline{u_i \frac{\partial \tau_{ij}}{\partial x_j}}$ (horizontal axis) with the modelled counterpart (vertical axis) evaluated a priori from the DNS with $\Delta \approx 8\eta$ and with scaling coefficient optimized for a best match of the PDF of energy transfer in Fig. 5</figcaption></figure>

filter width of size $\Delta \approx 8\eta$. However, the *a priori* analysis of the modelled subfilter stress contribution computed in this section was also done for $\Delta \approx 16\eta$ and the conclusions remain similar. The filter width remains relatively small so that the assumption that the dissipative scales are completely unresolved is not respected during the *a priori* validation. This is why a significant skewness is observed in the PDF of the energy transfer results in figures 5 and 6. Indeed, some of the dissipative scales affect the energy transfer of the resolved part of the model. This filter width limitation should not, however, affect the large tails of the PDF which are the main purpose of the resolved part of our model and should remain with a larger filter width.

## 2. Errors introduced with the different approximations

As our model is directly derived from the Navier Stokes equations, the errors introduced from the exact subfilter stress formula (23) to the approximated formulation (43) of our model can be analysed at each approximation step. The quantification of these approximation errors are important to better understand the performance of the model and to have the possibility to propose future improvements.

The analysis is only conducted with the first component of the subfilter stress gradient for simplicity. The joint PDF of the exact subfilter stress component ($\frac{\partial \tau_{1j}}{\partial X_j}$) with the corresponding approximated values at 3 steps of the
<figure><img src="image_15.png" /></figure>

(a) Exact formulation from (24):
$$ \iiint_{-\infty}^{\infty} \iiint_{-\infty}^{\infty} G(s-r)G(s+r)\nabla_r \cdot (u_{X1}\delta u) d^3r d^3s. $$

(b) After dissociation of scale and space integration as in (39):
$$ \iiint_{-\Delta}^{\Delta} \mathfrak{F}_{|r|=\Delta} u_{X1}\delta u \cdot nd^2rd^3s. $$

<figure><img src="image_16.png" /></figure>

(c) Real numerical implementation as in (44)

<caption>FIG. 9: Two-dimensional probability distribution function of $\frac{\partial \tau_{1j}}{\partial X_j}$ with its different approximations evaluated in the Taylor-Green DNS with $\Delta \approx 8\eta$. (a) Exact formulation from equation (24):</caption>

$\iiint_{-\infty}^{\infty} \iiint_{-\infty}^{\infty} G(s-r)G(s+r)\nabla_r \cdot (u_{X1}\delta u) d^3r d^3s$, (b) after dissociation of scale and space integration as in (39):

$\iiint_{-\Delta}^{\Delta} \mathfrak{F}_{|r|=\Delta} u_{X1}\delta u \cdot nd^2rd^3s$, (c) real numerical implementation of (44).

derivation corresponding to equations (31), (37) and (44) are evaluated a priori on the DNS dataset. Note that equations (31) and (37) must be discretized to be evaluated over our DNS dataset as opposed to equation (44) which is already discrete. The numerical implementation of the computation of the continuous terms is performed on the DNS grid without interpolation and discrete centered scheme are used to compute derivatives. As the DNS is highly resolved ($\Delta x = 0.7\eta$), the error introduced during the discretization of (31) and (37) are negligible. The results for the first step in figure 9a are presented to check the error introduced by the numerical implementation of the exact equation (31). The very good correlation demonstrates the low impact of this numerical implementation.

The error introduced when decoupling space and scale coordinates in equation (37) are presented in figure 9b. The elliptical shape indicates a high level of correlation meaning that the space-scale decoupling does not introduce significant error in the present case at moderate Reynolds number with a filter width $\Delta \approx 8\eta$. The results in figure 9c demonstrate that most of the numerical error is introduced when moving from equation (37) to equation (44) where the equation is discretized over the LES mesh and where the residual term is neglected. The error introduced by this approximation might be reduced with a better numerical implementation of the resolved term, using interpolation for example.
### C. *a posteriori* evaluation of the models

Large Eddy Simulations of a Taylor-Green flow are carried out with a same pseudo-spectral code as the DNS and with the same forcing. A 64 × 64 × 64 discretization is used (8 times coarser than the DNS) and the CFL is set to 0.5. The 12 × 12 × 12 top hat filter used for the *a priori* analysis corresponds to a -3dB cutoff wavenumber $k_{cutoff} \approx 19$. The resolution of the LES was chosen such that the truncation wavenumber ($\frac{2}{3}k_{max} \approx 21$) is close to this value which makes it easier to compare results. However, as the sharpness of the top hat and spectral cutoff filter are different, the comparison between the *a-priori* and a posteriori evaluation can only be partial. This may explain why the optimized model constants need to be adapted for each approach.

The LES are initialized with the DNS velocity field sub-sampled on a 64 × 64 × 64 mesh. The *a posteriori* analysis is conducted with the following models:

* The Smagorinsky model defined in equation (4) where $C_S = 0.18$. A second order central difference scheme is used for the derivatives so that the model is fully implemented in physical space.

* The model (26) implemented the same way as in [2] where the constant $C_b = 0.152$ is evaluated *a priori*. This model is tested alone and in combination with the Smagorinsky model with $C_S = 0.18$.

* The scale similarity model used in a mixed configuration (see equation (7)) with $C_3 = 1$. As the value of $C_3$, evaluated during the *a priori* analysis was leading to an over-filtered energy spectrum in the *a posteriori* simulation, different values of $C_3$ were tested based on both the energy spectrum and the distribution of energy transfer between resolved and residual scales (see figures 10b and 11a). The value of $C_3$ affects significantly the results but no good compromise could be found. However, the results with $C_3 = 1$ presented below are representative of the different simulations. For these results a top hat filter over 15 points is used in (7) for the explicit filter identified with a tilde symbol. The scale similarity model is tested alone and in association with the Smagorinsky model with $C_S = 0.18$.

* Our model (43) derived from the two point equations with the resolved part alone and associated with the Smagorinsky model in a mixed configuration. The value of the constant $C_1$ evaluated during the *a priori* analysis was leading to over-filtered energy spectrum in the *a posteriori* simulation. Therefore, the value of $C_1$ was optimized to have an energy spectrum close to the one of the simulation with the Smagorinsky model alone and the right energy transfer between filtered scales and residual scales. As a results of this optimization a constant $C_1 \approx 0.39$ was selected for the present simulation. In the mixed configuration, the Smagorinsky coefficient $C_S$ is adjusted to have the same amount of energy transfer $T^e = u_i \frac{\partial \tau_{ij}}{\partial x_j}$ as for the simulation with the Smagorinsky model alone. This method was also tested for the scale similarity model without success.

In a first stage, our model (43) with the resolved part only is evaluated by comparing the energy spectrum averaged over 24 turnover times with the two other models (the increment model, and the scale similarity model) and the reference DNS. The results in figure 10a suggest that the increment model (26), the scale similarity model (6) and the resolved part of our model do not introduce enough dissipation as evidenced by the level of energy at high wavenumbers ($k/k_\eta > 0.2$ with $k_\eta = \pi/\eta$) which is too high and close to the simulation without model. This confirms the suggestion of Cimarelli et al [14] to use a version of model (26) in combination with a dissipation model and our strategy to combine the resolved part of our model with a dissipation model as anticipated with the *a-priori* tests. However, some small scale dissipation is nevertheless introduced by the resolved part of our model as the energy at wavenumbers $k/k_\eta > 0.2$ is slightly smaller than for the simulation without model. This suggests that the results of the resolved model associated with the Smagorinsky model would be over-dissipative compared to the Smagorinsky model alone if we do not change the constant $C_S$. This result may come from the numerical resolution, which is chosen accordingly to the filtering width used in the *a priori* analysis. Indeed, the resolved quantities certainly contain some of the dissipative scales due to the relatively high numerical resolution as opposed to the assumptions of our model derivation.

In a second stage, all the models are tested in mixed configurations (combined with a Smakorinsky model) and compared to the DNS and to the coarse DNS (LES without model). The energy spectra are given in Fig. 10b. The results are highly dependent on the coefficients used for the different models. As the coefficients evaluated during the *a priori* analysis do not lead to good results, these coefficients are optimized to match the energy spectrum of the simulation with Smagorinsky model alone and the right distribution of energy transfer. Indeed, a small variation affects significantly the results so the analysis is more qualitative than quantitative. All the LES are over-dissipative compared to the filtered DNS results as observed with the excessive drop of energy at high wavenumbers ($k/k_\eta > 0.2$).
<figure><img src="image_17.png" /><figcaption>FIG. 10: Time averaged energy spectrum of large eddy simulation for several models with $k_{\eta} = \pi/\eta$ and $\eta/\Delta \approx 0.13$ are compared with the DNS. In (a), the model of Brun et al (eq. (26), purple), the scale similarity model (eq. (6), green) and the present model (eq. (43), light blue) are used without additional turbulent viscosity model and in (b) the same models are used in a mixed configuration with the Smagorinsky model.</figcaption></figure>

for the simulations with model in figure 10b. The shape of the energy spectra seems to be mainly driven by the eddy viscosity model as they are similar for the four simulations with models. The scale similarity model in a mixed configuration is more dissipative as compared to the simulation with the Smagorinsky model alone. For the scale similarity model, the use of a dynamic coefficient $C_s$, adjusted to have the same amount of energy transfer $T^e$ as for the simulation with the Smagorinsky model alone, was tested but did not improve the results. Therefore no dynamic method was used to evaluate $C_s$ for this model. However, this does not mean that the results cannot be improved with methods such as dynamics methods in [15] or [17] for the evaluation of the coefficient.

The energy spectrum is an average statistic which does not reflect the space-local properties of the model. This might explain why taking into account the local backscatter effect is not clearly visible in the *a posteriori* results. Indeed, the non-homogeneous contributions to the energy equation vanish in Fourier space. Therefore, the impact of the resolved model, which introduces some local non-homogeneity contributions, may only be visible in stationary non-homogeneous turbulence or in homogeneous turbulence but only locally. Thus, the PDFs of the energy transfer $T^e$ from resolved to subfilter scales are now analyzed for the different models to describe the local properties of the simulated turbulence.

These probability distribution functions of $T^e$ sampled in space and time over 24 turnover times are presented with linear axes in figure 11a and with logarithmic axes in figure 11b. As observed with the highly negatively skewed PDF with small tails, the Smagorinsky model alone produces mainly dissipation similarly to the previous *a priori* results meaning that energy is extracted from resolved scales and that the rare energetic transfer events are not well predicted. The increment model (26) and our model combined with the Smagorinsky model improve significantly the results by introducing rare energy transfer events (large backscatter and dissipative events) as observed with the near superposition of the PDFs of $T^e$ for these simulations and the exact results computed from DNS. The PDF for the scale similarity model in a mixed configuration (6) is slightly improved as compared to the one with Smagorinsky model alone but remains far from the DNS. As the coefficient $C_3$ in the scale similarity model (6) affects the results, larger values were also tested. With a larger coefficient, the tails of the PDF widen but their shape remains comparable to the results presented with $C_3 = 1$, while affecting significantly the energy spectrum described previously. Therefore, contrary to the conclusions on energy spectra where the results do not depend significantly on the model (whether model (26) of [2], scale similarity model (6) or the resolved part of our model), the space-local behaviour of the different simulations are significantly affected by these models. Indeed, the physics of the energy transfer, between filtered and residual scales, happening locally in the simulations with the increment model (26) combined with Smagorinsky or with our model is different to the physics present in the simulation with the Smagorinsky model alone. These models capture the local non-homogeneity of the flow with large energy transfers between filtered and residual scales
<figure><img src="image_18.png" /><figcaption>FIG. 11: PDF of $T = -\bar{u}_i \frac{\partial \tau_{ij}}{\partial x_j}$ during large eddy simulation with different models. In blue: a priori results from DNS, yellow: Smagorinsky model, purple: structure model in a mixed configuration, green: scale similarity model in a mixed configuration (7), light blue: our model in a mixed configuration.</figcaption></figure>

both positive and negative which do not seem to affect significantly global homogeneous quantities such as the energy spectrum in periodic/homogeneous turbulence. From this analysis, the increment model and our model, in mixed configurations, reproduce better the energy exchanged between filtered scales and residual scales compared to the simulations with Smagorinsky model alone and scale similarity model in a mixed configuration. The low performance of the scale similarity model in the *a posteriori* analysis is unexpected given the very good results obtained from the *a priori* analysis. These results may be improved with a dynamic method or by changing the explicit filter (·) in (7). Our model, which is derived with simple approximations from exact two-points equations has therefore some rather good performances even compared to existing LES model commonly used.

The increment model (26) is found to affect global quantities such as mean flow or velocity RMS in non-homogeneous flows as observed in channel flows in [2] and the other variants of this model in [42] and [14]. Therefore, our model is to be tested in non-homogeneous turbulence in future work to see if the physics of local energy transfers between resolved and filtered scales, captured both *a priori* and *a posteriori*, improve the predictions of global quantities in LES of strongly non-homogeneous turbulent flows.

V. CONCLUSION

The small-scale Kármán-Howarth-Monin-Hill two-point equation [28] and large scale Germano equation [1] are used to develop turbulence physical knowledge and understanding. This mathematically exact framework allows the analysis of the turbulent energy cascade in non-homogeneous turbulent flows.

A new LES model is proposed based on approximations of the exact subfilter stress formulation (24) obtained by Germano [22]. This model is designed and tested as a proof of concept that new LES models can be derived as approximations of the exact two-point Germano equation.

Knowledge about the small-scale interscale energy transfer rate [33] and the large-scale interscale energy transfer rate [1] is used to justify the expected behavior of the LES model. LES models should reproduce the large fluctuations, including both forward and backward energy transfers, of the small- and the large-scale interscale energy transfer rates. A mixed model is proposed with a resolved part able to capture these fluctuations combined with a conventional dissipation model. A priori analysis of the resolved part of our model is carried out in a forced Taylor-Green flow and compared with the a priori analysis of the Smagorinsky model and the increment model [2] used in a mixed configuration as in [14]. Rather good correlation with DNS of local energy transfers between filtered and residual scales is obtained with the resolved part of our model. Particularly, backscatter events are well predicted which is not the case with turbulent viscosity models. The correlation of the exact rate of energy exchange $T^e = -\bar{u}_i \frac{\partial \tau_{ij}}{\partial x_j}$ with the
modelled one is better with the resolved part of our model than with the model of [2] but not as good as for the scale similarity model. However, our model is encouraging because it is derived from very crude approximations applied to exact equations and returns rather good results. This paves the way for future improvements of our approximations.

In a posteriori simulations, the resolved part of our model under-predicts the mean dissipation. This result was anticipated during the formulation of the model and this is why the resolved part is combined with a turbulent viscosity model in a mixed model configuration. This mixed model captures well, a posteriori, the probability distribution function of the local energy transfers between filtered and residual scales, and particularly backscatter events. The energy transfer characteristics of our model are comparable to those of the increment model [2] used in a mixed configuration as in [14] and even better than those of the scale similarity model in a mixed configuration. However, our model's good prediction of local energy transfers is not enough to ensure proper spatially averaged quantities such as the energy spectrum which is shown to be mostly determined by the dissipative model.

Future detailed assessments of our model as it currently stands or with improvements in the approximations leading to it must cover several other flow configurations. Our model can be expected to deliver enhanced performance in non-homogeneous turbulence and especially in regions of the flow where mean backscatter is present such as regions with high mean flow gradients, for example in turbulent channel flow as in [14]. The model should also be tested at higher Reynolds numbers to assess the Reynolds dependence of the resolved part of our model.

ACKNOWLEDGMENTS

This work was directly supported by JCV's Chair of Excellence CoPreFlo unded by I-SITE-ULNE (grant number R-TALENT-19-001-VASSILICOS), MEL (grant number CONVENTION_219_ESR_06) and Region Hauts de France (grant number 20003862).

This project was provided with computing HPC and storage resources by GENCI (Grand Equipment National de Calcul Intensif) at IDRIS thanks to the grant 2022-021741 on the supercomputer Jean Zay.

Funded by the European Union (ERC, NoStaHo, 101054117). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

Appendix A: Galilean invariance

In this appendix section, the resolved part of the model is shown to be Galilean invariant. The resolved term has the same form as the exact term $(\delta u \cdot \nabla_r) u_X$ in equation (11) once integrated in scale space. Therefore, the invariance properties of the Navier-Stokes equation are expected to be also applicable for the resolved term. However, the Galilean invariance property of our model is derived for clarity. Note that the Smagorinsky model is Galilean invariant [48].

Introducing a new reference frame moving at a constant velocity $U$ compared to our previous reference system: $t' = t$, $x' = x - Ut$ and $u' = u - U$ the resolved part $\oint_{|r|_{\infty}=\Delta} (\bar{u}_{X_i}(x, r)\delta\bar{u}(x, r)).n d^2r$ of the model (42) must be evaluated in this new frame using

$$ \left\{ \begin{aligned} \delta\bar{u}(x', r', t') &= \delta\bar{u}(x - Ut, r, t) \\ \bar{u}_X(x', r', t') &= \bar{u}_X(x - Ut, r, t) - U \end{aligned} \right. \qquad (A1) $$

Using (A1) and the property that $\nabla_r \cdot \delta\bar{u} = 0$ so that the closed surface integral of $\delta\bar{u}$ vanishes, the resolved part of the model can be written as follows in the reference frame moving with velocity $U$ relative to the original frame.

$$ \begin{aligned} \oint_{|r|_{\infty}=\Delta} & (\bar{u}_{X_i}(x', r', t')) \delta\bar{u}(x', r', t') . n \, d^2 r' \\ &= \oint_{|r|_{\infty}=\Delta} (\bar{u}_{X_i}(x - Ut, r, t)) \delta\bar{u}(x - Ut, r, t) . n \, d^2 r - U_i \oint_{|r|_{\infty}=\Delta} \delta\bar{u}(x - Ut, r, t) . n \, d^2 r \quad (A2) \\ &= \oint_{|r|_{\infty}=\Delta} (\bar{u}_{X_i}(x - Ut, r, t)) \delta\bar{u}(x - Ut, r, t) . n \, d^2 r \end{aligned} $$
Therefore, the resolved model is Galilean invariant as expected.

[1] M. Germano, The elementary energy transfer between the two-point velocity mean and difference, Physics of Fluids 19 (2007).

[2] C. Brun, R. Friedrich, and C. B. Da Silva, A Non-Linear SGS Model Based On The Spatial Velocity Increment: Application to LES of fully developed pipe flow and round turbulent jet, Theoretical and Computational Fluid Dynamics 20 (2006).

[3] J. Bardina, J. H. Ferziger, and W. C. Reynolds, Improved subgrid models for large eddy simulation, AIAA paper 80, 1357 (1980).

[4] A. N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds number, C.R. Acad. Sci URSS 30, 301 (1941).

[5] L. Richardson, Weather prediction by numerical process, cambridge university press ed. (1922).

[6] W. C. Reynolds, The potential and limitations of direct and large eddy simulations, in Whither Turbulence? Turbulence at the Crossroads, Vol. 357, edited by J. L. Lumley (Springer Berlin Heidelberg, 1990) pp. 313–343.

[7] S. B. Pope, Turbulent Flows (Cambridge Univ. Press, 2000).

[8] T. Dairay, E. Lamballais, S. Laizet, and J. Vassilicos, Numerical dissipation vs. subgrid-scale modelling for large eddy simulation, Journal of Computational Physics 337, 252 (2017).

[9] A. Vela-Martín, Subgrid-scale models of isotropic turbulence need not produce energy backscatter, J. Fluid Mech. 937, A14 (2022).

[10] U. Piomelli, W. H. Cabot, P. Moin, and S. Lee, Subgrid-scale backscatter in turbulent and transitionnal flows, Phys. Fluids A 3, 1766 (1991).

[11] J. Smagorinsky, General circulation experiment with the primitive equations, part i: The basic experiment, Mon. Wea. Rev. 91, 99 (1963).

[12] R. Vicente Cruz and E. Lamballais, Physical/numerical duality of explicit/implicit subgrid-scale modelling, Journal of Turbulence 24, 235 (2023).

[13] C. Meneveau and J. Katz, Scale-inveriance and turbulence models for large eddy simulation, Annu. Rev. Fluid Mech. 32, 1 (2000).

[14] A. Cimarelli and E. De Angelis, The physics of energy transfer toward improved subgrid-scale models, Physics of Fluids 26 (2014).

[15] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, A dynamical subgrid-scale eddy viscosity model, Phys. Fluids A 3, 1760 (1991).

[16] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, Erratum: "A dynamic subgrid-scale eddy viscosity model" [Phys. Fluids A 3, 1760 (1991)], Physics of Fluids A: Fluid Dynamics 3, 3128 (1991).

[17] D. K. Lilly, A proposed modification of the Germano subgrid-scale closure method, Phys. Fluids A 4, 633 (1992).

[18] S. Liu, C. Meneveau, and J. Katz, On the properties of similarity subgrid-scale models as deduced from measurments in a turbulent jet, J. Fluid Mech. 275, 83 (1994).

[19] R. A. Clark, J. H. Ferziger, and W. C. Reynolds, Evaluation of subgrid-scale turbulence models using a fully simulated turbulent flow, Engng Dept., Stanford Univ. Rep. TF-9 (1977).

[20] B. Vreman, B. Geurts, and H. Kuerten, Large-eddy simulation of the temporal mixing layer using the Clark model, Theoretical and Computational Fluid Dynamics (1995).

[21] A. Cimarelli, A. Abbà, and M. Germano, General formalism for a reduced description and modelling of momentum and energy transfer in turbulence, Journal of Fluid Mechanics 866, 865 (2019).

[22] M. Germano, The two-point average and the related subgrid model, in Proceedings of the Fifth International Symposium on Turbulence and Shear Flow Phenomena (TSFP5) (27-29 August 2007, Munich, Germany, 2007) p. 6.

[23] P. Valente and J. Vassilicos, The energy cascade in grid-generated non-equilibrium decaying turbulence, Physics of Fluids 27, 045103 (2015).

[24] F. Alves Portela, G. Papadakis, and J. Vassilicos, The turbulence cascade in the near wake of a square prism, J. Fluid Mech. 825, 315 (2017).

[25] A. N. Knutsen, P. Baj, J. M. Lawson, E. Bodenschatz, J. R. Dawson, and N. A. Worth, The inter-scale energy budget in a von Kármán mixing flow, Journal of Fluid Mechanics 895 (2020).

[26] P. Beaumard, P. Braganca, C. Cuvier, K. Steiros, and J. Vassilicos, Scale-by-scale non-equilibrium with kolmogorov-like scalings in non-homogeneous stationary turbulence, J. fluid Mech. 884 (2024).

[27] R. J. Hill, Exact second-order structure-function relationships, Journal of Fluid Mechanics 468, 317 (2002).

[28] R. J. Hill, The approach of turbulence to the locally homogeneous asymptote as studied using exact structure-function equations, ArXiv:physics (2002).

[29] L. Danaila, J. Krawczynski, F. Thiesset, and B. Renou, Yaglom-like equation in axisymmetric anisotropic turbulence, Physica D: Nonlinear Phenomena 241, 216 (2012).

[30] J. Chen and J. Vassilicos, Scalings of scale-by-scale turbulence energy in non-homogeneous turbulence, J. Fluid Mech. 938, A7 (2022).

[31] H. Yao, M. Schnaubelt, A. S. Szalay, T. A. Zaki, and C. Meneveau, Comparing local energy cascade rates in isotropic turbulence using structure-function and filtering formulations, Journal of Fluid Mechanics 980, A42 (2024).
[32] J. C. Vassilicos and J.-P. Laval, Scale-by-scale non-equilibrium in turbulent flows (Cambridge University Press, 2025).

[33] R. J. Hill, Equations relating structure functions of all orders, Journal of Fluid Mechanics **434**, 379 (2001).

[34] T. Yasuda and J. Vassilicos, Spatio-temporal intermittency of the turbulent energy cascade, J. Fluid Mech. **853**, 235 (2018).

[35] J. Domaradzki and E. Saiki, Backscatter Models for Large-Eddy Simulations, Theoretical and Computational Fluid Dynamics **9**, 75 (1997).

[36] S. Cerutti and C. Meneveau, Intermittency and relative scaling of subgrid-scale energy dissipation in isotropic turbulence, Physics of Fluids **10** (1998).

[37] T. Aoyama, T. Ishihara, Y. Kaneda, M. Yokokawa, K. Itakura, and A. Uno, Statistics of Energy Transfer in High-Resolution Direct Numerical Simulation of Turbulence in a Periodic Box, Journal of the Physical Society of Japan **74** (2005).

[38] S. Goto, A physical mechanism of the energy cascade in homogeneous isotropic turbulence, Journal of Fluid Mechanics **605**, 355 (2008).

[39] T. Ishihara, T. Gotoh, and Y. Kaneda, Study of high-Reynolds number isotropic turbulence by direct numerical simulation, Ann. Rev. Fluid Mech. **41**, 165 (2005).

[40] M. Germano, A direct relation between the filtered subgrid stress and the second order structure function, Phys. Fluids **19**, 038102 (2007).

[41] M. Germano, The simplest decomposition of a turbulent field, Physica D **241**, 284 (2012).

[42] L. Fang, L. Shao, J.-P. Bertoglio, G. X. Cui, and Z. S. Xu, C. X. and Zhang, An improved velocity increment model based on Kolmogorov equation of filtered velocity, Phys. Fluids **21**, 1 (2009).

[43] C. Hartel, L. Kleiser, F. Unger, and R. Friedrich, Subgrid-scale energy transfer in the near wall region of turbulent flows, Phys. Fluids **6**, 3130 (1994).

[44] A. Leonard, Energy cascade in large eddy simulations of turbulent fluid flow, Advances in Geophysics **08**, 60464 (1974).

[45] A. Vincent and M. Meneguzzi, The spatial structure and statistical properties of homogenous turbulence, J. Fluid Mech. **225**, 1 (1991).

[46] M. E. Brachet, D. I. Meiron, S. A. Orszag, B. G. Nickel, R. H. Morf, and U. Frisch, Small-scale structure of the Taylor-Green vortex, Journal of Fluid Mechanics **130**, 411 (1983).

[47] H. Larssen and J. C. Vassilicos, Spatio-temporal fluctuations of interscale and interspace energy transfer dynamics in homogeneous turbulence, Journal of Fluid Mechanics **969**, A14 (2023).

[48] M. Oberlack, Invariant modeling in large-eddy simulation of turbulence, Center for Turbulence Research, Annual Research Briefs, 3 (1997).

[49] P. Beaumard, J.-P. Laval, and J. Vassilicos, Derivation of a new les model approximated from exact two-point equations and evaluation in a taylor-green flow, Physical Review Fluids (2025).

[50] R. A. Clark, J. H. Ferziger, and W. C. Reynolds, Evaluation of subgrid-scale models using an accurately simulated turbulent flow, Journal of Fluid Mechanics **91** (1979).

[51] B. Dubrulle, Beyond Kolmogorov cascades, Journal of Fluid Mechanics **867** (2019).

[52] S. Goto and J. C. Vassilicos, Local equilibrium hypothesis and Taylor's dissipation law, Fluid Dynamics Research **48**, 021402 (2016).