The Lagrangian transport in the laminar incompressible flow in a two-dimensional square cavity driven by a harmonic tangential oscillation of the lid is investigated numerically for a wide range of Reynolds and Strouhal numbers. The topology of fluid trajectories is analyzed by stroboscopic projections revealing the co-existence of chaotic trajectories and regular Kolmogorov–Arnold–Moser (KAM) tori. The pathline structure strongly depends on the Reynolds number and the oscillation frequency of the lid. Typically, most pathlines are chaotic when the oscillation frequency is small, with few KAM tori being strongly stretched along instantaneous streamlines of the flow. As the frequency is increased, the fluid motion becomes more regular and the size of the KAM tori grows until, at high frequencies, they resemble streamlines of a mean flow.

The steady lid-driven cavity flow serves as a popular computational benchmark, due to the combination of simple geometry and numerical challenges due to the discontinuous boundary conditions (Botella Peyret, 1998; Gelfgat, 2006; Hansen and Kelmanson, 1994; Gupta et al., 1981). At the same time, lid-driven cavity flows exhibit fundamental fluid-dynamical phenomena such as an infinite sequence of viscous corner eddies (Moffatt, 2019; Moffatt, 1964; Kalita et al., 2018), hydrodynamic instabilities (Albensoeder et al., 2001; Kuhlmann and Albensoeder, 2014), and finite-size Lagrangian coherent structures of advected particles (Romanò et al., 2019a; Wu et al., 2021; Wu et al., 2023). Comprehensive reviews of the lid-driven cavity problem are due to Shankar and Deshpande (2000) and Kuhlmann and Romanò (2018).

In a two-dimensional cavity flow, a harmonic modulation of the lid velocity introduces, moreover, a significant difference between separation points of streamlines and those of streaklines (Han et al., 2018), a spatiotemporal symmetry breaking (Vogel et al., 2003; Blackburn and Lopez, 2003), oscillatory boundary layers (O'Brien, 1975), and the coexistence of regular (periodic or quasi-periodic) and chaotic fluid trajectories (Ottino, 1989). The former arise on nested Kolmogorov–Arnold–Moser (KAM) tori in the space spanned by the two spatial coordinates and time. Each set of KAM tori houses a periodic pathline in its center. While quasi-periodic trajectories wind on KAM tori, chaotic trajectories fill the space between the KAM tori, called the chaotic sea.

The KAM tori are material surfaces, i.e., barriers for convective fluid transport. Thus, their structure (called Lagrangian topology) is closely related to the transport and mixing of passive scalars such as the concentration of chemical species or temperature in a forced convection flow (i.e., when the effect of the scalar on the fluid properties is neglected). The mixing in laminar flows is of practical interest, e.g., in microdevices (Petkovic et al., 2017) and several other industrial and natural flows, as summarized by Speetjens et al. (2021). In particular, the micro-fluid applications benefit from efficient mixing under low-Reynolds-number flow conditions using simple geometries, which are relatively easy to manufacture. Thus, the presence of chaotic trajectories in steady spatially periodic duct flows has been utilized for the design of continuous in-line mixers as, e.g., the rotated arc mixer (Speetjens et al., 2014). Qualitatively, the same physical phenomena that govern the mixing in such mixers can be also investigated in two-dimensional time-periodic flows. A recent practical application of an electroosmotic in-line homogenizer for the reduction of Taylor–Aris dispersion is due to Biagioni et al. (2022). Furthermore, the Lagrangian topology strongly rules the transport of small solid particles. Romanò et al. (2019b) have shown that density-matched particles can accumulate solely due to their finite size if they are advected in a flow with closed streamlines and KAM tori, which arrive sufficiently close to a tangentially moving boundary or to a moving free surface.

Although the Lagrangian topology in laminar flows and, in particular, in time-periodic two-dimensional lid-driven cavities has been thoroughly investigated in the literature, most studies (Chien et al., 1986; Ottino et al., 1988; Leong and Ottino, 1989; Franjione et al., 1989) considered a Stokes flow. The common feature of closed time-periodic Stokes flows with chaotic advection is that the net Lagrangian transport relies on a net tangential displacement of the driving wall after one complete period of the driving protocol (see also Jana et al., 1994, for a general dependence of the Lagrangian topology on the time-periodic perturbation of a steady driving velocity in Stokes flows). The emergence of chaotic trajectories due to inertia has been analyzed primarily in three-dimensional flows (Bajer and Moffatt, 1992; Castelain et al., 2001; Pouransari et al., 2010; Kang and Anderson, 2014; Contreras et al., 2017).

The two-dimensional inertial oscillatory cavity flow has been considered by Anderson et al. (2000) for a fixed mean Reynolds number Rem=UmL/ν=50, where Um, L, and ν are the mean velocity of the lid, the height of the cavity, and the kinematic viscosity of the fluid, respectively. They show that the extent of the chaotic sea increases with the amplitude of the lid-velocity oscillation. Takasaki et al. (1994) quantified mixing in the same setup for a range of Reynolds and Strouhal numbers. They found an optimal driving frequency for fast mixing 0.7Strm1, where Strm=ΩL/(2πUm) and Ω is the circular frequency. The optimal frequency was not very sensitive to the Reynolds number nor the amplitude of the velocity oscillation. The underlying Lagrangian topology was, however, not analyzed in detail. Horner et al. (2002) found a qualitatively similar dependence of the mixing quality on the driving frequency for an open shear-driven cavity, but the optimal driving frequency for mixing Str=ΩL/(2π(U+Um))0.1, where U is the amplitude of driving velocity oscillation, is one order of magnitude smaller compared to the closed lid-driven cavity. They show that there exists a central flow region consisting of KAM tori, which hinders mixing. However, the size of these KAM tori decreases when the Reynolds number Re=(U+Um)L/ν is increased.

Experimentally, a nearly two-dimensional time-periodic cavity flow can be realized in the central region between the two end walls of a three-dimensional periodically driven cavity, if it is sufficiently extended in the spanwise direction. This flow can become unstable with respect to a small three-dimensional perturbation. The linear stability boundary Rec of such a two-dimensional harmonically modulated cavity flow with zero mean, Rem=0, as a function of the Stokes number StΩL2/(2πν) has been determined experimentally by Vogel et al. (2003) and numerically by Blackburn and Lopez (2003) for a cavity aspect ratio ΓW/L=2, where W is the cavity width. In both investigations, the bifurcating flow was synchronous and was breaking the translational symmetry in the spanwise direction (z) of the two-dimensional basic flow. The critical Reynolds number reaches a minimum of Rec=532.5 at St=20. Within the investigated range St[10,60],Rec is higher than the critical Reynolds number of the steady two-dimensional lid-driven cavity flow. The latter has been determined for the same aspect ratio by Albensoeder et al. (2001) as Rem,c=353. While the linear stability boundary of the periodically driven flow in a square cavity is not known to the authors, the critical Reynolds number for a steady motion of the lid is Rem,c=786 (Albensoeder et al., 2001). Considering the qualitative observations of Vogel et al. (2003) and Blackburn and Lopez (2003), we thus expect the harmonically driven flow to remain stable to three-dimensional perturbations at least up to the largest Reynolds number Re=500 considered in our present study.

The aim of this article is to investigate the dependence of the Lagrangian topology in a two-dimensional lid-driven cavity flow with a harmonically modulated lid velocity on the Reynolds number Re=(U+Um)L/ν and the Strouhal number Str=ΩL/(2π(U+Um)). The mathematical problem is formulated in Sec. II. Section III describes the methods of analysis, and Sec. III B deals with the numerical methods employed. Results are presented in Sec. IV, where the dependence of the flow structure (Sec. IV A) and its Lagrangian topology (Secs. IV B and IV C) on the governing parameters is discussed. The observations are summarized in Sec. V.

Compared with Horner et al. (2002), we employ a different approach to isolate the effect of inertia. Namely, we consider the case of a zero mean velocity of the lid Um = 0 (Secs. IV A and IV B). For a zero mean lid velocity, any net transport of fluid after one period of driving relies on inertia. The effect of the mean lid velocity is illustrated for a selected case in Sec. IV C. The present study explores the Lagrangian topology in a wider range of Reynolds numbers Re[1,500] and driving frequencies Str[0.01,1] compared to previous works.

We consider the two-dimensional flow of an incompressible Newtonian fluid of density ρ and kinematic viscosity ν in a square cavity with equal height and width L = W (Fig. 1). The flow is induced by a lid moving tangentially with velocity U(t)=[Um+Ucos(Ωt)]ex, where U is the amplitude of the lid velocity oscillation, Ω its circular frequency, and Um its mean. Using the length, velocity, time, and pressure scales L, U+Um,2π/Ω, and ρ(U+Um)2, respectively, the problem is governed by the Navier–Stokes and continuity equations

(1)
Strut+u·u=p+1Re2u,
(1a)
·u=0.
(1b)

FIG. 1.

Sketch of a square cavity with unit width-to-length ratio W/L=1, which is driven by a tangentially moving lid with velocity consisting of a steady component Um and an oscillating component with amplitude U and circular frequency Ω. The origin of the coordinate system is located in the cavity center.

FIG. 1.

Sketch of a square cavity with unit width-to-length ratio W/L=1, which is driven by a tangentially moving lid with velocity consisting of a steady component Um and an oscillating component with amplitude U and circular frequency Ω. The origin of the coordinate system is located in the cavity center.

Close modal

The solution of Eq. (1) in the non-dimensional domain D=[0.5,0.5]×[0.5,0.5], corresponding to the physical domain sketched in Fig. 1, must satisfy the boundary conditions

u(x=±1/2)=0,
(2a)
u(y=1/2)=0,
(2b)
u(y=1/2)=S+cos(2πt)S+1ex,
(2c)

where S=Um/U is the ratio of the mean lid velocity to the amplitude of the lid velocity oscillation. The Reynolds and Strouhal numbers are defined using the maximum velocity

Re=(U+Um)Lν,Str=ΩL2π(U+Um).
(3)

The Strouhal number Str is used here to measure the driving frequency with respect to the convective timescale. Alternatively, the Stokes number St=ReStr=L2Ω/(2πν), which quantifies the frequency in viscous scaling, can be employed.

For small Reynolds numbers, the flow evolves, after an initial transient, to a time-periodic flow u(x,t)=u(x,t+k),k. The advection of an arbitrary fluid element with position vector x=(x,y)T in an incompressible two-dimensional (time-periodic) flow is described by the dynamical system

dxdt=u(x,t)=(ψ/yψ/x),
(4)

where ψ(x,y,t) is the stream function. For any initial condition X=x(t0)D, there exists a unique trajectory x(t)=Φt(X), where Φt is the formal solution x(t) of Eq. (4). It is well known (Aref et al., 2017) that in the case of a two-dimensional time-dependent incompressible flow, Eq. (4) can be interpreted as a Hamiltonian dynamical system with 1.5 degrees of freedom, where the time-periodic stream function ψ(x,y,t)=ψ(x,y,t+k) takes the role of the Hamiltonian. It allows the coexistence of chaotic and regular trajectories, the latter in form of Kolmogorov–Arnold–Moser (KAM) tori.

We imagine time as the third dimension. In the extended space (x, y, t), a closed streamline in a steady flow defines a straight tube parallel to the t axis. A fluid element then spirals on the tube defined by the closed streamline on which it moves in the two-dimensional physical space (x, y). When the two-dimensional flow becomes time-periodic, only subsets of the straight tubes survive and become wavy in the t direction with period nk, n, where k is the period of the driving. Therefore, only part of the fluid elements spiral on periodic tubes (see Fig. 8), while others move chaotically. A single period of such a tube is also called a period-n KAM torus. The KAM tori are densely nested, and the contour of a nested set of tori at a given instant of time is denoted a regular island [see Fig. 7(h)]. Non-periodic trajectories in the extended space (x, y, t) fill the space between the KAM tubes and belong to the chaotic sea.

For a discussion of the Lagrangian transport, we consider a strictly time-periodic flow with a non-dimensional period of 1. A complete picture of the transport, further called Lagrangian topology, can be obtained by taking the stroboscopic projection (or Poincaré section)

P(X)=k=0KΦk(X),k,
(5)

which is the collection of subsequent positions Φk(X) at discrete time levels t=t0+k of an initial position X at t = t0. The stroboscopic projection of a set of initial conditions distributed over the domain reveals the structure of fluid trajectories in terms of cross sections of KAM tori and the chaotic sea in the (x, y) Poincaré plane defined by t0. While the projection P depends on the initial time or the phase of the strobing ϕ0=2πt0, its choice only affects the positions and shapes of the cross sections of the KAM tori. Their cross-sectional areas are independent of ϕ0. Here, we fix the phase by defining ϕ0=0 as the phase at which the lid displacement is zero and its velocity has a positive maximum, i.e., u(y=1/2,ϕ0)=ex. At this phase, as well as for ϕ0=π, the KAM tori appear relatively compact and are not excessively stretched.

A hallmark of chaos is the existence of transverse intersections between stable and unstable manifolds of hyperbolic points, called transverse homoclinic/heteroclinic connections or homoclinic/heteroclinic points (Ottino, 1989). In order to reconstruct the stable and unstable manifolds (Fig. 13), we compute, respectively, the advection forward and backward in time of a small circle initialized around a hyperbolic point detected by the stroboscopic projection. The circle is numerically approximated by a set of connected points, which are advected according to Eq. (4). In the course of the temporal evolution, the circle is stretched into two graphically indistinguishable curves housing the unstable manifold, which can thus be graphically identified by the apparent (double) line into which the circle transforms. The (double) line, which emerges from the stretched circle, is locally reconstructed by adaptive cubic spline interpolation such that the maximum distance between neighboring points is bounded by Δl/(1+ακ), where Δl=5×103,α=15Δl/π and κ is the local curvature of the interpolation. The refinement algorithm is described in more detail by Meunier and Villermaux (2010).

1. Flow field

The velocity field u(x,y,t) is computed by solving (1) using the spectral-element code NEK5000 (Fischer et al., 2008). The Galerkin weak formulation of the governing equations is discretized on 20 × 20 square elements of a uniform size, where the variation of the velocity components across each element in both spatial directions is approximated by Lagrange polynomials of the order N = 7 on Gauss–Legendre–Lobatto nodes. The pressure is approximated by Lagrange polynomials of the order N−2 on Gauss–Legendre nodes (PNPN2 method). The non-linear convective term is over-integrated in the weak formulation with a 12-point (3(N+1)/2) Gauss–Legendre–Lobatto quadrature in each direction to avoid aliasing (Mengaldo et al., 2015). For the temporal discretization, a third-order backward-differencing scheme is used for the linear terms, and a third-order Adams–Bashforth extrapolation is applied to the convective terms, yielding an overall third-order accurate scheme. A constant time step of size Δt=103 was used for all computations, leading to a Courant number C0.3.

The solver was verified and the grid convergence was tested by computing a steady two-dimensional lid-driven cavity flow (S1=0) starting from rest with u(t=0)=0 and u(y=1/2)=ex. Steady state was signaled by the termination criterion

maxi,j|ui(xj,t)ui(xj,tΔt)|Δt107,
(6)

where j enumerates all nodal points. The solution was compared to the benchmark spectral data of Botella and Peyret (1998). As shown in Table I, the values of the velocity extrema along the centerlines of the cavity agree with the benchmark data up to O(106) for Re=100.

TABLE I.

Extrema [umin(x=0),vmax(y=0),vmin(y=0)] of the velocity components and their respective locations (ymin,xmax,xmin) along the centerlines for the steady flow in a lid-driven square cavity at Re=100 computed without (NR) and with (R) regularization of the lid velocity (see the text) according to Eq. (8). The reference data of Botella and Peyret (1998) are indicated by (B and P).

ReferenceGriduminyminvmaxxmaxvminxmin
NR 101 × 101 0.214 043 4 −0.0419 0.179 574 5 −0.2630 0.253 808 9 0.3104 
NR 141 × 141 0.214 042 3 −0.0419 0.179 572 7 −0.2630 0.253 803 1 0.3104 
141 × 141 0.214 041 3 −0.0419 0.179 568 1 −0.2630 0.253 800 2 0.3104 
B and P 96 × 97 0.214 042 4 −0.0419 0.179 572 8 −0.2630 0.253 803 0 0.3104 
ReferenceGriduminyminvmaxxmaxvminxmin
NR 101 × 101 0.214 043 4 −0.0419 0.179 574 5 −0.2630 0.253 808 9 0.3104 
NR 141 × 141 0.214 042 3 −0.0419 0.179 572 7 −0.2630 0.253 803 1 0.3104 
141 × 141 0.214 041 3 −0.0419 0.179 568 1 −0.2630 0.253 800 2 0.3104 
B and P 96 × 97 0.214 042 4 −0.0419 0.179 572 8 −0.2630 0.253 803 0 0.3104 

All calculations for periodic forcing were likewise initiated from a state of rest. The convergence of the flow field to a time-periodic state was checked by monitoring the magnitude of the maximum local increment of any velocity component after one period of forcing

Rk=maxi,j|ui(xj,k)ui(xj,k1)|,k.
(7)

Figure 2 shows the evolution of Rk as a function of the number of forcing periods k for S = 0. The number of periods after which a time-periodic state is reached, signaled by Rk<ϵ=107, increases with Re and Str. For Re=500, we find convergence after k = 1 for Str=0.02 and after k = 100 for Str=1. Figure 3 confirms that the converged flow is indeed periodic with the same period as the lid oscillation. Once Rk has dropped below the tolerance level of ϵ=107, the computation was terminated. It was then initiated again with initial conditions corresponding to the last time step, and 1000 snapshots of the velocity field were stored over one period of the forcing. The further temporal evolution of the flow was then constructed as a periodic continuation of this single period. A subharmonic response of the flow was not found.

FIG. 2.

Maximum local increment Rk of any velocity component after one period of forcing as a function of the number of forcing periods k passed since the lid oscillation was initiated.

FIG. 2.

Maximum local increment Rk of any velocity component after one period of forcing as a function of the number of forcing periods k passed since the lid oscillation was initiated.

Close modal
FIG. 3.

Fourier amplitudes û (black) and v̂ (blue) of the velocity components u and v, respectively, as functions of the dimensionless frequency f at the point (x,y)=(1/4,1/4) obtained from t[20,120] for Re=500 and Str=0.1.

FIG. 3.

Fourier amplitudes û (black) and v̂ (blue) of the velocity components u and v, respectively, as functions of the dimensionless frequency f at the point (x,y)=(1/4,1/4) obtained from t[20,120] for Re=500 and Str=0.1.

Close modal

For computing the flow under time-periodic driving with a non-zero mean (Sec. IV C), the velocity boundary condition at y = 1/2 (moving lid) was replaced by

u(x,y=1/2,t)=[12 cosh(2x/ξ)exp(1/ξ)]S+cos(2πt)S+1ex.
(8)

For the value of ξ=0.01 used, the residual velocity jump at x=±0.5 is insignificant. Table I shows that in the case of a steady driving (S1=0), this regularization has no significant effect on the flow in the bulk of the cavity.

2. Streamlines

Streamlines of the instantaneous flow field were obtained from isolines of the stream function ψ(x,y,t), which was computed by solving the Poisson equation

2ψ=ω
(9)

in Galerkin weak formulation and subject to Dirichlet boundary conditions ψ = 0 on all boundaries. For the solution of Eq. (9), the same grid of spectral elements and the same polynomial order were used as for the flow solver. The vorticity ω=v/xu/y was computed by the subroutine comp_vort3 of NEK5000. We also consider the mean flow

u¯(x)=01u(x,t)dt11000n=11000un(x),
(10)

where un(x)=u(x,tn) and n enumerates 1000 uniformly distributed time steps returned by the flow solver per one period. The mean flow becomes relevant for the Lagrangian transport at high forcing frequencies, as discussed in Sec. IV.

3. Trajectories of fluid elements

The trajectories of fluid elements were computed by solving (4) using the Dormand–Prince Runge–Kutta method (Dormand and Prince, 1980; Shampine and Reichelt, 1997) as implemented in the MATLAB function ode45. It is an explicit adaptive time-stepping algorithm that compares the results of fourth- and fifth-order Runge–Kutta formulas to estimate the error at each time step. An absolute and relative tolerance of 5×107 was used in the present study. The computation of the trajectories requires an interpolation of the velocity field in space and time, which is described in  Appendix A.

4. Stroboscopic projection

To numerically approximate the Lagrangian topology, 10 × 10 up to 25 × 25 fluid elements were initialized on a regular uniform grid of X. The marked fluid elements were then tracked over 100 up to 1000 periods of the forcing, depending on the sizes of the KAM tori to be visualized and the duration τ (further called winding period) required for a quasi-periodic trajectory to make one revolution about the central closed pathline. For some tori in the region of slow motion near the bottom wall (y=1/2), the winding period may be very large with τ1000, exceeding the time K employed for the stroboscopic projection (5). Therefore, the contours in the Poincaré plane of KAM tori near the bottom of the cavity might not be fully covered [see Figs. 7(a) and 7(c)]. Boundaries between the chaotic sea and regular islands (shown in Fig. 8) were obtained as α-shapes of trajectories belonging to the largest reconstructible torus of a set of nested KAM tori.

In order to test the sensitivity of the results to the number of spectral elements employed by the flow solver, we computed the stroboscopic projections for Re=50,Str=0.25, and S = 0 with two different meshes made of 10 × 10 and 20 × 20 spectral elements. The Lagrangian topologies revealed by the stroboscopic projections are shown in Fig. 4. The results are almost indistinguishable from each other. As indicated above, the finer mesh of 20 × 20 elements is employed throughout the present study.

FIG. 4.

Stroboscopic projections for Re=50,Str=0.25, and S = 0 computed using 10 × 10 (a) and 20 × 20 (b) spectral elements.

FIG. 4.

Stroboscopic projections for Re=50,Str=0.25, and S = 0 computed using 10 × 10 (a) and 20 × 20 (b) spectral elements.

Close modal

In the following, we shall primarily consider the case of harmonic modulation of the lid velocity with zero mean (S = 0) such that Re=UL/v. Only in Sec. IV C, the influence of a non-zero mean lid velocity is studied.

In the Stokes flow limit Re0, the inertial term u·u in Eq. (1) can be neglected and the Eulerian velocity field shares the temporal symmetries of the forcing function such that u(x,y,t)=ũ(x,y)cos(2πt). Under these conditions, the Lagrangian topology is trivial

Φt=1(X)=X,XD,
(11)

and every fluid element returns to its initial position after one period of forcing (Taylor, 1966). Furthermore, the streamlines exhibit a reflection symmetry with respect to x = 0 at any time, i.e.,

(ũṽ)(x,y)=(ũṽ)(x,y).
(12)

For inertial flow and S = 0, the symmetries (11) and (12) are reduced to a reflection–translation symmetry (RT symmetry, Vogel et al., 2003), i.e., reflections about the vertical centerline together with a temporal translation by half a period,

(uv)(x,y,t)=(uv)(x,y,t+1/2).
(13)

For vanishing lid frequency Str0, the partial derivative term Stru/t in Eq. (1a) becomes negligibly small and the flow is quasi-steady. It thus exhibits the known characteristics of classical cavity flow under stationary driving. For unit aspect ratio, these are the primary vortex occupying most of the cavity and an infinite sequence of self-similar vortices (Moffatt, 1964) near the corners, where two stationary walls meet.

As the lid acceleration increases, the velocity fields deviate from the quasi-steady ones and the complexity of the instantaneous flow increases with Str and Re (Zhu et al., 2020). Since the primary vortex created by the uni-directional lid motion during one half-period survives longer than a half-period, this vortex will coexist with the newly created counter-rotating primary vortex initiated near the lid when it changes the direction of motion. Figure 5 shows the evolution during the acceleration phase of the lid in the positive direction: A shear layer is created near the moving lid. This layer rolls up into a vortex near the downstream end of the cavity near (x,y)=(0.5,0.5) (a). The vortex grows and replaces the previous vortex (b), which eventually detaches from the wall near the upstream corner (x,y)=(0.5,0.5) of the cavity (c) and finally vanishes. The horizontal velocity profile along a vertical centerline (blue in Fig. 5) resembles the profile arising in Stokes' second problem (also pointed out by O'Brien, 1975). As the relative thickness of the Stokes layer λ=St1/2 decreases with increasing Stokes number, the velocity profile at x = 0 in the square cavity approaches the profile for the Stokes problem. The thickness of the Stokes layer is relevant for the mixing capability of the flow, as will be discussed in Sec. IV B 2.

FIG. 5.

Streamlines (black) and velocity profile u(x=0) on the vertical centerline (blue) for Str=1 and Re=St=500 shown during the acceleration phase of the lid in the positive x direction. Full and dashed lines indicate negative (clockwise rotation) and positive stream functions (counterclockwise rotation), respectively. The phase ϕ=2π(tmod1) is indicated for each subfigure. (a) ϕ = 1.625π, (b) ϕ = 1.75π, and (c) ϕ = 1.9π.

FIG. 5.

Streamlines (black) and velocity profile u(x=0) on the vertical centerline (blue) for Str=1 and Re=St=500 shown during the acceleration phase of the lid in the positive x direction. Full and dashed lines indicate negative (clockwise rotation) and positive stream functions (counterclockwise rotation), respectively. The phase ϕ=2π(tmod1) is indicated for each subfigure. (a) ϕ = 1.625π, (b) ϕ = 1.75π, and (c) ϕ = 1.9π.

Close modal

Even for a zero mean velocity of the lid S = 0, there exists a non-zero mean flow u¯ (10) for Re0. The RT symmetry implies that the mean flow is reflection symmetric. It is shown in Fig. 6 for different parameters. The mean flow exhibits a symmetric four-vortex structure with the two counter-rotating upper vortices being dominant in strength and size. With the increasing frequency of the lid (Str), above a certain threshold, the mean flow vortices in the upper half of the cavity grow in size at the expense of the mean flow vortices in the lower part of the cavity and all vortex centers tend to move toward the corners of the cavity.

FIG. 6.

Magnitude of the mean velocity u¯ (color) and the corresponding streamlines (black) for S = 0; Re and Str as indicated. (a) Str = 0.02, Re = 1, (b) Str = 0.02, Re = 200, (c) Str = 0.02, Re = 500, (d) Str = 1, Re = 1, (e) Str = 1, Re = 200, and (f) Str = 1, Re = 500.

FIG. 6.

Magnitude of the mean velocity u¯ (color) and the corresponding streamlines (black) for S = 0; Re and Str as indicated. (a) Str = 0.02, Re = 1, (b) Str = 0.02, Re = 200, (c) Str = 0.02, Re = 500, (d) Str = 1, Re = 1, (e) Str = 1, Re = 200, and (f) Str = 1, Re = 500.

Close modal

The Lagrangian topology has a well-known generic structure for two-dimensional time-periodic flows, consisting of KAM tori as well as chaotic trajectories. This section concerns the variation of the topology as a function of the parameters Re and Str for S = 0.

1. KAM tori

Figure 7 gives an overview of the Lagrangian topology in terms of stroboscopic projections onto the plane ϕ0=0 for S = 0 with Str and Re covering the range investigated. The continuation of the KAM tori in the three-dimensional phase space is shown in Fig. 8 for representative cases.

FIG. 7.

Lagrangian topologies for S = 0 and different combinations of Re and Str as indicated. The stroboscopic projection is carried out for K1000 and ϕ0=0. Black and magenta dots indicate regular trajectories, while gray dots indicate chaotic trajectories. Blue lines in (g) are streamlines of the velocity field at ϕ=0. For regular trajectories near the bottom of the cavity y=0.5, only every 20th to 100th intersection with the stroboscopic plane is plotted. (a) Str = 1, Re = 1, (b) Str = 1, Re = 10, (c) Str = 1, Re = 500, (d) Str = 0.05, Re = 1, (e) Str = 0.05, Re = 10, (f) Str = 0.05, Re = 500, (g) Str = 0.01, Re = 1, (h) Str = 0.01, Re = 10, and (i) Str = 0.01, Re = 500. (a) and (g) are reproduced with permission from L. Babor and H. C. Kuhlmann, Proc. Appl. Math. Mech. 20, e202000194 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution Licence.

FIG. 7.

Lagrangian topologies for S = 0 and different combinations of Re and Str as indicated. The stroboscopic projection is carried out for K1000 and ϕ0=0. Black and magenta dots indicate regular trajectories, while gray dots indicate chaotic trajectories. Blue lines in (g) are streamlines of the velocity field at ϕ=0. For regular trajectories near the bottom of the cavity y=0.5, only every 20th to 100th intersection with the stroboscopic plane is plotted. (a) Str = 1, Re = 1, (b) Str = 1, Re = 10, (c) Str = 1, Re = 500, (d) Str = 0.05, Re = 1, (e) Str = 0.05, Re = 10, (f) Str = 0.05, Re = 500, (g) Str = 0.01, Re = 1, (h) Str = 0.01, Re = 10, and (i) Str = 0.01, Re = 500. (a) and (g) are reproduced with permission from L. Babor and H. C. Kuhlmann, Proc. Appl. Math. Mech. 20, e202000194 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution Licence.

Close modal
FIG. 8.

Stroboscopic projections (dots) and outermost KAM tori (surfaces) of selected sets for S = 0 shown over one full period of the driving. The low-opacity surfaces in (d) represent isosurfaces of the instantaneous stream function for ψ=±0.4 and ±0.6. The phase at the front face of the domain corresponds to ϕ0=0 in (a)–(c) and to ϕ0=3π/2 in (d). (a) Str = 0.1, Re = 20, (b) Str = 0.25, Re = 50, (c) Str = 0.01, Re = 5, and (d) Str = 0.05, Re = 500. (d) is reproduced with permission from L. Babor and H. C. Kuhlmann, Proc. Appl. Math. Mech. 20, e202000194 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution Licence.

FIG. 8.

Stroboscopic projections (dots) and outermost KAM tori (surfaces) of selected sets for S = 0 shown over one full period of the driving. The low-opacity surfaces in (d) represent isosurfaces of the instantaneous stream function for ψ=±0.4 and ±0.6. The phase at the front face of the domain corresponds to ϕ0=0 in (a)–(c) and to ϕ0=3π/2 in (d). (a) Str = 0.1, Re = 20, (b) Str = 0.25, Re = 50, (c) Str = 0.01, Re = 5, and (d) Str = 0.05, Re = 500. (d) is reproduced with permission from L. Babor and H. C. Kuhlmann, Proc. Appl. Math. Mech. 20, e202000194 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution Licence.

Close modal

For the phase ϕ0=0, the KAM tori are arranged roughly reflection symmetrically with respect to the vertical midplane x = 0. This symmetry is not exact due to the flow inertia associated with Re>0. The breaking of the symmetry is more pronounced for higher values of Re and longer oscillation periods Str1. Due to the exact RT symmetry (13) of the flow when the lid velocity has a zero mean, the KAM tori must arise in pairs that satisfy this symmetry. This is visually confirmed in Fig. 8.

For the small Reynolds number Re=1 and Str[0.01,1] [first column of Fig. 7: Fig. 7(a), 7(d), and 7(g)], only regular trajectories are found. As Re increases the sea of chaotic advection emerges visibly and grows, with regular subharmonic islands appearing occasionally [Figs. 7(c), 7(e), 7(f), 7(h), and 7(i)].

At the highest frequency investigated [Str=1, Fig. 7(a)], three pairs of regular regions are observed when Re=1. The KAM tori resemble the streamlines of the mean flow [Fig. 6(c)] for the same parameters, except for the upper pair of tori, which are stretched along the moving lid. The symmetry of the KAM tori with respect to the midplane is only weakly broken. As the Reynolds number increases [Figs. 7(b) and 7(c)], the periodic trajectories in the centers of the two largest sets of tori are displaced toward the singular corners (x,y)=(±0.5,0.5), just like the streamlines of the mean flow (Fig. 6), while the smaller pairs of KAM tori with period 1 near the bottom corners and the moving lid shrink. Differences between KAM tori and the mean streamlines remain, however, in regions where the mean flow is too weak compared to the oscillatory part of the flow. As a result, chaotic and subharmonic KAM islands can be created in these regions, as shown by the inset in Fig. 7(c).

An interesting behavior is found for small Reynolds numbers (Re=1) when the Strouhal number is decreased [Figs. 7(a), 7(d), and 7(g)]. For Str=0.05, the regular regions next to the lid, present for Str=1, have vanished. Furthermore, for ϕ0, a free stagnation point that splits into two hyperbolic points is created near the center of the steady primary vortex when the Strouhal number drops below Str5.9×102. For steady Stokes flow, this center is located at (x,y)=(0,0.2600) (Shankar, 1993). As a result, a new pair of KAM tori is born. With a further decrease in Str, the newly created KAM tori grow in size [Fig. 7(d)] until another pair of tori is generated at the same location. This process of creation of separated regular regions repeats itself as Str becomes smaller. Those KAM tori that have come into existence at higher Strouhal numbers are displaced toward the walls and get stretched along the instantaneous streamlines, which are close to Stokesian for Re=1 and which do not change very much during the whole period. Already for Str=0.01 and Re=1, each separated regular region is almost perfectly bounded by two instantaneous streamlines [blue lines in Fig. 7(g)], except for the pair of KAM tori near bottom corners. The boundaries of these latter tori deviate from the instantaneous streamlines due to the large difference between Eulerian and Lagrangian separation points. For small Str and Re, the layered KAM tori become very thin and, for Str1, wind many times about the center of the quasi-steady vortex during half of the period. In the following half period, the tori wind in the opposite direction to reconstruct the initial state [see Fig. 8(c) for Re=5].

This behavior can be explained as follows. In the case of a Stokes flow with Re0, every point on a streamline is a fixed point of the Poincaré map, because of the time reversibility of Stokes flow. Furthermore, each fluid element moves on a straight KAM tube defined by a streamline of the Stokes flow. If the Stokes flow is slightly perturbed by a small non-zero Reynolds number Re=ϵ1, the points on the streamline are no longer fixed points, and the return to the Poincaré plane is slightly offset by η=(ηx,ηy). Furthermore, if the Strouhal number is very small, a fluid element initiated on a streamline of the main Stokes flow vortex makes many revolution about the vortex center, before its orbiting motion is reversed for the next half period of the forcing. In this situation, a weak nonlinear effect due to ϵ0 is expected to be mainly due to a runtime lag, which would manifest itself in an offset η of the Poincaré return, which is mainly tangential to the unperturbed streamline on which the fluid element started at t0=0.

A trajectory X0(t;x0=0,t0=0) initiated on the symmetry plane x=0 at the instant at which the harmonic lid velocity has an extremum (t0=0) exhibits the time inversion symmetry (X0,Y0)(t)=(X0,Y0)(t). This symmetry is broken for trajectories initiated at x00. If X+(t;x+,t0=0)=(X+,Y+)(t) is a trajectory initiated at x+>0 and t0=0, then the trajectory initiated mirror symmetrically at (x,y0)=(x+,y0) at t0=0 satisfies (X,Y)(t)=(X+,Y+)(t). Therefore, the runtime lag for X+(t) and X(t) should be antisymmetric with respect to x = 0 yielding an opposite tangential offset. This is confirmed by the numerical calculations, which yield (ηx,ηy)(ηx,ηy)+, mainly tangentially. Therefore, successive returns to the Poincaré plane t = 0 of trajectories initiated at x+ and x at t = 0 with (x,y0)=(x+,y0) approach each other in a mirror-symmetric fashion and approximately along the streamlines of the unperturbed flow (Stokes flow).

If the offset of the returns perpendicular to the unperturbed streamline vanishes, two hyperbolic fixed points are created in the Poincarè plane t = 0 at the intersections of this streamline with the symmetry plane x = 0. On the other hand, if the offset normal to the unperturbed streamline is nonzero, the sequence of the Poincaré points will rapidly turn just before reaching x = 0 and continue in opposite direction along a neighboring streamline of the unperturbed Stokes flow. This is a consequence of the stroboscopic projection of a 2D incompressible flow being an area-preserving map (Ottino, 1989). The Poincaré points then depart from each other along their new streamline, only to turn again near the other hyperbolic point on x = 0, and are transported back to the original streamline, closing the contours of two RT-symmetric KAM tori. Apparently, the offset normal to the unperturbed closed streamlines is smaller for longer periods of the driving. Thus, its sign must change more rapidly as function of the distance from the vortex center due to the area preservation of the mapping. Therefore, more hyperbolic points are created for smaller Str, which leads to the characteristic onion-like KAM tubes in the Poincaré plane t = 0 [see Fig. 7(g)] whose thickness becomes thinner as Str0. The slightly asymmetric locations of the hyperbolic points around x0 suggest heteroclinic tangling such that the KAM tori are separated by thin chaotic layers, which, however, could not be resolved numerically.

The organization of the periodic points with respect to the symmetries of the streamlines of creeping flow and to the temporal symmetries of the driving has also been noted by Leong and Ottino (1989) and was further explored by Franjione et al. (1989). The phenomenon is generic for Stokesian time-periodic flows with symmetries, also for non-zero mean velocities of the driving.

Figure 8 shows examples of the motion of the KAM tori. For Str=0.1 and Re=20 [Fig. 8(a)], the tori near the bottom corners occupy the whole Lagrangian separation zones. Since the outermost tori of these nested sets are attached to the walls, they are not reconstructible due to an infinite winding time of the limiting trajectories. Therefore, slightly smaller tori were used for visualization. For Re50 [Figs. 8(a)–8(c)], the tori of period 1 have an oscillatory (sinuous) helical shape and become stretched, while they approach the moving lid. As mentioned above, when the period Str11 is large [Fig. 8(c)], the tori that are not attached to stationary walls (red, yellow, and cyan) wind about the center of the main Eulerian vortex for multiple revolutions. The torus in gray near the left bottom corner is stretched along the left sidewall for the first half-period (part of the torus is hidden behind itself and the other tori shown), before returning to its original position near the bottom wall during the second half-period.

For Str=0.05 and large Reynolds number Re=500, the cross sections of the KAM tori at ϕ=0 are shown in Fig. 7(f), while the three-dimensional continuation of the largest reconstructible tori is represented by solid colors in Fig. 8(d). In addition, the Eulerian vortex is visualized by transparent isosurfaces of ψ. For this relatively large Reynolds number, the Eulerian vortex occupies mostly the downstream half of the cavity, depending on the instantaneous direction of the lid motion (alternating transparent blue and transparent brown isosurfaces). The two upper sets of KAM tori (blue and green) alternatingly spiral close to the center of the respective Eulerian vortex for half a period. As their respective Eulerian vortex vanishes during the second half-period, the tori are advected back to the position that they had one period earlier.

To better demonstrate the evolution of the KAM tori upon a variation of Re, detailed structures near the main KAM torus in the Poincaré plane ϕ0=0 are shown for Str=0.1 [Figs. 9(a)–9(c)] and Str=0.25 [Figs. 9(d)–9(f)]. Furthermore, the mean winding periods of quasi-periodic trajectories to complete a full revolution about their KAM torus were computed in the following way. Let us define a cumulative polar position angle θ(t) of a regular trajectory with respect to the elliptic trajectory at the center of the KAM torus. The number of windings within a given time interval t[0,K] is then

Nw(K)=θ(K)θ(0)2π.
(14)
FIG. 9.

Stroboscopic projections to ϕ0=0 showing nested KAM tori around a central periodic trajectory x(0) (blue dot). Subharmonic resonances are shown in magenta. The Reynolds and Strouhal numbers are indicated in the sub-captions and color-coded corresponding to the line color in Fig. 10. (a) Re = 5, Str = 0.1, (b) Re = 10, Str = 0.1, (c) Re = 20, Str = 0.1, (d) Re = 10, Str = 0.25, (e) Re = 20, Str = 0.25, and (f) Re = 29, Str = 0.25.

FIG. 9.

Stroboscopic projections to ϕ0=0 showing nested KAM tori around a central periodic trajectory x(0) (blue dot). Subharmonic resonances are shown in magenta. The Reynolds and Strouhal numbers are indicated in the sub-captions and color-coded corresponding to the line color in Fig. 10. (a) Re = 5, Str = 0.1, (b) Re = 10, Str = 0.1, (c) Re = 20, Str = 0.1, (d) Re = 10, Str = 0.25, (e) Re = 20, Str = 0.25, and (f) Re = 29, Str = 0.25.

Close modal

The number of windings per one period of driving is nw=Nw(K)/K. The inverse of nw gives the number of driving periods per winding

τ=1nw=2πKθ(K)θ(0).
(15)

The mean winding period τ converges with the number of driving periods K employed for the averaging. K should be sufficiently large, such that the regular trajectory completes several windings within the time interval considered. The value K = 500 was typically sufficient, except from KAM tori too close to stationary walls, for which the winding period diverges. Figure 10 displays the winding periods τ(r¯) on the main period 1 KAM tori for the same parameters as in Fig. 9 (parameters coded by color). The nested tori are parameterized by their mean distance

r¯=1500Km=1500K||x(tm)x(0)(tm)||2
(16)

from the central periodic trajectory x(0)(t), using 500 equidistant time steps tm=m/500 per period of the lid motion. From Fig. 10, the winding period τ of regular trajectories about the central periodic trajectory generally grows with r¯ but decreases with increasing Reynolds number. As r¯ gets larger, the computed values of τ develop wiggles because the trajectory approaches hyperbolic points or the chaotic sea.

FIG. 10.

Dependence of the winding time τ of a regular trajectory about the central closed streamline x(0)(t) as a function of its mean distance r¯ from x(0)(t). The parameters are given in the legend. The cases shown are distinguished by line type and color corresponding to the color of the sub-captions of Figs. 9(a)–9(e).

FIG. 10.

Dependence of the winding time τ of a regular trajectory about the central closed streamline x(0)(t) as a function of its mean distance r¯ from x(0)(t). The parameters are given in the legend. The cases shown are distinguished by line type and color corresponding to the color of the sub-captions of Figs. 9(a)–9(e).

Close modal

For Str=0.1,Re=5 [Fig. 9(a)], only synchronous (period-1) tori are observed and their winding periods are longer than 20. When Re is increased to 10, the winding is more rapid and resonant tori with τ = 11 and τ = 12 break into subharmonic tori [Fig. 9(b)]. These disappear in the chaotic sea when the Reynolds number is further increased, while synchronous tori of smaller and smaller winding period τ start to break [Fig. 9(c)].

For a slightly higher frequency with Str=0.25, the breaking of KAM tori is postponed to 10<Re<20 when several chains of subharmonic tori are created, with periodicities ranging from 5 to 14 [Fig. 9(e)]. Upon a further increase in the Reynolds number, the tori of high periodicity disappear, while some tori of lower periodicity (5) appreciably grow in size [Fig. 9(f)].

According to the Kolmogorov–Arnold–Moser theorem (Kolmogorov, 1954a; Kolmogorov, 1954b; Arnold, 1963; Möser, 1962; Rüssmann, 1970; Ottino, 1989), the tori consisting of regular trajectories with irrational values of τ survive a perturbation of the Hamiltonian system, while those with rational τ=p/q(p,q) break into subharmonic tori (Birkhoff, 1934; Helleman, 1980). These periodically return to their original locations in the Poincaré plane after p periods of the lid motion, winding q times about the central periodic trajectory. Such subharmonic tori must intersect the Poincaré plane at p distinct locations. Comparing the stroboscopic projections shown in Figs. 9(a)–9(f) with the corresponding winding times in Fig. 10, we observe that the first breaking of tori occurs when τ is an integer number [Figs. 9(b), 9(e), and 9(f)]. Afterward, also some tori with other rational values of τ start to break. Figure 9(c) shows the resonances τ=21/4 and 16/3 producing tori of periodicities 21 and 16, respectively.

2. Emergence of the chaotic sea

For low driving frequencies Str0.1 and 1<Re<10, the sea of chaotic trajectories comes into existence along heteroclinic connections between hyperbolic points of period one [synchronous trajectories, Figs. 7(a), 7(d), and 7(g)] if the Reynolds number is increased. As shown in Sec. IV B 2, KAM tori near the edges of the regular islands with resonant winding times τ=p/q break into chains of elliptic and hyperbolic periodic trajectories. These chains subsequently turn into chaotic layers containing subharmonic islands. On the other hand, for higher driving frequencies Str0.25 and Re10, the chaotic sea is initiated in the bulk of the regular regions [Fig. 7(c)] and not only near the synchronous heteroclinic connections.

Usually, the extent of the chaotic sea increases with Re, implying shrinkage and disappearance of regular regions. Only for one of the frequencies investigated, for Str=0.05, we observe the creation of new KAM tori within the chaotic sea upon an increase in Re above Re340. The evolution with Re is shown in Fig. 11.

FIG. 11.

Stroboscopic projections for Str=0.05 show the creation of regular islands within the chaotic sea as the Reynolds number is increased. The Reynolds number is specified in the sub-captions. Colors as in Fig. 7. (a) Re = 340 and (b) Re = 345.

FIG. 11.

Stroboscopic projections for Str=0.05 show the creation of regular islands within the chaotic sea as the Reynolds number is increased. The Reynolds number is specified in the sub-captions. Colors as in Fig. 7. (a) Re = 340 and (b) Re = 345.

Close modal

Less deformed islands persist up to higher Reynolds numbers than very stretched ones do. Typically, KAM tori are more stretched at lower driving frequencies and the sensitivity of the topology with respect to a variation of the Reynolds number is higher. This can be seen from Fig. 7 for Re=10: For Str=0.01 [Fig. 7(h)], the cavity is mainly occupied by chaotic trajectories, whereas for Str=1 [Fig. 7(b)], the topology is still regular. The character of the flow as a function of Re and Str is classified in Fig. 12: For low Reynolds numbers, all trajectories detected are regular (blue squares). The onset Reynolds number for chaotic streamlines seems to slightly increase with Str. For most parameters (Re,Str) investigated, we find regular and chaotic trajectories to coexist (green circles). For Re100 and Str0.03, regular trajectories are found only near the bottom corners (x,y)=(±0.5,0.5) occupying less than 1% of the domain (red diamonds).

FIG. 12.

(a) Regimes of the Lagrangian topology in which only regular trajectories exist (blue squares), regular and chaotic trajectories coexist (green circles), and in which regular regions occupy less than 1% of the domain (red diamonds). (b) Dependence of the convective mixing time τm on the Strouhal number for different Reynolds numbers Re=10 (), Re=20 (), Re=50 (), Re=100 (), Re=200 (), and Re=500 ().

FIG. 12.

(a) Regimes of the Lagrangian topology in which only regular trajectories exist (blue squares), regular and chaotic trajectories coexist (green circles), and in which regular regions occupy less than 1% of the domain (red diamonds). (b) Dependence of the convective mixing time τm on the Strouhal number for different Reynolds numbers Re=10 (), Re=20 (), Re=50 (), Re=100 (), Re=200 (), and Re=500 ().

Close modal

The mixing time in the convective scale τm, as defined in  Appendix B, is shown in Fig. 12(b). Within the range of Str considered and for Reynolds numbers Re=10 and 20, the mixing is the fastest (smallest τm) for Str=0.25. The minimum of τm moves to smaller driving frequencies when the Reynolds number increases. While for small Str the mixing is improved by a higher Reynolds number, this is not the case for the high frequency Str=1. When the Stokes layer thickness λ=(ReStr)1/2 is too small, the flow is restricted to the thin Stokes layer near the lid, while the flow in the rest of the cavity becomes too weak for efficient mixing. From a comparison of Fig. 12(a) with Fig. 12(b), it is seen that the largest extent of the chaotic sea does not warrant the shortest mixing time τm for a given Reynolds number.

Next, we consider the typical behavior near hyperbolic points. The stable and unstable manifolds of two hyperbolic points forming a heteroclinic connection coincide. This applies, e.g., to the hyperbolic points on the axis in axisymmetric vortex breakdown (e.g., see Sotiropoulos et al., 2001). For the apparently regular Lagrangian topology at Str=0.05 and Re=5 [Fig. 13(a)], the unstable manifold Wu of the hyperbolic synchronous periodic point P at (x,y)=(0.015,0.064) seems to coincide with the stable manifold Ws of the hyperbolic point Q at (x,y)=(0.011,0.405) and vice versa. Yet, due to the lack of symmetry, a thin chaotic layer is expected, which, however, could not be resolved numerically. On the other hand, for Str=0.02 and Re=10, the hyperbolic point is clearly located in the sea of chaotic trajectories [Fig. 13(b)] and the unstable manifold of P: (x,y)=(0.016,0.099) develops undulations, which become strongly amplified while approaching another hyperbolic point Q at (x,y)=(0.029,0.089) by stretching along its unstable manifold. These stretched undulations indicate transverse heteroclinic intersections, i.e., the hallmarks of chaos (Ottino, 1989; Wiggins, 1992). The increasingly rapid creation of new folds of material lines near different hyperbolic points leads to an exponential stretching rate, and the unstable manifold eventually penetrates the entire chaotic sea. Figure 13(c) shows another example for Str=0.25,Re=50, where the stable and unstable manifolds originate from a hyperbolic point of period 4.

FIG. 13.

Examples of heteroclinic connections. Dots are the stroboscopic projection of fluid trajectories to ϕ0=0. The blue lines in each figure represent an initially small circle about the hyperbolic point P [or P1 in (c)] being advected and shown after 42 (a), 6 (b), and 20 (c) periods of lid motion. Lines in cyan represent an initially small circle about the point Q [or P2 in (c)] advected backward in time. In (a), stable and unstable manifolds are indicated by green and red arrows, respectively. Black, gray, and magenta colors as in Fig. 7. (a) Str = 0.05, Re = 5, (b) Str = 0.02, Re = 10, and (c) Str = 0.25, Re = 50.

FIG. 13.

Examples of heteroclinic connections. Dots are the stroboscopic projection of fluid trajectories to ϕ0=0. The blue lines in each figure represent an initially small circle about the hyperbolic point P [or P1 in (c)] being advected and shown after 42 (a), 6 (b), and 20 (c) periods of lid motion. Lines in cyan represent an initially small circle about the point Q [or P2 in (c)] advected backward in time. In (a), stable and unstable manifolds are indicated by green and red arrows, respectively. Black, gray, and magenta colors as in Fig. 7. (a) Str = 0.05, Re = 5, (b) Str = 0.02, Re = 10, and (c) Str = 0.25, Re = 50.

Close modal

Finally, we compute the forward finite-time Lyapunov exponent (FTLE) field, as defined in Sec. 4.1 of Haller (2015), over one period of driving using the LCS Tool (Onu et al., 2015) on a grid of 100 × 100 initial positions. The FTLE is the average exponential stretching rate over the given finite time interval of an infinitesimal fluid filament aligned with the direction of maximum stretching, depending on the initial position X. The FTLE field over the time interval t[1/4,3/4] is compared to the stroboscopic projection to t0=1/4 (i.e., at the instance of zero lid velocity) for Re=500,Str=0.05, and S = 0 in Figs. 14(a) and 14(b). While KAM tori are typically associated with regions of low FTLEs, the distinction between regular and chaotic trajectories is not evident from the FTLE field over one driving period. However, ridges of the FTLE are clearly visible. Such ridges tend to outline normally repelling Lagrangian coherent structures (shrinklines) (Haller, 2015; Onu et al., 2015), which in our case correspond to stable manifolds of hyperbolic trajectories.

FIG. 14.

Stroboscopic projection to t0=1/4 (a) compared to the field of FTLE over the time interval t[1/4,3/4] (b) and t[1/4,51/4] (c) for Re=500,Str=0.05, and S = 0.

FIG. 14.

Stroboscopic projection to t0=1/4 (a) compared to the field of FTLE over the time interval t[1/4,3/4] (b) and t[1/4,51/4] (c) for Re=500,Str=0.05, and S = 0.

Close modal

In order to obtain a better approximation of the Lyapunov exponent field, we also computed the FTLE over a longer time interval of five driving periods t[1/4,51/4], but on a coarser grid of 50 × 50 initial points to reduce the computational effort. As shown in Fig. 14(c), the distinction between regular and chaotic regions now becomes much more pronounced. Yet, the time interval is still too short to reveal the chaotic nature of the trajectories near the bottom wall. The dependence of the FTLE on the time interval t1t0, where t0 and t1 are the initial and final times of the trajectories, is shown Fig. 15(a) for one regular and one chaotic trajectory. For the chaotic trajectory (red), the FTLE approaches a constant value, while for a regular trajectory (blue), it decays to zero. For long time intervals 10, a weak decay of FTLE of the chaotic trajectory is observed due to the finite initial distance between the diverging trajectories marking the endpoints of the filament and the finite size of the flow domain.

FIG. 15.

(a) Dependence of the FTLE on the time interval t1t0 for Re=500,Str=0.05, S = 0 and two initial points at t0=1/4 with (x,y)=(0.16,0.14) (blue) and (x,y)=(0.16,0.14) (red). (b) Dependence on Str of the space-averaged FTLE for (t1t0)/Str=100 and Reynolds numbers Re=1 (blue), 5 (red), 10 (yellow), 20 (purple), 50 (green), 100 (cyan), and 200 (dark red).

FIG. 15.

(a) Dependence of the FTLE on the time interval t1t0 for Re=500,Str=0.05, S = 0 and two initial points at t0=1/4 with (x,y)=(0.16,0.14) (blue) and (x,y)=(0.16,0.14) (red). (b) Dependence on Str of the space-averaged FTLE for (t1t0)/Str=100 and Reynolds numbers Re=1 (blue), 5 (red), 10 (yellow), 20 (purple), 50 (green), 100 (cyan), and 200 (dark red).

Close modal

Another measure of the mixing capability of the flow is the spatial average over the flow domain of the finite-time Lyapunov exponent FTLE¯. The result is shown in Fig. 15(b) for the time interval (t1t0)/Str=100. For low Reynolds numbers Re[1,20], the FTLE¯ decreases as the Strouhal number increases from Str=102. This effect is related to the large number of hyperbolic points existing for small Str. A comparison with Fig. 12(b) reveals, however, that for the case considered, the larger values of FTLE¯ do not guarantee a more efficient global mixing. On the other hand, for large Re[50,200], we find a maximum of FTLE¯ near Str=0.1, which is close to the minimum of τm.

In this section, we investigate how the structure of KAM tori responds to adding a constant velocity Um = SU (S > 0) to the oscillating lid velocity with amplitude U. This approach is complementary to the one of Anderson et al. (2000) who considered the steady lid-driven cavity as a reference and investigated the effect of adding an oscillating component to the lid velocity.

The change of the Lagrangian topology upon an increase in the mean velocity in the positive x direction from S = 0 to S = 3 is demonstrated in Fig. 16 by stroboscopic projections to ϕ0=0 for Re=(U+Um)L/ν=50 and Str=0.25. Even a mean lid velocity as small as S = 0.1 [Fig. 16(b)] leads to a visible breaking of the RT symmetry: While the KAM structures of the two main regular regions in Fig. 16(a) have the same qualitative structure, this is not the case for S = 0.1 shown in Fig. 16(b). In the left part of the cavity, the subharmonic tori of periods 3 and 4 are replaced by smaller subharmonic tori with higher periods. In the left part of the cavity, the synchronous tori have moved closer to the left wall. In the right part of the cavity, the tori of period 3 have moved closer to the synchronous periodic trajectory and have become embedded in the synchronous tori. Furthermore, the Lagrangian separation points on the walls (red dots) are displaced. In particular, the separation point on the moving lid and the left separation point on the bottom wall have moved to the left wall. The two upper separation points on the left wall merge and disappear together with the left system of KAM tori at S0.2 [Fig. 16(c)].

FIG. 16.

Stroboscopic projections showing the evolution the Lagrangian topology in the plane ϕ0=0 for Re=(U+Um)L/ν=50 and St/Re=0.25 as the relative magnitude S of the mean velocity of the lid increases from S = 0 to S = 3. The small red dots indicate the locations of degenerate hyperbolic points on the walls. The green and red arrows in figures (a) and (b) qualitatively symbolize their attracting and repelling directions, respectively. (a) S = 0, (b) S = 0.1, (c) S = 0.2, (d) S = 1, (e) S = 1.1, (f) S = 1.5, (g) S = 1.65, (h) S = 1.8, and (i) S = 3. (a) is reproduced with permission from L. Babor and H. C. Kuhlmann, Proc. Appl. Math. Mech. 20, e202000194 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution Licence.

FIG. 16.

Stroboscopic projections showing the evolution the Lagrangian topology in the plane ϕ0=0 for Re=(U+Um)L/ν=50 and St/Re=0.25 as the relative magnitude S of the mean velocity of the lid increases from S = 0 to S = 3. The small red dots indicate the locations of degenerate hyperbolic points on the walls. The green and red arrows in figures (a) and (b) qualitatively symbolize their attracting and repelling directions, respectively. (a) S = 0, (b) S = 0.1, (c) S = 0.2, (d) S = 1, (e) S = 1.1, (f) S = 1.5, (g) S = 1.65, (h) S = 1.8, and (i) S = 3. (a) is reproduced with permission from L. Babor and H. C. Kuhlmann, Proc. Appl. Math. Mech. 20, e202000194 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution Licence.

Close modal

As S is further increased, the remaining synchronous tori on the right side of the cavity move toward the vertical midplane, while subharmonic tori of various periodicities appear and disappear in the chaotic sea [Fig. 16(d)]. At S1, the synchronous elliptic point changes its character to a hyperbolic point of period one, associated with a bifurcation of a pair of elliptic points of period two [Fig. 16(e)]. These elliptic points depart from the hyperbolic point as S is further increased [Fig. 16(f)] until, at S1.65, the hyperbolic point is approached by a new pair of elliptic points [Fig. 16(g)]. The new elliptic points of period 2 merge with the hyperbolic point at S1.8, creating a single elliptic point of period one inside a large regular island of synchronous tori [Fig. 16(h)]. The Lagrangian topology is not significantly affected by S in the range 1.8<S<3 until, at S3, subharmonic tori of period 3 appear inside the synchronous regular island [Fig. 16(i)]. The largest extent of the chaotic sea is observed around S1.6, just before the new elliptic points are created.

The effect of S on the convective mixing time τm is shown in Fig. 17. The shrinking of the set of KAM tori in the upper left quarter of the cavity when S increases from 0 to 0.2 [Figs. 16(a)–16(c)] is accompanied by a steep increase in τm. The mixing time reaches a maximum at S = 0.3, soon after the upper left system of KAM tori vanishes. Upon a further increase in S from 0.3 to 1, the mixing time drops again to a similar value as for S = 0. Beyond S0.8, τm levels off with minimum values around S1.5.

FIG. 17.

Convective mixing time τm (dots) as function of the steady-to-oscillatory component of the lid velocity S=Um/U for Re=(U+Um)L/ν=50 and Str=0.25.

FIG. 17.

Convective mixing time τm (dots) as function of the steady-to-oscillatory component of the lid velocity S=Um/U for Re=(U+Um)L/ν=50 and Str=0.25.

Close modal

The limit S corresponds to the steady lid-driven cavity, for which the Hamiltonian system (4) is integrable (i.e., with exclusively regular trajectories) and the KAM tori coincide with the isosurfaces of the stream function ψ [Anderson et al., 2000, Fig. 3(a)]. For a finite S, the integrability is perturbed by the oscillating component of the lid motion. The periodic trajectory of period one in the upper central region of the cavity at S1.8 and the synchronous KAM tori in its neighborhood have also been observed by Anderson et al. (2000) and Horner et al. (2002). They evolve from the streamlines of the steady lid-driven cavity flow when the lid velocity is perturbed by the oscillating component. The resonance of the periodic trajectory for S = 3 [Fig. 16(i)] is briefly considered in  Appendix C.

The time-dependent incompressible flow in a square cavity driven by a harmonically moving lid has been revisited. We focused on the Lagrangian motion for a zero mean velocity of the lid where any net transport relies on fluid inertia. The properties of the Lagrangian structures were investigated in some detail, and their dependence on the governing parameters was established.

The existence and abundance of KAM tori in the Hamiltonian system (4) based on (1) depends sensitively on the Reynolds number. From Fig. 12(a), the amount of fluid occupied by chaotic streamlines tends to be larger for a given Reynolds number when the frequency is decreased. In terms of mixing efficiency, however, it is important to note that a frequency reduction also increases the timescale required to achieve the mixing. Therefore, an optimum driving frequency exists which grows with the Reynolds number. For the case considered, the mixing cannot significantly be improved by adding a steady component to the lid velocity. On the contrary, the mixing time doubles when the mean lid velocity is about 30% of the oscillation amplitude (S0.3) for Re=50 and Str=0 lid velocities S do not significantly improve on the mixing at S = 0. We expect that this behavior will also qualitatively apply to in-line mixers, where an optimal spatial period of transversal driving should exist for a given mean axial velocity of the flow.

The structure and location of KAM tori for S = 0 depends predominantly on the frequency of forcing. For relatively large Reynolds numbers and high frequencies (Str0.5), the periodic orbits inside of the remaining KAM tori are located close to the boundaries, near the singular corners of the moving lid. In such a situation, suspended particles may accumulate in or near these KAM tori due to their finite size, which causes repulsion from the moving boundary (see Romanò et al., 2019a; Wu et al., 2021; Wu et al., 2023). This property of the flow may be exploited to segregate neutrally buoyant particles by means of particle–wall interactions (Hofmann and Kuhlmann, 2011), which would be a promising extension of the present study.

Further interesting extensions of this work could be concerned with a systematic and quantitative study of the limiting cases of Str0 and Str for different Reynolds numbers. According to our observations, it is expected that for large frequencies Str1, the net Lagrangian transport becomes governed by the mean velocity field, similarly to streaming flows. On the other hand, decreasing the inertia of the flow (Str1 and Re1), we found an increasing number of hyperbolic points and KAM tori, which are extremely squeezed. In this case, the KAM tori are stretched along the streamlines of the (almost Stokesian) flow at all times. It would be interesting to investigate this limit by asymptotic perturbation methods.

Lukas Babor thanks AIC Androsch International Management Consulting GmbH for partial financial support of this work under Project No. 5168. The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Programme.

Furthermore, we thank Francesco Romanò for setting up the computations of the velocity field and sharing his Poisson solver and MATLAB packages for the computation of fluid trajectories. We are very grateful for his help and encouragement during the initial phase of this project and for stimulating discussions. We also thank Nelson Poumaëre for computing and analyzing the stroboscopic projections for the cases with non-zero mean lid velocity and Joris Marceaux for the computations of the finite-time Lyapunov exponents.

The authors have no conflicts to disclose.

Lukas Babor: Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Hendrik C. Kuhlmann: Conceptualization (lead); Project administration (lead); Resources (equal); Supervision (lead); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

To provide the velocity field at any point xD and at any phase ϕ=2πt, the nodal values uijn:=u(xi,yj,ϕn) returned by the flow solver at 1000 equidistant times tn were interpolated element-wise in space using Lagrange polynomials

un(N)(x,y)=i,j=1N+1uijnli(x)lj(y),li(x)=k=1,kiN+1xxkxixk
(A1)

of the same order N as used by the flow solver. Interpolation in time is then accomplished by a M-point Fourier interpolant in Lagrange form

u(N,M)(x,y,ϕ)=k=1Mu1000k/M(N)(x,y)Msin[M2(ϕϕ1000k/M)]cot[12(ϕϕ1000k/M)],
(A2)

as given by Kopriva (2009). The relative temporal interpolation error

ϵ(M)=||u(M)(xi,yj,ϕn)un(xi,yj)||2/||un(xi,yj)||2
(A3)

was of O(107) if the number of interpolation time-steps M per forcing period ranged from 20 for Re50 to 100 for Re=500 when Str=0.01, as shown in Fig. 18. For Str=1 and M > 20, the interpolation error was of O(106), but it was localized at the first time steps after the flow was restarted from the time-periodic state. At these time steps, NEK5000 applies a lower-order time-stepping scheme until sufficient time levels for the third-order multistep method have been generated. Furthermore, improvement of accuracy could, in this case, be achieved by restarting the computations from more history points, e.g., using the checkpointing routines of Peplinski (2016).

FIG. 18.

Relative error ϵ according to Eq. (A3) of the Fourier interpolation (A2) as a function of the number of temporal nodes M per period used for the interpolation of the flow field for Str=0.01(),Str=0.1 (red, ), Str=1 (blue, ×) and Re=500 (black solid line), Re=200 (dashed line), Re=50 (dotted line). S = 0.

FIG. 18.

Relative error ϵ according to Eq. (A3) of the Fourier interpolation (A2) as a function of the number of temporal nodes M per period used for the interpolation of the flow field for Str=0.01(),Str=0.1 (red, ), Str=1 (blue, ×) and Re=500 (black solid line), Re=200 (dashed line), Re=50 (dotted line). S = 0.

Close modal

Unless the prescribed velocity at the moving lid is regularized, the flow field exhibits very large gradients near the corners made by the lid and the stationary side walls, i.e., at x(1)=(1/2,1/2) and x(2)=(1/2,1/2), because the boundary conditions are discontinuous. The polynomial interpolation of the flow field returned by the solver exhibits Gibbs oscillations close to these corners, which makes the problem (4) stiff. This leads to prohibitively small time steps in the Dormand–Prince algorithm. The velocity fields for each time step are post-processed before fluid trajectories are computed to reduce the Gibbs phenomenon. To that end, the leading-order term of the asymptotic expansion of the singular flow field provided by Gupta et al. (1981) is employed

(A4)
ψ0(r,θ)=4ru(y=1/2)π24(π2+2πθ4sinθθcosθ),
(A4a)
where r and θ are polar coordinates

r(1,2)=||xx(1,2)||2,θ(1,2)=sign(x(1,2))arctan[(yy(1,2))/(xx(1,2))]
(A4b)

centered at and depending on the corner (1,2) to which (A4a) is applied. Within the spectral element containing the singularity, the asymptotic solution (A4) is subtracted from the nodal values obtained from the flow solver, and the residual velocity field is interpolated in space. The interpolated velocity field within the affected spectral elements is then constructed as a superposition of the interpolated residual field and the leading-order approximation (A4). Figure 19 shows this approach significantly reduces the Gibbs oscillations, such that the interpolated flow field computed on 20 × 20 elements approaches the nodal values computed on a finer grid of 40 × 40 elements. If either a non-Newtonian fluid or a leakage of fluid between the lid and the solid walls is considered, the solutions of Riedler and Schneider (1983) should be employed.

FIG. 19.

Interpolated profile of the horizontal velocity component u at x=0.497 in the top left spectral element, adjacent to the corner (x,y)=(1/2,1/2). Symbols denote nodal values computed on a uniform grid of 20 × 20 elements () and on a non-uniform grid of 40 × 40 elements refined toward the lid (×). The dashed line is the Lagrange interpolation of the nodal values (), while the full line represents the leading-order term of an asymptotic expansion of the flow near the singular corner at (x,y)=(1/2,1/2) (Gupta et al., 1981) plus the interpolated residual flow (see the text).

FIG. 19.

Interpolated profile of the horizontal velocity component u at x=0.497 in the top left spectral element, adjacent to the corner (x,y)=(1/2,1/2). Symbols denote nodal values computed on a uniform grid of 20 × 20 elements () and on a non-uniform grid of 40 × 40 elements refined toward the lid (×). The dashed line is the Lagrange interpolation of the nodal values (), while the full line represents the leading-order term of an asymptotic expansion of the flow near the singular corner at (x,y)=(1/2,1/2) (Gupta et al., 1981) plus the interpolated residual flow (see the text).

Close modal

To quantify the capability of the flow to homogenize the spatial distribution of initially unmixed portions of fluid, we use the concept of coarse-grained density as in Aref et al. (2017). To that end, we initialize a square of side length 1/2 of marked fluid in the center of the cavity and let it be advected by the flow [Fig. 20(a)]. Each side of the square is initially discretized into 500 fluid elements, which are advected according to Eq. (4) with the methods described in Sec. III B 3. The discretization is adaptively refined using the aforementioned algorithm of Meunier and Villermaux (2010). The computation is terminated either when one of the advected sides of the initial square is refined into more than 50 000 fluid elements or when the convective time t/Str from the start of the advection reaches 100. Moreover, we divide the cavity into Nδ=400 square cells of side length δ=0.05 and area Sδ=δ2.

FIG. 20.

Advection of a portion of fluid delimited by an initial square during eight periods of lid oscillation for Re=50 and Str=0.25. (a) Initial (black) and final (red) shape of the advected portion, (b) distribution of the coarse-grained density Dn corresponding to the final state, indicated by gray scale from 0 (white) to 1 (black). (c) Evolution of the intensity of segregation I (dots) and its least squares exponential fit (line).

FIG. 20.

Advection of a portion of fluid delimited by an initial square during eight periods of lid oscillation for Re=50 and Str=0.25. (a) Initial (black) and final (red) shape of the advected portion, (b) distribution of the coarse-grained density Dn corresponding to the final state, indicated by gray scale from 0 (white) to 1 (black). (c) Evolution of the intensity of segregation I (dots) and its least squares exponential fit (line).

Close modal

Let Sb(n) be the area of intersection of the advected shape of marked fluid with the nth cell. Then, the proportion Dn(t)=Sb(n)(t)/Sδ of marked fluid inside each bin is calculated once per half a period, i.e., with Δt=1/2 [Fig. 20(b)]. With the mean values

D(t)=1Nδn=1NδDn(t)andD2(t)=1Nδn=1NδDn2(t),
(B1)

one obtains the variance D2D2, which is normalized to obtain the coarse-grained intensity of segregation

I(t)=D2D2D(1D)[0,1].
(B2)

The evolution of I(t), which varies between 1 for a perfectly unmixed and 0 for a perfectly mixed state, was fitted with the function et/(τmStr) [Fig. 20(c)], where the parameter τm is a characteristic mixing time in convective scaling.

Gorodetskyi et al. (2012) have shown that the coarse-graining corresponds to a numerical diffusion acting on the length scale δ, which homogenizes the fluid composition within each cell. This numerical diffusion can be matched to an effective physical diffusion, depending on δ. τm is thus a good approximation of the characteristic mixing time of a diffusive scalar with such an effective diffusivity (it captures the effect of the dominant eigenvalue of the advection–diffusion operator).

For Re=50,Str=0.25, and S = 3, we observe in [Fig. 16(i)] a resonance τ(r¯=0)=p/q=3/1 of the synchronous periodic trajectory in the upper-central part of the cavity. Since the winding time τ(0) of a periodic trajectory cannot be computed by Eq. (15), because the winding angle θ is not well defined, it is extrapolated from the nearby trajectories

τ(0)=limr¯0τ(r¯).
(C1)

The structure of the KAM tori near the periodic trajectory at the resonance resembles Figs. 5–8 of Ishii et al. (2012) who also observe 3/1 resonances of closed streamlines in three-dimensional steady lid-driven cavities with spanwise aspect ratio Λ=(zmaxzmin)/W=1; aspect ratios Γ1=0.4, 0.6, 1, and 1.4; and Reynolds numbers Re=Rem=164, 144, 222, and 343, respectively. They have shown that the trajectories in a neighborhood of the closed trajectory at resonance can be qualitatively reproduced by three leading terms of the Hamiltonian

Hp,q(r,ϑ)=δr+Brp/2cos[p(ϑ+2πtqp+ϑ0)]+r2A(r)
(C2)

given in Chap. 8.4.3 of Arnold et al. (2007), where (r,ϑ) are scaled symplectic polar coordinates centered at the closed trajectory

xx(0)ϵ=2rcosϑ,yy(0)ϵ=2rsinϑ,
(C3)

where ϵ is the scale, A(r)=a0+a1r+a2r2+O(r3) is a generic power series in r, B and ϑ0 are free parameters, and δ=[qτ(0)]1p1 quantifies any mismatch between τ(0) and the resonant winding time p/q. For δ0, resonant tori at some distance r > 0 from the closed trajectory are broken into subharmonic tori, but small non-resonant synchronous tori remain in the vicinity of the closed trajectory [e.g., see the right system of tori in Fig. 16(b)]. On the other hand, for δ = 0, the outermost subharmonic tori meet at the closed trajectory. Following Ishii et al. (2012), selecting B = 1 and taking into account only the leading-order term A(r)=a0=1 we find a good qualitative agreement between the phase portrait of Eq. (C2) for δ=0.2,ϑ0=π/6,ϵ=0.1 (Fig. 21) and the neighborhood of the resonant periodic trajectory in Fig. 16(i). For S = 1.1 and S = 1.65 [Figs. 16(e) and 16(g)], we find another resonance τ(r¯=0)=p/q=2/1.

FIG. 21.

Isolines of the Hamiltonian (C2) for the parameters A(r) = 1, B = 1, δ=0.2,ϵ=0.1, and ϑ0=π/6.

FIG. 21.

Isolines of the Hamiltonian (C2) for the parameters A(r) = 1, B = 1, δ=0.2,ϵ=0.1, and ϑ0=π/6.

Close modal
1.
Albensoeder
,
S.
,
Kuhlmann
,
H. C.
, and
Rath
,
H. J.
, “
Three-dimensional centrifugal-flow instabilities in the lid-driven cavity problem
,”
Phys. Fluids
13
,
121
135
(
2001
).
2.
Anderson
,
P. D.
,
Galaktionov
,
O. S.
,
Peters
,
G. W. M.
,
van de Vosse
,
F. N.
, and
Meijer
,
H. E. H.
, “
Chaotic fluid mixing in non-quasi-static time-periodic cavity flows
,”
Int. J. Heat Fluid Flow
21
,
176
185
(
2000
).
3.
Aref
,
H.
,
Blake
,
J. R.
,
Budišić
,
M.
,
Cardoso
,
S. S. S.
,
Cartwright
,
J. H. E.
,
Clercx
,
H. J. H.
,
El Omari
,
K.
,
Feudel
,
U.
,
Golestanian
,
R.
,
Gouillart
,
E.
,
van Heijst
,
G. F.
,
Krasnopolskaya
,
T. S.
,
Le Guer
,
Y.
,
MacKay
,
R. S.
,
Meleshko
,
V. V.
,
Metcalfe
,
G.
,
Mezić
,
I.
,
de Moura
,
A. P. S.
,
Piro
,
O.
,
Speetjens
,
M. F. M.
,
Sturman
,
R.
,
Thiffeault
,
J.-L.
, and
Tuval
,
I.
, “
Frontiers of chaotic advection
,”
Rev. Mod. Phys.
89
,
025007
(
2017
).
4.
Arnold
,
V. I.
, “
Small denominators and problems of stability in classical and celestial mechanics
,”
Russ. Math. Surv.
18
,
85
191
(
1963
).
5.
Arnold
,
V. I.
,
Kozlov
,
V. V.
, and
Neishtadt
,
A. I.
,
Mathematical Aspects of Classical and Celestial Mechanics
(
Springer
,
2007
), Vol.
3
.
6.
Bajer
,
K.
and
Moffatt
,
H. K.
, “
Chaos associated with fluid inertia
,” in
Topological Aspects of the Dynamics of Fluids and Plasmas
(
Springer Netherlands
,
Dordrecht
,
1992
), pp.
517
534
.
7.
Biagioni
,
V.
,
Venditti
,
C.
,
Adrover
,
A.
,
Giona
,
M.
, and
Cerbelli
,
S.
, “
Taming Taylor–Aris dispersion through chaotic advection
,”
J. Chromatogr. A
1673
,
463110
(
2022
).
8.
Birkhoff
,
G. D.
,
Nouvelles recherches sur les systèmes dynamiques
(
Ex aedibus academicis in Civitate Vaticana
,
1934
).
9.
Blackburn
,
H. M.
and
Lopez
,
J. M.
, “
The onset of three-dimensional standing and modulated travelling waves in a periodically driven cavity flow
,”
J. Fluid Mech.
497
,
289
317
(
2003
).
10.
Botella
,
O.
and
Peyret
,
R.
, “
Benchmark spectral results on the lid-driven cavity flow
,”
Comput. Fluids
27
,
421
433
(
1998
).
11.
Castelain
,
C.
,
Mokrani
,
A.
,
Guer
,
Y. L.
, and
Peerhossaini
,
H.
, “
Experimental study of chaotic advection regime in a twisted duct flow
,”
Eur. J. Mech. B
20
,
205
232
(
2001
).
12.
Chien
,
W.-L.
,
Rising
,
H.
, and
Ottino
,
J. M.
, “
Laminar mixing and chaotic mixing in several cavity flows
,”
J. Fluid Mech.
170
,
355
377
(
1986
).
13.
Contreras
,
P. S.
,
Speetjens
,
M. F. M.
, and
Clercx
,
H. J. H.
, “
Lagrangian transport in a class of three-dimensional buoyancy-driven flows
,”
J. Fluid Mech.
832
,
5
40
(
2017
).
14.
Dormand
,
J. R.
and
Prince
,
P. J.
, “
A family of embedded Runge–Kutta formulae
,”
J. Comput. Appl. Math.
6
,
19
26
(
1980
).
15.
Fischer
,
P. F.
,
Lottes
,
J. W.
, and
Kerkemeier
,
S. G.
, see http://nek5000.mcs.anl.gov for the official Nek5000 web page (
2008
).
16.
Franjione
,
J. G.
,
Leong
,
C.-W.
, and
Ottino
,
J. M.
, “
Symmetries within chaos: A route to effective mixing
,”
Phys. Fluids A
1
,
1772
1783
(
1989
).
17.
Gelfgat
,
A. Y.
, “
Implementation of arbitrary inner product in the global Galerkin method for incompressible Navier–Stokes equations
,”
J. Comput. Phys.
211
,
513
530
(
2006
).
18.
Gorodetskyi
,
O.
,
Giona
,
M.
, and
Anderson
,
P. D.
, “
Spectral analysis of mixing in chaotic flows via the mapping matrix formalism: Inclusion of molecular diffusion and quantitative eigenvalue estimate in the purely convective limit
,”
Phys. Fluids
24
,
073603
(
2012
).
19.
Gupta
,
M. M.
,
Manohar
,
R. P.
, and
Noble
,
B.
, “
Nature of viscous flows near sharp corners
,”
Comput. Fluids
9
,
379
388
(
1981
).
20.
Haller
,
G.
, “
Lagrangian coherent structures
,”
Annu. Rev. Fluid Mech.
47
,
137
162
(
2015
).
21.
Han
,
S.
,
Zhang
,
S.
, and
Zhang
,
H.
, “
A Lagrangian criterion of unsteady flow separation for two-dimensional periodic flows
,”
Appl. Math. Mech.
39
,
1007
1018
(
2018
).
22.
Hansen
,
E. B.
and
Kelmanson
,
M. A.
, “
An integral equation justification of the boundary conditions of the driven-cavity problem
,”
Comput. Fluids
23
,
225
240
(
1994
).
23.
Helleman
,
R. H.
, “
Self generated chaotic behavior in nonlinear mechanics
,” in
Fundamental Problems in Statistical Mechanics
, 5th ed., edited by
E. G. D.
Cohen
(
North-Holland Publishing Company
,
Amsterdam
,
1980
), pp.
165
275
.
24.
Hofmann
,
E.
and
Kuhlmann
,
H. C.
, “
Particle accumulation on periodic orbits by repeated free surface collisions
,”
Phys. Fluids
23
,
072106
(
2011
).
25.
Horner
,
M.
,
Metcalfe
,
G.
,
Wiggins
,
S.
, and
Ottino
,
J. M.
, “
Transport enhancement mechanisms in open cavities
,”
J. Fluid Mech.
452
,
199
229
(
2002
).
26.
Ishii
,
K.
,
Ota
,
C.
, and
Adachi
,
S.
, “
Streamlines near a closed curve and chaotic streamlines in steady cavity flows
,”
Proc. IUTAM
5
,
173
186
(
2012
).
27.
Jana
,
S. C.
,
Metcalfe
,
G.
, and
Ottino
,
J. M.
, “
Experimental and computational studies of mixing in complex Stokes flows: The vortex mixing flow and multicellular cavity flows
,”
J. Fluid Mech.
269
,
199
246
(
1994
).
28.
Kalita
,
J. C.
,
Biswas
,
S.
, and
Panda
,
S.
, “
Finiteness of corner vortices
,”
Z. Angew. Math. Phys.
69
,
37
(
2018
).
29.
Kang
,
T. G.
and
Anderson
,
P. D.
, “
The effect of inertia on the flow and mixing characteristics of a chaotic serpentine mixer
,”
Micromachines
5
,
1270
1286
(
2014
).
30.
Kolmogorov
,
A. N.
, “
On conservation of conditionally periodic motions under small perturbations of the Hamiltonian
,”
Dokl. Akad. Nauk SSSR
98
,
527
530
(
1954a
).
31.
Kolmogorov
,
A. N.
, “
The general theory of dynamical systems and classical mechanics
,” in
Proceedings of the 1954 Congress in Mathematics
(North-Holland, Amsterdam,
1954b
), pp.
315
333
.
32.
Kopriva
,
D. A.
,
Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers
(
Springer Science & Business Media
,
2009
).
33.
Kuhlmann
,
H. C.
, and
Albensoeder
,
S.
, “
Stability of the steady three-dimensional lid-driven flow in a cube and the supercritical flow dynamics
,”
Phys. Fluids
26
,
024104
(
2014
).
34.
Kuhlmann
,
H. C.
and
Romanò
,
F.
, “
The lid-driven cavity
,” in
Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics
, Computational Methods in Applied Sciences, Vol. 50, edited by A. Gelfgat (
Springer
,
Berlin, Heidelberg
,
2018
), pp.
233
309
.
35.
Leong
,
C. W.
and
Ottino
,
J. M.
, “
Experiments on mixing due to chaotic advection in a cavity
,”
J. Fluid Mech.
209
,
463
499
(
1989
).
36.
Mengaldo
,
G.
,
De Grazia
,
D.
,
Moxey
,
D.
,
Vincent
,
P. E.
, and
Sherwin
,
S. J.
, “
Dealiasing techniques for high-order spectral element methods on regular and irregular grids
,”
J. Comput. Phys.
299
,
56
81
(
2015
).
37.
Meunier
,
P.
and
Villermaux
,
E.
, “
The diffusive strip method for scalar mixing in two dimensions
,”
J. Fluid Mech.
662
,
134
172
(
2010
).
38.
Moffatt
,
H. K.
, “
Viscous and resistive eddies near a sharp corner
,”
J. Fluid Mech.
18
,
1
18
(
1964
).
39.
Moffatt
,
H. K.
, “
Singularities in fluid mechanics
,”
Phys. Rev. Fluids
4
,
110502
(
2019
).
40.
Möser
,
J.
, “
On invariant curves of area-preserving mappings of an annulus
,” in
Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse
(Vandenhoeck & Ruprecht, Göttingen,
1962
), pp.
1
20
.
41.
O'Brien
,
V.
, “
Unsteady cavity flows: Oscillatory flat box flows
,”
J. Appl. Mech.
42
,
557
563
(
1975
).
42.
Onu
,
K.
,
Huhn
,
F.
, and
Haller
,
G.
, “
LCS tool: A computational platform for Lagrangian coherent structures
,”
J. Comput. Sci.
7
,
26
36
(
2015
).
43.
Ottino
,
J. M.
,
The Kinematics of Mixing: Stretching, Chaos, and Transport
, Cambridge Texts in Applied Mathematics (
Cambridge University Press
,
Cambridge
,
1989
).
44.
Ottino
,
J. M.
,
Leong
,
C. W.
,
Rising
,
H.
, and
Swanson
,
P. D.
, “
Morphological structures produced by mixing in chaotic flows
,”
Nature
333
,
419
425
(
1988
).
45.
Peplinski
,
A.
, see https://kth-nek5000.github.io/KTH_Framework/ for “
KTH framework for nek5000 toolboxes
” (
2016
).
46.
Petkovic
,
K.
,
Metcalfe
,
G.
,
Chen
,
H.
,
Gao
,
Y.
,
Best
,
M.
,
Lester
,
D.
, and
Zhu
,
Y.
, “
Rapid detection of Hendra virus antibodies: An integrated device with nanoparticle assay and chaotic micromixing
,”
Lab Chip
17
,
169
177
(
2017
).
47.
Pouransari
,
Z.
,
Speetjens
,
M. F. M.
, and
Clercx
,
H. J. H.
, “
Formation of coherent structures by fluid inertia in three-dimensional laminar flows
,”
J. Fluid Mech.
654
,
5
34
(
2010
).
48.
Riedler
,
J.
and
Schneider
,
W.
, “
Viscous flow in corner regions with a moving wall and leakage of fluid
,”
Acta Mech.
48
,
95
102
(
1983
).
49.
Romanò
,
F.
,
Kunchi Kannan
,
P.
, and
Kuhlmann
,
H. C.
, “
Finite-size Lagrangian coherent structures in a two-sided lid-driven cavity
,”
Phys. Rev. Fluids
4
,
024302
(
2019a
).
50.
Romanò
,
F.
,
Wu
,
H.
, and
Kuhlmann
,
H. C.
, “
A generic mechanism for finite-size coherent particle structures
,”
Int. J. Multiphase Flow
111
,
42
52
(
2019b
).
51.
Rüssmann
,
H.
, “
Kleine Nenner 1, Über invariante Kurven differenzierbarer Abbildung eines Kreisringes
,” in
Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse
(Vandenhoeck & Ruprecht, Göttingen,
1970
), pp.
67
105
.
52.
Shampine
,
L. F.
and
Reichelt
,
M. W.
, “
The MATLAB ODE suite
,”
SIAM J. Sci. Comput.
18
,
1
22
(
1997
).
53.
Shankar
,
P. N.
, “
The eddy structure in Stokes flow in a cavity
,”
J. Fluid Mech.
250
,
371
383
(
1993
).
54.
Shankar
,
P. N.
and
Deshpande
,
M. D.
, “
Fluid mechanics in the driven cavity
,”
Annu. Rev. Fluid Mech.
32
,
93
136
(
2000
).
55.
Sotiropoulos
,
F.
,
Ventikos
,
Y.
, and
Lackey
,
T. C.
, “
Chaotic advection in three-dimensional stationary vortex-breakdown bubbles: Šil'nikov's chaos and the devil's staircase
,”
J. Fluid Mech.
444
,
257
297
(
2001
).
56.
Speetjens
,
M.
,
Metcalfe
,
G.
, and
Rudman
,
M.
, “
Lagrangian transport and chaotic advection in three-dimensional laminar flows
,”
Appl. Mech. Rev.
73
,
030801
(
2021
).
57.
Speetjens
,
M. F. M.
,
Demissie
,
E. A.
,
Metcalfe
,
G.
, and
Clercx
,
H. J. H.
, “
Lagrangian transport characteristics of a class of three-dimensional inline-mixing flows with fluid inertia
,”
Phys. Fluids
26
,
113601
(
2014
).
58.
Takasaki
,
S.
,
Ogawara
,
K.
, and
Iida
,
S.
, “
A study on chaotic mixing in 2D cavity flows: Effects of Reynolds number and amplitude of lid velocity
,”
JSME Int. J., Ser. B
37
,
237
241
(
1994
).
59.
Taylor
,
G. I.
, see http://web.mit.edu/hml/ncfmf.html for “
Low-Reynolds-number flows
” (
1966
).
60.
Vogel
,
M. J.
,
Hirsa
,
A. H.
, and
Lopez
,
J. M.
, “
Spatio-temporal dynamics of a periodically driven cavity flow
,”
J. Fluid Mech.
478
,
197
226
(
2003
).
61.
Wiggins
,
S.
,
Chaotic Transport in Dynamical Systems
, Interdisciplinary Applied Mathematics Vol.
2
(
Springer
,
New York
,
1992
).
62.
Wu
,
H.
,
Romanò
,
F.
, and
Kuhlmann
,
H. C.
, “
Attractors for the motion of a finite-size particle in a two-sided lid-driven cavity
,”
J. Fluid Mech.
906
,
A4
(
2021
).
63.
Wu
,
H.
,
Romanò
,
F.
, and
Kuhlmann
,
H. C.
, “
Attractors for the motion of a finite-size particle in a cuboidal lid-driven cavity
,”
J. Fluid Mech.
955
,
A16
(
2023
).
64.
Zhu
,
J.
,
Holmedal
,
L. E.
,
Wang
,
H.
, and
Myrhaug
,
D.
, “
Vortex dynamics and flow patterns in a two-dimensional oscillatory lid-driven rectangular cavity
,”
Eur. J. Mech. B
79
,
255
269
(
2020
).