2.1. Problem set-up

We consider the natural convection flow of an air–$\textrm {H}_2\textrm {O}$–$\textrm {CO}_2$ mixture confined in a cubic Rayleigh–Bénard cell, with top and bottom isothermal walls at $T_{cold}$ and $T_{hot}$ and adiabatic vertical walls. The six walls are characterised by uniform grey emissivities $\varepsilon$, the horizontal isothermal walls being black ($\varepsilon =1$) and the vertical adiabatic walls being perfectly diffuse reflecting ($\varepsilon =0$). The parameter controlling the flow regime is the Rayleigh number defined by

(2.1)\begin{equation} Ra = \frac{g\beta\Delta T L^3}{\nu_f a}, \end{equation}

where $g$ is the gravitational acceleration, $\beta =1/T_0$ is the thermal expansion coefficient ($T_0$ being the mean temperature), $L$ is the size of the cavity, $\nu _f$ is the kinematic viscosity, $a$ is the thermal diffusivity and $\Delta T=T_{hot}-T_{cold}$ is the temperature difference between hot and cold walls. We investigate the Rayleigh range $Ra\in [10^6\text {--}10^8]$ in which reorientations are likely to be observed (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010). We vary the Rayleigh number by changing the temperature difference $\Delta T$ and other parameters are fixed. In order to make our study relevant for building applications, we consider an air–$\textrm {H}_2\textrm {O}$–$\textrm {CO}_2$ mixture at a mean temperature of $T_0=300\ \textrm {K}$, of molar composition $X_{{H}_2{O}}=0.02$ and $X_{{CO}_2}=0.001$, and the cavity size is set to $L=1\ \textrm {m}$, so that the temperature difference vary in the range $\Delta T\in [1.1\times 10^{-2}\text {--}1.1]\ \textrm {K}$. Thermophysical properties are assumed to be uniform (low temperature differences), not affected by the small amount of water vapour and carbon dioxide, and constant at all Rayleigh numbers (thermal conductivity $\lambda =0.0263\ \textrm {W}\,\textrm {m}^{-1}\,\textrm {K}^{-1}$, thermal diffusivity $a=2.25\times 10^{-5}\ \textrm {m}^2\,\textrm {s}^{-1}$, Prandtl number $Pr=\nu _f/a=0.707$).

Mass, momentum and energy balance are made dimensionless using the cavity size $L$, the reference time $L^2/(a\sqrt {Ra})$ and the reduced temperature $\theta =(T-T_0)/\Delta T$ and write under Boussinesq approximation as

(2.2)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} =0, \end{gather}

(2.3)\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla}p +Pr \theta \boldsymbol{e_z} + \frac{Pr}{\sqrt{Ra}}\nabla^2 \boldsymbol{u}, \end{gather}

(2.4)\begin{gather} \frac{\partial \theta}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta = \frac{1}{\sqrt{Ra}}\left(\nabla^2 \theta+\mathcal{P}^{rad}\right), \end{gather}

where $\boldsymbol {u}=(u,v,w)$ is the dimensionless velocity vector and $p$ the dimensionless motion pressure. The velocity is zero on all walls. The temperature is set to 0.5 and $-0.5$, respectively, on the bottom and top walls. The conductive flux $-\boldsymbol {\nabla }\theta \boldsymbol {\cdot } \boldsymbol {n}$ is zero on the four lateral walls, as these walls have a zero radiative emissivity.

The dimensionless radiative power $\mathcal {P}^{rad}$ in (2.4) accounts for absorption and emission of thermal radiation by the medium and is defined by

(2.5)\begin{equation} \mathcal{P}^{rad}(\boldsymbol{r}) = \frac{L^2}{\lambda \Delta T} \int_\nu \kappa_\nu \left( \int_{4{\rm \pi}} I_\nu(\boldsymbol{\varOmega},\boldsymbol{r})\,\textrm{d}\boldsymbol{\varOmega} - 4{\rm \pi} I_\nu^\circ (T(\boldsymbol{r})) \right)\,\textrm{d}\nu, \end{equation}

where $I_\nu (\boldsymbol {\varOmega },\boldsymbol {r})$ is the actual radiative intensity at frequency $\nu$, direction $\boldsymbol {\varOmega }$ and position $\boldsymbol {r}$. $I_\nu ^\circ (T(\boldsymbol {r}))$ is the Planck equilibrium intensity at temperature $T$ and $\kappa _\nu$ is the absorption coefficient of the medium. In accordance with the Boussinesq approximation, the absorption coefficient is assumed to be spatially uniform. The Planck equilibrium intensity is given by

(2.6)\begin{equation} I_\nu^\circ (T(\boldsymbol{r}))=\frac{2h\nu^3}{c_0^2}\frac{1}{\exp\left(\dfrac{h\nu}{k_BT(\boldsymbol{r})}\right)-1}, \end{equation}

where $h$ is the Planck constant, $k_B$ is the Boltzmann constant and $c_0$ the speed of light in vacuum. The radiative intensity field is obtained by solving the radiative transfer equation

(2.7)\begin{equation} \boldsymbol{\varOmega}\boldsymbol{\cdot}\boldsymbol{\nabla}I_\nu(\boldsymbol{\varOmega},\boldsymbol{r})=\kappa_\nu \left( I_\nu^\circ (T(\boldsymbol{r})) -I_\nu(\boldsymbol{\varOmega}, \boldsymbol{r})\right). \end{equation}

The associated boundary condition at boundary points $\boldsymbol {r}^w$ for grey diffuse reflecting walls writes as

(2.8)\begin{equation} I_\nu(\boldsymbol{\varOmega},\boldsymbol{r}^w)=\varepsilon I_\nu^\circ (T(\boldsymbol{r}^w))+\frac{1-\varepsilon}{\rm \pi} \int_{\boldsymbol{\varOmega}'\boldsymbol{\cdot} \boldsymbol{n}<0} I_\nu(\boldsymbol{\varOmega}',\boldsymbol{r}^w)\vert \boldsymbol{\varOmega}'\boldsymbol{\cdot} \boldsymbol{n} \vert\,\textrm{d}\boldsymbol{\varOmega}', \end{equation}

for directions $\boldsymbol {\varOmega }$ such that $\boldsymbol {\varOmega }\boldsymbol {\cdot } \boldsymbol {n}>0$, $\boldsymbol {n}$ being the wall normal directed towards the inside of the domain.

It is worth noting that the flow equations are written and solved in dimensionless form, while the radiative transfer equations are treated in dimensional form since we consider an actual molecular radiating gas. A key parameter for radiation is the optical thickness $\tau _\nu =\kappa _\nu L$ that varies over several orders of magnitude in our model. Considering a grey fluid (wavelength-independent absorption) would facilitate a parametric study of radiation effects but would fail to represent the behaviour of actual radiating gases.

2.2. Numerical methods

The numerical methods used for solving the coupled system of (2.2)–(2.4) and (2.7) are presented and validated in detail by Soucasse, Rivière & Soufiani (Reference Soucasse, Rivière and Soufiani2016) and Soucasse et al. (Reference Soucasse, Podvin, Rivière and Soufiani2020). We briefly mention here the main features of the coupled algorithm.

Navier–Stokes equations are solved using a Chebyshev collocation method (Xin & Le Quéré Reference Xin and Le Quéré2002). Domain decomposition along the vertical direction is carried out by the Schur complement method to make the computations parallel (Xin, Chergui & Le Quéré Reference Xin, Chergui and Le Quéré2008). Time integration is performed through a second-order semi-implicit scheme. The velocity divergence-free condition is enforced using a projection method. The radiative transfer equation is solved using a ray-tracing algorithm, made parallel by distributing the rays among the different processors. The $4{\rm \pi}$ angular domain is uniformly discretised using 900 rays from volume cell centres and 450 rays from boundary cell centres. The Absorption distribution function (known as ADF) model (Pierrot et al. Reference Pierrot, Rivière, Soufiani and Taine1999) is used to take into account the spectral variations of the absorption coefficient of the air–$\textrm {H}_2\textrm {O}$–$\textrm {CO}_2$ mixture: it consists in substituting the integration over the frequency with an integration over the values of the absorption coefficient, for which a coarse logarithmic discretisation is sufficient. In the present study, the values of the absorption coefficient have been logarithmically discretised in 16 classes and the accuracy of the model has been shown to be better than 1 % (Soucasse et al. Reference Soucasse, Rivière, Xin, Le Quéré and Soufiani2012). Model parameters and computational details for the considered mixture are given by Soucasse (Reference Soucasse2013) and Soucasse et al. (Reference Soucasse, Rivière, Xin, Le Quéré and Soufiani2012).

Direct numerical simulations have been performed, considering the air–$\textrm {H}_2\textrm {O}$–$\textrm {CO}_2$ mixture as radiating (coupled case) or transparent (uncoupled case, $X_{{H}_2{O}}=0$, $X_{{CO}_2}=0$ and thus $\mathcal {P}^{rad}=0$). Simulation parameters are given in table 1. Five different Rayleigh numbers have been considered: $Ra=\{ 10^6 ; 3\times 10^6 ; 10^7 ; 3\times 10^7 ; 10^8\}$. The convection mesh is built from Chebyshev collocation points. We have checked that the number of points in the boundary layers is sufficient regarding the criterion proposed by Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). Up to $Ra=10^7$, the radiation mesh is built from the convection mesh, coarsened by a factor of two in each direction of space. For $Ra=3\times 10^7$ and $Ra=10^8$, the radiation mesh is coarsened by a factor of four in each direction of space compared with the convection mesh and we use a radiation subgrid model (Soucasse, Rivière & Soufiani Reference Soucasse, Rivière and Soufiani2014b) to account for the radiation of small spatial scales. This subgrid model has been validated in various configurations and its accuracy is approximately a few per cent on radiative power and wall fluxes. It has been used for the simulation of coupled natural convection and radiation in a differentially heated cavity at $Ra=3\times 10^9$ (Soucasse et al. Reference Soucasse, Rivière and Soufiani2016). Finally, an explicit coupling is carried out between flow and radiation calculations and the radiation source term is updated every 10 convection time steps $\delta t$ (the flow time step is imposed by numerical stability constraints and does not correspond to significant variations of the temperature field). Time integration is carried out over a period $\Delta t$ once the asymptotic regime (statistically steady) is reached.

Table 1. Simulation parameters: convection mesh, radiation mesh, convection time step $\delta t$ and integration time interval $\Delta t$ in the asymptotic regime. For the convection mesh, numbers in parenthesis indicate the number of spatial domains times the number of collocation points in the vertical in each domain.

$^{a}$Radiation subgrid model is used; $^{b}\delta t$ is $2\times 10^{-3}$ in the uncoupled simulation and $1.5\times 10^{-3}$ in the coupled simulation.

It should be mentioned here that radiation calculations are much more computationally expensive than convection calculations (the CPU time is approximately 30 times greater in the coupled case).

2.3. Radiative transfer effects

Radiative transfer effects on RBC at $Ra=10^7$ have been discussed in a previous work (Soucasse et al. Reference Soucasse, Podvin, Rivière and Soufiani2020). When the fluid emits and absorbs radiation, heat transfer is no longer restricted to the boundary layer region. An energy exchange between convection and radiation in the core of the cavity leads to a significant increase of the convective flux compared with the uncoupled case. The LSC is strengthened and both mean and turbulent kinetic energies increase. Temperature fluctuations also increase but to a lesser extent because of radiative damping.

The same effects are observed at other Rayleigh numbers. The total kinetic energy $e_k=\int 0.5\left \langle \boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {u} \right \rangle \,\textrm {d}\boldsymbol {r}$, the total thermal energy $e_\theta =\int 0.5\left \langle \theta \theta \right \rangle d \boldsymbol {r}$, the conductive flux at the walls $q_{cond}=\int \left \langle \partial \theta /\partial z \right \rangle \,\textrm {d}x\,\textrm {d}y$ and the convective flux in the core ($z=0.5$) $q_{conv}=\sqrt {Ra}\int \left \langle w\theta \right \rangle \,\textrm {d}x\,\textrm {d}y$ are displayed in figure 1 as a function of the Rayleigh number (where $\left \langle \cdot \right \rangle$ denotes the time average). It can be noticed that the kinetic energy and the convective flux significantly increase in the presence of radiation while the thermal energy and the conductive flux are not much affected. However, radiation effects diminish with the Rayleigh number. This can be explained by the following scaling analysis. If the temperature dependence of the Planck equilibrium intensity is linearised around the mean temperature $T_0$, an order of magnitude for the dimensional radiative power is $16\kappa _P \sigma T_0^3\Delta T$, where $\kappa _P= \int _\nu \kappa _\nu I_\nu ^\circ (T_0)\,\textrm {d}\nu \times ({\rm \pi} /\sigma T_0^4)$ is the Planck mean absorption coefficient of the mixture, $\sigma$ is the Stefan–Boltzmann constant. Therefore, the radiative source term in the energy balance (2.4) roughly scales as

(2.9)\begin{equation} \frac{\mathcal{P}^{rad}}{\sqrt{Ra}} \thicksim O\left( \frac{1}{\sqrt{Ra}}\frac{16\kappa_P \sigma T_0^3 L^2}{\lambda}\right), \end{equation}

and decreases in $Ra^{-1/2}$ while the order of magnitude of the convective term $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\theta$ remains constant.

Figure 1. Total kinetic energy $e_k$, total thermal energy $e_\theta$, conductive flux at the wall $q_{cond}$ and convective flux in the core $q_{conv}$ as a function of the Rayleigh number for coupled (red squares) and uncoupled (black triangles) cases.

The unsteady dynamics of RBC in a cubic container is characterised by low-frequency reorientations of the LSC (Bai et al. Reference Bai, Ji and Brown2016; Foroozani et al. Reference Foroozani, Niemela, Armenio and Sreenivasan2017), that can be monitored using the time evolution of the $x$ and $y$ components of the angular momentum with respect to the cavity centre $\boldsymbol {r}_0$:

(2.10)\begin{equation} \boldsymbol{L}=\int (\boldsymbol{r} -\boldsymbol{r}_0)\times\boldsymbol{u}\,\textrm{d}\boldsymbol{r}. \end{equation}

Figure 2 shows the time evolution of components $L_x$ and $L_y$ at the different Rayleigh numbers for the coupled and uncoupled cases. Note that time integration is shorter for the two highest Rayleigh numbers.

Figure 2. Time evolution of $x$ and $y$ components of angular momentum $L_x$ (blue lines) and $L_y$ (red line).

The uncoupled case is characterised at each Rayleigh number by quasi-stable diagonal states, with abrupt reorientations between these states. A diagonal state means that the LSC lies in one of the two diagonal planes $x=y$ or $x=1-y$, with clockwise or anticlockwise motion (four diagonal states are available), and is characterised by a non-zero equilibrium value for both $L_x$ and $L_y$ components. Sudden reorientations between two diagonal states occur when either $L_x$ or $L_y$ changes sign, which corresponds to a rotation of ${\rm \pi} /2$ of the LSC around the vertical axis. A detailed description of the reorientation process in the cubic cell and associated flow patterns is provided by Foroozani et al. (Reference Foroozani, Niemela, Armenio and Sreenivasan2017) and Vasiliev et al. (Reference Vasiliev, Frick, Kumar, Stepanov, Sukhanovskii and Verma2019). Although the overall dynamics is similar over the Rayleigh number range, it can be noticed that the stability of the diagonal states increases with the Rayleigh number: the flow spends less time around zero angular momentum and reorientations are less frequent when the Rayleigh number increases. At $Ra=10^6$, the dynamics is more chaotic with quick passing around zero of either $L_x$ or $L_y$ components. At $Ra=10^8$, the oscillation amplitude of $L_x$ and $L_y$ around the equilibrium value is weaker and only one reorientation event is observed during the sequence (only two of the four diagonal states are visited in this case).

In the coupled case, a significant change in dynamics compared with the uncoupled case is noticeable at $Ra=10^6$, where quasi-stable planar states are observed. A planar state means that the LSC lies either in $x$ planes or $y$ planes, with clockwise or anticlockwise motion (four planar states are available), and is characterised by a zero equilibrium value for one of the two components $L_x$ or $L_y$. Figure 3 shows a snapshot of temperature and velocity fields at $Ra=10^6$. In the uncoupled case (figure 3a) a diagonal state is observed as the LSC lies in the diagonal plane $x=1-y$ with $L_x>0$, $L_y>0$. The fluid flows up along the right vertical edge $(x;y)=(0;1)$ and flows down along the left vertical edge $(x;y)=(1;0)$. This main diagonal roll is slightly tilted and small counter-rotating rolls are noticeable in the top right corner $(x;y;z)=(0;1;1)$ and in the bottom left corner $(x;y;z)=(1;0;0)$. In the coupled case (figure 3b) a planar state is observed as the LSC lies in $y$ planes with $L_x\simeq 0$, $L_y<0$. The fluid flows up along the front plane $x=1$ and flows down along the rear plane $x=0$. Counter-rotating structures are noticeable near the top horizontal edge $(x;z)=(1;1)$ and the bottom horizontal edge $(x;z)=(0;0)$. For $Ra\ge 3\times 10^7$, quasi-stable diagonal states are observed with radiation. However, for a given Rayleigh number, the dynamics is more chaotic and reorientations seem to be more frequent when the flow is coupled with radiation. It is worth noticing here that, although stable planar states are not observed in the uncoupled case, they have been reported at low Rayleigh numbers (Puigjaner et al. Reference Puigjaner, Herrero, Giralt and Simó2004).

Figure 3. Instantaneous flow field at $Ra=10^6$ and $t=5000$ for uncoupled (a) and coupled (b) cases. Streamlines and isotherms $\theta =\{ 0;{\pm }0.05; {\pm }0.1\}$. In the uncoupled case (panel (a)) a diagonal state is observed ($L_x>0$, $L_y>0$) while in the coupled case (panel (b)) a planar state is observed ($L_x\simeq 0$, $L_y<0$).

In order to quantify radiative transfer effects on the temporal dynamics we have computed two characteristic frequency scales: the circulation frequency $f_c$ (or circulation time $\tau _c=1/f_c$) and the reorientation frequency $f_r$ (or reorientation time $\tau _r=1/f_r$). The circulation frequency is a high frequency associated with the rotation frequency of the LSC roll. Frequencies are reported in table 2. The circulation frequency is estimated using a reference ellipsoid path length in the diagonal plane and a reference velocity. In the uncoupled case, the circulation frequency is nearly constant with the Rayleigh number as convective time units are used. In the coupled case, the increase of the kinetic energy leads to an increase of the circulation frequency. This increase compared with the uncoupled case diminishes with the Rayleigh number. The reorientation frequency is estimated by tracking the zeros of the filtered time evolution of $L_x$ and $L_y$ to avoid small-scale noise. Data in table 2 confirm the observations of figure 2: the reorientation frequency seems to decrease with the Rayleigh number and, at a given Rayleigh number, seems to be higher when radiation is taken into account.

Table 2. Circulation frequency $f_c$ and reorientation frequency $f_r$ for the uncoupled and coupled cases at different Rayleigh numbers. The value $N/A$ is indicated when it was not possible to obtain a value (less than two switches observed at $Ra=10^8$ – uncoupled case, many rapid switches at $Ra=3\times 10^6$ – coupled case) or when planar flow states are observed at $Ra=10^6$ – coupled case.

It should be noted here that the uncertainty associated with reorientation frequencies is significant, owing to the few reorientations. If we assume a Poisson distribution of reorientation events (as it has been observed for reversals in cylindrical cells Sreenivasan et al. (Reference Sreenivasan, Bershadski and Niemela2002)), the uncertainty on the number of reorientations $N_r$ would be $\pm \sqrt {N_r}$, and the uncertainty on the reorientation frequency $f_r$ would be $\pm \sqrt {(f_r/\Delta t)}$, $\Delta t$ being the integration time. Uncertainties on reorientation frequencies are given in table 2. Relative uncertainties range from 20 % ($Ra=3\times 10^6$, uncoupled case) to 70 % ($Ra=3\times 10^7$, uncoupled case). These uncertainties moderate the conclusions on radiative transfer effects on the reorientation frequency, especially at $Ra=10^7$ where the increase of the reorientation frequency with radiation is not statistically significant. However, the decrease of the reorientation frequency with the Rayleigh number in the uncoupled case seems to be statistically significant in the range $3\times 10^6 \le Ra \le 10^8$ and has been reported in other works (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010).

3.1. Methodology

The POD in fluid mechanics aims at finding an optimal basis of spatial eigenfunctions $\boldsymbol {\phi }_n(\boldsymbol {r})$ to represent an unsteady flow variable vector $\boldsymbol {U}(\boldsymbol {r},t)$ of size $M$ on a low-dimensional subspace. The POD eigenfunctions or POD modes $\boldsymbol {\phi }_n(\boldsymbol {r})$ are solutions of the following eigenvalue problem (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993):

(3.1)\begin{equation} \int \sum_{k=1}^M \left\langle U^m(\boldsymbol{r},t)U^k(\boldsymbol{r'},t) \right\rangle\phi_n^k(\boldsymbol{r'})\,d\boldsymbol{r'}=\lambda_n \phi_n^m(\boldsymbol{r}), \end{equation}

where $\left \langle \cdot \right \rangle$ denotes the time average. They form an orthonormal basis allowing the decomposition

(3.2)\begin{equation} \boldsymbol{U}(\boldsymbol{r},t)= \sum_{n=1}^\infty a_n(t) \boldsymbol{\phi}_n(\boldsymbol{r}), \end{equation}

where the projection coefficients $a_n(t)$ are statistically uncorrelated and their energy is equal to the eigenvalue such that $\left \langle a_n(t) a_m(t) \right \rangle =\delta _{nm}\lambda _n$, where $\delta _{nm}$ is the Kronecker symbol. The eigenvalue associated with a POD mode is thus a measure of its energy content and the objective is to restrict the decomposition (3.2) to a few modes with the largest eigenvalues so that the associated low-order subspace captures most of the energy of the field $\boldsymbol {U}(\boldsymbol {r},t)$. To take into account the coupling between velocity and temperature fields in thermal convection, we define the flow variable vector as $\boldsymbol {U}=\{ \boldsymbol {u}, \gamma \theta \}$ ($M=4$), where $\gamma$ is a rescaling factor

(3.3)\begin{equation} \gamma = \sqrt{\left\langle \frac{\int \boldsymbol{u}(\boldsymbol{r},t)\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{r},t) \,\textrm{d}\boldsymbol{r}}{\int \theta^2(\boldsymbol{r},t) \,\textrm{d}\boldsymbol{r}} \right\rangle}, \end{equation}

such that velocity and temperature fields have the same energy (Podvin & Le Quéré Reference Podvin and Le Quéré2001). The POD modes thus combine velocity and temperature: $\boldsymbol {\phi }_n=\{\boldsymbol {\phi }_n^u, \gamma \phi _n^\theta \}$.

Equation (3.1) is solved in practice using the method of snapshots (Sirovich Reference Sirovich1987). A snapshot set of 1000 samples $\boldsymbol {U}(\boldsymbol {r},t_i)$ is extracted from each simulation at discrete times $t_i$ with a fixed sampling period of 10 dimensionless time units for $10^6 \le Ra \le 10^7$ and a fixed sampling period of five dimensionless time units for $3\times 10^7 \le Ra \le 10^8$. This sampling period might seem large compared with the transition time between two states and the few reorientations, but according to Podvin & Sergent (Reference Podvin and Sergent2017), precursor events for reorientations are associated with large-scale interactions during relatively large time scales (of the order of the reorientation period). However, we have seen in the time evolution of the angular momentum components in figure 2 that the four possible quasi-stable flow states (planar or diagonal) were not necessarily visited or not equally represented during the time sequence, although there is no physical evidence suggesting that these states are not equiprobable. This is an artefact owing to the relatively short simulation time compared with the time scale separating two reorientations, especially at high Rayleigh number. In order to enforce an equal statistical weight for each flow state and to improve the convergence of the POD method, we have built enlarged snapshot sets, obtained by the action of the symmetry group of the problem on the original snapshot sets. The use of statistical symmetries of the flow in POD is discussed at length by Moin & Moser (Reference Moin and Moser1989) and Holmes, Lumley & Berkooz (Reference Holmes, Lumley and Berkooz1996).

In the uncoupled case, the problem satisfies four independent reflection symmetries $S_x$, $S_y$, $S_z$ and $S_d$ with respect to the planes $x=0.5$, $y=0.5$, $z=0.5$ and $x=y$ (Puigjaner et al. Reference Puigjaner, Herrero, Simó and Giralt2008). In the coupled case, radiative transfer should break the $S_z$ symmetry (radiative emission being proportional to $T^4$), but owing to the weak temperature gradients, nonlinear effects are negligible so that we can consider that the $S_z$ symmetry is still satisfied. These four elementary symmetries generate a symmetry group of 16 elements. This allows us to multiply the number of snapshot by a factor of 16. In conclusion, for each simulation (coupled, uncoupled, $10^6\le Ra \le 10^8$), we work on an enlarged snapshot set made of 16 000 samples (1000 original snapshots multiplied by 16 after applying all the symmetries).

3.2. Energy spectra

The POD spectra obtained in the uncoupled case at different Rayleigh numbers are shown in figure 4(a). Each spectrum is normalised such that $\sum _n\lambda _n=1$ but since the total POD energy is nearly constant with the Rayleigh number in the uncoupled case, it is relevant to compare the spectra between them. The mode ordering roughly corresponds to a ranking of the eigenfunctions in terms of a characteristic spatial scale, and we can therefore associate the low-order modes with the largest spatial scales and the high-order modes with the smallest spatial scales. Examination of figure 4 suggests that the POD spectra can be split into three parts: (i) the large scales for $n\lesssim 10$ which will correspond to the modes retain in the models in § 4; (ii) the intermediate scales for $10 \lesssim n \lesssim N$, with $N$ such that the $N$ first POD modes capture 95 % of the total energy ($N$ increases with the Rayleigh number, it goes from 250 at $Ra=10^6$ to approximately 5000 at $Ra=10^8$); (iii) the small scales for $n\gtrsim N$. At high Rayleigh number ($Ra=3\times 10^7$, $Ra=10^8$), there is a different behaviour of large scales and intermediate scales, with a very fast decay of the low-order modes followed by a very slow decay in the rest of the spectrum owing to the turbulent nature of the flow. On the contrary, at low Rayleigh number ($Ra=10^6$, $Ra=3\times 10^6$), the energy contained in the intermediate scales is proportionately more important and a fast decay of the energy of the small scales is observed. Interestingly, the spectrum decay is roughly log–linear in the intermediate range and this will be used in § 6 to model the variations of the POD spectrum with the Rayleigh number.

Figure 4. The POD eigenspectrum obtained from uncoupled simulations, normalised such that $\sum _n\lambda _n=1$ (a) and POD eigenspectrum ratio between coupled and uncoupled results (b).

Figure 4(b) shows the ratio between the coupled POD spectrum and the uncoupled POD spectrum at each Rayleigh number. This ratio is always greater than one (except for one mode at $Ra=10^6$) which confirms that the coupled POD spectrum captures the energy increase associated with radiation effects. However, this energy increase depends on the Rayleigh number and on the part of the spectrum. At low Rayleigh number, the energy increase is proportionately higher in the large-scale range, while at high Rayleigh number, the energy increase is proportionately higher in the intermediate scale range. In the small-scale range, the spectrum ratio is nearly constant and remains greater than one. As expected from the discussion in § 2.3, the ratio between coupled and uncoupled total POD energies decreases with the Rayleigh number as radiation effects weaken: it is equal to 1.59 at $Ra=10^6$ and equal to 1.19 at $Ra=10^8$.

A last parameter characterising the POD spectra is the factor $\gamma$, the values of which are reported in table 3. This factor $\gamma$ (corresponding to the ratio between mechanical energy and thermal energy, see (3.3)) increases with the Rayleigh number and is always greater in the coupled case at a given Rayleigh number, which is consistent with the observations of figure 1.

Table 3. The $\gamma$ factor for uncoupled and coupled cases at different Rayleigh numbers.

3.3. Spatial eigenfunctions

3.3.1. Uncoupled case, $Ra=10^6$

Figure 5 shows the first 12 POD eigenfunctions at $Ra=10^6$ in the uncoupled case, which correspond to nine spatial structures (three of the modes are doubly degenerated). The figure displays the contribution of each structure to the mean convective heat flux in the vertical direction which can be written as $\sqrt {Ra}\sum _n\lambda _n \phi ^\theta _n\phi ^w_n$ owing to the POD decomposition. Five spatial structures have been already highlighted in previous works (Soucasse et al. Reference Soucasse, Podvin, Rivière and Soufiani2019, Reference Soucasse, Podvin, Rivière and Soufiani2020) and correspond to the first seven POD modes at $Ra=10^7$. They are labelled $M$, $L_{x/y}$, $D$, $BL_{x/y}$ and $C$. We briefly mention their properties and physical meaning.

1. (i) The $M$ mode corresponds to the mean flow: it is made of two counter rotating torus and it is thermally stratified.

2. (ii) The $L_x$ and $L_y$ modes form a pair of degenerate modes (only the $L_x$ mode is shown in figure 5, the $L_y$ mode is its image by a rotation of ${\rm \pi} /2$ around the vertical axis). They correspond to a single large-scale roll lying in either $x$ planes or $y$ planes. When combined, the $L_x$ and $L_y$ modes from a single large-scale diagonal roll.

3. (iii) The $D$ mode is an 8-roll mode that transports fluid from one corner to the other and strengthens the circulation along the diagonal.

4. (iv) The $BL_x$ and $BL_y$ modes (pair of degenerate modes) are constituted of two longitudinal corotating structures around the $x$ axis or the $y$ axis. They connect the core of the cell with the horizontal boundary layers.

5. (v) The $C$ mode is a corner-roll mode which favours planar flow states and was found to promote reorientations between diagonal flow states.

Figure 5. First 12 modes at $Ra=10^6$ (uncoupled case). Streamlines and isosurfaces of the contribution to the mean convective heat flux $\phi ^\theta \phi ^w=0.25$, coloured by mode temperature. For degenerate modes, only the $x$-oriented one is displayed. Colour map for mode temperature ranges from $-0.5$ (blue) to 0.5 (red).

In addition, four other spatial structures are noticeable in figure 5. They are labelled $D^*$, $C^*$, $BL_{x/y}^*$ and $M^*$ because they share common features with modes $D$, $C$, $BL_{x/y}$ and $M$.

1. (i) The $D^*$ mode is a 4-roll mode that shares some similarities with the $D$ mode, as it transports fluid from one corner to the other. However, the rolls of the $D^*$ mode extend from bottom to top, while the rolls of the $D$ mode are confined in either the upper half or the lower half of the cell.

2. (ii) The $C^*$ mode is another corner-roll mode. Unlike the $C$-mode these corner rolls extend from top to bottom and are not confined in a half-cell.

3. (iii) The $BL_x^*$ and $BL_y^*$ modes (pair of degenerate modes) are made of two longitudinal counter-rotating structures around the $x$ or the $y$ axis. A two-dimensional version of these POD modes has been highlighted by Podvin & Sergent (Reference Podvin and Sergent2017) in the square cell (symmetry-breaking mode labelled $S$).

4. (iv) The $M^*$ mode is a single toroidal structure linking the top and bottom walls. The fluid flows up along the adiabatic vertical walls and flows down in the centre of the cell.

3.3.2. Uncoupled case, $Ra=10^8$

Figure 6 shows the first 12 POD eigenfunctions at $Ra=10^8$ in the uncoupled case. Despite the large difference in Rayleigh number, most of the modes observed at $Ra=10^6$ are retrieved. One can identify again modes $M$, $L_{x/y}$, $D$, $BL_{x/y}$, $C$, as well as modes $D^*$, $C^*$ and $BL_{x/y}^*$. However, the mode ordering is not exactly the same. In addition, streamlines and isosurfaces of convective heat flux in figure 6 are closer to the walls than those in figure 5, as boundary layers are thinner when the Rayleigh number increases. The $M^*$ mode is no longer observed, but a new mode, the $L_z$ mode, appears. This mode is constituted of a single roll lying in $z$ plane. This is a pure mechanical mode ($\Vert \boldsymbol {\phi }^u\Vert \simeq 1$) and one can note that the convective heat flux associated with the $L_z$ mode is much weaker than that of the other modes.

Figure 6. First 12 modes at $Ra=10^8$ (uncoupled case). Streamlines and isosurfaces of the contribution to the mean convective heat flux $\phi ^\theta \phi ^w=0.25$ (0.025 for the $L_z$ mode), coloured by mode temperature. For degenerate modes, only the $x$-oriented one is displayed. Colour map for mode temperature ranges from $-0.5$ (blue) to 0.5 (red).

3.3.3. Other coupled and uncoupled cases, similarities

For the other cases studied (other Rayleigh numbers, coupled with radiation), the first 12 POD modes belong to the set of the thirteen POD modes described above: $M$, $L_{x/y}$, $D$, $BL_{x/y}$, $C$, $D^*$, $C^*$, $BL_{x/y}^*$, $M^*$ and $L_z$. The list of the modes, ordered according to their eigenvalue, is given in table 4. A remarkable feature is that the eigenfunctions are globally preserved when radiation is taken into account, although the associated eigenvalues are higher than in the uncoupled case. In the range $10^7 \le Ra \le 10^8$, there are very few differences in the mode ordering, whatever the coupling conditions. It can be noted than the $L_z$ mode becomes more important when the Rayleigh number increases and when radiation is considered. In the range $10^6 \le Ra \le 3\times 10^6$, the mode ordering is less stable. Indeed, the eigenvalues of the modes are closer to each other, given the slower decay of the POD spectrum at low Rayleigh number. The ordering observed in the coupled case at $Ra=10^6$ is the most different from the others: the $D$ mode is no longer in fourth position, the $BL_{x/y}$ modes are also downshift far from the first modes and both $M^*$ and $L_z$ modes are present. This is not surprising as this case is the only one associated with quasi-stable planar flow states.

Table 4. List of the first 12 POD modes for each case studied, ordered from top to bottom. The last line of the table gives the error measure $e^{\mathcal {B}}$ in percentage between the POD basis and the reference POD basis ($Ra=10^7$, uncoupled) as defined in (3.4).

The question arising is: How close are two eigenfunctions of the same nature but belonging to different POD bases (for instance the $D$ mode in the uncoupled case at $Ra=10^6$ and the $D$ mode in the coupled case at $Ra=3\times 10^7$)? To answer this, we have taken the POD basis of the uncoupled case at $Ra=10^7$ as a reference and we have computed, for a given POD basis $\mathcal {B}$, the differences compared with this reference using the following error measure:

(3.4)\begin{equation} e^{\mathcal{B}}=\sqrt{\sum_{n=1}^{12}\frac{1}{12}\left( 1- \left\vert \int \boldsymbol{\phi}_n^{\mathcal{B}}(\boldsymbol{r}) \boldsymbol{\cdot} \boldsymbol{\phi}_n^{ref}(\boldsymbol{r}) \right\vert \,\textrm{d}\boldsymbol{r}\right)}, \end{equation}

where $\boldsymbol {\phi }_n^{ref}$ are the eigenfunctions of the reference POD basis and $\boldsymbol {\phi }_n^{\mathcal {B}}$ are the eigenfunctions of the POD basis $\mathcal {B}$, reordered according to the mode ranking of the reference case. For the $L_z$ mode, which is not contained in the reference basis, we have taken the $L_z$ mode obtained at $Ra=10^7$ in the coupled case as the reference. Values of indicator $e^{\mathcal {B}}$ are reported in table 4. The differences compared with the reference basis are rather small and are approximately 1 % or less, except for the coupled case at $Ra=10^6$. This result will be key for the development of a POD model across the Rayleigh range (see § 6) and for predicting radiative transfer effects.

A last comment on the 13 identified modes can be made regarding their contribution to the global angular momentum and their symmetry properties, given in table 5. The symmetry analysis sheds light on mode interactions in the low-order models and the angular momentum is used to monitor the flow reorientations. The modes contributing to the $x$ and $y$ components of the angular momentum are of course the LSC modes $L_x$ and $L_y$ but also, in a lesser extent the boundary layer modes $BL_x$ and $BL_y$. Furthermore, the $L_x$ and the $BL_x$ modes have the same symmetry properties, as well as the $L_y$ and the $BL_y$ modes. It denotes strong interactions between these two pairs of modes, corresponding to the connection between the LSC and the boundary layers. Interestingly, all the degenerate modes break the $S_d$ symmetry and each pair possesses opposite symmetry properties in $S_x$ and $S_y$. In addition, one can note that $D^*$, $C^*$, $BL_{x/y}^*$ and $M^*$ modes have the same symmetry properties in $S_x$, $S_y$ and $S_d$ as their companion mode $D$, $C$, $BL_{x/y}$ and $M$, but have opposite symmetry properties in $S_z$.

Table 5. Name, angular momentum and symmetry properties of the 13 identified modes. An X indicates a non-zero angular momentum. Here $S_{x/y/z/d}$ and $AS_{x/y/z/d}$ denote, respectively, a symmetry and an antisymmetry with respect to the planes $x=0.5$, $y=0.5$, $z=0.5$ and $x=y$.

4.1. Construction

Proper orthogonal decomposition low-order models are derived from Galerkin projection of Navier–Stokes equations ((2.3)–(2.4)) onto a POD basis set. Using the decomposition (3.2), this yields a system of ordinary differential equations for the mode amplitudes $a_n(t)$ of the form

(4.1)\begin{equation} \frac{\textrm{d} a_n(t)}{\textrm{d} t} = (L_{nm}^B+L_{nm}^D) a_m(t) + Q_{nmp} a_m(t) a_p(t) +T_n(t), \end{equation}

where $L_{nm}^B$ and $L_{nm}^D$ are linear coefficients associated with buoyancy and diffusion and $Q_{nmp}$ are quadratic coefficients associated with advection. They are directly determined from spatial modes $\boldsymbol {\phi }_n$, see Appendix A. Here $T_n(t)$ is a closure term, which models the effect of the truncation and corresponds to an equivalent dissipation term. Following Podvin & Sergent (Reference Podvin and Sergent2015), the closure term is expressed as a combination of linear and cubic terms, and by the addition of noise such that

(4.2)\begin{equation} T_n(t)= L_{nm}^A \left( 1 + \frac{ 1}{\left\langle k \right\rangle} \sum_{p \ge 2}|a_p(t)|^2 \right) a_m(t) +\sigma\epsilon_n(t), \end{equation}

where $L_{nm}^A$ are adjustable linear parameters, $\left \langle k \right \rangle$ is the time average of the energy of the fluctuating modes in the truncation and $\epsilon _n(t)$ is a random noise perturbation of amplitude one.

The adjustable coefficients $L_{nm}^A$ were determined for equilibria $a_n(t)=a_n^{eq}$ corresponding to diagonal states from the DNS ($\textrm {d} a_n^{eq}/\textrm {d}t=0$). We retain the following equilibrium: $a^{eq}_M=\sqrt {\lambda _M}$; $a^{eq}_{L_{x/y}}=\sqrt {\lambda _{L_{x/y}}}$; $a^{eq}_D=\sqrt {\lambda _D}$; $a^{eq}_{BL_{x/y}}=-\eta _{BL}\sqrt {\lambda _{BL_{x/y}}}$ and $a^{eq}_n=0$ for the other modes. Note that boundary layer modes $BL_{x/y}$ are negatively correlated with LSC modes $L_{x/y}$ and their equilibrium value in absolute value $\vert a^{eq}_{BL_{x/y}}\vert$ is smaller than the standard deviation of their fluctuations ($\eta _{BL}<1$). The coefficient $\eta _{BL}$ is approximately 0.15–0.2 except at $Ra=10^8$ where it is much weaker ($\eta _{BL}=0.05$). A large value of $\eta$ represents a strong connection between the boundary layers and the LSC in the diagonal state. The coefficients $L_{nm}^A$ can only be determined unambiguously for modes $M$, $L_{x,y}$, $D$ and $BL_{x,y}$, which are non-zero at equilibrium. For modes which are zero at equilibrium, the sum of linear contributions was adjusted to balance quadratic interactions with the mean mode, in order to avoid excessive energy. We checked that the model behaviour did not change when small changes were made in the adjustable parameters.

The noise level was estimated as follows:

(4.3)\begin{equation} \sigma = C\frac{\sum_{8}^{2000} \lambda_n}{ \sqrt{\sum_2^7 \lambda_n}}, \end{equation}

where $C=0.016$ is a constant, determined at $Ra=10^7$ from a scaling analysis (Soucasse et al. Reference Soucasse, Podvin, Rivière and Soufiani2019). The ratio in (4.3) gives a rough estimate of the energy transferred from the low-order truncation to the next higher-order scales. This ratio reaches a maximum at $Ra=3\times 10^6$ and is lowest for the highest Rayleigh numbers. The noise level is indicated in Appendix A.

4.2. Time integration

The dynamics of modes $L_{x/y}$,$D$, $C$ and $BL_{x/y}$ have already been discussed in previous work for the case $Ra=10^7$ (Soucasse et al. Reference Soucasse, Podvin, Rivière and Soufiani2019). We checked that the general features of the phase portraits held over the range of Rayleigh numbers considered. In particular, the following points were established.

1. (i) Quasi-steady states constituting the LSC are characterised by roll $L_{x/y}$ and boundary layer modes $BL_{x/y}$ connecting the boundary layer and the core region.

2. (ii) The mode $D$ indicates along which diagonal the LSC takes place and tends to stabilise the LSC.

3. (iii) The corner mode $C$ is essential to reproduce reorientations in the model. The frequency of reorientations also depends on the coupling between the boundary layer modes and the roll modes, as well as on the noise level. A higher intensity of the boundary layer mode makes the LSC less likely to become unstable.

Figure 7 shows histograms of modes $L_x$ and $L_y$ at different Rayleigh numbers obtained with the 11-D model and with the DNS. They have been computed from time series of mode amplitudes $a_{L_x}$ and $a_{L_y}$ and represent the probability density function of the flow in the phase space $L_x$–$L_y$. It can be seen that the model captures the change in dynamics with the Rayleigh number: the planar states are more frequently visited at low Rayleigh number and the stability of the diagonal states increases at high Rayleigh number.

Figure 7. Histograms of modes $L_x$ and $L_y$: DNS (a,c,e,g,i) and model (b,d,f,h,j). Uncoupled case with integration time for the model equal to $\Delta t$ (DNS). Rayleigh numbers are: (a,b) $Ra=10^6$; (c,d) $Ra=3\times 10^6$; (e,f) $Ra=10^7$; (g,h) $Ra=3\times 10^7$; (i,j) $Ra=10^8$.

This is further confirmed by the values of reorientation frequencies. Frequencies of zero crossings for modes $L_x$, $L_y$ and $D$, for the uncoupled DNS and model are reported in table 6. Frequency $f_D$ corresponds to the reorientation frequency $f_r$ defined in § 2.3, as the $D$ mode changes sign at each reorientation. Values for the model reported in the table were obtained for an integration time of 40 000 time units. Uncertainties on model frequencies are approximately four to eight times smaller than uncertainties on DNS frequencies, according to the respective integration times. It can be noted in table 6 that the model captures the decrease of reorientation frequencies with the Rayleigh number and fairly reproduces the frequencies observed in the DNS. Discrepancies between DNS and model results are of the order of the uncertainty on DNS frequencies (see table 2). At the higher Rayleigh numbers, we found that it is not always possible to determine reorientation frequencies in the DNS as the available simulation time is relatively small compared with the characteristic time between reorientations. However, it was possible to obtain values for the model since the integration time could be easily extended.

Table 6. Average frequencies $f_n=1/\tau _n$ in the uncoupled DNS and in the uncoupled model where $\tau _n$ is the average time between zeros of $a_n$ (restricted to times larger than $5 \tau _c$ ). The value $N/A$ is indicated when it was not possible to obtain a value (less than two switches observed at higher Rayleigh numbers). Integration time for the model is 40 000 time units.

Time spectra of roll mode $L_x$ and higher-order mode $D^{*}$ are represented in figure 8 and compared with the DNS at $Ra=3\times 10^6$. A good agreement is observed for frequencies up to $10^{-2}$. However, the model is not able to predict the local increase of the circulation frequency $f_c \thicksim 2\times 10^{-2}$ and the energetic content at higher frequencies is underpredicted, which is not surprising given the small size of the truncation.

Figure 8. The DNS and model spectra at $Ra=3\times 10^6$ for mode $L_x$ (a) and mode $D^*$ (b).

4.3. Effects of the additional modes

We study in this section the effects of the additional modes included in the 11-D truncation: modes $D^*$, $C^*$, $BL_{x/y}^*$, $M^*$ (for $Ra \le 10^7$) and $L_z$ (for $Ra \ge 3 \times 10^7$). Modes $L_z$ and $M^*$, respectively, correspond to a toroidal and poloidal circulation in the cube. These differences correspond to different interactions with other modes: the evolution of mode $L_{x,y}$ is determined by the interaction of $M^*$ with mode $BL_{x,y}^*$, as well as mode $L_z$ with mode $BL_{x,y}$. Moreover, the dynamics of $M^*$ depends on interactions between modes $L$ and $BL$, modes $BL$ and $BL^*$ with the same orientation ($x$ or $y$), as well as interactions of mode $C$ and $C^*$, and $D$ and $D^*$. In contrast, the evolution of $L_z$ depends on cross-interactions of modes $C$ and $D$, as well as $C^*$ and $D^*$, and modes $L$ and $BL$ corresponding to different directions ($x$ and $y$). Phase portraits in figure 9 show that the correlations observed in the DNS are reproduced by the model. More quantitatively, the correlation coefficient between modes $C$ and product $DL_z$ at $Ra = 10^8$ is equal to 0.66 in both the model and the DNS. At $Ra=10^7$, the correlation coefficient between modes $M^*$ and product $CC^*$ is equal to 0.5 in the model and 0.3 in the DNS.

Figure 9. (a,b) Phase portrait of $M^*$ and $CC^*$ in the DNS (a) and in the model (b) at $Ra=10^7$; (c,d) the phase portrait of $C$ and $D L_z$ in the DNS (c) and in the model (d) at $Ra=10^8$.

To get further insights on the effects of the additional modes, we have compared the dynamics of the model using either the 11-D truncation or the 6-D truncation. Generally speaking, the results obtained with the 6-D truncation are very similar to those obtained with the 11-D truncation, which means the additional modes do not alter the general dynamics of the reorientations. Reorientation frequencies given in table 6 are retrieved with the 6-D model within the same accuracy. This small influence of the additional modes is further confirmed by a linear stability analysis performed around diagonal equilibria for both truncations at each Rayleigh number. The resulting eigenvalue spectra are given in figure 10. First, it can be noted that the shape of the spectrum is globally preserved with the Rayleigh number. Comparison of the 6-D and the 11-D truncations shows that the 6-D reproduces the essential features of the 11-D truncation, that is: (i) a pair of very stable eigenvalues associated with an eigenvector with a strong component along the $D$ mode; (ii) two real eigenvalues associated with eigenvectors coupling $L_{x/y}$ and $BL_{x/y}$ modes; (iii) a pair of eigenvalues close to the stability limit associated with an eigenvector with a strong $C$ component. This last point confirms the key role of the corner flows associated with the $C$ mode in the reorientation process, even at low Rayleigh number where the energy of this mode is proportionately weaker. In addition, the 11-D spectrum contains one real eigenvalue, associated with either $L_z$ or $M^*$, and two additional pairs of complex eigenvalues associated with modes $C^*$ and $BL_{x/y}^*$, and $D^*$ and $BL_{x/y}^*$. These two pairs of eigenvalues are close to each other, except at $Ra=10^6$, where the real part of the eigenvalues associated with $D^*$ and modes $BL_{x/y}^*$ is much more stable than at all other Rayleigh numbers.

Figure 10. Linear stability analysis spectrum around the $(L_x,L_y)$ equilibria. Comparison between the 6-D and 11-D truncations of the uncoupled model. Rayleigh numbers are (a) $Ra=10^6$; (b) $Ra=3\times 10^6$; (c) $Ra=10^7$; (d) $Ra=3\times 10^7$; (e) $Ra=10^8$.

5.1. Construction

In the framework of the Boussinesq approximation (weak temperature differences), the temperature dependence of the radiation field can be assumed to be linear. We are therefore able to define a modal-radiative power $\mathcal {P}^{rad}_n(\boldsymbol {r})$, corresponding to the radiative response to a temperature eigenfunction $\phi _n^\theta (\boldsymbol {r})$ and such that $\mathcal {P}^{rad}(\boldsymbol {r},t)=\sum _n a_n(t) \mathcal {P}^{rad}_n(\boldsymbol {r})$. Using this decomposition in the energy balance and applying Galerkin projection, a linear radiation term is obtained and the coupled POD model takes the general form

(5.1)\begin{equation} \frac{\textrm{d}a_n(t)}{\textrm{d}t} = (L_{nm}^B+L_{nm}^D+L_{nm}^R) a_m(t) + Q_{nmp} a_m(t) a_p(t) +T_n(t), \end{equation}

where $L_{nm}^R$ are linear radiation coefficients. Definition and computation details for coefficients $L_{nm}^R$ and for the modal radiative power $\mathcal {P}^{rad}_n(\boldsymbol {r})$ are given in Appendix B and can be also found in Soucasse et al. (Reference Soucasse, Podvin, Rivière and Soufiani2020). We will assume that the form of the closure law (see (4.2)) remains valid in the presence of radiation.

Our goal is to build an a priori model of radiative transfer effects based on uncoupled DNS data and uncoupled POD eigenfunctions. Therefore, model coefficients are determined as follows:

1. (i) Coefficients $L_{nm}^B$, $L_{nm}^D$ and $Q_{nmp}$ are taken from the uncoupled model.

2. (ii) Coefficients $L_{nm}^R$ are added and computed from uncoupled temperature eigenfunctions $\phi ^\theta _n$ and uncoupled factor $\gamma$.

3. (iii) Adjustable coefficients $L_{nm}^A$ are taken from the uncoupled model. This means we assume that the energy transfer from the large scales to the small scales is unchanged.

However, equilibrium values of mode amplitude $a_n^{eq,rad}$ are required in the coupled model to estimate the contribution of the mean mode $M$, to estimate the fluctuating energy $\left \langle k \right \rangle$ in (4.2) and to estimate the noise level. They are not known a priori and need to be predicted. To do so, we use equilibrium relations ($da_n^{eq,rad}/dt=0$) for diagonal states, with the additional radiation terms, the adjustable coefficients being known. Predicted equilibria $\tilde {a}_n^{eq,rad}$ for modes $M$, $L_{x/y}$, $D$ and $BL_{x/y}$ are given in table 7 ($a_n^{eq,rad}=0$ for the other modes) and compared with actual coupled equilibria $a_n^{eq,rad}$ and uncoupled equilibria $a_n^{eq}$ extracted from the DNS. We can see that the effect of radiative coupling is to increase the energy of the equilibria, by lesser amounts as the Rayleigh number increases. Finally, the noise level in the coupled model $\sigma ^{rad}$ is determined from the energy increase of the large-scale modes compared with the uncoupled case

(5.2)\begin{equation} \frac{\sigma^{rad}}{\sigma} = \sqrt{\frac{\sum_n \left(\tilde{a}_n^{eq,rad}\right)^2}{\sum_n \left(a_n^{eq}\right)^2}}. \end{equation}

Noise level values are reported in Appendix A. As expected, the noise level with radiation $\sigma ^{rad}$ is always greater than the noise level without radiation due to the energy increase.

Table 7. Non-zero equilibrium values in the uncoupled ($a_n^{eq}$) and coupled ($a_n^{eq,rad}$) cases (extracted from the DNS) and predicted equilibrium values ($\thicksim$) with the predicted coupled model. Relative ratios compared with the uncoupled case are indicated in parentheses.

5.2. Mode energies

The predicted coupled model was then integrated in time, using an 11-D truncation, leading to predicted amplitudes $a_i^{rad}$. The predicted energy (or predicted eigenvalue) of the modes with radiative coupling $\tilde {\lambda }^{rad}$ was then determined using:

1. (a) the ratio $\eta _n$ determined above for the non-zero equilibrium values $\eta _n=(\tilde {a}_n^{eq,rad}/a_n^{eq})^2$ (with $n$ corresponding to modes $M,L_{x,y},D,BL_{x,y})$;

2. (b) the ratio $\tilde {\eta }_n$ predicted from the model $\tilde {\eta }_n = \left \langle (a_n^{rad})^2 \right \rangle /\left \langle a_n^2 \right \rangle$ for modes which are zero at equilibria (with $n$ corresponding to modes $D^*,C,C^*, BL_{x,y}^*$, $L_z$ or $M_z$).

Figure 11 compares for each mode and each Rayleigh number the predicted energy $\tilde {\lambda }_n^{rad}$ with the eigenvalues $\lambda _n$ and $\lambda _n^{rad}$ found in the uncoupled and coupled DNS. In addition, the values of the total energy contained in the first 12 modes are reported in table 8. Except for the lowest Rayleigh number, corresponding to a large radiative increase which is overpredicted by the model, the total energy increase is correctly predicted by the model (within approximately 10 % maximum). Overall, the model tends to overpredict the amplitude of the $BL_{x/y}$ modes (which are associated with the largest relative increase) and underpredict that of the $L_{x/y}$ modes. However, the agreement is generally good, although the ordering of the modes observed in the coupled DNS is not always captured.

Figure 11. Predicted energies compared with POD eigenvalues in the uncoupled and coupled DNS. Energies are ranked according to the mode ordering in the uncoupled simulation. Rayleigh numbers are (a) $Ra=10^6$; (b) $Ra=3\times 10^6$; (c) $Ra=10^7$; (d) $Ra=3\times 10^7$; (e) $Ra=10^8$.

Table 8. Predicted energy ($\thicksim$) of the first 12 modes compared with eigenvalues in the uncoupled simulation $\lambda _n$ and in the coupled simulation $\lambda ^{rad}_n$. Relative ratios compared with the uncoupled case are indicated in parentheses.

Figure 12 shows histograms of modes $L_x$ and $L_y$ at different Rayleigh numbers obtained with the coupled predicted model and the coupled DNS. At $Ra=10^6$, the radiative coupling makes diagonal states unstable, and the LSC in the simulation tends to become parallel to the cavity sides. This change in dynamics is qualitatively predicted by the model: it can be seen in the histograms of figure 12 that the system spends considerable time near the roll states ($L_x=0$ or $L_y=0$). This will be further confirmed by the linear stability analysis in § 5.4. At $Ra=3\times 10^6$, the coupled system seems to spend more time near the roll states although the diagonal states are still dominant. This is also correctly captured by the predicted model. Changes due to radiative coupling are less important at $Ra \ge 10^7$ in both the model and the DNS. At $Ra \ge 3\times 10^7$, very few reorientations are observed and the amplitude of the oscillations about the equilibrium states decrease.

Figure 12. Histograms of modes $L_x$ and $L_y$: DNS (a,c,e,g,i) and model (b,d,f,h,j). Coupled case with integration time for the model equal to $\Delta t$ (DNS). Rayleigh numbers are (a,b) $Ra=10^6$; (c,d) $Ra=3\times 10^6$; (e,f) $Ra=10^7$; (g,h) $Ra=3\times 10^7$; (i,j) $Ra=10^8$.

5.3. Reorientation frequencies

Table 9 compares frequencies of zero crossings for modes $L_x$, $L_y$ and $D$, for the coupled DNS and the predicted coupled model. Generally speaking, a good agreement is observed between the predicted model and the coupled DNS (discrepancies of 10 %–30 %, of the order of the uncertainty on DNS frequencies), except at $Ra=10^8$ where the reorientation frequency is underestimated. We note that results are not given for the coupled case at $Ra=10^6$ since the diagonal states are no longer observed. We note that at $Ra=3\times 10^6$ large and fast oscillations in the $D$ mode do not make it possible to determine a frequency in both the coupled simulation and the model. As in § 4.2, values for the model reported in the table were obtained for an integration time of 40 000 time units.

Table 9. Average frequencies $f_n=1/\tau _n$ in the coupled DNS and in the coupled model where $\tau _n$ is the average time between zeros of $a_n$ (restricted to times larger than $5 \tau _c$). The value $N/A$ is indicated when it was not possible to obtain a value (less than two switches observed at higher Rayleigh numbers, many rapid switches at $Ra=3\times 10^6$). Integration time for the model is 40 000 time units.

5.4. Model stability at $Ra=10^6$

A remarkable feature of the predicted coupled model is its ability to predict the loss of stability of the diagonal rolls for the benefit of the planar rolls at $Ra = 10^6$. In fact, linear stability analysis of the diagonal equilibria for both 6-D and 11-D truncations (see figure 13) shows that the least stable pair of eigenvalues associated with mode $C$ becomes markedly unstable with radiative coupling. Integration of the models without noise shows that in the presence of radiation the now unstable fixed point destabilises towards another fixed point in the $L_x=0$ or $L_y=0$ space (also in figure 13) corresponding to roll modes. We emphasise that both 6-D and 11-D truncations provide very similar results, indicating that the essential dynamics are captured by the 6-D truncation.

Figure 13. (a) Linear stability analysis spectrum around the $(L_x,L_y)$ diagonal equilibria for the predicted coupled model at $Ra=10^6$. (b) Trajectory in the $(L_x,L_y)$ space of the coupled predicted model without noise from unstable diagonal equilibria: 6-D truncation (black dashed lines) and 11-D truncation (red dotted lines). Symbols indicate stable fixed point for both coupled and uncoupled models.

At other Rayleigh numbers, we found that for both 6-D and 11-D truncations the pair of eigenvectors associated with the corner mode $C$ became less stable when radiation was taken into account, but their eigenvalue remained negative.

6.1. Uncoupled case

The basic assumption in building the generalised model is that the variations of the eigenfunctions can be neglected over the Rayleigh number range. The key ingredient is thus to estimate correctly the distribution of the POD spectrum in order to determine the equilibria, the adjustable coefficients and the noise level. The form of the uncoupled generalised model is the same as that of the uncoupled model ((4.1) together with closure law (4.2)). Model parameters are determined as follows.

1. (i) Linear and quadratic coefficients $L^B_{nm}$, $L^D_{nm}$ and $Q_{nmp}$ are directly determined from a set of POD eigenfunctions $\boldsymbol {\phi }_n=\{\boldsymbol {\phi }_n^u, \gamma \phi _n^\theta \}$. Here, $\boldsymbol {\phi }_n^u$ and $\phi _n^\theta$ are assumed to be constant with the Rayleigh number and taken at $Ra=10^7$. The coefficient $\gamma ^2$ is assumed to vary with the Rayleigh number according to the following law:

   (6.1)\begin{equation} \gamma^2 \thicksim 0.032 Ra^{1/4}, \end{equation}
   which is consistent with the scaling of the turbulent fluctuations in the bulk $\theta ^* \thicksim Ra^{-0.14}$ (Castaing et al. Reference Castaing, Gunaratne, Heslot, Kadanoff, Libchaber, Thomae, Wu, Zaleski and Zanetti1989).

2. (ii) The adjustable coefficients $L^A_{nm}$ are determined from equilibrium relations for a diagonal state. Equilibrium values $a^{eq}_n$ strongly depend on the Rayleigh number and are estimated using the following ansatz for the energy $\lambda _n$ of the large-scale modes $n \in \{M,L_{x/y},D,BL_{x,y},C \}$:

   (6.2)\begin{equation} \lambda_n = h_n \log(Ra) + g_n, \end{equation}
   which is a good approximation at all Rayleigh numbers except the lowest one ($Ra=10^6$), as it can be seen in figure 14. The lowest Rayleigh number was therefore not taken into account when fitting the coefficients $h_n$ and $g_n$. Note that the ratio between the equilibrium value of modes $BL_{x/y}$ and their amplitude was taken to be constant and equal to $\eta _{BL}=0.2$.

3. (iii) As in § 4, the noise level is determined using (4.3), from the energy of the large scales and the energy of the intermediate scales. We have highlighted in § 3.2 that the spectrum decay of the intermediate scales is log–linear and we adopt the following modelling for their energy:

   (6.3)\begin{equation} \lambda_n = C(Ra) n^{\alpha(Ra)} \quad \mathrm{with} \ \left\{ \begin{array}{@{}l} \log(C(Ra))= 1.9-0.5 \log(Ra), \\ \alpha(Ra)={-}2.0+0.07 \log(Ra). \end{array} \right. \end{equation}
   The fit was determined from the variations of the spectrum $\lambda _n$ with $8 \le n \le N/2$, with $N$ such that the $N$ first POD modes capture 95 % of the total energy. We have checked that the asymptotic decay was correctly captured. It can also be noted in figure 14 that the intermediate-scale fit provides a correct estimation of the energy of mode $C$, with a discrepancy between the large-scale and the intermediate-scale fit of 25 % maximum at the lowest Rayleigh number. Noise levels are reported in Appendix A.

Figure 14. Evolution of the energy $\lambda _n$ of the large-scale modes given by the DNS and the large-scale fit (6.2).

The uncoupled generalised model was integrated in time for the different Rayleigh numbers. The measured reorientation frequencies are given in table 10. Overall, the model correctly captures the amplitude of the modes and the associated time scales, except at $Ra=10^6$. At this Rayleigh number, the energy of the large-scale modes is underestimated by the large-scale fit, which leads to an overestimation of the noise level and an overestimation of the reorientation frequency. At $Ra=10^8$, no reorientation event has been detected by the model despite a large integration time.

Table 10. Reorientation frequencies obtained with the uncoupled generalised model. Integration time for the model is 40 000 time units.

6.2. Coupled case

Following the procedure described in § 5.1, radiation effects are incorporated in the generalised model. Model parameters for the coupled generalised model are thus determined as follows.

1. (i) Coefficients $L_{nm}^B$, $L_{nm}^D$ and $Q_{nmp}$ are taken from the uncoupled generalised model.

2. (ii) Coefficients $L_{nm}^R$ are added and computed from uncoupled temperature eigenfunctions $\phi ^\theta _n$ at $Ra=10^7$ and $\gamma$ factor modelled by (6.1).

3. (iii) Adjustable coefficients $L_{nm}^A$ are taken from the uncoupled generalised model. This allows us to estimate equilibrium values $\tilde {a}_n^{eq,rad}$ and energies $\tilde {\lambda }_n^{rad}$ for the modes in the truncation (except mode $C$ which is zero at equilibrium).

4. (iv) The noise level is determined using (5.2).

Differences between all the models (uncoupled model, predicted coupled model, uncoupled generalised model and coupled generalised model) are highlighted in Appendix A.

The total energy of the large-scale modes $\sum _{n=1}^7\tilde {\lambda }_n^{rad}$ predicted by the coupled generalised model is represented in figure 15 and compared with values given by the coupled DNS and the coupled predicted model. In addition, the energy increase for each large-scale mode $M, L_{x,y}, D, BL_{x,y}$ is given in figure 16. In terms of total large-scale energy, a very good agreement is observed for the generalised model (discrepancies approximately 5 %), even at $Ra=10^6$. In this case, as figure 16 shows, it is due to the compensation of estimation errors ($M$ is overpredicted by 25 %, while modes $L$ are underpredicted by 30 %). At other Rayleigh numbers, the evolution of the energy of the different modes is relatively well predicted (between 5 % and 20 %), except for an underprediction of modes $L$. Overall a good agreement is observed between the different models.

Figure 15. Total energy of the large-scale modes $\sum _{n=1}^7$ as a function of the Rayleigh number – comparison of the DNS, the coupled predicted model and the coupled generalised model.

Figure 16. Relative increase of the amplitude of the large-scale modes due to radiation with the Rayleigh number – comparison of the DNS, the coupled predicted model and the coupled generalised model.

The measured reorientation frequencies obtained with the coupled generalised model are given in table 11 and compared with the coupled DNS. Generally speaking, the model fairly reproduces the coupled DNS results (discrepancies of 10 %–30 %) and predicts the increase of the reorientation frequency with radiation, which tends to subside with the Rayleigh number. However, it should be pointed that the change in dynamics observed at $Ra=10^6$ (loss of stability of the diagonal rolls) is not predicted with the generalised coupled model, which is consistent with the lack of validity of the modelling assumptions at this Rayleigh number.

Table 11. Reorientation frequencies obtained with the coupled generalised model. Integration time for the model is 40 000 time units.