Lagrangian transport in the time-periodic two-dimensional lid-driven square cavity

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 U_m, 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.7≲Strm≲1$⁠, 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 U_m = 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+U cos (Ωt)]ex$⁠, where U is the amplitude of the lid velocity oscillation, Ω its circular frequency, and U_m 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

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

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

The Strouhal number $Str$ is used here to measure the driving frequency with respect to the convective timescale. Alternatively, the Stokes number $St=Re Str=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

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)

which is the collection of subsequent positions $Φk(X)$ at discrete time levels $t=t0+k$ of an initial position X at t = t₀. 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 t₀. 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×10−3,α=15Δl/π$ and κ is the local curvature of the interpolation. The refinement algorithm is described in more detail by Meunier and Villermaux (2010).

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 (⁠$PN−PN−2$ 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=10−3$ was used for all computations, leading to a Courant number $C≤0.3$⁠.

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

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(10−6)$ for $Re=100$⁠.

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

Figure 2 shows the evolution of R_k 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<ϵ=10−7$⁠, 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 R_k has dropped below the tolerance level of $ϵ=10−7$⁠, 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.

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

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 (⁠$S−1=0$⁠), this regularization has no significant effect on the flow in the bulk of the cavity.

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

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/∂x−∂u/∂y$ was computed by the subroutine comp_vort3 of NEK5000. We also consider the mean flow

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.

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×10−7$ 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.

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.

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 $Re→0$⁠, 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

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.,

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,

For vanishing lid frequency $Str→0$⁠, the partial derivative term $Str∂u/∂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 $λ=St−1/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.

Even for a zero mean velocity of the lid S = 0, there exists a non-zero mean flow $u¯$ (10) for $Re≠0$⁠. 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.

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.

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.

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 $Str−1$⁠. 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 $Str≈5.9×10−2$⁠. 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 $Str≪1$⁠, 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 $Re→0$⁠, 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 $x0≠0$⁠. 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 $Str→0$⁠. The slightly asymmetric locations of the hyperbolic points around $x≈0$ 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 $Re≲50$ [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 $Str−1≫1$ 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

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

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

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.

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.

For low driving frequencies $Str≲0.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 $Str≳0.25$ and $Re≳10$⁠, 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 $Re≳340$⁠. The evolution with $Re$ is shown in Fig. 11.

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 $Re≳100$ and $Str≲0.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).

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.

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.

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,5−1/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 $t1−t0$⁠, where t₀ and t₁ 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.

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 $(t1−t0)/Str=100$⁠. For low Reynolds numbers $Re∈[1,20]$⁠, the $FTLE¯$ decreases as the Strouhal number increases from $Str=10−2$⁠. 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 U_m = 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 $S≈0.2$ [Fig. 16(c)].

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 $S≈1$⁠, 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 $S≈1.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 $S≈1.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 $S≈3$⁠, subharmonic tori of period 3 appear inside the synchronous regular island [Fig. 16(i)]. The largest extent of the chaotic sea is observed around $S≈1.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 $S≈0.8$⁠, τ_m levels off with minimum values around $S≈1.5$⁠.

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 $S≳1.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 $(S≈0.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 (⁠$Str≳0.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 $Str→0$ and $Str→∞$ for different Reynolds numbers. According to our observations, it is expected that for large frequencies $Str≫1$⁠, 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 (⁠$Str≪1$ and $Re≪1$⁠), 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 $x∈D$ and at any phase $ϕ=2πt$⁠, the nodal values $uijn:=u(xi,yj,ϕn)$ returned by the flow solver at 1000 equidistant times t_n were interpolated element-wise in space using Lagrange polynomials

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

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

was of $O(10−7)$ if the number of interpolation time-steps M per forcing period ranged from 20 for $Re≤50$ 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(10−6)$⁠, 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).

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

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.

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$⁠.

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

one obtains the variance $⟨D2⟩−⟨D⟩2$⁠, which is normalized to obtain the coarse-grained intensity of segregation

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 $e−t/(τ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

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 $Λ=(zmax−zmin)/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

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

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)]−1−p−1$ 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$⁠.