Inertial effects on reactive particles advected by turbulence

We study the problem of the advection of passive particles with inertia in a two-dimensional, synthetic, and stationary turbulent flow. The asymptotic analytical result and numerical simulations show the importance of inertial bias in collecting the particles preferentially in certain regions of the flow, depending on their density relative to that of the flow. We also study how these aggregates are affected when a simple chemical reaction mechanism is introduced through a Eulerian scheme. We find that inertia can be responsible for maintaining a stationary concentration pattern even under nonfavorable reactive conditions or destroying it under favorable ones.


I. INTRODUCTION
Advection of particles by turbulent flows can be considered as one of the most interesting problems in fluid mechanics, and a basic issue in environmental and engineering fields.The mixing efficiency of such flows has been an object of study for years and the idea that turbulence is a good mixing mechanism is today inconsistent with experimental observations that show the existence of particle aggregation patterns even in the most turbulent flows.This is related to the fact that turbulent flows are clearly characterized by structure and underlying coherence ͓1͔.Parallel to this finding, numerical simulations of the dispersion of passively advected particles in two-dimensional ͑2D͒ kinematic turbulence ͓2,3͔ show a clustering feature proving the poor mixing ability of such flows.
This aggregation phenomenon is even more dramatic when considering nonpassively advected particles, i.e., taking into account inertial effects.In this case, particles do not follow the flow lines of the fluid exactly, and therefore, even under an incompressible turbulent flow, an initially homogenized condition segregates and particles tend to accumulate in certain regions of the advecting field.The existence and characterization of such regions ͓4͔, as well as the study of the motion of particles in a turbulent fluid, are keys to the understanding of a large variety of problems like the settling of aerosol particles ͓5,6͔, dispersion of pollutants in the atmosphere, distribution of planktonic organisms in the ocean ͓7͔, etc.
Many numerical and analytical approaches have been devoted to studying the motion of particles in turbulent flow fields.It is well known that the relation between the density of the fluid and that of the advected particles largely determines the zones of the flow where the particles accumulate.In this paper, we reconsider this problem using a synthetic, two-dimensional turbulent flow with some prescribed statistical properties.In Sec.II we present the flow and a brief introduction to the procedure through which we generate it.In Sec.III we present the Lagrangian scheme to study the segregation of heavy and light particles under turbulent flows, also including molecular diffusion.As well as checking the validity of our synthetic turbulent field in reproducing the already known behavior of inertial particles advected by turbulence, we also study the effects of inertia when the advected particles are able to react autocatalytically in some way.To study this problem we first move from a Lagrangian scheme to a Eulerian approach.The details, approximations, and numerical procedure of such an approach are described in Sec.IV.Once the Eulerian scheme is implemented we include a reaction term with creation and annihilation contributions, and we proceed to numerically study it for different particular situations in Sec.V.In particular, we show that the inertial drift can permanently sustain a nonzero mean concentration pattern under nonfavorable reaction-diffusion conditions for heavy particles and, contrarily, can wash out a system of light particles even under favorable reactiondiffusion conditions.

II. SYNTHETIC VELOCITY FIELD
We work with a statistically homogeneous, isotropic, and stationary two-dimensional velocity field which represents a ''synthetic'' or kinematic turbulent flow ͓8͔.This flow corresponds to the velocity of the elements of fluid where the particles are immersed, and it is not supposed to be affected by them.However, the presence of inertia makes the velocity of these particles differ from the velocity of the fluid at the particle position.In order to distinguish between the two velocities, we denote the fluid velocity by U i, j ͓since its generation procedure is discretized in space using a square lattice (i, j)], whereas when referring to the particles velocity we will use V i, j or V(X n ) ͑depending on whether we are using a Eulerian or a Lagrangian scheme, respectively͒.
The generation of the turbulent flow comes from the twodimensional simulation of a Langevin equation for the stream function (r,t), where 0 stands for the kinematic viscosity and Q͓ 2 ٌ 2 ͔ denotes an operator that controls the spatial structure of the flow, with its typical correlation length.Furthermore, (r,t) represents a Gaussian white-noise field with zero mean value and whose covariance is given by where ⑀ 0 is the parameter that determines the intensity of the noise and further on that of the mimicked turbulent flow.This Langevin equation can be formally integrated in Fourier space to get the temporal evolution of the stream function.
Our technique corresponds to building up the flow field from its independent Fourier modes, and this implies changing from an intrinsic randomness ͑associated with the complex behavior resulting from the nonlinearity of the Navier-Stokes equation͒ to a system of independent Fourier modes coupled to an external spatiotemporal noise with prescribed statistical characteristics.Therefore, the incompressible twodimensional velocity field U i, j follows from the stream function, This synthetic fluid flow is characterized by three basic well-defined statistical properties: u 0 2 , the intensity of the flow, l 0 , the length correlation of the flow, and t 0 , the time correlation of the flow ͓8͔.These flow parameters are expressed in terms of the velocity correlation function, which in turn is formally associated with the energy spectrum ͑which depends on the form of the operator Q), and depends explicitly on the input parameters 0 , ⑀ 0 , and .In this paper we consider a flow with the Kraichnan spectrum ͓9͔.A detailed presentation of the algorithm just introduced to simulate turbulent flows can be found in ͓8͔.
For the sake of simplicity, and in order to avoid the existence of a great number of temporal scales, we take a frozen flow, so that we have fixed t 0 ϭϱ in our flow simulations.The whole procedure is discretized in space using a square lattice of NϫN points and unit spacing ⌬r.In all the simulations in this paper, Nϭ128 and ⌬rϭ0.5.

III. LAGRANGIAN SCHEME
A visual and suitable technique to study particle mixing is by means of the Lagrangian scheme.This approach has been used in most of the literature on this topic.In this section we present the most relevant aspects of the problem of inertial particles advected by turbulence using the proposed synthetic turbulent field.

A. Model and equations
According to a Lagrangian description of the system, we follow the dynamics of a given set of particles by solving where the subscript n denotes the particle label, X n stands for its position, and V n its velocity.
Our starting point in order to determine the influence of the fluid velocity on the particle trajectory is the complete version of the equation of motion for a rigid sphere in a flow proposed in ͓10͔.This equation reads where m p is the mass of the particle, a corresponds to its diameter, m f is the mass of the fluid in the spherical volume displaced by the advected particle, and and are the dynamic and kinematic viscosities, respectively.U(X n ) is the flow velocity at the point where the particle n is located.
Since our fluid flow is generated using a discrete lattice U i, j , we have to interpolate these values using a bilinear form in order to get U(X n ) ͓3͔.
Notice that we distinguish between two different time derivatives: d/dt is used to denote a time derivative following the particle motion, whereas D/Dt corresponds to a time derivative following a fluid element, Equation ͑5͒ is derived under the assumption that the particle radius, its Reynolds number, and the velocity gradients around it are very small ͓10͔.It has also been assumed that the particles are sufficiently small so that Faxen corrections ͓11͔ have been neglected.The first term on the right is the Bernouilli term, and represents the force acting on the particle from the undisturbed flow.The second contribution stands for the viscous Stokes drag.The third term corresponds to the added mass contribution, and the last one to the Basset history force ͓12͔.In this paper we consider the cases where the drag force dominates over other forces, so that the Basset term can be neglected.
According to the intuitive idea that the added-mass term stems from the force that the particle has to exert in order to move the fluid out from its new position and since this effort depends locally on the state of the fluid at that position, we decide to change the ''particle derivative'' to the ''fluid derivative'' for the flow field in this term.This change was first proposed by Taylor ͓13͔, noticed by Auton et al. ͓14͔, and assumed as correct in much further work ͓15͔.Including all the former approximations and changes, Eq. ͑5͒ now reads

͑7͒
Then if we substitute Eq. ͑6͒ into Eq.͑7͒, taking into account that our flow field does not change in time, and we rescale space, time, and velocity by the scale factors l 0 , l 0 /u 0 , and u 0 , we obtain the final equation for the velocity of a particle, For the sake of simplicity we keep the same notation after the adimensionalization.We note here that all quantities plotted in the figures of this paper are dimensionless.A hereafter refers to the dimensionless inertia parameter defined as the ratio between the viscous and inertial forces, and ␤ is the relation between masses that determines the importance of the Bernouilli term, We are interested in the cases of advected particles for which Eq. ͑8͒ applies and the Stokes term dominates (A Ͼ1).In this case, as we have already anticipated, the density of the particles relative to that of the fluid is the relevant quantity that fixes the magnitude of the Bernouilli term and, in turn, controls the preferred directional motion of particles.Although other intermediate cases have been studied ͓12,15͔, we will focus here on the cases of either extremely heavy or extremely light advected particles.These two cases will show a very different aggregation behavior under turbulence.

B. Asymptotic approach for an effective particle velocity field
Equation ͑8͒ cannot be solved analytically but it is possible to get reliable information following an asymptotic analysis introduced by Maxey ͓5͔.In this reference, the Bernouilli term was not considered since only very heavy particles were studied.Here we include this term in order to study also the case of particles lighter than the fluid.We start with the formal integration of the particle motion following Eqs.͑4͒ and ͑8͒.As A is large, one eliminates exponential transients and by differentiation the following expression is obtained for the velocity of the particle: where the last term will be neglected because we consider only the first order in A Ϫ1 .
As the time derivative of U along the trajectory X(t) can be expressed as it can also be expanded in powers of 1/A, giving to first order Considering that we are using a frozen flow ‫ץ‬U/‫ץ‬tϭ0, the final expression for the particle velocity in Eq. ͑11͒ reads This is an important result of this paper and according to it the particle velocity is completely determined only by its position, so that we can define an effective particle velocity field v(r) as Before going ahead we have to realize the physical sense of the approximations considered to obtain Eq. ͑15͒.By eliminating the initial transients and taking only the first term in the expansion of dU/dt, one assumes that the particles rapidly adopt the velocity of the fluid around, since they arrive at their new position with a velocity not very different from the fluid velocity at that point.Under this condition, in some sense, we ''forget'' about the particle trajectory and therefore an effective particle velocity field v(r) can be defined, and, more important, it gives relevant information about the real motion of the particles.
As corroborated when performing numerical simulations, one of the most relevant effects of particle inertia consists in a drift of the particle motion away from the advecting flow lines and, as a consequence, the aggregation of particles in certain regions of the flow field.This important effect is also captured by the asymptotic approach developed above.Particles will accumulate in those zones of the effective particle velocity field with negative divergence, "•vϽ0.Notice that, since our turbulent flow is incompressible, "•Uϭ0, passively advected particles do not get collected in any part of the system.However, the second term of Eq. ͑15͒ gives a contribution that makes the particles follow a different path from that of the flow field lines.Starting from Eq. ͑15͒, the divergence of v"r… can be expressed as a function of some local properties of the original turbulent flow U,

͑16͒
where S 2 stands for the squared local strain rate of the flow U and ͉͉ 2 corresponds to the squared modulus of its local vorticity.
From Eq. ͑16͒ it is very easy to distinguish two different cases.On one hand, heavy particles, ␤Ͻ1, will accumulate in regions of the turbulence with S 2 Ͼ 2 /4, i.e., low vorticity and high strain rate ͓4 -7͔.On the other hand, light particles, ␤Ͼ1, collect in those zones of the turbulent flow with S 2 Ͻ 2 /4, i.e., high vorticity and low strain rate.These effects can also be conjectured from the observation of Fig. 1 where the particle velocity field v"r… in Eq. ͑15͒ is plotted for generic heavy and light particles.From this figure one easily realizes ͑look at the arrows of the velocity field͒ that inertia makes the heavy particles escape from inside the eddies until they get into a path of almost parallel flow lines.In contrast, light particles seem to be attracted by the eddies of the turbulent flow.These results coming from the asymptotic approach will be confirmed by a numerical study of particle distributions presented in Sec.III C.

C. Numerical results
In order to get numerical results for the motion of particles under turbulence, we integrate Eqs.͑4͒ and ͑8͒ using a second order Runge-Kutta method with a time step ⌬t.For the Lagrangian results in this paper we use a value of ⌬t ϭ0.005.We choose a spatially random initial distribution of particles all over the system and the initial velocity of the particles has been fixed to the velocity of the fluid at their original positions.Our Lagrangian simulations contain typically 1000 particles.Periodic boundary conditions apply for the velocity flow as well as for the motion of particles.
In most of the simulations with turbulence in this paper, the frozen flow in the second panel of Fig. 1 is used, with u 0 2 ϭ1, l 0 ϭ6, and an inertia parameter Aϭ6.This value for A has been checked as the optimum in order to observe strong effects of inertia in the motion of the advected particles.Much smaller values would imply that the particles respond very slowly to the changes in the flow during its motion.On the other hand, for much larger A, the particles rapidly adopt the velocity of the flow and behave as a fluid element.Under both these extreme situations, particles show no aggregation behavior.
The behavior of advected particles under these conditions is shown in Figs. 2 and 3, for extremely heavy particles (␤ ϭ0) and very light particles (␤ϭ2), respectively.Comparing the two figures, one notices how the aggregation behavior is evident in both cases, but the regions where the particles aggregate are quite different.For the case of denser particles, we observe a strong aggregation of moving particles along very narrow paths.In contrast, the case of lighter particles shows that they accumulate at rest inside the vortical structures of the flow.
In Fig. 4 we observe the transition from the pattern of circulation paths to the clustered pattern corresponding to light particles.In this figure, snapshots at tϭ180 are plotted for systems with Aϭ6 but different values of ␤.
At this point it is worth defining a variable which could give us quantitative information concerning the degree of particle aggregation.We denote this aggregation variable as ⌸ and we define it through In order to compute this variable we have divided our system into n cells square cells with linear size L c ͑in our results we use n cells ϭ32 2 ).N(k,t) stands for the number of particles inside the cell k at time t, and N part is the total number of particles.According to this definition, ⌸ gives the fraction of the system occupied.In particular, the limits are 1 for a completely homogenized system and 1/n cells for a system with all the particles collected in one cell.In our case, the initial configuration of the system follows a Poissonian distribution, and therefore the value for ⌸(0) will be smaller than 1.In our initial Poissonian distributions of 1000 particles it can easily be checked that ⌸(0)ϭ0.5.
We compute the ⌸ function for several values of ␤, and average over 100 realizations of the flow and initial particle distribution for each case.The results are shown in Fig. 5 and we can see how for extreme values of ␤ the particles show the largest tendency for preferential accumulation.In the case of ␤ϭ1, the system is rather homogenized, and this result is coherent with Eq. ͑16͒ of our asymptotic approach.We have to point out, however, that cases with ␤ near 1 also show aggregation but after long times ͓of the order of (1 Ϫm p /m f ) Ϫ1 ].In these situations, the Basset force and Faxen corrections, which have been neglected here, play an important role in the particle motion drift.A complete analysis of this situation can be found in ͓12͔.
What is also important to notice from Fig. 5 is that aggregation is quantitatively stronger, in general, for light particles than for heavy ones.This is caused by the fact that light advected particles are attracted to ''points'' of the flow ͑centers of vortical structures͒, whereas heavy ones circulate following paths across the whole system size, avoiding flow eddies.
The validity of the asymptotic approach in Eq. ͑16͒ is also quantitatively checked in Fig. 6, where the Lagrangian average ͑particle average͒ for the S 2 and ͉͉ 2 variables has been plotted for simulations of heavy and light particles advected by turbulence.We compare their behavior with a noninertial case.We can clearly see how the inertial heavy particles accumulate in regions of higher strain rate and lower vorticity than passively advected ones, whereas the light ones tend to aggregate in regions of lower strain rate and higher vorticity than the noninertial ones.

D. Diffusion in the Lagrangian scheme
Before moving to the Eulerian scheme, and in order to compare the Lagrangian results with the Eulerian ones that will be presented in the next section, we introduce molecular diffusion in this discrete scheme by adding a noise term in Eq. ͑4͒, where n corresponds to a zero-mean, Gaussian white noise, and D is the diffusion coefficient.In Fig. 7 the corresponding concentration maps at tϭ180 for the Lagrangian simulations of heavy and light particles are shown for Dϭ0.25.The effect of diffusion consists in making the aggregation zones wider when comparing panels in Fig. 7 with the particle plots of Figs. 2 and 3.The only difference between heavy and light particles is that, in the case of heavy particles, a larger diffusion makes possible the opening of new paths for circulating particles that were not accessible in the simple deterministic simulations ͑Fig.4͒.On the other hand, for light particles the final patterns are qualitatively the same but simply with more diffused clusters.This difference is important when introducing reaction in these systems.

IV. MOVING TO THE EULERIAN SCHEME
Up to this point the typical Lagrangian approach for inertial particles advected by a flow has been used.Nevertheless, our main interest in what follows will consist in incorporating reaction effects affecting the dispersed inertial particles.In order to include reaction we decided to change the particle description and move to a continuum scheme of particle concentrations.The main reasons for using this approach are the following.First, the inclusion of a reaction term is rather easy to implement as a part of an advection-reaction equation, simplifying the numerical and analytical procedure.Second, there is a clear computational advantage: field approaches usually need much less computational effort to obtain the same results than schemes based on particle dynamics, because smaller systems and shorter CPU times are needed.On the other hand we have to be aware of the limitations of this technique, since some approximations are implicitly assumed when working with such mean field ͑meso-scopic͒ approaches.

͑19͒
where c(r,t) corresponds to the local concentration, D is the diffusion constant, R(r,t) stands for the reaction term, and v(r,t) is the velocity of the particles at the position r at time t.This velocity is related to the particle velocity field as given in Eq. ͑8͒, dv͑r,t ͒ dt ϭA͓U͑r,t ͒Ϫv͑ r,t ͔͒ϩ␤͓U͑ r,t ͒•"͔U͑ r,t ͒.

͑20͒
In general, the total time derivative of the variables associated with the constituent particles of the system, and considering that matter flows with a velocity v, can be decomposed as usual, From Eqs. ͑20͒ and ͑21͒ one obtains the temporal evolution for the local particle velocity, ‫ץ‬v͑r,t ͒ ‫ץ‬t ϭA͓U͑r,t ͒Ϫv͑ r,t ͔͒Ϫ͓v͑ r,t ͒•"͔v͑ r,t ͒ ϩ␤͓U͑r,t ͒•"͔U͑ r,t ͒. ͑22͒ Both Eqs.͑19͒ and ͑22͒ will determine the evolution of our system.
Before going ahead we want to compare the Lagrangian and Eulerian prescriptions.The Eulerian scheme in Eq. ͑22͒ generates v as a function of the position and time.In any case, since we work here with a frozen fluid flow U(r), the evolution of v will reach a stationary state that makes this velocity depend only on r.At this point it is important to remark, however, that the Eulerian prescription implicitly assumes a very strong condition.In particular, the assumption of the very existence of a particle velocity field v implies that all the particles that reach the same position at the same time have the same velocity, no matter where they come from.The same assumption is implicitly made in the analytical approach in Sec.III B when obtaining the effective velocity particle field within a Lagrangian scheme, Eq. ͑15͒.Therefore, both approximations stem from similar assumptions, and we have checked numerically that the velocity field v obtained for the asymptotic Eulerian treatment in Eq. ͑22͒ and the one following from the asymptotic approximation in the Lagrangian scheme, Eq. ͑15͒, are rather similar.In other words, the Eulerian scheme should be considered as appropriate only when referring to large values of the inertia parameter A.
A technically difficult step at this point is to discretize Eqs.͑19͒ and ͑22͒.In general, the discretization of an advection term is rather complicated when trying to obtain good numerical accuracy.Specifically, this problem is even worse in systems leading to aggregation since in these cases we have to deal with high concentration gradients.We have chosen the two-step Lax-Wendroff scheme, which is a second order in time method that defines the intermediate values of concentrations at half time steps at the half mesh points.A complete description of this method in one dimension can be found in ͓16͔.This discretization scheme has been extensively checked and compared with other explicit methods, FIG. 7. Concentration map at tϭ180 for a diffusion coefficient Dϭ0.25 and Aϭ6.Left: heavy particles (␤ϭ0).Right: light particles (␤ϭ2).In these Lagrangian simulations 10 6 particles are advected and the concentration values have been obtained by making a cross-grained average.giving us the best numerical accuracy.The systems are spatially discretized using a square lattice of 128ϫ128 cells and unit spacing ⌬rϭ0.5.In all the Eulerian simulations ⌬t ϭ0.0002.
Before going ahead, we first reproduce the results of former sections for systems without reaction.Actually, although the Lax-Wendroff scheme gives good numerical accuracy, a diffusion term is needed in order to assure numerical convergence, at least when no reaction is considered in Eq. ͑19͒.
We start from a random initial distribution of concentrations by assigning a local concentration c i, j (tϭ0) around 1. Depending on the value of D the behavior of the system changes slightly.Figure 8 shows the concentration distributions for heavy particles ͑left͒ and light particles ͑right͒ for a given diffusion constant Dϭ0.25.Gray scales have been used here: black for the higher concentrations and white for the smaller.The results obtained from this scheme are very similar to the ones from the Lagrangian simulations with diffusion; compare Fig. 8 with Fig. 7.For small D the concentration distributions are rather similar to those obtained in the simple deterministic Lagrangian simulations, and, similarly to the Lagrangian simulations with diffusion, larger D makes the aggregation zones get wider.

V. DIFFUSION ¿ ADVECTION ¿ REACTION PROBLEM
In this section we consider the whole problem of Eq. ͑19͒; namely, a system with reactive particles advected by a quenched kinematic turbulent flow and subject also to diffusion.For the sake of simplicity, and as a formal example, we take a reaction term that reads R͑r,t ͒ϭKc͑ r,t ͓͒1Ϫc͑ r,t ͔͓͒c͑ r,t ͒Ϫ0.5͔, ͑23͒ where K is the reaction constant.This reaction term has three steady states, two of them stable attractors at concentrations 0 and 1.As a consequence of these bounds, the reaction avoids the large gradients of particle concentrations that appear in the nonreactive case.This fact allows us to obtain rather good numerical accuracy even with small or zero diffusion coefficient.
The effect of the reaction term in the distribution of concentrations consists of eliminating matter in those zones where cϽ0.5 ͑first attractor at cϭ0) and in building up and maintaining a concentration close to 1 in those zones where cϾ0.5 ͑second attractor at cϭ1).Obviously, inertial drift is going to accumulate particles in particular zones of the system, and in those zones the second attractor will dominate transiently.The other, nonpreferred, zones will fall into the first attractor and become empty.For this reason, patterns of concentration distributions will appear and will be similar, at least at intermediate times, to those appearing in the nonreactive case.Besides this general behavior, the different roles of diffusion, reaction constant, and the initial concentration are going to lead to some different dynamical effects in the temporal evolution of such systems.
In order to study the effects of diffusion and reaction we ran a set of simulations with initial local concentrations randomly distributed around 0.6 and different values of D and K.As expected, with an initial mean concentration larger than 0.5, in the absence of inertia, the system completely fills up and falls into the attractor at cϭ1.This is due to the fact that the zones with local concentration larger than 0.5 finally dominate whatever the values of K and D. However, when inertia is included, some of the regions of the system become empty, leading to new stationary states with smaller mean concentrations ͑in some cases quite close to zero values͒ than in the case of passively advected reactive particles without inertia.Diffusion and reaction will affect the temporal evolution of such aggregation patterns and, in turn, the newly appearing stationary states.In Fig. 9 the temporal evolution of the averaged concentration is presented for different values of D and K, for heavy particles ͑left͒ and light particles ͑right͒.Obviously, for a given K, the larger D is the more zones of the system completely fall into the first attractor (cϭ0), leading to smaller stationary concentrations.This can be checked by looking at the curves with Kϭ1 in both panels of Fig. 9. On the other hand, competition between diffusion D and reaction K is exhibited, and a larger mean concentration value can be maintained for a given D by fixing a larger value for the reaction constant K. Comparing both panels, one realizes that for light particles smaller stationary mean concentrations are achieved as compared to the situation for heavy particles, and some systems even get completely empty despite the initial concentration being larger than 0.5.This result is related to the fact that light particles accumulate in smaller regions of the system than heavy particles, as we pointed out in former sections.
The opposite effect of inertia, at least for heavy particles, is achieved when the initial conditions correspond to an initial mean concentration below 0.5.These cases are shown in Fig. 10 where simulations for different diffusion coefficients and reaction constants are plotted, again for heavy and light particles in the left and right panels, respectively.We observe in both cases that with such small values for the local concentrations, in the absence of inertia, the system falls into the first attractor and becomes totally empty ͑whatever the values of D and K).However, for the case of heavy particles, the presence of inertia makes the particles collect in some places where the concentration threshold at cϭ0.5 can be exceeded and consequently a segregation pattern can remain for ever ͑at least for not very large values of D).This idea is enforced by looking at the mean concentration curves plotted in the left panel of Fig. 10 for Kϭ1 and different values of D. Again, a large value of D can be compensated by fixing a large enough value for K ͑for instance, in the case of heavy particles, a higher mean concentration for Dϭ0.01 is attained by increasing K from 1 to 3).This effect was not obtained for any set of values (K,D) for light particles.In these latter cases particles are always attracted to such small regions of the system as can be easily annihilated by reaction if we start from a mean concentration below 0.5.
Finally, in Fig. 11 cases for heavy particles are shown with different initial mean concentrations c 0 but for the same inertia, flow, diffusion, and reaction parameter values.As we have already pointed out, in the absence of inertia, the cases with c 0 Ͻ0.5 fall into the first attractor and the system becomes totally empty, whereas the cases with c 0 Ͼ0.5 fall into the second attractor, corresponding to a local concentration equal to 1 in the whole system.However, the existence of inertia makes these heavy particles achieve intermediate states of a stationary concentration.The larger the initial mean concentration, the more mass is maintained at long time, but always smaller than cϭ1.On the other hand, when starting from a rather small c 0 the system cannot sustain even a very small concentration.In an interval of initial concentration around c 0 ϭ0.5 the final state does not depend on this value.
The main effect we want to stress from all the results presented in this section is that inertia can maintain a mean concentration of heavy particles in our system despite the fact that the reaction and diffusion conditions are nonfavorable and try to make the system empty.On the contrary, when the conditions are favorable in the absence of inertia, the inclusion of the Stokes force reduces and even makes vanish the final stationary concentration.For light particles, however, inertial drift always has a destructive effect when combined with reaction.For the most favorable cases, c 0 Ͼ0.5, a very low ͑near zero͒ stationary mean concentration is reached.In contrast to the cases with heavy particles, any nonzero final concentration can be reached for light particles when starting from c 0 Ͻ0.5.This is easily explained because heavy particles accumulate in paths along the system, but light particles aggregate in a few spots at the center of the flow eddies.So it is clearly shown how inertia will determine  the long-time behavior of reactive particles advected by turbulence, depending on their relative density.

VI. CONCLUSIONS
In this paper the problem of inertial particles advected by turbulence has been studied.First, we applied the Lagrangian prescription to our model to characterize the presence of an aggregation behavior for heavy and light particles.We also extended the asymptotic approach developed by Maxey ͓5͔ to the cases where the Bernouilli term has to be considered ͑light particles͒ and we checked the validity of this approximation for light and heavy particles with our numerical simulations.We analyzed how regions of the flow where particles accumulate are related to flow properties.Heavy particles escape from the flow eddies and circulate following thin paths between them.In contrast, light particles are attracted by the regions of turbulence with high vorticity and remain accumulated inside the vortical structures of the flow.
In the second part of the paper we extended our model of particles with inertia to a Eulerian prescription.We found good agreement between the results of both Lagrangian and Eulerian schemes, and we studied under what conditions the continuous method is valid.As an application of the Eulerian scheme, we considered chemical reaction effects between the advected particles.We examined the importance of diffusion and inertia in these systems and found that inertia is responsible for maintaining a stationary concentration pattern of heavy particles even under nonfavorable reactive situations.In contrast, for light particles, inertia can wash out the system even under favorable situations.

FIG. 1 .
FIG. 1. Turbulent velocity field U i, j with u 0 2 ϭ1 and l 0 ϭ6 used in most of the simulations of this paper ͑middle͒.Particle velocity field v(r) in Eq. ͑15͒, in the case of generic heavy particles ␤Ͻ1 ͑first panel͒ and generic light particles ␤Ͼ1 ͑third panel͒.

FIG. 10
FIG. 10.Temporal evolution of the mean concentration with the bistable reaction term.We fix Aϭ6 and a random initial distribution of concentrations below 0.5.Different values of D and K are used.These curves have been obtained by averaging over 10 realizations.Left panel: heavy particles, ␤ϭ0.Right panel: light particles, ␤ϭ2.

FIG. 11 .
FIG. 11.Temporal evolution of the mean concentration with the bistable reaction term.We use a random initial distribution of concentrations with mean values equal to 0.8, 0.6, 0.55, 0.45, and 0.4.The same values for the other parameters are used: Dϭ0.01 and Kϭ3.Thicker lines correspond to the cases with inertia (Aϭ6) and single lines to the cases without inertia.These curves have been obtained by averaging over 10 realizations.All the cases correspond to heavy particles with ␤ϭ0.