3.1. Global characteristics and spatial distribution of the dissipation fields

There is extensive visual material available on VC, highlighting key features such as the self-organisation of thermal plumes into a large-scale circulation (LSC), the stable thermal stratification in the bulk and the increase in small-scale turbulent fluctuations as $\textit {Pr}$ decreases. Figure 3 provides additional insight into the spatial distribution of dissipation fields by presenting contour plots of $\log _{10} \varepsilon _u$ and $\log _{10} \varepsilon _T$ for ${\textit {Pr}} = 0.032$ , including both instantaneous and time-averaged fields. The instantaneous kinetic dissipation fields reveal the intense nature of liquid metal convection. As $Ra$ increases, extreme dissipation transitions from a ‘sheet-like’ to a finer ‘vortex-like’ structure, with the bulk becoming more chaotic due to disturbances from boundary layer separation flows. The extreme region of the mean field evolves from diagonally aligned, approximately elliptical shapes to more rectangular structures that better conform to the cavity geometry, indirectly indicating the specific morphology of the LSC. Importantly, the distribution of $\varepsilon _u$ remains predominantly concentrated within the VBL (both vertical and horizontal). The instantaneous thermal dissipation fields, influenced by the large thermal diffusivity, exhibit coarser plume structures that are advected downstream by the LSC. As $Ra$ increases, more small-scale structures extend randomly from the hot and cold plates into the bulk, though the overall spatial scale is significantly larger than the ‘filamentary’ structures observed in RBC at intermediate $\textit {Pr}$ (Xu et al. Reference Xu, Wang, Wan and Sun2021; He, Bao & Chen Reference He, Bao and Chen2022). Similarly, in the time-averaged fields, the contribution of plumes is filtered out, with $\varepsilon _T$ concentrating near the TBL where the temperature gradients are the largest, consistent with the adiabatic condition imposed on the horizontal boundaries.

Figure 3.

Time-averaged and instantaneous kinetic dissipation ( $\log _{10} \varepsilon _u$ ) and thermal dissipation fields ( $\log _{10} \varepsilon _T$ ) for $Ra=10^7$ and $5\times 10^8$ from 2-D simulations.

Figure 4.

(a,b) P.d.f.s of $\log _{10} \varepsilon _T$ and $\log _{10} \varepsilon _u$ calculated over the whole 2-D cell. Here, $\mu$ and $\sigma$ are respectively the mean value and standard deviation. The dashed lines mark the log-normal distribution for comparison. (c,d) P.d.f.s of $\varepsilon _T$ and $\varepsilon _u$ , obtained over the whole cell and normalised with respect to their r.m.s. values. The curves are the best fits of stretched exponential, as given by (3.1).

Figure 4 presents the probability density function (p.d.f.) of the dissipation fields, sampled from at least 50 snapshots after the flow has reached a statistically steady state. It is important to note that the current p.d.f. statistics are based on simulations using a fine, uniform grid, with specific grid details provided in table 1. The original non-uniform grid, with its exponential distribution near the walls, caused data clustering in phase space, actually increasing the weighting of the boundary layer and introducing non-physical information (Xu, Zhang & Xia Reference Xu, Zhang and Xia2024). To avoid this, it is necessary to ensure uniform spatial resolution across the entire field.

Figures 4(a) and 4(b) illustrate the deviation of the log-normalised p.d.f.s from a log-normal distribution, visually depicting the transition of VC from a flow state concentrated at specific patterns at the lowest $Ra$ to a progressively more intermittent and random turbulent state. As $Ra$ increases, the distribution of low dissipation (left tail) gradually converges, and the central region becomes more aligned with the log-normal curve, similar to findings in existing RBC studies (Zhang, Zhou & Sun Reference Zhang, Zhou and Sun2017; Venugopal, De & Mishra Reference Venugopal, De and Mishra2022). A pronounced shoulder appears in the high-dissipation region (right tail), which is believed to be associated with high-frequency extreme events (plume emissions) in the boundary layer, indicating the developed intermittency of localised dissipation. This is more clearly demonstrated in figure 4(b), where the two inflection points on the right can be used to distinguish between the bulk flow, plumes/mixing layers and the conductive sublayers (Kaczorowski & Wagner Reference Kaczorowski and Wagner2009). In other words, the strong dissipation in the right tail originates from the steep mean gradients near the conducting walls, while the weak dissipation in the left tail is associated with local fluctuations exhibited by the bulk flow.

Figures 4(c) and 4(d) show the p.d.f.s normalised by the root-mean-square values. At the lowest $Ra$ , the tail follows an exponential distribution that appears closer to a straight line on a log-linear scale, without any stretching. As $Ra$ increases, the tail gradually broadens and extends towards higher amplitudes, corresponding to the small-scale intermittency caused by more active plume dynamics. Emran & Schumacher (Reference Emran and Schumacher2008) proposed fitting the p.d.f. from the most probable amplitude to the tail using a stretched exponential function:

(3.1) \begin{equation} p(X)=\frac {C}{\sqrt {X}}\exp(-mX^\alpha ). \end{equation}

For passive scalar dissipation in the two-dimensional limit of large Peclet and Prandtl numbers, the stretched exponent $\alpha$ converges to the derived threshold value of 1/3. This functional form has since been widely applied in describing p.d.f.s in various RBC studies. Here, we extend this approach to VC, which shares similar physical characteristics, and obtain a reasonably good fit, with the related parameters listed in the figures. This analysis suggests that, across different thermal turbulence systems (RBC or VC), the instantaneous dissipation fields dominated by plume activity exhibit certain universal properties.

With a basic understanding of the dissipation distribution across the entire field, we can further quantify the contribution from each subregion of the convective cell. The subregions for kinetic dissipation include the vertical boundary layer near the hot and cold walls, the horizontal boundary layer near the top and bottom walls, and the bulk. For thermal dissipation, the subregions consist of the thermal boundary layer and the central bulk. To reasonably distinguish, it is necessary to accurately calculate the boundary layer thickness, which varies depending on the specific thermal turbulence system and boundary conditions. In 2-D RBC, Zhang et al. (Reference Zhang, Zhou and Sun2017) defines the VBL as the location of the maximum r.m.s. velocity to avoid non-physical thickness caused by corner vortices. Zhang, Dong & Xia (Reference Zhang, Dong and Xia2022), in inclined convection, defines the VBL thickness as the $x$ -direction position of the maximum vertical velocity. Lim et al. (Reference Lim, Chong, Ding and Xia2019) defines the stress boundary layer thickness as the location of the maximum sum of the velocity normal gradient squares. In VC, since the large-scale turbulent ‘wind’ is easily identifiable, the VBL thickness generally follows the laminar Prandtl–Blasius (P-B) scaling (Ng et al. Reference Ng, Ooi, Lohse and Chung2015). It can be defined as the point where the tangent to the mean velocity profile at the wall equals the turbulent wind velocity ( $\delta _u = U / ({\rm d}\bar {u}/{\rm d}x|_w)$ ), a method also applied in most RBC configurations.

Relatively speaking, the definition of TBL thickness is more consistent, typically employing the slope method, which defines the thickness as the intercept of the tangent to the mean temperature profile at the wall with the horizontal axis ( $\delta _T = -(\Delta T/2)/({\rm d}\bar {T}/{\rm d}x|_w)$ ). Additionally, Howland et al. (Reference Howland, Ng and Lohse2022) defines the conductive TBL in high- $Pr$ VC as the wall-normal location where the conductive heat flux equals the turbulent heat flux ( $-\kappa {\rm d}\bar {T}/{\rm d}x = \overline {u'T'}$ ), with the heat flux being dominated by mean molecular diffusion. Ng et al. (Reference Ng, Ooi, Lohse and Chung2015) compared three boundary layer definitions: the tangent intercept method, the dissipation cross-over point and the displacement thickness, demonstrating that the scaling relations were largely consistent,

(3.2) \begin{equation} \delta _u /H \approx 0.43Re^{-1/2}, \quad \delta _T /H \approx 2.1Re^{-1/2}. \end{equation}

Following Ng et al. (Reference Ng, Ooi, Lohse and Chung2015), Zwirner et al. (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022) and Shishkina (Reference Shishkina2016), we adopt the slope method to define the VBL and TBL thickness.

Figure 5 illustrates the percentage of dissipation in each subregion for different $Ra$ in the present 2-D simulations. Specific boundary layer information can be found in table 7 in the Appendix. The results indicate that the relative contributions of the VBL and bulk to $\varepsilon _u$ remain nearly constant, staying around $60\, \%$ and $40\, \%$ , respectively. Notably, $\varepsilon _{u, BL}$ from the horizontal boundary layer is a significant source of kinetic dissipation, comparable in magnitude to that from the generally considered vertical boundary layer. Additionally, the fraction of both alone is less than the bulk dissipation. As $Ra$ increases, the proportion of dissipation in the horizontal boundary layer decreases, though it still exceeds half of that in the vertical boundary layer, even at the highest $Ra$ . Also in the current parameter range, the contribution to $\varepsilon _T$ is also dominated by the TBL. Although the bulk contribution slightly increases with $Ra$ , it remains significantly lower than that of the TBL.

Figure 5.

Percentages of the bulk and BL contributions to the (a) kinetic and (b) thermal dissipation rates for different $Ra$ at $Pr=0.032$ .

Table 3.

A summary of the dominant subregions for kinetic and thermal dissipation that can be clearly identified from recent publications. The checkmark means that the corresponding dissipation is dominated by this subregion. The sources are as follows: (i, ii) Zhang et al. (Reference Zhang, Zhou and Sun2017), (iii) He et al. (Reference He, Bao and Chen2022), (iv) Zhang & Zhou (Reference Zhang and Zhou2024), (v) Xu, Shi & Xi (Reference Xu, Shi and Xi2019), (vi) Guo et al. (Reference Guo, Wang, Wu, Chong and Zhou2022), (vii) Ng et al. (Reference Ng, Ooi, Lohse and Chung2015).

For ease of reference, table 3 summarises the contributions from the bulk and boundary layer to the kinetic (thermal) dissipation rate that can be obtained clearly from several recent studies on thermal convection. It can be observed that for 2-D thermal convection models with unit aspect ratio, either RBC or VC, both kinetic and thermal dissipation are dominated by the boundary layer across a wide range of parameters. Nonetheless, the local flow organisation in thermal convection systems is strongly correlated with the aspect ratio of the container, and velocity and temperature profiles may exhibit considerably complex behaviour, which has a significant impact on the boundary layer thickness and dissipation statistics in different subregions. For example, Ng et al.’s (Reference Ng, Ooi, Lohse and Chung2015) 3-D VC simulations for a computational domain of $8H\times 4H \times H$ found $\varepsilon _u$ to be dominated by the bulk and $\varepsilon _T$ to remain dominated by the TBL. It might be due to the use of periodic boundary conditions in both directions except for the wall normal phase that the statistics on $\varepsilon _u$ includes only the vertical VBL and does not contain the horizontal VBL. Our study indicates that both contribute equally to kinetic dissipation. Additionally, Grossmann & Lohse (Reference Grossmann and Lohse2002) proposed that the theoretical $Ra$ scaling of the VBL thickness in RBC aligns more closely with that of the sidewalls rather than the top and bottom plates, suggesting that the main contribution to $\varepsilon _u$ appears to come from the sidewalls. For instance, in cylindrical RBC with unit aspect ratio, the contribution from the sidewalls is twice that of the top and bottom combined. Of course, for simplicity, the scaling theory does not distinguish between the positions of the plates and walls in the flow, but instead assumes a characteristic effective thickness in an averaged sense.

3.2. Re dependence of bulk and BL dissipations

Based on the aforementioned statistics of dissipation rates in the bulk and boundary layer, we can examine whether the present low- $\textit {Pr}$ VC satisfies the scaling laws proposed by the GL theory. The key relations are briefly outlined as follows:

(3.3) \begin{align} \varepsilon _{u,bulk}\sim \frac {U^3}{H}=\frac {\nu ^3}{H^4}Re^3,\quad\qquad \end{align}

(3.4) \begin{align} \varepsilon _{u,BL}\sim \nu \frac {U^2}{\delta _u^2}\frac {\delta _u}{H}\sim \frac {\nu ^3}{H^4}Re^{5/2},\qquad \end{align}

(3.5) \begin{align} \varepsilon _{T,bulk}\sim \frac {U\Delta T^2}{H}=\kappa \frac {\Delta T^2}{H^2}RePr,\quad \end{align}

(3.6) \begin{align} \varepsilon _{T,BL}\sim \kappa \frac {\Delta ^2}{\delta _T^2}\frac {\delta _T}{H}\sim \kappa \frac {\Delta T^2}{H^2}Re^{1/2}Pr^{1/2}. \end{align}

Figure 6.

Dissipation scaling laws in the bulk and BL for (a) $\varepsilon _u$ and (b) $\varepsilon _T$ . A distinction is made between total dissipation ( $\varepsilon _{u(T), bulk}$ ) and fluctuations ( $\varepsilon _{u'(T'), bulk}$ ) in the bulk.

Existing research has provided robust validation of the GL theory for RBC and VC in medium- $\textit {Pr}$ fluids, but a common issue persists: while the kinetic and thermal dissipation scaling within the boundary layer can be perfectly reproduced, the contribution from the bulk significantly deviates from the predicted $Re$ dependence, with bulk and boundary layer having almost parallel trends. Grossmann & Lohse (Reference Grossmann and Lohse2004) attributed this to the contamination of bulk scaling by plumes detaching from the TBL. Ng et al. (Reference Ng, Ooi, Lohse and Chung2015), however, considered the possibility that turbulent scaling is obscured by the strong mean flow component. Zhang et al. (Reference Zhang, Zhou and Sun2017) proposed that the Kolmogorov energy cascade does not hold in 2-D RBC, as all calculated third-order velocity structure functions are positive.

Figure 6 depicts the trends from the boundary layer and bulk scalings for 2-D simulations. Analogous to the partitioning occupancies demonstrated in figure 5, the boundary layer contributions are in better agreement with the GL predictions, while $\varepsilon _{u, bulk}$ maintains a roughly constant proportionality with $\varepsilon _{u, BL}$ . Comparatively, the scaling exponent of $\varepsilon _{T, bulk}$ is greater than 1/2, indicating that it is less affected by plume contamination, as the higher thermal diffusivity leads to shorter lifetime of local plumes, with fewer chances of penetrating into the bulk. Nevertheless, it is still far from the predicted $Re$ dependence. Considering the inherent non-zero mean flow velocity in VC and following Ng et al.’s (Reference Ng, Ooi, Lohse and Chung2015) approach, we subtracted the mean field from the bulk dissipation to isolate the fluctuations for statistical analysis:

(3.7) \begin{equation} \langle \varepsilon _{u'}\rangle =\langle \varepsilon _{u}\rangle -\langle \varepsilon _{\bar {u}}\rangle , \quad \langle \varepsilon _{T'}\rangle =\langle \varepsilon _{T}\rangle -\langle \varepsilon _{\bar {T}}\rangle , \end{equation}

where

(3.8) \begin{equation} \langle \varepsilon _{\bar {u}}\rangle =\frac {2}{H}\int _{\delta _u}^{H/2}\nu \left(\frac {{\rm d}\bar {u}}{{\rm d}x}\right)^2{\rm d}x, \quad \langle \varepsilon _{\bar {T}}\rangle =\frac {2}{H}\int _{\delta _T}^{H/2}\kappa \left(\frac {{\rm{d}}{\bar {\rm{T}}}}{{\rm d}x}\right)^2\rm{d}x. \end{equation}

The results successfully reproduce the corresponding theoretical scalings, further confirming that in VC, the small-scale turbulent transport driven by the large-scale mean field is consistent with the Kolmogorov framework. Notably, as $Ra$ increases, the proportion of dissipation attributed to the mean field (the difference between the solid and hollow triangles in figure 6) gradually decreases, indicating that turbulent fluctuations progressively dominate. This observation mirrors Wang et al.’s (Reference Wang, Liu, Shishkina and Lohse2021) results for high- $\textit {Pr}$ VC, although the transition $Ra$ in their case is significantly larger.

Figure 7.

Time and surface averaged horizontal profiles of (a) the vertical velocity component and (b) the temperature near the heated wall. The dashed line represents the tangent at the centre of the profile.

It should be emphasised that the current treatment is not identical to that of Ng et al. (Reference Ng, Ooi, Lohse and Chung2015), as no additional correction has been made to the characteristic velocity using the midline-averaged gradient. The key distinction lies in the velocity profiles presented in the reference (figure 7 from Ng et al. Reference Ng, Ooi, Lohse and Chung2015), which show a well-developed velocity distribution extending from the near-wall region to the bulk centre up to $Ra = 2 \times 10^7$ . The mean velocity obtained from the central gradient is a direct value, closely approximating the maximum velocity of the turbulent ‘wind’. Figure 7 illustrates the surface-averaged velocity profiles acquired in this study. As $Ra$ increases, the peak velocity rises sharply and the overall jet becomes more confined near the wall. Meanwhile, the velocity outside the jet, closer to the bulk, drops drastically, with the vertical velocity and temperature in a substantial part of the bulk remaining close to zero. At this stage, the dashed line representing the central tangent in the figure fails to reflect the average velocity or temperature gradient in the bulk. The fundamental cause of this profile discrepancy might be explained by Zwirner et al.’s (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022) figure 8. It is evident that the closer the aspect ratio is to 1, the more fully mixed the bulk becomes, which is characterised by an isothermal profile close to $(T_h +T_c)/2$ and nearly zero vertical velocity. As $\Gamma$ increases, the mean gradient at the bulk centre increases. This difference can be confirmed to be weakly related to $\textit {Pr}$ , since the current low- $\textit {Pr}$ profiles exhibit similar behaviour to those in Wang et al.’s (Reference Wang, Liu, Shishkina and Lohse2021) high- $\textit {Pr}$ cases. Despite this, it cannot be determined to be independent of the boundary conditions, as both Ng et al. (Reference Ng, Ooi, Lohse and Chung2015) and Howland et al. (Reference Howland, Ng and Lohse2022) applied a periodic boundary in the streamwise and spanwise directions with an aspect ratio of $8\times 4 \times 1$ . In contrast, both the current study and that of Wang et al. (Reference Wang, Liu, Shishkina and Lohse2021) consider a 2-D enclosed model with horizontal no-slip boundaries $( \Gamma = 1)$ . The presence of horizontal plates has a crucial influence on VC, primarily through the generation of corner vortex structures. The temperature profiles intuitively display the so-called ‘undershoot’ phenomenon, where the temperature drops below the bulk temperature after descending from $T_h$ before oscillating back to the bulk value. This results in the central gradient, as indicated by the dashed line in the figure, being opposite to the mean gradient. The reason for this can be understood as the rate of heat transfer from the wall to the outer boundary layer is not sufficient to keep the temperature synchronised with that of the exterior as the fluid rises vertically (Ravi, Henkes & Hoogendoorn Reference Ravi, Henkes and Hoogendoorn1994). As a consequence, cooler fluid relative to the core appears, descending as it passes through the corners and separating from the horizontal wall, forming additional circulations.

In conclusion, by filtering out the obscuring effect of the strong mean field, we are able to reproduce the classical theoretical scalings. Nevertheless, the influence of the vessel aspect ratio, including boundary conditions, on the local flow structures, profile shapes and related statistics requires additional consideration.

3.3. Modelling of global dissipation and extension to classical VC

The starting point of the GL theory lies in two exact relations for kinetic dissipation and thermal dissipation, (1.2) and (1.3), which can be derived by multiplying the general form of the momentum equation by $ u$ , and the temperature transport equation by $ T$ , followed by applying Gaussian integration and spatiotemporal averaging under steady-state conditions (Siggia Reference Siggia1994). In fact, the more general relation obtained from the kinetic energy equation is

(3.9) \begin{equation} \varepsilon _u=\frac {\nu ^3}{H^4}\frac {\langle u_y T\rangle }{\kappa \Delta /H}RaPr^{-2}, \end{equation}

where $ u_y$ is the velocity component in the direction opposite to gravity. This indicates that the buoyant potential energy represented by the vertical heat flux $\langle u_yT\rangle$ is converted to kinetic energy and then dissipated by viscous force. In RBC, the heat flux aligns with the buoyancy direction ( $Nu=\langle u_y T-\kappa \partial _y T\rangle /(\kappa \Delta /H)$ ) and (3.9) can be rewritten as (1.2).

In contrast, in VC, the mean horizontal heat flux is orthogonal to the buoyancy direction ( $Nu=\langle u_x T-\kappa \partial _x T\rangle /(\kappa \Delta /H)$ ). Therefore, the exact relationship that links the global dissipation rate to dimensionless parameters no longer holds. Meanwhile, the thermal dissipation rate depends solely on the $Nu$ and retains the same form as in the classical case. Thus, designing an accurate global kinetic dissipation model is a key point in extending the GL theory to the context of VC. Currently, there have been some forward-looking and highly valuable studies seeking alternative definitions for this unclosed dissipation relation. For instance, Ng et al. (Reference Ng, Ooi, Lohse and Chung2015) discovered that the horizontal heat flux component satisfies a corresponding scaling and also exhibits a 1/2 power-law in homogeneous VC (Ng et al. Reference Ng, Ooi, Lohse and Chung2018). Zhang et al.’s (Reference Zhang, Ding and Xia2021) successful modelling of the effective horizontal buoyancy in inclined convection is even more encouraging. Despite this, a complete extension of the GL theory to VC and a comparative analysis with existing scaling theories are still lacking at present.

Equation (3.9) informs us that dissipation is proportional to the system’s volume-averaged vertical heat flux, $ \langle u_y T \rangle$ . After non-dimensionalising with characteristic velocity $ u_0 = (\beta g\Delta TH)^{1/2}$ and temperature $ T_0 = \Delta T$ , the kinetic dissipation rate can be written as

(3.10) \begin{equation} \varepsilon _u = \frac {\nu ^3}{H^4}\langle u_y^*T^*\rangle Ra^{3/2}Pr^{-3/2}, \end{equation}

where $u_y^*$ and $T^*$ represent dimensionless vertical velocity and temperature. Taking into account the non-zero mean streamwise velocity in VC, combined with a careful examination of the velocity and temperature profiles, it can be demonstrated that vertical heat flux is primarily concentrated within the TBL (Zwirner et al. Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022), as shown in figure 8(a). As previously mentioned, in the bulk, due to strong flow mixing and thermal stratification, both mean vertical velocity and mean temperature are nearly zero, contributing very little to the overall vertical heat transfer. This phenomenon of concentrated heat flux is remarkably similar to the ‘heat transport channels’ observed in inclined convection, where the temperature and vertical velocity are highly correlated, leading to a significant enhancement in local heat flux. The inset of figure 8(a) shows that the temperature profile converges to 0 much faster despite the larger jet width, i.e. the thickness of the TBL is representative of the actual width of the heat flux (in grey). The following question arises. Can we approximate the total vertical heat flux of the entire field by statistically analysing this portion alone?

Figure 8.

(a) Mean vertical heat flux profiles for different $Ra$ , where the shaded area of the inset represents the TBL used to approximate statistical kinetic dissipation. (b) Comparison of exact and approximate values of global kinetic dissipation from 3-D DNS, where the inset shows the relative errors. Specific data can be found in table 8 in the Appendix.

The modelling of vertical heat flux within the TBL in this context involves several assumptions or approximations.

1. (i) For low- $\textit {Pr}$ fluids, the VBL is nested within the TBL and the characteristic velocity has developed into the so-called ‘turbulent wind velocity’.

2. (ii) The temperature distribution within the TBL is approximately linear, $\bar {T}=T_h +{\rm d}\bar {T}/{\rm d}x|_wx$ .

3. (iii) The thickness of the TBL is directly defined by the temperature gradient at the plate, $\delta _T=-(\Delta T/2)/{\rm d}\bar {T}/{\rm d}x|_w$ .

Conditions (ii) and (iii) are evidently valid, as they can be sufficiently verified or analytically derived from the current profiles or published work (Grossmann & Lohse Reference Grossmann and Lohse2004; Ng et al. Reference Ng, Ooi, Lohse and Chung2015). For condition (i), although $ u_y$ within the TBL varies, considering the very thin VBL and the fact that the vertical velocity remains close to the ‘turbulent wind velocity’ within the statistical region, it is still feasible to approximate with a constant value. This will naturally lead to an overestimation, but subsequent quantitative comparisons demonstrate that the resulting error is acceptable. An independent approximation of the product of the vertical heat flux is then obtained:

(3.11) \begin{equation} \langle u_y^*T^*\rangle \sim u_{wind}^*\left \langle \, T^*\right \rangle \sim \frac {RePr^{1/2}}{Ra^{1/2}}\int _{0}^{\delta _T}\left(1-\frac {x}{\delta _T}\right){\rm{d}}x = \frac {RePr^{1/2}}{4NuRa^{1/2}}, \end{equation}

which results in a fairly reasonable estimate of the global kinetic dissipation in VC (substitute (3.11) into (3.10)),

(3.12) \begin{equation} \varepsilon _u=\frac {\nu ^3ReRa}{4H^4NuPr}. \end{equation}

To verify (3.12), the exact and approximate values of kinetic dissipation are examined in 3-D VC simulations ( $\Gamma =1$ ) and the expression (1.2) applied to the RBC is additionally calculated for comparison. As shown in figure 8(b), the exact values are very close to the estimated ones, with the relative error, displayed in the inset, ranging between 6 $\, \%$ and 28 $\, \%$ . In contrast, the classical formula (1.2) significantly overestimates the actual dissipation, with errors ranging from 65 $\, \%$ to 215 $\, \%$ . In fact, the ratios, $\varepsilon _u/(\nu ^3H^{-4}(Nu-1)Ra)$ , measured by Shishkina (Reference Shishkina2016) in a cylindrical container with equal height and diameter, range from 0.2 to 0.6 ( $Pr = 0.1, Ra = 10^5{-} 10^{10}$ ), which are weakly correlated with $Ra$ . Similarly, Wang et al. (Reference Wang, Liu, Shishkina and Lohse2021) obtained a ratio of approximately 0.4 in 2-D VC, corresponding to a relative error of approximately 150 $\, \%$ . All these findings suggest that the classical formula is no longer applicable in VC due to the orthogonality between buoyancy and heat flux.

It is worth cautioning of the error at the highest $Ra$ in the current approximation, which cannot be overlooked. This calls for careful consideration of the model’s applicable range. The starting point of the modelling is an approximate statistic of the vertical heat flux, and the corresponding error arises from two main sources. On the one hand, it stems from our choice of a constant Reynolds number: as indicated by the inspection of velocity profiles in figure 7, the jet width narrows sharply with increasing $Ra$ , leading to a very steep velocity distribution within the TBL. In this case, using the peak $Re$ undoubtedly overestimates the actual vertical heat flux. On the other hand, the concentration of heat flux in the TBL requires the velocity profile to exhibit a full peak near the wall, which should be representative of the overall flow. The 2-D VC velocity profiles with $Pr$ = 10 shown by Wang et al.’s (Reference Wang, Liu, Shishkina and Lohse2021) figures 3 and 4 indicate that after $Ra$ reaches $2\times 10^{11}$ , the horizontal flow, originally concentrated near the top and bottom walls, evolves into a strong horizontal velocity in the bulk as zonal flow. Meanwhile, the peak vertical velocity diminishes due to more active plume dynamics, attenuating the overall vertical motion. At this point, any statistical treatment of the vertical heat flux becomes inaccurate. Moreover, for low- $Pr$ fluids ( $Pr$ = 0.03), the transition $Ra$ will undoubtedly occur even earlier.

Figure 9.

(a,b) $Ra$ dependences of the compensated forms of $Nu/Ra^{1/4}$ and $Re/Ra^{1/2}$ from Shishkina (Reference Shishkina2016) and Zwirner et al. (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022). (c,d) Reproductions of Shishkina’s (Reference Shishkina2016) figures 3(b) and 3(d): the $Pr$ dependences of $Nu$ and $RePr$ . The solid lines represent the scalings from (3.14).

Nevertheless, the above approximate model provides a relatively accurate alternative to the global kinetic dissipation statistics of VC under the current parameter range. Below, we apply (3.12) to the extended GL theory and compare the results with existing studies and well-established scaling theories. As shown in figure 5, the dissipation fields in the present simulations are dominated by the boundary layer, allowing for the simplification of (1.4) and (1.5) to obtain an effective pure power-law scaling. The joint equations

(3.13) \begin{equation} \frac {ReRa}{NuPr} \sim Re^{5/2} \quad \text {and} \quad Nu \sim Re^{1/2}Pr^{1/2} \end{equation}

are solved as follows:

(3.14) \begin{equation} Nu \sim Ra^{1/4}Pr^{1/8}, \quad \quad Re\sim Ra^{1/2}Pr^{-3/4}. \end{equation}

Surprisingly, the scaling (3.14) is identical to the predicted scaling for the $I_l$ regime in the unified RBC theory (Grossmann & Lohse Reference Grossmann and Lohse2000), despite the expressions for $\varepsilon _u$ being quite different. In comparison to the theoretical scaling obtained by Shishkina (Reference Shishkina2016) based on the Prandtl method,

(3.15) \begin{equation} Nu \sim Ra^{1/4}Pr^{1/4}, \quad \quad Re\sim Ra^{1/2}Pr^{-1/2}, \end{equation}

the classical $Ra$ dependence can be reproduced, i.e. $Nu \sim Ra^{1/4}$ and $Re \sim Ra^{1/2}$ . As shown by the compensated quantities from Shishkina (Reference Shishkina2016) and Zwirner et al. (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022) in figures 9(a) and 9(b), the results agree well with the current derivation and previously established theory. Upon narrowing the vertical axis range, slight fluctuations in the trend become noticeable. In particular, the compensated $NuRa^{-1/4}$ increases with $Ra$ , which could be due to the geometric constraints of the convective cells affecting heat and mass transfer, or it might suggest a transition towards the $Ra^{1/3}$ scaling ( $Ra \gt 10^7$ , Ching Reference Ching2023). At present, there is insufficient evidence to provide a definitive answer.

It is noteworthy that the current $Pr$ exponent is slightly different from the laminar scaling proposed by Shishkina (Reference Shishkina2016). Here, we reproduce figures 3(b) and 3(d) of that article, i.e. the $Pr$ dependence within an equal height and diameter cylindrical cell for fixed $Ra = 10^6$ and $10^7$ in figure 9(c) and 9(d). Interestingly, the scalings $Nu \sim Pr^{1/4}$ and $Re \sim Pr^{-1/2}$ fit well in the small $Pr$ limit ( $Pr \lt 0.03$ ), with a slight overestimation between 0.03 and 1 (especially for $RePr$ ). In contrast, the current scalings $Nu \sim Pr^{1/8}$ and $Re \sim Pr^{-3/4}$ , are supported by the data over almost the full range where $Pr \lt 1$ . Let us revisit the origins of these two $Pr$ scalings: the former is based on an extended Prandtl approach, providing a similarity solution to the heat and momentum transport problem on a vertical heated plate. For the detailed derivation, please refer to the original work (Shishkina Reference Shishkina2016); we briefly outline the method here to facilitate the following analysis.

(3.16) \begin{align} u_x\partial _xu_x+u_z\partial _zu_x = \nu \partial _z^2u_x+\beta g(T-T_0), \end{align}

(3.17) \begin{align} u_x\partial _xT+u_z\partial _zT = \kappa \partial _z^2T,\qquad\quad \end{align}

(3.18) \begin{align} \partial _xu_x+\partial _zu_z=0.\qquad\qquad \end{align}

Note that the $x$ direction is along the plate and the $z$ direction is horizontal away from the plate. Construct the following similar variable $\xi$ , stream function $\psi$ and temperature for the Prandtl boundary layer equations (3.16)–(3.18):

(3.19) \begin{align} \xi = Pr^aRa^b(z/H)(x/H)^c, \end{align}

(3.20) \begin{align} \psi = Pr^dRa^e\nu (x/H)^f\phi (\xi ), \end{align}

(3.21) \begin{align} T = T_0+(\Delta T/2)\theta (\xi ). \quad\end{align}

The Nusselt number is obtained from the average temperature gradient at the plate, while the Reynolds number is determined based on the maximum mean velocity ( $u_x=\partial _z\psi$ ),

(3.22) \begin{align} Nu = \frac {-H^{-1}\int _{0}^{H}\kappa \partial _zT|_{z=0}{\rm{d}}x}{\kappa \Delta /H}=\frac {-\theta _\xi (0)}{2(c+1)}Pr^aRa^b, \qquad \end{align}

(3.23) \begin{align} Re = \underset {z}{\max }\ H^{-1}\int _{0}^{H}u_x{\rm{d}}x = Pr^{a+b}Ra^{b+e}(c+f+1)^{-1}\phi _\xi . \end{align}

Substitute similarity variables into the heat transport equation (3.17),

(3.24) \begin{equation} \theta _{\xi \xi }+fPr^{d-a+1}Ra^{e-b}(x/H)^{f-c-1}\phi \theta _\xi =0, \end{equation}

and eliminate half of the parameters by requiring independence from $Ra$ and $Pr$ ,

(3.25) \begin{equation} d=a-1,\quad e=b, \quad f=c+1. \end{equation}

Similarly, substituting (3.19)–(3.21) into the momentum equation (3.16) yields

(3.26) \begin{equation} Pr^{4a}\phi _{\xi \xi \xi }+Pr^{4a-1}\left[(c+1)\phi \phi _{\xi \xi }-(2c+1)(\phi _\xi )^2\right]-Ra^{1-4b}(x/H)^{-4c-1}\theta /2=0. \end{equation}

Considering the independence of the momentum equation from $Ra$ and the existence of a similarity solution with respect to $\xi$ ,

(3.27) \begin{equation} b=1/4, \quad c=-1/4. \end{equation}

Equation (3.26) can be simplified as

(3.28) \begin{equation} 4Pr^{4a}\phi _{\xi \xi \xi }+Pr^{4a-1}\left[3\phi \phi _{\xi \xi }-2(\phi _\xi )^2\right]-2\theta =0. \end{equation}

For $ Pr \ll 1$ , the first term on the left-hand side of the above equation can be neglected compared with the second term, yielding $ a = 1/4$ , with the corresponding scaling as

(3.29) \begin{equation} Nu \sim Ra^{1/4}Pr^{1/4}, \quad Re\sim Ra^{1/2}Pr^{-1/2}. \end{equation}

For $ Pr \gg 1$ , the first term in (3.28) dominates and $a=0$ . In this case,

(3.30) \begin{equation} Nu \sim Ra^{1/4}Pr^{0}, \quad Re\sim Ra^{1/2}Pr^{-1}. \end{equation}

The above scalings are well validated in figure 9, and the same similarity solution approach has been successfully applied in deriving the scaling laws for large- $Pr$ RBC, i.e. $Nu \sim Ra^{1/3}Pr^{0},\ Re\sim Ra^{2/3}Pr^{-1}$ (Shishkina et al. Reference Shishkina, Emran, Grossmann and Lohse2017). Despite this, as mentioned in the article, ‘Since the $Pr$ -dependence holds only for very large or small $Pr$ , it formally breaks down for intermediate values of $Pr$ , which strictly speaking means the absence of the similarity solution’. This deviation is particularly noticeable in the range of $Pr = 0.03{-}1$ in figure 9, even though the scaling law remains entirely valid in a broader sense. A careful examination of the data trends in figure 9(c) suggests that the $Pr$ exponent $ a$ should vary continuously (slow transition from 1/4 to 0 at $Pr\lt 1$ , the 1/4 scaling may be more accurately captured when $Pr\lt 10^{-2}$ ).

For the case of $ Pr \lt 1$ (but not significantly smaller), it may not be appropriate to simply neglect the subdominant terms in (3.28), and the exponent $ a$ could lie between 0 and 1/4. The $Pr$ exponent for $Re$ , $ a + d = 2a - 1$ , may be less than -1/2. This aligns with the accepted view that, while pure power-law scaling is sufficiently valid, it may not be the optimal description of transport relationships. Returning to the model presented in this paper, the core focus is on the exact dissipation relations in a thermal convection system and the spatial distribution of vertical heat flux, which are connected to the BL-dominated nature of $ \varepsilon _u$ and $ \varepsilon _T$ . By solving a set of simple coupled equations, the scaling relations can be derived. The $Pr$ exponent of 1/8 for $Nu$ and $-3/4$ for $Re$ correspond well with the results for $ Pr$ between $ 10^{-2}\ \rm{and}\ 10^0$ . This exponent of 1/8 is based on the physical picture of a BL-dominated system and the profile morphology, representing a transition or average between 0 and 1/4 in the intermediate regime.

In fact, our discussion extends beyond this, and further exploration may yield even more interesting results. Although the subject of this paper is on vertical convection in low- $Pr$ fluids, we can reasonably make some predictions about the scaling relationships for $ Pr \gt 1$ . As discussed by Grossmann & Lohse (Reference Grossmann and Lohse2001) and Zhang et al. (Reference Zhang, Ding and Xia2021) regarding the ‘lower’ and ‘upper’ subregions, for high- $Pr$ fluids, the TBL is nested within the VBL. In terms of the profiles, this implies that at the point where the temperature curve converges from $ T_h$ to $ (T_h+T_c)/2$ , the velocity has not yet developed to the ‘turbulent wind’ velocity, but is instead represented as $ U{\delta _T}/\delta _u$ . At this stage, the expression for dimensionless vertical heat flux can be approximately written as

(3.31) \begin{equation} \langle u_y^*T^*\rangle \sim \frac {RePr^{1/2}}{Ra^{1/2}}\int _{0}^{\delta _T}\frac {x}{\delta _u}\left(1-\frac {x}{\delta _T}\right){\rm{d}}x = \frac {Re^{3/2}Pr^{1/2}}{24aNu^2Ra^{1/2}}, \end{equation}

where $\delta _u \sim a/\sqrt {Re}$ . The scalings shown in figures 9(c) and 9(d) indicate that for fixed $Ra$ of $ 10^6$ and $ 10^7$ , when $ Pr \gt 1$ , $Nu$ is independent of $Pr$ , while $Re\sim Pr^{-1}$ . To ensure rigour, we conduct additional simulations using fine grids for fluids with $ Pr = 4.3$ (water) and $ Pr = 10$ (lubricating oil) under the same parameter conditions to evaluate the corresponding profiles and dissipation field distributions. The results demonstrate that both kinetic and thermal dissipation are also dominated by the VBL and TBL. To avoid distracting from the main focus of this paper, a detailed description of the high- $Pr$ VC simulations is provided in the Appendix, though we strongly encourage readers to refer to this section. Similarly, by solving the following equations:

(3.32) \begin{equation} \frac {Re^{3/2}Ra}{Nu^2 Pr} \sim Re^{5/2} \quad \text {and} \quad Nu \sim Re^{1/2}Pr^{1/3}, \end{equation}

the corresponding scalings are

(3.33) \begin{equation} Nu \sim Ra^{1/4}Pr^{-1/12}, \quad \quad Re\sim Ra^{1/2}Pr^{-5/6}. \end{equation}

Note that when $ Pr \gt 1$ , the form of $ \varepsilon _{T,BL}$ differs in (3.32) (the exponent of $Pr$ ). Equation (3.33) compared with (3.30) has the same $Ra$ exponent but exhibits a slight difference in the $Pr$ exponent. To clarify this further, table 4 lists several scalings, including the current derivation (3.33), the laminar VC theory (3.30) (Shishkina Reference Shishkina2016) and the effective power laws recently obtained by Howland et al. (Reference Howland, Ng and Lohse2022) from high- $Pr$ VC simulations. As shown, the current $Pr$ -dependence is more consistent with the results of Howland et al. (Reference Howland, Ng and Lohse2022), with $Nu$ showing weak dependence on $Pr$ , and the $Pr$ exponent of $Re$ being slightly greater than -1. Howland et al. (Reference Howland, Ng and Lohse2022) attributes the deviation from laminar VC to the strong turbulence from their set-up, but this can also be explained within the GL framework by the dominant dissipation relationship.

Table 4.

Observed and derived scalings laws for various dimensionless response parameters from large- $Pr$ VC. The effective power laws from Howland et al. (Reference Howland, Ng and Lohse2022) for $1\leqslant Pr\leqslant 100$ , theoretical scalings for laminar VC from Shishkina (Reference Shishkina2016) and the present scalings are provided for comparison. Note that the following relations exist: $Re_{\tau }=H\sqrt {d\bar {u}_y/dx|_{w}}/\nu ^{1/2}$ , $C_D=Re_\tau ^2/Re^2$ , $C_T=2Nu/(Re_\tau Pr)$ and the Blasius drag law $C_D\sim Re^{-1/2}$ .

Overall, whether one starts from Prandtl boundary layer theory or from the exact dissipation relations, both approaches converge on the similar conclusion, showing full harmonisation. The $Ra$ dependence of all the response parameters is consistent, while the $Pr$ dependence shows slight divergence in applicability. For $ Pr \lt 1$ , the current approximation, yielding $ a = 1/8$ , provides an averaged characterisation in the intermediate range, which requires the system to be in a VBL–TBL-dominated state, thus imposing certain parameter limitations. If $Pr\ll 1$ and $Ra$ is extremely large, on the one hand, the VBL becomes extremely thin, potentially leading to $\varepsilon _u$ being dominated by the bulk. On the other hand, the high turbulence intensity in the flow field causes the approximate statistic of vertical heat flux to break down. In such cases, the scaling $ a = 1/4$ holds true. Currently, research on thermal convection in extremely low- $Pr$ fluids is quite scarce (Zwirner et al. (Reference Zwirner, Khalilov, Kolesnichenko, Mamykin, Mandrykin, Pavlinov, Shestakov, Teimurazov, Frick and Shishkina2020), $Pr\approx 0.009$ ). Determining the critical regimes where different dissipation mechanisms dominate presents a significant challenge, both in terms of computational resources and statistical analysis. For $ Pr \gt 1$ , the two sets of scaling relations are also unified and $Nu$ exhibits a weak dependence on $Pr$ , as the VBL is constrained by the vessel geometry. As $Pr$ expands, the flow tends to converge (Grossmann & Lohse Reference Grossmann and Lohse2001). Even so, there are also limitations in terms of $Ra$ range, because both the current approximation and boundary layer theory derivations are based on a common core assumption: the buoyancy-driven turbulent ‘wind’ is not strong enough to sustain the turbulent boundary layer in the classical Prandtl and von Kármán sense. As validated by Wang et al. (Reference Wang, Liu, Shishkina and Lohse2021) in 2-D VC simulations with $ Pr = 10$ , this transition $Ra$ number is approximately $ 5 \times 10^{10}$ . Beyond this threshold, the scalings $ Nu \sim Ra^{1/3}$ and $ Re \sim Ra^{4/9}$ emerge, which are consistent with the predictions of the unified RBC theory in the $IV_u$ regime (bulk–bulk dominated). The consistency of the VC scaling with the unified theory suggests a deeper connection between the transport mechanisms in RBC and VC, one that far exceeds previous expectations. This also urges us to treat the conclusions drawn here with great caution.

In the spirit of GL theory, substituting the generalised expressions (3.3)–(3.6) for local dissipation and the approximate statistics of global dissipation into the conservation equations gives the complete functional dependence $Nu(Ra, Pr)$ and $Re(Ra, Pr)$ , rather than being limited to asymptotic power laws

(3.34) \begin{align} \frac {\nu ^3 ReRa}{4H^4NuPr}=c_1\frac {\nu ^3}{H^4}Re^{5/2}+c_2\frac {\nu ^3}{H^4}Re^3, \qquad\qquad \end{align}

(3.35) \begin{align} \kappa \frac {(\Delta T)^2}{H^2}(Nu-1)=c_3\kappa \frac {(\Delta T)^2}{H^2}Pr^{1/2}Re^{1/2}+c_4\kappa \frac {(\Delta T)^2}{H^2}PrRe. \end{align}

Compared with the original set of equations applicable to RBC (Stevens et al. Reference Stevens, Poel, Grossmann and Lohse2013), we have the following.

1. (i) The expression for global kinetic dissipation has been modified.

2. (ii) The transition function that represents the cross-over between the VBL and the TBL thickness, $ f(x) = (1 + x^4)^{-1/4}$ , has been omitted. In fact, integrating this function would be straightforward. Following Grossmann & Lohse (Reference Grossmann and Lohse2001) and Zhang et al.’s (Reference Zhang, Ding and Xia2021) application to smoothing function, dividing (3.31) by (3.11) gives $Re^{1/2}/(6aNu)=\delta _T/(3\delta _u)$ . Additionally, by setting $ x = 3\delta _u/\delta _T$ , when $ Pr \ll 1$ , $ f(x) \rightarrow 1$ , whereas for $ Pr \gg 1$ , it approaches $ x^{-1}$ . Therefore, for the full range of $Pr$ , $ \varepsilon _u$ can be expressed as

   (3.36) \begin{equation} \varepsilon _u \sim \frac {\nu ^3 ReRa}{4H^4NuPr}f(6aNuRe^{-1/2}). \end{equation}
   The transition in thermal dissipation balance between the bulk and TBL remains the same as in the initial model. Given that our target working fluid is a low- $Pr$ fluid (liquid metal) and that modelling for high- $Pr$ cases is not generally applicable to scenarios where a magnetic field is applied, the function $f$ is neglected here.

3. (iii) The constraint function $ g(x) = x(1 + x^4)^{-1/4}$ , which accounts for the saturation of the VBL thickness as it reaches the container scale, has been removed. For the high- $Pr$ regime, $Re$ of the large-scale wind decreases as $Pr$ increases and $ \delta _u$ gradually grows to its maximum value, comparable to the characteristic scale $H$ ( $\delta _{u,max}=H/(4\sqrt {Re_c})$ ), where $Re_c$ corresponds to the threshold beyond which a discernible large-scale wind exists in the system. In the case of VC, due to the non-zero mean streamwise velocity, the presence of a turbulent wind is easily identifiable, especially for low- $Pr$ fluids. Therefore, the function $ g$ is not necessary within the current parameter range.

Figure 10.

Comparison of (a) Nusselt number and (b) Reynolds number data from Shishkina (Reference Shishkina2016) and Zwirner et al. (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022) (markers) with theoretical predictions (3.34) and (3.35) (lines).

For the coupled equations (3.34) and (3.35), we can use the available data to determine the values of the four priors. Here, four data points from Shishkina (Reference Shishkina2016) (cylindrical cell, $ \Gamma = 1$ ) and Zwirner et al. (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022) (square cavity, $ \Gamma = 1$ ) are chosen to fit the model (the data obtained from the current simulations align closely with those of Zwirner et al. (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022), as listed in table 1), including $Ra=10^5$ , $10^8$ at $Pr=0.032$ and $Ra=10^6$ , $10^7$ at $Pr=0.1$ . The resulting GL parameters are $c_1=6.009, c_2=0.0565, c_3=0.586, c_4=0.00577$ . Note that the coefficients $ c_i$ are strictly constrained by the aspect ratio and are not universal. As shown in figure 10, the curves represent the functional relationship described by (3.34) and (3.35) and cover the data points quite well. Strictly speaking, the predicted $Re$ at $ Pr = 0.1$ slightly underestimates the values at very high $Ra$ , which may stem from different cell models (cylinder versus square cavity, and geometric constraints can significantly affect momentum transport) or the critical point for scaling transitions occurring around $ Pr = 0.1$ (Shishkina Reference Shishkina2016). Overall, the extended model captures the general trends of $Nu$ and $Re$ as functions of $Ra$ and $Pr$ . Additionally, the fitted parameters reveal more insights: the magnitude difference between the bulk and VBL prefactors in the kinetic dissipation balance equation is approximately $ 10^2$ . The subregion contributions shown in figure 5 indicate a bulk contribution of approximately 40 $\, \%$ and a VBL contribution of 60 $\, \%$ , which are of the same order of magnitude. Considering $Re$ values ranging from $ 10^3$ to $ 10^4$ , the order-of-magnitude difference in the kinetic dissipation expression for subregions ( $ Re^{1/2} \sim 10^2$ ) suggests that multiplying the coefficient by this approximation gives a result that roughly corresponds to the actual partitioning of dissipation. This also aligns with the core principles of the GL theory.

At the end of this section, we would like to briefly evaluate this part of the work. Clearly, the current understanding of transport relationships in classical VC at moderate to high $Ra$ is quite mature and comprehensive. The laminar $ 1/4$ -power law and the turbulent $ 1/3$ -power law are sufficiently effective for practical applications, though the influence of aspect ratio and boundary conditions should not be overlooked. As Ng et al. (Reference Ng, Ooi, Lohse and Chung2015) emphasised, ‘the global kinetic dissipation cannot be determined a priori because the relationship for the buoyancy flux is unclosed’, which can be partially solved by means of a reasonable approximation of the vertical heat flux. Equation (3.36) is essentially appropriate within the range of $ Pr = 0.03 {-} 100$ (for moderate $Ra$ ), and the scaling relations for dissipation in the VBL–TBL-dominant regime have been well validated. Despite this, we are not yet able to make predictions for the other three regimes: VBL–bulk, bulk–TBL and bulk–bulk, due to a lack of fundamental understanding of these issues: the parameter ranges in which they occur, the spatial distribution of dissipation and the behaviour of the profiles, among others. Equation (3.36) cannot be guaranteed to hold in these unexplored regimes. While the VBL–TBL-dominant regime accounts for the majority of results from current simulations and experiments, presenting a complete physical picture of VC under all parameter conditions, with different subregion-dominant regimes, will undoubtedly require much more extensive work.

4.1. Flow characteristics and global relations in a horizontal magnetic field

Figure 11.

Isosurface of $Q_{c}=0.5\times (|\!|A|\!|_F^2-|\!|S|\!|_F^2)$ (the second invariant of the velocity gradient tensor, $S$ is a symmetric tensor and $A$ is an antisymmetric tensor), represents the roll structures, coloured by local temperature: (a) $Ra=10^6$ and $Ha=0$ ; (b) $Ha=100$ ; (c) $Ha=400$ ; (d) $Ha=800$ .

Before expanding on the theory of MC, it is essential to first develop a basic understanding of the flow structure. Figure 11 illustrates the evolution of the vortex structure: when no magnetic field is applied ( $Ha$ = 0), the flow exhibits a certain degree of isotropic turbulence, although the global circulation can still be vaguely identified. The bulk contains various small-scale structures, with corner vortices breaking up and being stretched along the flow direction due to the entrainment imposed by the mean flow. This will further result in irregular deformation of the cross-section of the main vortex. After the initial magnetic field is applied ( $Ha$ = 100), the flow develops clear anisotropy. The side vortices are elongated along the field lines without significant interruptions. The vortex surface becomes smoother, and small-scale complex structures are constrained, though some 3-D features remain, such as weak local instabilities between side vortices, intrinsic deformation of the primary vortex and global oscillations induced by the entrainment of side vortices. With a further increase in the magnetic field ( $Ha \gt 400$ ), the flow transitions into a more pronounced Q2-D state. The primary and side vortices become independent and strictly align along the field lines, forming a more coherent and improved self-organised large-scale circulation. Of interest are the two modes exhibited by the flow oscillations: first, the local oscillation of the corner vortex (recirculation region) itself and, second, the more pronounced global elliptical deformation of the main vortex accompanied by the entrainment effect of the corner vortex. The latter diminishes as $Ha$ increases, eventually reaching an approximate steady state, but this does not imply the complete disappearance of the corner vortex. As confirmed by our previous study (Zhu et al. Reference Zhu, Chen and Ni2024), the phenomenon of flow separation in weakly nonlinear or laminar flows under the influence of a magnetic field is jointly regulated by magnetic damping effects and the curvature of boundary layer transitions.

The definition of the Q2-D state can be quantitatively described by the correlation of physical fields. In the electrode-driven flow (Potherat & Klein Reference Potherat and Klein2014; Chen et al. Reference Chen, Potherat, Ni and Moreau2021), which is a typical paradigm for the study of MHD turbulence, the strength of the imposed magnetic field can control the flow’s dimensionality (with vortices either partially or fully extending across the cavity). For a given vortex structure with a transverse scale ( $l_\perp$ ), the time required for it to be stretched by the magnetic field to a specific longitudinal scale ( $l_\parallel$ ) is

(4.1) \begin{equation} \tau _{2D}=\frac {\rho }{\sigma B^2}(\frac {l_\parallel }{l_\perp }). \end{equation}

Also consider the characteristic motion time of the vortex structure, i.e. the vortex turnover time,

(4.2) \begin{equation} \tau _{u}=l_\perp /u(l_\perp ). \end{equation}

When the two are equal, the characteristic scale of the given structure along the magnetic field direction can be determined,

(4.3) \begin{equation} l_\parallel = \sqrt {N}l_\perp , \end{equation}

where $N=\sigma B^2l_\perp /\rho v(l_\perp )$ denotes the interaction parameter based on the scale $l_\perp$ and the velocity of motion $v(l_\perp )$ . Considering the spatial scale of the flow region along $B$ to be $h$ :

1. (i) $l_\parallel \gg h$ , the flow structure is strictly Q2-D;

2. (ii) $l_\parallel \gtrsim h$ , the flow exhibits weak-3-D characteristics, where the 2-D ‘slices’ along the magnetic field direction maintain approximately equivalent flow topology, but the flow intensity varies significantly due to the influence of the wall;

3. (iii) $l_\parallel \ll h$ , the flow manifests itself as a strong 3-D effect at different locations along the field lines.

Therefore, large-scale structures can reach the Q2-D state more quickly compared with small-scale structures. This can be quantitatively assessed by measuring the correlation at the Hartmann layer, which is fully aligned with the field lines. The correlation functions are given by the following:

(4.4) \begin{align} C_1^{\prime}(x, y)=\frac {\sum _{t=0}^{T}\phi _f(x, y, t)\phi _b(x, y, t)}{\sqrt {\sum _{t=0}^{T}\phi _f^2(x, y, t)\sum _{t=0}^{T}\phi _b^2(x, y, t)}}, \end{align}

(4.5) \begin{align} C_2^{\prime}(x, y)=\frac {\sum _{t=0}^{T}\phi _f(x, y, t)\phi _b(x, y, t)}{\sum _{t=0}^{T}\phi _f^2(x, y, t)}, \qquad\end{align}

where $\phi$ represents the physical field used for statistics, which may include velocity, electric potential, and subscripts $f$ and $b$ denote the front and back walls. According to Potherat & Klein (Reference Potherat and Klein2014), when strong 3-D characteristics are present, $ C_1^{\prime}$ is not equal to 1, but it cannot identify weak 3-D characteristics. For $ C_2^{\prime}$ , as long as there is any 3-D effect – whether weak or strong – the value will deviate from 1. In the current cases, note that MC differs from a unilaterally driving of electrode injection flow, as buoyancy as a driving force actually acts over the entire field. In the absence of a magnetic field, instability is triggered by gravitational internal waves generated at separation points or strong shear in the vertical boundary layer (Khoubani et al. Reference Khoubani, Mohanan, Augier and Flór2024). Due to the constraint of the no-slip walls, the highest convective intensity occurs near the centre. The magnetic field, on one hand, eliminates any parallel gradients, and on the other, attempts to smooth the convective intensity along the field lines to minimise wall effects, confining viscous effects to the extremely thin Hartmann layer.

Figure 12.

(a) Time- and surface-averaged horizontal profiles of the vertical velocity for different $Ha$ near the heated wall. (b) Correlations of the velocity components between front and back Hartmann layers ( $z_{1, 2}$ ) or the Hartmann layer and the centre ( $z_{1, 3}$ ), averaged over several groups.

Thus, in figure 12(b), in addition to measuring the velocity correlations within the relative Hartmann layers ( $ z_1 = H/Ha,\ z_2 = H(1-1/Ha)$ ), it is also necessary to consider the correlation values at HBL relative to the central position ( $ z_3 = H/2$ ) to capture the presence of weak 3-D characteristic. As shown, after the initial magnetic field is applied, the flow structures along the field lines always exhibit $C_1^{\prime}\approx 1$ , indicating the absence of strong 3-D, and the convective vortex topology remains self-similar. However, the features of weak 3-D are quite apparent. The values of $C_2^{\prime}|_{z_{1, 3}}$ , measured between the near-wall and central region, increase from 0.3 to 0.4 and converge to 0.8 as $ Ha$ increases, suggesting non-uniform convective intensity along $ B$ and possible 3-D characteristics in small-scale structures. Additionally, due to the symmetry imposed by the no-slip walls, $C_2^{\prime}|_{z_{1, 2}}$ is also close to 1 (not shown in the present work). So, we can infer that before $C_2^{\prime}|_{z_{1, 3}}$ approaches saturation, the flow is in a ‘weak 3-D’ state, where large-scale Q2-D structures coexist with small-scale 3-D structures. After $ C_2^{\prime}$ reaches convergence, the flow transitions into a more complete ‘strong Q2-D’ state. At the present $ Ra = 1 \times 10^6$ , this critical $ Ha$ lies between 400 and 500. The vortex structures shown in figure 11 provide a more intuitive insight of this discussion on flow dimensionality, setting the stage for the dissipation scaling discussion in § 4.2.

However, in conjunction with the flow images, it is worth emphasising the ‘rectification’ effect of the horizontal magnetic field, i.e. a slight enhancement in convective strength and heat transfer rate at small $Ha$ ( $Nu$ increases by $4.93\, \%$ at $Ha=100$ , as listed in table 5). This process aligns with the inverse energy cascade observed in classical two-dimensional turbulence, where energy moves from the forcing scale to larger-scale structures (Baker et al. Reference Baker, Potherat, Davoust and Debray2018). Specifically, the magnetic field supports system-size flow by directing the energy of small vortices upward towards the primary vortex. From the perspective of flow coherence, this can be interpreted as an increase in horizontal heat transport per unit area by the plumes, essentially generating better-organised, higher-energy clusters of hot and cold plumes (Zhu et al. Reference Zhu, Chen and Ni2024). Figure 12(a) shows the mean vertical velocity profiles (near the heated wall) at $Ra=10^6$ under different $Ha$ . It can be seen that at $Ha=100$ , the velocity of the ‘turbulent wind’ increases by $10.28\, \%$ compared with the profile without magnetic field, indicating that the rectification effect enhances the main circulation by suppressing small-scale 3-D structures. As $Ha$ increases further, the Hartmann braking due to magnetic damping intensifies. If flow coherence cannot continue to improve, the convective intensity decreases as dissipation increases.

Table 5.

Response parameters for MC with a horizontal magnetic field. From the sixth column from the left to the last column on the right indicate the following: the average kinetic energy of the main circulation, the kinetic energy along the field lines, total dissipation, the approximate expression (3.12), the relative errors.

Recalling the momentum equation (2.2), the additional Joule heating effect contributes to the dissipation due to the Lorentz force. Apply a similar treatment as in the previous section to obtain the kinetic energy equation

(4.6) \begin{align} \frac {\partial }{\partial t}(\frac {u^2}{2})&=-\boldsymbol {\nabla }\cdot \left[\left(p+\frac {u^2}{2}\right)\boldsymbol {u}\right]+\nabla \cdot (u_i\tau _{ij})-2\sqrt {\frac {Pr}{Ra}}S_{ij}S_{ij}\notag\\&\quad+Ha^2\sqrt {\frac {Pr}{Ra}}(\boldsymbol {j}\times \boldsymbol {e}_b)\cdot \boldsymbol {u}+u_yT, \end{align}

where $S_{ij}$ is the strain-rate tensor and $\tau _{ij}$ is the viscous stress. Taking volume averaging under steady-state condition omits the non-constant term and the divergence term to obtain the equilibrium relationship between potential energy and dissipation,

(4.7) \begin{equation} u_yT=2\sqrt {\frac {Pr}{Ra}}S_{ij}S_{ij}-Ha^2\sqrt {\frac {Pr}{Ra}}(\boldsymbol {j}\times \boldsymbol {e}_b)\cdot \boldsymbol {u}. \end{equation}

The work done by the Lorentz force can be rewritten as

(4.8) \begin{equation} -Ha^2\sqrt {\frac {Pr}{Ra}}(\boldsymbol {j}\times \boldsymbol {e}_b)\cdot \boldsymbol {u}=Ha^2\sqrt {\frac {Pr}{Ra}}(\boldsymbol {u}\times \boldsymbol {e}_b)\cdot \boldsymbol {j}=Ha^2\sqrt {\frac {Pr}{Ra}}j^2-\boldsymbol {\nabla }\cdot (\phi \boldsymbol {j}). \end{equation}

By removing the divergence terms from the global perspective, we obtain the expression for Joule dissipation under the quasi-static assumption,

(4.9) \begin{equation} \varepsilon _b = Ha^2\sqrt {\frac {Pr}{Ra}}j^2, \end{equation}

which can also be written in general form as

(4.10) \begin{equation} \varepsilon _b = \frac {\eta }{2}(\partial _ib_j+\partial _jb_i)^2. \end{equation}

Figure 13.

$ \langle \varepsilon _u \rangle + \langle \varepsilon _b \rangle$ compared with the exact relation (3.9) and the approximate expression (3.12) for different $Ha$ , combined with the original equation (1.2).

Clearly, the exact relationship (4.7) reflects that in MC, the buoyant potential energy (left-hand side) is entirely consumed by kinetic dissipation (first term on the right-hand side) and the work done by the Lorentz force (second term on the right-hand side). This should hold rigorously. The modelling of the total potential energy on the left-hand side can still follow the approximate expression we used for the non-magnetic case (3.12), as the energy source remains unique and definite – the vertical heat flux. Furthermore, the topology of the velocity and temperature profiles is essentially equivalent to that in the absence of a magnetic field, meaning the core argument that ‘vertical heat flux is primarily concentrated within the thermal boundary layer’ remains valid under the current parameters. Consequently, the dissipation relationship (1.2) and (1.3) can be written as

(4.11) \begin{align} \varepsilon _u+\varepsilon _b=\frac {\nu ^3}{H^4}\frac {\langle u_y T\rangle }{\kappa \Delta /H}RaPr^{-2}\sim \frac {\nu ^3ReRa}{4H^4NuPr}, \end{align}

(4.12) \begin{align} \varepsilon _T=\kappa \frac {\Delta T^2}{H^2}Nu. \qquad\qquad\qquad \end{align}

The conservation relation for thermal dissipation has already been verified in table 2 of § 2. Figure 13 presents the validation of the balance equation and the approximate relation (4.11), with the original equation (1.2) also included for comparison. As expected, the exact balance between potential energy and total dissipation is reproduced, and the approximate expression (3.12) aligns reasonably well with the total potential energy, with an error of approximately $ 14.2 \, \%$ at $ Ha = 200$ , decreasing to $ 4.5\, \%$ as $ Ha$ increases. The error observed at low- $Ha$ might be attributed to the residual unsteady characteristics of the flow, which transitions to a steady laminar state under a strong magnetic field. In contrast, the expression applicable to RBC shows a significantly larger deviation.

Before examining and adjusting the scaling relations of each component in the bulk and boundary layer, a fundamental issue needs to be reconsidered: the definition of the Reynolds number. In non-magnetic VC, it is typically to define the Reynolds number by the ‘turbulent wind’ velocity ( $Re_{wind} = HU_{wind}/\nu$ ) because the non-zero mean flow is easily identifiable, aligning with the basic assumption in GL theory that a large-scale wind exists. However, in practical applications of thermal convection, the definition of $Re$ is not fixed, but tailored to the specific physical model and flow image. For instance, the r.m.s. definition of $Re=H\sqrt {\langle \boldsymbol {u\cdot u}\rangle _{V, t}}/\nu$ has been widely used in validating RBC scaling relations (Emran & Schumacher Reference Emran and Schumacher2015; Scheel & Schumacher Reference Scheel and Schumacher2016; Zhang et al. Reference Zhang, Zhou and Sun2017), and applied in extending GL theory to double-diffusive convection (Yang & Lohse Reference Yang and Lohse2017). In MC research, $Re_{rms}$ is commonly adopted to provide a measure of the overall flow intensity (Zürner et al. Reference Zürner, Liu, Krasnov and Schumacher2016; Lim et al. Reference Lim, Chong, Ding and Xia2019; Yan et al. Reference Yan, Calkins, Maffei, Julien, Tobias and Marti2019; Akhmedagaev et al. Reference Akhmedagaev, Zikanov, Krasnov and Schumacher2020; Zürner et al. Reference Zürner2020; Nicoski, Yan & Calkins Reference Nicoski, Yan and Calkins2022). Grossmann & Lohse (Reference Grossmann and Lohse2002) proposed that $Re$ can be rescaled by any factor without affecting the results for $Nu$ , as the shape of $Re$ as a function in phase space remains the same, aside from a constant shift. Despite this, Wang et al.’s (Reference Wang, Liu, Shishkina and Lohse2021) results suggest that different definitions of $Re$ lead to different scaling relations. To obtain a definition of $Re$ suitable for revising the scaling relations, we start with the subregion contributions derived from the conservation of thermal dissipation, which remains invariant in form, and test (3.5) and (3.6).

Figure 14.

$Re$ dependence of the normalised thermal dissipation rates in the TBL and bulk based on the different definitions: (a) $Re_{rms}$ and (b) $Re_{wind}$ , validated with scaling laws (3.5) and (3.6).

Figure 14 illustrates the $Re$ dependence of thermal dissipation in the bulk and TBL with respect to two different definitions. It can be observed that when the r.m.s. definition of $Re$ is used, the classical scaling relations (3.5) and (3.6) are reproduced. In contrast, using the ‘turbulent wind’ definition results in a stronger $Re$ dependence (with exponents of 1.52 and 1.94), since the reduction in $ Re_{wind}$ is much smaller than that of $ Re_{rms}$ as $Ha$ increases. This stronger dependence can be explained by the rectification effect of the horizontal magnetic field discussed earlier: the main circulation is enhanced due to the energy injection from small-scale structures, offsetting part of the damping effect from Joule dissipation. However, the overall weakening of flow intensity is evident, particularly reflected in the horizontal Reynolds number $ Re_{u_x}$ and the circulation kinetic energy $ E_{xy}$ (see table 5). The horizontal flow, as a product of the transition in the vertical boundary layer of the confined system, is highly sensitive to changes in the magnetic field (with $Ha$ increasing from 200 to 800, $ Re_{u_y}$ decreases by approximately $10\, \%$ , $ Re_{u_x}$ by $24\, \%$ , $ E_{xy}$ by $46\, \%$ and $ Re_{rms}$ by $26\, \%$ ). Given the overall weakening of flow intensity (where the compensation for the main circulation comes at the expense of an almost stagnant central region, as seen in the vortex structures in figure 11), the physical context of horizontal magnetic fields, and the current preference for $Re$ definitions in MC research, using the r.m.s. definition of $Re$ is more appropriate. This is further validated by figure 14 and continues to be followed in subsequent studies of the scaling relations.

4.2. Scaling of the subregion dissipation under a horizontal magnetic field

From (4.7), the total potential energy in the magnetic field environment consists of kinetic dissipation and Joule dissipation, with a constant thermal dissipation relationship. Decomposing these into contributions from the boundary layer and the bulk:

(4.13) \begin{align} \varepsilon _{u, bulk}+\varepsilon _{u, BL}+\varepsilon _{b, bulk}+\varepsilon _{b, BL}\sim \frac {\nu ^3ReRa}{4H^4NuPr}, \end{align}

(4.14) \begin{align} \varepsilon _{T, bulk}+\varepsilon _{T, BL}=\kappa \frac {\Delta T^2}{H^2}Nu, \qquad\qquad \end{align}

where the scalings of the thermal dissipation component have been validated in § 4.1 concerning the definition of $Re$ . The kinetic dissipation scalings without the magnetic fields are (3.3) and (3.4), which will be carefully adjusted in the following according to the actual physical picture of the horizontal magnetic field VC. Additionally, Zürner et al. (Reference Zürner, Liu, Krasnov and Schumacher2016) and Zürner (Reference Zürner2020) suggested that the scalings for the Joule dissipation component should follow $\varepsilon _{b, bulk}\sim Re^2Q$ and $\varepsilon _{b, BL}\sim Re^2Q^{3/2}$ . There are also some new findings in this area that are worth discussing.

The derivation of the $\varepsilon _{u, bulk}$ scaling in the original GL theory follows the Kolmogorov framework of the turbulent energy cascade, where kinetic energy is transferred from large to small scales. This means that the large vortex, occupying the entire volume of the container, stirs the bulk with a characteristic velocity $U$ and a characteristic time scale of $H/U$ . Thus, $\varepsilon _{u, bulk}$ can be modelled as (3.3), which, at the order of magnitude, balances the large-scale convective term in the energy equation. When a horizontal magnetic field is applied, the $Re$ dependence of $\varepsilon _{u, bulk}$ exhibits an entirely different trend. As shown in figure 15(a), $Ha$ and $Re$ are negatively correlated due to the magnetic damping effect. As $Ha$ increases, $Re$ decreases, but $\varepsilon _{u, bulk}$ increases instead, and the $Re^3$ scaling is no longer present.

Figure 15.

(a) $Re$ dependence of the kinetic dissipation in the bulk for different $Ha$ at $Ra=10^6$ , with the $Re^3$ scaling shown as a reference. (b) Comparison between the simulated $\varepsilon _{u, bulk}$ with (4.21). (c,d) $Ha$ and $Re$ dependence of the compensated data.

Let us extend the Kolmogorov energy decay law to low- $Rm$ MHD turbulence. Due to the presence of both viscous and Joule dissipation, the total energy decay (dominated by kinetic energy at low- $Rm$ , as the induced magnetic field $b\sim RmB_0$ is very small) satisfies (neglecting buoyancy for the moment)

(4.15) \begin{equation} \frac {{\rm{d}}E}{{\rm{d}}t}=-\frac {1}{\sigma }\int \boldsymbol {j}^2{\rm{d}}V-\rho \nu \int \boldsymbol {\omega }^2{\rm{d}}V, \end{equation}

where $\omega$ represents the vorticity. Compared with the decay law in classical turbulence, ${\rm{d}}u^2/{\rm{d}}t \sim -u^3/l$ , the kinetic energy decay law in MHD can be obtained as

(4.16) \begin{equation} \frac {{\rm{d}}u^2}{{\rm{d}}t}\sim -\frac {u^3}{l_\perp }-\left (\frac {l_\perp }{l_\parallel }\right ) ^2\frac {u^2}{\tau }. \end{equation}

Here, $\tau$ is the Joule time scale, and $l_\perp$ and $l_\parallel$ denote the integral length scales of turbulence in the directions perpendicular and parallel to the magnetic field, respectively. As discussed earlier in the context of Q2-D flow, $l_\perp$ remains essentially unchanged, while $l_\parallel$ is elongated along the field lines due to momentum diffusion. When $N$ is small, (4.16) reverts to the Kolmogorov picture, and when $N$ becomes sufficiently large, such that the inertial force is negligible compared with the Lorentz force, it reduces to

(4.17) \begin{equation} \frac{ {\rm{d}}u^2}{{\rm{d}}t}\sim -\left (\frac {l_\perp }{l_\parallel }\right )^2\frac {u^2}{\tau }. \end{equation}

Since $N$ is proportional to the integral scale of the turbulence, only the largest vortex can feel the so-called ‘electromagnetic diffusion’ effect (being stretched into columnar vortex aligned with $\boldsymbol {B}$ ). Most small-scale vortices have turnover times much shorter than $ \tau$ . The above is a brief description of statistically homogeneous isotropic turbulence under the influence of a magnetic field, neglecting boundary effects. Considering the current flow in a confined region, the presence of Hartmann walls perpendicular to the field lines introduces a ‘truncation’ effect on vortex stretching, which drives the discussion on the flow dimensionality. (How should the dissipation scaling in the bulk be considered for wall-bounded flows arranged in a Q2-D state along the magnetic field?) Revisiting the momentum equation (2.2), considering insulated walls and averaging along the magnetic field lines yields a two-dimensional equation (Sommeria & Moreau Reference Sommeria and Moreau1982),

(4.18) \begin{equation} \frac {{\rm{d}}\boldsymbol {u}_\perp }{{\rm{d}}t}=\frac {1}{\rho }\nabla _\perp P+\nu \nabla _\perp ^2\boldsymbol {u}_\perp -\frac {nB}{H}\left (\frac {\sigma \nu }{\rho }\right )^{1/2}\boldsymbol {u}_\perp . \end{equation}

The Lorentz force is expressed here as a linear term, viewed as friction indirectly induced by shaping the Hartmann layer. The above equation reflects the physical picture: the electric currents in the Hartmann layer must close within the bulk, thereby generating a resistance in the main flow that transmits the vorticity loss caused by viscous stress at the wall to the core. In other words, the magnetic field’s resistance to the flow is manifested in the shear stress produced by the strong velocity gradient near the Hartmann walls. Here, $n$ represents the number of Hartmann walls, indicating the number of resistive interfaces experienced by the core, which is typically 2 (or 1 in the case of a single free surface). The linear friction term in (4.18) can also be expressed as $ \boldsymbol {u}_\perp /\tau _H$ , where

(4.19) \begin{equation} \tau _H = \frac {H}{2B}\left (\frac {\rho }{\sigma \nu }\right )^{1/2}. \end{equation}

This corresponds to the linear damping time exerted by the magnetic field on the Q2-D flow structures, also known as the ‘Hartmann friction time’. It is important to note that $ \tau _H$ differs from the previously mentioned Joule time $ \tau _j$ , the latter representing the characteristic time of the process in which the velocity gradient and electric current between adjacent fluid layers with velocity differences dissipate along the magnetic field lines.

We choose $ \tau _H$ as the characteristic time for $\varepsilon _{u, bulk}$ and substitute it into (4.16). Considering the aspect ratio of the large-scale vortex in cavity, $l_\perp /l_\parallel \sim 1$ , the following general form of the kinetic energy decay law can be obtained:

(4.20) \begin{equation} \frac {{\rm{d}}u^2}{{\rm{d}}t}\sim -\alpha \frac {u^3}{l_\perp }-\beta \frac {u^2}{\tau _H}. \end{equation}

Quantifying the relative magnitude of the vortex turnover time and Hartmann friction time, $\tau _{u}/\tau _{H}\sim 2Ha/Re$ . Naturally, it can serve as a criterion for the dominant time scale. When $ 2Ha/Re \gt 1$ , under the strong magnetic field approximation, $ \tau _H$ is smaller than $ H/U$ , meaning the energy of the main vortex dissipates on the friction time scale, which is faster than the vortex turnover time. As shown in table 5, this critical regime occurs around $ Ha \sim 500$ at $Ra=10^6$ . At this point, the results of the correlation coefficients (figure 12 b) indicate that the flow has just reached a fully ‘strong Q2-D’ state, which is a necessary condition for $ \tau _H$ to dominate (as also reflected by the vortex structure diagrams, figure 11). Therefore, (4.20) can be simplified to

(4.21) \begin{equation} \varepsilon _{u, bulk}\sim \frac {U^2}{\tau _H}\sim \frac {\nu ^3}{H^4}Re^2Ha. \end{equation}

This indicates that in wall-bounded Q2-D flows, the large-scale vortex energy is dissipated mainly through boundary layer friction on a $\tau _H$ time scale. Figure 15(b) compares the simulated $\varepsilon _{u, bulk}$ with (4.21), showing a good proportional relationship between the two when the magnetic field increases to the predicted value. Before reaching the critical regime, the flow field still exhibits some three-dimensional characteristics, and the direct energy cascade remains non-negligible. The actual dissipation deviates from the theoretical scaling to some extent, but converges as $Ha$ increases, as also shown in the compensated data plots in figures 15(c) and 15(d). A more general form of the decay law might better capture dissipation scaling over a broader range, but this would likely require introducing two additional prefactors to quantify their respective contributions, as previous studies on independent scaling laws have shown that the coefficients are often not equal to 1. For example, in ${\rm{d}}u^2/{\rm{d}}t\sim \alpha u^3/l$ , experiments found $ \alpha = 1.0 {-} 1.2$ , and in isotropic turbulence, $ \beta = 2/3$ in ${\rm{d}}u^2/{\rm{d}}t\sim -\beta (l_\perp /l_\parallel )^2u^2/\tau$ (Davidson Reference Davidson2001b ). Thus, to minimise the introduction of extra parameters and reduce constraints on the final model equations, we continue to use (4.21) as a proxy for the statistical $\varepsilon _{u, bulk}$ , though this imposes limitations on the model’s applicability. It is worth noting that figure 15(b) exhibits not only a proportional relationship but also a surprisingly close agreement in absolute magnitude, with a ratio of approximately 1.3.

It should be emphasised that (4.21) shares the same form as Zürner et al.’s (Reference Zürner, Liu, Krasnov and Schumacher2016) modelling of $\varepsilon _{u, BL}$ for vertical magnetic field RBC, but these are two fundamentally different concepts. First, they apply to different regions: the current case concerns the bulk, while the latter applies to the HBL. Second, the physical nature is different: the current theory is derived from the kinetic energy decay caused by linear friction of Q2-D vortices at the wall, whereas Zürner et al.’s (Reference Zürner, Liu, Krasnov and Schumacher2016) model is based on the HBL thickness scaling as $ Ha^{-1}$ , obtained through order-of-magnitude approximation. Therefore, their intrinsic meanings differ and it is crucial to distinguish between the two.

The correction for $\varepsilon _{u, BL}$ is based on the traditional GL theory perspective: the large-scale ‘wind’ in VC is apparent, and thus $\varepsilon _{u, BL}$ is closely related to the volume fraction of the VBL and the magnitude of the velocity gradient within it. Unlike in vertical magnetic field RBC, where the velocity BL at the hot and cold plates serve as HBL, in a horizontal magnetic field, they are actually ‘side layers’, also referred to as ‘parallel layers’. Since the magnetic field eliminates any gradients along the field lines, $\varepsilon _{u, BL}$ primarily occurs in the plane of the main circulation perpendicular to the field, which includes both the vertical and horizontal BLs. At low $Ha$ , the thickness of the parallel layers does not follow the asymptotic form $\sim bH/\sqrt {Ha}$ , but should revert to the classical Prandtl–Blasius scaling (3.2), $\sim aH/\sqrt {Re}$ . Following Lim et al. (Reference Lim, Chong, Ding and Xia2019) and Zürner (Reference Zürner2020), we introduce a combination of both to estimate the thickness,

(4.22) \begin{equation} \delta \sim (\delta _u^{-2}+\delta _{s}^{-2})^{-1/2}\sim \frac {H}{(a^{-2}Re+b^{-2}Ha)^{1/2}}. \end{equation}

A plain idea would be to fit the model with the actual thickness, but this would introduce two additional factors $a$ and $b$ , increasing the model’s complexity and sensitivity to parameters. Subsequent results prove that this is unnecessary, as an exact match of numerical values is not required as long as the fundamental physical processes are captured. For practical purposes, achieving a good dissipation scaling through dimensional analysis is sufficient. A brief discussion of the choice of P-B scaling prefactor for VC is given here. Zwirner et al.’s (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022) DNS results for $ \Gamma = 5$ and $ Pr = 0.032$ yielded $ a = 0.38$ , which was subsequently well validated by experimental results for $ \Gamma = 3$ and 5. Ng et al.’s (Reference Ng, Ooi, Lohse and Chung2015) DNS results for $ Pr = 0.71$ , using periodic boundaries, gave $ a = 0.43$ , which encompassed different definitions of the $\delta _u$ scaling. We will use these two values, while the scaling coefficient $ b$ for the parallel layer is set to 1 without any further adjustments. Figure 16(a) illustrates the comparison between the actual VBL thickness and the calculated values from (4.22). The predicted values closely match $\delta _{u_y}$ , while $\delta _{u_x}$ is significantly larger. This is caused by the fact that horizontal VBL is driven by the transition of vertical VBL under geometric constraints, lacking the steep velocity gradient directly induced by buoyancy. Moreover, the difference between 0.43 and 0.38 is minimal, and the model is not particularly sensitive to this variation. It is also worth noting that there is still a deviation in trend between (4.22) and $\delta _{u_y}$ , which stems from setting the prefactor $b$ for $Ha$ to 1, thus somewhat underestimating the influence of the magnetic field. In reality, $\delta _{u_y}$ decreases more rapidly as $Ha$ increases (in the original two-factor fit, $ b$ was approximately 0.41, though this was not shown in the paper). Even so, the extrapolated dissipation expression,

(4.23) \begin{equation} \varepsilon _{u, BL}\sim \nu \frac {U^2}{\delta ^2}\frac {\delta }{H}\sim \frac {\nu ^3}{H^4}(a^{-2}Re^5+Re^4Ha)^{1/2}, \end{equation}

still maintains a relatively proportional relation with the actual values (figure 16 b). The above analysis reflects an important physical insight: the influence of the magnetic field on $\varepsilon _{u, BL}$ does not change the time scale of the process, as it does in the bulk. Instead, it regulates the magnitude of dissipation by affecting the boundary layer thickness, i.e. the essence that $\varepsilon _{u, BL}$ is determined by the velocity gradient remains unchanged. The transition between buoyancy-dominated and Lorentz-force-dominated regimes in vertical magnetic field RBC corresponds to the cross-over of boundary layer thickness (Teimurazov et al. Reference Teimurazov, McCormack, Linkmann and Shishkina2024).

Figure 16.

(a) Comparison of the parallel layer thickness $\delta$ , including the vertical VBL $\delta _{u_y}$ , and the horizontal VBL $\delta _{u_x}$ , with the conjugate expression (4.22). (b) Ratio of the simulated $\varepsilon _{u,BL}$ to the extrapolated expression (4.23).

After the above discussion on the subregion scalings of kinetic dissipation, it is now possible to turn our attention to the analysis of Joule dissipation. Following Zürner et al.’s (Reference Zürner, Liu, Krasnov and Schumacher2016) suggestion, the scalings of Joule dissipation in both the bulk and HBL can also be derived through basic order-of-magnitude analysis,

(4.24) \begin{align} \varepsilon _{b, bulk} \sim \eta \frac {Rm^2B_0^2}{H^2}\sim \frac {\nu ^3}{H^4}Re^2Ha^2, \,\,\,\end{align}

(4.25) \begin{align} \varepsilon _{b, BL} \sim \eta \frac {Rm^2B_0^2}{\delta _{Ha}^2}\frac {\delta _{Ha}}{H}\sim \frac {\nu ^3}{H^4}Re^2Ha^3. \end{align}

To accurately compute the subregion’s Joule dissipation in simulations, it is essential to obtain the actual HBL thickness. Typically, for simplification, the $ \delta _{Ha} \sim H/Ha$ approximation is often applied, a scaling derived from the 2-D isothermal pressure-driven Poiseuille flow subjected to a vertical magnetic field. However, in buoyancy-driven or 3-D confined flows, deviations are expected, as seen in the conjugate expressions from Lim et al. (Reference Lim, Chong, Ding and Xia2019) and Zürner (Reference Zürner2020). Therefore, we directly determine the actual BL thickness through the velocity profile along the field lines. Note that in the current Q2-D flow, the plane perpendicular to the field line cannot be surface averaged because the vertical velocity directions on both sides are opposite. Therefore, the profile along $z$ at the $x$ position in the region of the vertical velocity peak is taken, as shown in figure 17(a). Although the $ \delta _{Ha}$ values obtained using the slope method are local (imposing a somewhat stricter constraint on thickness calculation), they are still reasonable because the strongest dissipation occurs where the currents accumulate, i.e. near the velocity peak. Equation (4.10) is not statistically tractable in the simulations, so we turn to calculating the work done by the Lorentz force in the bulk and HBL (4.8).

Figure 17.

(a) Distribution of the mean vertical velocity along the $z$ -axis at the peak of the circulation velocity to define the actual Hartmann layer thickness. (b) Hartmann solution for pressure-driven flow between two parallel plates under a transverse magnetic field; when $Ha=0$ , it would return to a parabolic solution ( $u_{max}=1.5u_{mean}$ ). Note the topological shifts in the centre and near-wall of the profile after applying a magnetic field.

The results presented in figure 18(a) are particularly noteworthy: the Lorentz force performs negative work in the bulk, leading to positive dissipation, while it exerts positive work in the HBL, corresponding to negative dissipation. The absolute values of both increase with $Ha$ , and the sum represents a net negative work, which manifests as dissipation of potential energy. This should not be a surprising result, analogous to Hartmann flow between two parallel plates (the velocity profile is shown in figure 17 b), as described by Shercliff (Reference Shercliff1956): ‘The imposition of the transverse field generates currents which retard the fluid near the centre of the channel and accelerate the fluid near the walls, flattening the normally parabolic velocity profile until finally there is a core of uniform flow and an exponential boundary layer’. Closed loops generate opposite current directions at specific positions in the flow field, resulting in differences in the work done.

Returning to the current MC model, taking the heated wall as an example: the buoyancy-driven vertical motion (y-positive) interacts with $B$ to generate $x$ -negative direction currents, which flow through the wall due to the geometrical constraints to form loops ( $x$ -positive, shown in figure 1). At this point, the currents in the bulk and near the Hartmann walls flow in opposite directions, resulting in negative and positive work. The specific currents distribution and the profiles of opposing Lorentz forces within the VC system have been detailed by Zhu et al. (Reference Zhu, Chen and Ni2024), so they are not elaborated here. From a different perspective, this result can be understood in terms of the dissipation components distinguished in (4.8). When considering the global work done by the Lorentz force, the divergence term $\nabla \cdot (\phi \boldsymbol {j})$ is neglected because no current enters or exits the system. Despite this, after dividing the bulk and HBL, the virtual boundary still exists within the flow field and $\phi \boldsymbol {j}$ has a net flux at the boundary, meaning the divergence integral is no longer zero. As depicted in figure 18(b), $j^2$ consistently represents positive dissipation in both the bulk and the HBL, which is undeniable. The divergence of $\phi \boldsymbol {j}$ is positive in the bulk, corresponding to the previously mentioned buoyancy-driven flow near the centre generating current that flows from the bulk to the HBL, acting as a current ‘source’. The divergence of $\phi \boldsymbol {j}$ in the HBL is negative, representing the current flowing from the bulk into the HBL to form a loop, acting as a ‘sink’. The absolute values of the two are equal and their sum yields a net contribution of zero across the entire system.

Figure 18.

(a) Statistics of the work done by the Lorentz force in the bulk and HBL. (b) The actual magnitude of the Joule dissipation component from (4.8).

In brief, flow in the bulk generates electricity which drives flow in the Hartmann boundary layers, while in this confined system, energy loss caused by Joule heating always exists. If only the $j^2$ term is considered when calculating Joule dissipation in subregions, it will underestimate the kinetic energy consumption by the Lorentz force in the bulk and inaccurately assess its effect within the HBL. This brings to mind the relationship between a generator (bulk) and a motor (HBL), where, in addition to the inevitable generation of Joule heating, there is also a conversion between electrical and kinetic energy.

4.3. Establishment and validation of total model equations under horizontal magnetic field

By combining the scalings of dissipation components from the previous section, their general expressions are incorporated into the approximate balance equations (4.13) and (4.14) in an interpolated form. Upon eliminating the coefficients $\nu ^3/H^4$ and $\kappa \Delta ^2/H^2$ , respectively, the final result is derived:

(4.26) \begin{align} \frac {ReRa}{4NuPr} &= c_{1}(a^{-2}Re^5+Re^4Ha)^{1/2}+c_{2}Re^2Ha \notag\\ &\quad + c_{5}Re^2Ha^3+c_{6}Re^2Ha^2, \end{align}

(4.27) \begin{equation} Nu-1=c_{3}Pr^{1/2}Re^{1/2}+c_{4}PrRe. \end{equation}

Note that here, the Reynolds number is still defined according to the root-mean-square definition, $Re_{rms}$ . Compared with the non-magnetic VC equations (3.34) and (3.35), the MC model modifies two components of the kinetic dissipation rate, introduces Joule dissipation, and includes prefactors $ c_5$ and $ c_6$ . A brief description is as follows: the thermal dissipation relation in (4.14) remains unchanged because the Joule heating effect is negligible compared with the overall temperature gradient, and the global thermal dissipation is accurate. The energy dissipation relation in (4.13) is approximately matched, with the left-hand side representing global dissipation based on the fact that ‘the vertical heat flux providing buoyant potential energy is primarily concentrated within the boundary layer’. The expression for $\varepsilon _{u, BL}$ introduces no new coefficients and demonstrates good Prandtl number generalisation, while $\varepsilon _{u, bulk}$ takes into account the time scale of the main vortex dissipation once the flow reaches a Q2-D state under strong magnetic field.

It is important to note the two terms for Joule dissipation. Based on the analysis in the previous section, due to the differences in the work done by the Lorentz force in the subregions, $ c_5$ should be constrained to a negative value, while $ c_6$ should be positive. Additionally, Zürner’s (Reference Zürner2020) fit for $ c_5$ in vertical magnetic field RBC yielded a value of $ 7.52 \times 10^{-18}$ , which is effectively zero, suggesting that $\varepsilon _{b, BL}$ has no impact on the result. To test this conclusion and complete the model equations, we perform two sets of fits for the prefactors using the current dataset of $Ra$ and $Ha$ (in addition to the standard cases of $ Ra = 1 \times 10^6$ , results for $ Ra = 5 \times 10^5$ and $ 2 \times 10^6$ are also included): (i) setting $ c_5 = 0$ and ensuring all other coefficients are positive; (ii) allowing $ C_5$ to vary freely while ensuring all other coefficients remain positive. The two sets of GL parameters obtained are: $c_1=8.601$ , $c_2=0.597$ , $c_3=0.331$ , $c_4=0.0706$ , $c_5=0$ and $c_6=3.84\times 10^{-4}$ ; $C_1=8.488$ , $C_2=0.635$ , $C_3=0.331$ , $C_4=0.0706$ , $C_5=-1.731\times 10^{-7}$ and $C_6=4.61\times 10^{-4}$ . From the fitting results, when $ C_5$ is not constrained, the fitted parameter is a negative value, consistent with the earlier conclusion that the Lorentz force does positive work in the BL. The corresponding exponent of the $ Ha^3$ scaling is large, actually, this portion of the dissipation contributes only a small fraction (figure 18 b). Therefore, $ C_5$ must be a small quantity close to zero.

Figure 19.

Comparison of the Nusselt number and Renolds number ( $Re_{rms}$ ) data from the present simulations and Chen, Liu & Ni (Reference Chen, Liu and Ni2018) with theoretical predictions from (4.26)–(4.27) (lines). (a,b) Disregard $\varepsilon _{b, BL}$ and let $c_5$ be 0. (c,d) Consider $\varepsilon _{b, BL}$ but do not constrain $C_5$ .

The predicted curves from these two sets of model equations, along with the data points, are plotted in figure 19, with Chen et al.’s (Reference Chen, Liu and Ni2018) results for $ Pr = 0.025$ included as a supplement. It can be observed that the theoretical predictions and the numerical simulations of both $Nu$ and $Re$ show a generally agreement. Notably, we only use data for $ Pr = 0.032$ in the fitting process; however, the resulting model equations still capture subtle differences in $Pr$ -dependence and provide reasonably accurate predictions for $ Pr = 0.025$ . It must be acknowledged that these two Prandtl numbers are quite close and some overlap between the two datasets can be expected. Nevertheless, the smaller $Pr$ exhibits a gentler decline and overall lower $Nu$ values.

Finally, two aspects deserve discussing.

1. (i) Considering $\varepsilon _{b, BL}$ versus not considering it shows little difference in the overall transport relationship. Even when included in the model equation, it does not significantly improve the prediction, as the function without it already sufficiently describes the behaviour and the value of $C_5$ represents a very small weighting factor. It is also possible that as $Ha$ increases further, the high power growth of $ Ha^3$ leads to numerical instability. Therefore, we recommend following Zürner’s (Reference Zürner2020) suggestion in practical applications and setting $ c_5 = 0$ , which is consistent with the physical picture of viscous dissipation being the dominant mechanism of energy loss in horizontal magnetic field MC (Chen et al. Reference Chen, Wang and Ni2024).

2. (ii) The overall decreasing trend of $Nu$ and $Re$ slows slightly as $Ra$ increases, which was an important consideration when selecting the parameters for this study and also serves as the starting point for correcting viscous dissipation in the bulk. In comparison to vertical magnetic field RBC, where convection is almost completely suppressed at $ Ra = 6.3 \times 10^6$ and $ Ha\sim 1000$ , with $ Nu\rightarrow 1$ (Zürner Reference Zürner2020), the convective intensity in horizontal magnetic field VC under similar conditions is much higher. The energy loss under vertical magnetic field is basically dominated by Joule dissipation on a much faster time scale, while the damping effect of horizontal magnetic field on convective motion is not as pronounced. For very large $Ra$ , the influence of the magnetic field may not be clearly manifested, and this would not align with our correction of bulk scaling for a strong Q2-D flow regime. Therefore, we must exercise caution regarding the applicability of the current theoretical model: (i) global energy dissipation can be approximately described by (3.12); (ii) the magnetic damping effect on convective intensity can be clearly quantified.