Bifurcation analysis of two-dimensional Rayleigh--Bénard convection using deflation
BBifurcation analysis of two-dimensional Rayleigh–B´enard convection using deflation
N. Boull´e, ∗ V. Dallas, † and P. E. Farrell ‡ Mathematical Institute, University of Oxford, Woodstock Road, Oxford OX2 6GG, UK (Dated: March 1, 2021)We perform a bifurcation analysis of the steady state solutions of Rayleigh–B´enard convectionwith no-slip boundary conditions in two dimensions using a numerical method called deflated con-tinuation. By combining this method with an initialisation strategy based on the eigenmodes of theconducting state, we are able to discover multiple solutions to this non-linear problem, includingdisconnected branches of the bifurcation diagram, without the need of any prior knowledge of thedynamics. One of the disconnected branches we find contains a s-shape bifurcation with hysteresis,which is the origin of the flow pattern that may be related to the dynamics of flow reversals in theturbulent regime. Linear stability analysis is also performed to analyse the steady and unsteadyregimes of the solutions in the parameter space and to characterise the type of instabilities.
I. INTRODUCTION
This paper focuses on the bifurcations of the two-dimensional Rayleigh–B´enard convection, which supports a profu-sion of states. There have been numerous investigations trying to visualise solutions and characterise the steady-statesof Rayleigh–B´enard convection, as well as computing its bifurcation diagrams with respect to Ra , the Rayleigh num-ber, in various settings and for different geometries [1, 2]. Ouertatani et al. [3] performed numerical simulations ofRayleigh–B´enard convection in a square cavity using the finite volume method and computed fluid flow profiles andtemperature patterns at Ra = 10 , , . Numerical studies have been conducted to understand the existence ofbifurcating solutions to this problem in a two-dimensional domain with periodic boundary conditions in the horizon-tal direction [4, 5]. Mishra et al. [6] analysed the effect of low Prandtl number on the bifurcation structures, whilePeterson [7] used arclength continuation [8] to study the evolution of cell solutions with respect to the aspect ratioof the domain. Bifurcation structures of Rayleigh–B´enard convection have been extensively studied in cylindricalgeometry by Ma et al. [9] as well as by Boro´nska and Tuckerman [10, 11]. In particular, references [10, 11] adapted atime-dependent code to perform branch continuation and they analysed the linear stability of the computed solutionsusing Arnoldi iterations [12]. Finally, Puigjaner et al. [13, 14] computed bifurcation diagrams for Rayleigh–B´enardconvection in a cubical cavity over different intervals of Ra < . × .A common way to reconstruct bifurcation diagrams is to use arclength continuation and branch switching tech-niques [8, 15, 16] to continue known solutions in a given parameter. These solutions can be computed by applyingNewton’s method to suitably chosen initial guesses or by using time-evolution algorithms. Then, a standard numeri-cal technique in fluid dynamics for finding initial states for the continuation algorithm is to perform time-dependentsimulation of initial solutions until a steady-state is reached [17–19]. These numerical techniques have been widelyused and extremely successful at computing bifurcation diagrams of nonlinear PDEs.In this paper, we employ a recent algorithm called deflation for computing multiple solutions to nonlinear partialdifferential equations [20, 21], which has been successfully applied to compute bifurcations to a wide range of physicalproblems such as the deformation of a hyperelastic beam [21], Carrier’s problem [22], cholesteric liquid crystals [23],and the nonlinear Schr¨odinger equation in two and three dimensions [24–26]. We emphasise that deflated continuationand its extensions offer some advantages that could be combined with the standard bifurcation analysis tools. Themain strength of deflation is that it allows the detection of disconnected branches, in addition to those connected toknown branches. An additional advantage is that it does not require the solution of augmented problems to find newsolutions. Using perturbed solutions to the linearised equations as initial conditions, we are able to discover numeroussteady-states solutions of Rayleigh–B´enard convection. These solutions can be obtained without the need of priorknowledge of the dynamics and regardless of the stability of the solutions.While much progress in fluid dynamics has been made on the study of hydrodynamic instabilities, instabilities thatoccur when a control parameter is varied within the turbulent regime remain poorly understood [27]. The difficultyin studying such transitions arises from the underlying turbulent fluctuations which make analytical approachescumbersome. These types of transitions resemble more closely phase transitions in statistical mechanics becausethe instability occurs on a fluctuating background [28]. In Rayleigh–B´enard convection, such transitions have been ∗ Email: [email protected] † Email: [email protected] ‡ Email: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b observed in the form of reversals of the large scale circulation in an enclosed rectangular geometry, both in experimentsand numerical simulations [29–31]. We therefore investigate the bifurcations of Rayleigh–B´enard convection in anenclosed rectangular geometry to attempt to understand the genesis of the flow patterns that play the central role inflow reversals in the turbulent regime.The paper is organised as follows. We first discuss the problem set-up and the choice of parameters in Section II.Then, in Section III, we briefly review deflated continuation and present the methods used to improve its initialisationprocedure and compute the linear stability of the solutions. Next, we analyse the different branches discovered bydeflation and their stability in Section IV. Finally, in Section V, we recapitulate the main conclusions of this studyand propose potential extensions for future work. II. PROBLEM SET-UP
We consider Rayleigh–B´enard convection [32–34] of a confined fluid heated from below and maintained at a constanttemperature difference ∆ T = T − T across a two dimensional square cell of height d . For simplicity, we employthe Oberbeck-Boussinesq approximation [35–37] in which the kinematic viscosity ν and the thermal diffusivity κ do not depend on temperature while the fluid density ρ is assumed to be constant except in the buoyancy termof the momentum equation. In this case, we assume the density to depend linearly on the temperature, ρ ( T ) (cid:39) ρ ( T )[1 − α ∆ T ], where α is the thermal expansion coefficient. Then, the governing equations of the problem are ∇ · u = 0 , (1a) ∂ t u + u · ∇ u = −∇ p + P r ∇ u + P rRaT ˆ z , (1b) ∂ t T + u · ∇ T = ∇ T, (1c)where u is the velocity field, T is the temperature field, p is the pressure, and ˆ z is the buoyancy direction. The aboveequations have been nondimensionalized using d , d /κ and ∆ T as the relevant scales for length, time and temperature,respectively. The two dimensionless parameters of the problem are the Rayleigh and Prandtl numbers Ra = gα ∆ T d /νκ, P r = ν/κ, (2)where g is the gravitational acceleration. The Rayleigh number represents the ratio of the acceleration gα ∆ T relatedto buoyancy to the stabilizing effects of ν and κ . The bifurcation parameter in the problem we study is Ra by fixing P r = 1. The cell is assumed to have rigid walls, with thermally conducting horizontal walls and thermally insulatingside walls, i.e. u = w = 0 , x = z = 0 , , (3a) T = 1 , z = 0 , (3b) T = 0 , z = 1 , (3c) ∂ x T = 0 , x = 0 , . (3d)The trivial steady state solution is a motionless state with a negative thermal gradient upwards through the layer u = 0 , T = 1 − z, (4)called conducting state. This is because the fluid acts as a conducting material. The cooler fluid near the top of thelayer is denser than the warmer fluid underneath it. The classical stability analysis on this steady state for a criticalvalue of the Rayleigh number Ra c reveals the well-known Rayleigh–B´enard cellular instability which gives a patternof two counter-rotating convection rolls. Note that such a pattern is not realisable in small aspect ratio domains asit will also be shown in Section IV A.The symmetries of the problem are: a) the mirror symmetry with respect to the axis x = 1 / u, w, T ]( x, z ) → [ − u, w, T ](1 − x, z ) , (5)which leaves the velocity and temperature field invariant, b) the mirror symmetry with respect to the axis z = 1 / u, w, T ]( x, z ) → [ u, − w, − T ]( x, − z ) , (6)which leaves the velocity field invariant but transforms the temperature field into its opposite and c) the combinationof these two symmetries, which is essentially a rotation by π with respect to the center of the cell ( x, z ) = (1 / , / III. NUMERICAL METHODSA. Computation of multiple solutions with deflation
Steady-states solutions to Rayleigh–B´enard convection are computed using a recent numerical technique calleddeflation [20]. This algorithm is based on Newton’s method and allows the computation of multiple solutions to anonlinear system of equations F ( φ, λ ) = 0, where φ is the solution and λ ∈ R is a bifurcation parameter. G ( φ, λ ) := (cid:18) (cid:107) φ − φ (cid:107) + 1 (cid:19) F ( φ, λ ) , (7)such that G ( φ, λ ) does not converge to zero as φ → φ . The term M ( φ, φ ) := (cid:107) φ − φ (cid:107) − + 1 is called a deflationoperator and ensures that Newton’s method applied again to G will not converge to a previously computed solution.By adding the unity in the expression of the deflation operator, we impose that the new problem G behaves similarlyto F away from the root φ . This process can be repeated to obtain a set of solutions S ( λ ) to a nonlinear problem F ( φ, λ ) = 0 at the bifurcation parameter λ . The next step is the continuation of the previously discovered solutionsat a new bifurcation parameter λ − ∆ λ : each φ ∈ S ( λ ) is used as initial guess for the deflation procedure in order toconstruct a new set of solution S ( λ − ∆ λ ) at λ − ∆ λ .We apply this algorithm repeatedly to discover multiple solutions to the steady state Rayleigh–B´enard problem: F ( φ, Ra ) := − u · ∇ u − ∇ p + P r ∇ u + P rRaT ˆ z = 0 , ∇ · u = 0 , − u · ∇ T + ∇ T = 0 , (8)where φ = ( u , p, T ) and Ra is the bifurcation parameter. This equation is discretised with Taylor–Hood finiteelements for the velocity and pressure (piecewise biquadratic and piecewise bilinear respectively), together withpiecewise bilinear polynomials for the temperature, using the Firedrake finite element library [38]. The unit squaredomain is represented by a mesh with 50 ×
50 square cells to preserve the symmetries of the problem and avoidnumerical artefacts. A coarse discretization is needed due to the complexity of the problem and number of distinctsolutions discovered, leading to hundred of thousands of deflation steps and Newton iterations. Moreover, we have notobserved significant differences in selected experiments where we refined the mesh. A crucial question in the deflationtechnique is the choice of the deflation operator M because it can affect the convergence of Newton’s method. Weuse the following operator, M (( u , p, T ) , ( u , p , T )) = (cid:18) (cid:107) u − u (cid:107) + (cid:107)∇ ( u − u ) (cid:107) + (cid:107) T − T (cid:107) + 1 (cid:19) , (9)where (cid:107) · (cid:107) denotes the L -norm. It is important to note that the deflation operator M does not depend on thepressure p because we wish to deflate all the solutions of the form ( u , p + c, T ), where c ∈ R , to avoid discoveringsolutions trivially related to known ones.We perform bifurcation analysis of Rayleigh–B´enard convection and compute solutions to the nonlinear system ofEq. (8) for 0 ≤ Ra ≤ using deflated continuation [20, 21]. This method combines deflation with continuation ofthe solutions in the bifurcation parameter in order to continue and discover new branches of solutions as we decreasethe bifurcation parameter from Ra = 10 , − ∆ R a , . . . ,
0. Here, ∆ Ra denotes the continuation step size in thebifurcation parameter, which is chosen to be equal to ∆ Ra = 100. In a previous work on the computation of solutionsto the 3D nonlinear Schr¨odinger equation [26], we found that having good initial guesses (e.g. solutions of the linearisedproblem) for the initial deflation steps facilitates the convergence of Newton’s method and leads to more complex andinteresting states. However, contrary to [26] we do not have the analytical expressions of the linear solutions to theRayleigh–B´enard convection due to the no-slip boundary conditions. Therefore, we numerically solve an eigenvalueproblem (see Section III B) to obtain the first ten unstable eigenmodes ( u , T ) , . . . , ( u , T ) linearised around theconducting steady state Eq. (4) (see Fig. 2). Then, the sums of the conducting state with the respective normalizedperturbations: ( u , − z + T ) , . . . , ( u , − z + T ) are used as initial guesses for deflated continuation. B. Linear stability analysis
The stability analysis of a given steady state ( u , T ) is performed by considering the perturbation ansatz u ( x, z, t ) = u + v e λt , (10a) T ( x, z, t ) = T + θe λt , (10b)where the velocity perturbation v (cid:28) θ (cid:28) λ .When at least one of the eigenvalues has a positive real part R ( λ ) >
0, we can have two types of instabilities. Thestationary instability which occurs when the imaginary part of the eigenvalue I ( λ ) = 0 and the oscillatory instabilitywhich occurs when I ( λ ) (cid:54) = 0. If all the eigenvalues have negative real parts, then the perturbed steady state is stable.Inserting Eq. (10) into Eq. (1) and considering the linearised system of equations, we obtain the following generalisedeigenvalue problem at leading order A −∇ RaP r ˆ z ∇· −∇ T · ∇ − u · ∇ v pθ = λ I I v pθ , (11)where I is the relevant identity operator and A is the linear operator defined by A v = ∇ v − u · ∇ v − v · ∇ u . The eigenvalue problem described in Eq. (11) is solved with a Krylov–Schur method [39] using the SLEPc library [40].To obtain the critical value of the Rayleigh number Ra c at which the conducting state Eq. (4) becomes unstable,we modify the eigenvalue problem in Eq. (11). Inserting Eq. (4) into Eq. (11) we obtain the following generalizedeigenvalue problem ∇ −∇ RaP r ˆ z ∇· z · ∇ v pθ = λ I I v pθ . (12)A steady state becomes unstable when λ = 0 and this corresponds to the critical Rayleigh number Ra c which is thesmallest Ra satisfying Eq. (12) for λ = 0. So, to obtain the critical value of the Rayleigh number for different basestates we have reformulated the problem to a generalised eigenvalue problem for Ra c as follows ∇ −∇ ∇· z · ∇ v pθ = Ra c − P r ˆ z v pθ . (13)In the following sections this linear analysis will be used to study the stability of the conducting state Eq. (4) and ofthe nonlinear steady state solutions that we obtain from deflated continuation. IV. RESULTSA. Primary instabilities
In this section, we analyse the stability of the conducting steady state defined by Eq. (4). We found that the firstinstability of this state arises at Ra (1) c := Ra c ≈ Ra ∗ c ≈ .
762 for an unbounded domain with periodic streamwise boundariesand no-slip spanwise boundaries. This difference in the value of Ra c is due to the effect of the side walls on theinstability. This effect should decrease when x becomes large such that the behaviour predicted for an unboundeddomain is recovered. To demonstrate this, we systematically increase the aspect ratioΓ = Ld , (14)by fixing the height d = 1 and increasing the length L of the domain. Figure 1 (left) clearly shows that for Γ (cid:29) Ra (1) c converges to Ra ∗ c . Note, however, that this is not the case for the spatially extended system with stress-free boundary conditions in the limit of Γ (cid:29) d and decreasing the length L of the domain.For aspect ratios Γ < Ra c ∝ Γ − . indicating a more stable flow as Γ decreases. The critical Rayleighnumber for Rayleigh-B´enard convection with stress-free boundary conditions is Ra ∗ c ( k ) = ( π + k ) /k [43], where k is the horizontal wavenumber. For k (cid:29) Ra ∗ c ( k ) ∝ k ∝ L − ∝ Γ − usingEq. (14). This scaling relation is very close to the observed one, which implies that the exponent is not affected muchby the no-slip boundary conditions. The exponent we find is in agreement with [44]. − − − − Γ − . Aspect ratio Γ ( R a c − R a ∗ c ) / R a ∗ c . . . . Ra · R ( λ ) FIG. 1. (Color online) Left: Critical Rayleigh number with respect to the aspect ratio Γ (blue dots). The red dot indicates theaspect ratio Γ = 1, which is the focus of our study. Right: Real part of the growth rate of the eigenmodes that emanate fromthe conducting steady state.TABLE I. The ten primary bifurcations from the conductive state within the range 0 ≤ Ra ≤ . Here n enumerates theeigenmodes by order of appearance as Ra increases. n Ra ( n ) c k x , k z ) (2,1) (3,1) (4,1) (3,2) (2,2) (4,2) (5,1) (5,2) (3,3) (4,3) The case on which we focus for the rest of the paper is the square cell (Γ = 1), which is denoted with a red dot inFig. 1 (left). We first study the linear stability of the conducting state Eq. (4). In the range 0 ≤ Ra ≤ we observe10 supercritical stationary bifurcations that arise from the conducting state. Figure 1 (right) shows the real part ofthe growth rate of the unstable eigenmodes as a function of Ra . The onset of these bifurcations occur when R ( λ ) = 0.The critical values of the Rayleigh number at the onsets are listed in Table I. Figure 1 (right) demonstrates thatthe growth rates of the eigenmodes vary with Rayleigh number and they have different dependence. The eigenmodewhich bifurcated from the conducting state at Ra (2) c = 6746 ends up being the most unstable as Ra → .The velocity magnitude and temperature fields of the eigenmodes from these primary bifurcations of the flow areshown in Fig. 2 with the corresponding critical Rayleigh numbers indicated on the top of each flow pattern. In TableI we characterise the flow patterns by their Fourier modes ( k x , k z ), where k x is the streamwise wavenumber and k z is the spanwise wavenumber of the Fourier modes of the temperature field. Note that due to the small aspect ratioof the domain the flow pattern of the first bifurcation at Ra (1) c = 2586 is just a single convection roll (see Fig. 2)instead of two counter-rotating rolls as it is known from the classical result. In the following sections, we first analysethe branches of solutions that are created from some of these primary bifurcations before moving to disconnectedbranches. B. Bifurcation diagrams
The combination of deflated continuation and the initialisation strategy presented in Section III A identified 129solutions to the two-dimensional steady-state Rayleigh–B´enard convection at Ra = 10 including 43 unique solutionswith respect to the symmetries. Here, we analyse the evolution and the linear stability of some of the branches arisingfrom the first four excited states (see Fig. 2) as well as two disconnected branches. Moreover, we discard the mirrorsymmetric solutions with respect to the x = 1 / z = 1 / Ra , Ra (1) c = 2586 Ra (2) c = 6746 Ra (3) c = 19655 Ra (4) c = 23346 Ra (5) c = 25780 Ra (6) c = 41772 Ra (7) c = 47431 Ra (8) c = 74878 Ra (9) c = 86313 Ra (10) c = 94543FIG. 2. (Color online) The eigenmodes of the primary bifurcations emanating from the conducting steady state solution Eq.(4) in the Rayleigh number range 0 ≤ Ra ≤ . The top rows show the magnitude of the velocity v (red color indicates a highmagnitude and blue a zero velocity profile), while the bottom rows show the temperature θ (red and blue respectively displaypositive and negative temperature fields). are the kinetic energy (cid:107) u (cid:107) , the potential energy (cid:107) T (cid:107) , and the Nusselt number N u , which are defined as (cid:107) u (cid:107) = (cid:90) u d x, (cid:107) T (cid:107) = (cid:90) T d x, N u = (cid:90) |∇ T | d x. (15)The Nusselt number characterises the efficiency of convective heat transfer and is given by the ratio of the total heattransfer (i.e. both advective and diffusive) to the conductive heat transfer N u ≡ HL/κd = (cid:104) wT (cid:105) − κ∂ z (cid:104) T (cid:105) = (cid:104) ( ∇ T ) (cid:105) ,where H is the heat flux and (cid:104) . (cid:105) stands for the area average. A Nusselt number of one represents heat transfer bypure conduction. Here, the Nusselt number is essentially defined as the dissipation rate of the temperature variance(see second equality). This is obtained by integrating the equation for potential energy over the entire area [45].The bifurcation diagrams are presented in Fig. 3, where the numbers denote the different branches of steadysolutions, highlighting the values of kinetic energy, potential energy and Nusselt number that correspond to eachbranch on each diagram.The first steady state solutions that we obtain from deflated continuation as Ra is varied are four branches thatarise from the first four bifurcations of the conducting steady state: branch (1) at Ra (1) c = 2586, branch (2) at Ra (2) c = 6746, branch (3) at Ra (3) c = 19655 and branch (4) at Ra (4) c = 23346. Moreover a secondary branch (5)bifurcates at Ra ≈ Ra ≈ Ra ≈ ≤ Ra ≤ branch (1) has the highest kinetic energy while branch(6) has the highest potential energy (see Fig. 3(a) and 3(b)). The steady state solutions that are most effective inconvecting heat transfer for 0 ≤ Ra (cid:46) × are in branch (1) from the first instability of the conducting state whilefor 3 × (cid:46) Ra ≤ . × are in the branch (2), and in branch (7) for Ra (cid:38) . × (see Fig. 3(c)). . . . . Ra · , , ,
000 (1)(2) (5) (3) (4)(6) (9) k u k . . . . Ra · . . . . Ra · . . . . k T k . . . . Ra · N u (a) (b)(c)FIG. 3. (Color online) Bifurcation diagram of the kinetic energy (cid:107) u (cid:107) of the branches studied in the paper (a). The panels(b) and (c) represent the bifurcation diagrams of the potential energy (cid:107) T (cid:107) and the Nusselt number Nu , respectively. Thebranches that originate from the linear instabilities are highlighted in dashed lines while bifurcations and disconnected branchesare depicted with solid lines. C. Bifurcations from the conducting steady state
In this section we analyse the states that emanate from the conducting steady state and their evolution as afunction of the bifurcation parameter Ra . Our analysis involves the detailed evolution of the velocity magnitude andthe temperature fields on the branches of the bifurcation diagrams of kinetic and potential energy, respectively. Inaddition, we present the real and imaginary parts of the largest growth rates from the linear stability analysis we haveperformed on the steady state solutions.Figure 4 shows results of the aforementioned analysis for branch (1), which arises at Ra (1) c = 2586 from the firsteigenmode (see Fig. 2). This bifurcation breaks only the rotational symmetry by π that is satisfied by the eigenmode.As Ra increases, the magnitude of the velocity field evolves from a circular to a bent convection roll at Ra ≈ × (see Fig. 4(a)). On the other hand, the z-shaped interface of the temperature field is sheared (see Fig. 4(b)), enhancingthe heat transfer in the convection cell (see also branch (1) in Fig. 3(c)). Figure 4(c) demonstrates that the steadystate solution of this branch remains stable ( R ( λ ) <
0) under perturbations with a constant growth rate over the wholerange of Rayleigh numbers we consider. The imaginary part of the largest growth rate is zero from the stationarybifurcation at Ra (1) c = 2586 to the Ra = 10 (see Fig. 4(d)), while the subsequent growth rate has a strong oscillatorycomponent and its real part is R ( λ ) = const < Ra (2) c = 6746, which is depicted in Fig. 5(a)-(b). While the symmetry of the eigenmode with respect to the x -axis is rapidly broken, the velocity and temperatureprofiles remain symmetric with respect to the z -axis in the range of parameters that we consider. Then, the largestgrowth rate in Fig. 5(c) indicates that the perturbed steady state solutions of this branch is stable over the interval Ra ∈ [13500 , Ra ≈ R ( λ ) = 0, which . . . . Ra · , , , k u k . . . . Ra · − − R ( λ ) 0 0 . . . . Ra · . . . . . k T k . . . . Ra · I ( λ ) (a) (b)(c) (d)FIG. 4. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofbranch (1), which arises from the first eigenmode. The color scheme for the temperature field ranges from 0 (blue) to 1 (red).The real and positive imaginary parts of the largest growth rates are presented in (c) and (b), respectively. is analysed later and represented in Fig. 6. Another bifurcation is observed at Ra ≈ Ra increases. This symmetry breaking leads to a primary large scale circulationspanning the domain and a secondary vortex with smaller amplitude in one of the corners depending on whichsymmetric solution one is referring to. As Ra → , the difference in the flow pattern of the velocity magnitudebetween branch (1) and (5) is the secondary vortex in the corner of the convection cell. The symmetry breaking is alsoobvious in the temperature field (see Fig. 6(b)) with the potential energy of the flow to decrease slightly in contrastto the increase of the kinetic energy and the convective heat transfer (see also Fig. 3(c)). The largest growth rateshows that the states are unstable to perturbations without any oscillations (see Fig. 6(d)). The subsequent growthrates shown in Fig. 6(c) are negative and almost constant in this range of Rayleigh numbers.Branch (3) originates from the eigenmode of the conducting state at Ra (3) c = 19655, which is depicted in Fig. 7.Similarly to branch (5), the three symmetries (see Section II) satisfied by the third eigenmode are rapidly broken andthe pattern of the velocity magnitude becomes more complex as Ra increases. This complex pattern is one reason forthe meandering of the potential energy in this range of Rayleigh numbers and enhances the convective heat transfer(see Fig. 3(c)) by mixing the temperature field. Figure 7(c) shows that the steady-state solutions of branch (3) arestationary unstable under the effect of perturbations until Ra ≈ I ( λ ) = 0 or I ( λ ) (cid:54) = 0.Figure 8 shows bifurcation diagrams and the linear stability results of branch (4) which was originated from thefourth linear eigenmode at Ra (4) c = 23346. The flow pattern of the velocity magnitude has a symmetric form ofan array of four vortices, which is almost conserved over the whole range of the Rayleigh numbers we consider asFig. 8(a) illustrates. Similarly, the symmetric pattern of the temperature field remains almost unaffected in the regime23346 ≤ Ra ≤ (see Fig. 8(b)). Actually, the bifurcation breaks symmetry (6) and the rotational symmetry by π of the eigenmode. This is visible within this range of Ra by inspecting the flow patterns displayed on the diagram ofpotential energy.Figure 8(c) shows the real parts of the four largest growth rates as a function of Ra . At Rayleigh numbers close to Ra (4) c = 23346 we see that the three eigenmodes are stationary unstable, while the fourth one is stationary stable. As Ra increases the first with the third most unstable eigenmode (coloured in red and green, respectively, in Fig. 8(c)) . . . . Ra · , , , , k u k . . . . Ra · − − R ( λ ) 0 0 . . . . Ra · . . . k T k . . . . Ra · I ( λ ) (a) (b)(c) (d)FIG. 5. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofbranch (2), which arises from the second eigenmode. The real and positive imaginary parts of the largest growth rates arepresented in (c) and (b), respectively. The two largest growth rates are shown in red and green while the subsequent ones aredepicted in blue. . . . . Ra · , , , , k u k . . . . Ra · − − R ( λ ) 0 0 . . . . Ra · . . . . . . k T k . . . . Ra · I ( λ ) (a) (b)(c) (d)FIG. 6. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofbranch (5), which bifurcates from branch (2) at Ra ≈ coalesce at Ra ≈ . . . . Ra · , , k u k . . . . Ra · R ( λ ) 0 0 . . . . Ra · . . . . k T k . . . . Ra · I ( λ ) (a) (b)(c) (d)FIG. 7. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energyof branch (3), which arises from the third eigenmode. The real and positive imaginary parts of the largest growth rates arepresented in (c) and (b), respectively. . . . . Ra · , k u k . . . . Ra · − R ( λ ) 0 . . . . Ra · . . . . k T k . . . . Ra · I ( λ ) (a) (b)(c) (d)FIG. 8. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofbranch (4), which arises from the fourth eigenmode. The real and positive imaginary parts of the largest growth rates arepresented in (c) and (b), respectively. as Ra → . Similarly, the second and fourth most unstable mode coalesce at Ra ≈ < Ra ≤ .We close our discussion on the branches arising from the conducting steady state solution with branch (6), presentedin Fig. 9. This branch bifurcates from branch (4) at Ra ≈ . . . . Ra · , k u k . . . . Ra · R ( λ ) 0 . . . . Ra · . . . . k T k . . . . Ra · I ( λ ) (a) (b)(c) (d)FIG. 9. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofbranch (6), which bifurcates from branch (4) at Ra ≈ by π of the flow pattern with the four vortices. Close to Ra ≈ Ra ≈ Ra → . The linear stability analysis of this branch shows that the steadystates are very sensitive to perturbations and thus stationary unstable throughout the range of Ra considered in thisstudy. D. Disconnected branches
We now focus on the disconnected branch (7) (see Fig. 3), which comprises an upper and a lower branch, depictedin blue and red in Fig. 10, respectively. At the upper (blue) branch the flow pattern of the velocity magnitude is athree vortex state with one large titled vortex spanning the domain and two smaller vortices located at the corners. As Ra → we notice that the state looks similar to the branches (1) and (5). The difference is that branch (5) involvesonly one smaller vortex at the corner of the domain and branch (1) has no vortices at the corners, while branch (7)has two smaller vortices located at the corners. At Ra ≈ R ( λ ) increases and seems to saturate as Ra → . The imaginary part of the growth rateis zero in the whole range of Ra numbers suggesting that the perturbed steady state is stationary unstable in thisregime. The potential energy of the upper (blue) branch decreases with Ra in contrast to the kinetic energy. As Ra increases the flow pattern of the potential energy becomes more efficient at transferring heat due to the large scalecirculation in the flow (see Fig. 3(c)).The upper (blue) branch plays an important role on the nonlinear dynamics of Rayleigh-B´enard convection atvalues of Ra > . The flow pattern of the kinetic energy shown for Ra = 10 (see Fig. 10(a)), with the large scalevortex and the two smaller vortices at the corners of the domain, is reminiscent to the flow patterns from studieson the dynamics of random-in-time flow reversals [29–31, 47]. These reversals in the flow direction of the large scale2 . . . . Ra · , , , k u k . . . . Ra · − − R ( λ ) 0 . . . . Ra · − − R ( λ ) 0 . . . . Ra · . . . . . k T k . . . . Ra · I ( λ ) 0 . . . . Ra · I ( λ ) (a) (b)(c) (d)(e) (f)FIG. 10. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofthe disconnected branch (7). The upper and lower parts of the branch are coloured in blue and red, respectively. The real andpositive imaginary parts of the largest growth rates of the upper (resp. lower) part of the branch are presented in (c) and (b)(resp. (e) and (f)), respectively. circulation at irregular time intervals have been observed in domains with free-slip or no-slip boundaries and in anarrow band of Rayleigh numbers in the range 10 ≤ Ra ≤ depending on the Prandtl number [29–31]. In thisregime the system fluctuates in time between essentially the state we show for Ra = 10 (see Fig. 10(a)) and itsmirror symmetric state. Therefore, for values of the bifurcation parameter Ra > we expect the perturbed steadystates of this branch to gradually become stable as they turn out to be the statistically steady state solution in therange 10 ≤ Ra ≤ . So, as Ra increases from 10 to 10 the system transitions from one steady state solution toanother. This big change in the dynamics of the flow are caused by bifurcations over a turbulent background [27].Such bifurcations have been observed in many other types of flows [48–52] and understanding their dynamics is anopen question of great interest.At the lower (red) branch illustrated in Fig. 10(a) we observe the central vortex to be oriented more vertically,giving space to the smaller vortices at the corners to extend along the z axis and almost reach the height of thedomain. This flow pattern is an extension of the two-vortex pattern observed in branch (5) (see Fig. 6) and it hassimilar features to the three-vortex pattern found near the linear limit of branch (3) (see Fig. 7). The patterns ofthe kinetic and potential energy of this branch are almost invariant within the range 30000 (cid:46) Ra ≤ . The flowpattern of the potential energy clearly manages to effectively push the hot fluid toward the cold plate at the top andvice versa. This efficiency in convective heat transfer is clearly depicted in Fig. 3(c), where this part of branch (7)has the highest Nusselt number for Ra > R ( λ ) ≈ Ra = 10 (see Fig. 10(e)). On the otherhand, the imaginary part of the growth rate I ( λ ) = 0 at the threshold of the disconnected branch and as Ra increasesthe instability becomes progressively oscillatory (see Fig. 10(f)). The subsequent growth rates presented in Fig. 10(e)3remain negative as Ra increases suggesting that this part of branch (7) is another preferred steady state solution ofthe system particularly for Ra > .Branch (8) illustrated in Fig. 11 bifurcates from the lower part of the disconnected branch (7) at Ra ≈ . . . . . . . Ra · , , , k u k . . . . . . . Ra · − R ( λ ) 0 . . . . . . . Ra · . . . . . k T k . . . . . . . Ra · I ( λ ) (a) (b)(c) (d)FIG. 11. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofbranch (8), which bifurcates from the lower part of the disconnected branch (7) at Ra ≈ elongated vortex on its left and a smaller vortex at the bottom right corner (see Fig. 11(a)). From Fig. 3 we can seethat even though the kinetic energy of branch (8) is equal to that of the upper part of branch (7) for Ra > Ra increases (see Fig. 3(c)).However, it is not as efficient in convecting heat as the pattern of branch (7) (see lower branch in 10(b)) because twohot plumes are obviously more effective than one. The steady state solution of branch (8) is stationary unstable tofinite perturbations over the range 37000 (cid:46) Ra ≤ except for a small range around Ra = 50000 where it becomesoscillatory unstable. This can be seen from Fig. 10(c) and 10(d), where the two largest growth rates (coloured in redand blue) swap with each other in this small range of Rayleigh numbers.The second disconnected branch that we find in this range of Rayleigh numbers is branch (9) and arises at Ra ≈ R ( λ ) > I ( λ ) = 0, with very rangegrowth rates. This suggests that this branch is very unlikely to be observed in a laboratory experiment as it is verysensitive to finite perturbations and highlight the effectiveness of deflated continuation in identifying such unstabledisconnected branches and steady-state solutions.4 . .
75 0 . .
85 0 . .
95 1 Ra · k u k . .
75 0 . .
85 0 . .
95 1 Ra · R ( λ ) 0 . .
75 0 . .
85 0 . .
95 1 Ra · R ( λ ) 0 . .
75 0 . .
85 0 . .
95 1 Ra · . . . k T k . .
75 0 . .
85 0 . .
95 1 Ra · I ( λ ) 0 . .
75 0 . .
85 0 . .
95 1 Ra · I ( λ ) (a) (b)(c) (d)(e) (f)FIG. 12. (Color online) Evolution of the steady state flow patterns for (a) the kinetic energy and (b) the potential energy ofthe disconnected branch (9). The upper and lower parts of the branch are coloured in blue and red, respectively. The real andpositive imaginary parts of the largest growth rates of the upper (resp. lower) part of the branch are presented in (c) and (b)(resp. (e) and (f)), respectively. V. CONCLUSIONS
In this work, we have computed and analysed the linear stability of steady-state solutions of the two-dimensionalRayleigh–B´enard convection in a square cell with rigid walls. The combination of deflated continuation and a suitableinitialization procedure based on the linear stability of the motionless state revealed to be a powerful numericaltechnique for finding multiples solutions including disconnected branches, which are not easily accessible by standardbifurcation techniques such as arclength continuation and time-evolution algorithms. A highlight is the success ofdeflation to discover these branches in an automated manner without any prior knowledge of the dynamics.We classify the solutions based on the kinetic and potential energy as well as the Nusselt number. Our classificationprovides a clear view of the nonlinear dynamics of Rayleigh–B´enard convection up to a Rayleigh number of 10 . Weshow that the steady state solutions of branch (1) (which originates from the conducting state) dictate the dynamicsof the flow for Ra < < Ra ≤ the flow patterns of the disconnected branch (7) start to playa very important role. This is also indicated by the linear stability of solutions as Ra varies. Solutions on branch (7)are reminiscent of the flow pattern that plays the fundamental role in the turbulent dynamics of flow reversals thathave been reported in the literature for Rayleigh–B´enard convection at Ra ∼ [29–31]. This disconnected branchexhibits a s-shape bifurcation with hysteresis which can be a challenge to discover with other bifurcation techniques.There are several possible extensions to this work. An interesting but computationally challenging future study isto extend certain solutions of interest to higher values of Ra and determine whether they give the required theoreticalscalings in the limit of high Rayleigh number. A first attempt toward this direction has been done in Rayleigh–B´enardconvection with periodic boundary conditions in the horizontal direction and no-slip in the vertical direction [53, 54].Moreover, the extensions of branch (7) and its stability in the regime where reversals occur is also of great interest.5Another extension of interest is to study how the states presented in this study change with respect to otherparameters of the system such as the Prandtl number or the aspect ratio of the domain Γ. As observed in the leftpanel of Fig. 1, the critical Rayleigh number converges to Ra ∗ c ≈ .
762 for large aspect ratio. Using alternativeformulations of the Rayleigh–B´enard problem, with the aspect ratio playing the role of the bifurcation parameter,it is possible to analyse the evolution of the branches for large aspect ratio and connect them with the solutionsof the periodic set-up, which has been studied well in the literature, both analytically [41] and numerically [4, 5].This connection will be of great interest to theoreticians and hopefully will spark their interest to explore analyticalsolutions on more realistic convection set-ups with spanwise walls.A different future direction would be to perform bifurcation analysis on Rayleigh–B´enard convection in a three-dimensional rectangular domain. Such an approach will be of interest to experimentalists since direct comparisonswith experiments would be possible. However, while the generalisation of the bifurcation technique employed inthis paper to three dimensions is straightforward, high Rayleigh numbers require an efficient solver to perform theunderlying Newton iterations. One possible solution would be to construct an efficient preconditioner robust withrespect to the Rayleigh number to be able to reach high Ra regimes in the spirit of [55, 56]. ACKNOWLEDGMENTS
We would like to thank P. A. Gazca-Orozco and S. Fauve for helpful discussions that improved the quality of themanuscript. This work is supported by the EPSRC Centre For Doctoral Training in Industrially Focused MathematicalModelling (EP/L015803/1) in collaboration with Simula Research Laboratory, and by EPSRC grants EP/V001493/1and EP/R029423/1. [1] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. , 851 (1993).[2] E. Bodenschatz, W. Pesch, and G. Ahlers, Recent developments in Rayleigh-B´enard convection, Annu. Rev. Fluid Mech. , 709 (2000).[3] N. Ouertatani, N. B. Cheikh, B. B. Beya, and T. Lili, Numerical simulation of two-dimensional Rayleigh–B´enard convectionin an enclosure, C. R. Mec. , 464 (2008).[4] E. Zienicke, N. Seehafer, and F. Feudel, Bifurcations in two-dimensional Rayleigh-B´enard convection, Phys. Rev. E ,428 (1998).[5] S. Paul, M. K. Verma, P. Wahi, S. K. Reddy, and K. Kumar, Bifurcation analysis of the flow patterns in two-dimensionalRayleigh–B´enard convection, Int. J. Bifurcat. Chaos , 1230018 (2012).[6] P. K. Mishra, P. Wahi, and M. K. Verma, Patterns and bifurcations in low–Prandtl-number Rayleigh-B´enard convection,EPL , 44003 (2010).[7] J. W. Peterson, Parallel adaptive finite element methods for problems in natural convection , Ph.D. thesis, University ofTexas, Austin (2008).[8] H. B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems., in
Applications of Bifurcation Theory (Academic Press, 1977) pp. 359–384.[9] D.-J. Ma, D.-J. Sun, and X.-Y. Yin, Multiplicity of steady states in cylindrical Rayleigh-B´enard convection, Phys. Rev. E , 037302 (2006).[10] K. Boro´nska and L. S. Tuckerman, Extreme multiplicity in cylindrical Rayleigh-B´enard convection. I. Time dependenceand oscillations, Phys. Rev. E , 036320 (2010).[11] K. Boro´nska and L. S. Tuckerman, Extreme multiplicity in cylindrical Rayleigh-B´enard convection. II. Bifurcation diagramand symmetry classification, Phys. Rev. E , 036321 (2010).[12] W. E. Arnoldi, The principle of minimized iterations in the solution of the matrix eigenvalue problem, Q. Appl. Math. ,17 (1951).[13] D. Puigjaner, J. Herrero, F. Giralt, and C. Sim´o, Stability analysis of the flow in a cubical cavity heated from below, Phys.Fluids , 3639 (2004).[14] D. Puigjaner, J. Herrero, F. Giralt, and C. Sim´o, Bifurcation analysis of multiple steady flow patterns for Rayleigh-B´enardconvection in a cubical cavity at P r = 130, Phys. Rev. E , 046304 (2006).[15] E. J. Doedel, AUTO: A program for the automatic bifurcation analysis of autonomous systems, Congr. Numer , 25(1981).[16] H. Uecker, D. Wetzel, and J. D. M. Rademacher, pde2path-A Matlab package for continuation and bifurcation in 2Delliptic systems, Numer. Math.-Theory Me. , 58 (2014).[17] H. A. Dijkstra, F. W. Wubs, A. K. Cliffe, E. Doedel, I. F. Dragomirescu, B. Eckhardt, A. Y. Gelfgat, A. L. Hazel,V. Lucarini, A. G. Salinger, et al. , Numerical Bifurcation Methods and their Application to Fluid Dynamics: Analysisbeyond Simulation, Commun. Comput. Phys. , 1 (2014). [18] L. S. Tuckerman and D. Barkley, Bifurcation analysis for timesteppers, in Numerical methods for bifurcation problems andlarge-scale dynamical systems (Springer, 2000) pp. 453–466.[19] C. K. Mamun and L. S. Tuckerman, Asymmetry and Hopf bifurcation in spherical Couette flow, Phys. Fluids , 80 (1995).[20] P. E. Farrell, A. Birkisson, and S. W. Funke, Deflation techniques for finding distinct solutions of nonlinear partial differ-ential equations, SIAM J. Sci. Comput. , A2026 (2015).[21] P. E. Farrell, C. H. Beentjes, and ´A. Birkisson, The computation of disconnected bifurcation diagrams, arXiv preprintarXiv:1603.00809 (2016).[22] S. J. Chapman and P. E. Farrell, Analysis of Carrier’s problem, SIAM J. Appl. Math. , 924 (2017).[23] D. B. Emerson, P. E. Farrell, J. H. Adler, S. P. MacLachlan, and T. J. Atherton, Computing equilibrium states of cholestericliquid crystals in elliptical channels with deflation algorithms, Liquid Crystals , 341 (2018).[24] E. G. Charalampidis, P. G. Kevrekidis, and P. E. Farrell, Computing stationary solutions of the two-dimensional Gross–Pitaevskii equation with deflated continuation, Commun. Nonlinear Sci. Numer. Simulat. , 482 (2018).[25] E. G. Charalampidis, N. Boull´e, P. G. Kevrekidis, and P. E. Farrell, Bifurcation analysis of stationary solutions of two-dimensional coupled Gross–Pitaevskii equations using deflated continuation, Commun. Nonlinear Sci. Numer. Simulat. ,105255 (2020).[26] N. Boull´e, E. G. Charalampidis, P. E. Farrell, and P. G. Kevrekidis, Deflation-based identification of nonlinear excitationsof the three-dimensional Gross-Pitaevskii equation, Phys. Rev. A , 053307 (2020).[27] S. Fauve, J. Herault, G. Michel, and F. P´etr´elis, Instabilities on a turbulent background, J. Stat. Mech. Theory Exp. ,064001 (2017).[28] L. P. Kadanoff, Statistical physics: statics, dynamics and renormalization (World Scientific Publishing Company, 2000).[29] K. Sugiyama, R. Ni, R. J. Stevens, T. S. Chan, S.-Q. Zhou, H.-D. Xi, C. Sun, S. Grossmann, K.-Q. Xia, and D. Lohse,Flow reversals in thermally driven turbulence, Phys. Rev. Lett. , 034503 (2010).[30] M. Chandra and M. K. Verma, Flow reversals in turbulent convection via vortex reconnections, Phys. Rev. Lett. ,114503 (2013).[31] B. Podvin and A. Sergent, Precursor for wind reversal in a square Rayleigh-B´enard cell, Phys. Rev. E , 013112 (2017).[32] H. B´enard, Etude exp´erimentale du mouvement des liquides propageant de la chaleur par convection. R´egime permanent:tourbillons cellulaires, C. r. hebd. s´eances Acad. sci. Paris , 1004 (1900).[33] L. Rayleigh, On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side, Phil.Mag. S. , 529 (1916).[34] H. B´enard, Sur les tourbillons cellulaires et la th´eorie de Rayleigh, C. r. hebd. s´eances Acad. sci. Paris , 1109 (1927).[35] A. Oberbeck, ¨Uber die W¨armeleitung der Fl¨ussigkeiten bei Ber¨ucksichtigung der Str¨omungen infolge von Temperaturdif-ferenzen, Ann. Phys. Chem. , 271 (1879).[36] J. Boussinesq, Th´eorie analytique de la chaleur mise en harmonic avec la thermodynamique et avec la th´eorie m´ecaniquede la lumi`ere , Vol. II (Gauthier-Villars, 1903).[37] D. J. Tritton,
Physical fluid dynamics (Springer, 2012).[38] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. Mcrae, G.-T. Bercea, G. R. Markall, and P. H. J.Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Trans. Math. Softw. , 1 (2016).[39] G. W. Stewart, A Krylov–Schur algorithm for large eigenproblems, SIAM J. Matrix Anal. A. , 601 (2002).[40] V. Hernandez, J. E. Roman, and V. Vidal, SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems,ACM Trans. Math. Softw. , 351 (2005).[41] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Oxford University Press, 1961).[42] P. Hirschberg and E. Knobloch, Mode interactions in large aspect ratio convection, J. Nonlinear Sci. , 537 (1997).[43] S. Fauve, Henri B´enard and pattern-forming instabilities, C. R. Phys. , 531 (2017).[44] J. Mizushima, Onset of the Thermal Convection in a Finite Two-Dimensional Box, J. Phys. Soc. Japan , 2420 (1995).[45] E. D. Siggia, High Rayleigh number convection, Annu. Rev. Fluid Mech. , 137 (1994).[46] Y. A. Kuznetsov, Elements of applied bifurcation theory , Vol. 112 (Springer Science & Business Media, 2013).[47] L. P. Kadanoff, Turbulent heat flow: structures and scaling, Phys. Today , 34 (2001).[48] F. Ravelet, L. Mari´e, A. Chiffaudel, and F. Daviaud, Multistability and memory effect in a highly turbulent flow: Experi-mental evidence for a global bifurcation, Phys. Rev. Lett. , 164501 (2004).[49] M. Berhanu, R. Monchaux, S. Fauve, N. Mordant, F. P´etr´elis, A. Chiffaudel, F. Daviaud, B. Dubrulle, L. Mari´e, F. Ravelet, et al. , Magnetic field reversals in an experimental turbulent dynamo, EPL-Europhys. Lett. , 59001 (2007).[50] O. Cadot, A. Evrard, and L. Pastur, Imperfect supercritical bifurcation in a three-dimensional turbulent wake, Phys. Rev.E , 063005 (2015).[51] V. Dallas, K. Seshasayanan, and S. Fauve, Transitions between turbulent states in a two-dimensional shear flow, Phys.Rev. Fluids , 084610 (2020).[52] P. Winchester, V. Dallas, and P. D. Howell, Zonal flow reversals in two-dimensional Rayleigh-B´enard convection, arXivpreprint arXiv:2010.15642 (2020).[53] D. Sondak, L. M. Smith, and F. Waleffe, Optimal heat transport solutions for Rayleigh–B´enard convection, J. Fluid Mech. , 565 (2015).[54] F. Waleffe, A. Boonkasame, and L. M. Smith, Heat transport by coherent Rayleigh-B´enard convection, Phys. Fluids ,051702 (2015).[55] P. E. Farrell and P. A. Gazca-Orozco, An augmented Lagrangian preconditioner for implicitly-constituted non-Newtonianincompressible flow, SIAM J. Sci. Comput. , B1329 (2020). [56] P. E. Farrell, L. Mitchell, and F. Wechsung, An Augmented Lagrangian Preconditioner for the 3D Stationary IncompressibleNavier–Stokes Equations at High Reynolds Number, SIAM J. Sci. Comput.41