Coherent Solutions and Transition to Turbulence in Two-Dimensional Rayleigh-Bénard Convection

Thermal convection is heat transfer in liquids and gases resulting from buoyancy-driven fluid motion, and is fundamental for establishing circulations in the earth’s atmosphere and oceans, planetary mantle dynamics, stellar evolution, as well as a broad range of engineering applications [3]. Rayleigh-Bénard convection (RBC) refers to convection generated by simplified dynamical equations and idealized boundary conditions. The simplified equations are called the Oberbeck-Boussinesq approximation, describing viscous fluids that are mechanically incompressible but thermally expansible [see, e.g., 4, and references therein]. The domain is confined between two parallel plates, heated from below and cooled from above such that the bottom plate is always warmer than the top plate by a fixed temperature difference. The Rayleigh-Bénard setup presents a framework that is at the same time experimentally realizable, analytically tractable and numerically accessible [see 5, 6, for overviews]. Its prominence in the literature stems partly from its practical importance, and partly from the fact that Rayleigh-Bénard convection represents a rich nonlinear dynamics to test physical reasoning and mathematical techniques.

In the Rayleigh-Bénard convection problem, the Rayleigh number $Ra$ and the Prandtl number $Pr$ are nondimensional parameters that control the flow dynamics: $Ra$ characterizes the relative strength of buoyancy-driven inertial forces to viscous forces and $Pr$ is the ratio of the kinematic viscosity to the thermal diffusivity. Furthermore, a key nondimensional diagnostic parameter is the Nusselt number $Nu$, measuring the vertical heat transport for a given temperature difference between the plates. Over the last several decades, there has been an intense focus on the scaling behavior of $Nu$ on $Ra$ and $Pr$. This relationship has been explored via experiments [7, 8, 9, 5, 10, 11, 12], scaling laws [13, 14, 15, 16, 17, 18, 19, 20], rigorous upper bounds on the allowable heat transport [21, 22, 23, 24, 25], and numerical simulations [26, 27, 1, 2, 28, 29, 30, 31, 32, 33]. The two competing scaling behaviors that emerge from this body of work are $Nu\propto Ra^{1/3}$ and $Nu\propto Ra^{1/2}$. The $1/3$ scaling, often referred to as the classical scaling, emerges from physical arguments [13, 14] involving independence of top and bottom boundary layers [14] as well as marginal stability of the boundary layers [13, 17]. The $1/2$ scaling, recently referred to as the ultimate regime, emanates from mixing length theories of turbulence [16, 15]. There has been a flurry of recent activity on the possible transition from scaling $Nu\propto Ra^{1/3}$ to scaling $Nu\propto Ra^{1/2}$, e.g. experiments by [7, 34, 11, 12]. Several studies have also been aimed at accessing the regime for $Ra>10^{13}$, and at verifying (or not) the existence of the scaling $Nu\propto Ra^{1/2}$ [35, 12, 36, 37, 38, 39, 40, 30]. Such considerations are beyond the scope of our work in this paper. We note that the flow Reynolds number can be expressed in terms of $Ra$ and $Pr$, and that the precise relationship is also of interest to researchers [19, 41, 42]. Interestingly, rigorous upper bounds on heat transport are sensitive to the boundary conditions, being bounded by $1/2$ in the no-slip case [23] and by $5/12$ in the 2D free-slip case [43, 44]. We also mention that substantial effort has been devoted to determining the effect of $Pr$ on the heat transport [45, 46, 47, 48, 49, 50]. Finally, the Grossmann-Lohse (GL) framework [19, 20] introduces nonlinear relationships for $Nu\left(Ra,Pr\right)$ and $Re\left(Ra,Pr\right)$ with six free parameters that have been fit to experimental data [42].

In [1, 2], the 2D Oberbeck-Boussinesq equations were solved numerically to find steady solutions. All of the computed solutions consist of a symmetric pair of hot and cold plumes emanating from the bottom and top boundary layers, respectively, and extending nearly wall-to-wall. The ‘primary’ solution bifurcates from the conduction state at $Ra\approx 1708$, and has a characteristic aspect ratio of $\Gamma_{prim}=L_{prim}/{H_{prim}}\approx 2$, where $L_{prim}$ and $H_{prim}$ are the length and height of the domain, respectively. By adjusting the horizontal length of the domain, $L$, [1, 2] also looked for maximal heat transport solutions, with characteristic aspect ratio $\Gamma_{max}<2$, decreasing with increasing $Ra$. For each value of $Ra$ and $Pr$ considered, they computed solutions that maximized heat transport, referred to as the ‘optimal solution’ (with highest $Nu$). For higher $Ra$, two local maxima emerged with different aspect ratios, and $Pr$ determined which maximum was the global optimal heat transport solution [2]. Best fit scaling relations $Nu-1\propto Ra^{\beta}$ produce $\beta\approx 0.315$ ($\beta\approx 0.311$) for the optimal solution with $Pr=7$ ($Pr=100$), and $\beta\approx 0.28$ ($\beta\approx 0.227$) for the primary solution with $Pr=7$ ($Pr=100$).

The work presented herein aims to establish a direct link between the steady solutions [1, 2], and the structures observed in time-evolving Rayleigh-Bénard flows. As motivation for the study, Figure 1 shows the Rayleigh-Bénard simulation data in a two-dimensional (2D) domain of aspect ratio $\Gamma=10$, along with the best-fit scalings for the primary and optimal solutions. With the expection of the optimal fit at $Pr=7$, which starts at $Ra=10^{7}$, all fits start at $Ra\approx 10^{6}$ and use the available data for higher $Ra$. The fit for the primary solutions have been extended beyond the highest available $Ra\approx 5\times 10^{7}$ to aid the eye. For both $Pr=7,100$ at low $Ra<10^{7}$, the $Nu$ values for the simulation data are closer to the values corresponding to the primary solution. In the range $10^{7}<Ra<10^{9}$, the $Nu$ values for the data consistently become higher than for the primary solution, and the optimal solutions provide a tight upper bound on both sets of data. For $Pr=100$, it is especially evident that $Nu$ values for the data are in between primary and optimal $Nu$ values, with approximate scaling exponent $0.293$ closer to the optimal exponent $0.311$ than to the primary exponent $0.227$.

Furthermore, Figure 2 shows the appearance of all three steady solutions — primary, optimal, and local maximum — within a single time snapshot of a simulation with domain aspect ratio $\Gamma=10$, $Ra=1.1\times 10^{5}$, and $Pr=100$. The time $t=1742$ is relatively early in the evolution from initial conditions (12) with random amplitudes. Notice that these initial conditions do not select a scale in the horizontal direction, and that the horizontal scales of the primary, optimal and local max arise spontaneously from nonlinear interactions. At this early time $t=1742$, the flow is quasi-steady, and later transitions to a quasi-periodic, statistically steady state. In Figure 2, the middle left, middle right and bottom panels compare, respectively, the optimal, local maximum and primary solution to different simulation sub-boxes. In each case, the sub-box has the same horizontal scale as the corresponding steady solution, e.g., the primary solution has horizontal scale $L_{prim}=2l_{p}\approx 4$.

Motivated by Figures 1 and 2, the purpose of the current investigation is to explore how transition to turbulence in 2D Rayleigh-Bénard convection is influenced by the primary and optimal steady solutions. We focus on Rayleigh numbers $10^{5}<Ra<10^{9}$ well above the onset of convection and well below the ultimate regime. In this range, the flow is potentially described as a dynamical systems ‘repeller’ consisting of unstable solutions that are sampled by turbulence trajectories [51]. In such a description, solutions with a small number of unstable directions in phase space would be sampled most frequently. With increasing $Ra$, one would expect more and more unstable (and unsteady) solutions of increasing complexity in the phase space, thus complicating such an investigation at higher $Ra$. The possibility of describing turbulence as a repeller has been previously explored in the context of wall-bounded shear flows (see [52, 53, 51, 54, 55, 56, 57] and references therein).

We restrict our investigation to $Pr=7$ and $Pr=100$, in part because their optimal solutions are distinct from one another [2]. The optimal thermal fields of $Pr=7$ cases have coiling arms and a larger wavelength, as compared to those of $Pr=100$ which are armless and columnar [1, 2]. Furthermore, $Pr=7$ has the added benefit of having been investigated extensively in the literature. Despite the differences in the nature of the optimal solution for $Pr=7$ and $Pr=100$, the scaling behavior of $Nu$ with $Ra$ is almost identical (see Figure 1). We perform well-resolved numerical simulations and analyze the resulting datasets using two main analysis techniques. First we consider the spatial correlation between the computed steady-state solutions (primary or optimal) and ‘windows’ of the simulation data with the same aspect ratio as the steady-state solution. Second, we use the singular value decomposition on the windows with high correlation to compare the underlying structures. Our work fits within a broader class of data-driven approaches to explore structures and dynamics in Rayleigh-Bénard convection. For example, the dynamic mode decomposition has been used to search turbulence data for unstable periodic orbits [58]. Turbulent superstructures with scale separation were identified in [59]. In [29], a convolutional neural network was trained and used to analyze turbulent superstructures. Constrained generative adversarial neural networks have been used to learn the underlying distribution of a Rayleigh-Bénard convection dataset [60].

The rest of the paper is organized as follows. Section II presents the governing equations and the details of the numerical simulations. Section III describes the methodologies we used to detect the footprint of the steady solutions, and the rationales behind them. The results are presented in Section IV, followed by discussion in Section V.

When considering strategies to detect steady solutions in the time developing simulations, we first explored approaches to search for the signature of the optimal solutions. An obvious strategy is to monitor the global heat transfer and focus on high instantaneous $Nu_{-}\left(t\right)$ events. Then, extracting flow fields at individual times corresponding to high $Nu_{-}\left(t\right)$, one may inspect these fields for signs of structural similarity with the optimal solutions. Although this method is successful at low $Ra$, it has a number of drawbacks. First, the global $Nu$ is governed by the combined contributions of all the plumes present in the field at any given time, and therefore may not be the best indicator of a close match between the optimal solution and a local plume in the large-aspect-ratio box. Furthermore, we wish to employ a more general technique to detect different steady solutions (optimal, local max, and primary).

To address these concerns, we developed a windowing technique to search locally for all of the known steady solutions, at all times, keeping $Ra$ and $Pr$ fixed. For fixed time, Figure 5 illustrates the idea using a schematic of the optimal solution for $Pr=7,Ra=1.1\times 10^{5}$.

By sliding the window across the length of the box, one can compute the spatial correlation between a steady solution and the sub-field centered at $x$, given by

where $\left<\cdot,\cdot\right>$ is an $L^{2}$ inner product and $T_{steady}$ is the temperature field of the steady solution. The notation $T_{sub}$ refers to the temperature field in a portion of the computational domain, which is centered at location $x$ with horizontal extent $2l_{steady}$. Repeating this for every snapshot of the temperature field, one obtains the alignment $\cos\left(\theta\left(x,t\right)\right)$ as a function of position $x$ and time $t$. This spatial correlation data can be used to approximate a probability density function, from which we can quantify the likelihood of convective plumes that are highly correlated with a particular steady solution. For a given time regime, e.g. the spin-up period or statistically steady state regime, we counted the number of sub-fields with correlation $0.8\leq\cos\left(\theta\right)\leq 1.0$. We define the incidence parameter $\mathcal{I}$ as

where $M_{\left[0.8,1.0\right]}$ is the number of sub-fields within the specified correlation range, and $M_{\text{tot}}$ is the total number of sub-fields. A few representative values of the incidence parameter $\mathcal{I}$ for the primary solution are given in Table 2 and for the optimal solution in Table 3. These numbers will help to inform the discussion in Section IV.

The windowing technique allows us to identify the influence of multiple steady states on the structure and statistics of transitional flows, especially with respect to horizontal scale selection by nonlinear interactions. For this first study, beyond our lowest $Ra\approx 10^{5}$, we focus on the primary and optimal solutions. We observe that several of their features can remain intact for $10^{7}\lesssim Ra\lesssim 10^{8}$ at both $Pr$ considered in this work. It is remarkable that their signatures endure, since turbulence significantly disrupts their wall-to-wall nature, along with their upright plumes, and in some cases, their delicate interior structures. Once highly-correlated sub-fields have been identified, singular value decomposition [77] is used for more detailed comparison between steady solutions and instantaneous simulation sub-fields. We find that the first two SVD modes are associated with the boundary layer and plume structures, respectively, and provide a quantitative means to separately assess agreement for these two structural features.

To set the stage, we begin this section by reviewing some relevant literature describing the dynamics of similar $Ra-Pr$ regimes studied herein. In [78], the authors employed a low dimensional model using the most energetic modes from direct numerical simulations to study 2D flow regimes at $Pr=6.8$. According to [78], the flow is chaotic for $r_{c}>48.4$, where $r_{c}=Ra/Ra_{c}$ and $Ra_{c}\approx 1708$. Experimental studies in [79] investigated the many routes the flow takes to reach a turbulent state, and the dependence of these routes on aspect ratio $\Gamma$, $Pr$, and the presence of a mean flow. Their results for $\Gamma=3.5$ at $Pr=5$ show that transition to a non-periodic state occurs for $r_{c}>50$. For $Pr=7$, the lowest value of $r_{c}$ investigated in our simulations is $r_{c}=64.4$, expected to be well above the threshold for chaos by comparison to [78, 79]. In the experimental studies [34] with $\Gamma=1$ at $4\leq Pr\leq 1350$, the flows with $Ra\geq 2\times 10^{7}$ have been described as being turbulent.

We also conducted studies of our own to classify the nature of our simulation fields. For example, we examined the time traces of vertical velocity and temperature, collected from probes placed in the bulk and the boundary layer. With the exception of the lowest $Ra=1.1\times 10^{5}$ for $Pr=100$, all cases exhibited spectra with a broad wavenumber and frequency distribution in statistically steady state (see Figure 4).

We next present our results, organized into two separate sections for $Pr=7$ and $Pr=100$ (Section IV.1 and IV.2, respectively). As we will demonstrate, primary and optimal structures are readily observed in temperature fields corresponding to transitional values of $Ra$. Visualizations and data analysis of the simulations for higher $Ra$ suggest that a signature of the primary solution persists for $Pr=7$, consistent with agreement for $Nu$ vs. $Ra$ scaling seen in Figure 1. The optimal solution for higher $Ra$ is easier to detect in the case of $Pr=100$, perhaps in part because of its simpler structure compared to the optimal solution for $Pr=7$ (see Figure 13).

We investigated a transitional range of Rayleigh numbers $10^{5}<Ra\lesssim 10^{8}$ for $Pr=7$ and $Pr=100$ to find signatures of 2D steady solutions that were first described in [1] and [2]. These steady solutions satisfy no-slip boundary conditions in the wall-normal direction, periodic boundary conditions in the horizontal direction, and were found by numerical continuation as bifurcations from the conduction state at $Ra\approx 1708$. The primary solution has aspect ratio $\approx 2$, and becomes unstable at $Ra\approx 53,000$ to a time-periodic solution. The optimal solution maximizes heat transport over steady solutions, and has aspect ratio $<2$ that decreases with increasing $Ra$. A third type of steady solution corresponds to a local maximum of the heat transport, also with aspect ratio smaller than $2$. These solutions impose mirror symmetry about $x=0$, and generally speaking, their temperature fields consist of one hot/cold plume pair, with the hot plume emanating from the bottom boundary layer and centered at $x=0$. On the other hand, the three solution types differ with respect to details such as horizontal scale, boundary layer thickness and interior structure.

A domain aspect ratio $\Gamma=10$ was chosen for the simulations in the current work, allowing multiple copies of coherent structures reminiscent of the steady solutions to appear spontaneously from nonlinear interactions, with minimal constraint by the domain dimensions. The simulation data for $Nu$ vs. $Ra$ is almost coincident with the data for the primary solution at $Pr=7$ at these relatively low values $10^{5}<Ra<10^{8}$, especially at the lower end of the range (Figure 1). For $Pr=100$, the primary solution becomes less relevant with increasing $Ra$, and the optimal scaling appears to be more dominant. Thus, from the outset, the statistical $Nu$ data strongly suggested a structural relevance of the primary solution for transition to turbulence. At the same time, the tight upper bound on $Nu$ provided by the optimal solution inspired us ask if the optimal solution might also influence the features of transitional and turbulent data (Figures 1, 6 and 14).

Our objective was to establish further links between the simulation data and the steady solutions, beyond the $Nu$ vs. $Ra$ information, and for two regimes of $Pr$ represented by $Pr=7$ and $Pr=100$. In particular, we used a moving window technique to identify simulation sub-fields with high correction coefficient (13), considering all three types of steady solutions for the comparison at $Ra\approx 10^{5}$. As evidenced by Figures 2, 7, and 15, the results are stunning, with all three steady solutions appearing prominently during transition to statistically steady state. Then, as the flows enter the statistically steady time regime, a modulated version of the primary solution is dominant. For higher $Ra\gtrsim 10^{5}$, we focused on the primary and optimal solutions, leaving further study of local maximal transport solutions for future work.

For $Ra$ in the two decades $10^{7}\leq Ra<10^{9}$, Schlieren-type plots – Figure 9 ($Pr=7)$ and Figure  16 ($Pr=100$) – show a multi-scale horizontal structure in the simulation snapshots. Visual inspection suggests that the primary solution sets the large horizontal scale of turbulent plumes, and the intriguing possibility for the optimal solution to impact the small-scale plumelets that converge together in a primary-scale updraft. We investigated these large-scale and small-scale features using singular value decomposition on highly-correlated sub-fields and, respectively, their corresponding primary and optimal solutions.

Consistent with the $Nu$ vs. $Ra$ plots in Figure 1, we found subtle differences in the dominance of the primary solutions at large scales, and the persistence of the optimal solutions at small scales. For $Pr=7$, the windowing technique and SVD analyses clearly identify the primary solution as embedded within the turbulence at $Ra\approx 10^{7}$ (Figure 12), while the signature of the optimal solution is much less apparent. We note the lack of scale separation between primary and optimal solutions for $Pr=7$ at the moderate value $Ra\approx 10^{7}$ (Figure 13). On the other hand, there is a definitive scale separation between primary and optimal solutions for $Pr=100$, $Ra\approx 10^{7}$, and the analyses favor the optimal solutions (Figure 18). It is conceivable that larger scale separation between primary and optimal solutions will facilitate detection techniques aimed at higher $Ra$ regimes.

An interesting possibility is presented by the high-$Ra$ 2D simulation data of [38] in the range $Ra=[10^{8},10^{14}]$ In a box of aspect ratio $\Gamma=2$ and for $Pr=1$, [38] find an approximate scaling $Nu$ vs. $Ra^{\beta}$ with $\beta\approx 0.29$ for $Ra=[10^{8},10^{13}]$, beyond which they report a transition to $\beta\approx 0.35$ [38, 39, 81]. Their data is reproduced here in Figure 19, along with (i) our lower $Ra$ data for $\Gamma=10$, $Pr=7$ and $Pr=100$; (ii) the primary solution best-fit scaling at $Pr=7$ extended to the range $Ra=[10^{8},10^{15}]$; and (iii) the optimal solution best-fit scaling at $Pr=7$ extended to the range $Ra=[10^{8},10^{15}]$. One can see a seamless transition between our 2D data for $Ra=[10^{5},10^{8}]$ and the 2D data for $Ra\geq 10^{8}$ [38]. Furthermore, the continuation of the optimal scaling with $\beta\gtrsim 0.3$ begs the question as to whether or not the transition in 2D exponent at $Ra=10^{13}$ [38] could be related to boundary-layer structures with horizontal and vertical scales determined by the optimal solution. Thus, it would undoubtedly be revealing to pursue higher values of $Ra$ in 2D, along with application of other types of data analyses ([82, 83, 29, 84, 60, 85]), in the exploration of exact coherent solutions.