Tracking Particles in Flows near Invariant Manifolds via Balance Functions
aa r X i v : . [ m a t h . D S ] A ug Tracking Particles in Flows near Invariant Manifolds via BalanceFunctions
Christian Kuehn †¶ Francesco Roman`o ‡ Hendrik C. Kuhlmann ‡ September 1, 2016
Abstract
Particles moving inside a fluid near, and interacting with, invariant manifolds is a com-mon phenomenon in a wide variety of applications. One elementary question is whether wecan determine once a particle has entered a neighbourhood of an invariant manifold, when itleaves again. Here we approach this problem mathematically by introducing balance functions,which relate the entry and exit points of a particle by an integral variational formula. Wedefine, study, and compare different natural choices for balance functions and conclude thatan efficient compromise is to employ normal infinitesimal Lyapunov exponents. We apply ourresults to two different model flows: a regularized solid-body rotational flow and the asymmetricKuhlmann–Muldoon model developed in the context of liquid bridges. Furthermore, we employfull numerical simulations of the Navier-Stokes equations of a two-way coupled particle in ashear–stress-driven cavity to test balance functions for a particle moving near an invariant wall.In conclusion, our theoretically-developed framework seems to be applicable to models as wellas data to understand particle motion near invariant manifolds.
Keywords:
Invariant manifold, fluid dynamics, entry-exit function, perturbation theory, particle-surface interaction, fast-slow systems, nonautonomous dynamics.
Tracking particles inside a flow is a topic of general importance in a wide variety of applications,ranging from small scales in microfluidics [69], to mesoscale problems in sedimentation [71], to largeoceanic scales [47], and even astrophysical scales [7]. Many results for particle motion are imme-diately useful for very classical problems such as turbulent pipe flow [74] as well as newly arisingrecent challenges, for example in the search for parts of crashed airplanes [62]. The dynamics ofparticles within certain domain and flow geometries has been studied intensively, e.g., for parti-cles moving near walls [8], colliding with walls [28], with respect to rotation-induced forces [64] orchanging viscosity [18], or via quite general partial differential and integral equations of motion forparticles [46, 55]. Another important branch of current research is to study separation structures † Technical Unversity Munich, Fakult¨at f¨ur Mathematik, Boltzmannstr. 3, 85748 Garching bei M¨unchen, Germany ‡ Institute of Fluid Mechanics and Heat Transfer, Vienna University of Technology, Getreidemarkt 9, 1060 Vienna,Vienna, Austria ¶ CK acknowledges support via an APART fellowship of the Austrian Academy of Sciences ( ¨OAW) and via aLichtenberg Professorship of the VolkswagenFoundation. Furthermore, CK would like to thank Peter Szmolyanand Daniel Karrasch for general discussions regarding non-autonomous dynamical systems, and Armin Rainer for aclarifying remark regarding differentiability of certain parametrized eigenvalues.
1n flows, which act as transport barriers, i.e., particles cannot pass through them. The defini-tion [23, 19, 68], analysis [20, 14], and computation [16, 42] of these Lagrangian coherent structureshas received a lot of attention; see also the references in the recent review [21] and the focus is-sue [54]. A directly related topic are invariant manifolds for general dynamical systems [12, 25].The autonomous (time-independent vector field) has to be studied first but since transport barriersin flows are often viewed in terms of non-autonomous dynamical systems [31, 58], i.e., as time-dependent invariant manifolds [1, 6, 10], we provide a setup, which also works for non-autonomouscases. t = t t ∈ I t = T M MM EEE
Figure 1: Sketch of the basic situation studied in this paper. A particle (grey disc) moves near aninvariant manifold M inside a flow. The question is how to relate the entry and exit points into asmall tubular neighbourhood E of M ; note that M = M ( t ) and E = E ( t ) could be time-dependentfor general fluids. The key calculation is to relate the time t at entry, via a drift time near M inside some compact time domain I , to the final exit time T .Despite considerable progress, there seems to be no completely general mathematical frameworkavailable for some relatively elementary-looking problems involving particles. One such question is,how to determine the entry-exit relationship for a particle near an invariant structure; see Figure 1.In Section 2, we provide further background from experiments, which motivated our work on entry-exit problems.In this paper, we develop a mathematical framework based on balancing ideas, i.e., we ask thequestion which quantity is adequate to compute (or measure) the entry-exit relationship of a parti-cle moving near a time-dependent invariant manifold. The main idea of this work is motivated bya classical problem in multiple time scale systems [35]. In this context, the invariant manifolds ofinterest are slow manifolds [27, 13]. Fast transverse dynamics to slow manifolds may be attractingas well as repelling in certain directions. Of particular interest are trajectories which first get at-tracted to slow manifolds, pass near certain bifurcation points of the fast dynamics, and eventuallyget repelled from a neighbourhood of the slow manifold. This problem has been studied intensivelyduring recent years in the context of various bifurcations such as fold points [11, 33], Hopf bifurca-tion [51, 2], as well as more complicated singularities [72, 34]. One key idea in this context is theso-called entry-exit (or way-in/way-out, or input-output) map which can be used to calculate therelationship between entry and exit from a neighbourhood of a slow manifold involving a stabilitychange. Here we substantially extend this idea to the case of fluid flows, as well as to arbitrary di-mensions, and to various different possible maps involving entry-exit from time-dependent invariantmanifolds. Then we apply it to several examples motivated by recent experiments and models ofparticle-surface interaction [49, 70, 39] as well as to full two-dimensional Navier-Stokes simulationsfor a particle motion fully coupled to the fluid. Our main results (in non-technical form) are thefollowing: Section 3:
We introduce several natural notions for balance functions to determine entry-exitrelationships based upon instantaneous eigenvalues, finite-time Lyapunov exponents (FTLE) [42,19], fast-slow scaling [35] and normal infinitesimal Lyapunov exponents (NILE) [22]. We prove2uitable differentiability as well as continuous dependence on data for the balance functions anddevelop perturbation theory results.
Section 4:
To evaluate which concept is most promising, we investigate two test models (solid-body rotation and the Kuhlmann–Muldoon model), where the particle does not actively influencethe flow. Three concepts turn out to have some drawbacks while NILE seems sufficiently easy tocompute and conceptually the most adequate for our purpose.
Section 5:
We consider direct numerical simulation (DNS) of the Navier-Stokes equations for atwo-way coupled particle in a shear–stress-driven cavity. In this context, there is an invariant free-surface in the cavity, which is tracked by the particle for a transient period. We use the DNS datato compute the NILE balance function. The results show that, although the NILE balance functionis only an integrated locally linearized estimator near the boundary, it performs very similar to afull nonlinear analysis of the particle velocities.
Section 6:
We conclude with an outlook of potential applications and future challenges. Webelieve many interesting directions could be pursued beyond the fundamental groundwork presentedin this study. In particular, there are interesting questions from theoretical, modelling, and dataanalysis perspectives.
An example for an entry-exit problem of current interest [56, 38, 48, 40] are particle accumulationstructures (PAS) [70, 66]. The phenomenon during which particles suspended in a liquid are rapidlysegregating.
Remark:
The final version of this paper contains a figure from [66] illustrating this experiment,which we cannot reproduce here on the arXiv.A natural transport barrier for fully wetted particles is the liquid–gas interface of a cylindricaldroplet (liquid bridge) in which a flow is driven by tangential shear stresses created by temperature-induced axial variations of the surface tension [37]. According to [26] and [49] PAS are attractorsfor the motion of small particles in a three-dimensional steady flow which are created in a closeneighbourhood of the invariant manifold represented by the liquid–gas interface. The mechanisminvolves the particle–free-surface interaction during which a particle enters a certain neighbourhoodof the invariant manifold, experiences extra forces due to the presences of the interface, and exitsto the bulk. If the properties of the entry–exit process are such that particles enter from a regionof chaotic streamlines and exit to a region of regular streamlines of the steady three-dimensionalflow a periodic attractor is typically created, in particular, if the particle inertia is small.A direct numerical simulation of this process based on the Navier–Stokes equations would requireenormous computing resources, because all length scales must be properly resolved, ranging fromthe scale of the droplet over the scale of the particles down to the scale of the lubrication gap betweena particle’s surface and the liquid–gas interface. For that reason suitable entry-exit relations whichgo beyond inelastic collision models [26] would be extremely helpful for a realistic and economicsimulation of PAS using one-way coupling, e.g. in the framework of the Maxey–Riley equation [46];we refer to Section 4.2 for the analysis of a first model using balance functions in this direction.Analytical approaches would also be a first step to check whether there is a stretching mechanismnear the surface so some particles stay there a lot longer than others although they entered innearby regions. Such a stretching mechanism is key for the generation of chaotic dynamics inmany situations including problems involving the classical Smale horseshoe [17]. Furthermore,3nalytical entry-exit relationships could be used for perturbation theory (see Section 3.6) as wellas for extrapolation of particle locations based upon previous particle data (see Section 5).Further real-world examples for entry-exit problems are mixing processes, either in a rotatingdrum [53] or by use of impellers [15]. In these devices suspended particles enter and exit a neigh-bourhood of a moving wall. Since the fluid flow typically becomes ergodic from the boundaries, theentry-exit properties may play a key role for the mixing of suspended particles which move differentfrom the flow, in particular near the boundaries. Other systems of interest are fluidized beds [73]with wall effects becoming important in micro-fluidized beds [43]. Finally, dynamical entry-exitproblems should be of importance for the deposition of evaporating micro-droplets in the humanairways [75].
Consider the ordinary differential equation (ODE)d z d t =: z ′ = h ( z ( t ) , t ) , z ∈ R d , t ∈ R , (1)for d ≥ z ( t ) =: z ∈ K , where K is a compact subset of R d . We alwaysassume that h is smooth, i.e., C ∞ = C ∞ ( R d +1 , R d ); this assumption could be weakened but it willbe convenient for the applications we have in mind. Denote the flow associated with (1) by φ tt : z z ( t ) = z ( t ; t , z ) . (2)Let M ( t ) ⊂ R d be an m -dimensional smooth invariant manifold with φ tt ( M ( t )) ⊂ M ( t ) , ∀ t ∈ [ t , t ] =: I , for some t < t , with t , t ∈ I , and for t ≥ t ; we only consider the dynamics on the compact timeinterval I as we are interested in finite-time non-autonomous dynamics. Furthermore, we emphasizethat we do not require any particular type of invariant manifold such as normally hyperbolic [13] ornormally elliptic. However, one has to pick a relevant invariant manifold (e.g. a domain boundary, aspecial co-dimension one surface, an interface between two physically/biologically relevant regions).Here we do not discuss this choice as it does depend upon the application but we refer to the outlookin Section 6, where one goal would be to try our methods for different definitions of Lagrangiancoherent structures.Let T p M ( t ) and N p M ( t ) denote the tangent and normal spaces at p to M ( t ). Let γ = γ ( t ; t , x )denote a trajectory in M ( t ), i.e., we require γ ( t ) ⊂ M ( t ) for all t ∈ I . Suppose γ is smooth, i.e., γ ∈ C ∞ ( I × I × K , R d ). Then we define the function t
7→ F ( t ) = ˜ F ( γ ( t ; t , z )) , F : I → R , (3)where the choice for ˜ F respectively the construction of F will be discussed below. As yet, F isquite a general function as it just maps a time t to R . The idea is to associate the location of zerosof F with certain balancing properties of trajectories contained in M ( t ) , to relate the entrance andexit points for particles moving near M ( t ). We emphasize that we only need a reference trajectory γ ( t ) ⊂ M ( t ), not the entire manifold M ( t ) for our calculation.4 .1 Instantaneous Eigenvalues A naive first guess to study the dynamical properties near M ( t ) is to consider instantaneouseigenvalues. Consider the non-autonomous linear system Z ′ = [D z h ( γ ( t ; t , z ))] Z =: A ( t ; t , z ) Z, (4)for Z ∈ R d , A ( t ) = A ( t ; t , z ) ∈ R d × d . Let λ j = λ j ( t ; t , z ) for j ∈ { , , . . . , d } denote theeigenvalues of A ( t ). One guess, how to balance attraction and repulsion near M ( t ) is to consider F λ j ( t ) := Z tt ℜ ( λ j ( s )) d s, (5)where ℜ ( · ) denotes the real part of a complex number. The problem with just monitoring eigenval-ues is that they do not take into account direction, i.e., we just get undirected rates of growth anddecay. Furthermore, the eigenvalues for fixed times do usually not imply any stability statementsfor non-autonomous dynamical systems. Therefore, we expect that the direct use of eigenvalues isinsufficient in some many cases and we shall demonstrate these issues in Section 4.1. Another idea is to consider multiple time scale dynamics and use fast subsystem eigenvalues [35].Suppose (1) can be written in the standard fast-slow form ε d x d τ = ˜ f ( x, ˜ y, τ, ε ) , d˜ y d τ = g ( x, ˜ y, τ, ε ) , (6)where ( x, ˜ y ) ∈ R m +˜ n , the maps f : R m +˜ n +2 → R m and g : R m +˜ n +2 → R m are smooth (we againrequire C ∞ without further notice), let h = ( f, g ) ⊤ , and 0 < ε ≪
1. Note that we can makethe system autonomous by setting d τ d s = 1. This is equivalent to appending another slow variable˙˜ y ˜ n +1 = 1, so we shall now restrict to ε d x d s = ε ˙ x = f ( x, y, ε ) , d y d s = ˙ y = g ( x, y, ε ) , (7)where z := ( x, y ) ∈ R m + n for n = ˜ n + 1. The critical manifold of (7) is given by C := { ( x, y ) ∈ R m + n : f ( x, y,
0) = 0 } . (8)Suppose C is a smooth manifold which coincides with M ( t ) on the given domain for t , i.e., for y n . Recall that C is called normally hyperbolic if the eigenvalues of the matrix D x f ( p, ∈ R m × m have no zero real part [35]. Letting ε → f ( x, y, , d y d s = g ( x, y, , (9)for the dynamics on the critical manifold. As before, consider a trajectory γ ( s ) = γ ( s ; s , z )contained in C . Then we can consider the eigenvalues ρ j ( s ) for j ∈ { , , . . . , m } of the matrixD x f ( γ ( s ) ,
0) and define F ρ j ( t ) := Z tt ℜ ( ρ j ( s )) d s (10)5o study the local attraction and repulsion rates near C . Note that there is still a choice of theneighbourhood E of C , which can be selected to be of size O ( ε ) according to Fenichel’s Theorem [27]near the normally hyperbolic parts and has to be adapted to scalings of fast subsystem bifurcationpoints; see [35, Ch.7-8].The functions F ρ j are natural generalizations of the classical entry–exit maps developed firstin the context of delayed Hopf bifurcation [51, 52]; in fact, the delayed Hopf case may also havebuffer points so that the entry-exit map does not give the correct exit point by itself. Since we donot consider any oscillatory instabilities in this paper, we will not be further concerned with thisproblem but see Section 6. However, since rigorous proofs for branch points are available [65] (andreferences in [35]), it is important to have the fast-slow approach as a comparative benchmark.Furthermore, fast-slow techniques have been successfully applied to fluid dynamics of particles(see e.g. [63]). The potential disadvantage of the fast-slow construction is that the ODE (1) isusually not directly available in the form (7). This generically requires identifying a new parameter ε and moving the invariant manifold M ( t ) into a form where it becomes a slow manifold. Inaddition, not every system has a well-defined time-scale separation. Another classic concept to deal with spectral properties of dynamical systems defined on finite-time intervals are finite-time Lyapunov exponents (FTLEs); see e.g. [42, 19] for introductions. Theintuition of FTLEs is to consider a perturbation ζ ∈ R d to an initial condition for the flow φ tt ( z + ζ ) = φ tt ( z ) + [D φ tt ] ζ + O ( k ζ k ) , as k ζ k →
0, (11)where k · k always denotes the Euclidean norm, the linearized flow map isD φ tt : T R d → T R d , ( p, v ) ( φ tt ( p ) , [D z z ( t ; t , p )] v ) , (12)and T R d ≃ R d denotes the tangent bundle to R d . Hence, the expansion (11) shows that thelinearized flow map characterizes local separation and contraction properties. Then one may definethe maximum FTLE as l max ( t ; t , z ) = 1 | t − t | ln (cid:18) max k ζ k6 =0 k [D φ tt ] ζ kk ζ k (cid:19) . (13)In addition, one may also consider FTLEs associated to specific directions. This construction isconveniently expressed by considering the fundamental solution Φ( t ; t , z ) ∈ R d × d for Z ′ = A ( t ) Z and letting δ j = δ j ( t ; t , z ) denote the singular values of Φ( t ; t , z ). Then one may define theassociated FTLEs [9, 68] as l j = l j ( t ; t ; z ) := 1 t − t ln δ j ( t ; t , z ) . (14)Note carefully that FTLEs are dependent upon the final time t and do not constitute an infinitesimalnotion. Since FTLE are defined over two time points, t and t , one possibility would be to justdefine balance functions by F l j ( t ) := l j ( t ; t ; z ) . (15)However, as the definition already suggests, the computation of FTLEs is frequently impossible an-alytically and very challenging numerically. Furthermore, FTLEs do not directly take into accountthe existence of the invariant manifold M ( t ) as the calculation works for any trajectory in phasespace via the variational equation. 6 .4 NILE Exponent Here we briefly recall the theory of normal infinitesimal Lyapunov exponents (NILE) from [22] inthe context of (1). LetΠ tp : T p R d = T p M ( t ) ⊕ N p M ( t ) → N p M ( t ) , Π tp ( u, v ) = v (16)denote the natural projection onto the normal space at each p ∈ M ( t ). Then the NILE at p ∈ M ( t )is defined by σ ( p ; t ) := lim s → + s ln (cid:13)(cid:13)(cid:13)(cid:13) Π t + sφ t + st ( p ) [D φ t + st ] (cid:12)(cid:12)(cid:12) N p M ( t ) (cid:13)(cid:13)(cid:13)(cid:13) . (17)Although the definition may look complicated, σ ( p ; t ) is just the infinitesimal growth rate in thenormal direction to M ( t ) at a point p . Essentially, σ ( p ; t ) was designed in [22] to be a morecomputable measure of growth and decay rates for invariant manifolds in comparison to the classicalLyapunov-type numbers [12]. In particular, for a smooth unit normal vector field n ( p ; t ) to M ( t )one has σ ( p ; t ) = max n ( p,t ) ∈ N p M ( t ) n ( p, t ) ⊤ [D z h ( p, t )] n ( p, t ) (18)by [22, Thm.4]. Of course, one could use (18) as a definition instead of (17). In the case when M ( t ) can be written as a graph, one may consider z = ( x, y ) ∈ R m + n and the ODE x ′ ( t ) = f ( x ( t ) , y ( t ) , t ) ,y ′ ( t ) = g ( x ( t ) , y ( t ) , t ) . (19)Suppose M ( t ) = { x = m M ( y, t ) } and defineΓ( y, t ) := ∂f∂x ( m ( y, t ) , y, t ) − ∂m M ∂y ( y, t ) ∂g∂x ( m ( y, t ) , y, t ) , (20)then one may conclude [22, Thm.4] that σ ( z, t ) = λ max [Γ( y, t ) + Γ( y, t ) ⊤ ] / , (21)where λ max [ · ] is the largest eigenvalue of a symmetric matrix. Considering the fact that σ ( p, t ) < M ( t ) while σ ( p, t ) > F σ ( t ) := Z tt σ ( γ ( s ; t , z ) , s ) d s. (22)The advantage of NILE is that it does measure the generic (i.e., typical) maximum growth andminimum decay rates for repelling, respectively attracting, manifolds. Furthermore, it does takeinto account the geometry by focusing on the normal direction to the manifold M ( t ), which is theone relevant for entry and exit of neighbourhoods of M ( t ). We again remark that we have to selectthe neighbourhood E ( t ) so that the linearization approximation in (22) is sufficiently accurate. Incomparison to the fast-slow definition (10) where natural scales in ε appear, the scales for NILEare related to uniform bounds on the normal directions and the requirement that the leading linearnormal direction is not dominated by nonlinear terms in E ( t ). As a practical strategy it seemssafest to start with a neighbourhood E ( t ) of very small volume, where the approximation must bevalid and gradually increase the neighbourhood. In particular, one may extend it until the linearapproximation differs on several small compact test subsets from the full flow for a given tolerance.7quation (21) in conjunction with [9] shows that the NILE concept is related to D-hyperbolicity [4],i.e., the generalization of the existence of an exponential dichotomy for finite-time intervals. Thereby,NILE also relates to other notions of spectra and hyperbolicity for non-autonomous systems on finitetime intervals as discussed in [3, 9, 30]. In this context one usually obtains spectral intervals [59, 30]and one could aim to use these intervals to define new balance functions. However, for our pur-poses we aim to find a concept, which is analytically and conceptually simple in low-dimensionalexamples and stays numerically computable for more complicated flows. We shall see in severalexamples that a balance function based on NILE satisfies these requirements quite well. We develop some basic general theory for the class of balance functions F defined above. Eachfunction depends not only on its argument T ∈ R but also on the inital time t , and the initialpoint z . Therefore, we also consider G υ : I × I × K → R , G υ ( t, t , z ) = F υ ( t ) (23)using the G -notation to emphasize this dependence and with υ ∈ { λ j , ρ j , l j , σ } ; recall that λ j refers to instantaneous eigenvalues, ρ j to the fast-slow case, l j to finite-time Lyapunov exponents(FTLEs), and σ to the normal infinitesimal Lyapunov exponent (NILE). Proposition 3.1. (basic properties) The following results hold: • G υ ( t , t , z ) = 0 , for υ ∈ { λ j , ρ j , σ } ; • lim t → t G l j ( t, t , z ) exists; • G υ is bounded for υ ∈ { λ j , ρ j , σ, l j } .Proof. The first statement is trivial due to the integral definition of the balance function for { λ j , ρ j , σ } . For the second statement, it suffices to study t = 0 and the remaining cases willfollow by a shift. Since γ ∈ C it follows that A ( t ) is C in t so A ( t ) = A + tA + o ( t ) by Taylor’sTheorem. For the fundamental matrix we haveΦ( t, , z ) = e R t A ( s ) d s = e t [ A + o ( t )] . (24)So the singular values δ j of Φ( t, , z ) satisfy δ j = e t [ a + o (1)] and the second result follows. The laststatement can be deduced from γ ∈ C ; for example, consider υ = λ j , then A ( t ) is C and definedon the compact set I so the eigenvalues are bounded on I , and so are their real parts. Integratinga bounded function over a compact set again yields a bounded function. Boundedness with respectto the initial point follows from the compactness of K . The other cases are equally easy.We remark that the last proof only used continuous differentiability. Recall that we are primar-ily interested in zeros of balance functions t
7→ G ( t, t , z ) = F ( t ). The existence of zeros is a globaldynamical problem and has to be considered on a case-by-case basis. However, the non-degeneracyof a zero is a local property. Non-degeneracy of a nontrivial zero T > t yields a transversality con-dition for the motion near the exit point, i.e., we expect a true exit and not just sliding near the exitboundary. Therefore, we have to study differentiability properties of G , particularly with respectto t . To state the next result, we say that a family of matrices is NIC (no-infinite-contacts) [57]if none of its eigenvalues meet with an infinite order of contacts when t is varied, i.e., roots of thecharacteristic polynomial have well-defined finite multiplicities for all t .8 roposition 3.2. (differentiability) Consider t , t in the interior of I and z ∈ K . Then thefollowing hold:(D1) t
7→ G λ j ( t ; t , z ) is C ; if A ( s ; t , z ) is normal and NIC for all s ∈ I then ( t, t , z ) λ j ( t, t , z ) is C ∞ ;(D2) t
7→ G ρ j ( t ; t , z ) is C ; if D x f ( γ ( s ; t , z ) , is normal and NIC for all s ∈ I then ( t, t , z ) ρ j ( t, t , z ) is C ∞ ;(D3) if AA ⊤ is NIC then ( t, t , z )
7→ G l j ( t ; t , z ) is C ∞ ;(D4) suppose M can be written globally as a graph of a C ∞ -function, and Γ( y, s ) is NIC for all s ∈ I then ( t, t , z )
7→ G σ ( t ; t , z ) is C ∞ ;Proof. Regarding (D1), differentiability in t is immediate from the Leibniz integral formula, whichyields a continuous derivative. For smoothness in the initial data, consider first t . Note that A = A ( s ; t , z ) is a C ∞ function of t by assumption on the smoothness of the vector field andthe assumption that γ ∈ C ∞ . Since A is normal, NIC, and smooth, it follows that its eigenvalues λ j = λ j ( s ; t , z ) are C ∞ as functions of t [57, Thm.7.8]. Therefore, the corresponding real partsare C ∞ as well. Using the Leibniz integral formula, we obtain smoothness with respect to t . Thesame argument can now be applied to each of the coordinates of z and (D1) follows. (D2) isimmediate as the proof from (D1) carries over. For (D3), the first observation is that AA ⊤ is anormal matrix. Since AA ⊤ is also NIC, the singular values are C ∞ ; now we can just use the sameTaylor expansion argument as in the proof of Proposition 3.1 to get that l j is C ∞ . For (D4), notethat since M can be written as a C ∞ -graph, it follows that Γ( y, t ) is C ∞ , so using formula (21)the last result follows by the same arguments as (D1)-(D2).The differentiability results can be improved but selecting eigenvalue-type quantities smoothly isan extremely technical topic and many special cases may occur [57]. Here we are mainly interestedin the dependence of zeros of F on the input data. A direct application of the implicit functiontheorem yields: Corollary 3.3.
Suppose h ∈ C ∞ and γ ∈ C ∞ for each trajectory and the assumptions in (D1)-(D4) hold. Consider ( t, t , z )
7→ G υ ( t, t , z ) for υ ∈ { λ j , ρ j , σ, l j } and suppose t
7→ F υ ( t ) has annon-degenerate root at T , i.e., F ′ ( T ) = 0 . Then there exists a neighbourhood of ( T, t , z ) and alocally unique continuous one-dimensional curve of solutions to G υ ( t, t , z ) . The last result is not really surprising. It just states that if we find a zero of the balancefunction, and this zero is isolated, then the C -dependence of the balance function on its inputdata guarantees that the zero perturbs uniquely. However, for applications it is of paramountimportance to have at least some computable conditions for robustness such as F ′ ( T ) = 0. Remark:
Regularity results are also relevant to study the dependence upon initial conditionssuch as sets of the form S t := { z ∈ K : F ( t ) = 0 } , (25) i.e., considering the evolution of co-dimension one submanifolds of initial conditions. Fixing aneighbourhood size of M ( t ) , makes S t depend only upon z . However, if M ( t ) interacts withanother invariant manifold then an initially smooth connected set S t may break up. Furthermore,since S t are sets imposing balance conditions, it would be natural ask, how this relates to Lagrangiancoherent structures and we aim to pursue this analysis in future work. .6 Perturbation Expansion Having reduced the entry-exit problem to a root-finding problem, we remark that it is now possibleto apply standard perturbation techniques if suitable boundedness and regularity conditions aresatisfied. For example, consider the case when we have to solve0 = F ( t ; δ ) (26)for some small system parameter δ ≥
0. If F ( T ; 0) = 0 and ∂ F ∂t ( T ; 0) = 0 then the implicit functiontheorem yields that we can locally solve (26) via T = T ( δ ) and formally can make the perturbationansatz: T ( δ ) = T + δT + δ T + · · · , T j ∈ R , < δ ≪ . (27)Furthermore, if (26) depends upon a trajectory γ ( t ) ⊂ M ( t ) which has a perturbation expansion γ ( t ) = γ ( t ) + δγ ( t ) + δ γ ( t ) + · · · (28)then we can just plug all terms into (26) and formally expand in δ , collect terms of different orders,and aim to solve the problem perturbatively if the perturbation problem is regular. Here we shall consider two basic models that represent elementary flows to benchmark the previ-ously introduced concepts. ’Passive’ refers to the fact that particles inside these flows are viewedas passively transported while Section 5 presents an ’active’ case, where there is a fully-coupledparticle-fluid interaction.
We start with a basic example, which is nevertheless quite insightful as it is relatively easy from acomputational perspective and demonstrates the basic features of balance functions. Consider theODE z ′ = − ( z − β ) ,z ′ = z (1 − e − αz ) , (29)where ( z , z ) ∈ R × [0 , + ∞ ), and α, β > z , z ) = (0 , β ). M = { y = 0 } is a time-independent invariant manifold with dynamics z ′ = β . We are interested in the entry-exit relation of a point particle getting passively transported by the flow (29) in and out of aneighbourhood E ( χ ) := { ( z , z ) ∈ R : z ∈ [0 , χ ] } . (30)Assume that χ > M . We expect a symmetric entry-exit relation. Indeed,given the initial point z = ( − b, χ ) ∈ ∂ E ( χ ), we know already that the first exit of the point isgiven by z T = ( b, χ ) due to the symmetry( z , z , t ) ( − z , z , − t ) (31)of (29). We discard this fact temporarily and apply the different functions F . We start by justnaively considering the concept of instantaneous eigenvalues discussed in Section 3.1.10 z z Figure 2: Trajectories of the ODE (29) for parameters β = 1, α = 2. The invariant manifold M = { z = 0 } (thick black line) is a fixed wall with neighbourhood E ( χ ) (indicated by a dashedblack line). The center of rotation at (0 , β ) is marked by a black dot. A typical relationship betweenentry (grey dot) and exit (grey circle) of a particle transported by the flow is indicated as well. Proposition 4.1.
There exists an open set of parameters α, β, b > such that the exit computedby balancing F λ j is correct, while there also exists an open set of parameters where balancing F λ j gives an incorrect exit point.Proof. As a trajectory in M , we must take γ ( t ; t , z ) = ( β ( t − t ) − b, ⊤ . (32)Without loss of generality consider t = 0. Linearization yields that A = A ( t ; t , z ) = (cid:18) − α ( βt − b ) (cid:19) . (33)The two eigenvalues of A are λ ± = (cid:16) α ( βt − b ) ± p [ α ( βt − b )] − (cid:17) . Therefore, consider F λ ± ( T ) = Z T ℜ ( λ − ( s )) d s, = Z bααβ ∧ T α βs − b ) d s ± Z T bααβ ∧ T (cid:16) α ( βt − b ) − p [ α ( βt − b )] − (cid:17) d s, (34)where bααβ ∧ T = min { bααβ , T } . The second integral can be evaluated in certain cases but is alengthy expression. It is interesting to just consider certain cases. Suppose bα <
4, then F λ ± hasa zero at T = bβ since bβ < bααβ . This yields the correct exit point since β bβ − b = b . This provesthe first part of the proposition. If bα > β = 1, α = 2 and b = 3, then we have F λ ± ( T ) = − ± Z T t − − p [( t − − s (35)11he last integral can be evaluated and one checks that F λ ± (6) = 0. However, γ (6; 0 , ( − , − ,
0) = (3 ,
0) is the correct exit point. A direct application of Corollary 3.3 yields an open setof parameters where F λ j yields an incorrect exit point.Since monitoring the real parts of instantaneous eigenvalues and requiring an integrated balancefunction F λ j does not yield the correct result in relevant cases, we discard this option from nowon. Another option considered in Section 3.2 is to scale the problem differently. Proposition 4.2.
There exists a scaling of (29) converting it to the standard fast-slow form (7) .Furthermore, balancing the function F ρ yields the correct exit point.Proof. Consider the scaling ( z , z , t ) = ( y/ √ ε, x, s/ √ ε ) and observe that (29) then becomes afast-slow system ε ˙ x = y (1 − e − αx ) , ˙ y = − ( x − β ) . (36)The critical manifold is C = { y = 0 } ∪ { x = 0 } . Note that C is not a smooth manifold but thesubset { x = 0 } = M is an invariant manifold for (36) for any ε > M here.The slow subsystem is given by ˙ y = β so γ ( s ) = (0 , β ( s − s ) − b ). In the notation of Section 3.2we have the linearization of the fast vector field yields just one eigenvalue so that ρ = D x f ( γ,
0) = α [ β ( s − s ) − b ] . (37)Therefore, the balance function can be calculated, say for s = 0, F ρ ( T ) = α Z T βs − b d s = α (cid:20) β T − bT (cid:21) = αT (cid:20) β T − b (cid:21) . (38)Now the balance condition F ρ ( T ) = 0 provides the correct answer that the exit point is given, viathe zero T = 2 b/β , as ( x, y ) = (0 , b ).Furthermore, we easily check that the zero of the balance function is indeed isolated since F ′ ρ (2 b ) = α (2 b − b ) = αb > α, b > Z ′ = A ( T ) Z , we know that a fundamental solution can be writtenas Φ( t ; t , z ) = exp (cid:20)Z tt A ( r ) d r (cid:21) . (39)Working out the integral is easy but the algebraic form of matrix exponential is extremely lengthyin the general case. However, observe that Z tt A ( r ) d r = (cid:18) − tt t ( βt − b ) α (cid:19) (40)So if t = 2 b/β then we easily find that the singular values of Φ(2 b/a, , ( − b, F l (2 b/β ) = 0 = F l (2 b/β ) yielding the correct balance time forthis case. Proposition 4.3.
Balancing F l j yields the correct entry-exit relationship for (29) . t corresponds actually to a full solution of the non-autonomous linear system along theinvariant manifold. Hence, it seems worthwhile to search for simpler balance functions that takeinto account the existence of M and its geometry. A natural option seems to be the NILE asdiscussed in Section 3.4. Proposition 4.4.
Balancing F σ yields the correct entry-exit relationship for (29) .Proof. Consider (32) for t = 0 and the invariant manifold M = { z = 0 } for (29). A unit normalvector field is given by n ( p, t ) = (0 , ⊤ and formula (18) yields σ ( γ, t ) = (0 , (cid:18) − α ( βt − b ) (cid:19) (cid:18) (cid:19) = α ( βt − b ) . Hence, the balance function F σ also yields the correct result for the exit point by the same calcu-lation as in (38).No preliminary scaling by ε was necessary for the NILE case nor was it necessary to solve thefull non-autonomous system. Hence, the concept seems to be well-suited to take into account thegeometry of the invariant manifold. Furthermore, we may calculate eigenvalues directly from thelinearized problem without the need to scale by a small parameter. However, NILE is a singlerate so one should keep in mind that it would be useful to generalize the case from a unit normalvector field to M to other transverse frames when projection operators depend upon time andposition; cf. [22, p.614-615]. However, generically the entry (resp. departure) speed is governed bythe largest normal rates in the attracting (resp. repelling) regime so NILE already captures whatwe are interested in for most applications. An axisymmetric liquid bridge of length d , stabilized by surface tension, can be established betweentwo parallel and coaxial wetted disks of cylindrical rods of radius R . If the disks are heated differ-entially and the aspect ratio Γ = d/R is of order O (1), temperature gradients along the capillaryinterface induce a toroidal vortex via the thermocapillary effect [67], which can become modulatedazimuthally for large imposed temperature differences. Muldoon and Kuhlmann [49] proposed aclosed-form approximation to the three-dimensional time-dependent Navier–Stokes flow in a cylin-drical liquid bridge aiming at modeling the rapid de-mixing of particles into curious accumulationstructures found experimentally (see e.g. [66]); see also Section 2 for further background. Usinga separation ansatz and free-slip boundary conditions on the supporting disks, the axisymmetricpart of the flow field was represented by harmonic functions. The algebraic radial dependence ofthe axisymmetric flow is then dictated by the incompressibility constraint. The remaining integra-tion constants were used to fit the strength and the radial shape of the vortex to the numericallyobtained Navier–Stokes solution. The axisymmetric steady part of the Kuhlmann–Muldoon modelflow is symmetric with respect to a midplane { z = 0 } . Here we extend this flow by an antisym-metric term which may take care of the inertia-induced asymmetry of the real flow. This also aimsto test balance functions in a setup without symmetry. Let the fluid be confined to the cylindricalvolume ( z , φ, z ) ∈ [0 , / Γ] × [0 , π ] × [ − . , .
5] and discard the second component. Then we definethe planar model flow for an aspect ratio Γ = 2 / z ′ = πz η (1 − z )[sin( πz ) − α s cos(2 πz )] ,z ′ = [( η + 1) z η − − ( η + 2) z η ][cos( πz ) + α s sin(2 πz )] , (41)13here z ∈ (0 , ] and z ∈ ( − , ). The parameter α s controls the asymmetry of the flow, and theexponent η ∈ [4 ,
5] determines the radial shape of the vortex. Figure 3 shows examples of the flowfield. z z . . . − . (a) z z . . . − . (b) Figure 3: Streamlines and vector field of the model flow (41) for η = 4 .
74 and α s = 0 . α s = 0 . M = { z = , z ∈ ( − , ) } (capillarysurface). Note that the flow on M can be calculated from z ′ = − (cid:18) (cid:19) η − [cos( πz ) + α s sin(2 πz )] . (42)Unfortunately, there does not seem to be a closed form solution to the ODE (42). As before, let γ denote a trajectory contained inside M . The linearization A ( t ) = D h ( γ ( t )) with t = 0 startingfrom initial data ( z , z ) = ( − ξ, b ) along γ can be calculated (cid:0) (cid:1) η − (2 α s cos(2 πγ ) − sin( πγ )) 0 − (cid:0) (cid:1) − η (2 η + 1)(cos( πγ ) + A sin(2 πγ )) − (cid:0) (cid:1) − η π (2 α s cos(2 πγ ) − sin( πγ )) ! The full analytical computation of FTLEs is not possible, hence we start by considering NILE.Since a unit normal vector field to M is just given by n ( p, t ) = (0 , ⊤ , we can easily calculate(0 , A (cid:18) (cid:19) = − (cid:18) (cid:19) − η π (2 α s cos(2 πγ ) − sin( πγ )) (43)Now one can actually numerically integrate (42) and then insert the solution γ into the integral F σ ( t ) = − (cid:18) (cid:19) − η π Z t α s cos(2 πγ ) − sin( πγ ) d s. (44)We already see that if we are just interested in zeros of F σ we may discard the prefactorin (44) as η ∈ [4 , α s , η on the exit time T . We observe a very clear linearrelationship between the parameters and the exit time T . The initial condition has been fixed to( , ). If α s , η are decreased, we observe that this induces a quicker escape in terms of escape times.14 TT α s η (a) (b)Figure 4: Computation for the Kuhlmann–Muldoon model relating the initial time t = 0 to a finalexit time T using the integral formula (44) finding T by requiring F σ ( T ) = 0. (a) The parameter η = 4 .
74 and the initial condition ( , ) are fixed and the asymmetry parameter is varied. Thecomputed exit times are marked as dots; a linear fit to the computed points is shown. (b) Identicalsetup as in (a) except that the parameter α s = is fixed and the parameter η is varied to computethe exit time.Figure 5 converts the escape times into escape points by integrating the one-dimensional ODE (42).Interestingly the dependence on the two parameters now shows a completely different behaviour.For the asymmetry parameter α s , a nonlinear behaviour is observed while for η , we observe a similarexit point regardless of the parameter, i.e., the changing escape time is compensated by a differentspeed on M . z ( T ) z ( T ) α s η (a) (b)Figure 5: Computation for the Kuhlmann–Muldoon model relating the initial time t = 0 to a finalexit time T using the integral formula (44). In contrast to Figure 4 we show the exit point z ( T )on the vertical axis obtained from numerically integrating (42) up to time T . (a) The parameter η = 4 .
74 and the initial condition ( , ) are fixed and the asymmetry parameter is varied. Thecomputed exit points are marked as dots. (b) Identical setup as in (a) except that the parameter α s = is fixed and the parameter η is varied to compute the exit point.The analysis shows that it is very convenient to use balance functions numerically and todetermine the influence of different parameters on particle motion near M .15 z z L a Figure 6: Sketch of the two-dimensional shear–stress-driven cavity seeded with a particle of radius a = 0 . · L . Remark:
Another way to arrive at the entry-exit map is to scale the original ODE (41) by firstmoving M to the manifold { z = 0 , z ∈ ( − . , . } and then scaling z by a suitable power of asmall parameter ε to make the z -variable a fast x -variable. However, the problem already showsthat scaling becomes more and more involved once the manifold M is nontrivial. Particle-laden flows are a class of multiphase flows in which there is a continuously connected fluid-phase and a dispersed particle phase which do not mix. When the particles scales are very smallcompared to the fluid scales, neglecting the feedback effect of the particles on the fluid flow can be agood approximation in numerically calculating each particle trajectory. In this case we talk aboutone-way coupling simulations. However, if the particle passes very close to a wall/free-surface,the particle and the lubrication-gap scales must be solved and the feedback effect of the particleon the fluid phase is essential [60, 32], i.e., the non-autonomous dynamics really differs from theautonomous case. Simulating both phases can be very expensive computationally. Therefore, un-derstanding the physical mechanisms which play the main role in the particle–boundary interactioncontributes to an accurate modelling of this phenomenon. Taking into account lubrication effectscan lead to a significant improvement in predicting the particle trajectories. Hence, this setup isan excellent test case for the framework of balance functions presented above.Consider a two-dimensional square cavity filled with an incompressible Newtonian liquid ofdensity ρ f and kinematic viscosity ν as shown in Figure 6; we denote the horizontal and verticalcoordinates by z and z respectively similar to our notation above and set z = ( z , z ) ⊤ . Thecavity of linear length L is open from above and the fluid flow is driven by a constant shear stress, τ , in z -direction. In the limit of asymptotically large surface tension the interface is flat. The flowin the shear–stress-driven cavity can be modelled using the Navier–Stokes equations and employinga viscous scaling ˆ u = νL u, ˆ z = Lx, ˆ t = L ν t, ˆ p = ρ f ν L p, (45)where ’hats’ indicate the dimensional variables and ’no-hats’ the non-dimensional ones, one obtains16he non-dimensional Navier–Stokes system ∇ · u = 0 , (46a)( ∂ t + u · ∇ ) u = −∇ p + ∇ u, (46b)where u = u ( z, t ) ∈ R and p = p ( z, t ) ∈ R represent the flow velocity and pressure fields,respectively. The flow field must satisfy the no-slip and constant-stress boundary conditionson { z = ± / } , u = 0 = u ; (47a)on { z = − / } , u = 0 = u ; (47b)on { z = 1 / } , ∂ z u = − Re , u = 0; (47c)where Re is the Reynolds number Re = τ L ρ f ν . (48)and in the following it will be set equal to 1000. Considering now the cavity seeded with one rigidcircular particle of radius a , rigid body motion equations will define its trajectory in the fluid flow F = M d V d t , (49a) T = I dΩd t (49b)where F and T are the force and torque acting on the particle, V and Ω denote its translationaland rotational velocities, and M and I are the mass and the inertia tensor of the particle. Thecoupling between the particle and the fluid motions results from the no-slip and no-penetrationconditions on the surface of the particle. In the following the case of particle radius a = 0 . L is considered for a particle-to-fluid density ratio ̺ = ρ p /ρ f = 2. The initial particle position isset to ( z , z ) = (0 . , . × − in the z -direction. Along the z -direction a stretching function proportional to z (2 .
4) is adopted with minimum side length of9 × − for the elements in contact with the free-surface. The time step is fixed to 1 × − [61].Having the full simulation data of the flow available, we want to check the predictive capabilitiesof the NILE balance condition F σ = 0 by computing F σ from the data. To make the problem harder,suppose we did not even know a reference trajectory inside the invariant manifold/free-surface M = { ( z , z ) ∈ R : z ∈ [ − . , . , z = 0 . } . Suppose we can only track the particle itself and have the local flow normal to M . We project theparticle trajectory ˜ γ ( t ) along the direction normal to the wall γ ( t ) := (0 , ˜ γ ( t )) ⊤ ⇒ γ ( t ) ⊂ M . This means we can calculate F σ ( t ) (using t = 0) at fixed times according to formula (22). Clearlythis predictor is locally linearized and only takes into account the flow near M . To have a com-parison, we consider a fully nonlinear balance function that takes into account the local velocity of17 F σ F σ F σ F σ F v F v F v F v z z z z z z z z tt tt tt tt Figure 7: Each row of the figure corresponds to one time snapshot of the full numerical simulationat time t , i.e., we used data available in [0 , t ]. The left column shows the flow fields including thefree-surface for z = 0 . z = 0 . E := { z ∈ [ − . , . , z ∈ [0 . , . } where we calculate theintegrals for the particle; red indicates the calculation is not active while green indicates an activecalculation of the NILE predictor as well as the full nonlinear predictor. The two rightmost columnsshow the predictors (multiplied by a factor of 10 to increase the scale). The dashed blue line marksthe zeros, i.e., we look for the first non-trivial zero crossing. The vertical dashed magenta line in thelast time snapshot shows the true crossing of the particle from the active region past the thresholdat z = 0 . F v ( t ) := Z t v (˜ γ ( s )) d s, (50)18here v (˜ γ ( s )) is the actual particle velocity normal to the wall at each time point s . To compute F v and F σ , we only consider a region E := { z ∈ [0 . , . } as shown in Figure 7, i.e., outside of E , thereis no contribution to the balance functions. Figure 7 shows the main results for the balance functionsat four different time snapshots. F σ based upon NILE performs extremely well considering the factthat only local linearized flow near M is taken into account. It slightly underestimates the trueexit point while using F v slightly overestimates it. In particular, there seems to be no significantperformance difference between the two balance functions, although one uses locally linear and theother global fully nonlinear information.It is important to point out that one may even try to use F σ as a predictor. Consider the thirdrow in Figure 7 at time t = 0 . F σ vanishes.This adds to the potential flexibility of the balance function framework discussed here for manypractical applications. In this paper we have proposed a framework to address the entry-exit relationship for trajectories indynamical systems with a strong focus on applications in fluid dynamics. Particles travelling nearinvariant manifolds have a tremendous relevance for phenomena in fluids. We have demonstratedthat the concept of balance functions is well-suited to understand where particles are expected toexit neighbourhoods of time-dependent invariant structures in fluid dynamics. Two model problemshave been used to determine that NILE balance functions provide a very efficient way to mergetheoretical precision as well as practical computation for particle trajectories. For the Kuhlmann–Muldoon model, we showed how parameter dependencies can be uncovered using NILE balancefunctions. Then we tested the concept on fully-resolved direct numerical simulation of a particle ina shear–stress-driven cavity. Again we obtained excellent agreement with the abstract theory. Inparticular, local data near the invariant manifold and a projected particle trajectory were sufficientto gain insight in the particle exit point and provided a practical strategy to predict it based uponthe current value and/or extrapolation of the balance function.We stress that we only aimed here at demonstrating the main technical elements of balancefunctions, and tried to show their applicability to key practical issues in fluid dynamics. Severalopen questions remain after this study. On the one hand, further mathematical issues arise. Forexample, what happens for balance functions in non-autonomous systems in the case of oscillatoryinstabilities, such as delayed Hopf bifurcations [51, 35]. The question of random coefficients andchange of stability points [36] as well as stochastic forcing [5, 41] are problems to be tackled in futurework. Improving the link to the theory of Lagrangian coherent structures is another interestingdirection to pursue so is the connection to purely set-valued approaches to dynamical systems.Furthermore, one might check whether balance functions can equivalently be defined using theKoopman operator, which is a strategy that already worked in other contexts for isochrons andisostable sets [45].With balance functions, many practical questions can now be tackled. The last example of theNavier-Stokes equation shows that working directly with data can be successful. To observe theresults presented here in laboratory and field experiments would be of interest. Applications todata analysis in oceanography and engineering are also conceivable as measurements near invariantmanifolds may potentially be easier to obtain in comparison to full monitoring of flow fields, e.g.,in the context of flows in rivers as well as in flows near the coastline.19 eferences [1] B. Aulbach, M. Rasmussen, and S. Siegmund. Invariant manifolds as pullback attractors ofnonautonomous differential equations.
Discr. Cont. Dyn. Syst. , 15(2):579–597, 2006.[2] S.M. Baer, T. Erneux, and J. Rinzel. The slow passage through a Hopf bifurcation: delay,memory effects, and resonance.
SIAM J. Appl. Math. , 49(1):55–71, 1989.[3] A. Berger. On finite-time hyperbolicity.
Commun. Pure Appl. Anal. , 10(3):963–981, 2011.[4] A. Berger, T.S. Doan, and S. Siegmund. A definition of spectrum for differential equations onfinite time.
J. Differential Equat. , 246:1098–1118, 2009.[5] N. Berglund and B. Gentz.
Noise-Induced Phenomena in Slow-Fast Dynamical Systems .Springer, 2006.[6] M. Branicki and S. Wiggins. An adaptive method for computing invariant manifolds in non-autonomous, three-dimensional dynamical systems.
Physica D , 238(15):1625–1657, 2009.[7] C. Clarke and B. Carswell.
Principles of Astrophysical Fluid Dynamics . CUP, 2007.[8] R.G. Cox and S.K. Matthews. The lateral migration of solid particles in a laminar flow neara plane.
Int. J. Multiphase Flow , 3(3):201–222, 1977.[9] T.S. Doan, D. Karrasch, T.Y. Nguyen, and S. Siegmund. A unified approach to finite-timehyperbolicity which extends finite-time lyapunov exponents.
J. Differential Equat. , 252:5535–5554, 2012.[10] L.H. Duc and S. Siegmund. Hyperbolicity and invariant manifolds for planar nonautonomoussystems on finite time intervals.
Int. J. Bif. Chaos , 18(3):641–674, 2008.[11] F. Dumortier and R. Roussarie.
Canard Cycles and Center Manifolds , volume 121 of
MemoirsAmer. Math. Soc.
AMS, 1996.[12] N. Fenichel. Persistence and smoothness of invariant manifolds for flows.
Indiana U. Math.J. , 21:193–225, 1971.[13] N. Fenichel. Geometric singular perturbation theory for ordinary differential equations.
J.Differential Equat. , 31:53–98, 1979.[14] G. Froyland, K. Padberg, M.H. England, and A.M. Treguier. Detection of coherent oceanicstructures via transfer operators.
Phys. Rev. Lett. , 98(22):224503, 2007.[15] P.R. Gogate, A.C.M Beenackers, and A.B. Pandit. Multiple-impeller systems with a specialemphasis on bioreactors: a critical review.
Biochem. Eng. J. , 6(2):109–144, 2000.[16] M.A. Green, C.W. Rowley, and G. Haller. Detection of Lagrangian coherent structures inthree-dimensional turbulence.
J. Fluid Mech. , 572:111–120, 2007.[17] J. Guckenheimer and P. Holmes.
Nonlinear Oscillations, Dynamical Systems, and Bifurcationsof Vector Fields . Springer, New York, NY, 1983.[18] S. Haber. A spherical particle moving slowly in a fluid with a radially varying viscosity.
SIAMJ. Appl. Math. , 67(1):279–304, 2006.[19] G. Haller. Distinguished material surfaces and coherent structures in three-dimensional fluidflows.
Physica D , 149(4):248–277, 2001. 2020] G. Haller. A variational theory of hyperbolic Lagrangian coherent structures.
Physica D ,240(7):574–598, 2011.[21] G. Haller. Lagrangian coherent structures.
Ann. Rev. Fluid Mech. , 47:137–162, 2015.[22] G. Haller and T. Sapsis. Localized instability and attraction along invariant manifolds.
SIAMJ. Appl. Dyn. Syst. , 9(2):611–633, 2010.[23] G. Haller and G. Yuan. Lagrangian coherent structures and mixing in two-dimensional turbu-lence.
Physica D , 147(3):352–370, 2000.[24] J.S. Hesthaven and T. Warburton.
Nodal Discontinuous Galerkin Methods: Algorithms, Anal-ysis, and Applications . Springer, 2007.[25] M.W. Hirsch, C.C. Pugh, and M. Shub.
Invariant Manifolds . Springer, 1977.[26] E. Hofmann and H.C. Kuhlmann. Particle accumulation on periodic orbits by repeated freesurface collisions.
Phys. Fluids , 23:0721106, 2011.[27] C.K.R.T. Jones. Geometric singular perturbation theory. In
Dynamical Systems (MontecatiniTerme, 1994) , volume 1609 of
Lect. Notes Math. , pages 44–118. Springer, 1995.[28] G.G. Joseph, R. Zenit, M.L. Hunt, and A.M. Rosenwinkel. Particle-wall collisions in a viscousfluid.
J. Fluid Mech. , 433:329–346, 2001.[29] G.E. Karniadakis, M. Israeli, and S.A. Orszag. High-order splitting methods for the incom-pressible Navier-Stokes equations.
J. Comput. Phys. , 97(2):414–443, 1991.[30] D. Karrasch. Linearization of hyperbolic finite-time processes.
J. Differential Equat. , 254:256–282, 2013.[31] P.E. Kloeden and M. Rasmussen.
Nonautonomous Dynamical Systems . AMS, 2011.[32] P. Kosinski, A. Kosinska, and A.C. Hoffmann. Simulation of solid particles behaviour in adriven cavity flow.
Powder Technol. , 191(3):327–339, 2009.[33] M. Krupa and P. Szmolyan. Relaxation oscillation and canard explosion.
J. Differential Equat. ,174:312–368, 2001.[34] C. Kuehn. Normal hyperbolicity and unbounded critical manifolds.
Nonlinearity , 27(6):1351–1366, 2014.[35] C. Kuehn.
Multiple Time Scale Dynamics . Springer, 2015. 814 pp.[36] C. Kuehn. Uncertainty transformation via Hopf bifurcation in fast-slow systems. arXiv:1512.03002 , pages 1–19, 2015.[37] H.C. Kuhlmann.
Thermocapillary Convection in Models of Crystal Growth , volume 152 of
Springer Tracts in Modern Physics . Springer, Berlin, Heidelberg, 1999.[38] H.C. Kuhlmann and F.H. Muldoon. Comment on “Ordering of small particles in one-dimensional coherent structures by time-periodic flows”.
Phys. Rev. Lett. , 108:249401, 2012.[39] H.C. Kuhlmann and F.H. Muldoon. Particle-accumulation structures in periodic free-surfaceflows: Inertia versus surface collisions.
Phys. Rev. E , 85:046310, 2012.[40] H.C. Kuhlmann and F.H. Muldoon. Comment on “Synchronization of finite-size particlesby a traveling wave in a cylindrical flow” [Phys. Fluids 25, 092108 (2013)].
Phys. Fluids ,26(9):099101, 2014. 2141] R. Kuske. Probability densities for noisy delay bifurcation.
J. Stat. Phys. , 96(3):797–816, 1999.[42] F. Lekien and S.D. Ross. The computation of finite-time Lyapunov exponents on unstructuredmeshes and for non-Euclidean manifolds.
Chaos , 20(1):017505, 2010.[43] X. Liu, G. Xu, and S. Gao. Micro fluidized beds: Wall effect and operability.
Chem. Eng. J. ,137(2):302–307, 2008.[44] X. Luo, M.R. Maxey, and G.E. Karniadakis. Smoothed profile method for particulate flows:error analysis and simulations.
J. Comput. Phys. , 228(5):1750–1769, 2009.[45] A. Mauroy, I. Mezic, and J. Moehlis. Isostables, isochrons, and Koopman spectrum for theaction-angle reduction of stable fixed point dynamics.
Physica D , 261:19–30, 2013.[46] M.R. Maxey and J.J. Riley. Equation of motion for a small rigid sphere in a nonuniform flow.
Phys. Fluids , 26(4):883–889, 1983.[47] I.N. McCave. Size spectra and aggregation of suspended particles in the deep ocean.
Deep SeaRes. A , 31(4):329–352, 1984.[48] D.E. Melnikov, D.O. Pushkin, and V.M. Shevtsova. Synchronization of finite-size particles bya traveling wave in a cylindrical flow.
Phys. Fluids , 25(9):092108, 2013.[49] F.H. Muldoon and H.C. Kuhlmann. Coherent particulate structures by boundary interactionof small particles in confined periodic flows.
Physica D , 253:40–65, 2013.[50] Y. Nakayama and R. Yamamoto. Simulation method to resolve hydrodynamic interactions incolloidal dispersions.
Phys. Rev. E , 71:036707, 2005.[51] A.I. Neishtadt. Persistence of stability loss for dynamical bifurcations. I.
Differential EquationsTranslations , 23:1385–1391, 1987.[52] A.I. Neishtadt. Persistence of stability loss for dynamical bifurcations. II.
Differential Equa-tions Translations , 24:171–176, 1988.[53] J.M. Ottino and D.V. Khakhar. Mixing and segregation of granular materials.
Ann. Rev.Fluid Mech. , 32:55–91, 2000.[54] T. Peacock and J. Dabiri. Introduction to focus issue: Lagrangian coherent structures.
Chaos ,10(1):017501, 2010.[55] H. Power and B. Febres de Power. Second-kind integral equation formulation for the slowmotion of a particle of arbitrary shape near a plane wall in a viscous fluid.
SIAM J. Appl.Math. , 53(1):60–70, 1993.[56] D.O. Pushkin, D.E. Melnikov, and V.M. Shevtsova. Ordering of small particles in one-dimensional coherent structures by time-periodic flows.
Phys. Rev. Lett. , 106:234501, 2011.[57] A. Rainer. Differentiable roots, eigenvalues, and eigenvectors.
Israel J. Math. , 201(1):99–122,2014.[58] M. Rasmussen.
Attractivity and Bifurcation for Nonautonomous Dynamical Systems . Springer,2007.[59] M. Rasmussen. Finite-time attractivity and bifurcation for nonautonomous differential equa-tions.
Differential Equations Dynam. Systems , 18(1):57–78, 2010.[60] F. Roman`o and H.C. Kuhlmann. Interaction of a finitesize particle with the moving lid of acavity.
Proc. Appl. Math. Mech. , 15(1):519–520, 2015.2261] F. Roman`o and H.C. Kuhlmann. Smoothed-profile method for momentum and heat trans-fer in particulate flows.
Int. J. Num. Meth. Fluids , pages 1–28, 2016. early view, DOI:10.1002/fld.4279.[62] U. Rosebrock, P.R. Oke, and G. Carroll. An application framework for the rapid deploymentof ocean models in support of emergency services: application to the MH370 search. In
Environmental Software Systems. Infrastructures, Services and Applications , pages 235–241.Springer, 2015.[63] J. Rubin, C.K.R.T. Jones, and M. Maxey. Settling and asymptotic motion of aerosol particlesin a cellular flow field.
J. Nonlinear Sci. , 5:337–358, 1995.[64] S.I. Rubinow and J.B. Keller. The transverse force on a spinning sphere moving in a viscousfluid.
J. Fluid Mech. , 11(3):447–459, 1961.[65] S. Schecter. Persistent unstable equilibria and closed orbits of a singularly perturbed equation.
J. Differential Equat. , 60:131–141, 1985.[66] D. Schwabe, A.I. Mizev, M. Udhayasankar, and S. Tanaka. Formation of dynamic particleaccumulation structures in oscillatory thermocapillary flow in liquid bridges.
Phys. Fluids ,19:072102, 2007.[67] L.E. Scriven and C.V. Sternling. The Marangoni effects.
Nature. , 187:186–188, 1960.[68] S.C. Shadden, F. Lekien, and J.E. Marsden. Definition and properties of Lagrangian coherentstructures from finite-time Lyapunov exponents in two-dimensional aperiodic flows.
Phys. D ,212(3):271–304, 2005.[69] H.A. Stone, A.D. Stroock, and A. Ajdari. Engineering flows in small devices: microfluidicstoward a lab-on-a-chip.
Annu. Rev. Fluid Mech. , 36:381–411, 2004.[70] S. Tanaka, H. Kawamura, I. Ueno, and D. Schwabe. Flow structure and dynamic particleaccumulation in thermocapillary convection in a liquid bridge.
Phys. Fluids , 18:067103, 2006.[71] R. Vasseur and R.G. Cox. The lateral migration of spherical particles sedimenting in a stagnantbounded fluid.
J. Fluid Mech. , 80(3):561–591, 1977.[72] M. Wechselberger. A propos de canards (apropos canards).
Trans. Amer. Math. Soc. , 364:3289–3309, 2012.[73] B.H. Xu and A.B. Yu. Numerical simulation of the gas-solid flow in a fluidized bed by combiningdiscrete particle method with computational fluid dynamics.
Chem. Eng. Sci. , 52(16):2785–2809, 1997.[74] J. Young and A. Leeming. A theory of particle deposition in turbulent pipe flow.
J. FluidMech. , 340:129–159, 1997.[75] Z. Zhang, C. Kleinstreuer, C.S. Kim, and Y.S. Cheng. Vaporizing microdroplet inhalation,transport, and deposition in a human upper airway model.