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 [“

## I. INTRODUCTION

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 Corrsin^{3} 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 Patterson^{13} 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.

Lundgren^{11} 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 Pope^{6} 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 Lundgren^{11} 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}κ^{2} − *A*)*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 128^{3} periodic cubic box (of length *L*_{b} = 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 Lee^{1} 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.

## II. MATHEMATICAL DESCRIPTION

### A. Fluid equation

The momentum equation of the fluid with velocity *v*_{f, j} is

where *T*_{j} 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

where μ_{f} is the dynamic viscosity of the fluid.

### B. Formulation of two-way coupling

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

where the summation is for all particles and fractions of them in each computational cell, *V*_{cell} is the local volume of the grid cell which the particles are present, *V*_{p} is the particle volume, β_{(f, p)} is the drag function, *v*_{f@p, j} is the fluid velocity interpolated at the particle, *v*_{p, j} is the particle velocity, and *N*_{p} 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 Pope^{20} 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 *v*_{f@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 Yu^{19} as

where α_{f} represents the fluid volume fraction, α_{p} is the particle volume fraction, *d*_{p} is the particle diameter, and *C*_{D} represents the coefficient of drag for an individual particle^{15}

where *Re*_{p} 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.

### C. New methodology to sustain homogeneous and isotropic turbulence

The forcing term *T*_{j} specified in this work as

where Δ*t* is the fluid time-step,

The random fluctuations,

^{12}

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

^{n}is the phase of each wavenumber mode, which is obtained randomly (the range is 0 ⩽ ψ

^{n}⩽ 2π), and

_{F}≈

*T*

_{E, int}in order to obtain smooth forcing, where

*T*

_{E, 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.

## III. SIMULATIONS AND CONDITIONS

Simulations on a 128^{3} periodic box (*L*_{box} = 0.128 m) were performed with fluid kinematic viscosity ν_{f} = 1.47 × 10^{−5}m^{2} 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 *d*_{p} = 67.6 μm. This diameter ensures that the particles are much smaller than the Kolmogorov microscale, which is *d*_{p}/η_{κ} = 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} m^{2} 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*/*L*_{box}, where n is the wavenumber, and the wavelength spacing is 49.1m^{−1}.

## IV. RESULTS

### A. Influence of forcing on the nonlinear transfer function

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

where

^{′}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

*T*

_{i}on each component in wavenumber space. The Fourier space equations are

where

*v*

_{i}and

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

The forcing scheme of Lundgren^{11} 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

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 Lundgren^{11} can be incorporated into the two-point spatial velocity correlation equations as follows:

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

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

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

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 Truesdell^{5} and Lucci *et al.*^{9}

#### 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

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

Integrating over spherical shells yields

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

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} = *L*_{box}/(2π).

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

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

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

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.

### B. Coherencespectra and phase

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 Lumley^{10} [Chap. 4, p. 116] of the cross-spectrum between a signal A and signal B is *C*_{AB}(κ) is defined as

*G*_{AA}(κ) and *G*_{BB}(κ) 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 *G*_{AB}(κ) is defined as

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

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

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

_{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

#### 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 Lumley^{18} 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.

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.

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

^{5}and Ferrante and Elghobashi.

^{7}

### C. Temporal evolution of turbulence kinetic energy

This section investigates the equation proposed by Abdelsamie and Lee^{1} 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

where *F*(*t*) is the forcing term expressed as

_{p}(

*t*) is the total particle dissipation defined as

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:

Note that in the analysis of Abdelsamie and Lee^{1} the term

_{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

Integrating Eq. (27) results in

From Eq. (24) and at *t* = *t*_{o}, then ε_{p}(*t*_{o}) = −β⟨*v*_{f, i}(*t*_{o})*v*_{f, i}(*t*_{o})⟩ + β⟨*v*_{f, i}(*t*_{o})*v*_{p, i}(*t*_{o})⟩, hence Eq. (28) simplifies to

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, ⟨*v*_{f, i}(*t*)*v*_{p, 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* = *t*_{o}. With the current forcing scheme there is no effect of the forcing on temporal evolution of the two-way coupling.

## V. SUMMARY AND CONCLUSIONS

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

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 Pope^{6} 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.

## ACKNOWLEDGMENTS

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.