Inertial particle acceleration in strained turbulence
Chung-min Lee, Ármann Gylfason, Prasad Perlekar, Federico Toschi
aa r X i v : . [ phy s i c s . f l u - dyn ] A p r Inertial particle acceleration in strainedturbulence
Chung-min Lee , ´Armann Gylfason ∗ , Prasad Perlekar , and Federico Toschi Department of Mathematics and Statistics, California State University Long Beach, USA School of Science and Engineering, Reykjav´ık University, Menntavegur 1, 101 Iceland TIFR Centre for Interdisciplinary Sciences, Hyderabad 500075, India Department of Applied Physics, Eindhoven University of Technology, The Netherlands Department of Mathematics and Computer Science, Eindhoven University of Technology, The Netherlands Istituto per le Applicazioni del Calcolo CNR, Via dei Taurini 19, 00185 Rome, Italy
Abstract
The dynamics of inertial particles in turbulence is modelled and investigated by meansof direct numerical simulation of an axisymmetrically expanding homogeneous turbulentstrained flow. This flow can mimic the dynamics of particles close to stagnation points.The influence of mean straining flow is explored by varying the dimensionless strain rateparameter Sk / ǫ from 0.2 to 20. We report results relative to the acceleration variancesand probability density functions for both passive and inertial particles. A high mean strainis found to have a significant effect on the acceleration variance both directly, through anincrease in wave number magnitude, and indirectly, through the coupling of the fluctuat-ing velocity and the mean flow field. The influence of the strain on normalized particleacceleration pdfs is more subtle. For the case of passive particle we can approximate theacceleration variance with the aid of rapid distortion theory and obtain good agreementwith simulation data. For the case of inertial particles we can write a formal expressionsfor the accelerations. The magnitude changes in the inertial particle acceleration varianceand the effect on the probability density function are then discussed in a wider context forcomparable flows, where the effects of the mean flow geometry and of the anisotropy atthe small scales are present. The motions of small, passive and inertial particles in turbulence has been extensively studiedin recent years, both from the experimental and theoretical viewpoints. This was motivated bya broad range of applications such as the spread of pollutants in the atmosphere and oceans,the process of rain and ice formation in cloud, the transport of sediments in rivers and estu-aries, see e.g. the reviews of Shaw (2003) and Toschi & Bodenschatz (2009). Progress in theunderstanding was on the one hand made possible by recent improvements in Lagrangian mea-surements through particle tracking methodologies, resulting in part from rapid advances inhigh speed imaging (Virant & Dracos (1997); Ott & Mann (2000); Voth et al. (2002); Xu et al. ∗ Email address for correspondence: [email protected] et al. (2006); Ayyalasomayajula & Warhaft (2006); Gualtieri & Meneveau (2010)). A mean strainingflow naturally appears in the proximity of stagnation points. Flow straining is furthermore offundamental interest since it induces a scale dependent anisotropy; the smallest scales of the flowmay be nearly isotropic, whereas the largest scales are highly anisotropic (Biferale & Procaccia(2005)).Furthermore, many flows naturally combine straining geometries and inertial particles. Theflow geometry presented here, namely particle laden turbulent flow undergoing an axisymmetricexpansion, has similarities with combustor diffusers in jet engines (Klein (1995)), where liquidfuel is injected in an expanding flow, and with the flow in the combustion chamber in an internalcombustion engine during the compression stroke of the fuel air mixture (Han & Reitz (1995)).While significant attention has been given to the study of Lagrangian acceleration statisticsin isotropic turbulence, less attention has been paid to the implications of anisotropic largescale flow geometry on the Lagrangian dynamics. Recent experimental and numerical work onthe Lagrangian behavior of inertial particles in shear flows and in turbulent boundary layerhas shown pronounced effects on the inertial particle statistics (Gerashchenco et al. (2008);Lavezzo et al. (2010); Gualtieri et al. (2009, 2012)). The persistent small scale anisotropy hasbeen found to influence the geometry and alignment of particle clusters and relative particle pairvelocities. In addition, the combined effects of gravity and shear on particle acceleration varianceresults in an increase in magnitude with the Stokes number. As a consequence, the accelerationprobability distribution functions (pdfs) became increasingly narrow and skewed with inertia.Here, we address a related topic, namely the complexity introduced in the Lagrangian dynamicsof tracer and inertial particles due to flow straining.In an effort to realize the effects of anisotropy in the particle dynamics, we numericallysimulate axisymmetric expansion of initially isotropic turbulence. The flow is seeded withinfinitesimal tracer and inertial particles of varied Stokes numbers. We measure the particlevelocity and acceleration statistics, including variances and probability density functions fordifferent strain rates and Stokes numbers. Comparisons are made with predictions of rapiddistortion theory on tracer accelerations, and the solutions of the Stokes equations for inertialparticles in the straining flow.The paper is organized as follows. In section 2 we briefly introduce the numerical methods forsimulating an axisymmetric turbulence and particle movements. Parameters of the simulationsare listed. Section 3 presents the underlying flow field. We discuss our main findings of particleacceleration variances and probability density functions in simulation data with the support oftheoretical estimations in section 4. We also discuss our result in the context of previous workin shear flows. In section 5 we present our conclusions.2 zy Figure 1: A sketch of the deformation of the simulation domain under straining. The mean flow U = (− Sx, Sy, Sz ) corresponds to an ideal flow onto a flat plate. The deforming domain isinitially elongated in the x -direction but becomes wider in y and z directions with time. Arrowsindicate the directions of the stream lines of the induced mean flow. The equations describing the motion of an incompressible Newtonian fluid are the continuityequation and the Navier-Stokes equations, respectively: ∇ ⋅ ˜ u = ,∂ ˜ u ∂t + ˜ u ⋅ ∇ ˜ u + ∇ ˜ p = ν ∇ ˜ u . (1)Here, ˜ u and ˜ p are the instantaneous flow velocity and pressure, respectively, and ν is thekinematic viscosity of the fluid.In this paper we are concerned with a turbulent fluid undergoing an axisymmetric expansion,where the mean flow field is described by U = ( − Sx, Sy, Sz ) , (2)here S is the mean strain rate S = √ ( ¯ S ij ¯ S ij ) / , and ¯ S ij = ( ∂U i ∂x j + ∂U j ∂x i ) is the mean rateof strain tensor. The mean flow corresponds to an ideal flow onto a flat plate, and is realizedin the flow between contracting pistons or in the expanding flow through a diffuser. Figure 1shows a sketch of the mean flow field, namely streamlines of the mean field, and the deformingdomain.Applying the Reynolds decomposition, and expressing (1) in terms of the vector potential b , with u = ∇ × b , one obtains − ∂ t ∇ b − ∇ × ( u × ω ) + Sx∂ x ∇ b − Sy∂ y ∇ b − Sz∂ z ∇ b + S ∇ b − S ∇ b ˆ e = − ν ∇ b , (3)where u = ˜ u − U is the velocity fluctuation, ω is the vorticity, defined as ω ≡ ∇ × u = −∇ b . b is the first component of b and ˆ e is the unit vector in the x -direction. In the following webriefly outline the numerical algorithm. 3s in Gylfason et al. (2011), in order to solve (3) we apply a pseudo-spectral method withRogallo’s algorithm (Rogallo (1981)) where the following variable transformations are per-formed: x ′ = e St x, y ′ = e − St y, z ′ = e − St z, t ′ = t, and hence the vector potential of the velocity fluctuation satisfies ∂ t ′ ∇ ′ b − ∇ ′ × ( u × ω ) + S ∇ ′ b − S ∇ ′ b ˆ e = − ν ∇ ′ b , (4)where ∇ ′ = ( e St ∂ x ′ , e − St ∂ y ′ , e − St ∂ z ′ ) . By adopting this new coordinate system, the physicaldomain deforms with time while the computational lattice grid is time independent, and theflow equations become periodic. We then apply the pseudo-spectral method to equations in(4). More details about the numerical method can be found in Gylfason et al. (2011).Numerical simulations of this axisymmetric expansion flow were carried out on a Cartesiangrid with 1024 × ×
256 and 2048 × ×
512 grid points in the x , y and z directions, respectively.The initial configurations are derived from statistically independent homogeneous and isotropicflow simulations which have reached a stationary state after more than 5 large-eddy turnovertimes. The Reynolds numbers, based on the Taylor microscale λ were R λ =
117 and R λ = [ , π ] × [ , π ] × [ , π ] in x , y and z directions respectively, andthe simulation was terminated when the domain has reached [ , . π ] × [ , . π ] × [ , . π ] toprevent the physical domain to become too flattened.The top of Figure 2 shows snapshots of the fluctuating velocity magnitude at three timeinstants during the straining. From left to right, the non-dimensional times are S × t = . S × t = .
64, and at S × t = .
96 (just before thestrain simulation is terminated due the large deformation of the physical domain). Additionally,the figure shows the coordinate system adopted in the text, and the geometry of the simulationdomain selected and its deformation. Production of turbulence overwhelms dissipation duringthe straining, reflected in an increase in the turbulent kinetic energy, most notably in thecompressed component ( x ). This can be seen in the warmer colors in the rightmost plot. Thebottom of Figure 2 shows isosurfaces of non-dimensional vorticity ω /( ε / ν ) / = .
17 at the sametime instants as above, which is respectively 4.36, 2.86 and 2.43 standard deviations above themean vorticity magnitude at the three time instants. From the increased number and size ofthe filaments, we observe that the vorticity is intensified during straining, and the filaments arefound to gradually align with the y, z -plane due to the mean flow extension in the plane.
To study Lagrangian aspects of this flow we seed the flow with tracers and inertial particles.Here, we are concerned with particles that are small compared to the smallest length scalespresent in the flow and their densities considerably higher than the fluid density. The particlenumber densities are furthermore assumed to be sufficiently low so that particle-particle inter-actions can be ignored. Under the above approximations, the coupling of the particles with thecarrier fluid can be ignored.The Lagrangian equations of inertial particle motion is derived from Newton’s second law,and represents the balance between the forces acting on the particles (inertia and Stokes drag).The equations describing the motion of a particle of diameter d p and density ρ p , located at x p ω = ∣∣ ∇ × u ∣∣ (bottom) in the deforming domain in a realization ofthe axisymmetrically expanding flow. Left to right: from the onset of the straining simulationto the end of the straining simulation at time instants S × t = . , .
64, and 0 .
96. The size ofthe simulation domain is noted in the figure, and its deformation is displayed. The coordinatesystem adopted in the text is also shown. We simulated the axisymmetric turbulence with tworesolutions using 2048 × ×
512 and 1024 × ×
256 computational nodes. The illustrationshows one of the realizations of the 1024 × ×
256 simulations with S =
10. The isosurfacesplotted have non-dimensional vorticity ω /( ε / ν ) / = .
17, which is respectively 4.36, 2.86 and2.43 standard deviations above the mean vorticity magnitude at the three time instants andchosen to illustrate the flow structure. 5nd with instantaneous velocity ˜ v p are (Maxey & Riley (1983); Bec et al. (2006)): d x p dt = ˜ v p , (5) d ˜ v p dt = τ p ( ˜ u ( x p ) − ˜ v p ) , (6)where τ p = βd p / ν is the Stokes relaxation time for the particle and β = ( ρ p − ρ f )/ ρ f isthe relative density ratio between the particle and the fluid. The Stokes number St = τ p / τ η characterizes the inertia of a particle in the flow, where τ η is the Kolmogorov timescale of theflow.For the tracer particles (zero inertia), the particle velocity is the same as the fluid velocityat the particle location, which yields the kinematic relation: d x p dt = ˜ u ( x p ) . (7)The ordinary differential equations (5), (6), and (7) are solved numerically by the secondorder Adams-Bashforth method. In equations (6) and (7), the instantaneous flow velocity atthe particle location, x p , is evaluated as˜ u ( x p ) = U ( x p ) + u ( x p ) ;that is, the mean flow velocity is evaluated at the location of the particle through the formula U ( x p ) = ( − Sx p , Sy p , Sz p ) , and the flow velocity fluctuation is interpolated to the particleposition.We initialize the particles and the fluid velocity with steady state homogeneous isotropicsimulations. The particles are uniformly distributed over the domain prior to the forced homo-geneous isotropic simulation is carried out. At the beginning of the straining, at t =
0, we addthe mean flow velocity and acceleration component due to the strain geometry to the existingparticle velocity. We conduct simulations with 1024 × and 2048 × collocation points.For the lower resolution simulation we use 16 independent flow realizations with 5 × particlesof each type (6 different Stokes numbers) and for the higher resolution simulations we perform10 independent flow realizations with 4 × particles of each type.Table 1 shows the various flow parameters for the simulations performed. The range of strainrates selected is such that its effect on the smallest scales of the flow range from being negligibleto substantial. The higher strain rates are felt intensely by the large scale flow, whereas thelower strain rates have mild effect on the large scales. The value of the strain parameters, Sτ η and Sk / ǫ , which compare the strain time with the local timescales of the flow, indicate theimportance of the various terms in the evolution equation of the velocity field. The left hand side of Figure 3 shows the evolution of the Reynolds stresses ( ⟨ u i ⟩ ) normalizedby the initial turbulent kinetic energy ( k ). The component along the compressed direction( x ) grows rapidly in most cases, whereas the components along the expanding directions aresuppressed or remain roughly constant. At the lowest strain rate all component are suppressedduring the simulation time. For large times, rapid distortion theory (RDT) predicts an expo-nential growth of the Reynolds stresses, in the proportions ⟨ u ⟩ = ⟨ u ⟩ = ⟨ u ⟩ . The kinetic6able 1: Flow parameters in the Direct Numerical Simulations, based on the homogeneousisotropic simulation prior to the application of the straining. Here k ≡ (⟨ u ⟩ + ⟨ u ⟩ + ⟨ u ⟩) is the turbulent kinetic energy, ǫ ≡ ⟨ ∂u i ∂x j ∂u i ∂x j ⟩ is the energy dissipation rate, ℓ ≡ ( u rms ) / ǫ is the integral length scale, and the Kolmogorov length scale η ≡ ( ν / ǫ ) / . The subscript 0indicates the parameter values are taken prior to the straining. (Units are arbitrary)Simulation domain 1024 × ×
256 2048 × × R λ ≡ u rms λ / ν
117 193 k ≡ (⟨ u ⟩ + ⟨ u ⟩ + ⟨ u ⟩) λ / η ≡ u rms ( ν / ǫ ) / / η S . , . , , ,
10 1 , , τ η ≡ ( ν / ǫ ) / u rms = ( k / ) / ℓ / η ≡ ( u rms ) / ǫ / η η ≡ ( ν / ǫ ) / ± ± ǫ ≡ ⟨ ∂u i ∂x j ∂u i ∂x j ⟩ ± ± ν St = τ p / τ η , . , . , . , , , . , . , . , . , . Sτ η . , . , . , . , .
51 0 . , . , . Sk / ǫ . , . , . , . , . . , . , . b ij = ⟨ u i u j ⟩⟨ u i u i ⟩ − δ ij with time. The curves corresponding to thelowest strain rates markedly deviate from the others as the turbulent kinetic energy decreasesduring the straining, and the straining motions are fairly mild, even for the largest scales ofmotions.The short term RDT prediction plotted in Figure 3 is derived from the Reynolds stressequation (Pope (2000)) ddt ⟨ u i u j ⟩ = P ij + R ( r ) ij , (8)where P ij ≡ − ⟨ u i u k ⟩ ∂U j ∂x k − ⟨ u j u k ⟩ ∂U i ∂x k is the production rate of Reynolds stress, R ( r ) ij ≡ ⟨ p ( r ) ρ ( ∂u i ∂x j + ∂u j ∂x i )⟩ is the rapid pressure-rate-of-strain tensor, and p ( r ) is the rapid pressure that satisfies ρ ∇ p ( r ) = − ∂U i ∂x j ∂u j ∂x i . Right before the straining starts, the initial configuration of flow is isotropic, and R ( r ) ij = − P ij (Pope (2000)). In an axisymmetric expansion flow, the production rates are P = S ⟨( u ) ⟩ , P = − S ⟨( u ) ⟩ , and P = − S ⟨( u ) ⟩ . Therefore, in early times RDTpredicts ⟨( u ) ⟩ = ⟨( u ) ⟩ e St , ⟨( u ) ⟩ = ⟨( u ) ⟩ e − St , ⟨( u ) ⟩ = ⟨( u ) ⟩ e − St , (9)where ⟨( u i ) ⟩ represent the initial values of the Reynolds stresses ( ⟨( u i ) ⟩ , i = , , Sτ η ≫ Sk / ǫ ≫ × t h u i i / k -1 S × t b ii -0.3-0.2-0.100.10.20.30.40.50.6 Figure 3: Left: Normalized Reynolds stresses ⟨ u i ⟩/ k vs. strain time S × t . Right: The flowanisotropy tensor b ij = ⟨ u i u j ⟩⟨ u i u i ⟩ − δ ij vs. strain time. Solid symbols represent the i = i = k is the initial kinetic energy of the fluidprior to the straining. The diamond ( ◇ ), pentagram ( ☆ ), circle ( ◯ ), square ( ◻ ) and triangle( △ ) symbols represent data from S = . , . , , ,
10 in the R λ =
117 flow, respectively. Blue,green and red lines, indicate the data from S = , ,
10 in the R λ =
193 flow; solid and dashedlines represent the i = i = S = . , ,
10 areshown and are computed according to (10) with X j being ⟨ u i ⟩/ k and b ii in the j -th realization.8 / τ η -1 h ( ∂ x u ) i / h ( ∂ x u ) i S = 1 , R λ = 117 S = 4 , R λ = 117 S = 10 , R λ = 117 S = 1 , R λ = 193 S = 4 , R λ = 193 S = 10 , R λ = 193HIT Figure 4: Ratios of the variances of velocity derivatives ⟨( ∂ x u ) ⟩/⟨( ∂ x u ) ⟩ vs. normalized time t / τ η . The circle ( ◯ ), square ( ◻ ) and triangle ( △ ) symbols represent strain rates S = , ,
10 fromthe R λ =
117 set, and crosses ( × ), asterisks ( ∗ ) and pluses ( + ) represent strain rates S = , , R λ =
193 set. An estimate of statistical error bar is shown (within the symbols) andcomputed according to (10) with X j being ⟨( ∂ x u ) ⟩/⟨( ∂ x u ) ⟩ in the j -th realization. Thesolid line shows the theoretical prediction for this ratio in the isotropic turbulence.(Batchelor (1953)). Only at the highest rate of strain is the latter constraint weakly satisfied,and therefore one does not observe close matches between the predictions of RDT and theReynolds stresses in simulation data. The global anisotropy is much less sensitive, and shortterm RDT predicts the anisotropy well.The error bars in Figure 3 indicate the statistical error of the quantities estimated from thefinite number of realizations of the flow. That is, in N realizations of the turbulent flow with aparticular strain rate S , one obtains samples { X , X , . . . , X N } of a quantity X . The estimatedstandard error is the sample standard deviation of { X , X , . . . , X N } divided by √ N . That is, √∑ Nj = ( X j − X ) √ N − √ N (10)for data from N realizations; here X is the mean of { X j } Nj = . In this work, the length of thesymmetric error bars in Figures 3 to 8 is twice of the estimated statistical error.Since tracer and inertial particle accelerations occur primarily at the smallest scales ofmotions, it is useful to look at the effects of the straining on the small scales. Figure 4 showsa measure of the small scale anisotropy, the ratio of variances of longitudinal derivatives of thetransverse and longitudinal velocity components ⟨( ∂ x u ) ⟩/⟨( ∂ x u ) ⟩ with respect to time. Atthe highest strain rates, the anisotropy due to the straining is present at the smallest scales ofmotion, whereas for the lower strain rates, the flow appears nearly isotropic at the small scales.Note that the isotropic prediction for the ratio is ⟨( ∂ x u ) ⟩/⟨( ∂ x u ) ⟩ =
2. The small scaleanisotropy appears to become close to a constant after an initial transition period for the lowerstrain rates, but a stationary state is not reached at the highest rate of strain for this quantity.The left hand side of Figure 5 shows the time evolution of longitudinal derivative skew-9 / τ η -1 -0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.100.10.20.3 x,S = 1 ,R λ = 117 y,S = 1 ,R λ = 117 x,S = 4 ,R λ = 117 y,S = 4 ,R λ = 117 x,S = 10 ,R λ = 117 y,S = 10 ,R λ = 117 x,S = 1 ,R λ = 193 y,S = 1 ,R λ = 193 x,S = 4 ,R λ = 193 y,S = 4 ,R λ = 193 x,S = 10 ,R λ = 193 y,S = 10 ,R λ = 193 t/ τ η -1 Figure 5: Left: Skewness in the compression and expansion directions. Solid lines (—): ⟨( ∂ x u ) ⟩/⟨( ∂ x u ) ⟩ / . Dashed lines ( −− ): ⟨( ∂ y u ) ⟩/⟨( ∂ y u ) ⟩ / . Right: Kurtosis in the com-pression and expansion directions. Solid lines (—): ⟨( ∂ x u ) ⟩/⟨( ∂ x u ) ⟩ . Dashed lines ( −− ): ⟨( ∂ y u ) ⟩/⟨( ∂ y u ) ⟩ . In both plots: circles ( ◯ ), squares ( ◻ ) and triangles ( △ ): S = , , R λ =
117 set, and crosses ( × ), asterisks ( ∗ ) and pluses ( + ): S = , ,
10 from the R λ =
193 set. An estimate of statistical error bar is shown and computed according to (10)with X j being the skewness and kurtosis in the j -th realization.ness, along the directions of compression and expansion in the flow. Before the straining isapplied, the skewness has a value of about [ − . , − . ] as expected, but the straining causesa marked change in its value and becomes positive in the expanding direction. The effect inthe compressed direction is more subtle, but an increase in magnitude (larger negative values)appear to occur for all strain rates given that the simulation is run for a sufficiently long time.The sign change indicates a change in the small scale structure of turbulence, namely thatvortex structures dominate sheet-like structures resulting in an inhibition of the energy cascade(Ayyalasomayajula & Warhaft (2006)).The right hand side of Figure 5 shows the longitudinal kurtosis, a measure of the flow inter-mittency. Here the effects are milder in both directions, although a small increase in the kurto-sis is noted in the expanding directions (from the expected value of about 5-6, Gylfason et al. (2004)). The non-uniform mean velocity field has a significant effect on the dynamics of tracers andinertial particles, both directly via the mean flow velocity and indirectly through the strainedturbulent field. 10 .1 Tracers accelerations in straining flow
The full acceleration of tracer particles in our flow is given by the equation: ̃ a p i = D ˜ u p i Dt = Du p i Dt + u p j ∂U p i ∂x j + U p j ∂U p i ∂x j , i = , , . (11)The subscript p indicates that the full instantaneous flow velocity ˜ u p i , the fluctuating flowvelocity u p j , and the mean flow velocity U p i are taken at the location of the tracer. Here, wehave assumed that the mean flow is time-independent. The mean of the tracer accelerationis U p j ∂U pi ∂x j ; equal to the acceleration in a laminar flow of the same strain geometry. We areinterested in the statistics of the tracer acceleration fluctuation, that is, when the pure meanflow contribution to the acceleration has been subtracted: a p i = Du p i Dt + u p j ∂U p i ∂x j , i = , , . (12)The first term on the right hand side refers to the material derivative of the fluctuating velocityfield and represents the acceleration experienced by the fluid particle advected by the fluctuatingflow field, and the second term refers to the tracer acceleration induced by turbulent transportin the mean velocity field.In the variance of acceleration, the cross terms give rise to ⟨ Du pi Dt u p i ⟩ (no sum over i ),representing the time derivative of the kinetic energy in the i -th component of velocity. In ad-dition, the latter term describes contributions of velocity variances, notably ( S ) ⟨( u p ) ⟩ and S ⟨( u p ) ⟩ for the first and second component, respectively. These terms are easily evaluatedif the statistics of flow velocity fluctuations are available. When the straining is sufficiently rapid, the nonlinear and viscous forces vanish from the Navier-Stokes equations (e.g. see Pope (2000)), and therefore their solution is particularly convenient incomparison to solving the full Navier-Stokes equations. Below, we re-derive the RDT predictionsfor the evolution of the fluctuating velocity variances, as well as deriving the prediction of RDTon the evolution of the tracer acceleration variance.In RDT, each Fourier mode evolves independently. Let us consider a single mode of thefluctuating velocity u ( x , t ) = ˆ u κ ( t ) e i κ ( t ) ⋅ x The wave number and the Fourier coefficients evolve according to the following set of equations(Pope (2000)) dκ ℓ dt = − κ j ∂U j ∂x ℓ , (13) d ˆ u j dt = − ˆ u k ∂U ℓ ∂x k ( δ jℓ − κ j κ ℓ ∣ κ ∣ ) . (14)When the mean flow geometry, U = ( − Sx, Sy, Sz ) , has been applied, equation (13), results in κ ℓ ( t ) = κ ℓ e − S ℓ t , (15)11here S = − S, S = S = S , and ∣ κ ∣ = ( κ ) e St + ( κ ) e − St + ( κ ) e − St . Similarly, for ˆ u j ,equation (14), gives, d ˆ u j dt = − S ℓ ˆ u ℓ ( δ jℓ − κ j κ ℓ ∣ κ ∣ ) . In particular, the long term prediction givesˆ u ( t ) ≈ ˆ u e − St , ˆ u ( t ) ≈ ˆ u e − St , ˆ u ( t ) ≈ ˆ u e − St . Taking the material derivative of the single mode flow fluctuating velocity, we have Du Dt ≈ − Su , Du Dt ≈ − Su , Du Dt ≈ − Su . (16)Therefore the anisotropic contribution to ⟨( Du i / Dt ) ⟩ can be estimated as S i ⟨( u i ) ⟩ , where S = − S, S = S = S are the strain rate at different directions. Together with the normalizedacceleration variances a ≡ ( / )⟨ a i a i ⟩/(⟨ ǫ ⟩ / ν ) / at the beginning of the straining (Voth et al. (2002)), the magnitude of tracer fluctuating acceleration variance can be approximated as ⟨( a p i ) ⟩ ≈ a ( ⟨ ǫ ⟩ ν ) / + S i ⟨( u i ) ⟩ + S i d ⟨( u i ) ⟩ dt , i = , , , (17)where the first term represents the isotropic contribution to acceleration variances with depen-dence on turbulence intensity, and all the terms on the right-hand side are Eulerian quantitiesof the straining flow. Here ⟨ ǫ ⟩ is the mean energy dissipation rate in homogeneous isotropicturbulence across all realizations, and ⟨ ǫ ⟩ is the time-dependent mean energy dissipation ratein the straining turbulence across all realizations.Figure 6 shows the acceleration variances of passive tracers in the straining flow. Thevariance is higher in the compression direction than in the expanding directions due to the meanstraining geometry, and the effects are seen immediately after the strain has been imposed. Theapproximations derived above, applying RDT, fit nicely to the simulation data as shown. Thechange in the isotropic dissipation rate due to straining is accounted for in the first term of(17), and the rest of the terms involve the mean strain and the velocity fluctuations. The meanstrain causes a rapid increase in the acceleration variance as the strain is applied, particularlyat the higher rate of strain. The subsequent increase and the differences between the individualcomponents are partially due to the evolution of the velocity variances for each component;namely that the compressed velocity components are emphasized (increasing energy content)whereas the expanding velocity components are either suppressed or maintained. The differential equations for the particle position as a function of time, obtained by combiningthe mean flow components in equations (2) with the equations (5) and (6) are second orderlinear ordinary differential equations with constant coefficients: d x p i dt + τ p dx p i dt − S i τ p x p i = τ p u i , i = , , . The roots λ , = − ± √ + S i τ p τ p of the characteristic equations of these equations prescribethe behavior of particle movements in absence of turbulence. In the compression direction, the12 × t h a p i / p ǫ / ν , h a p i / p ǫ / ν -2 -1 S =0 . ,xS =0 . ,yS =1 ,xS =1 ,yS =4 ,xS =4 ,yS =10 ,xS =10 ,yS =0 . ,a ( h ǫ i / ν ) / (4 . ,i =1(4 . ,i =2 S × t h a p i / p ǫ / ν , h a p i / p ǫ / ν S = 4 ,x,R λ = 117 S = 4 ,y,R λ = 117 S = 4 ,x,R λ = 193 S = 4 ,y,R λ = 193(4 . ,i = 1(4 . ,i = 2 Figure 6: Normalized acceleration variances of passive tracers in the compressed and expandingdirections vs. strain time at various strain rates. Solid and empty symbols represent thecompressed and expanding directions, respectively. Solid (—) and dashed ( −− ) lines indicateRDT long term predictions of the tracer acceleration variances from the straining together withnormalized flow acceleration variances in HIT shown in (17). Dash-dot line ( − . ) indicates theterm a / (⟨ ǫ ⟩ / ν ) / for the S = . ◇ ), circle ( ◯ ), square ( ◻ ), triangle ( △ )mark normalized tracer acceleration variances in strain rates S = . , , R λ = S = R λ =
117 (squares ( ∎ , ◻ ): compressed and expanding directions, respectively)and 193 (asterisks ( ∗ ) and pluses ( + ): compressed and expanding directions, respectively).An estimate of statistical error bar is shown and computed according to (10) with X j being ⟨ a p i ⟩/√ ǫ / ν in the j -th realization. 13ombination of strain rates S and the Stokes relaxation time τ p in our simulations give rise tothe discriminant D = − Sτ p , which separates two possible movements: a overdamped-decaymotion toward the stagnation plane x = D , = + Sτ p is always positive, andparticles move away from the axis y = z = D > ∶ ̃ a p ( t ) = λ − λ [( λ ( λ + S ) e λ t − λ ( λ + S ) e λ t ) x p + ( λ e λ t − λ e λ t ) v p + τ p ∫ t u ( x p ( τ ) , τ )( λ e λ ( t − τ ) − λ e λ ( t − τ ) ) dτ ] + τ p u ( x p ( t ) , t ) . (18) D < ∶ ̃ a p ( t ) = e − t τp [ − τ p v p cos ( ωt ) + τ p ω (( − Sτ p ) v p + S τ p x p ) sin ( ωt ) − τ p ∫ t u ( x p ( τ ) , τ ) e τ τp ( cos ( ω ( t − τ )) + Sτ p − τ p ω sin ( ω ( t − τ )))] + τ p u ( x p ( t ) , t ) . (19)In the expressions, λ , = − τ p ( ± √ − Sτ p ) , ω = τ p √ Sτ p − τ p determine time scales for various contributions to inertialparticle accelerations. x p and v p are position and velocity of the particle right before thestraining starts (at t = − ). In the expanding direction y : (the expression in the z direction issimilar) ̃ a p ( t ) = λ − λ [( λ ( λ − S ) e λ t − λ ( λ − S ) e λ t ) y p + ( λ e λ t − λ e λ t ) v p + τ p ∫ t u ( x p ( τ ) , τ )( λ e λ ( t − τ ) − λ e λ ( t − τ ) ) dτ ] + τ p u ( x p ( t ) , t ) . (21)Similarly, y p and v p are position and velocity of the particle at t = − , and λ , = − τ p ( ± √ + Sτ p ) . (22)The mean flow influences the variances when full acceleration is considered, and so is the casefor other statistics involving full particle velocity or accelerations. Since the magnitude of themean flow velocity depends on the location in the domain, the magnitudes in the variances of14 × t h ˜ a p i , h ˜ a p i S × t S = 0 . , xS = 0 . , yS = 1 , xS = 1 , yS = 4 , xS = 4 , yS = 10 , xS = 10 , y Figure 7: Acceleration variances of inertial particles in the compressed (solid symbols) andexpanding (empty symbols) directions vs. strain time for two types of particles Left: τ p = . τ p = .
05. Symbols refer to different strain rates, diamond ( ◇ ) S = .
1; circle ( ◯ ) S =
1; square ( ◻ ) S =
4; triangle ( △ ) S =
10. The data are from the data set R λ = X j being particleacceleration variances in the j -th realization.acceleration components are characterized by the domain size in respective directions (throughthe initial positions x p , y p in acceleration expressions (18), (19) and (21)) in addition to therate of strain.In an attempt to minimize the influence of the mean flow, in addition to ensuring thatour statistics are deduced from sufficient many independent samples, we condition our inertialparticle analysis on particles that started in a thin layer parallel to and next to the x = x -component statistics and a thin layer parallel to and next to the y = y -component statistics. For the strain rates S = . S =
1, we use layers with thickness of8 lattice units in the R λ =
117 simulations. For S = S =
10 flows, we use layers withthicknesses of 4 and 2 lattice units. With the selected thicknesses, we ensure the differenceof the mean flow velocities within the layers do not exceed 55% of u rms in the compressiondirection, and do not exceed 77% of u rms in the expanding direction. The numbers of particlesavailable in the layers decreases with the thickness of the layers from about 61,500 in S = S =
10 flows in the compression direction. In the expanding direction, thenumber of particles used for the analysis decreases from 240,000 in S = S = x = y = τ p = .
015 and τ p = .
05, which correspond to St = . S =
1, correspond to D > D < S =
4, and have negative discriminants in S = × t h e a p i , h e a p i tracer, S = 4 ,x tracer, S = 4 ,y τ p = 0 . ,S = 4 ,x τ p = 0 . ,S = 4 ,y τ p = 0 . ,S = 4 ,x τ p = 0 . ,S = 4 ,y S × t tracer, S = 10 ,x tracer, S = 10 ,y τ p = 0 . ,S = 10 ,x τ p = 0 . ,S = 10 ,y τ p = 0 . ,S = 10 ,x τ p = 0 . ,S = 10 ,y Figure 8: Acceleration variances for tracers and inertial particles with τ p = . , .
1. Lines withcircles ( ◯ ): τ p = .
1; lines with triangles ( △ ): τ p = .
01. Solid lines with solid symbols: ⟨(̃ a p ) ⟩ ;dashed lines with empty symbols: ⟨(̃ a p ) ⟩ . Solid lines (—): tracers ⟨( a p ) ⟩ ; dashed lines ( −− ):tracers ⟨( a p ) ⟩ . Left: S =
4, Right: S =
10. The data are from the data set R λ = X j being particleacceleration variances in the j -th realization. The acceleration variances of the inertial particles show a transition period at the beginning ofthe straining that is not seen in the tracer accelerations. Equations (18) and (19) indicate thatthe initial position x p and velocity v p of an inertial particle affects its acceleration. The timescale of this influence depends on the exponents of the exponential terms in the accelerationexpressions. In the compression direction, if D < τ p . If D >
0, both λ and λ are negative so the decaying rate is min ( ∣ λ ∣ , ∣ λ ∣ ) = ∣ λ ∣ . This explains the differentlengths of the initial transition period for particles with different Stokes numbers in flows with afixed strain rate. In the left figure of Figure 8, we show acceleration variances of particles with τ p = .
01 and 0 . S = ⟨( u ) ⟩ decreases with straining, ⟨̃ a p ⟩ increases at the beginning and a maximumoccurs. This is because λ is positive in the expanding direction. Similar to the case discussed inthe x -direction, the exponents determine the time scale of the behavior of acceleration variances.From (21) one can estimate the time that the maximum takes place by ignoring the flowfluctuation, and obtain the ratio between the times of maximum acceleration variances of thetwo particles presented in the figure to be 1.52. It appears that our simulation data confirmssuch an estimation. For S =
10, our simulations were too short to display the post-transitionperiod.
Figure 7 shows that the acceleration variances of inertial particles increases with the strain rate.This is mainly due to the mean flow contributing to the acceleration in terms of the factors 2 S and S in the expressions (18), (19) and (21) and, in a smaller part, through the exponents thatregulate the influence of the initial conditions.16t the onset of straining, the acceleration variances in the expanding direction is higherthan that in the compression direction. This could be explained through the exponents of thekernels. In the compression direction, the coefficients λ , λ and − τ p in the exponents of thekernels in (18) and (19) are all negative and indicate decay with time. But for the expandingdirection, λ = − + √ + Sτ p τ p is positive, which leads to the increase of magnitude. For sufficientlylong time, the influence of the initial condition is diminished. The magnitude of the flow velocityfluctuation (in the last two terms of the acceleration expressions) in the expanding directionremains roughly constant, while in the compression direction the flow velocity fluctuation growsexponentially. As a result the acceleration variances in the compression direction overtakes thatin the expanding direction.In contrast to the isotropic homogeneous turbulence (Ayyalasomayajula et al. (2006); Bec et al. (2006)), the acceleration variances of inertial particles do not necessary have lower magnitudescompared with that of tracers in strained flow. The inertial particle with τ p = .
01 in theright plot of Figure 8 depicts such a situation. Recall that (17) provides an approximation oftracer acceleration variances, and in higher strain rates the terms 2 S i ⟨( u i ) ⟩ + S i d ⟨( u i ) ⟩/ dt are the main contributor of the variances. From Figure 3 we know that the variance of flowfluctuating velocity grows exponentially in time, so its time derivative can be estimated as aconstant multiple c S ( c is a constant) of the variance itself. Hence the main terms in thetracer acceleration variance increases in the order of S i ⟨( u i ) ⟩ with the strain rate. For theinertial particles, however, their long term variances depends on the term ⟨( u i ) ⟩/ τ p . When τ p is small enough, the factor 1 / τ p is larger than the magnitude of S i , and the inertial particleacceleration variances surpass the tracer acceleration variances. Figures 9 to 11 show the particle acceleration pdfs at various rates of strain ( S = , R λ =
117 and R λ = τ p = . τ p = . ⟨( a p i − a p i ) ⟩/⟨( a p i − a p i ) ⟩ (a measure of theintermittency of the acceleration) is shown in the bottom right figure for the lower Reynoldsnumber R λ = et al. (2006); Bec et al. (2006)), owing to the heavier particles passing through, or selectively filtering,the most rapid motions in the flow field. However, here the situation is more complex, sincethe acceleration variance is not necessarily smaller for inertial particles, due to the strong effectof the mean flow as discussed above. Figure 11 shows that the higher Reynolds number has17 p /a rmsp -20 -10 0 10 20 P d f ( a p ) -6 -5 -4 -3 -2 -1 HITS = 4 S = 10 a p /a rmsp , a p /a rmsp -20 -10 0 10 20 P d f ( a p ) , P d f ( a p ) -6 -5 -4 -3 -2 -1 HIT, xHIT, yS = 10 , xS = 10 , yS = 4 , xS = 4 , y S × t h ( a p i − a p i ) i / h ( a p i − a p i ) i HITS = 4 , xS = 4 , yS = 10 , xS = 10 , y a p /a rmsp -20 -10 0 10 20 P d f ( a p ) -6 -5 -4 -3 -2 -1 HIT, R λ = 117 HIT, R λ = 193 S = 10 , R λ = 117 S = 10 , R λ = 193 S Re Figure 9: Normalized pdfs of tracer accelerations and tracer acceleration flatness. Top left:Pdf( a p ) in HIT (circle • ), straining flow S = ∎ ), and S =
10 (triangle ▲ ). Topright: Pdf( a p ) in HIT (circle • : R λ = ○ ∶ R λ = S = ▲ ∶ R λ = △ ∶ R λ = a p ) and Pdf( a p ) in R λ = • , ○ ): HIT, squares ( ∎ , ◻ ): S =
4, triangles ( ▲ , △ ): S =
10. Solid sym-bols: i = i = ⟨( a p i − a p i ) ⟩/⟨( a p i − a p i ) ⟩ at R λ = ∎ , ◻ ): S = ▲ , △ ): S =
10. Solid symbols: i = i = a p / e a rmsp -15 -10 -5 0 5 10 15 P d f ( e a p ) -6 -5 -4 -3 -2 -1 St = 0 . R λ = 117) HITS = 4 S = 10 e a p / e a rmsp -10 -5 0 5 1010 -6 -5 -4 -3 -2 -1 St = 1( R λ = 117) e a p / e a rmsp -15 -10 -5 0 5 10 15 P d f ( e a p ) -6 -5 -4 -3 -2 -1 e a p / e a rmsp -10 -5 0 5 1010 -6 -5 -4 -3 -2 -1 SS SS
Figure 10: Normalized probability distribution function of particle acceleration components ̃ a p i for different strain rates, for particles originating in thin slices parallel to the x = i = y = i = St = . R λ = St = R λ = i = i = • ): HIT; squares ( ∎ ): S =
4; triangles ( ▲ ): S = St =
1) in higher strain S =
10, the effect is not as evident in the expanding direction.We note that the Stokes number of a particle increases slightly during the straining due toa decrease in the Kolmogorov timescale (less than 10% for the highest rate of strain). Thiseffect contributes to the narrowing of the acceleration pdf tails, but to a lesser effect than theincreased variances.It is interesting to consider our results in a wider context of flows with non-zero meancomponents, for example by comparing with the dynamics in shear flow. Lavezzo et al. (2010)and Gerashchenco et al. (2008) considered tracers and inertial particles in a non-uniform shear,namely a turbulent boundary layer. An increased rate of shear, closer to the wall, resulted inan increased acceleration variance in the lateral component but a milder effect in the transversecomponent, which appears to be consistent with the predictions of equation (12) for theirgeometry. As a result of the increased acceleration variances, attenuation was found in the tailsof the inertial particle acceleration pdfs as observed in this work.19 a p / e a rmsp -15 -10 -5 0 5 10 15 P d f ( e a p ) -6 -5 -4 -3 -2 -1 St = 0 . R λ = 117) , St = 0 . R λ = 193) e a p / e a rmsp -10 -5 0 5 1010 -6 -5 -4 -3 -2 -1 St = 1( R λ = 117) , St = 1 . R λ = 193) e a p / e a rmsp -15 -10 -5 0 5 10 15 P d f ( e a p ) -6 -5 -4 -3 -2 -1 e a p / e a rmsp -10 -5 0 5 1010 -6 -5 -4 -3 -2 -1 Re ReRe
S SS S
Figure 11: Normalized probability distribution function of particle acceleration components ̃ a p i for two Reynolds numbers, for particles originating in thin slices parallel to the x = i = y = i = St = . R λ =
117 and St = .
34 in R λ = St = R λ = St = .
12 in R λ = i = i = • , ○ ): HIT; triangles ( ▲ , △ ): S =
10. Solid symbols: R λ = R λ = Conclusions
Our results show a strong effect of the mean flow straining on the Lagrangian accelerationstatistics, both for the passive tracers and the inertial particles. The effect of straining isprimarily felt in acceleration variances and pdfs when the rate of strain is sufficiently high suchthat the strain timescale is of a comparable magnitude to the Kolmogorov timescale. For highrates of strain the magnitude of the acceleration variances are increased significantly and thetails of the normalized acceleration pdfs for tracers and inertial particles narrow. The formereffect is well explained by observing the predicted behavior of the acceleration variance byrapid distortion theory. RDT provides us a relation between the flow velocity variances and itsacceleration variances and illustrates the dependence of the acceleration variance on the rate ofstrain.However, the effect is complex, partly due to the connection, or lack thereof, between particleacceleration component in one direction and the fluid flow in the same direction; a particletrajectory around a strong vortex will result in large acceleration values in the directions normalto the axis of the vortex. But there is also a direct contribution from the mean straining andfluctuating velocity in the acceleration, resulting for example in an increased variance value ofacceleration in both compressing and expanding directions.For tracers, the narrowing of the normalized acceleration pdfs stems therefore in a complexmanner from both of these effects. The same effect is also felt by the lighter inertial particles,given that their inertia is sufficiently small to sample the small scale motions. Because of theirsmall inertia the interplay with the mean flow enable their acceleration variances to rise beyondthat of tracers.When the inertia is further increased the ballistic particle motion in the rapidly acceleratingmean flow becomes increasingly important leading to Gaussian pdf tails. Here, the particlesare swept through the fluid, and the slower, large scales of motions are more likely to influencethe particles. This could also be seen by the lower magnitude of acceleration variances of theheavier particles compared with the lighter inertial particles. We derive the formal expressionsfor inertial particle acceleration, and these expressions reveal the complex interplay betweenflow straining and particle inertia.Our findings emphasize the importance of the presence of strong mean motions and imposedsmall scale anisotropy in particle laden flows. It is our opinion that the results have relevanceto the understanding and modelling of a range of practical deforming or straining flows whereinertial particles are important aspects of the process. In particular, we believe that our findingsmay help in the development of sub-grid Lagrangian models for particles in the proximity ofstraining regime near stagnation points.We thank Zellman Warhaft, Paolo Gualtieri, and Michel van Hinsberg for their commentsand suggestions. C. Lee acknowledges the hospitality of the School of Science and Engineer-ing at Reykjav´ık University and the Department of Mathematics and Computer Science andthe Department of Applied Physics at Eindhoven University of Technology. This work wassupported by the Icelandic Research Fund and was partially supported by a research programof the Foundation for Fundamental Research on Matter (FOM), which is part of the Nether-lands Organisation for Scientific Research (NWO). Support from COST Actions MP0806 andMP1305 is also kindly acknowledged.
References
Ayyalasomayajula, S., Gylfason, A., Salazar, J., Warhaft, Z., Collins, L. R. Bodenschatz, E.
ASME Fluids Conference . Ayyalasomayajula, S. & Warhaft, Z.
J. Fluid Mech. , 273–307.
Batchelor, G.K
The Theory of Homogeneous Turbulence . Cambridge University Press.
Bec, J., Biferale, L., Boffetta, G., Celani, A., Cencini, M., Lanotte, A., Musac-chio, S. & Toschi, F.
J. Fluid.Mech. , 349.
Biferale, L. & Procaccia, I.
Physics Reports (2-3), 43 – 164.
Celani, A.
J.Turbul. , N34. Chen, J., Meneveau, C. & Katz, J.
J. Fluid Mech. , 123–150.
Gerashchenco, S., Sharp, N. S., Neuscamman, S. & Warhaft, Z.
J. Fluid Mech. , 255–281.
Gualtieri, P. & Meneveau, C.
Phys. Fluids , 065104. Gualtieri, P., Picano, F. & Caisciola, C.M.
J. Fluid Mech. , 25–39.
Gualtieri, P., Picano, F., Sardina, G. & Casciola, C.M.
Physica D , 245–250.
Gylfason, A., Ayyalasomayajula, S. & Warhaft, Z.
J. Fluid Mech. , 213–229.
Gylfason, A., Lee, C., Perlekar, P. & Toschi, F.
Journal of Physics: Conference Series ,052003.
Han, Z. & Reitz, R.D. κ − ǫ models. Combustion Science and Technology , 267–295.
Hunt, J.C.R.
J. Fluid.Mech. , 625–706. Hunt, J.C.R. & Carruthers, D.J.
J. Fluid Mech. , 497–532.
Klein, A.
Prog. Aerospace Sci. , 171–271.22 avezzo, V., Soldati, A., Gerashchenko, S., Warhaft, Z. & Collins, L.R. J. FluidMech. , 229–246.
Maxey, M.R. & Riley, J.R.
Phys. Fluids (4), 883–889. Ott, S. & Mann, J.
J. Fluid. Mech. , 207–223.
Pope, S. B.
Turbulent Flows . Cambridge University Press.
Rogallo, R. S.
Tech. Rep.
Shaw, R.
Ann. Rev. Fluid. Mech. , 183–227. Toschi, F. & Bodenschatz, E.
Ann.Rev. Fluid Mech. , 375–404. Virant, M. & Dracos, T.
Meas. Sci.Tech. , 1539–1552. Voth, G.A., Porta, A. La, Crawford, A.M., Alexander, J. & Bodenschatz, E.
J. Fluid Mech ,121–160.
Warhaft, Z.
J. Fluid Mech. , 545–573. Xu, H., Ouelette, N.T. & Bodenschatz, E.
New J. Phys. , 013012. Yeung, P.K. & Pope, S.B.
J. Comput. Phys.79