In some applications of large-eddy simulation (LES), in addition to providing a closure model for the subgrid-scale stress tensor, it is necessary to also provide means to approximate the subgrid-scale velocity field. In this work, we derive a new model for the subgrid-scale velocity that can be used in such LES applications. The model consists in solving a linearized form of the momentum equation for the subgrid-scale velocity using a truncated Fourier-series approach. Solving within a structured grid of statistically homogeneous sub-domains enables the treatment of inhomogeneous problems. It is shown that the generated subgrid-scale velocity emulates key properties of turbulent flows, such as the right kinetic energy spectrum, realistic strain–rotation relations, and intermittency. The model is also shown to predict the correct inhomogeneous and anisotropic velocity statistics in unbounded flows. The computational costs of the model are still of the same order as the costs of the LES.

## I. INTRODUCTION

The concept of large-eddy simulation (LES) has emerged as an enormously successful tool to describe turbulent flows in their many different varieties. The property of turbulence of rapidly decreasing kinetic energy with increasing wave numbers enables capturing the majority of the kinetic energy of a flow by only resolving its largest scales. In many applications, a decent approximation of the large scales of the flow already allows a sufficient understanding of the underlying processes governing the main features of the flow. The present paper is dedicated to problems where this is not the case.

The majority of applications that suffer from missing subgrid-scale velocity content in LES can be found in multiphase flows. The LES of dispersed multiphase flows performs notoriously poorly with regard to predicting multi-particle statistics.^{1–3} Furthermore, when droplets are deformable, the lack of access to the unresolved rotation and strain of the velocity field renders the modeling of their deformation challenging.^{4} Even for simulation methods that at least partially resolve interfaces, a subgrid-scale velocity closure is needed, for example, as in the dual-scale approach of Herrmann and coworkers.^{5,6} Apart from multiphase flows, information about the subgrid-scale velocity may for example also be used to accurately predict the aerodynamical loads on wind turbines that are immersed in the planetary boundary layer,^{7} the chemical reactions and location of the flame front in reactive flows,^{8} and combustion applications or heat transfer with significant turbulent mixing due to the subgrid-scale velocity.

A perfect model for the subgrid-scale velocity in the context of LES would satisfy the governing flow equations [i.e., the Navier–Stokes equations (NSE)], which can only be achieved with the computational cost of a direct numerical simulation (DNS). Characteristic properties indicating that a generated velocity field is similar to that of a real turbulent flow are, among others, the kinetic energy spectrum, the non-Gaussian distribution of derived flow quantities, and the relations between strain and rotation of the velocity. In addition to the physical properties of the generated subgrid-scale velocity, additional demands can be formulated so as to guarantee the practical applicability of the model: (i) the computational costs of the model have to be in the same order of magnitude as the costs of a LES; (ii) if the model contains tuning parameters, the solution should not sensitively depend on their choice; (iii) the model should not be limited to homogeneous and isotropic flows. Homogeneous isotropic turbulence (HIT) constitutes an important case for theoretical studies of turbulent flows, but rarely occurs in realistic flow configurations; and (iv) finally, it is desirable for the model to be implemented and parallelized straightforwardly, as well as to be appendable to any LES flow solver.

A variety of different approaches for the approximate reconstruction of the subgrid-scale velocity have emerged in the past decades. Scotti and Meneveau^{9} derived a fractal interpolation technique to reconstruct the subgrid-scale velocity. This purely geometrical reconstruction does not incorporate important physical properties of turbulence. A related idea, but with a realization in spectral space, is the so-called “kinematic simulation.”^{10–14} With this approach, a divergence-free velocity field that matches a given kinetic energy spectrum is generated. Apart from the assumption of a spatially constant kinetic energy spectrum, it is difficult to define a physically realistic temporal evolution of the modes. Domaradzki *et al.*^{15,16} proposed a reconstruction of the subgrid-scale velocity by a defiltering operation that incorporates information of the non-linear term. In the work of Park *et al.*,^{17} a deconvolution is derived that is dynamically adjusted such that the model is either kinetic energy or dissipation consistent with the subgrid-scale model. However, the model does not introduce length scales that are smaller than the cutoff length scale of the LES grid. A combination of the models of Domaradski *et al.* and Park *et al.* was proposed by Bassenne *et al.*^{18} The idea is a stepwise reconstruction of the subgrid-scale velocity based on alternating applications of spectral extrapolation and dynamical approximate deconvolution. Although key parameters, such as the second invariant of the velocity gradient tensor, can be predicted well, a severe drawback of the method is the projection operation that has to be carried out in real space on a grid comparable to the DNS grid. Thus, the computational costs are of the order of a DNS. The class of the multilevel methods relies on successive explicit filtering operations of the Navier–Stokes equations. Within that class, the variational multiscale method has emerged that uses a Galerkin method to solve the linearized residual velocity equations together with the filtered Navier–Stokes equations (FNSE). Apart from this linearization, no modeling assumptions are introduced.^{19,20} However, even the solution of a linear set of equations with a high resolution leads to a large computational overhead compared to a classical LES on a coarse grid. A second type of multilevel methods uses more than two levels of velocity scales and is based on temporarily frozen velocity of the respective subgrid scales (see, e.g., Terracol *et al.*^{21}). Another approach is to derive a set of auxiliary subfilter equations from the Navier–Stokes equation using localized modes that are extracted via the Gabor transform.^{22,23} The conducted studies have been of a rather theoretical nature or have only been performed in two dimensions. More recently, Ghate and Lele^{24,25} further developed these ideas and achieved excellent agreement of the kinetic energy spectrum and velocity correlations in a boundary layer with the generated velocity and DNS results. However, for long simulation times, the quenching of the modes, which constitutes an application-specific external intervention of relatively high complexity, is required.

As it is clear from above, existing models for the subgrid-scale velocity are still limited in computational efficiency, practical applicability, and accuracy of the predictions.

In this paper, we derive a new model for the enrichment of the LES with a subgrid-scale velocity including physical assumptions and details for practical realization in Sec. II. Subsequently, two simulation cases are introduced for evaluating the performance of the enrichment model compared to the respective DNS in Sec. III. In Sec. IV, the results of the simulation cases are presented and discussed. Finally, Sec. V concludes the present paper.

## II. THEORY

Despite the main focus of this section lying in the modeling of the subgrid-scale velocity, the modeling of the filtered velocity (i.e., the LES) is also briefly summarized for the sake of completeness. The aim of the newly proposed model is to approximate the velocity field of a DNS using only information that is available in a LES.

First, the physical background is provided in Sec. II A. The remaining sections address the assumptions, the discretization, and the realization of the model in detail.

### A. Scale decomposition and the assumptions of the rapid distortion theory

The foundation of LES and of the following discussions is that any physical quantity of interest, $\varphi $ (e.g., velocity and pressure), is decomposed in a filtered (later referred to as LES part) and a subfilter part (later also referred to as subgrid part). The filtered part of a physical quantity, $\varphi \u0303$, may be defined by the following spatial low-pass filtering operation, as proposed by Leonard:^{26}

where *G* denotes a filter kernel. This filtering kernel determines the shape of the spectrum of the resulting filtered quantity. The subfilter-scale contribution, $\varphi \u2032$, is the difference between the original quantity, $\varphi $, and the filtered quantity, $\varphi \u0303$. Note that in a LES, the filtering operation is typically not explicitly performed but rather arises naturally from a numerical bandwidth limitation such as a (too) coarse grid. Furthermore, space and time are coupled in a flow by the velocity which implicitly introduces a filtering in time. However, the definition in Eq. (1) is commonly used for derivations in the LES context and is therefore also adapted in the present paper. Even if not fully consistent, the notations $u\u0303$ and $u\u2032$ are used for referring to the LES velocity and the modeled subgrid-scale velocity, respectively.^{27}

The incompressible Navier–Stokes equations, with a fluid density *ρ _{f}*, and a kinematic viscosity

*ν*are given by

_{f}where *u _{i}* is the velocity,

*p*the pressure, and

*s*an optional source term. Applying the explicit filtering operation (1) to the Navier–Stokes equations (NSE) leads to the filtered Navier–Stokes equations (FNSE) that can be written independently of the filter kernel

_{i}*G,*

The additional term $\u2202\tau ij/\u2202xj$ appears because of the application of the filter to the non-linear advection term. The subgrid-scale stress tensor *τ _{ij}* can be expanded using the general relation $\varphi =\varphi \u0303+\varphi \u2032$ to become

Subtracting the FNSE from the NSE results in a set of equations for the subgrid-scale variables

The subgrid-scale equations contain three advection terms, to which we refer to as: non-linear relaxation (second term on the left-hand side), sweeping (third term on the left-hand side), and straining (fourth term on the left-hand side). The coupled system consisting of Eqs. (4)–(8) can be solved equivalently to the NSE. However, from a numerical point of view solving this coupled system is much more challenging than solving the NSE.

Lavel *et al.*^{28} performed a numerical analysis of the coupled system consisting of Eqs. (4)–(8) and omitted the non-linear relaxation. Taking into account only the linear terms in Eqs. (7) and (8) resembles the governing equations for the rapid distortion theory (RDT) that makes use of this linearization to obtain insights of the small-scale turbulence modification under large-scale distortion.^{29} Note that under very high mean strain the viscous term is often also neglected in analytical studies.^{30,31} It emanates from these studies that without the non-linear relaxation term, the turbulent kinetic energy spectrum is overestimated and the intermittency is increased (i.e., the probability distribution function of the velocity increments shows wider tails). With the addition of a turbulent viscosity term, these two effects can essentially be compensated.^{28}

Apart from the energetic and intermittency effects that have been observed as a consequence of the neglected non-linear relaxation term by Laval *et al.*,^{28} the strain–rotation relations of the velocity may also be influenced. This can be quantified with the probability distribution function of the second invariant of the velocity gradient tensor $Q=1/2(\Omega ij\Omega ij\u2212SijSij)$, or the angle between the eigenvectors of the rotation-rate tensor Ω_{ij} and the strain-rate tensor *S _{ij}* as for example investigated by Horiuti.

^{32}In Sec. IV A, evidence is provided that the non-linear relaxation term contributes to the strain–rotation relations. Furthermore, the term $uj\u2032\u2202ui\u2032/\u2202xj$ redistributes kinetic energy in spectral space, whereas a diffusive term $\u22022ui\u2032/(\u2202xj\u2202xj)$ exclusively removes kinetic energy if the (turbulent) viscosity is positive (see, e.g., Pope

^{30}). It can be concluded that the replacement of the non-linear relaxation term with a turbulent viscosity may approximate the effects of attenuation of intermittency and high wave numbers of the kinetic energy spectrum well, but does not reproduce all physical mechanisms of the original term.

### B. Modeling the filtered and subgrid-scale equations

In applications where only a coarse resolution is feasible, the FNSE can be solved instead of the NSE, using a suitable approximation for *τ _{ij}*. In such a case, all present flow structures can be resolved by the grid, because the high wave number content of the physical quantities is removed. Since in a LES the equations for the small scales are not solved, the subgrid-scale stresses are unknown and require modeling. Many approaches exist that attempt to close the equations for the large scales that typically rely on the use of a turbulent viscosity to mimic at least the energetic effect of additional dissipation by small scales.

^{33–36}This eddy-viscosity type of subgrid-scale models yields a modeled subgrid-scale stress tensor of the form

where *ν _{t}* is the model-specific turbulent viscosity,

*δ*is the Kronecker tensor, and $S\u0303ij$ is the strain-rate tensor of the large-scale velocity, given by

_{ij}Even with a subgrid-scale model that perfectly reproduces the subgrid-scale stress tensor, only one part of the coupled problem is solved, that is, only Eqs. (4) and (5) of the system of equations (4)–(8) are solved. Due to the rapid decrease in the kinetic energy spectrum of a turbulent flow for increasing wave numbers, the filtered velocity already possesses the majority of the kinetic energy, which is sufficient for many applications. However, there are applications for which knowledge of the subgrid-scale velocity field is required. In this paper, we propose a model for the subgrid-scale velocity that enables an approximate solution of the full problem consisting of Eqs. (4)–(8) (including the filtered velocity and the subgrid-scale velocity), but with a computational cost of the same order of magnitude as that of LES.

The proposed model exploits the previous findings of Lavel *et al.*,^{28} who suggest replacing the non-linear relaxation term in Eq. (8) by a turbulent viscosity $\nu t\u2032$. Based on this work, the following expression for the turbulent viscosity is used, which is derived from renormalization groups to be^{37}

where *k* is the wave number, *E*(*k*) is the kinetic energy spectrum, and $C\nu $ is a constant with an analytical value of $C\nu =2/5$. Equation (11) is advantageous for the resulting kinetic energy spectrum compared to a constant turbulent viscosity for the following reason: The omitted non-linear relaxation is not only responsible for attenuating intermittency, but also for redistributing kinetic energy (mainly toward smaller scales). The turbulent viscosity of Eq. (11) is larger for small wave numbers and converges toward zero for high wave numbers. This means that the turbulent viscosity removes more kinetic energy from larger scales than from the smaller scales comparable to the physical mechanism of energy transfer toward small scales.

In addition to the introduction of a turbulent viscosity, a second modeling step is required for coupling the filtered scales with the subgrid scales, that is, modeling the subgrid-scale stress tensor. Assuming the forward energy cascade to be dominant, the influence of the filtered scales on the subgrid scales is mainly characterized by the subgrid scales being supplied with kinetic energy from the filtered scales. A simple way of mimicking this energy supply is by adding an additional source term to the right-hand side of the subgrid-scale momentum equations (8). This is inspired by numerical forcing of turbulence for maintaining statistical steadiness in HIT, in particular by the forcing scheme of Mallouppas *et al.*^{38} It is important to understand that the forcing maintains a desired kinetic energy of the subgrid scales, which may be estimated using the model of Yoshizawa^{39} or the dynamical model of Moin *et al.*^{40} The amount of kinetic energy that is added to the subgrid scales is not necessarily equal to the kinetic energy that is removed from the filtered scales. The enrichment constitutes a more general process if it operates independently of the subgrid-scale model and its choices of constants. Furthermore, modeling errors of the subgrid-scale model weigh much more for the subgrid scales since their kinetic energy is only a small fraction of the total kinetic energy.

### C. Discretization of the subgrid-scale equations

Even with the modeling described in Sec. II B, a numerical solution of Eqs. (7) and (8) (e.g., with a finite volume solver) would still require a grid that is fine enough to resolve the smallest flow structures. Since this would exceed the acceptable amount of computational time of a LES, a different discretization approach is proposed.

We assume that the subgrid-scale velocity $u\u2032$ can be represented by the finite sum

where $Am(t)$ and $Bm(t)$ are vectorial coefficients that depend on time *t*. Equation (12) is equivalent to a truncated Fourier series expansion with *N _{m}* modes, so the coefficients can be interpreted as contributions to the velocity at the respective wave number $km=|km|$, where $km$ is the wave number vector. If the coefficients and the wave number vectors are known, the subgrid-scale velocity can be computed at arbitrary positions using the expansion given by Eq. (12).

The same approach is utilized in kinematic simulations.^{10–14} It consists of determining the coefficients, such that (i) the kinetic energy equals that of a given spectrum, and (ii) the random Gaussian directions of the coefficients lead to a divergence free velocity field, which yields the requirement $km\xb7Am(t)=km\xb7Bm(t)=0$. In kinematic simulations, a continuous velocity field that can be evaluated at every arbitrary position is obtained. The resulting velocity field can lead to decent predictions of second-order Lagrangian statistics of fluid particles.^{14} Since the coefficients do not change in space and the directions are chosen from a Gaussian distribution with fixed mean and standard deviation, the resulting velocity is statistically homogeneous and stationary. Furthermore, it is not clear how to maintain realistic time correlations when the kinetic energy spectrum varies in time. These issues are addressed and overcome with the present approach.

Consider the NSE for the subgrid scales with the turbulent viscosity, from Eq. (11), that replaces the non-linear relaxation term

where *f _{i}* is a forcing term that models the kinetic energy transfer from the filtered scales to the subgrid scales. Note that we assume the spectral turbulent viscosity constant in space, so it is treated similar to the molecular viscosity in the viscous term. In the next step, the series expansion of the subgrid velocity in Eq. (12) is inserted into the modeled momentum equations (14) for every term of the sum separately and the pressure term is dropped. Sorting the result by sine and cosine terms the following equations for the coefficients $Am(t)$ and $Bm(t)$ are obtained:

where $\Delta t$ indicates the numerical time step. The time level of the coefficients is indicated with the index *n*. For this and all the following equations, no implicit summation over the index *m* is carried out. Without the pressure term, only a preliminary solution (indicated with an asterisk) is obtained. These coefficients do not lead to a divergence-free velocity field. However, if the pressure is formally expanded in spectral space and the divergence of the difference of Eqs. (15) and (16) and the same equations with pressure is computed, the following projection operation (details in the Appendix) is obtained to generate a set of coefficients that lead to a divergence-free velocity field:

The computational cost of the solution of one time step consisting of Eqs (15)–(18) scales linear with the number of modes *N _{m}*. In fact, the coefficients of the new time step are obtained by the explicit solution of two sets of algebraic equations, which makes their solution affordable.

Finally, determining the forcing terms $fm,i$ and $gm,i$ is addressed. Their modeling consists of two components: the magnitude and the direction. Ideally, the direction of the forcing terms is equal to the direction of the divergence of the subgrid-scale stress tensor, as defined in Eq. (9). However, this non-linear expression consists of contributions from the subgrid-scale velocity, the filtered velocity, their interaction, and the explicit filtering operations. In the present model, the directions are chosen from a uniform distribution over a spherical shell

with

In order for the forcing directions to have realistic time correlation, the following algorithm is used to propagate the directions in time:

where $Vm,0$ is a random vector that is newly generated for every wave number and every time step, and $Vmn\u22121$ is the normalized direction of the previous time step. The vector $Vmn$ is the direction of the forcing term at the current time step. Note that both forcing terms have directions that are propagated independently. Similar to the divergence of the subgrid-scale stress tensor, the vector field $Vmn$ is not divergence free ($km\xb7Vmn\u22600$). The weights *α _{m}* and

*β*determine how fast the directions change in time, and are computed as follows:

_{m}The rate of change of the directions depends on the wave number because the weights are determined based on the eddy-turnover time of the respective mode $TE,m=|km|3E(|km|)$.^{10} The kinetic energy spectrum $E(|km|)$ can be directly computed from the coefficients as

where $\Delta km$ is the distance between consecutive wave numbers for the wave number *k _{m}*. Similar to the forcing scheme of Mallouppas

*et al.*,

^{38}the source term in the subgrid-scale equation is

where the trigger velocity of the respective mode follows the slope of the inertial range and is given by

The forcing term is zero if the kinetic energy of the subgrid-scale velocity field *K _{actual}* is equal the desired kinetic energy

*K*. The trigger velocity is scaled such that it has the same characteristic magnitude as the resulting subgrid-scale velocity. In order to do so, the energy of the source has to be determined as

_{desired}where $VA,m$ and $VB,m$ are the directions for the forcing of the respective equation of the coefficients $Am$ and $Bm$, which are determined from Eq. (21). The kinetic energy that is actually present in the subgrid scale is computed analogously

The forcing term $gi,m$ has the same magnitude as $fi,m$, but is based on a separate evolution of the direction $VB,m$. Since the expression for $gi,m$ is almost identical to $fi,m$, but with $VB,m$ instead of $VA,m$ in Eq. (25), the explicit formula is not shown here for conciseness. The desired kinetic energy of the subgrid-scale velocity is estimated with the approach of Yoshizawa^{39}

with the constant $CI=0.0826$ and the filter width Δ. The constant can also be computed dynamically with the model of Moin *et al.*^{40}

### D. Transition to a sub-domain-based discretization

Since the model described in Sec. II C is constructed in Fourier space [i.e., the coefficients $Am(t)$ and $Bm(t)$ are spatially constant], it can only generate subgrid-scale velocities with spatially constant statistics. However, many types of flows are inhomogeneous, which is the motivation for modifying the above approach, so that inhomogeneous statistics across a domain can be achieved. Instead of assuming that the coefficients obtained from Eqs. (15)–(18) are constant globally, they are defined on a structured grid of sub-domains, coarser than the LES grid. This induces motion of the flow field whose scale is of the order of the sub-domain size. However, considering that the energy of the corresponding induced velocity is small compared to the energy of the large scales, its effect is minor compared to other modeling assumptions in LES. Every sub-domain possesses a distinct set of coefficients. In each sub-domain, the statistics are assumed homogeneous, but each sub-domain can have different statistics. Thus, the coefficients formally do not only depend on time, but also vary between the sub-domains. Therefore, the coefficients have an additional index for the respective sub-domain and are denoted as $Ad,m$ and $Bd,m$. This implies that the spatial variations of scales of the order of the domain size are represented by spatially varying coefficients and the high wave number fluctuations by the trigonometric functions in Eq. (12). The quantities from the LES that are required to solve Eqs. (15)–(18) are averaged within the respective sub-domain, as indicated by $\u27e8\xb7\u27e9domain$. In particular, this applies to the LES velocity $\u27e8u\u0303i\u27e9domain$, the gradient of the LES velocity $\u27e8\u2202u\u0303i/\u2202xj\u27e9domain$, and theoretically also the estimation of the kinetic energy of the subgrid-scale velocity $\u27e8Kdesired\u27e9domain$. The latter has severe consequences for the resulting subgrid-scale field, which are discussed in Sec. II F. A domain average of the subgrid-scale kinetic energy is thus not applied in the present model. Instead, the subgrid-scale kinetic energy is estimated at the center of the sub-domain.

Another consequence of this sub-domain discretization is that the spatial derivatives also possess contributions of the coefficients. This leads to a modified predictor step for the preliminary coefficients

The temporal and spatial dependencies of the coefficients are omitted for the sake of simplicity.

Figure 1 summarizes the proposed velocity enrichment strategy based on a LES. The first step is to perform a time step of the LES and average the required quantities in the sub-domains. Furthermore, the subgrid-scale kinetic energy is estimated and the direction and magnitude of the forcing terms are computed. The choice of the LES solving framework is arbitrary. Even a LES on an unstructured grid is possible, in which case solely the averaging has to be adapted. In the next step, Eqs. (29), (30), (17), and (18) are solved in every sub-domain to obtain a distinct set of coefficients for every sub-domain. This step is rather inexpensive, since the equations are solved fully explicit and independent for every sub-domain. The parallelization is very simple and efficient if it is carried out in spectral space. If every processor solves a portion of the wave numbers of all sub-domains, no communication is necessary. The final step consists of adding up the LES velocity and the subgrid-scale velocity, which is obtained by computing Eq. (12) at the desired locations.

### E. Treatment of the discontinuities

In some applications of the enrichment model, a spatially continuous reconstructed subgrid-scale velocity field may be required. Due to the proposed sub-domain discretization, this cannot be guaranteed at the boundaries of the sub-domains. A naive interpolation of the subgrid-scale velocity or the coefficients of the series expansion in Eq. (12) would not conserve the property of a divergence-free velocity field. However, since the divergence is only non-zero at sub-domain boundaries, the portion of the influenced regions is small compared to the entire simulation domain.

In this section, we present a modified interpolation between sub-domains that leads to a continuous and divergence-free velocity field. The modeled momentum equations are, in general, not satisfied within the interpolation regions. A compact interpolation kernel should be preferred in order to keep the regions with the subgrid-scale velocity that satisfies the model equations as large as possible.

The general idea is to add a correction velocity in a direction tangential to the boundary for every mode. We assume that the interpolated subgrid-scale velocity $uint\u2032$ consists of a contribution of the lower sub-domain $ulow\u2032$, the upper sub-domain $uup\u2032$, and a correction velocity $ucorr\u2032$,

where *W* is the interpolation kernel function. In the following, the interpolation in x direction is derived. The other directions can be obtained analogously. For this case, the interpolation kernel function *W*(*x*) is exclusively a function of *x* and the correction velocity has only contributions in y and z directions. The required criterion for the interpolated velocity is

Furthermore, we assume that the velocity fields can be decomposed in the following way:

where $u\u2032,m$ is the velocity contribution corresponding to the wave number $|km|$. Since the divergence is a linear operation, Eq. (32) is satisfied if the divergence of every velocity contribution $u\u2032,m$ is zero

In the x direction, the divergence-free requirement is expressed as

At this point, a choice has to be made for how the divergence of the velocity field is distributed over the boundary tangential velocities. One satisfactory way is to either set $wcorr\u2032,m=0$, which leads to $vcorr\u2032,m=\Psi m(x)/km,2$ or set $vcorr\u2032,m=0$, which leads to $wcorr\u2032,m=\Psi m(x)/km,3$. In order to keep the correction velocity magnitude small, $min(|vcorr\u2032,m|,|wcorr\u2032,m|)$ decides which case is used. The correction velocity kernel is defined as

The correction theoretically works with every interpolation kernel *W*(*x*). However, we suggest to use a sigmoid function based on the signed distance to the sub-domain boundary *r,*

where *α* is a parameter to control the sharpness of the interpolation kernel.

An extension to three dimensions is straightforward as long as the interpolation kernels of different directions are not overlapping. This leads to complicated integrations that have to be carried out analytically, which is not possible in some cases. If the overlapping of the interpolation kernels is ignored, the edges and corners of the sub-domains are not divergence free. With a sharp interpolation kernel, these regions are tiny and most likely unimportant for the majority of the applications. If the integrand is approximated with a Taylor series of the interpolation kernel, a nearly divergence-free field can also be obtained in the regions with overlapping interpolation kernels.

The divergence of the modeled subgrid-scale velocity with and without interpolation between sub-domains is depicted in Fig. 2. The detailed setup and the definition of the reference time *T _{ref}*, which is used for normalization, are specified in Sec. III B 1. The divergence is computed with central differences for the spatial derivatives of the subgrid-scale velocity field that is sampled on a DNS-grid (i.e., a grid that can resolve the Kolmogorov length scale). In Fig. 2(a), where no interpolation between sub-domains is used, it can be observed that the divergence is significant between the sub-domains. With the proposed interpolation in Fig. 2(b), the divergence is approximately zero almost everywhere; solely at the edges between four sub-domains does non-zero divergence occur.

The divergence-free interpolation constitutes an option for applications where the discontinuity in the subgrid-scale velocity is not acceptable. The results of the present paper are obtained without the interpolation (i.e., with a discontinuous subgrid-scale velocity). However, it is shown in Sec. IV A that the influence of the interpolation on the presented statistics is negligible.

### F. The role of the estimated subgrid-scale kinetic energy

As mentioned in Sec. II D, the subgrid-scale kinetic energy is evaluated in one grid cell within the sub-domain instead of averaged over all grid cells in the sub-domain. This section justifies this choice.

The estimation of the kinetic energy of the subgrid-scale velocity strongly influences the predictions of the presented model, as it directly impacts the kinetic energy spectrum of the generated velocity field. For HIT, the estimation of the kinetic energy with Eq. (28) leads to average subgrid-scale kinetic energies that are close to the values obtained from a filtered DNS. More specifically, the estimated subgrid-scale kinetic energy lies within a ten percent margin of the subgrid-scale kinetic energy obtained from the DNS of the HIT case described in Sec. III B using a spectrally sharp filter.

Apart from the mean value of the kinetic energy of the subgrid-scale kinetic energy, its probability distribution is of importance. A characteristic aspect of turbulent flows is that the probability distribution function of the velocity gradients (and other quantities) exhibits a non-Gaussian behavior, which is widely known as “intermittency.”^{41} Typically, as the Reynolds number increases, the tails of the probability distribution functions (PDFs) get wider, and that is, high-intensity events become more likely. Since the effect of the subgrid-scale velocity on the filtered scales in a LES is commonly modeled using a turbulent viscosity, a spatially averaged LES may be interpreted as a DNS with smaller Reynolds number. Thus, the PDFs can be expected to tend to a Gaussian distribution. This is confirmed in Sec. IV A.

Due to the forcing of the subgrid-scale velocity, its distribution directly depends on the distribution of the subgrid-scale kinetic energy. Distributions of derived values, such as the subgrid-scale velocity gradients, thus also depend on the distribution of the subgrid-scale kinetic energy. This is why the PDF of the subgrid-scale kinetic energy is of great importance.

In Fig. 3, the probability distribution function of the estimated subgrid-scale kinetic energy is shown together with a squared Gaussian distribution. The simulation configuration again refers to the HIT case described in Sec. III B. The kinetic energy is centered with respect to the mean $\u27e8Kest\u27e9$ and normalized with standard deviation *σ _{k}*. As expected, the two distributions do not coincide. More practically, relevant is how the PDF of the estimated subgrid-kinetic energy averaged in sub-domains consisting of 2

^{3}LES grid cells compares to this. It is not surprising that the averaging leads to a distribution that is closer to a Gaussian distribution. As a consequence, derived quantities of the modeled subgrid-scale velocity field possess different PDFs, dependent on the averaging volume of the estimated subgrid-scale kinetic energy.

Since the averaging of the subgrid-scale kinetic energy within the sub-domains leads to artificial modification of its PDF, the averaging is avoided in the application of the model. Instead, the kinetic energy of the subgrid-scale velocity is chosen to be equal to the one computed for one arbitrary LES grid cell in the sub-domain (e.g., the cell in the sub-domain center).

Sections III–IV will show the features of the model and validate it for two flow configurations.

## III. SIMULATIONS

The performance of the model proposed in this paper is evaluated on two test cases, HIT and a turbulent unbounded shear flow. The reference DNS is carried out on a very fine grid that captures all present length and time scales down to the Kolmogorov scales. LES that are enriched with the presented model are compared with the DNS results. Note that wall boundary conditions significantly increase the complexity of the flow and already constitute a problem for LES without modeling of the subgrid-scale velocity. This is why the scope of this paper is limited to unbounded flow configurations.

### A. Numerical solution of the Navier–Stokes equations and LES equations

For the two DNS (HIT and turbulent shear flow), the NSE are solved in their form as given in Eqs. (2) and (3). Since all time and length scales down to the Kolmogorov scales are resolved, no further modeling is applied to the set of governing physical equations. The LES is performed on a grid that is too coarse to resolve the full turbulent spectrum and is based on the FNSE given by Eqs. (4) and (5). Due to the presence of the subgrid-scale stress tensor, these equations are not closed. That is why the following set of equations is solved in the case of a LES that are referred to as LES equations:

The modified pressure $P$ accounts for the trace of the subgrid-scale stress tensor and the turbulent viscosity, *ν _{t}*, models the effect of additional dissipation by the omitted subgrid scales. For the HIT test case, the turbulent viscosity is determined by the Smagorinsky model

^{33}

where *C _{S}* is the Smagorinsky constant, Δ is the filter width, and $|S\u0303ij|$ is the following norm of the filtered strain rate tensor:

In shear-dominated turbulent flows, it can be advantageous to replace the Smagorinsky model with the Vreman model,^{36} which is known to produce more physical results in transitional shear flows. With this model, the turbulent viscosity is given by

with the filtered velocity gradient tensor

and

In the case of HIT, a statistically steady state is desired for the evaluation of the model. To achieve this, the same amount of kinetic energy that is dissipated has to be inserted into the flow. In addition to the constant kinetic energy, further requirements are that the flow is statistically homogeneous and isotropic. In the present work, the forcing scheme of Mallouppas *et al.*^{38} is used. The applied procedure is very similar to the forcing of the small scales that is described in Sec. II C

However, there are two differences. Instead of estimating the desired kinetic energy, *K _{desired}*, its value is given as a parameter and determines the kinetic energy in the flow. Furthermore, the trigger velocity $vitrigger$ is determined from the following series expansion:

The wave number vector $km,j$ has a magnitude corresponding to the triggered wave number and has random directions. The phase *ψ _{m}* is chosen randomly from a uniform distribution in the range $0\u2264\psi \u22642\pi $. The velocity coefficients $Um,i$ are chosen to lead to a divergence-free trigger velocity and thus fulfills the condition $Um\xb7km=0$. Its magnitude corresponds to a von Kármán kinetic energy spectrum in the triggered wave number range and is given by

where *E _{vk}* is a von Kármán spectrum and $\Delta k$ the difference between two consecutive discrete wave numbers. Similar to the directions in the forcing of the small scales, $vitrigger$ is smoothly propagated in time by applying Eqs. (21) and (22). The trigger velocity consists only of modes within a predefined wave number range. This has the reason that isotropy of the generated velocity field cannot be achieved if the largest possible scales are triggered.

^{42}

In this work, the NSE and the LES equations are discretized with a finite-volume approach that converges with second order in space and time. The numerical solution is similar to the one described in Denner *et al.*,^{43} but is carried out on an equidistant structured Cartesian grid with the corresponding geometrical simplifications that result from constant cell volumes, cell face areas, and normal vectors that are aligned with the coordinate axes. The continuity equation and the three-momentum equations are solved in a single equation system for multiple iterations every time step. The continuity equation is coupled with the momentum equations using momentum-weighted interpolation, a concept that was originally introduced by Rhie and Chow.^{44} As a result, two distinct velocity fields exist in the numerical framework, a cell centered velocity that satisfies the momentum balance with high accuracy and a face centered velocity that conserves the mass. With increasing spatial resolution, the velocities converge toward each other. The results that are presented in this paper are based on the cell-centered velocity.

### B. Flow configurations

#### 1. Homogeneous isotropic turbulence

The first of the two considered test cases is forced HIT. The simulation domain is a cube with the edge length *L* and periodic boundary conditions in all three directions. A statistically steady turbulence is obtained by applying the forcing procedure described in Sec. III A. In order to obtain a homogeneous, isotropic velocity field, the forcing is exclusively applied between the wave number *k _{start}* and

*k*. Physical parameters of the resulting turbulent velocity field are the Taylor–Reynolds number $Re\lambda =\lambda urms/\nu f$, the turbulent Reynolds number $Rel=l11urms/\nu f$, the Kolmogorov length scale

_{end}*η*, the Taylor microscale

*λ,*and the longitudinal (

*l*

_{11}) and transverse (

*l*

_{22},

*l*

_{33}) integral length scales, and are summarized in Table I. Spatial quantities are presented normalized by the domain length and temporal quantities are normalized by a reference time $Tref=L/urms$, with the root mean square velocity $urms=2/3K$ and the total kinetic energy

*K*.

Parameter . | Value . |
---|---|

$Re\lambda $ | 75 |

Re _{l} | 205 |

$\eta /L$ | 0.0017 |

$\tau \eta /Tref$ | 0.0075 |

$\lambda /L$ | 0.029 |

$l11/L$ | 0.079 |

$l22/L$ | 0.49 |

$l33/L$ | 0.38 |

$kstartL/2\pi $ | 3 |

$kendL/2\pi $ | 6 |

Parameter . | Value . |
---|---|

$Re\lambda $ | 75 |

Re _{l} | 205 |

$\eta /L$ | 0.0017 |

$\tau \eta /Tref$ | 0.0075 |

$\lambda /L$ | 0.029 |

$l11/L$ | 0.079 |

$l22/L$ | 0.49 |

$l33/L$ | 0.38 |

$kstartL/2\pi $ | 3 |

$kendL/2\pi $ | 6 |

The solution of the DNS is obtained on a numerical grid with 256^{3} cells. This leads to a product of the largest resolved wave number *k _{max}* and the Kolmogorov length scale of $kmax\eta =1.37$.

The LES of the same case is solved on a grid consisting of 32^{3} cells. The influence of the omitted subgrid scales is accounted for by adding a turbulent viscosity that is obtained from the Smagorinsky model with a model constant $CS=0.1$ and a filter width Δ corresponding to the size of a computational cell.

The resulting kinetic energy spectra *E*(*k*) of the DNS and the LES are shown in Fig. 4. The DNS spectrum possesses a pronounced inertial range that is characterized by a slope of $|km|\u22125/3$. It can be observed that the kinetic energy spectrum of the LES resembles the one of the DNS well for the smaller wave numbers. The cutoff wave number of the LES lies within the inertial range.

For a homogeneous isotropic flow, the spatial autocorrelations $B\alpha \alpha (1)$ are defined as^{45}

where *r* denotes the distance between the evaluated velocities and $e1$ the unit vector in the first coordinate direction. Thus, $B11(1)$ indicates the longitudinal and $B22(1)$ and $B33(1)$ the transverse spatial autocorrelation functions, respectively. As already mentioned, isotropy is not achieved if the smallest wave numbers are included in the forcing. Based on the spatial autocorrelations, the integral length scales $l\alpha \alpha $ are defined as

#### 2. Turbulent shear flow

The second case that is investigated is a turbulent shear flow that is driven by a source term in the momentum equation. Figure 5 sketches the simulation setup. The domain is a cube with edge length L and periodic boundary conditions in all the three directions. The momentum source acts in the x direction and varies with the shape of one period of a sine-profile in the y direction

The amplitude is given by *s _{max}* and $e1$ indicates the unit vector in x direction. Assuming a shear velocity $ushear=smaxL$, the shear Reynolds number can be defined as $Reshear=ushearL/\nu f$. The value of this and other physical parameters is summarized in Table II. Similar to the shear velocity, a characteristic shear time can be defined for normalization $Tshear=L/ushear$.

Parameter . | Value . |
---|---|

Re _{shear} | 3115 |

$\eta /L$ | 0.0029 |

$\tau \eta /Tshear$ | 0.026 |

Parameter . | Value . |
---|---|

Re _{shear} | 3115 |

$\eta /L$ | 0.0029 |

$\tau \eta /Tshear$ | 0.026 |

A DNS with 128^{3} grid cells fully resolves the length scales down to the Kolmogorov scale of the shear flow ($kmax\eta =1.16$). In addition to the DNS, a LES with 16^{3} grid cells is performed. Since it is well known that the Smagorinsky model is inaccurate in the prediction of shear-dominated flows (see, e.g., Lévêque *et al.*^{46}), the eddy viscosity is computed with the Vreman model using a model constant with the value $CV=0.025$.

Figure 6(a) shows the resulting mean velocity profiles of the DNS and the LES averaged in time. For both simulations, the mean velocity profiles in the *y* and *z* directions are very close to zero. According to the profile of the volumetric force that drives the flow, the mean velocity profiles in *x* direction have the shape of one period of a sine-function. However, the amplitudes of the mean velocity profiles differ significantly.

Figure 6(b) shows the kinetic energy of the DNS and the LES of the turbulent shear flow over time. The turbulent shear flow behaves aperiodic and contains strong fluctuations of flow quantities, for example, the kinetic energy. Periods of pronounced turbulence slow down the main velocity stream, which reduces the production of new turbulent structures. The flow becomes less turbulent and the mainstream accelerates because the resistance due to the turbulence decreases. This causes a chaotic change between periods of high and low kinetic energy. As already pointed out, the mean velocity profile of the LES possesses a higher amplitude than the one of the DNS. The consequence is that the average kinetic energy of the LES is also higher.

## IV. RESULTS AND DISCUSSION

In Sec. IV A, predictions of the presented model are evaluated in HIT and the turbulent shear flow. We investigate the extent with which the actual intermittency and spatial structures can be reproduced by the proposed enrichment strategy as described in Sec. IV A. Furthermore, the ability of the model to generate the anisotropic and inhomogeneous subgrid-scale velocity of the turbulent shear flow is assessed in Sec. IV B.

### A. Enriched velocity of the HIT case

In order to study the influence of different parameters on the resulting predictions of the subgrid-scale velocity, a standard model setup for the enrichment of the HIT case is defined. The standard model setup contains a grid with eight sub-domains per direction, which is of a factor four coarser per direction than the LES grid. The number of modes in the series expansion (12) is *N _{m}* = 100 and the constant in the turbulent viscosity in the enrichment model $\nu t\u2032$ is $C\nu =0.4$, which corresponds to the analytical value obtained from renormalization group theory.

^{37}The kinetic energy of the subgrid-scale velocity is estimated using the model of Yoshizawa

^{39}$Kest=CI\Delta 22S\u0303ijS\u0303ij$, with the originally proposed value of the model constant $CI=0.0826$.

In Fig. 7, slices of the instantaneous velocity magnitudes of the DNS, the LES, and the LES enriched with the proposed model are depicted. The LES does not contain small velocity structures, and the maximum velocity magnitude is still somewhat smaller than the one of the DNS. With enrichment, the velocity magnitude of the DNS is better recovered and smaller velocity structures than in the LES without enrichment are observed. However, the DNS field consists of thin elongated structures, whereas the enrichment yields a less structured velocity field. This missing structure in the subgrid-scale velocity can be traced back to the forcing of the subgrid scales that is applied in randomly chosen directions instead of the actual directions of the subgrid-scale stress tensor. Furthermore, the non-linear relaxation term, which is not part of the model equations, may contribute to the formation of the observed structures in the DNS velocity field.

Figure 8 shows the kinetic energy spectrum that results from the subgrid-scale velocity. The kinetic energy spectrum of the modeled subgrid-scale velocity coincides well with the DNS. The LES and the modeled subgrid-scale spectrum together recover the entire range of scales of the DNS. Note that the subgrid-scale spectrum is obtained directly from the coefficients of the series expansion Eq. (12). If the discrete Fourier transform of the sampled subgrid-scale velocity field, which is discontinuous across the sub-domains, were to be computed, additional spurious spectral content would occur. This is known as the Gibbs phenomenon and is avoided here by computing the kinetic energy spectrum directly from the coefficients in Eq. (23).

The kinetic energy spectrum provides insight into the average kinetic energy belonging to a certain wave number. However, it does not contain information on how likely events with a specific intensity are. A characteristic of turbulent flows is that the PDFs exhibit a higher incidence of high-intensity events, leading to a non-Gaussian (more widely shaped) distribution of the velocity gradients.^{41} Figure 9 shows the PDFs of the longitudinal and transverse velocity gradients *A*_{11} and *A*_{12} normalized with their standard deviations $\sigma A11$ and $\sigma A11$. In both cases, the LES predicts the gradients with smaller magnitudes, and the PDF tends to be closer to a Gaussian distribution than the PDF of the DNS. With the LES including the newly proposed enrichment model with the standard model setup, a much wider distribution of the velocity gradients is achieved that increases the intermittency of the flow compared to the LES. Solely, the very rare events with the greatest magnitude (i.e., the ends of the tails of the PDF) cannot be fully reproduced. It may seem surprising, that despite the Gaussian forcing of the subgrid-scale velocity, a non-Gaussian PDF of the velocity gradients is accomplished. There are two main reasons for the observed non-Gaussian distribution of the velocity gradients of the enriched LES: first, the forcing is only Gaussian within a sub-domain. Since the estimated subgrid-scale kinetic energy shows a non-Gaussian distribution, the forcing is also non-Gaussian on a scale of the domain size, and second, there are other terms in the model equations (13) and (14) apart from the forcing term that verifiably cause intermittency, that is, the sweeping term and the straining term.^{28} It is also worth noting that the PDF of the longitudinal velocity gradient possesses a non-zero skewness. The modeled subgrid-scale velocity obtained from the enrichment model has difficulties in reproducing this asymmetric behavior.

In addition to the distribution of single entries in the velocity gradient tensor, a characteristic feature of turbulence is the distribution of the second invariant *Q* of the velocity gradient tensor, which is defined as

where Ω_{ij} is the rotation-rate tensor and *S _{ij}* is the strain-rate tensor. The second invariant of the velocity gradient tensor can be interpreted as a measure of how strong rotation is compared to straining. The right prediction of

*Q*is of practical importance, for instance in particle-laden flows, since the vorticity-strain relations of a turbulent flow can be directly related to the clustering behavior of inertial particles of very small Stokes numbers that are immersed in the flow.

^{47,48}The ability to predict this behavior is often viewed as important property for structural subgrid-scale models to predict the right amount of particle clustering.

^{13,18,49,50}Note however that the right prediction of

*Q*does not guarantee that the right particle clustering is obtained.

The PDFs of the second invariant of the velocity gradient tensor of the DNS, the LES, and multiple predictions of the enriched LES velocity are compared in Fig. 10. Four different parameters of the model are systematically varied with respect to the standard model setup in order to determine their impact: (i) the number of sub-domains, (ii) the kinetic energy of the subgrid-scale velocity, (iii) the number of modes, and (iv) the model constant $C\nu $ in the turbulent viscosity $\nu t\u2032$ of the subgrid-scale momentum equation.

The most obvious observation in Fig. 10 is that the LES without an enrichment model severely underpredicts the width of the PDF compared to the DNS. The previously defined standard model setup, which is depicted for comparison in all of the figures, improves the width and also the shape of the PDF of the second invariant of the velocity gradient tensor. However, the skewness is somewhat less pronounced than for the DNS.

An indication for the observed bias in the PDF of the second invariant of the velocity gradient tensor may be obtained by decomposing *Q* as follows:

In fact, the divergence of the advective term is proportional to the second invariant of the velocity gradient tensor. Neither the non-linearity in the subgrid-scale momentum equation (8) nor the non-linearities in the subgrid-scale stress tensor are represented correctly in the model. The second invariant of the velocity gradient tensor can be further expanded

In contrast to the local contributions $u\u0303iu\u0303j$ and $ui\u2032uj\u2032$, the non-local terms in Eq. (54) are represented in the model equations and contribute to the improved results of the PDF compared to the LES. This suggests that for further improvements of the predictions of the second invariant of the velocity gradient tensor, the non-linear relaxation term has to be included and the forcing has to be modified.

Next to the standard setup of the enrichment model, other configurations with variation of the four previously defined influence parameters are also shown. Figure 10(a) compares the obtained PDFs for the second invariant of the velocity gradient tensor with a varying number of sub-domains that are used for the solution of the model equations. Instead of the 8^{3} sub-domains of the standard model setup, two further simulations are performed with 16^{3} and 32^{3} sub-domains, respectively. The PDFs of the three cases almost coincide, suggesting that neither the smaller averaging volume for the LES quantities nor the increased accuracy of the spatial derivatives of the coefficients in Eq. (12) yields improved strain–rotation relations in HIT. However, the number of sub-domains becomes important if an inhomogeneous flow is considered, since the sub-domains have to be small enough to resolve the spatial changes of the mean or larger structures.

The next investigated model parameter is the estimated subgrid-scale kinetic energy. The effect of multiplying the estimated kinetic energy with a factor of 1.5 and 2.0 on the PDF of the second invariant of the velocity gradient tensor is shown in Fig. 10(b). A higher kinetic energy leads to wider tails of the PDF, because larger absolute velocity gradients become more likely. The shape of the PDF does not notably change with different kinetic energies. The lack of asymmetry as observed with the standard model setup still remains.

Another investigated parameter is the number of modes in the series expansion of the subgrid-scale velocity in Eq. (12). Its influence is shown in Fig. 10(c). For this parameter, no clear trend can be observed in the PDF. The second invariant of the velocity gradient tensor is relatively independent of the number of modes of the subgrid-scale velocity, for $Nm\u2265100$. Since the non-linear term is not present in the model equations for the subgrid-scale velocity, Eqs. (13) and (14), there is no information exchange between the coefficients of different wave numbers. This has an advantageous side effect for the implementation. The parallelization of the solution procedure for Eqs. (29), (30), (17), and (18) can be carried out straightforwardly in wave number space. However, it is expected that, in order to properly represent the kinetic energy spectrum and achieve a suitable complexity of the sampled subgrid-scale velocity field, a minimal amount of wave numbers is required.

The last model parameter, which is investigated with respect to the PDF of the second invariant of the velocity gradient tensor, is the constant in the expression of the subgrid-scale turbulent viscosity given in Eq. (11), $C\nu $. The PDFs of the second invariant of the velocity gradient tensor are shown for different values of $C\nu $ in Fig. 10(d). Higher values of this constant lead to slightly wider PDFs compared to smaller values. This can be explained by the fact that higher turbulent viscosities marginally increase the kinetic energy for the higher wave numbers relative to the small wave numbers, which typically yields larger velocity gradient magnitudes and hence wider tails of the PDF. The shape of the PDF is not affected. Even though the constant is varied by a factor of 16, the difference between the PDFs is very small. Thus, the results are not very sensitive to the choice of $C\nu $.

The specific choice of parameters in the model depends on the application of the enrichment with subgrid-scale velocity. General guidelines for the choice of parameter are thus not possible. However, it is shown that the subgrid-scale velocity statistics are rather insensitive to the change of parameters.

In order to verify that the interpolation proposed in Sec. II E does not significantly influence the results of this paper, the PDFs of *Q* with and without the proposed divergence-free interpolation between sub-domains is shown in Fig. 11. The PDFs almost coincide, suggesting that the interpolation does not affect the turbulence statistics on a global scale. The fact that the desirable properties are conserved and the resulting velocity field is divergence free makes the combination of the enrichment with the proposed interpolation between sub-domains advisable.

For the application of the model, it is also interesting to measure the increase in computational time due to the enrichment compared to the LES. In Table III, the computational times for the DNS, the LES, and the enriched LES with the standard model setup are compared with each other. The computational time is shown that it is needed to simulate a physical time of *T _{ref}*. The enrichment approximately doubles the required computational time of the LES. Compared to the computational time required for the DNS, this is an acceptable increase in computational time. The computational time of the DNS is approximately four orders of magnitude larger compared to the enriched LES. The presented statistics of the enriched LES are much closer to the DNS compared to the LES by only increasing the computational costs by a factor of two.

### B. Enriched velocity of the turbulent shear flow case

In order to verify the ability of the proposed enrichment model at simulating inhomogeneous and anisotropic flows, the enrichment is applied to the turbulent shear flow configuration that is described in Sec. III B 2. The fact that a distinct set of coefficients exists for every sub-domain makes the model capable of representing inhomogeneous subgrid-scale velocity fields. Furthermore, anisotropy of the subgrid-scale velocity [i.e., preferential alignment of the coefficients in the series expansion of Eq. (12)] may be induced by the sweeping and straining terms in the governing equations for the subgrid scales.

Due to the chaotic fluctuations of the kinetic energy of the shear flow, temporal averaging has to be performed over long times. The statistics are evaluated within a time span of $140Tshear$, starting from a flow field that converged toward a statistically steady state for at least the same time, in order to obtain meaningful statistics.

The generated subgrid-scale velocity is investigated in terms of its anisotropy and spatial variation. The results are assessed based on the averaged fluctuation velocity correlations

where the ensemble average $\u27e8.\u27e9$ in this case stands for averaging in time and along the homogeneous directions (i.e., x and z directions). Similarly, the filtered velocity $u\u0303$ or the subgrid-scale velocity $u\u2032$ can be used in Eq. (55) in order to obtain the velocity correlations for the filtered velocity and the subgrid-scale velocity, respectively.

In Fig. 12, the velocity correlations for the DNS and the LES [Fig. 12(a)] as well as for the subgrid-scale velocity obtained with the model [Fig. 12(b)] are depicted. Similar to the standard model setup of Sec. IV A, the model equations are solved for *N _{m}* = 100 modes, with the analytical value of the constant in the subgrid-scale turbulent viscosity $C\nu =0.4$ and a subgrid-scale kinetic energy estimation with the model of Yoshizawa and the analytical constant $CI=0.0826$. However, the number of sub-domains is 16

^{3}, which corresponds to the LES grid. Since the model equations are solved based on the domain-averaged filtered velocity, the average modeled subgrid-scale velocity is approximately constant in a domain, which is why a relatively high number of sub-domains is required to resolve the mean velocity profile.

The profiles of the correlations $C11,\u2009C22$, and $C33$ in Fig. 12(a) exhibit the same qualitative behavior. Over one domain length two periods of the correlations can be observed, which is the double frequency of the mean velocity profile in x direction. The maximum values of the correlations are in the regions of the highest mean shear. The cross-correlation $C12$ has the same frequency as the mean velocity in x direction. It is negative for positive values of $\u2202u1/\u2202y$ and positive for negative values of $\u2202u1/\u2202y$. The other two cross-correlations are zero, which is not shown in Fig. 12(a).

The velocity correlations obtained with the LES clearly deviate from the velocity correlations of the DNS, especially for $C11$ and $C22$. This can be explained with the too high amplitudes of the mean velocity profile. Adding the subgrid-scale velocity to the LES velocity does not correct the observed deviations. Compared to the total kinetic energy, the kinetic energy of the subgrid scales is very small in this configuration. However, note that the subgrid-scale velocity possesses properties that are important for some applications (e.g., increased turbulent mixing) and omitting them can lead to wrong results.

As it can be observed in Fig. 12(b), the generated subgrid-scale velocity field possesses a similar qualitative behavior as the velocity field of the DNS. Note that the magnitude of the velocity correlations of the subgrid-scale velocity is smaller, since the kinetic energy of the subgrid scale is only a small fraction of the total kinetic energy. Since the correlations vary in the y direction, an inhomogeneous field is generated. Furthermore, the regions of high and low values of the correlations coincide with those of the DNS. The reasons for the stepwise increase and decrease in the correlations are caused by the statistically homogeneous subgrid-scale velocity in the sub-domains. The ability of the model to generate an inhomogeneous subgrid-scale velocity field is a prerequisite to predict turbophoresis in a particle-laden flow.

In addition to the inhomogeneity of the subgrid-scale velocity, there is also evidence that the velocity field is anisotropic. In particular, the cross-correlation $C12$ demonstrates the capability of the model to produce an anisotropic velocity field. It has a similar trend as the cross-correlation of the DNS, whereas the other two cross-correlations $C13$ and $C23$ are zero (not shown in the figure).

It can be concluded that the concept of sub-domains indeed enables inhomogeneous subgrid-scale velocity fields. Moreover, the spatial change of the temporally averaged subgrid-scale velocity correlations behaves qualitatively the right way. Furthermore, anisotropy of the subgrid-scale velocity is obtained, because the correlations differ in function of the used subgrid-scale velocity components.

## V. CONCLUSIONS

A new model for enriching the velocity of a LES with a subgrid-scale velocity is introduced. The enrichment model enables closures for LES when these are required. It is shown that the generated subgrid-scale velocity possesses several important characteristic properties of real turbulent fields, such as the shape of the kinetic energy spectrum, non-Gaussian PDFs of velocity gradients, and similar strain–rotation relations. Compared to other models for the subgrid-scale velocity, common problems are mitigated with our proposed model: the influence of the small number of tuning parameters required in the model is shown to be insignificant, the computational cost of the enrichment model is of the order of the LES, and it is shown that inhomogeneous and anisotropic subgrid-scale velocities can also be generated. These properties make our enrichment model suitable for a variety of applications, in which a closure for the subgrid-scale velocity is required.

## ACKNOWLEDGMENTS

This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID No. 457509672.

F. Evrard has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 101026017.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Max Hausmann:** Conceptualization (equal); Methodology (lead); Writing – original draft (lead). **Fabien Evrard:** Conceptualization (equal); Funding acquisition (supporting); Supervision (equal); Writing – review & editing (equal). **Berend van Wachem:** Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Resources (lead); Software (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

### APPENDIX: PROJECTION OF THE MODEL EQUATIONS

Equation (14) can be written in a short form

using $A,\u2009D$, and $F$ for the convective term, diffusive term, and source term, respectively. Neglecting the pressure term and approximating the time derivative with an explicit Euler scheme yields a preliminary subgrid-scale velocity $ui\u2032,*$

The same equation but with the pressure term leads to the physically right (up to discretization error of the time derivative) subgrid-scale velocity of the next time level $ui\u2032,n+1$

Assuming the divergence of the subgrid-scale velocity of the next time level to vanish, the following Poisson equation can be derived:

Similar to the subgrid-scale velocity, the pressure can be written as a truncated Fourier series

with the gradient

and the Laplacian

Plugging the Fourier series expansions into Eq. (A5) yields

After sorting by $cos$- and $sin$-terms and evaluating each mode separately, the following relations for the coefficients $Pm(t)$ and $Qm(t)$ are obtained:

Finally, the divergence-free coefficients for the subgrid-scale velocity $Amn+1(t)$ and $Bmn+1(t)$ can be computing using Eq. (A4)

## References

_{4}-H

_{2}fuel blends