This paper evaluates the effect of forcing to sustain turbulence on the transfer function of the fluid with particles suspended in a homogeneous and isotropic flow. As mentioned by Lucci et al [“

Modulation of isotropic turbulence by particles of Taylor length-scale size
,” J. Fluid Mech.650, 5
55
(2010)], there are three limitations of forcing particle-laden homogeneous and isotropic turbulence: (a) large fluctuations on the temporal evolution of the kinetic energy are created when forcing is active at low wavenumbers, (b) the redistribution of energy is affected when forcing is performed over all wavenumbers, and (c) the nonlinear transfer function of the fluid due to the triadic interactions is affected when forcing is active over a wavenumber range. These limitations make the interpretation of the effects of particles on the energy spectrum of the fluid difficult. A new forcing scheme in physical space has been designed which avoids these limitations in wavenumber space, so the spectral effects of particles can be evaluated. The performance of this forcing scheme is tested using Direct Numerical Simulations. It is shown that the nonlinear transfer function of the fluid with the current forcing scheme is only affected at the wavenumbers it is acting, consistent with the theory. Even so, the spatial coherence and phase spectra between the two-way coupling and the fluid computed from the simulations show that the new forcing scheme is only moderately correlated even for the forced wavenumbers, with correlation coefficient typically about 10%.

Statistically, stationary Homogeneous Isotropic Turbulence (HIT) provides us with the opportunity to study turbulence in the absence of boundary conditions and turbulence transport. HIT is a valuable tool because it enables us to isolate and therefore study specific aspects of a turbulent flow. In addition, it can be used in the development of models for predicting the statistical characteristics of turbulence. It is important to note that HIT cannot be found naturally and it is even very difficult to be obtained via experiments. Nevertheless there are successful experiments such as the one performed by Comte-Bellot and Corrsin3 via the use of static grids. On the other hand, it is very easy to obtain a reasonable approximation to HIT via computer simulations, and since the pioneering study of Riley and Patterson13 has been a subject of study for over three decades. In the absence of any production mechanism, homogeneous and isotropic turbulence will decay due to the dissipative effect of viscosity. Therefore body forces, in terms of a forcing scheme, may be added in the momentum equation of the fluid if statistical stationarity is to be maintained.

Lundgren11 proposed a turbulence forcing scheme which balances the fluid dissipation rate and the turbulent kinetic energy in order to produce a stationary, homogeneous, and isotropic turbulent flow by adding a constant body force in physical space. This method triggers all wavenumber modes of the energy spectrum simultaneously. Other schemes avoid this by triggering specific wavenumbers, which are usually the low wavenumbers of the energy spectrum such as the one proposed by Eswaran and Pope.6 The method of Eswaran and Pope6 involves the use of a pseudo-spectral framework. The idea of their forcing scheme is to keep the total dissipation of the fluid constant by making use of statistically independent stochastic processes at each wavenumber. These random processes, which also satisfy continuity, supply energy to the first few wavenumber modes of the fluid by a prescribed variance and integral time scale of the turbulence.14 It is worthwhile mentioning that when energy is supplied at these first few low wavenumbers the temporal fluctuations are very high.9,14

Taking advantage of the retained statistics offered by forced homogeneous and isotropic turbulence, Squires and Eaton,16 Squires and Eaton,17 and Boivin et al.2 among others have performed studies involving two-way coupled forced turbulence in order to examine the dispersion of particles and the modulation of turbulence. Recently, it is argued by Lucci et al.,9 that existing forcing schemes may have a direct impact on the behaviour of the two-way coupling spectrum. Lucci et al.9 list three possible reasons that make the study of forced turbulence with current forcing schemes and two-way coupling “incorrect”: (a) the forcing of the first few low wavenumber modes in spectral space creates large fluctuations on the temporal evolution of the turbulence kinetic energy which can be of the same order of magnitude as the effect of the two-way coupling; (b) forcing performed in physical space, where all wavenumber modes are triggered as Lundgren11 proposes, distorts the redistribution of energy and hence the effects of particles on E(κ, t) cannot be quantified; and (c) due to the triadic interactions, it is not possible to investigate the modification of the nonlinear transfer function, T(κ, t), and consequently the fluid energy spectrum at the wavenumber modes where forcing is active. Their conclusions are later shown to be valid for the forcing scheme proposed by Lundgren,11 because the nonlinear transfer function is affected at all wavenumbers. In fact the nonlinear transfer function can be expressed as T(κ, t) = 2(νfκ2A)E(κ, t), where the factor A is the forcing constant. This constant therefore affects all the wavenumbers of the energy spectrum, but it is also shown that more general forcing schemes can avoid this difficulty, at least for the unforced wavenumbers.

The purpose of this paper is to present a new forcing scheme which avoids the three limitations suggested by Lucci et al.9 It avoids the influence of the forcing scheme on the nonlinear transfer function for all wavenumbers, but keeps the implementation in physical space. Using the new forcing scheme, Direct Numerical Simulations (DNS) of forced HIT on a 1283 periodic cubic box (of length Lb = 0.128 m) fully coupled with particles of Stokes numbers of 0.07 and 3.45 at a Taylor Reynolds number of 35.4 are performed, in order to evaluate the relative performance of the current forcing scheme both at a statistical and a theoretical level. To elucidate the effects of the new forcing scheme on the fluid and particles, this article investigates the coherence spectra between the current forcing scheme with both the fluid and the two-way coupling. The results examined are also used to determine whether the influence of the forcing on the fluid and the two-way coupling is statistically significant.

Furthermore, the analysis and results of Abdelsamie and Lee1 are revisited. They examine the temporal evolution of the turbulence kinetic energy in forced turbulence by using the spectral forcing scheme of Eswaran and Pope.6 They proposed that at times smaller than the forcing integral scale the one-point, one-time two-way coupling term is “substantially modified by the forcing scheme,” especially for light particles (Stokes numbers, St < 1). By contrast, it is shown that the new forcing scheme used in this work does not appear in the temporal evolution of the two-way coupling term. This is because the current forcing scheme exactly balances the amount of forcing with the fluid and particle dissipation and also keeps the turbulent kinetic energy constant at a pre-defined value.

This article is organised into four sections. Section II describes the framework of the fluid-phase and the forcing of HIT. Section III presents the simulations and conditions of this study. Section IV compares and assesses the relative performance of the new forcing scheme in terms of the three limitations discussed by Lucci et al.9 Finally, Sec. V summarises the main conclusions of this work.

The momentum equation of the fluid with velocity vf, j is

\begin{eqnarray}\frac{\partial {(v_{f,j}})}{\partial {t}} + \frac{\partial {(v_{f,j} v_{f,i})}}{\partial {x_i}} = - \frac{1}{\rho _f} \frac{\partial {p}}{\partial {x_j}} + \frac{1}{\rho _f}\frac{\partial {(\tau _{ij})}}{\partial {x_i}} + T_j + \Pi _j,\end{eqnarray}
(vf,j)t+(vf,jvf,i)xi=1ρfpxj+1ρf(τij)xi+Tj+Πj,
(1)

where Tj is the term representing the forcing to sustain HIT and the last term on the right hand side is the inter-phase momentum exchange between the fluid and suspended particles, further discussed in Sec. II B. τij is the viscous stress tensor of the fluid, given by

\begin{eqnarray}\tau _{ij} = \mu _f \bigg (\frac{\partial {v_{f,i}}}{\partial {x_j}} + \frac{\partial {v_{f,j}}}{\partial {x_i}}\bigg ),\end{eqnarray}
τij=μfvf,ixj+vf,jxi,
(2)

where μf is the dynamic viscosity of the fluid.

The particle effects are included in the fluid momentum equation (Eq. (1)) as the source term Πj. This source term is the two-way coupling. The PSI-Cell method of Crowe et al.4 is used and the force of the particles on the fluid in each computational cell is volume averaged, yielding

\begin{eqnarray}\Pi _j = -\frac{1}{\rho _f}\frac{1}{V_{cell}}\displaystyle \sum _{p=1}^{N_p}V_p\beta _{(f,p)}\big [v_{f@p,j} - v_{p,j}\big ],\end{eqnarray}
Πj=1ρf1Vcellp=1NpVpβ(f,p)vf@p,jvp,j,
(3)

where the summation is for all particles and fractions of them in each computational cell, Vcell is the local volume of the grid cell which the particles are present, Vp is the particle volume, β(f, p) is the drag function, vf@p, j is the fluid velocity interpolated at the particle, vp, j is the particle velocity, and Np is the total number of particles in the grid cell under consideration. In this way, the Lagrangian particle properties are transformed to Eulerian by volume averaging. Note that this method involves the interpolation of Lagrangian quantities to the fluid cell centres. Yeung and Pope20 perform a study on the interpolation schemes in homogeneous and isotropic turbulence and report that the cubic spline interpolation method has the least effect on the fluid energy spectrum. Therefore, cubic splines are used in this work to determine vf@p, j from the Eulerian velocity field.

The drag function, β, which is the reciprocal of the Eulerian fluid-particle timescale, is given by Wen and Yu19 as

\begin{eqnarray}\beta = \frac{3}{4} C_D \frac{\alpha _p \alpha _f \rho _f \big | v_{f@p,i} - v_{p,i} \big |}{d_p} \alpha ^{-2.65}_f,\end{eqnarray}
β=34CDαpαfρf|vf@p,ivp,i|dpαf2.65,
(4)

where αf represents the fluid volume fraction, αp is the particle volume fraction, dp is the particle diameter, and CD represents the coefficient of drag for an individual particle15 

\begin{eqnarray}C_D =\frac{24\bigg [1 + 0.15 \big ((1 - \alpha _p) Re_p\big )^{0.687}\bigg ]}{Re_p(1 - \alpha _p)},\end{eqnarray}
CD=241+0.15(1αp)Rep0.687Rep(1αp),
(5)

where Rep is the particle Reynolds number. Note that Eq. (5) takes into account the fact that the drag function of a single particle is influenced by the suspension of a large number of particles, hence it includes the particle volume fraction.

The forcing term Tj specified in this work as

\begin{eqnarray}T_j = \frac{\rho _f}{\Delta t} \frac{\sqrt{q^2_{f,wanted}} - \sqrt{q^2_{f,computed}}}{\sqrt{q^2_{f,wanted}}} v_{f,j}^{triggered},\end{eqnarray}
Tj=ρfΔtqf,wanted2qf,computed2qf,wanted2vf,jtriggered,
(6)

where Δt is the fluid time-step,

$q^2_{f,wanted}$
qf,wanted2 is the specified turbulence kinetic energy,
$q^2_{f,computed}$
qf,computed2
is the computed turbulence kinetic energy of the domain, and
$v_{f,j}^{triggered}$
vf,jtriggered
is a carefully chosen pseudo-velocity field, which is synthesised from a model spectrum. Specifically,
$v_{f,j}^{triggered}$
vf,jtriggered
is synthesised from the von-Karman spectrum sampled from pseudo-wavenumber vectors which inherently satisfy continuity. It is important to note that
$v_{f,j}^{triggered}$
vf,jtriggered
is not related to the velocity field obtained by solving the fluid momentum equation. Note that this forcing scheme has the added benefit of enabling the triggering of specific wavenumber modes.

The random fluctuations,

$v_{f,j}^{triggered}$
vf,jtriggered⁠, can be represented in discrete Fourier representation12 

\begin{equation}v_{f,j}^{triggered}= 2\sum ^{N}_{n=1} (\hat{u}^n \cos {\kappa ^n_ix_i + \psi ^n)}\sigma ^n_j,\end{equation}
vf,jtriggered=2n=1N(ûncosκinxi+ψn)σjn,
(7)

where n is the mode number (note that n are chosen to trigger specific wavenumber modes), N is the total number of modes,

$\hat{u}^n$
ûn is the amplitude of the fluctuations at each mode, which is determined from the model spectrum (i.e.,
$\hat{u} = \sqrt{E(\kappa ^n)\Delta \kappa }$
û=E(κn)Δκ
), ψn is the phase of each wavenumber mode, which is obtained randomly (the range is 0 ⩽ ψn ⩽ 2π), and
$\sigma ^n_j$
σjn
is the direction of each mode,
$\kappa ^n_i$
κin
are random wavenumber vectors in the three orthogonal directions. The aim of this scheme is to sustain the turbulence kinetic energy by replenishing the energy lost from any dissipation mechanism by only forcing a specified set of wavenumbers. A time correlation similar to the fluid integral time scale is used for the forcing τFTE, int in order to obtain smooth forcing, where TE, int is the Eulerian integral timescale. In addition, the forcing is performed only over a small range of wavenumber modes to avoid large fluctuations of the turbulence kinetic energy so that most wavenumbers are not affected by the forcing. This forcing scheme, therefore, overcomes the three issues raised by Lucci et al.9 discussed in Sec. I and provides an alternative forcing scheme in physical space.

Simulations on a 1283 periodic box (Lbox = 0.128 m) were performed with fluid kinematic viscosity νf = 1.47 × 10−5m2 s−1 and fluid density ρf = 1.17 kgm−3 with particle volume fraction of 1.37 × 10−5, corresponding to 223, 520 particles. The particles considered in this work are spherical and elastic with a fixed diameter of dp = 67.6 μm. This diameter ensures that the particles are much smaller than the Kolmogorov microscale, which is dpκ = 0.1. The density of the particles varied as follows: 150 kgm−3, 2500 kgm−3, 8000 kgm−3, and 12 000 kgm−3 in order to obtain particles of different Stokes numbers. The corresponding Stokes number (St) of the particles is St = 0.07, St = 1.12, St = 3.45, and St = 5.15. Note that the turbulence kinetic energy levels remained constant at about 6.52 × 10−3 m2 s−2. The Taylor Reynolds number, Reλ is 35.4. Post-processing statistics were gathered after 10 integral time scales and sampled for a time of 20 integral time scales. Note that κ is defined as 2πn/Lbox, where n is the wavenumber, and the wavelength spacing is 49.1m−1.

The influence of the forcing scheme on the nonlinear transfer function is examined by considering the single-phase two-point spatial velocity correlation arising from Eq. (1) for a homogeneous and isotropic flow

\begin{eqnarray}\frac{\partial \langle v^{\prime \prime }_i v^{\prime }_j \rangle }{\partial t} + \frac{\partial \Big (\langle v^{\prime }_j v^{\prime }_k v^{\prime \prime }_i \rangle - \langle v^{\prime \prime }_i v^{\prime \prime }_k v^{\prime }_j \rangle \Big )}{\partial r_k} =&& -\frac{1}{\rho }\Bigg [ -\frac{\partial \langle v^{\prime }_j P^{\prime \prime }\rangle }{\partial r_i} + \frac{\partial \langle v^{\prime \prime }_i P^{\prime }\rangle }{\partial r_j} \Bigg ] \nonumber \\&&+ 2 \nu _f \frac{\partial ^2\langle v^{\prime \prime }_i v^{\prime }_j\rangle }{\partial r^2_k} + \langle v^{\prime }_j T^{\prime \prime }_i\rangle + \langle v^{\prime \prime }_i T^{\prime }_j\rangle ,\end{eqnarray}
vivjt+vjvkvivivkvjrk=1ρvjPri+viPrj+2νf2vivjrk2+vjTi+viTj,
(8)

where

$r_i = x^\prime _i - x^{\prime \prime }_i$
ri=xixi and the superscripts and denote the fluctuating components at two different positions and ⟨.⟩ denotes spatial averaging. Equation (8) can be spatially Fourier transformed in order to determine the effects of the forcing scheme Ti on each component in wavenumber space. The Fourier space equations are

\begin{eqnarray}\frac{\partial \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_j \rangle }{\partial t} = && i \kappa _k \Big (\langle \widehat{v}^{\star \prime }_j \widehat{v}^{\star \prime }_k \widehat{v}^{\prime \prime }_i \rangle - \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\prime \prime }_k \widehat{v}^{\star \prime }_j \rangle \Big ) -\frac{1}{\rho }\Big [ i \kappa _i \langle \widehat{v}^{\star \prime }_j \widehat{P}^{\prime \prime }\rangle -i \kappa _j \langle \widehat{v}^{\prime \prime }_i \widehat{P}^{\star \prime }\rangle \Big ]\nonumber \\&&+ 2 \nu _f \kappa ^2_k \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_j\rangle + \langle \widehat{v}^{\star \prime }_j \widehat{T}^{\prime \prime }_i\rangle + \langle \widehat{v}^{\prime \prime }_i \widehat{T}^{\star \prime }_j\rangle ,\end{eqnarray}
v̂iv̂jt=iκkv̂jv̂kv̂iv̂iv̂kv̂j1ρiκiv̂jP̂iκjv̂iP̂+2νfκk2v̂iv̂j+v̂jT̂i+v̂iT̂j,
(9)

where

$\widehat{v}_i$
v̂i denotes the Fourier transform of vi and
$\widehat{v}^{\star }_i$
v̂i
is the complex conjugate of
$\widehat{v}_i$
v̂i
.

1. Influence of linear forcing scheme of Lundgren on the nonlinear transfer function

The forcing scheme of Lundgren11 is examined first to provide a historical perspective. He proposed a forcing scheme, which is in physical space, so the forcing term in Fourier space is

\begin{eqnarray}\widehat{T}_i = A \widehat{v}_i,\end{eqnarray}
T̂i=Av̂i,
(10)

where term A is a ratio of the production to dissipation, thus keeping the total turbulence kinetic energy constant. Substituting Eq. (10) into Eq. (9), the forcing scheme of Lundgren11 can be incorporated into the two-point spatial velocity correlation equations as follows:

\begin{eqnarray}\frac{\partial \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_j \rangle }{\partial t} = && i \kappa _k \Big (\langle \widehat{v}^{\star \prime }_j \widehat{v}^{\star \prime }_k \widehat{v}^{\prime \prime }_i \rangle - \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\prime \prime }_k \widehat{v}^{\star \prime }_j \rangle \Big ) -\frac{1}{\rho }\Big [ i \kappa _i \langle \widehat{v}^{\star \prime }_j \widehat{P}^{\prime \prime }\rangle -i \kappa _j \langle \widehat{v}^{\prime \prime }_i \widehat{P}^{\star \prime }\rangle \Big ] \nonumber \\&& - 2 \nu _f \kappa ^2_k \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_j\rangle + \langle \widehat{v}^{\star \prime }_j A\widehat{v}^{\prime \prime }_i\rangle + \langle \widehat{v}^{\prime \prime }_i A\widehat{v}^{\star \prime }_j\rangle .\end{eqnarray}
v̂iv̂jt=iκkv̂jv̂kv̂iv̂iv̂kv̂j1ρiκiv̂jP̂iκjv̂iP̂2νfκk2v̂iv̂j+v̂jAv̂i+v̂iAv̂j.
(11)

Setting i = j, it is immediately apparent that the forcing term,

$2A\langle \widehat{v}^{\star \prime }_i \widehat{v}^{\prime \prime }_i\rangle$
2Av̂iv̂i has the same form as the dissipation term,
$2 \nu _f \kappa ^2_k \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_i\rangle$
2νfκk2v̂iv̂i
, hence Eq. (11) simplifies to

\begin{eqnarray}\frac{\partial \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_i \rangle }{\partial t} = i \kappa _k \Big (\langle \widehat{v}^{\star \prime }_i \widehat{v}^{\star \prime }_k \widehat{v}^{\prime \prime }_i \rangle - \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\prime \prime }_k \widehat{v}^{\star \prime }_i \rangle \Big ) + 2(A - \nu _f \kappa ^2_k) \langle \widehat{v}^{\star \prime }_i \widehat{v}^{\prime \prime }_i\rangle ,\end{eqnarray}
v̂iv̂it=iκkv̂iv̂kv̂iv̂iv̂kv̂i+2(Aνfκk2)v̂iv̂i,
(12)

where the pressure term contribution drops out.8 Integrating over spherical shells of radius κi = |κi| yields the dynamical equation for the energy spectrum,E(κ), as

\begin{eqnarray}\frac{\partial E(\kappa ,t)}{\partial t} = T(\kappa , t) + 2(A -\nu _f \kappa ^2)E(\kappa ,t).\end{eqnarray}
E(κ,t)t=T(κ,t)+2(Aνfκ2)E(κ,t).
(13)

In stationary HIT, the term ∂E(κ, t)/∂t = 0, hence the transfer function, T(κ), is given by

\begin{eqnarray}T(\kappa , t) = 2(\nu _f\kappa ^2 - A)E(\kappa ,t).\end{eqnarray}
T(κ,t)=2(νfκ2A)E(κ,t).
(14)

Equation (14) is also given by Rosales and Meneveau.14 It is obvious that this forcing scheme affects all wavenumbers of the transfer function.

Figure 1 illustrates a schematic of the individual terms of the particle-laden dynamic energy equation in forced turbulence using the forcing scheme of Lundgren.11 As already shown previously via Eq. (14), this forcing scheme affects all wavenumbers. Thus any effects due to the particles on the transfer function and energy spectrum of the fluid cannot be separated from the forcing. This is exactly the point of Elghobashi and Truesdell5 and Lucci et al.9 

FIG. 1.

Schematic of transfer function (T(κ), continuous line), dissipation spectrum (νfκ2E(κ) sparsely dashed line), two-way coupling spectrum (Ψ(κ) densely dashed line), and the forcing scheme of Lundgren11 (AE(κ) dotted-dashed line). Note that the figure is an illustration and it is not obtained from DNS.

FIG. 1.

Schematic of transfer function (T(κ), continuous line), dissipation spectrum (νfκ2E(κ) sparsely dashed line), two-way coupling spectrum (Ψ(κ) densely dashed line), and the forcing scheme of Lundgren11 (AE(κ) dotted-dashed line). Note that the figure is an illustration and it is not obtained from DNS.

Close modal

2. Influence of new scheme on the nonlinear transfer function

Performing a similar procedure for the new forcing scheme, expressed by Eq. (6), the evolution of

$\langle \widehat{v^{\prime \prime }_i} \widehat{v^{\star \prime }_i} \rangle$
vîvî with time is

\begin{eqnarray}\frac{\partial \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_i \rangle }{\partial t} = && i \kappa _k \Big (\langle \widehat{v}^{\star \prime }_i \widehat{v}^{\star \prime }_k \widehat{v}^{\prime \prime }_i \rangle - \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\prime \prime }_k \widehat{v}^{\star \prime }_i \rangle \Big ) - 2 \nu _f \kappa ^2_k \langle \widehat{v}^{\star \prime }_i \widehat{v}^{\prime \prime }_i \rangle \nonumber \\&& + B(\kappa )\Big (\langle \widehat{v}^{\star \prime }_i \widehat{v}^{\prime \prime triggered}_i\rangle + \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime triggered}_i\rangle \Big ),\end{eqnarray}
v̂iv̂it=iκkv̂iv̂kv̂iv̂iv̂kv̂i2νfκk2v̂iv̂i+B(κ)v̂iv̂itriggered+v̂iv̂itriggered,
(15)

where B(κ, t) acts as a top hat function, where it is given as

\begin{eqnarray}B(\kappa ) = \left\lbrace {\begin{array}{c@{\quad}c}\rho _f \frac{\sqrt{q^2_{f,wanted}} - \sqrt{q^2_{f,computed}}}{\sqrt{q^2_{f,wanted}}}\Delta t&, \mbox{if } \kappa _{min} \le \kappa \le \kappa _{max}\\\\0 &, \mbox{if } \kappa < \kappa _{min} \mbox{ and } \kappa > \kappa _{max}. \end{array}}\right.\end{eqnarray}
B(κ)=ρfqf,wanted2qf,computed2qf,wanted2Δt,ifκminκκmax0,ifκ<κminandκ>κmax.
(16)

Integrating over spherical shells yields

\begin{eqnarray}\frac{\partial E(\kappa ,t)}{\partial t} = T(\kappa , t) - 2\nu _f \kappa ^2E(\kappa ,t) + B(\kappa ,t)F(\kappa ,t),\end{eqnarray}
E(κ,t)t=T(κ,t)2νfκ2E(κ,t)+B(κ,t)F(κ,t),
(17)

where F(κ, t) is the three-dimensional Fourier transform of the current forcing scheme. In stationary HIT, the transfer function is now expressed as

\begin{eqnarray}T(\kappa , t) = 2\nu _f\kappa ^2E(\kappa ,t) - B(\kappa ,t)F(\kappa ,t).\end{eqnarray}
T(κ,t)=2νfκ2E(κ,t)B(κ,t)F(κ,t).
(18)

It is obvious now that the forcing scheme affects the transfer function only in the range where the forcing scheme is active, i.e., κmin ⩽ κ ⩽ κmax. For the simulations described later, forcing is applied only in the range 3κo ⩽ κ ⩽ 7κo, where κo = Lbox/(2π).

When particles are added, two extra terms appear in Eq. (15), i.e.,

\begin{eqnarray}\frac{\partial \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime }_i \rangle }{\partial t} = && i \kappa _k \Big (\langle \widehat{v}^{\star \prime }_i \widehat{v}^{\star \prime }_k \widehat{v}^{\prime \prime }_i \rangle - \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\prime \prime }_k \widehat{v}^{\star \prime }_i \rangle \Big ) - 2 \nu _f \kappa ^2_k) \langle \widehat{v}^{\star \prime }_i \widehat{v}^{\prime \prime }_i \rangle \nonumber \\&& + B(\kappa )\Big (\langle \widehat{v}^{\star \prime }_i \widehat{v}^{\prime \prime triggered}_i \rangle + \langle \widehat{v}^{\prime \prime }_i \widehat{v}^{\star \prime triggered}_i\rangle \Big ) + \Big (\langle \widehat{v}^{\star \prime }_i \widehat{\Pi }^{\prime \prime }_i\rangle + \langle \widehat{v}^{\prime \prime }_i \widehat{\Pi }^{\star \prime }_i\rangle \Big ).\end{eqnarray}
v̂iv̂it=iκkv̂iv̂kv̂iv̂iv̂kv̂i2νfκk2)v̂iv̂i+B(κ)v̂iv̂itriggered+v̂iv̂itriggered+v̂iΠ̂i+v̂iΠ̂i.
(19)

Hence in stationary HIT, the transfer function including the effect of particles is

\begin{eqnarray}T(\kappa , t) = 2\nu _f\kappa ^2E(\kappa ,t) - B(\kappa ,t)F(\kappa ,t) - \Psi (\kappa ,t),\end{eqnarray}
T(κ,t)=2νfκ2E(κ,t)B(κ,t)F(κ,t)Ψ(κ,t),
(20)

where Ψ(κ, t) is the three-dimensional two-way coupling spectrum. From Eq. (20) it is clear that when B(κ, t) is zero, the transfer function, T(κ, t) can only be directly affected by the particles and the fluid viscosity, but not by the forcing. On the other hand, in the range when B(κ, t) is non-zero, however, the forcing can directly affect the transfer function of the fluid.

In the simulations described below the forcing is confined to the range of wavenumbers κmin ⩽ κ ⩽ κmax. A schematic for this forcing scheme is illustrated in Figure 2, which shows that the current forcing scheme is only present at the forced wavenumbers. Note, however, that when particles are added in the domain the current forcing scheme adjusts itself to a new value in order to keep the integral

$\int ^\infty _0 T(\kappa )d\kappa = 0$
0T(κ)dκ=0⁠.

FIG. 2.

Schematic of transfer function (T(κ), continuous line), dissipation spectrum (νfκ2E(κ) sparsely dashed line), two-way coupling spectrum (Ψ(κ) densely dashed line) and the new forcing scheme (B(κ)F(κ) dotted-dashed line), which are described by Eq. (20). Note that the figure is an illustration and it is not obtained from DNS.

FIG. 2.

Schematic of transfer function (T(κ), continuous line), dissipation spectrum (νfκ2E(κ) sparsely dashed line), two-way coupling spectrum (Ψ(κ) densely dashed line) and the new forcing scheme (B(κ)F(κ) dotted-dashed line), which are described by Eq. (20). Note that the figure is an illustration and it is not obtained from DNS.

Close modal

To obtain a clear picture on the effects of the forcing on the transfer function a statistical analysis has been performed in order to determine the statistical coherence between the fluid and the forcing. This is further described in Sec. IV B below where the spatial coherence spectra between the forcing and the particles and between the forcing and fluid is studied.

This section investigates both the temporal and spatial coherence spectra between the two-way coupling and the fluid velocity, between the forcing and the fluid velocity and between the forcing and the two-way coupling. The mean squared coherence, as given in Lumley10 [Chap. 4, p. 116] of the cross-spectrum between a signal A and signal B is CAB(κ) is defined as

\begin{eqnarray}C_{AB}(\kappa ) = \frac{\Big | G_{AB}(\kappa )\Big |^2}{G_{AA}(\kappa )G_{BB}(\kappa )}.\end{eqnarray}
CAB(κ)=|GAB(κ)|2GAA(κ)GBB(κ).
(21)

GAA(κ) and GBB(κ) are the spectra of signals A and B, respectively. The mean squared coherence is similar to the cross-correlation coefficients, with the only difference that a value of the correlation coefficient at each κ is obtained. In addition, the phase θAB(κ) of the cross-spectrum GAB(κ) is defined as

\begin{eqnarray}\theta _{AB}(\kappa ) = \frac{\Im \Big \lbrace G_{AB}(\kappa )\Big \rbrace }{\Re \Big \lbrace G_{AB}(\kappa )\Big \rbrace },\end{eqnarray}
θAB(κ)=GAB(κ)GAB(κ),
(22)

where ℑ{.} and ℜ{.} are the imaginary and real parts of the cross-spectrum GAB(κ), respectively. Note that θ(κ) varies from

$-180^\text{o}$
180o to
$+180^\text{o}$
+180o
, θ(κ) = 0 means that signals A and B are in-phase. A randomly distributed phase across κ means that signals A and B are incoherent and any value on the coherence spectrum is unimportant.

1. Spatial coherence and phase

This section investigates three spatial coherence spectra and the corresponding phase of the cross-spectra for forced turbulence with particles with St = 0.07 (light particles) and particles with St = 3.45 (heavy particles) at Reλ = 35.4. These are cross-spectra between: (a) the two-way coupling and the fluid velocity, (b) the forcing and the fluid velocity, and (c) the forcing and the two-way coupling. Note that the investigated spectra are one-dimensional in direction 1. Figures 3(a) and 4(a) compare the aforementioned spatial coherence spectra for light and heavy particles, respectively. The forcing and particles for both Stokes numbers are moderately correlated, C(κ) < 0.17, in the range where forcing is acting, 100 ⩽ κ ⩽ 300. Outside this range, the correlation quickly drops to very low values (less than 0.001) which confirms that the forcing does not influence the two-way coupling spectrum at these wavenumbers. Similarly, the forcing does not affect the fluid scales for κ > 300, because the correlation between the forcing and the fluid is C(κ) ≪ 0.01.

FIG. 3.

Comparison of (a) spatial coherence spectrum and (b) spatial phase between particles and fluid, between forcing and fluid and between forcing and particles for St = 0.07 and Reλ = 35.4.

FIG. 3.

Comparison of (a) spatial coherence spectrum and (b) spatial phase between particles and fluid, between forcing and fluid and between forcing and particles for St = 0.07 and Reλ = 35.4.

Close modal
FIG. 4.

Comparison of (a) spatial coherence spectrum and (b) spatial phase between particles and fluid, between forcing and fluid and between forcing and particles for St = 3.45 and Reλ = 35.4.

FIG. 4.

Comparison of (a) spatial coherence spectrum and (b) spatial phase between particles and fluid, between forcing and fluid and between forcing and particles for St = 3.45 and Reλ = 35.4.

Close modal

Figures 3(b) and 4(b) compare the phase of the aforementioned cross-spectra. Outside the forcing range, κ > 300, the phase fluctuates randomly for both fluid-forcing and particle-forcing cross-spectra. This means that the forcing at these wavenumbers is completely uncorrelated to both fluid and particles because the phase is incoherent. On the other hand, the phase of the cross-spectrum between the two-way coupling and the fluid are in-phase and then in anti-phase, i.e., alternating between

$-180^\text{o}$
180o and
$+180^\text{o}$
+180o
. In addition, the two-way coupling and fluid coherence spectrum for both particles appear to have double peaks. The second peak is after the signals are out of phase. Moreover, the phases shown for each individual cross-spectrum are related as θAB(κ) = −θBA(κ). This means that the combined cross-correlations are symmetric; which means in Fourier space, the summation of the last four terms in Eq. (15) is only real coefficients. As shown in Figure 5, which compares the combined phase of the cross-spectra, the combined phase between the forcing and the fluid and the forcing and the two-way coupling alternates between
$0^\text{o}$
0o
and
$180^\text{o}$
180o
for the unforced wavenumbers. This means that the coefficients of these spectra at these wavenumbers are incoherent. On the other hand, the combined phase between the two-way coupling and the fluid remains in phase at the forced wavenumbers and then switches to +180° and remains at this angle. This is because the cross-spectrum between the two-way coupling and the fluid at these wavenumbers is negative. This also explains the double peak of the corresponding coherence spectra of Figures 3(a) and 4(a). This behaviour of the two-way coupling and fluid spectrum could be a consequence of the particles generating fluid scales at high wavenumbers. This is because at low wavenumbers the phase between the fluid and the two-way coupling is zero, which means the fluid is “forcing” the particles. On the other hand, at high wavenumbers and as reflected in the phase shift, the particles are now “forcing” the fluid. This is further examined is Sec. ??? where particles of various Stokes numbers are examined.

FIG. 5.

Phase of combined cross-spectra between particles and fluid, between forcing and fluid, and between forcing and particles for St = 3.45 and Reλ = 35.4.

FIG. 5.

Phase of combined cross-spectra between particles and fluid, between forcing and fluid, and between forcing and particles for St = 3.45 and Reλ = 35.4.

Close modal

2. Spatial mean square coherence and phase of fluid velocity and two-way coupling

Figure 6 compares the different coherence functions resulting from various Stokes numbers at Reλ = 35.4. The value of C(0) decreases with decreasing Stokes number, which means that for light particles the fluid velocity is almost incoherent with the two-way coupling term. When the Stokes number increases, the fluid velocity and the two-way coupling term become more coherent, but only at low wavelengths (about κ < 700). At κ > 1000 there is a slight increase in the coherence spectrum. This means that: (a) the two-way coupling interacts with the intermediate fluid structures and interacts with a smaller effect at higher wavelengths. The reason the coherence spectrum peaks at κ = 0 is because one-dimensional spectra in the longitudinal directions are used. Tennekes and Lumley18 show that isotropic one-dimensional spectra in the longitudinal direction peak at κ = 0. One-dimensional spectra in both the longitudinal or lateral directions, by contrast to three-dimensional spectra, are non-zero at κ = 0 because they are proportional to an integral scale.8 (b) the amount of interaction decreases with decreasing Stokes number because the particles become less inertial and therefore lose their ability to affect the fluid.

FIG. 6.

Mean squared coherence of cross-spectrum

$G^{(1)}_{v_{f,1} \Pi }$
Gvf,1Π(1) with different Stokes numbers at Reλ = 35.4. Note that the superscript (1) indicates one-dimensional spectrum in direction 1.

FIG. 6.

Mean squared coherence of cross-spectrum

$G^{(1)}_{v_{f,1} \Pi }$
Gvf,1Π(1) with different Stokes numbers at Reλ = 35.4. Note that the superscript (1) indicates one-dimensional spectrum in direction 1.

Close modal

From Figure 7, it is apparent that the phase shift from 0° to +180° is Stokes number dependent. For very light particles the shift occurs at high wavenumbers, whereas for heavy particles (St = 3.45 and St = 5.15) the phase shift occurs at lower wavenumbers. However, this transition in phase shift is not linear with Stokes number. Figure 7 also compares the phase shift behaviour for St ≈ 1 particles. The phase shift in this case occurs at very low wavenumbers (even lower for St = 5.15). This is because, although not presented, St ≈ 1 particles preferentially concentrate, whereas particles of different Stokes numbers (St < 1 and St > 1) appear to be more randomly distributed in the domain. This behaviour is also physical and cannot be attributed to any effects of the artificial forcing.

FIG. 7.

Combined phase of cross-spectrum

$G^{(1)}_{v_{f,1}\Pi }$
Gvf,1Π(1) with different Stokes numbers at Reλ = 35.4. Note that the superscript (1) indicates one-dimensional spectrum in direction 1.

FIG. 7.

Combined phase of cross-spectrum

$G^{(1)}_{v_{f,1}\Pi }$
Gvf,1Π(1) with different Stokes numbers at Reλ = 35.4. Note that the superscript (1) indicates one-dimensional spectrum in direction 1.

Close modal

Figure 8 compares the modulation of the single-phase transfer function when particles of different Stokes numbers (St = 0.07, St = 1.12, and St = 3.45) are added in the domain at Reλ = 35.4. The figure clearly shows that when particles are added the transfer function at the forced wavenumbers is also modulated. But even though the effects of the particles at the forced wavenumbers cannot be clearly quantified in this range, the simulations are still correct at the unforced wavenumbers, only the separate effects cannot be distinguished. The net effect over the forcing range can be determined, however, by using the information obtained over the unforced range. This is simply because the integral

$\int ^\infty _0 T(\kappa )d\kappa = 0$
0T(κ)dκ=0⁠. In addition, since this integral of the transfer function is equal to zero for all the cases shown in Figure 8, it can be deduced that the forcing scheme responds differently to particles with different Stokes numbers. This is shown in Figure 8 where at the forced wavenumbers the peak of the transfer function becomes less negative with increasing Stokes number. Moreover, the modulation of the transfer function in this work shows similar trends to the works of Elghobashi and Truesdell5 and Ferrante and Elghobashi.7 

FIG. 8.

Comparison between single phase and particle-laden; St = 0.07, St = 1.12, and St = 3.45, transfer functions at Reλ = 35.4.

FIG. 8.

Comparison between single phase and particle-laden; St = 0.07, St = 1.12, and St = 3.45, transfer functions at Reλ = 35.4.

Close modal

This section investigates the equation proposed by Abdelsamie and Lee1 who study the temporal evolution of the two-way coupling in stationary HIT. In homogeneous turbulence, the temporal evolution of turbulence kinetic energy is expressed as

\begin{eqnarray}\frac{\partial {\frac{1}{2}\langle v_{f,i} v_{f,i} \rangle }}{\partial {t}} = F(t) - \epsilon _f(t) + \epsilon _p(t),\end{eqnarray}
12vf,ivf,it=F(t)εf(t)+εp(t),
(23)

where F(t) is the forcing term expressed as

$F(t) = \langle v_{f_,i}(t) T_{i}(t)\rangle$
F(t)=vf,i(t)Ti(t) and εp(t) is the total particle dissipation defined as

\begin{eqnarray}\epsilon _p(t) &&= \langle v_{f,i}(t) \Pi _i(t) \rangle \nonumber \\&&= -\beta \langle v_{f,i}(t) v_{f,i}(t) \rangle + \beta \langle v_{f,i}(t) v_{p,i}(t)\rangle ,\end{eqnarray}
εp(t)=vf,i(t)Πi(t)=βvf,i(t)vf,i(t)+βvf,i(t)vp,i(t),
(24)

where β is the drag function, which is defined by Eq. (4) and assumed to be constant for simplicity. Note that all quantities are expressed in the Eulerian framework. Differentiating Eq. (24) the particle dissipation with time, the following expression is obtained:

\begin{eqnarray}\frac{d\epsilon _p(t)}{dt} = - \beta \frac{d\langle v_{f,i}(t) v_{f,i}(t) \rangle }{dt} + \beta \frac{\langle v_{f,i}(t)v_{p,i}(t) \rangle }{dt}.\end{eqnarray}
dεp(t)dt=βdvf,i(t)vf,i(t)dt+βvf,i(t)vp,i(t)dt.
(25)

Substituting Eq. (23) into Eq. (25), a first order differential equation results

\begin{eqnarray}\frac{d \epsilon _p(t)}{dt} + 2\beta \epsilon _p(t) = -2\beta \Bigg ( F(t) - \epsilon _f(t) \Bigg ) + \beta \frac{d \langle v_{f,i}(t)v_{p,i}(t) \rangle }{dt}.\end{eqnarray}
dεp(t)dt+2βεp(t)=2βF(t)εf(t)+βdvf,i(t)vp,i(t)dt.
(26)

Note that in the analysis of Abdelsamie and Lee1 the term

$\beta \frac{d \langle v_{f,i}(t)v_{p,i}(t) \rangle }{dt}$
βdvf,i(t)vp,i(t)dt is neglected and the dissipation term 2βεf(t) is missing. They then continue their analysis and derive an equation which shows that the temporal evolution of the two-way coupling term is directly affected by the forcing.

In stationary HIT, however, Eq. (23) is exactly equal to zero because the forcing balances both the dissipation due to the fluid and due to the particles, i.e., F(t) = εf(t) − εp(t). Hence, Eq. (26) simplifies to

\begin{eqnarray}\frac{d \epsilon _p(t)}{dt} = \beta \frac{d \langle v_{f,i}(t)v_{p,i}(t) \rangle }{dt}.\end{eqnarray}
dεp(t)dt=βdvf,i(t)vp,i(t)dt.
(27)

Integrating Eq. (27) results in

\begin{eqnarray}\epsilon _p(t) = \epsilon _p(t_o) + \beta \Big ( \langle v_{f,i}(t)v_{p,i}(t)\rangle - \langle v_{f,i}(t_o)v_{p,i}(t_o)\rangle \Big ).\end{eqnarray}
εp(t)=εp(to)+βvf,i(t)vp,i(t)vf,i(to)vp,i(to).
(28)

From Eq. (24) and at t = to, then εp(to) = −β⟨vf, i(to)vf, i(to)⟩ + β⟨vf, i(to)vp, i(to)⟩, hence Eq. (28) simplifies to

\begin{eqnarray}\epsilon _p(t) = \beta \Big ( \langle v_{f,i}(t)v_{p,i}(t)\rangle - \langle v_{f,i}(t_o)v_{f,i}(t_o)\rangle \Big ).\end{eqnarray}
εp(t)=βvf,i(t)vp,i(t)vf,i(to)vf,i(to).
(29)

Equation (29) illustrates that the two-way coupling term is independent from the forcing term. In fact, it is dependent on the product between the drag function β and the covariance between the fluid and particle velocities, ⟨vf, i(t)vp, i(t)⟩, which acts as a transient term. The last term on the right hand side of Eq. (29) is the product between the drag function and the variance of the fluid at t = to. With the current forcing scheme there is no effect of the forcing on temporal evolution of the two-way coupling.

The purpose of this work was to test the suitability and evaluate the performance of a newly proposed forcing scheme for homogeneous and isotropic turbulence laden with particles. In view of the problems that previous simulations have faced the three limitations concerning forced particle HIT discussed by Lucci et al.9 were examined. The new scheme was evaluated using both theory and statistics from DNS simulations of forced turbulence simulations with particles. The nonlinear transfer function was examined in order to determine which wavenumber modes were influenced by the forcing scheme. The theoretical analysis showed that the newly proposed scheme only affects T(κ, t) at the wavenumber modes where the forcing is active. Thus the additional effect of particles can be quantified at all other wavenumbers. The effects on E(κ, t) due to the particles at intermediate and high wavenumbers is possible because the transfer function is not influenced directly outside of the range of the forcing. By contrast, by performing a similar analysis with the linear forcing of Lundgren,11 it is shown that the transfer function was affected at all wavenumbers so that the effects of the particles cannot be clearly quantified.

The spatial coherencespectra between the two-way coupling and the fluid velocity, between the forcing and the fluid velocity and between the forcing and the two-way coupling were investigated in order to determine the cross-correlation coefficients at each wavenumber mode. The spatial coherence spectra show that the forcing is moderately correlated with the particles at the wavenumber range it is acting (less than 0.17), but quickly drops to very low correlation values (less than 0.001) outside the forcing range; the latter is effectively zero to within the computational error. Similarly, it was shown that the forcing does not affect the fluid velocity outside its working range. Moreover, the phase of the combined cross-spectra between the forcing and the fluid and between the forcing and the particles was examined. The results show that outside the working range of the forcing, the phase of these spectra fluctuates randomly between

$\pm 180^\text{o}$
±180o⁠. This indicates that at these wavenumbers the forcing is uncorrelated with both the fluid and the particles. On the other hand, the phase of the cross-spectrum between the fluid and the particles changes from
$0^\text{o}$
0o
to
$180^\text{o}$
180o
with increasing wavenumbers. This means that the fluid is directly influencing the particles at low wavenumbers but at large wavenumbers the particles interact with the fluid structures and induce fluid motion.

The temporal evolution of the turbulence kinetic energy and two-way coupling term were investigated in order to determine analytically whether the newly proposed forcing affects the two-way coupling spectrum. The analysis showed that the forcing did not affect the two-way coupling term. The results of this work showed that the newly proposed forcing scheme can be used to examine forced homogeneous and isotropic turbulence with particles. The benefit of this forcing scheme over the forcing scheme of Eswaran and Pope6 is that the scheme is in physical space. Codes that are implemented in physical space do not need to Fourier transform quantities, thereby reducing the computational time of the simulation.

The authors are grateful to the Engineering and Physical Sciences Research Council (EPSRC) for their financial support (Grant No. EP/G049262/1). The simulations have been performed on the MECTOR computer cluster of Imperial College London.

1.
A. H.
Abdelsamie
and
C.
Lee
, “
Decaying versus stationary turbulence in particle-laden isotropic turbulence: Turbulence modulation mechanism
,”
Phys. Fluids
24
(
1
),
015106
(
2012
).
2.
M.
Boivin
,
O.
Simonin
, and
K. D.
Squires
, “
Direct numerical simulation of turbulence modulation by particles in isotropic turbulence
,”
J. Fluid Mech.
375
(
1
),
235
263
(
1998
).
3.
G.
Comte-Bellot
and
S.
Corrsin
, “
Simple Eulerian time correlation of full-and narrow-band velocity signals in grid-generated, isotropic turbulence
,”
J. Fluid Mech.
48
(
02
),
273
337
(
1971
).
4.
C. T.
Crowe
,
M. P.
Sharma
, and
D. E.
Stock
, “
The particle-source-in cell (PSI-CELL) model for gas-droplet flows
,”
J. Fluids Eng.
99
(
2
),
325
333
(
1977
).
5.
S.
Elghobashi
and
G. C.
Truesdell
, “
On the two-way interaction between homogeneous turbulence and dispersed solid particles. I: Turbulence modification
,”
Phys. Fluids A
5
(
7
),
1790
(
1993
).
6.
V.
Eswaran
and
S. B.
Pope
, “
An examination of forcing in direct numerical simulations of turbulence
,”
Comput. Fluids
16
(
3
),
257
278
(
1988
).
7.
A.
Ferrante
and
S.
Elghobashi
, “
On the physical mechanisms of two-way coupling in particle-laden isotropic turbulence
,”
Phys. Fluids
15
(
2
),
315
(
2003
).
8.
W. K.
George
, “
Lectures in turbulence for the 21st century
,” Lecture Notes (
2013
), http://turbulence-online.com/Publications/Lecture_Notes/Turbulence_Lille/TB_16January2013.pdf.
9.
F.
Lucci
,
A.
Ferrante
, and
S.
Elghobashi
, “
Modulation of isotropic turbulence by particles of Taylor length-scale size
,”
J. Fluid Mech.
650
(
2010
),
5
55
(
2010
).
10.
J. L.
Lumley
,
Stochastic Tools in Turbulence
(
Dover Publications Inc.
,
New York
,
1970
).
11.
T. S.
Lundgren
, “
Linearly forced isotropic turbulence
,” Technical Report 2, DTIC Document,
2003
.
12.
S. A.
Orszag
and
G. S.
Patterson
, “
Numerical simulation of three-dimensional homogeneous isotropic turbulence
,”
Phys. Rev. Lett.
28
(
2
),
76
79
(
1972
).
13.
J. J.
Riley
and
G. S.
Patterson
, “
Diffusion experiments with numerically integrated isotropic turbulence
,”
Phys. Fluids
17
(
2
),
292
(
1974
).
14.
C.
Rosales
and
C.
Meneveau
, “
Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties
,”
Phys. Fluids
17
,
095106
(
2005
).
15.
P. N.
Rowe
, “
Drag forces in a hydraulic model of a fluidized bed, part II
,”
Trans. Inst. Chem. Eng.
39
,
175
180
(
1961
).
16.
K. D.
Squires
and
J. K.
Eaton
, “
Particle response and turbulence modification in isotropic turbulence
,”
Phys. Fluids A
2
,
1191
(
1990
).
17.
K. D.
Squires
and
J. K.
Eaton
, “
Preferential concentration of particles by turbulence
,”
Phys. Fluids A
3
(
5
),
1169
(
1991
).
18.
H.
Tennekes
and
J. L.
Lumley
,
A First Course in Turbulence
, 3rd ed. (
The MIT Press
,
1972
).
19.
C. Y.
Wen
and
Y. H.
Yu
, “
A generalized method for predicting the minimum fluidization velocity
,”
AIChE J.
12
(
3
),
610
612
(
1966
).
20.
P. K.
Yeung
and
S. B.
Pope
, “
An algorithm for tracking fluid particles in numerical simulations of homogeneous turbulence
,”
J. Comput. Phys.
79
,
373
416
(
1988
).