Effect of turbulent fluctuations on the drag and lift forces on a towed sphere and its boundary layer

0.39 MB·English

by  Holger Homann

Immersed boundaries with a Fourier-spectral method
Holger Homanna, J´er´emie Beca, Rainer Grauerb
aUniversit´e de Nice-Sophia Antipolis, CNRS, Observatoire de la Coˆte d’Azur, BP4228, 06304 Nice, France.
bTheoretische Physik I, Ruhr-Universita¨t, 44780 Bochum, Germany.
0
1
0 Abstract
2
A novelscheme is presentedfor the direct numericalsimulation of two-waycoupling between a particle and
c
anincompressiblefluidflowmaintainedinahigh-Reynolds-numberturbulentregime. Themainideaconsists
e
D incombiningaFourierpseudo-spectralmethodforthefluidwithanimmersed-boundarytechniquetoimpose
theno-slipboundaryconditiononthesurfaceoftheparticle. Thisschemeisshowntoconvergeasthepower

3 3/2ofthespatialresolution. Thisbehaviorisexplainedbythe L convergenceofthe Fourierrepresentation
2
2
of a velocity field displaying discontinuities of its derivative. Benchmarking of the code is performed by

measuring the drag and lift coefficients and the torque-free rotation rate of a spherical particle in various
]
n configurations of the carrier flow. Such studies show a good agreement with experimental and numerical
y
measurements from other groups, and validate the code. A study of the turbulent wake downstream the
d
sphereis alsoreported. The meanvelocitydeficit is shownto behave asthe inverseofthe distance fromthe
-
u particle,as predicted from classicalsimilarity analysis. This law is reinterpretedin terms ofthe principle of
l
f “permanence of large eddies” that relates infrared asymptotic self-similarity to the law of decay of energy
.
s in homogeneous turbulence.
c
i Keywords: Fourier-spectral methods, penalty methods, flow-particle interaction
s
y
h
p
1. Introduction
[

1 The dynamics of suspended particles whose size is larger than the smallest active scales of the carrier
v turbulent flow is important for the understanding of many natural and industrial applications. Let us
5 mention two concrete problems borrowed from the study of planet formation and of rain processes. In the
0
first situation, the role of turbulence is still unsettled in the dynamical mechanisms that take place during
2
the origination of planets in a proto-planetary disk (see, e.g., de Pater and Lissauer, 2001). A key issue is
5
. tounderstandwhetherornotthe diskturbulencecangenerateanincreaseofthe rateatwhichalargebody
2
can accrete dust. In other terms, it requires identifying the turbulent fluctuation mechanisms that lead
1
to a snowball effect where the larger bodies have a growth speed all the more fast as their size increases.
0
1 The same kind of question arises in the precipitation process when a large raindrop falls under the effect
: of gravity and collides with much smaller droplets that are transported by the turbulent airflow in the
v
cloud (see Shaw, 2003). For this problem, another question is to know the efficiency of wet deposition (or
i
X scavenging) of particle pollutants that are present in the atmosphere (see, e.g., Seinfeld and Pandis, 1998).
r These various phenomena are still poorly understood and are subject to empirical ad-hoc modeling. Such
a
fundamental open questions, which are taken from natural sciences but also arise in a number of industrial
settings, reflect an inadequate current level of modeling for large particles in turbulent flows.
Current approaches for the dynamics of finite-size particles are commonly assuming that the particles
are infinitely small (point particles) and that the Reynolds number associated to their interaction with the
flow vanishes, so that their dynamics can be obtained by matching the solution to the Stokes equation
in their vicinity to an unperturbed fluid velocity. Recent experimental work showed that the dynamics
of finite-size particles is not well described by such models (see Qureshi et al., 2007; Xu and Bodenschatz,
2008;Volk et al.,2010). Significantsizeeffectshavebeenforinstancemeasuredintheaccelerationstatistics
of such particles. Despite the fact that the flatness of the acceleration distribution depends only weakly
Preprint submitted toElsevier December 24, 2010
on the diameter of the particle, its variance strongly decreases with diameter and temporal correlations of
the accelerationvary withthe particle size in a manner thatcannotbe predictedby standardpoint-particle
models. Improvedmodelsaccountingfortheso-called“Fax´encorrections”(spatialaveragesoffluidvelocities
andaccelerationsoverthevolumeoftheparticle)areinmuchbetteragreementwithsuchqualitativetrends
but hardly reproduce quantitatively the mentioned experimental data (see Calzavarini et al., 2009).
One of the mostchallenging requirementfor improvingsuch models is to geta preciseunderstanding on
howthecarrierflowismodifiedinthevicinityoftheparticleinsituationswheresucheffectsoccursonscales
which arecomparable to the active scales of the surrounding turbulence. For instance, it is clear that much
workisstillneededinordertogetabetterhandlingonhowthekineticturbulentenergy,thedissipationrate,
aswellasotherstandardstatisticalpropertiesofturbulence,variesintheturbulentboundarylayerenclosing
theparticle. Directnumericalsimulationsappeartobevaluablecandidatesastheygiveauniqueglobaland
instantaneous access to such quantities. However, the numerical problem is challenging for several reasons:
(i) the numerical approach has to resolve all active scales of a fully developed turbulent flow, which span
several decades; (ii) at small scales the boundary layer around the particle is of great importance and has
to be well resolved; (iii) if one considers a moving particle the algorithm has to follow the displacement of
the particle accurately andefficiently; (iv) for a sufficient separationof the differentscales involved(forcing
scales, inertial range of scales, dissipative scales) the algorithm has to parallelize efficiently on massively
parallel computers.
Two different strategies have been recently proposed in order to solve the incompressible Navier–Stokes
equationstogetherwithno-slipboundaryconditionsatthe surfaceofa sphericalparticle. Bothmakeuse of
ahomogeneousgridcoveringofthecompletephysicaldomain. Thefirst,“Physalis”,whichwasdevelopedby
Prosperetti and Oguz (2001) and extended and used by Naso and Prosperetti (2010), makes use of known
analytic solutions to the Stokes equation for the boundary layer in the vicinity of the particle, in order
to describe the velocity at the nearest surface points of the grid. Iteratively, the analytic solution and
the solution on the homogeneous grid are matched. An advantage of this strategy is that the no-slip
conditionis imposedinanexactmanner. However,improvementare requiredinorderto achievevery-large
Reynolds number carrier flows (the value Re = 20 of the Taylor-microscale Reynolds number was chosen
λ
in Naso and Prosperetti, 2010). The other recently developed approach consists in using and matching
two different grids. A spherical one encloses the particle surface and is designed to resolve the boundary
layer and a second homogeneous grid covers the entire domain. This “Overset Grid” technique proposed
by Burton and Eaton (2005) uses a third-order interpolation to transfer values from one grid to the other
in order to couple the two solutions. This method was shown to reproduce very accurately experimental
measurementsof the drag coefficientat a particle Reynolds number ofRe =20. In turbulent settings, this
p
method was used with a Reynolds number Re =32.
λ
Both of these methods reside on finite-difference schemes and are at the moment limited to rathersmall
carrier-flow Reynolds numbers. We propose here a new method that allows one to reach much higher
Reynolds numbers. The idea is to combine a standardpseudo-Fourier-spectralwith a penalty method. The
former is well adapted to incompressible homogeneous and isotropic turbulence and has been extensively
used in such settings during the past. Its advantages comprise its degree of accuracy and its performance
onmassiveparallelsupercomputerssuchasthe BlueGene/PsystemorlargeIntel/AMDclusters. Incombi-
nation with an immersed boundary and penalization strategy, high Reynolds numbers turbulent flow with
anextended inertial rangeof scales can be reached. The proposedmethod is a modificationof that usedby
Homann and Bec (2010) for investigating the dynamics of neutrally buoyant particles suspended in turbu-
lent flows. The former version of the method was well-adapted to particles displaying a small slip velocity
with respect to the fluid but happened to show some limitations when the slip is too large. We present
here a modification of this method that unravels such shortcomings. This paper is organized as follows:
in Section 2 we present in details the proposed method, together with some considerations concerning its
convergence. Section 3 presents measurements of drag and lift coefficients in various settings and compare
them to other numerical and experimental works. We consider the cases of a stationary sphere in a flow
at a constant speed, of a freely rotating sphere and a translating sphere in a flow maintained at rest. Sec-
tion 4 deals with the characterization of the turbulence decay in the particle wake. Section 5 encompasses
concluding remarks and open questions.
2
2. The numerical method
The dynamics of an incompressible flow in which is embedded a particle is given by the Navier-Stokes
equations
∂ u=L(u)=−u·∇u−∇p+ν∇2u+f , ∇·u=0, for x∈Ω\Ω (1)
t e p
together with the no-slip boundary condition on the surface of the particle
u(x,t)=v (x,t) for x∈∂Ω (t). (2)
p p
f denotes here external forces that are exerted on the flow, Ω is the full physical domain containing both
e
the fluid and the particle, ∂Ω (t) is the boundary of the physical location Ω (t)⊂Ω of the particle at time
p p
t, and v denotes the velocity of the particle surface. Note that in the case of a solidsphericalparticle with
p
diameter d, one has Ω = {x : |x−X (t)| ≤ d/2} and v (x,t) = V (t)+[x−X (t)]×Ω (t), where V
p p p p p p p
andΩ arethe translationaland angularvelocities ofthe particle,respectively. The particlemotionis then
p
obtainedby solvingNewtonequations wherethe forceexertedbythe fluid isobtainedby integratingonthe
surface of the particle the contribution from pressure and viscous stress to the fluid stress tensor.
The numerical method presented here consists in integrating the Navier–Stokes equation (1) by using
a strongly-stable, low-storage Runge-Kutta scheme of third order that was introduced by Shu and Osher
(1988). It reads
1 2 2 3 1 1
un+1 = un+ u′′+ ∆tL(u′′), where u′′ = un+ u′+ ∆tL(u′) and u′ =un+∆tL(u0). (3)
3 3 3 4 4 4
∆t denotes here the time step and un and un+1 are the velocity fields at time steps n and n+1 and the
operator L denotes the right-hand side (rhs) of (1). The later is integrated on a homogeneous rectangular
periodic spatial grid with a Fourier pseudo-spectral solver using the standard 2/3 rule to avoid aliasing
due to its quadratic nonlinearity (see, e.g., Patterson and Orszag, 1971). The pressure term is obtained by
solving the Poisson equation
∇2p=−∇·(u·∇u), (4)
directly in Fourier space. The viscous term is treated implicitly via an exponential factor.
Toimposetheno-slipboundaryconditions(2)atthesurfaceoftheparticle,wemakeuseofanimmersed
boundary technique. It consists in introducing in the right-hand side of the Navier-Stokes equation (1)
a penalty force f (x,t), which acts as a Lagrange multiplier associated to the constraint defined by the
b
boundary condition (2). The full problem (1)-(2) can then be rewritten as
∂ u=L(u)+f , ∇·u=0, for x∈Ω. (5)
t b
This system then has to be associated with the problem
find f such that f (x,t)=0 for all x∈Ω\Ω (t) and u(x,t)=v (x,t) for all x∈Ω (t). (6)
b b p p p
One can easily see that, when the solution is approximated by the Galerkin truncation associated with the
Fourier spectralmethod, these two conditions cannot be satisfied in an exactmanner: annon-zero function
thatvanishes ona full sub-domainis necessarilynon-analyticandcannotbe representedby a finite number
of Fourier modes. The idea of the method which is used here is to relax the problem of finding the penalty
force f to
b
find f such that |f |2d3x minimal and |u−v |2d3x minimal. (7)
b b p
ZΩ\Ωp(t) ZΩp(t)
All the difficulty of integrating this system is now in finding an efficient manner of solving the optimization
problem above, which involves local condition in the physical space, while spectral methods try to be as
local as possible in Fourier space. As this has to be done at every sub-time steps, the use of standard
3
minimization techniques would introduce a very large number of operations and result in a cost-inefficient
algorithm.
Efficientstrategiesconsistinguessingacceptablevalues ofthe fieldf . We here makeuse ofthe “direct-
b
forcing” method proposed by Fadlun et al. (2000). To explain this method, let us consider a simple Euler-
discretization in time of (5)
un+1 =un+∆t [L(un)+fn]. (8)
b
The idea ofthe direct-forcingmethod is to use the factthat we wantto impose un+1 =vn+1 inΩ to write
p p
from (8)
fn =−χ (x) L(un)+(un−vn+1)/∆t , (9)
b p
where χp(x) denotes the characteristic functio(cid:2)n of the particle: χp(x) =(cid:3) 1 if x ∈ Ωp(t) and 0 otherwise.
Other similar approaches consist in imposing the boundary conditions by means of a Darcy term, i.e.
f =−αχ (x)(u−v ) (see Angot et al., 1999). The particle volume is then seen as a porous media where
b p p
alargevalueoftheparameterαensurestheno-slipboundaryconditionsbutresultsinchoosingasufficiently
small value of the time step. The advantage of direct forcing is that it consists in combining the value of
the parameterα that is optimal for a givenchoiceof the time step ∆t, together with a compensationof the
overallevolution of the velocity field inside the particle due to L.
2.1. The smoothed particle
The expression(9) forf wouldbe correctifthe positionofthe gridpoints coincidedwiththe immersed
b
boundary. This is not true in the case of a spherical particle embedded in a cubic grid. An interpolation
procedureisthenneededinordertoprescribetheboundaryconditionsatthegrid-pointsnexttothesurface.
Forthis weuse a volume-fractionschemeproposedbyFadlun et al. (2000) andPasquetti et al.(2008). The
characteristic function χ describing the geometry of the particle is averagedover two grid cell:
p
1 h h h
χ¯ (x)= χ (x+y)d3y. (10)
p (2h)3 p
Z−hZ−hZ−h
This representationhastwo advantages. First,as atypicalvelocityfieldarounda particle iscontinuousbut
not differentiable at the surface, Gibbs oscillations are created in the vicinity of the boundary; replacing
the characteristic function by χ¯ reduces such oscillations. Second, this volume averaged χ¯ allows for a
p p
continuous displacement of the immersed boundary which is important in the case of a moving particle.
2.2. The pressure predictor
Another issue arises when combining a spectral with a penalty method. Although applying the force
f to the velocity field imposes the boundary conditions precisely, the divergence-freeness of u requires
b
a subsequent projection of the velocity field onto its solenoidal part which violates the formally correct
boundary conditions. A solution to this conflict is to solve the Poisson equation (4) for pressure with
modifiedboundary conditions as for example proposedby Taira and Colonius(2007)in combinationwith a
finite-difference scheme. However, this is not possible with a Fourier-spectral method as in this case (4) is
solved with generic periodic boundary conditions.
We attenuate this intrinsic problem of a Fourier-spectral approachby using a predictor for the pressure
gradientin(5),asconsideredbyBrown et al.(2001). Wewillnowexplainourprocedure. Theusualstrategy
in spectral codes without immersed boundaries consists in advancing first the velocity without taking into
account the pressure gradient in (1)
u˜n+1 =un+∆t[−(un·∇)un−ν∆un]. (11)
The intermediate velocity field u˜n+1 is compressible and is projected in a second step onto its solenoidal
part
un+1 =P(u˜n+1)≡u˜n+1−∇pn. (12)
4
100 #-3/2
p
1
0.8
10-1
u/Uc 0.6 Lu2B
0.4 10-2
0.2 predictor
no predictor
0 10-3
4 8 16 32
0 50 100 150 200 250 #
p
x/grid points
Figure1: (Left)Profileofthevelocity(normalizedtotheinflowspeed)inthedirectionofthemeanflowontheaxisofsymmetry.
The solid line corresponds to a computation with predictor pn−1 while the points correspond to a computation without the
predictor. (Right)L2 normofthevelocityinsidetheparticleasafunctionofthenumberofgridpointsonitsdiameter
ThisprojectionisequivalenttosubtractingthegradientofthepressurethatisgivenbythePoissonequation
(4). This splitting of the integration is called fractional-step method.
If the domain contains immersed boundaries with no-slip conditions, the penalty force f has to be
b
additionally applied. This can be seen as an intermediate step
u˜˜n+1 =B(u˜n+1)≡u˜n+1+fn (13)
b
leading to a compressible u˜˜n+1, which is again projected onto its solenoidal part as in (12). A problem
arises since the two operators P and B do not commute. The application of P spoils the result of B. Our
approach consists in attenuating this conflict by adding a predictor for the pressure gradient to (11). We
make of ∇pn−1, i.e. of the value of the pressure gradient at the last time step. This predictor already
includes the major part of the pressureneeded to keepthe velocity field very close to being divergence-free.
The predictor is updated during the projection step according to
pn =pn−1+Φn, where ∇2Φn =∇·u˜˜n+1/∆t. (14)
It is important to stress that we apply the operatorsP and B and also the update of the pressure predictor
during each sub-time step of the third-order temporal scheme (3).
TheimportanceofthepredictorforaFourier-spectralschemeisevidentfromFig.1wherethestreamwise
velocityprofileonalinepassingthroughtheparticleisshown. Theflowconditionscorrespondtoanumerical
wind tunnel experiment where a fixed spherical particle (v =0) is embedded in a homogeneous flow. The
p
velocity fluctuations introduced by the wake of the particle are eliminated in a zone at the right boundary
by an additional application of the penalty method. The accuracy of the upstream boundary condition
at the particle surface is strongly affected by the pressure predictor. With the predictor we achieve a no-
slip condition at the surface while the fluid significantly enters the particle when we neglect the pressure
predictor. Alsothesmoothingofthevelocityfluctuationattheexitofthewindtunnelismuchmoreefficient
with the pressure predictor.
A measure for the accuracy of the boundary conditions is given by the L -norm
2
1/2
uL2 = u(x)2dx (15)
B
(cid:18)ZB (cid:19)
of the velocity inside the fixed spherical particle. This norm is represented in Fig. 1 (Right) for a flow with
a particleReynolds number Re =20 andas afunction of the number of grid-points# alongthe diameter
p p
of the sphere. This number quantifies the resolution of the geometry of the particle by the underlying
5
homogeneous grid. The norm uL2 shows a power-law decrease with an exponent −3/2. This value can be
B
understoodbyrememberingthatafunctionf(x)inaFourierschemeisrepresentedbyaf (x), atruncated
K
series of K modes. Using a function f(x)=max(0,cos(x)) which resembles qualitatively the flow structure
in the vicinity of a particle at rest one can show that f2(x)dx∝K−3.
B K
R
2.3. The forces on the sphere
The force that is exerted by the fluid on the particle is given by the integral over the surface of the
particle of the fluid stress tensor
F = ∇·T d3x= T·dS, where T=−pI +ν(∇u+∇uT). (16)
3
ZΩp(t) Z∂Ωp(t)
This total force F can easily be computed from the penalty force f because
b
∇·T d3x= f d3x.
b
ZΩp(t) ZΩp(t)
In order to distinguish the pressure and viscous terms appearing in the stress tensor (16), our strategy
consists in computing the pressurecontributionby constructing a homogeneousgridof points at the sphere
surface. Then, we integrate the directed pressure pn, n being the normal vector field at the surface of the
particle,overthis grid. The viscouspartisobtainedbysubtractingthe pressurecontributionfromthe total
force.
3. Benchmarks
Now we are going to present several types of benchmarks of the proposed scheme. For this we turn
back to the numerical wind tunnel experiment already mentioned before. The principle flow configuration
is depicted in Fig. 2 above the critical Reynolds number where a non-stationary wake develops.
Figure2: TurbulentwakeatRe=400(volumerenderingofhigh-vorticityisosurfaces).
3.1. Drag coefficient
An important quantity of the dynamics of a solid sphere exposed to a homogeneous inflow is the drag-
coefficient
8F
D
C = , (17)
D πd2ρ U2
f c
whereF denotesthe totalforceexertedonthe particleinthe stream-wisedirection, U the inflowvelocity,
D c
dis the spherediameter, andρ the fluid massdensity (setto onein ournumericalsettings). The measured
f
coefficients are given in Fig. 3 in comparisonto the empirical formula
24
CSN = (1+0.15Re0.687) (18)
D Re p
p
proposed by Schiller and Naumann (1933) which fits well experimentally measured coefficients. For all
Reynolds numbers considered(fromRe =2 to Re=400)the DNS coefficients are in goodagreementwith
p
the fitting formula.
6
Re=2
p
Re=20
2 100 #p-3/2
P p
C
C/F 1
10 (16) SN
0.5 DNS CD
CD (16) Le Clair et al. (1970) SN)/D10-1
(16) 5 20Rep 100 400 C-CD
(32) (
1 (32) (48)
(32) (64)
Schiller-NaumDaNnSn 10-2
2 5 10 20 50 100 200 400 3 4 6 8 12 16 32
Re #
p p
Figure 3: (Left) Drag-coefficients CD as a function of the particle Reynolds number; in parenthesis is given the number of
internal grid points. Inset: ratio between the viscous contribution CF to the drag coefficient and the contribution CP from
pressure. (Right)Convergence rateasafunctionofthenumberofgridpoints#p insidethesphere.
The numbers in parenthesis in Fig. 3 (Left) denote the number of grid-points # on the diameter of
p
the sphere. The higher this number the better is the resolution of the boundary layer around the sphere.
The thickness of this layer decreases as Re−1/2 thus the higher the Reynolds number the more points #
p
are required in order to resolve the boundary layer. In this experiment the sphere is fixed and the inflow
velocity is one. The Reynolds number is varied by changing the viscosity of the fluid.
In order to check to performance of the scheme in case where the particle is moving we changed the
frame of reference of this experiment to the particle frame and thus move the particle in a flow at rest.
Also the zone where we remove the velocity fluctuations is moved according to the particle position. The
measured drag-coefficients and the velocity profile around the sphere are nearly indistinguishable from the
former wind tunnel experiments. We observe fluctuating deviations from the values given in Fig. 3 of less
than one percent for Reynolds number up to 50 and less than five percent up to 400.
TheinsetofFig.3(Left)showstheratiooftheviscousC tothepressureC contributiontothecomplete
F P
dragC =C +C . ThemeasuredratiosareinreasonableagreementwithresultsofLe-Clair et al.(1970).
D F P
Here we stress that the pressure contribution is not needed for the integrationof the particle dynamics and
might only be of interest for post-processing purposes.
Theconvergenceofthedrag-coefficientcanbequantifiedbylookingattherelativeerror(C −CSN)/CSN
D
asafunctionoftheinternalpoints# . FromFig.3(Right)weobservethatthecoefficientsconvergewithan
p
exponent of -3/2 for the two Reynolds numbers considered. This is the same rate of convergence as for the
L norm uL2 of the velocity inside the particle (see previous section) which suggests that the convergence-
2 B
rate of the proposed scheme is determined by the convergence-rate of the truncated representation of the
C0-velocity field.
3.2. Lift-coefficient
If the inflow velocity of the wind tunnel is not homogeneous but possesses a constant gradient
U =U +g(y−y ) (19)
y c 0
this linear shear flow creates a lift force F = F ·e on the particle perpendicular to the mean flow. The
L y
rate of shear is usually given in terms of a non-dimensional shear rate
gd
s= , (20)
U
c
where g is the dimensional shear in (19). Recently, the corresponding lift-coefficients
8F
L
C = , (21)
L πd2ρ U2
f c
7
where F represents the amplitude of the lift force, have been measured by three different groups (see
L
Dandy and Dwyer, 1990; Kurose and Komori, 1999; Bagchi and Balachandar, 2002).
0.1 1
1
DNS
0.08 a exp(-bx)
0.06 0.8 ΩΩ/stf
0.04
0.6
0.02 Ωf
CL 0 Ω/st 0.1 0 10 20 30 40
0.4 Re
G
-0.02
s=0.2
DNS
-0.04 DNS (double domain) 0.2
Bagchi-Balachandar 2002
-0.06 Kurose-Komori 1998 Bagchi and Balachandar 2002
Dandy-Dwyer 1989 DNS
-0.08 0
10 100 0 5 10 15 20 25 30 35 40
Re Re
p G
Figure4: (Left)Lift-coefficientsCL asafunctionoftheparticleReynoldsnumberfors=0.2. (Right)Stationaryrotationrate
(normalizedtotheambientrotationrateΩf =g/2)asafunctionoftheparticleReynoldsnumberdefinedwiththemeanfluid
velocitygradientReG=gd2/ν. Inset: Logarithmicrepresentation. Fortheexponential fitweuseda=0.8andb=0.043.
From Fig 4 (Left) one observesa quite largespreadof the lift coefficients measuredby the three groups.
Onereasonforthese differing values is the sensitivityofthe lift forceonthe stressesinthe flowandthus on
the size of the computational domain. It has been claimed by Bagchi and Balachandar (2002) that the box
size in the work of Dandy and Dwyer (1990) was too small, leading to an overestimation of the lift force.
We have computed the lift force for two different box sizes and observe a decreasing force with increasing
box-size. Our measurements are in good agreement with the two recent works and reproduce in particular
the transition to negative values of the lift force for Re ≈60−80.
p
3.3. Torque-free rotation rate
As a last benchmark we present in Fig. 4 (Right) the torque-free rotation rate in a linear shear flow
described in the precedent section. The sphere is fixed but allowed to rotate freely. After a relaxation
period we observe a rotation rate which is exponentially (see inset of Fig. 4 (Right)) decreasing with the
Reynoldsnumber. ThemeasuredvaluesagreewellwiththoseofBagchi and Balachandar(2002). Thesmall
deviations might again be a result of different box-sizes used.
4. Statistical properties of the turbulent wake at moderate Reynolds numbers
This section is dedicated to the study of turbulent velocity statistics in the wake of a sphericalobstacle.
ThedevelopmentofaFourier-spectralcodethatintegratesthefluidmotionandprescribesno-slipboundary
conditionsonthe body surfaceallowsoneto achieveprecisedirectnumericalsimulationsofturbulentwakes
at a reasonable computational time.
4.1. Mean flow and wake similarity
Most models for turbulent flows past an obstacle rely on studying solutions to the Reynolds-averaged
equations (see, e.g., Pope, 2000). This approach consists in writing the velocity field u(x,t) as a sum of
its stationary average U(x) = hu(x,t)i and turbulent fluctuations u′(x,t). The mean velocity is then a
solution to
(U ·∇)U =−∇·τ −∇hpi+ν∇2U, (22)
8
where τ designates the Reynolds stress tensor τ (x)=hu′(x,t)u′(x,t)i. In order to obtain the properties
ij i j
ofthe meanvelocityprofileU, Eq.(22) is generallyclosedby means ofthe eddy viscosityassumption. This
hypotheses consists in relating linearly the anisotropic part of the Reynolds stress to the mean rate, that is
2
τ (x)− kδ =ν (x)[∂ U +∂ U ], (23)
ij ij T i j j i
3
wherek ≡(1/2)traceτ =(1/2)hu′u′iisthe turbulentkinetic energyandthe positivequantityν (x)isthe
i i T
eddyviscosity. Whentheeddyviscosityν isassumedtobeconstantinspace,andinparticularindependent
T
ofthestreamwisedistancexfromtheparticlecenter,theReynolds-averagedequation(22)canbeintegrated
to obtain the following self-similar estimate for the streamwise mean velocity U(x,y,z) = hu (x,t)i (see,
x
e.g., Schlichting, 1979)
C d2 y2+z2
D
U(x,y,z)=U 1− exp − , (24)
c 32ν (x−x ) r2 (x)
" T 0 1/2 !#
where the wake half-width r is given by
1/2
r (x)=d 2U (x−x )/(ν ln2). (25)
1/2 c 0 T
p
Such formula, which are very commonly used in industrial applications, have two important consequences.
First, the mean velocity deficit ∆U = 1−U(x,0,0)/U should behave at large distances in the turbulent
c
wake as x−1. Secondly, the normalized velocity deficit profile f(x,y,z)=[U −U(x,y,z)]/[U −U(x,0,0)]
c c
should collapse to the similar profile f(x,y,z)=2−ξ2 where ξ =r/r (x) and r= y2+z2.
1/2
An analysis of the turbulent wake from a direct simulation performed at the moderate value of the
p
particle Reynolds number Re =400 shows that this similar eddy-viscosity theory is able to rightly predict
p
thedecayofthemeanvelocitydeficit. Figure5(Left)representsthebehaviorof∆U asafunctionofx/d−1.
Its displays a power-lawdecay with a exponent −1 over nearly half a decade. The value −1 is undoubtedly
obtained as seen from the inset of Fig. 5, which shows the value of the local slope (logarithmic derivative)
of the velocity deficit as a function of the distance in the wake. The oscillations around the value −1 are a
signature of the Gibbs phenomenon induced by the strip forcing at x≈12d.
1
x = d
100 0,0)] 0.8 xx == 23dd
1 − U (x,0,0)/Uc10−1 Local Slope−−01210 2 4 6 8 slope = −1 U−U (x,y,z)]/[U−U (x,cc 000...246 xxxxxf ====== 245678−dddddξ2
x / d − 1 = [ 0
f
10−1 x /1 d00 − 1 101 −0.20 0.5 1 1.5 ξ = r2 / r1/2 2.5 3 3.5 4
Figure 5: (Left) mean velocity deficit ∆U = 1−U(x,0,0)/Uc in the axis of symmetry y = 0 and z = 0 of the wake as
a function of the normalized streamwise distance x/d. As seen from the inset, which represents its logarithmic derivative
[dln∆U]/[dln(x/d−1)],itdisplaysapower-lawbehavior∝1/(x−x0)withx0≃d. (Right)meanvelocitydeficitprofilef as
afunctionofthereducedspanwisedistanceξ=r/r1/2 forvariousvaluesofx/daslabeled;thebolddashedlinerepresentsthe
Gaussianfunctionalformoff obtainedfromaconstant-eddy-viscosity similarprofile.
However, as seen from the right-hand panel of Fig. 5, the mean velocity deficit profile does not seem
to display any similarity. The half width r was estimated here as the value of the distance to the
1/2
axisymmetric axis of the wake such that the velocity deficit is exactly equal to 1/2. By definition, the
various curves should go through f = 1 at ξ = 0 and f = 1/2 at ξ = 1. Despite these two constraints,
9
it is seen that they do not collapse on the top of each other for values of ξ close to 1/2 and 3/2. The
discrepancy is even stronger for ξ >2, where f tends to negative values that depend on x. These deviation
at large ξ’s might be due to the periodicity of the spatial domain in the y and z directions that results
in a mean streamwise velocity of the fluid that is larger than U at y = ±π and z = ±π. However, the
c
deviation from the similar profile at moderate values of ξ are comparable to those obtained experimentally
by Uberoi and Freymuth (1970) and Wu and Faeth (1994). In particular, a striking observation is that in
allcaseswithmoderatelylargeReynoldsnumber,the meanvelocitydefectprofileseemstovanishforξ ≈2.
Such a feature is clearly reproducedin our simulation, as seen from Fig. 5 (Right) where the various curves
are pinched together before reaching negative values.
4.2. Two-point statistics and turbulence decay
We have seen above that the self-similar wake approach reproduces well the decay ∝ x−1 of the mean
velocity deficit. However, this approach cannot be fully true as no similarity is observed in the wake. In
order to get a better handling, we relate in this section this observed behavior to the decay of turbulent
kinetic energy in the wake of the spherical particle.
It is known since the seminal work of Kolmogorov (1941) that the law of decay of turbulent kinetic
energy is dictated by the “principle of permanence of large eddies”. We follow here the phenomenological
presentationof this law made by Frisch (1995). If we assume that the initial turbulent velocity correlations
follow a infrared asymptotic self-similarity of the form hu(x+r)u(x)i ∼ r−h at large scales r, the decay
behavior of all turbulent quantities is dictated by the value of the exponent h. Indeed, by continuity of
velocity statistics, the integral-scale velocity should then itself behave as U = k1/2 ∼ L−h, where k is the
turbulent kinetic energy and L is the integral scale. However, we know that the energy decays with a rate
that can be related to velocity by dk/dt = −ε ∼ −k3/2/L. Replacing k by L−2h in this equation, we
obtain that dL/dt∼L−h, so that L∼t1/(1+h). This evolution of the integral scale implies that the kinetic
turbulent energy behave as k ∼t−2h/(1+h).
Whenconsideringthedecayofwaketurbulence,alltheturbulentquantitieshavetobeofcourseevaluated
in planes perpendicular to the axisymmetric axis and time has to be replaced by x/U by invoking Taylor
c
hypothesis. Theaboveconsiderationshence implythatifthe velocitycorrelationshu(r)u(0)ibehaveasr−h
in a plane x = x located in the developing wake, then at any x > x , the large-scale correlations should
0 0
rememberthisinitialformandtheturbulentkineticenergyisexpectedtodecayask ∼x−2h/(1+h). Figure6
representsthe correlationofthe radialcomponentof the fluctuating velocityat variousvalues ofx. As seen
0.1
100
0)|
> C (
u(x,0) r 0.05 C (r) / 10−2
r) |
x, 0 0.5 1 1.5 2
u(r r / d
<
= 0
r)
C (
x = 2d x = 4d x = 6d x = 8d x = 10d
−0.05
0 0.5 1 1.5 2
r / d
Figure 6: Correlation of the velocity radial fluctuation as a function of the distance r to the axisymmetric axis for various
valuesofthestreamwisedistancexaslabeled. Inset: sameinsemi-logarithmiccoordinatesshowingthatthevariouscorrelation
fallonthetopofeachother andfollowanexponential decay.
from this figure, velocity correlations collapse at large scales on the top of each other onto an exponential
10