A method to calculate finite-time Lyapunov exponents for inertial particles in incompressible flows

The present study aims to improve the calculus of finite-time Lyapunov exponents (FTLEs) applied to describe the transport of inertial particles in a fluid flow. To this aim, the deformation tensor is modified to take into account that the stretching rate between particles separated by a certain distance is influenced by the initial velocity of the particles. Thus, the inertial FTLEs (iFTLEs) are defined in terms of the maximum stretching between infinitesimally close trajectories that have different initial velocities. The advantages of this improvement, if compared to the standard method (Shadden et al., 2005), are discussed for the double-gyre flow and the meandering jet flow. The new method allows one to identify the initial velocity that inertial particles must have in order to maximize their dispersion.


Introduction
The advection of finite-size or inertial particles in open, unsteady flows was initially assessed independently by Maxey and Riley (1983) and Gatignol (1983).In contrast to the idealized Lagrangian view, the two authors consider a number of additional forces acting upon inertial particles moving in a flow such as buoyancy, the Stokes drag, the added mass, or the Basset-Boussinesq memory force (see Michaelides (2003) for a review).In this view, the motion of a spherical small solid particle in an incompressible two-dimensional flow can be approximated by a Stokes flow (with Basset-Boussinesq memory force and Faxén terms neglected) governed by the equation (Babiano et al., 2000;Bec, 2003;Boffetta et al., 2004) where β = 3ρ f /(2ρ p + ρ f ), τ p = a 2 /(3νβ) is the Stokes time or particle viscous response time, v i = dx i /dt the velocity of the sphere and u i that of the fluid, ρ p (ρ f ) the density of the particle (fluid) it displaces, ν = µ/ρ f the kinematic viscosity, and a the radius of the sphere.The derivative D/Dt is taken along the path of a fluid element, whereas the derivative d/dt is taken along the trajectory of the particle.The dimensionless Stokes number St = τ p /τ f characterizes the effect of particle inertia and τ f accounts for the characteristic timescale of the fluid.A modified version of Eq. ( 1) has been recently studied in the context of compressible fluids (Pérez-Muñuzuri, 2015).
During the last years, the dynamics of inertial particles have been studied in many research fields such as sedimentation processes (Santamaria et al., 2013), turbulent flows (Boffetta et al., 2004), rain generation (Falkovich et al., 2002), composite materials (Segurado et al., 2003), volcanic ash transport (Haszpra and Tél, 2011), and the formation of planetesimals in the early Solar system (Bracco et al., 1999), and other processes as transport of dust from soil erosion, combustion or the mixing of sprays.
Although the transport of Lagrangian particles as described by finite-time Lyapunov exponents (FTLEs) has been extensively studied to date, to our knowledge only a few studies have been performed with inertial particles (Haller and Sapsis, 2008;Sapsis and Haller, 2009;Peng and Dabiri, 2009;Beron-Vera et al., 2015).Traditionally, FTLEs are calculated as the logarithm of the maximum eigenvalue of the right Cauchy-Green strain tensor.They measure the maximum stretching rate at a given location r 0 over the time Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.

572
D. Garaboa-Paz and V. Pérez-Muñuzuri: FTLE for inertial particles interval τ = t − t 0 of trajectories starting near the point r 0 at time t 0 (Shadden et al., 2005).However, the motion of inertial particles following the simplified model mentioned above requires not only the knowledge of their initial positions, but also that of their initial velocity fields.In this study, a modified procedure for calculating inertial FTLEs (iFTLEs) is proposed that takes into account this additional factor.Considering two idealized flows (the double-gyre and the meandering jet), the results obtained by the new method are illustrated and compared to those of the traditional one.

Methods
For the dynamical system described above (Eq.1), an initial fluid element with coordinates λ = (r 0 (t 0 ), v 0 (t 0 )) is advected by the flow u(r, t) during a finite-time period τ , = φ t 0 +τ t 0 λ until reaching the final position = (r(t 0 + τ ), v(t 0 + τ )), such that the Cauchy problem, holds.The gradient of the flow map (the deformation gradient tensor F) is F m,n = ∂ m /∂λ n , and the right Cauchy-Green strain tensor is then defined as At each initial condition λ, the tensor C(λ) is represented by a symmetric, positive definite, 4 × 4 matrix (for a two-dimensional flow), with four positive eigenvalues µ 1 > . . .> µ 4 > 0 and four eigenvectors (ξ k (λ)) 1≤k≤4 satisfying In order to characterize the transport of inertial particles, we introduce the iFTLEs σ that are computed along their trajectories in the flow as (Peacock and Dabiri, 2010) where √ µ 1 denotes the ratio of stretching between two particles that are initially located close to each other and have different initial velocities in the direction of the largest stretching.The iFTLE at a given location measures the maximum stretching rate of an infinitesimal fluid parcel over the interval [t 0 , t 0 + τ ] starting at the point λ at time t 0 .Repelling (attracting) coherent structures for τ > 0 (τ < 0) can be thought of as finite-time generalizations of the stable (unstable) manifolds of the system.These structures govern the stretching and folding mechanism that control the mixing of the flow of inertial particles.
For the computation of the iFTLE field, the integration time τ and the initial conditions λ at time t 0 must be predefined.Basically, τ has to be long enough to allow the trajectories to explore the coherent structures present in the flow.
For the initial state λ at time t 0 , we are interested in studying the iFTLEs as a function of the direction and intensity of the particle velocity with respect to the velocity of the flow u.To this end, instead of using the Cartesian components of the particle velocity, we will use from now on its polar coordinates defined by v 0 = R|u|(cos(α + α 0 ), sin (α + α 0 )), with R > 0 and α 0 the direction of the flow velocity u at t 0 .Then, the elements of the deformation gradient tensor are computed by a finite-difference method as and similarly for the rest of the terms.The indices i, j correspond to the directions x, y, respectively, and the superscripts 1 and 2 to two layers of particles with different initial velocities, M k (M k 0 ) and α k (α k 0 ) account for the modulus and direction of the particle velocities v k at t = τ (t = t 0 ), respectively, in both layers k = 1, 2. The advantage of defining the initial velocities as above allows one to use a constant initial gradient δα 0 for the deformation tensor, although some caution must be taken in order to select R in δM 0 (see below), with the distance between particles initially located in a regular grid.
The results obtained from using the improved iFTLE calculation method will be compared to those obtained from the standard method (Shadden et al., 2005), which only takes into account deformation due to the initial positioning (dim(F) = 2).It is recognized that initial velocity has been considered for this purpose in previous studies (see for example Peng and Dabiri, 2009).These studies, however, did not explicitly consider its effect on the deformation tensor F, as is done above.

Results
We investigate the effect of initial velocity v 0 on inertial particle motion in two chaotic flows: the double-gyre flow and the meandering jet flow (see Appendix A for equations).The dynamics of these flows are completely different.While the first model describes a closed flow, the second one describes an open one.In the double-gyre model, the flow is the combination of two vortices rotating in opposite directions that conversely expand and contract periodically in the x direction such that the rectangle enclosing the gyres remains in- variant.The periodic perturbation leads to mixing between both gyres.On the other hand, the jet flow is a simple kinematic model for the Gulf Stream and is frequently used to describe western boundary current extensions in the ocean.Three different regimes can be observed in this flow and the presence of temporal variability in the amplitude of the meander favors mixing and chaotic transport.The (x, y) projection of the iFTLE fields calculated with the method depicted above is shown in the upper rows of Figs.1-2, assuming two different initial angles α.For comparison, the FTLE fields calculated with the standard method (Shadden et al., 2005) are shown in the lower rows of these figures.The two methods are applied to identical initial velocities.Regarding the double-gyre flow (see Fig. 1), the pattern is sensitive to variations in α for both methods in the sense that the ridges are not equally distributed.Due to inertia, particle motion is affected by the distribution of the initial velocities, as was already noted by Peng and Dabiri (2009).The iFTLE patterns, however, provide new information as they indicate the regions where larger mixing is expected for a given angle α.The projection of the four-dimensional iF-TLE onto the (x, y) plane results in a stroboscopic pattern characterized by crossing ridges.
A comparison between Figs. 2 and 1 reveals that the method-related differences in the FTLE fields are larger for the meandering jet flow than for the double-gyre flow.For the meandering jet flow (see Fig. 2), the maximum dispersion (as displayed by the red areas) is attained within the central meandering jet in all cases.For α = 90 • , however, the fila- mentary structures obtained by our method are wider than those obtained by the standard method (compare upper and lower panels in Fig. 2).This indicates that the largest dispersion is achieved by those particles that (i) are initially located near the axis of the jet and (ii) have an initial velocity perpendicular to the flow (v 0 • u = 0).In this case, particles are trapped within the circulation regimes located above the meander troughs and below the meander crests where the mixing of the flow is larger than in the remaining areas.This result is consistent with physical intuition and demonstrates that the initial velocity of inertial particles contributes to their dispersion in the flow.For ε = 0, due to the periodic modulation of the meander amplitude, mixing is favored between the circulations and the peripheral currents and between the circulation cells and the jet (Prants et al., 2006).Particles are trapped in the eastward currents and in the jet where they may spend longer periods of time.For this flow, dispersion is largest in the circulation cells located above the meander troughs and below the meander crests.The new method allows one to identify those areas where maximum dispersion is attained.
Figure 3 shows the mean iFTLE as a function of α (on the x axis) and R (as represented by distinct lines) for heavy and light inertial particles (right and left columns, respectively) in the two considered flows.A modification of the initial velocity R leads to distinct local iFTLE maxima along the range of possible α values.The mean iFTLEs obtained by our method are systematically more sensitive to changes R and/or α than the corresponding values obtained by the standard method (which is insensitive to changes in R by definition).These differences are larger for the heavier particles than for the lighter ones (compare right and left columns) and are particularly pronounced for the case of a chaotic jet flow (ε = 0.5).
The mean FTLEs are larger for the heavier particles (β = 0.2) than for lighter ones (β = 3.0) because clustering diminishes with increasing ρ p .Heavier particles are less influenced by the flow (they move slowly) and they are less likely to be displaced towards the trapping areas.Without periodic modulation (ε = 0), light particles concentrate inside the inner cores of the vortices as they are driven inward due to the lesser inertial.On the contrary, heavy particles tend to move away from the inner core to regions characterized by low vorticity and high strain, a process known as preferential concentration (Eaton and Fessler, 1994).Thus, when u × v 0 and ∇ × u have the same sign, the dispersion within the vortices will decrease for lighter particles and increase for heavier ones.Next, the iFTLE-field response to changes in the initial conditions is assessed by means of the eigenvector, decomposed into the factors "initial position" and "initial velocity".To this end, we calculate the ratio between both, where < ξ 1 > (M,α) is the spatial average of the eigenvector ξ 1 corresponding to the largest eigenvalue, Eq. ( 4), projected onto the (M, α) plane.Figure 4 shows as a function of R and α.For some angle ranges, and depending on the initial modulus R, the initial velocity ( > 1) is the largest contributor to the main eigenvector ξ 1 ( ).Note, for example, that for light particles in the double-gyre flow there is nearly no dependence on the initial velocity parameters.Thus, is a parameter describing the sensitivity of the iFTLE to changes in the initial position and velocity.

Conclusions
In the present study, a new method for calculating iFTLEs used to describe the motion of inertial particles in a fluid flow was suggested, which, in contrast to the standard method (Shadden et al., 2005;Peng and Dabiri, 2009), takes into account the particle initial velocity.It is shown that the results obtained with the new method are sensitive to variations in initial velocity, indicating that the suggested modifications are meaningful from a practical point of view.The presence of different flow regimes leads to larger differences in the particle dispersion when taking into account the initial velocity, as particles may change between different regimes.Thus, for the jet flow, initial velocities perpendicular to the flow close to the jet stream regime lead to larger dispersion.More importantly, our method does not only hold for particle motions described by the Maxey-Riley equation (Eq.1), but can be applied to any type of equation used to describe the motion in an N-dimensional flow as well as the associated transport and mixing properties.It can be generalized to any type of particle motion, such as the dispersion of marine litter (micro and macro plastics) (Lebreton et al., 2012), or the motion of drifting buoys near the coast where inertial effects may be particularly relevant (Huhn et al., 2012).

Figure 1 .
Figure 1.Finite-time Lyapunov exponents σ (x, y) for heavy inertial particles in a gyre flow, assuming two different initial angles α.The iFTLE fields in the upper row were calculated using Eq.(5), while those in the lower row were calculated using the standard method.The model parameters are R = 0.6, A = 0.1, ε = 0.25, T = 10, τ p = 1, and β = 0.2, and the integration time is τ = 30.

Figure 3 .
Figure 3. Spatial average FTLE as a function of the initial velocity parameters R (displayed by distinct lines) and α (displayed along the x axis) for light (β = 3.0, left column panels) and heavy (β = 0.2, right column panels) inertial particles and the double-gyre and meandering jet flows.Color lines correspond to the results obtained by the new method Eq. (5), whereas black dashed lines correspond to those obtained by the standard method.Rest of the parameters as in Figs.1-2.

Figure 4 .
Figure 4. Contribution to the eigenvector ξ 1 corresponding to the largest eigenvalue due to the initial velocity and position , Eq. (9), as a function of the parameters R and α for light (β = 3.0, left column panels) and heavy (β = 0.2, right column panels) inertial particles for the double-gyre and meandering jet flows.For > 1, the main contribution to the eigenvector is due to the initial velocity.Rest of the parameters as in Figs.1-2.