2.1. Description of the flow

A two-dimensional differentially heated cavity is a closed container whose two vertical walls are subject to different homogeneous temperatures. The temperature of the hot and cold walls are denoted by ${\it\theta}_{h}$ and ${\it\theta}_{c}$ , respectively (see figure 1). We refer to $x$ and $z$ as the horizontal and vertical directions, respectively. The origin of the system $(x,z)=(0,0)$ is set as the bottom left corner of the cavity. The gravity $g$ is parallel to the heated walls and points downwards in the $z$ direction. This system is characterised by its height $H$ and width $W$ , and the aspect ratio of the cavity is defined as $A=H/W$ . The difference in temperature between the hot and the cold wall is defined as ${\rm\Delta}{\it\theta}={\it\theta}_{h}-{\it\theta}_{c}$ .

Figure 1.

The differentially heated cavity with geometrical aspect ratio $A=2$ for the steady state at $\mathit{Ra}=10^{8}$ : (a) isovalues of the temperature field ${\it\theta}$ from hot (top, red online) to cold (bottom, blue online); (b) streamlines. Boundary conditions for the velocity are $\boldsymbol{U}=0$ at each walls. Gravity points downwards in the $z$ direction.

Following the classical Boussinesq hypothesis, we assume that the thermophysical coefficients of the fluid are constant except for the density in the buoyancy force, which expands as

where ${\it\theta}_{0}=({\it\theta}_{h}+{\it\theta}_{c})/2$ and ${\it\rho}_{0}=\text{const}$ . Defining ${\it\kappa}$ as the thermal diffusivity of the fluid, ${\it\nu}$ as the kinematic viscosity and ${\it\beta}$ as the thermal expansion coefficient, a Rayleigh number based on the cavity height is defined as

and the Prandtl number as $\mathit{Pr}={\it\nu}/{\it\kappa}$ .

Following the nomenclature of Xin & Le Quéré (Reference Xin and Le Quéré2006), the reference length is chosen as $H$ , the reference velocity as $V_{ref}={\it\kappa}\mathit{Ra}^{1/2}/H$ and the reference pressure as $P_{ref}={\it\rho}_{0}V_{ref}^{2}{\it\nu}/{\it\kappa}={\it\rho}_{0}g{\it\beta}{\rm\Delta}{\it\theta}H$ . The convective time unit is defined by $t_{conv}=L_{ref}/V_{ref}=H^{2}/({\it\kappa}\mathit{Ra}^{0.5})$ . We use a reduced dimensionless temperature ${\it\Theta}=({\it\theta}-{\it\theta}_{0})/{\rm\Delta}{\it\theta}$ . Using dimensionless variables $\boldsymbol{U}$ for the velocity field, $P$ for the pressure, ${\it\Theta}$ for the temperature, $X$ for the horizontal direction and $Z$ for the vertical direction, the Navier–Stokes equations read

Boundary conditions for the temperature are set as adiabatic for the two horizontal walls (Neumann conditions) and isothermal (Dirichlet conditions) for each vertical wall:

The velocity field obeys the no-slip condition $\boldsymbol{U}=0$ at all walls. The bounded fluid domain defined by the cavity is referred to as $\mathscr{D}=\{\boldsymbol{X}=(X,Z),X\in [0,A^{-1}],Z\in [0,1]\}$ . For this study we fix $\mathit{Pr}=0.71$ , which corresponds to air, and $A=2$ . The Rayleigh number $\mathit{Ra}$ is the only control parameter in our study.

2.2. Numerical method

DNS was performed using a code written by S. Xin and P. Le Quéré. It is based on a pseudo-spectral Chebyshev collocation method for the spatial discretisation and a direct Uzawa method for the velocity–pressure coupling. The time-stepping algorithm combines a second-order Adams–Bashforth extrapolation for the convective terms with a second-order backward differentiation formula for the linear ones. The resulting Helmholtz equations are solved by full diagonalisation of the second-order partial derivatives. Based on previous works, for all $\mathit{Ra}$ values presented in this paper, we used 81 and 101 collocation points in the $X$ and $Z$ directions, respectively, with a timestep of ${\rm\Delta}t=10^{-3}$ . This resolution is sufficient to accurately predict the first Hopf bifurcations of the system (Xin & Le Quéré Reference Xin and Le Quéré2006). Accurate convergence towards steady states below and beyond $\mathit{Ra}_{c}$ (where the steady state is unstable) is achieved using the selective frequency damping algorithm (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006), with the convergence criterion $\Vert \boldsymbol{u}(t)-\boldsymbol{u}(t-{\rm\Delta}t)\Vert _{\infty }<10^{-10}$ and $\boldsymbol{u}=(U,W,{\it\Theta},P)$ .

2.3. Steady base flow

Due to buoyancy, fluid rises along the hot wall and sinks along the cold wall irrespectively of the value of $\mathit{Ra}$ . As $\mathit{Ra}$ increases, vertical boundary layers narrow while a homogeneous stratification appears in the core of the cavity for $\mathit{Ra}\gtrsim 10^{6}$ . At $\mathit{Ra}\approx 10^{7}$ detached recirculations form in the vicinity of two of the cavity corners. In addition, for all times the base flow $\boldsymbol{u}$ respects the $Z_{2}$ inversion symmetry

where the operators $R$ and ${\it\zeta}$ are defined by

The base flow is shown in figure 1 for $\mathit{Ra}=10^{8}$ . The stratified core of the base flow supports internal gravity waves. The relaxation time of the least-damped internal waves varies as $0.1\times \mathit{Ra}^{1/2}$ , which is typically $10^{3}$ convective time units for $\mathit{Ra}\approx 10^{8}$ (Le Quéré & Behnia Reference Le Quéré and Behnia1998).

2.4. Onset of unsteadiness

The first Hopf bifurcation of the Eulerian system occurs at $\mathit{Ra}_{c}=1.5865\times 10^{8}$ (Xin & Le Quéré Reference Xin and Le Quéré2006). The bifurcated state exhibits oscillations of the two detached areas in phase opposition, together with an internal mode oscillating less vigorously in the core of the cavity. It is characterised by an angular frequency ${\it\omega}_{0}=2{\rm\pi}/T_{0}\simeq 0.28225$ , with the period $T_{0}$ only weakly varying with $\mathit{Ra}$ . Four different phases ${\it\Phi}=2{\rm\pi}t/T_{0_{[2{\rm\pi}]}}$ of the corresponding global mode are shown in figure 2 for $\mathit{Ra}=1.6\times 10^{8}$ . As reported in Burroughs et al. (Reference Burroughs, Romero, Lehoucq and Salinger2002), the corresponding solution $\boldsymbol{v}(\boldsymbol{X},t)$ of the linearised problem exhibits the anti-symmetry property

while the full nonlinear field does not (Le Quéré & Behnia Reference Le Quéré and Behnia1998). For $\mathit{Ra}\simeq 2.07\times 10^{8}$ , the time-periodic state along this branch undergoes a secondary Hopf bifurcation. The new bifurcated state is characterised by the additional propagation of rolls confined to the vertical boundary layers, with an angular frequency ${\it\omega}_{2}\simeq 3.0919$ . The resulting regime is quasiperiodic in time and will not be investigated here. We however report the value of $\mathit{Ra}$ at which it bifurcates since it marks the upper limit in $\mathit{Ra}$ of the time-periodic regime.

Figure 2.

Departure from steady state of the velocity field for $\mathit{Ra}=1.6\times 10^{8}$ . The steady state has been obtained with the SFD algorithm: (a) $t=T_{0}$ ; (b)  $t=T_{0}+T_{0}/4$ ; (c) $t=T_{0}+2T_{0}/4$ ; (d) $t=T_{0}+3T_{0}/4$ .

3.1. Particle tracking: numerical method

Our Lagrangian approach is based on post-processing the velocity field obtained from the DNS computations. Whereas only one velocity field needs to be stored in the steady regime, in the oscillatory regime the period $T_{0}$ is split into approximatively $2\times 10^{4}$ snapshots. Usually, $T_{0}$ is a priori not an integer multiple of ${\rm\Delta}t$ , consequently, a fifth-order polynomial interpolation is performed to compute the last timestep of the period.

In the Lagrangian formulation, the motion of a non-diffusive tracer is described by

where $\boldsymbol{U}$ is the solution of (2.3). Here $\boldsymbol{U}$ fields are extracted from DNS and fifth-order Chebyshev polynomial interpolation is used to determine the velocity field $\boldsymbol{U}$ at any arbitrary spatial location. Owing to the post-processing approach, spatial interpolation leads to faster computations than a Galerkin evaluation of the velocity field. The position of the particle $\boldsymbol{X}(t)$ is updated using a fourth-order Runge–Kutta scheme using as a timestep ${\rm\Delta}t=10^{-3}$ (identical to the timestep in the DNS). The accuracy of the scheme is assessed by quantifying the closure of streamlines in the steady regime for $\mathit{Ra}=2.05\times 10^{8}$ . While streamlines should be closed periodic trajectories, with our parameters the average closure error is in practice below $5\times 10^{-7}$ and around $10^{-4}$ for streamlines around saddle points.

3.2. Tracers in steady velocity fields

In an incompressible two-dimensional flow, (3.1) combined with the divergence-free condition (2.5) leads to

where ${\it\Psi}$ is the streamfunction. The system (3.2) has a Hamiltonian structure. When ${\it\Psi}={\it\Psi}_{0}(X,Z)$ is time-independent, the system (3.2) is integrable. In this case, the trajectories of the tracers coincide with the isocontours of ${\it\Psi}$ (i.e. streamlines) and chaotic advection does not occur.

We first focus on the topology of the steady velocity field displayed in figure 3, which is organised around the fixed points of $\boldsymbol{U}$ . Streamlines originating from saddle fixed points of the velocity fields either form homoclinic loops, closing on the fixed point, or heteroclinic orbits connecting two different saddle fixed points together. Their period is infinite due to the slowing down in the vicinity of the fixed points. We identified in total five saddle fixed points (referred to as $\boldsymbol{X}_{S_{i}},\ i=1,\dots ,5$ ) and four centres. (Fixed points of the steady velocity field were found using the algorithm described in the appendix A.) Owing to the centro-symmetry, $\boldsymbol{X}_{S_{4}}=R\boldsymbol{X}_{S_{1}}$ and $\boldsymbol{X}_{S_{5}}=R\boldsymbol{X}_{S_{2}}$ . The three saddle points $\boldsymbol{X}_{S_{1}}$ , $\boldsymbol{X}_{S_{2}}$ and $\boldsymbol{X}_{S_{3}}$ lie near the corner vortex, at mid-height of the cavity near the hot wall and at the geometrical centre of the cavity, respectively. They are associated with three homoclinic and two heteroclinic connections. Those connections are referred to as ${\it\Gamma}_{i,j}$ where $i=1,2,3$ is the label of the saddle point $\boldsymbol{X}_{S_{i}}$ and $j=1,2$ refers to the direction considered in each unstable eigenspace.

Figure 3.

Homoclinic and heteroclinic connections of the (unstable) steady base flow for $\mathit{Ra}=1.625\times 10^{8}$ . Filled circles: saddle points; open circles: centres.

For any initial condition $\boldsymbol{X}_{0}$ belonging to a closed streamline, we can define the period $T_{SL}(\boldsymbol{X}_{0})$ of the corresponding orbit. A spatial distribution of $T_{SL}$ is obtained by seeding at $t=0$ the segment ${\bf\chi}=\{X_{0}\in ]0,1/(2A)],Z=A/2\}$ with $10^{5}$ tracers. Tracers are labelled by their initial horizontal position $X_{0}$ in ${\bf\chi}$ . For each of them, equation (3.1) is integrated up to time $T_{SL}$ after which the tracer has come back to its initial position. The distribution of $T_{SL}$ versus $X_{0}$ is shown in figure 4 for $\mathit{Ra}=1.625\times 10^{8}$ . It is smooth except for the presence of three singularities of $T_{SL}$ (referred to as $\mathit{Sg}_{i}$ , $i=1,2,3$ ). Here $\mathit{Sg}_{1}$ is related to the no-slip at the wall whereas $\mathit{Sg}_{2}$ and $\mathit{Sg}_{3}$ are the signatures of the heteroclinic and homoclinic connections $R{\it\Gamma}_{1,2}$ and ${\it\Gamma}_{2,2}$ , respectively. A last zone not displayed in figure 4 lies around $\boldsymbol{X}_{S_{5}}$ , the geometrical centre of the cavity, and is bounded by the heteroclinic connections ${\it\Gamma}_{2,1}$ and $R{\it\Gamma}_{2,1}$ . That zone features two homoclinic loops attached to $\boldsymbol{X}_{S_{5}}$ . It is associated with large values of $T_{SL}$ since the velocity magnitude is very weak.

Figure 4.

Distribution of return times for the (unstable) steady base flow $T_{SL}(X_{0})$ for $\mathit{Ra}=1.625\times 10^{8}$ near the hot wall. The left divergence ( $\mathit{Sg}_{1}$ ) is the signature of the non-slip boundary condition at the wall, whereas $\mathit{Sg}_{2}$ and $\mathit{Sg}_{3}$ are linked to a homoclinic or heteroclinic connection. Return times are computed by seeding the segment ${\bf\chi}$ with approximatively $4\times 10^{5}$ tracers.

3.3. Tracers in time-periodic velocity fields

Beyond the Hopf bifurcation, the dynamical system (3.1) becomes non-autonomous. Slightly above $Ra_{c}$ , the steady streamfunction is perturbed by a small time-dependent term, i.e.  ${\it\Psi}(X,Z,t)={\it\Psi}_{0}(X,Z)+{\it\varepsilon}{\it\delta}{\it\Psi}(X,Z,t)$ , with ${\it\varepsilon}\ll 1$ . Owing to time periodicity, ${\it\delta}{\it\Psi}(X,Z,t+T_{0})={\it\delta}{\it\Psi}(X,Z,t)$ . In contrast to the steady case, here particles initially in ${\bf\chi}$ are not expected to return exactly to their initial position. The mapping between the initial position and the position at the next crossing defines an input–output problem and can be seen as an example of scattering process. Such an approach is usually dedicated to open systems where particles leave the test domain never to return. The times of escape from the scatterer can be evaluated statistically in order to investigate, through the presence of fractal singularities, the properties of the underlying chaotic saddle (Lai & Tél Reference Lai and Tél2011). Chaotic saddles are built around an infinite set of periodic orbits around which particles can spend a long and arbitrary time before being ejected further downstream (Stolovitzky, Kaper & Sirovich Reference Stolovitzky, Kaper and Sirovich1995; Neufeld & Tél Reference Neufeld and Tél1998; Budyansky, Uleysky & Prants Reference Budyansky, Uleysky and Prants2004; Biemond et al. Reference Biemond, de Moura, Grebogi and Nijmeijer2008). While the concept of chaotic saddle does not necessarily make sense in a closed recirculating flow, we demonstrate here how the tools of chaotic scattering can shed light on the dynamics of tracers.

Figure 5.

(a) Return times $T_{P}$ as a function of $X_{0}$ and phase ${\it\Phi}$ , (b) projection on the $(X_{0},{\it\tau}=T_{P}/T_{SL})$ plane, (c) cut through $X_{0}=0.02$ (area $A_{2}$ ), for varying ${\it\Phi}$ , (d) cut through ${\it\Phi}=0$ for varying $X_{0}$ . Inset: zoom on the peaks in $A_{2}$ . The fractal distribution is obtained by seeding the segment $\{X\in [2.151\times 10^{-2},2.152\times 10^{-2}],Z={\textstyle \frac{1}{2}}\}$ with $10^{5}$ tracers. Here $\mathit{Ra}=1.625\times 10^{8}$ .

Return times are here computed for initial conditions on the segment ${\bf\chi}$ , now using the unsteady velocity field $\boldsymbol{U}(\boldsymbol{X},t)$ . However, since particle trajectories are no longer closed, except in the case of periodic orbits, the return time (now denoted $T_{P}$ ) is defined as the time required for a particle initially on ${\bf\chi}$ to come back through ${\bf\chi}$ with ${\dot{Z}}>0$ . Once back, the tracer provides a new initial condition from which a new return time can be extracted, and so forth until the next crossing. We record the return times associated with particles having crossed ${\bf\chi}$ again at a phase ${\it\Phi}$ from an initial position $X_{0}$ . Here $T_{P}$ is normalised by the return time $T_{SL}$ of the particle originating from $(X_{0},{\textstyle \frac{1}{2}})$ in the steady regime. This normalisation avoids the divergences resulting from the walls and from the homoclinic/heteroclinic connections. Such a procedure has been applied up to the 40th return through ${\bf\chi}$ for $10^{4}$ particles. Figure 5 shows the distribution of ${\it\tau}=T_{P}/T_{SL}$ as a function of ${\it\Phi}$ and $X_{0}$ . Three seemingly fractal distributions of ${\it\tau}$ are found in the regions referred to as $A_{i},\ i=1,2,3$ . They indicate the presence of chaotic scattering regions explored, be it transiently or not, by the tracers. Here $A_{1}$ corresponds to a region including the proximity of the walls. Here $A_{2}$ and $A_{3}$ correspond to the proximity of former homoclinic and heteroclinic loops ${\it\Gamma}_{1,j}$ and ${\it\Gamma}_{2,j}$ , $j=1,2$ . The regions $A_{i}$ are delimited by intermediate regions with a smooth dependence on ${\it\Phi}$ and smooth variations of ${\it\tau}$ . Such regions are referred to as $N_{i},\ i=1,2,3$ in figures 5 and 6. The central zone around the centre of the cavity is not treated here because of the huge return times associated with nearly vanishing velocities.

Figure 6.

Trajectories inside the $A_{i}$ areas (in black) associated with singularities in figure 5. Here $N_{i}$ corresponds to the trajectories inside the non-hyperbolic areas (in grey). Here $i=1,2,3$ and $\mathit{Ra}=1.625\times 10^{8}$ .

Figure 7 shows the cumulative probability distributions $P({\it\tau}>{\it\tau}_{\ast })$ for the normalised return time ${\it\tau}$ to be larger than a given value ${\it\tau}_{\ast }$ . Only events for which ${\it\tau}>2$ are taken into account. Two possible asymptotic scalings can be expected for $P({\it\tau}>{\it\tau}_{\ast })$ . An exponential decay of the form $P({\it\tau}>{\it\tau}_{\ast })\simeq \text{e}^{-{\it\alpha}{\it\tau}_{\ast }}$ would indicate hyperbolic scattering, while an algebraic decay of the form $P({\it\tau}>{\it\tau}_{\ast })\simeq {\it\tau}_{\ast }^{-{\it\lambda}}$ would indicate non-hyperbolic scattering (Ott & Tél Reference Ott and Tél1993). This distinction is crucial for a description of the long-time properties of the mixing process. In addition, it allows us to establish a classification of the different dynamics occurring in different spatial zones of the flow.

Figure 7.

Cumulative probabilities $P({\it\tau}>{\it\tau}_{\ast })$ where ${\it\tau}=T_{P}/T_{SL}$ (see the text) for different regions $A_{1}$ , $A_{2}$ , $A_{3}$ and $\mathit{Ra}=1.625\times 10^{8}$ : (a) $A_{1}$ , algebraic fit with exponent ${\it\lambda}=3.3$ ; (b) $A_{2}$ , exponential fit with ${\it\alpha}=1.15$ ; (c) $A_{3}$ , exponential fit with ${\it\alpha}=0.10$ ; (d) same as (c) in log–log plot, showing a better algebraic fit with ${\it\lambda}=4.11$ . The corresponding areas are shown in figure 5(b) and their representation in the physical space in figure 6.

In the region $A_{1}$ , $P({\it\tau}_{\ast })$ scales algebraically, see figure 7(a). This strongly suggests that the particles departing from $A_{1}$ are advected through an area where all invariant sets possess non-hyperbolic properties. The probability distribution in region $A_{2}$ almost perfectly fits an exponential decay (figure 7 b). Therefore, we can expect the dynamics to be hyperbolic. Finally, $P({\it\tau}_{\ast })$ in region $A_{3}$ shows a mixed exponential/algebraic scaling (figure 7 d). Since the algebraic fit is valid over a wider range of values of ${\it\tau}_{\ast }$ than the exponential one (see figure 7 c), it suggests that $A_{3}$ is dominated by non-hyperbolic dynamics for larger times but still displays signs of near-hyperbolic dynamics on shorter timescales, at least for the observation times considered here. As we shall see next, hyperbolic properties of $A_{2}$ and $A_{3}$ are related to the breakdown of the homoclinic and heteroclinic loops ${\it\Gamma}_{i,j}$ in the unsteady regime, whereas non-hyperbolic properties are associated with the survival of KAM tori and to a lesser degree with the no-slip condition at the walls. The next section is devoted to an explanation of all of these statistical properties from a deterministic point of view.

4.1. Definition of the stroboscopic section ${\it\Pi}$

Taking advantage of the time-periodicity of the flow field, the dynamics of the tracers can be studied using a stroboscopic map ${\rm\pi}$ transporting the position of any tracer in $\mathscr{D}$ onto its position one period $T_{0}$ later. In this representation, (3.1) is reduced to an autonomous two-dimensional map ${\rm\pi}:\mathscr{D}\rightarrow \mathscr{D}$ such that ${\rm\pi}(\boldsymbol{X}_{0})={\it\phi}^{T_{0}}(\boldsymbol{X}_{0})$ , ${\it\phi}^{T_{0}}$ being the propagator over one period $T_{0}$ . We define the stroboscopic section of a given trajectory passing through $\boldsymbol{X}_{0}$ as ${\it\Pi}=\{{\rm\pi}^{n}(\boldsymbol{X}_{0})\}_{n\in \mathbb{Z}}$ . Note that since the flow is incompressible, the continuous-time flow ${\it\phi}^{t}$ is area-preserving, hence the initial phase ${\it\Phi}$ of the velocity field can be chosen arbitrarily in $[0,2{\rm\pi}[$ to define ${\it\Pi}$ . Furthermore, the harmonic forcing is not external but intrinsic to the flow since it is determined by the Hopf bifurcation occurring at $\mathit{Ra}=\mathit{Ra}_{c}$ .

4.2. Numerical tools in the stroboscopic section

An orbit crossing the stroboscopic section $n$ times before returning to its original position defines $n$ periodic points of period $n$ in ${\it\Pi}$ . The numerical identification of such periodic points is based on two numerical algorithms: a standard Newton–Raphson algorithm, and an alternative algorithm described in the appendix A. The convergence criterion for both algorithms is $\Vert {\rm\pi}^{n}(\boldsymbol{X})-\boldsymbol{X}\Vert _{\infty }<10^{-10}$ . A Jacobian matrix is evaluated at order one by releasing two extra tracers in the neighbourhood of the reference periodic point and monitoring the separation versus time. Depending on the imaginary part of the spectrum of this matrix, we can distinguish whether the periodic point is a centre or a saddle. Therefore, eigenvalues are not expected to be accurate at all. However, it is expected that the nature of the identified periodic point (saddle or centre) will still be reflected in the imaginary part of the eigenvalues. Once an unstable periodic point is identified, a disk of radius $10^{-7}$ seeded with 100 tracers distributed on 10 concentric circles is placed around this point. Tracers are then tracked in ${\it\Pi}$ in backward and forward time in order to highlight the stable and unstable manifolds of each hyperbolic point.

4.3. Destabilisation of homoclinic and heteroclinic connections

As long as the Lagrangian system is autonomous, the stable ( $W^{s}$ ) and unstable ( $W^{u}$ ) manifolds of a homoclinic saddle point $\boldsymbol{X}_{S}$ coincide and form a homoclinic loop ${\it\Gamma}$ . When subject to harmonic forcing, the saddle point becomes a periodic orbit, i.e. a periodic point $\boldsymbol{X}_{P}$ in ${\it\Pi}$ whose manifolds $W^{u}(\boldsymbol{X}_{P})$ and $W^{s}(\boldsymbol{X}_{P})$ no longer coincide. However, if they cross transversally at least once in ${\it\Pi}$ it is known that they cross an infinite countable number of times (Poincaré Reference Poincaré1892; Guckenheimer & Holmes Reference Guckenheimer and Holmes1983; Wiggins Reference Wiggins1988). This gives rise to a so-called homoclinic tangle consisting of an infinity of lobes of equal area. The conservation of the areas is a direct consequence of the Hamiltonian structure of (3.2) and implies a complex foliation of the lobes. Periodic points are dense in the chaotic region and the dynamics restrained to this area displays hyperbolic mixing properties. The transport from lobe to lobe is responsible for the mixing properties in the vicinity of the original homoclinic loop ${\it\Gamma}$ (Ottino Reference Ottino1990; Rom-Kedar et al. Reference Rom-Kedar, Leonard and Wiggins1990). This is illustrated in figure 8 where the manifolds are depicted over an observation time long enough to display a few transversal crossings. The same scenario is valid for the destabilisation of heteroclinic orbits.

Figure 8.

(a) Sketch of a homoclinic connection ${\it\Gamma}$ from a fixed point $\boldsymbol{X}_{S}\ ({\it\varepsilon}=0)$ . (b) Transversal intersections of the stable ( $W^{s}$ ) and unstable ( $W^{u}$ ) manifolds of the associated periodic point $\boldsymbol{X}_{P}$ lead to a homoclinic tangle ( ${\it\varepsilon}>0$ ).

Melnikov functions are a convenient tool to predict whether the stable and unstable manifolds of a periodic point cross each other or not (Melnikov Reference Melnikov1963). Such a function is proportional to the signed $O({\it\varepsilon})$ distance between the stable and unstable manifolds of the hyperbolic point:

where $\boldsymbol{U}_{0}$ corresponds to the (unstable) steady velocity field evaluated along the homoclinic loop ${\it\Gamma}$ , $\boldsymbol{X}_{0}(t)$ corresponds to a trajectory along ${\it\Gamma}$ parametrised by time $t$ , and $\boldsymbol{U}_{P}$ is the finite-amplitude perturbation from $\boldsymbol{U}_{0}$ . If $M(t_{0})$ vanishes for some isolated values of $t_{0}$ , the stable and unstable manifolds intersect transversally and a homoclinic tangle exists in the vicinity of ${\it\Gamma}$ . The mass transport from lobe to lobe associated with the homoclinic tangle produces Lagrangian mixing. The integral of the unsigned Melnikov function provides a first order approximation of the area covered by the lobes (Rom-Kedar et al. Reference Rom-Kedar, Leonard and Wiggins1990; Hackborn, Ulucakli & Yuster Reference Hackborn, Ulucakli and Yuster1997).

The Melnikov function in (4.1) can be redefined as function of the phase variable ${\it\Phi}=2{\rm\pi}t_{0}/T_{0}\in [\!0,2{\rm\pi}\![$ . In order to account for the variation of $M$ with $\mathit{Ra}$ , we furthermore consider the normalisation of the Melnikov function given by

where the subscripts $i$ and $j$ refer to the index of the homoclinic or heteroclinic loop as defined in § 3.2. The Melnikov functions are computed numerically by direct evaluation of the integral along ${\it\Gamma}$ in (4.1). Figure 9 shows the normalised Melnikov functions for the homoclinic loops ${\it\Gamma}_{1,1}$ and ${\it\Gamma}_{2,2}$ and the heteroclinic connections ${\it\Gamma}_{1,2}$ and ${\it\Gamma}_{2,1}$ , each of them computed for four values of $\mathit{Ra}>\mathit{Ra}_{c}$ . In all cases ${\it\Lambda}({\it\Phi})$ vanishes at least twice, indicating transversal intersections of the manifolds, and hence the existence of chaotic tangles via the break-up of the connections.

Figure 9.

Normalised Melnikov functions ${\it\Lambda}_{i,j}$ associated with saddle points of the unstable steady state for $\mathit{Ra}=1.587\times 10^{8}$  (—),  $\mathit{Ra}=1.590\times 10^{8}$ ( $-$   $-$ ),  $\mathit{Ra}=1.600\times 10^{8}$ ( $-\cdot$ ) and $\mathit{Ra}=1.625\times 10^{8}$ ( $\cdot$ ). The presence of isolated zeros shows the existence of homoclinic/heteroclinic tangles.

Both unstable periodic points $\boldsymbol{X}_{P_{1}}$ and $\boldsymbol{X}_{P_{2}}$ , associated with $\boldsymbol{X}_{S_{1}}$ and $\boldsymbol{X}_{S_{2}}$ respectively, give rise to the tangles shown in figure 10 for the case $\mathit{Ra}=1.587\times 10^{8}$ . It is noteworthy from figure 10 that for $\mathit{Ra}$ closest to $\mathit{Ra}_{c}$ , the distance $|\boldsymbol{X}_{P_{i}}-\boldsymbol{X}_{S_{i}}|$ is larger for $i=1$ than for $i=2$ . This is due to the intensity of the oscillations of $\boldsymbol{U}_{\boldsymbol{P}}$ being much larger near $\boldsymbol{X}_{S_{1}}$ than near $\boldsymbol{X}_{S_{2}}$ , as is already clear from figure 3. Note the similarities between figure 10(b,c) and the sketch in figure 8(b). All destabilised homoclinic/heteroclinic connections identified in figure 5 induce fluid mixing, at least in a finite-size neighbourhood. Tracers seeded in the close vicinity of such hyperbolic points appear to spread over an area comparable to that of the lobes, so that mixing is only partial for $\mathit{Ra}=1.587\times 10^{8}$ (figure 10 a). This is typical of the existence of transport barriers in an $O({\it\varepsilon})$ vicinity of the homoclinic tangle, as investigated in the following subsections. The case of higher $\mathit{Ra}$ (when ${\it\varepsilon}$ is no longer small) will be treated in the last section of the paper.

Figure 10.

(a) Spreading of particles initially placed close to the hyperbolic fixed points $X_{P_{1}}$ and $X_{P_{2}}$ . Circles, fixed points associated with the unstable steady state at $\mathit{Ra}=1.587\times 10^{8}$ ; solid lines, homoclinic and heteroclinic orbits; black dots (red online), forward time integration; grey dots (green online), backward time integration. Both disks are composed of $10^{4}$ particles and are integrated over 50 periods in ${\it\Pi}$ following the procedure described in § 4.2. (b,c) Enlarged views of the crossing manifolds originating from $X_{P_{1}}$ (b) and $X_{P_{2}}$ (c) on shorter integration times. Note that $X_{P_{2}}$ and $X_{S_{2}}$ are nearly indistinguishable.

4.4. Periodic points in the stroboscopic section

Because of incompressibility, periodic points in ${\it\Pi}$ can only be saddles or centres, while sinks are precluded. Because mixing is related to stretching by the velocity field, local mixing can only occur in the immediate vicinity of hyperbolic saddle points in ${\it\Pi}$ , which in continuous time correspond to either saddle points or unstable periodic orbits (UPOs), but not around centres. As seen in the previous subsection, the periodic points corresponding to saddle points of the base flow all induce mixing in the unsteady regime. In fact, all hyperbolic points identified in ${\it\Pi}$ have the property to mix fluid over some (a priori unknown) finite-size area. Within the mixing zone corresponding to one given periodic point $\boldsymbol{X}_{P}$ , a dense aperiodic trajectory linking any two points should exist, as a consequence of the chaotic nature of mixing. Hence, any point within that zone should be arbitrarily close to both $W^{s}(\boldsymbol{X}_{P})$ and $W^{u}(\boldsymbol{X}_{P})$ . In other words the mixing area associated with a given periodic point $\boldsymbol{X}_{P}$ in ${\it\Pi}$ coincides with the closure of its stable and unstable manifolds, themselves being interpreted as the skeleton of the hyperbolic mixing zone. This point of view is prominent in the study by Feudel et al. (Reference Feudel, Witt, Gellert, Kurths, Gregobi and Sanjuán2005) who suggest that the total area occupied by mixing in the system corresponds to the closure of a bundle of manifolds.

Here we wish to determine numerically the mixing zone associated with each periodic point that can be identified in ${\it\Pi}$ , by tracing out its stable and unstable manifolds and measuring the area covered by their intertwining. This task is ambitious since (i) there are infinitely many periodic points of arbitrary high period and (ii) the infinitely convoluted geometry of these manifolds requires in practice a coarse-graining procedure. Moreover, each mixing zone contains new periodic points whose manifolds cover the same area as the periodic point originally considered to define the mixing zone. Hence redundancies are expected in practice. The benefit of this approach is an accurate a priori prediction of the various mixing zones of the flow from a deterministic point of view, rather than from the statistical a posteriori analysis of § 3. In addition, as we shall see, in practice a finite number of periodic points is enough to obtain a quantitatively satisfying coverage of all mixing zones, and hence to determine to which degree mixing in the cavity can be considered as complete or incomplete. It is moreover expected that the chaotic area increases with $\mathit{Ra}$ and we suggest in what follows a method to quantify it.

The method relies on the computation of spreading maps associated with every individual periodic point in ${\it\Pi}$ . We first begin by extracting numerically as many periodic points as possible using the method explained in § 4.2. We introduce the disk $\mathscr{B}_{i}$ of radius $d$ centred on a periodic point $P_{i}$ of period $T_{i}$ . A finite number of tracers ( $\mathscr{B}_{i}^{\prime }$ ) is chosen in $\mathscr{B}_{i}$ as in § 4.2. The partial spreading map associated with $P_{i}$ is defined as

We later define the total spreading map as the union of all partial spreading maps for a given value of $\mathit{Ra}$ :

where $\mathscr{R}_{n}$ is the set of all identified periodic points whose period is less than $nT_{0}$ . In the present investigation, we used $d=10^{-7}$ , $Q=50$ and $n=40$ , as a compromise between accuracy and computational effort (see § 5).

Figure 11.

Two examples of partial spreading maps for $\mathit{Ra}=1.625\times 10^{8}$ . (a,d) UPOs (lines) and associated periodic points (open dots) in ${\it\Pi}$ . (b,e) partial spreading map for the periodic points in ${\it\Pi}$ . (c,f) zoom on the manifolds of the periodic points. Forward iterations shown in black (red online), backward in grey (green online). (a–d) tangle associated with a periodic point with period $T=2T_{0}$ . (d–f) non-mixing island associated with resonant streamline of period $T=6T_{0}$ . ( f) is a close-up of the rectangle in (a).

Two examples of partial spreading maps are shown in figure 11. The corresponding UPO (a,d), the associated partial spreading map (b,e) and a zoom on its stable/unstable manifolds (c,f) are successively displayed.

1. (a) The first UPO (first row) has period $2T_{0}$ and is associated with two periodic points of period two. The UPO is entirely located in one of the two detached corners. Even if it occupies a small fraction of the cavity (figure 11 a), the mixing induced by the periodic point is considered as efficient since it covers a comparatively large area that coincides with $A_{2}$ (figure 11 b). It is clear from figure 11(c) that the stable and unstable manifolds of these points cross transversally and that the scenario of § 4.3 is reproduced even in the absence of homoclinic loop.

2. (b) The second UPO shown in figure 11(d) (solid line) has period $6T_{0}$ and is associated with six periodic points in ${\it\Pi}$ . As can be seen from figure 11(e), the mixing induced by this UPO is highly inefficient even though the diameter of the UPO is comparable to the cavity size. A zoom on the topological structure around two of these periodic points is shown in figure 11( f). The vicinity of the orbit exhibits heteroclinic connections between four periodic points and the appearance of several heteroclinic tangles. However, the structure defined by the periodic points composing the heteroclinic tangles does not spread particles efficiently.

Those two examples demonstrate that the diameter of a given UPO and its mixing properties are not necessarily correlated.

4.5. Resonant streamlines

In order to highlight the dynamics inside the structure displayed in figure 11( f), we have added lines of particles inside each region lying between the saddle points represented, and propagated these new tracers in ${\it\Pi}$ over a time of $50T_{0}$ . This computation shows that if a particle starts inside the structure, it cannot escape from it. Most of those particles trace out ellipses in ${\it\Pi}$ . Equivalently the continuous-time dynamics is quasiperiodic and the associated trajectories form toroidal surfaces in a space–time representation. These tori are precisely the KAM tori resulting from the non-integrability of the Hamiltonian system (3.1). These structures are not associated with an unstable manifold and, hence, are not repelling. Tracers in a neighbourhood of KAM tori are known to ‘stick’, which results in increased residence times. As a consequence of this stickiness effect, they are commonly interpreted as the origin for non-hyperbolicity in the statistics of return times. The KAM tori define strict material barriers in ${\it\Pi}$ (Ottino Reference Ottino1990), forbidding transport from one mixing zone to another one.

In the steady regime any streamline, parametrised by its closure time $T_{SL}$ , acts as a barrier to transport. As the flow becomes oscillatory, the periodic dynamics of tracers becomes in principle quasiperiodic with two periods, one close to $T_{SL}$ and the second being $T_{0}$ . The KAM theorem (see e.g. Guckenheimer & Holmes Reference Guckenheimer and Holmes1983) ensures that for $0<{\it\varepsilon}\ll 1$ a measurable set of toroidal trajectories are preserved, allowing for the presence of non-mixing islands. The condition for toroidal trajectories, i.e. material barriers to persist is that the ratio ${\it\delta}=T_{SL}/T_{0}$ be sufficiently far from any rational number $p/q$ (with $p$ and $q$ not too large in practice when dealing with computer arithmetics). This predictive argument associates transport barriers in the unsteady regime with certain streamlines of the steady base flow that are more robust to unsteady forcing than others. Streamlines for which ${\it\delta}$ is close to a rational number $p/q$ can lock on a periodic motion with $T_{SL}=pT_{0}/q$ . Such resonant streamlines undergo a saddle–node bifurcation giving rise to the alternation between saddles (hyperbolic orbits) and centres (elliptic orbits) illustrated in figure 11( f). Manifolds originating from the saddle points form homoclinic and heteroclinic tangles enclosing the elliptic points and are responsible for the emergence of a stochastic layer in which mixing occurs locally. Similarly, other resonances can occur inside each region surrounded by a heteroclinic tangle, resulting in a complex tapestry of barriers and stochastic layers. The width of these stochastic layers is usually smaller with increasing values of $q$ but grows with ${\it\varepsilon}$ . Several neighbouring layers predicted by this arithmetic criterion can hence overlap as ${\it\varepsilon}$ increases, hence producing a wider mixing zone in which diffusive motion of the tracers is possible (Cencini, Cecconi & Vulpiani Reference Cencini, Cecconi and Vulpiani2010; Mitchell & Grigoriev Reference Mitchell and Grigoriev2012).

Resonant streamlines, via the formation of heteroclinic tangles, are naturally associated with periodic points in ${\it\Pi}$ , be they hyperbolic or elliptic. Hyperbolic and elliptic points appear together via saddle–node bifurcations and hence come in pairs. Conversely, once a periodic point is identified in ${\it\Pi}$ , it is not trivial to distinguish whether it is originally associated with the resonance of a given streamline or with one of the tangles predicted by Melnikov analysis. The KAM theorem suggests that non-resonant tori should be present in the vicinity of most elliptic periodic points, implying non-mixing islands (Wiggins & Ottino Reference Wiggins and Ottino2004) and thus non-hyperbolic statistics. The question remains however non-trivial in the case of hyperbolic points. The distinction is further complicated by the fact that periodic points can undergo bifurcations as ${\it\varepsilon}$ is increased. The simplest approach to this question is to check, orbit by orbit, whether the resonance criterion

is verified or not. Considering again the segment ${\bf\chi}$ defined in § 3.2, we plot in figure 12 the period of all periodic orbits crossing ${\bf\chi}$ versus the crossing coordinate $X_{0}$ , together with the return times $T_{SL}$ associated with the steady base flow. Saddle points are shown as points, elliptic points as open circles. The periodic points whose period falls close to integer multiples of $T_{SL}(X_{0})$ should correspond to the thickest resonant tori with $q=1$ . The value of $\mathit{Ra}=1.625\times 10^{8}$ is chosen for comparison with the statistical analysis of § 3.3.

Figure 12.

Plot of the period of all periodic orbits crossing ${\bf\chi}$ versus the crossing coordinate $X_{0}$ , together with the return times $T_{SL}$ associated with the steady base flow. Solid lines: $T_{SL}$ as a function of $X_{0}$ with integer multiples $n\times T_{SL}$ for $\mathit{Ra}=1.625\times 10^{8}$ . Black dots (red online): UPOs crossing the segment ${\bf\chi}=\{X\in [0,0.05],Z={\textstyle \frac{1}{2}}\}$ . Circles (blue online): elliptic periodic orbits. Grey dots: distribution of return times $T_{P}({\it\Phi},X_{0})$ (from figure 5 b). Here $A_{i}$ and $N_{i}$ , $i=1,2,3$ correspond to the areas defined in figure 6.

The distribution of the periods is of course limited by the number of UPOs found. However even with 168 orbits this distribution reveals interesting features. A vast majority of the periodic orbits, both hyperbolic and elliptic, whose period is commensurate with $T_{SL}$ in figure 12, correspond to values of ${\it\delta}=T_{SL}/T_{0}=1$ to 5 and belong to the region denoted $N_{2}$ in § 3.3. Other elliptic points were also identified in the region $N_{1}$ , as well as in the region $A_{3}$ with ${\it\delta}=1$ and in $N_{3}$ with ${\it\delta}=2$ and 5. Neither hyperbolic orbits nor elliptic orbits have been detected in $A_{1}$ and $A_{2}$ using this criterion. Note that this approach is necessarily partial since it requires the periodic orbits to cross the segment ${\bf\chi}$ , while the choice of ${\bf\chi}$ is convenient but still arbitrary.

A direct interpretation of figure 12 in connection with the statistics of § 3.3 is now possible. The presence of elliptic orbits satisfying the criterion (4.5) implies the proximity of non-resonant orbits (KAM tori or non-mixing islands). By default, periodic orbits which according to (4.5) do not come from resonances, are most likely embedded in the tangles associated with the destabilisation of homo/heteroclinic loops as in § 4.3. Consequently the region $N_{2}$ is expected to be populated by many non-hyperbolic objects such as the non-mixing islands of figure 11( f). Here $N_{1}$ and $N_{3}$ contain a small but non-zero number of such elliptic orbits and, hence, should contain a few non-mixing islands too. Here $A_{3}$ contains very few low-period elliptic orbits, which explains that its non-hyperbolic statistical properties emerge for large times. Here $A_{2}$ does not show any such resonance, which by default is consistent so far with its hyperbolic nature. In addition, the only periodic points identified in $A_{2}$ are hyperbolic and occur close to the fractal distribution of the return times $T_{P}$ , which is consistent with the picture of a homo/heteroclinic tangle of figure 8. Yet it is fair to recall that the statistics in question are also dependent on the choice of the segment ${\bf\chi}$ , so that non-hyperbolicity is not strictly ruled out provided a few non-mixing regions exist that do not cut ${\bf\chi}$ . This case will be examined in further detail in § 5 using spreading maps. The case of $A_{1}$ , which does not feature resonances according to figure 12 but still shows clear non-hyperbolic statistics, seems contradictory. This case is however peculiar since this region also includes the streamlines closest to the walls. Because of the no-slip condition, the velocity $\boldsymbol{U}$ vanishes at the wall, and hence the period $T_{SL}$ diverges too, making resonances impossible. Here $A_{1}$ is thus expected without any preliminary analysis to display strongly non-hyperbolic statistics since the wall itself, as well as all neighbouring streamlines, can be interpreted as a non-resonant KAM torus.

We pay particular attention to some of the periodic points represented in figure 12 both as centres and saddles, mainly distributed in the region $N_{2}$ . We cannot certify the stability of those periodic points through the finite-difference computation of § 4.2. Because of the fractal configuration of the KAM tori alternating with mixing areas, the perturbed trajectories used for the linearisation can very well be located outside the zone under investigation. They can land on another KAM torus (the point is identified as a centre) or in a locally mixing area (it is detected as a saddle). This ambiguous property is not a numerical artefact but reflects the multiscale structure of ${\it\Pi}$ . Note that this is in any case a minor issue since these points never contribute significantly to the spreading of the tracers.

The fine-scale coexistence of hyperbolic and non-hyperbolic regions is finally illustrated for $\mathit{Ra}=1.625\times 10^{8}$ around the vortex in the middle-right of the cavity, i.e. around $X=Z\approx 0.45$ , see figure 13. This area involves mainly the regions $A_{3}$ and $N_{3}$ as well as the homoclinic loop $R{\it\Gamma}_{2,2}$ . Figure 13(a) shows all periodic orbits identified in the continuous-time representation, along with the homo/heteroclinic loops. Figure 13(b,c) show the spreading maps associated with all those orbits for a short and a long observation time, respectively. Note that the heteroclinic loop ${\it\Gamma}_{2,1}$ acts as a barrier isolating the central region of the cavity from the rest of the flow. It is visually clear from figure 13 that the complete or incomplete character of the mixing, as well as the associated statistics, result from a compromise between the hyperbolic tangle predicted by the Melnikov function ${\it\Lambda}_{2,2}$ and the many KAM tori (barriers) remaining in the flow for this value of $\mathit{Ra}$ .

Figure 13.

(a) Periodic orbits (grey lines) and corresponding periodic points in ${\it\Pi}$ (circles) identified numerically in the vicinity of $R{\it\Gamma}_{2,2}$ for $\mathit{Ra}=1.625\times 10^{8}$ . Black lines: ${\it\Gamma}_{21}$ , $R{\it\Gamma}_{2,2}$ and $R{\it\Gamma}_{2,1}$ . (b) Stable (grey, green online) and unstable (black, red online) manifolds of these periodic points. The main homoclinic tangle originates from $X_{P_{2}}$ and surrounds non-mixing islands and stochastic layers formed by resonant streamlines. (c) Same as (b) with partial spreading maps superimposed. Stochastic layers overlap leading to a complex mixing zone.