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.

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 Meneveau9 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 Lele24,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.

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.

The foundation of LES and of the following discussions is that any physical quantity of interest, ϕ (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, ϕ̃, may be defined by the following spatial low-pass filtering operation, as proposed by Leonard:26 

ϕ̃(x)=G(xξ)ϕ(ξ)dξ,
(1)

where G denotes a filter kernel. This filtering kernel determines the shape of the spectrum of the resulting filtered quantity. The subfilter-scale contribution, ϕ, is the difference between the original quantity, ϕ, and the filtered quantity, ϕ̃. 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 ũ and u 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 νf are given by

uixi=0,
(2)
uit+ujuixj=1ρfpxi+νf2uixjxj+si,
(3)

where ui is the velocity, p the pressure, and si 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 G,

ũixi=0,
(4)
ũit+ũjũixj=1ρfp̃xi+νf2ũixjxjτijxj+s̃i.
(5)

The additional term τij/xj 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 ϕ=ϕ̃+ϕ to become

τij=ũiũj̃+uiũj̃+ũiuj̃+uiuj̃ũiũj.
(6)

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

uixi=0,
(7)
uit+ujuixj+ũjuixj+ujũixj=1ρfpxi+νf2uixjxj+τijxj+si.
(8)

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(ΩijΩijSijSij), or the angle between the eigenvectors of the rotation-rate tensor Ωij and the strain-rate tensor Sij 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 ujui/xj redistributes kinetic energy in spectral space, whereas a diffusive term 2ui/(xjxj) exclusively removes kinetic energy if the (turbulent) viscosity is positive (see, e.g., Pope30). 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.

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

τij,mod=2νtS̃ij+13τkk,modδij,
(9)

where νt is the model-specific turbulent viscosity, δij is the Kronecker tensor, and S̃ij is the strain-rate tensor of the large-scale velocity, given by

S̃ij=12(ũixj+ũjxi).
(10)

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 νt. Based on this work, the following expression for the turbulent viscosity is used, which is derived from renormalization groups to be37 

νt(k)=(νf2Cνkq2E(q)dq)1/2νf,
(11)

where k is the wave number, E(k) is the kinetic energy spectrum, and Cν is a constant with an analytical value of Cν=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 Yoshizawa39 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.

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 can be represented by the finite sum

u=m=0Nm1(Am(t)cos(km·x)+Bm(t)sin(km·x)),
(12)

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 Nm 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·Am(t)=km·Bm(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

uixi=0,
(13)
uit+ũjuixj+ujũixj=1ρfpxi+(νf+νt)2uixjxj+fi,
(14)

where fi 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:

Am,i*(t)Am,in(t)Δt+ũjnkm,jBm,in(t)+Am,jn(t)ũinxj=(νf+νt)|km|2Am,in(t)+fm,i,
(15)
Bm,i*(t)Bm,in(t)Δtũjnkm,jAm,in(t)+Bm,jn(t)ũinxj=(νf+νt)|km|2Bm,in(t)+gm,i,
(16)

where Δ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:

Amn+1(t)=Am*(t)kmkm·Am*(t)|km|2,
(17)
Bmn+1(t)=Bm*(t)kmkm·Bm*(t)|km|2.
(18)

The computational cost of the solution of one time step consisting of Eqs (15)–(18) scales linear with the number of modes Nm. 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

Vm,0=(sin(θm)cos(ϕm),sin(θm)sin(ϕm),cos(θm))T
(19)

with

ϕm[0,2π],cos(θm)[1,1].
(20)

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

Vmn=αmVmn1+βmVm,0,
(21)

where Vm,0 is a random vector that is newly generated for every wave number and every time step, and Vmn1 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·Vmn0). The weights αm and βm determine how fast the directions change in time, and are computed as follows:

αm=exp(ΔtTE,m),βm=1αm2.
(22)

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

E(|km|)=14Δkm(|Am|2+|Bm|2),
(23)

where Δkm is the distance between consecutive wave numbers for the wave number km. Similar to the forcing scheme of Mallouppas et al.,38 the source term in the subgrid-scale equation is

fi,m=1ΔtKdesiredKactualKdesiredvi,mtrigger,
(24)

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

vi,mtrigger=KdesiredKsourceVA,mkm5/6.
(25)

The forcing term is zero if the kinetic energy of the subgrid-scale velocity field Kactual is equal the desired kinetic energy Kdesired. 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

Ksource=14m=0Nm1(|VA,mkm5/6|2+|VB,mkm5/6|2),
(26)

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

Kactual=14m=0Nm1(|Am|2+|Bm|2).
(27)

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 Yoshizawa39 

Kdesired=CIΔ22S̃ijS̃ij
(28)

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 

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 ·domain. In particular, this applies to the LES velocity ũidomain, the gradient of the LES velocity ũi/xjdomain, and theoretically also the estimation of the kinetic energy of the subgrid-scale velocity Kdesireddomain. 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

Am,i*Am,inΔt+ũjndomain(km,jBm,in+Am,inxj)+Am,jnũinxjdomain=(νf+νt)(|km|2Am,in+2Am,inxjxj+2km,jBm,inxj)+fm,i,
(29)
Bm,i*Bm,inΔt+ũjndomain(Bm,inxjkm,jAm,in)+Bm,jnũinxjdomain=(νf+νt)(|km|2Bm,in+2Bm,inxjxj2km,jAm,inxj)+gm,i.
(30)

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.

FIG. 1.

General enrichment strategy consisting of three steps. (I) Information on the large-scale flow dynamics is obtained by averaging within sub-domains that are coarser than the LES grid. (II) The model equations (29), (30), (17), and (18) are solved in every sub-domain. (III) The enriched velocity is approximated as the sum of the resolved LES velocity and the modeled subgrid-scale velocity.

FIG. 1.

General enrichment strategy consisting of three steps. (I) Information on the large-scale flow dynamics is obtained by averaging within sub-domains that are coarser than the LES grid. (II) The model equations (29), (30), (17), and (18) are solved in every sub-domain. (III) The enriched velocity is approximated as the sum of the resolved LES velocity and the modeled subgrid-scale velocity.

Close modal

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 consists of a contribution of the lower sub-domain ulow, the upper sub-domain uup, and a correction velocity ucorr,

uint=W(x)ulow+(1W(x))uup+ucorr,
(31)

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

·uint=0.
(32)

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

u=m=0Nm1u,m=m=0Nm1(Am(t)cos(km·x)+Bm(t)sin(km·x)),
(33)

where u,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,m is zero

·u,m=0.
(34)

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

dW(x)dx(ulow,muup,m)+vcorr,my+wcorr,mz=0.
(35)

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,m=0, which leads to vcorr,m=Ψm(x)/km,2 or set vcorr,m=0, which leads to wcorr,m=Ψm(x)/km,3. In order to keep the correction velocity magnitude small, min(|vcorr,m|,|wcorr,m|) decides which case is used. The correction velocity kernel is defined as

Ψ(x)=W(x)x[(Alow,m,1Aup,m,1)sin(km·x)(Blow,m,1Bup,m,1)cos(km·x)].
(36)

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,

W(r)=111+eαr,
(37)

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 Tref, 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.

FIG. 2.

Normalized divergence of the modeled subgrid-scale velocity computed using central differences. (a) is obtained without interpolation between the sub-domains and (b) is obtained using the proposed interpolation.

FIG. 2.

Normalized divergence of the modeled subgrid-scale velocity computed using central differences. (a) is obtained without interpolation between the sub-domains and (b) is obtained using the proposed interpolation.

Close modal

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.

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 Kest 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 23 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.

FIG. 3.

PDF of the centered and normalized estimated subgrid-scale kinetic energy. We show the difference between averaging the kinetic energy in multiple LES grid cells (average) and evaluating it for every LES grid cells separately (single). In order to be able to assess the resulting distributions, a squared Gaussian distribution is plotted.

FIG. 3.

PDF of the centered and normalized estimated subgrid-scale kinetic energy. We show the difference between averaging the kinetic energy in multiple LES grid cells (average) and evaluating it for every LES grid cells separately (single). In order to be able to assess the resulting distributions, a squared Gaussian distribution is plotted.

Close modal

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.

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.

First, the numerical solution of the NSE and LES equations is briefly summarized in Sec. III A. Then, that the exact configurations of the two test cases are explained in Sec. III B.

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:

ũixi=0,
(38)
ũit+ũjũixj=1ρfPxi+xj((νf+νt)2(ũixj+ũjxj))+Fi.
(39)

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 model33 

νt=(CSΔ)2|S̃ij|2,
(40)

where CS is the Smagorinsky constant, Δ is the filter width, and |S̃ij| is the following norm of the filtered strain rate tensor:

|S̃ij|=2S̃ijS̃ij.
(41)

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

νt=CVBβÃijÃij
(42)

with the filtered velocity gradient tensor

Ãij=ũixj
(43)

and

βij=Δ2ÃikÃjk,
(44)
Bβ=β11β22β122+β11β33β132+β22β33β232.
(45)

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

Fi=1ΔtKdesiredKactualKdesiredvitrigger.
(46)

However, there are two differences. Instead of estimating the desired kinetic energy, Kdesired, 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:

vitrigger=mUm,icos(km,jxj+ψm).
(47)

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ψ2π. The velocity coefficients Um,i are chosen to lead to a divergence-free trigger velocity and thus fulfills the condition Um·km=0. Its magnitude corresponds to a von Kármán kinetic energy spectrum in the triggered wave number range and is given by

|Um|=2Evk(|km|)Δk,
(48)

where Evk is a von Kármán spectrum and Δ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.

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 kstart and kend. Physical parameters of the resulting turbulent velocity field are the Taylor–Reynolds number Reλ=λurms/νf, the turbulent Reynolds number Rel=l11urms/νf, the Kolmogorov length scale η, the Taylor microscale λ, and the longitudinal (l11) and transverse (l22, l33) 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.

TABLE I.

Parameters of the HIT simulation configuration.

ParameterValue
Reλ 75 
Rel 205 
η/L 0.0017 
τη/Tref 0.0075 
λ/L 0.029 
l11/L 0.079 
l22/L 0.49 
l33/L 0.38 
kstartL/2π 
kendL/2π 
ParameterValue
Reλ 75 
Rel 205 
η/L 0.0017 
τη/Tref 0.0075 
λ/L 0.029 
l11/L 0.079 
l22/L 0.49 
l33/L 0.38 
kstartL/2π 
kendL/2π 

The solution of the DNS is obtained on a numerical grid with 2563 cells. This leads to a product of the largest resolved wave number kmax and the Kolmogorov length scale of kmaxη=1.37.

The LES of the same case is solved on a grid consisting of 323 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|5/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.

FIG. 4.

Kinetic energy spectrum of forced HIT. We compare the spectra of the HIT case of a DNS with 2563 grid cells and a LES with 323 grid cells.

FIG. 4.

Kinetic energy spectrum of forced HIT. We compare the spectra of the HIT case of a DNS with 2563 grid cells and a LES with 323 grid cells.

Close modal

For a homogeneous isotropic flow, the spatial autocorrelations Bαα(1) are defined as45 

Bαα(1)(r)=uα(x,t)uα(x+re1,t)uα(x,t)uα(x,t),
(49)

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αα are defined as

lαα=0Bαα(1)(r)dr.
(50)

Note that no implicit summation is carried over double occurring indices α. In isotropic turbulence, the transverse integral length scale is half the longitudinal integral length scale.30 It can be seen in Table I that this relation approximately holds for the generated turbulent field.

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

s=smaxsin(2πy/L)e1.
(51)
FIG. 5.

Sketch of the shear flow configuration. The momentum source (black line) varies in the shape of a sin-profile in y direction and acts in the x direction. The cubic domain is periodic in all directions. The colors show the velocity in x direction.

FIG. 5.

Sketch of the shear flow configuration. The momentum source (black line) varies in the shape of a sin-profile in y direction and acts in the x direction. The cubic domain is periodic in all directions. The colors show the velocity in x direction.

Close modal

The amplitude is given by smax 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/ν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.

TABLE II.

Parameters of the turbulent shear flow simulation configuration.

ParameterValue
Reshear 3115 
η/L 0.0029 
τη/Tshear 0.026 
ParameterValue
Reshear 3115 
η/L 0.0029 
τη/Tshear 0.026 

A DNS with 1283 grid cells fully resolves the length scales down to the Kolmogorov scale of the shear flow (kmaxη=1.16). In addition to the DNS, a LES with 163 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.

FIG. 6.

(a) Comparison of the mean velocity profiles over the y-coordinate in the turbulent shear flow of the DNS (solid line) and the LES with the Vreman model (dashed line) and (b) dimensionless kinetic energy over dimensionless time for the DNS and the LES of the turbulent shear flow. The dashed lines are the average kinetic energies.

FIG. 6.

(a) Comparison of the mean velocity profiles over the y-coordinate in the turbulent shear flow of the DNS (solid line) and the LES with the Vreman model (dashed line) and (b) dimensionless kinetic energy over dimensionless time for the DNS and the LES of the turbulent shear flow. The dashed lines are the average kinetic energies.

Close modal

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.

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.

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 Nm = 100 and the constant in the turbulent viscosity in the enrichment model νt is Cν=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 Yoshizawa39Kest=CIΔ22S̃ijS̃ij, 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.

FIG. 7.

Slices in the x–y plane of the velocity magnitudes of (a) the DNS velocity, (b) the LES velocity without enrichment, and (c) the LES velocity enriched with the proposed model for the HIT case.

FIG. 7.

Slices in the x–y plane of the velocity magnitudes of (a) the DNS velocity, (b) the LES velocity without enrichment, and (c) the LES velocity enriched with the proposed model for the HIT case.

Close modal

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).

FIG. 8.

Kinetic energy spectrum of forced HIT. We compare the results of the DNS with 2563 grid cells, the LES with 323 grid cells, and the modeled subgrid-scale velocity field.

FIG. 8.

Kinetic energy spectrum of forced HIT. We compare the results of the DNS with 2563 grid cells, the LES with 323 grid cells, and the modeled subgrid-scale velocity field.

Close modal

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.41Figure 9 shows the PDFs of the longitudinal and transverse velocity gradients A11 and A12 normalized with their standard deviations σA11 and σ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.

FIG. 9.

PDF of longitudinal (a) and transverse (b) velocity gradient for DNS, LES without enrichment, and enriched LES (LES-E) velocity field normalized with their standard deviations together with a Gaussian distribution for comparison. The results of the enriched LES are obtained with the standard model setup.

FIG. 9.

PDF of longitudinal (a) and transverse (b) velocity gradient for DNS, LES without enrichment, and enriched LES (LES-E) velocity field normalized with their standard deviations together with a Gaussian distribution for comparison. The results of the enriched LES are obtained with the standard model setup.

Close modal

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

Q=12(ΩijΩijSijSij),
(52)

where Ωij is the rotation-rate tensor and Sij 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ν in the turbulent viscosity νt of the subgrid-scale momentum equation.

FIG. 10.

PDF of the second invariant of the velocity gradient tensor for DNS, LES without enrichment, and the enriched LES (LES-E) for the HIT case. The influence of four parameters of the proposed enrichment model is compared to the previously defined standard model setup: (a) the number of sub-domains, (b) the estimated subgrid-scale kinetic energy, (c) the number of modes in the series expansion (12), and (d) the constant Cν in the subgrid-scale turbulent viscosity.

FIG. 10.

PDF of the second invariant of the velocity gradient tensor for DNS, LES without enrichment, and the enriched LES (LES-E) for the HIT case. The influence of four parameters of the proposed enrichment model is compared to the previously defined standard model setup: (a) the number of sub-domains, (b) the estimated subgrid-scale kinetic energy, (c) the number of modes in the series expansion (12), and (d) the constant Cν in the subgrid-scale turbulent viscosity.

Close modal

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:

Q=12(ΩijΩijSijSij)=12uixjujxi=12xi(ujuixj).
(53)

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

Q=xixj(ũiũj+ũiuj+uiũjnonlocalcontributions+uiuj).
(54)

In contrast to the local contributions ũiũj and uiuj, 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 83 sub-domains of the standard model setup, two further simulations are performed with 163 and 323 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 Nm100. 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ν. The PDFs of the second invariant of the velocity gradient tensor are shown for different values of Cν 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ν.

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.

FIG. 11.

PDF of the second invariant of the velocity gradient tensor of the enriched velocity of the HIT case. The velocity field of the enrichment with the standard model setup is used with discontinuities between the sub-domains and with the proposed divergence-free interpolation.

FIG. 11.

PDF of the second invariant of the velocity gradient tensor of the enriched velocity of the HIT case. The velocity field of the enrichment with the standard model setup is used with discontinuities between the sub-domains and with the proposed divergence-free interpolation.

Close modal

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 Tref. 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.

TABLE III.

Comparison of the computational times for the DNS, the LES, and the enriched LES for the simulation time Tref.

TypeComputational time per Tref
LES 32 × 103
LES + model (standard setup) 69 × 103
DNS 340 × 106
TypeComputational time per Tref
LES 32 × 103
LES + model (standard setup) 69 × 103
DNS 340 × 106

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

Cij=(uiui)(ujuj),
(55)

where the ensemble average . in this case stands for averaging in time and along the homogeneous directions (i.e., x and z directions). Similarly, the filtered velocity ũ or the subgrid-scale velocity u 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 Nm = 100 modes, with the analytical value of the constant in the subgrid-scale turbulent viscosity Cν=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 163, 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.

FIG. 12.

Dimensionless velocity correlations of the turbulent shear flow in y direction. (a) Velocity correlations of the DNS (solid lines) and the LES (dashed lines). (b) Velocity correlations of the generated subgrid-scale velocity.

FIG. 12.

Dimensionless velocity correlations of the turbulent shear flow in y direction. (a) Velocity correlations of the DNS (solid lines) and the LES (dashed lines). (b) Velocity correlations of the generated subgrid-scale velocity.

Close modal

The profiles of the correlations C11,C22, 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 u1/y and positive for negative values of u1/y. 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.

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.

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.

The authors have no conflicts to disclose.

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).

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

Equation (14) can be written in a short form

uit+Ai=1ρfpxi+Di+Fi
(A1)

using A,D, 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,*

ui,*ui,nΔt+Ain=Din+Fin.
(A2)

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,n+1

ui,n+1ui,nΔt+Ain=1ρfpxi+Din+Fin.
(A3)

The difference of ui,n+1 and ui,* is obtained by subtracting Eq. (A2) from (A3)

ui,n+1ui,*Δt=1ρfpxi.
(A4)

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

ui,*xi=Δtρf2pxixi.
(A5)

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

p=m=0Nm1(Pm(t)cos(km·x)+Qm(t)sin(km·x))
(A6)

with the gradient

p=m=0Nm1km(Pm(t)sin(km·x)+Qm(t)cos(km·x))
(A7)

and the Laplacian

2p=m=0Nm1|km|2(Pm(t)cos(km·x)+Qm(t)sin(km·x)).
(A8)

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

m=0Nm1km(Am*(t)sin(km·x)+Bm*(t)cos(km·x))=Δtρfm=0Nm1|km|2(Pm(t)cos(km·x)+Qm(t)sin(km·x)).
(A9)

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:

Pm(t)=ρfΔt|km|2km·Bm*(t),
(A10)
Qm(t)=ρfΔt|km|2km·Am*(t).
(A11)

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

Amn+1(t)=Am*(t)ΔtρfQm(t)km=Am*(t)kmkm·Am*(t)|km|2,
(A12)
Bmn+1(t)=Bm*(t)+ΔtρfPm(t)km=Bm*(t)kmkm·Bm*(t)|km|2.
(A13)
1.
P.
Fede
and
O.
Simonin
, “
Numerical study of the subgrid fluid turbulence effects on the statistics of heavy colliding particles
,”
Phys. Fluids
18
,
045103
(
2006
).
2.
B.
Ray
and
L. R.
Collins
, “
Preferential concentration and relative velocity statistics of inertial particles in Navier-Stokes turbulence with and without filtering
,”
J. Fluid Mech.
680
,
488
510
(
2011
).
3.
C.
Marchioli
, “
Large-eddy simulation of turbulent dispersed flows: A review of modelling approaches
,”
Acta Mech.
228
,
741
771
(
2017
).
4.
P. L.
Johnson
and
C.
Meneveau
, “
Predicting viscous-range velocity gradient dynamics in large-eddy simulations of turbulence
,”
J. Fluid Mech.
837
,
80
114
(
2018
).
5.
M.
Herrmann
, “
A sub-grid surface dynamics model for sub-filter surface tension induced interface dynamics
,”
Comput. Fluids
87
,
92
101
(
2013
).
6.
D.
Kedelty
,
J.
Uglietta
, and
M.
Herrmann
, “
A volume-of-fluid dual-scale approach for simulating turbulent liquid/gas interactions
,” in
Turbulence and Interactions
, edited by
M.
Deville
,
C.
Calvin
,
V.
Couaillier
,
M.
De La Llave Plata
,
J.-L.
Estivalèzes
,
T. H.
, and
S.
Vincent
(
Springer International Publishing
,
Cham
,
2021
), Vol.
149
, pp.
39
51
, Series Title: Notes on Numerical Fluid Mechanics and Multidisciplinary Design.
7.
Z.
Zhang
,
P.
Huang
,
G.
Bitsuamlak
, and
S.
Cao
, “
Large-eddy simulation of wind-turbine wakes over two-dimensional hills
,”
Phys. Fluids
34
,
065123
(
2022
).
8.
A. M.
Garcia
,
S. L.
Bras
,
J.
Prager
,
M.
Häringer
, and
W.
Polifke
, “
Large eddy simulation of the dynamics of lean premixed flames using global reaction mechanisms calibrated for CH4-H2 fuel blends
,”
Phys. Fluids
34
,
095105
(
2022
).
9.
A.
Scotti
and
C.
Meneveau
, “
A fractal model for large eddy simulation of turbulent flow
,”
Physica D
127
,
198
232
(
1999
).
10.
N. A.
Malik
and
J. C.
Vassilicos
, “
A Lagrangian model for turbulent dispersion with turbulent-like flow structure: Comparison with direct numerical simulation for two-particle statistics flow structure
,”
Phys. Fluids
11
,
1572
(
1999
).
11.
J. C. H.
Fung
,
J. C. R.
Hunt
,
N. A.
Malik
, and
R. J.
Perkins
, “
Kinematic simulation of homogeneous turbulence by unsteady random Fourier modes
,”
J. Fluid Mech.
236
,
281
318
(
1992
).
12.
R. H.
Kraichnan
, “
Diffusion by a random velocity field
,”
Phys. Fluids
13
,
22
(
1970
).
13.
B.
Ray
and
L. R.
Collins
, “
A subgrid model for clustering of high-inertia particles in large-eddy simulations of turbulence
,”
J. Turbul.
15
,
366
385
(
2014
).
14.
Z.
Zhou
,
S.
Wang
, and
G.
Jin
, “
A structural subgrid-scale model for relative dispersion in large-eddy simulation of isotropic turbulent flows by coupling kinematic simulation with approximate deconvolution method
,”
Phys. Fluids
30
,
105110
(
2018
).
15.
J. A.
Domaradzki
and
K.-C.
Loh
, “
The subgrid-scale estimation model in the physical space representation
,”
Phys. Fluids
11
,
2330
2342
(
1999
).
16.
J. A.
Domaradzki
and
E. M.
Saiki
, “
A subgrid-scale model based on the estimation of unresolved scales of turbulence
,”
Phys. Fluids
9
,
2148
2164
(
1997
).
17.
G. I.
Park
,
M.
Bassenne
,
J.
Urzay
, and
P.
Moin
, “
A simple dynamic subgrid-scale model for LES of particle-laden turbulence
,”
Phys. Rev. Fluids
2
,
044301
(
2017
).
18.
M.
Bassenne
,
M.
Esmaily
,
D.
Livescu
,
P.
Moin
, and
J.
Urzay
, “
A dynamic spectrally enriched subgrid-scale model for preferential concentration in particle-laden turbulence
,”
Int. J. Multiphase Flow
116
,
270
280
(
2019
).
19.
T. J.
Hughes
,
G. R.
Feijóo
,
L.
Mazzei
, and
J.-B.
Quincy
, “
The variational multiscale method—A paradigm for computational mechanics
,”
Comput. Methods Appl. Mech. Eng.
166
,
3
24
(
1998
).
20.
T. J.
Hughes
,
L.
Mazzei
, and
K. E.
Jansen
, “
Large eddy simulation and the variational multiscale method
,”
Comput. Visualization Sci.
3
,
47
59
(
2000
).
21.
M.
Terracol
,
P.
Sagaut
, and
C.
Basdevant
, “
A multilevel algorithm for large-eddy simulation of turbulent compressible flows
,”
J. Comput. Phys.
167
,
439
474
(
2001
).
22.
B.
Dubrulle
,
J. P.
Laval
,
S.
Nazarenko
, and
N. K.-R.
Kevlahan
, “
A dynamic subfilter-scale model for plane parallel flows
,”
Phys. Fluids
13
,
2045
2064
(
2001
).
23.
J.-P.
Laval
,
B.
Dubrulle
, and
S.
Nazarenko
, “
Fast numerical simulations of 2D turbulence using a dynamic model for subfilter motions
,”
J. Comput. Phys.
196
,
184
207
(
2004
).
24.
A. S.
Ghate
and
S. K.
Lele
, “
Subfilter-scale enrichment of planetary boundary layer large eddy simulation using discrete Fourier-Gabor modes
,”
J. Fluid Mech.
819
,
494
539
(
2017
).
25.
A. S.
Ghate
and
S. K.
Lele
, “
Gabor mode enrichment in large eddy simulations of turbulent flow
,”
J. Fluid Mech.
903
,
A13
(
2020
).
26.
A.
Leonard
, “
Energy cascade in large-eddy simulations of turbulent fluid flows
,” in
Advances in Geophysics
(
Elsevier
,
1975
), Vol.
18
, pp.
237
248
.
27.
P.
Sagaut
,
Large Eddy Simulation for Incompressible Flows
, 3rd ed. (
Springer
,
2005
).
28.
J.-P.
Laval
,
B.
Dubrulle
, and
S.
Nazarenko
, “
Nonlocality and intermittency in three-dimensional turbulence
,”
Phys. Fluids
13
,
1995
2012
(
2001
).
29.
J. C. R.
Hunt
and
D. J.
Carruthers
, “
Rapid distortion theory and the ‘problems’ of turbulence
,”
J. Fluid Mech.
212
,
497
(
1990
).
30.
S. B.
Pope
,
Turbulent Flows
, 6th ed. (
Cambridge University Press
,
Cambridge
,
2000
).
31.
C.
Cambon
and
J. F.
Scott
, “
Linear and nonlinear models of anisotropic turbulence
,”
Annu. Rev. Fluid Mech.
31
,
1
53
(
1999
).
32.
K.
Horiuti
, “
Roles of non-aligned eigenvectors of strain-rate and subgrid-scale stress tensors in turbulence generation
,”
J. Fluid Mech.
491
,
65
100
(
2003
).
33.
J.
Smagorinsky
, “
General circulation experiments with the primitive equations. I. The basic experiment
,”
Mon. Weather Rev.
91
,
99
165
(
1963
).
34.
M.
Germano
,
U.
Piomelli
,
P.
Moin
, and
W. H.
Cabot
, “
A dynamic subgrid-scale eddy viscosity model
,”
Phys. Fluids A
3
,
1760
1765
(
1991
).
35.
F.
Nicoud
and
F.
Ducros
, “
Subgrid-scale stress modelling based on the square of the velocity gradient tensor
,”
Flow Turbul. Combust.
62
,
183
200
(
1999
).
36.
A. W.
Vreman
, “
An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications
,”
Phys. Fluids
16
,
3670
3681
(
2004
).
37.
V. M.
Canuto
and
M. S.
Dubovikov
, “
A dynamical model for turbulence. I. General formalism
,”
Phys. Fluids
8
,
571
586
(
1996
).
38.
G.
Mallouppas
,
W. K.
George
, and
B.
van Wachem
, “
New forcing scheme to sustain particle-laden homogeneous and isotropic turbulence
,”
Phys. Fluids
25
,
083304
(
2013
).
39.
A.
Yoshizawa
, “
Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling
,”
Phys. Fluids
29
,
2152
(
1986
).
40.
P.
Moin
,
K.
Squires
,
W.
Cabot
, and
S.
Lee
, “
A dynamic subgrid-scale model for compressible turbulence and scalar transport
,”
Phys. Fluids A
3
,
2746
(
1991
).
41.
J.
Jimenez
, “
Intermittency in turbulence
,” in
Proceedings of the 15th “Aha Huliko” a Winter Workshop
(
2007
).
42.
G.
Mallouppas
, “
Modelling of turbulent gas-solid flows from DNS to LES
,” Ph.D. thesis (
Imperial College, London
,
2013
).
43.
F.
Denner
,
F.
Evrard
, and
B.
van Wachem
, “
Conservative finite-volume framework and pressure-based algorithm for flows of incompressible, ideal-gas and real-gas fluids at all speeds
,”
J. Comput. Phys.
409
,
109348
(
2020
).
44.
C. M.
Rhie
and
W. L.
Chow
, “
Numerical study of the turbulent flow past an airfoil with trailing edge separation
,”
AIAA J.
21
,
1525
1532
(
1983
).
45.
W. K.
George
,
Lectures in Turbulence for the 21st Century
(
Chalmers University of Technology
,
Gothenburg, Sweden
,
2013
), Vol. 550,
353
pp.
46.
E.
Lévêque
,
F.
Toschi
,
L.
Shao
, and
J.-P.
Bertoglio
, “
Shear-improved Smagorinsky model for large-eddy simulation of wall-bounded turbulent flows
,”
J. Fluid Mech.
570
,
491
502
(
2007
).
47.
S.
Balachandar
and
J. K.
Eaton
, “
Turbulent dispersed multiphase flow
,”
Annu. Rev. Fluid Mech.
42
,
111
133
(
2010
).
48.
M. R.
Maxey
, “
The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields
,”
J. Fluid Mech.
174
,
441
(
1987
).
49.
Z.
Zhou
,
S.
Wang
,
X.
Yang
, and
G.
Jin
, “
A structural subgrid-scale model for the collision-related statistics of inertial particles in large-eddy simulations of isotropic turbulent flows
,”
Phys. Fluids
32
,
095103
(
2020
).
50.
L.
Brandt
and
F.
Coletti
, “
Particle-laden turbulence: Progress and perspectives
,”
Annu. Rev. Fluid Mech.
54
,
159
189
(
2022
).