text:
# Effect of particle inertia on the turbulence in a suspension

## Introduction

The interaction of solid particles or liquid droplets with the turbulence in a gas controls the performance of various engineering devices and is important for many practical applications like the combustion of pulverized coal and liquid sprays and cyclone separation. This interaction plays also an important role in many areas of environmental science and physics of the atmosphere. Dust storms, rain triggering, dusting and spraying for agricultural or forestry purposes, preparation and processing of aerosols are typical examples. For a review of turbulent flows with particles and droplets see, e.g. the book by C.T. Crowe, M. Sommerfeld and Y.Tsuji B2 .

In dilute suspensions with small volume fractions of particles, $`C_\text{p}`$, the particle-particle interactions are negligible. Nevertheless, for $`\rho _\text{p}/\rho _\text{f}1`$ (the ratio of the solid particle material and the gas densities), the mass loading $`\varphi =C_\text{p}\rho _\text{p}/\rho _\text{f}`$ may exceed unity and the kinetic energies of the particles and the carrier gas may be comarable. Hence the “two-way coupling” effect of the fluid on the particles and vice versa must be accounted for. Current understanding of the turbulence in dilute suspensions is still at its infancy due to the highly nonlinear nature of the physically relevant interactions and a wide spectrum of governing parameters (the particle size $`a`$ vs. $`L`$ and $`\eta `$, the outer and inner scales of turbulence, the particle response time $`\tau _\text{p}`$ vs. $`\gamma __L`$ and $`\gamma _\eta `$, the turnover frequencies of $`L`$\- and $`\eta `$\- scale eddies).

Existing analytical studies of the problem are mainly based upon a two-fluid model description wherein both the carrying fluid and particle phases are treated as interpenetrating continua B2 ; 83EA ; 99DE ; 98BSS . This model deals with non-interacting solid spherical particles with a radius $`a`$ small enough such that:

1. One can neglect the effect of preferential concentration and may assume homogeneity of the particle space distribution. This is not always so. Above some critical radius $`a_{\text{cr}}`$ the space homogeneous distribution of particles becomes unstable. Resulting clustering instability leads to preferential concentration. For a detailed theory of this effect, see Ref. BeerSheva  and references therein. In the present paper we consider only particles with $`a<a_{\text{cr}}`$.
2. The Stokes viscous drag law for particle acceleration, $`d𝒖_\text{p}/dt=\left[𝒖_\text{f}𝒖_\text{p}\right]/\tau _\text{p}`$, is valid ($`𝒖_\text{f}`$ is the fluid velocity).

Unfortunately, the statistical description of two-fluid turbulence with closure procedures requires a set of additional questionable simplifications due to the lack of understanding of the relevant physics of the particle-fluid interactions. This makes closures of the two-fluid model highly qualitative at best 94Egl ; 98BSS ; 00BSS .

We think that the basic physics of the problem may be better described by a simpler one-fluid model for turbulent dilute suspensions, which uses standard closure relations of one-phase turbulence. The present paper suggests such a model and, as a first step, uses a properly modified simple closure, based on the Kolmogorov-Richardson cascade picture of turbulence. The resulting non-linear differential equation for the energy budget were solved analytically. This provides an economical and internally consistent analytical description of the turbulence modification by particles including the dependence of suppression or enhancement of the turbulence on the three governing parameters: $`\left(\tau _\text{p}\gamma __L\right)`$, $`\varphi `$ and the scale of eddies. These effects were previously observed in numerous experimental and numerical publications, see, e.g.the review by Crowe, Trout and Chung B3 . Many groups carried out experimental work; for an overview see Pietryga 01Pie . Other researchers studied the modification of turbulence by small particles using direct numerical simulations 90SE ; 93ET ; 98BSS ; 99SC ; 01Dru  or by large-eddy simulation B13 . Nevertheless the complicated physics of turbulently flowing suspensions in the two-way coupling regime still wait for a detailed analytical description.

Our analytical findings in this paper successfully correlate important features of turbulence modification observed in numerical simulations Refs. 98BSS ; 01Dru ; 99SC . We believe that the one-fluid model (together with more advanced closures of one-phase turbulence) offers an insight in basic physics of particle-laden turbulent flows. The next step in this development should include the effect of preferential concentrations, which was studied so far only for a given turbulent flow field of the carrier fluid BeerSheva .

The paper is organized as follows. In Sec. II we review, after a presentation of the notation (Sec. LABEL:ss:notations) and an evaluation of the characteristic time scales (Sec. LABEL:ss:timescales), some publications about DNS-simulations (Sec. LABEL:ss:DNS), about experimental work (Sec. LABEL:ss:lab-exp) and about some analytical models(Sec. LABEL:ss:models). A critical evaluation of the existing analytical models98BSS ; 71BP ; 77AT ; B5 ; B6 ; B7 ; B11  is made.

In Sec. III we suggest a new one-fluid equation of motion LABEL:NS-basic for turbulently flowing suspensions with small particles. This is a modified version of the Navier-Stokes equation with two wave-number-dependent parameters, $`\rho _{\text{eff}}\left(k\right)`$ and $`\gamma _\text{p}\left(k\right)`$:

* The $`k`$-dependent effective density of suspensions $`\rho _{\text{eff}}\left(k\right)`$ describes the different degree of involvement of heavy particles in turbulent fluctuations with different wave-numbers \[referred to below as $`k`$-eddies\]. For $`k`$-eddies with a turnover time $`1/\gamma \left(k\right)`$, which is much smaller than the particle response time $`\tau _\text{p}`$, the particles may be considered at rest and $`\rho _{\text{eff}}\left(k\right)`$ is about the density of the fluid itself, $`\rho _\text{f}`$. For $`k`$-eddies with $`\tau _\text{p}\gamma \left(k\right)1`$ the effect of the particle inertia may be neglected and particles may be considered as fully involved in the motion of eddies. Therefore for small enough $`k`$ the effective density $`\rho _{\text{eff}}\left(k\right)`$ is close to the mean density of the suspension (fluid plus particles), $`\rho _\text{s}=\rho _\text{f}\left(1+\varphi \right)`$. Our Eq. (16) reasonably describes $`\rho _{\text{eff}}\left(k\right)`$ for all values of $`k`$.
* The damping term $`\gamma _\text{p}\left(k\right)`$, given by Eq. (17) describes the fluid-particle viscous friction. The function $`\gamma _\text{p}\left(k\right)`$ saturates at the level $`1/\tau _\text{p}`$ for small scale eddies with $`\tau _\text{p}\gamma \left(k\right)1`$, when the particles may be considered to be almost at rest. In this regime the damping is $`k`$-independent, while the turnover frequency of $`k`$-eddies $`\gamma \left(k\right)`$ grows with $`k`$. Therefore for large $`k`$ $`\gamma _\text{p}\left(k\right)\gamma \left(k\right)`$ and the particle-induced damping of these $`k`$-eddies may be neglected with respect to their energy loss in the cascade process, which is determined by the frequency $`\gamma \left(k\right)`$. In contrast, for small enough $`k`$ \[when $`\tau _\text{p}\gamma \left(k\right)1`$\] the particles are almost completely involved in the motions of $`k`$-eddies and their contribution to $`\gamma _\text{p}\left(k\right)`$ is suppressed by the factor $`\left[\tau _\text{p}\gamma \left(k\right)\right]^21`$ with respect to $`1/\tau _\text{p}`$.

Our one-fluid model for turbulent suspensions LABEL:NS-basic is first postulated in Sec. LABEL:ss:adv. Its physical interpretation is discussed in Sec. LABEL:ss:NS-physics. A detailed derivation of Eq. (III) is given in Secs. LABEL:ss:ass, LABEL:ss:NSo-derivation and LABEL:ss:NSk-derivation. The most difficult problem here is how to account for the nonlinear effect of the interaction of $`k`$-eddies within the one-fluid model of turbulent suspensions. The suggested form of the nonlinear term LABEL:NL-basic is a modification of the standard Navier-Stokes nonlinearity and is based on:

* a rigorous description of eddy interactions in both limiting cases $`\tau _\text{p}\gamma \left(k\right)1`$ and $`\tau _\text{p}\gamma \left(k\right)1`$
* respect of the fundamental symmetries of the problem – Galilean invariance and conservation of energy.

Section IV deals with the budget of the kinetic energy in turbulently flowing suspensions. One has to account not only for the dissipation of energy due to the fluid-particle friction but also for the effect of particles on the energy redistribution in the system due to the eddy interaction. First we derive in LABEL:ss:energy the budget equation LABEL:budget which accounts for the energy pumping due to a stirring force, energy damping due to the kinematic viscosity and fluid-particle friction and also describes the flux of energy over the scales due to the nonlinearity of the problem. Equation LABEL:budget is exact but unfortunately is not closed. As usual it includes a 3rd order velocity correlation functions. As a first step in the analysis of turbulent suspensions in the framework of our one-fluid model Eq. (III) and the budget Eq. (IV), we use in this paper, sect. LABEL:ss:closure, a simple closure procedure based on the Richardson-Kolmogorov cascade picture of turbulence in which the energy flux is accounted for in a differential approximation. Needless to say that there are various closure procedures for the Navier-Stokes turbulence in the literature. They may be straightforwardly applied to our Eq. (III). This important part of the project will be done elsewhere.

The derived energy balance equations are summarized in Sec. LABEL:ss:basic-eq. They have a very simple and transparent analytical form LABEL:bud2LABEL:rel-1, allowing their effective analytical analysis, see Sects. V and VI. In particular in Sec. V.2 we found a simple solution for the case of micro-particles having a very small response time. In Sec. V.3 we found the iterative solution for the case of a suspension with heavy particles in the inertial interval and analyzed its accuracy in Sec. V.4.

In section VI we analytically describe a complicated interplay between two competitive effects: of the turbulence suppression and the turbulence enhancement in the inertial interval of scales, as well as in the viscous subrange. A brief comparison of our finding with DNS results is done in Sec. VI.5.

In the concluding Sec. VII we summarize the results of the paper and present our ideas for further work.

## I Notations and relevant timescales

SSnotationsNomenclature

* $`\rho _\text{f}`$, $`𝒖(t,𝒓)`$ – density and velocity of the fluid
* $`\stackrel{~}{𝒖}(\omega ,𝒓)`$, $`𝒖(t,𝒌)`$, $`\stackrel{~}{𝒖}(\omega ,𝒌)`$ – Fourier transform of $`𝒖(t,𝒓)`$ with respect to $`t`$ \[$`t\omega `$\], with respect to $`𝒓`$ \[$`𝒓𝒌`$\], and to both variables \[$`t\omega ,r𝒌`$\]
* $`F(t,𝒌)`$, $`\stackrel{~}{F}(\omega ,𝒌)`$ – pair correlation functions of fluid velocity in $`(t,𝒌)`$ and $`(\omega ,𝒌)`$ representation
* $`E\left(k\right)=\rho _\text{f}k^2F(0,k)/2\pi `$ – one dimensional spectrum of the turbulent kinetic energy of the pure fluid (fluid without particles)
* $`\left(k\right)`$ – one dimensional spectrum (of the turbulent kinetic energy) of the suspension
* $`\gamma \left(k\right)`$ – turnover frequency of $`k`$-eddies (turbulent fluctuations of the characteristic scale $`1/k`$). May be understood also as $`1/\tau \left(k\right)`$, where $`\tau \left(k\right)`$ is the life time of $`k`$-eddies. In the Kolmogorov 41 picture of turbulence $`\gamma \left(k\right)k\sqrt{kE\left(k\right)/\rho _\text{f}}`$.
* $`E=\text{d}kE\left(k\right)/2\pi `$, $`=\text{d}k\left(k\right)/2\pi `$ – total turbulent kinetic energy of respectively the pure fluid and the suspension
* $`a`$, $`\rho _\text{p}`$, $`m_\text{p}=4\pi a^3\rho _\text{p}/3`$ – radius, density and mass of the particles
* $`C_\text{p}`$, $`\mathrm{}^3=1/C_\text{p}`$ , volume fraction of particles and volume of suspension per particle
* $`\psi =\left[4\pi a^3/3\right]/\mathrm{}^3`$, $`\varphi =m_\text{p}/\rho _\text{f}\mathrm{}^3`$ – volume fraction and mass loading parameter
* $`\tau _\text{p}`$ – particle response time, also referred to as Stokes time scale
* $`\tau __L`$ – turnover time of the energy containing eddies (of scale $`L`$)
* $`\delta \tau _\text{p}/\tau __L`$ – the particle response time in the units of $`\tau __L`$.
* $`\eta `$; $`v_\eta `$, $`\tau _\eta =\eta /u_\eta `$ – Kolmogorov (viscous) microscale; characteristic velocity and time at scale $`\eta `$ of turbulence
* $`\rho _{\text{eff}}\left(k\right)`$ – effective density of the suspension for turbulent fluctuations of characteristic scale $`1/k`$ \[referred to as $`k`$-eddies\]
* $`\nu `$, $`\nu _{\text{eff}}\left(k\right)`$ – kinematic viscosity of the pure fluid, effective kinematic viscosity of $`k`$-eddies in the suspension
* $`\gamma _\text{p}\left(k\right)`$ – effective damping frequency in the suspension due to the fluid-particle friction
* $`\epsilon \left(k\right)`$ – (one dimensional) flux of the turbulent kinetic energy of the suspension via a sphere of radius $`k`$ in $`𝒌`$-space, also referred to as energy flux over scales.

SStimescalesEvaluation of time scales

The radius of the particles is supposed to be small enough, so that the particle Reynolds number $`e_\text{p}`$ is less than a critical value ($`e_{\text{cr}}`$). In this case we can apply the Stokes approximation (according to which the fluid-particle friction force is proportional to the difference between the particle velocity and the fluid velocity). Careful analysis by Lumley B10  shows that in a turbulent flow the condition for the validity of $`e_\text{p}e_{\text{cr}}`$ may be expressed via the particle radius $`a`$ and the Kolmogorov micro-scale $`\eta `$ in the following way

$`a2\eta \left(\rho _\text{f}/\rho _\text{p}\right)^{1/3}.`$ (1)
It is clear that one of the important parameters in the physics of turbulently flowing suspensions is the ratio of the inertial time scale of the particles (the Stokes time scale $`\tau _\text{p}`$) and the life time $`\tau _\eta `$ of eddies of the Kolmogorov micro-scale. The particle response time is given by

$`\tau _\text{p}`$ $`=`$ $`{\displaystyle \frac{m_\text{p}}{6\pi \nu \rho _\text{f}a}}={\displaystyle \frac{2\rho _\text{p}a^2}{9\rho _\text{f}\nu }},`$ (2)
where we use the expression for the particle mass $`m_\text{p}`$:

$`m_\text{p}={\displaystyle \frac{4\pi }{3}}a^3\rho _\text{p}.`$ (3)
As is well-known the Kolmogorov micro-scale $`\eta `$ is found from the condition that the Reynolds number for eddies of scale $`\eta `$ is equal to unity:

$`e_\eta =\eta v_\eta /\nu =1.`$ (4)
Here $`v_\eta `$ is the characteristic velocity of $`\eta `$-scale eddies. It is related to the turnover time of these eddies in the following manner $`\tau _\eta =\eta /v_\eta `$. This allows us to rewrite the requirement LABEL:ReK as follows

$`\tau _\eta =\eta ^2/\nu .`$ (5)
The ratio of the time-scales $`\tau _\text{p}`$ and $`\tau _\eta `$ immediately follows from Eqs. LABEL:Stime and LABEL:Ktime:

$`{\displaystyle \frac{\tau _\text{p}}{\tau _\eta }}={\displaystyle \frac{2}{9}}{\displaystyle \frac{\rho _\text{p}}{\rho _\text{f}}}{\displaystyle \frac{a^2}{\eta ^2}}.`$ (6)
Substituting the condition LABEL:St-apll for the validity of the Stokes approximation we find

$`{\displaystyle \frac{\tau _\text{p}}{\tau _\eta }}\left({\displaystyle \frac{\rho _\text{p}}{\rho _\text{f}}}\right)^{1/3},`$ (7)
where we neglected the difference between $`8/9`$ and 1. Equation LABEL:t-rat means that for “heavy” particles in a gas, that satisfy Stokes approximation, the particle response time scale may be about ten times larger than the Kolmogorov time scale: $`\tau _\text{p}10\tau _\eta `$. For such particles in a liquid the two time scales are about the same. So we may conclude that for “heavy” particles in a gas, that satisfy Stokes approximation, the inertia of the particles may be expected to be important in a considerable part of the energy spectrum. For particles in a liquid the particle inertia will only be significant for the smallest eddies, for which $`\tau _\text{p}\tau _\eta `$.

## II  Review of literature

This section is devoted to a review of the literature about the problem of a turbulently flowing suspension. We will review important findings from published numerical experiments, physical experiments and analytical models. SSDNSReview of some DNS simulations

To study the two-way coupling effect several groups have applied the direct-numerical-simulation technique (DNS) to particle-laden isotropic turbulence. A brief review of some of the publications is given below.

Squires and Eaton 90SE  considered the particle motion in the Stokes regime in which gravitational settling was neglected. They assumed statistically stationary isotropic turbulence. Mass loadings from zero to unity were considered for a series of particle response times varying from $`0.3\tau _\eta `$ to $`11\tau _\eta `$, where $`\tau _\eta `$ is the Kolmogorov time scale. They found that the overall reduction in turbulence kinetic energy for increasing mass loading was insensitive to the particle response time. They attributed the non-uniform distortion of the turbulent energy spectrum by particles to the preferential concentration of particles into regions of low vorticity and/or high strain rate.

Elghobashi and Truesdell 93ET  examined turbulence modulation by particles in decaying isotropic turbulence. They used the particle equation of motion derived by Maxey and Riley 83MR , and found that for the large density ratio considered in their simulations the particle motion was influenced mostly by drag and gravity. They found that the coupling between particles and fluid resulted in an increase in small-scale energy. This increase in the energy of the high-wave-number components of the velocity field resulted in a larger dissipation rate. They also found that the effect of gravity resulted in an anisotropic modulation of the turbulence and an enhancement of turbulence energy levels in the direction aligned with gravity.

Boivin, Simonin and Squires 98BSS  also made a very detailed DNS-study of the modulation of isotropic turbulence by particles. The focus of their work was on the class of dilute flows in which particle volume fractions and inter-particle collisions are negligible. Gravitational settling was also neglected and the particle motion was assumed to be governed by drag with particle response times ranging from the Kolmogorov scale to the Eulerian time scale of the turbulence and particle mass loadings up to unity. The velocity field was made statistically stationary by forcing the low wave-numbers of the flow. Like in 90SE ; 93ET  the effect of particles on the turbulence was included by using the point-force approximation. The DNS-results showed that particles increasingly dissipate fluid kinetic energy with increased mass loading, with the reduction in kinetic energy being relatively independent of the particle response time (as was already found in 90SE ). The viscous dissipation in the fluid decreases with increased mass loading and is larger for particles with smaller response times. The fluid energy spectra show that there is a non-uniform distortion of the turbulence spectrum with a relative increase in small-scale energy (as was found in 93ET ). They state that the fluid drags the particles at low wave-numbers while the converse is true at high wave-numbers for small particles.

Sundaram and Collins 99SC  performed DNS-simulations of particle-laden isotropic decaying turbulence. The particle response time was in the range $`1.6\tau _\eta \tau _p6.4\tau _\eta `$. The ratio of the particle density and fluid density was of the order $`10^3`$. The particle Reynolds number $`e_\text{p}`$ remained less than $`e_{\text{cr}}`$, and the drag force on the particles was described by Stokes law. The point-force approximation was employed to represent the two-way coupling force in the fluid momentum equation. The DNS-results showed that the particles reduce the turbulent kinetic energy as compared to the particle-free case, and this reduction is less pronounced for smaller response times $`\tau _p`$. The results also showed that the total turbulent energy dissipation is increased by the particles, and the increase is larger for smaller $`\tau _p`$. The turbulent energy spectrum is reduced at small wave-numbers and increased at high wave-numbers by the two-way coupling, and the location of the cross-over point is shifted towards larger wave-numbers for larger $`\tau _p`$.

Druzhinin 01Dru  examined the modulation of isotropic decaying turbulence by microparticles, for which $`2a<\eta `$, $`\tau _\text{p}<\tau _\eta `$ and $`e_\text{p}<e_{\text{cr}}`$. The gravitational settling is neglected. Due to the fact that $`\rho _\text{p}\rho _\text{f}`$, the mass loading may be large enough to modify the carrier flow. Druzhinin first derived an approximate analytical solution for the energy spectrum and then performed also DNS simulations. The results obtained for particles whose $`\tau _\text{p}0.4\tau _\eta `$ show that both the turbulence kinetic energy and the turbulence dissipation rate are increased by the two-way coupling effect as compared to the particle-free case. For particles with sufficiently high inertia ($`\tau _\text{p}0.5\tau _\eta `$) the two-way coupling effect caused a reduction in the turbulence kinetic energy as compared to the particle free case. Druzhinin, therefore, showed that there occurs a qualitative transition in the two-way coupling effect of particles on isotropic turbulence as the particle response time is increased from $`\tau _\text{p}\tau _\eta `$ (microparticles) to $`\tau _\text{p}\tau _\eta `$ (particles with finite inertia). For microparticles there is an increase of all wave-numbers in the energy spectrum. For particles with a higher inertia that is no longer the case.

SSlab-expReview of some laboratory experiments

Many experiments have been carried out to study the modulation of turbulence in the carrier phase by particles. An overview of the experimental work up to 1999 is given by Pietryga 01Pie . Experimental measurements in shear flows, e.g. particle-laden jets and boundary layers, have shown that turbulence velocity fluctuations may be either increased or decreased due to the modulation of the flow by (heavy) particles. However in turbulent shear flows it is often difficult to separate the direct modulation of the turbulence due to the momentum exchange with the particles from the indirect changes occurring through modification of turbulence production mechanisms via interaction with mean gradients. In grid-generated turbulence these production mechanisms are absent. It approximates in the best possible way the homogeneous, isotropic turbulence with particles that we study in this publication. We will, therefore, briefly review below some literature publications about experimental work devoted to the study of the modulation by particles of grid-generated turbulence.

Schreck and Kleis 93SK  studied the effect of almost neutrally buoyant plastic particles (density 1050 kg/m<sup>3</sup>) and heavy glass particles (density 2400 kg/m<sup>3</sup>) on grid-generated turbulence in a water flow facility. The average particle size was $`655\mu `$m. The particle Reynolds number of the plastic particles $`e_\text{p}8`$, for the glass particles $`e_\text{p}20`$. The particle volume fraction was varied between $`0.4\%`$ and $`1.5\%`$, so the system was very dilute. Mean velocity and velocity fluctuations of both phases were measured by a laser-Doppler velocimeter. The presence of the particles in sufficiently high concentration modified the turbulence downstream of the grid. The decay rate of the turbulence energy increased monotonically with particle concentration. The additional dissipation rate for the suspensions with the heavier glass particles was about double that of the almost neutrally buoyant plastic particles. A simple model based on the slip velocity between the phases under-predicted the measured increase in the dissipation rate. Schreck and Kleis, therefore, assumed that a large portion of the additional dissipation is associated with the measured modification of the spectral distribution of the turbulence energy. They speculate that the particles enhance the transfer of energy to smaller eddies extending the dissipation spectrum to smaller scale. Since only part of the high wave-number end of the spectrum could be resolved experimentally, this speculation could not be conclusively demonstrated by their experimental data.

Hussainov et al. 00HKR  studied the modulation of grid-generated turbulence by coarse glass particles in a vertical downward channel flow of air. Two different types of grids were used. Glass beads with an average diameter of $`700\mu `$m and a mass loading of $`10\%`$ were used. The particles were about 7 times larger than the Kolmogorov length scale $`\eta `$ and $`e_\text{p}70`$ or 93, dependent on the type of grid used. The particle response time scale of the particles $`\tau _\text{p}`$ was about 5000 to 7000 times larger than the Kolmogorov time scale $`\tau _\eta `$. The mean velociy and the turbulence intensity along the channel axis (and in some cross-sections) were measured by means of a laser-Doppler velocimeter. The decay curves of the turbulence intensity in the streamwise direction showed an attenuation of turbulence intensity of the flow by the particles. The particles caused an increase in the total dissipation rate of the turbulence. Hussainov et al. found that the presence of the particles decreased the energy spectra at high frequencies. This seems to be in contradiction with the speculation of Schreck and Kleis, that the particles enhance the transfer of energy to smaller eddies. SSmodelsAnalytical models

The starting point for analytical models, described in the literature, is often the Navier-Stokes (NS) equation for the velocity of the pure fluid (fluid without particles) $`𝒖(t,𝒓)`$

$`\rho _\text{f}\left[{\displaystyle \frac{}{t}}+\left(𝒖\mathbf{}\right)\nu ^2\right]𝒖+\mathbf{}p=𝒇_\text{p}+𝒇,`$ (8)
where $`p(t,𝒓)`$ is the pressure and $`\rho _\text{f}`$ is the fluid density. The random vector field $`𝒇(t,𝒓)`$ represents the stirring force responsible for the maintenance of the turbulent flow. Equation LABEL:NSf includes also the force $`𝒇_\text{p}(t,𝒓)`$ caused by the friction of the fluid with particles:

$`𝒇_\text{p}(t,𝒓)={\displaystyle \frac{\varphi \rho _\text{f}}{\tau _\text{p}}}\left[𝒗(t,𝒓)𝒖(t,𝒓)\right].`$ (9)
Here $`𝒗(t,𝒓)`$ is the velocity field of the particles, considered as a continuous medium with density $`m_\text{p}/\mathrm{}^3=\rho _\text{f}\varphi `$, where $`m_\text{p}`$ is the mass of a particle, $`\mathrm{}^3`$\- suspension volume per particle and $`\varphi `$ is the mass loading parameter

$`\varphi =m_\text{p}/\rho _\text{f}\mathrm{}^3.`$ (10)
The validity to represent $`𝒇_\text{p}(t,𝒓)`$ in the form of LABEL:friction is based on the assumption of space homogeneity of the particle distribution. It is also assumed that the particles are small enough for the Stokes drag law to be valid. The equation of motion, suggested in the literature, for the continuum phase of the particles does often not include the pressure and viscous terms

$`\left({\displaystyle \frac{m_\text{p}}{\mathrm{}^3}}\right)\left[{\displaystyle \frac{}{t}}+\left(𝒗\mathbf{}\right)\right]𝒗=𝒇_\text{p}.`$ (11)
Equations LABEL:NSf and LABEL:NSp were used by Baw and Peskin 71BP  to derive a set of “energy balance” equations for the following functions:

* $`E_{\text{ff}}\left(k\right)`$ – energy spectrum of the fluid turbulence, $`E\left(k\right)`$ in our nomenclature
* $`E_{\text{ff,p}}\left(k\right)`$ – energy spectrum of the fluid turbulence along a particle trajectory
* $`E_{\text{fp}}\left(k\right)`$ – fluid-particle covariance spectrum
* $`E_{\text{pp}}\left(k\right)`$ – particle energy spectrum

In the balance equations the following energy transfer functions occur

* $`T_{\text{ff,f}}\left(k\right)`$ – energy transfer in fluid turbulence
* $`T_{\text{fp,f}}\left(k\right)`$ – transfer of fluid-particle correlated motion by the fluid turbulence along the particle path
* $`T_{\text{fp,p}}\left(k\right)`$ – transfer of fluid-particle correlated motion by the particles
* $`T_{\text{pp}}\left(k\right)`$ – transfer of particle-particle correlated motion by the particle motion

* $`\mathrm{\Pi }_{\text{q,f}}\left(k\right)`$ – fluid-particle energy exchange rate.

Baw and Peskin 71BP  made a set of simplifying assumptions in order to be able to analyze the balance equations. First, they assumed that the particles do not respond to the fluid velocity fluctuations due to their (very large) inertia. Therefore

$$E_{\text{ff,p}}\left(k\right)=E_{\text{ff}}\left(k\right),$$
(12)
$$T_{\text{fp,f}}\left(k\right)=T_{\text{fp,p}}\left(k\right)=T_{pp,p}\left(k\right)=0.$$

This assumption is, of course, not realistic for particles satisfying the Stokes’ approximation. Their next assumption

$$\mathrm{\Pi }_{\text{q,f}}=\varphi \left[E_{\text{fp}}\left(k\right)E_{\text{ff,p}}\left(k\right)\right]/\tau _\text{p},$$
(13)
may be understood as a statement that the fluid-particle exchange rate is statistically the same for all scales characterized by a $`k`$-independent frequency $`\gamma _\text{p}=\varphi /\tau _\text{p}`$. This is reasonable for particles with very large inertia, but then Stokes law is not valid. For particles satisfying Stokes law, assumption LABEL:BP2 has to be replaced with a more realistic, $`k`$-dependent frequency $`\gamma _\text{p}\left(k\right)`$. We will come back to this point when discussing our new model.

A serious difficulty in the derivation of the energy balance equations is how to find a closure expression for third-order velocity correlation functions, responsible for the various energy transfer functions. Baw and Peskin assumed that $`T_{\text{ff,f}}\left(k\right)`$ can be expressed similarly as in the case of a pure (single phase) flow

$$T_{\text{ff,f}}\left(k\right)=\frac{\text{d}}{\text{d}k}\frac{ϵ_f^{1/3}k^{5/3}E_{\text{ff}}\left(k\right)}{\alpha },$$
(14)
where $`ϵ_\text{f}`$ is the viscous dissipation in the pure fluid (without particles) and $`\alpha `$ is the so-called Kolmogorov constant. This assumption seems questionable to us. According to the spirit of the Richardson-Kolmogorov cascade picture of turbulence one may express inertial range objects, like $`T_{\text{ff,f}}\left(k\right)`$ in terms of again inertial range quantities, like $`k`$, $`E_{\text{ff}}\left(k\right)`$ (which is done in Eq. (14)) and $`\epsilon \left(k\right)`$, the energy flux in $`k`$-space. In a single phase flow, indeed $`\epsilon \left(k\right)=ϵ_\text{f}`$. However this is not the case for a turbulent suspension due to the fluid-particle energy exchange, given by Eq. (13). We think that our closure (to be discussed later on) is an improvement in this respect.

With this simplified model Baw and Peskin predicted the following influences on the energy spectrum of the fluid turbulence due to the particles:

* a decrease of the energy in the energy-containing range of the spectrum
* an increase in the inertial range of the spectrum
* a decrease in the viscous dissipation range.

Boivin, Simonin and Squires 98BSS  used the same model as in Ref. 71BP, . They also applied assumptions similar to Eq. (13) and Eq. (14). Fortunately, they took into account the response of the particles to the turbulent velocity fluctuations by relaxing assumptions LABEL:BP1 and also accounted for the very important physical effect of the energy dissipation due to the drag around the particles. For that reason they approximated $`T_{\text{ff,f}}\left(k\right)`$ and $`T_{\text{fp,f}}\left(k\right)`$ as follows:

$$T_{\text{ff,f}}\left(k\right)=\frac{\text{d}}{\alpha \text{d}k}\left[ϵ_\text{f}\mathrm{\Pi }_{\text{q,f}}\left(k\right)\right]^{1/3}k^{5/3}E_{\text{ff}}\left(k\right),$$
(BSS1)
$$T_{\text{fp,f}}\left(k\right)=\frac{\text{d}}{\alpha \text{d}k}\left[ϵ_\text{f}\mathrm{\Pi }_{\text{q,f}}\left(k\right)\right]^{1/3}k^{5/3}E_{\text{fp}}\left(k\right).$$

Notice that this closure has the same weakness as Eq. (14), involving the dissipation range value $`ϵ_\text{f}`$. With the above described changes with respect to the model described in Ref. 71BP,  Boivin, Simonin and Squires found an increase in the viscous dissipation range of the fluid turbulence spectrum for small values of the particle response time $`\tau _\text{p}`$.

Al Taweel 77AT  calculated the rate of additional energy dissipation due to the presence of the particles. Because of their inertia the particles were assumed not to follow completely the turbulent velocity fluctuations of the carrier fluid. They expressed the additional dissipation in terms of the turbulent kinetic energy of the suspension. Then they added this term to the balance equation of the turbulent kinetic energy, making the (questionable) assumption that the energy flux across the spectrum has the same functional form as in a single-phase flow. Solving this equation they found an attenuation of the high-frequency fluctuations with a small alteration of the energy-containing low frequencies. Although there was an additional energy dissipation due to the particles, the total energy dissipation was reduced due to the reduction of viscous dissipation in the carrier fluid.

In a number of publications B5 ; B6 ; B7 ; B11  Felderhof, Ooms and Jansen developed an analytical model for the dynamics of a suspension of solid spherical particles in an incompressible fluid based on the linearized version of the Navier-Stokes equation. In particular they studied the effect of the particles-fluid interaction on the effective transport coefficients and on the turbulent energy spectrum of the suspension. Also the hydrodynamic interaction between the particles and the influence of the finite size of the particles were incorporated. However it is needless to say that the nonlinearity of the Navier-Stokes equation is of crucial importance in the problem of turbulence. Felderhof, Ooms and Jansen were well aware of this problem, but wanted to study in particular the influences of the particle-particle hydrodynamic interaction and of the finite particle size at a high particle volume concentration.

## III One-fluid model Navier-Stokes Equation for turbulent suspensions

In Sec. LABEL:ss:models we discussed the two-fluid model of suspensions consisting of the Navier-Stokes (NS) equation LABEL:NSf for the fluid and Eq. LABEL:NSp for the ”gaseous” phase of particles. This approximation is based on the assumptions of space homogeneity of the particle distribution and applicability of the Stokes drag law for the fluid-particle friction. We think that the basic physics of a turbulently flowing suspension with these assumptions may be described in the framework of the much more simple one-fluid equation. This model is presented in Sec. LABEL:ss:adv, discussed in Sec. LABEL:ss:NS-physics and “derived” in Secs. LABEL:ss:NSo-derivation and LABEL:ss:NSk-derivation.

SSadvThe model The following equation may be considered as a model equation for turbulently flowing suspensions:

$`\rho _{\text{eff}}\left(k\right)\left[{\displaystyle \frac{}{t}}+\gamma _\text{p}\left(k\right)+\gamma _0\left(k\right)\right]𝒖(t,𝒌)`$
$`=`$ $`𝓝\{𝒖,𝒖\}_{t,𝒌}+𝒇(t,𝒌),`$
The linear part of this equation involves:

$`\rho _{\text{eff}}\left(k\right)=\rho _\text{f}\left\{1\psi +\varphi {\displaystyle \frac{1+2\tau _\text{p}\gamma \left(k\right)}{\left[1+\tau _\text{p}\gamma \left(k\right)\right]^2}}\right\},`$ (16)
$`\gamma _\text{p}\left(k\right)={\displaystyle \frac{\varphi \tau _\text{p}\left[\gamma \left(k\right)\right]^2}{\left(1+\varphi \right)\left[1+2\tau _\text{p}\gamma \left(k\right)\right]+\left[\tau _\text{p}\gamma \left(k\right)\right]^2}},`$ (17)
$`\gamma _0\left(k\right)=\nu _{\text{eff}}\left(k\right)k^2,\nu _{\text{eff}}\left(k\right)={\displaystyle \frac{\nu \rho _\text{f}}{\rho _{\text{eff}}\left(k\right)}}.`$ (18)
The nonlinear term in Eq. (III) has the usual NS equation form:

$`𝒩\{𝒖,𝒖\}_{t,𝒌}^\alpha `$ $`=`$ $`{\displaystyle \frac{\text{d}^3k_1\text{d}^3k_2}{\left(2\pi \right)^3}\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }u_\beta ^{}(t,𝒌_1)u_\gamma ^{}(t,𝒌_2)}.`$
However the vertex $`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$ differs from the standard vertex $`\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$ of the NS equation (see e.g. Refs. B8, ; B4, ):

$`\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$ $`=`$ $`{\displaystyle \frac{\rho _\text{f}}{2}}\left[P^{\alpha \beta }\left(𝒌\right)k^\gamma +P^{\alpha \gamma }\left(𝒌\right)k^\beta \right]\delta \left(𝒌+𝒌_1+𝒌_2\right),`$
as follows:

$`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }=\rho _{\text{eff}}\left({\displaystyle \frac{2k_1k_2k_3}{k_1^2+k_2^2+k_3^2}}\right){\displaystyle \frac{\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }}{\rho _\text{f}}}.`$ (21)
Our model differs from the standard NS equation in three aspects:

a. Equation Eq. (III) involves the $`k`$-dependent effective density of suspensions $`\rho _{\text{eff}}\left(k\right)`$ given by Eq. (16). The function $`\rho _{\text{eff}}\left(k\right)`$ satisfies the inequality $`\rho _\text{f}\rho _{\text{eff}}\left(k\right)\rho _\text{f}\left(1+\varphi \right)`$. One could say that $`\rho _{\text{eff}}\left(k\right)\rho _\text{f}`$ represents the contribution of the particles involved in turbulent fluctuations with characteristic scale $`1/k`$ to the effective density of suspensions.

b. Equation LABEL:NS-basic includes the additional damping term $`\gamma _\text{p}\left(k\right)`$, Eq. (17), describing the loss of kinetic energy caused by the viscous fluid-particle friction.

c. In the absence of a stirring force $`𝒇(t,𝒓)`$ and both damping terms, Eq. (III) conserves the total kinetic energy of suspensions $``$ \[given by Eq. (56)\] which is different from the kinetic energy $`E`$ of the fluid itself.

The explicit form (LABEL:eq:NL-basic) of the nonlinear term is not necessary for the simple closure procedure that we applied in this publication. For the introduction of the energy flux in used closure procedure it is enough to use the fact that the modelled nonlinearity must be conservative. However, the explicit form is needed for more advanced closure procedures that we intend to use in future work. For this reason we include it in this publication.

SSNS-physicsPhysical interpretation of the one-fluid model

In a simplified fashion we may interpret $`\rho _{\text{eff}}\left(k\right)`$, the $`k`$-dependent density of suspension in our model equation LABEL:NS-basic as follows.

Denote as $`f_{\text{com}}\left(k\right)`$ the fraction of particles co-moving with the $`k`$-eddies (turbulent fluctuation with some wave-number $`k`$), in the sense that their velocity is the same as the velocity of $`k`$-eddies. These particles also participate in the motion of eddies with smaller wave-number $`k^{}<k`$ but not necessarily in the motion of $`k^{\prime \prime }`$-eddies with $`k^{\prime \prime }>k`$. For small $`k`$ the turnover frequency $`\gamma \left(k\right)`$ of $`k`$-eddies is small in the sense $`\gamma \left(k\right)\tau _\text{p}1`$. Therefore, in this region of $`k`$, the particle velocity is very close to that of the carrier fluid and we can describe the suspension as a single fluid with effective density $`\rho _{\text{eff}}\left(k\right)`$, which is close to the density of suspension:

$`\rho _\text{s}`$ $`=`$ $`\rho _\text{f}\left(1\psi \right)+C_\text{p}\rho _\text{p}=\rho _\text{f}\left(1\psi +\varphi \right),`$ (22)
$`\psi `$ $``$ $`C_\text{p}\left[4\pi a^3/3\right],\varphi =C_\text{p}\rho _\text{p}/\rho _\text{f}.`$
Here $`C_\text{p}`$ is the particle concentration, $`\psi `$ and $`\varphi `$ are the volume fraction and mass loading parameter. However, for large $`k`$, when $`\gamma \left(k\right)\tau _\text{p}1`$, the particles cannot follow these very fast motions and may be considered at rest. Thus the particles do not contribute to the effective density and $`\rho _{\text{eff}}\left(k\right)\rho _\text{f}`$. In the general case $`\rho _{\text{eff}}\left(k\right)`$ may be written as

$`\rho _{\text{eff}}\left(k\right)=\rho _\text{f}\left[1\psi +\varphi f_{\text{com}}\left(k\right)\right],`$ (23)
Here a statistical ensemble of all particles, partially involved in the motion of $`k`$-eddies, is replaced by two sub-ensembles of “fully co-moving” (fraction $`f_{\text{com}}\left(k\right)`$ ) and “fully at rest” (fraction $`f_{\text{rest}}\left(k\right)=1f_{\text{com}}\left(k\right)`$) particles, which does not contribute to $`\rho _{\text{eff}}\left(k\right)`$.

The particles at rest cause the fluid-particle friction. According to Newton’s third law, the damping frequency of a suspension $`\gamma _\text{p}\left(k\right)`$ may be related to the particle response time, $`\tau _\text{p}`$, via the ratio of total mass of particles $`M_\text{p}`$ at rest to the total effective mass of the suspension $`M_{\text{eff}}\left(k\right)`$:

$`\gamma _\text{p}\left(k\right)={\displaystyle \frac{M_\text{p}}{\tau _\text{p}M_{\text{eff}}\left(k\right)}}={\displaystyle \frac{C_\text{p}\rho _\text{p}f_{\text{rest}}\left(k\right)}{\tau _\text{p}\rho _{\text{eff}}\left(k\right)}}={\displaystyle \frac{\varphi \rho _\text{f}f_{\text{rest}}\left(k\right)}{\tau _\text{p}\rho _{\text{eff}}\left(k\right)}}.`$ (24)
As we mentioned, the fractions $`f_{\text{com}}\left(k\right)`$ and $`f_{\text{rest}}\left(k\right)`$ depend on $`\tau _\text{p}\gamma \left(k\right)`$. Moreover, the portion $`f_{\text{rest}}\left(k\right)`$ is independent on the sign of the velocity, therefore we expect $`f_{\text{rest}}\left(k\right)\left[\tau _\text{p}\gamma \left(k\right)\right]^2`$. In the opposite case, when $`1/\tau _\text{p}\gamma \left(k\right)`$ is small, $`f_{\text{com }}\left(k\right)`$ has corresponding smallness: $`f_{\text{com}}\left(k\right)1/\tau _\text{p}\gamma \left(k\right)`$. As a simple model of such a function we adopt

$$f_{\text{rest}}\left(k\right)=1f_{\text{com}}\left(k\right)=\left[\tau _\text{p}\gamma \left(k\right)/\left(1+\tau _\text{p}\gamma \left(k\right)\right)\right]^2.$$
(25)
Using Eq. (25), we rewrite Eqs. LABEL:eff-den and LABEL:gamma-p1 as Eqs. LABEL:r-eff and LABEL:gamma-p. Note that these equations, which follow from the physical reasoning described above, give the same expression for $`\gamma _\text{p}\left(k\right)`$ as Eq. (17) in our “derivation” in Sec. III.0.2. We consider this fact as a strong support of the physical relevance of our one-fluid model for a turbulently flowing suspension given by Eqs. LABEL:NS-basic-LABEL:nu-eff, with $`k`$-dependent effective density, fluid-particle damping frequency $`\gamma _\text{p}`$, and effective kinematic viscosity $`\nu _{\text{eff}}\left(k\right)`$.

SSassBasic assumptions

The theory developed in this paper is based on a number of assumptions and simplifications described below:

1. All particles in the suspension are spheres with the same density $`\rho _\text{p}`$ and the same radius $`a`$.
2. The radius of the particles is small enough, see Eq(1).
3. The particle-particle interaction will be neglected, assuming that the volume fraction $`\psi 1`$. Nevertheless, for the very heavy particles with $`\rho _\text{p}\rho _\text{f}`$, the mass loading $`\varphi `$ may be of the order of unity, leading to a significant modification of the turbulent flow by particles.
4. The turbulent flow is stationary, homogeneous and isotropic.
5. In our equations for the energy balance LABEL:budget we will use simple (but physically relevant) closure procedures based on our effective (one-fluid) Navier-Stokes equation for suspensions LABEL:NS-basic and on the Richardson-Kolmogorov cascade picture of turbulence.

SSNSo-derivationFormal derivation of the effective NS equation for suspensions

In the derivation we begin with the NS Eq. (8) for the fluid component. Instead of the averaged expression LABEL:friction for the fluid-particle friction force we will use the following detailed expression

$`𝒇_\text{p}(t,𝒓)={\displaystyle \underset{j}{}}𝑭_\text{p}(t,𝒓_j)\delta \left(𝒓𝒓_j\right),`$ (26)
in which $`𝑭_\text{p}(𝒓_j,t)`$ is the force between the fluid and $`j`$-particle positioned at $`𝒓=𝒓_j`$. Assume \[as in derivation of Eqs. LABEL:friction and LABEL:NSp\] that the statistics of particles is independent of the statistics of turbulence and, moreover, that their distribution is space homogeneous. In that case, we can replace the sum over the position of particles by a space integration:

$$\underset{j}{}\frac{1}{\mathrm{}^3}\text{d}𝒓_j,$$

where $`\mathrm{}^3`$ is the volume per particle. In this approximation

$`𝒇_\text{p}(𝒓,t)=𝑭_\text{p}(𝒓,t)/\mathrm{}^3.`$ (27)
We compute $`𝑭_\text{p}(t,𝒓)`$ for small enough particles with a radius $`a`$ satisfying inequality LABEL:St-apll, such that the fluid flow in the vicinity of a particle may be considered as laminar (assumption LABEL:ss:ass-2). Then one can apply Stokes law for the force $`𝑭_\text{p}(t,𝒓)`$:

$`𝑭_\text{p}(t,𝒓)=\zeta \left[𝒗_\text{p}\left(t\right)𝒖(t,𝒓)\right]`$ (28)
with the friction coefficient $`\zeta `$ for heavy particles (with the density $`\rho _\text{p}\rho _\text{f}`$) is given by

$`\zeta `$ $`=`$ $`6\pi \rho _\text{f}\nu a.`$ (29)
The Newton equation for a particle reads:

$`m_\text{p}{\displaystyle \frac{\text{d}𝒗_\text{p}\left(t\right)}{\text{d}t}}=𝑭_\text{p}(t,𝒓)=\zeta \left[𝒖(t,𝒓)𝒗_\text{p}\left(t\right)\right].`$ (30)
A formal solution of this equation

$`𝒗_\text{p}\left(t\right)=\left[\tau _\text{p}{\displaystyle \frac{\text{d}}{\text{d}t}}+1\right]^1𝒖(t,𝒓),`$ (31)
allows one to express the force $`𝑭_\text{p}(t,𝒓)`$ as follows

$`𝑭_\text{p}(t,𝒓)=m_\text{p}{\displaystyle \frac{\text{d}}{\text{d}t}}\left[\tau _\text{p}{\displaystyle \frac{\text{d}}{\text{d}t}}+1\right]^1𝒖(t,𝒓).`$ (32)
Here $`\tau _\text{p}`$ is the particle response time:

$`\tau _\text{p}=m_\text{p}/6\pi \nu \rho _\text{f}a.`$ (33)
The total time derivative $`\left(\text{d}/\text{d}t\right)`$ as usual takes into account the time dependence of the particle coordinate $`𝒓`$:

$`{\displaystyle \frac{\text{d}}{\text{d}t}}=\left[{\displaystyle \frac{}{t}}+𝒗_\text{p}\left(t\right)\mathbf{}\right].`$ (34)
Due to the (immersed) particle inertia they do not follow Lagrangian trajectories of fluid particles. Therefore, generally speaking, $`\left(\text{d}/\text{d}t\right)`$ does not coincide with the total Lagrangian time derivative in the fluid,

$`{\displaystyle \frac{D}{Dt}}=\left[{\displaystyle \frac{}{t}}+𝒖(t,𝒓)\mathbf{}\right].`$ (35)
Consider the relationship between $`\left(\text{d}/\text{d}t\right)`$ and $`\left(D/Dt\right)`$:

$`{\displaystyle \frac{D𝒖(t,𝒓)}{Dt}}={\displaystyle \frac{\text{d}𝒖(t,𝒓)}{\text{d}t}}+\left[v_\text{p}^\alpha u^\alpha (t,𝒓)\right]_\alpha 𝒖(t,𝒓)`$
$`=`$ $`{\displaystyle \frac{\text{d}𝒖(t,𝒓)}{\text{d}t}}{\displaystyle \frac{\text{d}}{\text{d}t}}{\displaystyle \frac{\tau _\text{p}}{1+\tau _\text{p}\frac{\text{d}}{\text{d}t}}}u(t,𝒓)^\alpha _1^\alpha 𝒖_1`$ (37)
$`=`$ $`{\displaystyle \frac{\text{d}}{\text{d}t}}{\displaystyle \frac{1}{1+\tau _\text{p}\frac{\text{d}}{\text{d}t}}}\left[\left(1+\tau _\text{p}{\displaystyle \frac{}{t}}\right)+\widehat{}\right]𝒖(t,𝒓),`$
$`\widehat{}𝒖(t,𝒓)\tau _\text{p}\left[\left(𝒗_{\text{p},1}\mathbf{}\right)𝒖(t,𝒓)\left(𝒖\mathbf{}_1\right)𝒖_1\right].`$
Here $`𝒖_1𝒖(t_1,𝒓_1)`$, $`\mathbf{}_1=\text{d}/\text{d}𝒓_1`$, and all derivatives with respect $`t_1`$ and $`𝒓_1`$ are taken at $`t_1=t`$ and $`𝒓_1=𝒓`$. Together with Eq. (35) this gives:

$`𝑭_\text{p}(t,𝒓)={\displaystyle \frac{D𝒖(t,𝒓)}{Dt}}{\displaystyle \frac{1}{1+\tau _\text{p}\frac{}{t}+\widehat{}}}𝒖(t,𝒓)`$ (38)
For particles with a small response time Ferry and Balachandar 01FB  show, that the particle velocity depends only on local fluid quantities (the velocity and its spatial and temporal derivatives). They derive an expansion of the particle velocity in terms of the particle response time which generalizes those of previous researchers. For large values of the ratio of the particle density and the fluid density and for small values of the particle response time our Eq. (38) for the force on a particle gives the same equation for the particle velocity as derived in Ref. 01FB .

Substitution of Eq. (38) into NS equation LABEL:NSf yields:

$`\rho _\text{f}[{\displaystyle \frac{}{t}}`$ $`+`$ $`(𝒖\mathbf{})\left]\right\{1+{\displaystyle \frac{\varphi }{1+\tau _\text{p}\left(\frac{}{t}+\widehat{}\right)}}\}𝒖+\mathbf{}p`$ (39)
$`=`$ $`\rho _\text{f}\nu ^2𝒖+𝒇,`$
where $`\varphi `$ is the mass loading parameter. For simplicity we consider here only the case of heavy particles with negligibly small volume loading parameter, $`\psi 1`$. However, the mass loading parameter may be of the order of unity. For example, for the water droplets in the air, $`\rho _\text{p}/\rho _\text{f}10^3`$ and for $`\varphi =1`$, the volume fraction $`\psi 10^3`$.

The inverse operator in Eq. LABEL:NS1 may be understood as a Taylor expansion with respect to the nonlinearity $`\left(𝒖\mathbf{}\right)`$:

$`{\displaystyle \frac{\varphi }{1+\tau _\text{p}\left(\frac{}{t}+\widehat{}\right)}}={\displaystyle \frac{1}{1+\tau _\text{p}\frac{}{t}}}{\displaystyle \frac{\widehat{}}{\left(1+\tau _\text{p}\frac{}{t}\right)^2}}+\mathrm{}`$
(40)
This expansion produces higher order \[ in $`\left(𝒖\mathbf{}\right)`$ \] nonlinear terms in the effective NS equation LABEL:NS1. These terms are not important for big eddies with $`\tau _\text{p}\gamma \left(k\right)1`$ for which the operator in the braces in the LHS of Eq. LABEL:NS1, $`\left\{\mathrm{}\right\}`$, is close to the factor $`1+\varphi `$. In the opposite case, for small scale eddies with $`\tau _\text{p}\gamma \left(k\right)1`$ the operator $`\left\{\mathrm{}\right\}=1`$. Both limiting cases one easily gets from the first term in the Taylor expansion LABEL:Taylor in which there is no contribution from $`\widehat{}`$. It means that only for intermediate scales with $`\tau _\text{p}\gamma \left(k\right)1`$ this operator may be quantitatively important. For a qualitative description of the “transient” process between these two regimes it is enough to account for the first term of expansion LABEL:Taylor. In this approximation the turbulent fluid velocity around the particle is approximated by the velocity at a fixed coordinate, which is reasonable in statistical sense and exact in the limit $`\tau _\text{p}\gamma \left(k\right)1`$. With this approximation Eq. LABEL:NS1 turns into

$`\widehat{\rho }_{\text{eff}}\left[{\displaystyle \frac{}{t}}+\left(𝒖\mathbf{}\right)\right]𝒖+\mathbf{}p=\rho _\text{f}\nu ^2𝒖+𝒇,`$ (41)
$`\widehat{\rho }_{\text{eff}}\rho _\text{f}\left\{1+{\displaystyle \frac{\varphi }{\tau _\text{p}\frac{}{t}+1}}\right\},`$ (42)
where $`\widehat{\rho }_{\text{eff}}`$ may be considered as an operator of effective density for a suspensions.

Since we are interested in the incompressible flows, we can project the potential components out of the equation of motion. This may be done by the projection operator $`𝓟`$, defined via its kernel $`𝒫^{\alpha \beta }\left(𝒓\right)`$:

$`𝒫^{\alpha \beta }\left(𝒓\right)`$ $``$ $`{\displaystyle \frac{\text{d}^3k}{\left(2\pi \right)^3}P^{\alpha \beta }\left(𝒌\right)\mathrm{exp}\left[i𝒌𝒓\right]},`$ (43)
$`P^{\alpha \beta }\left(𝒌\right)`$ $`=`$ $`\delta _{\alpha \beta }k_\alpha k_\beta /k^2.`$ (44)
The application of $`𝓟`$ to any given vector field $`𝒂\left(𝒓\right)`$ is nonlocal, and it has the form

$`\left[𝓟𝒂\left(𝒓\right)\right]_\alpha ={\displaystyle 𝑑𝒓_1𝒫^{\alpha \beta }\left(𝒓𝒓_1\right)a_\beta \left(𝒓_1\right)}.`$ (45)
Applying $`𝓟`$ to Eq. LABEL:NS2 we find

$`\widehat{\rho }_{\text{eff}}\left[{\displaystyle \frac{}{t}}+𝓟\left(𝒖\mathbf{}\right)\right]𝒖=\rho _\text{f}\nu ^2𝒖+𝒇,`$ (46)
This equation together with the definition (42) for the operator of the effective density constitutes a one-fluid description of a turbulently flowing suspension. However, the operator form of the effective density is not convenient for practical calculations. To overcome this inconvenience we will derive below another form for the effective parameters of this equation.

SSNSk-derivationNS equation for suspensions with $`k`$-dependent parameters

In our analytical description of space homogeneous, stationary turbulence it is convenient to consider Eq. (46) in the $`(𝒌,\omega )`$-representation:

$`\left\{\omega \left[\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)\right]\text{i}\rho _\text{f}\left[\nu k^2+\stackrel{~}{\gamma }_\text{p}\left(\omega \right)\right]\right\}\stackrel{~}{𝒖}(\omega ,𝒌)`$ (47)
$`=𝓝\{𝒖,𝒖\}_{\omega ,𝒌}+\stackrel{~}{𝒇}(\omega ,𝒌).`$
Here

$`\stackrel{~}{𝒖}(\omega ,𝒌)`$ $`=`$ $`{\displaystyle \text{d}t\text{d}𝒓𝒖(t,𝒓)\mathrm{exp}\left(\text{i}\omega t+\text{i}𝒌𝒓\right)},`$ (48)
$`\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)`$ $`=`$ $`\text{Re}\left\{\stackrel{~}{\rho }_{\text{eff}}\left(\omega \right)\right\}=\rho _\text{f}\left[1+{\displaystyle \frac{\varphi }{1+\left(\omega \tau _\text{p}\right)^2}}\right],`$ (49)
$`\stackrel{~}{\gamma }_\text{p}\left(\omega \right)`$ $`=`$ $`{\displaystyle \frac{\omega }{\rho _\text{f}}}\text{Im}\left[\left\{\stackrel{~}{\rho }_{\text{eff}}\left(\omega \right)\right]\right\}={\displaystyle \frac{\varphi \omega ^2\tau _\text{p}}{1+\left(\omega \tau _\text{p}\right)^2}},`$ (50)
$`\rho _{\text{eff}}\left(\omega \right)`$ $`=`$ $`\rho _\text{f}\left[1+{\displaystyle \frac{\varphi }{1\text{i}\omega \tau _\text{p}}}\right],`$
$`𝓝\{𝒖,𝒖\}_{\omega ,𝒌}\left[\stackrel{~}{\rho }_{\text{eff}}\left(\omega \right)𝓟\left(\stackrel{~}{𝒖}\mathbf{}\right)\stackrel{~}{𝒖}\right]_{\omega ,𝒌}.`$ (51)
$`𝓝\{𝒖,𝒖\}_{\omega ,𝒌}`$ denote the nonlinear term in $`(\omega ,𝒌)`$-representation and frequency $`\nu k^2`$ describes the viscous damping.

The Navier-Stokes equation for suspensions LABEL:NS-s involves a frequency-dependent effective density of suspensions, $`\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)`$ and $`\omega `$-dependent frequency $`\stackrel{~}{\gamma }_\text{p}\left(\omega \right)`$ responsible for the damping due to fluid-particle friction. To use standard closure procedures in the statistical description of turbulence one needs frequency independent coefficients in the basic equation of motion. On other hand, these procedures may be applied to equations with $`k`$-dependent coefficients. Therefore for further analysis it is much more convenient to deal with a $`k`$-dependent effective density $`\rho _{\text{eff}}\left(k\right)`$ of $`k`$-eddies. To relate these functions we note that the $`k`$-eddies have a characteristic frequency of motions, $`\gamma \left(k\right)`$ \[related to their life-time $`\tau \left(k\right)`$ by a simple relation $`\gamma \left(k\right)1/\tau \left(k\right)`$\].

#### III.0.1 $`k`$-dependent effective density of suspensions

In the inertial interval of scales Eq. (47) must preserve the total kinetic energy of a suspension $``$ if one neglects the fluid-particle friction $`\stackrel{~}{\gamma }_\text{p}\left(\omega \right)0`$. The equation for $``$ may be written in terms of the density $`[\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)`$ and $`\stackrel{~}{F}(\omega ,𝒌)`$, the pair correlation function of the $`(\omega ,𝒌)`$-Fourier components of velocity, $`\stackrel{~}{u}(\omega ,𝒌)`$. Namely

$`={\displaystyle \frac{\text{d}^3k}{\left(2\pi \right)^3}\frac{\text{d}\omega }{2\pi }\frac{\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)}{2}\stackrel{~}{F}(\omega ,𝒌)},`$ (52)
$`\left(2\pi \right)^4\delta \left(𝒌𝒌_1\right)\delta \left(\omega \omega _1\right)\stackrel{~}{F}(\omega ,𝒌)\stackrel{~}{𝒖}(\omega ,𝒌)\stackrel{~}{𝒖}(\omega _1,𝒌_1).`$
For isotropic turbulence $`\stackrel{~}{F}(\omega ,𝒌)=\stackrel{~}{F}(\omega ,k)`$ and Eq. (52) allows one to introduce the one-dimensional energy spectrum of suspension, $`\left(k\right)`$ according to:

$``$ $`=`$ $`{\displaystyle \frac{\text{d}k}{2\pi }\left(k\right)},`$ (53)
$`\left(k\right)`$ $`=`$ $`{\displaystyle \frac{\text{d}\omega }{\left(2\pi \right)}\frac{\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)}{2\pi }k^2\stackrel{~}{F}(\omega ,k)}.`$ (54)
Define a $`k`$-dependent effective density of suspension, which gives the same one-dimensional spectrum $`\left(k\right)`$ as the $`\omega `$-dependent effective density $`\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)`$:

$`\rho _{\text{eff}}\left(k\right)={\displaystyle \frac{\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)\stackrel{~}{F}(\omega ,k)\text{d}\omega }{\stackrel{~}{F}(\omega ,k)\text{d}\omega }}.`$ (55)
Then, Eq. (54) takes the form

$`\left(k\right)`$ $`=`$ $`{\displaystyle \frac{\rho _{\text{eff}}\left(k\right)}{2\pi }}k^2{\displaystyle \frac{\text{d}\omega }{\left(2\pi \right)}\stackrel{~}{F}(\omega ,k)}.`$ (56)
Notice that

$`{\displaystyle \frac{\text{d}\omega }{\left(2\pi \right)}\stackrel{~}{F}(\omega ,k)}=F\left(k\right),`$ (57)
is the simultaneous velocity pair correlation function. With this notations Eq. (56) may be written as

$`\left(k\right)`$ $`=`$ $`{\displaystyle \frac{\rho _{\text{eff}}\left(k\right)}{2\pi }}k^2F\left(k\right),`$ (58)
while the traditional notation for 1-dimensional spectrum of kinetic energy of fluid itself is $`E\left(k\right)`$:

$`E\left(k\right)`$ $`=`$ $`{\displaystyle \frac{\rho _\text{f}}{2\pi }}k^2F\left(k\right).`$ (59)
Formally speaking, in order to evaluate $`\rho _{\text{eff}}\left(k\right)`$ by Eq. (55) we need to know the $`\omega `$ dependence of $`F(\omega ,𝒌)`$. This is not a simple task. Instead we will use a few reasonable forms of $`F(\omega ,𝒌)`$ and compare the resulting functions $`\rho _{\text{eff}}\left(k\right)`$. One of the frequently used is the Lorentzian form

$`\stackrel{~}{F}(\omega ,k)=F\left(k\right){\displaystyle \frac{\gamma \left(k\right)/\pi }{\omega ^2+\gamma ^2\left(k\right)}},`$ (60)
which corresponds to the simplest “one-pole” approximation for the Greens functions. Using this $`\omega `$-dependence in Eq. LABEL:relation2 we have the following simple form for $`\rho _{\text{eff}}\left(k\right)`$.

$`\rho _{\text{eff}}\left(k\right)`$ $`=`$ $`\rho _\text{f}\left[1+{\displaystyle \frac{\varphi }{1+\tau _\text{p}\gamma \left(k\right)}}\right].`$ (61)
For small $`\tau _\text{p}\gamma \left(k\right)`$ this gives a correction linear in $`\tau _\text{p}\gamma \left(k\right)`$

$`\rho _{\text{eff}}\left(k\right)\rho _\text{f}\left[1\varphi \tau _\text{p}\gamma \left(k\right)\right],`$ (62)
which contradicts the physical intuition. Indeed, one may consider $`k`$-eddies as having a motion with the characteristic frequency $`\gamma \left(k\right)`$ and expect that $`\rho _{\text{eff}}\left(k\right)`$ may be obtained from $`\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)`$ with the substitution $`\omega \gamma \left(k\right)`$. This gives a correction quadratic in $`\tau _\text{p}\gamma \left(k\right)`$

$`\rho _{\text{eff}}\left(k\right)\rho _\text{f}\left(1+\varphi \right)\varphi \rho _\text{f}\tau _\text{p}^2\gamma ^2\left(k\right).`$ (63)
This contradiction follows from the fact that the model function Eq. (60) decays very slowly for $`\omega \mathrm{}`$, like $`1/\omega ^2`$. It is known from the diagrammatic analysis of the different time velocity correlation function $`F(\tau ,𝒌)`$ that for small $`\tau `$ the difference $`F(\tau ,𝒌)F(0,𝒌)`$ does not contains $`\left|\tau \right|`$ and decays not slower than $`\tau ^2`$. Therefore the Fourier transform of $`F(\tau ,𝒌)`$, $`\stackrel{~}{F}(\omega ,𝒌)`$ has to decay with $`\omega `$ faster than $`1/\omega ^2`$, at least as $`1/\omega ^4`$. To meet this requirement we consider instead of Eq. (60) the function:

$`\stackrel{~}{F}(\omega ,k)=F\left(k\right){\displaystyle \frac{2\gamma ^3\left(k\right)/\pi }{\left[\omega ^2+\gamma ^2\left(k\right)\right]^2}},`$ (64)
which gives instead of Eq. (61)

$`\rho _{\text{eff}}\left(k\right)`$ $`=`$ $`\rho _\text{f}\left\{1+{\displaystyle \frac{\varphi \left[1+2\tau _\text{p}\gamma \left(k\right)\right]}{\left[1+\tau _\text{p}\gamma \left(k\right)\right]^2}}\right\}`$
$`=`$ $`\rho _\text{f}\left(1+\varphi \right)\varphi \rho _\text{f}\left[{\displaystyle \frac{\tau _\text{p}\gamma \left(k\right)}{1+\tau _\text{p}\gamma \left(k\right)}}\right]^2.`$
Now the correction to $`\rho _{\text{eff}}\left(k\right)`$ is quadratic in $`\tau _\text{p}\gamma \left(k\right)`$ which agrees with the expectation LABEL:r-eff3. One observes the same agreement for any other model dependence $`\stackrel{~}{F}(\omega ,k)`$ decaying even faster than $`1/\omega ^4`$.

Therefore the linear part of Eq. (47) may be modelled as

$`\left\{\omega \rho _{\text{eff}}\left(k\right)\text{i}\rho _\text{f}\left[\nu k^2+\stackrel{~}{\gamma }_\text{p}\left(\omega \right)\right]\right\}\stackrel{~}{u}(\omega ,𝒌)=\mathrm{},`$ (66)
with $`\rho _{\text{eff}}\left(k\right)`$ given by Eq. (III.0.1).

#### III.0.2 Effective fluid-particle damping frequency $`\gamma _\text{p}(k)`$

Using Eq. (47) or Eq. (66) together with Eq. (54) we can compute the contribution of the fluid-particle friction to the damping of $`\left(k\right)`$:

$`{\displaystyle \frac{\left(k\right)}{t}}|_\text{p}=2\rho _\text{f}{\displaystyle \frac{\text{d}\omega }{\left(2\pi \right)}\frac{\stackrel{~}{\gamma }_\text{p}\left(\omega \right)}{2\pi }k^2\stackrel{~}{F}(\omega ,k)}.`$ (67)
Introduce an $`\omega `$-independent fluid-particle damping friction by a standard relation:

$`{\displaystyle \frac{\left(k\right)}{t}}|_\text{p}=2\gamma _\text{p}\left(k\right)\left(k\right).`$ (68)
Combining these two equations with Eq. (54) one gets

$`\gamma _\text{p}\left(k\right)={\displaystyle \frac{\rho _\text{f}\stackrel{~}{\gamma }_\text{p}\left(\omega \right)\stackrel{~}{F}(\omega ,k)\text{d}\omega }{\stackrel{~}{\rho }_{\text{eff}}^{}\left(\omega \right)\stackrel{~}{F}(\omega ,k)\text{d}\omega }}.`$ (69)
Substitution of Eqs. LABEL:gp-omega and LABEL:omega-dep into Eq. (69) gives Eq. (17) for $`\gamma _\text{p}\left(k\right)`$. With this knowledge, LABEL:NS-linear1 may be further simplified as follows

$`\rho _{\text{eff}}\left(k\right)\left\{\omega \text{i}\left[\nu _{\text{eff}}\left(k\right)k^2+\gamma _\text{p}\left(k\right)\right]\right\}\stackrel{~}{𝒖}(\omega ,𝒌)`$ (70)
$`=𝓝\{𝒖,𝒖\}_{\omega ,𝒌}+\stackrel{~}{𝒇}.`$
Here $`\nu _{\text{eff}}\left(k\right)`$ is given by Eq. (18). Notice that this equation gives the same dissipation rate LABEL:p-dam2 due to fluid-particle friction as Eq. (47) and the same dissipation rate

$`{\displaystyle \frac{\left(k\right)}{t}}|_\nu =2\nu _{\text{eff}}k^2\left(k\right)`$ (71)
due to the kinematic viscosity.

The suggested form of $`𝓝\{𝒖,𝒖\}_{\omega ,𝒌}`$ in terms of $`\rho _{\text{eff}}\left(k\right)`$ will be discussed in the following Sec. III.0.3.

#### III.0.3 $`\omega `$-independent nonlinearity of the effective NS equation

##### Nonlinearity in the usual NS equation.

Consider first the nonlinear term in the “usual” NS equation for single-phase flow. In $`(\omega ,𝒌)`$-representation it has the form (see e.g. Refs. B8, ; B4, ):

$`𝒩\{𝒖,𝒖\}_{\omega ,𝒌}^\alpha ={\displaystyle \frac{\text{d}^3k_1\text{d}^3k_2\text{d}\omega _1\text{d}\omega _2}{\left(2\pi \right)^4}}`$
$`\times `$ $`\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }u_\beta ^{}(\omega _1,𝒌_1)\stackrel{~}{u}_\gamma ^{}(\omega _2,𝒌_2)\delta \left(\omega +\omega _1+\omega _2\right).`$
Here $`\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$ is the so called vertex of interaction given by Eq. (LABEL:eq:vertex). It includes transversal projectors accounting for the incompressibility of the fluid, delta function of $`𝒌`$-vectors originating from the space homogeneity of the problem and is proportional to $`k`$ (as a reflection of $`\mathbf{}`$ operator in the nonlinear term in $`𝒓`$-representation).

The vertex $`\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$ satisfies so-called Jacobi identity

$`\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }+\gamma _{𝒌_2𝒌𝒌_1}^{\gamma \alpha \beta }+\gamma _{𝒌_1𝒌_2𝒌}^{\beta \gamma \alpha }=0,`$ (73)
as a consequence of the energy conservation by the Euler equation.

##### Nonlinearity in the effective NS Eq. (III).

A rigorous derivation of the nonlinear term in the effective NS Eq. (III) is a very delicate issue. For example, in Eq. (41) we used the operator of the effective density LABEL:rho-eff containing only the first term of expansion LABEL:Taylor. This approximation does not account for all terms, second orderin $`𝒖`$. This derivation is beyond the scope of this paper. Instead we present here physical arguments which allows us to propose a form of $`𝓝\{𝒖,𝒖\}_{\omega ,𝒌}`$ which satisfies all needed requirements.

By analogy with Eq. (III.0.3) we can write:

$`𝒩\{𝒖,𝒖\}_{\omega ,𝒌}^\alpha ={\displaystyle \frac{\text{d}^3k_1\text{d}^3k_2\text{d}\omega _1\text{d}\omega _2}{\left(2\pi \right)^4}}`$
$`\times `$ $`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }u_\beta ^{}(\omega _1,𝒌_1)\stackrel{~}{u}_\gamma ^{}(\omega _2,𝒌_2)\delta \left(\omega +\omega _1+\omega _2\right),`$
where the vertex $`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$ differs from $`\gamma _{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$, Eq. (LABEL:eq:vertex), because now $`\rho _\text{f}\rho _{\text{eff}}\left(k\right)`$.

The simplest possible generalization of the vertex, just a replacement $`\rho _\text{f}\rho _{\text{eff}}\left(k\right)`$ in Eq. (LABEL:eq:vertex) leads to a vertex $`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$, which does not satisfy the Jacobi identity

$`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }+\mathrm{\Gamma }_{𝒌_2𝒌𝒌_1}^{\gamma \alpha \beta }+\mathrm{\Gamma }_{𝒌_1𝒌_2𝒌}^{\beta \gamma \alpha }=0,`$ (75)
leading to violation of the conservation of the kinetic energy of suspension $``$. We suggest Eq. (21) for $`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$. Clearly, due to Eq. (73) this vertex satisfies requirement LABEL:Jac1 and consequently, Eq. (III) conserves the energy $``$.

Another physical requirement is Galilean invariance of the problem. This is the case for the standard NS equation with vertex LABEL:vertex in which $`\rho _\text{f}`$ is $`k`$-independent. For the $`k`$-dependent density in the vertex LABEL:vertex-basic Galilean invariance implies that in the limit, when one of the wave-numbers is much smaller then two others (say $`k_1k_2k_3`$), the effective density must depend on the smallest $`k`$-vector. Obviously, this is the case for the $`k`$-argument of $`\rho _{\text{eff}}\left(k\right)`$ in Eq. (21). This guarantees Galilean invariance of Eq. (III).

Now Eq. (70) involves only $`\omega `$-independent coefficients and may be rewritten in $`(t,𝒌)`$ representation, see Eq. (III).

## IV Budget of kinetic energy in turbulent suspensions

In this section we consider the budget of kinetic energy in turbulent suspension. In Sec. LABEL:ss:energy we will use the one-fluid model LABEL:NS-basic to derive (for homogeneous, isotropic turbulence) the following budget equation for the (1-dimensional) density of kinetic energy:

$`{\displaystyle \frac{(t,k)}{2t}}+\left[\gamma _0\left(k\right)+\gamma _\text{p}\left(k\right)\right](t,k)`$
$`=`$ $`𝒲(t,k)+𝒥(t,k).`$
The left hand side (LHS) of this equation includes two damping terms, $`\gamma _0\left(k\right)(t,k)`$, caused by the effective kinematic viscosity and $`\gamma _\text{p}\left(k\right)(t,k)`$ caused by the fluid-particle friction. The density $`(t,k)`$ is given by Eq. (82). The right hand side (RHS) of Eq. (IV) includes the source of energy $`𝒲(t,k)`$, localized in the energy containing interval, and the energy redistribution term, $`𝒥(t,k)`$.

The budget equation LABEL:budget is exact, but unfortunately not closed. Equations Eq. (18) for the effective kinematic viscosity and Eq. (17) for $`\gamma _\text{p}\left(k\right)`$ includes “turnover frequency” of $`k`$-eddies, $`\gamma \left(k\right)`$. Also $`𝒲(t,k)`$ contains yet unknown $`(u,f)`$ correlations, Eq. (81). And finally $`𝒥(t,k)`$ is given by Eqs. LABEL:CJk, LABEL:bal2 via 3-rd order velocity correlations $`F_3`$. There are many reasonable closure procedures for the approximation of high-order velocity correlations by lower-order ones. To elucidate the basic physics of the problem at hand, in this paper we will use the simplest possible closure. Application of more sophisticated closures will be done elsewhere.

SSenergyThe energy budget equation

In order to derive Eq. (IV) we multiply Eq. (III) by $`𝒖(t,𝒌^{})`$ and average:

$`\rho _{\text{eff}}\left(k\right)\left\{{\displaystyle \frac{F(t,𝒌)}{2t}}+\left[\gamma _0\left(k\right)+\gamma _\text{p}\left(k\right)\right]F(t,𝒌)\right\}`$
$`=J(t,𝒌)+W(t,𝒌),`$ (77)
$`J(t,𝒌){\displaystyle \frac{\text{d}^3k_1\text{d}^3k_2}{\left(2\pi \right)^3}\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }F_3^{\alpha \beta \gamma }(t;𝒌,𝒌_1,𝒌_2)}.`$ (78)
Here $`F(t,𝒌)`$, and $`F_3(t;\mathrm{})`$ are the 2nd and 3rd order simultaneous velocity correlation functions taken at overall time $`t`$:

$`\left(2\pi \right)^3\delta \left(𝒌+𝒌_1\right)F(t;𝒌)=𝒖(t,𝒌)𝒖(t,𝒌_1),`$ (79)
$`\left(2\pi \right)^3\delta \left(𝒌+𝒌_1+𝒌_2\right)F_3^{\alpha \beta \gamma }(t;𝒌,𝒌_1,𝒌_2)`$ (80)
$`=u^\alpha (t,𝒌)u^\beta (t,𝒌_1)u^\gamma (t,𝒌_2).`$
Note that the time $`t`$ in the argument of $`F\left(𝒌\right)`$ in Eqs. (57) – (60) is omitted: $`F(t,𝒌)=F\left(𝒌\right)`$.

In Eq. (77) we introduce also the simultaneous $`(u,f)`$-cross correlation functions

$`\left(2\pi \right)^3\delta \left(𝒌𝒌_1\right)W(t,𝒌)=𝒖(t,𝒌)𝒇(t,𝒌_1).`$ (81)
We can rewrite Eq. (56) for the density of the kinetic energy of suspension in terms of $`F(t,𝒌)=F(t,k)`$ (for isotropic turbulence):

$`(t,k)={\displaystyle \frac{\rho _{\text{eff}}\left(k\right)}{2\pi }}k^2F(t,k).`$ (82)
Multiplying Eq. (78) by $`k^2/2\pi `$ one gets finally the balance Eq. (IV) in which

$`𝒲(t,k)={\displaystyle \frac{k^2}{2\pi }}W(t,k),`$ (83)
$`𝒥(t,k)={\displaystyle \frac{k^2}{2\pi }}J(t,k).`$ (84)
Notice that effective vertex $`\mathrm{\Gamma }_{𝒌𝒌_1𝒌_2}^{\alpha \beta \gamma }`$ in Eq. (78) was constructed such that the total kinetic energy $``$ is the integral of motion (neglecting pumping and damping): $`_0^{\mathrm{}}𝒥(t,k)\text{d}k=0`$. Therefore the energy redistribution term $`𝒥(t,k)`$ may be written in the divergent form:

$`𝒥(t,k)`$ $`=`$ $`{\displaystyle \frac{\text{d}\epsilon (t,k)}{\text{d}k}},\text{where}`$ (85)
$`\epsilon (t,k)`$ $`=`$ $`{\displaystyle _k^{\mathrm{}}}\text{d}kJ(t,k)`$ (86)
is the (one-dimensional) energy flux over scales.

In the rest of the paper we will consider only stationary solutions of Eq. (IV). Omitting (here and below) time argument one finally has:

$`\left[\gamma _0\left(k\right)+\gamma _\text{p}\left(k\right)\right]\left(k\right)=𝒲\left(k\right){\displaystyle \frac{\text{d}\epsilon \left(k\right)}{\text{d}k}}.`$ (87)
SSclosureSimple closure for the energy budget equation

The effective density in our approach, Eq. (16), depends on the characteristic frequency $`k`$-eddies $`\gamma \left(k\right)`$. This object may be evaluated as the inverse life time of these eddies which is determined by their viscous damping and energy loss in the cascade processes. Accordingly, $`\gamma \left(k\right)`$ is a sum of two terms

$$\gamma \left(k\right)=\gamma _0\left(k\right)+\gamma _\text{c}\left(k\right),$$
(88)
where $`\gamma _0\left(k\right)`$, Eq. (18), is the viscous frequency and $`\gamma _\text{c}\left(k\right)`$ may be evaluated as the turnover frequency of $`k`$-eddies:

$`\gamma _\text{c}\left(k\right)kU_k,\text{where}U_k\sqrt{k\left(k\right)/\rho _{\text{eff}}\left(k\right)}`$ (89)
is the characteristic velocity of $`k`$-eddies. We therefore define

$`\gamma _\text{c}\left(k\right)=C_\gamma \sqrt{k^3\left(k\right)/\rho _{\text{eff}}\left(k\right)},`$ (90)
where $`C_\gamma `$ is some dimensional constant, presumably of the order of unity. Clearly, the same evaluation LABEL:gamma-k one gets from a dimensional reconstruction of $`\gamma \left(k\right)`$ in terms of the only relevant (in the K41 picture of turbulence) variables $`k`$, $`\left(k\right)`$ and $`\rho _{\text{eff}}\left(k\right)`$.

In the same manner, by the dimensional reasoning, one gets the following evaluation of the energy flux:

$`\epsilon \left(k\right)=C_\epsilon \sqrt{k^5^3\left(k\right)/\rho _{\text{eff}}\left(k\right)},`$ (91)
where $`C_\epsilon =O\left(1\right)`$. Notice that in pure fluids \[with $`\rho _{\text{eff}}\left(k\right)\rho _\text{f}`$\] Eqs. LABEL:gamma-k and LABEL:e-flux are nothing but the K41 evaluation of the corresponding objects. This become even more transparent if one rewrites Eq. LABEL:e-flux in the more familiar form:

$`\left(k\right)=C_1\left[\epsilon ^2\left(k\right)\rho _{\text{eff}}\left(k\right)\right]^{1/3}k^{5/3},C_1C_\epsilon ^{2/3}.`$ (92)
Together with Eq. (90) this gives a useful evaluation of $`\gamma _\text{c}\left(k\right)`$ via $`\epsilon \left(k\right)`$:

$`\gamma _\text{c}\left(k\right)=C_2\left[{\displaystyle \frac{\epsilon \left(k\right)k^2}{\rho _{\text{eff}}\left(k\right)}}\right]^{1/3},C_2{\displaystyle \frac{C_\gamma }{C_\epsilon ^{1/3}}}.`$ (93)
Lastly, we have to evaluate the energy input in the system $`𝒲\left(k\right)`$. It follows from Eq. (III) that $`u\left(k\right)`$ may be evaluated as $`f\left(k\right)/\gamma \left(k\right)\rho _{\text{eff}}\left(k\right)`$. Together with Eqs. LABEL:W, LABEL:CWk and LABEL:gamma-k this gives:

$`𝒲\left(k\right)=C_wf_k^2\sqrt{{\displaystyle \frac{k\left(k\right)}{\rho _{\text{eff}}\left(k\right)}}},`$ (94)
where $`C_w=O\left(1\right)`$ and $`f_k^2`$ is the pair correlation of the forcing \[which may be defined similarly to LABEL:W\].

SSbasic-eqDimensionless budget equation

For the convenience of the reader we present here the full set of equations which will be studied below. To non-dimensionalize this equation we define a dimensionless wave-number, $`\kappa `$, and the integral-scale related parameters

$$\kappa =kL,\epsilon __L=\epsilon \left(\frac{1}{L}\right),\gamma __L=\gamma \left(\frac{1}{L}\right),\rho __L=\rho _{\text{eff}}\left(\frac{1}{L}\right),$$
(95)
Define also the dimensionless functions

$`\epsilon _\kappa `$ $`=`$ $`{\displaystyle \frac{\epsilon \left(k\right)}{\epsilon __L}},\gamma _\kappa ={\displaystyle \frac{\gamma \left(k\right)}{\gamma __L}},`$ (96)
$`\rho _\kappa `$ $`=`$ $`{\displaystyle \frac{\rho _{\text{eff}}\left(k\right)}{\rho __L}},𝒲_\kappa ={\displaystyle \frac{𝒲\left(k\right)}{𝒲\left(1/L\right)}},`$
in which the argument $`\kappa `$ is written as a subscript to distinguish them from the corresponding dimensional functions of the dimensional argument $`k`$.

The resulting dimensionless budget equation reads

$`{\displaystyle \frac{\text{d}\epsilon _\kappa }{\text{d}y}}`$ $`+`$ $`{\displaystyle \frac{\epsilon _\kappa }{y}}CT_\kappa `$
$`+`$ $`{\displaystyle \frac{C_1}{e_\text{s}}}\left({\displaystyle \frac{y\epsilon _\kappa ^2}{\rho _\kappa ^2}}\right)^{1/3}\left(1+T_\kappa \right)=𝒲_\kappa ,`$
$`T_\kappa `$ $``$ $`{\displaystyle \frac{\varphi \delta \gamma _\kappa }{\left(1+\varphi \right)\left(1+2\delta \gamma _\kappa \right)+\left(\delta \gamma _\kappa \right)^2}},`$ (98)
$`C`$ $`=`$ $`C_1C_2,\delta =\gamma __L\tau _\text{p}.`$
Here we used Eq. (92) and defined the Reynolds numbers for the carrier fluid, $`e_\text{f}`$, and the effective Reynolds number for the suspension $`e_\text{s}`$

$`e_\text{f}`$ $`=`$ $`{\displaystyle \frac{Lv__L}{\nu }},v__L=\left({\displaystyle \frac{\epsilon __LL}{\rho __L}}\right)^{1/3},`$ (99)
$`Re_\text{s}`$ $`=`$ $`{\displaystyle \frac{Lv__L}{\nu __L}},\nu __L\nu _{\text{eff}}\left(L^1\right)=\nu {\displaystyle \frac{\rho __L}{\rho _\text{f}}},`$
$`\rho __L`$ $`=`$ $`\rho _\text{f}\left[1+\varphi {\displaystyle \frac{1+2\delta }{\left(1+\delta \right)^2}}\right].`$ (100)
in terms of the rms turbulent velocity, $`v__L`$, dominated by $`L`$-eddies. Obviously, $`e_\text{f}`$ involves the kinematic viscosity of the carrier fluid, $`\nu `$, while $`e_\text{s}`$ depends on the effective kinematic viscosity of the suspension, $`\nu _{\text{eff}}\left(L^1\right)`$, for the outer scale of turbulence, $`L`$.

Eq. (IV) has to be considered together with equations for $`\rho _\kappa `$ and $`\gamma _\kappa `$, which follows from Eqs. LABEL:r-eff, LABEL:sum-g and LABEL:gamma2:

$`\rho _\kappa `$ $`=`$ $`\left[1+\varphi {\displaystyle \frac{1+2\delta \gamma _\kappa }{\left(1+\delta \gamma _\kappa \right)^2}}\right]/\left[1+\varphi {\displaystyle \frac{1+2\delta }{\left(1+\delta \right)^2}}\right],`$
$`\gamma _\kappa `$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{C_2e_\text{s}\rho _\kappa }}+{\displaystyle \frac{\epsilon _\kappa ^{1/3}\kappa ^{2/3}}{\rho _\kappa ^{1/3}}}.`$ (101)
These two equations allow us to express the functions $`\gamma _\kappa `$ and $`\rho _\kappa `$ in terms of $`\epsilon _\kappa `$. With these solutions, Eq. (IV) becomes an ordinary differential equation for the only function $`\epsilon _\kappa `$.

The first line of Eq. (IV) describes the effect of particles in the inertial integral of scales. This part involves the mass loading, $`\varphi `$, the dimensionless particle response time (normalized by the eddy life time), $`\delta `$, and the parameter $`C`$, characterizing our version of the K41 closure.

The second line of Eq. (IV) represents the effect of the viscous friction, which is proportional to $`1/e_\text{s}`$, and the pumping term $`𝒲_\kappa `$, which we choose as follows:

$$𝒲_\kappa =\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left[\frac{\left(y1\right)^2}{2\sigma ^2}\right].$$
(102)
This function has a maximum at $`y=1`$ (the input of energy is largest at $`k=1/L`$), while the parameter $`\sigma `$ describes the characteristic width of the pumping region. In addition, the function $`𝒲_\kappa `$ satisfies the normalization constrain

$$\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝒲_\kappa \text{d}y=1,$$
(103)
which follows from Eq. (IV) in the limit $`\sigma 1`$.

## V Solution of the budget equation

### V.1 Simplification of the energy pumping term

First notice, that the turbulence statistics in the energy containing range, $`kL=y1`$ is not universal and depends on the type of energy pumping, in our case, on the function $`𝒲_\kappa `$. Therefore for general analysis, which is not aimed at the study of some particular type of turbulence generation, we can take the pumping of energy in a narrow shell in the $`𝒌`$-space. This means

$$\underset{\sigma \mathrm{\hspace{0.17em}0}}{lim}\left\{𝒲_\kappa \right\}=\delta \left(\kappa \right),$$
(104)
where $`\delta \left(\kappa \right)`$ is the Dirac $`\delta `$ function. In this limit and with zero boundary conditions for $`\epsilon _\kappa `$, $`\gamma _\kappa `$ at $`\kappa =0`$ (and, consequently, $`\rho _\kappa =1`$ at $`\kappa =0`$) , Eq. (IV) can be solved on the interval $`0\kappa 1`$. This gives

$`\epsilon _\kappa =1,\gamma _\kappa =1,\rho _\kappa =1,\text{at}\kappa =1.`$ (105)
In the limit LABEL:lim, Eq. (IV) has zero RHS for $`\kappa >1`$:

$$\frac{\text{d}\epsilon _\kappa }{\text{d}\kappa }+\frac{\epsilon _\kappa }{\kappa }CT_\kappa +\frac{C_1}{e_\text{s}}\left(\frac{\kappa \epsilon _\kappa ^2}{\rho _\kappa ^2}\right)^{1/3}\left(1+T_\kappa \right)=0,$$
(106)
Relations LABEL:bound can be considered as the boundary conditions for Eq. (106) at $`\kappa =1`$.

### V.2 Particle free case and limit of small particles

Consider now the particle-free case, $`\varphi =0`$, and the case of very small particles, $`\delta =0`$, for finite $`e_\text{s}`$. In both cases Eq. (106) becomes:

$`{\displaystyle \frac{\text{d}\epsilon _\kappa }{\text{d}y}}+{\displaystyle \frac{C_1\kappa ^{1/3}\epsilon ^{2/3}}{e_\text{s}}}=0.`$ (107)
We took here in account that according to Eq. (101) $`\rho _\kappa =1`$ for $`\varphi =0`$ and also for $`\delta =0`$. Notice, that for $`\varphi =0`$, $`\nu __L=\nu `$ and, consequently, $`e_\text{s}=e_\text{f}`$, while for $`\delta =0`$, $`\nu __L=\nu /\left(1+\varphi \right)`$ and $`e_\text{s}=e_\text{f}\left(1+\varphi \right)`$. The reason is that for $`\delta 0`$ all particles are fully involved in turbulent motions and one can consider a suspension as a single but heavier fluid with the density $`\rho _\text{f}\left(1+\varphi \right)`$.

The solution of Eq. (107) with the boundary condition $`\epsilon _1=1`$ is

$`\epsilon _\kappa =\left[1+{\displaystyle \frac{C_1}{4e_\text{s}}}\left(1\kappa ^{4/3}\right)\right]^3.`$ (108)
In the particle free case $`\varphi =0`$ and this solution turns into

$`\epsilon _\kappa ^{\left(0\right)}=\left[1+{\displaystyle \frac{C_1}{4e_\text{f}}}\left(1\kappa ^{4/3}\right)\right]^3,`$ (109)
where $`e_\text{s}e_\text{f}`$ as we discussed above. In the bulk of the inertial interval these solutions give a small viscous correction to the K41 solution with the constant energy flux $`\epsilon _\kappa =1`$. Namely

$$\epsilon _\kappa 1+3C_1\left(1\kappa ^{4/3}\right)/\left(4e_\text{s}\right),\text{for}\kappa 1/e_\text{s}.$$
(110)
The local in the $`k`$-space closure procedure, used in the paper, works reasonably well in the inertial interval, where the energy exchange between eddies is dominated by the eddies of compatible scales. However, it is violated in the bulk of the viscous subrange, where the dynamics of eddies of very small scales is dominated not by their self-interaction, but by their energy exchange with $`\eta `$-eddies of the Kolmogorov micro-scale $`\eta `$. Therefore we cannot expect Eq. (106) to reproduce the exponential decay of the energy flux in the viscous subrange. Nevertheless, this equation describes on a qualitative level the behavior of the energy flux until the very end of the inertial interval giving the crossover scale to the viscous subrange, i.e. the value of $`\eta `$. According to Eq. (108), the energy flux become zero at

$$\kappa _{\text{cr}}=\left(1+4e_\text{s}/C_1\right)^{3/4}.$$
(111)
It convenient to introduce here an effective Reynolds number of the carrier fluid and suspensions

$$e_\text{f}^{\text{eff}}=4e_\text{f}/C_1,e_\text{s}^{\text{eff}}=4e_\text{s}/C_1,$$
(112)
which enters in the corresponding Kolmogorov-41 evaluations of the viscous cutoff. For example for the fluid

$$\eta _\text{f}=L/\kappa _{\text{cr}}L/\left[e_\text{f}^{\text{eff}}\right]^{3/4}.$$
(113)
### V.3 Iterative solution in the inertial interval

In the bulk of the inertial interval, after neglecting the viscous terms (i.e. for $`e_\text{s}\mathrm{})`$, Eq. (106) becomes

$`{\displaystyle \frac{\text{d}\epsilon _\kappa }{\text{d}\kappa }}`$ $`+`$ $`{\displaystyle \frac{\epsilon _\kappa }{\kappa }}{\displaystyle \frac{C\varphi \delta \gamma _\kappa }{\left(1+\varphi \right)\left(1+2\delta \gamma _\kappa \right)+\left(\delta \gamma _\kappa \right)^2}}=0,`$ (114)
$`\gamma _\kappa ^3\rho _\kappa `$ $`=`$ $`\epsilon _\kappa \kappa ^2,`$ (115)
$`\rho _\kappa `$ $`=`$ $`\left[1+\varphi {\displaystyle \frac{1+2\delta \gamma _\kappa }{\left(1+\delta \gamma _\kappa \right)^2}}\right]/\left[1+\varphi {\displaystyle \frac{1+2\delta }{\left(1+\delta \right)^2}}\right].`$
#### V.3.1 Large scale solution of the budget equation

In region of large scales, $`\kappa 1`$, the functions $`\rho _\kappa 1`$ and we can simplify Eq. (115) by the replacement $`\rho _\kappa 1`$ in the equation for $`\gamma _\kappa `$, i.e. $`\gamma _\kappa \epsilon _\kappa ^{1/3}\kappa ^{2/3}`$. In the denominator of Eq. (114), where the $`\kappa `$-dependence of $`\gamma _\kappa `$ is less essential, we can make further simplification, replacing $`\gamma _\kappa \kappa ^{2/3}`$. The resulting equation allows separation of variables:

$`{\displaystyle \frac{2}{\delta }}{\displaystyle \frac{\text{d}\epsilon _\kappa ^{1/3}}{\text{d}\kappa ^{2/3}}}`$ $`=`$ $`C\mathrm{\Psi }_0\left(\kappa \right),`$ (116)
$`\mathrm{\Psi }_0\left(\kappa \right)`$ $`=`$ $`{\displaystyle \frac{\varphi }{\left(1+\varphi \right)\left(1+2\delta \kappa ^{2/3}\right)+\delta ^2\kappa ^{4/3}}}.`$
The solution of this equation with the boundary conditions $`\epsilon _1=1`$ is $`\epsilon _\kappa =\epsilon _{0,\kappa }`$, where

$`\epsilon _{0,\kappa }`$ $`=`$ $`{\displaystyle \frac{1}{\left[1+CJ_0\left(\kappa \right)\right]^3}},`$ (117)
$`J_0\left(\kappa \right)`$ $`=`$ $`{\displaystyle \frac{\delta }{2}}{\displaystyle \underset{1}{\overset{\kappa }{}}}\mathrm{\Psi }_0\left(x\right)\text{d}x^{2/3}`$
$`=`$ $`{\displaystyle \frac{\sqrt{\varphi }}{4\sqrt{1+\varphi }}}\{\mathrm{ln}\left[{\displaystyle \frac{\delta \kappa ^{2/3}+1+\varphi \sqrt{\varphi \left(1+\varphi \right)}}{\delta +1+\varphi \sqrt{\varphi \left(1+\varphi \right)}}}\right]`$
$`\mathrm{ln}\left[{\displaystyle \frac{\delta \kappa ^{2/3}+1+\varphi +\sqrt{\varphi \left(1+\varphi \right)}}{\delta +1+\varphi +\sqrt{\varphi \left(1+\varphi \right)}}}\right]\}.`$
Now with Eq. (115) we find the following approximations

$`\rho _{0,\kappa }`$ $`=`$ $`\left[1+\varphi {\displaystyle \frac{1+2\delta \epsilon _{0,\kappa }^{1/3}\kappa ^{2/3}}{\left(1+\delta \epsilon _{0,\kappa }^{1/3}\kappa ^{2/3}\right)^2}}\right]/\left[1+\varphi {\displaystyle \frac{1+2\delta }{\left(1+\delta \right)^2}}\right],`$
$`\gamma _{0,\kappa }`$ $`=`$ $`\left(\epsilon _{0,\kappa }\kappa ^2/\rho _{0,\kappa }\right)^{1/3}.`$ (118)
Formally speaking, the analytical solution (117118) is valid only for $`\kappa 1`$. To find the solution of the initial Eq. (114) in the whole interval of scales, we will iterate this equation, taking the analytical solution as a starting function. Comparing in Sec. V.4 this solution with the next order iterations and with the numerical solution of Eq. (114), we will find an actual region of applicability of the analytical solution (117118).

#### V.3.2 First improvement and subsequent iterations

With the analytical solution (117118), we can improve approximation LABEL:bud1 of Eq. (114) by $`\rho _\kappa \rho _{0,\kappa }`$ \[instead of $`\rho _\kappa 1`$\], which gives $`\gamma _\kappa \epsilon _\kappa ^{1/3}\kappa ^{2/3}/\rho _{0,\epsilon }`$ in the numerator of Eq. (114). In the denominator we replace $`\gamma _\kappa \epsilon _{0,\kappa }^{1/3}\kappa ^{2/3}/\rho _{0,\epsilon }`$. The improved simplification of Eq. (114) reads:

$`{\displaystyle \frac{2}{\delta }}{\displaystyle \frac{\text{d}\epsilon _\kappa ^{1/3}}{\text{d}\kappa ^{2/3}}}`$ $`=`$ $`C\mathrm{\Psi }_1\left(\kappa \right),`$ (119)
$`\mathrm{\Psi }_1\left(\kappa \right)`$ $`=`$ $`{\displaystyle \frac{\varphi }{\rho _{0,\kappa }^{1/3}\left[\left(1+\varphi \right)\left(1+2\delta \gamma _{0,\kappa }\right)+\left(\delta \gamma _{0,\kappa }\right)^2\right]}}.`$
Integration of this equation gives the first iterative solution of Eq. (114), $`\epsilon _\kappa =\epsilon _{1,\kappa }`$, where

$`\epsilon _{1,\kappa }`$ $`=`$ $`{\displaystyle \frac{1}{\left[1+CJ_1\left(\kappa \right)\right]^3}},`$ (120)
$`J_1\left(\kappa \right)`$ $`=`$ $`{\displaystyle \frac{\delta }{2}}{\displaystyle \underset{1}{\overset{\kappa }{}}}\mathrm{\Psi }_1\left(x\right)\text{d}x^{2/3}.`$
This allows further improvement of approximations LABEL:gamma0

$`\rho _{1,\kappa }`$ $`=`$ $`\left[1+\varphi {\displaystyle \frac{1+2\delta \gamma _{0,\kappa }}{\left(1+\delta \gamma _{0,\kappa }\right)^2}}\right]/\left[1+\varphi {\displaystyle \frac{1+2\delta }{\left(1+\delta \right)^2}}\right],`$
$`\gamma _{1,\kappa }`$ $`=`$ $`\left(\epsilon _{1,\kappa }\kappa ^2/\rho _{1,\kappa }\right)^{1/3}.`$ (121)
Now, the next iteration steps are obvious. The $`n`$-order solution is

$`\epsilon _{n,\kappa }`$ $`=`$ $`{\displaystyle \frac{1}{\left[1+CJ_n\left(\kappa \right)\right]^3}},`$ (122)
$`J_n\left(\kappa \right)`$ $`=`$ $`{\displaystyle \frac{\delta }{2}}{\displaystyle \underset{1}{\overset{\kappa }{}}}\mathrm{\Psi }_n\left(x\right)\text{d}x^{2/3},`$
$`\mathrm{\Psi }_n\left(\kappa \right)`$ $`=`$ $`{\displaystyle \frac{\varphi }{\rho _{n1,\kappa }^{1/3}\left[\left(1+\varphi \right)\left(1+2\delta \gamma _{n1,\kappa }\right)+\left(\delta \gamma _{n1,\kappa }\right)^2\right]}},`$
$`\rho _{n,\kappa }`$ $`=`$ $`\left[1+\varphi {\displaystyle \frac{1+2\delta \gamma _{n1,\kappa }}{\left(1+\delta \gamma _{n1,\kappa }\right)^2}}\right]/\left[1+\varphi {\displaystyle \frac{1+2\delta }{\left(1+\delta \right)^2}}\right],`$
$`\gamma _{n,\kappa }`$ $`=`$ $`\left(\epsilon _{n,\kappa }\kappa ^2/\rho _{n,\kappa }\right)^{1/3}.`$ (123)
### V.4 Accuracy of the iterative solutions

To get an understanding of the accuracy of the analytical solution $`\epsilon _{0,\kappa }`$, Eq. (117), and the first iterative solution $`\epsilon _{1,\kappa }`$, we compare them with the ”numerically exact” solutions of Eq. (114), $`\epsilon _\kappa `$ in the wide inertial range of four decades.

We found that for all values of $`\kappa `$ and $`\delta `$ the analytical function $`\epsilon _{0,\kappa }`$ works unexpectedly well for $`C\varphi 0.25`$. To illustrate this, we plot in Fig. 1 functions $`\epsilon _{0,\kappa }`$, $`\epsilon _{1,\kappa }`$ and $`\epsilon _\kappa `$ for $`C=0.25`$ (panel $`a`$) and $`C=1`$ (panel $`b`$) for $`\varphi =0.25,\mathrm{\hspace{0.17em}0.5}`$ and $`\varphi =1`$ with $`\delta =0.1`$. The relative difference between $`\epsilon _{0,\kappa }`$ and $`\epsilon _\kappa `$ is about a few percents for all three cases in panel $`a`$ and for the case $`\varphi =0.25`$ in panel $`b`$.

For larger values of the product $`C\varphi `$ the accuracy of few percents is achieved in the smaller region of $`\kappa `$, where $`\tau _\text{p}\gamma \left(k\right)=\delta \kappa ^{2/3}<1`$, i.e. approximately for $`\kappa \delta ^{3/2}`$. For example for $`\delta =0.01`$ this is three decades, $`\kappa <10^3`$, while for $`\delta =0.1`$ only for $`\kappa <30`$, as we show in Fig. 1$`b`$. Moreover, the first iterative solution, $`\epsilon _{1,\kappa }`$ gives a very good approximation to $`\epsilon _\kappa `$ for all reasonable values of parameters. This is illustrated in Fig. 1$`b`$ for $`C=1`$ and $`\varphi =1`$. Notice, that for $`C=0.25`$ and $`\varphi =1`$ (Fig. 1$`a`$) the plots of $`\epsilon _{1,\kappa }`$ $`\epsilon _\kappa `$ are undistinguished within the line width.

The conclusion is that for the qualitative and semi-quantitative description of the turbulence modification by particles in the inertial interval we can use the analytical solution (117118), corrected, if needed, by the first iteration, $`\epsilon _{1,\kappa }`$.

## VI Turbulence modification by particles

### VI.1 Preliminaries

Consider now separately the density of kinetic energy of the carrier fluid, $`E_\text{f}\left(k\right)`$, and that of the particle, $`E_\text{p}\left(k\right)`$ (i.e. the density of the kinetic energy of the particle velocity field). According to Eqs. LABEL:r-eff, LABEL:Eks3, LABEL:Ek2 and LABEL:energy3

$`E_\text{f}\left(k\right)`$ $`=`$ $`C_1\rho _\text{f}\left[\epsilon \left(k\right)/\rho _{\text{eff}}\left(k\right)\right]^{2/3}k^{5/3},`$ (124)
$`E_\text{p}\left(k\right)`$ $`=`$ $`\varphi {\displaystyle \frac{1+2\tau _\text{p}\gamma \left(k\right)}{\left[1+\tau _\text{p}\gamma \left(k\right)\right]^2}}E_\text{f}\left(k\right).`$ (125)
It is convenient to introduce the dimensionless functions of $`\kappa =kL`$, $`E_\kappa ^\text{f}`$ and $`E_\kappa ^\text{p}`$, both normalized by $`E_\text{f}\left(L^1\right)`$:

$`E_\kappa ^\text{f}={\displaystyle \frac{E_\text{f}\left(k\right)}{E_\text{f}\left(L^1\right)}},E_\kappa ^\text{p}={\displaystyle \frac{E_\text{p}\left(k\right)}{E_\text{f}\left(L^1\right)}},`$ (126)
which may be written as

$`E_\kappa ^\text{f}`$ $`=`$ $`\left(\epsilon _\kappa /\rho _\kappa \right)^{2/3}\kappa ^{5/3},`$ (127)
$`E_\kappa ^\text{p}`$ $`=`$ $`\varphi {\displaystyle \frac{1+2\delta \gamma _\kappa }{\left[1+\delta \gamma _\kappa \right]^2}}E_\kappa ^\text{f}.`$ (128)
Next, introduce the dimensionless ratio

$`R\left(\kappa \right)`$ $``$ $`{\displaystyle \frac{E_\kappa ^\text{f}}{E_\kappa ^{0,\mathrm{f}}}}=\left[{\displaystyle \frac{\epsilon _\kappa }{\epsilon _\kappa ^{\left(0\right)}\rho _\kappa }}\right]^{2/3},`$ (129)
where $`E_\kappa ^{0,\mathrm{f}}=\left[\epsilon _\kappa ^{\left(0\right)}\right]^{2/3}\kappa ^{5/3}`$ is the density of turbulent kinetic energy and $`\epsilon _\kappa ^{\left(0\right)}`$ is the energy flux in the particle-free case, Eq. (109). The ratio $`R\left(\kappa \right)`$ is larger (smaller) than unity in the case of enhancement (suppression) of the turbulent energy by particles.

### VI.2 Energy flux

Our model with local in $`k`$ space parametrization of the energy flux involves the parameter of the closure procedure $`C`$, which has to be considered as a fit parameter which may be evaluated, for example, by comparison with the direct numerical simulation. Generally speaking, it is expected to be of the order of unity. For simplification of the qualitative analysis of the effect of particles on the statistics of turbulence we choose usually $`C=1/4`$, for which we can use the analytical solution (117118). The effect of $`C`$ on the behavior of the function $`\epsilon _\kappa `$ is qualitativly described by the approximation $`\epsilon _\kappa \left(C\right)1+C\left[\epsilon _\kappa \left(C1\right)\right]`$. This is illustrated by comparison of plots $`\epsilon _{0,\kappa }`$ for $`C=0.25`$ and $`C=1`$ in Figs. 1$`a`$ and 1$`b`$.

The plots in Fig. 1 also demonstrate the expected fact that for small $`\varphi `$ the effect of the particles is proportional to $`\varphi `$. This may also clearly be seen from the balance Eq. (114). The balance equation itself also reflects a less trivial fact of saturation of the effect of the particles in the limit $`\varphi 1`$; the beginning of this saturation is clearly seen in Fig. 1. The physical reason for the saturation is that the main governing parameter in the problem is the ratio of the particle energy to the energy of the *suspension* but not to the energy of the *carrier fluid*. As an example of the quantitative description of the effect of saturation consider Eq. (117) for $`\epsilon _{0,\kappa }`$ in the limit $`\kappa \mathrm{}`$:

$$\epsilon _{0,\mathrm{}}=\left[1+\frac{C\sqrt{\varphi }}{4\sqrt{1+\varphi }}\mathrm{ln}\left(\frac{\delta +1+\varphi +\sqrt{\varphi \left(1+\varphi \right)}}{\delta +1+\varphi \sqrt{\varphi \left(1+\varphi \right)}}\right)\right]^3.$$
(130)
Evidently:

$`\epsilon _{0,\mathrm{}}`$ $``$ $`\left[1+{\displaystyle \frac{C\varphi }{2\left(1+\delta \right)}}\right]^3,\text{for}\varphi 1,`$ (131)
$`\epsilon _{0,\mathrm{}}`$ $``$ $`\left[1+{\displaystyle \frac{C}{4}}\mathrm{ln}\left({\displaystyle \frac{2\varphi }{\delta +1}}\right)\right]^3,\text{for}\varphi 1.`$
One sees in Fig. 1 that the increase in the mass loading $`\varphi `$, leads to the suppression of the energy flux for large $`\kappa `$ ( small scales). The onset of this suppression shifts to smaller $`\kappa `$ (larger scales) with increasing $`\delta `$, Fig. 2. To understand this, we note that $`\tau _\text{p}\gamma _\kappa \delta \kappa ^{2/3}`$ is an important governing parameter in the energy budget equation. Consequently, with increasing particle response time, the fluid-particle friction dissipate energy in the larger scale region. As is evident from this figure, the main dissipation of energy occurs in the region $`0.1<\delta \kappa ^{2/3}<10`$. Therefore, for $`\delta >0.1`$ the total loss in the energy decreases with further increase in $`\delta `$.

### VI.3 Suppression and enhancement of turbulence in the inertial interval

As we discussed in Sec. VI.1 the effect of particles on the energy distribution in suspension (with respect to the particle free case) may be characterized by the ratio $`R\left(\kappa \right)`$, Eq. (129). This effect is twofold: the fluid-particle friction leads to suppression of the energy flux with increasing $`\kappa `$. Accordingly, $`\epsilon _\kappa `$ in the numerator of Eq. LABEL:ratio decreases towards larger $`\kappa `$. On the other hand, for larger $`k`$ less particles are involved in the motion and the effective density, $`\rho _\kappa `$, decreases with $`\kappa `$. The factor $`\epsilon _\kappa ^{\left(0\right)}=1`$ in the limit $`e_\text{s}\mathrm{}`$. Therefore in the inertial interval of scales the behavior of $`R\left(\kappa \right)`$ is defined by the strongest $`\kappa `$-dependence either of $`\epsilon _\kappa `$ or of $`\rho _\kappa `$. As we discussed, for $`C<0.25`$ it is sufficient to use the analytical solutions for $`\epsilon _{0,\kappa }`$, Eq. (117), and $`\gamma _{0,\kappa }`$, Eq. (118). This gives:

$`R_0\left(\kappa \right)=\left\{{\displaystyle \frac{\epsilon _{0,\kappa }\left[1+\varphi \left(1+2\delta \right)/\left(1+\delta \right)^2\right]}{1+\varphi \left(1+2\delta \gamma _{0,\kappa }\right)/\left(1+\delta \gamma _{0,\kappa }\right)^2}}\right\}^{2/3}.`$
(132)
Fig. 3$`a`$ demonstrates how the ratio $`R\left(\kappa \right)`$ depends on the fit parameter of our model, $`C`$, which appeared in the budget Eq. (114) in the front of the term, responsible for the fluid-particle friction. Clearly, the relative importance of the fluid-particle friction (with respect to the effect of the density variation) increases with the value of $`C`$. In particular, for $`C=1`$ the friction dominates and $`R\left(\kappa \right)<1`$, for $`C=0.25`$ the density variation dominates and $`R\left(\kappa \right)>1`$. For $`0.25<C<1`$, the density of energy of the carrier fluid is suppressed for smaller $`\kappa `$ and enhanced towards larger $`\kappa `$. As is clearly seen in the figure, the function $`R\left(\kappa \right)`$ has a minimum around some critical value $`\kappa _{\text{cr}}`$ which depends on $`C`$ and $`\kappa `$. For $`C1`$ the value of $`\kappa _{\text{cr}}`$ agrees with the expected estimate $`\tau _\text{p}\gamma \left(k_{\text{cr}}\right)\delta \gamma _{\kappa _{\text{cr}}}1`$. With decreasing $`\delta `$ the position of the crossover ( and of the minimum) is shifted towards larger $`\kappa `$. It is evident that for $`\kappa <\kappa _{\text{cr}}`$ the effect of the fluid-particle friction wins, while for $`\kappa >\kappa _{\text{cr}}`$ the effect of the decrease in the effective density of the suspension is stronger. For $`\kappa \kappa _{\text{cr}}`$ the function $`R\left(\kappa \right)`$ saturates.

Notice that for $`C<0.5`$ the analytical prediction (dashed lines) is pretty close to the ”exact” numerical result (solid lines) indicating the qualitative validity of the analytical description of the effect of particles on the energy distribution in suspensions (within the model limitations). Therefore, one can find the limiting value of $`R_{\mathrm{}}R\left(\kappa \mathrm{}\right)`$ from Eq. (132):

$`R_{0,\mathrm{}}=\left[1+{\displaystyle \frac{\varphi \left(1+2\delta \right)}{\left(1+\delta \right)^2}}\right]^{2/3}`$ (133)
$`\times \left[1+{\displaystyle \frac{C}{2}}\sqrt{{\displaystyle \frac{\varphi }{1+\varphi }}}\mathrm{ln}\left({\displaystyle \frac{\delta +1+\varphi +\sqrt{\varphi \left(1+\varphi \right)}}{\delta +1+\varphi \sqrt{\varphi \left(1+\varphi \right)}}}\right)\right]^2.`$
The analysis of this equation shows that the largest possible enhancement of the turbulent energy in the inertial interval is achieved for $`\delta 1`$ and increases with $`\varphi `$.

### VI.4 Turbulence modification for finite $`e`$

In the previous section we discussed the mechanism of turbulent enhancement in the inertial interval caused by the density variation in the energy cascade processes. There is one more mechanism of the turbulent enhancement, near the viscous subrange, that may be even more important at moderate $`e`$. This effect is due to the renormalization of the kinematic viscosity in suspensions, $`\nu \nu _{\text{eff}}\left(k\right)`$, Eq. (18), caused again by the density variation. Since $`\rho _{\text{eff}}\left(k\right)`$ near the viscous cutoff, $`k1/\eta `$, is larger than $`\rho _\text{f}`$ ( and consequently $`\nu _{\text{eff}}\left(1/\eta \right)<\nu `$), the extent of the inertial interval in suspension is therefore larger than that in the particle free case for the same energy pumping to the system. Within our model this effect may be described for very small particles with a response time smaller than the turnover time of $`\eta `$-eddies. In this case the effective density is $`k`$-independent in the inertial subrange, $`\rho _\kappa =1`$, $`\epsilon _\kappa `$ and $`\epsilon _\kappa ^{\left(0\right)}`$ (for the particle free case) are given by Eqs. LABEL:sol1 and LABEL:sol0. Thus Eq. (129) yields

$$R\left(\kappa \right)=\left[1+\frac{C_1}{4e_\text{s}}\left(1\kappa ^{4/3}\right)\right]^2/\left[1+\frac{C_1}{4e_\text{f}}\left(1\kappa ^{4/3}\right)\right]^2.$$
(134)
Plots of $`R\left(\kappa \right)`$ for different $`e_\text{f}`$ are shown in Fig. 4$`a`$ by dashed lines together with the numerical results for a quite small $`\delta =10^3`$, solid lines. The numerical results for $`\delta =0.01`$ and $`\delta =0.1`$ are shown in Figs. 4$`b`$ and 4$`c`$. With $`e_\text{f}`$ growing above $`10^610^8`$ we return back to the situation in the inertial interval, described above, see Fig. 3. For the comparison we show the plots for $`e_\text{f}\mathrm{}`$ in Fig. 4.

For very small $`\delta `$ the effect of particles on the turbulent statistics in the inertial interval is negligible; as an illustration see Fig. 4$`a`$ for $`\delta =10^3`$. In this case there is only the viscous range enhancement. Clearly, with decreasing $`e_\text{f}`$ this effect is more pronounced. For the moderate values of $`\delta `$, there is a turbulence suppression in the beginning of the inertial interval, which turns into a turbulence enhancement in the bulk of the inertial interval, see, e.g. Fig. 3$`b`$ for $`\delta =0.01`$ and the line marked $`e_\text{f}=\mathrm{}`$ in Fig. 4$`b`$. One sees that already for $`e_\text{f}=10^4`$ the energy enhancement increases. This enlargement becomes more and more pronounced for even smaller $`e_\text{f}`$. For $`e_\text{f}=10^2`$ the turbulence suppression in the beginning of the inertial interval is negligible. Further development of these tendencies is illustrated in Fig. 4$`c`$ for $`\delta =0.1`$

### VI.5 Brief comparison with DNS

In order to get an analytical description of the main physical mechanisms of the particle effect on turbulence we used in this paper as simple as possible approximations, which nevertheless preserve the basic physics of the problem. In particular, we have used the differential approximation of the energy flux term, Eq. (86) with local in $`k`$-space closure procedure, which gives a reasonable approximation in the extended inertial intervals of several decades. However, in the direct numerical simulations of turbulence in suspensions, e.g. in Ref. 98BSS , there is almost no inertial interval, definitely smaller than one decade. Therefore the detailed comparison of our simple theoretical picture with DNS may be only qualitative.

For such a comparison with DNS by Boivin, Simonin and Squires 98BSS ) we re-plotted in Fig. 5 their Fig. 5b for the kinetic energy spectra $`E_\text{s}(k,\varphi )`$ of suspensions in the log-log coordinates (solid lines). The solid line, labelled by $`\varphi =0`$, describes the particle free case, in which the energy spectrum in the inertial interval should be scale invariant. The K41 dependence is shown in Fig. 5 by a dash-dotted line, labelled by $`\kappa ^{5/3}`$. One sees that only the first half of the decade may be considered as the inertial interval. The viscous corrections to this dependence may be accounted for with the help of Eq. (109). Using also Eqs. LABEL:Re-eff and LABEL:Efn, one gets:

$`E_\kappa ^\text{f}=\kappa ^{5/3}\left[1+{\displaystyle \frac{1}{e_\text{f}^{\text{eff}}}}\left(1\kappa ^{4/3}\right)\right]^2.`$ (135)
With an appropriate value of $`e_\text{f}^{\text{eff}}`$ this equation reasonably approximates the numerical data almost in the whole decade of $`\kappa `$, in which $`E_\kappa ^\text{f}`$ decays more than three orders of magnitude, see dashed line $`\varphi =0`$. The chosen value $`e_\text{f}^{\text{eff}}=40`$ agrees with parameters given in Ref. 98BSS  with an acceptable value of the closure parameter $`C_1`$, which enters in the definition (112) for $`e_\text{f}^{\text{eff}}`$. With $`C=13`$ the numerical solutions of Eq. (106) approximate well all the DNS energy spectra $`E_\kappa ^\text{f}\left(\varphi \right)`$ with $`\varphi =0.2,\mathrm{\hspace{0.17em}0.5},`$ and $`1`$ in the region, bounded from above by some value of $`\kappa `$ referred to as $`\kappa _{\text{max}}`$. In this region the spectra decrease from unity (at $`\kappa =1`$) to some values, smaller than $`10^3`$. The value of $`\kappa _{\text{max}}`$ decreases from $`\kappa _{\text{max}}14`$ for $`\varphi =0`$ spectrum to $`\kappa _{\text{max}}7`$ for the spectrum with $`\varphi =1`$.

For $`\kappa >\kappa _{\text{max}}`$ the solutions of Eq. (106) give too small values of the turbulent energy. As already discussed, this is due to the differential approximation for the energy flux, which is absolutely not realistic in the viscous subrange. Clearly, the larger the value of $`\varphi `$, the more energy is dissipated by the fluid-particle friction, diminishing the energy flux at the end of the inertial interval. Consequently, the Kolmogorov microscale $`\eta =\nu ^{3/4}/\epsilon ^{1/3}`$ increases. Since $`\kappa _{\text{max}}1/\eta `$, it shifts toward smaller values.

One observes also some deviations of the DNS data and our numerical solutions in the energy containing region $`\kappa 1`$. This is again related to the differential approximation for the energy flux. To improve the description of the particle effect on the turbulent statistics a better approximation for the energy transfer term is required. This calls for the more elaborated closure procedures, based on a proper analysis of the triad interactions.

## VII Summary

• In this paper we propose a one-fluid dynamical model of turbulently flowing dilute suspensions which differs from the usual Navier-Stokes equation in two aspects:

1. Instead of fluid density, $`\rho _\text{f}`$, our model involves a $`k`$ dependent effective density $`\rho _{\text{eff}}\left(k\right)`$ which varies between $`\rho _\text{f}`$ (for large $`k`$) and the mean density of suspension $`\rho _\text{s}=\rho _\text{f}\left(1+\varphi \right)`$ (for small $`k`$);
2. The model equation includes an additional damping term $`\gamma _\text{p}\left(k\right)`$ which describes the fluid-particle viscous friction.

• Our model may be considered as a mean-field approximation in which one uses a dynamical equation of motion with “effective” coefficients which depend on the statistics of the resulting stochastic solutions. In our case $`\rho _{\text{eff}}\left(k\right)`$ and $`\gamma _\text{p}\left(k\right)`$ are determined by the eddy turnover frequency $`\gamma \left(k\right)`$ which, in its turn, depends on the resulting energy distribution in the system.

• Our model is based on the same set of assumptions (applicability of the Stokes law for the fluid-particle friction and space homogeneity of the particle distribution) as widely used in two-fluid models for suspensions. We believe that the one-fluid description of turbulent suspensions contains the same physics as the essentially more complicated two-fluid models. Our feeling is that a possible minor difference in the level of accuracy between these two models is beyond a current level of understanding of the problem and is definitely smaller than the “absolute” accuracy of each model itself.

• In order to keep the description of the problem as simple and transparent as possible we used in this paper a closure procedure based on the Kolmogorov-41 dimensional reasoning with an additional simplification — the differential form of the energy transfer term in which the energy flux $`\epsilon \left(k\right)`$ is evaluated locally in $`k`$-space, via the spectrum $`\left(k\right)`$ taken at the same wave-number $`k`$. This allows us to derive the quite simple ordinary differential equation for the energy budget in the system LABEL:bud3.

• As a reward, our budget Eq. (106) allows an effective analytical analysis in various important limiting cases, i.e.:

1. In the particle free case, see Eq. (109);
2. For the micro-particles case ($`\delta <e^{3/4}`$), see Eq. (108);
3. For the first decades of the inertial interval (in the case $`\delta <1`$) or in the whole inertial interval (if $`C<1/4`$), see Eqs. LABEL:ve0 and LABEL:gamma0;
4. For any reasonable values of parameters at hand, see Eqs. LABEL:ve1 and LABEL:gamma1, involving one-dimensional integration.

In the general case the budget equation LABEL:bud3 may be easily solved numerically.

• We derived the analytical expression LABEL:Rkappa for the dimensionless ratio $`R_0\left(kL\right)`$, which describes the energy suppression and enhancement in the *inertial interval of scales*.

• In Sec. VI.4 we described the additional ”viscous” mechanism of the turbulence suppression and enhancement, caused by the particle effect on the *extent of the inertial interval*:

1. The decrease of the effective kinematic viscosity in suspensions (due to the increase in the effective density for small scale motions)  *elongate the inertial interval*.
2. The fluid-particle friction causes a decrease of the energy flux at the viscous end of the inertial interval and hence  *shorten the inertial interval*.

The winner of this competition depends mainly on the value of $`\delta `$, see, e.g. Fig. 4.

The complicated interplay of the inertial-range and the viscous-range mechanisms of the suppression and the enhancement of the turbulent activity in suspensions is the main topic of Sects. V and VI.

• Our model successfully correlates observed features of numerical simulations. These features are the following:

1. For a suspension with particles with a response time much larger than the Kolmogorov time the main effect of the particles is suppression of the turbulence energy of fluid eddies of all sizes (at the same energy input as for the particle-free case). See for instance Fig. 5, where a comparison with the DNS-results of Boivin, Simonin and Squires 98BSS  is shown.
2. For a suspension with particles with a response time comparable to or smaller than the Kolmogorov time the Kolmogorov length scale of the fluid eddies will decrease and the turbulence energy of eddies of (nearly) all sizes increases (at the same energy input as for the particle-free case). This result was also reported by Druzhinin 01Dru , who carried out DNS-simulations for the case of micro-particles.
3. For a suspension with particles with a response time in between the two limiting cases mentioned above the energy of the larger fluid eddies is suppressed whereas the energy of the smaller eddies is enhanced. The cross-over between suppression and enhancement depends on the ratio of the particle response time and the Kolmogorov time. The strength of the effect depends on the mass loading. This is in agreement with the DNS-results of Sundaram and Collins 99SC .

The more detailed comparison of our approach to turbulent suspensions with the physical and numerical experiments requires:

1. From DNS side more detailed analysis of joint statistics of the velocity field of the particle and the carried fluid;
2. From the theoretical side an application of the more advanced non-local closure procedures, explicitly accounting for the triad interactions.

• An additional advantage of our one-fluid approach is that one can use standard and well developed closures from analytical theory of one-phase turbulence. This fact and the relative simplicity and physical transparency of the one-fluid model equations may essentially help in the further progress towards a theory of turbulent suspensions for more realistic cases with space inhomogeneities, gravitational settling, etc.

### Acknowledgements

We thank T. Elperin, N. Kleeorin, and I. Rogachevskii for helpful discussions, which contributed to this paper. This work has been partially supported by the Israel Science Foundation and the Nitherlands Foundation of Applied Sciences. Two of us (V.L and A.P) acknowledge the hospitality at the Burgerscentrum, Delft, The Netherlands.