Asymmetry of Two-Dimensional Thermal Convection at High Rayleigh Numbers

Abstract

While thermal convection cells exhibit left–right and top–bottom symmetries at low Rayleigh numbers ($Ra$), the emergence of coherent flow structures, such as elliptical large-scale circulation in Rayleigh–Bénard convection (RBC), breaks these symmetries as the Rayleigh number increases. Recently, spatial double-reflection symmetry was proposed and verified for two-dimensional RBC at a Prandtl number of 6.5 and $Ra$ values up to ${10}^{10}$. In this study, we examined this new symmetry at a lower Prandtl number of 0.7 and across a wider range of Rayleigh numbers, from ${10}^7$ to ${10}^{13}$. Our findings reveal that the double-reflection symmetry is preserved for the mean profiles and flow fields of velocity and temperature for $Ra<{10}^9$, but it is broken at higher Rayleigh numbers. This asymmetry at high $Ra$ values is inferred to be induced by a flow-pattern transition at $Ra={10}^9$. Together with the previous study, our results demonstrate that the Prandtl number has an important influence on the symmetry preservation in RBC.

1. Introduction

Turbulence, which is prevalent in both nature and industrial systems [1,2], is a complex phenomenon in fluids that is characterized by chaotic changes in velocity and pressure. One example of these complex phenomena is the Rayleigh–Bénard convection (RBC). This phenomenon occurs in a system that is heated from the bottom and cooled from the top with two constant temperatures ($T_{\mathrm{bot}}$ and $T_{\mathrm{top}}$), while the sidewalls remain adiabatic. The fluid in this system is driven by buoyancy and gravity. The cold fluid sinks while the hot fluid rises. They encounter and mix in the center of the cell, leading to complex flow patterns. The geometry of such a system can be arbitrary, including cubic, rectangular, or cylindrical configurations. The characteristics of RBC are determined by three key parameters: the Rayleigh number $Ra={\displaystyle \frac{\alpha g\Delta T{H}^3}{\nu \kappa }}$, the Prandtl number $Pr={\displaystyle \frac{\nu }{\kappa }}$, and the aspect ratio $\Gamma ={\displaystyle \frac{D}{H}}$. Here, $\alpha$ is the thermal-expansion coefficient, $\nu$ is the viscosity, $\kappa$ is the thermal-diffusion coefficient, g is the gravitational acceleration, H is the system’s height, and D is the system’s width or diameter. The temperature difference between the horizontal conducting plates at the bottom and top of the system is represented by $\Delta T=T_{\mathrm{bot}}-T_{\mathrm{top}}$.

Over the years, research on RBC has encompassed various aspects, including the system’s heat transfer, the emergence of coherent structures, and fluid dynamics. In early studies, the instability and chaos in RBC systems attracted much attention [3,4,5,6]. After the moderately high-$Ra$ experiments reported by Castaing et al. in 1989 [7], the properties of turbulence in RBC were widely investigated by researchers, including considerations of velocity and thermal fluctuations [8,9], kinetic and thermal dissipation [10,11,12], boundary layers [13,14], and small-scale turbulence [15,16]. Turbulent RBC shows different flow patterns with different $Ra$–$Pr$ phases [17,18], with some coherent structures emerging, i.e., large-scale circulation (LSC) and plumes [12]. In particular, the flow patterns experience a transition with increasing $Ra$, leading to different statistical behaviors [18]. These properties are closely related to the heat transfer of RBC. For instance, thermal plumes are transported upward or downward by LSC, facilitating the transfer of heat flux; however, alterations to the boundary conditions or constraints on the flow motion can lead to significant changes in the system’s heat transfer behavior [19,20,21]. Furthermore, heat transfer is influenced by the cell’s aspect ratio and shape, both of which have garnered significant attention [22,23,24,25,26]. The heat transfer is quantified by the Nusselt number ($Nu$), and its scaling with $Ra$ and $Pr$ has been discussed in numerous reports [27,28,29,30,31,32,33]; however, no unified conclusions have yet been reached on this topic.

Symmetry in RBC systems has also been investigated in many studies [34,35,36,37]. Symmetry is an important concept in physics as it reflects the invariance within a system [38], and it is associated with a pattern that satisfies invariance under a specific transformation rule. Generally speaking, if there is an invariant quantity—one that remains unchanged under transformation—then a corresponding symmetry exists. With symmetric boundary conditions and symmetric geometry, symmetry is expected in the properties of the flow within an RBC system. At very low $Ra$ values, an RBC system will maintain a motionless state with a uniform vertical temperature gradient. Above a critical Rayleigh number $R{a}_{\mathrm{c}}$, a convective motion appears and swirl structures are formed [35]. Then, as the $Ra$ is increased, the flow state of the RBC experiences convection, oscillation, chaos, transition, and turbulence [39]. In these states, the symmetry of an RBC system also experiences different states, such as axisymmetric convection and non-axisymmetric motions [35]. With cryogenic helium gas as the working fluid, for moderately high $Ra$ values in the turbulent state, it was found that the flow is statistically symmetric for $Ra>{10}^{11}$ in 3D cylinders with aspect ratios ($\Gamma$) of 0.5 and 1.0 [34]. In experiments with water ($Pr=5.3$), the mean velocity field is also azimuthally symmetric in a cylinder with $\Gamma =0.5$, while no such symmetry is observed in a cylinder with $\Gamma =1.0$ [36,37]. For a large $Pr$, the local boundary layers have been found to be either symmetric [40] or asymmetric [41,42,43]. In a 2D RBC, elliptical LSC forms within the cell at moderately high $Ra$ values [11,44], as shown in Figure 1; therefore, the left–right and top–bottom symmetry of the system will be broken. These results indicate that the symmetry might be broken or restored in different cells or those with different $Ra$–$Pr$ phases; in other words, symmetry plays a crucial role in the potential transitions exhibited by this system. Recently, the x–y double-reflection symmetry in thermal convection at $Re$ values between ${10}^7$ and ${10}^{10}$ has been proposed, and the symmetry under this reflection was validated [45]. This symmetry for other $Pr$ values or higher $Ra$ values is unknown. For this reason, we sought to examine double-reflection symmetry for higher $Ra$ values at $Pr=0.7$.

In this work, the x–y double-reflection symmetry in a 2D RBC with $Ra$ values ranging from ${10}^7$ to ${10}^{13}$ and $Pr=0.7$ was investigated using direct numerical simulation (DNS). In the simulation model, the aspect ratio of the cell was unity, and it was confined by adiabatic walls. The $Ra$ range in the previous study was only up to ${10}^{10}$, and it was necessary to verify the symmetry at higher $Ra$. In addition, the $Pr$ effect in the double reflection symmetry was unknown. Therefore, the $Pr$ in this paper was set as 0.7, which is different from the $Pr=6.5$ set in the previous study. The remainder of this paper is structured as follows. The numerical setup of the DNS is described in Section 2, the temperature and velocity flow fields are presented in Section 3, the verification of the double-reflection symmetry using the data is discussed in Section 4, and the conclusions are presented in Section 5.

The RBC system is considered to be incompressible, with the temperature treated as an active scalar. In general, the Oberbeck–Boussinesq approximation is used to simplify the Navier–Stokes equations [14,18,46], i.e., $\rho \left(T\right)=\rho \left(T_0\right)\left[1-\alpha (T-T_0)\right]$, while various transport coefficients (e.g., $\alpha$, $\nu$, and $\kappa$) are treated as constants. Under this approximation, the effect of temperature on the velocity field is realized through the buoyancy term. Using the scaling factors H, $\Delta T$, $U=\sqrt{\alpha gH\Delta T}$ and $\tau =\sqrt{H/\alpha g\Delta T}$ [14], the governing equations can be written in a dimensionless form as follows:

$$\begin{array}{ccc}\hfill \nabla ·\mathbf{u}& =& 0\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\left(\mathbf{u}·\nabla \right)\mathbf{u}& =& -\nabla p+{\displaystyle \frac{1}{\sqrt{Ra/Pr}}}{\nabla }^2\mathbf{u}+\theta \widehat{y},\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \theta }{\partial t}}+\left(\mathbf{u}·\nabla \right)\theta & =& {\displaystyle \frac{1}{\sqrt{Ra·Pr}}}{\nabla }^2\theta ,\hfill \end{array}$$

(3)

where $\mathbf{u}=(u,v)$ is the dimensionless velocity, p is the dimensionless pressure, $\theta$ is the dimensionless temperature, and $\widehat{y}$ is a vertical unit vector whose direction is opposite to that of gravity. For the thermal boundary conditions, the temperatures of the two conducting plates were kept constant, while the two sidewalls were adiabatic:

$$\begin{array}{c}\hfill \theta (y=0)=0.5,\theta (y=1)=-0.5\end{array}$$

(4)

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \theta }{\partial x}}\left(x=0\right)=0,{\displaystyle \frac{\partial \theta }{\partial x}}\left(x=1\right)=0.\end{array}$$

(5)

No-slip velocity boundary conditions were imposed on all the conducting plates and sidewalls:

$$\begin{array}{c}\hfill u(y=0)=0,u(y=1)=0,u(x=0)=0,u(x=1)=0\end{array}$$

(6)

$$\begin{array}{c}\hfill v(y=0)=0,v(y=1)=0,v(x=0)=0,v(x=1)=0.\end{array}$$

(7)

2. Equations and Numerical Settings

The simulation scheme applied in this work uses a finite-difference method known as the parallel direct method of DNS [47]. This process for solving the Navier–Stokes equations are illustrated in Figure 2. Specifically, the process of solving the Poisson equation is also shown there. In the flowchart, $f({\mathbf{u}}^n,{\theta }^n)=-\left({\mathbf{u}}^n·\nabla \right){\mathbf{u}}^n+{\displaystyle \frac{1}{\sqrt{Ra/Pr}}}{\nabla }^2{\mathbf{u}}^n+{\theta }^n\widehat{y}$ and $g({\mathbf{u}}^{n+1},{\theta }^n)=-\left({\mathbf{u}}^{n+1}·\nabla \right){\theta }^n+{\displaystyle \frac{1}{\sqrt{Ra·Pr}}}{\nabla }^2{\theta }^n$. A fully explicit projection method, with second-order accuracy in both space and time, was used to solve the governing equations. In the first step, the pressure $p^{n+1}$ was unknown, and the Poisson equation for the pressure was solved in the third step. An explicit second order Runge–Kutta scheme was employed in the time-marching. The Poisson equation for pressure was solved using the parallel diagonal dominant (PDD) algorithm in combination with fast Fourier transforms (FFTs). Finally, the velocity ${\mathbf{u}}^{n+1}$ was corrected with $u^{*}$ and pressure $p^{n+1}$, and the temperature ${\theta }^{n+1}$ was calculated with ${\mathbf{u}}^{n+1}$ and ${\theta }^n$.

The code was written in Fortran using parallel technologies, including MPI and OpenMP. We used this code for the work we reported in previous publications that related to investigating flow patterns, the scaling transition of thermal dissipation, and boundary layers [14]. In our simulations, staggered grids were employed. The grids were uniform in the horizontal direction (x), while they were un-uniform in the vertical direction (y). To accurately capture the intense fluctuations near the top and bottom plates, the grid points were clustered closer to these regions.

The key parameters and grid for the cases in this work are outlined in

Table 1. Details of the numerical simulations employed for the present simulations. Here, $Nu$ is calculated as $Nu=\sqrt{RaPr}{〈v\theta 〉}_{\Omega ,t}-{〈{\partial }_y\theta 〉}_{\Omega ,t}$, where ${〈〉}_{\Omega ,t}$ represents the time–space average over the entire system; $N_x$ and $N_y$ represent the grid resolution in the horizontal and vertical directions, respectively; and $N_{\mathrm{BL}}$ is the number of grid points in the thermal boundary layer of a conducting plate.

$\mathit{Ra}$	$\mathit{Pr}$	$\Gamma$	${\mathit{N}}_{\mathit{x}}\times {\mathit{N}}_{\mathit{y}}$	${\mathit{t}}_{\mathbf{avg}}$	$\mathit{Nu}$	${\mathit{N}}_{\mathbf{BL}}$
$1\times {10}^7$	$0.70$	$1.0$	$512\times 576$	1000	$11.38$	33
$2\times {10}^7$	$0.70$	$1.0$	$512\times 576$	1000	$14.33$	27
$5\times {10}^7$	$0.70$	$1.0$	$512\times 576$	1000	$19.71$	20
$1\times {10}^8$	$0.70$	$1.0$	$512\times 576$	1000	$25.22$	16
$2\times {10}^8$	$0.70$	$1.0$	$512\times 576$	1000	$31.24$	13
$5\times {10}^8$	$0.70$	$1.0$	$1024\times 1152$	1000	$42.13$	26
$1\times {10}^9$	$0.70$	$1.0$	$1024\times 1152$	400	$50.53$	21
$2\times {10}^9$	$0.70$	$1.0$	$1024\times 1152$	400	$59.83$	18
$5\times {10}^9$	$0.70$	$1.0$	$1024\times 1152$	400	$78.25$	14
$1\times {10}^{10}$	$0.70$	$1.0$	$1536\times 1728$	600	$95.50$	30
$2\times {10}^{10}$	$0.70$	$1.0$	$2048\times 2304$	400	$118.23$	36
$5\times {10}^{10}$	$0.70$	$1.0$	$2048\times 2304$	400	$158.12$	28
$1\times {10}^{11}$	$0.70$	$1.0$	$2048\times 2304$	400	$198.43$	22
$2\times {10}^{11}$	$0.70$	$1.0$	$2560\times 2880$	200	$243.14$	40
$5\times {10}^{11}$	$0.70$	$1.0$	$2560\times 2880$	200	$326.14$	30
$1\times {10}^{12}$	$0.70$	$1.0$	$4096\times 4608$	200	$404.86$	51
$2\times {10}^{12}$	$0.70$	$1.0$	$4096\times 4608$	300	$502.11$	41
$5\times {10}^{12}$	$0.70$	$1.0$	$5120\times 5760$	200	$666.74$	57
$1\times {10}^{13}$	$0.70$	$1.0$	$5120\times 5760$	50	$839.71$	48

. To resolve the smallest turbulence scales—i.e., the Kolmogorov scale ${\eta }_K=HP{r}^{1/2}/{\left[Ra(Nu-1)\right]}^{1/4}$ and the Batchelor scale ${\eta }_B={\eta }_{K}P{r}^{-1/2}$—the grid in the vertical direction (y) was refined in the boundary layer, and the number of grid points ($N_{\mathrm{BL}}$ in Table 1) in the thermal boundary layer satisfied the condition derived by Shishkina et al. [48]:

$$\begin{array}{ccc}\hfill N_{\theta ,\mathrm{BL}}& \ge & \sqrt{2}aNuP{r}^{-0.5355+0.033logPr},\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill N_{v,\mathrm{BL}}& \ge & \sqrt{2}aN{u}^{1/2}P{r}^{-0.1785+0.011logPr},\hfill \end{array}$$

(9)

where $a\approx 0.482$ and $N_{\theta ,\mathrm{BL}}$ and $N_{v,\mathrm{BL}}$ are the minimum numbers of grid points in the thermal boundary layer and the viscous boundary layer, respectively. Note that Equations (8) and (9) were only valid for $3\times {10}^{-4}\le Pr\le 1$; the conditions for other $Pr$ values can be found in Ref. [48] and are not discussed here. The dimensionless time step $\Delta t$ was set to be less than $1/1000$ of the Kolmogorov time scale ${\tau }_K=\sqrt{Pr/(Nu-1)}$ to accurately capture the intense fluctuations. The simulations were executed on the Tianhe-2 supercomputer, requiring millions of core hours to complete.

For cases with low $Ra$ values (≤$2\times {10}^8$), the number of grid points was relatively small, allowing the statistical time $t_{\mathrm{avg}}$ to reach up to 1000 free-fall times ($t_{\mathrm{f}}$), with $t_{\mathrm{f}}$ calculated as $t_{\mathrm{f}}=\sqrt{H/\alpha g\Delta T}$. For higher $Ra$ values, the number of grid points rapidly increased, and the cost of the simulation also increased, resulting in a shorter statistical time. The value of $t_{\mathrm{avg}}$ was set to at least $200t_{\mathrm{f}}$ before $Ra={10}^{13}$ to ensure convergence of the data. It was noted that $Nu$ is a key response variable for assessing the state of the system and is commonly used to verify convergence. For the statistical analysis, $Nu$ was calculated as the average of data from two successive long-term statistical segments. It was found that the relative differences between the $Nu$ values from these segments were within 1%. For the highest $Ra$ case, $t_{\mathrm{avg}}$ was only $50t_{\mathrm{f}}$, but the relative difference was about 3.6%.

3. Flow Fields

3.1. Instantaneous Flow Fields

The instantaneous velocity and temperature fields are presented in Figure 3. In the top row of Figure 3, the high-speed fluid is represented in red and the low-speed fluid is represented in blue. In the bottom row of Figure 3, red is used to represent the high-temperature fluid, while blue represents the low-temperature fluid. Note that, in the 2D RBC, the direction of the LSC was arbitrary as it could be either clockwise or counter-clockwise. Similar to the approach taken in Ref. [13], the direction of the LSC was unified to be counter-clockwise in this work, and this did not affect the analysis.

In cases with low $Ra$ values, a large elliptical swirl (red on the outside and blue on the inside) fills almost the entire system, and there are usually two small swirls beside it (see Figure 3a,b top-right and bottom-left corners). Generally, this large elliptical swirl is referred to as LSC, and the small counter-rotating swirls are referred to as corner rolls. In contrast, the top-left and bottom-right corners were almost blue, meaning that the motion of the fluid in these regions was very slow. The four corners were filled by cold and hot fluid, as shown in Figure 3e,f, as a result of the emission of thermal plumes. The motion of the thermal plumes was primarily constrained by the LSC, resulting in an elliptical distribution of the mean-temperature fluid in the bulk region. This flow pattern has also been observed at lower $Ra$ values in some previous works [11,44]. Due to the presence of the LSC and corner rolls, the flow pattern remained almost steady for the low $Ra$ values in the current range up to $5\times {10}^8$.

For relatively high $Ra$ values (≥${10}^9$ in this paper), the red areas were less stable, and a large-scale swirl moved around in the system, as shown in Figure 3c,d. Additionally, the corner rolls became less prominent and were able to move away from their original positions. More small swirls (shown by vectors) appeared in the system, making the flow pattern more complex. These swirls varied in scale, which is related to the energy inverse cascade in 2D turbulence [49]. The small swirls were generated by the shearing action between the large-scale swirl and the boundary layers. As these swirls moved throughout the system with the large-scale swirl, they mixed with the larger ones. Meanwhile, the motion of the plumes was disorganized, as shown in Figure 3g,h, and the plumes were mainly carried by the small swirls. Moreover, the size of the plumes decreased as the $Ra$ increased, as has been previously reported [14]. The flow pattern was markedly different from that observed at lower $Ra$ values, as can be seen in Figure 3a,b,e,f. This transition in a flow pattern was also reported by Gao et al., who conducted a more detailed investigation using stability analysis [18]. This transition also lead to changes in the thermal dissipation [14] and fluctuations [9].

In these instantaneous fields, the top–bottom and left–right symmetry of the system was, evidently, broken. Next, we examined the time-averaged temperature and velocity fields.

3.2. Time-Averaged Flow Fields

The time-averaged velocity and temperature fields are shown in Figure 4. As shown in Figure 4a, there were only two corner rolls, in the top-right and bottom-left corners, while there were no swirls in the other two corners. The high-speed flow formed a large swirl, which is the LSC mentioned above. The maximum velocity was located in the region between the LSC and corner roll. As shown in Figure 4e, the time-averaged temperature field revealed a yellow ellipse, indicating a mean temperature distribution of zero in most parts of the cell. The thermal plumes were primarily carried by the LSC or were confined in the corner rolls, leading to the high- and low-temperature regions that correspond to the LSC shown in Figure 4a.

For the $Ra={10}^8$ in Figure 4b,f, the flow pattern and temperature fields were very similar to those for $Ra={10}^7$; however, the flow motion in top-left and bottom-right corners became more complex, suggesting the appearance of additional swirls in these regions. Moreover, the velocity magnitude remained small, indicating that these swirls had little influence on the heat transfer. For these two cases, the time-averaged velocity fields were similar to the instantaneous velocity fields, suggesting that the flow motion was relatively steady for $Ra$ values ranging from ${10}^7$ to $5\times {10}^8$.

As shown in Figure 4c,d, a large red annulus almost filled the system together with four low-speed corner rolls. At $Ra={10}^{10}$, the shapes of these corner rolls appeared different, and their sizes were smaller. This distribution was quite different from that of the instantaneous velocity fields that are shown in Figure 3. The time-averaged velocity field indicated that a large-scale swirl persisted in the cell while many smaller swirls were present, making it challenging to identify the largest swirl. Reviewing the fields in Figure 3c,d, the red region resembles an annulus and it was surrounded by several swirls. Therefore, the LSC still exists in the system at high $Ra$ values, but it was difficult to identify. For the temperature fields in Figure 4g,h, the fluctuations of the plumes were eliminated in the bulk region. The high- (or low-) temperature region in the bottom-left (or top-right) corner disappeared as the flow in these two corner rolls was not steady, unlike those in Figure 4a,b. The hot and cold layers near the conducting plates were thermal boundary layers. In previous studies [14], the boundary layer was divided into three regions (from left to right): the impacting region, the shearing region, and the ejecting region. In these three regions, the motions of the fluid were different, and their properties changed with different $Ra$ values.

The mean fields in Figure 4 indicate that the top–bottom and left–right symmetries were absent. The flow fields for $Ra\le 5\times {10}^8$ were similar to the results for $Pr=6.5$, which were also found in a recent study [45], although the effects of the flow transitions were not included. Therefore, we applied double reflection to the flow fields, as shown in the next section, to investigate the symmetry at higher $Ra$ values.

4. Double-Reflection Symmetry and Data Verification

In this section, we will examine the symmetry by applying an x–y reflection at the center of the system. If a reflection is performed over $x=1/2$ and is followed by another reflection over $y=1/2$, then

$$\begin{array}{c}\hfill t^{*}=t,x^{*}=1-x,y^{*}=1-y,{\mathbf{u}}^{*}=-\mathbf{u},p^{*}=p,{\theta }^{*}=-\theta ,\end{array}$$

(10)

where the superscript * indicates transformed variables. Equations (1)–(3) can then be written as follows:

$$\begin{array}{ccc}\hfill {\nabla }^{*}·{\mathbf{u}}^{*}& =& 0\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial {\mathbf{u}}^{*}}{\partial t^{*}}}+\left({\mathbf{u}}^{*}·\nabla \right){\mathbf{u}}^{*}& =& -\nabla p^{*}+{\displaystyle \frac{1}{\sqrt{Ra/Pr}}}{{\nabla }^{*}}^2{\mathbf{u}}^{*}+{\theta }^{*}{\widehat{y}}^{*}\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial {\theta }^{*}}{\partial t^{*}}}+\left({\mathbf{u}}^{*}·\nabla \right){\theta }^{*}& =& {\displaystyle \frac{1}{\sqrt{Ra·Pr}}}{\nabla }^2{\theta }^{*}.\hfill \end{array}$$

(13)

Comparing Equations (11)–(13) with Equations (1)–(3), it can be seen that their forms are the same; it is, thus, clear that the governing equations are symmetric with the reflections applied in both the x and y directions.

4.1. Flow Fields

As shown in this section, the time-averaged temperature ${\langle \theta (x,y)\rangle }_t$ and velocity magnitude fields $\mathbb{V}=\sqrt{{\langle u\rangle }_t^2+{\langle v\rangle }_t^2}$ were transformed using double reflection, where ${\langle \rangle }_t$ indicates time averaging. For brevity, we use $\Theta (x,y)$ to denote ${\langle \theta (x,y)\rangle }_t$; $U(x,y)$ to denote ${\langle u(x,y)\rangle }_t$; and $V(x,y)$ to denote ${\langle v(x,y)\rangle }_t^2$. Taking the temperature field $\Theta (x,y)$ as an example, the process of the double-reflection transformation is as follows:

$$\begin{array}{c}\hfill \Theta (x,y)\to \Theta (1-x,y)\to -\Theta (1-x,1-y)={\Theta }^{*}(x^{*},y^{*}).\end{array}$$

(14)

This process at $Ra={10}^8$ is directly shown in Figure 5. The process was similar for other physical fields, as shown in Equation (10). From Figure 5a,c, it can be seen that the differences in the temperature fields before and after the transformation were not visually apparent; the two plots are almost identical, showing a dominant LSC accompanied by two corner rolls. Therefore, the relative difference between the two fields was defined to quantify the quality of the double reflection as follows:

$$\begin{array}{c}\hfill D(\Theta )=\left|\Theta (x,y)-{\Theta }^{*}(x^{*},y^{*})\right|/\Delta \Theta \times 100\%,\end{array}$$

(15)

where $\Delta \Theta ={\Theta }_{\mathrm{bot}}-{\Theta }_{\mathrm{top}}=1$. The relative difference between the two temperature fields is shown clearly in Figure 5d. If the difference $D(\Theta )$ across most of the region is very small, the double-reflection symmetry is validated. The symmetry of the velocity magnitude field $\mathbb{V}$ will also be examined in a similar way, with the relative difference $D\left(\mathbb{V}\right)$ defined as follows:

$$\begin{array}{c}\hfill D\left(\mathbb{V}\right)=\left|\mathbb{V}-{\mathbb{V}}^{*}\right|/{\mathbb{V}}_{\max }\times 100\%,\end{array}$$

(16)

where ${\mathbb{V}}^{*}=\sqrt{{U^{*}}^2+{V^{*}}^2}$ and ${\mathbb{V}}_{\max }$ is the maximum velocity magnitude.

For $Ra<{10}^9$, the mean velocity fields were symmetrical after the double reflection, with the relative difference $D\left(\mathbb{V}\right)$ remaining below $1\%$, as shown in Figure 6a,b. The maximum value of $D\left(\mathbb{V}\right)$ was found in the region between the LSC and the corner rolls, where the velocity fluctuations were the most intense [50].

In contrast, for $Ra\ge {10}^9$, the distribution of the difference $D\left(\mathbb{V}\right)$ changed, as shown in Figure 6c,d. The value of $D\left(\mathbb{V}\right)$ in the center of the cell was much larger, exceeding $10\%$. The value of $D\left(\mathbb{V}\right)$ also increased in the middle of the conducting plates (the shearing regions), as well as in the bottom-left and top-right regions, which were dominated by impacting plumes. These distributions indicate that the double-reflection symmetry was broken in the mean velocity fields for $Ra\le {10}^9$ at $Pr=0.7$. This asymmetry corresponds to the flow-pattern transition mentioned above and reflects the change in the flow fields, as shown in Figure 3 and Figure 4. In the flow state for $Ra\le {10}^9$, more of the plumes moved around the center of the cell, inducing more intense fluctuations. Therefore, $D\left(\mathbb{V}\right)$ became larger in the center of the cell.

The mean temperature fields exhibited symmetry after a double reflection for $Ra<{10}^9$, with $D(\Theta )$ values remaining below 1% in Figure 6e,f; however, for $Ra\ge {10}^9$, the $D(\Theta )$ values increased throughout the system, which is similar to the behavior seen in $D\left(\mathbb{V}\right)$. The maximum $D(\Theta )$ value was around $5\%$, as shown in Figure 6g,h, and this was located in the ejecting region. In this region, thermal plumes were primarily emitted, leading to stronger thermal fluctuations [50]. For the mean temperature fields, the symmetry of double reflection held for $Ra<{10}^9$ within the current range, while the symmetry was broken for $Ra\ge {10}^9$. According to [18], the flow pattern transition occurred at about $Ra=1.1\times {10}^{9}P{r}^{1.41}$. In our previous study for $Pr=6.5$, the highest $Ra$ was ${10}^{10}$, while the transition $Ra$ for this $Pr$ was, approximately, $1.54\times {10}^{10}$. Under the transition $Ra$, the double-reflection symmetry held. In this paper, the highest $Ra$ was ${10}^{13}$, which was significantly higher than the transition $Ra$ for $Pr=0.7$ ($6.65\times {10}^8$ approximately). As a result, asymmetry was observed in this regime.

4.2. Mean Velocity and Temperature Profiles

In this section, the symmetry at the low $Ra$ values was examined in more detail using the horizontal profiles of the mean temperature ($\Theta$) and velocity components ($U,V$). Two horizontal sections at different heights, $y=0.01$ and $y=0.1$, were selected to study the double-reflection symmetry of the velocity and temperature profiles. The $y=0.01$ section was closer to the boundary layer, while the $y=0.1$ section was situated in the bulk region.

The mean profile at $y=0.99$ was transformed to correspond with the mean profile at $y=0.01$ in accordance with Equation (10). If the two profiles collapsed onto one another, the symmetry of the double reflection held; otherwise, the symmetry was considered broken. Similarly, the profile at $y=0.9$ was transformed to match the profile at $y=0.1$.

The profiles for $Ra={10}^7$ and ${10}^8$ were compared, as shown in Figure 7. It can be seen that the profile $U(x,y=0.01)$ closely matches $-U(x,y=0.99)$, confirming the symmetry after the double reflection. Similarly, the profiles of $V(x,y=0.01)$ and $\Theta (x,y=0.01)$ also aligned with their counterparts following double reflection (see Figure 7, left-hand panel). Furthermore, the profiles at $y=0.1$ and $y=0.9$ were nearly identical, validating the symmetry of the double reflection. Although the profiles for other $Ra$ values are not shown here and they varied with $Ra$, the double-reflection symmetry consistently held for all $Ra$ values below ${10}^9$ in this study.

It was also evident that these profiles exhibited asymmetry in the x direction. Using the profiles at $y=0.1$ as an example, we can explain the source of the left–right asymmetry. On the left-hand side, the corner roll generated a negative horizontal velocity U and a positive vertical velocity V (see Figure 4a,b). In contrast, on the right-hand side, the absence of a corner roll resulted in velocity components that differed from those on the left side of the cell. The plateaus observed in Figure 7b,f were induced by the shear from the elliptical LSC, and they were not symmetrical about the center axis at $x=0.5$. Therefore, the left–right symmetry was broken in these cases.

At high $Ra$ values (≥${10}^9$), the double-reflection symmetry was broken, as indicated in Figure 8. For the velocity components U and V, as well as the temperature $\Theta$, the trends of the profiles near the bottom plate ($y=0.01,0.1$) were similar to those near the top plate ($y=0.99,0.9$). However, the profiles near the bottom plate (represented by plus symbols) did not coincide with those near the top plate (shown by solid lines). We can see clear discrepancies in Figure 8, and these discrepancies increased with increasing $Ra$.

4.3. Boundary Layer Characteristics

As shown in this section, the wall units of the boundary layers were compared to examine the double-reflection symmetry. The wall units were calculated for two regions (as shown in Figure 1): the top-left region ($y=1$, $0\le x\le 0.5$) and the bottom-right region ($y=0$, $0.5\le x\le 1$).

Following Ref. [46], the friction velocity in the top-left (TL) region $u_{\tau }^T$ or in the bottom-right (BR) region $u_{\tau }^B$ were calculated as

$$u_{\tau }^T={\left({\displaystyle \frac{Pr}{Ra}}\right)}^{1/4}{〈{\left|{\displaystyle \frac{\partial u}{\partial y}}\right|}_{y=1}^{1/2}〉}_{t,0\le x\le 0.5},u_{\tau }^B={\left({\displaystyle \frac{Pr}{Ra}}\right)}^{1/4}{〈{\left|{\displaystyle \frac{\partial u}{\partial y}}\right|}_{y=0}^{1/2}〉}_{t,0.5\le x\le 1},$$

(17)

where ${〈〉}_{t,a\le x\le b}$ denotes the time–space average over the region $a\le x\le b$. The wall units of temperature ${\theta }_{\tau }^T$ and ${\theta }_{\tau }^B$ were also investigated, and these were calculated as follows [46]:

$${\theta }_{\tau }^T=-{\displaystyle \frac{1}{u_{\tau }\sqrt{Ra·Pr}}}{〈{\displaystyle \frac{\partial \theta }{\partial y}}{|}_{y=1}〉}_{t,0\le x\le 0.5},{\theta }_{\tau }^B=-{\displaystyle \frac{1}{u_{\tau }\sqrt{Ra·Pr}}}{〈{\displaystyle \frac{\partial \theta }{\partial y}}{|}_{y=0}〉}_{t,0.5\le x\le 1}.$$

(18)

As shown in Figure 9a,b, the wall units $u_{\tau }$ and ${\theta }_{\tau }$ in the TL and BR regions were identical before $Ra={10}^9$; however, after $Ra={10}^9$, the $u_{\tau }$ in the TL region decreased with increasing $Ra$, while the $u_{\tau }$ in the BR region was slightly smaller than that in the TL. Additionally, the temperature wall unit ${\theta }_{\tau }$ showed a clear discrepancy between the TL and BR regions. Since ${\theta }_{\tau }$ was calculated using $u_{\tau }$ in Equation (18), a decrease in $u_{\tau }$ led to an increase in ${\theta }_{\tau }$.

We further investigated the boundary layer thickness in the TL and BR regions. The boundary layer thickness ${\delta }_{*}^T$ was defined as the distance between the top wall and the maximum value of $\langle \left|U\right|\rangle$ in the domain $(0\le x\le 0.5$, $0.5\le y\le 1.0)$. Similarly, ${\delta }_{*}^B$ is the distance from the bottom wall to the maximum value of $\langle \left|U\right|\rangle$ in the domain $(0.5\le x\le 1.0$, $0\le y\le 0.5)$. Consequently, two friction Reynolds numbers $R{e}_{\tau }$ were defined as follows:

$$R{e}_{\tau }^T=u_{\tau }^T{\delta }_{*}^T\sqrt{Ra/Pr},R{e}_{\tau }^B=u_{\tau }^B{\delta }_{*}\sqrt{Ra/Pr}.$$

(19)

The boundary layer thickness ${\delta }_{*}$ remained almost constant for $Ra\le 5\times {10}^7$, but it began to diverge from $Ra=1\times {10}^8$, as shown in Figure 9c. A sharp increase appeared at $Ra={10}^9$, after which ${\delta }_{*}$ was consistently larger until $Ra={10}^{10}$. Although the variations in ${\delta }_{*}^T$ and ${\delta }_{*}^B$ with respect to $Ra$ were similar, they showed some differences for $Ra\ge {10}^{10}$. The values of $R{e}_{\tau }$ were almost the same for $Ra={10}^{11}$, as shown in Figure 9d, but, for the highest five cases, the discrepancy was clear, meaning that the double-reflection symmetry was broken.

Based on these results for the TL and BR regions, we can conclude that the flow in RBC remained symmetrical after an x–y reflection for low $Ra$ values in the current cases, i.e., $Pr=0.7$ and ${10}^7\le Ra\le 5\times {10}^8$. However, the symmetry was broken above $Ra={10}^9$. There was a slight deviation between the wall units of the TL and BR regions, indicating that the flow was different after double reflection. This transition at $Ra={10}^9$ was correlated with the flow-pattern transition discussed above [18].

5. Conclusions

In this study, we investigated the symmetry after an x–y double reflection using the DNS data of a 2D RBC in a confined square cell for $Ra$ values ranging from ${10}^7$ to ${10}^{13}$ at $Pr=0.7$. As illustrated by the instantaneous and mean flow fields, both the top–bottom and left–right symmetries were broken within the current range, as was expected. We examined the mean temperature and velocity fields, the horizontal profiles at two specific heights, and the boundary layer properties in the TL and BR regions, including the friction velocity, friction temperature, boundary layer thickness, and friction Reynolds number. To assess the validity of the symmetry, the relative differences between the original and transformed quantities were calculated. It was revealed that the double-reflection symmetry held for $Ra<{10}^9$ but was broken at higher $Ra$ values within the range considered. In contrast to a previous study [45], which found double-reflection symmetry for $Ra$ values up to ${10}^{10}$ at $Pr=6.5$, the present results indicate that the symmetry was broken at $Ra={10}^9$ when the $Pr$ value was decreased to 0.7. Therefore, the effect of $Pr$ was crucial for determining the symmetry properties in a RBC. The observed asymmetry may be related to the flow-pattern transitions previously reported [18], leading to stronger fluctuations in the flow fields. The current findings demonstrate the parameter range for the double-reflection symmetry. The application in the latter in the boundary settings could potentially accelerate the simulations of RBC cases and reduce computational costs, but this will have to be explored in the future.