History effects in the sedimentation of light aerosols in turbulence: the case of marine snow

3.4 MB·

by  Ksenia Guseva

History effects in the sedimentation of light aerosols in turbulence: the case of marine
snow
Ksenia Guseva,1, Anton Daitche,1 Ulrike Feudel,1,2, and Tamás Tél2,3,
∗ † ‡
1Theoretical Physics/Complex Systems, ICBM, University of Oldenburg, 26129 Oldenburg, Germany
2MTA-ELTE Theoretical Physics Research Group, Eötvös University,
Pázmány P. s. 1/A, H-1117, Budapest, Hungary
3Institute for Theoretical Physics Eötvös University,
Pázmány P. s. 1/A, H-1117, Budapest, Hungary
WeanalyzetheeffectoftheBassethistoryforceonthesedimentationofnearlyneutrallybuoyant
particles, exemplified by marine snow, in a three-dimensional turbulent flow. Particles are charac-
terized by Stokes numbers much smaller than unity, and still water settling velocities, measured in
unitsoftheKolmogorovvelocity,oforderone. ThepresenceofthehistoryforceintheMaxey-Riley
equation leads to individual trajectories which differ strongly from the dynamics of both inertial
particleswithoutthisforce,andidealsettlingtracers. Whenconsidering,however,alargeensemble
ofparticles,thestatisticalpropertiesofallthreedynamicsbecomemoresimilar. Themaineffectof
the history force is a rather slow, power-law type convergence to an asymptotic settling velocity of
6
the center of mass, which is found numerically to be the settling velocity in still fluid. The spatial
1
extension of the ensemble grows diffusively after an initial ballistic growth lasting up to ca. one
0
large eddy turnover time. We demonstrate that the settling of the center of mass for such light
2
aggregatesisbestapproximatedbythesettlingdynamicsinstillfluidfoundwiththehistoryforce,

n on top of which fluctuations appear which follow very closely those of the turbulent velocity field.
a
J

6 I. INTRODUCTION up to now – to our knowledge – only in stratified turbu-
lence [21] and for the plankton problem [22].

]
n There is an increasing evidence, both theoretical and The motivation for our particular range of parameters
y experimental, pointing out the relevance of memory ef- comes from recent studies of marine ecosystems which
d fects in the advection of inertial particles (see e.g. [1– emphasize the importance of marine snow. Marine snow
- 10]). Several further studies concerning these effects in plays a central role in the carbon cycle [23–25], and its
u
l turbulencearereviewedin[11]. Theequationsofmotion formation is mainly due to physical aggregation, a con-
f
. for small spherical inertial particles were formulated by sequence of particle-flow interactions. Sedimentation of
s Maxey and Riley [12] and Gatignol [13] with corrections marine snow is considered to account for a large frac-
c
i by Auton et al. [14] and are of integro-differential type tion of carbon sequestration into the deep ocean [26, 27]
s
in their full form. They contain an integral term which this net ocean sequestration flux is estimated to reach
y
h accountsforthediffusionofvorticityaroundtheparticle 1015 g Carbon/yr [27]. The physical and biological
∼
p throughoutitsentirehistory. Thisintegraltermiscalled properties of marine snow aggregates make it rather dif-
[ thehistory(orBasset)force[15],andithasbecomeclear ficult to estimate their sinking velocity. The techniques

by now that the often used approximation in which this employed to evaluate settling velocities vary across dif-
1
v term is neglected is improper, and the full Maxey-Riley ferent measurements [28, 29], and the results are diffi-
9 equation should be considered [1–11, 16, 17]. cult to compare due to the variation in density and sizes
2 In this work, we analyze the effect of the history force of the used aggregates. An additional difficulty for in
0 on sedimenting particles in turbulence in the presence of situexperimentsisthefactthatturbulentkineticenergy
1 gravity. Previous efforts to understand the importance varies with depth [30]. While some laboratory experi-
0
of the history force in the presence of gravity in smooth mentswithgridgeneratedturbulence[31,32]andinsitu
.
1 flows are due to Mordant and Pinton [7] and to Lohse measurements[33]findindicationsforaretardedsettling
0 and coworkers [18, 19] who also carried out experiments. insitucomparedtolaboratorymeasurementsinstillwa-
6
Their studies, however, concentrated on free sedimen- ter, other observations in coastal areas [34] and in the
1
tation, that is on the particle motion in a fluid at rest, laboratoryusingCouettedevices[35]reportanenhance-
:
v andonbubbledynamicsinastandingwave,respectively. ment of the sinking speed in turbulence.
i
X More recent papers investigate the problem in a station- Marine snow particles contain organic and inorganic
ary [20] and periodically changing cellular flow [4]. The components as primary particles which stick together in
r
a sedimentation problem in turbulent flows is considered afractal-likestructurepossessingarelativelyhighporos-
ity. This fact has been taken into account in concepts
working with an effective density [36] [37, 38], an effec-
tivediameter[39]oramodifiedStokeslaw[40]oftheag-
∗ [email protected] gregates. Moreover, sinking marine aggregates undergo
† [email protected] changes in size and density due to aggregation and frag-
‡ [email protected] mentationprocessesinfluencingthesettlingofthem [41].
2
Those biological properties are difficult to take into and coastal areas [43]. The sizes and the flow set the
account when modeling the sinking of marine aggregates Stokesnumbers(seeEq.(6)below)whichtakethevalues
asinertialparticlesusingtheMaxey-Rileyequation. The St = 0.083 and St = 0.03, respectively. The parameters
fractal shape can be taken into account by means of an characterizing the six cases are summarized in Table I.
effective density [42], but this has been studied so far
only neglecting the history force. Since the latter has 0.08
not yet been formulated for more complicated objects 0.07 Soulsby
than spheres, the study presented here will work exclu- 0.06 McCave
sively with spherical particles the properties of which ) 0.05
3 (IV) Alldredge
m
are based on the effective aggregate diameters and ef- 0.04
c
fective densities given in the literature. The effective g/ 0.03
( (V)
densities of marine aggregates are usually very close to ρ 0.02 (I)
∆ (VI)
the water density and a general property is that larger 0.01 (II)
aggregates have smaller density than small ones. The 0.00 (III)
relationship between the size (average effective radius 0.01
− 100 150 200 250 300 350 400 450 500 550
a) and effective excess density ∆ρ between particle and
fluid (∆ρ = ρ ρ ) for aggregates rich in biological a(µm)
p f
−
components was proposed by McCave et al [43] to be
∆ρ a 1.3 , and this relation was used to fit the ex- FIG. 1. Representation of the chosen parameters (see Table
−
perim∝ental results [44], see middle curve in Fig. 1. This I)ontheparticleradius–excessdensityplane. Thetwolower
curves representtherelationship betweenthe effectiveexcess
relation is close to the one obtained from in situ mea-
density∆ρandtheeffectiveradiusaforaggregateswithpre-
surements from Santa Barbara Channel by Alldredge
dominantlyorganiccomposition[33,43],asexpecetedforthe
et al [33], ∆ρ a 1.6 for marine snow characterized
∝ − openocean,whiletheuppermostcurveshowsthisrelationship
by a > 250 µm, also with predominately organic com-
for aggregates from coastal areas and estuaries [46], contain-
position (lowest curve Fig. 1). On the other hand, in ing a large fraction of inorganic components.
estuaries and coastal regions aggregate composition in-
cludes more inorganic components [45], therefore, they
As a preliminary qualitative analysis, let us concen-
are smaller and slightly denser than the ones formed in
trate here on the particle Reynolds number Re that
theocean. Thecorrespondingeffectivesizeeffectiveden- p
should not exceed a limit. It should be below or of
sity relationship was studied by Soulsby et al [46], and
the order of unity for the Stokesian drag to be valid,
assumes ∆ρ a 0.66 (see uppermost curve in Fig. 1).
− at least in a good approximation. Since gravity breaks
∝
For a review see [47]. Typical velocities in the ocean’s
the isotropy of the advection problem by preferring the
upper layer are strongly dependent on the wind and can
vertical (z) direction, and particles are nearly neutrally
reach up to 0.5 m/s [48]. The turbulent kinetic energy
bouyant, it is worth defining a vertical and a horizontal
(cid:15) typical for the open ocean is (cid:15) = 10 6 m2/s3 [23, 49],
− particle Reynolds number for spheres of radius a and of
which sets the size of the smallest possible eddies, the
typical slip velocity(cid:126)v (cid:126)u relative to the fluid
Kolmogorov length η to be 10 3 m. The size of ag- −
−
∼
gregates(macroaggregates)variesfrom0.1tolessthan1 av u av(cid:126) u(cid:126)
mm[24,50],howevertheaverageaggregatesizeisalways Re∗z = | zν− z|, Re∗h = | hν− h|, (1)
atmostη/2accordingto [24,50],thoughtherelationship
between the average aggregate size and the turbulent ki- where index h refers to the horizontal component, and
netic energy in the ocean is not well established due to ν is the fluid’s kinematic viscosity. The corresponding
the difficulty of in situ measurements. usual particle Reynolds number Re follows from the
p
identity Re2 = Re2 +Re2. The order of magnitude of
Sinceweareinterestedintheeffectofthehistoryforce p z h
the vertical slip velocity is the settling velocity in still
we are confined to a certain, yet realistic, set of param-
water which we write as W = Wu , where W is
settling η
eters for size and density of our marine snow particles
the dimensionless settling velocity taken in units of the
which maximize the impact of the history force. On the
Kolmogorov velocity u (all results for turbulent advec-
η
onehandweneedStokesnumbersthatarenottoosmall
tion will be given in Kolmogorov units). Measuring the
for the history term to play an important role. On the
particle radius in Kolmogorov length, η, we find
other hand, having a small Stokes number also decreases
the impact of preferential concentrations. To study the aWu η a
Re = η =W , (2)
problem,weselectsixdensity-sizepairstypicalofmarine ∗z η ν η
snow. The radii are 0.5 and 0.3 mm, since the strongest
impactofthehistoryforceisexpectedatthelargestsizes, sincethefluidReynoldsnumberontheKolmogorovscale
largest possible Stokes numbers [11]. To both of these u η/ν is by definition unity. The horizontal slip velocity
η
sizes we assign three different densities, see Fig 1. We is expected to vanish with the Stokes number. There-
also display the results of the relationship of the excess fore, the horizontal slip velocity should be proportional
density ∆ρ and particle diameter a for open ocean [46] toStu . Takingtheproportionalityfactortobeunity,we
η
3
case ∆ρ (g/cm3) β a (m) St W Re∗
z
(I) 0.015 0.9900 8.21 4.1
(II) 0.0075 0.9950 5·10−4 0.083 4.12 2
(III) 0.003 0.9980 1.65 0.8
(IV) 0.05 0.9677 9.67 2.9
(V) 0.025 0.9836 3·10−4 0.03 4.91 1.5
(VI) 0.01 0.9934 1.98 0.6
TABLE I. Parameters for six representative cases of marine aggregates in the ocean ((I), (II), (III)) and coastal areas ((IV),
(V), (VI)). Parameters β, St, W and Re∗ are defined in equations (5), (6), (7), and (1), respectively, and turbulence data are
z
taken from Table II.
find for the horizontal Reynolds number in an analogous
manner the estimate
d(cid:126)v 1 W D(cid:126)u 3β t d((cid:126)v−(cid:126)u)
a = ((cid:126)u (cid:126)v)+ (cid:126)n+β dτ dτ, (4)
Re∗h =Stη. (3) dt St − St Dt−(cid:114)πSt(cid:90)0 √t−τ
where (cid:126)n is the vertical unit vector pointing downwards.
This form of the equation holds when the particle is ini-
SinceSt W (seeTableI.),wefindthatRe Re . We
p z
(cid:28) ≈ tialized at time zero with a velocity coinciding with that
think that this is a central property of marine snow sed-
of the fluid, zero initial slip velocity. We have to dis-
imentation, which expresses that these particles behave
tinguish the full derivative along a fluid element and a
horizontally as nearly neutrally bouyant, but they sedi-
particle trajectory, given by
ment with a speed comparable to that of the small scale
fluctuations of the fluid (u ), i.e. they are not neutrally
η D ∂ d ∂
bouyant from the point of view of the vertical dynam- = +(cid:126)u and = +(cid:126)v ,
Dt ∂t ·∇ dt ∂t ·∇
ics. We shall in fact see that the instantaneous particle
Reynolds numbers converge in time towards Re∗z. The respectively. The velocity of the particle changes due to
characteristic numbers Re∗z are also indicated in Table I. the action of different forces. The forces in (4) represent
Thepaperisorganizedasfollows: InSec.IIwepresent from left to right: the Stokes drag, the gravity, the pres-
an overview of the equation of motion with the history sure force (which accounts for the force felt by a fluid
force. Next, we recall an infinite series solution of it in element together with the added mass force), and lastly
still fluid and find a simple analytic approximation to be the Basset history force. The equation is written in di-
valid after relatively short times, presented in Sec. III. mensionlessform,rescaledbytheKolmogorovtimeτ and
Then we summarize the approach used to compute the the Kolmogorov length scale η of the flow (u = η/τ ).
η η
historyforce, andtogeneratetheturbulentvelocityfield The ratio
inSectionsIV.InSec.Vournumericalresultsconcerning
3ρ 3ρ
thesedimentationdynamicsinspacearesummarized. In β = f = f (5)
ρ +2ρ 3ρ +2∆ρ
section VI we turn to results on velocities and accelera- f p f
tions. Sec. VII is devoted to estimating the relevance of
characterizes the excess density of the particle ∆ρ and
the Faxén corrections. Our final conclusions are given in
the density of the fluid ρ . For aerosols β <1 [51].
f
Sec. VIII.
Another dimensionless parameter in Eq. (4) is the
Stokes number
a2 τ
St= = p, (6)
II. EQUATION OF MOTION AND NOTATIONS 3νβτ τ
η η
which is the ratio of the particles’ relaxation time τ due
We analyze the advection of spherical, rigid particles p
to kinematic viscosity ν of the fluid to the Kolmogorov
withasmallparticleReynoldsnumberinanincompress-
time.
ible and viscous fluid. The Lagrangian trajectories of
Additionally, parameter W governs the dimensionless
suchparticlesareevaluatedaccordingtotheMaxey-Riley
settling velocity in still fluid. It can be written as
equation[12,13],includingthecorrectionsbyAutonand
coworkers [14]. In the full Maxey-Riley picture one de- gη
W =St(β 1) , (7)
scribes the dimensionless evolution of the particle posi- − u2
η
tion(cid:126)x(t) and velocity(cid:126)v(t)=d(cid:126)x/dt in a flow field(cid:126)u((cid:126)x,t).
Without Faxén corrections the equation of motion reads where the last factor corresponds to the reciprocal of a
as turbulentFroudenumber. ItistobeemphasizedthatW
4
cannotbevariedfreely: anad-hocchoiceofW toagiven term in the parenthesis it is very close to the exact so-
St could imply that, for a fixed density, the flow and/or lution for t > 22. Note that these forms do not depend
(cid:48)
the gravity g are changed. The W values given in Table ontheparticlesizesinceStokesnumberscanonlybede-
I.,usedbyus,aretheoneswhichfollowfromtheparticle fined in a moving fluid. Whether the particle Reynolds
properties and the characteristics of our turbulent flow. number Re remains small, i.e. whether equation (4) re-
p
Anyhow, in sedimentation the role of the dimensionless mains valid during the entire free fall, should be checked
settling velocity might be more relevant than that of the a posteriori in the knowledge of the dimensional settling
Stokes number. velocity, the particle size and the fluid’s kinematic vis-
We shall compare the Maxey-Riley equation (Eq.(4)) cosity.
to the approximation which does not take into account For comparison, we mention that the solution of the
the history force, widelyusedinertialdynamicsequation(Eq.(8))provides
forthesameproblemalineardifferentialequationwhose
d(cid:126)v 1 D(cid:126)u solution is with the same zero initial condition, and in
= ((cid:126)u (cid:126)v+W(cid:126)n)+β , (8)
dt St − Dt the same units:
often called the advective equation of inertial particles. v (t)=1 e t(cid:48). (12)
We emphasize that Eq.(8) does not follow from any ap- z(cid:48) (cid:48) − −
proximation of the Maxey-Riley equation for our sets of
This solution is of completely different character.
particle parameters, its use is motivated by mere numer-
The solution of the ideal tracer problem (Eq.(9)) is
ical convenience.
that the particle velocity jumps immediately from 0 to
We also carry out simulations with the equation
unity and remains there forever. Note that this behavior
follows from both formulas (11) and (12) in the limit of
(cid:126)v =(cid:126)u+W(cid:126)n, (9)
τ 0,whichisequivalenttotakingt intheseex-
p (cid:48)
→ →∞
pressions. Eq. (12),however,doesnotfollowasanylimit
valid for ideal non-inertial particles. Note that this case
of (10). Fig. 2 provides a comparison of these different
ariseswhenSt 0,andisthelimitofbothequations(4)
→ dynamics.
and (8), with different convergence properties, of course.
1.2
III. SETTLING IN STILL FLUID 1.0
0.8
Theexactsolutionforthesettlinginastillfluid((cid:126)u=0) v′z 0.6
was worked out by Belmonte and coworkers [2]. In
0.4
this case a natural velocity unit is the settling velocity
W ,andtimecanbemeasuredinunitsoftheparti- 0.2
settling
cle relaxation time τp. In these units, the dimensionless 0.0
verticalvelocityv (t)indimensionlesstimet canbeex- 0.0 0.5 1.0 1.5 2.0 2.5 3.0
z(cid:48) (cid:48) (cid:48)
pressed in terms of complementary error functions erfc. t
′
With zero initial velocity it reads in our notation as
FIG. 2. Short-term behaviorof the settling instill fluid ((cid:126)u=
0)inthedifferentdynamicsinvestigated. Withmemory(10):
√3β eα1t(cid:48)Erfc √α1t(cid:48) continuous line; without memory (12): dashed line, and non-
vz(cid:48)(t(cid:48))=1+ α1−α2(cid:20) √α(cid:0)1 (cid:1) idnoetrtteidallipnaer.tiNcloetse(tvhz(cid:48)e=rat1hefrordifft(cid:48)er>ent0,vealsocfiotlileoswpsrefdroicmted(9f)o)r:
− eα2t(cid:48)Er√fcα(cid:0)2√α2t(cid:48)(cid:1)(cid:21) (10) aτpnyantidmWeisnetsttlianng,t.reTsipmecetaivnedlyv.elocityaremeasuredinunitsof
whereα ,α aretherootsofthequadraticequationα2+
1 2
(2 3β)α + 1 = 0 depending only on the density via
−
parameter β. IV. TURBULENT FLOW AND NUMERICAL
By keeping only the leading terms of the power law SIMULATION
expansion of the function euErfc(√u) for large u (long
times t), we find Weconsiderherethecaseofparticlesmovinginstatis-
(cid:48)
ticallyhomogeneous,isotropicandstationaryturbulence
3β (3β 2) [52]. To this end we solve the vorticity equation, which
vz(cid:48)(t(cid:48))=1− πt 1− 2−t . (11) is equivalent to the incompressible Naiver-Stokes equa-
(cid:114) (cid:48) (cid:18) (cid:48) (cid:19) tion, on a grid in a triply-periodic box of size L . The
box
Thisformturnsouttoprovidearatheraccurateapprox- energy is injected by a large scale forcing, see [11]. For
imation for t > 2, and even by neglecting the second the integration of the flow we use a standard dealiased
(cid:48)
5
Reλ Lbox/η L/η λ/η ∆x/η Tsim/τη T/τη ∆t/τη urms/uη N3
112 633 156 20.9 1.24 1020 29.0 0.015 5.39 5123
TABLE II. Parameters of the simulated turbulent flow: Taylor Reynolds number Reλ = λurms/ν, size of the periodic box
L ,integralscaleL=u3 /(cid:15),Taylormicroscaleλ=u (cid:112)15ν/(cid:15),sizeofagridcell∆x,lengthofthewholesimulationT ,
box rms rms sim
(cid:112)
large-eddy turnover time T =L/u , time step ∆t, root-mean-square of the velocity u = (cid:104)(cid:126)u2(cid:105)/3, number of grid points
rms rms
N3. All dimensional quantities are given in multiples of the corresponding Kolmogorov units.
Fourier-pseudo-spectral method [53, 54] with a third- respectively.
order Runge-Kutta time-stepping scheme [55]. The val- Trajectories with the same initial condition, but fol-
ues of Eulerian quantities, which are available on a grid, lowing these three distinct dynamics, deviate from each
are obtained at the particle positions through tricubic other already after a short period of time. The distance
interpolation. The characteristics of the turbulent flow amongthesetrajectoriesincreasessignificantlywithtime
and the simulation parameters are depicted in Table II. in both the horizontal and the vertical directions. This
Since the Kolmogorov scale is η =(ν3/(cid:15))1/4 [52], a fixed resultsinstrongdifferencesinthepredictionsforthepo-
value of it can belong to any kinematic viscosity ν and sition of a particle, since the trajectories of different dy-
mean energy dissipation (cid:15), as long as the ratio ν3/(cid:15) is namics can be thousand Kolmogorov lengths away from
fixed. For the particular choice of η = 1 mm, which we each other after 500τ as Fig. 3a,d illustrates [56].
η
shall take as a typical value in our estimations, one finds
Although there are strong differences for the predic-
with the viscosity of water (cid:15) 10 6 m2/s3.
∼ − tions of the position of an individual particle for these
The presence of the history integral in (4) leads to
three dynamics Fig. 3d, these differences are smaller
two problems from the numerical point of view. First,
when an ensemble is considered Fig. 3e,f. For this anal-
the singularity of the history kernel impedes an accurate
ysis we initialize clouds of particles, one with smaller
numerical solution. This problem can be solved by the
and other with larger number of particles, with the ini-
use of a specialized integration scheme [17] which treats
tial condition mentioned above, and evolve them accord-
the history force appropriately. This third order scheme
ing to our three possible dynamics. The center of mass
has been adjusted for our purposes, see [11] for details.
of each cloud also follows a distinct trajectory, however
Second, it is necessary to recompute the history integral
the distances among the centers of mass do not grow as
foreverynewtimestep. Thisleadstohighcomputational
fast as that of the trajectories of individual particles, see
costsandahighdemandformemory(tostorethehistory
Fig. 3b,c where the final horizontal difference is of a few
of each particle). This second problem is inherent to
Kolmogorov lengths only. Moreover, in the x,y planes
thedynamicswithmemoryand,asaconsequence,limits
(upper panels a–c), it becomes clear that the total dis-
us to a moderate number of particles. For each case of
placementsareinratherdifferentdirectionswiththedif-
particleparameterswesimulatedN =1.5 105 particles.
p · ferent dynamics. The horizontal distances between the
The initial particle positions have been chosen randomly
centers of mass decrease with the number of particles,
andhomogeneouslydistributedinthetriple-periodicbox
which can be considered as a consequence of the law of
ofsizeL ofthesimulation; theinitialparticlevelocity
box large numbers. The total horizontal displacement in all
is that of the fluid at the particle’s position.
three dynamics is therefore expected to be zero in the
large particle number limit [57].
V. TURBULENCE: RESULTS ON THE After having seen the results for the center of mass
POSITION OF PARTICLES of the particle ensembles, we show in Fig. 4 their dis-
tribution in space at three different time instants. It
We start by comparing individual trajectories in the is clear that with large settling velocities the ensemble
Maxey-Riley equation (4), in the inertial equation (8) in blobs are well separated after 500 time units, this sepa-
which memory is neglected and in the non-inertial dy- ration decreases, however, with W, and with the small-
namics (9). Throughout the paper we will use the fol- est settling velocity there is hardly any separation, the
lowing notation for particles following the different dy- blobs strongly overlap, and some points are even above
namical equations (4), (8), and (9), and show their cor- thecubeofinitialconditionsafter1020timeunits. These
responding curves with particular line types: resultsareobtainedinthepresenceofmemoryeffectsbut
we generated the corresponding figures without memory
particles with memory, continuous line, computed and with non-inertial particles (governed by Eqs.(8) and
•
by (4), (9), respectively), too. No difference can be recognized
by naked eyes. This is the first hint to the fact that in-
particleswithoutmemory,dashedlinecomputedby
• spite of the difference in the individual and in the center
(8), and
of mass trajectories following from the different dynam-
non-inertial particle, dotted line, computed by (9), ics,thestatisticalpropertiesarequitesimilar. Toseethe
•
6
(a) (b) (c)
800 10 10
8
600
6 5
400 4
η[] 200 η[] 02 η[] 0
y 0 y 2 y
−4 5
−200 −6 −
400 −8 10
− 200 0 200 400 600 800 − 10 5 0 5 10 − 10 5 0 5 10
− − − − −
x[η] x[η] x[η]
(d) (e) (f)
withmemory (I)
2000 withoutmemory 2000 2000
(II)
non-inertial
(V)
0 0 0
2000 2000 2000
− − −
η]
[
z 4000 4000 4000
− − −
6000 6000 6000
− − −
8000 8000 8000
− − −
10000 10000 10000
− 200 0 200 400 600 − 10 5 0 5 10− 2 1 0 1 2 3
− − − − −
x[η] x[η] x[η]
FIG. 3. (left) Individual trajectories of the same particle started with the same initial conditions (with zero slip velocity)
with the parameter set of case (I), however, following three distinct equations of motion (x = 293.16 η, y = 214.18 η,
0 0
z = 105.71 η). (middle and right) Position of the center of mass of the ensemble containing 10% of our standard particle
0
number(Np =1.5·104),andthestandardnumberNp =1.5·105 ofparticles,respectively,forthecase(I),(II),and(V)evolved
with the different equations of motion up to 1020 τη. Here and in the following figures we use the convention that x and t
denote the dimensional space and time, respectively, and the dimensions are given in parentheses.
differences, more quantitative methods should be taken. blobsislargerthanthehorizontaloneatanyinstant. The
data are plotted on a log-log scale to enlighten the ap-
In Figure 5 we present a time-continuous plot of the
pearance of power law behavior. The two black straight
z-coordinate of the center of mass for the different cases
lines represent ballistic (σ2 t2) and diffusive (σ2 t)
with the three different dynamics. Here differences be- ∼ ∼
spreading. A crossover to the diffusive behavior can be
come visible, and are on the order of a few hundred Kol-
observed at t 20 30 time units. It is natural to un-
mogorovlengthsattheendofthesimulation. Thediffer- ∼ −
derstand that when the blobs are large, the ensembles
ence is, however, never larger than a few percent of the
become subjected to a diffusive spreading by wandering
instantaneous value of z . The long-term behavior is a
(cid:104) (cid:105) in-between the largest scale vortices.
roughly linear increase in all cases, indicating a rather
uniform settling. Note that the graphs for the largest
Inordertoexplaintheballisticbehavior, werecallthe
excess density (case (I) and (IV)) are close to each other
theory of Batchelor [58] for the separation of pairs of
in spite of the different Stokes numbers. The cases with
ideal tracers in three-dimensional homogeneous isotropic
intermediate and small excess densities behave also sim-
turbulence. Thistheoryclaimsthatthemeansquaresep-
ilarly.
aration should grow as t2 for times shorter than a char-
As seen from Fig. 4, the blob sizes are also impor- acteristic time t . For times larger than t the famous
0 0
tant. To monitor their time evolution, we determined Richardson scaling [52] should hold characterized by a
the standard deviation σ about the center of mass in the scaling proportional to t3. This regimes extends, how-
horizontal x direction, and in the vertical direction, as ever, only up to the time when the effect of the largest
Fig. 6 shows. The horizontal and the vertical behavior coherent structures becomes dominant, i.e. up to the
are different, reflecting again that gravity prefers a cer- eddyturnovertimeT. Thecharacteristictimet depends
0
tain direction. In fact, the vertical extension of all the ontheinitialspatialseparation(cid:126)r betweenthetwoparti-
0
7
(IV) (V) (VI)
1000.0 1000.0 1000.0
-1500.0 -1500.0 -1500.0
-4000.0 -4000.0 -4000.0
z[η] z[η] z[η]
-6500.0 -6500.0 -6500.0
-9000.0 -9000.0 -9000.0
-11500.0 -11500.0 -11500.0
-14000.0 -14000.0 -14000.0
-40x00-20[0η0] 020004000 -4000-2000020004000y[η] -40x00-20[0η0] 020004000 -4000-2000020004000y[η] -40x00-20[0η0] 020004000 -4000-2000020004000y[η]
FIG.4. Spatialdistributionofthesedimentingparticleensemblesforcases(IV),(V)and(VI)simulatedwiththeMaxey-Riley
equation (4). Black, blue, and red dots represent the location of the particles at times t=0, 510τη, and 1020τη, respectively.
A 2D histogram of the particle density is projected onto the (y,z) plane, curves represent isolines of densities.
10000 (a) 106 (b) 106
(I) (IV) 2] (I) (IV) 2]
8000 (II) (V) η[ 105 ((IIII)I) ((VVI)) η[ 105
6000 (III) (VI) 0) 0)
tη)[]i 4000 2σ(x−110034 t2 t1 2σ(z−110034 t2 t1
( ) )
zh 2000 t( t(
2x 2z
σ 102 σ 102
0 100 101 102 103 100 101 102 103
t[τ ] t[τ ]
2000 η η
− 0 200 400 600 800 1000
t[τη] FIG. 6. Time-dependence of the variances in the horizontal
(x,leftpanel)andvertical(zrightpanel)directionsinthedif-
FIG.5. Time-dependenceofthez-coordinateofthecenterof ferent cases (coloring and line types as in as in the previous
massinthedifferentdynamics,distinguishedbydifferentline figure)o√nlog-logscales. Forclarity,theinitialvarianceσx(0)
types, and for all six cases distinguished by different colors, (=633/ 12=183indimensionlessunits)issubtracted. The
dashedanddottedlinesarehardtodistinguishedinthisrep- graphsofdifferentdynamicsareoverlaidanddashedanddot-
resentation. Thetwo-sidedarrowsindicatethetypicalspatial tedlinesarehardtodistinguishedinthisrepresentation. The
extension in z dimensions of the ensemble for case (I) at the thin continuous black lines have slopes 2 and 1, respectively,
given instances. to guide the eye.
cles. In dimensional units t0 =( (cid:126)r0 2 /(cid:15))1/3. Hence, for √3L /√12 = √3σ (0)η, where σ (0) denotes the di-
| | box x x
an ensemble of particles with different initial distances
mensionless variance in x direction, used in the plots of
no unique t0 can be found, so that only a typical t0 can Fig. 6. To estimate the dimensionless t /τ we replace
0 η
be estimated. (cid:126)r 2 by 3σ2(0)η2 to find
Although the original theory applies to ideal, i.e. | 0 | x
non-settling tracers, it is worth estimating t . For
our initial ensemble a natural choice is the var0iance of t0 = 3σx2(0)η2 1/3 =(3σ2(0))1/3 =(3 1832)1/3 =46,
their positions in the initial cube of size L , what is τ (cid:15)τ3 x ×
box η (cid:18) η (cid:19)
8
whereweusedthat τ =(ν/(cid:15))1/2 and η =(ν3/(cid:15))1/4 [52].
η 0.0
This value turns out to be larger than the dimensionless
turnover time, T/τ which is about 30 (see Table II). 0.2
η −
Thus, there is no possibility for seeing the Richardson W
− 0.4
scaling due to the broad initial distribution of the par- piW −
ticles. The anisotropy of the problem is reflected in the vzh −0.6 (I) (V) ESqt.1=2,0.083
fact that the crossover to the diffusive behavior occurs (II) (VI) Eq.10,
somewhatlaterintheverticalthaninthehorizontal. We −0.8 ((IIIVI)) ESqt.1=0,0.083 SEqt.1=2,0.03
St=0.03
1.0
haveverifiedthattheevolutionofaninitiallystronglylo- − 0.0 0.2 0.4 0.6 0.8 1.0
calized ensemble leads to Richardson behaviour (see Ap- t[τη]
pendix A).
It is worth mentioning a related problem. Single par-
FIG. 7. Short-term behavior of the center of mass velocities
ticle dispersion was investigated in stratified turbulence expressed as ((cid:104)vz(cid:105)−W)/W: empty symbols: with memory
byvanAartrijkandClercxinthepresenceofthehistory (Eq. 4), and full symbols: without memory (Eq. 8). Black
force [21], but with larger typical excess densities and stars indicate the results for non-inertial particles but only
Stokes numbers compared to ours. They also found a up to t=0.3τη in order to avoid heavy overlap. Continuous
crossover between ballistic and diffusive spread, i.e. the (dashed) curve represents the still fluid result with memory
history force did not change the exponents. (without memory) as expressed by (10) ((12)).
exactly on the still-fluid curves.
VI. TURBULENCE: RESULTS ON VELOCITIES To see the long-term behavior, and check if the rela-
AND ACCELERATIONS
tionwiththestillfluidresultsholdalsoonthistimescale,
we show in Fig. 8 the difference between W and the en-
Inertial effects are more pronounced in the velocity sembleaveragedverticalvelocityupto1020Kolmogorov
data characterizing the ensemble. An investigation of times. Since the large eddy turnover time T in Table II
the short-term behavior, up to a single time unit (one isabout30τ ,t=1020correspondstoabout34turnover
η
Kolmogorov time), indicates clearly that particles for times, quite a considerable time span in turbulence. The
which the history force is neglected approach typically results of two cases (I and III) are shown in Fig. 8a,d
much faster the asymptotic settling velocity than those fortheMaxey-Rileydynamicswithmemory,inFig. 8b,e
for which the history force is taken into account. This for the dynamics without memory and in Fig. 8c,f for
can very well be seen in Fig. 7 which exhibits the z- the non-inertial dynamics. With memory, the ensemble
component of the particle velocities averaged over the averaged settling velocity is always below W, and has
ensemble, for all six cases, and for all three types of dy- not yet reached a steady value by the end of the inves-
namics. At t = 0.4 the particle dynamics without mem- tigated time interval. This is so even for averages taken
ory indicates a settling with W for all the cases, without over finite time windows of, say, one large eddy turnover
any further change, while the results following from the time. Such smoothed time series (not shown) are, how-
Maxey-RileyequationpredictasettlingwithaboutW/2, ever, remarkably close to the result valid in still fluid.
with a difference between the cases of different Stokes In order to check if this property is not a consequence
numbers,andamonotonousincreasefort>0.4. Onthis of the relatively large settling velocities, we carried out
scale no difference can be seen between the cases with additionalsimulationswith10timessmallerW-sbutthe
differentexcessdensities∆ρ. Itisinterestingtocompare sameStokesnumbersasinTableI.Theresultsinturbu-
thenumericaldatawiththeanalyticexpressionpresented lentflowarefoundtocorrelatewiththestillfluidsettling
for the free fall in still fluids in Section III. To this end, just as in Figs. 7, 8. In fact, the red lines represent the
we have to rescale equations (10) and (12) according to function:
theunitsusedfortheturbulentflow. Sincethetimeunit
instillfluidcanonlybeτ ,butinturbulenceitischosen 3St
p v (t)=W 1 , (13)
to be the Kolmogorov time τη, and the Stokes number is (cid:104) z(cid:105) (cid:34) −(cid:114) πt (cid:35)
exactly τ /τ (see Eq.(6)), the dimensionless time t of
p η (cid:48)
those equation should be transformed into a t/St, where which follows from (11) to be valid asymptotically, and
t is the dimensionless time used in all our equations. Si- the deviation of β from unity can be neglected since all
multaneously, v ofthestillfluidcaseshouldbereplaced our access densities are rather small. This indicates that
z(cid:48)
by Wv in order to be converted to our units. The dif- the decay towards the asymptotic settling velocity is of
z
ferentcurvesinFig.7representthestillfluidresults(10) power-law type, decaying as one over the square-root of
and (12) in these units. Since our six parameter sets are time. Since such functions are scale-free, no characteris-
grouped around two Stokes numbers, with which time is tic time can be associated with them (in contrast e.g. to
scaled, each type of solution appears with two curves. A exponential decays).
surprising observation is that all points representing the Regarding the right column of panels, note the very
ensemble averages (symbols) of the turbulent results fall small scale on the vertical axes. In all cases the average
9
(a) (b) 0.08
0.0 ]η 0.06 withoutmemory
u 0.04
[
W 0.02
0.2 − 0.00
−
u[]η 0.4 vzhi−00..0024
W − −
− (c) 0.08
i
vzh−0.6 ]η 0.06 non-inertial
u 0.04
[
W 0.02
0.8
− − 0.00
1.0 withmemory stillfluid vzhi−00..0024
− 0 200 400 600 800 1000 −
0 200 400 600 800 1000
t[τ ]
η t[τ ]
η
(d) (e) 0.08
0.0 ]η 0.06 withoutmemory
u 0.04
[
W 0.02
0.2 − 0.00
−
u[]η 0.4 vzhi−00..0024
W − −
− (f) 0.08
i
vzh−0.6 ]η 0.06 non-inertial
u 0.04
[
W 0.02
0.8
− − 0.00
1.0 withmemory stillfluid vzhi−00..0024
− 0 200 400 600 800 1000 −
0 200 400 600 800 1000
t[τ ]
η t[τ ]
η
FIG. 8. Long-term behavior of the velocity difference ((cid:104)vz(cid:105)−W) for cases I (a,b,c) and III (d,e,f). The left column (a,d)
showstheresultsoftheMaxey-Rileyequation,therightcolumnthosewithoutmemoryeffects. Continuousredlinesrepresent
(13) an approximate form of settling in still fluid, and fits nevertheless very well to the turbulent data.
is zero, meaning that the average settling velocity over muchsmallerthanunity. Theydependontheaccessden-
theinvestigatedtimeintervalisW,asinstillwater. The sity, and differ a little bit with and without memory. In
graphs of the non-inertial and memoryless dynamics are any case they happen to be close to the estimated value
somewhat different, but they basically represent a ran- Re = St a/η. The difference is changing with W, and
∗h
dom process around zero. Fluctuations in all panels are we can discover a simple relation
on the order of 0.05. These features also hold for the
Re Re
results of the other four cases not shown here. hw− hm =const. (14)
W
We also evaluate the slip velocities of the particle en- to hold, where index m and w stand for memory and
semble as time series with and without memory and, without memory, respectively. The value of the constant
based on these define an instantaneous vertical and hor- isfoundtobeabout0.005and0.0006forSt=0.083and
izontal Reynolds number Re (t) and Re (t) in an analo- St = 0.03, respectively. It reflects that for W 0 the
z h
→
gouswayasRe andRe aredefinedin(1),butthistime particles have smaller and smaller excess densities, and
∗z ∗h
withtheaverageofthemodulusoftheinstantaneousslip their dynamics approaches that of ideal fluid elements
velocity (cid:126)v(t) (cid:126)u(t) . The results obtained for cases V with the zero initial slip velocity condition used in this
(cid:104)| − |(cid:105)
andVIaresummarizedinFig.9. Withmemory,thever- paper.
tical Reynolds number converges according to a power It is worth also considering the distributions (pdf-s) of
law to a long-term limit, which is close to Re =Wa/η, the different types of accelerations. In Fig. 10 the accel-
∗z
avaluetheotherdynamicsreachpracticallyimmediately. eration due to the drag, pressure and history force are
TheorderofmagnitudeofthelimitingReynoldsnumber plotted, for case IV, at time instants t = 10.5 τ and
η
is unity in all the cases (the values coincide with those t=1020τ for the vertical, and only for the last instant
η
given in Table I). The horizontal Reynolds numbers are for the horizontal components. In the vertical, the drag
10
(a) (b) theaccelerationsthemselveswithouttakingthemodulus
0.018 would all be zero with the exception to the vertical drag
1.4 0.016
and vertical history acceleration.
0.014
1.2 0.012
ez 1.0 eh 0.010
R R 0.008 VII. ESTIMATING THE RELEVANCE OF THE
0.8 0.006 FAXÉN CORRECTIONS
0.004
0.6 0.002
0.000
TheFaxéncorrectionsarecorrectionsto(4)duetothe
0 100200300400500 0 100200300400500
finitesizeoftheparticleandtothecurvatureoftheflow.
t[τ ] t[τ ]
η η They appear as terms proportional to a2∆(cid:126)u which are
(V)withoutmemory (VI)withmemory added to the slip velocity in the Stokes drag and in the
(V)withmemory Re∗h nominator of the memory integral, as well as to the fluid
(VI)withoutmemory velocity in the added mass term [12, 13]. They appear
with a coefficient 1/6 and 1/10, respectively. We con-
centrate here on the correction to the slip velocity and
FIG. 9. Vertical (left panel) and horizontal (right panel) consider the ratio of the average modulus of the correc-
Reynolds numbers as a function of time for case V (upper tion to that of the slip velocity
curves) and VI (lower curves). The grey horizontal line rep-
resents the order of magnitude estimate Re∗h of (1). C = a2 (cid:104)|∆(cid:126)uj|(cid:105) , (15)
j
6 (cid:126)v (cid:126)u
j j
(cid:104)| − |(cid:105)
where index j stands for the Cartesian components x,y
dominates, and has a rather narrow distribution. This is or z in this correction factor. Because of the anisotropy
due to the fact that the slip velocity becomes quickly to due to gravity, it is worth treating the horizontal and
beoforderW. Ontheotherhand,thepdfoftheacceler- verticalcomponentsseparately. Theirdifferencebecomes
ation from the pressure term is rather broad and hardly clear from a simple estimation. Since the characteristic
changes with time after t = 10.5 τ . These two pdf-s length and velocity scale of the turbulent flow are η and
η
are found nearly identical with those in the memoryless u , respectively, the Laplacian in any component can be
η
equations. Onlythepdfofthehistoryforce(red)changes estimated as u /η2. The slip velocity in the vertical is
η
with time rather dramatically: at t = 10.5 τ it is sharp approximatelyWu ,whilethatinthehorizontalisStu ,
η η η
an has a much larger average than the pressure contri- asusedin(3). Wethusfindtheestimatesforthevertical
bution. Bytheendoftheobservationalperiod, however, and horizontal correction factors
the pdf broadens and becomes shifted towards smaller
2 2
1 a 1 1 a 1
values. Itsaverageremainsonlyslightlylargerthanthat C = , C = .
of the pressure. In the horizontal, the distributions are z∗ 6 η W h∗ 6 η St
(cid:18) (cid:19) (cid:18) (cid:19)
similar, do not change too much in time, the averages
Since W is larger than unity in our cases, but St < 1
are ordered as drag, pressure and history with not very
(see Table I.), the relative importance of the Faxén cor-
much differences. These pdf-s are rather different from
rectionsisexpectedtobemuchsmallerinverticalthanin
those obtained without gravity in [11]: all distributions
horizontal direction. For our largest particles a/η =1/2,
are broad there, the pressure contribution is the largest,
andSt=0.083,thustheestimateC amountstoavalue
those of drag and history are comparable, and the aver- h∗
0.5.
ages are not separated by several orders of magnitudes.
We numerically determine the correction factors C as
j
The closeness of the averages resembles Fig. 10c.
functions of time. The results in the presence of mem-
To gain insight into the full time dependence, we plot ory are shown in Fig. 12 for C and C for cases (V)
z x
inFig.11theensembleaverageofthepdf-sshownabove. and(VI)withandwithoutmemory. Theverticalcorrec-
Astrikingfeatureintheverticalcomponents(leftpanel) tions (Fig. 12a) appear to be at most 0.1 % consistently,
is the monotonous decay of the history force. Sooner or withhardlyanydifferencewithandwithoutmemory. For
later, the average of the history force is likely to become smallerexcessdensity(case(VI))thecorrectionislarger.
smaller than that of the pressure force. This internal The measured values are about a factor 5 smaller than
degradationofthehistoryforceseemstobespecifictothe theestimatesC . Inthehorizontal,thecorrectionsfactor
z∗
sedimentation dynamics with all the parameters investi- reachesnearly20%, butisaboutafactor3smallerthan
gated. Horizontally (right panel), however, everything is what the estimate C predicts. There is a measurable
h∗
stationary after t = 10.5. The relatively small values of difference in the correction with memory and without,
the history acceleration explain why no sign of a slow and the former one is consistently larger by about 10
convergence is seen in Fig. 9b. The stationarity of the %. The effect for lighter particles is here stronger again.
average pressure acceleration indicates the stationarity A comparison with Fig. 9 reveals that the tendencies in
of our turbulent flow. Its order of magnitude is indeed the Reynolds numbers and in the correction factors are
u /τ = η/τ2. It is worth noting that the averages of roughly the opposites.
η η η